Philippe 17/06/2016: change variable status INOUT->OUT
[MNH-git_open_source-lfs.git] / src / MNH / advec_weno_k_2_aux.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       MODULE MODI_ADVEC_WENO_K_2_AUX
7 !     ##############################
8 !
9 INTERFACE
10 !
11       SUBROUTINE ADVEC_WENO_K_2_UX(HLBCX,PSRC, PRUCT, PR, TPHALO2)
12 !
13 USE MODE_ll
14 USE MODD_LUNIT
15 USE MODD_CONF
16 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
17 !
18 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
19 !
20 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
21 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRUCT ! contrav. comp. on MASS GRID
22 !
23 ! output source term
24 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
25 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
26 !
27 END SUBROUTINE ADVEC_WENO_K_2_UX
28 !
29 !                    ----------------------------
30 !
31       SUBROUTINE ADVEC_WENO_K_2_MX(HLBCX,PSRC, PRUCT, PR, TPHALO2)
32 !
33 USE MODE_ll
34 USE MODD_LUNIT
35 USE MODD_CONF
36 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
37 !
38 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
39 !
40 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
41 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRUCT ! contrav. comp. on MASS GRID
42 !
43 ! output source term
44 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
45 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
46 !
47 END SUBROUTINE ADVEC_WENO_K_2_MX
48 !
49 !                     ---------------------------
50 !
51       SUBROUTINE ADVEC_WENO_K_2_VY(HLBCY,PSRC, PRVCT, PR, TPHALO2)
52 !
53 USE MODE_ll
54 USE MODD_LUNIT
55 USE MODD_CONF
56 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
57 !
58 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
59 !
60 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
61 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRVCT ! contrav. comp. on MASS GRID
62 !
63 !
64 ! output source term
65 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
66 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
67 !
68 END SUBROUTINE ADVEC_WENO_K_2_VY
69 !
70 !                  ------------------------------
71 !
72       SUBROUTINE ADVEC_WENO_K_2_MY(HLBCY,PSRC, PRVCT, PR, TPHALO2)
73 !
74 USE MODE_ll
75 USE MODD_LUNIT
76 USE MODD_CONF
77 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
78 !
79 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
80 !
81 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
82 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRVCT ! contrav. comp. on MASS GRID
83 !
84 ! output source term
85 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
86 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
87 !
88 END SUBROUTINE ADVEC_WENO_K_2_MY
89 !
90 !                     -------------------------------
91 !
92 FUNCTION WENO_K_2_WZ(PSRC, PRWCT) RESULT(PR)
93 !
94 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on W grid at t
95 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRWCT ! contrav. comp. on MASS GRID
96 !
97 ! output source term
98 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: PR
99 !
100 END FUNCTION WENO_K_2_WZ
101 !
102 !                      ------------------------------
103 !
104 FUNCTION WENO_K_2_MZ(PSRC, PRWCT) RESULT(PR)
105 !
106 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on MASS grid at t
107 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRWCT ! contrav. comp. on W grid
108 !
109 ! output source term
110 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: PR
111 !
112 END FUNCTION WENO_K_2_MZ
113 !
114 END INTERFACE
115 !
116 END MODULE MODI_ADVEC_WENO_K_2_AUX
117 !
118 !-----------------------------------------------------------------------------
119 !
120 !     ############################################################
121       SUBROUTINE ADVEC_WENO_K_2_UX(HLBCX,PSRC, PRUCT, PR, TPHALO2)
122 !     ############################################################
123 !!
124 !!**** Computes PRUCT * PUT. Upstream fluxes of U in X direction.  
125 !!     Input PUT is on U Grid 'ie' (i,j,k) based on UGRID reference
126 !!     Output PR is on mass Grid 'ie' (i+1/2,j,k) based on UGRID reference
127 !!              
128 !!    AUTHOR
129 !!    ------
130 !!    F. Visentin   *CNRS/LA*               
131 !!
132 !!    MODIFICATIONS
133 !!    -------------
134 !!      J.Escobar 21/03/2013: for HALOK comment all NHALO=1 test
135 !!
136 !-------------------------------------------------------------------------------
137 !
138 USE MODE_ll
139 USE MODD_LUNIT
140 USE MODD_CONF
141 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
142 USE MODI_GET_HALO
143 !
144 IMPLICIT NONE
145 !
146 !*       0.1   Declarations of dummy arguments :
147 !
148 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
149 !
150 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
151 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRUCT ! contrav. comp. on MASS GRID
152 !
153 ! output source term
154 !
155 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
156 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
157 !
158 !*       0.2   Declarations of local variables :
159 !
160 INTEGER :: IIB,IJB    ! Begining useful area in x,y,z directions
161 INTEGER :: IIE,IJE    ! End useful area in x,y,z directions
162 INTEGER:: IW,IE,IWF,IEF   ! Coordinate of third order diffusion area
163 !
164 INTEGER:: ILUOUT,IRESP   ! for prints
165 !
166 ! intermediate reconstruction fluxes for positive wind case
167 !
168 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFPOS1, ZFPOS2
169 !
170 ! intermediate reconstruction fluxes for negative wind case
171 ! we need only one since ZFNEG2 = ZFPOS2
172 !
173 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFNEG1, ZFNEG2
174 !
175 ! smoothness indicators for positive wind case
176 !
177 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBPOS1, ZBPOS2
178 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBNEG1, ZBNEG2
179 !
180 ! WENO weights
181 !
182 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMP1, ZOMP2
183 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMN1, ZOMN2
184 !
185 ! standard weights
186 !
187 REAL, PARAMETER :: ZGAMMA1 = 1./3.
188 REAL, PARAMETER :: ZGAMMA2 = 2./3.
189 !
190 REAL, PARAMETER :: ZEPS = 1.0E-15
191 !
192 !-----------------------------------------------------------------------------
193 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
194 !                 ------------------------------
195 !
196 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
197 !
198 !-------------------------------------------------------------------------------
199 !
200 !*       0.4.   INITIALIZE THE FIELD 
201 !               ---------------------
202 !
203 PR(:,:,:) = 0.0
204 !
205 ZFPOS1  = 0.0
206 ZFPOS2  = 0.0
207 ZFNEG1  = 0.0
208 ZFNEG2  = 0.0
209 ZBPOS1  = 0.0
210 ZBPOS2  = 0.0
211 ZBNEG1  = 0.0
212 ZBNEG2  = 0.0
213 ZOMP1   = 0.0
214 ZOMP2   = 0.0
215 ZOMN1   = 0.0
216 ZOMN2   = 0.0
217 !
218 !-------------------------------------------------------------------------------
219 !
220 SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
221 !
222 !*       1.1    CYCLIC CASE IN THE X DIRECTION:
223 !
224 CASE ('CYCL')          ! In that case one must have HLBCX(1) == HLBCX(2)
225 !
226 !!$  IF(NHALO == 1) THEN
227     IW=IIB
228     IE=IIE
229 !!$  ELSE
230 !!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
231 !!$    WRITE(ILUOUT,*) 'ERROR : 3rd order advection in CYCLic case '
232 !!$    WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
233 !!$    CALL ABORT
234 !!$    STOP
235 !!$  END IF
236 !
237 ! r: many left cells in regard to 'i' cell for each stencil
238 !
239 ! intermediate fluxes at the mass point on Ugrid u(i+1/2,j,k) for positive wind
240 ! (r=1 for the first stencil ZFPOS1, r=0 for the second ZFPOS2)
241 !
242    ZFPOS1(IW:IE+1,:,:) = 0.5 * (3.0*PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:))
243    ZFPOS1(IW-1,   :,:) = 0.5 * (3.0*PSRC(IW-1,   :,:) - TPHALO2%WEST(:,:))
244 !
245    ZFPOS2(IW-1:IE,:,:) = 0.5 * (PSRC(IW-1:IE,:,:) + PSRC(IW:IE+1,:,:))
246    ZFPOS2(IE+1,   :,:) = 0.5 * (PSRC(IE+1,   :,:) + TPHALO2%EAST(:,:))
247 !
248 ! intermediate flux at the mass point on Ugrid (i+1/2,j,k) for negative wind 
249 ! case (from the right to the left)
250 ! (r=0 for the second stencil ZFNEG2=ZFPOS2, r=-1 for the first ZFNEG1)  
251 !
252   ZFNEG1(IW-1:IE-1,:,:) = 0.5 * (3.0*PSRC(IW:IE,:,:) - PSRC(IW+1:IE+1,:,:))
253   ZFNEG1(IE,   :,:) = 0.5 * (3.0*PSRC(IE+1,   :,:) - TPHALO2%EAST(:,:))
254   ZFNEG2(IW-1:IE,:,:) = 0.5 * (PSRC(IW-1:IE,:,:) + PSRC(IW:IE+1,:,:))
255   ZFNEG2(IE+1,   :,:) = 0.5 * (PSRC(IE+1,   :,:) + TPHALO2%EAST(:,:))
256 !
257 ! smoothness indicators for positive wind case
258 !
259   ZBPOS1(IW:IE+1,:,:) = (PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:))**2
260   ZBPOS1(IW-1,   :,:) = (PSRC(IW-1,   :,:) - TPHALO2%WEST(:,:))**2
261 !
262   ZBPOS2(IW-1:IE,:,:) = (PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:))**2
263   ZBPOS2(IE+1,   :,:) = (TPHALO2%EAST(:,:) - PSRC(IE+1,   :,:))**2
264 !
265 ! smoothness indicators for negative wind case
266 !       
267   ZBNEG1(IW-1:IE-1,:,:) = (PSRC(IW:IE,:,:)   - PSRC(IW+1:IE+1,:,:))**2
268   ZBNEG1(IE,   :,:)     = (PSRC(IE+1,   :,:) - TPHALO2%EAST(:,:))**2
269   ZBNEG2(IW-1:IE,:,:)   = (PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:))**2
270   ZBNEG2(IE+1,   :,:)   = (PSRC(IE+1,   :,:) - TPHALO2%EAST(:,:))**2
271 !
272 ! WENO weights
273 !
274   ZOMP1 = ZGAMMA1 / (ZEPS + ZBPOS1)**2
275   ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2
276   ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2
277   ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2
278 !
279 ! WENO fluxes
280 !
281   PR = (ZOMN2/(ZOMN1+ZOMN2) * ZFNEG2 +                           &
282        (ZOMN1/(ZOMN1+ZOMN2) * ZFNEG1)) * (0.5-SIGN(0.5,PRUCT)) + &
283        (ZOMP2/(ZOMP1+ZOMP2) * ZFPOS2 +                           &
284        (ZOMP1/(ZOMP1+ZOMP2) * ZFPOS1)) * (0.5+SIGN(0.5,PRUCT))
285 !
286 !
287 !       OPEN, WALL, NEST CASE IN THE X DIRECTION
288 !
289 CASE ('OPEN','WALL','NEST')
290 !
291   IW=IIB
292   IE=IIE
293 !
294 !       USE A FIRST ORDER UPSTREAM SCHEME AT THE PHYSICAL BORDER
295 !
296   IF(LWEST_ll()) THEN
297     PR(IW-1,:,:) = PSRC(IW-1,:,:) * (0.5+SIGN(0.5,PRUCT(IW-1,:,:))) + PSRC(IW,:,:) * (0.5-SIGN(0.5,PRUCT(IW-1,:,:)))
298 !
299 !!$  ELSEIF (NHALO == 1) THEN
300   ELSE 
301     ZFPOS1(IW-1,:,:) = 0.5 * (3.0*PSRC(IW-1,:,:) - TPHALO2%WEST(:,:))
302     ZFPOS2(IW-1,:,:) = 0.5 * (PSRC(IW-1,    :,:) + PSRC(IW,:,:))
303     ZBPOS1(IW-1,:,:) = (PSRC(IW-1,:,:) - TPHALO2%WEST(:,:))**2
304     ZBPOS2(IW-1,:,:) = (PSRC(IW,  :,:) - PSRC(IW-1,:,:))**2
305 !
306     ZFNEG1(IW-1,:,:) = 0.5 * (3.0*PSRC(IW,:,:) - PSRC(IW+1,:,:))
307     ZFNEG2(IW-1,:,:) = 0.5 * (PSRC(IW-1,  :,:) + PSRC(IW,  :,:))
308     ZBNEG1(IW-1,:,:) = (PSRC(IW,  :,:) - PSRC(IW+1,:,:))**2
309     ZBNEG2(IW-1,:,:) = (PSRC(IW-1,:,:) - PSRC(IW,  :,:))**2
310 !
311     ZOMP1(IW-1,:,:) = ZGAMMA1 / (ZEPS + ZBPOS1(IW-1,:,:))**2
312     ZOMP2(IW-1,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW-1,:,:))**2
313     ZOMN1(IW-1,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW-1,:,:))**2
314     ZOMN2(IW-1,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW-1,:,:))**2
315 !
316     PR(IW-1,:,:) = (ZOMN2(IW-1,:,:)/(ZOMN1(IW-1,:,:)+ZOMN2(IW-1,:,:)) * ZFNEG2(IW-1,:,:) +   &
317                    (ZOMN1(IW-1,:,:)/(ZOMN1(IW-1,:,:)+ZOMN2(IW-1,:,:)) * ZFNEG1(IW-1,:,:))) * (0.5-SIGN(0.5,PRUCT(IW-1,:,:))) + &
318                    (ZOMP2(IW-1,:,:)/(ZOMP1(IW-1,:,:)+ZOMP2(IW-1,:,:)) * ZFPOS2(IW-1,:,:) +  &
319                    (ZOMP1(IW-1,:,:)/(ZOMP1(IW-1,:,:)+ZOMP2(IW-1,:,:)) * ZFPOS1(IW-1,:,:))) * (0.5+SIGN(0.5,PRUCT(IW-1,:,:)))
320 !
321   ENDIF
322 !
323   IF(LEAST_ll()) THEN
324     PR(IE,:,:) = PSRC(IE,:,:) * (0.5+SIGN(0.5,PRUCT(IE,:,:))) + PSRC(IE+1,:,:) * (0.5-SIGN(0.5,PRUCT(IE,:,:)))
325 !
326 !!$  ELSEIF (NHALO == 1) THEN
327   ELSE
328     ZFPOS1(IE,:,:) = 0.5 * (3.0*PSRC(IE,:,:) - PSRC(IE-1,:,:))
329     ZFPOS2(IE,:,:) = 0.5 * (PSRC(IE,    :,:) + PSRC(IE+1,:,:)) 
330     ZBPOS1(IE,:,:) = (PSRC(IE,:,:) - PSRC(IE-1,:,:))**2
331     ZBPOS2(IE,:,:) = (PSRC(IE+1,:,:) - PSRC(IE,:,:))**2
332 !
333     ZFNEG1(IE,:,:) = 0.5 * (3.0*PSRC(IE+1,:,:) - TPHALO2%EAST(:,:))
334     ZFNEG2(IE,:,:) = 0.5 * (PSRC(IE,:,:) + PSRC(IE+1,:,:))
335     ZBNEG1(IE,:,:) = (PSRC(IE+1,:,:) - TPHALO2%EAST(:,:))**2
336     ZBNEG2(IE,:,:) = (PSRC(IE,  :,:) - PSRC(IE+1,:,:))**2
337 !
338     ZOMP1(IE,:,:) = ZGAMMA1 / (ZEPS + ZBPOS1(IE,:,:))**2
339     ZOMP2(IE,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IE,:,:))**2
340     ZOMN1(IE,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IE,:,:))**2
341     ZOMN2(IE,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IE,:,:))**2
342 !
343     PR(IE,:,:) = (ZOMN2(IE,:,:)/(ZOMN1(IE,:,:)+ZOMN2(IE,:,:)) * ZFNEG2(IE,:,:) +  &
344        (ZOMN1(IE,:,:)/(ZOMN1(IE,:,:)+ZOMN2(IE,:,:)) * ZFNEG1(IE,:,:))) * (0.5-SIGN(0.5,PRUCT(IE,:,:))) + &
345                  (ZOMP2(IE,:,:)/(ZOMP1(IE,:,:)+ZOMP2(IE,:,:)) * ZFPOS2(IE,:,:) +  &
346        (ZOMP1(IE,:,:)/(ZOMP1(IE,:,:)+ZOMP2(IE,:,:)) * ZFPOS1(IE,:,:))) * (0.5+SIGN(0.5,PRUCT(IE,:,:)))
347 !
348   ENDIF
349 !
350 !      USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE 
351 !
352   ZFPOS1(IW:IE-1,:,:) = 0.5 * (3.0*PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:))
353   ZFPOS2(IW:IE-1,:,:) = 0.5 * (PSRC(IW:IE-1,    :,:) + PSRC(IW+1:IE,  :,:))
354   ZBPOS1(IW:IE-1,:,:) = (PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:))**2
355   ZBPOS2(IW:IE-1,:,:) = (PSRC(IW+1:IE,:,:) - PSRC(IW:IE-1,  :,:))**2
356 !
357   ZFNEG1(IW:IE-1,:,:) = 0.5 * (3.0*PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:))
358   ZFNEG2(IW:IE-1,:,:) = 0.5 * (PSRC(IW:IE-1,    :,:) + PSRC(IW+1:IE,  :,:))
359   ZBNEG1(IW:IE-1,:,:) = (PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:))**2
360   ZBNEG2(IW:IE-1,:,:)   = (PSRC(IW:IE-1,:,:) - PSRC(IW+1:IE,:,:))**2
361 !
362   ZOMP1(IW:IE-1,:,:) = ZGAMMA1 / (ZEPS + ZBPOS1(IW:IE-1,:,:))**2
363   ZOMP2(IW:IE-1,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW:IE-1,:,:))**2
364   ZOMN1(IW:IE-1,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW:IE-1,:,:))**2
365   ZOMN2(IW:IE-1,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW:IE-1,:,:))**2
366 !
367     PR(IW:IE-1,:,:) = (ZOMN2(IW:IE-1,:,:)/(ZOMN1(IW:IE-1,:,:)+ZOMN2(IW:IE-1,:,:)) * ZFNEG2(IW:IE-1,:,:) +   &
368        (ZOMN1(IW:IE-1,:,:)/(ZOMN1(IW:IE-1,:,:)+ZOMN2(IW:IE-1,:,:)) * ZFNEG1(IW:IE-1,:,:))) * (0.5-SIGN(0.5,PRUCT(IW:IE-1,:,:))) + &
369        (ZOMP2(IW:IE-1,:,:)/(ZOMP1(IW:IE-1,:,:)+ZOMP2(IW:IE-1,:,:)) * ZFPOS2(IW:IE-1,:,:) +   &
370        (ZOMP1(IW:IE-1,:,:)/(ZOMP1(IW:IE-1,:,:)+ZOMP2(IW:IE-1,:,:)) * ZFPOS1(IW:IE-1,:,:))) * (0.5+SIGN(0.5,PRUCT(IW:IE-1,:,:)))
371 !
372 END SELECT
373 !
374 PR = PR * PRUCT
375 CALL GET_HALO(PR)
376 !
377 END SUBROUTINE ADVEC_WENO_K_2_UX
378 !
379 !------------------------------------------------------------------------------
380 !
381 !     ############################################################
382       SUBROUTINE ADVEC_WENO_K_2_MX(HLBCX,PSRC, PRUCT, PR, TPHALO2)
383 !     ############################################################
384 !!
385 !!**** Computes PRUCT * PWT (or PRUCT * PVT). Upstream fluxes of W (or V) 
386 !!     variables in X direction.  
387 !!     Input PWT is on W Grid 'ie' (i,j,k) based on WGRID reference
388 !!     Output PR is on mass Grid 'ie' (i-1/2,j,k) based on WGRID reference  
389 !!
390 !!    AUTHOR
391 !!    ------
392 !!    F. Visentin   *CNRS/LA*                
393 !!
394 !!    MODIFICATIONS
395 !!    -------------
396 !!
397 !------------------------------------------------------------------------------
398 !
399 USE MODE_ll
400 USE MODD_LUNIT
401 USE MODD_CONF
402 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
403 USE MODI_GET_HALO
404 !
405 IMPLICIT NONE
406 !
407 !*       0.1   Declarations of dummy arguments :
408 !
409 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCX  ! X direction LBC type
410 !
411 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
412 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRUCT ! contrav. comp. on MASS GRID
413 !
414 ! output source term
415 !
416 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
417 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
418 !
419 !*       0.2   Declarations of local variables :
420 !
421 INTEGER :: IIB,IJB    ! Begining useful area in x,y,z directions
422 INTEGER :: IIE,IJE    ! End useful area in x,y,z directions
423 INTEGER::  IW,IE   ! Coordinate of third order diffusion area
424 !
425 INTEGER:: ILUOUT,IRESP   ! for prints
426 !
427 ! intermediate reconstruction fluxes for positive wind case
428 !
429 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFPOS1, ZFPOS2
430 !
431 ! intermediate reconstruction fluxes for negative wind case
432 ! we need only one since ZFNEG2 = ZFPOS2
433 !
434 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFNEG1, ZFNEG2
435 !
436 ! smoothness indicators for positive wind case
437 !
438 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBPOS1, ZBPOS2
439 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBNEG1, ZBNEG2
440 !
441 ! WENO weights
442 !
443 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMP1, ZOMP2
444 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMN1, ZOMN2
445 !
446 ! standard weights
447 !
448 REAL, PARAMETER :: ZGAMMA1 = 1./3.
449 REAL, PARAMETER :: ZGAMMA2 = 2./3.
450 !
451 REAL, PARAMETER :: ZEPS = 1.0E-15
452 !
453 !-------------------------------------------------------------------------------
454 !
455 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
456 !                 ------------------------------
457 !
458 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
459 !
460 !-----------------------------------------------------------------------------
461 !
462 !*       0.4.   INITIALIZE THE FIELD 
463 !               ---------------------
464 !
465 PR(:,:,:) = 0.0
466 !
467 ZFPOS1 = 0.0
468 ZFPOS2 = 0.0
469 ZFNEG1 = 0.0
470 ZFNEG2 = 0.0 
471 ZBPOS1 = 0.0
472 ZBPOS2 = 0.0
473 ZBNEG1 = 0.0
474 ZBNEG2 = 0.0
475 ZOMP1  = 0.0
476 ZOMP2  = 0.0
477 ZOMN1  = 0.0
478 ZOMN2  = 0.0 
479 !
480 !------------------------------------------------------------------------------
481 !
482 SELECT CASE ( HLBCX(1) ) ! X direction LBC type: (1) for left side
483 !
484 !*       1.1    CYCLIC CASE IN THE X DIRECTION:
485 !
486 CASE ('CYCL')          ! In that case one must have HLBCX(1) == HLBCX(2)
487 !
488 !!$  IF(NHALO == 1) THEN
489     IW=IIB
490     IE=IIE
491 !!$  ELSE
492 !!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
493 !!$    WRITE(ILUOUT,*) 'ERROR : 3rd order advection in CYCLic case '
494 !!$    WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
495 !!$    CALL ABORT
496 !!$    STOP
497 !!$  END IF  
498 !
499 ! intermediate fluxes for positive wind case
500 !
501   ZFPOS1(IW+1:IE+1,:,:) = 0.5 * (3.0*PSRC(IW:IE,:,:) - PSRC(IW-1:IE-1,:,:))
502   ZFPOS1(IW,       :,:) = 0.5 * (3.0*PSRC(IW-1, :,:) - TPHALO2%WEST(:,:))
503 !!  ZFPOS1(IW-1,     :,:) = - 999.
504 !
505   ZFPOS2(IW:IE+1,:,:) = 0.5 * (PSRC(IW-1:IE,:,:) + PSRC(IW:IE+1,:,:))
506   ZFPOS2(IW-1,   :,:) = 0.5 * (TPHALO2%WEST(:,:) + PSRC(IW-1,   :,:))
507 !
508 ! intermediate flux for negative wind case
509 !
510   ZFNEG1(IW-1:IE,:,:) = 0.5 * (3.0*PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:))
511   ZFNEG1(IE+1,   :,:) = 0.5 * (3.0*PSRC(IE+1,   :,:) - TPHALO2%EAST(:,:))
512 !
513   ZFNEG2(IW:IE+1,:,:) = 0.5 * (PSRC(IW:IE+1,:,:) + PSRC(IW-1:IE,:,:))
514   ZFNEG2(IW-1,       :,:) = 0.5 * (PSRC(IW-1, :,:) + TPHALO2%WEST(:,:))
515
516 ! smoothness indicators for positive wind case
517 !
518   ZBPOS1(IW+1:IE+1,:,:) = (PSRC(IW:IE,:,:) - PSRC(IW-1:IE-1,:,:))**2
519   ZBPOS1(IW,       :,:) = (PSRC(IW-1, :,:) - TPHALO2%WEST(:,:))**2
520 !!  ZBPOS1(IW-1,     :,:) = - 999.
521 !
522   ZBPOS2(IW:IE+1,:,:) = (PSRC(IW:IE+1,:,:) - PSRC(IW-1:IE,:,:))**2
523   ZBPOS2(IW-1,   :,:) = (PSRC(IW-1,   :,:) - TPHALO2%WEST(:,:))**2
524 !
525 ! smoothness indicators for negative wind case
526 !       
527   ZBNEG1(IW-1:IE,:,:) = (PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:))**2
528   ZBNEG1(IE+1,   :,:) = (PSRC(IE+1,   :,:) - TPHALO2%EAST(:,:))**2
529 !
530   ZBNEG2(IW:IE+1,:,:) = (PSRC(IW-1:IE,:,:) - PSRC(IW:IE+1,:,:))**2
531   ZBNEG2(IW-1,   :,:) = (TPHALO2%WEST(:,:) - PSRC(IW-1,:,:))**2
532 !
533 ! WENO weights
534 !
535   ZOMP1 = ZGAMMA1 / (ZEPS + ZBPOS1)**2
536   ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2
537   ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2
538   ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2
539 !
540 ! WENO fluxes
541 !
542   PR = (ZOMP2/(ZOMP1+ZOMP2) * ZFPOS2 +                            &
543        (ZOMP1/(ZOMP1+ZOMP2) * ZFPOS1)) * (0.5+SIGN(0.5,PRUCT )) + &
544        (ZOMN2/(ZOMN1+ZOMN2) * ZFNEG2 +                            &
545        (ZOMN1/(ZOMN1+ZOMN2) * ZFNEG1)) * (0.5-SIGN(0.5,PRUCT ))
546 !
547 !
548 !       OPEN, WALL, NEST CASE IN THE X DIRECTION
549 !
550 CASE ('OPEN','WALL','NEST')
551 !
552   IW=IIB
553   IE=IIE
554 !
555 !       USE A FIRST ORDER UPSTREAM SCHEME AT THE PHYSICAL BORDER
556 !
557   IF(LWEST_ll()) THEN
558     PR(IW,:,:) = PSRC(IW-1,:,:) * (0.5+SIGN(0.5,PRUCT(IW,:,:))) + PSRC(IW,:,:) * (0.5-SIGN(0.5,PRUCT(IW,:,:)))
559 !
560 !!$  ELSEIF (NHALO == 1) THEN
561   ELSE
562     ZFPOS1(IW,:,:) = 0.5 * (3.0*PSRC(IW-1, :,:) - TPHALO2%WEST(:,:))
563     ZFPOS2(IW,:,:) = 0.5 * (PSRC(IW-1,     :,:) + PSRC(IW,     :,:))
564     ZBPOS1(IW,:,:) = (PSRC(IW-1,:,:) - TPHALO2%WEST(:,:))**2
565     ZBPOS2(IW,:,:) = (PSRC(IW,  :,:) - PSRC(IW-1,:,:))**2
566 !
567     ZFNEG1(IW,:,:) = 0.5 * (3.0*PSRC(IW,:,:) - PSRC(IW+1,:,:))
568     ZFNEG2(IW,:,:) = 0.5 * (PSRC(IW,    :,:) + PSRC(IW-1,:,:))
569     ZBNEG1(IW,:,:) = (PSRC(IW,:,:) - PSRC(IW+1,:,:))**2
570     ZBNEG2(IW,:,:) = (PSRC(IW-1,:,:) - PSRC(IW,:,:))**2
571 !
572     ZOMP1(IW,:,:) = ZGAMMA1 / (ZEPS + ZBPOS1(IW,:,:))**2
573     ZOMP2(IW,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW,:,:))**2
574     ZOMN1(IW,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW,:,:))**2
575     ZOMN2(IW,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW,:,:))**2
576 !
577     PR(IW,:,:) = (ZOMP2(IW,:,:)/(ZOMP1(IW,:,:)+ZOMP2(IW,:,:)) * ZFPOS2(IW,:,:) +  &
578        (ZOMP1(IW,:,:)/(ZOMP1(IW,:,:)+ZOMP2(IW,:,:)) * ZFPOS1(IW,:,:))) * (0.5+SIGN(0.5,PRUCT(IW,:,:))) + &
579        (ZOMN2(IW,:,:)/(ZOMN1(IW,:,:)+ZOMN2(IW,:,:)) * ZFNEG2(IW,:,:) +   &
580        (ZOMN1(IW,:,:)/(ZOMN1(IW,:,:)+ZOMN2(IW,:,:)) * ZFNEG1(IW,:,:))) * (0.5-SIGN(0.5,PRUCT(IW,:,:)))
581 !
582   ENDIF
583 !
584   IF(LEAST_ll()) THEN
585     PR(IE+1,:,:) = PSRC(IE,:,:) * (0.5+SIGN(0.5,PRUCT(IE+1,:,:))) + PSRC(IE+1,:,:) * (0.5-SIGN(0.5,PRUCT(IE+1,:,:)))
586 !
587 !!$  ELSEIF (NHALO == 1) THEN
588   ELSE
589     ZFPOS1(IE+1,:,:) = 0.5 * (3.0*PSRC(IE,:,:) - PSRC(IE-1,:,:))
590     ZFPOS2(IE+1,:,:) = 0.5 * (PSRC(IE,    :,:) + PSRC(IE+1,:,:))
591     ZBPOS1(IE+1,:,:) = (PSRC(IE,:,:) - PSRC(IE-1,:,:))**2
592     ZBPOS2(IE+1,:,:) = (PSRC(IE+1,:,:) - PSRC(IE,:,:))**2
593 !
594     ZFNEG1(IE+1,:,:) = 0.5 * (3.0*PSRC(IE+1,:,:) - TPHALO2%EAST(:,:))
595     ZFNEG2(IE+1,:,:) = 0.5 * (PSRC(IE+1,    :,:) + PSRC(IE,:,:))
596     ZBNEG1(IE+1,:,:) = (PSRC(IE+1,:,:) - TPHALO2%EAST(:,:))**2
597     ZBNEG2(IE+1,:,:) = (PSRC(IE,  :,:) - PSRC(IE+1,:,:))**2
598 !
599     ZOMP1(IE+1,:,:) = ZGAMMA1 / (ZEPS + ZBPOS1(IE+1,:,:))**2
600     ZOMP2(IE+1,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IE+1,:,:))**2
601     ZOMN1(IE+1,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IE+1,:,:))**2
602     ZOMN2(IE+1,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IE+1,:,:))**2
603 !
604     PR(IE+1,:,:) = (ZOMP2(IE+1,:,:)/(ZOMP1(IE+1,:,:)+ZOMP2(IE+1,:,:)) * ZFPOS2(IE+1,:,:) +  &
605        (ZOMP1(IE+1,:,:)/(ZOMP1(IE+1,:,:)+ZOMP2(IE+1,:,:)) * ZFPOS1(IE+1,:,:))) * (0.5+SIGN(0.5,PRUCT(IE+1,:,:))) + &
606        (ZOMN2(IE+1,:,:)/(ZOMN1(IE+1,:,:)+ZOMN2(IE+1,:,:)) * ZFNEG2(IE+1,:,:) +   &
607        (ZOMN1(IE+1,:,:)/(ZOMN1(IE+1,:,:)+ZOMN2(IE+1,:,:)) * ZFNEG1(IE+1,:,:))) * (0.5-SIGN(0.5,PRUCT(IE+1,:,:)))
608 !
609   ENDIF
610 !
611 !      USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE 
612 !
613   ZFPOS1(IW+1:IE,:,:) = 0.5 * (3.0*PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:))
614   ZFPOS2(IW+1:IE,:,:) = 0.5 * (PSRC(IW:IE-1,    :,:) + PSRC(IW+1:IE,  :,:))
615   ZBPOS1(IW+1:IE,:,:) = (PSRC(IW:IE-1,:,:) - PSRC(IW-1:IE-2,:,:))**2
616   ZBPOS2(IW+1:IE,:,:) = (PSRC(IW+1:IE,:,:) - PSRC(IW:IE-1,:,:))**2
617 !
618   ZFNEG1(IW+1:IE,:,:) = 0.5 * (3.0*PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:))
619   ZFNEG2(IW+1:IE,:,:) = 0.5 * (PSRC(IW+1:IE,    :,:) + PSRC(IW:IE-1,  :,:))
620   ZBNEG1(IW+1:IE,:,:) = (PSRC(IW+1:IE,:,:) - PSRC(IW+2:IE+1,:,:))**2
621   ZBNEG2(IW+1:IE,:,:) = (PSRC(IW:IE-1,:,:) - PSRC(IW+1:IE,:,:))**2
622 !
623   ZOMP1(IW+1:IE,:,:) = ZGAMMA1 / (ZEPS + ZBPOS1(IW+1:IE,:,:))**2
624   ZOMP2(IW+1:IE,:,:) = ZGAMMA2 / (ZEPS + ZBPOS2(IW+1:IE,:,:))**2
625   ZOMN1(IW+1:IE,:,:) = ZGAMMA1 / (ZEPS + ZBNEG1(IW+1:IE,:,:))**2
626   ZOMN2(IW+1:IE,:,:) = ZGAMMA2 / (ZEPS + ZBNEG2(IW+1:IE,:,:))**2
627 !
628   PR(IW+1:IE,:,:) = (ZOMP2(IW+1:IE,:,:)/(ZOMP1(IW+1:IE,:,:)+ZOMP2(IW+1:IE,:,:)) * ZFPOS2(IW+1:IE,:,:) + &
629        (ZOMP1(IW+1:IE,:,:)/(ZOMP1(IW+1:IE,:,:)+ZOMP2(IW+1:IE,:,:)) * ZFPOS1(IW+1:IE,:,:))) * (0.5+SIGN(0.5,PRUCT(IW+1:IE,:,:))) + &
630        (ZOMN2(IW+1:IE,:,:)/(ZOMN1(IW+1:IE,:,:)+ZOMN2(IW+1:IE,:,:)) * ZFNEG2(IW+1:IE,:,:) +              &
631        (ZOMN1(IW+1:IE,:,:)/(ZOMN1(IW+1:IE,:,:)+ZOMN2(IW+1:IE,:,:)) * ZFNEG1(IW+1:IE,:,:))) * (0.5-SIGN(0.5,PRUCT(IW+1:IE,:,:)))
632 !
633 END SELECT
634 !
635 PR = PR * PRUCT
636 CALL GET_HALO(PR)
637 !
638 END SUBROUTINE ADVEC_WENO_K_2_MX
639 !
640 !-------------------------------------------------------------------------------
641 !
642 !     ############################################################
643       SUBROUTINE ADVEC_WENO_K_2_MY(HLBCY,PSRC, PRVCT, PR, TPHALO2)
644 !     ############################################################
645 !!
646 !!****  Computes PRVCT * PUT (or PRVCT * PWT). Upstream fluxes of U (or W) 
647 !!      variables in Y direction.  
648 !!      Input PUT is on U Grid 'ie' (i,j,k) based on UGRID reference
649 !!      Output PR is on mass Grid 'ie' (i,j-1/2,k) based on UGRID reference 
650 !!
651 !!    AUTHOR
652 !!    ------
653 !!    F. Visentin   *CNRS/LA*                 
654 !!
655 !!    MODIFICATIONS
656 !!    -------------
657 !!
658 !-------------------------------------------------------------------------------
659 !
660 USE MODE_ll
661 USE MODD_LUNIT
662 USE MODD_CONF
663 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
664 USE MODI_GET_HALO
665 !
666 IMPLICIT NONE
667 !
668 !*       0.1   Declarations of dummy arguments :
669 !
670 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
671 !
672 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
673 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRVCT ! contrav. comp. on MASS GRID
674 !
675 ! output source term
676 !
677 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
678 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2 ! halo2 for the field at t
679 !
680 !
681 !*       0.2   Declarations of local variables :
682 !
683 INTEGER :: IIB,IJB    ! Begining useful area in x,y,z directions
684 INTEGER :: IIE,IJE    ! End useful area in x,y,z directions
685 INTEGER::  IS,IN      ! Coordinate of third order diffusion area
686 !
687 INTEGER:: ILUOUT,IRESP   ! for prints
688 !
689 ! intermediate reconstruction fluxes for positive wind case
690 !
691 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFPOS1, ZFPOS2
692 !
693 ! intermediate reconstruction fluxes for negative wind case
694 ! we need only one since ZFNEG2 = ZFPOS2
695 !
696 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFNEG1, ZFNEG2
697 !
698 ! smoothness indicators for positive wind case
699 !
700 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBPOS1, ZBPOS2
701 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBNEG1, ZBNEG2
702 !
703 ! WENO weights
704 !
705 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMP1, ZOMP2
706 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMN1, ZOMN2
707 !
708 ! standard weights
709 !
710 REAL, PARAMETER :: ZGAMMA1 = 1./3.
711 REAL, PARAMETER :: ZGAMMA2 = 2./3.
712 !
713 REAL, PARAMETER :: ZEPS = 1.0E-15
714 !
715 !-----------------------------------------------------------------------------
716 !
717 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
718 !                 ------------------------------
719 !
720 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
721 !
722 !---------------------------------------------------------------------------
723 !
724 !*       0.4.   INITIALIZE THE FIELD 
725 !               ---------------------
726 !
727 PR(:,:,:) = 0.0
728 !
729 ZFPOS1 = 0.0
730 ZFPOS2 = 0.0
731 ZFNEG1 = 0.0
732 ZFNEG2 = 0.0
733 ZBPOS1 = 0.0
734 ZBPOS2 = 0.0
735 ZBNEG1 = 0.0
736 ZBNEG2 = 0.0
737 ZOMP1  = 0.0
738 ZOMP2  = 0.0
739 ZOMN1  = 0.0
740 ZOMN2  = 0.0 
741 !
742 !-------------------------------------------------------------------------------
743 !
744 SELECT CASE ( HLBCY(1) ) ! 
745 !
746 !*       1.1    CYCLIC CASE IN THE Y DIRECTION:
747 !
748 CASE ('CYCL')          ! In that case one must have HLBCY(1) == HLBCY(2)
749 !
750 !!$  IF(NHALO == 1) THEN
751     IS=IJB
752     IN=IJE
753 !!$  ELSE
754 !!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
755 !!$    WRITE(ILUOUT,*) 'ERROR : 4th order advection in CYCLic case '
756 !!$    WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
757 !!$    CALL ABORT
758 !!$    STOP
759 !!$  END IF
760 !
761 ! intermediate fluxes for positive wind case
762 !
763   ZFPOS1(:,IS+1:IN+1,:) = 0.5 * (3.0*PSRC(:,IS:IN,:) - PSRC(:,IS-1:IN-1,:))
764   ZFPOS1(:,IS,       :) = 0.5 * (3.0*PSRC(:,IS-1, :) - TPHALO2%SOUTH(:,:))
765 !!  ZFPOS1(:,IS-1,     :) = - 999.
766 !
767   ZFPOS2(:,IS:IN+1,:) = 0.5 * (PSRC(:,IS-1:IN,:) + PSRC(:,IS:IN+1,:))
768   ZFPOS2(:,IS-1,   :) = 0.5 * (TPHALO2%SOUTH(:,:) + PSRC(:,IS-1,   :))
769 !
770   ZFNEG1(:,IS-1:IN,:) = 0.5 * (3.0*PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:))
771   ZFNEG1(:,IN+1,   :) = 0.5 * (3.0*PSRC(:,IN+1,   :) - TPHALO2%NORTH(:,:))
772 !
773   ZFNEG2(:,IS:IN+1,:) = 0.5 * (PSRC(:,IS:IN+1,:) + PSRC(:,IS-1:IN,:))
774   ZFNEG2(:,IS-1,   :) = 0.5 * (PSRC(:,IS-1,   :) + TPHALO2%SOUTH(:,:))
775 !
776 ! smoothness indicators for positive wind case
777 !
778   ZBPOS1(:,IS+1:IN+1,:) = (PSRC(:,IS:IN,:) - PSRC(:,IS-1:IN-1,:))**2
779   ZBPOS1(:,IS,       :) = (PSRC(:,IS-1,   :) - TPHALO2%SOUTH(:,:))**2
780 !!  ZBPOS1(:,IS-1,     :) = - 999. 
781 !
782   ZBPOS2(:,IS:IN+1,:) = (PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:))**2
783   ZBPOS2(:,IS-1,   :) = (PSRC(:,IS-1,   :) - TPHALO2%SOUTH(:,:))**2
784 !
785 ! smoothness indicators for negative wind case
786 !
787   ZBNEG1(:,IS-1:IN,:) = (PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:))**2
788   ZBNEG1(:,IN+1,   :) = (PSRC(:,IN+1,   :) - TPHALO2%NORTH(:,:))**2
789 !
790   ZBNEG2(:,IS:IN+1,:) = (PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:))**2
791   ZBNEG2(:,IS-1,   :) = (TPHALO2%SOUTH(:,:) - PSRC(:,IS-1,:))**2
792 !
793 ! WENO weights
794 !
795   ZOMP1 = ZGAMMA1 / (ZEPS + ZBPOS1)**2
796   ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2
797   ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2
798   ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2
799 !
800 ! WENO fluxes
801 !
802   PR = (ZOMP2/(ZOMP1+ZOMP2) * ZFPOS2 +                           &
803        (ZOMP1/(ZOMP1+ZOMP2) * ZFPOS1)) * (0.5+SIGN(0.5,PRVCT)) + &
804        (ZOMN2/(ZOMN1+ZOMN2) * ZFNEG2 +                           &
805        (ZOMN1/(ZOMN1+ZOMN2) * ZFNEG1)) * (0.5-SIGN(0.5,PRVCT))
806 !
807 !
808 !       OPEN, WALL, NEST CASE IN THE Y DIRECTION
809 !
810 CASE ('OPEN','WALL','NEST')
811 !
812   IS=IJB
813   IN=IJE
814 !
815 !       USE A FIRST ORDER UPSTREAM SCHEME AT THE PHYSICAL BORDER
816 !
817   IF(LSOUTH_ll()) THEN
818     PR(:,IS,:) = PSRC(:,IS-1,:) * (0.5+SIGN(0.5,PRVCT(:,IS,:))) + PSRC(:,IS,:) * (0.5-SIGN(0.5,PRVCT(:,IS,:)))
819 !
820 !!$  ELSEIF (NHALO == 1) THEN
821   ELSE
822     ZFPOS1(:,IS,:) = 0.5 * (3.0*PSRC(:,IS-1,:) - TPHALO2%SOUTH(:,:))
823     ZFPOS2(:,IS,:) = 0.5 * (PSRC(:,IS-1,:) + PSRC(:,IS,:))
824     ZBPOS1(:,IS,:) = (PSRC(:,IS-1,:) - TPHALO2%SOUTH(:,:))**2
825     ZBPOS2(:,IS,:) = (PSRC(:,IS,  :) - PSRC(:,IS-1,:))**2
826 !
827     ZFNEG1(:,IS,:) = 0.5 * (3.0*PSRC(:,IS,:) - PSRC(:,IS+1,:))
828     ZFNEG2(:,IS,:) = 0.5 * (PSRC(:,IS,    :) + PSRC(:,IS-1,:))
829     ZBNEG1(:,IS,:) = (PSRC(:,IS,  :) - PSRC(:,IS+1,:))**2
830     ZBNEG2(:,IS,:) = (PSRC(:,IS-1,:) - PSRC(:,IS,  :))**2
831 !
832     ZOMP1(:,IS,:) = ZGAMMA1 / (ZEPS + ZBPOS1(:,IS,:))**2
833     ZOMP2(:,IS,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS,:))**2
834     ZOMN1(:,IS,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS,:))**2
835     ZOMN2(:,IS,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS,:))**2
836 !
837     PR(:,IS,:) = (ZOMP2(:,IS,:)/(ZOMP1(:,IS,:)+ZOMP2(:,IS,:)) * ZFPOS2(:,IS,:) + &
838        (ZOMP1(:,IS,:)/(ZOMP1(:,IS,:)+ZOMP2(:,IS,:)) * ZFPOS1(:,IS,:))) * (0.5+SIGN(0.5,PRVCT(:,IS,:))) + &
839        (ZOMN2(:,IS,:)/(ZOMN1(:,IS,:)+ZOMN2(:,IS,:)) * ZFNEG2(:,IS,:) +  &
840        (ZOMN1(:,IS,:)/(ZOMN1(:,IS,:)+ZOMN2(:,IS,:)) * ZFNEG1(:,IS,:))) * (0.5-SIGN(0.5,PRVCT(:,IS,:)))
841 !
842   ENDIF
843 !
844   IF(LNORTH_ll()) THEN
845     PR(:,IN+1,:) = PSRC(:,IN,:) * (0.5+SIGN(0.5,PRVCT(:,IN+1,:))) + PSRC(:,IN+1,:) * (0.5-SIGN(0.5,PRVCT(:,IN+1,:)))
846 !
847 !!$  ELSEIF (NHALO == 1) THEN
848   ELSE
849     ZFPOS1(:,IN+1,:) = 0.5 * (3.0*PSRC(:,IN,:) - PSRC(:,IN-1,:))
850     ZFPOS2(:,IN+1,:) = 0.5 * (PSRC(:,IN,    :) + PSRC(:,IN+1,:))
851     ZBPOS1(:,IN+1,:) = (PSRC(:,IN,:) - PSRC(:,IN-1,:))**2
852     ZBPOS2(:,IN+1,:) = (PSRC(:,IN+1,:) - PSRC(:,IN,:))**2
853 !
854     ZFNEG1(:,IN+1,:) = 0.5 * (3.0*PSRC(:,IN+1,:) - TPHALO2%NORTH(:,:))
855     ZFNEG2(:,IN+1,:) = 0.5 * (PSRC(:,IN+1,    :) + PSRC(:,IN,:))
856     ZBNEG1(:,IN+1,:) = (PSRC(:,IN+1,:) - TPHALO2%NORTH(:,:))**2
857     ZBNEG2(:,IN+1,:) = (PSRC(:,IN,  :) - PSRC(:,IN+1,:))**2
858 !
859     ZOMP1(:,IN+1,:) = ZGAMMA1 / (ZEPS + ZBPOS1(:,IN+1,:))**2
860     ZOMP2(:,IN+1,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IN+1,:))**2
861     ZOMN1(:,IN+1,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IN+1,:))**2
862     ZOMN2(:,IN+1,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IN+1,:))**2
863 !
864     PR(:,IN+1,:) = (ZOMP2(:,IN+1,:)/(ZOMP1(:,IN+1,:)+ZOMP2(:,IN+1,:)) * ZFPOS2(:,IN+1,:) +   &
865        (ZOMP1(:,IN+1,:)/(ZOMP1(:,IN+1,:)+ZOMP2(:,IN+1,:)) * ZFPOS1(:,IN+1,:))) * (0.5+SIGN(0.5,PRVCT(:,IN+1,:))) + &
866        (ZOMN2(:,IN+1,:)/(ZOMN1(:,IN+1,:)+ZOMN2(:,IN+1,:)) * ZFNEG2(:,IN+1,:) +     &
867        (ZOMN1(:,IN+1,:)/(ZOMN1(:,IN+1,:)+ZOMN2(:,IN+1,:)) * ZFNEG1(:,IN+1,:))) * (0.5-SIGN(0.5,PRVCT(:,IN+1,:)))
868 !
869   ENDIF
870 !
871 !      USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE 
872 !
873   ZFPOS1(:,IS+1:IN,:) = 0.5 * (3.0*PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:))
874   ZFPOS2(:,IS+1:IN,:) = 0.5 * (PSRC(:,IS:IN-1,    :) + PSRC(:,IS+1:IN,  :))
875   ZBPOS1(:,IS+1:IN,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:))**2
876   ZBPOS2(:,IS+1:IN,:) = (PSRC(:,IS+1:IN,:) - PSRC(:,IS:IN-1,  :))**2
877 !
878   ZFNEG1(:,IS+1:IN,:) = 0.5 * (3.0*PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:))
879   ZFNEG2(:,IS+1:IN,:) = 0.5 * (PSRC(:,IS+1:IN,    :) + PSRC(:,IS:IN-1,  :))
880   ZBNEG1(:,IS+1:IN,:) = (PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:))**2
881   ZBNEG2(:,IS+1:IN,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS+1:IN,:))**2
882 !
883   ZOMP1(:,IS+1:IN,:) = ZGAMMA1 / (ZEPS + ZBPOS1(:,IS+1:IN,:))**2
884   ZOMP2(:,IS+1:IN,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS+1:IN,:))**2
885   ZOMN1(:,IS+1:IN,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS+1:IN,:))**2
886   ZOMN2(:,IS+1:IN,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS+1:IN,:))**2
887 !
888   PR(:,IS+1:IN,:) = (ZOMP2(:,IS+1:IN,:)/(ZOMP1(:,IS+1:IN,:)+ZOMP2(:,IS+1:IN,:)) * ZFPOS2(:,IS+1:IN,:) +  &
889        (ZOMP1(:,IS+1:IN,:)/(ZOMP1(:,IS+1:IN,:)+ZOMP2(:,IS+1:IN,:)) * ZFPOS1(:,IS+1:IN,:))) * (0.5+SIGN(0.5,PRVCT(:,IS+1:IN,:))) + &
890        (ZOMN2(:,IS+1:IN,:)/(ZOMN1(:,IS+1:IN,:)+ZOMN2(:,IS+1:IN,:)) * ZFNEG2(:,IS+1:IN,:) + &
891        (ZOMN1(:,IS+1:IN,:)/(ZOMN1(:,IS+1:IN,:)+ZOMN2(:,IS+1:IN,:)) * ZFNEG1(:,IS+1:IN,:))) * (0.5-SIGN(0.5,PRVCT(:,IS+1:IN,:)))
892 !
893 END SELECT
894 !
895 PR = PR * PRVCT
896 CALL GET_HALO(PR)
897 !
898 END SUBROUTINE ADVEC_WENO_K_2_MY
899 !-------------------------------------------------------------------------------
900 !
901 !     #############################################################
902       SUBROUTINE ADVEC_WENO_K_2_VY(HLBCY, PSRC, PRVCT, PR, TPHALO2)
903 !     #############################################################
904 !!
905 !!**** Computes PRVCT * PVT. Upstream fluxes of V in Y direction.  
906 !!     Input PVT is on V Grid 'ie' (i,j,k) based on VGRID reference
907 !!     Output PR is on mass Grid 'ie' (i,j+1/2,k) based on VGRID reference
908 !!
909 !!    AUTHOR
910 !!    ------
911 !!    F. Visentin   *CNRS/LA*
912 !!
913 !!    MODIFICATIONS
914 !!    -------------
915 !!
916 !-------------------------------------------------------------------------------
917 !
918 USE MODE_ll
919 USE MODD_LUNIT
920 USE MODD_CONF
921 USE MODD_ARGSLIST_ll, ONLY : HALO2LIST_ll
922 USE MODI_GET_HALO
923 !
924 IMPLICIT NONE
925 !
926 !*       0.1   Declarations of dummy arguments :
927 !
928 CHARACTER (LEN=4), DIMENSION(2), INTENT(IN) :: HLBCY  ! Y direction LBC type
929 !
930 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on U grid at t
931 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRVCT ! contrav. comp. on MASS GRID
932 !
933 ! output source term
934 REAL, DIMENSION(:,:,:), INTENT(OUT) :: PR
935 TYPE(HALO2_ll), OPTIONAL, POINTER :: TPHALO2      ! halo2 for the field at t
936 !
937 !*       0.2   Declarations of local variables :
938 !
939 INTEGER :: IIB,IJB    ! Begining useful area in x,y,z directions
940 INTEGER :: IIE,IJE    ! End useful area in x,y,z directions
941 INTEGER::  IS,IN      ! Coordinate of third order diffusion area
942 !
943 INTEGER:: ILUOUT,IRESP   ! for prints
944 !
945 ! intermediate reconstruction fluxes for positive wind case
946 !
947 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFPOS1, ZFPOS2
948 !
949 ! intermediate reconstruction fluxes for negative wind case
950 ! we need only one since ZFNEG2 = ZFPOS2
951 !
952 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFNEG1, ZFNEG2
953 !
954 ! smoothness indicators for positive wind case
955 !
956 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBPOS1, ZBPOS2
957 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBNEG1, ZBNEG2
958 !
959 ! WENO weights
960 !
961 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMP1, ZOMP2
962 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMN1, ZOMN2
963 !
964 ! standard weights
965 !
966 REAL, PARAMETER :: ZGAMMA1 = 1./3.
967 REAL, PARAMETER :: ZGAMMA2 = 2./3.
968 !
969 REAL, PARAMETER :: ZEPS = 1.0E-15
970 !
971 !----------------------------------------------------------------------------
972 !
973 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
974 !                 ------------------------------
975 !
976 CALL GET_INDICE_ll(IIB,IJB,IIE,IJE)
977 !
978 !--------------------------------------------------------------------------
979 !
980 !*       0.4.   INITIALIZE THE FIELD 
981 !               ---------------------
982 !
983 PR(:,:,:) = 0.0
984 !
985 ZFPOS1 = 0.0
986 ZFPOS2 = 0.0
987 ZFNEG1 = 0.0
988 ZFNEG2 = 0.0
989 ZBPOS1 = 0.0
990 ZBPOS2 = 0.0
991 ZBNEG1 = 0.0
992 ZBNEG2 = 0.0
993 ZOMP1  = 0.0
994 ZOMP2  = 0.0
995 ZOMN1  = 0.0
996 ZOMN2  = 0.0
997 !
998 !-------------------------------------------------------------------------------
999 !
1000 SELECT CASE ( HLBCY(1) ) ! Y direction LBC type: (1) for left side
1001 !
1002 !*       1.1    CYCLIC CASE IN THE Y DIRECTION:
1003 !
1004 CASE ('CYCL')          ! In that case one must have HLBCX(1) == HLBCX(2)
1005 !
1006 !!$  IF(NHALO == 1) THEN
1007     IS=IJB
1008     IN=IJE
1009 !!$  ELSE
1010 !!$    CALL FMLOOK_ll(CLUOUT0,CLUOUT0,ILUOUT,IRESP)
1011 !!$    WRITE(ILUOUT,*) 'ERROR : 4th order advection in CYCLic case '
1012 !!$    WRITE(ILUOUT,*) 'cannot be used with NHALO=2'
1013 !!$    CALL ABORT
1014 !!$    STOP
1015 !!$  END IF
1016 !
1017 ! intermediate fluxes for positive wind case
1018 !
1019   ZFPOS1(:,IS:IN+1,:) = 0.5 * (3.0*PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:))
1020   ZFPOS1(:,IS-1,   :) = 0.5 * (3.0*PSRC(:,IS-1,   :) - TPHALO2%SOUTH(:,:))
1021 !
1022   ZFPOS2(:,IS-1:IN,:) = 0.5 * (PSRC(:,IS-1:IN,:) + PSRC(:,IS:IN+1,:))
1023   ZFPOS2(:,IN+1,   :) = 0.5 * (PSRC(:,IN+1,   :) + TPHALO2%NORTH(:,:))
1024 !
1025 ! intermediate flux for negative wind case
1026 !
1027   ZFNEG1(:,IS-1:IN-1,:) = 0.5 * (3.0*PSRC(:,IS:IN,:) - PSRC(:,IS+1:IN+1,:))
1028   ZFNEG1(:,IN,   :) = 0.5 * (3.0*PSRC(:,IN+1,   :) - TPHALO2%NORTH(:,:))
1029 !
1030   ZFNEG2(:,IS-1:IN,:) = 0.5 * (PSRC(:,IS-1:IN,:) + PSRC(:,IS:IN+1,:))
1031   ZFNEG2(:,IN+1,   :) = 0.5 * (PSRC(:,IN+1,   :) + TPHALO2%NORTH(:,:))
1032 !
1033 ! smoothness indicators for positive wind case
1034 !
1035   ZBPOS1(:,IS:IN+1,:) = (PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:))**2
1036   ZBPOS1(:,IS-1,   :) = (PSRC(:,IS-1,   :) - TPHALO2%SOUTH(:,:))**2
1037 !
1038   ZBPOS2(:,IS-1:IN,:) = (PSRC(:,IS:IN+1,:) - PSRC(:,IS-1:IN,:))**2
1039   ZBPOS2(:,IN+1,   :) = (TPHALO2%NORTH(:,:) - PSRC(:,IN+1,     :))**2
1040 !
1041 ! smoothness indicators for negative wind case
1042 !
1043   ZBNEG1(:,IS-1:IN-1,:) = (PSRC(:,IS:IN,:) - PSRC(:,IS+1:IN+1,:))**2
1044   ZBNEG1(:,IN,       :) = (PSRC(:,IN+1, :) - TPHALO2%NORTH(:,:))**2
1045 !
1046   ZBNEG2(:,IS-1:IN,:) = (PSRC(:,IS-1:IN,:) - PSRC(:,IS:IN+1,:))**2
1047   ZBNEG2(:,IN+1,   :) = (PSRC(:,IN+1,   :) - TPHALO2%NORTH(:,:))**2 
1048 !
1049 ! WENO weights
1050 !
1051   ZOMP1 = ZGAMMA1 / (ZEPS + ZBPOS1)**2
1052   ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2
1053   ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2
1054   ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2
1055 !
1056   PR = (ZOMP2/(ZOMP1+ZOMP2) * ZFPOS2 +                           &
1057        (ZOMP1/(ZOMP1+ZOMP2) * ZFPOS1)) * (0.5+SIGN(0.5,PRVCT)) + &
1058        (ZOMN2/(ZOMN1+ZOMN2) * ZFNEG2 +                           &
1059        (ZOMN1/(ZOMN1+ZOMN2) * ZFNEG1)) * (0.5-SIGN(0.5,PRVCT))
1060 !
1061 !
1062 !       OPEN, WALL, NEST CASE IN THE Y DIRECTION
1063 !
1064 CASE ('OPEN','WALL','NEST')
1065 !
1066   IS=IJB
1067   IN=IJE
1068 !
1069 !       USE A FIRST ORDER UPSTREAM SCHEME AT THE PHYSICAL BORDER
1070 !
1071   IF(LSOUTH_ll()) THEN
1072     PR(:,IS-1,:) = PSRC(:,IS-1,:) * (0.5+SIGN(0.5,PRVCT(:,IS-1,:))) + PSRC(:,IS,:) * (0.5-SIGN(0.5,PRVCT(:,IS-1,:)))
1073 !
1074 !!$  ELSEIF (NHALO == 1) THEN
1075   ELSE
1076     ZFPOS1(:,IS-1,:) = 0.5 * (3.0*PSRC(:,IS-1,:) - TPHALO2%SOUTH(:,:))
1077     ZFPOS2(:,IS-1,:) = 0.5 * (PSRC(:,IS-1,    :) + PSRC(:,IS,:))
1078     ZBPOS1(:,IS-1,:) = (PSRC(:,IS-1,:) - TPHALO2%SOUTH(:,:))**2
1079     ZBPOS2(:,IS-1,:) = (PSRC(:,IS,  :) - PSRC(:,IS-1,:))**2
1080 !
1081     ZFNEG1(:,IS-1,:) = 0.5 * (3.0*PSRC(:,IS,:) - PSRC(:,IS+1,:))
1082     ZFNEG2(:,IS-1,:) = 0.5 * (PSRC(:,IS-1,  :) + PSRC(:,IS,:))
1083     ZBNEG1(:,IS-1,:) = (PSRC(:,IS,:) - PSRC(:,IS+1,:))**2
1084     ZBNEG2(:,IS-1,:) = (PSRC(:,IS-1,:) - PSRC(:,IS,:))**2
1085 !
1086     ZOMP1(:,IS-1,:) = ZGAMMA1 / (ZEPS + ZBPOS1(:,IS-1,:))**2
1087     ZOMP2(:,IS-1,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS-1,:))**2
1088     ZOMN1(:,IS-1,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS-1,:))**2
1089     ZOMN2(:,IS-1,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS-1,:))**2
1090 !
1091     PR(:,IS-1,:) = (ZOMP2(:,IS-1,:)/(ZOMP1(:,IS-1,:)+ZOMP2(:,IS-1,:)) * ZFPOS2(:,IS-1,:) +   &
1092        (ZOMP1(:,IS-1,:)/(ZOMP1(:,IS-1,:)+ZOMP2(:,IS-1,:)) * ZFPOS1(:,IS-1,:))) * (0.5+SIGN(0.5,PRVCT(:,IS-1,:))) + &
1093        (ZOMN2(:,IS-1,:)/(ZOMN1(:,IS-1,:)+ZOMN2(:,IS-1,:)) * ZFNEG2(:,IS-1,:) +        &
1094        (ZOMN1(:,IS-1,:)/(ZOMN1(:,IS-1,:)+ZOMN2(:,IS-1,:)) * ZFNEG1(:,IS-1,:))) * (0.5-SIGN(0.5,PRVCT(:,IS-1,:)))
1095 !
1096   ENDIF
1097 !
1098   IF(LNORTH_ll()) THEN
1099     PR(:,IN,:) = PSRC(:,IN,:) * (0.5+SIGN(0.5,PRVCT(:,IN,:))) + PSRC(:,IN+1,:) * (0.5-SIGN(0.5,PRVCT(:,IN,:)))
1100 !
1101 !!$  ELSEIF (NHALO == 1) THEN
1102   ELSE
1103     ZFPOS1(:,IN,:) = 0.5 * (3.0*PSRC(:,IN,:) - PSRC(:,IN-1,:))
1104     ZFPOS2(:,IN,:) = 0.5 * (PSRC(:,IN,    :) + PSRC(:,IN+1,:))
1105     ZBPOS1(:,IN,:) = (PSRC(:,IN,  :) - PSRC(:,IN-1,:))**2
1106     ZBPOS2(:,IN,:) = (PSRC(:,IN+1,:) - PSRC(:,IN,  :))**2
1107 !
1108     ZFNEG1(:,IN,:) = 0.5 * (3.0*PSRC(:,IN+1,:) - TPHALO2%NORTH(:,:))
1109     ZFNEG2(:,IN,:) = 0.5 * (PSRC(:,IN,      :) + PSRC(:,IN+1,:))
1110     ZBNEG1(:,IN,:) = (PSRC(:,IN+1,:) - TPHALO2%NORTH(:,:))**2
1111     ZBNEG2(:,IN,:) = (PSRC(:,IN,  :) - PSRC(:,IN+1,:))**2
1112 !
1113     ZOMP1(:,IN,:) = ZGAMMA1 / (ZEPS + ZBPOS1(:,IN,:))**2
1114     ZOMP2(:,IN,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IN,:))**2
1115     ZOMN1(:,IN,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IN,:))**2
1116     ZOMN2(:,IN,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IN,:))**2
1117 !
1118     PR(:,IN,:) = (ZOMP2(:,IN,:)/(ZOMP1(:,IN,:)+ZOMP2(:,IN,:)) * ZFPOS2(:,IN,:) + &
1119        (ZOMP1(:,IN,:)/(ZOMP1(:,IN,:)+ZOMP2(:,IN,:)) * ZFPOS1(:,IN,:))) * (0.5+SIGN(0.5,PRVCT(:,IN,:))) + &
1120        (ZOMN2(:,IN,:)/(ZOMN1(:,IN,:)+ZOMN2(:,IN,:)) * ZFNEG2(:,IN,:) +  &
1121        (ZOMN1(:,IN,:)/(ZOMN1(:,IN,:)+ZOMN2(:,IN,:)) * ZFNEG1(:,IN,:))) * (0.5-SIGN(0.5,PRVCT(:,IN,:)))
1122 !
1123   ENDIF
1124 !
1125 !      USE A THIRD ORDER UPSTREAM WENO SCHEME ELSEWHERE 
1126 !
1127   ZFPOS1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:))
1128   ZFPOS2(:,IS:IN-1,:) = 0.5 * (PSRC(:,IS:IN-1,    :) + PSRC(:,IS+1:IN,  :))
1129   ZBPOS1(:,IS:IN-1,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS-1:IN-2,:))**2
1130   ZBPOS2(:,IS:IN-1,:) = (PSRC(:,IS+1:IN,:) - PSRC(:,IS:IN-1,  :))**2
1131 !  
1132   ZFNEG1(:,IS:IN-1,:) = 0.5 * (3.0*PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:))  
1133   ZFNEG2(:,IS:IN-1,:) = 0.5 * (PSRC(:,IS:IN-1,    :) + PSRC(:,IS+1:IN,  :))
1134   ZBNEG1(:,IS:IN-1,:) = (PSRC(:,IS+1:IN,:) - PSRC(:,IS+2:IN+1,:))**2
1135   ZBNEG2(:,IS:IN-1,:) = (PSRC(:,IS:IN-1,:) - PSRC(:,IS+1:IN,  :))**2
1136 !
1137   ZOMP1(:,IS:IN-1,:) = ZGAMMA1 / (ZEPS + ZBPOS1(:,IS:IN-1,:))**2
1138   ZOMP2(:,IS:IN-1,:) = ZGAMMA2 / (ZEPS + ZBPOS2(:,IS:IN-1,:))**2
1139   ZOMN1(:,IS:IN-1,:) = ZGAMMA1 / (ZEPS + ZBNEG1(:,IS:IN-1,:))**2
1140   ZOMN2(:,IS:IN-1,:) = ZGAMMA2 / (ZEPS + ZBNEG2(:,IS:IN-1,:))**2
1141 !
1142   PR(:,IS:IN-1,:) = (ZOMP2(:,IS:IN-1,:)/(ZOMP1(:,IS:IN-1,:)+ZOMP2(:,IS:IN-1,:)) * ZFPOS2(:,IS:IN-1,:) + &
1143        (ZOMP1(:,IS:IN-1,:)/(ZOMP1(:,IS:IN-1,:)+ZOMP2(:,IS:IN-1,:)) * ZFPOS1(:,IS:IN-1,:))) * (0.5+SIGN(0.5,PRVCT(:,IS:IN-1,:))) + &
1144        (ZOMN2(:,IS:IN-1,:)/(ZOMN1(:,IS:IN-1,:)+ZOMN2(:,IS:IN-1,:)) * ZFNEG2(:,IS:IN-1,:) + &
1145        (ZOMN1(:,IS:IN-1,:)/(ZOMN1(:,IS:IN-1,:)+ZOMN2(:,IS:IN-1,:)) * ZFNEG1(:,IS:IN-1,:))) * (0.5-SIGN(0.5,PRVCT(:,IS:IN-1,:)))
1146 !
1147 END SELECT
1148 !
1149 PR = PR * PRVCT
1150 CALL GET_HALO(PR)
1151 !
1152 END SUBROUTINE ADVEC_WENO_K_2_VY
1153 !
1154 !-------------------------------------------------------------------------------
1155 !
1156 !     ############################################
1157       FUNCTION WENO_K_2_WZ(PSRC, PRWCT) RESULT(PR)
1158 !     ############################################
1159 !!
1160 !!* Computes PRWCT * PWT. Upstream fluxes of W in Z direction.  
1161 !!  Input PWT is on W Grid 'ie' (i,j,k) based on WGRID reference
1162 !!  Output PR is on mass Grid 'ie' (i,j,k+1/2) based on WGRID reference
1163 !!
1164 !!    AUTHOR
1165 !!    ------
1166 !!    F. Visentin   *CNRS/LA*
1167 !!
1168 !!    MODIFICATIONS
1169 !!    -------------
1170 !!
1171 !-------------------------------------------------------------------------------
1172 !
1173 USE MODE_ll
1174 USE MODD_CONF
1175 USE MODD_PARAMETERS,ONLY: JPVEXT
1176 USE MODI_GET_HALO
1177 !
1178 IMPLICIT NONE
1179 !
1180 !*       0.1   Declarations of dummy arguments :
1181 !
1182 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on W grid at t
1183 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRWCT ! contrav. comp. on MASS GRID
1184 !
1185 ! output source term
1186 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: PR
1187 !
1188 !*       0.2   Declarations of local variables :
1189 !
1190 INTEGER :: IB    ! Begining useful area in x,y,z directions
1191 INTEGER :: IT    ! End useful area in x,y,z directions
1192 !
1193 ! WENO-related variables:
1194 ! intermediate reconstruction fluxes for positive wind case
1195 !
1196 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFPOS1, ZFPOS2
1197 !
1198 ! intermediate reconstruction fluxes for negative wind case
1199 ! we need only one since ZFNEG2 = ZFPOS2
1200 !
1201 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFNEG1, ZFNEG2
1202 !
1203 ! smoothness indicators for positive wind case
1204 !
1205 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBPOS1, ZBPOS2
1206 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBNEG1, ZBNEG2
1207 !
1208 ! WENO weights
1209 !
1210 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMP1, ZOMP2
1211 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMN1, ZOMN2
1212 !
1213 ! standard weights
1214 !
1215 REAL, PARAMETER :: ZGAMMA1 = 1./3.
1216 REAL, PARAMETER :: ZGAMMA2 = 2./3.
1217 !
1218 REAL, PARAMETER :: ZEPS = 1.0E-15
1219 !
1220 !-------------------------------------------------------------------------------
1221 !
1222 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
1223 !                 ------------------------------
1224 !
1225 IB = 1 + JPVEXT
1226 IT = SIZE(PSRC,3) - JPVEXT
1227 !
1228 PR(:,:,:) = 0.0
1229 !
1230 ZFPOS1 = 0.0
1231 ZFPOS2 = 0.0
1232 ZFNEG1 = 0.0
1233 ZFNEG2 = 0.0
1234 ZBPOS1 = 0.0
1235 ZBPOS2 = 0.0
1236 ZBNEG1 = 0.0
1237 ZBNEG2 = 0.0
1238 ZOMP1  = 0.0
1239 ZOMP2  = 0.0
1240 ZOMN1  = 0.0
1241 ZOMN2  = 0.0 
1242 !
1243 ! intermediate fluxes at the mass point on Wgrid w(i,j,k+1/2) for positive 
1244 ! wind case (L. to the R.)
1245 ! (r=1 for the first stencil ZFPOS1, r=0 for the second ZFPOS2)
1246 !
1247 ZFPOS1(:,:,IB:IT-1) = 0.5 * (3.0*PSRC(:,:,IB:IT-1) - PSRC(:,:,IB-1:IT-2))
1248 ZFPOS2(:,:,IB:IT-1) = 0.5 * (PSRC(:,:,IB:IT-1) + PSRC(:,:,IB+1:IT))
1249 !
1250 ! intermediate flux at the mass point on Wgrid w(i,j,k+1/2) for negative 
1251 ! wind case (R. to the L.)
1252 ! (r=-1 for the first stencil ZFNEG1, r=0 for the second ZFNEG2=ZFPOS2)
1253 !
1254 ZFNEG1(:,:,IB-1:IT-1) = 0.5 * (3.0*PSRC(:,:,IB:IT) - PSRC(:,:,IB+1:IT+1))
1255 ZFNEG2(:,:,IB-1:IT) = 0.5 * (PSRC(:,:,IB-1:IT) + PSRC(:,:,IB:IT+1))
1256 !
1257 ! smoothness indicators for positive wind case
1258 !
1259 ZBPOS1(:,:,IB:IT-1) = (PSRC(:,:,IB:IT-1) - PSRC(:,:,IB-1:IT-2))**2
1260 ZBPOS2(:,:,IB:IT-1) = (PSRC(:,:,IB+1:IT) - PSRC(:,:,IB:IT-1))**2
1261 !
1262 ! smoothness indicators for negative wind case
1263 !
1264 ZBNEG1(:,:,IB-1:IT-1) = (PSRC(:,:,IB:IT) - PSRC(:,:,IB+1:IT+1))**2
1265 ZBNEG2(:,:,IB-1:IT) = (PSRC(:,:,IB-1:IT) - PSRC(:,:,IB:IT+1))**2
1266 !
1267 ! WENO weights
1268 !
1269 ZOMP1 = ZGAMMA1 / (ZEPS + ZBPOS1)**2
1270 ZOMP2 = ZGAMMA2 / (ZEPS + ZBPOS2)**2
1271 ZOMN1 = ZGAMMA1 / (ZEPS + ZBNEG1)**2
1272 ZOMN2 = ZGAMMA2 / (ZEPS + ZBNEG2)**2
1273 !
1274 ! WENO fluxes
1275 !
1276 PR(:,:,IB:IT-1) = (ZOMP2(:,:,IB:IT-1)/(ZOMP1(:,:,IB:IT-1)+ZOMP2(:,:,IB:IT-1))* &
1277                                                          ZFPOS2(:,:,IB:IT-1) + &
1278                   (ZOMP1(:,:,IB:IT-1)/(ZOMP1(:,:,IB:IT-1)+ZOMP2(:,:,IB:IT-1))* &
1279                  ZFPOS1(:,:,IB:IT-1))) * (0.5+SIGN(0.5,PRWCT(:,:,IB:IT-1) )) + &
1280                   (ZOMN2(:,:,IB:IT-1)/(ZOMN1(:,:,IB:IT-1)+ZOMN2(:,:,IB:IT-1))* &
1281                                                          ZFNEG2(:,:,IB:IT-1) + &
1282                   (ZOMN1(:,:,IB:IT-1)/(ZOMN1(:,:,IB:IT-1)+ZOMN2(:,:,IB:IT-1))* &
1283                  ZFNEG1(:,:,IB:IT-1))) * (0.5-SIGN(0.5,PRWCT(:,:,IB:IT-1) ))
1284 !
1285 PR(:,:,IB-1) = PSRC(:,:,IB-1) * (0.5+SIGN(0.5,PRWCT(:,:,IB-1) )) + &
1286                PSRC(:,:,IB)   * (0.5-SIGN(0.5,PRWCT(:,:,IB-1) ))
1287 PR(:,:,IT)   = PSRC(:,:,IT)   * (0.5+SIGN(0.5,PRWCT(:,:,IT) ))   + &
1288                PSRC(:,:,IT+1) * (0.5-SIGN(0.5,PRWCT(:,:,IT) ))
1289 PR(:,:,IT+1) = -999.
1290 !
1291 PR = PR * PRWCT
1292 CALL GET_HALO(PR)
1293 !
1294 END FUNCTION WENO_K_2_WZ
1295 !
1296 !-----------------------------------------------------------------------------
1297 !
1298 !     ############################################
1299       FUNCTION WENO_K_2_MZ(PSRC, PRWCT) RESULT(PR)
1300 !     ############################################
1301 !!
1302 !!* Computes PRWCT * PUT (or PRWCT * PVT). Upstream fluxes of U (or V) 
1303 !!  variables in Z direction.  
1304 !!  Input PUT is on U Grid 'ie' (i,j,k) based on UGRID reference                
1305 !!  Output PR is on mass Grid 'ie' (i,j,k-1/2) based on UGRID reference
1306 !!
1307 !!    AUTHOR
1308 !!    ------
1309 !!    F. Visentin   *CNRS/LA*
1310 !!
1311 !!    MODIFICATIONS
1312 !!    -------------
1313 !!
1314 !-------------------------------------------------------------------------------
1315 !
1316 USE MODE_ll
1317 USE MODD_CONF
1318 USE MODD_PARAMETERS,ONLY: JPVEXT
1319 USE MODI_GET_HALO
1320 !
1321 IMPLICIT NONE
1322 !
1323 !*       0.1   Declarations of dummy arguments :
1324 !
1325 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PSRC  ! variable on MASS grid at t
1326 REAL, DIMENSION(:,:,:), INTENT(IN)  :: PRWCT ! contrav. comp. on W grid
1327 !
1328 ! output source term
1329 !
1330 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: PR
1331 !
1332 !*       0.2   Declarations of local variables :
1333 !
1334 INTEGER :: IB    ! Begining useful area in x,y,z directions
1335 INTEGER :: IT    ! End useful area in x,y,z directions
1336 !
1337 ! WENO-related variables:
1338 !
1339 ! intermediate reconstruction fluxes for positive wind case
1340 !
1341 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFPOS1, ZFPOS2
1342 !
1343 ! intermediate reconstruction fluxes for negative wind case
1344 ! we need only one since ZFNEG2 = ZFPOS2
1345 !
1346 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZFNEG1, ZFNEG2
1347 !
1348 ! smoothness indicators for positive wind case
1349 !
1350 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBPOS1, ZBPOS2
1351 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZBNEG1, ZBNEG2
1352 !
1353 ! WENO weights
1354 !
1355 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMP1, ZOMP2
1356 REAL, DIMENSION(SIZE(PSRC,1),SIZE(PSRC,2),SIZE(PSRC,3)) :: ZOMN1, ZOMN2
1357 !
1358 ! standard weights
1359 !
1360 REAL, PARAMETER :: ZGAMMA1 = 1./3.
1361 REAL, PARAMETER :: ZGAMMA2 = 2./3.
1362 !
1363 REAL, PARAMETER :: ZEPS = 1.0E-15
1364 !
1365 !-------------------------------------------------------------------------------
1366 !
1367 !*       0.3.     COMPUTES THE DOMAIN DIMENSIONS
1368 !                 ------------------------------
1369 !
1370 IB = 1 + JPVEXT
1371 IT = SIZE(PSRC,3) - JPVEXT
1372 !
1373 PR(:,:,:) = 0.0
1374 !
1375 ZFPOS1 = 0.0
1376 ZFPOS2 = 0.0
1377 ZFNEG1 = 0.0
1378 ZFNEG2 = 0.0 
1379 ZBPOS1 = 0.0
1380 ZBPOS2 = 0.0
1381 ZBNEG1 = 0.0
1382 ZBNEG2 = 0.0
1383 ZOMP1  = 0.0
1384 ZOMP2  = 0.0
1385 ZOMN1  = 0.0
1386 ZOMN2  = 0.0
1387 !
1388 ! intermediate fluxes at the flux point on the Wgrid u(i,j,k-1/2) for 
1389 ! positive wind case
1390 !
1391 ZFPOS1(:,:,IB+1:IT) = 0.5 * (3.0*PSRC(:,:,IB:IT-1) - PSRC(:,:,IB-1:IT-2))
1392 ZFPOS2(:,:,IB+1:IT) = 0.5 * (PSRC(:,:,IB:IT-1) + PSRC(:,:,IB+1:IT))
1393 !
1394 ! intermediate flux at the flux point on the Wgrid u(i,j,k-1/2) for 
1395 ! negative wind case
1396 !
1397 ZFNEG1(:,:,IB+1:IT) = 0.5 * (3.0*PSRC(:,:,IB+1:IT) - PSRC(:,:,IB+2:IT+1))
1398 ZFNEG2(:,:,IB+1:IT) = 0.5 * (PSRC(:,:,IB:IT-1) + PSRC(:,:,IB+1:IT))
1399 !
1400 ! smoothness indicators for positive wind case
1401 !
1402 ZBPOS1(:,:,IB+1:IT) = (PSRC(:,:,IB:IT-1) - PSRC(:,:,IB-1:IT-2))**2
1403 ZBPOS2(:,:,IB+1:IT) = (PSRC(:,:,IB+1:IT) - PSRC(:,:,IB:IT-1))**2
1404 !
1405 ! smoothness indicators for negative wind case
1406 !
1407 ZBNEG1(:,:,IB+1:IT) = (PSRC(:,:,IB+1:IT) - PSRC(:,:,IB+2:IT+1))**2
1408 ZBNEG2(:,:,IB+1:IT) = (PSRC(:,:,IB:IT-1) - PSRC(:,:,IB+1:IT))**2
1409 !
1410 ! WENO weights
1411 !
1412 ZOMP1(:,:,IB+1:IT) = ZGAMMA1 / (ZEPS + ZBPOS1(:,:,IB+1:IT))**2
1413 ZOMP2(:,:,IB+1:IT) = ZGAMMA2 / (ZEPS + ZBPOS2(:,:,IB+1:IT))**2
1414 ZOMN1(:,:,IB+1:IT) = ZGAMMA1 / (ZEPS + ZBNEG1(:,:,IB+1:IT))**2
1415 ZOMN2(:,:,IB+1:IT) = ZGAMMA2 / (ZEPS + ZBNEG2(:,:,IB+1:IT))**2
1416 !
1417 PR(:,:,IB+1:IT) = (ZOMP2(:,:,IB+1:IT)/(ZOMP1(:,:,IB+1:IT)+ZOMP2(:,:,IB+1:IT))* &
1418                                                          ZFPOS2(:,:,IB+1:IT) + &
1419                   (ZOMP1(:,:,IB+1:IT)/(ZOMP1(:,:,IB+1:IT)+ZOMP2(:,:,IB+1:IT))* &
1420                  ZFPOS1(:,:,IB+1:IT))) * (0.5+SIGN(0.5,PRWCT(:,:,IB+1:IT) )) + &
1421                   (ZOMN2(:,:,IB+1:IT)/(ZOMN1(:,:,IB+1:IT)+ZOMN2(:,:,IB+1:IT))* &
1422                                                          ZFNEG2(:,:,IB+1:IT) + &
1423                   (ZOMN1(:,:,IB+1:IT)/(ZOMN1(:,:,IB+1:IT)+ZOMN2(:,:,IB+1:IT))* &
1424                  ZFNEG1(:,:,IB+1:IT))) * (0.5-SIGN(0.5,PRWCT(:,:,IB+1:IT) ))
1425 !
1426 PR(:,:,IB)   = PSRC(:,:,IB-1) * (0.5+SIGN(0.5,PRWCT(:,:,IB) ))   + &
1427                PSRC(:,:,IB)   * (0.5-SIGN(0.5,PRWCT(:,:,IB) ))
1428 PR(:,:,IT+1) = PSRC(:,:,IT)   * (0.5+SIGN(0.5,PRWCT(:,:,IT+1) )) + &
1429                PSRC(:,:,IT+1) * (0.5-SIGN(0.5,PRWCT(:,:,IT+1) ))
1430 !
1431 PR = PR * PRWCT
1432 CALL GET_HALO(PR)
1433 !
1434 END FUNCTION WENO_K_2_MZ