FRZLVL2.f

Go to the documentation of this file.

!-------------------------------------------------------------------------------

      SUBROUTINE frzlvl2(ISOTHERM,ZFRZ,RHFRZ,PFRZL)

!     
      use vrbls3d, only: pint, t, zmid, pmid, q, zint, alpint
      use vrbls2d, only: fis, tshltr, pshltr, qz0, qs, qshltr
      use masks, only: lmh, sm
      use params_mod, only: gi, d00, capa, d0065, tfrz, pq0, a2, a3, a4, d50
      use ctlblk_mod, only: jsta, jend, spval, lm, modelname, im, ista, iend
      use physcons_post, only: con_rd, con_rv, con_eps, con_epsm1
      use upp_physics, only: fpvsnew

      implicit none

!- - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!      implicit none
!     
!     DECLARE VARIABLES.
!
      REAL,PARAMETER::pucap=300.0e2
      real,intent(in)                   ::  isotherm
      REAL,dimension(ista:iend,jsta:jend),intent(out) ::  rhfrz, zfrz, pfrzl
!jw
      integer i,j,l,lice,llmh
      real htsfc,psfc,qsfc,rhsfc,qw,qsat,delz,delt,delq,delalp,delzp,  &
           zl,zu,dzabv,qfrz,alpl,alph,alpfrz,pfrz,qsfrz,rhz,dzfr, &
           tsfc,es
!     
!*********************************************************************
!     START FRZLVL.
!
!     LOOP OVER HORIZONTAL GRID.
!     

      DO 20 j=jsta,jend
      DO 20 i=ista,iend
         IF(fis(i,j)<spval)THEN
         htsfc    = fis(i,j)*gi
         llmh     = nint(lmh(i,j))
         rhfrz(i,j) = d00
         zfrz(i,j)  = htsfc
         psfc     = pint(i,j,llmh)
         pfrzl(i,j) = psfc
!     
!        FIND THE HIGHEST LAYER WHERE THE TEMPERATURE
!        CHANGES FROM ABOVE TO BELOW ISOTHERM.
!     
!         TSFC = (SM(I,J)*THZ0(I,J)+(1.-SM(I,J))*THS(I,J))    &
!            *(PINT(I,J,NINT(LMH(I,J))+1)/P1000)**CAPA   
         IF(tshltr(i,j)/=spval .AND. pshltr(i,j)/=spval)THEN
          tsfc=tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa
         ELSE
! GFS analysis does not have flux file to retrieve TSFC from
          tsfc=t(i,j,lm)+d0065*(zmid(i,j,lm)-htsfc-2.0)
         END IF
         lice=llmh
! Per AWC's request, put a 300 mb cap for highest isothermal level so that it
! does not go into stratosphere
         DO l = llmh-1,1,-1
            IF (pmid(i,j,l)>=pucap .AND. &
        (t(i,j,l)<=isotherm.AND.t(i,j,l+1)>isotherm))lice=l
         ENDDO
!     
!        CHECK IF ISOTHERM LEVEL IS AT THE GROUND.
!     
         IF (lice==llmh.AND.tsfc<=isotherm) THEN
            zfrz(i,j) = htsfc+2.0+(tsfc-isotherm)/d0065
            qsfc    = sm(i,j)*qz0(i,j)+(1.-sm(i,j))*qs(i,j)
            IF(qshltr(i,j)/=spval)THEN
             psfc=pshltr(i,j)
             qsfc=qshltr(i,j)
            ELSE
             qsfc=q(i,j,lm)
             psfc=pmid(i,j,lm)
            END IF
            pfrzl(i,j) = psfc
!
            IF(modelname == 'GFS' .OR. modelname == 'RAPR')THEN
             es=fpvsnew(tsfc)
             es=min(es,psfc)
             qsat=con_eps*es/(psfc+con_epsm1*es)
            ELSE
             qsat=pq0/psfc  &
               *exp(a2*(tsfc-a3)/(tsfc-a4))
            END IF
!
            rhsfc   = qsfc/qsat
            rhsfc   = amax1(0.01,rhsfc)
            rhsfc   = amin1(rhsfc,1.0)
            rhfrz(i,j)= rhsfc
