CALMICT.f

Go to the documentation of this file.

      SUBROUTINE calmict_new(P1D,T1D,Q1D,C1D,FI1D,FR1D,FS1D,CUREFL,     &
                        qw1,qi1,qr1,qs1,dbz1,dbzr1,dbzi1,dbzc1,nlice1,nrain1)

!
      use params_mod, only: dbzmin, epsq, tfrz, eps, rd, d608
      use ctlblk_mod, only: jsta, jend, jsta_2l,jend_2u,im
      use cmassi_mod, only: t_ice, rqr_drmin, n0rmin, cn0r_dmrmin,      &
              mdrmin, rqr_drmax, cn0r_dmrmax, mdrmax, n0r0, xmrmin,     &
              xmrmax, massi, cn0r0, mdimin, xmimax, mdimax 
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
      INTEGER indexs, indexr
!aligo
      REAL, PARAMETER :: cice=1.634e13, cwet=1./.189, cboth=cice/.224,   &
     &  NLI_min=1.E3, RFmax=45.259, RQmix=0.1E-3,NSI_max=250.E3
!aligo
      real,dimension(IM,jsta_2l:jend_2u),intent(in)    :: p1d,t1d,q1d,c1d,fi1d,fr1d, &
                                                          fs1d,curefl
      real,dimension(IM,jsta_2l:jend_2u),intent(inout) :: qw1,qi1,qr1,qs1,dbz1,dbzr1,&
                                                          dbzi1,dbzc1,nlice1,nrain1
      
      integer i,j
      real :: tc,frain,fice,rimef,xsimass,qice,qsat,esat,wv,rho,rrho,    &
     &        RQR,DRmm,Qsigrd,WVQW,Dum,XLi,Qlice,WC,DLI,NLImax,NSImax,   &
     &        RQLICE, N0r,Ztot,Zrain,Zice,Zconv,Zmin,Zmix,NLICE,NSmICE,  &
     &        QSmICE,NRAIN,NMIX,Zsmice
      logical :: large_rf, hail
      real,external :: fpvs
!************************************************************************
!--- Determine composition of condensate in the form of cloud water, 
!    total ice (cloud ice & snow), & rain following GSMDRIVE in NMM model
!
      zmin=10.**(0.1*dbzmin)
      DO j=jsta,jend
        DO i=1,im
          qw1(i,j)=0.
          qi1(i,j)=0.
          qr1(i,j)=0.
          qs1(i,j)=0.
          nlice1(i,j)=0.
          nrain1(i,j)=0.
          dbz1(i,j)=dbzmin
          dbzr1(i,j)=dbzmin
          dbzi1(i,j)=dbzmin
          dbzc1(i,j)=dbzmin
        ENDDO
      ENDDO
      DO j=jsta,jend
        DO i=1,im
          ztot=0.             !--- Total radar reflectivity
          zrain=0.            !--- Radar reflectivity from rain
          zice=0.             !--- Radar reflectivity from ice
          zsmice=0.           !--- Radar reflectivity from small ice
          zconv=curefl(i,j)   !--- Radar reflectivity from convection
          IF (c1d(i,j) <= epsq) THEN
!
!--- Skip rest of calculatiions if no condensate is present
!
            go to 10
          ELSE
            wc=c1d(i,j)
          ENDIF
!
!--- Code below is from GSMDRIVE for determining:
!    QI1 - total ice (cloud ice & snow) mixing ratio
!    QW1 - cloud water mixing ratio
!    QR1 - rain mixing ratio
!
          tc=t1d(i,j)-tfrz
          fice=fi1d(i,j)
          frain=fr1d(i,j)
          IF (tc<=t_ice .OR. fice>=1.) THEN
!          IF (Fice>=1.) THEN  
            qi1(i,j)=wc
          ELSE IF (fice <= 0.) THEN
            qw1(i,j)=wc
          ELSE
            qi1(i,j)=fice*wc
            qw1(i,j)=wc-qi1(i,j)
          ENDIF   
          IF (qw1(i,j)>0. .AND. frain>0.) THEN
            IF (frain >= 1.) THEN 
              qr1(i,j)=qw1(i,j)
              qw1(i,j)=0.
            ELSE
              qr1(i,j)=frain*qw1(i,j)
              qw1(i,j)=qw1(i,j)-qr1(i,j)
            ENDIF 
          ENDIF
          wv=q1d(i,j)/(1.-q1d(i,j))
