	PROGRAM HURRICANE
C
C  ***   NRD is the maximum number of radial nodes; NTG is   ***
C  ***      the maximum number of time points in the         ***
C  ***               graphics time series                    ***
C
	PARAMETER(NRD=200, NTG=600)
C
C  ***  NZD is the number of nodes in the vertical used in ***
C  ***     constructing the two-dimensional plots          ***
C  ***  MI and NI are the maximum dimensions of the 2-D plot arrays ***
C
	PARAMETER(NZD=60,MI=100,NI=150,MI2=2*MI+1)
C
C  ***     Dimension arrays of dependent variables   ***
C
	REAL RBS1(NRD), RBS2(NRD), RBS3(NRD), RTS1(NRD), RTS2(NRD) 
	REAL X1(NRD), X2(NRD), X3(NRD), XS1(NRD), XS2(NRD), XS3(NRD)
	REAL XM1(NRD), XM2(NRD), XM3(NRD),RTS3(NRD)
	REAL MU1(NRD), MU2(NRD), MU3(NRD)
C
C  ***     Dimension ocean variables   ***
C
	REAL HMIX(NRD), UHMIX1(NRD), UHMIX2(NRD), UHMIX3(NRD), HM
	REAL HEDDY
C
C  ***      Dimension various diagnostic quantities     ***
C
	REAL P(NRD), PS0(NRD), PS2(NRD), PS3(NRD), GB(NRD), RMS2(NRD)
	REAL RB1(NRD), RB2(NRD), RT1(NRD)
C
C  ***   Dimension VISCOUS and working arrays   ***
C
	REAL VIS(NRD), VIS2(NRD), XVIS(NRD), DV(NRD), AF(NRD)
	REAL Q2(NRD), Q3(NRD), XMVIS(NRD), RMM2(NRD), EPRE(NRD)
C
C  ***      Dimension graphics arrays       ***
C
	REAL RG(NRD), RTG(NRD), RMG(NTG), RTMG(NTG), VBG(NRD)
	REAL RGP(NRD), PG(NRD), HUMG(NRD), HMG(NTG)
	REAL SSTT(NTG), SSTR(NRD), SST1(NRD),SST2(NRD)
	REAL SST3(NRD)
	REAL VTG(NRD), VMG(NTG,2), VTMG(NTG), UG(NRD,2), WG(NRD,2)
	REAL MUG(NTG), MDG(NTG), PCG(NTG,2),TIMEG(NTG), XG(NRD,3)
 	REAL PSIT(NRD,NZD),PRT(NRD,NZD),TTE(NRD,NZD),PPT(NRD,NZD)
 	REAL VVT(NRD,NZD),RRT(NRD,NZD),ZZT(NRD,NZD),XIT(NRD,NZD)
	REAL PBO(NZD),ZBO(NZD),ZTE(NZD),RRTI(NRD,NZD)
	REAL WIT(NRD,NZD),UIT(NRD,NZD),WORKIT(NRD,NZD),WORKITI(MI,NI)
	DIMENSION VVTI(MI,NI),PSITI(MI,NI),UITI(MI,NI),WITI(MI,NI)
	DIMENSION TTEI(MI,NI),PPTI(MI,NI),XITI(MI,NI),RRI(MI)
	DIMENSION RMT(NRD,NZD),RMTI(MI,NI),ZZI(NI),RAIN(MI,NI)
C
C   ***   Hovmuller diagrams
C
	REAL WHOV(MI,NTG),VHOV(MI,NTG),UHOV(MI,NTG)
C
C  ***       Other quantities    ***
C
	REAL H_A, MEQ, MF, MT, MUMAX, MDMIN
	REAL TIME,DTG,ROG,VM,RM,R0,TS,TO,ALAT,TLAND,TSHEAR,VEXT,UT
	REAL EDDYTIME,RWIDE,RO,AHM,PA,CD,CD1,CDCAP,CECD,PNU,TAUR
	REAL RADMAX,TAUC,EFRAC,DSST,GM,HS
	REAL*8 DT,TT,ATIME
	CHARACTER*4 DIM,FMT,OM,SURFACE
	CHARACTER*1 VDISP,RST
C
C   ***  Assign Unit 11 to the Input Parameter File   ***
C
	OPEN(UNIT=11,FILE='hurrparams.txt',STATUS='OLD')
C
C           ***  Read in PARAMETERS  ***
C
	READ(11,*)
	READ(11,*)
	READ(11,*)
	READ(11,*)TIME
	READ(11,*)DTG
	READ(11,*)ROG
	READ(11,*)
	READ(11,*)
	READ(11,*)RST
	READ(11,*)VM
	READ(11,*)RM
	READ(11,*)R0
	READ(11,*)
	READ(11,*)
	READ(11,*)TS
	READ(11,*)TO
	READ(11,*)H_A
	READ(11,*)ALAT
	READ(11,*)TSHEAR
	READ(11,*)VEXT
	READ(11,*)
	READ(11,*)
	READ(11,*)TLAND
	READ(11,*)SURFACE
	READ(11,*)HS
	READ(11,*)
	READ(11,*)
	READ(11,*)OM
	READ(11,*)UT
	READ(11,*)EDDYTIME
	READ(11,*)HEDDY
	READ(11,*)RWIDE
	READ(11,*)
	READ(11,*)
	READ(11,*)DIM
	READ(11,*)FMT
	READ(11,*)
	READ(11,*)
	READ(11,*)NR
	READ(11,*)DT
	READ(11,*)RO
	READ(11,*)
	READ(11,*)
	READ(11,*)AHM
	READ(11,*)AHMCORE
	READ(11,*)RCORE
	READ(11,*)PA
	READ(11,*)CD
	READ(11,*)CD1
	READ(11,*)CDCAP
	READ(11,*)CECD
	READ(11,*)PNU
	READ(11,*)TAUR
	READ(11,*)RADMAX
	READ(11,*)TAUC
	READ(11,*)EFRAC
	READ(11,*)DPB
	READ(11,*)
	READ(11,*)
	READ(11,*)HM
	READ(11,*)DSST
	READ(11,*)GM
C
	CLOSE(UNIT=11,STATUS='SAVE')
C
	IF(HS.LT.0.2.AND.SURFACE.EQ.'swmp')THEN
	 PRINT*, 'Swamp depth must be at least 0.2 meters'
	 PRINT*, 'Setting depth to 0.2 meters and proceeding'
	 HS=MAX(HS,0.2)
	END IF
C
C  ***  Open output files   ***
C
	IF(FMT.EQ.'tab'.OR.FMT.EQ.'bth')THEN
       OPEN(UNIT=12,FILE='output/time.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=14,FILE='output/radius.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=15,FILE='output/vcon.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=16,FILE='output/psicon.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=17,FILE='output/ucon.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=18,FILE='output/wcon.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=19,FILE='output/tcon.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=20,FILE='output/pcon.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=21,FILE='output/tecon.out',
     1  STATUS='UNKNOWN')
	 OPEN(UNIT=22,FILE='output/rgraph.out',
     1  STATUS='UNKNOWN')
	 OPEN(UNIT=23,FILE='output/zgraph.out',
     1  STATUS='UNKNOWN')
	 OPEN(UNIT=24,FILE='output/whov.out',
     1  STATUS='UNKNOWN')
	 OPEN(UNIT=25,FILE='output/vhov.out',
     1  STATUS='UNKNOWN')
	 OPEN(UNIT=26,FILE='output/uhov.out',
     1  STATUS='UNKNOWN')
       OPEN(UNIT=61,FILE='output/raincon.out',
     1  STATUS='UNKNOWN')
	END IF
	IF(FMT.EQ.'csv'.OR.FMT.EQ.'bth')THEN
         OPEN(UNIT=32,FILE='output/time.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=34,FILE='output/radius.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=35,FILE='output/vcon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=36,FILE='output/psicon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=37,FILE='output/ucon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=38,FILE='output/wcon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=39,FILE='output/tcon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=40,FILE='output/pcon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=41,FILE='output/tecon.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=42,FILE='output/whov.csv',
     1     STATUS='UNKNOWN')
         OPEN(UNIT=43,FILE='output/vhov.csv',
     1     STATUS='UNKNOWN')
	   OPEN(UNIT=44,FILE='output/uhov.csv',
     1     STATUS='UNKNOWN')
	   OPEN(UNIT=45,FILE='output/raincon.csv',
     1     STATUS='UNKNOWN')
	END IF
C
C      *** Non-dimensionalize Parameters and Initial Conditions  ***
C
C
C   ***  More external parameters  ***
C
	XMDIF=1500.0*FLOAT(NR)*0.01
	RIC=1.0
	ZGRAD=100.0
	TIMED=TIME
	TIME=TIME+0.1
