module windspeed_routines implicit NONE real, parameter, private :: pi = 3.14159265 real, parameter, private :: e = 2.71828183 real, parameter, private :: rtd = 57.2957802 real, parameter, private :: dtk = 111.1 real, parameter, private :: dtr = 0.01745329 contains subroutine ctorh(x,y,r,theta) ! This routine converts from Cartesian coordinates ! to radial coordinates, where theta is in ! degrees measured clockwise from ! the +y-axis (standard meteorological heading). ! Input / Output Variables real, intent(in):: x, y real, intent(out) :: r, theta !local variables r = sqrt(x*x + y*y) if (r .le. 0.0) then theta = 0.0 return endif theta = rtd*acos(x/r) if (y .lt. 0.0) theta = 360.0 - theta ! Convert theta to heading theta = 90.0 - theta if (theta .lt. 0.0) theta = theta + 360.0 return end subroutine ctorh subroutine rhtoc(r,thetah,x,y) ! This routine converts from radial coordinates ! to Cartesian coordinates, where theta is in ! degrees measured lockwise from ! the +y-axis (standard meteorological heading). ! Convert theta from heading to standard definition ! input / output variables real, intent(in) :: r, thetah real, intent(out) :: x, y ! local variable real theta theta = 90.0 - thetah if (theta .lt. 0.0) theta = theta + 360.0 x = r*cos(theta/rtd) y = r*sin(theta/rtd) return end subroutine rhtoc subroutine distk(rlon1,rlat1,rlon2,rlat2,dx,dy,rad) ! This routine calculates the distance in km (rad) between the ! points (rlon1,rlat1) and (rlon2,rlat2) using an approximate ! formula. The lon and lat are in deg E and N. The east and ! north components of the distance (dx,dy) are also calculated. real, intent(in) :: rlat1, rlon1, rlon2, rlat2 real, intent(out) :: dx, dy, rad ! local variables real :: cfac cfac = cos(0.5*dtr*(rlat1+rlat2)) dx = dtk*(rlon2-rlon1)*cfac dy = dtk*(rlat2-rlat1) rad = sqrt(dx*dx + dy*dy) return end subroutine distk subroutine distki(rlon1,rlat1,dx,dy,rlon2,rlat2) ! This routine performs the inverse of routine distk. ! The latitude and longitude (rlon2,rlat2) of the point a distance ! dx and dy (km) from the point (rlon1,rlat1) is calculated ! using an approximate formula. The lon and lat are in deg E ! and deg N. real, intent(in) :: rlat1, rlon1, dx, dy real, intent(out) :: rlon2, rlat2 ! local variables real :: cfac if (rlat1 .eq. 90.0) then rlat2 = rlat1 rlon2 = rlon1 return endif rlat2 = rlat1 + dy/dtk cfac = cos(0.5*dtr*(rlat1+rlat2)) rlon2 = rlon1 + dx/(dtk*cfac) return end subroutine distki subroutine rcal(rlat1,rlon1,rlat2,rlon2,rnm,azdeg) ! ********************************************************************* ! This routine calculates the radial distance in ! nm from the point (rlat1,rlon1) to (rlat2,rlon2) ! using an approximate formula. The azimuthal angle ! (azdeg) in degrees is also calculated, assuming ! the point (rlat1,rlon1) is at the origin, and the ! x-axis is east and y-axis is north. The angle azdeg ! is measured counter-clockwise in degrees from ! east (i.e., azdeg=0 if the point rlat2,rlon2 is ! due east of the point rlat1,rlon2, azdeg=90 if ! rlat2,rlon2 is due north of rlat1,rlon2, etc. ! The variables rlat,rlon are input in degrees N and ! E, respectively (degrees W is negative). ! ********************************************************************* real, intent(in) :: rlat1, rlon1, rlat2, rlon2 real, intent(out) :: rnm, azdeg real :: cfac, dx, dy cfac = cos(0.5*dtr*(rlat1+rlat2)) dx = cfac*60.0*(rlon2-rlon1) dy = 60.0*(rlat2-rlat1) rnm = sqrt(dx*dx + dy*dy) if (rnm .le. 0.0) then azdeg = 0.0 elseif (dy/rnm .gt. 1.0) then ! round off fix azdeg = 90.0 elseif (dy/rnm .lt. -1.0) then ! round off fix azdeg = 270.0 else azdeg = (asin(dy/rnm))/dtr if ( dx .lt. 0.0) azdeg = 180.0 - azdeg if (azdeg .lt. 0.0) azdeg = azdeg + 360.0 endif return end subroutine rcal subroutine shcal(flon,flat,ftime,fspd,fhead,mt,nt) !This routine calculates the speed (kt) and heading along !forecast storm track. It is assumed that routine linte !has already been called to fill in missing time periods !when possible. The input lat/lon are in deg N and deg E, !(deg W negative) respectively. implicit none real, dimension(0:mt), intent(in) :: flon,flat,ftime real, dimension(0:mt), intent(out) :: fspd,fhead integer, intent(in) :: nt, mt real :: dt, rlat1, rlat2, rlon1, rlon2, head, rad, dx, dy integer :: k do k=1,nt rlat2 = flat(k ) rlat1 = flat(k-1) rlon2 = flon(k ) rlon1 = flon(k-1) dt = ftime(k)-ftime(k-1) if(dt .gt. 0.0)then call distk(rlon1,rlat1,rlon2,rlat2,dx,dy,rad) fspd(k) = rad*60.0/(111.1*dt) call ctorh(dx,dy,rad,head) fhead(k) = head else fhead(k) = 0.0 fspd(k) = 0.0 endif end do fspd(0) = fspd(1) fhead(0) = fhead(1) return end subroutine shcal subroutine yint(x,f,n,iflag,lflag,xi,fi,ierr) !This routine applies a quadratic interpolation procedure !to f(x) between x(1) and x(n). f(x) is assumed to be !represented by quadratic polynomials between the points !x(i). The polynomials are chosen so that they equal f(i) !at the points x(i), the first derviatives on either !side of the interior x(i) match at x(i), and the second !derivative of the approximated function integrated !over the domain is minimized. !This verision is for interpolating latitude !Input: x(1),x(2) ... x(n) The x values (must be sequential) !f(1),f(2) ... f(n) The function values !n The number of x,f pairs !iflag Flag for initialization !=1 for coefficient calculation !=0 to use previous coefficients !lflag Flag for linear interpolation !=0 to perform linear interpolation !=1 to perform quadratic interpolation !xi The x value at which to interpolate f !Output: fi The interpolated function value !ierr Error flag !=0 Normal return !=1 Parameter nmax is too small !=2 The x values are not sequential !=3 Coefficient iteration did not !converge !=4 Mix-up finding coefficients !=5 if xi .gt. x(n) or .lt. x(1), !xi is set to nearest endpoint !before the interpolation !Note: fi is set to -99.9 if !ierr=1,2,3 or 4 integer, parameter :: nmax=100 real, dimension(n), intent(in) :: x,f integer, intent(in) :: n, iflag, lflag real, intent(inout) :: xi real, intent(out) :: fi integer, intent(out) :: ierr !Save variables real, dimension(nmax) :: ay, by, cy !Temporary local variables integer :: ii, i logical :: interval_found common /ysave/ ay,by,cy !Initialize error flag ierr = 0 fi = -99.9 interval_found = .false. !Check to make sure nmax is large enough, and n is .gt. 1 if (n .gt. nmax .or. n .lt. 2) then ierr=1 return endif if (iflag .eq. 1) then ! Initialize the function weights call init_int(x,f,n,lflag,fi,ierr, ay, by, cy) if (ierr .ne. 0) return endif !Interpolate the function !Check for xi out of bounds if (xi .lt. x(1)) then xi = x(1) ierr = 5 endif if (xi .gt. x(n)) then xi = x(n) ierr = 5 endif !Find the interval for the interpolation ii = 1 do i=2,n if (xi .le. x(i)) then ii = i-1 interval_found = .true. exit endif end do if (.not. interval_found) then ierr = 4 return endif fi = ay(ii) + by(ii)*xi + cy(ii)*xi*xi return end subroutine yint subroutine init_int(x,f,n,lflag,fi,ierr, ay, by, cy) integer, parameter :: nmax=100 real, dimension(n), intent(in) :: x,f integer, intent(in) :: n, lflag real, intent(out) :: fi integer, intent(out) :: ierr real, dimension(nmax), intent(out) :: ay, by, cy !Temporary local variables real, dimension(nmax) :: df, dx, gm, ct real :: thresh, slt, slo, rel, eps, dev, den real :: dsdc1, d, cft, cfsave, cf0, bb, aa integer :: nit, lutest, k, j, idbug, i logical :: converged !Perform the initialization for later interpolation !Specify unit number for debug write statements !and debug flag idbug = 0 lutest = 6 converged=.false. !Initialize error flag ierr = 0 !Specify minimum reduction in cost function for convergence thresh = 1.0e-10 !Check to make sure x is sequential do i=1,n-1 if (x(i) .ge. x(i+1)) then ierr=2 fi = -99.9 return endif end do !Check for special case where n=2. Only linear interpolation !is possible. if (n .eq. 2) then cy(1) = 0.0 by(1) = (f(2)-f(1))/(x(2)-x(1)) ay(1) = f(1) - by(1)*x(1) return endif !Calculate x and f differences do i=1,n-1 df(i) = f(i+1)-f(i) dx(i) = x(i+1)-x(i) end do !Calculate domain size d = x(n) - x(1) !Check for linearity of input points eps = 1.0e-10 bb = (f(2)-f(1))/(x(2)-x(1)) aa = f(1) - bb*x(1) dev = 0.0 do i=3,n dev = dev + abs(aa + bb*x(i) - f(i)) end do if (dev .lt. eps .or. lflag .eq. 0) then cy(:) = 0.0 converged = .true. endif !Iterate to find the c-coefficients cy(1) = 0.0 nit = 100 slt = 0.01 cfsave = 1.0e+10 do k=1,nit if (converged) exit !Calculate c values do i=2,n-1 cy(i) = -cy(i-1)*dx(i-1)/dx(i) & & -df(i-1)/(dx(i)*dx(i-1)) & & +df(i )/(dx(i)*dx(i )) end do !Calculate current value of cost function cf0 = 0.0 do i=1,n-1 cf0 = cf0 + cy(i)*cy(i)*dx(i) end do cf0 = 0.5*cf0/d if (idbug .ne. 0) then write(lutest,'(A, e13.6)') '/ cf0=', cf0 endif !Check for convergence rel = abs(cf0 - cfsave)/abs(cfsave) if (rel .lt. thresh)then converged=.true. exit endif cfsave = cf0 !Calculate values of Lagrange multipliers gm(n-1) = cy(n-1)*dx(n-1)/d if (n .gt. 3) then do i=n-2,2,-1 gm(i) = cy(i)*dx(i)/d - gm(i+1)*dx(i)/dx(i+1) end do endif !Calculate gradient of cost function with respect to c1 dsdc1 = dx(1)*(cy(1)/d - gm(2)/dx(2)) !Adjust cy(1) using trial step ct(1) = cy(1) - slt*dsdc1 !Calculate remaining c values at trial step do i=2,n-1 ct(i) = -ct(i-1)*dx(i-1)/dx(i) & & -df(i-1)/(dx(i)*dx(i-1)) & & +df(i )/(dx(i)*dx(i )) end do !Calculate cost function at trial step cft = 0.0 do i=1,n-1 cft = cft + ct(i)*ct(i)*dx(i) end do cft = 0.5*cft/d !write(6,*) 'dsdc1,cft,cf0',dsdc1,cft,cf0 !Calculate optimal step length and re-adjust cy(1) den = 2.0*((cft-cf0) + slt*dsdc1*dsdc1) if (den .ne. 0.0) then slo = dsdc1*dsdc1*slt*slt/den else slo =0.0 endif !Adjust slo if desired slo = 1.0*slo cy(1) = cy(1) - slo*dsdc1 if (idbug .ne. 0) then write(lutest,'(A, i4, 3(A,e11.4))') ' Iteration=',k, & &' cf1=', cft,' slt=',slt,' slo=',slo do j=1,n-1 write(lutest,'(A,i2, a, f8.4)') ' i=',j,' c=',cy(j) end do endif !Calculate trial step for next time step slt = 0.5*slo end do !Iteration did not converge if (.not. converged)then ierr=3 fi=-99.9 return endif if (idbug .ne. 0) then write(lutest,'(A)') '/ Iteration converged' endif !Calculate b and a coefficients do i=1,n-1 by(i) = df(i)/dx(i) - cy(i)*(x(i+1) + x(i)) ay(i) = f(i) - by(i)*x(i) - cy(i)*x(i)*x(i) end do return end subroutine init_int subroutine xint(x,f,n,iflag,lflag,xi,fi,ierr) !This routine applies a quadratic interpolation procedure !to f(x) between x(1) and x(n). f(x) is assumed to be !represented by quadratic polynomials between the points !x(i). The polynomials are chosen so that they equal f(i) !at the points x(i), the first derviatives on either !side of the interior x(i) match at x(i), and the second !derivative of the approximated function integrated !over the domain is minimized. !This verision is for interpolating latitude !Input: x(1),x(2) ... x(n) The x values (must be sequential) !f(1),f(2) ... f(n) The function values !n The number of x,f pairs !iflag Flag for initialization !=1 for coefficient calculation !=0 to use previous coefficients !lflag Flag for linear interpolation !=0 to perform linear interpolation !=1 to perform quadratic interpolation !xi The x value at which to interpolate f !Output: fi The interpolated function value !ierr Error flag !=0 Normal return !=1 Parameter nmax is too small !=2 The x values are not sequential !=3 Coefficient iteration did not !converge !=4 Mix-up finding coefficients !=5 if xi .gt. x(n) or .lt. x(1), !xi is set to nearest endpoint !before the interpolation !Note: fi is set to -99.9 if !ierr=1,2,3 or 4 integer, parameter :: nmax=100 real, dimension(n), intent(in) :: x,f integer, intent(in) :: n, iflag, lflag real, intent(inout) :: xi real, intent(out) :: fi integer, intent(out) :: ierr !Save variables real, dimension(nmax) :: ax, bx, cx !Temporary local variables integer :: ii, i logical :: interval_found common /xsave/ ax,bx,cx !Initialize error flag ierr = 0 fi = -99.9 interval_found = .false. !Check to make sure nmax is large enough, and n is .gt. 1 if (n .gt. nmax .or. n .lt. 2) then ierr=1 return endif if (iflag .eq. 1) then ! Initialize the function weights call init_int(x,f,n,lflag,fi,ierr, ax, bx, cx) if (ierr .ne. 0) return endif !Interpolate the function !Check for xi out of bounds if (xi .lt. x(1)) then xi = x(1) ierr = 5 endif if (xi .gt. x(n)) then xi = x(n) ierr = 5 endif !Find the interval for the interpolation ii = 1 do i=2,n if (xi .le. x(i)) then ii = i-1 interval_found = .true. exit endif end do if (.not. interval_found) then ierr = 4 return endif fi = ax(ii) + bx(ii)*xi + cx(ii)*xi*xi return end subroutine xint subroutine linfill(f_arr, f_len, missvalue, fillvalue, inmaxgap) integer, intent(in) :: f_len, inmaxgap real, intent(inout), dimension(f_len) :: f_arr real, intent(in) :: fillvalue, missvalue real :: srbc, erbc, progress integer :: j, k, space, maxgap logical extrapolate extrapolate = (inmaxgap < 0) maxgap = abs(inmaxgap) srbc = fillvalue do j=1, f_len if (f_arr(j) .eq. missvalue) then ! if missing erbc = fillvalue do k=j+1, f_len ! find next non missing within 36hrs space = (k-j) ! how big is the gap so far if (f_arr(k) .ne. missvalue) then erbc = f_arr(k) ! hold onto non-missing value exit endif enddo do k=j, j+space ! now interpolate across the gap progress = (k-(j-1)) / (real(space)+1.) f_arr(k) = (srbc * (1-progress)) + (erbc * progress) enddo ! if we started at the begining of the array, backfill value if (space .ge. maxgap .and. j .eq. 1 .and. extrapolate) then f_arr(:space) = erbc endif ! if we reached the end of the array then extrapolate or fill if (space .ge. maxgap .and. k .ge. f_len) then if (extrapolate) then f_arr(j:) = srbc ! propagate last non missing through end else f_arr(j+space:) = fillvalue endif exit endif else ! not missing, hold onto value and continue srbc = f_arr(j) endif enddo where(f_arr .eq. missvalue) f_arr = fillvalue end subroutine linfill subroutine lcase(charr, lenr) ! This routine converts all letters in the character charr ! of length lent to lower case implicit none ! input and input variables integer, intent(in) :: lenr character (len=*), intent(inout) :: charr ! local variables integer m, k Character(26), Parameter :: ucl = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' Character(26), Parameter :: lcl = 'abcdefghijklmnopqrstuvwxyz' ! Capitalize each letter if it is lowecase do m = 1, lenr k = INDEX(ucl, charr(m:m)) if (k > 0) charr(m:m) = lcl(k:k) end do return end subroutine lcase subroutine uppcase(charv, lenr) ! This routine converts all the letters in ! the character variable charv to upper case. nchar is the number ! of characters to convert. ! implicit none ! input and input variables integer, intent(in) :: lenr character (len=*), intent (inout) :: charv ! local variables integer m, k Character(26), Parameter :: ucl = 'ABCDEFGHIJKLMNOPQRSTUVWXYZ' Character(26), Parameter :: lcl = 'abcdefghijklmnopqrstuvwxyz' ! Capitalize each letter if it is lowecase do m = 1, lenr k = INDEX(lcl, charv(m:m)) if (k > 0) charv(m:m) = ucl(k:k) end do return end subroutine uppcase real function ran2(idum) implicit none integer idum,IM1,IM2,IMM1,IA1,IA2,IQ1,IQ2,IR1,IR2,NTAB,NDIV real AM,EPS,RNMX parameter (IM1=2147483563,IM2=2147483399,AM=1./IM1,IMM1=IM1-1, & & IA1=40014,IA2=40692,IQ1=53668,IQ2=52774,IR1=12211,IR2=3791, & & NTAB=32,NDIV=1+IMM1/NTAB,EPS=1.2e-7,RNMX=1.-EPS) integer idum2,j,k,iv(NTAB),iy,idumtmp save iv,iy,idum2 data idum2/123456789/, iv/NTAB*0/, iy/0/ idumtmp=idum if (idumtmp.le.0) then idumtmp = max (-idumtmp , 1) idum2 = idumtmp do j=NTAB+8,1,-1 k=idumtmp/IQ1 idumtmp=IA1*(idumtmp-k*IQ1)-k*IR1 if (idumtmp.lt.0) idumtmp=idumtmp+IM1 if (j.le.NTAB) iv(j)=idumtmp end do iy=iv(1) endif k=idumtmp/IQ1 idumtmp=IA1*(idumtmp-k*IQ1)-k*IR1 if (idumtmp.lt.0) idumtmp=idumtmp+IM1 k=idum2/IQ2 idum2=IA2*(idum2-k*IQ2)-k*IR2 if (idum2.lt.0) idum2=idum2+IM2 j=1+iy/NDIV iy=iv(j)-idum2 iv(j)=idumtmp if(iy.lt.1)iy=iy+IMM1 ran2=min(AM*iy,RNMX) return end function ran2 subroutine randn(irand,ndr) ! This routine calcuates a random number (irand) ! with ndr digits. On the first call to rand, ! irand should be initialized to -1. ! (ndr=1 gives a number from 0 to 9, ! (ndr=2 gives number from 0 to 99, etc) ! ndr is limited to a value from 1 to 6 implicit none integer, intent(out) :: irand integer, intent(in) ::ndr integer ndrt real ff, fmax ndrt = ndr if (ndr .lt. 1) ndrt = 1 if (ndr .gt. 6) ndrt = 6 ! maximum of the random number generator fmax=10.**ndrt ! calculate a random number between 1 and fmax ff = ran2(1) irand = int(ff*fmax) return end subroutine randn subroutine sindex(i, ni, j, nj) ! This routine scales the number i, which lies between 1 and ni, ! to a number j, which lies between 1 and nj. implicit none integer, intent(in) :: i, ni, nj integer, intent(out) :: j real dj, b, s if (ni .eq. nj) then j = i elseif (ni .lt. nj) then dj = float(nj-1)/float(ni) j = ifix( 1.499 + (float(i)-0.5)*dj ) elseif (ni .gt. nj) then b = float(ni-nj)/float(ni-1) s = float(nj- 1)/float(ni-1) j = ifix(0.499 + b + s*float(i)) endif return end subroutine sindex subroutine tolck(x1,x2,tol,idiff) ! This routine calculates the difference between ! x1 and x2. If the absolute value of the difference ! is less than or equal to tol, then idiff=0. If not, ! then idiff is set to 1 if x2-x1 is positive for -1 if ! x2-x1 is negative. implicit none ! In/Out Variables real, intent(in):: x1, x2, tol integer, intent(out) :: idiff ! Local variables real :: diff, adiff diff = x2-x1 adiff = abs(diff) if (adiff .le. tol) then idiff=0 return endif if (diff .gt. 0.0) then idiff = 1 else idiff = -1 endif return end subroutine tolck subroutine ofti(iego,nx,ny,i,j,ij) ! This routine finds the 1-D index (ij) from a regularly spaced ! 2-D where i goes from 1 to nx and j from 1 to ny. It is assumed that ! x-points were put in the 1-D grid first. Also, 4 options are allowed ! for the starting point of the 2-D grid, depending on the iego value, ! as follows: ! iego = 1 Lower Left ! = 2 Upper Left ! = 3 Upper Right ! = 4 Lower Right implicit NONE integer, intent(in) :: iego, nx, ny, i, j integer, intent(out) :: ij if (iego .eq. 1) then ij = i + (j-1)*nx elseif (iego .eq. 2) then ij = i + (ny-j)*nx elseif (iego .eq. 3) then ij = 1 + (nx-i) + (ny-j)*nx elseif (iego .eq. 4) then ij = 1 + (nx-i) + (j-1)*nx else return endif return end subroutine ofti subroutine rint(azdeg,rad,radt) ! ! The routine interpolates between the wind radii (ne,se,sw,nw) ! in the array rad to determine the radius at the angle azdeg. ! The angle azdeg is measured counter-clockwise from east. ! implicit none real, intent(in) :: azdeg, rad(4) real, intent(out) :: radt real :: aindx, wt1, wt2, radp(0:5), arel integer :: indx ! Check valid range of azdeg if (azdeg .lt. 0 .or. azdeg .gt. 360.) then write(6,'(A, e11.4)') & & '/ Unexpected value of azdeg in routine rint: ', azdeg stop endif ! Make a periodic rad array that increases counter-clockwise radp(0) = rad(2) radp(1) = rad(1) radp(2) = rad(4) radp(3) = rad(3) radp(4) = rad(2) radp(5) = rad(1) ! Convert azdeg to an array index aindx = 4.0*(azdeg+45.0)/360.0 indx = ifix(aindx) arel = aindx - float(indx) wt1 = 1.0 - arel wt2 = arel radt = wt1*radp(indx) + wt2*radp(indx+1) return end subroutine rint end module windspeed_routines module windspeed_dates implicit none contains subroutine jday(imon,iday,iyear,julday) ! This routine calculates the Julian day (julday) from ! the month (imon), day (iday), and year (iyear). The ! appropriate correction is made for leap year. implicit none integer, intent(in) :: imon, iday, iyear integer, intent(out) :: julday integer, dimension(12) :: ndmon integer i ! Specify the number of days in each month ndmon(1) = 31 ndmon(2) = 28 ndmon(3) = 31 ndmon(4) = 30 ndmon(5) = 31 ndmon(6) = 30 ndmon(7) = 31 ndmon(8) = 31 ndmon(9) = 30 ndmon(10) = 31 ndmon(11) = 30 ndmon(12) = 31 ! Correct for leap year if (mod(iyear,4) .eq. 0) ndmon(2)=29 ! Check for illegal input if (imon .lt. 1 .or. imon .gt. 12) then julday=-1 return endif if (iday .lt. 1 .or. iday .gt. ndmon(imon)) then julday=-1 return endif ! Calculate the Julian day julday = iday if (imon .gt. 1) then do i=2,imon julday = julday + ndmon(i-1) end do endif return end subroutine jday subroutine jdayi(julday,iyear,imon,iday) ! This routine calculates the month (imon) and day (iday) ! from the Julian day (julday) and year (iyear). ! The appropriate correction is made for leap year. integer, intent(in) :: julday, iyear integer, intent(out) :: imon, iday integer, dimension(12) :: ndmon integer, dimension(13) :: nsum integer mxjul, i ! Specify the number of days in each month ndmon(1) = 31 ndmon(2) = 28 ndmon(3) = 31 ndmon(4) = 30 ndmon(5) = 31 ndmon(6) = 30 ndmon(7) = 31 ndmon(8) = 31 ndmon(9) = 30 ndmon(10) = 31 ndmon(11) = 30 ndmon(12) = 31 ! Correct for leap year if (mod(iyear,4) .eq. 0) ndmon(2)=29 ! Check for illegal input if (mod(iyear,4) .eq. 0) then mxjul = 366 else mxjul = 365 endif if (julday .lt. 1 .or. julday .gt. mxjul) then imon = -1 iday = -1 return endif ! Calculate the month and day nsum(1) = 0 do i=1,12 nsum(i+1) = nsum(i) + ndmon(i) end do do i=2,13 if (julday .le. nsum(i)) then imon = i-1 exit endif end do iday = julday - nsum(imon) return end subroutine jdayi subroutine tadd(iyr,imon,iday,itime,ihra,iyra,imona,idaya,itimea) ! This subroutine calculates the year,month,day,time ! (iyra,imona,idaya,itimea) that are ihra hours after ! the input year,month,day,time (iyr,imon,iday,itime). ! ! ihra can not exceed the number of hours in 9 months. ! ! If ihra is negative, the year,month,day,time that are abs(ihra) ! before the input year,month,day,time are calculated. ! ! The Julian day utilities jday and jdayi are called ! by this routine. ! implicit none integer, intent(in) :: iyr, imon, iday, itime, ihra integer, intent(out) :: iyra, imona, idaya, itimea integer :: idayex, jdate, jmax iyra = iyr imona = imon idaya = iday ! if (ihra .gt. 0) then !** Start calculation for positive ihra ** !Add the hours to the input time itimea = itime+ihra if (itimea .lt. 24) return !Calculate the number of extra days in the itimea variable !and subtract them from itimea idayex = itimea/24 itimea = itimea - idayex*24 !Calculate Julian day of input date !and add extra days to it call jday(imon,iday,iyr,jdate) jdate = jdate + idayex !Check to see if year has changed call jday(12,31,iyr,jmax) if (jdate .gt. jmax) then iyra = iyr+1 jdate = jdate-jmax else iyra = iyr endif !Calculate month and day corresponding to the !increased jdate call jdayi(jdate,iyra,imona,idaya) else !** Start calculation for negative ihra ** !Add the hours to the input time itimea = itime+ihra if (itimea .ge. 0) return !Calculate adjusted time idayex = -1 + itimea/24 itimea = itimea - idayex*24 if (itimea .eq. 24) then itimea=0 idayex = idayex + 1 endif !Calculate Julian day of input date !and subtract extra days from it call jday(imon,iday,iyr,jdate) jdate = jdate + idayex !Check to see if year has changed if (jdate .lt. 1) then iyra = iyra-1 jdate = jdate + 365 if (mod(iyr-1,4) .eq. 0) jdate = jdate + 1 endif !Calculate month and day corresponding to the !decreased jdate call jdayi(jdate,iyra,imona,idaya) end if return end subroutine tadd end module windspeed_dates module windspeed_cliper implicit none real, parameter, private :: pi = 3.14159265 real, parameter, private :: dtr = 0.01745329 real, parameter, private :: lcor = 0.8 real, parameter, private, dimension(4) :: phi = (/45.0, 315.0, 225.0, 135.0/) contains subroutine radcal(rlat,vmax,asym,beta,th0,rmax,xt,& & r34,r50,r64,r100,ap,rmwp,thp) ! This routine calculates the wind radii from the parameteric vortex real, intent(out), dimension(4) :: r34,r50,r64,r100 real, intent(out) :: ap, rmwp, thp real, intent(in) :: rlat, vmax, asym, beta, th0, rmax, xt real :: x, xi, a, rm,v34, v50, v64, v100, cfac, rotang, sgnlat real :: theta integer :: k ap=0.0 rmwp=0.0 thp=0.0 x = xt if (x .lt. 0.2) x=0.2 if (x .gt. 1.2) x=1.2 xi = 1.0/x a = asym rm = rmax v34 = 34.0 v50 = 50.0 v64 = 64.0 v100 = 100.0 if (vmax .ge. 34.0 .and. vmax .le. 36. ) v34 = 33.0 if (vmax .ge. 49.0 .and. vmax .le. 51. ) v50 = 48.0 if (vmax .ge. 63.0 .and. vmax .le. 65. ) v64 = 63.0 if (vmax .ge. 99.0 .and. vmax .le. 101.) v100 = 98.0 if (rlat .lt. 0.0) then rotang = -90.0 sgnlat = -1.0 else rotang = 90.0 sgnlat = 1.0 endif do k=1,4 theta = (rotang + phi(k) - beta) cfac = cos(dtr*(theta-sgnlat*th0)) r34(k) = rm*( (vmax-a)/( v34-a*cfac) )**xi r50(k) = rm*( (vmax-a)/( v50-a*cfac) )**xi r64(k) = rm*( (vmax-a)/( v64-a*cfac) )**xi r100(k) = rm*( (vmax-a)/(v100-a*cfac) )**xi if ( r34(k) .lt. rm) r34(k) = 0.0 if ( r50(k) .lt. rm) r50(k) = 0.0 if ( r64(k) .lt. rm) r64(k) = 0.0 if (r100(k) .lt. rm) r100(k) = 0.0 ! added by J. Knaff 012904 (atlantic has a small rmw now) if (vmax .lt. 34.)then r34(k)=0.0 else if (vmax .lt. 48.)then r50(k)=0.0 else if (vmax .lt. 63.)then r64(k)=0.0 else if (vmax .lt. 98.)then r100(k)=0.0 endif enddo ! return th0, a and rmw for the calcuation of the r inside the radius of ! maximum winds ap = a rmwp = rm thp = th0 return end subroutine radcal subroutine radcal_land(rlat,vmax,asym,beta,th0,rmax,xt,& & r34,r50,r64,r100) ! This routine calculates the wind radii from the parameteric vortex ! assuming TC is over land (using reduction factor, lcor) real, dimension(4), intent(out) :: r34,r50,r64,r100 real :: rlat, vmax, asym, beta, th0, rmax, xt, ap, rmwp, thp real :: x, xi, a, rm,v34, v50, v64, v100, cfac, rotang, sgnlat real :: theta integer :: k ap=0.0 rmwp=0.0 thp=0.0 x = xt if (x .lt. 0.2) x=0.2 if (x .gt. 1.2) x=1.2 xi = 1.0/x a = asym rm = rmax v34 = 34.0/lcor v50 = 50.0/lcor v64 = 64.0/lcor v100 = 100.0/lcor ! if (vmax .ge. 34.0 .and. vmax .le. 36. ) v34 = 33.0 ! if (vmax .ge. 49.0 .and. vmax .le. 51. ) v50 = 48.0 ! if (vmax .ge. 63.0 .and. vmax .le. 65. ) v64 = 63.0 ! if (vmax .ge. 99.0 .and. vmax .le. 101.) v100 = 98.0 if (rlat .lt. 0.0) then rotang = -90.0 sgnlat = -1.0 else rotang = 90.0 sgnlat = 1.0 endif do k=1,4 theta = (rotang + phi(k) - beta) cfac = cos(dtr*(theta-sgnlat*th0)) r34(k) = rm*( (vmax-a)/( v34-a*cfac) )**xi r50(k) = rm*( (vmax-a)/( v50-a*cfac) )**xi r64(k) = rm*( (vmax-a)/( v64-a*cfac) )**xi r100(k) = rm*( (vmax-a)/(v100-a*cfac) )**xi if ( r34(k) .lt. rm) r34(k) = 0.0 if ( r50(k) .lt. rm) r50(k) = 0.0 if ( r64(k) .lt. rm) r64(k) = 0.0 if (r100(k) .lt. rm) r100(k) = 0.0 ! added by J. Knaff 012904 (atlantic has a small rmw now) if (vmax .lt. 34.)then r34(k)=0.0 else if (vmax .lt. 48.)then r50(k)=0.0 else if (vmax .lt. 63.)then r64(k)=0.0 else if (vmax .lt. 98.)then r100(k)=0.0 endif enddo ! return th0, a and rmw for the calcuation of the r inside the radius of ! maximum winds ap = a rmwp = rm thp=th0 return end subroutine radcal_land subroutine radcal1(rlat,rlon,vmx,spd,head,vmxth,radarr,& & ap,rmwp,thp,xp) !This routine calculates the wind radii (nmi) ne,nw,se,sw !of the storm center. These values are returned in array !radarr. This routine works in either the northern or southern !hemisphere. real, intent(in) :: rlat, rlon, vmx, spd, head, vmxth real, intent(out) :: ap, rmwp, thp, xp, radarr(4) real :: x, xi, rm, asf, quad, thetar, ct, rota integer :: m ap=0.0 rmwp=0.0 thp=0.0 !Check to see if requested threshold exceeds !current intensity if (vmx .lt. vmxth) then do m=1,4 radarr(m) = 0.0 enddo return endif !Calculate climatological modified Rankine vortex parameters !as a function of maximum winds x = 0.42 + .0025*vmx xi = 1.0/x rm= 54. - .27*vmx if (rm .lt. 15.0) rm = 15.0 !Calculate asymmetry factor if (spd .gt. 0.0) then asf = 0.9*(spd**0.63) else asf = 0.0 endif !Calculate wind radii !Rotation variable for southern hemisphere if (rlat .lt. 0.0) then rota = -90.0 else rota = 90.0 endif do m=1,4 quad = 45.0 + 90.0*float(m-1) thetar = dtr*(head + rota - quad) ct = cos(thetar) radarr(m) = rm*( (vmx-asf)/(vmxth-asf*ct) )**xi if (radarr(m) .lt. rm) radarr(m) = 0.0 enddo !return a and rmw associated with the vortex for the calculation !of the radii of significant winds inside the radius of maximum wind. !Also return x parameter. ap=asf rmwp=rm xp = x return end subroutine radcal1 end module windspeed_cliper module windspeed_model use windspeed_routines use windspeed_cliper, only : radcal, radcal1, radcal_land implicit none real, parameter :: pi=3.141593, dtr=pi/180.0 real, parameter, private :: reflat = 25 !Set cmta to convert the model wind radii to the average !in a quadrant rather than the maximum in a quadrant. !For no correction, set cmta=1.0 real, parameter, private:: cmta = 0.85 !Set rbuf to the minmum ratio of the rmw and smallest wind radii real, parameter :: rbuf=0.75 ! this will be overridden if the value is in the namelist character (len=4), private :: coeff_yrs='2024' real :: land_buffer_km = 0 integer :: ioper ! Climo Vortex Variables integer, parameter, private :: nc=13 real, dimension(nc) :: coef character (len=2), private :: lcoef(nc) contains subroutine vpcal(vmax,rlat,c,asym,th0,rmax,x) ! This routine calculates the parameters of the climatological ! vortex, given the max wind, latitutde and storm translational ! speed real, intent(in) :: vmax, rlat, c real, intent(out) :: asym, th0, rmax, x real t0, t1, t2, a0, a1, a2, a3, x0, x1, x2, r0, r1, r2 real rlatp, rm, a t0 = coef( 1) t1 = coef( 2) t2 = coef( 3) a0 = coef( 4) a1 = coef( 5) a2 = coef( 6) a3 = coef( 7) x0 = coef( 8) x1 = coef( 9) x2 = coef(10) r0 = coef(11) r1 = coef(12) r2 = coef(13) ! Calculate climatological vortex parameters rlatp = abs(rlat)-reflat rm = r0 + r1*vmax + r2*rlatp x = x0 + x1*vmax + x2*rlatp a = a0 + a1*c + a2*c*c + a3*rlatp th0 = t0 + t1*rlatp + t2*c ! if (rm .lt. 10.0) rm = 10.0 if (rm .lt. 15.0) rm = 15.0 if (x .lt. 0.2) x = 0.2 rmax = rm asym = a return end subroutine vpcal subroutine landcross(evmxs,edtls,nt,mti,dts,evmxs_new) ! Passed variables and arrays implicit none integer :: nt, mti real :: dts real, dimension(0:mti) :: evmxs real, dimension(0:mti) :: edtls real, dimension(0:mti) :: evmxs_new real :: rf, vb, al, m, vmxp integer :: k, kk, kk1, kk2, ilst, ilnd, tland ! Specify parameters for exponential inland decay model rf = 0.9 vb = 26.7 al = 0.095 do kk = 0, mti evmxs_new(kk) = evmxs(kk) enddo do k = 0, nt kk1 = 12*k kk2 = 12*(k+1) ! Find out if there is a land point between 2 over water points ilst=0 ilnd=0 if (edtls(kk1) .gt. 0. .and. edtls(kk2) .gt. 0. .and. & & minval(edtls(kk1+1:kk2-1)) .lt. 0.) then do kk = kk1+1, kk2-1 if (edtls(kk) .lt. 0) then ilst=kk exit endif enddo do kk = kk2-1,kk1+1,-1 if (edtls(kk) .lt. 0) then ilnd=kk exit endif enddo do kk = kk1, ilst-1 evmxs_new(kk) = evmxs(kk1) enddo vmxp = evmxs(ilst-1) do kk = ilst,ilnd tland = kk - ilst ! Adjust forecast intensity using the ! Kaplan and DeMaria (JAM, 1995) exponential decay model evmxs_new(kk) = vb + (rf*vmxp - vb)*exp(-al*tland) enddo m=(evmxs(kk2)-evmxs_new(ilnd))/real(kk2-ilnd) do kk = ilnd+1, kk2 evmxs_new(kk) = m*dts + evmxs_new(kk-1) enddo endif enddo return end subroutine landcross logical function is_point_in_windfield(point_lat, point_lon, & & storm_lat, storm_lon,ap, betat, thp, rmwp, dtl1, vmx1,vmxth, & & radarr_land, radarr ) real, intent(in) :: point_lat, point_lon, storm_lat, storm_lon real, intent(in) :: radarr_land(4), radarr(4) real, intent(in) :: ap, betat, thp, rmwp, dtl1, vmx1,vmxth real :: dtlv, rnm, azdeg, radt, radi, rlat2, rlon2, rlat1, rlon1 !Check to see if this grid point for this ensemble !member has already been counted. If so, there is no !need to check it again. !if (wprobs(ii,ntsave) .eq. 1.0) go to 1000 !Note: This check was disabled for the wprobi calculation. rlat2 = point_lat rlon2 = point_lon rlat1 = storm_lat rlon1 = storm_lon is_point_in_windfield = .false. !Calculate distance to land for lat/lon point call gdland_table(rlon2,rlat2,dtlv) !If lat/lon point is over land and storm center is over !water, use radii that take into account overland !reduction in wind speeds due to friction call rcal(rlat1,rlon1,rlat2,rlon2,rnm,azdeg) if (dtlv .le. land_buffer_km) then call rint(azdeg,radarr_land,radt) else call rint(azdeg,radarr,radt) endif !Check to see if interpolated radius is less than !the radius of maximum wind. If so, set it to slightly less than zero !so this point is excluded from the "yes" area. if (radt .lt. rmwp) return !Check if r is within outer radius if (rnm .gt. radt) return if (radt .gt. rmwp .and. vmx1 .gt. ap) then !Calculate the radius inside the radius of maximum wind radi=(vmxth-ap*cos((90.0+azdeg-betat-thp)*dtr)) & & *rmwp/(vmx1-ap) if (radi .gt. rmwp) radi=rbuf*rmwp else radi = 0.0 endif !check if r is within inner radius if (rnm .lt. radi) return is_point_in_windfield = .true. return end function is_point_in_windfield subroutine ptadcal(nts,nte,ntoap,itarad,ftimes,toap,icnt,ptd) ! This routine calculates the probability distributions of the time ! of arrival or departure from the times accumulated in itarad. implicit none ! Input Variables integer, intent(in) :: nts, nte, ntoap, itarad(0:nte) real, intent(in) :: ftimes(0:nte), toap(ntoap) ! Output Variables integer, intent(out) :: icnt real, intent(inout) :: ptd(ntoap) ! local variables integer kk, n, icumt real probt ! Calculate total number of realizations contributing to the ! probability calculation. icnt = 0 do kk=nts,nte icnt = icnt + itarad(kk) enddo ptd(:) = -99. if (icnt .lt. 1) then icnt = -99 return endif icumt = 0 do kk=nts,nte icumt = icumt + itarad(kk) probt = 100.0*float(icumt)/float(icnt) do n=1,ntoap if (ptd(n) .lt. 0.0 .and. probt .ge. toap(n)) then ptd(n) = ftimes(kk) endif enddo enddo return end subroutine ptadcal subroutine idfile(flat,flon,fvmx,ftime,elat,elon,evmx,& & mte,ntr,mt,nt,wtprob,ibasin) !**************************************************************** ! This routine calcuates intensity PDFs and writes related ! information to a file if requested. !**************************************************************** implicit none real, dimension(0:mt), intent(in) :: flat,flon,fvmx,ftime real, dimension(mte, 0:mt), intent(in) :: elat,elon,evmx real, dimension(7, 9), intent(out) :: wtprob integer, intent(in) :: mt, mte, ibasin, nt,ntr ! ++ Local variables for intensity PDFs integer, parameter :: mtl=50, mii=9 real :: picnt(mtl, mii), vi(mii), phurr(mtl), dvi, vs integer :: iicnt(mtl, mii),isum(mtl),iwtprob(7,9),iwptime(7),ivmx(7) integer :: iprtf, iprtd, iprtn, iprte, ihurr integer :: nv, nvi, m, kk, n, luidst, i, ii, k character (len=8) :: fmt_hour(7), tmpstr character (len=17) :: wplab(9) character (len=35) :: luidst_fmt ! Set iprtf=1 to print forecast track to mcidst.dat file iprtf=0 ! Set iprtd=1 to print intensity PDFs to mcidst.dat iprtd=0 ! Set iprtn=1 to print NHC wind speed probability table to mcidst.dat iprtn=0 ! Set iprte=1 to print all lat/lon/vmax of all realizations to mcidst.dat iprte=0 ! Intensity PDF calculation dvi = 10.0 vs = -150.0 nvi = mii nv = nvi-1 vi = (/-300.0, 14.5, 33.5, 63.5, 82.5, 95.5, 113.5, 135.5, 300.0/) ! Set ihurr to the index of the hurricane threshold ! to sum all probabilities for hurricane intensity. Set ! ihurr=0 to skip this option. ihurr = 4 iicnt = 0 picnt = 0.0 phurr = 0.0 isum = 0.0 do k=1,nt do i=1,nv do n=1,ntr if (evmx(n,k) .ge. vi(i) .and. evmx(n,k) .lt. vi(i+1)) then iicnt(k,i) = iicnt(k,i) + 1 isum(k) = isum(k) + 1 endif enddo enddo enddo do k=1,nt do i=1,nv picnt(k,i) = 100.0*float(iicnt(k,i))/float(ntr) enddo enddo if (ihurr .gt. 0) then do k=1,nt do i=ihurr,nv phurr(k) = phurr(k) + picnt(k,i) enddo enddo endif if (iprtf .eq. 1 .or. iprtd .eq. 1 .or. & & iprtn .eq. 1 .or. iprte .eq. 1 ) then luidst = 61 open(file='mcidst.dat',unit=luidst,form='formatted', status='replace') endif ! iyr2 = iyr - 100*(iyr/100) if (iprtf .eq. 1) then ! Write forecast track and intensity write(luidst,'(i6,1x,i3)') ntr,nt ! Write control forecast do k=0,nt write(luidst,'(f5.0,1x,f6.1,1x,f6.1,1x,f6.1)') & & ftime(k),flat(k),flon(k),fvmx(k) enddo endif if (iprtd .eq. 1) then ! Write the intensity PDFs write(luidst,'(/, A)') 'Intensity PDF info' do k=1,nt luidst_fmt = '(A,f6.1,A,f6.1,A,i3,A,i5)' write(luidst,luidst_fmt) 't=',ftime(k),' fcst v=',fvmx(k), & & ' nint=', nv,' isum=', isum(k) enddo write(luidst,'(/,18x,11(f6.1,1x))') (ftime(k),k=1,nt,1) do i=1,nv luidst_fmt = "(f6.1,A,f6.1, 11(1x,f6.1))" write(luidst,luidst_fmt) & & vi(i),' to ',vi(i+1),(picnt(k,i),k=1,nt,1) enddo write(luidst,'(A)') ' ' write(luidst,'(A)') ' ' endif ! Calculate the probabilities for the NHC wind probability table wplab(1) = 'DISSIPATED ' wplab(2) = 'TROP DEPRESSION ' wplab(3) = 'TROPICAL STORM ' wplab(4) = 'HURRICANE ' wplab(5) = 'HUR CAT 1 ' wplab(6) = 'HUR CAT 2 ' wplab(7) = 'HUR CAT 3 ' wplab(8) = 'HUR CAT 4 ' wplab(9) = 'HUR CAT 5 ' kk = 0 do k=1,10 if (k .eq. 5) cycle if (k .eq. 7) cycle if (k .eq. 9) cycle kk = kk + 1 iwptime(kk) = ifix(ftime(k)) ivmx(kk) = ifix(fvmx(k)) ii = 0 do i=1,8 ii = i if (i .ge. 4) ii = i + 1 wtprob(kk,ii) = picnt(k,i) iwtprob(kk,ii) = ifix(0.5 + picnt(k,i)) enddo wtprob(kk,4) = phurr(k) iwtprob(kk,4) = ifix(0.5 + phurr(k)) enddo if (iprtn .eq. 1) then ! Write the probabilities out in a format to match the NHC ! wind probability table write(luidst,'(A)') ' - - - MAXIMUM WIND SPEED (INTENSITY) PROBABILITIES - - -' write(luidst,'(A)') ' ' do k=1, 7 write(tmpstr, "(i3)") iwptime(k) write(fmt_hour(k),"(3A)") "(",trim(adjustl(tmpstr)),")" end do luidst_fmt = '(A, 7(8A))' write(luidst,luidst_fmt) 'FORECAST HOUR ',(adjustr(fmt_hour(k)),k=1,7) write(luidst,'(A,7(i3,5x))') 'FCST MAX WIND KT ', (ivmx(k),k=1,7) write(luidst,'(72A)') repeat('- ',36) do i=1,4 write(luidst,'(a17,1x,7(i3,5x))') wplab(i),(iwtprob(k,i),k=1,7) enddo write(luidst,'(72A)') repeat('- ',36) do i=5,9 write(luidst,'(a17,1x,7(i3,5x))') wplab(i),(iwtprob(k,i),k=1,7) enddo write(luidst,'(A)') repeat('- ',36) write(luidst,'(A)') ' ' write(luidst,'(A)') ' ' ! write(luidst,305) (wplab(i),i=1,9) ! 305 format('Time ',9(a4,2x)) ! do k=1,7 ! write(luidst,310) iwptime(k),(iwtprob(k,i),i=1,9) ! 310 format(1x,i3,2x,9(1x,i3,2x)) ! enddo endif if (iprte .eq. 1) then ! Write track and intensities of all ensemble members do k=1,nt write(luidst,'(A2,f6.1)') 't=', ftime(k) do m=1,ntr write(luidst,'(f7.2,1x,f7.2,1x,f7.2)') & & elat(m,k),elon(m,k),evmx(m,k) enddo enddo endif if (iprtf .eq. 1 .or. iprtd .eq. 1 .or. & & iprtn .eq. 1 .or. iprte .eq. 1 ) then close(luidst) endif return end subroutine idfile subroutine atfile(flat,flon,fvmx,elat,elon,evmx,& & iyr,imon,iday,itime,mte,ntr,mt,ibasin) !**************************************************************** !This routine writes the ensemble forecasts to a file in a !modified version of the old-ATCF A-deck format !**************************************************************** use windspeed_dates, only : tadd implicit none real, dimension(0:mt), intent(in) :: flat,flon,fvmx real, dimension(mte,0:mt), intent(in):: elat,elon,evmx integer, intent(in) :: iyr, imon, iday, itime, mte, ntr, mt, ibasin !++ Local variables character *1 echar character (len=50) fmtstr_data, fmtstr_lng, fmtstr_shr integer, dimension(10) :: ilat,ilon,ivmx real :: dlat, dlon integer :: luatcf, iyr2, ilat0, ilon0, ivmx0, ihra integer :: ilatm, ilonm, ivmxm, idaym, iyrm, imonm, itimem integer :: k, i, ii luatcf = 60 open(file='mcatcf.dat',unit=luatcf,form='formatted',status='replace') iyr2 = iyr - 100*(iyr/100) if (ibasin.eq.1)then ilat0 = ifix( 10.0*flat(0)) ilon0 = iabs(ifix(10.*(flon(0)-360.))) ivmx0 = ifix( fvmx(0)) elseif (ibasin.eq.2)then ilat0 = ifix( 10.0*flat(0)) ilon0 = iabs(ifix(10.*(flon(0)-360.))) ivmx0 = ifix( fvmx(0)) elseif (ibasin.eq.3)then ilat0 = ifix( 10.0*flat(0)) ilon0 = iabs(ifix(10.0*flon(0))) ivmx0 = ifix( fvmx(0)) elseif (ibasin.eq.4)then ilat0 = ifix(10.0*flat(0)) if (flon(0).gt.180.)then ilon0 = iabs(ifix(10.0*(flon(0)-360.))) else ilon0 = ifix(10.0*flon(0)) endif ivmx0 = ifix( fvmx(0)) endif !Estimate positions at t=-12 h ihra = -12 call tadd(iyr2,imon ,iday ,itime ,ihra,iyrm,imonm,idaym,itimem) dlat = flat(1)-flat(0) dlon = flon(1)-flon(0) ilatm = ifix( 10.0*(flat(0) - dlat)) ilonm = ifix(-10.0*(flon(0) - dlon)) ivmxm = ivmx0 k = 0 fmtstr_lng = '(a6, 4(i2.2),6(i4),6(i3))' fmtstr_shr = '(a6, 4(i2.2))' write(luatcf,fmtstr_lng) '00COMA', iyrm,imonm,idaym,itimem,ilatm,ilonm, & & k,k,k,k,ivmxm,k,k,k,k,k write(luatcf,fmtstr_shr) '00COMB',iyrm,imonm,idaym,itimem write(luatcf,fmtstr_lng) '00COMA',iyr2,imon,iday,itime,ilat0,ilon0, & & k,k,k,k,ivmx0,k,k,k,k,k write(luatcf,fmtstr_shr) '00COMB',iyr2,imon,iday,itime !Write control forecast do k=1,mt if (ibasin.eq.1)then ilat(k) = nint( 10.0*flat(k)) ilon(k) = iabs(nint(10.*(flon(k)-360.))) ivmx(k) = nint( fvmx(k)) elseif (ibasin.eq.2)then ilat(k) = nint( 10.0*flat(k)) ilon(k) = iabs(nint(10.*(flon(k)-360.))) ivmx(k) = nint( fvmx(k)) elseif (ibasin.eq.3)then ilat(k) = nint( 10.0*flat(k)) ilon(k) = nint(-10.0*flon(k)) ivmx(k) = nint( fvmx(k)) elseif (ibasin.eq.4)then ilat(k) = nint(10.0*flat(k)) if (flon(k).gt.180.)then ilon(k) = iabs(nint(10.0*(flon(k)-360.))) else ilon(k) = nint(10.0*flon(k)) endif ivmx(k) = nint( fvmx(k)) endif enddo echar = 'E' ii = 0 fmtstr_data = '(a2,a1,i3.3,4(i2.2),20(i5),1x,10(i5),a5,i2.2)' write(luatcf,fmtstr_data) '99', echar,ii,iyr2,imon,iday,itime, & & (ilat(k),ilon(k),k=1,mt), & & (ivmx(k),k=1,mt),' ALXX', iyr2 !Write ensemble members do i=1,ntr do k=1,mt if (ibasin.eq.1)then ilat(k) = nint( 10.0*elat(i,k)) ilon(k) = iabs(nint(10.*(elon(i,k)-360.))) ivmx(k) = nint( evmx(i,k)) elseif (ibasin.eq.2)then ilat(k) = nint( 10.0*elat(i,k)) ilon(k) = iabs(nint(10.*(elon(i,k)-360.))) ivmx(k) = nint( evmx(i,k)) elseif (ibasin.eq.3)then ilat(k) = nint( 10.0*elat(i,k)) ilon(k) = nint(-10.0*elon(i,k)) ivmx(k) = nint( evmx(i,k)) elseif (ibasin.eq.4)then ilat(k) = nint(10.0*elat(i,k)) if (elon(i,k).gt.180.)then ilon(k) = iabs(nint(10.0*(elon(i,k)-360.))) else ilon(k) = nint(10.0*elon(i,k)) endif ivmx(k) = nint( evmx(i,k)) endif end do if (i .lt. 1000) then echar = 'E' ii = i else echar = 'F' ii = i-1000 endif write(luatcf,fmtstr_data) '99', echar,ii,iyr2,imon,iday,itime, & & (ilat(k),ilon(k),k=1,mt), & & (ivmx(k),k=1,mt),' ALXX', iyr2 end do close(luatcf) return end subroutine atfile subroutine gopck(plon,plat,mxy,nxy) !This routine determines if the grid optimization procedure !can be implemented. The input lon/lat points in the 1-D !arrays must have been derived from a regularly spaced 2-D lon/lat !grid, where the lon values were put into the 1-D array first, followed !by the lat points. Four options are acceptable, where the 1-D points !started in the lower left, upper left, upper right or lower right corner !of the 2-D grid (iego is set to 1, 2, 3 or 4). If plon,plat were not !derived from a regular 2-D grid with x incremented first, iego is set to !zero. Several other grid parameters (dlone,dlate,plone,plate,nxe,nye) !are also calculated for use with the optimzed option. implicit NONE integer, parameter :: mxe=2000,mye=1000 integer :: mxy, nxy real, dimension(mxy) :: plon, plat !Arrays for code optimization for regular lat/lon grid real :: plone(mxe),plate(mye),dlone, dlate integer :: nxe, nye, iego real :: tol, dxf, dyf, z, nxy1, dlones, dlates integer idiff, idiffx, idiffy, i,idx, idy, j, ij save common /egopt/ dlone,dlate,plone,plate,nxe,nye,iego !Default values iego = 0 dlone = 0.0 dlate = 0.0 nxe = 0 nye = 0 !Note: Uncomment the line below to always skip the optimization !return if (nxy .lt. 100) then !Optimization is not worth the effort iego = 0 return endif !Set tolerance for difference checking of real variables tol = 0.001 !Calculate dx and dy of first two 1-D points dxf = plon(2)-plon(1) dyf = plat(2)-plat(1) z = 0.0 call tolck(dxf,z,tol,idiffx) call tolck(dyf,z,tol,idiffy) if (idiffx .eq. 0 .or. idiffy .ne. 0) then !Grid is not indexed in x first or is not equally spaced iego = 0 return endif !Find nxe and nye do i=2,nxy call tolck(plat(1),plat(i),tol,idiff) if (idiff .ne. 0) then nxe = i - 1 exit endif enddo !If you get to here, all the y values are the same, so !there is no underlying 2-D grid. if(nxe .eq. 0)then iego = 0 return endif nye = nxy/nxe !Check to see if grid really is rectangular nxy1 = nxe*nye if (nxy1 .ne. nxy) then iego = 0 return endif !Check to make sure mxe and mye are large enough if (nxe .gt. mxe .or. nye .gt. mye) then iego = 0 return endif !Calculate lat/lons for use in 2-D dlones = plon( 2)-plon(1) dlates = plat(nxe+1)-plat(1) call tolck(z,dlones,tol,idx) call tolck(z,dlates,tol,idy) if (idx .eq. 0 .or. idy .eq. 0) iego = 0 if (idx .gt. 0 .and. idy .gt. 0) iego = 1 if (idx .gt. 0 .and. idy .lt. 0) iego = 2 if (idx .lt. 0 .and. idy .lt. 0) iego = 3 if (idx .lt. 0 .and. idy .gt. 0) iego = 4 if (idx .gt. 0) then do i=1,nxe plone(i) = plon(i) enddo else do i=1,nxe plone(i) = plon(1+nxe-i) enddo endif if (idy .gt. 0) then do j=1,nye plate(j) = plat(1 + (j-1)*nxe) enddo else do j=1,nye plate(j) = plat(1 + nxe*(nye-j)) enddo endif dlone = abs(dlones) dlate = abs(dlates) !Perform final check on 2-D points to make sure !they match all the 1-D points do j=1,nye do i=1,nxe call ofti(iego,nxe,nye,i,j,ij) call tolck(plone(i),plon(ij),tol,idiffx) call tolck(plate(j),plat(ij),tol,idiffy) if (idiffx .ne. 0 .or. idiffy .ne. 0) then iego = 0 return endif enddo enddo return end subroutine gopck character (len=200) function get_parm_dir() implicit none if (ioper .eq. 1) then call getenv ( "ATCFPROBLTY", get_parm_dir ) ! ATCF elseif (ioper .eq. 2) then call getenv ( "PRBLTYINC", get_parm_dir ) ! IBM else get_parm_dir = "include" ! local endif return end function get_parm_dir subroutine xfit(vmax,rlat,c,beta,asym,rmax,xc,x, & & r34,r50,r64,r100,c34,c50,c64,c100,p34,p50,p64,p100, errmin) !This rouinte calculates the climatological !values of the asymmetry factor asym, the radius of !maximum wind (rmax), and the size parameter xc, for !the idealized surface wind model, given !the maximum wind vmax, the latitude rlat, !the storm motion vector and the angle ccw !relative to east (beta) of the storm motion !vector. !The value of x that best fits the observed !wind radii (r34, r50, r64, r100) is then calculated !along with the the predicted wind radii using the !climatological value of the size parameter xc !(c34, c50, c64, c100) and the radii using the optimal !size parameter x (p34, p50, p64, p100) real :: vmax, rlat, c, beta, asym, rmax, xc, x real, intent(in), dimension(4) :: r34,r50,r64,r100 real, intent(inout), dimension(4) :: c34,c50,c64,c100 real, intent(inout), dimension(4) :: p34,p50,p64,p100 real :: errmin ! local variables real, dimension(4) :: pt34,pt50,pt64,pt100 real :: a,ap, dx, err100, err34, err50, err64, errt real :: rlatkp, rm, rmwp, rsum, th0, thp, v100, v34, v50, v64 real :: xi, xl, xmin integer :: i, idbug,imin, k, nx character (len=60) :: dbug_fmt, dbug_array_fmt v34 = 34.0 v50 = 50.0 v64 = 64.0 v100 = 100.0 !Calculate climatological vortex parameters call vpcal(vmax,rlat,c,asym,th0,rmax,x) rlatkp = abs(rlat)-reflat xi = 1.0/x xc = x if (vmax .lt. 34.0) then p34(:) = 0.0 p50(:) = 0.0 p64(:) = 0.0 p100(:) = 0.0 c34(:) = 0.0 c50(:) = 0.0 c64(:) = 0.0 c100(:) = 0.0 x = xc return endif a = asym rm= rmax !Calculate wind radii from climatological vortex !do k=1,4 !theta = (90.0 + phi(k) - beta) !cfac = cos(dtr*(theta-th0)) !c34(k) = rm*( (vmax-a)/(v34-a*cfac) )**xi !c50(k) = rm*( (vmax-a)/(v50-a*cfac) )**xi !c100(k) = rm*( (vmax-a)/(v100-a*cfac) )**xi !if (c34(k) .lt. rm) c34(k) = 0.0 !if (c50(k) .lt. rm) c50(k) = 0.0 !if (c100(k) .lt. rm) c100(k) = 0.0 !enddo call radcal(rlat,vmax,a,beta,th0,rm,xc,c34,c50,c64,c100,ap,rmwp,thp) !Check for special case where all predicted radii are zero rsum = SUM(c34) + SUM(c50) + SUM(c64)+ SUM(c100) if (rsum .le. 0.0) then x = xc p34(:) = 0.0 p50(:) = 0.0 p64(:) = 0.0 p100(:) = 0.0 return endif !Find x that best fits the observed radii !Set search increment dx, starting x value and number of searches dx = 0.01 xl = 0.2 nx = 100 errmin = 1.0e+10 imin = 0 xmin = xc do i=0,nx err34 = 0.0 err50 = 0.0 err64 = 0.0 err100 = 0.0 x = xl + dx*float(i) xi = 1.0/x call radcal(rlat,vmax,a,beta,th0,rm,x,pt34,pt50,pt64,pt100,ap,rmwp,thp) ! do k=1,4 !theta = (90.0 + phi(k) - beta) !cfac = cos(dtr*(theta-th0)) !pt34(k) = rm*( (vmax-a)/(v34-a*cfac) )**xi !pt50(k) = rm*( (vmax-a)/(v50-a*cfac) )**xi !pt100(k) = rm*( (vmax-a)/(v100-a*cfac) )**xi !if (pt34(k) .lt. rm) pt34(k) = 0.0 !if (pt50(k) .lt. rm) pt50(k) = 0.0 !if (pt100(k) .lt. rm) pt100(k) = 0.0 ! end do if (vmax .ge. 34.0) then err34 = 0.25*sum(abs(pt34(:)-r34(:))) endif if (vmax .ge. 50.0) then err50 = 0.25*sum(abs(pt50(:)-r50(:))) endif if (vmax .ge. 64.0) then err64 = 0.25*sum(abs(pt64(:)-r64(:))) endif if (vmax .ge. 100.0) then err100 = 0.25*sum(abs(pt100(:)-r100(:))) endif errt = err34 + err50 + err64 + err100 if (errt .lt. errmin) then errmin = errt imin = i xmin = x p34(:) = pt34(:) p50(:) = pt50(:) p64(:) = pt64(:) p100(:) = pt100(:) endif idbug=0 dbug_fmt = '(A14, 4(f6.1,1x,f6.1,3x))' if (idbug .eq. 1) then write(6,dbug_fmt) ' p34, r34 : ',(pt34(k),r34(k),k=1,4) write(6,dbug_fmt) ' p50, r50 : ',(pt50(k),r50(k),k=1,4) write(6,dbug_fmt) ' p64, r64 : ',(pt64(k),r64(k),k=1,4) write(6,dbug_fmt) ' p100, r100 : ',(pt100(k),r100(k),k=1,4) dbug_array_fmt = '(i4,1x,A,(f7.4,1x,f7.4), A,5(f4.0,1x), A,5(f8.3,1x)' write(6,dbug_array_fmt) i,'xc,x=',xc,x, & & ' vm,rm,c,be,a=',vmax,rmax,c,beta,asym, & & ' err34,50,100,t: ',err34,err50,err64,err100,errt endif end do x = xmin return end subroutine xfit subroutine iens(fvmx,fdtl,ftime,elat,elon,evmx,mte,ntr,mt, nt,& & lulog,ibasin) !This routine calculates the intensity ensemble members by !radomly sampling from the offical intensity error distribtuions. !Separate distributions are used for forecasts over land and !over water. implicit NONE real :: fvmx(0:nt),fdtl(0:nt),ftime(0:nt) real :: elat(mte,0:nt),elon(mte,0:nt),evmx(mte,0:nt) integer :: mte, ntr, mt, nt, lulog, ibasin !Local variables integer, parameter :: mc=6000,mp=6 integer :: nev(mt), ncv(mt) real :: verr(mc,mt) real :: vcoef(mp,mt) real :: edtl(0:mt),fvmxt(0:mt),telon(0:mt),telat(0:mt) real ::vmxp, vmxk, vldiss, verrp, verro, veold, venew, vcap, vb real :: tland, rf,ftimep, ftimek, ftimepp, evold, evnew, edtlt real :: dtlpp, dtlp,dtlk, dldiss, deldtl, al integer :: nrand, net, ndr, nct, n, m, luidi integer :: kp, k, kk, km1, pwk, nwk integer :: iswtl,isltw,iskip, irand, iicap, i, indx, ind character (len=200) :: include_path, fnidi !initialize the coeficient arrays so they are zero if not filled vcoef = 0.0 ! Specify name of file with intensity error distributions include_path = get_parm_dir() ind = index( include_path, " " ) - 1 if ( ibasin .eq. 1 ) then ! Atlantic basin fnidi = trim(include_path)//"/mcidist_al_"//coeff_yrs//".dat" endif if (ibasin .eq. 2)then ! eastern North Pacific basin fnidi = trim(include_path)//"/mcidist_ep_"//coeff_yrs//".dat" endif if (ibasin .eq. 3)then ! western North Pacific and N. Indian basins fnidi = trim(include_path)//"/mcidist_wp_"//coeff_yrs//".dat" endif if (ibasin .eq. 4)then ! Southern Hemisphere basin fnidi = trim(include_path)//"/mcidist_sh_"//coeff_yrs//".dat" endif !Set iskip=1 to skip perturbations of intensity. !For this case, the ensemble intensities are set equal !to the forecast intensity. iskip = 0 !Set iswtl=1 to skip the intensity correction for the case !where the forecast track is over water but the ensemble !member track is over land iswtl=0 !Set isltw=1 to skip the intensity correction for the case !where the forecast track is over land but the ensemble !member track is over water isltw=0 !Set iicap=1 to apply max wind cap as a function of the distance !inland iicap=1 !If iicap=1, set the minimum max wind (kt) before a storm !is considered dissipated (vldiss). If the max wind of an ensemble !member over land drops below this value, the intensity !for the rest of the track will be set to zero. !Also, specify the distance inland (neg km) to implement the !max wind cap (dldiss) vldiss=15.0 dldiss=-20.0 !Initialize random number routine and !specify number of places for random number !variable !irand=-1 ndr=4 nrand= 10**ndr !Specify parameters for exponential inland decay model rf = 0.9 vb = 26.7 al = 0.095 if (iskip .eq. 1) then !Skip perturbation of intensity do i=1,ntr do k=0,nt evmx(i,k) = fvmx(k) enddo enddo return endif !Open and read the intensity error file luidi = 41 open(file=fnidi,unit=luidi,form='formatted',status='old',err=900) do k=1,mt read(luidi,*,err=905,end=905) read(luidi,*,err=905,end=905) nev(k),ncv(k) nct = ncv(k)+1 net = nev(k) read(luidi,'(10(e12.5,1x))',end=905,err=905) (vcoef(n,k),n=1,nct) read(luidi,'(10(f7.1,1x))',end=905,err=905) (verr(i,k),i=1,net) end do close(luidi) !Set t=0 ensemble intensity to forecast intensity do m=1,ntr evmx(m,0) = fvmx(0) enddo do m=1,ntr if (mod(m,1000) .eq. 0) write(6,*) 'iens no. ',m !Calculate the distance to land along the ensemble track !and copy forecast intensity do k=0,nt telon(k) = elon(m,k) telat(k) = elat(m,k) fvmxt(k) = fvmx(k) enddo do k=0,nt call gdland_table(telon(k),telat(k),edtl(k)) end do do k=1,nt !Check to see if intensity requires adjustment !for the case where the forecast track point is over water !but the ensemble track point is over land. if (fdtl(k) .gt. 0.0 .and. edtl(k) .lt. 0.0 & & .and. iswtl .eq. 0) then vmxk = fvmx(k) dtlk = edtl(k) ftimek = ftime(k) !Find first previous ensemble point that was over water !or use t=0 point if all ensemble points are over land kp = 0 vmxp = fvmx(kp) dtlp = edtl(kp) dtlpp = edtl(kp+1) ftimep = ftime(kp) ftimepp= ftime(kp+1) if (k .ne. 1) then km1 = k-1 do kk=km1,1,-1 if (edtl(kk) .gt. 0.0) then kp = kk vmxp = fvmx(kp) dtlp = edtl(kp) dtlpp = edtl(kp+1) ftimep = ftime(kp) ftimepp= ftime(kp+1) exit endif enddo end if !Find time since last landfall tland = (ftimek-ftimep) deldtl = dtlp - dtlpp if (dtlp .gt. 0.0 .and. deldtl .gt. 0.0) then tland = tland - (dtlp/deldtl)*(ftimepp-ftimep) endif !Adjust forecast intensity using the !Kaplan and DeMaria (JAM, 1995) exponential decay model fvmxt(k) = vb + (rf*vmxp - vb)*exp(-al*tland) endif !Check to see if intensity requires adjustment !for the case where the forecast track point is over land !but the ensemble track point is over water. if (fdtl(k) .lt. 0.0 .and. edtl(k) .gt. 0.0 & & .and. isltw .eq. 0) then !Find last forecast track point that was over water ! km1 = k-1 ! do kk=km1,0,-1 ! if (fdtl(kk) .gt. 0.0) then ! fvmxt(k) = fvmx(kk) ! exit ! endif ! enddo !! New Version That reduces search and searches forwards and backs pwk = -1 nwk = -1 km1 = k-1 do kk=k-1,max(k-9, 0),-1 if (fdtl(kk) .gt. 0.0) then pwk = kk fvmxt(k) = fvmx(kk) exit endif enddo do kk=k+1,min(k+3,nt) if (fdtl(kk) .gt. 0.0) then nwk = kk exit endif enddo if ((pwk .ge.0) .and. (nwk.ge.0)) then fvmxt(k) = fvmx(pwk) + (((fvmx(nwk) - fvmx(pwk)) / (nwk - pwk)) * (k-pwk)) elseif ((pwk .ge.0) .and. (nwk.lt.0)) then fvmxt(k) = fvmx(pwk) elseif ((pwk .lt. 0) .and. (nwk.ge.0)) then fvmxt(k) = fvmx(nwk) endif !! write(6,*) pwk, k, nwk, fvmxt(k) endif enddo ! end k loop do k=1,nt !Calculate the predicted value of !the intensity error if (k .eq. 1)then verrp = vcoef(1,k) + vcoef(2,k)*fvmx(k) + & & vcoef(3,k)*min(edtl(k),500.) elseif (k .gt. 1) then verrp = vcoef(1,k) & & + vcoef(2,k)*verro & & + vcoef(3,k)*fvmx(k) & & + vcoef(4,k)*min(edtl(k),500.) endif !Pick randomly from the intensity !error distributions, relative to the predicted values call randn(irand,ndr) irand = irand + 1 call sindex(irand,nrand,indx,nev(k)) verrp = verrp + verr(indx,k) evmx(m,k) = fvmxt(k) + verrp if (iicap .eq. 1) then !Cap the max wind of each realization as a function of distance !inland edtlt = edtl(k) if (edtlt .gt. 0.0) edtlt = 0.0 vcap = 20.0 + 120.0*exp(0.0035*edtlt) !write(6,997) m,k,edtl(k),evmx(m,k),vcap !997 format('m,k,dtl,v,vcap ',i4,1x,i4, !+ 1x,f8.0,1x,f7.1,1x,f7.1) if(edtl(k) .lt. dldiss .and. evmx(m,k) .gt. vcap) then evold = evmx(m,k) veold = verrp evmx(m,k) = vcap verrp = vcap - fvmxt(k) evnew = evmx(m,k) venew = verrp ! write(6,998) m,k,edtl(k),evold,evnew,veold,venew,vcap ! 998 format('m,k,dtl,eo,en,vo,vn,vcap: ',i4,1x,i4,1x,6(f8.2,1x)) endif ! /TODO: implment this edit for v1.6.0 ! I _think_ I came up with a better zombie check here. ! Instead of checking for any overland point below a threshold and ! zeroing everything subsequent. While generating the intensities, ! check if the current point is inland and if the previous vmx was ! below the threshold then zero it out. Any point that remains overland ! will continue to get zero'd, but a point that re-emerges will be able ! to regain the forecast intensity with pertubation. ! This was relevant during AL03 2020-06-03T12:00 ! if ( k.gt.1 .and. edtl(k) .lt. dldiss .and. evmx(m,k-1) .lt. vldiss) then ! evmx(m,k) = 0.0 ! end if endif !Save predicted error for next forecast time verro = verrp end do if (iicap .eq. 1) then !Check to see if an ensemble storm has dissipated over land. !If so, set the intensity for the rest of the track to !zero to prevent unrealistic regeneration over land from !the random intensity errors. do k=1,nt if (edtl(k) .lt. dldiss .and. evmx(m,k) .lt. vldiss) then do kk=k,nt evmx(m,kk) = 0.0 enddo exit endif enddo endif end do return !Error processing 900 continue write(lulog,950) fnidi 950 format(/,' Error opening intensity distribution file: ',a200,& & /,' in routine iens') stop 905 continue write(lulog,955) fnidi 955 format(/,' Error reading intensity distribution file: ',a200,& & /,' in routine iens') stop end subroutine iens subroutine tens(flat,flon,gpce,igpcef,indxg,demax,& & cal,sal,elat,elon,mte,ntr,mt,nt,lulog,ibasin) !This routine calculates the track member ensembles by !randomly sampling from the official track error distributions !divided into along and cross track components. The sampling !account for error correlations from one time period to the !next. implicit NONE real, dimension(0:nt) :: flat,flon,gpce,cal,sal real, dimension(mte,0:nt) :: elat,elon real :: demax integer :: igpcef, indxg(mt) integer :: mte,ntr,mt,nt,lulog,ibasin !Local variables integer, parameter :: mc=6000,mp=6 real, dimension(mp,mt) :: acoef,ccoef integer :: ng character (len=200) :: include_path, fntdi !Local variables for GPCE implementation integer, parameter :: ngcat=3 character (len=200) :: fntdg(ngcat), coef_fmt1, coef_fmt2 character (len=200) :: gpce_dbg_fmt, exp_fmt, flt_fmt integer, dimension(mt,0:ngcat) :: neag,necg,ncag,nccg real, dimension(mc,mt,0:ngcat) :: aerrg,cerrg real, dimension(mp,mt,0:ngcat) :: acoefg,ccoefg real :: gthresh(mt,ngcat,2) real :: rlon1, rlon2, rlat2, rlat1, rerr, dx, dy, aerro real :: aerrp, cerro, cerrp, daerr, dcerr, demin integer :: nct, n, ndr, net, nrand, m, lutdi, k, iskip, irand integer :: indx, igdist, i character (len=2) :: cbasin(4) !Initialize the coeficient arrays so they are zero if not filled ccoef(:,:) = 0 acoef(:,:) = 0 ccoefg(:,:,:) = 0 acoefg(:,:,:) = 0 !Set iskip=1 to skip perturbations of track. !For this case, the ensemble tracks are set equal !to the forecast track. iskip = 0 if (iskip .eq. 1) then !Skip perturbation of track do i=1,ntr do k=0,nt elat(i,k) = flat(k) elon(i,k) = flon(k) enddo enddo return endif !Specify name of file with track error distributions igdist=0 !Check to make sure GPCE distribution files are available if (ibasin .ge. 1 .and. ibasin.le.4) igdist = 1 if (igpcef .eq. 1 .and. igdist .eq. 0) then write(lulog, '(A)') 'GPCE tracks dist. not available in routine tens' stop end if include_path = get_parm_dir() cbasin = (/'al', 'ep', 'wp', 'sh'/) coef_fmt1 = '(A, A9,A2,A1,A4,A4)' write(fntdi, coef_fmt1) trim(include_path),"/mctdist_",cbasin(ibasin),'_',coeff_yrs,'.dat' coef_fmt2 = '(A, A9,A2,A1,A4,A1,i2.2,A4)' do n=1,ngcat write(fntdg(n), coef_fmt2) & & trim(include_path),"/mctdist_",cbasin(ibasin),'_',coeff_yrs,'_',n,'.dat' end do !Initialize random number routine and !specify number of places for random number !variable ndr=4 nrand= 10**ndr !Open and read the basin-wide track error file. !Note that the basin-wide info is read into the !the 0 element of various arrays. lutdi = 40 open(file=trim(fntdi),unit=lutdi,form='formatted',status='old',& & err=900) exp_fmt = '(10(e12.5,1x))' flt_fmt = '(10(f7.1,1x))' !Along track do k=1,mt read(lutdi,*,err=905,end=905) read(lutdi,*,err=905,end=905) neag(k,0),ncag(k,0) nct = ncag(k,0)+1 net = neag(k,0) read(lutdi,exp_fmt,end=905,err=905) (acoefg(n,k,0),n=1,nct) read(lutdi,flt_fmt,end=905,err=905) (aerrg(i,k,0),i=1,net) end do !Cross track do k=1,mt read(lutdi,*,err=905,end=905) read(lutdi,*,err=905,end=905) necg(k,0),nccg(k,0) nct = nccg(k,0)+1 net = necg(k,0) read(lutdi,exp_fmt,end=905,err=905) (ccoefg(n,k,0),n=1,nct) read(lutdi,flt_fmt,end=905,err=905) (cerrg(i,k,0),i=1,net) end do close(lutdi) if (igpcef .eq. 1) then !Open and read the track error file stratified by GPCE lutdi = 40 do m=1,ngcat open(file=fntdg(m),unit=lutdi,form='formatted',status='old',& & err=930) !Along track do k=1,mt read(lutdi,'(24x,f6.1,6x,f6.1)',err=935,end=935) & & gthresh(k,m,1), gthresh(k,m,2) read(lutdi,*,err=935,end=935) neag(k,m),ncag(k,m) nct = ncag(k,m)+1 net = neag(k,m) read(lutdi,exp_fmt,end=935,err=935) (acoefg(n,k,m),n=1,nct) read(lutdi,flt_fmt,end=935,err=935) (aerrg(i,k,m),i=1,net) enddo !Cross track do k=1,mt read(lutdi,*,err=935,end=935) read(lutdi,*,err=935,end=935) necg(k,m),nccg(k,m) nct = nccg(k,m)+1 net = necg(k,m) read(lutdi,exp_fmt,end=935,err=935) (ccoefg(n,k,m),n=1,nct) read(lutdi,flt_fmt,end=935,err=935) (cerrg(i,k,m),i=1,net) enddo close(lutdi) enddo endif !Calculcate track ensemble members !Set initial lat/lon equal to the forecast lat/lon do m=1,ntr elat(m,0) = flat(0) elon(m,0) = flon(0) enddo !Determine the index that determines !which GPCE distribution to use at each forecast !period to use. do k=1,nt if (igpcef .eq. 0 .or. gpce(k) .le. 0.0) then indxg(k) = 0 else indxg(k) = 0 do n=1,ngcat if (gpce(k) .gt. gthresh(k,n,1) .and.& & gpce(k) .le. gthresh(k,n,2) ) indxg(k) = n enddo endif gpce_dbg_fmt = '(A, i2,1x,i1,1x,f6.1,A,4(f6.1,1x))' write(6,gpce_dbg_fmt) 'k,indxg,gpce: ', k,indxg(k),gpce(k), & & ' gpce thresh: ', gthresh(k,1,1),gthresh(k,1,2), & & gthresh(k,3,1),gthresh(k,3,2) enddo do m=1,ntr do k=1,nt ng = indxg(k) !Calculate the predicted value of !the along track and cross track error !based upon the bias (t=12 h), or !the bias plus the serial correlation !(t > 12 hr) if (k .eq. 1) then aerrp = acoefg(1,k,ng) cerrp = ccoefg(1,k,ng) elseif (k .gt. 1) then aerrp = acoefg(1,k,ng) + acoefg(2,k,ng)*aerro cerrp = ccoefg(1,k,ng) + ccoefg(2,k,ng)*cerro endif !Pick randomly from the along/cross track !error distributions, relative to the predicted values call randn(irand,ndr) irand = irand + 1 call sindex(irand,nrand,indx,neag(k,ng)) aerrp = aerrp + aerrg(indx,k,ng) call randn(irand,ndr) irand = irand + 1 call sindex(irand,nrand,indx,necg(k,ng)) cerrp = cerrp + cerrg(indx,k,ng) if (k .gt. 1) then !Check to make sure along and cross track error change !does not exceed the upper limit demin = -demax daerr = aerrp - aerro dcerr = cerrp - cerro if (daerr .gt. demax) aerrp = aerro + demax if (daerr .lt. demin) aerrp = aerro + demin if (dcerr .gt. demax) cerrp = cerro + demax if (dcerr .lt. demin) cerrp = cerro + demin ! if (daerr .gt. demax .or. daerr .lt. demin .or.& ! & dcerr .gt. demax .or. dcerr .lt. demin ) then ! write(6,'(A,3(e11.4,1x),A,i3,A,i5)') & ! & 'AE or CE restricted, demax,daerr,dcerr: ',& ! & demax,daerr,dcerr,k,m, ' k=',i3,' niter=',i5 ! endif endif !Rotate along and cross track errors !into geographic coordinate system dx = cal(k)*aerrp + sal(k)*cerrp dy = sal(k)*aerrp - cal(k)*cerrp rlat1 = flat(k) rlon1 = flon(k) call distki(rlon1,rlat1,dx,dy,rlon2,rlat2) elat(m,k) = rlat2 elon(m,k) = rlon2 rerr = sqrt(dx*dx+dy*dy) !Save predicted error for next forecast time aerro = aerrp cerro = cerrp end do ! end k loop end do ! end m loop return !Error processing 900 continue write(lulog,950) fntdi 950 format(/,' Error opening track distribution file: ',a275,& & /,' in routine tens') stop 905 continue write(lulog,955) fntdi 955 format(/,' Error reading track distribution file: ',a275,& & /,' in routine tens') stop 930 continue write(lulog,980) fntdg(m) 980 format(/,' Error opening track distribution file: ',a275,& & /,' in routine tens') stop 935 continue write(lulog,985) fntdg(m) 985 format(/,' Error reading track distribution file: ',a275,& & /,' in routine tens') stop end subroutine tens subroutine rens(flat,fspd,fhead,fvmx,r34,r50,r64,r100,& & exdev,mte,ntr,mt, nt,lulog,ibasin) !This routine creates the wind radii ensembles. The size !of the cyclone ensemble is based on perturbation to the !x parameter. The routine returns 34, 50, 64 !and 100 kt wind radii based on track and intensity !perturbations implicit NONE !Passed arrays integer :: mte, ntr, mt, nt, lulog, ibasin real :: flat(0:nt),fspd(0:nt),fhead(0:nt),fvmx(0:nt) real :: r34(4),r50(4),r64(4),r100(4) real :: exdev(mte,0:nt) !Local variables integer, parameter :: mc=6000,mp=6 real :: dxerr(mc),dxcoef(mp), errmin real :: c34(4),c50(4),c64(4),c100(4) real :: p34(4),p50(4),p64(4),p100(4) real :: asym, beta, c, dxinit, dxnow, dxprev,rlat, rmax, vmax, x,xc integer :: i, ind, indx, irand, ityperp, k, lurdi, n integer :: nct, ndr, net, nrand character (len=200) :: include_path, fnrdi !initialize the coeficients so they are filled with zeros dxcoef(:)=0.0 !Set flag for type of wind perturbations !ityperp=0 for radii from climatology model !ityperp=1 for radii from climatology plus persistence model !ityperp=2 for radii from climatology plus persistence model !with random perturbations added ityperp = 2 !Initialize random number routine and !specify number of places for random number !variable ndr=4 nrand= 10**ndr ! Open file with climo wind model coefficients and xdev distributions include_path = get_parm_dir() ind = index( include_path, " " ) - 1 if (ibasin .eq. 1) then ! Atlantic Basin fnrdi = trim(include_path)//"/mcrdist_al_8803.dat" endif if (ibasin .eq. 2) then ! eastern North Pacific Basin fnrdi = trim(include_path)//"/mcrdist_ep_0103.dat" endif if (ibasin .eq. 3) then ! western North Pacific and N. Indian Basins fnrdi = trim(include_path)//"/mcrdist_wp_0102.dat" endif if (ibasin .eq. 4) then ! Southern Hemisphere Basin fnrdi = trim(include_path)//"/mcrdist_al_8803.dat" endif lurdi = 40 open(file=fnrdi,unit=lurdi,form='formatted',status='old',& & err=900) !Read climo wind model coefficients read(lurdi,*,err=905,end=905) nct if (nct .ne. nc) then write(lulog,'(A, i6, 1x, i6, A)') ' nc and nct inconsistent, nc,nct: ',& & nc, nct, ' in routine rens' stop end if !write(lulog,'(/, A)'), 'Climo wind model coefficients' do k=1,nct read(lurdi,'(f9.4,2x,a2)',err=905,end=905) coef(k),lcoef(k) ! write(lulog,190) coef(k),lcoef(k) enddo !Read dx error distribution read(lurdi,*,err=905,end=905) read(lurdi,*,err=905,end=905) net,nct nct = nct + 1 read(lurdi,'(10(E12.5,1x))',end=905,err=905) (dxcoef(n),n=1,nct) read(lurdi,'(10(f7.4,1x))',end=905,err=905) (dxerr(i),i=1,net) close(lurdi) !Fit initial radii to wind model vmax = fvmx(0) rlat = flat(0) c = fspd(0) beta = (90.0-fhead(0)) if (beta .lt. 0.) beta = beta + 360.0 !write(lulog,510) vmax,rlat,c,beta !510 format(/,'Input to initial vortex fit',/, !+ 'vmax=',f6.1,' rlat=',f6.1,' c=',f5.1,' beta=',f6.1) call xfit(vmax,rlat,c,beta,asym,rmax,xc,x,& & r34,r50,r64,r100,& & c34,c50,c64,c100,& & p34,p50,p64,p100,& & errmin) !write(lulog,520) rmax,asym,xc,x !520 format(/,'Output from initial vortex fit',/, !+ 'rmax=',f6.1,' asym=',f6.1,' xc=',f7.4,' x=',f7.4) !write(lulog,530) !530 format(/,'Wind radii from initial vortex fit') !write(lulog,531) (r34(k),k=1,4),(r50(k),k=1,4),(r64(k),k=1,4) !+ ,(r100(k),k=1,4) !531 format('r34,50,64, 100: ',4(4(f5.0,1x),2x)) !write(lulog,532) (c34(k),k=1,4),(c50(k),k=1,4),(c64(k),k=1,4) !+ ,(c100(k),k=1,4) !532 format('c34,50,64, 100: ',4(4(f5.0,1x),2x)) !write(lulog,533) (p34(k),k=1,4),(p50(k),k=1,4),(p64(k),k=1,4) !+ ,(p100(k),k=1,4) !533 format('p34,50,64, 100: ',4(4(f5.0,1x),2x)) if (ityperp .eq. 0) then !Set all the x perturbations to zero. x will then !be calculated from the climatology model in routine wpcal exdev(:,:) = 0.0 return endif !Set the t=0 x deviation of each ensemble to the difference !between the climo value and the fit value dxinit = x-xc do i=1,ntr exdev(i,0) = dxinit enddo !Calculate the x deviations at later times do i=1,ntr do k=1,nt dxprev = exdev(i,k-1) dxnow = dxcoef(1) + dxprev*dxcoef(2) if (ityperp .eq. 2) then !Add random component to x-deviation for t>0 call randn(irand,ndr) irand = irand + 1 call sindex(irand,nrand,indx,net) dxnow = dxnow + dxerr(indx) endif exdev(i,k) = dxnow enddo enddo return !Error processing 900 continue write(lulog,950) fnrdi 950 format(/,' Error opening radii distribution file: ',a200,& & /,' in routine rens') stop 905 continue write(lulog,955) fnrdi 955 format(/,' Error reading radii distribution file: ',a200,& & /,' in routine rens') stop end subroutine rens subroutine gdland_table(lon_in, lat_in, dtl) ! VERSION 1.0.1 ! Last updated: 02 Mar 2020 ! This program calculates the distance to land (dtl in km) ! Input: ! rlon,rlat - (deg W neg, and deg N pos) ! Output: ! dtl - Distance to nearest landmass in km ! Created 10 Feb 2020 (v1.0.0) (SS) ! ++ Modified from gfland.f - WSP version using lookup table for dtl (not fractional land) ! Updated 02 Mar 2020 (v1.0.1) (KM, AB) ! ++ IMPLICIT NONE and freeform format ! ++ adjusted ioper to include versions for SHIPS and WSP (uncomment needed version) ! ++ dtl returns missing (9999.) in case of error instead of 1000. ! ++ removed checks in subroutine for lsdiag length; moved to shipst7 IMPLICIT NONE ! Calling variables real, intent(in) :: lon_in, lat_in real, intent(out) :: dtl ! Local variables real, parameter :: rounding_factor = 100. ! round to 2 decimal places to reduce floating point diffs integer, parameter :: mx = 3601, my = 1801 integer, parameter :: iperiodic = 1 real, dimension(mx,my) :: gdtab character(len=200) :: include_path, fntab character(len=10), parameter :: gdtab_fmt='(10(f9.1))' integer :: iread !has table file been read integer :: ierr = 1 !initialize ierr to 1 integer :: lutab !file number integer :: i, j, ii, jj !counters ! Read in from table file header real :: slat, elat, dlat integer :: nlat real :: slon, elon, dlon integer :: nlon ! Temporary calculations real :: x, xlon, y, ylat, xy, rlon, rlat real :: w00, w01, w10, w11 ! Save table to only read file once during program execution save gdtab save slat,elat,dlat,nlat,slon,elon,dlon,nlon,iread data iread /1/ include_path = get_parm_dir() ! Set defaults dtl = 9999.!-999. if (iread .eq. 1) then ! Read the data tables lutab = 27 fntab = trim(include_path)//"/gdland_table.dat" open(file=trim(fntab),unit=lutab,form='formatted', & status='old',err=927) read(lutab,*,err=927,end=927) slon,elon,dlon,nlon, & slat,elat,dlat,nlat do j=1,nlat read(lutab,gdtab_fmt,err=927,end=927) (gdtab(i,j),i=1,nlon) enddo ! Add a column at lon=360 matching lon=0 if (iperiodic .eq. 1) then nlon=nlon+1 elon=elon+dlon do j=1,nlat gdtab(nlon,j)=gdtab(1,j) enddo endif iread=0 endif ! normalize input longitudes to 0-360 regardless of input rlon = mod(lon_in+360., 360.) rlat = lat_in ! Check input values if (rlat .lt. slat .or. rlat .gt. elat) then ierr=2 endif if (rlon .lt. slon .or. rlon .gt. elon) then ierr=3 endif if (ierr .ne. 1) then write(6,*) 'Error in subroutine gdland_table, ierr=',ierr write(6,*) 'slon,elon,slat,elat,rlon,rlat: ', & slon,elon,slat,elat,rlon,rlat return endif ! Calculate indices of table value closest but to the lower left ! of the requested lat/lon ii = 1 + ifix( (rlon-slon)/dlon ) jj = 1 + ifix( (rlat-slat)/dlat ) if (ii .eq. nlon) ii = nlon-1 if (jj .eq. nlat) jj = nlat-1 ! Calculate normalized x,y values relative to lower left table point xlon = slon + dlon*float(ii-1) x = (rlon-xlon)/dlon ylat = slat + dlat*float(jj-1) y = (rlat-ylat)/dlat xy = x*y ! Interpolate table values to find dtl and afraction w00 = 1.0+xy-(x+y) w10 = x-xy w01 = y-xy w11 = xy dtl = w00*gdtab(ii ,jj ) + & w10*gdtab(ii+1,jj ) + & w01*gdtab(ii ,jj+1) + & w11*gdtab(ii+1,jj+1) ! round result to reduce floating point differences dtl = real(nint(dtl*rounding_factor))/rounding_factor return 927 continue write(6,*) 'Error in subroutine gdland_table, ierr=',ierr write(6,*) 'while reading from file', trim(fntab), & slon,elon,slat,elat,rlon,rlat return end subroutine gdland_table subroutine calc_radii_bc(radarrt,mtetoa,mtitoa, mti,mt, ntr, & & flon, flat, eradt, ftimeh, rbc ) !Calculate mean of realization wind radii (from radii-CLIPER above) real, intent(in) :: radarrt(mtetoa,0:mtitoa,2,4,4) real, intent(in) :: eradt(4, mti, 4), flon(0:mt), flat(0:mt), ftimeh(0:mt) real, intent(out) :: rbc(4, mti) integer, intent(in) ::mtetoa,mtitoa, mti, mt, ntr real, dimension(4, 4) :: radmn, radmn_land real :: smooth_rbc(4,mti) real, dimension(4) :: radct, radct_land integer :: j, k,r, decayhr, iquad, kp1, kn1, land_ocean_flag real :: oradm(4, 4), fillvalue, min_rbc, max_rbc, trbc(4,4),trbc_ct(4) decayhr= 72 fillvalue = 1 min_rbc = 0.25 max_rbc = 5 land_ocean_flag = 1 ! always compare against ocean radii write(6, *) 'Starting Ensemble Bias Calculation' rbc(:,:) = 1 radmn = 0.0 radmn_land = 0.0 radct = 0.0 radct_land = 0.0 ! loop through forecast hours and calculate "correction" per quadrant ! save corrections at forecast hours do j = 1, mt-1 oradm = eradt(:,j,:) * cmta if (maxval(oradm) .eq. 0) then continue endif k = int(ftimeh(j-1)) ! write(6,*) 'k:', k, radarrt(1, k, i, 1, :) radmn = sum(radarrt(:, k, land_ocean_flag, :, :), 1) radct = count(sum(radarrt(:, k, land_ocean_flag, :, :), 3).gt.0, 1) * 1. do r=1,4 ! write(6,*) radmn(:,r), "|", radct where(radct .gt. 0) radmn(:,r) = radmn(:,r) / radct end do ! Calculate ratio between obs and forecast ! for each threshold and radial quadrant trbc = 0 where(oradm .gt. 0 .and. radmn .gt. 0) trbc = oradm/radmn trbc_ct = count(oradm .gt. 0 .and. radmn .gt. 0, 1) ! then average the correction across the tresholds, assuming that ! the windfield would increase somewhat consistently. where(trbc_ct .ge. 1) rbc(:,k+1) = sum(trbc,1)/trbc_ct end do ! apply min/max limits on correction where(rbc(:,:) .gt. max_rbc) rbc(:,:) = max_rbc where(rbc(:,:) .lt. min_rbc) rbc(:,:) = min_rbc ! smooth correction in time to reduce any spikes. ! utilises fill value of 1 if no correction was calculated. smooth_rbc = rbc do j=1, mt-2 kp1 = int(ftimeh(j+1)) + 1 k = int(ftimeh(j)) + 1 kn1 = int(ftimeh(j-1)) + 1 smooth_rbc(:,k) = (rbc(:,k) + rbc(:,kp1) + rbc(:,kn1))/3. end do rbc = smooth_rbc ! linear fill from forecast hours to hourly and extrapolate at end do iquad=1,4 call linfill(rbc(iquad,:), mti, fillvalue, fillvalue, decayhr) end do ! do j = 0, mt-1 ! k = int(ftimeh(j)) +1 ! write(6,*) 'rbc:', k, rbc(:,k+1) ! enddo end subroutine subroutine calc_ens_radii(mtetoa,mtitoa,nts, nte,mti,mt,i,ithresh, & & elats,etimes,elons,espds, eheads,evmxs,exdevs, & & r34i,r50i,r64i,r100i, vmxth, & & allradarr, rmwpt, apt,betatt,thpt, & & xtt, spdt) integer :: mtetoa,mtitoa,nts, nte real :: elats(0:mti),elons(0:mti),evmxs(0:mti),exdevs(0:mti) integer :: mti, mt, ithresh real :: r34i(0:mt,4),r50i(0:mt,4),r64i(0:mt,4),r100i(0:mt,4) real :: espds(0:mti),eheads(0:mti),etimes(0:mti), vmxth real :: rfix real :: allradarr(mtetoa,0:mtitoa,2,4,4) real :: rmwpt(mtetoa,0:mtitoa), apt(mtetoa,0:mtitoa) real :: betatt(mtetoa,0:mtitoa),thpt(mtetoa,0:mtitoa) real :: xtt(mtetoa,0:mtitoa), spdt(mtetoa,0:mtitoa) real :: Eo34(4),Eo50(4),Eo64(4),Eo100(4) real :: Eo34_land(4),Eo50_land(4),Eo64_land(4),Eo100_land(4) real :: Eo34t(4),Eo50t(4),Eo64t(4),Eo100t(4) real :: Eo34t_land(4),Eo50t_land(4),Eo64t_land(4),Eo100t_land(4) logical :: exist34,exist50,exist64,exist100 real :: r34t(4),r50t(4),r64t(4),r100t(4) real :: r34t_land(4),r50t_land(4),r64t_land(4),r100t_land(4) real, dimension(4) :: radarr integer kk, mm, i, iaper, ibias,ifixr, irtype, n real :: ap, thp, xt, asymt, betat, biasm, biasp real :: epsil, head real :: avg100p, avg34p, avg50p, avg64p real :: avg100p_land, avg34p_land, avg50p_land, avg64p_land real :: r34fix, r50fix, r64fix,r100fix,xtmn,xtmx real :: r34mx,r50mx,r64mx,r100mx,r34mn,r50mn,r64mn,r100mn real :: avg34i, avg50i, avg64i, avg100i real :: rlat1, rlon1, rmaxt, rmwp real :: sba64, spd, th0t, vmx1, wrmax,wt !++ Set the flag irtype to determine which method to use for !radii calculation. !irtype=0 for original climo vortex model, !irtype=1 for new climo vortex model plus perturbations. irtype=1 !++ Set ifixr=1 to fix the wind radii to constant values. !This option is normally used only for testing purposes. !Also set the fixed values (nmi). This option requires !that irtype=0. ifixr=0 r34fix = 125.0 r50fix = 75.0 r64fix = 50.0 r100fix= 25.0 if (ithresh .eq. 1) then rfix = r34fix elseif (ithresh .eq. 2) then rfix = r50fix elseif (ithresh .eq. 3) then rfix = r64fix elseif (ithresh .eq. 4) then rfix = r100fix else rfix = 0.0 endif !++ Specify minimum and maximum allowable 34,50,64 wind radii (nmi) !for radii for use when irtype=1. !(mx currently set to 99th percentile of observed Atlantic wind radii) ! r34mx = 500. ! Increased after Hurricane Sandy, formerly 350. (AS) r34mx = 350. r50mx = 200. r64mx = 130. r100mx = 95. r34mn = 15. r50mn = 10. r64mn = 8. r100mn = 8. !++ Specify minimum and maximum allowable x factor for Rankine vortex xtmn = 0.2 xtmx = 1.2 !Set ibias=1 for bias correction to r34, r64 and r100. This !adjustment corrects for systematic biases in the !Rankine vortex wind model used to estimate the radii. !This bias correction is applied to r34 (to increase r34) and !r64 (decrease r64) when the storm intensity is greater !than 94 kt. A constant 6% bias reduction is also applied !to r64 and r100 at all wind speeds. ibias=1 !Set iaper=1 to use the information in the initial observed !wind radii asymetries. iaper=1 !Determine if observed radii exist in any quadrant if (iaper .eq. 1) then avg34i=0. avg50i=0. avg64i=0. avg100i=0. do n=1,4 avg34i =avg34i + r34i(0,n)/4. avg50i =avg50i + r50i(0,n)/4. avg64i =avg64i + r64i(0,n)/4. avg100i=avg100i + r100i(0,n)/4. enddo exist34 = .false. exist50 = .false. exist64 = .false. exist100 = .false. if (avg34i .gt. 1.) exist34 = .true. if (avg50i .gt. 1.) exist50 = .true. if (avg64i .gt. 1.) exist64 = .true. if (avg100i .gt. 1.) exist100= .true. endif do kk=nts,nte rlat1 = elats(kk) rlon1 = elons(kk) vmx1 = evmxs(kk) vmx1 = max(vmx1, 5.0) spd = espds(kk) head = eheads(kk) betat = 90.0 - eheads(kk) if (betat .lt. 0.0) betat = betat + 360.0 if (irtype .eq. 0) then !Calculate radii from climatology vortex model call radcal1(rlat1,rlon1,vmx1,spd,head,vmxth,radarr, & & ap,rmaxt,thp,xt) if (ifixr .eq. 1) then !Set radii in each quadrant to a constant do mm=1,4 radarr(mm) = rfix enddo endif allradarr(i,kk,1,ithresh,:) = radarr allradarr(i,kk,2,ithresh,:) = radarr else !Calculate vortex model parameters call vpcal(vmx1,rlat1,spd,asymt,th0t,rmaxt,xt) !Add perturbation to x xt = xt + exdevs(kk) !Check bounds of perturbed x if (xt .gt. xtmx) xt = xtmx if (xt .lt. xtmn) xt = xtmn !Calculate wind radii from vortex model call radcal(rlat1,vmx1,asymt,betat,th0t,rmaxt,xt, & & r34t,r50t,r64t,r100t,ap,rmwp,thp) !Calculate wind radii from vortex model call radcal_land(rlat1,vmx1,asymt,betat,th0t,rmaxt,xt, & & r34t_land,r50t_land,r64t_land,r100t_land) if (ibias .eq. 1) then !Apply bias correction to vortex model radii !Calculate speed dependent vortex model bias factor if (vmx1 .ge. 95.) then epsil = 0.0025*(vmx1-94.) else epsil = 0.0 endif biasp = 1.0 + epsil biasm = 1.0 - epsil !Specify systematic r64 and r100 bias adjustment sba64 = 0.94 do n=1,4 r34t(n) = r34t(n) *(1.0 + epsil) r64t(n) = r64t(n) *(sba64 - epsil) r100t(n) = r100t(n)*(sba64 - epsil) r34t_land(n) = r34t_land(n) *(1.0 + epsil) r64t_land(n) = r64t_land(n) *(sba64 - epsil) r100t_land(n) = r100t_land(n)*(sba64 - epsil) enddo endif if (iaper .eq. 1) then !Radii will exactly match observed radii at t=0. !The influence of the initial radii decays exponentially !with time, with an e-folding time of 32 hours. !Calculate initial radii errors if (kk .eq. 0) then do n=1,4 Eo34(n) = r34i(0,n) - r34t(n) Eo50(n) = r50i(0,n) - r50t(n) Eo64(n) = r64i(0,n) - r64t(n) Eo100(n)= r100i(0,n)- r100t(n) Eo34_land(n) = r34i(0,n) - r34t_land(n) Eo50_land(n) = r50i(0,n) - r50t_land(n) Eo64_land(n) = r64i(0,n) - r64t_land(n) Eo100_land(n)= r100i(0,n)- r100t_land(n) end do endif !Average of predicted radii. avg34p=0. avg50p=0. avg64p=0. avg100p=0. do n=1,4 avg34p = avg34p + r34t(n) / 4. avg50p = avg50p + r50t(n) / 4. avg64p = avg64p + r64t(n) / 4. avg100p = avg100p + r100t(n) / 4. enddo !Average of predicted radii. avg34p_land=0. avg50p_land=0. avg64p_land=0. avg100p_land=0. do n=1,4 avg34p_land = avg34p_land + r34t_land(n) / 4. avg50p_land = avg50p_land + r50t_land(n) / 4. avg64p_land = avg64p_land + r64t_land(n) / 4. avg100p_land = avg100p_land + r100t_land(n) / 4. enddo !Make sure avg predicted radii are not zero. avg34p = max(avg34p, 1.0) avg50p = max(avg50p, 1.0) avg64p = max(avg64p, 1.0) avg100p = max(avg100p,1.0) avg34p_land = max(avg34p_land, 1.0) avg50p_land = max(avg50p_land, 1.0) avg64p_land = max(avg64p_land, 1.0) avg100p_land = max(avg100p_land,1.0) do n=1,4 Eo34t(n)=0. Eo50t(n)=0. Eo64t(n)=0. Eo100t(n)=0. Eo34t_land(n)=0. Eo50t_land(n)=0. Eo64t_land(n)=0. Eo100t_land(n)=0. end do !Calculate time weight wt=exp(-1.*etimes(kk)/32.) !Calculate time weighted average of t=0 observed !and current predicted radii do mm = 1,4 if (exist34) then r34t(mm) = r34t(mm)+wt*Eo34(mm) r34t_land(mm) = r34t_land(mm)+wt*Eo34_land(mm) else r34t(mm) = r34t(mm) r34t_land(mm) = r34t_land(mm) endif if (exist50) then r50t(mm) = r50t(mm) + wt*Eo50(mm) r50t_land(mm) = r50t_land(mm)+wt*Eo50_land(mm) elseif (exist34) then Eo50t(mm)=avg50p/avg34p*Eo34(mm) r50t(mm) = r50t(mm)+wt*Eo50t(mm) Eo50t_land(mm)=avg50p/avg34p*Eo34_land(mm) r50t_land(mm) = r50t_land(mm)+wt*Eo50t_land(mm) else r50t(mm) = r50t(mm) r50t_land(mm) = r50t_land(mm) endif if (exist64) then r64t(mm) = r64t(mm)+wt*Eo64(mm) r64t_land(mm) = r64t_land(mm)+wt*Eo64_land(mm) elseif (exist50) then Eo64t(mm)=avg64p/avg50p*Eo50(mm) r64t(mm) = r64t(mm)+wt*Eo64t(mm) Eo64t_land(mm)=avg64p/avg50p*Eo50_land(mm) r64t_land(mm) = r64t_land(mm)+wt*Eo64t_land(mm) elseif (exist34) then Eo64t(mm)=avg64p/avg34p*Eo34(mm) r64t(mm) = r64t(mm)+wt*Eo64t(mm) Eo64t_land(mm)=avg64p/avg34p*Eo34_land(mm) r64t_land(mm) = r64t_land(mm)+wt*Eo64t_land(mm) else r64t(mm) = r64t(mm) r64t_land(mm) = r64t_land(mm) endif if (exist100) then r100t(mm) = r100t(mm)+wt*Eo100(mm) r100t_land(mm) = r100t_land(mm)+wt*Eo100_land(mm) elseif (exist64 ) then Eo100t(mm)=avg100p/avg64p*Eo64(mm) r100t(mm) = r100t(mm)+wt*Eo100t(mm) Eo100t_land(mm)=avg100p/avg64p*Eo64_land(mm) r100t_land(mm) = r100t_land(mm)+wt*Eo100t_land(mm) elseif (exist50 ) then Eo100t(mm)=avg100p/avg50p*Eo50(mm) r100t(mm) = r100t(mm)+wt*Eo100t(mm) Eo100t_land(mm)=avg100p/avg50p*Eo50_land(mm) r100t_land(mm) = r100t_land(mm)+wt*Eo100t_land(mm) elseif (exist34 ) then Eo100t(mm)=avg100p/avg34p*Eo34(mm) r100t(mm) = r100t(mm)+wt*Eo100t(mm) Eo100t_land(mm)=avg100p/avg34p*Eo34_land(mm) r100t_land(mm) = r100t_land(mm)+wt*Eo100t_land(mm) else r100t(mm) = r100t(mm) r100t_land(mm) = r100t_land(mm) endif enddo endif ! end iaper nest do mm=1,4 !Check for wind radii vs vmax if (vmx1 .lt. 34.) r34t(mm )=0. if (vmx1 .lt. 50.) r50t(mm )=0. if (vmx1 .lt. 64.) r64t(mm )=0. if (vmx1 .lt. 100.) r100t(mm)=0. if (vmx1 .lt. 34.) r34t_land(mm )=0. if (vmx1 .lt. 50.) r50t_land(mm )=0. if (vmx1 .lt. 64.) r64t_land(mm )=0. if (vmx1 .lt. 100.) r100t_land(mm)=0. !Restrict radii to minimum observed values if (r34t(mm) .lt. r34mn ) r34t(mm)=0. if (r50t(mm) .lt. r50mn ) r50t(mm)=0. if (r64t(mm) .lt. r64mn ) r64t(mm)=0. if (r100t(mm) .lt. r100mn) r100t(mm)=0. if (r34t_land(mm) .lt. r34mn ) r34t_land(mm)=0. if (r50t_land(mm) .lt. r50mn ) r50t_land(mm)=0. if (r64t_land(mm) .lt. r64mn ) r64t_land(mm)=0. if (r100t_land(mm) .lt. r100mn) r100t_land(mm)=0. !Restrict radii to maximum observed values if (r34t(mm) .gt. r34mx ) r34t(mm ) = r34mx if (r50t(mm) .gt. r50mx ) r50t(mm ) = r50mx if (r64t(mm) .gt. r64mx ) r64t(mm ) = r64mx if (r100t(mm) .gt. r100mx) r100t(mm) = r100mx if (r34t_land(mm) .gt. r34mx ) r34t_land(mm ) = r34mx if (r50t_land(mm) .gt. r50mx ) r50t_land(mm ) = r50mx if (r64t_land(mm) .gt. r64mx ) r64t_land(mm ) = r64mx if (r100t_land(mm) .gt. r100mx) r100t_land(mm) = r100mx end do !Check thresholds and asign values to those radii !Also, apply the max to avg radii correction factor allradarr(i,kk,1,1,:) = r34t * cmta allradarr(i,kk,1,2,:) = r50t * cmta allradarr(i,kk,1,3,:) = r64t * cmta allradarr(i,kk,1,4,:) = r100t * cmta allradarr(i,kk,2,1,:) = r34t_land * cmta allradarr(i,kk,2,2,:) = r50t_land * cmta allradarr(i,kk,2,3,:) = r64t_land * cmta allradarr(i,kk,2,4,:) = r100t_land * cmta endif !Find the smallest non-zero 34, 50, 64 kt radius wrmax = 999.9 do mm=1,4 if (r34t(mm) .gt. 0.0 .and. r34t(mm) .lt. wrmax) wrmax = r34t(mm) if (r50t(mm) .gt. 0.0 .and. r50t(mm) .lt. wrmax) wrmax = r50t(mm) if (r64t(mm) .gt. 0.0 .and. r64t(mm) .lt. wrmax) wrmax = r64t(mm) enddo !Apply max to mean in quadradant adjustment to wrmax if (wrmax .lt. 999.0) wrmax = cmta*wrmax !Check to make sure the rmw is smaller than the smallest !non-zero wind radii. If not, set it to a fraction of the smallest !radius. if (rmwp .ge. wrmax .and. wrmax .lt. 100.0 & & .and. wrmax .gt. 0.0) rmwp=rbuf*wrmax !Save rmwp for TOA calculation rmwpt(i,kk) = rmwp !Save other vortex parameters for TOA calculation xtt(i,kk) = xt apt(i,kk) = ap betatt(i,kk) = betat thpt(i,kk) = thp spdt(i,kk) = spd enddo end subroutine calc_ens_radii subroutine mcprob ( cbasin, flatt, flont, fvmxt, igpcef, gpcet,& & indxg,sname,isnum,iyr,imon,iday,itime,& & mt, ktime, ntr, ithresh, its, ite, kstime,& & ntsave, iprall,r34, r50, r64, r100,& & nxy, mxy, plon, plat, wprobs, wprobi, wprob,& & wtprob, ierr ) ! c**************************************************************************** !Version 8.0.1 ! c**************************************************************************** !This routine generates a field of probabilities !of winds exceeding a certain threshold (determined by ithresh) !from time its to time ite. !This is an experimental version where the error distributions depend !on the GPCE values. !INPUT: !cbasin - character (len=2) TC basin (AL,al,EP,ep,CP,cp,WP,wp,SH, !sh,IO,io, etc...) !flatt(0:mt)- forecast lats (deg N, 0 for missing) at 12 hr intervals !flont(0:mt)- forecast lons (deg E, 0 for missing) at 12 hr intervals !Note: longitude is continuous (0.01-360.00 E) !fvmxt(0:mt)- forecast intensity (kt, 0 for missing) at 12 hr intervals !igpcef - Flag for using error distributions stratified by GPCE !- values (=1 to include the option, =0 to turn it off) !gpcet(0:mt)- Expected consensus model forecast error (nmi) from the !GPCE algorithm. Not used if igpcef=0. !indxg(mt) - GPCE error index (1-3) !sname - character (len=10) storm name !isnum - storm number (01-99) !iyr - year (eg, 2010) !imon - month (1-12) !iday - day of month (1-31) !itime - hour of day (0-23) !ktime(0:mt)- array contains forecast times (hrs) !ntr - number of MC realizations (ensembles) !Note: must be .le. mte !ithresh - Wind speed threshold indicator !(1 for 34 kt, 2 for 50 kt, 3 for 64 kt, 4 for 100 kt) !its - Starting time (hr) for calculating probabilities !ite - Ending time (hr) for calculating probabilities !kstime(0:ntsave) - Array containing times for saving cumulative and !incremental probabilities. The times must be sequential !in hours and must be at least 2 hours apart from its to ite. !The first element of kstime (kstime(0)) is always set to its. !and the last element (kstime(ntsave) is always set to ite. !If some other value is provided in kstime( 0) it !is replaced by its. !If some other value is provided in kstime(ntsave) it !is replaced by ite. !For probabilities at evenly spaced 12 hour time intervals from !0 to 120 hours, set its=0 ite=120 ntsave=10 and !kstime(n) for n=0,10 to 0, 12, 24, 36, 48, 60, 72, 84, 96, 108, 120 !The values in kstime do not need to be evenly spaced. !For example, to define intervals consistent with the times !of the NHC official forecast points, set its=0 ite=120 ntsave=7 !and kstime(n) for n=0,7 to 0,12,24,36,48,72,96,120 !ntsave - The number of times to save the probabilities !from its to ite. !in the array wprob at the end of each interval. !iprall - Flag for printing a file with all the realizations in it !=0 to skip printing !=1 to print in simple ASCII format !=2 to print in ATCF format !r34(0:mt,4) - radius of 34 kt winds (nmi) to NE, SE, SW, NW at t=0,...,mt*12 !r50(0:mt,4) - radius of 50 kt winds (nmi) to NE, SE, SW, NW at t=0,...,mt*12 !r64(0:mt,4) - radius of 64 kt winds (nmi) to NE, SE, SW, NW at t=0,...,mt*12 !r100(0:mt,4) - radius of 100 kt winds (nmi) to NE, SE, SW, NW at t=0,...,mt*12 !Note: Set radii to zero if missing !nxy - Number of lat/lon points to calculate probabilities !mxy - Maximum number of lat/lon probability points !plon(mxy) - 1-D longitude array for calculation of probabilities !(deg west negative) !plat(mxy) - 1-D latitude array for calculation of probabilities !(deg north) !wprobs - 2-D array (mxy,0:ntsave) temporary work array !cbasin - character (len=2) TC basin (AL,EP,CP,WP,SH,IO etc...) !OUTPUT: !wprob - 2-D array (mxy,0:ntsave) containing cumulative wind !probabilities at plon,plat !from its to ite, saved every (ite-its)/ntsave hours !wprobi - 2-D array (mxy,0:ntsave) containing interval wind !probabilities at plon,plat !from its to ite, saved every (ite-its)/ntsave hours !wtprob(7,9)- Probability values (%) for the the NHC wind speed probability !table. !wtprob(k,i) k=1,7 for t=12,24,36,48,72,96,120 !i=1,9 for Dissipated, TD, TS, Hurricane, !cat1, cat2, cat3, cat4, cat5 !ierr - error flag (=0 if okay, >0 if error) !mcrall... .dat - If iprall=1 then an output file is created !containing the track, intensity and structure info from !all the realizations. The file name is of the form !mcrall_bbnnyyyy_mmddtt.dat where bb is the basin ID, !nn is the ATCF storm number, yyyy is the year, mm is !month, dd is the day and tt is the UTC time. !If iprall=2, the file name will end with _adeck.dat. !mctoar... .dat If iptoar=1 then an output file is created !containing the time of arrival probabilities. the file !name is of the form mctoar_bbnnyyyy_mmddtt.dat, similar !to the mcrall... .dat file !***************************************************************************** !Written by Mark DeMaria and John Knaff (NOAA/NESDIS), and Buck Sampson (NRL) !Last updated: February 16, 2004 !Modified August 4, 2004 -- wind radii smaller than the climatological !radius of maximum wind fix, wind radii !cliper bias correction fixed !Modified April 20, 2005 - 1. Radii adjustment for max in quadrant versus !mean in quadrant. !2. t=0 hour probability added. !3. Extra output array for total probabilities !in each sub-interval !Modified April 6, 2007 - 1. variable ibasin is now a function of !location instead of basin of orgin. This !effects which error distributions and !wind radii climatology is used throughout !the relization generation process (JAK). !Modified March 18, 2008 - 1. New calling argument added (wtprob) to !return info for the wind speed probability !table (MD). !2. Subroutine idfile and call from main routine !added for WSPT calculations. !Modified May 30, 2008 - 1. ioper variable added to automatically !select proper include path !2. Option for new inland wind "capping" !function to restrict realization max !winds to maximum observed for given !distance inland (see iicap variable !in routine iens) !Modified Nov 2008 - GPCE option added !Modified May 2010 - Option to print file with all realizations added !Modified Mar 2011 - Correction to azimuthal interpolation of wind radii !so outer radii cannot be < RMW !Modified Apr 2011 1. Fixes for southern hemisphere added for BoM applications !2. Program uses new track,intensity distributions for !2011 season (based on 2006-2010 forecasts) !3. demax variable added to prevent unrealistically !large changes in the along and cross track errors. !4. Modification to make 34, 50, 64 probs consistent. !Fixes the problem identified in the Fay (2008) forecasts. !Modified Jun 2011 1. Check added to make sure rmw is always less !that smallest nonzero wind radii in any of the !4 quadrants. Fixes the problem identified during !Beatriz (2011) EPAC forecasts. !Modified Jul 11 2011 1. Format of last item in mcrall changes from f5.1 to f6.1 !to prevent overflow (MD) !Modified Feb 2012 - Fixed bug where incremental probabilities were being set to !0 at the last time step when a full 5-day official forecast !was not available (A. Schumacher, CIRA) !Modified May 2012 - Fixed problem experienced with Atlantic TCs that cross !0 E - interpolation between longitudes on either !side of 0 E were incorrect (e.g., interpolating !between 359.0 E and 1.0 as 180.0 E). For these cases, !360 is added to the forecast longitudes that cross 0 E. !(A. Schumacher, CIRA) !Modified Mar 2013 - Updated error distribution files (2008-2012) and !filenames in code, v3.03 (A. Schumacher, CIRA) !Modified Jul 2014 - Updated error distribution files (2009-2013) and !filenames in code, made consistent with ATCF version !provided by B. Sampson, v4.0 (A. Schumacher, CIRA) !Modified Nov 2014 - Added options to use spline interpolation instead of !linear interpolation using the "switches": !isintp (ofcl forecast track) = 0 lin int, = 1 spline !isintpe (realization tracks) = 0 lin int, = 1 spline !Updated code to version 4.1 (A. Schumacher, CIRA) !Modified Jan 2015 - Added time of arrival/departure calculation, which !outputs new mctoar files. Set iptoar parameter below, and !related wprobmin. Version changed to 5.0.0. Three digit !convention used to ahdere to NCO standards. !Modified Feb 2015 - v.5.0.1: Increased character size of variables !that hold filenames to allow for longer data paths !(Sampson and Mattocks). !Modified Jan 2016 - Updated error distribution files (2010-2014) !(A. Schumacher) !Modified Apr 2016 - v6.0.0 - This version is for the 2016 season. !AL and EP track and intensity distributions updated !to include 2011-2015. WP is still 2010-2014. !Several bug fixes: !1. iyr2 calculation removed in subroutine idfile !2. avg34p, ... 100p values checked to make sure they are not zero in wpcal !3. Special calculation for Yokosuka commented out in wpcal !5. mxly corrected to mxyl in wpcal !6. fntdi corrected to fnidi in 955 error processing in subroutine iens !7. Call to splintt with "tension" variable instead of numeric value !8. NCONWOB dimension increased to 900 in subroutine WLAND !Modified Apr 2016 - v6.0.1 - Replaced tab characters with spaces, !changed line endings of all files !from DOS to unix. !Modified Jan 2018 - v7.0.0 - This version is for the 2018 season. !1. Updated error distributions, AL/EP 2013-2017 !(2017 provisional BT used), WP2012-2016 (Schumacher) !2. Added hourly interpolation for cases crossing land !via subroutine landcross (Schumacher) !Uses subroutine afland and static file atable.dat !Uses subroutine wafland and static file wtable.dat !3. Reduced wind radii (vmax fac=0.8) used for realizations !over land using hourly interp. (Schumacher) !Uses subroutine radcal_land !Modified Jul 2018 - v7.0.1 - Extended overland updates (Jan 2018) to the part of !the code that runs on breakpoints. Solves !inconsistency between gridded and bpt probs. !(Schumacher) !Modified Feb 2019 - v8.0.1 - Replaced basin-specific subs (afland/wafland) calculating !dist to land with global version: gfland !Requires static global land file, gtable.dat, in include !subdirectory. !Changed "tension" for spline subroutines to 0.9. !Updated coefficients to 2014-2018 errors using !preliminary best tracks.(Schumacher) !Modified Feb 2020 - v9.0.1 - Updated to run out to 7 days if requested !- Replaced a/w/shland calls with gdland_table !- Read forecast R34/50/64/100 and pass to wpcal !- Reset r34mx = 350 !- Fixed potential bug in spline routines !- convert to free format and implicit none ! ***************************************************************************** implicit None !** Passed arrays !Specify maximum number of forecast times integer :: mt integer :: ibasin character (len=2) :: cbasin ! character (len=8) :: atcfid character (len=10) :: sname real, dimension(0:mt) :: flatt,flont,fvmxt, gpcet integer, dimension(mt) :: indxg integer :: ktime(0:mt),kstime(0:ntsave) real r34(0:mt,4),r50(0:mt,4),r64(0:mt,4),r100(0:mt,4) real :: plat(mxy),plon(mxy) real :: wprobs(mxy,0:ntsave),wprob(mxy,0:ntsave) real :: wprobi(mxy,0:ntsave),wtprob(7,9) integer ithresh !** Local arrays real :: ftime(0:mt) integer :: igpcef,isnum,iyr, imon, iday, itime, its, ite, ierr integer :: idprt, iatcf, ivmxth, k, kk, lulog, iprall integer nxy, ntr, mxy,ntsave, nt real :: ftemp, alpha !Additional forecast arrays real :: fdtl(0:mt),fspd(0:mt),fhead(0:mt) real :: flon(0:mt),flat(0:mt),fvmx(0:mt),gpce(0:mt) real :: cal(0:mt),sal(0:mt) !Spline interpolation arrays integer :: mts real :: spl_lon(1:mt+1),spl_lat(1:mt+1), spl_tim(1:mt+1) real :: spl_lon2(1:mt+1),spl_lat2(1:mt+1) !Arrays for ensemble members (track, intensity radii) integer, parameter :: mte=5000 real :: elat(mte,0:mt),elon(mte,0:mt) real :: evmx(mte,0:mt) real :: exdev(mte,0:mt) !Arrays for climo wind radii model character (len=200) :: include_path, namelist_filepath integer :: iptoar, isintp, isintpe, irbc, ipgrid real :: demax, tension, wprobmin, tensione NAMELIST /config/ ioper, iptoar, demax, isintp, tension, isintpe, & & wprobmin, tensione, irbc, coeff_yrs, iprall, ipgrid !goptc ++ !Arrays for code optimization for regular lat/lon grid integer, parameter :: mxe=2000,mye=1000 real :: plone(mxe),plate(mye) real :: dlone, dlate integer nxe, nye, iego common /egopt/ dlone,dlate,plone,plate,nxe,nye,iego !goptc ++ !**** Executable code **** !These are now defaults !!! Don't edit them here, change the respective namelist file !Set ioper=0 for CIRA test runs !or ioper=1 for runs on the ATCF ioper=2 ! for runs on WCOSS !or ioper=0 !Set demax to the maximum allowable change (km) in the along and cross !track errors between 12 h intervals. This prevents the random sampling !technique to change these faster than in the original error data. demax = 600.0 !Filling in forecast tracks: set =0 for linear interp, =1 for spline isintp=1 !If isintp=1, set the tension for the spline routine. !tension should be between 0.0 and 1.0 tension = 0.9 !Interpolating ensemble tracks: set =0 for linear interp, =1 for spline isintpe=1 !If isintpe=1, set the tension for the spline routine used for !interpolating ensemble tracks !tensione should be between 0.0 and 1.0 tensione = 0.9 !Set iptoar=1 to calculate and print time of arrival probabilities !Other options are iptoar=11 to only calculate TOA if ithresh=1 !iptoar=12 to only calculate TOA if ithresh=2 !iptoar=13 to only calculate TOA if ithresh=3 !iptoar=14 to only calculate TOA if ithresh=4 !iptoar= 0 to skip TOA calculation iptoar=0 ipgrid=1 ! irbc=1 will adjust forecast radii to average OFCL, issues with asymetries remain irbc = 0 !For iptoar .ne. 0 specify min value of cumulative probability !to perform toa calculation wprobmin = 5.0 ierr = 0 lulog = 6 call getenv ( "WSP_NAMELIST", namelist_filepath ) ! WCOSS if (len_trim(namelist_filepath) .gt. 0) then write(lulog, '(A, A)') 'Reading namelist from:', trim(namelist_filepath) OPEN(FILE=trim(namelist_filepath), UNIT=57, STATUS='OLD', IOSTAT=ierr) if (ierr .ne. 0) then write(0,'(A, i4, A)') 'Error:',ierr, 'Opening '//trim(namelist_filepath) stop end if READ(57, nml=config, IOSTAT=ierr) close(57) if (ierr .gt. 0) then write(0,'(A, i4, A)') 'Error:',ierr, 'Reading '//trim(namelist_filepath) stop end if else write(lulog, '(A)') 'Running on default configuration' end if write(lulog,'(A5, i2)') 'ioper=', ioper write(lulog,'(A6, i2)') 'iptoar=', iptoar write(lulog, '(A7, A)') 'cbasin=', cbasin write(lulog, '(A7, i2)') 'irbc=', irbc write(lulog, '(A10, A4)') 'coeff_yrs=', coeff_yrs !Re-set iptoar for theshold-dependent cases if (iptoar .gt. 1) then if (iptoar .eq. 11 .and. ithresh .eq. 1) iptoar=1 if (iptoar .eq. 12 .and. ithresh .eq. 2) iptoar=1 if (iptoar .eq. 13 .and. ithresh .eq. 3) iptoar=1 if (iptoar .eq. 14 .and. ithresh .eq. 4) iptoar=1 endif if (iptoar .gt. 1) iptoar=0 !Check input if (ithresh .eq. 1) then ivmxth = 34 elseif (ithresh .eq. 2) then ivmxth = 50 elseif (ithresh .eq. 3) then ivmxth = 64 elseif (ithresh .eq. 4) then ivmxth = 100 else write(lulog,'(a, i6)') 'Illegal value of ithresh in mcprob: ', ithresh ierr = 1 return endif if (its .lt. 0) then write(lulog,'(a,i6)') ' Starting hour can not be negative: ', its ierr = 2 return endif if (ite .gt. 168) then write(lulog,'(a,i6)') ' Ending hour can not be > 168 h: ', ite ierr = 3 return endif if (ite .le. its) then write(lulog,'(a)') ' Time interval for wind probs must be positive' ierr = 4 return endif if (ntr .gt. mte) then write(lulog,'(a)') 'Requested ntr too large, increase mte parameter' ierr = 5 return endif if (nxy .gt. mxy) then write(lulog,'(a,i6,1x,i6)') & & 'Requested nxy too large, increase mxy parameter nxy,mxy: ,', & & nxy,mxy ierr = 6 return endif if (ntsave .lt. 1) then write(lulog,'(a)') ' ntsave must be at least 1' ierr = 7 return endif kstime( 0) = its kstime(ntsave) = ite do k=1,ntsave if (kstime(k) .le. kstime(k-1)) then write(lulog,*) ' kstime values must be sequential' ierr=8 return endif if (kstime(k) .gt. ite) then write(lulog,*) ' Warning: kstime value is greater than ite' endif if (kstime(k) .lt. its) then write(lulog,*) ' Warning: kstime value is less than its' endif enddo !goptc ++ !Check for a regularly spaced lat/lon grid to determine if grid !optimization feature can be turned on. call gopck(plon,plat,mxy,nxy) print *, 'iego=',iego !goptc ++ !Convert forecast times to real ftime = 0. do k=0,mt ftime(k) = float(ktime(k)) enddo !Move forecast lat,lon,vmx from input arrays flat = 0. flon = 0. fvmx = 0. flat(:mt) = flatt(:mt) flon(:mt) = flont(:mt) fvmx(:mt) = fvmxt(:mt) !2012 Fix - when Atlantic TC's cross 0 E, sets lon to 360 + lon !so interpolation is done correctly (ABS) do k=0,mt if (cbasin .eq. 'AL' .or. cbasin .eq. 'al') then if (flon(k) .gt. 0.0 .and. flon(k) .lt. 100.0) then flon(k) = 360.0 + flon(k) endif endif enddo !Fill in missing forecast track positions using !linear interpolation (ili = 0) or spline (ili = 1) IF (isintp .EQ. 0) THEN call linte(flat,fdtl,0,mt) call linte(flon,fdtl,0,mt) ELSE kk = 0 spl_tim = 0. spl_lat = 0. spl_lon = 0. DO k = 0, mt if (flat(k) .eq. 0. .and. flon(k) .eq. 0.) CYCLE kk = kk + 1 spl_tim(kk) = ftime(k) spl_lat(kk) = flat(k) spl_lon(kk) = flon(k) ENDDO mts = kk IF (mts .GT. 1) THEN CALL spline(spl_tim(1:mts),spl_lat(1:mts),mts, & & 99e31,99e31,spl_lat2(1:mts)) CALL spline(spl_tim(1:mts),spl_lon(1:mts),mts, & & 99e31,99e31,spl_lon2(1:mts)) DO k=1,mt IF (ftime(k) .LE. spl_tim(mts) .AND. ftime(k) .GT. 0.) THEN CALL splintt(spl_tim(1:mts),spl_lat(1:mts), & & spl_lat2(1:mts),mts,ftime(k),tension,flat(k)) CALL splintt(spl_tim(1:mts),spl_lon(1:mts), & & spl_lon2(1:mts),mts,ftime(k),tension,flon(k)) ENDIF ENDDO ENDIF ENDIF !Write over flat and flon with filled-in arrays so interpolated tracks !are passed by to mcdriver for writing output files flatt = flat flont = flon !Find the last time with a valid forecast track nt = mt do k=0,mt if (flon(k) .le. 0.0) then nt = k-1 exit endif enddo if (nt .lt. 1) then ierr=10 return endif !Move and interpolate the GPCE values if (igpcef .eq. 1) then do k=0,mt gpce(k) = gpcet(k) enddo call linte(gpce,fdtl,0,mt) endif !Initialize (idum < 0) the random number generator ran2 !once a more random value of idum could be utilized !here in the future. !itemp=abs(nint(flatt(0)/100.)) !itemp=-1*(iday+imon+ihour+itemp) !ftemp=ran2(itemp) ftemp=ran2(-1) !Switch the sign of storm longitude to deg W negative !convention and calcualte distance to land along !forecast track. The longitude convention when !calling this routine is a value 0 to 360. ibasin = 0 do k=0,mt call gdland_table(flon(k), flat(k), fdtl(k)) end do ! convert basin to integer select case(cbasin) case ('AL', 'al') ! Atlanti ibasin = 1 case ('EP' ,'ep', 'CP', 'cp') ! E./C. North Pacifi if (flon(0).ge.180.0) then !NHC/CPHC AOR ibasin = 2 !if a storm that forms in the cp or ep and tracks across ! the dateline (e.g., Ioke 2006) then use JTWC errors and ! the correct wind radii CLIPER (JAK) else if (flon(0).lt.180.0) then ! JTWC AOR ibasin = 3 end if case ('WP', 'wp', 'IO', 'io') ! NW Pac or N. IO if (flon(0).lt.180.0) then ! JTWC AOR ibasin = 3 !if a storm that forms in the wp and tracks across the dateline ! (rare, but has happens) then use JTWC errors and the ! correct wind radii CLIPER (JAK) else if (flon(0).ge.180.0) then ! NHC/CPHC AOR ibasin = 2 end if case('SH', 'sh', 'sp', 'SP') ! Southern Hemisphere ibasin = 4 case default write(lulog,'(a)') 'ibasin must be an integer value equal to 1,2,3, or 4' ierr = 9 return end select !Convert cbasin to lower case call lcase(cbasin,2) !Linearly interpolate or extrapolate to fill in missing intensity !forecast positions call linte(fvmx,fdtl,1,nt) !Calculate the speed and heading along the forecast track call shcal(flon,flat,ftime,fspd,fhead,mt,nt) !Calculate sin and cos factors for along/cross track rotations do k=0,nt alpha = dtr*(90.0-fhead(k)) cal(k) = cos(alpha) sal(k) = sin(alpha) enddo !Write adjusted forecast track/intensity to log file ! cc do k=0,nt ! cc write(lulog,100) ktime(k),flat(k),flon(k),fvmx(k),fdtl(k), ! cc + fspd(k),fhead(k) ! ! 100 format(1x,i3,6(f7.1)) ! cc enddo !Set ensemble variables to zero elat(:,:) = 0.0 elon(:,:) = 0.0 evmx(:,:) = 0.0 exdev(:,:) = 0.0 !Calculate track ensembles call tens(flat,flon,gpce,igpcef,indxg,demax, & & cal,sal,elat,elon,mte,ntr,mt,nt, & & lulog,ibasin) !Calculate the intensity of the ensembles call iens(fvmx,fdtl,ftime,elat,elon,evmx,mte,ntr,mt,nt, lulog,ibasin) !Read the climo wind model coefficients and !calculate the size parameter x deviations !for the ensembles call rens(flat,fspd,fhead,fvmx,r34(0,1:4),r50(0,1:4), & & r64(0,1:4),r100(0,1:4), & & exdev,mte,ntr,mt,nt,lulog,ibasin) iatcf=0 if (iatcf .eq. 1) then !Write ensemble tracks and intensities to file !in modified-old ATCF format. This call is basin specifi call atfile(flat,flon,fvmx,elat,elon,evmx, & & iyr,imon,iday,itime,mte,ntr,mt,ibasin) endif idprt=1 if (idprt .eq. 1) then !Calculate intensity probability distributions and write !related info to a file if requested. call idfile(flat,flon,fvmx,ftime,elat,elon,evmx, & & mte,ntr,mt,nt,wtprob,ibasin) endif write(lulog,*) 'Start call of wpcal' !Calculate wind probabilities call wpcal(flat,flon,fvmx,elat,elon,evmx,exdev,igpcef,gpce, & & isintpe,cbasin,sname,isnum,iyr,imon,iday,itime, & & ktime,plat,plon,wprobs,wprobi,wprob,nxy,mxy, & & mte,ntr,mt,nt,its,ite,kstime,ntsave,iprall,iptoar,ipgrid,& & irbc,wprobmin,r34,r50,r64,r100,ithresh,lulog) return end subroutine mcprob subroutine wpcal(flat,flon,fvmx,elat,elon,evmx,exdev,igpcef,gpce, & & isintpe,cbasin,sname,isnum,iyr,imon,iday,itime, & & ktime,plat,plon,wprobs,wprobi,wprob,nxy,mxy, & & mte,ntr,mt,nt,its,ite,kstime,ntsave,iprall,iptoar,ipgrid,& & irbc, wprobmin,r34i,r50i,r64i,r100i,ithresh,lulog) ! c********************************************************************** !This routine calculates the wind probabilties !using the ensemble tracks, intensities and wind radii !Wind radii climatology and persistence follow the same !methodology employed for the operational wind radii cliper. !Coefficients used to calculate the wind radii are derived for and !applied in the North Atlantic, Eastern North Pacific (including the !Central Pacific Basin) and the Western North Pacific ( including the !North Indian Ocean). !This methodology has three unique components as follows: !1. The wind radii account for the initial wind radii including !asymmetries. The radii will match the observed at t=0. The influence !of the observed radii decays exponentially with time with !an efolding time of 32-h. Set iaper=1 or 0 to turn this option on/off. !2. For ensemble members that have vmax greater than 95 knots a !linear bias correction is applied to the R34 (increase) and R64 !(decrease) from the vortex model to better fit observations. !An additional bias correction (6% reduction) for R64 is applied !at all intensities. Set ibias=1 or 0 to turn this option on/off. !3. The vortex model was developed from the operational wind radii !database, where the radii represent the maximum value in each of the !four quadrants. However, the model radii are being applied as the !average in each quadrant. The model radii values are adjusted !by the factor cmta. If no conversion is desired, set cmta to 1.0. !A study of HRD surface wind analyses suggests that cmta should be !0.85. !See wind radii flags section. ! c********************************************************************** ! c** ! c** Written by John Knaff(CIRA/CSU) and Mark DeMaria (NOAA/NESDIS) ! c** ! c** ! c** Last Modified ! c** August 4, 2004 - wind radii smaller than the climatological ! c** radius of maximum wind allowed, bias correction ! c** to R64 corrected SEE !* ! c** ! c** January 2015 - Time of arrival calclations added (MD) ! c** April 2020 - Converted to free format, pass through forecast radii ! c********************************************************************** implicit None !** Passed variables real :: flat(0:mt),flon(0:mt),fvmx(0:mt) real :: elat(mte,0:mt),elon(mte,0:mt),evmx(mte,0:mt) real :: exdev(mte,0:mt),gpce(0:mt) integer :: ktime(0:mt),kstime(0:ntsave) real :: plat(mxy),plon(mxy) real :: wprob(mxy,0:ntsave),wprobs(mxy,0:ntsave) real :: wprobi(mxy,0:ntsave) real :: r34i(0:mt,4),r50i(0:mt,4),r64i(0:mt,4),r100i(0:mt,4) real :: wprobmin integer irbc, mt character *2 cbasin character *10 sname !** Local variables integer, parameter :: mti=500 real :: elatt(mti),elont(mti),evmxt(mti) real :: elatt2(mti),elont2(mti) real :: exdevt(mti),eradt(4, mti,4) real :: elats(0:mti),elons(0:mti),evmxs(0:mti) real :: exdevs(0:mti),edtls(0:mti),evmxs2(0:mti) real :: etimes(0:mti),espds(0:mti),eheads(0:mti) real :: erads(0:mti), ftimes(0:mti) real :: ap, cfac, azdeg, dtl1, dtlv, dts, dtsp real :: tkk, thp, vmx1, vmxth, rlat1, rlon1, rmwp, rktime real :: radt, betat, fi, rnm, tensione real :: radarrm, radi, rlat2, rlon2 integer :: iisave(0:mti), je1, je2, jstm, ie, ii,je real :: radarr(4) real :: radarr_land(4) real :: rbc(4, mti) character *100 tmp_path, fn character *200 fnall,fntoa !Variables for wprobi calculation !Note: mxyl must be .ge. nxy integer, parameter :: mxyl=200000 real :: wptemi(mxyl),wptemc(mxyl) !Variables for time of arrival calculation integer, parameter :: mtoap=10 real :: toap(mtoap) real :: ptarrive(mxyl,mtoap),ptdepart(mxyl,mtoap),ptadtem(mtoap) integer :: icnttoap(mxyl),icnttodp(mxyl) integer :: i,iday,ie1,ie2,ierrt,iflag,igpcef,ihit,iis,imon,ind,intx,inty integer :: ipint,iprall,iptoar,ipgrid,isintpe,iskipxy,isnum,istm,ite,itfirst integer :: ithresh,itime,itlast,its,ivmxth,iyr integer :: k, kk integer :: lflag, luall, lulog, lutoa integer :: m,mte,mxy integer :: n,nt,nte,ntesave,ntoap,ntr,nts,ntsave,ntt,nxy integer, parameter:: mtetoa=1000,mtitoa=168 !Note: mtetoa is the max number of realizations that can be !used with the toar is requested. This is much smaller !than the number allowed without the toar calculation. !Also, ntr must be .le. mtetoa. real :: ftimeh(21) real :: elatst(mtetoa,0:mtitoa), elonst(mtetoa,0:mtitoa) real :: evmxst(mtetoa,0:mtitoa),edtlst(mtetoa,0:mtitoa) real :: rmwpt(mtetoa,0:mtitoa),thpt(mtetoa,0:mtitoa) real :: xtt(mtetoa,0:mtitoa),spdt(mtetoa,0:mtitoa) real :: apt(mtetoa,0:mtitoa),betatt(mtetoa,0:mtitoa) real :: allradarr(mtetoa,0:mtitoa,2,4,4) integer :: ittar(0:mtitoa),ittarf(0:mtitoa),ittarl(0:mtitoa) real :: r34t(4), r50t(4), r64t(4), r100t(4) !Arrays for code optimization for regular lat/lon grid integer, parameter :: mxe=2000,mye=1000 real :: plone(mxe),plate(mye) real :: dlone, dlate integer nxe, nye, iego common /egopt/ dlone,dlate,plone,plate,nxe,nye,iego !For iprall=1, set ipint to the time step interval for printing the information !(ipint=1 to print every time step). ipint = 1 !For iptoar=1, specify number of probabilities and prob values ntoap=9 toap(1) = 1.0 toap(2) = 5.0 toap(3) = 10.0 toap(4) = 25.0 toap(5) = 50.0 toap(6) = 75.0 toap(7) = 90.0 toap(8) = 95.0 toap(9) = 99.0 !ntoap = 3 !toap(1) = 25.0 !toap(2) = 50.0 !toap(3) = 75.0 !Unit number for realization file luall = 78 !Unit number for time of arrival file lutoa = 79 !Initialize intermediate save array to zero do i=0,mti iisave(i) = 0 enddo !Determine requested wind threshold if (ithresh .eq. 1) then vmxth = 34.0 elseif (ithresh .eq. 2) then vmxth = 50.0 elseif (ithresh .eq. 3) then vmxth = 64.0 elseif (ithresh .eq. 4) then vmxth = 100.0 else write(lulog,931) ithresh 931 format(' Illegal value of ithresh in routine wpcal: ',i6) stop endif !Check to see if local arrays are large enough if (nxy .ge. mxyl) then write(lulog,932) mxyl,nxy 932 format(' mxyl must be .ge. nxy in wpcal, nxy,mxyl=', & & i6,1x,i6) stop endif if (ntr .gt. mtetoa .and. iptoar .eq. 1) then write(lulog,933) ntr,mtetoa 933 format(/,' ntr must be .le. mtetoa when time of arrival', & & ' calculation is requested',/, & & ' ntr,mtetoa: ',i6,1x,i6,' toa calcuation skipped') iptoar=0 endif !Specify time interval (hr) for subdividing !ensembles for probability calculation dts = 1.0 !Calculate starting and ending indices !for sub-time intervals nts = ifix(float(its)/dts) nte = ifix(float(ite)/dts) ntesave = nte !Time array at sub-time intervals do k=0,nte ftimes(k) = dts*float(k) enddo write(lulog,*) 'wpcal started' !Calculate variables for intermediate saving !of probabilities iisave(nte) = ntsave if (ntsave .gt. 1) then do i=1,ntsave-1 ii = nts + ifix(float(kstime(i))/dts + 0.0001) iisave(ii) = i enddo endif !Limit time to last period that ensembles !are available ntt = ifix(float(ktime(nt) )/dts) if (ntt .lt. nte) nte=ntt !Forecast time array for interpolation do k=0,nt ftimeh(k+1) = float(ktime(k)) enddo !Error checking if (nte .gt. mti) then write(lulog,910) nte 910 format(' mti too small, value=',i6) stop endif !Initialize the wind probabilties to zero do m=0,ntsave do i=1,mxy wprob(i,m) = 0.0 wprobi(i,m) = 0.0 end do end do call getenv ( "NHCMESSAGES", tmp_path ) ! NHC messages directory ind = index( tmp_path, " " ) - 1 if (iprall .eq. 1) then !Open file for saving track/intensity/radii of all realizations fn ='mcrall_al991955_093000.dat' fn(8:9) = cbasin write(fn(10:22),661) isnum,iyr,imon,iday,itime 661 format(i2.2,i4.4,'_',3(i2.2)) if (ind .ge. 1) then fnall = tmp_path(1:ind)//"/"//fn else fnall = fn endif open(file=fnall,unit=luall,form='formatted', status='replace') !Write header to mcrall file write(luall,663) cbasin,isnum,iyr,imon,iday,itime,sname,igpcef 663 format(a2,i2.2,i4.4,1x,3(i2.2),1x,a10,1x,i1) dtsp = dts*float(ipint) write(luall,664) ftimes(nts),ftimes(nte),dtsp,ntr 664 format(f5.1,1x,f5.1,1x,f5.1,1x,i6) do k=0,mt rktime = float(ktime(k)) write(luall,665) rktime,flon(k),flat(k),fvmx(k),gpce(k) 665 format(f5.1,1x,f7.2,1x,f7.2,1x,f6.1,1x,f6.1) enddo else if (iprall .eq. 2) then !Open file for saving adeck output of all realizations fn ='mcrall_al991955_093000.dat' fn(8:9) = cbasin write(fn(10:32),666) isnum,iyr,imon,iday,itime 666 format(i2.2,i4.4,'_',3(i2.2),'_adeck.dat') fnall = tmp_path(1:ind)//"/"//fn open(file=fnall,unit=luall,form='formatted', status='replace') endif if (iptoar .eq. 1) then !Open mctoar file for writing time of arrival output fn ='mctoar_al991955_093000_000.dat' fn(8:9) = cbasin write(fn(10:22),661) isnum,iyr,imon,iday,itime ivmxth = nint(vmxth) write(fn(24:26),674) ivmxth 674 format(i3.3) if (ind .ge. 1) then fntoa = tmp_path(1:ind)//"/"//fn else fntoa = fn endif open(file=fntoa,unit=lutoa,form='formatted', status='replace') !Write header to mctoar file write(lutoa,673) cbasin,isnum,iyr,imon,iday,itime,sname,igpcef, & & ithresh,nxy 673 format(a2,i2.2,i4.4,1x,3(i2.2),1x,a10,1x,i1,1x,i1,1x,i6) dtsp = dts*float(ipint) write(lutoa,664) ftimes(nts),ftimes(nte),dtsp,ntr do k=0,mt rktime = float(ktime(k)) write(lutoa,665) rktime,flon(k),flat(k),fvmx(k),gpce(k) enddo write(lutoa,675) ntoap,(toap(n),n=1,ntoap) 675 format(18x,i5,10(f5.1)) !Initialize the probability of time of arrival/departure arrays to zero do i=1,nxy icnttoap(i) = 0 icnttodp(i) = 0 do n=1,ntoap ptarrive(i,n) = 0.0 ptdepart(i,n) = 0.0 enddo enddo endif if (nts .gt. ntt) then !Starting time for probability calculation is after the !end of the forecast period. Thus, the wind probabilities !are all equal to zero because there is no storm during !the requested time interval. return endif !Initialize input to interpolation routines lflag=0 tkk =0.0 eradt = 0 ! initialize radii holding array for times after valid forecast !Start the main ensemble loop do 99 i=1,ntr if (ioper .eq. 0) then if (mod(i,50) .eq. 0) write(lulog,*) 'realization ',i endif !Extract ensemble member track, intensity and x-deviations !for this case for interpolation do k=0,nt elatt(k+1) = elat(i,k) elont(k+1) = elon(i,k) evmxt(k+1) = evmx(i,k) exdevt(k+1)= exdev(i,k) eradt(1,k+1,:) = r34i(k,:) eradt(2,k+1,:) = r50i(k,:) eradt(3,k+1,:) = r64i(k,:) eradt(4,k+1,:) = r100i(k,:) enddo !Interpolate lat/lon to sub-time intervals IF (isintpe .EQ. 0) THEN iflag = 1 call xint(ftimeh,elatt,nt+1,iflag,lflag,tkk,fi,ierrt) call yint(ftimeh,elont,nt+1,iflag,lflag,tkk,fi,ierrt) iflag = 0 do kk=0,nte tkk = dts*float(kk) etimes(kk) = tkk call xint(ftimeh,elatt,nt+1,iflag,lflag,tkk,elats(kk),ierrt) call yint(ftimeh,elont,nt+1,iflag,lflag,tkk,elons(kk),ierrt) enddo ELSE CALL spline(ftimeh(1:nt+1),elatt(1:nt+1),nt+1, & & 99e31,99e31,elatt2(1:nt+1)) CALL spline(ftimeh(1:nt+1),elont(1:nt+1),nt+1, & & 99e31,99e31,elont2(1:nt+1)) DO kk=0,nte tkk = dts*float(kk) etimes(kk) = tkk CALL splintt(ftimeh(1:nt+1),elatt(1:nt+1), & & elatt2(1:nt+1),nt+1,tkk,tensione,elats(kk)) CALL splintt(ftimeh(1:nt+1),elont(1:nt+1), & & elont2(1:nt+1),nt+1,tkk,tensione,elons(kk)) ENDDO ENDIF !Save lat,lon at small time intervals for TOA calculation do kk=0,nte elatst(i,kk) = elats(kk) elonst(i,kk) = elons(kk) enddo !Calculate speed and heading at sub-time intervals call shcal(elons,elats,etimes,espds,eheads,mti,nte) !Interpolate intensity to sub-time intervals iflag = 1 tkk = 0.0 call xint(ftimeh,evmxt,nt+1,iflag,lflag,tkk,fi,ierrt) iflag = 0 do kk=0,nte tkk = etimes(kk) call xint(ftimeh,evmxt,nt+1,iflag,lflag,tkk,evmxs(kk),ierrt) enddo !Calculate distance to land at interpolated hourly positions do kk=0,nte call gdland_table(elons(kk),elats(kk),edtls(kk)) !Save dtl at small time steps for TOA calculation edtlst(i,kk) = edtls(kk) enddo !Reduce hourly intensity when realization crosses land between !12-hr forecast times (e.g., TC crosses narrow !peninsula or large island) if (irbc .eq. 2 .and. i .eq. 1) then iflag = 1 else call landcross(evmxs,edtls,nt,mti,dts,evmxs2) evmxs = evmxs2 end if do kk=0,nte !Save vmax at small time steps for TOA calculation evmxst(i,kk) = evmxs(kk) enddo !Interpolate x-deviations to sub-time intervals iflag = 1 tkk = 0.0 call xint(ftimeh,exdevt,nt+1,iflag,lflag,tkk,fi,ierrt) iflag = 0 do kk=0,nte tkk = etimes(kk) call xint(ftimeh,exdevt,nt+1,iflag,lflag,tkk, exdevs(kk),ierrt) enddo !Start the main time loop for this realization !Calculate wind radii at each small time step and then !check to see if probability grid points fall within !the radii of the specified wind threshold. call calc_ens_radii(mtetoa,mtitoa,nts, nte,mti,mt,i,ithresh, & & elats,etimes,elons,espds,eheads,evmxs,exdevs, & & r34i,r50i,r64i,r100i, vmxth, & & allradarr, rmwpt, apt,betatt, thpt, & & xtt, spdt) 99 continue ! c---------------------------------------------------------------------- !If flag irbc=1, use the radii bias correction calculated above !If not, set bias correction to 0.0 rbc = 1.0 if (irbc .eq. 1) then ! c---------------------------------------------------------------------- ! Calculate mean of realization wind radii (from radii-CLIPER above) write(lulog,*) 'Starting Ensemble Bias Calculation: ', irbc call calc_radii_bc(allradarr,mtetoa,mtitoa,& & mti, mt, ntr, flon, flat, eradt, ftimeh, rbc) elseif(irbc.eq.2) then write(lulog,*) 'Starting MRCL Bias Calculation: ', irbc call calc_radii_bc(allradarr,mtetoa,mtitoa,& & mti, mt, 1, flon, flat, eradt, ftimeh, rbc) endif if(irbc .gt. 0) then write(6,*) 'Applying bias correction' do i=1, ntr do kk=nts, nte do ii=1,4 allradarr(i,kk,1,ii,:) = allradarr(i,kk,1,ii,:) * rbc(:,kk+1) allradarr(i,kk,2,ii,:) = allradarr(i,kk,2,ii,:) * rbc(:,kk+1) enddo enddo enddo endif ! c---------------------------------------------------------------------- 701 format(i4,1x,a2,f5.0,1x,f5.1,1x,f6.2,1x,f6.2,1x,f6.1, & & 1x,f6.3,4(2x,4(f5.0,1x)),3(f5.1,1x),f5.0 ) do i=1,ntr do kk=nts,nte if (iprall .eq. 1 .and. mod(kk,ipint) .eq. 0) then ! undo the max to mean reduction for writing out to adeck like file r34t = allradarr(i,kk,1,1,:) * 1./cmta r50t = allradarr(i,kk,1,2,:) * 1./cmta r64t = allradarr(i,kk,1,3,:) * 1./cmta r100t = allradarr(i,kk,1,4,:) * 1./cmta write(luall, 701) & & i, 't=',etimes(kk),max(5.,evmxst(i,kk)),elatst(i,kk),& & elonst(i,kk),rmwpt(i,kk),xtt(i,kk), & & r34t(1),r34t(2),r34t(3),r34t(4), & & r50t(1),r50t(2),r50t(3),r50t(4), & & r64t(1),r64t(2),r64t(3),r64t(4), & & r100t(1),r100t(2),r100t(3),r100t(4), & & spdt(i,kk),betatt(i,kk),apt(i,kk),thpt(i,kk) elseif (iprall .eq. 2 .and. mod(kk,12) .eq. 0) then write(lulog, *) ' doconvWriteAidData is not included in this version' ! call doconvWriteAidData(luall,i,kk,ntesave,nte, & ! & etimes(kk),max(5.,evmxst(i,kk)),elatst(i,kk), & ! & elonst(i,kk),r34t,r50t,r64t,r100t) endif enddo enddo if (iprall .eq. 1 .or. iprall .eq. 2) then close(luall) endif if (ipgrid .eq. 0 .and. iptoar .eq. 0) return !Start the main ensemble loop do i=1,ntr !Initialize prob arrays for this ensemble member to zero wprobs = 0.0 wptemc = 0.0 wptemi = 0.0 do kk=nts,nte radarr(:) = allradarr(i,kk,1,ithresh,:) radarr_land(:) = allradarr(i,kk,2,ithresh,:) !Restore vortex struture parameters rlat1 = elatst(i,kk) rlon1 = elonst(i,kk) vmx1 = max(5.0,evmxst(i,kk)) dtl1 = edtlst(i,kk) ! Reset rmwp to be the smallest non-zero quadrant or the existing rmwp calc if (count(allradarr(i,kk,1,:,:).ne.0) .gt. 0) then rmwp = minval(allradarr(i,kk,1,:,:), allradarr(i,kk,1,:,:).ne.0) rmwp = min(rbuf*rmwp, rmwpt(i,kk)) else rmwp = 0 endif ap = apt(i,kk) betat = betatt(i,kk) thp = thpt(i,kk) !Find max radius radarrm = maxval(radarr) if ( radarrm .gt. 0.0) then !At least one of the wind radii are greater than zero, ! so checking the probability grid. !Check the probability grid points to see if they !fall within the radius of the requested wind threshold if (iego .ne. 0) then !Take advantage of evenly spaced lat/lon probability grid !Calculate indices to search that are close enough to the storm center !to be within a distance radarrm intx = 1 + ifix( radarrm/(dlone*60.0*cos(dtr*rlat1)) ) inty = 1 + ifix( radarrm/(dlate*60.0) ) istm = 1 + ifix( 0.5+((rlon1-plone(1))/dlone) ) jstm = 1 + ifix( 0.5+((rlat1-plate(1))/dlate) ) ! dbg_fmt ='(A, f5.1,1x,f6.1,1x,i3,1x,i3,1x,f6.1,1x,i2,1x,i2)' ! write(lulog,dbg_fmt) rlat1,rlon1,istm,jstm,radm,intx,inty: ': ! & rlat1,rlon1,istm,jstm,radarrm,intx,inty ie1=istm-intx ie2=istm+intx je1=jstm-inty je2=jstm+inty if (ie1 .lt. 1) ie1=1 if (je1 .lt. 1) je1=1 if (ie2 .gt. nxe) ie2=nxe if (je2 .gt. nye) je2=nye do je=je1,je2 do ie=ie1,ie2 !Calculate 1-D index from 2-D index call ofti(iego,nxe,nye,ie,je,ii) !ii = (je-1)*nxe + ie !write(lulog,669) ie,je,ii,plone(ie),plate(je), !+ plon(ii),plat(ii) !669 format('i,j,ii,plone,plate,plat,plon: ', !+ i3,1x,i3,1x,i6,2(1x,f5.1,1x,f6.1)) rlat2 = plat(ii) rlon2 = plon(ii) if (is_point_in_windfield(rlat2, rlon2, & & rlat1, rlon1,ap, betat, thp, rmwp, dtl1, vmx1,vmxth, & & radarr_land, radarr )) then !If you get to here, current x,y,t point for this realization !is within the wind threshold region wprobs(ii,ntsave) = 1.0 wptemc(ii) = 1.0 if (kk .eq. nts) then !Accumulate "initial" probabilities (at first time) wprob( ii,0) = wprob( ii,0) + 1.0 wprobi(ii,0) = wprobi(ii,0) + 1.0 endif if (iisave(kk) .ne. 0) then !Save the probabilities at this time step wptemi(ii) = 1.0 endif end if enddo enddo else !Use code for arbitrary lat/lon probability points do ii = 1,nxy !Check to see if this grid point for this ensemble !member has already been counted. If so, there is no !need to check it again. !if (wprobs(ii,ntsave) .eq. 1.0) go to 1000 !Note: This check was disabled for the wprobi calculation. rlat2 = plat(ii) rlon2 = plon(ii) if (is_point_in_windfield(rlat2, rlon2, & & rlat1, rlon1,ap, betat, thp, rmwp, dtl1, vmx1,vmxth, & & radarr_land, radarr )) then !If you get to here, current x,y,t point for this realization !is within the wind threshold region wprobs(ii,ntsave) = 1.0 wptemc(ii) = 1.0 if (kk .eq. nts) then !Accumulate "initial" probabilities (at first time) wprob( ii,0) = wprob( ii,0) + 1.0 wprobi(ii,0) = wprobi(ii,0) + 1.0 endif if (iisave(kk) .ne. 0) then !Save the probabilities at this time step wptemi(ii) = 1.0 endif endif enddo endif endif !Save wprobs for intermediate summations if (iisave(kk) .ne. 0) then iis = iisave(kk) do ii=1,nxy wprobs(ii,iis) = wprobs(ii,ntsave) wprobi(ii,iis) = wprobi(ii,iis) + wptemc(ii) enddo !Reinitialize wprobc array do ii=1,nxy wptemc(ii) = wptemi(ii) wptemi(ii) = 0.0 enddo endif enddo !Debug printing (check for -1 to disable) if (i .eq. -1) then write(lulog,*) 'nts,nte: ',nts,nte do kk=nts,nte write(lulog,801) kk,etimes(kk),elats(kk),elons(kk), & & evmxs(kk), espds(kk),eheads(kk),erads(kk) 801 format(i3,1x,f5.1,1x,5(f7.2,1x),1x,4(f5.1,1x)) enddo endif do m=1,ntsave do ii=1,nxy wprob(ii,m) = wprob(ii,m) + wprobs(ii,m) enddo enddo enddo cfac = 100.0/float(ntr) do m=0,ntsave do ii=1,nxy wprob(ii,m) = cfac*wprob(ii,m) wprobi(ii,m) = cfac*wprobi(ii,m) enddo enddo if (ntesave .gt. nte) then !Requested time period for probabilities extended past the !end of the forecast track. Check to see if intermediate !probabilities need to be added for the period from the end !of the forecast track to the end probability period. !Feb 2011 - Changed k=nte, to k=nte+1 (AS) if (ntsave .gt. 1) then do k=nte+1,ntesave iis = iisave(k) if (iis .gt. 0) then do ii=1,nxy wprob(ii,iis) = wprob(ii,ntsave) wprobi(ii,iis) = 0.0 enddo endif enddo endif endif if (iptoar .eq. 0) return !Processing for time of arrival/depature calculation iskipxy=0 do ii=1,nxy if (mod(ii,500) .eq. 0) then write(lulog,*) 'TOA processing for grid point ',ii,' of ',nxy endif !Check final cumulative probability to determine if this !grid point has enough points for TOA calculation if (wprob(ii,ntsave) .lt. wprobmin) then iskipxy = iskipxy + 1 icnttoap(ii) = -99 icnttodp(ii) = -99 do n=1,ntoap ptarrive(ii,n) = -99. ptdepart(ii,n) = -99. enddo cycle ! cycle loop early we don't need to do the rest endif rlat2 = plat(ii) rlon2 = plon(ii) !Initialize temp time array sum arrays to zero before interation loop ittarf(:) = 0 ittarl(:) = 0 do i=1,ntr !Initialize temp time array to zero before time loop do kk= 0,nte ittar(kk) = 0 enddo ihit=0 do kk=nts,nte rlat1 = elatst(i,kk) rlon1 = elonst(i,kk) vmx1 = max(5.0, evmxst(i,kk)) ap = apt(i,kk) betat = betatt(i,kk) thp = thpt(i,kk) if (count(allradarr(i,kk,1,:,:).ne.0) .gt. 0) then rmwp = minval(allradarr(i,kk,1,:,:), allradarr(i,kk,1,:,:).ne.0) rmwp = min(rbuf*rmwp, rmwpt(i,kk)) else rmwp = 0 endif if (is_point_in_windfield(rlat2, rlon2, & & rlat1, rlon1,ap, betat, thp, rmwp, edtlst(i,kk), vmx1,vmxth, & & allradarr(i,kk,2,ithresh,:), allradarr(i,kk,1,ithresh,:) )) then !If you get to here, current x,y,t point for this realization !is within the wind threshold region. ittar(kk) = 1 ihit = 1 endif enddo if (ihit .eq. 1) then !Find time of arrival and departure for this iteration !at this grid point itfirst = -1 itlast = -1 do kk=nts,nte if (ittar(kk) .eq. 1 .and. itfirst .lt. 0) itfirst=kk if (ittar(kk) .eq. 1 ) itlast =kk enddo if (itfirst .ge. 0) then ittarf(itfirst) = ittarf(itfirst) + 1 endif if (itlast .ge. 0) then ittarl( itlast) = ittarl( itlast) + 1 endif endif enddo !Calculate time of arrival, departure probabilities !for this grid point call ptadcal(nts,nte,ntoap,ittarf,ftimes,toap,icnttoap(ii),ptadtem) do n=1,ntoap ptarrive(ii,n) = ptadtem(n) enddo call ptadcal(nts,nte,ntoap,ittarl,ftimes,toap,icnttodp(ii),ptadtem) do n=1,ntoap ptdepart(ii,n) = ptadtem(n) enddo enddo write(lulog,*) 'iskipxy=',iskipxy, 'nxy=',nxy !Print results to mctoar file do ii=1,nxy write(lutoa,677) plon(ii),plat(ii),icnttoap(ii), & ! TODO: is this a ioper thing? !+ (ptarrive(ii,n),n=1,ntoap) ! WCOSS & (ptarrive(ii,n),n=1,ntoap),wprob(ii,ntsave) ! graphics 677 format(f8.3,1x,f8.3,1x,i6,9(f5.0),3x,f5.1) write(lutoa,678) plon(ii),plat(ii),icnttodp(ii), & & (ptdepart(ii,n),n=1,ntoap) 678 format(f8.3,1x,f8.3,1x,i6,9(f5.0)) enddo close(lutoa) end subroutine wpcal end module windspeed_model