r4627 | gill | 2010-12-29 16:29:58 -0700 (Wed, 29 Dec 2010) | 5 lines
[wrffire.git] / wrfv2_fire / dyn_em / module_big_step_utilities_em.F
blob3585c99c0a941c3db604d987fc502e0c49ec0ac5
1 !wrf:MODEL_LAYER:DYNAMICS
4 #if (RWORDSIZE == 4)
5 #   define VPOWX vspowx
6 #   define VPOW  vspow
7 #else
8 #   define VPOWX vpowx
9 #   define VPOW  vpow
10 #endif
13 MODULE module_big_step_utilities_em
15    USE module_domain, ONLY : domain
16    USE module_model_constants
17    USE module_state_description
18    USE module_configure, ONLY : grid_config_rec_type
19    USE module_wrf_error
21 CONTAINS
23 !-------------------------------------------------------------------------------
25 SUBROUTINE calc_mu_uv ( config_flags,                 &
26                         mu, mub, muu, muv,            &
27                         ids, ide, jds, jde, kds, kde, &
28                         ims, ime, jms, jme, kms, kme, &
29                         its, ite, jts, jte, kts, kte )
31    IMPLICIT NONE
32    
33    ! Input data
35    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
37    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
38                                        ims, ime, jms, jme, kms, kme, &
39                                        its, ite, jts, jte, kts, kte 
41    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
42    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub
44    !  local stuff
46    INTEGER :: i, j, itf, jtf, im, jm
48 !<DESCRIPTION>
50 !  calc_mu_uv calculates the full column dry-air mass at the staggered
51 !  horizontal velocity points (u,v) and places the results in muu and muv.
52 !  This routine uses the reference state (mub) and perturbation state (mu)
54 !</DESCRIPTION>
57       itf=ite
58       jtf=MIN(jte,jde-1)
60       IF      ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
61          DO j=jts,jtf
62          DO i=its,itf
63             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
64          ENDDO
65          ENDDO
66       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
67          DO j=jts,jtf
68          DO i=its+1,itf
69             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
70          ENDDO
71          ENDDO
72          i=its
73          im = its
74          if(config_flags%periodic_x) im = its-1
75          DO j=jts,jtf
76 !            muu(i,j) =      mu(i,j)          +mub(i,j)
77 !  fix for periodic b.c., 13 march 2004, wcs
78             muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
79          ENDDO
80       ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
81          DO j=jts,jtf
82          DO i=its,itf-1
83             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
84          ENDDO
85          ENDDO
86          i=ite
87          im = ite-1
88          if(config_flags%periodic_x) im = ite
89          DO j=jts,jtf
90 !            muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
91 !  fix for periodic b.c., 13 march 2004, wcs
92             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
93          ENDDO
94       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
95          DO j=jts,jtf
96          DO i=its+1,itf-1
97             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
98          ENDDO
99          ENDDO
100          i=its
101          im = its
102          if(config_flags%periodic_x) im = its-1
103          DO j=jts,jtf
104 !            muu(i,j) =      mu(i,j)          +mub(i,j)
105 !  fix for periodic b.c., 13 march 2004, wcs
106             muu(i,j) = 0.5*(mu(i,j)+mu(im,j)+mub(i,j)+mub(im,j))
107          ENDDO
108          i=ite
109          im = ite-1
110          if(config_flags%periodic_x) im = ite
111          DO j=jts,jtf
112 !            muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
113 !  fix for periodic b.c., 13 march 2004, wcs
114             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j)+mub(i-1,j)+mub(im,j))
115          ENDDO
116       END IF
118       itf=MIN(ite,ide-1)
119       jtf=jte
121       IF      ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
122          DO j=jts,jtf
123          DO i=its,itf
124              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
125          ENDDO
126          ENDDO
127       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
128          DO j=jts+1,jtf
129          DO i=its,itf
130              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
131          ENDDO
132          ENDDO
133          j=jts
134          jm = jts
135          if(config_flags%periodic_y) jm = jts-1
136          DO i=its,itf
137 !             muv(i,j) =      mu(i,j)          +mub(i,j)
138 !  fix for periodic b.c., 13 march 2004, wcs
139              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
140          ENDDO
141       ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
142          DO j=jts,jtf-1
143          DO i=its,itf
144              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
145          ENDDO
146          ENDDO
147          j=jte
148          jm = jte-1
149          if(config_flags%periodic_y) jm = jte
150          DO i=its,itf
151              muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
152 !  fix for periodic b.c., 13 march 2004, wcs
153              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
154          ENDDO
155       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
156          DO j=jts+1,jtf-1
157          DO i=its,itf
158              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
159          ENDDO
160          ENDDO
161          j=jts
162          jm = jts
163          if(config_flags%periodic_y) jm = jts-1
164          DO i=its,itf
165 !             muv(i,j) =      mu(i,j)          +mub(i,j)
166 !  fix for periodic b.c., 13 march 2004, wcs
167              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm)+mub(i,j)+mub(i,jm))
168          ENDDO
169          j=jte
170          jm = jte-1
171          if(config_flags%periodic_y) jm = jte
172          DO i=its,itf
173 !             muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
174 !  fix for periodic b.c., 13 march 2004, wcs
175              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm)+mub(i,j-1)+mub(i,jm))
176          ENDDO
177       END IF
179 END SUBROUTINE calc_mu_uv
181 !-------------------------------------------------------------------------------
183 SUBROUTINE calc_mu_uv_1 ( config_flags,                 &
184                           mu, muu, muv,                 &
185                           ids, ide, jds, jde, kds, kde, &
186                           ims, ime, jms, jme, kms, kme, &
187                           its, ite, jts, jte, kts, kte )
189    IMPLICIT NONE
190    
191    ! Input data
193    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
195    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
196                                        ims, ime, jms, jme, kms, kme, &
197                                        its, ite, jts, jte, kts, kte 
199    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
200    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
202    !  local stuff
204    INTEGER :: i, j, itf, jtf, im, jm
206 !<DESCRIPTION>
208 !  calc_mu_uv calculates the full column dry-air mass at the staggered
209 !  horizontal velocity points (u,v) and places the results in muu and muv.
210 !  This routine uses the full state (mu)
212 !</DESCRIPTION>
213    
214       itf=ite
215       jtf=MIN(jte,jde-1)
217       IF      ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
218          DO j=jts,jtf
219          DO i=its,itf
220             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
221          ENDDO
222          ENDDO
223       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
224          DO j=jts,jtf
225          DO i=its+1,itf
226             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
227          ENDDO
228          ENDDO
229          i=its
230          im = its
231          if(config_flags%periodic_x) im = its-1
232          DO j=jts,jtf
233             muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
234          ENDDO
235       ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
236          DO j=jts,jtf
237          DO i=its,itf-1
238             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
239          ENDDO
240          ENDDO
241          i=ite
242          im = ite-1
243          if(config_flags%periodic_x) im = ite
244          DO j=jts,jtf
245             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
246          ENDDO
247       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
248          DO j=jts,jtf
249          DO i=its+1,itf-1
250             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j))
251          ENDDO
252          ENDDO
253          i=its
254          im = its
255          if(config_flags%periodic_x) im = its-1
256          DO j=jts,jtf
257             muu(i,j) = 0.5*(mu(i,j)+mu(im,j))
258          ENDDO
259          i=ite
260          im = ite-1
261          if(config_flags%periodic_x) im = ite
262          DO j=jts,jtf
263             muu(i,j) = 0.5*(mu(i-1,j)+mu(im,j))
264          ENDDO
265       END IF
267       itf=MIN(ite,ide-1)
268       jtf=jte
270       IF      ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
271          DO j=jts,jtf
272          DO i=its,itf
273              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
274          ENDDO
275          ENDDO
276       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
277          DO j=jts+1,jtf
278          DO i=its,itf
279              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
280          ENDDO
281          ENDDO
282          j=jts
283          jm = jts
284          if(config_flags%periodic_y) jm = jts-1
285          DO i=its,itf
286              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
287          ENDDO
288       ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
289          DO j=jts,jtf-1
290          DO i=its,itf
291              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
292          ENDDO
293          ENDDO
294          j=jte
295          jm = jte-1
296          if(config_flags%periodic_y) jm = jte
297          DO i=its,itf
298              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
299          ENDDO
300       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
301          DO j=jts+1,jtf-1
302          DO i=its,itf
303              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1))
304          ENDDO
305          ENDDO
306          j=jts
307          jm = jts
308          if(config_flags%periodic_y) jm = jts-1
309          DO i=its,itf
310              muv(i,j) = 0.5*(mu(i,j)+mu(i,jm))
311          ENDDO
312          j=jte
313          jm = jte-1
314          if(config_flags%periodic_y) jm = jte
315          DO i=its,itf
316              muv(i,j) = 0.5*(mu(i,j-1)+mu(i,jm))
317          ENDDO
318       END IF
320 END SUBROUTINE calc_mu_uv_1
322 !-------------------------------------------------------------------------------
324 ! Map scale factor comments for this routine:
325 ! Locally not changed, but sent the correct map scale factors
326 ! from module_em (msfuy, msfvx, msfty)
328 SUBROUTINE couple_momentum ( muu, ru, u, msfu,              &
329                              muv, rv, v, msfv, msfv_inv,    &
330                              mut, rw, w, msft,              &
331                              ids, ide, jds, jde, kds, kde,  &
332                              ims, ime, jms, jme, kms, kme,  &
333                              its, ite, jts, jte, kts, kte  )
335    IMPLICIT NONE
337    ! Input data
339    INTEGER ,             INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
340                                           ims, ime, jms, jme, kms, kme, &
341                                           its, ite, jts, jte, kts, kte
343    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: ru, rv, rw
345    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: muu, muv, mut
346    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: msfu, msfv, msft, msfv_inv
347    
348    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: u, v, w
349    
350    ! Local data
351    
352    INTEGER :: i, j, k, itf, jtf, ktf
353    
354 !<DESCRIPTION>
356 ! couple_momentum couples the velocities to the full column mass and
357 ! the map factors.
359 !</DESCRIPTION>
361    ktf=MIN(kte,kde-1)
362    
363       itf=ite
364       jtf=MIN(jte,jde-1)
366       DO j=jts,jtf
367       DO k=kts,ktf
368       DO i=its,itf
369          ru(i,k,j)=u(i,k,j)*muu(i,j)/msfu(i,j)
370       ENDDO
371       ENDDO
372       ENDDO
374       itf=MIN(ite,ide-1)
375       jtf=jte
377       DO j=jts,jtf
378       DO k=kts,ktf
379       DO i=its,itf
380            rv(i,k,j)=v(i,k,j)*muv(i,j)*msfv_inv(i,j)
381 !           rv(i,k,j)=v(i,k,j)*muv(i,j)/msfv(i,j)
382       ENDDO
383       ENDDO
384       ENDDO
386       itf=MIN(ite,ide-1)
387       jtf=MIN(jte,jde-1)
389       DO j=jts,jtf
390       DO k=kts,kte
391       DO i=its,itf
392          rw(i,k,j)=w(i,k,j)*mut(i,j)/msft(i,j)
393       ENDDO
394       ENDDO
395       ENDDO
397 END SUBROUTINE couple_momentum
399 !-------------------------------------------------------------------
401 SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv,            &
402                                   ids, ide, jds, jde, kds, kde, &
403                                   ims, ime, jms, jme, kms, kme, &
404                                   its, ite, jts, jte, kts, kte )
406    IMPLICIT NONE
407    
408    ! Input data
410    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
411                                        ims, ime, jms, jme, kms, kme, &
412                                        its, ite, jts, jte, kts, kte 
414    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
415    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub
417    !  local stuff
419    INTEGER :: i, j, itf, jtf
421 !<DESCRIPTION>
423 ! calc_mu_staggered calculates the full dry air mass at the staggered
424 ! velocity points (u,v).
426 !</DESCRIPTION>
427    
428       itf=ite
429       jtf=MIN(jte,jde-1)
431       IF      ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
432          DO j=jts,jtf
433          DO i=its,itf
434             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
435          ENDDO
436          ENDDO
437       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
438          DO j=jts,jtf
439          DO i=its+1,itf
440             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
441          ENDDO
442          ENDDO
443          i=its
444          DO j=jts,jtf
445             muu(i,j) =      mu(i,j)          +mub(i,j)
446          ENDDO
447       ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
448          DO j=jts,jtf
449          DO i=its,itf-1
450             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
451          ENDDO
452          ENDDO
453          i=ite
454          DO j=jts,jtf
455             muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
456          ENDDO
457       ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
458          DO j=jts,jtf
459          DO i=its+1,itf-1
460             muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j))
461          ENDDO
462          ENDDO
463          i=its
464          DO j=jts,jtf
465             muu(i,j) =      mu(i,j)          +mub(i,j)
466          ENDDO
467          i=ite
468          DO j=jts,jtf
469             muu(i,j) =      mu(i-1,j)        +mub(i-1,j)
470          ENDDO
471       END IF
473       itf=MIN(ite,ide-1)
474       jtf=jte
476       IF      ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
477          DO j=jts,jtf
478          DO i=its,itf
479              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
480          ENDDO
481          ENDDO
482       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
483          DO j=jts+1,jtf
484          DO i=its,itf
485              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
486          ENDDO
487          ENDDO
488          j=jts
489          DO i=its,itf
490              muv(i,j) =      mu(i,j)          +mub(i,j)
491          ENDDO
492       ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
493          DO j=jts,jtf-1
494          DO i=its,itf
495              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
496          ENDDO
497          ENDDO
498          j=jte
499          DO i=its,itf
500              muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
501          ENDDO
502       ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
503          DO j=jts+1,jtf-1
504          DO i=its,itf
505              muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1))
506          ENDDO
507          ENDDO
508          j=jts
509          DO i=its,itf
510              muv(i,j) =      mu(i,j)          +mub(i,j)
511          ENDDO
512          j=jte
513          DO i=its,itf
514              muv(i,j) =      mu(i,j-1)        +mub(i,j-1)
515          ENDDO
516       END IF
518 END SUBROUTINE calc_mu_staggered
520 !-------------------------------------------------------------------------------
522 SUBROUTINE couple ( mu, mub, rfield, field, name, &
523                     msf,                          &
524                     ids, ide, jds, jde, kds, kde, &
525                     ims, ime, jms, jme, kms, kme, &
526                     its, ite, jts, jte, kts, kte )
528    IMPLICIT NONE
530    ! Input data
532    INTEGER ,             INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
533                                           ims, ime, jms, jme, kms, kme, &
534                                           its, ite, jts, jte, kts, kte
536    CHARACTER(LEN=1) ,     INTENT(IN   ) :: name
538    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: rfield
540    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub, msf
541    
542    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field
543    
544    ! Local data
545    
546    INTEGER :: i, j, k, itf, jtf, ktf
547    REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
549 !<DESCRIPTION>
551 ! subroutine couple couples the input variable with the dry-air 
552 ! column mass (mu).  
554 !</DESCRIPTION>
556    
557    ktf=MIN(kte,kde-1)
558    
559    IF (name .EQ. 'u')THEN
561       CALL calc_mu_staggered ( mu, mub, muu, muv,            &
562                                   ids, ide, jds, jde, kds, kde, &
563                                   ims, ime, jms, jme, kms, kme, &
564                                   its, ite, jts, jte, kts, kte )
566       itf=ite
567       jtf=MIN(jte,jde-1)
569       DO j=jts,jtf
570       DO k=kts,ktf
571       DO i=its,itf
572          rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j)
573       ENDDO
574       ENDDO
575       ENDDO
577    ELSE IF (name .EQ. 'v')THEN
579       CALL calc_mu_staggered ( mu, mub, muu, muv,            &
580                                ids, ide, jds, jde, kds, kde, &
581                                ims, ime, jms, jme, kms, kme, &
582                                its, ite, jts, jte, kts, kte )
584       itf=ite
585       itf=MIN(ite,ide-1)
586       jtf=jte
588       DO j=jts,jtf
589       DO k=kts,ktf
590       DO i=its,itf
591            rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j)
592       ENDDO
593       ENDDO
594       ENDDO
596    ELSE IF (name .EQ. 'w')THEN
597       itf=MIN(ite,ide-1)
598       jtf=MIN(jte,jde-1)
599       DO j=jts,jtf
600       DO k=kts,kte
601       DO i=its,itf
602          rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j)
603       ENDDO
604       ENDDO
605       ENDDO
607    ELSE IF (name .EQ. 'h')THEN
608       itf=MIN(ite,ide-1)
609       jtf=MIN(jte,jde-1)
610       DO j=jts,jtf
611       DO k=kts,kte
612       DO i=its,itf
613          rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
614       ENDDO
615       ENDDO
616       ENDDO
618    ELSE 
619       itf=MIN(ite,ide-1)
620       jtf=MIN(jte,jde-1)
621       DO j=jts,jtf
622       DO k=kts,ktf
623       DO i=its,itf
624          rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))
625       ENDDO
626       ENDDO
627       ENDDO
628    
629    ENDIF
631 END SUBROUTINE couple
634 !-------------------------------------------------------------------------------
636 SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww,              &
637                         rdx, rdy, msftx, msfty,          &
638                         msfux, msfuy, msfvx, msfvx_inv,  &
639                         msfvy, dnw,                      &
640                         ids, ide, jds, jde, kds, kde,    &
641                         ims, ime, jms, jme, kms, kme,    &
642                         its, ite, jts, jte, kts, kte    )
644    IMPLICIT NONE
646    ! Input data
649    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
650                                  ims, ime, jms, jme, kms, kme, &
651                                  its, ite, jts, jte, kts, kte
653    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: u, v
654    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mup, mub, &
655                                                             msftx, msfty, &
656                                                             msfux, msfuy, &
657                                                             msfvx, msfvy, &
658                                                             msfvx_inv
659    REAL , DIMENSION( kms:kme ) , INTENT(IN   ) :: dnw
660    
661    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT  ) :: ww
662    REAL , INTENT(IN   )  :: rdx, rdy
663    
664    ! Local data
665    
666    INTEGER :: i, j, k, itf, jtf, ktf
667    REAL , DIMENSION( its:ite ) :: dmdt
668    REAL , DIMENSION( its:ite, kts:kte ) :: divv
669    REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
671 !<DESCRIPTION>
673 !  calc_ww calculates omega using the velocities (u,v) and the dry-air 
674 !  column mass (mup+mub).
675 !  The algorithm integrates the continuity equation through the column
676 !  followed by a diagnosis of omega.
678 !</DESCRIPTION>
680 !<DESCRIPTION>
682 !  calc_ww_cp calculates omega using the velocities (u,v) and the 
683 !  column mass mu.
685 !</DESCRIPTION>
687     jtf=MIN(jte,jde-1)
688     ktf=MIN(kte,kde-1)  
689     itf=MIN(ite,ide-1)
691 !  mu coupled with the appropriate map factor
693       DO j=jts,jtf
694       DO i=its,min(ite+1,ide)
695         ! u is always coupled with my
696         muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfuy(i,j)
697       ENDDO
698       ENDDO
700       DO j=jts,min(jte+1,jde)
701       DO i=its,itf
702        ! v is always coupled with mx
703 !        muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfvx(i,j)
704         muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))*msfvx_inv(i,j)
705       ENDDO
706       ENDDO
708       DO j=jts,jtf
710         DO i=its,ite
711           dmdt(i) = 0.
712           ww(i,1,j) = 0.
713           ww(i,kte,j) = 0.
714         ENDDO
716 !       Comments on the modifications for map scale factors
717 !       ADT eqn 47 / my (putting rho -> 'mu') is:
718 !       (1/my) partial d mu/dt = -mx partial d/dx(mu u/my) 
719 !                                -mx partial d/dy(mu v/mx)
720 !                                -partial d/dz(mu w/my)
722 !       Using nu instead of z the last term becomes:
723 !                                -partial d/dnu(mu (dnu/dt)/my)
725 !       Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
726 !       and bottom, the last term becomes = 0
728 !       Integral|bot->top[(1/my) partial d mu/dt]dnu = 
729 !       Integral|bot->top[-mx partial d/dx(mu u/my) 
730 !                         -mx partial d/dy(mu v/mx)]dnu
732 !       muu='mu'[on u]/my, muv='mu'[on v]/mx
733 !       (1/my) partial d mu/dt is independent of nu
734 !         => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
736 !         => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
737 !                                        partial d/dy(mu v/mx)]dnu
738 !         => dmdt = sum_bot->top[divv]
739 !       where
740 !         divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
742         DO k=kts,ktf
743         DO i=its,itf
745           divv(i,k) = msftx(i,j)*dnw(k)*( rdx*(muu(i+1,j)*u(i+1,k,j)-muu(i,j)*u(i,k,j))  &
746                                         +rdy*(muv(i,j+1)*v(i,k,j+1)-muv(i,j)*v(i,k,j))   )
748 !          dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j))  &
749 !                                       +rdy*(rv(i,k,j+1)-rv(i,k,j))   )
751           dmdt(i) = dmdt(i) + divv(i,k)
754         ENDDO
755         ENDDO
757 !       Further map scale factor notes:
758 !       Now integrate from bottom to top, level by level:
759 !       mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt 
760 !                           -mx partial d/dx(mu u/my)
761 !                           -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
762 !       ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
763 !                = ww [k] -dmdt * dnw[k] - divv[k]
765         DO k=2,ktf
766         DO i=its,itf
768 !           ww(i,k,j)=ww(i,k-1,j)                                       &
769 !                        - dnw(k-1)* ( dmdt(i)                          &
770 !                                     +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j))  &
771 !                                     +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
773            ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
775         ENDDO
776         ENDDO
777      ENDDO
780 END SUBROUTINE calc_ww_cp
783 !-------------------------------------------------------------------------------
785 SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
786                      ids, ide, jds, jde, kds, kde,  &
787                      ims, ime, jms, jme, kms, kme,  &
788                      its, ite, jts, jte, kts, kte  )
790    IMPLICIT NONE
791    
792    ! Input data
794    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
795                                        ims, ime, jms, jme, kms, kme, &
796                                        its, ite, jts, jte, kts, kte 
798    INTEGER ,          INTENT(IN   ) :: n_moist
799    
801    REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN   ) :: moist
802                                               
803    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: cqu, cqv, cqw
805    ! Local stuff
807    REAL :: qtot
808    
809    INTEGER :: i, j, k, itf, jtf, ktf, ispe
811 !<DESCRIPTION>
813 !  calc_cq calculates moist coefficients for the momentum equations.
815 !</DESCRIPTION>
817       itf=ite
818       jtf=MIN(jte,jde-1)
819       ktf=MIN(kte,kde-1)
821       IF(  n_moist >= PARAM_FIRST_SCALAR ) THEN
823         DO j=jts,jtf
824         DO k=kts,ktf
825         DO i=its,itf
826           qtot = 0.
827 !DEC$ loop count(3)
828           DO ispe=PARAM_FIRST_SCALAR,n_moist
829             qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
830           ENDDO
831 !           qtot = 0.5*( moist(i  ,k,j,1)+moist(i  ,k,j,2)+moist(i  ,k,j,3)+  &
832 !     &                  moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) )
833 !           cqu(i,k,j) = 1./(1.+qtot)
834            cqu(i,k,j) = 1./(1.+0.5*qtot)
835         ENDDO
836         ENDDO
837         ENDDO
839         itf=MIN(ite,ide-1)
840         jtf=jte
842         DO j=jts,jtf
843         DO k=kts,ktf
844         DO i=its,itf
845           qtot = 0.
846 !DEC$ loop count(3)
847           DO ispe=PARAM_FIRST_SCALAR,n_moist
848             qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
849           ENDDO
850 !           qtot = 0.5*( moist(i,k,j  ,1)+moist(i,k,j  ,2)+moist(i,k,j  ,3)+  &
851 !     &                  moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) )
852 !           cqv(i,k,j) = 1./(1.+qtot)
853            cqv(i,k,j) = 1./(1.+0.5*qtot)
854         ENDDO
855         ENDDO
856         ENDDO
858         itf=MIN(ite,ide-1)
859         jtf=MIN(jte,jde-1)
860         DO j=jts,jtf
861         DO k=kts+1,ktf
862         DO i=its,itf
863           qtot = 0.
864 !DEC$ loop count(3)
865           DO ispe=PARAM_FIRST_SCALAR,n_moist
866             qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
867           ENDDO
868 !           qtot = 0.5*( moist(i,k  ,j,1)+moist(i,k  ,j,2)+moist(i,k-1,j,3)+  &
869 !     &                  moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k  ,j,3) )
870 !           cqw(i,k,j) = qtot
871            cqw(i,k,j) = 0.5*qtot
872         ENDDO
873         ENDDO
874         ENDDO
876       ELSE
878         DO j=jts,jtf
879         DO k=kts,ktf
880         DO i=its,itf
881            cqu(i,k,j) = 1.
882         ENDDO
883         ENDDO
884         ENDDO
886         itf=MIN(ite,ide-1)
887         jtf=jte
889         DO j=jts,jtf
890         DO k=kts,ktf
891         DO i=its,itf
892            cqv(i,k,j) = 1.
893         ENDDO
894         ENDDO
895         ENDDO
897         itf=MIN(ite,ide-1)
898         jtf=MIN(jte,jde-1)
899         DO j=jts,jtf
900         DO k=kts+1,ktf
901         DO i=its,itf
902            cqw(i,k,j) = 0.
903         ENDDO
904         ENDDO
905         ENDDO
907       END IF
909 END SUBROUTINE calc_cq
911 !----------------------------------------------------------------------
913 SUBROUTINE calc_alt ( alt, al, alb,                  &
914                       ids, ide, jds, jde, kds, kde,  &
915                       ims, ime, jms, jme, kms, kme,  &
916                       its, ite, jts, jte, kts, kte  )
918    IMPLICIT NONE
919    
920    ! Input data
922    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
923                                        ims, ime, jms, jme, kms, kme, &
924                                        its, ite, jts, jte, kts, kte 
926    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alb, al
927    REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(  OUT) :: alt
929    ! Local stuff
931    INTEGER :: i, j, k, itf, jtf, ktf
933 !<DESCRIPTION>
935 ! calc_alt computes the full inverse density
937 !</DESCRIPTION>
939       itf=MIN(ite,ide-1)
940       jtf=MIN(jte,jde-1)
941       ktf=MIN(kte,kde-1)
943       DO j=jts,jtf
944       DO k=kts,ktf
945       DO i=its,itf
946         alt(i,k,j) = al(i,k,j)+alb(i,k,j)
947       ENDDO
948       ENDDO
949       ENDDO
952 END SUBROUTINE calc_alt
954 !----------------------------------------------------------------------
956 SUBROUTINE calc_p_rho_phi ( moist, n_moist,                &
957                             al, alb, mu, muts, ph, p, pb,  &
958                             t, p0, t0, znu, dnw, rdnw,     &
959                             rdn, non_hydrostatic,          &
960                             ids, ide, jds, jde, kds, kde,  &
961                             ims, ime, jms, jme, kms, kme,  &
962                             its, ite, jts, jte, kts, kte  )
964   IMPLICIT NONE
965    
966    ! Input data
968   LOGICAL ,          INTENT(IN   ) :: non_hydrostatic
970   INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
971                                       ims, ime, jms, jme, kms, kme, &
972                                       its, ite, jts, jte, kts, kte 
974   INTEGER ,          INTENT(IN   ) :: n_moist
976   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alb,  &
977                                                                    pb,   &
978                                                                    t
980   REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN   ) :: moist
982   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(  OUT) :: al, p
984   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph
986   REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN   ) :: mu, muts
988   REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: znu, dnw, rdnw, rdn
990   REAL,   INTENT(IN   ) :: t0, p0
992   ! Local stuff
994   INTEGER :: i, j, k, itf, jtf, ktf, ispe
995   REAL    :: qvf, qtot, qf1, qf2
996   REAL, DIMENSION( its:ite) :: temp,cpovcv_v
999 !<DESCRIPTION>
1001 ! For the nonhydrostatic option, calc_p_rho_phi calculates the
1002 ! diagnostic quantities pressure and (inverse) density from the
1003 ! prognostic variables using the equation of state.
1005 ! For the hydrostatic option, calc_p_rho_phi calculates the
1006 ! diagnostic quantities (inverse) density and geopotential from the
1007 ! prognostic variables using the equation of state and the hydrostatic 
1008 ! equation.
1010 !</DESCRIPTION>
1012   itf=MIN(ite,ide-1)
1013   jtf=MIN(jte,jde-1)
1014   ktf=MIN(kte,kde-1)
1016 #ifndef INTELMKL
1017   cpovcv_v = cpovcv
1018 #endif
1020   IF (non_hydrostatic) THEN
1022       IF (n_moist >= PARAM_FIRST_SCALAR ) THEN  
1024         DO j=jts,jtf
1025         DO k=kts,ktf
1026         DO i=its,itf
1027           qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1028           al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j)  &
1029                      +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1030           temp(i)=(r_d*(t0+t(i,k,j))*qvf)/                 &
1031                         (p0*(al(i,k,j)+alb(i,k,j)))
1032         ENDDO
1033 #ifdef INTELMKL
1034         CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1035 #else
1036 ! use vector version from libmassv or from compat lib in frame/libmassv.F
1037         CALL VPOW  ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1038 #endif
1039         DO i=its,itf
1040            p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1041         ENDDO
1042         ENDDO
1043         ENDDO
1045       ELSE
1047         DO j=jts,jtf
1048         DO k=kts,ktf
1049         DO i=its,itf
1050           al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j)  &
1051                      +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1052           p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/                     &
1053                         (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv  &
1054                            -pb(i,k,j)
1055         ENDDO
1056         ENDDO
1057         ENDDO
1059       END IF
1061    ELSE
1063 !  hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1066       IF (n_moist >= PARAM_FIRST_SCALAR ) THEN  
1068         DO j=jts,jtf
1070           k=ktf          ! top layer
1071           DO i=its,itf
1073             qtot = 0.
1074             DO ispe=PARAM_FIRST_SCALAR,n_moist
1075               qtot = qtot + moist(i,k,j,ispe)
1076             ENDDO
1077             qf2 = 1./(1.+qtot)
1078             qf1 = qtot*qf2
1080             p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1081             qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1082             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1083                 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1085           ENDDO
1087           DO k=ktf-1,kts,-1  ! remaining layers, integrate down
1088             DO i=its,itf
1090             qtot = 0.
1091             DO ispe=PARAM_FIRST_SCALAR,n_moist
1092               qtot = qtot + 0.5*(  moist(i,k  ,j,ispe) + moist(i,k+1,j,ispe) )
1093             ENDDO
1094             qf2 = 1./(1.+qtot)
1095             qf1 = qtot*qf2
1097             p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1098             qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1099             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1100                         (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1101             ENDDO
1102           ENDDO
1104           DO k=2,ktf+1  ! integrate hydrostatic equation for geopotential
1105             DO i=its,itf
1107 !              ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*(       &
1108 !                           (muts(i,j)+mu(i,j))*al(i,k-1,j)+    &
1109 !                            mu(i,j)*alb(i,k-1,j)  )
1110               ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*(           &
1111                            (muts(i,j))*al(i,k-1,j)+            &
1112                             mu(i,j)*alb(i,k-1,j)  )
1113                                                    
1115             ENDDO
1116           ENDDO
1118         ENDDO
1120       ELSE
1122         DO j=jts,jtf
1124           k=ktf          ! top layer
1125           DO i=its,itf
1127             qtot = 0.
1128             qf2 = 1./(1.+qtot)
1129             qf1 = qtot*qf2
1131             p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2
1132             qvf = 1.
1133             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1134                 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1136           ENDDO
1138           DO k=ktf-1,kts,-1  ! remaining layers, integrate down
1139             DO i=its,itf
1141             qtot = 0.
1142             qf2 = 1./(1.+qtot)
1143             qf1 = qtot*qf2
1145             p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1)
1146             qvf = 1.
1147             al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1148                         (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1149             ENDDO
1150           ENDDO
1152           DO k=2,ktf+1  ! integrate hydrostatic equation for geopotential
1153             DO i=its,itf
1155 !              ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*(       &
1156 !                           (muts(i,j)+mu(i,j))*al(i,k-1,j)+    &
1157 !                            mu(i,j)*alb(i,k-1,j)  )
1158               ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*(           &
1159                            (muts(i,j))*al(i,k-1,j)+            &
1160                             mu(i,j)*alb(i,k-1,j)  )
1161                                                    
1163             ENDDO
1164           ENDDO
1166         ENDDO
1168      END IF
1170    END IF
1172 END SUBROUTINE calc_p_rho_phi
1174 !----------------------------------------------------------------------
1176 SUBROUTINE calc_php ( php, ph, phb,                  &
1177                       ids, ide, jds, jde, kds, kde,  &
1178                       ims, ime, jms, jme, kms, kme,  &
1179                       its, ite, jts, jte, kts, kte  )
1181    IMPLICIT NONE
1182    
1183    ! Input data
1185    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1186                                        ims, ime, jms, jme, kms, kme, &
1187                                        its, ite, jts, jte, kts, kte 
1189    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) :: phb, ph
1190    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(  OUT) :: php
1192    ! Local stuff
1194    INTEGER :: i, j, k, itf, jtf, ktf
1196 !<DESCRIPTION>
1198 !  calc_php calculates the full geopotential from the reference state
1199 !  geopotential and the perturbation geopotential (phb_ph).
1201 !</DESCRIPTION>
1203       itf=MIN(ite,ide-1)
1204       jtf=MIN(jte,jde-1)
1205       ktf=MIN(kte,kde-1)
1207       DO j=jts,jtf
1208       DO k=kts,ktf
1209       DO i=its,itf
1210         php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1211       ENDDO
1212       ENDDO
1213       ENDDO
1215 END SUBROUTINE calc_php
1217 !-------------------------------------------------------------------------------
1219 SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt,  &
1220                        u, v, ht,                            &
1221                        cf1, cf2, cf3, rdx, rdy,             &
1222                        msftx, msfty,                        &
1223                        ids, ide, jds, jde, kds, kde,        &
1224                        ims, ime, jms, jme, kms, kme,        &
1225                        its, ite, jts, jte, kts, kte        )
1227    IMPLICIT NONE
1229    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1230                                        ims, ime, jms, jme, kms, kme, &
1231                                        its, ite, jts, jte, kts, kte 
1233    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   ph_tend, &
1234                                                                      ph_new,  &
1235                                                                      ph_old,  &
1236                                                                      u,       &
1237                                                                      v
1240    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(  OUT) :: w
1242    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mu, ht, msftx, msfty
1244    REAL, INTENT(IN   ) :: dt, cf1, cf2, cf3, rdx, rdy
1246    INTEGER :: i, j, k, itf, jtf
1248    itf=MIN(ite,ide-1)
1249    jtf=MIN(jte,jde-1)
1251 !<DESCRIPTION>
1253 ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1254 ! Used with the hydrostatic option.
1256 !</DESCRIPTION>
1258    DO j = jts, jtf
1260 !  lower b.c. on w
1262 !  Notes on map scale factors:
1263 !  Chain rule: if Z=Z(X,Y) [true at the surface] then
1264 !  dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1265 !  Using capitals to denote actual values
1266 !  In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1267 !    => w = dz/dt = mx u dz/dx + my v dz/dy
1268 !  [where dz/dx is just the surface height change between x
1269 !   gridpoints, and dz/dy is the change between y gridpoints]
1270 !  [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1271 !   nearest the surface]
1273 !  Previously msft multiplied by rdy and rdx terms.
1274 !  Now msfty multiplies rdy term, and msftx multiplies msftx term
1275      DO i = its, itf
1276          w(i,1,j)=  msfty(i,j)*.5*rdy*(                      &
1277                            (ht(i,j+1)-ht(i,j  ))             &
1278           *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1))    &
1279                           +(ht(i,j  )-ht(i,j-1))             &
1280           *(cf1*v(i,1,j  )+cf2*v(i,2,j  )+cf3*v(i,3,j  ))  ) &
1281                  +msftx(i,j)*.5*rdx*(                        &
1282                            (ht(i+1,j)-ht(i,j  ))             &
1283           *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j))    &
1284                           +(ht(i,j  )-ht(i-1,j))             &
1285           *(cf1*u(i  ,1,j)+cf2*u(i  ,2,j)+cf3*u(i  ,3,j))  )
1286      ENDDO
1288 !  use geopotential equation to diagnose w
1290 !  Further notes on map scale factors
1291 !  If ph_tend contains:  -mx partial d/dx(mu rho u/my) 
1292 !                        -mx partial d/dy(phi mu v/mx)
1293 !                        -partial d/dz(phi mu w/my)
1294 !  then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1295 !    => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
1297      DO k = 2, kte
1298      DO i = its, itf
1299        w(i,k,j) =  msfty(i,j)*(  (ph_new(i,k,j)-ph_old(i,k,j))/dt       &
1300                                - ph_tend(i,k,j)/mu(i,j)        )/g 
1302      ENDDO
1303      ENDDO
1305    ENDDO
1307 END SUBROUTINE diagnose_w
1309 !-------------------------------------------------------------------------------
1311 SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
1312                    ph, ph_old, phb, w,              &
1313                    mut, muu, muv,                   &
1314                    fnm, fnp,                        &
1315                    rdnw, cfn, cfn1, rdx, rdy,       &
1316                    msfux, msfuy, msfvx,             &
1317                    msfvx_inv, msfvy,                &
1318                    msftx, msfty,                    &
1319                    non_hydrostatic,                 &
1320                    config_flags,                    &
1321                    ids, ide, jds, jde, kds, kde,    &
1322                    ims, ime, jms, jme, kms, kme,    &
1323                    its, ite, jts, jte, kts, kte    )
1324    IMPLICIT NONE
1326    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
1328    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1329                                        ims, ime, jms, jme, kms, kme, &
1330                                        its, ite, jts, jte, kts, kte 
1332    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::        &
1333                                                                      u,   &
1334                                                                      v,   &
1335                                                                      ww,  &
1336                                                                      ph,  &
1337                                                                      ph_old, &
1338                                                                      phb, & 
1339                                                                     w
1341    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1343    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: muu, muv, mut,   &
1344                                                             msfux, msfuy, &
1345                                                             msfvx, msfvy, &
1346                                                             msftx, msfty, &
1347                                                             msfvx_inv
1349    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, fnm, fnp
1351    REAL,  INTENT(IN   ) :: cfn, cfn1, rdx, rdy
1353    LOGICAL,  INTENT(IN   )  ::  non_hydrostatic
1355    ! Local stuff
1357    INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1358    REAL    :: ur, ul, ub, vr, vl, vb
1359    REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1361    INTEGER :: advective_order
1363    LOGICAL :: specified
1365 !<DESCRIPTION>
1367 ! rhs_ph calculates the large-timestep tendency terms for the geopotential
1368 ! equation.  These terms include the advection and "gw".  The geopotential
1369 ! equation is cast in advective form, so we don't use the flux form advection
1370 ! algorithms here.
1372 !</DESCRIPTION>
1374    specified = .false.
1375    if(config_flags%specified .or. config_flags%nested) specified = .true.
1377    advective_order = config_flags%h_sca_adv_order 
1379    itf=MIN(ite,ide-1)
1380    jtf=MIN(jte,jde-1)
1381    ktf=MIN(kte,kde-1)
1383 !  Notes on map scale factors (WCS, 2 march 2008)
1384 !  phi equation is:   mu/my d/dt(phi) = -(1/my) mx mu u  d/dx(phi)
1385 !                                       -(1/my) my mu v d/dy(phi)
1386 !                                       - omega d/d_eta(phi)
1387 !                                               +mu g w/my
1389 !  A little further explanation...
1390 !  The tendency term we are computing here is for mu/my d/dt(phi).  It is advective form 
1391 !  but it is multiplied be mu/my.  It will be decoupled from (mu/my) when the implicit w-phi
1392 !  solution is computed in subourine advance_w.  The formulation dates from the early 
1393 !  days of the mass coordinate model when we were testing both a flux and an advective formulation
1394 !  for the geopotential equation and different forms of the vertical momentum equation and the 
1395 !  vertically implicit solver.
1397 ! advective form for the geopotential equation
1399    DO j = jts, jtf
1401      DO k = 2, kte
1402      DO i = its, itf
1403           wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1)               &
1404                         *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1405      ENDDO
1406      ENDDO
1408 !  RHS term 3 is: - omega partial d/dnu(phi)
1410      DO k = 2, kte-1
1411      DO i = its, itf
1412            ph_tend(i,k,j) = ph_tend(i,k,j)                           &
1413                              - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1414      ENDDO
1415      ENDDO
1417    ENDDO
1419    IF (non_hydrostatic) THEN  ! add in "gw" term.
1420    DO j = jts, jtf            ! in hydrostatic mode, "gw" will be diagnosed
1421                               ! after the timestep to give us "w"
1422      DO i = its, itf
1423         ph_tend(i,kde,j) = 0.
1424      ENDDO
1426      DO k = 2, kte
1427      DO i = its, itf
1428         ! phi equation RHS term 4
1429         ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
1430      ENDDO
1431      ENDDO
1433    ENDDO
1435    END IF
1437 !  Notes on map scale factors:
1438 !  RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1439 !                         -(1/my) my v mu partial d/dy(phi)
1441    IF (advective_order <= 2) THEN
1443 !  y (v) advection
1445    i_start = its
1446    j_start = jts
1447    itf=MIN(ite,ide-1)
1448    jtf=MIN(jte,jde-1)
1450    IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
1451    IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-2
1453    DO j = j_start, jtf
1455      DO k = 2, kte-1
1456      DO i = i_start, itf
1457         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
1458                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
1459                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
1460                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
1461                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
1462      ENDDO
1463      ENDDO
1465      k = kte
1466      DO i = i_start, itf
1467         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
1468                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
1469                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
1470                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
1471                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
1472      ENDDO
1474    ENDDO
1476 !  x (u) advection
1478    i_start = its
1479    j_start = jts
1480    itf=MIN(ite,ide-1)
1481    jtf=MIN(jte,jde-1)
1483    IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
1484    IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = itf-2
1486    DO j = j_start, jtf
1488      DO k = 2, kte-1
1489      DO i = i_start, itf
1490         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*         &
1491                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)*  &
1492                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))  &
1493                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j)*  &
1494                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))  )
1495      ENDDO
1496      ENDDO
1498      k = kte
1499      DO i = i_start, itf
1500         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*                        &
1501                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)*    &
1502                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))               &
1503                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(  i,j)*    &
1504                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))             )
1505      ENDDO
1507    ENDDO
1509    ELSE IF (advective_order <= 4) THEN
1511 !  y (v) advection
1513    i_start = its
1514    j_start = jts
1515    itf=MIN(ite,ide-1)
1516    jtf=MIN(jte,jde-1)
1518    IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
1519    IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jtf-3
1521    DO j = j_start, jtf
1523      DO k = 2, kte-1
1524      DO i = i_start, itf
1525         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*(                     &
1526                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)                &
1527                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  ))* (1./12.)*(   &
1528                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                    &
1529                       -(ph(i,k,j+2)-ph(i,k,j-2))                                    &
1530                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                  &
1531                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1534      ENDDO
1535      ENDDO
1537      k = kte
1538      DO i = i_start, itf
1539         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*(                                 &
1540                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)                &
1541                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  ))* (1./12.)*(   &
1542                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                               &
1543                       -(ph(i,k,j+2)-ph(i,k,j-2))                                               &
1544                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                             &
1545                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1547      ENDDO
1549    ENDDO
1551 !  pick up near boundary rows using 2nd order stencil for open and specified conditions
1553    IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 )  THEN
1555      j = jds+1
1556      DO k = 2, kte-1
1557      DO i = i_start, itf
1558         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
1559                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
1560                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
1561                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
1562                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
1563      ENDDO
1564      ENDDO
1566      k = kte
1567      DO i = i_start, itf
1568         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
1569                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
1570                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
1571                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
1572                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
1573      ENDDO
1575    END IF
1577    IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 )  THEN
1579      j = jde-2
1580      DO k = 2, kte-1
1581      DO i = i_start, itf
1582         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
1583                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
1584                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
1585                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
1586                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
1587      ENDDO
1588      ENDDO
1590      k = kte
1591      DO i = i_start, itf
1592         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
1593                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
1594                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
1595                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
1596                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
1597      ENDDO
1599    END IF
1601 !  x (u) advection
1603    i_start = its
1604    j_start = jts
1605    itf=MIN(ite,ide-1)
1606    jtf=MIN(jte,jde-1)
1608    IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
1609    IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-3
1611    DO j = j_start, jtf
1613      DO k = 2, kte-1
1614      DO i = i_start, itf
1615         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*(                    &
1616                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)               &
1617                   +muu(i,j  )*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j) )* (1./12.)*( &
1618                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                                   &
1619                       -(ph(i+2,k,j)-ph(i-2,k,j))                                   &
1620                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                                 &
1621                       -(phb(i+2,k,j)-phb(i-2,k,j))  )   )                
1622      ENDDO
1623      ENDDO
1625      k = kte
1626      DO i = i_start, itf
1627         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*(                                 &
1628                  ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)                &
1629                   +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(i  ,j) )* (1./12.)*(  &
1630                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                                               &
1631                       -(ph(i+2,k,j)-ph(i-2,k,j))                                               &
1632                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                                             &
1633                       -(phb(i+2,k,j)-phb(i-2,k,j))  )     )
1634      ENDDO
1636    ENDDO
1638 !  pick up near boundary rows using 2nd order stencil for open and specified conditions
1640    IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
1642      i = ids + 1
1644      DO j = j_start, jtf
1645      DO k = 2, kte-1
1646         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*         &
1647                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)*  &
1648                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))  &
1649                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j)*  &
1650                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))  )
1651      ENDDO
1652      ENDDO
1654      k = kte
1655      DO j = j_start, jtf
1656         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*                        &
1657                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)*    &
1658                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))               &
1659                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(  i,j)*    &
1660                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))             )
1661      ENDDO
1663    END IF
1665    IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
1667      i = ide-2
1668      DO j = j_start, jtf
1669      DO k = 2, kte-1
1670         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*         &
1671                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)*  &
1672                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))  &
1673                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j)*  &
1674                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))  )
1675      ENDDO
1676      ENDDO
1678      k = kte
1679      DO j = j_start, jtf
1680         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*                        &
1681                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)*    &
1682                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))               &
1683                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(  i,j)*    &
1684                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))             )
1685      ENDDO
1687    END IF
1689 !--------------------------
1691    ELSE IF (advective_order <= 6) THEN
1693 !!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
1694 !!!       the branches covering the other advective_order cases have not.  20090923. JM
1696 !  y (v) advection
1698    i_start = its
1699    j_start = jts
1700    itf=MIN(ite,ide-1)
1701    jtf=MIN(jte,jde-1)
1703    IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
1704    IF (config_flags%open_ye .or. specified ) jtf     = min(jtf,jde-4)
1706    DO j = j_start, jtf
1708      DO k = 2, kte-1
1709      DO i = i_start, itf
1710         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* (                    &
1711                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)                &
1712                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  ) )* (1./60.)*(  &
1713                    45.*(ph(i,k,j+1)-ph(i,k,j-1))                                    &
1714                    -9.*(ph(i,k,j+2)-ph(i,k,j-2))                                    &
1715                       +(ph(i,k,j+3)-ph(i,k,j-3))                                    &
1716                   +45.*(phb(i,k,j+1)-phb(i,k,j-1))                                  &
1717                    -9.*(phb(i,k,j+2)-phb(i,k,j-2))                                  &
1718                       +(phb(i,k,j+3)-phb(i,k,j-3))  )   )                
1721      ENDDO
1722      ENDDO
1724      k = kte
1725      DO i = i_start, itf
1726         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* (                                &
1727                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)                &
1728                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  ) )* (1./60.)*(  &
1729                    45.*(ph(i,k,j+1)-ph(i,k,j-1))                                               &
1730                    -9.*(ph(i,k,j+2)-ph(i,k,j-2))                                               &
1731                       +(ph(i,k,j+3)-ph(i,k,j-3))                                               &
1732                   +45.*(phb(i,k,j+1)-phb(i,k,j-1))                                             &
1733                    -9.*(phb(i,k,j+2)-phb(i,k,j-2))                                             &
1734                       +(phb(i,k,j+3)-phb(i,k,j-3))  )   )                
1736      ENDDO
1738    ENDDO
1740 !  4th order stencil for open or specified conditions two in form the boundary
1742    IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte )  THEN
1744      j = jds+2
1745      DO k = 2, kte-1
1746      DO i = i_start, itf
1747         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* (                     &
1748                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)                &
1749                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  ) )* (1./12.)*(  &
1750                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                    &
1751                       -(ph(i,k,j+2)-ph(i,k,j-2))                                    &
1752                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                  &
1753                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1756      ENDDO
1757      ENDDO
1759      k = kte
1760      DO i = i_start, itf
1761         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* (                              &
1762                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)              &
1763                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j) )* (1./12.)*(  &
1764                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                             &
1765                       -(ph(i,k,j+2)-ph(i,k,j-2))                                             &
1766                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                           &
1767                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1769      ENDDO
1771    END IF
1773    IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte )  THEN
1774      j = jde-3
1775      DO k = 2, kte-1
1776      DO i = i_start, itf
1777         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* (                  &
1778                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)              &
1779                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j) )* (1./12.)*(  &
1780                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                  &
1781                       -(ph(i,k,j+2)-ph(i,k,j-2))                                  &
1782                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                &
1783                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1786      ENDDO
1787      ENDDO
1789      k = kte
1790      DO i = i_start, itf
1791         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* (                              &
1792                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)              &
1793                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j) )* (1./12.)*(  &
1794                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                             &
1795                       -(ph(i,k,j+2)-ph(i,k,j-2))                                             &
1796                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                                           &
1797                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1799      ENDDO
1801    END IF
1803 !  2nd order stencil for open and specified conditions one row in from boundary
1805    IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte )  THEN
1807      j = jds+1
1808      DO k = 2, kte-1
1809      DO i = i_start, itf
1810         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
1811                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
1812                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
1813                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
1814                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
1815      ENDDO
1816      ENDDO
1818      k = kte
1819      DO i = i_start, itf
1820         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
1821                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
1822                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
1823                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
1824                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
1825      ENDDO
1827    END IF
1829    IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte )  THEN
1831      j = jde-2
1832      DO k = 2, kte-1
1833      DO i = i_start, itf
1834         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*          &
1835                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)*   &
1836                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))   &
1837                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvy(i,j  )*   &
1838                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1)) )
1839      ENDDO
1840      ENDDO
1842      k = kte
1843      DO i = i_start, itf
1844         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*                        &
1845                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)*    &
1846                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))               &
1847                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvy(i,j  )*    &
1848                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))              )
1849      ENDDO
1851    END IF
1853 !  x (u) advection
1855    i_start = its
1856    j_start = jts
1857    itf=MIN(ite,ide-1)
1858    jtf=MIN(jte,jde-1)
1860    IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
1861    IF (config_flags%open_xe .or. specified ) itf     = min(itf,ide-4)
1863    DO j = j_start, jtf
1865      DO k = 2, kte-1
1866      DO i = i_start, itf
1867         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*(                   &
1868                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)              &
1869                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))*msfux(i,j) )* (1./60.)*(  &
1870                    45.*(ph(i+1,k,j)-ph(i-1,k,j))                                  &
1871                    -9.*(ph(i+2,k,j)-ph(i-2,k,j))                                  &
1872                       +(ph(i+3,k,j)-ph(i-3,k,j))                                  &
1873                   +45.*(phb(i+1,k,j)-phb(i-1,k,j))                                &
1874                    -9.*(phb(i+2,k,j)-phb(i-2,k,j))                                &
1875                       +(phb(i+3,k,j)-phb(i-3,k,j))  )   )                
1876      ENDDO
1877      ENDDO
1879      k = kte
1880      DO i = i_start, itf
1881         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*(                             &
1882                  ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)            &
1883                   +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*(  &
1884                    45.*(ph(i+1,k,j)-ph(i-1,k,j))                                           &
1885                    -9.*(ph(i+2,k,j)-ph(i-2,k,j))                                           &
1886                       +(ph(i+3,k,j)-ph(i-3,k,j))                                           &
1887                   +45.*(phb(i+1,k,j)-phb(i-1,k,j))                                         &
1888                    -9.*(phb(i+2,k,j)-phb(i-2,k,j))                                         &
1889                       +(phb(i+3,k,j)-phb(i-3,k,j))  )     )
1890      ENDDO
1892    ENDDO
1894 !  4th order stencil two in from the boundary for open and specified conditions
1896    IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
1897      i = ids + 2
1898      DO j = j_start, jtf
1899        DO k = 2, kte-1
1900         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*(                   &
1901                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)              &
1902                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))*msfux(i,j) )* (1./12.)*(  &
1903                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                                  &
1904                       -(ph(i+2,k,j)-ph(i-2,k,j))                                  &
1905                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                                &
1906                       -(phb(i+2,k,j)-phb(i-2,k,j))  )   )                
1907        ENDDO
1908        k = kte
1909        ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*(                             &
1910                 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)            &
1911                  +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*(  &
1912                    8.*(ph(i+1,k,j)-ph(i-1,k,j))                                           &
1913                      -(ph(i+2,k,j)-ph(i-2,k,j))                                           &
1914                   +8.*(phb(i+1,k,j)-phb(i-1,k,j))                                         &
1915                      -(phb(i+2,k,j)-phb(i-2,k,j))  )     )
1917      ENDDO
1918    END IF
1920    IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
1921      i = ide-3
1922      DO j = j_start, jtf
1923        DO k = 2, kte-1
1924         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*(                   &
1925                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)              &
1926                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))*msfux(i,j) )* (1./12.)*(  &
1927                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                                  &
1928                       -(ph(i+2,k,j)-ph(i-2,k,j))                                  &
1929                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                                &
1930                       -(phb(i+2,k,j)-phb(i-2,k,j))  )   )                
1931        ENDDO
1932        k = kte
1933        ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*(                             &
1934                 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)            &
1935                  +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*(  &
1936                    8.*(ph(i+1,k,j)-ph(i-1,k,j))                                           &
1937                      -(ph(i+2,k,j)-ph(i-2,k,j))                                           &
1938                   +8.*(phb(i+1,k,j)-phb(i-1,k,j))                                         &
1939                      -(phb(i+2,k,j)-phb(i-2,k,j))  )     )
1941      ENDDO
1942    END IF
1944 !  2nd order stencil for open and specified conditions one in from the boundary
1946    IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
1948      i = ids + 1
1950      DO j = j_start, jtf
1951      DO k = 2, kte-1
1952         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*         &
1953                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)*  &
1954                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))  &
1955                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j)*  &
1956                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))  )
1957      ENDDO
1958      ENDDO
1960      k = kte
1961      DO j = j_start, jtf
1962         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*                        &
1963                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)*    &
1964                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))               &
1965                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(  i,j)*    &
1966                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))             )
1967      ENDDO
1969    END IF
1971    IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
1973      i = ide-2
1974      DO j = j_start, jtf
1975      DO k = 2, kte-1
1976         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*         &
1977                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)*  &
1978                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))  &
1979                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*msfux(i  ,j)*  &
1980                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))  )
1981      ENDDO
1982      ENDDO
1984      k = kte
1985      DO j = j_start, jtf
1986         ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*                        &
1987                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)*    &
1988                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))               &
1989                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*msfux(  i,j)*    &
1990                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))             )
1991      ENDDO
1993    END IF
1995    END IF  ! 6th order advection
1997 !  lateral open boundary conditions,
1998 !  start with north and south (y) boundaries
2000    i_start = its
2001    itf=MIN(ite,ide-1)
2003    !  south
2005    IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2007      j=jts
2009      DO k=2,kde
2010        kz = min(k,kde-1)
2011        DO i = its,itf
2012          vb =.5*( fnm(kz)*(v(i,kz  ,j+1)+v(i,kz  ,j  ))    &
2013                  +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j  )) )
2014          vl=amin1(vb,0.)
2015          ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*(      &
2016                               +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
2017        ENDDO
2018      ENDDO
2020    END IF
2022    ! north
2024    IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2026      j=jte-1
2028      DO k=2,kde
2029        kz = min(k,kde-1)
2030        DO i = its,itf
2031         vb=.5*( fnm(kz)*(v(i,kz  ,j+1)+v(i,kz  ,j))   &
2032                +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
2033         vr=amax1(vb,0.)
2034         ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*(      &
2035                    +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
2036        ENDDO
2037      ENDDO
2039    END IF
2041    !  now the east and west (y) boundaries
2043    j_start = its
2044    jtf=MIN(jte,jde-1)
2046    !  west
2048    IF ( (config_flags%open_xs) .and. its == ids ) THEN
2050      i=its
2052      DO j = jts,jtf
2053        DO k=2,kde-1
2054          kz = k
2055          ub =.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i  ,kz  ,j))     &
2056                  +fnp(kz)*(u(i+1,kz-1,j)+u(i  ,kz-1,j)) )
2057          ul=amin1(ub,0.)
2058          ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*(       &
2059                               +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2060        ENDDO
2062          k = kde
2063          kz = k
2064          ub =.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i  ,kz  ,j))     &
2065                  +fnp(kz)*(u(i+1,kz-1,j)+u(i  ,kz-1,j)) )
2066          ul=amin1(ub,0.)
2067          ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*(       &
2068                               +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2069      ENDDO
2071    END IF
2073    ! east
2075    IF ( (config_flags%open_xe) .and. ite == ide ) THEN
2077      i = ite-1
2079      DO j = jts,jtf
2080        DO k=2,kde-1
2081         kz = k
2082         ub=.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i,kz  ,j))  &
2083                +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2084         ur=amax1(ub,0.)
2085         ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
2086                    +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2087        ENDDO
2089         k = kde    
2090         kz = k-1
2091         ub=.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i,kz  ,j))   &
2092                +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2093         ur=amax1(ub,0.)
2094         ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*(  &
2095                    +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2097      ENDDO
2099    END IF
2101   END SUBROUTINE rhs_ph
2104 !-------------------------------------------------------------------------------
2106 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend,                &
2107                                          ph,alt,p,pb,al,php,cqu,cqv,     &
2108                                          muu,muv,mu,fnm,fnp,rdnw,        &
2109                                          cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
2110                                          msfvx,msfvy,msftx,msfty,        &
2111                                          config_flags, non_hydrostatic,  &
2112                                          top_lid,                        &
2113                                          ids, ide, jds, jde, kds, kde,   &
2114                                          ims, ime, jms, jme, kms, kme,   &
2115                                          its, ite, jts, jte, kts, kte   )
2117    IMPLICIT NONE
2118    
2119    ! Input data
2122    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2124    LOGICAL, INTENT (IN   ) :: non_hydrostatic, top_lid
2126    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2127                                        ims, ime, jms, jme, kms, kme, &
2128                                        its, ite, jts, jte, kts, kte 
2130    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::        &
2131                                                                      ph,  &
2132                                                                      alt, &
2133                                                                      al,  &
2134                                                                      p,   &
2135                                                                      pb,  &
2136                                                                      php, &
2137                                                                      cqu, &
2138                                                                      cqv
2141    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::           &
2142                                                                     ru_tend, &
2143                                                                     rv_tend
2145    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: muu, muv, mu,    &
2146                                                             msfux, msfuy, &
2147                                                             msfvx, msfvy, &
2148                                                             msftx, msfty
2150    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, fnm, fnp
2152    REAL,  INTENT(IN   ) :: rdx, rdy, cf1, cf2, cf3
2154    INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
2155    REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
2156    REAL :: dpx, dpy
2158    LOGICAL :: specified
2160 !<DESCRIPTION>
2162 !  horizontal_pressure_gradient calculates the 
2163 !  horizontal pressure gradient terms for the large-timestep tendency 
2164 !  in the horizontal momentum equations (u,v).
2166 !</DESCRIPTION>
2168    specified = .false.
2169    if(config_flags%specified .or. config_flags%nested) specified = .true.
2171 !  Notes on map scale factors:
2172 !  Calculates the pressure gradient terms in ADT eqns 44 and 45
2173 !  With upper rho -> 'mu', these are:
2174 !  Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
2175 !  Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
2177 !  As we are on nu, rather than height, surfaces:
2179 !  mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
2180 !           + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
2182 !  mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
2183 !           + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
2185 ! start with the north-south (y) pressure gradient
2187    itf=MIN(ite,ide-1)
2188    jtf=jte
2189    ktf=MIN(kte,kde-1)
2190    i_start = its
2191    j_start = jts
2192    IF ( (config_flags%open_ys .or. specified .or. &
2193          config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2194    IF ( (config_flags%open_ye .or. specified .or. &
2195          config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2197    DO j = j_start, jtf
2199      IF ( non_hydrostatic )  THEN
2201         k=1
2203         DO i = i_start, itf
2204           dpn(i,k) = .5*( cf1*(p(i,k  ,j-1)+p(i,k  ,j))   &
2205                          +cf2*(p(i,k+1,j-1)+p(i,k+1,j))   &
2206                          +cf3*(p(i,k+2,j-1)+p(i,k+2,j))  )
2207           dpn(i,kde) = 0.
2208         ENDDO
2209         IF (top_lid) THEN
2210           DO i = i_start, itf
2211             dpn(i,kde) = .5*( cf1*(p(i,kde-1,j-1)+p(i,kde-1,j))   &
2212                              +cf2*(p(i,kde-2,j-1)+p(i,kde-2,j))   &
2213                              +cf3*(p(i,kde-3,j-1)+p(i,kde-3,j))  )
2214           ENDDO
2215         ENDIF
2216                
2217         DO k=2,ktf
2218           DO i = i_start, itf
2219             dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j-1)+p(i,k  ,j))  &
2220                            +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2221           END DO
2222         END DO
2224 !       ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2225 !       [alt, al are 1/rho terms; muv, mu are NOT coupled]
2226         DO K=1,ktf
2227           DO i = i_start, itf
2228             ! Here are mu dp/dy terms 1-3 
2229             dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*(                 &
2230                      (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1))  &
2231                     +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p (i,k,j)-p (i,k,j-1))  &
2232                     +(al (i,k  ,j)+al (i,k  ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2233             ! Here is mu dp/dy term 4 
2234             dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2235                 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2236             rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2237           END DO
2238         END DO
2240      ELSE
2242 !       ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2243 !       [alt, al are 1/rho terms; muv, mu are NOT coupled]
2244         DO K=1,ktf
2245           DO i = i_start, itf
2246             ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2247             dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*(                 &
2248                      (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1))  &
2249                     +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p (i,k,j)-p (i,k,j-1))  &
2250                     +(al (i,k  ,j)+al (i,k  ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2251             rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2252           END DO
2253         END DO
2255      END IF
2257    ENDDO
2259 !  now the east-west (x) pressure gradient
2261    itf=ite
2262    jtf=MIN(jte,jde-1)
2263    ktf=MIN(kte,kde-1)
2264    i_start = its
2265    j_start = jts
2266    IF ( (config_flags%open_xs .or. specified .or. &
2267            config_flags%nested ) .and. its == ids ) i_start = its+1
2268    IF ( (config_flags%open_xe .or. specified .or. &
2269            config_flags%nested ) .and. ite == ide ) itf = itf-1
2270    IF ( config_flags%periodic_x ) i_start = its
2271    IF ( config_flags%periodic_x ) itf=ite
2273    DO j = j_start, jtf
2275      IF ( non_hydrostatic )  THEN
2277         k=1
2279         DO i = i_start, itf
2280           dpn(i,k) = .5*( cf1*(p(i-1,k  ,j)+p(i,k  ,j))   &
2281                          +cf2*(p(i-1,k+1,j)+p(i,k+1,j))   &
2282                          +cf3*(p(i-1,k+2,j)+p(i,k+2,j))  )
2283           dpn(i,kde) = 0.
2284         ENDDO
2285         IF (top_lid) THEN
2286           DO i = i_start, itf
2287             dpn(i,kde) = .5*( cf1*(p(i-1,kde-1,j)+p(i,kde-1,j))   &
2288                              +cf2*(p(i-1,kde-2,j)+p(i,kde-2,j))   &
2289                              +cf3*(p(i-1,kde-3,j)+p(i,kde-3,j))  )
2290           ENDDO
2291         ENDIF
2292                
2293         DO k=2,ktf
2294           DO i = i_start, itf
2295             dpn(i,k) = .5*( fnm(k)*(p(i-1,k  ,j)+p(i,k  ,j))  &
2296                            +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2297           END DO
2298         END DO
2300 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2301 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2302         DO K=1,ktf
2303           DO i = i_start, itf
2304             ! Here are mu dp/dy terms 1-3
2305             dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*(                    &
2306                         (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j))  &
2307                        +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p (i,k,j)-p (i-1,k,j))  &
2308                        +(al (i,k  ,j)+al (i-1,k  ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2309             ! Here is mu dp/dy term 4
2310             dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2311                 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2312             ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2313           END DO
2314         END DO
2316      ELSE
2318 !       ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2319 !       [alt, al are 1/rho terms; muu, mu are NOT coupled]
2320         DO K=1,ktf
2321           DO i = i_start, itf
2322             ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2323             dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*(                    &
2324                         (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j))  &
2325                        +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p (i,k,j)-p (i-1,k,j))  &
2326                        +(al (i,k  ,j)+al (i-1,k  ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2327             ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2328           END DO
2329         END DO
2331      END IF
2333    ENDDO
2335 END SUBROUTINE horizontal_pressure_gradient
2337 !-------------------------------------------------------------------------------
2339 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub,       &
2340                       rdnw, rdn, g, msftx, msfty,     &
2341                       ids, ide, jds, jde, kds, kde,   &
2342                       ims, ime, jms, jme, kms, kme,   &
2343                       its, ite, jts, jte, kts, kte   )
2345    IMPLICIT NONE
2346    
2347    ! Input data
2349    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2350                                        ims, ime, jms, jme, kms, kme, &
2351                                        its, ite, jts, jte, kts, kte 
2353    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   p
2354    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::   cqw
2357    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::  rw_tend
2359    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mub, mu, msftx, msfty
2361    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, rdn
2363    REAL,  INTENT(IN   ) :: g
2365    INTEGER :: itf, jtf, i, j, k
2366    REAL    :: cq1, cq2
2369 !<DESCRIPTION>
2371 !  pg_buoy_w calculates the 
2372 !  vertical pressure gradient and buoyancy terms for the large-timestep 
2373 !  tendency in the vertical momentum equation.
2375 !</DESCRIPTION>
2377 !  BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2379 !  Map scale factor notes
2380 !  ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2381 !  Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2382 !  term 6: +(g/my) partial dp'/dnu
2383 !  term 7: -(g/my) mu'
2385 !  For moisture-free atmosphere, cq1=1, cq2=0
2386 !  => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2388    itf=MIN(ite,ide-1)
2389    jtf=MIN(jte,jde-1)
2391    DO j = jts,jtf
2393      k=kde
2394      DO i=its,itf
2395        cq1 = 1./(1.+cqw(i,k-1,j))
2396        cq2 = cqw(i,k-1,j)*cq1
2397        rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*(      &
2398                         cq1*2.*rdnw(k-1)*(  -p(i,k-1,j))  &
2399                         -mu(i,j)-cq2*mub(i,j)            )
2400      END DO
2402      DO k = 2, kde-1
2403      DO i = its,itf
2404       cq1 = 1./(1.+cqw(i,k,j))
2405       cq2 = cqw(i,k,j)*cq1
2406       cqw(i,k,j) = cq1
2407       rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*(      &
2408                        cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j))  &
2409                        -mu(i,j)-cq2*mub(i,j)            )
2410      END DO
2411      ENDDO           
2414    ENDDO
2416 END SUBROUTINE pg_buoy_w
2418 !-------------------------------------------------------------------------------
2420 SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2421                       u, v, ww, w, mut, rdnw,         &
2422                       rdx, rdy, msfux, msfuy,         &
2423                       msfvx, msfvy, dt,               &
2424                       config_flags,                   &
2425                       ids, ide, jds, jde, kds, kde,   &
2426                       ims, ime, jms, jme, kms, kme,   &
2427                       its, ite, jts, jte, kts, kte   )
2429    USE module_llxy
2430    IMPLICIT NONE
2432    ! Input data
2434    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
2436    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2437                                        ims, ime, jms, jme, kms, kme, &
2438                                        its, ite, jts, jte, kts, kte
2440    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   u, v, ww, w
2442    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::  rw_tend
2444    REAL, INTENT(OUT) ::  max_vert_cfl
2445    REAL, INTENT(OUT) ::  max_horiz_cfl
2446    REAL              ::  horiz_cfl
2448    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mut
2450    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw
2452    REAL, INTENT(IN)    :: dt
2453    REAL, INTENT(IN)    :: rdx, rdy
2454    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux, msfuy
2455    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfvx, msfvy
2457    REAL                :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2459    INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2460    INTEGER :: some
2461    CHARACTER*512 :: temp
2463    CHARACTER (LEN=256) :: time_str
2464    CHARACTER (LEN=256) :: grid_str
2466    integer :: total
2467    REAL :: msfuxt , msfxffl
2468    
2469 !<DESCRIPTION>
2471 !  w_damp computes a damping term for the vertical velocity when the
2472 !  vertical Courant number is too large.  This was found to be preferable to 
2473 !  decreasing the timestep or increasing the diffusion in real-data applications
2474 !  that produced potentially-unstable large vertical velocities because of
2475 !  unphysically large heating rates coming from the cumulus parameterization 
2476 !  schemes run at moderately high resolutions (dx ~ O(10) km).
2478 !  Additionally, w_damp returns the maximum cfl values due to vertical motion and
2479 !  horizontal motion.  These values are returned via the max_vert_cfl and 
2480 !  max_horiz_cfl variables.  (Added by T. Hutchinson, WSI, 3/5/2007)
2482 !</DESCRIPTION>
2484    itf=MIN(ite,ide-1)
2485    jtf=MIN(jte,jde-1)
2487    some = 0
2488    max_vert_cfl = 0.
2489    max_horiz_cfl = 0.
2490    total = 0
2492    IF(config_flags%map_proj == PROJ_CASSINI ) then
2493      msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad) 
2494    END IF
2496    IF ( config_flags%w_damping == 1 ) THEN
2497      DO j = jts,jtf
2499      DO k = 2, kde-1
2500      DO i = its,itf
2501 #if 1
2502         IF(config_flags%map_proj == PROJ_CASSINI ) then
2503            msfuxt = MIN(msfux(i,j), msfxffl)
2504         ELSE
2505            msfuxt = msfux(i,j)
2506         END IF
2507         vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2509         IF ( vert_cfl > max_vert_cfl ) THEN
2510            max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k 
2511            maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2512         ENDIF
2513         
2514         horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt),                          &
2515              abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2516         if (horiz_cfl > max_horiz_cfl) then
2517            max_horiz_cfl = horiz_cfl
2518         endif
2519         
2520         if(vert_cfl .gt. w_beta)then
2521 #else
2522 ! restructure to get rid of divide
2524 ! This had been used for efficiency, but with the addition of returning the cfl values, 
2525 !   the old version (above) was reinstated.  (T. Hutchinson, 3/5/2007)
2527         cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2528         cf_d = abs(mut(i,j))
2529         if(cf_n .gt. cf_d*w_beta )then
2530 #endif
2532            WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2533            CALL wrf_debug ( 100 , TRIM(temp) )
2534            if ( vert_cfl > 2. ) some = some + 1
2535            rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(vert_cfl-w_beta)*mut(i,j)
2536         endif
2537      END DO
2538      ENDDO
2539      ENDDO
2540    ELSE
2541 ! just print
2542      DO j = jts,jtf
2544      DO k = 2, kde-1
2545      DO i = its,itf
2547 #if 1
2548         IF(config_flags%map_proj == PROJ_CASSINI ) then
2549            msfuxt = MIN(msfux(i,j), msfxffl)
2550         ELSE
2551            msfuxt = msfux(i,j)
2552         END IF
2553         vert_cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2554         
2555         IF ( vert_cfl > max_vert_cfl ) THEN
2556            max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k 
2557            maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2558         ENDIF
2559         
2560         horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt),                          &
2561              abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2563         if (horiz_cfl > max_horiz_cfl) then
2564            max_horiz_cfl = horiz_cfl
2565         endif
2566         
2567         if(vert_cfl .gt. w_beta)then
2568 #else
2569 ! restructure to get rid of divide
2571 ! This had been used for efficiency, but with the addition of returning the cfl values, 
2572 !   the old version (above) was reinstated.  (T. Hutchinson, 3/5/2007)
2574         cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2575         cf_d = abs(mut(i,j))
2576         if(cf_n .gt. cf_d*w_beta )then
2577 #endif
2578            WRITE(temp,*)i,j,k,' vert_cfl,w,d(eta)=',vert_cfl,w(i,k,j),-1./rdnw(k)
2579            CALL wrf_debug ( 100 , TRIM(temp) )
2580            if ( vert_cfl > 2. ) some = some + 1
2581         endif
2582      END DO
2583      ENDDO
2584      ENDDO
2585    ENDIF
2586    IF ( some .GT. 0 ) THEN
2587      CALL get_current_time_string( time_str )
2588      CALL get_current_grid_name( grid_str )
2589      WRITE(temp,*)some,                                            &
2590             ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2591      CALL wrf_debug ( 0 , TRIM(temp) )
2592      WRITE(temp,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' vert_cfl,w,d(eta)=',max_vert_cfl, &
2593                              maxdub,maxdeta
2594      CALL wrf_debug ( 0 , TRIM(temp) )
2595    ENDIF
2597 END SUBROUTINE w_damp
2599 !-------------------------------------------------------------------------------
2601 SUBROUTINE horizontal_diffusion ( name, field, tendency, mu,           &
2602                                   config_flags,                        &
2603                                   msfux, msfuy, msfvx, msfvx_inv,      &
2604                                   msfvy, msftx, msfty,                 &
2605                                   khdif, xkmhd, rdx, rdy,              &
2606                                   ids, ide, jds, jde, kds, kde,        &
2607                                   ims, ime, jms, jme, kms, kme,        &
2608                                   its, ite, jts, jte, kts, kte        )
2610    IMPLICIT NONE
2611    
2612    ! Input data
2614    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2616    INTEGER ,        INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2617                                      ims, ime, jms, jme, kms, kme, &
2618                                      its, ite, jts, jte, kts, kte
2620    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
2622    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field, xkmhd
2624    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2626    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
2628    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
2629                                                                     msfuy,      &
2630                                                                     msfvx,      &
2631                                                                     msfvx_inv,  &
2632                                                                     msfvy,      &
2633                                                                     msftx,      &
2634                                                                     msfty
2636    REAL ,                                      INTENT(IN   ) :: rdx,       &
2637                                                                 rdy,       &
2638                                                                 khdif
2640    ! Local data
2641    
2642    INTEGER :: i, j, k, itf, jtf, ktf
2644    INTEGER :: i_start, i_end, j_start, j_end
2646    REAL :: mrdx, mkrdxm, mkrdxp, &
2647            mrdy, mkrdym, mkrdyp
2649    LOGICAL :: specified
2651 !<DESCRIPTION>
2653 !  horizontal_diffusion computes the horizontal diffusion tendency
2654 !  on model horizontal coordinate surfaces.
2656 !</DESCRIPTION>
2658    specified = .false.
2659    if(config_flags%specified .or. config_flags%nested) specified = .true.
2661    ktf=MIN(kte,kde-1)
2662    
2663    IF (name .EQ. 'u') THEN
2665       i_start = its
2666       i_end   = ite
2667       j_start = jts
2668       j_end   = MIN(jte,jde-1)
2670       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2671       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-1,ite)
2672       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2673       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2674       IF ( config_flags%periodic_x ) i_start = its
2675       IF ( config_flags%periodic_x ) i_end = ite
2678       DO j = j_start, j_end
2679       DO k=kts,ktf
2680       DO i = i_start, i_end
2682          ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2683          ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2685          mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
2686          mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
2687          mrdx=msfux(i,j)*msfuy(i,j)*rdx 
2688          mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2689                 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2690                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2691          mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2692                 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
2693                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2694          ! need to do four-corners (t) for diffusion coefficient as there are
2695          ! no values at u,v points
2696          ! msfuy - has to be y as part of d/dY
2697          !         has to be u as we're at a u point
2698          mrdy=msfux(i,j)*msfuy(i,j)*rdy 
2700          ! correctly averaged version of rho~ * m^2 * 
2701          !    [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2702             tendency(i,k,j)=tendency(i,k,j)+( &
2703                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j))  &
2704                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2705                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  ))  &
2706                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2707       ENDDO
2708       ENDDO
2709       ENDDO
2710    
2711    ELSE IF (name .EQ. 'v')THEN
2713       i_start = its
2714       i_end   = MIN(ite,ide-1)
2715       j_start = jts
2716       j_end   = jte
2718       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2719       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2720       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2721       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-1,jte)
2722       IF ( config_flags%periodic_x ) i_start = its
2723       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2724       IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2725       IF ( config_flags%polar ) j_end   = MIN(jde-1,jte)
2727       DO j = j_start, j_end
2728       DO k=kts,ktf
2729       DO i = i_start, i_end
2731          mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )*    &
2732                 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2733                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2734          mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )*    &
2735                 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
2736                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2737          mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2738          mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2739          mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2740          mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2742             tendency(i,k,j)=tendency(i,k,j)+( &
2743                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j))  &
2744                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2745                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  ))  &
2746                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2747       ENDDO
2748       ENDDO
2749       ENDDO
2750    
2751    ELSE IF (name .EQ. 'w')THEN
2753       i_start = its
2754       i_end   = MIN(ite,ide-1)
2755       j_start = jts
2756       j_end   = MIN(jte,jde-1)
2758       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2759       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2760       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2761       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2762       IF ( config_flags%periodic_x ) i_start = its
2763       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2765       DO j = j_start, j_end
2766       DO k=kts+1,ktf
2767       DO i = i_start, i_end
2769          mkrdxm=(msfux(i,j)/msfuy(i,j))*   &
2770                 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
2771                 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2772          mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*   &
2773                 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
2774                 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2775          mrdx=msftx(i,j)*msfty(i,j)*rdx
2776 !         mkrdym=(msfvy(i,j)/msfvx(i,j))*   &
2777          mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*   &
2778                 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
2779                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2780 !         mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*   &
2781          mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*   &
2782                 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
2783                 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2784          mrdy=msftx(i,j)*msfty(i,j)*rdy
2786             tendency(i,k,j)=tendency(i,k,j)+( &
2787                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j)) &
2788                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2789                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  )) &
2790                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2791       ENDDO
2792       ENDDO
2793       ENDDO
2794    
2795    ELSE
2798       i_start = its
2799       i_end   = MIN(ite,ide-1)
2800       j_start = jts
2801       j_end   = MIN(jte,jde-1)
2803       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2804       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2805       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2806       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2807       IF ( config_flags%periodic_x ) i_start = its
2808       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2810       DO j = j_start, j_end
2811       DO k=kts,ktf
2812       DO i = i_start, i_end
2814          mkrdxm=(msfux(i,j)/msfuy(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*0.5*(mu(i,j)+mu(i-1,j))*rdx
2815          mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*0.5*(mu(i+1,j)+mu(i,j))*rdx
2816          mrdx=msftx(i,j)*msfty(i,j)*rdx
2817 !         mkrdym=(msfvy(i,j)/msfvx(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
2818          mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
2819 !         mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
2820          mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
2821          mrdy=msftx(i,j)*msfty(i,j)*rdy
2823             tendency(i,k,j)=tendency(i,k,j)+( &
2824                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j))  &
2825                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2826                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  ))  &
2827                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2828       ENDDO
2829       ENDDO
2830       ENDDO
2831            
2832    ENDIF
2834 END SUBROUTINE horizontal_diffusion
2836 !-----------------------------------------------------------------------------------------
2838 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu,           &
2839                                        config_flags, base_3d,               &
2840                                        msfux, msfuy, msfvx, msfvx_inv,      &
2841                                        msfvy, msftx, msfty,                 &
2842                                        khdif, xkmhd, rdx, rdy,              &
2843                                        ids, ide, jds, jde, kds, kde,        &
2844                                        ims, ime, jms, jme, kms, kme,        &
2845                                        its, ite, jts, jte, kts, kte        )
2847    IMPLICIT NONE
2848    
2849    ! Input data
2850    
2851    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2853    INTEGER ,        INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2854                                      ims, ime, jms, jme, kms, kme, &
2855                                      its, ite, jts, jte, kts, kte
2857    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
2859    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field, &
2860                                                                       xkmhd, &
2861                                                                       base_3d
2863    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2865    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
2867    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
2868                                                                     msfuy,      &
2869                                                                     msfvx,      &
2870                                                                     msfvx_inv,  &
2871                                                                     msfvy,      &
2872                                                                     msftx,      &
2873                                                                     msfty
2875    REAL ,                                      INTENT(IN   ) :: rdx,       &
2876                                                                 rdy,       &
2877                                                                 khdif
2879    ! Local data
2880    
2881    INTEGER :: i, j, k, itf, jtf, ktf
2883    INTEGER :: i_start, i_end, j_start, j_end
2885    REAL :: mrdx, mkrdxm, mkrdxp, &
2886            mrdy, mkrdym, mkrdyp
2888    LOGICAL :: specified
2890 !<DESCRIPTION>
2892 !  horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2893 !  on model horizontal coordinate surfaces.  This routine computes diffusion
2894 !  a perturbation scalar (field-base_3d).
2896 !</DESCRIPTION>
2898    specified = .false.
2899    if(config_flags%specified .or. config_flags%nested) specified = .true.
2901    ktf=MIN(kte,kde-1)
2902    
2903       i_start = its
2904       i_end   = MIN(ite,ide-1)
2905       j_start = jts
2906       j_end   = MIN(jte,jde-1)
2908       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2909       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2910       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2911       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2912       IF ( config_flags%periodic_x ) i_start = its
2913       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2915       DO j = j_start, j_end
2916       DO k=kts,ktf
2917       DO i = i_start, i_end
2919          mkrdxm=(msfux(i,j)/msfuy(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*0.5*(mu(i,j)+mu(i-1,j))*rdx
2920          mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*0.5*(mu(i+1,j)+mu(i,j))*rdx
2921          mrdx=msftx(i,j)*msfty(i,j)*rdx
2922 !         mkrdym=(msfvy(i,j)/msfvx(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
2923 !         mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
2924          mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*(mu(i,j)+mu(i,j-1))*rdy
2925          mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*(mu(i,j+1)+mu(i,j))*rdy
2926          mrdy=msftx(i,j)*msfty(i,j)*rdy
2928             tendency(i,k,j)=tendency(i,k,j)+(                        &
2929                     mrdx*( mkrdxp*(   field(i+1,k,j)  -field(i  ,k,j)      &
2930                                    -base_3d(i+1,k,j)+base_3d(i  ,k,j) )    &
2931                           -mkrdxm*(   field(i  ,k,j)  -field(i-1,k,j)      &
2932                                    -base_3d(i  ,k,j)+base_3d(i-1,k,j) )  ) &
2933                    +mrdy*( mkrdyp*(   field(i,k,j+1)  -field(i,k,j  )      &
2934                                    -base_3d(i,k,j+1)+base_3d(i,k,j  ) )    &
2935                           -mkrdym*(   field(i,k,j  )  -field(i,k,j-1)      &
2936                                    -base_3d(i,k,j  )+base_3d(i,k,j-1) )  ) &
2937                                                                          ) 
2938       ENDDO
2939       ENDDO
2940       ENDDO
2942 END SUBROUTINE horizontal_diffusion_3dmp
2944 !-----------------------------------------------------------------------------------------
2946 SUBROUTINE vertical_diffusion ( name, field, tendency,        &
2947                                 config_flags,                 &
2948                                 alt, mut, rdn, rdnw, kvdif,   &
2949                                 ids, ide, jds, jde, kds, kde, &
2950                                 ims, ime, jms, jme, kms, kme, &
2951                                 its, ite, jts, jte, kts, kte )
2954    IMPLICIT NONE
2955    
2956    ! Input data
2957    
2958    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2960    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2961                                  ims, ime, jms, jme, kms, kme, &
2962                                  its, ite, jts, jte, kts, kte
2964    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
2966    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
2967                                                INTENT(IN   ) :: field,    &
2968                                                                 alt
2970    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2972    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut
2974    REAL , DIMENSION( kms:kme ) ,                   INTENT(IN   ) :: rdn, rdnw
2976    REAL ,                                      INTENT(IN   ) :: kvdif
2977    
2978    ! Local data
2979    
2980    INTEGER :: i, j, k, itf, jtf, ktf
2981    INTEGER :: i_start, i_end, j_start, j_end
2983    REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
2984    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2986    REAL :: rdz
2988    LOGICAL :: specified
2990 !<DESCRIPTION>
2992 !  vertical_diffusion
2993 !  computes vertical diffusion tendency.
2995 !</DESCRIPTION>
2997    specified = .false.
2998    if(config_flags%specified .or. config_flags%nested) specified = .true.
3000    ktf=MIN(kte,kde-1)
3001    
3002    IF (name .EQ. 'w')THEN
3004    
3005    i_start = its
3006    i_end   = MIN(ite,ide-1)
3007    j_start = jts
3008    j_end   = MIN(jte,jde-1)
3010 j_loop_w : DO j = j_start, j_end
3012      DO k=kts,ktf-1
3013        DO i = i_start, i_end
3014           vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
3015        ENDDO
3016      ENDDO
3018      DO i = i_start, i_end
3019        vflux(i,ktf)=0.
3020      ENDDO
3022      DO k=kts+1,ktf
3023        DO i = i_start, i_end
3024             tendency(i,k,j)=tendency(i,k,j)                                         &
3025                               +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j)))  &
3026                                          *(vflux(i,k)-vflux(i,k-1))
3027        ENDDO
3028      ENDDO
3030     ENDDO j_loop_w
3032    ELSE IF(name .EQ. 'm')THEN
3034      i_start = its
3035      i_end   = MIN(ite,ide-1)
3036      j_start = jts
3037      j_end   = MIN(jte,jde-1)
3039 j_loop_s : DO j = j_start, j_end
3041      DO k=kts,ktf-1
3042        DO i = i_start, i_end
3043          vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j)))   &
3044                   *(field(i,k+1,j)-field(i,k,j))
3045        ENDDO
3046      ENDDO
3048      DO i = i_start, i_end
3049        vflux(i,0)=vflux(i,1)
3050      ENDDO
3052      DO i = i_start, i_end
3053        vflux(i,ktf)=0.
3054      ENDDO
3056      DO k=kts,ktf
3057        DO i = i_start, i_end
3058          tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j)  &
3059                 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3060        ENDDO
3061      ENDDO
3063  ENDDO j_loop_s
3065    ENDIF
3067 END SUBROUTINE vertical_diffusion
3070 !-------------------------------------------------------------------------------
3072 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
3073                                    base,                          &
3074                                    alt, mut, rdn, rdnw, kvdif,    &
3075                                    ids, ide, jds, jde, kds, kde,  &
3076                                    ims, ime, jms, jme, kms, kme,  &
3077                                    its, ite, jts, jte, kts, kte  )
3080    IMPLICIT NONE
3081    
3082    ! Input data
3083    
3084    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
3086    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3087                                  ims, ime, jms, jme, kms, kme, &
3088                                  its, ite, jts, jte, kts, kte
3090    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
3091                                                INTENT(IN   ) :: field,    &
3092                                                                 alt
3094    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3096    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut
3098    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn,  &
3099                                                                   rdnw, &
3100                                                                   base
3102    REAL ,                                      INTENT(IN   ) :: kvdif
3103    
3104    ! Local data
3105    
3106    INTEGER :: i, j, k, itf, jtf, ktf
3107    INTEGER :: i_start, i_end, j_start, j_end
3109    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3111    REAL :: rdz
3113    LOGICAL :: specified
3115 !<DESCRIPTION>
3117 !  vertical_diffusion_mp
3118 !  computes vertical diffusion tendency of a perturbation variable
3119 !  (field-base).  Note that base as a 1D (k) field.
3121 !</DESCRIPTION>
3123    specified = .false.
3124    if(config_flags%specified .or. config_flags%nested) specified = .true.
3126    ktf=MIN(kte,kde-1)
3127    
3128      i_start = its
3129      i_end   = MIN(ite,ide-1)
3130      j_start = jts
3131      j_end   = MIN(jte,jde-1)
3133 j_loop_s : DO j = j_start, j_end
3135      DO k=kts,ktf-1
3136        DO i = i_start, i_end
3137          vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j)))   &
3138                     *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
3139        ENDDO
3140      ENDDO
3142      DO i = i_start, i_end
3143        vflux(i,0)=vflux(i,1)
3144      ENDDO
3146      DO i = i_start, i_end
3147        vflux(i,ktf)=0.
3148      ENDDO
3150      DO k=kts,ktf
3151        DO i = i_start, i_end
3152          tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j)  &
3153                 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3154        ENDDO
3155      ENDDO
3157  ENDDO j_loop_s
3159 END SUBROUTINE vertical_diffusion_mp
3162 !-------------------------------------------------------------------------------
3164 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
3165                                      base_3d,                       &
3166                                      alt, mut, rdn, rdnw, kvdif,    &
3167                                      ids, ide, jds, jde, kds, kde,  &
3168                                      ims, ime, jms, jme, kms, kme,  &
3169                                      its, ite, jts, jte, kts, kte  )
3172    IMPLICIT NONE
3173    
3174    ! Input data
3175    
3176    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
3178    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3179                                  ims, ime, jms, jme, kms, kme, &
3180                                  its, ite, jts, jte, kts, kte
3182    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
3183                                                INTENT(IN   ) :: field,    &
3184                                                                 alt,      &
3185                                                                 base_3d
3187    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3189    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut
3191    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn,  &
3192                                                                   rdnw
3194    REAL ,                                      INTENT(IN   ) :: kvdif
3195    
3196    ! Local data
3197    
3198    INTEGER :: i, j, k, itf, jtf, ktf
3199    INTEGER :: i_start, i_end, j_start, j_end
3201    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3203    REAL :: rdz
3205    LOGICAL :: specified
3207 !<DESCRIPTION>
3209 !  vertical_diffusion_3dmp
3210 !  computes vertical diffusion tendency of a perturbation variable
3211 !  (field-base_3d).  
3213 !</DESCRIPTION>
3215    specified = .false.
3216    if(config_flags%specified .or. config_flags%nested) specified = .true.
3218    ktf=MIN(kte,kde-1)
3219    
3220      i_start = its
3221      i_end   = MIN(ite,ide-1)
3222      j_start = jts
3223      j_end   = MIN(jte,jde-1)
3225 j_loop_s : DO j = j_start, j_end
3227      DO k=kts,ktf-1
3228        DO i = i_start, i_end
3229          vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j)))   &
3230                     *(   field(i,k+1,j)  -field(i,k,j)               &
3231                       -base_3d(i,k+1,j)+base_3d(i,k,j) )
3232        ENDDO
3233      ENDDO
3235      DO i = i_start, i_end
3236        vflux(i,0)=vflux(i,1)
3237      ENDDO
3239      DO i = i_start, i_end
3240        vflux(i,ktf)=0.
3241      ENDDO
3243      DO k=kts,ktf
3244        DO i = i_start, i_end
3245          tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j)  &
3246                 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3247        ENDDO
3248      ENDDO
3250  ENDDO j_loop_s
3252 END SUBROUTINE vertical_diffusion_3dmp
3255 !-------------------------------------------------------------------------------
3258 SUBROUTINE vertical_diffusion_u ( field, tendency,              &
3259                                   config_flags, u_base,         &
3260                                   alt, muu, rdn, rdnw, kvdif,   &
3261                                   ids, ide, jds, jde, kds, kde, &
3262                                   ims, ime, jms, jme, kms, kme, &
3263                                   its, ite, jts, jte, kts, kte )
3266    IMPLICIT NONE
3267    
3268    ! Input data
3269    
3270    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
3272    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3273                                  ims, ime, jms, jme, kms, kme, &
3274                                  its, ite, jts, jte, kts, kte
3276    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
3277                                                INTENT(IN   ) :: field,    &
3278                                                                 alt
3280    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3282    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: muu
3284    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn, rdnw, u_base
3286    REAL ,                                      INTENT(IN   ) :: kvdif
3287    
3288    ! Local data
3289    
3290    INTEGER :: i, j, k, itf, jtf, ktf
3291    INTEGER :: i_start, i_end, j_start, j_end
3293    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3295    REAL :: rdz, zz
3297    LOGICAL :: specified
3299 !<DESCRIPTION>
3301 !  vertical_diffusion_u computes vertical diffusion tendency for 
3302 !  the u momentum equation.  This routine assumes a constant eddy
3303 !  viscosity kvdif.
3305 !</DESCRIPTION>
3307    specified = .false.
3308    if(config_flags%specified .or. config_flags%nested) specified = .true.
3310    ktf=MIN(kte,kde-1)
3312       i_start = its
3313       i_end   = ite
3314       j_start = jts
3315       j_end   = MIN(jte,jde-1)
3317       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3318       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-1,ite)
3319       IF ( config_flags%periodic_x ) i_start = its
3320       IF ( config_flags%periodic_x ) i_end = ite
3323 j_loop_u : DO j = j_start, j_end
3325      DO k=kts,ktf-1
3326        DO i = i_start, i_end
3327          vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i  ,k  ,j)      &
3328                                         +alt(i-1,k  ,j)      &
3329                                         +alt(i  ,k+1,j)      &
3330                                         +alt(i-1,k+1,j) ) )  &
3331                              *(field(i,k+1,j)-field(i,k,j)   &
3332                                -u_base(k+1)   +u_base(k)  )
3333        ENDDO
3334      ENDDO
3336      DO i = i_start, i_end
3337        vflux(i,0)=vflux(i,1)
3338      ENDDO
3340      DO i = i_start, i_end
3341        vflux(i,ktf)=0.
3342      ENDDO
3344      DO k=kts,ktf-1
3345        DO i = i_start, i_end
3346          tendency(i,k,j)=tendency(i,k,j)+                             &
3347                 g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3348                               (vflux(i,k)-vflux(i,k-1))
3349        ENDDO
3350      ENDDO
3352  ENDDO j_loop_u
3353    
3354 END SUBROUTINE vertical_diffusion_u
3356 !-------------------------------------------------------------------------------
3359 SUBROUTINE vertical_diffusion_v ( field, tendency,              &
3360                                   config_flags, v_base,         &
3361                                   alt, muv, rdn, rdnw, kvdif,   &
3362                                   ids, ide, jds, jde, kds, kde, &
3363                                   ims, ime, jms, jme, kms, kme, &
3364                                   its, ite, jts, jte, kts, kte )
3367    IMPLICIT NONE
3368    
3369    ! Input data
3370    
3371    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
3373    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3374                                  ims, ime, jms, jme, kms, kme, &
3375                                  its, ite, jts, jte, kts, kte
3377    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
3378                                                INTENT(IN   ) :: field,    &
3379                                                                 alt
3380    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn, rdnw, v_base
3382    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3384    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: muv
3386    REAL ,                                      INTENT(IN   ) :: kvdif
3387    
3388    ! Local data
3389    
3390    INTEGER :: i, j, k, itf, jtf, ktf, jm1
3391    INTEGER :: i_start, i_end, j_start, j_end
3393    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3395    REAL :: rdz, zz
3397    LOGICAL :: specified
3399 !<DESCRIPTION>
3401 !  vertical_diffusion_v computes vertical diffusion tendency for 
3402 !  the v momentum equation.  This routine assumes a constant eddy
3403 !  viscosity kvdif.
3405 !</DESCRIPTION>
3407    specified = .false.
3408    if(config_flags%specified .or. config_flags%nested) specified = .true.
3410    ktf=MIN(kte,kde-1)
3411    
3412       i_start = its
3413       i_end   = MIN(ite,ide-1)
3414       j_start = jts
3415       j_end   = MIN(jte,jde-1)
3417       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3418       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-1,jte)
3420 j_loop_v : DO j = j_start, j_end
3421 !     jm1 = max(j-1,1)
3422      jm1 = j-1
3424      DO k=kts,ktf-1
3425        DO i = i_start, i_end
3426          vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k  ,j  )      &
3427                                         +alt(i,k  ,jm1)      &
3428                                         +alt(i,k+1,j  )      &
3429                                         +alt(i,k+1,jm1) ) )  &
3430                              *(field(i,k+1,j)-field(i,k,j)   &
3431                                -v_base(k+1)   +v_base(k)  )
3432        ENDDO
3433      ENDDO
3435      DO i = i_start, i_end
3436        vflux(i,0)=vflux(i,1)
3437      ENDDO
3439      DO i = i_start, i_end
3440        vflux(i,ktf)=0.
3441      ENDDO
3443      DO k=kts,ktf-1
3444        DO i = i_start, i_end 
3445          tendency(i,k,j)=tendency(i,k,j)+                              &
3446                 g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))*  &
3447                               (vflux(i,k)-vflux(i,k-1))
3448        ENDDO
3449      ENDDO
3451  ENDDO j_loop_v
3452    
3453 END SUBROUTINE vertical_diffusion_v
3455 !***************  end new mass coordinate routines
3457 !-------------------------------------------------------------------------------
3459 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp,     &
3460                             ids, ide, jds, jde, kds, kde, &
3461                             ims, ime, jms, jme, kms, kme, &
3462                             its, ite, jts, jte, kts, kte )
3464    IMPLICIT NONE
3465    
3466    ! Input data
3467    
3468    INTEGER ,      INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3469                                    ims, ime, jms, jme, kms, kme, &
3470                                    its, ite, jts, jte, kts, kte 
3471    
3472    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: rfieldb, &
3473                                                                       rfieldp
3475    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT  ) :: rfield
3476    
3477    ! Local indices.
3478    
3479    INTEGER :: i, j, k, itf, jtf, ktf
3480    
3481 !<DESCRIPTION>
3483 !  calculate_full
3484 !  calculates full 3D field from pertubation and base field.
3486 !</DESCRIPTION>
3488    itf=MIN(ite,ide-1)
3489    jtf=MIN(jte,jde-1)
3490    ktf=MIN(kte,kde-1)
3492    DO j=jts,jtf
3493    DO k=kts,ktf
3494    DO i=its,itf
3495       rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3496    ENDDO
3497    ENDDO
3498    ENDDO
3500 END SUBROUTINE calculate_full
3502 !------------------------------------------------------------------------------
3504 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3505                       config_flags,                          &
3506                       msftx, msfty, msfux, msfuy,            &
3507                       msfvx, msfvy,                          &
3508                       f, e, sina, cosa, fzm, fzp,            &
3509                       ids, ide, jds, jde, kds, kde,          &
3510                       ims, ime, jms, jme, kms, kme,          &
3511                       its, ite, jts, jte, kts, kte          )
3513    IMPLICIT NONE
3514    
3515    ! Input data
3516    
3517    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
3519    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3520                                               ims, ime, jms, jme, kms, kme, &
3521                                               its, ite, jts, jte, kts, kte
3523    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3524                                                                 rv_tend, &
3525                                                                 rw_tend
3526    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: ru, &
3527                                                                 rv, &
3528                                                                 rw
3530    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
3531                                                                 msfuy,      &
3532                                                                 msfvx,      &
3533                                                                 msfvy,      &
3534                                                                 msftx,      &
3535                                                                 msfty
3537    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: f,    &
3538                                                                     e,    &
3539                                                                     sina, &
3540                                                                     cosa
3542    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm, &
3543                                                                   fzp
3544    
3545    ! Local indices.
3546    
3547    INTEGER :: i, j , k, ktf
3548    INTEGER :: i_start, i_end, j_start, j_end
3549    
3550    LOGICAL :: specified
3552 !<DESCRIPTION>
3554 !  coriolis calculates the large timestep tendency terms in the 
3555 !  u, v, and w momentum equations arise from the coriolis force.
3557 !</DESCRIPTION>
3559    specified = .false.
3560    if(config_flags%specified .or. config_flags%nested) specified = .true.
3562    ktf=MIN(kte,kde-1)
3564 ! coriolis for u-momentum equation
3566 !  Notes on map scale factor
3567 !  cosa, sina are related to rotating the coordinate frame if desired
3568 !  generally sina=0, cosa=1
3569 !  ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3570 !                                + 2 mu v omega sin(lat)/my
3571 !  Define f=2 omega sin(lat), e=2 omega cos(lat)
3572 !   => terms are: -e mu w / my + f mu v / my
3573 !  rv = mu v / mx ; rw = mu w / my
3574 !   => terms are: -e rw + f rv *mx / my
3576    i_start = its
3577    i_end   = ite
3578    IF ( config_flags%open_xs .or. specified .or. &
3579         config_flags%nested) i_start = MAX(ids+1,its)
3580    IF ( config_flags%open_xe .or. specified .or. &
3581         config_flags%nested) i_end   = MIN(ide-1,ite)
3582       IF ( config_flags%periodic_x ) i_start = its
3583       IF ( config_flags%periodic_x ) i_end = ite
3585    DO j = jts, MIN(jte,jde-1)
3587    DO k=kts,ktf
3588    DO i = i_start, i_end
3589    
3590      ru_tend(i,k,j)=ru_tend(i,k,j) + (msfux(i,j)/msfuy(i,j))*0.5*(f(i,j)+f(i-1,j)) &
3591        *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3592            - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3593        *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3595    ENDDO
3596    ENDDO
3598 ! boundary loops for coriolis not needed for open bdy  (commented out 20100611 JD)
3599 !  IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3601 !    DO k=kts,ktf
3602 !  
3603 !      ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j))   &
3604 !        *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3605 !            - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3606 !        *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3608 !    ENDDO
3610 !  ENDIF
3612 !  IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3614 !    DO k=kts,ktf
3615 !  
3616 !      ru_tend(ite,k,j)=ru_tend(ite,k,j) + (msfux(ite,j)/msfuy(ite,j))*0.5*(f(ite-1,j)+f(ite-1,j)) &
3617 !        *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3618 !            - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3619 !        *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3621 !    ENDDO
3623 !  ENDIF
3625    ENDDO
3627 !  coriolis term for v-momentum equation
3629 !  Notes on map scale factors
3630 !  ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3631 !  Define f=2 omega sin(lat), e=2 omega cos(lat)
3632 !   => terms are: -f mu u / mx
3633 !  ru = mu u / my ; rw = mu w / my
3634 !   => terms are: -f ru *my / mx + ?
3636    j_start = jts
3637    j_end   = jte
3639    IF ( config_flags%open_ys .or. specified .or. &
3640         config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3641    IF ( config_flags%open_ye .or. specified .or. &
3642         config_flags%nested .or. config_flags%polar) j_end   = MIN(jde-1,jte)
3644 ! boundary loops for coriolis not needed for open bdy  (commented out 20100611 JD)
3645 !  IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3647 !    DO k=kts,ktf
3648 !    DO i=its,MIN(ide-1,ite)
3649 !  
3650 !       rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts))    &
3651 !        *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts))   &
3652 !            + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts))   &
3653 !            *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts)) 
3655 !    ENDDO
3656 !    ENDDO
3658 !  ENDIF
3660    DO j=j_start, j_end
3661    DO k=kts,ktf
3662    DO i=its,MIN(ide-1,ite)
3663    
3664       rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*0.5*(f(i,j)+f(i,j-1))    &
3665        *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3666            + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3667            *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j)) 
3669    ENDDO
3670    ENDDO
3671    ENDDO
3674 ! boundary loops for coriolis not needed for open bdy  (commented out 20100611 JD)
3675 !  IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3677 !    DO k=kts,ktf
3678 !    DO i=its,MIN(ide-1,ite)
3679 !  
3680 !       rv_tend(i,k,jte)=rv_tend(i,k,jte) - (msfvy(i,jte)/msfvx(i,jte))*0.5*(f(i,jte-1)+f(i,jte-1))        &
3681 !        *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1))   &
3682 !            + (msfvy(i,jte)/msfvx(i,jte))*0.5*(e(i,jte-1)+e(i,jte-1))*0.5*(sina(i,jte-1)+sina(i,jte-1))   &
3683 !            *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1)) 
3685 !    ENDDO
3686 !    ENDDO
3688 !  ENDIF
3690 ! coriolis term for w-mometum 
3692 ! Notes on map scale factors
3693 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3694 ! Define e=2 omega cos(lat)
3695 !  => terms are: e mu u / my + ???
3696 ! ru = mu u / my ; ru = mu v / mx
3697 !  => terms are: e ru + ???
3699    DO j=jts,MIN(jte, jde-1)
3700    DO k=kts+1,ktf
3701    DO i=its,MIN(ite, ide-1)
3703        rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)*           &
3704           (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3705           +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j)))           &
3706           -(msftx(i,j)/msfty(i,j))*                      &
3707            sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3708           +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3710    ENDDO
3711    ENDDO
3712    ENDDO
3714 END SUBROUTINE coriolis
3716 !------------------------------------------------------------------------------
3718 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3719                                    config_flags,                                &
3720                                    u_base, v_base, z_base,                      &
3721                                    muu, muv, phb, ph,                           &
3722                                    msftx, msfty, msfux, msfuy, msfvx, msfvy,    &
3723                                    f, e, sina, cosa, fzm, fzp,                  &
3724                                    ids, ide, jds, jde, kds, kde,                &
3725                                    ims, ime, jms, jme, kms, kme,                &
3726                                    its, ite, jts, jte, kts, kte                )
3728    IMPLICIT NONE
3729    
3730    ! Input data
3731    
3732    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
3734    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3735                                               ims, ime, jms, jme, kms, kme, &
3736                                               its, ite, jts, jte, kts, kte
3738    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3739                                                                 rv_tend, &
3740                                                                 rw_tend
3741    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: ru_in, &
3742                                                                       rv_in, &
3743                                                                       rw,    &
3744                                                                       ph,    &
3745                                                                       phb
3748    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
3749                                                                 msfuy,      &
3750                                                                 msfvx,      &
3751                                                                 msfvy,      &
3752                                                                 msftx,      &
3753                                                                 msfty
3755    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: f,    &
3756                                                                     e,    &
3757                                                                     sina, &
3758                                                                     cosa
3760    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: muu, &
3761                                                                     muv
3762                                                                     
3764    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm, &
3765                                                                   fzp
3767    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: u_base,  &
3768                                                                   v_base,  &
3769                                                                   z_base
3770    
3771    ! Local storage
3773    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3774                                                       rv
3776    REAL  :: z_at_u, z_at_v, wkp1, wk, wkm1
3778    ! Local indices.
3779    
3780    INTEGER :: i, j , k, ktf
3781    INTEGER :: i_start, i_end, j_start, j_end
3782    
3783    LOGICAL :: specified
3785 !<DESCRIPTION>
3787 !  perturbation_coriolis calculates the large timestep tendency terms in the 
3788 !  u, v, and w momentum equations arise from the coriolis force.  This version
3789 !  subtracts off the horizontal velocities from the initial sounding when
3790 !  computing the forcing terms, hence "perturbation" coriolis.
3792 !</DESCRIPTION>
3794    specified = .false.
3795    if(config_flags%specified .or. config_flags%nested) specified = .true.
3797    ktf=MIN(kte,kde-1)
3799 ! coriolis for u-momentum equation
3801    i_start = its
3802    i_end   = ite
3803    IF ( config_flags%open_xs .or. specified .or. &
3804         config_flags%nested) i_start = MAX(ids+1,its)
3805    IF ( config_flags%open_xe .or. specified .or. &
3806         config_flags%nested) i_end   = MIN(ide-1,ite)
3807       IF ( config_flags%periodic_x ) i_start = its
3808       IF ( config_flags%periodic_x ) i_end = ite
3810 !  compute perturbation mu*v for use in u momentum equation
3812    DO j = jts, MIN(jte,jde-1)+1
3813    DO k=kts+1,ktf-1
3814    DO i = i_start-1, i_end
3815      z_at_v = 0.25*( phb(i,k,j  )+phb(i,k+1,j  )  &
3816                     +phb(i,k,j-1)+phb(i,k+1,j-1)  &
3817                     +ph(i,k,j  )+ph(i,k+1,j  )    &
3818                     +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3819      wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3820      wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3821      wk   = 1.-wkp1-wkm1
3822      rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*(            &
3823                                   wkm1*v_base(k-1)    &
3824                                  +wk  *v_base(k  )    &
3825                                  +wkp1*v_base(k+1)   )
3826    ENDDO
3827    ENDDO
3828    ENDDO
3831 !  pick up top and bottom v 
3833    DO j = jts, MIN(jte,jde-1)+1
3834    DO i = i_start-1, i_end
3836      k = kts
3837      z_at_v = 0.25*( phb(i,k,j  )+phb(i,k+1,j  )  &
3838                     +phb(i,k,j-1)+phb(i,k+1,j-1)  &
3839                     +ph(i,k,j  )+ph(i,k+1,j  )    &
3840                     +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3841      wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3842      wk   = 1.-wkp1
3843      rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*(            &
3844                                  +wk  *v_base(k  )    &
3845                                  +wkp1*v_base(k+1)   )
3847      k = ktf
3848      z_at_v = 0.25*( phb(i,k,j  )+phb(i,k+1,j  )  &
3849                     +phb(i,k,j-1)+phb(i,k+1,j-1)  &
3850                     +ph(i,k,j  )+ph(i,k+1,j  )    &
3851                     +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3852      wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3853      wk   = 1.-wkm1
3854      rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*(            &
3855                                   wkm1*v_base(k-1)    &
3856                                  +wk  *v_base(k  )   )
3858    ENDDO
3859    ENDDO
3861 !  compute coriolis forcing for u
3863 !  Map scale factors: see comments above for Coriolis
3865    DO j = jts, MIN(jte,jde-1)
3867    DO k=kts,ktf
3868      DO i = i_start, i_end
3869        ru_tend(i,k,j)=ru_tend(i,k,j) + (msfux(i,j)/msfuy(i,j))*0.5*(f(i,j)+f(i-1,j)) &
3870          *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3871              - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3872          *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3873      ENDDO
3874    ENDDO
3876 ! boundary loops for coriolis not needed for open bdy  (commented out 20100611 JD)
3877 !  IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3879 !    DO k=kts,ktf
3880 !  
3881 !      ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j))   &
3882 !        *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3883 !            - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3884 !        *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3886 !    ENDDO
3888 !  ENDIF
3890 !  IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3892 !    DO k=kts,ktf
3893 !  
3894 !      ru_tend(ite,k,j)=ru_tend(ite,k,j) + (msfux(ite,j)/msfuy(ite,j))*0.5*(f(ite-1,j)+f(ite-1,j)) &
3895 !        *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3896 !            - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3897 !        *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3899 !    ENDDO
3901 !  ENDIF
3903    ENDDO
3905 !  coriolis term for v-momentum equation
3906 !  Map scale factors: see comments above for Coriolis
3908    j_start = jts
3909    j_end   = jte
3911    IF ( config_flags%open_ys .or. specified .or. &
3912         config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3913    IF ( config_flags%open_ye .or. specified .or. &
3914         config_flags%nested .or. config_flags%polar) j_end   = MIN(jde-1,jte)
3916 !  compute perturbation mu*u for use in v momentum equation
3918    DO j = j_start-1,j_end
3919    DO k=kts+1,ktf-1
3920    DO i = its, MIN(ite,ide-1)+1
3921      z_at_u = 0.25*( phb(i  ,k,j)+phb(i  ,k+1,j)  &
3922                     +phb(i-1,k,j)+phb(i-1,k+1,j)  &
3923                     +ph(i  ,k,j)+ph(i  ,k+1,j)    &
3924                     +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3925      wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3926      wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3927      wk   = 1.-wkp1-wkm1
3928      ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*(            &
3929                                   wkm1*u_base(k-1)    &
3930                                  +wk  *u_base(k  )    &
3931                                  +wkp1*u_base(k+1)   )
3932    ENDDO
3933    ENDDO
3934    ENDDO
3936 !  pick up top and bottom u
3938    DO j = j_start-1,j_end
3939    DO i = its, MIN(ite,ide-1)+1
3941      k = kts
3942      z_at_u = 0.25*( phb(i  ,k,j)+phb(i  ,k+1,j)  &
3943                     +phb(i-1,k,j)+phb(i-1,k+1,j)  &
3944                     +ph(i  ,k,j)+ph(i  ,k+1,j)    &
3945                     +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3946      wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3947      wk   = 1.-wkp1
3948      ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*(            &
3949                                  +wk  *u_base(k  )    &
3950                                  +wkp1*u_base(k+1)   )
3953      k = ktf
3954      z_at_u = 0.25*( phb(i  ,k,j)+phb(i  ,k+1,j)  &
3955                     +phb(i-1,k,j)+phb(i-1,k+1,j)  &
3956                     +ph(i  ,k,j)+ph(i  ,k+1,j)    &
3957                     +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3958      wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3959      wk   = 1.-wkm1
3960      ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*(            &
3961                                   wkm1*u_base(k-1)    &
3962                                  +wk  *u_base(k  )   )
3964    ENDDO
3965    ENDDO
3967 !  compute coriolis forcing for v momentum equation
3968 !  Map scale factors: see comments above for Coriolis
3970 ! boundary loops for coriolis not needed for open bdy  (commented out 20100611 JD)
3971 !  IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3973 !    DO k=kts,ktf
3974 !    DO i=its,MIN(ide-1,ite)
3975 !  
3976 !       rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts))    &
3977 !        *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts))   &
3978 !            + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts))   &
3979 !            *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts)) 
3981 !    ENDDO
3982 !    ENDDO
3984 !  ENDIF
3986    DO j=j_start, j_end
3987    DO k=kts,ktf
3988    DO i=its,MIN(ide-1,ite)
3989    
3990       rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*0.5*(f(i,j)+f(i,j-1))    &
3991        *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3992            + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3993            *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j)) 
3995    ENDDO
3996    ENDDO
3997    ENDDO
4000 ! boundary loops for coriolis not needed for open bdy  (commented out 20100611 JD)
4001 !  IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
4003 !    DO k=kts,ktf
4004 !    DO i=its,MIN(ide-1,ite)
4005 !  
4006 !       rv_tend(i,k,jte)=rv_tend(i,k,jte) - (msfvy(i,jte)/msfvx(i,jte))*0.5*(f(i,jte-1)+f(i,jte-1))        &
4007 !        *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1))   &
4008 !            + (msfvy(i,jte)/msfvx(i,jte))*0.5*(e(i,jte-1)+e(i,jte-1))*0.5*(sina(i,jte-1)+sina(i,jte-1))   &
4009 !            *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1)) 
4011 !    ENDDO
4012 !    ENDDO
4014 !  ENDIF
4016 ! coriolis term for w-mometum 
4017 !  Map scale factors: see comments above for Coriolis
4019    DO j=jts,MIN(jte, jde-1)
4020    DO k=kts+1,ktf
4021    DO i=its,MIN(ite, ide-1)
4023        rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)*           &
4024           (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
4025           +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j)))           &
4026           -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
4027           +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
4029    ENDDO
4030    ENDDO
4031    ENDDO
4033 END SUBROUTINE perturbation_coriolis
4035 !------------------------------------------------------------------------------
4037 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
4038                         config_flags,                                       &
4039                         msfux, msfuy, msfvx, msfvy, msftx, msfty,       &
4040                         xlat, fzm, fzp, rdx, rdy,                       &
4041                         ids, ide, jds, jde, kds, kde,                   &
4042                         ims, ime, jms, jme, kms, kme,                   &
4043                         its, ite, jts, jte, kts, kte                   )
4046    IMPLICIT NONE
4047    
4048    ! Input data
4050    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
4052    INTEGER ,                  INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4053                                                ims, ime, jms, jme, kms, kme, &
4054                                                its, ite, jts, jte, kts, kte
4055    
4056    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                     &
4057                                                INTENT(INOUT) :: ru_tend, &
4058                                                                 rv_tend, &
4059                                                                 rw_tend
4061    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                     &
4062                                                INTENT(IN   ) :: ru,      &
4063                                                                 rv,      &
4064                                                                 rw,      &
4065                                                                 u,       &
4066                                                                 v,       &
4067                                                                 w
4069    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,    &
4070                                                                 msfuy,    &
4071                                                                 msfvx,    &
4072                                                                 msfvy,    &
4073                                                                 msftx,    &
4074                                                                 msfty,    &
4075                                                                 xlat
4077    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,     &
4078                                                                 fzp
4080    REAL ,                                      INTENT(IN   ) :: rdx,     &
4081                                                                 rdy
4082    
4083    ! Local data
4084    
4085 !   INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
4086    INTEGER :: i, j, k, itf, jtf, ktf
4087    INTEGER :: i_start, i_end, j_start, j_end
4088 !   INTEGER :: irmin, irmax, jrmin, jrmax
4090    REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
4092    LOGICAL :: specified
4094 !<DESCRIPTION>
4096 !  curvature calculates the large timestep tendency terms in the 
4097 !  u, v, and w momentum equations arise from the curvature terms.  
4099 !</DESCRIPTION>
4101    specified = .false.
4102    if(config_flags%specified .or. config_flags%nested) specified = .true.
4104       itf=MIN(ite,ide-1)
4105       jtf=MIN(jte,jde-1)
4106       ktf=MIN(kte,kde-1)
4108 !   irmin = ims
4109 !   irmax = ime
4110 !   jrmin = jms
4111 !   jrmax = jme
4112 !   IF ( config_flags%open_xs ) irmin = ids
4113 !   IF ( config_flags%open_xe ) irmax = ide-1
4114 !   IF ( config_flags%open_ys ) jrmin = jds
4115 !   IF ( config_flags%open_ye ) jrmax = jde-1
4116    
4117 ! Define v cross grad m at scalar points - vxgm(i,j)
4119    i_start = its-1
4120    i_end   = ite
4121    j_start = jts-1
4122    j_end   = jte
4124    IF ( ( config_flags%open_xs .or. specified .or. &
4125         config_flags%nested) .and. (its == ids) ) i_start = its
4126    IF ( ( config_flags%open_xe .or. specified .or. &
4127         config_flags%nested) .and. (ite == ide) ) i_end   = ite-1
4128    IF ( ( config_flags%open_ys .or. specified .or. &
4129         config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
4130    IF ( ( config_flags%open_ye .or. specified .or. &
4131         config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end   = jte-1
4132       IF ( config_flags%periodic_x ) i_start = its-1
4133       IF ( config_flags%periodic_x ) i_end = ite
4135    DO j=j_start, j_end
4136    DO k=kts,ktf
4137    DO i=i_start, i_end
4138 !     Map scale factor notes:
4139 !     msf...y is constant everywhere for cylindrical map projection
4140 !     msf...x varies with y only
4141 !     But we know that this is not = 0 for cylindrical,
4142 !     therefore use msfvX in 1st line
4143 !     which => by symmetry use msfuY in 2nd line - ???  
4144       vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
4145                   0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
4146    ENDDO
4147    ENDDO
4148    ENDDO
4150 !  Pick up the boundary rows for open (radiation) lateral b.c.
4151 !  Rather crude at present, we are assuming there is no
4152 !    variation in this term at the boundary.
4154    IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4155         config_flags%nested) .and. (its == ids) ) THEN
4157      DO j = jts, jte-1
4158      DO k = kts, ktf
4159        vxgm(its-1,k,j) =  vxgm(its,k,j)
4160      ENDDO
4161      ENDDO
4163    ENDIF
4165    IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4166         config_flags%nested) .and. (ite == ide) ) THEN
4168      DO j = jts, jte-1
4169      DO k = kts, ktf
4170        vxgm(ite,k,j) =  vxgm(ite-1,k,j)
4171      ENDDO
4172      ENDDO
4174    ENDIF
4176 !  Polar boundary condition:
4177 !  The following change is needed in case one tries using the vxgm route with
4178 !  polar B.C.'s in the future, but not needed if 'tan' used
4179    IF ( ( config_flags%open_ys .or. specified .or. &
4180         config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
4182      DO k = kts, ktf
4183      DO i = its-1, ite
4184        vxgm(i,k,jts-1) =  vxgm(i,k,jts)
4185      ENDDO
4186      ENDDO
4188    ENDIF
4190 !  Polar boundary condition:
4191 !  The following change is needed in case one tries using the vxgm route with
4192 !  polar B.C.'s in the future, but not needed if 'tan' used
4193    IF ( ( config_flags%open_ye .or. specified .or. &
4194         config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
4196      DO k = kts, ktf
4197      DO i = its-1, ite
4198        vxgm(i,k,jte) =  vxgm(i,k,jte-1)
4199      ENDDO
4200      ENDDO
4202    ENDIF
4204 !  curvature term for u momentum eqn.
4206 !  Map scale factor notes:
4207 !  ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4208 !                                               - mu u w /(a my)
4209 !  ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4210 !   => terms are:
4211 !  (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
4212 !  ru v tan(lat) / a - u rw / a
4213 !  xlat defined with end points half grid space from pole,
4214 !  hence are on u latitude points
4216    i_start = its
4217    IF ( config_flags%open_xs .or. specified .or. &
4218         config_flags%nested) i_start = MAX ( ids+1 , its )
4219    IF ( config_flags%open_xe .or. specified .or. &
4220         config_flags%nested) i_end   = MIN ( ide-1 , ite )
4221       IF ( config_flags%periodic_x ) i_start = its
4222       IF ( config_flags%periodic_x ) i_end = ite
4224 !  Polar boundary condition
4225    IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4227       DO j=jts,MIN(jde-1,jte)
4228       DO k=kts,ktf
4229       DO i=i_start,i_end
4231             ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius*                 ( &
4232                         (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4233                                     rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
4234                         - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4235       ENDDO
4236       ENDDO
4237       ENDDO
4239    ELSE  ! normal code
4242       DO j=jts,MIN(jde-1,jte)
4243       DO k=kts,ktf
4244       DO i=i_start,i_end
4246          ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4247                  *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4248                   - u(i,k,j)*reradius &
4249                  *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4251       ENDDO
4252       ENDDO
4253       ENDDO
4255    END IF
4257 !  curvature term for v momentum eqn.
4259 !  Map scale factor notes
4260 !  ADT eqn 45, RHS terms 4 and 5, in cylindrical:  - mu u*u tan(lat)/(a mx)
4261 !                                               - mu v w /(a mx)
4262 !  ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4263 !  terms are:
4264 !  - (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a 
4265 !  = - [my/(mx*a)]*[u ru tan(lat) + v rw]
4266 !  - (1/a)*[(my/mx)*u ru tan(lat) + w rv]
4267 !  xlat defined with end points half grid space from pole, hence are on
4268 !  u latitude points => av here
4270 !  in original wrf, there was a sign error for the rw contribution
4272    j_start = jts
4273    IF ( config_flags%open_ys .or. specified .or. &
4274         config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4275    IF ( config_flags%open_ye .or. specified .or. &
4276         config_flags%nested .or. config_flags%polar) j_end   = MIN ( jde-1 , jte )
4278    IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4280       DO j=j_start,j_end
4281       DO k=kts,ktf
4282       DO i=its,MIN(ite,ide-1)
4283             rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius*   (  &
4284                         0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))*     &
4285                         tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)*                &
4286                         0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1))  &
4287                        + v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+              &
4288                                                       rw(i,k+1,j)+rw(i,k,j))    )
4289       ENDDO
4290       ENDDO
4291       ENDDO
4293    ELSE  ! normal code
4295       DO j=j_start,j_end
4296       DO k=kts,ktf
4297       DO i=its,MIN(ite,ide-1)
4299          rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4300                  *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4301                        - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius       &
4302                  *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4304       ENDDO
4305       ENDDO
4306       ENDDO
4308    END IF
4310 !  curvature term for vertical momentum eqn.
4312 !  Notes on map scale factors:
4313 !  ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4314 !  ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4315 !  terms are: u ru / a + (mx/my)v rv / a
4317    DO j=jts,MIN(jte,jde-1)
4318    DO k=MAX(2,kts),ktf
4319    DO i=its,MIN(ite,ide-1)
4321       rw_tend(i,k,j)=rw_tend(i,k,j) + reradius*                              &
4322     (0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j))+fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
4323     *0.5*(fzm(k)*( u(i,k,j) +u(i+1,k,j))+fzp(k)*( u(i,k-1,j) +u(i+1,k-1,j)))     &
4324     +(msftx(i,j)/msfty(i,j))*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1))+fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))) &
4325     *0.5*(fzm(k)*( v(i,k,j) +v(i,k,j+1))+fzp(k)*( v(i,k-1,j) +v(i,k-1,j+1))))
4327    ENDDO
4328    ENDDO
4329    ENDDO
4331 END SUBROUTINE curvature
4333 !------------------------------------------------------------------------------
4335 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4336                       fzm, fzp,                          &
4337                       ids, ide, jds, jde, kds, kde,      &
4338                       ims, ime, jms, jme, kms, kme,      &
4339                       its, ite, jts, jte, kts, kte      )
4341    IMPLICIT NONE
4343    ! Input data
4345    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
4347    INTEGER ,                                   INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4348                                                                 ims, ime, jms, jme, kms, kme, &
4349                                                                 its, ite, jts, jte, kts, kte
4351    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
4353    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: rfield
4355    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: rr
4356    
4357    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT  ) :: field
4358    
4359    REAL , DIMENSION( kms:kme ) , INTENT(IN   ) :: fzm, fzp
4360    
4361    ! Local data
4362    
4363    INTEGER :: i, j, k, itf, jtf, ktf
4364    
4365 !<DESCRIPTION>
4367 !  decouple decouples a variable from the column dry-air mass.
4369 !</DESCRIPTION>
4371    ktf=MIN(kte,kde-1)
4372    
4373    IF (name .EQ. 'u')THEN
4374       itf=ite
4375       jtf=MIN(jte,jde-1)
4377       DO j=jts,jtf
4378       DO k=kts,ktf
4379       DO i=its,itf
4380          field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4381       ENDDO
4382       ENDDO
4383       ENDDO
4385    ELSE IF (name .EQ. 'v')THEN
4386       itf=MIN(ite,ide-1)
4387       jtf=jte
4389       DO j=jts,jtf
4390       DO k=kts,ktf
4391         DO i=its,itf
4392              field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4393         ENDDO
4394       ENDDO
4395       ENDDO
4397    ELSE IF (name .EQ. 'w')THEN
4398       itf=MIN(ite,ide-1)
4399       jtf=MIN(jte,jde-1)
4400       DO j=jts,jtf
4401       DO k=kts+1,ktf
4402       DO i=its,itf
4403          field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4404       ENDDO
4405       ENDDO
4406       ENDDO
4408       DO j=jts,jtf
4409       DO i=its,itf
4410         field(i,kte,j) = 0.
4411       ENDDO
4412       ENDDO
4414    ELSE 
4415       itf=MIN(ite,ide-1)
4416       jtf=MIN(jte,jde-1)
4417    ! For theta we will decouple tb and tp and add them to give t afterwards
4418       DO j=jts,jtf
4419       DO k=kts,ktf
4420       DO i=its,itf
4421          field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4422       ENDDO
4423       ENDDO
4424       ENDDO
4425    
4426    ENDIF
4428 END SUBROUTINE decouple
4430 !-------------------------------------------------------------------------------
4433 SUBROUTINE zero_tend ( tendency,                     &
4434                        ids, ide, jds, jde, kds, kde, &
4435                        ims, ime, jms, jme, kms, kme, &
4436                        its, ite, jts, jte, kts, kte )
4439    IMPLICIT NONE
4440    
4441    ! Input data
4442    
4443    INTEGER ,                                   INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4444                                                                 ims, ime, jms, jme, kms, kme, &
4445                                                                 its, ite, jts, jte, kts, kte
4447    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4449    ! Local data
4450    
4451    INTEGER :: i, j, k, itf, jtf, ktf
4453 !<DESCRIPTION>
4455 !  zero_tend sets the input tendency array to zero.
4457 !</DESCRIPTION>
4459       DO j = jts, jte
4460       DO k = kts, kte
4461       DO i = its, ite
4462         tendency(i,k,j) = 0.
4463       ENDDO
4464       ENDDO
4465       ENDDO
4467       END SUBROUTINE zero_tend
4469 !-------------------------------------------------------------------------------
4470 ! Sets the an array on the polar v point(s) to zero
4471 SUBROUTINE zero_pole ( field,                        &
4472                        ids, ide, jds, jde, kds, kde, &
4473                        ims, ime, jms, jme, kms, kme, &
4474                        its, ite, jts, jte, kts, kte )
4477   IMPLICIT NONE
4479   ! Input data
4480    
4481   INTEGER , INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4482                              ims, ime, jms, jme, kms, kme, &
4483                              its, ite, jts, jte, kts, kte
4485   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4487   ! Local data
4489   INTEGER :: i, k
4491   IF (jts == jds) THEN
4492      DO k = kts, kte
4493      DO i = its-1, ite+1
4494         field(i,k,jts) = 0.
4495      END DO
4496      END DO
4497   END IF
4498   IF (jte == jde) THEN
4499      DO k = kts, kte
4500      DO i = its-1, ite+1
4501         field(i,k,jte) = 0.
4502      END DO
4503      END DO
4504   END IF
4506 END SUBROUTINE zero_pole
4508 !-------------------------------------------------------------------------------
4509 ! Sets the an array on the polar v point(s)
4510 SUBROUTINE pole_point_bc ( field,                        &
4511                        ids, ide, jds, jde, kds, kde, &
4512                        ims, ime, jms, jme, kms, kme, &
4513                        its, ite, jts, jte, kts, kte )
4516   IMPLICIT NONE
4518   ! Input data
4519    
4520   INTEGER , INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4521                              ims, ime, jms, jme, kms, kme, &
4522                              its, ite, jts, jte, kts, kte
4524   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4526   ! Local data
4528   INTEGER :: i, k
4530   IF (jts == jds) THEN
4531      DO k = kts, kte
4532      DO i = its, ite
4533 !        field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4534         field(i,k,jts) = field(i,k,jts+1)
4535      END DO
4536      END DO
4537   END IF
4538   IF (jte == jde) THEN
4539      DO k = kts, kte
4540      DO i = its, ite
4541 !        field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4542         field(i,k,jte) = field(i,k,jte-1)
4543      END DO
4544      END DO
4545   END IF
4547 END SUBROUTINE pole_point_bc
4549 !======================================================================
4550 !   physics prep routines
4551 !======================================================================
4553    SUBROUTINE phy_prep ( config_flags,                                &  ! input
4554                          mu, muu, muv, u, v, p, pb, alt, ph,          &  ! input
4555                          phb, t, tsk, moist, n_moist,                 &  ! input
4556                          rho, th_phy, p_phy , pi_phy ,                &  ! output
4557                          u_phy, v_phy, p8w, t_phy, t8w,               &  ! output
4558                          z, z_at_w, dz8w,                             &  ! output
4559                          p_hyd, p_hyd_w,                              &  ! output
4560                          fzm, fzp, znw, p_top,                        &  ! params
4561                          RTHRATEN,                                    &
4562                          RTHBLTEN, RUBLTEN, RVBLTEN,                  &
4563                          RQVBLTEN, RQCBLTEN, RQIBLTEN,                &
4564                          RTHCUTEN, RQVCUTEN, RQCCUTEN,                &
4565                          RQRCUTEN, RQICUTEN, RQSCUTEN,                &
4566                          RTHFTEN,  RQVFTEN,                           &
4567                          RUCUTEN,  RVCUTEN,                           &  ! Tiedtke and NSAS
4568                          RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN,            &
4569                          RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN,           &
4570                          ids, ide, jds, jde, kds, kde,                &
4571                          ims, ime, jms, jme, kms, kme,                &
4572                          its, ite, jts, jte, kts, kte                )
4573 !----------------------------------------------------------------------
4574    IMPLICIT NONE
4575 !----------------------------------------------------------------------
4577    TYPE(grid_config_rec_type) ,     INTENT(IN   ) :: config_flags
4579    INTEGER ,        INTENT(IN   ) ::   ids, ide, jds, jde, kds, kde, &
4580                                        ims, ime, jms, jme, kms, kme, &
4581                                        its, ite, jts, jte, kts, kte
4582    INTEGER ,          INTENT(IN   ) :: n_moist
4584    REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4587    REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN   )   ::     TSK, mu, muu, muv
4589    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
4590           INTENT(  OUT)                                  ::   u_phy, &
4591                                                               v_phy, &
4592                                                              pi_phy, &
4593                                                               p_phy, &
4594                                                                 p8w, &
4595                                                               t_phy, &
4596                                                              th_phy, &
4597                                                                 t8w, &
4598                                                                 rho, &
4599                                                                   z, &
4600                                                                dz8w, &
4601                                                               z_at_w 
4603    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
4604           INTENT(  OUT)                                  ::   p_hyd, &
4605                                                               p_hyd_w
4607    REAL , INTENT(IN   )                                  ::   p_top
4609    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
4610           INTENT(IN   )                                  ::      pb, &
4611                                                                   p, &
4612                                                                   u, &
4613                                                                   v, &
4614                                                                 alt, &
4615                                                                  ph, &
4616                                                                 phb, &
4617                                                                   t
4620    REAL , DIMENSION( kms:kme ) ,           INTENT(IN   ) ::     fzm,   &
4621                                                                 fzp
4623    REAL , DIMENSION( kms:kme ) ,           INTENT(IN   ) ::     znw
4625    REAL,  DIMENSION( ims:ime , kms:kme, jms:jme ),                   &
4626           INTENT(INOUT)   ::                               RTHRATEN  
4628    REAL,  DIMENSION( ims:ime , kms:kme, jms:jme ),                   &
4629           INTENT(INOUT)   ::                               RTHCUTEN, &
4630                                                            RQVCUTEN, &
4631                                                            RQCCUTEN, &
4632                                                            RQRCUTEN, &
4633                                                            RQICUTEN, &
4634                                                            RQSCUTEN
4636    REAL,  DIMENSION( ims:ime , kms:kme, jms:jme ),                   &
4637           INTENT(INOUT)   ::                                RUCUTEN, &
4638                                                             RVCUTEN
4640    REAL,  DIMENSION( ims:ime, kms:kme, jms:jme )                   , &
4641           INTENT(INOUT)   ::                                RUBLTEN, &
4642                                                             RVBLTEN, &
4643                                                            RTHBLTEN, &
4644                                                            RQVBLTEN, &
4645                                                            RQCBLTEN, &
4646                                                            RQIBLTEN
4648    REAL,  DIMENSION( ims:ime, kms:kme, jms:jme )                   , &
4649           INTENT(INOUT)   ::                                RTHFTEN, &
4650                                                             RQVFTEN
4652    REAL,  DIMENSION( ims:ime, kms:kme, jms:jme )                   , &
4653           INTENT(INOUT)   ::                                RUNDGDTEN, &
4654                                                             RVNDGDTEN, &
4655                                                            RTHNDGDTEN, &
4656                                                            RPHNDGDTEN, &
4657                                                            RQVNDGDTEN
4659    REAL,  DIMENSION( ims:ime, jms:jme )                            , &
4660           INTENT(INOUT)   ::                               RMUNDGDTEN
4662    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4663    INTEGER :: i, j, k
4664    REAL    :: w1, w2, z0, z1, z2
4665    REAL    :: e_vapor
4667 !-----------------------------------------------------------------------
4669 !<DESCRIPTION>
4671 !  phys_prep calculates a number of diagnostic quantities needed by
4672 !  the physics routines.  It also decouples the physics tendencies from
4673 !  the column dry-air mass (the physics routines expect to see/update the
4674 !  uncoupled tendencies).
4676 !</DESCRIPTION>
4678 !  set up loop bounds for this grid's boundary conditions
4680     i_start = its
4681     i_end   = min( ite,ide-1 )
4682     j_start = jts
4683     j_end   = min( jte,jde-1 )
4685     k_start = kts
4686     k_end = min( kte, kde-1 )
4688 !  compute thermodynamics and velocities at pressure points (or half levels)
4690     do j = j_start,j_end
4691     do k = k_start, k_end
4692     do i = i_start, i_end
4694       th_phy(i,k,j) = t(i,k,j) + t0
4695       p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4696       pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4697       t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4698       rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4699       u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4700       v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4702     enddo
4703     enddo
4704     enddo
4706 !  compute z at w points
4708     do j = j_start,j_end
4709     do k = k_start, kte
4710     do i = i_start, i_end
4711       z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4712     enddo
4713     enddo
4714     enddo
4716     do j = j_start,j_end
4717     do k = k_start, kte-1
4718     do i = i_start, i_end
4719       dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4720     enddo
4721     enddo
4722     enddo
4724     do j = j_start,j_end
4725     do i = i_start, i_end
4726       dz8w(i,kte,j) = 0.
4727     enddo
4728     enddo
4730 !  compute z at p points or half levels (average of z at full levels)
4732     do j = j_start,j_end
4733     do k = k_start, k_end
4734     do i = i_start, i_end
4735       z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4736     enddo
4737     enddo
4738     enddo
4740 !  interp t and p to full levels
4742     do j = j_start,j_end
4743     do k = 2, k_end
4744     do i = i_start, i_end
4745       p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4746       t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4747     enddo
4748     enddo
4749     enddo
4751 !  extrapolate p and t to surface and top.
4752 !  we'll use an extrapolation in z for now
4754     do j = j_start,j_end
4755     do i = i_start, i_end
4757 ! bottom
4759       z0 = z_at_w(i,1,j)
4760       z1 = z(i,1,j)
4761       z2 = z(i,2,j)
4762       w1 = (z0 - z2)/(z1 - z2)
4763       w2 = 1. - w1
4764       p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4765       t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4767 ! top
4769       z0 = z_at_w(i,kte,j)
4770       z1 = z(i,k_end,j)
4771       z2 = z(i,k_end-1,j)
4772       w1 = (z0 - z2)/(z1 - z2)
4773       w2 = 1. - w1
4775 !      p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4776 !!!  bug fix      extrapolate ln(p) so p is positive definite
4777       p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4778       t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4780     enddo
4781     enddo
4783 ! calculate hydrostatic pressure at both full and half levels
4784 ! first, full level p: assuming dry over model top
4786     do j = j_start,j_end
4787     do i = i_start, i_end
4788        p_hyd_w(i,kte,j) = p_top
4789     enddo
4790     enddo
4792     e_vapor = 0.
4793     do j = j_start,j_end
4794     do k = kte-1, k_start, -1
4795     do i = i_start, i_end
4796 !      e_vapor = 1./alt(i,k,j)*moist(i,k,j,P_QV)*g*dz8w(i,k,j) 
4797 !      p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+e_vapor
4798        p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))*g*dz8w(i,k,j)
4799     enddo
4800     enddo
4801     enddo
4803 ! now calculate hydrostatic pressure at half levels
4805     do j = j_start,j_end
4806     do k = k_start, k_end
4807     do i = i_start, i_end
4808        p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
4809     enddo
4810     enddo
4811     enddo
4813 ! decouple all physics tendencies
4815    IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4817       DO J=j_start,j_end
4818       DO K=k_start,k_end
4819       DO I=i_start,i_end
4820          RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4821       ENDDO
4822       ENDDO
4823       ENDDO
4825    ENDIF
4827    IF (config_flags%cu_physics .gt. 0) THEN
4829       DO J=j_start,j_end
4830       DO I=i_start,i_end
4831       DO K=k_start,k_end
4832          RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4833       ENDDO
4834       ENDDO
4835       ENDDO
4837       IF (config_flags%cu_physics == TIEDTKESCHEME .or.        &
4838           config_flags%cu_physics == NSASSCHEME ) THEN
4839         DO J=j_start,j_end
4840         DO I=i_start,i_end
4841         DO K=k_start,k_end
4842            RUCUTEN(I,K,J)=RUCUTEN(I,K,J)/mu(I,J)
4843            RVCUTEN(I,K,J)=RVCUTEN(I,K,J)/mu(I,J)
4844         ENDDO
4845         ENDDO
4846         ENDDO
4847       ENDIF
4849       IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4850          DO J=j_start,j_end
4851          DO I=i_start,i_end
4852          DO K=k_start,k_end
4853             RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4854          ENDDO
4855          ENDDO
4856          ENDDO
4857       ENDIF
4859       IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4860          DO J=j_start,j_end
4861          DO I=i_start,i_end
4862          DO K=k_start,k_end
4863             RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4864          ENDDO
4865          ENDDO
4866          ENDDO
4867       ENDIF
4869       IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4870          DO J=j_start,j_end
4871          DO I=i_start,i_end
4872          DO K=k_start,k_end
4873             RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4874          ENDDO
4875          ENDDO
4876          ENDDO
4877       ENDIF
4879       IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4880          DO J=j_start,j_end
4881          DO I=i_start,i_end
4882          DO K=k_start,k_end
4883             RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4884          ENDDO
4885          ENDDO
4886          ENDDO
4887       ENDIF
4889       IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4890          DO J=j_start,j_end
4891          DO I=i_start,i_end
4892          DO K=k_start,k_end
4893             RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4894          ENDDO
4895          ENDDO
4896          ENDDO
4897       ENDIF
4899    ENDIF
4901    IF (config_flags%bl_pbl_physics .gt. 0) THEN
4903       DO J=j_start,j_end
4904       DO K=k_start,k_end
4905       DO I=i_start,i_end
4906          RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
4907          RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
4908          RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
4909       ENDDO
4910       ENDDO
4911       ENDDO
4913       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4914          DO J=j_start,j_end
4915          DO K=k_start,k_end
4916          DO I=i_start,i_end
4917             RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
4918          ENDDO
4919          ENDDO
4920          ENDDO
4921       ENDIF
4923       IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
4924          DO J=j_start,j_end
4925          DO K=k_start,k_end
4926          DO I=i_start,i_end
4927            RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
4928          ENDDO
4929          ENDDO
4930          ENDDO
4931       ENDIF
4933       IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
4934          DO J=j_start,j_end
4935          DO K=k_start,k_end
4936          DO I=i_start,i_end
4937             RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
4938          ENDDO
4939          ENDDO
4940          ENDDO
4941       ENDIF
4943     ENDIF
4945 !  decouple advective forcing required by Grell-Devenyi scheme
4947    if(( config_flags%cu_physics == GDSCHEME ) .OR.    &
4948       ( config_flags%cu_physics == G3SCHEME ) .OR.    &
4949       ( config_flags%cu_physics == KFETASCHEME ) .OR. &
4950       ( config_flags%cu_physics == TIEDTKESCHEME ) ) then  ! Tiedtke ZCX&YQW
4952       DO J=j_start,j_end
4953       DO I=i_start,i_end
4954          DO K=k_start,k_end
4955             RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
4956          ENDDO
4957       ENDDO
4958       ENDDO
4960       IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4961          DO J=j_start,j_end
4962          DO I=i_start,i_end
4963             DO K=k_start,k_end
4964                RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
4965             ENDDO
4966          ENDDO
4967          ENDDO
4968       ENDIF
4970    END IF
4972 ! fdda
4973 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
4974 !   so only decouple those
4976    IF (config_flags%grid_fdda .gt. 0) THEN
4978       i_startu=MAX(its,ids+1)
4979       j_startv=MAX(jts,jds+1)
4981       DO J=j_start,j_end
4982       DO K=k_start,k_end
4983       DO I=i_startu,i_end
4984          RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
4985       ENDDO
4986       ENDDO
4987       ENDDO
4988       DO J=j_startv,j_end
4989       DO K=k_start,k_end
4990       DO I=i_start,i_end
4991          RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
4992       ENDDO
4993       ENDDO
4994       ENDDO
4995       DO J=j_start,j_end
4996       DO K=k_start,k_end
4997       DO I=i_start,i_end
4998          RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
4999 !        RMUNDGDTEN(I,J) - no coupling
5000       ENDDO
5001       ENDDO
5002       ENDDO
5004       IF (config_flags%grid_fdda .EQ. 2) THEN
5005       DO J=j_start,j_end
5006       DO K=k_start,kte
5007       DO I=i_start,i_end
5008          RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/mu(I,J)
5009       ENDDO
5010       ENDDO
5011       ENDDO
5013       ELSE IF (config_flags%grid_fdda .EQ. 1) THEN
5014       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5015          DO J=j_start,j_end
5016          DO K=k_start,k_end
5017          DO I=i_start,i_end
5018             RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
5019          ENDDO
5020          ENDDO
5021          ENDDO
5022       ENDIF
5023       ENDIF
5025     ENDIF
5027 END SUBROUTINE phy_prep
5029 !------------------------------------------------------------
5031    SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
5032                                      p, p8w, p0, pb, ph, phb,        &
5033                                      th_phy, pii, pf,                &
5034                                      z, z_at_w, dz8w,                &
5035                                      dt,h_diabatic,                  &
5036                                      config_flags,fzm, fzp,          &
5037                                      ids,ide, jds,jde, kds,kde,      &
5038                                      ims,ime, jms,jme, kms,kme,      &
5039                                      its,ite, jts,jte, kts,kte      )
5041    IMPLICIT NONE
5043 ! Here we construct full fields
5044 ! needed by the microphysics
5046    TYPE(grid_config_rec_type),    INTENT(IN   )    :: config_flags
5048    INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
5049    INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
5050    INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
5052    REAL, INTENT(IN   )  ::  dt
5054    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
5055          INTENT(IN   ) ::                           al,  &
5056                                                     alb, &
5057                                                     p,   &
5058                                                     pb,  &
5059                                                     ph,  &
5060                                                     phb
5063    REAL , DIMENSION( kms:kme ) ,           INTENT(IN   ) ::   fzm, &
5064                                                               fzp
5066    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),       &
5067          INTENT(  OUT) ::                         rho,  &
5068                                                th_phy,  &
5069                                                   pii,  &
5070                                                   pf,   &
5071                                                     z,  &
5072                                                z_at_w,  &
5073                                                  dz8w,  &
5074                                                   p8w
5076    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),       &
5077          INTENT(INOUT) ::                         h_diabatic
5079    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
5080          INTENT(INOUT) ::                         t_new, &
5081                                                   t_old
5083    REAL, INTENT(IN   ) :: t0, p0
5084    REAL                :: z0,z1,z2,w1,w2
5086    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5087    INTEGER :: i, j, k
5089 !--------------------------------------------------------------------
5091 !<DESCRIPTION>
5093 !  moist_phys_prep_em calculates a number of diagnostic quantities needed by
5094 !  the microphysics routines.
5096 !</DESCRIPTION>
5098 !  set up loop bounds for this grid's boundary conditions
5100     i_start = its    
5101     i_end   = min( ite,ide-1 )
5102     j_start = jts    
5103     j_end   = min( jte,jde-1 )
5105     k_start = kts
5106     k_end = min( kte, kde-1 )
5108      DO j = j_start, j_end
5109      DO k = k_start, kte
5110      DO i = i_start, i_end
5111        z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
5112      ENDDO
5113      ENDDO
5114      ENDDO
5116     do j = j_start,j_end
5117     do k = k_start, kte-1
5118     do i = i_start, i_end
5119       dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
5120     enddo
5121     enddo
5122     enddo
5124     do j = j_start,j_end
5125     do i = i_start, i_end
5126       dz8w(i,kte,j) = 0.
5127     enddo
5128     enddo
5131            !  compute full pii, rho, and z at the new time-level
5132            !  (needed for physics).
5133            !  convert perturbation theta to full theta (th_phy)
5134            !  use h_diabatic to temporarily save pre-microphysics full theta
5136      DO j = j_start, j_end
5137      DO k = k_start, k_end
5138      DO i = i_start, i_end
5140 #ifdef REVERT
5141        t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
5142 #endif
5143        th_phy(i,k,j) = t_new(i,k,j) + t0
5144        h_diabatic(i,k,j) = th_phy(i,k,j)
5145        rho(i,k,j)  = 1./(al(i,k,j)+alb(i,k,j))
5146        pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
5147        z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
5148        pf(i,k,j) = p(i,k,j)+pb(i,k,j)
5150      ENDDO
5151      ENDDO
5152      ENDDO
5154 !  interp t and p at w points
5156     do j = j_start,j_end
5157     do k = 2, k_end
5158     do i = i_start, i_end
5159       p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
5160     enddo
5161     enddo
5162     enddo
5164 !  extrapolate p and t to surface and top.
5165 !  we'll use an extrapolation in z for now
5167     do j = j_start,j_end
5168     do i = i_start, i_end
5170 ! bottom
5172       z0 = z_at_w(i,1,j)
5173       z1 = z(i,1,j)
5174       z2 = z(i,2,j)
5175       w1 = (z0 - z2)/(z1 - z2)
5176       w2 = 1. - w1
5177       p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
5179 ! top
5181       z0 = z_at_w(i,kte,j)
5182       z1 = z(i,k_end,j)
5183       z2 = z(i,k_end-1,j)
5184       w1 = (z0 - z2)/(z1 - z2)
5185       w2 = 1. - w1
5186 !      p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
5187       p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
5189     enddo
5190     enddo
5192    END SUBROUTINE moist_physics_prep_em
5194 !------------------------------------------------------------------------------
5196    SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut,     &
5197                                        th_phy, h_diabatic, dt,    &
5198                                        config_flags,              &
5199 #if ( WRF_DFI_RADAR == 1 )
5200                                        dfi_tten_rad,dfi_stage,    &
5201 #endif
5202                                        ids,ide, jds,jde, kds,kde, &
5203                                        ims,ime, jms,jme, kms,kme, &
5204                                        its,ite, jts,jte, kts,kte )
5206    IMPLICIT NONE
5208 ! Here we construct full fields
5209 ! needed by the microphysics
5211    TYPE(grid_config_rec_type),    INTENT(IN   )    :: config_flags
5213    INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
5214    INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
5215    INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
5217    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
5218          INTENT(INOUT) ::                         t_new, &
5219                                                   t_old, &
5220                                                  th_phy, &
5221                                                   h_diabatic
5222 #if ( WRF_DFI_RADAR == 1 )
5223    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
5224          INTENT(IN), OPTIONAL ::               dfi_tten_rad
5225    INTEGER,      INTENT(IN   ) ,OPTIONAL   :: dfi_stage
5226    REAL :: dfi_tten_max, old_max
5227 #endif
5229    REAL mpten, mptenmax, mptenmin
5231    REAL, DIMENSION( ims:ime , jms:jme ),  INTENT(INOUT) ::  mut
5234    REAL, INTENT(IN   ) :: t0, dt
5236    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5237    INTEGER :: i, j, k, imax, jmax, imin, jmin
5239 !--------------------------------------------------------------------
5241 !<DESCRIPTION>
5243 !  moist_phys_finish_em resets theta to its perturbation value and
5244 !  computes and stores the microphysics diabatic heating term.
5246 !</DESCRIPTION>
5248 !  set up loop bounds for this grid's boundary conditions
5251     i_start = its    
5252     i_end   = min( ite,ide-1 )
5253     j_start = jts    
5254     j_end   = min( jte,jde-1 )
5255 !      i_start=max(its,ids+4)
5256 !      i_end=min(ite,ide-5)
5257 !      j_start=max(jts,jds+4)
5258 !      j_end=min(jte,jde-5)
5260     k_start = kts
5261     k_end = min( kte, kde-1 )
5263 #if ( WRF_DFI_RADAR == 1 )
5264          if ( PRESENT(dfi_stage) .and.  dfi_stage ==DFI_FWD            &
5265              .and. PRESENT(dfi_tten_rad) ) then
5266              WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
5267              CALL wrf_debug ( 50 , TRIM(wrf_err_message) )
5268          endif
5269      dfi_tten_max=-999
5270      old_max=-999
5271 #endif
5273 !  add microphysics theta diff to perturbation theta, set h_diabatic
5275      IF ( config_flags%no_mp_heating .eq. 0 ) THEN
5276        mptenmax = 0.
5277        mptenmin = 999.
5278      DO j = j_start, j_end
5279      DO k = k_start, k_end
5280      DO i = i_start, i_end
5281           mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
5282 #if ( WRF_DFI_RADAR == 1 )
5283        if(mpten.gt.mptenmax) then
5284           mptenmax=mpten
5285           imax=i
5286           jmax=j
5287        endif
5288        if(mpten.lt.mptenmin) then
5289           mptenmin=mpten
5290           imin=i
5291           jmin=j
5292        endif
5293           mpten=min(config_flags%mp_tend_lim*dt, mpten)
5294           mpten=max(-config_flags%mp_tend_lim*dt, mpten)
5296        if(k < k_end ) then
5297          if(dfi_tten_max < dfi_tten_rad(i,k,j) ) dfi_tten_max = dfi_tten_rad(i,k,j)
5298          if(old_max < (th_phy(i,k,j)-h_diabatic(i,k,j)) ) old_max=th_phy(i,k,j)-h_diabatic(i,k,j)
5299        endif
5301        if ( PRESENT(dfi_stage) .and.                                   &
5302           dfi_stage == DFI_FWD .and. PRESENT(dfi_tten_rad) .and.       &
5303           dfi_tten_rad(i,k,j) >= -0.1 .and. dfi_tten_rad(i,k,j) <= 0.1 &
5304           .and. k < k_end ) then
5305 ! add radar temp tendency
5306 ! there is radar coverage
5307            t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))*dt
5308        else
5309 ! no radar coverage
5310            t_new(i,k,j) = t_new(i,k,j) + mpten
5311        endif
5312 #else
5313          t_new(i,k,j) = t_new(i,k,j) + mpten
5314 #endif
5315          h_diabatic(i,k,j) =  mpten/dt
5316      ENDDO
5317      ENDDO
5318      ENDDO
5320      ELSE
5322      DO j = j_start, j_end
5323      DO k = k_start, k_end
5324      DO i = i_start, i_end
5325 !        t_new(i,k,j) = t_new(i,k,j)
5326          h_diabatic(i,k,j) = 0.
5327      ENDDO
5328      ENDDO
5329      ENDDO
5330      ENDIF
5332    END SUBROUTINE moist_physics_finish_em
5334 !----------------------------------------------------------------
5337    SUBROUTINE init_module_big_step
5338    END SUBROUTINE init_module_big_step
5340 SUBROUTINE set_tend ( field, field_adv_tend, msf,       &
5341                       ids, ide, jds, jde, kds, kde,     &
5342                       ims, ime, jms, jme, kms, kme,     &
5343                       its, ite, jts, jte, kts, kte       )
5345    IMPLICIT NONE
5347    ! Input data
5349    INTEGER ,  INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
5350                                ims, ime, jms, jme, kms, kme, &
5351                                its, ite, jts, jte, kts, kte
5353    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5355    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN)  :: field_adv_tend
5357    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN)  :: msf
5359    ! Local data
5361    INTEGER :: i, j, k, itf, jtf, ktf
5363 !<DESCRIPTION>
5365 !  set_tend copies the advective tendency array into the tendency array.
5367 !</DESCRIPTION>
5369       jtf = MIN(jte,jde-1)
5370       ktf = MIN(kte,kde-1)
5371       itf = MIN(ite,ide-1)
5372       DO j = jts, jtf
5373       DO k = kts, ktf
5374       DO i = its, itf
5375          field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5376       ENDDO
5377       ENDDO
5378       ENDDO
5380 END SUBROUTINE set_tend
5382 !------------------------------------------------------------------------------
5384     SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf,              &
5385                                  rw_tendf, t_tendf,               &
5386                                  u, v, w, t, t_init,              &
5387                                  mut, muu, muv, ph, phb,          &
5388                                  u_base, v_base, t_base, z_base,  &
5389                                  dampcoef, zdamp,                 &
5390                                  ids, ide, jds, jde, kds, kde,    &
5391                                  ims, ime, jms, jme, kms, kme,    &
5392                                  its, ite, jts, jte, kts, kte   )
5394 ! History:     Apr 2005  Modifications by George Bryan, NCAR:
5395 !                  - Generalized the code in a way that allows for
5396 !                    simulations with steep terrain.
5398 !              Jul 2004  Modifications by George Bryan, NCAR:
5399 !                  - Modified the code to use u_base, v_base, and t_base
5400 !                    arrays for the background state.  Removed the hard-wired
5401 !                    base-state values.
5402 !                  - Modified the code to use dampcoef, zdamp, and damp_opt,
5403 !                    i.e., the upper-level damper variables in namelist.input.
5404 !                    Removed the hard-wired variables in the older version.
5405 !                    This damper is used when damp_opt = 2.
5406 !                  - Modified the code to account for the movement of the
5407 !                    model surfaces with time.  The code now obtains a base-
5408 !                    state value by interpolation using the "_base" arrays.
5410 !              Nov 2003  Bug fix by Jason Knievel, NCAR
5412 !              Aug 2003  Meridional dimension, some comments, and
5413 !                        changes in layout of the code added by
5414 !                        Jason Knievel, NCAR
5416 !              Jul 2003  Original code by Bill Skamarock, NCAR
5418 ! Purpose:     This routine applies Rayleigh damping to a layer at top
5419 !              of the model domain.
5421 !-----------------------------------------------------------------------
5422 ! Begin declarations.
5424     IMPLICIT NONE
5426     INTEGER, INTENT( IN )  &
5427     :: ids, ide, jds, jde, kds, kde,  &
5428        ims, ime, jms, jme, kms, kme,  &
5429        its, ite, jts, jte, kts, kte
5431     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
5432     :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5434     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
5435     :: u, v, w, t, t_init, ph, phb
5437     REAL, DIMENSION( ims:ime, jms:jme ),  INTENT( IN )  &
5438     :: mut, muu, muv
5440     REAL, DIMENSION( kms:kme ) ,  INTENT(IN   )  &
5441     :: u_base, v_base, t_base, z_base
5443     REAL, INTENT(IN   )   &
5444     :: dampcoef, zdamp
5446 ! Local variables.
5448     INTEGER  &
5449     :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5451     REAL  &
5452     :: pii, dcoef, z, ztop
5454     REAL :: wkp1, wk, wkm1
5456     REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5458 ! End declarations.
5459 !-----------------------------------------------------------------------
5461     pii = 2.0 * asin(1.0)
5463     ktf = MIN( kte,   kde-1 )
5465 !-----------------------------------------------------------------------
5466 ! Adjust u to base state.
5468     DO j = jts, MIN( jte, jde-1 )
5469     DO i = its, MIN( ite, ide   )
5471       ! Get height at top of model
5472       ztop = 0.5*( phb(i  ,kde,j)+phb(i-1,kde,j)   &
5473                   +ph(i  ,kde,j)+ph(i-1,kde,j) )/g
5475       ! Find bottom of damping layer
5476       k1 = ktf
5477       z = ztop
5478       DO WHILE( z >= (ztop-zdamp) )
5479         z = 0.25*( phb(i  ,k1,j)+phb(i  ,k1+1,j)  &
5480                   +phb(i-1,k1,j)+phb(i-1,k1+1,j)  &
5481                   +ph(i  ,k1,j)+ph(i  ,k1+1,j)    &
5482                   +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5483         z00(k1) = z
5484         k1 = k1 - 1
5485       ENDDO
5486       k1 = k1 + 2
5488       ! Get reference state at model levels
5489       DO k = k1, ktf
5490         k2 = ktf
5491         DO WHILE( z_base(k2) .gt. z00(k) )
5492           k2 = k2 - 1
5493         ENDDO
5494         if(k2+1.gt.ktf)then
5495           u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) )   &
5496                               * (     z00(k) - z_base(k2)   )   &
5497                               / ( z_base(k2) - z_base(k2-1) )
5498         else
5499           u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) )   &
5500                               * (       z00(k) - z_base(k2) )   &
5501                               / ( z_base(k2+1) - z_base(k2) )
5502         endif
5503       ENDDO
5505       ! Apply the Rayleigh damper
5506       DO k = k1, ktf
5507         dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5508         dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5509         ru_tendf(i,k,j) = ru_tendf(i,k,j) -                    &
5510                           muu(i,j) * ( dcoef * dampcoef ) *    &
5511                           ( u(i,k,j) - u00(k) )
5512       END DO
5514     END DO
5515     END DO
5517 ! End adjustment of u.
5518 !-----------------------------------------------------------------------
5520 !-----------------------------------------------------------------------
5521 ! Adjust v to base state.
5523     DO j = jts, MIN( jte, jde   )
5524     DO i = its, MIN( ite, ide-1 )
5526       ! Get height at top of model
5527       ztop = 0.5*( phb(i,kde,j  )+phb(i,kde,j-1)   &
5528                   +ph(i,kde,j  )+ph(i,kde,j-1) )/g
5530       ! Find bottom of damping layer
5531       k1 = ktf
5532       z = ztop
5533       DO WHILE( z >= (ztop-zdamp) )
5534         z = 0.25*( phb(i,k1,j  )+phb(i,k1+1,j  )  &
5535                   +phb(i,k1,j-1)+phb(i,k1+1,j-1)  &
5536                   +ph(i,k1,j  )+ph(i,k1+1,j  )    &
5537                   +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5538         z00(k1) = z
5539         k1 = k1 - 1
5540       ENDDO
5541       k1 = k1 + 2
5543       ! Get reference state at model levels
5544       DO k = k1, ktf
5545         k2 = ktf
5546         DO WHILE( z_base(k2) .gt. z00(k) )
5547           k2 = k2 - 1
5548         ENDDO
5549         if(k2+1.gt.ktf)then
5550           v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) )   &
5551                               * (     z00(k) - z_base(k2)   )   &
5552                               / ( z_base(k2) - z_base(k2-1) )
5553         else
5554           v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) )   &
5555                               * (       z00(k) - z_base(k2) )   &
5556                               / ( z_base(k2+1) - z_base(k2) )
5557         endif
5558       ENDDO
5560       ! Apply the Rayleigh damper
5561       DO k = k1, ktf
5562         dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5563         dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5564         rv_tendf(i,k,j) = rv_tendf(i,k,j) -                    &
5565                           muv(i,j) * ( dcoef * dampcoef ) *    &
5566                           ( v(i,k,j) - v00(k) )
5567       END DO
5569     END DO
5570     END DO
5572 ! End adjustment of v.
5573 !-----------------------------------------------------------------------
5575 !-----------------------------------------------------------------------
5576 ! Adjust w to base state.
5578     DO j = jts, MIN( jte,   jde-1 )
5579     DO i = its, MIN( ite,   ide-1 )
5580       ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5581       DO k = kts, MIN( kte,   kde   )
5582         z = ( phb(i,k,j) + ph(i,k,j) ) / g
5583         IF ( z >= (ztop-zdamp) ) THEN
5584           dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5585           dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5586           rw_tendf(i,k,j) = rw_tendf(i,k,j) -  &
5587                             mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5588         END IF
5589       END DO
5590     END DO
5591     END DO
5593 ! End adjustment of w.
5594 !-----------------------------------------------------------------------
5596 !-----------------------------------------------------------------------
5597 ! Adjust potential temperature to base state.
5599     DO j = jts, MIN( jte,   jde-1 )
5600     DO i = its, MIN( ite,   ide-1 )
5602       ! Get height at top of model
5603       ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5605       ! Find bottom of damping layer
5606       k1 = ktf
5607       z = ztop
5608       DO WHILE( z >= (ztop-zdamp) )
5609         z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) +  &
5610                      ph(i,k1,j) +  ph(i,k1+1,j) ) / g
5611         z00(k1) = z
5612         k1 = k1 - 1
5613       ENDDO
5614       k1 = k1 + 2
5616       ! Get reference state at model levels
5617       DO k = k1, ktf
5618         k2 = ktf
5619         DO WHILE( z_base(k2) .gt. z00(k) )
5620           k2 = k2 - 1
5621         ENDDO
5622         if(k2+1.gt.ktf)then
5623           t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) )   &
5624                               * (     z00(k) - z_base(k2)   )   &
5625                               / ( z_base(k2) - z_base(k2-1) )
5626         else
5627           t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) )   &
5628                               * (       z00(k) - z_base(k2) )   &
5629                               / ( z_base(k2+1) - z_base(k2) )
5630         endif
5631       ENDDO
5633       ! Apply the Rayleigh damper
5634       DO k = k1, ktf
5635         dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5636         dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5637         t_tendf(i,k,j) = t_tendf(i,k,j) -                      &
5638                          mut(i,j) * ( dcoef * dampcoef )  *    &
5639                          ( t(i,k,j) - t00(k) )
5640       END DO
5642     END DO
5643     END DO
5645 ! End adjustment of potential temperature.
5646 !-----------------------------------------------------------------------
5648     END SUBROUTINE rk_rayleigh_damp
5650 !==============================================================================
5651 !==============================================================================
5652                                                                                 
5653       SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt,  &
5654                                         config_flags,                   &
5655                                         diff_6th_opt, diff_6th_factor,  &
5656                                         ids, ide, jds, jde, kds, kde,   &
5657                                         ims, ime, jms, jme, kms, kme,   &
5658                                         its, ite, jts, jte, kts, kte )
5659                                                                                 
5660 ! History:       14 Nov 2006   Name of variable changed by Jason Knievel
5661 !                07 Jun 2006   Revised and generalized by Jason Knievel  
5662 !                25 Apr 2005   Original code by Jason Knievel, NCAR
5663                                                                                 
5664 ! Purpose:       Apply 6th-order, monotonic (flux-limited), numerical
5665 !                diffusion to 3-d velocity and to scalars.
5666                                                                                 
5667 ! References:    Ming Xue (MWR Aug 2000)
5668 !                Durran ("Numerical Methods for Wave Equations..." 1999)
5669 !                George Bryan (personal communication)
5671 !------------------------------------------------------------------------------
5672 ! Begin: Declarations.
5674     IMPLICIT NONE
5676     INTEGER, INTENT(IN)  &
5677     :: ids, ide, jds, jde, kds, kde,   &
5678        ims, ime, jms, jme, kms, kme,   &
5679        its, ite, jts, jte, kts, kte
5681     TYPE(grid_config_rec_type), INTENT(IN)  &
5682     :: config_flags
5684     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT)  &
5685     :: tendency
5687     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN)  &
5688     :: field
5690     REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN)  &
5691     :: mu
5693     REAL, INTENT(IN)  &
5694     :: dt
5696     REAL, INTENT(IN)  &
5697     :: diff_6th_factor
5699     INTEGER, INTENT(IN)  &
5700     :: diff_6th_opt
5702     CHARACTER(LEN=1) , INTENT(IN)  &
5703     :: name
5705     INTEGER  &
5706     :: i, j, k,         &
5707        i_start, i_end,  &
5708        j_start, j_end,  &
5709        k_start, k_end,  &
5710        ktf
5712     REAL  &
5713     :: dflux_x_p0, dflux_y_p0,  &
5714        dflux_x_p1, dflux_y_p1,  &
5715        tendency_x, tendency_y,  &
5716        mu_avg_p0, mu_avg_p1,    &
5717        diff_6th_coef
5719     LOGICAL  &
5720     :: specified
5722 ! End: Declarations.
5723 !------------------------------------------------------------------------------
5725 !------------------------------------------------------------------------------
5726 ! Begin: Translate the diffusion factor into a diffusion coefficient.  See
5727 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5728 ! fourth) and for diffusion in two dimensions (not one).  For reference, a
5729 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5730 ! although application of the flux limiter reduces somewhat the effects of
5731 ! diffusion for a given coefficient.
5733     diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )  
5735 ! End: Translate diffusion factor.
5736 !------------------------------------------------------------------------------
5738 !------------------------------------------------------------------------------
5739 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
5740 ! The halo regions are already filled with values by the time this subroutine
5741 ! is called, which allows the stencil to extend beyond the domains' edges.
5743     ktf = MIN( kte, kde-1 )
5745     IF ( name .EQ. 'u' ) THEN
5747       i_start = its
5748       i_end   = ite
5749       j_start = jts
5750       j_end   = MIN(jde-1,jte)
5751       k_start = kts
5752       k_end   = ktf
5754     ELSE IF ( name .EQ. 'v' ) THEN
5756       i_start = its
5757       i_end   = MIN(ide-1,ite)
5758       j_start = jts
5759       j_end   = jte
5760       k_start = kts
5761       k_end   = ktf
5763     ELSE IF ( name .EQ. 'w' ) THEN
5765       i_start = its
5766       i_end   = MIN(ide-1,ite)
5767       j_start = jts
5768       j_end   = MIN(jde-1,jte)
5769       k_start = kts+1
5770       k_end   = ktf
5772     ELSE
5774       i_start = its
5775       i_end   = MIN(ide-1,ite)
5776       j_start = jts
5777       j_end   = MIN(jde-1,jte)
5778       k_start = kts
5779       k_end   = ktf
5781     ENDIF
5783 ! End: Assignment of limits of spatial loops.
5784 !------------------------------------------------------------------------------
5786 !------------------------------------------------------------------------------
5787 ! Begin: Loop across spatial dimensions.
5789     DO j = j_start, j_end
5790     DO k = k_start, k_end
5791     DO i = i_start, i_end
5793 !------------------------------------------------------------------------------
5794 ! Begin: Diffusion in x (i index).
5796 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
5798       dflux_x_p0 = (  10.0 * ( field(i,  k,j) - field(i-1,k,j) )    &
5799                      - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) )    &
5800                      +       ( field(i+2,k,j) - field(i-3,k,j) ) )
5802       dflux_x_p1 = (  10.0 * ( field(i+1,k,j) - field(i  ,k,j) )    &
5803                      - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) )    &
5804                      +       ( field(i+3,k,j) - field(i-2,k,j) ) )
5806 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5807 ! (variation on Xue's eq. 10).
5809       IF ( diff_6th_opt .EQ. 2 ) THEN
5811         IF ( dflux_x_p0 * ( field(i  ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
5812           dflux_x_p0 = 0.0
5813         END IF
5815         IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i  ,k,j) ) .LE. 0.0 ) THEN
5816           dflux_x_p1 = 0.0
5817         END IF
5819       END IF
5821 ! Apply 6th-order diffusion in x direction.
5823       IF      ( name .EQ. 'u' ) THEN
5824         mu_avg_p0 = mu(i-1,j)
5825         mu_avg_p1 = mu(i  ,j)
5826       ELSE IF ( name .EQ. 'v' ) THEN
5827         mu_avg_p0 = 0.25 * (       &
5828                     mu(i-1,j-1) +  &
5829                     mu(i  ,j-1) +  &
5830                     mu(i-1,j  ) +  &
5831                     mu(i  ,j  ) )
5832         mu_avg_p1 = 0.25 * (       &
5833                     mu(i  ,j-1) +  &
5834                     mu(i+1,j-1) +  &
5835                     mu(i  ,j  ) +  &
5836                     mu(i+1,j  ) )
5837       ELSE
5838         mu_avg_p0 = 0.5 * (        &
5839                     mu(i-1,j) +    &
5840                     mu(i  ,j) )
5841         mu_avg_p1 = 0.5 * (        &
5842                     mu(i  ,j) +    &
5843                     mu(i+1,j) )
5844       END IF
5846       tendency_x = diff_6th_coef *  &
5847                  ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
5849 ! End: Diffusion in x.
5850 !------------------------------------------------------------------------------
5852 !------------------------------------------------------------------------------
5853 ! Begin: Diffusion in y (j index).
5855 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
5857       dflux_y_p0 = (  10.0 * ( field(i,k,j  ) - field(i,k,j-1) )    &
5858                      - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) )    &
5859                      +       ( field(i,k,j+2) - field(i,k,j-3) ) )
5861       dflux_y_p1 = (  10.0 * ( field(i,k,j+1) - field(i,k,j  ) )    &
5862                      - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) )    &
5863                      +       ( field(i,k,j+3) - field(i,k,j-2) ) )
5865 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5866 ! (variation on Xue's eq. 10).
5868       IF ( diff_6th_opt .EQ. 2 ) THEN
5870         IF ( dflux_y_p0 * ( field(i,k,j  )-field(i,k,j-1) ) .LE. 0.0 ) THEN
5871           dflux_y_p0 = 0.0
5872         END IF
5874         IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j  ) ) .LE. 0.0 ) THEN
5875           dflux_y_p1 = 0.0
5876         END IF
5878       END IF
5880 ! Apply 6th-order diffusion in y direction.
5882       IF      ( name .EQ. 'u' ) THEN
5883         mu_avg_p0 = 0.25 * (       &
5884                     mu(i-1,j-1) +  &
5885                     mu(i  ,j-1) +  &
5886                     mu(i-1,j  ) +  &
5887                     mu(i  ,j  ) )
5888         mu_avg_p1 = 0.25 * (       &
5889                     mu(i-1,j  ) +  &
5890                     mu(i  ,j  ) +  &
5891                     mu(i-1,j+1) +  &
5892                     mu(i  ,j+1) )
5893       ELSE IF ( name .EQ. 'v' ) THEN
5894         mu_avg_p0 = mu(i,j-1)
5895         mu_avg_p1 = mu(i,j  )
5896       ELSE
5897         mu_avg_p0 = 0.5 * (      &
5898                     mu(i,j-1) +  &
5899                     mu(i,j  ) )
5900         mu_avg_p1 = 0.5 * (      &
5901                     mu(i,j  ) +  &
5902                     mu(i,j+1) )
5903       END IF
5905       tendency_y = diff_6th_coef *  &
5906                  ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
5908 ! End: Diffusion in y.
5909 !------------------------------------------------------------------------------
5911 !------------------------------------------------------------------------------
5912 ! Begin: Combine diffusion in x and y.
5913      
5914       tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
5916 ! End: Combine diffusion in x and y.
5917 !------------------------------------------------------------------------------
5919     ENDDO
5920     ENDDO
5921     ENDDO
5923 ! End: Loop across spatial dimensions.
5924 !------------------------------------------------------------------------------
5926     END SUBROUTINE sixth_order_diffusion
5928 !==============================================================================
5929 !==============================================================================
5931 END MODULE module_big_step_utilities_em