C
	H_A=0.01*H_A
	AHM=0.01*AHM
	AHMCORE=0.01*AHMCORE
	ES=6.112*EXP(17.67*TS/(243.5+TS))
	EA=H_A*ES
	QS=0.622*ES/PA
	TSA=TS+273.15
	TOA=TO+273.15
	EF=(TS-TO)/TSA
	CHI=2.5E6*EF*QS*(1.-H_A)
	SCHI=SQRT(CHI)
	FC=(3.14159/(12.*3600.))*SIN(3.14159*ALAT/180.)
	CD=CD*0.001
	CD1=CD1*1.0E-5
	ATHA=1000.*LOG(TSA)+2.5E6*QS*H_A/TSA-QS*H_A*
     1    461.*LOG(H_A)
	THETA_E=EXP((ATHA-287.*LOG(0.001*PA))*0.001)
	PT=1000.0*(TOA/THETA_E)**3.484
	DELP=0.5*(PA-PT)
	GAMMA=DPB/DELP
	TM=0.85*TS+0.15*TO
	ESM=6.112*EXP(17.67*TM/(243.5+TM))
	QSM=0.622*ESM/(PA-0.5*DELP)
	TMA=TM+273.15
	GRATB=(1.+2.5E6*QSM/(287.*TMA))/
     1   (1.+2.5E6*2.5E6*QSM/(461.*1005.*TMA*TMA))
	ALENGTH=SCHI/FC
	ATIME=287.*TSA*DELP/(CD*9.8*PA*SCHI)
	BETA=CHI/(287.*TSA)
	Q=0.5*H_A/(GRATB*(1.-H_A))
	TAUC=TAUC*3600./ATIME
	AL=1./TAUC
	TAUR=TAUR*3600./ATIME
	RAD=1./TAUR
	RADMAX=RADMAX/(3600.*24.)
	RADMAX=RADMAX*ATIME*1000.*(TS-TO)/(CHI*GRATB*320.0)
	R0=R0*1000./ALENGTH
	RM=RM*1000./ALENGTH
	RO=RO*1000./ALENGTH	
	ROGD=ROG
	ROG=ROG*1000./ALENGTH
	TIME=TIME*3600.*24./ATIME
	DT=DT/ATIME
	DTG=DTG*3600./ATIME
	TLAND=TLAND*24.*3600./ATIME
	TSHEAR=TSHEAR*24.*3600./ATIME
	nttmax=time/dtg+1
	VM=VM/SCHI
	VEXT1=2.0*VEXT*ATIME/ALENGTH
	XM0=2.5E6*(TS-TO)*QSM*(AHM-1.)/TMA
	XM0=XM0/CHI
	CDV=CD1*SCHI/CD
	RHAA=H_A
	GM=GM*0.01
	UT=UT*ATIME/ALENGTH
	AMIXFAC=0.5*1000.*EF*GM*(1.+2.5E6*2.5E6*QSM/(1000.*461.*
     1   TSA*TSA))/CHI
	HM=HM*AMIXFAC
        DSST=DSST*2.*AMIXFAC/GM
	HEDDY=HEDDY*AMIXFAC
	RWIDE=RWIDE*1000./ALENGTH
	REDDY=EDDYTIME*3600.*24.*UT/ATIME
	TWIDE=TWIDE*3600.*24./ATIME
	AMIX=(2.*RIC*CHI/(9.8*3.3E-4*GM))*1.0E-6*(287.*TSA/9.8)**2*
     1   (DELP/PA)**2*AMIXFAC**4
	AMIX=SQRT(AMIX)
	FACSST=0.002*CD*CECD*CHI*SCHI*ATIME*AMIXFAC/(HS*GM*4160.0)
      DRIG=ROG/FLOAT(MI-2)
C
C   ***   Find Theoretical Maximum Wind and Minimum Pressure
C
	CALL THEORY(CECD,RHAA,TS,TO,PA,VMTH,PCTH,IFAIL)
C
	VMTHI=10.*float(int((vmth+10.05)/10))
	AMTHI=FLOAT(INT(100.*CD*VMTH+1.05))
	if(dim.eq.'n'.or.dim.eq.'N')then
	 vmthi=vmth/schi
	end if
C
C  ***    Initialize certain arrays
C
	do i=1,ntg
	 timeg(i)=float(i-1)*dtg*atime/(3600.*24.)
	 vmg(i,1)=0.0
	 vmg(i,2)=vmth
	 sstt(i)=0.0
	 if(dim.eq.'n'.or.dim.eq.'N')then
	  timeg(i)=float(i-1)*dtg
	  vmg(i,2)=vmth/schi
	 end if
	 do j=1,mi
	  vhov(j,i)=0.0
	  whov(j,i)=0.0
	 end do
	end do
C
C   ***   SET CERTAIN CONSTANTS    ***
C	
	DR=RO/FLOAT(NR-2)
	DRI=1./DR
	TIME=TIME+0.2
	NT=TIME/DT
	DAMP=0.1
	DPR=1./GAMMA
	BFAC=1.-DR*SQRT(2./Q)
	PBM=1.0E-5
	SIXTI=1./16.
C
C              *** SET ACTUAL DIFFUSIVITIES PROPORTIONAL TO DR  ***
C
	PNU=PNU*PNU*DR*0.01744
C
C            ****   INITIALIZE FIELDS   ***
C
	IF(RST.EQ.'y'.OR.RST.EQ.'Y')THEN
C
	OPEN(UNIT=10,FILE='output/restart',STATUS='UNKNOWN')
	DO I=1,NR
	 READ(10,*)RBS1(I),RTS1(I),X1(I),XS1(I),XM1(I),MU1(I),
     1    RBS2(I),RTS2(I),X2(I),XS2(I),XM2(I),MU2(I),PS2(I),
     2	UHMIX1(I),UHMIX2(I),SST1(I),SST2(I),HMIX(I)
	END DO
	CLOSE(10)
C
	ELSE
C
	DO 60 I=2,NR
	 R=FLOAT(I-1)*DR
	 IF(R.GT.R0)GOTO 40
	 IF(R.GT.SQRT(RM*RM*(1.+2.*VM/RM)))GOTO 30
	 RBS1(I)=(R*R/(1.+2.*VM/RM))
	 RBS2(I)=RBS1(I)
	 GOTO 50
   30    RBS1(I)=(0.5*R*R-VM*RM/(1.-(RM/R0)**2))/
     1     (0.5-VM*RM/(R0*R0-RM*RM))
	 RBS2(I)=RBS1(I)
	 GOTO 50
   40    RBS1(I)=R*R
	 RBS2(I)=RBS1(I)
   50    CONTINUE
	 PS2(I)=0.0
	 PS3(I)=0.0
	 P(I)=0.0
	 RTS1(I)=R*R
	 RTS2(I)=RTS1(I)
	 RBS3(I)=RBS2(I)
	 RTS3(I)=RTS2(I)
	 X1(I)=-0.05
	 X2(I)=X1(I)
	 XM1(I)=XM0
	 XM2(I)=XM0
	 MU1(I)=0.0
	 MU2(I)=0.0
	 HMIX(I)=HM+HEDDY/(1.+((SQRT(RBS2(I))-REDDY)/RWIDE)**2)
         UHMIX1(I)=SQRT(HMIX(I)*HMIX(I)*HMIX(I)*DSST)/AMIX
	 UHMIX2(I)=UHMIX1(I)
	 SST1(I)=0.0
	 SST2(I)=SST1(I)
   60   CONTINUE
	XS1(2)=0.0
	HMIX(1)=HM
	DO 70 I=3,NR
	 R=(FLOAT(I-2))*DR
	 XS1(I)=XS1(I-1)+R*DR*0.5*(1.-R*R/RBS1(I-1))
   70 CONTINUE
	PS2(1)=0.0
	MU2(1)=0.0
	SST2(1)=0.0
	SST1(1)=0.0
C
	END IF
C
	DO 75 I=2,NR
	 RB1(I)=SQRT(RBS1(I))
	 RB2(I)=RB1(I)
	 RT1(I)=SQRT(RTS1(I))
	 XS1(I)=XS1(I)-XS1(NR)
	 XS2(I)=XS1(I)
	 IF((RB2(I)*0.001*ALENGTH).LT.RCORE)THEN
	  XM1(I)=XS1(I)+2.5E6*(TS-TO)*QSM*(AHMCORE-1.)/(TMA*CHI)
	  XM2(I)=XM1(I)
       END IF
	 XVIS(I)=0.0
	 XMVIS(I)=0.0
   75 CONTINUE
	PS3(1)=0.0
	PS0(1)=0.0
	VIS(1)=0.0
	GB(1)=0.0
	RMM2(1)=0.0
C
C             ***  SET TIME LOOPING PARAMETERS AND BEGIN TIME LOOP  ***
C
	TT=-DT
	TT1=0.0
	NTT=0
	DT1=DT
C
C           ***  PROGRAM RETURNS TO 77 AFTER EACH TIME STEP  ***
C
   77 CONTINUE
	TT=TT+DT
	IF(TT.GT.TIME)GOTO 705
C
C            *** SET BOUNDARY VALUES FOR AND ADVECTION TERMS ***
C
	RBS1(1)=0.0
	RTS1(1)=0.0
	RBS2(1)=0.0
	RBS3(1)=0.0
	RTS2(1)=0.0
	RTS3(1)=0.0
	RMS2(1)=0.0
	RB1(1)=0.0
	RB2(1)=0.0
	RT1(1)=0.0
	RTS2(NR)=RTS2(NR-1)+DR*(DR+2.*RO)
	RTS1(NR)=RTS1(NR-1)+DR*(DR+2.*RO)
	RBS2(NR)=RBS2(NR-1)+DR*(DR+2.*RO)
	RBS1(NR)=RBS1(NR-1)+DR*(DR+2.*RO)
	RB1(NR)=SQRT(RBS1(NR))
	RB2(NR)=SQRT(RBS2(NR))
	RT1(NR)=SQRT(RTS1(NR))
	X1(1)=X1(2)
	XS2(1)=XS2(2)
	XS1(1)=XS1(2)
	X2(1)=X2(2)
	XM2(1)=XM2(2)
	XM1(1)=XM1(2)
	MU2(1)=MU2(2)
	SST2(1)=SST2(2)
	SST2(NR)=0.0
	MU1(NR-1)=0.0
	MU2(NR-1)=0.0
	X1(NR)=X1(NR-1)
	X2(NR)=X2(NR-1)
	XS1(NR)=XS1(NR-1)
	XS2(NR)=XS2(NR-1)
	XM2(NR)=XM2(NR-1)
	XM1(NR)=XM1(NR-1)
	XVIS(NR)=XVIS(NR-1)
	XMVIS(NR)=XMVIS(NR-1)
	UHMIX1(1)=UHMIX1(2)
	UHMIX2(1)=UHMIX2(2)
	RCENT=REDDY-UT*TT
	TFAC=HM+HEDDY/(1.+((R0-RCENT)/RWIDE)**2)
      UHMIX1(NR)=SQRT(TFAC*TFAC*TFAC*DSST)/AMIX
	UHMIX2(NR)=UHMIX1(NR)
