Skip to content

Commit 9f54c41

Browse files
dustinswalesdustinswales
committed
Merge branch 'ufs/dev' of https://github.com/ufs-community/ccpp-physics into feature/mpas_nuopc_init
2 parents 8cd11b0 + 5d6bacf commit 9f54c41

3 files changed

Lines changed: 332 additions & 9 deletions

File tree

‎physics/PBL/SATMEDMF/satmedmfvdifq.F‎

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

0 commit comments

Comments
 (0)