Skip to content

Commit 092e520

Browse files
grantfirlgrantfirl
authored
Merge pull request #1154 from grantfirl/ufs-dev-PR253
UFS-dev PR#253
2 parents c9b47a4 + 8a79b16 commit 092e520

File tree

3 files changed

+330
-7
lines changed

3 files changed

+330
-7
lines changed

‎physics/PBL/SATMEDMF/satmedmfvdifq.F

Lines changed: 234 additions & 6 deletions
Original file line numberDiff line numberDiff line change
@@ -7,6 +7,9 @@ module satmedmfvdifq
77
use mfpbltq_mod
88
use tridi_mod
99
use mfscuq_mod
10+
!PCC_CANOPY_utilities
11+
use canopy_utils_mod
12+
1013
contains
1114

1215
!> \defgroup module_satmedmfvdifq GFS TKE-EDMF PBL Module
@@ -74,7 +77,7 @@ end subroutine satmedmfvdifq_init
7477
!! (mfscuq.f).
7578
!! \section detail_satmedmfvidfq GFS satmedmfvdifq Detailed Algorithm
7679
subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
77-
& ntiw,ntke,grav,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
80+
& ntiw,ntke,grav,pi,rd,cp,rv,hvap,hfus,fv,eps,epsm1, &
7881
!The following three variables are for SA-3D-TKE
7982
& def_1,def_2,def_3,sa3dtke,dku3d_h,dku3d_e, &
8083
& dv,du,tdt,rtg,u1,v1,t1,q1,usfco,vsfco,icplocn2atm, &
@@ -85,9 +88,15 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
8588
& dspheat,dusfc,dvsfc,dtsfc,dqsfc,hpbl,dkt,dku,tkeh, &
8689
& kinver,xkzm_m,xkzm_h,xkzm_s,dspfac,bl_upfr,bl_dnfr, &
8790
& rlmx,elmx,sfc_rlm,tc_pbl,use_lpt, &
91+
!IVAI: canopy inputs from AQM
92+
& do_canopy, cplaqm, claie, cfch, cfrt, cclu, cpopu, &
93+
!IVAI
8894
& ntqv,dtend,dtidx,index_of_temperature,index_of_x_wind, &
8995
& index_of_y_wind,index_of_process_pbl,gen_tend,ldiag3d, &
90-
& errmsg,errflg)
96+
& errmsg,errflg)
97+
!IVAI: aux arrays
98+
! & naux2d,naux3d,aux2d,aux3d)
99+
91100
!
92101
use machine , only : kind_phys
93102
use funcphys , only : fpvs
@@ -104,11 +113,17 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
104113
integer, intent(out) :: kpbl(:)
105114
logical, intent(in) :: gen_tend,ldiag3d
106115
!
107-
real(kind=kind_phys), intent(in) :: grav,rd,cp,rv,hvap,hfus,fv, &
116+
real(kind=kind_phys), intent(in) :: grav,pi,rd,cp,rv,hvap,hfus,fv, &
108117
& eps,epsm1
109118
real(kind=kind_phys), intent(in) :: delt, xkzm_m, xkzm_h, xkzm_s
110119
real(kind=kind_phys), intent(in) :: dspfac, bl_upfr, bl_dnfr
111120
real(kind=kind_phys), intent(in) :: rlmx, elmx
121+
!PCC CANOPY------------------------------------
122+
logical, intent(in) :: do_canopy, cplaqm
123+
!IVAI: canopy inputs
124+
real(kind=kind_phys), optional, intent(in) :: claie(:), cfch(:), &
125+
& cfrt(:), cclu(:), cpopu(:)
126+
!----------------------------------------------
112127
real(kind=kind_phys), intent(inout) :: dv(:,:), du(:,:), &
113128
& tdt(:,:), rtg(:,:,:), tkeh(:,:)
114129
real(kind=kind_phys), intent(in) :: &
@@ -289,8 +304,40 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
289304
real(kind=kind_phys) h1
290305