!
!--- Saturation vapor pressure w/r/t water ( >=0C ) or ice ( <0C )
!
          esat=1000.*fpvs(t1d(i,j))
          qsat=eps*esat/(p1d(i,j)-esat)
          rho=p1d(i,j)/(rd*t1d(i,j)*(1.+d608*q1d(i,j)))
          rrho=1./rho
  !
  !--- Based on code from GSMCOLUMN in model to determine reflectivity from rain
  !
          rqr=0.
          IF (qr1(i,j) > epsq) THEN
            rqr=rho*qr1(i,j)
            IF (rqr <= rqr_drmin) THEN
              n0r=max(n0rmin, cn0r_dmrmin*rqr)
              indexr=mdrmin
            ELSE IF (rqr >= rqr_drmax) THEN
              n0r=cn0r_dmrmax*rqr
              indexr=mdrmax
            ELSE
              n0r=n0r0
              indexr=max( xmrmin, min(cn0r0*rqr**.25, xmrmax) )
            ENDIF
  !
  !--- INDEXR is the mean drop size in microns; convert to mm
  !
            drmm=1.e-3*REAL(indexr)
  !
  !--- Number concentration of rain drops (convert INDEXR to m)
  !
            nrain=n0r*1.e-6*REAL(indexr)
            nrain1(i,j)=nrain
            zrain=0.72*n0r*drmm*drmm*drmm*drmm*drmm*drmm*drmm
          ENDIF        !--- End IF (QR1(I,J) > EPSQ) block
!
!--- Based on code from GSMCOLUMN in model to determine partition of 
!    total ice into cloud ice & snow (precipitation ice)
!
          rqlice=0.
          IF (qi1(i,j) > epsq) THEN
            qice=qi1(i,j)
!
!  ->  Small ice particles are assumed to have a mean diameter of 50 microns.
!  * INDEXS  - mean size of snow to the nearest micron (units of microns)
!  * RimeF   - Rime Factor, which is the mass ratio of total (unrimed &
!              rimed) ice mass to the unrimed ice mass (>=1)
!  * QTICE   - time-averaged mixing ratio of total ice
!  * QLICE   - mixing ratio of large ice
!  * RQLICE  - mass content of large ice
!  * NLICE1  - time-averaged number concentration of large ice
!
            IF (tc>=0.) THEN
   !
   !--- Eliminate small ice particle contributions for melting & sublimation
   !
              nsmice=0.
              qsmice=0.
            ELSE
!
!--- Max # conc of small ice crystals based on 10% of total ice content
!
!              NSImax=0.1*RHO*QICE/MASSI(MDImin)
!aligo
               nsimax=max(nsi_max,0.1*rho*qice/massi(mdimin) )
!aligo
!
!-- Specify Fletcher, Cooper, Meyers, etc. here for ice nuclei concentrations
!
              nsmice=min(0.01*exp(-0.6*tc), nsimax)       !- Fletcher (1962)
              dum=rrho*massi(mdimin)
              nsmice=min(nsmice, qice/dum)
              qsmice=nsmice*dum
            ENDIF            ! End IF (TC>=0.) THEN
            qlice=max(0., qice-qsmice)
            qs1(i,j)=qlice
            qi1(i,j)=qsmice
            rimef=amax1(1., fs1d(i,j) )
            rimef=min(rimef, rfmax)
            rqlice=rho*qlice
            dum=xmimax*exp(.0536*tc)
            indexs=min(mdimax, max(mdimin, int(dum) ) )
!
!-- Specify NLImax depending on presence of high density ice (rime factors >10)
!
            IF (rimef>10.) THEN
              large_rf=.true.         !-- For convective precipitation
              nlimax=1.e3
            ELSE
              large_rf=.false.        !-- For non-convective precipitation
              dum=max(tc, t_ice)
              nlimax=10.e3*exp(-0.017*dum)
            ENDIF
            nlice=rqlice/(rimef*massi(indexs))
            dum=nli_min*massi(mdimin)     !-- Minimum large ice mixing ratio
new_nlice:  IF (rqlice<dum) THEN
              nlice=rqlice/massi(mdimin)
            ELSE IF (nlice<nli_min .OR. nlice>nlimax) THEN  new_nlice
!
!--- Force NLICE to be between NLI_min and NLImax, but allow for exceptions
!
              hail=.false.
              nlice=max(nli_min, min(nlimax, nlice) )
              xli=rqlice/(nlice*rimef)