C
	VMAX=0.0
	DO I=2,NR-1
	 R=(FLOAT(I-1))*DR
	 V=0.5*(R*R-RBS2(I))/RB2(I)
	 VMAX=MAX(V,VMAX)
	END DO
	VEXT=0.0
	IF(TT.GE.TSHEAR)VEXT=0.02*VMAX*VMAX*VEXT1*VEXT1
C
C           ***  IN THIS LOOP WE CALCULATE THE INTEGRATED CUMULUS
C                MASS FLUX, THE TOTAL PBL STREAMFUNCTION, AND
C                THE DIFFUSIVE TERMS   ***
C
	DO 108 I=2,NR-1
	 R=(FLOAT(I-1))*DR
	 RP=R+0.5*DR
	 R2=R+DR
	 R1=R-DR
	 R10=R1-DR
c
C            *** CALCULATE INTEGRATED CUMULUS MASS FLUX ***
C
C            ***        Precipitation Efficiency        ***  
C
	 DENOM=MAX(0.001,(X1(I)-XM0))
	 FAC=(XM1(I)-XM0)/DENOM
	 FAC=MAX(FAC,0.0)
	 FAC=MIN(FAC,1.0)
	 EPRE(I)=FAC
C
	 MT=MU2(I)*EPRE(I)
	 GB(I)=GB(I-1)-MT*(RBS2(I)-RBS2(I-1))
C
C            *** CALCULATE EKMAN FLOW PS0(I) ***
C
	 PS0(I)=0.25*((RB1(I+1)-RB1(I-1))
     1     *(R*R-RBS1(I))*ABS(R*R-RBS1(I)))/(R*DR)
	 VREL=0.5*(R*R-RBS1(I))/RB1(I)
	 VABS=ABS(VREL)
	 CDFAC=1.0+CDV*VABS
	 CDFAC=MIN(CDFAC,CDCAP)
C
	IF(TT.GE.TLAND)THEN
	 IF(SURFACE.EQ.'swmp')THEN
	  CDFAC=1.0
	 ELSE IF(SURFACE.EQ.'pln')THEN
	  CDFAC=1.5
	 ELSE IF(SURFACE.EQ.'hill')THEN
	  CDFAC=2.5
	 ELSE IF(SURFACE.EQ.'mtn')THEN
	  CDFAC=4.0
	 END IF
	END IF
C
	 PS0(I)=PS0(I)*CDFAC
C
C            ***   CALCULATE VISCOUS TERM FOR LAYER 1  ***
C
	 C2=(R2*R2/RBS1(I+1)-R*R/RBS1(I))/(RB1(I+1)-RB1(I))
	 RM2=0.5*(RB1(I+1)+RB1(I))
	 RM2=RM2*RM2
	 RM2=RM2*RM2
	 IF(I.EQ.2)THEN
	   C1=(R2*R2/RBS1(I+1)-R*R/RBS1(I))/RB1(I)
	 ELSE
	   C1=(R*R/RBS1(I)-R1*R1/RBS1(I-1))/(RB1(I)-RB1(I-1))
	 END IF
	 RM1=0.5*(RB1(I)+RB1(I-1))
	 RM1=RM1*RM1
	 RM1=RM1*RM1
	 VIS(I)=-(PNU/(R*DR))*(RM2*C2*ABS(C2)-RM1*C1*ABS(C1))
C
C            ***   CALCULATE VISCOUS TERM FOR LAYER 2  ***
C
	 C2=(R2*R2/RTS1(I+1)-R*R/RTS1(I))/(RT1(I+1)-RT1(I))
	 RM2=0.5*(RT1(I+1)+RT1(I))
	 RM2=RM2*RM2
	 RM2=RM2*RM2
	 IF(I.EQ.2)THEN
	   C1=(R2*R2/RTS1(I+1)-R*R/RTS1(I))/RT1(I)
	 ELSE
	   C1=(R*R/RTS1(I)-R1*R1/RTS1(I-1))/(RT1(I)-RT1(I-1))
	 END IF
	 RM1=0.5*(RT1(I)+RT1(I-1))
	 RM1=RM1*RM1
	 RM1=RM1*RM1
	 VIS2(I)=-(PNU/(R*DR))*(RM2*C2*ABS(C2)-RM1*C1*ABS(C1))
C
C         ***  CALCULATE DIFFUSION TERM IN THERMODYNAMIC EQUATION  ***
C
	 IF(I.EQ.2)THEN
	  C2=2.*RBS1(I)*ABS(R2*R2/RBS1(I+1)-R*R/RBS1(I))/RBS1(I+1)
	  BS=C2
	  C1=0.0
	 ELSE
	  C2=RBS1(I)*ABS(R2*R2/RBS1(I+1)-R1*R1/RBS1(I-1))/
     1    ((RB1(I+1)-RB1(I-1))*(RB1(I+1)-RB1(I-1)))
	 END IF
	 IF(I.EQ.3)THEN
	  C1=BS
	 ELSE IF(I.GT.3)THEN
	  C1=RBS1(I-1)*ABS(R*R/RBS1(I)-R10*R10/RBS1(I-2))/
     1    ((RB1(I)-RB1(I-2))*(RB1(I)-RB1(I-2)))
	 END IF
	 XVIS(I)=4.*PNU*(C2*(XS1(I+1)-XS1(I))-C1*(XS1(I)-
     1     XS1(I-1)))/(RBS1(I)-RBS1(I-1))
C
C          *** CALCULATE r**2 AT THE MIDDLE LEVEL ***
C
	 RMS2(I)=RBS1(I)/(0.7+RBS1(I)*0.3/RTS1(I))
	 RMM2(I)=SQRT(RMS2(I))
C
  108   CONTINUE
C
C          ***  Various boundary values   ***
C
	RMS2(NR)=RMS2(NR-1)+DR*(DR+2.*RO)
	RMM2(NR)=SQRT(RMS2(NR))
	GB(NR)=GB(NR-1)
	PS0(NR)=PS0(NR-1)*BFAC
	EPRE(NR)=EPRE(NR-1)
	EPRE(1)=EPRE(2)
C
C   Viscous term for middle layer
C
        DO 109 I=2,NR-1
         R=(FLOAT(I-1))*DR
         RP=R+0.5*DR
         R2=R+DR
         R1=R-DR
         R10=R1-DR
         IF(I.EQ.2)THEN
          C2=2.*RMS2(I)*ABS(R2*R2/RMS2(I+1)-R*R/RMS2(I))/RMS2(I+1)
   	    C2B=RMM2(I)
          BS=C2
	    BSB=C2B
          C1=0.0
	    C1B=0.0
         ELSE
          C2=RMS2(I)*ABS(R2*R2/RMS2(I+1)-R1*R1/RMS2(I-1))/
     1    ((RMM2(I+1)-RMM2(I-1))*(RMM2(I+1)-RMM2(I-1)))
	    C2B=RMM2(I)
         END IF
         IF(I.EQ.3)THEN
          C1=BS
	    C1B=BSB
         ELSE IF(I.GT.3)THEN
          C1=RMS2(I-1)*ABS(R*R/RMS2(I)-R10*R10/RMS2(I-2))/
     1    ((RMM2(I)-RMM2(I-2))*(RMM2(I)-RMM2(I-2)))
	    C1B=RMM2(I-1)
         END IF
	 C1=C1B*4.*XMDIF*PNU
	 C2=C2B*4.*XMDIF*PNU
c	 C1=2.*C1*PNU
c	 C2=2.*C2*PNU
       XMVIS(I)=(C2*(XM1(I+1)-XM1(I))-C1*(XM1(I)-
     1     XM1(I-1)))/(RMS2(I)-RMS2(I-1))
  109	CONTINUE
C
C         *** CALCULATE ARRAYS FOR SOLVING ELLIPTIC EQUATION FOR PSI ***
C
	DO 150 I=2,NR-1
	 R=(FLOAT(I-1))*DR
	 RF2=RAD*XS1(I+1)
	 RF2=MIN(RF2,RADMAX)
	 RF1=RAD*XS1(I)
	 RF1=MIN(RF1,RADMAX)
	 GF1=R*R/RTS2(I)
	 GF2=R*R/RBS2(I)
	 GF=GF1*GF1+GF2*GF2
	 DF=0.5*GF*DR/R
	 Q3(I)=Q/(RMS2(I+1)-RMS2(I))
	 Q2(I)=Q/(RMS2(I)-RMS2(I-1))
	 RAT=(RBS2(I)/RTS2(I))
	 RAT=RAT*RAT
	 PF=(PS0(I)+VIS(I))/(1.+RAT)+GB(I)
	 PF=PF-VIS2(I)*RAT/(1.+RAT)
	 DV(I)=DF*PF+RF2-RF1-XVIS(I+1)+XVIS(I)
	 AF(I)=1./(Q2(I)+Q3(I)+DF)
  150 CONTINUE
