	PROGRAM WRT_AER
C
C  ***       This program uses ECMWF ensmebles to generate random           ***
C  ***       tracks from monthly 850 and 250 mb mean flow and variance      ***
C  ***       It then runs CHIPS to produce an intensity estimate            ***
C  ***       It can run off of models or NCEP reanalysis data, and can      ***
C  ***       alternatively draw from best track genesis pdfs.               ***
C
C  ***               Created 2/16/2010 by K. Emanuel                        ***
C  ***                Copyright WindRiskTech, 2010                          ***
C
C  ***     Modified June, 2010, to allow for TC motion to be influenced     ***
C  ***              by ECMWF track ensemble-derived motion.                 ***
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=100, NTG=600)
c
C  ***      Flow-related files  ***
C
	REAL*8, ALLOCATABLE :: U850(:,:,:), U250(:,:,:),V850(:,:,:)
	REAL*8, ALLOCATABLE :: V250(:,:,:)
	REAL*8, ALLOCATABLE :: UVAR250(:,:,:), VVAR250(:,:,:)
	REAL*8, ALLOCATABLE :: UVAR850(:,:,:),VVAR850(:,:,:)
	REAL*8, ALLOCATABLE :: U250V250(:,:,:),U250U850(:,:,:)
	REAL*8, ALLOCATABLE :: V250V850(:,:,:), V850U850(:,:,:)
	REAL*8, ALLOCATABLE :: UEC(:),VEC(:),U2EC(:),V2EC(:),UVEC(:)
        REAL, ALLOCATABLE :: VOR(:,:,:)
	REAL*8, ALLOCATABLE :: VTEM(:,:,:), RH(:,:,:), T600(:,:,:)
	REAL*8, ALLOCATABLE :: VPOTC(:,:,:)
	REAL, ALLOCATABLE :: X1VERT(:), Y1VERT(:)
	REAL, ALLOCATABLE :: X2VERT(:), Y2VERT(:)
c
	REAL AMP(15), PHASEU850(15), PHASEV850(15)
	REAL PHASEU250(15),PHASEV250(15)
	REAL A(4,4), PCH(4)
	REAL AMPEC(15),PHASEUEC(15),PHASEVEC(15)
	REAL AEC(2,2), PCHEC(2)
	REAL LONGINC, LATINC, LONGSTART, LATSTART
C
C  ***   Print quantities   **
C
	INTEGER MONTHPR(125),HOURPR(125),DAYPR(125)
        REAL LATPR(125),LONGPR(125),SHEARPR(125)
C
	INTEGER MONTH, DATE
C
C  ***       Other quantities    ***
C
	REAL*8 TT,DTT,TIME, ALAT, ALOG,PRINTIME,PIAMP
	REAL LATMAX, LOGMAX, LOGMIN, LATMING, LATMAXG
	real lat_init,long_init
	INTEGER DAY(12)
	INTEGER ISEED, year
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), XSURST(NRD)
	REAL MU1(NRD), MU2(NRD), MU3(NRD), RTS3(NRD)
C
C  ***     Dimension ocean variables   ***
C
	REAL HMIX(NRD), UHMIX1(NRD), UHMIX2(NRD), UHMIX3(NRD), HM
	REAL DISTANCE(600),UTIME(600), LATG(600), LOGG(600)
	REAL HMIX0(600),BATHTRACK(600),VPTRACK(600)
	REAL STRAT(600),BATHVAR(600),RHTRACK(600),T600TRACK(600)
C
C  ***     Climatological ocean and potential intensity files  ***
C
	REAL BATH(1440,721),MIXDEPTH(12,180,360)
	REAL THERMAL(12,180,360)
	REAL LATTE(200),LOGGTE(200),VSTE(200),VSTES(200)
	INTEGER MONTE(200),DATETE(200),WRITETIME(200)
	REAL DISTE(200),UTE(200),VSTEG(200),UTES(200)
C
C  ***      Shear-related files  ***
C
	REAL SHEARTE(200), SHEARTRACK(600)
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 RMG(NTG), RTMG(NTG), LATTIME(NTG), LOGTIME(NTG)
	REAL HMG(NTG), UTRANS(NTG), VTRANS(NTG)
	REAL SSTT(NTG), SSTR(NRD), HMG0(NTG)
	REAL VMG(NTG)
	REAL VSHG(NTG), VPOTG(NTG)
	REAL MUG(NTG), MDG(NTG), PCG(NTG,2)
        REAL*8 TIMEG(NTG)
C
C   ***   Character arrays    ***
C
	CHARACTER*2 AYEAR,AMON,ADAY,AHOUR,BASIN,ISTNUM,DATEC,MOND
	CHARACTER*2 ABAS,ABASU,ASNUM,ABAS2,SHOUR
	CHARACTER*3 MONC,DATETEMP
	CHARACTER*1 LATN,LOGGB,PSHEAR,ILATI,ILONGI,AK
	CHARACTER*10 ndate, DATEA
	CHARACTER*1 llat,llong
	CHARACTER*80 JFILE,FNAMEA,FNAMEB,PATHBUF
	CHARACTER*95 NAM
	CHARACTER*30 NAM2
	CHARACTER*20 ANAM,MODEL,STORM
	CHARACTER*11 SNAME
	CHARACTER*15 SDATE
	CHARACTER*80 DIREC3
	CHARACTER*120 DIRECTORY
	CHARACTER*5  FN
	CHARACTER*4 IDENT, COMPFILE, FCST, FN2, PIFILE, FCSTO,AYEAR4
	CHARACTER*4 FCST2, gmeth
	CHARACTER*3 lev
	CHARACTER*50 SFILENAME
        CHARACTER GENPOINTS*8
	CHARACTER AF1*5, CITYNAME*90, AMON1*1, BAS*2, BASlc*2
	CHARACTER ch1*1,ch2*2,ch3*3,ch4*4,ch5*5
C
C  ***       Other quantities    ***
C
	REAL H_A, MEQ, MF, MT, MUMAX, MDMIN
	REAL*8 DT,ATIME,DT1,DTG,TT1,DT0
	CHARACTER*4 FMT,OM
	integer*8 adate(3)
c
	data day/31,28,31,30,31,30,31,31,30,31,30,31/
C
C  ***   Read in city parameters  ***
C
	OPEN(UNIT=11,FILE='wrt_pred.in',STATUS='OLD')
	READ(11,*)BAS,MODEL,DATEA,STORM,imoni,
     1   idatei,ihouri,NUMTRACKS,VMCRIT,dtotr,rerandom,adate(1)
	CLOSE(11)
	adate(2)=1
	adate(3)=1
	if(long_init.lt.0.0)long_init=long_init+360.0
	gmeth='rand'
	clat=0.0
	clog=180.0
	crad=100000.0
	IF(BAS.EQ.'AL')BASlc='al'
        IF(BAS.EQ.'EP')BASlc='ep'
	IF(BAS.EQ.'WP')BASlc='wp'
	IF(BAS.EQ.'IO')BASlc='io'
	IF(BAS.EQ.'SH')BASlc='sh'
