PROGRAM BAMADV C C** THIS IS THE Beta-Advection model C C** MAIN PROGRAM DOCUMENTATION BLOCK C C MAIN PROGRAM: BAM DEEP (BETA-ADVECTION MODEL) C PRGMMR: MARKS ORG: NMC421 DATE:88-01-22 C GROSS WD80JG 90-01-01 C C ABSTRACT: PRODUCES A HURRICANE TRACK FORECAST BASED ON THE C WORK OF G. HOLLAND (1983). USES THE LARGE-SCALE FIELDS FROM C THE SPECTRAL MODEL AND SOLVES THE DIVERGENT, BAROTROPIC C VORTICITY EQUATION ON A BETA PLANE. C C PROGRAM HISTORY LOG: C 88-01-14 DON MARKS FIX S.H. FORECASTS AND DELETE C SOME EXTRANEOUS ERROR MSGS C 88-02-18 DON MARKS ADD CHECK FOR VERTICAL SHEAR IN SUBR GETLVL C 88-05-02 DON MARKS ADD PARAMETER STMTS TO USE GRIDDED FLDS C 88-05-09 DON MARKS CHANGE FT16F001 TO TOP-DOWN STACK C 88-11-03 DON MARKS CHANGE CALC OF C FOR S.H. STORMS C 89-05-10 DON MARKS INTERPOLATION OF MEAN WINDS OVER LARGER AREA C 89-07-05 DON MARKS NEW DATE-TIME MESSAGE C 90-05-31 JIM GROSS ATCF CHANGES C 90-06-07 JIM GROSS ADDITION OF 84 HOUR FORECAST C 91-05-22 JIM GROSS REMOVED VERTICAL SHEAR CHECK C 92-07-30 JIM GROSS CLEAN UP ROUTINE AND INSERT COMMENTS C 96-04-00 JIM GROSS REVISED FOR THE CRAY AND TO USE GRIB C 97-06-23 JIM GROSS REVISED AVIATION FIELDS ACCESS TO HANDLE C FOUR RUNS PER DAY TO 78 HOURS C 2012 C. SISKO REVISED FOR 168 HR FORECAST C C USAGE: C C INPUT FILES: C C OUTPUT FILES: C FT06F001 - OUTPUT PRINT FILE C C SUBPROGRAMS CALLED: C UNIQUE: - STHETA,RADIUS,VARIN,GETLVL C C LIBRARY: C W3LIB - C C EXIT STATES: C COND = 0 - SUCCESSFUL RUN C C C ATTRIBUTES: C LANGUAGE: CF77 FORTRAN C C REMARKS: THIS PROGRAM UNIT CONTAINS THE STATEMENTS WHICH DETERMINE C WHICH BAM MODEL IS WRITTEN TO THE ATCF FILES C include 'dataformats.inc' c CAS PARAMETER ( IPTS = 360*181, KF = 4, ntau = 10 ) PARAMETER ( IPTS = 360*181, KF = 4, ntau = 14 ) C COMMON /OBSV/ BETA,GAMMA,REFF,VS,XI,VB,ALPHA,GAMMAH,C, & SPD,VNB,VEB,XLATC COMMON /MNWIND/ UCOMP(IPTS),UTEND(IPTS),VCOMP(IPTS),VTEND(IPTS), & KGDSI(22),MXWAVE C CAS DIMENSION jNLEVS(7),LEVS(6),NLEVS(6),ibam(11,3) DIMENSION jNLEVS(7),LEVS(6),NLEVS(6),ibam(15,3) REAL THINIT(8),SPDINI(8),FPBAM(2,ntau),URAD(KF),VRAD(KF) REAL RLAT(KF),RLON(KF),CROT(KF),SROT(KF) INTEGER MBLEV(8),INLEVS(7),IPOPT(20),KGDSO(22) LOGICAL LGI(IPTS),LGO(KF),YORN C CHARACTER*1 EW CHARACTER*2 CYCLE CHARACTER*3 FTIME(ntau) character*4 bam_type CHARACTER*8 STRMID, YRMNDY CHARACTER*10 aYMDH,mYMDH character*50 input_file CHARACTER*80 FLNAME C type ( AID_DATA ) comRcd, tauData c DATA FTIME / '012', '024', '036', '048', '060', & '072', '084', '096', '108', '120', & '132', '144', '156', '168' / DATA MBLEV /100,200,300,400,500,700,850,1000/ DATA NLEVS /850,700,500,400,300,200/ C C** ZERO ARRAYS AND INITIALIZE CONSTANTS C DO I = 1, ntau FPBAM( 1, I ) = 0.0 FPBAM( 2, I ) = 0.0 enddo c do i = 1, ntau + 1 do j = 1, 3 ibam( i, j ) = 0 enddo enddo C IFP = 0 IHR = 0 IREFF = 0 ISTART = 0 GAMMA = 0.0 c c** Get the command line arguments c call getarg ( 1, bam_type ) call getarg ( 2, input_file) c strmid = input_file(1:8) c c** Open the input file c open ( 11, file=input_file, status='old', iostat=ios, err=1010 ) c c** Read in the compute data and close the file c call getARecord ( 11, "CARQ", comRcd, istat ) if ( istat .eq. 0 ) goto 1020 c close ( 11 ) c c** Read in the current data c call getSingleTAU ( comRcd, 0, tauData, istat ) if ( istat .ne. 1 ) goto 1030 c aymdh = tauData%aRecord(1)%DTG xlat = tauData%aRecord(1)%lat xlon = tauData%aRecord(1)%lon drn = tauData%aRecord(1)%dir speed = tauData%aRecord(1)%speed c if (tauData%aRecord(1)%windcode .eq. 'AAA' ) then rd34ne = tauData%aRecord(1)%radii(1) rd34se = tauData%aRecord(1)%radii(1) rd34sw = tauData%aRecord(1)%radii(1) rd34nw = tauData%aRecord(1)%radii(1) else rd34ne = tauData%aRecord(1)%radii(1) rd34se = tauData%aRecord(1)%radii(2) rd34sw = tauData%aRecord(1)%radii(3) rd34nw = tauData%aRecord(1)%radii(4) endif C C** GET THE INPUT PARAMETERS C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C V2 - WIND (M/S) AT OUTER EDGE OF STORM C R2 - RADIUS OF V2 C XLAT2 - LATITUDE OF STORM C XLON2 - LONGITUDE OF STORM C JHR - HOURS PAST FORECAST FILE (<12) AT WHICH TO START C DRN - DIRECTION OF STORM MOTION C SPD - SPEED OF STORM MOTION C LVLS - UP TO 6 VALUES (200,300,ETC) OF LEVELS TO BE USED C LHR - HOURS PAST INITIAL FORECAST DATA TIME C C INPUT FILES: (DELETE IF NO INPUT FILES IN SUBPROGRAM) C FT05F001 - INPUT PARAMETERS IN NAMELIST NAMED INPUT C - XLAT,XLON -- POSITION OF STORM C - IHR -- HOUR OF OBSERVATION C - DIR,SPEED --CURRENT DIRECTION&SPEED OF STORM C - OUTWND,RADIi -- SPEED OF OUTER WIND & ITS RADIUS C - LEVS -- STEERING LEVELS DESIRED c OUTWND = 34.0 radii = (RD34NE + RD34SE + RD34SW + RD34NW)/4.0 C C** IF THE RADIUS OF 34 KNOT WINDS IS LESS THAN OR EQUAL TO ZERO, C** SET WINDS TO 25 KNOTS AND RADIUS TO 50 NM. C IF ( radii .LE. 0.0 ) THEN OUTWND = 25.0 radii = 50.0 WRITE (*,'(/)') WRITE (*,*) ' THE RADIUS OF 34 KNOT WINDS IS LESS THAN ZERO' WRITE (*,*) ' SO THE OUTER WINDS WILL BE SET TO 25 KNOTS' WRITE (*,*) ' AND THE RADIUS TO 50 NM' WRITE (*,'(/)') ENDIF C C** ALLOW FOR SOUTHERN HEMISPHERE STORMS C IF ( XLAT .LT. 0.0 ) OUTWND = - OUTWND C V2 = OUTWND/1.94254 R2 = radii*111137.0/60.0 C XLAT2 = XLAT XLON2 = XLON IF (XLON2.LT.0.0) XLON2 = XLON2 + 360.0 C Dir = Drn SPD = SPEED*0.515 Dir = 360.0 - Dir C if ( bam_type .eq. 'BAMD' ) then c NLEV = 6 ! DEEP BAM (850MB TO 200MB) C elseif ( bam_type .eq. 'BAMM' ) then c NLEV = 4 ! MEDIUM BAM (850MB TO 400MB) C elseif ( bam_type .eq. 'BAMS' ) then c NLEV = 2 ! SHALLOW BAM (850MB TO 700MB) c else c stop ' ERROR - The type of BAM has not been defined' c endif C DO I = 1,6 C IF (I.LE.NLEV) THEN LEVS(I) = NLEVS(I) ELSE LEVS(I) = 0 ENDIF C jNLEVS(I) = LEVS(I) C enddo C jNLEVS(7) = 0 LVLf = 0 C DO 10 I = 1,6 C J = 7 - I LEV = jNLEVS(J)/100 IF (LEV.LE.1) GO TO 10 IF (LEV.LT.7) LVLf = LVLf*10 + LEV IF (LEV.EQ.7) LVLf = LVLf*10 + 6 IF (LEV.EQ.8) LVLf = LVLf*10 + 7 C 10 continue C WRITE (*,'(//,'' BAM FORECAST MODEL'',/)') WRITE (*,'('' INPUT PARAMETERS FOR TROPICAL CYCLONE '',A8, & '' ON YMDH = '',A10,/)') STRMID,aYMDH WRITE (*,'('' LA0 = '',F5.1,'' N LO0 = '',F6.1, & '' W DIR0 = '',F4.0,'' DEG SPD0 = '',F4.0,'' KT'',/, & '' AVERAGE RADIUS OF 34 KT WINDS = '',F5.1,'' NM'',//)') & XLAT,XLON,DRn,SPEED,RADIi C WRITE (6,1) XLAT,XLON,Drn,SPEED,OUTWND,RADIi,LEVS 1 FORMAT (/,' XLAT=',F5.1,' XLON=',F6.1,' DIR=',F4.0, & ' SPEED=',F4.1,' OUTWND=',F4.1,' RADIUS=',F4.0,6X,/, & 1X,'LEVS =',6I8) c c** end varin subroutine C C** CHECK STATUS OF THE CURRENT 78-HOUR AVIATION FORECAST GRIB FILE. C** THAT'S A 6 HOUR LAG. IF IT DOES NOT EXIST, CHECK THE EXISTENCE C** OF THE 78-HOUR FILE OF THE PREVIOUS AVIATION RUN. THAT'S A C** 12 HOUR LAG. IF THAT FILE DOES NOT EXIST, QUIT. C CALL dtgmod ( aYMDH, -6, mYMDH, istat ) C YRMNDY = mYMDH(1: 8) CYCLE = mYMDH(9:10) C if ( aymdh(9:10) .eq. '00' .or. aymdh(9:10) .eq. '06' .or. & aymdh(9:10) .eq. '12' .or. aymdh(9:10) .eq. '18' ) then c CAS FLNAME = '/com/gfs/prod/gfs.'//YRMNDY//'/gfs.t'//CYCLE// CAS & 'z.pgrbf126' FLNAME = '/tmp/ARK/TEST-DATA/gfs.'//YRMNDY// & '/gfs.t'//CYCLE// & 'z.pgrbf174' print *,FLNAME INQUIRE ( FILE=FLNAME, EXIST=YORN, IOSTAT=IOS, ERR=1040 ) C IF ( YORN ) THEN c WRITE (*,*) ' CURRENT AVIATION GRIB AVAILABLE, USING ', & FLNAME(1:40) IHR = 6 CAS IHREND = 126 IHREND = 174 IFHR = 6 c ELSE STOP ' CURRENT AVIATION RUN FILES UNAVAILABLE' ENDIF c ELSE STOP ' Current hour not equal to 00, 06, 12 or 18!' ENDIF C PRINT 2, mYMDH(1:4),mYMDH(5:6),mYMDH(7:8),mYMDH(9:10) 2 FORMAT (/,' DATE/TIME OF DATA IS ',A4,'/',A2,'/',A2,'/',A2, & ' FROM THE AVN MODEL',/) C INITHR = IHR XI = 0.5 C = V2*(R2**XI) C = (ABS(C) + 7460.0)/2.0 C IF (XLAT2.LT.0.0) C = - C C C** SET THE SPECTRAL FILTERING WAVE NUMBER MXWAVE = 25 C C** SET R34 = RADIUS OF 34 KNOT WIND, SFAC = DIAMETER IN DEGREES C R34 = (V2**2/17.5**2)*R2 SFAC = 2.0*R34/111137.0 C C** USE FIRST GUESS FOR EFFECTIVE RADIUS (REFF) AS 280000.0 KM C REFF = 280000.0 VS = C/(REFF**XI) REDEG = AMAX1(REFF,1.3*R2,4.0E5)/111137.0 C MXWAVE = 360.0/(4.0*AMAX1(REDEG,SFAC)) C cc PRINT *, ' # R2 = ',R2,' SFAC = ',SFAC,' REFF = ',REFF cc PRINT *, ' ## REDEG = ',REDEG,' MXWAVE = ',MXWAVE C LVLS = 2 C 20 CONTINUE cc PRINT *, ' * GOT TO 20 CONTINUE' C THETAM = DIR/57.3 XLATC = XLAT2 XLONC = XLON2 C 40 CONTINUE cc PRINT *, ' *** GOT TO 40 CONTINUE' C C** READ IN MEAN LAYER VALUE AS SPECIFIED IN LVLS C cc PRINT *, '** PRE-GETLVL **' cc PRINT *, 'FYMDH = ',FYMDH,' LVLS = ',LVLS,' IFHR = ',IFHR C CALL GETLVL ( mYMDH, LVLS, IFHR ) C CC IF ((UCOMP(1).NE.0.0).OR.(VCOMP(1).NE.0)) GO TO 50 CC IF (ISTART .EQ. 0) GO TO 60 CC GO TO 120 C 50 CONTINUE cc PRINT *, ' **** GOT TO 50 CONTINUE' C C** GET U,V COMPONENTS AT PROPER RADIUS FROM CENTER C ELONG = 360.0 - XLONC TLAT = XLATC + REDEG BLAT = XLATC - REDEG C SCL = COS(XLATC/57.3) EP = ELONG + REDEG/SCL WP = ELONG - REDEG/SCL C cc PRINT *, ' TLAT = ',TLAT,' BLAT = ',BLAT,' ELONG = ',ELONG cc PRINT *, ' XLATC = ',XLATC,' EP = ',EP,' WP = ',WP C RLAT(1) = TLAT RLON(1) = ELONG RLAT(2) = BLAT RLON(2) = ELONG RLAT(3) = XLATC RLON(3) = EP RLAT(4) = XLATC RLON(4) = WP C DO II = 1,22 KGDSO(II) = -1 enddo C IPOPT(1) = 0 CC IPOPT(1) = 1 CC IPOPT(2) = 32 CC IPOPT(2) = 25 IPOPT(2) = MXWAVE + 5 CC IPOPT(2) = 35 C DO II = 1,KF CROT(II) = 1.0 SROT(II) = 0.0 URAD(II) = -999 VRAD(II) = -999 enddo C CALL IPOLATEV (4,IPOPT,KGDSI,KGDSO,IPTS,KF,1,0,LGI,UCOMP,VCOMP, & KF,RLAT,RLON,CROT,SROT,IBO,LGO,URAD,VRAD,IRET) IF (IRET.NE.0) STOP ' *** ERROR *** INTERPOLATING' C DO II = 1,KF URAD(II) = URAD(II)*COS(RLAT(II)*3.141593/180.0) VRAD(II) = VRAD(II)*COS(RLAT(II)*3.141593/180.0) enddo C UC = (URAD(1) + URAD(2) + URAD(3) + URAD(4))*0.25 VC = (VRAD(1) + VRAD(2) + VRAD(3) + VRAD(4))*0.25 C cc WRITE (*,'( 2X,2F8.2 )') UC,VC cc WRITE (*,'(2(2X,2F8.2))') (URAD(K),VRAD(K), K = 1,KF) C UC = UC * 1.05 VC = VC * 1.05 C VNB = VC VEB = UC VB = SQRT(UC**2 + VC**2) ALPHA = ATAN2(-UC,VC) BETA = 2.0*7.292E-5*COS(XLATC/57.3)/(111137.0*57.3) C IF (IREFF.GE.1) GO TO 100 C IF (ISTART.EQ.1) GO TO 90 C C** FIND INITIAL MOVEMENT AT EACH LEVEL C CALL STHETA (THETAM,SPD) C THINIT(LVLS) = 360.0 - THETAM*57.3 THINIT(LVLS) = AMOD(THINIT(LVLS),360.0) SPDINI(LVLS) = SPD*1.94254 C 60 CONTINUE cc PRINT *, ' ***** GOT TO 60 CONTINUE' C LVLS = LVLS + 1 IF (LVLS .LE. 7) GO TO 20 C ISTART = 1 GAMMAH = 0.2 C C** PRINT OUT WINDS AT EACH LEVEL C cc PRINT *, ' I DO GET HERE' C DO LEV = 2,7 WRITE (*,3) MBLEV(LEV),THINIT(LEV),SPDINI(LEV) 3 FORMAT (/,' DIR & SPEED FOR LEVEL ',I3,'MB IS ',F6.1, & ' DEGREES ',F6.1,' KNOTS') enddo C IREFF = 0 LVLS = LVLF ILEV = LVLS*10 C DO 80 I = 1,6 ILEV = ILEV/10 IF (ILEV.EQ.0) GO TO 80 K = ILEV - (ILEV/10)*10 INLEVS(I) = MBLEV(K) 80 CONTINUE C WRITE (6,4) (INLEVS(I),I = 1,6) 4 FORMAT(/,' LEVELS TO BE USED ARE ',6I8) C GO TO 20 C 90 CONTINUE cc PRINT *, ' ****** GOT TO 90 CONTINUE' C C** CHECK THE OBSERVATION C VON = SPD * COS(THETAM) VOW = SPD * SIN(THETAM) IF (VS .GT. 0.0) THEN VON = AMAX1(VON,VC - ABS(VC*0.20)) VON = AMIN1(VON,VC + 1.5) ELSE VON = AMAX1(VON,VC + ABS(VC*0.20)) VON = AMIN1(VON,VC - 1.5) ENDIF C VOW = AMAX1(VOW,-UC - ABS(UC*0.15)) VOW = AMIN1(VOW,-UC + 2.0) SPD = SQRT(VOW**2 + VON**2) C THETAM = ATAN2(VOW,VON) DIR = THETAM*57.3 C CALL RADIUS (DIR) C THETAM = DIR/57.3 C REFF = AMAX1(REFF,60000.0) RX = AMAX1(R2,500000.0) REFF = AMIN1(REFF,RX) REDEG = AMAX1(REFF,R2,4.E5)/111137.0 C MXWAVE = 360.0/(4.0*AMAX1(REDEG,SFAC)) C cc PRINT *, ' ### REFF = ',REFF,' R2 = ',R2,' RX = ',RX cc PRINT *, ' #### REDEG = ',REDEG,' MXWAVE = ',MXWAVE C IREFF = IREFF + 1 C GO TO 50 C 100 CONTINUE cc PRINT *, ' ******* GOT TO 100 CONTINUE' C CALL STHETA (THETAM,SPEED) C XLONC = XLONC + SPEED*3600.0*SIN(THETAM)/(111137.0* & COS(XLATC/57.3)) XLATC = XLATC + SPEED*3600.0*COS(THETAM)/111137.0 XLATP = XLATC XLONP = XLONC IF (XLONP.GT.180.0) XLONP = XLONP - 360.0 C DO I = 1,IPTS UCOMP(I) = UCOMP(I) + UTEND(I) VCOMP(I) = VCOMP(I) + VTEND(I) enddo C IHR = IHR + 1 C 6 FORMAT (' FORECAST HOUR=',I8,'PREDICTED ANGLE, SPEED',2F8.2) C CC IF (MOD(IHR,12) .NE.0) GO TO 50 IF (MOD(IHR, 6) .NE.0) GO TO 50 KHR = IHR - INITHR C C** SAVE THE BAM FORECAST FOR ATCF FORMATING C IF (MOD(KHR,12).EQ.0) THEN IFP = IFP + 1 FPBAM(1,IFP) = XLATP FPBAM(2,IFP) = XLONP ENDIF C PRINT 7, IHR,KHR,XLATP,XLONP 7 FORMAT (/,' BAM FORECAST POSITION AT ',I3,' HRS AFTER MODEL OR ', & I3,' HRS AFTER INPUT POSITION',/,' LAT, LONG= ',F5.1,',',F6.1) C C** THE FORECAST FIELDS GO TO 78 HOURS C IF ( IHR .LT. IHREND ) THEN IF ( MOD( IHR, 6 ) .EQ. 0 ) THEN IFHR = IFHR + 6 IF ( INITHR .EQ. 6 ) GO TO 40 IF ( INITHR .EQ. 12 .AND. IFHR .LE. 120 ) GO TO 40 ENDIF GO TO 50 ENDIF C 120 CONTINUE C C** WRITE OUT THE ATCF FORECAST C WRITE (*,'(/,'' THE '',a4,'' FORECAST POSITIONS FOR '',a8, & '' ON '',a10,/)') bam_type, strmid, aymdh C WRITE (*,'('' 000 HR = '',F5.1,'' N '',F6.1,'' W'')') XLAT2,XLON2 C DO K = 1, ntau ibam( k + 1, 1 ) = int ( fpbam( 1, k )*10.0 + 0.5 ) ibam( k + 1, 2 ) = int ( fpbam( 2, k )*10.0 + 0.5 ) c EW = 'W' FPRITE = ABS(FPBAM(2,K)) C IF (FPBAM(2,K).LT.0.0) THEN EW = 'E' FPBAM(2,K) = 360.0 + FPBAM(2,K) ENDIF C WRITE (*,'(1X,A3,'' HR = '',F5.1,'' N '',F6.1,1X,A)') FTIME(K), & FPBAM(1,K),FPRITE,EW C enddo C C** OPEN, WRITE TO AND CLOSE THE ATCF OUTPUT FILE C OPEN ( 41, FILE='bamadv.dat', STATUS='UNKNOWN' ) C call newWriteAidRcd ( 41, strmid, aymdh, bam_type, ibam ) C CLOSE (41) C STOP '***** BAM HAS FINISHED *****' C C** ERROR MESSAGES C 1010 print *, ' ERROR - opening file = ', input_file, ' istat = ', istat STOP C 1020 print *, ' ERROR - reading file = ', input_file, ' istat = ', istat STOP C 1030 print *, ' ERROR - reading data = ', input_file, ' istat = ', istat STOP C 1040 print *, ' ERROR - inquiring file = ', flname, ' IOS = ', IOS STOP C END C***************************************************************** C DATA SET MULTIPHM AT LEVEL 042 AS OF 04/22/86 SUBROUTINE GETLVL ( FYMDH, LVLS, IFHR ) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: GETLVL GET LEVELS FOR BAM FORECAST C PRGMMR: DON MARKS ORG: W/NMC42 DATE: 87-06-25 C C ABSTRACT: READS IN U AND V-COMPONENT WINDS FOR LEVELS SPECIFIED C IN LEVS AND COMBINES THEM INTO A SINGLE WEIGHTED MEAN WIND FIELD C C PROGRAM HISTORY LOG: C 87-06-25 DON MARKS C 88-02-18 DON MARKS ADD CHECK FOR VERTICAL SHEAR C 91-05-22 JIM GROSS REMOVED VERTICAL SHEAR CHECK C 92-07-28 JIM GROSS - CLEAN UP ROUTINE C C USAGE: CALL GETLVL(LVLS,IFHR,LHR) C INPUT ARGUMENT LIST: C LVLS - ARRAY (6) LISTING LEVELS DESIRED C LHR - HOUR PAST FORECAST FILE AT WHICH TO START C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C UCOMP - ARRAYS CONTAINED IN COMMON-- MNWIND -- C VCOMP - CONTAINING U-COMPONENTS, V-COMPONENTS C UTEND - AND THEIR 1-HOUR TENDENCIES C VTEND - C C OUTPUT FILES: (DELETE IF NO OUTPUT FILES IN SUBPROGRAM) C FT06F001 - DIAGNOSTIC, ERROR PRINT C C REMARKS: C C ATTRIBUTES: C LANGUAGE: VS FORTRAN C MACHINE: NAS C C$$$ PARAMETER (IPTS=360*181) C COMMON /MNWIND/ UCOMP(IPTS),UTEND(IPTS),VCOMP(IPTS),VTEND(IPTS), & KGDSI(22),MXWAVE C REAL UPART(IPTS),VPART(IPTS) INTEGER JPDS(25),JGDS(22),KPDSI(25) INTEGER LEVEL(8),IUNITB(2),IUNITI(2) LOGICAL LGI(IPTS) C CHARACTER*2 CYCLE character*3 FHR00, FHRP6 CHARACTER*8 YRMNDY CHARACTER*10 FYMDH CHARACTER*80 FILE21, FILE22, FILE31, FILE32 C DATA IUNITB /21,22/ DATA IUNITI /31,32/ DATA LEVEL /100,200,300,400,500,700,850,1000/ DATA IROMB /1/ C C** SETUP THE FILE NAMES C YRMNDY = FYMDH(1: 8) CYCLE = FYMDH(9:10) C if ( ifhr .lt. 100 ) WRITE ( FHR00, '(I2.2)' ) IFHR if ( ifhr .ge. 100 ) WRITE ( FHR00, '(I3.3)' ) IFHR IFHR6 = IFHR + 6 if ( ifhr6 .lt. 100 ) WRITE ( FHRP6, '(I2.2)' ) IFHR6 if ( ifhr6 .ge. 100 ) WRITE ( FHRP6, '(I3.3)' ) IFHR6 C CARK FILE21 = '/com/gfs/prod/gfs.'//YRMNDY//'/gfs.t' CARK & //CYCLE//'z.pgrbf'//FHR00 CARK FILE31 = '/com/gfs/prod/gfs.'//YRMNDY//'/gfs.t' CARK & //CYCLE//'z.pgrbif'//FHR00 CARK FILE22 = '/com/gfs/prod/gfs.'//YRMNDY//'/gfs.t' CARK & //CYCLE//'z.pgrbf'//FHRP6 CARK FILE32 = '/com/gfs/prod/gfs.'//YRMNDY//'/gfs.t' CARK & //CYCLE//'z.pgrbif'//FHRP6 FILE21 = '/tmp/ARK/TEST-DATA/gfs.'//YRMNDY// & '/gfs.t'//CYCLE//'z.pgrbf'//FHR00 FILE31 = '/tmp/ARK/TEST-DATA/gfs.'//YRMNDY// & '/gfs.t'//CYCLE//'z.pgrbif'//FHR00 FILE22 = '/tmp/ARK/TEST-DATA/gfs.'//YRMNDY// & '/gfs.t'//CYCLE//'z.pgrbf'//FHRP6 FILE32 = '/tmp/ARK/TEST-DATA/gfs.'//YRMNDY// & '/gfs.t'//CYCLE//'z.pgrbif'//FHRP6 C PRINT *, FILE21 PRINT *, FILE31 PRINT *, FILE22 PRINT *, FILE32 C CALL BAOPENR (IUNITB(1),FILE21,IER) IF (IER.NE.0) GO TO 1010 CALL BAOPENR (IUNITI(1),FILE31,IER) IF (IER.NE.0) GO TO 1020 CALL BAOPENR (IUNITB(2),FILE22,IER) IF (IER.NE.0) GO TO 1030 CALL BAOPENR (IUNITI(2),FILE32,IER) IF (IER.NE.0) GO TO 1040 C C** ZERO THE ARRAYS AND INITIALIZE CONSTANTS C DO 10 I = 1,IPTS UCOMP(I) = 0.0 UTEND(I) = 0.0 VCOMP(I) = 0.0 VTEND(I) = 0.0 10 CONTINUE C WTOT1 = 0.0 WTOT2 = 0.0 C IRET = 0 JRET = 0 JREW = 0 C DO 15 II = 1,25 JPDS(II) = -1 15 CONTINUE JPDS(3) = 3 JPDS(6) = 100 C C** PROCESS THE DESIRED LEVELS C LVLS2 = LVLS 20 LEV = LVLS2 - (LVLS2/10)*10 LVLS2 = LVLS2/10 C IF (LEV.NE.0) THEN C WFACT = (LEVEL(LEV + 1) - LEVEL(LEV - 1))/100.0 IF (LEV.EQ.7) WFACT = WFACT/2.0 C cc PRINT *, ' *** IN-GETLVLS LEV = ',LEV,' LVLS2 = ',LVLS2 cc PRINT *, ' **** WFACT = ',WFACT,' WTOT1 = ',WTOT1 C JPDS(7) = LEVEL(LEV) C C** FIRST LOOP DOES THE CURRENT WIND COMPONENTS AND THE SECOND LOOP C** DOES THE 12 HOUR FORECAST WIND COMPONENTS C DO 40 LOOP = 1,2 C C** READ THE UPART FOR A LEVEL FOR U C JPDS(5) = 33 C PRINT *, LOOP,IUNITB(LOOP),IUNITI(LOOP),IPTS PRINT *, JPDS(3),JPDS(5),JPDS(6),JPDS(7),IRET C CALL GETGB (IUNITB(LOOP),IUNITI(LOOP),IPTS,0,JPDS, & JGDS,KF,JREW,KPDSI,KGDSI,LGI,UPART,IRET) IF (IRET.NE.0) GO TO 1050 print *, 'the number of points unpacked are ',KF print *, '1st u comp unpacked data element is ',UPART(1) C C** READ THE VPART FOR A LEVEL FOR V C JPDS(5) = 34 C CALL GETGB (IUNITB(LOOP),IUNITI(LOOP),IPTS,0,JPDS, & JGDS,KF,JREW,KPDSI,KGDSI,LGI,VPART,JRET) IF (JRET.NE.0) GO TO 1060 print *, 'the number of points unpacked are ',KF print *, '1st v comp unpacked data element is ',VPART(1) C C** SMOOTH THE U, V FIELDS TO WAVE 30 C IDRT = KGDSI(1) IMAX = KGDSI(2) JMAX = KGDSI(3) KM = 1 C CALL SPTRUNV (IROMB,MXWAVE,IDRT,IMAX,JMAX,IDRT,IMAX,JMAX, & KM,0,0,0,IPTS,0,0,IPTS,0, & UPART,VPART,.TRUE.,UPART,VPART, & .FALSE.,DUM,DUM,.FALSE.,DUM,DUM) C DO 30 I = 1,IPTS C IF (LOOP.EQ.1) THEN UCOMP(I) = (UCOMP(I)*WTOT1 + UPART(I)*WFACT)/(WTOT1 + WFACT) VCOMP(I) = (VCOMP(I)*WTOT1 + VPART(I)*WFACT)/(WTOT1 + WFACT) ELSE UTEND(I) = (UTEND(I)*WTOT1 + UPART(I)*WFACT)/(WTOT1 + WFACT) VTEND(I) = (VTEND(I)*WTOT1 + VPART(I)*WFACT)/(WTOT1 + WFACT) ENDIF C 30 CONTINUE 40 CONTINUE C WTOT1 = WTOT1 + WFACT C GO TO 20 C ENDIF C C** FIND THE VELOCITY TENDENCIES PER HOUR C DO 50 I = 1,IPTS UTEND(I) = (UTEND(I) - UCOMP(I))/6.0 VTEND(I) = (VTEND(I) - VCOMP(I))/6.0 50 CONTINUE C C** CLOSE THE FILE UNITS C CALL BACLOSE (IUNITB(1),IER) CALL BACLOSE (IUNITB(2),IER) CALL BACLOSE (IUNITI(1),IER) CALL BACLOSE (IUNITI(2),IER) C RETURN C C* ERROR MESSAGES C 1010 PRINT *, ' *** ERROR READING FIRST FILE, IER = ',IER STOP C 1020 PRINT *, ' *** ERROR READING FIRST FILE INDEX, IER = ',IER STOP C 1030 PRINT *, ' *** ERROR READING SECOND FILE, IER = ',IER STOP C 1040 PRINT *, ' *** ERROR READING SECOND FILE INDEX, IER = ',IER STOP C 1050 PRINT *, ' *** ERROR READING U GRIB, IRET = ',IRET,' LOOP = ',LOOP STOP C 1060 PRINT *, ' *** ERROR READING V GRIB, JRET = ',JRET,' LOOP = ',LOOP STOP C END C****************************************************************** C DATA SET PHMSOLVR AT LEVEL 010 AS OF 03/30/89 SUBROUTINE STHETA (TACT,SPEED) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: SOLVE FOR THETA (DIRECTION OF BAM FORECAST) C PRGMMR: DON MARKS ORG: W/NMC42 DATE: 88-04-29 C C ABSTRACT: SOLVES THE BETA AND ADVECTION MODEL EQUATIONS FOR C TROPICAL CYCLONE MOTION. GIVEN RADIUS AND INFLOW ANGLE C C PROGRAM HISTORY LOG: C 88-04-29 DON MARKS C 89-02-14 DON MARKS ADJUST FOR L-S CONVERGENCE AND C REDUCE NORTHWARD DRIFT C 92-07-16 JIM GROSS CLEANED UP THE ROUTINE AND TRIED TO C UNDERSTAND IT C C USAGE: CALL STHETA (TACT,SPEED) C C OUTPUT ARGUMENT LIST: (INCLUDING WORK ARRAYS) C DVDTX - CHANGE OF REL VORTICITY W.R.T. TIME IN DIR OF MOTION C TACT - ANGLE OF MOTION OF STORM C SPEED - SPEED OF MOTION, (IN M/S) C C C REMARKS: USES REFF, GAMMA FROM COMMON ARRAYS C THIS SUBROUTINE DETERMINES WHAT PHYSICAL PROCEDURES ARE C USED IN THE MODEL C C ATTRIBUTES: C LANGUAGE: VS FORTRAN C MACHINE: NAS C C$$$ C C FINDS THE VALUE OF THETA M, THE DIRECTION, AND C THE SPEED OF THE STORM MOTION C COMMON /OBSV/ BETA,GAMMA,REFF,VS,XI,VB,ALPHA,GAMMAH,C, & SPD,VNB,VEB,XLATC C DT = 0.01745 TACT = ALPHA VB1 = SQRT(VNB**2 + VEB**2) SPEED = VB1 DVDTX = 0.0 C C** **** **** **** DRIFT **** **** **** C GANGLE = ATAN(GAMMA) UDRFT = (REFF**2)*BETA C C** **** **** **** DRIFT **** **** C XLATR = ABS(XLATC/57.3) T1 = 0.0 BGS = BETA*GAMMAH/4.0 K1 = (VS**2)/(2.0*7.29E-5*6.37E6*SIN(XLATR)) C DO 20 I = 1,360 T1 = T1 + DT CT1 = COS(T1) ST1 = SIN(T1) ST = SIN(T1 + GANGLE) CT = COS(T1 + GANGLE) IF (VS.LE.0.0) THEN CT = -CT CT1 = -CT1 ENDIF C C** CHANGE COMPONENTS AND MULTIPLY BY 2/PI C VNS = AMIN1(ABS(VNB/VS),1.0) THDIV = ASIN(VNS) CFACT = COS(THDIV) C C** THIS LINE IS A DON MARKS TEST C CCDM WNORTH = UDRFT*0.4627*CT*CFACT + 0.06*COT(XLATR)*CT1 C C** THE OPERATIONAL ATCF BAM'S (BAMD, BAMM, BAMS) CONTAIN THE C** FOLLOWING TWO LINES TO INCLUDE THE BETA EFFECT AND C** CLIMATOLOGICAL UPWARD MOTION C WNORTH = UDRFT*0.4627*CT + 0.06*(1.0/tan(XLATR))*CT1 WWEST = UDRFT*0.6366*ST + 0.06*(1.0/tan(XLATR))*1.414*ST1 cc WNORTH = UDRFT*0.4627*CT + 0.06*COT(XLATR)*CT1 cc WWEST = UDRFT*0.6366*ST + 0.06*COT(XLATR)*1.414*ST1 C C** USE THESE TWO LINES TO INCLUDE THE BETA EFFECT ONLY (BAMB) C CC WNORTH = UDRFT*0.4627*CT CC WWEST = UDRFT*0.6366*ST C C** USE THESE TWO LINES TO INCLUDE NO BETA EFFECT (BAMA) C CC WNORTH = 0.0 CC WWEST = 0.0 C C** THE OPERATIONAL ATCF BAM'S (BAMD, BAMM, BAMS) CONTAIN THE C** FOLLOWING TWO LINES TO INCLUDE ? EFFECT (I DON'T KNOW) C** THE RESEARCH BAMS (BAMA, BAMB) DO NOT CONTAIN THESE LINES C WNORTH = WNORTH - (K1/3.1416)*CT1 WWEST = WWEST + (K1/2.0)*ST1 C DZSDR = VS*(1.0 - XI**2)/REFF**2 C DZDT = (VB1*COS(T1 - ALPHA) + WNORTH + WWEST)*DZSDR C C** THE OPERATIONAL ATCF BAM'S (BAMD, BAMM, BAMS) CONTAIN THE C** FOLLOWING LINE TO INCLUDE ? EFFECT (I DON'T KNOW) C** THE RESEARCH BAMS (BAMA, BAMB) DO NOT CONTAIN THIS LINE C DZDT = DZDT + BGS*(ST1 - CT1) C C** THIS LINE IS A DON MARKS TEST C CCDM DZDT = DZDT + BGS*(ST1 - CT1) - VB1*BETA*CT1*SIN(ALPHA - T1) C IF (VS.GT.0) GO TO 10 C C** SOUTHERN HEMISPHERE CASE C IF (DZDT.GT.DVDTX) GO TO 20 DVDTX = DZDT C C** INCLUDE FOR VORTEX ASYMMETRIES C C** THE OPERATIONAL ATCF BAM'S (BAMD, BAMM, BAMS) AND BAMB CONTAIN C** THE FOLLOWING LINE TO INCLUDE THE BETA EFFECT SPEED = DZDT/(DZSDR - BETA*ST*1.27) C C** USE THIS LINE TO INCLUDE NO BETA EFFECT (BAMA) C CC SPEED = DZDT/DZSDR C TF = T1 CF = CT1 GO TO 20 C C** NORTHERN HEMISPHERE CASE C 10 IF (DZDT.LT.DVDTX) GO TO 20 DVDTX = DZDT C C** INCLUDE FOR VORTEX ASYMMETRIES C C** THE OPERATIONAL ATCF BAM'S (BAMD, BAMM, BAMS) AND BAMB CONTAIN C** FOLLOWING LINE TO INCLUDE THE BETA EFFECT C SPEED = DZDT/(DZSDR + BETA*ST*1.27) C C** USE THIS LINE TO INCLUDE NO BETA EFFECT (BAMA) C CC SPEED = DZDT/DZSDR C TF = T1 CF = CT1 C 20 CONTINUE C TACT = TF VFN = SPEED*CF C C** INCLUDE FOR VORTEX ASYMMETRIES C C** THE OPERATIONAL ATCF BAM'S (BAMD, BAMM, BAMS) CONTAIN C** FOLLOWING LINE C AFACT = (BETA*REFF**2)/((1.0 - XI**2)*VS)*(VFN/VS) C C** THE RESEARCH ATCF BAM'S (BAMA, BAMB) CONTAIN C** FOLLOWING LINE C CC AFACT = 0.0 C VCN = VFN*(1.0/(1.0 + ABS(AFACT))) VCW = SPEED*SIN(TF) + AFACT*VCN C IF (VS.LT.0.0) VCN = -VCN C TACT = ATAN2(VCW,VCN) SPEED = SQRT(VCW**2 + VCN**2) C RETURN END C********************************************************************** SUBROUTINE RADIUS (DIR) C$$$ SUBPROGRAM DOCUMENTATION BLOCK C . . . . C SUBPROGRAM: RADIUS FINDER FOR BAM MODEL C PRGMMR: DON MARKS ORG: W/NMC42 DATE: 88-04-29 C C ABSTRACT: CALCULATES THE EFFECTIVE RADIUS FOR USE BY THE C BETA AND ADVECTION MODEL C C PROGRAM HISTORY LOG: C 88-04-29 DON MARKS C 92-05-14 JIM GROSS -- CLEANED UP THE CODE AND TRIED TO C UNDERSTAND IT C 97-06-25 JIM GROSS -- REMOVED XLAT2 AND XLON2 AS ARGUMENTS SINCE C THEY WERE NOT USED C C USAGE: CALL RADIUS (DIR) C INPUT ARGUMENT LIST: C DIR - DIRECTION OF CURRENT STORM MOTION (90 EQUALS TO WEST) C C C REMARKS: CHANGES REFF, AND GAMMA IN COMMON FOR USE LATER C C ATTRIBUTES: C LANGUAGE: VS FORTRAN C MACHINE: NAS C C$$$ C C** CALCULATES EFFECTIVE RADIUS FROM PAST MOVEMENT C COMMON /OBSV/ BETA,GAMMA,REFF,VS,XI,VB,ALPHA,GAMMAH,C, & SPD,VNB,VEB,XLATC C THETAM=DIR*0.01745 CC IDIAG = 0 C C** ESTIMATE ADVECTION USING BETA*REFF**2=1, XI=0.5 C REFF = SQRT(1.2/BETA) VX = ABS(C/SQRT(REFF)) C FACTW = (0.75*VX/(0.75*VX + 1.0)) IF (VEB.GT.0.0) FACTW = (0.75*VX/(0.75*VX - 1.0)) C FACTN = 1.0 - 2.0/(VX*3.1416) C X1 = SPD*SIN(THETAM)/FACTW + VEB X2 = SPD*COS(THETAM)/FACTN - VNB C GAMMA = 0.6285 - ATAN2(X1,X2) C C** TO CHANGE NORTH COMP -- USE 0.6285 C IF (C.LT.0.0) X2 = -X2 IF (C.LT.0.0) GAMMA = ATAN2(X1,X2) + 0.6285 C IF (GAMMA.GT. 3.1416) GAMMA = 6.283 - GAMMA IF (GAMMA.LT.-3.1416) GAMMA = 6.283 + GAMMA IF (GAMMA.GT. 1.57) GAMMA = 0.1 IF (GAMMA.LT.-1.57) GAMMA = -0.1 C X1 = AMAX1(X1,0.0) X2 = AMAX1(X2,0.0) X3 = SQRT(X1**2 + X2**2) C REFF = SQRT(X3/(BETA*0.63662)) C C** TO CHANGE NORTH COMP -- USE 0.61688 C REFF = AMAX1(REFF,110000.0) REFF = AMIN1(REFF,350000.0) C GAMMA = AMIN1(GAMMA, 0.5) GAMMA = AMAX1(GAMMA,-0.5) C IF (VS.GT.0.0) GAMMA = AMAX1(GAMMA,-0.3) IF (VS.LT.0.0) GAMMA = AMIN1(GAMMA, 0.3) C DELR = 50000.0 C DO 20 KR = 1,4 DELR = DELR/5.0 REFF = REFF - 5.0*DELR REFF = AMAX1(REFF,60000.0) RX = 1.0E10 RF = REFF C DO 10 I = 1,11 VS = C/(REFF**XI) C CALL STHETA (THETA1,SPD1) C CC IF (IDIAG.EQ.1) PRINT 1, REFF,THETA1,SPD1 CC 1 FORMAT (' RADIUS, THETA, SPEED',3F12.4) C TH1 = THETAM - THETA1 R1 = (SPD1*SIN(TH1))**2 + (SPD - SPD1*COS(TH1))**2 C IF (R1.LT.RX) THEN RX = R1 RF = REFF GF = GAMMA ENDIF C REFF = REFF + DELR C 10 CONTINUE C REFF = RF GF1 = GAMMA*0.75 GAMMA = GF1 C CALL STHETA (THETA1,SPD1) C TH1 = THETAM - THETA1 R1 = (SPD1*SIN(TH1))**2 + (SPD - SPD1*COS(TH1))**2 GF2 = GAMMA*1.25 GAMMA = GF2 C CALL STHETA (THETA1,SPD1) C TH1 = THETAM - THETA1 R2 = (SPD1*SIN(TH1))**2 + (SPD - SPD1*COS(TH1))**2 GAMMA = GF C IF (R1.LT.RX) GAMMA = GF1 IF ((R2.LT.RX).AND.(R2.LT.R1)) GAMMA = GF2 C GAMMA = AMIN1(GAMMA, 0.5) GAMMA = AMAX1(GAMMA,-0.5) C IF (VS.GT.0.0) GAMMA = AMAX1(GAMMA,-0.3) IF (VS.LT.0.0) GAMMA = AMIN1(GAMMA, 0.3) C VS = C/(REFF**XI) C 20 CONTINUE C RETURN END