C
C           ***  SOLVE ELLIPTIC EQUATION FOR PSI  ***
C
	DO 360 K=1,1000
	 KP=K
	 PMAX=0.0
	 PS3(NR)=PS3(NR-1)*BFAC+GB(NR-1)*(1.-BFAC)
	 DO 320 I=NR-1,2,-1
	  R=(FLOAT(I-1))*DR
	  PS3(I)=AF(I)*(DV(I)+Q3(I)*PS3(I+1)+Q2(I)*PS2(I-1))
	  A=PS3(I)-PS2(I)
	  PMAX=MAX(PMAX,ABS(A))
  320    CONTINUE
	 DO 340 I=2,NR-1
	  PS2(I)=PS3(I)
  340    CONTINUE
	 IF(PMAX.LT.PBM)GOTO 370
  360   CONTINUE
  370   CONTINUE
C
C            ***  SOLVE SURFACE PRESSURE EQUATION  ***
C
	P(2)=0.0
	PAV=0.0
	DO 200 I=3,NR-1
	 R=(FLOAT(I-1))*DR
	 RMI=R-DR
	 R4=R*R*R*R
	 RMI2=RMI*RMI
	 RMI3=RMI*RMI2
	 RMI4=RMI3*RMI
	 RMIM4=(RMI-DR)*(RMI-DR)
	 RMIM4=RMIM4*RMIM4
	 FM1=XS2(2)+SIXTI*(RMI4/RBS2(2)+RBS2(2))
	 IF(I.GT.3)FM1=SIXTI*((RMI4/RBS2(I-1))+RBS2(I-1)+(
     1    RMIM4/RBS2(I-2))+RBS2(I-2))+XS2(I-1)
	 FM2=SIXTI*((R4/RBS2(I))+RBS2(I)+(RMI4/RBS2(I-1))
     1     +RBS2(I-1))+XS2(I)
	 P(I)=P(I-1)+0.5*DR*RMI3/RTS2(I-1)+FM1-FM2
	 PAV=PAV+P(I)*(RBS2(I)-RBS2(I-1))
  200 CONTINUE
	PAV=PAV/RBS2(NR-1)
C
C           ***  RECTIFY PRESSURES  ***
C
	DO 250 I=2,NR-1
	  P(I)=P(I)-PAV
  250 CONTINUE
	P(1)=P(2)
	P(NR)=P(NR-1)
C
C            ***  IN THIS LOOP THE RIGHT-HAND SIDES OF THE
C                 TIME-DEPENDENT EQUATIONS ARE CALCULATED AND THE
C                 QUANTITIES ARE THEN INCREMENTED IN TIME ***
C
	DO 100 I=2,NR-1
	R=(FLOAT(I-1))*DR
	R2=R+DR
	R1=R-DR
C
C            ***  CALCULATE FORCING OF RT AND RB ***
C
	RTF=PS2(I)-GB(I)+VIS2(I)
	RBF=PS0(I)+VIS(I)-PS2(I)+GB(I)
C
C            ***  CALCULATE VERTICAL ADVECTION TERMS IN s* EQN   ***
C
	VADV=-Q2(I)*(PS2(I)-PS2(I-1))
C
C            *** CALCULATE V AT X POINTS AND PBL ADVECTIONS  ***
C
	VB2=0.5*(R*R-RBS2(I))/RB2(I)
c	UB1=PS0(I)/((RB2(I+1)-RB2(I-1)))
	UB1=(PS0(I)+GAMMA*VIS(I))/((RB2(I+1)-RB2(I-1)))
	IF(I.EQ.2)THEN
	 UB1=MAX(UB1,0.0)
	 UB2=0.0
	 VB1=0.0
	ELSE
	 VB1=0.5*((R-DR)*(R-DR)-RBS2(I-1))/RB2(I-1)
c	 UB2=PS0(I-1)/(RB2(I)-RB2(I-2))
	 UB2=(PS0(I-1)+GAMMA*VIS(I-1))/(RB2(I)-RB2(I-2))
	END IF
C
C  ***  Here we form a weighted average of centered and upstream  ***
C  ***      differences to provide some smoothing    ***
C
	UBI=2.*MIN(UB2,0.0)
	UBO=2.*MAX(UB1,0.0)
	UBI=0.8*UB2+0.2*UBI
	UBO=0.8*UB1+0.2*UBO
	ADV=(UBO*(X2(I+1)-X2(I))+UBI*(X2(I)-X2(I-1)))/
     1   (RB2(I)+RB2(I-1))
	ADVM=(UBO*(MU2(I+1)-MU2(I))+UBI*(MU2(I)-MU2(I-1)))/
     1   (RB2(I)+RB2(I-1))
C
C        *** CALCULATE VERTICAL VELOCITY AT TOP OF EKMAN LAYER ***
C
C	W0=(PS0(I)-PS0(I-1))/(RBS2(I)-RBS2(I-1))
	W0=(PS0(I)-PS0(I-1)+GAMMA*(VIS(I)-VIS(I-1)))/(RBS2(I)-RBS2(I-1))
C
C  ***    Calculate downdraft and vertical velocity in between clouds   ***
C
	AMD=-MU2(I)*(1.-EPRE(I))
	W0D=W0-MU2(I)-AMD
C
C  *** Calculate surface fluxes, radiation and equilibrium mass flux   ***
C  ***          Include dissipative heating term here                  ***
C
	VABS=0.5*ABS(VB1+VB2)
	CDFAC=1.0+CDV*VABS
	CDFAC=MIN(CDFAC,CDCAP)
	CKFAC=CDFAC
	IF(TT.GE.TLAND)THEN
	 IF(SURFACE.EQ.'swmp')THEN
	  CKFAC=1.0
	 ELSE
	  CKFAC=0.0
	 END IF
	 CDFAC=1.0
	 IF(SURFACE.EQ.'pln')THEN
	  CDFAC=1.5
	 ELSE IF(SURFACE.EQ.'hill')THEN
	  CDFAC=2.5
	 ELSE IF(SURFACE.EQ.'mtn')THEN
	  CDFAC=4.0
	 END IF
	END IF
C
C   ***  Calculate forcing of ocean temperature  ***
C
	SSTR(I)=0.0
	SEAF=0.0
	SSTF=0.0
      IF(OM.EQ.'Y'.OR.OM.EQ.'y'.AND.TT.LT.TLAND)THEN
	 AO1=(UHMIX2(I+1)-UHMIX2(I))/(RB2(I+1)-RB2(I))
	 AO2=(UHMIX2(I)-UHMIX2(I-1))/(RB2(I)-RB2(I-1))
	 ADVOM=(UT+0.5*RBF/RB2(I))*(0.8*AO1+0.2*AO2)
C	 SEAF=ADVOM+CDFAC*VTEN*VTEN+SPMFAC*VTEN*SPFUNC
	 SEAF=ADVOM+CDFAC*VABS*VABS
	 AHMLOCAL1=HM+HEDDY/(1.+((RB2(I)-RCENT)/RWIDE)**2)
	 AHMLOCAL2=HM+HEDDY/(1.+((RB2(2)-RCENT)/RWIDE)**2)
	 AMIXD1=MAX(HMIX(I),0.0001)
	 AMIXD2=MAX(HMIX(2),0.0001)
	 XSURM1=-((HMIX(I)-AHMLOCAL1)**2)/AMIXD1
	 XSURM2=-((HMIX(2)-AHMLOCAL2)**2)/AMIXD2
       XSURM1=XSURM1+DSST*(AHMLOCAL1-HMIX(I))/AMIXD1
       XSURM2=XSURM2+DSST*(AHMLOCAL2-HMIX(2))/AMIXD2
	 SSTR(I)=0.5*(XSURM1+XSURM2)
	END IF
C
	IF(TT.GE.TLAND.AND.SURFACE.EQ.'swmp')THEN
	 SSTR(I)=SST2(I)
	 XSUR=1.0+SSTR(I)-EF*P(I)+(EXP(-BETA*P(I))-1.)/(1.-H_A)
	 AS1=(SST2(I+1)-SST2(I))/(RB2(I+1)-RB2(I))
	 AS2=(SST2(I)-SST2(I-1))/(RB2(I)-RB2(I-1))
	 ADVOS=(UT+0.5*RBF/RB2(I))*(0.8*AS1+0.2*AS2)
	 SSTF=-FACSST*CKFAC*VABS*(XSUR-X1(I))+ADVOS
	ELSE
	 XSUR=1.0+SSTR(I)-EF*P(I)+(EXP(-BETA*P(I))-1.)/(1.-H_A)
	END IF
C
	FLUX=CKFAC*VABS*(XSUR-X1(I))*CECD
C	FLUX=CKFAC*min(VABS,0.5)*(XSUR-X1(I))*CECD
	FLUX=FLUX+EF*CDFAC*VABS*VABS*VABS
C
	RAD1=RAD
	IF((RAD1*XS1(I)).GT.RADMAX)RAD1=RADMAX/XS1(I)
	IF(W0D.LE.0.0)THEN
	 ACHIM=MAX(0.001,(X1(I)-XM1(I)))
	 MEQ=W0+(ADV+FLUX-GAMMA*(VADV+XVIS(I)-RAD1*XS1(I)))/
     1    ACHIM
	ELSE
	 ATEMP=(X1(I)-XM1(I))*(1.-EPRE(I))
	 ACHIM=MAX(0.001,ATEMP)
	 MEQ=(ADV+FLUX-GAMMA*(VADV+XVIS(I)-RAD1*XS1(I)))/
     1  ACHIM
	END IF
C
C  ***  Limit maximum value of M   ***
C
	MEQ=MIN(MEQ,600.0)
C
C  ***    If conditionally stable, set MEQ to 0    ***
C 
	IF(X2(I).LT.(XS2(I)-0.0001))MEQ=0.0