!     
!        OTHERWISE, LOCATE THE ISOTHERM LEVEL ALOFT.
!
         ELSE IF (lice<llmh) THEN
                  l=lice
                  delz = d50*(zint(i,j,l)-zint(i,j,l+2))
                  zl   = d50*(zint(i,j,l+1)+zint(i,j,l+2))
                  delt = t(i,j,l)-t(i,j,l+1)
                  zfrz(i,j) = zl+(isotherm-t(i,j,l+1))/delt*delz
!     
                  dzabv = zfrz(i,j)-zl
                  delq  = q(i,j,l)-q(i,j,l+1)
                  qfrz  = q(i,j,l+1) + delq/delz*dzabv
                  qfrz  = amax1(0.0,qfrz)
!     
                  alpl   = alpint(i,j,l+2)
                  alph   = alpint(i,j,l)
                  delalp = alph - alpl
                  delzp  = zint(i,j,l)-zint(i,j,l+2)
                  dzfr   = zfrz(i,j) - zint(i,j,l+2)
                  alpfrz = alpl + delalp/delzp*dzfr
                  pfrz   = exp(alpfrz)
                  pfrzl(i,j) = pfrz
                  IF(modelname == 'GFS'.OR.modelname == 'RAPR')THEN
                    es=fpvsnew(isotherm)
                    es=min(es,pfrz)
                    qsfrz=con_eps*es/(pfrz+con_epsm1*es)
                  ELSE
                    qsfrz=pq0/pfrz  & 
                     *exp(a2*(isotherm-a3)/(isotherm-a4))
                  END IF
!                  QSFRZ  = PQ0/PFRZ
!     
                  rhz      = qfrz/qsfrz
                  rhz      = amax1(0.01,rhz)
                  rhz      = amin1(rhz,1.0)
                  rhfrz(i,j) = rhz
!     
         ELSE
                  l=lice
                  zu      = zmid(i,j,l)
                  zl      = htsfc+2.0
                  delz    = zu-zl
                  IF(tshltr(i,j)/=spval .AND. pshltr(i,j)/=spval)THEN
                   tsfc=tshltr(i,j)*(pshltr(i,j)*1.e-5)**capa
                  ELSE
! GFS analysis does not have flux file to retrieve TSFC from
                   tsfc=t(i,j,lm)+d0065*(zmid(i,j,lm)-htsfc-2.0)
                  END IF
                  delt    = t(i,j,l)-tsfc
                  zfrz(i,j) = zl + (isotherm-tsfc)/delt*delz
!     
                  dzabv   = zfrz(i,j)-zl
! GFS does not output QS
                  IF(qshltr(i,j)/=spval)THEN
                   qsfc=qshltr(i,j)
                  ELSE
                   qsfc=q(i,j,lm)
                  END IF
                  delq    = q(i,j,l)-qsfc
                  qfrz    = qsfc + delq/delz*dzabv
                  qfrz    = amax1(0.0,qfrz)
!     
                  alph    = alpint(i,j,l)
                  alpl    = alog(psfc)
                  delalp  = alph-alpl
                  alpfrz  = alpl + delalp/delz*dzabv
                  pfrz    = exp(alpfrz)
                  pfrzl(i,j) = pfrz
                  IF(modelname == 'GFS'.OR.modelname == 'RAPR')THEN
                    es=fpvsnew(isotherm)
                    es=min(es,pfrz)
                    qsfrz=con_eps*es/(pfrz+con_epsm1*es)
                  ELSE
                    qsfrz=pq0/pfrz  &
                     *exp(a2*(isotherm-a3)/(isotherm-a4))
                  END IF
!
                  rhz     = qfrz/qsfrz
                  rhz     = amax1(0.01,rhz)
                  rhz     = amin1(rhz,1.0)
                  rhfrz(i,j)= rhz
         ENDIF
!     
!              BOUND ISOTHERM LEVEL RH.  ISOTHERM LEVEL HEIGHT IS
!              MEASURED WITH RESPECT TO MEAN SEA LEVEL.
!
               rhfrz(i,j) = amax1(0.01,rhfrz(i,j))
               rhfrz(i,j) = amin1(rhfrz(i,j),1.00)
               zfrz(i,j)  = amax1(0.0,zfrz(i,j))
         ELSE
               rhfrz(i,j) = spval
               zfrz(i,j)  = spval
         ENDIF
 20   CONTINUE
!     
!     END OF ROUTINE.
!     
      RETURN
      END