c
C  Read in B deck and find initial storm intensity and position
c
	v_init=0.0
	lv_init=0
	OPEN(UNIT=11,FILE='b'//BASlc//TRIM(STORM(3:4))//
     1   TRIM(DATEA(1:4))//'.dat',STATUS='OLD')
	do i=1,1000
	 v_init_old=float(lv_init)
	 read(11,18,end=3)ndate,llat_init,llat,llong_init,llong,lv_init
	 if(ndate.eq.datea)then
	  lat_init=0.1*float(llat_init)
	  long_init=0.1*float(llong_init)
	  v_init=float(lv_init)
	  if(llat.eq.'S')lat_init=-lat_init
	  if(llong.eq.'W')long_init=360.0-long_init
	  goto 3
	 end if
	end do
   3    continue
	close(11)
  18    FORMAT(8X,A10,17X,I3,A1,2X,I4,A1,2X,I3)
c
C  VMCRIT is the critical maximum wind (kts) within city radius below which track is ignored
C
	IF(CLOG.LT.0.0)CLOG=CLOG+360.0
c
	clog1=clog
	clat1=clat
	jint=0
c
c  ***    The follwing are attributes of particular models....also see PARAMETER statment   ***
c
	DIREC3='data/'//trim(model)//'/'//trim(datea)//'/'//trim(storm)
c
c       Grid and other specifications
c
c  ***  The parameter PIAMP is adjusted to bring the 20th century potential intensity more in line with obs  ***
c
	lev='250'
        NLO=181
        NLA=91 
	NN=21
        LONGINC=2.0    
        LATINC=2.0    
        LONGSTART=0.0
        LATSTART=-90.0
        PIAMP=0.8
c
	ALLOCATE(U850(NLO,NLA,NN),U250(NLO,NLA,NN),V850(NLO,NLA,NN))
        ALLOCATE(V250(NLO,NLA,NN),UVAR250(NLO,NLA,NN))
	ALLOCATE(UVAR850(NLO,NLA,NN),VVAR250(NLO,NLA,NN))
	ALLOCATE(VVAR850(NLO,NLA,NN),U250V250(NLO,NLA,NN))
	ALLOCATE(U250U850(NLO,NLA,NN),V250V850(NLO,NLA,NN))
	ALLOCATE(V850U850(NLO,NLA,NN),VOR(NLO,NLA,NN))
	ALLOCATE(VTEM(NLO,NLA,NN),RH(NLO,NLA,NN),VPOTC(NLO,NLA,NN))
	ALLOCATE(T600(NLO,NLA,NN))
	ALLOCATE(UEC(NN),VEC(NN),U2EC(NN),V2EC(NN),UVEC(NN))
c
	CITYNAME=TRIM(DATEA)//'/'//TRIM(STORM)
c
      OPEN(UNIT=11, FILE='input.in',status='old')
C
C           ***  Read in PARAMETERS  ***
C
      READ(11,13)FN, COMPFILE, PIFILE, PSHEAR
   13 FORMAT(/////,T60,A8,/,T60,A4,/,T60,A4,/,T60,A1,/)
      READ(11,15)RM,R0,VOFFSET,VINOFFSET,VSMIN
   15 FORMAT(//,4(T57,F9.4,/),/,T57,F9.4,/)
      READ(11,20)NR,DT,RO
   20 FORMAT(///,T62,I3,/,2(T57,F9.4,/))
      READ(11,22)OM
   22 FORMAT(T62,A1)
      READ(11,25)DTG,ROG,FMT
   25 FORMAT(///,2(T57,F9.4,/)/,T62,A3)
C
       CLOSE(UNIT=11)
	DTG=2.0D0
c
	IF(PSHEAR.EQ.'Y')PSHEAR='y'
	IF(PSHEAR.EQ.'N')PSHEAR='n'
C
C  ***   A few constants   ***
C
C  ***   Beta drift   ***
C
	UDRIFT=0.0 
        VDRIFT=3.0
c	VDRIFT=0.0
C
C  ***    Time scale in Fourier series wind fields, and weight given             ***
C  ***      to 250 and 850 hPa flow in setting storm motion                      ***
C 
        TAUMAX=15.0
	WEIGHT=0.75
c
c  ***    Global frequency factor (empirical, matched to present climate)       ***
c  ***             Used only for random seeding method                          ***
c
        globalfac=22.0
C
C  ***  Latitude and longitude boundaries on tracks (only ABS(LATMAX)) used      ***
C
c
c   The following bounds for genesis locations are used in all basins except global
c
        LATMING=3.0
        LATMAXG=45.0
c
c  ***      Definitions of basin bounds and frequency factors for climo. These have been             ***
c  ***  calibrated to give the observed frequencies averaged between 1970 and 2005 in all basins     ***
c  ***          for all storms whose maximum winds exceed 40 knots. Driven by NCEP reanalysis.       ***
c
	IF(BAS.EQ.'AL')THEN
	 LATMAX=40.0
	 LOGMAX=355.0
	 LOGMIN=260.0
	 LATMAXT=65.0
	 LOGMAXT=370.0
	 LOGMINT=240.0
	 climfac=46.23 
	ELSE IF(BAS.EQ.'EP')THEN
         LATMAX=35.0
         LOGMAX=285.0
         LOGMIN=200.0
	 LATMAXT=60.0
	 LOGMAXT=300.0
	 LOGMINT=100.0
	 climfac=65.93 
        ELSE IF(BAS.EQ.'WP')THEN
         LATMAX=55.0
         LOGMAX=200.0
         LOGMIN=100.0 
	 LATMAXT=65.0
	 LOGMAXT=260.0
	 LOGMINT=50.0
	 climfac=53.5 
        ELSE IF(BAS.EQ.'IO')THEN
         LATMAX=32.0
         LOGMAX=100.0
         LOGMIN=40.0 
	 LATMAXT=40.0
	 LOGMAXT=150.0
	 LOGMINT=20.0
	 climfac=15.08
        ELSE IF(BAS.EQ.'SH'.OR.BAS.EQ.'BM')THEN
         LATMAX=50.0
         LOGMAX=290.0
c	 LOGMAX=380.0
         LOGMIN=15.0 
	 LATMAXT=60.0
	 LOGMAXT=300.0
	 LOGMINT=5.0
	 climfac=91.8
	 if(BAS.EQ.'BM')climfac=39.93
        ELSE IF(BAS.EQ.'SA')THEN
         LATMAX=50.0
         LOGMAX=380.0
         LOGMIN=300.0
	 LATMAXT=60.0
	 LOGMAXT=380.0
	 LOGMINT=290.0
        ELSE IF(BAS.EQ.'GB')THEN
         LATMAX=75.0
	 LATMAXG=75.0
	 LATMING=-65.0
         LOGMAX=360.0
         LOGMIN=0.0
         LATMAXT=80.0
         LOGMAXT=360.0
         LOGMINT=0.0
        END IF
c
c  ***   Revise latitude bound     ***
c
        latmaxg=min(latmaxg,latmax)
c
c  ***   Calculate area of basin    ***
c
        PI=ACOS(-1.0)
        ADEG=PI/180.
        areafac=(sin(adeg*latmaxg)-sin(adeg*latming))*(logmax-logmin)
        IF(BAS.EQ.'AL')then
         areafac=areafac-(sin(adeg*15.0)-sin(adeg*latming))*
     1     (275.0-logmin)
        END IF
        IF(BAS.EQ.'EP')THEN
         areafac=areafac-(sin(adeg*latmaxg)-sin(adeg*16.0))*
     1    (logmax-262.0)
        END IF
C
        ENDTIME=10.*24.*3600.    
	DELPRINT=6.*3600.
	flfac=PI/(6.0*3600.0)
	AFAC2=ADEG*6.38E6
	DISFACTOR=1./(60.*1852.)
	COEFS=2.*PI/(TAUMAX*24.*3600.)
	DTT=30.0D0*60.0D0
C
	IFLAG=0
c	ICOUNT=0
C
C   ***  Set default values of certain parameters
C
        TS=27.0
        TO=-70.0
        H_A=0.8
        PA=1005.0
        FC=5.0
        CECD=1.0
        CD=0.8
        CD1=3.0
        PNU=0.03
        TAUR=8.0
        RADMAX=2.0
        TAUC=4.0
        EFRAC=0.5
        DPB=50.0
        AHM=0.6
        DSST=0.6
        GM=8.0
        TRANSFAC=0.7
        CDM=3.0
        CDMAX=1.8
        CKMAX=2.0
        VCD=30.0
        IVMAXMIN=15
        CKEXP=2.0
        IVSMIN=VSMIN*3600./1852.
C
C   ***  More external parameters  ***
C
	XMDIF=750.0
	RIC=1.0
	NRESTART=0
C
	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=FC*1.0E-5
	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)
	TM1=TM+3.0
	ESM1=6.112*EXP(17.67*TM1/(243.5+TM1))
	QSM1=0.622*ESM1/(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
c
	RMFAC=SQRT(20.0/v_init)
	RMFAC=MIN(RMFAC,1.0)
	RM=RM*RMFAC
c
	RM=RM*1000./ALENGTH
	RMINIT=RM
	RMBE=RM
	RO=RO*1000./ALENGTH 
	ROG=ROG*1000./ALENGTH
	DT=DT/ATIME
	DT0=DT
	DTG=DTG*3600./ATIME
	AKNOTFAC=1852.0/(3600.0*SCHI)
c
        VM=v_init*aknotfac
c
	VMINIT0=VM
	VMINITOLD=v_init_old*aknotfac
	VMINITNEW=2.0*VM-VMINITOLD
	VMINITNEW=MAX(VMINITNEW,20.0*AKNOTFAC)
	VOFFSETND=VOFFSET/SCHI
        VINOFFSET=VINOFFSET/(3600.*24.)
	XM0=2.5E6*(TS-TO)*QSM*(AHM-1.)/TMA
	XM0=XM0/CHI
	CDV=CD1*SCHI/CD
	RHAA=H_A
	GM=GM*0.01
	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
	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)
        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
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
	BFAC1=1.-DR*SQRT(2./Q)
	PBM=1.0E-5
	SIXTI=1./16.
C
C              *** SET ACTUAL DIFFUSIVITIES PROPORTIONAL TO DR**2    ***
C
	PNU=PNU*PNU*DR*0.01744
C
C   ***    Specify initial intensification rate  ***
C
        IF(gmeth.eq.'clim')THEN
c         VINTEN=1.5/(6.*3600.)
	 VINTEN=100.0
        ELSE
         VINTEN=100.0
        END IF
c
	VINTEN0=VINTEN
C
C  ***   Read in monthly mean wind data  ***
C
	OPEN(11,file=TRIM(DIREC3)//'/u850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)U850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u'//lev,
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)U250
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/v850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)V850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/v'//lev,
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)V250
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u850u850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)UVAR850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u'//lev//'u'//lev,
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)UVAR250
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/v850v850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)VVAR850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/v'//lev//'v'//lev,
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)VVAR250
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u'//lev//'v'//lev,
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)U250V250
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u'//lev//'u850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)U250U850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u850v850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)V850U850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/v'//lev//'v850',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)V250V850
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/uec',
     1	 status='old', access='direct',recl=8*NN)
	READ(11,REC=1)UEC 
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/vec',
     1	 status='old', access='direct',recl=8*NN)
	READ(11,REC=1)VEC 
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/u2ec',
     1	 status='old', access='direct',recl=8*NN)
	READ(11,REC=1)U2EC 
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/v2ec',
     1	 status='old', access='direct',recl=8*NN)
	READ(11,REC=1)V2EC 
	CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/uvec',
     1	 status='old', access='direct',recl=8*NN)
	READ(11,REC=1)UVEC 
	CLOSE(11)
C
C   ***   Calculate vorticity   ***
C
	do k=1,NN
	 do j=1,nla
	  jp=j+1
	  jp=min(jp,nla)
	  jm=j-1
	  jm=max(jm,1)
	  coslat=cos((-90.+latinc*float(j-1))*adeg)
	  coslatp=cos((-90.+latinc*float(jp-1))*adeg)
	  coslatm=cos((-90.+latinc*float(jm-1))*adeg)
	  do i=1,nlo
	   ip=i+1
	   if(i.eq.nlo)ip=1
	   im=i-1
	   if(i.eq.1)im=nlo
	   vor(i,j,k)=((u850(i,jm,k)*coslatm-u850(i,jp,k)*coslatp)/
     1       (coslat*latinc)
     2      +(v850(ip,j,k)-v850(im,j,k))/longinc)/(2.0*111120.) 
	  end do
	 end do
	end do