C
C  ***   Modification introduced 11/6/05 to limit excessive intensity   ***
C
	IF(VMAX.GT.(1.2*VMTH/SCHI))MEQ=0.5*MEQ	
C
C        *** Calculate forcing of mass flux  ***
C
	MF=0.1*ADVM+AL*(MEQ-MU1(I))
C
C        *** CALCULATE FORCING OF ENTROPY EQNS ***
C
	WXM=MAX(0.0,W0D)
c	XMF=(EFRAC*MU2(I)+WXM)*(X1(I)-XM1(I))-(XS1(I)-XM1(I))*
c     1   (PS2(I)-PS2(I-1))/(RMS2(I)-RMS2(I-1))
	XMF=EFRAC*(MU2(I)+WXM)*(X1(I)-XM1(I))
      XMF=XMF+XMVIS(I)-RAD1*XS1(I)*(0.75*GRATB+0.25)
      XMF=XMF-300.*VEXT*(XM1(I)-XM0)
	AMC=MIN(0.0,W0D)
	XF=ADV+FLUX+(AMD+AMC)*(X1(I)-XM1(I))
	XF=XF*DPR
	XSF=VADV+XVIS(I)-RAD1*XS1(I)
C
C           ***  INCREMENT IN TIME   ***
C
	RBS3(I)=RBS1(I)+2.*DT*RBF
	RTS3(I)=RTS1(I)+2.*DT*RTF
	XS3(I)=XS1(I)+2.*DT*XSF
	X3(I)=X1(I)+2.*DT*XF
	X3(I)=MIN(X3(I),XS3(I),XSUR)
	XM3(I)=XM1(I)+2.*DT*XMF
	XM3(I)=MAX(XM3(I),MIN(-1.0,XM0))
	XM3(I)=MIN(XM3(I),XS3(I))
	MU3(I)=MU1(I)+2.*DT*MF
	MU3(I)=MAX(MU3(I),0.0)
	UHMIX3(I)=UHMIX1(I)+2.*DT*SEAF
	UHMIX3(I)=MAX(UHMIX3(I),0.0)
	SST3(I)=SST1(I)+2.*DT*SSTF
C
C  ***    Bail out if r**2 is negative   ***
C
	IF(RBS3(I).LT.0.0)GOTO 710
  100 CONTINUE
C
	SSTR(NR)=SSTR(NR-1)
	SSTR(1)=SSTR(2)
	SST3(1)=SST3(2)
	SST3(NR)=0.0
C
C          ***  STORE CERTAIN GRAPHICS ARRAYS EVERY DTG TIME UNITS  ***
C
	TT1=TT1+DT
	IF(TT1.LT.DTG.AND.TT.GT.0)GOTO 500
	TT1=0
	NTT=NTT+1
	VMAX=0.0
	VTMIN=0.0
	MUMAX=0.0
	MDMIN=0.0
	PMIN=0.0
C
C           ***  CALCULATE EXTREME VALUES OF QUANTITIES    ***
C
	DO 430 I=2,NR-1
	 R=(FLOAT(I-1))*DR
	 V=0.5*(R*R-RBS2(I))/RB2(I)
	 VT=0.5*(R*R-RTS1(I))/RT1(I)
	 IF(V.GE.VMAX)THEN
	  VMAX=V
	  RMAX=RB2(I)
	  PRMAX=P(I)
	 END IF
	 IF(VT.LE.VTMIN)THEN
	  VTMIN=VT
	  RTMIN=RT1(I)
	 END IF
	 PMIN=MIN(PMIN,P(I))
	 MUMAX=MAX(MUMAX,MU3(I))
	 AMD=-MU3(I)*(1.-EPRE(I))
	 MDMIN=MIN(MDMIN,AMD)
  430   CONTINUE
C
C   ***   Create time series    ***
C
	IF(DIM.NE.'y')THEN
	 TIMEG(NTT)=TT
	 VMG(NTT,1)=VMAX
c	 VMG(NTT,2)=VMTH/SCHI
	 RMG(NTT)=RMAX
	 VTMG(NTT)=VTMIN
	 RTMG(NTT)=RTMIN
	 PCG(NTT,1)=PMIN
	 PCG(NTT,2)=LOG(PCTH/PA)/BETA
	 MUG(NTT)=MUMAX
	 MDG(NTT)=MDMIN
	 HMG(NTT)=-HMIX(2)
	 SSTT(NTT)=SSTR(2)
	 IF(TT.GE.TLAND)SSTT(NTT)=SST2(2)
	 DO 435 J=1,MI-1
	  RR=DRIG*FLOAT(J-1)
	  JI=NR
	  DO 434 I=2,NR
	   IF(RB2(I).GT.RR)JI=MIN(JI,I)
  434	  CONTINUE
	  WE1=RR-RB2(JI-1)
	  WE2=RB2(JI)-RR
	  WDEN=1./(RB2(JI)-RB2(JI-1))
	  WHOV(J,NTT)=(MU2(JI)*WE1+MU2(JI-1)*WE2)*WDEN
	  R=DR*FLOAT(JI-1)
	  V1=0.5*(R*R-RBS2(JI))/RB2(JI)
	  U1=(PS2(JI)-GB(JI))/RB2(JI)
	  IF(JI.GT.2)THEN
	   V2=0.5*((R-DR)*(R-DR)-RBS2(JI-1))/RB2(JI-1)
	   U2=(PS2(JI-1)-GB(JI-1))/RB2(JI-1)
	  ELSE
	   V2=0.0
	   U2=0.0
	  END IF
	  VHOV(J,NTT)=(V1*WE1+V2*WE2)*WDEN
	  UHOV(J,NTT)=-(U1*WE1+U2*WE2)*WDEN
  435	 CONTINUE
C
	ELSE
C
	 TIMEG(NTT)=TT*ATIME
	 TIMEG(NTT)=TIMEG(NTT)/(24.*3600.)
	 VMG(NTT,1)=VMAX*SCHI
c	 VMG(NTT,2)=VMTH
	 RMG(NTT)=RMAX*SCHI*0.001/FC
	 VTMG(NTT)=VTMIN*SCHI
	 RTMG(NTT)=RTMIN*SCHI*0.001/FC
c	 PCG(NTT,1)=PA*EXP(BETA*PMIN)
	 PEXP=0.5
	 PCG(NTT,1)=PA*EXP(BETA*PRMAX-VMG(NTT,1)**2/
     1	 (2.*PEXP*287.*TSA))
	 PCG(NTT,2)=PCTH
	 MUG(NTT)=MUMAX*CD*SCHI
	 MDG(NTT)=MDMIN*CD*SCHI
	 HMG(NTT)=-HMIX(2)/AMIXFAC
	 SSTT(NTT)=SSTR(2)*0.5*GM/AMIXFAC
	 IF(TT.GE.TLAND)SSTT(NTT)=SST2(2)*0.5*GM/AMIXFAC
         DO 445 J=1,MI-1
          RR=DRIG*FLOAT(J-1)
          JI=NR
          DO 444 I=2,NR
           IF(RB2(I).GT.RR)JI=MIN(JI,I)
  444     CONTINUE
          WE1=RR-RB2(JI-1)
          WE2=RB2(JI)-RR
          WDEN=1./(RB2(JI)-RB2(JI-1))
          WHOV(J,NTT)=(MU2(JI)*WE1+MU2(JI-1)*WE2)*WDEN
          R=DR*FLOAT(JI-1)
          V1=0.5*(R*R-RBS2(JI))/RB2(JI)
          U1=(PS2(JI)-GB(JI))/RB2(JI)
          IF(JI.GT.2)THEN
           V2=0.5*((R-DR)*(R-DR)-RBS2(JI-1))/RB2(JI-1)
           U2=(PS2(JI-1)-GB(JI-1))/RB2(JI-1)
          ELSE
           V2=0.0
           U2=0.0
          END IF
          VHOV(J,NTT)=(V1*WE1+V2*WE2)*WDEN
          UHOV(J,NTT)=-(U1*WE1+U2*WE2)*WDEN
	  WHOV(J,NTT)=WHOV(J,NTT)*CD*SCHI
	  VHOV(J,NTT)=VHOV(J,NTT)*SCHI
	  UHOV(J,NTT)=UHOV(J,NTT)*1.25*9.8*CD*BETA/FC
  445    CONTINUE
	END IF
C
C         *** IF VMAX EXCEEDS CERTAIN VALUES, CUT DOWN TIME STEP ***
C
	IF(VMAX.GT.1.5)DT=0.5*DT1
	IF(VMAX.GT.2.0)DT=0.2*DT1
C
  500 CONTINUE
C
C            ***  ADVANCE VARIABLES WITH TIME SMOOTHER ***
C
	DO 600 I=2,NR-1
	 RBS1(I)=RBS2(I)+DAMP*(RBS1(I)+RBS3(I)-2.*RBS2(I))
	 RTS1(I)=RTS2(I)+DAMP*(RTS1(I)+RTS3(I)-2.*RTS2(I))
	 RB1(I)=SQRT(RBS1(I))
	 RT1(I)=SQRT(RTS1(I))
	 X1(I)=X2(I)+DAMP*(X1(I)+X3(I)-2.*X2(I))
	 XS1(I)=XS2(I)+DAMP*(XS1(I)+XS3(I)-2.*XS2(I))
	 XM1(I)=XM2(I)+DAMP*(XM1(I)+XM3(I)-2.*XM2(I))
	 MU1(I)=MU2(I)+DAMP*(MU1(I)+MU3(I)-2.*MU2(I))
	 MU1(I)=MAX(MU1(I),0.0)
	 UHMIX1(I)=UHMIX2(I)+DAMP*(UHMIX1(I)+UHMIX3(I)-2.*UHMIX2(I))
	 SST1(I)=SST2(I)+DAMP*(SST1(I)+SST3(I)-2.*SST2(I))
