OTLIFT.f

Go to the documentation of this file.

      SUBROUTINE otlift(SLINDX)

!     
      use vrbls3d,    only: pmid, t, q
      use vrbls2d,    only: t500
      use masks,      only: lmh
      use lookup_mod, only: thl, rdth, jtb, qs0, sqs, rdq,itb, ptbl, pl, &
                            rdp, the0, sthe, rdthe, ttbl
      use ctlblk_mod, only: jsta, jend, im, spval
      use params_mod, only: d00,h10e5, capa, elocp, eps, oneps
      use upp_physics, only: fpvsnew
!

!
      implicit none
!
!     SET LOCAL PARAMETERS.
       real,PARAMETER :: d8202=.820231e0 , h5e4=5.e4 , p500=50000.

!     
!     DECLARE VARIABLES.
      real,intent(out) :: slindx(im,jsta:jend)
      REAL :: tvp, esatp, qsatp
      REAL :: tth, tp, apesp, partmp, thesp, tpsp
      REAL :: bqs00, sqs00, bqs10, sqs10, bq, sq, tq
      REAL :: p00, p10, p01, p11, t00, t10, t01, t11
      REAL :: bthe00, sthe00, bthe10, sthe10, bth, sth
      REAL :: tqq, qq, qbt, tthbt, tbt, apebt, ppq, pp
!
      INTEGER :: i, j, lbtm, ittbk, iq, it, iptbk, ith, ip, iqtb
      INTEGER :: ittb, iptb, ithtb
!     
!***********************************************************************
!     START OTLIFT HERE
!     
!     INTIALIZE LIFTED INDEX ARRAY TO ZERO.
!$omp parallel do private(i,j)
      DO j=jsta,jend
        DO i=1,im
          slindx(i,j) = d00
        ENDDO
      ENDDO
!--------------FIND EXNER AT LOWEST LEVEL-------------------------------
      DO j=jsta,jend
        DO i=1,im
          lbtm=nint(lmh(i,j))
          IF(t(i,j,lbtm)<spval .AND. q(i,j,lbtm)<spval) THEN
          tbt = t(i,j,lbtm)
          qbt = q(i,j,lbtm)
          apebt = (h10e5/pmid(i,j,lbtm))**capa
!--------------SCALING POTENTIAL TEMPERATURE & TABLE INDEX--------------
          tthbt = tbt*apebt
          tth = (tthbt-thl)*rdth
          tqq = tth-aint(tth)
          ittb = int(tth)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
          IF(ittb < 1)THEN
            ittb = 1
            tqq = d00
          ENDIF
            IF(ittb >= jtb)THEN
            ittb = jtb-1
            tqq = d00
          ENDIF
!--------------BASE AND SCALING FACTOR FOR SPEC. HUMIDITY---------------
          ittbk = ittb
          bqs00=qs0(ittbk)
          sqs00=sqs(ittbk)
          bqs10=qs0(ittbk+1)
          sqs10=sqs(ittbk+1)
!--------------SCALING SPEC. HUMIDITY & TABLE INDEX---------------------
          bq=(bqs10-bqs00)*tqq+bqs00
          sq=(sqs10-sqs00)*tqq+sqs00
          tq=(qbt-bq)/sq*rdq
          ppq = tq-aint(tq)
          iqtb = int(tq)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
          IF(iqtb < 1)THEN
            iqtb = 1
            ppq = d00
          ENDIF
          IF(iqtb >= itb)THEN
            iqtb = itb-1
            ppq = d00
          ENDIF
!--------------SATURATION PRESSURE AT FOUR SURROUNDING TABLE PTS.-------
          iq=iqtb
          it = ittb
          p00=ptbl(iq,it)
          p10=ptbl(iq+1,it)
          p01=ptbl(iq,it+1)
          p11=ptbl(iq+1,it+1)
!--------------SATURATION POINT VARIABLES AT THE BOTTOM-----------------
          tpsp = p00+(p10-p00)*ppq+(p01-p00)*tqq       &
                +(p00-p10-p01+p11)*ppq*tqq
          IF(tpsp <= d00) tpsp = h10e5
          apesp = (h10e5/tpsp)**capa
          thesp = tthbt*exp(elocp*qbt*apesp/tthbt)
!--------------SCALING PRESSURE & TT TABLE INDEX------------------------
          tp=(h5e4-pl)*rdp
          qq = tp-aint(tp)
          iptb = int(tp)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
          IF(iptb < 1)THEN
            iptb = 1
            qq = d00
          ENDIF
          IF(iptb >= itb)THEN
            iptb = itb-1
            qq = d00
          ENDIF
!--------------BASE AND SCALING FACTOR FOR THE--------------------------
          iptbk=iptb
          bthe00=the0(iptbk)
          sthe00=sthe(iptbk)
          bthe10=the0(iptbk+1)
          sthe10=sthe(iptbk+1)
!--------------SCALING THE & TT TABLE INDEX-----------------------------
          bth=(bthe10-bthe00)*qq+bthe00
          sth=(sthe10-sthe00)*qq+sthe00
          tth=(thesp-bth)/sth*rdthe
          pp = tth-aint(tth)
          ithtb = int(tth)+1
!--------------KEEPING INDICES WITHIN THE TABLE-------------------------
          IF(ithtb < 1)THEN
            ithtb = 1
            pp = d00
          ENDIF
          IF(ithtb >= jtb)THEN
            ithtb = jtb-1
            pp = d00
          ENDIF
!--------------TEMPERATURE AT FOUR SURROUNDING TT TABLE PTS.------------
          ith=ithtb
          ip=iptb
          t00=ttbl(ith,ip)
          t10=ttbl(ith+1,ip)
          t01=ttbl(ith,ip+1)
          t11=ttbl(ith+1,ip+1)
!--------------PARCEL TEMPERATURE AT 500MB----------------------------
          IF(tpsp >= h5e4)THEN
            partmp=(t00+(t10-t00)*pp+(t01-t00)*qq    &
                +(t00-t10-t01+t11)*pp*qq)
          ELSE
            partmp=tbt*apebt*d8202
          ENDIF
!--------------LIFTED INDEX---------------------------------------------
!
! GSM  THE PARCEL TEMPERATURE AT 500 MB HAS BEEN COMPUTED, AND WE
!       FIND THE MIXING RATIO AT THAT LEVEL WHICH WILL BE THE SATURATION
!       VALUE SINCE WE'RE FOLLOWING A MOIST ADIABAT.    NOTE THAT THE
!       AMBIENT 500 MB SHOULD PROBABLY BE VIRTUALIZED, BUT THE IMPACT
!       OF MOISTURE AT THAT LEVEL IS QUITE SMALL

         esatp=fpvsnew(partmp)
         qsatp=eps*esatp/(p500-esatp*oneps)
         tvp=partmp*(1+0.608*qsatp)
         slindx(i,j)=t500(i,j)-tvp
        ENDIF !end T(I,J,LBTM)<SPVAL 
        ENDDO
      ENDDO
!       write(*,*) ' in otlift t500 partmp ',t500(1,1),partmp(1,1)
!       write(*,*) ' in otlift tbt ',tbt(1,1)
!
      RETURN
      END