      PROGRAM SEMIE
      PARAMETER(NP=257,NL=33)
      DOUBLE PRECISION  DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +  F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,Z2(NP,NL),P(NP,NL),
     + Q(NP,NL),PSI(NP,NL),VORTTS(NP,NL),RNSTS(NP,NL),Z1(NP,NL),
     + OUTPUT(20),EPFRM,EPSRM,EPFFM,EPSFM,DAY,DAYS,RTAPE,
     + VORTGS(NP,NL),VGS(NP,NL),RTEST,KEOUT(20),TS(NP,NL)

C       
      INTEGER MX,NX,NZ,NXP1,NZM1,NXM1,LOUT(20),I,J,NT,KTSD,KTS,
     +NSTEPS,NOUT,MDUMM,KT,LF,LS,KOUNT,W0,KTSK,LKEOUT(20),
     +KET,NKEOUT,IQ,IK
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM
C      
      COMMON/PREV/Z2
      COMMON/INIT/P,Q,PSI
      COMMON/VORT2/VORTTS,RNSTS,VORTGS,
     :VGS,TS

      COMMON/ZZ/Z1

      COMMON/TAPE/RTAPE

      COMMON/UNIF/IQ

c      DIMENSION KEOUT(20),LKEOUT(20)
      PRINT*, ' ----------------------------------------------'
      PRINT*, ' THIS IS THE SEMI-GEOSTROPHIC EADY-TYPE PROBLEM'
      PRINT*, ' ----------------------------------------------'
      IQ=0
      IK=0
	RTAPE=1.
	PRINT*, ' RTAPE = ' , RTAPE
      RTEST=1.
      IF(RTEST.EQ.0) GOTO 985

            

c      OPEN(UNIT=15,FILE=
c     : '/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/vg3')
c      OPEN(UNIT=16,FILE=
c     : '/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/pg3')
c      OPEN(UNIT=17,FILE=
c     :'/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/vgg3')
c      OPEN(UNIT=18,FILE=
c     : '/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/tg3')
c      OPEN(UNIT=19,FILE=
c     : '/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/wg3')
c      OPEN(UNIT=20,FILE=
c     : '/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/rv3')
c      OPEN(UNIT=21,FILE=
c     :'/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/rvt3')
c      OPEN(UNIT=22,FILE=
c     :'/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/uag3')
c      OPEN(UNIT=23,FILE=
c     :'/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/rts3')
c      OPEN(UNIT=24,FILE=
c     :'/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/vts3')
c      OPEN(UNIT=25,FILE=
c     :'/export/tlaloc/tcyclone-01/swsthcri/sgmodel/dyster/results/ts3')



  985    CONTINUE
C        OPEN(UNIT=51,FILE='/home2/chris2/pg115',
C     +		FORM='FORMATTED')
C        OPEN(UNIT=52,FILE='/home2/chris2/tg115',
C     +		FORM='FORMATTED')
	
C
C  SOLUTION OF THE TWO DIMENSIONAL EADY PROBLEM
C  NX AND NSPECIFY THE NUMBER OF GRID POINTS IN THE HORIZONTAL AND
C  VERTICAL RESPECTIVELY NZ.LT.OR.EQ.40 AND NX=2**MX ,MX.LT.OR.EQ.7
C
C-----------------------------------------------------------------------
C
C  THE TIME INTERVAL IS TIMINC
C  FOR STABILITY WE REQUIRE THAT  TIMINC.LT.DELX/(SH*ZL*0.5)
      DO 914 J=1,NL
        DO 914 I=1,NP
            VORTTS(I,J)=0.
            RNSTS(I,J)=0.
 914  CONTINUE
      
      


	

	KTSK=0

C  SET UP PARAMETERS , INITIAL FIELDS AND VARIABLES
C     REWIND(77)
C     REWIND(78)
      CALL START
C  KOUNT IS USED TO COUNT THE NUMBER OF INTEGRATIONS
      KOUNT=0
      KTS=0
C
C  SPECIFY THE TOTAL INTEGRATION PERIOD (DAYS) OR NUMBER OF TIMESTEPS
C  (NSTEPS), AND OUTPUT TIMES AS AN ARRAY OUTPUT , TIMES IN DAYS .
C  EPSRM AND EPSFM ARE MEASURES OF ACCURACY  OF THE SOLUTION
C  IF NSTEPS IS SPECIFIED DAYS MUST BE SET TO ZERO
C


C SET INTEGRATION LENGTH:
      DAYS=0.




      NSTEPS=1
      print*, ' days = ' , days
      print*, ' nsteps = ' , nsteps
      NT=1
      KTSD=20
      KTSD=4

     
        KTSD=(DAYS*24.*60.*60.)/(TIMINC*(NL-1))

      EPFRM=2.0E-5
      EPSRM=5.0E-05
C      EPSRM=2.0E-08 IS THE USUAL VALUE , LARGE VALUE IS 5.0E-05
      EPFFM=1.0E-4
      EPSFM=1.0E-4
C
C  SPECIFY OUTPUT DAYS , OUTPUT(1) MUST BE SET TO ZERO
C
      OUTPUT(1)=0.
      OUTPUT(2)=1.
      OUTPUT(3)=2.
      OUTPUT(4)=3.
      OUTPUT(5)=4.
      OUTPUT(6)=5.
      OUTPUT(7)=6.
      OUTPUT(8)=6.
      DO 1 I=9,20
 1    OUTPUT(I) = 0.0


      KEOUT(1)=0.
      KEOUT(2)=1.
      KEOUT(3)=2.
      KEOUT(4)=3.
      KEOUT(5)=4.
      KEOUT(6)=5.
      KEOUT(7)=6.
      KEOUT(8)=3.5
      KEOUT(9)=4.
      KEOUT(10)=4.5
      KEOUT(11)=5.
      KEOUT(12)=5.5
      KEOUT(13)=6.
      DO 113 I=14,20
 113    KEOUT(I) = 0.0


C
      IF(DAYS.LT.0.001)GOTO 2
      NSTEPS=IFIX(DAYS*86400.0/TIMINC+1.5)
 2    LOUT(1)=1
      LKEOUT(1)=1
      DO 3 I=2,20
      LKEOUT(I)=IFIX(86400.0*KEOUT(I)/TIMINC+1.5)

 3    LOUT(I)=IFIX(86400.0*OUTPUT(I)/TIMINC+1.5)
      NOUT=1
      NKEOUT=1
      MDUMM=2
      KT=1
      KET=1
      DO 101 J=1,NZ
      DO 101 I=1,NXP1
      PSI(I,J)=0.
  101 Z2(I,J)=PSI(I,J)
C     WRITE(14)PSI
C     WRITE(14)PSI
C
C-----------------------------------------------------------------------
C
C  MARCHING LOOP
C
C  CALCULATE THE COEFICIENT OF PHIXX : U(X,Z)
 7    CALL UCALC
C
C
C   CALCULATE THE RIGHT HAND SIDE OF THE PHI EQUATION
       CALL RHSPHI
C
C  SOLVE THE P.D.E. FOR PHI    :
C      U(X,Z) D2PHI(X,Z)/DX**2 + D2PHI(X,Z)/DZ**2 = U(X,Z) F**2
C
      CALL SOLVEFG(EPFRM,EPFFM,LF)
C
C  CALCULATE THE RIGHT HAND SIDE OF THE PSI EQUATION
      CALL RHPSI(KOUNT)
C
C  SOLVE ELLIPTIC FOR PSI :
C RO(Z)D((1/RO(Z))DPSI(X,Z)/DZ)/D+ D(U(X,Z)DPSI(X,Z)/DX)/(DX*F**2)=...
C..     -(2.0*SH*D(VG)/DX*(1/F))
C
      DO 102 J=1,NZ
      DO 102 I=1,NXP1
  102 Z1(I,J)=PSI(I,J)
      CALL SOLVESG(EPSRM,EPSFM,LS,KOUNT)
      DO 103 J=1,NZ
      DO 103 I=1,NXP1
  103 Z2(I,J)=Z1(I,J)
C
C  CALCULATE Q
      CALL QCALC
C
C  CALCULATE THE VERTICAL VELOCITY
C
      CALL VERTVEL
C
C
           CALL FRICT(KOUNT)
C
      KOUNT=KOUNT+1
C      PRINT*, ' KOUNT = ' , KOUNT
      KTS=KTS+1
      IF(KOUNT.EQ.1) KTS=0
C
C-----------------------------------------------------------------------
C
C  OUTPUT IF REQUIRED
      IF(KOUNT.EQ.NSTEPS)GOTO 5

      IF(KOUNT.EQ.KET)GOTO 44
  66    IF(KOUNT.EQ.KT)GOTO 4

C
C
C  OUTPUT FIELDS AND STOP IF THE CONVERGENCE IS BAD
      IF(MAX0(LF,LS).GT.150) GO TO 5
C
C-----------------------------------------------------------------------
C
C  INTEGRATE FORWARD IN TIME
C
    6 CONTINUE
       DAY=(FLOAT(KOUNT-1)/86400.0)*TIMINC

             CALL ADVQFXDI(KOUNT,DAY)
               CALL ADVPHIZ(KOUNT)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C TIME SERIES ADDED TO
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C

c        print*, ' kts = ', kts
c        print*, ' ktsd = ', ktsd
c        print*, ' nt = ', nt
c        print*, ' nl = ', nl   
     
         IF((KTS.EQ.KTSD.OR.KTS.EQ.0).AND.NT.LE.NL) THEN
         CALL XCALC(KOUNT)
         PRINT*, ' NT BEFORE CALL TO INTVORT' , NT
c         print*, ' vortts before vortts ', vortts
         CALL INTVORT(NT)
c         print*, ' vortts after intvort ' , vortts
         NT=NT+1
	 KTSK=KTSK+1
         KTS=0
         ENDIF
C
C
C-----------------------------------------------------------------------
C
      GOTO 7
C
 4    NOUT=NOUT+1
      KT=LOUT(NOUT)

      MDUMM=1
 5    DAY=(FLOAT(KOUNT-1)/86400.0)*TIMINC
      WRITE(15,'(5E10.3)') DAY,RTAPE,DELX,DELZ,XA
      CALL OUT(KOUNT,DAY,MDUMM)
      CALL PVSUM2(DAY)

      IF(MDUMM.EQ.2)GOTO 900
      MDUMM=2
      GOTO 6

 44      NKEOUT=NKEOUT+1
      KET=LKEOUT(NKEOUT)
      DAY=(FLOAT(KOUNT-1)/86400.0)*TIMINC
      CALL PVSUM2(DAY)
      GOTO 66

C
 900   PRINT*, ' VORTICITY TIME SERIES AT SURFACE '
       WRITE(21,'(4E10.3)') TIMINC,KTSD,DELX,KTSK,RTAPE
       WRITE(23,'(4E10.3)') TIMINC,KTSD,DELX,KTSK,RTAPE
       WRITE(24,'(4E10.3)') TIMINC,KTSD,DELX,KTSK,RTAPE
       WRITE(25,'(4E10.3)') TIMINC,KTSD,DELX,KTSK,RTAPE

        CALL OUTP(VORTTS,DAY,8)
        PRINT*, ' P V       TIME SERIES IN MIDDLE '
        CALL OUTP(RNSTS,DAY,9)
        CALL OUTP(VGS,DAY,11)
        CALL OUTP(TS,DAY,12)

      STOP
      END
      SUBROUTINE START  
      INTEGER MX,NX,NZ,NXP1,NZM1,W0,NXM1,I,J,L,AN,
     +NP,NL,KB,JJ,JJJ,KE,J1,NOTROP
      PARAMETER(NP=257,NL=33)
      DOUBLE PRECISION  DELX,DELZ,DELXSQ,DELZSQ,CRHS,
     +  F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,
     +PZB(NP),PZT(NP),RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +QINT(NP,NL),QTR(NP,NL),AS(NP),
     +PVINT(NL),XDIST,ZTROP,ZTHICK,GTROP,CTROP,GTPAUSE,
     +GSTRAT,ZTT,ZC,CP,CV,C,THETAZ,ROE0,P0,PI,STABLTY,GAMMAM1,GAMMA,
     +Z,ZA,ZTEMP,XTEMP,QANOM2,XLxx,ZL2,XSIG,
     +qanom1,delta1,xl1,zl1,xsig1,zsig1,delta2,xlb,xlt,xsigb,xsigt,
     +thetpb,zan1,expz1,expx1,
     +XL2,XSIG2,ZL4,ZSIG2,QPERT1,QPERT2,QSUM,QAV,ROSUM,ZL3,
     +EXPX2,EXPZ2,EXPX,ROAV,THETP,EXPXB,PZBP,
     +SL,DPERT,WV,RR,TH,CTH,BA,X,TROPL,CX,QE(NP,NL),HTR(NL),FT,
     +FUN,ZAN2,THETPT,EXPXT,PZTP
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/QMOIST/QE,FT,HTR

      COMMON/QINIT/QINT
      COMMON/PVTEST/PVINT
      COMMON/SLOPE/SL
	COMMON/QFORCE/QTR
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC IF KE=1 THEN DO EADY PROBLEM                        C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      KE=0
      IF(KE.EQ.1) THEN
      PRINT*, ' EADY PROBLEM'
      WV=1.61
      PRINT*, ' WAVENUMBER = ' , WV
      DPERT=0.1
      PRINT*, ' PERTURBATION TEMPERATURE = ' , DPERT
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCC IF NOTROP=1 THEN INCLUDE TROPOPAUSE                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      NOTROP=0.
      IF(NOTROP.EQ.0) PRINT*, ' NO TROPOPAUSE'
CCC AN=1 MEANS ANALYTICAL TROPOPAUSE 
CCC AN=0 MEANS 1103 TROPOPAUSE

      AN=1

      IF(AN.EQ.1.AND.NOTROP.EQ.1) THEN
      PRINT*, ' WITH ANALYTICAL TROPOPAUSE'
      ENDIF
      IF(AN.EQ.0.AND.NOTROP.EQ.1) THEN
      PRINT*, ' WITH 1103-TYPE TROPOPAUSE'
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      IF KB=0 THEN BOUSSINESQ
      KB=0
      IF(KB.EQ.0) THEN
      PRINT*, ' BOUSSINESQ'
      ELSE
      PRINT*, ' NON-BOUSSINESQ'
      ENDIF

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C     IF W0 = 0 LID BOUNDARY CONDITION IS W = 0
      W0=0
      IF(W0.EQ.1) THEN
      PRINT*, ' SLOPING LID'
      ELSE
      PRINT*, ' W = 0 AT LID'
      ENDIF


C      HEIGHT OF LID ZL
      ZL=10000.

C     HORIZONTAL LENGTH OF DOMAIN XDIST
      XDIST=5.E6


C DETERMINES NUMBER OF TIME STEPS
      IF(NP.EQ.257) THEN
      MX=8
      TIMINC=450.
      TIMINC=225.
      ELSE
      MX=7
      TIMINC=900.

      ENDIF

      NZ=NL
      DELZ=ZL/(NL-1)


CCCC HEATING PARAM. IGNORE!
      HTR(1)=1.
      DO 185 J=2,NL
      HTR(J)=HTR(1)
 185  CONTINUE
      PRINT*, ' HTR IN HEATING PARAM = ' , HTR(1)
C     HTR =1 THEN DRY IF SMALLER THEN MOIST

      DO 15 J=1,NL
            ZTT=(J-1)*DELZ
            HTR(J)=HTR(1)*EXP((ZTT/5000.)**2)
            IF(HTR(J).GT.1.) HTR(J)=1.
 15   CONTINUE


C
C      CONSTANTS
C
C  SPECIFIC HEAT OF DRY AIR AT CONSTANT PRESSURE
      CP=1004.0
C
C  SPECIFIC HEAT OF DRY AIR AT CONSTANT VOLUME
      CV=717.0
C
C  RATIO OF SPECIFIC HEATS
      GAMMA=CP/CV
C   VERTICAL - SHEAR
      SH=3.E-03
C
C  HORIZONTAL DIFFUSION PARAMETER EP
      EP=4.
C  GRAVITATIONAL ACCELERATION
      G=9.8066
C
C  CORIOLIS PARAMETER
c          F=0.0001
         F=0.000059
C
C  STANDARD TEMPERATURE (SURFACE)
      THETA0=288.15
C
C  RATIO OF G TO A STANDARD TEMPERATURE
      C=G/THETA0
C
C VERTICAL TEMPERATURE GRADIENT
        THETAZ= (4.6E-3)*100.
c        THETAZ=(THETA0*1.E-4/G)
c       THETAZ=(THETA0*f/G)
      print*, ' thetaz = ', thetaz
C
C  DENSITY OF AIR AT ZERO ALTITUDE
      ROE0=1.225
C
C  SURFACE PRESSURE
      P0=1.01325E5
C
C
      PI=3.1415926535898

CCC EADY STUFF  CCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF(KE.EQ.1) THEN
        W0=0
        ZL=10000.
        RR=ZL*SQRT(C*THETAZ)/F
        PRINT*, ' ROSSBY RADIUS = ' ,RR
        XDIST=2.*PI*RR/WV
        PRINT*, ' XDIST = ' , XDIST
      ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C
      XA=-XDIST/2.0
      XB=XA+XDIST
C
C  HORIZONTAL GRID LENGTH
      DELX=XDIST/2**MX
C
C  VERTICAL GRID LENGTH
      DELZ=ZL/(NZ-1)
C
      NXP1=2**MX+1
C
C    FOR STABILITY WE REQUIRE THAT
C    SH.LT.(2*DELX)/(TIMINC*ZL)
C    STABILITY PARAMETER
      STABLTY=SH*ZL*0.5*TIMINC/DELX
C-----------------------------------------------------------------------
C
C  SET UP PARAMETERS
      CALL CONST
C
C-----------------------------------------------------------------------
C
C  PSEUDO DENSITY DISTRIBUTION : RO(Z)
C
      GAMMAM1=GAMMA-1.0
      ZA=P0*GAMMA/(ROE0*G*GAMMAM1)
C
      DO 1 J=1,NZ
      Z=1.0-DELZ*(J-1)/ZA
      RO(J)=ROE0*DEXP(DLOG(Z)/GAMMAM1)
      IF(KB.EQ.0)RO(J)=ROE0
C
C  ROZ(Z) = (1/RO(Z))*DRO(Z)/D
C
      ROZ(J)=-1.0/(ZA*GAMMAM1*Z)
      IF(KB.EQ.0)ROZ(J)=0.
 1    CONTINUE
C
      IF(KE.EQ.1) GOTO 400

      DO 902 J=1,NZ
      ZTEMP=(J-1)*DELZ
C ZTEMP IS HEIGHT OF CURRENT LEVEL
      DO 902 I=1,NXP1
       P(I,J)=1000.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      XTEMP=(I-1)*DELX