C
	 RBS2(I)=RBS3(I)
	 RTS2(I)=RTS3(I)
	 RB2(I)=SQRT(RBS2(I))
	 X2(I)=X3(I)
	 XS2(I)=XS3(I)
	 XM2(I)=XM3(I)
	 MU2(I)=MU3(I)
	 SST2(I)=SST3(I)
	 UHMIX2(I)=UHMIX3(I)
         IF(OM.EQ.'Y'.OR.OM.EQ.'y')THEN
	  AHM0=HM+HEDDY/(1.+((RB2(I)-RCENT)/RWIDE)**2)
          AHM02=AHM0*(AHM0-DSST)
	  HMIX(I)=SQRT(0.5*AHM02+SQRT(0.25*AHM02*AHM02+AMIX*AMIX*
     1     UHMIX2(I)*UHMIX2(I)))
	 END IF
  600   CONTINUE
C
  700   GOTO 77
  705   CONTINUE
C
	OPEN(UNIT=10,FILE='output/restart',STATUS='UNKNOWN')
	DO I=1,NR
	 WRITE(10,*)RBS1(I),RTS1(I),X1(I),XS1(I),XM1(I),MU1(I),
     1    RBS2(I),RTS2(I),X2(I),XS2(I),XM2(I),MU2(I),PS2(I),
     2    UHMIX1(I),UHMIX2(I),SST1(I),SST2(I),HMIX(I)
	END DO
	CLOSE(10)
C
C         ***  Fill Graphics Arrays     ***
C
	NG=NR
	DO 490 I=2,NR
	 R=(FLOAT(I-1))*DR
	 RG(I)=RB2(I)
	 RGP(I)=0.5*(RB2(I)+RB2(I-1))
	 RTG(I)=RT1(I)
	 IF(RG(I).LE.ROG)THEN
	  NG=I
	 END IF
	 VBG(I)=0.5*(R*R-RBS2(I))/RG(I)
	 VTG(I)=0.5*(R*R-RTS1(I))/RTG(I)
	 PG(I)=P(I)
	 IF(TT.GE.TLAND)SSTR(I)=SST2(I)
       XSUR=1.0+SSTR(I)-EF*P(I)+(EXP(-BETA*P(I))-1.)/(1.-H_A)
       HUMG(I)=1.-EXP(BETA*P(I))*(1.-H_A)*(XSUR-X2(I))
	 IF(I.GT.2) THEN
	  WG(I,1)=10.*(PS2(I)-PS2(I-1))/(RBS2(I)-RBS2(I-1))
	 ELSE
	  WG(I,1)=10.*PS2(2)/RBS2(2)
	 END IF
	 WG(I,2)=MU2(I)
	 UG(I,1)=(PS2(I)-GB(I))/RG(I)
	 UG(I,2)=PS0(I)/RG(I)
	 XG(I,1)=X2(I)
	 XG(I,2)=XS2(I)
	 XG(I,3)=XM2(I)
C
C   *** Re-dimensionalize Output, If Desired ***
C
	 HMIX(I)=-HMIX(I)
	 IF(DIM.EQ.'y')THEN
	  RG(I)=RG(I)*0.001*ALENGTH
	  RGP(I)=RGP(I)*0.001*ALENGTH
	  RTG(I)=RTG(I)*0.001*ALENGTH
	  VBG(I)=VBG(I)*SCHI
	  VTG(I)=VTG(I)*SCHI
	  HMIX(I)=HMIX(I)/AMIXFAC
	  SSTR(I)=SSTR(I)*0.5*GM/AMIXFAC
	  PG(I)=PA*EXP(BETA*P(I))
	  DO 488 J=1,2
	   WG(I,J)=WG(I,J)*CD*SCHI
	   UG(I,J)=UG(I,J)*1.25*9.8*CD*BETA/FC
  488     CONTINUE
	  DO 489 J=1,3
	   XG(I,J)=EXP((CHI*XG(I,J)/(TS-TO)+ATHA)*0.001)
  489     CONTINUE
	 END IF
  490   CONTINUE
	RG(1)=0.0
	RGP(1)=0.0
	VBG(1)=0.0
	RTG(1)=0.0
	VTG(1)=0.0
	UG(1,1)=0.0
	UG(1,2)=0.0
	WG(1,1)=WG(1,2)
	WG(1,2)=WG(2,2)
	XG(1,1)=XG(2,1)
	XG(1,2)=XG(2,2)
	XG(1,3)=XG(2,3)
      HUMG(1)=HUMG(2)
	HMIX(1)=HMIX(2)
	SSTR(1)=SSTR(2)
	PG(1)=PG(2)
	NRM=NR-1
	NGM=MIN(NG,NRM)
C
C   ***    Calculate Two-dimensional Distributions   ***
C
C   ***   First calculate 2-D arrays of r,V,p,z,p',T',and psi  ***
C
	DO 495 J=1,NZD-1
	 ZTEM=FLOAT(J-1)/FLOAT(NZD-2)
C	 ZTE(J)=0.5+1.5*(ZTEM-0.5)-2.*(ZTEM-0.5)**3
	 ZTE(J)=ZTEM
	 ZTE(J)=MAX(ZTE(J),0.0)
  495	CONTINUE
	ZTE(NZD)=2.-ZTE(NZD-2)
	DO 550 I=NRM,2,-1
	 R=FLOAT(I-1)*DR
C
C   ***   Pressure and azimuthal velocity at sea surface  ***
C
	 PR=PA*EXP(BETA*P(I))
	 PTEMP=P(I)
	 VB=0.5*(R*R-RBS2(I))/RB2(I)
C
C   ***   Physical radius and azimuthal velocity   ***
C
	 DO 530 J=1,NZD-1
	  RRT(I,J)=(RBS2(I)*RTS2(I)/(RBS2(I)*ZTE(J)+
     1     RTS2(I)*(1.-ZTE(J))))
          RRTI(I,J)=SQRT(MAX(RRT(I,J),0.0))
	  VVT(I,J)=0.5*(R*R-RRT(I,J))/RRTI(I,J) 
	  RMT(I,J)=R
C
C   ***    Iteration using Newton's method to find pressure     ***
C   ***               and q_s by inverting theta_es             ***
C
	  FAC=1.-EF*ZTE(J)
	  TTEMP=TSA*FAC
	  TTC=TTEMP-273.15
	  EST=6.112*EXP(17.67*TTC/(243.5+TTC))
	  AFAC=3.484*LOG(FAC)/BETA-XS2(I)/EF
	   DO 528 JK=1,15
	    QST=0.622*EST/PR
	    RES=AFAC+(QST/(QS*FAC)-H_A)/(EF*(1.-H_A))-PTEMP
	    DRDP=-QST/(QS*FAC*EF*(1.-H_A)*PR)-1./
     1        (BETA*PR)
	    PR=PR-RES/DRDP
	    PTEMP=LOG(PR/PA)/BETA
  528	   CONTINUE
	  IF(J.EQ.1)QSR=QST
	  PRT(I,J)=PR
	  TR=TSA-ZTE(J)*(TS-TO)
C
C   ***     Physical height by conservation of total energy   ***
C
	  ZZT(I,J)=0.5*(VB*VB-VVT(I,J)*VVT(I,J))+3.484*ZTE(J)*
     1     EF/BETA+(QSR-QST)/(EF*QS*(1.-H_A))
C
C   ***     Outer boundary values of pressure and altitude    ***
C
	  IF(I.EQ.NRM)THEN
	   PBO(J)=PTEMP
	   ZBO(J)=ZZT(I,J)
	  END IF
C
C   ***    Perturbation temperature and pressure              ***
C
        GMGD=(1.+8710.*QST/TR)/(1.+1.36E7*QST/(TR*TR))
	  TTE(I,J)=0.287*TSA*BETA*(ZZT(I,J)-ZBO(J))*GMGD
	  TMINT=(TS-TO)*(1.-ZTE(J))
	  TTE(I,J)=MIN(TTE(I,J),TMINT)
	  PPT(I,J)=PTEMP-PBO(J)+3.484*TTE(I,J)/(TSA*BETA*GMGD*FAC)
	  PCMAX=P(I)-P(NRM)
	  PPT(I,J)=MAX(PPT(I,J),PCMAX)
	  WORKIT(I,J)=(10.*CD*SCHI*CHI/(287.*TR*7.0))*QST*(-1.+
     1   GMGD*1555./TR)
	  IF(DIM.EQ.'y')THEN
	   PPT(I,J)=PR*(1.-EXP(-BETA*PPT(I,J)))
	  ELSE
	   PPT(I,J)=PPT(I,J)*EXP(BETA*PBO(J))   
	  END IF
  530	 CONTINUE
  550 CONTINUE
C
C   ***   Mass streamfunction      ***
C
	Z0=0.5*GAMMA
	DO 540 I=NRM,2,-1
	PSIT(I,1)=0.0
	 DO 538 J=2,NZD-1
	  DZ=ZTE(J+1)-ZTE(J)
	  RAT=RBS2(I)/RTS2(I)
	  RTF=PS2(I)-GB(I)+VIS2(I)
	  RBF=PS0(I)+VIS(I)-PS2(I)+GB(I)
	  IF(ZTE(J).LE.(Z0+0.5*DZ))THEN
	   PSIT(I,J)=(PS0(I)+VIS(I))*ZTE(J)*(2.*Z0-ZTE(J))/(Z0*Z0)
	  ELSE
	   UCF=(RRT(I,J)*RRT(I,J)/(RBS2(I)*RBS2(I)))*(RBF*(
     1      1.-ZTE(J))+ZTE(J)*RAT*RAT*RTF)/DELP
	   PSIT(I,J)=PSIT(I,J-1)+(PRT(I,J)-PRT(I,J-1))*UCF
	  END IF
  538	 CONTINUE
  540	CONTINUE
