114 SUBROUTINE jsfc(FLAG_ITER,ITER,ME &
115 & ,NTSD,EPSQ2,HT,DZ &
116 & ,PHMID,PHINT,TH,T,Q,QC,U,V,Q2 &
117 & ,TSK,QSFC,THZ0,QZ0,UZ0,VZ0 &
119 & ,USTAR,Z0,Z0BASE,PBLH,MAVAIL,RMOL &
120 & ,AKHS,AKMS,CHKLOWQ,HLFLX,RIB &
121 & ,CM,CH,STRESS,FFM,FFH,WIND,FM10,FH2 &
123 & ,IDS,IDE,JDS,JDE,KDS,KDE &
124 & ,IMS,IME,JMS,JME,KMS,KME &
125 & ,ITS,ITE,JTS,JTE,KTS,LM,errmsg,errflg)
147 INTEGER,
INTENT(IN) :: ids,ide,jds,jde,kds,kde &
148 & ,IMS,IME,JMS,JME,KMS,KME &
149 & ,ITS,ITE,JTS,JTE,KTS,LM
151 INTEGER,
INTENT(IN) :: ntsd,iter,me
152 LOGICAL,
DIMENSION(IMS:IME,JMS:JME),
INTENT(IN) :: flag_iter
153 real(kind=kfpt),
dimension(1:lm),
intent(in):: epsq2
155 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME),
INTENT(IN) :: ht,mavail,tsk &
159 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME,1:LM),
INTENT(IN) :: dz,phmid
161 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME,1:LM+1),
INTENT(IN) :: phint
163 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME,1:LM),
INTENT(IN) :: q,qc,u,v,q2,t,th
169 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME) :: flx_lh,hfx &
170 & ,qfx,q10,th10,t02 &
172 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME) :: pshltr,qshltr,tshltr
174 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME),
INTENT(INOUT) :: akhs,akms &
177 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME),
INTENT(INOUT) :: qz0,rmol,thz0 &
181 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME),
INTENT(OUT) :: hlflx,chklowq
182 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME),
INTENT(OUT) :: cm,ch,stress,ffm &
183 & ,ffh,wind,fm10,fh2 &
185 character(len=*),
intent(out) :: errmsg
186 integer,
intent(out) :: errflg
191 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME) :: chs,chs2,cqs2 &
193 REAL(kind=kfpt),
DIMENSION(IMS:IME,JMS:JME) :: qgh,cpm,ct
198 INTEGER :: i,j,k,lmh,lpbl
200 REAL(kind=kfpt) :: a,apesfc,b,btgx,cwmlow &
201 & ,dqdt,dtdif,dtdt,dudt,dvdt &
203 & ,p02p,p10p,plow,psfc,ptop,qlow,qs02,qs10 &
204 & ,rapa,rapa02,rapa10,ratiomx,rdz,seamask,sm &
205 & ,t02p,t10p,tem,th02p,th10p,thlow,thelow,thm &
206 & ,tlow,tz0,ulow,vlow,zsl
208 REAL(kind=kfpt),
DIMENSION(KTS:LM) :: cwmk,pk,q2k,qk,thek,thk,tk,uk,vk
210 REAL(kind=kfpt),
DIMENSION(KTS:LM-1) :: el,elm
212 REAL(kind=kfpt),
DIMENSION(KTS:LM+1) :: zhk
214 REAL(kind=kfpt),
DIMENSION(ITS:ITE,JTS:JTE) :: thsk
216 REAL(kind=kfpt),
DIMENSION(ITS:ITE,JTS:JTE,KTS:LM+1) :: zint
230 IF(flag_iter(i,j))
THEN
240 IF(flag_iter(i,j))
THEN
241 zint(i,j,lm+1)=ht(i,j)
256 IF(flag_iter(i,j))
THEN
258 zint(i,j,k)=zint(i,j,k+1)+dz(i,j,k)
264 IF(ntsd==0.and.iter==1)
then
267 IF(flag_iter(i,j))
THEN
294 setup_integration:
DO j=jts,jte
298 IF(flag_iter(i,j))
THEN
305 psfc=phint(i,j,lmh+1)
307 thsk(i,j)=tsk(i,j)/(psfc*1.e-5)**cappa
311 seamask=xland(i,j)-1.
321 thek(k)=(cwmk(k)*(-elocp/tk(k))+1.)*thk(k)
341 IF(q2k(k)<=epsq2(k)*fh)
THEN
353 110 pblh(i,j)=zhk(lpbl)-zhk(lmh+1)
356 IF(qc(i,j,lm).GT.epsq)
THEN
373 zsl=(zhk(lmh)-zhk(lmh+1))*0.5
375 apesfc=(psfc*1.e-5)**cappa
382 CALL sfcdif(ntsd,seamask,thsk(i,j),qsfc(i,j),psfc &
383 & ,uz0(i,j),vz0(i,j),tz0,thz0(i,j),qz0(i,j) &
384 & ,ustar(i,j),z0(i,j),z0base(i,j),ct(i,j),rmol(i,j) &
385 & ,akms(i,j),akhs(i,j),pblh(i,j),mavail(i,j) &
386 & ,chs(i,j),chs2(i,j),cqs2(i,j) &
387 & ,hfx(i,j),qfx(i,j),flx_lh(i,j) &
388 & ,flhc(i,j),flqc(i,j),qgh(i,j),cpm(i,j) &
389 & ,ulow,vlow,tlow,thlow,thelow,qlow,cwmlow &
390 & ,zsl,plow,hlflx(i,j) &
392 & ,u10(i,j),v10(i,j),tshltr(i,j),th10(i,j) &
393 & ,qshltr(i,j),q10(i,j),pshltr(i,j) &
394 & ,ffm(i,j),ffh(i,j),fm10(i,j),fh2(i,j) &
395 & ,a1u(i,j),a1t(i,j),a1q(i,j) &
396 & ,ids,ide,jds,jde,kds,kde &
397 & ,ims,ime,jms,jme,kms,kme &
398 & ,its,ite,jts,jte,kts,lm,i,j,zhk(lmh+1),rib(i,j) &
406 th02(i,j)=tshltr(i,j)
408 rapa02=rapa-gocp02/th02p
409 rapa10=rapa-gocp10/th10p
414 t02(i,j) = th02(i,j)*apesfc
416 p02p=(rapa02**rcap)*1.e5
417 p10p=(rapa10**rcap)*1.e5
419 qs02=pq0/p02p*exp(a2*(t02p-a3)/(t02p-a4))
420 qs10=pq0/p10p*exp(a2*(t10p-a3)/(t10p-a4))
422 IF(qshltr(i,j)>qs02)qshltr(i,j)=qs02
423 IF(q10(i,j)>qs10)q10(i,j)=qs10
424 q02(i,j)=qshltr(i,j)/(1.-qshltr(i,j))
427 wind(i,j)=max(ustar(i,j)*ffm(i,j)/vkarman,1.0)
428 cm(i,j)=vkarman*vkarman/(ffm(i,j)*ffm(i,j))
429 ch(i,j)=vkarman*vkarman/(ffm(i,j)*ffh(i,j))
430 tem=0.00001/dz(i,j,lm)
431 cm(i,j)=max(cm(i,j),tem)
432 ch(i,j)=max(ch(i,j),tem)
433 stress(i,j)=cm(i,j) * wind(i,j) * wind(i,j)
434 ustar(i,j)=sqrt(stress(i,j))
442 ENDDO setup_integration
452 SUBROUTINE sfcdif(NTSD,SEAMASK,THS,QS,PSFC &
453 & ,UZ0,VZ0,TZ0,THZ0,QZ0 &
454 & ,USTAR,Z0,Z0BASE,CT,RLMO,AKMS,AKHS,PBLH,WETM &
455 & ,CHS,CHS2,CQS2,HFX,QFX,FLX_LH,FLHC,FLQC,QGH,CPM &
456 & ,ULOW,VLOW,TLOW,THLOW,THELOW,QLOW,CWMLOW &
458! & ,VEGF,SNOC & !added 5/17/2013
459 & ,U10,V10,TH02,TH10,Q02,Q10,PSHLTR &
460 & ,FFM,FFH,FM10,FH2,A1U,A1T,A1Q &
461 & ,IDS,IDE,JDS,JDE,KDS,KDE &
462 & ,IMS,IME,JMS,JME,KMS,KME &
463 & ,ITS,ITE,JTS,JTE,KTS,LM,I,J,ZSFC,RIB & ! Added Bulk Richardson No.
475 INTEGER,
INTENT(IN) :: IDS,IDE,JDS,JDE,KDS,KDE &
476 & ,IMS,IME,JMS,JME,KMS,KME &
477 & ,ITS,ITE,JTS,JTE,KTS,LM,i,j
479 INTEGER,
INTENT(IN) :: NTSD
481 REAL(kind=kfpt),
INTENT(IN) :: cwmlow,pblh,plow,qlow,psfc,seamask,zsfc &
482 & ,thelow,thlow,ths,tlow,tz0,ulow,vlow,wetm,zsl &
486 REAL(kind=kfpt),
INTENT(OUT) :: chs,chs2,cpm,cqs2,ct,flhc,flqc,flx_lh,hfx &
487 & ,rib,pshltr,q02,q10,qfx,qgh,rlmo,th02,th10,u10,v10
488 REAL(kind=kfpt),
INTENT(OUT) :: ffm,ffh,fm10,fh2,a1u,a1t,a1q
490 REAL(kind=kfpt),
INTENT(INOUT) :: akhs,akms,qz0,thz0,ustar,uz0,vz0,z0,qs
491 character(len=*),
intent(out) :: errmsg
492 integer,
intent(out) :: errflg
499 REAL(kind=kfpt) :: a,b,btgh,btgx,cxchl,cxchs,dthv,du2,elfc,fct &
500 & ,hlflx,hsflx,hv,psh02,psh10,pshz,pshzl,psm10,psmz,psmzl &
501 & ,rdz,rdzt,rlma,rlmn,rlmp &
502 & ,rlogt,rlogu,rwgh,rz,rzst,rzsu,simh,simm,tem,thm &
503 & ,umflx,ustark,vmflx,wght,wghtt,wghtq,wstar2 &
504 & ,x,xlt,xlt4,xlu,xlu4,xt,xt4,xu,xu4,zetalt,zetalu &
505 & ,zetat,zetau,zq,zslt,zslu,zt,zu,topoterm,zzil,czetmax
510 REAL(kind=kfpt) :: akhs02,akhs10,akms02,akms10,ekms10,qsat10,qsat2 &
511 & ,rlnt02,rlnt10,rlnu10,simh02,simh10,simm10,t02,t10 &
512 & ,term1,rlow,u10e,v10e,wstar,xlt02,xlt024,xlt10 &
513 & ,xlt104,xlu10,xlu104,xu10,xu104,zt02,zt10,ztat02,ztat10 &
514 & ,ztau,ztau10,zu10,zuuz
550 z0=max(ustfc*ustar*ustar,1.59e-5)
567 zu=fzu1*sqrt(sqrt(z0*ustar*rvisc))/ustar
570 uz0=(ulow*rwgh+uz0)*0.5
571 vz0=(vlow*rwgh+vz0)*0.5
579 thz0=((wghtt*thlow+ths)/(wghtt+1.)+thz0)*0.5
580 qz0=((wghtq*qlow+qs)/(wghtq+1.)+qz0)*0.5
582 thz0=(wghtt*thlow+ths)/(wghtt+1.)
583 qz0=(wghtq*qlow+qs)/(wghtq+1.)
588 IF(ustar>=ustr.AND.ustar<ustc)
THEN
593 zt=fzt2*sqrt(sqrt(z0*ustar*rvisc))/ustar
599 thz0=((wghtt*thlow+ths)/(wghtt+1.)+thz0)*0.5
600 qz0=((wghtq*qlow+qs)/(wghtq+1.)+qz0)*0.5
602 thz0=(wghtt*thlow+ths)/(wghtt+1.)
603 qz0=(wghtq*qlow+qs)/(wghtq+1.)
623 thm=(thelow+thz0)*0.5
626 b=(elocp/tem-1.-p608)*thm
628 dthv=((thelow-thz0)*((qlow+qz0+cwmlow)*(0.5*p608)+1.) &
629 & +(qlow-qz0+cwmlow)*a+cwmlow*b)
631 du2=max((ulow-uz0)**2+(vlow-vz0)**2,epsu2)
632 rib=btgx*dthv*zsl/du2
654 rlmo=elfc*akhs*dthv/ustar**3
661 zetalu=min(max(zetalu,ztmin1),ztmax1)
662 zetalt=min(max(zetalt,ztmin1),ztmax1)
663 zetau=min(max(zetau,ztmin1/rzsu),ztmax1/rzsu)
664 zetat=min(max(zetat,ztmin1/rzst),ztmax1/rzst)
670 rz=(zetau-ztmin1)/dzeta1
675 psmz=(psim1(k+2)-psim1(k+1))*rdzt+psim1(k+1)
677 rz=(zetalu-ztmin1)/dzeta1
682 psmzl=(psim1(k+2)-psim1(k+1))*rdzt+psim1(k+1)
684 simm=psmzl-psmz+rlogu
686 rz=(zetat-ztmin1)/dzeta1
691 pshz=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
693 rz=(zetalt-ztmin1)/dzeta1
698 pshzl=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
700 simh=(pshzl-pshz+rlogt)*fh01
703 if(abs(simm)<1.e-10.or.abs(simh)<1.e-10)
then
704 write(0,*)
' simm=',simm,
' simh=',simh,
' at i=',i,
' j=',j
708 if(abs(simm).lt.1.e-10.or.abs(simh).lt.1.e-10)
then
709 print*,
'SIMM,PSMZL,PSMZ,RLOGU=',simm,psmzl,psmz,rlogu
710 print*,
'SIMH,PSHZL,PSHZ,RLOGT,FH01=',simh,pshzl,pshz,rlogt,fh01
711 print*,
'USTARK,CXCHS=',ustark,cxchs
712 print*,
'PSIM1(1,2),K=',psim1(k+1),psim1(k+2),k
713 print*,
'ZETAU,ZTMIN1,DZETA1=',zetau,ztmin1,dzeta1
714 print*,
'PSIH1(1,2),RDZT=',psih1(k+1),psih1(k+2),rdzt
715 print*,
'ZSLU,ZSLT,RLMO,ZU,ZT=',zslu,zslt,rlmo,zu,zt
716 print*,
'A,B,DTHV,DU2,RIB=',a,b,dthv,du2,rib
718 errmsg =
'ERROR(SFCDIF): '
724 akms=max(ustark/simm,cxchs)
725 akhs=max(ustark/simh,cxchs)
732 wstar2=wwst2*abs(btgh*akhs*dthv)**(2./3.)
736 ustar=max(sqrt(akms*sqrt(du2+wstar2)),epsust)
768 thm=(thelow+thz0)*0.5
771 b=(elocp/tem-1.-p608)*thm
773 dthv=((thelow-thz0)*((qlow+qz0+cwmlow)*(0.5*p608)+1.) &
774 & +(qlow-qz0+cwmlow)*a+cwmlow*b)
776 du2=max((ulow-uz0)**2+(vlow-vz0)**2,epsu2)
777 rib=btgx*dthv*zsl/du2
808 zzil=zilfc*(1.0+(rib/ric)*(rib/ric)*czetmax)
810 zzil=zilfc*(czetmax+1.0)
819 land_point_iteration:
DO itr=1,itrmx
825 zt=max(exp(zzil*sqrt(ustar*zu))*zu,epszt)
834 rlmo=elfc*akhs*dthv/ustar**3
840 zetalu=min(max(zetalu,ztmin2),ztmax2)
841 zetalt=min(max(zetalt,ztmin2),ztmax2)
842 zetau=min(max(zetau,ztmin2/rzsu),ztmax2/rzsu)
843 zetat=min(max(zetat,ztmin2/rzst),ztmax2/rzst)
849 rz=(zetau-ztmin2)/dzeta2
854 psmz=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
856 rz=(zetalu-ztmin2)/dzeta2
861 psmzl=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
863 simm=psmzl-psmz+rlogu
879 rz=(zetat-ztmin2)/dzeta2
884 pshz=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
886 rz=(zetalt-ztmin2)/dzeta2
891 pshzl=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
893 simh=(pshzl-pshz+rlogt)*fh02
896 akms=max(ustark/simm,cxchl)
897 akhs=max(ustark/simh,cxchl)
904 wstar2=wwst2*abs(btgh*akhs*dthv)**(2./3.)
909 ustar=max(sqrt(akms*sqrt(du2+wstar2)),epsust)
912 ENDDO land_point_iteration
948 wstar=sqrt(wstar2)/wwst
950 umflx=akms*(ulow -uz0 )
951 vmflx=akms*(vlow -vz0 )
952 hsflx=akhs*(thlow-thz0)
953 hlflx=akhs*(qlow -qz0 )
988 ztau10=min(max(ztau10,ztmin1),ztmax1)
989 ztat02=min(max(ztat02,ztmin1),ztmax1)
990 ztat10=min(max(ztat10,ztmin1),ztmax1)
992 rz=(ztau10-ztmin1)/dzeta1
997 psm10=(psim1(k+2)-psim1(k+1))*rdzt+psim1(k+1)
999 simm10=psm10-psmz+rlnu10
1001 rz=(ztat02-ztmin1)/dzeta1
1006 psh02=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
1008 simh02=(psh02-pshz+rlnt02)*fh01
1010 rz=(ztat10-ztmin1)/dzeta1
1015 psh10=(psih1(k+2)-psih1(k+1))*rdzt+psih1(k+1)
1017 simh10=(psh10-pshz+rlnt10)*fh01
1019 akms10=max(ustark/simm10,cxchs)
1020 akhs02=max(ustark/simh02,cxchs)
1021 akhs10=max(ustark/simh10,cxchs)
1030 ztau10=min(max(ztau10,ztmin2),ztmax2)
1031 ztat02=min(max(ztat02,ztmin2),ztmax2)
1032 ztat10=min(max(ztat10,ztmin2),ztmax2)
1034 rz=(ztau10-ztmin2)/dzeta2
1039 psm10=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
1041 simm10=psm10-psmz+rlnu10
1043 rz=(ztat02-ztmin2)/dzeta2
1048 psh02=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
1050 simh02=(psh02-pshz+rlnt02)*fh02
1052 rz=(ztat10-ztmin2)/dzeta2
1057 psh10=(psih2(k+2)-psih2(k+1))*rdzt+psih2(k+1)
1059 simh10=(psh10-pshz+rlnt10)*fh02
1061 akms10=ustark/simm10
1062 akhs02=ustark/simh02
1063 akhs10=ustark/simh10
1065 IF(akms10<=cxchl) akms10=akms
1066 IF(akhs02<=cxchl) akhs02=akhs
1067 IF(akhs10<=cxchl) akhs10=akhs
1075 u10 =umflx/akms10+uz0
1076 v10 =vmflx/akms10+vz0
1077 th02=hsflx/akhs02+thz0
1078 th10=hsflx/akhs10+thz0
1079 q02 =hlflx/akhs02+qz0
1080 q10 =hlflx/akhs10+qz0
1081 term1=-0.068283/tlow
1082 pshltr=psfc*exp(term1)
1098 zuuz=amin1(zu*0.50,0.18)
1099 zu=amax1(zu*0.35,zuuz)
1108 ztau10=min(max(ztau10,ztmin2),ztmax2)
1109 zetau=min(max(zetau,ztmin2/rzsu),ztmax2/rzsu)
1111 rz=(ztau10-ztmin2)/dzeta2
1116 psm10=(psim2(k+2)-psim2(k+1))*rdzt+psim2(k+1)
1117 simm10=psm10-psmz+rlnu10
1118 ekms10=max(ustark/simm10,cxchl)
1120 u10e=umflx/ekms10+uz0
1121 v10e=vmflx/ekms10+vz0
1132 rlow=plow/(r_d*tlow)
1137 qfx=-rlow*hlflx*wetm
1142 qgh=((1.-seamask)*pq0+seamask*pq0sea) &
1143 & /plow*exp(a2s*(tlow-a3s)/(tlow-a4s))
1144 cpm=cp*(1.+0.8*qlow)
1152 qs=qlow+qfx/(rlow*akhs)
1156 fm10=simm10+wght*simm
1157 fh2=simh02+wghtt*simh