[Fortran] syntaxhighlighter_viewsource syntaxhighlighter_copycode
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