MDL2SIGMA2.f

Go to the documentation of this file.

!
      SUBROUTINE mdl2sigma2

!
!
      use vrbls3d, only:  pint, pmid, t, zint, q
!      use vrbls2d, only:
      use masks, only: lmh
      use params_mod, only: pq0, a2, a3, a4, rgamog
      use ctlblk_mod, only: pt, jsta_2l, jend_2u, spval, lp1, lm, jsta, jend,&
                            grib, cfld, datapd, fld_info, im, jm, im_jm
      use rqstfld_mod, only: iget, lvls, id, iavblfld, lvlsxml
!     
      implicit none
!
      integer,PARAMETER :: lsig=5
!     
!     DECLARE VARIABLES.
!     
      LOGICAL readthk
!     REAL,dimension(im,jm) :: FSL, TSL, QSL, osl, usl, vsl, q2sl, fsl1,     &
      REAL,dimension(im,jsta_2l:jend_2u) :: tsl
      REAL,dimension(im,jsta_2l:jend_2u) :: grid1
      REAL sigo(lsig+1),dsigo(lsig),asigo(lsig)
!
!     INTEGER,dimension(im,jm) :: IHOLD,JHOLD,NL1X,NL1XF
      INTEGER,dimension(im,jsta_2l:jend_2u) :: nl1x
!
!
!--- Definition of the following 2D (horizontal) dummy variables
!
!  C1D   - total condensate
!  QW1   - cloud water mixing ratio
!  QI1   - cloud ice mixing ratio
!  QR1   - rain mixing ratio
!  QS1   - snow mixing ratio
!
!     REAL,dimension(im,jm) :: C1D,QW1,QI1,QR1,qs1,qg1,akh 
!
      integer i,j,l,ll,lp,llmh,nhold,ii,jj
      real ptsigo,psigo,apsigo,fact,ai,bi,tmt0,qsat,tvrl,  &
           tvrblo,tblo,ql,rhl,zl,pl,tl
!
!     
!******************************************************************************
!
!     START MDL2P. 
!     
!     SET TOTAL NUMBER OF POINTS ON OUTPUT GRID.
!
!---------------------------------------------------------------
!
!     *** PART I ***
!
!     VERTICAL INTERPOLATION OF EVERYTHING ELSE.  EXECUTE ONLY
!     IF THERE'S SOMETHING WE WANT.
!
      IF((iget(296)>0) ) THEN  !!Air Quality (Plee Oct2003)
!
!---------------------------------------------------------------------
!
!---  VERTICAL INTERPOLATION OF GEOPOTENTIAL, SPECIFIC HUMIDITY, TEMPERATURE, 
!     OMEGA, TKE, & CLOUD FIELDS.  START AT THE UPPERMOST TARGET SIGMA LEVEL.
!
        ptsigo=pt   
        readthk=.false.
        IF(readthk)THEN   ! EITHER READ DSG THICKNESS
          READ(41)dsigo  !DSIGO FROM TOP TO BOTTOM
!
          sigo(1)=0.0
          DO l=2,lsig+1
            sigo(l)=sigo(l-1)+dsigo(lsig-l+2)
          END DO 
          sigo(lsig+1)=1.0
          DO l=1,lsig
            asigo(l)=0.5*(sigo(l)+sigo(l+1))
          END DO
        ELSE  ! SPECIFY SIGO
          asigo( 1)=   0.7000
          asigo( 2)=   0.7500
          asigo( 3)=   0.8000
          asigo( 4)=   0.8500
          asigo( 5)=   0.9000
        END IF
!***
!***  BECAUSE SIGMA LAYERS DO NOT GO UNDERGROUND,  DO ALL
!***  INTERPOLATION ABOVE GROUND NOW.
!***
!

        DO 310 lp=1,lsig
          nhold=0
!
          DO j=jsta_2l,jend_2u
            DO i=1,im

!
              tsl(i,j)=spval
!
!***  LOCATE VERTICAL INDEX OF MODEL MIDLAYER JUST BELOW
!***  THE PRESSURE LEVEL TO WHICH WE ARE INTERPOLATING.
!
              nl1x(i,j)=lp1
              DO l=2,lm
                llmh = nint(lmh(i,j))
                psigo=ptsigo+asigo(lp)*(pint(i,j,llmh+1)-ptsigo)
                IF(nl1x(i,j)==lp1.AND.pmid(i,j,l)>psigo)THEN
                  nl1x(i,j)=l
                ENDIF
              ENDDO
!
!  IF THE PRESSURE LEVEL IS BELOW THE LOWEST MODEL MIDLAYER
!  BUT STILL ABOVE THE LOWEST MODEL BOTTOM INTERFACE,
!  WE WILL NOT CONSIDER IT UNDERGROUND AND THE INTERPOLATION
!  WILL EXTRAPOLATE TO THAT POINT
!
              IF(nl1x(i,j)==lp1.AND.pint(i,j,llmh+1)>=psigo)THEN
                nl1x(i,j)=lm
              ENDIF
!
!        if(NL1X(I,J)==LP1)print*,'Debug: NL1X=LP1 AT '
!     1 ,i,j,lp
            ENDDO
          ENDDO
!
!mptest        IF(NHOLD==0)GO TO 310
!
!!$omp  parallel do
!!$omp& private(nn,i,j,ll,fact,qsat,rhl)
!hc        DO 220 NN=1,NHOLD
!hc        I=IHOLD(NN)
!hc        J=JHOLD(NN)
!         DO 220 J=JSTA,JEND
!         DO 220 J=JSTA_2L,JEND_2U
          DO 220 j=jsta,jend           ! Moorthi on Nov 26, 2014
            DO 220 i=1,im
              ll=nl1x(i,j)
!---------------------------------------------------------------------
!***  VERTICAL INTERPOLATION OF GEOPOTENTIAL, TEMPERATURE, SPECIFIC
!***  HUMIDITY, CLOUD WATER/ICE, OMEGA, WINDS, AND TKE.
!---------------------------------------------------------------------
!
!HC           IF(NL1X(I,J)<=LM)THEN
              llmh = nint(lmh(i,j))
              psigo=ptsigo+asigo(lp)*(pint(i,j,llmh+1)-ptsigo) 
              apsigo=log(psigo)
              IF(nl1x(i,j)<=llmh.and. &
               (pmid(i,j,ll)-pmid(i,j,ll-1))/=0.)THEN
!
!---------------------------------------------------------------------
!          INTERPOLATE LINEARLY IN LOG(P)
!***  EXTRAPOLATE ABOVE THE TOPMOST MIDLAYER OF THE MODEL
!***  INTERPOLATION BETWEEN NORMAL LOWER AND UPPER BOUNDS
!***  EXTRAPOLATE BELOW LOWEST MODEL MIDLAYER (BUT STILL ABOVE GROUND)
!---------------------------------------------------------------------
!

                fact=(apsigo-log(pmid(i,j,ll)))/                        &
     &              (log(pmid(i,j,ll))-log(pmid(i,j,ll-1)))
                tsl(i,j)=t(i,j,ll)+(t(i,j,ll)-t(i,j,ll-1))*fact
! FOR UNDERGROUND PRESSURE LEVELS, ASSUME TEMPERATURE TO CHANGE 
! ADIABATICLY, RH TO BE THE SAME AS THE AVERAGE OF THE 2ND AND 3RD
! LAYERS FROM THE GOUND, WIND TO BE THE SAME AS THE LOWEST LEVEL ABOVE
! GOUND
              ELSE
                ii=91
                jj=13
!                if(i==ii.and.j==jj)                                 &
!                  print*,'Debug: underg extra at i,j,lp',i,j,lp
                pl = pint(i,j,lm-1)
                zl = zint(i,j,lm-1)
                tl = 0.5*(t(i,j,lm-2)+t(i,j,lm-1))
                ql = 0.5*(q(i,j,lm-2)+q(i,j,lm-1))
                tmt0  = tl - a3
                ai    = 0.008855
                bi    = 1.
                IF(tmt0<-20.)THEN
                  ai = 0.007225
                  bi = 0.9674
                ENDIF
                qsat = pq0/pl*exp(a2*(tl-a3)/(tl-a4))
!
                rhl = ql/qsat
!
                IF(rhl>1.)THEN
                  rhl = 1.
                  ql  = rhl*qsat
                ENDIF
!
                IF(rhl<0.01)THEN
                  rhl = 0.01
                  ql  = rhl*qsat
                ENDIF
!
!               print *,' tl=',tl,' ql=',ql,' i=',i,' j=',j,' pl=',pl
                tvrl   = tl*(1.+0.608*ql)
                tvrblo = tvrl*(psigo/pl)**rgamog
                tblo   = tvrblo/(1.+0.608*ql)
!     
                tmt0  = tblo-a3
                ai    = 0.008855
                bi    = 1.
                IF(tmt0<-20.)THEN
                  ai = 0.007225
                  bi = 0.9674
                ENDIF
                qsat = pq0/psigo*exp(a2*(tblo-a3)/(tblo-a4))
!
                tsl(i,j) = tblo
              END IF
  220       CONTINUE

!---------------------------------------------------------------------
!        *** PART II ***
!---------------------------------------------------------------------
!---------------------------------------------------------------------
!
!        OUTPUT SELECTED FIELDS.
!     
!***  TEMPERATURE
!
        IF(iget(296)>0) THEN
          IF(lvls(lp,iget(296))>0)THEN
             DO j=jsta,jend
               DO i=1,im
                 grid1(i,j)=tsl(i,j)
               ENDDO
             ENDDO
            if(grib=='grib2')then
             cfld=cfld+1
             fld_info(cfld)%ifld=iavblfld(iget(296))
             fld_info(cfld)%lvl=lvlsxml(lp,iget(296))
             datapd(1:im,1:jend-jsta+1,cfld)=grid1(1:im,jsta:jend)
            endif
          ENDIF
        ENDIF
!     
!***  END OF MAIN VERTICAL LOOP
!     
  310   CONTINUE
!***  ENDIF FOR IF TEST SEEING IF WE WANT ANY OTHER VARIABLES
!
      ENDIF
!
!     END OF ROUTINE.
!
      RETURN
      END