text:
c*****************************************************************************
c*****************************************************************************
c
c
c                    ROTUNNO-EMANUEL 1987 MODEL
c
c   Rotunno, R., and K.A. Emanuel, 1987:  An air-sea interaction theory for
c      tropical cyclones, Part II: Evolutionary study using axisymmetric
c      nonhydrostatic numerical model.  J. Atmos. Sci., 44, 542-561.
c
c   Modified where indicated by K. Emanuel; last modified July, 2020
c
c
c*****************************************************************************
c*****************************************************************************
c*****************************************************************************
c
      PROGRAM HURAXI
C
      PARAMETER(M=400,N=80,MP1=M+1,MM1=M-1,NP1=N+1,NM1=N-1)
      PARAMETER (NTM=5000)
C
      DIMENSION  QRAD(M,N), UB(N),
     $ U(MP1,NP1),V(MP1,NP1),W(MP1,NP1),P(MP1,NP1),
     $ T(MP1,NP1)  , QV(MP1,NP1) , QL(MP1,NP1),VTERM2D(M,N),
     $ DUM(M,N), RC2(N)  , ZC2(N), CPTDR(N), CPTDZ(N), QS(N),
     $ RHOTVB(N), RHOTVW(NP1), RHOT(N),RHOW(NP1),RTB(N),RQVB(N),
     $ DTSR(M),DTSZ(M,N),DTSRW(N),DTSZW(N),DTLR(M),DTLZ(N),DTSV(M),
     $ A(N), B(N), C(N), D(M,N), E(N), WS(M,NP1), PS(M,N)
     $,UAV(MP1,N),VAV(M,N),WAV(M,NP1),TAV(M,N),PAV(M,N),QVAV(M,N)
     $,QLAV(M,N), XKMAV(M,N), SST(M), ESS(M)
       REAL UMN(NTM),VMX(NTM),VMXS(NTM),WMX(NTM),
     $ RMX(NTM),PMN(NTM), TTP(NTM)
       REAL VMAXZ(NTM,N),THEHOV(NTM,M)
       REAL VBHOV(NTM,M),UBHOV(NTM,M),VTHOV(NTM,M),UTHOV(NTM,M)
       REAL RPL(M),UPL(M),VPL(M),PPL(M),THEPL(M)
       REAL ZPL(N), DUMT(M)
       character*18 chtime
       character*2 anp
       CHARACTER*1 PDEP,EVAP
C
      COMMON/DIFF/
     $ UA(MP1,N) , WA(M,NP1), VA(M,N)  , TA(M,N), QVA(M,N), QLA(M,N),
     $ U1(MP1,N) , V1(M,N)  , W1(M,NP1), T1(M,N), QV1(M,N), QL1(M,N),
     $ XKM(M,NP1), P1(M,N)  , R(MP1) ,RS(M),Z(NP1),ZS(N),TAD(M,N),
     $ RTBW(N),RQVBW(N),TAU(N),TB(N), QVB(N),PN(NP1), PD(NP1),
     $ TSURF(M),QSURF(M),RDR,RDZ,RDR2,RDZ2,DTS,DTL,DZ,
     $ XVL2, XHL2, CD, CE, XLDCP, XKAPPA, A1, G,
     $ CERS(M), CDRS(M), DISS, VCAP, IVMAXS
C
        DATA  CSTAR/30./,EP/0.1/,NZTR/1/,QVBT/0./
C
C   Open Parameter Input File and Read Parameters
C
        OPEN(UNIT=11,FILE='hurr.in',STATUS='old')
C
        READ(11,*)
        READ(11,*)
        READ(11,*)PDS
        READ(11,*)FC
        READ(11,*)CD
        READ(11,*)CE
        READ(11,*)CD1
        READ(11,*)CDCAP
        READ(11,*)CKCAP
        READ(11,*)VCAP
        READ(11,*)PDEP
        READ(11,*)XHL
        READ(11,*)XVL
        READ(11,*)DISS
        READ(11,*)TAUR
        READ(11,*)RADMAX
        READ(11,*)ACT
        READ(11,*)VTERM
        READ(11,*)VTERMSNOW
        READ(11,*)EVAP
        READ(11,*)TMIN
        READ(11,*)
        READ(11,*)
        READ(11,*)
        READ(11,*)RMAX
        READ(11,*)RO
        READ(11,*)VMAX
        READ(11,*)TMID
        READ(11,*)RSST
        READ(11,*)
        READ(11,*)
        READ(11,*)
        READ(11,*)DT
        READ(11,*)NS
        READ(11,*)EPS
        READ(11,*)ALPHA
        READ(11,*)RB
        READ(11,*)ZB
        READ(11,*)NSPONGE
        READ(11,*)ETIME
        READ(11,*)TAVE
        READ(11,*)PLTIME
        READ(11,*)TIMMAX
        READ(11,*)TIMEPL
        READ(11,*)ROG
        READ(11,*)ZOG
c
        CLOSE(11)
c
c-----------------------------------------------------------------------
c
c   Read input sounding file
c
        DZTEMP=ZB*1000./FLOAT(N)
        CALL INTERPOLATE(N,PDS,DZTEMP,TB,QVB,DISS,TBS,VMAXTH)
c
c  Modify VMAX according to Emanuel and Rotunno, J. Atmos. Sci., 2011
c
        CR1=CE/CD
        IF(CR1.EQ.2.0)CR1=2.01 ! Avoid singularity (singular limit)
        VMAXTH=VMAXTH*SQRT(2.0*CR1)*(0.5*CR1)**(CR1/(4.-2.*CR1))
c
        TBSW=TBS
        TBS=(TBS+273.)*(1000./PDS)**0.286
        TBT=2.0*TB(N)-TB(N-1)
c
c     open output file
c
        open(17,file='output/texout.txt',
     1   status='unknown',form='formatted')
c
c  Define and/or normalize parameters
c
        F=FC*1.0E-5
        CD=CD*0.001
        CD1=CD1*1.0E-5
        CE=CE*0.001
        CDCAP=CDCAP*0.001
        CKCAP=CKCAP*0.001
        RADT=1.0/(3600.*TAUR)
        RADMAX=2.*DT*RADMAX/(3600.*24.)
        RMAX=RMAX*1000.0
        RO=RO*1000.0
        RSST=RSST*1000.0
        MPL=(FLOAT(M)*ROG)/RB
        NPL=(FLOAT(N)*ZOG)/ZB
        RB=1000.*RB
        ZB=1000.*ZB
        ISTOP=ETIME*24.*3600./DT
        NZSP=N-NSPONGE+1
        IAVE=TAVE*3600.*24./DT
        IPLOT=PLTIME*3600.*24./DT
        ITMAX=TIMMAX*3600./DT
        IPRINT=ISTOP/2-1
        TIMEPL=TIMEPL*3600.*24.
        ISTART=1
        AVE=FLOAT(IAVE)
        IPL=ISTART-1+IPLOT
        PDEPA=1.0
        IF(PDEP.EQ.'N'.OR.PDEP.EQ.'n')PDEPA=0.0
        ACT=0.001*ACT ! Normalize autoconversion threshold (July, 2020)
c
c     Added  9/3/2004 by KE
c
c	iave=min(iave,ipl)
c
      DR  = RB / FLOAT(M)
      DZ  = ZB / FLOAT(N)
      DTS = DT / FLOAT(NS)
      DTL = DT
      DZ2 = DZ * DZ
      DR2 = DR * DR
      RDR = 1. / DR
      RDZ = 1. / DZ
      RDR2= 1. / DR2
      RDZ2= 1. / DZ2
      PI  = 4. * ATAN(1.)
      C2  = 90000.
      CP  = 1005.
      CV  = 718.
      RD  = 287.
      XKAPPA= RD / CP
      XLDCP = 2500.
      G     = 9.81
      A1    = 7.5 * LOG(10.)
      PNS   = (PDS/1000.) ** XKAPPA
      TAMB  = TBS*PNS
      ESSS   =  6.11 * EXP( A1* (TAMB-273.)/( TAMB-36.)  )
      QSP   = .622 * ESSS /( PDS -ESSS )
      QVBS  = QVB(1)
      TVBS  = TBS * ( 1. + .61 * QVBS )
c      CC1 = .61 * G * DZ /(2.* CP * TVBS )/(1. + .61 *QVBS   )
c      CC2 = 1. - G * DZ / (2. * CP * TVBS * PNS )
c      CC5 = - TBS / PNS
      CC3   =  G * DZ / CP
c
c   Modified by K. Emanuel, 11/4/2010
c
c     XHL2  = .04 * DR2
      XHL2  = XHL*XHL
c
      XVL2  = XVL * XVL
      TIME  = DT * ( FLOAT(ISTART) - 1. )
C
C        DEFINE GRID ARRAYS
C
      DO 10 I = 1 , MP1
      R(I) = ( FLOAT(I) - 1.0 ) * DR
      IF( I .NE. MP1 )  RS(I) = ( FLOAT(I) - 0.5 ) * DR
   10 CONTINUE
      DO 20 J = 1 , NP1
      Z(J) = ( FLOAT(J) - 1.0 ) * DZ
      IF( J .NE. NP1)  ZS(J) = ( FLOAT(J) - 0.5 ) * DZ
   20 CONTINUE
      TSTOP  = DT * FLOAT(ISTOP)
      write(17,920)  TIME, TSTOP, ISTART, ISTOP, IPLOT, IPRINT,
     $ DT, EPS, XVL, F, RB, ZB, M, N, CD, CE
  920 FORMAT(1H ,
     $  9X,'TSTART=',F10.3,/,10X,'TSTOP =',F10.3,/,
     $ 10X,'ISTART=',I5   ,/,10X,'ISTOP =',I5   ,/,
     $ 10X,'IPLOT =',I5   ,/,10X,'IPRINT=',I5   ,/,
     $ 10X,'DT    =',F10.3,/,10X,'EPS   =',F10.3,/,
     $ 10X,'XL    =',F10.3,/,10X,'F     =',F7.5 ,/,
     $ 10X,'RB    =',F10.1,/,10X,'ZB    =',F10.3,/,
     $ 10X,'M     =',I3   ,/,10X,'N     =',I3   ,/,
     $ 10X,'CD    =',F10.5,/,10X,'CE    =',F10.5   )
C
C        BASE STATE (HYDROSTATIC EQUATION FOR PN)
C
      PN(1) = PNS -.5*CC3/(.5*(TB(1)+TBS)*(1.+.61*.5*(QVB(1)+QVBS)))
      PD(1) = 1000. * PN(1) ** ( 1. / XKAPPA )
      DO 40 J = 2 , N
       PN(J) = PN(J-1) - CC3/(.5*(TB(J)+TB(J-1))*
     $                          (1. + .61*.5*(QVB(J)+QVB(J-1)) ) )
         TBAR=0.5*(PN(J)*TB(J)+PN(J-1)*TB(J-1))
         IF(TBAR.LE.TMIN)THEN
          PN(J)=PN(J-1)*EXP(-CC3/TMIN)
          TB(J)=TMIN/PN(J)
         END IF
       PD(J) = 1000. * PN(J) ** ( 1. / XKAPPA )
   40 CONTINUE
C   mod by KE  8/17/04
      PN(NP1)=2.*PN(N)-PN(N-1)
        PD(NP1)=2.*PD(N)-PD(N-1)
C  end mod
        TBAR=0.5*(PN(N-1)*TB(N-1)+PN(N)*TB(N))
      PNT = PN(N) -.5*CC3/(.5*(TB(N)+TBT)*(1.+.61*.5*(QVB(N)+QVBT)))
        IF(TBAR.LE.TMIN)THEN
          PNT=PN(N)*EXP(-.5*CC3/TMIN)
        END IF
      PDT = 1000. * PNT ** ( 1. / XKAPPA )
      write(17,945)
  945 FORMAT(1H ,'     Z          P          PI        PT       ',
     $'QVB         QVS')
      write(17,940) ZB, PDT, PNT, TBT, QVBT, 0.
      DO 50 J = N , 2, -1
      ES =6.11*EXP(A1*(PN(J)*TB(J)- 273.)/(PN(J)*TB(J)-36.))
      QS(J) = .622 * ES /MAX(( PD(J)-ES ),ES)
      write(17,940) ZS(J), PD(J), PN(J), TB(J), QVB(J), QS(J)
  940 FORMAT(1H , 6(F10.4,1X) )
   50 CONTINUE
      write(17,940) 0., PDS, PNS, TBS, QVBS, QSP
C
C       ARRAYS FOR SMALL TIME STEP AND BUOYANCY CALCULATION
C
      DO 60 J = 1 , N
      RC2(J)    =  RDR*DTS*C2/(CP*TB(J)*(1.+.61*QVB(J)))
      CPTDR(J)  =  DTS*RDR * CP * TB(J)* (1. + .61 * QVB(J) )
      RTB(J)    = .5  * G / TB(J)
c      RQVB(J)   = .61 * G * .5 / ( 1. + .61 * QVB(J) )
      RHOTVB(J) = (1000./RD) * PN(J) ** (CV /RD)
      RHOT(J) = 100.* RHOTVB(J)/( TB(J)*( 1.+.61*QVB(J) ) )
      ZC2(J)    = RDZ*RC2(J)/RHOTVB(J)/RDR
      IF( J.EQ.1 )  GO TO 60
      PNW = .5 * ( PN(J) + PN(J-1) )
      TVW = .5*(TB(J)+TB(J-1) )*(1.+.305*(QVB(J)+QVB(J-1)))
      RHOTVW(J)=  (1000./RD) * PNW ** (CV /RD)
      RHOW(J)  = 100. * RHOTVW(J) / TVW
      RTBW(J) = G / ( TB(J) + TB(J-1) )
      RQVBW(J)= .305*G/( 1. + .305*(QVB(J)+QVB(J-1)) )
      CPTDZ(J)=DTS*RDZ*CP*.5*(TB(J)+TB(J-1))*
     1   (1.+.305*(QVB(J)+QVB(J-1)))
   60 CONTINUE
      RHOTVW(1  )= (1000./RD) * PNS ** (CV /RD)
      RHOTVW(NP1)= (1000./RD) * PNT ** (CV /RD)
      RHOW(1  ) = 100. * RHOTVW(1  ) / TVBS
      RHOW(NP1)=  100. * RHOTVW(NP1) / TBT
C
C         SPONGE LAYER DAMPING COEFFICIENT
C
      DO 65 J = 1 , N
      ZSP = ( ZS(J) - ZS(NZSP) ) / ( ZS(N) - ZS(NZSP) )
      IF(ZSP.LT.0.) TAU(J) = 0.
      IF(ZSP.GE.0..AND.ZSP.LE..5) TAU(J)=-.5*ALPHA*(1.-COS(ZSP*PI) )
      IF(ZSP.GT..5) TAU(J)=-.5*ALPHA*( 1. + PI*(ZSP-.5) )
   65 CONTINUE
C
C         ARRAYS TO INCREASE EFFICIENCY
C
      DTSF = .5 * DTS * F
      DTSG = .5 * DTS * G
      DTL2 = .5 * DTL
      DO 70 I = 1 , M
      DTLR(I)  = .5  * DTL * RDR /RS(I)
      DTSV(I)  = .5  * DTS  / RS(I)
      IF( I .NE. 1 )  DTSR(I)  = .25 * DTS * RDR / R(I)
   70 CONTINUE
      DO 75 J = 1 , N
      DTLZ(J)  = .5  * DTL * RDZ / RHOT(J)
      DTSZW(J) = .25 * DTS * RDZ / RHOW(J)
      DTSRW(J) = .25 * DTS * RDR / RHOW(J)
   75 CONTINUE
      DO 80 J = 1 , N
      DO 80 I = 2 , M
      DTSZ(I,J)  = .25 * DTS * RDZ /( RHOT(J) * R(I) )
   80 CONTINUE
C
C        ARRAYS FOR SEMI - IMPLICIT SMALL TIME STEP
C
      CC4=.25*(1.+EP)**2
      DO 85 J = 2 , N
      A(J)=   CC4*CPTDZ(J)*RHOTVW(J+1)*    ZC2(J)
      B(J)=1.+CC4*CPTDZ(J)*RHOTVW(J  )*(ZC2(J)+ZC2(J-1))
      C(J)=   CC4*CPTDZ(J)*RHOTVW(J-1)*   ZC2(J-1)
   85 CONTINUE
      E(1)=0.
      DO 95 I = 1 , M
      D(I,1)=0.
95    CONTINUE
      DO 90 J = 2 , N
      E(J) = A(J)/(B(J)-C(J)*E(J-1))
   90 CONTINUE