C XTEMP IS HORIZONTAL DISTANCE IN SG SPACE OF CURRENT LOCATION
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  PV ANOMALY                                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C QANOM1 IS PVA ANOMALY VALUE
C DELTA1 DETERMINES VERTICAL TILT OF ANOMALY
C XL1 IS HORIZONTAL POSITION OF ANOMALY  IN METRES
C ZL1 IS VERTICAL POSITION OF ANOMALY IN METRES
C XSIG1 DETERMINES HORIZONTAL SCALE OF ANOMALY IN METRES
C ZSIG1 DETERMINES VERTICAL SCALE OF ANOMALY  IN METRES
      QANOM1=0.5E-6
      DELTA1=-0.
      XL1=XDIST/2.
      ZL1=4000.
C     ZL1=4.*DELZ
      XSIG1=200000.
      ZSIG1=1500.
C     IF(ZTEMP.GT.8000) ZSIG1=250.
      ZAN1 = ZL1-ZTEMP
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CAN INCLUDE A SECOND PV ANOMALY IF REQUIRED:
      QANOM2=0.E-6
C  IN KEYSER AND ZEHNDER PV ANOMALIES WRITTEN IN THESE TERMS?
      XL2=XL1+100000.
      XSIG2=200000.
      DELTA2=0.
      ZL2=7000.
      ZSIG2=1500.
      ZAN2 = ZL2-ZTEMP
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      EXPX1=((XTEMP-XL1-DELTA1*ZAN1)**2)/(XSIG1*XSIG1)
      EXPX2=((XTEMP-XL2-DELTA2*ZAN2)**2)/(XSIG2*XSIG2)
      EXPZ1=((ZTEMP-ZL1)**2)/(ZSIG1*ZSIG1)
      EXPZ2=((ZTEMP-ZL2)**2)/(ZSIG2*ZSIG2)
      
C DEFINE PV ANOMALY AS A WAVE IF REQUIRED 
C      LAMBDA=1000000.
C      TILT=1.
C      WAVE=SIN((2.*PI*XTEMP/LAMBDA)+TILT*PI*ZTEMP/ZL)
C      WAVE=1.
C      QPERT1=WAVE*QANOM1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC     



CCCCCCC   LOCALISED WIND ANOMALY IF REQUIRED CCCCCCC
C       XFUN=1.
C      XFUN=(XTEMP-XL)*SQRT(3.)/1000000.
C      FUN=1./(1+(XFUN**2))
C      FUN=FUN*FUN-4.*XFUN*XFUN*FUN*FUN*FUN
C      QPERT1=FUN*QANOM1*DEXP(-EXPZ)
C      IF(ZTEMP.EQ.4000.) THEN
C      PRINT*, ' XTEMP = ' , XTEMP, ' FUN = ' , FUN
C      ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      QPERT1=QANOM1*DEXP(-EXPX1)*DEXP(-EXPZ1)
      QPERT2=QANOM2*DEXP(-EXPX2)*DEXP(-EXPZ2)

C     IF(ZTEMP.GT.ZL2) QPERT1=0.
C
C PVINT(J) DEFINES BACKGROUND PV VALUE IN PV UNITS
       PVINT(J)=0.25

C Q IS ACTUTAL INITIAL PV DISTRIBUTION
       Q(I,J)=(PVINT(J))*1.E-6+QPERT1+QPERT2


C QTR AND QE USED IF DIABATIC HEATING INCLUDED
       QTR(I,J)=0.
       QE(I,J)=HTR(J)*PVINT(J)*1.E-6  
     
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC  TROPOPAUSE STUFF!                  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF(NOTROP.EQ.0) GOTO 943
      IF(AN.EQ.1) THEN
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ANALYTICAL TROPOPAUSE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        TROPL=0.5*DELZ
        TROPL=500.
        XSIG=250000.
        EXPX=(XTEMP-XL)*(XTEMP-XL)/(XSIG*XSIG)

        ZTROP=8000.-2000.*DEXP(-EXPX)
        ZTROP=8000.
        CX=-3.25*PI/(ATAN2(-ZTROP,TROPL))
        ZTT=(J-1)*DELZ
        ZL3=ZTT-ZTROP
        PVINT(J)=((ATAN2(ZL3,TROPL)*CX/3.14159)+3.5)
        Q(I,J)=PVINT(J)*1.E-6+QPERT1+QPERT2

        QE(I,J)=HTR(J)*PVINT(J)*1.E-6
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      ELSE


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  DO ANALYTICAL PV STRUCTURE LINEARLY INCREASING C
C  WITH HEIGHT BUT WITH DIFFERENT GRADIENTS IN    C
C  TROPOSPHERE TROPOPAUSE AND STRATOSPHERE        C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      XL=XDIST/2.
      XSIG=350000.
      EXPX=(XTEMP-XL)*(XTEMP-XL)/(XSIG*XSIG)

      ZTROP=10000.-1000.*DEXP(-EXPX)
      ZTROP=8000.

      ZTHICK=3000.

      GTROP=0.00E-3

      CTROP=0.25

      GTPAUSE=5.75/3000.

      GSTRAT=0.0E-3
      DO 541 J1=1,NL
         ZTT=(J1-1)*DELZ
         IF(ZTT.GT.ZTROP) GOTO 111
         PVINT(J1)=GTROP*ZTT+CTROP
 541    CONTINUE
 111    DO 121 JJ=J1,NL
         ZTT=(JJ-1)*DELZ-ZTROP
         ZC=(JJ-1)*DELZ
         IF(ZC.GT.(ZTROP+ZTHICK)) GOTO 131
         PVINT(JJ)=GTPAUSE*ZTT+PVINT(J1-1)
 121   CONTINUE
 131   DO 141 JJJ=JJ,NL
         ZTT=(JJJ-1)*DELZ-ZTROP-ZTHICK
         PVINT((JJJ))=GSTRAT*ZTT+PVINT(JJ-1)
 141   CONTINUE

        Q(I,J)=PVINT(J)*1.E-6+QPERT1+QPERT2

        QE(I,J)=HTR(J)*PVINT(J)*1.E-6

      ENDIF

 943  CONTINUE

  902 CONTINUE

      SL=F*F*SH/(C*RO(NZ)*PVINT(NZ)*1.E-6)
      PRINT*, ' SLOPE OF UPPER THETA SURFACE = ' , SL

      QSUM=0.
      QAV=0.
      ROSUM=0.
      DO 910 L=2,NZ-1
      DO 911 I=1,NX
      QSUM=QSUM+Q(I,L)*RO(L)
      QAV=QAV+Q(I,L)
      ROSUM=ROSUM+RO(L)
 911  CONTINUE
 910   CONTINUE
      ROAV=ROSUM/(NZ)
      QAV=QAV/NZ
      DO 9122 I=1,NX
      QSUM=QSUM+(Q(I,1)*RO(1)+Q(I,NZ)*RO(NZ))*0.5
 9122  CONTINUE
      THETZM=QSUM/((NZ-1)*F*NX)
      PRINT*, ' THETZM = ' , THETZM


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATES BOUNDARY TEMPERATURE AND ANOMALIES 
C IF REQUIRED
C THETPB= BOTTOM BOUNDARY  TEMP ANOMALY
C THETPT= TOP BOUNDARY TEMP ANOMALY
C XL* IS HORIZONTAL DISTANCE IN METRES OF ANOMALIES
C XSIG* DETERMINES HORIZONTAL SCALE OF ANOMALIES
C  C = G/THETAE0 NEEDED FOR CORRECT SCALING AND SET ABOVE
      THETPB=2.0
      THETPT=0.
      XLB=XDIST/2.+1000000.
      XLT=XL2-2000000.
      XSIGB=350000.
      XSIGT=350000.
      DO 903 I=1,NXP1
       PZB(I)=0.
       XTEMP=(I-1)*DELX

       EXPXB=(XTEMP-XLB)*(XTEMP-XLB)/(XSIGB*XSIGB)
       EXPXT=(XTEMP-XLT)*(XTEMP-XLT)/(XSIGT*XSIGT)
C 
C THESE ARE THE PERTURBATIONS
       PZBP=THETPB*DEXP(-EXPXB)*C
       PZTP=THETPT*DEXP(-EXPXT)*C

CCCCCCC   LOCALISED WIND ANOMALY CCCCCCCCCCCCCCCC
C      XFUN=(XTEMP-XL)*SQRT(3.)/1000000.
C      FUN=1./(1+(XFUN**2))
C      FUN=FUN*FUN-4.*XFUN*XFUN*FUN*FUN*FUN
C      PZBP=FUN*THETP*C
C      PRINT*, ' XTEMP = ' , XTEMP, ' FUN = ' , FUN
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


C THESE ARE THE ACTUAL BOUNDARY TEMPERTATURES INCLUDING
C MEAN AND ANOMALY:
      PZB(I)=PZBP+THETP
      PZT(I)=PZTP+PZB(1)+THETZM*ZL*C

  903 CONTINUE
  901 CONTINUE 



      GOTO 555



 400    CONTINUE



CCCCCCCCCCCCC EADY STUFF CCCCCCCC

      TH=TANH(WV*0.5)
      CTH=1./TH
      BA=(WV*0.5*TH-1.)/(WV*0.5*CTH-1.)
      IF(BA.LT.0.) THEN
      BA=-SQRT(-BA)
      ELSE
      BA=SQRT(BA)
      ENDIF
      DO 907 J=1,NZ
        Z=FLOAT(J-1)/FLOAT(NZM1)-0.5
        DO 907 I=1,NXP1
          X=FLOAT(I-1)*DELX/RR
          IF(BA.LT.0.) THEN
          FUN=SIN(WV*X)
          ELSE
          FUN=COS(WV*X)
          ENDIF
          P(I,J)=THETAZ*C*0.5*(FLOAT(J-1)*DELZ)**2+
     1           (DPERT*C*ZL/WV)*(COSH(WV*Z)*COS(WV*X)
     2           +BA*SINH(WV*Z)*FUN)
          Q(I,J)=F*THETAZ/RO(J)
          QE(I,J)=HTR(J)*Q(I,J)
 907  CONTINUE
      DO 908 I=1,NXP1
        X=FLOAT(I-1)*DELX/RR
        IF(BA.LT.0.) THEN
        FUN=SIN(WV*X)
        ELSE
        FUN=COS(WV*X)
        ENDIF
        PZB(I)=-DPERT*C*(SINH(0.5*WV)*COS(WV*X)-
     1         BA*COSH(0.5*WV)*FUN)
        PZT(I)=THETAZ*C*ZL+DPERT*C*(SINH(0.5*WV)*
     1         COS(WV*X)+BA*COSH(0.5*WV)*FUN)
 908  CONTINUE


 555  CONTINUE
C  AS AN INITIAL GUESS AT THE PSI EQUATION SET PSI(X,Z) TO ZERO
      DO 10 J=1,NZ
      DO 10 I=1,NXP1
 10   PSI(I,J)=0.0
C
C
C  OUTPUT MAIN PARAMETERS
C
      PRINT 100,NXP1
      PRINT 110,NZ
      PRINT 120,DELX
      PRINT 130,DELZ
      PRINT 140,TIMINC
      PRINT 150,SH
      PRINT 160,ZL
      PRINT 170,ZTROP
      PRINT 180,XDIST

C      PRINT 200,THETAZ
      PRINT 210,STABLTY
      IF(STABLTY.LE.1.0)GOTO 9
      PRINT 220
 9    CONTINUE
C9    PRINT 230
 100  FORMAT('THE NUMBER OF GRID POINTS IN THE HORIZONTAL IS' ,I4)
 110  FORMAT('THE NUMBER OF GRID POINTS IN THE VERTICAL IS ',I4)
 120  FORMAT('THE HORIZONTAL GRID LENGTH IS ',E12.5,' METERS')
 130  FORMAT('THE VERTICAL GRID LENGTH IS ',E12.5,' METERS')
 140  FORMAT(' THE TIME INCREMENT IS ',E12.5,' SECONDS')
 150  FORMAT(' THE SHEAR IS',E12.5, 'SECONDS-1')
 160  FORMAT(' THE HEIGHT OF THE LID IS ',E12.5,' METERS')
 170  FORMAT(' THE HEIGHT OF THE TROPOPAUSE IS ',E12.5,' METERS
     1')
 180  FORMAT(' THE HORIZONTAL DOMAIN IS OF LENGTH ',E12.5,' METERS')
 190  FORMAT(' THE INITIAL HORIZONTAL TEMPERATURE VARIATION IS ',E12
     1.5,' DEGREES')
 200  FORMAT(' THE INITIAL VERTICAL TEMPERATURE GRADIENT IS ',E12.5,
     1' DEGREES')
 210  FORMAT(' THE STABILITY PARAMETER : SH*ZL*0.5*TIMINC/DELX IS '
     1,F5.3)
 220  FORMAT(' WARNING ------ THE COMPUTATIONAL SCHEME IS UNSTABLE')
C230  FORMAT(1H1)
C
      PRINT*, ' NOW LETS GET RUNNING!!!!!!'
      RETURN
      END
      SUBROUTINE SOLVEFG(EPSRM,EPSFM,L)
C        IMPLICIT UNDEFINED(A-Z)      
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,L,J,I
      PARAMETER(NP=257,NL=33)


      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),DF(NP,NL),
     +DFB(NP),DFT(NP),RO(NL),ROZ(NL),F(NP,NL),U(NP,NL),PSI(NP,NL),
     +QINT(NP,NL),EPSRM,EPSFM,EPSR,EPSF

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/FI/DF
      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/F,U,PSI
      COMMON/QINIT/QINT
C
C-----------------------------------------------------------------------
C
C  THIS SUBROUTINE SOLVES THE P.D.E.
C      U(X,Z)*D2F(X,Z)/DX**2 + D2F(X,Z)/DZ**2 = RHS(X,Z)
C
C-----------------------------------------------------------------------
C
C      INITIAL F(I,J)
C
C      CALL SETF0
C      PRINT*, ' PRINTING F VALUES OUT AT TOP'
C       DO 987 I=1,NXP1
C         PRINT*, ' I = ' ,I , ' FTOP ' , F(I,NZ)
C 987    CONTINUE
C
C      CALCULATE V
C
      CALL SETV(U)
C
C      SET LAMBDAS FOR GAUSSIAN ELIMINATION
C
      CALL GAUSSA
C
      DO 1 L=1,200
C      IF(L.EQ.2) THEN
C      CALL SETF0
C      PRINT*, ' PHI AFTER CALL TO F0 FOR THE SECOND TIME'
C CC      CALL OUTP(F,99.00,6)
C      CALL GAUSSA
C      PRINT*, ' HERE WE ARE'
C      ENDIF
C
C      SET R.H.S.
C
      CALL SETRHS(EPSR)
C
C  SOLVE LINEARIZED EQUATION
C
      CALL SOLVEFL
C
C      ADD THE SOLUTION
C
      EPSF=0.0
      DO 10 J=1,NZ
      DO 10 I=1,NXP1
      EPSF=AMAX1(EPSF,ABS(DF(I,J)))
 10   F(I,J)=F(I,J)+DF(I,J)
C      PRINT *,'PHI AFTER AN ITERATION'
C      CALL OUTP(F,DAY,6)
C
C  TEST FOR ACCURAC
C        PRINT*, ' EPSR= ' , EPSR, ' EPSRM = ' , EPSRM
C        PRINT 100,L,EPSR,EPSF 
C       PRINT*, ' EPSF = ' , EPSF, ' EPSFM = ' , EPSFM
C
      IF((EPSF.LT.EPSFM).AND.(EPSR.LT.EPSRM))GOTO 2
 1    CONTINUE
C
 2    PRINT 100,L,EPSR,EPSF   
C      PRINT*, ' EPSR= ' , EPSR, ' EPSRM = ' , EPSRM
C C     PRINT*, ' EPSF = ' , EPSF, ' EPSFM = ' , EPSFM
 100  FORMAT(1X,' AFTER ',I4,' ITERATIONS THE RESIDUAL OF THE R.H.S. WAS
     1 ',E11.4,'  AND THE RESIDUAL OF THE FUNCTION WAS ',E11.4)
C
 9999   CONTINUE
      RETURN
      END
      SUBROUTINE SETF0
C        IMPLICIT UNDEFINED(A-Z)

      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,J,I,L
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),DFB(NP),
     +DFT(NP),RO(NL),ROZ(NL),F(NP,NL),U(NP,NL),PSI(NP,NL),RHS(NP,NL),
     + A(NP),R(NP),DSQ
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/F,U,PSI
      COMMON/RHS1/RHS
      COMMON/ARC/A,R
C
C-----------------------------------------------------------------------
C
C  THIS SUBROUTINE ADJUSTS TOP AND BOTTOM VALUES OF A TRIAL SOLUTION TO
C  MATCH THE NEUMANN BOUNDARY CONDITIONS
C DFB(I) AND DFT(I) ARE THE BOUNDARY DERIVATIVE DF(X,Z)/DFOR X AT THE
C  BOTTOM AND TOP OF THE REGION.
C
C-----------------------------------------------------------------------
C
      DSQ=2.0*(DELX/DELZ)**2
      DO 1 J=1,NZ,NZM1
C
C  CALCULATION OF THE R.H.S. AND A
C
      DO 11 I=2,NX
      R(I)=RHS(I,J)*DELXSQ/U(I,J)
 11   A(I)=-2.0-DSQ/U(I,J)
C
C  ADJUST FOR BOTTOM AND TOP
C
      IF(J.EQ.NZ)GOTO 13
      DO 12 I=2,NX
 12   R(I)=R(I)+(2.0+A(I))*(F(I,2)-DFB(I)*DELZ)
      GOTO 15
 13   DO 14 I=2,NX
 14   R(I)=R(I)+(2.0+A(I))*(F(I,NZM1)+DFT(I)*DELZ)
 15   CONTINUE
C  ADJUSTING FOR LEFT AND RIGHT
      R(2)=R(2)-F(1,J)
      R(NX)=R(NX)-F(NXP1,J)
C  GAUSSIAN ELIMINATION
      A(2)=1.0/A(2)
      R(2)=R(2)*A(2)
      DO 20 I=3,NX
      A(I)=1.0/(A(I)-A(I-1))
   20 R(I)=(R(I)-R(I-1))*A(I)
C   SOLUTION
      F(NX,J)=R(NX)
      L=NX-1
      DO 21 I=3,NX
      F(L,J)=R(L)-F(L+1,J)*A(L)
   21 L=L-1
 1    CONTINUE