C
C   ***  Add a layer of points at top to make graphs nicer   ***
C
	DO 541 I=2,NRM
	 ZZT(I,NZD)=ZZT(2,NZD-1)+3.0
	 PRT(I,NZD)=PRT(2,NZD-1)*EXP(-3.*BETA/(1.-EF))
	 RRT(I,NZD)=RRT(I,NZD-1)
	 RMT(I,NZD)=DR*FLOAT(I-1)
	 RRTI(I,NZD)=SQRT(MAX(RRT(I,NZD),0.0))
	 VVT(I,NZD)=0.0
	 TTE(I,NZD)=0.0
	 PPT(I,NZD)=PPT(I,NZD-1)
	 PSIT(I,NZD)=0.0
  541	CONTINUE	
C
C   ***   Boundary conditions at R=0    ***
C
	DO 542 J=1,NZD
	 ZZT(1,J)=ZZT(2,J)
	 RRT(1,J)=0.0
	 RMT(1,J)=0.0
	 RRTI(1,J)=0.0
	 RRT(NR,J)=(SQRT(RRT(NRM,J))+DR)**2
	 VVT(1,J)=0.0
	 TTE(1,J)=TTE(2,J)
	 PPT(1,J)=PPT(2,J)
	 PRT(1,J)=PRT(2,J)
	 PSIT(NR,J)=PSIT(NRM,J)
	 PSIT(1,J)=0.0
  542 CONTINUE
C
C   ***  Fill graphics arrays and calculate theta_e, u and w   ***
C
	ZM=0.2
	PI=3.14159
	DO 546 I=NRM,2,-1
	 DO 544 J=1,NZD
	  RN=SQRT(RRT(I,J))
C
C   ***  Fit a sine curve to X2 and XM to create 2-D theta_e plots  ***
C   ***     (ZM is the level of minimum theta_e)                    ***
C
	   AN=LOG(2.0)/LOG(1.0/ZM)
	   XIT(I,J)=X2(I)-(X2(I)-XM2(I))*SIN(0.5*PI*(ZTE(J)/ZM)**AN)
	   IF(J.EQ.NZD)THEN
c            XIT(I,J)=XS2(I)+EF*(ZZT(I,J)-ZZT(I,J-1))/(1.-EF)
	    XIT(I,J)=2.*XIT(I,J-1)-XIT(I,J-2)
	   END IF
C
C   ***     Calculate radial and vertical velocities       ***
C   ***  Mesh boundary layer with interior radial velocity ***
C
	   IF(J.GT.1)THEN
	    DENOM=(PRT(I,J-1)-PRT(I,J))
	    DENOM2=(RRT(I,J)-RRT(I-1,J))
	    DENOM3=1.-(RRT(I,J-1)-RRT(I,J))*(PRT(I,J)-
     1        PRT(I-1,J))/(DENOM*DENOM2)
	    DENOM3=MAX(DENOM3,1.0E-12)
	    UIT(I,J)=DELP*(PSIT(I,J-1)-PSIT(I,J))/
     1       DENOM-(PSIT(I,J)-PSIT(I-1,J))*
     2       DELP*(RRT(I,J-1)-RRT(I,J))/(DENOM*DENOM2)
	    UIT(I,J)=UIT(I,J)/(RN*DENOM3)
	    WIT(I,J)=(PSIT(I,J)-PSIT(I-1,J))/DENOM2-
     1       RN*UIT(I,J)*(PRT(I,J)-PRT(I-1,J))/(DELP*DENOM2)
	    WIT(I,J)=WIT(I,J)*(1.-EF*ZTE(J))*PA/PRT(I,J)
	   ELSE
	    UIT(I,J)=0.0
	    WIT(I,J)=0.0
	   END IF
	   IF(WIT(I,J).GT.0.0)THEN
	    WORKIT(I,J)=WORKIT(I,J)*WIT(I,J)/(1.-
     1	   MIN((0.2*CD*SCHI*WIT(I,J)),0.95))
	   ELSE
	    WORKIT(I,J)=-0.5*CHI*(X2(I)-XIT(I,J))/(9.8*1000.*(X2(I)-XM0))
	   END IF
  544	 CONTINUE
  546	CONTINUE	
C
C   ***   Boundary values for graphics arrays at r=0   ***
C
	DO 552 J=1,NZD
	 WIT(1,J)=WIT(2,J)
	 XIT(1,J)=XIT(2,J)
	 UIT(1,J)=0.0
	 WORKIT(1,J)=WORKIT(2,J)
  552	CONTINUE
c
c   ***   Linearly interpolate fields to regular grid   ***
c
	ZMAX=ZZT(3,NZD)-0.01
C
C   ***   DRI and DZI are the nondimensional radial and vertical   ***
C   ***            increments for the 2-d plots                    ***
C
	DRI=ROG/FLOAT(MI-2)
	DZI=ZMAX/FLOAT(NI-2)
C
	CALL INTERGRAPH(RMT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     RMTI,NRI,NZI)
      CALL INTERGRAPH(PSIT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     PSITI,NRI,NZI)
      CALL INTERGRAPH(TTE,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     TTEI,NRI,NZI)
      CALL INTERGRAPH(UIT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     UITI,NRI,NZI)
      CALL INTERGRAPH(PPT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     PPTI,NRI,NZI)
      CALL INTERGRAPH(WIT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     WITI,NRI,NZI)
      CALL INTERGRAPH(XIT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     XITI,NRI,NZI)
      CALL INTERGRAPH(WORKIT,RRTI,ZZT,NRD,NZD,ROG,ZMAX,DRI,DZI,
     1     WORKITI,NRI,NZI)
C
C  *** Calculate gridded azimuthal velocity from the gridded  ***
C  ***                 potential radius                       ***
C
	DO 558 J=1,NZI
	 VVTI(1,J)=0.0
	 DO 558 I=2,NRI
	  RTEMP=DRI*FLOAT(I-1)
	  VVTI(I,J)=0.5*(RMTI(I,J)*RMTI(I,J)-RTEMP*RTEMP)/RTEMP
	  PSITI(I,J)=PSITI(I,J)*1000.0
	  IF(WITI(I,J).LT.0.0)THEN
	   WITI(I,J)=100.*WITI(I,J)
	  END IF
  558	CONTINUE
C
C   ***  Calculate rainwater content from vertical advection equation  ***
C
	DO I=1,NRI
	 RAIN(I,NZI)=0.0
	 DO J=NZI-1,1,-1
	  IF(WORKITI(I,J).GT.0.0)THEN
	   RAIN(I,J)=RAIN(I,J+1)+WORKITI(I,J)*DZI
	  ELSE
	   RAIN(I,J)=RAIN(I,J+1)*(1.+DZI*WORKITI(I,J))/
     1    (1.-DZI*WORKITI(I,J))
	   RAIN(I,J)=MAX(RAIN(I,J),0.0)
	  END IF
	 END DO
	END DO
C
C   ***     Write to output files     ***
C
	IF(FMT.EQ.'csv'.OR.FMT.EQ.'bth')THEN
C
        WRITE(32,609)
  609   FORMAT(2X,'Time,','V_b, ','V_theory,','r_b, ','V_t ,',
     1   'r_t, ','p_c, ','p_theory,','M_u, ','M_d,','h_mix(1),',
     2   'SST')
        DO 601 I=1,NTT
         WRITE(32,611)TIMEG(I),VMG(I,1),VMG(I,2),RMG(I),VTMG(I),
     1    RTMG(I),PCG(I,1),PCG(I,2),MUG(I),MDG(I),HMG(I),SSTT(I)
  601   CONTINUE
  611   FORMAT(2X,11(F12.6,','),F12.6)
        WRITE(34,625)
  625   FORMAT(2X,'r_b, ','r_bi, ','r_t, ','v_b, ','v_t, ','p, ',
     1  'u_0, ','u, ','M_u, ','w_e, ','theta_e, ','theta_es, ',
     2  'theta_em,','h_mix,','SST')
        DO 620 I=1,NG
         WRITE(34,630)RG(I),RGP(I),RTG(I),VBG(I),VTG(I),PG(I),UG(I,1),
     1   UG(I,2),WG(I,1),WG(I,2),XG(I,1),XG(I,2),XG(I,3),HMIX(I),
     2   SSTR(I),HUMG(I)
  620   CONTINUE
  630   FORMAT(2X,15(F12.6,','),F12.6)
