| blob | 35267f246bb2907629c7d2f44055f5c23b8c48d1 |
1 !----------------------------------------------------------------------
2 !#define BIT_FOR_BIT
3 !----------------------------------------------------------------------
4 #include "nmm_loop_basemacros.h"
5 #include "nmm_loop_macros.h"
6 !----------------------------------------------------------------------
7 !
8 !NCEP_MESO:MODEL_LAYER: HORIZONTAL AND VERTICAL ADVECTION
9 !
10 !----------------------------------------------------------------------
11 !
12 MODULE MODULE_ADVECTION
13 !
14 !----------------------------------------------------------------------
15 USE MODULE_MODEL_CONSTANTS
16 USE MODULE_EXT_INTERNAL
17 !----------------------------------------------------------------------
18 #if defined(DM_PARALLEL) && !defined(STUBMPI)
19 INCLUDE "mpif.h"
20 #endif
21 !----------------------------------------------------------------------
22 !
23 REAL,PARAMETER :: FF2=-0.64813,FF3=0.24520,FF4=-0.12189
24 REAL,PARAMETER :: FFC=1.533,FBC=1.-FFC
25 REAL :: CONSERVE_MIN=0.9,CONSERVE_MAX=1.1
26 !
27 !----------------------------------------------------------------------
28 !*** CRANK-NICHOLSON OFF-CENTER WEIGHTS FOR CURRENT AND FUTURE
29 !*** TIME LEVELS.
30 !-----------------------------------------------------------------------
31 !
32 REAL,PARAMETER :: WGT1=0.90
33 REAL,PARAMETER :: WGT2=2.-WGT1
34 !
35 !*** FOR CRANK_NICHOLSON CHECK ONLY.
36 !
37 INTEGER :: ITEST=47,JTEST=70
38 REAL :: ADTP,ADUP,ADVP,TTLO,TTUP,TULO,TUUP,TVLO,TVUP
39 !
40 !----------------------------------------------------------------------
41 CONTAINS
42 !
43 !***********************************************************************
44 SUBROUTINE ADVE(NTSD,DT,DETA1,DETA2,PDTOP &
45 & ,CURV,F,FAD,F4D,EM_LOC,EMT_LOC,EN,ENT,DX,DY &
46 & ,HBM2,VBM2 &
47 & ,T,U,V,PDSLO,TOLD,UOLD,VOLD &
48 & ,PETDT,UPSTRM &
49 & ,FEW,FNS,FNE,FSE &
50 & ,ADT,ADU,ADV &
51 & ,N_IUP_H,N_IUP_V &
52 & ,N_IUP_ADH,N_IUP_ADV &
53 & ,IUP_H,IUP_V,IUP_ADH,IUP_ADV &
54 & ,IHE,IHW,IVE,IVW &
55 & ,IDS,IDE,JDS,JDE,KDS,KDE &
56 & ,IMS,IME,JMS,JME,KMS,KME &
57 & ,ITS,ITE,JTS,JTE,KTS,KTE)
58 !***********************************************************************
59 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
60 ! . . .
61 ! SUBPROGRAM: ADVE HORIZONTAL AND VERTICAL ADVECTION
62 ! PRGRMMR: JANJIC ORG: W/NP22 DATE: 93-10-28
63 !
64 ! ABSTRACT:
65 ! ADVE CALCULATES THE CONTRIBUTION OF THE HORIZONTAL AND VERTICAL
66 ! ADVECTION TO THE TENDENCIES OF TEMPERATURE AND WIND AND THEN
67 ! UPDATES THOSE VARIABLES.
68 ! THE JANJIC ADVECTION SCHEME FOR THE ARAKAWA E GRID IS USED
69 ! FOR ALL VARIABLES INSIDE THE FIFTH ROW. AN UPSTREAM SCHEME
70 ! IS USED ON ALL VARIABLES IN THE THIRD, FOURTH, AND FIFTH
71 ! OUTERMOST ROWS. THE ADAMS-BASHFORTH TIME SCHEME IS USED.
72 !
73 ! PROGRAM HISTORY LOG:
74 ! 87-06-?? JANJIC - ORIGINATOR
75 ! 95-03-25 BLACK - CONVERSION FROM 1-D TO 2-D IN HORIZONTAL
76 ! 96-03-28 BLACK - ADDED EXTERNAL EDGE
77 ! 98-10-30 BLACK - MODIFIED FOR DISTRIBUTED MEMORY
78 ! 99-07- JANJIC - CONVERTED TO ADAMS-BASHFORTH SCHEME
79 ! COMBINING HORIZONTAL AND VERTICAL ADVECTION
80 ! 02-02-04 BLACK - ADDED VERTICAL CFL CHECK
81 ! 02-02-05 BLACK - CONVERTED TO WRF FORMAT
82 ! 02-08-29 MICHALAKES - CONDITIONAL COMPILATION OF MPI
83 ! CONVERT TO GLOBAL INDEXING
84 ! 02-09-06 WOLFE - MORE CONVERSION TO GLOBAL INDEXING
85 ! 04-05-29 JANJIC,BLACK - CRANK-NICHOLSON VERTICAL ADVECTION
86 ! 04-11-23 BLACK - THREADED
87 ! 05-12-14 BLACK - CONVERTED FROM IKJ TO IJK
88 !
89 ! USAGE: CALL ADVE FROM SUBROUTINE SOLVE_NMM
90 ! INPUT ARGUMENT LIST:
91 !
92 ! OUTPUT ARGUMENT LIST:
93 !
94 ! OUTPUT FILES:
95 ! NONE
96 !
97 ! SUBPROGRAMS CALLED:
98 !
99 ! UNIQUE: NONE
100 !
101 ! LIBRARY: NONE
102 !
103 ! ATTRIBUTES:
104 ! LANGUAGE: FORTRAN 90
105 ! MACHINE : IBM SP
106 !$$$
107 !***********************************************************************
108 !-----------------------------------------------------------------------
109 !
110 IMPLICIT NONE
111 !
112 !-----------------------------------------------------------------------
113 !
114 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
115 & ,IMS,IME,JMS,JME,KMS,KME &
116 & ,ITS,ITE,JTS,JTE,KTS,KTE
117 !
118 INTEGER, DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW &
119 ,N_IUP_H,N_IUP_V &
120 & ,N_IUP_ADH,N_IUP_ADV
121 !
122 INTEGER, DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V &
123 & ,IUP_ADH,IUP_ADV
124 !
125 INTEGER,INTENT(IN) :: NTSD
126 !
127 REAL,INTENT(IN) :: DT,DY,EN,ENT,F4D,PDTOP
128 !
129 REAL,DIMENSION(NMM_MAX_DIM),INTENT(IN) :: EM_LOC,EMT_LOC
130 !
131 REAL,DIMENSION(KMS:KME),INTENT(IN) :: DETA1,DETA2
132 !
133 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: CURV,DX,F,FAD,HBM2 &
134 & ,PDSLO,VBM2
135 !
136 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(OUT) :: ADT,ADU,ADV
137 !
138 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: PETDT
139 !
140 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(INOUT) :: T,TOLD &
141 & ,U,UOLD &
142 & ,V,VOLD
143 !
144 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(OUT) :: FEW,FNE &
145 & ,FNS,FSE
146 !
147 !-----------------------------------------------------------------------
148 !*** LOCAL VARIABLES
149 !-----------------------------------------------------------------------
150 !
151 LOGICAL :: UPSTRM
152 !
153 INTEGER :: I,IEND,IFP,IFQ,II,IPQ,ISP,ISQ,ISTART &
154 & ,IUP_ADH_J,IVH,IVL &
155 & ,J,J1,JA,JAK,JEND,JGLOBAL,JJ,JKNT,JP2,JSTART &
156 & ,K,KNTI_ADH,KSTART,KSTOP &
157 & ,N,N_IUPH_J,N_IUPADH_J,N_IUPADV_J
158 !
159 INTEGER :: MY_IS_GLB,MY_IE_GLB,MY_JS_GLB,MY_JE_GLB
160 !
161 INTEGER,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ISPA,ISQA
162 !
163 REAL :: ADPDX,ADPDY,ARRAY3_X,CFL,CFT,CFU,CFV,CMT,CMU,CMV &
164 & ,DTE,DTQ,F0,F1,F2,F3,FEWP,FNEP,FNSP,FPP,FSEP,HM &
165 & ,PDOP,PDOPU,PDOPV,PP &
166 & ,PVVLO,PVVLOU,PVVLOV,PVVUP,PVVUPU,PVVUPV &
167 & ,QP,RDP,RDPU,RDPV &
168 & ,TEMPA,TEMPB,TTA,TTB,UDY &
169 & ,VDX,VM,VVLO,VVLOU,VVLOV,VVUP,VVUPU,VVUPV
170 !
171 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ARRAY0,ARRAY1 &
172 & ,ARRAY2,ARRAY3 &
173 & ,DPDE,RDPD,RDPDX,RDPDY &
174 & ,TEW,TNE,TNS,TSE,TST &
175 & ,UNE,UNED,UEW,UNS,USE &
176 & ,USED,UST &
177 & ,VEW,VNE,VNS,VSE &
178 & ,VST
179 !
180 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5,KTS:KTE) :: VAD_TEND_T &
181 & ,VAD_TEND_U &
182 & ,VAD_TEND_V
183 !
184 REAL,DIMENSION(KTS:KTE) :: CRT,CRU,CRV,DETA1_PDTOP &
185 & ,RCMT,RCMU,RCMV,RSTT,RSTU,RSTV &
186 & ,T_K,TN,U_K,UN,V_K,VN
187 !
188 !-----------------------------------------------------------------------
189 !***********************************************************************
190 !
191 ! DPDE ----- 3
192 ! | J Increasing
193 ! |
194 ! | ^
195 ! FNS ----- 2 |
196 ! | |
197 ! | |
198 ! | |
199 ! VNS ----- 1 |
200 ! |
201 ! |
202 ! |
203 ! ADV ----- 0 ------> Current J
204 ! |
205 ! |
206 ! |
207 ! VNS ----- -1
208 ! |
209 ! |
210 ! |
211 ! FNS ----- -2
212 ! |
213 ! |
214 ! |
215 ! DPDE ----- -3
216 !
217 !***********************************************************************
218 !-----------------------------------------------------------------------
219 DO J=JTS-5,JTE+5
220 DO I=ITS-5,ITE+5
221 ARRAY0(I,J)=0.0
222 ARRAY1(I,J)=0.0
223 ARRAY2(I,J)=0.0
224 ARRAY3(I,J)=0.0
225 DPDE(I,J)=0.0
226 RDPD(I,J)=0.0
227 RDPDX(I,J)=0.0
228 RDPDY(I,J)=0.0
229 TEW(I,J)=0.0
230 TNE(I,J)=0.0
231 TNS(I,J)=0.0
232 TSE(I,J)=0.0
233 TST(I,J)=0.0
234 UNE(I,J)=0.0
235 UNED(I,J)=0.0
236 UEW(I,J)=0.0
237 UNS(I,J)=0.0
238 USE(I,J)=0.0
239 USED(I,J)=0.0
240 UST(I,J)=0.0
241 VEW(I,J)=0.0
242 VNE(I,J)=0.0
243 VNS(I,J)=0.0
244 VSE(I,J)=0.0
245 VST(I,J)=0.0
246 ENDDO
247 ENDDO
248 !-----------------------------------------------------------------------
249 !
250 DTQ=DT*0.25
251 DTE=DT*(0.5*0.25)
252 !
253 !-----------------------------------------------------------------------
254 !***
255 !*** PRECOMPUTE DETA1 TIMES PDTOP.
256 !***
257 !-----------------------------------------------------------------------
258 !
259 DO K=KTS,KTE
260 DETA1_PDTOP(K)=DETA1(K)*PDTOP
261 ENDDO
262 !
263 !-----------------------------------------------------------------------
264 !***
265 !*** INITIALIZE SOME WORKING ARRAYS TO ZERO
266 !***
267 !
268 !-----------------------------------------------------------------------
269 !-----------------------------------------------------------------------
270 !
271 !*** COMPUTE VERTICAL ADVECTION TENDENCIES USING CRANK-NICHOLSON.
272 !
273 !-----------------------------------------------------------------------
274 !-----------------------------------------------------------------------
275 !
276 !-----------------------------------------------------------------------
277 !*** FIRST THE TEMPERATURE
278 !-----------------------------------------------------------------------
279 !$omp parallel do &
280 !$omp& private(cft,cfu,cfv,cmt,cmu,cmv,crt,cru,crv,i,k &
281 !$omp& ,pdop,pdopu,pdopv,pvvlo,pvvlou,pvvlov,pvvup,pvvupu,pvvupv &
282 !$omp& ,rcmt,rcmu,rcmv,rdp,rdpu,rdpv,rstt,rstu,rstv,t_k,tn &
283 !$omp& ,u_k,un,v_k,vn,vvlo,vvlou,vvlov,vvup,vvupu,vvupv)
284 !!$omp& private(adtp,adup,advp,ttlo,ttup,tulo,tuup,tvlo,tvup)
285 !-----------------------------------------------------------------------
286 !
287 main_vertical: DO J=MYJS2,MYJE2
288 !
289 !-----------------------------------------------------------------------
290 !
291 iloop_for_t: DO I=MYIS1,MYIE1
292 !
293 !-----------------------------------------------------------------------
294 !*** EXTRACT T FROM THE COLUMN
295 !-----------------------------------------------------------------------
296 !
297 DO K=KTS,KTE
298 T_K(K)=T(I,J,K)
299 ENDDO
300 !
301 !-----------------------------------------------------------------------
302 !
303 PDOP=PDSLO(I,J)
304 PVVLO=PETDT(I,J,KTE-1)*DTQ
305 VVLO=PVVLO/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOP)
306 CMT=-VVLO*WGT2+1.
307 RCMT(KTE)=1./CMT
308 CRT(KTE)=VVLO*WGT2
309 RSTT(KTE)=-VVLO*WGT1*(T_K(KTE-1)-T_K(KTE))+T_K(KTE)
310 !
311 !-----------------------------------------------------------------------
312 !
313 DO K=KTE-1,KTS+1,-1
314 RDP=1./(DETA1_PDTOP(K)+DETA2(K)*PDOP)
315 PVVUP=PVVLO
316 PVVLO=PETDT(I,J,K-1)*DTQ
317 VVUP=PVVUP*RDP
318 VVLO=PVVLO*RDP
319 CFT=-VVUP*WGT2*RCMT(K+1)
320 CMT=-CRT(K+1)*CFT+((VVUP-VVLO)*WGT2+1.)
321 RCMT(K)=1./CMT
322 CRT(K)=VVLO*WGT2
323 RSTT(K)=-RSTT(K+1)*CFT+T_K(K) &
324 & -(T_K(K)-T_K(K+1))*VVUP*WGT1 &
325 & -(T_K(K-1)-T_K(K))*VVLO*WGT1
326 ENDDO
327 !
328 !-----------------------------------------------------------------------
329 !
330 PVVUP=PVVLO
331 VVUP=PVVUP/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOP)
332 CFT=-VVUP*WGT2*RCMT(KTS+1)
333 CMT=-CRT(KTS+1)*CFT+VVUP*WGT2+1.
334 CRT(KTS)=0.
335 RSTT(KTS)=-(T_K(KTS)-T_K(KTS+1))*VVUP*WGT1 &
336 & -RSTT(KTS+1)*CFT+T_K(KTS)
337 TN(KTS)=RSTT(KTS)/CMT
338 VAD_TEND_T(I,J,KTS)=TN(KTS)-T_K(KTS)
339 !
340 DO K=KTS+1,KTE
341 TN(K)=(-CRT(K)*TN(K-1)+RSTT(K))*RCMT(K)
342 VAD_TEND_T(I,J,K)=TN(K)-T_K(K)
343 ENDDO
344 !
345 !-----------------------------------------------------------------------
346 !*** The following section is only for checking the implicit solution
347 !*** using back-substitution. Remove this section otherwise.
348 !-----------------------------------------------------------------------
349 ! if(ntsd<=10.or.ntsd>=6000)then
350 ! IF(I==ITEST.AND.J==JTEST)THEN
351 !!
352 ! PVVLO=PETDT(I,J,KTE-1)*DT*0.25
353 ! VVLO=PVVLO/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOP)
354 ! TTLO=VVLO*(T(I,J,KTE-1)-T(I,J,KTE) &
355 ! & +TN(KTE-1)-TN(KTE))
356 ! ADTP=TTLO+TN(KTE)-T(I,J,KTE)
357 ! WRITE(0,*)' NTSD=',NTSD,' I=',ITEST,' J=',JTEST,' K=',KTE &
358 ! &, ' ADTP=',ADTP
359 ! WRITE(0,*)' T=',T(I,J,KTE),' TN=',TN(KTE) &
360 ! &, ' VAD_TEND_T=',VAD_TEND_T(I,J,KTE)
361 ! WRITE(0,*)' '
362 !!
363 ! DO K=KTE-1,KTS+1,-1
364 ! RDP=1./(DETA1_PDTOP(K)+DETA2(K)*PDOP)
365 ! PVVUP=PVVLO
366 ! PVVLO=PETDT(I,J,K-1)*DT*0.25
367 ! VVUP=PVVUP*RDP
368 ! VVLO=PVVLO*RDP
369 ! TTUP=VVUP*(T(I,J,K)-T(I,J,K+1)+TN(K)-TN(K+1))
370 ! TTLO=VVLO*(T(I,J,K-1)-T(I,J,K)+TN(K-1)-TN(K))
371 ! ADTP=TTLO+TTUP+TN(K)-T(I,J,K)
372 ! WRITE(0,*)' NTSD=',NTSD,' I=',I,' J=',J,' K=',K &
373 ! &, ' ADTP=',ADTP
374 ! WRITE(0,*)' T=',T(I,J,K),' TN=',TN(K) &
375 ! &, ' VAD_TEND_T=',VAD_TEND_T(I,J,K)
376 ! WRITE(0,*)' '
377 ! ENDDO
378 !!
379 ! PVVUP=PVVLO
380 ! VVUP=PVVUP/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOP)
381 ! TTUP=VVUP*(T(I,J,KTS)-T(I,J,KTS+1)+TN(KTS)-TN(KTS+1))
382 ! ADTP=TTUP+TN(KTS)-T(I,J,KTS)
383 ! WRITE(0,*)' NTSD=',NTSD,' I=',I,' J=',J,' K=',KTS &
384 ! &, ' ADTP=',ADTP
385 ! WRITE(0,*)' T=',T(I,J,KTS),' TN=',TN(KTS) &
386 ! &, ' VAD_TEND_T=',VAD_TEND_T(I,J,KTS)
387 ! WRITE(0,*)' '
388 ! ENDIF
389 ! endif
390 !
391 !-----------------------------------------------------------------------
392 !*** End of check.
393 !-----------------------------------------------------------------------
394 !
395 ENDDO iloop_for_t
396 !
397 !-----------------------------------------------------------------------
398 !
399 !*** NOW VERTICAL ADVECTION OF WIND COMPONENTS
400 !
401 !-----------------------------------------------------------------------
402 !
403 iloop_for_uv: DO I=MYIS1,MYIE1
404 !
405 !-----------------------------------------------------------------------
406 !*** EXTRACT U AND V FROM THE COLUMN
407 !-----------------------------------------------------------------------
408 !
409 DO K=KTS,KTE
410 U_K(K)=U(I,J,K)
411 V_K(K)=V(I,J,K)
412 ENDDO
413 !
414 !-----------------------------------------------------------------------
415 !
416 PDOPU=(PDSLO(I+IVW(J),J)+PDSLO(I+IVE(J),J))*0.5
417 PDOPV=(PDSLO(I,J-1)+PDSLO(I,J+1))*0.5
418 PVVLOU=(PETDT(I+IVW(J),J,KTE-1)+PETDT(I+IVE(J),J,KTE-1))*DTE
419 PVVLOV=(PETDT(I,J-1,KTE-1)+PETDT(I,J+1,KTE-1))*DTE
420 VVLOU=PVVLOU/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPU)
421 VVLOV=PVVLOV/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPV)
422 CMU=-VVLOU*WGT2+1.
423 CMV=-VVLOV*WGT2+1.
424 RCMU(KTE)=1./CMU
425 RCMV(KTE)=1./CMV
426 CRU(KTE)=VVLOU*WGT2
427 CRV(KTE)=VVLOV*WGT2
428 RSTU(KTE)=-VVLOU*WGT1*(U_K(KTE-1)-U_K(KTE))+U_K(KTE)
429 RSTV(KTE)=-VVLOV*WGT1*(V_K(KTE-1)-V_K(KTE))+V_K(KTE)
430 !
431 !-----------------------------------------------------------------------
432 !
433 DO K=KTE-1,KTS+1,-1
434 RDPU=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPU)
435 RDPV=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPV)
436 PVVUPU=PVVLOU
437 PVVUPV=PVVLOV
438 PVVLOU=(PETDT(I+IVW(J),J,K-1)+PETDT(I+IVE(J),J,K-1))*DTE
439 PVVLOV=(PETDT(I,J-1,K-1)+PETDT(I,J+1,K-1))*DTE
440 VVUPU=PVVUPU*RDPU
441 VVUPV=PVVUPV*RDPV
442 VVLOU=PVVLOU*RDPU
443 VVLOV=PVVLOV*RDPV
444 CFU=-VVUPU*WGT2*RCMU(K+1)
445 CFV=-VVUPV*WGT2*RCMV(K+1)
446 CMU=-CRU(K+1)*CFU+(VVUPU-VVLOU)*WGT2+1.
447 CMV=-CRV(K+1)*CFV+(VVUPV-VVLOV)*WGT2+1.
448 RCMU(K)=1./CMU
449 RCMV(K)=1./CMV
450 CRU(K)=VVLOU*WGT2
451 CRV(K)=VVLOV*WGT2
452 RSTU(K)=-RSTU(K+1)*CFU+U_K(K) &
453 & -(U_K(K)-U_K(K+1))*VVUPU*WGT1 &
454 & -(U_K(K-1)-U_K(K))*VVLOU*WGT1
455 RSTV(K)=-RSTV(K+1)*CFV+V_K(K) &
456 & -(V_K(K)-V_K(K+1))*VVUPV*WGT1 &
457 & -(V_K(K-1)-V_K(K))*VVLOV*WGT1
458 ENDDO
459 !
460 !-----------------------------------------------------------------------
461 !
462 RDPU=1./(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPU)
463 RDPV=1./(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPV)
464 PVVUPU=PVVLOU
465 PVVUPV=PVVLOV
466 VVUPU=PVVUPU*RDPU
467 VVUPV=PVVUPV*RDPV
468 CFU=-VVUPU*WGT2*RCMU(KTS+1)
469 CFV=-VVUPV*WGT2*RCMV(KTS+1)
470 CMU=-CRU(KTS+1)*CFU+VVUPU*WGT2+1.
471 CMV=-CRV(KTS+1)*CFV+VVUPV*WGT2+1.
472 CRU(KTS)=0.
473 CRV(KTS)=0.
474 RSTU(KTS)=-(U_K(KTS)-U_K(KTS+1))*VVUPU*WGT1 &
475 & -RSTU(KTS+1)*CFU+U_K(KTS)
476 RSTV(KTS)=-(V_K(KTS)-V_K(KTS+1))*VVUPV*WGT1 &
477 & -RSTV(KTS+1)*CFV+V_K(KTS)
478 UN(KTS)=RSTU(KTS)/CMU
479 VN(KTS)=RSTV(KTS)/CMV
480 VAD_TEND_U(I,J,KTS)=UN(KTS)-U_K(KTS)
481 VAD_TEND_V(I,J,KTS)=VN(KTS)-V_K(KTS)
482 !
483 DO K=KTS+1,KTE
484 UN(K)=(-CRU(K)*UN(K-1)+RSTU(K))*RCMU(K)
485 VN(K)=(-CRV(K)*VN(K-1)+RSTV(K))*RCMV(K)
486 VAD_TEND_U(I,J,K)=UN(K)-U_K(K)
487 VAD_TEND_V(I,J,K)=VN(K)-V_K(K)
488 ENDDO
489 !
490 !-----------------------------------------------------------------------
491 !*** The following section is only for checking the implicit solution
492 !*** using back-substitution. Remove this section otherwise.
493 !-----------------------------------------------------------------------
494 !
495 ! if(ntsd<=10.or.ntsd>=6000)then
496 ! IF(I==ITEST.AND.J==JTEST)THEN
497 !!
498 ! PDOPU=(PDSLO(I+IVW(J),J)+PDSLO(I+IVE(J),J))*0.5
499 ! PDOPV=(PDSLO(I,J-1)+PDSLO(I,J+1))*0.5
500 ! PVVLOU=(PETDT(I+IVW(J),J,KTE-1) &
501 ! & +PETDT(I+IVE(J),J,KTE-1))*DTE
502 ! PVVLOV=(PETDT(I,J-1,KTE-1) &
503 ! & +PETDT(I,J+1,KTE-1))*DTE
504 ! VVLOU=PVVLOU/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPU)
505 ! VVLOV=PVVLOV/(DETA1_PDTOP(KTE)+DETA2(KTE)*PDOPV)
506 ! TULO=VVLOU*(U(I,J,KTE-1)-U(I,J,KTE)+UN(KTE-1)-UN(KTE))
507 ! TVLO=VVLOV*(V(I,J,KTE-1)-V(I,J,KTE)+VN(KTE-1)-VN(KTE))
508 ! ADUP=TULO+UN(KTE)-U(I,J,KTE)
509 ! ADVP=TVLO+VN(KTE)-V(I,J,KTE)
510 ! WRITE(0,*)' NTSD=',NTSD,' I=',I,' J=',J,' K=',KTE &
511 ! &, ' ADUP=',ADUP,' ADVP=',ADVP
512 ! WRITE(0,*)' U=',U(I,J,KTE),' UN=',UN(KTE) &
513 ! &, ' VAD_TEND_U=',VAD_TEND_U(I,KTE) &
514 ! &, ' V=',V(I,J,KTE),' VN=',VN(KTE) &
515 ! &, ' VAD_TEND_V=',VAD_TEND_V(I,KTE)
516 ! WRITE(0,*)' '
517 !!
518 ! DO K=KTE-1,KTS+1,-1
519 ! RDPU=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPU)
520 ! RDPV=1./(DETA1_PDTOP(K)+DETA2(K)*PDOPV)
521 ! PVVUPU=PVVLOU
522 ! PVVUPV=PVVLOV
523 ! PVVLOU=(PETDT(I+IVW(J),J,K-1) &
524 ! & +PETDT(I+IVE(J),J,K-1))*DTE
525 ! PVVLOV=(PETDT(I,J-1,K-1)+PETDT(I,J+1,K-1))*DTE
526 ! VVUPU=PVVUPU*RDPU
527 ! VVUPV=PVVUPV*RDPV
528 ! VVLOU=PVVLOU*RDPU
529 ! VVLOV=PVVLOV*RDPV
530 ! TUUP=VVUPU*(U(I,J,K)-U(I,J,K+1)+UN(K)-UN(K+1))
531 ! TVUP=VVUPV*(V(I,J,K)-V(I,J,K+1)+VN(K)-VN(K+1))
532 ! TULO=VVLOU*(U(I,J,K-1)-U(I,J,K)+UN(K-1)-UN(K))
533 ! TVLO=VVLOV*(V(I,J,K-1)-V(I,J,K)+VN(K-1)-VN(K))
534 ! ADUP=TUUP+TULO+UN(K)-U(I,J,K)
535 ! ADVP=TVUP+TVLO+VN(K)-V(I,J,K)
536 ! WRITE(0,*)' NTSD=',NTSD,' I=',ITEST,' J=',JTEST,' K=',K &
537 ! &, ' ADUP=',ADUP,' ADVP=',ADVP
538 ! WRITE(0,*)' U=',U(I,J,K),' UN=',UN(K) &
539 ! &, ' VAD_TEND_U=',VAD_TEND_U(I,K) &
540 ! &, ' V=',V(I,J,K),' VN=',VN(K) &
541 ! &, ' VAD_TEND_V=',VAD_TEND_V(I,K)
542 ! WRITE(0,*)' '
543 ! ENDDO
544 !!
545 ! PVVUPU=PVVLOU
546 ! PVVUPV=PVVLOV
547 ! VVUPU=PVVUPU/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPU)
548 ! VVUPV=PVVUPV/(DETA1_PDTOP(KTS)+DETA2(KTS)*PDOPV)
549 ! TUUP=VVUPU*(U(I,J,KTS)-U(I,J,KTS+1)+UN(KTS)-UN(KTS+1))
550 ! TVUP=VVUPV*(V(I,J,KTS)-V(I,J,KTS+1)+VN(KTS)-VN(KTS+1))
551 ! ADUP=TUUP+UN(KTS)-U(I,J,KTS)
552 ! ADVP=TVUP+VN(KTS)-V(I,J,KTS)
553 ! WRITE(0,*)' NTSD=',NTSD,' I=',ITEST,' J=',JTEST,' K=',KTS &
554 ! &, ' ADUP=',ADUP,' ADVP=',ADVP
555 ! WRITE(0,*)' U=',U(I,J,KTS),' UN=',UN(KTS) &
556 ! &, ' VAD_TEND_U=',VAD_TEND_U(I,KTS) &
557 ! &, ' V=',V(I,J,KTS),' VN=',VN(KTS) &
558 ! &, ' VAD_TEND_V=',VAD_TEND_V(I,KTS)
559 ! WRITE(0,*)' '
560 ! ENDIF
561 ! endif
562 !
563 !-----------------------------------------------------------------------
564 !*** End of check.
565 !-----------------------------------------------------------------------
566 !
567 ENDDO iloop_for_uv
568 !
569 !-----------------------------------------------------------------------
570 !
571 ENDDO main_vertical
572 !
573 !-----------------------------------------------------------------------
574 !-----------------------------------------------------------------------
575 !
576 !*** COMPUTE HORIZONTAL ADVECTION TENDENCIES.
577 !
578 !-----------------------------------------------------------------------
579 !-----------------------------------------------------------------------
580 !$omp parallel do &
581 !$omp& private(adpdx,adpdy,adt,adu,adv,array0,array1,array2,array3 &
582 !$omp& ,array3_x,dpde,f0,f1,f2,f3,fewp,fnep,fnsp,fpp,fsep,hm &
583 !$omp& ,i,ifp,ifq,ii,ipq,isp,ispa,isq,isqa,iup_adh_j,j,k &
584 !$omp& ,knti_adh,n_iupadh_j,n_iupadv_j,n_iuph_j,pp,qp &
585 !$omp& ,rdpd,rdpdx,rdpdy,tew,tne,tns,tse,tst,tta,ttb &
586 !$omp& ,uew,udy,une,uned,uns,use,used,ust &
587 !$omp& ,vdx,vew,vm,vne,vns,vse,vst)
588 !-----------------------------------------------------------------------
589 !
590 main_horizontal: DO K=KTS,KTE
591 !
592 !-----------------------------------------------------------------------
593 !
594 DO J=MYJS_P4,MYJE_P4
595 DO I=MYIS_P4,MYIE_P4
596 DPDE(I,J)=DETA1_PDTOP(K)+DETA2(K)*PDSLO(I,J)
597 RDPD(I,J)=1./DPDE(I,J)
598 TST(I,J)=T(I,J,K)*FFC+TOLD(I,J,K)*FBC
599 UST(I,J)=U(I,J,K)*FFC+UOLD(I,J,K)*FBC
600 VST(I,J)=V(I,J,K)*FFC+VOLD(I,J,K)*FBC
601 ENDDO
602 ENDDO
603 !
604 !-----------------------------------------------------------------------
605 !*** MASS FLUXES AND MASS POINT ADVECTION COMPONENTS
606 !*** THE NS AND EW FLUXES IN THE FOLLOWING LOOP ARE ON V POINTS
607 !*** FOR T.
608 !-----------------------------------------------------------------------
609 !
610 DO J=MYJS1_P3,MYJE1_P3
611 DO I=MYIS_P3,MYIE_P3
612 !
613 ADPDX=DPDE(I+IVW(J),J)+DPDE(I+IVE(J),J)
614 ADPDY=DPDE(I,J-1)+DPDE(I,J+1)
615 RDPDX(I,J)=1./ADPDX
616 RDPDY(I,J)=1./ADPDY
617 !
618 UDY=U(I,J,K)*DY
619 VDX=V(I,J,K)*DX(I,J)
620 !
621 FEWP=UDY*ADPDX
622 FNSP=VDX*ADPDY
623 !
624 FEW(I,J,K)=FEWP
625 FNS(I,J,K)=FNSP
626 !
627 TEW(I,J)=FEWP*(TST(I+IVE(J),J)-TST(I+IVW(J),J))
628 TNS(I,J)=FNSP*(TST(I,J+1)-TST(I,J-1))
629 !
630 UNED(I,J)=UDY+VDX
631 USED(I,J)=UDY-VDX
632 !
633 ENDDO
634 ENDDO
635 !
636 !-----------------------------------------------------------------------
637 !*** DIAGONAL FLUXES AND DIAGONALLY AVERAGED WIND
638 !*** THE NE AND SE FLUXES ARE ASSOCIATED WITH H POINTS
639 !*** (ACTUALLY JUST TO THE NE AND SE OF EACH H POINT).
640 !-----------------------------------------------------------------------
641 !
642 DO J=MYJS1_P2,MYJE2_P2
643 DO I=MYIS_P2,MYIE_P2
644 FNEP=(UNED(I+IHE(J),J)+UNED(I ,J+1)) &
645 & *(DPDE(I ,J)+DPDE(I+IHE(J),J+1))
646 FNE(I,J,K)=FNEP
647 TNE(I,J)=FNEP*(TST(I+IHE(J),J+1)-TST(I,J))
648 ENDDO
649 ENDDO
650 !
651 DO J=MYJS2_P2,MYJE1_P2
652 DO I=MYIS_P2,MYIE_P2
653 FSEP=(USED(I+IHE(J),J)+USED(I ,J-1)) &
654 & *(DPDE(I ,J)+DPDE(I+IHE(J),J-1))
655 FSE(I,J,K)=FSEP
656 TSE(I,J)=FSEP*(TST(I+IHE(J),J-1)-TST(I,J))
657 !
658 ENDDO
659 ENDDO
660 !
661 !-----------------------------------------------------------------------
662 !*** HORIZONTAL T ADVECTION TENDENCY ADT IS ON H POINTS OF COURSE.
663 !-----------------------------------------------------------------------
664 !
665 DO J=MYJS5,MYJE5
666 DO I=MYIS2,MYIE2
667 ADT(I,J)=(TEW(I+IHW(J),J)+TEW(I+IHE(J),J) &
668 & +TNS(I,J-1)+TNS(I,J+1) &
669 & +TNE(I+IHW(J),J-1)+TNE(I,J) &
670 & +TSE(I,J)+TSE(I+IHW(J),J+1)) &
671 & *RDPD(I,J)*FAD(I,J)
672 ENDDO
673 ENDDO
674 !
675 !
676 !-----------------------------------------------------------------------
677 !*** CALCULATION OF MOMENTUM ADVECTION COMPONENTS.
678 !-----------------------------------------------------------------------
679 !
680 DO J=MYJS4_P1,MYJE4_P1
681 DO I=MYIS_P1,MYIE_P1
682 !
683 !-----------------------------------------------------------------------
684 !*** THE NS AND EW FLUXES ARE ON H POINTS FOR U AND V.
685 !-----------------------------------------------------------------------
686 !
687 UEW(I,J)=(FEW(I+IHW(J),J,K)+FEW(I+IHE(J),J,K)) &
688 & *(UST(I+IHE(J),J)-UST(I+IHW(J),J))
689 UNS(I,J)=(FNS(I+IHW(J),J,K)+FNS(I+IHE(J),J,K)) &
690 & *(UST(I,J+1)-UST(I,J-1))
691 VEW(I,J)=(FEW(I,J-1,K)+FEW(I,J+1,K)) &
692 & *(VST(I+IHE(J),J)-VST(I+IHW(J),J))
693 VNS(I,J)=(FNS(I,J-1,K)+FNS(I,J+1,K)) &
694 & *(VST(I,J+1)-VST(I,J-1))
695 !
696 !-----------------------------------------------------------------------
697 !*** THE FOLLOWING NE AND SE FLUXES ARE TIED TO V POINTS AND ARE
698 !*** LOCATED JUST TO THE NE AND SE OF THE GIVEN I,J.
699 !-----------------------------------------------------------------------
700 !
701 UNE(I,J)=(FNE(I+IVW(J),J,K)+FNE(I+IVE(J),J,K)) &
702 & *(UST(I+IVE(J),J+1)-UST(I,J))
703 USE(I,J)=(FSE(I+IVW(J),J,K)+FSE(I+IVE(J),J,K)) &
704 & *(UST(I+IVE(J),J-1)-UST(I,J))
705 VNE(I,J)=(FNE(I,J-1,K)+FNE(I,J+1,K)) &
706 & *(VST(I+IVE(J),J+1)-VST(I,J))
707 VSE(I,J)=(FSE(I,J-1,K)+FSE(I,J+1,K)) &
708 & *(VST(I+IVE(J),J-1)-VST(I,J))
709 !
710 !-----------------------------------------------------------------------
711 !
712 ENDDO
713 ENDDO
714 !
715 !-----------------------------------------------------------------------
716 !*** COMPUTE THE ADVECTION TENDENCIES FOR U AND V.
717 !*** THE AD ARRAYS ARE ON THE VELOCITY POINTS.
718 !-----------------------------------------------------------------------
719 !
720 DO J=MYJS5,MYJE5
721 DO I=MYIS2,MYIE2
722 ADU(I,J)=(UEW(I+IVW(J),J)+UEW(I+IVE(J),J) &
723 & +UNS(I,J-1)+UNS(I,J+1) &
724 & +UNE(I+IVW(J),J-1)+UNE(I,J) &
725 & +USE(I,J)+USE(I+IVW(J),J+1)) &
726 & *RDPDX(I,J)*FAD(I+IVW(J),J)
727 !
728 ADV(I,J)=(VEW(I+IVW(J),J)+VEW(I+IVE(J),J) &
729 & +VNS(I,J-1)+VNS(I,J+1) &
730 & +VNE(I+IVW(J),J-1)+VNE(I,J) &
731 & +VSE(I,J)+VSE(I+IVW(J),J+1)) &
732 & *RDPDY(I,J)*FAD(I+IVW(J),J)
733 ENDDO
734 ENDDO
735 !
736 !-----------------------------------------------------------------------
737 !
738 !*** END OF JANJIC HORIZONTAL ADVECTION
739 !
740 !-----------------------------------------------------------------------
741 !
742 !*** UPSTREAM ADVECTION OF T
743 !
744 !-----------------------------------------------------------------------
745 !
746 upstream: IF(UPSTRM)THEN
747 !
748 !-----------------------------------------------------------------------
749 !***
750 !*** COMPUTE UPSTREAM COMPUTATIONS ON THIS TASK'S ROWS.
751 !***
752 !-----------------------------------------------------------------------
753 !
754 jloop_upstream: DO J=MYJS2,MYJE2
755 !
756 N_IUPH_J=N_IUP_H(J) ! See explanation in START_DOMAIN_NMM
757 DO II=0,N_IUPH_J-1
758 !
759 I=IUP_H(IMS+II,J)
760 TTA=EMT_LOC(J)*(UST(I,J-1)+UST(I+IHW(J),J) &
761 & +UST(I+IHE(J),J)+UST(I,J+1))
762 TTB=ENT *(VST(I,J-1)+VST(I+IHW(J),J) &
763 & +VST(I+IHE(J),J)+VST(I,J+1))
764 PP=-TTA-TTB
765 QP= TTA-TTB
766 !
767 IF(PP<0.)THEN
768 ISPA(I,J)=-1
769 ELSE
770 ISPA(I,J)= 1
771 ENDIF
772 !
773 IF(QP<0.)THEN
774 ISQA(I,J)=-1
775 ELSE
776 ISQA(I,J)= 1
777 ENDIF
778 !
779 PP=ABS(PP)
780 QP=ABS(QP)
781 ARRAY3_X=PP*QP
782 ARRAY0(I,J)=ARRAY3_X-PP-QP
783 ARRAY1(I,J)=PP-ARRAY3_X
784 ARRAY2(I,J)=QP-ARRAY3_X
785 ARRAY3(I,J)=ARRAY3_X
786 ENDDO
787 !
788 !-----------------------------------------------------------------------
789 !
790 N_IUPADH_J=N_IUP_ADH(J)
791 KNTI_ADH=1
792 IUP_ADH_J=IUP_ADH(IMS,J)
793 !
794 iloop_T: DO II=0,N_IUPH_J-1
795 !
796 I=IUP_H(IMS+II,J)
797 !
798 ISP=ISPA(I,J)
799 ISQ=ISQA(I,J)
800 IFP=(ISP-1)/2
801 IFQ=(-ISQ-1)/2
802 IPQ=(ISP-ISQ)/2
803 !
804 !-----------------------------------------------------------------------
805 !
806 IF(I==IUP_ADH_J)THEN ! Upstream advection T tendencies
807 !
808 ISP=ISPA(I,J)
809 ISQ=ISQA(I,J)
810 IFP=(ISP-1)/2
811 IFQ=(-ISQ-1)/2
812 IPQ=(ISP-ISQ)/2
813 !
814 F0=ARRAY0(I,J)
815 F1=ARRAY1(I,J)
816 F2=ARRAY2(I,J)
817 F3=ARRAY3(I,J)
818 !
819 ADT(I,J)=F0*T(I,J,K) &
820 & +F1*T(I+IHE(J)+IFP,J+ISP,K) &
821 & +F2*T(I+IHE(J)+IFQ,J+ISQ,K) &
822 +F3*T(I+IPQ,J+ISP+ISQ,K)
823 !
824 !-----------------------------------------------------------------------
825 !
826 IF(KNTI_ADH<N_IUPADH_J)THEN
827 IUP_ADH_J=IUP_ADH(IMS+KNTI_ADH,J)
828 KNTI_ADH=KNTI_ADH+1
829 ENDIF
830 !
831 ENDIF ! End of upstream advection T tendency IF block
832 !
833 ENDDO iloop_T
834 !
835 !-----------------------------------------------------------------------
836 !
837 !*** UPSTREAM ADVECTION OF VELOCITY COMPONENTS
838 !
839 !-----------------------------------------------------------------------
840 !
841 N_IUPADV_J=N_IUP_ADV(J)
842 !
843 DO II=0,N_IUPADV_J-1
844 I=IUP_ADV(IMS+II,J)
845 !
846 TTA=EM_LOC(J)*UST(I,J)
847 TTB=EN *VST(I,J)
848 PP=-TTA-TTB
849 QP=TTA-TTB
850 !
851 IF(PP<0.)THEN
852 ISP=-1
853 ELSE
854 ISP= 1
855 ENDIF
856 !
857 IF(QP<0.)THEN
858 ISQ=-1
859 ELSE
860 ISQ= 1
861 ENDIF
862 !
863 IFP=(ISP-1)/2
864 IFQ=(-ISQ-1)/2
865 IPQ=(ISP-ISQ)/2
866 PP=ABS(PP)
867 QP=ABS(QP)
868 F3=PP*QP
869 F0=F3-PP-QP
870 F1=PP-F3
871 F2=QP-F3
872 !
873 ADU(I,J)=F0*U(I,J,K) &
874 & +F1*U(I+IVE(J)+IFP,J+ISP,K) &
875 & +F2*U(I+IVE(J)+IFQ,J+ISQ,K) &
876 & +F3*U(I+IPQ,J+ISP+ISQ,K)
877 !
878 ADV(I,J)=F0*V(I,J,K) &
879 & +F1*V(I+IVE(J)+IFP,J+ISP,K) &
880 & +F2*V(I+IVE(J)+IFQ,J+ISQ,K) &
881 & +F3*V(I+IPQ,J+ISP+ISQ,K)
882 !
883 ENDDO
884 !
885 ENDDO jloop_upstream
886 !
887 !-----------------------------------------------------------------------
888 !
889 ENDIF upstream
890 !
891 !-----------------------------------------------------------------------
892 !
893 !*** END OF HORIZONTAL ADVECTION
894 !
895 !-----------------------------------------------------------------------
896 !
897 !*** NOW SUM THE VERTICAL AND HORIZONTAL TENDENCIES,
898 !*** CURVATURE AND CORIOLIS TERMS.
899 !
900 !-----------------------------------------------------------------------
901 !
902 DO J=MYJS2,MYJE2
903 DO I=MYIS1,MYIE1
904 HM=HBM2(I,J)
905 VM=VBM2(I,J)
906 ADT(I,J)=(VAD_TEND_T(I,J,K)+2.*ADT(I,J))*HM
907 !
908 FPP=CURV(I,J)*2.*UST(I,J)+F(I,J)*2.
909 ADU(I,J)=(VAD_TEND_U(I,J,K)+2.*ADU(I,J)+VST(I,J)*FPP)*VM
910 ADV(I,J)=(VAD_TEND_V(I,J,K)+2.*ADV(I,J)-UST(I,J)*FPP)*VM
911 ENDDO
912 ENDDO
913 !
914 !-----------------------------------------------------------------------
915 !*** SAVE THE OLD VALUES FOR TIMESTEPPING
916 !-----------------------------------------------------------------------
917 !
918 DO J=MYJS_P4,MYJE_P4
919 DO I=MYIS_P4,MYIE_P4
920 TOLD(I,J,K)=T(I,J,K)
921 UOLD(I,J,K)=U(I,J,K)
922 VOLD(I,J,K)=V(I,J,K)
923 ENDDO
924 ENDDO
925 !
926 !-----------------------------------------------------------------------
927 !*** FINALLY UPDATE THE PROGNOSTIC VARIABLES
928 !-----------------------------------------------------------------------
929 !
930 DO J=MYJS2,MYJE2
931 DO I=MYIS1,MYIE1
932 T(I,J,K)=ADT(I,J)+T(I,J,K)
933 U(I,J,K)=ADU(I,J)+U(I,J,K)
934 V(I,J,K)=ADV(I,J)+V(I,J,K)
935 ENDDO
936 ENDDO
937 !
938 !-----------------------------------------------------------------------
939 !
940 ENDDO main_horizontal
941 !
942 !-----------------------------------------------------------------------
943 !
944 END SUBROUTINE ADVE
945 !
946 !-----------------------------------------------------------------------
947 !
948 !***********************************************************************
949 SUBROUTINE VAD2(NTSD,DT,IDTAD,DX,DY &
950 & ,AETA1,AETA2,DETA1,DETA2,PDSL,PDTOP,HBM2 &
951 & ,Q,Q2,CWM,PETDT &
952 & ,N_IUP_H,N_IUP_V &
953 & ,N_IUP_ADH,N_IUP_ADV &
954 & ,IUP_H,IUP_V,IUP_ADH,IUP_ADV &
955 & ,IHE,IHW,IVE,IVW &
956 & ,IDS,IDE,JDS,JDE,KDS,KDE &
957 & ,IMS,IME,JMS,JME,KMS,KME &
958 & ,ITS,ITE,JTS,JTE,KTS,KTE)
959 !***********************************************************************
960 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
961 ! . . .
962 ! SUBPROGRAM: VAD2 VERTICAL ADVECTION OF H2O SUBSTANCE AND TKE
963 ! PRGRMMR: JANJIC ORG: W/NP22 DATE: 96-07-19
964 !
965 ! ABSTRACT:
966 ! VAD2 CALCULATES THE CONTRIBUTION OF THE VERTICAL ADVECTION
967 ! TO THE TENDENCIES OF WATER SUBSTANCE AND TKE AND THEN UPDATES
968 ! THOSE VARIABLES. AN ANTI-FILTERING TECHNIQUE IS USED.
969 !
970 ! PROGRAM HISTORY LOG:
971 ! 96-07-19 JANJIC - ORIGINATOR
972 ! 98-11-02 BLACK - MODIFIED FOR DISTRIBUTED MEMORY
973 ! 99-03-17 TUCCILLO - INCORPORATED MPI_ALLREDUCE FOR GLOBAL SUM
974 ! 02-02-06 BLACK - CONVERTED TO WRF FORMAT
975 ! 02-09-06 WOLFE - MORE CONVERSION TO GLOBAL INDEXING
976 ! 04-11-23 BLACK - THREADED
977 ! 05-12-14 BLACK - CONVERTED FROM IKJ TO IJK
978 ! 07-08-14 janjic - bc & no conservation in the advection step
979 !
980 ! USAGE: CALL VAD2 FROM SUBROUTINE SOLVE_NMM
981 ! INPUT ARGUMENT LIST:
982 !
983 ! OUTPUT ARGUMENT LIST
984 !
985 ! OUTPUT FILES:
986 ! NONE
987 ! SUBPROGRAMS CALLED:
988 !
989 ! UNIQUE: NONE
990 !
991 ! LIBRARY: NONE
992 !
993 ! ATTRIBUTES:
994 ! LANGUAGE: FORTRAN 90
995 ! MACHINE : IBM SP
996 !$$$
997 !***********************************************************************
998 !----------------------------------------------------------------------
999 !
1000 IMPLICIT NONE
1001 !
1002 !----------------------------------------------------------------------
1003 !
1004 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1005 & ,IMS,IME,JMS,JME,KMS,KME &
1006 ,ITS,ITE,JTS,JTE,KTS,KTE
1007 !
1008 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
1009 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: N_IUP_H,N_IUP_V &
1010 & ,N_IUP_ADH,N_IUP_ADV
1011 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V &
1012 & ,IUP_ADH,IUP_ADV
1013 !
1014 INTEGER,INTENT(IN) :: IDTAD,NTSD
1015 !
1016 REAL,INTENT(IN) :: DT,DY,PDTOP
1017 !
1018 REAL,DIMENSION(KMS:KME),INTENT(IN) :: AETA1,AETA2,DETA1,DETA2
1019 !
1020 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DX,HBM2,PDSL
1021 !
1022 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: PETDT
1023 !
1024 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(INOUT) :: CWM,Q,Q2
1025 !
1026 !*** LOCAL VARIABLES
1027 !----------------------------------------------------------------------
1028 !
1029 REAL,PARAMETER :: FF1=0.500
1030 !
1031 LOGICAL,SAVE :: TRADITIONAL=.TRUE.
1032 !
1033 INTEGER :: I,IRECV,J,JFP,JFQ,K,LAP,LLAP
1034 !
1035 INTEGER,DIMENSION(KTS:KTE) :: LA
1036 !
1037 REAL*8 :: ADDT,AFRP,D2PQE,D2PQQ,D2PQW,DEP,DETAP,DQP &
1038 & ,DWP,E00,E4P,EP,EP0,HADDT,HBM2IJ &
1039 & ,Q00,Q4P,QP,QP0 &
1040 & ,rdpdn,rdpup,sfacek,sfacqk,sfacwk,RFC,RR &
1041 & ,SUMNE,SUMNQ,SUMNW,SUMPE,SUMPQ,SUMPW &
1042 & ,W00,W4P,WP,WP0
1043 !
1044 REAL,DIMENSION(KTS:KTE) :: AFR,DEL,DQL,DWL,E3,E4,PETDTK &
1045 & ,RFACE,RFACQ,RFACW,Q3,Q4,W3,W4
1046 !
1047 !-----------------------------------------------------------------------
1048 !***********************************************************************
1049 !-----------------------------------------------------------------------
1050 !
1051 ADDT=REAL(IDTAD)*DT
1052 !
1053 !-----------------------------------------------------------------------
1054 !$omp parallel do &
1055 !$omp& private(afr,afrp,bot,d2pqe,d2pqq,d2pqw,del,dep,detap,dpdn,dpup &
1056 !$omp& ,dql,dqp,dwl,dwp,e00,e3,e4,e4p,ep,ep0,haddt,i,j,k &
1057 !$omp& ,la,lap,llap,petdtk,q00,q3,q4,q4p,qp,qp0,rfacek,rfacqk &
1058 !$omp& ,rfacwk,rfc,rr,sumne,sumnq,sumnw,sumpe,sumpq,sumpw,top &
1059 !$omp& ,w00,w3,w4,w4p,wp,wp0)
1060 !-----------------------------------------------------------------------
1061 !
1062 main_integration: DO J=MYJS2,MYJE2
1063 !
1064 !-----------------------------------------------------------------------
1065 !
1066 main_iloop: DO I=MYIS1_P1,MYIE1_P1
1067 !
1068 !-----------------------------------------------------------------------
1069 !
1070 E3(KTE)=Q2(I,J,KTE)*0.5
1071 !
1072 DO K=KTE-1,KTS,-1
1073 E3(K)=MAX((Q2(I,J,K+1)+Q2(I,J,K))*0.5,EPSQ2)
1074 ENDDO
1075 !
1076 DO K=KTS,KTE
1077 Q3(K)=MAX(Q(I,J,K),EPSQ)
1078 W3(K)=MAX(CWM(I,J,K),CLIMIT)
1079 E4(K)=E3(K)
1080 Q4(K)=Q3(K)
1081 W4(K)=W3(K)
1082 ENDDO
1083 !
1084 IF(TRADITIONAL)THEN
1085 PETDTK(KTE)=PETDT(I,J,KTE-1)*0.5
1086 !
1087 DO K=KTE-1,KTS+1,-1
1088 PETDTK(K)=(PETDT(I,J,K)+PETDT(I,J,K-1))*0.5
1089 ENDDO
1090 !
1091 PETDTK(KTS)=PETDT(I,J,KTS)*0.5
1092 !
1093 ELSE
1094 !
1095 !-----------------------------------------------------------------------
1096 !*** PERFORM HORIZONTAL AVERAGING OF VERTICAL VELOCITY
1097 !-----------------------------------------------------------------------
1098 !
1099 PETDTK(KTE)=(PETDT(I+IHW(J-1),J-1,KTE-1) &
1100 & +PETDT(I+IHE(J-1),J-1,KTE-1) &
1101 & +PETDT(I+IHW(J+1),J+1,KTE-1) &
1102 & +PETDT(I+IHE(J+1),J+1,KTE-1) &
1103 & +PETDT(I,J,KTE-1)*4. )*0.0625
1104 !
1105 DO K=KTE-1,KTS+1,-1
1106 PETDTK(K)=(PETDT(I+IHW(J-1),J-1,K-1) &
1107 +PETDT(I+IHE(J-1),J-1,K-1) &
1108 & +PETDT(I+IHW(J+1),J+1,K-1) &
1109 & +PETDT(I+IHE(J+1),J+1,K-1) &
1110 & +PETDT(I+IHW(J-1),J-1,K ) &
1111 & +PETDT(I+IHE(J-1),J-1,K ) &
1112 & +PETDT(I+IHW(J+1),J+1,K ) &
1113 & +PETDT(I+IHE(J+1),J+1,K ) &
1114 & +(PETDT(I,J,K-1)+PETDT(I,J,K))*4. &
1115 & )*0.0625
1116 ENDDO
1117 !
1118 PETDTK(KTS)=(PETDT(I+IHW(J-1),J-1,KTS) &
1119 & +PETDT(I+IHE(J-1),J-1,KTS) &
1120 & +PETDT(I+IHW(J+1),J+1,KTS) &
1121 & +PETDT(I+IHE(J+1),J+1,KTS) &
1122 & +PETDT(I,J,KTS)*4. )*0.0625
1124 ENDIF
1125 !
1126 !-----------------------------------------------------------------------
1127 !
1128 HADDT=-ADDT*HBM2(I,J)
1129 !
1130 DO K=KTE,KTS,-1
1131 RR=PETDTK(K)*HADDT
1132 !
1133 IF(RR<0.)THEN
1134 LAP=1
1135 ELSE
1136 LAP=-1
1137 ENDIF
1138 !
1139 LA(K)=LAP
1140 LLAP=K+LAP
1141 !
1142 if(llap.gt.kts-1.and.llap.lt.kte+1) then ! internal and outflow pts.
1143 rr=abs(rr &
1144 & /((aeta1(llap)-aeta1(k))*pdtop &
1145 & +(aeta2(llap)-aeta2(k))*pdsl(i,j)))
1146 if(rr.gt.0.999) rr=0.999
1147 !
1148 AFR(K)=(((FF4*RR+FF3)*RR+FF2)*RR+FF1)*RR
1149 dql(k)=(q3(llap)-q3(k))*rr
1150 dwl(k)=(w3(llap)-w3(k))*rr
1151 del(k)=(e3(llap)-e3(k))*rr
1152 elseif(llap.eq.kts-1) then
1153 !
1154 !chem rr=abs(rr &
1155 !chem /((1.-aeta2(kts))*pdsl(i,j)))
1156 !chem afr(kts)=0.
1157 !chem dql(kts)=(epsq -q3(kts))*rr
1158 !chem dwl(kts)=(climit-w3(kts))*rr
1159 !chem del(kts)=(epsq2 -e3(kts))*rr
1160 !
1161 rr=0.
1162 afr(kts)=0.
1163 dql(kts)=0.
1164 dwl(kts)=0.
1165 del(kts)=0.
1166 else
1167 rr=abs(rr &
1168 /(aeta1(kte)*pdtop))
1169 afr(kte)=0.
1170 dql(kte)=(epsq -q3(kte))*rr
1171 dwl(kte)=(climit-w3(kte))*rr
1172 del(kte)=(epsq2 -e3(kte))*rr
1173 endif
1174 ENDDO
1175 !
1176 !-----------------------------------------------------------------------
1177 !
1178 DO K=KTS,KTE
1179 Q4(K)=Q3(K)+DQL(K)
1180 W4(K)=W3(K)+DWL(K)
1181 E4(K)=E3(K)+DEL(K)
1182 ENDDO
1183 !
1184 !-----------------------------------------------------------------------
1185 !*** ANTI-FILTERING STEP
1186 !-----------------------------------------------------------------------
1187 !
1188 SUMPQ=0.
1189 SUMNQ=0.
1190 SUMPW=0.
1191 SUMNW=0.
1192 SUMPE=0.
1193 SUMNE=0.
1194 !
1195 !*** ANTI-FILTERING LIMITERS
1196 !
1197 antifilter: DO K=KTE-1,KTS+1,-1
1198 !
1199 DETAP=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J)
1200 !
1201 DQL(K)=0.
1202 DWL(K)=0.
1203 DEL(K)=0.
1204 !
1205 Q4P=Q4(K)
1206 W4P=W4(K)
1207 E4P=E4(K)
1208 !
1209 LAP=LA(K)
1210 !
1211 if(lap.ne.0)then
1212 rdpdn=1./((aeta1(k+lap)-aeta1(k))*pdtop &
1213 & +(aeta2(k+lap)-aeta2(k))*pdsl(i,j))
1214 rdpup=1./((aeta1(k)-aeta1(k-lap))*pdtop &
1215 & +(aeta2(k)-aeta2(k-lap))*pdsl(i,j))
1216 !
1217 afrp=afr(k)*detap
1218 !
1219 d2pqq=((q4(k+lap)-q4p)*rdpdn &
1220 & -(q4p-q4(k-lap))*rdpup)*afrp
1221 d2pqw=((w4(k+lap)-w4p)*rdpdn &
1222 & -(w4p-w4(k-lap))*rdpup)*afrp
1223 d2pqe=((e4(k+lap)-e4p)*rdpdn &
1224 & -(e4p-e4(k-lap))*rdpup)*afrp
1225 ELSE
1226 D2PQQ=0.
1227 D2PQW=0.
1228 D2PQE=0.
1229 ENDIF
1230 !
1231 QP=Q4P-D2PQQ
1232 WP=W4P-D2PQW
1233 EP=E4P-D2PQE
1234 !
1235 Q00=Q3(K)
1236 QP0=Q3(K+LAP)
1237 !
1238 W00=W3(K)
1239 WP0=W3(K+LAP)
1240 !
1241 E00=E3(K)
1242 EP0=E3(K+LAP)
1243 !
1244 IF(LAP/=0)THEN
1245 QP=MAX(QP,MIN(Q00,QP0))
1246 QP=MIN(QP,MAX(Q00,QP0))
1247 WP=MAX(WP,MIN(W00,WP0))
1248 WP=MIN(WP,MAX(W00,WP0))
1249 EP=MAX(EP,MIN(E00,EP0))
1250 EP=MIN(EP,MAX(E00,EP0))
1251 ENDIF
1252 !
1253 dqp=qp-q4p
1254 dwp=wp-w4p
1255 dep=ep-e4p
1256 !
1257 DQL(K)=DQP
1258 DWL(K)=DWP
1259 DEL(K)=DEP
1260 !
1261 DQP=DQP*DETAP
1262 DWP=DWP*DETAP
1263 DEP=DEP*DETAP
1264 !
1265 IF(DQP>0.)THEN
1266 SUMPQ=SUMPQ+DQP
1267 ELSE
1268 SUMNQ=SUMNQ+DQP
1269 ENDIF
1270 !
1271 IF(DWP>0.)THEN
1272 SUMPW=SUMPW+DWP
1273 ELSE
1274 SUMNW=SUMNW+DWP
1275 ENDIF
1276 !
1277 IF(DEP>0.)THEN
1278 SUMPE=SUMPE+DEP
1279 ELSE
1280 SUMNE=SUMNE+DEP
1281 ENDIF
1282 !
1283 ENDDO antifilter
1284 !
1285 !-----------------------------------------------------------------------
1286 !
1287 DQL(KTS)=0.
1288 DWL(KTS)=0.
1289 DEL(KTS)=0.
1290 !
1291 DQL(KTE)=0.
1292 DWL(KTE)=0.
1293 DEL(KTE)=0.
1294 !
1295 !-----------------------------------------------------------------------
1296 !*** FIRST MOMENT CONSERVING FACTOR
1297 !-----------------------------------------------------------------------
1298 !
1299 if(sumpq*(-sumnq).gt.1.e-9) then
1300 sfacqk=-sumnq/sumpq
1301 else
1302 sfacqk=0.
1303 endif
1304 !
1305 if(sumpw*(-sumnw).gt.1.e-9) then
1306 sfacwk=-sumnw/sumpw
1307 else
1308 sfacwk=0.
1309 endif
1310 !
1311 if(sumpe*(-sumne).gt.1.e-9) then
1312 sfacek=-sumne/sumpe
1313 else
1314 sfacek=0.
1315 endif
1316 !
1317 !-----------------------------------------------------------------------
1318 !*** IMPOSE CONSERVATION ON ANTI-FILTERING
1319 !-----------------------------------------------------------------------
1320 !
1321 DO K=KTE,KTS,-1
1322 !
1323 dqp=dql(k)
1324 if(sfacqk.gt.0.) then
1325 if(sfacqk.ge.1.) then
1326 if(dqp.lt.0.) dqp=dqp/sfacqk
1327 else
1328 if(dqp.gt.0.) dqp=dqp*sfacqk
1329 endif
1330 else
1331 dqp=0.
1332 endif
1333 q (i,j,k)=q4(k)+dqp
1334 !
1335 dwp=dwl(k)
1336 if(sfacwk.gt.0.) then
1337 if(sfacwk.ge.1.) then
1338 if(dwp.lt.0.) dwp=dwp/sfacwk
1339 else
1340 if(dwp.gt.0.) dwp=dwp*sfacwk
1341 endif
1342 else
1343 dwp=0.
1344 endif
1345 cwm(i,j,k)=w4(k)+dwp
1346 !
1347 dep=del(k)
1348 if(sfacek.gt.0.) then
1349 if(sfacek.ge.1.) then
1350 if(dep.lt.0.) dep=dep/sfacek
1351 else
1352 if(dep.gt.0.) dep=dep*sfacek
1353 endif
1354 else
1355 dep=0.
1356 endif
1357 e3 ( k)=e4(k)+dep
1358 !
1359 ENDDO
1360 !-----------------------------------------------------------------------
1361 HBM2IJ=HBM2(I,J)
1362 Q2(I,J,KTE)=MAX(E3(KTE)+E3(KTE)-EPSQ2,EPSQ2)*HBM2IJ &
1363 & +Q2(I,J,KTE)*(1.-HBM2IJ)
1364 DO K=KTE-1,KTS+1,-1
1365 Q2(I,J,K)=MAX(E3(K)+E3(K)-Q2(I,J,K+1),EPSQ2)*HBM2IJ &
1366 & +Q2(I,J,K)*(1.-HBM2IJ)
1367 ENDDO
1368 !-----------------------------------------------------------------------
1369 !
1370 ENDDO main_iloop
1371 !
1372 !-----------------------------------------------------------------------
1373 !
1374 ENDDO main_integration
1375 !
1376 !-----------------------------------------------------------------------
1377 !
1378 END SUBROUTINE VAD2
1379 !
1381 !-----------------------------------------------------------------------
1382 !
1383 !-----------------------------------------------------------------------
1384 !***********************************************************************
1385 SUBROUTINE HAD2( &
1386 #if defined(DM_PARALLEL)
1387 & domdesc , &
1388 #endif
1389 & NTSD,DT,IDTAD,DX,DY &
1390 & ,AETA1,AETA2,DETA1,DETA2,PDSL,PDTOP &
1391 & ,HBM2,HBM3 &
1392 & ,Q,Q2,CWM,U,V,Z,HYDRO &
1393 & ,N_IUP_H,N_IUP_V &
1394 & ,N_IUP_ADH,N_IUP_ADV &
1395 & ,IUP_H,IUP_V,IUP_ADH,IUP_ADV &
1396 & ,IHE,IHW,IVE,IVW &
1397 & ,IDS,IDE,JDS,JDE,KDS,KDE &
1398 & ,IMS,IME,JMS,JME,KMS,KME &
1399 & ,ITS,ITE,JTS,JTE,KTS,KTE)
1400 !***********************************************************************
1401 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
1402 ! . . .
1403 ! SUBPROGRAM: HAD2 HORIZONTAL ADVECTION OF H2O AND TKE
1404 ! PRGRMMR: JANJIC ORG: W/NP22 DATE: 96-07-19
1405 !
1406 ! ABSTRACT:
1407 ! HAD2 CALCULATES THE CONTRIBUTION OF THE HORIZONTAL ADVECTION
1408 ! TO THE TENDENCIES OF WATER SUBSTANCE AND TKE AND THEN
1409 ! UPDATES THOSE VARIABLES. AN ANTI-FILTERING TECHNIQUE IS USED.
1410 !
1411 ! PROGRAM HISTORY LOG:
1412 ! 96-07-19 JANJIC - ORIGINATOR
1413 ! 98-11-02 BLACK - MODIFIED FOR DISTRIBUTED MEMORY
1414 ! 99-03-17 TUCCILLO - INCORPORATED MPI_ALLREDUCE FOR GLOBAL SUM
1415 ! 02-02-06 BLACK - CONVERTED TO WRF FORMAT
1416 ! 02-09-06 WOLFE - MORE CONVERSION TO GLOBAL INDEXING
1417 ! 03-05-23 JANJIC - ADDED SLOPE FACTOR
1418 ! 04-11-23 BLACK - THREADED
1419 ! 05-12-14 BLACK - CONVERTED FROM IKJ TO IJK
1420 ! 07-08-14 janjic - no conservation in advection step
1421 !
1422 ! USAGE: CALL HAD2 FROM SUBROUTINE SOLVE_NMM
1423 ! INPUT ARGUMENT LIST:
1424 !
1425 ! OUTPUT ARGUMENT LIST
1426 !
1427 ! OUTPUT FILES:
1428 ! NONE
1429 ! SUBPROGRAMS CALLED:
1430 !
1431 ! UNIQUE: NONE
1432 !
1433 ! LIBRARY: NONE
1434 !
1435 ! ATTRIBUTES:
1436 ! LANGUAGE: FORTRAN 90
1437 ! MACHINE : IBM SP
1438 !$$$
1439 !***********************************************************************
1440 !-----------------------------------------------------------------------
1441 !
1442 IMPLICIT NONE
1443 !
1444 !-----------------------------------------------------------------------
1445 !
1446 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
1447 & ,IMS,IME,JMS,JME,KMS,KME &
1448 & ,ITS,ITE,JTS,JTE,KTS,KTE
1449 !
1450 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
1451 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: N_IUP_H,N_IUP_V &
1452 & ,N_IUP_ADH,N_IUP_ADV
1453 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V &
1454 & ,IUP_ADH,IUP_ADV
1455 !
1456 !-----------------------------------------------------------------------
1457 !
1458 INTEGER,INTENT(IN) :: IDTAD,NTSD
1459 !
1460 REAL,INTENT(IN) :: DT,DY,PDTOP
1461 !
1462 REAL,DIMENSION(KMS:KME),INTENT(IN) :: AETA1,AETA2,DETA1,DETA2
1463 !
1464 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DX,HBM2,HBM3,PDSL
1465 !
1466 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: U,V,Z
1467 !
1468 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(INOUT) :: CWM,Q,Q2
1469 !
1470 LOGICAL,INTENT(IN) :: HYDRO
1471 !
1472 !-----------------------------------------------------------------------
1473 !*** LOCAL VARIABLES
1474 !-----------------------------------------------------------------------
1475 !
1476 REAL,PARAMETER :: FF1=0.530
1477 !
1478 #ifdef DM_PARALLEL
1479 INTEGER :: DOMDESC
1480 #endif
1481 !
1482 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
1483 LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
1484 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME,6) :: XSUMS_L
1485 REAL,DIMENSION(IDS:IDE,JDS:JDE,KDS:KDE,6) :: XSUMS_G
1486 #endif
1487 !
1488 LOGICAL :: BOT,TOP
1489 !
1490 INTEGER :: I,IRECV,J,JFP,JFQ,K,LAP,LLAP,MPI_COMM_COMP
1491 INTEGER :: N
1492 !
1493 INTEGER,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5,KTS:KTE) :: IFPA,IFPF &
1494 & ,IFQA,IFQF &
1495 & ,JFPA,JFPF &
1496 & ,JFQA,JFQF
1497 !
1498 REAL :: ADDT,AFRP,CRIT,D2PQE,D2PQQ,D2PQW,DEP,DESTIJ,DQP,DQSTIJ &
1499 & ,DVOLP,DWP,DWSTIJ,DZA,DZB,E00,E0Q,E1X,E2IJ,E4P,ENH,EP,EP0 &
1500 & ,ESTIJ,FPQ,HAFP,HAFQ,HBM2IJ,HM,PP,PPQ00,Q00,Q0Q &
1501 & ,Q1IJ,Q4P,QP,QP0,QSTIJ,RDY,RFACE,RFACQ,RFACW,RFC &
1502 & ,RFEIJ,RFQIJ,RFWIJ,RR,SLOPAC,SPP,SQP,SSA,SSB,SUMNE,SUMNQ &
1503 & ,SUMNW,SUMPE,SUMPQ,SUMPW,TTA,TTB,W00,W0Q,W1IJ,W4P,WP,WP0 &
1504 & ,WSTIJ
1505 !
1506 DOUBLE PRECISION,DIMENSION(6,KTS:KTE) :: GSUMS,XSUMS
1507 !
1508 REAL,DIMENSION(KTS:KTE) :: AFR,DEL,DQL,DWL,E3,E4 &
1509 & ,Q3,Q4,W3,W4
1510 !
1511 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: DARE,EMH
1512 !
1513 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5,KTS:KTE) :: AFP,AFQ,DEST &
1514 & ,DQST,DVOL,DWST &
1515 & ,E1,E2,Q1,W1
1516 !-----------------------------------------------------------------------
1517 integer :: nunit,ier
1518 save nunit
1519 !-----------------------------------------------------------------------
1520 !***********************************************************************
1521 !-----------------------------------------------------------------------
1522 !
1523 RDY=1./DY
1524 SLOPAC=SLOPHT*SQRT(2.)*0.5*50.
1525 CRIT=SLOPAC*REAL(IDTAD)*DT*RDY*1000.
1526 !
1527 ADDT=REAL(IDTAD)*DT
1528 ENH=ADDT/(08.*DY)
1529 !
1530 !-----------------------------------------------------------------------
1531 !$omp parallel do &
1532 !$omp& private(i,j)
1533 DO J=MYJS_P3,MYJE_P3
1534 DO I=MYIS_P2,MYIE_P2
1535 EMH (I,J)=ADDT/(08.*DX(I,J))
1536 DARE(I,J)=HBM3(I,J)*DX(I,J)*DY
1537 E1(I,J,KTE)=MAX(Q2(I,J,KTE)*0.5,EPSQ2)
1538 E2(I,J,KTE)=E1(I,J,KTE)
1539 ENDDO
1540 ENDDO
1541 !-----------------------------------------------------------------------
1542 !$omp parallel do &
1543 !$omp& private(dza,dzb,e1x,fpq,hm,i,j,jfp,jfq,k,pp,qp,ssa,ssb,spp,sqp &
1544 !$omp& ,tta,ttb)
1545 !-----------------------------------------------------------------------
1546 !
1547 vertical_1: DO K=KTS,KTE
1548 !
1549 !-----------------------------------------------------------------------
1550 !
1551 DO J=MYJS_P3,MYJE_P3
1552 DO I=MYIS_P2,MYIE_P2
1553 DVOL(I,J,K)=DARE(I,J)*(DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J))
1554 Q (I,J,K)=MAX(Q (I,J,K),EPSQ)
1555 CWM(I,J,K)=MAX(CWM(I,J,K),CLIMIT)
1556 Q1 (I,J,K)=Q (I,J,K)
1557 W1 (I,J,K)=CWM(I,J,K)
1558 ENDDO
1559 ENDDO
1560 !
1561 IF(K<KTE)THEN
1562 DO J=MYJS_P3,MYJE_P3
1563 DO I=MYIS_P2,MYIE_P2
1564 E1X=(Q2(I,J,K+1)+Q2(I,J,K))*0.5
1565 E1(I,J,K)=MAX(E1X,EPSQ2)
1566 E2(I,J,K)=E1(I,J,K)
1567 ENDDO
1568 ENDDO
1569 ENDIF
1570 !
1571 !-----------------------------------------------------------------------
1572 !
1573 DO J=MYJS2_P1,MYJE2_P1
1574 DO I=MYIS1_P1,MYIE1_P1
1575 !
1576 HM=HBM2(I,J)
1577 TTA=(U(I,J-1,K)+U(I+IHW(J),J,K)+U(I+IHE(J),J,K)+U(I,J+1,K)) &
1578 & *EMH(I,J)*HM
1579 TTB=(V(I,J-1,K)+V(I+IHW(J),J,K)+V(I+IHE(J),J,K)+V(I,J+1,K)) &
1580 & *ENH*HBM2(I,J)
1581 !
1582 SPP=-TTA-TTB
1583 SQP= TTA-TTB
1584 !
1585 IF(SPP<0.)THEN
1586 JFP=-1
1587 ELSE
1588 JFP=1
1589 ENDIF
1590 IF(SQP<0.)THEN
1591 JFQ=-1
1592 ELSE
1593 JFQ=1
1594 ENDIF
1595 !
1596 IFPA(I,J,K)=IHE(J)+I+( JFP-1)/2
1597 IFQA(I,J,K)=IHE(J)+I+(-JFQ-1)/2
1598 !
1599 JFPA(I,J,K)=J+JFP
1600 JFQA(I,J,K)=J+JFQ
1601 !
1602 IFPF(I,J,K)=IHE(J)+I+(-JFP-1)/2
1603 IFQF(I,J,K)=IHE(J)+I+( JFQ-1)/2
1604 !
1605 JFPF(I,J,K)=J-JFP
1606 JFQF(I,J,K)=J-JFQ
1607 ! if(i==111.and.j==438.and.k==1)then
1608 ! endif
1609 !
1610 !-----------------------------------------------------------------------
1611 IF(.NOT.HYDRO)THEN ! z currently not available for hydro=.true.
1612 DZA=(Z(IFPA(I,J,K),JFPA(I,J,K),K)-Z(I,J,K))*RDY
1613 DZB=(Z(IFQA(I,J,K),JFQA(I,J,K),K)-Z(I,J,K))*RDY
1614 !
1615 IF(ABS(DZA)>SLOPAC)THEN
1616 SSA=DZA*SPP
1617 IF(SSA>CRIT)THEN
1618 SPP=0. !spp*.1
1619 ENDIF
1620 ENDIF
1621 !
1622 IF(ABS(DZB)>SLOPAC)THEN
1623 SSB=DZB*SQP
1624 IF(SSB>CRIT)THEN
1625 SQP=0. !sqp*.1
1626 ENDIF
1627 ENDIF
1628 !
1629 ENDIF
1630 !
1631 !-----------------------------------------------------------------------
1632 !
1633 FPQ=SPP*SQP*0.25
1634 PP=ABS(SPP)
1635 QP=ABS(SQP)
1636 !
1637 AFP(I,J,K)=(((FF4*PP+FF3)*PP+FF2)*PP+FF1)*PP
1638 AFQ(I,J,K)=(((FF4*QP+FF3)*QP+FF2)*QP+FF1)*QP
1639 !
1640 Q1(I,J,K)=(Q (IFPA(I,J,K),JFPA(I,J,K),K)-Q (I,J,K))*PP &
1641 & +(Q (IFQA(I,J,K),JFQA(I,J,K),K)-Q (I,J,K))*QP &
1642 & +(Q (I,J-2,K)+Q (I,J+2,K) &
1643 & -Q (I-1,J,K)-Q (I+1,J,K))*FPQ &
1644 & +Q(I,J,K)
1645 !
1646 W1(I,J,K)=(CWM(IFPA(I,J,K),JFPA(I,J,K),K)-CWM(I,J,K))*PP &
1647 & +(CWM(IFQA(I,J,K),JFQA(I,J,K),K)-CWM(I,J,K))*QP &
1648 & +(CWM(I,J-2,K)+CWM(I,J+2,K) &
1649 & -CWM(I-1,J,K)-CWM(I+1,J,K))*FPQ &
1650 & +CWM(I,J,K)
1651 !
1652 E2(I,J,K)=(E1 (IFPA(I,J,K),JFPA(I,J,K),K)-E1 (I,J,K))*PP &
1653 & +(E1 (IFQA(I,J,K),JFQA(I,J,K),K)-E1 (I,J,K))*QP &
1654 & +(E1 (I,J-2,K)+E1 (I,J+2,K) &
1655 & -E1 (I-1,J,K)-E1 (I+1,J,K))*FPQ &
1656 & +E1(I,J,K)
1657 !
1658 ENDDO
1659 ENDDO
1660 !
1661 !-----------------------------------------------------------------------
1662 !
1663 ENDDO vertical_1
1664 !
1665 !-----------------------------------------------------------------------
1666 !*** ANTI-FILTERING STEP
1667 !-----------------------------------------------------------------------
1668 !
1669 DO K=KTS,KTE
1670 XSUMS(1,K)=0.
1671 XSUMS(2,K)=0.
1672 XSUMS(3,K)=0.
1673 XSUMS(4,K)=0.
1674 XSUMS(5,K)=0.
1675 XSUMS(6,K)=0.
1676 ENDDO
1677 !-----------------------------------------------------------------------
1678 !
1679 !*** ANTI-FILTERING LIMITERS
1680 !
1681 !-----------------------------------------------------------------------
1682 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
1683 DO N=1,6
1684 !
1685 !$omp parallel do &
1686 !$omp& private(i,j,k)
1687 DO K=KMS,KME
1688 DO J=JMS,JME
1689 DO I=IMS,IME
1690 XSUMS_L(I,J,K,N)=0.
1691 ENDDO
1692 ENDDO
1693 ENDDO
1694 !
1695 !$omp parallel do &
1696 !$omp& private(i,j,k)
1697 DO K=KDS,KDE
1698 DO J=JDS,JDE
1699 DO I=IDS,IDE
1700 XSUMS_G(I,J,K,N)=0.
1701 ENDDO
1702 ENDDO
1703 ENDDO
1704 !
1705 ENDDO
1706 !
1707 #endif
1708 !-----------------------------------------------------------------------
1709 !$omp parallel do &
1710 !$omp& private(d2pqe,d2pqq,d2pqw,destij,dqstij,dvolp,dwstij &
1711 !$omp& ,e00,e0q,e2ij,ep0,estij,hafp,hafq,i,j,k &
1712 !$omp& ,q00,q0q,q1ij,qp0,qstij,w00,w0q,w1ij,wp0,wstij)
1713 !-----------------------------------------------------------------------
1714 !
1715 vertical_2: DO K=KTS,KTE
1716 !
1717 !-----------------------------------------------------------------------
1718 !
1719 DO J=MYJS2,MYJE2
1720 DO I=MYIS1,MYIE1
1721 !
1722 DVOLP=DVOL(I,J,K)
1723 Q1IJ =Q1(I,J,K)
1724 W1IJ =W1(I,J,K)
1725 E2IJ =E2(I,J,K)
1726 !
1727 HAFP=AFP(I,J,K)
1728 HAFQ=AFQ(I,J,K)
1729 !
1730 D2PQQ=(Q1(IFPA(I,J,K),JFPA(I,J,K),K)-Q1IJ &
1731 & -Q1IJ+Q1(IFPF(I,J,K),JFPF(I,J,K),K)) &
1732 & *HAFP &
1733 & +(Q1(IFQA(I,J,K),JFQA(I,J,K),K)-Q1IJ &
1734 & -Q1IJ+Q1(IFQF(I,J,K),JFQF(I,J,K),K)) &
1735 & *HAFQ
1736 !
1737 D2PQW=(W1(IFPA(I,J,K),JFPA(I,J,K),K)-W1IJ &
1738 & -W1IJ+W1(IFPF(I,J,K),JFPF(I,J,K),K)) &
1739 & *HAFP &
1740 & +(W1(IFQA(I,J,K),JFQA(I,J,K),K)-W1IJ &
1741 & -W1IJ+W1(IFQF(I,J,K),JFQF(I,J,K),K)) &
1742 & *HAFQ
1743 !
1744 D2PQE=(E2(IFPA(I,J,K),JFPA(I,J,K),K)-E2IJ &
1745 & -E2IJ+E2(IFPF(I,J,K),JFPF(I,J,K),K)) &
1746 & *HAFP &
1747 & +(E2(IFQA(I,J,K),JFQA(I,J,K),K)-E2IJ &
1748 & -E2IJ+E2(IFQF(I,J,K),JFQF(I,J,K),K)) &
1749 & *HAFQ
1750 !
1751 QSTIJ=Q1IJ-D2PQQ
1752 WSTIJ=W1IJ-D2PQW
1753 ESTIJ=E2IJ-D2PQE
1754 !
1755 Q00=Q (I ,J ,K)
1756 QP0=Q (IFPA(I,J,K),JFPA(I,J,K),K)
1757 Q0Q=Q (IFQA(I,J,K),JFQA(I,J,K),K)
1758 !
1759 W00=CWM(I ,J ,K)
1760 WP0=CWM(IFPA(I,J,K),JFPA(I,J,K),K)
1761 W0Q=CWM(IFQA(I,J,K),JFQA(I,J,K),K)
1762 !
1763 E00=E1 (I ,J ,K)
1764 EP0=E1 (IFPA(I,J,K),JFPA(I,J,K),K)
1765 E0Q=E1 (IFQA(I,J,K),JFQA(I,J,K),K)
1766 !
1767 QSTIJ=MAX(QSTIJ,MIN(Q00,QP0,Q0Q))
1768 QSTIJ=MIN(QSTIJ,MAX(Q00,QP0,Q0Q))
1769 WSTIJ=MAX(WSTIJ,MIN(W00,WP0,W0Q))
1770 WSTIJ=MIN(WSTIJ,MAX(W00,WP0,W0Q))
1771 ESTIJ=MAX(ESTIJ,MIN(E00,EP0,E0Q))
1772 ESTIJ=MIN(ESTIJ,MAX(E00,EP0,E0Q))
1773 !
1774 ! DQSTIJ=QSTIJ-Q(I,J,K)
1775 ! DWSTIJ=WSTIJ-CWM(I,J,K)
1776 ! DESTIJ=ESTIJ-E1(I,J,K)
1777 !
1778 dqstij=qstij-q1(i,j,k)
1779 dwstij=wstij-w1(i,j,k)
1780 destij=estij-e2(i,j,k)
1781 !
1782 DQST(I,J,K)=DQSTIJ
1783 DWST(I,J,K)=DWSTIJ
1784 DEST(I,J,K)=DESTIJ
1785 !
1786 DQSTIJ=DQSTIJ*DVOLP
1787 DWSTIJ=DWSTIJ*DVOLP
1788 DESTIJ=DESTIJ*DVOLP
1789 !
1790 !-----------------------------------------------------------------------
1791 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
1792 !-----------------------------------------------------------------------
1793 DO N=1,6
1794 XSUMS_L(I,J,K,N)=0.
1795 ENDDO
1796 !
1797 IF(DQSTIJ>0.)THEN
1798 XSUMS_L(I,J,K,1)=DQSTIJ
1799 ELSE
1800 XSUMS_L(I,J,K,2)=DQSTIJ
1801 ENDIF
1802 !
1803 IF(DWSTIJ>0.)THEN
1804 XSUMS_L(I,J,K,3)=DWSTIJ
1805 ELSE
1806 XSUMS_L(I,J,K,4)=DWSTIJ
1807 ENDIF
1808 !
1809 IF(DESTIJ>0.)THEN
1810 XSUMS_L(I,J,K,5)=DESTIJ
1811 ELSE
1812 XSUMS_L(I,J,K,6)=DESTIJ
1813 ENDIF
1814 !-----------------------------------------------------------------------
1815 #else
1816 !-----------------------------------------------------------------------
1817 IF(DQSTIJ>0.)THEN
1818 XSUMS(1,K)=XSUMS(1,K)+DQSTIJ
1819 ELSE
1820 XSUMS(2,K)=XSUMS(2,K)+DQSTIJ
1821 ENDIF
1822 !
1823 IF(DWSTIJ>0.)THEN
1824 XSUMS(3,K)=XSUMS(3,K)+DWSTIJ
1825 ELSE
1826 XSUMS(4,K)=XSUMS(4,K)+DWSTIJ
1827 ENDIF
1828 !
1829 IF(DESTIJ>0.)THEN
1830 XSUMS(5,K)=XSUMS(5,K)+DESTIJ
1831 ELSE
1832 XSUMS(6,K)=XSUMS(6,K)+DESTIJ
1833 ENDIF
1834 !-----------------------------------------------------------------------
1835 #endif
1836 !-----------------------------------------------------------------------
1837 !
1838 ENDDO
1839 ENDDO
1840 !
1841 !-----------------------------------------------------------------------
1842 !
1843 ENDDO vertical_2
1844 !
1845 !-----------------------------------------------------------------------
1846 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
1847 !-----------------------------------------------------------------------
1848 DO N=1,6
1849 CALL WRF_PATCH_TO_GLOBAL_REAL(XSUMS_L(IMS,JMS,KMS,N) &
1850 &, XSUMS_G(1,1,1,N),DOMDESC &
1851 &, 'xyz','xyz' &
1852 &, IDS,IDE,JDS,JDE,KDS,KDE &
1853 &, IMS,IME,JMS,JME,KMS,KME &
1854 &, ITS,ITE,JTS,JTE,KTS,KTE)
1855 ENDDO
1856 !
1857 DO K=KTS,KTE
1858 DO N=1,6
1859 GSUMS(N,K)=0.
1860 ENDDO
1861 ENDDO
1862 !
1863 IF(WRF_DM_ON_MONITOR())THEN
1864 DO N=1,6
1865 !$omp parallel do &
1866 !$omp& private(i,j,k)
1867 DO K=KTS,KTE
1868 DO J=JDS,JDE
1869 DO I=IDS,IDE
1870 GSUMS(N,K)=GSUMS(N,K)+XSUMS_G(I,J,K,N)
1871 ENDDO
1872 ENDDO
1873 ENDDO
1874 ENDDO
1875 ENDIF
1877 CALL WRF_DM_BCAST_BYTES(GSUMS,2*RWORDSIZE*6*(KTE-KTS+1) )
1879 !-----------------------------------------------------------------------
1880 #else
1881 !-----------------------------------------------------------------------
1882 !
1883 !-----------------------------------------------------------------------
1884 !*** GLOBAL REDUCTION
1885 !-----------------------------------------------------------------------
1886 !
1887 # if defined(DM_PARALLEL) && !defined(STUBMPI)
1888 CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
1889 CALL MPI_ALLREDUCE(XSUMS,GSUMS,6*(KTE-KTS+1) &
1890 & ,MPI_DOUBLE_PRECISION,MPI_SUM &
1891 & ,MPI_COMM_COMP,IRECV)
1892 # else
1893 DO K=KTS,KTE
1894 DO N=1,6
1895 GSUMS(N,K)=XSUMS(N,K)
1896 ENDDO
1897 ENDDO
1898 # endif
1899 !
1900 !-----------------------------------------------------------------------
1901 #endif
1902 !-----------------------------------------------------------------------
1903 !
1904 !-----------------------------------------------------------------------
1905 !*** END OF GLOBAL REDUCTION
1906 !-----------------------------------------------------------------------
1907 !
1908 ! if(mype==0)then
1909 !!! if(ntsd==0)then
1910 !!! call int_get_fresh_handle(nunit)
1911 !!! close(nunit)
1912 ! nunit=56
1913 !!! open(unit=nunit,file='gsums',form='unformatted',iostat=ier)
1914 !!! endif
1915 ! endif
1916 !-----------------------------------------------------------------------
1917 !$omp parallel do &
1918 !$omp& private(destij,dqstij,dwstij,i,j,k,rface,rfacq,rfacw &
1919 !$omp& ,rfeij,rfqij,rfwij,sumne,sumnq,sumnw,sumpe,sumpq,sumpw)
1920 !-----------------------------------------------------------------------
1921 !
1922 vertical_3: DO K=KTS,KTE
1923 !
1924 !-----------------------------------------------------------------------
1925 ! if(mype==0)then
1926 ! write(nunit)(gsums(i,k),i=1,6)
1927 ! endif
1928 !!! read(nunit)(gsums(i,k),i=1,6)
1929 !-----------------------------------------------------------------------
1930 !
1931 SUMPQ=GSUMS(1,K)
1932 SUMNQ=GSUMS(2,K)
1933 SUMPW=GSUMS(3,K)
1934 SUMNW=GSUMS(4,K)
1935 SUMPE=GSUMS(5,K)
1936 SUMNE=GSUMS(6,K)
1937 !
1938 !-----------------------------------------------------------------------
1939 !*** FIRST MOMENT CONSERVING FACTOR
1940 !-----------------------------------------------------------------------
1941 !
1942 IF(SUMPQ>1.)THEN
1943 RFACQ=-SUMNQ/SUMPQ
1944 ELSE
1945 RFACQ=1.
1946 ENDIF
1947 !
1948 IF(SUMPW>1.)THEN
1949 RFACW=-SUMNW/SUMPW
1950 ELSE
1951 RFACW=1.
1952 ENDIF
1953 !
1954 IF(SUMPE>1.)THEN
1955 RFACE=-SUMNE/SUMPE
1956 ELSE
1957 RFACE=1.
1958 ENDIF
1959 !
1960 IF(RFACQ<CONSERVE_MIN.OR.RFACQ>CONSERVE_MAX)RFACQ=1.
1961 IF(RFACW<CONSERVE_MIN.OR.RFACW>CONSERVE_MAX)RFACW=1.
1962 IF(RFACE<CONSERVE_MIN.OR.RFACE>CONSERVE_MAX)RFACE=1.
1963 !
1964 !-----------------------------------------------------------------------
1965 ! if(mype==0.and.ntsd==181)close(nunit)
1966 !-----------------------------------------------------------------------
1967 !
1968 !-----------------------------------------------------------------------
1969 !*** IMPOSE CONSERVATION ON ANTI-FILTERING
1970 !-----------------------------------------------------------------------
1971 !
1972 if(rfacq<1.)then
1973 do j=MYJS2,MYJE2
1974 DO I=MYIS1,MYIE1
1975 DQSTIJ=DQST(I,J,K)
1976 RFQIJ=HBM2(I,J)*(RFACQ-1.)+1.
1977 IF(DQSTIJ>=0.)DQSTIJ=DQSTIJ*RFQIJ
1978 q(i,j,k)=q1(i,j,k)+dqstij
1979 ENDDO
1980 enddo
1981 else
1982 do j=MYJS2,MYJE2
1983 DO I=MYIS1,MYIE1
1984 DQSTIJ=DQST(I,J,K)
1985 RFQIJ=HBM2(I,J)*(RFACQ-1.)+1.
1986 IF(DQSTIJ<0.)DQSTIJ=DQSTIJ/RFQIJ
1987 q(i,j,k)=q1(i,j,k)+dqstij
1988 ENDDO
1989 enddo
1990 endif
1991 !
1992 !-----------------------------------------------------------------------
1993 !
1994 if(rfacw<1.)then
1995 do j=MYJS2,MYJE2
1996 DO I=MYIS1,MYIE1
1997 DWSTIJ=DWST(I,J,K)
1998 RFWIJ=HBM2(I,J)*(RFACW-1.)+1.
1999 IF(DWSTIJ>=0.)DWSTIJ=DWSTIJ*RFWIJ
2000 cwm(i,j,k)=w1(i,j,k)+dwstij
2001 ENDDO
2002 enddo
2003 else
2004 do j=MYJS2,MYJE2
2005 DO I=MYIS1,MYIE1
2006 DWSTIJ=DWST(I,J,K)
2007 RFWIJ=HBM2(I,J)*(RFACW-1.)+1.
2008 IF(DWSTIJ<0.)DWSTIJ=DWSTIJ/RFWIJ
2009 cwm(i,j,k)=w1(i,j,k)+dwstij
2010 ENDDO
2011 enddo
2012 endif
2014 !-----------------------------------------------------------------------
2015 !
2016 if(rface<1.)then
2017 do j=MYJS2,MYJE2
2018 DO I=MYIS1,MYIE1
2019 DESTIJ=DEST(I,J,K)
2020 RFEIJ=HBM2(I,J)*(RFACE-1.)+1.
2021 IF(DESTIJ>=0.)DESTIJ=DESTIJ*RFEIJ
2022 e1(i,j,k)=e2(i,j,k)+destij
2023 ENDDO
2024 enddo
2025 else
2026 do j=MYJS2,MYJE2
2027 DO I=MYIS1,MYIE1
2028 DESTIJ=DEST(I,J,K)
2029 RFEIJ=HBM2(I,J)*(RFACE-1.)+1.
2030 IF(DESTIJ<0.)DESTIJ=DESTIJ/RFEIJ
2031 e1(i,j,k)=e2(i,j,k)+destij
2032 ENDDO
2033 enddo
2034 endif
2035 !
2036 !-----------------------------------------------------------------------
2037 !
2038 DO J=MYJS,MYJE
2039 DO I=MYIS,MYIE
2040 Q (I,J,K)=MAX(Q (I,J,K),EPSQ)
2041 CWM(I,J,K)=MAX(CWM(I,J,K),CLIMIT)
2042 ENDDO
2043 ENDDO
2044 !
2045 !-----------------------------------------------------------------------
2046 !
2047 ENDDO vertical_3
2048 !
2049 !-----------------------------------------------------------------------
2050 !
2051 !$omp parallel do &
2052 !$omp& private(i,j)
2053 DO J=MYJS,MYJE
2054 DO I=MYIS,MYIE
2055 Q2(I,J,KTE)=MAX(E1(I,J,KTE)+E1(I,J,KTE)-EPSQ2,EPSQ2)
2056 ENDDO
2057 ENDDO
2058 !
2059 !-----------------------------------------------------------------------
2060 !
2061 DO K=KTE-1,KTS+1,-1
2062 !$omp parallel do &
2063 !$omp& private(i,j)
2064 DO J=MYJS,MYJE
2065 DO I=MYIS,MYIE
2066 IF(K>KTS)THEN
2067 Q2(I,J,K)=MAX(E1(I,J,K)+E1(I,J,K)-Q2(I,J,K+1),EPSQ2)
2068 ELSE
2069 Q2(I,J,K)=Q2(I,J,K+1)
2070 ENDIF
2071 ENDDO
2072 ENDDO
2073 ENDDO
2074 !-----------------------------------------------------------------------
2075 !
2076 END SUBROUTINE HAD2
2077 !
2078 !-----------------------------------------------------------------------
2079 !***********************************************************************
2080 !-----------------------------------------------------------------------
2081 ! New routines added by Georg Grell to handle advection more like ARW
2082 ! core. Instead of VAD2/HAD2 that advect TKE, specific humidity, and
2083 ! condensed water species all in one routine, we call VAD2/HAD2_SCAL
2084 ! with multidimensioned arrays to advect each variable. For purposes
2085 ! here, solve_nmm.F calls this routine once for TKE, then again for
2086 ! all the species held in the moist array (qv, qc, qi, qr, qs, qg),
2087 ! then call again for number concentrations held in scalar array (qni).
2088 ! The dummy argument lstart is the starting index of the multidimensioned
2089 ! array for starting the advection since the 1st index of moist and
2090 ! scalar are actually empty placeholders (and the 2nd element is vapor,
2091 ! then qc, etc.) When calling with single 3D array (like TKE), just
2092 ! set NUM_SCAL=1 and lstart=1. The variable to advect is called SCAL
2093 ! herein.
2094 !***********************************************************************
2095 SUBROUTINE VAD2_SCAL(NTSD,DT,IDTAD,DX,DY &
2096 & ,AETA1,AETA2,DETA1,DETA2,PDSL,PDTOP &
2097 & ,HBM2 &
2098 & ,SCAL,PETDT &
2099 & ,N_IUP_H,N_IUP_V &
2100 & ,N_IUP_ADH,N_IUP_ADV &
2101 & ,IUP_H,IUP_V,IUP_ADH,IUP_ADV &
2102 & ,IHE,IHW,IVE,IVW &
2103 & ,NUM_SCAL,LSTART &
2104 & ,IDS,IDE,JDS,JDE,KDS,KDE &
2105 & ,IMS,IME,JMS,JME,KMS,KME &
2106 & ,ITS,ITE,JTS,JTE,KTS,KTE)
2107 !***********************************************************************
2108 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2109 ! . . .
2110 ! SUBPROGRAM: VAD2_SCAL VERTICAL ADVECTION OF SCALARS
2111 !
2112 ! PRGRMMR: JANJIC ORG: W/NP22 DATE: 96-07-19
2113 ! GRELL,PECKHAM ORG: NOAA/FSL DATE: 05-02-03
2114 !
2115 ! ABSTRACT:
2116 ! VAD2_SCAL CALCULATES THE CONTRIBUTION OF THE VERTICAL ADVECTION
2117 ! TO THE TENDENCIES OF SCALAR SUBSTANCES AND THEN UPDATES
2118 ! THOSE VARIABLES. AN ANTI-FILTERING TECHNIQUE IS USED.
2119 !
2120 ! PROGRAM HISTORY LOG:
2121 ! 96-07-19 JANJIC - ORIGINATOR
2122 ! 05-02-03 GRELL,PECKHAM - MODIFIED FOR SCALARS
2123 !
2124 ! USAGE: CALL VAD2_SCAL FROM SUBROUTINE SOLVE_NMM
2125 ! INPUT ARGUMENT LIST:
2126 !
2127 ! OUTPUT ARGUMENT LIST
2128 !
2129 ! OUTPUT FILES:
2130 ! NONE
2131 ! SUBPROGRAMS CALLED:
2132 !
2133 ! UNIQUE: NONE
2134 !
2135 ! LIBRARY: NONE
2136 !
2137 ! ATTRIBUTES:
2138 ! LANGUAGE: FORTRAN 90
2139 ! MACHINE : IBM
2140 !$$$
2141 !***********************************************************************
2142 !----------------------------------------------------------------------
2143 !
2144 IMPLICIT NONE
2145 !
2146 !----------------------------------------------------------------------
2147 !
2148 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2149 & ,IMS,IME,JMS,JME,KMS,KME &
2150 ,ITS,ITE,JTS,JTE,KTS,KTE
2151 !
2152 INTEGER,INTENT(IN) :: LSTART,NUM_SCAL
2153 !
2154 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
2155 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: N_IUP_H,N_IUP_V &
2156 & ,N_IUP_ADH,N_IUP_ADV
2157 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V &
2158 & ,IUP_ADH,IUP_ADV
2159 !
2160 INTEGER,INTENT(IN) :: IDTAD,NTSD
2161 !
2162 REAL,INTENT(IN) :: DT,DY,PDTOP
2163 !
2164 REAL,DIMENSION(KMS:KME),INTENT(IN) :: AETA1,AETA2,DETA1,DETA2
2165 !
2166 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DX,HBM2,PDSL
2167 !
2168 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: PETDT
2169 !
2170 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME,1:NUM_SCAL) &
2171 ,INTENT(INOUT) :: SCAL
2172 !
2173 !----------------------------------------------------------------------
2174 !*** LOCAL VARIABLES
2175 !----------------------------------------------------------------------
2176 !
2177 REAL,PARAMETER :: FF1=0.500
2178 !
2179 LOGICAL,SAVE :: TRADITIONAL=.TRUE.
2180 !
2181 INTEGER :: I,IRECV,J,JFP,JFQ,K,L,LAP,LLAP
2182 !
2183 INTEGER,DIMENSION(KTS:KTE) :: LA
2184 !
2185 REAL*8 :: ADDT,AFRP,D2PQQ,DETAP,DPDN,DPUP,DQP &
2186 & ,HADDT,HBM2IJ &
2187 & ,Q00,Q4P,QP,QP0 &
2188 & ,RFACQK,RFC,RR &
2189 & ,SUMNQ,SUMPQ
2190 !
2191 REAL :: SFACQK
2192 !
2193 REAL,DIMENSION(KTS:KTE) :: AFR,DEL,DQL,DWL,E3,E4,PETDTK &
2194 & ,RFACE,RFACQ,RFACW,Q3,Q4,W3,W4
2195 !
2196 !-----------------------------------------------------------------------
2197 !***********************************************************************
2198 !-----------------------------------------------------------------------
2199 !
2200 ADDT=REAL(IDTAD)*DT
2201 !
2202 !-----------------------------------------------------------------------
2203 !$omp parallel do &
2204 !$omp& private(afr,afrp,d2pqq,detap,dpdn,dpup &
2205 !$omp& ,dql,dqp,haddt,i,j,k &
2206 !$omp& ,la,lap,llap,petdtk,q00,q3,q4,q4p,qp,qp0,rfacqk &
2207 !$omp& ,rfc,rr,sfacqk,sumnq,sumpq)
2208 !-----------------------------------------------------------------------
2209 !
2210 scalar_loop: DO L=LSTART,NUM_SCAL
2211 !
2212 main_integration: DO J=MYJS2,MYJE2
2213 !
2214 !-----------------------------------------------------------------------
2215 !
2216 main_iloop: DO I=MYIS1_P1,MYIE1_P1
2217 !
2218 !-----------------------------------------------------------------------
2219 !
2220 DO K=KTS,KTE
2221 Q3(K)=SCAL(I,J,K,L)
2222 Q4(K)=Q3(K)
2223 ENDDO
2224 !
2225 IF(TRADITIONAL)THEN
2226 PETDTK(KTE)=PETDT(I,J,KTE-1)*0.5
2227 !
2228 DO K=KTE-1,KTS+1,-1
2229 PETDTK(K)=(PETDT(I,J,K)+PETDT(I,J,K-1))*0.5
2230 ENDDO
2231 !
2232 PETDTK(KTS)=PETDT(I,J,KTS)*0.5
2233 !
2234 ELSE
2235 !
2236 !-----------------------------------------------------------------------
2237 !*** PERFORM HORIZONTAL AVERAGING OF VERTICAL VELOCITY
2238 !-----------------------------------------------------------------------
2239 !
2240 PETDTK(KTE)=(PETDT(I+IHW(J-1),J-1,KTE-1) &
2241 & +PETDT(I+IHE(J-1),J-1,KTE-1) &
2242 & +PETDT(I+IHW(J+1),J+1,KTE-1) &
2243 & +PETDT(I+IHE(J+1),J+1,KTE-1) &
2244 & +PETDT(I,J,KTE-1)*4. )*0.0625
2245 !
2246 DO K=KTE-1,KTS+1,-1
2247 PETDTK(K)=(PETDT(I+IHW(J-1),J-1,K-1) &
2248 +PETDT(I+IHE(J-1),J-1,K-1) &
2249 & +PETDT(I+IHW(J+1),J+1,K-1) &
2250 & +PETDT(I+IHE(J+1),J+1,K-1) &
2251 & +PETDT(I+IHW(J-1),J-1,K ) &
2252 & +PETDT(I+IHE(J-1),J-1,K ) &
2253 & +PETDT(I+IHW(J+1),J+1,K ) &
2254 & +PETDT(I+IHE(J+1),J+1,K ) &
2255 & +(PETDT(I,J,K-1)+PETDT(I,J,K))*4. &
2256 & )*0.0625
2257 ENDDO
2258 !
2259 PETDTK(KTS)=(PETDT(I+IHW(J-1),J-1,KTS) &
2260 & +PETDT(I+IHE(J-1),J-1,KTS) &
2261 & +PETDT(I+IHW(J+1),J+1,KTS) &
2262 & +PETDT(I+IHE(J+1),J+1,KTS) &
2263 & +PETDT(I,J,KTS)*4. )*0.0625
2265 ENDIF
2266 !
2267 !-----------------------------------------------------------------------
2268 !
2269 HADDT=-ADDT*HBM2(I,J)
2270 !
2271 DO K=KTE,KTS,-1
2272 RR=PETDTK(K)*HADDT
2273 !
2274 IF(RR<0.)THEN
2275 LAP=1
2276 ELSE
2277 LAP=-1
2278 ENDIF
2279 !
2280 LA(K)=LAP
2281 LLAP=K+LAP
2282 !
2283 IF(LLAP>KTS-1.AND.LLAP<KTE+1)THEN
2284 RR=ABS(RR/((AETA1(LLAP)-AETA1(K))*PDTOP &
2285 & +(AETA2(LLAP)-AETA2(K))*PDSL(I,J)))
2286 IF(RR>0.9)RR=0.9
2287 !
2288 AFR(K)=(((FF4*RR+FF3)*RR+FF2)*RR+FF1)*RR
2289 DQP=(Q3(LLAP)-Q3(K))*RR
2290 DQL(K)=DQP
2291 ELSE
2292 RR=0.
2293 AFR(K)=0.
2294 DQL(K)=0.
2295 ENDIF
2296 ENDDO
2297 !
2298 !-----------------------------------------------------------------------
2299 !
2300 IF(LA(KTE-1)>0)THEN
2301 RFC=(DETA1(KTE-1)*PDTOP+DETA2(KTE-1)*PDSL(I,J)) &
2302 & /(DETA1(KTE )*PDTOP+DETA2(KTE )*PDSL(I,J))
2303 DQL(KTE)=-DQL(KTE-1)*RFC
2304 ENDIF
2305 !
2306 IF(LA(KTS+1)<0)THEN
2307 RFC=(DETA1(KTS+1)*PDTOP+DETA2(KTS+1)*PDSL(I,J)) &
2308 & /(DETA1(KTS )*PDTOP+DETA2(KTS )*PDSL(I,J))
2309 DQL(KTS)=-DQL(KTS+1)*RFC
2310 ENDIF
2311 !
2312 DO K=KTS,KTE
2313 Q4(K)=Q3(K)+DQL(K)
2314 ENDDO
2315 !
2316 !-----------------------------------------------------------------------
2317 !*** ANTI-FILTERING STEP
2318 !-----------------------------------------------------------------------
2319 !
2320 SUMPQ=0.
2321 SUMNQ=0.
2322 !
2323 !*** ANTI-FILTERING LIMITERS
2324 !
2325 antifilter: DO K=KTE-1,KTS+1,-1
2326 !
2327 DETAP=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J)
2328 !
2329 Q4P=Q4(K)
2330 !
2331 LAP=LA(K)
2332 !
2333 DPDN=(AETA1(K+LAP)-AETA1(K))*PDTOP &
2334 & +(AETA2(K+LAP)-AETA2(K))*PDSL(I,J)
2335 DPUP=(AETA1(K)-AETA1(K-LAP))*PDTOP &
2336 & +(AETA2(K)-AETA2(K-LAP))*PDSL(I,J)
2337 !
2338 AFRP=2.*AFR(K)*DPDN*DPDN/(DPDN+DPUP)
2339 D2PQQ=((Q4(K+LAP)-Q4P)/DPDN &
2340 & -(Q4P-Q4(K-LAP))/DPUP)*AFRP
2341 !
2342 QP=Q4P-D2PQQ
2343 !
2344 Q00=Q3(K)
2345 QP0=Q3(K+LAP)
2346 !
2347 QP=MAX(QP,MIN(Q00,QP0))
2348 QP=MIN(QP,MAX(Q00,QP0))
2349 !
2350 DQP=QP-Q00
2351 !
2352 DQL(K)=DQP
2353 !
2354 ENDDO antifilter
2355 !
2356 !-----------------------------------------------------------------------
2357 !
2358 IF(LA(KTE-1)>0)THEN
2359 RFC=(DETA1(KTE-1)*PDTOP+DETA2(KTE-1)*PDSL(I,J)) &
2360 & /(DETA1(KTE )*PDTOP+DETA2(KTE )*PDSL(I,J))
2361 DQL(KTE)=-DQL(KTE-1)*RFC+DQL(KTE)
2362 ENDIF
2363 !
2364 IF(LA(KTS+1)<0)THEN
2365 RFC=(DETA1(KTS+1)*PDTOP+DETA2(KTS+1)*PDSL(I,J)) &
2366 & /(DETA1(KTS )*PDTOP+DETA2(KTS )*PDSL(I,J))
2367 DQL(KTS)=-DQL(KTS+1)*RFC+DQL(KTS)
2368 ENDIF
2369 !
2370 DO K=KTS,KTE
2371 DETAP=DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J)
2372 DQP=DQL(K)*DETAP
2373 !
2374 IF(DQP>0.)THEN
2375 SUMPQ=SUMPQ+DQP
2376 ELSE
2377 SUMNQ=SUMNQ+DQP
2378 ENDIF
2379 ENDDO
2380 !
2381 !-----------------------------------------------------------------------
2382 !*** FIRST MOMENT CONSERVING FACTOR
2383 !-----------------------------------------------------------------------
2384 !
2385 IF(SUMPQ>1.E-9)THEN
2386 SFACQK=-SUMNQ/SUMPQ
2387 ELSE
2388 SFACQK=1.
2389 ENDIF
2390 !
2391 IF(SFACQK<CONSERVE_MIN.OR.SFACQK>CONSERVE_MAX)SFACQK=1.
2392 !
2393 RFACQK=1./SFACQK
2394 !
2395 !-----------------------------------------------------------------------
2396 !*** IMPOSE CONSERVATION ON ANTI-FILTERING
2397 !-----------------------------------------------------------------------
2398 !
2399 DO K=KTE,KTS,-1
2400 DQP=DQL(K)
2401 IF(SFACQK>=1.)THEN
2402 IF(DQP<0.)DQP=DQP*RFACQK
2403 ELSE
2404 IF(DQP>0.)DQP=DQP*SFACQK
2405 ENDIF
2406 SCAL(I,J,K,L)=Q3(K)+DQP
2407 ENDDO
2408 !
2409 !-----------------------------------------------------------------------
2410 !
2411 ENDDO main_iloop
2412 !
2413 !-----------------------------------------------------------------------
2414 !
2415 ENDDO main_integration
2416 !
2417 !-----------------------------------------------------------------------
2418 !
2419 ENDDO scalar_loop
2420 !
2421 !-----------------------------------------------------------------------
2422 !
2423 END SUBROUTINE VAD2_SCAL
2424 !
2425 !-----------------------------------------------------------------------
2426 !
2427 !***********************************************************************
2428 SUBROUTINE HAD2_SCAL( &
2429 #if defined(DM_PARALLEL)
2430 & DOMDESC , &
2431 #endif
2432 & NTSD,DT,IDTAD,DX,DY &
2433 & ,AETA1,AETA2,DETA1,DETA2,PDSL,PDTOP &
2434 & ,HBM2,HBM3 &
2435 & ,SCAL,U,V,Z,HYDRO &
2436 & ,N_IUP_H,N_IUP_V &
2437 & ,N_IUP_ADH,N_IUP_ADV &
2438 & ,IUP_H,IUP_V,IUP_ADH,IUP_ADV &
2439 & ,IHE,IHW,IVE,IVW &
2440 & ,NUM_SCAL,LSTART &
2441 & ,IDS,IDE,JDS,JDE,KDS,KDE &
2442 & ,IMS,IME,JMS,JME,KMS,KME &
2443 & ,ITS,ITE,JTS,JTE,KTS,KTE)
2444 !***********************************************************************
2445 !$$$ SUBPROGRAM DOCUMENTATION BLOCK
2446 ! . . .
2447 ! SUBPROGRAM: HAD2_SCAL HORIZONTAL ADVECTION OF SCALAR
2448 ! PRGRMMR: JANJIC ORG: W/NP22 DATE: 96-07-19
2449 ! GRELL,PECKHAM ORG: NOAA/FSL DATE: 05-02-03
2450 !
2451 ! ABSTRACT:
2452 ! HAD2_SCAL CALCULATES THE CONTRIBUTION OF THE HORIZONTAL ADVECTION
2453 ! TO THE TENDENCIES OF SCALAR SUBSTANCES AND THEN
2454 ! UPDATES THOSE VARIABLES. AN ANTI-FILTERING TECHNIQUE IS USED.
2455 !
2456 ! PROGRAM HISTORY LOG:
2457 ! 96-07-19 JANJIC - ORIGINATOR
2458 ! 05-01-03 GRELL,PECKHAM - MODIFIED FOR SCALAR
2459 !
2460 ! USAGE: CALL HAD2_SCAL FROM SUBROUTINE SOLVE_NMM
2461 ! INPUT ARGUMENT LIST:
2462 !
2463 ! OUTPUT ARGUMENT LIST
2464 !
2465 ! OUTPUT FILES:
2466 ! NONE
2467 ! SUBPROGRAMS CALLED:
2468 !
2469 ! UNIQUE: NONE
2470 !
2471 ! LIBRARY: NONE
2472 !
2473 ! ATTRIBUTES:
2474 ! LANGUAGE: FORTRAN 90
2475 ! MACHINE : IBM
2476 !$$$
2477 !***********************************************************************
2478 !-----------------------------------------------------------------------
2479 !
2480 IMPLICIT NONE
2481 !
2482 !-----------------------------------------------------------------------
2483 !
2484 INTEGER,INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
2485 & ,IMS,IME,JMS,JME,KMS,KME &
2486 & ,ITS,ITE,JTS,JTE,KTS,KTE
2487 !
2488 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: IHE,IHW,IVE,IVW
2489 INTEGER,DIMENSION(JMS:JME),INTENT(IN) :: N_IUP_H,N_IUP_V &
2490 & ,N_IUP_ADH,N_IUP_ADV
2491 INTEGER,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: IUP_H,IUP_V &
2492 & ,IUP_ADH,IUP_ADV
2493 !
2494 !-----------------------------------------------------------------------
2495 !
2496 INTEGER,INTENT(IN) :: IDTAD,LSTART,NTSD,NUM_SCAL
2497 !
2498 REAL,INTENT(IN) :: DT,DY,PDTOP
2499 !
2500 REAL,DIMENSION(KMS:KME),INTENT(IN) :: AETA1,AETA2,DETA1,DETA2
2501 !
2502 REAL,DIMENSION(IMS:IME,JMS:JME),INTENT(IN) :: DX,HBM2,HBM3,PDSL
2503 !
2504 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME),INTENT(IN) :: U,V,Z
2505 !
2506 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME,1:NUM_SCAL) &
2507 ,INTENT(INOUT) :: SCAL
2508 !
2509 LOGICAL,INTENT(IN) :: HYDRO
2510 !
2511 !-----------------------------------------------------------------------
2512 !*** LOCAL VARIABLES
2513 !-----------------------------------------------------------------------
2514 !
2515 REAL,PARAMETER :: FF1=0.530
2516 !
2517 #ifdef DM_PARALLEL
2518 INTEGER :: DOMDESC
2519 #endif
2520 !
2521 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
2522 LOGICAL,EXTERNAL :: WRF_DM_ON_MONITOR
2523 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME,2) :: XSUMS_L
2524 REAL,DIMENSION(IDS:IDE,JDS:JDE,KDS:KDE,2) :: XSUMS_G
2525 #endif
2526 !
2527 LOGICAL :: BOT,TOP
2528 !
2529 INTEGER :: I,IRECV,J,JFP,JFQ,K,L,LAP,LLAP,MPI_COMM_COMP
2530 INTEGER :: N
2531 !
2532 INTEGER,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5,KTS:KTE) :: IFPA,IFPF &
2533 & ,IFQA,IFQF &
2534 & ,JFPA,JFPF &
2535 & ,JFQA,JFQF
2536 !
2537 REAL :: ADDT,AFRP,CRIT,D2PQE,D2PQQ,D2PQW,DEP,DESTIJ,DQP,DQSTIJ &
2538 & ,DVOLP,DWP,DWSTIJ,DZA,DZB,E00,E0Q,E1X,E2IJ,E4P,ENH,EP,EP0 &
2539 & ,ESTIJ,FPQ,HAFP,HAFQ,HBM2IJ,HM,PP,PPQ00,Q00,Q0Q &
2540 & ,Q1IJ,Q4P,QP,QP0,QSTIJ,RDY,RFACE,RFACQ,RFACW,RFC &
2541 & ,RFEIJ,RFQIJ,RFWIJ,RR,SLOPAC,SPP,SQP,SSA,SSB,SUMNQ,SUMPQ &
2542 & ,TTA,TTB,W00,W0Q,W1IJ,W4P,WP,WP0,WSTIJ
2543 !
2544 DOUBLE PRECISION,DIMENSION(2,KTS:KTE) :: GSUMS,XSUMS
2545 !
2546 REAL,DIMENSION(KTS:KTE) :: AFR,DEL,DQL,DWL,Q3,Q4
2547 !
2548 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: DARE,EMH
2549 !
2550 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5,KTS:KTE) :: AFP,AFQ,DEST &
2551 & ,DQST,DVOL,DWST &
2552 & ,Q1
2553 !
2554 REAL,DIMENSION(IMS:IME,JMS:JME,KMS:KME) :: Q
2555 !
2556 !-----------------------------------------------------------------------
2557 integer :: nunit,ier
2558 save nunit
2559 !-----------------------------------------------------------------------
2560 !***********************************************************************
2561 !-----------------------------------------------------------------------
2562 !
2563 RDY=1./DY
2564 SLOPAC=SLOPHT*SQRT(2.)*0.5*50.
2565 CRIT=SLOPAC*REAL(IDTAD)*DT*RDY*1000.
2566 !
2567 ADDT=REAL(IDTAD)*DT
2568 ENH=ADDT/(08.*DY)
2569 !
2570 !-----------------------------------------------------------------------
2571 !$omp parallel do &
2572 !$omp& private(i,j)
2573 DO J=MYJS_P3,MYJE_P3
2574 DO I=MYIS_P2,MYIE_P2
2575 EMH (I,J)=ADDT/(08.*DX(I,J))
2576 DARE(I,J)=HBM3(I,J)*DX(I,J)*DY
2577 ENDDO
2578 ENDDO
2579 !-----------------------------------------------------------------------
2580 !
2581 scalar_loop: DO L=LSTART,NUM_SCAL
2582 !
2583 !-----------------------------------------------------------------------
2584 !$omp parallel do &
2585 !$omp& private(dza,dzb,e1x,fpq,hm,i,j,jfp,jfq,k,pp,qp,ssa,ssb,spp,sqp &
2586 !$omp& ,tta,ttb)
2587 !-----------------------------------------------------------------------
2588 !
2589 vertical_1: DO K=KTS,KTE
2590 !
2591 !-----------------------------------------------------------------------
2592 !
2593 DO J=MYJS_P3,MYJE_P3
2594 DO I=MYIS_P2,MYIE_P2
2595 DVOL(I,J,K)=DARE(I,J)*(DETA1(K)*PDTOP+DETA2(K)*PDSL(I,J))
2596 Q (I,J,K)=SCAL(I,J,K,L)
2597 Q1(I,J,K)=Q(I,J,K)
2598 ENDDO
2599 ENDDO
2600 !
2601 !-----------------------------------------------------------------------
2602 !
2603 DO J=MYJS2_P1,MYJE2_P1
2604 DO I=MYIS1_P1,MYIE1_P1
2605 !
2606 HM=HBM2(I,J)
2607 TTA=(U(I,J-1,K)+U(I+IHW(J),J,K)+U(I+IHE(J),J,K)+U(I,J+1,K)) &
2608 & *EMH(I,J)*HM
2609 TTB=(V(I,J-1,K)+V(I+IHW(J),J,K)+V(I+IHE(J),J,K)+V(I,J+1,K)) &
2610 & *ENH*HBM2(I,J)
2611 !
2612 SPP=-TTA-TTB
2613 SQP= TTA-TTB
2614 !
2615 IF(SPP<0.)THEN
2616 JFP=-1
2617 ELSE
2618 JFP=1
2619 ENDIF
2620 IF(SQP<0.)THEN
2621 JFQ=-1
2622 ELSE
2623 JFQ=1
2624 ENDIF
2625 !
2626 IFPA(I,J,K)=IHE(J)+I+( JFP-1)/2
2627 IFQA(I,J,K)=IHE(J)+I+(-JFQ-1)/2
2628 !
2629 JFPA(I,J,K)=J+JFP
2630 JFQA(I,J,K)=J+JFQ
2631 !
2632 IFPF(I,J,K)=IHE(J)+I+(-JFP-1)/2
2633 IFQF(I,J,K)=IHE(J)+I+( JFQ-1)/2
2634 !
2635 JFPF(I,J,K)=J-JFP
2636 JFQF(I,J,K)=J-JFQ
2637 !
2638 !-----------------------------------------------------------------------
2639 IF(.NOT.HYDRO)THEN ! z currently not available for hydro=.true.
2640 DZA=(Z(IFPA(I,J,K),JFPA(I,J,K),K)-Z(I,J,K))*RDY
2641 DZB=(Z(IFQA(I,J,K),JFQA(I,J,K),K)-Z(I,J,K))*RDY
2642 !
2643 IF(ABS(DZA)>SLOPAC)THEN
2644 SSA=DZA*SPP
2645 IF(SSA>CRIT)THEN
2646 SPP=0. !spp*.1
2647 ENDIF
2648 ENDIF
2649 !
2650 IF(ABS(DZB)>SLOPAC)THEN
2651 SSB=DZB*SQP
2652 IF(SSB>CRIT)THEN
2653 SQP=0. !sqp*.1
2654 ENDIF
2655 ENDIF
2656 !
2657 ENDIF
2658 !
2659 !-----------------------------------------------------------------------
2660 !
2661 FPQ=SPP*SQP*0.25
2662 PP=ABS(SPP)
2663 QP=ABS(SQP)
2664 !
2665 AFP(I,J,K)=(((FF4*PP+FF3)*PP+FF2)*PP+FF1)*PP
2666 AFQ(I,J,K)=(((FF4*QP+FF3)*QP+FF2)*QP+FF1)*QP
2667 !
2668 Q1(I,J,K)=(Q (IFPA(I,J,K),JFPA(I,J,K),K)-Q (I,J,K))*PP &
2669 & +(Q (IFQA(I,J,K),JFQA(I,J,K),K)-Q (I,J,K))*QP &
2670 & +(Q (I,J-2,K)+Q (I,J+2,K) &
2671 & -Q (I-1,J,K)-Q (I+1,J,K))*FPQ &
2672 & +Q(I,J,K)
2673 !
2674 ENDDO
2675 ENDDO
2676 !
2677 !-----------------------------------------------------------------------
2678 !
2679 ENDDO vertical_1
2680 !
2681 !-----------------------------------------------------------------------
2682 !*** ANTI-FILTERING STEP
2683 !-----------------------------------------------------------------------
2684 !
2685 DO K=KTS,KTE
2686 XSUMS(1,K)=0.
2687 XSUMS(2,K)=0.
2688 ENDDO
2689 !
2690 !-----------------------------------------------------------------------
2691 !
2692 !*** ANTI-FILTERING LIMITERS
2693 !
2694 !-----------------------------------------------------------------------
2695 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
2696 DO N=1,2
2697 !
2698 !$omp parallel do &
2699 !$omp& private(i,j,k)
2700 DO K=KMS,KME
2701 DO J=JMS,JME
2702 DO I=IMS,IME
2703 XSUMS_L(I,J,K,N)=0.
2704 ENDDO
2705 ENDDO
2706 ENDDO
2707 !
2708 !$omp parallel do &
2709 !$omp& private(i,j,k)
2710 DO K=KDS,KDE
2711 DO J=JDS,JDE
2712 DO I=IDS,IDE
2713 XSUMS_G(I,J,K,N)=0.
2714 ENDDO
2715 ENDDO
2716 ENDDO
2717 !
2718 ENDDO
2719 !
2720 #endif
2721 !-----------------------------------------------------------------------
2722 !$omp parallel do &
2723 !$omp& private(d2pqe,d2pqq,d2pqw,destij,dqstij,dvolp,dwstij &
2724 !$omp& ,e00,e0q,ep0,estij,hafp,hafq,i,j,k &
2725 !$omp& ,q00,q0q,q1ij,qp0,qstij,w00,w0q,wp0,wstij)
2726 !-----------------------------------------------------------------------
2727 !
2728 vertical_2: DO K=KTS,KTE
2729 !
2730 !-----------------------------------------------------------------------
2731 !
2732 DO J=MYJS2,MYJE2
2733 DO I=MYIS1,MYIE1
2734 !
2735 DVOLP=DVOL(I,J,K)
2736 Q1IJ =Q1(I,J,K)
2737 !
2738 HAFP=AFP(I,J,K)
2739 HAFQ=AFQ(I,J,K)
2740 !
2741 D2PQQ=(Q1(IFPA(I,J,K),JFPA(I,J,K),K)-Q1IJ &
2742 & -Q1IJ+Q1(IFPF(I,J,K),JFPF(I,J,K),K)) &
2743 & *HAFP &
2744 & +(Q1(IFQA(I,J,K),JFQA(I,J,K),K)-Q1IJ &
2745 & -Q1IJ+Q1(IFQF(I,J,K),JFQF(I,J,K),K)) &
2746 & *HAFQ
2747 !
2748 QSTIJ=Q1IJ-D2PQQ
2749 !
2750 Q00=Q (I ,J ,K)
2751 QP0=Q (IFPA(I,J,K),JFPA(I,J,K),K)
2752 Q0Q=Q (IFQA(I,J,K),JFQA(I,J,K),K)
2753 !
2754 QSTIJ=MAX(QSTIJ,MIN(Q00,QP0,Q0Q))
2755 QSTIJ=MIN(QSTIJ,MAX(Q00,QP0,Q0Q))
2756 !
2757 DQSTIJ=QSTIJ-Q(I,J,K)
2758 !
2759 DQST(I,J,K)=DQSTIJ
2760 !
2761 DQSTIJ=DQSTIJ*DVOLP
2762 !
2763 !-----------------------------------------------------------------------
2764 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
2765 !-----------------------------------------------------------------------
2766 DO N=1,2
2767 XSUMS_L(I,J,K,N)=0.
2768 ENDDO
2769 !
2770 IF(DQSTIJ>0.)THEN
2771 XSUMS_L(I,J,K,1)=DQSTIJ
2772 ELSE
2773 XSUMS_L(I,J,K,2)=DQSTIJ
2774 ENDIF
2775 !
2776 !-----------------------------------------------------------------------
2777 #else
2778 !-----------------------------------------------------------------------
2779 IF(DQSTIJ>0.)THEN
2780 XSUMS(1,K)=XSUMS(1,K)+DQSTIJ
2781 ELSE
2782 XSUMS(2,K)=XSUMS(2,K)+DQSTIJ
2783 ENDIF
2784 !
2785 !-----------------------------------------------------------------------
2786 #endif
2787 !-----------------------------------------------------------------------
2788 !
2789 ENDDO
2790 ENDDO
2791 !
2792 !-----------------------------------------------------------------------
2793 !
2794 ENDDO vertical_2
2795 !
2796 !-----------------------------------------------------------------------
2797 #if defined(BIT_FOR_BIT) && defined(DM_PARALLEL) && !defined(STUBMPI)
2798 !-----------------------------------------------------------------------
2799 DO N=1,2
2800 CALL WRF_PATCH_TO_GLOBAL_REAL(XSUMS_L(IMS,JMS,KMS,N) &
2801 &, XSUMS_G(1,1,1,N),DOMDESC &
2802 &, 'xyz','xzy' &
2803 &, IDS,IDE,KDS,KDE,JDS,JDE &
2804 &, IMS,IME,KMS,KME,JMS,JME &
2805 &, ITS,ITE,KTS,KTE,JTS,JTE)
2806 ENDDO
2807 !
2808 DO K=KTS,KTE
2809 DO N=1,2
2810 GSUMS(N,K)=0.
2811 ENDDO
2812 ENDDO
2813 !
2814 IF(WRF_DM_ON_MONITOR())THEN
2815 DO N=1,2
2816 !$omp parallel do &
2817 !$omp& private(i,j,k)
2818 DO K=KDS,KDE
2819 DO J=JDS,JDE
2820 DO I=IDS,IDE
2821 GSUMS(N,K)=GSUMS(N,K)+XSUMS_G(I,J,K,N)
2822 ENDDO
2823 ENDDO
2824 ENDDO
2825 ENDDO
2826 ENDIF
2828 CALL WRF_DM_BCAST_BYTES(GSUMS,2*RWORDSIZE*2*(KDE-KDS+1) )
2830 !-----------------------------------------------------------------------
2831 #else
2832 !-----------------------------------------------------------------------
2833 !
2834 !-----------------------------------------------------------------------
2835 !*** GLOBAL REDUCTION
2836 !-----------------------------------------------------------------------
2837 !
2838 # if defined(DM_PARALLEL) && !defined(STUBMPI)
2839 CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
2840 CALL MPI_ALLREDUCE(XSUMS,GSUMS,2*(KTE-KTS+1) &
2841 & ,MPI_DOUBLE_PRECISION,MPI_SUM &
2842 & ,MPI_COMM_COMP,IRECV)
2843 # else
2844 DO K=KTS,KTE
2845 DO N=1,2
2846 GSUMS(N,K)=XSUMS(N,K)
2847 ENDDO
2848 ENDDO
2849 # endif
2850 !
2851 !-----------------------------------------------------------------------
2852 #endif
2853 !-----------------------------------------------------------------------
2854 !
2855 !-----------------------------------------------------------------------
2856 !*** END OF GLOBAL REDUCTION
2857 !-----------------------------------------------------------------------
2858 !
2859 ! if(mype==0)then
2860 !!! if(ntsd==0)then
2861 !!! call int_get_fresh_handle(nunit)
2862 !!! close(nunit)
2863 ! nunit=56
2864 !!! open(unit=nunit,file='gsums',form='unformatted',iostat=ier)
2865 !!! endif
2866 ! endif
2867 !-----------------------------------------------------------------------
2868 !$omp parallel do &
2869 !$omp& private(destij,dqstij,dwstij,i,j,k,rface,rfacq,rfacw &
2870 !$omp& ,rfeij,rfqij,rfwij,sumnq,sumpq)
2871 !-----------------------------------------------------------------------
2872 !
2873 vertical_3: DO K=KTS,KTE
2874 !
2875 !-----------------------------------------------------------------------
2876 ! if(mype==0)then
2877 ! write(nunit)(gsums(i,k),i=1,6)
2878 ! endif
2879 !!! read(nunit)(gsums(i,k),i=1,6)
2880 !-----------------------------------------------------------------------
2881 !
2882 SUMPQ=GSUMS(1,K)
2883 SUMNQ=GSUMS(2,K)
2884 !
2885 !-----------------------------------------------------------------------
2886 !*** FIRST MOMENT CONSERVING FACTOR
2887 !-----------------------------------------------------------------------
2888 !
2889 IF(SUMPQ>1.)THEN
2890 RFACQ=-SUMNQ/SUMPQ
2891 ELSE
2892 RFACQ=1.
2893 ENDIF
2894 !
2895 IF(RFACQ<CONSERVE_MIN.OR.RFACQ>CONSERVE_MAX)RFACQ=1.
2896 !
2897 !-----------------------------------------------------------------------
2898 ! if(mype==0.and.ntsd==181)close(nunit)
2899 !-----------------------------------------------------------------------
2900 !
2901 !-----------------------------------------------------------------------
2902 !*** IMPOSE CONSERVATION ON ANTI-FILTERING
2903 !-----------------------------------------------------------------------
2904 !
2905 DO J=MYJS2,MYJE2
2906 IF(RFACQ<1.)THEN
2907 DO I=MYIS1,MYIE1
2908 DQSTIJ=DQST(I,J,K)
2909 RFQIJ=HBM2(I,J)*(RFACQ-1.)+1.
2910 IF(DQSTIJ>=0.)DQSTIJ=DQSTIJ*RFQIJ
2911 Q(I,J,K)=Q(I,J,K)+DQSTIJ
2912 ENDDO
2913 ELSE
2914 DO I=MYIS1,MYIE1
2915 DQSTIJ=DQST(I,J,K)
2916 RFQIJ=HBM2(I,J)*(RFACQ-1.)+1.
2917 IF(DQSTIJ<0.)DQSTIJ=DQSTIJ/RFQIJ
2918 Q(I,J,K)=Q(I,J,K)+DQSTIJ
2919 ENDDO
2920 ENDIF
2921 ENDDO
2922 !
2923 !-----------------------------------------------------------------------
2924 !
2925 DO J=MYJS,MYJE
2926 DO I=MYIS,MYIE
2927 SCAL(I,J,K,L)=Q(I,J,K)
2928 ENDDO
2929 ENDDO
2930 !
2931 !-----------------------------------------------------------------------
2932 !
2933 ENDDO vertical_3
2934 !
2935 !-----------------------------------------------------------------------
2936 !
2937 ENDDO scalar_loop
2938 !
2939 !-----------------------------------------------------------------------
2940 !
2941 END SUBROUTINE HAD2_SCAL
2942 !
2943 !-----------------------------------------------------------------------
2944 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2945 subroutine adv2 &
2946 (UPSTRM &
2947 ,mype,kss,kse &
2948 ,ids,ide,jds,jde,kds,kde &
2949 ,ims,ime,jms,jme,kms,kme &
2950 ,its,ite,jts,jte,kts,kte &
2951 ,N_IUP_H &
2952 ,N_IUP_ADH &
2953 ,IUP_H,IUP_ADH &
2954 ,ENT &
2955 ,idtad &
2956 ,dt,pdtop &
2957 ,ihe,ihw,ive,ivw &
2958 ,deta1,deta2 &
2959 ,EMT_LOC &
2960 ,fad,hbm2,pdsl,pdslo &
2961 ,petdt &
2962 ,UOLD,VOLD &
2963 ,s,sp &
2964 !---temporary arguments-------------------------------------------------
2965 ,fne,fse,few,fns,s1,tcs)
2966 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
2967 implicit none
2968 !-----------------------------------------------------------------------
2969 real,parameter:: &
2970 cfc=1.533 & ! adams-bashforth positioning in time
2971 ,bfc=1.-cfc & ! adams bashforth positioning in time
2972 ,cflc=9.005 & !
2973 ,epsq=1.e-20 & ! floor value for specific humidity
2974 ,epsq2=0.2 & ! floor value for 2tke
2975 ,epscm=2.e-6 & ! a floor value (not used)
2976 ,w1=1.0 & ! crank-nicholson uncentering
2977 !,w1=-1.00 & ! crank-nicholson uncentering
2978 ,w2=2.-w1 ! crank-nicholson uncentering
2980 logical,intent(in):: &
2981 upstrm
2983 integer,intent(in):: &
2984 idtad & ! time step multiplier
2985 ,kse & ! terminal species index
2986 ,kss & ! initial species index
2987 ,mype & !
2988 ,ids,ide,jds,jde,kds,kde &
2989 ,ims,ime,jms,jme,kms,kme &
2990 ,its,ite,jts,jte,kts,kte
2992 real,intent(in):: &
2993 ent & !
2994 ,dt & ! dynamics time step
2995 ,pdtop !
2997 real,dimension(kts:kte),intent(in):: &
2998 deta1 & ! delta sigmas
2999 ,deta2 ! delta pressures
3001 integer,dimension(jms:jme),intent(in):: &
3002 ihe,ihw,ive,ivw &
3003 ,n_iup_adh,n_iup_h
3005 integer,dimension(ims:ime,jms:jme),intent(in):: &
3006 iup_h,iup_adh
3008 real,dimension(2600),intent(in):: & !!!zj see nmm_max_dim in adve !!!zj
3009 emt_loc
3011 real,dimension(ims:ime,jms:jme),intent(in):: &
3012 fad & !
3013 ,hbm2 & !
3014 ,pdsl & ! sigma range pressure difference
3015 ,pdslo ! sigma range pressure difference
3017 real,dimension(ims:ime,jms:jme,kms:kme),intent(in):: &
3018 petdt & ! vertical mass flux
3019 ,uold,vold
3021 real,dimension(ims:ime,jms:jme,kms:kme,kss:kse),intent(inout):: &
3022 s ! tracers
3024 real,dimension(ims:ime,jms:jme,kms:kme,kss:kse),intent(inout):: &
3025 sp ! s at previous time level
3027 !---temporary arguments-------------------------------------------------
3028 real,dimension(ims:ime,jms:jme,kms:kme),intent(in):: &
3029 fne & ! mass flux, ne direction
3030 ,fse & ! mass flux, se direction
3031 ,few & ! mass flux, x direction
3032 ,fns ! mass flux, y direction
3034 real,dimension(ims:ime,jms:jme,kms:kme,kss:kse),intent(inout):: &
3035 s1 & ! intermediate value of s
3036 ,tcs ! timechange of s
3038 !--local variables------------------------------------------------------
3039 integer:: &
3040 i & !
3041 ,j & !
3042 ,k & !
3043 ,ks !
3045 INTEGER :: IEND,IFP,IFQ,II,IPQ,ISP,ISQ,ISTART &
3046 & ,IUP_ADH_J &
3047 & ,J1,JA,JAK,JEND,JGLOBAL,JJ,JKNT,JP2,JSTART &
3048 & ,KNTI_ADH,KSTART,KSTOP &
3049 & ,N,N_IUPH_J,N_IUPADH_J,N_IUPADV_J
3051 INTEGER :: MY_IS_GLB,MY_IE_GLB,MY_JS_GLB,MY_JE_GLB
3053 real:: &
3054 cf & ! temporary
3055 ,cms & ! temporary
3056 ,dtq & ! dt/4
3057 ,fahp & ! temporary grid factor
3058 ,sn & !
3059 ,rdp & ! 1/deltap
3060 ,vvlo & ! vertical velocity, lower interface
3061 ,vvup & ! vertical velocity, upper interface
3062 ,pvvup ! vertical mass flux, upper interface
3064 REAL :: ARRAY3_X &
3065 & ,F0,F1,F2,F3 &
3066 & ,PP &
3067 & ,QP &
3068 & ,TEMPA,TEMPB,TTA,TTB
3070 real,dimension(kts:kte):: &
3071 deta1_pdtop !
3073 INTEGER,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ISPA,ISQA
3075 real,dimension(its-5:ite+5,jts-5:jte+5):: &
3076 pdop & ! hydrostatic pressure difference at h points
3077 ,pvvlo & ! vertical mass flux, lower interface
3078 ,ss1 & ! extrapolated species between time levels
3079 ,ssne & ! flux, ne direction
3080 ,ssse & ! flux, nw direction
3081 ,ssx & ! flux, x direction
3082 ,ssy ! flux, y direction
3084 REAL,DIMENSION(ITS-5:ITE+5,JTS-5:JTE+5) :: ARRAY0,ARRAY1 &
3085 & ,ARRAY2,ARRAY3
3087 real,dimension(its-5:ite+5,jts-5:jte+5,kts:kte):: &
3088 crs & ! vertical advection temporary
3089 ,rcms ! vertical advection temporary
3091 real,dimension(its-5:ite+5,jts-5:jte+5,kts:kte,kss:kse):: &
3092 rsts ! vertical advection temporary
3093 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3094 DO J=JTS-5,JTE+5
3095 DO I=ITS-5,ITE+5
3096 pdop (i,j)=0.
3097 pvvlo(i,j)=0.
3098 ss1 (i,j)=0.
3099 ssne (i,j)=0.
3100 ssse (i,j)=0.
3101 enddo
3102 enddo
3103 !
3104 DO K=KTS,KTE
3105 DO J=JTS-5,JTE+5
3106 DO I=ITS-5,ITE+5
3107 crs (i,j,k)=0.
3108 rcms(i,j,k)=0.
3109 enddo
3110 enddo
3111 enddo
3112 !
3113 do ks=kss,kse
3114 DO K=KTS,KTE
3115 DO J=JTS-5,JTE+5
3116 DO I=ITS-5,ITE+5
3117 rsts(i,j,k,ks)=0.
3118 s1 (i,j,k,ks)=0.
3119 enddo
3120 enddo
3121 enddo
3122 enddo
3123 !
3124 do ks=kss,kse
3125 DO K=KMS,KME
3126 DO J=JMS,JME
3127 DO I=IMS,IME
3128 s1 (i,j,k,ks)=0.
3129 tcs(i,j,k,ks)=0.
3130 enddo
3131 enddo
3132 enddo
3133 enddo
3134 !-----------------------------------------------------------------------
3135 do k=kts,kte
3136 deta1_pdtop(k)=deta1(k)*pdtop
3137 enddo
3138 !-----------------------------------------------------------------------
3139 do ks=kss,kse ! loop by species
3140 !-----------------------------------------------------------------------
3141 DO K=KTS,KTE
3142 DO J=MYJS_P4,MYJE_P4
3143 DO I=MYIS_P4,MYIE_P4
3144 s1(i,j,k,ks)=sqrt(s(i,j,k,ks))
3145 enddo
3146 enddo
3147 enddo
3148 !-----------------------------------------------------------------------
3149 enddo ! end of the loop by species
3150 !-----------------------------------------------------------------------
3151 DO J=MYJS_P4,MYJE_P4
3152 DO I=MYIS_P4,MYIE_P4
3153 pdop(i,j)=(pdslo(i,j)+pdsl(i,j))*0.5
3154 enddo
3155 enddo
3156 !---crank-nicholson vertical advection----------------------------------
3157 dtq=dt*idtad*0.25
3159 DO J=MYJS2,MYJE2
3160 DO I=MYIS1,MYIE1
3161 pvvlo(i,j)=petdt(i,j,kte-1)*dtq
3162 vvlo=pvvlo(i,j)/(deta2(kte)*pdop(i,j)+deta1_pdtop(kte))
3163 !
3164 cms=-vvlo*w2+1.
3165 rcms(i,j,kte)=1./cms
3166 crs(i,j,kte)=vvlo*w2
3167 !
3168 do ks=kss,kse
3169 rsts(i,j,kte,ks)=(-vvlo*w1) &
3170 *(s1(i,j,kte-1,ks)-s1(i,j,kte,ks)) &
3171 +s1(i,j,kte,ks)
3172 enddo
3173 enddo
3174 enddo
3175 DO K=KTE-1,KTS+1,-1
3176 DO J=MYJS2,MYJE2
3177 DO I=MYIS1,MYIE1
3178 rdp=1./(deta2(k)*pdop(i,j)+deta1_pdtop(k))
3179 pvvup=pvvlo(i,j)
3180 pvvlo(i,j)=petdt(i,j,k-1)*dtq
3181 !
3182 vvup=pvvup*rdp
3183 vvlo=pvvlo(i,j)*rdp
3184 !
3185 ! if(abs(vvlo).gt.cflc) then
3186 ! if(vvlo.lt.0.) then
3187 ! vvlo=-cflc
3188 ! else
3189 ! vvlo= cflc
3190 ! endif
3191 ! endif
3192 !
3193 cf=-vvup*w2*rcms(i,j,k+1)
3194 cms=-crs(i,j,k+1)*cf+((vvup-vvlo)*w2+1.)
3195 rcms(i,j,k)=1./cms
3196 crs(i,j,k)=vvlo*w2
3197 !
3198 do ks=kss,kse
3199 rsts(i,j,k,ks)=-rsts(i,j,k+1,ks)*cf+s1(i,j,k,ks) &
3200 -(s1(i,j,k ,ks)-s1(i,j,k+1,ks))*vvup*w1 &
3201 -(s1(i,j,k-1,ks)-s1(i,j,k ,ks))*vvlo*w1
3202 enddo
3203 enddo
3204 enddo
3205 enddo
3206 DO J=MYJS2,MYJE2
3207 DO I=MYIS1,MYIE1
3208 pvvup=pvvlo(i,j)
3209 vvup=pvvup/(deta2(kts)*pdop(i,j)+deta1_pdtop(kts))
3210 !
3211 cf=-vvup*w2*rcms(i,j,kts+1)
3212 cms=-crs(i,j,kts+1)*cf+(vvup*w2+1.)
3213 rcms(i,j,kts)=1./cms
3214 crs(i,j,kts)=0.
3215 !
3216 do ks=kss,kse
3217 rsts(i,j,kts,ks)=-rsts(i,j,kts+1,ks)*cf+s1(i,j,kts,ks) &
3218 -(s1(i,j,kts,ks)-s1(i,j,kts+1,ks))*vvup*w1
3219 !
3220 tcs(i,j,kts,ks)=rsts(i,j,kts,ks)*rcms(i,j,kts)-s1(i,j,kts,ks)
3221 enddo
3222 enddo
3223 enddo
3224 do ks=kss,kse
3225 DO K=KTS+1,KTE
3226 DO J=MYJS2,MYJE2
3227 DO I=MYIS1,MYIE1
3228 tcs(i,j,k,ks)=(-crs(i,j,k)*(s1(i,j,k-1,ks)+tcs(i,j,k-1,ks)) &
3229 +rsts(i,j,k,ks)) &
3230 *rcms(i,j,k)-s1(i,j,k,ks)
3231 enddo
3232 enddo
3233 enddo
3234 enddo
3235 !-----------------------------------------------------------------------
3236 do ks=kss,kse ! loop by species
3237 !-----------------------------------------------------------------------
3238 DO K=KTS,KTE
3239 DO J=MYJS_P5,MYJE_P5
3240 DO I=MYIS_P5,MYIE_P5
3241 ss1(i,j)=s1(i,j,k,ks)*cfc+sp(i,j,k,ks)*bfc
3242 sp(i,j,k,ks)=s1(i,j,k,ks)
3243 enddo
3244 enddo
3245 !---fluxes--------------------------------------------------------------
3246 DO J=MYJS1_P2,MYJE1_P2
3247 DO I=MYIS_P2,MYIE_P3
3248 ssx(i,j)=(ss1(i+ive(j),j )-ss1(i+ivw(j),j ))*few(i,j,k) &
3249 *hbm2(i,j)
3250 ssy(i,j)=(ss1(i ,j+1)-ss1(i ,j-1))*fns(i,j,k) &
3251 *hbm2(i,j)
3252 enddo
3253 enddo
3254 DO J=MYJS1_P2,MYJE2_P2
3255 DO I=MYIS_P2,MYIE_P2
3256 ssne(i,j)=(ss1(i+ihe(j),j+1)-ss1(i,j))*fne(i,j,k)*hbm2(i,j)
3257 enddo
3258 enddo
3259 DO J=MYJS2_P2,MYJE1_P2
3260 DO I=MYIS_P2,MYIE_P2
3261 ssse(i,j)=(ss1(i+ihe(j),j-1)-ss1(i,j))*fse(i,j,k)*hbm2(i,j)
3262 enddo
3263 enddo
3264 !---advection of species------------------------------------------------
3265 DO J=MYJS5,MYJE5
3266 DO I=MYIS2,MYIE2
3267 tcs(i,j,k,ks)=((ssx (i+ihw(j),j )+ssx (i+ihe(j),j ) &
3268 +ssy (i ,j-1)+ssy (i ,j+1) &
3269 +ssne(i+ihw(j),j-1)+ssne(i ,j ) &
3270 +ssse(i ,j )+ssse(i+ihw(j),j+1)) &
3271 *fad(i,j)*2.0*idtad & !! 2.0 compensates for fad
3272 /(deta2(k)*pdop(i,j)+deta1_pdtop(k)) &
3273 +tcs(i,j,k,ks))*hbm2(i,j)
3274 enddo
3275 enddo
3276 !-----------------------------------------------------------------------
3277 !
3278 !*** upstream advection
3279 !
3280 !-----------------------------------------------------------------------
3281 !
3282 upstream: IF(UPSTRM)THEN
3283 !
3284 !-----------------------------------------------------------------------
3285 !***
3286 !*** COMPUTE UPSTREAM COMPUTATIONS ON THIS TASK'S ROWS.
3287 !***
3288 !-----------------------------------------------------------------------
3289 !
3290 jloop_upstream: DO J=MYJS2,MYJE2
3291 !
3292 N_IUPH_J=N_IUP_H(J) ! See explanation in START_DOMAIN_NMM
3293 DO II=0,N_IUPH_J-1
3294 !
3295 I=IUP_H(IMS+II,J)
3296 tta=emt_loc(j) &
3297 *(uold(i ,j-1,k)+uold(i+ihw(j),j ,k) &
3298 +uold(i+ihe(j),j ,k)+uold(i ,j+1,k))
3299 ttb=ent &
3300 *(vold(i ,j-1,k)+vold(i+ihw(j),j ,k) &
3301 +vold(i+ihe(j),j ,k)+vold(i ,j+1,k))
3302 PP=-TTA-TTB
3303 QP= TTA-TTB
3304 !
3305 IF(PP<0.)THEN
3306 ISPA(I,J)=-1
3307 ELSE
3308 ISPA(I,J)= 1
3309 ENDIF
3310 !
3311 IF(QP<0.)THEN
3312 ISQA(I,J)=-1
3313 ELSE
3314 ISQA(I,J)= 1
3315 ENDIF
3316 !
3317 PP=ABS(PP)
3318 QP=ABS(QP)
3319 ARRAY3_X=PP*QP
3320 ARRAY0(I,J)=ARRAY3_X-PP-QP
3321 ARRAY1(I,J)=PP-ARRAY3_X
3322 ARRAY2(I,J)=QP-ARRAY3_X
3323 ARRAY3(I,J)=ARRAY3_X
3324 ENDDO
3325 !
3326 !-----------------------------------------------------------------------
3327 !
3328 N_IUPADH_J=N_IUP_ADH(J)
3329 KNTI_ADH=1
3330 IUP_ADH_J=IUP_ADH(IMS,J)
3331 !
3332 iloop_T: DO II=0,N_IUPH_J-1
3333 !
3334 I=IUP_H(IMS+II,J)
3335 !
3336 ISP=ISPA(I,J)
3337 ISQ=ISQA(I,J)
3338 IFP=(ISP-1)/2
3339 IFQ=(-ISQ-1)/2
3340 IPQ=(ISP-ISQ)/2
3341 !
3342 !-----------------------------------------------------------------------
3343 !
3344 IF(I==IUP_ADH_J)THEN ! Upstream advection T tendencies
3345 !
3346 ISP=ISPA(I,J)
3347 ISQ=ISQA(I,J)
3348 IFP=(ISP-1)/2
3349 IFQ=(-ISQ-1)/2
3350 IPQ=(ISP-ISQ)/2
3351 !
3352 F0=ARRAY0(I,J)
3353 F1=ARRAY1(I,J)
3354 F2=ARRAY2(I,J)
3355 F3=ARRAY3(I,J)
3356 !
3357 tcs(i,j,k,ks)=(f0*s1(i,j,k,ks) &
3358 & +f1*s1(i+ihe(j)+ifp,j+isp,k,ks) &
3359 & +f2*s1(i+ihe(j)+ifq,j+isq,k,ks) &
3360 & +f3*s1(i+ipq,j+isp+isq,k,ks))*2.0 &
3361 & *idtad &
3362 & +tcs(i,j,k,ks)*hbm2(i,j)
3363 !
3364 !-----------------------------------------------------------------------
3365 !
3366 IF(KNTI_ADH<N_IUPADH_J)THEN
3367 IUP_ADH_J=IUP_ADH(IMS+KNTI_ADH,J)
3368 KNTI_ADH=KNTI_ADH+1
3369 ENDIF
3370 !
3371 ENDIF ! End of upstream advection tendency IF block
3372 !
3373 ENDDO iloop_T
3374 !
3375 !-----------------------------------------------------------------------
3376 enddo jloop_upstream
3377 !-----------------------------------------------------------------------
3378 endif upstream
3379 !-----------------------------------------------------------------------
3380 enddo ! kts,kte
3381 !-----------------------------------------------------------------------
3382 enddo ! end of the loop by the species
3383 !-----------------------------------------------------------------------
3384 endsubroutine adv2
3385 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3386 subroutine mono &
3387 ( &
3388 #if defined(DM_PARALLEL)
3389 domdesc, &
3390 #endif
3391 mype,ntsd,hours,kss,kse &
3392 ,ids,ide,jds,jde,kds,kde &
3393 ,ims,ime,jms,jme,kms,kme &
3394 ,its,ite,jts,jte,kts,kte &
3395 ,idtad &
3396 ,dy,pdtop &
3397 ,sumdrrw &
3398 ,ihe,ihw &
3399 ,deta1,deta2 &
3400 ,dx,hbm2,pd &
3401 ,s &
3402 !---temporary arguments-------------------------------------------------
3403 ,s1,tcs)
3404 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3405 implicit none
3406 !-----------------------------------------------------------------------
3407 real,parameter:: &
3408 epsq=1.e-20 & ! floor value for specific humidity
3409 ,epsq2=0.02 ! floor value for 2tke
3411 #ifdef DM_PARALLEL
3412 INTEGER :: DOMDESC
3413 #endif
3415 integer,intent(in):: &
3416 idtad & ! time step multiplier
3417 ,kse & ! terminal species index
3418 ,kss & ! initial species index
3419 ,mype & !
3420 ,ntsd & !
3421 ,ids,ide,jds,jde,kds,kde &
3422 ,ims,ime,jms,jme,kms,kme &
3423 ,its,ite,jts,jte,kts,kte
3425 real,intent(in):: &
3426 dy & !
3427 ,pdtop & !
3428 ,hours !
3430 real,intent(inout):: &
3431 sumdrrw
3433 integer,dimension(jms:jme),intent(in):: &
3434 ihe,ihw
3436 real,dimension(kts:kte),intent(in):: &
3437 deta1 & ! delta sigmas
3438 ,deta2 ! delta sigmas
3440 real,dimension(ims:ime,jms:jme),intent(in):: &
3441 dx & !
3442 ,hbm2 & !
3443 ,pd ! sigma range pressure difference
3445 real,dimension(ims:ime,jms:jme,kms:kme,kss:kse),intent(in):: &
3446 s ! s
3447 !---temporary arguments-------------------------------------------------
3448 real,dimension(ims:ime,jms:jme,kms:kme,kss:kse),intent(inout):: &
3449 s1 & ! intermediate value of s
3450 ,tcs ! timechange of s
3451 !--local variables------------------------------------------------------
3452 integer:: &
3453 i & !
3454 ,ierr & !
3455 ,irecv & !
3456 ,j & !
3457 ,k & !
3458 ,ks & !
3459 ,lngth & !
3460 ,mpi_comm_comp !
3462 real:: &
3463 dsks & !
3464 ,dvolp & !
3465 ,rfacs & !
3466 ,s1p & !
3467 ,sfacs & !
3468 ,smax & ! local maximum
3469 ,smin & ! local minimum
3470 ,smaxh & ! horizontal local maximum
3471 ,sminh & ! horizontal local minimum
3472 ,smaxv & ! vertical local maximum
3473 ,sminv & ! vertical local minimum
3474 ,sn & !
3475 ,sumns & !
3476 ,sumps & !
3477 ,steep
3479 double precision:: &
3480 xsmp,gsmp
3482 double precision,dimension(2*kss-1:2*kse):: &
3483 vgsms
3485 real,dimension(kts:kte):: &
3486 deta1_pdtop !
3488 real,dimension(2*kss-1:2*kse,kts:kte):: &
3489 gsms_single !
3491 double precision,dimension(2*kss-1:2*kse,kts:kte):: &
3492 gsms & ! sum of neg/pos changes all global fields
3493 ,xsms ! sum of neg/pos changes all local fields
3495 double precision,dimension(2*kss-1:2*kse):: &
3496 vgsums !
3498 real,dimension(its-5:ite+5,jts-5:jte+5,kts:kte):: &
3499 dvol & ! grid box volume
3500 ,rdvol ! 1./grid box volume
3502 logical, save :: first=.true.
3503 real, save :: sumrrw_first, summass_first
3504 real :: summass
3505 double precision, save :: gsmp_first
3506 !+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
3508 ! steep=1.-0.040*idtad
3509 steep=1.
3511 DO K=KTS,KTE
3512 deta1_pdtop(k)=deta1(k)*pdtop
3513 enddo
3514 !
3515 DO K=KTS,KTE
3516 DO J=MYJS2,MYJE2
3517 DO I=MYIS1,MYIE1
3518 dvolp=(deta2(k)*pd(i,j)+deta1_pdtop(k))*dx(i,j)*dy
3519 rdvol(i,j,k)=hbm2(i,j)/dvolp
3520 dvol (i,j,k)=hbm2(i,j)*dvolp
3521 enddo
3522 enddo
3523 enddo
3525 ! go to 109
3526 !---monotonization------------------------------------------------------
3527 do ks=kss,kse ! loop by species
3528 !-----------------------------------------------------------------------
3529 DO K=KTS,KTE
3530 xsms(2*ks-1,k)=0.
3531 xsms(2*ks ,k)=0.
3532 gsms(2*ks-1,k)=0.
3533 gsms(2*ks ,k)=0.
3534 !
3535 DO J=MYJS2,MYJE2
3536 DO I=MYIS1,MYIE1
3537 !
3538 s1p=(s1(i,j,k,ks)+tcs(i,j,k,ks))**2
3539 tcs(i,j,k,ks)=s1p-s(i,j,k,ks)
3540 !
3541 sminh=min(s(i ,j-2,k,ks) &
3542 ,s(i+ihw(j),j-1,k,ks) &
3543 ,s(i+ihe(j),j-1,k,ks) &
3544 ,s(i-1 ,j ,k,ks) &
3545 ,s(i ,j ,k,ks) &
3546 ,s(i+1 ,j ,k,ks) &
3547 ,s(i+ihw(j),j+1,k,ks) &
3548 ,s(i+ihe(j),j+1,k,ks) &
3549 ,s(i ,j+2,k,ks))
3550 smaxh=max(s(i ,j-2,k,ks) &
3551 ,s(i+ihw(j),j-1,k,ks) &
3552 ,s(i+ihe(j),j-1,k,ks) &
3553 ,s(i-1 ,j ,k,ks) &
3554 ,s(i ,j ,k,ks) &
3555 ,s(i+1 ,j ,k,ks) &
3556 ,s(i+ihw(j),j+1,k,ks) &
3557 ,s(i+ihe(j),j+1,k,ks) &
3558 ,s(i ,j+2,k,ks))
3559 !
3560 if(k.gt.kts.and.k.lt.kte) then
3561 sminv=min(s(i,j,k-1,ks),s(i,j,k ,ks),s(i,j,k+1,ks))
3562 smaxv=max(s(i,j,k-1,ks),s(i,j,k ,ks),s(i,j,k+1,ks))
3563 elseif(k.eq.kts) then
3564 sminv=min(s(i,j,k ,ks),s(i,j,k+1,ks))
3565 smaxv=max(s(i,j,k ,ks),s(i,j,k+1,ks))
3566 elseif(k.eq.kte) then
3567 sminv=min(s(i,j,k-1,ks),s(i,j,k ,ks))
3568 smaxv=max(s(i,j,k-1,ks),s(i,j,k ,ks))
3569 endif
3570 !
3571 smin=min(sminh,sminv)
3572 smax=max(smaxh,smaxv)
3573 !
3574 sn=s1p
3575 if(sn.gt.steep*smax) sn=smax
3576 if(sn.lt. smin) sn=smin
3577 !
3578 dsks=(sn-s1p)*dvol(i,j,k)
3579 s1(i,j,k,ks)=dsks
3580 if(dsks.gt.0.) then
3581 xsms(2*ks-1,k)=xsms(2*ks-1,k)+dsks
3582 else
3583 xsms(2*ks ,k)=xsms(2*ks ,k)+dsks
3584 endif
3585 !
3586 enddo
3587 enddo
3588 enddo
3589 !-----------------------------------------------------------------------
3590 enddo ! end of the loop by species
3591 !-----------------------------------------------------------------------
3592 !-----------------------------------------------------------------------
3593 !*** GLOBAL REDUCTION
3594 !-----------------------------------------------------------------------
3595 !
3596 # ifdef DM_PARALLEL
3597 CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
3598 lngth=(2*kse-2*kss+2)*(KTE-KTS+1)
3599 CALL MPI_ALLREDUCE(XSMS,GSMS,lngth &
3600 & ,MPI_DOUBLE_PRECISION,MPI_SUM &
3601 & ,MPI_COMM_COMP,IRECV)
3602 # else
3603 DO K=KTS,KTE
3604 do ks=kss,kse
3605 gsms(2*ks-1,k)=xsms(2*ks-1,k)
3606 gsms(2*ks ,k)=xsms(2*ks ,k)
3607 enddo
3608 enddo
3609 # endif
3610 !
3611 DO K=KTS,KTE
3612 do ks=kss,kse
3613 gsms_single(2*ks-1,k)=gsms(2*ks-1,k)
3614 gsms_single(2*ks ,k)=gsms(2*ks ,k)
3615 enddo
3616 enddo
3617 !
3618 !-----------------------------------------------------------------------
3619 !*** END OF GLOBAL REDUCTION
3620 !-----------------------------------------------------------------------
3621 !
3622 !dusan!---forced conservation after monotonization----------------------------
3623 !dusan do ks=kss,kse
3624 !dusan DO K=KTS,KTE
3625 !dusan sumps=gsms_single(2*ks-1,k)
3626 !dusan sumns=gsms_single(2*ks ,k)
3627 !dusan!
3628 !dusan if(sumps*(-sumns).gt.1.) then
3629 !dusan sfacs=-sumns/sumps
3630 !dusan rfacs=1./sfacs
3631 !dusan else
3632 !dusan sfacs=0.
3633 !dusan rfacs=0.
3634 !dusan endif
3635 !dusan!
3636 !dusan DO J=MYJS2,MYJE2
3637 !dusan DO I=MYIS1,MYIE1
3638 !dusan dsks=s1(i,j,k,ks)*rdvol(i,j,k)
3639 !dusan if(sfacs.gt.1.) then
3640 !dusan if(dsks.lt.0.) dsks=dsks*rfacs
3641 !dusan!zjtest if(dsks.lt.0.) dsks=dsks
3642 !dusan else
3643 !dusan if(dsks.ge.0.) dsks=dsks*sfacs
3644 !dusan endif
3645 !dusan tcs(i,j,k,ks)=tcs(i,j,k,ks)+dsks
3646 !dusan enddo
3647 !dusan enddo
3648 !dusan enddo
3649 !dusan!-----------------------------------------------------------------------
3650 !dusan enddo ! end of the loop by species
3652 !---forced conservation after monotonization----------------------------
3653 do ks=kss,kse
3654 vgsums(2*ks-1)=0.
3655 vgsums(2*ks )=0.
3656 DO K=KTS,KTE
3657 vgsums(2*ks-1)=gsms(2*ks-1,k)+vgsums(2*ks-1)
3658 vgsums(2*ks )=gsms(2*ks ,k)+vgsums(2*ks )
3659 enddo
3660 enddo
3662 do ks=kss,kse
3663 sumps=vgsums(2*ks-1)
3664 sumns=vgsums(2*ks )
3665 !
3666 if(sumps*(-sumns).gt.1.) then
3667 sfacs=-sumns/sumps
3668 rfacs=1./sfacs
3669 else
3670 sfacs=0.
3671 rfacs=0.
3672 endif
3673 !
3674 DO K=KTS,KTE
3675 DO J=MYJS2,MYJE2
3676 DO I=MYIS1,MYIE1
3677 dsks=s1(i,j,k,ks)*rdvol(i,j,k)
3678 if(sfacs.lt.1.) then
3679 if(dsks.gt.0.) dsks=dsks*sfacs
3680 endif
3681 tcs(i,j,k,ks)=tcs(i,j,k,ks)+dsks
3682 enddo
3683 enddo
3684 enddo
3685 !-----------------------------------------------------------------------
3686 enddo ! end of the loop by species
3687 109 continue
3688 !-----------------------------------------------------------------------
3689 !zjwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwwww
3690 !-----------------------------------------------------------------------
3691 do ks=kss,kse ! loop by species
3692 DO K=KTS,KTE
3693 xsms(2*ks-1,k)=0.
3694 xsms(2*ks ,k)=0.
3695 gsms(2*ks-1,k)=0.
3696 gsms(2*ks ,k)=0.
3697 !
3698 DO J=MYJS2,MYJE2
3699 DO I=MYIS1,MYIE1
3700 xsms(2*ks-1,k)=xsms(2*ks-1,k)+s (i,j,k,ks)*dvol(i,j,k)
3701 xsms(2*ks ,k)=xsms(2*ks ,k)+tcs(i,j,k,ks)*dvol(i,j,k)
3702 enddo
3703 enddo
3704 enddo
3705 enddo
3706 !
3707 xsmp=0.
3708 DO J=MYJS2,MYJE2
3709 DO I=MYIS1,MYIE1
3710 xsmp=pd(i,j)*dx(i,j)*dy*hbm2(i,j)+xsmp
3711 enddo
3712 enddo
3713 !
3714 !-----------------------------------------------------------------------
3715 !*** GLOBAL REDUCTION
3716 !-----------------------------------------------------------------------
3717 !
3718 # ifdef DM_PARALLEL
3719 CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
3720 lngth=1
3721 CALL MPI_ALLREDUCE(xsmp,gsmp,lngth &
3722 & ,MPI_DOUBLE_PRECISION,MPI_SUM &
3723 & ,MPI_COMM_COMP,IRECV)
3724 # else
3725 gsmp=xsmp
3726 # endif
3727 !
3728 !-----------------------------------------------------------------------
3729 !*** END OF GLOBAL REDUCTION
3730 !-----------------------------------------------------------------------
3731 !
3732 !
3733 !-----------------------------------------------------------------------
3734 !*** GLOBAL REDUCTION
3735 !-----------------------------------------------------------------------
3736 !
3737 # ifdef DM_PARALLEL
3738 CALL WRF_GET_DM_COMMUNICATOR(MPI_COMM_COMP)
3739 lngth=(2*kse-2*kss+2)*(KTE-KTS+1)
3740 CALL MPI_ALLREDUCE(XSMS,GSMS,lngth &
3741 & ,MPI_DOUBLE_PRECISION,MPI_SUM &
3742 & ,MPI_COMM_COMP,IRECV)
3743 # else
3744 DO K=KTS,KTE
3745 do ks=kss,kse
3746 gsms(2*ks-1,k)=xsms(2*ks-1,k)
3747 gsms(2*ks ,k)=xsms(2*ks ,k)
3748 enddo
3749 enddo
3750 # endif
3751 !
3752 DO K=KTS,KTE
3753 do ks=kss,kse
3754 gsms_single(2*ks-1,k)=gsms(2*ks-1,k)
3755 gsms_single(2*ks ,k)=gsms(2*ks ,k)
3756 enddo
3757 enddo
3758 !
3759 !-----------------------------------------------------------------------
3760 !*** END OF GLOBAL REDUCTION
3761 !-----------------------------------------------------------------------
3762 !
3764 do ks=kss,kse
3765 vgsms(2*ks-1)=0.
3766 vgsms(2*ks )=0.
3767 enddo
3768 !
3769 do ks=kss,kse
3770 DO K=KTS,KTE
3771 vgsms(2*ks-1)=gsms(2*ks-1,k)+vgsms(2*ks-1)
3772 vgsms(2*ks )=gsms(2*ks ,k)+vgsms(2*ks )
3773 enddo
3774 enddo
3775 !
3776 sumdrrw=vgsms(6)+sumdrrw
3777 !
3778 summass=0.0
3779 DO K=KTS,KTE
3780 DO J=MYJS2,MYJE2
3781 DO I=MYIS1,MYIE1
3782 summass=summass+dvol(i,j,k)
3783 ENDDO
3784 ENDDO
3785 ENDDO
3786 !
3787 if (first) then
3788 sumrrw_first = vgsms(5)
3789 gsmp_first = gsmp
3790 summass_first = summass
3791 first=.false.
3792 end if
3793 !
3794 ! write(0,1000) ntsd,hours, & ! 4,5
3795 ! (vgsms(ks),ks=2*kss-1,2*kse), & ! 6-13
3796 ! gsmp,gsmp_first,gsmp/gsmp_first, & ! 14-16
3797 ! sumdrrw, & ! 17
3798 ! vgsms(5)/sumrrw_first,sumrrw_first, & ! 18-19
3799 ! summass,summass_first,summass/summass_first ! 20-22
3800 1000 format('global vol sums ',i6,f8.3,30d13.5)
3801 !zjmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmmm
3802 endsubroutine mono
3803 !-----------------------------------------------------------------------
3804 !-----------------------------------------------------------------------
3805 !
3806 END MODULE MODULE_ADVECTION
3807 !
3808 !-----------------------------------------------------------------------