new_size:     IF (xli <= massi(mdimin) ) THEN
                indexs=mdimin
              ELSE IF (xli <= massi(450) ) THEN   new_size
                dli=9.5885e5*xli**.42066         ! DLI in microns
                indexs=min(mdimax, max(mdimin, int(dli) ) )
              ELSE IF (xli <= massi(mdimax) ) THEN   new_size
                dli=3.9751e6*xli**.49870         ! DLI in microns
                indexs=min(mdimax, max(mdimin, int(dli) ) )
              ELSE   new_size
                indexs=mdimax
                IF (large_rf)  hail=.true.
              ENDIF    new_size
no_hail:      IF (.NOT. hail) THEN
                nlice=rqlice/(rimef*massi(indexs))   !-- NLICE > NLImax
              ENDIF    no_hail
            ENDIF      new_nlice
            nlice1(i,j)=nlice
   !
   !--- Equation (C.8) in Ferrier (1994, JAS, p. 272), which when
   !    converted from cgs units to mks units results in the same
   !    value for Cice, which is equal to the {} term below:
   !
   !    Zi={.224*720*(10**18)/[(PI*RHOL)**2]}*(RHO*QLICE)**2/NLICE1(I,J),
   !    where RHOL=1000 kg/m**3 is the density of liquid water
   !
   !--- Valid only for exponential ice distributions
   !
            IF (nsmice > 0.) THEN
               zsmice=cice*rho*rho*qsmice*qsmice/nsmice
            ENDIF
            if (nlice1(i,j) /= 0.0) zice=cice*rqlice*rqlice/nlice1(i,j)
            IF (tc>=0.) zice=cwet*zice      ! increased for wet ice
          ENDIF                 ! End IF (QI1(I,J) > 0.) THEN
!
!--- Assumed enhanced radar reflectivity when rain and ice coexist
!    above an assumed threshold mass content, RQmix
!
dbz_mix:  IF (rqr>rqmix .AND. rqlice>rqmix) THEN
            IF (rqr>rqlice) THEN
              nmix=nrain
            ELSE
              nmix=nlice1(i,j)
            ENDIF
            dum=rqr+rqlice
            zmix=cboth*dum*dum/nmix
            IF (zmix > zrain+zice) THEN
              IF (rqr>rqlice) THEN
                zrain=zmix-zice
              ELSE
                zice=zmix-zrain
              ENDIF
            ENDIF
          ENDIF  dbz_mix
!
!---  Calculate total (convective + grid-scale) radar reflectivity
!
10        zice=zice+zsmice
          ztot=zrain+zice+zconv
          IF (ztot > zmin)  dbz1(i,j)= 10.*alog10(ztot)
          IF (zrain > zmin) dbzr1(i,j)=10.*alog10(zrain)
          IF (zice > zmin)  dbzi1(i,j)=10.*alog10(zice)
!          IF (Zconv > Zmin) DBZC1(I,J)=10.*ALOG10(Zsmice)
          IF (zconv > zmin) dbzc1(i,j)=10.*alog10(zconv)
        ENDDO
      ENDDO
!
      RETURN
      END
!
!-- For the old version of the microphysics:
!
      SUBROUTINE calmict_old(P1D,T1D,Q1D,C1D,FI1D,FR1D,FS1D,CUREFL,     &
                        qw1,qi1,qr1,qs1,dbz1,dbzr1,dbzi1,dbzc1,nlice1,nrain1)
      use params_mod, only: dbzmin, epsq, tfrz, eps, rd, d608, oneps, nlimin
      use ctlblk_mod, only: jsta, jend, jsta_2l, jend_2u, im
      use rhgrd_mod, only: rhgrd
      use cmassi_mod, only: t_ice, rqr_drmin, n0rmin, cn0r_dmrmin, mdrmin, &
              rqr_drmax,cn0r_dmrmax, mdrmax, n0r0, xmrmin, xmrmax,flarge2, &
              massi, cn0r0, mdimin, xmimax, mdimax,nlimax 
!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      implicit none
!
      INTEGER indexs, indexr
      REAL, PARAMETER :: cice=1.634e13

      real,dimension(IM,jsta_2l:jend_2u),intent(in)    :: p1d,t1d,q1d,c1d,fi1d,fr1d, &
                                                          fs1d,curefl
      real,dimension(IM,jsta_2l:jend_2u),intent(inout) :: qw1,qi1,qr1,qs1,dbz1,dbzr1,&
                                                          dbzi1,dbzc1,nlice1,nrain1

      REAL n0r,ztot,zrain,zice,zconv,zmin
      integer i,j
      real tc, frain,fice,flimass,flarge,          &
           fsmall,rimef,xsimass,qice,qsat,esat,wv,rho,rrho,rqr,          &
           drmm,qsigrd,wvqw,dum,xli,qlice,wc,dli,xlimass
      real,external :: fpvs
