module_big_step_utilities_em.F

References to this file elsewhere.
1 !WRF:MODEL_LAYER:DYNAMICS
2 !
3 
4 #if (RWORDSIZE == 4)
5 #   define VPOWX vspowx
6 #   define VPOW  vspow
7 #else
8 #   define VPOWX vpowx
9 #   define VPOW  vpow
10 #endif
11 
12 
13 MODULE module_big_step_utilities_em
14 
15    USE module_domain
16    USE module_model_constants
17    USE module_state_description
18    USE module_configure
19    USE module_wrf_error
20 
21 CONTAINS
22 
23 !-------------------------------------------------------------------------------
24 
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 )
30 
31    IMPLICIT NONE
32    
33    ! Input data
34 
35    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
36 
37    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
38                                        ims, ime, jms, jme, kms, kme, &
39                                        its, ite, jts, jte, kts, kte 
40 
41    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
42    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub
43 
44    !  local stuff
45 
46    INTEGER :: i, j, itf, jtf, im, jm
47 
48 !<DESCRIPTION>
49 !
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)
53 !
54 !</DESCRIPTION>
55 
56 
57       itf=ite
58       jtf=MIN(jte,jde-1)
59 
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
117 
118       itf=MIN(ite,ide-1)
119       jtf=jte
120 
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
178 
179 END SUBROUTINE calc_mu_uv
180 
181 !-------------------------------------------------------------------------------
182 
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 )
188 
189    IMPLICIT NONE
190    
191    ! Input data
192 
193    TYPE(grid_config_rec_type   ) ,   INTENT(IN   ) :: config_flags
194 
195    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
196                                        ims, ime, jms, jme, kms, kme, &
197                                        its, ite, jts, jte, kts, kte 
198 
199    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
200    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
201 
202    !  local stuff
203 
204    INTEGER :: i, j, itf, jtf, im, jm
205 
206 !<DESCRIPTION>
207 !
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)
211 !
212 !</DESCRIPTION>
213    
214       itf=ite
215       jtf=MIN(jte,jde-1)
216 
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
266 
267       itf=MIN(ite,ide-1)
268       jtf=jte
269 
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
319 
320 END SUBROUTINE calc_mu_uv_1
321 
322 !-------------------------------------------------------------------------------
323 
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)
327 
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  )
334 
335    IMPLICIT NONE
336 
337    ! Input data
338 
339    INTEGER ,             INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
340                                           ims, ime, jms, jme, kms, kme, &
341                                           its, ite, jts, jte, kts, kte
342 
343    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: ru, rv, rw
344 
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>
355 !
356 ! couple_momentum couples the velocities to the full column mass and
357 ! the map factors.
358 !
359 !</DESCRIPTION>
360 
361    ktf=MIN(kte,kde-1)
362    
363       itf=ite
364       jtf=MIN(jte,jde-1)
365 
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
373 
374       itf=MIN(ite,ide-1)
375       jtf=jte
376 
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
385 
386       itf=MIN(ite,ide-1)
387       jtf=MIN(jte,jde-1)
388 
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
396 
397 END SUBROUTINE couple_momentum
398 
399 !-------------------------------------------------------------------
400 
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 )
405 
406    IMPLICIT NONE
407    
408    ! Input data
409 
410    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
411                                        ims, ime, jms, jme, kms, kme, &
412                                        its, ite, jts, jte, kts, kte 
413 
414    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(  OUT) :: muu, muv
415    REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu, mub
416 
417    !  local stuff
418 
419    INTEGER :: i, j, itf, jtf
420 
421 !<DESCRIPTION>
422 !
423 ! calc_mu_staggered calculates the full dry air mass at the staggered
424 ! velocity points (u,v).
425 !
426 !</DESCRIPTION>
427    
428       itf=ite
429       jtf=MIN(jte,jde-1)
430 
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
472 
473       itf=MIN(ite,ide-1)
474       jtf=jte
475 
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
517 
518 END SUBROUTINE calc_mu_staggered
519 
520 !-------------------------------------------------------------------------------
521 
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 )
527 
528    IMPLICIT NONE
529 
530    ! Input data
531 
532    INTEGER ,             INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
533                                           ims, ime, jms, jme, kms, kme, &
534                                           its, ite, jts, jte, kts, kte
535 
536    CHARACTER(LEN=1) ,     INTENT(IN   ) :: name
537 
538    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(  OUT) :: rfield
539 
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
548 
549 !<DESCRIPTION>
550 !
551 ! subroutine couple couples the input variable with the dry-air 
552 ! column mass (mu).  
553 !
554 !</DESCRIPTION>
555 
556    
557    ktf=MIN(kte,kde-1)
558    
559    IF (name .EQ. 'u')THEN
560 
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 )
565 
566       itf=ite
567       jtf=MIN(jte,jde-1)
568 
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
576 
577    ELSE IF (name .EQ. 'v')THEN
578 
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 )
583 
584       itf=ite
585       itf=MIN(ite,ide-1)
586       jtf=jte
587 
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
595 
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
606 
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
617 
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
630 
631 END SUBROUTINE couple
632 
633 
634 !-------------------------------------------------------------------------------
635 
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    )
643 
644    IMPLICIT NONE
645 
646    ! Input data
647 
648 
649    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
650                                  ims, ime, jms, jme, kms, kme, &
651                                  its, ite, jts, jte, kts, kte
652 
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
670 
671 !<DESCRIPTION>
672 !
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.
677 !
678 !</DESCRIPTION>
679 
680 !<DESCRIPTION>
681 !
682 !  calc_ww_cp calculates omega using the velocities (u,v) and the 
683 !  column mass mu.
684 !
685 !</DESCRIPTION>
686 
687     jtf=MIN(jte,jde-1)
688     ktf=MIN(kte,kde-1)  
689     itf=MIN(ite,ide-1)
690 
691 !  mu coupled with the appropriate map factor
692 
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
699 
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
707 
708       DO j=jts,jtf
709 
710         DO i=its,ite
711           dmdt(i) = 0.
712           ww(i,1,j) = 0.
713           ww(i,kte,j) = 0.
714         ENDDO
715 
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)
721 !
722 !       Using nu instead of z the last term becomes:
723 !                                -partial d/dnu(mu (dnu/dt)/my)
724 !
725 !       Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
726 !       and bottom, the last term becomes = 0
727 !
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
731 !
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
735 !
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
741 
742         DO k=kts,ktf
743         DO i=its,itf
744 
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))   )
747 
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))   )
750 
751           dmdt(i) = dmdt(i) + divv(i,k)
752 
753 
754         ENDDO
755         ENDDO
756 
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]
764 
765         DO k=2,ktf
766         DO i=its,itf
767 
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)) )
772 
773            ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1)
774 
775         ENDDO
776         ENDDO
777      ENDDO
778 
779 
780 END SUBROUTINE calc_ww_cp
781 
782 
783 !-------------------------------------------------------------------------------
784  
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  )
789 
790    IMPLICIT NONE
791    
792    ! Input data
793 
794    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
795                                        ims, ime, jms, jme, kms, kme, &
796                                        its, ite, jts, jte, kts, kte 
797 
798    INTEGER ,          INTENT(IN   ) :: n_moist
799    
800 
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
804 
805    ! Local stuff
806 
807    REAL :: qtot
808    
809    INTEGER :: i, j, k, itf, jtf, ktf, ispe
810 
811 !<DESCRIPTION>
812 !
813 !  calc_cq calculates moist coefficients for the momentum equations.
814 !
815 !</DESCRIPTION>
816 
817       itf=ite
818       jtf=MIN(jte,jde-1)
819       ktf=MIN(kte,kde-1)
820 
821       IF(  n_moist >= PARAM_FIRST_SCALAR ) THEN
822 
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
838 
839         itf=MIN(ite,ide-1)
840         jtf=jte
841 
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
857 
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
875 
876       ELSE
877 
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
885 
886         itf=MIN(ite,ide-1)
887         jtf=jte
888 
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
896 
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
906 
907       END IF
908 
909 END SUBROUTINE calc_cq
910 
911 !----------------------------------------------------------------------
912 
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  )
917 
918    IMPLICIT NONE
919    
920    ! Input data
921 
922    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
923                                        ims, ime, jms, jme, kms, kme, &
924                                        its, ite, jts, jte, kts, kte 
925 
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
928 
929    ! Local stuff
930 
931    INTEGER :: i, j, k, itf, jtf, ktf
932 
933 !<DESCRIPTION>
934 !
935 ! calc_alt computes the full inverse density
936 !
937 !</DESCRIPTION>
938 
939       itf=MIN(ite,ide-1)
940       jtf=MIN(jte,jde-1)
941       ktf=MIN(kte,kde-1)
942 
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
950 
951 
952 END SUBROUTINE calc_alt
953 
954 !----------------------------------------------------------------------
955 
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  )
963 
964   IMPLICIT NONE
965    
966    ! Input data
967 
968   LOGICAL ,          INTENT(IN   ) :: non_hydrostatic
969 
970   INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
971                                       ims, ime, jms, jme, kms, kme, &
972                                       its, ite, jts, jte, kts, kte 
973 
974   INTEGER ,          INTENT(IN   ) :: n_moist
975 
976   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN   ) :: alb,  &
977                                                                    pb,   &
978                                                                    t
979 
980   REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN   ) :: moist
981 
982   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(  OUT) :: al, p
983 
984   REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph
985 
986   REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN   ) :: mu, muts
987 
988   REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: znu, dnw, rdnw, rdn
989 
990   REAL,   INTENT(IN   ) :: t0, p0
991 
992   ! Local stuff
993 
994   INTEGER :: i, j, k, itf, jtf, ktf, ispe
995   REAL    :: qvf, qtot, qf1, qf2
996   REAL, DIMENSION( its:ite) :: temp,cpovcv_v
997 
998 
999 !<DESCRIPTION>
1000 !
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.
1004 !
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.
1009 !
1010 !</DESCRIPTION>
1011 
1012   itf=MIN(ite,ide-1)
1013   jtf=MIN(jte,jde-1)
1014   ktf=MIN(kte,kde-1)
1015 
1016 #ifndef INTELMKL
1017   cpovcv_v = cpovcv
1018 #endif
1019 
1020   IF (non_hydrostatic) THEN
1021 
1022       IF (n_moist >= PARAM_FIRST_SCALAR ) THEN  
1023 
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
1044 
1045       ELSE
1046 
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
1058 
1059       END IF
1060 
1061    ELSE
1062 
1063 !  hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1064 
1065 
1066       IF (n_moist >= PARAM_FIRST_SCALAR ) THEN  
1067 
1068         DO j=jts,jtf
1069 
1070           k=ktf          ! top layer
1071           DO i=its,itf
1072 
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
1079 
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)
1084 
1085           ENDDO
1086 
1087           DO k=ktf-1,kts,-1  ! remaining layers, integrate down
1088             DO i=its,itf
1089 
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
1096 
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
1103 
1104           DO k=2,ktf+1  ! integrate hydrostatic equation for geopotential
1105             DO i=its,itf
1106 
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                                                    
1114 
1115             ENDDO
1116           ENDDO
1117 
1118         ENDDO
1119 
1120       ELSE
1121 
1122         DO j=jts,jtf
1123 
1124           k=ktf          ! top layer
1125           DO i=its,itf
1126 
1127             qtot = 0.
1128             qf2 = 1./(1.+qtot)
1129             qf1 = qtot*qf2
1130 
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)
1135 
1136           ENDDO
1137 
1138           DO k=ktf-1,kts,-1  ! remaining layers, integrate down
1139             DO i=its,itf
1140 
1141             qtot = 0.
1142             qf2 = 1./(1.+qtot)
1143             qf1 = qtot*qf2
1144 
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
1151 
1152           DO k=2,ktf+1  ! integrate hydrostatic equation for geopotential
1153             DO i=its,itf
1154 
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                                                    
1162 
1163             ENDDO
1164           ENDDO
1165 
1166         ENDDO
1167 
1168      END IF
1169 
1170    END IF
1171 
1172 END SUBROUTINE calc_p_rho_phi
1173 
1174 !----------------------------------------------------------------------
1175 
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  )
1180 
1181    IMPLICIT NONE
1182    
1183    ! Input data
1184 
1185    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1186                                        ims, ime, jms, jme, kms, kme, &
1187                                        its, ite, jts, jte, kts, kte 
1188 
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
1191 
1192    ! Local stuff
1193 
1194    INTEGER :: i, j, k, itf, jtf, ktf
1195 
1196 !<DESCRIPTION>
1197 !
1198 !  calc_php calculates the full geopotential from the reference state
1199 !  geopotential and the perturbation geopotential (phb_ph).
1200 !
1201 !</DESCRIPTION>
1202 
1203       itf=MIN(ite,ide-1)
1204       jtf=MIN(jte,jde-1)
1205       ktf=MIN(kte,kde-1)
1206 
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
1214 
1215 END SUBROUTINE calc_php
1216 
1217 !-------------------------------------------------------------------------------
1218 
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        )
1226 
1227    IMPLICIT NONE
1228 
1229    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1230                                        ims, ime, jms, jme, kms, kme, &
1231                                        its, ite, jts, jte, kts, kte 
1232 
1233    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   ph_tend, &
1234                                                                      ph_new,  &
1235                                                                      ph_old,  &
1236                                                                      u,       &
1237                                                                      v
1238 
1239 
1240    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(  OUT) :: w
1241 
1242    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mu, ht, msftx, msfty
1243 
1244    REAL, INTENT(IN   ) :: dt, cf1, cf2, cf3, rdx, rdy
1245 
1246    INTEGER :: i, j, k, itf, jtf
1247 
1248    itf=MIN(ite,ide-1)
1249    jtf=MIN(jte,jde-1)
1250 
1251 !<DESCRIPTION>
1252 !
1253 ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1254 ! Used with the hydrostatic option.
1255 !
1256 !</DESCRIPTION>
1257 
1258    DO j = jts, jtf
1259 
1260 !  lower b.c. on w
1261 
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]
1272 
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
1287 
1288 !  use geopotential equation to diagnose w
1289 
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]
1296 !  THE FOLLOWING SEEMS TO IMPLY THAT ph_new, old ARE ALREADY divided by mu ???
1297 
1298      DO k = 2, kte
1299      DO i = its, itf
1300        w(i,k,j) =  msfty(i,j)*(  (ph_new(i,k,j)-ph_old(i,k,j))/dt       &
1301                                - ph_tend(i,k,j)/mu(i,j)        )/g 
1302 
1303      ENDDO
1304      ENDDO
1305 
1306    ENDDO
1307 
1308 END SUBROUTINE diagnose_w
1309 
1310 !-------------------------------------------------------------------------------
1311 
1312 SUBROUTINE rhs_ph( ph_tend, u, v, ww,               &
1313                    ph, ph_old, phb, w,              &
1314                    mut, muu, muv,                   &
1315                    fnm, fnp,                        &
1316                    rdnw, cfn, cfn1, rdx, rdy,       &
1317                    msfux, msfuy, msfvx,             &
1318                    msfvx_inv, msfvy,                &
1319                    msftx, msfty,                    &
1320                    non_hydrostatic,                 &
1321                    config_flags,                    &
1322                    ids, ide, jds, jde, kds, kde,    &
1323                    ims, ime, jms, jme, kms, kme,    &
1324                    its, ite, jts, jte, kts, kte    )
1325    IMPLICIT NONE
1326 
1327    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
1328 
1329    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1330                                        ims, ime, jms, jme, kms, kme, &
1331                                        its, ite, jts, jte, kts, kte 
1332 
1333    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::        &
1334                                                                      u,   &
1335                                                                      v,   &
1336                                                                      ww,  &
1337                                                                      ph,  &
1338                                                                      ph_old, &
1339                                                                      phb, & 
1340                                                                     w
1341 
1342 ! pjj/cray
1343 !  REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(  OUT) :: ph_tend
1344    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1345 
1346    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: muu, muv, mut,   &
1347                                                             msfux, msfuy, &
1348                                                             msfvx, msfvy, &
1349                                                             msftx, msfty, &
1350                                                             msfvx_inv
1351 
1352    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, fnm, fnp
1353 
1354    REAL,  INTENT(IN   ) :: cfn, cfn1, rdx, rdy
1355 
1356    LOGICAL,  INTENT(IN   )  ::  non_hydrostatic
1357 
1358    ! Local stuff
1359 
1360    INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1361    REAL    :: ur, ul, ub, vr, vl, vb
1362    REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1363 
1364    INTEGER :: advective_order
1365 
1366    LOGICAL :: specified
1367 
1368 !<DESCRIPTION>
1369 !
1370 ! rhs_ph calculates the large-timestep tendency terms for the geopotential
1371 ! equation.  These terms include the advection and "gw".  The geopotential
1372 ! equation is cast in advective form, so we don't use the flux form advection
1373 ! algorithms here.
1374 !
1375 !</DESCRIPTION>
1376 
1377    specified = .false.
1378    if(config_flags%specified .or. config_flags%nested) specified = .true.
1379 
1380    advective_order = config_flags%h_sca_adv_order 
1381 !   advective_order = 2  !  original configuration (pre Oct 2001)
1382 
1383    itf=MIN(ite,ide-1)
1384    jtf=MIN(jte,jde-1)
1385    ktf=MIN(kte,kde-1)
1386 
1387 !  Notes on map scale factors:
1388 !  phi equation is:   partial d/dt(mu phi/my) = -mx partial d/dx(mu rho u/my)
1389 !                                               -mx partial d/dy(phi mu v/mx)
1390 !                                               -partial d/dz(phi mu w/my)
1391 !                                               +mu g w/my
1392 
1393 ! advective form for the geopotential equation
1394 
1395    DO j = jts, jtf
1396 
1397      DO k = 2, kte
1398      DO i = its, itf
1399           wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1)               &
1400                         *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1401      ENDDO
1402      ENDDO
1403 
1404 !  Further notes on map scale factors
1405 !  Using nu not z, RHS term 3 is: - partial d/dnu(phi mu dnu/dt/my)
1406 !                               = - partial d/dnu(phi ww)
1407 
1408      DO k = 2, kte-1
1409      DO i = its, itf
1410            ph_tend(i,k,j) = ph_tend(i,k,j)                           &
1411                              - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1412      ENDDO
1413      ENDDO
1414 
1415    ENDDO
1416 
1417    IF (non_hydrostatic) THEN  ! add in "gw" term.
1418    DO j = jts, jtf            ! in hydrostatic mode, "gw" will be diagnosed
1419                               ! after the timestep to give us "w"
1420      DO i = its, itf
1421         ph_tend(i,kde,j) = 0.
1422      ENDDO
1423 
1424      DO k = 2, kte
1425      DO i = its, itf
1426         ! phi equation RHS term 4
1427         ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msfty(i,j)
1428      ENDDO
1429      ENDDO
1430 
1431    ENDDO
1432 
1433    END IF
1434 
1435 !  Notes on map scale factors:
1436 !  RHS terms 1 and 2 are: -mx partial d/dx(mu rho u/my)
1437 !                         -mx partial d/dy(phi mu v/mx)
1438 !  Here v, u, phb, ph are not coupled with mu or map scale factors
1439 !  ?? shouldn't there be 1/m inside by each v, u and then mx outside ??
1440 !  ?? maybe it's equivalent to using 'avgeraged' m outside (as given by msft)??
1441 
1442    IF (advective_order <= 2) THEN
1443 
1444 !  y (v) advection
1445 
1446    i_start = its
1447    j_start = jts
1448    itf=MIN(ite,ide-1)
1449    jtf=MIN(jte,jde-1)
1450 
1451    IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1452    IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1453 
1454    DO j = j_start, jtf
1455 
1456      DO k = 2, kte-1
1457      DO i = i_start, itf
1458         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdy*                       &
1459                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*               &
1460                   (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))*msfvx_inv(i,j+1)  &
1461                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*               &
1462                   (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))*msfvx_inv(i,j) )
1463      ENDDO
1464      ENDDO
1465 
1466      k = kte
1467      DO i = i_start, itf
1468         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdy*                         &
1469                   ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*    &
1470                    (phb(i,k,j+1)-phb(i,k,j  )+ph(i,k,j+1)-ph(i,k,j  ))*msfvx_inv(i,j+1)  &
1471                    +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*    &
1472                    (phb(i,k,j  )-phb(i,k,j-1)+ph(i,k,j  )-ph(i,k,j-1))*msfvx_inv(i,j) )
1473      ENDDO
1474 
1475    ENDDO
1476 
1477 !  x (u) advection
1478 
1479    i_start = its
1480    j_start = jts
1481    itf=MIN(ite,ide-1)
1482    jtf=MIN(jte,jde-1)
1483 
1484    IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1485    IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1486 
1487    DO j = j_start, jtf
1488 
1489      DO k = 2, kte-1
1490      DO i = i_start, itf
1491         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdx*                        &
1492                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))*                &
1493                   (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))/msfuy(i+1,j)   &
1494                   +muu(i  ,j)*(u(i  ,k,j)+u(i  ,k-1,j))*                &
1495                   (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))/msfuy(i,j) )
1496      ENDDO
1497      ENDDO
1498  
1499      k = kte
1500      DO i = i_start, itf
1501         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdx*                         &
1502                   ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*    &
1503                    (phb(i+1,k,j)-phb(i  ,k,j)+ph(i+1,k,j)-ph(i  ,k,j))/msfuy(i+1,j)  &
1504                    +muu(i  ,j)*(cfn*u(i  ,k-1,j)+cfn1*u(i  ,k-2,j))*    &
1505                    (phb(i  ,k,j)-phb(i-1,k,j)+ph(i  ,k,j)-ph(i-1,k,j))/msfuy(i,j) )
1506      ENDDO
1507 
1508    ENDDO
1509 
1510    ELSE IF (advective_order <= 4) THEN
1511 
1512 !  y (v) advection
1513 
1514    i_start = its
1515    j_start = jts
1516    itf=MIN(ite,ide-1)
1517    jtf=MIN(jte,jde-1)
1518 
1519    IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1520    IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1521 
1522    DO j = j_start, jtf
1523 
1524      DO k = 2, kte-1
1525      DO i = i_start, itf
1526         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdy*           (          &
1527                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvx_inv(i,j+1)               &
1528                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvx_inv(i,j) )* (1./12.)*( &
1529                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                      &
1530                       -(ph(i,k,j+2)-ph(i,k,j-2))                      &
1531                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                    &
1532                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1533 
1534 
1535      ENDDO
1536      ENDDO
1537 
1538      k = kte
1539      DO i = i_start, itf
1540         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdy*           (                      &
1541                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvx_inv(i,j+1)               &
1542                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvx_inv(i,j) )* (1./12.)*( &
1543                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                 &
1544                       -(ph(i,k,j+2)-ph(i,k,j-2))                                 &
1545                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                               &
1546                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1547 
1548      ENDDO
1549 
1550    ENDDO
1551 
1552 
1553 !  x (u) advection
1554 
1555    i_start = its
1556    j_start = jts
1557    itf=MIN(ite,ide-1)
1558    jtf=MIN(jte,jde-1)
1559 
1560    IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1
1561    IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1
1562 
1563    DO j = j_start, jtf
1564 
1565      DO k = 2, kte-1
1566      DO i = i_start, itf
1567         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdx*(                     &
1568                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))/msfuy(i+1,j)               &
1569                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))/msfuy(i,j) )* (1./12.)*( &
1570                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                      &
1571                       -(ph(i+2,k,j)-ph(i-2,k,j))                      &
1572                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                    &
1573                       -(phb(i+2,k,j)-phb(i-2,k,j))  )   )                
1574      ENDDO
1575      ENDDO
1576  
1577      k = kte
1578      DO i = i_start, itf
1579         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdx*(                                 &
1580                  ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))/msfuy(i+1,j)               &
1581                   +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))/msfuy(i,j) )* (1./12.)*(   &
1582                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                                 &
1583                       -(ph(i+2,k,j)-ph(i-2,k,j))                                 &
1584                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                               &
1585                       -(phb(i+2,k,j)-phb(i-2,k,j))  )     )
1586      ENDDO
1587 
1588    ENDDO
1589 
1590    ELSE IF (advective_order <= 6) THEN
1591 
1592 !  y (v) advection
1593 
1594    i_start = its
1595    j_start = jts
1596    itf=MIN(ite,ide-1)
1597    jtf=MIN(jte,jde-1)
1598 
1599 !   IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1
1600 !   IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1
1601 
1602    IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+2)
1603    IF (config_flags%open_ye .or. specified ) jtf     = min(jtf,jde-3)
1604 
1605    DO j = j_start, jtf
1606 
1607      DO k = 2, kte-1
1608      DO i = i_start, itf
1609         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdy* (                    &
1610                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvx_inv(i,j+1)               &
1611                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvx_inv(i,j) )* (1./60.)*( &
1612                    45.*(ph(i,k,j+1)-ph(i,k,j-1))                      &
1613                    -9.*(ph(i,k,j+2)-ph(i,k,j-2))                      &
1614                       +(ph(i,k,j+3)-ph(i,k,j-3))                      &
1615                   +45.*(phb(i,k,j+1)-phb(i,k,j-1))                    &
1616                    -9.*(phb(i,k,j+2)-phb(i,k,j-2))                    &
1617                       +(phb(i,k,j+3)-phb(i,k,j-3))  )   )                
1618 
1619 
1620      ENDDO
1621      ENDDO
1622 
1623      k = kte
1624      DO i = i_start, itf
1625         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdy* (                                &
1626                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvx_inv(i,j+1)               &
1627                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvx_inv(i,j) )* (1./60.)*( &
1628                    45.*(ph(i,k,j+1)-ph(i,k,j-1))                                 &
1629                    -9.*(ph(i,k,j+2)-ph(i,k,j-2))                                 &
1630                       +(ph(i,k,j+3)-ph(i,k,j-3))                                 &
1631                   +45.*(phb(i,k,j+1)-phb(i,k,j-1))                               &
1632                    -9.*(phb(i,k,j+2)-phb(i,k,j-2))                               &
1633                       +(phb(i,k,j+3)-phb(i,k,j-3))  )   )                
1634 
1635      ENDDO
1636 
1637    ENDDO
1638 
1639 
1640 !  pick up near boundary rows using 4th order stencil 
1641 !  (open bc copy only goes out to jds-1 and jde, hence 4rth is ok but 6th is too big)
1642 
1643    IF ( (config_flags%open_ys) .and. jts <= jds+1 )  THEN
1644 
1645      j = jds+1
1646      DO k = 2, kte-1
1647      DO i = i_start, itf
1648         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdy* (                    &
1649                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvx_inv(i,j+1)               &
1650                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvx_inv(i,j) )* (1./12.)*( &
1651                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                      &
1652                       -(ph(i,k,j+2)-ph(i,k,j-2))                      &
1653                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                    &
1654                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1655 
1656 
1657      ENDDO
1658      ENDDO
1659 
1660      k = kte
1661      DO i = i_start, itf
1662         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdy* (                                &
1663                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvx_inv(i,j+1)               &
1664                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvx_inv(i,j) )* (1./12.)*( &
1665                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                 &
1666                       -(ph(i,k,j+2)-ph(i,k,j-2))                                 &
1667                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                               &
1668                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1669 
1670      ENDDO
1671 
1672    END IF
1673 
1674    IF ( (config_flags%open_ye) .and. jte >= jde-2 )  THEN
1675 
1676      j = jde-2
1677      DO k = 2, kte-1
1678      DO i = i_start, itf
1679         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdy* (                    &
1680                  ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))*msfvx_inv(i,j+1)               &
1681                   +muv(i,j  )*(v(i,k,j  )+v(i,k-1,j  ))*msfvx_inv(i,j) )* (1./12.)*( &
1682                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                      &
1683                       -(ph(i,k,j+2)-ph(i,k,j-2))                      &
1684                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                    &
1685                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1686 
1687 
1688      ENDDO
1689      ENDDO
1690 
1691      k = kte
1692      DO i = i_start, itf
1693         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdy* (                                &
1694                  ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvx_inv(i,j+1)               &
1695                   +muv(i,j  )*(cfn*v(i,k-1,j  )+cfn1*v(i,k-2,j  ))*msfvx_inv(i,j) )* (1./12.)*( &
1696                     8.*(ph(i,k,j+1)-ph(i,k,j-1))                                 &
1697                       -(ph(i,k,j+2)-ph(i,k,j-2))                                 &
1698                    +8.*(phb(i,k,j+1)-phb(i,k,j-1))                               &
1699                       -(phb(i,k,j+2)-phb(i,k,j-2))  )   )                
1700 
1701      ENDDO
1702 
1703    END IF
1704 
1705 !  x (u) advection
1706 
1707    i_start = its
1708    j_start = jts
1709    itf=MIN(ite,ide-1)
1710    jtf=MIN(jte,jde-1)
1711 
1712    IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+2)
1713    IF (config_flags%open_xe .or. specified ) itf     = min(itf,ide-3)
1714    IF ( config_flags%periodic_x ) i_start = its
1715    IF ( config_flags%periodic_x ) itf=MIN(ite,ide-1)
1716 
1717    DO j = j_start, jtf
1718 
1719      DO k = 2, kte-1
1720      DO i = i_start, itf
1721         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdx*(                     &
1722                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))/msfuy(i+1,j)               &
1723                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))/msfuy(i,j) )* (1./60.)*( &
1724                    45.*(ph(i+1,k,j)-ph(i-1,k,j))                      &
1725                    -9.*(ph(i+2,k,j)-ph(i-2,k,j))                      &
1726                       +(ph(i+3,k,j)-ph(i-3,k,j))                      &
1727                   +45.*(phb(i+1,k,j)-phb(i-1,k,j))                    &
1728                    -9.*(phb(i+2,k,j)-phb(i-2,k,j))                    &
1729                       +(phb(i+3,k,j)-phb(i-3,k,j))  )   )                
1730      ENDDO
1731      ENDDO
1732  
1733      k = kte
1734      DO i = i_start, itf
1735         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdx*(                                 &
1736                  ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))/msfuy(i+1,j)               &
1737                   +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))/msfuy(i,j) )* (1./60.)*(   &
1738                    45.*(ph(i+1,k,j)-ph(i-1,k,j))                                 &
1739                    -9.*(ph(i+2,k,j)-ph(i-2,k,j))                                 &
1740                       +(ph(i+3,k,j)-ph(i-3,k,j))                                 &
1741                   +45.*(phb(i+1,k,j)-phb(i-1,k,j))                               &
1742                    -9.*(phb(i+2,k,j)-phb(i-2,k,j))                               &
1743                       +(phb(i+3,k,j)-phb(i-3,k,j))  )     )
1744      ENDDO
1745 
1746    ENDDO
1747 
1748    IF ( (config_flags%open_xs) .and. its <= ids+1 ) THEN
1749      i = ids + 1
1750      DO j = j_start, jtf
1751        DO k = 2, kte-1
1752         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdx*(                     &
1753                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))/msfuy(i+1,j)               &
1754                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))/msfuy(i,j) )* (1./12.)*( &
1755                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                      &
1756                       -(ph(i+2,k,j)-ph(i-2,k,j))                      &
1757                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                    &
1758                       -(phb(i+2,k,j)-phb(i-2,k,j))  )   )                
1759        ENDDO
1760        k = kte
1761        ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdx*(                                 &
1762                 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))/msfuy(i+1,j)               &
1763                  +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))/msfuy(i,j) )* (1./12.)*(   &
1764                    8.*(ph(i+1,k,j)-ph(i-1,k,j))                                 &
1765                      -(ph(i+2,k,j)-ph(i-2,k,j))                                 &
1766                   +8.*(phb(i+1,k,j)-phb(i-1,k,j))                               &
1767                      -(phb(i+2,k,j)-phb(i-2,k,j))  )     )
1768 
1769      ENDDO
1770    END IF
1771 
1772    IF ( (config_flags%open_xe) .and. ite >= ide-2 ) THEN
1773      i = ide-2
1774      DO j = j_start, jtf
1775        DO k = 2, kte-1
1776         ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.25*rdx*(                     &
1777                  ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))/msfuy(i+1,j)               &
1778                   +muu(i,j  )*(u(i,k,j  )+u(i,k-1,j  ))/msfuy(i,j) )* (1./12.)*( &
1779                     8.*(ph(i+1,k,j)-ph(i-1,k,j))                      &
1780                       -(ph(i+2,k,j)-ph(i-2,k,j))                      &
1781                    +8.*(phb(i+1,k,j)-phb(i-1,k,j))                    &
1782                       -(phb(i+2,k,j)-phb(i-2,k,j))  )   )                
1783        ENDDO
1784        k = kte
1785        ph_tend(i,k,j)=ph_tend(i,k,j) - msftx(i,j)*.5*rdx*(                                 &
1786                 ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))/msfuy(i+1,j)               &
1787                  +muu(i,j  )*(cfn*u(i  ,k-1,j)+cfn1*u(i,k-2,j))/msfuy(i,j) )* (1./12.)*(   &
1788                    8.*(ph(i+1,k,j)-ph(i-1,k,j))                                 &
1789                      -(ph(i+2,k,j)-ph(i-2,k,j))                                 &
1790                   +8.*(phb(i+1,k,j)-phb(i-1,k,j))                               &
1791                      -(phb(i+2,k,j)-phb(i-2,k,j))  )     )
1792 
1793      ENDDO
1794    END IF
1795 
1796    END IF
1797 
1798 !  lateral open boundary conditions,
1799 !  start with north and south (y) boundaries
1800 
1801    i_start = its
1802    itf=MIN(ite,ide-1)
1803 
1804    !  south
1805 
1806    IF ( (config_flags%open_ys) .and. jts == jds ) THEN
1807    ! Notes on map scale factors:
1808    ! The following use only one 'mut' value, so maybe using one 1/mx value
1809    ! would be correct in which case it would just cancel with the *mx on top.
1810    ! So no change???
1811 
1812      j=jts
1813 
1814      DO k=2,kde
1815        kz = min(k,kde-1)
1816        DO i = its,itf
1817 !         vb =.5*( fnm(kz)*(v(i,kz  ,j+1)*msfvx_inv(i,j+1)+v(i,kz  ,j  )*msfvx_inv(i,j))    &
1818 !                 +fnp(kz)*(v(i,kz-1,j+1)*msfvx_inv(i,j+1)+v(i,kz-1,j  )*msfvx_inv(i,j)) )
1819          vb =.5*( fnm(kz)*(v(i,kz  ,j+1)+v(i,kz  ,j  ))    &
1820                  +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j  )) )
1821          vl=amin1(vb,0.)
1822 !         ph_tend(i,k,j)=ph_tend(i,k,j)-msftx(i,j)*rdy*mut(i,j)*(      &
1823          ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*(      &
1824                               +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
1825        ENDDO
1826      ENDDO
1827 
1828    END IF
1829 
1830    ! north
1831 
1832    IF ( (config_flags%open_ye) .and. jte == jde ) THEN
1833 
1834      j=jte-1
1835 
1836      DO k=2,kde
1837        kz = min(k,kde-1)
1838        DO i = its,itf
1839 !        vb=.5*( fnm(kz)*(v(i,kz  ,j+1)*msfvx_inv(i,j+1)+v(i,kz  ,j)*msfvx_inv(i,j))    &
1840 !               +fnp(kz)*(v(i,kz-1,j+1)*msfvx_inv(i,j+1)+v(i,kz-1,j)*msfvx_inv(i,j)) )
1841         vb=.5*( fnm(kz)*(v(i,kz  ,j+1)+v(i,kz  ,j))   &
1842                +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
1843         vr=amax1(vb,0.)
1844 !        ph_tend(i,k,j)=ph_tend(i,k,j)-msftx(i,j)*rdy*mut(i,j)*(      &
1845         ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*(      &
1846                    +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
1847        ENDDO
1848      ENDDO
1849 
1850    END IF
1851 
1852    !  now the east and west (y) boundaries
1853 
1854    j_start = its
1855    jtf=MIN(jte,jde-1)
1856 
1857    !  west
1858 
1859    IF ( (config_flags%open_xs) .and. its == ids ) THEN
1860    ! Notes on map scale factors
1861    ! The following use only one 'mut' value, so maybe using one 1/my value
1862    ! would be correct in which case we just need *(mx/my) on top???
1863 
1864      i=its
1865 
1866      DO j = jts,jtf
1867        DO k=2,kde-1
1868          kz = k
1869          ub =.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i  ,kz  ,j))     &
1870                  +fnp(kz)*(u(i+1,kz-1,j)+u(i  ,kz-1,j)) )
1871          ul=amin1(ub,0.)
1872          ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*(       &
1873                               +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1874        ENDDO
1875 
1876          k = kde
1877          kz = k
1878          ub =.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i  ,kz  ,j))     &
1879                  +fnp(kz)*(u(i+1,kz-1,j)+u(i  ,kz-1,j)) )
1880          ul=amin1(ub,0.)
1881          ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*(       &
1882                               +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
1883      ENDDO
1884 
1885    END IF
1886 
1887    ! east
1888 
1889    IF ( (config_flags%open_xe) .and. ite == ide ) THEN
1890 
1891      i = ite-1
1892 
1893      DO j = jts,jtf
1894        DO k=2,kde-1
1895         kz = k
1896         ub=.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i,kz  ,j))  &
1897                +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1898         ur=amax1(ub,0.)
1899         ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*( &
1900                    +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1901        ENDDO
1902 
1903         k = kde    
1904         kz = k-1
1905         ub=.5*( fnm(kz)*(u(i+1,kz  ,j)+u(i,kz  ,j))   &
1906                +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
1907         ur=amax1(ub,0.)
1908         ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*mut(i,j)*(  &
1909                    +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
1910 
1911      ENDDO
1912 
1913    END IF
1914 
1915   END SUBROUTINE rhs_ph
1916 
1917 !-------------------------------------------------------------------------------
1918 
1919 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend,                &
1920                                          ph,alt,p,pb,al,php,cqu,cqv,     &
1921                                          muu,muv,mu,fnm,fnp,rdnw,        &
1922                                          cf1,cf2,cf3,rdx,rdy,msfux,msfuy,&
1923                                          msfvx,msfvy,msftx,msfty,        &
1924                                          config_flags, non_hydrostatic,  &
1925                                          ids, ide, jds, jde, kds, kde,   &
1926                                          ims, ime, jms, jme, kms, kme,   &
1927                                          its, ite, jts, jte, kts, kte   )
1928 
1929    IMPLICIT NONE
1930    
1931    ! Input data
1932 
1933 
1934    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
1935 
1936    LOGICAL, INTENT (IN   ) :: non_hydrostatic
1937 
1938    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
1939                                        ims, ime, jms, jme, kms, kme, &
1940                                        its, ite, jts, jte, kts, kte 
1941 
1942    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::        &
1943                                                                      ph,  &
1944                                                                      alt, &
1945                                                                      al,  &
1946                                                                      p,   &
1947                                                                      pb,  &
1948                                                                      php, &
1949                                                                      cqu, &
1950                                                                      cqv
1951 
1952 
1953    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::           &
1954                                                                     ru_tend, &
1955                                                                     rv_tend
1956 
1957    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: muu, muv, mu,    &
1958                                                             msfux, msfuy, &
1959                                                             msfvx, msfvy, &
1960                                                             msftx, msfty
1961 
1962    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, fnm, fnp
1963 
1964    REAL,  INTENT(IN   ) :: rdx, rdy, cf1, cf2, cf3
1965 
1966    INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
1967    REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
1968    REAL :: dpx, dpy
1969 
1970    LOGICAL :: specified
1971 
1972 !<DESCRIPTION>
1973 !
1974 !  horizontal_pressure_gradient calculates the 
1975 !  horizontal pressure gradient terms for the large-timestep tendency 
1976 !  in the horizontal momentum equations (u,v).
1977 !
1978 !</DESCRIPTION>
1979 
1980    specified = .false.
1981    if(config_flags%specified .or. config_flags%nested) specified = .true.
1982 
1983 !  Notes on map scale factors:
1984 !  Calculates the pressure gradient terms in ADT eqns 44 and 45
1985 !  With upper rho -> 'mu', these are:
1986 !  Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
1987 !  Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
1988 !
1989 !  As we are on nu, rather than height, surfaces:
1990 !
1991 !  mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
1992 !           + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
1993 !
1994 !  mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
1995 !           + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
1996 
1997 ! start with the north-south (y) pressure gradient
1998 
1999    itf=MIN(ite,ide-1)
2000    jtf=jte
2001    ktf=MIN(kte,kde-1)
2002    i_start = its
2003    j_start = jts
2004    IF ( (config_flags%open_ys .or. specified .or. &
2005          config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2006    IF ( (config_flags%open_ye .or. specified .or. &
2007          config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2008 
2009    DO j = j_start, jtf
2010 
2011      IF ( non_hydrostatic )  THEN
2012 
2013         k=1
2014 
2015         DO i = i_start, itf
2016           dpn(i,k) = .5*( cf1*(p(i,k  ,j-1)+p(i,k  ,j))   &
2017                          +cf2*(p(i,k+1,j-1)+p(i,k+1,j))   &
2018                          +cf3*(p(i,k+2,j-1)+p(i,k+2,j))  )
2019           dpn(i,kde) = 0.
2020         ENDDO
2021                
2022         DO k=2,ktf
2023           DO i = i_start, itf
2024             dpn(i,k) = .5*( fnm(k)*(p(i,k  ,j-1)+p(i,k  ,j))  &
2025                            +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2026           END DO
2027         END DO
2028 
2029 !       ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2030 !       [alt, al are 1/rho terms; muv, mu are NOT coupled]
2031         DO K=1,ktf
2032           DO i = i_start, itf
2033             ! Here are mu dp/dy terms 1-3 
2034             dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*(                 &
2035                      (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1))  &
2036                     +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p (i,k,j)-p (i,k,j-1))  &
2037                     +(al (i,k  ,j)+al (i,k  ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2038             ! Here is mu dp/dy term 4 
2039             dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2040                 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j)))
2041             rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2042           END DO
2043         END DO
2044 
2045      ELSE
2046 
2047 !       ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2048 !       [alt, al are 1/rho terms; muv, mu are NOT coupled]
2049         DO K=1,ktf
2050           DO i = i_start, itf
2051             ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2052             dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*muv(i,j)*(                 &
2053                      (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1))  &
2054                     +(alt(i,k  ,j)+alt(i,k  ,j-1))*(p (i,k,j)-p (i,k,j-1))  &
2055                     +(al (i,k  ,j)+al (i,k  ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2056             rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2057           END DO
2058         END DO
2059 
2060      END IF
2061 
2062    ENDDO
2063 
2064 !  now the east-west (x) pressure gradient
2065 
2066    itf=ite
2067    jtf=MIN(jte,jde-1)
2068    ktf=MIN(kte,kde-1)
2069    i_start = its
2070    j_start = jts
2071    IF ( (config_flags%open_xs .or. specified .or. &
2072            config_flags%nested ) .and. its == ids ) i_start = its+1
2073    IF ( (config_flags%open_xe .or. specified .or. &
2074            config_flags%nested ) .and. ite == ide ) itf = itf-1
2075    IF ( config_flags%periodic_x ) i_start = its
2076    IF ( config_flags%periodic_x ) itf=ite
2077 
2078    DO j = j_start, jtf
2079 
2080      IF ( non_hydrostatic )  THEN
2081 
2082         k=1
2083 
2084         DO i = i_start, itf
2085           dpn(i,k) = .5*( cf1*(p(i-1,k  ,j)+p(i,k  ,j))   &
2086                          +cf2*(p(i-1,k+1,j)+p(i,k+1,j))   &
2087                          +cf3*(p(i-1,k+2,j)+p(i,k+2,j))  )
2088           dpn(i,kde) = 0.
2089         ENDDO
2090                
2091         DO k=2,ktf
2092           DO i = i_start, itf
2093             dpn(i,k) = .5*( fnm(k)*(p(i-1,k  ,j)+p(i,k  ,j))  &
2094                            +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2095           END DO
2096         END DO
2097 
2098 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2099 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2100         DO K=1,ktf
2101           DO i = i_start, itf
2102             ! Here are mu dp/dy terms 1-3
2103             dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*(                    &
2104                         (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j))  &
2105                        +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p (i,k,j)-p (i-1,k,j))  &
2106                        +(al (i,k  ,j)+al (i-1,k  ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2107             ! Here is mu dp/dy term 4
2108             dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2109                 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j)))
2110             ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2111           END DO
2112         END DO
2113 
2114      ELSE
2115 
2116 !       ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2117 !       [alt, al are 1/rho terms; muu, mu are NOT coupled]
2118         DO K=1,ktf
2119           DO i = i_start, itf
2120             ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2121             dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*muu(i,j)*(                    &
2122                         (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j))  &
2123                        +(alt(i,k  ,j)+alt(i-1,k  ,j))*(p (i,k,j)-p (i-1,k,j))  &
2124                        +(al (i,k  ,j)+al (i-1,k  ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2125             ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2126           END DO
2127         END DO
2128 
2129      END IF
2130 
2131    ENDDO
2132 
2133 END SUBROUTINE horizontal_pressure_gradient
2134 
2135 !-------------------------------------------------------------------------------
2136 
2137 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub,       &
2138                       rdnw, rdn, g, msftx, msfty,     &
2139                       ids, ide, jds, jde, kds, kde,   &
2140                       ims, ime, jms, jme, kms, kme,   &
2141                       its, ite, jts, jte, kts, kte   )
2142 
2143    IMPLICIT NONE
2144    
2145    ! Input data
2146 
2147    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2148                                        ims, ime, jms, jme, kms, kme, &
2149                                        its, ite, jts, jte, kts, kte 
2150 
2151    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   p
2152    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::   cqw
2153 
2154 
2155    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::  rw_tend
2156 
2157    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mub, mu, msftx, msfty
2158 
2159    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw, rdn
2160 
2161    REAL,  INTENT(IN   ) :: g
2162 
2163    INTEGER :: itf, jtf, i, j, k
2164    REAL    :: cq1, cq2
2165 
2166 
2167 !<DESCRIPTION>
2168 !
2169 !  pg_buoy_w calculates the 
2170 !  vertical pressure gradient and buoyancy terms for the large-timestep 
2171 !  tendency in the vertical momentum equation.
2172 !
2173 !</DESCRIPTION>
2174 
2175 !  BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2176 
2177 !  Map scale factor notes
2178 !  ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2179 !  Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2180 !  term 6: +(g/my) partial dp'/dnu
2181 !  term 7: -(g/my) mu'
2182 !
2183 !  For moisture-free atmosphere, cq1=1, cq2=0
2184 !  => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-mu(i,j)]
2185 
2186    itf=MIN(ite,ide-1)
2187    jtf=MIN(jte,jde-1)
2188 
2189    DO j = jts,jtf
2190 
2191      k=kde
2192      DO i=its,itf
2193        cq1 = 1./(1.+cqw(i,k-1,j))
2194        cq2 = cqw(i,k-1,j)*cq1
2195        rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*(      &
2196                         cq1*2.*rdnw(k-1)*(  -p(i,k-1,j))  &
2197                         -mu(i,j)-cq2*mub(i,j)            )
2198      END DO
2199 
2200      DO k = 2, kde-1
2201      DO i = its,itf
2202       cq1 = 1./(1.+cqw(i,k,j))
2203       cq2 = cqw(i,k,j)*cq1
2204       cqw(i,k,j) = cq1
2205       rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*(      &
2206                        cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j))  &
2207                        -mu(i,j)-cq2*mub(i,j)            )
2208      END DO
2209      ENDDO           
2210 
2211 
2212    ENDDO
2213 
2214 END SUBROUTINE pg_buoy_w
2215 
2216 !-------------------------------------------------------------------------------
2217 
2218 SUBROUTINE w_damp( rw_tend, ww, w, mut, rdnw, dt,     &
2219                       w_damping,                      &
2220                       ids, ide, jds, jde, kds, kde,   &
2221                       ims, ime, jms, jme, kms, kme,   &
2222                       its, ite, jts, jte, kts, kte   )
2223 
2224    IMPLICIT NONE
2225 
2226    ! Input data
2227 
2228    INTEGER ,          INTENT(IN   ) :: w_damping
2229 
2230    INTEGER ,          INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2231                                        ims, ime, jms, jme, kms, kme, &
2232                                        its, ite, jts, jte, kts, kte
2233 
2234    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN   ) ::   ww, w
2235 
2236    REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) ::  rw_tend
2237 
2238    REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN   ) :: mut
2239 
2240    REAL, DIMENSION( kms:kme ), INTENT(IN   ) :: rdnw
2241 
2242    REAL, INTENT(IN)    :: dt
2243    REAL                :: cfl, cf_n, cf_d, maxcfl, maxdub, maxdeta
2244 
2245    INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2246    INTEGER :: some
2247    CHARACTER*512 :: temp
2248    CHARACTER (LEN=256) :: time_str
2249    CHARACTER (LEN=256) :: grid_str
2250 
2251 !<DESCRIPTION>
2252 !
2253 !  w_damp computes a damping term for the vertical velocity when the
2254 !  vertical Courant number is too large.  This was found to be preferable to 
2255 !  decreasing the timestep or increasing the diffusion in real-data applications
2256 !  that produced potentially-unstable large vertical velocities because of
2257 !  unphysically large heating rates coming from the cumulus parameterization 
2258 !  schemes run at moderately high resolutions (dx ~ O(10) km).
2259 !
2260 !</DESCRIPTION>
2261 
2262    itf=MIN(ite,ide-1)
2263    jtf=MIN(jte,jde-1)
2264 
2265    some = 0
2266    maxcfl = 0.
2267 
2268    IF ( w_damping == 1 ) THEN
2269      DO j = jts,jtf
2270 
2271      DO k = 2, kde-1
2272      DO i = its,itf
2273 #if 0
2274         cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2275         if(cfl .gt. w_beta)then
2276 #else
2277 ! restructure to get rid of divide
2278         cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2279         cf_d = abs(mut(i,j))
2280         if(cf_n .gt. cf_d*w_beta )then
2281 #endif
2282            cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2283            IF ( cfl > maxcfl ) THEN
2284              maxcfl = cfl ; maxi = i ; maxj = j ; maxk = k 
2285              maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2286            ENDIF
2287            WRITE(temp,*)i,j,k,' cfl,w,d(eta)=',cfl,w(i,k,j),-1./rdnw(k)
2288            CALL wrf_debug ( 100 , TRIM(temp) )
2289            if ( cfl > 2. ) some = some + 1
2290            rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(cfl-w_beta)*mut(i,j)
2291         endif
2292      END DO
2293      ENDDO
2294      ENDDO
2295    ELSE
2296 ! just print
2297      DO j = jts,jtf
2298 
2299      DO k = 2, kde-1
2300      DO i = its,itf
2301         cf_n = abs(ww(i,k,j)*rdnw(k)*dt)
2302         cf_d = abs(mut(i,j))
2303         if(cf_n .gt. cf_d*w_beta )then
2304            cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt)
2305            IF ( cfl > maxcfl ) THEN
2306              maxcfl = cfl ; maxi = i ; maxj = j ; maxk = k 
2307              maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2308            ENDIF
2309            WRITE(temp,*)i,j,k,' cfl,w,d(eta)=',cfl,w(i,k,j),-1./rdnw(k)
2310            CALL wrf_debug ( 100 , TRIM(temp) )
2311            if ( cfl > 2. ) some = some + 1
2312         endif
2313      END DO
2314      ENDDO
2315      ENDDO
2316    ENDIF
2317    IF ( some .GT. 0 ) THEN
2318      CALL get_current_time_string( time_str )
2319      CALL get_current_grid_name( grid_str )
2320      WRITE(wrf_err_message,*)some,                                            &
2321             ' points exceeded cfl=2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2322      CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2323      WRITE(wrf_err_message,*)'MAX AT i,j,k: ',maxi,maxj,maxk,' cfl,w,d(eta)=',maxcfl, &
2324                              maxdub,maxdeta
2325      CALL wrf_debug ( 0 , TRIM(wrf_err_message) )
2326    ENDIF
2327 
2328 END SUBROUTINE w_damp
2329 
2330 !-------------------------------------------------------------------------------
2331 
2332 SUBROUTINE horizontal_diffusion ( name, field, tendency, mu,           &
2333                                   config_flags,                        &
2334                                   msfux, msfuy, msfvx, msfvx_inv,      &
2335                                   msfvy, msftx, msfty,                 &
2336                                   khdif, xkmhd, rdx, rdy,              &
2337                                   ids, ide, jds, jde, kds, kde,        &
2338                                   ims, ime, jms, jme, kms, kme,        &
2339                                   its, ite, jts, jte, kts, kte        )
2340 
2341    IMPLICIT NONE
2342    
2343    ! Input data
2344 
2345    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2346 
2347    INTEGER ,        INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2348                                      ims, ime, jms, jme, kms, kme, &
2349                                      its, ite, jts, jte, kts, kte
2350 
2351    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
2352 
2353    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field, xkmhd
2354 
2355    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2356 
2357    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
2358 
2359    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
2360                                                                     msfuy,      &
2361                                                                     msfvx,      &
2362                                                                     msfvx_inv,  &
2363                                                                     msfvy,      &
2364                                                                     msftx,      &
2365                                                                     msfty
2366 
2367    REAL ,                                      INTENT(IN   ) :: rdx,       &
2368                                                                 rdy,       &
2369                                                                 khdif
2370 
2371    ! Local data
2372    
2373    INTEGER :: i, j, k, itf, jtf, ktf
2374 
2375    INTEGER :: i_start, i_end, j_start, j_end
2376 
2377    REAL :: mrdx, mkrdxm, mkrdxp, &
2378            mrdy, mkrdym, mkrdyp
2379 
2380    LOGICAL :: specified
2381 
2382 !<DESCRIPTION>
2383 !
2384 !  horizontal_diffusion computes the horizontal diffusion tendency
2385 !  on model horizontal coordinate surfaces.
2386 !
2387 !</DESCRIPTION>
2388 
2389    specified = .false.
2390    if(config_flags%specified .or. config_flags%nested) specified = .true.
2391 
2392    ktf=MIN(kte,kde-1)
2393    
2394    IF (name .EQ. 'u') THEN
2395 
2396       i_start = its
2397       i_end   = ite
2398       j_start = jts
2399       j_end   = MIN(jte,jde-1)
2400 
2401       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2402       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-1,ite)
2403       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2404       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2405       IF ( config_flags%periodic_x ) i_start = its
2406       IF ( config_flags%periodic_x ) i_end = ite
2407 
2408 
2409       DO j = j_start, j_end
2410       DO k=kts,ktf
2411       DO i = i_start, i_end
2412 
2413          ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2414          ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2415 
2416          mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*mu(i-1,j)*xkmhd(i-1,k,j)*rdx
2417          mkrdxp=(msftx(i,j)/msfty(i,j))*mu(i,j)*xkmhd(i,k,j)*rdx
2418          mrdx=msfux(i,j)*msfuy(i,j)*rdx 
2419          mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2420                 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2421                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2422          mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2423                 0.25*(mu(i,j)+mu(i,j+1)+mu(i-1,j+1)+mu(i-1,j))* &
2424                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2425          ! need to do four-corners (t) for diffusion coefficient as there are
2426          ! no values at u,v points
2427          ! msfuy - has to be y as part of d/dY
2428          !         has to be u as we're at a u point
2429          mrdy=msfux(i,j)*msfuy(i,j)*rdy 
2430 
2431          ! correctly averaged version of rho~ * m^2 * 
2432          !    [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2433             tendency(i,k,j)=tendency(i,k,j)+( &
2434                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j))  &
2435                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2436                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  ))  &
2437                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2438       ENDDO
2439       ENDDO
2440       ENDDO
2441    
2442    ELSE IF (name .EQ. 'v')THEN
2443 
2444       i_start = its
2445       i_end   = MIN(ite,ide-1)
2446       j_start = jts
2447       j_end   = jte
2448 
2449       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2450       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2451       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2452       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-1,jte)
2453       IF ( config_flags%periodic_x ) i_start = its
2454       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2455       IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2456       IF ( config_flags%polar ) j_end   = MIN(jde-1,jte)
2457 
2458       DO j = j_start, j_end
2459       DO k=kts,ktf
2460       DO i = i_start, i_end
2461 
2462          mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )*    &
2463                 0.25*(mu(i,j)+mu(i,j-1)+mu(i-1,j-1)+mu(i-1,j))* &
2464                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2465          mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )*    &
2466                 0.25*(mu(i,j)+mu(i,j-1)+mu(i+1,j-1)+mu(i+1,j))* &
2467                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2468          mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2469          mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2470          mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2471          mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2472 
2473             tendency(i,k,j)=tendency(i,k,j)+( &
2474                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j))  &
2475                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2476                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  ))  &
2477                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2478       ENDDO
2479       ENDDO
2480       ENDDO
2481    
2482    ELSE IF (name .EQ. 'w')THEN
2483 
2484       i_start = its
2485       i_end   = MIN(ite,ide-1)
2486       j_start = jts
2487       j_end   = MIN(jte,jde-1)
2488 
2489       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2490       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2491       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2492       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2493       IF ( config_flags%periodic_x ) i_start = its
2494       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2495 
2496       DO j = j_start, j_end
2497       DO k=kts+1,ktf
2498       DO i = i_start, i_end
2499 
2500          mkrdxm=(msfux(i,j)/msfuy(i,j))*   &
2501                 0.25*(mu(i,j)+mu(i-1,j)+mu(i,j)+mu(i-1,j))* &
2502                 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2503          mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*   &
2504                 0.25*(mu(i+1,j)+mu(i,j)+mu(i+1,j)+mu(i,j))* &
2505                 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2506          mrdx=msftx(i,j)*msfty(i,j)*rdx
2507 !         mkrdym=(msfvy(i,j)/msfvx(i,j))*   &
2508          mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*   &
2509                 0.25*(mu(i,j)+mu(i,j-1)+mu(i,j)+mu(i,j-1))* &
2510                 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2511 !         mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*   &
2512          mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*   &
2513                 0.25*(mu(i,j+1)+mu(i,j)+mu(i,j+1)+mu(i,j))* &
2514                 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2515          mrdy=msftx(i,j)*msfty(i,j)*rdy
2516 
2517             tendency(i,k,j)=tendency(i,k,j)+( &
2518                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j)) &
2519                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2520                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  )) &
2521                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2522       ENDDO
2523       ENDDO
2524       ENDDO
2525    
2526    ELSE
2527 
2528 
2529       i_start = its
2530       i_end   = MIN(ite,ide-1)
2531       j_start = jts
2532       j_end   = MIN(jte,jde-1)
2533 
2534       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2535       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2536       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2537       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2538       IF ( config_flags%periodic_x ) i_start = its
2539       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2540 
2541       DO j = j_start, j_end
2542       DO k=kts,ktf
2543       DO i = i_start, i_end
2544 
2545          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
2546          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
2547          mrdx=msftx(i,j)*msfty(i,j)*rdx
2548 !         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
2549          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
2550 !         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
2551          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
2552          mrdy=msftx(i,j)*msfty(i,j)*rdy
2553 
2554             tendency(i,k,j)=tendency(i,k,j)+( &
2555                             mrdx*(mkrdxp*(field(i+1,k,j)-field(i  ,k,j))  &
2556                                  -mkrdxm*(field(i  ,k,j)-field(i-1,k,j))) &
2557                            +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j  ))  &
2558                                  -mkrdym*(field(i,k,j  )-field(i,k,j-1))))
2559       ENDDO
2560       ENDDO
2561       ENDDO
2562            
2563    ENDIF
2564 
2565 END SUBROUTINE horizontal_diffusion
2566 
2567 !-----------------------------------------------------------------------------------------
2568 
2569 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu,           &
2570                                        config_flags, base_3d,               &
2571                                        msfux, msfuy, msfvx, msfvx_inv,      &
2572                                        msfvy, msftx, msfty,                 &
2573                                        khdif, xkmhd, rdx, rdy,              &
2574                                        ids, ide, jds, jde, kds, kde,        &
2575                                        ims, ime, jms, jme, kms, kme,        &
2576                                        its, ite, jts, jte, kts, kte        )
2577 
2578    IMPLICIT NONE
2579    
2580    ! Input data
2581    
2582    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2583 
2584    INTEGER ,        INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2585                                      ims, ime, jms, jme, kms, kme, &
2586                                      its, ite, jts, jte, kts, kte
2587 
2588    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
2589 
2590    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: field, &
2591                                                                       xkmhd, &
2592                                                                       base_3d
2593 
2594    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2595 
2596    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN   ) :: mu
2597 
2598    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
2599                                                                     msfuy,      &
2600                                                                     msfvx,      &
2601                                                                     msfvx_inv,  &
2602                                                                     msfvy,      &
2603                                                                     msftx,      &
2604                                                                     msfty
2605 
2606    REAL ,                                      INTENT(IN   ) :: rdx,       &
2607                                                                 rdy,       &
2608                                                                 khdif
2609 
2610    ! Local data
2611    
2612    INTEGER :: i, j, k, itf, jtf, ktf
2613 
2614    INTEGER :: i_start, i_end, j_start, j_end
2615 
2616    REAL :: mrdx, mkrdxm, mkrdxp, &
2617            mrdy, mkrdym, mkrdyp
2618 
2619    LOGICAL :: specified
2620 
2621 !<DESCRIPTION>
2622 !
2623 !  horizontal_diffusion_3dmp computes the horizontal diffusion tendency
2624 !  on model horizontal coordinate surfaces.  This routine computes diffusion
2625 !  a perturbation scalar (field-base_3d).
2626 !
2627 !</DESCRIPTION>
2628 
2629    specified = .false.
2630    if(config_flags%specified .or. config_flags%nested) specified = .true.
2631 
2632    ktf=MIN(kte,kde-1)
2633    
2634       i_start = its
2635       i_end   = MIN(ite,ide-1)
2636       j_start = jts
2637       j_end   = MIN(jte,jde-1)
2638 
2639       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2640       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-2,ite)
2641       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2642       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-2,jte)
2643       IF ( config_flags%periodic_x ) i_start = its
2644       IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2645 
2646       DO j = j_start, j_end
2647       DO k=kts,ktf
2648       DO i = i_start, i_end
2649 
2650          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
2651          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
2652          mrdx=msftx(i,j)*msfty(i,j)*rdx
2653 !         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
2654 !         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
2655          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
2656          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
2657          mrdy=msftx(i,j)*msfty(i,j)*rdy
2658 
2659             tendency(i,k,j)=tendency(i,k,j)+(                        &
2660                     mrdx*( mkrdxp*(   field(i+1,k,j)  -field(i  ,k,j)      &
2661                                    -base_3d(i+1,k,j)+base_3d(i  ,k,j) )    &
2662                           -mkrdxm*(   field(i  ,k,j)  -field(i-1,k,j)      &
2663                                    -base_3d(i  ,k,j)+base_3d(i-1,k,j) )  ) &
2664                    +mrdy*( mkrdyp*(   field(i,k,j+1)  -field(i,k,j  )      &
2665                                    -base_3d(i,k,j+1)+base_3d(i,k,j  ) )    &
2666                           -mkrdym*(   field(i,k,j  )  -field(i,k,j-1)      &
2667                                    -base_3d(i,k,j  )+base_3d(i,k,j-1) )  ) &
2668                                                                          ) 
2669       ENDDO
2670       ENDDO
2671       ENDDO
2672 
2673 END SUBROUTINE horizontal_diffusion_3dmp
2674 
2675 !-----------------------------------------------------------------------------------------
2676 
2677 SUBROUTINE vertical_diffusion ( name, field, tendency,        &
2678                                 config_flags,                 &
2679                                 alt, mut, rdn, rdnw, kvdif,   &
2680                                 ids, ide, jds, jde, kds, kde, &
2681                                 ims, ime, jms, jme, kms, kme, &
2682                                 its, ite, jts, jte, kts, kte )
2683 
2684 
2685    IMPLICIT NONE
2686    
2687    ! Input data
2688    
2689    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2690 
2691    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2692                                  ims, ime, jms, jme, kms, kme, &
2693                                  its, ite, jts, jte, kts, kte
2694 
2695    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
2696 
2697    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
2698                                                INTENT(IN   ) :: field,    &
2699                                                                 alt
2700 
2701    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2702 
2703    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut
2704 
2705    REAL , DIMENSION( kms:kme ) ,                   INTENT(IN   ) :: rdn, rdnw
2706 
2707    REAL ,                                      INTENT(IN   ) :: kvdif
2708    
2709    ! Local data
2710    
2711    INTEGER :: i, j, k, itf, jtf, ktf
2712    INTEGER :: i_start, i_end, j_start, j_end
2713 
2714    REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
2715    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2716 
2717    REAL :: rdz
2718 
2719    LOGICAL :: specified
2720 
2721 !<DESCRIPTION>
2722 !
2723 !  vertical_diffusion
2724 !  computes vertical diffusion tendency.
2725 !
2726 !</DESCRIPTION>
2727 
2728    specified = .false.
2729    if(config_flags%specified .or. config_flags%nested) specified = .true.
2730 
2731    ktf=MIN(kte,kde-1)
2732    
2733    IF (name .EQ. 'w')THEN
2734 
2735    
2736    i_start = its
2737    i_end   = MIN(ite,ide-1)
2738    j_start = jts
2739    j_end   = MIN(jte,jde-1)
2740 
2741 j_loop_w : DO j = j_start, j_end
2742 
2743      DO k=kts,ktf-1
2744        DO i = i_start, i_end
2745           vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
2746        ENDDO
2747      ENDDO
2748 
2749      DO i = i_start, i_end
2750        vflux(i,ktf)=0.
2751      ENDDO
2752 
2753      DO k=kts+1,ktf
2754        DO i = i_start, i_end
2755             tendency(i,k,j)=tendency(i,k,j)                                         &
2756                               +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j)))  &
2757                                          *(vflux(i,k)-vflux(i,k-1))
2758        ENDDO
2759      ENDDO
2760 
2761     ENDDO j_loop_w
2762 
2763    ELSE IF(name .EQ. 'm')THEN
2764 
2765      i_start = its
2766      i_end   = MIN(ite,ide-1)
2767      j_start = jts
2768      j_end   = MIN(jte,jde-1)
2769 
2770 j_loop_s : DO j = j_start, j_end
2771 
2772      DO k=kts,ktf-1
2773        DO i = i_start, i_end
2774          vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j)))   &
2775                   *(field(i,k+1,j)-field(i,k,j))
2776        ENDDO
2777      ENDDO
2778 
2779      DO i = i_start, i_end
2780        vflux(i,0)=vflux(i,1)
2781      ENDDO
2782 
2783      DO i = i_start, i_end
2784        vflux(i,ktf)=0.
2785      ENDDO
2786 
2787      DO k=kts,ktf
2788        DO i = i_start, i_end
2789          tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j)  &
2790                 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
2791        ENDDO
2792      ENDDO
2793 
2794  ENDDO j_loop_s
2795 
2796    ENDIF
2797 
2798 END SUBROUTINE vertical_diffusion
2799 
2800 
2801 !-------------------------------------------------------------------------------
2802 
2803 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
2804                                    base,                          &
2805                                    alt, mut, rdn, rdnw, kvdif,    &
2806                                    ids, ide, jds, jde, kds, kde,  &
2807                                    ims, ime, jms, jme, kms, kme,  &
2808                                    its, ite, jts, jte, kts, kte  )
2809 
2810 
2811    IMPLICIT NONE
2812    
2813    ! Input data
2814    
2815    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2816 
2817    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2818                                  ims, ime, jms, jme, kms, kme, &
2819                                  its, ite, jts, jte, kts, kte
2820 
2821    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
2822                                                INTENT(IN   ) :: field,    &
2823                                                                 alt
2824 
2825    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2826 
2827    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut
2828 
2829    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn,  &
2830                                                                   rdnw, &
2831                                                                   base
2832 
2833    REAL ,                                      INTENT(IN   ) :: kvdif
2834    
2835    ! Local data
2836    
2837    INTEGER :: i, j, k, itf, jtf, ktf
2838    INTEGER :: i_start, i_end, j_start, j_end
2839 
2840    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2841 
2842    REAL :: rdz
2843 
2844    LOGICAL :: specified
2845 
2846 !<DESCRIPTION>
2847 !
2848 !  vertical_diffusion_mp
2849 !  computes vertical diffusion tendency of a perturbation variable
2850 !  (field-base).  Note that base as a 1D (k) field.
2851 !
2852 !</DESCRIPTION>
2853 
2854    specified = .false.
2855    if(config_flags%specified .or. config_flags%nested) specified = .true.
2856 
2857    ktf=MIN(kte,kde-1)
2858    
2859      i_start = its
2860      i_end   = MIN(ite,ide-1)
2861      j_start = jts
2862      j_end   = MIN(jte,jde-1)
2863 
2864 j_loop_s : DO j = j_start, j_end
2865 
2866      DO k=kts,ktf-1
2867        DO i = i_start, i_end
2868          vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j)))   &
2869                     *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
2870        ENDDO
2871      ENDDO
2872 
2873      DO i = i_start, i_end
2874        vflux(i,0)=vflux(i,1)
2875      ENDDO
2876 
2877      DO i = i_start, i_end
2878        vflux(i,ktf)=0.
2879      ENDDO
2880 
2881      DO k=kts,ktf
2882        DO i = i_start, i_end
2883          tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j)  &
2884                 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
2885        ENDDO
2886      ENDDO
2887 
2888  ENDDO j_loop_s
2889 
2890 END SUBROUTINE vertical_diffusion_mp
2891 
2892 
2893 !-------------------------------------------------------------------------------
2894 
2895 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
2896                                      base_3d,                       &
2897                                      alt, mut, rdn, rdnw, kvdif,    &
2898                                      ids, ide, jds, jde, kds, kde,  &
2899                                      ims, ime, jms, jme, kms, kme,  &
2900                                      its, ite, jts, jte, kts, kte  )
2901 
2902 
2903    IMPLICIT NONE
2904    
2905    ! Input data
2906    
2907    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
2908 
2909    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
2910                                  ims, ime, jms, jme, kms, kme, &
2911                                  its, ite, jts, jte, kts, kte
2912 
2913    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
2914                                                INTENT(IN   ) :: field,    &
2915                                                                 alt,      &
2916                                                                 base_3d
2917 
2918    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2919 
2920    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: mut
2921 
2922    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn,  &
2923                                                                   rdnw
2924 
2925    REAL ,                                      INTENT(IN   ) :: kvdif
2926    
2927    ! Local data
2928    
2929    INTEGER :: i, j, k, itf, jtf, ktf
2930    INTEGER :: i_start, i_end, j_start, j_end
2931 
2932    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
2933 
2934    REAL :: rdz
2935 
2936    LOGICAL :: specified
2937 
2938 !<DESCRIPTION>
2939 !
2940 !  vertical_diffusion_3dmp
2941 !  computes vertical diffusion tendency of a perturbation variable
2942 !  (field-base_3d).  
2943 !
2944 !</DESCRIPTION>
2945 
2946    specified = .false.
2947    if(config_flags%specified .or. config_flags%nested) specified = .true.
2948 
2949    ktf=MIN(kte,kde-1)
2950    
2951      i_start = its
2952      i_end   = MIN(ite,ide-1)
2953      j_start = jts
2954      j_end   = MIN(jte,jde-1)
2955 
2956 j_loop_s : DO j = j_start, j_end
2957 
2958      DO k=kts,ktf-1
2959        DO i = i_start, i_end
2960          vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j)))   &
2961                     *(   field(i,k+1,j)  -field(i,k,j)               &
2962                       -base_3d(i,k+1,j)+base_3d(i,k,j) )
2963        ENDDO
2964      ENDDO
2965 
2966      DO i = i_start, i_end
2967        vflux(i,0)=vflux(i,1)
2968      ENDDO
2969 
2970      DO i = i_start, i_end
2971        vflux(i,ktf)=0.
2972      ENDDO
2973 
2974      DO k=kts,ktf
2975        DO i = i_start, i_end
2976          tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j)  &
2977                 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
2978        ENDDO
2979      ENDDO
2980 
2981  ENDDO j_loop_s
2982 
2983 END SUBROUTINE vertical_diffusion_3dmp
2984 
2985 
2986 !-------------------------------------------------------------------------------
2987 
2988 
2989 SUBROUTINE vertical_diffusion_u ( field, tendency,              &
2990                                   config_flags, u_base,         &
2991                                   alt, muu, rdn, rdnw, kvdif,   &
2992                                   ids, ide, jds, jde, kds, kde, &
2993                                   ims, ime, jms, jme, kms, kme, &
2994                                   its, ite, jts, jte, kts, kte )
2995 
2996 
2997    IMPLICIT NONE
2998    
2999    ! Input data
3000    
3001    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
3002 
3003    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3004                                  ims, ime, jms, jme, kms, kme, &
3005                                  its, ite, jts, jte, kts, kte
3006 
3007    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
3008                                                INTENT(IN   ) :: field,    &
3009                                                                 alt
3010 
3011    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3012 
3013    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: muu
3014 
3015    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn, rdnw, u_base
3016 
3017    REAL ,                                      INTENT(IN   ) :: kvdif
3018    
3019    ! Local data
3020    
3021    INTEGER :: i, j, k, itf, jtf, ktf
3022    INTEGER :: i_start, i_end, j_start, j_end
3023 
3024    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3025 
3026    REAL :: rdz, zz
3027 
3028    LOGICAL :: specified
3029 
3030 !<DESCRIPTION>
3031 !
3032 !  vertical_diffusion_u computes vertical diffusion tendency for 
3033 !  the u momentum equation.  This routine assumes a constant eddy
3034 !  viscosity kvdif.
3035 !
3036 !</DESCRIPTION>
3037 
3038    specified = .false.
3039    if(config_flags%specified .or. config_flags%nested) specified = .true.
3040 
3041    ktf=MIN(kte,kde-1)
3042 
3043       i_start = its
3044       i_end   = ite
3045       j_start = jts
3046       j_end   = MIN(jte,jde-1)
3047 
3048       IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3049       IF ( config_flags%open_xe .or. specified ) i_end   = MIN(ide-1,ite)
3050       IF ( config_flags%periodic_x ) i_start = its
3051       IF ( config_flags%periodic_x ) i_end = ite
3052 
3053 
3054 j_loop_u : DO j = j_start, j_end
3055 
3056      DO k=kts,ktf-1
3057        DO i = i_start, i_end
3058          vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i  ,k  ,j)      &
3059                                         +alt(i-1,k  ,j)      &
3060                                         +alt(i  ,k+1,j)      &
3061                                         +alt(i-1,k+1,j) ) )  &
3062                              *(field(i,k+1,j)-field(i,k,j)   &
3063                                -u_base(k+1)   +u_base(k)  )
3064        ENDDO
3065      ENDDO
3066 
3067      DO i = i_start, i_end
3068        vflux(i,0)=vflux(i,1)
3069      ENDDO
3070 
3071      DO i = i_start, i_end
3072        vflux(i,ktf)=0.
3073      ENDDO
3074 
3075      DO k=kts,ktf-1
3076        DO i = i_start, i_end
3077          tendency(i,k,j)=tendency(i,k,j)+                             &
3078                 g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3079                               (vflux(i,k)-vflux(i,k-1))
3080        ENDDO
3081      ENDDO
3082 
3083  ENDDO j_loop_u
3084    
3085 END SUBROUTINE vertical_diffusion_u
3086 
3087 !-------------------------------------------------------------------------------
3088 
3089 
3090 SUBROUTINE vertical_diffusion_v ( field, tendency,              &
3091                                   config_flags, v_base,         &
3092                                   alt, muv, rdn, rdnw, kvdif,   &
3093                                   ids, ide, jds, jde, kds, kde, &
3094                                   ims, ime, jms, jme, kms, kme, &
3095                                   its, ite, jts, jte, kts, kte )
3096 
3097 
3098    IMPLICIT NONE
3099    
3100    ! Input data
3101    
3102    TYPE(grid_config_rec_type), INTENT(IN   ) :: config_flags
3103 
3104    INTEGER ,    INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3105                                  ims, ime, jms, jme, kms, kme, &
3106                                  its, ite, jts, jte, kts, kte
3107 
3108    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                      &
3109                                                INTENT(IN   ) :: field,    &
3110                                                                 alt
3111    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: rdn, rdnw, v_base
3112 
3113    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3114 
3115    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: muv
3116 
3117    REAL ,                                      INTENT(IN   ) :: kvdif
3118    
3119    ! Local data
3120    
3121    INTEGER :: i, j, k, itf, jtf, ktf, jm1
3122    INTEGER :: i_start, i_end, j_start, j_end
3123 
3124    REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3125 
3126    REAL :: rdz, zz
3127 
3128    LOGICAL :: specified
3129 
3130 !<DESCRIPTION>
3131 !
3132 !  vertical_diffusion_v computes vertical diffusion tendency for 
3133 !  the v momentum equation.  This routine assumes a constant eddy
3134 !  viscosity kvdif.
3135 !
3136 !</DESCRIPTION>
3137 
3138    specified = .false.
3139    if(config_flags%specified .or. config_flags%nested) specified = .true.
3140 
3141    ktf=MIN(kte,kde-1)
3142    
3143       i_start = its
3144       i_end   = MIN(ite,ide-1)
3145       j_start = jts
3146       j_end   = MIN(jte,jde-1)
3147 
3148       IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3149       IF ( config_flags%open_ye .or. specified ) j_end   = MIN(jde-1,jte)
3150 
3151 j_loop_v : DO j = j_start, j_end
3152 !     jm1 = max(j-1,1)
3153      jm1 = j-1
3154 
3155      DO k=kts,ktf-1
3156        DO i = i_start, i_end
3157          vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k  ,j  )      &
3158                                         +alt(i,k  ,jm1)      &
3159                                         +alt(i,k+1,j  )      &
3160                                         +alt(i,k+1,jm1) ) )  &
3161                              *(field(i,k+1,j)-field(i,k,j)   &
3162                                -v_base(k+1)   +v_base(k)  )
3163        ENDDO
3164      ENDDO
3165 
3166      DO i = i_start, i_end
3167        vflux(i,0)=vflux(i,1)
3168      ENDDO
3169 
3170      DO i = i_start, i_end
3171        vflux(i,ktf)=0.
3172      ENDDO
3173 
3174      DO k=kts,ktf-1
3175        DO i = i_start, i_end 
3176          tendency(i,k,j)=tendency(i,k,j)+                              &
3177                 g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))*  &
3178                               (vflux(i,k)-vflux(i,k-1))
3179        ENDDO
3180      ENDDO
3181 
3182  ENDDO j_loop_v
3183    
3184 END SUBROUTINE vertical_diffusion_v
3185 
3186 !***************  end new mass coordinate routines
3187 
3188 !-------------------------------------------------------------------------------
3189 
3190 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp,     &
3191                             ids, ide, jds, jde, kds, kde, &
3192                             ims, ime, jms, jme, kms, kme, &
3193                             its, ite, jts, jte, kts, kte )
3194 
3195    IMPLICIT NONE
3196    
3197    ! Input data
3198    
3199    INTEGER ,      INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3200                                    ims, ime, jms, jme, kms, kme, &
3201                                    its, ite, jts, jte, kts, kte 
3202    
3203    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: rfieldb, &
3204                                                                       rfieldp
3205 
3206    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT  ) :: rfield
3207    
3208    ! Local indices.
3209    
3210    INTEGER :: i, j, k, itf, jtf, ktf
3211    
3212 !<DESCRIPTION>
3213 !
3214 !  calculate_full
3215 !  calculates full 3D field from pertubation and base field.
3216 !
3217 !</DESCRIPTION>
3218 
3219    itf=MIN(ite,ide-1)
3220    jtf=MIN(jte,jde-1)
3221    ktf=MIN(kte,kde-1)
3222 
3223    DO j=jts,jtf
3224    DO k=kts,ktf
3225    DO i=its,itf
3226       rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3227    ENDDO
3228    ENDDO
3229    ENDDO
3230 
3231 END SUBROUTINE calculate_full
3232 
3233 !------------------------------------------------------------------------------
3234 
3235 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3236                       config_flags,                          &
3237                       msftx, msfty, msfux, msfuy,            &
3238                       msfvx, msfvy,                          &
3239                       f, e, sina, cosa, fzm, fzp,            &
3240                       ids, ide, jds, jde, kds, kde,          &
3241                       ims, ime, jms, jme, kms, kme,          &
3242                       its, ite, jts, jte, kts, kte          )
3243 
3244    IMPLICIT NONE
3245    
3246    ! Input data
3247    
3248    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
3249 
3250    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3251                                               ims, ime, jms, jme, kms, kme, &
3252                                               its, ite, jts, jte, kts, kte
3253 
3254    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3255                                                                 rv_tend, &
3256                                                                 rw_tend
3257    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: ru, &
3258                                                                 rv, &
3259                                                                 rw
3260 
3261    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
3262                                                                 msfuy,      &
3263                                                                 msfvx,      &
3264                                                                 msfvy,      &
3265                                                                 msftx,      &
3266                                                                 msfty
3267 
3268    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: f,    &
3269                                                                     e,    &
3270                                                                     sina, &
3271                                                                     cosa
3272 
3273    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm, &
3274                                                                   fzp
3275    
3276    ! Local indices.
3277    
3278    INTEGER :: i, j , k, ktf
3279    INTEGER :: i_start, i_end, j_start, j_end
3280    
3281    LOGICAL :: specified
3282 
3283 !<DESCRIPTION>
3284 !
3285 !  coriolis calculates the large timestep tendency terms in the 
3286 !  u, v, and w momentum equations arise from the coriolis force.
3287 !
3288 !</DESCRIPTION>
3289 
3290    specified = .false.
3291    if(config_flags%specified .or. config_flags%nested) specified = .true.
3292 
3293    ktf=MIN(kte,kde-1)
3294 
3295 ! coriolis for u-momentum equation
3296 
3297 !  Notes on map scale factor
3298 !  cosa, sina are related to rotating the coordinate frame if desired
3299 !  generally sina=0, cosa=1
3300 !  ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3301 !                                + 2 mu v omega sin(lat)/my
3302 !  Define f=2 omega sin(lat), e=2 omega cos(lat)
3303 !   => terms are: -e mu w / my + f mu v / my
3304 !  rv = mu v / mx ; rw = mu w / my
3305 !   => terms are: -e rw + f rv *mx / my
3306 
3307    i_start = its
3308    i_end   = ite
3309    IF ( config_flags%open_xs .or. specified .or. &
3310         config_flags%nested) i_start = MAX(ids+1,its)
3311    IF ( config_flags%open_xe .or. specified .or. &
3312         config_flags%nested) i_end   = MIN(ide-1,ite)
3313       IF ( config_flags%periodic_x ) i_start = its
3314       IF ( config_flags%periodic_x ) i_end = ite
3315 
3316    DO j = jts, MIN(jte,jde-1)
3317 
3318    DO k=kts,ktf
3319    DO i = i_start, i_end
3320    
3321      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)) &
3322        *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3323            - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3324        *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3325 
3326    ENDDO
3327    ENDDO
3328 
3329    IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3330 
3331      DO k=kts,ktf
3332    
3333        ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j))   &
3334          *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3335              - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3336          *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3337 
3338      ENDDO
3339 
3340    ENDIF
3341 
3342    IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3343 
3344      DO k=kts,ktf
3345    
3346        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)) &
3347          *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3348              - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3349          *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3350 
3351      ENDDO
3352 
3353    ENDIF
3354 
3355    ENDDO
3356 
3357 !  coriolis term for v-momentum equation
3358 
3359 !  Notes on map scale factors
3360 !  ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3361 !  Define f=2 omega sin(lat), e=2 omega cos(lat)
3362 !   => terms are: -f mu u / mx
3363 !  ru = mu u / my ; rw = mu w / my
3364 !   => terms are: -f ru *my / mx + ?
3365 
3366    j_start = jts
3367    j_end   = jte
3368 
3369    IF ( config_flags%open_ys .or. specified .or. &
3370         config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3371    IF ( config_flags%open_ye .or. specified .or. &
3372         config_flags%nested .or. config_flags%polar) j_end   = MIN(jde-1,jte)
3373 
3374    IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3375 
3376      DO k=kts,ktf
3377      DO i=its,MIN(ide-1,ite)
3378    
3379         rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts))    &
3380          *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts))   &
3381              + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts))   &
3382              *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts)) 
3383 
3384      ENDDO
3385      ENDDO
3386 
3387    ENDIF
3388 
3389    DO j=j_start, j_end
3390    DO k=kts,ktf
3391    DO i=its,MIN(ide-1,ite)
3392    
3393       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))    &
3394        *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3395            + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3396            *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j)) 
3397 
3398    ENDDO
3399    ENDDO
3400    ENDDO
3401 
3402 
3403    IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3404 
3405      DO k=kts,ktf
3406      DO i=its,MIN(ide-1,ite)
3407    
3408         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))        &
3409          *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1))   &
3410              + (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))   &
3411              *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1)) 
3412 
3413      ENDDO
3414      ENDDO
3415 
3416    ENDIF
3417 
3418 ! coriolis term for w-mometum 
3419 
3420 ! Notes on map scale factors
3421 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3422 ! Define e=2 omega cos(lat)
3423 !  => terms are: e mu u / my + ???
3424 ! ru = mu u / my ; ru = mu v / mx
3425 !  => terms are: e ru + ???
3426 
3427    DO j=jts,MIN(jte, jde-1)
3428    DO k=kts+1,ktf
3429    DO i=its,MIN(ite, ide-1)
3430 
3431        rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)*           &
3432           (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3433           +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j)))           &
3434           -(msftx(i,j)/msfty(i,j))*                      &
3435            sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3436           +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3437 
3438    ENDDO
3439    ENDDO
3440    ENDDO
3441 
3442 END SUBROUTINE coriolis
3443 
3444 !------------------------------------------------------------------------------
3445 
3446 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3447                                    config_flags,                                &
3448                                    u_base, v_base, z_base,                      &
3449                                    muu, muv, phb, ph,                           &
3450                                    msftx, msfty, msfux, msfuy, msfvx, msfvy,    &
3451                                    f, e, sina, cosa, fzm, fzp,                  &
3452                                    ids, ide, jds, jde, kds, kde,                &
3453                                    ims, ime, jms, jme, kms, kme,                &
3454                                    its, ite, jts, jte, kts, kte                )
3455 
3456    IMPLICIT NONE
3457    
3458    ! Input data
3459    
3460    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
3461 
3462    INTEGER ,                 INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3463                                               ims, ime, jms, jme, kms, kme, &
3464                                               its, ite, jts, jte, kts, kte
3465 
3466    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3467                                                                 rv_tend, &
3468                                                                 rw_tend
3469    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: ru_in, &
3470                                                                       rv_in, &
3471                                                                       rw,    &
3472                                                                       ph,    &
3473                                                                       phb
3474 
3475 
3476    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,      &
3477                                                                 msfuy,      &
3478                                                                 msfvx,      &
3479                                                                 msfvy,      &
3480                                                                 msftx,      &
3481                                                                 msfty
3482 
3483    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: f,    &
3484                                                                     e,    &
3485                                                                     sina, &
3486                                                                     cosa
3487 
3488    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: muu, &
3489                                                                     muv
3490                                                                     
3491 
3492    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm, &
3493                                                                   fzp
3494 
3495    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: u_base,  &
3496                                                                   v_base,  &
3497                                                                   z_base
3498    
3499    ! Local storage
3500 
3501    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3502                                                       rv
3503 
3504    REAL  :: z_at_u, z_at_v, wkp1, wk, wkm1
3505 
3506    ! Local indices.
3507    
3508    INTEGER :: i, j , k, ktf
3509    INTEGER :: i_start, i_end, j_start, j_end
3510    
3511    LOGICAL :: specified
3512 
3513 !<DESCRIPTION>
3514 !
3515 !  perturbation_coriolis calculates the large timestep tendency terms in the 
3516 !  u, v, and w momentum equations arise from the coriolis force.  This version
3517 !  subtracts off the horizontal velocities from the initial sounding when
3518 !  computing the forcing terms, hence "perturbation" coriolis.
3519 !
3520 !</DESCRIPTION>
3521 
3522    specified = .false.
3523    if(config_flags%specified .or. config_flags%nested) specified = .true.
3524 
3525    ktf=MIN(kte,kde-1)
3526 
3527 ! coriolis for u-momentum equation
3528 
3529    i_start = its
3530    i_end   = ite
3531    IF ( config_flags%open_xs .or. specified .or. &
3532         config_flags%nested) i_start = MAX(ids+1,its)
3533    IF ( config_flags%open_xe .or. specified .or. &
3534         config_flags%nested) i_end   = MIN(ide-1,ite)
3535       IF ( config_flags%periodic_x ) i_start = its
3536       IF ( config_flags%periodic_x ) i_end = ite
3537 
3538 !  compute perturbation mu*v for use in u momentum equation
3539 
3540    DO j = jts, MIN(jte,jde-1)+1
3541    DO k=kts+1,ktf-1
3542    DO i = i_start-1, i_end
3543      z_at_v = 0.25*( phb(i,k,j  )+phb(i,k+1,j  )  &
3544                     +phb(i,k,j-1)+phb(i,k+1,j-1)  &
3545                     +ph(i,k,j  )+ph(i,k+1,j  )    &
3546                     +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3547      wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3548      wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3549      wk   = 1.-wkp1-wkm1
3550      rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*(            &
3551                                   wkm1*v_base(k-1)    &
3552                                  +wk  *v_base(k  )    &
3553                                  +wkp1*v_base(k+1)   )
3554    ENDDO
3555    ENDDO
3556    ENDDO
3557 
3558 
3559 !  pick up top and bottom v 
3560 
3561    DO j = jts, MIN(jte,jde-1)+1
3562    DO i = i_start-1, i_end
3563 
3564      k = kts
3565      z_at_v = 0.25*( phb(i,k,j  )+phb(i,k+1,j  )  &
3566                     +phb(i,k,j-1)+phb(i,k+1,j-1)  &
3567                     +ph(i,k,j  )+ph(i,k+1,j  )    &
3568                     +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3569      wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3570      wk   = 1.-wkp1
3571      rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*(            &
3572                                  +wk  *v_base(k  )    &
3573                                  +wkp1*v_base(k+1)   )
3574 
3575      k = ktf
3576      z_at_v = 0.25*( phb(i,k,j  )+phb(i,k+1,j  )  &
3577                     +phb(i,k,j-1)+phb(i,k+1,j-1)  &
3578                     +ph(i,k,j  )+ph(i,k+1,j  )    &
3579                     +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3580      wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3581      wk   = 1.-wkm1
3582      rv(i,k,j) = rv_in(i,k,j) - muv(i,j)*(            &
3583                                   wkm1*v_base(k-1)    &
3584                                  +wk  *v_base(k  )   )
3585 
3586    ENDDO
3587    ENDDO
3588 
3589 !  compute coriolis forcing for u
3590 
3591 !  Map scale factors: see comments above for Coriolis
3592 
3593    DO j = jts, MIN(jte,jde-1)
3594 
3595    DO k=kts,ktf
3596      DO i = i_start, i_end
3597        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)) &
3598          *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3599              - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3600          *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3601      ENDDO
3602    ENDDO
3603 
3604    IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3605 
3606      DO k=kts,ktf
3607    
3608        ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j))   &
3609          *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3610              - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3611          *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3612 
3613      ENDDO
3614 
3615    ENDIF
3616 
3617    IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3618 
3619      DO k=kts,ktf
3620    
3621        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)) &
3622          *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3623              - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3624          *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3625 
3626      ENDDO
3627 
3628    ENDIF
3629 
3630    ENDDO
3631 
3632 !  coriolis term for v-momentum equation
3633 !  Map scale factors: see comments above for Coriolis
3634 
3635    j_start = jts
3636    j_end   = jte
3637 
3638    IF ( config_flags%open_ys .or. specified .or. &
3639         config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3640    IF ( config_flags%open_ye .or. specified .or. &
3641         config_flags%nested .or. config_flags%polar) j_end   = MIN(jde-1,jte)
3642 
3643 !  compute perturbation mu*u for use in v momentum equation
3644 
3645    DO j = j_start-1,j_end
3646    DO k=kts+1,ktf-1
3647    DO i = its, MIN(ite,ide-1)+1
3648      z_at_u = 0.25*( phb(i  ,k,j)+phb(i  ,k+1,j)  &
3649                     +phb(i-1,k,j)+phb(i-1,k+1,j)  &
3650                     +ph(i  ,k,j)+ph(i  ,k+1,j)    &
3651                     +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3652      wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3653      wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3654      wk   = 1.-wkp1-wkm1
3655      ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*(            &
3656                                   wkm1*u_base(k-1)    &
3657                                  +wk  *u_base(k  )    &
3658                                  +wkp1*u_base(k+1)   )
3659    ENDDO
3660    ENDDO
3661    ENDDO
3662 
3663 !  pick up top and bottom u
3664 
3665    DO j = j_start-1,j_end
3666    DO i = its, MIN(ite,ide-1)+1
3667 
3668      k = kts
3669      z_at_u = 0.25*( phb(i  ,k,j)+phb(i  ,k+1,j)  &
3670                     +phb(i-1,k,j)+phb(i-1,k+1,j)  &
3671                     +ph(i  ,k,j)+ph(i  ,k+1,j)    &
3672                     +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3673      wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
3674      wk   = 1.-wkp1
3675      ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*(            &
3676                                  +wk  *u_base(k  )    &
3677                                  +wkp1*u_base(k+1)   )
3678 
3679 
3680      k = ktf
3681      z_at_u = 0.25*( phb(i  ,k,j)+phb(i  ,k+1,j)  &
3682                     +phb(i-1,k,j)+phb(i-1,k+1,j)  &
3683                     +ph(i  ,k,j)+ph(i  ,k+1,j)    &
3684                     +ph(i-1,k,j)+ph(i-1,k+1,j))/g
3685      wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
3686      wk   = 1.-wkm1
3687      ru(i,k,j) = ru_in(i,k,j) - muu(i,j)*(            &
3688                                   wkm1*u_base(k-1)    &
3689                                  +wk  *u_base(k  )   )
3690 
3691    ENDDO
3692    ENDDO
3693 
3694 !  compute coriolis forcing for v momentum equation
3695 !  Map scale factors: see comments above for Coriolis
3696 
3697    IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3698 
3699      DO k=kts,ktf
3700      DO i=its,MIN(ide-1,ite)
3701    
3702         rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts))    &
3703          *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts))   &
3704              + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts))   &
3705              *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts)) 
3706 
3707      ENDDO
3708      ENDDO
3709 
3710    ENDIF
3711 
3712    DO j=j_start, j_end
3713    DO k=kts,ktf
3714    DO i=its,MIN(ide-1,ite)
3715    
3716       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))    &
3717        *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3718            + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3719            *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j)) 
3720 
3721    ENDDO
3722    ENDDO
3723    ENDDO
3724 
3725 
3726    IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3727 
3728      DO k=kts,ktf
3729      DO i=its,MIN(ide-1,ite)
3730    
3731         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))        &
3732          *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1))   &
3733              + (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))   &
3734              *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1)) 
3735 
3736      ENDDO
3737      ENDDO
3738 
3739    ENDIF
3740 
3741 ! coriolis term for w-mometum 
3742 !  Map scale factors: see comments above for Coriolis
3743 
3744    DO j=jts,MIN(jte, jde-1)
3745    DO k=kts+1,ktf
3746    DO i=its,MIN(ite, ide-1)
3747 
3748        rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)*           &
3749           (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3750           +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j)))           &
3751           -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3752           +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3753 
3754    ENDDO
3755    ENDDO
3756    ENDDO
3757 
3758 END SUBROUTINE perturbation_coriolis
3759 
3760 !------------------------------------------------------------------------------
3761 
3762 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
3763                         config_flags,                                       &
3764                         msfux, msfuy, msfvx, msfvy, msftx, msfty,       &
3765                         xlat, fzm, fzp, rdx, rdy,                       &
3766                         ids, ide, jds, jde, kds, kde,                   &
3767                         ims, ime, jms, jme, kms, kme,                   &
3768                         its, ite, jts, jte, kts, kte                   )
3769 
3770 
3771    IMPLICIT NONE
3772    
3773    ! Input data
3774 
3775    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
3776 
3777    INTEGER ,                  INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
3778                                                ims, ime, jms, jme, kms, kme, &
3779                                                its, ite, jts, jte, kts, kte
3780    
3781    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                     &
3782                                                INTENT(INOUT) :: ru_tend, &
3783                                                                 rv_tend, &
3784                                                                 rw_tend
3785 
3786    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                     &
3787                                                INTENT(IN   ) :: ru,      &
3788                                                                 rv,      &
3789                                                                 rw,      &
3790                                                                 u,       &
3791                                                                 v,       &
3792                                                                 w
3793 
3794    REAL , DIMENSION( ims:ime , jms:jme ) ,         INTENT(IN   ) :: msfux,    &
3795                                                                 msfuy,    &
3796                                                                 msfvx,    &
3797                                                                 msfvy,    &
3798                                                                 msftx,    &
3799                                                                 msfty,    &
3800                                                                 xlat
3801 
3802    REAL , DIMENSION( kms:kme ) ,                 INTENT(IN   ) :: fzm,     &
3803                                                                 fzp
3804 
3805    REAL ,                                      INTENT(IN   ) :: rdx,     &
3806                                                                 rdy
3807    
3808    ! Local data
3809    
3810 !   INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
3811    INTEGER :: i, j, k, itf, jtf, ktf
3812    INTEGER :: i_start, i_end, j_start, j_end
3813 !   INTEGER :: irmin, irmax, jrmin, jrmax
3814 
3815    REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
3816 
3817    LOGICAL :: specified
3818 
3819 !<DESCRIPTION>
3820 !
3821 !  curvature calculates the large timestep tendency terms in the 
3822 !  u, v, and w momentum equations arise from the curvature terms.  
3823 !
3824 !</DESCRIPTION>
3825 
3826    specified = .false.
3827    if(config_flags%specified .or. config_flags%nested) specified = .true.
3828 
3829       itf=MIN(ite,ide-1)
3830       jtf=MIN(jte,jde-1)
3831       ktf=MIN(kte,kde-1)
3832 
3833 !   irmin = ims
3834 !   irmax = ime
3835 !   jrmin = jms
3836 !   jrmax = jme
3837 !   IF ( config_flags%open_xs ) irmin = ids
3838 !   IF ( config_flags%open_xe ) irmax = ide-1
3839 !   IF ( config_flags%open_ys ) jrmin = jds
3840 !   IF ( config_flags%open_ye ) jrmax = jde-1
3841    
3842 ! Define v cross grad m at scalar points - vxgm(i,j)
3843 
3844    i_start = its-1
3845    i_end   = ite
3846    j_start = jts-1
3847    j_end   = jte
3848 
3849    IF ( ( config_flags%open_xs .or. specified .or. &
3850         config_flags%nested) .and. (its == ids) ) i_start = its
3851    IF ( ( config_flags%open_xe .or. specified .or. &
3852         config_flags%nested) .and. (ite == ide) ) i_end   = ite-1
3853    IF ( ( config_flags%open_ys .or. specified .or. &
3854         config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
3855    IF ( ( config_flags%open_ye .or. specified .or. &
3856         config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end   = jte-1
3857       IF ( config_flags%periodic_x ) i_start = its-1
3858       IF ( config_flags%periodic_x ) i_end = ite
3859 
3860    DO j=j_start, j_end
3861    DO k=kts,ktf
3862    DO i=i_start, i_end
3863 !     Map scale factor notes:
3864 !     msf...y is constant everywhere for cylindrical map projection
3865 !     msf...x varies with y only
3866 !     But we know that this is not = 0 for cylindrical,
3867 !     therefore use msfvX in 1st line
3868 !     which => by symmetry use msfuY in 2nd line - ???  
3869       vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
3870                   0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
3871    ENDDO
3872    ENDDO
3873    ENDDO
3874 
3875 !  Pick up the boundary rows for open (radiation) lateral b.c.
3876 !  Rather crude at present, we are assuming there is no
3877 !    variation in this term at the boundary.
3878 
3879    IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3880         config_flags%nested) .and. (its == ids) ) THEN
3881 
3882      DO j = jts, jte-1
3883      DO k = kts, ktf
3884        vxgm(its-1,k,j) =  vxgm(its,k,j)
3885      ENDDO
3886      ENDDO
3887 
3888    ENDIF
3889 
3890    IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
3891         config_flags%nested) .and. (ite == ide) ) THEN
3892 
3893      DO j = jts, jte-1
3894      DO k = kts, ktf
3895        vxgm(ite,k,j) =  vxgm(ite-1,k,j)
3896      ENDDO
3897      ENDDO
3898 
3899    ENDIF
3900 
3901 !  Polar boundary condition:
3902 !  The following change is needed in case one tries using the vxgm route with
3903 !  polar B.C.'s in the future, but not needed if 'tan' used
3904    IF ( ( config_flags%open_ys .or. specified .or. &
3905         config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
3906 
3907      DO k = kts, ktf
3908      DO i = its-1, ite
3909        vxgm(i,k,jts-1) =  vxgm(i,k,jts)
3910      ENDDO
3911      ENDDO
3912 
3913    ENDIF
3914 
3915 !  Polar boundary condition:
3916 !  The following change is needed in case one tries using the vxgm route with
3917 !  polar B.C.'s in the future, but not needed if 'tan' used
3918    IF ( ( config_flags%open_ye .or. specified .or. &
3919         config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
3920 
3921      DO k = kts, ktf
3922      DO i = its-1, ite
3923        vxgm(i,k,jte) =  vxgm(i,k,jte-1)
3924      ENDDO
3925      ENDDO
3926 
3927    ENDIF
3928 
3929 !  curvature term for u momentum eqn.
3930 
3931 !  Map scale factor notes:
3932 !  ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
3933 !                                               - mu u w /(a my)
3934 !  ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
3935 !   => terms are:
3936 !  (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
3937 !  ru v tan(lat) / a - u rw / a
3938 !  xlat defined with end points half grid space from pole,
3939 !  hence are on u latitude points
3940 
3941    i_start = its
3942    IF ( config_flags%open_xs .or. specified .or. &
3943         config_flags%nested) i_start = MAX ( ids+1 , its )
3944    IF ( config_flags%open_xe .or. specified .or. &
3945         config_flags%nested) i_end   = MIN ( ide-1 , ite )
3946       IF ( config_flags%periodic_x ) i_start = its
3947       IF ( config_flags%periodic_x ) i_end = ite
3948 
3949 !  Polar boundary condition
3950    IF (config_flags%polar) THEN
3951 
3952       DO j=jts,MIN(jde-1,jte)
3953       DO k=kts,ktf
3954       DO i=i_start,i_end
3955 
3956             ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius*                 ( &
3957                         (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
3958                                     rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
3959                         - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
3960       ENDDO
3961       ENDDO
3962       ENDDO
3963 
3964    ELSE  ! normal code
3965 
3966 
3967    DO j=jts,MIN(jde-1,jte)
3968    DO k=kts,ktf
3969    DO i=i_start,i_end
3970 
3971       ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
3972               *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3973                - u(i,k,j)*reradius &
3974               *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3975 
3976    ENDDO
3977    ENDDO
3978    ENDDO
3979 
3980    END IF
3981 
3982 !  curvature term for v momentum eqn.
3983 
3984 !  Map scale factor notes
3985 !  ADT eqn 45, RHS terms 4 and 5, in cylindrical: mu u*u tan(lat)/(a mx)
3986 !                                               - mu v w /(a mx)
3987 !  ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
3988 !  terms are:
3989 !  (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a 
3990 !  = [my/(mx*a)]*[u ru tan(lat) - v rw]
3991 !  (1/a)*[(my/mx)*u ru tan(lat) - w rv]
3992 !  xlat defined with end points half grid space from pole, hence are on
3993 !  u latitude points => av here
3994 !
3995 !  in original wrf, there was a sign error for the rw contribution
3996 
3997    j_start = jts
3998    IF ( config_flags%open_ys .or. specified .or. &
3999         config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4000    IF ( config_flags%open_ye .or. specified .or. &
4001         config_flags%nested .or. config_flags%polar) j_end   = MIN ( jde-1 , jte )
4002 
4003    IF (config_flags%polar) THEN
4004 
4005       DO j=j_start,j_end
4006       DO k=kts,ktf
4007       DO i=its,MIN(ite,ide-1)
4008             rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius*   (  &
4009                         0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))*     &
4010                         tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)*                &
4011                         0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1))  &
4012                        - v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+              &
4013                                                       rw(i,k+1,j)+rw(i,k,j))    )
4014       ENDDO
4015       ENDDO
4016       ENDDO
4017 
4018    ELSE  ! normal code
4019 
4020    DO j=j_start,j_end
4021    DO k=kts,ktf
4022    DO i=its,MIN(ite,ide-1)
4023 
4024       rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4025               *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4026                     - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius       &
4027               *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4028 
4029    ENDDO
4030    ENDDO
4031    ENDDO
4032 
4033    END IF
4034 
4035 !  curvature term for vertical momentum eqn.
4036 
4037 !  Notes on map scale factors:
4038 !  ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4039 !  ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4040 !  terms are: u ru / a + (mx/my)v rv / a
4041 
4042    DO j=jts,MIN(jte,jde-1)
4043    DO k=MAX(2,kts),ktf
4044    DO i=its,MIN(ite,ide-1)
4045 
4046       rw_tend(i,k,j)=rw_tend(i,k,j) + reradius*                              &
4047     (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))) &
4048     *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)))     &
4049     +(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))) &
4050     *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))))
4051 
4052    ENDDO
4053    ENDDO
4054    ENDDO
4055 
4056 END SUBROUTINE curvature
4057 
4058 !------------------------------------------------------------------------------
4059 
4060 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4061                       fzm, fzp,                          &
4062                       ids, ide, jds, jde, kds, kde,      &
4063                       ims, ime, jms, jme, kms, kme,      &
4064                       its, ite, jts, jte, kts, kte      )
4065 
4066    IMPLICIT NONE
4067 
4068    ! Input data
4069 
4070    TYPE(grid_config_rec_type) ,           INTENT(IN   ) :: config_flags   
4071 
4072    INTEGER ,                                   INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4073                                                                 ims, ime, jms, jme, kms, kme, &
4074                                                                 its, ite, jts, jte, kts, kte
4075 
4076    CHARACTER(LEN=1) ,                          INTENT(IN   ) :: name
4077 
4078    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: rfield
4079 
4080    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN   ) :: rr
4081    
4082    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT  ) :: field
4083    
4084    REAL , DIMENSION( kms:kme ) , INTENT(IN   ) :: fzm, fzp
4085    
4086    ! Local data
4087    
4088    INTEGER :: i, j, k, itf, jtf, ktf
4089    
4090 !<DESCRIPTION>
4091 !
4092 !  decouple decouples a variable from the column dry-air mass.
4093 !
4094 !</DESCRIPTION>
4095 
4096    ktf=MIN(kte,kde-1)
4097    
4098    IF (name .EQ. 'u')THEN
4099       itf=ite
4100       jtf=MIN(jte,jde-1)
4101 
4102       DO j=jts,jtf
4103       DO k=kts,ktf
4104       DO i=its,itf
4105          field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4106       ENDDO
4107       ENDDO
4108       ENDDO
4109 
4110    ELSE IF (name .EQ. 'v')THEN
4111       itf=MIN(ite,ide-1)
4112       jtf=jte
4113 
4114       DO j=jts,jtf
4115       DO k=kts,ktf
4116         DO i=its,itf
4117              field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4118         ENDDO
4119       ENDDO
4120       ENDDO
4121 
4122    ELSE IF (name .EQ. 'w')THEN
4123       itf=MIN(ite,ide-1)
4124       jtf=MIN(jte,jde-1)
4125       DO j=jts,jtf
4126       DO k=kts+1,ktf
4127       DO i=its,itf
4128          field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4129       ENDDO
4130       ENDDO
4131       ENDDO
4132 
4133       DO j=jts,jtf
4134       DO i=its,itf
4135         field(i,kte,j) = 0.
4136       ENDDO
4137       ENDDO
4138 
4139    ELSE 
4140       itf=MIN(ite,ide-1)
4141       jtf=MIN(jte,jde-1)
4142    ! For theta we will decouple tb and tp and add them to give t afterwards
4143       DO j=jts,jtf
4144       DO k=kts,ktf
4145       DO i=its,itf
4146          field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4147       ENDDO
4148       ENDDO
4149       ENDDO
4150    
4151    ENDIF
4152 
4153 END SUBROUTINE decouple
4154 
4155 !-------------------------------------------------------------------------------
4156 
4157 
4158 SUBROUTINE zero_tend ( tendency,                     &
4159                        ids, ide, jds, jde, kds, kde, &
4160                        ims, ime, jms, jme, kms, kme, &
4161                        its, ite, jts, jte, kts, kte )
4162 
4163 
4164    IMPLICIT NONE
4165    
4166    ! Input data
4167    
4168    INTEGER ,                                   INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4169                                                                 ims, ime, jms, jme, kms, kme, &
4170                                                                 its, ite, jts, jte, kts, kte
4171 
4172    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4173 
4174    ! Local data
4175    
4176    INTEGER :: i, j, k, itf, jtf, ktf
4177 
4178 !<DESCRIPTION>
4179 !
4180 !  zero_tend sets the input tendency array to zero.
4181 !
4182 !</DESCRIPTION>
4183 
4184       DO j = jts, jte
4185       DO k = kts, kte
4186       DO i = its, ite
4187         tendency(i,k,j) = 0.
4188       ENDDO
4189       ENDDO
4190       ENDDO
4191 
4192       END SUBROUTINE zero_tend
4193 
4194 !-------------------------------------------------------------------------------
4195 ! Sets the an array on the polar v point(s) to zero
4196 SUBROUTINE zero_pole ( field,                        &
4197                        ids, ide, jds, jde, kds, kde, &
4198                        ims, ime, jms, jme, kms, kme, &
4199                        its, ite, jts, jte, kts, kte )
4200 
4201 
4202   IMPLICIT NONE
4203 
4204   ! Input data
4205    
4206   INTEGER , INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4207                              ims, ime, jms, jme, kms, kme, &
4208                              its, ite, jts, jte, kts, kte
4209 
4210   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4211 
4212   ! Local data
4213 
4214   INTEGER :: i, k
4215 
4216   IF (jts == jds) THEN
4217      DO k = kts, kte
4218      DO i = its, ite
4219         field(i,k,jts) = 0.
4220      END DO
4221      END DO
4222   END IF
4223   IF (jte == jde) THEN
4224      DO k = kts, kte
4225      DO i = its, ite
4226         field(i,k,jte) = 0.
4227      END DO
4228      END DO
4229   END IF
4230 
4231 END SUBROUTINE zero_pole
4232 
4233 !-------------------------------------------------------------------------------
4234 ! Sets the an array on the polar v point(s)
4235 SUBROUTINE pole_point_bc ( field,                        &
4236                        ids, ide, jds, jde, kds, kde, &
4237                        ims, ime, jms, jme, kms, kme, &
4238                        its, ite, jts, jte, kts, kte )
4239 
4240 
4241   IMPLICIT NONE
4242 
4243   ! Input data
4244    
4245   INTEGER , INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4246                              ims, ime, jms, jme, kms, kme, &
4247                              its, ite, jts, jte, kts, kte
4248 
4249   REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4250 
4251   ! Local data
4252 
4253   INTEGER :: i, k
4254 
4255   IF (jts == jds) THEN
4256      DO k = kts, kte
4257      DO i = its, ite
4258 !        field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4259         field(i,k,jts) = field(i,k,jts+1)
4260      END DO
4261      END DO
4262   END IF
4263   IF (jte == jde) THEN
4264      DO k = kts, kte
4265      DO i = its, ite
4266 !        field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4267         field(i,k,jte) = field(i,k,jte-1)
4268      END DO
4269      END DO
4270   END IF
4271 
4272 END SUBROUTINE pole_point_bc
4273 
4274 !======================================================================
4275 !   physics prep routines
4276 !======================================================================
4277 
4278    SUBROUTINE phy_prep ( config_flags,                                &  ! input
4279                          mu, muu, muv, u, v, p, pb, alt, ph,          &  ! input
4280                          phb, t, tsk, moist, n_moist,                 &  ! input
4281                          mu_3d, rho, th_phy, p_phy , pi_phy ,         &  ! output
4282                          u_phy, v_phy, p8w, t_phy, t8w,               &  ! output
4283                          z, z_at_w, dz8w,                             &  ! output
4284                          fzm, fzp,                                    &  ! params
4285                          RTHRATEN,                                    &
4286                          RTHBLTEN, RUBLTEN, RVBLTEN,                  &
4287                          RQVBLTEN, RQCBLTEN, RQIBLTEN,                &
4288                          RTHCUTEN, RQVCUTEN, RQCCUTEN,                &
4289                          RQRCUTEN, RQICUTEN, RQSCUTEN,                &
4290                          RTHFTEN,  RQVFTEN,                           &
4291                          RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN,            &
4292                          RQVNDGDTEN, RMUNDGDTEN,                      &
4293                          ids, ide, jds, jde, kds, kde,                &
4294                          ims, ime, jms, jme, kms, kme,                &
4295                          its, ite, jts, jte, kts, kte                )
4296 !----------------------------------------------------------------------
4297    IMPLICIT NONE
4298 !----------------------------------------------------------------------
4299 
4300    TYPE(grid_config_rec_type) ,     INTENT(IN   ) :: config_flags
4301 
4302    INTEGER ,        INTENT(IN   ) ::   ids, ide, jds, jde, kds, kde, &
4303                                        ims, ime, jms, jme, kms, kme, &
4304                                        its, ite, jts, jte, kts, kte
4305    INTEGER ,          INTENT(IN   ) :: n_moist
4306 
4307    REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4308 
4309 
4310    REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN   )   ::     TSK, mu, muu, muv
4311 
4312    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
4313           INTENT(  OUT)                                  ::   u_phy, &
4314                                                               v_phy, &
4315                                                              pi_phy, &
4316                                                               p_phy, &
4317                                                                 p8w, &
4318                                                               t_phy, &
4319                                                              th_phy, &
4320                                                                 t8w, &
4321                                                               mu_3d, &
4322                                                                 rho, &
4323                                                                   z, &
4324                                                                dz8w, &
4325                                                               z_at_w 
4326 
4327    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) ,                 &
4328           INTENT(IN   )                                  ::      pb, &
4329                                                                   p, &
4330                                                                   u, &
4331                                                                   v, &
4332                                                                 alt, &
4333                                                                  ph, &
4334                                                                 phb, &
4335                                                                   t
4336 
4337 
4338    REAL , DIMENSION( kms:kme ) ,           INTENT(IN   ) ::     fzm,   &
4339                                                                 fzp
4340 
4341    REAL,  DIMENSION( ims:ime , kms:kme, jms:jme ),                   &
4342           INTENT(INOUT)   ::                               RTHRATEN  
4343 
4344    REAL,  DIMENSION( ims:ime , kms:kme, jms:jme ),                   &
4345           INTENT(INOUT)   ::                               RTHCUTEN, &
4346                                                            RQVCUTEN, &
4347                                                            RQCCUTEN, &
4348                                                            RQRCUTEN, &
4349                                                            RQICUTEN, &
4350                                                            RQSCUTEN
4351 
4352    REAL,  DIMENSION( ims:ime, kms:kme, jms:jme )                   , &
4353           INTENT(INOUT)   ::                                RUBLTEN, &
4354                                                             RVBLTEN, &
4355                                                            RTHBLTEN, &
4356                                                            RQVBLTEN, &
4357                                                            RQCBLTEN, &
4358                                                            RQIBLTEN
4359 
4360    REAL,  DIMENSION( ims:ime, kms:kme, jms:jme )                   , &
4361           INTENT(INOUT)   ::                                RTHFTEN, &
4362                                                             RQVFTEN
4363 
4364    REAL,  DIMENSION( ims:ime, kms:kme, jms:jme )                   , &
4365           INTENT(INOUT)   ::                                RUNDGDTEN, &
4366                                                             RVNDGDTEN, &
4367                                                            RTHNDGDTEN, &
4368                                                            RQVNDGDTEN, &
4369                                                            RMUNDGDTEN
4370 
4371    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4372    INTEGER :: i, j, k
4373    REAL    :: w1, w2, z0, z1, z2
4374 
4375 !-----------------------------------------------------------------------
4376 
4377 !<DESCRIPTION>
4378 !
4379 !  phys_prep calculates a number of diagnostic quantities needed by
4380 !  the physics routines.  It also decouples the physics tendencies from
4381 !  the column dry-air mass (the physics routines expect to see/update the
4382 !  uncoupled tendencies).
4383 !
4384 !</DESCRIPTION>
4385 
4386 !  set up loop bounds for this grid's boundary conditions
4387 
4388     i_start = its
4389     i_end   = min( ite,ide-1 )
4390     j_start = jts
4391     j_end   = min( jte,jde-1 )
4392 
4393     k_start = kts
4394     k_end = min( kte, kde-1 )
4395 
4396 !  compute thermodynamics and velocities at pressure points
4397 
4398     do j = j_start,j_end
4399     do k = k_start, k_end
4400     do i = i_start, i_end
4401 
4402       th_phy(i,k,j) = t(i,k,j) + t0
4403       p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4404       pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4405       t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4406       rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4407       mu_3d(i,k,j) = mu(i,j)
4408       u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4409       v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4410 
4411     enddo
4412     enddo
4413     enddo
4414 
4415 !  compute z at w points
4416 
4417     do j = j_start,j_end
4418     do k = k_start, kte
4419     do i = i_start, i_end
4420       z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4421     enddo
4422     enddo
4423     enddo
4424 
4425     do j = j_start,j_end
4426     do k = k_start, kte-1
4427     do i = i_start, i_end
4428       dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4429     enddo
4430     enddo
4431     enddo
4432 
4433     do j = j_start,j_end
4434     do i = i_start, i_end
4435       dz8w(i,kte,j) = 0.
4436     enddo
4437     enddo
4438 
4439 !  compute z at p points (average of z at w points)
4440 
4441     do j = j_start,j_end
4442     do k = k_start, k_end
4443     do i = i_start, i_end
4444       z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4445     enddo
4446     enddo
4447     enddo
4448 
4449 !  interp t and p at w points
4450 
4451     do j = j_start,j_end
4452     do k = 2, k_end
4453     do i = i_start, i_end
4454       p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4455       t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4456     enddo
4457     enddo
4458     enddo
4459 
4460 !  extrapolate p and t to surface and top.
4461 !  we'll use an extrapolation in z for now
4462 
4463     do j = j_start,j_end
4464     do i = i_start, i_end
4465 
4466 ! bottom
4467 
4468       z0 = z_at_w(i,1,j)
4469       z1 = z(i,1,j)
4470       z2 = z(i,2,j)
4471       w1 = (z0 - z2)/(z1 - z2)
4472       w2 = 1. - w1
4473       p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4474       t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4475 
4476 ! top
4477 
4478       z0 = z_at_w(i,kte,j)
4479       z1 = z(i,k_end,j)
4480       z2 = z(i,k_end-1,j)
4481       w1 = (z0 - z2)/(z1 - z2)
4482       w2 = 1. - w1
4483 
4484 !      p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4485 !!!  bug fix      extrapolate ln(p) so p is positive definite
4486       p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4487       t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4488 
4489     enddo
4490     enddo
4491 
4492 ! decouple all physics tendencies
4493 
4494    IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
4495 
4496       DO J=j_start,j_end
4497       DO K=k_start,k_end
4498       DO I=i_start,i_end
4499          RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/mu(I,J)
4500       ENDDO
4501       ENDDO
4502       ENDDO
4503 
4504    ENDIF
4505 
4506    IF (config_flags%cu_physics .gt. 0) THEN
4507 
4508       DO J=j_start,j_end
4509       DO I=i_start,i_end
4510       DO K=k_start,k_end
4511          RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/mu(I,J)
4512       ENDDO
4513       ENDDO
4514       ENDDO
4515 
4516       IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4517          DO J=j_start,j_end
4518          DO I=i_start,i_end
4519          DO K=k_start,k_end
4520             RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/mu(I,J)
4521          ENDDO
4522          ENDDO
4523          ENDDO
4524       ENDIF
4525 
4526       IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
4527          DO J=j_start,j_end
4528          DO I=i_start,i_end
4529          DO K=k_start,k_end
4530             RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/mu(I,J)
4531          ENDDO
4532          ENDDO
4533          ENDDO
4534       ENDIF
4535 
4536       IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
4537          DO J=j_start,j_end
4538          DO I=i_start,i_end
4539          DO K=k_start,k_end
4540             RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/mu(I,J)
4541          ENDDO
4542          ENDDO
4543          ENDDO
4544       ENDIF
4545 
4546       IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
4547          DO J=j_start,j_end
4548          DO I=i_start,i_end
4549          DO K=k_start,k_end
4550             RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/mu(I,J)
4551          ENDDO
4552          ENDDO
4553          ENDDO
4554       ENDIF
4555 
4556       IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
4557          DO J=j_start,j_end
4558          DO I=i_start,i_end
4559          DO K=k_start,k_end
4560             RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/mu(I,J)
4561          ENDDO
4562          ENDDO
4563          ENDDO
4564       ENDIF
4565 
4566    ENDIF
4567 
4568    IF (config_flags%bl_pbl_physics .gt. 0) THEN
4569 
4570       DO J=j_start,j_end
4571       DO K=k_start,k_end
4572       DO I=i_start,i_end
4573          RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/mu(I,J)
4574          RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/mu(I,J)
4575          RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/mu(I,J)
4576       ENDDO
4577       ENDDO
4578       ENDDO
4579 
4580       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4581          DO J=j_start,j_end
4582          DO K=k_start,k_end
4583          DO I=i_start,i_end
4584             RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/mu(I,J)
4585          ENDDO
4586          ENDDO
4587          ENDDO
4588       ENDIF
4589 
4590       IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
4591          DO J=j_start,j_end
4592          DO K=k_start,k_end
4593          DO I=i_start,i_end
4594            RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/mu(I,J)
4595          ENDDO
4596          ENDDO
4597          ENDDO
4598       ENDIF
4599 
4600       IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
4601          DO J=j_start,j_end
4602          DO K=k_start,k_end
4603          DO I=i_start,i_end
4604             RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/mu(I,J)
4605          ENDDO
4606          ENDDO
4607          ENDDO
4608       ENDIF
4609 
4610     ENDIF
4611 
4612 !  decouple advective forcing required by Grell-Devenyi scheme
4613 
4614    if ( config_flags%cu_physics == GDSCHEME ) then
4615 
4616       DO J=j_start,j_end
4617       DO I=i_start,i_end
4618          DO K=k_start,k_end
4619             RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/mu(I,J)
4620          ENDDO
4621       ENDDO
4622       ENDDO
4623 
4624       IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
4625          DO J=j_start,j_end
4626          DO I=i_start,i_end
4627             DO K=k_start,k_end
4628                RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/mu(I,J)
4629             ENDDO
4630          ENDDO
4631          ENDDO
4632       ENDIF
4633 
4634    END IF
4635 
4636 ! fdda
4637 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
4638 !   so only decouple those
4639 
4640    IF (config_flags%grid_fdda .gt. 0) THEN
4641 
4642       i_startu=MAX(its,ids+1)
4643       j_startv=MAX(jts,jds+1)
4644 
4645       DO J=j_start,j_end
4646       DO K=k_start,k_end
4647       DO I=i_startu,i_end
4648          RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/muu(I,J)
4649       ENDDO
4650       ENDDO
4651       ENDDO
4652       DO J=j_startv,j_end
4653       DO K=k_start,k_end
4654       DO I=i_start,i_end
4655          RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/muv(I,J)
4656       ENDDO
4657       ENDDO
4658       ENDDO
4659       DO J=j_start,j_end
4660       DO K=k_start,k_end
4661       DO I=i_start,i_end
4662          RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/mu(I,J)
4663 !        RMUNDGDTEN(I,J) - no coupling
4664       ENDDO
4665       ENDDO
4666       ENDDO
4667       IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
4668          DO J=j_start,j_end
4669          DO K=k_start,k_end
4670          DO I=i_start,i_end
4671             RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/mu(I,J)
4672          ENDDO
4673          ENDDO
4674          ENDDO
4675       ENDIF
4676 
4677     ENDIF
4678 
4679 END SUBROUTINE phy_prep
4680 
4681 !------------------------------------------------------------
4682 
4683    SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
4684                                      p, p8w, p0, pb, ph, phb,        &
4685                                      th_phy, pii, pf,                &
4686                                      z, z_at_w, dz8w,                &
4687                                      dt,h_diabatic,                  &
4688                                      config_flags,fzm, fzp,          &
4689                                      ids,ide, jds,jde, kds,kde,      &
4690                                      ims,ime, jms,jme, kms,kme,      &
4691                                      its,ite, jts,jte, kts,kte      )
4692 
4693    IMPLICIT NONE
4694 
4695 ! Here we construct full fields
4696 ! needed by the microphysics
4697 
4698    TYPE(grid_config_rec_type),    INTENT(IN   )    :: config_flags
4699 
4700    INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
4701    INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
4702    INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
4703 
4704    REAL, INTENT(IN   )  ::  dt
4705 
4706    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
4707          INTENT(IN   ) ::                           al,  &
4708                                                     alb, &
4709                                                     p,   &
4710                                                     pb,  &
4711                                                     ph,  &
4712                                                     phb
4713 
4714 
4715    REAL , DIMENSION( kms:kme ) ,           INTENT(IN   ) ::   fzm, &
4716                                                               fzp
4717 
4718    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),       &
4719          INTENT(  OUT) ::                         rho,  &
4720                                                th_phy,  &
4721                                                   pii,  &
4722                                                   pf,   &
4723                                                     z,  &
4724                                                z_at_w,  &
4725                                                  dz8w,  &
4726                                                   p8w
4727 ! pjj/cray
4728 !                                                 p8w,  &
4729 !                                          h_diabatic
4730 
4731    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),       &
4732          INTENT(INOUT) ::                         h_diabatic
4733 
4734    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
4735          INTENT(INOUT) ::                         t_new, &
4736                                                   t_old
4737 
4738    REAL, INTENT(IN   ) :: t0, p0
4739    REAL                :: z0,z1,z2,w1,w2
4740 
4741    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4742    INTEGER :: i, j, k
4743 
4744 !--------------------------------------------------------------------
4745 
4746 !<DESCRIPTION>
4747 !
4748 !  moist_phys_prep_em calculates a number of diagnostic quantities needed by
4749 !  the microphysics routines.
4750 !
4751 !</DESCRIPTION>
4752 
4753 !  set up loop bounds for this grid's boundary conditions
4754 
4755     i_start = its    
4756     i_end   = min( ite,ide-1 )
4757     j_start = jts    
4758     j_end   = min( jte,jde-1 )
4759 
4760     k_start = kts
4761     k_end = min( kte, kde-1 )
4762 
4763      DO j = j_start, j_end
4764      DO k = k_start, kte
4765      DO i = i_start, i_end
4766        z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
4767      ENDDO
4768      ENDDO
4769      ENDDO
4770 
4771     do j = j_start,j_end
4772     do k = k_start, kte-1
4773     do i = i_start, i_end
4774       dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4775     enddo
4776     enddo
4777     enddo
4778 
4779     do j = j_start,j_end
4780     do i = i_start, i_end
4781       dz8w(i,kte,j) = 0.
4782     enddo
4783     enddo
4784 
4785 
4786            !  compute full pii, rho, and z at the new time-level
4787            !  (needed for physics).
4788            !  convert perturbation theta to full theta (th_phy)
4789            !  use h_diabatic to temporarily save pre-microphysics full theta
4790 
4791      DO j = j_start, j_end
4792      DO k = k_start, k_end
4793      DO i = i_start, i_end
4794 
4795 #ifdef REVERT
4796        t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt
4797 #endif
4798        th_phy(i,k,j) = t_new(i,k,j) + t0
4799        h_diabatic(i,k,j) = th_phy(i,k,j)
4800        rho(i,k,j)  = 1./(al(i,k,j)+alb(i,k,j))
4801        pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
4802        z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4803        pf(i,k,j) = p(i,k,j)+pb(i,k,j)
4804 
4805      ENDDO
4806      ENDDO
4807      ENDDO
4808 
4809 !  interp t and p at w points
4810 
4811     do j = j_start,j_end
4812     do k = 2, k_end
4813     do i = i_start, i_end
4814       p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
4815     enddo
4816     enddo
4817     enddo
4818 
4819 !  extrapolate p and t to surface and top.
4820 !  we'll use an extrapolation in z for now
4821 
4822     do j = j_start,j_end
4823     do i = i_start, i_end
4824 
4825 ! bottom
4826 
4827       z0 = z_at_w(i,1,j)
4828       z1 = z(i,1,j)
4829       z2 = z(i,2,j)
4830       w1 = (z0 - z2)/(z1 - z2)
4831       w2 = 1. - w1
4832       p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
4833 
4834 ! top
4835 
4836       z0 = z_at_w(i,kte,j)
4837       z1 = z(i,k_end,j)
4838       z2 = z(i,k_end-1,j)
4839       w1 = (z0 - z2)/(z1 - z2)
4840       w2 = 1. - w1
4841 !      p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
4842       p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
4843 
4844     enddo
4845     enddo
4846 
4847    END SUBROUTINE moist_physics_prep_em
4848 
4849 !------------------------------------------------------------------------------
4850 
4851    SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut,     &
4852                                        th_phy, h_diabatic, dt,    &
4853                                        config_flags,              &
4854                                        ids,ide, jds,jde, kds,kde, &
4855                                        ims,ime, jms,jme, kms,kme, &
4856                                        its,ite, jts,jte, kts,kte )
4857 
4858    IMPLICIT NONE
4859 
4860 ! Here we construct full fields
4861 ! needed by the microphysics
4862 
4863    TYPE(grid_config_rec_type),    INTENT(IN   )    :: config_flags
4864 
4865    INTEGER,      INTENT(IN   )    :: ids,ide, jds,jde, kds,kde
4866    INTEGER,      INTENT(IN   )    :: ims,ime, jms,jme, kms,kme
4867    INTEGER,      INTENT(IN   )    :: its,ite, jts,jte, kts,kte
4868 
4869    REAL, DIMENSION( ims:ime , kms:kme, jms:jme ),        &
4870          INTENT(INOUT) ::                         t_new, &
4871                                                   t_old, &
4872                                                  th_phy, &
4873                                                   h_diabatic
4874 
4875    REAL, DIMENSION( ims:ime , jms:jme ),  INTENT(INOUT) ::  mut
4876 
4877 
4878    REAL, INTENT(IN   ) :: t0, dt
4879 
4880    INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
4881    INTEGER :: i, j, k
4882 
4883 !--------------------------------------------------------------------
4884 
4885 !<DESCRIPTION>
4886 !
4887 !  moist_phys_finish_em resets theta to its perturbation value and
4888 !  computes and stores the microphysics diabatic heating term.
4889 !
4890 !</DESCRIPTION>
4891 
4892 !  set up loop bounds for this grid's boundary conditions
4893 
4894 
4895     i_start = its    
4896     i_end   = min( ite,ide-1 )
4897     j_start = jts    
4898     j_end   = min( jte,jde-1 )
4899 
4900     k_start = kts
4901     k_end = min( kte, kde-1 )
4902 
4903 !  add microphysics theta diff to perturbation theta, set h_diabatic
4904 
4905      DO j = j_start, j_end
4906      DO k = k_start, k_end
4907      DO i = i_start, i_end
4908 
4909        t_new(i,k,j) = t_new(i,k,j) + (th_phy(i,k,j)-h_diabatic(i,k,j))
4910        h_diabatic(i,k,j) = (th_phy(i,k,j)-h_diabatic(i,k,j))/dt
4911 !       h_diabatic(i,k,j) = 0.
4912 
4913      ENDDO
4914      ENDDO
4915      ENDDO
4916 
4917    END SUBROUTINE moist_physics_finish_em
4918 
4919 !----------------------------------------------------------------
4920 
4921 
4922    SUBROUTINE init_module_big_step
4923    END SUBROUTINE init_module_big_step
4924 
4925 SUBROUTINE set_tend ( field, field_adv_tend, msf,       &
4926                       ids, ide, jds, jde, kds, kde,     &
4927                       ims, ime, jms, jme, kms, kme,     &
4928                       its, ite, jts, jte, kts, kte       )
4929 
4930    IMPLICIT NONE
4931 
4932    ! Input data
4933 
4934    INTEGER ,  INTENT(IN   ) :: ids, ide, jds, jde, kds, kde, &
4935                                ims, ime, jms, jme, kms, kme, &
4936                                its, ite, jts, jte, kts, kte
4937 
4938    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
4939 
4940    REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN)  :: field_adv_tend
4941 
4942    REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN)  :: msf
4943 
4944    ! Local data
4945 
4946    INTEGER :: i, j, k, itf, jtf, ktf
4947 
4948 !<DESCRIPTION>
4949 !
4950 !  set_tend copies the advective tendency array into the tendency array.
4951 !
4952 !</DESCRIPTION>
4953 
4954       jtf = MIN(jte,jde-1)
4955       ktf = MIN(kte,kde-1)
4956       itf = MIN(ite,ide-1)
4957       DO j = jts, jtf
4958       DO k = kts, ktf
4959       DO i = its, itf
4960          field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
4961       ENDDO
4962       ENDDO
4963       ENDDO
4964 
4965 END SUBROUTINE set_tend
4966 
4967 !------------------------------------------------------------------------------
4968 
4969     SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf,              &
4970                                  rw_tendf, t_tendf,               &
4971                                  u, v, w, t, t_init,              &
4972                                  mut, muu, muv, ph, phb,          &
4973                                  u_base, v_base, t_base, z_base,  &
4974                                  dampcoef, zdamp,                 &
4975                                  ids, ide, jds, jde, kds, kde,    &
4976                                  ims, ime, jms, jme, kms, kme,    &
4977                                  its, ite, jts, jte, kts, kte   )
4978 
4979 ! History:     Apr 2005  Modifications by George Bryan, NCAR:
4980 !                  - Generalized the code in a way that allows for
4981 !                    simulations with steep terrain.
4982 !
4983 !              Jul 2004  Modifications by George Bryan, NCAR:
4984 !                  - Modified the code to use u_base, v_base, and t_base
4985 !                    arrays for the background state.  Removed the hard-wired
4986 !                    base-state values.
4987 !                  - Modified the code to use dampcoef, zdamp, and damp_opt,
4988 !                    i.e., the upper-level damper variables in namelist.input.
4989 !                    Removed the hard-wired variables in the older version.
4990 !                    This damper is used when damp_opt = 2.
4991 !                  - Modified the code to account for the movement of the
4992 !                    model surfaces with time.  The code now obtains a base-
4993 !                    state value by interpolation using the "_base" arrays.
4994 
4995 !              Nov 2003  Bug fix by Jason Knievel, NCAR
4996 
4997 !              Aug 2003  Meridional dimension, some comments, and
4998 !                        changes in layout of the code added by
4999 !                        Jason Knievel, NCAR
5000 
5001 !              Jul 2003  Original code by Bill Skamarock, NCAR
5002 
5003 ! Purpose:     This routine applies Rayleigh damping to a layer at top
5004 !              of the model domain.
5005 
5006 !-----------------------------------------------------------------------
5007 ! Begin declarations.
5008 
5009     IMPLICIT NONE
5010 
5011     INTEGER, INTENT( IN )  &
5012     :: ids, ide, jds, jde, kds, kde,  &
5013        ims, ime, jms, jme, kms, kme,  &
5014        its, ite, jts, jte, kts, kte
5015 
5016     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT )  &
5017     :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5018 
5019     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN )  &
5020     :: u, v, w, t, t_init, ph, phb
5021 
5022     REAL, DIMENSION( ims:ime, jms:jme ),  INTENT( IN )  &
5023     :: mut, muu, muv
5024 
5025     REAL, DIMENSION( kms:kme ) ,  INTENT(IN   )  &
5026     :: u_base, v_base, t_base, z_base
5027 
5028     REAL, INTENT(IN   )   &
5029     :: dampcoef, zdamp
5030 
5031 ! Local variables.
5032 
5033     INTEGER  &
5034     :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5035 
5036     REAL  &
5037     :: pii, dcoef, z, ztop
5038 
5039     REAL :: wkp1, wk, wkm1
5040 
5041     REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5042 
5043 ! End declarations.
5044 !-----------------------------------------------------------------------
5045 
5046     pii = 2.0 * asin(1.0)
5047 
5048     ktf = MIN( kte,   kde-1 )
5049 
5050 !-----------------------------------------------------------------------
5051 ! Adjust u to base state.
5052 
5053     DO j = jts, MIN( jte, jde-1 )
5054     DO i = its, MIN( ite, ide   )
5055 
5056       ! Get height at top of model
5057       ztop = 0.5*( phb(i  ,kde,j)+phb(i-1,kde,j)   &
5058                   +ph(i  ,kde,j)+ph(i-1,kde,j) )/g
5059 
5060       ! Find bottom of damping layer
5061       k1 = ktf
5062       z = ztop
5063       DO WHILE( z >= (ztop-zdamp) )
5064         z = 0.25*( phb(i  ,k1,j)+phb(i  ,k1+1,j)  &
5065                   +phb(i-1,k1,j)+phb(i-1,k1+1,j)  &
5066                   +ph(i  ,k1,j)+ph(i  ,k1+1,j)    &
5067                   +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5068         z00(k1) = z
5069         k1 = k1 - 1
5070       ENDDO
5071       k1 = k1 + 2
5072 
5073       ! Get reference state at model levels
5074       DO k = k1, ktf
5075         k2 = ktf
5076         DO WHILE( z_base(k2) .gt. z00(k) )
5077           k2 = k2 - 1
5078         ENDDO
5079         if(k2+1.gt.ktf)then
5080           u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) )   &
5081                               * (     z00(k) - z_base(k2)   )   &
5082                               / ( z_base(k2) - z_base(k2-1) )
5083         else
5084           u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) )   &
5085                               * (       z00(k) - z_base(k2) )   &
5086                               / ( z_base(k2+1) - z_base(k2) )
5087         endif
5088       ENDDO
5089 
5090       ! Apply the Rayleigh damper
5091       DO k = k1, ktf
5092         dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5093         dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5094         ru_tendf(i,k,j) = ru_tendf(i,k,j) -                    &
5095                           muu(i,j) * ( dcoef * dampcoef ) *    &
5096                           ( u(i,k,j) - u00(k) )
5097       END DO
5098 
5099     END DO
5100     END DO
5101 
5102 ! End adjustment of u.
5103 !-----------------------------------------------------------------------
5104 
5105 !-----------------------------------------------------------------------
5106 ! Adjust v to base state.
5107 
5108     DO j = jts, MIN( jte, jde   )
5109     DO i = its, MIN( ite, ide-1 )
5110 
5111       ! Get height at top of model
5112       ztop = 0.5*( phb(i,kde,j  )+phb(i,kde,j-1)   &
5113                   +ph(i,kde,j  )+ph(i,kde,j-1) )/g
5114 
5115       ! Find bottom of damping layer
5116       k1 = ktf
5117       z = ztop
5118       DO WHILE( z >= (ztop-zdamp) )
5119         z = 0.25*( phb(i,k1,j  )+phb(i,k1+1,j  )  &
5120                   +phb(i,k1,j-1)+phb(i,k1+1,j-1)  &
5121                   +ph(i,k1,j  )+ph(i,k1+1,j  )    &
5122                   +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
5123         z00(k1) = z
5124         k1 = k1 - 1
5125       ENDDO
5126       k1 = k1 + 2
5127 
5128       ! Get reference state at model levels
5129       DO k = k1, ktf
5130         k2 = ktf
5131         DO WHILE( z_base(k2) .gt. z00(k) )
5132           k2 = k2 - 1
5133         ENDDO
5134         if(k2+1.gt.ktf)then
5135           v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) )   &
5136                               * (     z00(k) - z_base(k2)   )   &
5137                               / ( z_base(k2) - z_base(k2-1) )
5138         else
5139           v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) )   &
5140                               * (       z00(k) - z_base(k2) )   &
5141                               / ( z_base(k2+1) - z_base(k2) )
5142         endif
5143       ENDDO
5144 
5145       ! Apply the Rayleigh damper
5146       DO k = k1, ktf
5147         dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5148         dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5149         rv_tendf(i,k,j) = rv_tendf(i,k,j) -                    &
5150                           muv(i,j) * ( dcoef * dampcoef ) *    &
5151                           ( v(i,k,j) - v00(k) )
5152       END DO
5153 
5154     END DO
5155     END DO
5156 
5157 ! End adjustment of v.
5158 !-----------------------------------------------------------------------
5159 
5160 !-----------------------------------------------------------------------
5161 ! Adjust w to base state.
5162 
5163     DO j = jts, MIN( jte,   jde-1 )
5164     DO i = its, MIN( ite,   ide-1 )
5165       ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5166       DO k = kts, MIN( kte,   kde   )
5167         z = ( phb(i,k,j) + ph(i,k,j) ) / g
5168         IF ( z >= (ztop-zdamp) ) THEN
5169           dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
5170           dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
5171           rw_tendf(i,k,j) = rw_tendf(i,k,j) -  &
5172                             mut(i,j) * ( dcoef * dampcoef ) * w(i,k,j)
5173         END IF
5174       END DO
5175     END DO
5176     END DO
5177 
5178 ! End adjustment of w.
5179 !-----------------------------------------------------------------------
5180 
5181 !-----------------------------------------------------------------------
5182 ! Adjust potential temperature to base state.
5183 
5184     DO j = jts, MIN( jte,   jde-1 )
5185     DO i = its, MIN( ite,   ide-1 )
5186 
5187       ! Get height at top of model
5188       ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
5189 
5190       ! Find bottom of damping layer
5191       k1 = ktf
5192       z = ztop
5193       DO WHILE( z >= (ztop-zdamp) )
5194         z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) +  &
5195                      ph(i,k1,j) +  ph(i,k1+1,j) ) / g
5196         z00(k1) = z
5197         k1 = k1 - 1
5198       ENDDO
5199       k1 = k1 + 2
5200 
5201       ! Get reference state at model levels
5202       DO k = k1, ktf
5203         k2 = ktf
5204         DO WHILE( z_base(k2) .gt. z00(k) )
5205           k2 = k2 - 1
5206         ENDDO
5207         if(k2+1.gt.ktf)then
5208           t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) )   &
5209                               * (     z00(k) - z_base(k2)   )   &
5210                               / ( z_base(k2) - z_base(k2-1) )
5211         else
5212           t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) )   &
5213                               * (       z00(k) - z_base(k2) )   &
5214                               / ( z_base(k2+1) - z_base(k2) )
5215         endif
5216       ENDDO
5217 
5218       ! Apply the Rayleigh damper
5219       DO k = k1, ktf
5220         dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5221         dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5222         t_tendf(i,k,j) = t_tendf(i,k,j) -                      &
5223                          mut(i,j) * ( dcoef * dampcoef )  *    &
5224                          ( t(i,k,j) - t00(k) )
5225       END DO
5226 
5227     END DO
5228     END DO
5229 
5230 ! End adjustment of potential temperature.
5231 !-----------------------------------------------------------------------
5232 
5233     END SUBROUTINE rk_rayleigh_damp
5234 
5235 !==============================================================================
5236 !==============================================================================
5237                                                                                 
5238       SUBROUTINE sixth_order_diffusion( name, field, tendency, mu, dt,  &
5239                                         config_flags,                   &
5240                                         diff_6th_opt, diff_6th_factor,  &
5241                                         ids, ide, jds, jde, kds, kde,   &
5242                                         ims, ime, jms, jme, kms, kme,   &
5243                                         its, ite, jts, jte, kts, kte )
5244                                                                                 
5245 ! History:       14 Nov 2006   Name of variable changed by Jason Knievel
5246 !                07 Jun 2006   Revised and generalized by Jason Knievel  
5247 !                25 Apr 2005   Original code by Jason Knievel, NCAR
5248                                                                                 
5249 ! Purpose:       Apply 6th-order, monotonic (flux-limited), numerical
5250 !                diffusion to 3-d velocity and to scalars.
5251                                                                                 
5252 ! References:    Ming Xue (MWR Aug 2000)
5253 !                Durran ("Numerical Methods for Wave Equations..." 1999)
5254 !                George Bryan (personal communication)
5255  
5256 !------------------------------------------------------------------------------
5257 ! Begin: Declarations.
5258 
5259     IMPLICIT NONE
5260 
5261     INTEGER, INTENT(IN)  &
5262     :: ids, ide, jds, jde, kds, kde,   &
5263        ims, ime, jms, jme, kms, kme,   &
5264        its, ite, jts, jte, kts, kte
5265  
5266     TYPE(grid_config_rec_type), INTENT(IN)  &
5267     :: config_flags
5268  
5269     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT)  &
5270     :: tendency
5271  
5272     REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN)  &
5273     :: field
5274  
5275     REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN)  &
5276     :: mu
5277  
5278     REAL, INTENT(IN)  &
5279     :: dt
5280 
5281     REAL, INTENT(IN)  &
5282     :: diff_6th_factor
5283 
5284     INTEGER, INTENT(IN)  &
5285     :: diff_6th_opt
5286 
5287     CHARACTER(LEN=1) , INTENT(IN)  &
5288     :: name
5289 
5290     INTEGER  &
5291     :: i, j, k,         &
5292        i_start, i_end,  &
5293        j_start, j_end,  &
5294        k_start, k_end,  &
5295        ktf
5296  
5297     REAL  &
5298     :: dflux_x_p0, dflux_y_p0,  &
5299        dflux_x_p1, dflux_y_p1,  &
5300        tendency_x, tendency_y,  &
5301        mu_avg_p0, mu_avg_p1,    &
5302        diff_6th_coef
5303 
5304     LOGICAL  &
5305     :: specified
5306  
5307 ! End: Declarations.
5308 !------------------------------------------------------------------------------
5309 
5310 !------------------------------------------------------------------------------
5311 ! Begin: Translate the diffusion factor into a diffusion coefficient.  See
5312 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
5313 ! fourth) and for diffusion in two dimensions (not one).  For reference, a
5314 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
5315 ! although application of the flux limiter reduces somewhat the effects of
5316 ! diffusion for a given coefficient.
5317 
5318     diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )  
5319 
5320 ! End: Translate diffusion factor.
5321 !------------------------------------------------------------------------------
5322 
5323 !------------------------------------------------------------------------------
5324 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
5325 ! The halo regions are already filled with values by the time this subroutine
5326 ! is called, which allows the stencil to extend beyond the domains' edges.
5327 
5328     ktf = MIN( kte, kde-1 )
5329 
5330     IF ( name .EQ. 'u' ) THEN
5331 
5332       i_start = its
5333       i_end   = ite
5334       j_start = jts
5335       j_end   = MIN(jde-1,jte)
5336       k_start = kts
5337       k_end   = ktf
5338 
5339     ELSE IF ( name .EQ. 'v' ) THEN
5340  
5341       i_start = its
5342       i_end   = MIN(ide-1,ite)
5343       j_start = jts
5344       j_end   = jte
5345       k_start = kts
5346       k_end   = ktf
5347  
5348     ELSE IF ( name .EQ. 'w' ) THEN
5349 
5350       i_start = its
5351       i_end   = MIN(ide-1,ite)
5352       j_start = jts
5353       j_end   = MIN(jde-1,jte)
5354       k_start = kts+1
5355       k_end   = ktf
5356 
5357     ELSE
5358 
5359       i_start = its
5360       i_end   = MIN(ide-1,ite)
5361       j_start = jts
5362       j_end   = MIN(jde-1,jte)
5363       k_start = kts
5364       k_end   = ktf
5365  
5366     ENDIF
5367 
5368 ! End: Assignment of limits of spatial loops.
5369 !------------------------------------------------------------------------------
5370 
5371 !------------------------------------------------------------------------------
5372 ! Begin: Loop across spatial dimensions.
5373 
5374     DO j = j_start, j_end
5375     DO k = k_start, k_end
5376     DO i = i_start, i_end
5377 
5378 !------------------------------------------------------------------------------
5379 ! Begin: Diffusion in x (i index).
5380  
5381 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
5382  
5383       dflux_x_p0 = (  10.0 * ( field(i,  k,j) - field(i-1,k,j) )    &
5384                      - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) )    &
5385                      +       ( field(i+2,k,j) - field(i-3,k,j) ) )
5386  
5387       dflux_x_p1 = (  10.0 * ( field(i+1,k,j) - field(i  ,k,j) )    &
5388                      - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) )    &
5389                      +       ( field(i+3,k,j) - field(i-2,k,j) ) )
5390  
5391 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5392 ! (variation on Xue's eq. 10).
5393 
5394       IF ( diff_6th_opt .EQ. 2 ) THEN
5395  
5396         IF ( dflux_x_p0 * ( field(i  ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
5397           dflux_x_p0 = 0.0
5398         END IF
5399  
5400         IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i  ,k,j) ) .LE. 0.0 ) THEN
5401           dflux_x_p1 = 0.0
5402         END IF
5403 
5404       END IF
5405 
5406 ! Apply 6th-order diffusion in x direction.
5407  
5408       IF      ( name .EQ. 'u' ) THEN
5409         mu_avg_p0 = mu(i-1,j)
5410         mu_avg_p1 = mu(i  ,j)
5411       ELSE IF ( name .EQ. 'v' ) THEN
5412         mu_avg_p0 = 0.25 * (       &
5413                     mu(i-1,j-1) +  &
5414                     mu(i  ,j-1) +  &
5415                     mu(i-1,j  ) +  &
5416                     mu(i  ,j  ) )
5417         mu_avg_p1 = 0.25 * (       &
5418                     mu(i  ,j-1) +  &
5419                     mu(i+1,j-1) +  &
5420                     mu(i  ,j  ) +  &
5421                     mu(i+1,j  ) )
5422       ELSE
5423         mu_avg_p0 = 0.5 * (        &
5424                     mu(i-1,j) +    &
5425                     mu(i  ,j) )
5426         mu_avg_p1 = 0.5 * (        &
5427                     mu(i  ,j) +    &
5428                     mu(i+1,j) )
5429       END IF
5430  
5431       tendency_x = diff_6th_coef *  &
5432                  ( ( mu_avg_p1 * dflux_x_p1 ) - ( mu_avg_p0 * dflux_x_p0 ) )
5433  
5434 ! End: Diffusion in x.
5435 !------------------------------------------------------------------------------
5436  
5437 !------------------------------------------------------------------------------
5438 ! Begin: Diffusion in y (j index).
5439  
5440 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
5441  
5442       dflux_y_p0 = (  10.0 * ( field(i,k,j  ) - field(i,k,j-1) )    &
5443                      - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) )    &
5444                      +       ( field(i,k,j+2) - field(i,k,j-3) ) )
5445  
5446       dflux_y_p1 = (  10.0 * ( field(i,k,j+1) - field(i,k,j  ) )    &
5447                      - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) )    &
5448                      +       ( field(i,k,j+3) - field(i,k,j-2) ) )
5449  
5450 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
5451 ! (variation on Xue's eq. 10).
5452 
5453       IF ( diff_6th_opt .EQ. 2 ) THEN
5454  
5455         IF ( dflux_y_p0 * ( field(i,k,j  )-field(i,k,j-1) ) .LE. 0.0 ) THEN
5456           dflux_y_p0 = 0.0
5457         END IF
5458  
5459         IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j  ) ) .LE. 0.0 ) THEN
5460           dflux_y_p1 = 0.0
5461         END IF
5462 
5463       END IF
5464  
5465 ! Apply 6th-order diffusion in y direction.
5466  
5467       IF      ( name .EQ. 'u' ) THEN
5468         mu_avg_p0 = 0.25 * (       &
5469                     mu(i-1,j-1) +  &
5470                     mu(i  ,j-1) +  &
5471                     mu(i-1,j  ) +  &
5472                     mu(i  ,j  ) )
5473         mu_avg_p1 = 0.25 * (       &
5474                     mu(i-1,j  ) +  &
5475                     mu(i  ,j  ) +  &
5476                     mu(i-1,j+1) +  &
5477                     mu(i  ,j+1) )
5478       ELSE IF ( name .EQ. 'v' ) THEN
5479         mu_avg_p0 = mu(i,j-1)
5480         mu_avg_p1 = mu(i,j  )
5481       ELSE
5482         mu_avg_p0 = 0.5 * (      &
5483                     mu(i,j-1) +  &
5484                     mu(i,j  ) )
5485         mu_avg_p1 = 0.5 * (      &
5486                     mu(i,j  ) +  &
5487                     mu(i,j+1) )
5488       END IF
5489  
5490       tendency_y = diff_6th_coef *  &
5491                  ( ( mu_avg_p1 * dflux_y_p1 ) - ( mu_avg_p0 * dflux_y_p0 ) )
5492  
5493 ! End: Diffusion in y.
5494 !------------------------------------------------------------------------------
5495  
5496 !------------------------------------------------------------------------------
5497 ! Begin: Combine diffusion in x and y.
5498      
5499       tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
5500  
5501 ! End: Combine diffusion in x and y.
5502 !------------------------------------------------------------------------------
5503 
5504     ENDDO
5505     ENDDO
5506     ENDDO
5507 
5508 ! End: Loop across spatial dimensions.
5509 !------------------------------------------------------------------------------
5510  
5511     END SUBROUTINE sixth_order_diffusion
5512  
5513 !==============================================================================
5514 !==============================================================================
5515 
5516 END MODULE module_big_step_utilities_em