!WRF:MODEL_LAYER:DYNAMICS !MODULE module_big_step_utilities_em 6 USE module_domain
USE module_model_constants
USE module_state_description
USE module_configure
CONTAINS !------------------------------------------------------------------------------- !*************** new mass coordinate routines
SUBROUTINE calc_mu_uv ( mu, mub, muu, muv, & 3 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime , jms:jme ) , INTENT( OUT) :: muu, muv REAL, DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub ! local stuff INTEGER :: i, j, itf, jtf itf=ite jtf=MIN(jte,jde-1) IF ( ( its .NE. ids ) .AND. ( ite .NE. ide ) ) THEN DO j=jts,jtf DO i=its,itf muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j)) ENDDO ENDDO ELSE IF ( ( its .EQ. ids ) .AND. ( ite .NE. ide ) ) THEN DO j=jts,jtf DO i=its+1,itf muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j)) ENDDO ENDDO i=its DO j=jts,jtf muu(i,j) = mu(i,j) +mub(i,j) ENDDO ELSE IF ( ( its .NE. ids ) .AND. ( ite .EQ. ide ) ) THEN DO j=jts,jtf DO i=its,itf-1 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j)) ENDDO ENDDO i=ite DO j=jts,jtf muu(i,j) = mu(i-1,j) +mub(i-1,j) ENDDO ELSE IF ( ( its .EQ. ids ) .AND. ( ite .EQ. ide ) ) THEN DO j=jts,jtf DO i=its+1,itf-1 muu(i,j) = 0.5*(mu(i,j)+mu(i-1,j)+mub(i,j)+mub(i-1,j)) ENDDO ENDDO i=its DO j=jts,jtf muu(i,j) = mu(i,j) +mub(i,j) ENDDO i=ite DO j=jts,jtf muu(i,j) = mu(i-1,j) +mub(i-1,j) ENDDO END IF itf=MIN(ite,ide-1) jtf=jte IF ( ( jts .NE. jds ) .AND. ( jte .NE. jde ) ) THEN DO j=jts,jtf DO i=its,itf muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1)) ENDDO ENDDO ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .NE. jde ) ) THEN DO j=jts+1,jtf DO i=its,itf muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1)) ENDDO ENDDO j=jts DO i=its,itf muv(i,j) = mu(i,j) +mub(i,j) ENDDO ELSE IF ( ( jts .NE. jds ) .AND. ( jte .EQ. jde ) ) THEN DO j=jts,jtf-1 DO i=its,itf muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1)) ENDDO ENDDO j=jte DO i=its,itf muv(i,j) = mu(i,j-1) +mub(i,j-1) ENDDO ELSE IF ( ( jts .EQ. jds ) .AND. ( jte .EQ. jde ) ) THEN DO j=jts+1,jtf-1 DO i=its,itf muv(i,j) = 0.5*(mu(i,j)+mu(i,j-1)+mub(i,j)+mub(i,j-1)) ENDDO ENDDO j=jts DO i=its,itf muv(i,j) = mu(i,j) +mub(i,j) ENDDO j=jte DO i=its,itf muv(i,j) = mu(i,j-1) +mub(i,j-1) ENDDO END IF END SUBROUTINE calc_mu_uv !-------------------------------------------------------------------------------
SUBROUTINE couple ( mu, mub, rfield, field, name, & 38,2 msf, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte CHARACTER(LEN=1) , INTENT(IN ) :: name REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: rfield REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, mub, msf REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field ! Local data INTEGER :: i, j, k, itf, jtf, ktf REAL , DIMENSION(ims:ime,jms:jme) :: muu , muv ktf=MIN(kte,kde-1) IF (name .EQ. 'u')THEN CALL calc_mu_uv
( mu, mub, muu, muv, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) itf=ite jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf rfield(i,k,j)=field(i,k,j)*muu(i,j)/msf(i,j) ENDDO ENDDO ENDDO ELSE IF (name .EQ. 'v')THEN CALL calc_mu_uv
( mu, mub, muu, muv, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) itf=ite itf=MIN(ite,ide-1) jtf=jte DO j=jts,jtf DO k=kts,ktf DO i=its,itf rfield(i,k,j)=field(i,k,j)*muv(i,j)/msf(i,j) ENDDO ENDDO ENDDO ELSE IF (name .EQ. 'w')THEN itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts,kte DO i=its,itf rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j))/msf(i,j) ENDDO ENDDO ENDDO ELSE itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf rfield(i,k,j)=field(i,k,j)*(mu(i,j)+mub(i,j)) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE couple !-----------------------------------------------------------------------
SUBROUTINE calc_ww ( mu, ru, rv, ww, & rdx, rdy, msft, dnw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, rv REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu, msft REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww REAL , INTENT(IN ) :: rdx, rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf REAL , DIMENSION( its:ite ) :: dmdt jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) itf=MIN(ite,ide-1) DO j=jts,jtf DO i=its,ite dmdt(i) = 0. ww(i,1,j) = 0. ww(i,kte,j) = 0. ENDDO !! DO k=kts,ktf+1 DO k=kts,ktf DO i=its,itf dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) & +rdy*(rv(i,k,j+1)-rv(i,k,j)) ) ENDDO ENDDO ! DO K=2,NZ-1 ! ww(K,I)=ww(K-1,I)-DNW(K-1)* ! & (DMDT+RDX*( xmu(i )*u(K,I ) ! & -xmu(im1)*u(k,im1)) ) ! END DO DO k=2,ktf DO i=its,itf ww(i,k,j)=ww(i,k-1,j) & - dnw(k-1)* ( dmdt(i) & +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) & +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) ) ENDDO ENDDO ENDDO END SUBROUTINE calc_ww !-------------------------------------------------------------------------------
SUBROUTINE calc_ww_cp ( u, v, mup, mub, ww, & 1 rdx, rdy, msft, msfu, msfv, dnw, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: u, v REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mup, mub, & msft, msfu, msfv REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: ww REAL , INTENT(IN ) :: rdx, rdy ! Local data INTEGER :: i, j, k, itf, jtf, ktf REAL , DIMENSION( its:ite ) :: dmdt REAL , DIMENSION( its:ite, kts:kte ) :: divv REAL , DIMENSION( its:ite+1, jts:jte+1 ) :: muu, muv jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) itf=MIN(ite,ide-1) ! mu coupled with the appropriate map factor DO j=jts,jtf DO i=its,min(ite+1,ide) muu(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i-1,j)+mub(i-1,j))/msfu(i,j) ENDDO ENDDO DO j=jts,min(jte+1,jde) DO i=its,itf muv(i,j) = 0.5*(mup(i,j)+mub(i,j)+mup(i,j-1)+mub(i,j-1))/msfv(i,j) ENDDO ENDDO DO j=jts,jtf DO i=its,ite dmdt(i) = 0. ww(i,1,j) = 0. ww(i,kte,j) = 0. ENDDO DO k=kts,ktf DO i=its,itf divv(i,k) = msft(i,j)*dnw(k)*( rdx*(muu(i+1,j)*u(i+1,k,j)-muu(i,j)*u(i,k,j)) & +rdy*(muv(i,j+1)*v(i,k,j+1)-muv(i,j)*v(i,k,j)) ) ! dmdt(i) = dmdt(i) + dnw(k)* ( rdx*(ru(i+1,k,j)-ru(i,k,j)) & ! +rdy*(rv(i,k,j+1)-rv(i,k,j)) ) dmdt(i) = dmdt(i) + divv(i,k) ENDDO ENDDO DO k=2,ktf DO i=its,itf ! ww(i,k,j)=ww(i,k-1,j) & ! - dnw(k-1)* ( dmdt(i) & ! +rdx*(ru(i+1,k-1,j)-ru(i,k-1,j)) & ! +rdy*(rv(i,k-1,j+1)-rv(i,k-1,j)) ) ww(i,k,j)=ww(i,k-1,j) - dnw(k-1)*dmdt(i) - divv(i,k-1) ENDDO ENDDO ENDDO END SUBROUTINE calc_ww_cp !-------------------------------------------------------------------------------
SUBROUTINE calc_cq ( moist, cqu, cqv, cqw, n_moist, & 1 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN ) :: n_moist REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN ) :: moist REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT( OUT) :: cqu, cqv, cqw ! Local stuff REAL :: qtot INTEGER :: i, j, k, itf, jtf, ktf, ispe itf=ite jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) IF( n_moist >= PARAM_FIRST_SCALAR ) THEN DO j=jts,jtf DO k=kts,ktf DO i=its,itf qtot = 0. DO ispe=PARAM_FIRST_SCALAR,n_moist qtot = qtot + moist(i,k,j,ispe) + moist(i-1,k,j,ispe) ENDDO ! qtot = 0.5*( moist(i ,k,j,1)+moist(i ,k,j,2)+moist(i ,k,j,3)+ & ! & moist(i-1,k,j,1)+moist(i-1,k,j,2)+moist(i-1,k,j,3) ) ! cqu(i,k,j) = 1./(1.+qtot) cqu(i,k,j) = 1./(1.+0.5*qtot) ENDDO ENDDO ENDDO itf=MIN(ite,ide-1) jtf=jte DO j=jts,jtf DO k=kts,ktf DO i=its,itf qtot = 0. DO ispe=PARAM_FIRST_SCALAR,n_moist qtot = qtot + moist(i,k,j,ispe) + moist(i,k,j-1,ispe) ENDDO ! qtot = 0.5*( moist(i,k,j ,1)+moist(i,k,j ,2)+moist(i,k,j ,3)+ & ! & moist(i,k,j-1,1)+moist(i,k,j-1,2)+moist(i,k,j-1,3) ) ! cqv(i,k,j) = 1./(1.+qtot) cqv(i,k,j) = 1./(1.+0.5*qtot) ENDDO ENDDO ENDDO itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts+1,ktf DO i=its,itf qtot = 0. DO ispe=PARAM_FIRST_SCALAR,n_moist qtot = qtot + moist(i,k,j,ispe) + moist(i,k-1,j,ispe) ENDDO ! qtot = 0.5*( moist(i,k ,j,1)+moist(i,k ,j,2)+moist(i,k-1,j,3)+ & ! & moist(i,k-1,j,1)+moist(i,k-1,j,2)+moist(i,k ,j,3) ) ! cqw(i,k,j) = qtot cqw(i,k,j) = 0.5*qtot ENDDO ENDDO ENDDO ELSE DO j=jts,jtf DO k=kts,ktf DO i=its,itf cqu(i,k,j) = 1. ENDDO ENDDO ENDDO itf=MIN(ite,ide-1) jtf=jte DO j=jts,jtf DO k=kts,ktf DO i=its,itf cqv(i,k,j) = 1. ENDDO ENDDO ENDDO itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts+1,ktf DO i=its,itf cqw(i,k,j) = 0. ENDDO ENDDO ENDDO END IF END SUBROUTINE calc_cq !----------------------------------------------------------------------
SUBROUTINE calc_alt ( alt, al, alb, & 1 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, al REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: alt ! Local stuff INTEGER :: i, j, k, itf, jtf, ktf itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf alt(i,k,j) = al(i,k,j)+alb(i,k,j) ENDDO ENDDO ENDDO END SUBROUTINE calc_alt !----------------------------------------------------------------------
SUBROUTINE calc_p_rho_phi ( moist, n_moist, & 2 al, alb, mu, muts, ph, p, pb, & t, p0, t0, znu, dnw, rdnw, & rdn, non_hydrostatic, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data LOGICAL , INTENT(IN ) :: non_hydrostatic INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN ) :: n_moist REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(IN ) :: alb, & pb, & t REAL, DIMENSION( ims:ime , kms:kme , jms:jme, n_moist ), INTENT(IN ) :: moist REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT( OUT) :: al, p REAL, DIMENSION( ims:ime , kms:kme , jms:jme ), INTENT(INOUT) :: ph REAL, DIMENSION( ims:ime , jms:jme ), INTENT(IN ) :: mu, muts REAL, DIMENSION( kms:kme ), INTENT(IN ) :: znu, dnw, rdnw, rdn REAL, INTENT(IN ) :: t0, p0 ! Local stuff INTEGER :: i, j, k, itf, jtf, ktf, ispe REAL :: qvf, qtot, qf1, qf2 itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) IF (non_hydrostatic) THEN IF (n_moist >= PARAM_FIRST_SCALAR ) THEN DO j=jts,jtf DO k=kts,ktf DO i=its,itf qvf = 1.+rvovrd*moist(i,k,j,P_QV) al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j) & +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j))) p(i,k,j)= ( (r_d*(t0+t(i,k,j))*qvf)/ & (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv & *p0 -pb(i,k,j) ENDDO ENDDO ENDDO ELSE DO j=jts,jtf DO k=kts,ktf DO i=its,itf al(i,k,j)=-1./muts(i,j)*(alb(i,k,j)*mu(i,j) & +rdnw(k)*(ph(i,k+1,j)-ph(i,k,j))) p(i,k,j)=p0*( (r_d*(t0+t(i,k,j)))/ & (p0*(al(i,k,j)+alb(i,k,j))) )**cpovcv & -pb(i,k,j) ENDDO ENDDO ENDDO END IF ELSE ! hydrostatic pressure, al, and ph1 calc; WCS, 5 sept 2001 IF (n_moist >= PARAM_FIRST_SCALAR ) THEN DO j=jts,jtf k=ktf ! top layer DO i=its,itf qtot = 0. DO ispe=PARAM_FIRST_SCALAR,n_moist qtot = qtot + moist(i,k,j,ispe) ENDDO qf2 = 1./(1.+qtot) qf1 = qtot*qf2 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2 qvf = 1.+rvovrd*moist(i,k,j,P_QV) al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* & (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j) ENDDO DO k=ktf-1,kts,-1 ! remaining layers, integrate down DO i=its,itf qtot = 0. DO ispe=PARAM_FIRST_SCALAR,n_moist qtot = qtot + 0.5*( moist(i,k ,j,ispe) + moist(i,k+1,j,ispe) ) ENDDO qf2 = 1./(1.+qtot) qf1 = qtot*qf2 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1) qvf = 1.+rvovrd*moist(i,k,j,P_QV) al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* & (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j) ENDDO ENDDO DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential DO i=its,itf ! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( & ! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ & ! mu(i,j)*alb(i,k-1,j) ) ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( & (muts(i,j))*al(i,k-1,j)+ & mu(i,j)*alb(i,k-1,j) ) ENDDO ENDDO ENDDO ELSE DO j=jts,jtf k=ktf ! top layer DO i=its,itf qtot = 0. qf2 = 1./(1.+qtot) qf1 = qtot*qf2 p(i,k,j) = - 0.5*(mu(i,j)+qf1*muts(i,j))/rdnw(k)/qf2 qvf = 1. al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* & (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j) ENDDO DO k=ktf-1,kts,-1 ! remaining layers, integrate down DO i=its,itf qtot = 0. qf2 = 1./(1.+qtot) qf1 = qtot*qf2 p(i,k,j) = p(i,k+1,j) - (mu(i,j) + qf1*muts(i,j))/qf2/rdn(k+1) qvf = 1. al(i,k,j) = (r_d/p1000mb)*(t(i,k,j)+t0)*qvf* & (((p(i,k,j)+pb(i,k,j))/p1000mb)**cvpm) - alb(i,k,j) ENDDO ENDDO DO k=2,ktf+1 ! integrate hydrostatic equation for geopotential DO i=its,itf ! ph(i,k,j) = ph(i,k-1,j) - (1./rdnw(k-1))*( & ! (muts(i,j)+mu(i,j))*al(i,k-1,j)+ & ! mu(i,j)*alb(i,k-1,j) ) ph(i,k,j) = ph(i,k-1,j) - (dnw(k-1))*( & (muts(i,j))*al(i,k-1,j)+ & mu(i,j)*alb(i,k-1,j) ) ENDDO ENDDO ENDDO END IF END IF END SUBROUTINE calc_p_rho_phi !----------------------------------------------------------------------
SUBROUTINE calc_php ( php, ph, phb, & 1 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: phb, ph REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: php ! Local stuff INTEGER :: i, j, k, itf, jtf, ktf itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf php(i,k,j) = 0.5*(phb(i,k,j)+phb(i,k+1,j)+ph(i,k,j)+ph(i,k+1,j)) ENDDO ENDDO ENDDO END SUBROUTINE calc_php !-------------------------------------------------------------------------------
SUBROUTINE diagnose_w( ph_tend, ph_new, ph_old, w, mu, dt, & 2 u, v, ht, & cf1, cf2, cf3, rdx, rdy, msft, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ph_tend, & ph_new, & ph_old, & u, & v REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: w REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mu, ht, msft REAL, INTENT(IN ) :: dt, cf1, cf2, cf3, rdx, rdy INTEGER :: i, j, k, itf, jtf itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j = jts, jtf ! lower b.c. on w DO i = its, itf w(i,1,j)= msft(i,j)*( & .5*rdy*( & (ht(i,j+1)-ht(i,j )) & *(cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1)) & +(ht(i,j )-ht(i,j-1)) & *(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )) ) & +.5*rdx*( & (ht(i+1,j)-ht(i,j )) & *(cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j)) & +(ht(i,j )-ht(i,j-1)) & *(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)) ) & ) ENDDO ! use geopotential equation to diagnose w DO k = 2, kte DO i = its, itf w(i,k,j) = msft(i,j)*( (ph_new(i,k,j)-ph_old(i,k,j))/dt & - ph_tend(i,k,j)/mu(i,j) )/g ENDDO ENDDO ENDDO END SUBROUTINE diagnose_w !-------------------------------------------------------------------------------
SUBROUTINE rhs_ph( ph_tend, u, v, ww, & 2 ph, ph_old, phb, w, & mut, muu, muv, & fnm, fnp, & rdnw, cfn, cfn1, rdx, rdy, msft, & non_hydrostatic, & config_flags, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: & u, & v, & ww, & ph, & ph_old, & phb, & w REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT( OUT) :: ph_tend REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mut, msft REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp REAL, INTENT(IN ) :: cfn, cfn1, rdx, rdy LOGICAL, INTENT(IN ) :: non_hydrostatic ! Local stuff INTEGER :: i, j, k, itf, jtf, ktf, kz, i_start, j_start REAL :: ur, ul, ub, vr, vl, vb REAL, DIMENSION(its:ite,kts:kte) :: wdwn INTEGER :: advective_order LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. advective_order = config_flags%h_sca_adv_order ! advective_order = 2 ! original configuration (pre Oct 2001) itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) ! advective form for the geopotential equation DO j = jts, jtf DO k = 2, kte DO i = its, itf wdwn(i,k) = .5*(ww(i,k,j)+ww(i,k-1,j))*rdnw(k-1) & *(ph(i,k,j)-ph(i,k-1,j)+phb(i,k,j)-phb(i,k-1,j)) ENDDO ENDDO DO k = 2, kte-1 DO i = its, itf ph_tend(i,k,j) = ph_tend(i,k,j) & - (fnm(k)*wdwn(i,k+1)+fnp(k)*wdwn(i,k)) ENDDO ENDDO ENDDO IF (non_hydrostatic) THEN ! add in "gw" term. DO j = jts, jtf ! in hydrostatic mode, "gw" will be diagnosed ! after the timestep to give us "w" DO i = its, itf ph_tend(i,kde,j) = 0. ENDDO DO k = 2, kte DO i = its, itf ph_tend(i,k,j) = ph_tend(i,k,j) + mut(i,j)*g*w(i,k,j)/msft(i,j) ENDDO ENDDO ENDDO END IF IF (advective_order <= 2) THEN ! y (v) advection i_start = its j_start = jts itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1 IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1 DO j = j_start, jtf DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdy* & ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1))* & (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) & +muv(i,j )*(v(i,k,j )+v(i,k-1,j ))* & (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdy* & ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1))* & (phb(i,k,j+1)-phb(i,k,j )+ph(i,k,j+1)-ph(i,k,j )) & +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j ))* & (phb(i,k,j )-phb(i,k,j-1)+ph(i,k,j )-ph(i,k,j-1)) ) ENDDO ENDDO ! x (u) advection i_start = its j_start = jts itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1 DO j = j_start, jtf DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdx* & ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j))* & (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) & +muu(i ,j)*(u(i ,k,j)+u(i ,k-1,j))* & (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdx* & ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j))* & (phb(i+1,k,j)-phb(i ,k,j)+ph(i+1,k,j)-ph(i ,k,j)) & +muu(i ,j)*(cfn*u(i ,k-1,j)+cfn1*u(i ,k-2,j))* & (phb(i ,k,j)-phb(i-1,k,j)+ph(i ,k,j)-ph(i-1,k,j)) ) ENDDO ENDDO ELSE IF (advective_order <= 4) THEN ! y (v) advection i_start = its j_start = jts itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1 IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1 DO j = j_start, jtf DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdy* ( & ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) & +muv(i,j )*(v(i,k,j )+v(i,k-1,j )) )* (1./12.)*( & 8.*(ph(i,k,j+1)-ph(i,k,j-1)) & -(ph(i,k,j+2)-ph(i,k,j-2)) & +8.*(phb(i,k,j+1)-phb(i,k,j-1)) & -(phb(i,k,j+2)-phb(i,k,j-2)) ) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdy* ( & ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) & +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j )) )* (1./12.)*( & 8.*(ph(i,k,j+1)-ph(i,k,j-1)) & -(ph(i,k,j+2)-ph(i,k,j-2)) & +8.*(phb(i,k,j+1)-phb(i,k,j-1)) & -(phb(i,k,j+2)-phb(i,k,j-2)) ) ) ENDDO ENDDO ! x (u) advection i_start = its j_start = jts itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) IF ( (config_flags%open_xs) .and. its == ids ) i_start = its+1 IF ( (config_flags%open_xe) .and. ite == ide ) itf = itf-1 DO j = j_start, jtf DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdx*( & ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) & +muu(i,j )*(u(i,k,j )+u(i,k-1,j )) )* (1./12.)*( & 8.*(ph(i+1,k,j)-ph(i-1,k,j)) & -(ph(i+2,k,j)-ph(i-2,k,j)) & +8.*(phb(i+1,k,j)-phb(i-1,k,j)) & -(phb(i+2,k,j)-phb(i-2,k,j)) ) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdx*( & ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) & +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j)) )* (1./12.)*( & 8.*(ph(i+1,k,j)-ph(i-1,k,j)) & -(ph(i+2,k,j)-ph(i-2,k,j)) & +8.*(phb(i+1,k,j)-phb(i-1,k,j)) & -(phb(i+2,k,j)-phb(i-2,k,j)) ) ) ENDDO ENDDO ELSE IF (advective_order <= 6) THEN ! y (v) advection i_start = its j_start = jts itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ! IF ( (config_flags%open_ys) .and. jts == jds ) j_start = jts+1 ! IF ( (config_flags%open_ye) .and. jte == jde ) jtf = jtf-1 IF (config_flags%open_ys .or. specified ) j_start = max(jts,jds+2) IF (config_flags%open_ye .or. specified ) jtf = min(jtf,jde-3) DO j = j_start, jtf DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdy* ( & ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) & +muv(i,j )*(v(i,k,j )+v(i,k-1,j )) )* (1./60.)*( & 45.*(ph(i,k,j+1)-ph(i,k,j-1)) & -9.*(ph(i,k,j+2)-ph(i,k,j-2)) & +(ph(i,k,j+3)-ph(i,k,j-3)) & +45.*(phb(i,k,j+1)-phb(i,k,j-1)) & -9.*(phb(i,k,j+2)-phb(i,k,j-2)) & +(phb(i,k,j+3)-phb(i,k,j-3)) ) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdy* ( & ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) & +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j )) )* (1./60.)*( & 45.*(ph(i,k,j+1)-ph(i,k,j-1)) & -9.*(ph(i,k,j+2)-ph(i,k,j-2)) & +(ph(i,k,j+3)-ph(i,k,j-3)) & +45.*(phb(i,k,j+1)-phb(i,k,j-1)) & -9.*(phb(i,k,j+2)-phb(i,k,j-2)) & +(phb(i,k,j+3)-phb(i,k,j-3)) ) ) ENDDO ENDDO ! pick up near boundary rows using 4th order stencil ! (open bc copy only goes out to jds-1 and jde, hence 4rth is ok but 6th is too big) IF ( (config_flags%open_ys) .and. jts <= jds+1 ) THEN j = jds+1 DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdy* ( & ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) & +muv(i,j )*(v(i,k,j )+v(i,k-1,j )) )* (1./12.)*( & 8.*(ph(i,k,j+1)-ph(i,k,j-1)) & -(ph(i,k,j+2)-ph(i,k,j-2)) & +8.*(phb(i,k,j+1)-phb(i,k,j-1)) & -(phb(i,k,j+2)-phb(i,k,j-2)) ) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdy* ( & ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) & +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j )) )* (1./12.)*( & 8.*(ph(i,k,j+1)-ph(i,k,j-1)) & -(ph(i,k,j+2)-ph(i,k,j-2)) & +8.*(phb(i,k,j+1)-phb(i,k,j-1)) & -(phb(i,k,j+2)-phb(i,k,j-2)) ) ) ENDDO END IF IF ( (config_flags%open_ye) .and. jte >= jde-2 ) THEN j = jde-2 DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdy* ( & ( muv(i,j+1)*(v(i,k,j+1)+v(i,k-1,j+1)) & +muv(i,j )*(v(i,k,j )+v(i,k-1,j )) )* (1./12.)*( & 8.*(ph(i,k,j+1)-ph(i,k,j-1)) & -(ph(i,k,j+2)-ph(i,k,j-2)) & +8.*(phb(i,k,j+1)-phb(i,k,j-1)) & -(phb(i,k,j+2)-phb(i,k,j-2)) ) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdy* ( & ( muv(i,j+1)*(cfn*v(i,k-1,j+1)+cfn1*v(i,k-2,j+1)) & +muv(i,j )*(cfn*v(i,k-1,j )+cfn1*v(i,k-2,j )) )* (1./12.)*( & 8.*(ph(i,k,j+1)-ph(i,k,j-1)) & -(ph(i,k,j+2)-ph(i,k,j-2)) & +8.*(phb(i,k,j+1)-phb(i,k,j-1)) & -(phb(i,k,j+2)-phb(i,k,j-2)) ) ) ENDDO END IF ! x (u) advection i_start = its j_start = jts itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) IF (config_flags%open_xs .or. specified ) i_start = max(its,ids+2) IF (config_flags%open_xe .or. specified ) itf = min(itf,ide-3) DO j = j_start, jtf DO k = 2, kte-1 DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdx*( & ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) & +muu(i,j )*(u(i,k,j )+u(i,k-1,j )) )* (1./60.)*( & 45.*(ph(i+1,k,j)-ph(i-1,k,j)) & -9.*(ph(i+2,k,j)-ph(i-2,k,j)) & +(ph(i+3,k,j)-ph(i-3,k,j)) & +45.*(phb(i+1,k,j)-phb(i-1,k,j)) & -9.*(phb(i+2,k,j)-phb(i-2,k,j)) & +(phb(i+3,k,j)-phb(i-3,k,j)) ) ) ENDDO ENDDO k = kte DO i = i_start, itf ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdx*( & ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) & +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j)) )* (1./60.)*( & 45.*(ph(i+1,k,j)-ph(i-1,k,j)) & -9.*(ph(i+2,k,j)-ph(i-2,k,j)) & +(ph(i+3,k,j)-ph(i-3,k,j)) & +45.*(phb(i+1,k,j)-phb(i-1,k,j)) & -9.*(phb(i+2,k,j)-phb(i-2,k,j)) & +(phb(i+3,k,j)-phb(i-3,k,j)) ) ) ENDDO ENDDO IF ( (config_flags%open_xs) .and. its <= ids+1 ) THEN i = ids + 1 DO j = j_start, jtf DO k = 2, kte-1 ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdx*( & ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) & +muu(i,j )*(u(i,k,j )+u(i,k-1,j )) )* (1./12.)*( & 8.*(ph(i+1,k,j)-ph(i-1,k,j)) & -(ph(i+2,k,j)-ph(i-2,k,j)) & +8.*(phb(i+1,k,j)-phb(i-1,k,j)) & -(phb(i+2,k,j)-phb(i-2,k,j)) ) ) ENDDO k = kte ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdx*( & ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) & +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j)) )* (1./12.)*( & 8.*(ph(i+1,k,j)-ph(i-1,k,j)) & -(ph(i+2,k,j)-ph(i-2,k,j)) & +8.*(phb(i+1,k,j)-phb(i-1,k,j)) & -(phb(i+2,k,j)-phb(i-2,k,j)) ) ) ENDDO END IF IF ( (config_flags%open_xe) .and. ite >= ide-2 ) THEN i = ide-2 DO j = j_start, jtf DO k = 2, kte-1 ph_tend(i,k,j)=ph_tend(i,k,j) - .25*rdx*( & ( muu(i+1,j)*(u(i+1,k,j)+u(i+1,k-1,j)) & +muu(i,j )*(u(i,k,j )+u(i,k-1,j )) )* (1./12.)*( & 8.*(ph(i+1,k,j)-ph(i-1,k,j)) & -(ph(i+2,k,j)-ph(i-2,k,j)) & +8.*(phb(i+1,k,j)-phb(i-1,k,j)) & -(phb(i+2,k,j)-phb(i-2,k,j)) ) ) ENDDO k = kte ph_tend(i,k,j)=ph_tend(i,k,j) - .5*rdx*( & ( muu(i+1,j)*(cfn*u(i+1,k-1,j)+cfn1*u(i+1,k-2,j)) & +muu(i,j )*(cfn*u(i ,k-1,j)+cfn1*u(i,k-2,j)) )* (1./12.)*( & 8.*(ph(i+1,k,j)-ph(i-1,k,j)) & -(ph(i+2,k,j)-ph(i-2,k,j)) & +8.*(phb(i+1,k,j)-phb(i-1,k,j)) & -(phb(i+2,k,j)-phb(i-2,k,j)) ) ) ENDDO END IF END IF ! lateral open boundary conditions, ! start with north and south (y) boundaries i_start = its itf=MIN(ite,ide-1) ! south IF ( (config_flags%open_ys) .and. jts == jds ) THEN j=jts DO k=2,kde kz = min(k,kde-1) DO i = its,itf vb =.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j )) & +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j )) ) vl=amin1(vb,0.) ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( & +vl*(ph_old(i,k,j+1)-ph_old(i,k,j))) ENDDO ENDDO END IF ! north IF ( (config_flags%open_ye) .and. jte == jde ) THEN j=jte-1 DO k=2,kde kz = min(k,kde-1) DO i = its,itf vb=.5*( fnm(kz)*(v(i,kz ,j+1)+v(i,kz ,j)) & +fnp(kz)*(v(i,kz-1,j+1)+v(i,kz-1,j)) ) vr=amax1(vb,0.) ph_tend(i,k,j)=ph_tend(i,k,j)-rdy*mut(i,j)*( & +vr*(ph_old(i,k,j)-ph_old(i,k,j-1))) ENDDO ENDDO END IF ! now the east and west (y) boundaries j_start = its jtf=MIN(jte,jde-1) ! west IF ( (config_flags%open_xs) .and. its == ids ) THEN i=its DO j = jts,jtf DO k=2,kde-1 kz = k ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) & +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) ) ul=amin1(ub,0.) ph_tend(i,k,j)=ph_tend(i,k,j)-rdx*mut(i,j)*( & +ul*(ph_old(i+1,k,j)-ph_old(i,k,j))) ENDDO k = kde kz = k ub =.5*( fnm(kz)*(u(i+1,kz ,j)+u(i ,kz ,j)) & +fnp(kz)*(u(i+1,kz-1,j)+u(i ,kz-1,j)) ) ul=amin1(ub,0.) ph_tend(i,k,j)=ph_tend(i,k,j)-rdx*mut(i,j)*( & +ul*(ph_old(i+1,k,j)-ph_old(i,k,j))) ENDDO END IF ! east IF ( (config_flags%open_xe) .and. ite == ide ) THEN i = ite-1 DO j = jts,jtf DO k=2,kde-1 kz = k ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) & +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) ) ur=amax1(ub,0.) ph_tend(i,k,j)=ph_tend(i,k,j)-rdx*mut(i,j)*( & +ur*(ph_old(i,k,j)-ph_old(i-1,k,j))) ENDDO k = kde kz = k-1 ub=.5*( fnm(kz)*(u(i+1,kz ,j)+u(i,kz ,j)) & +fnp(kz)*(u(i+1,kz-1,j)+u(i,kz-1,j)) ) ur=amax1(ub,0.) ph_tend(i,k,j)=ph_tend(i,k,j)-rdx*mut(i,j)*( & +ur*(ph_old(i,k,j)-ph_old(i-1,k,j))) ENDDO END IF END SUBROUTINE rhs_ph !-------------------------------------------------------------------------------
SUBROUTINE horizontal_pressure_gradient( ru_tend,rv_tend, & 1 ph,alt,p,pb,al,php,cqu,cqv, & muu,muv,mu,fnm,fnp,rdnw, & cf1,cf2,cf3,rdx,rdy,msft, & config_flags, non_hydrostatic, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags LOGICAL, INTENT (IN ) :: non_hydrostatic INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: & ph, & alt, & al, & p, & pb, & php, & cqu, & cqv REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: & ru_tend, & rv_tend REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: muu, muv, mu, msft REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, fnm, fnp REAL, INTENT(IN ) :: rdx, rdy, cf1, cf2, cf3 INTEGER :: i,j,k, itf, jtf, ktf, i_start, j_start REAL, DIMENSION( ims:ime, kms:kme ) :: dpn REAL :: dpx, dpy LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ! start with the north-south (y) pressure gradient itf=MIN(ite,ide-1) jtf=jte ktf=MIN(kte,kde-1) i_start = its j_start = jts IF ( (config_flags%open_ys .or. specified .or. & config_flags%nested ) .and. jts == jds ) j_start = jts+1 IF ( (config_flags%open_ye .or. specified .or. & config_flags%nested ) .and. jte == jde ) jtf = jtf-1 DO j = j_start, jtf IF ( non_hydrostatic ) THEN k=1 DO i = i_start, itf dpn(i,k) = .5*( cf1*(p(i,k ,j-1)+p(i,k ,j)) & +cf2*(p(i,k+1,j-1)+p(i,k+1,j)) & +cf3*(p(i,k+2,j-1)+p(i,k+2,j)) ) dpn(i,kde) = 0. ENDDO DO k=2,ktf DO i = i_start, itf dpn(i,k) = .5*( fnm(k)*(p(i,k ,j-1)+p(i,k ,j)) & +fnp(k)*(p(i,k-1,j-1)+p(i,k-1,j)) ) END DO END DO ELSE DO k=1,kde DO i = i_start, itf dpn(i,k) = 0. END DO END DO END IF DO K=1,ktf DO i = i_start, itf dpy = .5*rdy*muv(i,j)*( & (ph (i,k+1,j)-ph (i,k+1,j-1) + ph(i,k,j)-ph(i,k,j-1)) & +(alt(i,k ,j)+alt(i,k ,j-1))*(p (i,k,j)-p (i,k,j-1)) & +(al (i,k ,j)+al (i,k ,j-1))*(pb(i,k,j)-pb(i,k,j-1)) ) dpy = dpy + rdy*(php(i,k,j)-php(i,k,j-1))* & (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i,j-1)+mu(i,j))) rv_tend(i,k,j) = rv_tend(i,k,j)-cqv(i,k,j)*dpy END DO END DO ENDDO ! now the east-west (x) pressure gradient itf=ite jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) i_start = its j_start = jts IF ( (config_flags%open_xs .or. specified .or. & config_flags%nested ) .and. its == ids ) i_start = its+1 IF ( (config_flags%open_xe .or. specified .or. & config_flags%nested ) .and. ite == ide ) itf = itf-1 DO j = j_start, jtf IF ( non_hydrostatic ) THEN k=1 DO i = i_start, itf dpn(i,k) = .5*( cf1*(p(i-1,k ,j)+p(i,k ,j)) & +cf2*(p(i-1,k+1,j)+p(i,k+1,j)) & +cf3*(p(i-1,k+2,j)+p(i,k+2,j)) ) dpn(i,kde) = 0. ENDDO DO k=2,ktf DO i = i_start, itf dpn(i,k) = .5*( fnm(k)*(p(i-1,k ,j)+p(i,k ,j)) & +fnp(k)*(p(i-1,k-1,j)+p(i,k-1,j)) ) END DO END DO ELSE DO k=1,kde DO i = i_start, itf dpn(i,k) = 0. END DO END DO END IF DO K=1,ktf DO i = i_start, itf dpx = .5*rdx*muu(i,j)*( & (ph (i,k+1,j)-ph (i-1,k+1,j) + ph(i,k,j)-ph(i-1,k,j)) & +(alt(i,k ,j)+alt(i-1,k ,j))*(p (i,k,j)-p (i-1,k,j)) & +(al (i,k ,j)+al (i-1,k ,j))*(pb(i,k,j)-pb(i-1,k,j)) ) dpx = dpx + rdx*(php(i,k,j)-php(i-1,k,j))* & (rdnw(k)*(dpn(i,k+1)-dpn(i,k))-.5*(mu(i-1,j)+mu(i,j))) ru_tend(i,k,j) = ru_tend(i,k,j)-cqu(i,k,j)*dpx END DO END DO ENDDO END SUBROUTINE horizontal_pressure_gradient !-------------------------------------------------------------------------------
SUBROUTINE pg_buoy_w( rw_tend, p, cqw, mu, mub, & 1 rdnw, rdn, g, msft, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: p REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: cqw REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mub, mu, msft REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw, rdn REAL, INTENT(IN ) :: g INTEGER :: itf, jtf, i, j, k REAL :: cq1, cq2 ! BUOYANCY AND PRESSURE GRADIENT TERM IN W EQUATION AT TIME T itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j = jts,jtf k=kde DO i=its,itf cq1 = 1./(1.+cqw(i,k-1,j)) cq2 = cqw(i,k-1,j)*cq1 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msft(i,j))*g*( & cq1*2.*rdnw(k-1)*( -p(i,k-1,j)) & -mu(i,j)-cq2*mub(i,j) ) END DO DO k = 2, kde-1 DO i = its,itf cq1 = 1./(1.+cqw(i,k,j)) cq2 = cqw(i,k,j)*cq1 cqw(i,k,j) = cq1 rw_tend(i,k,j) = rw_tend(i,k,j)+(1./msft(i,j))*g*( & cq1*rdn(k)*(p(i,k,j)-p(i,k-1,j)) & -mu(i,j)-cq2*mub(i,j) ) END DO ENDDO ENDDO END SUBROUTINE pg_buoy_w !-------------------------------------------------------------------------------
SUBROUTINE w_damp( rw_tend, ww, w, mut, rdnw, dt, & 1 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(IN ) :: ww, w REAL, DIMENSION( ims:ime, kms:kme , jms:jme ), INTENT(INOUT) :: rw_tend REAL, DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: mut REAL, DIMENSION( kms:kme ), INTENT(IN ) :: rdnw REAL, INTENT(IN) :: dt REAL :: cfl INTEGER :: itf, jtf, i, j, k itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j = jts,jtf DO k = 2, kde-1 DO i = its,itf cfl = abs(ww(i,k,j)/mut(i,j)*rdnw(k)*dt) if(cfl .gt. w_beta)then print *,i,j,k,'cfl,w,d(eta)=',cfl,w(i,k,j),-1./rdnw(k) rw_tend(i,k,j) = rw_tend(i,k,j)-sign(1.,w(i,k,j))*w_alpha*(cfl-w_beta)*mut(i,j) endif END DO ENDDO ENDDO END SUBROUTINE w_damp !-------------------------------------------------------------------------------
SUBROUTINE horizontal_diffusion ( name, field, tendency, mu, & 8 config_flags, & msfu, msfv, msft, khdif, xkmhd, rdx, rdy, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte CHARACTER(LEN=1) , INTENT(IN ) :: name REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, xkmhd REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , INTENT(IN ) :: rdx, & rdy, & khdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end REAL :: mrdx, mkrdxm, mkrdxp, & mrdy, mkrdym, mkrdyp, & rcoup REAL :: pr = 3. LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) IF (name .EQ. 'u') THEN i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mkrdxm=msft(i-1,j)*xkmhd(i-1,k,j)*rdx mkrdxp=msft(i,j)*xkmhd(i,k,j)*rdx mrdx=msfu(i,j)*rdx mkrdym=0.5*(msfu(i,j)+msfu(i,j-1))* & 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdy mkrdyp=0.5*(msfu(i,j)+msfu(i,j+1))* & 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j+1)+xkmhd(i-1,k,j+1)+xkmhd(i-1,k,j))*rdy mrdy=msfu(i,j)*rdy rcoup=0.5*(mu(i,j)+mu(i-1,j)) tendency(i,k,j)=tendency(i,k,j)+rcoup*( & mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) & -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) & +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) & -mkrdym*(field(i,k,j )-field(i,k,j-1)))) ENDDO ENDDO ENDDO ELSE IF (name .EQ. 'v')THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = jte IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mkrdxm=0.5*(msfv(i,j)+msfv(i-1,j))* & 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i-1,k,j-1)+xkmhd(i-1,k,j))*rdx mkrdxp=0.5*(msfv(i,j)+msfv(i+1,j))* & 0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i+1,k,j-1)+xkmhd(i+1,k,j))*rdx mrdx=msfv(i,j)*rdx mkrdym=msft(i,j-1)*xkmhd(i,k,j-1)*rdy mkrdyp=msft(i,j)*xkmhd(i,k,j)*rdy mrdy=msfv(i,j)*rdy rcoup=0.5*(mu(i,j)+mu(i,j-1)) tendency(i,k,j)=tendency(i,k,j)+rcoup*( & mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) & -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) & +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) & -mkrdym*(field(i,k,j )-field(i,k,j-1)))) ENDDO ENDDO ENDDO ELSE IF (name .EQ. 'w')THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts+1,ktf DO i = i_start, i_end mkrdxm=msfu(i,j)*0.25*(xkmhd(i,k,j)+xkmhd(i-1,k,j)+xkmhd(i,k-1,j)+xkmhd(i-1,k-1,j))*rdx mkrdxp=msfu(i+1,j)*0.25*(xkmhd(i+1,k,j)+xkmhd(i,k,j)+xkmhd(i+1,k-1,j)+xkmhd(i,k-1,j))*rdx mrdx=msft(i,j)*rdx mkrdym=msfv(i,j)*0.25*(xkmhd(i,k,j)+xkmhd(i,k,j-1)+xkmhd(i,k-1,j)+xkmhd(i,k-1,j-1))*rdy mkrdyp=msfv(i,j+1)*0.25*(xkmhd(i,k,j+1)+xkmhd(i,k,j)+xkmhd(i,k-1,j+1)+xkmhd(i,k-1,j))*rdy mrdy=msft(i,j)*rdy rcoup=0.5*(mu(i,j)+mu(i,j)) tendency(i,k,j)=tendency(i,k,j)+rcoup*( & mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) & -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) & +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) & -mkrdym*(field(i,k,j )-field(i,k,j-1)))) ENDDO ENDDO ENDDO ELSE i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mkrdxm=msfu(i,j)*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*rdx*pr mkrdxp=msfu(i+1,j)*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*rdx*pr mrdx=msft(i,j)*rdx mkrdym=msfv(i,j)*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*rdy*pr mkrdyp=msfv(i,j+1)*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*rdy*pr mrdy=msft(i,j)*rdy rcoup=mu(i,j) tendency(i,k,j)=tendency(i,k,j)+rcoup*( & mrdx*(mkrdxp*(field(i+1,k,j)-field(i ,k,j)) & -mkrdxm*(field(i ,k,j)-field(i-1,k,j))) & +mrdy*(mkrdyp*(field(i,k,j+1)-field(i,k,j )) & -mkrdym*(field(i,k,j )-field(i,k,j-1)))) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE horizontal_diffusion !-----------------------------------------------------------------------------------------
SUBROUTINE horizontal_diffusion_3dmp ( name, field, tendency, mu, & 2 config_flags, base_3d, & msfu, msfv, msft, khdif, xkmhd, rdx, rdy, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte CHARACTER(LEN=1) , INTENT(IN ) :: name REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: field, & xkmhd, & base_3d REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv, & msft REAL , INTENT(IN ) :: rdx, & rdy, & khdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end REAL :: mrdx, mkrdxm, mkrdxp, & mrdy, mkrdym, mkrdyp, & rcoup REAL :: pr = 3. LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-2,ite) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-2,jte) DO j = j_start, j_end DO k=kts,ktf DO i = i_start, i_end mkrdxm=msfu(i,j)*0.5*(xkmhd(i,k,j)+xkmhd(i-1,k,j))*rdx*pr mkrdxp=msfu(i+1,j)*0.5*(xkmhd(i+1,k,j)+xkmhd(i,k,j))*rdx*pr mrdx=msft(i,j)*rdx mkrdym=msfv(i,j)*0.5*(xkmhd(i,k,j)+xkmhd(i,k,j-1))*rdy*pr mkrdyp=msfv(i,j+1)*0.5*(xkmhd(i,k,j+1)+xkmhd(i,k,j))*rdy*pr mrdy=msft(i,j)*rdy rcoup=mu(i,j) tendency(i,k,j)=tendency(i,k,j)+rcoup*( & mrdx*( mkrdxp*( field(i+1,k,j) -field(i ,k,j) & -base_3d(i+1,k,j)+base_3d(i ,k,j) ) & -mkrdxm*( field(i ,k,j) -field(i-1,k,j) & -base_3d(i ,k,j)+base_3d(i-1,k,j) ) ) & +mrdy*( mkrdyp*( field(i,k,j+1) -field(i,k,j ) & -base_3d(i,k,j+1)+base_3d(i,k,j ) ) & -mkrdym*( field(i,k,j ) -field(i,k,j-1) & -base_3d(i,k,j )+base_3d(i,k,j-1) ) ) & ) ENDDO ENDDO ENDDO END SUBROUTINE horizontal_diffusion_3dmp !-----------------------------------------------------------------------------------------
SUBROUTINE vertical_diffusion ( name, field, tendency, & 4 config_flags, & alt, mut, rdn, rdnw, kvdif, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte CHARACTER(LEN=1) , INTENT(IN ) :: name REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: field, & alt REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw REAL , INTENT(IN ) :: kvdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION(its:ite, jts:jte) :: vfluxm, vfluxp, zz REAL , DIMENSION(its:ite, 0:kte+1) :: vflux REAL :: rdz, rcoup LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) IF (name .EQ. 'w')THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) j_loop_w : DO j = j_start, j_end DO k=kts,ktf-1 DO i = i_start, i_end vflux(i,k)= (kvdif/alt(i,k,j))*rdnw(k)*(field(i,k+1,j)-field(i,k,j)) ENDDO ENDDO DO i = i_start, i_end vflux(i,ktf)=0. ENDDO DO k=kts+1,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j) & +rdn(k)*g*g/mut(i,j)/(0.5*(alt(i,k,j)+alt(i,k-1,j))) & *(vflux(i,k)-vflux(i,k-1)) ENDDO ENDDO ENDDO j_loop_w ELSE IF(name .EQ. 'm')THEN i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) j_loop_s : DO j = j_start, j_end DO k=kts,ktf-1 DO i = i_start, i_end vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) & *(field(i,k+1,j)-field(i,k,j)) ENDDO ENDDO DO i = i_start, i_end vflux(i,0)=vflux(i,1) ENDDO DO i = i_start, i_end vflux(i,ktf)=0. ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) & *rdnw(k)*(vflux(i,k)-vflux(i,k-1)) ENDDO ENDDO ENDDO j_loop_s ENDIF END SUBROUTINE vertical_diffusion !-------------------------------------------------------------------------------
SUBROUTINE vertical_diffusion_mp ( field, tendency, config_flags, & 2 base, & alt, mut, rdn, rdnw, kvdif, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: field, & alt REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, & rdnw, & base REAL , INTENT(IN ) :: kvdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION(its:ite, 0:kte+1) :: vflux REAL :: rdz, rcoup LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) j_loop_s : DO j = j_start, j_end DO k=kts,ktf-1 DO i = i_start, i_end vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) & *(field(i,k+1,j)-field(i,k,j)-base(k+1)+base(k)) ENDDO ENDDO DO i = i_start, i_end vflux(i,0)=vflux(i,1) ENDDO DO i = i_start, i_end vflux(i,ktf)=0. ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) & *rdnw(k)*(vflux(i,k)-vflux(i,k-1)) ENDDO ENDDO ENDDO j_loop_s END SUBROUTINE vertical_diffusion_mp !-------------------------------------------------------------------------------
SUBROUTINE vertical_diffusion_3dmp ( field, tendency, config_flags, & 2 base_3d, & alt, mut, rdn, rdnw, kvdif, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: field, & alt, & base_3d REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mut REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, & rdnw REAL , INTENT(IN ) :: kvdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION(its:ite, 0:kte+1) :: vflux REAL :: rdz, rcoup LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) j_loop_s : DO j = j_start, j_end DO k=kts,ktf-1 DO i = i_start, i_end vflux(i,k)=kvdif*rdn(k+1)/(0.5*(alt(i,k,j)+alt(i,k+1,j))) & *( field(i,k+1,j) -field(i,k,j) & -base_3d(i,k+1,j)+base_3d(i,k,j) ) ENDDO ENDDO DO i = i_start, i_end vflux(i,0)=vflux(i,1) ENDDO DO i = i_start, i_end vflux(i,ktf)=0. ENDDO DO k=kts,ktf DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+g*g/mut(i,j)/alt(i,k,j) & *rdnw(k)*(vflux(i,k)-vflux(i,k-1)) ENDDO ENDDO ENDDO j_loop_s END SUBROUTINE vertical_diffusion_3dmp !-------------------------------------------------------------------------------
SUBROUTINE vertical_diffusion_u ( field, tendency, & 2 config_flags, u_base, & alt, muu, rdn, rdnw, kvdif, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: field, & alt REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muu REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, u_base REAL , INTENT(IN ) :: kvdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION(its:ite, 0:kte+1) :: vflux REAL :: rdz, rcoup, zz LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) i_start = its i_end = ite j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_xs .or. specified ) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified ) i_end = MIN(ide-1,ite) j_loop_u : DO j = j_start, j_end DO k=kts,ktf-1 DO i = i_start, i_end vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i ,k ,j) & +alt(i-1,k ,j) & +alt(i ,k+1,j) & +alt(i-1,k+1,j) ) ) & *(field(i,k+1,j)-field(i,k,j) & -u_base(k+1) +u_base(k) ) ENDDO ENDDO DO i = i_start, i_end vflux(i,0)=vflux(i,1) ENDDO DO i = i_start, i_end vflux(i,ktf)=0. ENDDO DO k=kts,ktf-1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+ & g*g*rdnw(k)/muu(i,j)/(0.5*(alt(i-1,k,j)+alt(i,k,j)))* & (vflux(i,k)-vflux(i,k-1)) ENDDO ENDDO ENDDO j_loop_u END SUBROUTINE vertical_diffusion_u !-------------------------------------------------------------------------------
SUBROUTINE vertical_diffusion_v ( field, tendency, & 2 config_flags, v_base, & alt, muv, rdn, rdnw, kvdif, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: field, & alt REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: rdn, rdnw, v_base REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: muv REAL , INTENT(IN ) :: kvdif ! Local data INTEGER :: i, j, k, itf, jtf, ktf, jm1 INTEGER :: i_start, i_end, j_start, j_end REAL , DIMENSION(its:ite, 0:kte+1) :: vflux REAL :: rdz, rcoup, zz LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) i_start = its i_end = MIN(ite,ide-1) j_start = jts j_end = MIN(jte,jde-1) IF ( config_flags%open_ys .or. specified ) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified ) j_end = MIN(jde-1,jte) j_loop_v : DO j = j_start, j_end ! jm1 = max(j-1,1) jm1 = j-1 DO k=kts,ktf-1 DO i = i_start, i_end vflux(i,k)=kvdif*rdn(k+1)/(0.25*( alt(i,k ,j ) & +alt(i,k ,jm1) & +alt(i,k+1,j ) & +alt(i,k+1,jm1) ) ) & *(field(i,k+1,j)-field(i,k,j) & -v_base(k+1) +v_base(k) ) ENDDO ENDDO DO i = i_start, i_end vflux(i,0)=vflux(i,1) ENDDO DO i = i_start, i_end vflux(i,ktf)=0. ENDDO DO k=kts,ktf-1 DO i = i_start, i_end tendency(i,k,j)=tendency(i,k,j)+ & g*g*rdnw(k)/muv(i,j)/(0.5*(alt(i,k,jm1)+alt(i,k,j)))* & (vflux(i,k)-vflux(i,k-1)) ENDDO ENDDO ENDDO j_loop_v END SUBROUTINE vertical_diffusion_v !*************** end new mass coordinate routines !-------------------------------------------------------------------------------
SUBROUTINE calculate_full ( rfield, rfieldb, rfieldp, & 1 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfieldb, & rfieldp REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: rfield ! Local indices. INTEGER :: i, j, k, itf, jtf, ktf itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf rfield(i,k,j)=rfieldb(i,k,j)+rfieldp(i,k,j) ENDDO ENDDO ENDDO END SUBROUTINE calculate_full !------------------------------------------------------------------------------
SUBROUTINE coriolis ( ru, rv, rw, ru_tend, rv_tend, rw_tend, & 1 config_flags, & f, e, sina, cosa, fzm, fzp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: ru_tend, & rv_tend, & rw_tend REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: ru, & rv, & rw REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: f, & e, & sina, & cosa REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp ! Local indices. INTEGER :: i, j , k, ktf INTEGER :: i_start, i_end, j_start, j_end LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. ktf=MIN(kte,kde-1) ! coriolis for u-momentum equation i_start = its i_end = ite IF ( config_flags%open_xs .or. specified .or. & config_flags%nested) i_start = MAX(ids+1,its) IF ( config_flags%open_xe .or. specified .or. & config_flags%nested) i_end = MIN(ide-1,ite) DO j = jts, MIN(jte,jde-1) DO k=kts,ktf DO i = i_start, i_end ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(f(i,j)+f(i-1,j)) & *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) & - 0.5*(e(i,j)+e(i-1,j))*0.5*(cosa(i,j)+cosa(i-1,j)) & *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) ENDDO ENDDO IF ( (config_flags%open_xs) .and. (its == ids) ) THEN DO k=kts,ktf ru_tend(its,k,j)=ru_tend(its,k,j) + 0.5*(f(its,j)+f(its,j)) & *0.25*(rv(its,k,j+1)+rv(its,k,j+1)+rv(its,k,j)+rv(its,k,j)) & - 0.5*(e(its,j)+e(its,j))*0.5*(cosa(its,j)+cosa(its,j)) & *0.25*(rw(its,k+1,j)+rw(its,k,j)+rw(its,k+1,j)+rw(its,k,j)) ENDDO ENDIF IF ( (config_flags%open_xe) .and. (ite == ide) ) THEN DO k=kts,ktf ru_tend(ite,k,j)=ru_tend(ite,k,j) + 0.5*(f(ite-1,j)+f(ite-1,j)) & *0.25*(rv(ite-1,k,j+1)+rv(ite-1,k,j+1)+rv(ite-1,k,j)+rv(ite-1,k,j)) & - 0.5*(e(ite-1,j)+e(ite-1,j))*0.5*(cosa(ite-1,j)+cosa(ite-1,j)) & *0.25*(rw(ite-1,k+1,j)+rw(ite-1,k,j)+rw(ite-1,k+1,j)+rw(ite-1,k,j)) ENDDO ENDIF ENDDO ! coriolis term for v-momentum equation j_start = jts j_end = jte IF ( config_flags%open_ys .or. specified .or. & config_flags%nested) j_start = MAX(jds+1,jts) IF ( config_flags%open_ye .or. specified .or. & config_flags%nested) j_end = MIN(jde-1,jte) IF ( (config_flags%open_ys) .and. (jts == jds) ) THEN DO k=kts,ktf DO i=its,MIN(ide-1,ite) rv_tend(i,k,jts)=rv_tend(i,k,jts) - 0.5*(f(i,jts)+f(i,jts)) & *0.25*(ru(i,k,jts)+ru(i+1,k,jts)+ru(i,k,jts)+ru(i+1,k,jts)) & + 0.5*(e(i,jts)+e(i,jts))*0.5*(sina(i,jts)+sina(i,jts)) & *0.25*(rw(i,k+1,jts)+rw(i,k,jts)+rw(i,k+1,jts)+rw(i,k,jts)) ENDDO ENDDO ENDIF DO j=j_start, j_end DO k=kts,ktf DO i=its,MIN(ide-1,ite) rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(f(i,j)+f(i,j-1)) & *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) & + 0.5*(e(i,j)+e(i,j-1))*0.5*(sina(i,j)+sina(i,j-1)) & *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j)) ENDDO ENDDO ENDDO IF ( (config_flags%open_ye) .and. (jte == jde) ) THEN DO k=kts,ktf DO i=its,MIN(ide-1,ite) rv_tend(i,k,jte)=rv_tend(i,k,jte) - 0.5*(f(i,jte-1)+f(i,jte-1)) & *0.25*(ru(i,k,jte-1)+ru(i+1,k,jte-1)+ru(i,k,jte-1)+ru(i+1,k,jte-1)) & + 0.5*(e(i,jte-1)+e(i,jte-1))*0.5*(sina(i,jte-1)+sina(i,jte-1)) & *0.25*(rw(i,k+1,jte-1)+rw(i,k,jte-1)+rw(i,k+1,jte-1)+rw(i,k,jte-1)) ENDDO ENDDO ENDIF ! coriolis term for w-mometum DO j=jts,MIN(jte, jde-1) DO k=kts+1,ktf DO i=its,MIN(ite, ide-1) rw_tend(i,k,j)=rw_tend(i,k,j) + e(i,j)* & (cosa(i,j)*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))) & -sina(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)))) ENDDO ENDDO ENDDO END SUBROUTINE coriolis !------------------------------------------------------------------------------
SUBROUTINE curvature ( ru, rv, rw, u, v, w, ru_tend, rv_tend, rw_tend, & 1 config_flags, & msfu, msfv, fzm, fzp, rdx, rdy, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(INOUT) :: ru_tend, & rv_tend, & rw_tend REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: ru, & rv, & rw, & u, & v, & w REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, & msfv REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp REAL , INTENT(IN ) :: rdx, & rdy ! Local data ! INTEGER :: i, j, k, itf, jtf, ktf, kp1, im, ip, jm, jp INTEGER :: i, j, k, itf, jtf, ktf INTEGER :: i_start, i_end, j_start, j_end ! INTEGER :: irmin, irmax, jrmin, jrmax REAL , DIMENSION( its-1:ite , kts:kte, jts-1:jte ) :: vxgm LOGICAL :: specified specified = .false. if(config_flags%specified .or. config_flags%nested) specified = .true. itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ktf=MIN(kte,kde-1) ! irmin = ims ! irmax = ime ! jrmin = jms ! jrmax = jme ! IF ( config_flags%open_xs ) irmin = ids ! IF ( config_flags%open_xe ) irmax = ide-1 ! IF ( config_flags%open_ys ) jrmin = jds ! IF ( config_flags%open_ye ) jrmax = jde-1 ! Define v cross grad m at scalar points - vxgm(i,j) i_start = its-1 i_end = ite j_start = jts-1 j_end = jte IF ( ( config_flags%open_xs .or. specified .or. & config_flags%nested) .and. (its == ids) ) i_start = its IF ( ( config_flags%open_xe .or. specified .or. & config_flags%nested) .and. (ite == ide) ) i_end = ite-1 IF ( ( config_flags%open_ys .or. specified .or. & config_flags%nested) .and. (jts == jds) ) j_start = jts IF ( ( config_flags%open_ye .or. specified .or. & config_flags%nested) .and. (jte == jde) ) j_end = jte-1 DO j=j_start, j_end DO k=kts,ktf DO i=i_start, i_end vxgm(i,k,j)=0.5*(u(i,k,j)+u(i+1,k,j))*(msfv(i,j+1)-msfv(i,j))*rdy - & 0.5*(v(i,k,j)+v(i,k,j+1))*(msfu(i+1,j)-msfu(i,j))*rdx ENDDO ENDDO ENDDO ! Pick up the boundary rows for open (radiation) lateral b.c. ! Rather crude at present, we are assuming there is no ! variation in this term at the boundary. IF ( ( config_flags%open_xs .or. specified .or. & config_flags%nested) .and. (its == ids) ) THEN DO j = jts-1, jte DO k = kts, ktf vxgm(its-1,k,j) = vxgm(its,k,j) ENDDO ENDDO ENDIF IF ( ( config_flags%open_xe .or. specified .or. & config_flags%nested) .and. (ite == ide) ) THEN DO j = jts-1, jte DO k = kts, ktf vxgm(ite,k,j) = vxgm(ite-1,k,j) ENDDO ENDDO ENDIF IF ( ( config_flags%open_ys .or. specified .or. & config_flags%nested) .and. (jts == jds) ) THEN DO k = kts, ktf DO i = its-1, ite vxgm(i,k,jts-1) = vxgm(i,k,jts) ENDDO ENDDO ENDIF IF ( ( config_flags%open_ye .or. specified .or. & config_flags%nested) .and. (jte == jde) ) THEN DO k = kts, ktf DO i = its-1, ite vxgm(i,k,jte) = vxgm(i,k,jte-1) ENDDO ENDDO ENDIF ! curvature term for u momentum eqn. i_start = its IF ( config_flags%open_xs .or. specified .or. & config_flags%nested) i_start = MAX ( ids+1 , its ) IF ( config_flags%open_xe .or. specified .or. & config_flags%nested) i_end = MIN ( ide-1 , ite ) DO j=jts,MIN(jde-1,jte) DO k=kts,ktf DO i=i_start,i_end ru_tend(i,k,j)=ru_tend(i,k,j) + 0.5*(vxgm(i,k,j)+vxgm(i-1,k,j)) & *0.25*(rv(i-1,k,j+1)+rv(i,k,j+1)+rv(i-1,k,j)+rv(i,k,j)) & - u(i,k,j)*reradius & *0.25*(rw(i-1,k+1,j)+rw(i-1,k,j)+rw(i,k+1,j)+rw(i,k,j)) ENDDO ENDDO ENDDO ! curvature term for v momentum eqn. j_start = jts IF ( config_flags%open_ys .or. specified .or. & config_flags%nested) j_start = MAX ( jds+1 , jts ) IF ( config_flags%open_ye .or. specified .or. & config_flags%nested) j_end = MIN ( jde-1 , jte ) DO j=j_start,j_end DO k=kts,ktf DO i=its,MIN(ite,ide-1) rv_tend(i,k,j)=rv_tend(i,k,j) - 0.5*(vxgm(i,k,j)+vxgm(i,k,j-1)) & *0.25*(ru(i,k,j)+ru(i+1,k,j)+ru(i,k,j-1)+ru(i+1,k,j-1)) & + v(i,k,j)*reradius & *0.25*(rw(i,k+1,j-1)+rw(i,k,j-1)+rw(i,k+1,j)+rw(i,k,j)) ENDDO ENDDO ENDDO ! curvature term for vertical momentum eqn. DO j=jts,MIN(jte,jde-1) DO k=MAX(2,kts),ktf DO i=its,MIN(ite,ide-1) rw_tend(i,k,j)=rw_tend(i,k,j) + reradius* & (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))) & *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))) & +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))) & *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)))) ENDDO ENDDO ENDDO END SUBROUTINE curvature !------------------------------------------------------------------------------
SUBROUTINE decouple ( rr, rfield, field, name, config_flags, & fzm, fzp, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte CHARACTER(LEN=1) , INTENT(IN ) :: name REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rfield REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rr REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(OUT ) :: field REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, fzp ! Local data INTEGER :: i, j, k, itf, jtf, ktf ktf=MIN(kte,kde-1) IF (name .EQ. 'u')THEN itf=ite jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i-1,k,j))) ENDDO ENDDO ENDDO ELSE IF (name .EQ. 'v')THEN itf=MIN(ite,ide-1) jtf=jte DO j=jts,jtf DO k=kts,ktf DO i=its,itf field(i,k,j)=rfield(i,k,j)/(0.5*(rr(i,k,j)+rr(i,k,j-1))) ENDDO ENDDO ENDDO ELSE IF (name .EQ. 'w')THEN itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts+1,ktf DO i=its,itf field(i,k,j)=rfield(i,k,j)/(fzm(k)*rr(i,k,j)+fzp(k)*rr(i,k-1,j)) ENDDO ENDDO ENDDO DO j=jts,jtf DO i=its,itf field(i,kte,j) = 0. ENDDO ENDDO ELSE itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) ! For theta we will decouple tb and tp and add them to give t afterwards DO j=jts,jtf DO k=kts,ktf DO i=its,itf field(i,k,j)=rfield(i,k,j)/rr(i,k,j) ENDDO ENDDO ENDDO ENDIF END SUBROUTINE decouple !-------------------------------------------------------------------------------
SUBROUTINE mix_qv ( qv_tend, qv_mix, rho, & ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime, kms:kme , jms:jme ) , INTENT(IN ) :: qv_mix, & rho REAL , DIMENSION( ims:ime, kms:kme , jms:jme ) , INTENT(INOUT) :: qv_tend ! Local data INTEGER :: i, j, k, itf, jtf, ktf ktf=MIN(kte,kde-1) itf=MIN(ite,ide-1) jtf=MIN(jte,jde-1) DO j=jts,jtf DO k=kts,ktf DO i=its,itf qv_tend(i,k,j) = qv_tend(i,k,j) + qv_mix(i,k,j) ENDDO ENDDO ENDDO END SUBROUTINE mix_qv !-------------------------------------------------------------------------------
SUBROUTINE zero_tend ( tendency, & 17 ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) IMPLICIT NONE ! Input data INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tendency ! Local data INTEGER :: i, j, k, itf, jtf, ktf DO j = jts, jte DO k = kts, kte DO i = its, ite tendency(i,k,j) = 0. ENDDO ENDDO ENDDO END SUBROUTINE zero_tend !====================================================================== ! physics prep routines !======================================================================
SUBROUTINE phy_prep ( config_flags, & ! input 1 mu, u, v, p, pb, alt, ph, & ! input phb, t, tsk, moist, n_moist, & ! input mu_3d, rho, th_phy, p_phy , pi_phy , & ! output u_phy, v_phy, p8w, t_phy, t8w, & ! output z, z_at_w, dz8w, & ! output fzm, fzp, & ! params ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte ) !---------------------------------------------------------------------- IMPLICIT NONE !---------------------------------------------------------------------- TYPE(grid_config_rec_type) , INTENT(IN ) :: config_flags INTEGER , INTENT(IN ) :: ids, ide, jds, jde, kds, kde, & ims, ime, jms, jme, kms, kme, & its, ite, jts, jte, kts, kte INTEGER , INTENT(IN ) :: n_moist REAL, DIMENSION( ims:ime, kms:kme , jms:jme , n_moist ), INTENT(IN) :: moist REAL , DIMENSION( ims:ime, jms:jme ), INTENT(IN ) :: TSK, mu REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT( OUT) :: u_phy, & v_phy, & pi_phy, & p_phy, & p8w, & t_phy, & th_phy, & t8w, & mu_3d, & rho, & z, & dz8w, & z_at_w REAL , DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: pb, & p, & u, & v, & alt, & ph, & phb, & t REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end INTEGER :: i, j, k REAL :: w1, w2, z0, z1, z2 !----------------------------------------------------------------------- ! set up loop bounds for this grid's boundary conditions i_start = its i_end = min( ite,ide-1 ) j_start = jts j_end = min( jte,jde-1 ) k_start = kts k_end = min( kte, kde-1 ) ! compute thermodynamics and velocities at pressure points do j = j_start,j_end do k = k_start, k_end do i = i_start, i_end th_phy(i,k,j) = t(i,k,j) + t0 p_phy(i,k,j) = p(i,k,j) + pb(i,k,j) pi_phy(i,k,j) = (p_phy(i,k,j)/p1000mb)**rcp t_phy(i,k,j) = th_phy(i,k,j)*pi_phy(i,k,j) rho(i,k,j) = 1./alt(i,k,j)*(1.+moist(i,k,j,P_QV)) mu_3d(i,k,j) = mu(i,j) u_phy(i,k,j) = 0.5*(u(i,k,j)+u(i+1,k,j)) v_phy(i,k,j) = 0.5*(v(i,k,j)+v(i,k,j+1)) enddo enddo enddo ! compute z at w points do j = j_start,j_end do k = k_start, kte do i = i_start, i_end z_at_w(i,k,j) = (phb(i,k,j)+ph(i,k,j))/g enddo enddo enddo do j = j_start,j_end do k = k_start, kte-1 do i = i_start, i_end dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j) enddo enddo enddo do j = j_start,j_end do i = i_start, i_end dz8w(i,kte,j) = 0. enddo enddo ! compute z at p points (average of z at w points) do j = j_start,j_end do k = k_start, k_end do i = i_start, i_end z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) ) enddo enddo enddo ! interp t and p at w points do j = j_start,j_end do k = 2, k_end do i = i_start, i_end p8w(i,k,j) = fzm(k)*p_phy(i,k,j)+fzp(k)*p_phy(i,k-1,j) t8w(i,k,j) = fzm(k)*t_phy(i,k,j)+fzp(k)*t_phy(i,k-1,j) enddo enddo enddo ! extrapolate p and t to surface and top. ! we'll use an extrapolation in z for now do j = j_start,j_end do i = i_start, i_end ! bottom z0 = z_at_w(i,1,j) z1 = z(i,1,j) z2 = z(i,2,j) w1 = (z0 - z2)/(z1 - z2) w2 = 1. - w1 p8w(i,1,j) = w1*p_phy(i,1,j)+w2*p_phy(i,2,j) t8w(i,1,j) = w1*t_phy(i,1,j)+w2*t_phy(i,2,j) ! top z0 = z_at_w(i,kte,j) z1 = z(i,k_end,j) z2 = z(i,k_end-1,j) w1 = (z0 - z2)/(z1 - z2) w2 = 1. - w1 ! p8w(i,kde,j) = w1*p_phy(i,kde-1,j)+w2*p_phy(i,kde-2,j) !!! bug fix extrapolate ln(p) so p is positive definite p8w(i,kde,j) = exp(w1*log(p_phy(i,kde-1,j))+w2*log(p_phy(i,kde-2,j))) t8w(i,kde,j) = w1*t_phy(i,kde-1,j)+w2*t_phy(i,kde-2,j) enddo enddo END SUBROUTINE phy_prep !------------------------------------------------------------
SUBROUTINE moist_physics_prep_em( t_new, t_old, t0, rho, al, alb, & 1 p, p8w, p0, pb, ph, phb, pii, pf, & z, z_at_w, dz8w, & dt,h_diabatic, & config_flags,fzm, fzp, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) IMPLICIT NONE ! Here we construct full fields ! needed by the microphysics TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte REAL, INTENT(IN ) :: dt REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), & INTENT(IN ) :: al, & alb, & p, & pb, & ph, & phb REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fzm, & fzp REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), & INTENT( OUT) :: rho, & pii, & pf, & z, & z_at_w, & dz8w, & p8w, & h_diabatic REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), & INTENT(INOUT) :: t_new, & t_old REAL, INTENT(IN ) :: t0, p0 REAL :: z0,z1,z2,w1,w2 INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end INTEGER :: i, j, k !-------------------------------------------------------------------- ! set up loop bounds for this grid's boundary conditions i_start = its i_end = min( ite,ide-1 ) j_start = jts j_end = min( jte,jde-1 ) k_start = kts k_end = min( kte, kde-1 ) DO j = j_start, j_end DO k = k_start, kte DO i = i_start, i_end z_at_w(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g ENDDO ENDDO ENDDO do j = j_start,j_end do k = k_start, kte-1 do i = i_start, i_end dz8w(i,k,j) = z_at_w(i,k+1,j)-z_at_w(i,k,j) enddo enddo enddo do j = j_start,j_end do i = i_start, i_end dz8w(i,kte,j) = 0. enddo enddo ! compute full pii, rho, and z at the new time-level ! (needed for physics). ! convert perturbation theta to full theta DO j = j_start, j_end DO k = k_start, k_end DO i = i_start, i_end t_new(i,k,j) = t_new(i,k,j)-h_diabatic(i,k,j)*dt h_diabatic(i,k,j) = t_new(i,k,j) t_new(i,k,j) = t_new(i,k,j) + t0 t_old(i,k,j) = t_old(i,k,j) + t0 rho(i,k,j) = 1./(al(i,k,j)+alb(i,k,j)) pii(i,k,j) = ((p(i,k,j)+pb(i,k,j))/p0)**rcp z(i,k,j) = 0.5*(z_at_w(i,k,j) +z_at_w(i,k+1,j) ) pf(i,k,j) = p(i,k,j)+pb(i,k,j) ENDDO ENDDO ENDDO ! interp t and p at w points do j = j_start,j_end do k = 2, k_end do i = i_start, i_end p8w(i,k,j) = fzm(k)*pf(i,k,j)+fzp(k)*pf(i,k-1,j) enddo enddo enddo ! extrapolate p and t to surface and top. ! we'll use an extrapolation in z for now do j = j_start,j_end do i = i_start, i_end ! bottom z0 = z_at_w(i,1,j) z1 = z(i,1,j) z2 = z(i,2,j) w1 = (z0 - z2)/(z1 - z2) w2 = 1. - w1 p8w(i,1,j) = w1*pf(i,1,j)+w2*pf(i,2,j) ! top z0 = z_at_w(i,kte,j) z1 = z(i,k_end,j) z2 = z(i,k_end-1,j) w1 = (z0 - z2)/(z1 - z2) w2 = 1. - w1 ! p8w(i,kde,j) = w1*pf(i,kde-1,j)+w2*pf(i,kde-2,j) p8w(i,kde,j) = exp(w1*log(pf(i,kde-1,j))+w2*log(pf(i,kde-2,j))) enddo enddo END SUBROUTINE moist_physics_prep_em !------------------------------------------------------------------------------
SUBROUTINE moist_physics_finish_em( t_new, t_old, t0, mut, & 1 h_diabatic, dt, & config_flags, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) IMPLICIT NONE ! Here we construct full fields ! needed by the microphysics TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags INTEGER, INTENT(IN ) :: ids,ide, jds,jde, kds,kde INTEGER, INTENT(IN ) :: ims,ime, jms,jme, kms,kme INTEGER, INTENT(IN ) :: its,ite, jts,jte, kts,kte REAL, DIMENSION( ims:ime , kms:kme, jms:jme ), & INTENT(INOUT) :: t_new, & t_old, & h_diabatic REAL, DIMENSION( ims:ime , jms:jme ), INTENT(INOUT) :: mut REAL, INTENT(IN ) :: t0, dt INTEGER :: i_start, i_end, j_start, j_end, k_start, k_end INTEGER :: i, j, k !-------------------------------------------------------------------- ! set up loop bounds for this grid's boundary conditions i_start = its i_end = min( ite,ide-1 ) j_start = jts j_end = min( jte,jde-1 ) k_start = kts k_end = min( kte, kde-1 ) ! return theta to perturbation theta, calc new p, rho DO j = j_start, j_end DO k = k_start, k_end DO i = i_start, i_end t_new(i,k,j) = t_new(i,k,j) - t0 t_old(i,k,j) = t_old(i,k,j) - t0 h_diabatic(i,k,j) = (t_new(i,k,j)-h_diabatic(i,k,j))/dt ! h_diabatic(i,k,j) = 0. ENDDO ENDDO ENDDO END SUBROUTINE moist_physics_finish_em !----------------------------------------------------------------
SUBROUTINE init_module_big_step 1 END SUBROUTINE init_module_big_step END MODULE module_big_step_utilities_em