SUBROUTINE INIT1 CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC C In this revised version (Feb. 1998), the vertical profile of C C the target tangential wind is expressed in an empirical C C analytical form. The target vortex extends to SCAP level. C C SCAP level is assumed to be higher for stronger low level C C maximum wind but upper-bounded at a level which is slightly C C lower than the latitude dependent tropopause level. C C Time integration period for vortex generation is set to C C 60 hours, including one hour free-up time. C CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCC PARAMETER (IMX=250, JMX=1, KMAX=42, LGI=27 ,NP=26) parameter (nwat = 6, nwat3 = nwat + 3, nvar=nwat3) PARAMETER (KMAXP = KMAX+1, KMAXM = KMAX - 1, nnst=3) PARAMETER (IBLKMX=LGI*IMX+4*KMAX*IMX, IMXM=IMX-1,nmx=24) c include 'GDINFO.h' include 'COMMUNICATE.h' include 'CONMLEV.h' include 'CONSLEV.h' c COMMON /FILEC/ TYP1C(IMX,LGI),TYP2C(KMAX,IMX,nwat3) COMMON /SURFC/ TSMEAN,WTMEAN,WSIN,PC,PMIN,IOFF CC COMMON /NING/ TNGAS,NWRA COMMON /NSTP/ NSTEP1,NSTEP2,NRUN,INWRA,INWRB,NFMFLG,NWRB COMMON /FACT/ HURR,DEEP COMMON /BKINFO/ NSTEP,NNEST COMMON /COMPND/ FILEPC(IBLKMX),JUNK(IBLKMX,5) COMMON /BARF/ QQ(KMAX,IMX),QQ1(KMAX),QQ2(KMAX),PCHK COMMON /TANG/ VT(KMAX,IMX),TBAR(KMAX),RBAR(KMAX),PBAR,FCOR * ,PMXMN,IENV COMMON /SM/ CRW(NP),RC(NP),VC(NP),R(IMX),VCR(IMX),WC(IMX) COMMON /XXX/ XF(IMX,JMX),XC,YC,DXX,DYY COMMON /POSIT/ XOLD,YOLD,XCORN,YCORN,RO,XV,YV,ROMM,TSIZE,RB CCC COMMON /INGCUT/ TNGCUT,VMXCUT COMMON /MAXVAL/ RMAX, VMAX, SGMST, SCAPCT COMMON /STRMSZ/ TCHGHT, TCRDUS CCC DIMENSION FK(KMAX,IMX) ,ro(nmx) DIMENSION PSPR(IMX),AW(IMX) DIMENSION FR(IMX),TABS(8),TABL(8),ITYP1C(IMX,LGI) DIMENSION FILC(IBLKMX),WIND(IMX,4),NNQ(4),WINDA(IMX,8) CCC DIMENSION SGMS(KMAX),VPARA(IMX),SCAP(IMX),EXSS(IMX) DIMENSION FKX(KMAX) C DIMENSION RGALM(4) C REAL*8 A,B,RMAX EQUIVALENCE (FILC(1),TYP1C(1,1)),(TYP1C,ITYP1C) CHARACTER*12 IBLOCK DATA IBLOCK /'AXMD'/ CC DEFINE USEFUL CONSTANTS CCCCCC C DATA RBAR /.391E-5,.504E-5,.602E-5,.483E-4, C * .198E-3,.517E-3,.752E-3,.162E-2,.299E-2, C * .558E-2,.902E-2,.881E-2,.883E-2,.126E-1, C * .174E-1,.190E-1,195E-1,.202E-1/ C DATA TBAR / 219.23,205.18,198.08,214.85,230.14, C * 239.90,246.68,254.03,262.93,270.90, C * 278.28,281.95,283.92,289.58,294.30, C * 297.66,299.50,300.53/ CCCCCCCCCCCCCCCCCCC DATA TABL /.0,.35,.65,.82,.88,.97,.99,1./ DATA TABS /.15,.2,.3,.4,.5,.7,.85,1.0/ c CCC KTOP = 15 < level 18 ktop C C morristest BTOP = .88 BTOP = .84 DO 2284 K=1,KMAX IF(Q(K).GT.BTOP)GO TO 2011 2284 CONTINUE 2011 CONTINUE KTOP = K-1 print*,'This is the level above ktop: ',KTOP C PI = 4.*ATAN(1.0) PI180 = PI/180. C CCCCCCCCCCCCCCCCCCCCCCCCCC C&E CCCCCCCCCCCCCCCCCCCCCC C C set up the carr and elsberry formulation C REDEFINE RB USING CARR & ELSBERRY FORMULATION (1996): C NTOT = 0 VGALE = 0.0 RGALE = 0.0 C GF08 X = .4 C C increase x to make storms smaller C X = .415 READ(98)XV,YV READ(98)IQD READ(98)WNE,WSE,WSW,WNW READ(98)FATX,FATY,PC READ(98)HURR,DEEP READ(98)PLAST READ(98)IHOUR READ(98)TSIZE READ(98)XOLD,YOLD,XCORN,YCORN READ(98)RO,RB READ(98)NNQ(1),NNQ(2),NNQ(3),NNQ(4) DO 4044 IWQ = 1 , IQD IF(NNQ(IWQ).EQ.0)THEN PRINT*,'NO OBSERVATIONS FOR QUADRANT: ',IWQ ELSE NNN = NNQ(IWQ) + 1 C READ(98)(RC(I), I = 2 , NNN) READ(98)(WC(I), I = 2 , NNN) C if(YV .LT. 11.0)YV=11.0 NTOT = NTOT + 1 CC gf07 VGALE = 18.1 VGALE = VGALE + WC(NNN)*.01 RGALE = RGALE + RC(NNN)*.01 ENDIF CCC 4044 CONTINUE VGALE = VGALE/FLOAT(NTOT) RGALE = RGALE/FLOAT(NTOT) REWIND 98 C print*,'rgale, vgale, lat: ',RGALE,VGALE,YV OMEGA = 7.292E-5 FCOR = 2.*OMEGA*SIN(YV*PI180) C angm = rgale*vgale + .5*fcor*rgale*rgale angmp = angm/rgale**(1 - x) print*,'this is angular momentum: ',angm,angmp c alogf = alog(2.*angmp/fcor) expf = exp(alogf/1.4) RB = expf/1000.0 IF(RB.LT.500.0)THEN RB = 500.0 ENDIF RUSS = RB print* print* print*,'this is the new RB in km: ',RUSS C CCCCCCCCCCCCCCCCCCCCCCCCCC C&E CCCCCCCCCCCCCCCCCCCCCC C IOFF = 0 KUP = 4 DELC= 0.0 ARAD = 6.371E8 DEN = 1.15E-3 READ(3,3456)ILOW,IHIGH 3456 FORMAT(I2,1x,I2) TLOW = FLOAT(ILOW) THIGH = FLOAT(IHIGH) PRINT* PRINT*,'THESE ARE THE HIGH AND LOW LIMITS: ',IHIGH,ILOW CC C IFIRST = 1 C RESOL = 6.0 c biju RESOL = 12.0 C C C RC(1) = 0.0 WC(1) = 0.0 CC READ(98)XV,YV CC TVAR = 1.0 IF(YV.LT.0.0)TVAR=-1.0 CC YLAT = YV XLON = XV CC PRINT*,'XLON,YLAT IN INIT:',XLON,YLAT FACT = 1.0 IF(YLAT.LT.0.0)FACT = -1.0 CC READ(98)IQD CC CC CC C C*********************************************** C IF CONDITION IF READING IN 4 QUADRENTS OF DATA..... C IF(IQD.EQ.4)THEN READ(98)WNE,WSE,WSW,WNW ENDIF C C******* END OF IF CONDITION **************** C CC CC C READ(98)FATX,FATY,PC WRITE(6,3306)FATX,FATY,PC 3306 FORMAT(/,2X,'FINAL FACTORS:',3F8.2) CC CC CC CC********************************************************************* CC READ(98)HURR,DEEP IF(HURR.EQ.1.0)THEN WRITE(6,3505)DEEP 3505 FORMAT(2X,'HURRICANE WITH DEPTH FACTOR :',F7.1) ENDIF IF(HURR.EQ.0.0)THEN WRITE(6,3506)DEEP 3506 FORMAT(2X,'TROPICAL STORM WITH DEPTH FACTOR :',F7.1) ENDIF CC CC IF(DEEP.EQ.1.)THEN TABL(1) = .0 TABL(2) = .0 TABL(3) = .0 TABL(4) = .0 TABL(5) = .0 TABL(6) = .6 TABL(7) = .9 TABL(8) = 1.0 ENDIF IF(DEEP.EQ.2.)THEN TABL(1) = .0 TABL(2) = .0 TABL(3) = .0 TABL(4) = .0 TABL(5) = .5 TABL(6) = .85 TABL(7) = .95 TABL(8) = 1.0 ENDIF IF(DEEP.EQ.3.)THEN TABL(1) = .0 TABL(2) = .35 TABL(3) = .65 TABL(4) = .82 TABL(5) = .88 TABL(6) = .97 TABL(7) = .99 TABL(8) = 1.0 ENDIF CC CC**************************************************************** CC CC CC ITURN = 1 IF WE WANT TO TAKE TBAR AND RBAR FROM THE CENTER POINT VALUE CC CC ITURN = 1 CC C C*************************************************************** C C DETERMINE THE CONSTANTS THAT DECIDE WHEN WE SHOULD SHUT OFF THE C INTEGRATION UNITS OF (hPa): C CCC CCC set values of AAG, BBG, PBIG CCC values reset for 2006 operations C C AAG = 2.0 BBG = .10 PBIG = 5.5 C C READ(98)PLAST C C C PRINT*,'INITIAL MAX-MIN PRESSURES:',PLAST,PC C PDIFF = (PLAST - PC) PCORR = AAG + BBG*PDIFF IF(PCORR.GT.PBIG)THEN PCORR = PBIG ENDIF PMXMN = 1.0E3 * (PCORR + PLAST - PC) C PLAST = PLAST*1.0E3 C PRINT* PRINT*,'PLEASE NOTE THE FOLLOWING:' PRINT* PRINT* WRITE(6,3404)PMXMN,PCORR 3404 FORMAT(2X,' MAX-MIN PSTAR FROM THE PCORR: ',2X,2E12.6) PRINT* PRINT* CC C C************************************************************** C READ(98)IHOUR CC CC READ(98)TSIZE WRITE(6,3503)TSIZE 3503 FORMAT(/,2X,'PARAMETER FOR VORTEX SIZE:',F7.1) CC CC READ(98)XOLD,YOLD,XCORN,YCORN CC WRITE(6,3302)XV,YV,XOLD,YOLD,XCORN,YCORN 3302 FORMAT(/,2X,'COORDIANTES:',6E11.5) C READ(98)RO,RBOLD C WRITE(6,3303) WRITE(6,3304)RO,RBOLD 3303 FORMAT(/,2X,'RO AND RB RADIUS:') 4304 FORMAT(2X,5E11.5) CC RBMIN = RBOLD * 111.194 RBMAX = RBOLD * 1.33 * 111.194 IF(RUSS .LT. RBMAX)then RB = RUSS print*,'RB is not redefined : ',RUSS,RBMAX else RB = RBMAX print*,'RB is double lci : ',RUSS,RBMAX endif IF(RB .LT. RBMIN)then RSS=RB RB = RBMIN print*,'RB too small set to RBMIN : ',RSS,RBMIN endif CC RB = RBOLD * 1.175 * 111.194 if(RB .GT. 1000.)then RB = 1000. endif CCC RB = RBOLD * 1.16666 * 111.194 print*,'test :','RB set to 1.175*RBOLD' print*,'new rb: ',rb CALL PHASE CC W1 = .5 W1 = 0.0 C C PBAR = W1*PLAST + (1.-W1)*PBAR PRINT*,'THIS IS THE INITIAL SURFACE PRESSURE: ',PBAR CC PMAX = PBAR/1.0E3 C CC WRITE(6,4440)PBAR 4440 FORMAT(2X,'THE PRESSURE OF OUTER REGION :',E12.6) CC NNEST = 1 NSTEP = 0 CC**************************************** CC FCOR = .59512E-4 OMEGA = 7.292E-5 FCOR = 2.*OMEGA*SIN(YLAT*PI180) write(6,3340)FCOR 3340 format(2x,'fcor',e12.6) C C DRR = PI180/RESOL C C DR = COS(YLAT*PI180)*ARAD*PI180/RESOL c c R1 = RB*PI180*ARAD c c Modify for C & E c R1 = RB*1.0E5 c RSIZ = R1/DR IENV = IFIX(RSIZ) IMAX = IENV + 3 IIMAX = IMAX IE = IMAX+1 WRITE(6,3494)IENV,IMAX 3494 FORMAT(2X,'THIS IS IMAX FOR THE INTEGRATION: ',2I5) C IF(IMAX.GT.IMX)THEN WRITE(6,3420)IMAX,IMX 3420 FORMAT(2X,'IMAX :',I4,' IS LARGER THAN IMX: ',I5) WRITE(6,3427) 3427 FORMAT(2X,'WE MUST STOP THE INTEGRATION!!!!!!!!!') STOP 4301 ENDIF NGD = 1 NGR = 1 NTR = 1 DELTAT = 120./RESOL c biju cc DELTAT = 90./RESOL JS = 2 JN = 0 IE = IMAX+1 IW = 0 IIMAX = IIMAX IMAX = IMAX JJMAX = 1 JMAX = 1 NSTFLG = 0 ICX = 0 ICY = 0 IHX = 0 IHY = 0 DFTX = 0. DFTY = 0. CC CC************************************* CC FAR = 1./3. FACTR = 1./(Q(KTOP)-Q(3)) ROMM = R1 TLS = (.5*R1)**2 R1S = R1*R1 RLS = -R1S/TLS PN = PMAX CC CC************************************* CC C CC CC TFACT= 1.2 CC IF(IQD.EQ.4)THEN READ(98)NNQ(1),NNQ(2),NNQ(3),NNQ(4) ELSE READ(98)NNT NNN = NNT + 1 CCCC NNN = NNT N = NNN+10 ENDIF C DO 3044 IWQ = 1 , IQD C IF(NNQ(IWQ).EQ.0)THEN PRINT*,'NO OBSERVATIONS FOR QUADRANT: ',IWQ C RC(2) = 0.0 WC(2) = 0.0 C DO 3429 I = 1 , IMX WIND(I,IWQ) = 0.0 3429 CONTINUE ELSE C IF(IQD.EQ.4)THEN NNN = NNQ(IWQ) + 1 CCCC NNN = NNQ(IWQ) N = NNN+10 ENDIF C C C READ(98)(RC(I), I = 2 , NNN) READ(98)(WC(I), I = 2 , NNN) if(proc_num .EQ. 1) then WRITE(6,3304)(RC(I), I = 1 , NNN) WRITE(6,3305)(WC(I), I = 1 , NNN) 3304 FORMAT(/,2X,'RADIUSS:',(/,2X,10E11.5) ) 3305 FORMAT(/,2X,'WINDSS:',(/,2X,10E11.5) ) endif CC C RMAX = RC(2) VMAX = WC(2) CC CCC Set the integration length TNGAS to 60 hours. TNGCUT = 60.0 VMXCUT = 50.0*100.0 TNGAS = TNGCUT C CCC NSTEP2 = TNGAS*3600./DELTAT NRUN = NSTEP2 - NSTEP1 NWRA = NSTEP2 C CCC PRINT* PRINT*,' IMPORTANT PLEASE OBSERVE: ' PRINT* PRINT*,' THIS WILL BE HOW LONG WE INTEGRATE THE MODEL' PRINT* WRITE(6,4303)NSTEP2,TNGAS 4303 FORMAT(2X,'THE FINAL STEP: ',I5,' HOURS INTEGRATION: ',F4.1) PRINT* CC CCCC ENDIF RATIO = VMAX/RMAX CC FACQ = .3 C C C CC MEASURED WINDS WITH THE MEAN MOTION VECTOR REMOVED CC CRW(1) = 0.0 CRW(2) = 0.0 CRW(3) = 0.0 CRW(N) = 0.0 CC RC(N) = R1 WC(N) = 0.0 CC CC CC********************************************************************* CC RH = RC(N-10) CC WRH = WC(N-10) CC R1RO = R1 - RH ROV = WC(N-10)*RH ANGO = WC(N-10)*RH/R1RO C C C IF(IFIRST.EQ.1)THEN C C################################################################## C PRINT*,' SETTING UP THE INITIAL CONDITION' DELC = XV*PI180 THAC = YV*PI180 CC cc PNN = (PN-PC)*1.0E3 PNN=PDIFF*1.0E3 E = EXP(1.0) B = VMAX*VMAX*E*DEN/PNN A = DEXP(B*DLOG(RMAX)) C AKM = EXP(B*ALOG(RMAX/1.0E5)) WRITE(6,2300)A,B,VMAX,PC,RMAX 2300 FORMAT(2X,'THIS IS A,B FOR VMAX OF:',E11.5,'P-STAR MIN OF:', * E11.5,1X,'RADIUS OF MAXIMUM WIND:',E11.5,2X,'A AND B IS',2E11.5) C C R(1) = DR/2. TYP1C(1,12) = R(1) DO 17 I = 2 , IMX R(I) = R(I-1) + DR TYP1C(I,12) = R(I-1) + DR 17 CONTINUE DO 165 I = 2 , IMXM R2P = .5*(R(I+1)+R(I) ) R1P = .5*(R(I-1)+R(I) ) TYP1C(I,4) = PI*(R2P*R2P-R1P*R1P) 165 CONTINUE R2P = .5*(R(1) + R(2) ) TYP1C(1,4) = PI*R2P*R2P R2P = R(IMX) + .5*DR R1P = .5*(R(IMX) + R(IMX-1) ) TYP1C(IMX,4) = PI*(R2P*R2P - R1P*R1P) CC DO 170 I = 1 , IMX TYP1C(I,1) = 0.0 TYP1C(I,2) = DR TYP1C(I,3) = DR TYP1C(I,5) = XLON*PI180 + DRR * .5 + (I-1) * DRR TYP1C(I,6) = YLAT*PI180 ITYP1C(I,7) = I ITYP1C(I,8) = 1 TYP1C(I,9) = FCOR TYP1C(I,10) = TSMEAN TYP1C(I,11) = 0.0 TYP1C(I,13) = 0.0 TYP1C(I,14) = 0.0 TYP1C(I,15) = -99.0 TYP1C(I,16) = 1.0 TYP1C(I,17) = 0.0 TYP1C(I,18) = 0.0 TYP1C(I,19) = .06 170 CONTINUE C ENDIF C C################################################################ C DO 19 I = 1 , IMAX IF(R(I).GE.1.0E4)THEN C IF(TVAR.GT.0.0)THEN AW(I) = SQRT(A*B*PNN*EXP(-A/(R(I)**B))/(DEN*(R(I)**B)) + * R(I)*R(I)*TYP1C(I,9)*TYP1C(I,9)*.25) - .5*R(I)*TYP1C(I,9) CC ELSE CCCC AW(I) = SQRT(A*B*PNN*EXP(-A/(R(I)**B))/(DEN*(R(I)**B)) ) AW(I) = SQRT(A*B*PNN*EXP(-A/(R(I)**B))/(DEN*(R(I)**B)) + * R(I)*R(I)*TYP1C(I,9)*TYP1C(I,9)*.25) + .5*R(I)*TYP1C(I,9) ENDIF C C ELSE AW(I) = 0.0 ENDIF 19 CONTINUE CC CC C C CC******* DO THE WIND CORRECTION AT OUTER REGIONS******* CC NMAX = MAXMAG(RC(1),NNN) + 1 C PRINT*,'MAXIMUM POINT ',NMAX,RC(NMAX),WC(NMAX) NMN = NMAX CC C RFRT = 1.0E7 C C THE SIZE RADIUS OF MAXIMUM WIND WHERE WE C STOP FORCING THE HOLLAND PROFILE AND USE THE RANKINE VORTEX C INSIDE OF RMAX C RHOL = .60E7 IF(RC(NMAX).LT.RHOL)THEN FATTX = .6 FATTY = .7 ELSE FATTX = FATX FATTY = FATY ENDIF C WRITE(6,3353)FATTX,FATTY 3353 FORMAT(/,2X,'FACTORS REDEFINED:',2F9.3) CC C DO 780 NR = 1 , 9 NN = NR + N - 10 RC(NN) = RC(NMN) + FLOAT(NR)*(R1 - RC(NMN))/10. 780 CONTINUE DO 790 NR = 1 , 9 NN = NR + N - 10 CC WC(NN) = ROV/RC(NN) - (RC(NN) - RH)*(ANGO/RC(NN)) R1F = (RH/RC(NN) ) ** FATTX R2F = ( (R1 - RC(NN))/R1RO ) ** FATTY WC(NN) = WRH*R1F*R2F 790 CONTINUE CC VC(1) = 0.0 CC DO 777 NN = 2 , N IF(TVAR.GT.0.0)THEN VC(NN) = SQRT(A*B*PNN*EXP(-A/(RC(NN)**B))/(DEN*(RC(NN)**B)) + *RC(NN)*RC(NN)*TYP1C(5,9)*TYP1C(5,9)*.25)-.5*RC(NN)*TYP1C(5,9) ELSE CCCC VC(NN)= SQRT(A*B*PNN*EXP(-A/(RC(NN)**B))/(DEN*(RC(NN)**B))) VC(NN)= SQRT(A*B*PNN*EXP(-A/(RC(NN)**B))/(DEN*(RC(NN)**B)) + *RC(NN)*RC(NN)*TYP1C(5,9)*TYP1C(5,9)*.25)+.5*RC(NN)*TYP1C(5,9) ENDIF 777 CONTINUE c rayr: don't calculate difference vector if we have at least 2 obs. do 778 nn = 1 , n if (nnq(iwq) .ge. 2) then crw(nn)=wc(nn) else crw(nn) = wc(nn) - vc(nn) end if 778 continue c DO 778 NN = 1 , N c CRW(NN) = WC(NN) - VC(NN) c 778 CONTINUE CC if(proc_num .EQ. 1) then WRITE(6,2345)(VC(I) , I = 1 , N) WRITE(6,2345)(RC(I) , I = 1 , N) WRITE(6,2345)(WC(I) , I = 1 , N) 2345 FORMAT(2X,'BOGUS VORTEX AT THE MEASURED REGIONS:',(/,2X,2E11.5)) WRITE(6,2346)(CRW(I) , I = 1 , N) 2346 FORMAT(2X,'CORRECTION TO WINDS :',(/,2X,2E11.5)) endif CC do i = 1 , imx vcr(i) = 0.0 enddo c rayr: call monotonic spline routine if we have at least 2 obs. if (nnq(iwq) .ge. 2) then call akim2(n-1) else call akim(n) end if c if(proc_num .EQ. 1) then WRITE(6,3113)( R(I), I = 1 , IMX) c WRITE(6,3111)(AW(I), I = 1 , IMX) c WRITE(6,3111)(VCR(I), I = 1 , IMX) endif c WRITE(6,3113)( R(I), I = 1 , IMX) 3113 FORMAT(2X,'RR',(/,2X,10E11.5)) 3111 FORMAT(2X,'AW',(/,2X,10E11.5)) CC C C*******************CREATE AW(I): THE TARGET PROFILE******* C DO 779 I = 1 , IMAX CC IF(R(I).GE.RMAX)THEN c rayr: don't add back in Holland profile if we did monotonic spline if (nnq(iwq) .ge. 2) then aw(i)=vcr(i) * tvar else aw(i) = (vcr(i) + aw(i)) * tvar end if C c AW(I) = (VCR(I) + AW(I) ) * TVAR C C C ELSE C CCC IF(RMAX.GT.RHOL)THEN AW(I) = RATIO*R(I) * TVAR CCC ELSE CCC AW(I) = AW(I) CCC ENDIF C ENDIF CC 779 CONTINUE C C DO 269 I = 1 , IMX IF( R(I).GT.R1)AW(I) = .1 269 CONTINUE C if(proc_num .EQ. 1) then WRITE(6,3112)(AW(I), I = 1 , IMX) 3112 FORMAT(2X,'FINAL AW',(/,2X,10E11.5)) endif C C************************************************************ C DO 219 I = 1 , IMX TRR = (R1-R(I))**2 FR(I) = ( EXP(-TRR/TLS) - EXP(RLS) )/ * (1. - EXP(RLS)) 219 CONTINUE CC C C C DO 3187 I = 1 , IMX WIND(I,IWQ) = AW(I) 3187 CONTINUE C IFIRST = 0 C ENDIF C 3044 CONTINUE C DO 3993 I = 1 , IMX WINDA(I,1) = .5*(WIND(I,1)+WIND(I,2)) WINDA(I,2) = WIND(I,1) WINDA(I,3) = .5*(WIND(I,1)+WIND(I,4)) WINDA(I,4) = WIND(I,4) WINDA(I,5) = .5*(WIND(I,3)+WIND(I,4)) WINDA(I,6) = WIND(I,3) WINDA(I,7) = .5*(WIND(I,2)+WIND(I,3)) WINDA(I,8) = WIND(I,2) 3993 CONTINUE C DO 3994 ID = 1 , 8 WRITE(93)(WINDA(I,ID), I = 1 , IMX) 3994 CONTINUE REWIND 93 C IF(IQD.EQ.1)THEN WNE = 1.0 WSE = 0. WSW = 0. WNW = 0. ENDIF C DO 3197 I = 1 , IMX AW(I) = WNE*WIND(I,1) + WSE*WIND(I,2) + * WSW*WIND(I,3) + WNW*WIND(I,4) 3197 CONTINUE C WRITE(6,3111)( R(I), I = 1 , IMX) WRITE(6,3111)(AW(I), I = 1 , IMX) C FACT1 = 1.0 FACT2 = 1.0 FACT3 = 1.0 FACT4 = 1.0 IF(WNE.LT..05)FACT1 = 0.0 IF(WSE.LT..05)FACT2 = 0.0 IF(WSW.LT..05)FACT3 = 0.0 IF(WNW.LT..05)FACT4 = 0.0 C IF(IQD.EQ.4)THEN DO 3049 I = 1 , IMX A1 = FACT1*(WIND(I,1)/100. ) A2 = FACT2*(WIND(I,2)/100. ) A3 = FACT3*(WIND(I,3)/100. ) A4 = FACT4*(WIND(I,4)/100. ) AWW = FACT*(AW(I)/100. ) RWW = R(I)/1.0E5 write(75,2373)RWW,A1 write(76,2373)RWW,A2 write(77,2373)RWW,A3 write(78,2373)RWW,A4 write(79,2373)RWW,AWW 2373 FORMAT(2F8.3) 3049 CONTINUE ELSE DO 3051 I = 1 , IMX AWW = ABS(AW(I)/100.) RWW = R(I)/1.0E5 WRITE(79,*)RWW,AWW 3051 CONTINUE ENDIF C CCC *********************************************** CC CC CALCULATE THE CHANGE WITH HEIGHT OF THE VORTEX CCC -----------PROFILE FUNCTION----------- CC CCC DO 5100 K = 1 , KMAX SGMS(K) = Q(K) 5100 CONTINUE C Specify tropopause (at equator, pole, and YLAT) STEQ = 0.10 STPL = 0.25 STDIF = STPL - STEQ SGMST = STEQ + STDIF * 0.5 * (1.0 - COS(2.0*YLAT * *FACT*PI180)) C TC height (SCAPX) computation. SCAPX is bouded by C the upper-bound SCAPCT and the lower-bound SCAP00. C It increases with V, from V=0 to VCUT, and is set C to SCAPCT for V>VCUT. DELCT = 0.03 SCAPCT = SGMST + DELCT KPBLTP = KTOP - 1 SGMSB = Q(KPBLTP) SCAP00 = SGMSB C VCUT = THIGH*100.0 if(VMAX.LT.TLOW*100.0)THEN VMAX2 = TLOW*100.0 if(proc_num .EQ. 1) then print* print*,'caution, the vmax will be modified:' print* print*,'current value: ',VMAX print* print*,'new value: ',VMAX2 endif ELSE VMAX2 = VMAX if(proc_num .EQ. 1) then print* print*,'this is the current vmax: ',VMAX2 print* endif ENDIF C SCAPX by an empirical formula. IF(VMAX.GT.VCUT)THEN VPARAX = 1.0 ELSE VPARAX = VMAX2/VCUT ENDIF SANGL = PI*VPARAX SCAPX = SCAPCT + (SCAP00 - SCAPCT) * 0.5 * (1.0 + * COS(SANGL)) C Write TC size in /STRMSZ/ for use in SUBR. UPDATE. TCHGHT = SCAPX TCRDUS = R1 C Now, obtain the wind profile FK(K,I). C First, disposable parameter PWRP is specified. PWRP = 0.4 DO 5120 K = 1, KMAX IF(SGMS(K).LE.SCAPX)THEN FKX(K) = 0.0 ELSE SGX = (1.0 - SGMS(K))/(1.0 - SCAPX) FKX(K) = (1.0 - SGX * SGX) ** PWRP ENDIF 5120 CONTINUE DO 5230 I = 1 , IMX DO 5220 K = 1 , KMAX FK(K,I) = FKX(K) 5220 CONTINUE 5230 CONTINUE 5250 CONTINUE C DO 21 I = 1 , IMX TRR = (R1-R(I))**2 FR(I) = ( EXP(-TRR/TLS) - EXP(RLS) )/ * (1. - EXP(RLS)) 21 CONTINUE DO 26 I = 1 , IMX IF(R(I).GT.R1)FR(I) = 1.0 26 CONTINUE C CCC**************************************************** CCC Compute target wind VT(K,I). G = 0.0 DO 919 I = 1 , IMX DO 919 K = 1 , KMAX VT(K,I) = FK(K,I)*AW(I) TYP2C(K,I,2) = FK(K,I)*AW(I)*G TYP2C(K,I,1) = 0.0 TYP2C(K,I,3) = TBAR(K) TYP2C(K,I,4) = RBAR(K) 919 CONTINUE DO 920 I = 1 , IMX TYP1C(I,LGI) = PBAR 920 CONTINUE CC if(proc_num .EQ. 1) then WRITE(6,2340)(RBAR(K) , K = 1 , KMAX) WRITE(6,2340)(TBAR(K) , K = 1 , KMAX) 2340 FORMAT(2X,'INIT FIELD:',(/,2X,10E11.5)) WRITE(6,2341)PN 2341 FORMAT(2X,'INIT PRESSURE:',2X,E11.5) endif C DO 888 I = 1 , IBLKMX FILEPC(I) = FILC(I) 888 CONTINUE CC SGX = (1.0 - SGMS(K))/(1.0 - SCAPX) CC write(6,4050)(VT(30,i), i = 1 , imax) 4050 format(2x,'leaving init ',(10e11.5) ) CC CC RETURN END