program poe_initial_prblty ! This program creates arrays of demaria probability for single storm ! Modified 05/30/2013 by M. DeMaria to replace old read_edeck ! with new one that fixes the problem with the more general ! e-deck files. Also, the read_edeck subroutine was deleted from ! this file, and is compiled separately. implicit none integer maxtau, alltau, ntr, last_tau, igpcef parameter ( maxtau = 10, alltau = 120 ) character (len=1) poe_num character (len=5) ctr character (len=8) stormid character (len=10) cur_ymdh character (len=200) storm_fst, storm_edeck ! storms directory character (len=100) outfilename real gpcet(0:maxtau) common /stormstuff/ stormid, cur_ymdh ctr=' ' ! Read in the current yymmddhh and forecast file name and ! process number the command line call getarg ( 1, cur_ymdh ) call getarg ( 2, storm_fst ) ! call getarg ( 2, stormid ) ! poe_num is not required in ATCF call getarg ( 3, poe_num ) ! number of iterations read from command line (optional) ! call getarg ( 3, ctr ) call getarg ( 4, storm_edeck ) call getarg ( 5, stormid ) ! call getenv("ATCFSTRMS",storms) ! ind=index(storms," ")-1 ! storm_fst = storms(1:ind)//'/'//stormid//'.fst' ! storm_edeck = storms(1:ind)//'/e'//stormid//'.dat' ! if (ctr .ne. ' ') then ! read(ctr,*)ntr ! else ! number of iterations, always 1k for now ntr = 1000 ! endif if (cur_ymdh .ne. ' ') then print *, ' poe_initial_prblty' print *, ' current dtg = ', cur_ymdh print *, ' fst file = ', storm_fst print *, ' iterations= ', ntr else print *, ' USAGE: poe_initial_prblty yyyymmddhh stormid' stop endif ! Open the probability output file ! output_file = 'prblty_'//poe_num//'.dat' ! call getenv("ATCFSTRMS",storms) ! ind=index(storms," ")-1 ! l=index(storm_fst," ")-1 outfilename = 'prblty_'//poe_num//'.dat' print *, outfilename ! Read the official forecast data for the storm call read_fcst ( storm_fst, cur_ymdh, last_tau ) ! If no OFCL forecast for the current date/time then if ( last_tau .lt. 1 ) then write(0,'(a,x,a10,x,a,x,a1)') 'WARNING: There is no official '// & & 'forecast for this date/time ', cur_ymdh, trim(storm_fst), poe_num stop endif ! Read the GPCE radii, if they exist call read_edeck( storm_edeck, cur_ymdh, last_tau, igpcef, gpcet) if (igpcef .eq. 0) then write(0,'(a,a)') 'WARNING: GPCE data not found in ',trim(storm_edeck) end if ! pass data into main processing code call dm_probability ( stormid, cur_ymdh, ntr, gpcet, outfilename) stop ' ***** The wind probability program has finished! *****' end !********1*********2*********3*********4*********5*********6*********7** SUBROUTINE read_fcst ( storm_fst, cur_ymdh, last_tau ) implicit none ! THIS SUBROUTINE reads the official forecast data INCLUDE 'dataformats.inc' integer maxtau, alltau parameter ( maxtau = 10, alltau = 120 ) integer tau real flat, flon, vmax,radius, alat180, alon360 common /fcst_data/ tau(0:maxtau), flat(0:maxtau), flon(0:maxtau), & & vmax(0:maxtau), radius(0:maxtau,4,4), basin integer istat, ios, itau, numcy, last_tau, num, m character (len=2) basin, cynum character (len=8) strmid character (len=10) fst_ymdh, cur_ymdh character (len=100) storm_fst type (AID_DATA) fstRcd, tauData !Zero the arrays tau( : ) = 0 flat( : ) = 0.0 flon( : ) = 0.0 vmax( : ) = 0.0 radius( :, :, : ) = 0.0 !Open, read all the official forecast records and close the ! strmid.fst file. !Determine the storms directory path open ( 21, file=storm_fst, status='old', iostat=ios, err=1010 ) call readARecord ( 21, fstRcd, istat ) if ( istat .le. 0 ) then last_tau = -1 close ( 21 ) return endif !Get the official forecast data by searching for all possible TAUs do itau = 0, alltau, 12 call getSingleTAU ( fstRcd, itau, tauData, istat ) if ( istat .ne. 1) cycle !Read the current year, month, day and hour if ( itau .eq. 0 ) then basin = tauData%atcfRcd(1)%basin cynum = tauData%atcfRcd(1)%cyNum numcy = tauData%aRecord(1)%cyNum strmid = basin//cynum//cur_ymdh(1:4) !Tropical cyclone with storm ID's greater then 89 are not processed !Changed from 79........ sampson if ( numcy .gt. 89 ) then last_tau = -1 print *, 'Storm number 80 or greater, strmid = ',strmid print *, 'Forecast file not processed!' return endif fst_ymdh = tauData%aRecord(1)%DTG print *, 'Storm ID = ', strmid print *, 'Forecast ymdh = ', fst_ymdh print *, 'Current ymdh = ', cur_ymdh ! If the current date/time is not the initial forecast time, do not process if ( fst_ymdh .ne. cur_ymdh ) then last_tau = -1 print *, 'Forecast ymdh not equal the current ymdh!' print *, 'Forecast file not processed!' return endif endif num = (tauData%aRecord(1)%tau)/12 tau(num) = tauData%aRecord(1)%tau ! Read the position and maximum wind flat(num)=alat180(tauData%aRecord(1)%lat,tauData%aRecord(1)%NS) flon(num)=alon360(tauData%aRecord(1)%lon,tauData%aRecord(1)%EW) vmax(num)= float( tauData%aRecord(1)%vmax ) ! Read the 34, 50 and 64-knot wind radii do m = 1, tauData%numrcrds if ( tauData%aRecord(m)%windcode .eq. 'AAA' ) then radius( num, m, 1 ) = float( tauData%aRecord(m)%radii(1)) radius( num, m, 2 ) = float( tauData%aRecord(m)%radii(1)) radius( num, m, 3 ) = float( tauData%aRecord(m)%radii(1)) radius( num, m, 4 ) = float( tauData%aRecord(m)%radii(1)) elseif ( tauData%aRecord(m)%windcode .eq. 'NEQ' ) then radius( num, m, 1 ) = float( tauData%aRecord(m)%radii(1)) radius( num, m, 2 ) = float( tauData%aRecord(m)%radii(2)) radius( num, m, 3 ) = float( tauData%aRecord(m)%radii(3)) radius( num, m, 4 ) = float( tauData%aRecord(m)%radii(4)) endif enddo end do last_tau = num !Print out the values read in print *, ' ' print *, ' last tau = ', last_tau print *, ' tau(last tau) = ', tau(last_tau) print *, ' ' do num = 0, last_tau write ( *, '( '' i = '', i3, '' tau = '', i3 ) ') num, tau(num) write ( *, '( '' lat = '', f5.1, '' lon = '', f5.1, '' vmax = '',f4.0 ) ') & & flat(num), flon(num), vmax(num) write ( *, '('' 34kt radii = '', 4f6.0 )' ) radius(num,1,1),& & radius(num,1,2), radius(num,1,3), radius(num,1,4) write ( *, '('' 50kt radii = '', 4f6.0 )' ) radius(num,2,1),& & radius(num,2,2), radius(num,2,3), radius(num,2,4) write ( *, '('' 64kt radii = '', 4f6.0 )' ) radius(num,3,1),& & radius(num,3,2), radius(num,3,3), radius(num,3,4) enddo RETURN 1010 print *, ' ERROR - opening file = ', storm_fst, ' ios = ', ios last_tau = -1 return END !********1*********2*********3*********4*********5*********6*********7** subroutine dm_probability ( stormid, cur_ymdh, ntr, gpcet, & & outfilename) use windspeed_model, only : mcprob implicit none ! This program computes the Demaria probabilities !Parameter statement for 1.0 degree resolution for a total area of: ! latitude: 1.0N to 60.0N longitude: 15.0W to 110.0W ! parameter (imin=15,imax=110,ioffset=-14,iarray=96,deltax=1.0) ! parameter (jmin= 1,jmax= 60,joffset= 0,jarray=60,deltay=1.0) !Parameter statement for 0.5 degree resolution for a total area of: ! latitude: 1.0N to 60.0N longitude: 15.0W to 110.0W ! parameter (imin=30,imax=220,ioffset=-29,iarray=191,deltax=0.5) ! parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5) ! parameter (ix=22,jx=12) !Parameter statement for 0.5 degree resolution for a total area of: ! latitude: 1.0N to 60.0N longitude: 15.0W to 100.0E ! parameter (imin=30,imax=520,ioffset=-29,iarray=491,deltax=0.5) ! parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5) ! parameter (ix=52,jx=12) !Parameter statement for 0.5 degree resolution for a total area of: ! latitude: 1.0N to 60.0N longitude: 1.0W to 100.0E integer imin, imax, ioffset, iarray, jmin, jmax, joffset, jarray integer ix, jx real deltax, deltay parameter (imin=2,imax=520,ioffset=-1,iarray=519,deltax=0.5) parameter (jmin= 2,jmax=120,joffset= -1,jarray=119,deltay=0.5) parameter (ix=52,jx=12) !Parameter statement for 0.1 degree resolution for a total area of: ! latitude: 1.0N to 60.0N longitude: 15.0W to 110.0W ! parameter (imin=150,imax=1100,ioffset=-149,iarray=951,deltax=0.1) ! parameter (jmin= 10,jmax= 600,joffset= -9,jarray=591,deltay=0.1) integer maxtau, numits, totpts integer igpcef, k, its, ite, ik, ij, isnum, ntr, ierr parameter ( maxtau=10, numits=20, totpts=iarray*jarray ) ! parameter ( maxtau=10, numits=20, totpts=iarray*jarray ) integer, dimension(0:maxtau) :: tau real, dimension(0:maxtau) :: flat, flon,vmax real radius(0:maxtau,4,4) character (len=2) basin common /fcst_data/ tau, flat, flon, vmax, radius, basin integer kstime_short(0:1), kstime_long(0:numits) logical diagnostic real plat(totpts), plon(totpts) real wprob(totpts, 0:numits), wprobs(totpts, 0:numits) real wprobi(totpts, 0:numits), wtprob(7,9) integer ktime(0:maxtau), ivmaxth(3) integer isave2(0:imax,0:jmax) real, dimension(0:maxtau) :: flatt, flont integer iyear, imonth, iday, ihour, ithresh, i, j, ntsave, iprall integer indxg(maxtau) real gpcet(0:maxtau) ! character (len=3) shortcut character (len=8) stormid character (len=10) cur_ymdh ! character (len=50) file_name character (len=50) outfilename data ivmaxth / 34, 50, 64 / data diagnostic / .false. / ! Flag to set whether or not to use GPCE radii in wind prob computation igpcef = 1 !igpcef = 0 !Set the TAUs and the initial Monte Carlo parameters do k = 0, maxtau ktime(k) = k*12 enddo do k = 0, numits ! kstime_long(k) = k*12 kstime_long(k) = k*6 enddo its = 0 ite = 120 ! Put longitude in a 0 to 360 coordinated system by basin do ik = 0, maxtau flatt(ik) = flat(ik) flont(ik) = 360.0 - flon(ik) if ( flont(ik) .eq. 360.0 ) flont(ik) = 0.0 enddo read (stormid, '(2x,i2)' ) isnum read (cur_ymdh, '(i4,3i2)' ) iyear, imonth, iday, ihour write(6,*) 'writing output to:', outfilename open ( 31, file=trim(outfilename), status='unknown',& & form='unformatted') do ithresh = 1, 3 ! if ( ithresh .eq. 1 ) then ! ! first pass through, use coarse grid to define max domain ! ! Calculate at probs at 5deg resolution to then mask higher ! ! resolution points on subsequent iterations ! ! shortcut = 'yes' ! ! if ( shortcut .eq. 'yes' ) then ! ! ** Specify lat/lon points (Input to mcprob is 1-D string ! ! ** of lat/lon point, with deg west negative). ! ij = 0 ! do j = jmax, 0, int(-5.0/deltay) ! do i = imax, 0, int(-5.0/deltax) ! ij = ij + 1 ! plat(ij) = deltay*float( j ) ! plon(ij) = 360.0 - deltax*float( i ) ! enddo ! enddo ! print *, ' totpts, ij = ', totpts, ij ! ntsave = 1 ! kstime_short(0) = its ! kstime_short(1) = ite ! ! iprall is option to print ensemble members. This is done in the text product. ! iprall = 1 ! call mcprob ( basin, flatt, flont, vmax,& ! & igpcef, gpcet,& ! & indxg, 'NONAME ', isnum, iyear, imonth, iday, ihour,& ! & maxtau, ktime, ntr,& ! & ithresh, its, ite, kstime_short, ntsave,& ! & iprall, radius(:,1,:), radius(:,2,:), radius(:,3,:), & ! & radius(:,4,:),& ! & ij, totpts, plon, plat,& ! & wprobs, wprobi, wprob, wtprob, ierr ) ! if ( ierr .ne. 0 ) then ! print *, ' ierr = ', ierr ! stop ' *** ERROR1 in mcprblty*** ' ! endif ! call define_wsp_mask(wprob(:, ntsave),plon,plat,ix,jx,imax,& ! & jmax,deltax, deltay,ivmaxth(ithresh),diagnostic,isave2) ! endif ! ithresh.eq.1 for isave2 definition ! if ( shortcut .eq. 'yes' ) stop ' Stop after shortcut!' ij = 0 ierr = 0 isave2 = 1 ! temporary fix to mask issues. Utilise all points ! Specify lat/lon points (Input to mcprob is 1-D string ! of lat/lon point, with deg west negative). do j = jmax, jmin, -1 do i = imax, imin, -1 if ( isave2( i, j ) .ne. 0 ) then ij = ij + 1 plat(ij) = deltay*float( j ) plon(ij) = 360.0 - deltax*float( i ) endif enddo enddo write(6,*) ' totpts, ij = ', totpts, ij ! iprall is option to print ensemble members. ! This is done in the text product. iprall = 0 call mcprob ( basin, flatt, flont, vmax, igpcef, gpcet,& & indxg, 'NONAME ', isnum, iyear, imonth, iday, ihour,& & maxtau,& & ktime, ntr, ithresh, its, ite, kstime_long, numits,& & iprall, radius(:,1,:), radius(:,2,:), radius(:,3,:), & & radius(:,4,:),& & ij, totpts, plon, plat,& & wprobs, wprobi, wprob, wtprob, ierr ) if ( ierr .ne. 0 ) then print *, ' ierr = ', ierr stop ' *** ERROR2 in mcprblty*** ' endif !Probabilities write(6,*) 'write out array to file' call write_probability_array(31, wprob, wprobi,& & numits,totpts, diagnostic, plon, plat, ivmaxth(ithresh),& & isave2, imax, imin, jmax, jmin) end do close(31) return end subroutine !********1*********2*********3*********4*********5*********6*********7** subroutine write_probability_array(filehandle, wprob, wprobi, & & ntsave,totpts, diagnostic, plon, plat, ivmaxth,isave2,& & imax, imin, jmax, jmin) implicit none integer, intent(in) :: filehandle, isave2(0:imax, 0:jmax), ivmaxth integer, intent(in) :: ntsave, totpts, imax, jmax, imin, jmin real, intent(in), dimension(totpts,0:ntsave) :: wprob, wprobi real, intent(in), dimension(totpts) :: plon, plat logical, intent(in) :: diagnostic !local: character (len=100) :: diag_fmt integer :: k, ij, j, i, jarray, iarray, ioffset, joffset parameter (jarray=119, iarray=519, ioffset=-1,joffset= -1) real, dimension( iarray, jarray) :: parray, parrayi parray = 0 parrayi = 0 diag_fmt = '(f8.3, 1x, f7.3, 1x, f8.3, 1x, i5 )' do k = 0, ntsave ij = 0 do j = jmax, jmin, -1 do i = imax, imin, -1 if ( isave2( i, j ) .ne. 0 ) then ij = ij + 1 parray(i+ioffset,j+joffset) = wprob(ij,k) parrayi(i+ioffset,j+joffset) = wprobi(ij,k) if ( diagnostic ) then if ( wprob( ij, k ) .gt. 0.1 ) then write( *, diag_fmt ) plon(ij), plat(ij), & & wprob(ij, k), ivmaxth endif endif endif write (filehandle) parray( i+ioffset, j+joffset ),& & parrayi( i+ioffset, j+joffset) enddo ! end i loop enddo ! end j loop enddo ! end k loop end subroutine write_probability_array !********1*********2*********3*********4*********5*********6*********7** subroutine write_diagnostic_data(file_name, array, ix, jx, fmt) integer ix, jx, i, j integer array(0:ix, 0:jx) character (len=*) fmt character (len=16) file_name write(6,*) 'dumping diagnostic to:', file_name open ( 35, file=trim(file_name), status='unknown' ) DO j = jx, 0, -1 WRITE ( 35, fmt ) ( array(i, j), i = ix, 0, -1 ) enddo CLOSE ( 35 ) end subroutine write_diagnostic_data !********1*********2*********3*********4*********5*********6*********7** real function alat180 ( alat90, ns ) implicit none ! Transforms given latitude from 0 to 90 degrees to -90 to 90 degrees real :: alat90 character (len=1) ns alat180 = alat90 if ( ns .eq. "S" ) alat180 = - alat90 return end function alat180 !********1*********2*********3*********4*********5*********6*********7** real function alon360 ( alon180, ew ) implicit none ! Transforms given longitude from 0 to 180 degrees to 0 to 360 degrees ! ( East > 180 ) real :: alon180 character (len=1) ew alon360 = alon180 if ( ew .eq. "E" ) alon360 = 360.0 - alon180 return end function alon360 !********1*********2*********3*********4*********5*********6*********7** subroutine define_irregular_bitmap(wprob, plon, plat, ivmaxth, ix,& & jx, isave1) implicit none ! input: integer, intent(in) :: jx, ix, ivmaxth real, intent(in), dimension(ix*jx) :: wprob, plon, plat integer, intent(out) :: isave1(ix, jx) ! local integer :: j, i, ij character (len=50) :: bmpfmt isave1 = 0 do j = jx, 0, -1 do i = ix, 0, -1 ij = ij + 1 ! write( *, 10 ) ij, plon( ij ), plat( ij ), ! & wprob( ij, 1 ), ivmaxth( ithresh ) if ( wprob( ij) .gt. 0.01 ) then bmpfmt = '( i5, f8.3, 1x, f7.3, 1x, f6.1, 1x, i5 )' write( *, bmpfmt ) ij, plon( ij ), plat( ij ), & & wprob( ij), ivmaxth if ( j.ne.jx ) isave1( i , j+1 ) = 1 if ( i.ne.ix .and. j.ne.jx ) isave1( i+1, j+1 ) = 1 isave1( i , j ) = 1 if ( i.ne.ix ) isave1( i+1, j ) = 1 endif enddo enddo end subroutine define_irregular_bitmap !********1*********2*********3*********4*********5*********6*********7** subroutine define_square_bitmap(wprob, plon, plat, ix, jx, isave1) implicit none ! input: integer, intent(in) :: jx, ix real, dimension(ix*jx), intent(in) :: wprob, plon, plat ! output: integer, intent(out):: isave1(0:ix, 0:jx) ! local real :: plonmin, plonmax, platmin, platmax integer :: j, i, ij plonmin = 180.0 plonmax = -180.0 platmin = 90.0 platmax = -90.0 ij = 0 do j = jx, 0, -1 do i = ix, 0, -1 ij = ij + 1 if ( wprob( ij) .gt. 0.01 ) then if (plon(ij) .lt. plonmin) plonmin=plon(ij) if (plon(ij) .gt. plonmax) plonmax=plon(ij) if (plat(ij) .lt. platmin) platmin=plat(ij) if (plat(ij) .gt. platmax) platmax=plat(ij) endif enddo enddo ij = 0 do j = jx, 0, -1 do i = ix, 0, -1 ij = ij +1 if ( plon(ij) .ge. plonmin .and. & & plon(ij) .le. plonmax .and.& & plat(ij) .ge. platmin .and.& & plat(ij) .le. platmax) then if ( j.ne.jx ) isave1( i , j+1 ) = 1 if ( i.ne.ix .and. j.ne.jx ) isave1( i+1, j+1 ) = 1 isave1( i , j ) = 1 if ( i.ne.ix ) isave1( i+1, j ) = 1 endif enddo enddo end subroutine define_square_bitmap !********1*********2*********3*********4*********5*********6*********7** subroutine upscale_bitmap(isave1,jmax, imax, ix, jx, deltax, & & deltay, diagnostic, isave2) implicit none ! input: integer, intent(in) :: jmax, imax, ix, jx, isave1(0:ix, 0:jx) real, intent(in) :: deltax, deltay logical, intent(in) :: diagnostic ! output: integer, intent(out) :: isave2(0:imax, 0:jmax) ! local: integer :: j, i, itest, jabv, jblo, ilft, irgt do j = jmax, 0, -1 jabv = min0( jx, int( ( float(j)*deltay + 4.0 )/5.0 )) jblo = max0( 0, int( ( float(j)*deltay - 1.0 )/5.0 )) do i = imax, 0, -1 ilft = min0( ix, int( ( float(i)*deltax + 4.0 )/5.0 )) irgt = max0( 0, int( ( float(i)*deltax - 1.0 )/5.0 )) ! Test to see if any of the four corner from the isave array bit ! are greater than 1 itest = isave1( ilft, jabv ) + isave1( ilft, jblo ) + & & isave1( irgt, jblo ) + isave1( irgt, jabv ) if ( itest .ne. 0 ) then if ( diagnostic )then write( *,'( 3i3,2x,3i3,4x,5i3 )')& & j, jabv, jblo, i, ilft, irgt,& & isave1( ilft, jabv ), isave1( ilft, jblo ),& & isave1( irgt, jblo ), isave1( irgt, jabv ),& & itest end if isave2( i, j ) = 1 else isave2( i, j ) = 0 endif enddo enddo end subroutine upscale_bitmap !********1*********2*********3*********4*********5*********6*********7** subroutine define_wsp_mask(wprob, plon, plat, ix, jx, imax, jmax,& & deltax, deltay, ivmaxth, diagnostic, isave2) implicit none ! input: integer, intent(in) :: jmax, imax, ix, jx, ivmaxth real, dimension(ix*jx), intent(in) :: wprob, plon, plat real, intent(in) :: deltax, deltay logical, intent(in) :: diagnostic ! output: integer, intent(out) :: isave2(0:imax, 0:jmax) ! local logical :: square_output integer :: isave1(0:ix, 0:jx) ! Use the calculated 5 degree probabilities to set the bit map isave1 = 0 square_output = .true. if (.not. square_output) then ! Jim's code to assign points. Allows irregular shape call define_irregular_bitmap(wprob, plon, plat,& & ivmaxth, ix, jx, isave1) else ! Replacement constructs rectangular grid for internal optimization call define_square_bitmap(wprob, plon, plat, ix, jx, isave1) endif ! Open the bitmap file and write it out. if ( diagnostic ) then call write_diagnostic_data('prblty_tbl_5.dat', isave1, & & ix, jx, '(53i2)') endif ! Create the large bit map call upscale_bitmap(isave1,jmax, imax, ix, jx,deltax,& & deltay, diagnostic,isave2) if ( diagnostic ) then ! Open the large bitmap file and write it out. call write_diagnostic_data('prblty_tbl_1.dat', isave2,& & imax, jmax, '(111i2)') endif end subroutine define_wsp_mask