[Fortran] 纯文本查看 复制代码
subroutine calcvis3(q,clw,rnw,ice,snow,t,prs,ista_num,vis3)
implicit none
integer :: ista_num
integer :: k
real,dimension(ista_num) :: t,q,clw,ice,rnw,snow,prs,vis3
real,dimension(ista_num) :: vovermd,celkel,tice,coeflc,coeflp,coeffc,coeffp,exponlc, &
exponlp,exponfp,const1,rhoice,rhowat,tv,rhoair,beta,qprc, &
qcld,qrain,qsnow,qclw,qclice,rgas,exponfc,concfc,concfp, &
conclp,conclc
CHARACTER METH*1
----------------------------------------
celkel=273.15 !K
rgas=287.04 !J/K/kg
tice=celkel-10.
coeflc=144.7
coeflp=2.24
coeffc=327.8
coeffp=10.36
exponlc=0.8800
exponlp=0.7500
exponfc=1.0000
exponfp=0.7776
const1=-LOG(.02)
rhoice=917.
rhowat=1000.
DO k=1,ista_num
qprc(k)=rnw(k)+snow(k)
qcld(k)=clw(k)+ice(k)
qrain(k)=rnw(k)
qsnow(k)=snow(k)
qclw(k)=clw(k)
qclice(k)=ice(k)
rhoair(k)=prs(k)/(287.04*t(k)*(1.+0.608*q(k))) !air density 空气通量密度
vovermd(k)=(1.+q(k)/1000.)/rhoair(k)+(qclw(k)+qrain(k))&
/rhowat+(qclice(k)+qsnow(k))/rhoice
conclc(k)=qclw(k)/vovermd(k)*1000.
conclp(k)=qrain(k)/vovermd(k)*1000.
concfc(k)=qclice(k)/vovermd(k)*1000.
concfp(k)=qsnow(k)/vovermd(k)*1000.
beta=coeffc*concfc(k)**exponfc+coeffp*concfp(k)**exponfp &
+coeflc*conclc(k)**exponlc+coeflp*conclp(k)**exponlp &
+1.E-10
vis3(k)=min(25.,const1/beta)
ENDDO
RETURN
END