291306
real(kind=kind_phys) bfac, mffac
292-
307+
293308
real(kind=kind_phys) qice(im,km),qliq(im,km)
309+
310+
!PCC_CANOPY------------------------------------
311+
integer COUNTCAN,KCAN
312+
integer kount !IVAI
313+
real(kind=kind_phys) FCH, MOL, HOL, TLCAN,
314+
& SIGMACAN, RRCAN, BBCAN,
315+
& AACAN, ZCAN, ZFL, BOTCAN,
316+
& EDDYVEST1, EDDYVEST_INT
317+
318+
! in canopy eddy diffusivity [ m**2/s ]
319+
real(kind=kind_phys), allocatable :: EDDYVESTX ( : )
320+
! in canopy layer [m]
321+
real(kind=kind_phys), allocatable :: ZCANX ( : )
322+
! Declare local maximum canopy layers
323+
integer, parameter :: MAXCAN = 1000
324+
integer, parameter :: mvt = 30 ! use 30 instead of 27
325+
!Based on MODIS IGBP 20 Category Dataset
326+
real :: fch_table(mvt) !< top of canopy (m)
327+
data ( fch_table(i),i=1,mvt) /
328+
& 20.0, 20.0, 18.0, 16.0, 16.0, 1.10,
329+
& 1.10, 13.0, 10.0, 1.00, 5.00, 2.00,
330+
& 15.0, 1.50, 0.00, 0.00, 0.00, 4.00,
331+
& 2.00, 0.50, 0.00, 0.00, 0.00, 0.00,
332+
& 0.00, 0.00, 0.00, 0.00, 0.00, 0.00 /
333+
!----------------------------------------------
334+
335+
!IVAI
336+
! integer, intent(in) :: naux2d,naux3d
337+
! real(kind_phys), intent(inout) :: aux2d(:,:)
338+
! real(kind_phys), intent(inout) :: aux3d(:,:,:)
339+
!IVAI
340+
294341
!!
295342
parameter(bfac=100.)
296343
parameter(wfac=7.0,cfac=4.5)
@@ -314,7 +361,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
314361
parameter(ck1=0.15,ch1=0.15)
315362
parameter(cs0=0.4,csmf=0.5)
316363
parameter(rchck=1.5,ndt=20)
317-
!The following variables are for SA-3D-TKE
364+
!The following variables are for SA-3D-TKE
318365
parameter(cpl1=0.280,cpl2=0.870,cpl3=0.913)
319366
parameter(cpl4=0.153,cpl5=0.278,cpl6=0.720)
320367
parameter(cpnl1=0.243,cpnl2=0.936,cpnl3=1.11)
@@ -325,6 +372,14 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
325372
parameter(inv3=0.33333333)
326373
! parameter(aa1=0.92,aa2=0.649,bb1=17.7,bb2=9.5,cc1=0.08)
327374