!************************************************************************
!--- Determine composition of condensate in the form of cloud water,
!    total ice (cloud ice & snow), & rain following GSMDRIVE in NMM model
!
      zmin=10.**(0.1*dbzmin)
      DO j=jsta,jend
        DO i=1,im
          qw1(i,j)=0.
          qi1(i,j)=0.
          qr1(i,j)=0.
          qs1(i,j)=0.
          nlice1(i,j)=0.
          nrain1(i,j)=0.
          dbz1(i,j)=dbzmin
          dbzr1(i,j)=dbzmin
          dbzi1(i,j)=dbzmin
          dbzc1(i,j)=dbzmin
        ENDDO
      ENDDO
      DO j=jsta,jend
        DO i=1,im
          zrain=0.            !--- Radar reflectivity from rain
          zice=0.             !--- Radar reflectivity from ice
          zconv=curefl(i,j)   !--- Radar reflectivity from convection
          IF (c1d(i,j) <= epsq) THEN
!
!--- Skip rest of calculatiions if no condensate is present
!
            go to 10
          ELSE
            wc=c1d(i,j)
          ENDIF
!
!--- Code below is from GSMDRIVE for determining:
!    QI1 - total ice (cloud ice & snow) mixing ratio
!    QW1 - cloud water mixing ratio
!    QR1 - rain mixing ratio
!
          tc=t1d(i,j)-tfrz
          fice=fi1d(i,j)
          frain=fr1d(i,j)
          IF (tc<=t_ice .OR. fice>=1.) THEN
!          IF (Fice>=1.) THEN
            qi1(i,j)=wc
          ELSE IF (fice <= 0.) THEN
            qw1(i,j)=wc
          ELSE
            qi1(i,j)=fice*wc
            qw1(i,j)=wc-qi1(i,j)
          ENDIF
          IF (qw1(i,j)>0. .AND. frain>0.) THEN
            IF (frain >= 1.) THEN
              qr1(i,j)=qw1(i,j)
              qw1(i,j)=0.
            ELSE
              qr1(i,j)=frain*qw1(i,j)
              qw1(i,j)=qw1(i,j)-qr1(i,j)
            ENDIF
          ENDIF
          wv=q1d(i,j)/(1.-q1d(i,j))
!
!--- Saturation vapor pressure w/r/t water ( >=0C ) or ice ( <0C )
!
          esat=1000.*fpvs(t1d(i,j))
          qsat=eps*esat/(p1d(i,j)-esat)
          rho=p1d(i,j)/(rd*t1d(i,j)*(1.+d608*q1d(i,j)))
          rrho=1./rho
  !
  !--- Based on code from GSMCOLUMN in model to determine reflectivity from rain
  !
          IF (qr1(i,j) > epsq) THEN
            rqr=rho*qr1(i,j)
            IF (rqr <= rqr_drmin) THEN
              n0r=max(n0rmin, cn0r_dmrmin*rqr)
              indexr=mdrmin
            ELSE IF (rqr >= rqr_drmax) THEN
              n0r=cn0r_dmrmax*rqr
              indexr=mdrmax
            ELSE
              n0r=n0r0
              indexr=max( xmrmin, min(cn0r0*rqr**.25, xmrmax) )
            ENDIF
  !
  !--- INDEXR is the mean drop size in microns; convert to mm
  !
            drmm=1.e-3*REAL(indexr)
            zrain=0.72*n0r*drmm*drmm*drmm*drmm*drmm*drmm*drmm
  !
  !--- Number concentration of rain drops (convert INDEXR to m)
  !
            nrain1(i,j)=n0r*1.e-6*REAL(indexr)
          ENDIF        !--- End IF (QR1(I,J) > EPSQ) block
!
!--- Based on code from GSMCOLUMN in model to determine partition of
!    total ice into cloud ice & snow (precipitation ice)
!
          IF (qi1(i,j) > epsq) THEN
            qice=qi1(i,j)
            rho=p1d(i,j)/(rd*t1d(i,j)*(1.+oneps*q1d(i,j)))
            rrho=1./rho
            qsigrd=rhgrd*qsat
            wvqw=wv+qw1(i,j)
!
! * FLARGE  - ratio of number of large ice to total (large & small) ice
! * FSMALL  - ratio of number of small ice crystals to large ice particles
!  ->  Small ice particles are assumed to have a mean diameter of 50 microns.
!  * XSIMASS - used for calculating small ice mixing ratio
!  * XLIMASS - used for calculating large ice mixing ratio
!  * INDEXS  - mean size of snow to the nearest micron (units of microns)
!  * RimeF   - Rime Factor, which is the mass ratio of total (unrimed &
!              rimed) ice mass to the unrimed ice mass (>=1)
!  * FLIMASS - mass fraction of large ice
!  * QTICE   - time-averaged mixing ratio of total ice
!  * QLICE   - time-averaged mixing ratio of large ice
!  * NLICE1   - time-averaged number concentration of large ice
!
            IF (tc>=0. .OR. wvqw<qsigrd) THEN
              flarge=1.
            ELSE
              flarge=flarge2    !-- specified in MICROINIT.f
!!              IF (TC>=-8. .AND. TC<=-3.) FLARGE=.5*FLARGE
            ENDIF
            fsmall=(1.-flarge)/flarge
            xsimass=rrho*massi(mdimin)*fsmall
            dum=xmimax*exp(.0536*tc)
            indexs=min(mdimax, max(mdimin, int(dum) ) )
            rimef=amax1(1., fs1d(i,j) )
            xlimass=rrho*rimef*massi(indexs)
            flimass=xlimass/(xlimass+xsimass)
            qlice=flimass*qice
            nlice1(i,j)=qlice/xlimass
            IF (nlice1(i,j)<nlimin .OR. nlice1(i,j)>nlimax) THEN
!
!--- Force NLICE1 to be between NLImin and NLImax
!
              dum=max(nlimin, min(nlimax, nlice1(i,j)) )
              xli=rho*(qice/dum-xsimass)/rimef
              IF (xli <= massi(mdimin) ) THEN
                indexs=mdimin
              ELSE IF (xli <= massi(450) ) THEN
                dli=9.5885e5*xli**.42066         ! DLI in microns
                indexs=min(mdimax, max(mdimin, int(dli) ) )
              ELSE IF (xli <= massi(mdimax) ) THEN
                dli=3.9751e6*xli**.49870         ! DLI in microns
                indexs=min(mdimax, max(mdimin, int(dli) ) )
              ELSE
                indexs=mdimax
!
!--- 8/22/01: Increase density of large ice if maximum limits
!    are reached for number concentration (NLImax) and mean size
!    (MDImax).  Done to increase fall out of ice.
!
                IF (dum >= nlimax)                             &
                  rimef=rho*(qice/nlimax-xsimass)/massi(indexs)
              ENDIF             ! End IF (XLI <= MASSI(MDImin) )
              xlimass=rrho*rimef*massi(indexs)
              flimass=xlimass/(xlimass+xsimass)
              qlice=flimass*qice
              nlice1(i,j)=qlice/xlimass
            ENDIF               ! End IF (NLICE<NLImin ...
            qs1(i,j)=amin1(qi1(i,j), qlice)
            qi1(i,j)=amax1(0., qi1(i,j)-qs1(i,j))
   !
   !--- Equation (C.8) in Ferrier (1994, JAS, p. 272), which when
   !    converted from cgs units to mks units results in the same
   !    value for Cice, which is equal to the {} term below:
   !
   !    Zi={.224*720*(10**18)/[(PI*RHOL)**2]}*(RHO*QLICE)**2/NLICE1(I,J),
   !    where RHOL=1000 kg/m**3 is the density of liquid water
   !
   !--- Valid only for exponential ice distributions
   !
            zice=cice*rho*rho*qlice*qlice/nlice1(i,j)
          ENDIF                 ! End IF (QI1(I,J) > 0.) THEN
!
!---  Calculate total (convective + grid-scale) radar reflectivity
10        ztot=zrain+zice+zconv
          IF (ztot > zmin)  dbz1(i,j)= 10.*alog10(ztot)
          IF (zrain > zmin) dbzr1(i,j)=10.*alog10(zrain)
          IF (zice > zmin)  dbzi1(i,j)=10.*alog10(zice)
          IF (zconv > zmin) dbzc1(i,j)=10.*alog10(zconv)
        ENDDO
      ENDDO
!
      RETURN
      END