program iships
c This program reads in data from the istormcard.dat, itrack.dat and
c lsdiag.dat files which contain the information for the SHIPS
c intensity forecast. The forecast is written to the file
c ships.dat for the new ATCF format, and to ships.tst in a modified
c version of the old ATCF format.
c
c Ver 23.0.2
c
c This version can be run for any basin (AL, EP, CP, WP, IO, SH)
c
c The over-ocean forecasts are modified to take into
c account the effect of land.
c
c This version also includes forecasts out to 5 days,
c and writes the output in the new ATCF format.
c
c Modified: Jul 27, 2001 to add Atlantic rapid intensity index
c Modified: Apr 10, 2002 to add new SHIPS predictors for 2002 season
c Modified: Jul 16, 2002 to add new RII routine with 24-hr input parameters
c and for parallel forecast with OHC predictors
c Modified: May 9, 2003 to add new SHIPS predictors for 2003 season
c Modified: May 21, 2004 for 2004 season (new coefficients, RI index)
c Modified: May 17, 2005 for 2005 season (new coefficients, RI index,
c 6 hour data)
c Modified: Jul 14, 2005 to correct persistence for storms moving
c back over water
c Modified: May 15, 2006 for 2006 season
c 1) Speed adjusted MPI
c 2) New T250 predictor
c 3) Forecast lat/speed used in Cione SST cooling
c 4) 36 N restriction removed in Cione cooling
c 5) Forecast divergence field used
c 6) Order of predictors adjusted
c 7) New coefficients
c Modified: May 19, 2006 to add Logistic Growth Equation (LGE) forecast
c Modified: Jul 3, 2006 to add ATCF number to SHIPS output
c Modified: Aug 17, 2006 to add additional GOES error checking
c Modified: Aug 29, 2006 Still more GOES error checking (std check)
c Modified: Sep 21, 2006 Experimental shear version and annular
c hurricane index added
c Modified: Mar 30, 2007 for 2007 season
c 1) SHRD replaced by SHDC
c 2) New TWAC predictor added
c 3) RHHI replaced by RHMD
c 4) New coefficients (based on 1982-2006)
c 5) VMAX added to LGEM forecast
c 6) Experimental run turned off
c Modified: May 22, 2008 for 2008 season
c 1) Some variable names changes to avoid compiler warnings
c 2) New RII routines
c 3) New coefficients (based on 1982-2007)
c Modified: May 15, 2009 for the 2009 season
c 1) New coefficients (based on 1982-2008)
c 2) Perturbation model eliminated
c (satellite input included with main model)
c 3) OHC added to east Pacific
c 4) New quality control for GOES and OHC data
c (proxy or sample mean used when missing)
c 5) New predictor based on shear direction
c added to SHIPS and LGEM
c 6) Satellite data added to LGEM (GOES and OHC)
c 7) Code for "experimental" SHIPS removed
c Modified: Mar 11, 2010 to add PrSEFoNe model (prob of secndry eywll frmtn)
c 1) Small block of code to fill feature matrix required
c by model added to main program body.
c 2) PrSEFoNe driver subroutine added below annular hurricane
c index code. All required subroutines and functions
c (PrSEFoNe module) added below AHI module.
c Modified: May 25,2010 for 2010 season
c 1) New predictor SHGC - SHCGBAR(SHDC)
c Modified: Apr 19, 2011 for 2011 season
c 1) New predictor TADV added
c 2) Print modified to exclude track for OFCP and OFPI
c 3) Extratropical classifcation routine added
c
c Modified: Dec 19, 2011 to add W. Pacific option
c 1) Lon option changed to deg E pos, deg W neg
c 2) dtland.f used for all land calculations
c 3) Revised decay.f and fland.f routines used
c 4) Some subroutines removed from iships.f file:
c lgeim.f, tspdcal.f
c 5) lcheck modified to use dtland instead of aland
c
c Modified: June 22, 2012 for 2012 season (MD)
c 1) new coefficient files added
c 2) Old perturbation SHIPS code removed
c 3) TCLIP routine added
c 4) Storm type written to ATCF file
c 5) 6 hr forecast interval to ATCF file instead of 12 hr
c
c Modified: Oct 2013 for global version (MD) that includes IO and SH options (MD).
c 1) gbland replaces dtland calls
c 2) Options for IO and SH basins
c 3) Basin determined from initial storm location rather than ATCF ID (except AL and SH)
c 4) MPI routines removed from iships.f file
c 5) jdate routine replaced by jday utility routine
c 6) sgnlon removed from lcheck input variable list
c
c Modified: Jul 2014 (MD)
c 1) New RII predictor in LGEM (required prelim call to ric2o)
c 2) Patch routine separated from code
c Modified: May 2015 (MD) for version 15.0.0
c 1) New G200 predictor added to SHIPS and LGEM
c 2) Bug fix for GFS vortex tendency calculation for LGEM
c Modified: Jun 2015 (MD) for version 15.0.1
c 1) Call to patch0 for vsstl added before
c prelim ric2o call
c 2) Typo corrected in ann hurr index code
c (idiv replaced by div)
c Modified: May 2016 (MD and GC) for version 16.2.0
c 1) New RII routines from J. Kaplan and C. Rozoff
c 2) Module for picking best SST and OHC fields from G. Chirokova
c 3) Updated coefficients for 2016
c 4) Model TPW added
c 5) New IPRC processing
c Modified: Nov 2016 (MD) for version 16.3.0 for Cray conversion
c 1) tcliper forecast eliminated since that is standalone on the Cray
c Modified: Apr 2017 (MD) for version 17.0.0. Mods for 2017
c Hurricane season. New RII routines, new
c SHIPS/LGEM coefficient files. RII
c predictor for LGEM 55 kt in 48 hr
c instead of 30 kt in 24 hr.
c Modified: Feb 2018 (GC) Added option to use multiple SSTs simultaneosly.
c (MD) John Kaplan's 2018 RII routines added.
c Chris Rozoff's 2018 RII routines added.
c Some legacy code removed.
c SHIPS:
c New predictors RIIP, PW06
c RSST replaced by DSTA (area avg daily SST)
c RHCN replaced by NOHC
c LGEM:
c Same as for SHIPS, except RIIP already included
c Modified: Jun 2018 (MD/SS) Changed XSST to CSST
c Fixed SST printed to SHIPS txt file
c PW06 removed from SHIPS
c Modified: Nov 2019 (GC) Modified way to use multiple SST and OHC variables.
c All SST variables were renamed and redefined
c New logic was added to use different SST
c OHC variables independently for all different models, inluding SHIPS.
c LGEM, etclass, ric2o. RIIs, PrSEF, and AHI
c Modified: Mar 2019 (GC) Modified to remove DAVT variables. This is the starting version for making
c operational 2019 version; version for J Kaplan, and version for
c JHT reruns with structure and eye predcitors
c Modified: Apr 2019 (GC) 1) grouped SST and model_sst variables and parameters in structures to
c simplify subroutine calls. See
c sst_module.f90 for definition of structures
c 2) moved annular hurricane index into a separate file ahi.f
c 3) ahi_sst is now put together in
c iships.f and passed to ahi.f
c Note: Search for SELECT SST to find code location for setting SST for each model
c
c Modified: Apr 2019 (GC) added option to use updated ric2o.f where potmin
c is defined in the calling code and Vmax**2/70.
c
c Modifed May 7, 2019 (MD) v19.1.0: 2019 SHIPS/LGEM coefficients
c New Kaplan RII routines
c Small bug fixes
c
c Modified May 13, 2019 (SS) v19.1.1: Changed SHIPS txt header to 2019
c
c Modified Sep 20, 2019 (SS) v19.1.2: Fixed etclass input sst to be an integer (isst_etcl)
c Modified Apr 14, 2020 (SS) v20.0.0: 2020 SHIPS/LGEM coefficients, new Kaplan routines, updated for 7-day
c forecast, reformatted text file (inc. scrubbing day 5-7 track), changed
c DIS to N/A in txt file
c Modified Mar 31, 2021 (MD) v21.0.0 Coefficient names updated for 2021 season.
c Calls to new RII routines for 2021 season.
c Modified Aug 1, 2021 (MD) v22.0.0 Test version for 2022 CIRA parallel runs
c
c Modified Feb 5, 2023 (MD) v23.0.0 (1) Updated for 2023 season.
c (1) J. Cione cooling removed from AL SHIPS/LGEM
c (2) SST and LAT added as SHIPS predictors
c (3) LAT added as LGEM predictors
c (4) SHIPS/LGEM coefficients updated with developmental data thru 2022
c (5) J. Kaplan RII routines updated with developmental data thru 2022
c (6) Generic names (rapidga and rapidge) for J. Kaplan RII routines
c Modified Feb 6, 2023 (MD) v23.0.1 (1) Predictor table changed to add LAT/SST contributions (added with MPI term)
c (2) Code changes for NCO standards
c Modified Feb 14, 2023 (MD) v23.0.2. (1) Remove writing predictor contributions to log file and ships.txt file
c for times with no track forecast
c (2) Remove SHIPS forecasts for times with no track forecast
c
ccccccccccccccccccccccccccccccc
c structure HOWTO
c
c 4 data structureas are decalred in th emodule sst_module.f90
c structures of
c
c type data_sst_params
c type data_sst_set
c type model_sst_params
c type model_sst_set
c
c to use structure in iships:
c
c 1) use sst_module, and get acess to structures and
c subroutines
c
c use sst_module, only : select_best_sst,
c + get_model_sst, get_cooling_mpi,
c + data_sst_set, model_sst_set,
c + data_sst_params, model_sst_params,
c + init_data_sst_set,
c + init_model_sst_set,
c + get_sst_print_variables
c
c 2) declare structure:
c reynolds wsst
c type(data_sst_set),dimension(0:mft) :: wsst_set
c
c
c
c get single variable from the structure:
c test_arr = wsst_set%isst_inp
c
c get single variable values from the structure:
c NOTE: the structure contains a single element from each SST array; thus place index after
c structure name
c test_arr = wsst_set(k)%isst_inp
ccccccccccccccccccccccccccccccc
c inlude sst_module and
c list structures and subroutines to be used
use sst_module, only : select_best_sst,
+ get_model_sst, get_cooling_mpi,
+ data_sst_set, model_sst_set,
+ data_sst_params, model_sst_params,
+ init_data_sst_set,
+ init_model_sst_set,
+ get_sst_print_variables
c inlude ohc_module, and
c list subroutines to be used
use ohc_module, only : select_best_ohc,get_model_ohc
include 'dataformats.inc'
c
parameter (nvar=28)
parameter (mft=28)
dimension ilat(-2:mft),ilon(-2:mft),ishr(0:mft)
dimension ishtd(0:mft),ivmaxb(0:mft)
dimension idist(0:mft)
dimension iu200(0:mft)
c min potential to use for ric2o
c - setup separately for each basin
real potmin
c input SST data - these are variables read from lsdiag
c reynolds wsst
type(data_sst_set),dimension(0:mft) :: wsst_set
c reynolds daily
type(data_sst_set),dimension(0:mft) :: dsst_set
c reynolds daily averaged
type(data_sst_set),dimension(0:mft) :: dsta_set
c ncoda daily
type(data_sst_set),dimension(0:mft) :: nsst_set
c ncoda daily averaged
type(data_sst_set),dimension(0:mft) :: nsta_set
c input SST data parameters
c reynolds wsst
type(data_sst_params) :: wsst_params
c reynolds daily
type(data_sst_params) :: dsst_params
c reynolds daily averaged
type(data_sst_params) :: dsta_params
c ncoda daily
type(data_sst_params) :: nsst_params
c ncoda daily averaged
type(data_sst_params) :: nsta_params
c model sst sets
c for each piece there is an option for a) SST and b) MPI (vsst)
c each structure contains multiple elements including arrays
c arrays sre defined inside structure
c SHIPS
type(model_sst_set),dimension(0:mft) :: ships_sst_set
c LGEM
type(model_sst_set),dimension(0:mft) :: lgem_sst_set
c RIC2O
type(model_sst_set),dimension(0:mft) :: ric2o_sst_set
c etclass sub
type(model_sst_set),dimension(0:mft) :: etcl_sst_set
integer,dimension(0:mft) :: isst_etcl
c J Kaplan RII
type(model_sst_set),dimension(0:mft) :: kaplan_sst_set
c C Rozof RII
type(model_sst_set),dimension(0:mft) :: rozof_sst_set
c PrSEF
type(model_sst_set),dimension(0:mft) :: prsef_sst_set
c AHI
type(model_sst_set),dimension(0:mft) :: ahi_sst_set
c lgem
c For lgem keep sst and vsst variables in addition to structure to
c not mess with common
dimension sst_lgem(0:mft)
dimension vsst_lgem(0:mft)
c Array for SHIPS SST predictor
dimension sst_ships(0:mft)
c model sst parameters
type(model_sst_params) :: ships_sst_params
c LGEM
type(model_sst_params) :: lgem_sst_params
c RIC2O
type(model_sst_params) :: ric2o_sst_params
c etclass sub
type(model_sst_params) :: etcl_sst_params
c J Kaplan RII
type(model_sst_params) :: kaplan_sst_params
c C Rozof RII
type(model_sst_params) :: rozof_sst_params
c PrSEF
type(model_sst_params) :: prsef_sst_params
c AHI
type(model_sst_params) :: ahi_sst_params
dimension it150(0:mft),it200(0:mft),it250(0:mft),ig200(0:mft)
dimension irefc(0:mft),iz850(0:mft)
dimension id200(0:mft),iepos(0:mft),itadv(0:mft)
dimension iv500(0:mft),iv300(0:mft),iv20c(0:mft),itgrd(0:mft)
dimension irhlo(0:mft),ipslv(0:mft)
dimension irhhi(0:mft),irhmd(0:mft)
c
dimension igoes(0:mft),igoesm1(0:mft),igoesm3(0:mft)
dimension igoesxx(0:mft)
dimension ipw01(0:mft),ipw03(0:mft),ipw04(0:mft)
dimension ipw06(0:mft),ipw14(0:mft),ipw19(0:mft)
dimension pw06(0:mft)
dimension iohc_inp_ncoda(0:mft)
dimension iohc_climo_ncoda(0:mft)
dimension iohc_best_ncoda(0:mft)
dimension ohc_final_ncoda(0:mft)
dimension iohc_ships(0:mft)
dimension iohc_lgem(0:mft)
dimension iohc_ric2(0:mft)
dimension iohc_kaplan(0:mft)
dimension iohc_rozof(0:mft)
dimension iohc_prse(0:mft)
dimension iohc_etcl(0:mft)
c
c PARMOD+ For RII with TPW, IRPC and CFLUX (TIC) input
c Also for Lightning input (LI)
parameter (mxh=100)
dimension iyrll(0:mxh),imonll(0:mxh),idayll(0:mxh),itimell(0:mxh)
dimension ivmaxll(0:mxh)
dimension clr01(0:mxh),clr23(0:mxh)
dimension istpw(0:mxh),imtpw(0:mxh),iaipc(0:mxh),icflx(0:mxh)
dimension iusetpw(0:mxh)
dimension ipc00(0:mft),ipcm1(0:mft),ipcm3(0:mft)
c PARMOD-
c
dimension rlat(-2:mft),rlon(-2:mft)
dimension ftimec(-2:mft),cxt(-2:mft),cyt(-2:mft),cmagt(-2:mft)
c
dimension var(nvar),coef(nvar,mft),vart(nvar,mft)
dimension xbar(nvar,mft),xsig(nvar,mft)
dimension ybar(mft),ysig(mft)
dimension v(0:mft),delv(0:mft)
dimension iv(0:mft)
c
c PARMOD+ For RII with multipe times and consensus (MTC)
parameter (mxmul=8)
dimension idivc(0:mft),ishrg(0:mft),ienss(0:mft)
dimension iepss(0:mft),ishrd(0:mft)
dimension probridisc(mxmul),probribay(mxmul)
dimension probrilog(mxmul), probricon(mxmul)
dimension idipc(0:mxh)
dimension ishrs(0:mft), ie000(0:mft), ieneg(0:mft)
c
c PARMOD-
c
dimension ilatt(0:mft),ilont(0:mft)
c
c Arrays for printing intensity change per variable
dimension dvvar(nvar,mft)
character *20 cdvlab(nvar)
c
character *4 cv(0:mft),cshr(0:mft)
character *4 cshtd(0:mft)
character *4 cdist(0:mft),crefc(0:mft)
character *4 crhlo(0:mft),cu200(0:mft)
character *4 cz850(0:mft),cd200(0:mft),ctadv(0:mft)
character *4 crhhi(0:mft),cepos(0:mft)
character *4 crhmd(0:mft)
character *5 clat(0:mft),clon(0:mft),ct200(0:mft),ct250(0:mft)
character *5 cg200(0:mft)
character *4 cmagc(0:mft)
character *4 cnohc(0:mft)
c
character *4 sname4
character *4 tem4
character *4 tmname
character *5 varlab(0:nvar)
character *6 slab,natcf
character *8 natcf8
character *10 sname
character *256 fnis,fnls,fnit,fnco,fnat,fnsh,fnlg,fnts
character *256 fncop
character *256 coef_location
character *80 iline80
character *180 iline
c
c Arrays for decay program
parameter(idmx=mft+1)
dimension rlatd(idmx),rlond(idmx),ftime(idmx),vmaxo(idmx)
dimension vmaxd(idmx),ddland(idmx)
dimension ivd(0:mft)
character *4 cvd(0:mft)
c
c Variables for new ATCF output
dimension ifship(11,3)
character *8 strmid
character *10 aymdh
c
c Variables for LGE model
parameter (ilmx=mft+1)
dimension vmaxl(ilmx)
dimension t200(0:mft),t250(0:mft),g200(0:mft)
dimension ivl(0:mft)
character *4 cvl(0:mft)
c
c Variables for modified shear version and pre-operational tests
dimension itwac(0:mft),ishdc(0:mft),isddc(0:mft)
dimension twac(0:mft), shdc(0:mft),sdir(0:mft)
character*2 runid,coyear
character*4 cshdc(0:mft),cshgc(0:mft)
character*4 ctwac(0:mft)
c
c Variables for SHGC and transformed shdc
dimension ishgc(0:mft)
dimension shdct(0:mft),shgc(0:mft)
c
c Variables for annular hurricane index
character *8 tcfid
character *10 stname
character *34 fn_ships
character *34 fn_ird1
character *34 fn_ird2
character *34 fn_iri1
character *34 fn_iri2
c
c Variables for ET classification routine
character *4 stype(0:mft)
c
c Variables for TCLIPER
parameter(ndt=40)
dimension tcplat(0:ndt), tcplon(0:ndt), tcpvmx(0:ndt)
dimension itcplat(0:ndt),itcplon(0:ndt),itcpvmx(0:ndt)
dimension itimet(0:ndt)
character *4 tcpstype(0:ndt)
c
c Variables for 7 day forecasts
parameter (mft5=20,mft7=28)
dimension ilat7(mft7),ilon7(mft7)
dimension iv7(mft7),ivd7(mft7)
c
c Variables for RII
dimension probri(mxmul)
dimension iprobri(mxmul)
character *256 fned
integer lued
dimension ithres(mxmul)
data ithres /20, 25, 30, 35, 40, 45, 55, 65/
dimension iendtau(mxmul)
data iendtau /12, 24, 24, 24, 24, 36, 48, 72/
c
c Variables for global version
character *1 latcon,loncon
character *50 temchar
c variable determining which SST variables shoudl be used
c
character *4 cships_ohc_type
character *4 clgem_ohc_type
character *4 cric2_ohc_type
character *4 cetcl_ohc_type
character *4 ckaplan_ohc_type
character *4 crozof_ohc_type
character *4 cprse_ohc_type
c value to use as missing for cvhar variable
character *4 cmiss
c ************************************************************************
c Variables for PrSEFoNe (probability of secondary eyewall formation)
c NATL: 11 features from lsdiag. Include DTL for 12. PC4 is
c calculated later. Data needed for 4 forecast times.
dimension psef_feat_vec(12,4) ! NATL feature set
dimension iflag_feat(4)
dimension ipenc(0:mft),ivmpi2(0:mft) ! PrSEFoNe model needs these
integer imiss_feat,iflag_feat_ir
c ************************************************************************
c
c Array for climatological goes predictors
parameter (ngoes=16)
dimension igoessm(0:mft)
data igoessm /0,-365,145,-266,191,
+ 71,65,58,50,40,25,
+ -364,-402,15,
+ -481,-405,56,
+ 9999,9999,9999,9999,
+ 9999,9999,9999,9999,
+ 9999,9999,9999,9999/
c
c Max lon (deg W neg) for calling annular routine
data rlonam /-35.0/
c
c ++ Common blocks for main iships.f program for lgecal routine
common /lgestr/ aday,spdx,pslv,per,vmx,pc20,ohc_thresh_lgem,
+ pprobric2
common /lgetdi/ ishr,iepos,iz850,id200,ishdc,itadv,
+ ipw06,iohc_best_ncoda,iohc_lgem
common /lgetdr/ sst_lgem,vsst_lgem,rlat,rlon,t200,t250,twac,sdir,
+ shdct,shgc,g200
common /lgepar/ rmiss,imiss,ioper
common /lgetst/ coyear
c
c Initialize various arrays
data delv /29*0.0/
data v /29*0.0/
c
data dtr /0.017453/
data luis,luls,luit,luco,luat,lush,lulg,luts,lued
+ /30,31,32,33,34,35,36,37,39/
c
data lucop /38/
c
type ( AID_DATA ) comRcd, tauData
c flag printed to the log to describe which OHC was used
integer ipohc
real rmpi,rmpic,rmpip
c
c Set ioper=1 for operational run where coefficient files use the getenv command
c or ioper=0 for non-operational runs where coefficient files are copied
c to the local directory where the code is run
c or ioper=-1 to specify new SHIP/DSHP/LGEM model names, but otherwise the same as ioper=0
ioper= 1
c
c Set run id number for test runs
runid='08'
c
c Set coefficient year for test runs
coyear='24'
c
c Set ipf7 = 1 to print 7 day forecasts to ships.tst file
c or ipf7 = 0 to print 5 day forecasts to ships.tst file
ipf7 = 1
c
c get COEF Env Variable
call getenv( "SHIPS_COEF", coef_location )
c PARMOD+ For RII
c flag for use observed or model itpwtype
c 0 - observed; 1 - model
itpwtype = 1
c
c Set ipcnew = 1 to use new IRPC variables (from K. Musgrave code)
c to set needed variables for RII routines. Set ipcnew=0 to skip.
ipcnew = 1
c
c PARMOD-
c TODO - move to config file
isst_max_lag_days = 13
iohc_max_lag_days = 13
c Missing Values to use
! integer
imiss=9999
! real
rmiss=9999.
!character
cmiss = 'none'
c initialize values for all sst sets
call init_data_sst_set(mft, wsst_set,wsst_params,
+ rmiss, imiss, cmiss)
call init_data_sst_set(mft, dsst_set,dsst_params,
+ rmiss, imiss, cmiss)
call init_data_sst_set(mft, dsta_set,dsta_params,
+ rmiss, imiss, cmiss)
call init_data_sst_set(mft, nsst_set,nsst_params,
+ rmiss, imiss, cmiss)
call init_data_sst_set(mft, nsta_set,nsta_params,
+ rmiss, imiss, cmiss)
c
c initialize values for all model sst sets
call init_model_sst_set(mft, ships_sst_set,ships_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, lgem_sst_set,lgem_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, ric2o_sst_set,ric2o_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, etcl_sst_set,etcl_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, kaplan_sst_set,kaplan_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, rozof_sst_set,rozof_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, prsef_sst_set,prsef_sst_params,
+ rmiss, imiss, cmiss)
call init_model_sst_set(mft, ahi_sst_set,ahi_sst_params,
+ rmiss, imiss, cmiss)
c
c SST to use -----------------SELECT SST HERE - START
c - wsst
c - dsst
c - nsst
c - dsta
c - nsta
c SST to use for each model as SST
ships_sst_params%csst_type = 'nsta'
lgem_sst_params%csst_type = 'nsta'
ric2o_sst_params%csst_type = 'nsta'
etcl_sst_params%csst_type = 'wsst'
kaplan_sst_params%csst_type = 'nsta'
rozof_sst_params%csst_type = 'wsst'
prsef_sst_params%csst_type = 'wsst'
ahi_sst_params%csst_type = ships_sst_params%csst_type
c SST to use for each model as VSST
ships_sst_params%cvsst_type = 'nsta'
lgem_sst_params%cvsst_type = 'nsta'
ric2o_sst_params%cvsst_type = 'nsta'
etcl_sst_params%cvsst_type = 'wsst'
kaplan_sst_params%cvsst_type = 'nsta'
rozof_sst_params%cvsst_type = 'wsst'
prsef_sst_params%cvsst_type = 'wsst'
ahi_sst_params%cvsst_type = ships_sst_params%cvsst_type
c SST to use -----------------SELECT SST HERE - END
c OHC to use
c 0 - nohc
cships_ohc_type = 'nohc'
clgem_ohc_type = 'nohc'
cric2_ohc_type = 'nohc'
cetcl_ohc_type = 'nohc'
ckaplan_ohc_type = 'nohc'
crozof_ohc_type = 'nohc'
cprse_ohc_type = 'nohc'
c
c Specify intensity (kt) for which storm is considered dissipated
vdiss = 15.0
ivdiss = ifix(vdiss)
c
c Forecast time interval (hr)
delt = 6.0
idelt = ifix(delt)
c
do k=0,ndt
itimet(k) = k*idelt
enddo
c
c Specify interval for writing ATCF output (minc=1 for 6 hr, =2 for 12 hr)
minc = 2
c
c Open the input and output files
fnis = 'istormcard.dat'
fnls = 'lsdiag.dat'
fnit = 'itrack.dat'
c
fnat = 'ships.dat'
fnts = 'ships.tst'
fnsh = 'ships.txt'
fnlg = 'ships.log'
c added ships output edeck file (RI)
fned = 'ships.edk'
c
open(unit=luis,file=fnis,form='formatted',status='old',err=900)
open(unit=luls,file=fnls,form='formatted',status='old',err=900)
print*, 'OPENED LSDIAG ======================'
open(unit=luit,file=fnit,form='formatted',status='old',err=900)
c
open(unit=lulg,file=fnlg,form='formatted',status='replace',
+ err=900)
open(unit=luat,file=fnat,form='formatted',status='replace',
+ err=900)
open(unit=luts,file=fnts,form='formatted',status='replace',
+ err=900)
open(unit=lush,file=fnsh,form='formatted',status='replace',
+ err=900)
c added edeck file
open(unit=lued,file=fned,form='formatted',status='replace',
+ err=900)
c
print*,"runid=", runid
c Read the stormcard file
read(luis,100) isyr,ismon,isday
read(luis,102) istime
read(luis,103) ilat0,ilatm12,latcon
read(luis,103) ilon0,ilonm12,loncon
read(luis,104) iline80
read(luis,104) iline80
read(luis,*) ivmx0,ivmxm12
read(luis,*) ihead
read(luis,*) ispeed
read(luis,107) sname
read(luis,101) natcf
101 format(a6)
100 format(3(i2))
102 format(i2)
103 format(i4,1x,i4,1x,a1)
104 format(a80)
105 format(i3,2x,i3)
106 format(i3)
107 format(1x,a10)
c
if (latcon .eq. 'S') then
ilat0 = -iabs(ilat0)
ilatm12 = -iabs(ilatm12)
endif
c
c Convert lon to 0 to 360 convention for now
if (loncon .eq. 'W') then
ilon0 = 3600-iabs(ilon0)
ilonm12 = 3600-iabs(ilonm12)
endif
c
c Calculate 4 digit year
if (isyr .gt. 50) then
isyr4 = isyr + 1900
else
isyr4 = isyr + 2000
endif
c
c Calculate 8 character ATCF ID (with 4 digit year)
natcf8(1:4) = natcf(1:4)
c
read(natcf(5:6),102) iatcfyr2
if (iatcfyr2 .gt. 50) then
iatcfyr4 = iatcfyr2 + 1900
else
iatcfyr4 = iatcfyr2 + 2000
endif
c
write(natcf8(5:8),109) iatcfyr4
109 format(i4.4)
c
c Create other variables for output
write(aymdh,600) isyr4,ismon,isday,istime
600 format(i4.4,3(i2.2))
c
strmid(1:4)=natcf(1:4)
write(strmid(5:8),601) iatcfyr4
601 format(i4.4)
c
c Read the track model name from the track file
read(luit,108) tem4
108 format(24x,a4)
tmname = tem4
c
c Determine the basin from a combination of ATCF ID and initial position.
c AL and SH cases do not change basins but EP, CP, WP and IO can. Decide
c between EP, CP, WP and IO using initial position.
glat0 = 0.1*float(ilat0)
glon0 = 0.1*float(ilon0)
if (glon0 .lt. 0) glon0 = glon0 + 360.0
c
if (natcf(1:1) .eq. 'A') then
ibasin=1
elseif (natcf(1:1) .eq. 'S') then
ibasin=5
elseif (natcf(1:1) .eq. 'I' .or. natcf(1:1) .eq. 'W' .or.
+ natcf(1:1) .eq. 'C' .or. natcf(1:1) .eq. 'E') then
if (glon0 .ge. 15.0 .and. glon0 .lt. 105.0) then
ibasin = 4
elseif (glon0 .ge. 105.0 .and. glon0 .lt. 180.0) then
ibasin = 3
else
ibasin = 2
endif
else
write(lulg,911) natcf
911 format('Unrecognized basin in iships, natcf=',a8)
stop
endif
c
c Set basin dependent parameters
c fnco = Coefficient file name
c jdmax = Julian Day of peak hurricane season
c ishrx = The maximum no. of forcast times to use the
c shear (1=12 hr, 2=24 hr, 3=36 hr, etc)
c irefcx= Similar to ishrx for refc
c it200x= Similar to ishrx for t200
c it250x= Similar to ishrx for t250
c t250th= Threshold value for t250 predictor
c icione= flag for Joe Cione SST cooling algorithm
c 0 to skip cooling algorithm
c 1 to include it with initial value of lat and storm speed
c 2 to include it with forecast lat and storm speed values
c mpispd= flag to adjust MPI for storm translational speed
c 0 to use mean speed adjusted mpi routine (mpical)
c 1 to use original mpi routine (mpicalo) and add
c fraction of storm speed from Schewerdt formula a = 1.5*c**0.63
c iohcsm= Sample mean value of OHC
c rthresh=OHC threshold
c rsdir = Reference direction for shear direction predictor (sdp)
c th1,th2,f1,f2 = factors for scaling sdp by latitude
c sna,snb,snc,snd: Coefficients for log transformation of shear
c ga0,ga1,ga2: Coefficients for mean SHGC from SHDC
c
if (ibasin .eq. 5) then
c Southern Hemisphere parameters
if (ioper .eq. 1) then
fnco = trim( coef_location ) // 'ships13_coef_shem.dat'
else
fnco = 'ships13_coef_shem.dat'
fnco(6:7) = coyear
endif
sgnlon = 1.0
sgnlat = -1.0
jdmax = 31
raefld = 30.0
ishrx =mft
irefcx=mft
it200x=mft
it250x=mft
t250th= 0.0
c start icione ------------
icione=0
c
ships_sst_params%icione_sst = 0
ships_sst_params%icione_vsst = 0
c
lgem_sst_params%icione_sst = 0
lgem_sst_params%icione_vsst = 0
c
ric2o_sst_params%icione_sst = 0
ric2o_sst_params%icione_vsst = 0
c
etcl_sst_params%icione_sst = 0
etcl_sst_params%icione_vsst = 0
c
kaplan_sst_params%icione_sst = 0
kaplan_sst_params%icione_vsst = 0
c
rozof_sst_params%icione_sst = 0
rozof_sst_params%icione_vsst = 0
c
prsef_sst_params%icione_sst = 0
prsef_sst_params%icione_vsst = 0
ahi_sst_params%icione_sst = 0
ahi_sst_params%icione_vsst = 0
c end icione ------------
mpispd=1
iohcsm=20
rthresh= 25.0
rsdir = 120.0
th1 = -10.0
th2 = -30.0
f1 = 1.0
f2 = -0.75
c
ga0 = 10.12
ga1 = 0.8348
ga2 = 0.0041
c
sna = 20.0
snb = 0.5
snc = 5.0
snd = 33.5
potmin = 20.0
elseif (ibasin .eq. 4) then
c North Indian Ocean parameters
if (ioper .eq. 1) then
fnco = trim( coef_location ) // 'ships13_coef_iocn.dat'
else
fnco = 'ships13_coef_iocn.dat'
fnco(6:7) = coyear
endif
sgnlon = 1.0
sgnlat = 1.0
jdmax =180
raefld = 60.0
ishrx =mft
irefcx=mft
it200x=mft
it250x=mft
t250th= -39.0
c start icione ------------
icione=0
c
ships_sst_params%icione_sst = 0
ships_sst_params%icione_vsst = 0
c
lgem_sst_params%icione_sst = 0
lgem_sst_params%icione_vsst = 0
c
ric2o_sst_params%icione_sst = 0
ric2o_sst_params%icione_vsst = 0
c
etcl_sst_params%icione_sst = 0
etcl_sst_params%icione_vsst = 0
c
kaplan_sst_params%icione_sst = 0
kaplan_sst_params%icione_vsst = 0
c
rozof_sst_params%icione_sst = 0
rozof_sst_params%icione_vsst = 0
c
prsef_sst_params%icione_sst = 0
prsef_sst_params%icione_vsst = 0
c
ahi_sst_params%icione_sst = 0
ahi_sst_params%icione_vsst = 0
c end icione ------------
mpispd=1
iohcsm=55
rthresh= 0.0
rsdir = 60.0
th1 = 10.0
th2 = 20.0
f1 = 1.0
f2 = 0.25
c
ga0 = 10.17
ga1 = 0.9933
ga2 = 0.0035
c
sna = 20.0
snb = 0.5
snc = 5.0
snd = 33.5
potmin = 20.0
elseif (ibasin .eq. 3) then
c Western North Pacific parameters
if (ioper .eq. 1) then
fnco = trim( coef_location ) // 'ships13_coef_wpac.dat'
else
fnco = 'ships13_coef_wpac.dat'
fnco(6:7) = coyear
endif
sgnlon = 1.0
sgnlat = 1.0
jdmax =248
raefld = 30.0
ishrx =mft
irefcx=mft
it200x=mft
it250x=mft
t250th= -39.0
c start icione ------------
icione=0
c
ships_sst_params%icione_sst = 0
ships_sst_params%icione_vsst = 0
c
lgem_sst_params%icione_sst = 0
lgem_sst_params%icione_vsst = 0
c
ric2o_sst_params%icione_sst = 0
ric2o_sst_params%icione_vsst = 0
c
etcl_sst_params%icione_sst = 0
etcl_sst_params%icione_vsst = 0
c
kaplan_sst_params%icione_sst = 0
kaplan_sst_params%icione_vsst = 0
c
rozof_sst_params%icione_sst = 0
rozof_sst_params%icione_vsst = 0
c
prsef_sst_params%icione_sst = 0
prsef_sst_params%icione_vsst = 0
c
ahi_sst_params%icione_sst = 0
ahi_sst_params%icione_vsst = 0
c end icione ------------
mpispd=1
iohcsm=45
rthresh= 0.0
rsdir = 60.0
th1 = 10.0
th2 = 20.0
f1 = 1.0
f2 = 0.25
c
ga0 = 10.17
ga1 = 0.9933
ga2 = 0.0035
c
sna = 20.0
snb = 0.5
snc = 5.0
snd = 33.5
potmin = 20.0
elseif (ibasin .eq. 2) then
c East/Central North Pacific parameters
fnco = trim( coef_location ) // 'ships24_coef_epac.dat'
sgnlon = -1.0
sgnlat = 1.0
jdmax =238
raefld = 30.0
ishrx =mft
irefcx=mft
it200x=mft
it250x=mft
t250th= 0.0
c start icione ------------
icione=0
c
ships_sst_params%icione_sst = 0
ships_sst_params%icione_vsst = 0
c
lgem_sst_params%icione_sst = 0
lgem_sst_params%icione_vsst = 0
c
ric2o_sst_params%icione_sst = 0
ric2o_sst_params%icione_vsst = 0
c
etcl_sst_params%icione_sst = 0
etcl_sst_params%icione_vsst = 0
c
kaplan_sst_params%icione_sst = 0
kaplan_sst_params%icione_vsst = 0
c
rozof_sst_params%icione_sst = 0
rozof_sst_params%icione_vsst = 0
c
prsef_sst_params%icione_sst = 0
prsef_sst_params%icione_vsst = 0
c
ahi_sst_params%icione_sst = 0
ahi_sst_params%icione_vsst = 0
c end icione ------------
mpispd=1
iohcsm=35
rthresh=0.0
c
rsdir = 70.0
th1 = 10.0
th2 = 30.0
f1 = 1.0
f2 = -1.0
c
ga0 = 12.74
ga1 = 0.6572
ga2 = 0.0070
c
sna = 20.0
snb = 0.5
snc = 5.0
snd = 33.5
potmin = 20.0
else
c Atlantic parameters
fnco = trim( coef_location ) // 'ships24_coef_atlc.dat'
sgnlon = -1.0
sgnlat = 1.0
jdmax =253
raefld = 30.0
ishrx =mft
irefcx=mft
it200x=mft
it250x=mft
t250th= -44.0
c start icione ------------
icione=0
c
ships_sst_params%icione_sst = 0
ships_sst_params%icione_vsst = 0
c
lgem_sst_params%icione_sst = 0
lgem_sst_params%icione_vsst = 0
c
ric2o_sst_params%icione_sst = 0
ric2o_sst_params%icione_vsst = 0
c
etcl_sst_params%icione_sst = 2
etcl_sst_params%icione_vsst = 2
c
kaplan_sst_params%icione_sst = 2
kaplan_sst_params%icione_vsst = 2
c
c !TODO - shoudl that be 0 or 2 ??? Need to aks Chris
rozof_sst_params%icione_sst = 0
rozof_sst_params%icione_vsst = 2
c
prsef_sst_params%icione_sst = 2
prsef_sst_params%icione_vsst = 2
c
ahi_sst_params%icione_sst = 0
ahi_sst_params%icione_vsst = 0
c end icione ------------
mpispd=1
iohcsm=45
rthresh=60.0
rsdir = 55.0
th1 = 10.0
th2 = 40.0
f1 = 1.0
f2 = -0.75
c
ga0 = 10.12
ga1 = 0.8348
ga2 = 0.0041
c
sna = 20.0
snb = 0.5
snc = 5.0
snd = 33.5
potmin = 20.0
endif
c
isgnlon = nint(sgnlon)
c
c Initialize predictors to missing values
do k=-2,mft
ilat(k) = imiss
ilon(k) = imiss
enddo
c
do k=0,mft
ishr(k) = imiss
ishtd(k) = imiss
c iohcc(k) = imiss
c ibohc(k) = imiss
c irhcn(k) = imiss
igoes(k) = imiss
igoesm1(k)= imiss
igoesm3(k)= imiss
igoesxx(k)= imiss
ipw01(k) = imiss
ipw03(k) = imiss
ipw04(k) = imiss
ipw06(k) = imiss
ipw14(k) = imiss
ipw19(k) = imiss
sst_lgem(k) = rmiss
vsst_lgem(k) = rmiss
sst_ships(k) = rmiss
iohc_inp_ncoda(k) = imiss
iohc_climo_ncoda(k) = imiss
iohc_best_ncoda(k) = imiss
ohc_final_ncoda(k) = rmiss
iohc_ships(k) = imiss
iohc_lgem(k) = imiss
iohc_ric2(k) = imiss
iohc_kaplan(k) = imiss
iohc_rozof(k) = imiss
iohc_prse(k) = imiss
iohc_etcl(k) = imiss
idist(k) = imiss
iu200(k) = imiss
it150(k) = imiss
it200(k) = imiss
it250(k) = imiss
ig200(k) = imiss
iz850(k) = imiss
id200(k) = imiss
itadv(k) = imiss
itgrd(k) = imiss
iv500(k) = imiss
iv300(k) = imiss
iv20c(k) = imiss
iepos(k) = imiss
irefc(k) = imiss
irhlo(k) = imiss
irhmd(k) = imiss
irhhi(k) = imiss
ipslv(k) = imiss
ishdc(k) = imiss
ishgc(k) = imiss
ivmaxb(k)= imiss
isddc(k) = imiss
itwac(k) = imiss
icflx(k) = imiss
iaipc(k) = imiss
istpw(k) = imiss
ipc00(k) = imiss
ipcm1(k) = imiss
ipcm3(k) = imiss
c
imtpw(k) = imiss
iusetpw(k) = imiss
idivc(k) = imiss
ishrg(k) = imiss
iepss(k) = imiss
ienss(k) = imiss
ishrd(k) = imiss
idipc(k) = imiss
ishrs(k) = imiss
ie000(k) = imiss
ieneg(k) = imiss
enddo
c
print*, 'WILL NOW READ LSDIAG >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
c ++ Begin reading the lsdiag.dat file
read(luls,110) sname4,ilyr,ilmon,ilday,iltime,ivmx
print*, 'LSDIAG HEADER >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
print*, 'LSDIAG HEADER >>>>>',sname4,ilyr,ilmon,ilday,iltime,ivmx
print*, 'LSDIAG HEADER >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>'
110 format(1x,a4,1x,3(i2),1x,i2,1x,i4)
c
c Set max number of lines in an individual record of the lsdiag.dat
mxl = 1100
lsloop: do i=1,mxl
read(luls,111) iline
111 format(a180)
c
if (iline(157:160) .eq. 'LAST' ) then
exit lsloop
endif
c
if (iline(157:160) .eq. 'DELV' )then
read(iline,112) iper
112 format(1x,i4)
endif
c
if (iline(157:160) .eq. 'LAT ' )then
read(iline,114) (ilat(k),k=-2,mft)
114 format( 1x,31(i4,1x))
endif
c
if (iline(157:160) .eq. 'LON ' )then
read(iline,115) (ilon(k),k=-2,mft),loncon
115 format( 1x,31(i4,1x),8x,a1)
c
c Convert to 0 to 360 convention for now
if (loncon .eq. 'W') then
do k=-2,mft
if (ilon(k) .lt. imiss) then
ilon(k) = 3600 - ilon(k)
endif
enddo
endif
endif
c
if (iline(157:160) .eq. 'SHRD') then
read(iline,116) (ishr(k),k= 0,mft)
116 format(11x,29(i4,1x))
endif
c
if (iline(157:160) .eq. 'SHRS') then
read(iline,116) (ishrs(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'SDDC') then
read(iline,116) (isddc(k),k= 0,mft)
c
c Convert ishtd from heading to direction
do k=0,mft
if (isddc(k) .lt. imiss) then
ishtd(k) = 180 + isddc(k)
if (ishtd(k) .gt. 360) ishtd(k) = ishtd(k)-360
endif
enddo
endif
c
c---------------------------------------------
c if (iline(117:120) .eq. 'XOHC') then
c read(iline,116) (iohcc(k),k= 0,mft)
c endif
c
c if (iline(117:120) .eq. 'NOHC') then
c read(iline,117) (irhcn(k),k= 0,mft),iohclag
c endif
c
if (iline(157:160) .eq. 'XOHC') then
read(iline,116) (iohc_climo_ncoda(k),k= 0,mft)
endif
if (iline(157:160) .eq. 'NOHC') then
read(iline,117) (iohc_inp_ncoda(k),k= 0,mft),iohc_lag_days
endif
c---------------------------------------------
c
if (iline(157:160) .eq. 'IR00') then
read(iline,116) (igoes(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'IRM1') then
read(iline,116) (igoesm1(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'IRM3') then
read(iline,116) (igoesm3(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'IRXX') then
read(iline,116) (igoesxx(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PW01') then
read(iline,116) (ipw01(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PW03') then
read(iline,116) (ipw03(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PW04') then
read(iline,116) (ipw04(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PW06') then
read(iline,116) (ipw06(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PW14') then
read(iline,116) (ipw14(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PW19') then
read(iline,116) (ipw19(k),k= 0,mft)
endif
c ------------CLIMO SST --BEGIN ---------------------------
c wsst Reynolds
if (iline(157:160) .eq. 'CSST') then
read(iline,116) (wsst_set(k)%isst_climo,k= 0,mft)
endif
c
c daily Reynolds
c read XDST into 2 climo variables: for dsst and for dsta
if (iline(157:160) .eq. 'XDST') then
c for nsst climo
read(iline,116) (dsst_set(k)%isst_climo,k= 0,mft)
c for dsta climo
do k = 0,mft
dsta_set(k)%isst_climo = dsst_set(k)%isst_climo
enddo
endif
c
c daily NCODA
c read XNST into 2 climo variables: for nsst and for nsta
if (iline(157:160) .eq. 'XNST') then
c for nsst climo
read(iline,116) (nsst_set(k)%isst_climo,k= 0,mft)
c for nsta climo
do k = 0,mft
nsta_set(k)%isst_climo = nsst_set(k)%isst_climo
enddo
endif
c
c ------------DATA SST --BEGIN ---------------------------
c special format for reading SST, OHC variables that have
c extra number for age in days
117 format(11x,29(i4,1x),5x,i4)
cc
c wsst Reynolds
if (iline(157:160) .eq. 'RSST') then
read(iline,117) (wsst_set(k)%isst_inp,k= 0,mft),
+ wsst_params%isst_lag
endif
c
c daily Reynolds DSST
if (iline(157:160) .eq. 'DSST') then
read(iline,117) (dsst_set(k)%isst_inp,k= 0,mft),
+ dsst_params%isst_lag
endif
c
c daily NCODA NSST
if (iline(157:160) .eq. 'NSST') then
read(iline,117) (nsst_set(k)%isst_inp,k= 0,mft),
+ nsst_params%isst_lag
endif
c
c daily Reynolds DSTA
if (iline(157:160) .eq. 'DSTA') then
read(iline,117) (dsta_set(k)%isst_inp,k= 0,mft),
+ dsta_params%isst_lag
endif
c
c daily NCODA NSTA
if (iline(157:160) .eq. 'NSTA') then
read(iline,117) (nsta_set(k)%isst_inp,k= 0,mft),
+ nsta_params%isst_lag
endif
c
c ------------DATA SST --END ---------------------------
c
if (iline(157:160) .eq. 'DTL') then
read(iline,116) (idist(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'U200') then
read(iline,116) (iu200(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'T150') then
read(iline,116) (it150(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'T200') then
read(iline,116) (it200(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'T250') then
read(iline,116) (it250(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'G200') then
read(iline,116) (ig200(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'Z850') then
read(iline,116) (iz850(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'D200') then
read(iline,116) (id200(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'TADV') then
read(iline,116) (itadv(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'TGRD') then
read(iline,116) (itgrd(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'V500') then
read(iline,116) (iv500(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'V300') then
read(iline,116) (iv300(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'V20C') then
read(iline,116) (iv20c(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'E000') then
read(iline,116) (ie000(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'EPOS') then
read(iline,116) (iepos(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'ENEG') then
read(iline,116) (ieneg(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'REFC') then
read(iline,116) (irefc(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'RHLO') then
read(iline,116) (irhlo(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'RHMD') then
read(iline,116) (irhmd(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'RHHI') then
read(iline,116) (irhhi(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'PSLV') then
read(iline,116) (ipslv(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'SHDC') then
read(iline,116) (ishdc(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'SHGC') then
read(iline,116) (ishgc(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'VMAX') then
read(iline,116) (ivmaxb(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'TWAC') then
read(iline,116) (itwac(k),k= 0,mft)
endif
if (itwac(0) .eq. imiss) itwac(0) = 0
c
c PARMOD+ For TIC and LI RII
c
if (iline(157:160) .eq. 'CFLX') then
read(iline,116) (icflx(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'AIPC') then
read(iline,116) (iaipc(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'STPW') then
read(iline,116) (istpw(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'PC00') then
read(iline,116) (ipc00(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'PCM1') then
read(iline,116) (ipcm1(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'PCM3') then
read(iline,116) (ipcm3(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'MTPW') then
read(iline,116) (imtpw(k),k=0,mft)
endif
c PARMOD-
c
c PARMOD+ For MTC RII
if (iline(157:160) .eq. 'DIVC') then
read(iline,116) (idivc(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'SHRG') then
read(iline,116) (ishrg(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'EPSS') then
read(iline,116) (iepss(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'ENSS') then
read(iline,116) (ienss(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'SHRD') then
read(iline,116) (ishrd(k),k=0,mft)
endif
c
if (iline(157:160) .eq. 'DIPC') then
read(iline,116) (idipc(k),k=0,mft)
endif
c
c*************** start add code, PrSEFoNe (1 FEB 2010, jpk)
if (iline(157:160) .eq. 'PENC') then
read(iline,116) (ipenc(k),k= 0,mft)
endif
c
if (iline(157:160) .eq. 'VMPI') then
read(iline,116) (ivmpi2(k),k= 0,mft)
endif
c*************** end add code, PrSEFoNe
c
enddo lsloop
close(luls)
print*, 'DONE READING LSDIAG <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<,'
c
c ++ End reading the lsdiag.dat file
c
do k=0,mft
if (ivmaxb(k) .ge. imiss) ivmaxb(k) = 0
enddo
ccccccccccccccccccccccccccc START select best OHC ccccccccccccccccccccc
c Check OHC data. If missing or older than iohclagx,
c replace with PHCN.
c If that is missing, use XOHC
c call ohcclimbyobs(irhcn,iohcc,iohclag,
c + ipohc,mft,imiss,ibohc)
call select_best_ohc(iohc_max_lag_days,
+ imiss,
+ iohc_inp_ncoda,iohc_climo_ncoda,
+ iohc_lag_days,mft,
+ iohc_best_ncoda, iohc_best_ncoda_lag, ipohc)
c do k=0,mft
c print*, "DEBUG OHC2 ", (k), iohc_inp_ncoda(k),
c + iohc_climo_ncoda(k), iohc_best_ncoda(k)
c enddo
ccccccccccccccccccccccccccc END select best OHC ccccccccccccccccccccc
c
c ++ Begin GOES data quality control
c
vmx = float(ivmx)
c
if (igoes(4) .lt. imiss) then
c Check GOES00 standard deviation variable to see if
c it unreasonably large (+/- 3sigma)
stdck = abs(vmx*0.1*float(igoes(4))-1000.0)/450.0
c
write(lulg,824) vmx,igoes(4),stdck
824 format(/,' t=0 GOES std chk: vmx, igoes4, stdck=',
+ f5.0,1x,i4,1x,f6.2)
c
if (stdck .gt. 3.0) then
do ii=0,mft
igoes(ii) = imiss
enddo
endif
endif
c
if (igoesm3(4) .lt. imiss) then
c Check GOESM3 standard deviation variable to see if
c it unreasonably large (+/- 3sigma)
stdck = abs(vmx*0.1*float(igoesm3(4))-1000.0)/450.0
c
write(lulg,825) vmx,igoes(4),stdck
825 format(/,' t=-3 GOES std chk: vmx, igoes4, stdck=',
+ f5.0,1x,i4,1x,f6.2)
c
if (stdck .gt. 3.0) then
do ii=0,mft
igoes(ii) = imiss
enddo
endif
endif
c
c Perform buddy check of GOES data if possible
it1 = igoes(1)
is1 = igoes(2)
it2 = igoesm3(1)
is2 = igoesm3(2)
c
if (it1 .lt. imiss .and. it2 .lt. imiss .and.
+ is1 .lt. imiss .and. is2 .lt. imiss ) then
adt = 0.1*abs(float(it2-it1))
ads = 0.1*abs(float(is2-is1))
c
write(lulg,826) it1,it2,is1,is2,adt,ads
826 format(/,' GOES buddy check input: ',4(i4,1x),2(f6.1,1x))
c
c If dt or ds too large, one of the GOES profiles is bad
if (adt .gt. 40.0 .or. ads .gt. 20.0) then
do ii=0,mft
igoes(ii) = imiss
enddo
endif
endif
c
c If t=0 GOES data missing, replace with t=-1.5 data.
c If still missing, replace with t=-1.5 data.
c If still missing, replace with proxy GOES.
c If still missing, replace with sample mean values
ipgoes = 0
do k=0,ngoes
if (igoes(k) .ge. imiss .and. igoesm1(k) .lt. imiss) then
igoes(k) = igoesm1(k)
endif
enddo
c
do k=0,ngoes
if (igoes(k) .ge. imiss .and. igoesm3(k) .lt. imiss) then
igoes(k) = igoesm3(k)
endif
enddo
c
do k=0,ngoes
if (igoes(k) .ge. imiss .and. igoesxx(k) .lt. imiss) then
igoes(k) = igoesxx(k)
ipgoes=1
endif
enddo
c
do k=0,ngoes
if (igoes(k) .ge. imiss) then
igoes(k) = igoessm(k)
ipgoes=2
endif
enddo
c
c ++ End GOES data quality control
c
c Check to make sure at least the initial value of each atmospheric
c predictor is available
istop=0
if (ishdc(0) .eq. imiss) istop=1
if (ishgc(0) .eq. imiss) istop=1
if (isddc(0) .eq. imiss) istop=1
if (it200(0) .eq. imiss) istop=1
if (it250(0) .eq. imiss) istop=1
if (ig200(0) .eq. imiss) istop=1
if (iepos(0) .eq. imiss) istop=1
if (irhmd(0) .eq. imiss) istop=1
if (id200(0) .eq. imiss) istop=1
if (itadv(0) .eq. imiss) istop=1
if (iz850(0) .eq. imiss) istop=1
if (ipw06(0) .eq. imiss) istop=1
c
if (istop .eq. 1) then
write(lulg,990)
990 format(/' Some t=0 predictors missing')
stop '***** SHIPS HALTED DUE TO INPUT DATA ERROR *****'
endif
c
c Convert t200, t250 and g200 to real, and apply 250 threshold
do k=0,mft
if (it200(k) .lt. imiss) then
t200(k) = 0.1*float(it200(k))
else
t200(k) = rmiss
endif
c
if (ig200(k) .lt. imiss) then
g200(k) = 0.1*float(ig200(k))
else
g200(k) = rmiss
endif
c
if (it250(k) .lt. imiss) then
t250t = 0.1*float(it250(k))
t250t = t250t - t250th
if (t250t .gt. 0.0) t250t = 0.0
t250(k) = t250t
else
t250(k) = rmiss
endif
enddo
c
c Convert time dependent TPW variables to reals
do k=0,mft
if (ipw06(k) .lt. imiss) then
pw06(k) = float(ipw06(k))
else
pw06(k) = rmiss
endif
enddo
c
c Convert twac to real and fill out array
do k=0,mft
if (itwac(k) .eq. imiss) then
twac(k) = rmiss
else
twac(k) = 0.1*float(itwac(k))
endif
enddo
call patch0(twac,rmiss,mft)
cccccccccc START select best SST ccccccccccccccccccccccccccccccccccc
c Replace the climatological SST with the Reynolds SST
c data if it is available
c select best available sst for each input sst
c
c --START EXAMPLE -------------------
c Reynolds weekly WSST - WSST - WSST - CSST
c subroutine select_best_sst will select best sst
c for the 2nd argument
c for example:
c
c call select_best_sst(isst_max_lag_days,
c +nsta_set,nsta_params,
c +dsta_set,dsta_params,
c +wsst_set,wsst_params)
c
c - will select best sst for NSTA
c - will use NSTA input sst (read from lsdiag) if valid
c - will next use DSTA input sst (read from lsdiag) if valid
c - will next use WSST input sst (read from lsdiag) if valid
c - will next use NSTA climo sst - climo from 2nd argument
c
c --END EXAMPLE -------------------
c Reynolds weekly WSST - WSST - WSST - CSST
call select_best_sst(mft, isst_max_lag_days,
+wsst_set,wsst_params,
+wsst_set,wsst_params,
+wsst_set,wsst_params,lulg )
c Reynolds daily - DSST - DSST - WSST - XDST
call select_best_sst(mft, isst_max_lag_days,
+dsst_set,dsst_params,
+dsst_set,dsst_params,
+wsst_set,wsst_params,lulg )
c Reynolds daily spatial average - DSTA - DSTA - WSST - XDST
call select_best_sst(mft, isst_max_lag_days,
+dsta_set,dsta_params,
+dsta_set,dsta_params,
+wsst_set,wsst_params,lulg )
c NCODA daily - NSST - DSST - WSST - XNST
call select_best_sst(mft, isst_max_lag_days,
+nsst_set,nsta_params,
+dsta_set,dsta_params,
+wsst_set,wsst_params,lulg )
c
c NCODA daily spatial average - NSTA - DSTA - WSST - XNST
call select_best_sst(mft, isst_max_lag_days,
+nsta_set,nsta_params,
+dsta_set,dsta_params,
+wsst_set,wsst_params,lulg )
c
cccccccccc END select best SST ccccccccccccccccccccccccccccccccccc
c ===============================================================
c Calculate shear direction predictor
do k=0,mft
if (ilat(k) .lt. imiss .and. isddc(k) .lt. imiss) then
tlat = 0.1*float(ilat(k))
sddct= float(isddc(k))
c
call sdpar(sddct,rsdir,tlat,th1,th2,f1,f2,sdir(k))
c
c adif = abs(rsdir-sddct)
c if (adif .gt. 180.0) adif = abs(360.0-adif)
c if (adif .gt. 180.0) adif = 0.0
c adif = adif - 90.0
c
c if (tlat .le. th1) then
c fscale = f1
c elseif (tlat .ge. th2) then
c fscale = f2
c else
c z = (tlat-th1)/(th2-th1)
c fscale = f1 + 3.0*(f2-f1)*z*z - 2.0*(f2-f1)*z*z*z
c endif
c
c sdir(k) = fscale*adif
else
sdir(k) = 9999.
endif
c write(6,731) k,tlat,rsdir,sddct,adif,sdir(k)
c 731 format('k,lat,rsdir,sddc,adif,sdir: ',i2,6(1x,f7.1))
enddo
c
c Convert shgc to a deviation from value based on shdc and calculate
c log-transformed shdc if ilogt=1
ilogt = 0
do k=0,mft
if (ishdc(k) .lt. imiss .and. ishgc(k) .lt. imiss) then
sdtem = 0.1*float(ishdc(k))
sgtem = 0.1*float(ishgc(k))
c
if (ilogt .eq. 1) then
shdct(k) = sna*alog(snb*sdtem + snc) - snd
else
shdct(k) = sdtem
endif
sgbar = ga0 + ga1*sdtem + ga2*sdtem*sdtem
shgc(k) = sgtem - sgbar
else
shdct(k) = rmiss
shgc(k) = rmiss
endif
enddo
c
c Read the SHIPS coefficient file
print*, 'WILL NOW OPEN COEFF FILE ======================', fnco
open(unit=luco,file=fnco,form='formatted',status='old',err=900)
print*, 'OPENED SHIPS COEFF FILE ======================'
c
do 5 k=1,mft
read(luco,135) ybar(k),ysig(k)
135 format(14x,e12.5,1x,e12.5)
c
do 6 n=1,nvar
read(luco,136) coef(n,k),xbar(n,k),xsig(n,k)
136 format(1x,3(e12.5,1x))
6 continue
5 continue
close(luco)
c
c Write basic storm information to the log file
write(lulg,302) sname,isyr,ismon,isday,istime
302 format(' SHIPS forecast for ',a10,3(i2.2),1x,i2.2)
c
write(lulg,303) iobslag
303 format(/,' Reynolds SST is ',i4,' hrs old')
c
do 12 k=-2,mft
if (ilat(k) .ge. imiss) then
rlat(k) = rmiss
else
rlat(k) = 0.1*float(ilat(k))
endif
c
if (ilon(k) .ge. imiss) then
rlon(k) = rmiss
else
rlon(k) = 0.1*float(ilon(k))
if (sgnlon .lt. 0.0) then
rlon(k) = rlon(k) - 360.0
c rlon(k) = sgnlon*rlon(k)
endif
endif
12 continue
c
c Calculate the storm translational speed along the forecast track
do k=-2,mft
ftimec(k) = delt*float(k)
enddo
call tspdcale(rlat,rlon,ftimec,mft,rmiss,cxt,cyt,cmagt)
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c----BEGIN------------------GET MPI and cooled SST -----------------
c for each SST get 4 new variables:
c sst - float veriosn of SST
c sstc - sst cooled with icione
c vsst - MPI based on SST
c vsstc - MPI based on cooled SST
c EXAMPLE call:
c input - ISST
c output - SST,SSTC,VSST,VSSTC
c ----------------------------------
c ONLY use wsst_set in the call - the subroutine
c will read wsst_set%isst_inp
c and calculate
c wsst_set%sst_final,
c wsst_set%sstc_final,
c wsst_set%vsst_final,
c wsst_set%vsstc_final,
c ----------------------------------
c call get_cooling_mpi(mft,rlat,ispeed,cmagt,
c + mpispd,ibasin,rmiss,ivmpi2,
c + wsst_set )
c WSST
call get_cooling_mpi(mft,rlat,ispeed,cmagt,
+ mpispd,ibasin,rmiss,ivmpi2,
+ wsst_set )
c Reynolds daily
call get_cooling_mpi(mft,rlat,ispeed,cmagt,
+ mpispd,ibasin,rmiss,ivmpi2,
+ dsst_set )
c Reynolds daily spatial average
call get_cooling_mpi(mft,rlat,ispeed,cmagt,
+ mpispd,ibasin,rmiss,ivmpi2,
+ dsta_set )
c NCODA_daily
call get_cooling_mpi(mft,rlat,ispeed,cmagt,
+ mpispd,ibasin,rmiss,ivmpi2,
+ nsst_set )
c NCODA_daily spatial average
call get_cooling_mpi(mft,rlat,ispeed,cmagt,
+ mpispd,ibasin,rmiss,ivmpi2,
+ nsta_set )
c----END------------------GET MPI and cooled SST -----------------
c ===============================================================
c ===============================================================
c ----get OHC for each application
do k=0,mft
ohc_final_ncoda(k) = float(iohc_best_ncoda(k))
enddo
c =======================================
call get_model_ohc(mft,
+ cships_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_ships,iohc_sm_ships,ohc_thresh_ships,lulg)
c do k=0,mft
c print*, "DEBUG OHC1 ", iohc_ships(k),iohc_best_ncoda(k)
c enddo
call get_model_ohc(mft,
+ clgem_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_lgem,iohc_sm_lgem,ohc_thresh_lgem,lulg)
c
call get_model_ohc(mft,
+ cric2_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_ric2,iohc_sm_ric2,ohc_thresh_ric2,lulg)
call get_model_ohc(mft,
+ ckaplan_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_kaplan,iohc_sm_kaplan,ohc_thresh_kaplan,lulg)
call get_model_ohc(mft,
+ crozof_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_rozof,iohc_sm_rozof,ohc_thresh_rozof,lulg)
call get_model_ohc(mft,
+ cetcl_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_etcl,iohc_sm_etcl,ohc_thresh_etcl,lulg)
call get_model_ohc(mft,
+ cprse_ohc_type,
+ rthresh,
+ iohcsm,
+ iohc_best_ncoda,
+ iohc_prse,iohc_sm_prse,ohc_thresh_prse,lulg)
c ===============================================================
c ===============================================================
c ===============================================================
c
c Calculate non-time dependent predictors and
c other related variables
vmx = float(ivmx)
v(0) = vmx
c
if (iper .lt. 90) then
per = float(iper)
else
per = 0.0
endif
c
c Check to see if storm is now over the ocean but was
c recently over land. If so, modify the persistence
c accordingly.
call lcheck(ilat0,ilon0,ilatm12,ilonm12,12.0,tland,lmod)
if (lmod .eq. 1 .and. per .lt. 0.0) per = 0.0
c
vper = per*vmx
pslv = float(ipslv(0))
d200 = float(id200(0))
c
c Calculate the eastward component of the storm motion at t=0
speed = float(ispeed)
head = float(ihead)
call rhtoc(speed,head,spdx,spdy)
c
c Calculate the Julian day predictor
call jday(ismon,isday,isyr,julday)
aday = abs(float(julday-jdmax))
if (aday .gt. 182.5) aday = 365.0-aday
c
c Calculate new ADAY parameter
raday = (aday/raefld)
aday = exp(-raday*raday)
c
c GOES predictors
pc20 = float(igoes(6))
pc20m = float(igoessm(6))
btstd = 0.1*float(igoes(2))
btstdm = 0.1*float(igoessm(2))
gstd = vmx*btstd
c
ccc c ----BEGIN-------SELECT SST FOR each part of the code -------
c SHIPS
call get_model_sst(mft,
+ ships_sst_set,
+ ships_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
c LGEM
call get_model_sst(mft,
+ lgem_sst_set,
+ lgem_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
c For lgem, use additional variables to not mess with common
do k = 0,mft
sst_lgem(k) = lgem_sst_set(k)%sst_final
vsst_lgem(k) = lgem_sst_set(k)%vsst_final
enddo
c For ships, extract final SST
do k=0,mft
sst_ships(k) = ships_sst_set(k)%sst_final
enddo
c RIC2O
call get_model_sst(mft,
+ ric2o_sst_set,
+ ric2o_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
c ETCL
call get_model_sst(mft,
+ etcl_sst_set,
+ etcl_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
c KAPLAN
call get_model_sst(mft,
+ kaplan_sst_set,
+ kaplan_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
c ROZOF
call get_model_sst(mft,
+ rozof_sst_set,
+ rozof_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
c AHI
call get_model_sst(mft,
+ ahi_sst_set,
+ ahi_sst_params,
+ wsst_set,
+ dsst_set,
+ dsta_set,
+ nsst_set,
+ nsta_set,lulg)
ccc c outputs single verions for sst and vsst for each
ccc c application - cooled or non-cooled depending on 2 icione
ccc c flags for each application
c --- END-------SELECT SST FOR each part of the code -------
c ++ Begin variables for input to the rapid intensity index
perri = per
c
if (ipgoes .eq. 0) then
pcri30 = float(igoes(7))
sbtri = 0.1*float(igoes(2))
sir2ri = 0.1 * float(igoes(1))
sir4ri = 0.1 * float(igoes(3))
sir5ri = 0.1 * float(igoes(4))
sir6ri = float(igoes(5))
sir7ri = float(igoes(6))
sir9ri = float(igoes(8))
sir10ri = float(igoes(9))
sir11ri = float(igoes(10))
sir12ri = 0.1 * float(igoes(11))
sir13ri = 0.1 * float(igoes(12))
sir14ri = float(igoes(13))
else
pcri30 = 999.
sbtri = 999.
sir2ri = 999.
sir4ri = 999.
sir5ri = 999.
sir6ri = 999.
sir7ri = 999.
sir9ri = 999.
sir10ri = 999.
sir11ri = 999.
sir12ri = 999.
sir13ri = 999.
sir14ri = 999.
endif
c
shdcri = 0.1*float(ishdc(0))
rhlori = float(irhlo(0))
d200ri = d200
c
c this is for RII
vssttm = kaplan_sst_set(0)%vsst_final
potri = vssttm-vmx
rhcri = float(iohc_kaplan(0))
c
c Perform time average of RI variables if necessary
ic2=1
ic3=1
ic4=1
ic5=1
ic6=1
ic7=1
ic8=1
do i=1,4
if (ishdc(i) .lt. 9999) then
ic7 = ic7+1
shdcri = shdcri + 0.1*float(ishdc(i))
endif
c
vssttm = kaplan_sst_set(0)%vsst_final
c
if (vssttm .lt. 200.0) then
ic3 = ic3+1
ic4 = ic4+1
sstri = sstri + kaplan_sst_set(i)%vsst_final
potri = potri + (vssttm-vmx)
endif
c
if (irhlo(i) .lt. 9999) then
ic5 = ic5+1
rhlori = rhlori + float(irhlo(i))
endif
c
if (id200(i) .lt. 9999) then
ic6 = ic6+1
d200ri = d200ri + float(id200(i))
endif
c
if (iohc_kaplan(i) .lt. 9999) then
ic8 = ic8 + 1
rhcri = rhcri + float(iohc_kaplan(i))
endif
enddo
c
irimin=1
if (ic2 .ge. irimin) then
shrdri = shrdri/float(ic2)
else
shrdri = 999.
endif
c
if (ic7 .ge. irimin) then
shdcri = shdcri/float(ic7)
else
shdcri = 999.
endif
c
if (ic5 .ge. irimin) then
rhlori= rhlori/float(ic5)
else
rhlori = 999.
endif
c
if (ic6 .ge. irimin) then
d200ri=d200ri/float(ic6)
else
d200ri = 999.
endif
c
sstri = sstri/float(ic3)
potri = potri/float(ic4)
rhcri = rhcri/float(ic8)
c
c ++ End RI variable input
c
c **** Storm classification routine
if (isgnlat .lt. 0.0) then
do k=0,mft
iv500(k) = -iv500(k)
iv300(k) = -iv300(k)
iv20c(k) = -iv20c(k)
enddo
endif
c
c Fix etcl sst to be an integer for etclass routine
isst_etcl = nint(etcl_sst_set%sst_final*10.0)
do k=0,mft
if (isst_etcl(k) .gt. 9000) then
isst_etcl(k)=9999
endif
enddo
c
c call etclass(mft,etcl_sst_set%sst_final,it150,iepos,itgrd,iv500,
call etclass(mft,isst_etcl,it150,iepos,itgrd,iv500,
+ iv300,ishdc,
+ iv20c,cyt,cmagt,igoes,imiss,rmiss,stype)
c
c **** Preliminary call to global RII in case RI probability is used
c as predictor in SHIPS or LGEM.
c
c Set unit number to negative to skip writing in ric2o routine
lutemp = -99
c
c interpolate missing values in vsstl
rmiss1 = rmiss/10.0
call patch0(ric2o_sst_set%vsst_final,rmiss1,mft)
call patch0(sst_lgem,rmiss1,mft)
call patch0(vsst_lgem,rmiss1,mft)
call ric2o(mft,ibasin,perri,vmx,ric2o_sst_set%vsst_final,
+ ishdc,id200,igoes,
+ iohc_ric2,irhlo,itwac,imiss,rmiss,
+ sname,natcf8,ismon,isday,isyr4,istime,lutemp,
+ pprobric2,potmin,ierr)
c
c Check for missing ri probability
if (pprobric2 .gt. 99.0 .or. pprobric2 .lt. 0.0) then
pprobric2 = 5.0
endif
c
write(lulg,310) pprobric2,ibasin
write( 6,310) pprobric2,ibasin
310 format(/,'Preliminary RII probability = ',f5.1,' ibasin=',i2,/)
c
c **** End of Preliminary RII call ****
c
c **** Start of SHIPS Model ****
c
c Specify the SHIPS independent variable names
varlab(1) = 'VMAX '
varlab(2) = 'PER '
varlab(3) = 'ADAY '
varlab(4) = 'SPDX '
varlab(5) = 'PSLV '
varlab(6) = 'VPER '
varlab(7) = 'PC20 '
varlab(8) = 'GSTD '
varlab(9) = 'RIIP '
varlab(10) = 'POT '
varlab(11) = 'SHDC '
varlab(12) = 'T200 '
varlab(13) = 'T250P'
varlab(14) = 'EPOS '
varlab(15) = 'RHMD '
varlab(16) = 'TWAT '
varlab(17) = 'Z850 '
varlab(18) = 'D200 '
varlab(19) = 'LSHDC'
varlab(20) = 'VSHDC'
varlab(21) = 'POT2 '
varlab(22) = 'NOHC '
varlab(23) = 'SDIR '
varlab(24) = 'SHGC '
varlab(25) = 'TADV '
varlab(26) = 'G200 '
varlab(27) = 'RLAT '
varlab(28) = 'SST '
c varlab(27) = 'PW06 '
c
do 99 kk=1,mft
c
c Calculate the average mpi
vssta = 0.0
ccount = 0.0
do 20 k=0,kk
c this is for SHIPS
vssttm = ships_sst_set(k)%vsst_final
c
if (vssttm .lt. 200.0) then
vssta = vssta + vssttm
ccount = ccount + 1.0
endif
20 continue
c
if (ccount .gt. 0.9) then
vssta = vssta/ccount
else
vssta = 100.0
endif
tapot = vssta-vmx
c
c Calculate time averaged SST
tasst = 0.0
ccount = 0.0
do k=0,kk
if (sst_ships(k) .lt. 50.0) then
tasst = tasst + sst_ships(k)
ccount = ccount + 1.0
endif
enddo
if (ccount .gt. 0.9) then
tasst = tasst/ccount
else
tasst = 27.5
endif
c
c Calculate time averaged OHC
taohc = 0.0
ccount= 0.0
do k=0,kk
c print*, "DEBUG iohc_ships ",iohc_ships(k),ships_sst_set(k)%vsst_final
if (iohc_ships(k) .lt. imiss) then
taohc = taohc + float(iohc_ships(k))
ccount= ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
taohc = taohc/ccount
else
c taohc = float(iohcsm)
taohc = float(iohc_sm_ships)
endif
c
c taohc = taohc - rthresh
taohc = taohc - ohc_thresh_ships
if (taohc .lt. 0.0) taohc = 0.0
c
c print*,"DEBUG TEST OHC ",taohc,ohc_thresh_ships,iships_ohc_sm
c print*, "DEBUG TEST IOHC ", iohc_ships
c Calculate time-averaged rlat
tarlat = 0.0
ccount = 0.0
do k=0,kk
if (abs(rlat(k)) .lt. 900.0) then
tarlat = tarlat + abs(rlat(k))
ccount = ccount + 1.0
endif
enddo
if (ccount .gt. 0.9) then
tarlat = tarlat/ccount
endif
c Calculate the time-averaged shear and talshr
tashr = 0.0
talshr = 0.0
ccount = 0.0
c
do 25 k=0,kk
if (ishr(k) .lt. 9000 .and. rlat(k) .lt. 900.0) then
shrtem = 0.1*float(ishr(k))
c
tashr = tashr + shrtem
talshr = talshr + shrtem*sin(dtr*rlat(k))
ccount = ccount + 1.0
endif
25 continue
c
if (ccount .gt. 0.9) then
tashr = tashr/ccount
talshr = talshr/ccount
else
tashr = 18.0
talshr = 7.5
endif
c
c Calculate the time-averaged modified shear and lshr
tashdc = 0.0
talshdc = 0.0
ccount = 0.0
c
do k=0,kk
if (ishdc(k) .lt. 9000 .and. rlat(k) .lt. 900.0) then
shrtem = 0.1*float(ishdc(k))
c
tashdc = tashdc + shrtem
talshdc = talshdc + shrtem*sin(dtr*rlat(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tashdc = tashdc/ccount
talshdc = talshdc/ccount
else
tashdc = 18.0
talshdc = 7.5
endif
c
c Calculate the time-averaged log-transformed shear and lshr
tashdct = 0.0
talshdct = 0.0
ccount = 0.0
c
do k=0,kk
if (shdct(k) .lt. rmiss .and. rlat(k) .lt. 900.0) then
c
tashdct = tashdct + shdct(k)
talshdct = talshdct + shdct(k)*sin(dtr*rlat(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tashdct = tashdct/ccount
talshdct = talshdct/ccount
else
tashdct = 18.0
talshdct = 7.5
endif
c
c Calculate time-averaged shear direction variable
tasdir = 0.0
ccount = 0.0
c
do k=0,kk
if (sdir(k) .lt. rmiss) then
tasdir = tasdir + sdir(k)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tasdir = tasdir/ccount
else
tasdir = 0.0
endif
c
c Calculate time-averaged generalized shear
tashgc = 0.0
ccount = 0.0
c
do k=0,kk
if (shgc(k) .lt. rmiss) then
tashgc = tashgc + shgc(k)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tashgc = tashgc/ccount
else
tashgc = 0.0
endif
c
c Calculate the time-averaged eddy fluxes
tarefc = 0.0
ccount = 0.0
c
do 26 k=0,kk
if (irefc(k) .lt. 9000) then
tarefc = tarefc + float(irefc(k))
ccount = ccount + 1.0
endif
26 continue
c
if (ccount .gt. 0.9) then
tarefc = tarefc/ccount
else
tarefc = 0.0
endif
c
c Calculate the initial value of eddy fluxes
if (irefc(0) .lt. 9000) then
refc0 = float(irefc(0))
else
refc0 = 0.0
endif
c
c Calculate the time-averaged t200
tat200 = 0.0
ccount = 0.0
c
do k=0,kk
if (it200(k) .lt. 9000) then
tat200 = tat200 + 0.1*float(it200(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tat200 = tat200/ccount
else
tat200 = -53.2
endif
c
c Calculate the time-averaged g200
tag200 = 0.0
ccount = 0.0
c
do k=0,kk
if (ig200(k) .lt. 9000) then
tag200 = tag200 + 0.1*float(ig200(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tag200 = tag200/ccount
else
tag200 = 0.0
endif
c
c Calculate the time-averaged t250 threshold variable
tat250 = 0.0
ccount = 0.0
c
do 27 k=0,kk
if (it250(k) .lt. 9000) then
t250t = 0.1*float(it250(k))
t250t = t250t - t250th
if (t250t .gt. 0.0) t250t = 0.0
tat250 = tat250 + t250t
ccount = ccount + 1.0
endif
27 continue
c
if (ccount .gt. 0.9) then
tat250 = tat250/ccount
else
tat250 = 0.0
endif
c
c Calculate the time-averaged u200
tau200 = 0.0
ccount = 0.0
c
do 28 k=0,kk
if (iu200(k) .lt. 9000) then
tau200 = tau200 + 0.1*float(iu200(k))
ccount = ccount + 1.0
endif
28 continue
c
if (ccount .gt. 0.9) then
tau200 = tau200/ccount
else
tau200 = 8.0
endif
c
c Calculate the time-averaged z850
taz850 = 0.0
ccount = 0.0
c
do k=0,kk
if (iz850(k) .lt. 9000) then
taz850 = taz850 + float(iz850(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
taz850 = taz850/ccount
else
taz850 = 25.0
endif
c
c Calculate the time-averaged d200
tad200 = 0.0
ccount = 0.0
c
do k=0,kk
if (id200(k) .lt. 9000) then
tad200 = tad200 + float(id200(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tad200 = tad200/ccount
else
tad200 = 0.0
endif
c
c Calculate the time-averaged tadv
tatadv = 0.0
ccount = 0.0
c
do k=0,kk
if (id200(k) .lt. 9000) then
tatadv = tatadv + float(itadv(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tatadv = tatadv/ccount
else
tatadv = 1.5
endif
c
c Calculate the time-averaged pw06
tapw06 = 0.0
ccount = 0.0
c
do k=0,kk
if (ipw06(k) .lt. 9000) then
tapw06 = tapw06 + float(ipw06(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tapw06 = tapw06/ccount
else
tapw06 = 50.0
endif
c
c Calculate the time-averaged rhlo,rhmd and rhhi
tarhlo = 0.0
tarhmd = 0.0
tarhhi = 0.0
ccount = 0.0
c
do k=0,kk
if (irhlo(k) .lt. 9000 .and. irhhi(k) .lt. 9000 .and.
+ irhmd(k) .lt. 9000) then
tarhlo = tarhlo + float(irhlo(k))
tarhmd = tarhmd + float(irhmd(k))
tarhhi = tarhhi + float(irhhi(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
tarhlo = tarhlo/ccount
tarhmd = tarhmd/ccount
tarhhi = tarhhi/ccount
else
tarhlo = 65.0
tarhmd = 53.0
tarhhi = 45.0
endif
c
c Calculate the time-averaged thetae variable (epos)
taepos = 0.0
ccount = 0.0
c
do k=0,kk
if (iepos(k) .lt. 9000) then
taepos = taepos + 0.1*float(iepos(k))
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.9) then
taepos = taepos/ccount
else
taepos = 12.0
endif
c
c Calculate time tendency of GFS mean tangential wind
twat = twac(kk)-twac(0)
c
c Specify the values of the independent variables
var(1) = vmx
var(2) = per
var(3) = aday
var(4) = spdx
var(5) = pslv
var(6) = vper
var(7) = pc20
var(8) = gstd
var(9) = pprobric2
var(10) = tapot
var(11) = tashdct
var(12) = tat200
var(13) = tat250
var(14) = taepos
var(15) = tarhmd
var(16) = twat
var(17) = taz850
var(18) = tad200
var(19) = talshdct
var(20) = vmx*tashdct
var(21) = tapot*tapot
var(22) = taohc
var(23) = tasdir
var(24) = tashgc
var(25) = tatadv
var(26) = tag200
var(27) = tarlat
var(28) = tasst
c var(27) = tapw06
c
c Normalize the independent variables
do 29 j=1,nvar
var(j) = (var(j) - xbar(j,kk))/xsig(j,kk)
29 continue
c
c Save variable values for later output
do 31 j=1,nvar
vart(j,kk) = var(j)
31 continue
c
c Calculate the intensity change
delv(kk) = 0.0
do 30 n=1,nvar
delv(kk) = delv(kk) + var(n)*coef(n,kk)
30 continue
delv(kk) = ybar(kk) + delv(kk)*ysig(kk)
v(kk) = vmx + delv(kk)
if (rlat(kk) .ge. rmiss) then
v(kk) = 0.0
delv(kk) = 0.0
endif
c
if (rlat(kk) .lt. rmiss) then
write(lulg,300) kk*idelt,vmx,v(kk),delv(kk)
300 format(//,1x,i3,' hr forecast, v(0)=',f5.1,' v=',f6.1,
+ ' dv=',f6.1,/)
do 35 n=1,nvar
write(lulg,305) n,varlab(n),var(n),coef(n,kk),
+ var(n)*coef(n,kk),
+ var(n)*coef(n,kk)*ysig(kk)
305 format(1x,i2,1x,a5,' var=',f8.2,' coef=',f9.4,
+ ' dvn=',f6.1,' dv=',f6.1)
35 continue
endif
99 continue
c
c Check the forecast for dissipation
dissloop: do k=1,mft
if (v(k) .lt. vdiss) then
do kk=k,mft
v(kk) = 0.0
enddo
exit dissloop
endif
enddo dissloop
c
c Write the forecast to integer array for output
iv(0) = ifix(v(0))
do 45 k=1,mft
iv(k) = ifix(v(k) + 0.49)
45 continue
c
c Write required predictors
c to character array for later output
do 50 k=0,mft
if (ishr(k) .lt. 0 .or. ishr(k) .gt. 9000
+ .or. k .gt. ishrx) then
cshr(k) = ' N/A'
else
shrt = 0.1*float(ishr(k))
ishrt = ifix(shrt +0.49)
write(cshr(k), 200) ishrt
200 format(i4)
endif
c
if (ishdc(k) .lt. 0 .or. ishdc(k) .gt. 9000
+ .or. k .gt. ishrx) then
cshdc(k) = ' N/A'
else
shrt = 0.1*float(ishdc(k))
ishrt = ifix(shrt +0.49)
write(cshdc(k), 200) ishrt
endif
c
if (ishgc(k) .gt. 9000
+ .or. k .gt. ishrx) then
cshgc(k) = ' N/A'
else
shrt = shgc(k)
ishrt = ifix(shrt +0.49)
write(cshgc(k), 200) ishrt
endif
c
if (ishtd(k) .gt. 9000) then
cshtd(k) = ' N/A'
else
write(cshtd(k), 200) ishtd(k)
endif
c
if (irefc(k) .gt. 9000 .or. k .gt. irefcx) then
crefc(k)= ' N/A'
else
refct = float(irefc(k))
irefct = ifix(refct+0.49)
write(crefc(k), 200) irefct
endif
c
if (it200(k) .gt. 9000 .or. k .gt. it200x) then
ct200(k)= ' N/A'
else
t200t = 0.1*float(it200(k))
write(ct200(k),210) t200t
endif
c
if (ig200(k) .gt. 9000) then
cg200(k) = ' N/A'
else
write(cg200(k),210) g200(k)
endif
c
if (iu200(k) .gt. 9000) then
cu200(k) = ' N/A'
else
iu200t = ifix(0.1*iu200(k))
write(cu200(k), 200) iu200t
endif
c
if (irhlo(k) .gt. 9000) then
crhlo(k) = ' N/A'
else
write(crhlo(k), 200) irhlo(k)
endif
c
if (irhhi(k) .gt. 9000) then
crhhi(k) = ' N/A'
else
write(crhhi(k), 200) irhhi(k)
endif
c
if (irhmd(k) .gt. 9000) then
crhmd(k) = ' N/A'
else
write(crhmd(k), 200) irhmd(k)
endif
c
if (iepos(k) .gt. 9000) then
cepos(k) = ' N/A'
else
iepost = ifix(0.1*iepos(k))
write(cepos(k), 200) iepost
endif
c
if (iz850(k) .gt. 9000) then
cz850(k) = ' N/A'
else
write(cz850(k), 200) iz850(k)
endif
c
if (id200(k) .gt. 9000) then
cd200(k) = ' N/A'
else
write(cd200(k), 200) id200(k)
endif
c
if (itadv(k) .gt. 9000) then
ctadv(k) = ' N/A'
else
write(ctadv(k), 200) itadv(k)
endif
c
if (iohc_ships(k) .gt. 9000) then
cnohc(k) = ' N/A'
else
write(cnohc(k), 200) iohc_ships(k)
endif
c
50 continue
c
c Write the sst and mpi to character arrays for output.
c Similarly for dist,lat,lon
do 55 k=0,mft
c---------------SST and VSST char variables----BEGIN-------
c WSST
call get_sst_print_variables(k,mft,wsst_set)
c SST DSST
call get_sst_print_variables(k,mft,dsst_set)
c SST DSTA
call get_sst_print_variables(k,mft,dsta_set)
c SST NSST
call get_sst_print_variables(k,mft,nsst_set)
c SST NSTA
call get_sst_print_variables(k,mft,nsta_set)
c---------------SST and VSST char variables---END-------
c
if (idist(k) .gt. 9000) then
cdist(k) = ' N/A'
else
write(cdist(k),200) idist(k)
endif
c
if (rlat(k) .gt. 900.0 .or. rlon(k) .gt. 900.0) then
clat(k) = ' N/A'
clon(k) = ' N/A'
c remove track points in text file for 126-168 h
else if (k .ge. 21 .and. k .le. 28) then
clat(k) = 'xx.x'
clon(k) = 'xxx.x'
else
write(clat(k),210) rlat(k)
write(clon(k),210) sgnlon*rlon(k)
210 format(f5.1)
endif
c
if (cmagt(k) .ge. rmiss) then
cmagc(k) = ' N/A'
else
icmag = ifix(cmagt(k)+0.49)
write(cmagc(k),200) icmag
endif
55 continue
c
c Check for special track model name and exclude printing of lat/lons
if (tmname .eq. 'OFCP' .or. tmname .eq. 'OFPI') then
do k=1,mft
clat(k) = 'xx.x'
clon(k) = 'xxx.x'
enddo
endif
c
c Put track in an array for inclusion in atcf file
do 56 i=0,mft
ilatt(i) = ifix(10.0*rlat(i))
ilont(i) = ifix(10.0*rlon(i))
c
if (ilat(i) .gt. 9000 .or. ilon(i) .gt. 9000) then
ilatt(i) = 0
ilont(i) = 0
endif
56 continue
c
c Run the decay model to include land effects
do i=1,mft+1
ftime(i) = 0.0
rlatd(i) = 0.0
rlond(i) = 0.0
vmaxo(i) = 0.0
vmaxd(i) = 0.0
ddland(i) = 0.0
enddo
c
do i=0,mft
ftime(i+1) = delt*float(i)
c
if (v(i) .ge. vdiss) then
vmaxo(i+1) = v(i)
else
vmaxo(i+1) = 0.0
endif
c
if (rlat(i) .lt. 999. .and. rlon(i) .lt. 999.) then
rlatd(i+1) = rlat(i)
rlond(i+1) = rlon(i)
else
rlatd(i+1) = 0.0
rlond(i+1) = 0.0
endif
enddo
c
c Use persistence for cases with v, but with missing lat/lon
perloop: do k=2,mft+1
if ( vmaxo(k) .gt. 0.0 .and.
+ (rlat(k-1) .ge. 999.0 .or. rlon(k-1) .ge. 999.0) ) then
do kk=k,mft+1
rlatd(kk) = rlatd(k-1)
rlond(kk) = rlond(k-1)
enddo
c
exit perloop
endif
enddo perloop
c
call decay(ftime,rlatd,rlond,vmaxo,vmaxd,ddland,lulg)
c
do i=1,mft+1
ivd(i-1) = ifix(vmaxd(i) + 0.49)
if (i .ne. 1 .and. ivd(i-1) .lt. ivdiss) ivd(i-1)=0
enddo
c
c Write SHIPS forecasts in a modified form of the old ATCF format
c to ships.tst file
c if (ioper .eq. 1) then
if (ioper .eq. 2) then
slab = '98SHIP'
elseif (ioper .eq. -1) then
slab = '92SHPR'
else
slab = '92SHxx'
slab(5:6) = runid
endif
c
do k=1,mft7
ilat7(k) = 0
ilon7(k) = 0
iv7(k) = 0
enddo
c
do k=1,mft
ilat7(k) = ilatt(k)
ilon7(k) = isgnlon*ilont(k)
iv7(k) = iv(k)
enddo
c
if (ipf7 .eq. 1) then
write(luts,401) slab,isyr,ismon,isday,istime,
+ (ilat7(i),ilon7(i),i=2,mft7,2),
+ (iv7(k),k=2,mft7,2),natcf
401 format(a6,4(i2.2),28(i4),14(i3),1x,a6)
else
write(luts,400) slab,isyr,ismon,isday,istime,
+ (ilat7(i),ilon7(i),i=2,mft5,2),
+ (iv7(k),k=2,mft5,2),natcf
400 format(a6,4(i2.2),20(i4),10(i3),1x,a6)
endif
c
c if (ioper .eq. 1) then
if (ioper .eq. 2) then
slab = '98DSHP'
elseif (ioper .eq. -1) then
slab = '92DSHR'
else
slab = '92DSxx'
slab(5:6) = runid
endif
c
do k=1,mft7
iv7(k) = 0
enddo
c
do k=1,mft
iv7(k) = ivd(k)
enddo
c
if (ipf7 .eq. 1) then
write(luts,401) slab,isyr,ismon,isday,istime,
+ (ilat7(i),ilon7(i),i=2,mft7,2),
+ (iv7(k),k=2,mft7,2),natcf
else
write(luts,400) slab,isyr,ismon,isday,istime,
+ (ilat7(i),ilon7(i),i=2,mft5,2),
+ (iv7(k),k=2,mft5,2),natcf
endif
c
c
c close(luts)
c
c Write SHIPS forecasts in the new ATCF format to ships.dat file
c SHIPS
kkount=0
do k=1,mft+1,2
kkount=kkount+1
ifship(kkount,1) = ilatt(k-1)
ifship(kkount,2) = ilont(k-1)*ifix(sgnlon)
ifship(kkount,3) = iv(k-1)
enddo
c call newWriteAidRcd(luat,strmid,aymdh,"SHIP",ifship)
call writeaidlocal2(luat,strmid,aymdh,'SHIP',
+ ilatt,ilont,iv,stype,itimet,mft,minc)
c
c Decay SHIPS
kkount=0
do k=1,mft+1,2
kkount=kkount+1
ifship(kkount,3) = ivd(k-1)
enddo
c call newWriteAidRcd(luat,strmid,aymdh,"DSHP",ifship)
call writeaidlocal2(luat,strmid,aymdh,'DSHP',
+ ilatt,ilont,ivd,stype,itimet,mft,minc)
c
c Write intensities to character array for later output
write( cv(0),200) iv(0)
write(cvd(0),200) ivd(0)
c
do k=1,mft
if (iv(k) .ge. ivdiss) then
write(cv(k),200) iv(k)
else
cv(k) = ' N/A'
endif
c
if (ivd(k) .ge. ivdiss) then
write(cvd(k),200) ivd(k)
else
cvd(k) = ' N/A'
endif
enddo
c
c **** End of SHIPS Model ****
c
c **** Start LGE Model ****
c
c Initialize LGE forecast to zero
nftl = mft+1
do k=1,nftl
vmaxl(k) = 0.0
enddo
c
call lgecal(ibasin,nftl,ftime,vdiss,vmaxl,ierr,lulg)
c
c ++ LGE model output
do i=1,mft+1
ivl(i-1) = ifix(vmaxl(i) + 0.49)
enddo
c
c Write LGE forecast to the ATCF
kkount=0
do k=1,mft+1,2
kkount=kkount+1
ifship(kkount,3) = ivl(k-1)
enddo
c call newWriteAidRcd(luat,strmid,aymdh,"LGEM",ifship)
call writeaidlocal2(luat,strmid,aymdh,'LGEM',
+ ilatt,ilont,ivl,stype,itimet,mft,minc)
c
c Write LGE forecast in a modified form
c of the old ATCF format to the ships.tst file
c if (ioper .eq. 1) then
if (ioper .eq. 2) then
slab = '98LGEM'
elseif (ioper .eq. -1) then
slab = '92LGER'
else
slab = '92LGxx'
slab(5:6) = runid
endif
c
do k=1,mft7
iv7(k) = 0
enddo
c
do k=1,mft
iv7(k) = ivl(k)
enddo
c
if (ipf7 .eq. 1) then
write(luts,401) slab,isyr,ismon,isday,istime,
+ (ilat7(i),ilon7(i),i=2,mft7,2),
+ (iv7(k),k=2,mft7,2),natcf
else
write(luts,400) slab,isyr,ismon,isday,istime,
+ (ilat7(i),ilon7(i),i=2,mft5,2),
+ (iv7(k),k=2,mft5,2),natcf
endif
c
close(luts)
c
c Write intensities to character array for later output
write(cvl(0),200) ivl(0)
c
do k=1,mft
c
if (ivl(k) .gt. 0) then
write(cvl(k),200) ivl(k)
else
cvl(k) = ' N/A'
endif
enddo
c
c **** End LGE Model ****
c
c ++ Write storm classification to separate output file if needed
if (ioper .ne. 1) then
ietw=1
else
ietw=0
endif
c
if (ietw .eq. 1) then
c Write storm classification info to an output file
open(file='stclass.dat',unit=40,form='formatted',
+ status='replace')
write(40,640) strmid,aymdh,(stype(k),k=0,mft)
640 format(a8,1x,a10,1x,21(a4,1x))
close(40)
endif
c
c **** Begin SHIPS text file section ****
c
c Write GFS Vortex mean tangential wind to an array for printing
do k=0,mft
if (itwac(k) .eq. imiss) then
ctwac(k) = 'LOST'
else
ttem = 1.944*0.1*float(itwac(k))
ittem= ifix(ttem)
write(ctwac(k),701) ittem
701 format(i4)
endif
enddo
c
if (ibasin .eq. 1) then
c temchar= '* ATLANTIC SHIPS INTENSITY FORECAST *'
temchar= '* ATLANTIC 2024 SHIPS INTENSITY FORECAST *'
elseif (ibasin .eq. 2) then
c temchar= '* EAST PACIFIC SHIPS INTENSITY FORECAST *'
temchar= '* EAST PACIFIC 2024 SHIPS INTENSITY FORECAST *'
elseif (ibasin .eq. 3) then
temchar= '* WEST PACIFIC SHIPS INTENSITY FORECAST *'
elseif (ibasin .eq. 4) then
temchar= '* N. INDIAN OCEAN SHIPS INTENSITY FORECAST *'
elseif (ibasin .eq. 5) then
temchar= '* S. HEMISPHERE SHIPS INTENSITY FORECAST *'
else
temchar = '* UNRECOGNIZED BASIN IN SHIPS FORECAST *'
endif
write(lush,500) temchar
500 format(33x,a50)
c
if (ipgoes .eq. 0 .and. ipohc .le. 1) then
temchar = '* IR SAT DATA AVAILABLE, OHC AVAILABLE *'
elseif (ipgoes .ne. 0 .and. ipohc .le. 1) then
temchar = '* IR SAT DATA PROXY USED, OHC AVAILABLE *'
elseif (ipgoes .eq. 0 .and. ipohc .gt. 1) then
temchar = '* IR SAT DATA AVAILABLE, OHC PROXY USED *'
else
temchar = '* IR SAT DATA PROXY USED, OHC PROXY USED *'
endif
write(lush,500) temchar
c
write(lush,508) sname,natcf8,ismon,isday,isyr,istime
508 format(33x,'*',2x,a10,2x,a8,2x,i2.2,'/',i2.2,'/',i2.2,2x,
+ i2.2,' UTC *')
c
write(lush,510) (i*6,i=0,4),(i*6,i=6,mft,2)
510 format(/,
+ 'TIME (HR) ',21(3x,i3))
write(lush,512) (cv(i),i=0,4),(cv(i),i=6,mft,2)
512 format('V (KT) NO LAND',21(2x,a4))
write(lush,514) (cvd(i),i=0,4),(cvd(i),i=6,mft,2)
514 format('V (KT) LAND ',21(2x,a4))
write(lush,515) (cvl(i),i=0,4),(cvl(i),i=6,mft,2)
515 format('V (KT) LGEM ',21(2x,a4))
write(lush,615) (stype(i),i=0,4),(stype(i),i=6,mft,2)
615 format('Storm Type ',21(2x,a4))
c
write(lush,516) (cshdc(i),i=0,4),(cshdc(i),i=6,mft,2)
516 format(/,
+ 'SHEAR (KT) ',21(2x,a4))
write(lush,519) (cshgc(i),i=0,4),(cshgc(i),i=6,mft,2)
519 format('SHEAR ADJ (KT)',21(2x,a4))
write(lush,517) (cshtd(i),i=0,4),(cshtd(i),i=6,mft,2)
517 format('SHEAR DIR ',21(2x,a4))
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c START writing SST to SHIPS output file
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c ------------------- LEGACY LINES---------START-----
c write(lush,6201) (csst_dsst(i),i=0,4),(csst_dsst(i),i=6,mft,2)
c 6201 format('SST (C) ',21(2x,a4))
c vsst dsst
c write(lush,6221) (cvsst_dsst(i),i=0,4),
c + (cvsst_dsst(i),i=6,mft,2)
c 6221 format('POT. INT. (KT)',21(2x,a4))
c if (icione .ne. 0) then
c write(lush,6231) (cvsstc_dsst(i),i=0,4),
c + (cvsstc_dsst(i),i=6,mft,2)
c 6231 format('ADJ. POT. INT.',21(2x,a4))
c endif
c ------------------- LEGACY LINES---------END-----
c ------------------- UPDATED LINES---------START-----
write(lush,6201) (ships_sst_set(i)%csstn,i=0,4),
+ (ships_sst_set(i)%csstn,i=6,mft,2)
6201 format('SST (C) ',29(2x,a4))
c vsst dsst
write(lush,6221) (ships_sst_set(i)%cvsstn,i=0,4),
+ (ships_sst_set(i)%cvsstn,i=6,mft,2)
6221 format('POT. INT. (KT)',29(2x,a4))
if (ships_sst_params%icione_vsst .ne. 0) then
write(lush,6231) (ships_sst_set(i)%cvsstc,i=0,4),
+ (ships_sst_set(i)%cvsstc,i=6,mft,2)
6231 format('ADJ. POT. INT.',29(2x,a4))
endif
c ------------------- UPDATED LINES---------END-----
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c END writing SST to SHIPS output file
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c endif
c ===============================================
write(lush,524) (ct200(i),i=0,4),(ct200(i),i=6,mft,2)
524 format('200 MB T (C) ',29(1x,a5))
write(lush,525) (cg200(i),i=0,4),(cg200(i),i=6,mft,2)
525 format('200 MB VXT (C)',29(1x,a5))
write(lush,526) (cepos(i),i=0,4),(cepos(i),i=6,mft,2)
526 format('TH_E DEV (C) ',29(2x,a4))
write(lush,528) (crhmd(i),i=0,4),(crhmd(i),i=6,mft,2)
528 format('700-500 MB RH ',29(2x,a4))
c write(lush,518) (crefc(i),i=0,4),(crefc(i),i=6,mft,2)
c 518 format('MO FLX (M/S/D)',21(2x,a4))
write(lush,529) (ctwac(i),i=0,4),(ctwac(i),i=6,mft,2)
529 format('MODEL VTX (KT)',29(2x,a4))
write(lush,533) (cz850(i),i=0,4),(cz850(i),i=6,mft,2)
533 format('850 MB ENV VOR',29(2x,a4))
write(lush,531) (cd200(i),i=0,4),(cd200(i),i=6,mft,2)
531 format('200 MB DIV ',29(2x,a4))
write(lush,537) (ctadv(i),i=0,4),(ctadv(i),i=6,mft,2)
537 format('700-850 TADV ',29(2x,a4))
write(lush,530) (cdist(i),i=0,4),(cdist(i),i=6,mft,2)
530 format('LAND (KM) ',29(2x,a4))
write(lush,532) (clat(i),i=0,4),(clat(i),i=6,mft,2)
532 format('LAT (DEG N) ',29(1x,a5))
c
if (sgnlon .lt. 0.0) then
write(lush,534) (clon(i),i=0,4),(clon(i),i=6,mft,2)
534 format('LONG(DEG W) ',29(1x,a5))
else
write(lush,538) (clon(i),i=0,4),(clon(i),i=6,mft,2)
538 format('LONG(DEG E) ',29(1x,a5))
endif
c
write(lush,535) (cmagc(i),i=0,4),(cmagc(i),i=6,mft,2)
535 format('STM SPEED (KT)',29(1x,a5))
write(lush,716) (cnohc(i),i=0,4),
+ (cnohc(i),i=6,mft,2)
716 format('HEAT CONTENT ',29(1x,a5))
c
ispdx = ifix(spdx+0.5)
ispdy = ifix(spdy+0.5)
ipslvm = ifix(xbar(5,1))
c
write(lush,536) tmname,ihead,ispeed,ispdx,ispdy
536 format(/,' FORECAST TRACK FROM ',a4,
+ ' INITIAL HEADING/SPEED (DEG/KT):',I3,'/',I3,
+ ' CX,CY: ',i3,'/',i3)
c
write(lush,542) ivmxm12,ipslv(0),ipslvm
542 format(' T-12 MAX WIND: ',i3,
+ ' PRESSURE OF STEERING LEVEL (MB): ',i4,
+ ' (MEAN=',i3,')')
c
write(lush,718) btstd,btstdm
718 format(' GOES IR BRIGHTNESS TEMP. STD DEV. 50-200 KM RAD: ',
+ f5.1,' (MEAN=',f4.1,')')
write(lush,719) pc20,pc20m
719 format(' % GOES IR PIXELS WITH T < -20 C 50-200 KM RAD: ',
+ f5.1,' (MEAN=',f4.1,')')
write(lush,720) pprobric2
720 format(' PRELIM RI PROB (DV .GE. 35 KT IN 36 HR): ',f6.1)
c
c Set iarr=1 to print array showing contributions to
c intensity change from individual variables,
c else set iarr=0
iarr=1
c
if (iarr .eq. 1) then
c Calculate contributions to intensity change from various variables
c or groups of variables.
c Note: This section of code depends on specific order of
c independent variables
c
c Specify labels for variables
ngrp = 19
cdvlab(1) = 'SAMPLE MEAN CHANGE '
cdvlab(2) = 'SST POTENTIAL '
cdvlab(3) = 'VERTICAL SHEAR MAG '
cdvlab(4) = 'VERTICAL SHEAR ADJ '
cdvlab(5) = 'VERTICAL SHEAR DIR '
cdvlab(6) = 'PERSISTENCE '
cdvlab(7) = '200/250 MB TEMP. '
cdvlab(8) = 'THETA_E EXCESS '
cdvlab(9) = '700-500 MB RH '
cdvlab(10)= 'MODEL VTX TENDENCY '
cdvlab(11)= '850 MB ENV VORTICITY'
cdvlab(12)= '200 MB DIVERGENCE '
cdvlab(13)= '850-700 T ADVEC '
cdvlab(14)= 'ZONAL STORM MOTION '
cdvlab(15)= 'STEERING LEVEL PRES '
cdvlab(16)= 'DAYS FROM CLIM. PEAK'
cdvlab(17)= 'GOES PREDICTORS '
cdvlab(18)= 'OCEAN HEAT CONTENT '
cdvlab(19)= 'RI POTENTIAL '
c
c Calculate intensity change contributions
do 60 kk=1,mft
dvvar( 1,kk) = ybar(kk)
dvvar( 2,kk) = vart( 1,kk)*coef( 1,kk)*ysig(kk) +
+ vart(10,kk)*coef(10,kk)*ysig(kk) +
+ vart(21,kk)*coef(21,kk)*ysig(kk) +
+ vart(27,kk)*coef(27,kk)*ysig(kk) +
+ vart(28,kk)*coef(28,kk)*ysig(kk)
dvvar( 3,kk) = vart(11,kk)*coef(11,kk)*ysig(kk) +
+ vart(19,kk)*coef(19,kk)*ysig(kk) +
+ vart(20,kk)*coef(20,kk)*ysig(kk)
dvvar( 4,kk) = vart(24,kk)*coef(24,kk)*ysig(kk)
dvvar( 5,kk) = vart(23,kk)*coef(23,kk)*ysig(kk)
dvvar( 6,kk) = vart( 2,kk)*coef( 2,kk)*ysig(kk) +
+ vart( 6,kk)*coef( 6,kk)*ysig(kk)
dvvar( 7,kk) = vart(12,kk)*coef(12,kk)*ysig(kk) +
+ vart(13,kk)*coef(13,kk)*ysig(kk) +
+ vart(26,kk)*coef(26,kk)*ysig(kk)
dvvar( 8,kk) = vart(14,kk)*coef(14,kk)*ysig(kk)
dvvar( 9,kk) = vart(15,kk)*coef(15,kk)*ysig(kk)
dvvar(10,kk) = vart(16,kk)*coef(16,kk)*ysig(kk)
dvvar(11,kk) = vart(17,kk)*coef(17,kk)*ysig(kk)
dvvar(12,kk) = vart(18,kk)*coef(18,kk)*ysig(kk)
dvvar(13,kk) = vart(25,kk)*coef(25,kk)*ysig(kk)
dvvar(14,kk) = vart( 4,kk)*coef( 4,kk)*ysig(kk)
dvvar(15,kk) = vart( 5,kk)*coef( 5,kk)*ysig(kk)
dvvar(16,kk) = vart( 3,kk)*coef( 3,kk)*ysig(kk)
dvvar(17,kk) = vart( 7,kk)*coef( 7,kk)*ysig(kk)+
+ vart( 8,kk)*coef( 8,kk)*ysig(kk)
dvvar(18,kk) = vart(22,kk)*coef(22,kk)*ysig(kk)
dvvar(19,kk) = vart( 9,kk)*coef( 9,kk)*ysig(kk)
c
if (rlat(kk) .ge. rmiss) then
do m=1,ngrp
dvvar(m,kk) = 0.0
enddo
endif
c do n=1,nvar
c write(lulg,305) n,varlab(n),var(n),coef(n,kk),
c + var(n)*coef(n,kk),
c + var(n)*coef(n,kk)*ysig(kk)
c enddo
60 continue
c
c Print IC array
write(lush,570)
570 format(/,24x,'INDIVIDUAL CONTRIBUTIONS TO INTENSITY CHANGE')
c
write(lush,572) (kk*6,kk=1,4),(kk*6,kk=6,mft,2)
572 format(22x,20(1x,i3,2x))
write(lush,574)
574 format(24x,
+ '--------------------------------------------------------------
+------------------------------')
do 62 i=1,ngrp
write(lush,576) cdvlab(i),(dvvar(i,kk),kk=1,4),
+ (dvvar(i,kk),kk=6,mft,2)
576 format(2x,a20,20(f5.0,1x))
62 continue
c
write(lush,574)
write(lush,579) (delv(kk),kk=1,4),(delv(kk),kk=6,mft,2)
579 format(2x,'TOTAL CHANGE ',4x,28(f5.0,1x))
endif
c **** End of SHIPS text file section ****
c
c **** Begin rapid intensity indices ****
c
c PARMOD+ 2012 JHT Enhanced RII
rlont = -rlon(0)
rlatt = rlat(0)
write(lush,789) vmx,rlatt,rlont
789 format(/,15x,' CURRENT MAX WIND (KT):',f6.0,
+ ' LAT, LON: ',f6.1,1x,f7.1)
c
c PARMOD+ JHT RII
c ++ Section 3, RII with multiple times and consensue
c
c select which TPW to use
if (itpwtype .eq. 0) then
c observed
do k = 0,mft
iusetpw(k) = istpw(k)
enddo
elseif (itpwtype .eq. 1) then
c model
do k = 0,mft
iusetpw(k) = imtpw(k)
enddo
else
write(lulg,1051)
1051 format(/,'Unknown TPW type; will use observed')
do k = 0,mft
iusetpw(k) = istpw(k)
enddo
endif
c
c Three RI models are used to create a consensus product:
c 1. SHIPS-RII (2012 Operational Disc Analysis version but with TIC)
c 2. Bayesian model (Rozoff and Kossin 2011; WAF)
c 3. logistic regression model (Rozoff and Kossin 2011; WAF)
c
c Extract extra predictors for TIC RII
c
if (iusetpw(18) .ne. imiss) then
tpwri = float(iusetpw(18))/10.0
else
tpwri = rmiss
endif
c
if (iaipc( 1) .ne. imiss .and. iaipc( 1) .gt. -900) then
rirpcri = float(iaipc(1))/100.0
else
rirpcri = rmiss
endif
c
icount = 0
cflxri = 0.0
do k=0,4
if (icflx(k) .ne. imiss) then
cflxri = cflxri + float(icflx(k))
icount = icount + 1
endif
enddo
c
if (icount .ge. 1) then
cflxri = cflxri/float(icount)
else
cflxri = rmiss
endif
c
vmaxri = float(ivmx)
c
c
c PARMOD+ For RII with multipe times and consensus (MTC)
if (iusetpw(0) .ne. imiss) then
tpw1ri = float(iusetpw(0))
else
tpw1ri = rmiss
endif
c
if (iusetpw(1) .ne. imiss) then
tpw2ri = float(iusetpw(1))
else
tpw2ri = rmiss
endif
c
if (iusetpw(3) .ne. imiss) then
tpw4ri = float(iusetpw(3))
else
tpw4ri = rmiss
endif
c
if (iusetpw(5) .ne. imiss) then
tpw6ri = float(iusetpw(5))
else
tpw6ri = rmiss
endif
c
if (iusetpw(6) .ne. imiss) then
tpw7ri = float(iusetpw(6))
else
tpw7ri = rmiss
endif
c
if (iusetpw(11) .ne. imiss) then
tpw12ri = float(iusetpw(11))
else
tpw12ri = rmiss
endif
c
if (iusetpw(12) .ne. imiss) then
tpw13ri = float(iusetpw(12))
else
tpw13ri = rmiss
endif
c
if (iusetpw(17) .ne. imiss) then
tpw18ri = float(iusetpw(17))
else
tpw18ri = rmiss
endif
c
if (iusetpw(18) .ne. imiss) then
tpw19ri = float(iusetpw(18))
else
tpw19ri = rmiss
endif
c
if (iaipc(0) .ne. imiss .and. iaipc(0) .gt. -900) then
aipc1ri = float(iaipc(0))
else
aipc1ri = rmiss
endif
c
if (iaipc(1) .ne. imiss .and. iaipc(1) .gt. -900) then
aipc2ri = float(iaipc(1))
else
aipc2ri = rmiss
endif
c
if (idipc(1) .ne. imiss .and. idipc(1) .gt. -900) then
dipc2ri = float(idipc(1))
else
dipc2ri = rmiss
endif
c
if (ipcnew .eq. 1) then
c Use new PC variables to set input to RII routines
call riipcset(ipc00,ipcm1,ipcm3,mft,imiss,rmiss,
+ rirpcri,aipc1ri,aipc2ri,dipc2ri)
endif
pc1ri = aipc1ri
pc2ri = aipc2ri
c
c PARMOD-
c
c
c write(lush,801)
c 801 format(/,'++++++++ SECTION 3, RII WITH MULTIPLE TIMES ++++++++',
c + /,' AND CONSENSUS FOR JHT ')
c
if (ibasin .eq. 1) then
c
c Atlantic RII Models
c
c Discriminant Analysis Version
c call rapidga_exp_v1_jht_2023(perri,ishdc,id200,
call rapidga(perri,ishgc,id200,
+ kaplan_sst_set%vsst_final,
+ ipw19,sbtri,rirpcri,iohc_kaplan,icflx,vmaxri,sname,
+ natcf8,ismon,isday,isyr,istime,lush,probridisc)
c
c Logistic Regression Version
call rapidloga_2018( perri, rozof_sst_set%vsst_final,
+ rozof_sst_set%sst_final, vmx,
+ iohc_rozof,
+ iu200, ienss, id200, ishgc, ishrg, sir2ri,
+ sbtri, sir5ri, sir10ri, sir12ri, sir13ri,
+ sir14ri, ipw04, ipw19, pc2ri, probrilog)
c
c Bayesian Version
call rapidbaya_2018(perri, ilat, iohc_rozof,
+ iu200, iepss,
+ irhlo, irhmd, irhhi, id200, ishrg,
+ sbtri, sir4ri, sir5ri, sir10ri,
+ sir12ri, sir13ri, sir14ri,
+ rozof_sst_set%vsst_final,
+ rozof_sst_set%sst_final,
+ vmx, ipw03, ipw04, pc2ri, probribay)
c
c Print out a matrix of all RI probabilities & consensus
call rapidprint_2018(probridisc, probrilog, probribay,
+ lush, probricon)
elseif (ibasin .eq. 2) then
c East/Central Pacific RII models
c
c Discriminant Analysis Version
c call rapidge_exp_v1_jht_2023(perri,ishdc,id200,
call rapidge(perri,ishgc,id200,
+ kaplan_sst_set%vsst_final,
+ ipw19,sbtri,rirpcri,iohc_kaplan,icflx,vmaxri,sname,
+ natcf8,ismon,isday,isyr,istime,lush,probridisc)
c
c Logistic Regression Version
call rapidloge_2018( perri, ilat, iepss, ienss, irhlo,
+ irhmd, id200, ishdc, ishgc,
+ rozof_sst_set%vsst_final,
+ vmx, sir5ri, sir6ri, sir10ri, sir12ri,
+ ipw01, ipw04, ipw06, icflx, probrilog)
c
c Bayesian Version
call rapidbaye_2018(perri, ilat, iohc_rozof,
+ iepos, iepss,
+ irhhi, id200, ishdc, ishgc, sbtri, sir5ri,
+ sir7ri, sir9ri, sir10ri, sir11ri, sir12ri,
+ sir14ri,
+ rozof_sst_set%vsst_final,
+ rozof_sst_set%sst_final,
+ vmx, ipw06, probribay)
c
c Print out a matrix of all RI probabilities & consensus
call rapidprint_2018(probridisc, probrilog, probribay,
+ lush, probricon)
c
c PARMOD- JHT RII
elseif (ibasin .ge. 3 .and. ibasin .le. 5) then
call ric2o(mft,ibasin,perri,vmx,ric2o_sst_set%vsst_final,
+ ishdc,id200,igoes,
+ iohc_ric2,irhlo,itwac,imiss,rmiss,
+ sname,natcf8,ismon,isday,isyr4,istime,lush,
+ probric2,potmin,ierr)
endif
c
iprobri=999
if (ibasin .eq. 1 .or. ibasin .eq. 2) then
c
c RI aid
probri=probridisc
do kri=1,mxmul
iprobri(kri)=nint(probri(kri))
enddo
c
else
iprobri(3) = nint(probric2)
endif
c
call rapid_aid(natcf8,ismon,isday,isyr,istime,
+ iprobri(2),iprobri(3),iprobri(4),iprobri(5),
+ luat)
c
if (ibasin .le. 2) then
c Write RI info to log file
write(lulg,860) sname,natcf8,ismon,isday,isyr,istime,vmx,
+ pprobric2,(probri(k),k=1,8)
860 format(/,a10,1x,a8,1x,i2.2,i2.2,i2.2,1x,i2.2,' vmx=',f5.0,
+ ' Pre Prob=',f6.1,
+ ' Oper Probs: ',8(f6.1,1x))
endif
c Write RI to edeck
ibegtau=0
ibri=0
c SHIPS-RII
do kri=1,mxmul
if (iprobri(kri) .lt. 999) then
call writeedecklocal(lued,natcf8,aymdh,"RIOD",
+ iendtau(kri),ilatt(ibri),
+ ilont(ibri),iprobri(kri),
+ ithres(kri),ivmx,
+ ibegtau,iendtau(kri))
endif
enddo
c Logistic Regression Version
do kri=1,mxmul
iprobri(kri) = nint(probrilog(kri))
if (iprobri(kri) .lt. 999) then
call writeedecklocal(lued,natcf8,aymdh,"RIOL",
+ iendtau(kri),ilatt(ibri),
+ ilont(ibri),iprobri(kri),
+ ithres(kri),ivmx,
+ ibegtau,iendtau(kri))
endif
enddo
c Bayesian Version
do kri=1,mxmul
iprobri(kri) = nint(probribay(kri))
if (iprobri(kri) .lt. 999) then
call writeedecklocal(lued,natcf8,aymdh,"RIOB",
+ iendtau(kri),ilatt(ibri),
+ ilont(ibri),iprobri(kri),
+ ithres(kri),ivmx,
+ ibegtau,iendtau(kri))
endif
enddo
c Consensus Version
do kri=1,mxmul
iprobri(kri) = nint(probricon(kri))
if (iprobri(kri) .lt. 999) then
call writeedecklocal(lued,natcf8,aymdh,"RIOC",
+ iendtau(kri),ilatt(ibri),
+ ilont(ibri),iprobri(kri),
+ ithres(kri),ivmx,
+ ibegtau,iendtau(kri))
endif
enddo
close(lued)
c
c **** End rapid intensity index 2017****
c
c **** Begin Annular hurricane index ****
luout = lush
stname = sname
tcfid = natcf8
fn_ships='lsdiag.dat'
fn_ird1='IRRP1.dat'
fn_ird2='IRRP2.dat'
fn_iri1='IRRP1.inf'
fn_iri2='IRRP2.inf'
c
rlon00 = rlon(0)
c write(6,*) 'Call annular'
if (rlon00 .lt. rlonam .and. ibasin .lt. 3) then
call id_annular_op(luout,stname,tcfid,fn_ships,fn_ird1,
+ fn_ird2,fn_iri1,fn_iri2,ianmon,ianday,
+ ianyr,iantime,dfval,ann_prob,
+ ahi_sst_set%sst_final)
endif
c
c **** End Annular hurricane index ****
c
if (ibasin .eq. 1) then
c **** Begin PrSEFoNe (probability of Secondary Eyewall Formation) model ****
c
c At 4 forcast times 0,12,24,36h, 12 features are input into the model:
c 1. Current intensity (kt)
c 2. Latitude of center fix (degrees N) [SHIPS -> LAT / 10]
c 3. Climatological ocean heat content [SHIPS -> PHCN]
c 4. 200mb zonal wind (kt) [SHIPS -> U200 / 10]
c 5. Relative humidity between 500-300mb (%) [SHIPS -> RHHI]
c 6. Symmetric tangential wind (m/s) [SHIPS -> TWAC / 10]
c 7. Surface pressure at outer vortex edge (mb - 1000) [SHIPS -> PENC / 10]
c 8. Vertical shear (kt) [SHIPS -> SHRD / 10]
c 9. Potential intensity (kt) [SHIPS -> VMPI]
c 10. Standard deviation of Tb between 100-300km (C) [SHIPS -> IR00-05 / 10]
c 11. Mean Tb between 20-120km (C) [SHIPS -> IR00-16 / 10]
c 12. IR Tb Principal component #4 (added within subroutine psef_driver)
c
c fill feature matrix for PrSEFoNe model. DTL is included to identify cases
c where storm center is over land, but it's not a feature of the model.
c
c If any one of the standard ships features are missing from lsdiag,
c then flag that fix.
c
imiss_feat=9999
do i=1,4
j=(i-1)*2
iflag_feat(i)=0
if (ivd(j).eq.imiss_feat) iflag_feat(i)=1
if (ilat(j).eq.imiss_feat) iflag_feat(i)=1
if (iohc_prse(j).eq.imiss_feat) iflag_feat(i)=1
if (iu200(j).eq.imiss_feat) iflag_feat(i)=1
if (irhhi(j).eq.imiss_feat) iflag_feat(i)=1
if (itwac(j).eq.imiss_feat) iflag_feat(i)=1
if (ipenc(j).eq.imiss_feat) iflag_feat(i)=1
if (ishr(j).eq.imiss_feat) iflag_feat(i)=1
if (ivmpi2(j).eq.imiss_feat) iflag_feat(i)=1
enddo
c
c flag missing IR00 features
c
iflag_feat_ir=0
if (igoes(4).eq.imiss_feat) iflag_feat_ir=1
if (igoes(15).eq.imiss_feat) iflag_feat_ir=1
do i=1,4
j=(i-1)*2
psef_feat_vec(1,i)=ivd(j) ! vmx at 0,12,24,36h
psef_feat_vec(2,i)=ilat(j)/10. ! lat at 0,12,24,36h
psef_feat_vec(3,i)=iohc_prse(j) ! phcn at 0,12,24,36h
psef_feat_vec(4,i)=iu200(j)/10. ! u200 at 0,12,24,36h
psef_feat_vec(5,i)=irhhi(j) ! rhhi at 0,12,24,36h
psef_feat_vec(6,i)=itwac(j)/10. ! twac at 0,12,24,36h
psef_feat_vec(7,i)=ipenc(j)/10. ! penc at 0,12,24,36h
psef_feat_vec(8,i)=ishr(j)/10. ! shrd at 0,12,24,36h
psef_feat_vec(9,i)=ivmpi2(j) ! vmpi at 0,12,24,36h
psef_feat_vec(10,i)=igoes(4)/10. ! ir00-5 (fixed through all times)
psef_feat_vec(11,i)=igoes(15)/10. ! ir00-16 (fixed through all times)
psef_feat_vec(12,i)=idist(j) ! dtl at 0,12,24,36h
enddo
c call PrSEFoNe driver subroutine:
c luout is for writing to ships output file (ships.txt)
c stname, tcfid, isyr4, ilmon, ilday, itime is for output file header
c psef_feat_vec, iflag_feat, iflag_feat_ir are needed to calculate PrSEFoNe
call psef_driver(luout,stname,tcfid,isyr4,ilmon,ilday,
. iltime,psef_feat_vec,iflag_feat,iflag_feat_ir)
c Print adjusted forecast with climatological mean intensity change
c Adjustments to DSHIPS forecast
call erc_clim_driver(luout,ivd,mft)
c
c **** End PrSEFoNe (probability of Secondary Eyewall Formation) ****
endif
stop '***** SHIPS COMPLETED *****'
c
900 continue
write(lulg,950)
950 format(' Error during file open')
c
stop '***** SHIPS HALTED DUE TO FILE OPEN ERROR *****'
c
c
end
c
c ****** END SHIPS MAIN PROGRAM, BEGIN SUBROUTINES ******
c
subroutine lcheck(ilat0,ilon0,ilatm,ilonm,dt,tland,lmod)
c This routine checks to see if a storm is now over water,
c but was recently over land. If so, lmod is set to 1.
c
c Specify maximum allowable time over land (hr)
tlmax = 2.0
c
lmod = 0
c
c Check and convert lat/lon values
if (ilat0 .eq. 0 .or. ilatm .eq. 0) return
if (ilon0 .eq. 0 .or. ilonm .eq. 0) return
c
rlat0 = 0.1*float(ilat0)
rlatm = 0.1*float(ilatm)
rlon0 = 0.1*float(ilon0)
rlonm = 0.1*float(ilonm)
c
dlat = rlat0-rlatm
dlon = rlon0-rlonm
c
c Check to make sure lat/lon pair are physically close
dmax = dt*50.0/60.0
c
if (abs(dlat) .gt. dmax .or.
+ abs(dlon) .gt. dmax ) return
c
c If storm is currently over land, no adjustment is needed
call gdland_table(rlon0,rlat0,dland0)
if (dland0 .lt. 0.0) return
c
c Divide positions into smaller sub-intervals
dts = 1.0
nts = ifix(0.0001 + dt/dts)
c
dlati = dlat/float(nts)
dloni = dlon/float(nts)
tland = 0.0
do k=0,nts
if (k .eq. 0 .or. k .eq. nts) then
wt = 0.5*dts
else
wt = dts
endif
c
rlatt = rlatm + dlati*float(k)
rlont = rlonm + dloni*float(k)
call gdland_table(rlont,rlatt,dlandt)
c
if (dlandt .lt. 0.0) tland = tland + wt
c write(6,888) rlatt,rlont,dlandt,wt,tland
c 888 format(5(f8.1,1x))
enddo
c
if (tland .ge. tlmax) lmod = 1
c
return
end
c
c ++++BEGIN LGEM MODULE++++
c
subroutine lgecal(ibasin,nft,ftime,vdiss,vmaxl,ierr,lulg)
c This routine makes an intensity forecast using the
c Logistic Growth Equation. This routine constructs the cappa
c array using statistical relationship , and then calls routine
c lgeim for the numerical solution of the LGE equation. Most of the
c required input for this routine is passed through common blocks.
c
c ++ Calling argument variables
dimension ftime(nft),vmaxl(nft)
c
c ++ Common blocks from main iships.f for lgecal routine
common /lgestr/ aday,spdx,pslv,per,vmx,pc20,ohc_thresh_lgem,
+ pprobric2
common /lgetdi/ ishr,iepos,iz850,id200,ishdc,itadv,
+ ipw06,iohc_best_ncoda,iohc_lgem
common /lgetdr/ sst_lgem,vsst_lgem,rlat,rlon,t200,t250,twac,sdir,
+ shdct,shgc,g200
common /lgepar/ rmiss,imiss,ioper
common /lgetst/ coyear
c
parameter (mft=28)
dimension ishr(0:mft),iepos(0:mft),iz850(0:mft),id200(0:mft)
dimension ishdc(0:mft),itadv(0:mft)
dimension ipw06(0:mft)
dimension iohc_best_ncoda(0:mft)
dimension twac(0:mft),sdir(0:mft)
dimension shdct(0:mft),shgc(0:mft)
dimension t200(0:mft),t250(0:mft),g200(0:mft)
dimension sst_lgem(0:mft),vsst_lgem(0:mft)
dimension iohc_lgem(0:mft)
dimension ohc_lgem(0:mft)
dimension rlat(-2:mft),rlon(-2:mft)
character *2 coyear
c
c ++ Local variables for variables derived from common variables
dimension shr(0:mft),epos(0:mft),z850(0:mft),d200(0:mft)
dimension shrl(0:mft),sst_lgem_local(0:mft)
dimension tadv(0:mft)
dimension shdctl(0:mft),pw06(0:mft)
c
c ++ Local variables for coefficients
parameter (nvar=22)
dimension var(nvar),coef(nvar,0:mft)
dimension xbar(nvar,0:mft),xsig(nvar,0:mft)
dimension ybar(0:mft),ysig(0:mft)
character *5 varlab(nvar)
c
c ++ Local variables for lgeim call
parameter (nmx=mft+1)
dimension tlat(nmx),tlon(nmx),vmpi(nmx),cappa(nmx),dland(nmx)
intrinsic trim
c
character *256 fncol
character *256 coef_location
c
c get SHIPS_COEF env variable
call getenv( "SHIPS_COEF", coef_location )
c
dtr = 0.017453
c
c Set the number of variables by basin
nvarb = nvar
if (ibasin .ge. 3) nvarb = nvar -1
c **** Default forecast ****
vmaxl(1) = vmx
do k=2,nft
vmaxl(k) = 0.0
enddo
ierr=1
c
c **** Check of input variables ****
if (nft .lt. 2 ) return
if (vmx .le. 0.0 .or. vmx .ge. rmiss ) return
if (aday .ge. rmiss .or. spdx .ge. rmiss) return
if (pslv .ge. rmiss ) return
if (per .ge. rmiss ) return
if (pc20 .ge. rmiss ) return
c if (ishr(0) .ge. imiss ) return
if (ishdc(0) .ge. imiss ) return
if (iepos(0) .ge. imiss ) return
if (iz850(0) .ge. imiss ) return
if (id200(0) .ge. imiss ) return
if (iohc_lgem(0) .ge. imiss ) return
if (itadv(0) .ge. imiss ) return
if (t200(0) .ge. rmiss ) return
if (t250(0) .ge. rmiss ) return
if (g200(0) .ge. rmiss ) return
if (rlat(0) .ge. rmiss ) return
if (rlon(0) .ge. rmiss ) return
if (vsst_lgem(0) .ge. rmiss ) return
if (sst_lgem(0) .ge. rmiss ) return
if (sdir(0) .ge. rmiss ) return
if (shdct(0) .ge. rmiss ) return
if (shgc(0) .ge. rmiss ) return
if (ipw06(0) .ge. imiss ) return
c
c **** Specify basin-specific variables ****
if (ibasin .eq. 1) then
fncol=trim( coef_location )//'shipl24_coef_atlc.dat'
elseif (ibasin .eq. 2) then
fncol=trim( coef_location )//'shipl24_coef_epac.dat'
elseif (ibasin .eq. 3) then
if (ioper .eq. 1) then
fncol=trim( coef_location )//'shipl13_coef_wpac.dat'
else
fncol='shipl13_coef_wpac.dat'
fncol(6:7) = coyear
endif
elseif (ibasin .eq. 4) then
if (ioper .eq. 1) then
fncol=trim( coef_location )//'shipl13_coef_iocn.dat'
else
fncol='shipl13_coef_iocn.dat'
fncol(6:7) = coyear
endif
elseif (ibasin .eq. 5) then
if (ioper .eq. 1) then
fncol=trim( coef_location )//'shipl13_coef_shem.dat'
else
fncol='shipl13_coef_shem.dat'
fncol(6:7) = coyear
endif
else
write(lulg,*) 'Unexpected basin, ibasin=',ibasin
stop
endif
c
c Open and read coefficient file
lucol=45
open(file=fncol,unit=lucol,form='formatted',status='old',
+ err=900)
print*, "-- opened coeff file fncol ", fncol
c
ierr=2
do 5 k=0,mft
read(lucol,135,end=900,err=900) ybar(k),ysig(k)
135 format(14x,e12.5,1x,e12.5)
c
do 6 n=1,nvarb
read(lucol,136,end=900,err=900)
+ coef(n,k),xbar(n,k),xsig(n,k)
136 format(1x,3(e12.5,1x))
6 continue
5 continue
c
close(lucol)
c **** Variable initialization and specification ****
c
c ++Specify the beta (1/hr) and nn parameters for the MPI term
beta = 1.0/24.0
rnn = 2.5
c
c ++Specify the wind scale and exponent for the alternate LGE model
c where the growth term is k(v/vscale)**rmm. Set vscale=rmm=1.0 ot
c use the standard LGE model equation.
vscale= 1.0
rmm = 1.0
c
c ++Convert integer common block variables to real
do k=0,mft
if (ishdc(k) .lt. imiss) then
shr(k) = 0.1*float(ishdc(k))
if (rlat(k) .lt. rmiss) then
shrl(k) = shr(k)*sin(dtr*rlat(k))
else
shrl(k) = rmiss
endif
else
shr(k) = rmiss
shrl(k) = rmiss
endif
c
if (iepos(k) .lt. imiss) then
epos(k) = 0.1*float(iepos(k))
else
epos(k) = rmiss
endif
c
if (iz850(k) .lt. imiss) then
z850(k) = float(iz850(k))
else
z850(k) = rmiss
endif
c
if (id200(k) .lt. imiss) then
d200(k) = 0.1*float(id200(k))
else
d200(k) = rmiss
endif
c
c if (ibohc(k) .lt. imiss) then
c rhcn(k) = float(ibohc(k))
c else
c rhcn(k) = rmiss
c endif
if (iohc_lgem(k) .lt. imiss) then
ohc_lgem(k) = float(iohc_lgem(k))
else
ohc_lgem(k) = rmiss
endif
c
if (itadv(k) .lt. imiss) then
tadv(k) = float(itadv(k))
else
tadv(k) = rmiss
endif
c
if (ipw06(k) .lt. imiss) then
pw06(k) = float(ipw06(k))
else
pw06(k) = rmiss
endif
c
c sstl(k) = gen_sst(k)
sst_lgem_local(k) = sst_lgem(k)
enddo
c
c Calculate transformed shear times sin(lat)
do k=0,mft
if (shdct(k) .lt. rmiss .and. rlat(k) .lt. rmiss) then
shdctl(k) = shdct(k)*sin(dtr*rlat(k))
else
shdctl(k) = rmiss
endif
enddo
c
c ++ Fill in missing values in time arrays
rmiss1 = rmiss/10.0
call patch0(shr ,rmiss,mft)
call patch0(shrl,rmiss,mft)
call patch0(shdct,rmiss,mft)
call patch0(shdctl,rmiss,mft)
call patch0(shgc,rmiss,mft)
call patch0(sdir,rmiss,mft)
call patch0(epos,rmiss,mft)
call patch0(z850,rmiss,mft)
call patch0(d200,rmiss,mft)
call patch0(t200,rmiss,mft)
call patch0(t250,rmiss,mft)
call patch0(g200,rmiss,mft)
call patch0(sst_lgem_local,rmiss1,mft)
call patch0(vsst_lgem,rmiss1,mft)
call patch0(ohc_lgem,rmiss,mft)
call patch0(tadv,rmiss,mft)
call patch0(pw06,rmiss,mft)
c
c ++ Calculate initial value of cappa from previous 12 hour intensity change
cappa00 = per/(12.0*vmx) + beta*(vmx/vsst_lgem(0))**rnn
c
c **** Calculate cappa ****
c
c ++ Main time loop
do 99 k=0,mft
ktime = k*6
c
c ++ GFS vortex time tendency
klag = 2
klead = 2
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
deltg = float(ke-ks)
if (deltg .le. 0.0) then
twacten = 0.0
else
twacten = 4.0*(twac(ke)-twac(ks))/deltg
endif
c
c ++ time averaged shear and lshear
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tashr = 0.0
tashrl= 0.0
do kk=ks,ke
if (shr(kk) .lt. rmiss .and. shrl(kk) .lt. rmiss) then
tashr = tashr + shr(kk)
tashrl= tashrl+ shrl(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tashr = tashr/ccount
tashrl= tashrl/ccount
else
tashr = rmiss
tashrl= rmiss
endif
c
c ++ time averaged transformed shear, lshear
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tashdct = 0.0
tashdctl= 0.0
do kk=ks,ke
if (shdct(kk) .lt. rmiss .and. shdctl(kk) .lt. rmiss) then
tashdct = tashdct + shdct(kk)
tashdctl= tashdctl+ shdctl(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tashdct = tashdct/ccount
tashdctl= tashdctl/ccount
else
tashdct = rmiss
tashdctl= rmiss
endif
c
c ++ time averaged shear direction parameter
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tasdir = 0.0
do kk=ks,ke
if (sdir(kk) .lt. rmiss) then
tasdir = tasdir + sdir(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tasdir = tasdir/ccount
else
tasdir = rmiss
endif
c
c ++ time averaged generalized shear
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tashgc = 0.0
do kk=ks,ke
if (shgc(kk) .lt. rmiss) then
tashgc = tashgc + shgc(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tashgc = tashgc/ccount
else
tashgc = rmiss
endif
c
c ++ time averaged epos
klag=1
klead=1
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
taepos = 0.0
do kk=ks,ke
if (epos(kk) .lt. rmiss) then
taepos = taepos + epos(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
taepos = taepos/ccount
else
taepos = rmiss
endif
c
c ++ time averaged d200 and z850
klag=2
klead=2
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tad200 = 0.0
taz850 = 0.0
do kk=ks,ke
if (d200(kk) .lt. rmiss .and. z850(kk) .lt. rmiss) then
tad200 = tad200 + d200(kk)
taz850 = taz850 + z850(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tad200 = tad200/ccount
taz850 = taz850/ccount
else
tad200 = rmiss
taz850= rmiss
endif
c
c ++ time averaged t200, t250
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tat200 = 0.0
tat250 = 0.0
do kk=ks,ke
if (t200(kk) .lt. rmiss .and. t250(kk) .lt. rmiss) then
tat200 = tat200 + t200(kk)
tat250 = tat250 + t250(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tat200 = tat200/ccount
tat250 = tat250/ccount
else
tat200 = rmiss
tat250= rmiss
endif
c
c ++ time averaged g200
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tag200 = 0.0
do kk=ks,ke
if (g200(kk) .lt. rmiss) then
tag200 = tag200 + g200(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tag200 = tag200/ccount
else
tag200 = rmiss
endif
c
c ++ time averaged rlat
klag=4
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tarlat = 0.0
do kk=ks,ke
if (rlat(kk) .lt. rmiss) then
tarlat = tarlat + rlat(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tarlat = tarlat/ccount
else
tarlat = rmiss
endif
c
c ++ time averaged pw06
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tapw06 = 0.0
do kk=ks,ke
if (pw06(kk) .lt. rmiss) then
tapw06 = tapw06 + pw06(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tapw06 = tapw06/ccount
else
tapw06 = rmiss
endif
c
c ++ time averaged sst
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tarsst = 0.0
do kk=ks,ke
c if (sstl(kk) .lt. 200.0) then
c tarsst = tarsst + sstl(kk)
c ccount = ccount + 1.0
c endif
if (sst_lgem_local(kk) .lt. 200.0) then
tarsst = tarsst + sst_lgem_local(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tarsst = tarsst/ccount
else
tarsst = rmiss
endif
c
c ++ time averaged ohc
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
taohc = 0.0
do kk=ks,ke
c if (rhcn(kk) .lt. rmiss) then
c taohc = taohc + rhcn(kk)
c ccount = ccount + 1.0
c endif
if (ohc_lgem(kk) .lt. rmiss) then
taohc = taohc + ohc_lgem(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
taohc = taohc/ccount
c taohc = taohc-rthresh
taohc = taohc-ohc_thresh_lgem
if (taohc .lt. 0.0) taohc=0.0
else
taohc = rmiss
endif
c
c ++ time averaged 850-700 hPa temp advection
klag=5
klead=0
ks = k-klag
ke = k+klead
if (ks .lt. 0) ks = 0
if (ke .gt. mft) ke = mft
if (ke .lt. ks) ke = ks
ccount = 0.0
tatadv = 0.0
do kk=ks,ke
if (tadv(kk) .lt. rmiss) then
tatadv = tatadv + tadv(kk)
ccount = ccount + 1.0
endif
enddo
c
if (ccount .gt. 0.0) then
tatadv = tatadv/ccount
else
tatadv = rmiss
endif
c
c ++ Specify the values of the independent variables for cappa
var( 1) = cappa00
var( 2) = aday
var( 3) = spdx
var( 4) = pslv
var( 5) = tashdct
var( 6) = tashdctl
var( 7) = taepos
var( 8) = twacten
var( 9) = taz850
var(10) = tad200
var(11) = tat200
var(12) = tat250
var(13) = tarsst
var(14) = vmx
var(15) = taohc
var(16) = pc20
var(17) = tasdir
var(18) = tashgc
var(19) = tatadv
var(20) = pprobric2
var(21) = tag200
var(22) = tarlat
c var(22) = tapw06
c
varlab( 1) = 'CAP00'
varlab( 2) = 'ADAY '
varlab( 3) = 'SPDX '
varlab( 4) = 'PSLV '
varlab( 5) = 'SHDC '
varlab( 6) = 'LSHDC'
varlab( 7) = 'EPOS '
varlab( 8) = 'TWAT '
varlab( 9) = 'Z850 '
varlab(10) = 'D200 '
varlab(11) = 'T200 '
varlab(12) = 'T250P'
varlab(13) = 'SST '
varlab(14) = 'VMAX'
varlab(15) = 'NOHC '
varlab(16) = 'PC20 '
varlab(17) = 'SDIR '
varlab(18) = 'SHGC '
varlab(19) = 'TADV '
varlab(20) = 'PBRI '
varlab(21) = 'G200 '
varlab(22) = 'RLAT '
c varlab(22) = 'PW06 '
c
c ++ Normalize the independent variables
do j=1,nvarb
var(j) = (var(j) - xbar(j,k))/xsig(j,k)
enddo
c
c ++ calculate cappa
cappa(k+1) = 0.0
do j=1,nvarb
cappa(k+1) = cappa(k+1) + var(j)*coef(j,k)
enddo
cappa(k+1) = ybar(k) + cappa(k+1)*ysig(k)
c
if (rlat(k) .lt. rmiss) then
write(lulg,830) ktime,cappa(k+1)
830 format(/,'t= ',i4,' cappa=',f8.4)
c
do j=1,nvarb
write(lulg,831) j,varlab(j),var(j),coef(j,k),
+ var(j)*coef(j,k)
831 format(i3,1x,a5,1x,3(f8.4,1x))
c write(6,882) ktime,j,varlab(j),var(j)
c 882 format('t=',i3,' var=',i2,1x,a6,' value=',e11.4)
enddo
endif
c
c write(6,828) ktime,tashr,tashrl,tat200,tat250,tarsst,
c + tad200,taz850,taepos
c828 format('ktime,shr,l, t200,250, rsst,d200,z850,epos ',
c + i3,1x,10(f6.1,1x))
99 continue
c
c ++ Integrate LGE model ++
c
c Copy rlat,rlon,vsst_lgem to temporary arrays for the lgeim call
do k=0,mft
tlat(k+1) = rlat(k)
tlon(k+1) = rlon(k)
vmpi(k+1) = vsst_lgem(k)
enddo
c
call lgeim(nft,ftime,tlat,tlon,vmpi,cappa,vmx,
+ rmm,rnn,beta,vscale,lulg,vmaxl,dland,ierrl)
c
c ++ Check integration error
if (ierrl .ne. 0) then
vmaxl(1) = vmx
do i=2,nft
vmaxl(i) = 0.0
enddo
ierr=3
endif
c
c ++ Check for dissipation
dissloop: do i=2,nft
if (vmaxl(i) .lt. vdiss) then
do ii=i,nft
vmaxl(ii) = 0.0
enddo
exit dissloop
endif
enddo dissloop
c
c write(6,599) cappa00,aday,spdx,pslv
c 599 format('cappa00,aday,spdx,pslv: ',e11.4,1x,3(f6.1,1x))
c
c do i=1,nft
c write(6,600) ftime(i),tlat(i),tlon(i),sstl(i-1),
c + vmpi(i),shr(i-1),cappa(i),dland(i),vmaxl(i)
c 600 format(f6.0,1x,f6.1,1x,f7.1,1x,f6.1,1x,
c + f6.1,1x,f6.1,1x,e11.4,1x,f6.0,1x,f6.1)
c enddo
c
c Normal exit
write(lulg,400)
400 format(/,'LGE model completed normally')
return
c
c Error exit
900 continue
write(lulg,950) ierr
950 format(/,'Error in LGE model, no forecast made, ierr=',i2)
c
return
end
c
c ++++END LGEM MODULE++++
c
subroutine riipcset(ipc00,ipcm1,ipcm3,mft,imiss,rmiss,
+ rirpcri,aipc1ri,aipc2ri,dipc2ri)
c This routine calculates the IRPC variables for the RII
c routines using the new form of the IRPC input.
c
c Note that the IRPC info in ipc00,m1,m2 start at index 1, instead of
c index 0 as was the convention in the old IRPC variables. The 0 index
c in the new files contains the time of the image relative to the synoptic time.
c
dimension ipc00(0:mft),ipcm1(0:mft),ipcm3(0:mft)
c
c Initialize the output variables to missing.
c
rirpcri = rmiss
aipc1ri = rmiss
aipc2ri = rmiss
c
c The new IRPC files do not contain the PC tendency, so set that to zero,
c which is very close to the sample mean.
dipc2ri = 0.0
c
c Average as many of the times as are available
c
c PC1
ncount = 0
apc = 0.0
c
if (ipc00(1) .lt. imiss .and. ipc00(1) .gt. -900) then
apc = apc + float(ipc00(1))
ncount = ncount + 1
endif
c
if (ipcm1(1) .lt. imiss .and. ipcm1(1) .gt. -900) then
apc = apc + float(ipcm1(1))
ncount = ncount + 1
endif
c
if (ipcm3(1) .lt. imiss .and. ipcm3(1) .gt. -900) then
apc = apc + float(ipcm3(1))
ncount = ncount + 1
endif
c
if (ncount .gt. 0) then
aipc1ri = apc/float(ncount)
c rirpcri = aipc1ri/100.0
endif
c
c PC2
ncount = 0
apc = 0.0
c
if (ipc00(2) .lt. imiss .and. ipc00(2) .gt. -900) then
apc = apc + float(ipc00(2))
ncount = ncount + 1
endif
c
if (ipcm1(2) .lt. imiss .and. ipcm1(2) .gt. -900) then
apc = apc + float(ipcm1(2))
ncount = ncount + 1
endif
c
if (ipcm3(2) .lt. imiss .and. ipcm3(2) .gt. -900) then
apc = apc + float(ipcm3(2))
ncount = ncount + 1
endif
c
if (ncount .gt. 0) then
aipc2ri = apc/float(ncount)
rirpcri = aipc2ri/100.0
endif
c
return
end