C
C   ***  Read in potential intensity, humidity, and mixed layer files   ***
C
	OPEN(11,file=TRIM(DIREC3)//'/rh',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)RH
	CLOSE(11)
C
        OPEN(11,file=TRIM(DIREC3)//'/T600',
     1   status='old', access='direct',recl=8*NLO*NLA*NN)
        READ(11,REC=1)T600
        CLOSE(11)
C
	OPEN(11,file=TRIM(DIREC3)//'/pi',
     1	 status='old', access='direct',recl=8*NLO*NLA*NN)
	READ(11,REC=1)VTEM
	CLOSE(11)
c
c  ***   Added factor of 1.1 on 2/20/2007 to bring histograms more in line with obs ***
c  ***          If statement added to deal with NaNs in VPOTC file                  ***
c

      DO 1380 K=1,NN
c
c  Following statements added 9/30/2020 to lag potential intensity fields
c
	IF(K.GE.5)THEN
	 KLAG=K-4
	ELSE
	 KLAG=K+4
	END IF
c
        DO 1350 M=1,NLO
         DO 1350 N1=1,NLA
	    IF(VTEM(M,N1,K).GT.0.0.AND.VTEM(M,N1,K).LE.200.0)THEN
             VPOTC(M,N1,K)=PIAMP*1.17*MAX(VTEM(M,N1,KLAG),0.0D0)
c	     VPOTC(M,N1,K)=VPOTC(M,N1,K)+0.1682*
c     1        FLOAT(MAX((year-1976),0))
	    ELSE
	     VPOTC(M,N1,K)=0.0
	    END IF
 1350   CONTINUE
C
 1380	CONTINUE
C
       DO 1400 K=1,12
         WRITE(AMON,1310)K
	 OPEN(11,FILE='data/ocean/mix'//AMON, STATUS='OLD')
	 OPEN(12,FILE='data/ocean/stratnorm'//AMON, STATUS='OLD')
	 DO 1370 I=1,180
	  READ(11,*)(MIXDEPTH(K,I,J),J=1,360)
c
c     Variants on mixed layer depth
c
	 IF(DATEA(1:2).EQ.'eo'.OR.STORM(1:2).EQ.'pl')THEN
	  DO J=1,360
c	   MIXDEPTH(K,I,J)=MIXDEPTH(K,I,J)*1.1
	   MIXDEPTH(K,I,J)=100.0
	  END DO
	 END IF
c
	  READ(12,*)(THERMAL(K,I,J),J=1,360)
 1370	 CONTINUE 
 1400 CONTINUE
	close(11)
	close(12)
 1310 FORMAT(I2.2)
 1320 FORMAT(I2)
c
C   ***  Read in bathymetric file   ***
C
        NAM='data/ocean/bathsmooth2.dat'
        OPEN(11,FILE=NAM,status='old',
     1    access='direct',recl=4*1440*721)
        READ(11,REC=1)BATH
	close(11)
C
C           ***  Re-initialize random numbers  ***
C
        if(rerandom.eq.1)then
          call random_seed ( PUT = adate(1:3))
        end if
C
C           ***  Randomly generate genesis data  ***
C
c        iseed=time8 ()
c	 iseed=1365079
c	 iseed=5869487
	 iseed=9735681
c        iseed1=iseed-10*int(1.0e-1*float(iseed))
c        iseed2=iseed-1e2*int(1.0e-2*float(iseed))
c	iseed=100*iseed2+10*iseed1+iseed2+0.1*iseed1
c	iseed=iseed2
c	iseed=iseed-100*(iseed/100)
c	iseed=1e8*iseed+1e6*iseed+1e4*iseed
c
c	open(unit=11,file='irandom',status='old')
c	 READ(11,*) ise
c	 ise=1e5*ise
c	CLOSE(11)
c
c	gdum=rand(1)
c        call srand(iseed)
c	gdum2=rand(0)
c	call random_seed()
	call random_number(gdum)
	call random_number(gdum2)
c
	M=0
	NTOTAL=0
	NUMSTORMS=0
c
	IF(gmeth.eq.'clim')THEN
	  call random_seed()
	END IF
C
    5   CONTINUE
C
    7   continue
C
	NTOTAL=NTOTAL+1
	IFLAG=0
c
	call random_number(arand)
	AILAT=lat_init+0.3*(arand-0.5)
	call random_number(arand)
	AILOG=long_init+0.3*(arand-0.5)
c
         JLAT=1+(AILAT-LATSTART)/LATINC
         IF(AILOG.GT.0.0)THEN
          JLOGG=1+(AILOG-LONGSTART)/LONGINC
         ELSE
          JLOGG=1+(360.+AILOG-LONGSTART)/LONGINC
         END IF
         JLATP=JLAT+1
         JLOGGP=JLOGG+1
         IF(JLOGGP.GT.NLO)JLOGGP=1
C
	SIGN=1.0
	IF(AILAT.LT.0.0)SIGN=-1.0
	VDRIFTNEW=VDRIFT*SIGN
	NSTEPS=0
C
	ALAT=AILAT
	ALOG=AILOG
	IMON=IMONI
	IDATE=IDATEI
	IHOUR=IHOURI
C
	TT=0.
	RMIN=1.0E14
	JTIME=0
	PRINTIME=-DTT
	TIME=3600.*FLOAT(IHOUR)
C
	SUMM=0.0
	DO I=1,15
	 BI=FLOAT(I)
	 SUMM=SUMM+1./(BI**3)
	 AMP(I)=1./(BI**1.5)
c	 PHPS=RAND(0) 
	 call random_number(PHPS)
	 PHASEU850(I)=2.*PI*(PHPS-0.5)
c	 PHPS=RAND(0) 
	 call random_number(PHPS)
	 PHASEV850(I)=2.*PI*(PHPS-0.5)
c         PHPS=RAND(0)
	 call random_number(PHPS)
         PHASEU250(I)=2.*PI*(PHPS-0.5)
c         PHPS=RAND(0)
	 call random_number(PHPS)
         PHASEV250(I)=2.*PI*(PHPS-0.5)
	 call random_number(PHPS)
         PHASEUEC(I)=2.*PI*(PHPS-0.5)
	 call random_number(PHPS)
         PHASEVEC(I)=2.*PI*(PHPS-0.5)
	END DO
	SUMM=SQRT(2./SUMM)
C
   10 FORMAT(I2.2)
   11 FORMAT(I1.1)
c
	do i=1,4
	 do jj=1,4
	  A(i,jj)=0.0
	 end do
	end do
c
	do i=1,2
	 do jj=1,2
	  AEC(i,jj)=0.0
	 end do
	end do
C
C      ***  Time march storm position   ***
C
  101 CONTINUE
C
	TT=TT+DTT
	PRINTIME=PRINTIME+DTT
C
	TIME=TIME+DTT
	IF(TIME.GT.24.*3600.)THEN
	 TIME=TIME-24.*3600.
	 IDATE=IDATE+1
	    IF((IMON.EQ.1.OR.IMON.EQ.3.OR.IMON.EQ.5.OR.IMON.EQ.7.OR.
     1     IMON.EQ.8.OR.IMON.EQ.10).AND.IDATE.GT.31)THEN
		  IDATE=1
	      IMON=IMON+1
	    END IF
	    IF((IMON.EQ.4.OR.IMON.EQ.6.OR.IMON.EQ.9.OR.IMON.EQ.11)
     1		.AND.IDATE.GT.30)THEN
		  IDATE=1
	      IMON=IMON+1
	    END IF
	    IF(IMON.EQ.2.AND.IDATE.GT.28)THEN
	      IDATE=1
	      IMON=IMON+1
	    END IF
	    IF(IMON.GT.12)IMON=1
	END IF
	IHOUR=TIME/3600.
	IMIN=TIME/60.
C
	 JLAT=1+(ALAT-LATSTART)/LATINC
         IF(ALOG.GT.0.0)THEN
          JLOGG=1+(ALOG-LONGSTART)/LONGINC
         ELSE
          JLOGG=1+(360.+ALOG-LONGSTART)/LONGINC
         END IF
	 JLATP=JLAT+1
	 JLOGGP=JLOGG+1
	 IF(JLOGGP.GT.NLO)JLOGGP=1
C
          KTIME=1+INT(2.*dtotr+TT/(12.*3600.))
          KTIMEP=KTIME+1
	  KTIME=MIN(KTIME,NN)
	  KTIMEP=MIN(KTIMEP,NN)
          FAC=TT/(12.*3600.)-INT(TT/(12.*3600.))
c
c	  if(ktime.lt.4)then
c	   ktime=4
c	   ktimep=4
c	  end if
C
           DELY=1.+(ALAT-LATSTART)/LATINC-FLOAT(JLAT)
           DELX=1.+(360.+ALOG-LONGSTART)/LONGINC-
     1      INT(1.+(360.+ALOG-LONGSTART)/LONGINC)
           D1=(1.-DELX)*(1.-DELY)
           D2=DELY*(1.-DELX)
           D3=DELX*(1.-DELY)
           D4=DELX*DELY
C
	   SU1=(1.-FAC)*U850(JLOGG,JLAT,KTIME)+
     1	    FAC*U850(JLOGG,JLAT,KTIMEP)
	   SU2=(1.-FAC)*U850(JLOGG,JLATP,KTIME)+FAC*
     1      U850(JLOGG,JLATP,KTIMEP)
	   SU3=(1.-FAC)*U850(JLOGGP,JLAT,KTIME)+FAC*
     1      U850(JLOGGP,JLAT,KTIMEP)
	   SU4=(1.-FAC)*U850(JLOGGP,JLATP,KTIME)+FAC*
     1      U850(JLOGGP,JLATP,KTIMEP)
	   SV1=(1.-FAC)*V850(JLOGG,JLAT,KTIME)+
     1	    FAC*V850(JLOGG,JLAT,KTIMEP)
	   SV2=(1.-FAC)*V850(JLOGG,JLATP,KTIME)+FAC*
     1      V850(JLOGG,JLATP,KTIMEP)
	   SV3=(1.-FAC)*V850(JLOGGP,JLAT,KTIME)+FAC*
     1      V850(JLOGGP,JLAT,KTIMEP)
	   SV4=(1.-FAC)*V850(JLOGGP,JLATP,KTIME)+FAC*
     1      V850(JLOGGP,JLATP,KTIMEP)
C
	  U850TRACK=D1*SU1+D2*SU2+D3*SU3+D4*SU4
	  V850TRACK=D1*SV1+D2*SV2+D3*SV3+D4*SV4
C
	   SU1=(1.-FAC)*U250(JLOGG,JLAT,KTIME)+
     1	    FAC*U250(JLOGG,JLAT,KTIMEP)
	   SU2=(1.-FAC)*U250(JLOGG,JLATP,KTIME)+FAC*
     1      U250(JLOGG,JLATP,KTIMEP)
	   SU3=(1.-FAC)*U250(JLOGGP,JLAT,KTIME)+FAC*
     1      U250(JLOGGP,JLAT,KTIMEP)
	   SU4=(1.-FAC)*U250(JLOGGP,JLATP,KTIME)+FAC*
     1      U250(JLOGGP,JLATP,KTIMEP)
	   SV1=(1.-FAC)*V250(JLOGG,JLAT,KTIME)+
     1	    FAC*V250(JLOGG,JLAT,KTIMEP)
	   SV2=(1.-FAC)*V250(JLOGG,JLATP,KTIME)+FAC*
     1      V250(JLOGG,JLATP,KTIMEP)
	   SV3=(1.-FAC)*V250(JLOGGP,JLAT,KTIME)+FAC*
     1      V250(JLOGGP,JLAT,KTIMEP)
	   SV4=(1.-FAC)*V250(JLOGGP,JLATP,KTIME)+FAC*
     1      V250(JLOGGP,JLATP,KTIMEP)
C
	  U250TRACK=D1*SU1+D2*SU2+D3*SU3+D4*SU4
	  V250TRACK=D1*SV1+D2*SV2+D3*SV3+D4*SV4
C
           SU1=(1.-FAC)*U250V250(JLOGG,JLAT,KTIME)+
     1      FAC*U250V250(JLOGG,JLAT,KTIMEP)
           SU2=(1.-FAC)*U250V250(JLOGG,JLATP,KTIME)+FAC*
     1      U250V250(JLOGG,JLATP,KTIMEP)
           SU3=(1.-FAC)*U250V250(JLOGGP,JLAT,KTIME)+FAC*
     1      U250V250(JLOGGP,JLAT,KTIMEP)
           SU4=(1.-FAC)*U250V250(JLOGGP,JLATP,KTIME)+FAC*
     1      U250V250(JLOGGP,JLATP,KTIMEP)
C
          U250V250TRACK=D1*SU1+D2*SU2+D3*SU3+D4*SU4
C
           SU1=(1.-FAC)*U250U850(JLOGG,JLAT,KTIME)+
     1      FAC*U250U850(JLOGG,JLAT,KTIMEP)
           SU2=(1.-FAC)*U250U850(JLOGG,JLATP,KTIME)+FAC*
     1      U250U850(JLOGG,JLATP,KTIMEP)
           SU3=(1.-FAC)*U250U850(JLOGGP,JLAT,KTIME)+FAC*
     1      U250U850(JLOGGP,JLAT,KTIMEP)
           SU4=(1.-FAC)*U250U850(JLOGGP,JLATP,KTIME)+FAC*
     1      U250U850(JLOGGP,JLATP,KTIMEP)
C    
          U250U850TRACK=D1*SU1+D2*SU2+D3*SU3+D4*SU4
C    
           SU1=(1.-FAC)*V250V850(JLOGG,JLAT,KTIME)+
     1      FAC*V250V850(JLOGG,JLAT,KTIMEP)
           SU2=(1.-FAC)*V250V850(JLOGG,JLATP,KTIME)+FAC*
     1      V250V850(JLOGG,JLATP,KTIMEP)
           SU3=(1.-FAC)*V250V850(JLOGGP,JLAT,KTIME)+FAC*
     1      V250V850(JLOGGP,JLAT,KTIMEP)
           SU4=(1.-FAC)*V250V850(JLOGGP,JLATP,KTIME)+FAC*
     1      V250V850(JLOGGP,JLATP,KTIMEP)
C   
          V250V850TRACK=D1*SU1+D2*SU2+D3*SU3+D4*SU4
C
           SU1=(1.-FAC)*V850U850(JLOGG,JLAT,KTIME)+
     1      FAC*V850U850(JLOGG,JLAT,KTIMEP)
           SU2=(1.-FAC)*V850U850(JLOGG,JLATP,KTIME)+FAC*
     1      V850U850(JLOGG,JLATP,KTIMEP)
           SU3=(1.-FAC)*V850U850(JLOGGP,JLAT,KTIME)+FAC*
     1      V850U850(JLOGGP,JLAT,KTIMEP)
           SU4=(1.-FAC)*V850U850(JLOGGP,JLATP,KTIME)+FAC*
     1      V850U850(JLOGGP,JLATP,KTIMEP)
C   
          V850U850TRACK=D1*SU1+D2*SU2+D3*SU3+D4*SU4
C
	   UV1=(1.-FAC)*UVAR850(JLOGG,JLAT,KTIME)+
     1	    FAC*UVAR850(JLOGG,JLAT,KTIMEP)
	   UV2=(1.-FAC)*UVAR850(JLOGG,JLATP,KTIME)+FAC*
     1      UVAR850(JLOGG,JLATP,KTIMEP)
	   UV3=(1.-FAC)*UVAR850(JLOGGP,JLAT,KTIME)+FAC*
     1      UVAR850(JLOGGP,JLAT,KTIMEP)
	   UV4=(1.-FAC)*UVAR850(JLOGGP,JLATP,KTIME)+FAC*
     1      UVAR850(JLOGGP,JLATP,KTIMEP)
	   VV1=(1.-FAC)*VVAR850(JLOGG,JLAT,KTIME)+
     1	    FAC*VVAR850(JLOGG,JLAT,KTIMEP)
	   VV2=(1.-FAC)*VVAR850(JLOGG,JLATP,KTIME)+FAC*
     1      VVAR850(JLOGG,JLATP,KTIMEP)
	   VV3=(1.-FAC)*VVAR850(JLOGGP,JLAT,KTIME)+FAC*
     1      VVAR850(JLOGGP,JLAT,KTIMEP)
	   VV4=(1.-FAC)*VVAR850(JLOGGP,JLATP,KTIME)+FAC*
     1      VVAR850(JLOGGP,JLATP,KTIMEP)
C
	  UVAR850TRACK=D1*UV1+D2*UV2+D3*UV3+D4*UV4
	  VVAR850TRACK=D1*VV1+D2*VV2+D3*VV3+D4*VV4
C
	   UV1=(1.-FAC)*UVAR250(JLOGG,JLAT,KTIME)+
     1	    FAC*UVAR250(JLOGG,JLAT,KTIMEP)
	   UV2=(1.-FAC)*UVAR250(JLOGG,JLATP,KTIME)+FAC*
     1      UVAR250(JLOGG,JLATP,KTIMEP)
	   UV3=(1.-FAC)*UVAR250(JLOGGP,JLAT,KTIME)+FAC*
     1      UVAR250(JLOGGP,JLAT,KTIMEP)
	   UV4=(1.-FAC)*UVAR250(JLOGGP,JLATP,KTIME)+FAC*
     1      UVAR250(JLOGGP,JLATP,KTIMEP)
	   VV1=(1.-FAC)*VVAR250(JLOGG,JLAT,KTIME)+
     1	    FAC*VVAR250(JLOGG,JLAT,KTIMEP)
	   VV2=(1.-FAC)*VVAR250(JLOGG,JLATP,KTIME)+FAC*
     1      VVAR250(JLOGG,JLATP,KTIMEP)
	   VV3=(1.-FAC)*VVAR250(JLOGGP,JLAT,KTIME)+FAC*
     1      VVAR250(JLOGGP,JLAT,KTIMEP)
	   VV4=(1.-FAC)*VVAR250(JLOGGP,JLATP,KTIME)+FAC*
     1      VVAR250(JLOGGP,JLATP,KTIMEP)
C
	  UVAR250TRACK=D1*UV1+D2*UV2+D3*UV3+D4*UV4
	  VVAR250TRACK=D1*VV1+D2*VV2+D3*VV3+D4*VV4
C
	  UECTRACK=(1.-FAC)*UEC(KTIME)+FAC*UEC(KTIMEP)
	  VECTRACK=(1.-FAC)*VEC(KTIME)+FAC*VEC(KTIMEP)
	  U2ECTRACK=(1.-FAC)*U2EC(KTIME)+FAC*U2EC(KTIMEP)
	  V2ECTRACK=(1.-FAC)*V2EC(KTIME)+FAC*V2EC(KTIMEP)
	  UVECTRACK=(1.-FAC)*UVEC(KTIME)+FAC*UVEC(KTIMEP)
C
	U850SUM=0.0
	V850SUM=0.0
        U250SUM=0.0
        V250SUM=0.0
	UECSUM=0.0
	VECSUM=0.0
	DO I=1,15
	 U850SUM=U850SUM+AMP(I)*SIN(COEFS*FLOAT(I)*TT+PHASEU850(I))
	 V850SUM=V850SUM+AMP(I)*SIN(COEFS*FLOAT(I)*TT+PHASEV850(I))
         U250SUM=U250SUM+AMP(I)*SIN(COEFS*FLOAT(I)*TT+PHASEU250(I))
         V250SUM=V250SUM+AMP(I)*SIN(COEFS*FLOAT(I)*TT+PHASEV250(I))
	 UECSUM=UECSUM+AMP(I)*SIN(COEFS*FLOAT(I)*TT+PHASEUEC(I))
	 VECSUM=VECSUM+AMP(I)*SIN(COEFS*FLOAT(I)*TT+PHASEVEC(I))
	END DO
C
C   This is new and accounts for covariances between u850 and v850.  5/19/05
C      Uses cholesky decomposition from numerical recipes
C
        A(1,1)=UVAR250TRACK
        A(1,2)=U250V250TRACK
        A(1,3)=U250U850TRACK
        A(1,4)=0.0
        A(2,2)=VVAR250TRACK
        A(2,3)=0.0
        A(2,4)=V250V850TRACK
        A(3,3)=UVAR850TRACK
        A(3,4)=V850U850TRACK
        A(4,4)=VVAR850TRACK
c
	AEC(1,1)=U2ECTRACK
	AEC(1,2)=UVECTRACK
	AEC(2,2)=V2ECTRACK
c
        CALL choldc(A,4,4,PCH)
c
	CALL choldc(AEC,2,2,PCHEC)
C
        U250M=U250TRACK+SUMM*PCH(1)*U250SUM
        V250M=V250TRACK+SUMM*(A(2,1)*U250SUM+PCH(2)*V250SUM)
        U850M=U850TRACK+SUMM*(A(3,1)*U250SUM+PCH(3)*U850SUM)
        V850M=V850TRACK+SUMM*(A(4,2)*V250SUM+A(4,3)*U850SUM+
     1    PCH(4)*V850SUM)
c
	UECMEAN=UECTRACK+SUMM*PCHEC(1)*UECSUM
	VECMEAN=VECTRACK+SUMM*(AEC(2,1)*UECSUM+PCHEC(2)*VECSUM)
C
	UMEAN=WEIGHT*U850M+(1.-WEIGHT)*U250M+UDRIFT*COS(ALAT*ADEG)
c	UMEAN=WEIGHT*U850M+(1.-WEIGHT)*U250M-0.5*COS(ALAT*ADEG)
c
c  Added 6/16/2010
c
	VDRIFTMOD=VDRIFTNEW*MIN((TT/(3.0*24.0*3600.0)),1.0D0)
c
	VMEAN=WEIGHT*V850M+(1.-WEIGHT)*V250M+VDRIFTMOD*COS(ALAT*ADEG)
	UNET=U250M-U850M
	VNET=V250M-V850M
	VSHEAR=SQRT(UNET**2+VNET**2)
C
C  ***  Variables for adding inertia to storm tracks, added 2/25/2010
C
	IF(TT.LE.(1.5*DTT))THEN
	  UM1=UMEAN
	  VM1=VMEAN
	  UM1=-0.0
	  VM1=3.0
	  UM2=UM1
	  VM2=VM1
	  TAUDAMP=1.0/(12.*3600.)
	END IF
c
	UM3=UM1+2.*DTT*(UMEAN-UM1)*TAUDAMP
	VM3=VM1+2.*DTT*(VMEAN-VM1)*TAUDAMP
	UM1=UM2+0.1*(UM1+UM3-2.*UM2)
	VM1=VM2+0.1*(VM1+VM3-2.*VM2)
	UM2=UM3
	VM2=VM3
c
c   Added 4/9/2007 to weight genesis points by variance of 850 hPa wind
c
c        VMINIT=(3.0+3.0*SQRT(UVAR850TRACK+VVAR850TRACK))/SCHI
C
c	ALAT=ALAT+DTT*VMEAN*DISFACTOR
c	ALAT=ALAT+DTT*VM2*DISFACTOR
	COSFAC=COS(ALAT*ADEG)
c	ALOG=ALOG+DTT*UMEAN*DISFACTOR/COSFAC
c	ALOG=ALOG+DTT*UM2*DISFACTOR/COSFAC
c
	WFAC=TT/(7.*24.*3600.)
	WFAC=MIN(WFAC,1.0)
	IF(ABS(UEC(KTIMEP)).LT.0.02.AND.
     1   ABS(VEC(KTIMEP)).LT.0.02)THEN
	  WFAC=1.0
	END IF
c
	UMEANW=(1.-WFAC)*UECMEAN+WFAC*UMEAN
	VMEANW=(1.-WFAC)*VECMEAN+WFAC*VMEAN
c
	ALAT=ALAT+DTT*VMEANW*DISFACTOR
	ALOG=ALOG+DTT*UMEANW*DISFACTOR/COSFAC
c
	dlogg=alog
c	IF(ALOG.GT.360.)ALOG=ALOG-360.
C
	RADIUS=((ALAT-CLAT1)**2+(COSFAC*(dlogg-CLOG1))**2)
	RMIN=MIN(RADIUS,RMIN)
C
	IF(JTIME.EQ.0.OR.PRINTIME.GE.DELPRINT)THEN
	 PRINTIME=0.0
	 JTIME=JTIME+1
	 MONTHPR(JTIME)=IMON
	 DAYPR(JTIME)=IDATE
	 HOURPR(JTIME)=IHOUR
	 LATPR(JTIME)=ALAT
	 LONGPR(JTIME)=ALOG
	 SHEARPR(JTIME)=VSHEAR
	END IF
C
	IF(TT.GT.ENDTIME.OR.ABS(ALAT).GT.LATMAXT.OR.ALOG.GT.LOGMAXT
     1    .OR.ALOG.LT.LOGMINT.OR.ABS(ALAT).LT.2.0)THEN
	   GOTO 900
	ELSE
	   GOTO 101
	END IF
  900 CONTINUE
C
	jint=1
	RMIN=0.001*SQRT(RMIN)/DISFACTOR
        M=M+1
c
 1000 CONTINUE
  980 FORMAT(2X,I2,2X,I2,2X,I2,3(2X,F7.3))
	CLOSE(11)
C
	UTE(1)=0.0
	DISTE(1)=0.0
	NST=0
	IMONI=MONTHPR(1)
	IDATEI=DAYPR(1)
	ITIMEI=HOURPR(1)
c
	DO 28 I=1,JTIME
	 IMON=MONTHPR(I)
	 IDATE=DAYPR(I)
	 ITIME=HOURPR(I)
	 ALAT=LATPR(I)
	 ALOGG=LONGPR(I)
	 VSH=SHEARPR(I)
c	 IF(ALOGG.LT.0.0)ALOGG=360.0+ALOGG
	 IF(I.NE.1.AND.ITIME.EQ.ITIMEOLD.AND.IDATE.EQ.IDATEOLD)GOTO 28
	 NST=NST+1 
	 IF(NST.GT.1)THEN
	  DELX=AFAC2*COS(ALAT*ADEG)*(ALOGG-ALOGGOLD)
	  DELY=AFAC2*(ALAT-ALATOLD)
	  DELTIME=3600.*FLOAT(ITIME-ITIMEOLD)
	  IF(IDATE.NE.IDATEOLD)DELTIME=DELTIME+3600*24
	  DELDIS=SQRT(DELX*DELX+DELY*DELY)
	  DISTE(NST)=DISTE(NST-1)+DELDIS
	  UTE(NST)=UTE(NST-1)+DELTIME
	 END IF
	 LATTE(NST)=ALAT
	 LOGGTE(NST)=ALOGG
	 MONTE(NST)=IMON
	 DATETE(NST)=IDATE
	 WRITETIME(NST)=ITIME
c
c  Increase shear for sensitivity experiment
c
c	 VSH=1.2*VSH
c
	 SHEARTE(NST)=VSH
	 ALATOLD=ALAT
	 ALOGGOLD=ALOGG
	 ITIMEOLD=ITIME
	 ITIMEB=ITIME
	 IMONOLD=IMON
	 IDATEOLD=IDATE
   28	CONTINUE
   29	CONTINUE
   27	FORMAT(12X,I2,I2,I2,17X,I3,A1,2X,I4,A1,2X,I3)
        NST0=NST
C
	DO I=1,NST
         UTES(I)=UTE(I)
	END DO
C
c	PRINT*, 'Program has read all files and is running'
C
	 WRITE(MONC,1801)MONTE(NST0)
	 WRITE(MOND,1801)MONTE(NST0)
	 WRITE(DATEC,1801)DATETE(NST0)
	 WRITE(DATETEMP,1805)WRITETIME(NST0)
 1801   FORMAT(I2.2)
 1805	FORMAT(I2.2,'Z')
C
C  ***  Interpolate positions to 2 hours intervals for accurate landfall  ***
C
c	NSTORM=(UTE(NST)+1.)/(2.*3600.)+1
        NSTORM=(UTE(NST))/(2.*3600.)
	DISTANCE(1)=0.0
	UTIME(1)=0.0
	JNST0=NSTORM
	NSTMAX=NST+1
c
	DO J=1,NSTORM
	 UTIME(J)=2.*3600.*FLOAT(J-1)
	 IF(UTIME(J).GE.UTE(NST0))JNST0=MIN(JNST0,J)
	 KI=NST
	 KM=NSTMAX
	 DO I=2,NST
	  IF(UTE(I).GT.(UTIME(J)+1.))KI=MIN(KI,I)
	 END DO
	 DO I=2,NSTMAX
	  IF(UTES(I).GT.(UTIME(J)+1.))KM=MIN(KM,I)
	 END DO
	 UTINV=1./(UTE(KI)-UTE(KI-1))
	 DISTANCE(J)=(DISTE(KI-1)*(UTE(KI)-UTIME(J))+DISTE(KI)*
     1   (UTIME(J)-UTE(KI-1)))*UTINV
         SHEARTRACK(J)=(SHEARTE(KI-1)*(UTE(KI)-UTIME(J))+SHEARTE(KI)*
     1   (UTIME(J)-UTE(KI-1)))*UTINV
	 ALAT=(LATTE(KI-1)*(UTE(KI)-UTIME(J))+LATTE(KI)*
     1   (UTIME(J)-UTE(KI-1)))*UTINV
	 ALOGG=(LOGGTE(KI-1)*(UTE(KI)-UTIME(J))+LOGGTE(KI)*
     1   (UTIME(J)-UTE(KI-1)))*UTINV
	 LATG(J)=ALAT
	 LOGG(J)=ALOGG
	 IMON=MONTE(KI-1)
	 IDATE=DATETE(KI-1)
C
	 ILAT=ALAT+90.501
	 ILOGG=ALOGG
	 IF(ILOGG.LT.1)ILOGG=360
	 ILATP=ILAT+1
	 ILOGGP=ILOGG+1
	 IF(ILOGGP.GT.360)ILOGGP=ILOGGP-360
C
	 IVLAT=1+(90.-ALAT)
	 IVLATP=IVLAT+1
	 IVLATP=MIN(IVLATP,181)
	 IVLOGG=ALOGG+1.
	 IVLOGGP=IVLOGG+1
	 IF(IVLOGGP.GT.360)IVLOGGP=IVLOGGP-360
C
         JLAT=1+(ALAT-LATSTART)/LATINC
         IF(ALOGG.GT.0.0)THEN
          JLOGG=1+(ALOGG-LONGSTART)/LONGINC
         ELSE
          JLOGG=1+(360.+ALOGG-LONGSTART)/LONGINC
         END IF
         JLATP=JLAT+1
         JLOGGP=JLOGG+1
         IF(JLOGGP.GT.NLO)JLOGGP=1
c
c  Define grid points further removed from track for temperature and humidity (added 2/7/2011)
c
	 JLATR=JLAT-1
	 JLATR=MAX(JLATR,1)
	 JLOGR=JLOGG-1
	 IF(JLOGR.LT.1)JLOGR=JLOGR+NLO
	 JLATRP=JLATP+1
	 JLATRP=MIN(JLATRP,NLA)
	 JLOGRP=JLOGGP+1
	 IF(JLOGRP.GT.NLO)JLOGRP=1
C
       KLAT=1+4.*(ALAT+90.)
       IF(ALOGG.GT.0.0)THEN
        KLOGG=1+4.*(ALOGG+0.01)
       ELSE
        KLOGG=1+4.*(360.01+ALOGG)
       END IF
       KLATP=KLAT+1
       KLOGGP=KLOGG+1
       IF(KLOGGP.GT.1440)KLOGGP=1
c
          IF(IDATE.GE.15)THEN
           IMONK=IMON
           IMONKP=IMON+1
           IF(IMONKP.EQ.13)IMONKP=1
           FAC=FLOAT(IDATE-15)/30.
          ELSE
           IMONK=IMON-1
           IF(IMONK.EQ.0)IMONK=12
           IMONKP=IMON
           FAC=FLOAT(15+IDATE)/30.
          END IF
c
	  DE1=(1.-FAC)*MIXDEPTH(IMONK,ILAT,ILOGG)+FAC*MIXDEPTH(IMONKP,
     1    ILAT,ILOGG)
	  DE2=(1.-FAC)*MIXDEPTH(IMONK,ILATP,ILOGG)+FAC*MIXDEPTH(IMONKP,
     1    ILATP,ILOGG)
	  DE3=(1.-FAC)*MIXDEPTH(IMONK,ILAT,ILOGGP)+FAC*MIXDEPTH(IMONKP,
     1    ILAT,ILOGGP)
	  DE4=(1.-FAC)*MIXDEPTH(IMONK,ILATP,ILOGGP)+FAC*MIXDEPTH(
     1    IMONKP,ILATP,ILOGGP)
	  DTH1=(1.-FAC)*THERMAL(IMONK,ILAT,ILOGG)+FAC*THERMAL(IMONKP,
     1    ILAT,ILOGG)
	  DTH2=(1.-FAC)*THERMAL(IMONK,ILATP,ILOGG)+FAC*THERMAL(IMONKP,
     1    ILATP,ILOGG)
	  DTH3=(1.-FAC)*THERMAL(IMONK,ILAT,ILOGGP)+FAC*THERMAL(IMONKP,
     1    ILAT,ILOGGP)
	  DTH4=(1.-FAC)*THERMAL(IMONK,ILATP,ILOGGP)+FAC*THERMAL(
     1    IMONKP,ILATP,ILOGGP)
c
	   JMOD=1+INT(J/6)
	   JMOD=MIN(JMOD,NN)
	   JMODP=JMOD+1
	   JMODP=MIN(JMODP,NN)
	   FAC=FLOAT(J)/6.0-INT(FLOAT(J)/6.0)
c
	   V1=(1.-FAC)*VPOTC(JLOGG,JLAT,JMOD)+FAC*
     1      VPOTC(JLOGG,JLAT,JMODP)
	   V2=(1.-FAC)*VPOTC(JLOGG,JLATP,JMOD)+FAC*
     1      VPOTC(JLOGG,JLATP,JMODP)
	   V3=(1.-FAC)*VPOTC(JLOGGP,JLAT,JMOD)+FAC*
     1      VPOTC(JLOGGP,JLAT,JMODP)
	   V4=(1.-FAC)*VPOTC(JLOGGP,JLATP,JMOD)+FAC*
     1      VPOTC(JLOGGP,JLATP,JMODP)
c
c  Changed 2/7/2011 to off-track values
c
	   RH1=(1.-FAC)*RH(JLOGR,JLATR,JMOD)+FAC*
     1      RH(JLOGR,JLATR,JMODP)    
	   RH2=(1.-FAC)*RH(JLOGR,JLATRP,JMOD)+FAC*
     1      RH(JLOGR,JLATRP,JMODP)
	   RH3=(1.-FAC)*RH(JLOGRP,JLATR,JMOD)+FAC*
     1      RH(JLOGRP,JLATR,JMODP)
	   RH4=(1.-FAC)*RH(JLOGRP,JLATRP,JMOD)+FAC*
     1      RH(JLOGRP,JLATRP,JMODP)
           T6001=(1.-FAC)*T600(JLOGR,JLATR,JMOD)+FAC*
     1      T600(JLOGR,JLATR,JMODP)
           T6002=(1.-FAC)*T600(JLOGR,JLATRP,JMOD)+FAC*
     1      T600(JLOGR,JLATRP,JMODP)
           T6003=(1.-FAC)*T600(JLOGRP,JLATR,JMOD)+FAC*
     1      T600(JLOGRP,JLATR,JMODP)
           T6004=(1.-FAC)*T600(JLOGRP,JLATRP,JMOD)+FAC*
     1      T600(JLOGRP,JLATRP,JMODP)
c
         VMIN=300.0
         IF(V1.GT.1.0)VMIN=MIN(VMIN,V1)
         IF(V2.GT.1.0)VMIN=MIN(VMIN,V2)
         IF(V3.GT.1.0)VMIN=MIN(VMIN,V3)
         IF(V4.GT.1.0)VMIN=MIN(VMIN,V4)
         IF(VMIN.LT.280.0)THEN
          V1=MAX(V1,VMIN)
          V2=MAX(V2,VMIN)
          V3=MAX(V3,VMIN)
          V4=MAX(V4,VMIN)
         END IF
         B1=BATH(KLOGG,KLAT)
         B2=BATH(KLOGG,KLATP)
         B3=BATH(KLOGGP,KLAT)
         B4=BATH(KLOGGP,KLATP)
c
c  ***   These 8 lines added 2/20/07 after land mask removed from calculating VPOT  ***
c  ***      Now VPOT uses skin temperature over land. These lines zero out VPOT     ***
c  ***                                over land                                     ***
c
c  ***   Forgot that bathymetry points here need to coincide with potential intensity points.  ***
c  ***                    More lines added 12/11/2008 to accomplish this                       ***
c
	 DUMLAT=LATSTART+LATINC*(JLAT-1)
	 DUMLONG=LONGSTART+LONGINC*(JLOGG-1)
         DUMLATP=LATSTART+LATINC*JLAT
         DUMLONGP=LONGSTART+LONGINC*JLOGG
	 KLATA=1+4*INT(DUMLAT+90.01)
	 KLOGGA=1+4*INT(DUMLONG+0.01)
         KLATAP=1+4*INT(DUMLATP+90.01)
         KLOGGAP=1+4*INT(DUMLONGP+0.01)
c
	 B1A=BATH(KLOGGA,KLATA)
         B2A=BATH(KLOGGA,KLATAP)
         B3A=BATH(KLOGGAP,KLATA)
         B4A=BATH(KLOGGAP,KLATAP)
c
         B1D=-MAX(MIN(B1A,0.0),-1.0)
         B2D=-MAX(MIN(B2A,0.0),-1.0)
         B3D=-MAX(MIN(B3A,0.0),-1.0)
         B4D=-MAX(MIN(B4A,0.0),-1.0)
c
	 V1=B1D*V1
	 V2=B2D*V2
	 V3=B3D*V3
	 V4=B4D*V4
c
	  ALATI=ALAT+90.
	  DELY=ALATI-INT(ALATI)
	  COSFAC=COS(ALAT*ADEG)
	  ALOGGI=ALOGG+180.
          delx=ALOGGI-int(ALOGGI)
          D1=(1.-DELX)*(1.-DELY)
          D2=DELY*(1.-DELX)
          D3=DELX*(1.-DELY)
          D4=DELX*DELY
          HMIX0(J)=D1*DE1+D2*DE2+D3*DE3+D4*DE4
          STRAT(J)=D1*DTH1+D2*DTH2+D3*DTH3+D4*DTH4
C
          DELY=1.+4.*(ALAT+90.)-FLOAT(KLAT)
          DELX=1.+4.*(360.+ALOGG)-INT(1.+4.*(360.+ALOGG))
          D1=(1.-DELX)*(1.-DELY)
          D2=DELY*(1.-DELX)
          D3=DELX*(1.-DELY)
          D4=DELX*DELY
          BATHTRACK(J)=D1*B1+D2*B2+D3*B3+D4*B4
c
          B1M=MAX(0.0,B1)
          B2M=MAX(0.0,B2)
          B3M=MAX(0.0,B3)
          B4M=MAX(0.0,B4)
          BMEAN=0.25*(B1M+B2M+B3M+B4M)
	  HMIX0(J)=MIN(HMIX0(J),-BATHTRACK(J))
	  HMIX0(J)=MAX(HMIX0(J),20.0)
           DELY=1.+(ALAT-LATSTART)/LATINC-FLOAT(JLAT)
           DELX=1.+(360.+ALOGG-LONGSTART)/LONGINC-
     1      INT(1.+(360.+ALOGG-LONGSTART)/LONGINC)
           D1=(1.-DELX)*(1.-DELY)
           D2=DELY*(1.-DELX)
           D3=DELX*(1.-DELY)
           D4=DELX*DELY
          VPTRACK(J)=D1*V1+D2*V2+D3*V3+D4*V4
c
c   Changed weights to off-track values, 2/7/2011
c
           DELY=0.25+0.5*DELY
           DELX=0.25+0.5*DELX
           D1=(1.-DELX)*(1.-DELY)
           D2=DELY*(1.-DELX)
           D3=DELX*(1.-DELY)
           D4=DELX*DELY
c
          RHTRACK(J)=D1*RH1+D2*RH2+D3*RH3+D4*RH4
          T600TRACK(J)=D1*T6001+D2*T6002+D3*T6003+D4*T6004
c
	  IF(V1.LT.1.0.OR.V2.LT.1.0.OR.V3.LT.1.0.OR.V4.LT.1.0)THEN
	   VPTRACK(J)=MAX(V1,V2,V3,V4)
	  END IF
	  UTEMP=UTIME(J)/(24.*3600.)
	END DO
C
        IF(gmeth.ne.'clim')THEN
C
	IF(BATHTRACK(1).GE.0.0.OR.VPTRACK(1).LT.40.0)THEN
c	IF(BATHTRACK(1).GE.0.0.OR.VPTRACK(1).LT.40.0.OR.
c     1    SHEARTRACK(1).GE.14.0.OR.RHTRACK(1).LE.40)THEN
c     1    SHEARTRACK(1).GE.10.0)THEN
c         OPEN(UNIT=9,FILE='flag3',STATUS='UNKNOWN')
c         WRITE(9,*)IFLAG
c         CLOSE(9)
	 GOTO 5
	END IF
C
	END IF
C
C      *** Non-dimensionalize Parameters and Initial Conditions  ***
c
	TIME=UTIME(NSTORM)/ATIME
C
	UT=(DISTANCE(2)-DISTANCE(1))/(UTIME(2)-UTIME(1))
	UTINIT=UT
c
	call random_number(arand)
	VINTEN0=(0.0+2.0*(arand-0.5))/(6.*3600.)
c
        VINTEN=2.*(VINTEN0+VINOFFSET)*ATIME/(VPTRACK(1))
        VINTEN=VINTEN*2.*SCHI/(MAX(VPTRACK(1),1.0))
c        DWELL=0.8+1.5*VINTEN
c
c   Changed DWELL to 0.5  2/7/2011
c
	call random_number(arand)
	DWELL=0.85+0.3*(arand-0.5)
	DWELL=1.-DWELL
	DWELL=MAX(DWELL,0.0)
	DWELL=MIN(DWELL,1.0)
c
	DO 33 I=1,NSTORM
	 DISTANCE(I)=DISTANCE(I)/ALENGTH
	 UTIME(I)=UTIME(I)/ATIME
	 HMIX0(I)=HMIX0(I)*AMIXFAC
	 BATHTRACK(I)=BATHTRACK(I)*AMIXFAC
	 STRAT(I)=STRAT(I)*0.01/GM
	 SHEARTRACK(I)=SHEARTRACK(I)*ATIME/ALENGTH
   33	CONTINUE
C
C            ****   INITIALIZE FIELDS   ***
c
c  ***  Make random draw from a normal distribution
c
c 9056       av1=2.*rand(0)-1.
c            av2=2.*rand(0)-1.
 9056	    call random_number(arand)
            av1=2.*arand-1.
	    call random_number(arand)
            av2=2.*arand-1.
            avr=av1**2+av2**2
            if(avr.ge.1.0.or.avr.eq.0.)goto 9056
	    avr=max(avr,1.0e-8)
            avfac=sqrt(-2.*log(avr)/avr)
            Gasdev=av2*avfac
c
	VMINIT=VMINIT0+5.0*AKNOTFAC*Gasdev-0.7*UTINIT/SCHI
c	VMINIT=VMINIT0
	VMINIT=MAX(VMINIT,7.0/SCHI)
	RTEMP=SQRT(RMINIT*RMINIT+2.*RMINIT*VMINIT)
	IMAX=RTEMP/DR
	RTEMP=DR*FLOAT(IMAX)
	RM=-VMINIT+SQRT(VMINIT*VMINIT+RTEMP*RTEMP)
c
	DO 60 I=2,NR
	 R=FLOAT(I-1)*DR
	 IF(R.GT.R0)GOTO 40
	 IF(I.GT.IMAX)GOTO 30
	 RBS1(I)=R*R/(1.+2.*VMINIT/RM)
	 RBS2(I)=RBS1(I)
	 GOTO 50
   30  RBS1(I)=(0.5*R*R-VMINIT*RM/(1.-(RM/R0)**2))/
     1   (0.5-VMINIT*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)
	 RTEMP=SQRT(RBS2(I))
	  KH=NSTORM
	  DO 55 J=2,NSTORM
	   IF(DISTANCE(J).GT.RTEMP)KH=MIN(KH,J)
   55	  CONTINUE
	  HMIX(I)=HMIX0(NSTORM)
	  IF(KH.LT.NSTORM)THEN
	   HMIX(I)=(HMIX0(KH)*(RTEMP-DISTANCE(KH-1))+HMIX0(KH-1)*
     1      (DISTANCE(KH)-RTEMP))/(DISTANCE(KH)-DISTANCE(KH-1))
	  END IF
          UHMIX1(I)=SQRT(HMIX(I)*HMIX(I)*HMIX(I)*DSST)/AMIX
	  UHMIX2(I)=UHMIX1(I)
   60 CONTINUE
	XS1(2)=0.0
	HMIX(1)=HMIX(2)
	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
	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)
	 X1(I)=XS1(I)-0.05
	 X2(I)=X1(I)
	 XS2(I)=XS1(I)
	 XM1(I)=XM0
       IF(RB1(I).LT.1.3*RM)XM1(I)=(1.-DWELL)*XS1(I)+DWELL*XM0
	 XM2(I)=XM1(I)
	 XSURST(I)=XS1(I)
	 XVIS(I)=0.0
	 XMVIS(I)=0.0
	 MU1(I)=0.0
	 MU2(I)=0.0
   75 CONTINUE
	PS2(1)=0.0
	PS3(1)=0.0
	PS0(1)=0.0
	VIS(1)=0.0
	MU2(1)=0.0
	GB(1)=0.0
	RMM2(1)=0.0
	XSURST(1)=XS1(1)
	VEXT1=0.0
	VEXT2=VEXT1
	TT=0.0D0
	TT1=0.0D0
	NTT=0
C
C             ***  SET TIME LOOPING PARAMETERS AND BEGIN TIME LOOP  ***
C
	ADIST=0.0
	DT=DT0
	DT1=DT
	TIMEPRINT=12.0
	NEXT=0
	NTCOUNT=0
c	VTHRESH=13./SCHI
c
c  Switched to 8.0 on 5/31/05
c
        VTHRESH=15./SCHI
	VTHRESH2=0.5
	TTHRESH=2.*24.*3600./ATIME
	VMAX3=2.*VTHRESH
	VMAX2=VMAX3
C
C           ***  PROGRAM RETURNS TO 77 AFTER EACH TIME STEP  ***
C
   77 CONTINUE
C
	TT=TT+DT
c
c   Switched from VMAX3 to VMAX2 in next statement, 5/20/05
c
	IF(TT.GT.TIME.OR.(VMAX2.LT.VTHRESH.AND.TT.GT.TTHRESH))GOTO 705
C
C
C            *** SET BOUNDARY VALUES FOR AND ADVECTION TERMS ***
C
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)
	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)