C
C       INITIAL CONDITIONS
C
C
C     TEMPERATURE AND VAPOR PRESSURE AT SEA SURFACE
C
      DO 5 I = 1 , M
      SST(I)=TAMB+ TMID/(1.+(RS(I)/RSST)**2)
      ESS(I)=  6.11 * EXP( A1* (SST(I)-273.)/(SST(I)-36.)  )
 5    CONTINUE
      zd=zs(nzsp)
      E1 = 3.
      D2 = 2. * RMAX / (    RO + RMAX )
      DO 100 J = 1 , N
      DO 100 I = 1 , M
      D1 = 2. * RMAX / ( RS(I) + RMAX )
      IF( RS(I) .LT. RO )THEN
       VR = SQRT(VMAX**2 *(RS(I)/RMAX)**2 *
     $ ( D1 ** E1 - D2 ** E1 )+.25*F*F*RS(I)**2) - .5 * F * RS(I)
      ELSE
       VR = 0.
      END IF
      V(I,J )   = VR * ZS(J)/ZS(NZTR)
      IF( ZS(J).GT.ZS(NZTR) ) V(I,J)=VR*(ZS(J)-ZD)/(ZS(NZTR)-ZD)
      IF( ZS(J).GT.ZS(NZSP) ) V(I,J) = 0.0
      V1(I,J)   = V(I,J)
      QL(I,J)   = 0.
      QL1(I,J)  = 0.
  100 CONTINUE
      DO 105 J = 1 , N
      P(M,J) = 0.
      P1(M,J)=P(M,J)
      DO 105 I = M , 2 , -1
      P(I-1,J) = P(I,J) - ( DTS / CPTDR(J) ) * .5 *
     $ ( V(I,J)*V(I,J)/RS(I) + V(I-1,J)*V(I-1,J)/RS(I-1)
     $   + F * ( V(I,J) + V(I-1,J) ))
      P1(I-1,J)= P(I-1,J)
  105 CONTINUE
      DO 106 I = 1 , M
      DO 107 J = 2 , NM1
      T(I,J) = TB(J) + .25*( CPTDZ(J+1)*(P(I,J+1)-P(I,J  ))
     $  +CPTDZ(J  )*(P(I,J  )-P(I,J-1)) )/(DTS*RTB(J))
      T1(I,J)=T(I,J)
  107 CONTINUE
      T(I,1) = TB(1) + .5*CPTDZ(2)*(P(I,2)-P(I,1  ))/(DTS*RTB(1))
      T(I,N) = TB(N) + .5*CPTDZ(N)*(P(I,N)-P(I,NM1))/(DTS*RTB(N))
      T1(I,1)=T(I,1)
      T1(I,N)=T(I,N)
  106 CONTINUE
C
C    GIVE 0'S TO THOSE ARRAY ELEMENTS THAT ARE OUTSIDE PHYSICAL BOUNDARIES (IE
C    TO ELEMENTS THAT WERE ADDED TO ARRAY DEFINITIONS TO AVOID INDEX OVERFLOW)
C
      DO 141 I=1,MP1
         V(I,NP1)=0.0
         P(I,NP1)=0.0
         T(I,NP1)=0.0
         QV(I,NP1)=0.0
         QL(I,NP1)=0.0
         U(I,NP1)=0.0
  141  CONTINUE
      W(MP1,NP1)=0.0
      DO 142 J=1,NP1
         V(MP1,J)=0.0
         P(MP1,J)=0.0
         T(MP1,J)=0.0
         QV(MP1,J)=0.0
         QL(MP1,J)=0.0
         W(MP1,J)=0.0
  142  CONTINUE
       U(MP1,NP1)=0.0
c
c
c    ADJUST Qv SO THAT THETAe= constant across vortex
c
c
      DO 108 J = 1 , N
      DO 108 I = 1 , M
c      QV(I,J)   = QVB(J)+PN(J)*(TB(J)-T(I,J))/XLDCP
      QV(I,J)   = QVB(J)
c       IF(R(I).LT.150000.0.AND.J.GT.1)THEN
c         PNST= PN(J) + P1(I,J) - P1(M,1)
c         PDST= 1000. * PNST ** (1./XKAPPA)
c         TEMP= PNST * T1(I,J)
c         ES  = 6.11 * EXP(A1 * ( TEMP  - 273. ) / ( TEMP  - 36. ) )
c         QV(I,J) = 0.98*.622 * ES /(PDST-ES)
c       END IF
      QV1(I,J)  = QV(I,J)
 108  CONTINUE
      DO 110 J = 1 , N
      DO 110 I = 1 , MP1
      U(I,J)  = 0.0
      U1(I,J) = U(I,J)
  110 CONTINUE
      DO 120 J = 1 , NP1
      DO 120 I = 1 , M
      W(I,J)  = 0.0
      W1(I,J) = W(I,J)
  120 CONTINUE
      DO 505 J = 1 , N
      UA(1,J) = 0.
  505 CONTINUE
        TIPL=0.0
        NPT=0
C
C   Write some arrays
C
        DO I=1,MPL
         RPL(I)=R(I)*0.001
        END DO
        DO 821 I=1,NPL
         ZPL(I)=0.001*DZ*FLOAT(I-1)
  821   CONTINUE
        OPEN(UNIT=22,FILE='output/rgraph.out',STATUS='UNKNOWN')
        OPEN(UNIT=23,FILE='output/zgraph.out',STATUS='UNKNOWN')
        WRITE(22,1025)(RPL(I),I=1,MPL)
 1025   FORMAT(2X,800(F12.7,' '))
        WRITE(23,1030)(ZPL(I),I=1,NPL)
 1030   FORMAT(2X,800(F8.4,1X))
        CLOSE(22)
        CLOSE(23)
        NPLOT=0
C
C        BEGIN TIME MARCH
C
      DO 500 ITIME = ISTART , ISTOP
      TIME = TIME + DT
      TIPL=TIPL + DT
C
C       SURFACE TEMP AND WATER-VAPOR MIXING RATIO
C
      DO 515 I = 1 , M
      PNSURF = PDEPA*(P1(I,1)-P1(M,1))/(1.-.5*G*DZ/(CP*SST(I)))
      PDST   = 1000. * (PNS +PNSURF) ** (1./XKAPPA)
      QSURF(I) = .622 * ESS(I) /(PDST -ESS(I))
      TSURF(I) = SST(I)/(PNS+PNSURF)
  515 CONTINUE
      DO 6 I = 1 , M
        CDRS(I)=CD+CD1
     $    *SQRT(.25*(U1(I,1)+U1(I+1,1))**2+V1(I,1)**2)
        CDRS(I)=MIN(CDRS(I),CDCAP)
        CERS(I)=CE+CD1
     $    *SQRT(.25*(U1(I,1)+U1(I+1,1))**2+V1(I,1)**2)
        CERS(I)=MIN(CERS(I),CKCAP)
C
    6 CONTINUE
C
      CALL DIFFUSE
C
C     RADIATION (R < 1 DEG/DAY)
C
      DO 150 J = 1 , N
      DO 150 I = 1 , M
      TDIF=T1(I,J)-TB(J)
        QRAD(I,J) = -DTL*TDIF*RADT
        QRAD(I,J)=MIN(QRAD(I,J),RADMAX)
        QRAD(I,J)=MAX(QRAD(I,J),-RADMAX)
  150 CONTINUE
C
C        FORCING FOR U EQUATION
C
      DO 510 J = 1 , N
      DO 510 I = 2 , M
      IF(J.EQ.1)THEN
        UMINUS=U(I,J)
      ELSE
        UMINUS=U(I,J-1)
      END IF
      UA(I,J)    = UA(I,J)
     1 - DTSR(I)   *(
     2 ( R(I+1)*U(I+1,J)+R(I  )*U(I  ,J) )*( U(I+1,J)-U(I  ,J) )
     3+( R(I  )*U(I  ,J)+R(I-1)*U(I-1,J) )*( U(I  ,J)-U(I-1,J) ) )
     4 - DTSZ(I,J) *(
     5 RHOW(J+1)*(RS(I)*W(I,J+1)+RS(I-1)*W(I-1,J+1))*(U(I,J+1)-U(I,J))
     6+RHOW(J  )*(RS(I)*W(I,J)+RS(I-1)*W(I-1,J))*(U(I,J)-UMINUS))
     6+ DTSV(I) * V(I,J)*V(I,J) +  DTSV(I-1) * V(I-1,J)*V(I-1,J)
     7  +  DTSF    * ( V(I,J) + V(I-1,J) )
  510 CONTINUE
C
C       OUTER BOUNDARY
C
      DO 525 J = 1 , N
      UA(MP1,J) =
     $ -MAX( U(MP1,J) + CSTAR , 0. )*DTL*RDR*(U1(MP1,J)-U1(M,J))
     $+ DTL * (V(M,J)*V(M,J)/RS(M) +  F * V(M,J)   )
  525 CONTINUE
C
C       DRAG AT J = 1
C
      UA(MP1,1) = UA(MP1,1)
     $ - CD *DTL*RDZ*U1(MP1,1)*SQRT( U1(MP1,1)**2+V1(M,1)**2 )
C
C       FORCING FOR W EQUATION
C
      DO 550 J = 2 , N
      DO  I = 2 , MM1
      WA(I,J)    = WA(I,J)
     1 -DTSRW(J)*(
     2 (RHOT(J  )*U(I+1,J)+RHOT(J-1)*U(I+1,J-1))*(W(I+1,J)-W(I  ,J))
     3+(RHOT(J  )*U(I  ,J)+RHOT(J-1)*U(I  ,J-1))*(W(I  ,J)-W(I-1,J)))
     4 -DTSZW(J)*(
     5 (RHOW(J+1)*W(I,J+1)+RHOW(J  )*W(I  ,J  ))*(W(I,J+1)-W(I,J  ))
     6+(RHOW(J-1)*W(I,J-1)+RHOW(J  )*W(I  ,J  ))*(W(I,J  )-W(I,J-1)))
     7   + DTS*RTBW(J)  * ( T(I,J) - TB(J)  + T(I,J-1)  - TB(J-1)  )
     8   + DTS*RQVBW(J) * (QV(I,J) - QVB(J) + QV(I,J-1) - QVB(J-1) )
     9   - DTSG * ( QL(I,J) +QL(I,J-1) )
      END DO
      I=1
      WA(I,J)    = WA(I,J)
     1 -DTSRW(J)*(
     2 (RHOT(J  )*U(I+1,J)+RHOT(J-1)*U(I+1,J-1))*(W(I+1,J)-W(I  ,J)))
     4 -DTSZW(J)*(
     5 (RHOW(J+1)*W(I,J+1)+RHOW(J  )*W(I  ,J  ))*(W(I,J+1)-W(I,J  ))
     6+(RHOW(J-1)*W(I,J-1)+RHOW(J  )*W(I  ,J  ))*(W(I,J  )-W(I,J-1)))
     7   + DTS*RTBW(J)  * ( T(I,J) - TB(J)  + T(I,J-1)  - TB(J-1)  )
     8   + DTS*RQVBW(J) * (QV(I,J) - QVB(J) + QV(I,J-1) - QVB(J-1) )
     9   - DTSG * ( QL(I,J) +QL(I,J-1) )
  550 CONTINUE
C
C       OUTER BOUNDARY
C
      DO 520 J = 2 , N
      UB(J) =      MAX( RHOT(J  ) * ( U(MP1,J) + U(M,J  ) )
     $                   +RHOT(J-1) * ( U(M,J-1) + U(M,J-1) )  , 0.)
  520 CONTINUE
      I = M
      DO 555 J = 2 , N
      WA(I,J)    = WA(I,J)
     1 -DTSRW(J) * UB(J) * ( W1(I  ,J) - W1(I-1,J) )
     4 -DTSZW(J)*(
     5 (RHOW(J+1)*W(I,J+1)+RHOW(J  )*W(I  ,J  ))*(W(I,J+1)-W(I,J  ))
     6+(RHOW(J-1)*W(I,J-1)+RHOW(J  )*W(I  ,J  ))*(W(I,J  )-W(I,J-1)))
     7   + DTS*RTBW(J)  * ( T(I,J) - TB(J)  + T(I,J-1)  - TB(J-1)  )
     8   + DTS*RQVBW(J) * (QV(I,J) - QVB(J) + QV(I,J-1) - QVB(J-1) )
     9   - DTSG * ( QL(I,J) +QL(I,J-1) )
  555 CONTINUE
