Philippe 07/03/2019: IO bugfix: io_set_mnhversion must be called by all the processes
[MNH-git_open_source-lfs.git] / src / MNH / advecuvw_rk.f90
1 !MNH_LIC Copyright 1994-2014 CNRS, Meteo-France and Universite Paul Sabatier
2 !MNH_LIC This is part of the Meso-NH software governed by the CeCILL-C licence
3 !MNH_LIC version 1. See LICENSE, CeCILL-C_V1-en.txt and CeCILL-C_V1-fr.txt  
4 !MNH_LIC for details. version 1.
5 !-----------------------------------------------------------------
6 !     #####################
7       MODULE MODI_ADVECUVW_RK
8 !     #####################
9 !
10 INTERFACE
11       SUBROUTINE ADVECUVW_RK (HUVW_ADV_SCHEME,                         &
12                     HTEMP_SCHEME, KWENO_ORDER,                         &
13                     HLBCX, HLBCY, PTSTEP,                              &
14                     PU, PV, PW,                                        &
15                     PUT, PVT, PWT,                                     &
16                     PMXM_RHODJ, PMYM_RHODJ, PMZM_RHODJ,                &
17                     PRUCT, PRVCT, PRWCT,                               &
18                     PRUS_ADV, PRVS_ADV, PRWS_ADV,                      &
19                     PRUS_OTHER, PRVS_OTHER, PRWS_OTHER                 )
20 !
21 CHARACTER(LEN=6),         INTENT(IN)    :: HUVW_ADV_SCHEME! to the selected
22 CHARACTER(LEN=4),         INTENT(IN)    :: HTEMP_SCHEME   ! Temporal scheme
23 !
24 INTEGER,                  INTENT(IN)    :: KWENO_ORDER   ! Order of the WENO
25                                                          ! scheme (3 or 5)
26 !
27 CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
28 !
29 REAL,                     INTENT(IN)    :: PTSTEP
30 !
31 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PU , PV  , PW
32                                                   ! Variables to advect
33 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT, PVT , PWT
34                                                   ! Variables at t
35 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMXM_RHODJ
36 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMYM_RHODJ
37 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMZM_RHODJ
38                                                   !  metric coefficients
39 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUCT , PRVCT, PRWCT
40                                                   ! Contravariant wind components
41 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRUS_ADV , PRVS_ADV, PRWS_ADV
42                                                   ! Tendency due to advection
43 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER
44 !                                                 ! tendencies from other processes
45 !
46 END SUBROUTINE ADVECUVW_RK
47 !
48 END INTERFACE
49 !
50 END MODULE MODI_ADVECUVW_RK
51 !     ##########################################################################
52       SUBROUTINE ADVECUVW_RK (HUVW_ADV_SCHEME,                         &
53                     HTEMP_SCHEME, KWENO_ORDER,                         &
54                     HLBCX, HLBCY, PTSTEP,                              &
55                     PU, PV, PW,                                        &
56                     PUT, PVT, PWT,                                     &
57                     PMXM_RHODJ, PMYM_RHODJ, PMZM_RHODJ,                &
58                     PRUCT, PRVCT, PRWCT,                               &
59                     PRUS_ADV, PRVS_ADV, PRWS_ADV,                      &
60                     PRUS_OTHER, PRVS_OTHER, PRWS_OTHER                 )
61 !     ##########################################################################
62 !
63 !!****  *ADVECUVW_RK * - routine to call the specialized advection routines for wind
64 !!(:,:,IKE+
65 !!    PURPOSE
66 !!    -------
67 !!
68 !!**  METHOD
69 !!    ------
70 !!
71 !!    EXTERNAL
72 !!    --------
73 !!
74 !!    IMPLICIT ARGUMENTS
75 !!    ------------------
76 !!      NONE
77 !!
78 !!    REFERENCE
79 !!    ---------
80 !!      Book1 and book2 ( routine ADVECTION )
81 !!
82 !!    AUTHOR
83 !!    ------
84 !!      J.-P. Pinty      * Laboratoire d'Aerologie*
85 !!      J.-P. Lafore     * Meteo France *
86 !!
87 !!    MODIFICATIONS
88 !!    -------------
89 !!      Original    06/07/94 
90 !!                  01/04/95 (Ph. Hereil J. Nicolau) add the model number
91 !!                  23/10/95 (J. Vila and JP Lafore) advection schemes scalar
92 !!                  16/01/97 (JP Pinty)              change presentation 
93 !!                  30/04/98 (J. Stein P Jabouille)  extrapolation for the cyclic
94 !!                                                   case and parallelisation
95 !!                  24/06/99 (P Jabouille)           case of NHALO>1
96 !!                  25/10/05 (JP Pinty)              4th order scheme
97 !!                  24/04/06 (C.Lac)                 Split scalar and passive
98 !!                                                   tracer routines
99 !!                  08/06    (T.Maric)               PPM scheme
100 !!                  04/2011  (V. Masson & C. Lac)    splits the routine and adds
101 !!                                                   time splitting
102 !!                  J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
103 !!                  J.Escobar : 15/09/2015 : WENO5 & JPHEXT <> 1
104 !!                  F.Auguste and C.Lac : 08/16 : CEN4TH with RKC4
105 !!                  C.Lac   10/16 : Correction on RK loop
106 !!
107 !-------------------------------------------------------------------------------
108 !
109 !*       0.    DECLARATIONS
110 !              ------------
111 !
112 USE MODE_ll
113 USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll
114 USE MODD_PARAMETERS,  ONLY : JPVEXT
115 USE MODD_CONF,        ONLY : NHALO
116 !
117 USE MODI_SHUMAN
118 USE MODI_ADVECUVW_WENO_K
119 USE MODI_ADV_BOUNDARIES
120 USE MODI_GET_HALO
121 USE MODE_MPPDB
122 !
123 USE MODI_ADVECUVW_4TH
124 !
125 !-------------------------------------------------------------------------------
126 !
127 IMPLICIT NONE
128 !
129 !*       0.1   Declarations of dummy arguments :
130 !
131 CHARACTER(LEN=6),         INTENT(IN)    :: HUVW_ADV_SCHEME! to the selected
132 CHARACTER(LEN=4),         INTENT(IN)    :: HTEMP_SCHEME   ! Temporal scheme
133 !
134 INTEGER,                  INTENT(IN)    :: KWENO_ORDER   ! Order of the WENO
135                                                          ! scheme (3 or 5)
136 !
137 CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
138 !
139 REAL,                     INTENT(IN)    :: PTSTEP
140 !
141 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PU , PV  , PW
142                                                   ! Variables to advect
143 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT, PVT , PWT
144                                                   ! Variables at t
145 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMXM_RHODJ
146 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMYM_RHODJ
147 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMZM_RHODJ
148                                                   !  metric coefficients
149 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUCT , PRVCT, PRWCT
150                                                   ! Contravariant wind components
151 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRUS_ADV , PRVS_ADV, PRWS_ADV
152                                                   ! Tendency due to advection
153 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER
154 !                                                 ! tendencies from other processes
155 !
156 !
157 !
158 !*       0.2   declarations of local variables
159 !
160 !
161 !  
162 INTEGER             :: IKE       ! indice K End       in z direction
163 !
164 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZUT, ZVT, ZWT
165 ! Intermediate Guesses inside the RK loop              
166 !
167 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZRUS,ZRVS,ZRWS
168 ! Momentum tendencies due to advection
169 REAL, DIMENSION(:,:), ALLOCATABLE :: ZBUT ! Butcher array coefficients
170                                           ! at the RK sub time step
171 REAL, DIMENSION(:),   ALLOCATABLE :: ZBUTS! Butcher array coefficients
172                                           ! at the end of the RK loop
173
174 !JUAN
175 TYPE(LIST_ll), POINTER      :: TZFIELDMT_ll ! list of fields to exchange
176 TYPE(HALO2LIST_ll), POINTER :: TZHALO2MT_ll ! momentum variables
177 INTEGER                     :: INBVAR
178 INTEGER :: IIU, IJU, IKU ! array sizes
179 !JUAN
180
181 ! Momentum tendencies due to advection
182 INTEGER :: ISPL                ! Number of RK splitting loops
183 INTEGER :: JI, JS              ! Loop index
184 !
185 INTEGER                     :: IINFO_ll    ! return code of parallel routine
186 TYPE(LIST_ll), POINTER      :: TZFIELD_ll  ! list of fields to exchange
187 TYPE(LIST_ll), POINTER      :: TZFIELDS_ll ! list of fields to exchange
188 TYPE(LIST_ll), POINTER      :: TZFIELDS0_ll ! list of fields to exchange
189 TYPE(LIST_ll), POINTER      :: TZFIELDS4_ll ! list of fields to exchange
190 !
191 !
192 REAL          :: XPRECISION
193 !-------------------------------------------------------------------------------
194 !
195 !*       0.     INITIALIZATION                        
196 !               --------------
197 !
198 IKE = SIZE(PWT,3) - JPVEXT
199 IIU=SIZE(PUT,1)
200 IJU=SIZE(PUT,2)
201 IKU=SIZE(PUT,3)
202 !
203 SELECT CASE (HTEMP_SCHEME)
204  CASE('RK11')
205   ISPL = 1
206  CASE('RK21')
207   ISPL = 2
208  CASE('NP32')
209   ISPL = 3
210  CASE('SP32')
211   ISPL = 3
212  CASE('RK33')
213   ISPL = 3
214  CASE('RKC4')
215   ISPL = 4
216  CASE('RK4B')
217   ISPL = 4
218  CASE('RK53')
219   ISPL = 5
220  CASE('RK62')
221   ISPL = 6
222  CASE('RK65')
223   ISPL = 6
224  CASE DEFAULT
225   PRINT *,'ERROR: UNKNOWN HTEMP_SCHEME'
226   CALL ABORT()
227 END SELECT
228 !
229 !
230 ALLOCATE(ZBUT(ISPL-1,ISPL-1))
231 ALLOCATE(ZBUTS(ISPL))
232 !
233 SELECT CASE (HTEMP_SCHEME)
234   CASE('RK11')
235     ZBUTS = (/ 1. /)
236   CASE('RK21')
237     ZBUTS     = (/ 0. , 1. /)
238     ZBUT(1,1)   = 3./4.
239   CASE('RK33')
240     ZBUTS     = (/ 1./6. , 1./6. , 2./3. /)
241     ZBUT(1,1) = 1.
242     ZBUT(1,2) = 0.
243     ZBUT(2,1) = 1./4.
244     ZBUT(2,2) = 1./4.
245   CASE('NP32')
246     ZBUTS     = (/ 1./2. , 0., 1./2. /)
247     ZBUT(1,1) = 1./3.
248     ZBUT(1,2) = 0.
249     ZBUT(2,1) = 0.
250     ZBUT(2,2) = 1.
251   CASE('SP32')
252     ZBUTS     = (/ 1./3. , 1./3. , 1./3. /)
253     ZBUT(1,1) = 1./2.
254     ZBUT(1,2) = 0.
255     ZBUT(2,1) = 1./2.
256     ZBUT(2,2) = 1./2.
257   CASE('RKC4')
258     ZBUTS     = (/ 1./6. , 1./3. , 1./3. , 1./6./)
259     ZBUT      = 0.
260     ZBUT(1,1) = 1./2.
261     ZBUT(2,2) = 1./2.
262     ZBUT(3,3) = 1.
263   CASE('RK4B')
264     ZBUTS     = (/ 1./8. , 3./8. , 3./8. , 1./8./)
265     ZBUT      = 0.
266     ZBUT(1,1) = 1./3.
267     ZBUT(2,1) = -1./3.
268     ZBUT(2,2) = 1.
269     ZBUT(3,1) = 1.
270     ZBUT(3,2) = -1.
271     ZBUT(3,3) = 1.
272   CASE('RK53')
273     ZBUTS     = (/ 1./4. , 0. , 0. , 0. , 3./4. /)
274     ZBUT      = 0.
275     ZBUT(1,1) = 1./7.
276     ZBUT(2,2) = 3./16.
277     ZBUT(3,3) = 1./3.
278     ZBUT(4,4) = 2./3.
279   CASE('RK62')
280     ZBUTS     = (/ 1./6. , 1./6. , 1./6. , 1./6. , 1./6. , 1./6. /)
281     ZBUT      = 0.
282     ZBUT(1,1) = 1./5.
283     ZBUT(2,1) = 1./5.
284     ZBUT(2,2) = 1./5.
285     ZBUT(3,1) = 1./5.
286     ZBUT(3,2) = 1./5.
287     ZBUT(3,3) = 1./5.
288     ZBUT(4,1) = 1./5.
289     ZBUT(4,2) = 1./5.
290     ZBUT(4,3) = 1./5.
291     ZBUT(4,4) = 1./5.
292     ZBUT(5,1) = 1./5.
293     ZBUT(5,2) = 1./5.
294     ZBUT(5,3) = 1./5.
295     ZBUT(5,4) = 1./5.
296     ZBUT(5,5) = 1./5.
297 CASE('RK65')
298     ZBUTS= (/ 7./90. , 0. , 16./45. , 2./15. , 16./45. , 7./90. /)
299     ZBUT= 0.
300     ZBUT(1,1) = 1./4.
301     ZBUT(2,1) = 1./8.
302     ZBUT(2,2) = 1./8.
303     ZBUT(3,1) = 0
304     ZBUT(3,2) = -1./2.
305     ZBUT(3,3) = 1
306     ZBUT(4,1) = 3./16.
307     ZBUT(4,2) = 0
308     ZBUT(4,3) = 0
309     ZBUT(4,4) = 9./16.
310     ZBUT(5,1) = -3./7.
311     ZBUT(5,2) = 2./7.
312     ZBUT(5,3) = 12./7.
313     ZBUT(5,4) = -12./7.
314     ZBUT(5,5) = 8./7.
315 END SELECT
316 !
317 ALLOCATE(ZRUS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
318 ALLOCATE(ZRVS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
319 ALLOCATE(ZRWS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
320 !
321 PRUS_ADV = 0.
322 PRVS_ADV = 0.
323 PRWS_ADV = 0.
324 !
325 !-------------------------------------------------------------------------------
326 !
327 !*       2.     Wind guess before RK loop
328 !               -------------------------
329 !
330 ZUT = PU
331 ZVT = PV
332 ZWT = PW
333 !
334 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
335 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
336 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
337
338 NULLIFY(TZFIELDMT_ll)    
339 CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZUT)
340 CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZVT)
341 CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZWT)
342 INBVAR = 3
343 CALL INIT_HALO2_ll(TZHALO2MT_ll,INBVAR,SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3))
344 !
345 ZRUS = 0.
346 ZRVS = 0.
347 ZRWS = 0.
348 !-------------------------------------------------------------------------------
349 !
350 !*       3.     BEGINNING of Runge-Kutta loop
351 !               -----------------------------
352 !
353  DO JS = 1, ISPL
354 !               
355     CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
356     CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
357     CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
358 !    
359     CALL UPDATE_HALO_ll(TZFIELDMT_ll,IINFO_ll)        
360     CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll)
361 !
362 !*       4.     Advection with WENO
363 !               -------------------
364 !
365
366   IF (HUVW_ADV_SCHEME=='WENO_K') THEN                                                                         
367     CALL ADVECUVW_WENO_K (HLBCX, HLBCY, KWENO_ORDER, ZUT, ZVT, ZWT,     &
368                         PRUCT, PRVCT, PRWCT,                            &
369                         ZRUS(:,:,:,JS), ZRVS(:,:,:,JS), ZRWS(:,:,:,JS), &
370                         TZHALO2MT_ll                                    )
371   ELSE IF ((HUVW_ADV_SCHEME=='CEN4TH') .AND. (HTEMP_SCHEME=='RKC4')) THEN 
372     CALL ADVECUVW_4TH (HLBCX, HLBCY, PRUCT, PRVCT, PRWCT,               &
373                        ZUT, ZVT, ZWT,                                   &
374                        ZRUS(:,:,:,JS), ZRVS(:,:,:,JS), ZRWS(:,:,:,JS),  &
375                        TZHALO2MT_ll )
376   ENDIF
377 !
378      NULLIFY(TZFIELDS4_ll)
379 !
380     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRUS(:,:,:,JS))
381     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRVS(:,:,:,JS))
382     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRWS(:,:,:,JS))
383     CALL UPDATE_HALO_ll(TZFIELDS4_ll,IINFO_ll)
384     CALL CLEANLIST_ll(TZFIELDS4_ll)
385 !               
386   IF ( JS /= ISPL ) THEN
387 !
388       ZUT = PU
389       ZVT = PV
390       ZWT = PW
391 !
392    DO JI = 1, JS
393 !
394 ! Intermediate guesses inside the RK loop
395 !
396       ZUT(:,:,:) = ZUT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
397        ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ
398       ZVT(:,:,:) = ZVT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
399        ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ
400       ZWT(:,:,:) = ZWT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
401        ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ
402 !
403     END DO
404 !
405   ELSE  
406 !
407 ! Guesses at the end of the RK loop
408 !
409     DO JI = 1, ISPL
410      PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JI) * ZRUS(:,:,:,JI) 
411      PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JI) * ZRVS(:,:,:,JI) 
412      PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JI) * ZRWS(:,:,:,JI) 
413     END DO
414 !
415   END IF
416 !
417 ! End of the RK loop
418  END DO
419 !
420 !
421 DEALLOCATE(ZBUT, ZBUTS, ZRUS, ZRVS, ZRWS)
422 CALL CLEANLIST_ll(TZFIELDMT_ll)
423 CALL  DEL_HALO2_ll(TZHALO2MT_ll)
424 !-------------------------------------------------------------------------------
425 !
426 END SUBROUTINE ADVECUVW_RK