C
C  **  Calculate certain along-track quantities
C
	KU=NSTORM
	DO 80 J=2,NSTORM
	 IF(UTIME(J).GT.TT)KU=MIN(KU,J) 
   80	CONTINUE
	ADIST=DISTANCE(NSTORM)
	VPOTM=VPTRACK(NSTORM)
	RHM=RHTRACK(NSTORM)
	T600M=T600TRACK(NSTORM)
	VSHEAR=SHEARTRACK(NSTORM)
	BATHY=BATHTRACK(NSTORM)
        BATHV=BATHVAR(NSTORM)
	ACCEL=0.0
	DUTINV=1./(UTIME(KU)-UTIME(KU-1))
	IF(KU.LE.NSTORM)THEN
	 ADIST=(DISTANCE(KU)*(TT-UTIME(KU-1))+DISTANCE(KU-1)*(UTIME(KU)-
     1    TT))*DUTINV
         VPOTM=(VPTRACK(KU)*(TT-UTIME(KU-1))+VPTRACK(KU-1)*(UTIME(KU)-
     1    TT))*DUTINV
         RHM=(RHTRACK(KU)*(TT-UTIME(KU-1))+RHTRACK(KU-1)*(UTIME(KU)-
     1    TT))*DUTINV
         T600M=(T600TRACK(KU)*(TT-UTIME(KU-1))+T600TRACK(KU-1)*
     1    (UTIME(KU)-TT))*DUTINV
         VSHEAR=(SHEARTRACK(KU)*(TT-UTIME(KU-1))+SHEARTRACK(KU-1)*
     1    (UTIME(KU)-TT))*DUTINV
         BATHY=(BATHTRACK(KU)*(TT-UTIME(KU-1))+BATHTRACK(KU-1)*
     1    (UTIME(KU)-TT))/(UTIME(KU)-UTIME(KU-1))
         BATHV=(BATHVAR(KU)*(TT-UTIME(KU-1))+BATHVAR(KU-1)*
     1    (UTIME(KU)-TT))/(UTIME(KU)-UTIME(KU-1))
         IF(KU.GT.2.AND.KU.LE.(NSTORM-1))THEN
          ACCEL=((DISTANCE(KU+1)-DISTANCE(KU-1))/(UTIME(KU+1)-
     1     UTIME(KU-1))-(DISTANCE(KU)-DISTANCE(KU-2))/(UTIME(KU)-
     2     UTIME(KU-2)))*DUTINV
         END IF
	END IF
