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