FRZLVL.f

Go to the documentation of this file.

!-----------------------------------------------------------------------------
      SUBROUTINE frzlvl(ZFRZ,RHFRZ,PFRZL)

!     
!     
      use vrbls3d, only: pint, t, zmid, q, pmid
      use vrbls2d, only: fis, tshltr, pshltr, qshltr
      use masks, only: lmh
      use params_mod, only: gi, d00, capa, d0065, tfrz, pq0, a2, a3, a4
      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
!
!     DECLARE VARIABLES.
!     
      REAL,dimension(ista:iend,jsta:jend) :: rhfrz, zfrz, pfrzl
      integer i,j,llmh,l
      real htsfc,psfc,tsfc,qsfc,qsat,rhsfc,delz,delt,delq,delalp,     &
           delzp,zl,dzabv,qfrz,alpl,alph,alpfrz,pfrz,qsfrz,rhz,zu,    &
           dzfr,es
!     
!*********************************************************************
!     START FRZLVL.
!
!
!     
!     LOOP OVER HORIZONTAL GRID.
!     
!!$omp  parallel do                                                   &
!    & private(i,j,alpfrz,alph,alpl,delalp,delq,delt,delz,            &
!    &         delzp,dzabv,dzfr,htsfc,l,llmh,psfc,qfrz,               &
!    &         qsat,qsfc,qsfrz,rhsfc,rhz,tsfc,                        &
!    &         zl,zu)

       DO 20 j=jsta,jend
       DO 20 i=ista,iend
         htsfc    = fis(i,j)*gi
         llmh     = nint(lmh(i,j))
         rhfrz(i,j) = d00
         zfrz(i,j)  = htsfc
         psfc    = pint(i,j,llmh+1)
         pfrzl(i,j) = psfc
!     
!        CHECK IF FREEZING LEVEL IS AT THE GROUND.
!     
!         IF(SM(I,J)/=SPVAL .AND. THZ0(I,J)/=SPVAL .AND.        &
!          THS(I,J)/=SPVAL)THEN
!          TSFC = (SM(I,J)*THZ0(I,J)+(1.-SM(I,J))*THS(I,J))     &
!           *(PINT(I,J,NINT(LMH(I,J))+1)/P1000)**CAPA
! Per AWC's request, use 2m T instead of skin T so that freezing level
! would be above ground more often
         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  
         IF (tsfc<=tfrz) THEN
!            ZFRZ(I,J) = HTSFC+(TSFC-TFRZ)/D0065
        zfrz(i,j) = htsfc+2.0+(tsfc-tfrz)/d0065
!       IF(SM(I,J)/=SPVAL .AND. QZ0(I,J)/=SPVAL .AND.      &
!             QS(I,J)/=SPVAL)THEN
!             QSFC    = SM(I,J)*QZ0(I,J)+(1.-SM(I,J))*QS(I,J)
! GFS does not output QS           
!            ELSE IF(QSHLTR(I,J)/=SPVAL)THEN
        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  
!
            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
            pfrzl(i,j)= psfc
            cycle 
         ENDIF
!     
!        OTHERWISE, LOCATE THE FREEZING LEVEL ALOFT.
!
         DO 10 l = llmh,1,-1
            IF (t(i,j,l)<=tfrz) THEN
               IF (l<llmh) THEN
                  delz = zmid(i,j,l)-zmid(i,j,l+1)
                  zl   = zmid(i,j,l+1)
                  delt = t(i,j,l)-t(i,j,l+1)
                  zfrz(i,j) = zl + (tfrz-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   = alog(pmid(i,j,l+1))
                  alph   = alog(pmid(i,j,l))
                  alpfrz = alpl + (alph-alpl)/delz*dzabv
                  pfrz   = exp(alpfrz)
                  pfrzl(i,j)  = pfrz
          IF(modelname == 'GFS' .OR.modelname == 'RAPR')THEN
                es=fpvsnew(tfrz)
                es=min(es,pfrz)
                qsfrz=con_eps*es/(pfrz+con_epsm1*es)
              ELSE 
                    qsfrz=pq0/pfrz    &
                     *exp(a2*(tfrz-a3)/(tfrz-a4))
                  END IF
!     
                  rhz      = qfrz/qsfrz
                  rhz      = amax1(0.01,rhz)
                  rhz      = amin1(rhz,1.0)
                  rhfrz(i,j) = rhz
!     
               ELSE
                  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 + (tfrz-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    = alog(pmid(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(tfrz)
                es=min(es,pfrz)
                qsfrz=con_eps*es/(pfrz+con_epsm1*es)
              ELSE 
                    qsfrz=pq0/pfrz   &
                     *exp(a2*(tfrz-a3)/(tfrz-a4))
                  END IF
!
                  rhz     = qfrz/qsfrz
                  rhz     = amax1(0.01,rhz)
                  rhz     = amin1(rhz,1.0)
                  rhfrz(i,j)= rhz
               ENDIF
!     
!              BOUND FREEZING LEVEL RH.  FREEZING 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))
               EXIT             
            ENDIF
 10      CONTINUE
20   CONTINUE
!     
!     END OF ROUTINE.
!     
      RETURN
      END