c       CDBATH=1.0+MAX(0.0,BATHV)/120.
	CDBATH=1.0
c
	RHM=MIN(RHM,90.0)
	RHM=MAX(RHM,1.0)
	TM=T600M-273.15
	TMA=T600M
	ESM=6.112*EXP(17.67*TM/(243.5+TM))
	QSM=0.622*ESM/600.0
	SMEX=2.5E6*QSM*(0.01*RHM-1.)-461.0*TMA*QSM*0.01*RHM*
     1    LOG(0.01*RHM)
	XM0=(TS-TO)*SMEX/(CHI*TMA)
	XM1(NR)=XM0
	XM1(NR-1)=XM0
	XM2(NR)=XM0
	XM2(NR-1)=XM0
	IF(TT.LT.(1.5*DT))THEN
	 DO I=1,NR-2
	  XM1(I)=XM0
          IF(RB1(I).LT.1.6*RM)XM1(I)=(1.-DWELL)*XS1(I)+DWELL*XM0
	  XM2(I)=XM1(I)
	 END DO
	END IF
c
        GRATB1=(1.+2.5E6*QSM/(287.*TMA))/
     1   (1.+2.5E6*2.5E6*QSM/(461.*1005.*TMA*TMA))
	Q=2.0*2.5E6*QSM*EF/(CHI*GRATB1)
	BFAC=1.-DR*SQRT(2./Q)
