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