module nnic_module ! ! v1.3.0 ! written by M. DeMaria June 2019 ! ! Modified 03/27/2023 by MD (v1.3.0) ! 1. Option to replace HWFI with HFAI in routine mread ! (controlled by ihswap set in mread) ! 2. Removed hard-coded model index values for filling 60 hr ! HWFI and AVNI forecasts ! 3. Header comment in reverse chronological order ! ! Modified 2/06/2022 by MD (v1.2.0) to help meet NCO standards ! ! Modified 4/28/2020 by SS (v1.1.1) to limit forecast values to 15kt ! ! Modified 4/19/2020 by MD and SS (v1.1.0) to add new option to fill in ! missing HWFI and other forecasts using ! the consensus tendency. Also, the ! intensity consensus is also written to ! the a-deck formated output file. ! ! Modified 04/13/2020 by SS (v1.0.5) to update for 7-day input ! to add TCLP & new lsdiag vars for 2020 ! to add ls0check, lstart subroutines ! ! Modified 07/02/2019 by MD (v1.0.4) to add writing track to output files ! ! This module defines variables and contains subroutines ! needed by the Neural Network Intensity Consensus (NNIC) model integer, parameter:: mxi=20 ! Max number of input nodes integer, parameter:: mxh=20 ! Max number of output nodes integer, parameter:: nft=14 ! Number of nnic forecast times integer, parameter:: idt=12 ! Time step (hr) for model forecasts integer, parameter:: nftl=28 ! Number of lsdiag.dat variable forecast times integer, parameter:: idtl=6 ! Time step (hr) for lsdiag.dat variables integer, parameter:: ntab=4 ! Number of hours to go back for time averaging lsdiag variables ! Models to include in the consensus integer, parameter:: ncm=4 ! Number of models integer, parameter:: nmin=2 ! Min number of model for a consensus forecast integer, parameter:: vmin=15.0 ! Min nn intensity forecast character(len=4),dimension(0:ncm):: mname ! Array for model names integer:: mreptype ! Method for replacing missing model forecasts ! =0 to replace with consensus model ! =0 to replace with consensus model tendency character(len=4) mnamecp ! CLIPER model name character(len=4) mnamecpe! Model name to extend using mnamecp (usually HWFI) ! Note: cliper model extension ! skipped if mreptype=1 ! Model input variables integer,dimension(0:nft) :: mftime !Model forecast time array integer,dimension(0:nft) :: ncon !No. of models available for the consensus real,dimension(0:ncm,0:nft):: vmaxf !Model intensity forecast array real,dimension(0:nft) :: vmaxfcp !Model intensity forecast array for CLIPER real,dimension(0:nft) :: vmaxfha !Model intensity forecast array for HFAI real,dimension(0:nft) :: vmaxfhb !Model intensity forecast array for HFBI real,dimension(0:ncm,0:nft):: vmaxd !Deviations from equally weighted consensus real,dimension(0:nft):: vmaxnn !Unequally weighted NN forecast real,dimension(0:nft):: vmaxnc !Equally weighted consensus of NN input models integer,dimension(0:nft):: ivmaxnn !Integer version of vmaxnn integer,dimension(0:nft):: ivmaxnc !Integer version of vmaxnc integer mindexcp !Index for cliper model integer mindexcpe !Index for model to extend integer kindexcpe !Time index to start extension of mnamecpe real,parameter:: vmaxlo= 10.0 !Minimum valid max wind real,parameter:: vmaxhi=185.0 !Maximum valid max wind ! I/O variables character(len=30),parameter:: fnminp='nn_input.dat' !Input model forecasts character(len=30),parameter:: fnlinp='lsdiag.dat' !Input lsdiag.dat file character(len=30) fncoef !Input file with NN weight and normalization factors ! Note: fncoef name set in routine cread character(len=30),parameter:: fnlog ='nnic.log' !Output log file character(len=30),parameter:: fndat ='nnic.dat' !Output forecast file (ATCF format) character(len=30),parameter:: fntxt ='nnic.txt' !Output forecast file (old ATCF format) integer,parameter:: runid=28 !Model ID number for old ATCF file integer,dimension(0:nft):: latdum,londum !Dummy arrays for lat/lon for old ATCF file integer,parameter:: luminp=10 integer,parameter:: lulinp=11 integer,parameter:: lucoef=12 integer,parameter:: lulog =13 integer,parameter:: ludat =14 integer,parameter:: lutxt =15 ! Variables for lsdiag.dat file integer:: ivmax0,iyr2l,imonl,idayl,itimel integer:: iyr4l,idtgl character(len=8):: atcfidl real:: vmax0,dvm12 real,dimension(0:nftl):: slat,slon,dtl,shdc,sstn real,dimension(0:nftl):: slatta,slonta,sstnta,shdcta !Time-averaged lsdiag.dat variables integer,parameter:: imiss=9999 real,parameter:: rmiss=9999. ! Variables for NN calculation integer:: nftc,idtc,nxi,nxh,nxo character(len=4):: xalab(mxi) real:: sigmani,sigmano real,dimension(mxi,nft):: xa,xabar,xasig real,dimension( nft):: ya,yabar,yasig real,dimension(mxi,mxh,nft):: w real,dimension(mxh,nft):: x real,dimension(mxh,nft):: rn !Input to hidden nodes real,dimension(mxh,nft):: rh !Output from hidden nodes real,dimension(nft):: ra !Input to output nodes real,dimension(nft):: ro !Output from output nodes contains subroutine ioset ! This routine opens the I/O files with generic names implicit none open(file=fnlog, unit=lulog ,status='replace') open(file=fnminp,unit=luminp,status='old',err=900) open(file=fnlinp,unit=lulinp,status='old',err=900) open(file=fndat, unit=ludat ,status='replace',err=900) open(file=fntxt, unit=lutxt ,status='replace',err=900) return 900 continue write(lulog,950) 950 format(/,'Error opening files in ioset, fatal error') stop end subroutine ioset subroutine lread ! This routine reads the lsdiag.dat file implicit none ! Local variables character(len=4):: tlab character(len=180):: tline integer:: i,j,k integer:: mxll !Max lines in an lsdiag.dat file integer:: idvm12 integer,dimension(0:100):: itemp ! Initialize input variables to missing dvm12=rmiss do k=0,nftl slat(k) = rmiss slon(k) = rmiss dtl(k) = rmiss shdc(k) = rmiss sstn(k) = rmiss enddo ! Read lsdiag.dat file header line tline='XXXX' read(lulinp,100,end=1000,err=1000) tline 100 format(a180) 1000 continue i=1 if (tline(157:160) .ne. 'HEAD') then write(lulog,900) tline 900 format(/,'Fatal error reading lsdiag.dat, expecting header line',/,' tline=',a180) endif ! Extract info from the header line read(tline,200) iyr2l,imonl,idayl,itimel,ivmax0,atcfidl 200 format(6x,3(i2),1x,i2,1x,i4,20x,a8) ! Convert header line info vmax0 = float(ivmax0) if (iyr2l .lt. 50) then iyr4l = iyr2l + 2000 else iyr4l = iyr2l + 1900 endif idtgl = 1000000*iyr4l + 10000*imonl + 100*idayl + itimel ! Read interior lines of lsdiag.dat file mxll = 1000 lsloop: do j=1,mxll read(lulinp,100,end=3000) tline i = i + 1 if (tline(157:160) .eq. 'DELV') then read(tline,105) idvm12 105 format(i5) if (idvm12 .le. imiss) then dvm12 = float(idvm12) else dvm12 = 0.0 endif endif if (tline(157:160) .eq. 'LAT ') then read(tline,110) (itemp(k),k=0,nftl) 110 format(10x,29(i5)) do k=0,nftl if (itemp(k) .eq. imiss) then slat(k) = rmiss else slat(k) = 0.1*float(itemp(k)) endif enddo endif if (tline(157:160) .eq. 'LON ') then read(tline,110) (itemp(k),k=0,nftl) do k=0,nftl if (itemp(k) .eq. imiss) then slon(k) = rmiss else slon(k) = -0.1*float(itemp(k)) endif enddo endif if (tline(157:160) .eq. 'DTL ') then read(tline,110) (itemp(k),k=0,nftl) do k=0,nftl if (itemp(k) .eq. imiss) then dtl(k) = rmiss else dtl(k) = float(itemp(k)) endif enddo endif if (tline(157:160) .eq. 'SHDC') then read(tline,110) (itemp(k),k=0,nftl) do k=0,nftl if (itemp(k) .eq. imiss) then shdc(k) = rmiss else shdc(k) = 0.1*float(itemp(k)) endif enddo endif if (tline(157:160) .eq. 'NSTA') then read(tline,110) (itemp(k),k=0,nftl) do k=0,nftl if (itemp(k) .eq. imiss) then sstn(k) = rmiss else sstn(k) = 0.1*float(itemp(k)) endif enddo endif if (tline(157:160) .eq. 'LAST') then exit lsloop endif enddo lsloop 3000 continue ! Check to make sure at least the t=0 lsdiag values are available call ls0check(sstn(0),rmiss,atcfidl,imonl,idayl,itimel,lulog) call ls0check(shdc(0),rmiss,atcfidl,imonl,idayl,itimel,lulog) ! Fill in missing lsdiag.dat data call patch0(slat,rmiss,nftl) call patch0(slon,rmiss,nftl) call patch0(shdc,rmiss,nftl) call patch0(sstn,rmiss,nftl) ! Time average lsdiag inputs call lstart(nftl,nftl,slat,slatta) call lstart(nftl,nftl,slon,slonta) call lstart(nftl,nftl,shdc,shdcta) call lstart(nftl,nftl,sstn,sstnta) write(lulog,300) i 300 format('lsdiag.dat read completed, no. lines=',i4) write(lulog,310) atcfidl,iyr4l,imonl,idayl,itimel,vmax0,dvm12 310 format(/,'lsdiag.dat case info: ATCFID=',a8,2x,i4.4,3(i2.2),' v0=',f5.0,' DVM12=',f5.0) write(lulog,315) (k*idtl,k=0,nftl) 315 format('TIME ',29(3x,i4)) tlab = 'SLAT' write(lulog,320) tlab,(slat(k),k=0,nftl) tlab = 'SLON' write(lulog,320) tlab,(slon(k),k=0,nftl) tlab = 'DTL ' write(lulog,320) tlab,(dtl(k),k=0,nftl) tlab = 'SHDC' write(lulog,320) tlab,(shdc(k),k=0,nftl) tlab = 'NSTA' write(lulog,320) tlab,(sstn(k),k=0,nftl) 320 format(a4,1x,29(f7.1)) return end subroutine lread subroutine mselect ! This routine defines the input models implicit none ! Local variables integer:: m mname(0) = 'NNIB' mname(1) = 'DSHP' mname(2) = 'LGEM' mname(3) = 'HWFI' mname(4) = 'AVNI' ! Specify the method for replacing missing models ! Note: mreptype=0 or 1 to use consensus or consensus tendency mreptype=1 ! CLIPER model anem mnamecp = 'TCLP' ! Specific model name to use CLIPER to extend the forecast length mnamecpe = 'HWFI' ! Specify time index to start extending mnamecpe kindexcpe = 12 write(lulog,200) ncm,(mname(m),m=1,ncm) 200 format(/,i4,' models included: ',10(a4,1x)) if (mreptype .eq. 0) then write(lulog,201) mreptype 201 format(' mreptype=',i1,' missing model forecasts replaced by consensus',/) else write(lulog,202) mreptype 202 format(' mreptype=',i1,' missing model forecasts replaced by consensus tendency',/) endif end subroutine mselect subroutine mread ! This routine reads the models to be used in the consensus forecast implicit none ! Local variables integer:: i,j,k,kk,l,m,mm,n,tt integer:: idum1 integer:: isnum,idtg,mftimet,ivmaxt integer:: mxml !Max number of model lines character(len=2):: cbasin,cdum1 character(len=4):: mnamet character(len=5):: clat,clon integer :: mindx_hwfi !Model index for HWFI integer :: mindx_avni !Model index for AVNI integer :: ihswap !Parameter for swapping HWFI with HFAI or HFBI character(len=4) :: swap_name !Tech ID of HWFI replacement ! Set ihswap = 1 to swap HWFI with HFAI when available ! or ihswap = 2 to swap HWFI with HFBI when available ! or ihswap = 0 to leave HWFI alone ihswap = 1 if (ihswap .eq. 1) then swap_name='HFAI' elseif (ihswap .eq. 2) then swap_name='HFBI' else swap_name='XXXX' endif if (ihswap .gt. 0) then write(lulog,301) swap_name 301 format(/,'HWFI replaced by ',a4,' when possible') else write(lulog,303) 303 format(/,'HWFI not replaced by ',a4) endif ! Calculate time array do k=1,nft mftime(k) = k*idt enddo ! Initialize model forecasts to missing do m=0,ncm do k=0,nft vmaxf(m,k) = rmiss enddo enddo do k=0,nft vmaxfcp(k) = rmiss vmaxfha(k) = rmiss vmaxfhb(k) = rmiss enddo ! Read model file i=0 mindexcpe = -99 mxml=1000 modloop: do j=1,mxml read(luminp,*,end=2000,err=2000) cbasin,isnum,idtg,cdum1,mnamet,mftimet,clat,clon,ivmaxt i = i + 1 ! Search for required model forecasts if (idtg .ne. idtgl) cycle modloop if (cbasin .ne. atcfidl(1:2)) cycle modloop kk = -99 do k=0,nft if (mftimet .eq. mftime(k)) then kk = k endif enddo if (kk .lt. 0) cycle modloop mm = -99 do m=1,ncm if (mnamet .eq. mname(m)) then mm = m ! Identify model (HWFI) index to replace with TCLP if (mnamet .eq. mnamecpe) mindexcpe = m endif enddo ! Save HFAI model forecast if (mnamet .eq. 'HFAI') then vmaxfha(kk) = float(ivmaxt) endif ! Save HFBI model forecast if (mnamet .eq. 'HFBI') then vmaxfhb(kk) = float(ivmaxt) endif ! Save TCLP model forecast if (mnamet .eq. mnamecp) then vmaxfcp(kk) = float(ivmaxt) endif if (mm .lt. 0) cycle modloop ! If you get to here, a required model forecast is found, so save it vmaxf(mm,kk) = float(ivmaxt) if (vmaxf(mm,kk) .lt. vmin) vmaxf(mm,kk) = rmiss enddo modloop 2000 continue write(lulog,200) i,fnminp 200 format(i4,' lines read from model input file ',a30) ! Find model index of HWFI mindx_hwfi = -99 do m=1,ncm if (mname(m) .eq. 'HWFI') then mindx_hwfi = m endif enddo if (ihswap .eq. 1) then if (vmaxfha(0) .lt. rmiss .and. vmaxfha(0) .gt. 0.0) then write(lulog,305) (vmaxf(mindx_hwfi,k),k=0,nft) 305 format(/,'HWFI before replacement: ',31(f6.1,1x)) do k=0,nft vmaxf(mindx_hwfi,k) = vmaxfha(k) enddo write(lulog,306) (vmaxf(mindx_hwfi,k),k=0,nft) 306 format( 'HWFI after replacement: ',31(f6.1,1x)) endif elseif (ihswap .eq. 2) then if (vmaxfhb(0) .lt. rmiss .and. vmaxfhb(0) .gt. 0.0) then write(lulog,305) (vmaxf(mindx_hwfi,k),k=0,nft) do k=0,nft vmaxf(mindx_hwfi,k) = vmaxfha(k) enddo write(lulog,306) (vmaxf(mindx_hwfi,k),k=0,nft) endif endif ! Temporary repair of 60 hr point ! HWFI if (mindx_hwfi .gt. 0) then if (vmaxf(mindx_hwfi,4) .lt. rmiss .and. & vmaxf(mindx_hwfi,6) .lt. rmiss .and. & vmaxf(mindx_hwfi,5) .ge. rmiss) then vmaxf(mindx_hwfi,5) = 0.5*(vmaxf(mindx_hwfi,4)+vmaxf(mindx_hwfi,6)) endif endif ! AVNI mindx_avni = -99 do m=1,ncm if (mname(m) .eq. 'AVNI') then mindx_avni = m endif enddo if (vmaxf(mindx_avni,4) .lt. rmiss .and. & vmaxf(mindx_avni,6) .lt. rmiss .and. & vmaxf(mindx_avni,5) .ge. rmiss) then vmaxf(mindx_avni,5) = 0.5*(vmaxf(mindx_avni,4)+vmaxf(mindx_avni,6)) endif ! Extend selected forecast using cliper if (kindexcpe .gt. 0 .and. kindexcpe .le. nft .and. mreptype .eq. 0) then do k=kindexcpe,nft if (vmaxf(mindexcpe,kindexcpe-1) .gt. vmaxlo .and.& &vmaxf(mindexcpe,kindexcpe-1) .le. vmaxhi .and.& &vmaxf(mindexcpe,k) .ge. rmiss ) then vmaxf(mindexcpe,k) = vmaxfcp(k) endif enddo endif write(lulog,205) (mftime(k),k=0,nft) 205 format(/,'Input model max wind forecasts',/,'TIME ',29(3x,i4)) do m=1,ncm write(lulog,210) mname(m),(vmaxf(m,k),k=0,nft) 210 format(a4,2x,29(f7.0)) enddo mnamet = 'HFAI' write(lulog,210) mnamet,(vmaxfha(k),k=0,nft) mnamet = 'HFBI' write(lulog,210) mnamet,(vmaxfhb(k),k=0,nft) write(lulog,211) mnamecpe, mnamecp, (vmaxfcp(k),k=0,nft) 211 format(/,'CLIPER model used to fill in long range ',a4,' if mreptype=0 ',/,a4,2x,29(f7.0)) return end subroutine mread subroutine cread(atcfidc) ! This routine reads the file with the NN coefficients in it implicit none ! Calling arguments character(len=8):: atcfidc ! Local variables integer i,j,k,l,m,n integer ktime,idum1 character(len=2) cbasin2 character(len=10) cdum1 character(len=256) fnco,coef_location ! Create coefficent file name fncoef = 'nnfit_coef_bb.dat' fncoef(12:13) = atcfidc(1:2) call getenv( "SHIPS_COEF",coef_location ) fnco = trim( coef_location) // fncoef ! Open the coefficient file open(file=fnco,unit=lucoef,form='formatted',status='old',err=900) write(lulog,200) fnco(1:100) 200 format(/,'NN coef file opened: ',a100) ! Read the header line read(lucoef,*,end=910,err=910) cbasin2,nftc,idtc,nxi,nxh,nxo,sigmani,sigmano ! Check forecast time info for consistency if (nftc .ne. nft .or. idtc .ne. idt) then write(lulog,970) 970 format(/,'Inconsistent time info in NN coef file: ',a100) stop endif write(lulog,210) nxi,nxh 210 format(i2,' input nodes, ',i2,' hidden nodes') ! Read the NN parameters for each forecast hour do k=1,nft read(lucoef,*,err=910,end=910) ktime read(lucoef,*,err=910,end=910) idum1,yabar(k),yasig(k) do i=1,nxi read(lucoef,*,err=910,end=910) idum1,xabar(i,k),xasig(i,k),xalab(i) enddo read(lucoef,*,err=910,end=910) cdum1 do i=1,nxi read(lucoef,*,err=910,end=910) (w(i,j,k),j=1,nxh) enddo read(lucoef,*,err=910,end=910) cdum1 read(lucoef,*,err=910,end=910) (x(j,k),j=1,nxh) enddo do i=1,nxi write(lulog,220) i,xalab(i) 220 format('Input ',i2,2x,a4) enddo close(lucoef) return 900 continue write(lulog,950) fnco 950 format(/,'Error opening NN coef file: ',a100) stop 910 continue write(lulog,960) fnco 960 format(/,'Error reading NN coef file: ',a100) stop end subroutine cread subroutine concal ! This routine calculates the equally weighted consensus ! and the deviations from the consensus implicit none ! Local variables integer:: i,j,k,l,m,mm,ncont,kconlast real:: vmaxc,dvcon ! Calculate consensus forecast and put it in element m=0 of vmaxf m=0 k=0 vmaxf(m,k) = vmax0 ncon(k) = ncm do k=1,nft ncont = 0 vmaxc = 0.0 do m=1,ncm if (vmaxf(m,k) .lt. rmiss) then vmaxc = vmaxc + vmaxf(m,k) ncont = ncont + 1 endif enddo if (ncont .ge. nmin) then vmaxf(0,k) = vmaxc/float(ncont) ncon(k) = ncont else vmaxf(0,k) = rmiss ncon(k) = 0 endif enddo if (mreptype .eq. 1) then ! Fill in missing model forecasts with consensus tendency ! Check t=0 hr model values first to make sure there is at least one good point do m=1,ncm if (vmaxf(0,0) .lt. rmiss .and. vmaxf(m,0) .ge. rmiss) then vmaxf(m,0) =vmaxf(0,0) endif enddo ! Find last forecast time with valid consesnus forecasts up to that time kconlast=nft do k=1,nft if (vmaxf(0,k) .ge. rmiss) then kconlast = k-1 exit endif enddo ! if (kconlast .ge. 1) then ! Check for missing models and fill in with consensus tendency do k=1,kconlast dvcon = (vmaxf(0,k)-vmaxf(0,k-1)) do m=1,ncm if (vmaxf(m,k) .ge. rmiss) then vmaxf(m,k) = vmaxf(m,k-1) + dvcon endif enddo enddo endif write(lulog,201) mftime(kconlast) 201 format(/,'Last time with valid consensus forecast=',i4,' hr') write(lulog,202) (mftime(k),k=0,nft) 202 format(/,'Corr. model max wind forecasts',/,'TIME ',29(3x,i4)) do m=1,ncm write(lulog,203) mname(m),(vmaxf(m,k),k=0,nft) 203 format(a4,2x,29(f7.0)) enddo endif ! Calculate deviations from equally weighted consensus do k=0,nft do m=0,ncm if (vmaxf(0,k) .lt. rmiss .and. vmaxf(m,k) .lt. rmiss) then vmaxd(m,k) = vmaxf(m,k) - vmaxf(0,k) else vmaxd(m,k) = 0.0 endif enddo enddo ! Write equally weighted consensus forecast to the log file write(lulog,205) (mftime(k),k=0,nft) 205 format(/,'Equally weighted consensus forecasts',/,'TIME ',29(3x,i4)) write(lulog,210) mname(0),(vmaxf(0,k),k=0,nft) 210 format(a4,2x,29(f7.0)) write(lulog,215) (ncon(k),k=0,nft) 215 format('NCON ',29(3x,i4)) ! Write deviation model forecasts to the log file write(lulog,220) (mftime(k),k=0,nft) 220 format(/,'Deviations from consensus',/,'TIME ',29(3x,i4)) do m=0,ncm write(lulog,225) mname(m),(vmaxd(m,k),k=0,nft) 225 format(a4,2x,29(f7.0)) enddo return end subroutine concal subroutine nniccal ! This routine uses the neural net to calculate the unequally ! weighted consensus intensity forecast implicit none ! Local variables integer:: i,j,k,l,m,n real ftemp ! Specify the normalized model input variables do i=1,ncm do k=1,nft if (vmaxf(0,k) .lt. rmiss) then xa(i,k) = (vmaxd(i,k)-xabar(i,k))/xasig(i,k) else xa(i,k) = 0.0 endif enddo enddo ! Specify the normalized non-model input variables if (nxi .gt. ncm) then do k=1,nft if (vmaxf(0,k) .lt. rmiss) then xa(ncm+1,k) = (vmaxf(0,k) -xabar(ncm+1,k))/xasig(ncm+1,k) xa(ncm+2,k) = (dvm12 -xabar(ncm+2,k))/xasig(ncm+2,k) xa(ncm+3,k) = (slatta(2*k)-xabar(ncm+3,k))/xasig(ncm+3,k) xa(ncm+4,k) = (sstnta(2*k)-xabar(ncm+4,k))/xasig(ncm+4,k) xa(ncm+5,k) = (shdcta(2*k)-xabar(ncm+5,k))/xasig(ncm+5,k) ! write(lulog,*) 'k,slat(2k): ',k,slat(2*k) else xa(ncm+1,k) = 0.0 xa(ncm+2,k) = 0.0 xa(ncm+3,k) = 0.0 xa(ncm+4,k) = 0.0 xa(ncm+5,k) = 0.0 endif enddo endif ! Scale the input variables to lie between 0 and 1 do i=1,nxi do k=1,nft ftemp = 0.5 + xa(i,k)/(2.0*sigmani) if (ftemp .gt. 1.0) ftemp = 1.0 if (ftemp .lt. 0.0) ftemp = 0.0 xa(i,k) = ftemp enddo enddo ! Calculate the input to the hidden nodes do k=1,nft do j=1,nxh rn(j,k) = 0.0 do i=1,nxi rn(j,k) = rn(j,k) + w(i,j,k)*xa(i,k) enddo enddo enddo ! Calcuate output from hidden nodes do k=1,nft do j=1,nxh rh(j,k) = tfunc(rn(j,k)) enddo enddo ! Calculate input to output nodes do k=1,nft ra(k) = 0.0 do j=1,nxh ra(k) = ra(k) + x(j,k)*rh(j,k) enddo enddo ! Calculate output from output nodes do k=1,nft ya(k) = tfunc(ra(k)) enddo write(lulog,870) 870 format(10x) do k=1,nft write(lulog,880) k*12,(xa(i,k),i=1,nxi) 880 format(i3,' hr NN input ',28(f8.3,1x)) enddo ! do k=1,nft ! write(lulog,860) k*12 ! 860 format(/,'w(i,j), t=',i3,' hr') ! do i=1,nxi ! write(lulog,861) (w(i,j,k),j=1,nxh) ! 861 format(20(f8.3,1x)) ! enddo ! write(lulog,*) 'rn(j)' ! write(lulog,861) (rn(j,k),j=1,nxh) ! write(lulog,*) 'rh(j)' ! write(lulog,861) (rh(j,k),j=1,nxh) ! write(lulog,*) 'x(j)' ! write(lulog,861) (x(j,k),j=1,nxh) ! write(lulog,866) ra(k),ya(k) ! 866 format('ra=',f8.3,' ya=',f8.3) ! enddo write(lulog,883) (k*12,k=1,nft) 883 format(/,'TIME ',28(6x,i3)) write(lulog,885) (ya(k),k=1,nft) 885 format('NN output ',28(f8.3,1x)) ! Calculate NNIC forecast by un-normalizing output from the neural network. ! Also, save un-weighted consensus of NN input models vmaxnc(0) = vmax0 vmaxnn(0) = vmax0 ivmaxnc(0) = nint(vmax0) ivmaxnn(0) = nint(vmax0) do k=1,nft if (vmaxf(0,k) .lt. rmiss) then vmaxnc(k) = vmaxf(0,k) ftemp = 2.0*sigmano*(ya(k)-0.5) vmaxnn(k) = vmaxnc(k) + yasig(k)*ftemp + yabar(k) ! vmaxnn(k) = vmaxf(0,k) else vmaxnc(k) = 0.0 vmaxnn(k) = 0.0 endif ! if (vmaxnn(k) .ge. 5.0 .and. vmaxnn(k) .lt. 15.0) vmaxnn(k)=15.0 if (vmaxnn(k) .lt. 15.0) vmaxnn(k)=0.0 if (vmaxnc(k) .lt. 15.0) vmaxnc(k)=0.0 ivmaxnc(k) = nint(vmaxnc(k)) ivmaxnn(k) = nint(vmaxnn(k)) enddo ! Write NN consensus forecast to the log file write(lulog,205) (mftime(k),k=0,nft) 205 format(/,'Neural Net Intensity Consensus forecasts',/,'TIME ',29(3x,i4)) write(lulog,211) (vmaxnc(k),k=0,nft) 211 format('NNIB',2x,29(f7.0)) write(lulog,210) (vmaxnn(k),k=0,nft) 210 format('NNIC',2x,29(f7.0)) return end subroutine nniccal subroutine write_nnicoa ! This routine writes the NNIC forecast in the old ATCF format implicit none ! Local variables integer:: i,j,k,l,m,n ! Fill integer lat/lon array for printing do k=0,nft if (slat(k) .lt. rmiss) then latdum(k) = nint(10.0*slat(k)) else latdum(k) = 0 endif if (slon(k) .lt. rmiss) then londum(k) = nint(-10.0*slon(k)) else londum(k) = 0 endif enddo write(lutxt,200) runid,iyr2l,imonl,idayl,itimel,(latdum(k),londum(k),k=1,nft),& &(ivmaxnc(k),k=1,nft),atcfidl(1:4),atcfidl(7:8) 200 format('92NC',5(i2.2),28(i4),14(i3),1x,a4,a2) write(lutxt,210) runid,iyr2l,imonl,idayl,itimel,(latdum(k),londum(k),k=1,nft),& &(ivmaxnn(k),k=1,nft),atcfidl(1:4),atcfidl(7:8) 210 format('92NN',5(i2.2),28(i4),14(i3),1x,a4,a2) return end subroutine write_nnicoa subroutine write_nnicna ! This routine writes the NNIC forecast in ATCF format implicit none ! Local variables integer:: i,j,k,l,m,n,minc integer:: iftime(0:100) character(len=4) :: aidlab character(len=4),dimension(0:100):: stype character(len=3):: rdum character(len=10):: aymdh ! Set time interval for writing forecast minc=1 ! Forecast time array do k=0,nft iftime(k) = k*idt enddo ! Set storm type to default tropical do k=0,nft stype(k) = 'TROP' enddo ! Copy date/time group to character format write(aymdh,100) idtgl 100 format(i10) ! Fill integer lat/lon array for printing do k=0,nft if (slat(k) .lt. rmiss) then latdum(k) = nint(10.0*slat(k)) else latdum(k) = 0 endif if (slon(k) .lt. rmiss) then londum(k) = nint(10.0*slon(k)) else londum(k) = 0 endif enddo ! Write NN baseline consensus forecast to ATCF file aidlab='NNIB' call writeaidlocal2(ludat,atcfidl,aymdh,aidlab,& &latdum,londum,ivmaxnc,stype,iftime,nft,minc) ! Write NN forecast to ATCF file aidlab='NNIC' call writeaidlocal2(ludat,atcfidl,aymdh,aidlab,& &latdum,londum,ivmaxnn,stype,iftime,nft,minc) return end subroutine write_nnicna subroutine ls0check(var0,rmiss0,atcfid0,imon0,iday0,itime0,lulog0) ! This routine checks to make sure var0 is not missing ! Passed variables integer:: imon0,iday0,itime0,lulog0 character(len=8):: atcfid0 real:: var0,rmiss0 if (var0 .ge. rmiss0) then write(lulog0,980) var0,atcfid0,imon0,iday0,itime0 980 format(/,'Missing t=0 lsdiag.dat, val=',f9.2,2x,a8,1x,2(i2.2),1x,i2& &/,'Stop processing') stop endif return end subroutine ls0check subroutine lstart(nftmla,nftma,varlc,varta) ! This routine time averages the lsdiag variables. The time interval ! of the input lsdiag array (varlc) can be different than the time-averaged ! output lsdiag array (varta). This is a real-time version. implicit none ! Calling arguments integer:: nftmla, nftma real,dimension(0:nftmla):: varlc real,dimension(0:nftma):: varta ! Local variables integer:: j,k,kk,m,n integer:: ktime,kktime integer:: kklo,kkhi real:: varsum,vcount ! Copy t=0 values varta(0) = varlc(0) do k=1,nftma ktime = idtl*k kkhi = ktime/idtl kklo = kkhi - ntab if (kklo .lt. 0) kklo = 0 varsum = 0.0 vcount = 0.0 do kk=kklo,kkhi if (varlc(kk) .lt. rmiss) then varsum = varsum + varlc(kk) vcount = vcount + 1.0 endif enddo if (vcount .ge. 1.0) then varta(k) = varsum/vcount else varta(k) = rmiss endif enddo end subroutine lstart function tfunc(a) ! This routine evaluates the transfer function used in the NN implicit none ! Calling arguments real:: a,tfunc tfunc = 1.0/(1.0 + exp(-a)) return end function tfunc end module nnic_module