C
	ACCEL=0.0
      UT0=(DISTANCE(KU)-DISTANCE(KU-1))*DUTINV
      UT=UT0+ACCEL*(TT-0.5*(UTIME(KU)+UTIME(KU-1)))
C
	VMAX=0.0
	AMMAX=0.0
	DO I=2,NR-1
	 R=(FLOAT(I-1))*DR
	 V=0.5*(R*R-RBS2(I))/RB2(I)
	 VMAX=MAX(VMAX,V)
	 AMMAX=MAX(AMMAX,MU2(I))
	END DO
	AMMAX=MAX(AMMAX,20.0)
	VMAX2=VMAX
	VMAX4=MAX(VMAX,0.4)
	VMAX=VMAX+TRANSFAC*UT*ALENGTH/(SCHI*ATIME)
	VMAX3=VMAX
	wf=tt*atime/(12.*3600.)
	wf=min(wf,1.0)
	vmtarget=(1.-wf)*vminit0+wf*vminitnew
	VMINITR=vmtarget-0.7*UTINIT/SCHI
C
c  Changed coefficient in line below from 0.007 to 0.005 on 05/21/05
c    Changed to 0.0045 on 06/01/05
c
c  The following lines were modified on 9/28/2010 to flatten initial rate of change of intensity
c
	IF((TT*ATIME).LE.(12.*3600.))THEN