CC      F(1,1)=F(2,1)
C      F(NXP1,1)=F(1,1)
C      F(1,NZ)=F(2,NZ)
C      F(NXP1,NZ)=F(1,NZ)
C
      RETURN
      END
      SUBROUTINE SETRHS(EPS)
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,J,I,W0
C
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),R(NP,NL),
     +DFB(NP),DFT(NP),RO(NL),ROZ(NL),F(NP,NL),U(NP,NL),PSI(NP,NL),
     +RHS(NP,NL),EPS
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS

      COMMON/FI/R
      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/F,U,PSI
      COMMON/RHS1/RHS
      EPS=0.0
      DO 11 J=2,NZM1
      R(1,J)=RHS(1,J)-U(1,J)*(F(NX,J)-2.0*F(1,J)+F(2,J))
     1       /DELXSQ-(F(1,J-1)-2.0*F(1,J)+F(1,J+1))/DELZSQ
      R(NXP1,J)=R(1,J)
      DO 11 I=2,NX
      R(I,J)=RHS(I,J)-U(I,J)*(F(I-1,J)-2.0*F(I,J)+F(I+1,J))/DELXSQ
     1-(F(I,J-1)-2.0*F(I,J)+F(I,J+1))/DELZSQ
      EPS=AMAX1(EPS,ABS(R(I,J)))
 11   CONTINUE
      R(1,1)=RHS(1,1)-U(1,1)*(F(NX,1)-2.0*F(1,1)+F(2,1))
     1       /DELXSQ-2.0*(F(1,2)-F(1,1)-DELZ*DFB(1))/DELZSQ
      R(1,NZ)=RHS(1,NZ)-U(1,NZ)*(F(NX,NZ)-2.0*F(1,NZ)+F(2,NZ))
     1        /DELXSQ-2.0*(F(1,NZM1)-F(1,NZ)+DELZ*DFT(1))/DELZSQ
      R(NXP1,1)=R(1,1)+2.0*(DFB(NXP1)-DFB(1))/DELZ
      R(NXP1,NZ)=R(1,NZ)+2.0*(DFT(1)-DFT(NXP1))/DELZ
      DO 2 I=2,NX
      R(I,1)=RHS(I,1)-U(I,1)*(F(I-1,1)-2.0*F(I,1)+F(I+1,1))/DELXSQ-2.0*(
     1F(I,2)-F(I,1)-DELZ*DFB(I))/DELZSQ
      R(I,NZ)=RHS(I,NZ)-U(I,NZ)*(F(I-1,NZ)-2.0*F(I,NZ)+F(I+1,NZ))/DELXSQ
     1-2.0*(F(I,NZM1)-F(I,NZ)+DELZ*DFT(I))/DELZSQ
 2    CONTINUE
      RETURN
      END
      SUBROUTINE SETV(A)
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,J,I
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),XL(NP,NL),
     +V(NL),DFB(NP),DFT(NP),RO(NL),ROZ(NL),F(NP,NL),U(NP,NL),PSI(NP,NL),
     +QINT(NP,NL),A(NP,NL),DELTA1,DELTA2
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS

c      DIMENSION A(NP,NL)

      COMMON/GAUSS/XL,V
      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/F,U,PSI
      COMMON/QINIT/QINT

      DO 1 J=1,NZ
      DELTA1=A(1,J)
      DELTA2=DELTA1
      DO 10 I=2,NX
      DELTA1=AMAX1(DELTA1,A(I,J))
 10   DELTA2=AMIN1(DELTA2,A(I,J))
 1    V(J)=0.5*(DELTA1+DELTA2)
      RETURN
      END
      SUBROUTINE SOLVEFL
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER  NP,NL,MX,NX,NZ,NXP1,NZM1,W0,I,K,J,L
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),XL(NP,NL),
     +V(NL),R(NP,NL)
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS

      COMMON/GAUSS/XL,V
      COMMON/FI/R
C
C  SOLUTION OF P.D.E.
C        V(Z) D2F(X,Z)/DX**2 + D2F(X,Z)/DZ**2 = R(X,Z)
C
C BY PERIODIC FOURIER TRANSFORM IN X AND GAUSSIAN ELIMINATION IN DIRECTION
C  BOUNDARY CONDITIONS
C      F(I,J) = F(NX+I,J)
C      F(I,2)-F(I,0) = F(I,NZ+1) - F(I,NZ-1) = 0
C      1.LE.I.LE.NX+1          X - DIECTION
C      1.LE.J.LE.N           - DIRECTION
C
C-----------------------------------------------------------------------
C
C  PREPARE R.H.S.
C     PRINT *,'XL AFTER GAUSS'
C     CALL OUTP(XL)
      DO 1 I=1,NXP1
      R(I,1)=0.5*CRHS*R(I,1)
      DO 11 J=2,NZM1
 11   R(I,J)=CRHS*R(I,J)
 1    R(I,NZ)=0.5*CRHS*R(I,NZ)
C     PRINT *,'R IS PREPARED RHS'
C     CALL OUTP(R)
C-----------------------------------------------------------------------
C
C      PERIODIC FOURIER TRANSFORM OF R.H.S.
C
C     PRINT*, ' JUST BEFORE CALL TO PFT-1 '
      CALL PFT(-1)
C     PRINT*, ' JUST AFTER CALL TO PFT-1 '
C      PRINT *,'R AFTER PERIODIC FOURIER TRANSFORM'
C     CALL OUTP(R)
C
C-----------------------------------------------------------------------
C      GAUSSIAN ELIMINATION
C      PREPARATION OF THE R.H.S.
C
      DO 2 K=1,NX
      R(K,1)=R(K,1)*XL(K,1)
      DO 21 J=2,NZ
 21   R(K,J)=(R(K,J)-R(K,J-1))*XL(K,J)
C      BACK SUBSTITUTION
      L=NZM1
      DO 22 J=2,NZ
      R(K,L)=R(K,L)-R(K,L+1)*XL(K,L)
 22   L=L-1
 2    CONTINUE
C      PRINT *,'R JUST BEFORE INVERSE TRANSFORM'
C     CALL OUTP(R)
C
C-----------------------------------------------------------------------
C
C      INVERSE PERIODIC FOURIER TRANSFORM
C
C     PRINT*,' JUST BEFORE CALL TO PFT+1'
      CALL PFT(+1)
CC    PRINT*,' JUST AFTER CALL TO PFT+1'
      DO 90 J=1,NZ
   90 R(NXP1,J)=R(1,J)
C     PRINT *,'R AFTER INVERSE TRANSFORM'
C     CALL OUTP(R)
C
C-----------------------------------------------------------------------
      RETURN
      END
      SUBROUTINE PFT(IBC)
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,N2,IBC,I0,I,J,MGPP,INC,
     +IFAX(20),IC,IS,II
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),R(NP,NL),
     +Y(NP,NL),X(44000),WORK(44000),TRIGS(800)
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/FI/R
      COMMON/AMP/Y
c     DIMENSION X(44000),WORK(44000),TRIGS(800),IFAX(20)
C
      N2=NX/2+1
      MGPP=NX+3
      INC=1
C     CALL FAX(IFAX,NX,3)
C     CALL FFTRIG(TRIGS,NX,3)
C
      IF(IBC.EQ.+1)GOTO 200
C
      I0=0
      DO 1 J=1,NZ
      DO 11 I=1,NX
      X(I+I0) = R(I,J)
 11   CONTINUE
      I0=I0+MGPP
 1    CONTINUE
C     PRINT*, ' JUST BEFORE FFT-1 CALL '
      CALL FFT99A(X,WORK,TRIGS,IFAX,INC,MGPP,NX,NZ,-1)
C     PRINT*, ' JUST AFTER FFT-1 CALL '
      I0=0
      DO 9 J=1,NZ
      IC=1
      IS=N2
      R(1,J)=X(I0+1)
      DO 111 I=4,NX,2
      IS=IS+1
      II=I-1
      IC=IC+1
      R(IS,J)=X(I+I0)
  111 R(IC,J)=X(I0+II)
      R(N2,J)=X(I0+NXP1)
      I0=I0+MGPP
 9    CONTINUE
C
      RETURN
C
 200  CONTINUE
      I0=0
      DO 2 J=1,NZ
      IS=N2
      IC=1
      X(1+I0)=R(1,J)
      X(2+I0)=0.
      DO 21 I=4,NX,2
      IS=IS+1
      IC=IC+1
      II=I-1
      X(I+I0)=R(IS,J)
 21   X(II+I0) = R(IC,J)
      X(NXP1+I0)=R(N2,J)
      I0=I0+MGPP
 2    CONTINUE
C     PRINT*, ' JUST BEFORE FFT+1 CALL '
      CALL FFT99A(X,WORK,TRIGS,IFAX,INC,MGPP,NX,NZ,+1)
C     PRINT*, ' JUST AFTER FFT+1 CALL '
      I0=0
      DO 115 J=1,NZ
      DO 223 I=1,NX
 223  R(I,J) = X(I+I0)
      I0=I0+MGPP
 115  CONTINUE
      RETURN
      END
      SUBROUTINE GAUSSA
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,K,J
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),XL(NP,NL),
     +V(NL)
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS

      COMMON/GAUSS/XL,V
      DO 1 K=1,NX
      XL(K,1)=-1.0-0.5*AS(K)*V(1)
      DO 11 J=2,NZM1
 11   XL(K,J)=-2.0-AS(K)*V(J)
      XL(K,NZ)=-1.0-0.5*AS(K)*V(NZ)
    1 IF(K.EQ.1)XL(K,NZ)=0.
C
C      CALCULATION OF LAMBDAS
C
      DO 2 K=1,NX
      XL(K,1)=1.0/XL(K,1)
      DO 21 J=2,NZ
 21   XL(K,J)=1.0/(XL(K,J)-XL(K,J-1))
 2    CONTINUE
      XL(1,NZ)=0.
      DO 23 J=1,NZ
   23 XL(NXP1,J)=XL(1,J)
      RETURN
      END
      SUBROUTINE CONST
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER  NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,NXH,K
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),F,THETA0,
     +G,TIMINC,XA,XB,C1,C2,SH,ZL,EP,THETZM
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM
      W0=W0*1
      NX=2**MX
      NXP1=NX+1
      NXM1=NX-1
      NZM1=NZ-1
      DELXSQ=DELX**2
      DELZSQ=DELZ**2
      CRHS=DELZSQ
      C1=4.0*(DELZ/DELX)**2
      NXH=NX/2+1
      C2=3.1415926535898D0/NX
      DO 1 K=1,NXP1
      IF(K-NXH)2,2,3
    2 AS(K)=C1*DSIN(C2*FLOAT(K-1))**2
      GOTO 1
    3 AS(K)=C1*DSIN(C2*FLOAT(K-NXH))**2
    1 CONTINUE

      RETURN
      END
      SUBROUTINE SOLVESG(EPSRM,EPSFM,L,KOUNT)
C        IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,
     +NXM1,NZM2,I,J,L,K,KOUNT
C
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F0,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,DF(NP,NL),DFB(NP),
     +DFT(NP),RO(NL),ROZ(NL),F(NP,NL),U(NP,NL),S(NP,NL),QINT(NP,NL),
     +Z(NP,NL),Z2(NP,NL),WTOP(NP),Z1(NP,NL),VGT(NP),CZ,CV,CT,
     +EPSR,EPSF,EPSFM,EPSRM,SL
C
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F0,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM
 

      COMMON/FI/DF
      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/F,U,S
      COMMON/QINIT/QINT
      COMMON/PSINT/Z
      COMMON/PREV/Z2
      COMMON/WUPPER/WTOP
      COMMON/ZZ/Z1
      COMMON/SLOPE/SL
c added (np) jan 29 2003
c      DIMENSION VGT(NP)
      NZM2=NZM1-1
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATE W AT TOP LEVEL ASSUMING      C
C ADIABATIC ASCENT                       C
      CALL QCALC
      CZ=0.5/DELZ
      CV = 1.0/(F0*2.0*DELX)
      CT=THETA0/G
      IF(W0.EQ.0) CT=0.
C
C          CALCULATE  THE  LONG-FRONT   VELOCITY   V
C          AT TOP  VGT(Z) = (DPHI(X,Z)/DX)/F0
C
C   WTOP IS F0*W*THZ/Q  WHICH EQUALS D(PSI)/DX
      VGT(1) = CV*(F(2,NZ) - F(NX,NZ))
      WTOP(1)=-VGT(1)*CT*SH*RO(NZ)*F0*F0/U(1,NZ)

C     WTOP(1)=-VGT(1)*SL

      DO 803 I=2,NX
      VGT(I) = CV*(F(I+1,NZ) - F(I-1,NZ))
      WTOP(I)=-VGT(I)*CT*SH*F0*RO(NZ)*F0/U(I,NZ)

C       WTOP(I)=-VGT(I)*SL

 803   CONTINUE
      VGT(NXP1)=VGT(1)
      WTOP(NXP1)=WTOP(1)
      S(1,NZ)=0.0
      S(NXP1,NZ)=0.
      Z2(1,NZ)=S(1,NZ)
      Z2(NXP1,NZ)=S(1,NZ)
      Z1(1,NZ)=S(1,NZ)
      Z1(NXP1,NZ)=S(1,NZ)
      DO 225 I=2,NX
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC   PSI AT TOP
      S(I,NZ)=0.5*(WTOP(I)+WTOP(I-1))*DELX+S(I-1,NZ)
      IF(KOUNT.EQ.0)  Z2(I,NZ)=S(I,NZ)
      Z1(I,NZ)=S(I,NZ)
 225  CONTINUE
      CALL UCALC
      DO 111 J=2,NZM1
      Z(NXP1,J)=U(NXP1,J)
      Z(1,J)=Z(NXP1,J)
      DO 111 I=2,NX
      Z(I,J)=U(I,J)
  111 CONTINUE
      DO 112 I=1,NXP1
      Z(I,1)=Z(I,2)
  112 Z(I,NZ)=Z(I,NZM1)
C
C-----------------------------------------------------------------------
C
C  THIS SUBROUTINE SOLVES THE P.D.E.
C  RO(Z)D((1/RO(Z))DS(X,Z)/DZ)/D+ D(U(X,Z)DS(X,Z)/DX)/DX = RHS(X,Z)
C
C-----------------------------------------------------------------------
C
C      CALCULATE V
C
      CALL SETV(U)
C
C      SET LAMBDAS FOR GAUSSIAN ELIMINATION
C
      CALL GAUSSB
      GOTO 98
C     DO 99 J=1,NZ
C     DO 99 I=1,NXP1
C  99 S(I,J)=2.*S(I,J)-Z2(I,J)
   98 CONTINUE
      K=2
C
      DO 1 L=1,200
C
C      SET R.H.S.
C
      CALL SETRHS2(EPSR,K)
C
C  SOLVE LINEARIZED EQUATION
C
      CALL SOLVESL
C
C      ADD THE SOLUTION
C
      EPSF=0.0
      DO 10 J=2,NZM1
      DO 10 I=1,NXP1
      EPSF=AMAX1(EPSF,ABS(DF(I,J)))
 10   S(I,J)=S(I,J)+DF(I,J)
C
C  TEST FOR ACCURACY
C
C     IF((EPSF.LT.EPSFM).AND.(EPSR.LT.EPSRM))GOTO 2
       IF(EPSR.LT.EPSRM.AND.L.GE.10)GOTO 2
C     IF(EPSR.LT.EPSRM)GOTO2
 1    CONTINUE
C
 2    PRINT 100,L,EPSR,EPSF
 100  FORMAT(1X,' AFTER ',I4,' ITERATIONS THE RESIDUAL OF THE R.H.S. WAS
     1 ',E11.4,'  AND THE RESIDUAL OF THE FUNCTION WAS ',E11.4)
C
      K=2
C      CALL SETRHS2(EPSR,K)
C      CALL SOLVESL
      EPSF=0.
C      DO 11 J=2,NZM1
C      DO 11 I=1,NXP1
C   11 EPSF=AMAX1(EPSF,ABS(DF(I,J)))
C      L=L+1
C      PRINT 100,L,EPSR,EPSF
      RETURN
      END
      SUBROUTINE SETRHS2(EPS,K)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,
     +K,J,I,M,I0,IJ,NXM1



      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F0,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,R(NP,NL),DFB(NP),
     +DFT(NP),RO(NL),ROZ(NL),F(NP,NL),U(NP,NL),S(NP,NL),QINT(NP,NL),
     +RHS(NP,NL),Z1(NP,NL),Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),
     +RQMT(NP,NL),RQT(NP,NL),RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),
     +R1(NP,NL),FF(NL),Z(NP,NL),V0,VJ,EPS,DEL4SQ,RZ

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F0,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/FI/R
      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/F,U,  S
      COMMON/QINIT/QINT
      COMMON/RHS1/RHS
      
      COMMON/PSINT/Z1
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
c      DIMENSION R1(NP,NL),FF(NL),Z
      GO TO (222,225)K
  225 CONTINUE
      DO 111 J=2,NZM1
      Z(NXP1,J)=U(NXP1,J)
      Z(1,J)=Z(NXP1,J)
      DO 111 I=2,NX
      Z(I,J)=U(I,J)
  111 CONTINUE
  222 CONTINUE
C   IF M=1 THEN UPRIGHT HEATING
      M=0
      IF(M.EQ.0)THEN
      DO 14 J=1,NZ
      DO 14 I=1,NXP1
   14 R1(I,J)=RHS(I,J)
      GOTO 13
      ENDIF
      DO 12 J=1,NZ
   12 FF(J)=0.1*3.4E10*DSIN(FLOAT(J-1)*3.14159D0/FLOAT(NZM1))
      CALL VERTVEL
      DO 10 I=1,NXP1
      I0=I
      IF(I.EQ.1.OR.I.EQ.NXP1)THEN
      V0=(F(2,1)-F(NX,1))*0.5/(DELX*1.0E-04)
      ELSE
      V0=(F(I+1,1)-F(I-1,1))*0.5/(DELX*1.0E-04)
      ENDIF
      DO 10 J=1,NZ
      Z1(I,J)=0.
      Z(I,J)=U(I,J)
      IF(I.EQ.1.OR.I.EQ.NXP1)THEN
      VJ=(F(2,J)-F(NX,J))*0.5/(DELX*1.0E-04)
      ELSE
      VJ=(F(I+1,J)-F(I-1,J))*0.5/(DELX*1.0E-04)
      ENDIF
      IJ=I0+INT((VJ-V0)/(1.0E-04*DELX))
      IF(IJ.GT.NP)IJ=IJ-(NP-1)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC BEWARE 256 
      IF(IJ.LT.1)IJ=IJ+256
      IF(W(I,2).GE.0.)THEN
      Z1(IJ,J)=W(I,2)*FF(J)/DELZ
      ENDIF
   10 CONTINUE
      DO 17 J=1,NZ
      DO 17 I=1,NXP1
      IF(I.EQ.1.OR.I.EQ.NXP1)THEN
      R1(I,J)=RHS(I,J)-(Z1(2,J)-Z1(NX,J))*0.5/DELX
      ELSE
      R1(I,J)=RHS(I,J)-(Z1(I+1,J)-Z1(I-1,J))*0.5/DELX
      ENDIF
   17 CONTINUE
   13 CONTINUE
      EPS=0.0
      DEL4SQ=4.0*DELXSQ
      DO 1 J=2,NZM1
      RZ=ROZ(J)*0.5/DELZ
      R(1,J)=R1(1,J)-(S(1,J-1)-2.0*S(1,J)+S(1,J+1))
     1       /DELZSQ+RZ*(S(1,J+1)-S(1,J-1))-((Z(2,J)+
     2       Z(1,J))*(S(2,J)-S(1,J))-(Z(1,J)+Z(NX,J))*
     3       (S(1,J)-S(NX,J)))*0.5/DELXSQ
      DO 11 I=2,NX
      R(I,J)=R1(I,J)-(S(I,J-1)-2.0*S(I,J)+S(I,J+1))/DELZSQ+RZ*(S(I,J+1)
     1   -S(I,J-1))
     2              -((Z(I+1,J)+Z(I,J))*(S(I+1,J)-S(I,J))-
     3    (Z(I,J)+Z(I-1,J))*(S(I,J)-S(I-1,J)))*0.5/DELXSQ
      EPS=AMAX1(EPS,ABS(R(I,J)))
 11   CONTINUE
 1    R(NXP1,J)=R(1,J)
      R(1,1)=0.0
      R(1,NZ)=0.0
      R(NXP1,1)=0.
      R(NXP1,NZ)=0.
      DO 2 I=2,NX
      R(I,NZ)=0.0
      R(I,1)=0.0
 2    CONTINUE
      RETURN
      END
      SUBROUTINE SOLVESL
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,
     +I,J,K,L

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +XL(NP,NL),V(NL),AL(NL),XNU(NP,NL),R(NP,NL)


      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/GAUSS/XL,V
      COMMON/GAUSS2/AL,XNU
      COMMON/FI/R
C
C  SOLUTION OF P.D.E.
C      V(Z) D2S(X,Z)/DX**2 - ROZ(Z) DS(X,Z)/D+ D2S(X,Z)/DZ**2 = R(X,Z)
C
C BY PERIODIC FOURIER TRANSFORM IN X AND GAUSSIAN ELIMINATION IN DIRECTION
C  BOUNDARY CONDITIONS
C      S(I,J) = S(NX+I,J)
C      S(I,1) = S(I,NZ) = 0.0 OR W TERM
C      1.LE.I.LE.NX+1          X - DIRECTION
C      1.LE.J.LE.N           - DIRECTION
C-----------------------------------------------------------------------
C
C  PREPARE R.H.S.
      DO 1 I=1,NX
      DO 1 J=2,NZM1
 1    R(I,J)=CRHS*R(I,J)
C-----------------------------------------------------------------------
C
C      PERIODIC FOURIER TRANSFORM OF R.H.S.
C
      CALL PFT(-1)
C
C-----------------------------------------------------------------------
C      GAUSSIAN ELIMINATION
C      PREPARATION OF THE R.H.S.
C
      DO 2 K=1,NX
      R(K,2)=R(K,2)*XL(K,2)
      DO 21 J=3,NZM1
 21   R(K,J)=(R(K,J)-(1.0+AL(J))*R(K,J-1))*XL(K,J)
C      BACK SUBSTITUTION
      L=NZ-2
      DO 22 J=3,NZM1
      R(K,L)=R(K,L)-XNU(K,L)*R(K,L+1)
 22   L=L-1
 2    CONTINUE
C
C-----------------------------------------------------------------------
C
C      INVERSE PERIODIC FOURIER TRANSFORM
C
      CALL PFT(+1)
      DO 90 J=1,NZ
   90 R(NXP1,J)=R(1,J)
C
C-----------------------------------------------------------------------
      RETURN
      END
      SUBROUTINE GAUSSB
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,J,K

      PARAMETER(NP=257,NL=33)
C        IMPLICIT UNDEFINED(A-Z)
      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +XL(NP,NL),V(NL),AL(NL),XNU(NP,NL),PZB(NP),PZT(NP),RO(NL),
     +ROZ(NL),DELZ2

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS

      COMMON/GAUSS/XL,V
      COMMON/GAUSS2/AL,XNU
      COMMON/INTS/PZB,PZT,RO,ROZ
      DELZ2=0.5*DELZ
      DO 10 J=2,NZM1
 10   AL(J)=ROZ(J)*DELZ2
C
      DO 1 K=1,NX
      DO 1 J=2,NZM1
 1    XL(K,J)=-2.0-AS(K)*V(J)
C
C      CALCULATION OF LAMBDAS AND NUS
C
      DO 2 K=1,NX
      XL(K,2)=1.0/XL(K,2)
      XNU(K,2)=(1.0-AL(2))*XL(K,2)
      DO 21 J=3,NZM1
      XL(K,J)=1.0/(XL(K,J)-(1.0+AL(J))*XNU(K,J-1))
 21   XNU(K,J)=XL(K,J)*(1.0-AL(J))
 2    CONTINUE
      RETURN
      END
      SUBROUTINE RHSPHI
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,I,J

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),RHS(NP,NL),F2

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/RHS1/RHS
      F2=F**2
      DO 1 J=1,NZ
      DO 1 I=1,NXP1
 1    RHS(I,J)=F2*Q(I,J)
      RETURN
      END
      SUBROUTINE RHPSI(KOUNT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +KOUNT,J,I

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),F1(NP,NL),
     +FT1(NP,NL),A(NP,NL),B(NP,NL),RHS(NP,NL),C1,C2

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/FRIC1/A,B
      COMMON/RHS1/RHS
C
C     IF(KOUNT.EQ.0)GO TO 31
C     GOTO 31
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C DON'T THINK WE WANT THIS IN INVISCID
C PROBLEM. INCLUDES FRICTION EFFECTS
C IN Q VECTOR
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C     ------------------------------------------------------------------
C
C     C1 = -2.0*SH/F
C     CZ= 1.0/(F*DELXSQ)
C     REWIND 13
C     READ (13)A
C     PRINT*, ' A = ' , A
CC    DO 5 J=2,NZM1
C     DVZ=(P(2,J+1)+P(NX,J-1)-P(NX,J+1)-P(2,J-1))/
C    1     (4.*DELX*DELZ*F)
C     RHS(1,J)=((P(2,J)+P(NX,J)-2.*P(1,J))*CZ*C1+
C    1     DVZ*A(1,J)/(F**2))*RO(J)
C     RHS(NXP1,J)=RHS(1,J)
C     DO 5 I=2,NX
C     DVZ=(P(I+1,J+1)+P(I-1,J-1)-P(I-1,J+1)-P(I+1,J-1))/
C    1        (4.*DELX*DELZ*F)
C     RHS(I,J)=((P(I+1,J)+P(I-1,J)-2.*P(I,J))*CZ*C1+
C    1      DVZ*A(I,J)/(F**2))*RO(J)
C   5 CONTINUE
C
C-----------------------------------------------------------------------
C
C     C6 = 1.0/(DELX*F)
C     C6 = C6**2
C     C3 = G/(2.0*DELX*THETA0*F**2)
C     READ (13)A
C
C     DO 7 J=2,NZM1
C     C5 = C3*RO(J)
C     C4=RO(J)*A(1,J)/F
C     RHS(1,J)=RHS(1,J)-C5*(FT1(2,J)-FT1(NX,J))+
C    1         (1.-C6*(P(2,J)+P(NX,J)-2.*P(1,J)))*C4
C     RHS(NXP1,J)=RHS(1,J)
C     DO 7 I=2,NX
C     C4 = RO(J)*A(I,J)/F
C     RHS(I,J) = RHS(I,J) -C5*(FT1(I+1,J)
C    *  - FT1(I-1,J))
C    *  + (1.0 - C6*(P(I+1,J) - 2.0*P(I,J) + P(I-1,J)))*C4
C   7 CONTINUE
C     GO TO 32
C
C     ------------------------------------------------------------------
C
   31 C1 = -2.0*SH/(DELXSQ*F**2)
      DO 1 J=2,NZM1
      C2 = C1*RO(J)
      RHS(1,J)=C2*(P(2,J)+P(NX,J)-2.*P(1,J))
      RHS(NXP1,J)=RHS(1,J)
      DO 11 I=2,NX
   11 RHS(I,J)=C2*(P(I+1,J)+P(I-1,J)-2.*P(I,J))
    1 CONTINUE
C
C     ------------------------------------------------------------------
C
   32 CONTINUE
      RETURN
      END
      SUBROUTINE UCALC
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +I,J

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),UCOEF,UCOEFR,
     +SH,ZL,EP,THETZM

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM


      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      UCOEF=G/(THETA0*F**3)
      DO 1 J=1,NZ
      UCOEFR=UCOEF*RO(J)
      DO 11 I=1,NXP1
 11   Q(I,J)=UCOEFR*Q(I,J)
 1    CONTINUE
      RETURN
      END
      SUBROUTINE QCALC
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +I,J

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),QCOEF,QCOEFR,
     +SH,ZL,EP,THETZM

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM


      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      QCOEF=(THETA0*F**3)/G
      DO 1 J=1,NZ
      QCOEFR=QCOEF/RO(J)
      DO 11 I=1,NXP1
 11   Q(I,J)=QCOEFR*Q(I,J)
 1    CONTINUE
      RETURN
      END
      SUBROUTINE VERTVEL
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +NZM2,I,J

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),
     +RQT(NP,NL),RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),PZZ(NP),WC,
     +WA,WB,SH,ZL,EP,THETZM

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM
 

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
c     DIMENSION PZZ(NP)
C-----------------------------------------------------------------------
C
C  CALCULATE THE VERTICAL VELOCITY FIELD
C
C-----------------------------------------------------------------------
C
      WC=-G*DELZ**2/(2.0*DELX*THETA0*F)
C
C  PSI HAS PREIODIC BOUNDARY CONDITIONS
C
      DO 1 J=2,NZM1
      WA=WC*(PSI(2,J)-PSI(NX,J))*Q(1,J)
      WB=P(1,J+1)-2.0*P(1,J)+P(1,J-1)
      W(1,J)=WA/WB
      DO 11 I=2,NX
      WA=WC*(PSI(I+1,J)-PSI(I-1,J))*Q(I,J)
      WB=P(I,J+1)-2.0*P(I,J)+P(I,J-1)
      W(I,J)=WA/WB
 11   CONTINUE
      W(NXP1,J)=W(1,J)
 1    CONTINUE
C
C
C-----------------------------------------------------------------------
C
C  SET W TO ZERO ON TOP AND BOTTOM BOUNDARIES
      NZM2=NZM1-1
C
C-----------------------------------------------------------------------
C
      DO 3 I=1,NXP1
      W(I,1)=0.0
      W(I,NZ)=0.0
      PZZ(I)=(P(I,NZM2)-2.*P(I,NZM1)+P(I,NZ))
 3    CONTINUE
C
      IF(W0.EQ.0) GOTO 118
      WA=WC*(PSI(2,NZ)-PSI(NX,NZ))*Q(1,NZ)
      WB=PZZ(1)
      W(1,NZ)=WA/WB
      DO 112 I=2,NX
      WA=WC*(PSI(I+1,NZ)-PSI(I-1,NZ))*Q(I,NZ)
      WB=PZZ(I)
      W(I,NZ)=WA/WB
 112  CONTINUE
      W(NXP1,NZ)=W(1,NZ)
C     PRINT*, ' VERT VEL FROM SUBR VERT'
C      CALL OUTP(W,1,4)
 118  CONTINUE
      RETURN
      END
      SUBROUTINE ADVPHIZ(KOUNT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +   KOUNT,I,IK

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +FC,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),
     +RQT(NP,NL),RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),F1(NP,NL),
     +FT1(NP,NL),PZB1(NP),PZB2(NP),PZT1(NP),PZT2(NP),
     +QINT1(NP,NL),VGINT(NP,NL),PSIINT(NP,NL),XSEMI(NP,NL),
     +WI(NP,NL),PZBI(NP),PZTI(NP),PV1,PV2,DPZT,DELT,T1,TF,CONST,
     +CONSTT

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/FC,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/STORE/PZB1,PZB2,PZT1,PZT2
      COMMON/INT/QINT1,VGINT,PSIINT,XSEMI
     +,WI,PZBI,PZTI
      COMMON/PVADJ/PV1,PV2,DPZT

      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL


      IK=IK+1
C
      DELT=2.*TIMINC
      IF(KOUNT.EQ.1)GOTO 30
C-----------------------------------------------------------------------
C
C  ADVECT PHION TOP AND BOTTOM BOUNDARIES
C  USING SECOND ORDER CENTERED DIFFERENCES
C  FOR STABILITY WE REQUIRE THAT TIMINC.LT.DELX/(SH*ZL*0.5)
C
C     DELT=2.0*TIMINC
 5    T1 = DELT/(2.0*DELX)
      TF = G*DELT/THETA0
      CONST=SH*ZL*0.5
        CONSTT=CONST
       CONST=0.
       CONSTT=SH*ZL


C	CONST=-CONSTT/2.

C	CONSTT=-CONST

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ZERO WIND AT 3KM 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      CONST=-SH*3000.
      CONSTT=CONST+SH*ZL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C 0M/S WIND AT 0KM 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       CONST=-0.
       CONSTT=CONST+SH*ZL
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC


      DO 10 I=2,NX
      PZB2(I)=PZB1(I)-(CONST)*(PZB(I+1)-PZB(I-1))*T1
     1        +SH*(P(I+1,1)-P(I-1,1))*T1
      PZT2(I)=PZT1(I)-CONSTT*(PZT(I+1)-PZT(I-1))*T1
     1        +SH*(P(I+1,NZ)-P(I-1,NZ))*T1   


   10 CONTINUE
      PZB2(1)=PZB1(1)-(CONST)*(PZB(2)-PZB(NX))*T1
     1        +SH*(P(2,1)-P(NX,1))*T1 


      PZT2(1)=PZT1(1)-CONSTT*(PZT(2)-PZT(NX))*T1
     1        +SH*(P(2,NZ)-P(NX,NZ))*T1
      PZB2(NXP1)=PZB2(1)

      PZT2(NXP1)=PZT2(1)
      DO 20 I=1,NXP1
      PZB1(I) = PZB(I)
        PZB(I)=PZB2(I)
      PZT1(I) = PZT(I)



        IF(W0.EQ.0) PZT(I)=PZT2(I)
        IF(W0.EQ.1) PZT(I)=PZT1(I)

 20   CONTINUE
C
C-----------------------------------------------------------------------
C
      RETURN
C
C-----------------------------------------------------------------------
C
C  ADJUST SUCH THAT THE FIRST TIME STEP IS A FORWARD DIFFERENCE
C
 30   DELT=TIMINC
      DO 40 I=1,NXP1
      PZB1(I)=PZB(I)
      PZT1(I)=PZT(I)
 40   CONTINUE
C
      GOTO 5
C
C-----------------------------------------------------------------------
 
 
C
      END
      SUBROUTINE ADVQFXDI(KOUNT,DAY)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +NZM2,I,J,KOUNT,N,IQ

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),DAY,
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),DX(NP,NL),FRIC(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),QINT1(NP,NL),VGINT(NP,NL),
     +PSIINT(NP,NL),XSEMI(NP,NL),WI(NP,NL),PZBI(NP),PZTI(NP),
     +DELSQ(NP,NL),ZMT(NP,NL),VMT(NP,NL),THMTT(NP),THMTB(NP),
     +CT,CTZ,CV,DELT,RD,RDZ,UB,
     +AXZ,FIPH,FIMH,FJPH,FJMH,RQP,U,PVTEN,QTR(NP,NL)


      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM


      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/GAUSS2/AL,T
      COMMON/FRIC1/DX,FRIC
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/INT/QINT1,VGINT,PSIINT,XSEMI
     +,WI,PZBI,PZTI
	COMMON/QFORCE/QTR
      COMMON/UNIF/IQ


      COMMON/DIFF/DELSQ,ZMT,VMT,THMTT,THMTB
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C   THIS SUBROUTINE ADVECTS PV USING THE C
C   FLUX FORM OF THE ADVECTION EQUATION  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C                                        C
C   MULTIPLY PV BY DENSITY               C
C                                        C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      NZM2=NZM1-1
      CT = THETA0/G
      CTZ= CT/(2.0*DELZ)
      CV = 1.0/(F*2.0*DELX)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      PRINT*, ' PV BEFORE ADVECTION '
C      CALL OUTP(Q)
C      IF(RD.GT.1.OR.RDZ.GT.1) PRINT*, ' STOP: RD RDZ ' , RD,RDZ
C      IF(RD.GT.1.OR.RDZ.GT.1)  STOP  124
      DO 10 J=1,NL
        DO 11 I=1,NP
C         RQT(I,J)=QINT1(I,J)*RO(J)
          RQT(I,J)=Q(I,J)*RO(J)
          WI(I,J)=W(I,J)
 11     CONTINUE
 10   CONTINUE
      DELT=TIMINC
      NZM2=NZM1-1
      IF(KOUNT.EQ.1) THEN
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C IF KOUNT EQUALS 1 DO A FORWARD TIME STEP C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 12 J=1,NL
        DO 13 I=1,NP
          RQMT(I,J)=RQT(I,J)
C      RQMT(I,J)=RO(J)
 13     CONTINUE
 12   CONTINUE
      DELT=DELT/2.
      ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCC  SET TO UNIFORM Q CCCCC
C      IF(IQ.EQ.9.AND.DAY.GE.1) THEN
C         IQ=1
C         DO 145 J=1,NL
C           DO 146 I=1,NP
C             RQMT(I,J)=(0.25E-6)*RO(J)
C             RQT(I,J)=RQMT(I,J)
C 146       CONTINUE
C 145  CONTINUE
C      ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      CALL DIFFUSE(KOUNT)
      
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ADVECTION ON BOTTOM BOUNDARY       C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C DIFFUSION COEFFICEIENTS             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      RD=EP*DELT/(DELX*DELX)
      RDZ=EP*DELT/(DELZ*DELZ)
      UB=-SH*ZL*0.5
       UB=0.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC ZERO WIND AT 3KM
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      UB=-SH*3000.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC 0M/S WIND AT 0KM
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
       UB=0.
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      AXZ=1./(DELX*DELZ)
      FIPH=0.5*((UB)*RQT(2,1)+(UB)*RQT(1,1))
      FIMH=0.5*((UB)*RQT(1,1)+(UB)*RQT(NX,1))
      FJPH=WI(1,2)*RQT(1,2)
      FJMH=0.
      RQP=RQMT(1,1)-2.*DELT*(((FIPH-FIMH)*DELZ)+
     +   ((FJPH-FJMH)*DELX))*AXZ
      RQP=RQP+2.*DELT*Q(1,1)*RO(1)*WI(1,2)/DELZ
      RQP=RQP+EP*DELSQ(1,1)*2.*DELT
      RQS(1,1)=RQP
      RQS(NXP1,1)=RQS(1,1)
      DO 20 I=2,NX
        FIPH=0.5*((UB)*RQT(I+1,1)+(UB)*RQT(I,1))
        FIMH=0.5*((UB)*RQT(I,1)+(UB)*RQT(I-1,1))
        FJPH=WI(I,2)*RQT(I,2)
        FJMH=0.
        RQP=RQMT(I,1)-2.*DELT*(((FIPH-FIMH)*DELZ)+
     +      ((FJPH-FJMH)*DELX))*AXZ
        RQP=RQP+2.*DELT*Q(I,1)*RO(1)*WI(I,2)/DELZ
        RQP=RQP+EP*DELSQ(I,1)*2.*DELT
        RQS(I,1)=RQP
 20   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  ADVECTION OF Q AT INTERMEDIATE LEVELS      C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C CALCULATE RQS AT LEVEL 2
C
      U=UB+DELZ*SH
      FIPH=0.5*((U)*RQT(2,2)+(U)*RQT(1,2))
      FIMH=0.5*((U)*RQT(1,2)+(U)*RQT(NX,2))
      FJPH=0.5*(WI(1,3)*RQT(1,3)+WI(1,2)*RQT(1,2))
      FJMH=0.5*WI(1,2)*RQT(1,2)
      RQP=RQMT(1,2)-2.*DELT*(((FIPH-FIMH)*DELZ)
     ++((FJPH-FJMH)*DELX))*AXZ
      RQP=RQP+2.*DELT*Q(1,2)*(RO(3)*WI(1,3)-RO(1)*WI(1,1))/(2.*DELZ)
      RQP=RQP+EP*DELSQ(1,2)*2.*DELT
      RQS(1,2)=RQP
      RQS(NXP1,2)=RQS(1,2)
      DO 40 I=2,NX
        FIPH=0.5*((U)*RQT(I+1,2)+(U)*RQT(I,2))
        FIMH=0.5*((U)*RQT(I,2)+(U)*RQT(I-1,2))
        FJPH=0.5*(WI(I,3)*RQT(I,3)+WI(I,2)*RQT(I,2))
        FJMH=0.5*WI(I,2)*RQT(I,2)
        RQP=RQMT(I,2)-2.*DELT*(((FIPH-FIMH)*DELZ)
     +  +((FJPH-FJMH)*DELX))*AXZ
        RQP=RQP+2.*DELT*Q(I,2)*(RO(3)*WI(I,3)-RO(1)*WI(I,1))/(2.*DELZ)
        RQP=RQP+EP*DELSQ(I,2)*2.*DELT
        RQS(I,2)=RQP
 40   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATE RQS FOR LEVELS 3 TO NZM2             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 50 J=3,NZM2
        U=UB+DELZ*(J-1)*SH
        FIPH=0.5*((U)*RQT(2,J)+(U)*RQT(1,J))
        FIMH=0.5*((U)*RQT(1,J)+(U)*RQT(NX,J))
        FJPH=0.5*(WI(1,J+1)*RQT(1,J+1)+WI(1,J)*RQT(1,J))
        FJMH=0.5*(WI(1,J)*RQT(1,J)+WI(1,J-1)*RQT(1,J-1))
        RQP=RQMT(1,J)-2.*DELT*(((FIPH-FIMH)*DELZ)
     +      +((FJPH-FJMH)*DELX))*AXZ
        RQP=RQP+Q(1,J)*(RO(J+1)*WI(1,J+1)-RO(J-1)*WI(1,J-1))*2.*DELT
     +      /(2.*DELZ)
        RQP=RQP+EP*DELSQ(1,J)*2.*DELT
        RQS(1,J)=RQP
      RQS(NXP1,J)=RQS(1,J)
        DO 51 I=2,NX
          FIPH=0.5*((U)*RQT(I+1,J)+(U)*RQT(I,J))
          FIMH=0.5*((U)*RQT(I,J)+(U)*RQT(I-1,J))
          FJPH=0.5*(WI(I,J+1)*RQT(I,J+1)+WI(I,J)*RQT(I,J))
          FJMH=0.5*(WI(I,J)*RQT(I,J)+WI(I,J-1)*RQT(I,J-1))
          RQP=RQMT(I,J)-2.*DELT*(((FIPH-FIMH)*DELZ)
     +        +((FJPH-FJMH)*DELX))*AXZ
          RQP=RQP+Q(I,J)*(RO(J+1)*WI(I,J+1)-RO(J-1)*WI(I,J-1))*2.*DELT
     +       /(2.*DELZ)
          RQP=RQP+EP*DELSQ(I,J)*2.*DELT
          RQS(I,J)=RQP
      IF(J.EQ.20) THEN
         PVTEN=EP*DELSQ(I,J)
          PVTEN=PVTEN*60.*60.*24.*1.E6
C        PRINT*, ' PV TENDENCY AT I = ', I, 'IS ' ,PVTEN , 'PV PER DAY '
         ENDIF
 51     CONTINUE
 50   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATE RQS AT LEVEL NZM1               C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C LEFT BOUNDARY
C
      U=UB+(NZM1-1)*DELZ*SH
      N=NZM1
      FIPH=0.5*((U)*RQT(2,N)+(U)*RQT(1,N))
      FIMH=0.5*((U)*RQT(1,N)+(U)*RQT(NX,N))
      FJPH=0.5*(WI(1,NZM1)*RQT(1,NZM1)+WI(1,NZ)*RQT(1,NZ))
      FJMH=0.5*(WI(1,NZM1)*RQT(1,NZM1)+WI(1,NZM2)*RQT(1,NZM2))
      RQP=RQMT(1,NZM1)-2.*DELT*(((FIPH-FIMH)*DELZ)
     +    +((FJPH-FJMH)*DELX))*AXZ
      RQP=RQP+Q(1,N)*(RO(N+1)*WI(1,N+1)-RO(N-1)*WI(1,N-1))*2.*DELT
     +    /(2.*DELZ)
      RQP=RQP+EP*DELSQ(1,NZM1)*2.*DELT
      RQS(1,NZM1)=RQP
      RQS(NXP1,NZM1)=RQS(1,NZM1)
      DO 60 I=2,NX
        FIPH=0.5*((U)*RQT(I+1,N)+(U)*RQT(I,N))
        FIMH=0.5*((U)*RQT(I,N)+(U)*RQT(I-1,N))
        FJPH=0.5*(WI(I,NZM1)*RQT(I,NZM1)+WI(I,NZ)*RQT(I,NZ))
        FJMH=0.5*(WI(I,NZM1)*RQT(I,NZM1)+WI(I,NZM2)*RQT(I,NZM2))
        RQP=RQMT(I,NZM1)-2.*DELT*(((FIPH-FIMH)*DELZ)
     +      +((FJPH-FJMH)*DELX))*AXZ
        RQP=RQP+Q(I,N)*(RO(N+1)*WI(I,N+1)-RO(N-1)*WI(I,N-1))*2.*DELT
     +    /(2.*DELZ)
        RQP=RQP+EP*DELSQ(I,NZM1)*2.*DELT
        RQS(I,NZM1)=RQP
 60   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATE RQS ON TOP BOUNDARY                 C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C
C LEFT BOUNDARY
C
      U=UB+(NZ-1)*SH*DELZ
      FIPH=0.5*((U)*RQT(2,NZ)+(U)*RQT(1,NZ))
      FIMH=0.5*((U)*RQT(1,NZ)+(U)*RQT(NX,NZ))
      FJPH=0.
      FJPH=WI(1,NZ)*RQT(1,NZ)
      FJMH=WI(1,NZM1)*RQT(1,NZM1)
      RQP=RQMT(1,NZ)-2.*DELT*(((FIPH-FIMH)*DELZ)+
     +    ((FJPH-FJMH)*DELX))*AXZ
      RQP=RQP-2.*DELT*Q(1,NZ)*WI(1,NZM1)*RO(NZM1)/DELZ
      RQP=RQP+2.*DELT*Q(1,NZ)*WI(1,NZ)*RO(NZ)/DELZ
      RQP=RQP+EP*DELSQ(1,NZ)*2.*DELT
      RQS(1,NZ)=RQP
      RQS(NXP1,NZ)=RQS(1,NZ)
      DO 70 I=2,NX
        FIPH=0.5*((U)*RQT(I+1,NZ)+(U)*RQT(I,NZ))
        FIMH=0.5*((U)*RQT(I,NZ)+(U)*RQT(I-1,NZ))
        FJPH=0.
        FJPH=WI(I,NZ)*RQT(I,NZ)
        FJMH=WI(I,NZM1)*RQT(I,NZM1)
        RQP=RQMT(I,NZ)-2.*DELT*(((FIPH-FIMH)*DELZ)+
     +      ((FJPH-FJMH)*DELX))*AXZ
        RQP=RQP-2.*DELT*Q(I,NZ)*WI(I,NZM1)*RO(NZM1)/DELZ
        RQP=RQP+2.*DELT*Q(I,NZ)*WI(I,NZ)*RO(NZ)/DELZ
        RQP=RQP+EP*DELSQ(I,NZ)*2.*DELT
        RQS(I,NZ)=RQP
 70   CONTINUE

CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ADD IN Q FORCING LIKE A TROUGH
C  0.0069 IS EQUIVALENT TO 1PVU/DAY 
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
	IF(DAY.GE.100) THEN
	DO 734 J=1,NL
		DO 735 I=1,NP
		RQS(I,J)=RQS(I,J)+QTR(I,J)*RO(J)*0.0069
 735		CONTINUE
 734	CONTINUE
	ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

      DO 800 J=1,NL
        DO 801 I=1,NP
          RQMT(I,J)=RQT(I,J)
          QINT1(I,J)=RQS(I,J)/RO(J)
          Q(I,J)=RQS(I,J)/RO(J)
 801    CONTINUE
 800  CONTINUE
C           CALL BETATERM(KOUNT)
C      PRINT*, ' POT VORT AFTER TIME STEP'
C      CALL OUTP(Q)
      RETURN
      END
      SUBROUTINE DIFFUSE(KOUNT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +I,J,NZM2,KOUNT


      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),DX(NP,NL),FRIC(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),QINT1(NP,NL),VGINT(NP,NL),
     +PSIINT(NP,NL),XSEMI(NP,NL),WI(NP,NL),PZBI(NP),PZTI(NP),
     +DELSQ(NP,NL),ZMT(NP,NL),VMT(NP,NL),THMTT(NP),THMTB(NP),UAG(NP,NL),
     +VGG(NP,NL),VORT(NP,NL)
     +,RQX(NP,NL),RQZ(NP,NL),ZX(NP,NL),ZZ(NP,NL)
     +,VZ(NP,NL),RQZZ(NP,NL),VZZ(NP,NL),RQXX(NP,NL),VZX(NP,NL)
     +,RQZX(NP,NL)
     +,THT(NP),THB(NP),THXT(NP),THXB(NP),CT,F0,
     +T1,T2,T3,T4,T5,T6,T7,T8,T9,T10

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/GAUSS2/AL,T
      COMMON/FRIC1/DX,FRIC
      COMMON/UAGEO/UAG
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/INT/QINT1,VGINT,PSIINT,XSEMI
     +,WI,PZBI,PZTI
      COMMON/DIFF/DELSQ,ZMT,VMT,THMTT,THMTB
c      DIMENSION VGG(NP,NL),VORT(NP,NL)
c     +,RQX(NP,NL),RQZ(NP,NL),ZX(NP,NL),ZZ(NP,NL)
c     +,VZ(NP,NL),RQZZ(NP,NL),VZZ(NP,NL),RQXX(NP,NL),VZX(NP,NL)
c     +,RQZX(NP,NL)
c     +,THT(NP),THB(NP),THXT(NP),THXB(NP)
      CT = THETA0/G
      DO 111 I=1,NXP1
      THB(I) = CT*PZB(I)
      THT(I)=CT*PZT(I)
  111 CONTINUE
      DO 2 J=1,NL
      VGG(1,J)=(P(2,J)-P(NX,J))/(2.0E-04*DELX)
      VGG(NXP1,J)=VGG(1,J)
      DO 2 I=2,(NP-1)
      VGG(I,J)=(P(I+1,J)-P(I-1,J))/(2.0E-4*DELX)
 2    CONTINUE
      F0=1.E-4
      NZM2=NZM1-1
        VORT(1,1)=1./(1.-(VGG(2,1)-VGG(NX,1))/(F0*2.*DELX))
        VORT(1,NL)=1./(1.-(VGG(2,NL)-VGG(NX,NL))/(F0*2.*DELX))
        VORT(NXP1,1)=VORT(1,1)
        VORT(NXP1,NL)=VORT(1,NL)
        DO 3 I=2,NX
          VORT(I,1)=1./(1.-(VGG(I+1,1)-VGG(I-1,1))/(F0*2.*DELX))
          VORT(I,NL)=1./(1.-(VGG(I+1,NL)-VGG(I-1,NL))/(F0*2.*DELX))
 3    CONTINUE
      DO 4 J=2,NL-1
        VORT(1,J)=1./(1.-(VGG(2,J)-VGG(NX,J))/(F0*2.*DELX))
        VORT(NXP1,J)=VORT(1,J)
        DO 4 I=2,NX
          VORT(I,J)=1./(1.-(VGG(I+1,J)-VGG(I-1,J))/(F0*2.*DELX))
 4    CONTINUE
      IF(KOUNT.EQ.1) THEN
         DO 5 J=1,NL
           DO 6 I=1,NXP1
               ZMT(I,J)=VORT(I,J)
               VMT(I,J)=VGG(I,J)
               THMTB(I)=THB(I)
               THMTT(I)=THT(I)
 6         CONTINUE
 5      CONTINUE
      ENDIF
C
C
C FIRST VERTICAL DERIVATIVES
C      PRINT*, ' FIRST VERT'
C
      DO 10 J=2,NZM1
        DO 11 I=1,NXP1
           RQZ(I,J)=(RQMT(I,J+1)-RQMT(I,J-1))/(2.*DELZ)
           ZZ(I,J)=(ZMT(I,J+1)-ZMT(I,J-1))/(2.*DELZ)
           VZ(I,J)=(VMT(I,J+1)-VMT(I,J-1))/(2.*DELZ)
 11     CONTINUE
 10   CONTINUE
C
C  FIRST X DERIVATIVES
C      PRINT*, ' FIRST X DER'
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
        THXB(1)=(THMTB(2)-THMTB(NX))/(2.*DELX)
        VZ(1,1)=THXB(1)/(CT*F0)
        THXB(NXP1)=THXB(1)
        VZ(NXP1,1)=VZ(1,1)
        THXT(1)=(THMTT(2)-THMTT(NX))/(2.*DELX)
        VZ(1,NL)=THXT(1)/(CT*F0)
        THXT(NXP1)=THXT(1)
        VZ(NXP1,NL)=VZ(1,NL)
      DO     700 I=2,NX
           THXB(I)=(THMTB(I+1)-THMTB(I-1))/(2.*DELX)
           VZ(I,1)=THXB(I)/(CT*F0)
           THXT(I)=(THMTT(I+1)-THMTT(I-1))/(2.*DELX)
           VZ(I,NL)=THXT(I)/(CT*F0)
 700  CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      PRINT*, ' HI THERE'
      DO 7 J=1,NL
        RQX(1,J)=(RQMT(2,J)-RQMT(NX,J))/(2.*DELX)
        RQX(NXP1,J)=RQX(1,J)
        ZX(1,J)=(ZMT(2,J)-ZMT(NX,J))/(2.*DELX)
        ZX(NXP1,J)=ZX(1,J)
        DO 8 I=2,NX
           RQX(I,J)=(RQMT(I+1,J)-RQMT(I-1,J))/(2.*DELX)
           ZX(I,J)=(ZMT(I+1,J)-ZMT(I-1,J))/(2.*DELX)
 8      CONTINUE
 7    CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      PRINT*, ' HI THERE 2'
C   CALC RQZX INCORPORATING BOUNDARY CONDITIONS  C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 300 I=1,NXP1
          RQZ(I,1)=-ZMT(I,1)*VZ(I,1)*RQX(I,1)/F0
          RQZ(I,NL)=-ZMT(I,NL)*VZ(I,NL)*RQX(I,NL)/F0
 300  CONTINUE
      DO 301 J=2,NZM1
        RQZX(1,J)=(RQZ(2,J)-RQZ(NX,J))/(2.*DELX)
        VZX(1,J)=(VZ(2,J)-VZ(NX,J))/(2.*DELX)
        VZX(NXP1,J)=VZX(1,J)
        RQZX(NXP1,J)=RQZX(1,J)
        DO 302 I=2,NX
           VZX(I,J)=(VZ(I+1,J)-VZ(I-1,J))/(2.*DELX)
          RQZX(I,J)=(RQZ(I+1,J)-RQZ(I-1,J))/(2.*DELX)
 302    CONTINUE
 301  CONTINUE
CC
CC SECOND VERTICAL DERIVATIVES
CC
C      PRINT*, ' SECOND VERTICAL '
      DO 13 J=2,NZM1
        DO 14 I=1,NXP1
          RQZZ(I,J)=(RQMT(I,J+1)-2.*RQMT(I,J)+RQMT(I,J-1))/(DELZ**2)
           VZZ(I,J)=(VMT(I,J+1)-2.*VMT(I,J)+VMT(I,J-1))/(DELZ**2)
 14     CONTINUE
 13   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    RQZZ AND VZZ AT TOP AND BOTTOM                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 400 I=1,NXP1
         VZZ(I,1)=2.*VZZ(I,2)-VZZ(I,3)
         VZZ(I,NL)=2.*VZZ(I,NZM1)-VZZ(I,NZM2)
         ZZ(I,1)=2.*ZZ(I,2)-ZZ(I,3)
         ZZ(I,NL)=2.*ZZ(I,NZM1)-ZZ(I,NZM2)
         RQZZ(I,1)=(2.*(RQMT(I,2)-RQMT(I,1))+
     +      2.*DELZ*(ZMT(I,1)*THXB(I)*RQX(I,1)/(CT*F0*F0)))/(DELZ**2)
         RQZZ(I,NL)=(2.*(RQMT(I,NZM1)-RQMT(I,NL))+
     +      2.*DELZ*(ZMT(I,NL)*THXT(I)*RQX(I,NL)/(CT*F0*F0)))/(DELZ**2)
 400  CONTINUE
C
C SECOND HORIZONTAL DERIVATIVES
C      PRINT*, ' SECOND HORIZ '
      DO 16 J=1,NL
        RQXX(1,J)=(RQMT(2,J)-2.*RQMT(1,J)+RQMT(NX,J))/(DELX**2)
        RQXX(NXP1,J)=RQXX(1,J)
        DO 17 I=2,NX
          RQXX(I,J)=(RQMT(I+1,J)-2.*RQMT(I,J)+RQMT(I-1,J))/(DELX**2)
 17   CONTINUE
 16   CONTINUE
      VZX(1,1)=(THMTB(2)-2.*THMTB(1)+THMTB(NX))/(DELX*DELX*CT*F0)
      VZX(NXP1,1)=VZX(1,1)
      VZX(1,NL)=(THMTT(2)-2.*THMTT(1)+THMTT(NX))/(DELX*DELX*CT*F0)
      VZX(NXP1,NL)=VZX(1,NL)
      DO 555 I=2,NX
         VZX(I,1)=(THMTB(I+1)-2.*THMTB(I)+THMTB(I-1))
     +    /(DELX*DELX*CT*F0)
         VZX(I,NL)=(THMTT(I+1)-2.*THMTT(I)+THMTT(I-1))
     +  /(DELX*DELX*CT*F0)
 555  CONTINUE
C      PRINT*, ' HI THERE 4'
      DO 303 I=1,NXP1
      RQZX(I,1)=(-ZX(I,1)*VZ(I,1)*RQX(I,1)-ZMT(I,1)*VZX(I,1)*RQX(I,1)
     +           -ZMT(I,1)*VZ(I,1)*RQXX(I,1))/F0
      RQZX(I,NL)=(-ZX(I,NL)*VZ(I,1)*RQX(I,NL)
     +-ZMT(I,NL)*VZX(I,NL)*RQX(I,NL)-ZMT(I,NL)*VZ(I,NL)*RQXX(I,NL))/F0
 303  CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C INTERIOR DIFFUSION TERMS   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 100 J=2,NZM1
         DO 101 I=1,NXP1
            T1=ZMT(I,J)*ZX(I,J)*RQX(I,J)
            T2=ZMT(I,J)*ZMT(I,J)*RQXX(I,J)
            T3=ZMT(I,J)*VZ(I,J)/F0
            T4=ZX(I,J)*VZ(I,J)*RQX(I,J)/F0
            T5=ZMT(I,J)*VZX(I,J)*RQX(I,J)/F0
            T6=ZMT(I,J)*VZ(I,J)*RQXX(I,J)/F0
            T7=2.*RQZX(I,J)
            T8=ZZ(I,J)*VZ(I,J)*RQX(I,J)/F0
            T9=ZMT(I,J)*VZZ(I,J)*RQX(I,J)/F0
            T10=RQZZ(I,J)
            DELSQ(I,J)=T1+T2+T3*(T4+T5+T6+T7)+T8+T9+T10
 101     CONTINUE
 100   CONTINUE
         DO 104 I=1,NXP1
            T1=ZMT(I,1)*ZX(I,1)*RQX(I,1)
            T2=ZMT(I,1)*ZMT(I,1)*RQXX(I,1)
            T3=ZMT(I,1)*VZ(I,1)/F0
            T4=ZX(I,1)*VZ(I,1)*RQX(I,1)/F0
            T5=ZMT(I,1)*VZX(I,1)*RQX(I,1)/F0
            T6=ZMT(I,1)*VZ(I,1)*RQXX(I,1)/F0
            T7=2.*RQZX(I,1)
            T8=ZZ(I,1)*VZ(I,1)*RQX(I,1)/F0
            T9=ZMT(I,1)*VZZ(I,1)*RQX(I,1)/F0
            T10=RQZZ(I,1)
            DELSQ(I,1)=T1+T2-T3*(T4+T5+T6)+T8+T9+T10
            T1=ZMT(I,NL)*ZX(I,NL)*RQX(I,NL)
            T2=ZMT(I,NL)*ZMT(I,NL)*RQXX(I,NL)
            T3=ZMT(I,NL)*VZ(I,NL)/F0
            T4=ZX(I,NL)*VZ(I,NL)*RQX(I,NL)/F0
            T5=ZMT(I,NL)*VZX(I,NL)*RQX(I,NL)/F0
            T6=ZMT(I,NL)*VZ(I,NL)*RQXX(I,NL)/F0
            T7=2.*RQZX(I,NL)
            T8=ZZ(I,NL)*VZ(I,NL)*RQX(I,NL)/F0
            T9=ZMT(I,NL)*VZZ(I,NL)*RQX(I,NL)/F0
            T10=RQZZ(I,NL)
            DELSQ(I,NL)=T1+T2-T3*(T4+T5+T6)+T8+T9+T10
 104     CONTINUE
C        DO 103 I=1,NXP1
C           T1=ZMT(I,1)*ZX(I,1)*RQX(I,1)
C           T2=ZMT(I,1)*ZMT(I,1)*RQXX(I,1)
C           DELSQ(I,1)=T1+T2
C           T1=ZMT(I,NL)*ZX(I,NL)*RQX(I,NL)
C           T2=ZMT(I,NL)*ZMT(I,NL)*RQXX(I,NL)
C           DELSQ(I,NL)=T1+T2
C103     CONTINUE
CCCCCCCCCC SAVE CURRENT FIELDS FOR DIFFUSION IN NEXT TIME STEP
         DO 200 J=1,NL
           DO 210 I=1,NXP1
               ZMT(I,J)=VORT(I,J)
               VMT(I,J)=VGG(I,J)
               THMTB(I)=THT(I)
               THMTT(I)=THB(I)
 210        CONTINUE
 200    CONTINUE
C      IF(KOUNT.EQ.50.OR.KOUNT.EQ.100) THEN
C          PRINT*, ' DIFFUSION FIELD '
C         CALL OUTP(DELSQ,100,7)
C Cc     ENDIF
      RETURN
      END
      SUBROUTINE BETATERM(KOUNT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +NZM2,J,KOUNT,I,L


      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),DX(NP,NL),FRIC(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),QINT1(NP,NL),VGINT(NP,NL),
     +PSIINT(NP,NL),XSEMI(NP,NL),WI(NP,NL),PZBI(NP),PZTI(NP),
     +DELSQ(NP,NL),ZMT(NP,NL),VMT(NP,NL),THMTT(NP),THMTB(NP),
     +VG(NP,NL),THZ(NP,NL),BETA,QMT(NP,NL),BETA,CT,CTZ,CV,DELT

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/GAUSS2/AL,T
      COMMON/FRIC1/DX,FRIC
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/INT/QINT1,VGINT,PSIINT,XSEMI
     +,WI,PZBI,PZTI
      COMMON/DIFF/DELSQ,ZMT,VMT,THMTT,THMTB
C     DIMENSION VG(NP,NL),THZ(NP,NL),BETA,QMT(NP,NL)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C   THIS SUBROUTINE ADDS THE BETA TERM   C
C   TO THE PV EQUATION                   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      NZM2=NZM1-1
      BETA=1.02E-11
      CT = THETA0/G
      CTZ= CT/(2.0*DELZ)
      CV = 1.0/(F*2.0*DELX)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DELT=TIMINC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C IF KOUNT EQUALS 1 DO A FORWARD TIME STEP C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 12 J=1,NL
        DO 13 I=1,NP
           QMT(I,J)=Q(I,J)
 13     CONTINUE
 12   CONTINUE
      IF(KOUNT.EQ.1) THEN
      DELT=DELT/2.
      ENDIF
      DO 2 J=1,NL
      VG(1,J)=(P(2,J)-P(NX,J))/(2.0E-04*DELX)
      VG(NXP1,J)=VG(1,J)
      DO 2 I=2,(NP-1)
      VG(I,J)=(P(I+1,J)-P(I-1,J))/(2.0E-4*DELX)
 2    CONTINUE
CC SECOND VERTICAL DERIVATIVES
CC
C      PRINT*, ' SECOND VERTICAL '
      DO 23 J=2,NZM1
        DO 14 I=1,NXP1
          THZ(I,J)=(P(I,J+1)-2.*P(I,J)+P(I,J-1))/(DELZ**2)
     +              *THETA0/G
 14     CONTINUE
 23   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C    THZ          AT TOP AND BOTTOM                C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 3 I=1,NXP1
      THZ(I,NZ)=(P(I,NZM2)-2.*P(I,NZM1)+P(I,NZ))/(DELZ**2)
     +              *THETA0/G
      THZ(I,1)=(P(I,1)-2.*P(I,2)+P(I,3))/(DELZ**2)
     +              *THETA0/G
 3    CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C ADVECTION OF PLANETARY VORTICITY   C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 40 L=1,NL
      RQS(1,L)= QMT(1,L)-2.*DELT*BETA*THZ(1,L)*VG(1,L)/RO(L)
      RQS(NXP1,L)=RQS(1,L)
      DO 20 I=2,NX
      RQS(I,L)= QMT(I,L)-2.*DELT*BETA*THZ(I,L)*VG(I,L)/RO(L)
 20   CONTINUE
 40   CONTINUE
      DO 800 J=1,NL
        DO 801 I=1,NP
          Q(I,J)=RQS(I,J)
 801    CONTINUE
 800  CONTINUE
C      PRINT*, ' POT VORT AFTER TIME STEP'
C      CALL OUTP(Q)
      RETURN
      END
      SUBROUTINE XCALC(KOUNT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +J,I,KOUNT

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),DX(NP,NL),V(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),VG(NP,NL),X

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ               
       COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/GAUSS2/AL,T
      COMMON/FRIC1/DX,V
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/VGEOS/VG
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C THIS SUBROUTINE CALCULATES X IN PHYSICAL SPACE C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 1 J=1,NL
      XREAL(1,J)=XA-VG(1,J)/F
      DO 2 I=2,NX
        X=DELX*(I-1)+XA
        XREAL(I,J)=X-VG(I,J)/F
 2    CONTINUE
C     XREAL(NXP1,J)=XB-VG(NXP1,J)/F
 1    CONTINUE
C      CALL OUTP(XREAL)
      RETURN
      END
      SUBROUTINE INTVORT(NT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +J,NZM2,I,II,IR,IN,NT

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),DX(NP,NL),V(NP,NL),
     +Q1(NP,NL),Q2(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),VG(NP,NL),
     +QINT1(NP,NL),VGINT(NP,NL),PSIINT(NP,NL),XSEMI(NP,NL)
     +,WI(NP,NL),PZBI(NP),PZTI(NP),
     +VORTR(NP,NL),VORT(NP,NL)
     +,RNSR(NP,NL),RNSAV(NL),VORTTS(NP,NL),RNSTS(NP,NL),F0,XREG,
     +WT1,WT2,DELXR,VORTGS(NP,NL),VGS(NP,NL),TS(NP,NL),TR(NP,NL),
     +CT

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM
      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/GAUSS2/AL,T
      COMMON/FRIC1/DX,V
      COMMON/STOREQW/Q1,Q2,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/INT/QINT1,VGINT,PSIINT,XSEMI
     +,WI,PZBI,PZTI
      COMMON/VGEOS/VG
c      DIMENSION  VORTR(NP,NL),VORT
c     +,RNSR(NP,NL),RNSAV(NL),TR(NP
      COMMON/VORT2/VORTTS,RNSTS,VORTGS,
     :VGS,TS


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C THIS SUBROUTINE INTERPOLATES Q,VG AND PSI TO REGULAR GRID   C
C IN PHYSICAL SPACE                                           C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      DO 133 J=1,NL
      RNSAV(J)=0.
 133  CONTINUE
      NZM2=NZ-2
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATE N-SQUARED          C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      DO 112 I=1,NX
C        DO 113 J=2,NZM1
C          RNS(I,J)=(P(I,J+1)-2.*P(I,J)+P(I,J-1))/(DELZ**2)
C     RNSAV(J)=RNSAV(J)+RNS(I,J)
C 113   CONTINUE
C       RNS(I,1)=RNS(I,2)
C       RNS(I,NZ)=RNS(I,NZM1)
C THIS IS DONE SINCE FITTING A QUADRATIC TO P GIVES THIS RESULT
C     RNSAV(1)=RNSAV(1)+RNS(I,1)
C     RNSAV(NL)=RNSAV(NL)+RNS(I,NL)
C 112  CONTINUE
C     DO 9991 J=1,NL
C     DO 9992 I=1,NX
C      RNS(I,J)=RNS(I,J)-RNSAV(J)/NX
C9992 CONTINUE
C     RNS(NXP1,J)=RNS(1,J)
C9991 CONTINUE
C             PRINT*, ' STATIC STABILITY '
C             CALL OUTP(RNS,1,9)
C     CALCULATE  THE  TEMPERATURE
C
      CT = THETA0/G
      DO 111 I=1,NXP1
      T(I,1) = CT*PZB(I)-THETA0
C      PRINT*, ' T = ' , T(I,1)
  111 CONTINUE

      F0=1.E-4
      NZM2=NZM1-1
        VORT(1,1)=1./(1.-(VG(2,1)-VG(NX,1))/(F0*2.*DELX))
        VORT(1,NL)=1./(1.-(VG(2,NL)-VG(NX,NL))/(F0*2.*DELX))
         VORTGS(1,NL)=(VG(2,NL)-VG(NX,NL))/(2.*DELX)
         VORTGS(1,1)=(VG(2,1)-VG(NX,1))/(2.*DELX)
         VORTGS(NXP1,1)=VORTGS(1,1)
         VORTGS(NXP1,NL)=VORTGS(1,NL)
        VORT(NXP1,1)=VORT(1,1)
        VORT(NXP1,NL)=VORT(1,NL)
        DO 31 I=2,NX
          VORT(I,1)=1./(1.-(VG(I+1,1)-VG(I-1,1))/(F0*2.*DELX))
          VORT(I,NL)=1./(1.-(VG(I+1,NL)-VG(I-1,NL))/(F0*2.*DELX))
          VORTGS(I,1)=(VG(I+1,1)-VG(I-1,1))/(2.*DELX)
          VORTGS(I,NL)=(VG(I+1,NL)-VG(I-1,NL))/(2.*DELX)
 31   CONTINUE
      DO 41 J=2,NL-1
        VORT(1,J)=1./(1.-(VG(2,J)-VG(NX,J))/(F0*2.*DELX))
        VORT(NXP1,J)=VORT(1,J)
        VORTGS(1,J)=(VG(2,J)-VG(NX,J))/(2.*DELX)
        VORTGS(NXP1,J)=VORTGS(1,J)
        DO 41 I=2,NX
          VORT(I,J)=1./(1.-(VG(I+1,J)-VG(I-1,J))/(F0*2.*DELX))
          VORTGS(I,J)=(VG(I+1,J)-VG(I-1,J))/(2.*DELX)
 41   CONTINUE



        DO 2 J=1,NL
      II=1
        XREG=XA
 999    IF(XREG.LT.XREAL(1,J)) THEN
C          PRINT*, ' II= ' , II
          DELXR=DELX+XB-XREAL(NX,J)-(XREG-XREAL(1,J))
          WT1=(XREAL(1,J)-XREG)/DELXR
          WT2=1.-WT1
          VORTR(II,J)=WT1*VORT(NX,J)+WT2*VORT(1,J)
      IF(J.EQ.1) THEN
          TR(II,J)=WT1*T(NX,J)+WT2*T(1,J)
      ENDIF
          RNSR(II,J)=WT1*Q(NX,J)+WT2*Q(1,J)
          II=II+1
        ENDIF
        XREG=XREG+DELX
        IF(XREG.LT.XREAL(1,J)) GOTO 999
        IR=0
        DO 1 I=II,NXM1
        XREG=XA+(I-1)*DELX
         IF(XREG.GT.XREAL(NX,J)) GOTO 888
        DO 3 IR=1,NXM1
        IF(XREG.LT.XREAL(IR+1,J)) THEN
          DELXR=XREAL(IR+1,J)-XREAL(IR,J)
          WT1=(XREAL(IR+1,J)-XREG)/DELXR
          WT2=1.-WT1
          VORTR(I,J)=WT1*VORT(IR,J)+WT2*VORT(IR+1,J)
      IF(J.EQ.1) THEN
          TR(I,J)=WT1*T(IR,J)+WT2*T(IR+1,J)
      ENDIF
          RNSR(I,J)=WT1*Q(IR,J)+WT2*Q(IR+1,J)
           GOTO 1
          ENDIF
 3    CONTINUE
 1    CONTINUE
      I=NX
 888  VORTR(NXP1,J)=VORTR(1,J)
      IF(J.EQ.1) THEN
      TR(NXP1,J)=TR(1,J)
      ENDIF
      RNSR(NXP1,J)=RNSR(1,J)
        DO 889 IN=I,NX
        XREG=XA+(IN-1)*DELX
        IF(XREG.LT.XREAL(NX,J)) THEN
          DELXR=XREAL(NX,J)-XREAL(NXM1,J)
          WT1=(XREAL(NX,J)-XREG)/DELXR
          WT2=1.-WT1
          VORTR(IN,J)=WT1*VORT(NXM1,J)+WT2*VORT(NX,J)
      IF(J.EQ.1) THEN
          TR(IN,J)=WT1*T(NXM1,J)+WT2*T(NX,J)
      ENDIF
          RNSR(IN,J)=WT1*Q(NXM1,J)+WT2*Q(NX,J)
        ELSE
          DELXR=DELX+XB-XREAL(NX,J)-(XA-XREAL(1,J))
          WT1=(XREG-XREAL(NX,J))/DELXR
          WT2=1.-WT1
          VORTR(IN,J)=WT1*VORT(1,J)+WT2*VORT(NX,J)   
      IF(J.EQ.1) THEN
          TR(IN,J)=WT1*T(1,J)+WT2*T(NX,J)
      ENDIF
          RNSR(IN,J)=WT1*Q(1,J)+WT2*Q(NX,J)
        ENDIF
 889   CONTINUE
 2    CONTINUE
      DO 13 J=1,NL
        DO 14 I=1,NXP1

           VORTR(I,J)=((VORTR(I,J)*F0)-F0)/F0
           VORT(I,J)=((VORT(I,J)*F0)-F0)/F0
           RNSR(I,J)=RNSR(I,J)*1.E6
 14     CONTINUE
 13   CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C TIME SERIES OF SURFACE RELATIVE VORTICITY C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      IF(NT.LE.NL) THEN
      DO 15 I=1,NXP1
         VORTTS(I,NT)=VORT(I,1)
         TS(I,NT)=TR(I,1)
c       PRINT*, ' VORTTS = ' ,  VORTTS(I,NT)
         RNSTS(I,NT)=VORTGS(I,1)/F0
         VGS(I,NT)=VG(I,1)
 15   CONTINUE
      ENDIF
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      RETURN
      END
      SUBROUTINE OUTP(A,DAY,IOUT)
      INTEGER NP,NL,IOUT,I,J,IMA,IMI,JMA,JMI,A1,JJ
      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION A(NP,NL),RFAC(12),
     +RTAPE,RDAY,DAY,AMA,AMI

      COMMON/TAPE/RTAPE


c      DIMENSION A(NP,NL)
c      DIMENSION RFAC(12)
CCCCCCCCCCCCCCCCCCCCC
C SCALINGS       C
CCCCCCCCCCCCCCCCCCCC
      print*, day,iout
C
C PV:
      RFAC(1)=1.E6
C VG:
      RFAC(2)=1.
C THETA:
      RFAC(3)=1.
C W:
      RFAC(4)=100.
C VORTICITY:
      RFAC(5)=1.E4
C PSI:
      RFAC(6)=1.E-3
      RFAC(7)=1.E13
      RFAC(8)=1.0
      RFAC(9)=1.
	RFAC(10)=1.
	RFAC(11)=1.
	RFAC(12)=1.
C      REWIND(77)
      RDAY=DAY
C     WRITE(77) RDAY,RTAPE
      PRINT*, ' DAY = ' ,DAY, ' RTAPE = ' ,RTAPE
C     WRITE(77,'(8(1P,E10.3))') A
	IF(IOUT.EQ.1) THEN
c------------------------------------alter for write
c       DO 787 J=1,NL
c	JJ=(J-1)*257+1
c       WRITE(16,'(257E10.3)') (A(I,J),I=1,257)
c 787	CONTINUE
        DO 787 j = 1,33
        DO 1787 i = 1, 257
        write (16,'(E10.3)' ) a(i,j)
 1787   continue
 787    continue
       ENDIF 
	IF(IOUT.EQ.2) THEN
C	DO 788 J=1,NL
C	JJ=(J-1)*257+1
C        WRITE(17,'(257E10.3)') (A(I,J),I=1,257)
C 788	CONTINUE
        count=0
        DO 788 j = 1,33
        DO  1788 i =  1, 257
        count=count+1
         write (17,'(E10.3)' ) a(i,j)
 1788 continue
 788  continue
       write (999,*) count

       ENDIF 
 	IF(IOUT.EQ.3) THEN
c	DO 789 J=1,NL
c	JJ=(J-1)*257+1
c       WRITE(18,'(257E10.3)') (A(I,J),I=1,257)
c 789	CONTINUE
        DO 789 j = 1,33
        DO  1789 i =  1, 257
        count=count+1
         write (18,'(E10.3)' ) a(i,j)
 1789 continue
 789  continue

       ENDIF 
	IF(IOUT.EQ.4) THEN
c	DO 790 J=1,NL
c	JJ=(J-1)*257+1
c        WRITE(19,'(257E10.3)') (A(I,J),I=1,257)
c 790	CONTINUE
       DO 790 j = 1,33
       DO 1790 i = 1,257
             write (19, '(E10.3)' ) a(i,j)
 1790       continue
 790        continue

       ENDIF 
 	IF(IOUT.EQ.5) THEN

c	DO 791 J=1,NL
c	JJ=(J-1)*257+1
c        WRITE(20,'(257E10.3)') (A(I,J),I=1,257)
c 791	CONTINUE
       DO 791 j = 1,33
       DO 1791 i = 1,257
             write (20, '(E10.3)' ) a(i,j)
 1791     continue
 791   continue

       ENDIF 
 	IF(IOUT.EQ.8) THEN
	DO 792 J=1,NL
	JJ=(J-1)*257+1
       print*, ' here in outp ', a(i,j)
       WRITE(21,'(257E10.3)') (A(I,J),I=1,257)
 792	CONTINUE
       ENDIF 
 	IF(IOUT.EQ.9) THEN
	DO 799 J=1,NL
	JJ=(J-1)*257+1
       WRITE(23,'(257E10.3)') (A(I,J),I=1,257)
 799	CONTINUE
       ENDIF 
 	IF(IOUT.EQ.10) THEN
	DO 793 J=1,NL
	JJ=(J-1)*257+1
       WRITE(22,'(257E10.3)') (A(I,J),I=1,257)
 793	CONTINUE
       ENDIF 
 	IF(IOUT.EQ.11) THEN
	DO 780 J=1,NL
	JJ=(J-1)*257+1
       WRITE(24,'(257E10.3)') (A(I,J),I=1,257)
 780	CONTINUE
       ENDIF 
 	IF(IOUT.EQ.12) THEN
	DO 781 J=1,NL
	JJ=(J-1)*257+1
       WRITE(25,'(257E10.3)') (A(I,J),I=1,257)
 781	CONTINUE
       ENDIF 

      IF(IOUT.EQ.8.OR.IOUT.EQ.9) THEN
C     WRITE(78) RDAY,RTAPE
C     WRITE(78,'(8(1P,E10.3))') A
      ENDIF
      AMA=-1.0E30
      AMI=1.0E30
      IMA=2000
      IMI=2000
      JMA=2000
      JMI=2000
      DO 1 I=1,NP
      DO 1 J=1,NL
      IF(A(I,J).GT.AMA) THEN
      AMA=A(I,J)
      IMA=I
      JMA=J
      ENDIF
      IF(A(I,J).LT.AMI) THEN
      AMI=A(I,J)
      IMI=I
      JMI=J
      ENDIF
    1 CONTINUE
      IF(ABS(AMA).GT.ABS(AMI))THEN
      A1=AMA
      ELSE
      A1=AMI
      ENDIF
      PRINT *,A1
      PRINT*, ' MAX = ' , AMA, ' MIN = ' , AMI
      PRINT*, ' I,J MAX  = ' , IMA,JMA
      PRINT*, ' I,J MIN ' , IMI,JMI
      PRINT*, ' IOUT = ' , IOUT
C      RAT=AMA/AMI
C      PRINT*, ' RATIO = ' , RAT
      DO 3 J=NL,1,-1
      IF(IOUT.EQ.100) THEN
      PRINT 4,((A(I,J)/(A1+0.0000000001)),I=1,NP,12)
      ELSE
      PRINT 4,((A(I,J)*RFAC(IOUT)),I=1,NP,8)
      IF(J.EQ.1) THEN
      PRINT*, ' SFC VALUES '
      PRINT 4,((A(I,J)*RFAC(IOUT)),I=1,NP)
      ENDIF
      ENDIF
    3 CONTINUE
  4    FORMAT(1X,20(F6.1))
      RETURN
      END
      SUBROUTINE MOVIE(A,DAY,IOUT)
C      IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,IOUT,I,J,IMA,IMI,JMA,JMI,A1,MX,NX,NZ,NXP1,NZM1,
     +W0,J2,I2,I1,IGRID(256*256),FRAMIN(256/4)


      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION A(NP,NL),RFAC(9),
     +RTAPE,RDAY,DAY,AMA,AMI,DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +GRID(NP,256),Z1,Z2,WT1,WT2,DELZG,HG,X1,X2,DELXG,GRID1(256*256),
     +GMAX,GMIN,GRAD,CUT,GRID2(256,256)

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS


  
      CHARACTER*1 FRAME(256*256)


c      INTEGER FRAMIN(256/4)

c      DIMENSION A(NP,NL),GRID(NP,256),GRID1(256*256),GRID2(256,256)
c      DIMENSION RFAC(9),IGRID(256*256)   
      EQUIVALENCE (FRAME(1),FRAMIN(1))


CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
CC INTERPOLATING GRID TO 256X256 GRID FOR MOVIETOOL GRAPHICS CC
CC NB THIS IS FOR NL AND NP LESS THAN 256  AT THE MOMENT     CC
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C  SET TOP AND BOTTOM C
CC    
C      OPEN(10,FILE='pvgrid1',FORM='UNFORMATTED')
C      OPEN(UNIT=15,FILE='pvgrid1',ACCESS='APPEND',
C     :     FORM='UNFORMATTED',STATUS='NEW')
      DO 101 I=1,NP
      GRID(I,1)=A(I,1)
      GRID(I,256)=A(I,NZ)
      Z1=0.
      Z2=DELZ
      J2=2
      DELZG=(NL-1)*DELZ/255.
C      PRINT*, ' A = ' ,A(I,1)
      DO 100 J=2,255
        HG=(J-1)*DELZG
        IF(HG.GT.Z2) THEN
          J2=J2+1
          Z1=Z1+DELZ
          Z2=Z2+DELZ
C      PRINT*, ' A = ' , A(I,J2-1)
        ENDIF
        WT1=(Z2-HG)/DELZ
        WT2=(HG-Z1)/DELZ
        GRID(I,J)=WT1*A(I,J2-1)+WT2*A(I,J2)
C      PRINT*, ' GRID = ' , GRID(I,J)
 100  CONTINUE
 101  CONTINUE
      DO 102 J=1,256
C      PRINT*, ' J = ' ,J
       GRID2(1,J)=GRID(1,J)
       GRID2(256,J)=GRID(NP,J)
      X1=0.
      X2=DELX
      I2=2
      DELXG=(NP-1)*DELX/255.
C      PRINT*, ' A = ' ,A(1,J)
      DO 103 I=2,255
        HG=(I-1)*DELXG
        IF(HG.GT.X2) THEN
          I2=I2+1
          X1=X1+DELX
          X2=X2+DELX
C        PRINT*, ' A = ' , A(I2-1,J)
        ENDIF
        WT1=(X2-HG)/DELX
        WT2=(HG-X1)/DELX
        GRID2(I,J)=WT1*GRID(I2-1,J)+WT2*GRID(I2,J)
C C     PRINT*, ' GRID = ' , GRID(I,J)
 103  CONTINUE
 102  CONTINUE
      I1=0   
      GMIN=10000000.
      GMAX=-1000000.
      DO 1000 J=1,256
      DO 1001 I=1,256
      IF(GRID2(I,J).GT.GMAX) GMAX=GRID2(I,J)
      IF(GRID2(I,J).LT.GMIN) GMIN=GRID2(I,J)
      I1=I1+1
      GRID1(I1)=GRID2(I,J)
 1001 CONTINUE
 1000 CONTINUE
CCCCCCCCCCCCCCCCC CALC GRADIENT FOR DIGITALISING DATA
      GRAD=(GMAX-GMIN)/249.
      CUT=GMIN-GRAD
      PRINT*, ' GMAX,GMIN,CUT,GRAD', GMAX,GMIN,CUT,GRAD
      DO 900 I=1,256*256
      IGRID(I)=NINT((GRID1(I)-CUT)/GRAD)
      FRAME(I)=CHAR(IFIX(IGRID(I)))
C        PRINT*, ' IGRID = ' , IGRID(I)
 900  CONTINUE
Cc      DO 901 I=1,256*256,256
C      PRINT*, ' IGRID= ' , IGRID(I)
C 901  CONTINUE 
       DO 333 J=1,256
       I1=(J-1)*256+1
       PRINT 43,(IGRID(I),I=I1,I1+256,20)
 333    CONTINUE
  43    FORMAT(1X,20(I4))
       WRITE(15) FRAMIN
C
CC
C
C
CCCCCCCCCCCCCCCCCC
C SCALINGS       C
CCCCCCCCCCCCCCCCCC
C
C PV:
      RFAC(1)=1.E6
C VG:
      RFAC(2)=1.
C THETA:
      RFAC(3)=1.
C W:
      RFAC(4)=100.
C VORTICITY:
      RFAC(5)=1.E4
C PSI:
      RFAC(6)=1.E-3
      RFAC(7)=1.E13
      RFAC(8)=1.0
      RFAC(9)=1.
C      REWIND(77)
      RTAPE=23.1
      RDAY=DAY
C     WRITE(77) RDAY,RTAPE
      PRINT*, ' DAY = ' ,DAY, ' RTAPE = ' ,RTAPE
C     WRITE(77,'(8(1P,E10.3))') A
      IF(IOUT.EQ.8.OR.IOUT.EQ.9) THEN
C     WRITE(78) RDAY,RTAPE
C     WRITE(78,'(8(1P,E10.3))') A
      ENDIF
      AMA=-1.0E30
      AMI=1.0E30
      IMA=2000
      IMI=2000
      JMA=2000
      JMI=2000
      DO 1 I=1,NP
      DO 1 J=1,NL
      IF(A(I,J).GT.AMA) THEN
      AMA=A(I,J)
      IMA=I
      JMA=J
      ENDIF
      IF(A(I,J).LT.AMI) THEN
      AMI=A(I,J)
      IMI=I
      JMI=J
      ENDIF
    1 CONTINUE
      IF(ABS(AMA).GT.ABS(AMI))THEN
      A1=AMA
      ELSE
      A1=AMI
      ENDIF
      PRINT *,A1
      PRINT*, ' MAX = ' , AMA, ' MIN = ' , AMI
      PRINT*, ' I,J MAX  = ' , IMA,JMA
      PRINT*, ' I,J MIN ' , IMI,JMI
      PRINT*, ' IOUT = ' , IOUT
C      RAT=AMA/AMI
C      PRINT*, ' RATIO = ' , RAT
      DO 3 J=NL,1,-1
      IF(IOUT.EQ.100) THEN
      PRINT 4,((A(I,J)/(A1+0.0000000001)),I=1,NP,12)
      ELSE
      PRINT 4,((A(I,J)*RFAC(IOUT)),I=1,NP,8)
      ENDIF
    3 CONTINUE
  4    FORMAT(1X,20(F6.1))
      RETURN
      END
      SUBROUTINE OUT(KOUNT,DAY,MDUMM)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +NZM2,J,I,KOUNT,DAY,MDUMM,JJ

      PARAMETER(NP=257,NL=33)

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,DFB(NP),DFT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),
     +Q1(NP,NL),RCOND(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),
     +PZB1(NP),PZB2(NP),PZT1(NP),PZT2(NP),
     +VGG(NP,NL),VORT(NP,NL),CT,RTAPE,F0,QSUM,QAV,
     +VORTTS(NP,NL),RNSTS(NP,NL),UAGEO(NP,NL),VGZ,PSIZ,VORTGS(NP,NL),
     :VGS(NP,NL),TS(NP,NL)

      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/STOREQW/Q1,RCOND,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/STORE/PZB1,PZB2,PZT1,PZT2
      COMMON/GAUSS2/AL,T
      COMMON/VORT2/VORTTS,RNSTS,VORTGS,
     :VGS,TS

	COMMON/UAGEOS/UAGEO
c      DIMENSION VORT(NP,NL),VGG(NP,NL)
    
C
C-----------------------------------------------------------------------
C
C  WRITE THE FIELDS IN UNFORMATED FORM ONTO MAGNETIC TAPE
C
C-----------------------------------------------------------------------
C
      CT=THETA0/G
      NZM2=NZM1-1
C     WRITE (7)KOUNT,DAY,MDUMM
C     WRITE (7)PZB1
C     WRITE (7)PZT1
C     WRITE (7)DFB
C     WRITE (7)DFT
C     WRITE (7)P
C     WRITE (7)Q
      DO 345 J=1,NL
      VGG(1,J)=(P(2,J)-P(NX,J))/(2.0E-04*DELX)
      VGG(NXP1,J)=VGG(1,J)
      DO 345 I=2,(NP-1)
      VGG(I,J)=(P(I+1,J)-P(I-1,J))/(2.0E-4*DELX)
C      IF(J.EQ.1) PRINT*, ' V = ' , VGG(I,J)
 345  CONTINUE
      RTAPE=23.1
c-----------------------------------------VG change for write 	
c	DO 787 J=1,NL
c	JJ=(J-1)*257+1
c       WRITE(15,'(257E10.3)') (VGG(I,J),I=1,257)
c 787	CONTINUE
      DO 787 J = 1,33
      DO 1787 i = 1, 257
         write (15,'(E10.3)') VGG(i,j)
 1787    continue
 787     continue

C     WRITE(77) DAY,RTAPE,DELX,DELZ,XA
C     WRITE(77,'(8(1P,E10.3))') VGG
      PRINT *,'POTENTIAL VORTICITY'
C     CALL MOVIE(Q,DAY,1)
      CALL OUTP(Q,DAY,1)
      PRINT *,'VG'
      CALL OUTP(VGG,DAY,2)
      PRINT*, ' TEMPERATURE '
C      DO 346 I=1,NP
C        DO 346 J=1,NL
C            T(I,J)=T(I,J)-THETZM*(J-1)*DELZ
C
C       W(I,J)=W(I,J)*T(I,J)
C 346  CONTINUE
      CALL OUTP(T,DAY,3)
      PRINT*, ' VERTICAL VELOCITY'
      CALL OUTP(W,DAY,4)
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C VORTICITY CALCULATION             C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
      F0=1.E-4
      QSUM=0.
      DO 445 J=1,NL
        VORT(1,J)=1./(1.-(VGG(2,J)-VGG(NX,J))/(F0*2.*DELX))
C       VORT(1,J)=VORT(1,J)-F0
        QSUM=QSUM+(Q(1,J)/VORT(1,J))
        VORT(NXP1,J)=VORT(1,J)
        DO 445 I=2,NX
          VORT(I,J)=1./(1.-(VGG(I+1,J)-VGG(I-1,J))/(F0*2.*DELX))
          QSUM=QSUM+(Q(I,J)/VORT(I,J))
C         VORT(I,J)=VORT(I,J)-F0
 445  CONTINUE
      QAV=QSUM/(NX*NL)
C      PRINT*, 'AVERAGE Q = ' , QAV
      DO 999 J=1,NL
      DO 998 I=1,NXP1
      VORT(I,J)=VORT(I,J)*F0-F0
 998  CONTINUE
 999  CONTINUE
      PRINT*, ' RELATIVE VORTICITY'
      CALL OUTP(VORT,DAY,5)
      PRINT *,'PSI'
       CALL OUTP(PSI,DAY,6)

CCCCCCCCCCCCCCC UAG CALCULATION CCCCCCCCCCCCCCCCCCCC
C
C     UAG = -r W VGZ /F      +     PSIZ             
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	DO 981 J=2,NZM1
	  DO 982 I=1,NXP1
	   	VGZ=(VGG(I,J+1)-VGG(I,J-1))/(2.*DELZ)
		PSIZ=(PSI(I,J+1)-PSI(I,J-1))/(2.*DELZ)
		UAGEO(I,J)=-(W(I,J)*VGZ/F)+(PSIZ/RO(J))
 982	  CONTINUE
 981	CONTINUE
	DO 983 I=1,NP
	UAGEO(I,1)=(PSI(I,2)-PSI(I,1))/(RO(1)*DELZ)
	UAGEO(I,NL)=(PSI(I,NZ)-PSI(I,NZM1))/(RO(NZ)*DELZ)
 983	CONTINUE
	PRINT*, ' AGEOSTROPHIC U '
	CALL OUTP(UAGEO,DAY,10)
C     WRITE (7)PSI
C     WRITE (7)W
C      WRITE (7)F1
C      WRITE (7)FT1
C
C
C
C
C
C-----------------------------------------------------------------------
C
      RETURN
      END
 
      SUBROUTINE PVSUM2(DAY)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +J,I
 
      PARAMETER(NP=257,NL=33) 

      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,DFB(NP),DFT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),
     +Q1(NP,NL),RCOND(NP,NL),W(NP,NL),VM1(NP,NL),RQMT(NP,NL),RQT(NP,NL),
     +RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),
     +PZB1(NP),PZB2(NP),PZT1(NP),PZT2(NP),PV1,PV2,DPZT,
     +VGG(NP,NL),VORT(NP,NL),F0,ROSUM,QSUM,RQSUM,ROAV,QAV,RQAV,
     +VSUM,VAV,UAGEO(NP,NL),USUM,UAV,WSUM,WAV,DAY,VGZ,PSIZ,UB,UUSUM,
     +UUAV
      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/F,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM
     

      COMMON/INTS/DFB,DFT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/STOREQW/Q1,RCOND,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/STORE/PZB1,PZB2,PZT1,PZT2
      COMMON/GAUSS2/AL,T
      COMMON/PVADJ/PV1,PV2,DPZT
	COMMON/UAGEOS/UAGEO

c     DIMENSION VGG(NP,NL),VORT(NP,NL)
      DO 345 J=1,NL
      VGG(1,J)=(P(2,J)-P(NX,J))/(2.0E-04*DELX)
      VGG(NXP1,J)=VGG(1,J)
      DO 345 I=2,(NP-1)
      VGG(I,J)=(P(I+1,J)-P(I-1,J))/(2.0E-4*DELX)
 345  CONTINUE
CCCCCCCCCCCCCCC UAG CALCULATION CCCCCCCCCCCCCCCCCCCC
C
C     UAG = - W VGZ /F      +     PSIZ/r             
C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC

	DO 981 J=2,NZM1
	  DO 982 I=1,NXP1
	   	VGZ=(VGG(I,J+1)-VGG(I,J-1))/(2.*DELZ)
		PSIZ=(PSI(I,J+1)-PSI(I,J-1))/(2.*DELZ)
		UAGEO(I,J)=-(W(I,J)*VGZ/F)+(PSIZ/RO(J))
 982	  CONTINUE
 981	CONTINUE
	DO 983 I=1,NXP1
	UAGEO(I,1)=(PSI(I,2)-PSI(I,1))/(DELZ*RO(1))
	UAGEO(I,NL)=(PSI(I,NZ)-PSI(I,NZM1))/(RO(NZ)*DELZ)
 983	CONTINUE

      F0=1.E-4
      ROSUM=0.
      QSUM=0.
      RQSUM=0.
      VSUM=0.
      USUM=0.
      UUSUM=0.
      WSUM=0.
        VORT(1,1)=1./(1.-(VGG(2,1)-VGG(NX,1))/(F0*2.*DELX))
        VORT(1,NL)=1./(1.-(VGG(2,NL)-VGG(NX,NL))/(F0*2.*DELX))
          RQSUM=RQSUM+(Q(1,1)/VORT(1,1))*RO(1)*0.5
          RQSUM=RQSUM+(Q(1,NL)/VORT(1,NL))*RO(NL)*0.5
          VSUM=VSUM+(VGG(1,1)*VGG(1,1)/VORT(1,1))*RO(1)*0.5
          VSUM=VSUM+(VGG(1,NL)*VGG(1,NL)/VORT(1,NL))*RO(NL)*0.5
          USUM=USUM+(UAGEO(1,1)*UAGEO(1,1)/VORT(1,1))*RO(1)*0.5
          USUM=USUM+(UAGEO(1,NL)*UAGEO(1,NL)/VORT(1,NL))*RO(NL)*0.5
          WSUM=WSUM+(W(1,1)*W(1,1)/VORT(1,1))*RO(1)*0.5
          WSUM=WSUM+(W(1,NL)*W(1,NL)/VORT(1,NL))*RO(NL)*0.5
      UB=0.
          UUSUM=UUSUM+(UAGEO(1,1)*UB/VORT(1,1))*RO(1)*0.5
      UB=SH*(NZ-1)*DELZ
          UUSUM=UUSUM+(UAGEO(1,NL)*UB/VORT(1,NL))*RO(NL)*0.5

      QSUM=QSUM+0.5*(Q(1,1)/VORT(1,1))
      QSUM=QSUM+0.5*(Q(1,NL)/VORT(1,NL))
        DO 444 I=2,NX
          VORT(I,1)=1./(1.-(VGG(I+1,1)-VGG(I-1,1))/(F0*2.*DELX))
          VORT(I,NL)=1./(1.-(VGG(I+1,NL)-VGG(I-1,NL))/(F0*2.*DELX))
          RQSUM=RQSUM+(Q(I,1)/VORT(I,1))*RO(1)*0.5
          RQSUM=RQSUM+(Q(I,NL)/VORT(I,NL))*RO(NL)*0.5
          VSUM=VSUM+(VGG(I,1)*VGG(I,1)/VORT(I,1))*RO(1)*0.5
          VSUM=VSUM+(VGG(I,NL)*VGG(I,NL)/VORT(I,NL))*RO(NL)*0.5
          USUM=USUM+(UAGEO(I,1)*UAGEO(I,1)/VORT(I,1))*RO(1)*0.5
          USUM=USUM+(UAGEO(I,NL)*UAGEO(I,NL)/VORT(I,NL))*RO(NL)*0.5
          WSUM=WSUM+(W(I,1)*W(I,1)/VORT(I,1))*RO(1)*0.5
          WSUM=WSUM+(W(I,NL)*W(I,NL)/VORT(I,NL))*RO(NL)*0.5
      UB=0.
          UUSUM=UUSUM+(UAGEO(I,1)*UB/VORT(I,1))*RO(1)*0.5

      UB=SH*(NZ-1)*DELZ

          UUSUM=UUSUM+(UAGEO(I,NL)*UB/VORT(I,NL))*RO(NL)*0.5
      QSUM=QSUM+0.5*(Q(I,1)/VORT(I,1))
      QSUM=QSUM+0.5*(Q(I,NL)/VORT(I,NL))
 444  CONTINUE
      DO 445 J=2,NL-1
        VORT(1,J)=1./(1.-(VGG(2,J)-VGG(NX,J))/(F0*2.*DELX))
      ROSUM=ROSUM+RO(J)
        VORT(NXP1,J)=VORT(1,J)
      RQSUM=RQSUM+(Q(1,J)/VORT(1,J))*RO(J)
      VSUM=VSUM+(VGG(1,J)*VGG(1,J)/VORT(1,J))*RO(J)  
      USUM=USUM+(UAGEO(1,J)*UAGEO(1,J)/VORT(1,J))*RO(J)   
      WSUM=WSUM+(W(1,J)*W(1,J)/VORT(1,J))*RO(J)
      UB=SH*(J-1)*DELZ
      UUSUM=UUSUM+(UAGEO(1,J)*UB/VORT(1,J))*RO(J)

      QSUM=QSUM+(Q(1,J)/VORT(1,J))
        DO 445 I=2,NX
          VORT(I,J)=1./(1.-(VGG(I+1,J)-VGG(I-1,J))/(F0*2.*DELX))
          RQSUM=RQSUM+(Q(I,J)/VORT(I,J))*RO(J)
          VSUM=VSUM+(VGG(I,J)*VGG(I,J)/VORT(I,J))*RO(J)
          USUM=USUM+(UAGEO(I,J)*UAGEO(I,J)/VORT(I,J))*RO(J)
          WSUM=WSUM+(W(I,J)*W(I,J)/VORT(I,J))*RO(J)
          UUSUM=UUSUM+(UAGEO(I,J)*UB/VORT(I,J))*RO(J)
      QSUM=QSUM+(Q(I,J)/VORT(I,J))
 445  CONTINUE
      ROAV=ROSUM/NL
      QAV=QSUM*1.E6/(NX*(NL-1))
      RQAV=RQSUM*1.E6/(NX*(NL-1))  
      VAV=VSUM/(NX*(NL-1))
      UAV=USUM/(NX*(NL-1))
      WAV=WSUM/(NX*(NL-1))
      UUAV=UUSUM*2./(NX*(NL-1))
C      RQAV=RQAV/RO(1)
      PRINT*, '*************************'

      PRINT*, ' DAY = ' , DAY
      PRINT*, '*************************'
      PRINT*, ' AVERAGE  Q = ' , QAV
      PRINT*, ' AVERAGE RQ = ' , RQAV

      PRINT*, ' TOTAL KE = ' , VAV

      PRINT*, '*************************'

      RETURN
      END
      SUBROUTINE FRICT(KOUNT)
C       IMPLICIT UNDEFINED(A-Z)
      INTEGER NP,NL,MX,NX,NZ,NXP1,NZM1,W0,NXM1,
     +NZM2,KOUNT,J,I

      
      PARAMETER(NP=257,NL=33)


      DOUBLE PRECISION DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS(NP),
     +FC,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM,PZB(NP),PZT(NP),
     +RO(NL),ROZ(NL),P(NP,NL),Q(NP,NL),PSI(NP,NL),
     +F1(NP,NL),FT1(NP,NL),AL(NL),T(NP,NL),DFLUQ(NP,NL),
     +V(NP,NL),Q1(NP,NL),UAG(NP,NL),W(NP,NL),VM1(NP,NL),
     +RQMT(NP,NL),RQT(NP,NL),RQPT(NP,NL),RQS(NP,NL),XREAL(NP,NL),
     +TX2(NP),TTH2(NP),ST(NP),VG(NP,NL),NS(NP,NL),CV,CZ,CT,CTZ


      COMMON/CONSTA/MX,NX,NZ,NXP1,NZM1,W0
      COMMON/CONSTB/DELX,DELZ,DELXSQ,DELZSQ,CRHS,AS
      COMMON/CONST1/NXM1
      COMMON/CONST2/FC,THETA0,G,TIMINC,XA,XB,SH,ZL,EP,THETZM

      COMMON/INTS/PZB,PZT,RO,ROZ
      COMMON/INIT/P,Q,PSI
      COMMON/FRIC/F1,FT1
      COMMON/GAUSS2/AL,T
      COMMON/FRIC1/DFLUQ,V
      COMMON/STOREQW/Q1,UAG,W,VM1
     +,RQMT,RQT,RQPT,RQS,XREAL
      COMMON/STOREF/TX2,TTH2,ST
      COMMON/VGEOS/VG
c      DIMENSION NS(NP,NL)
C
      CV = 1.0/(FC*2.0*DELX)
      CZ=0.5/DELZ
      NZM2 = NZM1 - 1
C
C          CALCULATE  THE  LONG-FRONT   VELOCITY   V
C                  V(X,Z) = (DPHI(X,Z)/DX)/F
C
      DO 2 J=1,NL
      VG(1,J) = CV*(P(2,J) - P(NX,J))
      DO 1 I=2,NX
      VG(I,J) = CV*(P(I+1,J) - P(I-1,J))
    1 CONTINUE
      VG(NXP1,J)=VG(1,J)
 2    CONTINUE
C
C     CALCULATE  THE  TEMPERATURE
C
      CT = THETA0/G
      CTZ= CT/(2.0*DELZ)
      DO 111 I=1,NXP1
      T(I,1) = CT*PZB(I)
      DO 222 J=2,NZM1
      T(I,J)=(P(I,J+1)-P(I,J-1))*CTZ
 222  CONTINUE
      T(I,NZ)=CT*PZT(I)
  111 CONTINUE
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C CALCULATE N-SQUARED          C
CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC
C      PRINT*, ' TEMP1 '
C      CALL OUTP(T)
       DO 112 I=1,NXP1
         DO 113 J=2,NZM1
           NS(I,J)=(T(I,J+1)-T(I,J-1))/(2.*DELZ)
      NS(I,J)=NS(I,J)/CT
C      NS(I,J)=NS(I,J)-THETZM
C      IF(I.EQ.1) PRINT*, ' LEVEL: ' , J , '  ' , NS(I,J)
  113   CONTINUE
        NS(I,1)=2.*NS(I,2)-NS(I,3)
       NS(I,NZ)=NS(I,NZM1)
  112  CONTINUE
C           PRINT*, ' STATIC STABILITY '
C          CALL OUTP(NS)
C
C      PRINT*, ' PHI FIELD '
C      CALL OUTP(P)
C
C      PRINT*, ' VGG = '
C      CALL OUTP(V)
C      PRINT*, ' THETA'
C      CALL OUTP(T)
C      GOTO 999
C        CALCULATE  THE  CROSS  FRONT  VELOCITY  U(X,Z)
C
      CZ = .5/DELZ
      RETURN
C
      END
C     FUNCTION DELAY(N)v
C     DELAY=UNIT(N)
C     IF(DELAY)1,1,2
C   1 RETURN
C   2 PRINT 100,NZ
C 100 FORMAT(' PARITY  ERROR  ON  CHANNEL  ',I4)
C
C     RETURN
C     END
      SUBROUTINE FFT99A(GRID,WORK,TRIGS,IFAX,INC,JUMP,N,LOT,ISIGN)
C       IMPLICIT UNDEFINED(A-Z)
      DOUBLE PRECISION GRID(*),WORK(*),TRIGS(*),IFAX(*),RN

      INTEGER N2,INC,JUMP,N,LOT,ISIGN,L,IOFF,I

c      DIMENSION GRID(*),WORK(*),TRIGS(*),IFAX(*)
      N2=N/2
      IF(ISIGN.EQ.-1)THEN
      RN=1./N
      DO 10 L=1,LOT
10    CALL REALFT(GRID(1+(L-1)*JUMP),N2,1)
      DO 20 L=1,LOT
      IOFF=(L-1)*JUMP
      DO 20 I=1,N
20      GRID(I+IOFF)=GRID(I+IOFF)*RN
      DO 25 L=1,LOT
      IOFF=(L-1)*JUMP
       DO 25 I=2,N,2
25      GRID(I+IOFF)=-GRID(I+IOFF)
      ELSE
      DO 30 L=1,LOT
      IOFF=(L-1)*JUMP
      DO 30 I=1,N
30      GRID(I+IOFF)=2.*GRID(I+IOFF)
      DO 35 L=1,LOT
      IOFF=(L-1)*JUMP
      DO 35 I=2,N,2
35      GRID(I+IOFF)=-GRID(I+IOFF)
      DO 40 L=1,LOT
40      CALL REALFT(GRID(1+(L-1)*JUMP),N2,-1)
      ENDIF
      RETURN
      END
      SUBROUTINE REALFT(DATA,N,ISIGN)
C       IMPLICIT UNDEFINED(A-Z)
      DOUBLE PRECISION DATA(*),WR,WI,WPR,WTEMP,THETA,C1,C2,
     +H1R,H1I,H2R,H2I,N2P3,WPI

      INTEGER N,ISIGN,I,I1,I2,I3,I4

      REAL WRS,WIS

C      REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
c      DIMENSION DATA(*)
      THETA=6.28318530717959D0/2.0D0/DBLE(N)
      C1=0.5
      IF (ISIGN.EQ.1) THEN
      C2=-0.5
      CALL FOUR1(DATA,N,+1)
      ELSE
      C2=0.5
      THETA=-THETA
      ENDIF
      WPR=-2.0D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      WR=1.0D0+WPR
      WI=WPI
      N2P3=2*N+3
      DO 11 I=2,N/2+1
      I1=2*I-1
      I2=I1+1
      I3=N2P3-I2
      I4=I3+1
      WRS=SNGL(WR)
      WIS=SNGL(WI)
      H1R=C1*(DATA(I1)+DATA(I3))
      H1I=C1*(DATA(I2)-DATA(I4))
      H2R=-C2*(DATA(I2)+DATA(I4))
      H2I=C2*(DATA(I1)-DATA(I3))
      DATA(I1)=H1R+WRS*H2R-WIS*H2I
      DATA(I2)=H1I+WRS*H2I+WIS*H2R
      DATA(I3)=H1R-WRS*H2R+WIS*H2I
      DATA(I4)=-H1I+WRS*H2I+WIS*H2R
      WTEMP=WR
      WR=WR*WPR-WI*WPI+WR
      WI=WI*WPR+WTEMP*WPI+WI
11    CONTINUE
      IF (ISIGN.EQ.1) THEN
      H1R=DATA(1)
      DATA(1)=H1R+DATA(2)
      DATA(2)=H1R-DATA(2)
      ELSE
      H1R=DATA(1)
      DATA(1)=C1*(H1R+DATA(2))
      DATA(2)=C1*(H1R-DATA(2))
      CALL FOUR1(DATA,N,-1)
      ENDIF
      RETURN
      END
      SUBROUTINE FOUR1(DATA,NN,ISIGN)
      DOUBLE PRECISION DATA(*), WR,WI,WPR,WPI,WTEMP,THETA,TEMPR,
     +TEMPI,WTEMP

      INTEGER N,NN,ISIGN,I,J,M,ISTEP,MMAX

C C     REAL*8 WR,WI,WPR,WPI,WTEMP,THETA
c      DIMENSION DATA(*)
      N=2*NN
      J=1
      DO 11 I=1,N,2
      IF(J.GT.I)THEN
      TEMPR=DATA(J)
      TEMPI=DATA(J+1)
       DATA(J)=DATA(I)
       DATA(J+1)=DATA(I+1)
       DATA(I)=TEMPR
       DATA(I+1)=TEMPI
       ENDIF
       M=N/2
1       IF ((M.GE.2).AND.(J.GT.M)) THEN
        J=J-M
        M=M/2
      GO TO 1
      ENDIF
      J=J+M
11    CONTINUE
      MMAX=2
2     IF (N.GT.MMAX) THEN
      ISTEP=2*MMAX
      THETA=6.28318530717959D0/(ISIGN*MMAX)
      WPR=-2.D0*DSIN(0.5D0*THETA)**2
      WPI=DSIN(THETA)
      WR=1.D0
      WI=0.D0
      DO 13 M=1,MMAX,2
      DO 12 I=M,N,ISTEP
      J=I+MMAX
      TEMPR=SNGL(WR)*DATA(J)-SNGL(WI)*DATA(J+1)
      TEMPI=SNGL(WR)*DATA(J+1)+SNGL(WI)*DATA(J)
      DATA(J)=DATA(I)-TEMPR
      DATA(J+1)=DATA(I+1)-TEMPI
      DATA(I)=DATA(I)+TEMPR
      DATA(I+1)=DATA(I+1)+TEMPI
12        CONTINUE
       WTEMP=WR
       WR=WR*WPR-WI*WPI+WR
       WI=WI*WPR+WTEMP*WPI+WI
13      CONTINUE
       MMAX=ISTEP
      GO TO 2
      ENDIF
      RETURN
      END

	SUBROUTINE RANDOM(RT)
	DOUBLE PRECISION RM,RT
	INTEGER M,IA,IC,IFF,IDUM,J,IR(97)
	PARAMETER (M=714025,IA=1366,IC=150989,RM=1./M)
C        DIMENSION IR(97)
	DATA IFF/0/
	DO 99 II = 1,100
	IF(IDUM.LT.0.OR.IFF.EQ.0) THEN
	IFF=1
	IDUM=MOD(IC-IDUM,M)
	DO 11 J=1,97
		IDUM=MOD(IA*IDUM+IC,M)
		IR(J)=IDUM
 11	CONTINUE
	IDUM=MOD(IA*IDUM+IC,M)
	IY=IDUM
	ENDIF
	J=1+(97*IY)/M
	IF(J.GT.97.OR.J.LT.1) PAUSE
	IY=IR(J)
	RAN2=IY*RM
	IDUM=MOD(IA*IDUM+IC,M)
	IR(J)=IDUM
	RT=IR(J)*RM
 99	CONTINUE
	RETURN
	END
