program lsdiag c This program performs diagnostic calculations using c large-scale analyses in packed form. c c Version 16.6.1 c c Modified Apr 2005 - 6 hour data increment c Modified 20 Apr 2006 - for 2006 season c Modified 26 Aug 2006 - Legacy svn routines removed c 13 Sep 2006 - Experimental shear and related calculation c added c 16 Sep 2006 - New storm-relative D200 calculation c 2008 - lcmod calculations added c 25 Mar 2009 - U200C,V200C output added c 02 Apr 2010 - RH bias correction for reanalsysis fields c added c 07 Apr 2011 - New predictors added (V000,V850,V500,V300, c TGRD,TADV) c 07 Jun 2012 - Modifications for global pack files c Legacy tavg and ascal routines removed c 05 Feb 2013 - Option for reducing computational grid size c added (set by igcrop) c 01 Apr 2013 - Option to calculate just the t=0 predictors c (set by izero) c 10 Oct 2013 - sgnlon modified if E/W explictly indicated c on LON line of input file c 05 Nov 2013 - New sv7cal routine for predictors for 7 day c models c 12 Jan 2014 - New sv7cal2 routine with 925 hPa added and c option to print entire sounding c 30 Apr 2015 - New variables for balanced temperature field c in T averaging area in sv6cal (requires new c faavg routine and distk utility routine) c Legacy ssfcal routine removed. c New storm relative helicity predictors added c 26 Dec 2015 - Joe Cione downdraft variable added (CFLX) c Pressure vertical velocity added c (O500, O700) c 13 Jan 2016 - Error flag added to septoz routine in lcmod c so lsdiag.x processing continues when c sounding cannot be interpolated. c 15 Mar 2016 - TPW variables added c 05 Nov 2016 - Utility routines removed (jday, jdayi, tadd) c 10 Jan 2017 - ir00param3() replaced with ir00param4() c for v16.3.0 (KM) c 15 Jan 2018 - Model TPW predictors calculated and written c to lsdiag.dat file for all forecast times c instead of just t=0 hr. New variables are c labeled PW01, PW02, ... for 16.4.0 (MD) c 18 Feb 2018 - pgoes dimension bug fix put back in for c 16.4.1 (MD) c 18 Mar 2020 - Legacy sv6cal and sv7cal2 removed, c ilaghr increased to 12, added ntpw (SS) c 02 Feb 2023 - NCO specifications addressed (KM) c 06 Feb 2023 - (1) O500, O700, HE05, HE07 overflow fixed c (MD) (v16.6.0) c (2) Options to replace weekly Reynolds SST c with c (1st choice) Daily Reynolds SST c (2nd choice) NCODA SST c (3rd choice) Climo SST c because weekly Reynolds is no longer c available after Nov 2022. c 15 Feb 2023 (1) Bug fix for selecting SST option (MD) c (v16.6.1) c (2) writing the SST choice to log file. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) c C.DUP lsdiag,prima common /primv/ u,v,z,t,r dimension u(ixmx,iymx,ipmx),v(ixmx,iymx,ipmx),z(ixmx,iymx,ipmx) dimension t(ixmx,iymx,ipmx),r(ixmx,iymx,ipmx) c C.DUP lsdiag,index common /indexa/ indexp,indexv common /indexn/ indexvn parameter (nvar=5) dimension indexp(ipmx),indexv(nvar) character *1 indexvn(nvar) c Specify allowable variable names on data input file data indexvn /'U','V','Z','T','R'/ c C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 c C.DUP lsdiag,plevs common /plevs/ plev,plevr dimension plev(ipmx),plevr(ipmx) C.DUP lsdiag,plevsd DATA PLEV / 1070.,1000.,950.,925.,900.,850.,700.,600.,500., + 400., 350.,300.,250.,200.,150.,100.,70.,50., + 25., 10., 5./ DATA PLEVR / 1013.,1000.,950.,925.,900.,850.,700.,600.,500., + 400., 350.,300.,250.,200.,150.,100.,70.,50., + 25., 10., 5./ c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c C.DUP lsdiag,rawdfn common /rawdf/ rawdfn character *40 rawdfn c C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd C.DUP lsdiag,unnumd data lulog,luinp,luout,lurawd /21,22,23,24/ c c Local variables character *20 logname,inname,outname c c ******************** c c Specify log file name logname = 'lsdiag.log' c c Specify input file name inname = 'lsdiag.inp' c c Specify data output file name outname = 'lsdiag.dat' c c Open the log file open(unit=lulog,file=logname,status='replace',form='formatted') write(lulog,300) 300 format(' LSDIAG BEGINS, LOG FILE OPENED') c c Open the data input and output files open(unit=luinp,file=inname,status='old',form='formatted') open(unit=luout,file=outname,status='replace',form='formatted') c c Calculate constants call pncal c c Calcualte SHIPS variables (7-day input files with 6 hr data) call sv7cal c stop END program lsdiag c subroutine aunpack(ierr,sgnlon) c This routine unpacks the data file and puts the variables c in three-dimensional arrays. If there is a problem opening c or reading the data file, ierr=1, otherwise ierr=0. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,prima common /primv/ u,v,z,t,r dimension u(ixmx,iymx,ipmx),v(ixmx,iymx,ipmx),z(ixmx,iymx,ipmx) dimension t(ixmx,iymx,ipmx),r(ixmx,iymx,ipmx) c C.DUP lsdiag,index common /indexa/ indexp,indexv common /indexn/ indexvn parameter (nvar=5) dimension indexp(ipmx),indexv(nvar) character *1 indexvn(nvar) C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,rawdfn common /rawdf/ rawdfn character *40 rawdfn C.DUP lsdiag,plevs common /plevs/ plev,plevr dimension plev(ipmx),plevr(ipmx) C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c c Local variables parameter(imax=100000) c dimension tra(imax) character *2 code(imax) character *1 type character *10 vlabel c c ******************** c c Open data input file open(unit=lurawd,file=rawdfn,status='old',form='formatted', + err=900) c write(lulog,300) rawdfn 300 format(/,' Begin unpacking rawd file: ',a40) c c Initialize all index variables to zero do i=1,ipmx indexp(i) = 0 enddo c do i=1,nvar indexv(i) = 0 enddo c C READ MAIN HEADER ON PACKED DATA FILE READ(lurawd,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 FOR WEST PACIFIC RAWD FILES, SWITCH RLONMN AND RLONMX C TO ACCOUNT FOR DEG EAST LONGITUDE c IF (rawdfn(1:2) .EQ. 'WP') THEN c RLMIN = RLONMN c RLMAX = RLONMX c RLONMN = -RLMAX c RLONMX = -RLMIN c ENDIF C C WRITE HEADER INFORMATION WRITE(LULOG,200) DAYX,UTCX 200 FORMAT(/,' DATA DATE AND TIME: ',F7.0,1X,F5.0) C WRITE(LULOG,205) RLATMN,RLATMX,RLONMN,RLONMX 205 FORMAT(/,' MIN,MAX LAT: ',F6.1,1X,F6.1, + ' MIN,MAX LON: ',F6.1,1X,F6.1) C WRITE(LULOG,210) DLAT,DLON 210 FORMAT(/,' DLAT,DLON= ',F4.1,1X,F4.1) C C CALCULATE LAT AND LON 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 C C UNPACK THE DATA NROW = 1 + (IPTS-1)/36 c packloop: do c loop exits when read error, end-of-file, or lev/var error READ(lurawd,105,ERR=900,END=600) TYPE,PTEM,BSUB,SMPY 105 FORMAT(1X,A1,1X,F6.1,2(1X,G15.9)) c C DETERMINE THE PRESSURE LEVEL DO I=1,ipmx IF (PTEM .EQ. PLEV(I)) EXIT IF (I .EQ. ipmx) THEN WRITE(lulog,960) PTEM 960 FORMAT(2X,'PRESSURE= ',E11.4,' NOT FOUND IN PLEV ARRAY') STOP ENDIF ENDDO C indexp(i) = 1 lnow = i c C READ IN PACKED DATA DO N=1,NROW IS = 1 + (N-1)*36 IE = IS+35 READ(lurawd,110,END=902,ERR=902) (CODE(I),I=IS,IE) 110 FORMAT(36(A2)) ENDDO 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 indexv(1) = 1 DO I = 0,NLATI DO J = 0,NLONI IJ = (J+1) + I*(NLONI+1) u(J+1,I+1,LNOW) = TRA(IJ) ENDDO ENDDO elseif (type .eq. 'V') then indexv(2) = 1 DO I = 0,NLATI DO J = 0,NLONI IJ = (J+1) + I*(NLONI+1) v(J+1,I+1,LNOW) = TRA(IJ) ENDDO ENDDO elseif (type .eq. 'Z') then indexv(3) = 1 DO I = 0,NLATI DO J = 0,NLONI IJ = (J+1) + I*(NLONI+1) z(J+1,I+1,LNOW) = TRA(IJ) ENDDO ENDDO elseif (type .eq. 'T') then indexv(4) = 1 DO I = 0,NLATI DO J = 0,NLONI IJ = (J+1) + I*(NLONI+1) t(J+1,I+1,LNOW) = TRA(IJ) ENDDO ENDDO elseif (type .eq. 'R') then indexv(5) = 1 DO I = 0,NLATI DO J = 0,NLONI IJ = (J+1) + I*(NLONI+1) r(J+1,I+1,LNOW) = TRA(IJ) ENDDO ENDDO else write(lulog,970) type 970 format(/,' Unrecognized variable type: ',a1) c stop endif enddo packloop 600 continue c c Calculate latitude (deg) and longitudes (deg W neg.) c the of grid points do j=1,nlat1 rlatd(j) = rlatmn + dlat*float(j-1) enddo c do i=1,nlon1 rlond(i) = -(rlonmx - dlon*float(i-1)) enddo c if (sgnlon .lt. 0.0 .and. rlond(1) .ge. 0.0) then c Recalcuate longitudes for neg lon case rlonleft = rlonmx-360.0 do i=1,nlon1 rlond(i) = rlonleft + dlon*float(i-1) enddo endif c c List pressure levels and variable types found on rawd file write(lulog,305) 305 format(/, +' Pressure levels (mb) and variables found on rawd file:') do i=1,ipmx if (indexp(i) .eq. 1) then write(lulog,310) plevr(i) 310 format(1x,f6.1) endif enddo c vlabel = ' ' ii = 0 do i=1,nvar if (indexv(i) .eq. 1) then ii = ii+1 ii2 = 2*ii vlabel(ii2:ii2) = indexvn(i) endif enddo write(lulog,315) vlabel 315 format(1x,a10) c c Close the rawd file and set the error flag close(lurawd) ierr=0 return c 900 write(lulog,950) rawdfn 950 format(/,' Error opening rawd file: ',a40) close(lurawd) ierr=1 return c 902 write(lulog,952) rawdfn 952 format(/,' Error reading rawd file: ',a40) close(lurawd) ierr=1 return c end subroutine aunpack c c integer function idecod (code) c c parameter (base = 32) c character*(*) code c character dgtb*(base-1) c parameter (dgtb = '123456789ABCDEFGHIJKLMNOPQRSTUV') c c idecod = index (dgtb,code(1:1)) * base c 1 + index (dgtb,code(2:2)) c c return c END function idecod c integer function idecod (code) c hp version c character*(*) code character dgtb*(31) parameter (dgtb = '123456789ABCDEFGHIJKLMNOPQRSTUV') c idecod = index (dgtb,code(1:1)) * 32 1 + index (dgtb,code(2:2)) c return END function idecod c subroutine sv7cal c This routine calculates the variables needed for c the SHIPS intensity forecast model c for each case in the input file. c c This version uses the input file with 7-day forecasts. c c The lsdiag input file includes data at 6 hour intervals. c parameter (itmax=28) c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,rawdfn common /rawdf/ rawdfn character *40 rawdfn C.DUP lsdiag,plevh parameter (nplevh=10) common /plevhs/ up,vp,zp,tp,rp,tbarp,plevh dimension up(nplevh),vp(nplevh),zp(nplevh) dimension tp(nplevh),rp(nplevh) dimension tbarp(nplevh) dimension plevh(nplevh) C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Local variables dimension islat(-2:itmax),islon(-2:itmax) dimension slat(-2:itmax), slon(-2:itmax) character *1 loncon c dimension iu200(0:itmax),it200(0:itmax) dimension iu200c(0:itmax),iv200c(0:itmax) dimension it150(0:itmax),it250(0:itmax) dimension ishrd(0:itmax),ishrs(0:itmax) dimension ishtd(0:itmax),ishts(0:itmax) dimension ishrg(0:itmax) dimension ishxu(0:itmax),ishxl(0:itmax) dimension ivorb(0:itmax),idivb(0:itmax),irefc(0:itmax) dimension ipefc(0:itmax) dimension irhlo(0:itmax),irhhi(0:itmax),irhmd(0:itmax) dimension ipslv(0:itmax) dimension iz000(0:itmax),it000(0:itmax),ir000(0:itmax) c dimension ie000(0:itmax) dimension iepos(0:itmax),ieneg(0:itmax) dimension iepss(0:itmax),ienss(0:itmax) c dimension fxy(ixmx,iymx) c character *2 cbasin,csnum character *4 sname,tlabel character *170 hline,tline,iline(60) c c Variables for vertical profiles parameter (nvp=10) c dimension pvp(nvp),uvp(nvp),vvp(nvp),tvp(nvp),rvp(nvp) dimension qvp(nvp) dimension wvp(nvp),wt1(nvp),wt2(nvp),tevp(nvp),tevps(nvp) dimension zvp(nvp) c c Variables for J. Cione flux variable dimension icflx(0:itmax) c c Variables for omega calculation dimension divbp(nvp),omegap(nvp) dimension io500(0:itmax),io700(0:itmax) c c Vertical profile arrays for Kerry Emanuel MPI code and lcmod dimension pkmpi(nvp+1),tkmpi(nvp+1),qkmpi(nvp+1),rkmpi(nvp+1) dimension ivkmpi(0:itmax),ipkmpi(0:itmax) c SST arrays dimension irsst(-2:itmax),rsst(-2:itmax) dimension idsst(-2:itmax) dimension insst(-2:itmax) dimension icsst(-2:itmax) character *4 ssttype c c Variables for TPW calculation dimension psxy(ixmx,iymx),rsxy(ixmx,iymx),tpwxy(ixmx,iymx) dimension rptpw(nvp),tptpw(nvp),pptpw(nvp) parameter (mxtpw=itmax) c set maximum number of PW## outputs in lsdiag parameter (ntpw=20) c dimension itpw(0:mxtpw),itpw0(0:mxtpw) dimension itpwa(0:mxtpw,0:itmax) c c Variables for tracker and vortex removal code dimension uxyp(ixmx,iymx,nvp),vxyp(ixmx,iymx,nvp) dimension txyp(ixmx,iymx,nvp),rxyp(ixmx,iymx,nvp) dimension zxyp(ixmx,iymx,nvp) c dimension itlat(0:itmax),itlon(0:itmax) dimension ishdc(0:itmax),isddc(0:itmax),ishgc(0:itmax) dimension it15c(0:itmax),it20c(0:itmax),it25c(0:itmax) dimension iepoc(0:itmax),idivc(0:itmax) dimension itgrd(0:itmax),itadv(0:itmax) c dimension irtcc(0:itmax),irtxc(0:itmax) dimension itwac(0:itmax),itwxc(0:itmax),itwbc(0:itmax) dimension ipsfc(0:itmax),ipsfcc(0:itmax) dimension iv000(0:itmax),iv850(0:itmax) dimension iv500(0:itmax),iv300(0:itmax) c dimension ivvavg(0:itmax),ivmflux(0:itmax) dimension ivvavc(0:itmax) c parameter (mrcf=100,mtcf=16) dimension rcf(0:mrcf),thetacf(0:mtcf) dimension stheta(0:mtcf),ctheta(0:mtcf) dimension uraa(0:mrcf),vtaa(0:mrcf) dimension zaa(0:mrcf), taa(0:mrcf) dimension uraap(0:mrcf,nvp),vtaap(0:mrcf,nvp) dimension taap(0:mrcf,nvp), zaap(0:mrcf,nvp) dimension gtaap(0:mrcf,nvp) dimension ig250(0:itmax),ig200(0:itmax),ig150(0:itmax) c dimension zcrt(0:mrcf,0:mtcf),tcrt(0:mrcf,0:mtcf) dimension zraa(0:mrcf),traa(0:mrcf) c dimension isrh07(0:itmax),isrh05(0:itmax) c c Variables for proxy GOES data dimension pgoes(19) dimension ipgoes(0:itmax) c ******************** c c Set iatype=0 to use analyses to calculate predictors, c or iatype=1 to use forecasts for predictors iatype=1 c c If iatype=0, set last year (yyyy) to use re-analysis files c iray = 1999 iray = 1099 c c If iatype=0, set ibcrh=1 to apply the bias correction to the c RH from the reanalysis fields ibcrh=0 c c Set izero = 1 to only calculate t=0 variables, and fill in other c times with zeros izero = 0 c c Set igcrop=1 to crop the model lat/lon grid around the storm c center to increase efficiency. Set the distance (km) for cropping c the fields. igcrop=0 bufkm = 2500.00 c c Specify time interval between data points idelt=6 delt = float(idelt) c c Specify radii (km) for vertical profile variables c (These define averaging area for shear and other predictors) radinn = 200.0 radout = 800.0 c c Specify radius (km) for averaging tangential wind vortex tracking radtwa=600.0 c c Specify inner and outer radii (km) for azimuthal average c after filtering and vortex removal radinnc = 0.0 radoutc = 500.0 c c Specify radius (km) for area average relative vorticity radvor = 1000.0 c c Specify radius (km) for area average divergence raddiv = 1000.0 c c Specify radius (km) for radially averaging eddy momentum fluxes radrefc = 600.0 c c Specify radius (km) for radially averaging planetary momentum c fluxes radpefc = 600.0 c c Specify minimum points for area average calculations minpts=5 c c Specify vertical pressure levels (mb) for analysis pvp( 1) = 100.0 pvp( 2) = 150.0 pvp( 3) = 200.0 pvp( 4) = 250.0 pvp( 5) = 300.0 pvp( 6) = 400.0 pvp( 7) = 500.0 pvp( 8) = 700.0 pvp( 9) = 850.0 pvp(10) = 1000.0 c c Specify indices of pressures in pvp array i100 = 1 i150 = 2 i200 = 3 i250 = 4 i300 = 5 i400 = 6 i500 = 7 i700 = 8 i850 = 9 i000 = 10 c c Set imxshr=1 to include max shear variables imxshr=0 c c Missing values rmiss = 9999. imiss = 9999 c c Calculate deep layer mean weights for each pressure level dpvp = pvp(nvp)-pvp(1) c wvp( 1) = 0.5*(pvp( 2)-pvp( 1))/dpvp wvp(nvp) = 0.5*(pvp(nvp)-pvp(nvp-1))/dpvp c do n=2,nvp-1 wvp(n) = ( 0.5*(pvp(n+1)-pvp(n-1)) )/dpvp enddo c c Variables for center finding and vortex removal arrays nrcf = 24 drcf = 50.0 do k=0,nrcf rcf(k) = drcf*float(k) enddo c ntcf = 7 dthetacf = 360.0/float(ntcf+1) do m=0,ntcf thetacf(m) = dthetacf*float(m) ctheta(m) = cos(dtr*thetacf(m)) stheta(m) = sin(dtr*thetacf(m)) enddo c c Start the case by case loop nrawd=0 nskip=0 c caseloop: do c caseloop is the main loop for entire subroutine c c Read header line for this case read(luinp,100,err=751,end=2000) hline 100 format(a170) c c Extract storm name and date information from header line read(hline,102) sname,imyear,immon,imday,imtime,ihvmax 102 format(1x,a4,1x,3(i2),1x,i2,1x,i4) cbasin = hline(41:42) csnum = hline(43:44) c c Specify sign of longitude depending on the basin if (cbasin .eq. 'WP' .or. cbasin .eq. 'IO' + .or. cbasin .eq. 'SH') then sgnlon = 1.0 else sgnlon = -1.0 endif c write(lulog,748) sname,imyear,immon,imday,imtime,cbasin,csnum 748 format(/,' HEADER: ',a4,1x,3(i2.2),1x,i2.2,1x,a2,a2) c if (imyear .lt. 50) then imyear4 = 2000 + imyear else imyear4 = 1900 + imyear endif C iline(1) = hline c do k=-2,itmax irsst(k) = 9999 idsst(k) = 9999 insst(k) = 9999 icsst(k) = 9999 enddo idelvm = 9999 c c Read the rest of the data lines for this case c and look for the LAT and LON lines. lcount=2 readlsdloop: do c readlsdloop will fail for read err/end; exit for 'LAST' read(luinp,100,err=751,end=751) tline iline(lcount) = tline c if (tline(157:160) .eq. 'LAT ') then read(tline,105) (islat(k),k=-2,itmax) 105 format(31(1x,i4)) endif C if (tline(157:160) .eq. 'LON ') then read(tline,106) (islon(k),k=-2,itmax),loncon 106 format(31(1x,i4),9x,a1) c c Over-ride sgnlon if longitude convention explictly c appears on the LON line if (loncon .eq. 'W') then sgnlon = -1.0 elseif (loncon .eq. 'E') then sgnlon = 1.0 endif endif c if (tline(157:160) .eq. 'RSST') then read(tline,105) (irsst(k),k=-2,itmax) endif c if (tline(157:160) .eq. 'DSST') then read(tline,105) (idsst(k),k=-2,itmax) endif c if (tline(157:160) .eq. 'NSST') then read(tline,105) (insst(k),k=-2,itmax) endif c if (tline(157:160) .eq. 'CSST') then read(tline,105) (icsst(k),k=-2,itmax) endif c if (tline(157:160) .eq. 'DELV') then read(tline,105) idelvm endif c if (tline(157:160) .eq. 'LAST') exit readlsdloop c lcount=lcount+1 enddo readlsdloop c c Heirarchy of SSTs ssttype='RSST' if (irsst(0) .eq. 9999 .and. idsst(0) .lt. 9999) then do k=-2,itmax irsst(k) = idsst(k) ssttype='DSST' enddo endif if (irsst(0) .eq. 9999 .and. insst(0) .lt. 9999) then do k=-2,itmax irsst(k) = insst(k) ssttype='NSST' enddo endif if (irsst(0) .eq. 9999 .and. icsst(0) .lt. 9999) then do k=-2,itmax irsst(k) = icsst(k) ssttype='CSST' enddo endif write(lulog,302) ssttype 302 format(/,' SST used for KMPI, CFLX, lcmod=',a4) nrawd = nrawd + 1 write(6,200) nrawd,hline 200 format(/,' START CASE',I4,' FIRST LINE OF LUINP HEADER:',/,A80) C do k=-2,itmax slat(k) = float(islat(k))/10.0 slon(k) = sgnlon*float(islon(k))/10.0 c if (irsst(k) .lt. 9000) then rsst(k) = 0.1*float(irsst(k)) else rsst(k) = 9999. endif enddo c c Assign basin if (islat(0) .ne. imiss) then glat = 0.1*float(islat(0)) else glat = rmiss endif c if (islon(0) .ne. imiss) then glon = 0.1*sgnlon*float(islon(0)) else glon = rmiss endif if (glon .lt. 0.0) glon = glon + 360.0 c if (cbasin .eq. 'AL') then ibasin = 1 elseif (cbasin .eq. 'EP' .or. cbasin .eq. 'CP') then ibasin = 2 elseif (cbasin .eq. 'WP') then ibasin = 3 elseif (cbasin .eq. 'IO') then ibasin = 4 elseif (cbasin .eq. 'SH') then ibasin = 5 else call bassel(glat,glon,ibasin) endif c write(6,*) 'cbasin,ibasin: ',cbasin,ibasin c c Initialize synoptic variables to missing initloop: do kt=0,itmax iu200(kt) = 9999 iu200c(kt)= 9999 iv200c(kt)= 9999 it150(kt) = 9999 it200(kt) = 9999 it250(kt) = 9999 it000(kt) = 9999 iz000(kt) = 9999 ir000(kt) = 9999 ie000(kt) = 9999 iepos(kt) = 9999 ieneg(kt) = 9999 iepss(kt) = 9999 ienss(kt) = 9999 ishrd(kt) = 9999 ishtd(kt) = 9999 ishrs(kt) = 9999 ishts(kt) = 9999 ishrg(kt) = 9999 ishxu(kt) = 9999 ishxl(kt) = 9999 ivorb(kt) = 9999 idivb(kt) = 9999 irefc(kt) = 9999 ipefc(kt) = 9999 irhlo(kt) = 9999 irhhi(kt) = 9999 irhmd(kt) = 9999 ipslv(kt) = 9999 c itlat(kt) = 9999 itlon(kt) = 9999 itwac(kt) = 9999 itwbc(kt) = 9999 itwxc(kt) = 9999 irtcc(kt) = 9999 irtxc(kt) = 9999 ipsfc(kt) = 9999 ipsfcc(kt)= 9999 iv000(kt) = 9999 iv850(kt) = 9999 iv500(kt) = 9999 iv300(kt) = 9999 c ishdc(kt) = 9999 isddc(kt) = 9999 ishgc(kt) = 9999 it15c(kt) = 9999 it20c(kt) = 9999 it25c(kt) = 9999 iepoc(kt) = 9999 idivc(kt) = 9999 c itgrd(kt) = 9999 itadv(kt) = 9999 c ivkmpi(kt) = 9999 ipkmpi(kt) = 9999 c ivvavg(kt) = 9999 ivvavc(kt) = 9999 ivmflux(kt)= 9999 c isrh07(kt) = 9999 isrh05(kt) = 9999 c ig250(kt) = 9999 ig200(kt) = 9999 ig150(kt) = 9999 c icflx(kt) = 9999 io500(kt) = 9999 io700(kt) = 9999 c itpw(kt) = imiss itpw0(kt) = imiss do ii=0,mxtpw itpwa(ii,kt) = imiss enddo enddo initloop c ilost=0 tlast=-6.0 c c Start time loop ntsteps = itmax if (izero .eq. 1) ntsteps=0 c timeloop: do kt=0,ntsteps tnow = 6.0*float(kt) c c Check for missing track positions if (slat(kt) .gt. 900.0) cycle timeloop c c Construct rawd file if (iatype .eq. 0) then c Use analyses (re-analysis or operational) for predictors c ktime = kt*idelt c c Calculate date information for next analysis file call tadd(imyear,immon,imday,imtime,ktime, + iayear,iamon,iaday,iatime) c if (imyear4 .le. iray) then rawdfn = 'R00055_X0000_PACK.DAT' else if (cbasin .eq. 'WP' .or. cbasin .eq. 'IO' + .or. cbasin .eq. 'SH') then if (imyear4 .ge. 2005) then rawdfn = 'G00055_X0000_PACK.DAT' else rawdfn = 'H00055_X0000_PACK.DAT' endif else if (imyear4 .ge. 2001) then rawdfn = 'A00055_X0000_PACK.DAT' else rawdfn = 'B00055_X0000_PACK.DAT' endif endif endif c if (iatime .eq. 6 .or. iatime .eq. 18) then rawdfn(8:8) = 'Y' endif c write(rawdfn(5:6),115) iayear 115 format(i2.2) c write(rawdfn(9:10),115) iamon c if (iatime .ge. 12) then iadayt = iaday + 50 else iadayt = iaday endif write(rawdfn(11:12),115) iadayt else c Use forecast fields for the predictors. Try the c forecast initialized at the time of this case. c If that is not available, try a 6-hour old forecast. c if (kt .eq. 0) then ilaghr = 0 ityear = imyear itmon = immon itday = imday ittime = imtime endif c packexistloop: do if (cbasin .eq. 'WP' .or. cbasin .eq. 'IO' + .or. cbasin .eq. 'SH') then rawdfn = 'G00055_X0000_PACK.DAT' else rawdfn = 'A00055_X0000_PACK.DAT' endif c if (ittime .eq. 6 .or. ittime .eq. 18) then rawdfn(8:8) = 'Y' endif c ktime = ilaghr + kt*idelt c write(rawdfn(2:4),110) ktime 110 format(i3.3) c write(rawdfn(5:6),115) ityear c write(rawdfn(9:10),115) itmon c if (ittime .ge. 12) then itdayt = itday + 50 else itdayt = itday endif write(rawdfn(11:12),115) itdayt c if (kt .eq. 0) then c Check to see if first file is available. If not, c look for an earlier model run c call aunpack(ierr,sgnlon) if (ierr .eq. 0) exit packexistloop c ilaghr = ilaghr + 6 c c Set the limit on how far back to look for previous c model runs. if (ilaghr .gt. 12) exit packexistloop c call tadd(imyear,immon,imday,imtime,-ilaghr, + ityear,itmon,itday,ittime) cycle packexistloop else exit packexistloop endif enddo packexistloop endif c c Unpack the current rawd file call aunpack(ierr,sgnlon) c if (ierr .eq. 1 .and. kt .eq. 0) then c If first data file is missing, skip entire case write(lulog,300) 300 format(' First data file must exist, skip case') nskip=nskip+1 cycle caseloop endif c if (ierr .eq. 1) then c If subsequent data files are missing, skip c calculation only for individual time periods cycle timeloop endif c if (igcrop .eq. 1) then call gcrop(slon(kt),slat(kt),bufkm) endif c c Initialize grid parameters at the current c storm location call xyracal(slon(kt),slat(kt)) c c Calculate x,y components of storm motion (m/s) icok = 0 c if (slat(kt-2) .lt. 900.0 .and. + slat(kt ) .lt. 900.0) then t2lat = slat(kt) t2lon = slon(kt) t1lat = slat(kt-2) t1lon = slon(kt-2) deltc = 2.0*delt icok = 1 endif c if (kt+2 .le. itmax .and. icok .eq. 0) then if (slat(kt+2) .lt. 900.0 .and. + slat(kt ) .lt. 900.0 ) then t2lat = slat(kt+2) t2lon = slon(kt+2) t1lat = slat(kt) t1lon = slon(kt) deltc = 2.0*delt icok = 1 endif endif c if (kt+1 .le. itmax .and. icok .eq. 0) then if (slat(kt+1) .lt. 900.0 .and. + slat(kt ) .lt. 900.0 ) then t2lat = slat(kt+1) t2lon = slon(kt+1) t1lat = slat(kt) t1lon = slon(kt) deltc = 1.0*delt icok = 1 endif endif c if (icok .eq. 0) then t2lat = slat(kt) t2lon = slon(kt) t1lat = slat(kt) t1lon = slon(kt) deltc = 1.0*delt endif c call lltoxy(olons,olats,colat,t1lon,t1lat,temx1,temy1) call lltoxy(olons,olats,colat,t2lon,t2lat,temx2,temy2) c cx = 1000.0*(temx2-temx1)/(deltc*3600.0) cy = 1000.0*(temy2-temy1)/(deltc*3600.0) c if (kt .eq. 0) then cx0 = cx cy0 = cy endif c c Get vertical profiles of required variables do n=1,nvp ptem=pvp(n) call varget('U',ptem,fxy,1,iexist) call fcopy(fxy,uxyp(1,1,n)) call aavg(fxy,radinn,radout,minpts,1,ierr,uvp(n)) c call varget('V',ptem,fxy,1,iexist) call fcopy(fxy,vxyp(1,1,n)) call aavg(fxy,radinn,radout,minpts,1,ierr,vvp(n)) c call varget('T',ptem,fxy,1,iexist) call fcopy(fxy,txyp(1,1,n)) call aavg(fxy,radinn,radout,minpts,1,ierr,tvp(n)) c call varget('Z',ptem,fxy,1,iexist) call fcopy(fxy,zxyp(1,1,n)) call aavg(fxy,radinn,radout,minpts,1,ierr,zvp(n)) c if (ptem .ge. 300.0) then call varget('R',ptem,fxy,1,iexist) c if (iatype .eq. 0 .and. ibcrh .eq. 1 + .and. imyear4 .le. iray) then c Apply bias correction to RH call rhbc(fxy,ptem) endif c call fcopy(fxy,rxyp(1,1,n)) call aavg(fxy,radinn,radout,minpts,1,ierr,rvp(n)) else rvp(n) = 999.9 do j=1,ny do i=1,nx rxyp(i,j,n) = 999.9 enddo enddo endif enddo c c c Calculate thetae profile do n=1,nvp tkel = tvp(n) pmb = pvp(n) rh = rvp(n) if (rh .le. 0.0 .or. rh .gt. 99.9) rh = 50.0 call thetae(tkel,pmb,rh,plcl,tlcl,wmr,tevp(n)) qvp(n) = wmr rvp(n) = rh c rhs= 100.0 call thetae(tkel,pmb,rhs,plcl,tlcl,wmr,tevps(n)) enddo c c Save annular average u at 200 mb (in knots) iu200(kt) = ifix(10.0*1.944*uvp(i200)) c c Save annular average T at 150,200 and 250 mb (in deg C) it150(kt) = ifix( 10.0*(tvp(i150)-273.15) ) it200(kt) = ifix( 10.0*(tvp(i200)-273.15) ) it250(kt) = ifix( 10.0*(tvp(i250)-273.15) ) c c Save annular average T,RH and Z at 1000 mb (in deg C) it000(kt) = ifix( 10.0*(tvp(i000)-273.15) ) ir000(kt) = ifix(rvp(i000)) iz000(kt) = ifix(zvp(i000)) c c Calculate and save sea level pressure call psext(tvp(i000),zvp(i000),psfc) ipsfc(kt) = ifix(10.0*(psfc-1000.0)) c c Save annular average low-level RH irhlo(kt) = ifix( (rvp(i700)+rvp(i850))/2.0 ) c c Save annular average mid-level RH irhmd(kt) = ifix( (rvp(i500)+rvp(i700))/2.0 ) c c Save annular average high-level RH irhhi(kt) = ifix( (rvp(i300)+rvp(i400)+rvp(i500))/3.0 ) c c Save annular average thetae at 1000 mb ie000(kt) = ifix( 10.0*tevp(i000) ) c c Calculate and save the thetae positive area epos = 0.0 eneg = 0.0 epss = 0.0 enss = 0.0 do n=1,nvp delte = tevp(i000) - tevp(n) deltes= tevp(i000) - tevps(n) if (delte .gt. 0.0) epos = epos + wvp(n)*delte if (delte .lt. 0.0) eneg = eneg + wvp(n)*abs(delte) if (deltes .gt. 0.0) epss = epss + wvp(n)*deltes if (deltes .lt. 0.0) enss = enss + wvp(n)*abs(deltes) enddo iepos(kt) = ifix( 10.0*epos ) ieneg(kt) = ifix( 10.0*eneg ) iepss(kt) = ifix( 10.0*epss ) ienss(kt) = ifix( 10.0*enss ) c c Calculate and save the 850-200 mb vertical shear shearx = uvp(i200)-uvp(i850) sheary = vvp(i200)-vvp(i850) call ctorh(shearx,sheary,shrd,shtd) ishrd(kt) = ifix(10.0*1.944*shrd) ishtd(kt) = ifix(shtd) c c Calculate and save the 850-500 mb vertical shear shearx = uvp(i500)-uvp(i850) sheary = vvp(i500)-vvp(i850) call ctorh(shearx,sheary,shrs,shts) ishrs(kt) = ifix(10.0*1.944*shrs) ishts(kt) = ifix(shts) c c Calculate and save the generalized shear parameter call gshear(uvp,vvp,wvp,nvp,shrg) ishrg(kt) = ifix(10.0*1.944*shrg) c c Calculate and save the maximum shear in the c upper and lower part of the sounding call xshear(uvp,vvp,pvp,nvp,i150,i500,shxu) ishxu(kt) = ifix(10.0*1.944*shxu) c call xshear(uvp,vvp,pvp,nvp,i500,i000,shxl) ishxl(kt) = ifix(10.0*1.944*shxl) c c Calculate the storm relative helicity 1000-700 hPa call srhel(uvp,vvp,pvp,nvp,cx,cy,i700,i000,isrh07(kt)) if (isrh07(kt) .gt. 9999 .or. isrh07(kt) .lt. -999) then isrh07(kt) = 9999 endif c c Calculate the storm relative helicity 1000-500 hPa call srhel(uvp,vvp,pvp,nvp,cx,cy,i500,i000,isrh05(kt)) if (isrh05(kt) .gt. 9999 .or. isrh05(kt) .lt. -999) then isrh05(kt) = 9999 endif c ipsprt=1 if (kt .eq. 0) then c Calculate and save the steering pressure level alpha = 0.4 if (icok .eq. 1) then call splcal(pvp,uvp,vvp,cx,cy,alpha,nvp,wt1,wt2, + ubard,vbard,pbard,ubars,vbars,pbars) c if (ipsprt .eq. 1) then write(lulog,749) alpha,cx,cy,ubard,vbard,pbard, + ubars,vbars,pbars 749 format(/,' Steering pressure level output', + /,' alpha=',f5.1,' cx,cy: ',2(f5.1,1x), + /,' ub,vb,pb: ',3(f6.1,1x), + /,' us,vs,ps: ',3(f6.1,1x), + //,' p u v wd ws ws/wd', + ' T RH ThE ThEs') c do ll=1,nvp tvptem = tvp(ll)-273.15 write(lulog,750) pvp(ll),uvp(ll),vvp(ll), + wt1(ll),wt2(ll),wt2(ll)/wt1(ll) + ,tvptem,rvp(ll),tevp(ll), + tevps(ll) 750 format(1x,f5.0,1x,f5.1,1x,f5.1,1x, + f6.3,1x,f6.3,1x,f6.2,1x,f5.1,1x,f5.1, + 1x,f5.1,1x,f5.1) enddo endif else pbars = 525.0 endif ipslv( 0) = ifix(pbars) c if (icok .eq. 1) then ipslv( 1) = ifix(sgnf( cx)*0.49 + cx*10.0) ipslv( 2) = ifix(sgnf( cy)*0.49 + cy*10.0) ipslv( 3) = ifix(sgnf(ubard)*0.49 + ubard*10.0) ipslv( 4) = ifix(sgnf(vbard)*0.49 + vbard*10.0) ipslv( 5) = ifix(sgnf(ubars)*0.49 + ubars*10.0) ipslv( 6) = ifix(sgnf(vbars)*0.49 + vbars*10.0) ipslv( 7) = ifix(sgnf(alpha)*0.49 + alpha*100.0) do ll=1,nvp ipslv(7+ll) = ifix(sgnf(wt2(ll))*0.49 + + 1000.0*wt2(ll)) enddo endif c endif c c Calculate area average 850 mb vorticity ptem = 850.0 call voravg(ptem,slon(kt),slat(kt),radvor,vorbar,ierr) c if (ierr .eq. 0) then ivorb(kt) = ifix( (1.0e+7)*vorbar ) else ivorb(kt) = 9999 endif c c Calculate area average 200 mb divergence ptem = 200.0 call divavg(ptem,slon(kt),slat(kt),raddiv,divbar,ierr) c if (ierr .eq. 0) then idivb(kt) = ifix( (1.0e+7)*divbar ) else idivb(kt) = 9999 endif c c Calculate omega (hPa/day) at 500 and 700 hPa ierrp = 0 do ll=1,nvp ptem = pvp(ll) call divavg(ptem,slon(kt),slat(kt),raddiv,divbp(ll),ierr) if (ierr .ne. 0) then divbp(ll) = 9999. ierrp = 1 endif enddo c if (ierrp .eq. 0) then omegap(1) = 0.0 do ll=2,nvp omegap(ll) = omegap(ll-1) - + 0.5*(divbp(ll)+divbp(ll-1))* + 24.0*3600.0*( pvp(ll)- pvp(ll-1)) enddo else do ll=1,nvp omegap(ll) = 9999. enddo endif c io500(kt) = nint(omegap(i500)) if (io500(kt) .gt. 9999 .or. io500(kt) .lt. -999) then io500(k) = 9999 endif io700(kt) = nint(omegap(i700)) if (io700(kt) .gt. 9999 .or. io700(kt) .lt. -999) then io700(k) = 9999 endif c c write(6,768) kt*6 c 768 format(/,'Divergence profile at t=',i3) c c do ll=1,nvp c write(6,769) pvp(ll),divbp(ll),omegap(ll) c 769 format('p=',f6.0,' div=',e11.4,' omega=',e11.4) c enddo c c Calculate relative eddy flux convergence of angular c momentum ptem = 200.0 call refcal(ptem,slon(kt),slat(kt),cx,cy,radrefc,refc) irefc(kt) = ifix(refc) c c Calculate planetary eddy flux convergence of angular c momentum ptem = 200.0 radpefc=1000.0 call pefcal(ptem,slon(kt),slat(kt),radpefc,pefc) ipefc(kt) = ifix(pefc) c c Calculate Joe Cione flux parameter (CFLX) pcore = 9999.0 vmaxjc = float(ihvmax) if (it000(kt) .lt. 9000) then t000jc = 0.1*float(it000(kt)) else t000jc = 9999. endif z000jc = float(iz000(kt)) r000jc = float(ir000(kt)) rlojc = float(irhlo(kt)) rmdjc = float(irhmd(kt)) rhijc = float(irhhi(kt)) rlatjc = abs(slat(kt)) sstjc = rsst(kt) c sspdjc = 1.944*sqrt(cx*cx + cy*cy) if (sspdjc .gt. 80.0) sspdjc = 11.0 c ibiasc = 0 iyearjc = imyear4 c icfcal = 1 imiss =9999 if (it000(kt) .ge. imiss) icfcal = 0 c if (iz000(kt) .ge. imiss) icfcal = 0 if (irhlo(kt) .ge. imiss) icfcal = 0 if (irhmd(kt) .ge. imiss) icfcal = 0 if (irhhi(kt) .ge. imiss) icfcal = 0 if (rlatjc .ge. 100.0) icfcal = 0 if (sstjc .ge. 35.0) icfcal = 0 c if (icfcal .eq. 1) then call cflux(t000jc,z000jc,r000jc,rlojc,rmdjc,rhijc,sstjc, + vmaxjc,pcore,rlatjc,sspdjc,ibiasc,iyearjc, + ibasin,lulog,flux) else flux = 9999. endif c icflx(kt) = nint(flux) c c write(6,767) kt*6,vmaxjc,t000jc,z000jc,r000jc, c + rlojc,rmdjc,rhijc,rlatjc,sspdjc,sstjc, c + flux,icflx(kt) 767 format('t=',i3,' vmax=',f4.0,' t000=',f4.1,' z000=',f5.1, + ' r000=',f4.0,' rhlo=',f4.0,' rhmd=',f4.0,' rhhi=',f4.0, + ' lat=',f5.1,' sspd=',f5.1,' sstjc=',f5.1, + ' cflux=',f6.1,' icflux=',i4) c c Calculate Kerry Emanual MPI sstkmpi = rsst(kt) pskmpi = psfc c c Reverse order of P-level for KMPI routine do n=1,nvp nn = nvp - (n-1) pkmpi(n) = pvp(nn) tkmpi(n) = tvp(nn) - 273.15 qkmpi(n) = qvp(nn) rkmpi(n) = rvp(nn) enddo c ikflag=0 if (sstkmpi .lt. 35.0) then call pcmin(sstkmpi,pskmpi,pkmpi,tkmpi,qkmpi,nvp,nvp, + pmpi,vmpi,ikflag) endif c if (ikflag .eq. 1) then ivkmpi(kt) = ifix(vmpi*1.944 + 0.5) ipkmpi(kt) = ifix(pmpi+0.5) else ivkmpi(kt) = 9999 ipkmpi(kt) = 9999 endif c if (kt .eq. 900) then write(6,*) 'sst,psfc= ',sstkmpi,pskmpi do n=1,nvp write(6,821) pkmpi(n),tkmpi(n),qkmpi(n),rkmpi(n) 821 format(' p,t,q,rh: ',4(f7.1,1x)) enddo c write(6,*) 'pmpi,vmpi: ',ipkmpi(kt),ivkmpi(kt) endif c c Calculate convective mass flux from entraining parcel c cloud model do n=1,nvp pkmpi(n) = pvp(n) tkmpi(n) = tvp(n) - 273.15 rkmpi(n) = rvp(n) enddo c if (sstkmpi .lt. 35.0) then nepx = nvp + 1 pkmpi(nepx) = psfc tkmpi(nepx) = sstkmpi rkmpi(nepx) = rkmpi(nvp) else nepx = nvp endif c c do n=1,nepx c write(6,333) pkmpi(n),tkmpi(n),rkmpi(n) c 333 format('P,T,RH: ',3(f7.1,1x)) c enddo c upx = 15.0 zpx = 0.0 radx = 500.0 c radx = 1000.0 cex = 0.10 alphax = 1.0/600.0 iicex = 1 icweightx = 1 iprtx = 0 call lcmod(pkmpi,tkmpi,rkmpi,nepx,upx,zpx,radx, + cex,alphax,iicex,icweightx,iprtx, + vvavg,vmflux) c if (vvavg .lt. 50.0) then ivvavg(kt) = ifix(0.5 + 100.0*vvavg) else ivvavg(kt) = 9999 endif c if (vmflux .lt. 50.0) then ivmflux(kt) = ifix(0.5 + 100.0*vmflux) else ivmflux(kt) = 9999 endif c if (kt .eq. 0) then write(lulog,822) sstkmpi,pskmpi, + ipkmpi(kt),ivkmpi(kt),ihvmax, + vvavg,slat(kt),slon(kt) 822 format(/,'SST= ',f6.1,' PSFC=',f6.1, + ' PMPI=',i4,' VMPI=',i4,' VMAX=',i4, + ' VVAV=',f7.1,' LAT=',f5.1,' LON=',f6.1) endif c c ++ Start new code for vortex removal and adjusted predictors c c Specify radius (km) for storm-relative avg divergence raddivc = 1000.0 indxdivc = ifix(0.0001 + raddivc/drcf) c c Specific radius (km) for azimuthal average of tangential c wind radaav = 500.0 indxaav = ifix(0.0001 + radaav/drcf) c c Determine if this is a southern hemisphere case if (olats .lt. 0.0) then sgnlat = -1.0 else sgnlat = 1.0 endif c c First guess location for vortex search if (kt .eq. 0) then slatcf0 = slat(kt) sloncf0 = slon(kt) else slatcf0 = slatcf sloncf0 = sloncf endif c call lltoxy(olons,olats,colat,sloncf0,slatcf0,xcf0,ycf0) c c Find 850 mb circulation center call cfindxy(uxyp(1,1,i850),vxyp(1,1,i850),rcf,thetacf, + mrcf,mtcf,nrcf,ntcf, + radtwa,xcf0,ycf0,xcf,ycf,twa,rwa) c call xytoll(olons,olats,colat,xcf,ycf,sloncf,slatcf) slongw = sloncf slatgw = slatcf c itlat(kt) = ifix( 10.0*slatcf) itlon(kt) = ifix(sgnlon*10.0*sloncf) c c Check for lost vortex rdist=sqrt( (xcf-xcf0)**2 + (ycf-ycf0)**2 ) spdkt = (60.0/111.1)*rdist/(tnow-tlast) tlast = tnow c call spdlim(slatcf0,spdmax) c if (ilost .eq. 0) then if (spdkt .gt. spdmax) ilost=1 if (twa .lt. 1.0) ilost=1 endif c if (ilost .eq. 0) then itwac(kt) = ifix( 10.0*twa*sgnlat) else itwac(kt) = 9999 endif c c write(6,871) kt*6,slat(kt),slon(kt), c + slatcf0,sloncf0,slatcf,sloncf,twa,rwa, c + spdkt,spdmax,ilost c 871 format('t=',i3,' blat,blon: ',f6.2,1x,f7.2, c + ' lat0,lon0: ',f6.2,1x,f7.2, c + ' lat,lon: ' ,f6.2,1x,f7.2, c + ' t,rwa: ',f6.1,1x,f6.1, c + ' s,sm,ilost: ',f5.1,1x,f5.1,1x,i1) c if (ilost .eq. 0) then c Filter thermodynamic variables and c remove vortex from wind field c fradlf = 400.0 irmflag = 1 c do ll=1,nvp c call lfilter(txyp(1,1,ll),xcf,ycf,fradlf) c call lfilter(rxyp(1,1,ll),xcf,ycf,fradlf) c call lfilter(zxyp(1,1,ll),xcf,ycf,fradlf) c call vremove(uxyp(1,1,ll),vxyp(1,1,ll),zxyp(1,1,ll), + txyp(1,1,ll),rcf,thetacf, + xcf,ycf,mrcf,mtcf,nrcf,ntcf, + uraa,vtaa,zaa,taa,irmflag) c c Save azimuthally averaged radial and tangential winds, c t and z do i=0,nrcf uraap(i,ll) = uraa(i) vtaap(i,ll) = vtaa(i) taap(i,ll) = taa(i) zaap(i,ll) = zaa(i) enddo c if (ll .eq. i200) then c Calculate storm-relative 200 mb area average c divergence divc = 2.0*uraa(indxdivc)/(1000.0*raddivc) idivc(kt) = ifix(divc*1.0e+7) endif c if (ll .eq. i850) then c Calculate maximum and average cyclonic tangential c wind c c Find last radial point where vtaa .ge. 1 starting c at n=4 do kk=4,nrcf nrlast = kk vtaatem = vtaa(kk)*sgnlat if (vtaatem .lt. 1.0) then nrlast = kk-1 exit endif enddo c rtlast = rcf(nrlast) c vtmax = -1000.0 rtmax = 1000.0 vtavg = 0.0 do kk=0,nrlast vtaatem = vtaa(kk)*sgnlat vtavg = vtavg + vtaa(kk) if (vtaatem .gt. vtmax) then vtmax = vtaatem rtmax = rcf(kk) endif enddo vtavg = vtavg/(float(nrlast+1)) vtmax = vtmax*sgnlat c irtcc(kt) = ifix(rtlast) irtxc(kt) = ifix(rtmax) itwxc(kt) = ifix(vtmax*10.0) itwbc(kt) = ifix(vtavg*10.0) endif c if (ll .eq. i000) then c Save azimuthally averaged 1000 hPa tangential wind c at r=radaav iv000(kt) = ifix(vtaa(indxaav)*10.0) endif c if (ll .eq. i850) then c Save azimuthally averaged 850 hPa tangential wind c at r=radaav iv850(kt) = ifix(vtaa(indxaav)*10.0) endif c if (ll .eq. i500) then c Save azimuthally averaged 500 hPa tangential wind c at r=radaav iv500(kt) = ifix(vtaa(indxaav)*10.0) endif c if (ll .eq. i300) then c Save azimuthally averaged 300 hPa tangential wind c at r=radaav iv300(kt) = ifix(vtaa(indxaav)*10.0) endif c c do kk=0,nrcf c write(6,873) tnow,pvp(ll),rcf(kk),uraa(kk),vtaa(kk) c 873 format(f5.0,1x,f5.0,1x,f5.0,1x,f6.1,1x,f6.1) c enddo c enddo else c Vortex lost, so skip the filtering, but still calculate c the iv000, iv850, iv500 and iv300 variables. c irmftem = 0 call lltoxy(olons,olats,colat,slon(kt),slat(kt), + xcftem,ycftem) slongw = slon(kt) slatgw = slat(kt) c do ll=1,nvp call vremove(uxyp(1,1,ll),vxyp(1,1,ll),zxyp(1,1,ll), + txyp(1,1,ll),rcf,thetacf, + xcftem,ycftem,mrcf,mtcf,nrcf,ntcf, + uraa,vtaa,zaa,taa,irmftem) c Save azimuthally averaged radial and tangential winds, c T and Z do i=0,nrcf uraap(i,ll) = uraa(i) vtaap(i,ll) = vtaa(i) taap(i,ll) = taa(i) zaap(i,ll) = zaa(i) enddo c if (ll .eq. i000) then c Save azimuthally averaged 1000 hPa tangential wind c at r=radaav iv000(kt) = ifix(vtaa(indxaav)*10.0) endif c if (ll .eq. i850) then c Save azimuthally averaged 850 hPa tangential wind c at r=radaav iv850(kt) = ifix(vtaa(indxaav)*10.0) endif c if (ll .eq. i500) then c Save azimuthally averaged 500 hPa tangential wind c at r=radaav iv500(kt) = ifix(vtaa(indxaav)*10.0) endif c if (ll .eq. i300) then c Save azimuthally averaged 300 hPa tangential wind c at r=radaav iv300(kt) = ifix(vtaa(indxaav)*10.0) endif c enddo endif c c Calculate temperautre associated with symmetric tangential c wind using gradient thermal wind equation. call gtwind(vtaap,taap,gtaap,rcf,pvp,slatgw,mrcf,nrcf,nvp) c iuvprt = 0 if (iuvprt .eq. 1 .and. kt .eq. 0) then write(6,655) kt*6 655 format(/,'Azimuthal Avg Radial Wind (m/s) at t=',i4) write(6,656) (pvp(ll),ll=nvp,1,-1) 656 format('r (km) p=',10(f8.0)) c do i=0,nrcf write(6,657) rcf(i),(uraap(i,ll),ll=nvp,1,-1) 657 format(f6.0,5x,15(f8.1)) enddo c write(6,658) kt*6 658 format(/,'Azimuthal Avg Tangential Winds (m/s) at t=',i4) write(6,656) (pvp(ll),ll=nvp,1,-1) c do i=0,nrcf write(6,657) rcf(i),(vtaap(i,ll),ll=nvp,1,-1) enddo c write(6,659) 659 format(/,'Azimuthal Avg Z Pert (m) at t=',i4) write(6,656) (pvp(ll),ll=nvp,1,-1) c do i=0,nrcf write(6,657) rcf(i),((zaap( i,ll)- + zaap(nrcf,ll)),ll=nvp,1,-1) enddo c write(6,660) kt*6 660 format(/,'Azimuthal Avg Temp Pert (K) at t=',i4) write(6,656) (pvp(ll),ll=nvp,1,-1) c do i=0,nrcf write(6,657) rcf(i),((taap( i,ll)- + taap(nrcf,ll)),ll=nvp,1,-1) enddo c write(6,661) kt*6 661 format(/,'Azimuthal Avg GTW Temp Pert (K) at t=',i4) write(6,656) (pvp(ll),ll=nvp,1,-1) c do i=0,nrcf c write(6,657) rcf(i),(gtaap(i,ll),ll=nvp,1,-1) write(6,657) rcf(i),((gtaap( i,ll)- + gtaap(nrcf,ll)),ll=nvp,1,-1) enddo c write(6,*) endif c c Average T of associated with symmetric vortex averaged c over area used for the T calculation ipert = 20 call faavg(slon(kt),slat(kt),slongw,slatgw,radinn,radout, + rcf,gtaap(1,i250),mrcf,ipert,tgavg250,ierrg) c call faavg(slon(kt),slat(kt),slongw,slatgw,radinn,radout, + rcf,gtaap(1,i200),mrcf,ipert,tgavg200,ierrg) c call faavg(slon(kt),slat(kt),slongw,slatgw,radinn,radout, + rcf,gtaap(1,i150),mrcf,ipert,tgavg150,ierrg) c if (iuvprt .eq. 1) then write(6,876) ktime,ilost,slat(kt),slon(kt),slatgw,slongw, + tgavg250,tgavg200,tgavg150 876 format(i3,' lost=',i1,' slat/lon',f6.1,1x,f7.1, + ' vlat/lon',f6.1,1x,f7.1, + ' tg250,200,150',3(f6.1,1x)) endif c if (ilost .ne. 0) then tgavg250 = 0.0 tgavg200 = 0.0 tgavg150 = 0.0 endif c ig250(kt) = nint(10.0*tgavg250) ig200(kt) = nint(10.0*tgavg200) ig150(kt) = nint(10.0*tgavg150) c c Calculate surface pressure at edge of tangential wind area dxtemp = xgk(2,1)-xgk(1,1) dytemp = ygk(1,2)-ygk(1,1) x1temp = xgk(1,1) y1temp = ygk(1,1) c if (ilost .eq. 0) then xcftemp = xcf ycftemp = ycf nrlastt = nrlast else xcftemp = 0.0 ycftemp = 0.0 nrlastt = 16 endif c do j=0,ntcf do i=0,nrcf xtemp = xcftemp + rcf(i)*ctheta(j) ytemp = ycftemp + rcf(i)*stheta(j) c call lintcf(zxyp(1,1,i000),x1temp,dxtemp, + y1temp,dytemp, + ixmx,iymx,nx,ny, + xtemp,ytemp,zcrt(i,j)) call lintcf(txyp(1,1,i000),x1temp,dxtemp, + y1temp,dytemp, + ixmx,iymx,nx,ny, + xtemp,ytemp,tcrt(i,j)) enddo enddo c call azavgcf(zcrt,zraa,mrcf,mtcf,nrcf,ntcf) call azavgcf(tcrt,traa,mrcf,mtcf,nrcf,ntcf) c c do i=0,nrcf c write(6,*) kt,rcf(i),zraa(i),traa(i) c enddo c z00 = zraa(nrlastt) t00 = traa(nrlastt) call psext(t00,z00,psfc) ipsfcc(kt) = ifix(10.0*(psfc-1000.0)) c ipsfcc(kt) = ifix(psfc+0.5) c c write(6,889) z00,t00,ipsfcc(kt),ipsfc(kt),nrlastt c 889 format('z00,t00,psfcc,psfc,nrlastt: ', c + 2(f7.1,1x),i4,1x,i4,1x,i2) c c Azimuthally average u,v,t,r,z c do ll=1,nvp call aavg(uxyp(1,1,ll),radinnc,radoutc,minpts,1, + ierr,uvp(ll)) call aavg(vxyp(1,1,ll),radinnc,radoutc,minpts,1, + ierr,vvp(ll)) call aavg(txyp(1,1,ll),radinnc,radoutc,minpts,1, + ierr,tvp(ll)) call aavg(zxyp(1,1,ll),radinnc,radoutc,minpts,1, + ierr,zvp(ll)) if (pvp(ll) .ge. 300.0) then call aavg(rxyp(1,1,ll),radinnc,radoutc,minpts,1, + ierr,rvp(ll)) else rvp(ll) = 999.9 endif c write(6,875) pvp(ll),uvp(ll),vvp(ll),tvp(ll),rvp(ll) c 875 format('p,u,v,t,r: ',5(f7.1,1x)) enddo c c Calculate and save the 850-500 mb vertical shear shearx = uvp(i200)-uvp(i850) sheary = vvp(i200)-vvp(i850) call ctorh(shearx,sheary,shdc,sddc) ishdc(kt) = ifix(10.0*1.944*shdc) isddc(kt) = ifix(sddc) c c Calculate and save the generalized shear parameter call gshear(uvp,vvp,wvp,nvp,shgc) ishgc(kt) = ifix(10.0*1.944*shgc) c c Calculate and save the 850 to 700 hPa temperature gradient c and advection call tgrad(slat(kt),pvp(i700),uvp(i700),vvp(i700), + pvp(i850),uvp(i850),vvp(i850), + tgrd,tadv,ierr) c if (ierr .eq. 0) then itgrd(kt) = ifix(tgrd*1.0e+7) itadv(kt) = ifix(tadv*1.0e+6) else itgrd(kt) = 9999 itadv(kt) = 9999 endif c c Calculate convective mass flux from entraining parcel c cloud model with corrected T/RH profiles do n=1,nvp pkmpi(n) = pvp(n) tkmpi(n) = tvp(n) - 273.15 rkmpi(n) = rvp(n) enddo c if (sstkmpi .lt. 35.0) then nepx = nvp + 1 pkmpi(nepx) = psfc tkmpi(nepx) = sstkmpi rkmpi(nepx) = rkmpi(nvp) else nepx = nvp endif c call lcmod(pkmpi,tkmpi,rkmpi,nepx,upx,zpx,radx, + cex,alphax,iicex,icweightx,iprtx, + vvavg,vmflux) c if (vvavg .lt. 50.0) then ivvavc(kt) = ifix(0.5 + 100.0*vvavg) else ivvavc(kt) = 9999 endif c c Save annular average T at 150,200 and 250 mb (in deg C) it15c(kt) = ifix( 10.0*(tvp(i150)-273.15) ) it20c(kt) = ifix( 10.0*(tvp(i200)-273.15) ) it25c(kt) = ifix( 10.0*(tvp(i250)-273.15) ) c c Save u200c and v200c iu200c(kt) = ifix(10.0*1.944*uvp(i200)) iv200c(kt) = ifix(10.0*1.944*vvp(i200)) c c Calculate thetae profile do n=1,nvp tkel = tvp(n) pmb = pvp(n) rh = rvp(n) if (rh .le. 0.0 .or. rh .gt. 99.9) rh = 50.0 call thetae(tkel,pmb,rh,plcl,tlcl,wmr,tevp(n)) rhs= 100.0 call thetae(tkel,pmb,rhs,plcl,tlcl,wmr,tevps(n)) enddo c c Calculate and save the thetae positive area epos = 0.0 eneg = 0.0 epss = 0.0 enss = 0.0 do n=1,nvp delte = tevp(i000) - tevp(n) deltes= tevp(i000) - tevps(n) if (delte .gt. 0.0) epos = epos + wvp(n)*delte if (delte .lt. 0.0) eneg = eneg + wvp(n)*abs(delte) if (deltes .gt. 0.0) epss = epss + wvp(n)*deltes if (deltes .lt. 0.0) enss = enss + wvp(n)*abs(deltes) enddo iepoc(kt) = ifix( 10.0*epos ) c c ++ End new code for vortex removal and adjusted predictors c c ++ Start TPW calculation ++ c c if (kt .eq. 0) then if (kt .le. itmax) then c Fill surface RH array (assume it equals RH at 1000 hPa) c and calculate surface pressure from z(p=1000 hPa) do j=1,ny do i=1,nx rsxy(i,j) = rxyp(i,j,nvp) call psext(txyp(i,j,i000),zxyp(i,j,i000),psxy(i,j)) c do k=1,nvp tptpw(k) = txyp(i,j,k) rptpw(k) = rxyp(i,j,k) pptpw(k) = pvp(k) enddo c tstpw = txyp(i,j,i000) + 0.0065*(zxyp(i,j,i000)+111.0) rstpw = rsxy(i,j) pstpw = psxy(i,j) c call tpwcal(pptpw,tptpw,rptpw,nvp, + pstpw,tstpw,rstpw,tpwxy(i,j)) c c if (i .eq. 10 .and. j .eq. 10) then c write(6,*) ' t, z,ps: ',txyp(i,j,i000), c + zxyp(i,j,i000),psxy(i,j) c write(6,*) 'ts,rs,ps,tpw: ',tstpw,rstpw,pstpw, c + tpwxy(i,j) c endif enddo enddo c c Calculate TPW parameters call ctor(cx,cy,temp,cdir00) c rtemp = float(isddc(kt)) if (rtemp .ge. rmiss) then sdir00 = 0.0 else if (rtemp .gt. 180.0) rtemp = rtemp - 360.0 sdir00 = 90.0-rtemp if (sdir00 .gt. 360.0) sdir00 = sdir00 - 360.0 if (sdir00 .lt. 0.0) sdir00 = sdir00 + 360.0 endif c c write(6,*) 'cdir00,sdir00: ',cdir00,sdir00 call tpwpcal(tpwxy,ixmx,iymx,nx,ny, + rlond(1),rlond(nx),dlon, + rlatd(1),rlatd(ny),dlat, + slat(kt),slon(kt),cdir00,sdir00, + itpw,itmax) c if (kt .eq. 0) then do ii=0,mxtpw itpw0(ii) = itpw(ii) enddo endif c do ii=0,mxtpw itpwa(ii,kt) = itpw(ii) enddo endif c c ++ End TPW calculation ++ enddo timeloop c c ++ Start code to calculate proxy GOES predictors c for use when GOES data c is missing. c grsst = rsst(0) gvmax = float(ihvmax) c if (ishdc(0) .ne. imiss) then gshdc = 0.1*float(ishdc(0)) else gshdc = rmiss endif c if (idelvm .ne. imiss) then gdelv = float(idelvm) else gdelv = 0.0 endif c if (it200(0) .ne. imiss) then gt200 = 0.1*it200(0) else gt200 = rmiss endif c if (idivb(0) .ne. imiss) then gd200 = idivb(0) else gd200 = rmiss endif c gspd = 1.944*sqrt(cx0*cx0 + cy0*cy0) c call ir00param4(ibasin,glat,grsst,gvmax,gshdc,gdelv,gt200, + gd200,gspd,rmiss,pgoes,ierr) c if (ierr .ne. 0) then do i=0,itmax ipgoes(i) = imiss enddo else ipgoes(0) = 0 do i=1,19 ipgoes(i) = ifix(pgoes(i)) enddo do i=20,itmax ipgoes(i) = imiss enddo endif c c write(6,817) glat,glon,grsst,gvmax,gshdc,gdelv,gt200,gd200, c + gspd,ibasin c 817 format('ginput ',9(f7.1),1x,i2) c c ++ End code to calculate proxy GOES predictors c c Write input lines to output file (except the last line) do i=1,lcount-1 write(luout,100) iline(i) enddo c c Set iskip=1 to skip all but a few variables iskip=0 c c Write out synoptic variables for this case if (iskip .eq. 0) then write(luout,210) (iu200(k),k=0,itmax) 210 format(10x,29(1x,i4),1x,'U200') write(luout,208) (iu200c(k),k=0,itmax) 208 format(10x,29(1x,i4),1x,'U20C') write(luout,209) (iv200c(k),k=0,itmax) 209 format(10x,29(1x,i4),1x,'V20C') write(luout,212) (ie000(k),k=0,itmax) 212 format(10x,29(1x,i4),1x,'E000') write(luout,213) (iepos(k),k=0,itmax) 213 format(10x,29(1x,i4),1x,'EPOS') write(luout,214) (ieneg(k),k=0,itmax) 214 format(10x,29(1x,i4),1x,'ENEG') write(luout,215) (iepss(k),k=0,itmax) 215 format(10x,29(1x,i4),1x,'EPSS') write(luout,216) (ienss(k),k=0,itmax) 216 format(10x,29(1x,i4),1x,'ENSS') write(luout,217) (irhlo(k),k=0,itmax) 217 format(10x,29(1x,i4),1x,'RHLO') write(luout,218) (irhmd(k),k=0,itmax) 218 format(10x,29(1x,i4),1x,'RHMD') write(luout,219) (irhhi(k),k=0,itmax) 219 format(10x,29(1x,i4),1x,'RHHI') c write(luout,231) (ipslv(k),k=0,itmax) 231 format(10x,29(1x,i4),1x,'PSLV') write(luout,232) (ivorb(k),k=0,itmax) 232 format(10x,29(1x,i4),1x,'Z850') write(luout,233) (idivb(k),k=0,itmax) 233 format(10x,29(1x,i4),1x,'D200') c tlabel='REFC' write(luout,236) (irefc(k),k=0,itmax),tlabel tlabel='PEFC' write(luout,236) (ipefc(k),k=0,itmax),tlabel c tlabel='T000' write(luout,236) (it000(k),k=0,itmax),tlabel tlabel='R000' write(luout,236) (ir000(k),k=0,itmax),tlabel tlabel='Z000' write(luout,236) (iz000(k),k=0,itmax),tlabel c write(luout,239) (itlat(k),k=0,itmax) 239 format(10x,29(1x,i4),1x,'TLAT') write(luout,240) (itlon(k),k=0,itmax) 240 format(10x,29(1x,i4),1x,'TLON') write(luout,241) (itwac(k),k=0,itmax) 241 format(10x,29(1x,i4),1x,'TWAC') c write(luout,260) (itwbc(k),k=0,itmax) c 260 format(10x,29(1x,i4),1x,'TWBC') write(luout,261) (itwxc(k),k=0,itmax) 261 format(10x,29(1x,i4),1x,'TWXC') c tlabel = 'G150' write(luout,236) (ig150(k),k=0,itmax),tlabel tlabel = 'G200' write(luout,236) (ig200(k),k=0,itmax),tlabel tlabel = 'G250' write(luout,236) (ig250(k),k=0,itmax),tlabel c tlabel='V000' write(luout,236) (iv000(k),k=0,itmax),tlabel tlabel='V850' write(luout,236) (iv850(k),k=0,itmax),tlabel tlabel='V500' write(luout,236) (iv500(k),k=0,itmax),tlabel tlabel='V300' write(luout,236) (iv300(k),k=0,itmax),tlabel tlabel='TGRD' write(luout,236) (itgrd(k),k=0,itmax),tlabel tlabel='TADV' write(luout,236) (itadv(k),k=0,itmax),tlabel c c write(luout,262) (irtcc(k),k=0,itmax) c 262 format(10x,29(1x,i4),1x,'RTCC') c write(luout,263) (irtxc(k),k=0,itmax) c 263 format(10x,21(1x,i4),1x,'RTXC') write(luout,264) (ipsfcc(k),k=0,itmax) 264 format(10x,29(1x,i4),1x,'PENC') write(luout,242) (ishdc(k),k=0,itmax) 242 format(10x,29(1x,i4),1x,'SHDC') write(luout,243) (isddc(k),k=0,itmax) 243 format(10x,29(1x,i4),1x,'SDDC') write(luout,244) (ishgc(k),k=0,itmax) 244 format(10x,29(1x,i4),1x,'SHGC') c write(luout,245) (it15c(k),k=0,itmax) c 245 format(10x,29(1x,i4),1x,'T15C') c write(luout,246) (it20c(k),k=0,itmax) c 246 format(10x,29(1x,i4),1x,'T20C') c write(luout,247) (it25c(k),k=0,itmax) c 247 format(10x,29(1x,i4),1x,'T25C') c write(luout,248) (iepoc(k),k=0,itmax) c 248 format(10x,29(1x,i4),1x,'EPOC') write(luout,249) (idivc(k),k=0,itmax) 249 format(10x,29(1x,i4),1x,'DIVC') c endif c write(luout,311) (it150(k),k=0,itmax) 311 format(10x,29(1x,i4),1x,'T150') write(luout,312) (it200(k),k=0,itmax) 312 format(10x,29(1x,i4),1x,'T200') write(luout,313) (it250(k),k=0,itmax) 313 format(10x,29(1x,i4),1x,'T250') write(luout,220) (ishrd(k),k=0,itmax) 220 format(10x,29(1x,i4),1x,'SHRD') write(luout,222) (ishtd(k),k=0,itmax) 222 format(10x,29(1x,i4),1x,'SHTD') write(luout,225) (ishrs(k),k=0,itmax) 225 format(10x,29(1x,i4),1x,'SHRS') write(luout,227) (ishts(k),k=0,itmax) 227 format(10x,29(1x,i4),1x,'SHTS') write(luout,228) (ishrg(k),k=0,itmax) 228 format(10x,29(1x,i4),1x,'SHRG') c if (imxshr .eq. 1) then write(luout,229) (ishxu(k),k=0,itmax) 229 format(10x,29(1x,i4),1x,'SHXU') write(luout,230) (ishxl(k),k=0,itmax) 230 format(10x,29(1x,i4),1x,'SHXL') endif c tlabel='PENV' write(luout,236) (ipsfc(k),k=0,itmax),tlabel 236 format(10x,29(1x,i4),1x,a4) c tlabel='VMPI' write(luout,236) (ivkmpi(k),k=0,itmax),tlabel c tlabel='PMPI' c write(luout,236) (ipkmpi(k),k=0,itmax),tlabel c tlabel='VVAV' write(luout,236) (ivvavg(k),k=0,itmax),tlabel tlabel='VMFX' write(luout,236) (ivmflux(k),k=0,itmax),tlabel tlabel='VVAC' write(luout,236) (ivvavc(k),k=0,itmax),tlabel c tlabel='HE07' write(luout,236) (isrh07(k),k=0,itmax),tlabel tlabel='HE05' write(luout,236) (isrh05(k),k=0,itmax),tlabel c tlabel= 'O500' write(luout,236) (io500(k),k=0,itmax),tlabel tlabel= 'O700' write(luout,236) (io700(k),k=0,itmax),tlabel c tlabel= 'CFLX' write(luout,236) (icflx(k),k=0,itmax),tlabel c tlabel='MTPW' write(luout,236) (itpw0(k),k=0,itmax),tlabel c do ii=0,ntpw write(tlabel,237) ii+1 237 format('PW',i2.2) write(luout,236) (itpwa(ii,k),k=0,itmax),tlabel enddo c tlabel='IRXX' write(luout,236) (ipgoes(k),k=0,itmax),tlabel c c Write out last line write(luout,100) iline(lcount) c C Do next case enddo caseloop C 2000 continue C WRITE(LULOG,350) nrawd,nskip WRITE(6,350) nrawd,nskip 350 FORMAT(/,' NORMAL COMPLETION OF LSDIAG, ',I4,' CASES PROCESSED ', + I4,' SKIPPED') return c c Error processing 751 write(lulog,*) ' Error reading input file (luinp)' stop C return END subroutine sv7cal c subroutine psext(t00,z00,psfc) c This routine calculates the surface pressure in hPa c by extrapolating from the 1000 mb height surface. c c Input: t00 = temperature (K) at 1000 mb c z00 = height deviation (m) of 1000 mb from the height of c the standard atmosphere. c gamma = 0.0065 aa = 9.81/(287.0*gamma) p00 = 1000.0 call stndz(p00,zb00,tb00,thb00) z00t= zb00+z00 t11 = t00 + gamma*z00t c psfc= p00*( (t11/t00)**aa ) c return END subroutine psext c subroutine pncal c This routine specifies physical and numerical constants c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c pi = 3.14159265 dtr = pi/180.0 c dtk = 111.1 c erad = 6371.0 erot = 7.292E-5 c rd = 287.0 c return END subroutine pncal c subroutine refcal(ptem,rlon,rlat,cx,cy,radmax,refc) c This routine calculates the relative eddy c flux convergence (m/s/day). The refc is calculated c in storm relative coordinates, and is radially c averaged. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Local variables dimension u(ixmx,iymx),v(ixmx,iymx),ur(ixmx,iymx),vt(ixmx,iymx) dimension vtp(irmx,iamx),urp(irmx,iamx) dimension upvp(irmx,iamx) dimension urpb(irmx),vtpb(irmx) dimension upvpb(irmx),dupvpb(irmx),rrefc(irmx) c ierr=0 c c Get u and v ifatal=0 call varget('U',ptem,u,ifatal,iexu) call varget('V',ptem,v,ifatal,iexv) c if (iexu .ne. 1 .or. iexv .ne. 1) then ierr=1 refc = 0.0 return endif c c Convert u and v to storm relative coordinates do j=1,ny do i=1,nx u(i,j) = u(i,j) - cx v(i,j) = v(i,j) - cy enddo enddo c c Calculate radial and tangential winds call xyracal(rlon,rlat) call uvtort(u,v,ur,vt) c c Interpolate radial and tangential winds to polar grid call ctop(ur,urp) call ctop(vt,vtp) c c Azimuthally average ur and vt call azavg(urp,urpb,0) call azavg(vtp,vtpb,0) c c Calculate U'V' do j=1,ntheta do i=1,nrad upvp(i,j) = (urp(i,j) - urpb(i))* + (vtp(i,j) - vtpb(i)) enddo enddo c c Azimuthally average U'V' call azavg(upvp,upvpb,0) c c Take radial derivative of averaged U'V'. call ddrcal(upvpb,dupvpb) c c Calculate relative flux convergence in m/s/day cf = -24.0*3600.0 rrefc(1) = cf*dupvpb(1) do i=2,nrad rrefc(i) = cf*(2.*upvpb(i)/(1000.*radk(i)) + dupvpb(i)) enddo c refc = 0.0 icount = 0 do i=2,nrad if (radk(i) .gt. radmax) exit refc = refc + rrefc(i) icount = icount + 1 enddo c if (icount .ge. 1) then refc = refc/float(icount) else refc = 0.0 endif c return END subroutine refcal c subroutine pefcal(ptem,rlon,rlat,radmax,pefc) c This routine calculates the planetary eddy c flux convergence (m/s/day). The pefc is radially averaged. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c c Local variables dimension u(ixmx,iymx),v(ixmx,iymx),ur(ixmx,iymx),vt(ixmx,iymx) dimension urp(irmx,iamx) dimension urpb(irmx) c c Local variables for Coriolis calculation dimension f(irmx,iamx),fb(irmx) dimension t1ra(irmx,iamx),t1r(irmx) c ierr=0 c c Get u and v ifatal=0 call varget('U',ptem,u,ifatal,iexu) call varget('V',ptem,v,ifatal,iexv) c if (iexu .ne. 1 .or. iexv .ne. 1) then ierr=1 pefc = 0.0 return endif c c Calculate radial and tangential winds call xyracal(rlon,rlat) call uvtort(u,v,ur,vt) c c Interpolate radial winds to polar grid call ctop(ur,urp) c c Azimuthally average ur call azavg(urp,urpb,0) c c Calculate Coriolis parameter at all of the polar grid points do j=1,ntheta do i=1,nrad xtemp = xra(i,j) ytemp = yra(i,j) call xytoll(olons,olats,colat,xtemp,ytemp,tlon,tlat) f(i,j) = 2.0*erot*sin(dtr*tlat) enddo enddo c c Azimuthally average Coriolis parameter call azavg(f,fb,1) c c Calculate planetary eddy flux momentum convergence c in units of m/s/day do j=1,ntheta do i=1,nrad t1ra(i,j) = -(urp(i,j)-urpb(i))*(f(i,j)-fb(i)) enddo enddo call azavg(t1ra,t1r,0) c do i=1,nrad t1r(i) = t1r(i)*24.0*3600.0 enddo c pefc = 0.0 icount = 0 do i=2,nrad if (radk(i) .gt. radmax) exit pefc = pefc + t1r(i) icount = icount + 1 enddo c if (icount .ge. 1) then pefc = pefc/float(icount) else pefc = 0.0 endif c return END subroutine pefcal c subroutine divavg(ptem,rlon,rlat,radtem,divbar,ierr) c This routine calculates the area averaged c divergence at level ptem by azimuthally averaging c the radial wind. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) c C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Local variables dimension u(ixmx,iymx),v(ixmx,iymx),ur(ixmx,iymx),vt(ixmx,iymx) dimension urp(irmx,iamx),urpb(irmx) c ierr=0 c c Get u and v ifatal=0 call varget('U',ptem,u,ifatal,iexu) call varget('V',ptem,v,ifatal,iexv) c if (iexu .ne. 1 .or. iexv .ne. 1) then ierr=1 divbar=0.0 return endif c c Calculate radial and tangential winds call xyracal(rlon,rlat) call uvtort(u,v,ur,vt) c c Interpolate radial wind to polar grid call ctop(ur,urp) c c Azimuthally average vtp call azavg(urp,urpb,0) c c Find radial grid point closest to requested radius irad = 0 do i=3,nrad if (radk(i) .gt. radtem) then irad = i-1 exit endif enddo c if (irad .le. 0) then ierr=1 divbar=0.0 return endif c divbar = 2.0*urpb(irad)/(1000.0*radk(irad)) c return END subroutine divavg c subroutine voravg(ptem,rlon,rlat,radtem,vorbar,ierr) c This routine calculates the area averaged relative c vorticity at level ptem by azimuthally averaging c the tangential wind. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) c C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Local variables dimension u(ixmx,iymx),v(ixmx,iymx),ur(ixmx,iymx),vt(ixmx,iymx) dimension vtp(irmx,iamx),vtpb(irmx) c ierr=0 c c Get u and v ifatal=0 call varget('U',ptem,u,ifatal,iexu) call varget('V',ptem,v,ifatal,iexv) c if (iexu .ne. 1 .or. iexv .ne. 1) then ierr=1 vorbar=0.0 return endif c c Calculate radial and tangential winds call xyracal(rlon,rlat) call uvtort(u,v,ur,vt) c c Interpolate tangential wind to polar grid call ctop(vt,vtp) c c Azimuthally average vtp call azavg(vtp,vtpb,0) c c Find radial grid point closest to requested radius irad = 0 do i=3,nrad if (radk(i) .gt. radtem) then irad = i-1 exit endif enddo c if (irad .le. 0) then ierr=1 vorbar=0.0 return endif c vorbar = 2.0*vtpb(irad)/(1000.0*radk(irad)) c return END subroutine voravg c subroutine defcal(ptem,defxy,ierr) c This routine calculates the horizontal deformation at pressure c level ptem. The grid is defined by the last call to xyracal. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension defxy(ixmx,iymx) c c Local variables dimension uxy(ixmx,iymx),vxy(ixmx,iymx) dimension dudx(ixmx,iymx),dudy(ixmx,iymx) dimension dvdx(ixmx,iymx),dvdy(ixmx,iymx) c c Get u and v call varget('U',ptem,uxy,0,iexu) call varget('V',ptem,vxy,0,iexv) if (iexu .ne. 1 .or. iexv .ne. 1) then ierr = 1 return endif c c calculate required derivatives call dxcal(uxy,dudx) call dycal(uxy,dudy) call dxcal(vxy,dvdx) call dycal(vxy,dvdy) c c Calculate the magnitude of the horizontal deformation vector do j=1,ny do i=1,nx defxy(i,j) = sqrt( (dudy(i,j)+dvdx(i,j))**2 + + (dudx(i,j)-dvdy(i,j))**2 ) enddo enddo c ierr = 0 c return end subroutine defcal c subroutine dxcal(f,dfdx) c This routine calculates the x derivative of f using centered c differences, and one-sided differences on the boundaries. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension f(ixmx,iymx),dfdx(ixmx,iymx) c c Calculate dx in meters dx = 1000.0*(xgk(2,1) - xgk(1,1)) tdx = 2.0*dx c c Calculate centered differences do j=1,ny do i=2,nx-1 dfdx(i,j) = (f(i+1,j) - f(i-1,j))/tdx enddo enddo c c Calculate one-sided differences do j=1,ny dfdx( 1,j) = (f( 2,j) - f( 1,j))/dx dfdx(nx,j) = (f(nx,j) - f(nx-1,j))/dx enddo c return end subroutine dxcal c subroutine dycal(f,dfdy) c This routine calculates the y derivative of f using centered c differences, and one-sided differences on the boundaries. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension f(ixmx,iymx),dfdy(ixmx,iymx) c c Calculate dy in meters dy = 1000.0*(ygk(1,2) - ygk(1,1)) tdy = 2.0*dy c c Calculate centered differences do i=1,nx do j=2,ny-1 dfdy(i,j) = (f(i,j+1) - f(i,j-1))/tdy enddo enddo c c Calculate one-sided differences do i=1,nx dfdy(i, 1) = (f(i, 2) - f(i, 1))/dy dfdy(i,ny) = (f(i,ny) - f(i,ny-1))/dy enddo c return end subroutine dycal c subroutine gshear(uvp,vvp,wvp,nvp,shrg) c This routine calculates a generalized shear parameter, which c includes contributes from nvp pressure levels c dimension uvp(nvp),vvp(nvp),wvp(nvp) c c Calculate the vertically averaged wind ubar = 0.0 vbar = 0.0 do n=1,nvp ubar = ubar + wvp(n)*uvp(n) vbar = vbar + wvp(n)*vvp(n) enddo c c Calculate generarlized shear shrg = 0.0 do n=1,nvp shrg = shrg + wvp(n)*sqrt( (uvp(n)-ubar)**2 + + (vvp(n)-vbar)**2 ) enddo c c Add factor to make generalized shear equal to two level shear c for linear shear profiles shrg = 4.0*shrg c c do n=1,nvp c write(6,100) n,uvp(n),vvp(n),wvp(n),ubar,vbar,shrg c 100 format(1x,'n,u,v,w,ub,vb,sg: ',i2,6(f6.2,1x)) c enddo c return end subroutine gshear c subroutine tgrad(slat,ptop,utop,vtop, + pbot,ubot,vbot,tgrd,tadv,ierr) c This routine calculates the magnitude of the temperature c gradient (tgrd in deg C per m) and temperature advection c (tadv in deg C per sec) between two pressure levels c using the geostrophic thermal wind equation. c c Input: slat (deg N), p,u,v (hPa, m/s, m/s) at the top and bottom c of the layer. c ierr = 0 c Check input if (slat .lt. -90.0 .or. slat .gt. 90.0) ierr=1 if (pbot .le. 0.0 .or. pbot .gt. 1100.0) ierr=1 if (ptop .le. 0.0 .or. ptop .gt. 1100.0) ierr=1 if (utop .lt. -200.0 .or. utop .gt. 200.0) ierr=1 if (vtop .lt. -200.0 .or. vtop .gt. 200.0) ierr=1 if (ubot .lt. -200.0 .or. ubot .gt. 200.0) ierr=1 if (vbot .lt. -200.0 .or. vbot .gt. 200.0) ierr=1 c if (ierr .ne. 0) return c f = 2.0*7.292e-5*sin(slat*0.017453) r = 287.0 plog = alog(pbot/ptop) cfac = f/(r*plog) ubar = 0.5*(utop+ubot) vbar = 0.5*(vtop+vbot) c dtdx = -cfac*(vbot-vtop) dtdy = cfac*(ubot-utop) tgrd = sqrt(dtdx**2 + dtdy**2) tadv = -ubar*dtdx - vbar*dtdy c return end subroutine tgrad c subroutine xshear(uvp,vvp,pvp,nvp,i1,i2,shrx) c This routine calculates the maximum shear in the layer c pvp(i1) to pvp(i2). The max shear magnitude is normalized c to the 850 to 200 mb layer. c dimension uvp(nvp),vvp(nvp),pvp(nvp) c c Check input if (i2 .gt. i1) then n1 = i1 n2 = i2 elseif (i2 .lt. i1) then n1 = i2 n2 = i1 else shrx = 999.9 return endif c if (n1 .gt. nvp .or. n2 .gt. nvp) then shrx = 999.9 return endif c if (n1 .lt. 1 .or. n2 .lt. 1) then shrx = 999.9 return endif c c Calculate the vertically averaged wind shrx = -999.9 dpn = (850.0-200.0) do n=n1,n2-1 du = (uvp(n+1)-uvp(n)) dv = (vvp(n+1)-vvp(n)) dp = (pvp(n+1)-pvp(n)) shrt = sqrt(du*du + dv*dv)*dpn/dp c if (shrt .gt. shrx) shrx = shrt enddo c return end subroutine xshear c subroutine vshear(ps1,ps2,rlon,rlat,radkt,minpts,shear,sheard) c This routine calculates the magnitude (shear) and direction c (sheard) of the vertical shear between pressure c levels ps1,ps2 (mb). The horizontal velocities are first averaged c over a circle of radius radkt (km), centered at rlon,rlat. The c shear is set to 999.9 if there are less than minpts included c in the horizontal average. The direction of the shear is in c degrees clockwise from north. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,prima common /primv/ u,v,z,t,r dimension u(ixmx,iymx,ipmx),v(ixmx,iymx,ipmx),z(ixmx,iymx,ipmx) dimension t(ixmx,iymx,ipmx),r(ixmx,iymx,ipmx) c C.DUP lsdiag,index common /indexa/ indexp,indexv common /indexn/ indexvn parameter (nvar=5) dimension indexp(ipmx),indexv(nvar) character *1 indexvn(nvar) C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,plevs common /plevs/ plev,plevr dimension plev(ipmx),plevr(ipmx) C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Local variables dimension rad(ixmx,iymx) c c Find indicies of required pressure levels ips1 = 0 ips2 = 0 do k=1,ipmx if (plevr(k) .eq. ps1) ips1=k if (plevr(k) .eq. ps2) ips2=k enddo c if (ips1 .eq. 0 .or. ips2 .eq. 0) then write(lulog,900) ps1,ps2 900 format(/,' Pressure level not found in routine vshear', + ' ps1,ps2: ' ,f6.1,1x,f6.1) stop endif c if (indexp(ips1) .eq. 0 .or. indexp(ips2) .eq. 0) then write(lulog,902) ips1,ips2 902 format(/,' Pressure level not available on data file', + ' ps1,ps2: ',f6.1,1x,f6.1) stop endif c if (indexv(1) .eq. 0 .or. indexv(2) .eq. 0) then write(lulog,904) 904 format(/,' Wind data not available for shear calculation.') stop endif c c Calculate the radius of each grid point call xyracal(rlon,rlat) c c Apply horizontal average to wind components ubar1=0.0 ubar2=0.0 vbar1=0.0 vbar2=0.0 nbar =0 do j=1,nlat1 do i=1,nlon1 if (rgk(i,j) .le. radkt) then ubar1 = ubar1 + u(i,j,ips1) ubar2 = ubar2 + u(i,j,ips2) vbar1 = vbar1 + v(i,j,ips1) vbar2 = vbar2 + v(i,j,ips2) nbar = nbar + 1 endif enddo enddo c if (nbar .lt. minpts) then write(lulog,800) nbar,minpts 800 format(/,' Not enough data points for shear calculation.', + ' nbar,minpts: ',i4,1x,i4) shear = 999.9 return endif c rnbar = float(nbar) ubar1 = ubar1/rnbar ubar2 = ubar2/rnbar vbar1 = vbar1/rnbar vbar2 = vbar2/rnbar c shearx = ubar2 - ubar1 sheary = vbar2 - vbar1 call ctorh(shearx,sheary,shear,sheard) c c Convert shear to knots shear = 1.944*shear c return END subroutine vshear c subroutine xyracal(olon,olat) c This routine calculates the x and y coordinates and radius c in km and azimuth in deg relative to the origin olon,olat c (olon is deg W negative, olat is deg N positive) c at all of the Cartesian grid points. c Then, the radii and azimuths of a polar coordinate system c centered at olon,olat are calculated. Finally, the array pmask c is calculated which contains 0.0 or 1.0 if the polar grid point c is outside or inside of the Cartesian grid. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c ************************* c c Save lat,lon of origin olons = olon olats = olat colat = cos(dtr*olat) c c Determine number of x,y grid points nx = nlon1 ny = nlat1 c c Calculate x,y,r and azimuth of the Cartesian grid do j=1,ny do i=1,nx tlon = rlond(i) tlat = rlatd(j) call lltoxy(olon,olat,colat,tlon,tlat,xtem,ytem) xgk(i,j) = xtem ygk(i,j) = ytem enddo enddo c do j=1,ny do i=1,nx rgk(i,j) = sqrt(xgk(i,j)**2 + ygk(i,j)**2) if (rgk(i,j) .gt. 0.0) then agd(i,j) = ( acos(xgk(i,j)/rgk(i,j)) )/dtr else agd(i,j) = 0.0 endif enddo enddo c xmax = xgk(nx,1) xmin = xgk(1,1) ymax = ygk(1,ny) ymin = ygk(1,1) c c Calculate azimuthal coordinates ntheta = iamx dtheta = 360.0/float(ntheta) do j=1,ntheta thetad(j) = dtheta*float(j-1) enddo c c Specify radial interval in km drad = 100.0 c c Find maximum radius value on Cartesian grid radmax = -1000.0 do j=1,ny do i=1,nx if (rgk(i,j) .gt. radmax) radmax=rgk(i,j) enddo enddo c c Calculate radial coordinates nrad = 2 + ifix(radmax/drad) if (nrad .gt. irmx) then write(lulog,900) nrad,irmx 900 format(/,' nrad too big in routine xyracal, increase irmx.', + /,' nrad,irmx: ',i5,1x,i5) stop endif c do i=1,nrad radk(i) = drad*float(i-1) enddo c c Calculate x,y values at polar grid points do j=1,ntheta do i=1,nrad xra(i,j) = radk(i)*cos(dtr*thetad(j)) yra(i,j) = radk(i)*sin(dtr*thetad(j)) enddo enddo c xmaxm1 = xgk(nx-1,1) ymaxm1 = ygk(1,ny-1) c Calculate polar coordinate mask do j=1,ntheta do i=1,nrad xtem = xra(i,j) ytem = yra(i,j) pmask(i,j) = 1.0 if (xtem .ge. xmaxm1 .or. xtem .le. xmin) pmask(i,j)=0.0 if (ytem .ge. ymaxm1 .or. ytem .le. ymin) pmask(i,j)=0.0 enddo enddo c return END subroutine xyracal c subroutine lltoxy(olon,olat,colat,tlon,tlat,xtem,ytem) c This routine calculates the x,y coordinates (xtem,ytem) of the c point with longitude and latitude tlon,tlat relative to a c coordinate system with origin at olon,olat. In the current c version, a simple tangent plane approximation is used for the c transformation. c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c xtem = dtk*(tlon-olon)*colat ytem = dtk*(tlat-olat) c return END subroutine lltoxy c subroutine xytoll(olon,olat,colat,xtem,ytem,tlon,tlat) c This routine calculates the longitude and latitude tlon,tlat c of the point with x,y coordinates xtem,tem c relative to a coordinate system with origin at c olon,olat. In the current version, a simple tangent plane c approximation is used for the transformation. c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c tlon = olon + xtem/(dtk*colat) tlat = olat + ytem/(dtk) c return END subroutine xytoll c subroutine uvtort(u,v,ur,vt) c This routine calculates the radial (ur) and tangential (vt) c winds at the Cartesian grid points from the Cartesian c velocity components u,v. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Other variables dimension u(ixmx,iymx),v(ixmx,iymx),ur(ixmx,iymx),vt(ixmx,iymx) c c ******************** c do j=1,ny do i=1,nx radtem = amax1(0.1,rgk(i,j)) ctem = xgk(i,j)/radtem stem = ygk(i,j)/radtem ur(i,j) = ctem*u(i,j) + stem*v(i,j) vt(i,j) = -stem*u(i,j) + ctem*v(i,j) enddo enddo c return END subroutine uvtort c subroutine rttouv(ur,vt,u,v) c This routine calculates the Cartesian velocity components u,v c from the radial (ur) and tangential (vt) winds at the c Cartesian grid points. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c c Other variables dimension u(ixmx,iymx),v(ixmx,iymx),ur(ixmx,iymx),vt(ixmx,iymx) c c ******************** c do j=1,ny do i=1,nx radtem = amax1(0.1,rgk(i,j)) ctem = xgk(i,j)/radtem stem = ygk(i,j)/radtem u(i,j) = ctem*ur(i,j) - stem*vt(i,j) v(i,j) = stem*ur(i,j) + ctem*vt(i,j) enddo enddo c return END subroutine rttouv c subroutine varget(vtype,press,fxy,ifatal,iexist) c This routine gets the variable vtype at pressure level press c and puts it in the array fxy. c c If ifatal=1 then execution is terminated if the requested c field is not available. c If ifatal=0 then iexist is set to 1 (0) if the requested c field exists (does not exist). Execution is not c terminated if ifatal=0. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,prima common /primv/ u,v,z,t,r dimension u(ixmx,iymx,ipmx),v(ixmx,iymx,ipmx),z(ixmx,iymx,ipmx) dimension t(ixmx,iymx,ipmx),r(ixmx,iymx,ipmx) C.DUP lsdiag,index common /indexa/ indexp,indexv common /indexn/ indexvn parameter (nvar=5) dimension indexp(ipmx),indexv(nvar) character *1 indexvn(nvar) C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,plevs common /plevs/ plev,plevr dimension plev(ipmx),plevr(ipmx) C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c c Other variables character *1 vtype dimension fxy(ixmx,iymx) c c ******************** c c Search for variable type do n=1,nvar if (indexvn(n) .eq. vtype) then idv = n exit endif if (n .eq. nvar) then write(lulog,900) vtype 900 format(/,' Unrecognized variable type: ',a1, + ' in routine varget.') stop endif enddo c c Check to make sure variable is available if (indexv(idv) .eq. 0) then if (ifatal .eq. 1) then write(lulog,902) vtype 902 format(/,' Variable ',a1, + ' not available on data input file.', + /,' Processing halted in routine varget.') stop else write(lulog,802) vtype 802 format(/,' Variable ',a1, + ' not available on data input file.') iexist=0 return endif endif c c Search for pressure level do k=1,ipmx if (plevr(k) .eq. press) then idp = k if (indexp(k) .ne. 0) exit endif if (k .eq. ipmx) then write(lulog,904) press 904 format(/,' Pressure level: ',e11.4, + ' not found in routine varget.') if (ifatal .eq. 1) then stop else iexist=0 return endif endif enddo c iexist=1 c c Extract variable if (idv .eq. 1) then do j=1,nlat1 do i=1,nlon1 fxy(i,j) = u(i,j,idp) enddo enddo elseif (idv .eq. 2) then do j=1,nlat1 do i=1,nlon1 fxy(i,j) = v(i,j,idp) enddo enddo elseif (idv .eq. 3) then do j=1,nlat1 do i=1,nlon1 fxy(i,j) = z(i,j,idp) enddo enddo elseif (idv .eq. 4) then do j=1,nlat1 do i=1,nlon1 fxy(i,j) = t(i,j,idp) enddo enddo elseif (idv .eq. 5) then do j=1,nlat1 do i=1,nlon1 fxy(i,j) = r(i,j,idp) enddo enddo else write(lulog,908) idv 908 format(/,' Illegal variable number ',i2,' in routine varget.') stop endif c return END subroutine varget c subroutine ctop(fxy,fra) c This routine linearly interpolates the function fxy on c the Cartesian grid to the function fra on the polar grid. c fra is set to zero for points on the polar grid outside c of the Cartesian grid. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c c Other variables dimension fxy(ixmx,iymx),fra(irmx,iamx) c nlonck = nlon1-1 nlatck = nlat1-1 c do j=1,ntheta do i=1,nrad if (pmask(i,j) .lt. 0.5) then fra(i,j) = 0.0 else xtem = xra(i,j) ytem = yra(i,j) call xytoll(olons,olats,colat,xtem,ytem,tlon,tlat) c c Calculate x,y indices of the point closet to, and to the c lower left of the current polar coordinate grid point. idx = 1 + ifix( (tlon-rlond(1))/dlon ) idy = 1 + ifix( (tlat-rlatd(1))/dlat ) c c Check bounds of x,y indices if (idx .lt. 1 .or. idx .gt. nlonck) then write(lulog,900) idx 900 format(/,' idx=',i4,' out of bounds in routine ctop.') stop endif c if (idy .lt. 1 .or. idy .gt. nlatck) then write(lulog,902) idy 902 format(/,' idy=',i4,' out of bounds in routine ctop.') stop endif c x11 = xgk(idx ,idy ) x21 = xgk(idx+1,idy ) x12 = xgk(idx ,idy+1) x22 = xgk(idx+1,idy+1) c y11 = ygk(idx ,idy ) y21 = ygk(idx+1,idy ) y12 = ygk(idx ,idy+1) y22 = ygk(idx+1,idy+1) c f11 = fxy(idx ,idy ) f21 = fxy(idx+1,idy ) f12 = fxy(idx ,idy+1) f22 = fxy(idx+1,idy+1) c c Normalize x,y of polar grid point xnorm = (xtem-x11)/(x21-x11) ynorm = (ytem-y11)/(y12-y11) c a = f11 b = f21 - f11 c = f12 - f11 d = f11 + f22 - (f21 + f12) c fra(i,j) = a + b*xnorm + c*ynorm + d*xnorm*ynorm endif enddo enddo c return END subroutine ctop c subroutine aavg(fxy,radinn,radout,minpts,ifatal,ierr,fxybar) c This routine horizontally averages fxy over an annular c region with r between radinn and radout. The center of the circle c is determined by the last call to routine xyracal. c minpts is the minimum number of points allowed in the c circular area. If ifatal=1, execution stops if there is c not enough points for the calculation. If ifatal=0, execution c continues, but ierr is set to 1 if there are not enough points. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c dimension fxy(ixmx,iymx) c fxybar= 0.0 nbar = 0 do j=1,ny do i=1,nx if (rgk(i,j) .le. radout .and. rgk(i,j) .ge. radinn) then fxybar = fxybar + fxy(i,j) nbar = nbar + 1 endif enddo enddo c if (nbar .lt. minpts) then write(lulog,900) nbar,minpts 900 format(/,' Not enough data points for aavg calculation.', + ' nbar,minpts: ',i4,1x,i4) if (ifatal .eq. 1) then stop else ierr=1 return endif endif c fxybar = fxybar/float(nbar) ierr=0 c return END subroutine aavg c subroutine havg(fxy,radtem,minpts,ifatal,ierr,fxybar) c This routine horizontally averages fxy over a circular c region with radius radtem. The center of the circle c is determined by the last call to routine xyracal. c minpts is the minimum number of points allowed in the c circular area. If ifatal=1, execution stops if there is c not enough points for the calculation. If ifatal=0, execution c continues, but ierr is set to 1 if there are not enough points. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c dimension fxy(ixmx,iymx) c fxybar= 0.0 nbar = 0 do j=1,ny do i=1,nx if (rgk(i,j) .le. radtem) then fxybar = fxybar + fxy(i,j) nbar = nbar + 1 endif enddo enddo c if (nbar .lt. minpts) then write(lulog,900) nbar,minpts 900 format(/,' Not enough data points for havg calculation.', + ' nbar,minpts: ',i4,1x,i4) if (ifatal .eq. 1) then stop else ierr=1 return endif endif c fxybar = fxybar/float(nbar) ierr=0 c return END subroutine havg c subroutine hcavg(fxy,xtem,ytem,minpts,ifatal,ierr,fxybar) c This routine horizontally averages fxy over a rectangular c region from -xtem,ytem to +xtem,ytem. The center of the rectangle c is determined by the last call to routine xyracal. c minpts is the minimum number of points allowed in the c rectangular area. If ifatal=1, execution stops if there is c not enough points for the calculation. If ifatal=0, execution c continues, but ierr is set to 1 if there are not enough points. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c dimension fxy(ixmx,iymx) c fxybar= 0.0 nbar = 0 do j=1,ny do i=1,nx if (xgk(i,j) .le. xtem .and. xgk(i,j) .ge. -xtem .and. + ygk(i,j) .le. ytem .and. ygk(i,j) .ge. -ytem ) then fxybar = fxybar + fxy(i,j) nbar = nbar + 1 endif enddo enddo c if (nbar .lt. minpts) then write(lulog,900) nbar,minpts 900 format(/,' Not enough data points for hcavg calculation.', + ' nbar,minpts: ',i4,1x,i4) if (ifatal .eq. 1) then stop else ierr=1 return endif endif c fxybar = fxybar/float(nbar) ierr=0 c return END subroutine hcavg c subroutine azavg(fra,fr,iorz) c This routine azimuthally averages fra on the polar grid. c If iorz=0 then fr is set to zero at the origin. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension fra(irmx,iamx),fr(irmx) c do i=1,nrad jcount = 0 ftemp = 0.0 do j=1,ntheta if (pmask(i,j) .gt. 0.5) then ftemp = ftemp + fra(i,j) jcount = jcount + 1 endif enddo if (jcount .gt. 0) then fr(i) = ftemp/float(jcount) else fr(i) = 0.0 endif enddo c if (iorz .eq. 0) fr(1)=0.0 c return END subroutine azavg c subroutine patoc(fr,fxy) c This routine linearly interpolates the azimuthally averaged c function fr on the polar grid to fxy on the Cartesian grid. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension fr(irmx),fxy(ixmx,iymx) c do j=1,ny do i=1,nx radtem = rgk(i,j) c c Find index of radial point closest to, but c less than current Cartesian grid point idr = 1 + ifix(radtem/drad) c c Scale the radial distance rnorm = (radtem-radk(idr))/drad c c Perform linear interpolation fxy(i,j) = fr(idr) + (fr(idr+1)-fr(idr))*rnorm enddo enddo c return END subroutine patoc c subroutine ddrcal(fr,dfdr) c This routine calculates the radial derivative of fr c using centered derivatives. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension fr(irmx),dfdr(irmx) c c ******************** c drm = 1000.0*drad drm2 = 2.0*drm c c Interior points do i=3,nrad-1 dfdr(i) = (fr(i+1)-fr(i-1))/drm2 enddo c c End points dfdr(1) = 0.0 dfdr(2) = (fr( 3)-fr( 2))/drm dfdr(nrad) = (fr(nrad)-fr(nrad-1))/drm c return END subroutine ddrcal c subroutine vsound(rlon,rlat,sradk) c This routine horizontally averages variables at c specified pressure levels and puts the results pressure c arrays. The average is over a circular area of radius sradk c (km) centered at rlon,rlat. c c If a variable is not available on the data input c file, the value is set to 999.9 c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,plevh parameter (nplevh=10) common /plevhs/ up,vp,zp,tp,rp,tbarp,plevh dimension up(nplevh),vp(nplevh),zp(nplevh) dimension tp(nplevh),rp(nplevh) dimension tbarp(nplevh) dimension plevh(nplevh) C.DUP lsdiag,plevhd data plevh /1000.0, 850.0, 700.0, 500.0, 400.0, + 300.0, 250.0, 200.0, 150.0, 100.0/ data tbarp /287.43,278.68,268.57,251.92,241.44, + 228.58,220.79,216.65,216.65,216.65/ C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c dimension fxy(ixmx,iymx) c c Specify minimum number of points for the horizontal average minpts=1 c c Set itsub=1 to subtract standard atmosphere temperatures c from observed temperatures or c itsub=2 to convert from Kelvin to Celsius. c itsub=2 c c Initialize grid variables call xyracal(rlon,rlat) c c Start pressure level loop do k=1,nplevh ptem = plevh(k) c call varget('U',ptem,fxy,0,iexist) if (iexist .eq. 1) then call havg(fxy,sradk,minpts,1,ierr,fbar) up(k) = 1.944*fbar else up(k) = 999.9 endif c call varget('V',ptem,fxy,0,iexist) if (iexist .eq. 1) then call havg(fxy,sradk,minpts,1,ierr,fbar) vp(k) = 1.944*fbar else vp(k) = 999.9 endif c call varget('Z',ptem,fxy,0,iexist) if (iexist .eq. 1) then call havg(fxy,sradk,minpts,1,ierr,fbar) zp(k) = fbar else zp(k) = 999.9 endif c call varget('T',ptem,fxy,0,iexist) if (iexist .eq. 1) then call havg(fxy,sradk,minpts,1,ierr,fbar) if (itsub .eq. 1) then tp(k) = fbar - tbarp(k) elseif (itsub .eq. 2) then tp(k) = fbar - 273.15 else tp(k) = fbar endif else tp(k) = 999.9 endif c call varget('R',ptem,fxy,0,iexist) if (iexist .eq. 1 .and. ptem .gt. 275.0) then call havg(fxy,sradk,minpts,1,ierr,fbar) rp(k) = fbar else rp(k) = 999.9 endif enddo c write(lulog,300) 300 format(/,' P U V Z T RH') do k=1,nplevh write(lulog,302) plevh(k),up(k),vp(k),zp(k),tp(k),rp(k) 302 format(1x,6(f6.1,1x)) enddo c return END subroutine vsound c subroutine vcsound(rlon,rlat,xradk,yradk) c This routine horizontally averages variables at c specified pressure levels and puts the results pressure c arrays. The average is over a rectangulare area c (+/-xradk,+/-yradk) (km) centered at rlon,rlat. c c If a variable is not available on the data input c file, the value is set to 999.9 c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,plevh parameter (nplevh=10) common /plevhs/ up,vp,zp,tp,rp,tbarp,plevh dimension up(nplevh),vp(nplevh),zp(nplevh) dimension tp(nplevh),rp(nplevh) dimension tbarp(nplevh) dimension plevh(nplevh) C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c dimension fxy(ixmx,iymx) c c Specify minimum number of points for the horizontal average minpts=1 c c Set itsub=1 to subtract standard atmosphere temperatures c from observed temperatures or c itsub=2 to convert from Kelvin to Celsius. c itsub=2 c c Initialize grid variables call xyracal(rlon,rlat) c c Start pressure level loop do k=1,nplevh ptem = plevh(k) c call varget('U',ptem,fxy,0,iexist) if (iexist .eq. 1) then call hcavg(fxy,xradk,yradk,minpts,1,ierr,fbar) up(k) = 1.944*fbar else up(k) = 999.9 endif c call varget('V',ptem,fxy,0,iexist) if (iexist .eq. 1) then call hcavg(fxy,xradk,yradk,minpts,1,ierr,fbar) vp(k) = 1.944*fbar else vp(k) = 999.9 endif c call varget('Z',ptem,fxy,0,iexist) if (iexist .eq. 1) then call hcavg(fxy,xradk,yradk,minpts,1,ierr,fbar) zp(k) = fbar else zp(k) = 999.9 endif c call varget('T',ptem,fxy,0,iexist) if (iexist .eq. 1) then call hcavg(fxy,xradk,yradk,minpts,1,ierr,fbar) if (itsub .eq. 1) then tp(k) = fbar - tbarp(k) elseif (itsub .eq. 2) then tp(k) = fbar - 273.15 else tp(k) = fbar endif else tp(k) = 999.9 endif c call varget('R',ptem,fxy,0,iexist) if (iexist .eq. 1 .and. ptem .gt. 275.0) then call hcavg(fxy,xradk,yradk,minpts,1,ierr,fbar) rp(k) = fbar else rp(k) = 999.9 endif enddo c write(lulog,300) 300 format(/,' P U V Z T RH') do k=1,nplevh write(lulog,302) plevh(k),up(k),vp(k),zp(k),tp(k),rp(k) 302 format(1x,6(f6.1,1x)) enddo c return END subroutine vcsound c subroutine fcopy(fxy,fxyc) c This routine copies the 2-D array fxy to fxyc c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c dimension fxy(ixmx,iymx),fxyc(ixmx,iymx) c do j=1,ny do i=1,nx fxyc(i,j) = fxy(i,j) enddo enddo c return end subroutine fcopy c c ** cfindxy code marker (Start) ** c c This group of routines is for finding the circulation center c and removing the storm influence from various fields. c c Routines include c lfilter c cfindxy c lltxycf c xytllcf c lintcf c azavgcf c rdavgcf c rdavgcf c vremove c ctor c rtoc c subroutine lfilter(fxy,sxc,syc,frad) c This routine applies a Laplacian filter to all points within c a distance frad (km) to the point sxc, syc c c ++ common start in lsplot code ++ c common /xyrai/ nx,ny c common /xyrar/ xgk,ygk c c parameter (ixmx=26,iymx=21) c dimension xgk(ixmx,iymx),ygk(ixmx,iymx) c ++ common end in lsplot code ++ c c ++ common start in lsdiag code ++ C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c ++ common end in lsdiag code ++ c dimension fxy(ixmx,iymx) dimension mask(ixmx,iymx) c do j=1,ny do i=1,nx mask(i,j) = 0 enddo enddo c c Find maximum search indices dx = xgk(2,1)-xgk(1,1) dy = ygk(1,2)-ygk(1,1) kx = 1 + ifix(frad/dx) ky = 1 + ifix(frad/dx) c kxo = 1 + ifix( 0.5 + (sxc-xgk(1,1))/dx ) kyo = 1 + ifix( 0.5 + (syc-ygk(1,1))/dy ) c i1 = kxo - kx i2 = kxo + kx c j1 = kyo - ky j2 = kyo + ky c if (i1 .lt. 2) i1 = 2 if (i2 .gt. nx-1) i2 = nx-1 c if (j1 .lt. 2) j1 = 2 if (j2 .gt. ny-1) j2 = ny-1 c do j=j1,j2 do i=i1,i2 rad = sqrt( (xgk(i,j)-sxc)**2 + (ygk(i,j)-syc)**2 ) if (rad .le. frad) mask(i,j) = 1 enddo enddo c c do j=ny,1,-1 c write(6,800) (mask(i,j),i=1,nx) c 800 format(26(i2)) c enddo c c Perform relaxation for Laplacian filter solution nit=20 c r2 = (dy/dx)**2 c yfac = 0.5/(r2+1.0) xfac = r2*yfac c do n=1,nit do j=j1,j2 do i=i1,i2 if (mask(i,j) .eq. 1) then fxy(i,j) = xfac*(fxy(i+1,j)+fxy(i-1,j)) + + yfac*(fxy(i,j+1)+fxy(i,j-1)) endif enddo enddo enddo c return end subroutine lfilter c subroutine cfindxy(uxy,vxy,rcf,thetacf, + mrcf,mtcf,nrcf,ntcf, + radtwa,x0,y0,xcen,ycen,twa,rwa) c c This routine finds the coordinates of the center location c (xcen,ycen) in km that maximizes the mean tangential wind averaged c from 0 to radtwm (km), given the first guess coordinates (x0,y0) c and the horizontal wind components uxy,vxy. The average tangential c wind relative to the center location (twa) is also calculated. c c Input: uxy,vxy - horizontal wind components as a function of x,y c rcf(0:mrcf) - radius array for storm-rel coord's c thetacf(0:mtcf) - azimuth array for storm-rel coord's c mrcf,mtcf - max dimensions of r,theta arrays c nrcf,ntcf - working dimensions of r,theta arrays c radtwa - Maximum radius for radially averaged c tangential wind (km) c x0,y0 - first guess coordinates of center c c Output: xcen,ycen - coordinates of storm center (km) c twa,rwa - azimutally and radially averaged tangential c and radial winds relative to (xcen,ycen) c c Note: The x,y coordinates of uxy,vxy are passed through common c c ++ common start in lsplot code ++ c common /xyrai/ nx,ny c common /xyrar/ xgk,ygk c c parameter (ixmx=26,iymx=21) c dimension xgk(ixmx,iymx),ygk(ixmx,iymx) c ++ common end in lsplot code ++ c c ++ common start in lsdiag code ++ C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c ++ common end in lsdiag code ++ c c ++ Passed variables dimension uxy(ixmx,iymx),vxy(ixmx,iymx) dimension rcf(0:mrcf),thetacf(0:mtcf) c c ++ Temporary cylindrical grid variables parameter (irxcf=100,itxcf=32) dimension xrt(0:irxcf,0:itxcf), yrt(0:irxcf,0:itxcf) dimension ucrt(0:irxcf,0:itxcf),vcrt(0:irxcf,0:itxcf) dimension urrt(0:irxcf,0:itxcf),vtrt(0:irxcf,0:itxcf) dimension urr(0:irxcf),vtr(0:irxcf) dimension stheta(0:itxcf),ctheta(0:itxcf) c c Specify number of interations niter=8 c c Specify initial dx,dy (km) for the search dxy = 100.0 c c Specify the fraction to reduce dx,dy after each iteration dxyrf = 0.6 c c Specify the number of points on either side of the c first guess lat/lon for the search nsearch = 2 c c Calculate needed geometric info dtr = 3.14159/180.0 dx = xgk(2,1)-xgk(1,1) dy = ygk(1,2)-ygk(1,1) x1 = xgk(1,1) y1 = ygk(1,1) c do j=0,ntcf stheta(j) = sin(dtr*thetacf(j)) ctheta(j) = cos(dtr*thetacf(j)) enddo c c First guess center xcen = x0 ycen = y0 xcent= x0 ycent= y0 twa = -1.0e10 c c Determine if this is a southern hemisphere case if (olats .lt. 0.0) then sgnlat = -1.0 else sgnlat = 1.0 endif c c ++ Start iteration findloop: do n=1,niter c do jj=-nsearch,nsearch do ii=-nsearch,nsearch c c Calculate x,y coordinates of center-relative c r,theta grid points c xcenij = xcen + dxy*float(ii) ycenij = ycen + dxy*float(jj) c do j=0,ntcf do i=0,nrcf xrt(i,j) = xcenij + rcf(i)*ctheta(j) yrt(i,j) = ycenij + rcf(i)*stheta(j) enddo enddo c c Interpolate u,v from x,y to r,t do j=0,ntcf do i=0,nrcf call lintcf(uxy,x1,dx,y1,dy,ixmx,iymx,nx,ny, + xrt(i,j),yrt(i,j),ucrt(i,j)) call lintcf(vxy,x1,dx,y1,dy,ixmx,iymx,nx,ny, + xrt(i,j),yrt(i,j),vcrt(i,j)) enddo enddo c c Convert Cartesian u,v to radial and tangential winds do j=0,ntcf do i=0,nrcf urrt(i,j) = ucrt(i,j)*ctheta(j) + + vcrt(i,j)*stheta(j) vtrt(i,j) = vcrt(i,j)*ctheta(j) - + ucrt(i,j)*stheta(j) enddo enddo c c ++ Debug code c if (ii .eq. 0 .and. jj .eq. 0) then c do jk=10,6,-1 c if (jk .eq. 10) write(6,849) (ik,ik=8,12) c 849 format(/,'uxy',/,4x,20(1x,i4,1x)) c c write(6,850) jk,(uxy(ik,jk),ik=8,12) c 850 format(1x,i2,1x,20(f5.1,1x)) c enddo c c do jk=10,6,-1 c if (jk .eq. 10) write(6,851) (ik,ik=8,12) c 851 format(/,'vxy',/,4x,20(1x,i4,1x)) c c write(6,850) jk,(vxy(ik,jk),ik=8,12) c enddo c c do jk=0,ntcf c if (jk .eq. 0) write(6,852) (ik,ik=0,10) c 852 format(/,'ucrt',/,4x,20(1x,i4,1x)) c c write(6,850) jk,(ucrt(ik,jk),ik=0,10) c enddo c c do jk=0,ntcf c if (jk .eq. 0) write(6,853) (ik,ik=0,10) c 853 format(/,'vcrt',/,4x,20(1x,i4,1x)) c c write(6,850) jk,(vcrt(ik,jk),ik=0,10) c enddo c endif c c Azimuthally average ur,vt call azavgcf(urrt,urr,mrcf,mtcf,nrcf,ntcf) call azavgcf(vtrt,vtr,mrcf,mtcf,nrcf,ntcf) c c do k=0,nrcf c write(6,704) rcf(k),urr(k),vtr(k) c 704 format('r,ur,vt: ',3(f6.1,1x)) c enddo c c Radially average ur,vt call rdavgcf(urr,rcf,radtwa,ubar,mrcf,nrcf) call rdavgcf(vtr,rcf,radtwa,vbar,mrcf,nrcf) vbar = vbar*sgnlat c if (vbar .gt. twa) then xcent = xcenij ycent = ycenij twa = vbar rwa = ubar endif enddo enddo c xcen = xcent ycen = ycent c c write(6,605) n,dxy,xcen,ycen,twa,rwa c 605 format('n,dxy: ',i2,1x,f6.2,' xc,yc: ',f6.1,1x,f6.1, c + ' twa,ura: ',f6.2,1x,f6.2) c dxy = dxy*dxyrf c enddo findloop c c do k=0,nrcf c write(6,704) rcf(k),urr(k),vtr(k) c 704 format('r,ur,vt: ',3(f6.1,1x)) c enddo c return end subroutine cfindxy c subroutine lltxycf(olat,olon,colat,dtk,rlat,rlon,x,y) c This routine calculates the x,y coordinates of the c point (rlon,rlat) where the origin is at olon,olat and c colat=cos(olat). All lat/lons are in deg N/E and x,y c are in km, and dtk is the deg lat to km conversion factor. c All parameters except x,y are input parameters. c x = dtk*(rlon-olon)*colat y = dtk*(rlat-olat) c return end subroutine lltxycf c subroutine xytllcf(olat,olon,colat,dtk,x,y,rlat,rlon) c This routine calculates the lat,lon coordinates of the c point (x,y) where the origin is at olon,olat and c colat=cos(olat). All lat/lons are in deg N/E and x,y c are in km, and dtk is the deg lat to km conversion factor. c All parameters except rlat,rlon are input parameters. c rlon = olon + x/(dtk*colat) rlat = olat + y/(dtk ) c return end subroutine xytllcf c subroutine lintcf(fxy,x1,dx,y1,dy,mx,my,nx,ny,xi,yi,fxyii) c This routine bi-linearly interpolates fxy to the point c (xi,yi) to give fxyii. c c Input: fxy(mx,my) - array to be interpolated c dx,dy - spacing of evenly spaced x,y points c x1,y1 - x,y coorindates of lower-left point c mx,my - max dimensions of fxy c nx,ny - working dimensions of fxy c xi,yi - coordinates of point to be interpolated c c Output: fxyii - value of interpolated function at (xi,yi) c dimension fxy(mx,my) c c Find the indices of the lower-left point of the c c grid box containing the point (xi,yi) i0 = 1 + ifix( (xi-x1)/dx ) j0 = 1 + ifix( (yi-y1)/dy ) c c Check index bounds if (i0 .lt. 1) i0= 1 if (i0 .gt. nx-1) i0=nx-1 if (j0 .lt. 1) j0= 1 if (j0 .gt. ny-1) j0=ny-1 c i1 = i0+1 j1 = j0+1 c c Calculate normalized x,y distances xn = ( (xi-x1) - dx*float(i0-1) )/dx yn = ( (yi-y1) - dy*float(j0-1) )/dx c if (xn .lt. 0.0) xn = 0.0 if (xn .gt. 1.0) xn = 1.0 if (yn .lt. 0.0) yn = 0.0 if (yn .gt. 1.0) yn = 1.0 c c Calculate coefficients for interpolation function f00 = fxy(i0,j0) f10 = fxy(i1,j0) f01 = fxy(i0,j1) f11 = fxy(i1,j1) c a = f00 b = f10-f00 c = f01-f00 d = f00+f11-f10-f01 c fxyii = a + b*xn + c*yn + d*xn*yn c return end subroutine lintcf c subroutine azavgcf(frt,fr,mr,mt,nr,nt) c This routine performs and azimuthal average of frt to give ft dimension frt(0:mr,0:mt),fr(0:mr) c cf = 1.0/float(nt+1) c do i=0,nr fr(i) = frt(i,0) do j=1,nt fr(i) = fr(i) + frt(i,j) enddo c fr(i) = cf*fr(i) enddo c return end subroutine azavgcf c subroutine rdavgcf(fr,rr,radmax,fbar,mr,nr) c This routine radially averages fr from r=0 to r=radmax c dimension fr(0:mr),rr(0:mr) c count = 0.0 fbar = 0.0 do i=0,nr if (rr(i) .gt. radmax) exit count = count + 1.0 fbar = fbar + fr(i) enddo c if (count .ge. 1.0) then fbar = fbar/count else fbar = 0.0 endif c return end subroutine rdavgcf c subroutine vremove(uxy,vxy,zxy,txy,rcf,thetacf,xcen,ycen, + mrcf,mtcf,nrcf,ntcf,uraa,vtaa,zaa,taa,irflag) c c This routine calculates the azimuthally averaged radial c and tangential winds (uraa and vtaa) relative to the point c (xcen,ycen) from the Cartesian wind components uxy,vxy. c c The azimumthally averaged temperature (taa) is also calculated. c c If irflag=1, then the symmetric winds are removed from uxy,vxy c in regions of cyclonic rotation. c c ++ common start in lsplot code ++ c common /xyrai/ nx,ny c common /xyrar/ xgk,ygk c c parameter (ixmx=26,iymx=21) c dimension xgk(ixmx,iymx),ygk(ixmx,iymx) c ++ common end in lsplot code ++ c c ++ common start in lsdiag code ++ C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,xyra common /xyra/ xgk,ygk,rgk,agk,xmin,xmax,ymin,ymax, + olons,olats,colat dimension xgk(ixmx,iymx),ygk(ixmx,iymx) dimension rgk(ixmx,iymx),agd(ixmx,iymx) common /ixyra/ nx,ny common /raxy/ radk,thetad,xra,yra,pmask,drad,dtheta dimension radk(irmx),thetad(iamx) dimension xra(irmx,iamx),yra(irmx,iamx) dimension pmask(irmx,iamx) common /iraxy/ nrad,ntheta c ++ common end in lsdiag code ++ c c ++ Passed variables dimension uxy(ixmx,iymx),vxy(ixmx,iymx) dimension zxy(ixmx,iymx),txy(ixmx,iymx) dimension rcf(0:mrcf),thetacf(0:mtcf) dimension uraa(0:mrcf),vtaa(0:mrcf) dimension zaa(0:mrcf), taa(0:mrcf) c c ++ Temporary cylindrical grid variables parameter (irxcf=100,itxcf=32) dimension xrt(0:irxcf,0:itxcf), yrt(0:irxcf,0:itxcf) dimension ucrt(0:irxcf,0:itxcf),vcrt(0:irxcf,0:itxcf) dimension urrt(0:irxcf,0:itxcf),vtrt(0:irxcf,0:itxcf) dimension zrt(0:irxcf,0:itxcf), trt(0:irxcf,0:itxcf) dimension urr(0:irxcf),vtr(0:irxcf) dimension stheta(0:itxcf),ctheta(0:itxcf) c c Specify minimum radius (km) for vortex removal rminvr = 301.0 c c Calculate needed geometric info dtr = 3.14159/180.0 dx = xgk(2,1)-xgk(1,1) dy = ygk(1,2)-ygk(1,1) x1 = xgk(1,1) y1 = ygk(1,1) dr = rcf(2)-rcf(1) c c Determine if this is a southern hemisphere case if (olats .lt. 0.0) then sgnlat = -1.0 else sgnlat = 1.0 endif c do j=0,ntcf stheta(j) = sin(dtr*thetacf(j)) ctheta(j) = cos(dtr*thetacf(j)) enddo c c ++ Find symmetric average ++ c c Find x,y coordinates of storm-centered radial grid points do j=0,ntcf do i=0,nrcf xrt(i,j) = xcen + rcf(i)*ctheta(j) yrt(i,j) = ycen + rcf(i)*stheta(j) enddo enddo c c Interpolate u,v,t from x,y to r,t do j=0,ntcf do i=0,nrcf call lintcf(uxy,x1,dx,y1,dy,ixmx,iymx,nx,ny, + xrt(i,j),yrt(i,j),ucrt(i,j)) call lintcf(vxy,x1,dx,y1,dy,ixmx,iymx,nx,ny, + xrt(i,j),yrt(i,j),vcrt(i,j)) call lintcf(zxy,x1,dx,y1,dy,ixmx,iymx,nx,ny, + xrt(i,j),yrt(i,j),zrt(i,j)) call lintcf(txy,x1,dx,y1,dy,ixmx,iymx,nx,ny, + xrt(i,j),yrt(i,j),trt(i,j)) enddo enddo c c Convert Cartesian u,v to radial and tangential winds do j=0,ntcf do i=0,nrcf urrt(i,j) = ucrt(i,j)*ctheta(j) + vcrt(i,j)*stheta(j) vtrt(i,j) = vcrt(i,j)*ctheta(j) - ucrt(i,j)*stheta(j) enddo enddo c c Azimuthally average ur,vt, z and t call azavgcf(urrt,urr,mrcf,mtcf,nrcf,ntcf) call azavgcf(vtrt,vtr,mrcf,mtcf,nrcf,ntcf) call azavgcf( zrt,zaa,mrcf,mtcf,nrcf,ntcf) call azavgcf( trt,taa,mrcf,mtcf,nrcf,ntcf) c uraa(0) = 0.0 vtaa(0) = 0.0 do k=1,nrcf uraa(k) = urr(k) vtaa(k) = vtr(k) enddo c if (irflag .ne. 1) return c c ++ Vortex removal ++ c c Find the first radius beyond minimum where tangential wind c becomes negative (positive for SH) kstart = ifix(rminvr/dr) kend = nrcf c do k=kstart+1,nrcf vtrtem = vtr(k)*sgnlat if (vtrtem .lt. 0.0) then c write(6,*) 'k,vtr(k)',k,vtr(k) kend = k-1 exit endif enddo c rend = rcf(kend) vend = vtr(kend) uend = urr(kend) c if (kend .lt. nrcf) then c Set symmetric wind outside rend to zero do k=kend+1,nrcf vtr(k) = 0.0 urr(k) = 0.0 enddo endif c if (abs(vend) .gt. 1.0 .or. abs(uend) .gt. 1.0) then c Taper the last three points of symmetric circulation vtr(kend ) = 0.1*vtr(kend ) urr(kend ) = 0.1*urr(kend ) c vtr(kend-1) = 0.4*vtr(kend-1) urr(kend-1) = 0.4*urr(kend-1) c vtr(kend-2) = 0.7*vtr(kend-2) urr(kend-2) = 0.7*urr(kend-2) endif c c do k=0,nrcf c write(61,800) rcf(k),uraa(k),urr(k),vtaa(k),vtr(k) c 800 format(5(f6.1,1x)) c enddo c c Determine maximum range of x,y points for vortex removal xp = xcen + rend xm = xcen - rend yp = ycen + rend ym = ycen - rend c im = 1 + ifix( (xm-x1)/dx ) ip = 2 + ifix( (xp-x1)/dx ) jm = 1 + ifix( (ym-y1)/dy ) jp = 2 + ifix( (yp-y1)/dy ) c if (im .lt. 1) im = 1 if (ip .gt. nx) ip = nx if (jm .lt. 1) jm = 1 if (jp .gt. ny) jp = ny c c Set ipert=1 to put vortex only into uxy,vxy instead of c subtracting it ipert=0 if (ipert .eq. 1) then do j=1,ny do i=1,nx uxy(i,j) = 0.0 vxy(i,j) = 0.0 enddo enddo endif c c Loop through all x,y points that may need vortex removal do j=jm,jp do i=im,ip c c Calculate storm-relative r,theta for this x,y point xrel = xgk(i,j)-xcen yrel = ygk(i,j)-ycen c call ctor(xrel,yrel,rrel,thetarel) if (rrel .gt. rend) cycle !inner loop c c Linearly interpolate vortex radial and tangetial wind c to model x,y grid point i1 = ifix(rrel/dr) if (i1 .ge. nrcf-1) i1=nrcf-1 i2 = i1+1 rnorm = (rrel - dr*float(i1))/dr w1 = (1.0-rnorm) w2 = rnorm c urel = w1*urr(i1) + w2*urr(i2) vrel = w1*vtr(i1) + w2*vtr(i2) c c write(6,915) xrel,yrel,rrel,thetarel,rnorm,urel,vrel, c + vtr(i1),vtr(i2) c 915 format('x,y,r,t,rnorm,urel,vrel,vt1,vt2: ',10(f7.1,1x)) c c Convert urel,vrel to Cartesian coordinates sinth = sin(dtr*thetarel) costh = cos(dtr*thetarel) ucrel = urel*costh - vrel*sinth vcrel = vrel*costh + urel*sinth c c Subtract vortex if (ipert .eq. 1) then uxy(i,j) = ucrel vxy(i,j) = vcrel else uxy(i,j) = uxy(i,j)-ucrel vxy(i,j) = vxy(i,j)-vcrel endif c enddo enddo c return end subroutine vremove c subroutine ctor(x,y,r,theta) c This routine converts from Cartesion coordinates c to radial coordinates, where theta is in c degrees measured counter-clockwise from c the +x-axis. c r = sqrt(x*x + y*y) c if (r .le. 0.0) then theta = 0.0 return endif c rtd = 57.296 theta = rtd*acos(x/r) if (y .lt. 0.0) theta = 360.0 - theta c return end subroutine ctor c subroutine rtoc(r,theta,x,y) c This routine converts from radial coordinates c to Cartesian coordinates, where theta is in c degrees measured counter-clockwise from c the +x-axis. c rtd = 57.296 x = r*cos(theta/rtd) y = r*sin(theta/rtd) c return end subroutine rtoc c subroutine spdlim(rlat,spdmax) c This routine estimates the maximum expected translational c speed (kt) of an Atlantic or east Pacific tropical cyclone c at a given latitude. c alat = abs(rlat) c if (alat .lt. 25.0) then spdmax = 30.0 elseif (alat .ge. 25.0 .and. alat .lt. 40.0) then spdmax = 30.0 + (alat-25.0)*2.0 else spdmax = 60.0 endif c return end subroutine spdlim c c ** cfindxy code marker (end) ** c subroutine rhbc(rh,ptem) c This routine applies a bias correction the relative humidity c variable from GFS reanalysis fields to make them comparable c to those from the operational GFS analyses. c c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 C.DUP lsdiag,plevh parameter (nplevh=10) common /plevhs/ up,vp,zp,tp,rp,tbarp,plevh dimension up(nplevh),vp(nplevh),zp(nplevh) dimension tp(nplevh),rp(nplevh) dimension tbarp(nplevh) dimension plevh(nplevh) c dimension rh(ixmx,iymx) c c Local variables dimension a0(nplevh),a1(nplevh),a2(nplevh) c c Coefficients for the RH bias correction functions c (quadratic polynomials) data a0 /1.3886,6.0946,1.7537,4.1882,9.1137,13.624, + 0.0,0.0,0.0,0.0/ data a1 /1.1992,1.5182,1.6303,1.1531,0.9546,0.6662, + 1.0,1.0,1.0,1.0/ data a2 /-.0037,-.0072,-.0084,-.0045,-.0025,-.0004, + 0.0,0.0,0.0,0.0/ c c Find the pressure level eps = 0.0001 do k=1,nplevh if (abs(ptem-plevh(k)) .lt. eps) then kplev = k exit endif if (k .eq. nplevh) then c Pressure level not found, skip bias correction return endif enddo c b0 = a0(kplev) b1 = a1(kplev) b2 = a2(kplev) c do j=1,nlat1 do i=1,nlon1 rhx = rh(i,j) rh(i,j) = b0 + b1*rhx + b2*rhx*rhx enddo enddo c return end subroutine rhbc c function sgnf(x) if (x .ge. 0.0) then sgnf=1.0 else sgnf=-1.0 endif c return end function sgnf c subroutine gcrop(slon,slat,bufkm) c This routine crops the model grid to keep only c those points within bufkm of the specified point slat,slon. c C.DUP lsdiag,primd parameter(ixmx=401,iymx=201,ipmx=21) parameter(irmx=800,iamx=24) c C.DUP lsdiag,prima common /primv/ u,v,z,t,r dimension u(ixmx,iymx,ipmx),v(ixmx,iymx,ipmx),z(ixmx,iymx,ipmx) dimension t(ixmx,iymx,ipmx),r(ixmx,iymx,ipmx) c C.DUP lsdiag,latlon common /latlon/ RLATMN,RLATMX,RLONMN,RLONMX,DLAT,DLON, + RLATD,RLOND dimension rlatd(iymx),rlond(ixmx) common /ilatlon/ NLAT1,NLON1 c C.DUP lsdiag,plevh parameter (nplevh=10) common /plevhs/ up,vp,zp,tp,rp,tbarp,plevh dimension up(nplevh),vp(nplevh),zp(nplevh) dimension tp(nplevh),rp(nplevh) dimension tbarp(nplevh) dimension plevh(nplevh) c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c C.DUP lsdiag,index common /indexa/ indexp,indexv common /indexn/ indexvn parameter (nvar=5) dimension indexp(ipmx),indexv(nvar) character *1 indexvn(nvar) c c Local variables dimension workxyp(ixmx,iymx,ipmx),workx(ixmx),worky(iymx) c c write(6,*) 'gcrop with slon,slat,bufkm: ',slon,slat,bufkm c write(6,*) 'rlonf,l, nlon1: ',rlond(1),rlond(nlon1),nlon1 c write(6,*) 'rlatf,l, nlat1: ',rlatd(1),rlatd(nlat1),nlat1 c c Find indices of the model grid point to the lower left of c the storm center i1 = 1 + ifix( (slon-rlond(1))/dlon ) j1 = 1 + ifix( (slat-rlatd(1))/dlat ) c ibuf = 1 + ifix( bufkm/(dlon*dtk*cos(dtr*slat)) ) ilft = i1 - ibuf irgt = i1 + ibuf c if (ilft .lt. 1) ilft = 1 if (irgt .gt. nlon1) irgt = nlon1 c jbuf = 1 + ifix( bufkm/(dlat*dtk) ) jbot = j1 - jbuf jtop = j1 + jbuf c if (jbot .lt. 1) jbot = 1 if (jtop .gt. nlat1) jtop = nlat1 c c write(6,*) 'ilft,irgt,i1,ibuf,',ilft,irgt,i1,ibuf c write(6,*) 'jbot,jtop,j1,jbuf ',jbot,jtop,j1,jbuf c c Create cropped grid and copy variables c np = ipmx c c ++ Copy longitudes icntx = 1 do i=ilft,irgt workx(icntx) = rlond(i) icntx = icntx + 1 enddo nglon1 = icntx-1 nlon1 = nglon1 c do i=1,nglon1 rlond(i) = workx(i) enddo c c ++ Copy latitudes icnty = 1 do j=jbot,jtop worky(icnty) = rlatd(j) icnty = icnty + 1 enddo nglat1 = icnty-1 nlat1 = nglat1 c do j=1,nglat1 rlatd(j) = worky(j) enddo c c ++ Copy u(x,y,p) do k=1,np if (indexp(k) .eq. 1) then icnty = 1 do j=jbot,jtop icntx = 1 do i=ilft,irgt workxyp(icntx,icnty,k) = u(i,j,k) icntx = icntx + 1 enddo icnty = icnty + 1 enddo endif enddo c do k=1,np if (indexp(k) .eq. 1) then do j=1,nglat1 do i=1,nglon1 u(i,j,k) = workxyp(i,j,k) enddo enddo endif enddo c c ++ Copy v(x,y,p) do k=1,np if (indexp(k) .eq. 1) then icnty = 1 do j=jbot,jtop icntx = 1 do i=ilft,irgt workxyp(icntx,icnty,k) = v(i,j,k) icntx = icntx + 1 enddo icnty = icnty + 1 enddo endif enddo c do k=1,np if (indexp(k) .eq. 1) then do j=1,nglat1 do i=1,nglon1 v(i,j,k) = workxyp(i,j,k) enddo enddo endif enddo c c ++ Copy z(x,y,p) do k=1,np if (indexp(k) .eq. 1) then icnty = 1 do j=jbot,jtop icntx = 1 do i=ilft,irgt workxyp(icntx,icnty,k) = z(i,j,k) icntx = icntx + 1 enddo icnty = icnty + 1 enddo endif enddo c do k=1,np if (indexp(k) .eq. 1) then do j=1,nglat1 do i=1,nglon1 z(i,j,k) = workxyp(i,j,k) enddo enddo endif enddo c c ++ Copy T(x,y,p) do k=1,np if (indexp(k) .eq. 1) then icnty = 1 do j=jbot,jtop icntx = 1 do i=ilft,irgt workxyp(icntx,icnty,k) = t(i,j,k) icntx = icntx + 1 enddo icnty = icnty + 1 enddo endif enddo c do k=1,np if (indexp(k) .eq. 1) then do j=1,nglat1 do i=1,nglon1 t(i,j,k) = workxyp(i,j,k) enddo enddo endif enddo c c ++ Copy RH(x,y,p) do k=1,np if (indexp(k) .eq. 1) then icnty = 1 do j=jbot,jtop icntx = 1 do i=ilft,irgt workxyp(icntx,icnty,k) = r(i,j,k) icntx = icntx + 1 enddo icnty = icnty + 1 enddo endif enddo c do k=1,np if (indexp(k) .eq. 1) then do j=1,nglat1 do i=1,nglon1 r(i,j,k) = workxyp(i,j,k) enddo enddo endif enddo c c write(6,*) 'nglon1,lon1,lon2: ',nglon1,workx(1),workx(nglon1) c write(6,*) 'nglat1,lat1,lat2: ',nglat1,worky(1),worky(nglat1) c return END subroutine gcrop c subroutine gtwind(vtaap,taap,gtaap,rcf,pvp,slatgw,mrcf,nrcf,nvp) c This routine calculates the temperature associated with the c symmetric tangential wind using the gradient thermal wind equation c c Passed variables dimension vtaap(0:mrcf,nvp),taap(0:mrcf,nvp),gtaap(0:mrcf,nvp) dimension rcf(0:mrcf),pvp(nvp) c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot,rd c c Calculate Coriolis parameter f = 2.0*erot*sin(dtr*slatgw) c c Initialize gtw T with observed T at outer radius do ll=1,nvp gtaap(nrcf,ll) = taap(nrcf,ll) enddo c do ll=1,nvp c c Calculate pressure level differences in Pascals if (ll .eq. 1) then dpu = 0.0 dpl = (pvp(ll+1) - pvp(ll ))*100.0 elseif (ll. eq. nvp) then dpu = (pvp(ll ) - pvp(ll-1))*100.0 dpl = 0.0 else dpu = (pvp(ll )-pvp(ll-1))*100.0 dpl = (pvp(ll+1)-pvp(ll ))*100.0 endif c c Integrate radially inward at each pressure level do i=nrcf-1,0,-1 rdel = ( rcf(i) - rcf(i+1) )*1000.0 rbar = 0.5*( rcf(i) + rcf(i+1) )*1000.0 vbar = 0.5*( vtaap(i,ll)+vtaap(i+1,ll) ) vfac = -(pvp(ll)*100.0/rd)*(2.0*vbar/rbar + f) c c Calculate dv/dp if (ll .eq. 1) then vu = 0.0 vm = 0.5*(vtaap(i,ll ) + vtaap(i+1,ll )) vl = 0.5*(vtaap(i,ll+1) + vtaap(i+1,ll+1)) dvdp = (vl-vm)/dpl elseif (ll. eq. nvp) then vu = 0.5*(vtaap(i,ll-1) + vtaap(i+1,ll-1)) vm = 0.5*(vtaap(i,ll ) + vtaap(i+1,ll )) vl = 0.0 dvdp = (vm-vu)/dpu else vu = 0.5*(vtaap(i,ll-1) + vtaap(i+1,ll-1)) vm = 0.5*(vtaap(i,ll ) + vtaap(i+1,ll )) vl = 0.5*(vtaap(i,ll+1) + vtaap(i+1,ll+1)) dvdpu = (vm-vu)/dpu dvdpl = (vl-vm)/dpl wl = dpu/(dpu+dpl) wu = dpl/(dpu+dpl) dvdp = wl*dvdpl + wu*dvdpu endif c gtaap(i,ll) = gtaap(i+1,ll) + rdel*dvdp*vfac enddo enddo c return END subroutine gtwind c subroutine tpwcal(pptpw,tptpw,rptpw,nvp,pstpw,tstpw,rstpw,tpw) c This routine calculates the total precipitable water in mm c from a sounding c dimension pptpw(nvp),tptpw(nvp),rptpw(nvp) c c Local variables parameter (mxp=200) dimension pl(mxp),tl(mxp),rl(mxp) dimension esl(mxp),wsl(mxp),wl(mxp) c c Constants eps = 0.622 g = 9.81 rhow = 1000.0 c c Assemble sounding do k=1,nvp pl(k) = pptpw(k) tl(k) = tptpw(k) rl(k) = rptpw(k) enddo c c Put surface values in their proper place nl = 0 do k=1,nvp-1 if (pstpw .gt. pptpw(k ) .and. + pstpw .le. pptpw(k+1) ) then nl = k+1 pl(k+1) = pstpw tl(k+1) = tstpw rl(k+1) = rstpw c exit endif if (k .eq. nvp-1) then c If you get to here, the surface pressure is higher c than the highest pressure in the sounding, c so tack the surface values on to the end. nl = nvp + 1 pl(nl) = pstpw tl(nl) = tstpw rl(nl) = rstpw endif enddo c c Check RH values (some GFS analysis don't have values above 300 hPa) c and convert P values to Pa do k=1,nl if (rl(k) .gt. 100.0 .or. rl(k) .le. 0.0) rl(k) = 40.0 pl(k) = 100.0*pl(k) enddo c do k=1,nl c Calculate saturation vapor pressure (pa) call teten(tl(k),esl(k)) c c Calculate saturation mixing ratio and mixing ratio wsl(k) = esl(k)*eps/(pl(k)-esl(k)) wl(k) = rl(k)*wsl(k)/100.0 enddo c c Integrate sounding to find TPW tpw = 0.0 do k=2,nl tpw = tpw + 0.5*(wl(k)+wl(k-1))*(pl(k)-pl(k-1)) enddo c tpw = 1000.0*tpw/(rhow*g) c return END subroutine tpwcal