subroutine tab_model(rlon00,rlat00,rlonm12,rlatm12,iftype, + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt, + glon,glat,gp,gt,ugrid,vgrid,pbot,ptop, + flon,flat,lulog,ierr) c c v3.1.0 for 2018 hurricane season c v3.1.1 for 2019 hurricane season c c This routine makes a track forecast by following a trajectory c where the horizontal motion is a weighted sum of persistence c and a motion vector from the vertically and horizontally averaged flow c from the input 3-D u,v velocity field (usually from the GFS). c c This version has the bottom and top pressure levels as input variables. c C Modifications: c v3.0.0 - pninit constants from DAZ c v3.1.0 - Increased averaging radius in pninit to 450 km c to reduce track oscillations. MD, 11/13/2017. c v3.1.1 - Increased averaging radius in pninit to 600 km c SS, 05/16/2019. c v3.1.2 - NCO compliance. KM, 02/02/2023 c v3.1.3 - Time integration scheme rearranged to get rid of c last remaining go to statement. Function idecod c removed from end of file. It is already contained c in the generalutils library. JFD 02/06/2023 c c Input: rlon00 - Longitude at t=0 (0 to 360 deg) c rlat00 - Latitude at t=0 (-90 to 90 deg) c rlonm12 - same as rlon00 at t=-12 hr (not used if ipf=1) c rlatm12 - same as rlat 00 at t=-12 hr (not used if ipf=1) c iftype - Flag for combining persistence and GFS motion vectors c =0 for persistence, =1 for GFS, =2 for time weighted mean of each c mxx - Maximum number of longitudes c mxy - Maximum number of latitudes c mxp - Maximum number of pressure levels c mxt - Maximum number of forecast times c c nxx, nxy, nxp, nxt - Actual number of longitudes, latitudes, etc c glon - Longitudes of the 3-D wind field (0 to 360 deg) c glat - Latitudes of the 3-D wind field (-90 to 90 deg) c gp - Pressure levels of the 3-D wind field (hPa) ordered from low c to high pressure. gp should be set to -9999. if data at a pressure c level is missing. c gt - Times (hr, 0 to nxt) of the 3-D wind field c ugrid - u (m/s) as a function of lon,lat,p,t c vgrid - v (m/s) as a function of lon,lat,p,t c pbot - Bottom pressure level (hPa) for the vertical average c ptop - Top pressure level (hPa) for the vertical average c c Output: flon(0:nxt) - Forecast longitudes (0 to 360 deg) c flat(0:nxt) - Forecast latitudes (-90 to 90 deg) c lulog - Unit number for writing messages c ierr - Error variable (=0 for normal completion) c (=1 if not enough GFS data for a forecast) c (=2 if problems with pressure levels) c c Passed variables dimension glon(mxx),glat(mxy),gp(mxp),gt(0:mxt) dimension ugrid(mxx,mxy,mxp,0:mxt),vgrid(mxx,mxy,mxp,0:mxt) dimension flat(0:mxt),flon(0:mxt),fvmax(0:mxt) c c Local variables parameter (mdt=3000,mxpt=100) dimension flons(0:mdt),flats(0:mdt),vmaxs(0:mdt) dimension cx(0:mdt),cy(0:mdt) dimension cxc(0:mdt),cyc(0:mdt) dimension cxp(0:mdt),cyp(0:mdt) dimension wxc(0:mdt),wyc(0:mdt) dimension wxp(0:mdt),wyp(0:mdt) dimension tf(0:mdt) dimension wcc(0:mdt),wcp(0:mdt),vwts(mxpt) c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /pncon2/ auc,avc,efup,efvp,aup,avp,betax,betay,bcfac common /pncon3/ boxhx,boxhy,boxhr c ierr = 0 c c Check to make sure input fields are available at a minimum of 2 times if (nxt .lt. 1) then ierr=1 return endif c c Set the maximum time step (hr) to use for the integration dtmax = 1.0 c c Calculate the time step of the input u,v wind field and forecast times dthr = gt(1)-gt(0) c Set debug write flag idbug = 0 c c Calcuate vertical weights call vwcal(gp,mxp,nxp,mxpt,pbot,ptop,vwts,ierr) if (ierr .ne. 0) return c c Initialize forecast to zero ndt = nxt do k=0,mxt flat(k) = 0.0 flon(k) = 0.0 enddo c do i=0,mdt flons(i) = 0.0 flats(i) = 0.0 enddo c c Initialize numerical and physical constants call pninit c c Write newbam2 parameters to output file efupi = 1.0/efup efvpi = 1.0/efvp write(lulog,300) iftype,pbot,ptop,boxhx,boxhy,boxhr,betax,betay, + bcfac,aup,avp,efupi,efvpi 300 format(/,'TAB forecast with iftype=',i1,/, + 'Vertical average pbot,ptop: ',f6.0,1x,f6.0,/, + 'Horizontal average x,y half widths:',f6.0,1x,f6.0,/, + 'Radius for case with circular area:',f6.0 ,/, + 'Beta drift, x,y and cosine factor: ',3(f6.2,1x) ,/, + 'Persistence weights, x,y: ',(f6.2,1x,f6.2),/, + 'Persistence e-foldng times, x,y: ',(f6.0,1x,f6.0)) write(lulog,302) glon(1),glon(nxx),glon(2)-glon(1), + glat(1),glat(nxy),glat(2)-glat(1) 302 format('Lon domain, dlon: ',3(f6.1,1x),/, + 'Lat domain, dlat: ',3(f6.1,1x)) c c Initialize the position flon(0) = rlon00 flat(0) = rlat00 flons(0) = rlon00 flats(0) = rlat00 c c Calculate initial storm motion vector cfac = cos(0.5*dtr*(rlat00+rlatm12)) dlon00 = (rlon00-rlonm12)*cfac dlat00 = (rlat00-rlatm12) c c Check to see if Greenwich meridian was crossed if (dlon00 .gt. 180.0) dlon00 = dlon00 - 360.0 if (dlon00 .lt. -180.0) dlon00 = dlon00 + 360.0 c cxp(0) = (dlon00*deg2m)/(12*3600.0) cyp(0) = (dlat00*deg2m)/(12*3600.0) c c Compare dtmax and output time step and compute the c interval for saving the output (intsave) and the c total number of time steps to integrate the model if (dthr .lt. dtmax) then intsave = 1 ndts = ndt dts = dthr else intsave = ifix(0.1 + dthr/dtmax) ndts = ndt*intsave dts = dthr/float(intsave) endif dtss = 3600.0*dts c c Calcuate time (hr) at each time step, c and persistence vector at each time step tf(0) = 0.0 c do k=1,ndts tf(k) = tf(0) + dts*float(k) c cxp(k) = cxp(0) cyp(k) = cyp(0) enddo c if (idbug .eq. 1) then write(lulog,700) cxp(0),cyp(0),dthr,dts,ndts 700 format('Initial motion vector (m/s): ',f7.2,1x,f7.2, + /,'Output time step=',f6.1, + /,'Integration time step=',f6.1, + /,'No. of time steps= ',i6) endif c c *** Adams-Bashforth time steps to desired time c do k=1,ndts c ++ Track if(k .eq. 1) then call uvcal(flons(0),flats(0),tf(0),iftype, + cxp(0),cyp(0),cxc(0),cyc(0), + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt, + glon,glat,gp,gt,ugrid,vgrid,vwts, + wxp(0),wyp(0),wxc(0),wyc(0), + cx(0),cy(0)) c fxm1 = cx(0)/(dtr*eradm*cos(dtr*flats(0))) fym1 = cy(0)/(dtr*eradm) flons(k) = flons(k-1) + dtss*fxm1 flats(k) = flats(k-1) + dtss*fym1 call uvcal(flons(1),flats(1),tf(1),iftype, + cxp(1),cyp(1),cxc(1),cyc(1), + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt, + glon,glat,gp,gt,ugrid,vgrid,vwts, + wxp(1),wyp(1),wxc(1),wyc(1), + cx(1),cy(1)) c fxm0 = cx(1)/(dtr*eradm*cos(dtr*flats(1))) fym0 = cy(1)/(dtr*eradm) else x1 = -0.5 x0 = 1.5 flons(k) = flons(k-1) + x1*dtss*fxm1 + x0*dtss*fxm0 flats(k) = flats(k-1) + x1*dtss*fym1 + x0*dtss*fym0 fxm1 = fxm0 fym1 = fym0 c c ++ Prepare for next time step c call uvcal(flons(k),flats(k),tf(k),iftype, + cxp(k),cyp(k),cxc(k),cyc(k), + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt, + glon,glat,gp,gt,ugrid,vgrid,vwts, + wxp(k),wyp(k),wxc(k),wyc(k), + cx(k),cy(k)) fxm0 = cx(k)/(dtr*eradm*cos(dtr*flats(k))) fym0 = cy(k)/(dtr*eradm) endif c Check to make sure track does not go outside of c model domain. If it does, stop the time integration. if (flats(k) .le. glat(1) .or. + flats(k) .ge. glat(nxy) .or. + flons(k) .le. glon(1) .or. + flons(k) .ge. glon(nxx) ) then exit endif enddo c c *** Extract forecast at output time steps kk = 0 do k=1,ndts if (mod(k,intsave) .eq. 0) then kk = kk + 1 flon(kk) = flons(k) flat(kk) = flats(k) fvmax(kk)= vmaxs(k) endif enddo c ierr = 0 c return end subroutine pninit c This routine initializes needed constants c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /pncon2/ aup,avp,efup,efvp,auc,avc,betax,betay,bcfac common /pncon3/ boxhx,boxhy,boxhr c c Numerical constants pi = 3.14159265 dtr = pi/180.0 c c Physical constants eradkm = 6371.0 eradm = eradkm*1000.0 c c Conversion factors deg2m = 2.0*pi*eradm/360.0 ckt2ms = 1.0/1.944 c c Missing data flag rmiss = -999. c c Parameters for weighting persistence track forecast c for the case where iftype=1. c c aup,avp are t=0 weights (non-dimensional). aup = 0.85 avp = 0.85 c c efup,efvp are exponential time decay factors (1/hr). c Final persistence weights for up,vp are wup = aup*exp(-t*efup) and c wvp = avp*exp(-1*efvp) c GFS motion vector weights are 1-wup, 1-wvp. efup = 1.0/6.0 efvp = 1.0/6.0 c c Scaling GFS motion vector. Normally keep these 1.0 unless c you do not want the persistence and GFS weightst to add up to 1. auc = 1.0 avc = 1.0 c c Specify half width (km) of the lat,lon rectangular area for c the horizontal average of the wind field c boxhx = 500.0 c boxhy = 500.0 boxhx = 750.0 boxhy = 750.0 c c Specify radius of circle (boxhr) in km for the horizontal average. c If boxhr > 0.0 this over-rides boxhx and boxhy, which uses a rectangular area. c boxhr = 350.0 c boxhr = 450.0 boxhr = 600.0 c c Beta drift speeds (m/s) and cosine factor (bcfac). The beta drift speeds c decrease with latitude if bcfac=1.0, or are constant if bcfac=0.0. betax = -0.325 betay = 0.250 bcfac = 1.0 c c betax =-1.0 c betay = 1.0 c return end subroutine vwcal(gp,mxp,nxp,mxpt,pbot,ptop,vwts,ierr) c The routine calculates the vertical weights for calculating the c the storm motion vector from the GFS analyses. Missing pressure c levels in the input pressure levels (gp) are accounted for, provided c gp .le. 0.0 for missing levels. The weights are for a mass-weighted c average from pressure levels pbot to ptop. c dimension gp(mxp),vwts(mxpt) c c Initialize all the weights to zero do m=1,nxp vwts(m) = 0.0 enddo ierr = 2 c c Find indices closest to ptop and pbot iptop = 0 do m=1,nxp if (gp(m) .gt. 0.0 .and. gp(m) .ge. ptop) then iptop = m exit endif enddo c if (iptop .eq. 0) return c ipbot = 0 do m=nxp,1,-1 if (gp(m) .gt. 0.0 .and. gp(m) .le. pbot) then ipbot = m exit endif enddo c if (ipbot .eq. 0) return c c nlev = 1 + ipbot-iptop if (nlev .lt. 1) then return elseif (nlev .eq. 1) then vwts(iptop) = 1.0 elseif (nlev .eq. 2) then vwts(iptop) = 0.5 vwts(ipbot) = 0.5 else delp = gp(ipbot)-gp(iptop) do m=iptop,ipbot if (gp(m) .gt. 0.0) then c Find nearest good level above this level up to ptop if (m .eq. iptop) then m1 = iptop else do mm=m-1,iptop,-1 if (gp(mm) .gt. 0.0) then m1 = mm exit !inner loop endif enddo endif c if (m .eq. ipbot) then m2 = ipbot else do mm=m+1,ipbot if (gp(mm) .gt. 0.0) then m2 = mm exit !inner loop endif enddo endif vwts(m) = 0.5*(gp(m2)-gp(m1))/delp endif c c write(6,601) gp(m),m1,m2,gp(m1),gp(m2) c 601 format('p,m1,m2,p1,p2: ',f8.1,2(1x,i4),2(1x,f8.1)) enddo endif c ierr = 0 c c write(6,603) iptop,ipbot,gp(iptop),gp(ipbot) c 603 format(/,'iptop,ipbot,ptop,pbot: ',2(1x,i3),2(1x,f8.1)) c do m=1,nxp c write(6,602) gp(m),vwts(m) c 602 format(' p=',f8.1,' vwt=',f8.4) c enddo c return end subroutine uvcal(slon,slat,t,iftype, + cxp,cyp,cxc,cyc, + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt, + glon,glat,gp,gt,ugrid,vgrid,vwts, + wxp,wyp,wxc,wyc, + cx,cy) c c This routine estimates the storm motion vector (cx,cy) c as a weighted average of GFS motion vector, persistence c and a beta drift. c c Input: c slon: storm longitude ( 0 to 360) c slat: storm latitude (-90 to 90) c vmax: max wind (kt) c t : forecast time (hr) c iftype: Method for averaging persistence and climatology c (=0 for persistence, =1 for model, =2 for time weighted average) c cxp,cyp: t=0 motion vector (m/s) c c mxx - Maximum number of longitudes c mxy - Maximum number of latitudes c mxp - Maximum number of pressure levels c mxt - Maximum number of forecast times c c nxx, nxy, nxp, nxt - Actual number of longitudes, latitudes, etc c glon - Longitudes of the 3-D wind field (0 to 360 deg) c glat - Latitudes of the 3-D wind field (-90 to 90 deg) c gp - Pressure levels of the 3-D wind field (hPa) ordered from low c to high pressure c gt - Times (hr, 0 to nxt) of the 3-D wind field c ugrid - u (m/s) as a function of lon,lat,p,t c vgrid - v (m/s) as a function of lon,lat,p,t c vwts - Vertical pressure weights c c Output: c wxp: The weight on the x-component of persistence c wyp: The weight on the y-component of persistence c wxc: The weight on the x-component of the GFS model c wyc: The weight on the y-component of the GFS model c cx: The x-component of the motion vector (m/s) c cy: The y-component of the motion vector (m/s) c dimension glon(mxx),glat(mxy),gp(mxp),gt(0:mxt) dimension ugrid(mxx,mxy,mxp,0:mxt),vgrid(mxx,mxy,mxp,0:mxt) dimension vwts(mxp) c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /pncon2/ aup,avp,efup,efvp,auc,avc,betax,betay,bcfac c c Get GFS motion vector call gfsuv(slon,slat,t,cxc,cyc, + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt,vwts, + glon,glat,gp,gt,ugrid,vgrid) c Calcuate the GFS and persistence weights if (iftype .eq. 0) then wxp = 1.0 wyp = 1.0 wxc = 0.0 wyc = 0.0 elseif (iftype .eq. 1) then wxp = 0.0 wyp = 0.0 wxc = 1.0 wyc = 1.0 else wxp = aup*exp(-t*efup) wyp = avp*exp(-t*efvp) wxc = auc*(1.0-wxp) wyc = avc*(1.0-wyp) endif c cx = wxc*cxc + wxp*cxp cy = wyc*cyc + wyp*cyp c c Add beta drift cslat = cos(dtr*slat*bcfac) sslat = 1.0 if (slat .lt. 0.0) sslat=-1.0 c cx = cx + cslat*betax cy = cy + cslat*betay*sslat c return end subroutine gfsuv(slon,slat,t,u,v, + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt,vwts, + glon,glat,gp,gt,ugrid,vgrid) c c This routine calculate the values of u,v (m/s) from c vertically and horizontally averaged horizontal winds, centered c on the position (slon,slat). c dimension glon(mxx),glat(mxy),gp(mxp),gt(0:mxt) dimension ugrid(mxx,mxy,mxp,0:mxt),vgrid(mxx,mxy,mxp,0:mxt) dimension vwts(mxp) c common /pncon1/ pi,dtr,ckt2ms,deg2m,eradkm,eradm,rmiss common /pncon3/ boxhx,boxhy,boxhr c c Local variables parameter (mxpp=100) dimension ubarp(mxpp),vbarp(mxpp) c c Calculate lon,lat,t increments dlon = glon(2)-glon(1) dlat = glat(2)-glat(1) dt = gt(1)-gt(0) c c Find the time index of the time level just before t k0 = ifix( (t-gt(0))/dt ) if (k0 .eq. nxt) k0 = nxt-1 c c Find the time weights k1 = k0+1 tnorm = (t-gt(k0))/dt wtk0 = 1.0-tnorm wtk1 = tnorm c c Find the indices of the point to the lower left of the point (slon,slat) i0 = 1 + nint( (slon-glon(1))/dlon ) j0 = 1 + nint( (slat-glat(1))/dlat ) c c Do not let index be outside the domain if (i0 .gt. nxx) i0 = nxx if (i0 .lt. 1) i0 = 1 if (j0 .gt. nxy) j0 = nxy if (j0 .lt. 1) j0 = 1 c c Horizontal average of u,v cfactem = cos(dtr*slat) c if (boxhr .gt. 0.0) then jinc = 1 + nint(boxhr/(111.1*dlat)) iinc = 1 + nint(boxhr/(111.1*dlon*cfactem)) else jinc = 1 + nint(boxhy/(111.1*dlat)) iinc = 1 + nint(boxhx/(111.1*dlon*cfactem)) endif c ibgn = i0 - iinc iend = i0 + iinc jbgn = j0 - jinc jend = j0 + jinc c if (ibgn .lt. 1) ibgn = 1 if (jbgn .lt. 1) jbgn = 1 if (iend .gt. nxx) iend = nxx if (jend .gt. nxy) jend = nxy c c u = ugrid(i0,j0,8,k0) c v = vgrid(i0,j0,8,k0) c c do m=1,nxp if (gp(m) .gt. 0.0) then ubarp(m) = 0.0 vbarp(m) = 0.0 acount = 0.0 c if (boxhr .gt. 0.0) then do j=jbgn,jend do i=ibgn,iend dy = 111.1*(slat-glat(j)) dx = 111.1*(slon-glon(i))*cfactem rtem = sqrt(dx*dx+dy*dy) if (rtem .le. boxhr) then acount = acount + 1.0 ubarp(m) = ubarp(m) + wtk0*ugrid(i,j,m,k0) + + wtk1*ugrid(i,j,m,k1) vbarp(m) = vbarp(m) + wtk0*vgrid(i,j,m,k0) + + wtk1*vgrid(i,j,m,k1) endif enddo enddo else do j=jbgn,jend do i=ibgn,iend acount = acount + 1.0 ubarp(m) = ubarp(m) + wtk0*ugrid(i,j,m,k0) + + wtk1*ugrid(i,j,m,k1) vbarp(m) = vbarp(m) + wtk0*vgrid(i,j,m,k0) + + wtk1*vgrid(i,j,m,k1) enddo enddo endif c if (acount .ge. 1.0) then ubarp(m) = ubarp(m)/acount vbarp(m) = vbarp(m)/acount else ubarp(m) = 0.0 vbarp(m) = 0.0 endif else ubarp(m) = 0.0 vbarp(m) = 0.0 endif enddo c c Calculate vertical average u = 0.0 v = 0.0 do m=1,nxp u = u + vwts(m)*ubarp(m) v = v + vwts(m)*vbarp(m) enddo c c write(6,500) t,i0,j0,k0,u,v,jinc,iinc,tnorm,wtk0,wtk1 c 500 format('t,i0,j0,k0,u,v,jinc,iinc,tnorm,wtk0,wtk1: ', c + f6.1,3(1x,i4),2(1x,f6.1),2(i4),1x,3(1xf6.3)) c return end subroutine getgfsl(iyr,imon,iday,itime,idthr,ndt, + mxx,mxy,mxp,mxt,nxx,nxy,nxp,nxt, + glon,glat,gp,gt,ugrid,vgrid,tgrid,zgrid,rgrid) c c This routine gets the GFS files and puts them into lon,lat,p,t arrays. c The data is assumed to be in packed ASCII files that have the same c lat,lon,pressure grids. c c This version looks for fields from a 6 hr old forecast if c the file is missing. This option is only valid if forecast c fields are used. If analysis fields are used, that option is c skipped. c c 06/15/18 (SNS) - added capability to get 12 hr old forecast to c allow for ECMWF version of tab model to run c c Passed variables dimension ugrid(mxx,mxy,mxp,0:mxt),vgrid(mxx,mxy,mxp,0:mxt) dimension tgrid(mxx,mxy,mxp,0:mxt),zgrid(mxx,mxy,mxp,0:mxt) dimension rgrid(mxx,mxy,mxp,0:mxt) c dimension glon(mxx),glat(mxy),gp(mxp),gt(0:mxt) c c Local variables parameter (mxtt=100,mxpt=100) character *1 cxy character *21 fnpack(0:mxtt) ,fntemp character *21 fnpackl(0:mxtt),fntempl character *21 fnpackl12(0:mxtt),fntempl12 dimension ipfound(mxpt) c parameter (imax=100000) dimension tra(imax) character *2 code(imax) character *1 type c character *50 dpath character *80 ffntemp c c Set ipath=1 if packed ASCII files are in a different directory than c where the gfsget is being run. If so, also specify the full path to c the files. If the packed ASCII files are in the same directory, then c set ipath=0 ipath = 0 dpath='/mnt/hfipnas1/gfspack/' c c Set ktmin to the minimum time (hr) the forecast files c must exist for. If the forecast does not extend to that time, c try to get the fields from the 6 hour old run. ktmin = 120 ktmini = ktmin/idthr c c Set the pressure levels to look for c nxp = 11 nxp = 7 c gp( 1) = 100.0 c gp( 2) = 150.0 gp( 1) = 200.0 gp( 2) = 250.0 gp( 3) = 300.0 gp( 4) = 400.0 gp( 5) = 500.0 gp( 6) = 700.0 gp( 7) = 850.0 c gp(10) = 1000.0 c gp(10) = 925.0 c gp(11) = 1000.0 c do m=1,nxp ipfound(m) = 0 enddo c c Set ifcst=1 to use forecast files. Otherwise analysis fields will be used. ifcst=1 c c Unit number for reading the pack files lupack = 24 c c Create t=0 packed file name fntemp='A00011_X0901_PACK.DAT' c iyr2 = iyr - 100*(iyr/100) idayt = iday cxy = 'X' if (itime .eq. 12 .or. itime .eq. 18) idayt = idayt + 50 if (itime .eq. 6 .or. itime .eq. 18) cxy='Y' write(fntemp(5:6),300) iyr2 300 format(i2.2) fntemp(8:8) = cxy write(fntemp(9:12),305) imon,idayt 305 format(i2.2,i2.2) fnpack(0) = fntemp c c Create t=0 packed file name from 6 hr old forecast fntempl='A00611_X0901_PACK.DAT' c itlag = -6 call tadd(iyr ,imon ,iday ,itime ,itlag, + iyrl,imonl,idayl,itimel) c iyrl2 = iyrl - 100*(iyrl/100) idayt = idayl cxy = 'X' if (itimel .eq. 12 .or. itimel .eq. 18) idayt = idayt + 50 if (itimel .eq. 6 .or. itimel .eq. 18) cxy='Y' c write(fntempl(5:6),300) iyrl2 fntempl(8:8) = cxy write(fntempl(9:12),305) imonl,idayt fnpackl(0) = fntempl c c Create t=0 packed file name from 12 hr old forecast fntempl12='A00611_X0901_PACK.DAT' c itlag12 = -12 call tadd(iyr ,imon ,iday ,itime ,itlag12, + iyrl12,imonl12,idayl12,itimel12) c iyrl212 = iyrl12 - 100*(iyrl12/100) idayt = idayl12 cxy = 'X' if (itimel12 .eq. 12 .or. itimel12 .eq. 18) idayt = idayt + 50 if (itimel12 .eq. 6 .or. itimel12 .eq. 18) cxy='Y' c write(fntempl12(5:6),300) iyrl212 fntempl12(8:8) = cxy write(fntempl12(9:12),305) imonl12,idayt fnpackl12(0) = fntempl12 c c Calculate the rest of the packed file names do k=1,ndt itadd = k*idthr c if (ifcst .eq. 1) then write(fntemp(2:4),310) itadd 310 format(i3.3) fnpack(k) = fntemp c write(fntempl(2:4),310) itadd - itlag fnpackl(k) = fntempl c write(fntempl12(2:4),310) itadd - itlag12 fnpackl12(k) = fntempl12 else call tadd(iyr ,imon ,iday ,itime ,itadd, + iyra,imona,idaya,itimea) c iyr2 = iyra - 100*(iyra/100) idayt = idaya cxy = 'X' if (itimea .eq. 12 .or. itimea .eq. 18) idayt = idayt + 50 if (itimea .eq. 6 .or. itimea .eq. 18) cxy='Y' write(fntemp(5:6),300) iyr2 fntemp(8:8) = cxy write(fntemp(9:12),305) imon,idayt fnpack(k) = fntemp endif enddo c nxt = -1 c Process each file do k=0,ndt gt(k) = float(k)*float(idthr) fntemp = fnpack(k) if (ipath .eq. 0) then ffntemp = fntemp else ffntemp=trim(dpath) // trim(fntemp) endif c if (ifcst .eq. 1) then open(file=ffntemp,unit=lupack,form='formatted', + status='old',iostat=istat) if (istat .ne. 0) then c If you get to here, there was a problem opening the model field. c Try field from 6 hr old forecast instead. c close(lupack) fntemp = fnpackl(k) if (ipath .eq. 0) then ffntemp = fntemp else ffntemp=trim(dpath) // trim(fntemp) endif open(file=ffntemp,unit=lupack,form='formatted', + status='old',iostat=istat) if (istat .ne. 0) then c If you get to here, there was problem opening the 6 hr old c model field. Try field from 12 hr old forecast instead. c close(lupack) fntemp = fnpackl12(k) if (ipath .eq. 0) then ffntemp = fntemp else ffntemp=trim(dpath) // trim(fntemp) endif open(file=ffntemp,unit=lupack,form='formatted', + status='old',err=9000) endif endif endif c c write(6,700) k,ffntemp c 700 format('k= ',i3,' packfile= ',a80) c c Read main header on packed data file read(lupack,100) wx,dayx,utcx,rlatmn,rlatmx,rlonmn,rlonmx, + dlat,dlon 100 format(1x,f3.0,f7.0,f5.0,4f8.3,2f4.1) c C Calculate lat (-90 to 90) and lon (0 to 360) points epsil=0.001 nlati = ifix((rlatmx-rlatmn)/dlat + epsil) nloni = ifix((rlonmx-rlonmn)/dlon + epsil) ipts = (nlati+1)*(nloni+1) nlat1 = nlati + 1 nlon1 = nloni + 1 nrow = 1 + (ipts-1)/36 c nxx = nlon1 nxy = nlat1 c if (fntemp(1:1) .eq. 'A' .or. fntemp(1:1) .eq. 'R') then rlonleft = 180.0 - (rlonmx-180.0) else rlonlef = rlonmn endif c do i=1,nxx glon(i) = rlonleft + dlon*float(i-1) if (glon(i) .lt. 0.0) glon(i) = glon(i) + 360.0 if (fntemp(1:1) .eq. 'A' .or. fntemp(1:1) .eq. 'R') then c Let right boundary go past 360. for non-global domains else if (glon(i) .gt. 360.0) glon(i) = glon(i) - 360.0 endif enddo c do j=1,nxy glat(j) = rlatmn + dlat*float(j-1) enddo c c Read in each block of packed data packloop: do read(lupack,105,err=9000,end=600) type,ptem,bsub,smpy 105 format(1x,a1,1x,f6.1,2(1x,g15.9)) c do n=1,nrow is = 1 + (n-1)*36 ie = is + 35 read(lupack,110,end=9000,err=9000) (code(i),i=is,ie) 110 format(36(a2)) enddo c c Determine the pressure level index ipindex = 0 do m=1,nxp delp = abs(gp(m)-ptem) if (delp .lt. epsil) then ipindex = m ipfound(m) = 1 endif enddo if (ipindex .le. 0) cycle c c Unpack the data and put it in the proper array do i=1,ipts ix = idecod(code(i)) tra(i) = ix * smpy - bsub enddo c if (type .eq. 'U') then do i=0,nlati do j=0,nloni ij = (j+1) + i*(nloni+1) ugrid(j+1,i+1,ipindex,k) = tra(ij) enddo enddo elseif (type .eq. 'V') then do i=0,nlati do j=0,nloni ij = (j+1) + i*(nloni+1) vgrid(j+1,i+1,ipindex,k) = tra(ij) enddo enddo elseif (type .eq. 'Z') then do i=0,nlati do j=0,nloni ij = (j+1) + i*(nloni+1) zgrid(j+1,i+1,ipindex,k) = tra(ij) enddo enddo elseif (type .eq. 'T') then do i=0,nlati do j=0,nloni ij = (j+1) + i*(nloni+1) tgrid(j+1,i+1,ipindex,k) = tra(ij) enddo enddo elseif (type .eq. 'R') then do i=0,nlati do j=0,nloni ij = (j+1) + i*(nloni+1) rgrid(j+1,i+1,ipindex,k) = tra(ij) enddo enddo endif c enddo packloop 600 continue c close(lupack) nxt = k enddo c c Normal return c rmiss = -9999. do m=1,nxp if (ipfound(m) .eq. 0) gp(m) = rmiss enddo c return c c Return from here if some requested files are missing 9000 continue c c Set pressure value to rmiss if data was not found at that level rmiss = -9999. do m=1,nxp if (ipfound(m) .eq. 0) gp(m) = rmiss enddo c return c end subroutine tab_err_handling(ierrtype,luerr,ierr) c This routine handles error processing for tcliper. c Types of errors handled (ierrtype=): c -1: input file failed to open c -2: error reading input file c -3: TAB not available in basin c -4: error in tab() c Input c ierrtype - error type to process c luerr - unit number to write to (assumes file is already open) c ierr - error code to report c IMPLICIT NONE c c ++ Passed variables integer, intent(in) :: ierrtype, luerr, ierr c if (ierrtype .eq. 1) then write(luerr,950) 950 format(/,'Error opening tab.com input file') stop elseif (ierrtype .eq. 2) then write(luerr,970) ierr 970 format(/,'Error reading tab.com input file, istat=',i4) stop elseif (ierrtype .eq. 3) then write(luerr,960) 960 format(/,'TAB not available for this basin') stop elseif (ierrtype .eq. 4) then write(luerr,980) ierr 980 format(/,'Error in tab routine, ierr=',i2) stop else write(luerr,999) 999 format(/,'Unrecognized error in TAB') stop endif c return end subroutine tab_err_handling