c
c  Changed VEXT2 from 30 to 15, 2/7/2011
c
c	 VEXT2=0.01*(VMAX-VMINITR)*
c     1      (abs(VMAX-VMINITR))**2
c	 VEXT2=0.03*(VMAX-VMINITR)
	 VEXT2=0.001*(VMAX-VMINITR)
	ELSE
c
c   Changed coefficient from 0.02 to 0.03  2/6/2011
c
         VEXT2=0.06*0.0045*VMAX4*VMAX4*VSHEAR*VSHEAR*AMMAX
	END IF
	VEXT1=VEXT2
        VRAT=VPOTM/VMTH
	VRAT2=VRAT*VRAT
	DISR0=ADIST+R0
	KHR0=NSTORM
	DO 81 J=2,NSTORM
	 IF(DISTANCE(J).GT.DISR0)KHR0=MIN(KHR0,J)
   81   CONTINUE
	TFAC=HMIX0(NSTORM)
	IF(KHR0.LT.NSTORM)THEN
	 TFAC=(HMIX0(KHR0)*(DISR0-DISTANCE(KHR0-1))+HMIX0(KHR0-1)*
     1    (DISTANCE(KHR0)-DISR0))/(DISTANCE(KHR0)-DISTANCE(KHR0-1))
	END IF
      UHMIX1(NR)=SQRT(TFAC*TFAC*TFAC*DSST)/AMIX
	UHMIX2(NR)=UHMIX1(NR)
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)
	 VABS=0.5*ABS(R*R-RBS1(I))/RB1(I)
	 VABS=MAX(VABS,1.0E-9)
	 UR=UT/(VABS*SCHI)
       PS0(I)=PS0(I)*CDBATH*(1.+MIN(CDV*VABS,CDMAX))
c
c added 5/24/05
c
c	PS0(I)=PS0(I)*(1.0+1.5*UR*UR)**(2./3.)
C
C   ***   Reduce drag coefficient under stable conditions over water ***
C
	IF(BATHY.LT.0.0.AND.XS1(I).GT.XSURST(I))THEN
	 PS0(I)=PS0(I)/(1.+2000.*(XS1(I)-XSURST(I)))
	END IF
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)
c	  C2B=RMS2(I)/(RMM2(I+1)-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)
c	  C2B=RMS2(I)*DRI
         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)
c	  C1B=RMS2(I-1)*DRI
         END IF
	 C1=C1B*4.*XMDIF*PNU
	 C2=C2B*4.*XMDIF*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))
	  A1=PS3(I)-PS2(I)
	  PMAX=MAX(PMAX,ABS(A1))
  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
	VMTEMP=0.0
	DO 99 I=2,NR-1
         R=(FLOAT(I-1))*DR
	 VMTEM1=(R*R-RBS2(I))/RB2(I)
	 IF(VMTEM1.GT.VMTEMP)THEN
	  VMTEMP=VMTEM1
	  IVMAX=I
	 END IF
   99	CONTINUE
C
	KH1=NSTORM
	DO 991 J=2,NSTORM
	 IF(DISTANCE(J).GT.ADIST)KH1=MIN(KH1,J)
  991	CONTINUE
	AHMLOCAL2=HMIX0(NSTORM)
	SLOCAL=STRAT(NSTORM)
	IF(KH1.LT.NSTORM)THEN
	 AHMLOCAL2=(HMIX0(KH1)*(ADIST-DISTANCE(KH1-1))+HMIX0(KH1-1)*
     1    (DISTANCE(KH1)-ADIST))/(DISTANCE(KH1)-DISTANCE(KH1-1))
	 SLOCAL=(STRAT(KH1)*(ADIST-DISTANCE(KH1-1))+STRAT(KH1-1)*
     1    (DISTANCE(KH1)-ADIST))/(DISTANCE(KH1)-DISTANCE(KH1-1))
	END IF
	BMIX2=HMIX(2)
	BMIX2=MAX(BMIX2,0.0001,AHMLOCAL2)
	XSURM2=-SLOCAL*((BMIX2-AHMLOCAL2)**2)/BMIX2
      XSURM2=XSURM2+DSST*(AHMLOCAL2-BMIX2)/BMIX2
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
	W0=(PS0(I)-PS0(I-1)+GAMMA*(VIS(I)-VIS(I-1)))/(RBS2(I)-RBS2(I-1))
c	W0=(PS0(I)-PS0(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)
	VABS=MAX(VABS,1.0E-9)
C
C   ***  Calculate forcing of ocean temperature  ***
C
	SSTR(I)=0.0
	SEAF=0.0
      KHI=NSTORM
	DO 251 J=2,NSTORM
	 IF(DISTANCE(J).GT.(ADIST+RB2(I)))KHI=MIN(KHI,J)
  251	CONTINUE
      IF(OM.EQ.'Y'.OR.OM.EQ.'y')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)
	 CVABS=VABS*(1.+MIN(CDV*VABS,CDMAX))
c	 CVABS=VABS*(1.+CDM*VABS/VCD)/(1.+VABS/VCD)
	 SEAF=ADVOM+CVABS*VABS
	 AHMLOCAL1=HMIX0(NSTORM)
	 SLOCAL2=STRAT(NSTORM)
	 IF(KHI.LT.NSTORM)THEN
	  AHMLOCAL1=(HMIX0(KHI)*(ADIST+RB2(I)-DISTANCE(KHI-1))+
     1     HMIX0(KHI-1)*(DISTANCE(KHI)-ADIST-RB2(I)))/
     2     (DISTANCE(KHI)-DISTANCE(KHI-1))
	  SLOCAL2=(STRAT(KHI)*(ADIST+RB2(I)-DISTANCE(KHI-1))+
     1     STRAT(KHI-1)*(DISTANCE(KHI)-ADIST-RB2(I)))/
     2     (DISTANCE(KHI)-DISTANCE(KHI-1))
	 END IF
	 BMIXI=HMIX(I)
	 BMIXI=MAX(BMIXI,0.0001,AHMLOCAL1)
	 XSURM1=-SLOCAL2*((BMIXI-AHMLOCAL1)**2)/BMIXI
       XSURM1=XSURM1+DSST*(AHMLOCAL1-BMIXI)/BMIXI
	 SSTFAC=2.*RB2(IVMAX)/RB2(I)
	 SSTFAC=MIN(SSTFAC,1.0)
	 SSTR(I)=SSTFAC*XSURM2+(1.-SSTFAC)*XSURM1
	END IF
C
	CECD1=CECD
	IF(BATHY.GT.0.0)THEN
        CECD1=CECD*(1.-BATHY/(4.*AMIXFAC))
c	 CECD1=0.0
	 CECD1=MAX(CECD1,0.0)
	 VEXTF=0.0
	END IF
C
c	IF(TT.GT.TTHRESH.AND.VMAX2.LE.VTHRESH2)THEN
c	 CECD1=CECD1*VMAX2/VTHRESH2
c	END IF
C
	XSUR=1.+SSTR(I)-EF*P(I)+(EXP(-BETA*P(I))-1.)/(1.-H_A)
	XSUR=VRAT2*(XSUR+EF*P(I))-EF*P(I)
	XSURST(I)=XSUR
	CVABS=VABS*(1.+MIN(CDV*VABS,CDMAX))
	CKABS=VABS*(1.+MIN(CDV*VABS,CKMAX))
c	IF((SCHI*VABS).GT.33.0)CKABS=CKABS*1.13
	FLUX=CKABS*(XSUR-X1(I))
	FLUX=FLUX*CECD1+EF*CVABS*CDBATH*VABS*VABS
C
	RADB=RAD
	IF((RADB*XS1(I)).GT.RADMAX)RADB=RADMAX/
     1	(XS1(I))
      IF(W0D.LE.0.0)THEN
       ACHIM=MAX(0.001,(X1(I)-XM1(I)))
       MEQ=W0+(ADV+FLUX-GAMMA*(VADV+XVIS(I)-RADB*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)-RADB*XS1(I)))/
     1  ACHIM
      END IF
C
C  ***  Limit maximum value of M   ***
C
	MEQ=MIN(MEQ,1200.0)
c	MEQ=MIN(MEQ,1000.0)
C
C  ***    If conditionally stable, set MEQ to 0    ***
C
	IF(X2(I).LT.(XS2(I)-0.0001))MEQ=0.0
c
	IF(VMAX2.GT.(1.04*VPOTM/SCHI))MEQ=0.0
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)
      XMF=EFRAC*(MU2(I)+WXM)*(X1(I)-XM1(I))
      XMF=XMF+XMVIS(I)-RADB*XS1(I)*(0.75*GRATB+0.25)
      VEXTFAC=1.0
      IF(VEXT2.GE.0.0)THEN
       XMF=XMF-500.*VEXT2*VEXTFAC*(XM1(I)-XM0)
	ELSE
	 IF(I.LE.IVMAX)THEN
	  AMF=-500.*VEXT2*(XS1(I)-XM1(I))
c	  AMF=-1500.*VEXT2*(XS1(I)-XM1(I))
	  AMF=MIN(AMF,(-100.*VEXT2))
	  XMF=XMF+AMF
	 END IF
	END IF
	AMC=MIN(0.0,W0D)
	XF=ADV+FLUX+(AMD+AMC)*(X1(I)-XM1(I))
	XF=XF*DPR
	XSF=VADV+XVIS(I)-RADB*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),X3(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)
C
C  ***    Bail out if r**2 is negative   ***
C
	IF(RBS3(I).LT.0.0)THEN
c	 PRINT*, 'Program failed'
	 goto 5
	END IF
  100 CONTINUE
C
C  ***   Rectify saturation entropy to be consistent with thermal wind
C
c	DO I=NR-1,2,-1
c	 R=FLOAT(I-1)*DR
c	 XS3(I)=XS3(I+1)+0.5*R*R*R*DR*(1./RBS3(I)-1./RTS3(I))
c	END DO
c	XS3(1)=XS3(2)
C
	SSTR(NR)=SSTR(NR-1)
	SSTR(1)=SSTR(2)
C
C  ***  Print progress report to screen  ***
C
	DIMTIME=TT*ATIME/(3600.)
	IF(DIMTIME.GT.TIMEPRINT)THEN
c	 PRINT*, ' Finished  ',TIMEPRINT,'  hours'
	 TIMEPRINT=TIMEPRINT+12.0
	END IF
C
C  Interpolate latitude and longitude to current time
C
        ALATS=LATG(NSTORM)
        ALOGS=LOGG(NSTORM)
        IF(KU.LE.NSTORM)THEN
         ALATS=(LATG(KU)*(TT-UTIME(KU-1))+LATG(KU-1)*(UTIME(KU)-
     1    TT))*DUTINV
         ALOGS=(LOGG(KU)*(TT-UTIME(KU-1))+LOGG(KU-1)*(UTIME(KU)-
     1    TT))*DUTINV
        END IF
c
c   Find maximum wind within city radius
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(VMAX,V)
        END DO
        VMTEMP1=VMAX*SCHI+0.7*UT*ALENGTH/ATIME
        VMTEMP2=VMAX*SCHI*3600./1852.
        VMTEMP1=VMTEMP1*3600./1852.
c
        DELX=AFAC2*COS(ALATS*ADEG)*(ALOGS-CLOG1)
        DELY=AFAC2*(ALATS-CLAT1)
        DISTA=0.001*SQRT(DELX*DELX+DELY*DELY)
