program dtops_realtime !----------------------------------------------------------------! !--- ---! !--- Fortran program to compute the DTOPS percent chance of ---! !--- rapid intensification. This program expects 10 command ---! !--- line arguments: ---! !--- ---! !--- 1) cfile : file containing the logistic regression ---! !--- coefficients, means, and standard ---! !--- deviations for each dt and RI threshold ---! !--- 2) dt : number of hours over which del-V occurs ---! !--- 3) thresh : del-V threshold to consider for RI ---! !--- 4) vmax_init : TC max sustained winds at t = 0 h ---! !--- 5) tclat : TC latitude at t = 0 h ---! !--- 6) avni_delv : del-V from AVNI ---! !--- 7) hfai_delv : del-V from HFAI ---! !--- 8) dshp_delv : del-V from DSHP ---! !--- 9) lgem_delv : del-V from LGEM ---! !--- 10) emx_delp : del-P from EMX ---! !--- ---! !--- Usage: ./dtops_realtime.exe ${cfile} ${vmax_init} ${dt}\ ---! !--- ${thresh} ${tclat} ${avni_delv} ${hfai_delv} \ ---! !--- ${dshp_delv} ${lgem_delv} ${emx_delp} ---! !--- ---! !--- Compile statement: ---! !--- ---! !--- gfortran -O2 -o dtops_realtime.exe dtops_realtime.f ---! !--- ---! !--- Edit history: ---! !--- Matt Onderlinde - May 2017 - wrote original code ---! !--- Matt Onderlinde - Mar 2023 - updated for HAFS ---! !--- ---! !----------------------------------------------------------------! implicit none !--- Declare variables ------------------------------------------! !--- Declare input arguments character*900 :: cfile, vmax_initin, tclatin, dshp_delvin, dtin character*900 :: threshin, lgem_delvin, hfai_delvin character*900 :: avni_delvin, emx_delpin real :: dt, thresh, vmax_init, tclat, dshp_delv real :: lgem_delv, hfai_delv, avni_delv, emx_delp !--- Declare RI time period and threshold variables to be read from !--- each line of the coefficient file real :: curdt, curthresh !--- Declare logistic coefficient, mean, and standard devlation arrays real :: coef(14), mean(14), std(14) !--- Declare normalized predictor 1D array real :: pred(14) !--- Declare standard deviation cutoff factor real :: std_cutoff !--- Declare output percent chance of RI real :: rip !--- Declare other variables integer :: i character*15 :: dum !--- Done declaring variables -----------------------------------! !--- Get command line arguments ---------------------------------! call get_command_argument(1,cfile) call get_command_argument(2,dtin) call get_command_argument(3,threshin) call get_command_argument(4,vmax_initin) call get_command_argument(5,tclatin) call get_command_argument(6,avni_delvin) call get_command_argument(7,hfai_delvin) call get_command_argument(8,dshp_delvin) call get_command_argument(9,lgem_delvin) call get_command_argument(10,emx_delpin) read(dtin,*) dt read(threshin,*) thresh read(vmax_initin,*) vmax_init read(tclatin,*) tclat read(avni_delvin,*) avni_delv read(hfai_delvin,*) hfai_delv read(dshp_delvin,*) dshp_delv read(lgem_delvin,*) lgem_delv read(emx_delpin,*) emx_delp !--- Done getting arguments -------------------------------------! !--- Read in current coefficients from file ---------------------! open(10, file=trim(cfile), status='old') read(10,*) dum read(10,*) dum, std_cutoff do i = 1, 8 read(10,*) curdt, curthresh, coef(:) read(10,*) curdt, curthresh, mean(:) read(10,*) curdt, curthresh, std(:) if ( curdt .eq. dt .and. curthresh .eq. thresh ) then exit end if end do close(10) !--- Done reading coefficients ----------------------------------! !--- Normalize by mean and standard deviation -------------------! !--- (e.g., compute normalized predictor array) pred(1) = coef(1) !--- Intercept term pred(2) = (avni_delv - mean(2))/(std(2)+0.0000001) !--- AVNI term pred(3) = (hfai_delv - mean(3))/(std(3)+0.0000001) !--- HFAI term pred(4) = (dshp_delv - mean(4))/(std(4)+0.0000001) !--- DSHP term pred(5) = (lgem_delv - mean(5))/(std(5)+0.0000001) !--- LGEM term pred(6) = (emx_delp - mean(6))/(std(6)+0.0000001) !--- EMX term pred(7) = (vmax_init*cos(tclat*0.0174533) & - mean(7))/(std(7)+0.0000001) !--- VMAX*cosine(TC_latitude*deg2rad) term pred(8) = ((vmax_init**2.0) - mean(8))/(std(8)+0.0000001) !--- VMAX squared term pred(9) = ((dshp_delv*hfai_delv) - mean(9))/(std(9)+0.0000001) !--- DSHP*HFAI term if ( avni_delv .gt. (std_cutoff*std(2)) ) then pred(10) = 1.0 else pred(10) = 0.0 end if if ( hfai_delv .gt. (std_cutoff*std(3)) ) then pred(11) = 1.0 else pred(11) = 0.0 end if if ( dshp_delv .gt. (std_cutoff*std(4)) ) then pred(12) = 1.0 else pred(12) = 0.0 end if if ( lgem_delv .gt. (std_cutoff*std(5)) ) then pred(13) = 1.0 else pred(13) = 0.0 end if if ( emx_delp .lt. (-1.0*std_cutoff*std(6)) ) then !--- Notice for EMX we need to know if del-P pred(14) = 1.0 !--- is more negative than 1.25 STDs else pred(14) = 0.0 end if pred(10) = (pred(10) - mean(10))/(std(10)+0.0000001) pred(11) = (pred(11) - mean(11))/(std(11)+0.0000001) pred(12) = (pred(12) - mean(12))/(std(12)+0.0000001) pred(13) = (pred(13) - mean(13))/(std(13)+0.0000001) pred(14) = (pred(14) - mean(14))/(std(14)+0.0000001) !--- Done normalizing -------------------------------------------! !--- Compute and report DTOPS RI percentage ---------------------! rip = pred(1) + sum(pred(2:14)*coef(2:14)) !--- Sum of betas in the logistic equation rip = 1.0 / (1.0 + exp(-1.0 * rip)) !--- Final logistic RI probability write(*,'(i4)') nint(rip*100.0) !--- Done computing RI percentage -------------------------------! end program