C
	DO 635 I=1,NRI
	 RRI(I)=DRI*FLOAT(I-1)
	 IF(DIM.EQ.'y')RRI(I)=RRI(I)*0.001*ALENGTH
  635	CONTINUE
	  WRITE(35,636)(RRI(I),I=1,NRI)
        WRITE(36,636)(RRI(I),I=1,NRI)
        WRITE(37,636)(RRI(I),I=1,NRI)
        WRITE(38,636)(RRI(I),I=1,NRI)
        WRITE(39,636)(RRI(I),I=1,NRI)
        WRITE(40,636)(RRI(I),I=1,NRI)
        WRITE(41,636)(RRI(I),I=1,NRI)
        WRITE(42,636)(RRI(I),I=1,NRI)
        WRITE(43,636)(RRI(I),I=1,NRI)
        WRITE(44,636)(RRI(I),I=1,NRI)
        WRITE(45,636)(RRI(I),I=1,NRI)
  636	FORMAT(2X,F13.8,200(',',F13.8))
	 DO 640 J=NZI,1,-1
	  ZZI(J)=DZI*FLOAT(J-1)
	  IF(DIM.EQ.'y')THEN
	   ZZI(J)=ZZI(J)*0.001*CHI/9.8
	   DO 637 I=1,NRI
	   VVTI(I,J)=VVTI(I,J)*SCHI
	   PSITI(I,J)=PSITI(I,J)*0.59*9.8*CD*SCHI*ALENGTH*ALENGTH
	   PSITI(I,J)=PSITI(I,J)*1.0E-9 
	   WITI(I,J)=WITI(I,J)*CD*SCHI
	   UITI(I,J)=UITI(I,J)*PA*9.8*CD*BETA/(FC*2.*DELP)
	   XITI(I,J)=EXP((CHI*XITI(I,J)/(TS-TO)+ATHA)*0.001)
  637	   CONTINUE
	  END IF
          WRITE(35,645)ZZI(J),(VVTI(I,J),I=1,NRI)
          WRITE(36,645)ZZI(J),(PSITI(I,J),I=1,NRI)
          WRITE(37,645)ZZI(J),(UITI(I,J),I=1,NRI)
          WRITE(38,645)ZZI(J),(WITI(I,J),I=1,NRI)
          WRITE(39,645)ZZI(J),(TTEI(I,J),I=1,NRI)
          WRITE(40,645)ZZI(J),(PPTI(I,J),I=1,NRI)
          WRITE(41,645)ZZI(J),(XITI(I,J),I=1,NRI)
	    WRITE(45,645)ZZI(J),(RAIN(I,J),I=1,NRI)
  640	CONTINUE
  645   FORMAT(2X,F12.7,',',200(F10.4,','))
C
	DO 648 J=1,NTT
	 WRITE(42,645)TIMEG(J),(WHOV(I,J),I=1,NRI)
	 WRITE(43,645)TIMEG(J),(VHOV(I,J),I=1,NRI)
	 WRITE(44,645)TIMEG(J),(UHOV(I,J),I=1,NRI)
  648	CONTINUE
C
	END IF
	IF(FMT.EQ.'tab'.OR.FMT.EQ.'bth')THEN
C
        DO 721 I=1,NTT
         WRITE(12,722)TIMEG(I),VMG(I,1),VMG(I,2),RMG(I),VTMG(I),
     1    RTMG(I),PCG(I,1),PCG(I,2),MUG(I),MDG(I),HMG(I),SSTT(I)
  721   CONTINUE
  722   FORMAT(2X,11(F12.6,' '),F12.6)
        DO 730 I=1,NG
         WRITE(14,740)RG(I),RGP(I),RTG(I),VBG(I),VTG(I),PG(I),UG(I,1),
     1   UG(I,2),WG(I,1),WG(I,2),XG(I,1),XG(I,2),XG(I,3),HMIX(I),
     2   SSTR(I),HUMG(I)
  730   CONTINUE
  740   FORMAT(2X,15(F12.6,' '),F12.6)
	DO 750 I=1,NRI
	 RRI(I)=DRI*FLOAT(I-1)
	 IF(DIM.EQ.'y')RRI(I)=RRI(I)*0.001*ALENGTH
  750	CONTINUE
	WRITE(22,755)(RRI(I),I=1,NRI)
  755	FORMAT(2X,800(F8.3,' '))
	 DO 760 J=1,NZI     
	  ZZI(J)=DZI*FLOAT(J-1)
	  IF(DIM.EQ.'y')THEN
	   ZZI(J)=ZZI(J)*0.001*CHI/9.8
	  END IF
  760	 CONTINUE
	 WRITE(23,765)(ZZI(J),J=1,NZI)
  765	 FORMAT(2X,200(F8.4,1X))
	 DO 780 J=1,NZI
          IF(DIM.EQ.'y')THEN
	   DO 770 I=1,NRI
	   VVTI(I,J)=VVTI(I,J)*SCHI
	   PSITI(I,J)=PSITI(I,J)*0.59*9.8*CD*SCHI*ALENGTH*ALENGTH
	   PSITI(I,J)=PSITI(I,J)*1.0E-9 
	   WITI(I,J)=WITI(I,J)*CD*SCHI
	   UITI(I,J)=UITI(I,J)*PA*9.8*CD*BETA/(FC*2.*DELP)
	   XITI(I,J)=EXP((CHI*XITI(I,J)/(TS-TO)+ATHA)*0.001)
  770	   CONTINUE
	  END IF
          WRITE(15,790)(VVTI(I,J),I=1,NRI)
          WRITE(16,790)(PSITI(I,J),I=1,NRI)
          WRITE(17,790)(UITI(I,J),I=1,NRI)
          WRITE(18,790)(WITI(I,J),I=1,NRI)
          WRITE(19,790)(TTEI(I,J),I=1,NRI)
          WRITE(20,790)(PPTI(I,J),I=1,NRI)
          WRITE(21,790)(XITI(I,J),I=1,NRI)
	    WRITE(61,790)(RAIN(I,J),I=1,NRI)
  780	 CONTINUE
  790   FORMAT(2X,F10.4,' ',800(F10.4,' '))
C
	DO 850 I=1,NTT
	 WRITE(24,790)(WHOV(J,I),J=1,MI-1)
	 WRITE(25,790)(VHOV(J,I),J=1,MI-1)
	 WRITE(26,790)(UHOV(J,I),J=1,MI-1)
  850	CONTINUE
C
	END IF
  710 CONTINUE
  800 CONTINUE
  880	continue
	STOP
	END
C
C--------------------------------------------------------------------
C
       SUBROUTINE THEORY(CKCD,H,TS,TO,PA,VM,PC,IFAIL)
C
C       This subroutine calculates the theoretical maximum wind
C         speed and minimum central pressure.
C
       REAL H, LV, PA, PC, PM
       DELTAT=0.0    
       LV=2.5E6
       RD=287.0
       RV=461.0 
       IFAIL=0
C
C      AN is the assumed power dependence of v on radius inside the
C        radius of maximum winds (i.e. v~r**an) used to calculate PC
C
       AN=1.0
C       
        ES=6.112*EXP(17.67*TS/(243.5+TS))
        EP=(TS-TO)/(TO+273.15)
        COEF1=EP*LV*ES/(RV*(TS+273.15))
        COEF2=0.5*CKCD*(1.-H)*(EP+1.)
        COEF3=0.5*CKCD*EP*1000.*DELTAT/(RD*(TS+273.15)*(1.-EP*H))
        PM=PA   
        N=0
   10   CONTINUE
        N=N+1
        PG=PA*EXP(-COEF1*(COEF2/PM+H*(1./PM-1./PA))-COEF3)
        IF(ABS(PG-PM).LT.0.1)THEN
         PM=0.5*(PM+PG)
         VM=SQRT(EP*CKCD*(LV*0.622*ES*(1.-H)/PM+1000.*DELTAT))
         GOTO 20
        END IF 
        IF(N.GT.1000.OR.PG.LE.1.0)THEN
          IFAIL=1
          GOTO 20
        END IF
        PM=0.5*(PM+PG)  
        GOTO 10
   20   CONTINUE
	IF(IFAIL.EQ.1)THEN
	 PC=PA
	 VM=0.0
	ELSE
         PC=PM*EXP(-VM*VM/(2.*AN*RD*(TS+273.15)))
	END IF
C
      RETURN
      END         
C      
	SUBROUTINE INTERGRAPH(A,R,Z,M,N,RMAX,ZMAX,DR,DZ,B,NR,NZ)
	PARAMETER(MI=100,NI=150)
	DIMENSION A(M,N),Z(M,N),R(M,N),B(MI,NI),C(M,NI),RC(M,NI)
	NR=RMAX/DR+1.1
	NZ=ZMAX/DZ+1.1
	DO 100 I=1,M 
	 C(I,1)=A(I,1)
	 RC(I,1)=R(I,1)
	 DO 50 J=2,NZ
	  ZI=DZ*FLOAT(J-1)
	  KI=2
	  DO 30 K=2,N
	   IF(Z(I,K).GT.ZI)THEN
	    KI=K
	    GOTO 40
	   END IF
   30	  CONTINUE
   40	  CONTINUE
	  DELZ1=Z(I,KI)-ZI
          DELZ1=MAX(DELZ1,0.00001)
	  DELZ2=ZI-Z(I,KI-1)
          DELZ2=MAX(DELZ2,0.00001)
	  C(I,J)=(A(I,KI)*DELZ2+A(I,KI-1)*DELZ1)/(DELZ1+DELZ2)
	  RC(I,J)=(R(I,KI)*DELZ2+R(I,KI-1)*DELZ1)/(DELZ1+DELZ2)
   50	 CONTINUE
  100	CONTINUE
	DO 200 J=1,NZ
	 B(1,J)=C(1,J)
	 DO 150 I=2,NR
	  RI=DR*FLOAT(I-1)
	  LI=2
	  DO 130 K=2,M
	   IF(RC(K,J).GT.RI)THEN
	    LI=K
	    GOTO 140
	   END IF
  130	  CONTINUE
  140	  CONTINUE
	  DELR1=RC(LI,J)-RI
        DELR1=MAX(DELR1,0.00001)
	  DELR2=RI-RC(LI-1,J)
	  DELR2=MAX(DELR2,0.00001)
	  B(I,J)=(C(LI,J)*DELR2+C(LI-1,J)*DELR1)/(DELR1+DELR2)
  150	 CONTINUE
  200	CONTINUE
	DO 300 I=1,NR
	 B(I,1)=2.*B(I,2)-B(I,3)
  300	CONTINUE
	RETURN
	END
