Merge branch 'LIBTOOLS-master' into MNH-52X
[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 !!
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 !!
105 !-------------------------------------------------------------------------------
106 !
107 !*       0.    DECLARATIONS
108 !              ------------
109 !
110 USE MODE_ll
111 USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll
112 USE MODD_PARAMETERS,  ONLY : JPVEXT
113 USE MODD_CONF,        ONLY : NHALO
114 !
115 USE MODI_SHUMAN
116 USE MODI_ADVECUVW_WENO_K
117 USE MODI_ADV_BOUNDARIES
118 USE MODI_GET_HALO
119 USE MODE_MPPDB
120 !
121 !-------------------------------------------------------------------------------
122 !
123 IMPLICIT NONE
124 !
125 !*       0.1   Declarations of dummy arguments :
126 !
127 CHARACTER(LEN=6),         INTENT(IN)    :: HUVW_ADV_SCHEME! to the selected
128 CHARACTER(LEN=4),         INTENT(IN)    :: HTEMP_SCHEME   ! Temporal scheme
129 !
130 INTEGER,                  INTENT(IN)    :: KWENO_ORDER   ! Order of the WENO
131                                                          ! scheme (3 or 5)
132 !
133 CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
134 !
135 REAL,                     INTENT(IN)    :: PTSTEP
136 !
137 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PU , PV  , PW
138                                                   ! Variables to advect
139 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT, PVT , PWT
140                                                   ! Variables at t
141 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMXM_RHODJ
142 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMYM_RHODJ
143 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMZM_RHODJ
144                                                   !  metric coefficients
145 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUCT , PRVCT, PRWCT
146                                                   ! Contravariant wind components
147 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRUS_ADV , PRVS_ADV, PRWS_ADV
148                                                   ! Tendency due to advection
149 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER
150 !                                                 ! tendencies from other processes
151 !
152 !
153 !
154 !*       0.2   declarations of local variables
155 !
156 !
157 !  
158 INTEGER             :: IKE       ! indice K End       in z direction
159 !
160 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZU, ZV, ZW
161 ! Guesses at the beginning of the RK loop
162 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZUT, ZVT, ZWT
163 ! Intermediate Guesses inside the RK loop              
164 !
165 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZRUS,ZRVS,ZRWS
166 ! Momentum tendencies due to advection
167 REAL, DIMENSION(:,:), ALLOCATABLE :: ZBUT ! Butcher array coefficients
168                                           ! at the RK sub time step
169 REAL, DIMENSION(:),   ALLOCATABLE :: ZBUTS! Butcher array coefficients
170                                           ! at the end of the RK loop
171
172 !JUAN
173 TYPE(LIST_ll), POINTER      :: TZFIELDMT_ll ! list of fields to exchange
174 TYPE(HALO2LIST_ll), POINTER :: TZHALO2MT_ll ! momentum variables
175 INTEGER                     :: INBVAR
176 INTEGER :: IIU, IJU, IKU ! array sizes
177 !JUAN
178
179 ! Momentum tendencies due to advection
180 INTEGER :: ISPL                ! Number of RK splitting loops
181 INTEGER :: JI, JS              ! Loop index
182 !
183 INTEGER                     :: IINFO_ll    ! return code of parallel routine
184 TYPE(LIST_ll), POINTER      :: TZFIELD_ll  ! list of fields to exchange
185 TYPE(LIST_ll), POINTER      :: TZFIELDS_ll ! list of fields to exchange
186 TYPE(LIST_ll), POINTER      :: TZFIELDS0_ll ! list of fields to exchange
187 TYPE(LIST_ll), POINTER      :: TZFIELDS4_ll ! list of fields to exchange
188 !
189 !
190 REAL          :: XPRECISION
191 !-------------------------------------------------------------------------------
192 !
193 !*       0.     INITIALIZATION                        
194 !               --------------
195 !
196 IKE = SIZE(PWT,3) - JPVEXT
197 IIU=SIZE(PUT,1)
198 IJU=SIZE(PUT,2)
199 IKU=SIZE(PUT,3)
200 !
201 SELECT CASE (HTEMP_SCHEME)
202  CASE('RK11')
203   ISPL = 1
204  CASE('RK21')
205   ISPL = 2
206  CASE('NP32')
207   ISPL = 3
208  CASE('SP32')
209   ISPL = 3
210  CASE('RK33')
211   ISPL = 3
212  CASE('RKC4')
213   ISPL = 4
214  CASE('RK4B')
215   ISPL = 4
216  CASE('RK53')
217   ISPL = 5
218  CASE('RK62')
219   ISPL = 6
220  CASE('RK65')
221   ISPL = 6
222 END SELECT
223 !
224 !
225 ALLOCATE(ZBUT(ISPL-1,ISPL-1))
226 ALLOCATE(ZBUTS(ISPL))
227 !
228 SELECT CASE (HTEMP_SCHEME)
229   CASE('RK11')
230     ZBUTS = (/ 1. /)
231   CASE('RK21')
232     ZBUTS     = (/ 0. , 1. /)
233     ZBUT(1,1)   = 3./4.
234   CASE('RK33')
235     ZBUTS     = (/ 1./6. , 1./6. , 2./3. /)
236     ZBUT(1,1) = 1.
237     ZBUT(1,2) = 0.
238     ZBUT(2,1) = 1./4.
239     ZBUT(2,2) = 1./4.
240   CASE('NP32')
241     ZBUTS     = (/ 1./2. , 0., 1./2. /)
242     ZBUT(1,1) = 1./3.
243     ZBUT(1,2) = 0.
244     ZBUT(2,1) = 0.
245     ZBUT(2,2) = 1.
246   CASE('SP32')
247     ZBUTS     = (/ 1./3. , 1./3. , 1./3. /)
248     ZBUT(1,1) = 1./2.
249     ZBUT(1,2) = 0.
250     ZBUT(2,1) = 1./2.
251     ZBUT(2,2) = 1./2.
252   CASE('RKC4')
253     ZBUTS     = (/ 1./6. , 1./3. , 1./3. , 1./6./)
254     ZBUT      = 0.
255     ZBUT(1,1) = 1./2.
256     ZBUT(2,2) = 1./2.
257     ZBUT(3,3) = 1.
258   CASE('RK4B')
259     ZBUTS     = (/ 1./8. , 3./8. , 3./8. , 1./8./)
260     ZBUT      = 0.
261     ZBUT(1,1) = 1./3.
262     ZBUT(2,1) = -1./3.
263     ZBUT(2,2) = 1.
264     ZBUT(3,1) = 1.
265     ZBUT(3,2) = -1.
266     ZBUT(3,3) = 1.
267   CASE('RK53')
268     ZBUTS     = (/ 1./4. , 0. , 0. , 0. , 3./4. /)
269     ZBUT      = 0.
270     ZBUT(1,1) = 1./7.
271     ZBUT(2,2) = 3./16.
272     ZBUT(3,3) = 1./3.
273     ZBUT(4,4) = 2./3.
274   CASE('RK62')
275     ZBUTS     = (/ 1./6. , 1./6. , 1./6. , 1./6. , 1./6. , 1./6. /)
276     ZBUT      = 0.
277     ZBUT(1,1) = 1./5.
278     ZBUT(2,1) = 1./5.
279     ZBUT(2,2) = 1./5.
280     ZBUT(3,1) = 1./5.
281     ZBUT(3,2) = 1./5.
282     ZBUT(3,3) = 1./5.
283     ZBUT(4,1) = 1./5.
284     ZBUT(4,2) = 1./5.
285     ZBUT(4,3) = 1./5.
286     ZBUT(4,4) = 1./5.
287     ZBUT(5,1) = 1./5.
288     ZBUT(5,2) = 1./5.
289     ZBUT(5,3) = 1./5.
290     ZBUT(5,4) = 1./5.
291     ZBUT(5,5) = 1./5.
292 CASE('RK65')
293     ZBUTS= (/ 7./90. , 0. , 16./45. , 2./15. , 16./45. , 7./90. /)
294     ZBUT= 0.
295     ZBUT(1,1) = 1./4.
296     ZBUT(2,1) = 1./8.
297     ZBUT(2,2) = 1./8.
298     ZBUT(3,1) = 0
299     ZBUT(3,2) = -1./2.
300     ZBUT(3,3) = 1
301     ZBUT(4,1) = 3./16.
302     ZBUT(4,2) = 0
303     ZBUT(4,3) = 0
304     ZBUT(4,4) = 9./16.
305     ZBUT(5,1) = -3./7.
306     ZBUT(5,2) = 2./7.
307     ZBUT(5,3) = 12./7.
308     ZBUT(5,4) = -12./7.
309     ZBUT(5,5) = 8./7.
310 END SELECT
311 !
312 ALLOCATE(ZRUS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
313 ALLOCATE(ZRVS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
314 ALLOCATE(ZRWS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
315 !
316 PRUS_ADV = 0.
317 PRVS_ADV = 0.
318 PRWS_ADV = 0.
319 !
320 !-------------------------------------------------------------------------------
321 !
322 !*       2.     Wind guess before RK loop
323 !               -------------------------
324 !
325 ZUT = PU
326 ZVT = PV
327 ZWT = PW
328 !
329 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
330 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
331 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
332 ZWT (:,:,IKE+1 ) = 0.
333
334 ZU = PU
335 ZV = PV
336 ZW = PW
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     ZW (:,:,IKE+1 ) = 0.
360
361     CALL UPDATE_HALO_ll(TZFIELDMT_ll,IINFO_ll)        
362     CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll)
363 !
364 !*       4.     Advection with WENO
365 !               -------------------
366 !
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 !                       
372      NULLIFY(TZFIELDS4_ll)
373 !
374     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRUS(:,:,:,JS))
375     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRVS(:,:,:,JS))
376     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRWS(:,:,:,JS))
377     CALL UPDATE_HALO_ll(TZFIELDS4_ll,IINFO_ll)
378     CALL CLEANLIST_ll(TZFIELDS4_ll)
379 !               
380   IF ( JS /= ISPL ) THEN
381 !
382    DO JI = 1, JS
383
384       ZUT = ZU
385       ZVT = ZV
386       ZWT = ZW
387 !
388 ! Intermediate guesses inside the RK loop
389 !
390       ZUT(:,:,:) = ZUT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
391        ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ
392       ZVT(:,:,:) = ZVT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
393        ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ
394       ZWT(:,:,:) = ZWT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
395        ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ
396 !
397     END DO
398 !
399   ELSE  
400 !
401 ! Guesses at the end of the RK loop
402 !
403     DO JI = 1, ISPL
404      PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JI) * ZRUS(:,:,:,JI) 
405      PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JI) * ZRVS(:,:,:,JI) 
406      PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JI) * ZRWS(:,:,:,JI) 
407     END DO
408 !
409   END IF
410 !
411 ! End of the RK loop
412  END DO
413 !
414 !
415 DEALLOCATE(ZBUT, ZBUTS, ZRUS, ZRVS, ZRWS)
416 CALL CLEANLIST_ll(TZFIELDMT_ll)
417 CALL  DEL_HALO2_ll(TZHALO2MT_ll)
418 !-------------------------------------------------------------------------------
419 !
420 END SUBROUTINE ADVECUVW_RK