subroutine tcBogus(zgaIO,tgIO,qgIO,ugIO,vgIO,dim, & centralPressure,radius,oldPosition_ll,newPosition_ll) use Dimensions_Type use Constants , except=>pi ! (PACKAGE: Dynamics ) ! TROPICAL CYCLONE BOGUS BASED ON THE JMA METHOD OF UENO, 1987. ! MODIFIED AND ENHANCED FOR USE IN THE TASP SYSTEM ! VERSION AUTHOR/DATE COMMENTS ! 0 UENO 1987 JMA ORIGINAL CODE ! 1 K.KURIHARA,J.WADSLEY, MODIFIED FOR TASP SYSTEM ! N.DAVIDSON 1992 ! 2 R.BOWEN SEP 1992 TIDIED CODE ! 3 R.BOWEN JUN 1993 ADDED INPUT PARAMETERS ! 4 R.MCTAGGART-COWAN, ! S.WEESE 2003 FIXED CODE - UPGRD TO F90 - INTRO TO SPA ! KM = MAX NO OF ANALYSIS LEVELS USED !******************************************************************** ! COMPILING INSTRUCTIONS: USE "f90 -fixedform code.f" ! REFER TO "NOTE"'S IN THIS SCRIPT TO IDENTIFY PARAMETERS THAT ! CAN BE MANUALLY CHANGED FOR SUCCESSFUL OPERATION ! ******************************************************************* ! NOTE the inverted grid for j; this is to agree with the ! jma grid specification, which assumes the j index starts at one ! in the north and proceeds to the south ! The following units are used by the bogusing routine: ! tg .... temperature at each level in deg. K ! ug .... zonal wind in m/s ! vg .... meridional wind in m/s ! zga .... geopotential height (m) ! qg .... specific humidity g/kg ! pmsl ... mean sea level pressure in hPa ! !******************************************************************** ! Input variables type(dimensions), intent(in) :: dim !grid specifications real, intent(in), optional :: centralPressure, & radius !bogus hurricane central pressure,radius (defaults 970.hPa and 1000.km) integer, dimension(2), intent(in), optional :: & oldPosition_ll,newPosition_ll !original hurricane location (lat,lon), bogus hurricane position (lat,lon - default to old) ! Input/output variables real, dimension(:,:,:), intent(inout) :: & zgaIO,tgIO,qgIO,ugIO,vgIO !heights(dam), temperature(degC), specific humidity(kg/kg), u,v winds(m/s) ! End_Header# ! Internal variables ! NLVM = NO LEVS FOR MOIS IN BOXS PARAMETER ( NLVM=6 ) ! MAXBOG = MAX NO OF OBS PER PRESS LEV PER ELEMENT PARAMETER ( MAXBOG=500 ) ! MXTY = MAXIMUM NUM OF TC'S THAT CAN BE PROCESSED PARAMETER ( MXTY=10 ) PARAMETER ( MNSF=50, MNMF=50 ) PARAMETER ( LU1=5, LU2=6, LU3=10, LU4=20, LU5=30, LU6=40 ) PARAMETER ( PI=3.14159265358979 ) ! INTEGER PRESS(dim%nk), LSG(dim%ni,dim%nj) INTEGER INDSL(MNMF), INDML(MNMF,dim%nk), SFLD(MNSF) INTEGER MFLD(MNMF) REAL SPLA(dim%nk), DELSPL(dim%nk) REAL DATA(dim%ni,dim%nj), PHISG(dim%ni,dim%nj) REAL PSEAG(dim%ni,dim%nj), WORK(dim%ni,dim%nj,dim%nk) REAL SSTG(dim%ni,dim%nj), FLATG(dim%ni,dim%nj) REAL FLONG(dim%ni,dim%nj), RSMG(dim%ni,dim%nj) REAL BOGLAT(MAXBOG), BOGLON(MAXBOG), WOBS(dim%nk) REAL PAIG(dim%ni,dim%nj), WG1(dim%ni,dim%nj), WG2(dim%ni,dim%nj) REAL USERU(MXTY), USERV(MXTY) real pmsl(dim%ni,dim%nj), SLAT real, dimension(dim%ni,dim%nj,dim%nk) :: zga,tg,qg,ug,vg CHARACTER NPRO*4, DDN*8, MOTTYP*4, LOCBY*4 real :: myCentralPressure,myRadius integer, dimension(2) :: myOldPosition,myNewPosition integer, dimension(2) :: oldPosition,newPosition real, dimension(2) :: oldLatLong,newLatLong real, dimension(dim%ni,dim%nj) :: myLong,myLat EQUIVALENCE ( IOBDAY, OBSDAY ) EQUIVALENCE ( IOBTIM, OBSTIM ) ! EQUIVALENCE ( IDIREC, DIREC ) ! COMMON /TYDATA/ NTY, NUM(MXTY), TPC(MXTY), T15(MXTY), TLO(MXTY) & ,TLA(MXTY), UTY(MXTY), VTY(MXTY), XTY(MXTY) & ,YTY(MXTY), TCPC(MXTY) ! CHARACTER*4 SGMA, MSLP, TEMP, UCMP, VCMP, HGHT, TDPT, MIXR DATA SGMA / 'SGMA' /, MSLP / 'MSLP' /, TEMP / 'TEMP' / DATA UCMP / 'UCMP' /, VCMP / 'VCMP' /, HGHT / 'HGHT' / DATA TDPT / 'TDPT' /, MIXR / 'MIXR' / ! Bogus parameter character(len=132), parameter :: myFile="boguspar.cfg" character(len=70) :: string real :: basicp,crmtnp,prfrac,cpdif,rmv integer, parameter :: inUnit=1 integer :: err ! Get bogus data open (inUnit,file=myFile,status='OLD',iostat=err) if (err.ne.0) then print*, "FATAL: Could not find settings file ",myFile stop endif ! REDEFINE INPUTS im = dim%ni jm = dim%nj km = dim%nk dels = dim%dx slat = -99. ! initialize select case (dim%mapProj) case ('E') npro = 'MER' slat = dim%lat(dim%ni/2,dim%nj/2) ! map scale == 1 case ('M') npro = 'MER' slat = dim%lat(dim%ni/2,dim%nj/2) ! map scale == 1 case ('B') npro = 'LL' case ('A') npro = 'LL' case ('N') npro = 'PSN' slat = 60. ! map scale == 1 case ('S') npro = 'PSS' case DEFAULT print*,'WARNING: Map projection not recognized by bogus ', & 'program' print*,'Setting to default Mercator projection' print*,'mapProj = ',dim%mapProj npro = 'MER' end select !!Z In mc2 settings mapProj is "L", but here it is "E" !!Z So before this problem be resolved, let mapProj="L" !!Z npro = 'LL' ! print *, 'npro ',npro ! pause !!! press = dim%levs do j=1,dim%nj zga(:,j,:) = zgaIO(:,dim%nj-j+1,:)*damtom tg(:,j,:) = tgIO(:,dim%nj-j+1,:)+to ug(:,j,:) = ugIO(:,dim%nj-j+1,:)*kttoms vg(:,j,:) = vgIO(:,dim%nj-j+1,:)*kttoms qg(:,j,:) = qgIO(:,dim%nj-j+1,:)*1000. myLat(:,j) = dim%lat(:,dim%nj-j+1) myLong(:,j) = dim%long(:,dim%nj-j+1) enddo !!z !!Z Recalculating the ij of the position, !!z oldPosition(1) = int(abs((myLong(1,1) - oldPosition_ll(1)/10.))/(DELS/111000.00))+1 oldPosition(2) = dim%nj-int(abs((myLat(1,dim%nj) - oldPosition_ll(2)/10.))/(DELS/111000.00))-1 newPosition(1) = int(abs((myLong(1,1) - newPosition_ll(1)/10.))/(DELS/111000.00))+1 newPosition(2) = dim%nj-int(abs((myLat(1,dim%nj) - newPosition_ll(2)/10.))/(DELS/111000.00))-1 write(*,*) 'Old position = ', oldPosition(1),oldPosition(2) write(*,*) 'New position = ', newPosition(1),newPosition(2) write(*,*) 'grid res = ',DELS write(*,*) 'Projection type (npro)= ',npro pause do j=1,dim%nj do i=1,dim%ni if (myLong(i,j) < 0.) myLong=360.+myLong enddo enddo pmsl = 1000.*exp(9.8*zga(:,:,1)/(281.*tg(:,:,1))) if (present(centralPressure)) then myCentralPressure = centralPressure else myCentralPressure = 970. endif if (present(radius)) then myRadius = radius else myRadius = 1.e3 endif if (present(oldPosition_ll)) then myOldPosition = oldPosition else myOldPosition(1) = dim%ni/2; myOldPosition(2) = dim%nj/2 endif oldLatLong(1) = myLat(myOldPosition(1),myOldPosition(2)) oldLatLong(2) = myLong(myOldPosition(1),myOldPosition(2)) if (present(newPosition_ll)) then myNewPosition = newPosition newLatLong(1) = myLat(newPosition(1),newPosition(2)) newLatLong(2) = myLong(newPosition(1),newPosition(2)) else myNewPosition = -1 newLatLong(1) = 0.; newLatLong(2) = 0. endif ! END REDEFINE INPUT NL = KM NBGFLS=1 ! ! NCIRC = NO OF CIRCLES ON WHICH TO GENERATE BOGUS OBS ! NBAD = DISTANCE (DEG LAT) BETWEEN CIRCLES ! NDIR = NO OF DIRECTIONS IN CIRCLE FOR OBS NCIRC=3 BRAD=1.25 NDIR=8 IPP=1 ! READ(LU1,*) NCIRC, BRAD, NDIR, IPP ! print*,NCIRC, BRAD, NDIR, IPP IF ( NCIRC*NDIR .GT. MAXBOG ) THEN WRITE(LU2,*) '*** FATAL ERROR : NO OF BOGUS OBS IS TOO GREAT', & ' MAXBOG = ', MAXBOG WRITE(LU2,*) ' REDUCE NCIRC AND/OR NDIR', & ' IN INPUT OR EDIT PROGRAM !!!' stop ENDIF ! GENERATE LAT/LON POSITIONS OF DATA ON CIRCLES AND DIRECTIONS IF ( IPP .EQ. 1 ) THEN NBOGS = 1 BOGLAT(1) = 0.0 BOGLON(1) = 0.0 ELSEIF ( IPP .EQ. 2 ) THEN NBOGS = 0 ENDIF DO 10 NC = 1, NCIRC DO 10 NX = 1, NDIR NBOGS = NBOGS + 1 ANG = ( NX - 1 ) * PI * 0.25 BOGLAT(NBOGS) = NC * BRAD * COS ( ANG ) BOGLON(NBOGS) = NC * BRAD * SIN ( ANG ) ! PRINT *,' LT,LN RLTV TO CNTR ',NBOGS,BOGLAT(NBOGS),BOGLON(NBOGS) 10 CONTINUE ! MOTION DESCRIPTION MOTTYP='12HR' do i=1,mxty useru(i)=0. userv(i)=0. enddo ! CYCLONE LOCATOR ('WIND' OR 'PRES') LOCBY='WIND' ! BASICP, CRMTNP = 0.0 TO 1.0 OF BASIC FLOW & MOTION FLOW TO USE ! basicp=0.75 ! original setup of Ron ! crmtnp=0.5 ! PRFRAC = FRACTION OF ANAL P TO USE ! prfrac=0.25 ! PRESSURE DIFFERENCE / RADIUS OF MAX WIND (km) ! cpdif=50 ! rmw=95 read(inUnit,*) string, basicp,crmtnp read(inUnit,*) string, prfrac read(inUnit,*) string, cpdif,rmw write(6,*) basicp,crmtnp write(6,*) prfrac write(6,*) cpdif,rmw close(inUnit) ! Fill spla array for use as a pressure index DO 30 K = 1, NL SPLA(K) = PRESS(K) 30 CONTINUE DO 40 K = 2, NL - 1 DELSPL(K) = ( - SPLA(K+1) + SPLA(K-1) ) * 0.5 40 CONTINUE DELSPL(1) = - SPLA(2) + SPLA(1) DELSPL(NL) = - SPLA(NL) + SPLA(NL-1) ! SET UP DOMAIN PARAMETERS - PROJECTION ! *************************** ! ** PROJECTION PARAMETERS ** ! *************************** ! NOTE: here specify the reference point for the lat/long grid. ! reference point: give an arbitrary grid point and the ! corresponding latitude and longitude at that point. ! standard lat/lon: specify the latitude and longitude ! where the map scale factor is 1.0, i.e. for a rotated ! mercator projector, the center of the grid ! if (npro .ne. 'LL') then XI = dim%ni/2; XJ = dim%nj/2 XLAT = myLat(dim%ni/2,dim%nj/2) XLON = myLong(dim%ni/2,dim%nj/2) SLON = myLong(dim%ni/2,dim%nj/2) !map scale factor=1. at centre?... ok else XI = 1 XJ = jm XLAT = myLat (1,jm) XLON = myLong(1,jm) SLON = myLong(1,jm) endif call qqrelh ( qg, qg, tg, spla, im, jm, nl, 'G/KG', 'QTOH' ) DO 50 J = 1, JM DO 50 I = 1, IM PSEAG(I,J) = pmsl(i,j) 50 CONTINUE DO 45 I = 1, IM DO 45 J = 1, JM PAIG(I,J) = PSEAG(I,J) SSTG(I,J) = TG(I,J,1) LSG(I,J) = 0 PHISG(I,J) = 0 DO 45 K = NLVM+1, NL IF ( SPLA(K) .LT. 100.0 ) THEN QG(I,J,K) = 5 ELSE QG(I,J,K) = QG(I,J,NLVM) ENDIF 45 continue PTOP = 0.0 ! convert rel hum. to spec. hum. call qqrelh ( qg, qg, tg, spla, im, jm, nl, 'G/KG', 'HTOQ' ) ! SETTING FLATG FLONG CALL LATLON ( FLATG, FLONG, IM, JM, SLAT, & NPRO, DELS, SLON, XI, XJ, XLAT, XLON ) ! ******************************* ! ** READ IN BOGUS INFORMATION ** ! ******************************* CALL BOGSRD ( NBGFLS, NTY, NUM, TPC, T15, TLO, TLA, UTY, VTY, & XTY, YTY, NPRO, DELS, SLON, XI, XJ, XLAT, XLON, & IDATE, ITIME, UCOMP, VCOMP, MOTTYP, USERU, USERV, & CRMTNP, IPP, TCPC, SLAT ,oldLatLong(1),oldLatLong(2), & myCentralPressure,myRadius) CALL MAPFCT ( WG1, IM, JM, FLATG, NPRO, SLAT ) DO 55 J = 1, JM DO 55 I = 1, IM RSMG(I,J) = 1.0 / WG1(I,J) 55 CONTINUE ! >> TYPHOON TRANSPLANTATION << CALL TYPLNT ( PSEAG, PAIG, UG, VG, TG, QG, SSTG, & FLATG, FLONG, PTOP, IM, JM, NL, & PHISG, ZGA, DELSPL, SPLA, DELS, & XI, XJ, XLAT, XLON, SLON, NPRO, & CPDIF, RMW, LOCBY, & SLAT, BASICP, CRMTNP, PRFRAC, WORK, IPP , & myNewPosition, newLatLong, dim) WRITE(LU2,9000) 9000 FORMAT(//,10X,'BOGUS COMPLETED SUCCESSFULLY'/ & 10X,'----------------------------') !!!!! CONVERT TO OUTPUT FORMAT do j=1,dim%nj zgaIO(:,j,:) = zga(:,dim%nj-j+1,:)*mtodam tgIO(:,j,:) = tg(:,dim%nj-j+1,:)-to ugIO(:,j,:) = ug(:,dim%nj-j+1,:)*mstokt vgIO(:,j,:) = vg(:,dim%nj-j+1,:)*mstokt qgIO(:,j,:) = qg(:,dim%nj-j+1,:)/1000. enddo !!!!! END CONVERT ! STOP ! END RETURN END SUBROUTINE TCBOGUS !------------------------------------------------------------------- !******************************************************************** !------------------------------------------------------------------- SUBROUTINE AXISYM ( PSEA, T, GRAD, RH, POUTLG, IRM, KPM, DR, & R0, R3, R15, P3, PC, T3, TSEA, RH3, PIN, KM, & W1A, W1B, KM1, DZ, TP, TPV, DTC, TCL, TCLV, & QCL, RHP, POUT, IPP, TCPC, RMW ) ! +----------------------------------------------------------------+ ! + -------------------------------------------------------------- + ! + | CONSTRUCTION | + ! + | OF | + ! + | 2-DIMENSIONAL TYPHOON MASS FIELD ON ISOBARIC SURFACES | + ! + | BASED ON FUJITA'S FORMULA AND SUGI'S FORMULATION | + ! + | CREATED BY T. IWASAKI 1985/10/01 | + ! + | REFORMED BY M. UENO 1987/12/10 | + ! + -------------------------------------------------------------- + ! + + ! + INPUT + ! + PIN ** INPUT P-LEVEL GIVEN BY MOTHER PROGRAM + ! + + ! + PC : CENTRAL PRESSURE + ! + R0 : CHARACTERISTIC RADIUS IN FUJITA'S FORMULA + ! + TSEA : SEA SURFACE TEMPERATURE (REAL) + ! + R15 : 15 M/S WIND RADIUS + ! + P3 : ENVIRONMENT PSEA + ! + R3 : MARGINAL RADIUS OF 2-D FRAME + ! + T3 : ENVIROMENT TEMPERATURE (REAL) + ! + RH3 : ENVIROMENT RELATIVE HUMIDITY + ! + + ! + OUTPUT + ! + POUTLG ** OUTPUT LOG(P)-LEVEL DEFINED IN THIS PROGRAM + ! + + ! + PSEA: SURFACE PRESSURE (1-D) + ! + T : TEMPERATURE (2-D) (REAL) + ! + RH : RELATIVE HUMIDITY (2-D) + ! + GRAD: PRESSURE GRADIENT (2-D) =DZ/DR + ! + + ! + WORK ! + W1A, W1B, ................ , POUT + ! + + ! + + ! + + ! + + ! + PTOP ----------------------- + ! + |-|-|-|-|-|-|-|-|-|-|-| + ! + (KPM) |-|-|-|-|-|-|-|-|-|-|-| + ! + |-|-|-|-|-|-|-|-|-|-|-| + ! + P3 ----------------------- + ! + 0 (IRM) R3 + ! + + ! + HORIZONTAL GRID SPACING : DR=R3/(IRM-1) + ! + VERTICAL SPACING : DLOGP = (ALOG(P3/PTOP))/(KPM-1) + ! + + ! -----------------------------------------------------------------+ ! DIMENSION PSEA(IRM), T(IRM,KPM), RH(IRM,KPM) & ,GRAD(IRM,KPM), POUT(KPM) & ,T3(KM), RH3(KM), PIN(KM) & ,W1A(KM1+1), W1B(KM1+1) & ,DZ(IRM,KPM) & ,TP(KPM), TPV(KPM), DTC(KPM) & ,TCL(KPM), TCLV(KPM), QCL(KPM) & ,RHP(KPM), POUTLG(KPM) ! ! ! ***** PHYSICAL PARAMETERS ***** GRAV = 9.8 RGAS = 287.04 RBYG = RGAS/GRAV GBYR = 1./RBYG EPSL = 0.608E-3 EPS = 0.608 ! ! ! ***** MODEL ARBITRARY PARAMETERS ***** ! ----------------------------------------------------------------- ! RATIO OF DZ(TROPOPAUSE) TO DZ(SURFACE) FWC = -0.5 ! |------ WARM CORE FACTOR ! DZ/DR=0 AT PMID (MB) PMID = 20. ! |------ MID-STRATOSPHERE LEVEL ! ---------------------------------------------------------------- ! PRES. (MB) OF TOP-LEVEL OF 2-D FRAME PTOP = PMID DR = R3/(IRM-1) ! CONNECTING RADII OF FORMULAE FOR DZ(LCLT) HORIZONTAL PROFILE R1 = R0*3. R2 = R3/1.5 ! ***** DEFINE OUTPUT P-LEVEL ***** ! LOG(P) LINEAR SPACING FROM SEA SURFACE TO PTOP DLOGP = (ALOG(P3/PTOP))/(KPM-1) DO 1 K=1,KPM POUTLG(K)= ALOG(P3)-DLOGP*(K-1) POUT(K) = EXP(POUTLG(K)) 1 CONTINUE ! ***** VERTICAL INTERPOLATION ***** ! ---- TEMPERATURE ---------- W1A(1) = ALOG(P3) W1B(1) = TSEA DO 2 K=2,KM+1 W1A(K) = ALOG(PIN(K-1)) W1B(K) = T3(K-1) 2 CONTINUE !!$ print*, 'w1a: ',w1a !!$ print*, 'w1b: ',w1b !!$ print*, 'poutlg: ',poutlg !!$ print*, 'kpm: ',kpm !!$ print*, 'km1: ',km1 CALL SPLINE (TP,POUTLG,KPM, W1B,W1A,KM1,2) !!$ print*, 'tp: ',tp ! ---- RELATIVE HUMIDITY ------ W1B(1) = 90. DO 3 K=2,KM+1 W1B(K) = RH3(K-1) 3 CONTINUE CALL SPLINE (RHP,POUTLG,KPM, W1B,W1A,KM1,2) ! ********************************************** ! VERTICAL TEMPERATURE PROFILE AT TYPHOON CENTER ! ********************************************** PINV = 1./PC ! MIXING RATIO G/G CALL TETENS(TSEA,PINV,QSAT,DQSAT) AXX = RBYG*TSEA*(1.+0.608*QSAT) DZSRF =-AXX*ALOG(P3/PC) TCLB = TSEA 115 CONTINUE ! ------ CONSTRUCT MOIST ADIABAT ------------- CALL TCLOUD (TCL,QCL,POUT,KPM,PC,TCLB,KPSEA) ! ! KPSEA : LOWEST LEVEL OVER SEA SURFACE ! ------ FIND CLOUD TOP LEVEL (LCLT) --------- ! (( CLOUD-TOP-LEVEL MUST BE ABOVE 400 MB )) DO 110 K=KPSEA,KPM LCLT=K IF(POUT(K).LT.400 .AND. TCL(K).LT.TP(K)) GO TO 111 !!$ print*, 'in axisym',k,pout(k),tcl(k),tp(k) 110 CONTINUE 111 CONTINUE PCLT = POUT(LCLT) ! WRITE(6,800) PCLT ! 800 FORMAT( / ,' CLOUD-TOP-LEVEL =',F10.1,'(MB)' / ) IF (LCLT.GE.KPM-2) THEN PRINT *,'CLOUD TOP LEVEL IS TOO HIGH' write(*,*)'LCLT=',LCLT,'KPM=',KPM stop ENDIF ! ----- CALCULATE DZZMAX AT TYPHOON CENTER ----- ! ........... DZZ=DZ(LCLT)-DZSRF ................ DO 112 K=1,KPM TK = TP(K) ES = 6.11*10.**(7.5*(TK-273.2)/(TK-35.9)) ! G/KG QQ = 622.*ES*0.01*RHP(K)/POUT(K) TPV(K) = TK * (1. + EPSL*QQ) TCLV(K) = TCL(K) * (1. + EPS *QCL(K)) 112 CONTINUE DZZMAX = (TCLV(KPSEA) - TPV(KPSEA)) * 0.5 & * (POUTLG(KPSEA)-ALOG(PC)) * RBYG ! PRINT *, ' TP(1) NEAR DZZMAX ', tp(1),tpv(1) DO 120 K=KPSEA,LCLT-1 DZZMAX = DZZMAX+(TCLV(K)-TPV(K)+TCLV(K+1)-TPV(K+1)) * 0.5 & * (-DLOGP) * RBYG 120 CONTINUE ! -------- SET CTROP AND CSTRA ------------------------------------ ! CTROP : DISTRIB. CONST. IN TROPOSPHERE ! CSTRA : DISTRIB. CONST. IN STRATOSPHERE CTROP = (1.-FWC) * DZSRF/DZZMAX DLG = ALOG(PMID/PCLT) CSTRA = 6.*GBYR*FWC*DZSRF / ( DLG**3 ) ! WRITE(6,801) DZSRF,DZZMAX,CTROP,CSTRA 801 FORMAT( / ' ( DZSRF,DZZMAX,CTROP,CSTRA ) =',5X,4F10.2 / ) ! ------ CHECK CTROP -------------------- IF(CTROP.LE.0.0 .OR. CTROP.GE.0.8) THEN ! CLOUD-BASE-TEMP. IS TOO LOW TO CONSTRUCT WARM-CORE CCC TCLB = TCLB+0.1 ! WRITE(6,802) TSEA,TCLB 802 FORMAT( / ,' CLOUD BASE TEMP. IS INCREASED FROM',F10.2,'K' & ,'TO',F10.2,'K' / ) GO TO 115 ENDIF ! ------ VIRTUAL TEMP. DEVIATION AT CENTER ----------------- DO 10 K=1,KPM DTC(K) = 0.0 10 CONTINUE DO 130 K=KPSEA,LCLT DTC(K) = CTROP*(TCLV(K)-TPV(K)) 130 CONTINUE DO 131 K=LCLT+1,KPM DTC(K) = CSTRA*(POUTLG(K) - ALOG(PCLT)) & *(ALOG(PMID) - POUTLG(K)) 131 CONTINUE ! ////////// MONITOR ////////// DO 150 K=1,KPM QCL(K) = TPV(K)+DTC(K) 150 CONTINUE ! WRITE(6, *)'REFERENCE LEVEL PRES. IN SUBR.((AXISYM))' ! WRITE(6,900) POUT ! WRITE(6, *)'TV DEVIATION AT CENTER' ! WRITE(6,900) DTC ! WRITE(6, *)'TV SPECIFIED AT CENTER' ! WRITE(6,900) QCL 900 FORMAT((5X,20F6.1)) ! /////////////////////////////// ! ************ ! CALCULATE DZ ! ************ ! ***** VERTICAL DZ PROFILE AT TYPHOON CENTER ***** DZ(1,1) = DZSRF DO 140 K=2,KPM DZ(1,K) = DZ(1,K-1) & - RBYG*0.5*(DTC(K)+DTC(K-1))*(POUTLG(K)-POUTLG(K-1)) 140 CONTINUE ! ***** HORIZONTAL DZ PROFILE ***** ! --------- SURFACE -------------------------------- CALL HPROF1 ( DZ, PSEA, IRM, KPM, PC, P3, R0, R3, DR, AXX, & IPP, TCPC, RMW ) ! --------- CLOUD-TOP-LEVEL --------------------------- CALL HPROF2 (DZ, IRM,KPM, LCLT,R1,R2,R3,DR,FWC,DZSRF) ! ***** HORIZONTAL/VERTICAL INTERPOLATION ***** ! ------- PSEA > P > PCLT -------------------------- DO 300 I=2,IRM DFACT = 1. / (DZ(1,1)-DZ(1,LCLT)) AA = (DZ(I,1)-DZ(I,LCLT)) * DFACT BB = - (DZ(I,1)*DZ(1,LCLT) - DZ(I,LCLT)*DZ(1,1)) * DFACT DO 301 K=1,LCLT DZ(I,K) = AA*DZ(1,K) + BB 301 CONTINUE 300 CONTINUE ! ------ PCLT > P > PMID --------------------------- DO 302 I=2,IRM AA = DZ(I,LCLT)/DZ(1,LCLT) DO 303 K=LCLT+1,KPM DZ(I,K) = AA*DZ(1,K) 303 CONTINUE 302 CONTINUE ! ************************* ! VIRTUAL TEMPERATURE FIELD ! ************************* DINV = 1./DLOGP DINV2 = DINV * 0.5 DO 310 I=1,IRM DZDLGP = (-1.5*DZ(I,1) + 2.*DZ(I,2) - 0.5*DZ(I,3)) * DINV T(I,1) = GBYR*DZDLGP +TPV(1) DZDLGP = (1.5*DZ(I,KPM) - 2.*DZ(I,KPM-1) +0.5*DZ(I,KPM-2)) & * DINV T(I,KPM) = GBYR*DZDLGP +TPV(KPM) DO 311 K=2,KPM-1 DZDLGP = (DZ(I,K+1) - DZ(I,K-1)) * DINV2 T(I,K) = GBYR*DZDLGP + TPV(K) 311 CONTINUE 310 CONTINUE ! WRITE(6, *)'TV AT CENTER DERIVED FROM HEIGHT FIELD' ! WRITE(6,900) (T(1,K),K=1,KPM) ! **************** ! GRADIENT (DZ/DR) ! **************** DR2 = 1./(2.*DR) DO 320 K=1,KPM GRAD(1, K) = 0. GRAD(IRM,K) = 0. DO 320 I=2,IRM-1 GRAD(I,K) = (DZ(I+1,K) - DZ(I-1,K)) * DR2 * GRAV 320 CONTINUE ! *********** ! WATER VAPOUR ! *********** ! +++++ 3(VERTICAL) X 3(HORIZONTAL) REGIMES +++++ ! PSEA <----> PCLT <----> PTOPW <----> PMID ! CENTER <----> RW1 <----> RW2 <----> R3 ! MILLIBAR PTOPW = 50. ! METER RWMOD = 200.E+3 ! PERCENT RH0 = 90. ! PERCENT RH1 = 85. ! PERCENT RH2 = 5. PLGW = ALOG(PTOPW) PLGINV = 1./(POUTLG(LCLT)-PLGW) RW2 = MAX(R15,RWMOD) RW1 = RW2*0.5 RWINV = 1./(RW2-RW1) DO 330 K=1,KPM IF (K.LE.LCLT) THEN QCL(K) = RH0 ELSE QCL(K) = RH1*(POUTLG(K)-PLGW)*PLGINV + RH2 IF (POUTLG(K).LT.PLGW) QCL(K) = RH2 ENDIF 330 CONTINUE DO 335 I=1,IRM R = DR*(I-1.) IF(R.LT.RW1) THEN DO 331 K=1,KPM RH(I,K) = QCL(K) 331 CONTINUE ELSEIF ( R .LT. RW2 ) THEN DO 332 K=1,KPM RH(I,K) = (QCL(K)*(RW2-R)+RHP(K)*(R-RW1)) * RWINV 332 CONTINUE ELSE DO 333 K=1,KPM RH(I,K) = RHP(K) 333 CONTINUE ENDIF 335 CONTINUE ! ------------ RH --> WV , TV --> T ------------------------- DO 340 K=1,KPM DO 340 I=1,IRM TK = T(I,K) ES = 6.11*10.**(7.5*(TK-273.2)/(TK-35.9)) QQ1 = 622.*ES*0.01*RH(I,K)/POUT(K) TR = TK/(1.+EPSL*QQ1) ES = 6.11*10.**(7.5*(TR-273.2)/(TR-35.9)) QQ2 = 622.*ES*0.01*RH(I,K)/POUT(K) T(I,K)= TK/(1.+EPSL*QQ2) 340 CONTINUE RETURN END !------------------------------------------------------------------ !****************************************************************** !------------------------------------------------------------------ SUBROUTINE BFLOW ( UA, VA, U, V, IM, JM, KM, CLAT, CLON, RTG, & ANG, RO, DELS, XI, XJ, XLAT, XLON, SLON, & NPRO, ACCT, VT, VR, SLAT ) ! ///-----------------------------------------------/// ! /// EXTRACT AXISYMMETRIC WIND COMPONENT /// ! /// FROM TYPHOON AREA /// ! ///-----------------------------------------------/// PARAMETER ( NBAND = 2 , NSECT = 16 , NBM = 120 ) ! OUT DIMENSION UA (IM,JM,KM), VA (IM,JM,KM) & ! IN ,U (IM,JM,KM), V (IM,JM,KM) & ! IN ,RTG (IM,JM), ANG (IM,JM) & ! WORK ,ACCT(IM,JM) & ! WORK ,VT (NBM,KM), VR (NBM,KM) & ! WORK ,UF (28), VF (28) real, dimension(im,jm,km) :: saveU,saveV CHARACTER * 4 NPRO saveU = u saveV = v RAD = ASIN(1.)/90. RE = RO DO 10 I = 1, IM IF (RTG(I, 1) .LT. RE) RE = RTG(I, 1) IF (RTG(I,JM) .LT. RE) RE = RTG(I,JM) 10 CONTINUE DO 20 J = 1, JM IF (RTG( 1,J) .LT. RE) RE = RTG( 1,J) IF (RTG(IM,J) .LT. RE) RE = RTG(IM,J) 20 CONTINUE DISU = DELS/NBAND NB = RE/DISU + 1 ANGU = 360./NSECT ! ARBITRARY VIRTUAL DISPLACEMENT RINCRE = 50000. DO 30 K = 1, KM DO 30 N = 1, NBM VT(N,K) = 0.0 VR(N,K) = 0.0 30 CONTINUE ! -------------------------------------------------------------------- DO 1000 N = 2, NB ! -------------------------------------------------------------------- DIST = (N-1) * DISU ! ---------------------------- DO 500 L = 1, NSECT ! ---------------------------- SANG = (L-1) * ANGU CALL CLY2LL (ALAT,ALON,SANG,DIST,CLAT,CLON) CALL RLTLN (FI,FJ,ALAT,ALON,SLAT, & NPRO,DELS,SLON, XI,XJ,XLAT,XLON ) CALL INTPLZ (UF, FI,FJ, U,IM,JM,KM) CALL INTPLZ (VF, FI,FJ, V,IM,JM,KM) R = DIST + RINCRE CALL CLY2LL (ALAT,ALON,SANG,R,CLAT,CLON) CALL RLTLN (VI,VJ,ALAT,ALON,SLAT, & NPRO,DELS,SLON, XI,XJ,XLAT,XLON ) DJ = VJ - FJ DI = VI - FI Y = max(SQRT(DI**2+DJ**2),epsilonR()) ZSIN = DJ / Y ZCOS = DI / Y ! Ron(corresponding to SUB. RLTLN) ! ZCOS = -DI / Y ! original DO 100 K = 1, KM VT(N,K) = VT(N,K) + UF(K)*ZSIN - VF(K)*ZCOS VR(N,K) = VR(N,K) - UF(K)*ZCOS - VF(K)*ZSIN 100 CONTINUE ! ---------------------------- 500 CONTINUE ! ---------------------------- DO 150 K = 1, KM VT(N,K) = VT(N,K)/NSECT VR(N,K) = VR(N,K)/NSECT 150 CONTINUE ! -------------------------------------------------------------------- 1000 CONTINUE ! -------------------------------------------------------------------- CALL ROTANG (ACCT,IM,JM,RTG,ANG,CLAT,CLON, & NPRO,DELS,SLON,XI,XJ,XLAT,XLON,SLAT ) DO 40 K = 1, KM DO 40 J = 1, JM DO 40 I = 1, IM UA(I,J,K) = U(I,J,K) VA(I,J,K) = V(I,J,K) 40 CONTINUE DO 200 J = 1, JM DO 200 I = 1, IM ! ------------------------------ IF (RTG(I,J) .LT. RE) THEN ! ------------------------------ FX = RTG(I,J)/DISU + 1. CALL INTPL1 (UF, FX, VT, NBM,KM, NB) CALL INTPL1 (VF, FX, VR, NBM,KM, NB) ANGL = ACCT(I,J) * RAD DO 210 K = 1, KM UA(I,J,K) = U(I,J,K) - UF(K)*SIN(ANGL)+VF(K)*COS(ANGL) VA(I,J,K) = V(I,J,K) + UF(K)*COS(ANGL)+VF(K)*SIN(ANGL) 210 CONTINUE ! ------------------------------ ENDIF ! ------------------------------ 200 CONTINUE u = saveU v = saveV RETURN END !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- SUBROUTINE BOGSRD ( NBGFLS, NTY, NUM, TPC, T15, TLO, TLA, & UTY, VTY, XTY, YTY, NPRO, DELS, SLON, XI, XJ, & XLAT, XLON, IDATE, ITIME, UCOMP, VCOMP, & MOTTYP, USERU, USERV, CRMTNP, IPP, TCPC,SLAT, & LATC,LONC,CPT,ROCI) PARAMETER ( MXTY=10 ) PARAMETER ( LU3=19 ) INTEGER NUM(*) REAL TPC(*), T15(*), TLO(*), TLA(*) REAL UTY(*), VTY(*), XTY(*), YTY(*) REAL TCPC(*), USERU(*), USERV(*) REAL LATC, LONC, ROCI, CPT, LAT6, LON6, LAT12, LON12, ALATC, ALONC CHARACTER*4 NPRO, MOTTYP CHARACTER FUNC*3,DTYP*4 NTY = 0 ! DO 100 LU = LU3, LU3-1+NBGFLS ! READ(LU,1010,END=100) NDATE, NTIME ! 1010 FORMAT(I6,1X,I4) ! 55 CONTINUE ! LOOP READING BOGUS CYCLONES UNTIL FUNC='END' ! NOTE: READS INFORMATION FROM "fort.19", MUCH OF INFORMATION ! READ HAS NO APPARENT VALUE; SHOULD SET NDATE AND NTIME ! TO MATCH DATE AND TIME SPECIFIED IN PROGRAM BODY, SET ! ROCI TO NECESSARY VALUE (SEE CONDITIONS IN "IF" ! IN THIS ROUTINE, CPT AND VALU ARE DESIRED CENTRAL ! PRESSURE OF BOGUSED VORTEX, DTYP MUST BE "TCYC", ! LATC AND LONC SET CENTER OF LAT/LON BOX IN WHICH SEARCH ! FOR HURRICANE WILL BE CONDUCTED OVER ! MOST READ IN VARIABLES ARE CONFINED TO THIS ROUTINE !!$ NOTE: old READ statement for centre location (LATC,LONC) were here !!$ NOTE: old READ statement for central pressure and radius (CPT,ROCI) !!$ were here ! READ(LU,1005,END=100) ROCI, CPT, LAT12, LON12, LAT6, LON6, ! > ALATC, ALONC ! 1005 FORMAT(8F10.5) ! PRINT *,' BOGUS, DATE: ',NDATE,' TIME ',NTIME ! IF ( IDATE .NE. NDATE .OR. ITIME .NE. NTIME ) THEN ! PRINT *,' SKIP SINCE DOES NOT MATCH WITH ANAL DATE TIME' ! ELSEIF ( DTYP .EQ. 'TCYC' ) THEN ! PRINT *,' CYCLONE! TCYC FOR CORRECT DATE/TIME' NTY = NTY + 1 ! IF ( NTY .GT. MXTY ) THEN ! PRINT *,' YOU HAVE ASKED FOR TOO MUCH = ', NTY, ', CAN', ! > ' ONLY DO ',MXTY,' TYPHOONS ' ! GOTO 100 ! ENDIF ! PRINT 16,FUNC,ITYP,LATC,LONC,DTYP,LEVL,VALU ! 16 FORMAT(' FUNC ITYP LL DTYP LEVL VALU ', ! > A3,1X,I1,2(1X,F8.1),1X,A4,1X,I5,1X,F6.1) ! PRINT 26, ROCI, CPT, LAT12, LON12, LAT6, LON6, ALATC, ALONC ! 26 FORMAT(' ROCI CPT LL12 LL6 LL0 ',8F10.2) ! CHECK FOR DUPLICATE TC'S ! IF ( NTY .GT. 1 ) THEN ! DO 22 NT = 1, NTY-1 ! IF ( ABS(LATC-TLA(NT)) .LT. 1.0 .AND. ! > ABS(LONC-TLO(NT)) .LT. 1.0 ) THEN ! PRINT *, ' >> DUPLICATE OBS IN TC BOGUS DATA SO IGNORE' ! NTY = NTY - 1 ! GO BACK AND READ IN NEXT PAOB ! GOTO 55 ! ENDIF ! 22 CONTINUE ! ENDIF IF ( CPT .LT. 900.0 ) THEN PRINT *,' >>BOGUS PRESSURE TOO LOW, MUST BE >900 hPa' NTY = NTY - 1 ! GO BACK AND READ IN NEXT PAOB ! GOTO 55 ENDIF IF ( ROCI .LT. 200.0 ) THEN PRINT *,' >>TOO SMALL, ROCI RESET TO 200KM FROM ',ROCI,'KM' ROCI = 200.0 ELSEIF ( ROCI .GT. 1200.0 ) THEN PRINT *,' >>TOO LARGE, ROCI RESET TO 800KM FROM ',ROCI,'KM' ROCI = 1200.0 ENDIF MEANF = 0 UCOMP = 0.0 VCOMP = 0.0 ! IF ( ABS(LATC-LAT6) .GE. 10 .OR. ! > ABS(LONC-LON6) .GE. 10 ) THEN ! PRINT *,' 6 HR PAST MOTION IS BAD!' ! MEANF = 6 ! ENDIF ! IF ( ABS(LATC-LAT12) .GE. 20 .OR. ! > ABS(LONC-LON12) .GE. 20 ) THEN ! PRINT *,' 12 HR PAST MOTION IS BAD!' ! IF ( MEANF .EQ. 0 ) THEN ! MEANF = 12 ! ELSE ! IF ( CRMTNP .NE. 0.0 .AND. MOTTYP .NE. 'USER' ) THEN ! PRINT *,' FAILED CHECK ON 0, 6, 12 HOUR TC LOCATION' ! PRINT *,' SO NO BOGUS FOR THIS TC' ! NTY = NTY - 1 ! GOTO 55 ! ENDIF ! MEANF = 1 ! ENDIF ! ENDIF NUM(NTY) = NTY TLA(NTY) = LATC TLO(NTY) = LONC TPC(NTY) = CPT TCPC(NTY) = CPT CALL RLTLN ( FI, FJ, LATC, LONC, SLAT, & NPRO, DELS, SLON, XI, XJ, XLAT, XLON ) XTY(NTY) = FI YTY(NTY) = FJ T15(NTY) = ROCI ! << INITIAL 6-HOUR MOVEMENT >> ! ! IF ( MEANF .NE. 1 ) THEN ! CALL TYMOVE ( UCOMP, VCOMP, NPRO, DELS, SLON, XI, XJ, ! > XLAT, XLON, LATC, LONC, LAT6, LON6, ! > LAT12, LON12, UCMP2, VCMP2, UCMP3, VCMP3 ) ! ENDIF ! IF ( MOTTYP .EQ. '12HR' .AND. MEANF .EQ. 0 ) THEN ! PRINT *,' # USING 12HR FOR PAST MOTION:' ! UCOMP = UCMP2 ! VCOMP = VCMP2 ! ELSEIF ( MOTTYP .EQ. 'USER' ) THEN ! PRINT *,' # USING USER SPEC. MOTION',USERU(NTY),USERV(NTY) ! UCOMP = USERU(NTY) ! VCOMP = USERV(NTY) ! ELSEIF ( MEANF .EQ. 0 .OR. MEANF .EQ. 12 ) THEN ! PRINT *,' #USING 6HR ONLY FOR PAST MOTION:' ! ELSEIF ( MEANF .EQ. 6 ) THEN ! PRINT *,' # USING 12HR AND NOT 6HR (BAD) FOR PAST MOTION' ! UCOMP = UCMP3 ! VCOMP = VCMP3 ! ELSE ! PRINT *,' # PAST MOTION INVALID...' ! IF ( CRMTNP .GT. 0 ) THEN ! PRINT *,' FAILED CHECK ON 0, 6, 12 HOUR TC LOCATION' ! PRINT *,' NO BOGUS FOR THIS TC' ! NTY = NTY - 1 ! GOTO 55 ! ENDIF ! PRINT *,'... NOT USED ANYWAY - KEEP GOING' ! ENDIF ! PRINT *, 'UCOMP VCOMP ',UCOMP,VCOMP UTY(NTY) = UCOMP VTY(NTY) = VCOMP ! ENDIF ! RELOOP FOR MORE TC'S IF NOT 'END' ! IF ( FUNC .NE. 'END' ) GOTO 55 ! 100 CONTINUE ! IF ( NTY .EQ. 0 ) THEN ! PRINT *,' NO TROPICAL CYCLONES FOR DATE,TIME OVER DOMAIN ' ! STOP ! ELSEIF ( NTY .GT. MXTY ) THEN ! PRINT *,' **TOO MANY TYPHOONS ',NTY,' CAN ONLY DO ',MXTY ! NTY = MXTY ! ENDIF RETURN END !--------------------------------------------------------------------- !********************************************************************* !--------------------------------------------------------------------- SUBROUTINE CALCPS ( PNEW, PS, TV, Z, P2, P3, IL, JL, KL ) DIMENSION PS(IL,JL),TV(IL,JL,KL),Z(IL,JL,KL) DIMENSION PNEW(IL,JL) GRAV=9.81 RGAS=287. ! RECALC MSLP *****DICEY DO 720 I=1,IL DO 720 J=1,JL P1 = PS(I,J) TL = TV(I,J,1) TU = TV(I,J,2) TBAR = TL + (P1 - P2)*(TU - TL)/(P3-P2) ! IF(I.EQ.30.AND.J.EQ.30)PRINT *,' TBAR SRF-1000 ',TBAR PNEW(I,J) = P2*EXP(GRAV*Z(I,J,1)/(RGAS*TBAR)) ! IF(I.EQ.30.AND.J.EQ.30)PRINT *,' NEW MSLP ',PNEW(I,J) 720 CONTINUE BIGD = 0. RMSP = 0. PAV = 0. DO 771 I=1,IL DO 771 J=1,JL ! IF(I.EQ.30.AND.J.EQ.30)PRINT *,' ANL,CALC MSLP ',PS(I,J),PNEW(I,J) PAV = PAV + PNEW(I,J)/(IL*JL) RMSP = RMSP + (PNEW(I,J) - PS(I,J))**2 ADIF = ABS(PNEW(I,J) - PS(I,J)) ! IF(ADIF.GT.5.)PRINT *,' LARGE PDIF ',I,J,PNEW(I,J),PS(I,J) IF(ADIF .GT. BIGD)BIGD = ADIF 771 CONTINUE RMSP = SQRT(RMSP/(IL*JL)) ! PRINT *,' RMS P ERR,BIGD,PAV ',RMSP,BIGD,PAV RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE CENTR2 ( PAI, IM, JM, XCNTR, YCNTR, PCNTR, IJRNGX ) PARAMETER (NSRCHX=15,MAX=(2*NSRCHX+1)*(2*NSRCHX+1)) DIMENSION PAI(IM,JM),IR(MAX),JR(MAX) DATA IJCAL / 0 / ! --------- DEFINE (IR,JR) ---------------------- IF (IJRNGX.GT.NSRCHX) THEN WRITE(6,*)'IJRNGX IS BEYOND THE LIMITS IN SUBR.((CENTR2))' stop ENDIF IF (IJCAL.NE.999) THEN N = 0 MM = 0 6 CONTINUE DO 5 I=-NSRCHX,NSRCHX DO 5 J=-NSRCHX,NSRCHX IF (I**2+J**2.EQ.N) THEN MM = MM+1 IR(MM) = I JR(MM) = J ENDIF 5 CONTINUE N = N+1 IF (N.LE.2*NSRCHX**2) GOTO 6 IJCAL = 999 ENDIF ! --------- NO SEARCH --------------------------- ! ... MISSED IN PAST STEP ! PRINT *,'IS PCNTR < 0 ?' ! write(*,*)'PCNTR=',PCNTR IF(PCNTR.LE.0.) RETURN PCNTR = -1. ! PRINT *,'IS XCNTR/YCNTR OUT OF DOMAIN ?' IF(XCNTR.LT.1 .OR. XCNTR.GT.IM .OR. & ! ... OUT OF DOMAIN YCNTR.LT.1 .OR. YCNTR.GT.JM) then PRINT*, 'CENTRE POINT OUTSIDE DOMAIN (I,J): ', & XCNTR,YCNTR RETURN ENDIF ! ---------- SEARCH CENTER (I) ------------------- ICP = XCNTR JCP = YCNTR DO 10 M=1,(2*IJRNGX+1)*(2*IJRNGX+1) JC = JCP+JR(M) IC = ICP+IR(M) IF (IC-1.LT.1 .OR. JC-1.LT.1 .OR. & ! ... OUT OF DOMAIN IC+1.GT.IM .OR. JC+1.GT.JM) RETURN DO 11 J=JC-1,JC+1 DO 11 I=IC-1,IC+1 IF (PAI(IC,JC).GT.PAI(I,J)) GOTO 10 11 CONTINUE GOTO 20 10 CONTINUE PRINT *,'NO POINT IS A LOCAL MINIMUM!' ! ... OUT OF SEARCHING RANGE RETURN 20 CONTINUE ! ---------- SEARCH CENTER (II) ----------------- P3=PAI(IC+1,JC) P2=PAI(IC-1,JC) P1=PAI(IC,JC) P5=PAI(IC,JC+1) P4=PAI(IC,JC-1) ! WRITE (6,600) P1,P2,P3,P4,P5 ! 600 FORMAT(1H ,'(P1 P2 P3 P4 P5) = ',5(F7.2,1X)) F = P1 C = (P3-P2)/2. D = (P5-P4)/2. ALPHA = (P2+P3+P4+P5)/4.-P1 ! ... CANNOT BE DEFINED ! PRINT *,'IS ALPHA OK?...' IF (ALPHA.EQ.0) RETURN DX =-0.5*C/ALPHA DY =-0.5*D/ALPHA PCNTR = F - ALPHA*(DX*DX+DY*DY) YCNTR = JC+DY XCNTR = IC+DX PRINT *,'FINISHED OK, LOCATED HURRICANE CENTER',PCNTR RETURN END !------------------------------------------------------------ !************************************************************ !------------------------------------------------------------ SUBROUTINE CIDENT ( ITYP, IDENT ) INTEGER CHRINC(4) CHARACTER*12 IDENT CHARACTER*36 CHARS DATA CHRINC/1,1,1,1/ DATA CHARS/'0123456789ABCDEFGHIJKLMNOPQRSTUVWXYZ'/ DATA KT/0/ KT = KT + 1 IF ( ITYP .EQ. 1 ) IDENT = '99999 31' IF ( ITYP .EQ. 2 ) IDENT = '99999 11' IF ( ITYP .EQ. 3 ) IDENT = '99999 12' DO 10 I = 1, 4 CHRINC(I) = CHRINC(I) + 1 IF ( CHRINC(I) .LE. 36 ) GOTO 20 CHRINC(I) = 1 10 CONTINUE 20 CONTINUE DO 30 I=1,4 IDENT(10-I:10-I) = CHARS(CHRINC(I):CHRINC(I)) 30 CONTINUE ! IF(KT.LE.5)PRINT *,' IDENT ',IDENT RETURN END !----------------------------------------------------------------- !***************************************************************** !----------------------------------------------------------------- SUBROUTINE CLY2LL ( ALAT, ALON, ANG, RD, CLAT, CLON ) ! ------------------------------------------------------------+ ! +< OUTPUT > + ! + ALAT : LATITUDE OF THE POINT + ! + ALON : LONGITUDE OF THE POINT + ! + | + ! +< INPUT > + ! + RD : DISTANCE ON GREAT CIRCLE + ! + ANG : IANGLE ON GREAT CIRCLE + ! + SHOWN IN THE SCHEMATIC DIAGRAM + ! + (N) + ! + | + ! + 270 < ANG < 360 | 180 < ANG < 270 + ! + | + ! + --------------- (C) ---------------- + ! + | + ! + 0 < ANG < 90 | 90 < ANG < 180 + ! + | + ! + (S) + ! + CLON : LONGITUDE OF CYLINDRICAL CENTRE + ! + CLAT : LATITUDE OF CYLINDRICAL CENTRE + ! +-----------------------------------------------------------+ REAL X1, Y1, Y1C, Y1S, Y2C, Y2S, ALPHA, RAD, DEG, & Z, ZC, ZS, ANGD, DS R0 = 6371.E+3 RAD = ASIN(1.0)/90.0 DEG = 1.0/RAD ANGD= ANG * RAD Z = RD / R0 X1 = CLON * RAD Y1 = CLAT * RAD Y1C = COS(Y1) Y1S = SIN(Y1) ZC = COS(Z) ZS = SIN(Z) DS =- SIN(ANGD) Y2S = Y1C*DS*ZS + Y1S*ZC ALPHA = ASIN(Y2S) Y2C = COS(ALPHA) ALAT = ALPHA * DEG ALPHA = (ZC-Y1S*Y2S)/(Y1C*Y2C) IF ( ABS ( ALPHA ) .GT. 1.0 ) THEN IF ( ABS ( ALPHA ) .GT. 1.1 ) & WRITE(6,*) ' ** WARNING CLY2LL - ALPHA = ', ALPHA ALPHA = SIGN ( 1.0, ALPHA ) ENDIF ALPHA = ACOS(ALPHA) IF (ANG.LT.90..OR. ANG.GE.270.) ALON = (X1-ALPHA) * DEG IF (ANG.GE.90..AND.ANG.LT.270.) ALON = (X1+ALPHA) * DEG RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE COMBI ( FCMB, FTYM, IM, JM, KM, R1, R2, RTG, JD, W2A ) !U, UM, IM,JM,KM,R1,R2,RTG,NT+10,W2C ! ----TRANSPLANTATION OF MODEL TYPHOON---------------- ! FCMB : COMBINED FIELD FTYM : MODEL TYPHOON FIELD ! MADE (1987 1 5) DIMENSION FCMB(IM,JM, *), FTYM(IM,JM, *), RTG(IM,JM), W2A(IM,JM) DATA ID/0/ IF(ID.NE.JD) THEN ID = JD DO 10 I=1,IM DO 10 J=1,JM IF(RTG(I,J).LT.R1) THEN W2A(I,J) = 1. ELSEIF ( RTG(I,J) .LT. R2 ) THEN W2A(I,J) = (R2 - RTG(I,J))/(R2-R1) ELSE W2A(I,J) = 0. ENDIF 10 CONTINUE ENDIF DO 100 I=1, IM DO 100 J=1, JM C1=W2A(I,J) C2=(1.-W2A(I,J)) DO 101 K=1, KM FCMB(I,J,K) = FTYM(I,J,K) * C1 + FCMB(I,J,K) * C2 101 CONTINUE 100 CONTINUE RETURN END !--------------------------------------------------------------------- !********************************************************************* !--------------------------------------------------------------------- SUBROUTINE DEVSYM ( VAL, CLAT, CLON, RMAX, RTG, VALR, V, & IM, JM, KM, IRM, IRM2, NPRO, DELS, SLON, & XI, XJ, XLAT, XLON,SLAT ) PARAMETER (NA=64) DIMENSION VAL(IM,JM,KM) DIMENSION RTG(IM,JM) DIMENSION VALR(IRM2,KM), V(KM) CHARACTER*4 NPRO PI=3.14159 CALL LLTOXY ( CI,CJ,CLAT,CLON,SLAT,DELS,SLON,XI,XJ,XLAT,XLON ) GRDMAX=RMAX DO 105 I=1,IM IF (RTG(I,1).LT.GRDMAX) GRDMAX=RTG(I,1) IF (RTG(I,JM).LT.GRDMAX) GRDMAX=RTG(I,JM) 105 CONTINUE DO 106 J=1,JM IF (RTG(1,J).LT.GRDMAX) GRDMAX=RTG(1,J) IF (RTG(IM,J).LT.GRDMAX) GRDMAX=RTG(IM,J) 106 CONTINUE GRDMAX=GRDMAX/DELS ! Calculate azimuthal mean DO 101 IR=1,IRM2 RAD=GRDMAX*IR/IRM DO 150 K=1,KM VALR(IR,K)=0 150 CONTINUE DO 102 IA=1,NA ANG=2.*PI*(IA-1)/NA FI=RAD*COS(ANG)+CI FJ=RAD*SIN(ANG)+CJ CALL INTPLZ (V, FI,FJ, VAL,IM,JM,KM) DO 151 K=1,KM VALR(IR,K)=VALR(IR,K)+V(K) 151 CONTINUE 102 CONTINUE DO 152 K=1,KM VALR(IR,K)=VALR(IR,K)/NA 152 CONTINUE ! IF (IR/10.0.EQ.INT(IR/10)) ! > PRINT *,IR,RAD,VALR(IR,1) 101 CONTINUE DO 110 I=1,IM DO 110 J=1,JM IF (RTG(I,J).GE.DELS.AND.RTG(I,J).LE.GRDMAX*DELS) THEN IR1=INT(RTG(I,J)/GRDMAX/DELS*IRM) W2=RTG(I,J)/GRDMAX/DELS*IRM-IR1 W1=1.-W2 DO 111 K=1,KM VAL(I,J,K)=VAL(I,J,K)-(W1*VALR(IR1,K)+W2*VALR(IR1+1,K)) !calculates deviation from azimuthal mean (valr(ir1,k)) 111 CONTINUE ELSE DO 113 K=1,KM VAL(I,J,K)=0 113 CONTINUE ENDIF 110 CONTINUE ! WRITE(6,400) (I,I=CI-5,CI+5) 400 FORMAT(' ',5I6,'|',I5,'|',I5,4I6) DO 201 J=CJ-7,CJ+7 ! WRITE(6,401) J,(VAL(I,J,1),I=CI-5,CI+5) 401 FORMAT(I4,' ',11F6.0) 201 CONTINUE ! PRINT *,' ## dev SYM DONE ' RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE HPROF1 ( DZ, PSEA, IRM, KPM, PC, P3, R0, R3, DR, AXX, & IPP, TCPC, RMW ) !OUTPUT> ! DZ(I,1) : DEVIATION OF POUT(1)-ISOBARIC SURFACE HEIGHT ! PSEA(I) : SEA-LEVEL PRESSURE ! ======= FUJITA'S FORMULA IS MODIFIED TO BE DP/DR=0 AT R=R3 ======= ! ! AA,BB,CC IN " RFACT=CC*FX-AA*X**2-BB " ! RFACT...MODIFICATION FACTOR ! FX......TERM INCLUDED IN FUJITA'S FORMULA ! X.......R/R0 ! ! (1) D(RFACT)/DX=0 AT X=X3 ! (2) RFACT=1 AT X=0 ! (3) RFACT=0 AT X=X3 DIMENSION DZ(IRM,KPM), PSEA(IRM) DELP = P3-PC X3=R3/R0 FX3=1./SQRT(1.+X3**2) FDX3=-X3*(1.+X3**2)**(-1.5) CC=1./(1.-FX3+0.5*X3*FDX3) BB=CC-1. AA=FDX3*CC/(2.*X3) IF ( IPP .EQ. 1 ) THEN DO 200 I = 1, IRM X = DR * ( I - 1 ) / R0 FX = 1.0 / SQRT ( 1.0 + X**2 ) RFACT = CC * FX - AA * X**2 - BB PSEA(I) = P3 - DELP * RFACT DZ(I,1) = -AXX * ALOG ( P3 / PSEA(I) ) ! write(*,*)'HPROF1 dz=',DZ(I,1),' I=',I 200 CONTINUE ELSEIF ( IPP .EQ. 2 ) THEN ! USE HOLLAND'S MODIFIED RANKINE VORTEX DO 100 I = 1, IRM X = DR * ( I - 1 ) / R0 IF ( TCPC .LT. 995.0 ) THEN BTC = ( 995.0 - TCPC ) * 0.6 / 80.0 + 1.2 ELSE BTC = 1.2 ENDIF ATC = RMW ** BTC DKM = X * R0 / 1000.0 IF ( DKM .LT. 1.0 ) DKM = 1.0 PSEA(I) = TCPC + ( P3 - TCPC ) * EXP ( -ATC / ( DKM ** BTC ) ) DZ(I,1) = -AXX * ALOG ( P3 / PSEA(I) ) ! PRINT *,' DKM PSEA P3 TCPC PC ', DKM, PSEA, P3, TCPC, PC 100 CONTINUE ELSE WRITE(6,*) '*** ERROR IN HPROF1 - IPP MUST BE 1 OR 2, NOT ', IPP stop ENDIF RETURN END !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- SUBROUTINE HPROF2 ( DZ, IRM, KPM, LCLT, R1, R2, R3, DR, FWC, & DZSRF ) ! ! DZ(I,LCLT) : DEVIATION OF POUT(LCLT)-ISOBARIC SURFACE HEIGHT ! ======= THREE FORMULAE ARE ASSUMED FOR DZ-PROFILE ======== ! ! AA,BB,CC,DD,EE ! WHERE ! DZ=AA*R**2+BB (R.LE.R1) ! DZ=CC*R+DD (R.LE.R2) ! DZ=EE*(R-R3)**2 (R.GT.R2) ! ! (1) BB=FWC*DZSRF ! (2) AA*R1**2+BB=CC*R1+DD ! (3) 2*AA*R1=CC " DERIVATIVE ! (4) CC*R2+DD=EE*(R2-R3)**2 ! (5) CC=2*EE*(R2-R3) " DERIVATIVE DIMENSION DZ(IRM,KPM) BB=FWC*DZSRF CC=-2.*BB/(R2+R3-R1) DD=-(R2+R3)*0.5*CC EE=CC*0.5/(R2-R3) AA=CC*0.5/R1 DO 210 I=1,IRM R=(I-1)*DR IF(R.LE.R1) THEN DZ(I,LCLT) = AA*R**2 + BB ELSEIF ( R .GT. R1 .AND. R .LE. R2 ) THEN DZ(I,LCLT) = CC*R + DD ELSE DZ(I,LCLT) = EE*(R3-R)**2 ENDIF 210 CONTINUE RETURN END !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- SUBROUTINE INTPL1 ( V, FI, GD, NM, KM, NI ) REAL V(*), GD(NM,*) I =INT(FI) DX=FI-FLOAT(I) DX1 = 1.0 - DX IF (I.LE.1.OR.I.GE.NI-1) GOTO 1 A=0.25*DX*DX1 DO 100 K=1,KM V(K) = (A+DX)*GD(I+1,K) + (A+DX1)*GD(I,K) & - A*(GD(I-1,K) + GD(I+2,K)) 100 CONTINUE RETURN 1 IF (I.LT.1.OR.I.GE.NI) GOTO 2 DO 110 K=1,KM V(K) = DX*GD(I+1,K) + DX1*GD(I,K) 110 CONTINUE RETURN 2 IF (I.LT. 1) NL=1 IF (I.GE.NI) NL=NI DO 120 K=1,KM V(K) = GD(NL,K) 120 CONTINUE RETURN END !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- SUBROUTINE INTPLZ ( V, FI, FJ, GD, NI, NJ, KM ) ! INTERPOLATION BY THE GRID DATA; 12 POINTS ARE USED ! APPROXIMATION FORMULA IS "Z=A*X**2+B*X*Y+C*Y**2+D*X+E*Y+F" ! INNER 4 POINTS ARE ON THE APPROXIMATION SURFACE ! AND OUTER 8 POINTS ARE NEAREST TO THIS IN THE SENSE OF ! LEAST SQUARE MEAN ERROR REAL V(*), GD(NI,NJ,*) J =INT(FJ) I =INT(FI) DX=FI-FLOAT(I) DX1=1.0-DX DY=FJ-FLOAT(J) DY1=1.0-DY IF (I.LE.1.OR.I.GE.NI-1.OR.J.LE.1.OR.J.GE.NJ-1) GOTO 1 A=0.125*DX*DX1 B=0.125*DY*DY1 C=A+B DO 100 K=1,KM V(K) = (C+DX *DY )*GD(I+1,J+1,K) + (C+DX *DY1)*GD(I+1,J ,K) & +(C+DX1*DY )*GD(I ,J+1,K) + (C+DX1*DY1)*GD(I ,J ,K) & -A*(GD(I-1,J ,K)+GD(I+2,J ,K)+GD(I-1,J+1,K)+GD(I+2,J+1,K)) & -B*(GD(I ,J-1,K)+GD(I+1,J-1,K)+GD(I ,J+2,K)+GD(I+1,J+2,K)) 100 CONTINUE RETURN 1 IF (I.LT.1.OR.I.GE.NI.OR.J.LT.1.OR.J.GE.NJ) GOTO 2 DO 110 K=1,KM V(K) = DX *DY *GD(I+1,J+1,K)+DX *DY1*GD(I+1,J ,K) & +DX1*DY *GD(I ,J+1,K) +DX1*DY1*GD(I ,J ,K) 110 CONTINUE RETURN 2 IF (J.LT.1.OR.J.GE.NJ) GOTO 3 IF (I.LT. 1) NL=1 IF (I.GE.NI) NL=NI DO 120 K=1,KM V(K) = DY*GD(NL,J+1,K)+DY1*GD(NL,J,K) 120 CONTINUE RETURN 3 IF (I.LT.1.OR.I.GE.NI) GOTO 4 IF (J.LT. 1) NL=1 IF (J.GE.NJ) NL=NJ DO 130 K=1,KM V(K) = DX*GD(I+1,NL,K)+DX1*GD(I,NL,K) 130 CONTINUE RETURN 4 CONTINUE DO 140 K=1,KM IF (I.LT. 1.AND.J.LT. 1) V(K) = GD( 1, 1,K) IF (I.GE.NI.AND.J.LT. 1) V(K) = GD(NI, 1,K) IF (I.LT. 1.AND.J.GE.NJ) V(K) = GD( 1,NJ,K) IF (I.GE.NI.AND.J.GE.NJ) V(K) = GD(NI,NJ,K) 140 CONTINUE RETURN END !----------------------------------------------------------------------- !*********************************************************************** !----------------------------------------------------------------------- SUBROUTINE INVBAL ( U, V, PZ, SF, VP, FF, F, WORK, IL, JL, KL, & NPRO, DELS, SLON, XI, XJ, XLAT, XLON ) ! ! THIS ROUTINE CALCULATES GEOPOTENTIAL USING ! INVERSE BALANCE EQUATION ! PARAMETER ( PI=3.14159265358979, DTR=PI/180.0, RTD=180.0/PI ) PARAMETER ( G=9.81, LIMIT=75, DELTA=1.E-3 ) ! DIMENSION U(IL,JL,KL),V(IL,JL,KL),PZ(IL,JL,KL),FF(IL,JL) DIMENSION SF(IL,JL),WORK(IL,JL),VP(IL,JL),F(IL,JL) CHARACTER NPRO*4 DATA ISTART/0/ real, dimension(il,jl) :: junk ! !!$ ier=fnom(50,'bal.fst','STD',0) !!$ ier=fstouv(50,'RND') CALL OPTMUM (IL, JL, ALFA, 0.0 ) ! PRINT *,'ALFA = ',ALFA EMBAR=0.0 STLAT=0. DO 30 I=1,IL DO 30 J=1,JL ! TO COMPENSATE FOR INV JMA GRID--v CALL XYTOLL (RLA,RLO,FLOAT(I),FLOAT(JL+1-J),60. & ,NPRO,DELS,SLON,XI,XJ,XLAT,XLON) ZLAT = RLA*DTR EMBAR=EMBAR+ COS(STLAT*DTR)/COS(ZLAT) F(I,J) = 2.*7.292E-5 * SIN(ZLAT) ! IF(I.EQ.1.AND.J.EQ.1)PRINT *,' CORI ',I,J,RLA,F(I,J) 30 CONTINUE EMBAR=EMBAR/(IL*JL) JY2 = 2 IX2 = 2 JY3 = 3 IX3 = 3 JLM = JL-1 ILM = IL-1 ! ! LEVELS LOOP: ! DO 999 KK=1,KL ! ! SETUP FIRST GUESS FOR VP ! DS2=2.*DELS ! DO 88 J=JY2,JLM DO 88 I=IX2,ILM WORK(I,J)=(U(I+1,J,KK)-U(I-1,J,KK))/DS2 + & (V(I,J+1,KK)-V(I,J-1,KK))/DS2 ! WORK(I,J)=WORK(I,J) * 1.E4 88 CONTINUE DO 881 I=1,IL WORK(I,1) = 2.*WORK(I,JY2) - WORK(I,JY3) WORK(I,JL) = 2.*WORK(I,JLM) - WORK(I,JLM-1) 881 CONTINUE DO 882 J=1,JL WORK(1,J) = 2.*WORK(IX2,J) - WORK(IX3,J) WORK(IL,J) = 2.*WORK(ILM,J) - WORK(ILM-1,J) 882 CONTINUE ! CORNERS WORK(1,1)=0.5*(2.*WORK(1,JY2)-WORK(1,JY3)+ & 2.*WORK(IX2,1)-WORK(IX3,1)) WORK(1,JL)=0.5*(2.*WORK(1,JLM)-WORK(1,JLM-1)+ & 2.*WORK(IX2,JL)-WORK(IX3,JL)) WORK(IL,1)=0.5*(2.*WORK(IL,JY2)-WORK(IL,JY3)+ & 2.*WORK(ILM,1)-WORK(ILM-1,1)) WORK(IL,JL)=0.5*(2.*WORK(IL,JLM)-WORK(IL,JLM-1)+ & 2.*WORK(ILM,JL)-WORK(ILM-1,JL)) ! !!$ ier=fstecr(work,junk,-16,50,0,0,0,il,jl,1, & !!$ kk,0,0,'a','W0',' ','X',0,0,0,0,1,.true.) DO 8821 J=1,JL DO 8821 I=1,IL 8821 VP(I,J)=0. !!$ ier=fstecr(VP,junk,-16,50,0,0,0,il,jl,1, & !!$ kk,0,0,'a','VB',' ','X',0,0,0,0,1,.true.) ! SOLVE FOR VP CALL LIEBH ( VP, WORK, ALFA, IL, JL, DELS ) !!$ ier=fstecr(vp,junk,-16,50,0,0,0,il,jl,1, & !!$ kk,0,0,'a','VP',' ','X',0,0,0,0,1,.true.) ! VORTICITY AND STREAM FUNCTION DO 89 J=JY2,JLM DO 89 I=IX2,ILM WORK(I,J)=((V(I+1,J,KK)-V(I-1,J,KK))/DS2 - & (U(I,J+1,KK)-U(I,J-1,KK))/DS2 ) 89 CONTINUE !!$ ier=fstecr(work,junk,-16,50,0,0,0,il,jl,1, & !!$ kk,0,0,'a','W1',' ','X',0,0,0,0,1,.true.) DO 891 I=1,IL WORK(I,1) = 2.*WORK(I,JY2) - WORK(I,JY3) WORK(I,JL) = 2.*WORK(I,JLM) - WORK(I,JLM-1) 891 CONTINUE DO 892 J=1,JL WORK(1,J) = 2.*WORK(IX2,J) - WORK(IX3,J) WORK(IL,J) = 2.*WORK(ILM,J) - WORK(ILM-1,J) 892 CONTINUE ! CORNERS WORK(1,1)=0.5*(2.*WORK(1,JY2)-WORK(1,JY3)+ & 2.*WORK(IX2,1)-WORK(IX3,1)) WORK(1,JL)=0.5*(2.*WORK(1,JLM)-WORK(1,JLM-1)+ & 2.*WORK(IX2,JL)-WORK(IX3,JL)) WORK(IL,1)=0.5*(2.*WORK(IL,JY2)-WORK(IL,JY3)+ & 2.*WORK(ILM,1)-WORK(ILM-1,1)) WORK(IL,JL)=0.5*(2.*WORK(IL,JLM)-WORK(IL,JLM-1)+ & 2.*WORK(ILM,JL)-WORK(ILM-1,JL)) !!$ ier=fstecr(work,junk,-16,50,0,0,0,il,jl,1, & !!$ kk,0,0,'a','W2',' ','X',0,0,0,0,1,.true.) ! SET UP FIRST GUESS FOR STREAM FUNCTION !!Z following is the original scheme !!Z DO 893 J = 1, JL !!Z DO 893 I = 1, IL !!Z CORPM=0.5E-05 !!Z SF(I,J)=PZ(I,J,KK)/(CORPM) !!Z THE UNITS OF SF FIRST GUESS ARE NOW M**2 SEC-1.....I HOPE !!Z SF = ORDER(10**12),,VP=ORDER(10**11)??????? !!Z IN ROUTINE LIEBH, UNITS REQUIRED ARE CM AND SEC<<< (XLAT,XLON) * ! STANDARD POINT * !********************************************************************** DIMENSION RLAT(IM,*), RLON(IM,*) CHARACTER * 4 NPROJC INTEGER ICHOICE ! RADIUS OF EARTH" A1=6371.E3 PI=3.14159 POI=PI /180. RPOI=180./ PI !-----------------------------------POLAR STEREO----------------------- IF(NPROJC(1:3).EQ.'PSN') THEN ! ! MAP FACTOR=1 ! STANDARD LATITUDE if (slat < -98.) slat=60. ! default map scale == 1. SLAT1=SLAT*POI SLON1=SLON*POI XLAT1=XLAT*POI XLON1=XLON*POI RL0=A1*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PI-XLAT1)) X0=RL0*SIN(XLON1-SLON1) Y0=RL0*COS(XLON1-SLON1) DO 10 I=1,IM DO 10 J=1,JM XP=X0+(I-XI)*DELS YP=Y0+(J-XJ)*DELS RLAT(I,J)=90.-RPOI*2.*ATAN(SQRT(XP*XP+YP*YP)/ & (A1*(1.+SIN(SLAT1)))) RLON(I,J)=SLON+RPOI*ATAN2(XP,YP) 10 CONTINUE ENDIF IF(NPROJC(1:3).EQ.'PSS') THEN ! MAP FACTOR=1 ! STANDARD LATITUDE if (slat < -98) slat = -60. ! default map scale == 1. SLAT1=-SLAT*POI SLON1= SLON*POI XLAT1=-XLAT*POI XLON1= XLON*POI RL0=A1*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PI-XLAT1)) X0=RL0*SIN(XLON1-SLON1) Y0=RL0*COS(XLON1-SLON1) DO 15 I=1,IM DO 15 J=1,JM XP=X0+(I-XI)*DELS YP=Y0+(XJ-J)*DELS RLAT(I,J)=-90.+RPOI*2.*ATAN(SQRT(XP*XP+YP*YP)/(A1*(1.+SIN(SLAT1 & )))) RLON(I,J)=SLON+RPOI*ATAN2(XP,YP) 15 CONTINUE ENDIF !-----------------------------------MERCATOR--------------------------- IF(NPROJC(1:3).EQ.'MER') THEN ! MAP FACTOR=1 ! STANDARD LATITUDE NOTE: STANDARD LAT IS SUBJECT TO GRID if (slat < -98.) slat = 0. ! default map scale == 1. SLAT1=SLAT*POI SLON1=SLON*POI XLAT1=XLAT*POI XLON1=XLON*POI AC=A1*COS(SLAT1) DO 20 I=1,IM DO 20 J=1,JM RLON(I,J)=XLON+RPOI*(I-XI)*DELS/AC T1=EXP((XJ-J)*DELS/AC)*(1.+SIN(XLAT1))/COS(XLAT1) RLAT(I,J)=RPOI*ASIN((T1*T1-1.)/(T1*T1+1.)) 20 CONTINUE ENDIF !-----------------------------------LAMBERT---------------------------- IF(NPROJC(1:3).EQ.'LMN') THEN ! MAP FACTOR=1 ! STANDARD LATITUDE SLATA=30. SLATB=60. PI4=PI*0.25 SLATA1=SLATA*POI SLATB1=SLATB*POI SLAT1=PI4-SLATA*POI*0.5 SLAT2=PI4-SLATB*POI*0.5 SLON1=SLON*POI XLAT1=XLAT*POI XLON1=XLON*POI CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) RCK=1./CK ACN=A1*COS(SLATA1)*RCK RL0=ACN*(TAN(PI4-XLAT1*0.5)/TAN(SLAT1))**CK ! ---------------------------------------- XSLON = XLON - SLON XSLON = MOD (XSLON+900., 360.) - 180. XSLON = XSLON * POI * CK X0=RL0*SIN(XSLON) Y0=RL0*COS(XSLON) ! ---------------------------------------- ! ---- X0=RL0*SIN((XLON1-SLON1)*CK) ! ---- Y0=RL0*COS((XLON1-SLON1)*CK) DO 30 I=1,IM DO 30 J=1,JM XP=X0+(I-XI)*DELS YP=Y0+(J-XJ)*DELS RLAT(I,J)=90.-2.*RPOI*ATAN((SQRT(XP*XP+YP*YP)/ACN)**RCK & *TAN(SLAT1)) RLON(I,J)=SLON+RCK*RPOI*ATAN2(XP,YP) ! ---------------------------------------- RRLON = RLON(I,J) RRLON = RRLON - XLON RRLON = MOD (RRLON+900., 360.) - 180. RLON(I,J) = RRLON + XLON ! ---------------------------------------- 30 CONTINUE ENDIF IF(NPROJC(1:3).EQ.'LMS') THEN ! MAP FACTOR=1 ! STANDARD LATITUDE SLATA=-30. SLATB=-60. PI4=PI*0.25 SLATA1=-SLATA*POI SLATB1=-SLATB*POI SLAT1=PI4+SLATA*POI*0.5 SLAT2=PI4+SLATB*POI*0.5 SLON1=SLON*POI XLAT1=-XLAT*POI XLON1=XLON*POI CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) RCK=1./CK ACN=A1*COS(SLATA1)*RCK RL0=ACN*(TAN(PI4-XLAT1*0.5)/TAN(SLAT1))**CK ! ---------------------------------------- XSLON = XLON - SLON XSLON = MOD (XSLON+900., 360.) - 180. XSLON = XSLON * POI * CK X0=RL0*SIN(XSLON) Y0=RL0*COS(XSLON) ! ---------------------------------------- ! ---- X0=RL0*SIN((XLON1-SLON1)*CK) ! ---- Y0=RL0*COS((XLON1-SLON1)*CK) DO 35 I=1,IM DO 35 J=1,JM XP=X0+(I-XI)*DELS YP=Y0+(XJ-J)*DELS RLAT(I,J)=-90.+2.*RPOI*ATAN((SQRT(XP*XP+YP*YP)/ACN)**RCK & *TAN(SLAT1)) RLON(I,J)=SLON+RCK*RPOI*ATAN2(XP,YP) ! ---------------------------------------- RRLON = RLON(I,J) RRLON = RRLON - XLON RRLON = MOD (RRLON+900., 360.) - 180. RLON(I,J) = RRLON + XLON ! ---------------------------------------- 35 CONTINUE ENDIF !-----------------------------------LATITUDE LONGITUDE----------------- IF(NPROJC(1:2).EQ.'LL') THEN DO 50 J=1,JM DO 50 I=1,IM ! RLON(I,J) = (I-1)*(DELS/111000.) ! RLAT(I,J) = 90.-(J-1)*(DELS/111000.) RLON(I,J) = XLON + (I-XI)*(DELS/111000.) RLAT(I,J) = XLAT - (J-XJ)*(DELS/111000.) 50 CONTINUE ENDIF !--------------------------------------------------------------------- DO 90 J = 1, JM DO 90 I = 1, IM IF ( RLAT(I,J) .GT. 90.0 ) RLAT(I,J) = 90.0 IF ( RLAT(I,J) .LT. -90.0 ) RLAT(I,J) = -90.0 90 CONTINUE DO 95 KK=1,2 DO 95 J = 1, JM DO 95 I = 1, IM IF ( RLON(I,J) .LT. 0.0 ) RLON(I,J) = RLON(I,J) + 360.0 IF ( RLON(I,J) .GT. 360.0 ) RLON(I,J) = RLON(I,J) - 360.0 95 CONTINUE RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE LIEBH_ORI ( SAL, FORC, ALFA, IL, JL, DELS ) ! THIS S/R SOLVES THE STREAM FUNCTION - VORTICITY EQUATION PARAMETER ( RESTOL=1.E-5) REAL SAL(IL,JL), FORC(IL,JL) DSSQ = DELS ** 2 IP=0 60 CONTINUE KP=0 DO 20 J = 2, JL-1 DO 20 I = 2, IL-1 dsol = -dssq*forc(i,j) - 4.*sal(i,j) + & (sal(i+1,j)+sal(i-1,j)+sal(i,j+1)+sal(i,j-1)) sal(i,j) = sal(i,j) + alfa*dsol IF ( ABS(DSOL) .GT. ABS(RESTOL*SAL(I,J)) ) KP = 1 20 CONTINUE IP = IP + 1 IF ( KP .EQ. 1 .AND. IP .LT. 151 ) GOTO 60 RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE LIEBH ( SAL, FORC, ALFA, IL, JL, DELS ) ! THIS S/R SOLVES THE STREAM FUNCTION - VORTICITY EQUATION PARAMETER ( RESTOL=2.) REAL SAL(IL,JL), FORC(IL,JL) real meanFOR,sdSAL,dsolMax DSSQ = DELS ** 2 IP=0 60 CONTINUE meanFOR = sum(forc)/float(il*jl) dsolMax = 0. KP=0 DO 20 J = 2, JL-1 DO 20 I = 2, IL-1 dsol = dssq*forc(i,j) + 4.*sal(i,j) - & (sal(i+1,j)+sal(i-1,j)+sal(i,j+1)+sal(i,j-1)) sal(i,j) = sal(i,j) - alfa*dsol if (abs(dsol) > dsolMax) dsolMax = dsol ! IF ( ABS(DSOL) .GT. ABS(RESTOL*SAL(I,J)) ) KP = 1 20 CONTINUE if (abs(dsolMax) > abs(restol*dssq*meanFOR)) kp=1 IP = IP + 1 IF ( KP .EQ. 1 .AND. IP .LT. 151 ) GOTO 60 RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE LLTOXY ( GI, GJ, RLAT, RLON,SLAT, & DELS, SLON, XI, XJ, XLAT, XLON ) !********************************************************************** ! GRID POINT OF LATITUDE AND LONGITUDE * ! * !OUT-PUT (GI,GJ) GRID POINT * ! * !IN-PUT RLAT LATITUDE (DEGREE) -90.< < 90. * ! RLON LONGITUDE (DEGREE) 0.< <360. * ! 'PSN ' POLAR STEREO SLAT= 60 N * ! 'PSS ' POLAR STEREO SLAT=-60 N * ! 'MER ' MERCATOR SLAT= 0. N * ! 'LMN ' LAMBERT SLAT= 30 N , 60 N * ! 'LMS ' LAMBERT SLAT=-30 N ,-60 N * ! DELS GRID PARAM ( M) * ! Y-COORDINATE 0.< <360.* * ! SLON STANDARD LONGITUDE* ! * ! (XI,XJ) <----------> (XLAT,XLON) * ! STANDARD POINT * !********************************************************************** ! RADIUS OF EARTH" RA=6371.E3 PAI=3.14159 RAD=PAI /180. !-----------------------------------MERCATOR--------------------------- ! STANDARD LATITUDE ;MAP FACTOR=1" ! NOTE: SET SLAT ACCORDING TO SPECIFIC GRID ! SLAT=42.3687 SLAT1=SLAT*RAD SLON1=SLON*RAD XLAT1=XLAT*RAD XLON1=XLON*RAD AC=RA*COS(SLAT1) X0=0. Y0=AC*LOG((1.+SIN(XLAT1))/COS(XLAT1)) RLAT1=RLAT*RAD RLON1=RLON*RAD ! ---------------------------- DI = RLON - XLON DI = MOD (DI+900., 360.) - 180. DI = AC*(DI*RAD) ! ---------------------------- ! ---- DI=AC*(RLON1-XLON1) DJ=AC*LOG((1.+SIN(RLAT1))/COS(RLAT1)) GI=(DI-X0)/DELS+XI GJ=(Y0-DJ)/DELS+XJ ! ENDIF !---------------------------------------------------------------------- RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE LOCCEN ( U, V, KUSE, GLAT1, GLON1, GLAT2, GLON2, & MAXDIS, CLAT, CLON, IM, JM, KM, NPRO, DELS, & SLON, XI, XJ, XLAT, XLON, SLAT ) ! THIS ROUTINE LOCATES THE TYPHOON CENTRE USING THE ! APPROACH THAT THE CENTRE OF THE TYPHOON IS THE POINT WHERE ! MOST OF THE SURROUNDING WINDS ARE GOING CLOCKWISE (SH) ! AROUND THAT POINT, AND THERE IS MINIMAL RADIAL WIND COMPONENT ! PARAMETERS: ! U,V =WIND COMPONENTS (IM,JM,KM) ! KUSE =LEVEL TO LOCATE WIND CENTRE ON ! GLAT1,GLON1-GLAT2,GLON2=BOX TO LOOK FOR CENTRE IN ! MAXDIS =MAXIMUM DIFF FROM BOX CENTRE (DEGREES) ! DIM IM,JM,KM, PROJ: NPRO,DELS,SLON,XI,XJ,XLAT,XLON ! OUTPUT: ! CLAT,CLON =PROPOSED TYPHOON CENTRE PARAMETER (IRM=12,NA=64) PARAMETER (IBOX=5,JBOX=5,NBOX=IBOX*JBOX,MAXLEV=25) PARAMETER ( PI=3.14159265358979 ) DIMENSION U(IM,JM,KM),V(IM,JM,KM) DIMENSION BI(NBOX),BJ(NBOX),CIRC(NBOX),NB(NBOX) DIMENSION RD(IRM),VELU(30),VELV(30) REAL MAXDIS CHARACTER*4 NPRO DATA RD/.5,1,1.5,2,2.5,3,3.5,4,4.5,5,5.5,6/ ! PRINT *,'LOCCEN####FIND CENTRE ON PRESS LEVEL NO. (K=) ',KUSE ! >>> CHECK INPUT PARAMETERS IF (GLAT1*GLAT2.LT.0) THEN ! PRINT *,'SEARCH AREA OVERLAPS 2 HEMISPHERES SO USE S/R CENTR2 ', ! > 'TO FIND TC CENTRE' CLON = -9999 RETURN ELSEIF ( GLAT1+GLAT2 .LT. 0 ) THEN PRINT *,'SEARCHING SOUTHERN HEMISPHERE...' IHEMI=1 ELSE PRINT *,'SEARCHING NORTHERN HEMISPHERE...' IHEMI=-1 ENDIF ! PRINT *,'LOCCEN##START - ' ! PRINT *,'IM,JM ',IM,JM,' DELS ',DELS ! PRINT *,'NPRO ',NPRO,' SLON ',SLON ! PRINT *,'XI,XJ,XLAT,XLON ',XI,XJ,XLAT,XLON ! PRINT *,'BOX L/L ',GLAT1,GLON1,GLAT2,GLON2 ! CALCULATE CENTRE USING TANGENTIAL VELOCITY GLATM=(GLAT1+GLAT2)*.5 GLONM=(GLON1+GLON2)*.5 DLAT=ABS(GLAT2-GLAT1) DLON=ABS(GLON2-GLON1) DO 5 IL=1,MAXLEV ! PRINT *,'*** LEVEL ',IL,' ***' N=0 DO 100 IB=1,IBOX DO 100 JB=1,JBOX N=N+1 GLAT=GLATM+((IB-.5-IBOX*.5)/IBOX)*DLAT GLON=GLONM+((JB-.5-JBOX*.5)/JBOX)*DLON ! PRINT *,' LAT LON BOX ',N,' CENTRE ',GLAT,GLON CALL RLTLN(BI(N),BJ(N),GLAT,GLON,SLAT, & NPRO,DELS,SLON,XI,XJ,XLAT,XLON) ! PRINT *,' GI GJ BOX ',N,' CENTRE ',BI(N),BJ(N) 100 CONTINUE DO 6 IB=1,NBOX CIRC(IB)=0. DO 10 IR=1,IRM DO 20 JA=1,NA ANG=2.*PI*(JA-1)/NA FI=BI(IB)+RD(IR)*COS(ANG) FJ=BJ(IB)-RD(IR)*SIN(ANG) CALL INTPLZ(VELU, FI,FJ, U, IM,JM,KM) CALL INTPLZ(VELV, FI,FJ, V, IM,JM,KM) VT=VELU(KUSE)*SIN(ANG)-VELV(KUSE)*COS(ANG) VR=VELU(KUSE)*COS(ANG)+VELV(KUSE)*SIN(ANG) CIRC(IB)=CIRC(IB)+IHEMI*ABS(VT)*VT/(VR*VR+VT*VT+.001) 20 CONTINUE 10 CONTINUE ! WRITE(6,1000) IB,BI(IB),BJ(IB),CIRC(IB) 1000 FORMAT(' BOX ',I2,' GI,GJ ',F7.3,F7.3,' *CIRC* ',F17.10) 6 CONTINUE DO 101 IB=1,NBOX NB(IB)=IB 101 CONTINUE DO 7 I1=1,NBOX-1 DO 7 I2=I1+1,NBOX IF (CIRC(NB(I2)).GT.CIRC(NB(I1))) THEN IB=NB(I2) NB(I2)=NB(I1) NB(I1)=IB ENDIF 7 CONTINUE ! PRINT *,' BEST GUESS BOX ',NB(1),BI(NB(1)),BJ(NB(1)) CALL XYTOLL (GLAT,GLON,BI(NB(1)),BJ(NB(1)),SLAT,NPRO,DELS,SLON, & XI,XJ,XLAT,XLON) ! PRINT *,' T.C. AT L/L ',GLAT,GLON GLATM=GLAT GLONM=GLON DLAT=DLAT*.8 DLON=DLON*.8 5 CONTINUE CIRCI=BI(NB(1)) CIRCJ=BJ(NB(2)) ! PRINT *,' TANGENTIAL WIND GUESS GI,GJ ',CIRCI,CIRCJ CALL XYTOLL (CLAT,CLON,CIRCI,CIRCJ,SLAT,NPRO,DELS,SLON, & XI,XJ,XLAT,XLON) PRINT *,' => HURRICANE AT LAT/LONG: ',CLAT,CLON ! >>>>> CHECK THAT LOCATED TYPHOON WITHIN MAXDIS OF L/L BOX CENTRE IF ((ABS(CLAT-0.5*(GLAT1+GLAT2)).GE.MAXDIS).OR. & (ABS(CLON-0.5*(GLON1+GLON2)).GE.MAXDIS)) THEN PRINT *,'##CENTRE LOCATED TO FAR AWAY FROM LAT/LONG BOX CENTRE' PRINT *,'## MAXIMUM DISTANCE (DEGREES) = ',MAXDIS PRINT *,'##BOX CENTRE LAT/LONG: ',0.5*(GLAT1+GLAT2), & 0.5*(GLON1+GLON2) ! PRINT *,'##SETTING CYCLONE LONG TO -999' CLON=-999 RETURN ENDIF DO 200 J=-3,3 ! WRITE(6,201) 'U',J, ! > (U(I+INT(BI(NB(1))),J+INT(BJ(NB(1))),1),I=-3,3) ! WRITE(6,201) 'V',J, ! > (V(I+INT(BI(NB(1))),J+INT(BJ(NB(1))),1),I=-3,3) 201 FORMAT(' ',A1,I3,3F6.1,'|',F6.1,'|',3F6.1) 200 CONTINUE ! PRINT *,'LOCCEN####FINISHED O.K. ' RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE MAPFCT ( SM, IM, JM, RLAT, NPROJC, SLAT ) !********************************************************************** ! MAP FACTOR * ! * !OUT-PUT SM(IM,JM) MAP FACTOR * ! * !IN-PUT RLAT(IM,JM) LATITUDE (DEGREE) 0.< <360. * ! | NPROJ! MAP PROJECTION * ! | 'PSN' POLAR STEREO SLAT= 60 N * ! | 'PSS' POLAR STEREO SLAT=-60 N * ! | 'MER' MERCATOR SLAT= 0. N * ! | 'LMN' LAMBERT SLAT= 30 N , 60 N * ! | 'LMS' LAMBERT SLAT=-30 N ,-60 N * ! | 'LL' LAT/LONG SLAT= 0. N ! * !********************************************************************** DIMENSION RLAT(IM,*),SM(IM,*) CHARACTER*4 NPROJC PI=3.14159 POI=PI /180. !-----------------------------------LAT/LONG--------------------------- IF(NPROJC(1:2).EQ.'LL') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=0. SLAT1=SLAT*POI DO I=1,IM DO J=1,JM ! SM(I,J)=COS(RLAT(I,J)*POI) !POSSIBLE? SM(I,J)=1. ENDDO ENDDO ENDIF !-----------------------------------POLAR STEREO----------------------- IF(NPROJC(1:3).EQ.'PSN') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=60. SLAT1=SLAT*POI DO 10 I=1,IM DO 10 J=1,JM SM(I,J)=(1.+SIN(SLAT1))/(1.+SIN(RLAT(I,J)*POI)) 10 CONTINUE ENDIF IF(NPROJC(1:3).EQ.'PSS') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=-60. SLAT1=SLAT*POI DO 15 I=1,IM DO 15 J=1,JM SM(I,J)=(1.-SIN(SLAT1))/(1.-SIN(RLAT(I,J)*POI)) 15 CONTINUE ENDIF !-----------------------------------MERCATOR--------------------------- IF(NPROJC(1:3).EQ.'MER') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" ! NOTE: AGAIN SET SLAT ! SLAT=0. ! SLAT=42.3687 SLAT1=SLAT*POI DO 20 I=1,IM DO 20 J=1,JM SM(I,J)=COS(SLAT1)/COS(RLAT(I,J)*POI) 20 CONTINUE ENDIF !-----------------------------------LAMBERT---------------------------- IF(NPROJC(1:3).EQ.'LMN') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLATA=30. SLATB=60. PI4=PI*0.25 SLATA1=SLATA*POI SLATB1=SLATB*POI SLAT1=PI4-SLATA*POI*0.5 SLAT2=PI4-SLATB*POI*0.5 CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) CNS=COS(SLATA1)/(TAN(SLAT1))**CK DO 30 I=1,IM DO 30 J=1,JM SM(I,J)=CNS*(TAN(PI4-RLAT(I,J)*POI*0.5))**CK/COS(RLAT(I,J)*POI) 30 CONTINUE ENDIF IF(NPROJC(1:3).EQ.'LMS') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLATA=-30. SLATB=-60. PI4=PI*0.25 SLATA1=-SLATA*POI SLATB1=-SLATB*POI SLAT1=PI4-SLATA1*0.5 SLAT2=PI4-SLATB1*0.5 CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) CNS=COS(SLATA1)/(TAN(SLAT1))**CK DO 35 I=1,IM DO 35 J=1,JM SM(I,J)=CNS*(TAN(PI4+RLAT(I,J)*POI*0.5))**CK/COS(RLAT(I,J)*POI) 35 CONTINUE ENDIF !---------------------------------------------------------------------- RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE OUTY ( PSEA, PAI, T, WV, U, V, IM, JM, KM, IC, JC, & INCR, COMENT ) DIMENSION PSEA(IM,JM), PAI(IM,JM), T(IM,JM,KM) & , WV(IM,JM,KM), U(IM,JM,KM), V(IM,JM,KM) & , LABI(15),LABJ(15) CHARACTER*(*) COMENT IF (IC.LE.1.OR.IC.GE.IM.OR.JC.LE.1.OR.JC.GE.JM) RETURN DO 10 N=1,15 LABI(N) = IC + (N-8)*INCR LABJ(N) = JC + (N-8)*INCR 10 CONTINUE ! WRITE(6,*) '+++++++++++++++++++++++ ' ! WRITE(6,*) '+++ CROSS + SECTION +++ ',COMENT ! WRITE(6,*) '+++++++++++++++++++++++ ' IF ( LABI(1).LT.1 .OR. LABI(15).GT.IM ) GOTO 500 ! WRITE(6,600) LABI ! WRITE(6,610) (PSEA(I,JC),I=LABI(1),LABI(15)) ! WRITE(6,610) (PAI (I,JC),I=LABI(1),LABI(15)) ! WRITE(6,610) ((T (I,JC,K),I=LABI(1),LABI(15)),K=1,KM) ! WRITE(6,610) ((WV (I,JC,K),I=LABI(1),LABI(15)),K=1,KM) ! WRITE(6,610) ((U (I,JC,K),I=LABI(1),LABI(15)),K=1,KM) ! WRITE(6,610) ((V (I,JC,K),I=LABI(1),LABI(15)),K=1,KM) 500 CONTINUE IF ( LABJ(1).LT.1 .OR. LABJ(15).GT.JM ) GOTO 1000 ! WRITE(6,605) LABJ ! WRITE(6,610) (PSEA(IC,J),J=LABJ(1),LABJ(15)) ! WRITE(6,610) (PAI (IC,J),J=LABJ(1),LABJ(15)) ! WRITE(6,610) ((T (IC,J,K),J=LABJ(1),LABJ(15)),K=1,KM) ! WRITE(6,610) ((WV (IC,J,K),J=LABJ(1),LABJ(15)),K=1,KM) ! WRITE(6,610) ((U (IC,J,K),J=LABJ(1),LABJ(15)),K=1,KM) ! WRITE(6,610) ((V (IC,J,K),J=LABJ(1),LABJ(15)),K=1,KM) 1000 CONTINUE 600 FORMAT(1H ,'I =',7I7,' |',I7,' | ',7I7) 605 FORMAT(1H ,'J =',7I7,' |',I7,' | ',7I7) 610 FORMAT((1H ,3X,7F7.1,' |',F7.1,' | ',7F7.1)) RETURN END !-------------------------------------------------------------- !************************************************************** !-------------------------------------------------------------- SUBROUTINE PHICAL ( PHI, T, PHIS, PSEA, SPLA, MAX, KM ) ! CALL PHICAL(W3B, T,PHIS,W2G,SPLA, IM*JM,KM) DIMENSION & PHI(MAX,KM), T(MAX,KM), PHIS(MAX), PSEA(MAX), SPLA(KM) REAL * 8 DFKAPA ! --------------------------------------------------------- CP = 1004.6 R = 287.04 DFKAPA = 0.2857D0 ! --------------------------------------------------------- DO 100 K=1,KM IF (K .EQ. 1) THEN DO 210 M=1,MAX PHI(M,1) = PHIS(M) & - (R * LOG(SPLA(1)/PSEA(M))) * T(M,1) 210 CONTINUE ELSE AA = 0.5*CP * ((SPLA(K-1)/SPLA(K ))**DFKAPA - 1.) BB = - 0.5*CP * ((SPLA(K )/SPLA(K-1))**DFKAPA - 1.) DO 220 M=1,MAX PHI(M,K) = PHI(M,K-1) + AA*T(M,K ) & + BB*T(M,K-1) 220 CONTINUE ENDIF 100 CONTINUE RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE PSEA2D ( DATA2, IM, JM, RIJ, DATA1, IRM, DR ) ! ----------------------------------------------------------------- ! PROJECTION OF SCALAR DATA FROM AXISYMMETRI! 1-D FIELD ! INTO 2-D HORIZONTAL FIELD ! ----------------------------------------------------------------- DIMENSION DATA2(IM,JM), DATA1(IRM), RIJ(IM,JM) DRIN=1./DR DO 100 I=1,IM DO 100 J=1,JM IR = INT(RIJ(I,J)*DRIN) + 1 IF(IR.GE.IRM) THEN DATA2(I,J) = DATA1(IRM) GO TO 100 ENDIF DDR = RIJ(I,J) - DR*(IR-1) DATA2(I,J) = DATA1(IR) + (DATA1(IR+1)-DATA1(IR))*DDR*DRIN 100 CONTINUE RETURN END !--------------------------------------------------------------------- !********************************************************************* !--------------------------------------------------------------------- SUBROUTINE QQRELH ( WV, RH, T, SPLA, IM, JM, KM, WVFCT, IQQ ) ! IQQ = 'QTOH' : (SPECIFIC) TO (RELATIVE) HUMIDITY CONVERSION ! = 'HTOQ' : (RELATIVE) TO (SPECIFIC) HUMIDITY CONVERSION ! WVFCT = 'G/KG' OR 'G/G ' ! RH -----------> % DIMENSION WV(IM*JM,KM), T(IM*JM,KM), RH(IM*JM,KM), SPLA(KM) CHARACTER*4 IQQ,WVFCT EFCT = 0.622 IF (WVFCT(3:3) .EQ. 'K') EFCT = 1000.*EFCT AEFCT = 1./EFCT ! ------------------------- IF (IQQ .EQ. 'HTOQ') THEN ! ------------------------- DO 10 K=1,KM DO 10 I=1,IM*JM X = RH(I,K) * 0.01 TC = T(I,K)-273.2 ES = 6.11*10.**(7.5*TC/(237.3+TC)) E = ES*X WV(I,K) = EFCT*E/SPLA(K) 10 CONTINUE ! ------------------------- ELSE ! ------------------------- DO 20 K=1,KM DO 20 I=1,IM*JM TC = T(I,K)-273.15 ES = 6.11*10.**(7.5*TC/(237.3+TC)) E = AEFCT*WV(I,K)*SPLA(K) RH(I,K) = E/ES * 100. 20 CONTINUE ! ------------------------- ENDIF ! ------------------------- RETURN END !---------------------------------------------------------------- !**************************************************************** !---------------------------------------------------------------- SUBROUTINE RCONS ( R0, DP, PC, R15, F0, PE, TS, DELS ) ! DETERMINE THE CHARACTERISTIC RADIUS OF TYPHOON ! BASED ON FUJITA'S FORMULA ! (( P(R)=PE-DP/(1+(R/R0)**2.0)**0.5 )) ! UNDER GRADIENT WIND BALANCE ! ! F0 : CORIOLIS PARAMETER AT CENTER (/SEC) ! PE : ENVIRONMENTAL SEA-LEVEL PRESSURE (HPA) ! TS : SEA SURF. TEMPERATURE (K) ! DELS : GRID SIZE (M) ! ! PC : CENTRAL PRESSURE (HPA) ! R15 : 15(M/S) RADIUS (M) ! ! R0 : CHARACTERISTIC RADIUS (M) ! DP : = PE-PC ! ===== (R15,R0,DP) ARE NOT INDEPENDENT ======== PARAMETER ( ITEMAX=15, RGAS=287.05 ) ! SEA-LEVEL AIR DENSITY ROH = PE/RGAS/TS ! <<<1-ST>>>.....FUNDAMENTAL CHECK OF "DP" DP = PE - PC ! write(*,*)'DP=',DP ! write(*,*)'PC=',PC,' PE=',PE IF(DP.LT.0.)THEN WRITE(6,*) 'ERR CENT. PRES. =',PC,'ENVIRON. PRES. =',PE stop ENDIF ! <<<2-ND>>>.....FUNDAMENTAL CHECK OF "R15" ! ARBITRARY PARAM. R0MIN = DELS ! ARBITRARY PARAM. R15MIN = SQRT(2.)*R0MIN IF(R15.LT.R15MIN) THEN WRITE(6,*) '*** R15 LESS THAN R15MIN IN SUBR.((RCONS)) ***' WRITE(6,*) ' (R15,R15MIN) =',R15,R15MIN R15 = R15MIN ENDIF ! <<<3-RD>>>.....MUTUAL(LOGICAL) CHECK OF "R15 & DP" ITE = 0 VT = F0/ABS(F0)*15.0 GRF = (F0*VT+VT*VT/R15)*ROH ! MINIMUM OF DP WHICH CAN CAUSE THE WIND OF 15(M/S) AT R15 ! DPMIN <------ R0=R15/SQRT(2.) <------ D(DP)/D(R0)=0 DPMIN = 3.**1.5 * R15/2. * GRF IF(DP.LT.DPMIN) THEN WRITE(6,*) '*** DP LESS THAN DPMIN IN SUBR.((RCONS)) ***' WRITE(6,*) ' (DP,DPMIN) =',DP,DPMIN ! (I) MODIFY R15 ! AA = 3.**1.5/2.*ROH ! XR = (DP-AA*VT*VT)/(AA*F0*VT) ! IF(XR.GT.R15MIN) THEN ! R15 = XR ! ELSE ! DP = DPMIN ! ENDIF ! (II) MODIFY DP DP = DPMIN write(*,*)'DP2=',DP R0 = R15/SQRT(2.) GO TO 500 ENDIF ! CALCULATE R0 BY WEGSTEIN'S ITERATION XR0 = 0. ! FIRST ESTIMATE YR0 = GRF*(R15*R15+XR0*XR0)**1.5/DP/R15 XR1 = YR0 ! SECOND ESTIMATE YR1 = GRF*(R15*R15+XR1*XR1)**1.5/DP/R15 20 CONTINUE ITE = ITE + 1 WEG = (YR1-YR0)/(XR1-XR0) XR = (YR1-XR1*WEG)/(1.-WEG) IF (ABS(XR-XR1) .LT. 100.) GO TO 30 XR0 = XR1 XR1 = XR YR0 = YR1 YR1 = GRF*(R15*R15+XR1*XR1)**1.5/DP/R15 IF (ITE .LT. ITEMAX) GO TO 20 30 CONTINUE R0 = XR 500 CONTINUE ! WRITE(6,600) R0,ITE ! 600 FORMAT(1H ,' CHARACTERISTIC RADIUS (R0) =',F10.1 ! > ,' ITE =',I3 ) ! <<<4-TH>>>.....FUNDAMENTAL CHECK OF "R0" ! MINIMUM OF R0 WHICH CAN RESOLVE INNER STRUCTURE IF(R0.LT.R0MIN) THEN WRITE(6,*) '*** R0 LESS THAN R0MIN IN SUBR.((RCONS)) ***' WRITE(6,*) ' (R0,R0MIN) =',R0,R0MIN ! F20N = 4.97E-5 " CORIOLIS PARAMETER AT 20N ! IF(F0.GT.F20N) THEN ! (I) MODIFY DP R0 = R0MIN DP = GRF*(R0*R0+R15*R15)**1.5/R0/R15 ! ELSE ! (II) MODIFY R15 ! R0X = R0MIN ! GRF = (F0*VT+VT*VT/R15)*ROH ! BB = ROH/R0X/DP ! CC = (VT/F0)**2+4*DP/ROH*R0X/F0/VT ! XR0 = 0.5*(SQRT(CC)-VT/F0) ! YR0 = GRF*(R0X*R0X+XR0*XR0)**1.5/R0X/DP ! XR1 = YR0 ! YR1 = BB*(F0*VT+VT*VT/XR1)*(R0X*R0X+XR1*XR1)**1.5 ! ITE = 0 !25 CONTINUE ! ITE = ITE + 1 ! WEG = (YR1-YR0)/(XR1-XR0) ! XR = (YR1-XR1*WEG)/(1.-WEG) ! IF (ABS(XR-XR1) .LT. 1000.) GO TO 35 ! XR0 = XR1 ! XR1 = XR ! YR0 = YR1 ! YR1 = BB*(F0*VT+VT*VT/XR1)*(R0X*R0X+XR1*XR1)**1.5 ! IF (ITE .LT. ITEMAX) GO TO 25 !35 CONTINUE ! IF(XR.GT.R15) THEN ! R15 = XR ! R0 = R0MIN ! WRITE(6,605) R15,ITE !605 FORMAT(1H ,' MODIFIED 15M/S RADIUS (R15) =',F10.1 ! > ,' ITE =',I3 ) ! ENDIF ! ENDIF ENDIF PC = PE - DP ! CALCULATE MAXIMUM GRADIENT WIND SPEED FOR /* MONITOR */ ! WEGSTEIN'S ITERATION BB = 4.*DP*R0/(ROH*F0*F0) R02 = R0*R0 XR0 = 0. CC = R02+XR0 YR0 = 2./BB*(1.-SQRT(1.+BB/CC**1.5))*CC**2.5+2.*R02 XR1 = YR0 CC = R02+XR1 YR1 = 2./BB*(1.-SQRT(1.+BB/CC**1.5))*CC**2.5+2.*R02 ITE = 0 40 CONTINUE ITE = ITE + 1 WEG = (YR1-YR0)/(XR1-XR0) XR = (YR1-XR1*WEG)/(1.-WEG) IF (ABS(XR-XR1) .LT. 2.E6) GO TO 50 XR0 = XR1 XR1 = XR YR0 = YR1 CC = R02+XR1 YR1 = 2./BB*(1.-SQRT(1.+BB/CC**1.5))*CC**2.5+2.*R02 IF (ITE .LT. ITEMAX) GO TO 40 50 CONTINUE RMX = SQRT(XR) VMX = F0*RMX/2.*(-1.+SQRT(1.+BB/CC**1.5)) ! WRITE(6,610) RMX,VMX,ITE ! 610 FORMAT(1H ,' MAXIMUM GRADIENT WIND (RMX,VMX) =' ! > ,2F10.1,' ITE =',I3) RETURN END !---------------------------------------------------------------- !**************************************************************** !---------------------------------------------------------------- subroutine rdmean1d(rm1,rm2,zout,zin,im,jm,kz,rtg) dimension rtg(*),zin(im*jm,*),aa(50) real zout do K=1,KZ AA(K) = 0.0 enddo NN=0 do I=1,IM*JM if (RM2.GT.RTG(I) .AND. RTG(I).GT.RM1) THEN NN = NN + 1 do K=1,KZ AA(K) = AA(K) + ZIN(I,K) enddo endif enddo ANN = FLOAT(NN) do K=1,KZ ZOUT = AA(K) / ANN enddo return end !--------------------------------------------------------------- !*************************************************************** !--------------------------------------------------------------- SUBROUTINE RDMEAN ( RM1, RM2, ZOUT, ZIN, IM, JM, KZ, RTG ) DIMENSION RTG ( *), & ! IN " ZIN (IM*JM, *), & ! OUT " ZOUT( *), & ! WORK" AA (50) ! ------------------------------------------------ ! ------------------------------------------------ ! "ZOUT" = RADIAL MEAN OF "ZIN" OVER RM1 < R < RM2 ! ------------------------------------------------ ! ------------------------------------------------ DO 10 K=1,KZ AA(K) = 0.0 10 CONTINUE NN=0 DO 20 I=1,IM*JM IF (RM2.GT.RTG(I) .AND. RTG(I).GT.RM1) THEN NN = NN + 1 DO 100 K=1,KZ AA(K) = AA(K) + ZIN(I,K) 100 CONTINUE ENDIF 20 CONTINUE ANN = FLOAT(NN) DO 30 K=1,KZ ZOUT(K) = AA(K) / ANN 30 CONTINUE RETURN END !-------------------------------------------------------------- !************************************************************** !-------------------------------------------------------------- SUBROUTINE RH2MR ( WV, IM, JM, KM, RH, TP, PSL ) DIMENSION WV(IM,JM,KM),RH(IM,JM,KM),TP(IM,JM,KM),PSL(KM) ! ----------------------------------- ! REL. HUM. --> MIXING RATIO WV ! ----------------------------------- ! OUTPUT DATA WV MIXING RATIO (G/KG) ! IMPUT DATA RH RELATIVE HUMD.(%) ! TP TEMP (K) ! PSL PRESS (MB) TETEN(X) = 6.11 * EXP (17.2694 * X / (237.3 + X)) TZERO = 273.16 DO 20 K=1,KM DO 20 J=1,JM DO 20 I=1,IM WV(I,J,K)=RH(I,J,K) IF(WV(I,J,K) .GT. 95.)WV(I,J,K) = 95. IF(WV(I,J,K).LT. 0.) WV(I,J,K)= 0. 20 CONTINUE DO 30 K=1,KM DO 30 J=1,JM DO 30 I=1,IM TT=TP(I,J,K)-TZERO SVAP=TETEN(TT) VAP=SVAP * WV(I,J,K) * 0.01 WV(I,J,K)=(622.*VAP)/(PSL(K)-0.378*VAP) 30 CONTINUE RETURN END !--------------------------------------------------------------- !*************************************************************** !--------------------------------------------------------------- SUBROUTINE RH2TD ( TD, IM, JM, KM, RH, TT, PSL ) DIMENSION TT(IM,JM,KM),RH(IM,JM,KM),PSL(KM),TD(IM,JM,KM) ! ----------------------------------- ! REL. HUM. --> DEW POINT (TD :K) ! ----------------------------------- ! OUTPUT DATA TD DEW POINT (K) ! IMPUT DATA RH RELATIVE HUMD.(%) ! TT TEMP (K) ! PSL PRESS (MB) CALL RH2MR (TD, IM,JM,KM,RH, TT,PSL) DO 10 K=1,KM DO 10 J=1,JM DO 10 I=1,IM TDXX=273.16+237.3/(17.2694/ALOG(PSL(K)*TD(I,J,K)/(3800.42 & +2.31*TD(I,J,K)))-1.0) TD(I,J,K)=TDXX 10 CONTINUE RETURN END !--------------------------------------------------------------- !*************************************************************** !--------------------------------------------------------------- SUBROUTINE RLTLN ( GI, GJ, RLAT, RLON,SLAT, & NPROJC, DELS, SLON, XI, XJ, XLAT, XLON ) !********************************************************************** ! GRID POINT OF LATITUDE AND LONGITUDE * ! * !OUT-PUT (GI,GJ) GRID POINT * ! * !IN-PUT RLAT LATITUDE (DEGREE) -90.< < 90. * ! RLON LONGITUDE (DEGREE) 0.< <360. * ! NPROJC MAP PROJECTION * ! 'PSN ' POLAR STEREO SLAT= 60 N * ! 'PSS ' POLAR STEREO SLAT=-60 N * ! 'MER ' MERCATOR SLAT= 0. N * ! 'LMN ' LAMBERT SLAT= 30 N , 60 N * ! 'LMS ' LAMBERT SLAT=-30 N ,-60 N * ! DELS GRID PARAM ( M) * ! Y-COORDINATE 0.< <360.* * ! SLON STANDARD LONGITUDE* ! * ! (XI,XJ) <----------> (XLAT,XLON) * ! STANDARD POINT * !********************************************************************** CHARACTER*4 NPROJC ! RADIUS OF EARTH" RA=6371.E3 PAI=3.14159 RAD=PAI /180. !-----------------------------------POLAR STEREO----------------------- IF(NPROJC(1:3).EQ.'PSN') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=60. SLAT1=SLAT*RAD SLON1=SLON*RAD XLAT1=XLAT*RAD XLON1=XLON*RAD AL0=RA*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PAI-XLAT1)) X0=AL0*SIN(XLON1-SLON1) Y0=AL0*COS(XLON1-SLON1) RLAT1=RLAT*RAD RLON1=RLON*RAD ! call fij(GI,GJ,RLAT1,RLON1) ALI=RA*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PAI-RLAT1)) DI=ALI*SIN(RLON1-SLON1) DJ=ALI*COS(RLON1-SLON1) GI=(DI-X0)/DELS+XI GJ=(DJ-Y0)/DELS+XJ ENDIF IF(NPROJC(1:3).EQ.'PSS') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=-60. SLAT1=-SLAT*RAD SLON1= SLON*RAD XLAT1=-XLAT*RAD XLON1= XLON*RAD AL0=RA*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PAI-XLAT1)) X0=AL0*SIN(XLON1-SLON1) Y0=AL0*COS(XLON1-SLON1) RLAT1=-RLAT*RAD RLON1= RLON*RAD ALI=RA*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PAI-RLAT1)) DI=ALI*SIN(RLON1-SLON1) DJ=ALI*COS(RLON1-SLON1) GI=(DI-X0)/DELS+XI GJ=(Y0-DJ)/DELS+XJ ENDIF !-----------------------------------MERCATOR--------------------------- IF(NPROJC(1:3).EQ.'MER') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" ! SLAT=0. ! NOTE: AGAIN SET SLAT ! SLAT=42.3687 SLAT1=SLAT*RAD SLON1=SLON*RAD XLAT1=XLAT*RAD XLON1=XLON*RAD AC=RA*COS(SLAT1) X0=0. Y0=AC*LOG((1.+SIN(XLAT1))/COS(XLAT1)) RLAT1=RLAT*RAD RLON1=RLON*RAD ! ---------------------------- DI = RLON - XLON DI = MOD (DI+900., 360.) - 180. DI = AC*(DI*RAD) ! ---------------------------- ! ---- DI=AC*(RLON1-XLON1) DJ=AC*LOG((1.+SIN(RLAT1))/COS(RLAT1)) GI=(DI-X0)/DELS+XI GJ=(Y0-DJ)/DELS+XJ ENDIF !-----------------------------------LAMBERT---------------------------- IF(NPROJC(1:3).EQ.'LMN') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLATA=30. SLATB=60. PI4=PAI*0.25 SLATA1=SLATA*RAD SLATB1=SLATB*RAD SLAT1=PI4-SLATA1*0.5 SLAT2=PI4-SLATB1*0.5 SLON1=SLON*RAD XLAT1=XLAT*RAD XLON1=XLON*RAD CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) ACN=RA*COS(SLATA1)/CK R0=ACN/(TAN(SLAT1))**CK AL0=R0*(TAN(PI4-XLAT1*0.5))**CK ! ----------------------------------- XSLON = XLON-SLON XSLON = MOD (XSLON+900., 360.) - 180. XSLON = XSLON * RAD * CK X0=AL0*SIN(XSLON) Y0=AL0*COS(XSLON) ! ----------------------------------- ! ---- X0=AL0*SIN((XLON1-SLON1)*CK) ! ---- Y0=AL0*COS((XLON1-SLON1)*CK) RLAT1=RLAT*RAD RLON1=RLON*RAD ALI=R0*(TAN(PI4-RLAT1*0.5))**CK ! ----------------------------------- RXLON = RLON-XLON RXLON = MOD (RXLON+900., 360.) - 180. RXLON = RXLON * RAD * CK DI=ALI*SIN(RXLON + XSLON) DJ=ALI*COS(RXLON + XSLON) ! ----------------------------------- ! ---- DI=ALI*SIN((RLON1-SLON1)*CK) ! ---- DJ=ALI*COS((RLON1-SLON1)*CK) GI=(DI-X0)/DELS+XI GJ=(DJ-Y0)/DELS+XJ ENDIF IF(NPROJC(1:3).EQ.'LMS') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLATA=-30. SLATB=-60. PI4=PAI*0.25 SLATA1=-SLATA*RAD SLATB1=-SLATB*RAD SLAT1=PI4-SLATA1*0.5 SLAT2=PI4-SLATB1*0.5 SLON1= SLON*RAD XLAT1=-XLAT*RAD XLON1= XLON*RAD CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) ACN=RA*COS(SLATA1)/CK R0=ACN/(TAN(SLAT1))**CK AL0=R0*(TAN(PI4-XLAT1*0.5))**CK ! ----------------------------------- XSLON = XLON-SLON XSLON = MOD (XSLON+900., 360.) - 180. XSLON = XSLON * RAD * CK X0 = AL0*SIN(XSLON) Y0 = AL0*COS(XSLON) ! ----------------------------------- !--- X0=AL0*SIN((XLON1-SLON1)*CK) !--- Y0=AL0*COS((XLON1-SLON1)*CK) RLAT1=-RLAT*RAD RLON1= RLON*RAD ALI=R0*(TAN(PI4-RLAT1*0.5))**CK ! --------------------------- RXLON = RLON-XLON RXLON = MOD (RXLON+900., 360.) - 180. RXLON = RXLON * RAD * CK DI=ALI*SIN(RXLON + XSLON) DJ=ALI*COS(RXLON + XSLON) ! --------------------------- !---- DI=ALI*SIN((RLON1-SLON1)*CK) !---- DJ=ALI*COS((RLON1-SLON1)*CK) GI=(DI-X0)/DELS+XI GJ=(Y0-DJ)/DELS+XJ ENDIF !----------------------------LAT-LONG 'LL')--------------------------- IF(NPROJC(1:3).EQ.'LL') THEN GI = XI+INT((RLON-XLON)/(DELS/111000.0)) GJ = XJ-INT((RLAT-XLAT)/(DELS/111000.)+1) ! write(*,*)'RLON and RLAT=',RLON,RLAT ! write(*,*)'GI and GJ=',GI,GJ ENDIF !--------------------------------------------------------------------- RETURN END !--------------------------------------------------------------------- !********************************************************************* !--------------------------------------------------------------------- SUBROUTINE ROTANG ( ACCT, IM, JM, RTG, ANG, CLAT, CLON, & NPRO, DELS, SLON, XI, XJ, XLAT, XLON,SLAT ) ! CALL ROTANG (W2C, IM,JM, RTG,ANG, CLAT,CLON ! > ,NPRO, DELS, SLON, XI,XJ,XLAT,XLON ) ! ----------------------------------------------------------- ! CALCULATE IANGLE(DEG) FOR CYLINDRICAL-TO-CARTESIAN ! WIND COMPONENT TRANSFORM ! ACCT=0 CORRESPONDS TO (-X) CARTESIAN AXIS ! ----------------------------------------------------------- ! OUT DIMENSION ACCT(IM, *) & ! IN ,RTG (IM, *) ,ANG(IM, *) CHARACTER * 4 NPRO RAD = ASIN(1.0)/90.0 DEG = 1.0/RAD DEG = 180/3.1415926 ! ARBITRARY VIRTUAL DISPLACEMENT RINCRE = 50000. DO 100 J=1,JM DO 100 I=1,IM R = RTG(I,J) + RINCRE CALL CLY2LL (ALAT,ALON,ANG(I,J),R,CLAT,CLON) CALL RLTLN (FI,FJ,ALAT,ALON,SLAT, & NPRO,DELS,SLON, XI,XJ,XLAT,XLON ) DJ = FJ - J DI = FI - I Z = ACOS ( DI/max(SQRT(DI**2+DJ**2),epsilonR()) ) * DEG ! IF (DJ .GE. 0.) Z = 180. + Z ! IF (DJ .LT. 0.) Z = 180. - Z !recommented(Ron) IF (DJ .LT. 0.) Z = 180. + Z !original IF (DJ .GE. 0.) Z = 180. - Z ACCT(I,J) = Z 100 CONTINUE RETURN END !--------------------------------------------------------------------- !********************************************************************* !--------------------------------------------------------------------- SUBROUTINE SETDPS ( DPSEA, T, PHIS, PAI, SPLA1, PTOP, IM, JM ) ! ---------------------------------------------------------------- ! CALCULATE THE DIFFERENCE BETWEEN SURFACE AND SEA-LEVEL PRESSURES ! BY INTEGRATING AIR DENSITY UNDER CONSTANT LAPSE RATE ! ---------------------------------------------------------------- DIMENSION DPSEA(*), T(IM*JM, *), PHIS(*), PAI(*) RGAS = 287.04 GRAV = 9.8 TGRAD = 0.0065 AA = GRAV/(RGAS*TGRAD) RAA = 1./AA DO 100 I=1,IM*JM PS = PAI(I) + PTOP ! P1 = SIG1*PAI(I) + PTOP TS = T(I,1) + RAA*T(I,1) * ALOG(PS/SPLA1) ZS = PHIS(I)/GRAV BB = TS + TGRAD*ZS CC =-AA * ALOG(TS/BB) PP = PS * EXP(CC) DPSEA(I) = PP-PS 100 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE SETRDA ( RTG, ANG, IM, JM, ALAT, ALON, CLAT, CLON ) ! ------------------------------------------------------------+ ! +< OUTPUT > + ! + RTG : DISTANCE ON GREAT CIRCLE + ! + BETWEEN EACH GRID AND TYPHOON CENTRE + ! + ANG : IANGLE SHOWN IN THE SCHEMATI! DIAGRAM + ! + (N) + ! + | + ! + 270 < ANG < 360 | 180 < ANG < 270 + ! + | + ! + --- XD<0 ---------- (C) -------- XD>0 --------- + ! + | + ! + 0 < ANG < 90 | 90 < ANG < 180 + ! + | + ! + (S) + ! +< INPUT > + ! + ALON : LONGITUDE OF GRID + ! + ALAT : LATITUDE OF GRID + ! + CLON : LONGITUDE OF TYPHOON CENTRE + ! + CLAT : LATITUDE OF TYPHOON CENTRE + ! +-----------------------------------------------------------+ DIMENSION RTG(IM*JM), ANG(IM*JM), ALON(IM*JM), ALAT(IM*JM) REAL X1, Y1, Y1C, Y1S, X2, Y2, Y2C, Y2S, XD, ALPHA, RAD ,DEG, & Z, ZC, ZS R0 = 6371.E+3 RAD = ASIN(1.0)/90.0 DEG = 1.0/RAD X1 = CLON * RAD Y1 = CLAT * RAD Y1C = COS(Y1) Y1S = SIN(Y1) DO 100 I=1,IM*JM X2 = ALON(I) * RAD Y2 = ALAT(I) * RAD Y2C = COS(Y2) Y2S = SIN(Y2) XD = X2-X1 ZC = Y2C*Y1C* COS(XD) + Y2S*Y1S Z = ACOS(ZC) ZS = SIN(Z) RTG(I) = Z * R0 IF (XD.NE.0.) THEN ALPHA = min((Y2S - Y1S*ZC)/(Y1C*ZS),1.-epsilonR()) ALPHA = max(ALPHA,-1.+epsilonR()) ALPHA = ACOS(ALPHA) * DEG IF (XD.LT.0) THEN IF (ALPHA.LT.90.) ANG(I) = ALPHA + 270. IF (ALPHA.GE.90.) ANG(I) = ALPHA - 90. ELSE ANG(I) = 270. - ALPHA ENDIF ELSE ANG(I) = 90. IF (Y2.GT.Y1) ANG(I) = 270. ENDIF 100 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE SPLINE ( ZI, PILN, LMAX, Z, PLN, IMAX, IDX ) DIMENSION ZI(*), PILN(*), Z(*), PLN(*), & SM(90),H(90),AL(90),AM(90),AP(90),C(90) ! INTERPOLATION USING A CUBIC NONPERIODIC SPLINE ! ZI(I)........INTERPOLATED VALUE AT PILN(I) ! PILN(I)......LN(P) COORDINATE IN DECREASING ORDER ! LMAX.........NUMBER OF INTERPOLATION POINTS ! Z(I).........DATA AT PLN(I) ! PLN(I).......LN(P) COORDINATE IN DECREASING ORDER ! IMAX.........NUMBER OF DATA POINTS ! SM(I)........SECOND DERIVATIVES AT DATA POINTS ! ******* 1 INTERPOLATION FOR HEIGHT ! * IDX * 2 INTERPOLATION FOR WIND ! ******* 3 INTERPOLATION FOR TEMPERATURE ( CUBIC SPLINE FOR ! HEIGHT IS DIFFERENCIATED WITH RESPECT TO LOG P ! --- HYDROSTATIC RELATION ) IM1=IMAX-1 R=287.04 G=9.8 GR=-G/R DO 110 I=2,IMAX !!$ print*, 'pln: ',i,pln(i),pln(i-1),PLN(I)-PLN(I-1) 110 H(I)=PLN(I)-PLN(I-1) DO 120 I=2,IM1 AL(I)=0.5*H(I+1)/(H(I)+H(I+1)) 120 AM(I)=0.5-AL(I) IF( IDX.EQ.2 ) GO TO 80 ! END CONDITION FOR HEIGHT AND TEMPERATURE (LAPSE RATE IS CONSTANT) ! SM(1)=SM(2) ; SM(IMAX-1)=SM(IMAX) AL(1)=-1. AM(IMAX)=-1. GO TO 90 80 CONTINUE ! END CONDITION FOR WIND ( WIND SHEAR IS CONSTANT ) ! SM(1)=0. ; SM(IMAX)=0. AL(1)=0. AM(IMAX)=0. 90 CONTINUE AL(IMAX)=0.0 DO 130 I=2,IMAX AP(I)=1.0/(1.0-AL(I-1)*AM(I)) 130 AL(I)=AL(I)*AP(I) C(1)=0. C(IMAX)=0. DO 160 I=2,IM1 160 C(I)=3.0*((Z(I+1)-Z(I))/H(I+1)-(Z(I)-Z(I-1))/H(I)) & / (H(I)+H(I+1)) !!$160 print*, i,z(i+1),z(i),z(i-1),h(i+1),h(i),c(i) ! FORWARD SUBSTITUTION DO 200 I=2,IMAX !!$ print*, i,c(i),c(i-1),am(i),ap(i) 200 C(I)=(C(I)-C(I-1)*AM(I))*AP(I) !!$ print*, 'imax: ',imax,c(imax) !!$ print*, 'c: ',c !!$ print*, 'am: ',am(imax),am !!$ print*, 'ap: ',ap(imax),ap SM(IMAX)=C(IMAX) ! BACKWARD SUBSTUTUTION DO 220 K=1,IM1 I=IMAX-K 220 SM(I)=C(I)-AL(I)*SM(I+1) ! INTERPOLATION IB=2 DO 500 L=1,LMAX X=PILN(L) DO 300 I=IB,IMAX IF(X.GE.PLN(I)) GOTO 310 300 CONTINUE I=IMAX 310 IB=I IF(IDX.EQ.3) GO TO 400 ! print*, 'gate .5: ',i,pln(i),pln(i-1),x,h(i),z(i),sm(i),sm(i-1) ZI(L)=(PLN(I)-X)/H(I)* & (Z(I-1)-SM(I-1)/6.*(X-PLN(I-1))*(H(I)+PLN(I)-X)) & +(X-PLN(I-1))/H(I)* & (Z(I)-SM(I)/6.*(PLN(I)-X)*(H(I)+X-PLN(I-1))) ! print*, 'gate 1: ',l,zi(l) GO TO 500 400 CONTINUE ! DIFFERENTIAL CALCULUS OF CUBIC SPLINE FOR HEIGHT ! GR=-G/R ... COEFFICIENT IN HYDROSTATIC EQUATION ZI(L)=SM(I-1)*(-(PLN(I)-X)**2/(2.*H(I))+H(I)/6.) & +SM(I)*((X-PLN(I-1))**2/(2.*H(I))-H(I)/6.) & +(Z(I)-Z(I-1))/H(I) ZI(L)=ZI(L)*GR 500 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TCLOUD ( TCL, QCL, POUT, KM, PSTART, TSTART, KMIN ) ! TCLB = TSEA (sst) ! MOIST ADIABATIC LAPSE RATE FROM TSTART AND PSTART DIMENSION TCL(KM), QCL(KM), POUT(KM) RGAS=287.04 CL=2.5E6 CP=1004.6 CLBYCP=CL/CP AKAPA=RGAS/CP ! --------------------------------------------- ! BELOW THE PSTART LEVEL ! TCL AND QCL ARE NOT GIVEN IN THIS ROUTINE ! --------------------------------------------- KMIN = 1 DO 100 K=1,KM KMIN=K+1 IF(POUT(K).GE.PSTART) THEN ! write(*,*)'checking' ! write(*,*)'POUT',POUT(K),' PSTART',PSTART GO TO 101 ELSE ENDIF 100 CONTINUE 101 CONTINUE PINV=1./PSTART CALL TETENS(TSTART,PINV,QSAT,DQSAT) IF (KMIN .GT. 1) THEN DO 110 K=1,KMIN-1 TCL(K)=TSTART QCL(K)=QSAT 110 CONTINUE ENDIF K = KMIN TCL(K)=TSTART*(POUT(K)/PSTART)**AKAPA QCL(K)=QSAT PINV=1./POUT(K) DO 120 ITR=1,2 CALL TETENS(TCL(K),PINV,QSAT,DQSAT) COND=(QCL(K)-QSAT)/(1.+DQSAT) TCL(K)=TCL(K)+CLBYCP*COND ! write(*,*)'TCL in TCLOUD',TCL(K) QCL(K)=QCL(K)-COND 120 CONTINUE DO 130 K=KMIN+1,KM TCL(K)=TCL(K-1)*(POUT(K)/POUT(K-1))**AKAPA QCL(K)=QCL(K-1) PINV=1./POUT(K) DO 131 ITR=1,2 CALL TETENS(TCL(K),PINV,QSAT,DQSAT) COND=(QCL(K)-QSAT)/(1.+DQSAT) TCL(K)=TCL(K)+CLBYCP*COND QCL(K)=QCL(K)-COND 131 CONTINUE 130 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TDTORH ( RH, IM, KM, TD, TT, W1, W2 ) ! CONVERSION FROM DEWPOINT TO REL.HUM(%) / TETEN'S EQUATION DIMENSION RH(IM,KM),TD(IM,KM),TT(IM,KM),W1(IM),W2(IM) ! OUT----- IN---------------- WORK-------- T0=273.16 ! PRINT *,'T DEW T ',TD(1,1),TT(1,1) ! PRINT *,'T DEW T ',TD(10,2 ),TT(10,2 ) ! PRINT *,'T DEW T ',TD(20,3 ),TT(20,3 ) ! PRINT *,'T DEW T ',TD(30,4 ),TT(30,4 ) DO 50 K=1,KM DO 30 I=1,IM W1(I) = TD(I,K) -T0 30 CONTINUE CALL TETENZ (W2, W1, IM) DO 35 I=1,IM W1(I) = TT (I,K) - T0 35 CONTINUE CALL TETENZ (RH(1,K), W1, IM) DO 40 I=1,IM RH(I,K) = 100.0 * W2(I) / RH(I,K) 40 CONTINUE 50 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TECORR ( T00, SPLA, KM, PSEA0, TSEA0 ) ! ///////////////////////////////////// ! ADJUST SURF. TEMP. TO SEA SURF. TEMP. ! ///////////////////////////////////// PARAMETER ( LCOR = 2 ) DIMENSION T00(KM), SPLA(KM) & ,PLN0X(20),TM0X(20),PLNXX(LCOR),TMXX(LCOR) TM0X(1) = TSEA0 PLN0X(1) = LOG(PSEA0) DO 305 K=2,4 PLN0X(K) = LOG(SPLA(K+LCOR-1)) TM0X (K) = T00(K+LCOR-1) 305 CONTINUE DO 310 L=1,LCOR PLNXX(L) = LOG(SPLA(L)) 310 CONTINUE CALL SPLINE(TMXX,PLNXX,LCOR,TM0X,PLN0X,4,2) DO 320 L=1,LCOR T00(L) = MAX(T00(L),TMXX(L)) 320 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TEMP3D ( DATA3S, IM, JM, KM, RIJ, SPLA, & DATA2P, IRM, KPM, PLG, DR, W1A, W1B, W1C ) ! ----------------------------------------------------------------- ! PROJECTION OF SCALAR DATA FROM AXISYMMETRIC 2-D FIELD (P-LEVEL) ! INTO 3-D FIELD (SIGMA-LEVEL) ! ----------------------------------------------------------------- DIMENSION DATA3S(IM,JM,KM), W1A(KM), W1B(KM), & DATA2P(IRM,KPM), W1C(KPM), & SPLA(KM), PLG(KPM), RIJ(IM,JM) ICOPY=0 DRIN=1./DR DO 100 I=1,IM DO 100 J=1,JM IR = INT(RIJ(I,J)*DRIN) + 1 DDR = RIJ(I,J) - DR*(IR-1) IF(IR.GE.IRM) THEN IF(ICOPY.EQ.0)THEN ICOPY = 1 III = I JJJ = J IR = IRM-1 DDR = DR ELSE DO 101 K=1,KM DATA3S(I,J,K) = DATA3S(III,JJJ,K) 101 CONTINUE GO TO 100 ENDIF ENDIF DO 102 K=1,KM ! POUT = SPLA(K) W1A(K) = LOG(SPLA(K)) 102 CONTINUE DO 103 IP=1,KPM W1C(IP) = DATA2P(IR,IP) + & (DATA2P(IR+1,IP)-DATA2P(IR,IP))*DDR*DRIN 103 CONTINUE CALL SPLINE (W1B,W1A,KM,W1C,PLG,KPM,2) DO 110 K=1,KM DATA3S(I,J,K) = W1B(K) 110 CONTINUE 100 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TETENS ( TMP, PRSIV, SATQ, DSATQ ) ! COMPUTES SATURATION Q PARAMETER ( CA=6.11, CC=7.5, CE=0.622, CF=237.3, CL=2.5E6 ) PARAMETER ( CP=1004.6, TZERO=273.2 ) PARAMETER ( CSAT=CE*CA, CDSAT=CA*CC*CE*CF*CL/CP ) TEMP=TMP-TZERO X1=10.**(CC*TEMP/(CF+TEMP)) XX=PRSIV*X1 SATQ=CSAT*XX DSATQ = CDSAT * ALOG(10.0) * XX / ( (TEMP+CF)*(TEMP+CF) ) RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TETENZ ( VAP, TDC, IM ) ! --- OUTPUT VAP(IM) : VAPOUR PRESSURE (MB) ! --- INPUT TDC(IM) : DEWPOINT ( C) DIMENSION VAP(*), TDC(*) CLG = 7.5*LOG(10.0) DO 10 I=1,IM VAP(I) = CLG * TDC(I) / (237.3 + TDC(I)) 10 CONTINUE DO 20 I=1,IM VAP(I) = EXP (VAP(I)) 20 CONTINUE DO 30 I=1,IM VAP(I) = 6.11 * VAP(I) 30 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TYMOVE ( UCOMP, VCOMP, NPROJC, DELS, SLON, XI, XJ, & XLAT, XLON, TLAT0, TLON0, TLAT6, TLON6, & TLAT12, TLON12, UA, VA, UCMP12, VCMP12 ) PARAMETER ( INTOBS=06 ) DIMENSION SM(1), RLAT(1) CHARACTER*4 NPROJC ! CALCULATE TYPHOON MOVEMENT (M/S) FROM TWO-TIME OBSERVATIONS ! UA,VA ARE CALCULATED USING 12 HR POSITION AS WELL CALL RLTLN (FI0,FJ0,TLAT0,TLON0,SLAT, & NPROJC,DELS,SLON,XI,XJ,XLAT,XLON) CALL RLTLN (FI6,FJ6,TLAT6,TLON6,SLAT, & NPROJC,DELS,SLON,XI,XJ,XLAT,XLON) CALL RLTLN (FI12,FJ12,TLAT12,TLON12,SLAT, & NPROJC,DELS,SLON,XI,XJ,XLAT,XLON) RLAT(1) = 0.5*(TLAT0+TLAT6) CALL MAPFCT (SM,1,1,RLAT,NPROJC,SLAT) DELS0 = DELS / SM(1) XCOMP = (FI6-FI0)*DELS0 YCOMP = (FJ6-FJ0)*DELS0 U3HR = XCOMP/(60.*60.*INTOBS) V3HR = YCOMP/(60.*60.*INTOBS) RLAT(1) = 0.5*(TLAT6+TLAT12) CALL MAPFCT (SM,1,1,RLAT,NPROJC,SLAT) DELS0 = DELS / SM(1) XCOMP = (FI12-FI6)*DELS0 YCOMP = (FJ12-FJ6)*DELS0 U9HR = XCOMP/(60.*60.*INTOBS) V9HR = YCOMP/(60.*60.*INTOBS) RLAT(1) = 0.5*(TLAT0+TLAT12) CALL MAPFCT (SM,1,1,RLAT,NPROJC,SLAT) DELS0 = DELS / SM(1) XCOMP = (FI12-FI0)*DELS0 YCOMP = (FJ12-FJ0)*DELS0 UCMP12 = XCOMP/(60.*60.*INTOBS*2) VCMP12 = YCOMP/(60.*60.*INTOBS*2) ! UCOMP,VCOMP USE THE APPROX. U0HR=U3HR, V0HR=V3HR UCOMP=U3HR VCOMP=V3HR ! U0HR=U3HR+(ACC.U)3HR*(T=3HRS) ! APPROX. (ACC.U)3HR BY (ACC.U)6HR = (U3HR-U9HR)/(T=6HRS) UA=U3HR*1.5-U9HR*.5 VA=V3HR*1.5-V9HR*.5 ! PRINT *,'#### SUBR. TYMOVE:' ! PRINT *,'CURRENT MOTION BY VEL(0)=VEL(-3)=(X(0)-X(-6))/6HRS:' ! PRINT *,'UCOMP=',UCOMP,' VCOMP=',VCOMP ! PRINT *,'CURRENT MOTION BY VEL(0)=VEL(-3)*1.5-VEL(-9)*.5:' ! PRINT *,'UCOMP=',UA,' VCOMP=',VA UA = 0.5*(U3HR + U9HR) VA = 0.5*(V3HR + V9HR) ! UA = 0.0 ! VA = 0.0 ! PRINT *,' AV. MOTION OVER 12HRS = 0.5*(U3HR + U9HR) ' ! PRINT *,' UCOMP = ',UA,' VCOMP = ',VA RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE TYPLNT ( PSEA, PAI, U, V, T, WV, TSEA, & ALAT, ALON, PTOP, IM, JM, KM, & PHIS, ZP, DELSPL, SPLA, DELS, & XI, XJ, XLAT, XLON, SLON, NPRO, & CPDIF, RMW, LOCBY, & SLAT, BASICP, CRMTNP, PRFRAC, W3C, IPP, & myNewPosition, newLatLong, dim) ! ************************************ ! * TYPHOON VORTEX TRANSPLANTATION * ! * INTO * ! * GLOBAL ANALYSIS FIELDS * ! ************************************ ! CREATED JAN 12 1987 BY T.IWASAKI ! REVISED NOV 16 1987 BY M.UENO ! OCT 14 1988 BY M.UENO ! include 'dim.h' PARAMETER ( KPM=41, IRM=051 ) PARAMETER ( MXTY=10 ) ! ======== ANALYSIS ============== DIMENSION & PSEA(IM,JM), PAI(IM,JM), U(IM,JM,KM), V(IM,JM,KM) & ,T(IM,JM,KM), WV(IM,JM,KM), TSEA(IM,JM), PHIS(IM,JM) & ,ZP(IM,JM,KM), SPLA(KM) & ,ALAT(IM,JM), ALON(IM,JM) & ,DELSPL(KM) ! ======== TYPHOON =============== DIMENSION & PSEAM(IM,JM), PAIM(IM,JM), UM(IM,JM,KM) & ,VM(IM,JM,KM), TM(IM,JM,KM), WVM(IM,JM,KM) & ,T00(KM), RH00(KM), DPSEA(IM,JM),ZPM(IM,JM,KM) & ,RTG(IM,JM), ANG(IM,JM), PIN(KM) DIMENSION & PSEAR(IRM), TR(IRM,KPM), RHR(IRM,KPM) & ,GRADR(IRM,KPM), POUTLG(KPM) ! ======== WORK AREA ============= DIMENSION W1A(KM), W1B(KM), W1C(KM), W1D(KM) & ,W1E(KPM), W1F(KPM), W1G(KPM), W1H(KPM) & ,W1I(KPM), W1J(KPM), W1K(KPM), W1L(KPM) & ,W1N(IRM), W1O(IRM) & ,W2A(IRM,KM), W2B(IRM,KM), W2H(IRM,KPM) & ,W2C(IM,JM), W2D(IM,JM), W2E(IM,JM) & ,W2F(IM,JM), W2G(IM,JM) & ,W3A(IM,JM,KM), W3B(IM,JM,KM) & ,W3C(IM,JM,KM) ! EQUIVALENCE ! > (TM(1,1,1),W2D), (TM(1,1,2),W2E) ! > ,(TM(1,1,3),W2F), (TM(1,1,4),W2G) ! ======== LAT/LONG INFO ========== integer, dimension(2) :: myNewPosition,centrePosition real, dimension(2) :: newLatLong COMMON /TYDATA/ NTY, NUM(MXTY), TPC(MXTY), T15(MXTY), TLO(MXTY) & ,TLA(MXTY), UTY(MXTY), VTY(MXTY), XTY(MXTY) & ,YTY(MXTY), TCPC(MXTY) CHARACTER*4 NPRO, LOCBY !!$ interface !!$ SUBROUTINE INVBAL ( U, V, PZ, SF, VP, FF, F, WORK, IL, JL, KL, & !!$ NPRO, DELS, SLON, XI, XJ, XLAT, XLON, SLAT ) !!$ real, dimension(:,:,:) :: u,v,pz !!$ real, dimension(:,:) :: sf,ff,work,vp,f !!$ integer :: il,jl,kl !!$ real :: dels,slon,xi,xj,xlat,xlon,SLAT !!$ character(len=*) :: npro !!$ end subroutine invbal !!$ end interface ! set kw1 value kw1=km ! initialize arrays do j=1,jm do i=1,im pseam(i,j)=0. paim(i,j)=0. dpsea(i,j)=0. rtg(i,j)=0. ang(i,j)=0. do k=1,km um(i,j,k)=0. vm(i,j,k)=0. tm(i,j,k)=0. wvm(i,j,k)=0. zpm(i,j,k)=0. enddo enddo enddo ! ================================================ ! >>>> PREPROCESS <<<< ! ================================================ ! 0: WITHOUT OR 1: WITH SURFACE MODIFICATION ISURF = 1 ! ---- SAVE SEA-LEVEL CORRECTION OF SURFACE PRESSURE ------ DO 120 I=1,IM DO 120 J=1,JM DPSEA(I,J) = PSEA(I,J) - PAI(I,J) 120 CONTINUE ! -------- ( TV ---> T , WV ---> RH ) ------------------------ ! CALL VRTEMP (T, T,WV, IM,JM,KM, 'G/KG', 'VTOR') CALL QQRELH ( WV, WV, T, SPLA, IM, JM, KM, 'G/KG', 'QTOH' ) ! ============================================== ! >>>> <<<< ! >>>> TYPHOON LOOP START <<<< ! >>>> <<<< ! ============================================== DO 1999 NT=1,NTY ! WRITE(6,6000) NUM(NT),TPC(NT),T15(NT),UTY(NT),VTY(NT) ! > ,TLA(NT),TLO(NT) 6000 FORMAT(1H0,'T',I4,2X,'PCNTR =',F7.1,2X,'R15 =',F7.1,2X & ,'INITIAL MOVEMENT(U,V) =',2F6.1,2X & ,'POSITION (LAT,LON) =',2F7.1 ) PCNT = TPC(NT) IF (LOCBY.EQ.'WIND') THEN ! >>> CALCULATE ANAL. CENTRE USING LOCCEN: ! - USE LEVEL 2 ( 850 mb ) WINDS TO FIND CENTRE ! - SEARCH BOX IS 14*14 DEGREES, CENTRE MUST BE WITHIN ! 4 DEGREES LAT/LON OF TLA/TLO ( MAXDIS=4. ) CALL LOCCEN ( U, V, 4, & TLA(NT)-7, TLO(NT)-7, TLA(NT)+7, TLO(NT)+7, & 4.0, BLAT, BLON, IM, JM, KM, NPRO, DELS, SLON, & XI, XJ, XLAT, XLON, SLAT ) IF (BLON.LT.-99) THEN PRINT *,'##LOCCEN FAILED... USING SLP MINIMUM LOCATION' LOCBY='PRES' ENDIF ENDIF IF (LOCBY.NE.'WIND') THEN ! >>> CALCULATE ANAL. CENTRE USING MSLP: YCNT = YTY(NT) XCNT = XTY(NT) ! CALL CENTR2 (PSEA, IM,JM, XCNT,YCNT,PCNT,7) centrePosition = minloc(psea) xcnt = centrePosition(1) ycnt = centrePosition(2) pcnt = minval(psea) WRITE(6,*)'CENTR2### USING MSLP GIVES:' WRITE(6,*)'(X GRIDPT ,Y GRIDPT ,CENTRAL VALUE)=',XCNT,YCNT,PCNT CALL XYTOLL (BLAT,BLON, XCNT,YCNT,SLAT & ,NPRO, DELS, SLON, XI,XJ,XLAT,XLON ) WRITE(6,*)'ANALYSIS HURRICANE LAT,LON =',BLAT,BLON ENDIF CALL RLTLN (CI,CJ,BLAT,BLON,SLAT, & NPRO,DELS,SLON, XI,XJ,XLAT,XLON ) CALL INTPLZ (W1A, CI,CJ, PSEA,IM,JM,1) ANLCP = W1A(1) PRINT *,'## PRESSURE AT ANALYSIS CENTER: ',ANLCP PRINT *,'## BOGUS CENTRAL PRESSURE: ',TPC(NT) PC = ANLCP*(PRFRAC)+TPC(NT)*(1.0-PRFRAC) IF(PRFRAC .LT. 0.) PC = ANLCP + PRFRAC ! ONLY ALLOW PRESSURES CPDIF hPa LOWER THEN ANALYSED IF( ANLCP-PC .GT. CPDIF ) PC = ANLCP - CPDIF TPC(NT) = PC ! PRINT *,'##COMBINED (FINAL) C.P. ',TPC(NT) ! ----------- ! RADIUS REDUCTION PARAMETER ! T15 = RADIUS OF OUTER CLOSED ISOBAR ! R15 = ROCI APPROX. RPARM = 1.0 ! ----------- R15 = T15 (NT) * 1000. * RPARM CLON = TLO (NT) CLAT = TLA (NT) ! NOTE: ! To set a new location for the bogused vortex, set CLAT and CLONG ! to the necessary latitude and longitude location ! MUST PROVIDE A BACKGROUND FIELD TO IMPLANT NEW VORTEX INTO IF MOVING VORTEX ONLY if (myNewPosition(1) > 0) then CLAT = newLatLong(1) CLON = newLatLong(2) PRINT*,'WARNING: DISPLACEMENT OF VORTEX FROM ORIGINAL LOCATION ', & 'SHOULD BE MERGED INTO SMOOTHED ANALYSIS HURRICANE FIELDS ', & 'PROVIDED BY USER' endif ! ------------------------------------------------ ! ==== PREPARATION FOR TRANSPLANTATION ==== ! ------------------------------------------------ ! ----- CORIOLIS PARAMETER AT TYPHOON CENTER ----- ! ANGULAR VEL. OF THE EARTH 00011500 OMEGA = 2.*(2.*ASIN(1.))/24./60./60. RAD = ASIN(1.)/90. F0 = 2. * OMEGA * SIN(RAD*CLAT) ! WRITE(6,620) F0 620 FORMAT(1H ,7X,'CORIOLIS PARAM. AT TYPHOON CENTER =',E12.5) ! ----- DISTANCE BETWEEN EACH GRID AND TYPHOON CENTER ----- CALL SETRDA (RTG,ANG,IM,JM,ALAT,ALON,CLAT,CLON) ! ----- NEAREST GRID POINT BY CENTER ----- ! USED IN SUBR.((OUTY)) AND ((CNCALL)) ! PRINT *, 'GNI,GNJ,CLAT,CLON,NPRO,DELS,SLON,XI,XJ,XLAT,XLON', ! > GNI,GNJ,CLAT,CLON,NPRO, DELS, SLON, XI,XJ,XLAT,XLON CALL RLTLN (GNI,GNJ, CLAT,CLON,SLAT, & NPRO, DELS, SLON, XI,XJ,XLAT,XLON ) ! PRINT *, 'GNI,GNJ,CLAT,CLON,NPRO,DELS,SLON,XI,XJ,XLAT,XLON', ! > GNI,GNJ,CLAT,CLON,NPRO, DELS, SLON, XI,XJ,XLAT,XLON ICNT = GNI JCNT = GNJ XTY(NT) = GNI YTY(NT) = GNJ INCR = 1 CALL OUTY (PSEA,PAI,T,WV,U,V,IM,JM,KM,ICNT,JCNT,INCR, & 'BEFORE TRANSPLANTATION ') ! ------------------------------------------------------------ ! | IN CASE OF TYPHOON LOCATED CLOSE TO THE LATERAL BOUNDARY | ! | AND LOCATED OUT OF THE LATERAL BOUNDARY | ! | TRANSPLANTATION IS NOT PERFORMED | ! ------------------------------------------------------------ ! =============== RBMIN = 400.E+3 ! =============== IF(ICNT.LE.1.OR.ICNT.GE.IM.OR.JCNT.LE.1.OR.JCNT.GE.JM) THEN PRINT *,'##OFF THE MAP - NO TRANSPLANT ',ICNT,JCNT TPC(NT) = -1. ! NOT TRANSPLANT AND TRACK IN THE FORECAST GO TO 1000 ENDIF DO 260 I=1,IM IF(RTG(I,1).LT.RBMIN.OR.RTG(I,JM).LT.RBMIN) THEN PRINT *,'##NEAR N OR S EDGE - NO TRANSPLANT ',ICNT,JCNT TPC(NT) = -2. ! NOT TRANSPLANT AND TRACK IN THE FORECAST GO TO 1000 ENDIF 260 CONTINUE DO 261 J=1,JM IF(RTG(1,J).LT.RBMIN.OR.RTG(IM,J).LT.RBMIN) THEN PRINT *,'##NEAR E OR W EDGE - NO TRANSPLANT ',ICNT,JCNT TPC(NT) = -3. ! NOT TRANSPLANT AND TRACK IN THE FORECAST GO TO 1000 ENDIF 261 CONTINUE ! -------------------------------------- ! ==== SET VARIOUS PARAMETERS ==== ! -------------------------------------- ! --- SET BOUND RADII OF AVERAGING ZONE FOR ENVIRONMENTAL VALUE --- ! RM1 < R < RM2 ! RM1 AND RM2 DEFINE THE ANNULUS OVER WHICH THE MEAN LARGE SCALE ! ENVIRONMENT IS OBTAINED RM2 = 2.*R15 RM2 = MAX(RM2,500.E3) RM1 = 0.6*RM2 ! 86.08.06 RM1 = MIN(RM1,500.E3) ! --- ENVIRONMENTAL SURF. PRES. & SEA SURF. TEMP. --- ! CALL RDMEAN1d (RM1,RM2,PSEA0,PSEA,IM,JM,1,RTG) !PSEA0 error CALL RDMEAN1d (RM1,RM2,TSEA0,TSEA,IM,JM,1,RTG) !TSEA0 error ! --- SET PC,R0 (PARAMETERS IN FUJITA'S FORMULA) --- ! WRITE(6, *) ! WRITE(6, *) '************************************' ! WRITE(6, *) '*** TRANSPLANTATION-PARAMETERS ***' ! WRITE(6, *) '************************************' ! WRITE(6, *) ' BEFORE SUBR.((RCONS))' ! PRINT *,'PC,R15,F0,DELS',PC,R15,F0,DELS ! WRITE(6,625) R15,PC,PSEA0,TSEA0 CALL RCONS (R0,DP,PC,R15, F0,PSEA0,TSEA0,DELS) ! TPC(NT) = PC ! T15(NT) = R15*0.001 WRITE(6, *) ' AFTER SUBR.((RCONS))' WRITE(6,625) R15,PC,PSEA0,TSEA0 WRITE(6, *) 'DP =',DP WRITE(6, *) 'R0 = ',R0 WRITE(6, *) 'RMW = ',RMW IF ( IPP .EQ. 2 ) THEN ! RESET CENTRAL PRESSURE BACK TO ORIGINAL VALUE FOR HOLLAND ! SCHEME ! PC = TCPC(NT) ENDIF ! --- SET R1,R2 (BOUND RADII OF TRANSITION ZONE IN SUPERPOSITION) R2 = MIN(2.*R15,R15+300.E+3) R1 = 0.5 * R2 ! --- SET R3 ( DP/DR=0 AT R3 ) ------------------------------- R3 = 2.0 * R15 ! WRITE(6,626) RM1,RM2,R0,R1,R2,R3 625 FORMAT(1H ,6X,'R15 =',F12.1,4X,'PC =',F12.1,4x,'PSEA0 =',F12.1, & 4x,'TSEA0 =',F12.1) 626 FORMAT(1H ,6X,'RM1 =',F12.1,4X,'RM2 =',F12.1,4X,'R0 =',F12.1,4X, & 'R1 =',F12.1,4X,'R2 =',F12.1,4X,'R3 =',F12.1) CALL RDMEAN (RM1,RM2,T00 ,T ,IM,JM,KM,RTG) CALL RDMEAN (RM1,RM2,RH00,WV,IM,JM,KM,RTG) ! WRITE (6, *) ' BEFORE SUBR.((TECORR))' ! WRITE (6,641) T00,RH00,PSEA0,TSEA0 ! --- HUMIDITY MODIFICATION --- RH00(1) = 90. RH00(KM) = 5. DO 300 K=2,KM-1 IF (RH00(K).GT.90.) RH00(K) = 90. IF (RH00(K).LT.5. ) RH00(K) = 5. 300 CONTINUE ! --- LOW-LEVEL TEMPERATURE CORRECTION --- CALL TECORR ( T00, SPLA, KM, PSEA0, TSEA0 ) ! WRITE (6, *) ' AFTER SUBR.((TECORR))' ! WRITE(6,641) T00,RH00,PSEA0,TSEA0 641 FORMAT(1H , 6X,'T00 =',15F10.1 & / 7X,'RH00 =',15F10.1 & / 7X,'PSEA0,TSEA0 =',2F10.1) ! WRITE(6, *) '************************************' ! --------------------------------------------------------------- ! ==== CREATION AND TLANSPLANTATION OF 2-D IDEAL TYPHOON ==== ! --------------------------------------------------------------- DO 310 K=1,KM PIN(K) = SPLA(K) 310 CONTINUE CALL AXISYM ( PSEAR, TR, GRADR, RHR, POUTLG, IRM, KPM, DR, & R0, R3, R15, PSEA0, PC, T00, TSEA0, RH00, PIN, KM, & W1C, W1D, KW1, W2H, & W1E, W1F, W1G, W1H, W1I, W1J, W1K, W1L, & IPP, TCPC(NT), RMW ) ! ----------------------------------------------------- ! ==== PROJECT INTO EXTENDED DIMENSIONAL SPACE ==== ! ----------------------------------------------------- ! <<< PAI >>> CALL PSEA2D (PSEAM,IM,JM,RTG,PSEAR,IRM,DR) DO 320 I=1,IM DO 320 J=1,JM PAIM(I,J) = PSEAM(I,J) - DPSEA(I,J) - PTOP 320 CONTINUE ! <<< WIND >>> !!!!! U AND V ARE FINE HERE ! ==== EXTRACT BASIC FLOW FROM INITIAL WIND FIELD ==== VCOMP = VTY(NT) UCOMP = UTY(NT) ! WRITE(6,*)' (UCOMP,VCOMP) =',UCOMP,VCOMP IF (PCNT .GT. 800.) THEN CALL SETRDA (W2F,W2G,IM,JM,ALAT,ALON,BLAT,BLON) ! >>>> BLAT,BLON MARK THE CENTRE OF THE ANALYSIS TYPHOON ... CLAT and CLON ! are the centre of the bogus typhoon - not the same for vortex replacement. ! U AND V ARE SKREWED COMING OUT OF HERE... BUT WHY? THEY'RE NOT SUPPOSED TO BE REDEFINED IN BFLOW ! CALL BFLOW (W3A,W3B, U, V, IM,JM,KM, BLAT,BLON, W2F,W2G & ! ,R2, DELS, XI,XJ,XLAT,XLON, SLON, NPRO & ! ,W2C,W2D,W2E,SLAT) !original call to bflow CALL BFLOW (W3A,W3B, U, V, IM,JM,KM, CLAT,CLON, RTG,ANG & ,R2, DELS, XI,XJ,XLAT,XLON, SLON, NPRO & ,W2C,W2D,W2E,SLAT) ! McGill ELSE PRINT *,'##**## BFLOW NOT CALLED! ##**##' pause DO 20 K = 1, KM DO 20 J = 1, JM DO 20 I = 1, IM W3A(I,J,K) = U(I,J,K) W3B(I,J,K) = V(I,J,K) 20 CONTINUE ENDIF IF (BASICP.LT.1.0) THEN ! PRINT *,'##MULTIPLYING BASIC FLOW BY FACTOR:',BASICP DO 7865 I=1,IM DO 7865 J=1,JM DO 7865 K=1,KM W3A(I,J,K)=W3A(I,J,K)*BASICP W3B(I,J,K)=W3B(I,J,K)*BASICP 7865 CONTINUE ENDIF IF (CRMTNP.GT.0.0) THEN ! PRINT *,'##DOING UV CORRECTION USING CURRENT MOTION' ! PRINT *,' FACTOR=',CRMTNP,' UEFF,VEFF ',UCOMP*CRMTNP, ! > VCOMP*CRMTNP DO 40 K = 1, KM DO 40 J = 1, JM DO 40 I = 1, IM W3A(I,J,K) = W3A(I,J,K) + UCOMP * CRMTNP W3B(I,J,K) = W3B(I,J,K) + VCOMP * CRMTNP 40 CONTINUE ENDIF CALL ROTANG (W2C, IM,JM, RTG,ANG, CLAT,CLON & ,NPRO, DELS, SLON, XI,XJ,XLAT,XLON,SLAT ) CALL WIND3D ( UM, VM, IM, JM, KM, DELSPL, SPLA, RTG, ANG, & GRADR, TR, POUTLG, KPM, IRM, DR, F0, & W2A, W2B, W1A, W1B, W1E, W1N, W1O, ISURF ) ! PRINT *,'##DOING FULL IMPLANTATION' DO 10 K = 1, KM DO 10 J = 1, JM DO 10 I = 1, IM UM(I,J,K) = UM(I,J,K) + W3A(I,J,K) !add background flow VM(I,J,K) = VM(I,J,K) + W3B(I,J,K) !add background flow 10 CONTINUE ! <<< TEMP. >>> CALL TEMP3D ( TM, IM, JM, KM, RTG, SPLA, TR, IRM, KPM, POUTLG, & DR, W1A, W1B, W1E ) ! <<< HUMD. >>> CALL TEMP3D ( WVM, IM, JM, KM, RTG, SPLA, RHR, IRM, KPM, POUTLG, & DR, W1A, W1B, W1E ) ! -------------------------------------------------- ! ==== TRANSPLANT INTO ENVIRONMENTAL FIELDS ==== ! -------------------------------------------------- CALL COMBI (PAI, PAIM, IM,JM, 1,R1,R2,RTG,NT,W2C) incr = 1 call outy (psea,pai,t,wv,u,v,im,jm,km,icnt,jcnt,incr, & 'before combi ') CALL COMBI (T, TM, IM,JM,KM,R1,R2,RTG,NT,W2C) incr = 1 call outy (psea,pai,t,wv,u,v,im,jm,km,icnt,jcnt,incr, & 'after combi ') CALL COMBI (WV, WVM, IM,JM,KM,R1,R2,RTG,NT,W2C) CALL COMBI (U, UM, IM,JM,KM,R1,R2,RTG,NT+10,W2C) CALL COMBI (V, VM, IM,JM,KM,R1,R2,RTG,NT+10,W2C) ! ****>< CALCULATE HGHT ><**** ! PRINT *,'##B4 HGHTCALC PSEAM AROUND ',CLAT,CLON CALL RLTLN (CI,CJ,CLAT,CLON,SLAT, & NPRO,DELS,SLON, XI,XJ,XLAT,XLON ) ! WRITE(6,3100) (I,I=CI-5,CI+5) DO 3203 J=CJ-7,CJ+7 ! WRITE(6,3101) J,(PSEAM(I,J),I=CI-5,CI+5) 3203 CONTINUE ! PRINT *,'##INVBAL HGHT' DO 888 K=1,KM DO 888 I=1,IM DO 888 J=1,JM W3C(I,J,K)=ZP(I,J,K) W3A(I,J,K)=U(I,J,K) W3B(I,J,K)=V(I,J,K) 888 CONTINUE !!$ w3? are all ok here. In original program, the j of w3? is inversed !!$ CALL INVBAL ( W3A, W3B, W3C, W2C, W2D, W2E, W2F, W2G, & !!$ IM, JM, KM,NPRO,DELS,SLON,XI,XJ,XLAT,XLON,SLAT ) CALL INVBAL ( W3A, W3B, W3C, W2C, W2D, W2E, W2F, W2G, & IM, JM, KM,NPRO,DELS,SLON,XI,XJ,XLAT,XLON ) ! w3c is stuffed (quad) here DO 981 K=1,KM DO 981 I=1,IM DO 981 J=1,JM W3A(I,J,K)=W3C(I,J,K) 981 CONTINUE ! w3a is stuffed (quad) here ! **** CALCULATING HGHT DEVIATIONS (FROM AXISYMMETRIC ): CALL DEVSYM ( W3A, CLAT, CLON, R2, RTG, W2C, W2D, IM, JM, KM, & 100, 104,NPRO,DELS,SLON,XI,XJ,XLAT,XLON,SLAT ) ! ** PUT SYMMETRIC IMPLANTED PSEA IN W2G DO 30 J = 1, JM DO 30 I = 1, IM W2G(I,J) = PSEA(I,J) 30 CONTINUE CALL COMBI(W2G,PSEAM,IM,JM, 1,R1,R2,RTG,NT,W2C) CALL PHICAL(W3B, T,PHIS,W2G,SPLA, IM*JM,KM) ! w3b is symmetric here ! w3a is stuffed (quad) at lowest level only DO 400 I=1,IM DO 400 J=1,JM DO 400 K=1,KM ! W3A(I,J,K) = 0. IF ( W3A(I,J,K) .GT. 10.0 ) THEN ! PRINT *,' ASYM Z,DZ ',I,J,K,W3B(I,J,K)/9.8,W3A(I,J,K) ENDIF ! W3C(I,J,K) = ALOG(SPLA(K)) ! *** ADDING INVBAL DEVIATIONS TO PHICAL (SPCONV) HGHTS: ZPM(I,J,K) = W3B(I,J,K)/9.8+W3A(I,J,K) 400 CONTINUE ! *** TRANSPLANT HGHT *** ! REMOVE BIAS IN HGHT DATA DO 720 K=1,KM ZERR = 0. AZP = 0. AZPM = 0. NAV = 0 DO 721 I=1,IM DO 721 J=1,JM IF(ABS(ZP(I,J,K) - ZPM(I,J,K))/ZP(I,J,K) .GT. 0.01)GOTO 721 NAV = NAV + 1 AZP = AZP + ZP(I,J,K) AZPM = AZPM + ZPM(I,J,K) 721 CONTINUE IF(NAV .GT. 500)AZP = AZP/NAV IF(NAV .GT. 500)AZPM = AZPM/NAV IF(NAV .GT. 500)ZERR = AZP - AZPM ! PRINT *,' MEAN Z ERR ',K,NAV,AZP,AZPM,ZERR DO 722 I=1,IM DO 722 J=1,JM 722 ZPM(I,J,K) = ZPM(I,J,K) + ZERR 720 CONTINUE !zpm(1) is oblate here CALL COMBI (ZP, ZPM, IM,JM,KM,R1,R2,RTG,NT+15,W2C) ! ** RECALCULATE SURFACE PRESSURE AND USE FOR DEVIATIONS ** CALL QQRELH ( W3A, WV, T, SPLA, IM, JM, KM, 'G/KG', 'HTOQ' ) incr = 1 call outy (psea,pai,t,wv,u,v,im,jm,km,icnt,jcnt,incr, & 'before vrtemp ') CALL VRTEMP (W3B, T,W3A, IM,JM,KM, 'G/KG', 'RTOV') incr = 1 call outy (psea,pai,t,wv,u,v,im,jm,km,icnt,jcnt,incr, & 'after vrtemp ') CALL CALCPS(W2C, W2G,W3B,ZP,SPLA(1),SPLA(2), IM,JM,KM) CALL DEVSYM ( W2C, CLAT, CLON, R2, RTG, W2D, W2E, IM, JM, 1, & 100, 104, NPRO,DELS,SLON,XI,XJ,XLAT,XLON,SLAT ) PRINT *,'##PSEA AROUND ',CLAT,CLON WRITE(6,3100) (I,I=CI-5,CI+5) DO 3205 J=CJ-7,CJ+7 WRITE(6,3101) J,(PSEA(I,J),I=CI-5,CI+5) 3205 CONTINUE PRINT *,'##PSEAM AROUND ',CLAT,CLON WRITE(6,3100) (I,I=CI-25,CI+25) 3100 FORMAT(' ',5I6,'|',I5,'|',I5,4I6) DO 3201 J=CJ-7,CJ+7 WRITE(6,3101) J,(PSEAM(I,J),I=CI-5,CI+5) 3101 FORMAT(I4,' ',11F6.0) 3201 CONTINUE PRINT *,'##HGHT BASED DEV. IN PSEAM ',CLAT,CLON WRITE(6,3100) (I,I=CI-5,CI+5) DO 3202 J=CJ-7,CJ+7 WRITE(6,3101) J,(W2C(I,J),I=CI-5,CI+5) 3202 CONTINUE DO 402 I=1,IM DO 402 J=1,JM ! *** ADDING CALCPS DEVIATIONS TO JMA MSLPS: PSEAM(I,J) = PSEAM(I,J)+W2C(I,J) 402 CONTINUE ! *** TRANSPLANT MSLP *** CALL COMBI (PSEA,PSEAM,IM,JM, 1,R1,R2,RTG,NT+20,W2C) INCR = 1 CALL OUTY (PSEA,PAI,T,WV,U,V,IM,JM,KM,ICNT,JCNT,INCR, & 'AFTER TRANSPLANTATION ') 1000 CONTINUE ! PRINT *,'TRANSPLANT ? ',TPC(NT) 1999 CONTINUE ! ============================================== ! >>>> <<<< ! >>>> TYPHOON LOOP END <<<< ! >>>> <<<< ! ============================================== ! ============================================== ! >>>> POST-PROCESS <<<< ! ============================================== ! ----------- ( RH --> WV , T --> TV ) ----- CALL QQRELH ( WV, WV, T, SPLA, IM, JM, KM, 'G/KG', 'HTOQ' ) ! CALL VRTEMP (T, T,WV, IM,JM,KM, 'G/KG', 'RTOV') ! --- PAI CORRECTION DUE TO TEMPERATURE MODIFICATION ---- ! CALL SETDPS (DPSEA,T(1,1,1),PHIS,PAI,SPLA(1),PTOP,IM,JM) ! DO 401 I=1,IM ! DO 401 J=1,JM ! PAI(I,J) = PSEA(I,J) - DPSEA(I,J) - PTOP ! 401 CONTINUE) RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE VRTEMP ( TV, T, WV, IM, JM, KM, WVFCT, ITEMP ) ! ITEMP = 'VTOR' : (VIRTUAL TEMP) TO (REAL TEMP) CONVERSION ! = 'RTOV' : (REAL TEMP) TO (VIRTUAL TEMP) CONVERSION ! WVFCT = 'G/KG' OR 'G/G ' DIMENSION TV(IM*JM,KM), T(IM*JM,KM), WV(IM*JM,KM) CHARACTER*4 ITEMP,WVFCT EPSL = 0.608 IF (WVFCT(3:3) .EQ. 'K') EPSL = EPSL * 0.001 ! --------------------------- IF (ITEMP .EQ. 'RTOV') THEN ! --------------------------- DO 210 K=1,KM DO 210 I=1,IM*JM TV(I,K) = T (I,K) * (1. + EPSL*WV(I,K)) 210 CONTINUE ! --------------------------- ELSE ! --------------------------- DO 220 K=1,KM DO 220 I=1,IM*JM T (I,K) = TV(I,K) / (1. + EPSL*WV(I,K)) 220 CONTINUE ! --------------------------- ENDIF ! --------------------------- RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE WIND3D ( U, V, IM, JM, KM, DELSPL, SPLA, RIJ, ANGIJ, & GRAD, T, PLG, KPM, IRM, DR, F0, & VT, VR, W1A, W1B, W1C, W1D, W1E, ISURF ) DIMENSION U(IM,JM,KM), V(IM,JM,KM) & ,RIJ(IM,JM), ANGIJ(IM,JM) & ,GRAD(IRM,KPM), T(IRM,KPM), PLG(KPM) & ,VT(IRM,KM), VR(IRM,KM) & ,W1A(KM), W1B(KM), W1C(KPM), W1D(IRM) & ,W1E(IRM), DELSPL(KM), SPLA(KM) ! +-------------------------------------------------- ! + 3-DIMENSIONAL WIND PROFILE + ! + DETERMINED FROM 2-D GRADIENT FORCE + ! + , TEMPERATURE + ! + CREARED BY T.IWASAKI 1985/10/01 + ! + REFORMED BY M.UENO 1987/12/10 + ! + + ! + + ! + ISURF 1 : INCLUDING PBL MODIFICATION + ! + 0 : ONLY GRADIENT WIND + ! +-------------------------------------------------- ! ***************************** ! 2-D GRADIENT WIND ( VT(I,K) ) ! ***************************** ! ! IF GRAD---->0 , THEN VT------>0 CPAI = 3.1416 GRAV = 9.8 RGAS = 287.04 RBYG = RGAS/GRAV ! MONITOR AT EVERY (LSTP) POINT LSTP = IRM/20 LPMX = 1+(20-1)*LSTP DO 30 K = 1, KM DO 30 I = 1, IRM VT(I,K) = 0.0 VR(I,K) = 0.0 30 CONTINUE DO 40 K = 1, KM DO 40 J = 1, JM DO 40 I = 1, IM U(I,J,K) = 0.0 V(I,J,K) = 0.0 40 CONTINUE DO 100 I=2,IRM DO 110 K=1,KM W1B(K) = ALOG(SPLA(K)) 110 CONTINUE DO 111 KP=1,KPM W1C(KP) = GRAD(I,KP) 111 CONTINUE CALL SPLINE(W1A,W1B,KM,W1C,PLG,KPM,2) R = (I-1)*DR FRIN = 1./(F0**2*R) DO 120 K=1,KM GRFR = 1.+4.*W1A(K)*FRIN IF (GRFR.LT.0.) GRFR=0. SQRTGR = SQRT(GRFR) ! GRADIENT WIND VR(I,K) = 0. VT(I,K) = 0.5*F0*R*(-1.+SQRTGR) 120 CONTINUE ! STORE FOR SURFACE PROCESS W1E(I) = W1A(1) 100 CONTINUE DO 125 K=1,KM ! WRITE(6, *) 'GRADIENT WIND AT LEVEL =',K ! WRITE(6, *) 'I POINTS =',LPMX,'interval=',LSTP ! WRITE(6,600) (VT(I,K),I=1,LPMX,LSTP) 125 CONTINUE 600 FORMAT((5X,20F6.1)) ! --------------------------------------------------------------------- ! --------------------------------------------------------------------- IF (ISURF .EQ. 0) GO TO 1000 ! --------------------------------------------------------------------- ! --------------------------------------------------------------------- ! ********************************************************************* ! CALCULATE TANGENTIAL & RADIAL WIND ( VT(I,1), VR(I,1) ) ! UNDER ! NEUTRAL STRATIFICATION ( LOGARITHMIC LAW ) IN THE CONSTANT FLUX LAYER ! ********************************************************************* ! ! (1) VR*D(VT)/D(R)+W*D(VT)/D(P)+F0*VR+VT*VR/R=0 ! (2) VR*D(VR)/D(R)+W*D(VR)/D(P)-F0*VT-VT**2/R+GRAD=0 ! .....ADVECTION TERMS IN THE FORMULAE ARE NEGLECTED..... H00= RBYG * T(1,1) PSIGM1 = SPLA(1) PSIGM2 = SPLA(2) DZ = H00 * LOG(PSIGM1/PSIGM2) ! 350. : ARBITRARY DZ = MAX( DZ, 350.) Z1 = H00 * LOG(1000./SPLA(1)) ! Z0 = 0.001 " ROUGHNESS PARM. ON THE SEA SURF. ! F10 = LOG(10./Z0)/LOG(Z1/Z0) F10 = 0.8 ! WRITE(6,800) F10, DZ 800 FORMAT( / ,' 10M-LEVEL WIND FACTOR =',F10.3, & 5X,'LOWEST-LAYER DEPTH =',F10.0, / ) VTMAX = 0. DO 200 I=1,IRM W1D(I) = VT(I,1) VTMAX = MAX(VTMAX,W1D(I)) 200 CONTINUE ITER = 0 ! FOR ITERATION CHECK VTMAX = 1.5 * VTMAX ! FOR ITERATION CONVERGENCE CITE = 0.5 ! ------ ITERATION START ------------------------------- 250 CONTINUE DO 201 I=2,IRM-1 ! --- VR --- R = (I-1)*DR ! 10M-LEVEL WIND V0 = F10 * SQRT(VT(I,1)**2+VR(I,1)**2) ! DRAG COEF. CDG = 0.0012 + 0.000025 * V0 FCD = CDG * V0 / DZ VR1 = VR(I,1) VR2 = - FCD * F10*VT(I,1) / (F0 + VT(I,1)/R) VR(I,1) = (VR1 + VR2) * CITE 201 CONTINUE ! --- VT --- DO 202 I = 2,IRM-1 R = (I-1)*DR V0 = F10 * SQRT(VT(I,1)**2+VR(I,1)**2) CDG = 0.0012 + 0.000025 * V0 FCD = CDG * V0 / DZ GRFR = 1. + 4. * (W1E(I) + FCD * F10*VR(I,1)) / (F0**2*R) IF (GRFR.LT.0.) GRFR=0. SQRTGR = SQRT(GRFR) VT1 = VT(I,1) VT2 = 0.5*F0*R*(-1.+SQRTGR) VT(I,1)= (VT1 + VT2) * CITE 202 CONTINUE ITER = ITER + 1 ! --- ITERATION CHECK --- DO 210 I=1,IRM IF(ABS(VT(I,1)).GT.VTMAX .OR. ABS(VR(I,1)).GT.VTMAX) THEN DO 211 II=1,IRM VT(II,1) = W1D(II) VR(II,1) = 0. 211 CONTINUE WRITE(6,604) 604 FORMAT(1H0, 5X,'ILL CONVERGENCE IN SUBR.((WIND3D))', & '-------------------------->GRADIENT WIND ADOPTED') ! ---------- GO TO 1000 ! ---------- ENDIF 210 CONTINUE IF (ITER.LT.15) GO TO 250 ! -------- ITERATION END ------------- ! ****************************************************** ! MODIFY ( VR(I,2),,VR(I,KMPBL) ; VT(I,2),,VT(I,KMPBL) ) ! UNDER ! APPROXIMATION OF EKMN SPIRAL ! ****************************************************** ! VISCOSITY COEF. (M**2/SEC.) " VCOE = 3. DDEKM = SQRT(2.*VCOE/F0) HHEKM = CPAI*DDEKM KPBL = 1 P1 = SPLA(1) ! WRITE(6,801) HHEKM 801 FORMAT( / ,' EKMN-LAYER HEIGHT =',F10.0, / ) DO 300 K=2,KM PK = SPLA(K) H00 = RBYG * (T(1,1)+T(1,K))*0.5 DZ = H00*LOG(P1/PK) ZK = DZ+Z1 IF (ZK .GT. HHEKM) GO TO 350 KPBL = K FAC= EXP(-DZ/DDEKM)*SIN(ZK/DDEKM) DO 301 I = 2,IRM VR(I,K) = VR(I,1)*FAC 301 CONTINUE ! WRITE(6,802) K,ZK,FAC 802 FORMAT( / ,' IN-EKMN LEVEL =', I3,' HEIGHT =',F10.0, & 5X,'WIND FACTOR =',F10.3, / ) 300 CONTINUE 350 CONTINUE ! ***************************** ! OUTFLOW NEAR TROPOPAUSE LEVEL ! ***************************** ! --- DETERMINE THE VERTICAL WEIGHT ---- ! EXTENT (M) OF OUTFLOW LAYER DZOUT = 4000. ! DEPTH (M) BETWEEN TROPO. AND MAX OUTFLOW LEVEL DEP = 2000. TTROPO = T(IRM,1) DO 400 KP=1,KPM IF(T(IRM,KP).LT.TTROPO) THEN TTROPO = T(IRM,KP) PTROPO = EXP(PLG(KP)) ENDIF 400 CONTINUE CSUM = 0. DO 50 K = 1, KM W1A(K) = 0.0 50 CONTINUE DO 401 K=KPBL+1,KM P = SPLA(K) H00 = RBYG * (TTROPO+T(IRM,K))*0.5 DZ = ABS ( H00 * LOG(P/PTROPO) - DEP ) ! WEIGHTING FUNC. W1A(K) = EXP(-(DZ/DZOUT)**2.) CSUM = W1A(K)*DELSPL(K)/1000. + CSUM 401 CONTINUE DO 402 K=KPBL+1,KM W1A(K) = W1A(K)/CSUM 402 CONTINUE ! WRITE(6,803) KPBL,PTROPO,TTROPO 803 FORMAT( / ,' (KPBL,PTROPO,TTROPO) =',I3,5X,F10.0,5X,F10.0, / ) ! WRITE(6,804) W1A 804 FORMAT( / ,' OUTFLOW WEIGHT AT EACH LEVEL =',12F10.3, / ) DO 410 IR =2,IRM FLXMSS = 0 DO 411 K=1,KPBL FLXMSS = DELSPL(K)/1000.*VR(IR,K) + FLXMSS 411 CONTINUE DO 412 K=KPBL+1,KM VR(IR,K) = -FLXMSS*W1A(K) 412 CONTINUE 410 CONTINUE DO 420 K=1,KM ! WRITE(6, *) 'TANGENTIAL WIND AT LEVEL =',K ! WRITE(6,600) (VT(I,K),I=1,LPMX,LSTP) 420 CONTINUE DO 425 K=1,KM ! WRITE(6, *) 'RADIAL WIND AT LEVEL =',K ! WRITE(6,600) (VR(I,K),I=1,LPMX,LSTP) 425 CONTINUE !---------------------------------------------------------------------- !---------------------------------------------------------------------- 1000 CONTINUE !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! **************************************************** ! INTERPOLATION OF ( VT(I,K), VR(I,K) ) INTO 3-D SPACE ! **************************************************** RAD = ACOS(0.)/90. DRIN = 1./DR isp=118 jsp=30 DO 500 I=1,IM DO 500 J=1,JM FR = RIJ(I,J)*DRIN + 1 IR = FR DDR = FR - IR IF(IR.GE.IRM) GO TO 500 RADANG = RAD * ANGIJ(I,J) ANGCOS = COS(RADANG) ANGSIN = SIN(RADANG) ad=90. ad1=ad*3.1415926/180. DO 501 K=1,KM VRG = VR(IR,K) + (VR(IR+1,K)-VR(IR,K))*DDR VTG = VT(IR,K) + (VT(IR+1,K)-VT(IR,K))*DDR U(I,J,K) = -VRG*ANGCOS + VTG*ANGSIN V(I,J,K) = -VRG*ANGSIN - VTG*ANGCOS 501 CONTINUE 500 CONTINUE RETURN END !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- SUBROUTINE XYTOLL ( RLAT, RLON, FI, FJ,SLAT, & NPROJC, DELS, SLON, XI, XJ, XLAT, XLON ) !********************************************************************** ! LATITUDE AND LONGITUDE OF GRID POINT ARE CALUCULATED * ! * !OUT-PUT RLAT LATITUDE (DEGREE) -90.< < 90. * ! RLON LONGITUDE (DEGREE) 0.< <360. * !IN-PUT DELS GRID LENGTH ( M) * ! | IF NPROJC='LL ' (DEG) * ! | NPROJC MAP PROJECTION * ! | 'PSN ' POLAR STEREO SLAT=60 N NORTH * ! | 'PSS ' POLAR STEREO SLAT=-60 N SOUTH * ! | 'MER ' MERCATOR SLAT=0. N * ! | 'LMN ' LAMBERT SLAT=30 N , 60 N NORTH * ! | 'LMS ' LAMBERT SLAT=-30 N ,-60 N SOUTH * ! | 'LL ' LATITUDE-LONGITUDE * ! Y-COORDINATE 0.< <360. * ! | SLON STANDARD LONGITUDE ! | * ! | (XI,XJ) <----------> (XLAT,XLON) * ! STANDARD POINT * !********************************************************************** CHARACTER * 4 NPROJC ! RADIUS OF EARTH" A1=6371.E3 PI=3.14159 POI=PI /180. RPOI=180./ PI !-----------------------------------POLAR STEREO----------------------- IF(NPROJC(1:3).EQ.'PSN') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=60. SLAT1=SLAT*POI SLON1=SLON*POI XLAT1=XLAT*POI XLON1=XLON*POI RL0=A1*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PI-XLAT1)) X0=RL0*SIN(XLON1-SLON1) Y0=RL0*COS(XLON1-SLON1) XP=X0+(FI-XI)*DELS YP=Y0+(FJ-XJ)*DELS RLAT=90.-RPOI*2.*ATAN(SQRT(XP*XP+YP*YP)/(A1*(1.+SIN(SLAT1)))) RLON=SLON+RPOI*ATAN2(XP,YP) ENDIF IF(NPROJC(1:3).EQ.'PSS') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLAT=-60. SLAT1=-SLAT*POI SLON1= SLON*POI XLAT1=-XLAT*POI XLON1= XLON*POI RL0=A1*(1.+SIN(SLAT1))*TAN(0.5*(0.5*PI-XLAT1)) X0=RL0*SIN(XLON1-SLON1) Y0=RL0*COS(XLON1-SLON1) XP=X0+(FI-XI)*DELS YP=Y0+(XJ-FJ)*DELS RLAT=-90.+RPOI*2.*ATAN(SQRT(XP*XP+YP*YP)/(A1*(1.+SIN(SLAT1)))) RLON=SLON+RPOI*ATAN2(XP,YP) ENDIF !-----------------------------------MERCATOR--------------------------- IF(NPROJC(1:3).EQ.'MER') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" ! NOTE: AGAIN SET SLAT ! SLAT=0. ! SLAT=42.3687 SLAT1=SLAT*POI SLON1=SLON*POI XLAT1=XLAT*POI XLON1=XLON*POI AC=A1*COS(SLAT1) RLON=XLON+RPOI*(FI-XI)*DELS/AC T1=EXP((XJ-FJ)*DELS/AC)*(1.+SIN(XLAT1))/COS(XLAT1) RLAT=RPOI*ASIN((T1*T1-1.)/(T1*T1+1.)) ENDIF !-----------------------------------LAMBERT---------------------------- IF(NPROJC(1:3).EQ.'LMN') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLATA=30. SLATB=60. PI4=PI*0.25 SLATA1=SLATA*POI SLATB1=SLATB*POI SLAT1=PI4-SLATA*POI*0.5 SLAT2=PI4-SLATB*POI*0.5 SLON1=SLON*POI XLAT1=XLAT*POI XLON1=XLON*POI CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) RCK=1./CK ACN=A1*COS(SLATA1)*RCK RL0=ACN*(TAN(PI4-XLAT1*0.5)/TAN(SLAT1))**CK ! ---------------------------------------- XSLON = XLON - SLON XSLON = MOD (XSLON+900., 360.) - 180. XSLON = XSLON * POI * CK X0=RL0*SIN(XSLON) Y0=RL0*COS(XSLON) ! ---------------------------------------- ! ---- X0=RL0*SIN((XLON1-SLON1)*CK) ! ---- Y0=RL0*COS((XLON1-SLON1)*CK) XP=X0+(FI-XI)*DELS YP=Y0+(FJ-XJ)*DELS RLAT=90.-2.*RPOI*ATAN((SQRT(XP*XP+YP*YP)/ACN)**RCK*TAN(SLAT1)) RLON=SLON+RCK*RPOI*ATAN2(XP,YP) ! ---------------------------------------- RRLON = RLON RRLON = RRLON - XLON RRLON = MOD (RRLON+900., 360.) - 180. RLON = RRLON + XLON ! ---------------------------------------- ENDIF IF(NPROJC(1:3).EQ.'LMS') THEN ! STANDARD LATITUDE ;MAP FACTOR=1" SLATA=-30. SLATB=-60. PI4=PI*0.25 SLATA1=-SLATA*POI SLATB1=-SLATB*POI SLAT1=PI4+SLATA*POI*0.5 SLAT2=PI4+SLATB*POI*0.5 SLON1=SLON*POI XLAT1=-XLAT*POI XLON1=XLON*POI CK=LOG(COS(SLATA1)/COS(SLATB1))/LOG(TAN(SLAT1)/TAN(SLAT2)) RCK=1./CK ACN=A1*COS(SLATA1)*RCK RL0=ACN*(TAN(PI4-XLAT1*0.5)/TAN(SLAT1))**CK ! ---------------------------------------- XSLON = XLON - SLON XSLON = MOD (XSLON+900., 360.) - 180. XSLON = XSLON * POI * CK X0=RL0*SIN(XSLON) Y0=RL0*COS(XSLON) ! ---------------------------------------- ! ---- X0=RL0*SIN((XLON1-SLON1)*CK) ! ---- Y0=RL0*COS((XLON1-SLON1)*CK) XP=X0+(FI-XI)*DELS YP=Y0+(XJ-FJ)*DELS RLAT=-90.+2.*RPOI*ATAN((SQRT(XP*XP+YP*YP)/ACN)**RCK*TAN(SLAT1)) RLON=SLON+RCK*RPOI*ATAN2(XP,YP) ! ---------------------------------------- RRLON = RLON RRLON = RRLON - XLON RRLON = MOD (RRLON+900., 360.) - 180. RLON = RRLON + XLON ! ---------------------------------------- ENDIF !-----------------------------------LATITUDE LONGITUDE----------------- IF(NPROJC(1:2).EQ.'LL') THEN ! RLON = (FI-1)*(DELS/111000.) ! RLAT = 90.-(FJ-1)*(DELS/111000.) RLON = XLON + (FI-XI)*(DELS/111000.) RLAT = XLAT - (FJ-XJ)*(DELS/111000.) ENDIF !---------------------------------------------------------------------- IF(RLAT .GT. 90.) RLAT = 90. IF(RLAT .LT.-90.) RLAT =-90. DO 100 ITR=1,2 IF(RLON .LT. 0.) RLON =RLON + 360. IF(RLON .GT.360.) RLON =RLON - 360. 100 CONTINUE RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE OPTMUM(IDIM,JDIM,ALFA,SIGMA) ! THIS ROUTINE OPTIMISES THE OVER-RELAXATION COEFFICIENT ! FOR LIEBH PARAMETER ( PI=3.14159265358979 ) T = ( 4.0 + SIGMA ) / ( COS ( PI / IDIM ) + COS ( PI / JDIM ) ) C = 0.5 * T * T - 1.0 IF ( C .LE. 100.0 ) THEN V = C - SQRT ( C * C - 1 ) ELSE IF ( C .LE. 1.0E10 ) THEN T2 = 0.125 / ( C * C * C ) ELSE T2 = 0.0 ENDIF V = 0.5 / C + T2 ENDIF ALFA = ( 1.0 + V ) / ( 4.0 + SIGMA ) RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- SUBROUTINE TIMDAT ( ITIM, IDAT, IT, NTIM, NDAT ) ! THIS SUBROUTINE DETERMINES THE TIME AND DATE (NTIM,NDAT) AT 'IT' HOURS ! FROM (ITIM,IDAT). 'IT' MAY BE POSITIVE, ZERO OR NEGATIVE ! TIMES ARE IN UNITS OF HOURS DIMENSION IDPM(12) DATA IDPM / 31, 28, 31, 30, 31, 30, 31, 31, 30, 31, 30, 31 / IDY = MOD ( IDAT, 100 ) MTH = MOD ( IDAT, 10000 ) / 100 IYR = IDAT / 10000 ID = IABS ( IT ) / 24 IF ( IT .EQ. 0 ) THEN ! NO CHANGE TO DATE/TIME NTIM = ITIM NDAT = IDAT RETURN ELSEIF ( IT .LT. 0 ) THEN ! DECREMENT DATE/TIME BY IT HOURS IX = MOD ( -IT, 24 ) NTIM = ITIM - IX * 100 IF ( NTIM .LT. 0 ) THEN NTIM = NTIM + 2400 ID = ID + 1 ENDIF IF ( ID .NE. 0 ) THEN DO 100 I = 1, ID IDY = IDY - 1 IF ( IDY .EQ. 0 ) THEN MTH = MTH - 1 IF ( MTH .EQ. 0 ) THEN MTH = 12 IYR = IYR - 1 ENDIF IDY = IDPM ( MTH ) IF ( MOD ( IYR, 4 ) .EQ. 0 .AND. MTH .EQ. 2 ) IDY = 29 ENDIF 100 CONTINUE ENDIF ELSEIF ( IT .GT. 0 ) THEN ! INCREMENT DATE/TIME BY IT HOURS IX = MOD ( IT, 24 ) NTIM = ITIM + IX * 100 IF ( NTIM .GE. 2400 ) THEN NTIM = NTIM - 2400 ID = ID + 1 ENDIF IF ( ID .NE. 0 ) THEN DO 180 I = 1, ID IDY = IDY + 1 IF ( IDY .GT. IDPM ( MTH ) .AND. & ( MOD ( IYR, 4 ) .NE. 0 .OR. MTH .NE. 2 .OR. & IDY .NE. 29 ) ) THEN IDY = 1 MTH = MTH + 1 IF ( MTH .EQ. 13 ) THEN MTH = 1 IYR = IYR + 1 ENDIF ENDIF 180 CONTINUE ENDIF ENDIF IF ( IYR .LT. 0 ) IYR = IYR + 100 IF ( IYR .GT. 99 ) IYR = IYR - 100 NDAT = IYR * 10000 + MTH * 100 + IDY RETURN END !---------------------------------------------------------------------- !********************************************************************** !---------------------------------------------------------------------- !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! subroutine unicon(pres,iyear,imonth,iday,ihr, & aslat,anlat,awlong,aelong,tg,ug, & vg,zga,qg,nxx,nyy,nplevs) ! sub now writes out bogussed analyses in NCEP Reanalysis pressure ! level format ! include 'dim.h' real tg(nxx,nyy,nplevs),ug(nxx,nyy,nplevs),vg(nxx,nyy,nplevs) real zga(nxx,nyy,nplevs),qg(nxx,nyy,nplevs) real temp(nxx,nyy,nplevs) real zga1(nxx,nyy),qg1(nxx,nyy) character*4 grid,name integer pres(nplevs) ! Open the output files for bogussed data ! open(70,file='geohgt.out') ! open(71,file='airtemp.out') ! open(72,file='spechum.out') ! open(73,file='uwind.out') ! open(74,file='vwind.out') ! Read in the data, reversing the JMA grid specification needed ! for this program, back to FST format do k=1,nplevs do j=1,nyy do i=1,nxx write(70,*)zga(i,nyy-j+1,k) temp(i,j,k)=zga(i,nyy-j+1,k) write(71,*)tg(i,nyy-j+1,k) write(72,*)qg(i,nyy-j+1,k) write(73,*)ug(i,nyy-j+1,k) write(74,*)vg(i,nyy-j+1,k) enddo enddo enddo ! Close the output files close(70) close(71) close(72) close(73) close(74) return end !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- subroutine clatlon(deglat,deglon,ni,nj,clat,clon,xais,d60) !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc ! Given x, y ,to calculate the latitude and longitude of grid. ! In Polar sterographic projection. ! !CCc --------------------------------------------------------- real deglat(ni,nj),deglon(ni,nj) a=6.37122e+6 pi=4.*atan(1.) cnt=180./pi phi0=60. smap0=(1.+sin(phi0/cnt))/(1.+sin(clat/cnt)) g=smap0*a*cos(clat/cnt)*cos(clon/cnt) h=smap0*a*cos(clat/cnt)*sin(clon/cnt) ! print*,'g= ,h= ',g,h alpha=360.-xais do j=1,nj do i=1,ni !cccccccc !! 1. dx=(i-(ni+1)/2)*d60 dy=(j-(nj+1)/2)*d60 !ccc !! 2. x=g+dx*cos(alpha/cnt)-dy*sin(alpha/cnt) y=h+dx*sin(alpha/cnt)+dy*cos(alpha/cnt) !cc !c 3. g1=sqrt((x/a)**2+(y/a)**2)/(1.+sin(phi0/cnt)) deglat(i,j)=2.*atan((1.-g1)/(1.+g1)) if(x.eq.0.0.and.y.eq.0.) then deglon(i,j)=0. else if(x.eq.0.) then if(y.gt.0.) deglon(i,j)=pi/2. if(y.lt.0.) deglon(i,j)=2.*pi-pi/2. else if(y.eq.0.) then if(x.gt.0.) deglon(i,j)=0. if(x.lt.0.) deglon(i,j)=pi else f=abs(y/x) if(x.gt.0.0.and.y.gt.0.) then deglon(i,j)=atan(f) else if(x.gt.0..and.y.lt.0.) then deglon(i,j)=2.*pi -atan(f) else if(x.lt.0..and.y.gt.0.) then deglon(i,j)=pi-atan(f) else deglon(i,j)=pi+atan(f) endif endif aa=deglat(i,j)*180./3.1415 bb=deglon(i,j)*180./3.1415 ! if(i.eq.1)print*,'lat',aa,'lon',bb enddo enddo do i=1,ni do j=1,nj deglat(i,j)=deglat(i,j)*180./3.1415 deglon(i,j)=deglon(i,j)*180./3.1415 enddo enddo return end !--------------------------------------------------------- !********************************************************* !--------------------------------------------------------- subroutine fij(Gii0,Gjj0,lat1,long1) real lat1,long1 data pi/3.1415926/ cnt=pi/180. clat=39.*cnt clon=-60.9*cnt phi0=60. xais=330.9 nx=207 ny=185 alpha=360-xais alpha=alpha*cnt x0=fx(clat,clon) y0=fy(clat,clon) dx=30. dy=30. call xycoord(x1,y1,lat1,long1,x0,y0,alpha) GII0=(x1/dx +float(nx+1)/2.) GJJ0=(y1/dy+float(ny+1)/2.) return end !---------------------------------------------------------- !********************************************************** !---------------------------------------------------------- real function fx(lat,lon) real lat,lon real m a=6371.22 phi0=60. pi=3.1415926 rad=pi/180. phi0=phi0*rad m=(1.+sin(phi0))/(1.+sin(lat)) fx=m*a*cos(lat)*cos(lon) return end !---------------------------------------------------------- !********************************************************** !---------------------------------------------------------- real function fy(lat,lon) real lat,lon real m a=6371.22 phi0=60. pi=3.1415926 rad=pi/180. phi0=phi0*rad m=(1.+sin(phi0))/(1.+sin(lat)) fy=m*a*cos(lat)*sin(lon) return end !---------------------------------------------------------- !********************************************************** !---------------------------------------------------------- subroutine xycoord(x,y,lat,lon,x0,y0,alpha) real lat,lon x1=fx(lat,lon) y1=fy(lat,lon) x1=x1-x0 y1=y1-y0 x=x1*cos(alpha)+y1*sin(alpha) y=-x1*sin(alpha)+y1*cos(alpha) return end !-------------------------------------------------------------------- !******************************************************************** !-------------------------------------------------------------------- subroutine clatlon1(deglat,deglon,fi,fj,clat,clon,xais,d60) !CCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCCc ! Given x, y ,to calculate the latitude and longitude of grid. ! In Polar sterographic projection. ! !!! --------------------------------------------------------- ! real deglat(ni,nj),deglon(ni,nj) ni=207 nj=185 a=6.37122e+6 pi=4.*atan(1.) cnt=180./pi phi0=60. smap0=(1.+sin(phi0/cnt))/(1.+sin(clat/cnt)) g=smap0*a*cos(clat/cnt)*cos(clon/cnt) h=smap0*a*cos(clat/cnt)*sin(clon/cnt) ! print*,'g= ,h= ',g,h alpha=360.-xais !cccccccc ! 1. dx=(fi-(ni+1)/2)*d60 dy=(fj-(nj+1)/2)*d60 !ccc ! 2. x=g+dx*cos(alpha/cnt)-dy*sin(alpha/cnt) y=h+dx*sin(alpha/cnt)+dy*cos(alpha/cnt) !cc !c n. g1=sqrt((x/a)**2+(y/a)**2)/(1.+sin(phi0/cnt)) deglat=2.*atan((1.-g1)/(1.+g1)) if(x.eq.0.0.and.y.eq.0.) then deglon=0. else if(x.eq.0.) then if(y.gt.0.) deglon=pi/2. if(y.lt.0.) deglon=2.*pi-pi/2. else if(y.eq.0.) then if(x.gt.0.) deglon=0. if(x.lt.0.) deglon=pi else f=abs(y/x) if(x.gt.0.0.and.y.gt.0.) then deglon=atan(f) else if(x.gt.0..and.y.lt.0.) then deglon=2.*pi -atan(f) else if(x.lt.0..and.y.gt.0.) then deglon=pi-atan(f) else deglon=pi+atan(f) endif endif deglat=deglat*180./3.1415 deglon=deglon*180./3.1415 return end !---------------------------------------------------------- !********************************************************** !---------------------------------------------------------- real function epsilonR() ! This function returns machine epsilon for single- ! precision real values. real d1 d1=1. do while( (1.+d1) .gt. 1. ) d1=d1*0.5 enddo epsilonR=d1*2. return end !---------------------------------------------------------- !********************************************************** !----------------------------------------------------------