Juan 9/01/2014: add header MNH_LIC to all MNH/*.f* files
[MNH-git_open_source-lfs.git] / src / MNH / advecuvw_rk.f90
1 !MNH_LIC Copyright 1994-2013 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 LICENCE, 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 for boundary
35                                                   ! conditions
36 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMXM_RHODJ
37 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMYM_RHODJ
38 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMZM_RHODJ
39                                                   !  metric coefficients
40 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUCT , PRVCT, PRWCT
41                                                   ! Contravariant wind components
42 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRUS_ADV , PRVS_ADV, PRWS_ADV
43                                                   ! Tendency due to advection
44 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER
45 !                                                 ! tendencies from other processes
46 !
47 END SUBROUTINE ADVECUVW_RK
48 !
49 END INTERFACE
50 !
51 END MODULE MODI_ADVECUVW_RK
52 !     ##########################################################################
53       SUBROUTINE ADVECUVW_RK (HUVW_ADV_SCHEME,                         &
54                     HTEMP_SCHEME, KWENO_ORDER,                         &
55                     HLBCX, HLBCY, PTSTEP,                              &
56                     PU, PV, PW,                                        &
57                     PUT, PVT, PWT,                                     &
58                     PMXM_RHODJ, PMYM_RHODJ, PMZM_RHODJ,                &
59                     PRUCT, PRVCT, PRWCT,                               &
60                     PRUS_ADV, PRVS_ADV, PRWS_ADV,                      &
61                     PRUS_OTHER, PRVS_OTHER, PRWS_OTHER                 )
62 !     ##########################################################################
63 !
64 !!****  *ADVECUVW_RK * - routine to call the specialized advection routines for wind
65 !!
66 !!    PURPOSE
67 !!    -------
68 !!
69 !!**  METHOD
70 !!    ------
71 !!
72 !!    EXTERNAL
73 !!    --------
74 !!
75 !!    IMPLICIT ARGUMENTS
76 !!    ------------------
77 !!      NONE
78 !!
79 !!    REFERENCE
80 !!    ---------
81 !!      Book1 and book2 ( routine ADVECTION )
82 !!
83 !!    AUTHOR
84 !!    ------
85 !!      J.-P. Pinty      * Laboratoire d'Aerologie*
86 !!      J.-P. Lafore     * Meteo France *
87 !!
88 !!    MODIFICATIONS
89 !!    -------------
90 !!      Original    06/07/94 
91 !!                  01/04/95 (Ph. Hereil J. Nicolau) add the model number
92 !!                  23/10/95 (J. Vila and JP Lafore) advection schemes scalar
93 !!                  16/01/97 (JP Pinty)              change presentation 
94 !!                  30/04/98 (J. Stein P Jabouille)  extrapolation for the cyclic
95 !!                                                   case and parallelisation
96 !!                  24/06/99 (P Jabouille)           case of NHALO>1
97 !!                  25/10/05 (JP Pinty)              4th order scheme
98 !!                  24/04/06 (C.Lac)                 Split scalar and passive
99 !!                                                   tracer routines
100 !!                  08/06    (T.Maric)               PPM scheme
101 !!                  04/2011  (V. Masson & C. Lac)    splits the routine and adds
102 !!                                                   time splitting
103 !!
104 !-------------------------------------------------------------------------------
105 !
106 !*       0.    DECLARATIONS
107 !              ------------
108 !
109 USE MODE_ll
110 USE MODD_ARGSLIST_ll, ONLY : LIST_ll, HALO2LIST_ll
111 USE MODD_PARAMETERS,  ONLY : JPVEXT
112 USE MODD_CONF,        ONLY : NHALO
113 !
114 USE MODI_SHUMAN
115 USE MODI_ADVECUVW_WENO_K
116 USE MODI_ADV_BOUNDARIES
117 USE MODI_GET_HALO
118 !
119 !-------------------------------------------------------------------------------
120 !
121 IMPLICIT NONE
122 !
123 !*       0.1   Declarations of dummy arguments :
124 !
125 CHARACTER(LEN=6),         INTENT(IN)    :: HUVW_ADV_SCHEME! to the selected
126 CHARACTER(LEN=4),         INTENT(IN)    :: HTEMP_SCHEME   ! Temporal scheme
127 !
128 INTEGER,                  INTENT(IN)    :: KWENO_ORDER   ! Order of the WENO
129                                                          ! scheme (3 or 5)
130 !
131 CHARACTER(LEN=4),DIMENSION(2),INTENT(IN):: HLBCX, HLBCY  ! X- and Y-direc LBC
132 !
133 REAL,                     INTENT(IN)    :: PTSTEP
134 !
135 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PU , PV  , PW
136                                                   ! Variables to advect
137 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PUT, PVT , PWT
138                                                   ! Variables for boundary
139                                                   ! conditions
140 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMXM_RHODJ
141 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMYM_RHODJ
142 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PMZM_RHODJ
143                                                   !  metric coefficients
144 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUCT , PRVCT, PRWCT
145                                                   ! Contravariant wind components
146 REAL, DIMENSION(:,:,:),   INTENT(OUT)   :: PRUS_ADV , PRVS_ADV, PRWS_ADV
147                                                   ! Tendency due to advection
148 REAL, DIMENSION(:,:,:),   INTENT(IN)    :: PRUS_OTHER , PRVS_OTHER, PRWS_OTHER
149 !                                                 ! tendencies from other processes
150 !
151 !
152 !
153 !*       0.2   declarations of local variables
154 !
155 !
156 !  
157 INTEGER             :: IKE       ! indice K End       in z direction
158 !
159 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZU, ZV, ZW
160 ! Guesses at the beginning of the RK loop
161 REAL, DIMENSION(SIZE(PUT,1),SIZE(PUT,2),SIZE(PUT,3)) :: ZUT, ZVT, ZWT
162 ! Intermediate Guesses inside the RK loop              
163 !
164 REAL, DIMENSION(:,:,:,:), ALLOCATABLE :: ZRUS,ZRVS,ZRWS
165 ! Momentum tendencies due to advection
166 REAL, DIMENSION(:,:), ALLOCATABLE :: ZBUT ! Butcher array coefficients
167                                           ! at the RK sub time step
168 REAL, DIMENSION(:),   ALLOCATABLE :: ZBUTS! Butcher array coefficients
169                                           ! at the end of the RK loop
170
171 !JUAN
172 TYPE(LIST_ll), POINTER      :: TZFIELDMT_ll ! list of fields to exchange
173 TYPE(HALO2LIST_ll), POINTER :: TZHALO2MT_ll ! momentum variables
174 INTEGER                     :: INBVAR
175 INTEGER :: IIU, IJU, IKU ! array sizes
176 !JUAN
177
178 ! Momentum tendencies due to advection
179 INTEGER :: ISPL                ! Number of RK splitting loops
180 INTEGER :: JI, JS              ! Loop index
181 !
182 INTEGER                     :: IINFO_ll    ! return code of parallel routine
183 TYPE(LIST_ll), POINTER      :: TZFIELD_ll  ! list of fields to exchange
184 TYPE(LIST_ll), POINTER      :: TZFIELDS_ll ! list of fields to exchange
185 TYPE(LIST_ll), POINTER      :: TZFIELDS0_ll ! list of fields to exchange
186 TYPE(LIST_ll), POINTER      :: TZFIELDS4_ll ! list of fields to exchange
187 !
188 !
189 REAL          :: XPRECISION
190 !-------------------------------------------------------------------------------
191 !
192 !*       0.     INITIALIZATION                        
193 !               --------------
194 !
195 IKE = SIZE(PWT,3) - JPVEXT
196 IIU=SIZE(PUT,1)
197 IJU=SIZE(PUT,2)
198 IKU=SIZE(PUT,3)
199 !
200 SELECT CASE (HTEMP_SCHEME)
201  CASE('RK11')
202   ISPL = 1
203  CASE('RK21')
204   ISPL = 2
205  CASE('RK33')
206   ISPL = 3
207  CASE('RK53')
208   ISPL = 5
209 END SELECT
210 !
211 !
212 ALLOCATE(ZBUT(ISPL-1,ISPL-1))
213 ALLOCATE(ZBUTS(ISPL))
214 !
215 IF (ISPL == 1 ) ZBUTS = (/ 1. /)
216 IF (ISPL == 2 ) THEN
217    ZBUTS     = (/ 0. , 1. /)
218    ZBUT(1,1)   = 3./4.
219 END IF
220 IF (ISPL == 3 ) THEN
221    ZBUTS     = (/ 1./6. , 1./6. , 2./3. /)
222    ZBUT(1,1) = 1.
223    ZBUT(1,2) = 0.
224    ZBUT(2,:) = 1./4.
225 END IF
226 IF (ISPL == 5 ) THEN
227    ZBUTS     = (/ 1./4. , 0., 0., 0., 3./4. /)
228    ZBUT      = 0.
229    ZBUT(1,1) = 1./7.
230    ZBUT(2,2) = 3./16.
231    ZBUT(3,3) = 1./3.
232    ZBUT(4,4) = 2./3.
233 END IF
234 !
235 ALLOCATE(ZRUS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
236 ALLOCATE(ZRVS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
237 ALLOCATE(ZRWS(SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3),ISPL))
238 !
239 PRUS_ADV = 0.
240 PRVS_ADV = 0.
241 PRWS_ADV = 0.
242 !
243 !-------------------------------------------------------------------------------
244 !
245 !*       2.     Wind guess before RK loop
246 !               -------------------------
247 !
248 ZUT = PU
249 ZVT = PV
250 ZWT = PW
251 !
252 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
253 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
254 CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
255 ZWT (:,:,IKE+1 ) = 0.
256
257 ZU = PU
258 ZV = PV
259 ZW = PW
260 !
261 NULLIFY(TZFIELDMT_ll)
262 IF( NHALO==1 ) THEN      
263 !
264    CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZUT)
265    CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZVT)
266    CALL ADD3DFIELD_ll(TZFIELDMT_ll, ZWT)
267 !
268    INBVAR = 3
269    IF( NHALO==1 ) CALL INIT_HALO2_ll(TZHALO2MT_ll,INBVAR,SIZE(PUT,1),SIZE(PUT,2),SIZE(PWT,3))
270 !
271  END IF
272 !
273 ZRUS = 0.
274 ZRVS = 0.
275 ZRWS = 0.
276 !-------------------------------------------------------------------------------
277 !
278 !*       3.     BEGINNING of Runge-Kutta loop
279 !               -----------------------------
280 !
281  DO JS = 1, ISPL
282 !
283 !
284       CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZUT, PUT, 'U' )    
285       CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZVT, PVT, 'V' )    
286       CALL ADV_BOUNDARIES (HLBCX, HLBCY, ZWT, PWT, 'W' )
287       ZW (:,:,IKE+1 ) = 0.
288      !JUAN
289      IF ( NHALO == 1 ) THEN   
290         CALL UPDATE_HALO_ll(TZFIELDMT_ll,IINFO_ll)        
291         CALL UPDATE_HALO2_ll(TZFIELDMT_ll, TZHALO2MT_ll, IINFO_ll)
292      ENDIF
293      !JUAN
294 !
295 !*       4.     Advection with WENO
296 !               -------------------
297 !
298      CALL ADVECUVW_WENO_K (HLBCX, HLBCY, KWENO_ORDER, ZUT, ZVT, ZWT,      &
299                           PRUCT, PRVCT, PRWCT,                            &
300                           ZRUS(:,:,:,JS), ZRVS(:,:,:,JS), ZRWS(:,:,:,JS), &
301                           TZHALO2MT_ll                                    )
302 !
303 !
304 ! ==> verifier si c'est utile !
305 !
306      NULLIFY(TZFIELDS4_ll)
307      IF(NHALO == 1) THEN
308           CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRUS(:,:,:,JS))
309           CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRVS(:,:,:,JS))
310           CALL ADD3DFIELD_ll(TZFIELDS4_ll, ZRWS(:,:,:,JS))
311           CALL UPDATE_HALO_ll(TZFIELDS4_ll,IINFO_ll)
312           CALL CLEANLIST_ll(TZFIELDS4_ll)
313      END IF
314
315     IF ( JS /= ISPL ) THEN
316 !
317
318      DO JI = 1, JS
319
320 !
321 ! Intermediate guesses inside the RK loop
322 !
323         ZUT(:,:,:) = ZU(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
324          ( ZRUS(:,:,:,JI) + PRUS_OTHER(:,:,:) ) / PMXM_RHODJ
325         ZVT(:,:,:) = ZV(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
326          ( ZRVS(:,:,:,JI) + PRVS_OTHER(:,:,:) ) / PMYM_RHODJ
327         ZWT(:,:,:) = ZW(:,:,:) + ZBUT(JS,JI) *  PTSTEP *  &
328          ( ZRWS(:,:,:,JI) + PRWS_OTHER(:,:,:) ) / PMZM_RHODJ
329 !
330       END DO
331 !
332     ELSE  
333 !
334 ! Guesses at the end of the RK loop
335 !
336       DO JI = 1, ISPL
337        PRUS_ADV(:,:,:) = PRUS_ADV(:,:,:) + ZBUTS(JI) * ZRUS(:,:,:,JI) 
338        PRVS_ADV(:,:,:) = PRVS_ADV(:,:,:) + ZBUTS(JI) * ZRVS(:,:,:,JI) 
339        PRWS_ADV(:,:,:) = PRWS_ADV(:,:,:) + ZBUTS(JI) * ZRWS(:,:,:,JI) 
340       END DO
341 !
342     END IF
343 !
344 ! End of the RK loop
345  END DO
346
347 !
348 !
349 DEALLOCATE(ZBUT, ZBUTS, ZRUS, ZRVS, ZRWS)
350 CALL CLEANLIST_ll(TZFIELDMT_ll)
351 CALL  DEL_HALO2_ll(TZHALO2MT_ll)
352 !-------------------------------------------------------------------------------
353 !
354 END SUBROUTINE ADVECUVW_RK