Philippe 17/06/2016: minor: added default case to abort if unknown HTEMP_SCHEME
[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  CASE DEFAULT
223   PRINT *,'ERROR: UNKNOWN HTEMP_SCHEME'
224   CALL ABORT()
225 END SELECT
226 !
227 !
228 ALLOCATE(ZBUT(ISPL-1,ISPL-1))
229 ALLOCATE(ZBUTS(ISPL))
230 !
231 SELECT CASE (HTEMP_SCHEME)
232   CASE('RK11')
233     ZBUTS = (/ 1. /)
234   CASE('RK21')
235     ZBUTS     = (/ 0. , 1. /)
236     ZBUT(1,1)   = 3./4.
237   CASE('RK33')
238     ZBUTS     = (/ 1./6. , 1./6. , 2./3. /)
239     ZBUT(1,1) = 1.
240     ZBUT(1,2) = 0.
241     ZBUT(2,1) = 1./4.
242     ZBUT(2,2) = 1./4.
243   CASE('NP32')
244     ZBUTS     = (/ 1./2. , 0., 1./2. /)
245     ZBUT(1,1) = 1./3.
246     ZBUT(1,2) = 0.
247     ZBUT(2,1) = 0.
248     ZBUT(2,2) = 1.
249   CASE('SP32')
250     ZBUTS     = (/ 1./3. , 1./3. , 1./3. /)
251     ZBUT(1,1) = 1./2.
252     ZBUT(1,2) = 0.
253     ZBUT(2,1) = 1./2.
254     ZBUT(2,2) = 1./2.
255   CASE('RKC4')
256     ZBUTS     = (/ 1./6. , 1./3. , 1./3. , 1./6./)
257     ZBUT      = 0.
258     ZBUT(1,1) = 1./2.
259     ZBUT(2,2) = 1./2.
260     ZBUT(3,3) = 1.
261   CASE('RK4B')
262     ZBUTS     = (/ 1./8. , 3./8. , 3./8. , 1./8./)
263     ZBUT      = 0.
264     ZBUT(1,1) = 1./3.
265     ZBUT(2,1) = -1./3.
266     ZBUT(2,2) = 1.
267     ZBUT(3,1) = 1.
268     ZBUT(3,2) = -1.
269     ZBUT(3,3) = 1.
270   CASE('RK53')
271     ZBUTS     = (/ 1./4. , 0. , 0. , 0. , 3./4. /)
272     ZBUT      = 0.
273     ZBUT(1,1) = 1./7.
274     ZBUT(2,2) = 3./16.
275     ZBUT(3,3) = 1./3.
276     ZBUT(4,4) = 2./3.
277   CASE('RK62')
278     ZBUTS     = (/ 1./6. , 1./6. , 1./6. , 1./6. , 1./6. , 1./6. /)
279     ZBUT      = 0.
280     ZBUT(1,1) = 1./5.
281     ZBUT(2,1) = 1./5.
282     ZBUT(2,2) = 1./5.
283     ZBUT(3,1) = 1./5.
284     ZBUT(3,2) = 1./5.
285     ZBUT(3,3) = 1./5.
286     ZBUT(4,1) = 1./5.
287     ZBUT(4,2) = 1./5.
288     ZBUT(4,3) = 1./5.
289     ZBUT(4,4) = 1./5.
290     ZBUT(5,1) = 1./5.
291     ZBUT(5,2) = 1./5.
292     ZBUT(5,3) = 1./5.
293     ZBUT(5,4) = 1./5.
294     ZBUT(5,5) = 1./5.
295 CASE('RK65')
296     ZBUTS= (/ 7./90. , 0. , 16./45. , 2./15. , 16./45. , 7./90. /)
297     ZBUT= 0.
298     ZBUT(1,1) = 1./4.
299     ZBUT(2,1) = 1./8.
300     ZBUT(2,2) = 1./8.
301     ZBUT(3,1) = 0
302     ZBUT(3,2) = -1./2.
303     ZBUT(3,3) = 1
304     ZBUT(4,1) = 3./16.
305     ZBUT(4,2) = 0
306     ZBUT(4,3) = 0
307     ZBUT(4,4) = 9./16.
308     ZBUT(5,1) = -3./7.
309     ZBUT(5,2) = 2./7.
310     ZBUT(5,3) = 12./7.
311     ZBUT(5,4) = -12./7.
312     ZBUT(5,5) = 8./7.
313 END SELECT
314 !
315 ALLOCATE(ZRUS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
316 ALLOCATE(ZRVS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
317 ALLOCATE(ZRWS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
318 !
319 PRUS_ADV = 0.
320 PRVS_ADV = 0.
321 PRWS_ADV = 0.
322 !
323 !-------------------------------------------------------------------------------
324 !
325 !*       2.     Wind guess before RK loop
326 !               -------------------------
327 !
328 ZUT = PU
329 ZVT = PV
330 ZWT = PW
331 !
332 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
333 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
334 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
335 ZWT (:,:,IKE+1 ) = 0.
336
337 ZU = PU
338 ZV = PV
339 ZW = PW
340 !
341 NULLIFY(TZFIELDMT_ll)    
342 CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZUT)
343 CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZVT)
344 CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZWT)
345 INBVAR = 3
346 CALL INIT_HALO2_ll(TZHALO2MT_ll,INBVAR,SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3))
347 !
348 ZRUS = 0.
349 ZRVS = 0.
350 ZRWS = 0.
351 !-------------------------------------------------------------------------------
352 !
353 !*       3.     BEGINNING of Runge-Kutta loop
354 !               -----------------------------
355 !
356  DO JS = 1, ISPL
357 !               
358     CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
359     CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
360     CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
361 !    
362     ZW (:,:,IKE+1 ) = 0.
363
364     CALL UPDATE_HALO_ll(TZFIELDMT_ll,IINFO_ll)        
365     CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll)
366 !
367 !*       4.     Advection with WENO
368 !               -------------------
369 !
370     CALL ADVECUVW_WENO_K (HLBCX, HLBCY, KWENO_ORDER, ZUT, ZVT, ZWT,     &
371                         PRUCT, PRVCT, PRWCT,                            &
372                         ZRUS(:,:,:,JS), ZRVS(:,:,:,JS), ZRWS(:,:,:,JS), &
373                         TZHALO2MT_ll                                    )
374 !                       
375      NULLIFY(TZFIELDS4_ll)
376 !
377     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRUS(:,:,:,JS))
378     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRVS(:,:,:,JS))
379     CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRWS(:,:,:,JS))
380     CALL UPDATE_HALO_ll(TZFIELDS4_ll,IINFO_ll)
381     CALL CLEANLIST_ll(TZFIELDS4_ll)
382 !               
383   IF ( JS /= ISPL ) THEN
384 !
385    DO JI = 1, JS
386
387       ZUT = ZU
388       ZVT = ZV
389       ZWT = ZW
390 !
391 ! Intermediate guesses inside the RK loop
392 !
393       ZUT(:,:,:) = ZUT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
394        ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ
395       ZVT(:,:,:) = ZVT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
396        ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ
397       ZWT(:,:,:) = ZWT(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
398        ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ
399 !
400     END DO
401 !
402   ELSE  
403 !
404 ! Guesses at the end of the RK loop
405 !
406     DO JI = 1, ISPL
407      PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JI) * ZRUS(:,:,:,JI) 
408      PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JI) * ZRVS(:,:,:,JI) 
409      PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JI) * ZRWS(:,:,:,JI) 
410     END DO
411 !
412   END IF
413 !
414 ! End of the RK loop
415  END DO
416 !
417 !
418 DEALLOCATE(ZBUT, ZBUTS, ZRUS, ZRVS, ZRWS)
419 CALL CLEANLIST_ll(TZFIELDMT_ll)
420 CALL  DEL_HALO2_ll(TZHALO2MT_ll)
421 !-------------------------------------------------------------------------------
422 !
423 END SUBROUTINE ADVECUVW_RK