<p><b>laura@ucar.edu</b> 2012-10-04 10:01:06 -0600 (Thu, 04 Oct 2012)</p><p>updated the surface layer scheme and the YSU planetary boundary layer scheme to WRF versions 3.4.1<br>
</p><hr noshade><pre><font color="gray">Modified: branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_bl_ysu.F
===================================================================
--- branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_bl_ysu.F 2012-10-04 15:36:29 UTC (rev 2184)
+++ branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_bl_ysu.F 2012-10-04 16:01:06 UTC (rev 2185)
@@ -1,4 +1,4 @@
-!WRf:model_layer:physics
+!WRF:model_layer:physics
!
!
!
@@ -19,10 +19,11 @@
dz8w,psfc, &
znu,znw,mut,p_top, &
znt,ust,hpbl,psim,psih, &
- xland,hfx,qfx,gz1oz0,wspd,br, &
+ xland,hfx,qfx,wspd,br, &
dt,kpbl2d, &
exch_h, &
u10,v10, &
+ ctopo,ctopo2, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
@@ -74,7 +75,6 @@
!-- xland land mask (1 for land, 2 for water)
!-- hfx upward heat flux at the surface (w/m^2)
!-- qfx upward moisture flux at the surface (kg/m^2/s)
-!-- gz1oz0 log(z/z0) where z0 is roughness length
!-- wspd wind speed at lowest model level (m/s)
!-- u10 u-wind speed at 10 m (m/s)
!-- v10 v-wind speed at 10 m (m/s)
@@ -117,40 +117,39 @@
!
real, dimension( ims:ime, kms:kme, jms:jme ) , &
intent(in ) :: qv3d, &
- qc3d, &
- qi3d, &
- p3d, &
- pi3d, &
- th3d, &
- t3d, &
- dz8w
+ qc3d, &
+ qi3d, &
+ p3d, &
+ pi3d, &
+ th3d, &
+ t3d, &
+ dz8w
real, dimension( ims:ime, kms:kme, jms:jme ) , &
intent(in ) :: p3di
!
real, dimension( ims:ime, kms:kme, jms:jme ) , &
intent(inout) :: rublten, &
- rvblten, &
- rthblten, &
- rqvblten, &
+ rvblten, &
+ rthblten, &
+ rqvblten, &
rqcblten
!
real, dimension( ims:ime, kms:kme, jms:jme ) , &
intent(inout) :: exch_h
real, dimension( ims:ime, jms:jme ) , &
- intent(in ) :: u10, &
+ intent(inout) :: u10, &
v10
!
real, dimension( ims:ime, jms:jme ) , &
intent(in ) :: xland, &
- hfx, &
+ hfx, &
qfx, &
br, &
psfc
real, dimension( ims:ime, jms:jme ) , &
intent(in ) :: &
psim, &
- psih, &
- gz1oz0
+ psih
real, dimension( ims:ime, jms:jme ) , &
intent(inout) :: znt, &
ust, &
@@ -184,6 +183,10 @@
!
real, optional, intent(in ) :: p_top
!
+ real, dimension( ims:ime, jms:jme ) , &
+ optional , &
+ intent(in ) :: ctopo, &
+ ctopo2
!local
integer :: i,j,k
real, dimension( its:ite, kts:kte*ndiff ) :: rqvbl2dt, &
@@ -195,11 +198,10 @@
dvsfc, &
dtsfc, &
dqsfc
-
-!MPAS specific for debugging (Laura D. Fowler):
- real,intent(in) ,dimension(ims:ime,kms:kme,jms:jme),optional:: rho
+!MPAS specific (Laura D. Fowler):
+ real,intent(in),dimension(ims:ime,kms:kme,jms:jme),optional:: rho
+ real:: rho_d
real,intent(out),dimension(ims:ime,kms:kme,jms:jme),optional:: kzhout,kzmout,kzqout
- real:: rho_d
do j = jts,jte
do k = kts,kte
do i = its,ite
@@ -279,10 +281,12 @@
,dt=dt,rcl=1.0,kpbl1d=kpbl2d(ims,j) &
,exch_hx=exch_h(ims,kms,j) &
,u10=u10(ims,j),v10=v10(ims,j) &
- ,gz1oz0=gz1oz0(ims,j) &
,kzh=kzhout(ims,kms,j) &
,kzm=kzmout(ims,kms,j) &
,kzq=kzqout(ims,kms,j) &
+#if ( ! NMM_CORE == 1 )
+ ,ctopo=ctopo(ims,j),ctopo2=ctopo2(ims,j) &
+#endif
,ids=ids,ide=ide, jds=jds,jde=jde, kds=kds,kde=kde &
,ims=ims,ime=ime, jms=jms,jme=jme, kms=kms,kme=kme &
,its=its,ite=ite, jts=jts,jte=jte, kts=kts,kte=kte )
@@ -311,7 +315,7 @@
dt,rcl,kpbl1d, &
exch_hx, &
u10,v10, &
- gz1oz0, &
+ ctopo,ctopo2, &
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
@@ -349,11 +353,21 @@
! pressure-level diffusion, april 2009
! ==> negligible differences
! implicit forcing for momentum with clean up, july 2009
-! ==> prevents model blownup when sfc layer is too low
-! increase of lamda, 30 < 0.1 x del z < 300, feb 2010
+! ==> prevents model blowup when sfc layer is too low
+! incresea of lamda, maximum (30, 0.1 x del z) feb 2010
! ==> prevents model blowup when delz is extremely large
! revised prandtl number at surface, peggy lemone, feb 2010
! ==> increase kh, decrease mixing due to counter-gradient term
+! revised thermal, shin et al. mon. wea. rev. , songyou hong, aug 2011
+! ==> reduce the thermal strength when z1 < 0.1 h
+! revised prandtl number for free convection, dudhia, mar 2012
+! ==> pr0 = 1 + bke (=0.272) when newtral, kh is reduced
+! minimum kzo = 0.01, lo = min (30m,delz), hong, mar 2012
+! ==> weaker mixing when stable, and les resolution in vertical
+! gz1oz0 is removed, and phim phih are ln(z1/z0)-phim,h, hong, mar 2012
+! ==> consider thermal z0 when differs from mechanical z0
+! a bug fix in wscale computation in stable bl, sukanta basu, jun 2012
+! ==> wscale becomes small with height, and less mixing in stable bl
!
! references:
!
@@ -425,7 +439,7 @@
!
real, dimension( ims:ime ), intent(in ) :: psim, &
psih
- real, dimension( ims:ime ), intent(in ) :: gz1oz0
+
!
real, dimension( ims:ime ), intent(in ) :: psfcpa
integer, dimension( ims:ime ), intent(out ) :: kpbl1d
@@ -433,9 +447,12 @@
real, dimension( ims:ime, kms:kme ) , &
intent(in ) :: ux, &
vx
-!optional
real, dimension( ims:ime ) , &
optional , &
+ intent(in ) :: ctopo, &
+ ctopo2
+ real, dimension( ims:ime ) , &
+ optional , &
intent(inout) :: regime
!
! local vars
@@ -443,29 +460,30 @@
real, dimension( its:ite ) :: hol
real, dimension( its:ite, kts:kte+1 ) :: zq
!
- real, dimension( its:ite, kts:kte ) :: &
+ real, dimension( its:ite, kts:kte ) :: &
thx,thvx, &
del, &
dza, &
dzq, &
+ xkzo, &
za
!
real, dimension( its:ite ) :: &
rhox, &
govrth, &
zl1,thermal, &
- wscale,hgamt, &
- hgamq,brdn, &
- brup,phim, &
- phih, &
+ wscale, &
+ hgamt,hgamq, &
+ brdn,brup, &
+ phim,phih, &
dusfc,dvsfc, &
dtsfc,dqsfc, &
prpbl, &
wspd1
!
real, dimension( its:ite, kts:kte ) :: xkzm,xkzh, &
- f1,f2, &
- r1,r2, &
+ f1,f2, &
+ r1,r2, &
ad,au, &
cu, &
al, &
@@ -476,8 +494,8 @@
real, dimension( ims:ime, kms:kme ) , &
intent(inout) :: exch_hx
!
- real, dimension( ims:ime ) , &
- intent(in ) :: u10, &
+ real, dimension( ims:ime ) , &
+ intent(inout) :: u10, &
v10
real, dimension( its:ite ) :: &
brcr, &
@@ -496,14 +514,15 @@
!
!
real :: dt2,rdt,spdk2,fm,fh,hol1,gamfac,vpert,prnum,prnum0
- real :: xkzo,ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
+ real :: ss,ri,qmean,tmean,alph,chi,zk,rl2,dk,sri
real :: brint,dtodsd,dtodsu,rdz,dsdzt,dsdzq,dsdz2,rlamdz
real :: utend,vtend,ttend,qtend
real :: dtstep,govrthv
real :: cont, conq, conw, conwrc
!
- real, dimension( its:ite, kts:kte ) :: wscalek, &
- xkzml,xkzhl, &
+ real, dimension( its:ite, kts:kte ) :: wscalek
+ real, dimension( its:ite ) :: delta
+ real, dimension( its:ite, kts:kte ) :: xkzml,xkzhl, &
zfacent,entfac
real, dimension( its:ite ) :: ust3, &
wstar3,wstar, &
@@ -512,12 +531,15 @@
bfxpbl, &
hfxpbl,qfxpbl, &
ufxpbl,vfxpbl, &
- delta,dthvx
+ dthvx, &
+ zol1
real :: prnumfac,bfx0,hfx0,qfx0,delb,dux,dvx, &
- dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr,prfac
+ dsdzu,dsdzv,wm3,dthx,dqx,wspd10,ross,tem1,dsig,tvcon,conpr, &
+ prfac,prfac2,phim8z
!
-!ldf:
+!MPAS specific begin (Laura Fowler - 2012-08-24):
real,intent(out),dimension(ims:ime,kms:kme),optional::kzh,kzm,kzq
+!MPAS specific end.
!----------------------------------------------------------------------
!
@@ -628,6 +650,11 @@
zfac(i,k) = 0.0
enddo
enddo
+ do k = kts,klpbl-1
+ do i = its,ite
+ xkzo(i,k) = ckz*dza(i,k+1)
+ enddo
+ enddo
!
do i = its,ite
dusfc(i) = 0.
@@ -685,16 +712,15 @@
enddo
!
do i = its,ite
- fm = gz1oz0(i)-psim(i)
- fh = gz1oz0(i)-psih(i)
- hol(i) = max(br(i)*fm*fm/fh,rimin)
+ fm = psim(i)
+ fh = psih(i)
+ zol1(i) = max(br(i)*fm*fm/fh,rimin)
if(sfcflg(i))then
- hol(i) = min(hol(i),-zfmin)
+ zol1(i) = min(zol1(i),-zfmin)
else
- hol(i) = max(hol(i),zfmin)
+ zol1(i) = max(zol1(i),zfmin)
endif
- hol1 = hol(i)*hpbl(i)/zl1(i)*sfcfrac
- hol(i) = -hol(i)*hpbl(i)/zl1(i)
+ hol1 = zol1(i)*hpbl(i)/zl1(i)*sfcfrac
if(sfcflg(i))then
phim(i) = (1.-aphi16*hol1)**(-1./4.)
phih(i) = (1.-aphi16*hol1)**(-1./2.)
@@ -724,7 +750,7 @@
hgamt(i) = min(gamfac*hfx(i)/cp,gamcrt)
hgamq(i) = min(gamfac*qfx(i),gamcrq)
vpert = (hgamt(i)+ep1*thx(i,1)*hgamq(i))/bfac*afac
- thermal(i) = thermal(i)+max(vpert,0.)
+ thermal(i) = thermal(i)+max(vpert,0.)*min(za(i,1)/(sfcfrac*hpbl(i)),1.0)
hgamt(i) = max(hgamt(i),0.0)
hgamq(i) = max(hgamq(i),0.0)
brint = -15.9*ust(i)*ust(i)/wspd(i)*wstar3(i)/(wscale(i)**4.)
@@ -890,31 +916,34 @@
do i = its,ite
if(k.lt.kpbl(i)) then
zfac(i,k) = min(max((1.-(zq(i,k+1)-zl1(i))/(hpbl(i)-zl1(i))),zfmin),1.)
- xkzo = ckz*dza(i,k+1)
zfacent(i,k) = (1.-zfac(i,k))**3.
+ wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
if(sfcflg(i)) then
- prfac = conpr/phim(i)/(1.+4.*karman*wstar3(i)/ust3(i))
+ prfac = conpr
+ prfac2 = 15.9*wstar3(i)/ust3(i)/(1.+4.*karman*wstar3(i)/ust3(i))
prnumfac = -3.*(max(zq(i,k+1)-sfcfrac*hpbl(i),0.))**2./hpbl(i)**2.
else
prfac = 0.
+ prfac2 = 0.
prnumfac = 0.
+ phim8z = 1.+aphi5*zol1(i)*zq(i,k+1)/zl1(i)
+ wscalek(i,k) = ust(i)/phim8z
+ wscalek(i,k) = max(wscalek(i,k),0.001)
endif
prnum0 = (phih(i)/phim(i)+prfac)
- prnum0 = min(prnum0,prmax)
- prnum0 = max(prnum0,prmin)
- wscalek(i,k) = (ust3(i)+phifac*karman*wstar3(i)*(1.-zfac(i,k)))**h1
- xkzm(i,k) = xkzo+wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac
+ prnum0 = max(min(prnum0,prmax),prmin)
+ xkzm(i,k) = wscalek(i,k)*karman*zq(i,k+1)*zfac(i,k)**pfac
prnum = 1. + (prnum0-1.)*exp(prnumfac)
xkzq(i,k) = xkzm(i,k)/prnum*zfac(i,k)**(pfac_q-pfac)
- prnum0 = prnum0/(1.+prfac)
+ prnum0 = prnum0/(1.+prfac2*karman*sfcfrac)
prnum = 1. + (prnum0-1.)*exp(prnumfac)
xkzh(i,k) = xkzm(i,k)/prnum
xkzm(i,k) = min(xkzm(i,k),xkzmax)
- xkzm(i,k) = max(xkzm(i,k),xkzmin)
+ xkzm(i,k) = max(xkzm(i,k),xkzo(i,k))
xkzh(i,k) = min(xkzh(i,k),xkzmax)
- xkzh(i,k) = max(xkzh(i,k),xkzmin)
+ xkzh(i,k) = max(xkzh(i,k),xkzo(i,k))
xkzq(i,k) = min(xkzq(i,k),xkzmax)
- xkzq(i,k) = max(xkzq(i,k),xkzmin)
+ xkzq(i,k) = max(xkzq(i,k),xkzo(i,k))
endif
enddo
enddo
@@ -923,7 +952,6 @@
!
do k = kts,kte-1
do i = its,ite
- xkzo = ckz*dza(i,k+1)
if(k.ge.kpbl(i)) then
ss = ((ux(i,k+1)-ux(i,k))*(ux(i,k+1)-ux(i,k)) &
+(vx(i,k+1)-vx(i,k))*(vx(i,k+1)-vx(i,k))) &
@@ -943,25 +971,26 @@
endif
zk = karman*zq(i,k+1)
rlamdz = min(max(0.1*dza(i,k+1),rlam),300.)
+ rlamdz = min(dza(i,k+1),rlamdz)
rl2 = (zk*rlamdz/(rlamdz+zk))**2
dk = rl2*sqrt(ss)
if(ri.lt.0.)then
! unstable regime
sri = sqrt(-ri)
- xkzm(i,k) = xkzo+dk*(1+8.*(-ri)/(1+1.746*sri))
- xkzh(i,k) = xkzo+dk*(1+8.*(-ri)/(1+1.286*sri))
+ xkzm(i,k) = dk*(1+8.*(-ri)/(1+1.746*sri))
+ xkzh(i,k) = dk*(1+8.*(-ri)/(1+1.286*sri))
else
! stable regime
- xkzh(i,k) = xkzo+dk/(1+5.*ri)**2
+ xkzh(i,k) = dk/(1+5.*ri)**2
prnum = 1.0+2.1*ri
prnum = min(prnum,prmax)
- xkzm(i,k) = (xkzh(i,k)-xkzo)*prnum+xkzo
+ xkzm(i,k) = xkzh(i,k)*prnum
endif
!
xkzm(i,k) = min(xkzm(i,k),xkzmax)
- xkzm(i,k) = max(xkzm(i,k),xkzmin)
+ xkzm(i,k) = max(xkzm(i,k),xkzo(i,k))
xkzh(i,k) = min(xkzh(i,k),xkzmax)
- xkzh(i,k) = max(xkzh(i,k),xkzmin)
+ xkzh(i,k) = max(xkzh(i,k),xkzo(i,k))
xkzml(i,k) = xkzm(i,k)
xkzhl(i,k) = xkzh(i,k)
endif
@@ -999,7 +1028,7 @@
xkzh(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
xkzh(i,k) = sqrt(xkzh(i,k)*xkzhl(i,k))
xkzh(i,k) = min(xkzh(i,k),xkzmax)
- xkzh(i,k) = max(xkzh(i,k),xkzmin)
+ xkzh(i,k) = max(xkzh(i,k),xkzo(i,k))
f1(i,k+1) = thx(i,k+1)-300.
else
f1(i,k+1) = thx(i,k+1)-300.
@@ -1035,7 +1064,7 @@
enddo
enddo
!
-! compute tridiagonal matrix elements for moisture, clouds, and tracers
+! compute tridiagonal matrix elements for moisture, clouds, and gases
!
do k = kts,kte
do i = its,ite
@@ -1090,7 +1119,7 @@
xkzq(i,k) = -we(i)*dza(i,kpbl(i))*exp(-entfac(i,k))
xkzq(i,k) = sqrt(xkzq(i,k)*xkzhl(i,k))
xkzq(i,k) = min(xkzq(i,k),xkzmax)
- xkzq(i,k) = max(xkzq(i,k),xkzmin)
+ xkzq(i,k) = max(xkzq(i,k),xkzo(i,k))
f3(i,k+1,1) = qx(i,k+1)
else
f3(i,k+1,1) = qx(i,k+1)
@@ -1132,7 +1161,7 @@
enddo
enddo
!
-! solve tridiagonal problem for moisture, clouds, and tracers
+! solve tridiagonal problem for moisture, clouds, and gases
!
call tridin_ysu(al,ad,cu,r3,au,f3,its,ite,kts,kte,ndiff)
!
@@ -1171,8 +1200,15 @@
enddo
!
do i = its,ite
- ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
- *(wspd1(i)/wspd(i))**2
+! paj: ctopo=1 if topo_wind=0 (default)
+! mchen add this line to make sure NMM can still work with YSU PBL
+ if(present(ctopo)) then
+ ad(i,1) = 1.+ctopo(i)*ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
+ *(wspd1(i)/wspd(i))**2
+ else
+ ad(i,1) = 1.+ust(i)**2/wspd1(i)*rhox(i)*g/del(i,1)*dt2 &
+ *(wspd1(i)/wspd(i))**2
+ endif
f1(i,1) = ux(i,1)
f2(i,1) = vx(i,1)
enddo
@@ -1195,7 +1231,7 @@
xkzm(i,k) = prpbl(i)*xkzh(i,k)
xkzm(i,k) = sqrt(xkzm(i,k)*xkzml(i,k))
xkzm(i,k) = min(xkzm(i,k),xkzmax)
- xkzm(i,k) = max(xkzm(i,k),xkzmin)
+ xkzm(i,k) = max(xkzm(i,k),xkzo(i,k))
f1(i,k+1) = ux(i,k+1)
f2(i,k+1) = vx(i,k+1)
else
@@ -1238,13 +1274,22 @@
enddo
enddo
!
+! paj: ctopo2=1 if topo_wind=0 (default)
+!
+ do i = its,ite
+ if(present(ctopo).and.present(ctopo2)) then ! mchen for NMM
+ u10(i) = ctopo2(i)*u10(i)+(1-ctopo2(i))*ux(i,1)
+ v10(i) = ctopo2(i)*v10(i)+(1-ctopo2(i))*vx(i,1)
+ endif !mchen
+ enddo
+!
!---- end of vertical diffusion
!
do i = its,ite
kpbl1d(i) = kpbl(i)
enddo
!
-!MPAS specific (Laura D. Fowler)::
+!MPAS specific begin (Laura D. Fowler - 2012-08-24)::
if(present(kzh) .and. present(kzm) .and. present(kzq)) then
do i = its,ite
do k = kts,kte
Modified: branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_sf_sfclay.F
===================================================================
--- branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_sf_sfclay.F 2012-10-04 15:36:29 UTC (rev 2184)
+++ branches/atmos_physics/src/core_atmos_physics/physics_wrf/module_sf_sfclay.F 2012-10-04 16:01:06 UTC (rev 2185)
@@ -14,6 +14,7 @@
SUBROUTINE SFCLAY(U3D,V3D,T3D,QV3D,P3D,dz8w, &
CP,G,ROVCP,R,XLV,PSFC,CHS,CHS2,CQS2,CPM, &
ZNT,UST,PBLH,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
+ FM,FH, &
XLAND,HFX,QFX,LH,TSK,FLHC,FLQC,QGH,QSFC,RMOL, &
U10,V10,TH2,T2,Q2, &
GZ1OZ0,WSPD,BR,ISFFLX,DX, &
@@ -23,7 +24,8 @@
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
- ustm,ck,cka,cd,cda,isftcflx,iz0tlnd,areaCell )
+ ustm,ck,cka,cd,cda,isftcflx,iz0tlnd, &
+ scm_force_flux,areaCell)
!-------------------------------------------------------------------
IMPLICIT NONE
!-------------------------------------------------------------------
@@ -50,6 +52,8 @@
!-- REGIME flag indicating PBL regime (stable, unstable, etc.)
!-- PSIM similarity stability function for momentum
!-- PSIH similarity stability function for heat
+!-- FM integrated stability function for momentum
+!-- FH integrated stability function for heat
!-- XLAND land mask (1 for land, 2 for water)
!-- HFX upward heat flux at the surface (W/m^2)
!-- QFX upward moisture flux at the surface (kg/m^2/s)
@@ -146,7 +150,7 @@
!
REAL, DIMENSION( ims:ime, jms:jme ) , &
INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
- PSIM,PSIH
+ PSIM,PSIH,FM,FH
REAL, DIMENSION( ims:ime, kms:kme, jms:jme ) , &
INTENT(IN ) :: U3D, &
@@ -184,6 +188,7 @@
!MPAS specific end.
INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
+ INTEGER, OPTIONAL, INTENT(IN ) :: SCM_FORCE_FLUX
! LOCAL VARS
@@ -218,6 +223,7 @@
CQS2(ims,j),CPM(ims,j),PBLH(ims,j), RMOL(ims,j), &
ZNT(ims,j),UST(ims,j),MAVAIL(ims,j),ZOL(ims,j), &
MOL(ims,j),REGIME(ims,j),PSIM(ims,j),PSIH(ims,j), &
+ FM(ims,j),FH(ims,j), &
XLAND(ims,j),HFX(ims,j),QFX(ims,j),TSK(ims,j), &
U10(ims,j),V10(ims,j),TH2(ims,j),T2(ims,j), &
Q2(ims,j),FLHC(ims,j),FLQC(ims,j),QGH(ims,j), &
@@ -230,25 +236,25 @@
its,ite, jts,jte, kts,kte &
#if defined(non_hydrostatic_core) || defined(hydrostatic_core)
!MPAS specific (Laura D. Fowler):
- ,isftcflx,iz0tlnd, &
+ ,isftcflx,iz0tlnd,scm_force_flux, &
USTM(ims,j),CK(ims,j),CKA(ims,j), &
CD(ims,j),CDA(ims,j),areaCell(ims,j) &
-!#elseif ( EM_CORE == 1 )
-! ,isftcflx,iz0tlnd, &
-! USTM(ims,j),CK(ims,j),CKA(ims,j), &
-! CD(ims,j),CDA(ims,j) &
+#elif ( EM_CORE == 1 )
+ ,isftcflx,iz0tlnd,scm_force_flux, &
+ USTM(ims,j),CK(ims,j),CKA(ims,j), &
+ CD(ims,j),CDA(ims,j) &
#endif
)
ENDDO
-
+
END SUBROUTINE SFCLAY
!-------------------------------------------------------------------
SUBROUTINE SFCLAY1D(J,UX,VX,T1D,QV1D,P1D,dz8w1d, &
CP,G,ROVCP,R,XLV,PSFCPA,CHS,CHS2,CQS2,CPM,PBLH,RMOL, &
- ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH, &
+ ZNT,UST,MAVAIL,ZOL,MOL,REGIME,PSIM,PSIH,FM,FH,&
XLAND,HFX,QFX,TSK, &
U10,V10,TH2,T2,Q2,FLHC,FLQC,QGH, &
QSFC,LH,GZ1OZ0,WSPD,BR,ISFFLX,DX, &
@@ -258,7 +264,7 @@
ids,ide, jds,jde, kds,kde, &
ims,ime, jms,jme, kms,kme, &
its,ite, jts,jte, kts,kte, &
- isftcflx, iz0tlnd, &
+ isftcflx, iz0tlnd, scm_force_flux, &
ustm,ck,cka,cd,cda, &
areaCell)
!-------------------------------------------------------------------
@@ -296,7 +302,7 @@
!
REAL, DIMENSION( ims:ime ) , &
INTENT(INOUT) :: GZ1OZ0,WSPD,BR, &
- PSIM,PSIH
+ PSIM,PSIH,FM,FH
REAL, DIMENSION( ims:ime ) , &
INTENT(INOUT) :: ZNT, &
@@ -334,6 +340,7 @@
INTENT(OUT) :: ck,cka,cd,cda,ustm
INTEGER, OPTIONAL, INTENT(IN ) :: ISFTCFLX, IZ0TLND
+ INTEGER, OPTIONAL, INTENT(IN ) :: SCM_FORCE_FLUX
!MPAS specific (Laura D. Fowler): We take into accound the actual size of individual
!grid-boxes:
@@ -372,6 +379,7 @@
REAL :: ZL,TSKV,DTHVDZ,DTHVM,VCONV,RZOL,RZOL2,RZOL10,ZOL2,ZOL10
REAL :: DTG,PSIX,DTTHX,PSIX10,PSIT,PSIT2,PSIQ,PSIQ2,PSIQ10
REAL :: FLUXC,VSGD,Z0Q,VISC,RESTAR,CZIL,RESTAR2
+ REAL :: ZW, ZN1, ZN2
!-------------------------------------------------------------------
KL=kte
@@ -511,8 +519,6 @@
else
VSGD = 0.32 * (max(dx/5000.-1.,0.))**.33
endif
-! write(0,201) i,areaCell(i),vsgd
-! 201 format(i8,2(1x,e15.8))
!MPAS specific end.
WSPD(I)=SQRT(WSPD(I)*WSPD(I)+VCONV*VCONV+vsgd*vsgd)
WSPD(I)=AMAX1(WSPD(I),0.1)
@@ -523,7 +529,7 @@
RMOL(I)=-GOVRTH(I)*DTHVDZ*ZA(I)*KARMAN
!jdf
- 260 CONTINUE
+ 260 CONTINUE
!
!-----DIAGNOSE BASIC PARAMETERS FOR THE APPROPRIATED STABILITY CLASS:
@@ -782,20 +788,25 @@
DENOMQ(I)=PSIQ
DENOMQ2(I)=PSIQ2
DENOMT2(I)=PSIT2
+ FM(I)=PSIX
+ FH(I)=PSIT
330 CONTINUE
!
335 CONTINUE
!-----COMPUTE THE SURFACE SENSIBLE AND LATENT HEAT FLUXES:
-
+ IF ( PRESENT(SCM_FORCE_FLUX) ) THEN
+ IF (SCM_FORCE_FLUX.EQ.1) GOTO 350
+ ENDIF
DO i=its,ite
QFX(i)=0.
HFX(i)=0.
ENDDO
+ 350 CONTINUE
IF (ISFFLX.EQ.0) GOTO 410
-!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
+!-----OVER WATER, ALTER ROUGHNESS LENGTH (ZNT) ACCORDING TO WIND (UST).
DO 360 I=its,ite
IF((XLAND(I)-1.5).GE.0)THEN
@@ -804,8 +815,16 @@
IF ( PRESENT(ISFTCFLX) ) THEN
IF ( ISFTCFLX.NE.0 ) THEN
! ZNT(I)=10.*exp(-9.*UST(I)**(-.3333))
- ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333))
- ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
+! ZNT(I)=10.*exp(-9.5*UST(I)**(-.3333))
+! ZNT(I)=ZNT(I) + 0.11*1.5E-5/AMAX1(UST(I),0.01)
+! ZNT(I)=0.011*UST(I)*UST(I)/G+OZO
+! ZNT(I)=MAX(ZNT(I),3.50e-5)
+! AHW 2012:
+ ZW = MIN((UST(I)/1.06)**(0.3),1.0)
+ ZN1 = 0.011*UST(I)*UST(I)/G + OZO
+ ZN2 = 10.*exp(-9.5*UST(I)**(-.3333)) + &
+ 0.11*1.5E-5/AMAX1(UST(I),0.01)
+ ZNT(I)=(1.0-ZW) * ZN1 + ZW * ZN2
ZNT(I)=MIN(ZNT(I),2.85e-3)
ZNT(I)=MAX(ZNT(I),1.27e-7)
ENDIF
@@ -825,13 +844,16 @@
ELSE
FLHC(I)=0.
ENDIF
- 360 CONTINUE
+ 360 CONTINUE
!
-!-----COMPUTE SURFACE MOIST FLUX:
-!
-! IF(IDRY.EQ.1)GOTO 390
-!
+!-----COMPUTE SURFACE MOIST FLUX:
+!
+! IF(IDRY.EQ.1)GOTO 390
+ IF ( PRESENT(SCM_FORCE_FLUX) ) THEN
+ IF (SCM_FORCE_FLUX.EQ.1) GOTO 405
+ ENDIF
+!
DO 370 I=its,ite
QFX(I)=FLQC(I)*(QSFC(I)-QX(I))
QFX(I)=AMAX1(QFX(I),0.)
@@ -855,6 +877,8 @@
HFX(I)=AMAX1(HFX(I),-250.)
ENDIF
400 CONTINUE
+
+ 405 CONTINUE
DO I=its,ite
IF((XLAND(I)-1.5).GE.0)THEN
</font>
</pre>