375+
!PCC_CANOPY------------------------------------
376+
if (do_canopy) then
377+
if(.not.allocated(EDDYVESTX))
378+
& allocate( EDDYVESTX ( MAXCAN ) )
379+
if(.not.allocated(ZCANX))
380+
& allocate( ZCANX ( MAXCAN ) )
381+
endif
382+
!----------------------------------------------
328383
if (tc_pbl == 0) then
329384
ck0 = 0.4
330385
ch0 = 0.4
@@ -1562,7 +1617,6 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
15621617
endif
15631618
enddo
15641619
enddo
1565-
15661620
do k = kmscu, 1, -1
15671621
do i = 1, im
15681622
if(scuflg(i) .and.
@@ -1574,6 +1628,179 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
15741628
15751629
endif !sa3dtke
15761630
1631+
!PCC_CANOPY------------------------------------
1632+
kount=0 !IVAI
1633+
if (do_canopy .and. cplaqm) then
1634+
1635+
!IVAI
1636+
! print*, 'SATMEDMFVDIFQ_RUN: CLAIE = ', claie(:)
1637+
! print*, 'SATMEDMFVDIFQ_RUN: CFCH = ' , cfch (:)
1638+
! print*, 'SATMEDMFVDIFQ_RUN: CFRT = ' , cfrt (:)
1639+
! print*, 'SATMEDMFVDIFQ_RUN: CCLU = ' , cclu (:)
1640+
! print*, 'SATMEDMFVDIFQ_RUN: CPOPU= ' , cpopu(:)
1641+
! 2D aux arrays: canopy data in diffusion
1642+
! aux2d(:,1) = cfch (:)
1643+
! aux2d(:,2) = claie(:)
1644+
! aux2d(:,3) = cfrt(:)
1645+
1646+
! 3D aux arrays: before canopy correction
1647+
! aux3d(:,:,1) = dkq(:,:)
1648+
! aux3d(:,:,2) = dkt(:,:)
1649+
! aux3d(:,:,3) = dku(:,:)
1650+
!IVAI
1651+
do k = 1, km1-1
1652+
do i = 1, im
1653+
1654+
!IVAI: AQM canopy Inputs
1655+
! FCH = fch_table(vegtype(i)) !top of canopy from look-up table
1656+
FCH = cfch(i) !top of canopy from AQM canopy inputs
1657+
IF (k .EQ. 1) THEN !use model layer interfaces
1658+
KCAN = 1
1659+
ELSE
1660+
IF ( cfch(i) .GT. zi(i,k)
1661+
& .AND. cfch(i) .LE. zi(i,k+1) ) THEN
1662+
KCAN = 1
1663+
ELSE
1664+
KCAN = 0
1665+
END IF
1666+
END IF
1667+
1668+
IF (KCAN .EQ. 1) THEN !canopy inside model layer
1669+
! Check for other Contiguous Canopy Grid Cell Conditions
1670+
1671+
! Not a contigous canopy cell
1672+
IF ( claie(i) .LT. 0.1
1673+
& .OR. cfch (i) .LT. 0.5
1674+
!IVAI: modified contiguous canopy condition
1675+
! & .OR. MAX(0.0, 1.0 - cfrt(i)) .GT. 0.5
1676+
& .OR. MAX(0.0, 1.0 - cfrt(i)) .GT. 0.75
1677+
!IVAI
1678+
& .OR. cpopu(i) .GT. 10000.0
1679+
& .OR. (EXP(-0.5*claie(i)*cclu(i)) .GT. 0.45
1680+
& .AND. cfch(i) .LT. 18.) ) THEN
1681+
1682+
1683+
!TODO: Canopy Inputs
1684+
! IF ( XCANOPYLAI .LT. 0.1 !from canopy inputs
1685+
! IF ( lai(i) .LT. 0.1 !from LSM
1686+
! & .OR. FCH .LT. 0.5 ) THEN
1687+
! & .OR. MAX(0.0, 1.0 - XCANOPYFRT) .GT. 0.5
1688+
! & .OR. POPU .GT. 10000.0
1689+
! & .OR. EXP(-0.5*XCANOPYLAI*XCANOPYCLU).GT. 0.45
1690+
! & .AND. FCH .LT. 18.0 ) THEN
1691+
1692+
dkt(i,k)= dkt(i,k)
1693+
dkq(i,k)= dkq(i,k)
1694+
dku(i,k)= dku(i,k)
1695+
1696+
ELSE ! There is a contiguous forest canopy, apply correction over canopy layers
1697+
1698+
! Output contiguous canopy mask
1699+
! if (kount .EQ. 0 ) aux2d(i,5) = aux2d(i,5) + 1
1700+
1701+
!Raupauch M. R. A Practical Lagrangian method for relating scalar
1702+
!concentrations to
1703+
! source distributions in vegetation canopies. Q. J. R. Meteor. Soc.
1704+
! (1989), 115, pp 609-632
1705+
MOL = zol(i)/zl(i,k) !Monin-Obukhov Length in layer
1706+
HOL = FCH/MOL !local canopy stability parameter (hc/MOL)
1707+
ZCAN = zi(i,k+1) ! Initialize each model layer top that contains canopy (m)
1708+
! Integrate across total model interface
1709+
ZFL = ZCAN ! Set ZFL = ZCAN
1710+
COUNTCAN = 0 ! Initialize canopy layers
1711+
1712+
IF (k .EQ. 1) THEN !Find bottom in each model layer
1713+
BOTCAN = 0.5
1714+
ELSE
1715+
BOTCAN = zi(i,k)
1716+
END IF
1717+
1718+
DO WHILE (ZCAN.GE.BOTCAN)
1719+
! TLCAN = Lagrangian timescale
1720+
TLCAN = (FCH/ustar(i)) * (
1721+
& (0.256 * (ZCAN-(0.75*FCH))/FCH ) +
1722+
& (0.492*EXP((-0.256*ZCAN/FCH)/0.492)) )
1723+
IF ( HOL .LT. -0.1 ) THEN !STRONG UNSTABLE
1724+
IF ( ZCAN/FCH .GT. 1.25 ) THEN !SIGMACAN = Eulerian vertical velocity variance
1725+
SIGMACAN = 1.25*ustar(i)
1726+
END IF
1727+
IF ( ZCAN/FCH .GE. 0.175
1728+
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
1729+
SIGMACAN = ustar(i) * ( 0.75 +
1730+
& (0.5 * COS((PI/1.06818) *
1731+
& (1.25 - (ZCAN/FCH)))) )
1732+
END IF
1733+
IF ( ZCAN/FCH .LT. 0.175 ) THEN
1734+
SIGMACAN = 0.25*ustar(i)
1735+
END IF
1736+
END IF
1737+
IF ( HOL .GE. -0.1 .AND. HOL .LT. 0.1 ) THEN !WEAKLY UNSTABLE to NEUTRAL
1738+
IF ( ZCAN/FCH .GT. 1.25 ) THEN
1739+
SIGMACAN = 1.0*ustar(i)
1740+
END IF
1741+
IF ( ZCAN/FCH .GE. 0.175
1742+
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
1743+
SIGMACAN = ustar(i) * ( 0.625 +
1744+
& (0.375* COS((PI/1.06818) *
1745+
& (1.25 - (ZCAN/FCH)))) )
1746+
END IF
1747+
IF ( ZCAN/FCH .LT. 0.175 ) THEN
1748+
SIGMACAN = 0.25*ustar(i)
1749+
END IF
1750+
END IF
1751+
IF ( HOL .GE. 0.1 .AND. HOL .LT. 0.9 ) THEN !STABLE
1752+
IF ( ZCAN/FCH .GT. 1.25 ) THEN
1753+
SIGMACAN = 0.25*(4.375 - (3.75*HOL))*ustar(i)
1754+
END IF
1755+
IF ( ZCAN/FCH .GE. 0.175
1756+
& .AND. ZCAN/FCH .LE. 1.25 ) THEN
1757+
RRCAN=4.375-(3.75*HOL)
1758+
AACAN=(0.125*RRCAN) + 0.125
1759+
BBCAN=(0.125*RRCAN) - 0.125
1760+
SIGMACAN = ustar(i) * ( AACAN +
1761+
& (BBCAN * COS((PI/1.06818) *
1762+
& (1.25 - (ZCAN/FCH)))) )
1763+
END IF
1764+
IF ( ZCAN/FCH .LT. 0.175 ) THEN
1765+
SIGMACAN = 0.25*ustar(i)
1766+
END IF
1767+
END IF
1768+
IF ( HOL .GE. 0.9 ) THEN !VERY STABLE
1769+
SIGMACAN = 0.25*ustar(i)
1770+
END IF
1771+
IF ( ZCAN .EQ. ZFL ) THEN ! Each model layer that includes canopy
1772+
EDDYVEST1 = (SIGMACAN*SIGMACAN)*TLCAN
1773+
ELSE IF ( ZCAN .LE. FCH ) THEN !in-canopy layers and set arrays
1774+
COUNTCAN = COUNTCAN + 1
1775+
ZCANX (COUNTCAN) = ZCAN
1776+
EDDYVESTX (COUNTCAN) = (SIGMACAN*SIGMACAN)*TLCAN
1777+
END IF
1778+
ZCAN = ZCAN-0.5 !step down in-canopy resolution of 0.5m
1779+
END DO !end loop on canopy layers
1780+
EDDYVEST_INT = IntegrateTrapezoid((ZCANX(COUNTCAN:1:-1)
1781+
& ),EDDYVESTX(COUNTCAN:1:-1)) / ZFL
1782+
dkt(i,k)= (dkt(i,k)/EDDYVEST1) * EDDYVEST_INT !Scale dkt to resolved eddy diffusivity
1783+
dkq(i,k)= (dkq(i,k)/EDDYVEST1) * EDDYVEST_INT !Scale dkq to resolved eddy diffusivity
1784+
dku(i,k)= (dku(i,k)/EDDYVEST1) * EDDYVEST_INT !Scale dku to resolved eddy diffusivity
1785+
1786+
!IVAI: Output contiguos canopy correction bottom layer and 3D
1787+
! if ( kount .EQ. 0)
1788+
! & aux2d(i,4) = 1./EDDYVEST1 * EDDYVEST_INT
1789+
! aux3d(i,k,4) = 1./EDDYVEST1 * EDDYVEST_INT
1790+
!IVAI
1791+
1792+
END IF ! contigous canopy conditions
1793+
1794+
END IF ! (KCAN .EQ. 1) model layer(s) containing canopy
1795+
1796+
enddo !i
1797+
1798+
kount = kount + 1 !IVAI
1799+
1800+
enddo !k
1801+
1802+
endif !do_canopy .and. cplaqm
1803+
15771804
!> ## Compute TKE.
15781805
!! - Compute a minimum TKE deduced from background diffusivity for momentum.
15791806
!
@@ -1995,6 +2222,7 @@ subroutine satmedmfvdifq_run(im,km,ntrac,ntcw,ntrw, &
19952222
tkei(i,k) = 0.5 * (tke(i,k)+tke(i,k+1))
19962223
enddo
19972224
enddo
2225+
19982226
do k=1,kps
19992227
do i=1,im
20002228
e_diff(i,k) = tke(i,k) - tke(i,k+1)

0 commit comments

Comments
 (0)