2 !wrf:MODEL_LAYER:DYNAMICS
14 MODULE module_big_step_utilities_em
16 USE module_model_constants
17 USE module_state_description, only: p_qg, p_qs, p_qi, gdscheme, tiedtkescheme, ntiedtkescheme, kfetascheme, mskfscheme, &
18 g3scheme, gfscheme,p_qv, param_first_scalar, p_qr, p_qc, DFI_FWD
19 USE module_configure, ONLY : grid_config_rec_type
24 !-------------------------------------------------------------------------------
26 SUBROUTINE calc_mu_uv ( config_flags, &
28 ids, ide, jds, jde, kds, kde, &
29 ims, ime, jms, jme, kms, kme, &
30 its, ite, jts, jte, kts, kte )
36 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
38 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
39 ims, ime, jms, jme, kms, kme, &
40 its, ite, jts, jte, kts, kte
42 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
43 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
47 INTEGER :: i, j, itf, jtf, im, jm
51 ! calc_mu_uv calculates the full column dry-air mass at the staggered
52 ! horizontal velocity points (u,v) and places the results in muu and muv.
53 ! This routine uses the reference state (mub) and perturbation state (mu)
61 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
64 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
67 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
70 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
75 if(config_flags%periodic_x) im = its-1
77 ! MUU(i,j) = MU(i,j) +MUB(i,j)
78 ! fix for periodic b.c., 13 march 2004, wcs
79 MUU(i,j) = 0.5*(MU(i,j)+MU(im,j)+MUB(i,j)+MUB(im,j))
81 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
84 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
89 if(config_flags%periodic_x) im = ite
91 ! MUU(i,j) = MU(i-1,j) +MUB(i-1,j)
92 ! fix for periodic b.c., 13 march 2004, wcs
93 MUU(i,j) = 0.5*(MU(i-1,j)+MU(im,j)+MUB(i-1,j)+MUB(im,j))
95 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
98 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
103 if(config_flags%periodic_x) im = its-1
105 ! MUU(i,j) = MU(i,j) +MUB(i,j)
106 ! fix for periodic b.c., 13 march 2004, wcs
107 MUU(i,j) = 0.5*(MU(i,j)+MU(im,j)+MUB(i,j)+MUB(im,j))
111 if(config_flags%periodic_x) im = ite
113 ! MUU(i,j) = MU(i-1,j) +MUB(i-1,j)
114 ! fix for periodic b.c., 13 march 2004, wcs
115 MUU(i,j) = 0.5*(MU(i-1,j)+MU(im,j)+MUB(i-1,j)+MUB(im,j))
122 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
125 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
128 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
131 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
136 if(config_flags%periodic_y) jm = jts-1
138 ! MUV(i,j) = MU(i,j) +MUB(i,j)
139 ! fix for periodic b.c., 13 march 2004, wcs
140 MUV(i,j) = 0.5*(MU(i,j)+MU(i,jm)+MUB(i,j)+MUB(i,jm))
142 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
145 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
150 if(config_flags%periodic_y) jm = jte
152 ! MUV(i,j) = MU(i,j-1) +MUB(i,j-1)
153 ! fix for periodic b.c., 13 march 2004, wcs
154 MUV(i,j) = 0.5*(MU(i,j-1)+MU(i,jm)+MUB(i,j-1)+MUB(i,jm))
156 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
159 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
164 if(config_flags%periodic_y) jm = jts-1
166 ! MUV(i,j) = MU(i,j) +MUB(i,j)
167 ! fix for periodic b.c., 13 march 2004, wcs
168 MUV(i,j) = 0.5*(MU(i,j)+MU(i,jm)+MUB(i,j)+MUB(i,jm))
172 if(config_flags%periodic_y) jm = jte
174 ! MUV(i,j) = MU(i,j-1) +MUB(i,j-1)
175 ! fix for periodic b.c., 13 march 2004, wcs
176 MUV(i,j) = 0.5*(MU(i,j-1)+MU(i,jm)+MUB(i,j-1)+MUB(i,jm))
180 END SUBROUTINE calc_mu_uv
182 !-------------------------------------------------------------------------------
184 SUBROUTINE calc_mu_uv_1 ( config_flags, &
186 ids, ide, jds, jde, kds, kde, &
187 ims, ime, jms, jme, kms, kme, &
188 its, ite, jts, jte, kts, kte )
194 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
196 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
197 ims, ime, jms, jme, kms, kme, &
198 its, ite, jts, jte, kts, kte
200 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
201 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
205 INTEGER :: i, j, itf, jtf, im, jm
209 ! calc_mu_uv calculates the full column dry-air mass at the staggered
210 ! horizontal velocity points (u,v) and places the results in muu and muv.
211 ! This routine uses the full state (mu)
218 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
221 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j))
224 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
227 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j))
232 if(config_flags%periodic_x) im = its-1
234 MUU(i,j) = 0.5*(MU(i,j)+MU(im,j))
236 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
239 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j))
244 if(config_flags%periodic_x) im = ite
246 MUU(i,j) = 0.5*(MU(i-1,j)+MU(im,j))
248 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
251 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j))
256 if(config_flags%periodic_x) im = its-1
258 MUU(i,j) = 0.5*(MU(i,j)+MU(im,j))
262 if(config_flags%periodic_x) im = ite
264 MUU(i,j) = 0.5*(MU(i-1,j)+MU(im,j))
271 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
274 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1))
277 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
280 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1))
285 if(config_flags%periodic_y) jm = jts-1
287 MUV(i,j) = 0.5*(MU(i,j)+MU(i,jm))
289 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
292 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1))
297 if(config_flags%periodic_y) jm = jte
299 MUV(i,j) = 0.5*(MU(i,j-1)+MU(i,jm))
301 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
304 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1))
309 if(config_flags%periodic_y) jm = jts-1
311 MUV(i,j) = 0.5*(MU(i,j)+MU(i,jm))
315 if(config_flags%periodic_y) jm = jte
317 MUV(i,j) = 0.5*(MU(i,j-1)+MU(i,jm))
321 END SUBROUTINE calc_mu_uv_1
323 !-------------------------------------------------------------------------------
325 ! Map scale factor comments for this routine:
326 ! Locally not changed, but sent the correct map scale factors
327 ! from module_em (msfuy, msfvx, msfty)
329 SUBROUTINE couple_momentum ( muu, ru, u, msfu, &
330 muv, rv, v, msfv, msfv_inv, &
332 c1h, c2h, c1f, c2f, &
333 ids, ide, jds, jde, kds, kde, &
334 ims, ime, jms, jme, kms, kme, &
335 its, ite, jts, jte, kts, kte )
341 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
342 ims, ime, jms, jme, kms, kme, &
343 its, ite, jts, jte, kts, kte
345 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: ru, rv, rw
347 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, muv, mut
348 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, msfv, msft, msfv_inv
350 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v, w
351 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1f, c2f
355 INTEGER :: i, j, k, itf, jtf, ktf
359 ! couple_momentum couples the velocities to the full column mass and
372 ru(i,k,j)=u(i,k,j)*(c1h(k)*muu(i,j)+c2h(k))/msfu(i,j)
383 rv(i,k,j)=v(i,k,j)*(c1h(k)*muv(i,j)+c2h(k))*msfv_inv(i,j)
394 rw(i,k,j)=w(i,k,j)*(c1f(k)*mut(i,j)+c2f(k))/msft(i,j)
399 END SUBROUTINE couple_momentum
401 !-------------------------------------------------------------------
403 SUBROUTINE calc_mu_staggered ( mu, mub, muu, muv, &
404 ids, ide, jds, jde, kds, kde, &
405 ims, ime, jms, jme, kms, kme, &
406 its, ite, jts, jte, kts, kte )
412 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
413 ims, ime, jms, jme, kms, kme, &
414 its, ite, jts, jte, kts, kte
416 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv
417 REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub
421 INTEGER :: i, j, itf, jtf
425 ! calc_mu_staggered calculates the full dry air mass at the staggered
426 ! velocity points (u,v).
433 IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN
436 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
439 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN
442 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
447 MUU(i,j) = MU(i,j) +MUB(i,j)
449 ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN
452 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
457 MUU(i,j) = MU(i-1,j) +MUB(i-1,j)
459 ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN
462 MUU(i,j) = 0.5*(MU(i,j)+MU(i-1,j)+MUB(i,j)+MUB(i-1,j))
467 MUU(i,j) = MU(i,j) +MUB(i,j)
471 MUU(i,j) = MU(i-1,j) +MUB(i-1,j)
478 IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN
481 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
484 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN
487 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
492 MUV(i,j) = MU(i,j) +MUB(i,j)
494 ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN
497 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
502 MUV(i,j) = MU(i,j-1) +MUB(i,j-1)
504 ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN
507 MUV(i,j) = 0.5*(MU(i,j)+MU(i,j-1)+MUB(i,j)+MUB(i,j-1))
512 MUV(i,j) = MU(i,j) +MUB(i,j)
516 MUV(i,j) = MU(i,j-1) +MUB(i,j-1)
520 END SUBROUTINE calc_mu_staggered
522 !-------------------------------------------------------------------------------
524 SUBROUTINE couple ( mu, mub, rfield, field, name, &
525 msf, c1h, c2h, c1, c2, &
526 ids, ide, jds, jde, kds, kde, &
527 ims, ime, jms, jme, kms, kme, &
528 its, ite, jts, jte, kts, kte )
534 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
535 ims, ime, jms, jme, kms, kme, &
536 its, ite, jts, jte, kts, kte
538 CHARACTER(LEN=1) , INTENT(IN ) :: name
540 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield
542 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf
544 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field
546 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1, c2
550 INTEGER :: i, j, k, itf, jtf, ktf
551 REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv
555 ! subroutine couple couples the input variable with the dry-air
563 IF (name .EQ. 'u')THEN
565 CALL calc_mu_staggered ( mu, mub, muu, muv, &
566 ids, ide, jds, jde, kds, kde, &
567 ims, ime, jms, jme, kms, kme, &
568 its, ite, jts, jte, kts, kte )
576 rfield(i,k,j)=field(i,k,j)*(c1h(k)*muu(i,j)+c2h(k))/msf(i,j)
581 ELSE IF (name .EQ. 'v')THEN
583 CALL calc_mu_staggered ( mu, mub, muu, muv, &
584 ids, ide, jds, jde, kds, kde, &
585 ims, ime, jms, jme, kms, kme, &
586 its, ite, jts, jte, kts, kte )
595 rfield(i,k,j)=field(i,k,j)*(c1h(k)*muv(i,j)+c2h(k))/msf(i,j)
600 ELSE IF (name .EQ. 'w')THEN
606 rfield(i,k,j)=field(i,k,j)*((c1(k)*mu(i,j))+(c1(k)*mub(i,j)+c2(k)))/msf(i,j)
611 ELSE IF (name .EQ. 'h')THEN
617 rfield(i,k,j)=field(i,k,j)*((c1(k)*mu(i,j))+(c1(k)*mub(i,j)+c2(k)))
628 rfield(i,k,j)=field(i,k,j)*((c1(k)*mu(i,j))+(c1(k)*mub(i,j)+c2(k)))
635 END SUBROUTINE couple
638 !-------------------------------------------------------------------------------
640 SUBROUTINE calc_ww_cp ( u, v, mup, mub, c1h, c2h, ww, &
641 rdx, rdy, msftx, msfty, &
642 msfux, msfuy, msfvx, msfvx_inv, &
644 ids, ide, jds, jde, kds, kde, &
645 ims, ime, jms, jme, kms, kme, &
646 its, ite, jts, jte, kts, kte )
653 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
654 ims, ime, jms, jme, kms, kme, &
655 its, ite, jts, jte, kts, kte
657 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v
658 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, &
663 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
664 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h
666 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww
667 REAL , INTENT(IN ) :: rdx, rdy
671 INTEGER :: i, j, k, itf, jtf, ktf
672 REAL , DIMENSION( its:ite ) :: dmdt
673 REAL , DIMENSION( its:ite, kts:kte ) :: divv
674 REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv
678 ! calc_ww calculates omega using the velocities (u,v) and the dry-air
679 ! column mass (mup+mub).
680 ! The algorithm integrates the continuity equation through the column
681 ! followed by a diagnosis of omega.
687 ! calc_ww_cp calculates omega using the velocities (u,v) and the
696 ! mu coupled with the appropriate map factor
699 DO i=its,min(ite+1,ide)
700 MUU(i,j) = 0.5*(MUP(i,j)+MUB(i,j)+MUP(i-1,j)+MUB(i-1,j))
704 DO j=jts,min(jte+1,jde)
706 MUV(i,j) = 0.5*(MUP(i,j)+MUB(i,j)+MUP(i,j-1)+MUB(i,j-1))
718 ! Comments on the modifications for map scale factors
719 ! ADT eqn 47 / my (putting rho -> 'mu') is:
720 ! (1/my) partial d mu/dt = -mx partial d/dx(mu u/my)
721 ! -mx partial d/dy(mu v/mx)
722 ! -partial d/dz(mu w/my)
724 ! Using nu instead of z the last term becomes:
725 ! -partial d/dnu((c1(k)*mu(dnu/dt))/my)
727 ! Integrating with respect to nu over ALL levels, with dnu/dt=0 at top
728 ! and bottom, the last term becomes = 0
730 ! Integral|bot->top[(1/my) partial d mu/dt]dnu =
731 ! Integral|bot->top[-mx partial d/dx(mu u/my)
732 ! -mx partial d/dy(mu v/mx)]dnu
734 ! muu='mu'[on u]/my, muv='mu'[on v]/mx
735 ! (1/my) partial d mu/dt is independent of nu
736 ! => LHS = Integral|bot->top[con]dnu = conservation*(-1) = -dmdt
738 ! => dmdt = mx*Integral|bot->top[partial d/dx(mu u/my) +
739 ! partial d/dy(mu v/mx)]dnu
740 ! => dmdt = sum_bot->top[divv]
742 ! divv=mx*[partial d/dx(mu u/my) + partial d/dy(mu v/mx)]*delta nu
747 divv(i,k) = msftx(i,j)*dnw(k)*( rdx*((c1h(k)*muu(i+1,j)+c2h(k))*u(i+1,k,j)/msfuy(i+1,j)-(c1h(k)*muu(i,j)+c2h(k))*u(i,k,j)/msfuy(i,j)) &
748 +rdy*((c1h(k)*muv(i,j+1)+c2h(k))*v(i,k,j+1)*msfvx_inv(i,j+1)-(c1h(k)*muv(i,j)+c2h(k))*v(i,k,j)*msfvx_inv(i,j)) )
750 ! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) &
751 ! +rdy*(rv(i,k,j+1)-rv(i,k,j)) )
753 dmdt(i) = dmdt(i) + divv(i,k)
759 ! Further map scale factor notes:
760 ! Now integrate from bottom to top, level by level:
761 ! mu dnu/dt/my [k+1] = mu dnu/dt/my [k] + [-(1/my) partial d mu/dt
762 ! -mx partial d/dx(mu u/my)
763 ! -mx partial d/dy(mu v/mx)]*dnu[k->k+1]
764 ! ww [k+1] = ww [k] -(1/my) partial d mu/dt * dnu[k->k+1] - divv[k]
765 ! = ww [k] -dmdt * dnw[k] - divv[k]
770 ! ww(i,k,j)=ww(i,k-1,j) &
771 ! - dnw(k-1)* ( dmdt(i) &
772 ! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) &
773 ! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) )
775 ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*c1h(k-1)*dmdt(i) - divv(i,k-1)
782 END SUBROUTINE calc_ww_cp
785 !-------------------------------------------------------------------------------
787 SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, &
788 ids, ide, jds, jde, kds, kde, &
789 ims, ime, jms, jme, kms, kme, &
790 its, ite, jts, jte, kts, kte )
796 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
797 ims, ime, jms, jme, kms, kme, &
798 its, ite, jts, jte, kts, kte
800 INTEGER , INTENT(IN ) :: n_moist
803 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist
805 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw
809 ! Changes from Larry Meadows, Intel Corp. Improve vectorization of this routine
810 REAL :: qtot(its:ite)
812 INTEGER :: i, j, k, itf, jtf, ktf, ispe
816 ! calc_cq calculates moist coefficients for the momentum equations.
822 IF( n_moist >= PARAM_FIRST_SCALAR ) THEN
829 DO ispe=PARAM_FIRST_SCALAR,n_moist
831 qtot(i) = qtot(i) + moist(i,k,j,ispe) + moist(i-1,k,j,ispe)
835 cqu(i,k,j) = 1./(1.+0.5*qtot(i))
845 DO ispe=PARAM_FIRST_SCALAR,n_moist
847 qtot(i) = qtot(i) + moist(i,k,j,ispe) + moist(i,k,j-1,ispe)
851 cqv(i,k,j) = 1./(1.+0.5*qtot(i))
861 DO ispe=PARAM_FIRST_SCALAR,n_moist
863 qtot(i) = qtot(i) + moist(i,k,j,ispe) + moist(i,k-1,j,ispe)
867 cqw(i,k,j) = 0.5*qtot(i)
906 END SUBROUTINE calc_cq
908 !----------------------------------------------------------------------
910 SUBROUTINE calc_alt ( alt, al, alb, &
911 ids, ide, jds, jde, kds, kde, &
912 ims, ime, jms, jme, kms, kme, &
913 its, ite, jts, jte, kts, kte )
919 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
920 ims, ime, jms, jme, kms, kme, &
921 its, ite, jts, jte, kts, kte
923 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al
924 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt
928 INTEGER :: i, j, k, itf, jtf, ktf
932 ! calc_alt computes the full inverse density
943 alt(i,k,j) = al(i,k,j)+alb(i,k,j)
949 END SUBROUTINE calc_alt
951 !----------------------------------------------------------------------
953 SUBROUTINE calc_p_rho_phi ( moist, n_moist, hypsometric_opt, &
955 c1, c2, c3h, c4h, c3f, c4f, &
957 t, p0, t0, ptop, znu, znw, dnw, rdnw, &
958 rdn, non_hydrostatic, use_theta_m, &
959 ids, ide, jds, jde, kds, kde, &
960 ims, ime, jms, jme, kms, kme, &
961 its, ite, jts, jte, kts, kte )
967 LOGICAL , INTENT(IN ) :: non_hydrostatic
968 INTEGER , INTENT(IN ) :: use_theta_m
970 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
971 ims, ime, jms, jme, kms, kme, &
972 its, ite, jts, jte, kts, kte
974 INTEGER , INTENT(IN ) :: n_moist
975 INTEGER , INTENT(IN ) :: hypsometric_opt
977 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, &
981 REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist
983 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p
985 REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph, phb
987 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts
989 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, znw, dnw, rdnw, rdn
991 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: c1, c2, c3h, c4h, c3f, c4f
993 REAL, INTENT(IN ) :: t0, p0, ptop
997 INTEGER :: i, j, k, kk, itf, jtf, ktf, ispe
998 REAL :: qtot, qf1, qf2, qvf
999 REAL, DIMENSION( its:ite) :: temp,cpovcv_v
1000 REAL :: pfu, phm, pfd
1004 ! For the nonhydrostatic option, calc_p_rho_phi calculates the
1005 ! diagnostic quantities pressure and (inverse) density from the
1006 ! prognostic variables using the equation of state.
1008 ! For the hydrostatic option, calc_p_rho_phi calculates the
1009 ! diagnostic quantities (inverse) density and geopotential from the
1010 ! prognostic variables using the equation of state and the hydrostatic
1023 nonhydro : IF (non_hydrostatic) THEN
1025 hypso_nonhydro : IF (hypsometric_opt == 1) THEN
1029 al(i,k,j)=-1./(c1(k)*muts(i,j)+c2(k))*(alb(i,k,j)*(c1(k)*mu(i,j)) + rdnw(k)*(ph(i,k+1,j)-ph(i,k,j)))
1033 ELSE IF (hypsometric_opt == 2) THEN hypso_nonhydro
1035 ! The relation used to get specific volume, al, is: al = -dZ/dp,
1036 ! where dp = mut * d(eta). The pressure depth, dp, is replaced with
1037 ! p*(dp/p) ~ p*LOG((p+0.5dp)/(p-0.5dp)). Difference between dp and p*dLOG(p)
1038 ! is as follows: p*dLOG(p) - dp = 1/12*(dp/p)**3 + 1/90*(dp/p)**5 + ...
1039 ! Therefore, p*dLOG(p) is always larger than dp and the difference is
1040 ! in proportion to dp/p. TKW, 02/16/2010
1045 pfu = c3f(k+1)*MUTS(i,j) + c4f(k+1) + ptop
1046 pfd = c3f(k )*MUTS(i,j) + c4f(k ) + ptop
1047 phm = c3h(k )*MUTS(i,j) + c4h(k ) + ptop
1048 al(i,k,j) = (ph(i,k+1,j)-ph(i,k,j)+phb(i,k+1,j)-phb(i,k,j))/phm/LOG(pfd/pfu)-alb(i,k,j)
1053 CALL wrf_error_fatal ( 'calc_p_rho_phi: hypsometric_opt should be 1 or 2' )
1054 END IF hypso_nonhydro
1056 moist_nonhydro : IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1061 IF ( use_theta_m .EQ. 1 ) THEN
1062 temp(i)=(r_d*(t0+t(i,k,j)))/(p0*(al(i,k,j)+alb(i,k,j)))
1064 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1065 temp(i)=(r_d*(t0+t(i,k,j))*qvf)/(p0*(al(i,k,j)+alb(i,k,j)))
1069 CALL VPOWX ( itf-its+1, temp(its), cpovcv, p(its,k,j) )
1071 ! use vector version from libmassv or from compat lib in frame/libmassv.F
1072 CALL VPOW ( p(its,k,j), temp(its), cpovcv_v(its), itf-its+1 )
1075 p(i,k,j)= p(i,k,j)*p0-pb(i,k,j)
1085 p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/ &
1086 (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv &
1092 END IF moist_nonhydro
1096 ! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001
1099 moist_hydro : IF (n_moist >= PARAM_FIRST_SCALAR ) THEN
1107 DO ispe=PARAM_FIRST_SCALAR,n_moist
1108 qtot = qtot + moist(i,k,j,ispe)
1113 p(i,k,j) = - 0.5*((c1(k)*mu(i,j))+qf1*(c1(k)*muts(i,j)+c2(k)))/rdnw(k)/qf2
1114 IF ( use_theta_m .EQ. 1 ) THEN
1115 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)* &
1116 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1118 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1119 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1120 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1125 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1129 DO ispe=PARAM_FIRST_SCALAR,n_moist
1130 qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) )
1135 p(i,k,j) = p(i,k+1,j) - ((c1(k)*mu(i,j)) + qf1*(c1(k)*muts(i,j)+c2(k)))/qf2/rdn(k+1)
1136 IF ( use_theta_m .EQ. 1 ) THEN
1137 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)* &
1138 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1140 qvf = 1.+rvovrd*moist(i,k,j,P_QV)
1141 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* &
1142 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1160 p(i,k,j) = - 0.5*((c1(k)*mu(i,j))+qf1*(c1(k)*muts(i,j)+c2(k)))/rdnw(k)/qf2
1161 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)* &
1162 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1166 DO k=ktf-1,kts,-1 ! remaining layers, integrate down
1173 p(i,k,j) = p(i,k+1,j) - ((c1(k)*mu(i,j)) + qf1*(c1(k)*muts(i,j)+c2(k)))/qf2/rdn(k+1)
1174 al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)* &
1175 (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j)
1183 hypso_hydro : IF (hypsometric_opt == 1) THEN
1185 DO kk=2,ktf+1 ! integrate hydrostatic equation for geopotential
1188 ph(i,k+1,j) = ph(i,k,j) - (dnw(k))*( &
1189 ((c1(k)*muts(i,j)+c2(k)))*al(i,k,j)+ &
1190 (c1(k)*mu(i,j))*alb(i,k,j) )
1196 ! Revised hypsometric eq.: dZ=-al*p*dLOG(p), where p is dry pressure
1200 ph(i,kts,j) = phb(i,kts,j)
1205 pfu = c3f(k )*MUTS(i,j) + c4f(k ) + ptop
1206 pfd = c3f(k-1)*MUTS(i,j) + c4f(k-1) + ptop
1207 phm = c3h(k-1)*MUTS(i,j) + c4h(k-1) + ptop
1208 ph(i,k,j) = ph(i,k-1,j) + (al(i,k-1,j)+alb(i,k-1,j))*phm*LOG(pfd/pfu)
1214 ph(i,k,j) = ph(i,k,j) - phb(i,k,j)
1223 END SUBROUTINE calc_p_rho_phi
1225 !----------------------------------------------------------------------
1227 SUBROUTINE calc_php ( php, ph, phb, &
1228 ids, ide, jds, jde, kds, kde, &
1229 ims, ime, jms, jme, kms, kme, &
1230 its, ite, jts, jte, kts, kte )
1236 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1237 ims, ime, jms, jme, kms, kme, &
1238 its, ite, jts, jte, kts, kte
1240 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph
1241 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php
1245 INTEGER :: i, j, k, itf, jtf, ktf
1249 ! calc_php calculates the full geopotential from the reference state
1250 ! geopotential and the perturbation geopotential (phb_ph).
1261 php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j))
1266 END SUBROUTINE calc_php
1268 !-------------------------------------------------------------------------------
1270 SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mut, &
1273 cf1, cf2, cf3, rdx, rdy, &
1275 ids, ide, jds, jde, kds, kde, &
1276 ims, ime, jms, jme, kms, kme, &
1277 its, ite, jts, jte, kts, kte )
1281 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1282 ims, ime, jms, jme, kms, kme, &
1283 its, ite, jts, jte, kts, kte
1285 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, &
1292 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w
1294 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut, ht, msftx, msfty
1296 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: c1f, c2f
1298 REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy
1300 INTEGER :: i, j, k, itf, jtf
1307 ! diagnose_w diagnoses the vertical velocity from the geopoential equation.
1308 ! Used with the hydrostatic option.
1316 ! Notes on map scale factors:
1317 ! Chain rule: if Z=Z(X,Y) [true at the surface] then
1318 ! dZ/dt = dZ/dX * dX/dt + dZ/dY * dY/dt, U=dX/dt, V=dY/dt
1319 ! Using capitals to denote actual values
1320 ! In mapped values, u=U, v=V, z=Z, 1/dX=mx/dx, 1/dY=my/dy
1321 ! => w = dz/dt = mx u dz/dx + my v dz/dy
1322 ! [where dz/dx is just the surface height change between x
1323 ! gridpoints, and dz/dy is the change between y gridpoints]
1324 ! [NB: cf1, cf2 and cf3 do vertical weighting of u or v values
1325 ! nearest the surface]
1327 ! Previously msft multiplied by rdy and rdx terms.
1328 ! Now msfty multiplies rdy term, and msftx multiplies msftx term
1330 w(i,1,j)= msfty(i,j)*.5*rdy*( &
1331 (ht(i,j+1)-ht(i,j )) &
1332 *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1)) &
1333 +(ht(i,j )-ht(i,j-1)) &
1334 *(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )) ) &
1335 +msftx(i,j)*.5*rdx*( &
1336 (ht(i+1,j)-ht(i,j )) &
1337 *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j)) &
1338 +(ht(i,j )-ht(i-1,j)) &
1339 *(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)) )
1342 ! use geopotential equation to diagnose w
1344 ! Further notes on map scale factors
1345 ! If ph_tend contains: -mx partial d/dx(mu rho u/my)
1346 ! -mx partial d/dy(phi mu v/mx)
1347 ! -partial d/dz(phi mu w/my)
1348 ! then phi eqn is: partial d/dt(mu phi/my) = ph_tend + mu g w/my
1349 ! => w = [my/(mu*g)]*[partial d/dt(mu phi/my) - ph_tend]
1353 w(i,k,j) = msfty(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt &
1354 - ph_tend(i,k,j)/(c1f(k)*mut(i,j)+c2f(k)) )/g
1361 END SUBROUTINE diagnose_w
1363 !-------------------------------------------------------------------------------
1365 SUBROUTINE rhs_ph( ph_tend, u, v, ww, &
1366 ph, ph_old, phb, w, &
1370 rdnw, cfn, cfn1, rdx, rdy, &
1371 msfux, msfuy, msfvx, &
1376 ids, ide, jds, jde, kds, kde, &
1377 ims, ime, jms, jme, kms, kme, &
1378 its, ite, jts, jte, kts, kte )
1381 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
1383 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
1384 ims, ime, jms, jme, kms, kme, &
1385 its, ite, jts, jte, kts, kte
1387 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
1396 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: ph_tend
1398 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muuf, muvf, mut, &
1404 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
1406 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: c1f, c2f
1409 REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy
1411 LOGICAL, INTENT(IN ) :: non_hydrostatic
1415 INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start
1416 REAL :: ur, ul, ub, vr, vl, vb
1417 REAL, DIMENSION(its:ite,kts:kte) :: wdwn
1419 INTEGER :: advective_order
1421 LOGICAL :: specified
1425 ! rhs_ph calculates the large-timestep tendency terms for the geopotential
1426 ! equation. These terms include the advection and "gw". The geopotential
1427 ! equation is cast in advective form, so we don't use the flux form advection
1433 if(config_flags%specified .or. config_flags%nested) specified = .true.
1435 advective_order = config_flags%h_sca_adv_order
1441 ! Notes on map scale factors (WCS, 2 march 2008)
1442 ! phi equation is: mu/my d/dt(phi) = -(1/my) mx mu u d/dx(phi)
1443 ! -(1/my) my mu v d/dy(phi)
1444 ! - omega d/d_eta(phi)
1447 ! A little further explanation...
1448 ! The tendency term we are computing here is for mu/my d/dt(phi). It is advective form
1449 ! but it is multiplied be mu/my. It will be decoupled from (mu/my) when the implicit w-phi
1450 ! solution is computed in subourine advance_w. The formulation dates from the early
1451 ! days of the mass coordinate model when we were testing both a flux and an advective formulation
1452 ! for the geopotential equation and different forms of the vertical momentum equation and the
1453 ! vertically implicit solver.
1455 ! advective form for the geopotential equation
1457 ! RHS term 3 is: - omega partial d/dnu(phi)
1460 IF (config_flags%phi_adv_z == 2) THEN
1461 ! First get staggered partial d/dnu(phi) and then multiply with omega
1462 ! to avoid double staggering of omega
1466 wdwn(i,k) = rdnw(k-1)*(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1472 ph_tend(i,k,j) = ph_tend(i,k,j) &
1473 - ww(i,k,j)*(fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1477 ! First destagger omega and multiply with partial d/dnu(phi), then stagger the product
1481 wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) &
1482 *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j))
1488 ph_tend(i,k,j) = ph_tend(i,k,j) &
1489 - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k))
1496 IF (non_hydrostatic) THEN ! add in "gw" term.
1497 DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed
1498 ! after the timestep to give us "w"
1500 ph_tend(i,kde,j) = 0.
1505 ! phi equation RHS term 4
1506 ph_tend(i,k,j) = ph_tend(i,k,j) + (c1f(k)*mut(i,j)+c2f(k))*g*w(i,k,j)/msfty(i,j)
1514 ! Notes on map scale factors:
1515 ! RHS terms 1 and 2 are: -(1/my) mx u mu partial d/dx(phi)
1516 ! -(1/my) my v mu partial d/dy(phi)
1518 IF (advective_order <= 2) THEN
1527 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+1
1528 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jte-2
1534 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1535 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1536 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1537 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1538 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1544 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1545 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1546 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1547 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1548 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1560 IF ( (config_flags%open_xs .or. specified) .and. its == ids ) i_start = its+1
1561 IF ( (config_flags%open_xe .or. specified) .and. ite == ide ) itf = ite-2
1567 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1568 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1569 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1570 +(c1f(k)*muuf(i ,j)+c2f(k))*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1571 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1577 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1578 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1579 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1580 +(c1f(k)*muuf(i ,j)+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1581 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1586 ELSE IF (advective_order <= 4) THEN
1595 IF ( (config_flags%open_ys .or. specified) .and. jts == jds ) j_start = jts+2
1596 IF ( (config_flags%open_ye .or. specified) .and. jte == jde ) jtf = jte-3
1602 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))*( &
1603 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1604 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ))* (1./12.)*( &
1605 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1606 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1607 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1608 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1616 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))*( &
1617 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1618 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ))* (1./12.)*( &
1619 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1620 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1621 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1622 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1628 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1630 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 ) THEN
1635 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1636 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1637 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1638 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1639 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1645 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1646 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1647 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1648 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1649 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1654 IF ( (config_flags%open_ye .or. specified) .and. jte >= jde-2 ) THEN
1659 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1660 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1661 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1662 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1663 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1669 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1670 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1671 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1672 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1673 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1685 IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+2
1686 IF ( (config_flags%open_xe) .and. ite == ide ) itf = ite-3
1692 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1693 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1694 +(c1f(k)*muuf(i,j )+c2f(k))*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j) )* (1./12.)*( &
1695 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1696 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1697 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1698 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1704 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1705 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1706 +(c1f(k)*muuf(i,j )+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux(i ,j) )* (1./12.)*( &
1707 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1708 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1709 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1710 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1715 ! pick up near boundary rows using 2nd order stencil for open and specified conditions
1717 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 ) THEN
1723 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1724 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1725 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1726 +(c1f(k)*muuf(i ,j)+c2f(k))*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1727 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1733 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1734 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1735 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1736 +(c1f(k)*muuf(i ,j)+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1737 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1742 IF ( (config_flags%open_xe .or. specified) .and. ite >= ide-2 ) THEN
1747 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
1748 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
1749 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1750 +(c1f(k)*muuf(i ,j)+c2f(k))*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
1751 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1757 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
1758 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
1759 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
1760 +(c1f(k)*muuf(i ,j)+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
1761 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
1766 !--------------------------
1768 ELSE IF (advective_order <= 6) THEN
1770 !!! NOTE: this branch has been looked at and fixed with changes for overdecomposition
1771 !!! the branches covering the other advective_order cases have not. 20090923. JM
1780 IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+3)
1781 IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-4)
1787 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1788 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1789 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./60.)*( &
1790 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1791 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1792 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1793 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1794 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1795 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1803 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1804 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1805 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j ) )* (1./60.)*( &
1806 45.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1807 -9.*(ph(i,k,j+2)-ph(i,k,j-2)) &
1808 +(ph(i,k,j+3)-ph(i,k,j-3)) &
1809 +45.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1810 -9.*(phb(i,k,j+2)-phb(i,k,j-2)) &
1811 +(phb(i,k,j+3)-phb(i,k,j-3)) ) )
1817 ! 4th order stencil for open or specified conditions two in form the boundary
1819 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+2 .and. jds+2 <= jte ) THEN
1824 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1825 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1826 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j ) )* (1./12.)*( &
1827 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1828 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1829 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1830 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1838 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1839 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1840 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1841 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1842 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1843 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1844 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1850 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-3 .and. jde-3 <= jte ) THEN
1854 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* ( &
1855 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1) &
1856 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j) )* (1./12.)*( &
1857 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1858 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1859 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1860 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1868 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* ( &
1869 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1) &
1870 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j) )* (1./12.)*( &
1871 8.*(ph(i,k,j+1)-ph(i,k,j-1)) &
1872 -(ph(i,k,j+2)-ph(i,k,j-2)) &
1873 +8.*(phb(i,k,j+1)-phb(i,k,j-1)) &
1874 -(phb(i,k,j+2)-phb(i,k,j-2)) ) )
1880 ! 2nd order stencil for open and specified conditions one row in from boundary
1882 IF ( (config_flags%open_ys .or. specified) .and. jts <= jds+1 .and. jds+1 <= jte ) THEN
1887 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1888 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1889 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1890 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1891 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1897 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1898 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1899 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1900 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1901 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1906 IF ( (config_flags%open_ye .or. specified) .and. jts <= jde-2 .and. jde-2 <= jte ) THEN
1911 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdy/msfty(i,j))* &
1912 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(v(i,k,j+1)+v(i,k-1,j+1))*msfvy(i,j+1)* &
1913 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1914 +(c1f(k)*muvf(i,j )+c2f(k))*(v(i,k,j )+v(i,k-1,j ))*msfvy(i,j )* &
1915 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1921 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdy/msfty(i,j))* &
1922 ( (c1f(k)*muvf(i,j+1)+c2f(k))*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))*msfvy(i,j+1)* &
1923 (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) &
1924 +(c1f(k)*muvf(i,j )+c2f(k))*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))*msfvy(i,j )* &
1925 (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) )
1937 IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+3)
1938 IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-4)
1944 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1945 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1946 +(c1f(k)*muuf(i,j )+c2f(k))*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./60.)*( &
1947 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1948 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1949 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1950 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1951 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1952 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1958 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1959 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1960 +(c1f(k)*muuf(i,j )+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./60.)*( &
1961 45.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1962 -9.*(ph(i+2,k,j)-ph(i-2,k,j)) &
1963 +(ph(i+3,k,j)-ph(i-3,k,j)) &
1964 +45.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1965 -9.*(phb(i+2,k,j)-phb(i-2,k,j)) &
1966 +(phb(i+3,k,j)-phb(i-3,k,j)) ) )
1971 ! 4th order stencil two in from the boundary for open and specified conditions
1973 IF ( (config_flags%open_xs) .and. its <= ids+2 .and. ids+2 <= ite ) THEN
1977 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
1978 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
1979 +(c1f(k)*muuf(i,j )+c2f(k))*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
1980 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1981 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1982 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1983 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1986 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
1987 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
1988 +(c1f(k)*muuf(i,j )+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
1989 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
1990 -(ph(i+2,k,j)-ph(i-2,k,j)) &
1991 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
1992 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
1997 IF ( (config_flags%open_xe) .and. its <= ide-3 .and. ide-3 <= ite ) THEN
2001 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))*( &
2002 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j) &
2003 +(c1f(k)*muuf(i,j )+c2f(k))*(u(i,k,j )+u(i,k-1,j ))*msfux(i,j) )* (1./12.)*( &
2004 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
2005 -(ph(i+2,k,j)-ph(i-2,k,j)) &
2006 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
2007 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
2010 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))*( &
2011 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j) &
2012 +(c1f(k)*muuf(i,j )+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j))*msfux(i,j) )* (1./12.)*( &
2013 8.*(ph(i+1,k,j)-ph(i-1,k,j)) &
2014 -(ph(i+2,k,j)-ph(i-2,k,j)) &
2015 +8.*(phb(i+1,k,j)-phb(i-1,k,j)) &
2016 -(phb(i+2,k,j)-phb(i-2,k,j)) ) )
2021 ! 2nd order stencil for open and specified conditions one in from the boundary
2023 IF ( (config_flags%open_xs .or. specified) .and. its <= ids+1 .and. ids+1 <= ite ) THEN
2029 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
2030 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
2031 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2032 +(c1f(k)*muuf(i ,j)+c2f(k))*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
2033 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2039 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
2040 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
2041 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2042 +(c1f(k)*muuf(i ,j)+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
2043 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2048 IF ( (config_flags%open_xe .or. specified) .and. its <= ide-2 .and. ide-2 <= ite ) THEN
2053 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.25*rdx/msfty(i,j))* &
2054 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(u(i+1,k,j)+u(i+1,k-1,j))*msfux(i+1,j)* &
2055 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2056 +(c1f(k)*muuf(i ,j)+c2f(k))*(u(i ,k,j)+u(i ,k-1,j))*msfux(i ,j)* &
2057 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2063 ph_tend(i,k,j)=ph_tend(i,k,j) - (0.5*rdx/msfty(i,j))* &
2064 ( (c1f(k)*muuf(i+1,j)+c2f(k))*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))*msfux(i+1,j)* &
2065 (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) &
2066 +(c1f(k)*muuf(i ,j)+c2f(k))*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))*msfux( i,j)* &
2067 (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) )
2072 END IF ! 6th order advection
2074 ! lateral open boundary conditions,
2075 ! start with north and south (y) boundaries
2082 IF ( (config_flags%open_ys) .and. jts == jds ) THEN
2089 vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) &
2090 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) )
2092 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*(c1f(k)*mut(i,j)+c2f(k))*( &
2093 +vl*(ph_old(i,k,j+1)-ph_old(i,k,j)))
2101 IF ( (config_flags%open_ye) .and. jte == jde ) THEN
2108 vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) &
2109 +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) )
2111 ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*(c1f(k)*mut(i,j)+c2f(k))*( &
2112 +vr*(ph_old(i,k,j)-ph_old(i,k,j-1)))
2118 ! now the east and west (y) boundaries
2125 IF ( (config_flags%open_xs) .and. its == ids ) THEN
2132 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2133 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2135 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*(c1f(k)*mut(i,j)+c2f(k))*( &
2136 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2141 ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) &
2142 +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) )
2144 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*(c1f(k)*mut(i,j)+c2f(k))*( &
2145 +ul*(ph_old(i+1,k,j)-ph_old(i,k,j)))
2152 IF ( (config_flags%open_xe) .and. ite == ide ) THEN
2159 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2160 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2162 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*(c1f(k)*mut(i,j)+c2f(k))*( &
2163 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2168 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) &
2169 +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) )
2171 ph_tend(i,k,j)=ph_tend(i,k,j)-(msftx(i,j)/msfty(i,j))*rdx*(c1f(k)*mut(i,j)+c2f(k))*( &
2172 +ur*(ph_old(i,k,j)-ph_old(i-1,k,j)))
2178 END SUBROUTINE rhs_ph
2181 !-------------------------------------------------------------------------------
2183 SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, &
2184 ph,alt,p,pb,al,php,cqu,cqv, &
2185 muu,muv,mu,c1h,c2h,fnm,fnp,rdnw,&
2186 cf1,cf2,cf3,cfn,cfn1, &
2187 rdx,rdy,msfux,msfuy,&
2188 msfvx,msfvy,msftx,msfty, &
2189 config_flags, non_hydrostatic, &
2191 ids, ide, jds, jde, kds, kde, &
2192 ims, ime, jms, jme, kms, kme, &
2193 its, ite, jts, jte, kts, kte )
2200 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2202 LOGICAL, INTENT (IN ) :: non_hydrostatic, top_lid
2204 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2205 ims, ime, jms, jme, kms, kme, &
2206 its, ite, jts, jte, kts, kte
2208 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: &
2219 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: &
2223 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, &
2228 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp
2230 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: c1h, c2h
2232 REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3, cfn, cfn1
2234 INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start
2235 REAL, DIMENSION( ims:ime, kms:kme ) :: dpn
2237 REAL, DIMENSION( kms:kme ) :: c1
2239 LOGICAL :: specified
2243 ! horizontal_pressure_gradient calculates the
2244 ! horizontal pressure gradient terms for the large-timestep tendency
2245 ! in the horizontal momentum equations (u,v).
2251 if(config_flags%specified .or. config_flags%nested) specified = .true.
2253 ! Notes on map scale factors:
2254 ! Calculates the pressure gradient terms in ADT eqns 44 and 45
2255 ! With upper rho -> 'mu', these are:
2256 ! Eqn 30: -mu*(mx/my)*(1/rho)*partial dp/dx
2257 ! Eqn 31: -mu*(my/mx)*(1/rho)*partial dp/dy
2259 ! As we are on nu, rather than height, surfaces:
2261 ! mu dp/dx = mu alpha partial dp'/dx + (nu mu partial dmubar/dx) alpha'
2262 ! + mu partial dphi'/dx + (partial dphi/dx)*(partial dp'/dnu - mu')
2264 ! mu dp/dy = mu alpha partial dp'/dy + (nu mu partial dmubar/dy) alpha'
2265 ! + mu partial dphi'/dy + (partial dphi/dy)*(partial dp'/dnu - mu')
2267 ! start with the north-south (y) pressure gradient
2274 IF ( (config_flags%open_ys .or. specified .or. &
2275 config_flags%nested .or. config_flags%polar ) .and. jts == jds ) j_start = jts+1
2276 IF ( (config_flags%open_ye .or. specified .or. &
2277 config_flags%nested .or. config_flags%polar ) .and. jte == jde ) jtf = jtf-1
2281 IF ( non_hydrostatic ) THEN
2286 dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) &
2287 +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) &
2288 +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) )
2293 dpn(i,kde) = .5*( cfn *(p(i,kde-1,j-1)+p(i,kde-1,j)) &
2294 +cfn1*(p(i,kde-2,j-1)+p(i,kde-2,j)) )
2300 dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) &
2301 +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) )
2305 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2306 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2309 ! Here are mu dp/dy terms 1-3
2310 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*(c1h(k)*muv(i,j)+c2h(k))*( &
2311 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2312 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2313 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2314 ! Here is mu dp/dy term 4
2315 dpy = dpy + (msfvy(i,j)/msfvx(i,j))*rdy*(php(i,k,j)-php(i,k,j-1))* &
2316 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*((c1(k)*mu(i,j-1))+(c1(k)*mu(i,j))))
2317 rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy
2323 ! ADT eqn 45: -mu*(my/mx)*(1/rho)*partial dp/dy
2324 ! [alt, al are 1/rho terms; muv, mu are NOT coupled]
2327 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2328 dpy = (msfvy(i,j)/msfvx(i,j))*.5*rdy*(c1h(k)*muv(i,j)+c2h(k))*( &
2329 (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) &
2330 +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) &
2331 +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) )
2332 rv_tend(i,k,j) = rv_tend(i,k,j)-dpy
2340 ! now the east-west (x) pressure gradient
2347 IF ( (config_flags%open_xs .or. specified .or. &
2348 config_flags%nested ) .and. its == ids ) i_start = its+1
2349 IF ( (config_flags%open_xe .or. specified .or. &
2350 config_flags%nested ) .and. ite == ide ) itf = itf-1
2351 IF ( config_flags%periodic_x ) i_start = its
2352 IF ( config_flags%periodic_x ) itf=ite
2356 IF ( non_hydrostatic ) THEN
2361 dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) &
2362 +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) &
2363 +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) )
2368 dpn(i,kde) = .5*( cfn *(p(i-1,kde-1,j)+p(i,kde-1,j)) &
2369 +cfn1*(p(i-1,kde-2,j)+p(i,kde-2,j)) )
2375 dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) &
2376 +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) )
2380 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2381 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2384 ! Here are mu dp/dy terms 1-3
2385 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*(c1h(k)*muu(i,j)+c2h(k))*( &
2386 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2387 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2388 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2389 ! Here is mu dp/dy term 4
2390 dpx = dpx + (msfux(i,j)/msfuy(i,j))*rdx*(php(i,k,j)-php(i-1,k,j))* &
2391 (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*((c1(k)*mu(i-1,j))+(c1(k)*mu(i,j))))
2392 ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx
2398 ! ADT eqn 44: -mu*(mx/my)*(1/rho)*partial dp/dx
2399 ! [alt, al are 1/rho terms; muu, mu are NOT coupled]
2402 ! Here are mu dp/dy terms 1-3; term 4 not needed if hydrostatic
2403 dpx = (msfux(i,j)/msfuy(i,j))*.5*rdx*(c1h(k)*muu(i,j)+c2h(k))*( &
2404 (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) &
2405 +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) &
2406 +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) )
2407 ru_tend(i,k,j) = ru_tend(i,k,j)-dpx
2415 END SUBROUTINE horizontal_pressure_gradient
2417 !-------------------------------------------------------------------------------
2419 SUBROUTINE pg_buoy_w( rw_tend, p, cqw, muf, mubf, &
2421 rdnw, rdn, g, msftx, msfty, &
2422 ids, ide, jds, jde, kds, kde, &
2423 ims, ime, jms, jme, kms, kme, &
2424 its, ite, jts, jte, kts, kte )
2430 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2431 ims, ime, jms, jme, kms, kme, &
2432 its, ite, jts, jte, kts, kte
2434 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p
2435 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw
2438 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2440 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mubf, muf, msftx, msfty
2442 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn
2444 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: c1f, c2f
2446 REAL, INTENT(IN ) :: g
2448 INTEGER :: itf, jtf, i, j, k
2454 ! pg_buoy_w calculates the
2455 ! vertical pressure gradient and buoyancy terms for the large-timestep
2456 ! tendency in the vertical momentum equation.
2460 ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T
2462 ! Map scale factor notes
2463 ! ADT eqn 46 RHS terms 6 and 7 (where 7 is "-rho g")
2464 ! Dividing by my, and using mu and nu (see Klemp et al. eqns 32, 40)
2465 ! term 6: +(g/my) partial dp'/dnu
2466 ! term 7: -(g/my) mu'
2468 ! For moisture-free atmosphere, cq1=1, cq2=0
2469 ! => (1./msft(i,j)) * g * [rdn(k)*{p(i,k,j)-p(i,k-1,j)}-(c1(k)*mu(i,j))]
2478 cq1 = 1./(1.+cqw(i,k-1,j))
2479 cq2 = cqw(i,k-1,j)*cq1
2480 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2481 cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) &
2482 -(c1f(k)*muf(i,j))-cq2*(c1f(k)*mubf(i,j)+c2f(k)) )
2487 cq1 = 1./(1.+cqw(i,k,j))
2488 cq2 = cqw(i,k,j)*cq1
2490 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msfty(i,j))*g*( &
2491 cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) &
2492 -(c1f(k)*muf(i,j))-cq2*(c1f(k)*mubf(i,j)+c2f(k)) )
2499 END SUBROUTINE pg_buoy_w
2501 !-------------------------------------------------------------------------------
2503 SUBROUTINE w_damp( rw_tend, max_vert_cfl,max_horiz_cfl, &
2504 u, v, ww, w, mut, c1f, c2f, rdnw, &
2505 rdx, rdy, msfux, msfuy, &
2508 ids, ide, jds, jde, kds, kde, &
2509 ims, ime, jms, jme, kms, kme, &
2510 its, ite, jts, jte, kts, kte )
2517 TYPE(grid_config_rec_type ) , INTENT(IN ) :: config_flags
2519 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2520 ims, ime, jms, jme, kms, kme, &
2521 its, ite, jts, jte, kts, kte
2523 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: u, v, ww, w
2525 REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend
2527 REAL, INTENT(OUT) :: max_vert_cfl
2528 REAL, INTENT(OUT) :: max_horiz_cfl
2531 REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut
2533 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw
2535 REAL, DIMENSION( kms:kme ), INTENT(IN ) :: c1f, c2f
2537 REAL, INTENT(IN) :: dt
2538 REAL, INTENT(IN) :: rdx, rdy
2539 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, msfuy
2540 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfvx, msfvy
2542 REAL :: vert_cfl, cf_n, cf_d, maxdub, maxdeta
2544 INTEGER :: itf, jtf, i, j, k, maxi, maxj, maxk
2545 CHARACTER*512 :: temp
2547 CHARACTER (LEN=256) :: time_str
2548 CHARACTER (LEN=256) :: grid_str
2551 REAL :: msfuxt , msfxffl
2554 REAL :: w_crit_cfl = 2.0
2555 REAL :: w_flag_cfl = 1.2
2556 LOGICAL :: print_flag = .true.
2558 INTEGER :: some1 ! Now have two catagories of CFL information, hence some1 & some2
2564 ! W_damp computes a damping term for the vertical velocity when the
2565 ! vertical Courant number is too large. This was found to be preferable to
2566 ! decreasing the timestep or increasing the diffusion in real-data applications
2567 ! that produced potentially-unstable large vertical velocities because of
2568 ! unphysically large heating rates coming from the cumulus parameterization
2569 ! schemes run at moderately high resolutions (dx ~ O(10) km).
2571 ! Additionally, w_damp returns the maximum cfl values due to vertical motion and
2572 ! horizontal motion. These values are returned via the max_vert_cfl and
2573 ! max_horiz_cfl variables. (Added by T. Hutchinson, WSI, 3/5/2007)
2577 ! W_damp modified to be more flexible with IEVA capability.
2578 ! When IEVA is on, the value of the W-CFL where damping is turned on can now be set
2579 ! in the namelist using "w_crit_cfl".
2581 ! Typical settings for non-IEVA use: w_crit_cfl = 1.0
2582 ! Typical settings for IEVA use: w_crit_cfl = 2.0
2584 ! I also commented out a lot of code put in 8-13 years ago for timing assuming that
2585 ! this is no longer needed....there is some pretty hacky stuff. HTH.
2587 ! (Added by L. Wicker, NSSL, 5/4/2020)
2601 w_crit_cfl = config_flags%w_crit_cfl
2602 ieva = config_flags%zadvect_implicit
2604 IF( ieva .gt. 0 ) THEN
2605 w_damp_on = w_crit_cfl
2610 IF( print_flag ) THEN
2611 write(wrf_err_message,*) '----------------------------------------'
2612 CALL wrf_debug( 0, wrf_err_message )
2613 WRITE(temp,*) 'W-DAMPING BEGINS AT W-COURANT NUMBER = ',w_damp_on
2614 CALL wrf_debug ( 0 , TRIM(temp) )
2615 write(wrf_err_message,*) '----------------------------------------'
2616 CALL wrf_debug( 0, wrf_err_message )
2617 print_flag = .false.
2620 IF(config_flags%polar ) then
2621 msfxffl = 1.0/COS(config_flags%fft_filter_lat*degrad)
2624 ! Routine has been reorganized to hopefully reduce redundant code while maintaining efficiency
2626 #ifdef OPTIMIZE_CFL_TEST
2627 ! 20121025, L. Meadows vector optimization does not include special case for Cassini
2629 IF(config_flags%polar ) then
2630 CALL wrf_error_fatal('module_big_step_utilities_em.F: -DOPTIMIZE_CFL_TEST option does not support global domains')
2639 vert_cfl = abs(ww(i,k,j)/(c1f(k)*mut(i,j)+c2f(k))*rdnw(k)*dt)
2641 # ifdef OPTIMIZE_CFL_TEST
2642 ! L. Meadows, Intel, MIC optimization, 20121025
2646 IF ( vert_cfl > max_vert_cfl ) THEN
2647 max_vert_cfl = vert_cfl
2650 IF(config_flags%polar ) THEN
2651 msfuxt = MIN(msfux(i,j), msfxffl)
2656 IF ( vert_cfl > max_vert_cfl ) THEN
2657 max_vert_cfl = vert_cfl ; maxi = i ; maxj = j ; maxk = k
2658 maxdub = w(i,k,j) ; maxdeta = -1./rdnw(k)
2662 horiz_cfl = max( abs(u(i,k,j) * rdx * msfuxt * dt), &
2663 abs(v(i,k,j) * rdy * msfvy(i,j) * dt) )
2665 IF (horiz_cfl > max_horiz_cfl) THEN
2666 max_horiz_cfl = horiz_cfl
2669 ! Dump out more information without affecting performance
2671 #ifndef OPTIMIZE_CFL_TEST
2672 ! This internal write is costly on newer Xeon processors because it breaks
2673 ! vectorization. (J. Michalakes for L. Meadows at Intel, 12/13/2012)
2675 IF (vert_cfl .gt. 2.) THEN
2679 WRITE(temp,FMT="(3(1x,i5,1x),'w-crit_cfl: ',f7.4,2x,'w-cfl: ',f7.4,2x,'dETA: ',f7.4)") &
2680 i,j,k,vert_cfl,w(i,k,j),-1./rdnw(k)
2681 CALL wrf_debug ( 100 , TRIM(temp) )
2687 IF ((vert_cfl .gt. w_damp_on) .and. (config_flags%w_damping == 1) ) THEN
2689 rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(vert_cfl-w_crit_cfl)*(c1f(k)*mut(i,j)+c2f(k))
2697 IF ( some1 .GT. 0 ) THEN
2698 CALL get_current_time_string( time_str )
2699 CALL get_current_grid_name( grid_str )
2700 WRITE(temp,*) some1, &
2701 ' points exceeded v_cfl = 2 in domain '//TRIM(grid_str)//' at time '//TRIM(time_str)//' hours'
2702 CALL wrf_debug ( 0 , TRIM(temp) )
2703 WRITE(temp,FMT="('Max W: ',3(1x,i5,1x),'W: ',f7.2,2x,'w-cfl: ',f7.2,2x,'dETA: ',f7.2)") &
2704 maxi, maxj, maxk, maxdub, max_vert_cfl, maxdeta
2705 CALL wrf_debug ( 0 , TRIM(temp) )
2706 ! WRITE(temp,FMT="('Max U/V: ',3(1x,i5,1x),'U: ',f7.2,2x,'u-cfl: ',f7.2,2x,'V: ',f7.2,2x,'v-cfl: ',f7.2)") &
2707 ! maxi, maxj, maxk, u(maxi,maxk,maxj), dt*u(maxi,maxk,maxj)*rdx, v(maxi,maxk,maxj), dt*v(maxi,maxk,maxj)*rdy
2708 ! CALL wrf_debug ( 0 , TRIM(temp) )
2711 END SUBROUTINE W_DAMP
2713 !-------------------------------------------------------------------------------
2715 SUBROUTINE horizontal_diffusion ( name, field, tendency, MUT, c1, c2, &
2717 msfux, msfuy, msfvx, msfvx_inv, &
2718 msfvy, msftx, msfty, &
2719 khdif, xkmhd, rdx, rdy, &
2720 ids, ide, jds, jde, kds, kde, &
2721 ims, ime, jms, jme, kms, kme, &
2722 its, ite, jts, jte, kts, kte )
2728 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2730 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2731 ims, ime, jms, jme, kms, kme, &
2732 its, ite, jts, jte, kts, kte
2734 CHARACTER(LEN=1) , INTENT(IN ) :: name
2736 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd
2738 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2740 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: MUT
2742 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2750 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
2752 REAL , INTENT(IN ) :: rdx, &
2758 INTEGER :: i, j, k, itf, jtf, ktf
2760 INTEGER :: i_start, i_end, j_start, j_end
2762 REAL :: mrdx, mkrdxm, mkrdxp, &
2763 mrdy, mkrdym, mkrdyp
2765 LOGICAL :: specified
2769 ! horizontal_diffusion computes the horizontal diffusion tendency
2770 ! on model horizontal coordinate surfaces.
2775 if(config_flags%specified .or. config_flags%nested) specified = .true.
2779 IF (name .EQ. 'u') THEN
2784 j_end = MIN(jte,jde-1)
2786 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2787 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
2788 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2789 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2790 IF ( config_flags%periodic_x ) i_start = its
2791 IF ( config_flags%periodic_x ) i_end = ite
2794 DO j = j_start, j_end
2796 DO i = i_start, i_end
2798 ! The interior is grad: (m_x*d/dx), the exterior is div: (m_x*m_y*d/dx(/m_y))
2799 ! setting up different averagings of m^2 partial d/dX and m^2 partial d/dY
2801 mkrdxm=(msftx(i-1,j)/msfty(i-1,j))*(c1(k)*MUT(i-1,j)+c2(k))*xkmhd(i-1,k,j)*rdx
2802 mkrdxp=(msftx(i,j)/msfty(i,j))*(c1(k)*MUT(i,j)+c2(k))*xkmhd(i,k,j)*rdx
2803 mrdx=msfux(i,j)*msfuy(i,j)*rdx
2804 mkrdym=( (msfuy(i,j)+msfuy(i,j-1))/(msfux(i,j)+msfux(i,j-1)) )* &
2805 0.25*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k))+(c1(k)*MUT(i-1,j-1)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k)))* &
2806 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy
2807 mkrdyp=( (msfuy(i,j)+msfuy(i,j+1))/(msfux(i,j)+msfux(i,j+1)) )* &
2808 0.25*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i-1,j+1)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k)))* &
2809 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy
2810 ! need to do four-corners (t) for diffusion coefficient as there are
2811 ! no values at u,v points
2812 ! msfuy - has to be y as part of d/dY
2813 ! has to be u as we're at a u point
2814 mrdy=msfux(i,j)*msfuy(i,j)*rdy
2816 ! correctly averaged version of rho~ * m^2 *
2817 ! [partial d/dX(partial du^/dX) + partial d/dY(partial du^/dY)]
2818 tendency(i,k,j)=tendency(i,k,j)+( &
2819 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2820 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2821 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2822 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2827 ELSE IF (name .EQ. 'v')THEN
2830 i_end = MIN(ite,ide-1)
2834 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2835 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2836 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2837 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
2838 IF ( config_flags%periodic_x ) i_start = its
2839 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2840 IF ( config_flags%polar ) j_start = MAX(jds+1,jts)
2841 IF ( config_flags%polar ) j_end = MIN(jde-1,jte)
2843 DO j = j_start, j_end
2845 DO i = i_start, i_end
2847 mkrdxm=( (msfvx(i,j)+msfvx(i-1,j))/(msfvy(i,j)+msfvy(i-1,j)) )* &
2848 0.25*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k))+(c1(k)*MUT(i-1,j-1)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k)))* &
2849 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx
2850 mkrdxp=( (msfvx(i,j)+msfvx(i+1,j))/(msfvy(i,j)+msfvy(i+1,j)) )* &
2851 0.25*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k))+(c1(k)*MUT(i+1,j-1)+c2(k))+(c1(k)*MUT(i+1,j)+c2(k)))* &
2852 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx
2853 mrdx=msfvx(i,j)*msfvy(i,j)*rdx
2854 mkrdym=(msfty(i,j-1)/msftx(i,j-1))*xkmhd(i,k,j-1)*rdy
2855 mkrdyp=(msfty(i,j)/msftx(i,j))*xkmhd(i,k,j)*rdy
2856 mrdy=msfvx(i,j)*msfvy(i,j)*rdy
2858 tendency(i,k,j)=tendency(i,k,j)+( &
2859 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2860 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2861 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2862 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2867 ELSE IF (name .EQ. 'w')THEN
2870 i_end = MIN(ite,ide-1)
2872 j_end = MIN(jte,jde-1)
2874 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2875 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2876 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2877 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2878 IF ( config_flags%periodic_x ) i_start = its
2879 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2881 DO j = j_start, j_end
2883 DO i = i_start, i_end
2885 mkrdxm=(msfux(i,j)/msfuy(i,j))* &
2886 0.25*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k))+(c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k)))* &
2887 0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx
2888 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))* &
2889 0.25*((c1(k)*MUT(i+1,j)+c2(k))+(c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i+1,j)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))* &
2890 0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx
2891 mrdx=msftx(i,j)*msfty(i,j)*rdx
2892 ! mkrdym=(msfvy(i,j)/msfvx(i,j))* &
2893 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))* &
2894 0.25*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k))+(c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k)))* &
2895 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy
2896 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))* &
2897 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))* &
2898 0.25*((c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))* &
2899 0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy
2900 mrdy=msftx(i,j)*msfty(i,j)*rdy
2902 tendency(i,k,j)=tendency(i,k,j)+( &
2903 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2904 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2905 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2906 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2915 i_end = MIN(ite,ide-1)
2917 j_end = MIN(jte,jde-1)
2919 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
2920 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
2921 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
2922 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
2923 IF ( config_flags%periodic_x ) i_start = its
2924 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
2926 DO j = j_start, j_end
2928 DO i = i_start, i_end
2930 mkrdxm=(msfux(i,j)/msfuy(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*0.5*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k)))*rdx
2931 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*0.5*((c1(k)*MUT(i+1,j)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))*rdx
2932 mrdx=msftx(i,j)*msfty(i,j)*rdx
2933 ! mkrdym=(msfvy(i,j)/msfvx(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k)))*rdy
2934 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k)))*rdy
2935 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*((c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))*rdy
2936 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*((c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))*rdy
2937 mrdy=msftx(i,j)*msfty(i,j)*rdy
2939 tendency(i,k,j)=tendency(i,k,j)+( &
2940 mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) &
2941 -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) &
2942 +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) &
2943 -mkrdym*(field(i,k,j )-field(i,k,j-1))))
2950 END SUBROUTINE horizontal_diffusion
2952 !-----------------------------------------------------------------------------------------
2954 SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, MUT, c1, c2, &
2955 config_flags, base_3d, &
2956 msfux, msfuy, msfvx, msfvx_inv, &
2957 msfvy, msftx, msfty, &
2958 khdif, xkmhd, rdx, rdy, &
2959 ids, ide, jds, jde, kds, kde, &
2960 ims, ime, jms, jme, kms, kme, &
2961 its, ite, jts, jte, kts, kte )
2967 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
2969 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
2970 ims, ime, jms, jme, kms, kme, &
2971 its, ite, jts, jte, kts, kte
2973 CHARACTER(LEN=1) , INTENT(IN ) :: name
2975 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, &
2979 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
2981 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: MUT
2983 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
2991 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
2993 REAL , INTENT(IN ) :: rdx, &
2999 INTEGER :: i, j, k, itf, jtf, ktf
3001 INTEGER :: i_start, i_end, j_start, j_end
3003 REAL :: mrdx, mkrdxm, mkrdxp, &
3004 mrdy, mkrdym, mkrdyp
3006 LOGICAL :: specified
3010 ! horizontal_diffusion_3dmp computes the horizontal diffusion tendency
3011 ! on model horizontal coordinate surfaces. This routine computes diffusion
3012 ! a perturbation scalar (field-base_3d).
3017 if(config_flags%specified .or. config_flags%nested) specified = .true.
3022 i_end = MIN(ite,ide-1)
3024 j_end = MIN(jte,jde-1)
3026 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3027 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite)
3028 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3029 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte)
3030 IF ( config_flags%periodic_x ) i_start = its
3031 IF ( config_flags%periodic_x ) i_end = MIN(ite,ide-1)
3033 DO j = j_start, j_end
3035 DO i = i_start, i_end
3037 mkrdxm=(msfux(i,j)/msfuy(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*0.5*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i-1,j)+c2(k)))*rdx
3038 mkrdxp=(msfux(i+1,j)/msfuy(i+1,j))*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*0.5*((c1(k)*MUT(i+1,j)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))*rdx
3039 mrdx=msftx(i,j)*msfty(i,j)*rdx
3040 ! mkrdym=(msfvy(i,j)/msfvx(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k)))*rdy
3041 ! mkrdyp=(msfvy(i,j+1)/msfvx(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*((c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))*rdy
3042 mkrdym=(msfvy(i,j)*msfvx_inv(i,j))*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*0.5*((c1(k)*MUT(i,j)+c2(k))+(c1(k)*MUT(i,j-1)+c2(k)))*rdy
3043 mkrdyp=(msfvy(i,j+1)*msfvx_inv(i,j+1))*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*0.5*((c1(k)*MUT(i,j+1)+c2(k))+(c1(k)*MUT(i,j)+c2(k)))*rdy
3044 mrdy=msftx(i,j)*msfty(i,j)*rdy
3046 tendency(i,k,j)=tendency(i,k,j)+( &
3047 mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) &
3048 -base_3d(i+1,k,j)+base_3d(i ,k,j) ) &
3049 -mkrdxm*( field(i ,k,j) -field(i-1,k,j) &
3050 -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) &
3051 +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) &
3052 -base_3d(i,k,j+1)+base_3d(i,k,j ) ) &
3053 -mkrdym*( field(i,k,j ) -field(i,k,j-1) &
3054 -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) &
3060 END SUBROUTINE horizontal_diffusion_3dmp
3062 !-----------------------------------------------------------------------------------------
3064 SUBROUTINE vertical_diffusion ( name, field, tendency, &
3065 config_flags, c1, c2, &
3066 alt, MUT, rdn, rdnw, kvdif, &
3067 ids, ide, jds, jde, kds, kde, &
3068 ims, ime, jms, jme, kms, kme, &
3069 its, ite, jts, jte, kts, kte )
3076 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3078 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3079 ims, ime, jms, jme, kms, kme, &
3080 its, ite, jts, jte, kts, kte
3082 CHARACTER(LEN=1) , INTENT(IN ) :: name
3084 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3085 INTENT(IN ) :: field, &
3088 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3090 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: MUT
3092 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw
3094 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
3096 REAL , INTENT(IN ) :: kvdif
3100 INTEGER :: i, j, k, itf, jtf, ktf
3101 INTEGER :: i_start, i_end, j_start, j_end
3103 REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz
3104 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3108 LOGICAL :: specified
3112 ! vertical_diffusion
3113 ! computes vertical diffusion tendency.
3118 if(config_flags%specified .or. config_flags%nested) specified = .true.
3122 IF (name .EQ. 'w')THEN
3126 i_end = MIN(ite,ide-1)
3128 j_end = MIN(jte,jde-1)
3130 j_loop_w : DO j = j_start, j_end
3133 DO i = i_start, i_end
3134 vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j))
3138 DO i = i_start, i_end
3143 DO i = i_start, i_end
3144 tendency(i,k,j)=tendency(i,k,j) &
3145 +rdn(k)*g*g/(c1(k)*MUT(i,j)+c2(k))/(0.5*(alt(i,k,j)+alt(i,k-1,j))) &
3146 *(vflux(i,k)-vflux(i,k-1))
3152 ELSE IF(name .EQ. 'm')THEN
3155 i_end = MIN(ite,ide-1)
3157 j_end = MIN(jte,jde-1)
3159 j_loop_s : DO j = j_start, j_end
3162 DO i = i_start, i_end
3163 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3164 *(field(i,k+1,j)-field(i,k,j))
3168 DO i = i_start, i_end
3169 vflux(i,0)=vflux(i,1)
3172 DO i = i_start, i_end
3177 DO i = i_start, i_end
3178 tendency(i,k,j)=tendency(i,k,j)+g*g/(c1(k)*MUT(i,j)+c2(k))/alt(i,k,j) &
3179 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3187 END SUBROUTINE vertical_diffusion
3190 !-------------------------------------------------------------------------------
3192 SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, &
3194 alt, MUT, rdn, rdnw, kvdif, &
3195 ids, ide, jds, jde, kds, kde, &
3196 ims, ime, jms, jme, kms, kme, &
3197 its, ite, jts, jte, kts, kte )
3204 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3206 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3207 ims, ime, jms, jme, kms, kme, &
3208 its, ite, jts, jte, kts, kte
3210 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3211 INTENT(IN ) :: field, &
3214 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3216 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: MUT
3218 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3222 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
3224 REAL , INTENT(IN ) :: kvdif
3228 INTEGER :: i, j, k, itf, jtf, ktf
3229 INTEGER :: i_start, i_end, j_start, j_end
3231 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3235 LOGICAL :: specified
3239 ! vertical_diffusion_mp
3240 ! computes vertical diffusion tendency of a perturbation variable
3241 ! (field-base). Note that base as a 1D (k) field.
3246 if(config_flags%specified .or. config_flags%nested) specified = .true.
3251 i_end = MIN(ite,ide-1)
3253 j_end = MIN(jte,jde-1)
3255 j_loop_s : DO j = j_start, j_end
3258 DO i = i_start, i_end
3259 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3260 *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k))
3264 DO i = i_start, i_end
3265 vflux(i,0)=vflux(i,1)
3268 DO i = i_start, i_end
3273 DO i = i_start, i_end
3274 tendency(i,k,j)=tendency(i,k,j)+g*g/(c1(k)*MUT(i,j)+c2(k))/alt(i,k,j) &
3275 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3281 END SUBROUTINE vertical_diffusion_mp
3284 !-------------------------------------------------------------------------------
3286 SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, &
3288 alt, MUT, rdn, rdnw, kvdif, &
3289 ids, ide, jds, jde, kds, kde, &
3290 ims, ime, jms, jme, kms, kme, &
3291 its, ite, jts, jte, kts, kte )
3298 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3300 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3301 ims, ime, jms, jme, kms, kme, &
3302 its, ite, jts, jte, kts, kte
3304 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3305 INTENT(IN ) :: field, &
3309 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3311 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: MUT
3313 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, &
3316 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1, c2
3318 REAL , INTENT(IN ) :: kvdif
3322 INTEGER :: i, j, k, itf, jtf, ktf
3323 INTEGER :: i_start, i_end, j_start, j_end
3325 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3329 LOGICAL :: specified
3333 ! vertical_diffusion_3dmp
3334 ! computes vertical diffusion tendency of a perturbation variable
3340 if(config_flags%specified .or. config_flags%nested) specified = .true.
3345 i_end = MIN(ite,ide-1)
3347 j_end = MIN(jte,jde-1)
3349 j_loop_s : DO j = j_start, j_end
3352 DO i = i_start, i_end
3353 vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) &
3354 *( field(i,k+1,j) -field(i,k,j) &
3355 -base_3d(i,k+1,j)+base_3d(i,k,j) )
3359 DO i = i_start, i_end
3360 vflux(i,0)=vflux(i,1)
3363 DO i = i_start, i_end
3368 DO i = i_start, i_end
3369 tendency(i,k,j)=tendency(i,k,j)+g*g/(c1(k)*MUT(i,j)+c2(k))/alt(i,k,j) &
3370 *rdnw(k)*(vflux(i,k)-vflux(i,k-1))
3376 END SUBROUTINE vertical_diffusion_3dmp
3379 !-------------------------------------------------------------------------------
3382 SUBROUTINE vertical_diffusion_u ( field, tendency, &
3383 config_flags, u_base, c1h,c2h,&
3384 alt, muu, rdn, rdnw, kvdif, &
3385 ids, ide, jds, jde, kds, kde, &
3386 ims, ime, jms, jme, kms, kme, &
3387 its, ite, jts, jte, kts, kte )
3394 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3396 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3397 ims, ime, jms, jme, kms, kme, &
3398 its, ite, jts, jte, kts, kte
3400 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3401 INTENT(IN ) :: field, &
3404 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3406 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu
3408 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base
3410 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h
3412 REAL , INTENT(IN ) :: kvdif
3416 INTEGER :: i, j, k, itf, jtf, ktf
3417 INTEGER :: i_start, i_end, j_start, j_end
3419 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3423 LOGICAL :: specified
3427 ! vertical_diffusion_u computes vertical diffusion tendency for
3428 ! the u momentum equation. This routine assumes a constant eddy
3434 if(config_flags%specified .or. config_flags%nested) specified = .true.
3441 j_end = MIN(jte,jde-1)
3443 IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its)
3444 IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite)
3445 IF ( config_flags%periodic_x ) i_start = its
3446 IF ( config_flags%periodic_x ) i_end = ite
3449 j_loop_u : DO j = j_start, j_end
3452 DO i = i_start, i_end
3453 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) &
3456 +alt(i-1,k+1,j) ) ) &
3457 *(field(i,k+1,j)-field(i,k,j) &
3458 -u_base(k+1) +u_base(k) )
3462 DO i = i_start, i_end
3463 vflux(i,0)=vflux(i,1)
3466 DO i = i_start, i_end
3471 DO i = i_start, i_end
3472 tendency(i,k,j)=tendency(i,k,j)+ &
3473 g*g*rdnw(k)/(c1h(k)*muu(i,j)+c2h(k))/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* &
3474 (vflux(i,k)-vflux(i,k-1))
3480 END SUBROUTINE vertical_diffusion_u
3482 !-------------------------------------------------------------------------------
3485 SUBROUTINE vertical_diffusion_v ( field, tendency, &
3486 config_flags, v_base, c1h,c2h,&
3487 alt, muv, rdn, rdnw, kvdif, &
3488 ids, ide, jds, jde, kds, kde, &
3489 ims, ime, jms, jme, kms, kme, &
3490 its, ite, jts, jte, kts, kte )
3497 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
3499 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3500 ims, ime, jms, jme, kms, kme, &
3501 its, ite, jts, jte, kts, kte
3503 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
3504 INTENT(IN ) :: field, &
3506 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base
3508 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h
3510 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
3512 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv
3514 REAL , INTENT(IN ) :: kvdif
3518 INTEGER :: i, j, k, itf, jtf, ktf, jm1
3519 INTEGER :: i_start, i_end, j_start, j_end
3521 REAL , DIMENSION(its:ite, 0:kte+1) :: vflux
3525 LOGICAL :: specified
3529 ! vertical_diffusion_v computes vertical diffusion tendency for
3530 ! the v momentum equation. This routine assumes a constant eddy
3536 if(config_flags%specified .or. config_flags%nested) specified = .true.
3541 i_end = MIN(ite,ide-1)
3543 j_end = MIN(jte,jde-1)
3545 IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts)
3546 IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte)
3548 j_loop_v : DO j = j_start, j_end
3553 DO i = i_start, i_end
3554 vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) &
3557 +alt(i,k+1,jm1) ) ) &
3558 *(field(i,k+1,j)-field(i,k,j) &
3559 -v_base(k+1) +v_base(k) )
3563 DO i = i_start, i_end
3564 vflux(i,0)=vflux(i,1)
3567 DO i = i_start, i_end
3572 DO i = i_start, i_end
3573 tendency(i,k,j)=tendency(i,k,j)+ &
3574 g*g*rdnw(k)/(c1h(k)*muv(i,j)+c2h(k))/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* &
3575 (vflux(i,k)-vflux(i,k-1))
3581 END SUBROUTINE vertical_diffusion_v
3583 !*************** end new mass coordinate routines
3585 !-------------------------------------------------------------------------------
3587 SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, &
3588 ids, ide, jds, jde, kds, kde, &
3589 ims, ime, jms, jme, kms, kme, &
3590 its, ite, jts, jte, kts, kte )
3596 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3597 ims, ime, jms, jme, kms, kme, &
3598 its, ite, jts, jte, kts, kte
3600 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, &
3603 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield
3607 INTEGER :: i, j, k, itf, jtf, ktf, i_start, j_start
3612 ! calculates full 3D field from pertubation and base field.
3613 ! The input fields (mu base and perturbation mu) are communicated prior
3614 ! to this call. That fills the total field into the halo region. That
3615 ! extra row/column is only used by the IEVA scheme. Having the extra
3616 ! row/column filled with valid values is not a problem for the rest of
3631 rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j)
3636 END SUBROUTINE calculate_full
3638 !------------------------------------------------------------------------------
3640 SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, &
3642 msftx, msfty, msfux, msfuy, &
3644 f, e, sina, cosa, fzm, fzp, &
3645 ids, ide, jds, jde, kds, kde, &
3646 ims, ime, jms, jme, kms, kme, &
3647 its, ite, jts, jte, kts, kte )
3653 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3655 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3656 ims, ime, jms, jme, kms, kme, &
3657 its, ite, jts, jte, kts, kte
3659 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3662 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, &
3666 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3673 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3678 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3683 INTEGER :: i, j , k, ktf
3684 INTEGER :: i_start, i_end, j_start, j_end
3686 LOGICAL :: specified
3690 ! coriolis calculates the large timestep tendency terms in the
3691 ! u, v, and w momentum equations arise from the coriolis force.
3696 if(config_flags%specified .or. config_flags%nested) specified = .true.
3700 ! coriolis for u-momentum equation
3702 ! Notes on map scale factor
3703 ! cosa, sina are related to rotating the coordinate frame if desired
3704 ! generally sina=0, cosa=1
3705 ! ADT eqn 44, RHS terms 6 and 7: -2 mu w omega cos(lat)/my
3706 ! + 2 mu v omega sin(lat)/my
3707 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3708 ! => terms are: -e mu w / my + f mu v / my
3709 ! rv = mu v / mx ; rw = mu w / my
3710 ! => terms are: -e rw + f rv *mx / my
3714 IF ( config_flags%open_xs .or. specified .or. &
3715 config_flags%nested) i_start = MAX(ids+1,its)
3716 IF ( config_flags%open_xe .or. specified .or. &
3717 config_flags%nested) i_end = MIN(ide-1,ite)
3718 IF ( config_flags%periodic_x ) i_start = its
3719 IF ( config_flags%periodic_x ) i_end = ite
3721 DO j = jts, MIN(jte,jde-1)
3724 DO i = i_start, i_end
3726 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)) &
3727 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
3728 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
3729 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
3734 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3735 ! IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
3739 ! ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
3740 ! *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
3741 ! - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
3742 ! *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
3748 ! IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
3752 ! 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)) &
3753 ! *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
3754 ! - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
3755 ! *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
3763 ! coriolis term for v-momentum equation
3765 ! Notes on map scale factors
3766 ! ADT eqn 45, RHS terms 6 and 6b [0 for sina=0]: -2 mu u omega sin(lat)/mx + ?
3767 ! Define f=2 omega sin(lat), e=2 omega cos(lat)
3768 ! => terms are: -f mu u / mx
3769 ! ru = mu u / my ; rw = mu w / my
3770 ! => terms are: -f ru *my / mx + ?
3775 IF ( config_flags%open_ys .or. specified .or. &
3776 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
3777 IF ( config_flags%open_ye .or. specified .or. &
3778 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
3780 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3781 ! IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
3784 ! DO i=its,MIN(ide-1,ite)
3786 ! rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
3787 ! *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
3788 ! + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
3789 ! *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
3798 DO i=its,MIN(ide-1,ite)
3800 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)) &
3801 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
3802 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
3803 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
3810 ! boundary loops for coriolis not needed for open bdy (commented out 20100611 JD)
3811 ! IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
3814 ! DO i=its,MIN(ide-1,ite)
3816 ! 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)) &
3817 ! *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
3818 ! + (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)) &
3819 ! *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
3826 ! coriolis term for w-mometum
3828 ! Notes on map scale factors
3829 ! ADT eqn 46/my, RHS terms 5 and 5b [0 for sina=0]: 2 mu u omega cos(lat)/my +?
3830 ! Define e=2 omega cos(lat)
3831 ! => terms are: e mu u / my + ???
3832 ! ru = mu u / my ; ru = mu v / mx
3833 ! => terms are: e ru + ???
3835 DO j=jts,MIN(jte, jde-1)
3837 DO i=its,MIN(ite, ide-1)
3839 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
3840 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
3841 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
3842 -(msftx(i,j)/msfty(i,j))* &
3843 sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
3844 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
3850 END SUBROUTINE coriolis
3852 !------------------------------------------------------------------------------
3854 SUBROUTINE perturbation_coriolis ( ru_in, rv_in, rw, ru_tend, rv_tend, rw_tend, &
3856 u_base, v_base, z_base, &
3857 muu, muv, c1h, c2h, phb, ph, &
3858 msftx, msfty, msfux, msfuy, msfvx, msfvy, &
3859 f, e, sina, cosa, fzm, fzp, &
3860 ids, ide, jds, jde, kds, kde, &
3861 ims, ime, jms, jme, kms, kme, &
3862 its, ite, jts, jte, kts, kte )
3868 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
3870 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
3871 ims, ime, jms, jme, kms, kme, &
3872 its, ite, jts, jte, kts, kte
3874 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, &
3877 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru_in, &
3884 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
3891 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, &
3896 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu, &
3900 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
3903 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base, &
3907 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h
3911 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) :: ru, &
3914 REAL :: z_at_u, z_at_v, wkp1, wk, wkm1
3918 INTEGER :: i, j , k, ktf
3919 INTEGER :: i_start, i_end, j_start, j_end
3921 LOGICAL :: specified
3925 ! perturbation_coriolis calculates the large timestep tendency terms in the
3926 ! u, v, and w momentum equations arise from the coriolis force. This version
3927 ! subtracts off the horizontal velocities from the initial sounding when
3928 ! computing the forcing terms, hence "perturbation" coriolis.
3933 if(config_flags%specified .or. config_flags%nested) specified = .true.
3937 ! coriolis for u-momentum equation
3941 IF ( config_flags%open_xs .or. specified .or. &
3942 config_flags%nested) i_start = MAX(ids+1,its)
3943 IF ( config_flags%open_xe .or. specified .or. &
3944 config_flags%nested) i_end = MIN(ide-1,ite)
3945 IF ( config_flags%periodic_x ) i_start = its
3946 IF ( config_flags%periodic_x ) i_end = ite
3948 ! compute perturbation mu*v for use in u momentum equation
3950 DO j = jts, MIN(jte,jde-1)+1
3952 DO i = i_start-1, i_end
3953 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3954 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3955 +ph(i,k,j )+ph(i,k+1,j ) &
3956 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3957 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3958 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3960 rv(i,k,j) = rv_in(i,k,j) - (c1h(k)*muv(i,j)+c2h(k))*( &
3969 ! pick up top and bottom v
3971 DO j = jts, MIN(jte,jde-1)+1
3972 DO i = i_start-1, i_end
3975 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3976 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3977 +ph(i,k,j )+ph(i,k+1,j ) &
3978 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3979 wkp1 = min(1.,max(0.,z_at_v-z_base(k))/(z_base(k+1)-z_base(k)))
3981 rv(i,k,j) = rv_in(i,k,j) - (c1h(k)*muv(i,j)+c2h(k))*( &
3986 z_at_v = 0.25*( phb(i,k,j )+phb(i,k+1,j ) &
3987 +phb(i,k,j-1)+phb(i,k+1,j-1) &
3988 +ph(i,k,j )+ph(i,k+1,j ) &
3989 +ph(i,k,j-1)+ph(i,k+1,j-1))/g
3990 wkm1 = min(1.,max(0.,z_base(k)-z_at_v)/(z_base(k)-z_base(k-1)))
3992 rv(i,k,j) = rv_in(i,k,j) - (c1h(k)*muv(i,j)+c2h(k))*( &
3999 ! compute coriolis forcing for u
4001 ! Map scale factors: see comments above for Coriolis
4003 DO j = jts, MIN(jte,jde-1)
4006 DO i = i_start, i_end
4007 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)) &
4008 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4009 - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) &
4010 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4014 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4015 IF ( (config_flags%open_xs) .and. (its == ids) ) THEN
4019 ru_tend(its,k,j)=ru_tend(its,k,j) + (msfux(its,j)/msfuy(its,j))*0.5*(f(its,j)+f(its,j)) &
4020 *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) &
4021 - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) &
4022 *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j))
4028 IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN
4032 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)) &
4033 *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) &
4034 - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) &
4035 *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j))
4043 ! coriolis term for v-momentum equation
4044 ! Map scale factors: see comments above for Coriolis
4049 IF ( config_flags%open_ys .or. specified .or. &
4050 config_flags%nested .or. config_flags%polar) j_start = MAX(jds+1,jts)
4051 IF ( config_flags%open_ye .or. specified .or. &
4052 config_flags%nested .or. config_flags%polar) j_end = MIN(jde-1,jte)
4054 ! compute perturbation mu*u for use in v momentum equation
4056 DO j = j_start-1,j_end
4058 DO i = its, MIN(ite,ide-1)+1
4059 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
4060 +phb(i-1,k,j)+phb(i-1,k+1,j) &
4061 +ph(i ,k,j)+ph(i ,k+1,j) &
4062 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
4063 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
4064 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
4066 ru(i,k,j) = ru_in(i,k,j) - (c1h(k)*muu(i,j)+c2h(k))*( &
4074 ! pick up top and bottom u
4076 DO j = j_start-1,j_end
4077 DO i = its, MIN(ite,ide-1)+1
4080 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
4081 +phb(i-1,k,j)+phb(i-1,k+1,j) &
4082 +ph(i ,k,j)+ph(i ,k+1,j) &
4083 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
4084 wkp1 = min(1.,max(0.,z_at_u-z_base(k))/(z_base(k+1)-z_base(k)))
4086 ru(i,k,j) = ru_in(i,k,j) - (c1h(k)*muu(i,j)+c2h(k))*( &
4092 z_at_u = 0.25*( phb(i ,k,j)+phb(i ,k+1,j) &
4093 +phb(i-1,k,j)+phb(i-1,k+1,j) &
4094 +ph(i ,k,j)+ph(i ,k+1,j) &
4095 +ph(i-1,k,j)+ph(i-1,k+1,j))/g
4096 wkm1 = min(1.,max(0.,z_base(k)-z_at_u)/(z_base(k)-z_base(k-1)))
4098 ru(i,k,j) = ru_in(i,k,j) - (c1h(k)*muu(i,j)+c2h(k))*( &
4105 ! compute coriolis forcing for v momentum equation
4106 ! Map scale factors: see comments above for Coriolis
4108 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4109 IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN
4112 DO i=its,MIN(ide-1,ite)
4114 rv_tend(i,k,jts)=rv_tend(i,k,jts) - (msfvy(i,jts)/msfvx(i,jts))*0.5*(f(i,jts)+f(i,jts)) &
4115 *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) &
4116 + (msfvy(i,jts)/msfvx(i,jts))*0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) &
4117 *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts))
4126 DO i=its,MIN(ide-1,ite)
4128 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)) &
4129 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4130 + (msfvy(i,j)/msfvx(i,j))*0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) &
4131 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4138 ! boundary loops for perturbation coriolis is needed for open bdy (20110225 JD)
4139 IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN
4142 DO i=its,MIN(ide-1,ite)
4144 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)) &
4145 *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) &
4146 + (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)) &
4147 *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1))
4154 ! coriolis term for w-mometum
4155 ! Map scale factors: see comments above for Coriolis
4157 DO j=jts,MIN(jte, jde-1)
4159 DO i=its,MIN(ite, ide-1)
4161 rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* &
4162 (cosa(i,j)*0.5*(fzm(k)*(ru(i,k,j)+ru(i+1,k,j)) &
4163 +fzp(k)*(ru(i,k-1,j)+ru(i+1,k-1,j))) &
4164 -(msftx(i,j)/msfty(i,j))*sina(i,j)*0.5*(fzm(k)*(rv(i,k,j)+rv(i,k,j+1)) &
4165 +fzp(k)*(rv(i,k-1,j)+rv(i,k-1,j+1))))
4171 END SUBROUTINE perturbation_coriolis
4173 !------------------------------------------------------------------------------
4175 SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, &
4177 msfux, msfuy, msfvx, msfvy, msftx, msfty, &
4178 xlat, fzm, fzp, rdx, rdy, &
4179 ids, ide, jds, jde, kds, kde, &
4180 ims, ime, jms, jme, kms, kme, &
4181 its, ite, jts, jte, kts, kte )
4188 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4190 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4191 ims, ime, jms, jme, kms, kme, &
4192 its, ite, jts, jte, kts, kte
4194 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4195 INTENT(INOUT) :: ru_tend, &
4199 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4200 INTENT(IN ) :: ru, &
4207 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfux, &
4215 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4218 REAL , INTENT(IN ) :: rdx, &
4223 ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp
4224 INTEGER :: i, j, k, itf, jtf, ktf
4225 INTEGER :: i_start, i_end, j_start, j_end
4226 ! INTEGER :: irmin, irmax, jrmin, jrmax
4228 REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm
4230 LOGICAL :: specified
4234 ! curvature calculates the large timestep tendency terms in the
4235 ! u, v, and w momentum equations arise from the curvature terms.
4240 if(config_flags%specified .or. config_flags%nested) specified = .true.
4250 ! IF ( config_flags%open_xs ) irmin = ids
4251 ! IF ( config_flags%open_xe ) irmax = ide-1
4252 ! IF ( config_flags%open_ys ) jrmin = jds
4253 ! IF ( config_flags%open_ye ) jrmax = jde-1
4255 ! Define v cross grad m at scalar points - vxgm(i,j)
4262 IF ( ( config_flags%open_xs .or. specified .or. &
4263 config_flags%nested) .and. (its == ids) ) i_start = its
4264 IF ( ( config_flags%open_xe .or. specified .or. &
4265 config_flags%nested) .and. (ite == ide) ) i_end = ite-1
4266 IF ( ( config_flags%open_ys .or. specified .or. &
4267 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) j_start = jts
4268 IF ( ( config_flags%open_ye .or. specified .or. &
4269 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) j_end = jte-1
4270 IF ( config_flags%periodic_x ) i_start = its-1
4271 IF ( config_flags%periodic_x ) i_end = ite
4276 ! Map scale factor notes:
4277 ! msf...y is constant everywhere for cylindrical map projection
4278 ! msf...x varies with y only
4279 ! But we know that this is not = 0 for cylindrical,
4280 ! therefore use msfvX in 1st line
4281 ! which => by symmetry use msfuY in 2nd line - ???
4282 vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfvx(i,j+1)-msfvx(i,j))*rdy - &
4283 0.5*(v(i,k,j)+v(i,k,j+1))*(msfuy(i+1,j)-msfuy(i,j))*rdx
4288 ! Pick up the boundary rows for open (radiation) lateral b.c.
4289 ! Rather crude at present, we are assuming there is no
4290 ! variation in this term at the boundary.
4292 IF ( ( config_flags%open_xs .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4293 config_flags%nested) .and. (its == ids) ) THEN
4297 vxgm(its-1,k,j) = vxgm(its,k,j)
4303 IF ( ( config_flags%open_xe .or. (specified .AND. .NOT. config_flags%periodic_x) .or. &
4304 config_flags%nested) .and. (ite == ide) ) THEN
4308 vxgm(ite,k,j) = vxgm(ite-1,k,j)
4314 ! Polar boundary condition:
4315 ! The following change is needed in case one tries using the vxgm route with
4316 ! polar B.C.'s in the future, but not needed if 'tan' used
4317 IF ( ( config_flags%open_ys .or. specified .or. &
4318 config_flags%nested .or. config_flags%polar) .and. (jts == jds) ) THEN
4322 vxgm(i,k,jts-1) = vxgm(i,k,jts)
4328 ! Polar boundary condition:
4329 ! The following change is needed in case one tries using the vxgm route with
4330 ! polar B.C.'s in the future, but not needed if 'tan' used
4331 IF ( ( config_flags%open_ye .or. specified .or. &
4332 config_flags%nested .or. config_flags%polar) .and. (jte == jde) ) THEN
4336 vxgm(i,k,jte) = vxgm(i,k,jte-1)
4342 ! curvature term for u momentum eqn.
4344 ! Map scale factor notes:
4345 ! ADT eqn 44, RHS terms 4 and 5, in cylindrical: mu u v tan(lat)/(a my)
4347 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4349 ! (mx/my)*u rv tan(lat) / a - u rw / a = (u/a)*[(mx/my) rv tan(lat) - rw]
4350 ! ru v tan(lat) / a - u rw / a
4351 ! xlat defined with end points half grid space from pole,
4352 ! hence are on u latitude points
4355 IF ( config_flags%open_xs .or. specified .or. &
4356 config_flags%nested) i_start = MAX ( ids+1 , its )
4357 IF ( config_flags%open_xe .or. specified .or. &
4358 config_flags%nested) i_end = MIN ( ide-1 , ite )
4359 IF ( config_flags%periodic_x ) i_start = its
4360 IF ( config_flags%periodic_x ) i_end = ite
4362 ! Polar boundary condition
4363 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4365 DO j=jts,MIN(jde-1,jte)
4369 ru_tend(i,k,j)=ru_tend(i,k,j) + u(i,k,j)*reradius* ( &
4370 (msfux(i,j)/msfuy(i,j))*0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+ &
4371 rv(i-1,k,j)+rv(i,k,j))*tan(xlat(i,j)*degrad) &
4372 - 0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) )
4380 DO j=jts,MIN(jde-1,jte)
4384 ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) &
4385 *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) &
4386 - u(i,k,j)*reradius &
4387 *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j))
4395 ! curvature term for v momentum eqn.
4397 ! Map scale factor notes
4398 ! ADT eqn 45, RHS terms 4 and 5, in cylindrical: - mu u*u tan(lat)/(a mx)
4400 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4402 ! - (my/mx)*u ru tan(lat) / a - (my/mx)*v rw / a
4403 ! = - [my/(mx*a)]*[u ru tan(lat) + v rw]
4404 ! - (1/a)*[(my/mx)*u ru tan(lat) + w rv]
4405 ! xlat defined with end points half grid space from pole, hence are on
4406 ! u latitude points => av here
4408 ! in original wrf, there was a sign error for the rw contribution
4411 IF ( config_flags%open_ys .or. specified .or. &
4412 config_flags%nested .or. config_flags%polar) j_start = MAX ( jds+1 , jts )
4413 IF ( config_flags%open_ye .or. specified .or. &
4414 config_flags%nested .or. config_flags%polar) j_end = MIN ( jde-1 , jte )
4416 IF ((config_flags%map_proj == 6) .OR. (config_flags%polar)) THEN
4420 DO i=its,MIN(ite,ide-1)
4421 rv_tend(i,k,j)=rv_tend(i,k,j) - (msfvy(i,j)/msfvx(i,j))*reradius* ( &
4422 0.25*(u(i,k,j)+u(i+1,k,j)+u(i,k,j-1)+u(i+1,k,j-1))* &
4423 tan((xlat(i,j)+xlat(i,j-1))*0.5*degrad)* &
4424 0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4425 + v(i,k,j)*0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+ &
4426 rw(i,k+1,j)+rw(i,k,j)) )
4435 DO i=its,MIN(ite,ide-1)
4437 rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) &
4438 *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) &
4439 - (msfvy(i,j)/msfvx(i,j))*v(i,k,j)*reradius &
4440 *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j))
4448 ! curvature term for vertical momentum eqn.
4450 ! Notes on map scale factors:
4451 ! ADT eqn 46, RHS term 4: [mu/(a my)]*[u*u + v*v]
4452 ! ru = mu u / my ; rw = mu w / my ; rv = mu v / mx
4453 ! terms are: u ru / a + (mx/my)v rv / a
4455 DO j=jts,MIN(jte,jde-1)
4457 DO i=its,MIN(ite,ide-1)
4459 rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* &
4460 (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))) &
4461 *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))) &
4462 +(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))) &
4463 *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))))
4469 END SUBROUTINE curvature
4472 DANGER - this is a bad routine to have laying around - someone could use it
4473 !------------------------------------------------------------------------------
4475 SUBROUTINE decouple ( rr, rfield, field, name, config_flags, &
4477 ids, ide, jds, jde, kds, kde, &
4478 ims, ime, jms, jme, kms, kme, &
4479 its, ite, jts, jte, kts, kte )
4485 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4487 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4488 ims, ime, jms, jme, kms, kme, &
4489 its, ite, jts, jte, kts, kte
4491 CHARACTER(LEN=1) , INTENT(IN ) :: name
4493 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield
4495 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr
4497 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field
4499 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp
4503 INTEGER :: i, j, k, itf, jtf, ktf
4507 ! decouple decouples a variable from the column dry-air mass.
4513 IF (name .EQ. 'u')THEN
4520 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j)))
4525 ELSE IF (name .EQ. 'v')THEN
4532 field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1)))
4537 ELSE IF (name .EQ. 'w')THEN
4543 field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j))
4557 ! For theta we will decouple tb and tp and add them to give t afterwards
4561 field(i,k,j)=rfield(i,k,j)/rr(i,k,j)
4568 END SUBROUTINE decouple
4571 !-------------------------------------------------------------------------------
4573 SUBROUTINE zero_tend ( tendency, &
4574 ids, ide, jds, jde, kds, kde, &
4575 ims, ime, jms, jme, kms, kme, &
4576 its, ite, jts, jte, kts, kte )
4583 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4584 ims, ime, jms, jme, kms, kme, &
4585 its, ite, jts, jte, kts, kte
4587 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency
4591 INTEGER :: i, j, k, itf, jtf, ktf
4595 ! zero_tend sets the input tendency array to zero.
4602 tendency(i,k,j) = 0.
4607 END SUBROUTINE zero_tend
4609 !-------------------------------------------------------------------------------
4611 SUBROUTINE zero_tend2d( tendency, &
4612 ids, ide, jds, jde, kds, kde, &
4613 ims, ime, jms, jme, kms, kme, &
4614 its, ite, jts, jte, kts, kte )
4621 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4622 ims, ime, jms, jme, kms, kme, &
4623 its, ite, jts, jte, kts, kte
4625 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(INOUT) :: tendency
4629 INTEGER :: i, j, k, itf, jtf, ktf
4633 ! zero_tend sets the input tendency array to zero.
4643 END SUBROUTINE zero_tend2d
4645 !-------------------------------------------------------------------------------
4647 ! Sets the an array on the polar v point(s) to zero
4648 SUBROUTINE zero_pole ( field, &
4649 ids, ide, jds, jde, kds, kde, &
4650 ims, ime, jms, jme, kms, kme, &
4651 its, ite, jts, jte, kts, kte )
4658 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4659 ims, ime, jms, jme, kms, kme, &
4660 its, ite, jts, jte, kts, kte
4662 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4668 IF (jts == jds) THEN
4675 IF (jte == jde) THEN
4683 END SUBROUTINE zero_pole
4685 !-------------------------------------------------------------------------------
4686 ! Sets the an array on the polar v point(s)
4687 SUBROUTINE pole_point_bc ( field, &
4688 ids, ide, jds, jde, kds, kde, &
4689 ims, ime, jms, jme, kms, kme, &
4690 its, ite, jts, jte, kts, kte )
4697 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4698 ims, ime, jms, jme, kms, kme, &
4699 its, ite, jts, jte, kts, kte
4701 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: field
4707 IF (jts == jds) THEN
4710 ! field(i,k,jts) = 2*field(i,k,jts+1) - field(i,k,jts+2)
4711 field(i,k,jts) = field(i,k,jts+1)
4715 IF (jte == jde) THEN
4718 ! field(i,k,jte) = 2*field(i,k,jte-1) - field(i,k,jte-2)
4719 field(i,k,jte) = field(i,k,jte-1)
4724 END SUBROUTINE pole_point_bc
4726 !======================================================================
4727 ! physics prep routines
4728 !======================================================================
4730 SUBROUTINE phy_prep ( config_flags, & ! input
4732 c1h, c2h, c1f, c2f, &
4733 u, v, p, pb, alt, ph, & ! input
4734 phb, t, moist, n_moist, & ! input
4735 rho, th_phy, th_phy_m_t0, & ! output
4736 p_phy , pi_phy , & ! output
4737 u_phy, v_phy, p8w, t_phy, t8w, & ! output
4738 z, z_at_w, dz8w, & ! output
4739 p_hyd, p_hyd_w, dnw, & ! output
4740 fzm, fzp, znw, p_top, & ! params
4741 ids, ide, jds, jde, kds, kde, &
4742 ims, ime, jms, jme, kms, kme, &
4743 its, ite, jts, jte, kts, kte )
4744 !----------------------------------------------------------------------
4746 !----------------------------------------------------------------------
4748 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
4750 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
4751 ims, ime, jms, jme, kms, kme, &
4752 its, ite, jts, jte, kts, kte
4753 INTEGER , INTENT(IN ) :: n_moist
4755 REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist
4758 REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut, muu, muv
4760 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4761 INTENT( OUT) :: u_phy, &
4775 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4776 INTENT( OUT) :: p_hyd, &
4779 REAL , INTENT(IN ) :: p_top
4781 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , &
4782 INTENT(IN ) :: pb, &
4792 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
4795 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: znw, &
4798 REAL, DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1f, c2f
4800 REAL, DIMENSION( kms:kme ) :: c1, c2
4801 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
4803 REAL :: w1, w2, z0, z1, z2
4807 !-----------------------------------------------------------------------
4811 ! phys_prep calculates a number of diagnostic quantities needed by
4812 ! the physics routines.
4818 ! set up loop bounds for this grid's boundary conditions
4821 i_end = min( ite,ide-1 )
4823 j_end = min( jte,jde-1 )
4826 k_end = min( kte, kde-1 )
4828 ! compute thermodynamics and velocities at pressure points (or half levels)
4830 IF ( ( config_flags%use_theta_m .EQ. 1 ) .AND. (P_Qv .GE. PARAM_FIRST_SCALAR) ) THEN
4831 do j = j_start,j_end
4832 do k = k_start, k_end
4833 do i = i_start, i_end
4834 th_phy(i,k,j) = (t(i,k,j)+t0)/(1.+R_v/R_d*moist(i,k,j,P_QV))
4839 do j = j_start,j_end
4840 do k = k_start, k_end
4841 do i = i_start, i_end
4842 th_phy(i,k,j) = t(i,k,j)+t0
4848 do j = j_start,j_end
4849 do k = k_start, k_end
4850 do i = i_start, i_end
4852 th_phy_m_t0(i,k,j) = th_phy(i,k,j)-t0
4853 p_phy(i,k,j) = p(i,k,j) + pb(i,k,j)
4854 pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp
4855 t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j)
4856 rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))
4857 u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j))
4858 v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1))
4864 ! compute z at w points
4866 do j = j_start,j_end
4868 do i = i_start, i_end
4869 z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g
4874 do j = j_start,j_end
4875 do k = k_start, kte-1
4876 do i = i_start, i_end
4877 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
4882 do j = j_start,j_end
4883 do i = i_start, i_end
4888 ! compute z at p points or half levels (average of z at full levels)
4890 do j = j_start,j_end
4891 do k = k_start, k_end
4892 do i = i_start, i_end
4893 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
4898 ! interp t and p to full levels
4900 do j = j_start,j_end
4902 do i = i_start, i_end
4903 p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j)
4904 t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j)
4909 ! extrapolate p and t to surface and top.
4910 ! we'll use an extrapolation in z for now
4912 do j = j_start,j_end
4913 do i = i_start, i_end
4920 w1 = (z0 - z2)/(z1 - z2)
4922 p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j)
4923 t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j)
4927 z0 = z_at_w(i,kte,j)
4930 w1 = (z0 - z2)/(z1 - z2)
4933 ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j)
4934 !!! bug fix extrapolate ln(p) so p is positive definite
4935 p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j)))
4936 t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j)
4941 ! calculate hydrostatic pressure at both full and half levels
4942 ! first, full level p: assuming dry over model top
4944 do j = j_start,j_end
4945 do i = i_start, i_end
4946 p_hyd_w(i,kte,j) = p_top
4950 do j = j_start,j_end
4951 do k = kte-1, k_start, -1
4952 do i = i_start, i_end
4954 do n = PARAM_FIRST_SCALAR,n_moist
4955 qtot = qtot + moist(i,k,j,n)
4957 p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j) - (1.+qtot)*(c1(k)*MUT(i,j)+c2(k))*dnw(k)
4958 ! p_hyd_w(i,k,j) = p_hyd_w(i,k+1,j)+1./alt(i,k,j)*(1.+moist(i,k,j,P_QV))*g*dz8w(i,k,j)
4963 ! now calculate hydrostatic pressure at half levels
4965 do j = j_start,j_end
4966 do k = k_start, k_end
4967 do i = i_start, i_end
4968 p_hyd(i,k,j) = 0.5*(p_hyd_w(i,k,j)+p_hyd_w(i,k+1,j))
4973 END SUBROUTINE phy_prep
4976 !======================================================================
4977 ! routine to decouple physics tendencies
4978 !======================================================================
4980 SUBROUTINE phy_prep_part2 ( config_flags, &
4982 c1h, c2h, c1f, c2f, &
4984 RTHBLTEN, RUBLTEN, RVBLTEN, &
4985 RQVBLTEN, RQCBLTEN, RQIBLTEN, &
4986 RUCUTEN, RVCUTEN, RTHCUTEN, &
4987 RQVCUTEN, RQCCUTEN, RQRCUTEN, &
4988 RQICUTEN, RQSCUTEN, &
4989 RUSHTEN, RVSHTEN, RTHSHTEN, &
4990 RQVSHTEN, RQCSHTEN, RQRSHTEN, &
4991 RQISHTEN, RQSSHTEN, RQGSHTEN, &
4993 RUNDGDTEN, RVNDGDTEN, RTHNDGDTEN, &
4994 RPHNDGDTEN,RQVNDGDTEN, RMUNDGDTEN, &
4995 t_new, th_phy, qv, &
4996 ids, ide, jds, jde, kds, kde, &
4997 ims, ime, jms, jme, kms, kme, &
4998 its, ite, jts, jte, kts, kte )
4999 !----------------------------------------------------------------------
5001 !----------------------------------------------------------------------
5003 TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags
5005 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5006 ims, ime, jms, jme, kms, kme, &
5007 its, ite, jts, jte, kts, kte
5009 REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut, muu, muv
5011 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5012 INTENT(INOUT) :: RTHRATEN
5014 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5015 INTENT(INOUT) :: RUCUTEN, &
5033 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
5034 INTENT(INOUT) :: RUBLTEN, &
5041 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
5042 INTENT(INOUT) :: RTHFTEN, &
5045 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
5046 INTENT(INOUT) :: RUNDGDTEN, &
5052 REAL, DIMENSION( ims:ime, jms:jme ) , &
5053 INTENT(INOUT) :: RMUNDGDTEN
5055 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
5056 INTENT(IN ) :: t_new, qv
5058 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
5059 INTENT(INOUT) :: th_phy
5061 REAL, DIMENSION( kms:kme ) , INTENT(IN ) :: c1h, c2h, c1f, c2f
5063 REAL, DIMENSION( kms:kme ) :: c1, c2
5065 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end, i_startu, j_startv
5068 !-----------------------------------------------------------------------
5072 ! It decouples the physics tendencies from
5073 ! the column dry-air mass (the physics routines expect to see/update the
5074 ! uncoupled tendencies).
5078 ! set up loop bounds for this grid's boundary conditions
5081 i_end = min( ite,ide-1 )
5083 j_end = min( jte,jde-1 )
5086 k_end = min( kte, kde-1 )
5091 ! decouple all physics tendencies
5093 IF (config_flags%ra_lw_physics .gt. 0 .or. config_flags%ra_sw_physics .gt. 0) THEN
5098 RTHRATEN(I,K,J)=RTHRATEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5105 IF (config_flags%cu_physics .gt. 0) THEN
5110 RUCUTEN(I,K,J) =RUCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5111 RVCUTEN(I,K,J) =RVCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5112 RTHCUTEN(I,K,J)=RTHCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5117 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5121 RQVCUTEN(I,K,J)=RQVCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5127 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
5131 RQCCUTEN(I,K,J)=RQCCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5137 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
5141 RQRCUTEN(I,K,J)=RQRCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5147 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
5151 RQICUTEN(I,K,J)=RQICUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5157 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
5161 RQSCUTEN(I,K,J)=RQSCUTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5169 IF (config_flags%shcu_physics .gt. 0) THEN
5174 RUSHTEN(I,K,J) =RUSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5175 RVSHTEN(I,K,J) =RVSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5176 RTHSHTEN(I,K,J)=RTHSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5181 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5185 RQVSHTEN(I,K,J)=RQVSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5191 IF (P_QC .ge. PARAM_FIRST_SCALAR)THEN
5195 RQCSHTEN(I,K,J)=RQCSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5201 IF (P_QR .ge. PARAM_FIRST_SCALAR)THEN
5205 RQRSHTEN(I,K,J)=RQRSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5211 IF (P_QI .ge. PARAM_FIRST_SCALAR)THEN
5215 RQISHTEN(I,K,J)=RQISHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5221 IF(P_QS .ge. PARAM_FIRST_SCALAR)THEN
5225 RQSSHTEN(I,K,J)=RQSSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5231 IF(P_QG .ge. PARAM_FIRST_SCALAR)THEN
5235 RQGSHTEN(I,K,J)=RQGSHTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5243 IF (config_flags%bl_pbl_physics .gt. 0) THEN
5248 RUBLTEN(I,K,J) =RUBLTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5249 RVBLTEN(I,K,J) =RVBLTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5250 RTHBLTEN(I,K,J)=RTHBLTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5255 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5259 RQVBLTEN(I,K,J)=RQVBLTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5265 IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
5269 RQCBLTEN(I,K,J)=RQCBLTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5275 IF (P_QI .ge. PARAM_FIRST_SCALAR) THEN
5279 RQIBLTEN(I,K,J)=RQIBLTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5287 ! decouple advective forcing required by a few CPS schemes
5289 IF(( config_flags%cu_physics == GDSCHEME ) .OR. &
5290 ( config_flags%cu_physics == GFSCHEME ) .OR. &
5291 ( config_flags%cu_physics == G3SCHEME ) .OR. &
5292 ( config_flags%cu_physics == KFETASCHEME ) .OR. &
5293 ( config_flags%cu_physics == MSKFSCHEME ) .OR. &
5294 ( config_flags%cu_physics == TIEDTKESCHEME ) .OR. &
5295 ( config_flags%cu_physics == NTIEDTKESCHEME )) THEN
5297 IF (P_QV .ge. PARAM_FIRST_SCALAR)THEN
5301 RQVFTEN(I,K,J)=RQVFTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5309 IF(( config_flags%cu_physics == GDSCHEME ) .OR. &
5310 ( config_flags%cu_physics == GFSCHEME ) .OR. &
5311 ( config_flags%cu_physics == G3SCHEME ) .OR. &
5312 ( config_flags%cu_physics == NTIEDTKESCHEME )) THEN
5317 RTHFTEN(I,K,J)=RTHFTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5322 ! If using moist theta, get dry theta tendency for CPSs
5323 IF ( config_flags%use_theta_m == 1 ) THEN
5327 th_phy(i,k,j) = (t_new(i,k,j) + t0) / (1. + (R_v/R_d) * qv(i,k,j))
5328 rthften(i,k,j) = th_phy(i,k,j)/(t_new(i,k,j)+t0) * &
5329 (rthften(i,k,j) - (R_v/R_d) * th_phy(i,k,j) * rqvften(i,k,j))
5337 ! note fdda u and v tendencies are staggered, also only interior points have muu/muv,
5338 ! so only decouple those
5340 IF (config_flags%grid_fdda .gt. 0) THEN
5342 i_startu=MAX(its,ids+1)
5343 j_startv=MAX(jts,jds+1)
5348 RUNDGDTEN(I,K,J) =RUNDGDTEN(I,K,J)/(c1h(k)*muu(I,J)+c2h(k))
5355 RVNDGDTEN(I,K,J) =RVNDGDTEN(I,K,J)/(c1h(k)*muv(I,J)+c2h(k))
5362 RTHNDGDTEN(I,K,J)=RTHNDGDTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5363 ! RMUNDGDTEN(I,J) - no coupling
5368 IF (config_flags%grid_fdda .EQ. 2) THEN
5372 RPHNDGDTEN(I,K,J)=RPHNDGDTEN(I,K,J)/(c1f(k)*mut(I,J)+c2f(k))
5377 ELSE IF (config_flags%grid_fdda .GE. 1) THEN
5378 IF (P_QV .ge. PARAM_FIRST_SCALAR) THEN
5382 RQVNDGDTEN(I,K,J)=RQVNDGDTEN(I,K,J)/(c1(k)*MUT(I,J)+c2(k))
5391 END SUBROUTINE phy_prep_part2
5392 !------------------------------------------------------------
5394 SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, &
5395 p, p8w, p0, pb, ph, phb, &
5401 config_flags,fzm, fzp, &
5402 ids,ide, jds,jde, kds,kde, &
5403 ims,ime, jms,jme, kms,kme, &
5404 its,ite, jts,jte, kts,kte )
5408 ! Here we construct full fields
5409 ! needed by the microphysics
5411 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5413 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5414 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5415 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5417 REAL, INTENT(IN ) :: dt
5419 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5420 INTENT(IN ) :: al, &
5430 REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, &
5433 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5434 INTENT( OUT) :: rho, &
5443 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5444 INTENT(INOUT) :: h_diabatic, &
5448 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5449 INTENT(INOUT) :: t_new, &
5452 REAL, INTENT(IN ) :: t0, p0
5453 REAL :: z0,z1,z2,w1,w2
5455 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5458 !--------------------------------------------------------------------
5462 ! moist_phys_prep_em calculates a number of diagnostic quantities needed by
5463 ! the microphysics routines.
5467 ! set up loop bounds for this grid's boundary conditions
5470 i_end = min( ite,ide-1 )
5472 j_end = min( jte,jde-1 )
5475 k_end = min( kte, kde-1 )
5477 DO j = j_start, j_end
5479 DO i = i_start, i_end
5480 z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
5485 do j = j_start,j_end
5486 do k = k_start, kte-1
5487 do i = i_start, i_end
5488 dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j)
5493 do j = j_start,j_end
5494 do i = i_start, i_end
5500 ! compute full pii, rho, and z at the new time-level
5501 ! (needed for physics).
5502 ! convert perturbation theta to full theta (th_phy)
5503 ! use h_diabatic to temporarily save pre-microphysics full theta
5505 IF ( ( config_flags%use_theta_m .EQ. 1 ) .AND. (P_Qv .GE. PARAM_FIRST_SCALAR) ) THEN
5506 DO j = j_start, j_end
5507 DO k = k_start, k_end
5508 DO i = i_start, i_end
5509 th_phy(i,k,j) = (t_new(i,k,j) + t0) / (1. + (R_v/R_d) * qv(i,k,j))
5514 DO j = j_start, j_end
5515 DO k = k_start, k_end
5516 DO i = i_start, i_end
5517 th_phy(i,k,j) = t_new(i,k,j) + t0
5523 DO j = j_start, j_end
5524 DO k = k_start, k_end
5525 DO i = i_start, i_end
5526 h_diabatic(i,k,j) = th_phy(i,k,j)
5527 #if ( WRFPLUS == 1 )
5528 if ( P_QV >= PARAM_FIRST_SCALAR ) then
5529 qv_diabatic(i,k,j) = qv(i,k,j)
5531 qv_diabatic(i,k,j) = 0.0
5533 if ( P_QC >= PARAM_FIRST_SCALAR ) then
5534 qc_diabatic(i,k,j) = qc(i,k,j)
5536 qc_diabatic(i,k,j) = 0.0
5539 qv_diabatic(i,k,j) = qv(i,k,j)
5540 qc_diabatic(i,k,j) = qc(i,k,j)
5542 rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j))
5543 pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp
5544 z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) )
5545 pf(i,k,j) = p(i,k,j)+pb(i,k,j)
5551 ! interp p at w points
5553 do j = j_start,j_end
5555 do i = i_start, i_end
5556 p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j)
5561 ! extrapolate p to surface and top.
5562 ! we'll use an extrapolation in z for now
5564 do j = j_start,j_end
5565 do i = i_start, i_end
5572 w1 = (z0 - z2)/(z1 - z2)
5574 p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j)
5578 z0 = z_at_w(i,kte,j)
5581 w1 = (z0 - z2)/(z1 - z2)
5583 ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j)
5584 p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j)))
5589 END SUBROUTINE moist_physics_prep_em
5591 !------------------------------------------------------------------------------
5593 SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, &
5594 th_phy, h_diabatic, dt, &
5599 #if ( WRF_DFI_RADAR == 1 )
5600 dfi_tten_rad,dfi_stage, &
5602 ids,ide, jds,jde, kds,kde, &
5603 ims,ime, jms,jme, kms,kme, &
5604 its,ite, jts,jte, kts,kte )
5608 ! Here we construct full fields
5609 ! needed by the microphysics
5611 TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags
5613 INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde
5614 INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme
5615 INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte
5617 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5618 INTENT(INOUT) :: t_new, &
5625 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5626 INTENT( OUT) :: th_phy_m_t0
5628 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5629 INTENT(IN ) :: qv, &
5632 #if ( WRF_DFI_RADAR == 1 )
5633 REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), &
5634 INTENT(IN), OPTIONAL :: dfi_tten_rad
5635 INTEGER, INTENT(IN ) ,OPTIONAL :: dfi_stage
5636 REAL :: dfi_tten_max, old_max
5639 REAL mpten, mptenmax, mptenmin
5642 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut
5645 REAL, INTENT(IN ) :: t0, dt
5647 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end
5648 INTEGER :: i, j, k, imax, jmax, imin, jmin
5650 !--------------------------------------------------------------------
5654 ! moist_phys_finish_em resets theta to its perturbation value and
5655 ! computes and stores the microphysics diabatic heating term.
5659 ! set up loop bounds for this grid's boundary conditions
5663 i_end = min( ite,ide-1 )
5665 j_end = min( jte,jde-1 )
5667 k_end = min( kte, kde-1 )
5669 #if ( WRF_DFI_RADAR == 1 )
5670 IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5671 IF ( dfi_stage ==DFI_FWD ) THEN
5672 WRITE(wrf_err_message,*)'Add radar tendency: i_start,j_start: ', i_start, j_start
5673 CALL wrf_debug ( 100 , TRIM(wrf_err_message) )
5680 ! add microphysics theta diff to perturbation theta, set h_diabatic
5682 IF ( config_flags%no_mp_heating .eq. 0 ) THEN
5685 DO j = j_start, j_end
5686 DO k = k_start, k_end
5687 DO i = i_start, i_end
5688 mpten = th_phy(i,k,j)-h_diabatic(i,k,j)
5689 #if ( WRFPLUS == 1 )
5690 if ( P_QV >= PARAM_FIRST_SCALAR ) qvten = qv(i,k,j)-qv_diabatic(i,k,j)
5691 if ( P_QC >= PARAM_FIRST_SCALAR ) qcten = qc(i,k,j)-qc_diabatic(i,k,j)
5693 qvten = qv(i,k,j)-qv_diabatic(i,k,j)
5694 qcten = qc(i,k,j)-qc_diabatic(i,k,j)
5696 if(mpten.gt.mptenmax) then
5701 if(mpten.lt.mptenmin) then
5706 mpten=min(config_flags%mp_tend_lim*dt, mpten)
5707 mpten=max(-config_flags%mp_tend_lim*dt, mpten)
5709 #if ( WRF_DFI_RADAR == 1 )
5710 !compiler error, not handled yet
5712 if(dfi_tten_max < dfi_tten_rad(i,k,j) ) dfi_tten_max = dfi_tten_rad(i,k,j)
5713 if(old_max < (th_phy(i,k,j)-h_diabatic(i,k,j)) ) old_max=th_phy(i,k,j)-h_diabatic(i,k,j)
5716 IF ( PRESENT(dfi_stage) .and. PRESENT(dfi_tten_rad) ) THEN
5717 IF ( dfi_stage == DFI_FWD .and. dfi_tten_rad(i,k,j) >= -0.1 .and. &
5718 dfi_tten_rad(i,k,j) <= 0.1 .and. k < k_end ) THEN
5719 ! add radar temp tendency
5720 ! there is radar coverage
5721 t_new(i,k,j) = t_new(i,k,j) + (dfi_tten_rad(i,k,j))*dt
5724 t_new(i,k,j) = t_new(i,k,j) + mpten
5727 th_phy_m_t0(i,k,j) = t_new(i,k,j)
5728 h_diabatic(i,k,j) = mpten / dt
5730 ! pertubation theta_moist(new) = theta_dry(old)*(1+rv/rd qv(old)) + ! term 1
5731 ! delta theta_dry * (1+rv/rd qv(new)) + ! term 2
5732 ! (rv/rd)*delta qv * theta_dry(new) - ! term 3
5735 IF ( ( config_flags%use_theta_m .EQ. 1 ) .AND. (P_Qv .GE. PARAM_FIRST_SCALAR) ) THEN
5736 t_new(i,k,j) = h_diabatic(i,k,j)*(1. + (R_v/R_d)*qv_diabatic(i,k,j)) + &
5737 mpten*(1. + (R_v/R_d)*qv(i,k,j)) + &
5738 (R_v/R_d)*qvten*th_phy(i,k,j) - T0
5739 th_phy_m_t0(i,k,j) = (t_new(i,k,j)+T0)/(1.+(R_v/R_d)*qv(i,k,j)) - T0
5740 h_diabatic(i,k,j) = ( mpten*(1. + (R_v/R_d)*qv(i,k,j)) + &
5741 (R_v/R_d)*qvten*th_phy(i,k,j) ) / dt
5743 t_new(i,k,j) = t_new(i,k,j) + mpten
5744 th_phy_m_t0(i,k,j) = t_new(i,k,j)
5745 h_diabatic(i,k,j) = mpten/dt
5751 #if ( WRFPLUS == 1 )
5752 if ( P_QV >= PARAM_FIRST_SCALAR ) then
5753 qv_diabatic(i,k,j) = qvten/dt
5755 qv_diabatic(i,k,j) = 0.0
5757 if ( P_QC >= PARAM_FIRST_SCALAR ) then
5758 qc_diabatic(i,k,j) = qcten/dt
5760 qc_diabatic(i,k,j) = 0.0
5763 qv_diabatic(i,k,j) = qvten/dt
5764 qc_diabatic(i,k,j) = qcten/dt
5772 DO j = j_start, j_end
5773 DO k = k_start, k_end
5774 DO i = i_start, i_end
5775 ! t_new(i,k,j) = t_new(i,k,j)
5776 h_diabatic(i,k,j) = 0.
5777 qv_diabatic(i,k,j) = 0.
5778 qc_diabatic(i,k,j) = 0.
5784 END SUBROUTINE moist_physics_finish_em
5786 !----------------------------------------------------------------
5789 SUBROUTINE init_module_big_step
5790 END SUBROUTINE init_module_big_step
5792 SUBROUTINE set_tend ( field, field_adv_tend, msf, &
5793 ids, ide, jds, jde, kds, kde, &
5794 ims, ime, jms, jme, kms, kme, &
5795 its, ite, jts, jte, kts, kte )
5801 INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, &
5802 ims, ime, jms, jme, kms, kme, &
5803 its, ite, jts, jte, kts, kte
5805 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT) :: field
5807 REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN) :: field_adv_tend
5809 REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: msf
5813 INTEGER :: i, j, k, itf, jtf, ktf
5817 ! set_tend copies the advective tendency array into the tendency array.
5821 jtf = MIN(jte,jde-1)
5822 ktf = MIN(kte,kde-1)
5823 itf = MIN(ite,ide-1)
5827 field(i,k,j) = field_adv_tend(i,k,j)*msf(i,j)
5832 END SUBROUTINE set_tend
5834 !------------------------------------------------------------------------------
5836 SUBROUTINE rk_rayleigh_damp( ru_tendf, rv_tendf, &
5837 rw_tendf, t_tendf, &
5838 u, v, w, t, t_init, &
5839 c1h, c2h, c1f, c2f, &
5840 mut, muu, muv, ph, phb, &
5841 u_base, v_base, t_base, z_base, &
5843 ids, ide, jds, jde, kds, kde, &
5844 ims, ime, jms, jme, kms, kme, &
5845 its, ite, jts, jte, kts, kte )
5847 ! History: Apr 2005 Modifications by George Bryan, NCAR:
5848 ! - Generalized the code in a way that allows for
5849 ! simulations with steep terrain.
5851 ! Jul 2004 Modifications by George Bryan, NCAR:
5852 ! - Modified the code to use u_base, v_base, and t_base
5853 ! arrays for the background state. Removed the hard-wired
5854 ! base-state values.
5855 ! - Modified the code to use dampcoef, zdamp, and damp_opt,
5856 ! i.e., the upper-level damper variables in namelist.input.
5857 ! Removed the hard-wired variables in the older version.
5858 ! This damper is used when damp_opt = 2.
5859 ! - Modified the code to account for the movement of the
5860 ! model surfaces with time. The code now obtains a base-
5861 ! state value by interpolation using the "_base" arrays.
5863 ! Nov 2003 Bug fix by Jason Knievel, NCAR
5865 ! Aug 2003 Meridional dimension, some comments, and
5866 ! changes in layout of the code added by
5867 ! Jason Knievel, NCAR
5869 ! Jul 2003 Original code by Bill Skamarock, NCAR
5871 ! Purpose: This routine applies Rayleigh damping to a layer at top
5872 ! of the model domain.
5874 !-----------------------------------------------------------------------
5875 ! Begin declarations.
5879 INTEGER, INTENT( IN ) &
5880 :: ids, ide, jds, jde, kds, kde, &
5881 ims, ime, jms, jme, kms, kme, &
5882 its, ite, jts, jte, kts, kte
5884 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
5885 :: ru_tendf, rv_tendf, rw_tendf, t_tendf
5887 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
5888 :: u, v, w, t, t_init, ph, phb
5890 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
5893 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5894 :: u_base, v_base, t_base, z_base
5896 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
5897 :: c1h, c2h, c1f, c2f
5905 :: i_start, i_end, j_start, j_end, k_start, k_end, i, j, k, ktf, k1, k2
5907 REAL, DIMENSION( kms:kme ) &
5911 :: pii, dcoef, z, ztop
5913 REAL :: wkp1, wk, wkm1
5915 REAL, DIMENSION( kms:kme ) :: z00, u00, v00, t00
5918 !-----------------------------------------------------------------------
5923 pii = 2.0 * asin(1.0)
5925 ktf = MIN( kte, kde-1 )
5927 !-----------------------------------------------------------------------
5928 ! Adjust u to base state.
5930 DO j = jts, MIN( jte, jde-1 )
5931 DO i = its, MIN( ite, ide )
5933 ! Get height at top of model
5934 ztop = 0.5*( phb(i ,kde,j)+phb(i-1,kde,j) &
5935 +ph(i ,kde,j)+ph(i-1,kde,j) )/g
5937 ! Find bottom of damping layer
5940 DO WHILE( z >= (ztop-zdamp) )
5941 z = 0.25*( phb(i ,k1,j)+phb(i ,k1+1,j) &
5942 +phb(i-1,k1,j)+phb(i-1,k1+1,j) &
5943 +ph(i ,k1,j)+ph(i ,k1+1,j) &
5944 +ph(i-1,k1,j)+ph(i-1,k1+1,j))/g
5950 ! Get reference state at model levels
5953 DO WHILE( z_base(k2) .gt. z00(k) )
5957 u00(k) = u_base(k2) + ( u_base(k2) - u_base(k2-1) ) &
5958 * ( z00(k) - z_base(k2) ) &
5959 / ( z_base(k2) - z_base(k2-1) )
5961 u00(k) = u_base(k2) + ( u_base(k2+1) - u_base(k2) ) &
5962 * ( z00(k) - z_base(k2) ) &
5963 / ( z_base(k2+1) - z_base(k2) )
5967 ! Apply the Rayleigh damper
5969 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
5970 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
5971 ru_tendf(i,k,j) = ru_tendf(i,k,j) - &
5972 (c1h(k)*muu(i,j)+c2h(k)) * ( dcoef * dampcoef ) * &
5973 ( u(i,k,j) - u00(k) )
5979 ! End adjustment of u.
5980 !-----------------------------------------------------------------------
5982 !-----------------------------------------------------------------------
5983 ! Adjust v to base state.
5985 DO j = jts, MIN( jte, jde )
5986 DO i = its, MIN( ite, ide-1 )
5988 ! Get height at top of model
5989 ztop = 0.5*( phb(i,kde,j )+phb(i,kde,j-1) &
5990 +ph(i,kde,j )+ph(i,kde,j-1) )/g
5992 ! Find bottom of damping layer
5995 DO WHILE( z >= (ztop-zdamp) )
5996 z = 0.25*( phb(i,k1,j )+phb(i,k1+1,j ) &
5997 +phb(i,k1,j-1)+phb(i,k1+1,j-1) &
5998 +ph(i,k1,j )+ph(i,k1+1,j ) &
5999 +ph(i,k1,j-1)+ph(i,k1+1,j-1))/g
6005 ! Get reference state at model levels
6008 DO WHILE( z_base(k2) .gt. z00(k) )
6012 v00(k) = v_base(k2) + ( v_base(k2) - v_base(k2-1) ) &
6013 * ( z00(k) - z_base(k2) ) &
6014 / ( z_base(k2) - z_base(k2-1) )
6016 v00(k) = v_base(k2) + ( v_base(k2+1) - v_base(k2) ) &
6017 * ( z00(k) - z_base(k2) ) &
6018 / ( z_base(k2+1) - z_base(k2) )
6022 ! Apply the Rayleigh damper
6024 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
6025 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
6026 rv_tendf(i,k,j) = rv_tendf(i,k,j) - &
6027 (c1h(k)*muv(i,j)+c2h(k)) * ( dcoef * dampcoef ) * &
6028 ( v(i,k,j) - v00(k) )
6034 ! End adjustment of v.
6035 !-----------------------------------------------------------------------
6037 !-----------------------------------------------------------------------
6038 ! Adjust w to base state.
6040 DO j = jts, MIN( jte, jde-1 )
6041 DO i = its, MIN( ite, ide-1 )
6042 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
6043 DO k = kts, MIN( kte, kde )
6044 z = ( phb(i,k,j) + ph(i,k,j) ) / g
6045 IF ( z >= (ztop-zdamp) ) THEN
6046 dcoef = 1.0 - MIN( 1.0, ( ztop - z ) / zdamp )
6047 dcoef = ( SIN( 0.5 * pii * dcoef ) )**2
6048 rw_tendf(i,k,j) = rw_tendf(i,k,j) - &
6049 (c1f(k)*mut(i,j)+c2f(k)) * ( dcoef * dampcoef ) * w(i,k,j)
6055 ! End adjustment of w.
6056 !-----------------------------------------------------------------------
6058 !-----------------------------------------------------------------------
6059 ! Adjust potential temperature to base state.
6061 DO j = jts, MIN( jte, jde-1 )
6062 DO i = its, MIN( ite, ide-1 )
6064 ! Get height at top of model
6065 ztop = ( phb(i,kde,j) + ph(i,kde,j) ) / g
6067 ! Find bottom of damping layer
6070 DO WHILE( z >= (ztop-zdamp) )
6071 z = 0.5 * ( phb(i,k1,j) + phb(i,k1+1,j) + &
6072 ph(i,k1,j) + ph(i,k1+1,j) ) / g
6078 ! Get reference state at model levels
6081 DO WHILE( z_base(k2) .gt. z00(k) )
6085 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
6086 * ( z00(k) - z_base(k2) ) &
6087 / ( z_base(k2) - z_base(k2-1) )
6089 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
6090 * ( z00(k) - z_base(k2) ) &
6091 / ( z_base(k2+1) - z_base(k2) )
6095 ! Apply the Rayleigh damper
6097 dcoef = 1.0 - MIN( 1.0, ( ztop - z00(k) ) / zdamp )
6098 dcoef = (SIN( 0.5 * pii * dcoef ) )**2
6099 t_tendf(i,k,j) = t_tendf(i,k,j) - &
6100 (c1(k)*MUT(i,j)+c2(k)) * ( dcoef * dampcoef ) * &
6101 ( t(i,k,j) - t00(k) )
6107 ! End adjustment of potential temperature.
6108 !-----------------------------------------------------------------------
6110 END SUBROUTINE rk_rayleigh_damp
6112 !==============================================================================
6113 !==============================================================================
6115 SUBROUTINE theta_relaxation( t_tendf, t, t_init, &
6116 MUT, c1, c2, ph, phb, &
6118 ids, ide, jds, jde, kds, kde, &
6119 ims, ime, jms, jme, kms, kme, &
6120 its, ite, jts, jte, kts, kte )
6122 ! Purpose: Newtonian relaxation on potential temperature. Serves two
6123 ! purposes: 1) to mimic atmospheric radiation in a simple
6124 ! manner, and 2) to keep the vertical profile of temperature
6125 ! close to the initial (base-state) profile, which is useful
6126 ! for certain idealized applications.
6128 ! Reference: Rotunno and Emanuel, 1987, JAS, p. 546
6130 !-----------------------------------------------------------------------
6131 ! Begin declarations.
6135 INTEGER, INTENT( IN ) &
6136 :: ids, ide, jds, jde, kds, kde, &
6137 ims, ime, jms, jme, kms, kme, &
6138 its, ite, jts, jte, kts, kte
6140 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( INOUT ) &
6143 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT( IN ) &
6144 :: t, t_init, ph, phb
6146 REAL, DIMENSION( ims:ime, jms:jme ), INTENT( IN ) &
6149 REAL, DIMENSION( kms:kme), INTENT( IN ) &
6152 REAL, DIMENSION( kms:kme ) , INTENT(IN ) &
6157 INTEGER :: i, j, k, ktf, k2
6158 REAL :: tau_r , rmax , rmin , inv_tau_r , inv_g , rterm
6159 REAL, DIMENSION( kms:kme ) :: z00,t00
6162 !-----------------------------------------------------------------------
6164 ! set tau_r to 12 h, following RE87
6167 ! limit rterm to +/- 2 K/day
6171 ktf = MIN( kte, kde-1 )
6172 inv_tau_r = 1.0/tau_r
6175 !-----------------------------------------------------------------------
6176 ! Adjust potential temperature to base state.
6178 DO j = jts, MIN( jte, jde-1 )
6179 DO i = its, MIN( ite, ide-1 )
6181 ! Get height of model levels:
6183 z00(k) = 0.5 * ( phb(i,k,j) + phb(i,k+1,j) + &
6184 ph(i,k,j) + ph(i,k+1,j) ) * inv_g
6187 ! Get reference state:
6190 DO WHILE( z_base(k2) .gt. z00(k) .and. k2 .gt. 1 )
6194 t00(k) = t_base(k2) + ( t_base(k2) - t_base(k2-1) ) &
6195 * ( z00(k) - z_base(k2) ) &
6196 / ( z_base(k2) - z_base(k2-1) )
6198 t00(k) = t_base(k2) + ( t_base(k2+1) - t_base(k2) ) &
6199 * ( z00(k) - z_base(k2) ) &
6200 / ( z_base(k2+1) - z_base(k2) )
6204 ! Apply the RE87 R term:
6206 rterm = -( t(i,k,j) - t00(k) )*inv_tau_r
6208 rterm = min( rterm , rmax )
6209 rterm = max( rterm , rmin )
6210 t_tendf(i,k,j) = t_tendf(i,k,j) + (c1(k)*MUT(i,j)+c2(k))*rterm
6216 END SUBROUTINE theta_relaxation
6218 !==============================================================================
6220 SUBROUTINE sixth_order_diffusion( name, field, tendency, MUT, dt, &
6221 config_flags, c1, c2, &
6222 diff_6th_opt, diff_6th_factor, &
6228 ids, ide, jds, jde, kds, kde, &
6229 ims, ime, jms, jme, kms, kme, &
6230 its, ite, jts, jte, kts, kte )
6232 ! History: 14 Nov 2006 Name of variable changed by Jason Knievel
6233 ! 07 Jun 2006 Revised and generalized by Jason Knievel
6234 ! 25 Apr 2005 Original code by Jason Knievel, NCAR
6236 ! Purpose: Apply 6th-order, monotonic (flux-limited), numerical
6237 ! diffusion to 3-d velocity and to scalars.
6239 ! References: Ming Xue (MWR Aug 2000)
6240 ! Durran ("Numerical Methods for Wave Equations..." 1999)
6241 ! George Bryan (personal communication)
6243 !------------------------------------------------------------------------------
6244 ! Begin: Declarations.
6248 INTEGER, INTENT(IN) &
6249 :: ids, ide, jds, jde, kds, kde, &
6250 ims, ime, jms, jme, kms, kme, &
6251 its, ite, jts, jte, kts, kte
6253 TYPE(grid_config_rec_type), INTENT(IN) &
6256 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(INOUT) &
6259 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
6262 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), INTENT(IN) &
6265 REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN) &
6268 REAL, DIMENSION( kms:kme ), INTENT(IN) &
6277 INTEGER, INTENT(IN) &
6280 CHARACTER(LEN=1) , INTENT(IN) &
6290 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfux
6291 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfuy
6292 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvx
6293 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfvy
6294 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msftx
6295 REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfty
6298 :: dflux_x_p0, dflux_y_p0, &
6299 dflux_x_p1, dflux_y_p1, &
6300 tendency_x, tendency_y, &
6301 mu_avg_p0, mu_avg_p1, &
6305 :: dzmax_p0, dzmax_p1, dx, dy, slopedamp_p0, slopedamp_p1, dzthresh
6310 ! End: Declarations.
6311 !------------------------------------------------------------------------------
6313 !------------------------------------------------------------------------------
6314 ! Begin: Translate the diffusion factor into a diffusion coefficient. See
6315 ! Durran's text, section 2.4.3, then adjust for sixth-order diffusion (not
6316 ! fourth) and for diffusion in two dimensions (not one). For reference, a
6317 ! factor of 1.0 would mean complete diffusion of a 2dx wave in one time step,
6318 ! although application of the flux limiter reduces somewhat the effects of
6319 ! diffusion for a given coefficient.
6321 diff_6th_coef = diff_6th_factor * 0.015625 / ( 2.0 * dt )
6323 ! End: Translate diffusion factor.
6324 !------------------------------------------------------------------------------
6326 !------------------------------------------------------------------------------
6327 ! Begin: Assign limits of spatial loops depending on variable to be diffused.
6328 ! The halo regions outside the physical domain are not defined, hence the loop
6329 ! indices should start 3 rows and columns in for specified boundary conditions.
6330 ! Also modified to work with idealized boundary conditions.
6332 ktf = MIN( kte, kde-1 )
6334 if(config_flags%specified .or. config_flags%nested) specified = .true.
6336 IF ( name .EQ. 'u' ) THEN
6341 j_end = MIN(jde-1,jte)
6342 IF ( config_flags%open_xs ) THEN
6343 i_start = MAX(its,ids+3)
6345 IF ( config_flags%open_ys ) THEN
6346 j_start = MAX(jts,jds+3)
6348 IF ( config_flags%open_xe ) THEN
6349 i_end = MIN(ite,ide-3)
6351 IF ( config_flags%open_ye ) THEN
6352 j_end = MIN(jte,jde-4)
6354 IF ( specified ) THEN
6355 i_start = MAX(its,ids+3)
6356 i_end = MIN(ite,ide-3)
6357 j_start = MAX(jts,jds+3)
6358 j_end = MIN(jte,jde-4)
6363 ELSE IF ( name .EQ. 'v' ) THEN
6366 i_end = MIN(ide-1,ite)
6369 IF ( config_flags%open_xs ) THEN
6370 i_start = MAX(its,ids+3)
6372 IF ( config_flags%open_ys ) THEN
6373 j_start = MAX(jts,jds+3)
6375 IF ( config_flags%open_xe ) THEN
6376 i_end = MIN(ite,ide-4)
6378 IF ( config_flags%open_ye ) THEN
6379 j_end = MIN(jte,jde-3)
6381 IF ( config_flags%open_xs .or. specified ) THEN
6382 i_start = MAX(its,ids+3)
6383 i_end = MIN(ite,ide-4)
6384 j_start = MAX(jts,jds+3)
6385 j_end = MIN(jte,jde-3)
6390 ELSE IF ( name .EQ. 'w' ) THEN
6393 i_end = MIN(ide-1,ite)
6395 j_end = MIN(jde-1,jte)
6396 IF ( config_flags%open_xs ) THEN
6397 i_start = MAX(its,ids+3)
6399 IF ( config_flags%open_ys ) THEN
6400 j_start = MAX(jts,jds+3)
6402 IF ( config_flags%open_xe ) THEN
6403 i_end = MIN(ite,ide-4)
6405 IF ( config_flags%open_ye ) THEN
6406 j_end = MIN(jte,jde-4)
6408 IF ( specified ) THEN
6409 i_start = MAX(its,ids+3)
6410 i_end = MIN(ide-4,ite)
6411 j_start = MAX(jts,jds+3)
6412 j_end = MIN(jde-4,jte)
6420 i_end = MIN(ide-1,ite)
6422 j_end = MIN(jde-1,jte)
6423 IF ( config_flags%open_xs ) THEN
6424 i_start = MAX(its,ids+3)
6426 IF ( config_flags%open_ys ) THEN
6427 j_start = MAX(jts,jds+3)
6429 IF ( config_flags%open_xe ) THEN
6430 i_end = MIN(ite,ide-4)
6432 IF ( config_flags%open_ye ) THEN
6433 j_end = MIN(jte,jde-4)
6435 IF ( specified ) THEN
6436 i_start = MAX(its,ids+3)
6437 i_end = MIN(ide-4,ite)
6438 j_start = MAX(jts,jds+3)
6439 j_end = MIN(jde-4,jte)
6446 ! End: Assignment of limits of spatial loops.
6447 !------------------------------------------------------------------------------
6449 !------------------------------------------------------------------------------
6450 ! Begin: Loop across spatial dimensions.
6452 DO j = j_start, j_end
6453 DO k = k_start, k_end
6454 DO i = i_start, i_end
6456 !------------------------------------------------------------------------------
6457 ! Begin: Diffusion in x (i index).
6459 ! Calculate the diffusive flux in x direction (from Xue's eq. 3).
6461 dflux_x_p0 = ( 10.0 * ( field(i, k,j) - field(i-1,k,j) ) &
6462 - 5.0 * ( field(i+1,k,j) - field(i-2,k,j) ) &
6463 + ( field(i+2,k,j) - field(i-3,k,j) ) )
6465 dflux_x_p1 = ( 10.0 * ( field(i+1,k,j) - field(i ,k,j) ) &
6466 - 5.0 * ( field(i+2,k,j) - field(i-1,k,j) ) &
6467 + ( field(i+3,k,j) - field(i-2,k,j) ) )
6469 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6470 ! (variation on Xue's eq. 10).
6472 IF ( diff_6th_opt .EQ. 2 ) THEN
6474 IF ( dflux_x_p0 * ( field(i ,k,j)-field(i-1,k,j) ) .LE. 0.0 ) THEN
6478 IF ( dflux_x_p1 * ( field(i+1,k,j)-field(i ,k,j) ) .LE. 0.0 ) THEN
6487 IF (config_flags%diff_6th_slopeopt .ge. 1) THEN
6489 dzthresh = config_flags%diff_6th_thresh*9.81*dx
6490 IF ( name .EQ. 'u' ) THEN
6491 dzmax_p0 = MAX( ABS(phb(i ,k,j)-phb(i-1,k,j))*msfux(i ,j), ABS(phb(i-1,k,j)-phb(i-2,k,j))*msfux(i-1,j) )
6492 dzmax_p1 = MAX( ABS(phb(i+1,k,j)-phb(i ,k,j))*msfux(i+1,j), ABS(phb(i ,k,j)-phb(i-1,k,j))*msfux(i ,j) )
6493 ELSE IF ( name .EQ. 'v' ) THEN
6494 dzmax_p0 = max( ABS(phb(i ,k,j) - phb(i-1,k,j))*msfux(i ,j), ABS(phb(i ,k,j-1) - phb(i-1,k,j-1))*msfux(i ,j-1) )
6495 dzmax_p1 = max( ABS(phb(i+1,k,j) - phb(i ,k,j))*msfux(i+1,j), ABS(phb(i+1,k,j-1) - phb(i ,k,j-1))*msfux(i+1,j-1) )
6497 dzmax_p0 = ABS(phb(i ,k,j) - phb(i-1,k,j))*msfux(i,j)
6498 dzmax_p1 = ABS(phb(i+1,k,j) - phb(i ,k,j))*msfux(i+1,j)
6500 slopedamp_p0 = MAX(1.- (dzmax_p0/dzthresh) , 0.)
6501 slopedamp_p1 = MAX(1.- (dzmax_p1/dzthresh) , 0.)
6504 ! Apply 6th-order diffusion in x direction.
6506 IF ( name .EQ. 'u' ) THEN
6507 mu_avg_p0 = (c1(k)*MUT(i-1,j)+c2(k))
6508 mu_avg_p1 = (c1(k)*MUT(i ,j)+c2(k))
6509 tendency_x = diff_6th_coef * msfux(i,j) * &
6510 ( ( slopedamp_p1 * mu_avg_p1 * dflux_x_p1 ) - ( slopedamp_p0 * mu_avg_p0 * dflux_x_p0 ) )
6511 ELSE IF ( name .EQ. 'v' ) THEN
6512 mu_avg_p0 = 0.25 * ( &
6513 (c1(k)*MUT(i-1,j-1)+c2(k)) + &
6514 (c1(k)*MUT(i ,j-1)+c2(k)) + &
6515 (c1(k)*MUT(i-1,j )+c2(k)) + &
6516 (c1(k)*MUT(i ,j )+c2(k)) )
6517 mu_avg_p1 = 0.25 * ( &
6518 (c1(k)*MUT(i ,j-1)+c2(k)) + &
6519 (c1(k)*MUT(i+1,j-1)+c2(k)) + &
6520 (c1(k)*MUT(i ,j )+c2(k)) + &
6521 (c1(k)*MUT(i+1,j )+c2(k)) )
6522 tendency_x = diff_6th_coef * msfvx(i,j) * &
6523 ( ( slopedamp_p1 * mu_avg_p1 * dflux_x_p1 ) - ( slopedamp_p0 * mu_avg_p0 * dflux_x_p0 ) )
6525 mu_avg_p0 = 0.5 * ( &
6526 (c1(k)*MUT(i-1,j)+c2(k)) + &
6527 (c1(k)*MUT(i ,j)+c2(k)) )
6528 mu_avg_p1 = 0.5 * ( &
6529 (c1(k)*MUT(i ,j)+c2(k)) + &
6530 (c1(k)*MUT(i+1,j)+c2(k)) )
6531 tendency_x = diff_6th_coef * msftx(i,j) * &
6532 ( ( slopedamp_p1 * mu_avg_p1 * dflux_x_p1 ) - ( slopedamp_p0 * mu_avg_p0 * dflux_x_p0 ) )
6535 ! End: Diffusion in x.
6536 !------------------------------------------------------------------------------
6538 !------------------------------------------------------------------------------
6539 ! Begin: Diffusion in y (j index).
6541 ! Calculate the diffusive flux in y direction (from Xue's eq. 3).
6543 dflux_y_p0 = ( 10.0 * ( field(i,k,j ) - field(i,k,j-1) ) &
6544 - 5.0 * ( field(i,k,j+1) - field(i,k,j-2) ) &
6545 + ( field(i,k,j+2) - field(i,k,j-3) ) )
6547 dflux_y_p1 = ( 10.0 * ( field(i,k,j+1) - field(i,k,j ) ) &
6548 - 5.0 * ( field(i,k,j+2) - field(i,k,j-1) ) &
6549 + ( field(i,k,j+3) - field(i,k,j-2) ) )
6551 ! If requested in the namelist (diff_6th_opt=2), prohibit up-gradient diffusion
6552 ! (variation on Xue's eq. 10).
6554 IF ( diff_6th_opt .EQ. 2 ) THEN
6556 IF ( dflux_y_p0 * ( field(i,k,j )-field(i,k,j-1) ) .LE. 0.0 ) THEN
6560 IF ( dflux_y_p1 * ( field(i,k,j+1)-field(i,k,j ) ) .LE. 0.0 ) THEN
6569 IF (config_flags%diff_6th_slopeopt .ge. 1) THEN
6571 dzthresh = config_flags%diff_6th_thresh*9.81*dy
6572 IF ( name .EQ. 'u' ) THEN
6573 dzmax_p0 = max( ABS(phb(i,k,j ) - phb(i,k,j-1))*msfvy(i,j ), ABS(phb(i-1,k,j ) - phb(i-1,k,j-1))*msfvy(i-1,j ) )
6574 dzmax_p1 = max( ABS(phb(i,k,j+1) - phb(i,k,j ))*msfvy(i,j+1), ABS(phb(i-1,k,j+1) - phb(i-1,k,j ))*msfvy(i-1,j+1) )
6575 ELSE IF ( name .EQ. 'v' ) THEN
6576 dzmax_p0 = MAX( ABS(phb(i,k,j )-phb(i,k,j-1))*msfvy(i,j ), ABS(phb(i,k,j-1)-phb(i,k,j-2))*msfvy(i,j-1) )
6577 dzmax_p1 = MAX( ABS(phb(i,k,j+1)-phb(i,k,j ))*msfvy(i,j+1), ABS(phb(i,k,j )-phb(i,k,j-1))*msfvy(i,j ) )
6579 dzmax_p0 = ABS(phb(i,k,j ) - phb(i,k,j-1))*msfvy(i,j )
6580 dzmax_p1 = ABS(phb(i,k,j+1) - phb(i,k,j ))*msfvy(i,j+1)
6582 slopedamp_p0 = MAX(1.- (dzmax_p0/dzthresh) , 0.)
6583 slopedamp_p1 = MAX(1.- (dzmax_p1/dzthresh) , 0.)
6586 ! Apply 6th-order diffusion in y direction.
6588 IF ( name .EQ. 'u' ) THEN
6589 mu_avg_p0 = 0.25 * ( &
6590 (c1(k)*MUT(i-1,j-1)+c2(k)) + &
6591 (c1(k)*MUT(i ,j-1)+c2(k)) + &
6592 (c1(k)*MUT(i-1,j )+c2(k)) + &
6593 (c1(k)*MUT(i ,j )+c2(k)) )
6594 mu_avg_p1 = 0.25 * ( &
6595 (c1(k)*MUT(i-1,j )+c2(k)) + &
6596 (c1(k)*MUT(i ,j )+c2(k)) + &
6597 (c1(k)*MUT(i-1,j+1)+c2(k)) + &
6598 (c1(k)*MUT(i ,j+1)+c2(k)) )
6599 tendency_y = diff_6th_coef * msfuy(i,j) * &
6600 ( ( slopedamp_p1 * mu_avg_p1 * dflux_y_p1 ) - ( slopedamp_p0 * mu_avg_p0 * dflux_y_p0 ) )
6602 ELSE IF ( name .EQ. 'v' ) THEN
6603 mu_avg_p0 = (c1(k)*MUT(i,j-1)+c2(k))
6604 mu_avg_p1 = (c1(k)*MUT(i,j )+c2(k))
6605 tendency_y = diff_6th_coef * msfvy(i,j) * &
6606 ( ( slopedamp_p1 * mu_avg_p1 * dflux_y_p1 ) - ( slopedamp_p0 * mu_avg_p0 * dflux_y_p0 ) )
6608 mu_avg_p0 = 0.5 * ( &
6609 (c1(k)*MUT(i,j-1)+c2(k)) + &
6610 (c1(k)*MUT(i,j )+c2(k)) )
6611 mu_avg_p1 = 0.5 * ( &
6612 (c1(k)*MUT(i,j )+c2(k)) + &
6613 (c1(k)*MUT(i,j+1)+c2(k)) )
6614 tendency_y = diff_6th_coef * msfty(i,j) * &
6615 ( ( slopedamp_p1 * mu_avg_p1 * dflux_y_p1 ) - ( slopedamp_p0 * mu_avg_p0 * dflux_y_p0 ) )
6618 ! End: Diffusion in y.
6619 !------------------------------------------------------------------------------
6621 !------------------------------------------------------------------------------
6622 ! Begin: Combine diffusion in x and y.
6624 tendency(i,k,j) = tendency(i,k,j) + tendency_x + tendency_y
6626 ! End: Combine diffusion in x and y.
6627 !------------------------------------------------------------------------------
6633 ! End: Loop across spatial dimensions.
6634 !------------------------------------------------------------------------------
6636 END SUBROUTINE sixth_order_diffusion
6638 !==============================================================================
6640 SUBROUTINE initialize_moist_old ( moist_old , moist , &
6641 ids, ide, jds, jde, kds, kde , &
6642 ims, ime, jms, jme, kms, kme , &
6643 its, ite, jts, jte, kts, kte )
6645 ! For the theta_m option, the moist_old variable is uninitialized
6646 ! at the beginning of EACH of the RK steps. So, just set the
6647 ! starting value of moist_old as the final value of moist from the
6648 ! previous time step. Here "moist" is only the P_Qv index.
6652 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde , &
6653 ims, ime, jms, jme, kms, kme , &
6654 its, ite, jts, jte, kts, kte
6655 REAL , INTENT(IN ) , DIMENSION(ims:ime,kms:kme,jms:jme) :: moist
6656 REAL , INTENT( OUT) , DIMENSION(ims:ime,kms:kme,jms:jme) :: moist_old
6660 INTEGER :: i , j , k
6662 DO j = jts , MIN(jte,jde-1)
6664 DO i = its , MIN(ite,ide-1)
6665 moist_old(i,k,j) = moist(i,k,j)
6670 END SUBROUTINE initialize_moist_old
6672 !==============================================================================
6674 SUBROUTINE conv_t_tendf_to_moist ( t_1 , moist_old , &
6675 t_tendf , moist_tend , &
6676 ids, ide, jds, jde, kds, kde , &
6677 ims, ime, jms, jme, kms, kme , &
6678 its, ite, jts, jte, kts, kte )
6680 ! Convert dry potential temperature to "moist" theta:
6681 ! theta_m = theta ( 1 + Rv/Rd Qv )
6685 INTEGER , INTENT(IN) :: ids, ide, jds, jde, kds, kde , &
6686 ims, ime, jms, jme, kms, kme , &
6687 its, ite, jts, jte, kts, kte
6688 REAL , INTENT(IN ) , DIMENSION(ims:ime,kms:kme,jms:jme) :: moist_tend, moist_old, t_1
6689 REAL , INTENT(INOUT) , DIMENSION(ims:ime,kms:kme,jms:jme) :: t_tendf
6693 INTEGER :: i , j , k
6695 ! This dry tendency is from the physics packages. It is modified immediately after the
6696 ! call to the physics schemes, and the remains constant for the remainder of the RK loops.
6698 DO j = jts , MIN(jte,jde-1)
6700 DO i = its , MIN(ite,ide-1)
6701 t_tendf(i,k,j) = (1. + (R_v/R_d) * moist_old(i,k,j))*t_tendf(i,k,j) + &
6702 (R_v/R_d)*(t_1(i,k,j)+T0)/(1.+R_v/R_d*moist_old(i,k,j))*moist_tend(i,k,j)
6707 END SUBROUTINE conv_t_tendf_to_moist
6709 !==============================================================================
6711 subroutine Cloud_tracer_nudge (dt, dt_relax, dt_nudge, xtime, qc, qi, qs, &
6712 tr_qc, tr_qi, tr_qs, ids, ide, jds, jde, kds, kde, ims, ime, jms, jme, &
6713 kms, kme, its, ite, jts, jte, kts, kte)
6717 real, intent(in) :: dt, dt_relax, dt_nudge, xtime
6718 real, dimension (ims:ime, kms:kme, jms:jme), intent (inout) :: qc, qi, qs
6719 real, dimension (ims:ime, kms:kme, jms:jme), intent(in) :: tr_qc, tr_qi, tr_qs
6720 integer, intent(in) :: ids, ide, jds, jde, kds, kde, ims,ime, jms, jme, &
6721 kms, kme, its, ite, jts, jte, kts, kte
6724 integer :: i, j, k, i_start, i_end, j_start, j_end, k_start, k_end
6726 !--------------------------------------------------------------------
6730 ! Cloud_tracer_nudge nudges qc, qi, and qs towards advected tracer values
6731 ! tr_qc, tr_qi, tr_qs
6735 ! Set up loop bounds
6737 i_end = min (ite, ide - 1)
6739 j_end = min (jte, jde - 1)
6741 k_end = min (kte, kde - 1)
6744 if (xtime < dt_nudge) then
6745 do j = j_start, j_end
6746 do k = k_start, k_end
6747 do i = i_start, i_end
6748 qc(i, k, j) = qc(i, k, j) + (tr_qc(i, k, j) - qc(i, k, j)) * dt / dt_relax
6749 qi(i, k, j) = qi(i, k, j) + (tr_qi(i, k, j) - qi(i, k, j)) * dt / dt_relax
6750 qs(i, k, j) = qs(i, k, j) + (tr_qs(i, k, j) - qs(i, k, j)) * dt / dt_relax
6756 end subroutine Cloud_tracer_nudge
6758 END MODULE module_big_step_utilities_em