blsimp.F90

Go to the documentation of this file.

MODULE blsimp

!This subroutine takes the model prognostic variables as its inputs and produces
!the three components of the tri-diagonal matrix used to compute BL diffusion. 
!The same diffussivity Kh is applied to all model variables above the surface. 
!Diffusion coefficient is computed using a simplified Louis type approximation.

IMPLICIT NONE

PRIVATE
PUBLIC bl_simp

CONTAINS

subroutine bl_simp( IM, JM, LM, DT, UT, VT, PTT, QVT, QLT, QIT, PET, PKT, FROCEAN, &
                    MAPL_GRAV, MAPL_VIREPS, MAPL_KAPPA, MAPL_CP, MAPL_RGAS, &
                    AKQ, AKS, AKV, BKQ, BKS, BKV, CKQ, CKS, CKV &
                  )

IMPLICIT NONE


!INPUTS
INTEGER, INTENT(IN)                        :: im, jm, lm
REAL,    INTENT(IN)                        :: dt
REAL,    INTENT(IN), DIMENSION(IM,JM,LM)   :: ut, vt, ptt, qvt, qlt, qit
REAL,    INTENT(IN), DIMENSION(IM,JM,0:LM) :: pet
REAL,    INTENT(IN), DIMENSION(IM,JM,LM)   :: pkt
REAL,    INTENT(IN), DIMENSION(IM,JM)      :: frocean

REAL,    INTENT(IN)                        :: mapl_grav, mapl_vireps, mapl_kappa, mapl_cp, mapl_rgas

!OUTPUTS
REAL, INTENT(OUT), DIMENSION(IM,JM,LM)     :: akq, aks, akv, bkq, bks, bkv, ckq, cks, ckv

!LOCALS
INTEGER :: i, j, l
REAL :: ws, ri, rin
REAL :: ell2, cdrag_land, cdrag_sea, cdrag, tcoef
REAL :: dmi, pkh, dz, ckx, tvb, tve, tvt

REAL, DIMENSION(IM,JM) :: bksq, bksv, bkst
REAL, DIMENSION(LM)    :: kh

 !Initialize outputs
 aks = 0.0
 bks = 0.0
 cks = 0.0
 akq = 0.0
 bkq = 0.0
 ckq = 0.0
 akv = 0.0
 bkv = 0.0
 ckv = 0.0

 !These two are not used
 aks(:,:,1)  = 0.0
 cks(:,:,lm) = 0.0

 ell2       = 30.*30.
 cdrag_land = 0.002
 cdrag_sea  = 0.0015

 ! the following lines solves  ------------------------------------
 ! dQ/dT=-G/dP*Rho*CDrag*WS*Q, Q_(t+dt)=Q_t+dt/(Rho_t*dZ_t)*dTAU
 !                             dTAU=TAU_(z) - TAU_(z-1),   
 !                             TAU_(z)= Rho_t*CDrag*WS*Q_(t+dt), 
 !                             TAU_(z-1)=Rho_t*K_(z-1)/dZ*dQ_(t+dt)
 !                             dZ=RT*dP/(P*G), DMI=(G*dt)/dP
 ! ----------------------------------------------------------------
 
 do j=1,jm
    do i=1,im

       dmi    = (mapl_grav*dt)/(pet(i,j,1)-pet(i,j,0))
       tvt    = ptt(i,j,1)*pkt(i,j,1) &
                 * (1.0 + mapl_vireps *qvt(i,j,1) - qlt(i,j,1) - qit(i,j,1) )
       do l=2,lm
          pkh      = pet(i,j,l)**mapl_kappa
          dz       = mapl_cp*(ptt(i,j,l-1)*(pkh-pkt(i,j,l-1)) + ptt(i,j,l  )*(pkt(i,j,l )-pkh))
          ws       = (ut(i,j,l-1)-ut(i,j,l))**2+(vt(i,j,l-1)-vt(i,j,l))**2 + 0.01
          ri       = mapl_grav*((ptt(i,j,l-1)-ptt(i,j,l)) / (0.5*(ptt(i,j,l-1)+ptt(i,j,l))))*dz/ws
          rin      = ell2*sqrt(ws)/dz
          IF (ri < 0.) THEN
             kh(l) = max(0.01, rin*sqrt(1.-18.*ri))
          ELSE
             kh(l) = max(0.01, rin/(1.+10.*ri*(1.+8.*ri)))
          ENDIF
     
          tvb      = ptt(i,j,l)*pkt(i,j,l)*(1.0 + mapl_vireps *qvt(i,j,l) - qlt(i,j,l) - qit(i,j,l))
          tve      = 0.5*(tvt+tvb)

          tvt      = tvb 
          ckx      = -kh(l)*pet(i,j,l)/( mapl_rgas * tve )/dz
          cks(i,j,l-1) = ckx * dmi
          dmi      = (mapl_grav*dt)/(pet(i,j,l)-pet(i,j,l-1))
          aks(i,j,l)   = ckx * dmi
          bks(i,j,l-1) = 1.0 - (aks(i,j,l-1)+cks(i,j,l-1))

          if (l==lm) then
             bks(i,j,l) = 1.0 - (aks(i,j,l)+cks(i,j,l))
             bksq(i,j)   = bks(i,j,l)
             ws     = sqrt(ut(i,j,l)**2 + vt(i,j,l)**2 + 1.0)
             IF (frocean(i,j).eq.1.0) then
                cdrag=cdrag_sea
                tcoef=1.! CDRAG       ! assume sea surface T unperturbed
             else
                cdrag=cdrag_land
                tcoef=0.          ! assume land surface T same as air T
             endif
             kh(l)  = -cdrag*dmi*ws*pet(i,j,l)/(mapl_rgas*tvb)
             bksv(i,j)   = 1.0 - (aks(i,j,l)+cks(i,j,l)+kh(l))
             bkst(i,j)   = 1.0 - (aks(i,j,l)+cks(i,j,l)+kh(l)*tcoef)
          endif

       end do

    end do
 end do

 !Make copy for wind and tracers as not doing differently for this
 !simplified boundary layer.
 akv = aks
 bkv = bks
 ckv = cks
 akq = aks
 bkq = bks
 ckq = cks

 !Apply surface layer diffusivities
 bkv(:,:,lm) = bksv(:,:)
 bks(:,:,lm) = bkst(:,:)
 bkq(:,:,lm) = bksq(:,:)

endsubroutine bl_simp


ENDMODULE blsimp