C
C       FORCING FOR  V, T, QV, QL  EQUATIONS
C
      DO 540 J = 1 , N
      DO  I = 2 , MM1
      IF(J.EQ.1)THEN
       VMINUS=V(I,J)
       TMINUS=T(I,J)
       TBMINUS=TB(J)
       QVMINUS=QV(I,J)
       QLMINUS=QL(I,J)
      ELSE
       VMINUS=V(I,J-1)
       TMINUS=T(I,J-1)
       TBMINUS=TB(J-1)
       QVMINUS=QV(I,J-1)
       QLMINUS=QL(I,J-1)
      END IF
      VA(I,J)    = VA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( V(I+1,J  ) - V(I  ,J  ) )
     2            +R(I  ) * U(I  ,J) * ( V(I  ,J  ) - V(I-1,J  ) ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( V(I  ,J+1) - V(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( V(I  ,J  ) - VMINUS ) )
     5 -DTL2*(F+ V(I,J)/RS(I))*(R(I+1)*U(I+1,J)+R(I)*U(I,J))/RS(I)
      TA(I,J)    = TA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( T(I+1,J  ) - T(I  ,J  ) )
     2            +R(I  ) * U(I  ,J) * ( T(I  ,J  ) - T(I-1,J  ) ) )
c     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* (
c     & T(I  ,J+1) -TB(J+1) - T(I  ,J  ) +TB(J) )
c     4            +RHOW(J  )*W(I,J  )* (
c     & T(I  ,J  ) -TB(J)   - TMINUS +TBMINUS ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( T(I  ,J+1) - T(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( T(I  ,J  ) - TMINUS ) )
     5 +QRAD(I,J)
      QVA(I,J)    = QVA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( QV(I+1,J  ) - QV(I  ,J  ) )
     2            +R(I  ) * U(I  ,J) * ( QV(I  ,J  ) - QV(I-1,J  ) ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( QV(I  ,J+1) - QV(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( QV(I  ,J  ) - QVMINUS ) )
      QLA(I,J)    = QLA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( QL(I+1,J  ) - QL(I  ,J  ) )
     2            +R(I  ) * U(I  ,J) * ( QL(I  ,J  ) - QL(I-1,J  ) ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( QL(I  ,J+1) - QL(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( QL(I  ,J  ) - QLMINUS ) )
      END DO
c
c  Inner boundary
c
      I=1
c
      IF(J.EQ.1)THEN
       VMINUS=V(I,J)
       TMINUS=T(I,J)
       TBMINUS=TB(J)
       QVMINUS=QV(I,J)
       QLMINUS=QL(I,J)
      ELSE
       VMINUS=V(I,J-1)
       TMINUS=T(I,J-1)
       TBMINUS=TB(J-1)
       QVMINUS=QV(I,J-1)
       QLMINUS=QL(I,J-1)
      END IF
      VA(I,J)    = VA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( V(I+1,J  ) - V(I  ,J  ) ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( V(I  ,J+1) - V(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( V(I  ,J  ) - VMINUS ) )
     5 -DTL2*(F+ V(I,J)/RS(I))*(R(I+1)*U(I+1,J)+R(I)*U(I,J))/RS(I)
      TA(I,J)    = TA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( T(I+1,J  ) - T(I  ,J  ) ) )
c     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* (
c     & T(I  ,J+1) -TB(J+1) - T(I  ,J  ) +TB(J) )
c     4            +RHOW(J  )*W(I,J  )* (
c     & T(I  ,J  ) -TB(J)   - TMINUS +TBMINUS ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( T(I  ,J+1) - T(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( T(I  ,J  ) - TMINUS ) )
     5 +QRAD(I,J)
      QVA(I,J)    = QVA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( QV(I+1,J  ) - QV(I  ,J  ) ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( QV(I  ,J+1) - QV(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( QV(I  ,J  ) - QVMINUS ) )
      QLA(I,J)    = QLA(I,J)
     1 -DTLR(I)* ( R(I+1) * U(I+1,J) * ( QL(I+1,J  ) - QL(I  ,J  ) ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( QL(I  ,J+1) - QL(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( QL(I  ,J  ) - QLMINUS ) )
c
  540 CONTINUE
C
C       OUTER BOUNDARY
C
      DO 530 J = 1 , N
      UB(J) =  MAX( R(MP1) * U(MP1,J) + R(M) * U(M,J) , 0. )
  530 CONTINUE
      I = M
      DO 545 J = 1 , N
      IF(J.EQ.1)THEN
       VMINUS=V(I,J)
       TMINUS=T(I,J)
       TBMINUS=TB(J)
       QVMINUS=QV(I,J)
       QLMINUS=QL(I,J)
      ELSE
       VMINUS=V(I,J-1)
       TMINUS=T(I,J-1)
       TBMINUS=TB(J-1)
       QVMINUS=QV(I,J-1)
       QLMINUS=QL(I,J-1)
      END IF
      VA(I,J)    = VA(I,J)
     1 -DTLR(I) * UB(J) * ( V1(I  ,J  ) - V1(I-1,J  ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( V(I  ,J+1) - V(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( V(I  ,J  ) - VMINUS ) )
     5 -DTL2*(F+ V(I,J)/RS(I))*(R(I+1)*U(I+1,J)+R(I)*U(I,J))/RS(I)
      TA(I,J)    = TA(I,J)
     1 -DTLR(I) * UB(J) * ( T1(I  ,J  ) - T1(I-1,J  ) )
c     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* (
c     & T(I  ,J+1) -TB(J+1) - T(I  ,J  ) +TB(J) )
c     4            +RHOW(J  )*W(I,J  )* (
c     & T(I  ,J  ) -TB(J)   - TMINUS +TBMINUS ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( T(I  ,J+1) - T(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( T(I  ,J  ) - TMINUS ) )
     5 +QRAD(I,J)
      QVA(I,J)    = QVA(I,J)
     1 -DTLR(I) * UB(J) * ( QV1(I  ,J  ) - QV1(I-1,J  ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( QV(I  ,J+1) - QV(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( QV(I  ,J  ) - QVMINUS ) )
      QLA(I,J)    = QLA(I,J)
     1 -DTLR(I) * UB(J) * ( QL1(I  ,J  ) - QL1(I-1,J  ) )
     3 -DTLZ(J)* ( RHOW(J+1)*W(I,J+1)* ( QL(I  ,J+1) - QL(I  ,J  ) )
     4            +RHOW(J  )*W(I,J  )* ( QL(I  ,J  ) - QLMINUS ) )
  545 CONTINUE
C
C         RAIN (AND SNOW) FALL
C
      DO 570 J = 1 , N
      DO 570 I = 1 , M
c
      VTERM2D(I,J)=VTERM
c
      IF(VTERM.LT.1.0E-6)THEN ! If VTERM = 0 calculate fall speed
c
       VTERM2D(I,J)=36.34*SQRT(1.1/RHOT(J))*(RHOT(J)*QL(I,J))**0.1364
       PNST= PN(J) + P1(I,J) - P1(M,1)
       PDST= 1000. * PNST ** (1./XKAPPA)
       TEMP= PNST * T1(I,J) -273.0
       IF(TEMP.LT.0.0)THEN  !  Weight toward set fall speed of snow ast T approaches -20C
         SWEIGHT=(TEMP+20.0)/20.0
         SWEIGHT=MAX(SWEIGHT,0.0)
         VTERM2D(I,J)=SWEIGHT*VTERM2D(I,J)+VTERMSNOW*(1.-SWEIGHT)
       END IF
       VTERM2D(I,J)=MIN(VTERM2D(I,J),12.0)  ! Put cap on fall speed of rain
c
      END IF
c
      DUM(I,J)=0.
      IF((QL(I,J) - ACT).GT.0.0) DUM(I,J)=VTERM2D(I,J)
  570 CONTINUE
      DO 575 J = 2 , NM1
      DO 575 I = 1 , M
      QLA(I,J) = QLA(I,J) + DTLZ(J) *
     $(DUM(I,J+1)*RHOT(J+1)*QL(I,J+1)-DUM(I,J-1)*RHOT(J-1)*QL(I,J-1))
  575 CONTINUE
c
      DO 580 I = 1 , M
      DUM(I,2)=0.
c
      VTERM2D(I,J)=VTERM
c
      IF(VTERM.LT.1.0E-6)THEN ! If VTERM = 0 calculate fall speed
c
      VTERM2D(I,2)=36.34*SQRT(1.1/RHOT(2))*(RHOT(2)*QL1(I,2))**0.1364
      VTERM2D(I,1)=36.34*SQRT(1.1/RHOT(1))*(RHOT(2)*QL1(I,1))**0.1364
c
      END IF
c
      DUM(I,2)=0.0
      IF((QL1(I,2) - ACT).GT.0.0) DUM(I,2)=VTERM2D(I,2)
      DUM(I,1)=0.
      IF((QL1(I,1) - ACT).GT.0.0) DUM(I,1)=VTERM2D(I,1)
c
      QLA(I,1) = QLA(I,1) + 2.*DTLZ(1)*
     $( DUM(I,2)*RHOT(2)*QL1(I,2)  - DUM(I,1)*RHOT(1)*QL1(I,1 ) )
c
  580 CONTINUE
C
C        TIME SMOOTHER
C
c      cfnew= dts/dtl

      DO 160 J = 1 , N
      DO 160 I = 1 , M
      U(I+1,J)= U(I+1,J)+ EPS * ( U1(I+1,J)-2.*U(I+1,J) )
      V(I,J)  = V(I,J)  + EPS * ( V1(I,J)  -2.*V(I,J)   )
      W(I,J+1)= W(I,J+1)+ EPS * ( W1(I,J+1)-2.*W(I,J+1) )
      P(I,J)  = P(I,J)  + EPS * ( P1(I,J)  -2.*P(I,J)   )
      T(I,J)  = T(I,J)  + EPS * ( T1(I,J)  -2.*T(I,J)   )
      QV(I,J) = QV(I,J) + EPS * ( QV1(I,J) -2.*QV(I,J)  )
      QL(I,J) = QL(I,J) + EPS * ( QL1(I,J) -2.*QL(I,J)  )
  160 CONTINUE
C
C         SMALL TIME STEP
C
      DO 170 NSMALL = 1 , NS
      DO 171 J = 1 , N
      DO 171 I = 2 , M
c  Modified by K.E. 4/15/02
c  171 U1(I,J)=U1(I,J)-CPTDR(J)*(P1(I,J)-P1(I-1,J))+UA(I,J)
        CFAC=DTS*RDR*CP*0.5*(T1(I,J)*(1.+0.61*QV(I,J))+T1(I-1,J)*
     1 (1.+0.61*QV(I-1,J)))
  171 U1(I,J)=U1(I,J)-CFAC*(P1(I,J)-P1(I-1,J))+UA(I,J)
c      do 1761, j=1,N
c      do 1761, i=1,M
c        t1(i,j)= t1(i,j) +cfnew*( ta(i,j)
c     & -DTLZ(J)* ( RHOW(J+1)*W1(I,J+1)* ( tb(j+1) - tb(j) )
c     &            +RHOW(J  )*W1(I,J  )* ( tb(j) - tb(j-1) ) ) )
c 1761 continue
      DO 172 J = 2 , N
      DO 172 I = 1 , M
c  Modified by K.E. 4/15/02
C  172 WS(I,J)=W1(I,J)-.5*(1.-EP)*CPTDZ(J)*(P1(I,J)-P1(I,J-1))+WA(I,J)
        CFAC=DTS*RDZ*CP*0.5*(T1(I,J)+T1(I,J-1))*
     1   (1.+0.305*(QV(I,J)+QV(I,J-1)))
  172 WS(I,J)=W1(I,J)-.5*(1.-EP)*CFAC*(P1(I,J)-P1(I,J-1))+WA(I,J)
c     7   + DTS*RTBW(J)  * ( T1(I,J) - TB(J)  + T1(I,J-1)  - TB(J-1))
      DO 173 J = 1 , N
      DO 173 I = 1 , M
      PS(I,J)=P1(I,J)-RC2(J) *
     1   ( R(I+1)*U1(I+1,J)  -R(I)*U1(I,J)  )/RS(I)
     2               -.5*(1.-EP)*ZC2(J) *
     3    ( RHOTVW(J+1)*W1(I,J+1)-RHOTVW(J)*W1(I,J))
  173 CONTINUE
      DO 175 J = 2 , N
      DO 175 I = 1 , M
c  Modified by K.E. 4/15/02
C      D(I,J) = (WS(I,J) - CPTDZ(J)*.5*(1.+EP)*(PS(I,J)-PS(I,J-1))
C     $                + C(J) * D(I,J-1) ) * E(J) /A(J)
       CFAC=DTS*RDZ*CP*0.5*(T1(I,J)+T1(I,J-1))*
     1  (1.+0.305*(QV(I,J)+QV(I,J-1)))
      D(I,J) = (WS(I,J) - CFAC*.5*(1.+EP)*(PS(I,J)-PS(I,J-1))
     $                + C(J) * D(I,J-1) ) * E(J) /A(J)
  175 CONTINUE
      DO 174 J = N , 1 , -1
       DO 174 I = 1 , M
        W1(I,J) = E(J) * W1(I,J+1) + D(I,J)
        P1(I,J) = PS(I,J) -.5*(1.+EP)*ZC2(J)*
     $    ( RHOTVW(J+1)*W1(I,J+1)-RHOTVW(J)*W1(I,J))
  174  CONTINUE
  170 CONTINUE
C
C        OUTER BOUNDARY
C
      DO 195 J = 1 , N
       U1(MP1,J) = U1(MP1,J) + UA(MP1,J)
  195 CONTINUE
C
C        ADVANCE V , T , QV , QL
C
      DO 490 J=1,N
      DO 490 I=1,M
      V1(I,J)  = V1(I,J)  + VA(I,J)
      T1(I,J)  = T1(I,J)  + TA(I,J)
      QV1(I,J) = QV1(I,J) + QVA(I,J)
      QL1(I,J) = QL1(I,J) + QLA(I,J)
      IF(QV1(I,J).LT.0.   ) QV1(I,J)=0.
      IF(QL1(I,J).LT.0.   ) QL1(I,J)=0.
C
C        CONDENSATION / EVAPORATION
C
      PNST= PN(J) + P1(I,J) - P1(M,1)
      PDST= 1000. * PNST ** (1./XKAPPA)
      TEMP= PNST * T1(I,J)
      ES  = 6.11 * EXP(A1 * ( TEMP  - 273. ) / ( TEMP  - 36. ) )
      QSS = .622 * ES /MAX((PDST-ES),ES)
      IF( QV1(I,J) .LT. QSS .AND. QL1(I,J).LE. 1.E-8 ) GO TO 480
      R1 = 1./(1. + XLDCP * 237. * A1 * QSS /( TEMP - 36. )**2 )
      QVD = R1 * ( QV1(I,J) - QSS )
c
      IF(QVD.LT.0.0)THEN
C
      IF( EVAP .EQ. 'N' .OR. EVAP .EQ. 'n')THEN ! No Evaporation
C
        QVD=0.
C
      ELSEIF(  EVAP .EQ. 'I' .OR. EVAP .EQ. 'i' )THEN ! Instantaneous evaporation (per RE(1987))
c
      IF ( (QL1(I,J)   + QVD) .LT. 0. ) THEN
C
C         EVAPORATE
C
      T1 (I,J) = T1(I,J)  - XLDCP * QL1(I,J) /PNST
      QV1(I,J) = QV1(I,J) + QL1(I,J)
      QL1(I,J) = 0.
c
      ENDIF
c
      ELSE
C
C         EVAPORATE
C
c  Added July 3 2020
c
      eamount=0.5*DTL*(1.-QV1(I,J)/QSS)*sqrt(max(QL1(I,J),0.0))/ ! Evaporation scheme (Emanuel, 1991)
     1  (2.0e3+1.0e4/(PDST*QSS))
      eamount=min(eamount,QL1(I,J))
c
      T1 (I,J) = T1(I,J)  - XLDCP *eamount  /PNST
      QV1(I,J) = QV1(I,J) + eamount
      QL1(I,J) = QL1(I,J) - eamount
c
      END IF
c
      ELSE
c
C         CONDENSE
C
      T1(I,J)   = T1(I,J)  + XLDCP * QVD / PNST
      QV1(I,J)  = QV1(I,J)         - QVD
      QL1(I,J)  = QL1(I,J)         + QVD
c
      END IF
c
  480 CONTINUE
c
      IF(QV1(I,J).LT.0.   ) QV1(I,J)=0.
      IF(QL1(I,J).LT.0.   ) QL1(I,J)=0.
  490 CONTINUE
C
C        TIME FLIP
C
      DO 190 J=1,N
      DO 190 I=1,M
      UTEMP    = U(I+1,J) + EPS * U1(I+1,J)
      VTEMP    = V(I,J)   + EPS * V1(I,J)
      WTEMP    = W(I,J+1) + EPS * W1(I,J+1)
      PTEMP    = P(I,J)   + EPS * P1(I,J)
      TTEMP    = T(I,J)   + EPS * T1(I,J)
      QVTEMP   = QV(I,J)  + EPS * QV1(I,J)
      QLTEMP   = QL(I,J)  + EPS * QL1(I,J)

      U(I+1,J) = U1(I+1,J)
      V(I,J)   = V1(I,J)
      W(I,J+1) = W1(I,J+1)
      P(I,J)   = P1(I,J)
      T(I,J)   = T1(I,J)
      QV(I,J)  = QV1(I,J)
      QL(I,J)  = QL1(I,J)

      U1(I+1,J)= UTEMP
      V1(I,J)  = VTEMP
      W1(I,J+1)= WTEMP
      T1(I,J)  = TTEMP
      P1(I,J)  = PTEMP
      QV1(I,J) = QVTEMP
      QL1(I,J) = QLTEMP
  190 CONTINUE
      IF( mod(itime,ipl-iave+1).eq.0) then
                NPLOT=NPLOT+1
        IF(NPLOT.NE.1)THEN
            nmplot=nplot-1
            WRITE(ANP,191)NMPLOT
  191 FORMAT(I2.2)
      OPEN(UNIT=15,FILE='output/vcon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=16,FILE='output/ucon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=18,FILE='output/wcon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=19,FILE='output/tcon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=20,FILE='output/pcon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=21,FILE='output/tecon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=25,FILE='output/liqcon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=28,FILE='output/tescon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=29,FILE='output/qcon'//anp//'.out',STATUS='UNKNOWN')
      OPEN(UNIT=30,FILE='output/tfcon'//anp//'.out',STATUS='UNKNOWN')
        TT = TAV(1,1) * ( PN(1) + PAV(1,1) - PAV(1,1) )
        TEMAX=TAV(1,1) * EXP( XLDCP * QVAV(1,1) / TT )
        TEMAX=TEMAX+10.0
        DO J=1,NPL
         WRITE(15,1035)(VAV(I,J),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=.5 * ( UAV(I+1,J) + UAV(I,J) )
         END DO
         WRITE(16,1035)(DUMT(I),I=1,MPL)
       DO I=1,MPL
          DUMT(I)=.5 * ( WAV(I,J+1) + WAV(I,J) )
         END DO
         WRITE(18,1035)(DUMT(I),I=1,MPL)
c       DO I=1,MPL
c          DUMT(I)=TAV(I,J)*(PN(J)+PAV(I,J)-PAV(M,1)) - TB(J) * PN(J)
c          DUMT(I)=TAV(I,J)*(PN(J)+PAV(I,J)-PAV(M,1))
c         END DO
c         WRITE(19,1035)(DUMT(I),I=1,MPL)
       DO I=1,MPL
          DUMT(I)=1000.*QLAV(I,J)
         END DO
         WRITE(25,1035)(DUMT(I),I=1,MPL)
c       DO I=1,MPL
c          DUMT(I)=1000.*PN(J)**(1./XKAPPA) * (PAV(I,J)-PAV(M,1))/XKAPPA
c         END DO
c         WRITE(20,1035)(DUMT(I),I=1,MPL)
       DO I=1,MPL
        TT = TAV(I,J) * ( PN(J) + PAV(I,J) - PAV(M,1) )
          DUMT(I)=TAV(I,J) * EXP( XLDCP * QVAV(I,J) / TT )
c          DUMT(I)=MIN(DUMT(I),TEMAX)
         END DO
       WRITE(21,1035)(DUMT(I),I=1,MPL)
       DO I=1,MPL
        TT = TAV(I,J) * ( PN(J) + PAV(I,J) - PAV(M,1) )
          TTC=TT-273.15
          ESTT=6.112*EXP(17.67*TTC/(243.5+TTC))
          PFULL=1000.*( PN(J) + PAV(I,J)- PAV(M,1))**(1./XKAPPA)
          QSTT=0.622*ESTT/MAX(1.0,(PFULL-ESTT))
          DUMT(I)=TAV(I,J) * EXP( XLDCP * QSTT / TT )
c          DUMT(I)=MIN(DUMT(I),TEMAX)
         END DO
         WRITE(28,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          TT = TAV(I,J) * ( PN(J) + PAV(I,J) - PAV(M,1) )
          DUMT(I)=TT
         END DO
         WRITE(30,1035)(DUMT(I),I=1,MPL)
        END DO
c
        DO I=1,M
         DO J=1,N
          DUMT(J)=TAV(I,J)*(PN(J)+PAV(I,J)-PAV(M,1))
         END DO
         WRITE(19,1035)(DUMT(J),J=1,N)
         DO J=1,N
          DUMT(J)=1000.*( PN(J) + PAV(I,J)- PAV(M,1))**(1./XKAPPA)
         END DO
         WRITE(20,1035)(DUMT(J),J=1,N)
         DO J=1,N
          DUMT(J)=1000.*QVAV(I,J)
         END DO
         WRITE(29,1035)(DUMT(J),J=1,N)
        END DO
c
        CLOSE(15)
        CLOSE(16)
        CLOSE(18)
        CLOSE(19)
        CLOSE(20)
        CLOSE(21)
        CLOSE(25)
        CLOSE(28)
        CLOSE(29)
        CLOSE(30)
        END IF
c
      DO 375 J = 1 , N
      DO 375 I = 1 , M
      UAV(I+1,J) = 0.0
      VAV(I,J)   = 0.0
      WAV(I,J+1) = 0.0
      TAV(I,J)   = 0.0
      PAV(I,J)   = 0.0
      QVAV(I,J)  = 0.0
      QLAV(I,J)  = 0.0
      XKMAV(I,J) = 0.0
 375  CONTINUE
      ipl=ipl+iplot
      end if
      if(itime.ge.ipl-iplot-iave+1.and.itime.le.ipl-iplot) then
      DO 370 I = 1 , M
      DO 370 J = 1 , N
      UAV(I+1,J) = UAV(I+1,J) + U1(I+1,J)/AVE
      VAV(I,J)   = VAV(I,J)   + V1(I,J)/AVE
      WAV(I,J+1) = WAV(I,J+1) + W1(I,J+1)/AVE
      TAV(I,J)   = TAV(I,J)   + T1(I,J)/AVE
      PAV(I,J)   = PAV(I,J)   + P1(I,J)/AVE
      QVAV(I,J)  = QVAV(I,J)  + QV1(I,J)/AVE
      QLAV(I,J)  = QLAV(I,J)  + QL1(I,J)/AVE
      XKMAV(I,J) = XKMAV(I,J) + XKM(I,J)/AVE
  370 CONTINUE
      end if
C
C      COMPUTE WMAX , VMAX
C
      IWMAX=1
      JWMAX=1
      WMAX=0.
      IVMAX=1
      IVMAXS=1
      JVMAX=1
      VMAX=0.
      UMIN=0.0
       VMAXS=0.0
       DO I=1,M
         IF(V(I,1).GT.VMAXS)THEN
        VMAXS=V(I,1)
         IVMAXS=I
         END IF
       END DO
      DO 210 J = 2 , N
      DO 210 I = 1 , M
      IF(W(I,J).GT.WMAX) THEN
      WMAX = W(I,J)
      IWMAX = I
      JWMAX = J
      END IF
      IF(V(I,J).GT.VMAX) THEN
      VMAX = V(I,J)
      IVMAX = I
      JVMAX = J
      END IF
      UMIN=MIN(UMIN,U(I,2))
  210 CONTINUE
      TINHRS = TIME / 3600.
      IF( ITIME.EQ.ISTART ) GO TO 301
      IF(MOD(ITIME,ITMAX).EQ.0)
     $ write(17,925) TINHRS,WMAX,IWMAX,JWMAX,VMAX,IVMAX,JVMAX
  925 FORMAT(1H ,'TIME = ',F9.1,2X,'WMAX = ',E11.4,2X,'IW= ',I3,
     $2X,'JW= ',I2,2X,'VMAX = ',E11.4,2X,'IV= ',I3,2X,'JV= ',I3)
C
C        PRINT FIELDS
C
      IF(MOD(ITIME,IPRINT).NE.0) GO TO 200
  301 CONTINUE
      TINHRS = TIME / 3600.
      WRITE(17, 915) TINHRS
      WRITE(17, 900) 'U'
      DO 300 J=N,1,-1
  300 WRITE(17, 910) (U(I,J)  ,I=1,20,2)
      WRITE(17, 900) 'W'
      DO 305 J=N,2,-1
  305 WRITE(17, 910) (W(I,J)  ,I=1,20,2)
      WRITE(17, 900) 'V'
      DO 306 J =N,1,-1
  306 WRITE(17, 910) (V(I,J)  ,I=1,20,2)
      WRITE(17, 900) 'THTA'
      DO 310 J=N,1,-1
  310 WRITE(17, 910) ( T(I,J) ,I=1,20,2 )
      WRITE(17, 900) 'TSUR'
      WRITE(17, 910) (TSURF(I),I=1,20,2)
      WRITE(17, 900) 'TEMP'
      DO 311 J=N,1,-1
  311 WRITE(17, 910) ( T(I,J)*(PN(J)+P(I,J)-P(M,1))-TB(J)*PN(J)
     &                ,I=1,20,2 )
      WRITE(17, 900) 'QV'
      DO 315 J=N,1,-1
  315 WRITE(17, 910) (1000.*QV(I,J)  ,I=1,20,2)
      WRITE(17, 900) 'QSUR'
      WRITE(17, 910) (1000.*QSURF(I) ,I=1,20,2)
      WRITE(17, 900) 'QL'
      DO 316 J=N,1,-1
  316 WRITE(17, 910) (1000.*QL(I,J)  ,I=1,20,2)
      WRITE(17, 900) 'P'
      DO 320 J=N,1,-1
      FAC = PN(J)**1./XKAPPA
  320 WRITE(17, 910) ( 1000.*FAC*(P(I,J)-P(M,1))/XKAPPA ,I=1,20,2)
      WRITE(17, 900) 'KM'
      DO 325 J=N,1,-1
  325 WRITE(17, 910) ( .01*XKM(I,J) ,I=1,20,2)
      WRITE(17, 900) 'TE'
      DO 330 J =N,1,-1
      TT = T(I,J) * ( PN(J) + P(I,J) - P(M,1) )
  330 WRITE(17, 910) (T(I,J) * (1. + XLDCP * QV(I,J) / TT ),I=1,20,2)
  200 CONTINUE
C
      TINHRS = TIME/3600.
      write(chtime,812) tinhrs
  812 format('   TIME=',F10.4)
      IF(ITIME.GT.1) GO TO 499
      DTL = 2. * DT
      NS  = 2  * NS
      DTL2 = .5 * DTL
      DO 17 I = 1 , M
   17 DTLR(I)  = .5  * DTL * RDR /RS(I)
      DO 27 J = 1 , N
   27 DTLZ(J)  = .5  * DTL * RDZ / RHOT(J)
  499 CONTINUE
        IF(TIPL.GE.TIMEPL)THEN
         NPT=NPT+1
         TTP(NPT)=TIME/3600.0
         VMX(NPT)=VMAX
         VMXS(NPT)=VMAXS
         UMN(NPT)=UMIN
         WMX(NPT)=WMAX*100.
         PMN(NPT)= PDS+1000.*PN(1)**(1./XKAPPA) *
     1    (P(1,1)-P(M,1))/XKAPPA
         RMX(NPT)=0.001*DR*FLOAT(IVMAX-1)
          DO J=1,N
           VMAXZ(NPT,J)=0.0
           DO I=1,M
            VMAXZ(NPT,J)=MAX(VMAXZ(NPT,J),V(I,J))
           END DO
          END DO
           JOUT=15000.0/DZ
          DO I=1,MPL
           TT = T(I,1) * ( PN(1) + P(I,1) - P(M,1) )
           THEHOV(NPT,I)= T(I,1) * EXP( XLDCP * QV(I,1) / TT )
           VBHOV(NPT,I)=V(I,1)
           UBHOV(NPT,I)=U(I,1)
           VTHOV(NPT,I)=V(I,JOUT)
           UTHOV(NPT,I)=U(I,JOUT)
          END DO
         TIPL=0.0
        END IF
C
  500 CONTINUE
C
        DO 820 I=1,MPL
         UPL(I)=UAV(I,1)
         VPL(I)=VAV(I,1)
         PPL(I)=PDS+1000.*PN(1)**(1./XKAPPA)*
     1    (PAV(I,1)-PAV(M,1))/XKAPPA
         TT = TAV(I,1) * ( PN(1) + PAV(I,1) - PAV(M,1) )
         THEPL(I)= TAV(I,1) * EXP( XLDCP * QVAV(I,1) / TT )
  820   CONTINUE
  900 FORMAT(/,1X,A4,/)
  910 FORMAT(1H ,16F7.2)
  915 FORMAT(//,1H ,'TIME=',F10.3)
         OPEN(UNIT=13,FILE='output/time.out',STATUS='UNKNOWN')
         OPEN(UNIT=14,FILE='output/radius.out',STATUS='UNKNOWN')
         OPEN(UNIT=15,FILE='output/vcon.out',STATUS='UNKNOWN')
         OPEN(UNIT=16,FILE='output/ucon.out',STATUS='UNKNOWN')
         OPEN(UNIT=18,FILE='output/wcon.out',STATUS='UNKNOWN')
         OPEN(UNIT=19,FILE='output/tcon.out',STATUS='UNKNOWN')
         OPEN(UNIT=20,FILE='output/pcon.out',STATUS='UNKNOWN')
         OPEN(UNIT=21,FILE='output/tecon.out',STATUS='UNKNOWN')
         OPEN(UNIT=24,FILE='output/vmaxz.out',STATUS='UNKNOWN')
         OPEN(UNIT=25,FILE='output/liqcon.out',STATUS='UNKNOWN')
         OPEN(UNIT=27,FILE='output/thehov.out',STATUS='UNKNOWN')
         OPEN(UNIT=28,FILE='output/tescon.out',STATUS='UNKNOWN')
         OPEN(UNIT=29,FILE='output/s.out',STATUS='UNKNOWN')
         OPEN(UNIT=30,FILE='output/tfcon.out',STATUS='UNKNOWN')
         OPEN(UNIT=31,FILE='output/xkcon.out',STATUS='UNKNOWN')
          OPEN(UNIT=32,FILE='output/vbhov.out',STATUS='UNKNOWN')
          OPEN(UNIT=33,FILE='output/vthov.out',STATUS='UNKNOWN')
          OPEN(UNIT=34,FILE='output/ubhov.out',STATUS='UNKNOWN')
          OPEN(UNIT=35,FILE='output/uthov.out',STATUS='UNKNOWN')
        DO 1010 I=1,NPT
         WRITE(13,1005)TTP(I),VMX(I),UMN(I),WMX(I),PMN(I),RMX(I),
     1    VMAXTH,VMXS(I)
          WRITE(24,2006)(VMAXZ(I,J),J=1,NPL)
 1005  FORMAT(2X,7(F12.6,' '),F12.6)
 1010        CONTINUE
 2006 FORMAT(50(1X,F12.6))
        DO 1020 I=1,MPL
         WRITE(14,1015)RPL(I),VPL(I),UPL(I),PPL(I),THEPL(I)
 1015         FORMAT(2X,4(F12.6,' '),F12.6)
 1020        CONTINUE
        TT = TAV(1,1) * ( PN(1) + PAV(1,1) - PAV(1,1) )
        TEMAX=TAV(1,1) * EXP( XLDCP * QVAV(1,1) / TT )
        TEMAX=TEMAX+10.0
c
        WRITE(29,1031)
 1031   FORMAT(9X,'Sea Surface Temperature (C)')
        WRITE(29,1032)
 1032   FORMAT(9X,'---------------------------')
        WRITE(29,1033)TBSW
 1033   FORMAT(19X,F5.2)
        WRITE(29,*)
        WRITE(29,1034)
 1034   FORMAT(4X,'p (hPa)         T (C)            q (g/Kg)')
        WRITE(29,1036)
 1036   FORMAT(4X,'-----------------------------------------')
        WRITE(29,*)
C
        DO 1040 J=1,NPL
         WRITE(15,1035)(VAV(I,J),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=.5 * ( UAV(I+1,J) + UAV(I,J) )
         END DO
         WRITE(16,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=.5 * ( WAV(I,J+1) + WAV(I,J) )
         END DO
         WRITE(18,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=.5 * ( XKMAV(I,J+1) + XKMAV(I,J) )
         END DO
         WRITE(31,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=TAV(I,J)*(PN(J)+PAV(I,J)-PAV(M,1)) - TB(J) * PN(J)
         END DO
         WRITE(19,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=1000.*QLAV(I,J)
         END DO
         WRITE(25,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          DUMT(I)=1000.*PN(J)**(1./XKAPPA)*(PAV(I,J)-PAV(M,1))/XKAPPA
         END DO
         WRITE(20,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          TT = TAV(I,J) * ( PN(J) + PAV(I,J) - PAV(M,1) )
          DUMT(I)=TAV(I,J) * EXP( XLDCP * QVAV(I,J) / TT )
c          DUMT(I)=MIN(DUMT(I),TEMAX)
         END DO
         WRITE(21,1035)(DUMT(I),I=1,MPL)
         DO I=1,MPL
          TT = TAV(I,J) * ( PN(J) + PAV(I,J) - PAV(M,1) )
          TTC=TT-273.15
          ESTT=6.112*EXP(17.67*TTC/(243.5+TTC))
          PFULL=1000.*( PN(J) + PAV(I,J)- PAV(M,1))**(1./XKAPPA)
          QSTT=0.622*ESTT/MAX(1.0,(PFULL-ESTT))
          DUMT(I)=TAV(I,J) * EXP( XLDCP * QSTT / TT )
c          DUMT(I)=MIN(DUMT(I),TEMAX)
         END DO
         WRITE(28,1035)(DUMT(I),I=1,MPL)
         TT = TAV(M,J) * ( PN(J) + PAV(M,J) - PAV(M,1) )-273.15
         PFULL=1000.*( PN(J) + PAV(M,J)- PAV(M,1))**(1./XKAPPA)
          WRITE(29,1041) PFULL,TT,1000.*QVAV(M,J)
         DO I=1,MPL
          TT = TAV(I,J) * ( PN(J) + PAV(I,J) - PAV(M,1) )
          DUMT(I)=TT
         END DO
         WRITE(30,1035)(DUMT(I),I=1,MPL)
 1040    CONTINUE
        DO I=1,NPT
         WRITE(27,1035)(THEHOV(I,J),J=1,MPL)
         WRITE(32,1035)(vbhov(I,J),J=1,MPL)
         WRITE(33,1035)(vthov(I,J),J=1,MPL)
         WRITE(34,1035)(ubhov(I,J),J=1,MPL)
         WRITE(35,1035)(uthov(I,J),J=1,MPL)
        END DO
 1035   FORMAT(2X,F12.7,' ',800(F8.2,' '))
 1041   FORMAT(4X,F8.4,5X,F9.4,11X,F8.5)
C
      STOP
      END
C
      SUBROUTINE DIFFUSE
C
      PARAMETER(M=400,N=80,MP1=M+1,MM1=M-1,NP1=N+1,NM1=N-1)
      DIMENSION
     $ TRR(M,N  )  , TTT(M,N)  , TZZ(M,NP1),
     $ TRZ(MP1,NP1), TRT(MP1,N), TZT(M,NP1),
     $ TR(MP1, N)  , TZ(M,NP1) ,
     $ QVR(MP1, N) , QVZ(M,NP1), QLR(MP1,N), QLZ(M,NP1),
     $ XKMH(M,NP1) , TE(M,N)   , UDR(M,N)  , WDZ(M,N)  , VDR(MP1,N),
     $ UZWR(MP1,N) , VDZ(M,N)  , DEFH2(M,N), DEF2(M,N) , DTDZ(M,N)
      COMMON/DIFF/
     $ UA(MP1,N),WA(M,NP1),VA(M,N),TA(M,N),QVA(M,N),QLA(M,N),
     $ U1(MP1,N),V1(M,N),W1(M,NP1),T1(M,N), QV1(M,N), QL1(M,N),
     $ XKM(M,NP1),P1(M,N),R(MP1),RS(M),Z(NP1),ZS(N),TAD(M,N),
     $ RTBW(N),RQVBW(N),TAU(N),TB(N),QVB(N),PN(NP1),PD(NP1),
     $ TSURF(M),QSURF(M),RDR,RDZ,RDR2,RDZ2,DTS,DTL,DZ,
     $ XVL2, XHL2, CD, CE, XLDCP, XKAPPA, A1, G,
     $ CERS(M), CDRS(M), DISS, VCAP, IVMAXS
      DO 500 J = 1 , N
      DO 500 I = 1 , M
      TE(I,J)= T1(I,J) * (1.+ XLDCP * QV1(I,J) /
     $        ( T1(I,J) * ( PN(J) + P1(I,J) - P1(M,1) )  )  )
  500 CONTINUE
      DO 110 J = 1 , N
      DO 110 I = 1 , M
      UDR(I,J) = RDR * ( U1(I+1,J) - U1(I,J) )
      WDZ(I,J) = RDZ * ( W1(I,J+1) - W1(I,J) )
  110 CONTINUE
      DO 120 J = 1 , N
      DO 125 I = 2 , M
      VDR(I,J) = R(I) * RDR * ( V1(I,J)/RS(I) - V1(I-1,J)/RS(I-1) )
  125 CONTINUE
      VDR(1  ,J) = 0.
      VDR(MP1,J) = VDR(M,J)
  120 CONTINUE
      DO 130 J = 2 , N
      DO 135 I = 2 , M
      UZWR(I,J) = RDZ*(U1(I,J)-U1(I,J-1)) + RDR*(W1(I,J)-W1(I-1,J))
  135 CONTINUE
      UZWR(1  ,J) = 0.
      UZWR(MP1,J) = UZWR(M,J)
  130 CONTINUE
      DO 140 J = 2 , N
      DO 140 I = 1 , M
      VDZ(I,J) = RDZ * ( V1(I,J) - V1(I,J-1) )
  140 CONTINUE
      DO 150 J = 2 , N
      DO 150 I = 1 , M
      DEFH2(I,J) = UDR(I,J)*UDR(I,J) + UDR(I,J-1)*UDR(I,J-1)
     $ + ( .25/(RS(I)*RS(I)) )*(
     $  (U1(I+1,J  )+U1(I,J  )) *  (U1(I+1,J  )+U1(I,J  ))
     $ +(U1(I+1,J-1)+U1(I,J-1)) *  (U1(I+1,J-1)+U1(I,J-1)) )
     $ + .25*( VDR(I+1,J  )*VDR(I+1,J  ) + VDR(I,J  )*VDR(I,J  )
     $        +VDR(I+1,J-1)*VDR(I+1,J-1) + VDR(I,J-1)*VDR(I,J-1) )
  150 CONTINUE
      DO 160 J = 2 , N
      DO 160 I = 1 , M
      DEF2(I,J) = DEFH2(I,J) +
     $    WDZ(I,J)*WDZ(I,J) + WDZ(I,J-1)*WDZ(I,J-1)
     $ +.5*( UZWR(I+1,J)*UZWR(I+1,J) + UZWR(I,J)*UZWR(I,J) )
     $ + VDZ(I,J) * VDZ(I,J)
c
c   Line below added 9/29/2010 to change effective Richardson Number
c
c	DEF2(I,J)=2.*DEF2(I,J)
c
  160 CONTINUE
      DO 170 J = 2 , N
      DO 170 I = 1 , M
      DTDZ(I,J) = RDZ*(
     $ 2.*RTBW(J) * (T1(I,J) - T1(I,J-1)  )
     $+2.*RQVBW(J)* (QV1(I,J)- QV1(I,J-1) ) )
      IF( (QL1(I,J) * QL1(I,J-1)) .GE. 1.E-8 ) THEN
      TT = .5 * (  T1(I,J  ) * ( PN(J  )+P1(I,J  )-P1(M,1) )
     $           + T1(I,J-1) * ( PN(J-1)+P1(I,J-1)-P1(M,1) )  )
      AA=(1. + 4362.   *(QV1(I,J)+QV1(I,J-1)) / TT      )
     $  /(1. + 6738953.*(QV1(I,J)+QV1(I,J-1))/(TT*TT)   )
     $  * (2. / (TB(J)+TB(J-1)) )
      DTDZ(I,J) = RDZ* G * ( AA * ( TE(I,J) - TE(I,J-1) )
     $ - ( QL1(I,J) + QV1(I,J) - QL1(I,J-1) - QV1(I,J-1) ) )
      END IF
  170 CONTINUE
      DO 180 J = 2 , N
      DO 180 I = 1 , M
      XKM(I,J) = 0.
c
c  Factor of 0.5 inserted here to mimic larger Richardson Number, 9/29/2010
c
      IF(DEF2(I,J).GT.DTDZ(I,J))
     $  XKM(I,J)= XVL2 * SQRT( DEF2(I,J) - DTDZ(I,J))
      XKMH(I,J) = XHL2 * SQRT( DEFH2(I,J) )
      IF(XKM(I,J).GE. .4 * DZ *DZ/DTL) XKM(I,J) = .4 *DZ*DZ/DTL
      IF(XKMH(I,J).LT.XKM(I,J))   XKMH(I,J) = XKM(I,J)
  180 CONTINUE
      DO 190 I = 1 , M
      XKM(I,  1) = XKM(I,2)
      XKM(I,NP1) = XKM(I,N)
      XKMH(I,  1) = XKMH(I,2)
      XKMH(I,NP1) = XKMH(I,N)
  190 CONTINUE
C
C        CALCULATE STRESS
C
C        TRR(M,N)
C
      DO 200 J = 1 , N
      DO 200 I = 1 , M
      TRR(I,J) =  ( XKMH(I,J+1) + XKMH(I,J) ) * UDR(I,J)
  200 CONTINUE
C
C        TTT(M,N)
C
      DO 210 J = 1 , N
      TTT(1,J) = 0.
      DO 210 I = 2 , M
      TTT(I,J) =
     $.5*(XKMH(I,J+1)+XKMH(I,J)+XKMH(I-1,J+1)+XKMH(I-1,J))*U1(I,J)/R(I)
  210 CONTINUE
C
C        TRZ(M,NP1)
C
      DO 220 J=2,N
      DO 220 I=2,M
      TRZ(I,J)=.5* ( XKM(I-1,J) + XKM(I,J) ) *  UZWR(I,J)
  220 CONTINUE
      DO 230 I =2,M
      TRZ(I,1  )=.5*(CDRS(I)+CDRS(I-1))
     $  *U1(I,1)*SQRT(U1(I,1)**2+.25*(V1(I,1)+V1(I-1,1))**2)
      TRZ(I,NP1) = 0.
  230 CONTINUE
      DO 235 J = 1 , NP1
      TRZ(1  ,J) = 0.
      TRZ(MP1,J) = TRZ(M,J) * R(M) / R(MP1)
  235 CONTINUE
C
C        TRT(M,N)
C
      DO 245 J=1,N
      DO 240 I=2,M
      TRT(I,J) =
     $.25*( XKMH(I,J+1)+XKMH(I,J)+XKMH(I-1,J+1)+XKMH(I-1,J))* VDR(I,J)
  240 CONTINUE
      TRT(1  ,J) = 0.
      TRT(MP1,J) = TRT(M,J) * R(M)**2 / R(MP1)**2
  245 CONTINUE
C
C        TZT(M,NP1)
C
      DO 250 I = 1 , M
      DO 255 J = 2 , N
      TZT(I,J) = XKM(I,J) * VDZ(I,J)
  255 CONTINUE
      TZT(I,1  )=CDRS(I)
     $   *V1(I,1)*SQRT(.25*(U1(I,1)+U1(I+1,1))**2+V1(I,1)**2)
      TZT(I,NP1)= 0.
  250 CONTINUE
C
C          TZZ(M,N)
C
      DO 260 J = 1 , N
      DO 260 I = 1 , M
      TZZ(I,J) = ( XKM(I,J+1) + XKM(I,J) )*WDZ(I,J)
  260 CONTINUE
C
C       TEMPERATURE FLUX
C
C          TR(MP1,N)
C
      DO 270 J=1,N
      DO 275 I=2,M
      TR(I,J) = .25*( XKMH(I,J+1)+XKMH(I,J)+XKMH(I-1,J+1)+XKMH(I-1,J))
     $                          *RDR*(T1(I,J)-T1(I-1,J))
  275 CONTINUE
      TR(1  ,J)=0.0
      TR(MP1,J)= TR(M,J) * R(M) / R(MP1)
  270 CONTINUE
C
C         TZ(M,NP1)
C
      DO 280 I=1,M
      DO 285 J=2,N
      TZ(I,J) = XKM(I,J) * RDZ * (T1(I,J) - T1(I,J-1))
c   Mod by K.E.  4/12/02
        TZ(I,J)=TZ(I,J)*PD(J)
  285 CONTINUE
        VABS=SQRT(.25 * ( U1(I+1,1) + U1(I,1) ) **2 + V1(I,1)**2 )
c
        VABS=MIN(VABS, VCAP)
        TZ(I,1  ) =( T1(I,1)-TSURF(I))*CERS(I)*VABS
c
c
c   Experiment of 2/4/2016
c
c       IM=MAX(I-1,1)
c       VABSM=SQRT(.25 * ( U1(I,1) + U1(IM,1) ) **2 + V1(IM,1)**2 )
c       TZ(I,1  ) =( T1(IM,1)-TSURF(IM))*CERS(IM)*VABSM
c
c   Mod by K.E.  4/12/02
        TZ(I,1)=TZ(I,1)*PD(1)
      TZ(I,NP1) = 0.
  280 CONTINUE
C
C       QV FLUX
C
C         QVR(MP1,N)
C
      DO 370 J=1,N
      DO 375 I=2,M
      QVR(I,J)= .25*(XKMH(I,J+1)+XKMH(I,J)+XKMH(I-1,J+1)+XKMH(I-1,J))
     $                            *RDR*(QV1(I,J)-QV1(I-1,J))
  375 CONTINUE
      QVR(1  ,J)=0.0
      QVR(MP1,J)= QVR(M,J) * R(M) / R(MP1)
  370 CONTINUE
C
C         QVZ(M,NP1)
C
      DO 380 I=1,M
      DO 385 J=2,N
      QVZ(I,J) = XKM(I,J) * RDZ * (QV1(I,J) - QV1(I,J-1))
c   Mod by K.E.  4/12/02
        QVZ(I,J)=QVZ(I,J)*PD(J)/PN(J)
  385 CONTINUE
      VABS=SQRT(.25 * ( U1(I+1,1) + U1(I,1) ) **2 + V1(I,1)**2 )
c
      VABS=MIN(VABS, VCAP)
      QVZ(I,1  ) =( QV1(I,1)-QSURF(I))*CERS(I)*VABS
c
c
c
c   Experiment of 2/4/2016
c
c       IM=MAX(I-1,1)
c       VABSM=SQRT(.25 * ( U1(IM,1) + U1(IM,1) ) **2 + V1(IM,1)**2 )
c       QVZ(I,1  ) =( QV1(IM,1)-QSURF(IM))*CERS(IM)*VABSM
c
c   Mod by K.E.  4/12/02
      QVZ(I,1)=QVZ(I,1)*PD(1)/PN(1)
      QVZ(I,NP1) = 0.
  380 CONTINUE
C
C       QL FLUX
C
C        QLR(MP1,N)
C
      DO 570 J=1,N
      DO 575 I=2,M
      QLR(I,J)= .25*(XKMH(I,J+1)+XKMH(I,J)+XKMH(I-1,J+1)+XKMH(I-1,J))
     $                          *RDR*(QL1(I,J)-QL1(I-1,J))
  575 CONTINUE
      QLR(1  ,J)=0.0
      QLR(MP1,J)= QLR(M,J) * R(M) / R(MP1)
  570 CONTINUE
C
C         QLZ(M,NP1)
C
      DO 580 I=1,M
      DO 585 J=2,N
      QLZ(I,J) = XKM(I,J) * RDZ * (QL1(I,J) - QL1(I,J-1))
  585 CONTINUE
      QLZ(I,1  ) = 0.
      QLZ(I,NP1) = 0.
  580 CONTINUE
      DO 400 J = 1 , N
      DO 400 I = 2 , M
      UA(I,J)= -DTS*(
     1 - RDR * ( RS(I) * TRR(I,J) - RS(I-1) * TRR(I-1,J) ) / R(I)
     2 + TTT(I,J) / R(I)
     3 - RDZ * ( TRZ(I,J+1) - TRZ(I,J) )   )
     4 + DTS * TAU(J) * U1(I,J)
  400 CONTINUE
      DO 410  J = 2 , N
      DO 410  I = 1 , M
      WA(I,J) = -DTS*(
     1 - RDR * ( R(I+1)*TRZ(I+1,J) - R(I)*TRZ(I,J) ) / RS(I)
     2 - RDZ * ( TZZ(I,J) - TZZ(I,J-1) )   )
     3 + DTS * .5 * ( TAU(J) + TAU(J-1) ) * W1(I,J)
  410 CONTINUE
      DO 420  J = 1 , N
      DO 420  I = 1 , M
      VA(I,J) = -DTL * (
     1 - RDR * (R(I+1)*R(I+1)*TRT(I+1,J)-R(I)*R(I)*TRT(I,J))
     2                                           /(RS(I)*RS(I))
     3 - RDZ * ( TZT(I,J+1) - TZT(I,J) )   )
     4 +DTL * TAU(J) * V1(I,J)
c   Mod by K.E.  4/12/02
      TA(I,J) = -DTL * (
     1 - RDR * ( R(I+1) * TR(I+1,J) - R(I) * TR(I,J) ) / RS(I)
c     2 - RDZ * (  TZ(I,J+1) - TZ(I,J) )    )
     2 - RDZ * (  TZ(I,J+1) - TZ(I,J) )/(0.5*(PD(J)+PD(J+1))))
     3 + DTL * TAU(J) * ( T1(I,J) - TB(J) )
c   Mod by K.E.  4/12/02
      QVA(I,J) = -DTL * (
     1 - RDR * ( R(I+1) * QVR(I+1,J) - R(I) * QVR(I,J) ) / RS(I)
C     2 - RDZ * (  QVZ(I,J+1) - QVZ(I,J) )    )
     2 - RDZ * (  QVZ(I,J+1) - QVZ(I,J) )/
     * (0.5*(PD(J)/PN(J)+PD(J+1)/PN(J+1)))   )
     3 + DTL * TAU(J) *  QV1(I,J)
      QLA(I,J) = -DTL * (
     1 - RDR * ( R(I+1) * QLR(I+1,J) - R(I) * QLR(I,J) ) / RS(I)
     2 - RDZ * (  QLZ(I,J+1) - QLZ(I,J) )    )
     3 + DTL * TAU(J) *  QL1(I,J)
  420 CONTINUE
c
c       Dissipative Heating
c
        DO 430 I=1,M-1
         DO 425 J=2,N-1
C         TA(I,J)=TA(I,J)-(0.25*(U1(I+1,J)+U1(I,J))*(UA(I+1,J)+UA(I,J))
C     1    +V1(I,J)*VA(I,J))/(PN(J)*1005.0)
         TAD(I,J)=DTL*(TZT(I,J+1)*VDZ(I,J+1)+XKM(I,J)*0.25*(UZWR
     1    (I,J)+UZWR(I+1,J))**2+XKMH(I,J+1)*(1./16.)*(VDR(I,J)+VDR(
     2   I+1,J)+VDR(I,J+1)+VDR(I+1,J+1))**2+XKMH(I,J+1)*0.25*(UDR(I,J)+
     3    UDR(I,J+1))**2)/(PN(J)*1005.0)
         TA(I,J)=TA(I,J)+DISS*TAD(I,J)
  425        CONTINUE
         TAD(I,1)=DTL*(RDZ*CDRS(I)*(0.25*(U1(I,1)+U1(I+1,1))**2+
     1    V1(I,1)**2)**1.5+XKMH(I,2)*(1./16.)*(VDR(I,1)+VDR(
     2    I+1,1)+VDR(I,2)+VDR(I+1,2))**2+XKMH(I,2)*0.25*(UDR(I,1)+
     3    UDR(I,2))**2)/(PN(1)*1005.0)
         TA(I,1)=TA(I,1)+DISS*TAD(I,1)
  430	CONTINUE
      RETURN
      END
c-----------------------------------------------------------
C
        SUBROUTINE INTERPOLATE(NLEVELS,PDS,DZ,TZ,QZ,DISS,SST,VMAX)
        PARAMETER(NA=600)
        REAL TZ(NLEVELS),QZ(NLEVELS)
        REAL PZ(NLEVELS), TAZ(NLEVELS),QAZ(NLEVELS)
        REAL Z(NA),T(NA),Q(NA),P(NA)
        OPEN(UNIT=101,FILE='s.in',STATUS='OLD')
        RD=287.
        EPSI=1./0.622
        G=9.8
        READ(101,*)
        READ(101,*)
        READ(101,*)SST
        READ(101,21)
   21   FORMAT(3X,///)
c
        DO 40 I=1,2000
          READ(101,*,END=41)P(I),T(I),Q(I)
          NP=I
   40	CONTINUE
   41   CONTINUE
        CLOSE(101)
        DO I=1,NP
         T(I)=T(I)+273.15
         Q(I)=0.001*Q(I)
        END DO
        Z(1)=0.0
        DO 50 I=2,NP
         TV1=T(I-1)*(1.+EPSI*Q(I-1)-Q(I-1))
         TV2=T(I)*(1.+EPSI*Q(I)-Q(I))
         TVBAR=0.5*(TV1+TV2)
         PBAR=0.5*(P(I-1)+P(I))
         Z(I)=Z(I-1)+(RD*TVBAR/(G*PBAR))*(P(I-1)-P(I))
   50   CONTINUE
        DO 100 I=1,NLEVELS
         AZ=0.5*DZ+DZ*FLOAT(I-1)
         DO 70 J=1,NP
         K=J
         IF(Z(J).GE.AZ)GOTO 80
   70   CONTINUE
   80   CONTINUE
         TM1=T(K-1)*(1000./P(K-1))**(0.286)
         TM=T(K)*(1000./P(K))**(0.286)
         TZ(I)=TM1+(TM-TM1)*(AZ-Z(K-1))/(Z(K)-Z(K-1))
         QZ(I)=Q(K-1)+(Q(K)-Q(K-1))*(AZ-Z(K-1))/(Z(K)-Z(K-1))
C
         QAZ(I)=1000.*QZ(I)
         TAZ(I)=T(K-1)+(T(K)-T(K-1))*(AZ-Z(K-1))/(Z(K)-Z(K-1))-273.15
         PZ(I)=P(K-1)+(P(K)-P(K-1))*(AZ-Z(K-1))/(Z(K)-Z(K-1))
C
  100   CONTINUE
C
        CALL PCMIN(SST,PDS,PZ,TAZ,QAZ,NLEVELS,NLEVELS-1,DISS,
     1        PMIN,VMAX,IFL)
        RETURN
        END
C
       SUBROUTINE PCMIN(SST,PSL,P,T,R,NA,N,DISS,PMIN,VMAX,IFL)
C
C   ***   This subroutine calculates the maximum wind speed        ***
C   ***             and mimimum central pressure                   ***
C   ***    achievable in tropical cyclones, given a sounding       ***
C   ***             and a sea surface temperature.                 ***
C
C  INPUT:   SST: Sea surface temperature in C
C
C           PSL: Sea level pressure (mb)
C
C           P,T,R: One-dimensional arrays of dimension NA
C             containing pressure (mb), temperature (C),
C             and mixing ratio (g/kg). The arrays MUST be
C             arranged so that the lowest index corresponds
C             to the lowest model level, with increasing index
C             corresponding to decreasing pressure. The temperature
C             sounding should extend to at least the tropopause and
C             preferably to the lower stratosphere, however the
C             mixing ratios are not important above the boundary
C             layer. Missing mixing ratios can be replaced by zeros.
C
C           NA: The dimension of P,T and R
C
C           N:  The actual number of points in the sounding
C                (N is less than or equal to NA)
C
C  OUTPUT:  PMIN is the minimum central pressure, in mb
C
C           VMAX is the maximum surface wind speed, in m/s
C                  (reduced to reflect surface drag)
C
C           IFL is a flag: A value of 1 means OK; a value of 0
C              indicates no convergence (hypercane); a value of 2
C              means that the CAPE routine failed.
C
C-----------------------------------------------------------------------------
       REAL T(NA), P(NA), R(NA)
C
C   ***   Adjustable constant: Ratio of C_k to C_D    ***
C
       CKCD=1.0
C
C   ***   Adjustable constant for buoyancy of displaced parcels:  ***
C   ***    0=Reversible ascent;  1=Pseudo-adiabatic ascent        ***
C
      SIG=0.0
C
C   ***  Exponent, b, in assumed profile of azimuthal velocity in eye,   ***
C   ***   V=V_m(r/r_m)^b. Used only in calculation of central pressure   ***
C
       b=2.0
C
C   *** Set level from which parcels lifted   ***
C
       NK=1
C
C   *** Factor to reduce gradient wind to 10 m wind
C
       VREDUC=1.0
C
C   ***   Normalize certain quantities   ***
C
       SSTK=SST+273.15
       ES0=6.112*EXP(17.67*SST/(243.5+SST))
       DO 40 I=1,N
        R(I)=R(I)*0.001
        T(I)=T(I)+273.15
   40       CONTINUE
C
C   ***   Default values   ***
C
      VMAX=0.0
       PMIN=PSL
       IFL=1
C
       NP=0
       PM=950.0
C
C   ***   Find environmental CAPE ***
C
      TP=T(NK)
      RP=R(NK)
      PP=P(NK)
      CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEA,TOA,IFLAG)
      IF(IFLAG.NE.1)IFL=2
C
C   ***   Begin iteration to find mimimum pressure   ***
C
  100 CONTINUE
C
C   ***  Find CAPE at radius of maximum winds   ***
C
      TP=T(NK)
      PP=PM
      RP=0.622*R(NK)*PSL/(PM*(0.622+R(NK))-R(NK)*PSL)
      CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEM,TOM,IFLAG)
      IF(IFLAG.NE.1)IFL=2
      RAT=SSTK/TOM
      RAT=RAT/(DISS+RAT*(1.-DISS))
C
C  ***  Find saturation CAPE at radius of maximum winds   ***
C
      TP=SSTK
      PP=PM
      RP=0.622*ES0/(PM-ES0)
      CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEMS,TOMS,IFLAG)
      IF(IFLAG.NE.1)IFL=2
C
C  ***  Initial estimate of minimum pressure   ***
C
      RS0=RP
      TV1=T(1)*(1.+R(1)/0.622)/(1.+R(1))
       TVAV=0.5*(TV1+SSTK*(1.+RS0/0.622)/(1.+RS0))
       CAT=CAPEM-CAPEA+0.5*CKCD*RAT*(CAPEMS-CAPEM)
C       CAT=0.5*CKCD*RAT*(CAPEMS-CAPEM)
       PNEW=PSL*EXP(-CAT/(287.04*TVAV))
C
C   ***  Test for convergence   ***
C
       IF(ABS(PNEW-PM).GT.0.2)THEN
        PM=PNEW
        NP=NP+1
        IF(NP.GT.1000.OR.PM.LT.400.0)THEN
         PMIN=400.0
         IFL=0
         GOTO 900
        END IF
        GOTO 100
       ELSE
        CATFAC=0.5*(1.+1./b)
        CAT=CAPEM-CAPEA+CKCD*RAT*CATFAC*(CAPEMS-CAPEM)
C        CAT=CKCD*RAT*CATFAC*(CAPEMS-CAPEM)
        PMIN=PSL*EXP(-CAT/(287.04*TVAV))
       END IF
  900       CONTINUE
       FAC=MAX(0.0,(CAPEMS-CAPEM))
       VMAX=VREDUC*SQRT(CKCD*RAT*FAC)
C
C   ***  Renormalize sounding arrays   ***
C
       DO 910 I=1,N
        R(I)=R(I)*1000.0
        T(I)=T(I)-273.15
  910       CONTINUE
C
       RETURN
       END
C
      SUBROUTINE CAPE(TP,RP,PP,T,R,P,ND,N,SIG,CAPED,TO,IFLAG)
C
C     This subroutine calculates the CAPE of a parcel with pressure PP (mb),
C       temperature TP (K) and mixing ratio RP (gm/gm) given a sounding
C       of temperature (T in K) and mixing ratio (R in gm/gm) as a function
C       of pressure (P in mb). ND is the dimension of the arrays T,R and P,
C       while N is the actual number of points in the sounding. CAPED is
C       the calculated value of CAPE and TO is the temperature at the
C       level of neutral buoyancy.  IFLAG is a flag
C       integer. If IFLAG = 1, routine is successful; if it is 0, routine did
C       not run owing to improper sounding (e.g.no water vapor at parcel level).
C       IFLAG=2 indicates that routine did not converge.
C
      REAL T(ND),R(ND),P(ND),TVRDIF(100)
      REAL NA
C
C   ***   Default values   ***
C
      CAPED=0.0
      TO=T(1)
      IFLAG=1
C
C   ***   Check that sounding is suitable    ***
C
      IF(RP.LT.1.0E-6.OR.TP.LT.200.0)THEN
       IFLAG=0
       RETURN
      END IF
C
C   ***   Assign values of thermodynamic constants     ***
C
      CPD=1005.7
      CPV=1870.0
C      CL=4190.0
      CL=2500.0
      CPVMCL=CPV-CL
      RV=461.5
      RD=287.04
      EPS=RD/RV
      ALV0=2.501E6
C
C   ***  Define various parcel quantities, including reversible   ***
C   ***                       entropy, S.                         ***
C
      TPC=TP-273.15
      ESP=6.112*EXP(17.67*TPC/(243.5+TPC))
      EVP=RP*PP/(EPS+RP)
      RH=EVP/ESP
      ALV=ALV0-CPVMCL*TPC
      S=(CPD+RP*CL)*LOG(TP)-RD*LOG(PP-EVP)+
     1   ALV*RP/TP-RP*RV*LOG(RH)
C
C   ***  Find lifted condensation pressure, PLCL   ***
C
       CHI=TP/(1669.0-122.0*RH-TP)
       PLCL=PP*(RH**CHI)
C
C   ***  Begin updraft loop   ***
C
       NCMAX=0
        DO J=1,N
         TVRDIF(J)=0.0
        END DO
       DO 200 J=2,N
C
C    ***   Don't bother lifting parcel above 60 mb    ***
C
        IF(P(J).LT.59.0)GOTO 200
C
C    ***  Parcel quantities below lifted condensation level   ***
C
        IF(P(J).GE.PLCL)THEN
         TG=TP*(P(J)/PP)**(RD/CPD)
         RG=RP
C
C   ***   Calculate buoyancy   ***
C
         TLVR=TG*(1.+RG/EPS)/(1.+RG)
         TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J))
        ELSE
C
C   ***  Parcel quantities above lifted condensation level  ***
C
         TG=T(J)
         TJC=T(J)-273.15
         ES=6.112*EXP(17.67*TJC/(243.5+TJC))
         RG=EPS*ES/(P(J)-ES)
C
C   ***  Iteratively calculate lifted parcel temperature and mixing   ***
C   ***                ratio for reversible ascent                    ***
C
         NC=0
  120         CONTINUE
         NC=NC+1
C
C   ***  Calculate estimates of the rates of change of the entropy    ***
C   ***           with temperature at constant pressure               ***
C
         ALV=ALV0-CPVMCL*(TG-273.15)
         SL=(CPD+RP*CL+ALV*ALV*RG/(RV*TG*TG))/TG
         EM=RG*P(J)/(EPS+RG)
         SG=(CPD+RP*CL)*LOG(TG)-RD*LOG(P(J)-EM)+
     1      ALV*RG/TG
         IF(NC.LT.3)THEN
          AP=0.3
         ELSE
          AP=1.0
         END IF
         TGNEW=TG+AP*(S-SG)/SL
C
C   ***   Test for convergence   ***
C
         IF(ABS(TGNEW-TG).GT.0.01)THEN
          TG=TGNEW
          TC=TG-273.15
          ENEW=6.112*EXP(17.67*TC/(243.5+TC))
C
C   ***   Bail out if things get out of hand   ***
C
          IF(NC.GT.500.OR.ENEW.GT.(P(J)-1.0))THEN
            IFLAG=2
            RETURN
          END IF
          RG=EPS*ENEW/(P(J)-ENEW)
          GOTO 120
         END IF
         NCMAX=MAX(NC,NCMAX)
C
C   *** Calculate buoyancy   ***
C
          RMEAN=SIG*RG+(1.-SIG)*RP
         TLVR=TG*(1.+RG/EPS)/(1.+RMEAN)
         TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J))
        END IF
  200       CONTINUE
C
C  ***  Begin loop to find NA, PA, and CAPE from reversible ascent ***
C
       NA=0.0
       PA=0.0
C
C   ***  Find maximum level of positive buoyancy, INB    ***
C
       INB=1
       DO 550 J=N,1,-1
        IF(TVRDIF(J).GT.0.0)INB=MAX(INB,J)
  550       CONTINUE
       IF(INB.EQ.1)RETURN
C
C   ***  Find positive and negative areas and CAPE  ***
C
       IF(INB.GT.1)THEN
        DO 600 J=2,INB
         TVM=0.5*(TVRDIF(J)+TVRDIF(J-1))
         PMA=0.5*(P(J)+P(J-1))
         IF(TVM.LE.0.0)THEN
          NA=NA-RD*TVM*(P(J-1)-P(J))/PMA
         ELSE
          PA=PA+RD*TVM*(P(J-1)-P(J))/PMA
         END IF
  600        CONTINUE
C
C   ***   Find residual positive area above INB and TO  ***
C
       PAT=0.0
       TO=T(INB)
       IF(INB.LT.N)THEN
        PINB=(P(INB+1)*TVRDIF(INB)-P(INB)*TVRDIF(INB+1))/
     1   (TVRDIF(INB)-TVRDIF(INB+1))
        PAT=RD*TVRDIF(INB)*(P(INB)-PINB)/(P(INB)+PINB)
       TO=(T(INB)*(PINB-P(INB+1))+T(INB+1)*(P(INB)-PINB))/
     1    (P(INB)-P(INB+1))
       END IF
C
C   ***   Find CAPE  ***
C
        CAPED=PA+PAT-NA
        CAPED=MAX(CAPED,0.0)
       END IF
C
       RETURN
       END


text:
      SUBROUTINE order_input(extension, nl, bsurf, ierr)
      USE stel_kinds
      USE read_wout_mod
      USE normalize_data
      USE readin_data
      USE fmesh_quantities, ONLY: rmncf, zmnsf, mercierf
     1  , rmnsf, zmncf                                                        ! 110909 RS: Asymmetric input
      USE vparams, ONLY : mu0
      IMPLICIT NONE
!-----------------------------------------------
!   D u m m y   A r g u m e n t s
!-----------------------------------------------
      CHARACTER*(*), INTENT(IN) :: extension
      INTEGER, INTENT(IN) :: nl
      INTEGER, INTENT(OUT) :: ierr
      INTEGER, DIMENSION(nl), INTENT(IN):: bsurf
!-----------------------------------------------
!   L o c a l   P a r a m e t e r s
!-----------------------------------------------
      REAL(rprec), PARAMETER :: zero = 0
!-----------------------------------------------
!   L o c a l   V a r i a b l e s
!-----------------------------------------------
      INTEGER :: lk, i, j, k
      INTEGER, DIMENSION(:), ALLOCATABLE :: plist
      REAL(rprec) :: aspect_v, r0max_v, r0min_v, betaxis_v
      LOGICAL, ALLOCATABLE, DIMENSION(:):: lballoon
!-----------------------------------------------
!
!      Variables allocated in this routine (deallocated elsewhere)
!         hiota, hpres, hphip, xm_v, xn_v, rmncf, zmnsf, mercierf,
!         lmnsh, bmnch, list, bsupumnch, bsupvmnch
!
!      110909 RS: IF LASYM_V = TRUE, also allocate arrays for asymmetric input
!         lmnch, bmnsh, bsupumnsh, bsupvmnsh, rmnsf, zmncf
!
!      All half mesh quantities have first surface (js =1) equal to zero
!

      ! SAL - PPPL Modification to avoid bad allocations
       IF (.not. lwout_opened) THEN
          call read_wout_file('wout_'//trim(extension), ierr)
          if (ierr .ne. 0) then
             ierr = 0
             call read_wout_file('wout.'//trim(extension), ierr)
             if (ierr .ne. 0) then
                print *,' ierr = ', ierr,
     1          'error in read_wout_file called from COBRA'
             end if
          end if
       END IF
!      IF (.not. lwout_opened) THEN
!         CALL read_wout_file(extension, ierr)                                    !...   Read surface quantities and normalization data
!         IF (ierr .ne. 0) THEN
!            PRINT *,' ierr = ', ierr,
!     1              'error in read_wout_file called from COBRA'
!         END IF
!      END IF
      IF (ierr_vmec.ne.0 .or. ierr.ne.0) GOTO 1000

      nfp_v = nfp
      ns_cob = ns
      aspect_v = aspect
      r0max_v = rmax_surf
      r0min_v = rmin_surf
      betaxis_v = betaxis
      mnmax_v = mnmax
      lasym_v = lasym                                                         ! 111009: RS, Logical for Asymmetric input

!...   IOTA, PRES, PHIP are ALL on HALF-MESH!!!

       k=0
       IF (.not. ALLOCATED(hiota)) ALLOCATE(hiota(ns_cob),stat=k)
       if (k .ne. 0) stop 'Allocation error 1 in cobra order_input'
       IF (.not. ALLOCATED(hpres)) ALLOCATE(hpres(ns_cob),stat=k)
       if (k .ne. 0) stop 'Allocation error 1 in cobra order_input'
       IF (.not. ALLOCATED(hphip)) ALLOCATE(hphip(ns_cob),stat=k)
       if (k .ne. 0) stop 'Allocation error 1 in cobra order_input'
       IF (.not. ALLOCATED(mercierf)) ALLOCATE(mercierf(ns_cob),stat=k)
       if (k .ne. 0) stop 'Allocation error 1 in cobra order_input'
!       allocate (hiota(ns_cob), hpres(ns_cob), hphip(ns_cob),
!     1   mercierf(ns_cob), stat=k)
!       if (k .ne. 0) stop 'Allocation error 1 in cobra order_input'

      hiota = iotas(1:ns_cob)
      IF (version_ .gt. 6.0) THEN
          hpres = mu0*pres(1:ns_cob)                                          ! in VMEC versions > 6.00 pressure is given in pascals
      ELSE
         hpres = pres(1:ns_cob)                                               ! in VMEC versions  .LE. 6.00 pressure is p= mu_o * p(pascals)
      ENDIF
      hphip = phip(1:ns_cob)

      mercierf = Dmerc(1:ns_cob)

      r0 = (r0max_v+r0min_v)/2
      amin = r0/aspect_v
      beta0 = betaxis_v
      IF (beta0 .le. zero) THEN
          beta0 = EPSILON(beta0)
          hpres = beta0*( /(1 - (j - 1)/REAL(ns - 1,rprec), j=1,ns)/ )
      END IF
      b0_v = SQRT((2._dp/beta0)*(1.5_dp*hpres(2)-.5_dp*hpres(3)))

       k=0
       IF (.not. ALLOCATED(xm_v)) ALLOCATE(xm_v(mnmax_v),stat=k)
       if (k .ne. 0) stop 'Allocation error 2 in cobra order_input'
       IF (.not. ALLOCATED(xn_v)) ALLOCATE(xn_v(mnmax_v),stat=k)
       if (k .ne. 0) stop 'Allocation error 2 in cobra order_input'
       IF (.not. ALLOCATED(plist)) ALLOCATE(plist(ns_cob),stat=k)
       if (k .ne. 0) stop 'Allocation error 2 in cobra order_input'
!       allocate (xm_v(mnmax_v), xn_v(mnmax_v), plist(ns_cob), stat=k)
!       if (k .ne. 0) stop 'Allocation error 2 in cobra order_input'
      xm_v = xm(1:mnmax_v)
      xn_v = xn(1:mnmax_v)

      mndim = ns_cob*mnmax_v

!...   RMN, ZMN are on FULL-MESH; LMNS, BSUPUMN (c), BSUPVMN (c),
!      and BMNC are on HALF-MESH!!!

       k=0
       IF (.not. ALLOCATED(rmncf)) ALLOCATE(rmncf(mndim),stat=k)
       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
       IF (.not. ALLOCATED(zmnsf)) ALLOCATE(zmnsf(mndim),stat=k)
       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
       IF (.not. ALLOCATED(lmnsh)) ALLOCATE(lmnsh(mndim),stat=k)
       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
       IF (.not. ALLOCATED(bmnch)) ALLOCATE(bmnch(mndim),stat=k)
       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
       IF (.not. ALLOCATED(bsupumnch)) ALLOCATE(bsupumnch(mndim),stat=k)
       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
       IF (.not. ALLOCATED(bsupvmnch)) ALLOCATE(bsupvmnch(mndim),stat=k)
       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
!       allocate (rmncf(mndim), zmnsf(mndim), lmnsh(mndim),
!     1    bmnch(mndim),  bsupumnch(mndim), bsupvmnch(mndim),
!     2    stat=k)
!       if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'

      zmnsf=0; rmncf=0; lmnsh=0; bmnch=0; bsupvmnch=0; bsupumnch=0

      IF (lasym_v) THEN                                                      !  110909 RS: For asymmetric input
         k=0
         IF (.not. ALLOCATED(rmnsf)) ALLOCATE(rmnsf(mndim),stat=k)
         if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
         IF (.not. ALLOCATED(zmncf)) ALLOCATE(zmncf(mndim),stat=k)
         if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
         IF (.not. ALLOCATED(lmnch)) ALLOCATE(lmnch(mndim),stat=k)
         if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
         IF (.not. ALLOCATED(bmnsh)) ALLOCATE(bmnsh(mndim),stat=k)
         if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
         IF (.not. ALLOCATED(bsupumnsh)) ALLOCATE(bsupumnsh(mndim),
     1       stat=k)
         if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
         IF (.not. ALLOCATED(bsupvmnsh)) ALLOCATE(bsupvmnsh(mndim),
     1       stat=k)
         if (k .ne. 0) stop 'Allocation error 3 in cobra order_input'
!         ALLOCATE (rmnsf(mndim), zmncf(mndim), lmnch(mndim),
!     1    bmnsh(mndim),  bsupumnsh(mndim), bsupvmnsh(mndim), stat=k)
!         IF (k .ne. 0) STOP 'Allocation error 3 in cobra order_input'
         zmncf=0; rmnsf=0; lmnch=0; bmnsh=0; bsupvmnsh=0; bsupumnsh=0
      ENDIF

      DO k=1, ns_cob
         lk = (k-1)*mnmax_v
         rmncf(lk+1:lk+mnmax_v) = rmnc(1:mnmax_v,k)
         zmnsf(lk+1:lk+mnmax_v) = zmns(1:mnmax_v,k)
         bmnch(lk+1:lk+mnmax_v) = bmnc(1:mnmax_v,k)
         lmnsh(lk+1:lk+mnmax_v) = lmns(1:mnmax_v,k)
         bsupvmnch(lk+1:lk+mnmax_v) = bsupvmnc(1:mnmax_v,k)
         bsupumnch(lk+1:lk+mnmax_v) = bsupumnc(1:mnmax_v,k)
         IF (lasym_v) THEN                                                   !  110909 RS: For asymmetric input
           rmnsf(lk+1:lk+mnmax_v) = rmns(1:mnmax_v,k)
           zmncf(lk+1:lk+mnmax_v) = zmnc(1:mnmax_v,k)
           bmnsh(lk+1:lk+mnmax_v) = bmns(1:mnmax_v,k)
           lmnch(lk+1:lk+mnmax_v) = lmnc(1:mnmax_v,k)
           bsupvmnsh(lk+1:lk+mnmax_v) = bsupvmns(1:mnmax_v,k)
           bsupumnsh(lk+1:lk+mnmax_v) = bsupumns(1:mnmax_v,k)
         ENDIF

      ENDDO

!...   identify wanted surfaces

      ALLOCATE (lballoon(ns_cob), stat=k)
      lballoon=.false.
      DO i=1,nl
         lballoon(bsurf(i))=.true.
      ENDDO

!...   check for feasibility of ballooning calculation

      nlist=0
      DO i = 2, ns_cob - 1                                                  ! exclude Boundary (i=ns)
         IF (lballoon(i)) THEN
            nlist=nlist+1
            plist(nlist)=i                                                  ! i-th surface on full is between i-th and (i+1)-th on half
         ENDIF
      ENDDO

!...   store NLIST surfaces WHERE calculation can be done in LIST(1:NLIST) vector

      IF (.not.ALLOCATED(list)) allocate (list(nlist), stat=k)
      list(1:nlist) = plist(1:nlist)
      DEALLOCATE (plist,lballoon)
!1000  CALL read_wout_deallocate
1000  IF (ierr .eq. 0) THEN
         ierr = -ierr_vmec
      ELSE
         PRINT *,' ierr = ', ierr, ' writing in READ_WOUT_BOOZ'
      END IF

      END SUBROUTINE order_input

ext: f

max_stars_repo_name: joseluisvelasco/STELLOPT

max_issues_repo_name: joseluisvelasco/STELLOPT

max_forks_repo_name: joseluisvelasco/STELLOPT


text:
      SUBROUTINE DROTM(N,DX,INCX,DY,INCY,DPARAM)
*     .. Scalar Arguments ..
      INTEGER INCX,INCY,N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION DPARAM(5),DX(*),DY(*)
*     ..
*
*  Purpose
*  =======
*
*     APPLY THE MODIFIED GIVENS TRANSFORMATION, H, TO THE 2 BY N MATRIX
*
*     (DX**T) , WHERE **T INDICATES TRANSPOSE. THE ELEMENTS OF DX ARE IN
*     (DY**T)
*
*     DX(LX+I*INCX), I = 0 TO N-1, WHERE LX = 1 IF INCX .GE. 0, ELSE
*     LX = (-INCX)*N, AND SIMILARLY FOR SY USING LY AND INCY.
*     WITH DPARAM(1)=DFLAG, H HAS ONE OF THE FOLLOWING FORMS..
*
*     DFLAG=-1.D0     DFLAG=0.D0        DFLAG=1.D0     DFLAG=-2.D0
*
*       (DH11  DH12)    (1.D0  DH12)    (DH11  1.D0)    (1.D0  0.D0)
*     H=(          )    (          )    (          )    (          )
*       (DH21  DH22),   (DH21  1.D0),   (-1.D0 DH22),   (0.D0  1.D0).
*     SEE DROTMG FOR A DESCRIPTION OF DATA STORAGE IN DPARAM.
*
*  Arguments
*  =========
*
*  N      (input) INTEGER
*         number of elements in input vector(s)
*
*  DX     (input/output) DOUBLE PRECISION array, dimension N
*         double precision vector with N elements
*
*  INCX   (input) INTEGER
*         storage spacing between elements of DX
*
*  DY     (input/output) DOUBLE PRECISION array, dimension N
*         double precision vector with N elements
*
*  INCY   (input) INTEGER
*         storage spacing between elements of DY
*
*  DPARAM (input/output)  DOUBLE PRECISION array, dimension 5
*     DPARAM(1)=DFLAG
*     DPARAM(2)=DH11
*     DPARAM(3)=DH21
*     DPARAM(4)=DH12
*     DPARAM(5)=DH22
*
*  =====================================================================
*
*     .. Local Scalars ..
      DOUBLE PRECISION DFLAG,DH11,DH12,DH21,DH22,TWO,W,Z,ZERO
      INTEGER I,KX,KY,NSTEPS
*     ..
*     .. Data statements ..
      DATA ZERO,TWO/0.D0,2.D0/
*     ..
*
      DFLAG = DPARAM(1)
      IF (N.LE.0 .OR. (DFLAG+TWO.EQ.ZERO)) GO TO 140
      IF (.NOT. (INCX.EQ.INCY.AND.INCX.GT.0)) GO TO 70
*
      NSTEPS = N*INCX
      IF (DFLAG) 50,10,30
   10 CONTINUE
      DH12 = DPARAM(4)
      DH21 = DPARAM(3)
      DO 20 I = 1,NSTEPS,INCX
          W = DX(I)
          Z = DY(I)
          DX(I) = W + Z*DH12
          DY(I) = W*DH21 + Z
   20 CONTINUE
      GO TO 140
   30 CONTINUE
      DH11 = DPARAM(2)
      DH22 = DPARAM(5)
      DO 40 I = 1,NSTEPS,INCX
          W = DX(I)
          Z = DY(I)
          DX(I) = W*DH11 + Z
          DY(I) = -W + DH22*Z
   40 CONTINUE
      GO TO 140
   50 CONTINUE
      DH11 = DPARAM(2)
      DH12 = DPARAM(4)
      DH21 = DPARAM(3)
      DH22 = DPARAM(5)
      DO 60 I = 1,NSTEPS,INCX
          W = DX(I)
          Z = DY(I)
          DX(I) = W*DH11 + Z*DH12
          DY(I) = W*DH21 + Z*DH22
   60 CONTINUE
      GO TO 140
   70 CONTINUE
      KX = 1
      KY = 1
      IF (INCX.LT.0) KX = 1 + (1-N)*INCX
      IF (INCY.LT.0) KY = 1 + (1-N)*INCY
*
      IF (DFLAG) 120,80,100
   80 CONTINUE
      DH12 = DPARAM(4)
      DH21 = DPARAM(3)
      DO 90 I = 1,N
          W = DX(KX)
          Z = DY(KY)
          DX(KX) = W + Z*DH12
          DY(KY) = W*DH21 + Z
          KX = KX + INCX
          KY = KY + INCY
   90 CONTINUE
      GO TO 140
  100 CONTINUE
      DH11 = DPARAM(2)
      DH22 = DPARAM(5)
      DO 110 I = 1,N
          W = DX(KX)
          Z = DY(KY)
          DX(KX) = W*DH11 + Z
          DY(KY) = -W + DH22*Z
          KX = KX + INCX
          KY = KY + INCY
  110 CONTINUE
      GO TO 140
  120 CONTINUE
      DH11 = DPARAM(2)
      DH12 = DPARAM(4)
      DH21 = DPARAM(3)
      DH22 = DPARAM(5)
      DO 130 I = 1,N
          W = DX(KX)
          Z = DY(KY)
          DX(KX) = W*DH11 + Z*DH12
          DY(KY) = W*DH21 + Z*DH22
          KX = KX + INCX
          KY = KY + INCY
  130 CONTINUE
  140 CONTINUE
      RETURN
      END


text:
Amber graduated from UCD in 2004 with a BS in Human Development and minor in Psychology.  She currently resides in Sacramento and works for the California Mortgage Bankers Association in downtown Sac but is looking to apply to grad schools for a masters in Sports Psychology.  In her spare time she can be found playing soccer, mainly DARSL, playing most other sports, hanging out at the race track (check out http://www.420racing.net 420 Racing), or listening to her cousins band, http://www.notoriousrocks.com Notorious play in Sacramento.

While at Davis she was very active in the Davis College Republicans, serving as Membership Director for the 20032004 school year, and intramural sports.  She was also a member of the Undergraduate Psychology Association.

Although she no longer lives in Davis, she journeys to Davis several times a week to coach an U14 girls http://www.davisayso.org AYSO team.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
Retro Film Revival used to be held at Regal Cinemas Davis Holiday 6 Holiday 6 on Thursdays at 10:00PM. Tickets cost $5, making it a pretty good cheap date. They have a fun atmosphere before the showing. Sometimes there is trivia contest (with prizes!) and they also toss out free candy. Unfortunately, the film schedule is not on their website, so this is THE place to see it. You can find a list of upcoming movies on Fandango under the Holiday 6 listing.

Upcoming shows:

None scheduled.

Heres the list of what they played before:

Retro Film Revival:
8 Films with Tongue in Cheek

20051103 16:12:01 nbsp The were checking bags for illicit food and drink when I went to see Scarface Users/RyanCastellucci

  Open your bag. Lessee... knife... gun... camera... smoke bomb... ok, you pass... WAIT! whats this here, Extreme Milk Duds Challenge Milk Duds?! You are so busted.

  20060104 14:30:40 nbsp hey i went to that scarface show too and the person in front of me had tried to bring in some candy in her bag and they made her dump it.  I however secreted a beer into the movie in my pants. Users/NathanielJarrett
    thats hot
      I bet the fact that you can secrete beer in your pants makes you a real hit with the ladies.

20060107 08:45:47 nbsp Added the new schedule here.  Please tell your friends as this is a really cool schedule of films. Users/DavidPetrie

20060208 22:44:12 nbsp The Big Lebowski is Retro now? I feel old. Users/DavidReid

20060216 22:01:18 nbsp UA Arden Faire 6 is also doing the same line of movies, normally a week behind Davis: http://movies.yahoo.com/showtimes/theater?id5784 Users/TarZxf

20060619 18:41:04 nbsp Have they been getting actual film prints of these movies, or are they shown on those video projectors that usually show the ads and crap? Users/AlanSmithee
  20060620 01:40:31 nbsp Theyre real prints Users/ArlenAbraham


!
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!     *                                                               *
!     *                  copyright (c) 1998 by UCAR                   *
!     *                                                               *
!     *       University Corporation for Atmospheric Research         *
!     *                                                               *
!     *                      all rights reserved                      *
!     *                                                               *
!     *                          Spherepack                           *
!     *                                                               *
!     *       A Package of Fortran77 Subroutines and Programs         *
!     *                                                               *
!     *              for Modeling Geophysical Processes               *
!     *                                                               *
!     *                             by                                *
!     *                                                               *
!     *                  John Adams and Paul Swarztrauber             *
!     *                                                               *
!     *                             of                                *
!     *                                                               *
!     *         the National Center for Atmospheric Research          *
!     *                                                               *
!     *                Boulder, Colorado  (80307)  U.S.A.             *
!     *                                                               *
!     *                   which is sponsored by                       *
!     *                                                               *
!     *              the National Science Foundation                  *
!     *                                                               *
!     * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * *
!
!
!
!
!     11/96
!
!     a program for testing all vector analysis and synthesis subroutines
!
!     (1) first a scalar stream function and a scalar velocity potential function
!         are set in st, sv by restricting polys in x, y, z to the sphere surface
!
!     (2) the vector vield (v, w) is set by analytically differenting the scalar fields in (1)
!         using the standard formula relating a vector field and the stream and velocity
!         potential scalar fields in colatitude X longitude spherical coordinates
!
!          v = -1/sin(theta)*d(st)/dphi + d(sv)/dtheta
!
!          w =  1/sin(theta)*d(sv)/dphi + d(st)/dtheta
!
!     (3) a vector analysis is performed on (v, w)
!
!     (4) a vector synthesis is performed using coeffs from (3)
!
!     (5) the synthesized vector field from (4) is compared with the vector field from (2)
!
!     note:  vhaec, vhaes, vhagc, vhags, vhsec, vhses, vhsgc, vhsgs are all tested!
!
program test_all_vector_analysis_and_synthesis_routines

    use spherepack

    ! Explicit typing only
    implicit none

    ! Dictionary
    integer(ip), parameter              :: NLAT= 25, NLON= 19, NT = 2
    integer(ip), parameter              :: ITYP = 0
    real(wp), dimension(NLAT, NLAT, NT) :: br, bi, cr, ci
    real(wp), dimension(NLAT, NLON, NT) :: st, sv, v, w
    real(wp), dimension(NLAT)           :: gaussian_latitudes, gaussian_weights
    real(wp), allocatable               :: wavetable(:)
    real(wp)                            :: cosp, cost, dlat, dphi
    real(wp)                            :: dstdp, dstdt, dsvdp, dsvdt
    real(wp)                            :: dxdp, dxdt, dydp, dydt, dzdp, dzdt
    real(wp)                            :: err2v, err2w
    integer(ip)                         :: i, j, k, icase
    integer(ip)                         :: error_flag
    real(wp)                            :: phi, sinp, sint, theta
    real(wp)                            :: ve, we, x, y, z

    call name("Testing all vector analysis and synthesis procedures")

    ! Set dimension variables
    call iout(NLAT, "nlat")
    call iout(NLON, "nlon")
    call iout(NT, "  nt")

    ! Set equally spaced colatitude and longitude increments
    dphi = TWO_PI/NLON
    dlat = PI/(NLAT-1)

    ! Compute nlat-many gaussian latitudinal points
    call compute_gaussian_latitudes_and_weights(NLAT, gaussian_latitudes, gaussian_weights, error_flag)

    call name("gaqd")
    call iout(error_flag, " ier")
    call vecout(gaussian_latitudes, "thtg", NLAT)

    ! Test all analysis and synthesis subroutines
    do icase=1, 4

        call name("*******************************")

        ! Set scalar stream and velocity potential fields as polys in x, y, z
        ! and then set v, w from st, sv scalar fields
        do k=1, NT
            do j=1, NLON
                phi = real(j - 1, kind=wp) * dphi
                sinp = sin(phi)
                cosp = cos(phi)
                do i=1, NLAT
                    select case (icase)
                        case (0:2)
                            theta = real(i-1, kind=wp) * dlat
                        case default
                            theta = gaussian_latitudes(i)
                    end select
                    cost = cos(theta)
                    sint = sin(theta)
                    x = sint*cosp
                    y = sint*sinp
                    z = cost
                    dxdt = cost*cosp
                    dxdp = -sint*sinp
                    dydt = cost*sinp
                    dydp = sint*cosp
                    dzdt = -sint
                    dzdp = 0.0_wp
                    select case (k)
                        case (1)
                            st(i, j, k) = x*y
                            sv(i, j, k) = y*z
                            dstdt = x*dydt+y*dxdt
                            dstdp = x*dydp+y*dxdp
                            dsvdt = y*dzdt+z*dydt
                            dsvdp = y*dzdp+z*dydp
                            v(i, j, k) = -(cosp*dydp+sinp*dxdp) + dsvdt
                            w(i, j, k) = sinp*dzdp + cost*cosp + dstdt
                        case (2)
                            st(i, j, k) = x*z
                            sv(i, j, k) = x*y
                            dstdp = x*dzdp+z*dxdp
                            dstdt = x*dzdt+z*dxdt
                            dsvdp = x*dydp+y*dxdp
                            dsvdt = x*dydt+y*dxdt

                            ! v = -1/sin(theta)*d(st)/dphi + d(sv)/dtheta
                            v(i, j, k) = z*sinp + dsvdt

                            ! w =  1/sin(theta)*d(sv)/dphi + d(st)/dtheta
                            w(i, j, k) = cosp*dydp+ sinp*dxdp + dstdt
                    end select
                end do
            end do
        end do

        select case (icase)
            case (1)

                call name("testing vhaec and vhsec")
                call initialize_vhaec(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhaec")
                call iout(error_flag, "error_flag = ")
                call vhaec(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhaec")
                call iout(error_flag, "error_flag = ")

                !  Now synthesize v, w from br, bi, cr, ci and compare with original
                call initialize_vhsec(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhsec")
                call iout(error_flag, "error_flag = ")
                call vhsec(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhsec")
                call iout(error_flag, "error_flag = ")
            case (2)

                call name("testing vhaes and vhses")
                call initialize_vhaes(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhaes")
                call iout(error_flag, "error_flag = ")
                call vhaes(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhaes")
                call iout(error_flag, "error_flag = ")

                ! Now synthesize v, w from br, bi, cr, ci and compare with original
                call initialize_vhses(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhses")
                call iout(error_flag, "error_flag = ")
                call vhses(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhses")
                call iout(error_flag, "error_flag = ")
            case (3)

                call name("testing vhagc and vhsgc")
                call iout(NLAT, "nlat")
                call initialize_vhagc(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhagc")
                call iout(error_flag, "error_flag = ")
                call vhagc(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhagc")
                call iout(error_flag, "error_flag = ")

                ! Now synthesize v, w from br, bi, cr, ci and compare with original
                call initialize_vhsgc(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhsgc")
                call iout(error_flag, "error_flag = ")

                call vhsgc(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhsgc")
                call iout(error_flag, "error_flag = ")
            case (4)

                call name("testing vhags and vhsgs")
                call initialize_vhags(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhags")
                call iout(error_flag, "error_flag = ")
                call vhags(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhags")
                call iout(error_flag, "error_flag = ")

                ! Now synthesize v, w from br, bi, cr, ci and compare with original
                call initialize_vhsgs(NLAT, NLON, wavetable, error_flag)
                call name("initialize_vhsgs")
                call iout(error_flag, "error_flag = ")
                call vhsgs(NLAT, NLON, ITYP, NT, v, w, NLAT, NLON, br, bi, cr, ci, NLAT, &
                    NLAT, wavetable, error_flag)
                call name("vhsgs")
                call iout(error_flag, "error_flag = ")
        end select

        err2v = 0.0_wp
        err2w = 0.0_wp
        do k=1, NT
            do j=1, NLON
                phi = real(j - 1, kind=wp)*dphi
                sinp = sin(phi)
                cosp = cos(phi)
                do i=1, NLAT
                    select case (icase)
                        case (0:2)
                            theta = real(i-1, kind=wp) * dlat
                        case default
                            theta = gaussian_latitudes(i)
                    end select
                    cost = cos(theta)
                    sint = sin(theta)
                    x = sint*cosp
                    y = sint*sinp
                    z = cost
                    dxdt = cost*cosp
                    dxdp = -sint*sinp
                    dydt = cost*sinp
                    dydp = sint*cosp
                    dzdt = -sint
                    dzdp = 0.0
                    select case (k)
                        case (1)
                            st(i, j, k) = x*y
                            sv(i, j, k) = y*z
                            dstdt = x*dydt+y*dxdt
                            dstdp = x*dydp+y*dxdp
                            dsvdt = y*dzdt+z*dydt
                            dsvdp = y*dzdp+z*dydp
                            ve = -(cosp*dydp+sinp*dxdp) + dsvdt
                            we = sinp*dzdp + cost*cosp + dstdt
                        case (2)
                            st(i, j, k) = x*z
                            sv(i, j, k) = x*y
                            dstdp = x*dzdp+z*dxdp
                            dstdt = x*dzdt+z*dxdt
                            dsvdp =  x*dydp+y*dxdp
                            dsvdt = x*dydt+y*dxdt
                            ve = z*sinp + dsvdt
                            we = cosp*dydp+ sinp*dxdp + dstdt
                    end select
                    err2v = err2v + (v(i, j, k) - ve)**2
                    err2w = err2w + (w(i, j, k) - we)**2
                end do
            end do
        end do

        ! Set and print least squares error in v, w
        err2v = sqrt(err2v/(NT*NLAT*NLON))
        err2w = sqrt(err2w/(NT*NLAT*NLON))
        call vout(err2v, "errv")
        call vout(err2w, "errw")
    end do

    ! Release memory
    deallocate (wavetable)

end program test_all_vector_analysis_and_synthesis_routines


text:
      integer             MAXTOK
      parameter           (MAXTOK = 10)
      character*(MAXTOK)
*     intervening comment
*     intervening comment
     x   toklst
      end

max_stars_repo_name: kurtsansom/ftnchek


text:
I graduated from UC Davis Davis in the fall but Im still hanging around. Im in Sacramento Sac during the day and out and about (probably drunk) at NIght Time night.

My best friend is PamelaYoo and I have 7 hamsters (a fucking hamster farm) and one adopted cat.

Pam, Me, and Dana~ drunk on Pams 23rd BDay


text:
      module timedata

      integer ntspts, freq, iterat, varcod
      integer iblkts
      real*8  tolpt
      logical exts

      integer,  allocatable :: statptts(:,:)
      real*8,  allocatable :: ptts(:,:)
      real*8,  allocatable :: parptts(:,:)
      real*8,  allocatable :: varts(:,:)

      end module

      module pp_data

      integer numppnodes, ppfreq

      integer,  allocatable :: ppnodes(:,:)

      end module

c-----------------------------------------------------------------------
c allocate the arrays
c-----------------------------------------------------------------------

      subroutine sTD

      use timedata
      include "common.h"

      allocate (statptts(ntspts,2))
      allocate (ptts(ntspts,nsd))
      allocate (parptts(ntspts,nsd))
      allocate (varts(ntspts,ndof))

      return
      end
c-----------------------------------------------------------------------
c delete the arrays
c-----------------------------------------------------------------------

      subroutine dTD

      use timedata

      deallocate (ptts)
      deallocate (varts)

      return
      end


text: #redirect North Davis Elementary School


text:
      SUBROUTINE CDTFBS (DX,DY,IOBUF,FILEU,NROW)
C
C     CDTFBS IS A SPECIAL VERSION OF GFBS USED BY COMPLEX DETERMINATE
C     METHOD
C
C     DEFINITION OF INPUT PARAMETERS
C
C     FILEU = MATRIX CONTROL BLOCK FOR THE UPPER TRIANGLE U
C     DX    = THE LOAD VECTOR B
C     DY    = THE SOLUTION VECTOR X
C     IOBUF = THE INPUT BUFFER
C
      INTEGER          FILEU(7)  ,TYPEAR   ,CDP      ,PARM(4)  ,
     1                 IOBUF(1)
      DOUBLE PRECISION DX(1)     ,DY(1)    ,DTEMP
      COMMON /NAMES /  RD        ,RDREW    ,WRT      ,WRTREW   ,
     1                 REW       ,NOREW    ,EOFNRW   ,RSP      ,
     2                 RDP       ,CSP      ,CDP      ,SQR      ,
     3                 RECT      ,DIAG     ,LOWTRI   ,UPRTRI   ,
     4                 SYM       ,ROW      ,IDENTY
      COMMON /UNPAKX/  ITYPEX    ,IUNPAK   ,JUNPAK   ,INCRX
      DATA    PARM(3), PARM(4)   /4HDETF,4HBS  /
C
      INCRX  = 1
      ITYPEX = CDP
      TYPEAR = CDP
C
C     BEGIN BACKWARD PASS
C
      IOFF = FILEU(7) - 1
      PARM(2) = FILEU(1)
      CALL OPEN (*80,FILEU,IOBUF,RD)
      DO 70 I = 1,NROW
      IUNPAK = 0
      J  = NROW - I + 1
      JJ = J + J
      CALL BCKREC (FILEU)
      CALL UNPACK (*90,FILEU,DY)
      CALL BCKREC (FILEU)
      ISING = 0
      K = (JUNPAK-IUNPAK+1)*2
      JU = JUNPAK + JUNPAK
      GO TO 30
   10 ISING = 1
C
C     DIVIDE BY THE DIAGONAL
C
      DTEMP    = (DX(JJ)*DY(K-1)-DX(JJ-1)*DY(K))/(DY(K)**2+DY(K-1)**2)
      DX(JJ-1) = (DX(JJ-1)*DY(K-1)+DX(JJ)*DY(K))/(DY(K)**2+DY(K-1)**2)
      DX(JJ  ) = DTEMP
   20 K  = K  - 2
      JU = JU - 2
      JUNPAK = JUNPAK - 1
      IF (K .EQ. 0) GO TO 60
      IF (DY(K).EQ.0.D0 .AND. DY(K-1).EQ.0.D0) GO TO 20
   30 IF (JUNPAK-J) 50,10,40
   40 JK = (J-IOFF)*2
      DX(JK-1) = DX(JK-1) - DX(JU-1)*DY(K-1) + DX(JU  )*DY(K)
      DX(JK  ) = DX(JK  ) - DX(JU  )*DY(K-1) - DX(JU-1)*DY(K)
      GO TO 20
   50 CONTINUE
      DX(JU-1) = DX(JU-1) - DX(JJ-1)*DY(K-1) + DX(JJ  )*DY(K)
      DX(JU  ) = DX(JU  ) - DX(JJ  )*DY(K-1) - DX(JJ-1)*DY(K)
      GO TO 20
   60 IF (ISING .EQ. 0) GO TO 90
   70 CONTINUE
      CALL CLOSE (FILEU,REW)
      RETURN
C
   80 PARM(1) = -1
      GO TO 100
   90 PARM(1) = -5
  100 CALL MESAGE (PARM(1),PARM(2),PARM(3))
      RETURN
      END


text:
! { dg-do compile }
! { dg-options "-std=f2003" }
! Test the patch for PR25098, where passing a variable as an
! actual argument to a formal argument that is a procedure
! went undiagnosed.
!
! Based on contribution by Joost VandeVondele  <jv244@cam.ac.uk>
!
integer function y()
  y = 1
end
integer function z()
  z = 1
end

module m1
contains
  subroutine s1(f)
    interface
      function f()
        integer f
      end function f
    end interface
  end subroutine s1
  subroutine s2(x)
    integer :: x
  end subroutine
end module m1

  use m1
  external y
  interface
   function x()
     integer x
   end function x
  end interface

  integer :: i, y, z
  i=1
  call s1(i) ! { dg-error "Expected a procedure for argument" }
  call s1(w) ! { dg-error "used as actual argument" }
  call s1(x) ! explicit interface
  call s1(y) ! declared external
  call s1(z) ! { dg-error "Expected a procedure for argument" }
  call s2(x) ! { dg-error "Invalid procedure argument" }
contains
  integer function w()
    w = 1
  end function w
end

! { dg-final { cleanup-modules "m1" } }

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
! { dg-do compile }
! Tests the fix for PR24325 in which the lack of any declaration
! that foo is a function or even a procedure was not detected.
!
! Contributed by Jakub Jelinek <jakub@gcc.gnu.org>
!
  integer foo
  call test
contains
  subroutine test
    integer :: i
    i = foo () ! { dg-error "is not a function" }
  end subroutine test
end


text:
	SUBROUTINE GDTXEV ( rlvlsi, nlvlsi, icord, rlvlso, nlvlso,
     +			    iret )
C************************************************************************
C* GDTXEV								*
C*									*
C* This subroutine gets evenly spaced levels for interpolation		*
C*									*
C* GDTXEV ( RLVLSI, NLVLSI, ICORD, RLVLSO, NLVLSO, IRET )		*
C*									*
C* Input parameters:							*
C*	RLVLSI (NLVLSI)	REAL		input levels			*
C*	NLVLSI		INTEGER		# of input levels		*
C*	ICORD		INTEGER		coordinate type:		*
C*					  1 = linear			*
C*					  2 = log			*
C*					  3 = **RKAPPA			*
C*									*
C* Output parameters:							*
C*	RLVLSO (NLVLSO)	REAL		output levels			*
C*	NLVLSO		INTEGER		# of output levels		*
C*	IRET		INTEGER		return code			*
C*					  0 = normal return		*
C*					 -4 = invalid vert coord	*
C**									*
C* Log:									*
C* K. F. Brill/GSC       7/89  						*
C* K. Brill/GSC          8/89           Interpolate to exact levels	*
C* K. Brill/GSC          1/90           Set lower/upper lvls for int	*
C* K. Brill/NMC		 1/91		Check for equality between 	*
C*        				input and output levels		*
C* T.W.Barker		 9/91		Created from gdtxgd		*
C* K. Brill/NMC		11/92		Fix calls to ER_WMSG		*
C************************************************************************
	INCLUDE 	'GEMPRM.PRM'
C*
        REAL 		rlvlsi(*),rlvlso(*)
C------------------------------------------------------------------------
	iret = 0
C*
C*
	rkinv = 1. / RKAPPA
C
C*	  Compute the vertical increment for the output grid using the
C*        minimum vertical differential distance in the transformed
C*        coordinate.
C
	difmin = 1.e32
	DO k = 1, nlvlsi
	    IF ( icord .eq. 1 ) THEN
		qk = rlvlsi ( k )
	    ELSE IF ( icord .eq. 2 ) THEN
		IF ( rlvlsi ( k ) .gt. 0.0 ) THEN
		    qk = ALOG ( rlvlsi ( k ) )
		ELSE
		    iret = -7
		    CALL ER_WMSG ( 'GDTHGT', -4, ' ', ier )
		    RETURN
		END IF
	    ELSE IF ( icord .eq. 3 ) THEN
		IF ( rlvlsi ( k ) .gt. 0.0 ) THEN
		    qk = rlvlsi ( k ) ** RKAPPA
		ELSE
		    iret = -7
		    CALL ER_WMSG ( 'GDTHGT', -4, ' ', ier )
		    RETURN
		END IF
	    ELSE
		iret = -7
		CALL ER_WMSG ( 'GDTHGT', -4, ' ', ier )
		RETURN
	    END IF
	    IF ( k .gt. 1 ) THEN
		dif = ABS ( qk - qkm1 )
		IF ( dif .lt. difmin ) difmin = dif
	    END IF
	    qkm1 = qk
	END DO
C
C*      Compute the bounds in output coordinate system.
C
	ybot = rlvlsi(1)
	ytop = rlvlsi(nlvlsi)
        IF ( icord .eq. 1 ) THEN
	    q1 = ybot
	    q2 = ytop
        ELSE IF ( icord .eq. 2 ) THEN
	    IF ( ybot .gt. 0.0 .and. ytop .gt. 0.0 ) THEN
		q1 = ALOG ( ybot )
		q2 = ALOG ( ytop )
	    ELSE
		iret = -7
		CALL ER_WMSG ( 'GDTHGT', -4, ' ', ier )
		RETURN
	    END IF
	ELSE IF ( icord .eq. 3 ) THEN
	    IF ( ybot .gt. 0.0 .and. ytop .gt. 0.0 ) THEN
		q1 = ybot ** RKAPPA
		q2 = ytop ** RKAPPA
	    ELSE
		iret = -7
		CALL ER_WMSG ( 'GDTHGT', -4, ' ', ier )
		RETURN
	    END IF
        ELSE
	    iret = -7
	    CALL ER_WMSG ( 'GDTHGT', -4, ' ', ier )
	    RETURN
	END IF

C*      Determine the number of levels to which to interpolate.
C
	rnolv = ABS ( q1 - q2 ) / difmin
	nvout = INT ( rnolv + .5 ) + 1
	IF ( nvout .gt. LLMXLV ) nvout = LLMXLV
C*
 	nvm1 = nvout - 1
	dz = ( q2 - q1 ) / FLOAT ( nvm1 )
C
C*	  Create array of output levels.
C
	z = q1
	DO k = 2, nvm1
	   z = z + dz
	   IF ( icord .eq. 1 ) rlvlso ( k ) = z
           IF ( icord .eq. 2 ) rlvlso ( k ) = EXP ( z )
           IF ( icord .eq. 3 ) rlvlso ( k ) = z ** rkinv
	END DO
C
C*	Set lower and upper most levels for the interpolation.
C
	rlvlso ( 1 ) = rlvlsi ( 1 )
	rlvlso ( nvout ) = rlvlsi ( nlvlsi )
	nlvlso = nvout
C
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak