!WRF:MODEL_LAYER:PHYSICS
!
MODULE module_diffusion_em 1
USE module_configure
USE module_bc
USE module_state_description
USE module_big_step_utilities_em
USE module_model_constants
USE module_wrf_error
CONTAINS
!==========================================================================
SUBROUTINE cal_deform_and_div ( config_flags,u,v,w,div,defor11,defor22, & 1
defor33,defor12,defor13,defor23, &
u_base, v_base,msfu, msfv, msft, &
rdx, rdy, dn, dnw, rdz, rdzw, &
fnm,fnp,cf1,cf2,cf3,zx,zy, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! Mass coordinate model tke code, ported from height October 2001, WCS
!--------------------------------------------------------------------------
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 , INTENT(IN ) :: rdx, rdy
REAL , INTENT(IN ) :: cf1, cf2, cf3
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: msfu, &
msfv, &
msft
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
INTENT(IN ) :: u, &
v, &
w, &
zx, &
zy, &
rdz, &
rdzw
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), &
INTENT(INOUT) ::defor11, &
defor22, &
defor33, &
defor12, &
defor13, &
defor23, &
div
! Local vars
INTEGER :: i, j, k, ktf, ktes1, ktes2
INTEGER :: i_start, i_end, j_start, j_end
REAL :: tmp
REAL , DIMENSION( its:ite, jts:jte) :: mm
REAL , DIMENSION( its:ite, jts:jte) :: zzavg
! new
! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: zxavg
REAL , DIMENSION( its:ite, jts:jte) :: zeta_zd12
REAL , DIMENSION( its-2:ite+2, kts:kte, jts-2:jte+2) :: tmp1, &
hat, &
hatavg
! new
REAL :: tmpzx, tmpzy, tmpzeta_z, cft1, cft2
! zeta_dir_avg needs memory size
!-------------------------------------------------------------------------
! deforxx = deformation
ktes1=kte-1
ktes2=kte-2
cft2 = -0.5*dnw(ktes1)/dn(ktes1)
cft1 = 1.-cft2
ktf=MIN(kte,kde-1)
!*********************************************************
!* calculate defor at p points with larger domain *
!* defor11 *
!* defor22 *
!* defor33 *
!*********************************************************
i_start = its
i_end = min(ite,ide-1)
j_start = jts
j_end = min(jte,jde-1)
do j=j_start,j_end
do i=i_start,i_end
mm(i,j)=msft(i,j)*msft(i,j)
enddo
enddo
!---------------------------
! partial u / partial x
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end+1
hat(i,k,j)=u(i,k,j)/msfu(i,j)
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
hatavg(i,k,j)=0.5*(fnm(k)*(hat(i,k ,j)+hat(i+1, k,j))+ &
fnp(k)*(hat(i,k-1,j)+hat(i+1,k-1,j)))
enddo
enddo
enddo
! extrapolations to top and bottom of domain (to w levels)
do j=j_start,j_end
do i=i_start,i_end
hatavg(i,1,j) =0.5*(cf1*hat(i ,1,j)+cf2*hat(i ,2,j)+cf3*hat(i ,3,j)+ &
cf1*hat(i+1,1,j)+cf2*hat(i+1,2,j)+cf3*hat(i+1,3,j))
hatavg(i,kte,j)=0.5*( cft1*(hat(i,ktes1,j)+hat(i+1,ktes1,j)) &
+cft2*(hat(i,ktes2,j)+hat(i+1,ktes2,j)) )
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
! new
tmpzx =0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ &
zx(i,k+1,j)+zx(i+1,k+1,j) )
! tmp1(i,k,j)= (hatavg(i,k+1,j)-hatavg(i,k,j))*tmpzx*zeta_z(i,j)/dnw(k)
tmp1(i,k,j)= (hatavg(i,k+1,j)-hatavg(i,k,j))*tmpzx*rdzw(i,k,j)
enddo
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
tmp1(i,k,j)= mm(i,j)*(rdx*(hat(i+1,k,j)-hat(i,k,j))-tmp1(i,k,j))
enddo
enddo
enddo
!-----------------------------------
! defor11 : 2 partial u / partial x
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
defor11(i,k,j)= 2.*tmp1(i,k,j)
enddo
enddo
enddo
! defor11 done
!-----------------------------------
! Div : add partial u / partial x
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
div(i,k,j)=tmp1(i,k,j)
enddo
enddo
enddo
!---------------------------
! partial v / partial y
do j=j_start,j_end+1
do k=kts,ktf
do i=i_start,i_end
hat(i,k,j)=v(i,k,j)/msfv(i,j)
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
hatavg(i,k,j)=0.5*(fnm(k)*(hat(i,k ,j)+hat(i,k ,j+1))+ &
fnp(k)*(hat(i,k-1,j)+hat(i,k-1,j+1)))
enddo
enddo
enddo
! extrapolations to top and bottom of domain (to w levels)
do j=j_start,j_end
do i=i_start,i_end
hatavg(i,1,j) =0.5*(cf1*hat(i,1,j )+cf2*hat(i,2,j )+cf3*hat(i,3,j )+ &
cf1*hat(i,1,j+1)+cf2*hat(i,2,j+1)+cf3*hat(i,3,j+1))
hatavg(i,kte,j)=0.5*( cft1*(hat(i,ktes1,j)+hat(i,ktes1,j+1)) &
+cft2*(hat(i,ktes2,j)+hat(i,ktes2,j+1)) )
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
! new
tmpzy =0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ &
zy(i,k+1,j)+zy(i,k+1,j+1) )
! tmp1(i,k,j)= (hatavg(i,k+1,j)-hatavg(i,k,j))*tmpzy*zeta_z(i,j)/dnw(k)
tmp1(i,k,j)= (hatavg(i,k+1,j)-hatavg(i,k,j))*tmpzy*rdzw(i,k,j)
enddo
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
tmp1(i,k,j)= mm(i,j)*(rdy*(hat(i,k,j+1)-hat(i,k,j))-tmp1(i,k,j))
enddo
enddo
enddo
!-------------------------------
! defor22 : 2 partial v / partial y
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
defor22(i,k,j)= 2.*tmp1(i,k,j)
enddo
enddo
enddo
! defor22 done
!-------------------------------
! Div : add partial v / partial y
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
div(i,k,j)=div(i,k,j)+tmp1(i,k,j)
enddo
enddo
enddo
!-------------------------
! partial w / partial z
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
! tmp1(i,k,j)= (w(i,k+1,j)-w(i,k,j))*zeta_z(i,j)/dnw(k)
tmp1(i,k,j)= (w(i,k+1,j)-w(i,k,j))*rdzw(i,k,j)
enddo
enddo
enddo
!----------------------------------
! defor33 : 2 partial w / parital z
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
defor33(i,k,j)= 2.*tmp1(i,k,j)
enddo
enddo
enddo
! defor33 done
!----------------------------------
! Div : add partial w / parital z
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
div(i,k,j)=div(i,k,j)+tmp1(i,k,j)
enddo
enddo
enddo
! Div done
!*****************************************************************
!* calculate defor at vorticity points with larger domain *
!* defor12 *
!*****************************************************************
! calculate defor at vorticity points
i_start = its
i_end = ite
j_start = jts
j_end = jte
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-1,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-1,jte)
!------------------------
! partial u / partial y
do j=j_start,j_end
do i=i_start, i_end
mm(i,j)=0.25*(msfu(i,j-1)+msfu(i,j))*(msfv(i-1,j)+msfv(i,j))
enddo
enddo
do j=j_start-1,j_end
do k=kts,ktf
do i=i_start,i_end
hat(i,k,j)=u(i,k,j)/msfu(i,j)
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
hatavg(i,k,j)=0.5*(fnm(k)*(hat(i,k ,j-1)+hat(i,k ,j))+ &
fnp(k)*(hat(i,k-1,j-1)+hat(i,k-1,j)))
enddo
enddo
enddo
! extrapolations to top and bottom of domain (to w levels)
do j=j_start,j_end
do i=i_start,i_end
hatavg(i,1,j) =0.5*(cf1*hat(i,1,j-1)+cf2*hat(i,2,j-1)+cf3*hat(i,3,j-1)+ &
cf1*hat(i,1,j)+cf2*hat(i,2,j)+cf3*hat(i,3,j))
hatavg(i,kte,j)=0.5*( cft1*(hat(i,ktes1,j-1)+hat(i,ktes1,j)) &
+cft2*(hat(i,ktes2,j-1)+hat(i,ktes2,j)) )
enddo
enddo
! do j=j_start, j_end
! do k=kts,ktf
! do i=i_start, i_end
! new
! zeta_zd12(i,j)=0.25*( zeta_z(i-1,j )+zeta_z(i,j )+ &
! zeta_z(i-1,j-1)+zeta_z(i,j-1) )
! enddo
! enddo
! enddo
do j=j_start, j_end
do k=kts,ktf
do i=i_start, i_end
! new
tmpzy =0.25*( zy(i-1,k ,j)+zy(i,k ,j)+&
zy(i-1,k+1,j)+zy(i,k+1,j) )
tmp1(i,k,j)= (hatavg(i,k+1,j)-hatavg(i,k,j))* &
tmpzy*0.5*(rdzw(i,k,j)+rdzw(i-1,k,j))
enddo
enddo
enddo
!-------------------------------------
! deform12: Add partial u / partial y
do j=j_start, j_end
do k=kts,ktf
do i=i_start, i_end
defor12(i,k,j)= mm(i,j)*(rdy*(hat(i,k,j)-hat(i,k,j-1))-tmp1(i,k,j))
enddo
enddo
enddo
!-------------------------------------
! partial v over partial x
do j=j_start,j_end
do k=kts,ktf
do i=i_start-1,i_end
hat(i,k,j)=v(i,k,j)/msfv(i,j)
enddo
enddo
enddo
do j=j_start, j_end
do k=kts+1,ktf
do i=i_start, i_end
hatavg(i,k,j)=0.5*(fnm(k)*(hat(i-1,k ,j)+hat(i,k ,j))+ &
fnp(k)*(hat(i-1,k-1,j)+hat(i,k-1,j)))
enddo
enddo
enddo
! extrapolations to top and bottom of domain (to w levels)
do j=j_start,j_end
do i=i_start,i_end
hatavg(i,1,j) =0.5*(cf1*hat(i-1,1,j)+cf2*hat(i-1,2,j)+cf3*hat(i-1,3,j)+ &
cf1*hat(i,1,j)+cf2*hat(i,2,j)+cf3*hat(i,3,j))
hatavg(i,kte,j)=0.5*( cft1*(hat(i,ktes1,j)+hat(i-1,ktes1,j)) &
+cft2*(hat(i,ktes2,j)+hat(i-1,ktes2,j)) )
enddo
enddo
do j=j_start, j_end
do k=kts,ktf
do i=i_start, i_end
! new
tmpzx =0.25*( zx(i,k ,j-1)+zx(i,k ,j)+ &
zx(i,k+1,j-1)+zx(i,k+1,j) )
tmp1(i,k,j)= (hatavg(i,k+1,j)-hatavg(i,k,j))* &
tmpzx*0.5*(rdzw(i,k,j)+rdzw(i,k,j-1))
! tmpzx*zeta_zd12(i,j)/dnw(k)
enddo
enddo
enddo
!------------------------------------
! deform12: Add partial v / partial x
do j=j_start, j_end
do k=kts,ktf
do i=i_start, i_end
defor12(i,k,j)= defor12(i,k,j) + &
mm(i,j)*(rdx*(hat(i,k,j)-hat(i-1,k,j))-tmp1(i,k,j))
enddo
enddo
enddo
!------------------------------------------------
! update boundary !!! might need to change later
!
if (.not. config_flags%periodic_x .and. i_start .eq. ids+1) then
do j=jts,jte
do k=kts,kte
defor12(ids,k,j)= defor12(ids+1,k,j)
enddo
enddo
endif
!
if (.not. config_flags%periodic_y .and. j_start .eq. jds+1) then
do k=kts,kte
do i=its,ite
defor12(i,k,jds)= defor12(i,k,jds+1)
enddo
enddo
endif
!
if (.not. config_flags%periodic_x .and. i_end .eq. ide-1) then
do j=jts,jte
do k=kts,kte
defor12(ide,k,j)= defor12(ide-1,k,j)
enddo
enddo
endif
if (.not. config_flags%periodic_y .and. j_end .eq. jde-1) then
do k=kts,kte
do i=its,ite
defor12(i,k,jde)= defor12(i,k,jde-1)
enddo
enddo
endif
! defor12 done
!*******************************************
!* calculate defor at u points *
!* defor13 *
!*******************************************
i_start = its
i_end = min(ite,ide-1)
j_start = jts
j_end = min(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%periodic_x ) i_end = min(ite,ide)
IF ( config_flags%periodic_y ) j_end = min(jte,jde)
do j=jts,jte
do i=its,ite
mm(i,j)=msfu(i,j)*msfu(i,j)
enddo
enddo
! for both defor13 and defor23
! do j=j_start-1,min(jte,jde-1)
!! do j=j_start,jte
do j=j_start,j_end
do k=kts,kte
! do i=i_start,min(ite,ide-1)
!! do i=i_start,ite
do i=i_start,i_end
hat(i,k,j)=w(i,k,j)/msft(i,j)
enddo
enddo
enddo
i=i_start-1
do j=j_start,min(jte,jde-1)
do k=kts,kte
hat(i,k,j)=w(i,k,j)/msft(i,j)
enddo
enddo
j=j_start-1
do k=kts,kte
do i=i_start,min(ite,ide-1)
hat(i,k,j)=w(i,k,j)/msft(i,j)
enddo
enddo
!----------------------------------
! defor13
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
hatavg(i,k,j)=0.25*(hat(i,k,j)+hat(i,k+1,j)+hat(i-1,k,j)+hat(i-1,k+1,j))
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
tmp1(i,k,j)= (hatavg(i,k,j)-hatavg(i,k-1,j))*zx(i,k,j)* &
0.5*(rdz(i,k,j)+rdz(i-1,k,j))
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
defor13(i,k,j)= mm(i,j)*(rdx*(hat(i,k,j)-hat(i-1,k,j))-tmp1(i,k,j))
enddo
enddo
enddo
do j=j_start,j_end
do i=i_start,i_end
defor13(i,kts,j )= 0.
defor13(i,ktf+1,j)= 0.
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
tmp1(i,k,j)= (u(i,k,j)-u_base(k)-u(i,k-1,j)+u_base(k-1))* &
0.5*(rdz(i,k,j)+rdz(i-1,k,j))
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
defor13(i,k,j)=defor13(i,k,j)+tmp1(i,k,j)
enddo
enddo
enddo
! defor13 done
!-------------
! defor23
i_start = its
i_end = min(ite,ide-1)
j_start = jts
j_end = min(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%periodic_y ) j_end = min(jte,jde)
do j=jts,jte
do i=its,ite
mm(i,j)=msfv(i,j)*msfv(i,j)
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
hatavg(i,k,j)=0.25*(hat(i,k,j)+hat(i,k+1,j)+hat(i,k,j-1)+hat(i,k+1,j-1))
! new
! zxavg(i,k,j)=0.5*(zy(i,k,j)+zy(i,k,j-1))
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
! new
! tmpzeta_z = 0.5*(zeta_z(i,j-1)+zeta_z(i,j))
! tmp1(i,k,j)= (hatavg(i,k,j)-hatavg(i,k-1,j))*zy(i,k,j)*tmpzeta_z/dn(k)
tmp1(i,k,j)= (hatavg(i,k,j)-hatavg(i,k-1,j)) &
*zy(i,k,j)*0.5*(rdz(i,k,j)+rdz(i,k,j-1))
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
defor23(i,k,j)= mm(i,j)*(rdy*(hat(i,k,j)-hat(i,k,j-1))-tmp1(i,k,j))
enddo
enddo
enddo
do j=j_start,j_end
do i=i_start,i_end
defor23(i,kts,j )= 0.
defor23(i,ktf+1,j)= 0.
enddo
enddo
! do j=j_start,j_end
! do i=i_start,i_end
! zzavg(i,j)=0.5*(zeta_z(i,j)+zeta_z(i,j-1))
! enddo
! enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
tmp1(i,k,j)= (v(i,k,j)-v_base(k)-v(i,k-1,j)+v_base(k-1))* &
0.5*(rdz(i,k,j)+rdz(i,k,j-1))
! zzavg(i,j)/dn(k)
enddo
enddo
enddo
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
defor23(i,k,j)= defor23(i,k,j) + tmp1(i,k,j)
enddo
enddo
enddo
! defor23 done
!------------------------------------------------
! update boundary !!! might need to change later
!
if (.not. config_flags%periodic_x .and. i_start .eq. ids+1) then
do j=jts,jte
do k=kts,kte
defor13(ids,k,j)= defor13(ids+1,k,j)
defor23(ids,k,j)= defor23(ids+1,k,j)
enddo
enddo
endif
!
if (.not. config_flags%periodic_y .and. j_start .eq. jds+1) then
do k=kts,kte
do i=its,ite
defor13(i,k,jds)= defor13(i,k,jds+1)
defor23(i,k,jds)= defor23(i,k,jds+1)
enddo
enddo
endif
!
if (.not. config_flags%periodic_x .and. i_end .eq. ide-1) then
do j=jts,jte
do k=kts,kte
defor13(ide,k,j)= defor13(ide-1,k,j)
defor23(ide,k,j)= defor23(ide-1,k,j)
enddo
enddo
endif
if (.not. config_flags%periodic_y .and. j_end .eq. jde-1) then
do k=kts,kte
do i=its,ite
defor13(i,k,jde)= defor13(i,k,jde-1)
defor23(i,k,jde)= defor23(i,k,jde-1)
enddo
enddo
endif
END SUBROUTINE cal_deform_and_div
!===============================================================================
SUBROUTINE calculate_km_kh ( config_flags,dt,dampcoef,zdamp,damp_opt, & 1,7
xkmh,xkmhd,xkmv,xkhh,xkhv,BN2, &
khdif,kvdif,div,defor11,defor22,defor33,defor12, &
defor13,defor23,tke,p8w,t8w,theta,t,p,moist, &
dn,dnw,dx,dy,rdz,rdzw,cr_len,n_moist, &
cf1, cf2, cf3, warm_rain, &
kh_tke_upper_bound, kv_tke_upper_bound, &
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, damp_opt
LOGICAL , INTENT(IN ) :: warm_rain
REAL , INTENT(IN ) :: cr_len,dx,dy,zdamp,dt,dampcoef
REAL, INTENT(IN ) :: cf1, cf2, cf3
REAL, INTENT(IN ) :: khdif,kvdif
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT(INOUT) &
:: moist
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmv, &
xkmh, &
xkmhd, &
xkhv, &
xkhh, &
BN2
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
defor11, &
defor22, &
defor33, &
defor12, &
defor13, &
defor23, &
div, &
rdz, &
rdzw
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: p8w, &
t8w, &
theta, &
t, &
p
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: tke
REAL, INTENT( IN) :: kh_tke_upper_bound, kv_tke_upper_bound
! LOCAL VAR
INTEGER:: i_start,i_end,j_start,j_end,ktf,i,j,k
!-------------------------------------------------------------------------------
! calculate deformation and divergence
ktf=min(kte,kde-1)
i_start = its
i_end = min(ite,ide-1)
j_start = jts
j_end = min(jte,jde-1)
CALL calculate_N2( config_flags,BN2,moist,theta,t,p, &
p8w,t8w,dnw,dn,rdz,rdzw, &
n_moist,cf1,cf2,cf3, &
warm_rain, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
! choice scheme for calculating diffusion coefficients
km_coef: SELECT CASE(config_flags%km_opt)
CASE (1)
CALL isotropic_km( config_flags,xkmh,xkmhd,xkmv, &
xkhh,xkhv,khdif,kvdif, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CASE (2)
CALL tke_km( config_flags,xkmh,xkmhd,xkmv, &
xkhh,xkhv,tke,p8w,t8w,theta, &
rdz, rdzw, dx,dy,cr_len, &
kh_tke_upper_bound, kv_tke_upper_bound, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CASE (3)
CALL smag_km( config_flags,xkmh,xkmhd,xkmv,xkhh,xkhv,BN2, &
div,defor11,defor22,defor33,defor12, &
defor13,defor23, &
rdzw,dx,dy,cr_len, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CASE (4)
CALL smag2d_km( config_flags,xkmh,xkmhd,xkmv,xkhh,xkhv, &
defor11,defor22,defor12, &
rdzw,dx,dy, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
CASE DEFAULT
CALL wrf_error_fatal( 'Please choose a diffusion coefficient scheme' )
END SELECT km_coef
IF (damp_opt .eq. 1) THEN
CALL cal_dampkm( config_flags,xkmhd,xkhh,xkmv,xkhv, &
dx,dy,dt,dampcoef, &
rdz, rdzw, zdamp, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte )
ENDIF
END SUBROUTINE calculate_km_kh
!===============================================================================
SUBROUTINE cal_dampkm( config_flags,xkmhd,xkhh,xkmv,xkhv, & 1
dx,dy,dt,dampcoef, &
rdz, rdzw ,zdamp, &
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 , INTENT(IN ) :: zdamp,dx,dy,dt,dampcoef
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: xkmhd, &
xkhh , &
xkmv , &
xkhv
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdz, &
rdzw
! LOCAL VARS
INTEGER :: i_start, i_end, j_start, j_end, ktf, ktfm1, i, j, k
REAL :: kmmax,kmmvmax,degrad90,dz,tmp
REAL , DIMENSION( its:ite ) :: deltaz
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: dampk,dampkv
ktf = min(kte,kde-1)
ktfm1 = ktf-1
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
kmmax=dx*dx/dt
degrad90=DEGRAD*90.
DO j = j_start, j_end
k=ktf
DO i = i_start, i_end
! deltaz(i)=0.5*dnw(k)/zeta_z(i,j)
! dz=dnw(k)/zeta_z(i,j)
dz = 1./rdzw(i,k,j)
deltaz(i) = 0.5*dz
kmmvmax=dz*dz/dt
tmp=min(deltaz(i)/zdamp,1.)
dampk(i,k,j)=cos(degrad90*tmp)*kmmax*dampcoef
dampkv(i,k,j)=cos(degrad90*tmp)*kmmvmax*dampcoef
ENDDO
DO k = ktfm1,kts,-1
DO i = i_start, i_end
! deltaz(i)=deltaz(i)+dn(k)/zeta_z(i,j)
! dz=dnw(k)/zeta_z(i,j)
dz = 1./rdz(i,k,j)
deltaz(i) = deltaz(i) + dz
dz = 1./rdzw(i,k,j)
kmmvmax=dz*dz/dt
tmp=min(deltaz(i)/zdamp,1.)
dampk(i,k,j)=cos(degrad90*tmp)*kmmax*dampcoef
dampkv(i,k,j)=cos(degrad90*tmp)*kmmvmax*dampcoef
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
xkmhd(i,k,j)=max(xkmhd(i,k,j),dampk(i,k,j))
xkhh(i,k,j)=max(xkhh(i,k,j),dampk(i,k,j))
xkmv(i,k,j)=max(xkmv(i,k,j),dampkv(i,k,j))
xkhv(i,k,j)=max(xkhv(i,k,j),dampkv(i,k,j))
ENDDO
ENDDO
ENDDO
END SUBROUTINE cal_dampkm
!===============================================================================
SUBROUTINE calculate_N2( config_flags,BN2,moist,theta,t,p, & 1
p8w,t8w,dnw,dn,rdz,rdzw, &
n_moist,cf1,cf2,cf3, &
warm_rain, &
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
LOGICAL , INTENT(IN ) :: warm_rain
REAL, INTENT(IN ) :: cf1, cf2, cf3
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: BN2
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdz
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdzw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), INTENT(INOUT) &
:: moist
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: theta, &
t, &
p, &
p8w, &
t8w
! Local VAR
INTEGER :: i, j, k, ktf, ispe, ktes1, ktes2
INTEGER :: i_start, i_end, j_start, j_end
REAL :: coefa,thetaep1,thetaem1,qc_cr,es,tc,qlpqi,qsw,qsi, &
tmpdz,xlvqv,thetaesfc,thetasfc,qvtop,qvsfc,thetatop, &
thetaetop
REAL , DIMENSION( its:ite, jts:jte) :: tmp1sfc, &
tmp1top
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmp1, &
qvs, &
qctmp
!------------------------------------------------------------------------------
! qc_cr is in Kg/Kg
qc_cr = 0.00001
ktf=MIN(kte,kde-1)
ktes1=kte-1
ktes2=kte-2
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
!
IF (P_QC .ge. PARAM_FIRST_SCALAR) THEN
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
qctmp(i,k,j)=moist(i,k,j,P_QC)
ENDDO
ENDDO
ENDDO
ELSE
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
qctmp(i,k,j)=0.
ENDDO
ENDDO
ENDDO
ENDIF
!
DO j = jts,jte
DO k = kts,kte
DO i = its,ite
tmp1(i,k,j)=0.
ENDDO
ENDDO
ENDDO
!
DO j = jts,jte
DO i = its,ite
tmp1sfc(i,j)=0.
tmp1top(i,j)=0.
ENDDO
ENDDO
!
DO ispe=PARAM_FIRST_SCALAR,n_moist
IF (ispe .eq. P_QV .or. ispe .eq. P_QC .or. ispe .eq. P_QI) THEN
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tmp1(i,k,j)=tmp1(i,k,j)+moist(i,k,j,ispe)
ENDDO
ENDDO
ENDDO
!
DO j = j_start, j_end
DO i = i_start, i_end
tmp1sfc(i,j)=tmp1sfc(i,j)+cf1*moist(i,1,j,ispe)+cf2*moist(i,2,j,ispe)+&
cf3*moist(i,3,j,ispe)
tmp1top(i,j)=tmp1top(i,j)+moist(i,ktes1,j,ispe)+ &
(moist(i,ktes1,j,ispe)-moist(i,ktes2,j,ispe)) &
*0.5*dnw(ktes1)/dn(ktes1)
ENDDO
ENDDO
ENDIF
ENDDO
! calculate saturation mixing ratio
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tc=t(i,k,j)-SVPT0
es=1000.*SVP1*exp(SVP2*tc/(t(i,k,j)-SVP3) )
qvs(i,k,j)=EP_2*es/(p(i,k,j)-es)
ENDDO
ENDDO
ENDDO
!
DO j = j_start, j_end
DO k = kts+1, ktf-1
DO i = i_start, i_end
! tmpdz=dn(k+1)+dn(k)
tmpdz=1./rdz(i,k,j)+1./rdz(i,k+1,j)
IF (moist(i,k,j,P_QV) .ge. qvs(i,k,j) .or. qctmp(i,k,j) .ge. qc_cr) THEN
xlvqv=XLV*moist(i,k,j,P_QV)
coefa=(1.+xlvqv/R_d/t(i,k,j))/ &
(1.+XLV*xlvqv/Cp/R_v/t(i,k,j)/t(i,k,j))/theta(i,k,j)
thetaep1=theta(i,k+1,j)*exp(XLV*moist(i,k+1,j,P_QV)/Cp/t(i,k+1,j))
thetaem1=theta(i,k-1,j)*exp(XLV*moist(i,k-1,j,P_QV)/Cp/t(i,k-1,j))
! BN2(i,k,j)=g*zeta_z(i,j)*(coefa*(thetaep1-thetaem1)/tmpdz-&
! (tmp1(i,k+1,j)-tmp1(i,k-1,j))/tmpdz &
BN2(i,k,j)=g*(coefa*(thetaep1-thetaem1)/tmpdz-&
(tmp1(i,k+1,j)-tmp1(i,k-1,j))/tmpdz &
)
ELSE
BN2(i,k,j)=g*( (theta(i,k+1,j)-theta(i,k-1,j))/theta(i,k,j)/tmpdz + &
0.61*(moist(i,k+1,j,P_QV)-moist(i,k-1,j,P_QV))/tmpdz &
)
! BN2(i,k,j)=g*zeta_z(i,j)* &
! ( (theta(i,k+1,j)-theta(i,k-1,j))/theta(i,k,j)/tmpdz + &
! 0.61*(moist(i,k+1,j,P_QV)-moist(i,k-1,j,P_QV))/tmpdz &
! )
ENDIF
ENDDO
ENDDO
ENDDO
k=kts
DO j = j_start, j_end
DO i = i_start, i_end
! tmpdz=dn(k+1)+0.5*dnw(k)
tmpdz=1./rdz(i,k+1,j)+0.5/rdzw(i,k,j)
thetasfc=T8w(i,kts,j)/(p8w(i,k,j)/p1000mb)**(R_d/Cp)
qvsfc=cf1*moist(i,1,j,P_QV)+cf2*moist(i,2,j,P_QV)+cf3*moist(i,3,j,P_QV)
IF (moist(i,k,j,P_QV) .ge. qvs(i,k,j) .or. qctmp(i,k,j) .ge. qc_cr) THEN
xlvqv=XLV*moist(i,k,j,P_QV)
coefa=(1.+xlvqv/R_d/t(i,k,j))/ &
(1.+XLV*xlvqv/Cp/R_v/t(i,k,j)/t(i,k,j))/theta(i,k,j)
thetaep1=theta(i,k+1,j)*exp(XLV*moist(i,k+1,j,P_QV)/Cp/t(i,k+1,j))
thetaesfc=thetasfc*exp(XLV*qvsfc/Cp/t8w(i,kts,j))
BN2(i,k,j)=g*(coefa*(thetaep1-thetaesfc)/tmpdz- &
(tmp1(i,k+1,j)-tmp1sfc(i,j))/tmpdz &
)
! BN2(i,k,j)=g*zeta_z(i,j)*(coefa*(thetaep1-thetaesfc)/tmpdz- &
! (tmp1(i,k+1,j)-tmp1sfc(i,j))/tmpdz &
! )
ELSE
BN2(i,k,j)=g* ( (theta(i,k+1,j)-thetasfc)/theta(i,k,j)/tmpdz + &
0.61*(moist(i,k+1,j,P_QV)-qvsfc)/tmpdz &
)
! BN2(i,k,j)=g*zeta_z(i,j)* &
! ( (theta(i,k+1,j)-thetasfc)/theta(i,k,j)/tmpdz + &
! 0.61*(moist(i,k+1,j,P_QV)-qvsfc)/tmpdz &
! )
ENDIF
ENDDO
ENDDO
!
k=ktf
DO j = j_start, j_end
DO i = i_start, i_end
! tmpdz=dn(k)+0.5*dnw(k)
tmpdz=1./rdz(i,k,j)+0.5/rdzw(i,k,j)
thetatop=T8w(i,kde,j)/(p8w(i,kde,j)/p1000mb)**(R_d/Cp)
qvtop=moist(i,ktes1,j,P_QV)+(moist(i,ktes1,j,P_QV)-moist(i,ktes2,j,P_QV)) &
*0.5*dnw(ktes1)/dn(ktes1)
IF (moist(i,k,j,P_QV) .ge. qvs(i,k,j) .or. qctmp(i,k,j) .ge. qc_cr) THEN
xlvqv=XLV*moist(i,k,j,P_QV)
coefa=(1.+xlvqv/R_d/t(i,k,j))/ &
(1.+XLV*xlvqv/Cp/R_v/t(i,k,j)/t(i,k,j))/theta(i,k,j)
thetaetop=thetatop*exp(XLV*qvtop/Cp/t8w(i,kde,j))
thetaem1=theta(i,k-1,j)*exp(XLV*moist(i,k-1,j,P_QV)/Cp/t(i,k-1,j))
BN2(i,k,j)=g*(coefa*(thetaetop-thetaem1)/tmpdz-&
(tmp1top(i,j)-tmp1(i,k-1,j))/tmpdz &
)
! BN2(i,k,j)=g*zeta_z(i,j)*(coefa*(thetaetop-thetaem1)/tmpdz-&
! (tmp1top(i,j)-tmp1(i,k-1,j))/tmpdz &
! )
ELSE
BN2(i,k,j)=g*( (thetatop-theta(i,k-1,j))/theta(i,k,j)/tmpdz + &
0.61*(qvtop-moist(i,k-1,j,P_QV))/tmpdz &
)
! BN2(i,k,j)=g*zeta_z(i,j)* &
! ( (thetatop-theta(i,k-1,j))/theta(i,k,j)/tmpdz + &
! 0.61*(qvtop-moist(i,k-1,j,P_QV))/tmpdz &
! )
ENDIF
ENDDO
ENDDO
END SUBROUTINE calculate_N2
!===============================================================================
SUBROUTINE isotropic_km( config_flags, & 1
xkmh,xkmhd,xkmv,xkhh,xkhv,khdif,kvdif, &
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 , INTENT(IN ) :: khdif,kvdif
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
xkmhd, &
xkmv, &
xkhh, &
xkhv
! LOCAL VARS
INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
REAL :: khdif3,kvdif3
! ktf = min(kte,kde-1)
ktf = kte
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
! khdif3=khdif*3.
! kvdif3=kvdif*3.
khdif3=khdif*prandtl
kvdif3=kvdif*prandtl
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
xkmh(i,k,j)=khdif
xkmhd(i,k,j)=khdif
xkmv(i,k,j)=kvdif
xkhh(i,k,j)=khdif3
xkhv(i,k,j)=kvdif3
ENDDO
ENDDO
ENDDO
END SUBROUTINE isotropic_km
!===============================================================================
SUBROUTINE smag_km( config_flags,xkmh,xkmhd,xkmv,xkhh,xkhv,BN2, & 1
div,defor11,defor22,defor33,defor12, &
defor13,defor23, &
rdzw,dx,dy,cr_len_in, &
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 , INTENT(IN ) :: cr_len_in, dx, dy
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: BN2, &
rdzw
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
xkmhd, &
xkmv, &
xkhh, &
xkhv
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
defor11, &
defor22, &
defor33, &
defor12, &
defor13, &
defor23, &
div
! LOCAL VARS
INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
REAL :: deltas, tmp, pr, mlen_h, mlen_v, cr_len
REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
!-------------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
pr=1./3.
cr_len = cr_len_in
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
def2(i,k,j)=0.5*(defor11(i,k,j)*defor11(i,k,j) + &
defor22(i,k,j)*defor22(i,k,j) + &
defor33(i,k,j)*defor33(i,k,j))
enddo
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
tmp=0.25*(defor12(i ,k,j)+defor12(i ,k,j+1)+ &
defor12(i+1,k,j)+defor12(i+1,k,j+1))
def2(i,k,j)=def2(i,k,j)+0.5*tmp*tmp
enddo
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
tmp=0.25*(defor13(i ,k+1,j)+defor13(i ,k,j)+ &
defor13(i+1,k+1,j)+defor13(i+1,k,j))
def2(i,k,j)=def2(i,k,j)+0.5*tmp*tmp
enddo
enddo
enddo
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
tmp=0.25*(defor23(i,k+1,j )+defor23(i,k,j )+ &
defor23(i,k+1,j+1)+defor23(i,k,j+1))
def2(i,k,j)=def2(i,k,j)+0.5*tmp*tmp
enddo
enddo
enddo
!
cr_len = dx + 1. ! hardwire for mixing length = (dx*dy*dz)**(1/3).
! remove this for the alternate formulation
IF (dx .gt. cr_len) THEN
mlen_h=sqrt(dx*dy)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
! mlen_v=dnw(k)/zeta_z(i,j)
mlen_v= 1./rdzw(i,k,j)
tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
tmp=tmp**0.5
xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
xkmhd(i,k,j)=xkmh(i,k,j)
xkmv(i,k,j)=max(c_s*c_s*mlen_v*mlen_v*tmp, 1.0E-6*mlen_v*mlen_v )
xkmv(i,k,j)=min(xkmv(i,k,j), 10.*mlen_v )
xkhh(i,k,j)=xkmh(i,k,j)/pr
xkhv(i,k,j)=xkmv(i,k,j)/pr
ENDDO
ENDDO
ENDDO
ELSE
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
! deltas=(dx*dy*dnw(k)/zeta_z(i,j))**0.33333333
deltas=(dx*dy/rdzw(i,k,j))**0.33333333
tmp=max(0.,def2(i,k,j)-BN2(i,k,j)/pr)
tmp=tmp**0.5
xkmh(i,k,j)=max(c_s*c_s*deltas*deltas*tmp, 1.0E-6*deltas*deltas )
xkmh(i,k,j)=min(xkmh(i,k,j), 10.*deltas )
xkmhd(i,k,j)=xkmh(i,k,j)
xkmv(i,k,j)=xkmh(i,k,j)
xkhh(i,k,j)=xkmh(i,k,j)/pr
xkhv(i,k,j)=xkmv(i,k,j)/pr
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE smag_km
!===============================================================================
SUBROUTINE smag2d_km( config_flags,xkmh,xkmhd,xkmv,xkhh,xkhv, & 1
defor11,defor22,defor12, &
rdzw,dx,dy, &
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 , INTENT(IN ) :: dx, dy
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: rdzw
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
xkmhd, &
xkmv, &
xkhh, &
xkhv
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ), INTENT(IN ) :: &
defor11, &
defor22, &
defor12
! LOCAL VARS
INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
REAL :: deltas, tmp, pr, mlen_h
REAL, DIMENSION( its:ite , kts:kte , jts:jte ) :: def2
!-------------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
pr=1./3.
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
def2(i,k,j)=0.25*(defor11(i,k,j)-defor22(i,k,j))**2 + &
defor12(i,k,j)*defor12(i,k,j)
enddo
enddo
enddo
!
mlen_h=sqrt(dx*dy)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tmp=def2(i,k,j)**0.5
! xkmh(i,k,j)=max(c_s*c_s*mlen_h*mlen_h*tmp, 1.0E-6*mlen_h*mlen_h )
xkmh(i,k,j)=c_s*c_s*mlen_h*mlen_h*tmp
xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
xkmhd(i,k,j)=xkmh(i,k,j)
xkmv(i,k,j)=0.
xkhh(i,k,j)=xkmh(i,k,j)/pr
xkhv(i,k,j)=0.
ENDDO
ENDDO
ENDDO
END SUBROUTINE smag2d_km
!===============================================================================
SUBROUTINE tke_km( config_flags,xkmh,xkmhd,xkmv,xkhh,xkhv,tke,p8w,t8w,theta, & 1
rdz, rdzw, dx,dy,cr_len_in, &
kh_tke_upper_bound, kv_tke_upper_bound, &
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 , INTENT(IN ) :: cr_len_in, dx, dy
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(IN ) :: tke, &
p8w, &
t8w, &
theta, &
rdz, &
rdzw
REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , INTENT(INOUT) :: xkmh, &
xkmhd, &
xkmv, &
xkhh, &
xkhv
REAL, INTENT( IN) :: kh_tke_upper_bound, kv_tke_upper_bound
! LOCAL VARS
REAL , DIMENSION( its:ite, kts:kte, jts:jte) ::dthrdn
REAL :: deltas, tmp, mlen_s, mlen_h, mlen_v, mlen, tmpdz, thetasfc, &
thetatop, minkx, pr, pr_h, pr_v, cr_len
INTEGER :: i_start, i_end, j_start, j_end, ktf, i, j, k
! REAL , DIMENSION( its:ite ,kts:kte, jts:jte) :: pr
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
cr_len = cr_len_in
DO j = j_start, j_end
DO k = kts+1, ktf-1
DO i = i_start, i_end
! tmpdz=dn(k+1)+dn(k)
tmpdz=1./(rdz(i,k+1,j)+rdz(i,k,j))
dthrdn(i,k,j)=(theta(i,k+1,j)-theta(i,k-1,j))/tmpdz
ENDDO
ENDDO
ENDDO
k=kts
DO j = j_start, j_end
DO i = i_start, i_end
! tmpdz=dn(k+1)+0.5*dnw(k)
tmpdz=1./(rdzw(i,k+1,j)+rdzw(i,k,j))
thetasfc=T8w(i,kts,j)/(p8w(i,k,j)/p1000mb)**(R_d/Cp)
dthrdn(i,k,j)=(theta(i,k+1,j)-thetasfc)/tmpdz
ENDDO
ENDDO
k=ktf
DO j = j_start, j_end
DO i = i_start, i_end
! tmpdz=dn(k)+0.5*dnw(k)
tmpdz=1./rdz(i,k,j) + 0.5/rdzw(i,k,j)
thetatop=T8w(i,kde,j)/(p8w(i,kde,j)/p1000mb)**(R_d/Cp)
dthrdn(i,k,j)=(thetatop-theta(i,k-1,j))/tmpdz
ENDDO
ENDDO
cr_len = dx + 1. ! hardwire for mixing length = (dx*dy*dz)**(1/3).
! remove this for the alternate formulation
IF (dx .gt. cr_len) THEN
mlen_h=sqrt(dx*dy)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tmp=sqrt(max(tke(i,k,j),0.0))
! deltas=dnw(k)/zeta_z(i,j)
deltas=1./rdzw(i,k,j)
mlen_v=deltas
IF (dthrdn(i,k,j) .gt. 0.) THEN
! mlen_s=0.76*tmp/ &
! abs(g/theta(i,k,j)*zeta_z(i,j)*dthrdn(i,k,j))
mlen_s=0.76*tmp/(abs(g/theta(i,k,j)*dthrdn(i,k,j)))**0.5
mlen_v=min(mlen_v,mlen_s)
ENDIF
xkmh(i,k,j)=max(c_k*tmp*mlen_h,1.0E-6*mlen_h*mlen_h)
xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen_h )
xkmhd(i,k,j)=xkmh(i,k,j)
xkmv(i,k,j)=max(c_k*tmp*mlen_v,1.0E-6*deltas*deltas)
xkmv(i,k,j)=min(xkmv(i,k,j), 10.*deltas )
pr_h=1./3.
pr_v=1.+2.*mlen_v/deltas
xkhh(i,k,j)=xkmh(i,k,j)/pr_h
xkhv(i,k,j)=xkmv(i,k,j)*pr_v
ENDDO
ENDDO
ENDDO
ELSE
DO j = j_start, j_end
DO k = kts+1, ktf-1
DO i = i_start, i_end
tmp=sqrt(max(tke(i,k,j),0.0))
! deltas=(dx*dy*dnw(k)/zeta_z(i,j))**0.33333333
deltas=(dx*dy/rdzw(i,k,j))**0.33333333
mlen=deltas
IF (dthrdn(i,k,j) .gt. 0.) THEN
! mlen_s=0.76*tmp/ &
! abs(g/theta(i,k,j)*zeta_z(i,j)*dthrdn(i,k,j))
!! mlen_s=0.76*tmp/(abs(g/theta(i,k,j)*dthrdn(i,k,j)))**0.5
!! mlen=min(mlen,mlen_s)
ENDIF
minkx=1.0E-6*deltas*deltas
xkmh(i,k,j)=max(0.1*tmp*mlen,minkx)
!! xkmh(i,k,j)=min(xkmh(i,k,j), 10.*mlen )
xkmh(i,k,j)=min(kh_tke_upper_bound,min(xkmh(i,k,j), 10.*mlen ))
xkmhd(i,k,j)=xkmh(i,k,j)
xkmv(i,k,j)=max(0.1*tmp*mlen,minkx)
!! xkmv(i,k,j)=min(xkmv(i,k,j), 10.*mlen )
xkmv(i,k,j)=min(kv_tke_upper_bound,min(xkmv(i,k,j), 10.*mlen ))
pr=1.+2.*mlen/deltas
!! xkhh(i,k,j)=xkmh(i,k,j)/pr
!! xkhv(i,k,j)=xkmv(i,k,j)/pr
xkhh(i,k,j)=min(kh_tke_upper_bound,xkmh(i,k,j)*pr)
xkhv(i,k,j)=min(kv_tke_upper_bound,xkmv(i,k,j)*pr)
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE tke_km
!=================================================================================
SUBROUTINE horizontal_diffusion_2 ( rt_tendf, ru_tendf, rv_tendf, rw_tendf, & 1,7
tke_tendf, &
moist_tendf, n_moist, &
scalar_tendf, n_scalar, &
thp, theta, mu, tke, config_flags, &
defor11, defor22, defor12, &
defor13, defor23, div, &
moist, scalar, &
msfu, msfv, msft, xkmh, xkhh,km_opt, &
rdx, rdy, rdz, rdzw, fnm, fnp, &
cf1, cf2, cf3, zx, zy, dn, dnw, &
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, n_scalar, km_opt
REAL , INTENT(IN ) :: cf1, cf2, cf3
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfu, &
msfv, &
msft, &
mu
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::rt_tendf,&
ru_tendf,&
rv_tendf,&
rw_tendf,&
tke_tendf
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
INTENT(INOUT) :: moist_tendf
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar), &
INTENT(INOUT) :: scalar_tendf
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
INTENT(IN ) :: moist
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
INTENT(IN ) :: scalar
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
defor22, &
defor12, &
defor13, &
defor23, &
div, &
xkmh, &
xkhh, &
zx, &
zy, &
theta, &
thp, &
tke, &
rdz, &
rdzw
REAL , INTENT(IN ) :: rdx, &
rdy
! LOCAL VARS
INTEGER :: im
! REAL , DIMENSION(its-1:ite+1, kts:kte, jts-1:jte+1) :: xkhh
!--------------------------------------------------------------------------
CALL horizontal_diffusion_u_2( ru_tendf, mu, config_flags, &
defor11, defor12, div, &
tke(ims,kms,jms), &
msfu, xkmh, rdx, rdy, fnm, fnp, &
zx, zy, rdzw, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
CALL horizontal_diffusion_v_2( rv_tendf, mu, config_flags, &
defor12, defor22, div, &
tke(ims,kms,jms), &
msfv, xkmh, rdx, rdy, fnm, fnp, &
zx, zy, rdzw, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
CALL horizontal_diffusion_w_2( rw_tendf, mu, config_flags, &
defor13, defor23, div, &
tke(ims,kms,jms), &
msft, xkmh, rdx, rdy, fnm, fnp, &
zx, zy, rdz, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! calculate khh
CALL horizontal_diffusion_s ( rt_tendf, mu, config_flags, thp, &
msft, msfu, msfv, xkhh, rdx, rdy, &
fnm, fnp, cf1, cf2, cf3, &
zx, zy, rdz, rdzw, dnw, dn, &
.false., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
IF (km_opt .eq. 2) &
CALL horizontal_diffusion_s ( tke_tendf(ims,kms,jms), &
mu, config_flags, &
tke(ims,kms,jms), &
msft, msfu, msfv, xkhh, rdx, rdy, &
fnm, fnp, cf1, cf2, cf3, &
zx, zy, rdz, rdzw, dnw, dn, &
.true., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
CALL horizontal_diffusion_s( moist_tendf(ims,kms,jms,im), &
mu, config_flags, &
moist(ims,kms,jms,im), &
msft, msfu, msfv, xkhh, rdx, rdy, &
fnm, fnp, cf1, cf2, cf3, &
zx, zy, rdz, rdzw, dnw, dn, &
.false., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
ENDDO moist_loop
ENDIF
IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar
CALL horizontal_diffusion_s( scalar_tendf(ims,kms,jms,im), &
mu, config_flags, &
scalar(ims,kms,jms,im), &
msft, msfu, msfv, xkhh, rdx, rdy, &
fnm, fnp, cf1, cf2, cf3, &
zx, zy, rdz, rdzw, dnw, dn, &
.false., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
ENDDO scalar_loop
ENDIF
END SUBROUTINE horizontal_diffusion_2
!==========================================================================
SUBROUTINE horizontal_diffusion_u_2( tendency, mu, config_flags, & 1,2
defor11, defor12, div, tke, &
msfu, xkmh, rdx, rdy, fnm, fnp, &
zx, zy, rdzw, &
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( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfu, &
mu
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: rdzw
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
defor12, &
div, &
tke, &
xkmh, &
zx, &
zy
REAL , INTENT(IN ) :: rdx, &
rdy
! Local data
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
titau2avg, &
titau1, &
titau2, &
xkxavg, &
rravg
! new
! zxavg, &
! zyavg
REAL :: mrdx, mrdy, rcoup
! new
REAL :: tmpzy, tmpzeta_z
REAL :: term1, term2, term3
ktf=MIN(kte,kde-1)
!----------------------------------------------------------------------
! u : p (.), u(|), w(-)
!
! p u p u u u
!
! p | . | . | . | k+1 | . | . | . | k+1
!
! w - 13 - - k+1 13 k+1
!
! p | 11 O 11 | . | k | 12 O 12 | . | k
!
! w - 13 - - k 13 k
!
! p | . | . | . | k-1 | . | . | . | k-1
!
! i-1 i i i+1 j-1 j j j+1 j+1
!
i_start = its
i_end = ite
j_start = jts
j_end = MIN(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-1,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! titau1 = titau11
is_ext=1
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_11_22_33( config_flags,titau1,mu,tke, &
xkmh,defor11,div, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau2 = titau12
is_ext=0
ie_ext=0
js_ext=0
je_ext=1
CALL cal_titau_12_21( config_flags,titau2,mu,xkmh,defor12, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau1avg = titau11avg
! titau2avg = titau12avg
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i-1,k ,j)+titau1(i,k ,j))+ &
fnp(k)*(titau1(i-1,k-1,j)+titau1(i,k-1,j)))
titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k ,j+1)+titau2(i,k ,j))+ &
fnp(k)*(titau2(i,k-1,j+1)+titau2(i,k-1,j)))
tmpzy = 0.25*( zy(i-1,k,j )+zy(i,k,j )+ &
zy(i-1,k,j+1)+zy(i,k,j+1) )
! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
! titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)*tmpzeta_z
! titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy *tmpzeta_z
titau1avg(i,k,j)=titau1avg(i,k,j)*zx(i,k,j)
titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
ENDDO
ENDDO
ENDDO
!
DO j = j_start, j_end
DO i = i_start, i_end
titau1avg(i,kts,j)=0.
titau1avg(i,ktf+1,j)=0.
titau2avg(i,kts,j)=0.
titau2avg(i,ktf+1,j)=0.
ENDDO
ENDDO
!
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
mrdx=msfu(i,j)*rdx
mrdy=msfu(i,j)*rdy
tendency(i,k,j)=tendency(i,k,j)- &
(mrdx*(titau1(i,k,j )-titau1(i-1,k,j))+ &
mrdy*(titau2(i,k,j+1)-titau2(i,k,j ))- &
msfu(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
(titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
) )
ENDDO
ENDDO
ENDDO
END SUBROUTINE horizontal_diffusion_u_2
!==========================================================================
SUBROUTINE horizontal_diffusion_v_2( tendency, mu, config_flags, & 1,2
defor12, defor22, div, tke, &
msfv, xkmh, rdx, rdy, fnm, fnp, &
zx, zy, rdzw, &
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( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfv, &
mu
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor12, &
defor22, &
div, &
tke, &
xkmh, &
zx, &
zy, &
rdzw
REAL , INTENT(IN ) :: rdx, &
rdy
! Local data
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
titau2avg, &
titau1, &
titau2, &
xkxavg, &
rravg
! new
! zxavg, &
! zyavg
REAL :: mrdx, mrdy, rcoup
! new
REAL :: tmpzx, tmpzeta_z
ktf=MIN(kte,kde-1)
!----------------------------------------------------------------------
! v : p (.), v(+), w(-)
!
! p v p v v v
!
! p + . + . + . + k+1 + . + . + . + k+1
!
! w - 23 - - k+1 23 k+1
!
! p + 22 O 22 + . + k + 21 O 21 + . + k
!
! w - 23 - - k 23 k
!
! p + . + . + . + k-1 + . + . + . + k-1
!
! j-1 j j j+1 i-1 i i i+1 i+1
!
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = jte
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-1,jte)
! titau1 = titau21
is_ext=0
ie_ext=1
js_ext=0
je_ext=0
CALL cal_titau_12_21( config_flags,titau1,mu,xkmh,defor12, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau2 = titau22
is_ext=0
ie_ext=0
js_ext=1
je_ext=0
CALL cal_titau_11_22_33( config_flags,titau2,mu,tke, &
xkmh,defor22,div, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
titau1avg(i,k,j)=0.5*(fnm(k)*(titau1(i+1,k ,j)+titau1(i,k ,j))+ &
fnp(k)*(titau1(i+1,k-1,j)+titau1(i,k-1,j)))
titau2avg(i,k,j)=0.5*(fnm(k)*(titau2(i,k ,j-1)+titau2(i,k ,j))+ &
fnp(k)*(titau2(i,k-1,j-1)+titau2(i,k-1,j)))
tmpzx = 0.25*( zx(i,k,j )+zx(i+1,k,j )+ &
zx(i,k,j-1)+zx(i+1,k,j-1) )
titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
titau2avg(i,k,j)=titau2avg(i,k,j)*zy(i,k,j)
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO i = i_start, i_end
titau1avg(i,kts,j)=0.
titau1avg(i,ktf+1,j)=0.
titau2avg(i,kts,j)=0.
titau2avg(i,ktf+1,j)=0.
ENDDO
ENDDO
!
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
mrdx=msfv(i,j)*rdx
mrdy=msfv(i,j)*rdy
tendency(i,k,j)=tendency(i,k,j)- &
(mrdy*(titau2(i ,k,j)-titau2(i,k,j-1))+ &
mrdx*(titau1(i+1,k,j)-titau1(i,k,j ))- &
msfv(i,j)*rdzw(i,k,j)*((titau1avg(i,k+1,j)-titau1avg(i,k,j))+ &
(titau2avg(i,k+1,j)-titau2avg(i,k,j)) &
) &
)
ENDDO
ENDDO
ENDDO
END SUBROUTINE horizontal_diffusion_v_2
!==========================================================================
SUBROUTINE horizontal_diffusion_w_2( tendency, mu, config_flags, & 1,2
defor13, defor23, div, tke, &
msft, xkmh, rdx, rdy, fnm, fnp, &
zx, zy, rdz, &
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( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msft, &
mu
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
defor23, &
div, &
tke, &
xkmh, &
zx, &
zy, &
rdz
REAL , INTENT(IN ) :: rdx, &
rdy
! Local data
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau1avg, &
titau2avg, &
titau1, &
titau2, &
xkxavg, &
rravg
! new
! zxavg, &
! zyavg
REAL :: mrdx, mrdy, rcoup
! new
REAL :: tmpzx, tmpzy, tmpzeta_z
ktf=MIN(kte,kde-1)
!----------------------------------------------------------------------
! w : p (.), u(|), v(+), w(-)
!
! p u p u p v p v
!
! w - - - k+1 w - - - k+1
!
! p . | 33 | . k p . + 33 + . k
!
! w - 31 O 31 - k w - 32 O 32 - k
!
! p . | 33 | . k-1 p . | 33 | . k-1
!
! w - - - k-1 w - - - k-1
!
! i-1 i i i+1 j-1 j j j+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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! titau1 = titau31
is_ext=0
ie_ext=1
js_ext=0
je_ext=0
CALL cal_titau_13_31( config_flags,titau1,defor13, &
mu,xkmh,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau2 = titau32
is_ext=0
ie_ext=0
js_ext=0
je_ext=1
CALL cal_titau_23_32( config_flags,titau2,defor23, &
mu,xkmh,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau1avg = titau31avg * zx * zeta_z = titau13avg * zx * zeta_z
! titau2avg = titau32avg * zy * zeta_z = titau23avg * zy * zeta_z
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
titau1avg(i,k,j)=0.25*(titau1(i+1,k+1,j)+titau1(i,k+1,j)+ &
titau1(i+1,k ,j)+titau1(i,k ,j))
titau2avg(i,k,j)=0.25*(titau2(i,k+1,j+1)+titau2(i,k+1,j)+ &
titau2(i,k ,j+1)+titau2(i,k ,j))
! new
tmpzx =0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ &
zx(i,k+1,j)+zx(i+1,k+1,j) )
tmpzy =0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ &
zy(i,k+1,j)+zy(i,k+1,j+1) )
titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx
titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy
! titau1avg(i,k,j)=titau1avg(i,k,j)*tmpzx*zeta_z(i,j)
! titau2avg(i,k,j)=titau2avg(i,k,j)*tmpzy*zeta_z(i,j)
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO i = i_start, i_end
titau1avg(i,kts ,j)=0.
titau2avg(i,kts ,j)=0.
titau1avg(i,ktf+1,j)=0.
titau2avg(i,ktf+1,j)=0.
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
mrdx=msft(i,j)*rdx
mrdy=msft(i,j)*rdy
tendency(i,k,j)=tendency(i,k,j)- &
(mrdx*(titau1(i+1,k,j)-titau1(i,k,j))+ &
mrdy*(titau2(i,k,j+1)-titau2(i,k,j))- &
msft(i,j)*rdz(i,k,j)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
titau2avg(i,k,j)-titau2avg(i,k-1,j) &
) &
)
! msft(i,j)/dn(k)*(titau1avg(i,k,j)-titau1avg(i,k-1,j)+ &
! titau2avg(i,k,j)-titau2avg(i,k-1,j) &
! ) &
! )
ENDDO
ENDDO
ENDDO
END SUBROUTINE horizontal_diffusion_w_2
!==========================================================================
SUBROUTINE horizontal_diffusion_s (tendency, mu, config_flags, var, & 4
msft, msfu, msfv, xkhh, rdx, rdy, &
fnm, fnp, cf1, cf2, cf3, &
zx, zy, rdz, rdzw, dn, dnw, &
doing_tke, &
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
LOGICAL, INTENT(IN ) :: doing_tke
REAL , INTENT(IN ) :: cf1, cf2, cf3
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfu
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msfv
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msft
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: mu
! REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
! INTENT(IN ) :: xkhh
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: &
xkhh, &
rdz, &
rdzw
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: var, &
zx, &
zy
REAL , INTENT(IN ) :: rdx, &
rdy
! Local data
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: H1avg, &
H2avg, &
H1, &
H2, &
xkxavg
! new
! zxavg, &
! zyavg
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
REAL :: mrdx, mrdy, rcoup
! new
REAL :: tmpzx, tmpzy, tmpzeta_z
INTEGER :: ktes1,ktes2
!-------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
!------------------------------------------------------------------------
! scalars: t (.), u(|), v(+), w(-)
!
! t u t u t v t v
!
! w - 3 - k+1 w - 3 - k+1
!
! t . 1 O 1 . k t . 2 O 2 . k
!
! w - 3 - k w - 3 - k
!
! t . | . | . k-1 t . + . + . k-1
!
! w - - - k-1 w - - - k-1
!
! t i-1 i i i+1 j-1 j j j+1
!
ktes1=kte-1
ktes2=kte-2
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! diffusion of the TKE needs mutiple 2
IF ( doing_tke ) THEN
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
tmptendf(i,k,j)=tendency(i,k,j)
ENDDO
ENDDO
ENDDO
ENDIF
! H1 = partial var over partial x
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end + 1
! new
! zxavg(i,k,j) =0.5*( zx(i-1,k,j)+ zx(i,k,j))
xkxavg(i,k,j)=0.5*(xkhh(i-1,k,j)+xkhh(i,k,j))
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts+1, ktf
DO i = i_start, i_end + 1
H1avg(i,k,j)=0.5*(fnm(k)*(var(i-1,k ,j)+var(i,k ,j))+ &
fnp(k)*(var(i-1,k-1,j)+var(i,k-1,j)))
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO i = i_start, i_end + 1
H1avg(i,kts ,j)=0.5*(cf1*var(i ,1,j)+cf2*var(i ,2,j)+ &
cf3*var(i ,3,j)+cf1*var(i-1,1,j)+ &
cf2*var(i-1,2,j)+cf3*var(i-1,3,j))
H1avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
var(i-1,ktes1,j)+(var(i-1,ktes1,j)- &
var(i-1,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1))
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end + 1
! new
tmpzx = 0.5*( zx(i,k,j)+ zx(i,k+1,j))
H1(i,k,j)=-msfu(i,j)*xkxavg(i,k,j)*( &
rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx* &
(H1avg(i,k+1,j)-H1avg(i,k,j))*rdzw(i,k,j) )
! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i-1,j))
! H1(i,k,j)=-msfu(i,j)*xkxavg(i,k,j)*( &
! rdx*(var(i,k,j)-var(i-1,k,j)) - tmpzx*tmpzeta_z* &
! (H1avg(i,k+1,j)-H1avg(i,k,j))/dnw(k))
ENDDO
ENDDO
ENDDO
! H2 = partial var over partial y
DO j = j_start, j_end + 1
DO k = kts, ktf
DO i = i_start, i_end
! new
! zyavg(i,k,j) =0.5*( zy(i,k,j-1)+ zy(i,k,j))
xkxavg(i,k,j)=0.5*(xkhh(i,k,j-1)+xkhh(i,k,j))
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end + 1
DO k = kts+1, ktf
DO i = i_start, i_end
! new
H2avg(i,k,j)=0.5*(fnm(k)*(var(i,k ,j-1)+var(i,k ,j))+ &
fnp(k)*(var(i,k-1,j-1)+var(i,k-1,j)))
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end + 1
DO i = i_start, i_end
H2avg(i,kts ,j)=0.5*(cf1*var(i,1,j )+cf2*var(i ,2,j)+ &
cf3*var(i,3,j )+cf1*var(i,1,j-1)+ &
cf2*var(i,2,j-1)+cf3*var(i,3,j-1))
H2avg(i,ktf+1,j)=0.5*(var(i,ktes1,j)+(var(i,ktes1,j)- &
var(i,ktes2,j))*0.5*dnw(ktes1)/dn(ktes1)+ &
var(i,ktes1,j-1)+(var(i,ktes1,j-1)- &
var(i,ktes2,j-1))*0.5*dnw(ktes1)/dn(ktes1))
ENDDO
ENDDO
DO j = j_start, j_end + 1
DO k = kts, ktf
DO i = i_start, i_end
! new
tmpzy = 0.5*( zy(i,k,j)+ zy(i,k+1,j))
H2(i,k,j)=-msfv(i,j)*xkxavg(i,k,j)*( &
rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy* &
(H2avg(i ,k+1,j)-H2avg(i,k,j))*rdzw(i,k,j))
! tmpzeta_z = 0.5*(zeta_z(i,j)+zeta_z(i,j-1))
! H2(i,k,j)=-msfv(i,j)*xkxavg(i,k,j)*( &
! rdy*(var(i,k,j)-var(i,k,j-1)) - tmpzy*tmpzeta_z* &
! (H2avg(i ,k+1,j)-H2avg(i,k,j))/dnw(k))
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts+1, ktf
DO i = i_start, i_end
H1avg(i,k,j)=0.5*(fnm(k)*(H1(i+1,k ,j)+H1(i,k ,j))+ &
fnp(k)*(H1(i+1,k-1,j)+H1(i,k-1,j)))
H2avg(i,k,j)=0.5*(fnm(k)*(H2(i,k ,j+1)+H2(i,k ,j))+ &
fnp(k)*(H2(i,k-1,j+1)+H2(i,k-1,j)))
! new
! zxavg(i,k,j)=fnm(k)*zx(i,k,j)+fnp(k)*zx(i,k-1,j)
! zyavg(i,k,j)=fnm(k)*zy(i,k,j)+fnp(k)*zy(i,k-1,j)
! H1avg(i,k,j)=zx*H1avg*zeta_z
! H2avg(i,k,j)=zy*H2avg*zeta_z
tmpzx = 0.5*( zx(i,k,j)+ zx(i+1,k,j ))
tmpzy = 0.5*( zy(i,k,j)+ zy(i ,k,j+1))
H1avg(i,k,j)=H1avg(i,k,j)*tmpzx
H2avg(i,k,j)=H2avg(i,k,j)*tmpzy
! H1avg(i,k,j)=H1avg(i,k,j)*tmpzx*zeta_z(i,j)
! H2avg(i,k,j)=H2avg(i,k,j)*tmpzy*zeta_z(i,j)
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO i = i_start, i_end
H1avg(i,kts ,j)=0.
H1avg(i,ktf+1,j)=0.
H2avg(i,kts ,j)=0.
H2avg(i,ktf+1,j)=0.
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
mrdx=msft(i,j)*rdx
mrdy=msft(i,j)*rdy
tendency(i,k,j)=tendency(i,k,j)- &
(mrdx*0.5*((mu(i+1,j)+mu(i,j))*H1(i+1,k,j)- &
(mu(i-1,j)+mu(i,j))*H1(i ,k,j))+ &
mrdy*0.5*((mu(i,j+1)+mu(i,j))*H2(i,k,j+1)- &
(mu(i,j-1)+mu(i,j))*H2(i,k,j ))- &
msft(i,j)*(H1avg(i,k+1,j)-H1avg(i,k,j)+ &
H2avg(i,k+1,j)-H2avg(i,k,j) &
)*rdzw(i,k,j) &
)
ENDDO
ENDDO
ENDDO
IF ( doing_tke ) THEN
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
tendency(i,k,j)=tmptendf(i,k,j)+2.* &
(tendency(i,k,j)-tmptendf(i,k,j))
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE horizontal_diffusion_s
!==============================================================================
SUBROUTINE vertical_diffusion_2 ( ru_tendf, rv_tendf, rw_tendf, rt_tendf, & 1,7
tke_tendf, moist_tendf, n_moist, &
scalar_tendf, n_scalar, &
u_2, v_2, &
thp,u_base,v_base,t_base,qv_base,mu,tke, &
config_flags,defor13,defor23,defor33,div, &
moist, scalar, xkmv, xkhv,km_opt, &
fnm, fnp, dn, dnw, rdz, rdzw, &
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, n_scalar, km_opt
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: qv_base
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: u_base
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: v_base
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: t_base
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::ru_tendf,&
rv_tendf,&
rw_tendf,&
tke_tendf,&
rt_tendf
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
INTENT(INOUT) :: moist_tendf
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
INTENT(INOUT) :: scalar_tendf
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_moist), &
INTENT(INOUT) :: moist
REAL , DIMENSION( ims:ime, kms:kme, jms:jme, n_scalar) , &
INTENT(IN ) :: scalar
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor13, &
defor23, &
defor33, &
div, &
xkmv, &
xkhv, &
tke, &
rdz, &
u_2, &
v_2, &
rdzw
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: thp
! LOCAL VAR
LOGICAL, PARAMETER :: filter_perturbations = .true.
REAL , DIMENSION( ims:ime, kms:kme, jms:jme) :: var_mix
INTEGER :: im, i,j,k
INTEGER :: i_start, i_end, j_start, j_end
! REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: xkhv
!***************************************************************************
!***************************************************************************
!MODIFICA VARIABILI PER I FLUSSI
!
! REAL , DIMENSION( its:ite, jts:jte) :: Cd
! REAL :: V0_u,V0_v,tao_xz,tao_yz,ustar
!
!FINE MODIFICA VARIABILI PER I FLUSSI
!***************************************************************************
!***************************************************************************
!
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
!
!--------------------------------------------------------------------------
CALL vertical_diffusion_u_2( ru_tendf, config_flags, mu, &
defor13, xkmv, &
dnw, rdzw, fnm, fnp, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
CALL vertical_diffusion_v_2( rv_tendf, config_flags, mu, &
defor23, xkmv, &
dnw, rdzw, fnm, fnp, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
CALL vertical_diffusion_w_2( rw_tendf, config_flags, mu, &
defor33, tke(ims,kms,jms), &
div, xkmv, &
dn, rdz, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
!*****************************************
!*****************************************
! MODIFICA al flusso di momento alla parete
!
! DO j = j_start, j_end+1
! DO i = i_start, i_end+1
! Cd(i,j)=0.007
! ENDDO
! ENDDO
!
!calcolo del modulo della velocita
! DO j = j_start, j_end
! DO i = i_start, i_end+1
! V0_u=0.
! tao_xz=0.
! V0_u= sqrt((u_2(i,kts,j)**2) + &
! (((v_2(i ,kts,j )+ &
! v_2(i ,kts,j+1)+ &
! v_2(i-1,kts,j )+ &
! v_2(i-1,kts,j+1))/4)**2))+epsilon
! tao_xz=Cd(i,j)*V0_u*u_2(i,kts,j)
! ru_tendf(i,kts,j)=ru_tendf(i,kts,j) &
! -0.25*(mu(i,j)+mu(i-1,j))*tao_xz*(rdzw(i,kts,j)+rdzw(i-1,kts,j))
! ENDDO
! ENDDO
!
!
! DO j = j_start, j_end+1
! DO i = i_start, i_end
! V0_v=0.
! tao_yz=0.
! V0_v= sqrt((v_2(i,kts,j)**2) + &
! (((u_2(i ,kts,j )+ &
! u_2(i ,kts,j-1)+ &
! u_2(i+1,kts,j )+ &
! u_2(i+1,kts,j-1))/4)**2))+epsilon
! tao_yz=Cd(i,j)*V0_v*v_2(i,kts,j)
! rv_tendf(i,kts,j)=rv_tendf(i,kts,j) &
! -0.25*(mu(i,j)+mu(i,j-1))*tao_yz*(rdzw(i,kts,j)+rdzw(i,kts,j-1))
! ENDDO
! ENDDO
!
! FINE MODIFICA al flusso di momento alla parete
!*****************************************
!*****************************************
IF (filter_perturbations) THEN
DO j=jts,min(jte,jde-1)
DO k=kts,kte-1
DO i=its,min(ite,ide-1)
var_mix(i,k,j) = thp(i,k,j) - t_base(k)
ENDDO
ENDDO
ENDDO
ELSE
DO j=jts,min(jte,jde-1)
DO k=kts,kte-1
DO i=its,min(ite,ide-1)
var_mix(i,k,j) = thp(i,k,j)
ENDDO
ENDDO
ENDDO
END IF
CALL vertical_diffusion_s( rt_tendf, config_flags, var_mix, mu, xkhv, &
dn, dnw, rdz, rdzw, fnm, fnp, &
.false., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
!*****************************************
!*****************************************
!MODIFICA al flusso di calore
!
!
! DO j = j_start, j_end
! DO i = i_start, i_end
!
! rt_tendf(i,kts,j)=rt_tendf(i,kts,j)+mu(i,j)*0.06*rdzw(i,kts,j)
!
! ENDDO
! ENDDO
!
! FINE MODIFICA al flusso di calore
!*****************************************
!*****************************************
If (km_opt .eq. 2) then
CALL vertical_diffusion_s( tke_tendf(ims,kms,jms), &
config_flags, tke(ims,kms,jms), &
mu, xkhv, &
dn, dnw, rdz, rdzw, fnm, fnp, &
.true., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
endif
IF (n_moist .ge. PARAM_FIRST_SCALAR) THEN
moist_loop: do im = PARAM_FIRST_SCALAR, n_moist
IF (filter_perturbations .and. (im == P_QV)) THEN
DO j=jts,min(jte,jde-1)
DO k=kts,kte-1
DO i=its,min(ite,ide-1)
var_mix(i,k,j) = moist(i,k,j,im) - qv_base(k)
ENDDO
ENDDO
ENDDO
ELSE
DO j=jts,min(jte,jde-1)
DO k=kts,kte-1
DO i=its,min(ite,ide-1)
var_mix(i,k,j) = moist(i,k,j,im)
ENDDO
ENDDO
ENDDO
END IF
CALL vertical_diffusion_s( moist_tendf(ims,kms,jms,im), &
config_flags, var_mix, &
mu, xkhv, &
dn, dnw, rdz, rdzw, fnm, fnp, &
.false., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
ENDDO moist_loop
ENDIF
IF (n_scalar .ge. PARAM_FIRST_SCALAR) THEN
scalar_loop: do im = PARAM_FIRST_SCALAR, n_scalar
CALL vertical_diffusion_s( scalar_tendf(ims,kms,jms,im), &
config_flags, scalar(ims,kms,jms,im), &
mu, xkhv, &
dn, dnw, rdz, rdzw, fnm, fnp, &
.false., &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
ENDDO scalar_loop
ENDIF
END SUBROUTINE vertical_diffusion_2
!==========================================================================
SUBROUTINE vertical_diffusion_u_2( tendency, config_flags, mu, & 1,1
defor13, xkmv, &
dnw, rdzw, fnm, fnp, &
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( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
INTENT(IN ) ::defor13, &
xkmv, &
rdzw
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
! LOCAL VARS
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
REAL , DIMENSION( its:ite, jts:jte) :: zzavg
REAL :: rdzu
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-1,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! titau3 = titau13
is_ext=0
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_13_31( config_flags,titau3,defor13,mu,xkmv,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
!
DO j = j_start, j_end
DO k=kts+1,ktf
DO i = i_start, i_end
rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j)-titau3(i,k,j))
ENDDO
ENDDO
ENDDO
! ******** MODIF...
! we will pick up the surface drag (titau3(i,kts,j)) later
!
! DO j = j_start, j_end
! k=kts
! DO i = i_start, i_end
!
! rdzu = 2./(1./rdzw(i,k,j) + 1./rdzw(i-1,k,j))
! tendency(i,k,j)=tendency(i,k,j)-rdzu*(titau3(i,k+1,j))
! ENDDO
! ENDDO
! ******** MODIF...
END SUBROUTINE vertical_diffusion_u_2
!==========================================================================
SUBROUTINE vertical_diffusion_v_2( tendency, config_flags, mu, & 1,1
defor23, xkmv, &
dnw, rdzw, fnm, fnp, &
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( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
! REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: zeta_z
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
INTENT(IN ) ::defor23, &
xkmv, &
rdzw
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN ) :: mu
! LOCAL VARS
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
REAL , DIMENSION( its:ite, jts:jte) :: zzavg
REAL :: rdzv
!--------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = jte
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-1,jte)
! titau3 = titau23
is_ext=0
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_23_32( config_flags,titau3,defor23,mu,xkmv,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j)-titau3(i,k,j))
ENDDO
ENDDO
ENDDO
! ******** MODIF...
! we will pick up the surface drag (titau3(i,kts,j)) later
!
! DO j = j_start, j_end
! k=kts
! DO i = i_start, i_end
!
! rdzv = 2./(1./rdzw(i,k,j) + 1./rdzw(i,k,j-1))
! tendency(i,k,j)=tendency(i,k,j)-rdzv*(titau3(i,k+1,j))
!
! ENDDO
! ENDDO
! ******** MODIF...
END SUBROUTINE vertical_diffusion_v_2
!==========================================================================
SUBROUTINE vertical_diffusion_w_2(tendency, config_flags, mu, & 1,1
defor33, tke, div, xkmv, &
dn, rdz, &
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( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
INTENT(IN ) ::defor33, &
tke, &
div, &
xkmv, &
rdz
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! LOCAL VARS
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: titau3
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! titau3 = titau33
is_ext=0
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_11_22_33( config_flags,titau3,mu,tke, &
xkmv,defor33,div, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! DO j = j_start, j_end
! DO k = kts+1, ktf
! DO i = i_start, i_end
! titau3(i,k,j)=titau3(i,k,j)*zeta_z(i,j)
! ENDDO
! ENDDO
! ENDDO
DO j = j_start, j_end
DO k = kts+1, ktf
DO i = i_start, i_end
tendency(i,k,j)=tendency(i,k,j)-rdz(i,k,j)*(titau3(i,k,j)-titau3(i,k-1,j))
ENDDO
ENDDO
ENDDO
END SUBROUTINE vertical_diffusion_w_2
!==========================================================================
SUBROUTINE vertical_diffusion_s( tendency, config_flags, var, mu, xkhv, & 4
dn, dnw, rdz, rdzw, fnm, fnp, &
doing_tke, &
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
LOGICAL, INTENT(IN ) :: doing_tke
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
INTENT(IN ) :: var, &
rdz, &
rdzw
! LOCAL VARS
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: H3, &
xkxavg, &
rravg
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: tmptendf
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
IF (doing_tke) THEN
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
tmptendf(i,k,j)=tendency(i,k,j)
ENDDO
ENDDO
ENDDO
ENDIF
! H3
xkxavg = 0.
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
xkxavg(i,k,j)=fnm(k)*xkhv(i,k,j)+fnp(k)*xkhv(i,k-1,j)
H3(i,k,j)=-xkxavg(i,k,j)*(var(i,k,j)-var(i,k-1,j))*rdz(i,k,j)
! H3(i,k,j)=-xkxavg(i,k,j)*zeta_z(i,j)* &
! (var(i,k,j)-var(i,k-1,j))/dn(k)
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO i = i_start, i_end
H3(i,kts,j)=0.
! H3(i,ktf+1,j)=0.
! H3(i,kts,j)=H3(i,kts+1,j)
H3(i,ktf+1,j)=H3(i,ktf,j)
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
tendency(i,k,j)=tendency(i,k,j) &
-mu(i,j)*(H3(i,k+1,j)-H3(i,k,j))*rdzw(i,k,j)
ENDDO
ENDDO
ENDDO
IF (doing_tke) THEN
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
tendency(i,k,j)=tmptendf(i,k,j)+2.* &
(tendency(i,k,j)-tmptendf(i,k,j))
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE vertical_diffusion_s
!==========================================================================
SUBROUTINE vertical_diffusion_s_qv( tendency, config_flags, var, &
var_base, mu, xkhv, &
rdz, rdzw, fnm, fnp, &
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( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) ::var_base
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , INTENT(IN) :: xkhv
REAL , DIMENSION( ims:ime , jms:jme ) , INTENT(IN) :: mu
REAL , DIMENSION( ims:ime , kms:kme, jms:jme ) , &
INTENT(IN ) :: var, &
rdz, &
rdzw
! LOCAL VARS
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its:ite, kts:kte+1, jts:jte) :: H3, &
xkxavg
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! H3
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
xkxavg(i,k,j)=fnm(k)*xkhv(i,k,j)+fnp(k)*xkhv(i,k-1,j)
! H3(i,k,j)=-xkxavg(i,k,j)*zeta_z(i,j)* &
! (var(i,k ,j)-var_base(k )- &
! var(i,k-1,j)+var_base(k-1))/dn(k)
H3(i,k,j)=-xkxavg(i,k,j)* &
(var(i,k ,j)-var_base(k )- &
var(i,k-1,j)+var_base(k-1))*rdz(i,k,j)
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO i = i_start, i_end
H3(i,kts,j)=0.
H3(i,ktf+1,j)=0.
ENDDO
ENDDO
! H3 = rr * zeta_z * H3
! DO j = j_start, j_end
! DO k = kts+1,ktf
! DO i = i_start, i_end
! rravg(i,k,j)=fnm(k)*rr(i,k,j)+fnp(k)*rr(i,k-1,j)
! H3(i,k,j)=rravg(i,k,j)*zeta_z(i,j)*H3(i,k,j)
! ENDDO
! ENDDO
! ENDDO
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
! tendency(i,k,j)=tendency(i,k,j)- &
! (H3(i,k+1,j)-H3(i,k,j))/dnw(k)
tendency(i,k,j)=tendency(i,k,j)-mu(i,j)*(H3(i,k+1,j)-H3(i,k,j))*rdzw(i,k,j)
ENDDO
ENDDO
ENDDO
END SUBROUTINE vertical_diffusion_s_qv
!==========================================================================
SUBROUTINE calculate_khh( xkmh, xkhh, 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( its-1:ite+1, kts:kte, jts-1:jte+1), &
INTENT(OUT ) :: xkhh
REAL , DIMENSION( ims:ime ,kms:kme, jms:jme), INTENT(IN ) :: xkmh
!LOCAL VAR
INTEGER :: i, j, k, ktf, i_start, i_end, j_start, j_end
! REAL , DIMENSION( its-1:ite+1 ,kts:kte, jts-1:jte+1) :: pr
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! DO j = j_start-1, j_end+1
! DO k = kts,ktf
! DO i = i_start-1, i_end+1
! pr(i,k,j)=1./3.
! ENDDO
! ENDDO
! ENDDO
! we have hardwired the prandtl number to three, probably should be
! a model constant.
DO j = j_start-1, j_end+1
DO k = kts,ktf
DO i = i_start-1, i_end+1
xkhh(i,k,j)=xkmh(i,k,j)*prandtl
ENDDO
ENDDO
ENDDO
END SUBROUTINE calculate_khh
!==========================================================================
SUBROUTINE calculate_khv( xkmv, xkhv, 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( its:ite, kts:kte, jts:jte), INTENT(OUT ) :: xkhv
REAL , DIMENSION( ims:ime ,kms:kme, jms:jme), INTENT(IN ) :: xkmv
!LOCAL VAR
INTEGER :: i, j, k, ktf, i_start, i_end, j_start, j_end
! REAL , DIMENSION( its-1:ite+1,kts:kte,jts-1:jte+1) :: pr
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
! DO j = j_start, j_end
! DO k = kts,ktf
! DO i = i_start, i_end
! pr(i,k,j)=1./3.
! ENDDO
! ENDDO
! ENDDO
! we have hardwired the prandtl number to three, probably should be
! a model constant.
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
xkhv(i,k,j)=xkmv(i,k,j)*prandtl
ENDDO
ENDDO
ENDDO
END SUBROUTINE calculate_khv
!==========================================================================
SUBROUTINE cal_titau_11_22_33( config_flags,titau,mu,tke,xkx,defor,div, & 6
is_ext,ie_ext,js_ext,je_ext, &
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 ) :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
INTENT(INOUT) :: titau
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: defor, &
xkx, &
div, &
tke
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! LOCAL VAR
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL :: x2r3
!--------------------------------------------------------------------------
x2r3 = 2./3.
ktf=MIN(kte,kde-1)
i_start = its
i_end = ite
j_start = jts
j_end = jte
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-1,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-1,jte)
i_start = i_start-is_ext
i_end = i_end +ie_ext
j_start = j_start-js_ext
j_end = j_end +je_ext
IF ( config_flags%km_opt .eq. 2) THEN
DO j = j_start,j_end
DO k = kts,ktf
DO i = i_start,i_end
titau(i,k,j)=mu(i,j)*(x2r3*max(tke(i,k,j),0.0)-xkx(i,k,j)* &
(defor(i,k,j)-x2r3*div(i,k,j)))
ENDDO
ENDDO
ENDDO
ELSE
DO j = j_start,j_end
DO k = kts,ktf
DO i = i_start,i_end
titau(i,k,j)=-mu(i,j)*xkx(i,k,j)*defor(i,k,j)
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE cal_titau_11_22_33
!==========================================================================
SUBROUTINE cal_titau_12_21( config_flags,titau,mu,xkx,defor, & 3
is_ext,ie_ext,js_ext,je_ext, &
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 ) :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
INTENT(INOUT) :: titau
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: defor, &
xkx
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! LOCAL VAR
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: xkxavg
REAL , DIMENSION( its-1:ite+1, jts-1:jte+1) :: muavg
!--------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
! needs one more point in the x and y directions
i_start = its
i_end = ite
j_start = jts
j_end = jte
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-1,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-1,jte)
i_start = i_start-is_ext
i_end = i_end +ie_ext
j_start = j_start-js_ext
j_end = j_end +je_ext
DO j = j_start,j_end
DO k = kts,ktf
DO i = i_start,i_end
xkxavg(i,k,j)=0.25*(xkx(i-1,k,j )+xkx(i,k,j )+ &
xkx(i-1,k,j-1)+xkx(i,k,j-1))
ENDDO
ENDDO
ENDDO
DO j = j_start,j_end
DO i = i_start,i_end
muavg(i,j) =0.25*(mu(i-1,j )+mu(i,j )+ &
mu(i-1,j-1)+mu(i,j-1))
ENDDO
ENDDO
! titau12 or titau21
DO j = j_start,j_end
DO k = kts,ktf
DO i = i_start,i_end
titau(i,k,j)=-muavg(i,j)*xkxavg(i,k,j)*defor(i,k,j)
ENDDO
ENDDO
ENDDO
END SUBROUTINE cal_titau_12_21
!==========================================================================
SUBROUTINE cal_titau_13_31( config_flags,titau,defor,mu,xkx,fnm,fnp, & 4
is_ext,ie_ext,js_ext,je_ext, &
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 ) :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
INTENT(INOUT) :: titau
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: defor, &
xkx
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! LOCAL VAR
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: xkxavg
REAL , DIMENSION( its-1:ite+1, jts-1:jte+1 ) :: muavg
!--------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
! need ide-1 and jde-1 for averaging to p point
i_start = its
i_end = ite
j_start = jts
j_end = MIN(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-1,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
i_start = i_start-is_ext
i_end = i_end +ie_ext
j_start = j_start-js_ext
j_end = j_end +je_ext
DO j = j_start,j_end
DO k = kts+1,ktf
DO i = i_start,i_end
xkxavg(i,k,j)=0.5*(fnm(k)*(xkx(i,k ,j)+xkx(i-1,k ,j))+ &
fnp(k)*(xkx(i,k-1,j)+xkx(i-1,k-1,j)))
ENDDO
ENDDO
ENDDO
DO j = j_start,j_end
DO i = i_start,i_end
muavg(i,j) =0.5*(mu(i,j)+mu(i-1,j))
ENDDO
ENDDO
DO j = j_start,j_end
DO k = kts+1,ktf
DO i = i_start,i_end
titau(i,k,j)=-muavg(i,j)*xkxavg(i,k,j)*defor(i,k,j)
ENDDO
ENDDO
ENDDO
DO j = j_start,j_end
DO i = i_start,i_end
titau(i,kts,j)=0.
titau(i,ktf+1,j)=0.
ENDDO
ENDDO
END SUBROUTINE cal_titau_13_31
!==========================================================================
SUBROUTINE cal_titau_23_32(config_flags,titau,defor,mu,xkx,fnm,fnp, & 4
is_ext,ie_ext,js_ext,je_ext, &
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 ) :: is_ext,ie_ext,js_ext,je_ext
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1), &
INTENT(INOUT) :: titau
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: defor, &
xkx
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! LOCAL VAR
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: xkxavg
REAL , DIMENSION( its-1:ite+1, jts-1:jte+1) :: muavg
!--------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
! need ide-1 and jde-1 for averaging to p point
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = jte
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-1,jte)
i_start = i_start-is_ext
i_end = i_end +ie_ext
j_start = j_start-js_ext
j_end = j_end +je_ext
DO j = j_start,j_end
DO k = kts+1, ktf
DO i = i_start,i_end
xkxavg(i,k,j)=0.5*(fnm(k)*(xkx(i,k ,j)+xkx(i,k ,j-1))+ &
fnp(k)*(xkx(i,k-1,j)+xkx(i,k-1,j-1)))
ENDDO
ENDDO
ENDDO
DO j = j_start,j_end
DO i = i_start,i_end
muavg(i,j) =0.5*(mu(i,j)+mu(i,j-1))
ENDDO
ENDDO
DO j = j_start,j_end
DO k = kts+1, ktf
DO i = i_start,i_end
titau(i,k,j)=-muavg(i,j)*xkxavg(i,k,j)*defor(i,k,j)
ENDDO
ENDDO
ENDDO
DO j = j_start,j_end
DO i = i_start,i_end
titau(i,kts,j)=0.
titau(i,ktf+1,j)=0.
ENDDO
ENDDO
END SUBROUTINE cal_titau_23_32
!===========================================================================
SUBROUTINE phy_bc ( config_flags,div,defor11,defor22,defor33, & 1,15
defor12,defor13,defor23,xkmh,xkmhd,xkmv,xkhh,xkhv,tke, &
RUBLTEN, RVBLTEN, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
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, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::RUBLTEN, &
RVBLTEN, &
defor11, &
defor22, &
defor33, &
defor12, &
defor13, &
defor23, &
xkmh, &
xkmhd, &
xkmv, &
xkhh, &
xkhv, &
tke, &
div
IF(config_flags%bl_pbl_physics .GT. 0) THEN
CALL set_physical_bc3d( RUBLTEN , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( RVBLTEN , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
ENDIF
IF(config_flags%diff_opt .eq. 2) THEN
CALL set_physical_bc3d( xkmh , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( xkmhd , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( xkmv , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( xkhh , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( xkhv , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( tke , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( div , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( defor11 , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( defor22 , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( defor33 , 't', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( defor12 , 'd', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( defor13 , 'e', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
CALL set_physical_bc3d( defor23 , 'f', config_flags, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
ips, ipe, jps, jpe, kps, kpe, &
its, ite, jts, jte, kts, kte )
ENDIF
END SUBROUTINE phy_bc
!==========================================================================
SUBROUTINE tke_rhs ( tendency,BN2,config_flags, & 1,3
defor11,defor22,defor33, &
defor12,defor13,defor23,u,v,w,div,tke,mu, &
theta,p,p8w,t8w,z,fnm,fnp,cf1,cf2,cf3, &
msft,xkmh,xkmv,xkhv,rdx,rdy,dx,dy,dt, &
zx,zy,rdz,rdzw,dn,dnw,cr_len, &
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 , INTENT(IN ) :: cf1, cf2, cf3, dt
REAL , INTENT(IN ) :: rdx, rdy, dx, dy, cr_len
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msft
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
defor22, &
defor33, &
defor12, &
defor13, &
defor23, &
div, &
BN2, &
tke, &
xkmh, &
xkmv, &
xkhv, &
zx, &
zy, &
u, &
v, &
w, &
theta, &
p, &
p8w, &
t8w, &
z, &
rdz, &
rdzw
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
!--------------------------------------------------------------------------
! call shear stress production term
CALL tke_shear( tendency,config_flags,defor11,defor22,defor33, &
defor12,defor13,defor23,u,v,w,div,tke,mu, &
fnm,fnp,cf1,cf2,cf3,msft,xkmh,xkmv,rdx,rdy, &
zx,zy,rdz, rdzw, dnw, dn, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
CALL tke_dissip( tendency,config_flags,mu,tke,theta,p8w,t8w,z, &
dx, dy,rdz, rdzw, cr_len, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
CALL tke_buoyancy( tendency,config_flags,mu,tke,xkhv,BN2,dt, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
END SUBROUTINE tke_rhs
!==========================================================================
SUBROUTINE tke_buoyancy( tendency,config_flags,mu,tke,xkhv,BN2,dt, & 1
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 , INTENT(IN ) :: dt
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: xkhv, &
tke, &
BN2
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! Local VAR
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
!-------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
!
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tendency(i,k,j)= tendency(i,k,j) - &
mu(i,j)*xkhv(i,k,j)*BN2(i,k,j)
tendency(i,k,j)= max(tendency(i,k,j),-mu(i,j)*tke(i,k,j)/dt)
ENDDO
ENDDO
ENDDO
END SUBROUTINE tke_buoyancy
!==========================================================================
SUBROUTINE tke_dissip ( tendency,config_flags,mu,tke,theta,p8w,t8w,z, & 1
dx, dy, rdz, rdzw ,cr_len_in, &
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 , INTENT(IN ) :: dx, dy, cr_len_in
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) ::tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) :: tke, &
theta, &
p8w, &
t8w, &
z, &
rdz, &
rdzw
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! Local VAR
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: dthrdn
REAL , DIMENSION( its:ite ) :: sumtke, &
sumtkez
INTEGER :: i, j, k, ktf
INTEGER :: i_start, i_end, j_start, j_end
REAL :: disp_len, deltas, coefc, tmpdz, len_s, thetasfc, &
thetatop, dissip_l, len_0, tketmp, cr_len
!--------------------------------------------------------------------------
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. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
cr_len = cr_len_in
cr_len = dx + 1. ! hardwire for mixing length = (dx*dy*dz)**(1/3).
! remove this for the alternate formulation
IF (dx .gt. cr_len) THEN
DO j = j_start, j_end
DO i = i_start, i_end
sumtke(i)=0.
sumtkez(i)=0.
ENDDO
DO k = kts, ktf
DO i = i_start, i_end
tketmp=max(tke(i,k,j),0.)
! sumtke(i)=sumtke(i)+sqrt(tketmp)*dnw(k)
sumtke(i)=sumtke(i)+sqrt(tketmp)/rdzw(i,k,j)
sumtkez(i)=sumtkez(i)+sumtke(i)*z(i,k,j)
IF (abs(sumtke(i)) .gt. 0.01) THEN
len_0=0.2*sumtkez(i)/sumtke(i)
ELSE
len_0=80.
ENDIF
len_0=min(80.,len_0)
dissip_l=KARMAN*z(i,k,j)/(1.+KARMAN*z(i,k,j)/len_0)
tendency(i,k,j)= tendency(i,k,j) - &
mu(i,j)*2.*sqrt(2.)/15.* &
tketmp**1.5/dissip_l
ENDDO
ENDDO
ENDDO
ELSE
DO j = j_start, j_end
DO k = kts+1, ktf-1
DO i = i_start, i_end
! tmpdz=dn(k+1)+dn(k)
tmpdz = 1./rdz(i,k+1,j)+1./rdz(i,k,j)
dthrdn(i,k,j)=(theta(i,k+1,j)-theta(i,k-1,j))/tmpdz
ENDDO
ENDDO
ENDDO
k=kts
DO j = j_start, j_end
DO i = i_start, i_end
! tmpdz=dn(k+1)+0.5*dnw(k)
tmpdz = 1./rdz(i,k+1,j)+0.5/rdzw(i,k,j)
thetasfc=T8w(i,kts,j)/(p8w(i,k,j)/p1000mb)**(R_d/Cp)
dthrdn(i,k,j)=(theta(i,k+1,j)-thetasfc)/tmpdz
ENDDO
ENDDO
k=ktf
DO j = j_start, j_end
DO i = i_start, i_end
! tmpdz=dn(k)+0.5*dnw(k)
tmpdz = 1./rdz(i,k,j)+0.5/rdzw(i,k,j)
thetatop=T8w(i,kde,j)/(p8w(i,kde,j)/p1000mb)**(R_d/Cp)
dthrdn(i,k,j)=(thetatop-theta(i,k-1,j))/tmpdz
ENDDO
ENDDO
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
! deltas=(dx*dy*dnw(k)/zeta_z(i,j))**0.33333333
deltas=(dx*dy/rdzw(i,k,j))**0.33333333
dissip_l=deltas
tketmp=max(tke(i,k,j),0.0)
IF (dthrdn(i,k,j) .gt. 0.) THEN
len_s=max(0.76*sqrt(tketmp)/ &
(abs(g/theta(i,k,j)*dthrdn(i,k,j)))**0.5, 1.0E-5)
! here we are setting coefc (c_e) = .93, uncomment the following line to
! retrieve the alternate formulation
! dissip_l=min(dissip_l,len_s)
ENDIF
coefc=0.19+0.74*dissip_l/deltas
! coefc=.20
tendency(i,k,j)= tendency(i,k,j) - &
mu(i,j)*coefc*tketmp**1.5/dissip_l
ENDDO
ENDDO
ENDDO
ENDIF
END SUBROUTINE tke_dissip
!==========================================================================
SUBROUTINE tke_shear ( tendency,config_flags,defor11,defor22,defor33, & 1,8
defor12,defor13,defor23,u,v,w,div,tke,mu, &
fnm,fnp,cf1,cf2,cf3,msft,xkmh,xkmv,rdx,rdy, &
zx,zy,rdz,rdzw,dn,dnw, &
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 , INTENT(IN ) :: cf1, cf2, cf3
REAL , INTENT(IN ) :: rdx, rdy
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnm
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: fnp
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dn
REAL , DIMENSION( kms:kme ) , INTENT(IN ) :: dnw
REAL , DIMENSION( ims:ime, jms:jme) , INTENT(IN ) :: msft
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(INOUT) :: tendency
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT(IN ) ::defor11, &
defor22, &
defor33, &
defor12, &
defor13, &
defor23, &
div, &
tke, &
xkmh, &
xkmv, &
zx, &
zy, &
u, &
v, &
w, &
rdz, &
rdzw
REAL , DIMENSION( ims:ime, jms:jme), INTENT(IN ) :: mu
! Local VAR
INTEGER :: i, j, k, ktf, ktes1, ktes2
INTEGER :: i_start, i_end, j_start, j_end
INTEGER :: is_ext,ie_ext,js_ext,je_ext
REAL :: mtau
REAL , DIMENSION( its-1:ite+1, kts:kte, jts-1:jte+1) :: avg, &
titau
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: titau12, &
tmp1
! new
REAL , DIMENSION( its:ite, kts:kte, jts:jte) :: zxavg, &
zyavg
!-------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
ktes1=kte-1
ktes2=kte-2
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
IF ( config_flags%open_xs .or. config_flags%specified .or. &
config_flags%nested) i_start = MAX(ids+1,its)
IF ( config_flags%open_xe .or. config_flags%specified .or. &
config_flags%nested) i_end = MIN(ide-2,ite)
IF ( config_flags%open_ys .or. config_flags%specified .or. &
config_flags%nested) j_start = MAX(jds+1,jts)
IF ( config_flags%open_ye .or. config_flags%specified .or. &
config_flags%nested) j_end = MIN(jde-2,jte)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
! new
zxavg(i,k,j)=0.25*( zx(i,k ,j)+zx(i+1,k ,j)+ &
zx(i,k+1,j)+zx(i+1,k+1,j) )
zyavg(i,k,j)=0.25*( zy(i,k ,j)+zy(i,k ,j+1)+ &
zy(i,k+1,j)+zy(i,k+1,j+1) )
ENDDO
ENDDO
ENDDO
!----------------
! avg = u at w point
do j=j_start,j_end
do k=kts+1, ktf
do i=i_start,i_end
avg(i,k,j)=0.5*(fnm(k)*(u(i,k ,j)+u(i+1,k ,j))+ &
fnp(k)*(u(i,k-1,j)+u(i+1,k-1,j)))
enddo
enddo
enddo
do j=j_start,j_end
do i=i_start,i_end
avg(i,1,j) =0.5*(cf1*u(i ,1,j)+cf2*u(i ,2,j)+cf3*u(i ,3,j)+ &
cf1*u(i+1,1,j)+cf2*u(i+1,2,j)+cf3*u(i+1,3,j))
avg(i,kte,j)=0.5*(u(i,ktes1,j)+(u(i,ktes1,j)-u(i,ktes2,j)) &
*0.5*dnw(ktes1)/dn(ktes1)+ &
u(i+1,ktes1,j)+(u(i+1,ktes1,j)-u(i+1,ktes2,j)) &
*0.5*dnw(ktes1)/dn(ktes1))
enddo
enddo
! tmp1 = delta u at p point
do j=j_start,j_end
do k=kts, ktf
do i=i_start,i_end
tmp1(i,k,j)= (avg(i,k+1,j)-avg(i,k,j))
enddo
enddo
enddo
!----------------
! titau = titau11
is_ext=0
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_11_22_33(config_flags,titau,mu,tke, &
xkmh,defor11,div, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! update tendency with tau11 ( partial u / partial x)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
mtau=msft(i,j)*titau(i,k,j)
! new
tendency(i,k,j)= tendency(i,k,j) &
-mtau*( rdx*(u(i+1,k,j)-u(i,k,j))- &
tmp1(i,k,j)*zxavg(i,k,j)*rdzw(i,k,j) )
! tendency(i,k,j)= tendency(i,k,j) &
! -mtau*( rdx*(u(i+1,k,j)-u(i,k,j))- &
! tmp1(i,k,j)*zxavg(i,k,j)*zeta_z(i,j) &
! /dnw(k) )
ENDDO
ENDDO
ENDDO
!----------------
! titau = titau12
is_ext=0
ie_ext=1
js_ext=0
je_ext=1
CALL cal_titau_12_21( config_flags,titau,mu,xkmh,defor12, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau12 = titau12 at p point (keep for later use)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
titau12(i,k,j)=0.25*(titau(i ,k,j)+titau(i ,k,j+1)+ &
titau(i+1,k,j)+titau(i+1,k,j+1))
ENDDO
ENDDO
ENDDO
! avg = u at v point
DO j = j_start, j_end+1
DO k = kts, ktf
DO i = i_start, i_end
avg(i,k,j)=0.25*(u(i,k,j )+u(i+1,k,j )+ &
u(i,k,j-1)+u(i+1,k,j-1))
ENDDO
ENDDO
ENDDO
! update tendency with tau12 ( partial u / partial y)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
mtau=msft(i,j)*titau12(i,k,j)
! new
tendency(i,k,j)= tendency(i,k,j) &
-mtau*( rdy*(avg(i,k,j+1)-avg(i,k,j)) - &
tmp1(i,k,j)*zyavg(i,k,j)*rdzw(i,k,j) )
! tendency(i,k,j)= tendency(i,k,j) &
! -mtau*( rdy*(avg(i,k,j+1)-avg(i,k,j)) - &
! tmp1(i,k,j)*zyavg(i,k,j)*zeta_z(i,j) &
! /dnw(k) )
ENDDO
ENDDO
ENDDO
!---------------
! avg = titau13
is_ext=0
ie_ext=1
js_ext=0
je_ext=0
CALL cal_titau_13_31(config_flags,avg,defor13,mu,xkmv,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau = titau 13 at p point
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
titau(i,k,j)=0.25*(avg(i ,k+1,j)+avg(i ,k,j)+ &
avg(i+1,k+1,j)+avg(i+1,k,j))
ENDDO
ENDDO
ENDDO
! update tendency with tau13 ( partial u / partial z)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tendency(i,k,j)= tendency(i,k,j) - &
titau(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j)
! tendency(i,k,j)= tendency(i,k,j) - &
! titau(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k)
ENDDO
ENDDO
ENDDO
!-------------------------------------
! avg = v at w point
do j=j_start,j_end
do k=kts+1,ktf
do i=i_start,i_end
avg(i,k,j)=0.5*(fnm(k)*(v(i,k ,j)+v(i,k ,j+1))+ &
fnp(k)*(v(i,k-1,j)+v(i,k-1,j+1)))
enddo
enddo
enddo
do j=j_start,j_end
do i=i_start,i_end
avg(i,1,j) =0.5*(cf1*v(i,1,j )+cf2*v(i,2,j )+cf3*v(i,3,j )+ &
cf1*v(i,1,j+1)+cf2*v(i,2,j+1)+cf3*v(i,3,j+1))
avg(i,kte,j)=0.5*(v(i,ktes1,j)+(v(i,ktes1,j)-v(i,ktes2,j)) &
*0.5*dnw(ktes1)/dn(ktes1)+ &
v(i,ktes1,j+1)+(v(i,ktes1,j+1)-v(i,ktes2,j+1)) &
*0.5*dnw(ktes1)/dn(ktes1))
enddo
enddo
! tmp1 = delta v at p point
do j=j_start,j_end
do k=kts, ktf
do i=i_start,i_end
tmp1(i,k,j)= (avg(i,k+1,j)-avg(i,k,j))
enddo
enddo
enddo
!-------------------
! avg = v at u point
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end+1
avg(i,k,j)=0.25*(v(i ,k,j)+v(i ,k,j+1)+ &
v(i-1,k,j)+v(i-1,k,j+1))
ENDDO
ENDDO
ENDDO
! update tendency with tau21 ( partial v / partial x)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
mtau=msft(i,j)*titau12(i,k,j)
! new
tendency(i,k,j)= tendency(i,k,j) &
-mtau*(rdx*(avg(i+1,k,j)-avg(i,k,j)) - &
zxavg(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j) )
! tendency(i,k,j)= tendency(i,k,j) &
! -mtau*(rdx*(avg(i+1,k,j)-avg(i,k,j)) - &
! zxavg(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k))
ENDDO
ENDDO
ENDDO
!-------------------------------------
! titau = titau22
is_ext=0
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_11_22_33(config_flags,titau,mu,tke, &
xkmh,defor22,div, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! update tendency with tau 22 ( partial v / partial y)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
mtau=msft(i,j)*titau(i,k,j)
! new
tendency(i,k,j)= tendency(i,k,j) &
-mtau*(rdy*(v(i,k,j+1)-v(i,k,j)) - &
zyavg(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j) )
! tendency(i,k,j)= tendency(i,k,j) &
! -mtau*(rdy*(v(i,k,j+1)-v(i,k,j)) - &
! zyavg(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k))
ENDDO
ENDDO
ENDDO
!-------------------------------------
! avg = titau23
is_ext=0
ie_ext=0
js_ext=0
je_ext=1
CALL cal_titau_23_32( config_flags,avg,defor23,mu,xkmv,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau = titau 23 at p point
DO j = j_start, j_end
DO k = kts+1,ktf
DO i = i_start, i_end
titau(i,k,j)=0.25*(avg(i,k+1,j )+avg(i,k,j )+ &
avg(i,k+1,j+1)+avg(i,k,j+1))
ENDDO
ENDDO
ENDDO
! approximation for grid with constant dz
DO j = j_start, j_end
k = kts
DO i = i_start, i_end
titau(i,k,j)= 0.5*( 1.5*(avg(i,k+1,j)+avg(i,k+1,j+1)) &
-0.5*(avg(i,k+2,j)+avg(i,k+2,j+1)) )
ENDDO
ENDDO
! update tendency with tau23 ( partial v / partial z)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tendency(i,k,j)= tendency(i,k,j) - &
titau(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j)
! tendency(i,k,j)= tendency(i,k,j) - &
! titau(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k)
ENDDO
ENDDO
ENDDO
!-------------------------------------
! tmp1 = partial w over partial zeta at p point
do j=j_start,j_end
do k=kts,ktf
do i=i_start,i_end
tmp1(i,k,j)= (w(i,k+1,j)-w(i,k,j))
enddo
enddo
enddo
!-------------------------------------
! avg = titau31
is_ext=0
ie_ext=1
js_ext=0
je_ext=0
CALL cal_titau_13_31(config_flags,avg,defor13,mu,xkmh,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau = titau 31 at p point
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
titau(i,k,j)=0.25*(avg(i ,k+1,j)+avg(i ,k,j)+ &
avg(i+1,k+1,j)+avg(i+1,k,j))
ENDDO
ENDDO
ENDDO
! avg = w at u point
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end+1
avg(i,k,j)=0.25*(w(i ,k,j)+w(i ,k+1,j)+ &
w(i-1,k,j)+w(i-1,k+1,j))
ENDDO
ENDDO
ENDDO
! update tendency with tau31 ( partial w / partial x)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
mtau=msft(i,j)*titau(i,k,j)
! new
tendency(i,k,j)= tendency(i,k,j) &
-mtau*(rdx*(avg(i+1,k,j)-avg(i,k,j)) - &
zxavg(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j) )
! tendency(i,k,j)= tendency(i,k,j) &
! -mtau*(rdx*(avg(i+1,k,j)-avg(i,k,j)) - &
! zxavg(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k))
ENDDO
ENDDO
ENDDO
!-------------------------------------
! avg = titau32
is_ext=0
ie_ext=0
js_ext=0
je_ext=1
CALL cal_titau_23_32( config_flags,avg,defor23,mu,xkmh,fnm,fnp, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! titau = titau 32 at p point
DO j = j_start, j_end
DO k = kts,ktf
DO i = i_start, i_end
titau(i,k,j)=0.25*(avg(i,k+1,j )+avg(i,k,j )+ &
avg(i,k+1,j+1)+avg(i,k,j+1))
ENDDO
ENDDO
ENDDO
! avg = w at v point
DO j = j_start, j_end+1
DO k = kts, ktf
DO i = i_start, i_end
avg(i,k,j)=0.25*(w(i,k,j )+w(i,k+1,j )+ &
w(i,k,j-1)+w(i,k+1,j-1))
ENDDO
ENDDO
ENDDO
! update tendency with tau32 ( partial w / partial y)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
mtau=msft(i,j)*titau(i,k,j)
! new
tendency(i,k,j)= tendency(i,k,j) &
-mtau*(rdy*(avg(i,k,j+1)-avg(i,k,j)) - &
zyavg(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j) )
! tendency(i,k,j)= tendency(i,k,j) &
! -mtau*(rdy*(avg(i,k,j+1)-avg(i,k,j)) - &
! zyavg(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k))
ENDDO
ENDDO
ENDDO
!-------------------------------------
! titau = titau33
is_ext=0
ie_ext=0
js_ext=0
je_ext=0
CALL cal_titau_11_22_33(config_flags,titau,mu,tke, &
xkmv,defor33,div, &
is_ext,ie_ext,js_ext,je_ext, &
ids, ide, jds, jde, kds, kde, &
ims, ime, jms, jme, kms, kme, &
its, ite, jts, jte, kts, kte )
! update tendency with tau33 ( partial w / partial z)
DO j = j_start, j_end
DO k = kts, ktf
DO i = i_start, i_end
tendency(i,k,j)= tendency(i,k,j) - &
titau(i,k,j)*tmp1(i,k,j)*rdzw(i,k,j)
! tendency(i,k,j)= tendency(i,k,j) - &
! titau(i,k,j)*zeta_z(i,j)*tmp1(i,k,j)/dnw(k)
ENDDO
ENDDO
ENDDO
END SUBROUTINE tke_shear
!==========================================================================
SUBROUTINE compute_diff_metrics ( config_flags, ph, phb, z, rdz, rdzw, & 1
zx, zy, rdx, rdy, &
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 ) :: ph, &
phb
REAL , DIMENSION( ims:ime, kms:kme, jms:jme), INTENT( OUT) :: rdz, &
rdzw, &
zx, &
zy, &
z
REAL , INTENT(IN ) :: rdx, rdy
! Local VAR
INTEGER :: i, j, k
INTEGER :: i_start, i_end, j_start, j_end, ktf
!-------------------------------------------------------------------------
ktf=MIN(kte,kde-1)
!!! bug fix, wcs, 22 april 2002
!!! we need rdzw in halo for average to u and v points
! i_start = its
! i_end = MIN(ite,ide-1)
! j_start = jts
! j_end = MIN(jte,jde-1)
j_start = jts-1
j_end = jte
! begin with dz computations
DO j = j_start, j_end
if( (j_start >= jts) .and. (j_end <= MIN(jte,jde-1)) ) then
i_start = its-1
i_end = ite
else
i_start = its
i_end = MIN(ite,ide-1)
end if
! compute z at w points for rdz and rdzw computations,
! we'll switch z to z at p points before returning
DO k = 1, kte
!!! bug fix, wcs, 22 april 2002
!!! DO i = 1, i_end
DO i = i_start, i_end
z(i,k,j) = (ph(i,k,j)+phb(i,k,j))/g
ENDDO
ENDDO
DO k = 1, ktf
DO i = i_start, i_end
rdzw(i,k,j) = 1./(z(i,k+1,j)-z(i,k,j))
ENDDO
ENDDO
DO k = 2, ktf
DO i = i_start, i_end
rdz(i,k,j) = 2./(z(i,k+1,j)-z(i,k-1,j))
ENDDO
ENDDO
!!! bug fix, wcs, 22 april 2002
!!! added the following code
DO i = i_start, i_end
rdz(i,1,j) = 2./(z(i,2,j)-z(i,1,j))
ENDDO
ENDDO
!!! end bug fix
! now compute zx and zy - we'll assume that
! the halo for ph and phb is properly filled.
i_start = its
i_end = MIN(ite,ide-1)
j_start = jts
j_end = MIN(jte,jde-1)
DO j = j_start, j_end
DO k = 1, kte
DO i = max(ids+1,its), i_end
zx(i,k,j) = rdx*(phb(i,k,j)-phb(i-1,k,j))/g
ENDDO
ENDDO
ENDDO
DO j = j_start, j_end
DO k = 1, kte
DO i = max(ids+1,its), i_end
zx(i,k,j) = zx(i,k,j)+ rdx*(ph(i,k,j)-ph(i-1,k,j))/g
ENDDO
ENDDO
ENDDO
DO j = max(jds+1,jts), j_end
DO k = 1, kte
DO i = i_start, i_end
zy(i,k,j) = rdy*(phb(i,k,j)-phb(i,k,j-1))/g
ENDDO
ENDDO
ENDDO
DO j = max(jds+1,jts), j_end
DO k = 1, kte
DO i = i_start, i_end
zy(i,k,j) = zy(i,k,j)+ rdy*(ph(i,k,j)-ph(i,k,j-1))/g
ENDDO
ENDDO
ENDDO
! DO j = j_start, j_end
! DO k = 1, kte
! DO i = max(ids+1,its), i_end
! zx(i,k,j) = rdx*(ph(i,k,j)+phb(i,k,j)-ph(i-1,k,j)-phb(i-1,k,j))/g
! ENDDO
! ENDDO
! ENDDO
! DO j = max(jds+1,jts), j_end
! DO k = 1, kte
! DO i = i_start, i_end
! zy(i,k,j) = rdy*(ph(i,k,j)+phb(i,k,j)-ph(i,k,j-1)-phb(i,k,j-1))/g
! ENDDO
! ENDDO
! ENDDO
! some b.c on zx and zy
IF (.not. config_flags%periodic_x) THEN
IF ( ite == ide ) THEN
DO j=j_start,j_end
DO k=1,ktf
zx(ide,k,j) = 0.
ENDDO
ENDDO
ENDIF
IF ( its == ids ) THEN
DO j=j_start,j_end
DO k=1,ktf
zx(ids,k,j) = 0.
ENDDO
ENDDO
ENDIF
ELSE
IF ( ite == ide ) THEN
DO j=j_start,j_end
DO k=1,ktf
zx(ide,k,j) = rdx*(phb(ide,k,j)-phb(ide-1,k,j))/g
ENDDO
ENDDO
DO j=j_start,j_end
DO k=1,ktf
zx(ide,k,j) = zx(ide,k,j) + rdx*(ph(ide,k,j)-ph(ide-1,k,j))/g
ENDDO
ENDDO
ENDIF
IF ( its == ids ) THEN
DO j=j_start,j_end
DO k=1,ktf
zx(ids,k,j) = rdx*(phb(ids,k,j)-phb(ids-1,k,j))/g
ENDDO
ENDDO
DO j=j_start,j_end
DO k=1,ktf
zx(ids,k,j) = zx(ids,k,j) + rdx*(ph(ids,k,j)-ph(ids-1,k,j))/g
ENDDO
ENDDO
ENDIF
END IF
IF (.not. config_flags%periodic_y) THEN
IF ( jte == jde ) THEN
DO k=1,ktf
DO i=i_start,i_end
zy(i,k,jde) = 0.
ENDDO
ENDDO
ENDIF
IF ( jts == jds ) THEN
DO k=1,ktf
DO i=i_start,i_end
zy(i,k,jds) = 0.
ENDDO
ENDDO
ENDIF
ELSE
IF ( jte == jde ) THEN
DO j=j_start,j_end
DO k=1,ktf
zy(i,k,jde) = rdy*(phb(i,k,jde)-phb(i,k,jde-1))/g
ENDDO
ENDDO
DO j=j_start,j_end
DO k=1,ktf
zy(i,k,jde) = zy(i,k,jde) + rdy*(ph(i,k,jde)-ph(i,k,jde-1))/g
ENDDO
ENDDO
ENDIF
IF ( jts == jds ) THEN
DO j=j_start,j_end
DO k=1,ktf
zy(i,k,jds) = rdy*(phb(i,k,jds)-phb(i,k,jds-1))/g
ENDDO
ENDDO
DO j=j_start,j_end
DO k=1,ktf
zy(i,k,jds) = zy(i,k,jds) + rdy*(ph(i,k,jds)-ph(i,k,jds-1))/g
ENDDO
ENDDO
ENDIF
END IF
! output z at p points
DO j = j_start, j_end
DO k = 1, ktf
DO i = i_start, i_end
z(i,k,j) = 0.5*(ph(i,k,j)+phb(i,k,j)+ph(i,k+1,j)+phb(i,k+1,j))/g
ENDDO
ENDDO
ENDDO
END SUBROUTINE compute_diff_metrics
!==========================================================================
END MODULE module_diffusion_em