c
c    Statment below added 8/1/2009
c
c	IF(VMTEMP1.GE.VMCRIT.AND.VMTEMP2.GE.25.0)IFLAGB=1
c
        IF((DISTA.LE.CRAD.or.jint.eq.-1).AND.
     1    VMTEMP1.GE.VMCRIT.AND.VMTEMP2.GE.25.0)THEN
         IFLAG=1
        END IF
C
C          ***  STORE CERTAIN GRAPHICS ARRAYS EVERY DTG TIME UNITS  ***
C
	TT1=TT1+DT
        IF(ABS(TT1-DTG).GT.(0.5*DT).AND.TT.GT.(1.5*DT))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
	VMAX=0.0
	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
	 TIMEG(NTT)=TT*ATIME
	 TIMEG(NTT)=TIMEG(NTT)/(24.*3600.)
	 LATTIME(NTT)=ALATS
	 LOGTIME(NTT)=ALOGS
	 IF(LOGTIME(NTT).LT.0.0)LOGTIME(NTT)=LOGTIME(NTT)+360.0
	 IF(LOGTIME(NTT).GT.360.0)LOGTIME(NTT)=LOGTIME(NTT)-360.0
c
	 UTRANS(NTT)=AFAC2*COS(ADEG*LATG(KU))*(LOGG(KU)-
     1     LOGG(KU-1))*DUTINV/ATIME
	 VTRANS(NTT)=AFAC2*(LATG(KU)-LATG(KU-1))*DUTINV/ATIME
C
	 VMG(NTT)=VMAX*SCHI
	 VMG(NTT)=3600.*VMG(NTT)/1852.
	 VMGP=VMAX*SCHI
C
	 VSHG(NTT)=VSHEAR*ALENGTH/ATIME
	 VPOTG(NTT)=VPOTM*3600./1852.
c
c  Factor of 2 added 1/6/2008 to be consistent with previous implementations of 'citystore'
c  Factor of 2 replaced with 1.9 on 4/4/2009 based on comparison with DeMaria and Kossin data
c
	 RMG(NTT)=1.9*RMAX*SCHI*0.001/FC
	 PCG(NTT,1)=PA*EXP(BETA*PMIN)
	 PEXP=0.5
	 PCG(NTT,2)=PA*EXP(BETA*PRMAX-VMGP**2/
     1	 (2.*PEXP*287.*TSA))
	 MUG(NTT)=MUMAX*CD*SCHI
	 MDG(NTT)=MDMIN*CD*SCHI
	 HMG(NTT)=-HMIX(2)/AMIXFAC
	 HMG0(NTT)=-AHMLOCAL2/AMIXFAC
	 SSTT(NTT)=SSTR(2)*0.5*GM/AMIXFAC
  445  CONTINUE
C
C         *** IF VMAX EXCEEDS CERTAIN VALUES, CUT DOWN TIME STEP ***
C
	IF(VMAX.GT.1.5)THEN
	   DT=0.5*DT1
	END IF
	IF(VMAX.GT.2.0)THEN
           DT=0.2*DT1
	END IF
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))
	 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)
	 UHMIX2(I)=UHMIX3(I)
	 ADIST1=ADIST+UT*DT
         IF(OM.EQ.'Y'.OR.OM.EQ.'y')THEN
          KHI=NSTORM
          DO 501 J=2,NSTORM
           IF(DISTANCE(J).GT.(ADIST1+RB2(I)))KHI=MIN(KHI,J)
  501     CONTINUE
          AHM0=HMIX0(NSTORM)
	  SLOCAL=STRAT(NSTORM)
	  BATHY=BATHTRACK(NSTORM)
          IF(KHI.LT.NSTORM)THEN
           AHM0=(HMIX0(KHI)*(ADIST1+RB2(I)-DISTANCE(KHI-1))+
     1      HMIX0(KHI-1)*(DISTANCE(KHI)-ADIST1-RB2(I)))/
     2      (DISTANCE(KHI)-DISTANCE(KHI-1))
           SLOCAL=(STRAT(KHI)*(ADIST1+RB2(I)-DISTANCE(KHI-1))+
     1      STRAT(KHI-1)*(DISTANCE(KHI)-ADIST1-RB2(I)))/
     2      (DISTANCE(KHI)-DISTANCE(KHI-1))
	   BATHY=(BATHTRACK(KHI)*(ADIST1+RB2(I)-DISTANCE(KHI-1))+
     1      BATHTRACK(KHI-1)*(DISTANCE(KHI)-ADIST1-RB2(I)))/
     2      (DISTANCE(KHI)-DISTANCE(KHI-1))
          END IF
          AHM02=AHM0*(AHM0-DSST/MAX(0.1,SLOCAL))	
	  TEMP=0.25*AHM02*AHM02+AMIX*AMIX*
     1     UHMIX2(I)*UHMIX2(I)/MAX(0.1,SLOCAL)
	  HMIX(I)=SQRT(ABS(0.5*AHM02+SQRT(ABS(TEMP))))
	  HMIX(I)=MIN(HMIX(I),-BATHY)
	  HMIX(I)=MAX(HMIX(I),AHM0)
	 END IF
  600 CONTINUE
C
  700 GOTO 77
  705 CONTINUE
c
	IF(IFLAG.EQ.0)THEN
c	 iseed=time8 ()
c	 gdum=rand(1)
c	 call srand(iseed)
c
c  Statements added 8/1/2009
c
c	 IF(IFLAGB.EQ.0.AND.ICOUNT.LE.20)THEN
c           BACKSPACE(18)
c	   ICOUNT=ICOUNT+1
c	 ELSE
c	   ICOUNT=0
c	 END IF
         GOTO 7
	END IF
C
C   ***     Write to output files     ***
C
	IMON=IMONI
	IDATE=IDATEI
	ITIME=ITIMEI
c
	NUMSTORMS=NUMSTORMS+1
c
	if(numstorms.lt.10)then
	 write(ch1,1201)numstorms
         DIRECTORY='ascii/'//TRIM(CITYNAME)//'/hurr'//ch1//'.out'
c
	elseif(numstorms.ge.10.and.numstorms.lt.100)then
	 write(ch2,1202)numstorms
         DIRECTORY='ascii/'//TRIM(CITYNAME)//'/hurr'//ch2//'.out'
c
        elseif(numstorms.ge.100.and.numstorms.lt.1000)then
         write(ch3,1203)numstorms
         DIRECTORY='ascii/'//TRIM(CITYNAME)//'/hurr'//ch3//'.out'
c
        elseif(numstorms.ge.1000.and.numstorms.lt.10000)then
         write(ch4,1204)numstorms
         DIRECTORY='ascii/'//TRIM(CITYNAME)//'/hurr'//ch4//'.out'
c
        elseif(numstorms.ge.10000.and.numstorms.lt.100000)then
         write(ch5,1205)numstorms
         DIRECTORY='ascii/'//TRIM(CITYNAME)//'/hurr'//ch5//'.out'
c
	end if
c
 1201   format(I1.1)
 1202   format(I2.2)
 1203   format(I3.3)
 1204   format(I4.4)
 1205   format(I5.5)
c
	OPEN(UNIT=12,FILE=TRIM(DIRECTORY),STATUS='UNKNOWN')
C
        DO I=1,NTT
C
	  IF(I.GT.1)THEN
	   IDATE=IDATEI+TIMEG(I)+FLOAT(ITIMEI)/24.
	   ITIME=ITIMEI+24*(TIMEG(I)-INT(TIMEG(I)))
	   IF(ITIME.GE.24)ITIME=ITIME-24
	   IMON=IMONI
	   IMOB=IMON
	   IDATB=IDATE
	   IF((IMOB.EQ.1.OR.IMOB.EQ.3.OR.IMOB.EQ.5.OR.IMOB.EQ.7.OR.
     1     IMOB.EQ.8.OR.IMOB.EQ.10.OR.IMOB.EQ.12).AND.IDATB.GT.31)THEN
		  IDATE=IDATE-31
	      IMON=IMON+1
	    END IF
	    IF((IMOB.EQ.4.OR.IMOB.EQ.6.OR.IMOB.EQ.9.OR.IMOB.EQ.11)
     1		.AND.IDATB.GT.30)THEN
		  IDATE=IDATE-30
	      IMON=IMON+1
	    END IF
	    IF(IMOB.EQ.2.AND.IDATB.GT.28)THEN
	      IDATE=IDATE-28
	      IMON=IMON+1
	    END IF
	    IF(IMON.GT.12)IMON=1
         END IF
C
	  WRITE(12,810)IMON,IDATE,ITIME,LATTIME(I),LOGTIME(I),
     1    VMG(I),PCG(I,2),RMG(I),VSHG(I),VPOTG(I)
c
	  END DO
c
	CLOSE(12)
C
  710 CONTINUE
  800 CONTINUE
  810 FORMAT(3(1X,I2),2(1X,F8.3),5(1X,F7.2))
c
        if(gmeth.ne.'clim')then
c
         SFREQ=globalfac*areafac*float(numstorms)/float(ntotal)
c
        else
c
	 SFREQ=1.1*climfac*float(numstorms)/float(ntotal)
c
        end if
c
	 OPEN(UNIT=12,FILE='ascii/'//TRIM(CITYNAME)//'/stats.out',
     1    status='unknown')
c	 write(12,1222)BAS,' ',gmeth,' ',' ',TRIM(CITYNAME),CLAT,
c     1    CLOG,CRAD,numstorms,ntotal,SFREQ
	 write(12,1222)BAS,gmeth,TRIM(CITYNAME),CLAT,
     1    CLOG,CRAD,numstorms,ntotal,SFREQ
	 CLOSE(12)
 1222    FORMAT(1X,A2,1X,A4,1X,A40,1X,F9.4,1X,F9.4,1X,F9.1,1X,
     1    I5,1X,I9,1X,F12.6)
C
	IF(NUMSTORMS.GE.NUMTRACKS)THEN
c
         OPEN(UNIT=12,FILE='ascii/'//TRIM(CITYNAME)//'/stats.txt',
     1    status='unknown')
	 write(12,'(1X,A35,T45,A)')'Basin:',BAS
	 write(12,'(1X,A35,T45,A)')'Filename:',TRIM(CITYNAME)
	 write(12,'(1X,A35,T45,A)')'Seeding method (random or climo):'
     1     , gmeth
	 write(12,'(1X,a35,T41,F9.4)')'POI Latitude:',CLAT
         write(12,'(1X,a35,T41,F9.4)')'POI Longitude:',CLOG
	 write(12,'(1X,A35,T42,F8.1)')'Search radius (km):',CRAD
         write(12,'(1X,A35,T45,F5.1)')'Critical wind speed (kts):',
     1     VMCRIT
	 write(12,'(1X,A35,T44,I6)')'Number of tracks:',numstorms
	 write(12,'(1X,A35,T44,I8)')'Number of unfiltered tracks:',
     1     ntotal
	 write(12,'(1X,A35,T44,F8.2)')'Annual frequency:',sfreq
         CLOSE(12)
c
	 STOP
c
	ELSE
	 GOTO 5
	END IF
c
	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.5
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 choldc(a,n,np,p)
	INTEGER n,np
	REAL a(np,np),p(n)
c
c	Given a positive-definite symmetric matrix a(1:n,1:n), with physical dimension np, this
c	routine constructs its Cholesky decomposition, A = L·LT . On input, only the upper triangle
c	of a need be given; it is not modified. The Cholesky factor L is returned in the lower triangle
c	of a, except for its diagonal elements which are returned in p(1:n).
c
	INTEGER i,j,k
	REAL sum
	do  i=1,n
	 do  j=i,n
	  sum=a(i,j)
	  do  k=i-1,1,-1
	   sum=sum-a(i,k)*a(j,k)
	  enddo 
	  if(i.eq.j)then
c	   sum=max(sum,0.0)
	   if(sum.le.0.)return
	   p(i)=sqrt(sum)
	  else
c	   if(p(i).eq.0.0)then
c 	    a(j,i)=0.0
c	   else
	    a(j,i)=sum/p(i)
c	   endif
	  endif
	 enddo 
	enddo 
	return
	END
c
