program lsdiag c This program performs diagnostic calculations using c large-scale analyses in packed form. 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 a c 16 Sep 2006 - New storm-relative D200 calculation c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,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.,500., + 400., 350.,300.,250.,200.,150.,100.,70.,50., + 25., 10., 5., 1./ DATA PLEVR / 1013.,1000.,1000.,925.,850.,850.,700.,500., + 400., 400.,300.,250.,200.,150.,100.,70.,50., + 25., 10., 5., 1./ c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot 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 Calculate vertical shear from aviation model fields c call ascal c c Calculate SHIPS variables (2005-2006 version with 5-day input files c 6 hr data) call sv6cal c c Calculate time averages of specified variables c call tavg c stop END subroutine aunpack(ierr) 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=121,iymx=91,ipmx=21) parameter(irmx=300,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=10000) 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 10 i=1,ipmx indexp(i) = 0 10 continue c do 15 i=1,nvar indexv(i) = 0 15 continue 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 IF (rawdfn(1:2) .EQ. 'WP') THEN RLMIN = RLONMN RLMAX = RLONMX RLONMN = -RLMAX RLONMX = -RLMIN 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 500 CONTINUE 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 25 I=1,ipmx IF (PTEM .EQ. PLEV(I)) GO TO 1000 25 CONTINUE C WRITE(lulog,960) PTEM 960 FORMAT(2X,'PRESSURE= ',E11.4,' NOT FOUND IN PLEV ARRAY') STOP c 1000 CONTINUE indexp(i) = 1 lnow = i c C READ IN PACKED DATA DO 30 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)) 30 CONTINUE c C UNPACK THE DATA AND PUT IT IN THE PROPER ARRAY DO 35 I=1,IPTS IX = IDECOD(CODE(I)) TRA(I) = IX * SMPY - BSUB 35 CONTINUE C if (type .eq. 'U') then indexv(1) = 1 DO 40 I = 0,NLATI DO 40 J = 0,NLONI IJ = (J+1) + I*(NLONI+1) u(J+1,I+1,LNOW) = TRA(IJ) 40 CONTINUE elseif (type .eq. 'V') then indexv(2) = 1 DO 45 I = 0,NLATI DO 45 J = 0,NLONI IJ = (J+1) + I*(NLONI+1) v(J+1,I+1,LNOW) = TRA(IJ) 45 CONTINUE elseif (type .eq. 'Z') then indexv(3) = 1 DO 50 I = 0,NLATI DO 50 J = 0,NLONI IJ = (J+1) + I*(NLONI+1) z(J+1,I+1,LNOW) = TRA(IJ) 50 CONTINUE elseif (type .eq. 'T') then indexv(4) = 1 DO 55 I = 0,NLATI DO 55 J = 0,NLONI IJ = (J+1) + I*(NLONI+1) t(J+1,I+1,LNOW) = TRA(IJ) 55 CONTINUE elseif (type .eq. 'R') then indexv(5) = 1 DO 60 I = 0,NLATI DO 60 J = 0,NLONI IJ = (J+1) + I*(NLONI+1) r(J+1,I+1,LNOW) = TRA(IJ) 60 CONTINUE else write(lulog,970) type 970 format(/,' Unrecognized variable type: ',a1) c stop endif go to 500 600 continue c c Calculate latitude (deg) and longitudes (deg W neg.) c the of grid points do 65 j=1,nlat1 rlatd(j) = rlatmn + dlat*float(j-1) 65 continue c do 70 i=1,nlon1 rlond(i) = -(rlonmx - dlon*float(i-1)) 70 continue 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 75 i=1,ipmx if (indexp(i) .eq. 1) then write(lulog,310) plevr(i) 310 format(1x,f6.1) endif 75 continue c vlabel = ' ' ii = 0 do 80 i=1,nvar if (indexv(i) .eq. 1) then ii = ii+1 ii2 = 2*ii vlabel(ii2:ii2) = indexvn(i) endif 80 continue 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 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 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 subroutine ascal c This routine calculates the vertical shear from the Aviation c model fields for each case in the input file. c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,iamx=24) 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,rawdfn common /rawdf/ rawdfn character *40 rawdfn C.DUP lsdiag,unnum common /unitno/ lulog,luinp,luout,lurawd c c Local common character *72 hline,tline,iline(12) c dimension islat(-2:7),islon(-2:7) c c ******************** c c Generic rawd file name for Aviation files rawdfn = 'AF0000_Z0000_PACK.DAT' c c Start the case by case loop nrawd=0 2001 continue c C READ BEST TRACK FILE FOR CASE SELECTION READ(LUINP,100,ERR=751,END=2000) HLINE 100 FORMAT(A72) c GO TO 5100 751 WRITE(LULOG,*) ' ERROR READING INPUT FILE (LUINP)' STOP 5100 CONTINUE C iline(1) = hline DO 10 I=2,12 READ(LUINP,100,ERR=751,END=752) ILINE(I) GO TO 5200 752 WRITE(LULOG,*) ' UNEXPECTED END ON INPUT FILE (LUINP)' STOP 5200 CONTINUE 10 CONTINUE C nrawd = nrawd + 1 WRITE(6,200) NRAWD,HLINE 200 FORMAT(/,' START CASE',I4,' FIRST LINE OF LUINP HEADER:',/,A72) C C READ STORM POSITIONS OFF OF ILINE FROM LUINP FILE READ(ILINE(3),105) (ISLAT(K),K=-2,6) READ(ILINE(4),105) (ISLON(K),K=-2,6) 105 FORMAT(9(1X,I4)) C c Calculate vertical shear if necessary do 15 k=0,6 if (islat(k) .lt. 9000 .and. islon(k) .lt. 9000) then slat = float(islat(k))/10.0 slon = -1.0*float(islon(k))/10.0 c c Create Aviation file name k12 = k*12 write(rawdfn(3:4),405) k12 400 format(i2) c rawdfn(5:6) = HLINE(7:8) rawdfn(9:10) = HLINE(9:10) c read(hline(11:12),400) iday read(hline(14:15),400) itime if (itime .eq. 12) iday = iday+50 write(rawdfn(11:12),405) iday 405 format(i2.2) c c Unpack the rawd file call aunpack(ierr) c if (ierr .ne. 0) then write(6,205) 205 format(' Error opening or reading data input file') stop endif c c Calculate the vertical shear ps1 = 850.0 ps2 = 200.0 radkt = 600.0 minpts = 8 call vshear(ps1,ps2,slon,slat,radkt,minpts,shear,sheard) c ishear = ifix(0.49 + 10.0*shear) c tline = iline(5) ks = 12 + 5*k ke = ks + 3 write(tline(ks:ke),410) ishear 410 format(i4) iline(5) = tline endif 15 continue c do 30 i=1,12 write(luout,100) iline(i) 30 continue c C DO NEXT CASE GO TO 2001 C 2000 CONTINUE C WRITE(LULOG,350) NRAWD 350 FORMAT(/,' NORMAL COMPLETION OF LSDIAG, ',I4,' CASES PROCESSED') C RETURN END subroutine tavg c This routine finds the time average of various variables c from the files listed in lsdiag.in c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,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 c Other variables dimension upavg(nplevh),vpavg(nplevh),tpavg(nplevh), + zpavg(nplevh),rpavg(nplevh) c dimension tppavg(nplevh),deltp(nplevh) parameter (mxd=300) dimension pps(mxd),tts(mxd),tvs(mxd) c ******************** c c Generic rawd file name for the input rawd files rawdfn = 'BACK00_X0000_PACK.DAT' c c Specify center lat,lon for sounding clat = 12.0 clon = -40.0 c clat = 13.0 clon = -40.0 c c Set istavg=1 to average sounding over circular area of radius radk c or istavg=2 to avarage sounding over rectangular area +/- xradk,yr istavg=2 c radk = 666.0 xradk=1085.0 yradk=555.6 c write(luout,203) clat,clon 203 format(/,' Center of area average, lat,lon: ', + f5.1,1x,f6.1) c if (istavg .eq. 1) then write(luout,201) radk 201 format(/,' Radius of area=',f6.1,' km') else write(luout,202) xradk,yradk 202 format(/,' Rectangular area +/- x=',f6.1,' km ', + '+/- y=',f6.1,' km ') endif c c Initialize average variables to zero do 5 i=1,nplevh upavg(i) = 0.0 vpavg(i) = 0.0 tpavg(i) = 0.0 zpavg(i) = 0.0 rpavg(i) = 0.0 5 continue c c Start the case by case loop nrawd=0 2001 continue c C READ FILE FOR CASE SELECTION READ(LUINP,100,ERR=751,END=2000) rawdfn 100 FORMAT(A40) c GO TO 5100 751 WRITE(LULOG,*) ' ERROR READING INPUT FILE (LUINP)' STOP 5100 CONTINUE C nrawd = nrawd + 1 WRITE(6,200) NRAWD,rawdfn WRITE(LUOUT,200) NRAWD,rawdfn 200 FORMAT(/,' START CASE',I4,' FILE: ',A40) C c Unpack the rawd file call aunpack(ierr) c if (ierr .ne. 0) then write(6,205) 205 format(' Error opening or reading data input file') nrawd = nrawd-1 go to 2001 endif c c Calculate vertical soundings if (istavg .eq. 1) then call vsound(clon,clat,radk) else call vcsound(clon,clat,xradk,yradk) endif c if (rawdfn(1:4) .eq. 'BACK') then c Convert virtual temperature to temperature do 8 i=1,nplevh if (rp(i) .le. 100.0) then call virti(tp(i),plevh(i),rp(i),tt) tp(i) = tt endif 8 continue endif c c Convert deviation heights to total heights do 9 i=1,nplevh ptemp = plevh(i) call stndz(ptemp,ztemp,ttemp,thtemp) zp(i) = zp(i) + ztemp 9 continue c c Add current case for time average do 10 i=1,nplevh upavg(i) = upavg(i) + up(i) vpavg(i) = vpavg(i) + vp(i) tpavg(i) = tpavg(i) + tp(i) zpavg(i) = zpavg(i) + zp(i) rpavg(i) = rpavg(i) + rp(i) 10 continue c c Write current sounding to the output file write(luout,210) (up(i),i=1,nplevh) write(luout,215) (vp(i),i=1,nplevh) write(luout,220) (tp(i),i=1,nplevh) write(luout,224) (zp(i),i=1,nplevh) write(luout,222) (rp(i),i=1,nplevh) write(luout,225) (plevh(i),i=1,nplevh) 210 format(' U ',10(f6.1)) 215 format(' V ',10(f6.1)) 220 format(' T ',10(f6.1)) 221 format(' TP',10(f6.1)) 223 format(' DT',10(f6.1)) 224 format(' Z ',10(f6.0)) 222 format(' R ',10(f6.1)) 225 format(' P ',10(f6.1)) c C Do next case go to 2001 C 2000 continue c c Calculate time averages rnrawd = float(nrawd) if (nrawd .le. 0) rnrawd = 1.0 c do 20 i=1,nplevh upavg(i) = upavg(i)/rnrawd vpavg(i) = vpavg(i)/rnrawd tpavg(i) = tpavg(i)/rnrawd zpavg(i) = zpavg(i)/rnrawd rpavg(i) = rpavg(i)/rnrawd 20 continue c c Calculate surface pressure and temperature z1 = zpavg(1) z2 = zpavg(2) p1 = plevh(1)*100.0 p2 = plevh(2)*100.0 t1 = tpavg(1) + 273.15 t2 = tpavg(2) + 273.15 call pstcal(z1,z2,t1,t2,p1,psfc,tsfc) c c Calculate lifted parcel temperature from sfc to 100 mb p1000 = 1000.0e+2 t1000 = tpavg(1)+273.13 r1000 = rpavg(1) dpp = 10.0e+2 ptop = 100.0e+2 iadj = 1 c if (nrawd .gt. 0) then if (psfc .lt. p1000) then call parsub(p1000,t1000,r1000,dpp,ptop,iadj,mxd, + psat,pps,tts,tvs,npp) c else call parsub(psfc,tsfc,r1000,dpp,ptop,iadj,mxd, + psat,pps,tts,tvs,npp) endif endif c do kk=1,nplevh do ii=1,npp if (pps(ii) .eq. 100.0*plevh(kk)) then tppavg(kk) = tts(ii)-273.15 endif enddo deltp(kk) = tppavg(kk) - tpavg(kk) enddo c c Sum up the temperature differences from level 1 to npmax npmax = 8 dptot = plevh(1)-plevh(npmax) dttp = 0.0 do kk=1,npmax if (kk .eq. 1) then wtt = 0.5*(plevh( 1)-plevh( 2))/dptot elseif (kk .eq. npmax) then wtt = 0.5*(plevh(npmax-1)-plevh(npmax))/dptot else wtt = 0.5*(plevh( kk-1)-plevh( kk+1))/dptot endif c dttp = dttp + wtt*(tppavg(kk) - tpavg(kk)) enddo c delu = upavg(8) - upavg(2) delt = tppavg(8) - tpavg(8) c write(luout,230) nrawd 230 format(//,' Time averaged variables, no. cases= ',i4) write(luout,210) (upavg(i),i=1,nplevh) write(luout,215) (vpavg(i),i=1,nplevh) write(luout,220) (tpavg(i),i=1,nplevh) write(luout,221) (tppavg(i),i=1,nplevh) write(luout,223) (deltp(i),i=1,nplevh) write(luout,224) (zpavg(i),i=1,nplevh) write(luout,222) (rpavg(i),i=1,nplevh) write(luout,225) (plevh(i),i=1,nplevh) C write(luout,227) upavg(2),upavg(8),delu,rpavg(1),delt,tsfc-273.15, + tpavg(8),psfc/100.0,psat/100.0,dttp, + nrawd,rawdfn(5:12) 227 format(/,'u2=',f4.0,' u8=',f4.0,' du=',f5.1, + ' rh1=',f3.0,' dT=',f4.1, + ' Ts=',f4.1,' T8=',f5.1,' Ps=',f6.1,' Pl=',f5.0, + ' S=',f4.1,1x,i2,1x,a8) c WRITE(LULOG,350) NRAWD 350 FORMAT(/,' NORMAL COMPLETION OF LSDIAG, ',I4,' CASES PROCESSED') C RETURN END subroutine sv6cal c This routine calculates the variables needed for c the SHIPS intensity forecast model (2005-2006 version) c for each case in the input file. c c This version uses the input file with 5-day forecasts. c c The lsdiag input file includes data at 6 hour intervals c rather than 12 hour intervals for this version c parameter (itmax=20) c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,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 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) c dimension iu200(0:itmax),it200(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 *4 sname,tlabel character *130 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 Vertical profile arrays for Kerry Emanuel MPI code dimension pkmpi(nvp),tkmpi(nvp),qkmpi(nvp),rkmpi(nvp) dimension irsst(-2:itmax),rsst(-2:itmax) dimension ivkmpi(0:itmax),ipkmpi(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) 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) 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) c dimension zcrt(0:mrcf,0:mtcf),tcrt(0:mrcf,0:mtcf) dimension zraa(0:mrcf),traa(0:mrcf) c 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 iray = 2000 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 fluxe 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 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 5 n=2,nvp-1 wvp(n) = ( 0.5*(pvp(n+1)-pvp(n-1)) )/dpvp 5 continue 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 2001 continue c c Read header line for this case read(luinp,100,err=751,end=2000) hline 100 format(a130) 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) c write(lulog,748) sname,imyear,immon,imday,imtime 748 format(/,' HEADER: ',a4,1x,3(i2.2),1x,i2.2) c if (imyear .lt. 50) then imyear4 = 2000 + imyear else imyear4 = 1900 + imyear endif c go to 5100 751 write(lulog,*) ' Error reading input file (luinp)' stop 5100 continue C iline(1) = hline c do k=-2,itmax irsst(k) = 9999 enddo c c Read the rest of the data lines for this case c and look for the LAT and LON lines. lcount=2 5200 continue read(luinp,100,err=751,end=751) tline iline(lcount) = tline c if (tline(117:120) .eq. 'LAT ') then read(tline,105) (islat(k),k=-2,itmax) 105 format(23(1x,i4)) endif C if (tline(117:120) .eq. 'LON ') then read(tline,105) (islon(k),k=-2,itmax) endif c if (tline(117:120) .eq. 'RSST') then read(tline,105) (irsst(k),k=-2,itmax) endif c if (tline(117:120) .eq. 'LAST') go to 5201 c lcount=lcount+1 go to 5200 5201 continue c nrawd = nrawd + 1 write(6,200) nrawd,hline 200 format(/,' START CASE',I4,' FIRST LINE OF LUINP HEADER:',/,A80) C do 15 k=-2,itmax slat(k) = float(islat(k))/10.0 slon(k) = -1.0*float(islon(k))/10.0 c if (irsst(k) .lt. 9000) then rsst(k) = 0.1*float(irsst(k)) else rsst(k) = 9999. endif 15 continue c c Initialize synoptic variables to missing do 20 kt=0,itmax iu200(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 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 ivkmpi(kt) = 9999 ipkmpi(kt) = 9999 20 continue c ilost=0 tlast=-6.0 c c Start time loop do 99 kt=0,itmax tnow = 6.0*float(kt) c c Check for missing track positions if (slat(kt) .gt. 900.0) go to 99 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 rawdfn = 'A00055_X0000_PACK.DAT' endif c if (iatime .eq. 6 .or. iatime .eq. 18) then rawdfn(8:8) = 'Y' endif c write(rawdfn(5:6),115) imyear 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 3001 continue rawdfn = 'A00055_X0000_PACK.DAT' 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 call aunpack(ierr) if (ierr .eq. 0) go to 3002 c ilaghr = ilaghr + 6 c c Set the limit on how far back to look for previous mode if (ilaghr .gt. 6) go to 3002 c call tadd(imyear,immon,imday,imtime,-ilaghr, + ityear,itmon,itday,ittime) go to 3001 endif endif c 3002 continue c c Unpack the current rawd file call aunpack(ierr) 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 go to 2001 endif c if (ierr .eq. 1) then c If subsequent data files are missing, skip c calculation only for individual time periods go to 99 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) 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 else if (slat(kt+2) .lt. 900.0 .and. + slat(kt ) .lt. 900.0 .and. + kt+2 .le. itmax) then t2lat = slat(kt+2) t2lon = slon(kt+2) t1lat = slat(kt) t1lon = slon(kt) deltc = 2.0*delt icok = 1 elseif (slat(kt+1) .lt. 900.0 .and. + slat(kt ) .lt. 900.0 .and. + kt+1 .le. itmax) then t2lat = slat(kt+1) t2lon = slon(kt+1) t1lat = slat(kt) t1lon = slon(kt) deltc = 1.0*delt icok = 1 else t2lat = slat(kt) t2lon = slon(kt) t1lat = slat(kt) t1lon = slon(kt) deltc = 1.0*delt icok = 0 endif 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 c Get vertical profiles of required variables do 40 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) 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 40 continue 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 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(kt) = ifix(pbars) 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 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 Kery 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 if (kt .eq. 0) then write(lulog,822) sstkmpi,pskmpi, + ipkmpi(kt),ivkmpi(kt),ihvmax 822 format(/,'SST= ',f6.1,' PSFC=',f6.1, + ' PMPI=',i4,' VMPI=',i4,' VMAX=',i4) 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 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) c itlat(kt) = ifix( 10.0*slatcf) itlon(kt) = ifix(-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) 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 call lfilter(txyp(1,1,ll),xcf,ycf,fradlf) call lfilter(rxyp(1,1,ll),xcf,ycf,fradlf) call lfilter(zxyp(1,1,ll),xcf,ycf,fradlf) c call vremove(uxyp(1,1,ll),vxyp(1,1,ll),rcf,thetacf, + xcf,ycf,mrcf,mtcf,nrcf,ntcf, + uraa,vtaa,irmflag) c if (ll .eq. i200) then c Calculate storm-relative 200 mb area average diverge 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 wi c c Find last radial point where vtaa .ge. 1 starting at do kk=4,nrcf nrlast = kk if (vtaa(kk) .lt. 1.0) then nrlast = kk-1 go to 1125 endif enddo 1125 continue rtlast = rcf(nrlast) c vtmax = -1000.0 rtmax = 1000.0 vtavg = 0.0 do kk=0,nrlast vtavg = vtavg + vtaa(kk) if (vtaa(kk) .gt. vtmax) then vtmax = vtaa(kk) rtmax = rcf(kk) endif enddo vtavg = vtavg/(float(nrlast+1)) c irtcc(kt) = ifix(rtlast) irtxc(kt) = ifix(rtmax) itwxc(kt) = ifix(vtmax*10.0) itwbc(kt) = ifix(vtavg*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 endif 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 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 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 99 continue c c Write input lines to output file (except the last line) do 16 i=1,lcount-1 write(luout,100) iline(i) 16 continue 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,21(1x,i4),1x,'U200') write(luout,212) (ie000(k),k=0,itmax) 212 format(10x,21(1x,i4),1x,'E000') write(luout,213) (iepos(k),k=0,itmax) 213 format(10x,21(1x,i4),1x,'EPOS') write(luout,214) (ieneg(k),k=0,itmax) 214 format(10x,21(1x,i4),1x,'ENEG') write(luout,215) (iepss(k),k=0,itmax) 215 format(10x,21(1x,i4),1x,'EPSS') write(luout,216) (ienss(k),k=0,itmax) 216 format(10x,21(1x,i4),1x,'ENSS') write(luout,217) (irhlo(k),k=0,itmax) 217 format(10x,21(1x,i4),1x,'RHLO') write(luout,218) (irhmd(k),k=0,itmax) 218 format(10x,21(1x,i4),1x,'RHMD') write(luout,219) (irhhi(k),k=0,itmax) 219 format(10x,21(1x,i4),1x,'RHHI') c write(luout,231) (ipslv(k),k=0,itmax) 231 format(10x,21(1x,i4),1x,'PSLV') write(luout,232) (ivorb(k),k=0,itmax) 232 format(10x,21(1x,i4),1x,'Z850') write(luout,233) (idivb(k),k=0,itmax) 233 format(10x,21(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,21(1x,i4),1x,'TLAT') write(luout,240) (itlon(k),k=0,itmax) 240 format(10x,21(1x,i4),1x,'TLON') write(luout,241) (itwac(k),k=0,itmax) 241 format(10x,21(1x,i4),1x,'TWAC') c write(luout,260) (itwbc(k),k=0,itmax) c 260 format(10x,21(1x,i4),1x,'TWBC') write(luout,261) (itwxc(k),k=0,itmax) 261 format(10x,21(1x,i4),1x,'TWXC') c write(luout,262) (irtcc(k),k=0,itmax) c 262 format(10x,21(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,21(1x,i4),1x,'PENC') write(luout,242) (ishdc(k),k=0,itmax) 242 format(10x,21(1x,i4),1x,'SHDC') write(luout,243) (isddc(k),k=0,itmax) 243 format(10x,21(1x,i4),1x,'SDDC') write(luout,244) (ishgc(k),k=0,itmax) 244 format(10x,21(1x,i4),1x,'SHGC') c write(luout,245) (it15c(k),k=0,itmax) c 245 format(10x,21(1x,i4),1x,'T15C') c write(luout,246) (it20c(k),k=0,itmax) c 246 format(10x,21(1x,i4),1x,'T20C') c write(luout,247) (it25c(k),k=0,itmax) c 247 format(10x,21(1x,i4),1x,'T25C') c write(luout,248) (iepoc(k),k=0,itmax) c 248 format(10x,21(1x,i4),1x,'EPOC') write(luout,249) (idivc(k),k=0,itmax) 249 format(10x,21(1x,i4),1x,'DIVC') c endif c write(luout,311) (it150(k),k=0,itmax) 311 format(10x,21(1x,i4),1x,'T150') write(luout,312) (it200(k),k=0,itmax) 312 format(10x,21(1x,i4),1x,'T200') write(luout,313) (it250(k),k=0,itmax) 313 format(10x,21(1x,i4),1x,'T250') write(luout,220) (ishrd(k),k=0,itmax) 220 format(10x,21(1x,i4),1x,'SHRD') write(luout,222) (ishtd(k),k=0,itmax) 222 format(10x,21(1x,i4),1x,'SHTD') write(luout,225) (ishrs(k),k=0,itmax) 225 format(10x,21(1x,i4),1x,'SHRS') write(luout,227) (ishts(k),k=0,itmax) 227 format(10x,21(1x,i4),1x,'SHTS') write(luout,228) (ishrg(k),k=0,itmax) 228 format(10x,21(1x,i4),1x,'SHRG') c if (imxshr .eq. 1) then write(luout,229) (ishxu(k),k=0,itmax) 229 format(10x,21(1x,i4),1x,'SHXU') write(luout,230) (ishxl(k),k=0,itmax) 230 format(10x,21(1x,i4),1x,'SHXL') endif c tlabel='PENV' write(luout,236) (ipsfc(k),k=0,itmax),tlabel 236 format(10x,21(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 c Write out last line write(luout,100) iline(lcount) c C Do next case go to 2001 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') C return END 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 the c 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 ssfcal(slon,slat,ssx,ssy,ssize,shear,refc,pefc) c This routine calculates the storm motion vector (ssx,ssy), c storm size (ssize), vertical shear along the storm track c (shear), and the relative (refc) and planetary (pefc) eddy c flux convergence of angular momentum. c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,iamx=24) 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,pncons common /pncons/ pi,dtr,dtk,erad,erot 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 slat(-2:7), slon(-2:7) dimension shear(0:6) c c f(x,y) variables dimension u200(ixmx,iymx),v200(ixmx,iymx) dimension u850(ixmx,iymx),v850(ixmx,iymx) c dimension us200(ixmx,iymx),vs200(ixmx,iymx) dimension us850(ixmx,iymx),vs850(ixmx,iymx) c dimension ur850(ixmx,iymx),vt850(ixmx,iymx) dimension ur200(ixmx,iymx),vt200(ixmx,iymx) c dimension ur850cb(ixmx,iymx),vt850cb(ixmx,iymx) dimension ur200cb(ixmx,iymx),vt200cb(ixmx,iymx) c dimension u850b(ixmx,iymx),v850b(ixmx,iymx) dimension u200b(ixmx,iymx),v200b(ixmx,iymx) c c f(r,a) variables dimension ur850p(irmx,iamx),vt850p(irmx,iamx) dimension ur200p(irmx,iamx),vt200p(irmx,iamx) c dimension f(irmx,iamx),t1ra(irmx,iamx),t2ra(irmx,iamx) c c f(r) variables dimension ur850b(irmx),vt850b(irmx) dimension ur200b(irmx),vt200b(irmx) c dimension fb(irmx),pefc(irmx),refc(irmx) dimension t1r(irmx),t2r(irmx) c c ******************** c c ++ Shear calculation parameters ++ c c Specify radius (km) for shear calculation radshr = 600.0 c c Specify minimum points for shear calculation minpts = 10 c c ++ Size calculation parameters ++ rsmin = 395.0 rsmax = 805.0 vnorm = 5.0 c c ++ End parameter specification c slat0 = slat(0) slon0 = slon(0) c c Calculate grid information relative to initial storm position call xyracal(slon0,slat0) c c Calculate x,y components of storm motion IF (SLAT(-1) .LT. 900.0 .AND. SLAT(0) .LT. 900.0) THEN T2LAT = SLAT(0) T2LON = SLON(0) T1LAT = SLAT(-1) T1LON = SLON(-1) ELSEIF (SLAT(0) .LT. 900.0 .AND. SLAT(1) .LT. 900.0) THEN T2LAT = SLAT(1) T2LON = SLON(1) T1LAT = SLAT(0) T1LON = SLON(0) ELSE WRITE(LULOG,904) 904 FORMAT(/,' NOT ENOUGH VALID TRACK POINTS FOR STORM SPEED') STOP ENDIF C call lltoxy(olons,olats,colat,t1lon,t1lat,temx1,temy1) call lltoxy(olons,olats,colat,t2lon,t2lat,temx2,temy2) c ssx = 1000.0*(temx2-temx1)/(12.0*3600.0) ssy = 1000.0*(temy2-temy1)/(12.0*3600.0) c c Calculate Coriolis parameter at all of the polar grid points do 9 j=1,ntheta do 9 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) 9 continue c c Get 200 and 850 winds call varget('U',200.0,u200,1,iexist) call varget('V',200.0,v200,1,iexist) call varget('U',850.0,u850,1,iexist) call varget('V',850.0,v850,1,iexist) c c Calculate radial and tangential winds at 850,200 call uvtort(u200,v200,ur200,vt200) call uvtort(u850,v850,ur850,vt850) c c Interpolate 850,200 mb radial,tangential wind to polar grid call ctop(ur200,ur200p) call ctop(vt200,vt200p) call ctop(ur850,ur850p) call ctop(vt850,vt850p) c c Azimuthally average polar grid winds call azavg(ur200p,ur200b,0) call azavg(vt200p,vt200b,0) call azavg(ur850p,ur850b,0) call azavg(vt850p,vt850b,0) c c Azimuthally average Coriolis parameter call azavg(f,fb,1) c c Calculate storm size anorm = 0.0 aram = 0.0 do 20 i=1,nrad if (radk(i) .gt. rsmin .and. radk(i) .lt. rsmax) then aram = aram + vt850b(i)*radk(i)*radk(i) anorm = anorm + vnorm*radk(i)*radk(i) endif 20 continue if (anorm .gt. 0.0) then ssize = aram/anorm else ssize=0.0 endif c write(lulog,301) ssize 301 format(/,' Size= ',f6.2) c c Calculate planetary eddy flux momentum convergence c in units of m/s/day do 22 j=1,ntheta do 22 i=1,nrad t1ra(i,j) = -(ur200p(i,j)-ur200b(i))*(f(i,j)-fb(i)) 22 continue call azavg(t1ra,pefc,0) c do 23 i=1,nrad pefc(i) = pefc(i)*24.0*3600.0 23 continue c c Calculate relative eddy flux momentum convergence c c Calculate -U'V' in storm relative coordinates do 26 j=1,ntheta ctheta = cos(dtr*thetad(j)) stheta = sin(dtr*thetad(j)) ssrad = ssx*ctheta + ssy*stheta sstan =-ssx*stheta + ssy*ctheta do 28 i=1,nrad t1ra(i,j) = -(ur200p(i,j) - (ssrad + ur200b(i)))* + (vt200p(i,j) - (sstan + vt200b(i))) 28 continue 26 continue c c Azimuthally average -U'V' call azavg(t1ra,t1r,0) c c Take radial derivative of averaged -U'V'. call ddrcal(t1r,t2r) c c Calculate relative flux convergence in m/s/day refc(1) = 24.*3600.*t2r(1) do 29 i=2,nrad refc(i) = 24.*3600.*(2.*t1r(i)/(1000.*radk(i)) + t2r(i)) 29 continue c write(lulog,302) 302 format(/,' r(km) ur850b vt850b ur200b vt200b refc pefc') do 15 i=1,nrad write(lulog,304) radk(i),ur850b(i),vt850b(i), + ur200b(i),vt200b(i), + refc(i),pefc(i) 304 format(1x,f7.1,6(f6.1,1x)) 15 continue c c Interpolate azimuthally averaged radial,tangential c wind back to Cartesian grid. call patoc(ur200b,ur200cb) call patoc(vt200b,vt200cb) call patoc(ur850b,ur850cb) call patoc(vt850b,vt850cb) c c Convert radial and tangetial winds to Cartesian components call rttouv(ur200cb,vt200cb,u200b,v200b) call rttouv(ur850cb,vt850cb,u850b,v850b) c c Subtract azimuthally averaged flow from u,v do 25 j=1,ny do 25 i=1,nx us200(i,j) = u200(i,j)-u200b(i,j) vs200(i,j) = v200(i,j)-v200b(i,j) us850(i,j) = u850(i,j)-u850b(i,j) vs850(i,j) = v850(i,j)-v850b(i,j) 25 continue c c Calculate vertical shear do 30 k=0,6 c Calculate grid information slatt = slat(k) slont = slon(k) if (slatt .gt. 900.0 .or. slont .gt. 900.0) then ua200 = 999.9 va200 = 999.9 ua850 = 999.9 va850 = 999.9 shrtem = 999.9 go to 1030 endif c if (k .gt. 0) call xyracal(slont,slatt) c c Horizontally average 850,200 mb winds call havg(us200,radshr,minpts,0,ierr1,ua200) call havg(vs200,radshr,minpts,0,ierr2,va200) call havg(us850,radshr,minpts,0,ierr3,ua850) call havg(vs850,radshr,minpts,0,ierr4,va850) c c Check for missing data ierrs = ierr1 + ierr2 + ierr3 + ierr4 if (ierrs .gt. 0) then ua200 = 999.9 va200 = 999.9 ua850 = 999.9 va850 = 999.9 shrtem = 999.9 go to 1030 endif c c Convert average winds to kts ua200 = 1.944*ua200 va200 = 1.944*va200 ua850 = 1.944*ua850 va850 = 1.944*va850 c shrtem = sqrt( (ua200-ua850)**2 + (va200-va850)**2 ) 1030 continue c shear(k) = shrtem ttime = 12.0*float(k) c write(lulog,303) ttime,slatt,slont 303 format(/' time,lat,lon: ',3(f6.1,1x)) c write(lulog,300) ua850,va850,ua200,va200,shrtem 300 format(' u850,v850,u200,v200,shr: ',5(f6.1,1x)) 30 continue c return END subroutine pncal c This routine specifies physical and numerical constants c C.DUP lsdiag,pncons common /pncons/ pi,dtr,dtk,erad,erot c pi = 3.14159265 dtr = pi/180.0 c dtk = 111.1 c erad = 6371.0 erot = 7.292E-5 c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 i=1,nx u(i,j) = u(i,j) - cx v(i,j) = v(i,j) - cy 10 continue 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 20 j=1,ntheta do 20 i=1,nrad upvp(i,j) = (urp(i,j) - urpb(i))* + (vtp(i,j) - vtpb(i)) 20 continue 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 30 i=2,nrad rrefc(i) = cf*(2.*upvpb(i)/(1000.*radk(i)) + dupvpb(i)) 30 continue c refc = 0.0 icount = 0 do 40 i=2,nrad if (radk(i) .gt. radmax) go to 2040 refc = refc + rrefc(i) icount = icount + 1 40 continue c 2040 continue if (icount .ge. 1) then refc = refc/float(icount) else refc = 0.0 endif c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 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) go to 2040 pefc = pefc + t1r(i) icount = icount + 1 enddo c 2040 continue if (icount .ge. 1) then pefc = pefc/float(icount) else pefc = 0.0 endif c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 i=3,nrad if (radk(i) .gt. radtem) then irad = i-1 go to 1000 endif 10 continue c 1000 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 i=3,nrad if (radk(i) .gt. radtem) then irad = i-1 go to 1000 endif 10 continue c 1000 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 i=1,nx defxy(i,j) = sqrt( (dudy(i,j)+dvdx(i,j))**2 + + (dudx(i,j)-dvdy(i,j))**2 ) 10 continue c ierr = 0 c return end 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 i=2,nx-1 dfdx(i,j) = (f(i+1,j) - f(i-1,j))/tdx 10 continue c c Calculate one-sided differences do 20 j=1,ny dfdx( 1,j) = (f( 2,j) - f( 1,j))/dx dfdx(nx,j) = (f(nx,j) - f(nx-1,j))/dx 20 continue c return end 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 i=1,nx do 10 j=2,ny-1 dfdy(i,j) = (f(i,j+1) - f(i,j-1))/tdy 10 continue c c Calculate one-sided differences do 20 i=1,nx dfdy(i, 1) = (f(i, 2) - f(i, 1))/dy dfdy(i,ny) = (f(i,ny) - f(i,ny-1))/dy 20 continue c return end 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 10 n=1,nvp ubar = ubar + wvp(n)*uvp(n) vbar = vbar + wvp(n)*vvp(n) 10 continue c c Calculate generarlized shear shrg = 0.0 do 20 n=1,nvp shrg = shrg + wvp(n)*sqrt( (uvp(n)-ubar)**2 + + (vvp(n)-vbar)**2 ) 20 continue 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 30 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 30 continue c return end 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 10 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 10 continue c return end 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 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 10 k=1,ipmx if (plevr(k) .eq. ps1) ips1=k if (plevr(k) .eq. ps2) ips2=k 10 continue 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 20 j=1,nlat1 do 20 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 20 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 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 10 j=1,ny do 10 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 10 continue c do 15 j=1,ny do 15 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 15 continue 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 20 j=1,ntheta thetad(j) = dtheta*float(j-1) 20 continue c c Specify radial interval in km drad = 100.0 c c Find maximum radius value on Cartesian grid radmax = -1000.0 do 25 j=1,ny do 25 i=1,nx if (rgk(i,j) .gt. radmax) radmax=rgk(i,j) 25 continue 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 30 i=1,nrad radk(i) = drad*float(i-1) 30 continue c c Calculate x,y values at polar grid points do 35 j=1,ntheta do 35 i=1,nrad xra(i,j) = radk(i)*cos(dtr*thetad(j)) yra(i,j) = radk(i)*sin(dtr*thetad(j)) 35 continue c c Calculate polar coordinate mask do 40 j=1,ntheta do 40 i=1,nrad xtem = xra(i,j) ytem = yra(i,j) pmask(i,j) = 1.0 if (xtem .ge. xmax .or. xtem .le. xmin) pmask(i,j)=0.0 if (ytem .ge. ymax .or. ytem .le. ymin) pmask(i,j)=0.0 40 continue c return END 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 c xtem = dtk*(tlon-olon)*colat ytem = dtk*(tlat-olat) c return END 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 c tlon = olon + xtem/(dtk*colat) tlat = olat + ytem/(dtk) c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 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) 10 continue c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 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) 10 continue c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 n=1,nvar if (indexvn(n) .eq. vtype) then idv = n go to 1000 endif 10 continue c write(lulog,900) vtype 900 format(/,' Unrecognized variable type: ',a1,' in routine varget.') stop 1000 continue 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 15 k=1,ipmx if (plevr(k) .eq. press) then idp = k if (indexp(k) .eq. 0) go to 1495 go to 1500 endif 1495 continue 15 continue c 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 c 1500 continue iexist=1 c c Extract variable if (idv .eq. 1) then do 20 j=1,nlat1 do 20 i=1,nlon1 fxy(i,j) = u(i,j,idp) 20 continue elseif (idv .eq. 2) then do 25 j=1,nlat1 do 25 i=1,nlon1 fxy(i,j) = v(i,j,idp) 25 continue elseif (idv .eq. 3) then do 30 j=1,nlat1 do 30 i=1,nlon1 fxy(i,j) = z(i,j,idp) 30 continue elseif (idv .eq. 4) then do 35 j=1,nlat1 do 35 i=1,nlon1 fxy(i,j) = t(i,j,idp) 35 continue elseif (idv .eq. 5) then do 40 j=1,nlat1 do 40 i=1,nlon1 fxy(i,j) = r(i,j,idp) 40 continue else write(lulog,908) idv 908 format(/,' Illegal variable number ',i2,' in routine varget.') stop endif c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ntheta do 10 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 10 continue c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 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 10 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 i=1,nx if (rgk(i,j) .le. radtem) then fxybar = fxybar + fxy(i,j) nbar = nbar + 1 endif 10 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 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 10 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 i=1,nrad jcount = 0 ftemp = 0.0 do 20 j=1,ntheta if (pmask(i,j) .gt. 0.5) then ftemp = ftemp + fra(i,j) jcount = jcount + 1 endif 20 continue if (jcount .gt. 0) then fr(i) = ftemp/float(jcount) else fr(i) = 0.0 endif 10 continue c if (iorz .eq. 0) fr(1)=0.0 c return END 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 j=1,ny do 10 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 10 continue c return END subroutine ddrcal(fr,dfdr) c This routine calculates the radial derivative of fr c using centered derivatives. c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 i=3,nrad-1 dfdr(i) = (fr(i+1)-fr(i-1))/drm2 10 continue 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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 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 10 continue c write(lulog,300) 300 format(/,' P U V Z T RH') do 15 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)) 15 continue c return END 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 (+/-xradk,+/-yradk) 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 10 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 10 continue c write(lulog,300) 300 format(/,' P U V Z T RH') do 15 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)) 15 continue c return END subroutine tadd(iyr,imon,iday,itime,ihra,iyra,imona,idaya,itimea) c This routine calculates the year,month,day,time c (iyra,imona,idaya,itimea) that are ihra hours after c the input year,month,day,time (iyr,imon,iday,itime). c c ihra can not exceed the number of hours in 9 months. c c If ihra is negative, the year,month,day,time that are abs(ihra) c before the input year,month,day,time are calculated. c c The Julian day utilities jday and jdayi are called c by this routine. c iyra = iyr imona = imon idaya = iday c if (ihra .lt. 0) go to 1000 c c ** Start calculation for positive ihra ** c c Add the hours to the input time itimea = itime+ihra c if (itimea .lt. 24) return c c Calculate the number of extra days in the itimea variable c and subtract them from itimea idayex = itimea/24 itimea = itimea - idayex*24 c c Calculate Julian day of input date c and add extra days to it call jday(imon,iday,iyr,jdate) jdate = jdate + idayex c c Check to see if year has changed call jday(12,31,iyr,jmax) c if (jdate .gt. jmax) then iyra = iyr+1 jdate = jdate-jmax else iyra = iyr endif c c Calculate month and day corresponding to the c increased jdate call jdayi(jdate,iyra,imona,idaya) c return c 1000 continue c ** Start calculation for negative ihra ** c Add the hours to the input time itimea = itime+ihra c if (itimea .ge. 0) return c c Calculate adjusted time idayex = -1 + itimea/24 itimea = itimea - idayex*24 if (itimea .eq. 24) then itimea=0 idayex = idayex + 1 endif c c Calculate Julian day of input date c and subtract extra days from it call jday(imon,iday,iyr,jdate) jdate = jdate + idayex c c Check to see if year has changed if (jdate .lt. 1) then iyra = iyra-1 jdate = jdate + 365 if (mod(iyr-1,4) .eq. 0) jdate = jdate + 1 endif c c Calculate month and day corresponding to the c decreased jdate call jdayi(jdate,iyra,imona,idaya) c return end subroutine jday(imon,iday,iyear,julday) c This routine calculates the Julian day (julday) from c the month (imon), day (iday), and year (iyear). The c appropriate correction is made for leap year. c dimension ndmon(12) c c Specify the number of days in each month ndmon(1) = 31 ndmon(2) = 28 ndmon(3) = 31 ndmon(4) = 30 ndmon(5) = 31 ndmon(6) = 30 ndmon(7) = 31 ndmon(8) = 31 ndmon(9) = 30 ndmon(10) = 31 ndmon(11) = 30 ndmon(12) = 31 c c Correct for leap year if (mod(iyear,4) .eq. 0) ndmon(2)=29 c c Check for illegal input if (imon .lt. 1 .or. imon .gt. 12) then julday=-1 return endif c if (iday .lt. 1 .or. iday .gt. ndmon(imon)) then julday=-1 return endif c c Calculate the Julian day julday = iday if (imon .gt. 1) then do 10 i=2,imon julday = julday + ndmon(i-1) 10 continue endif c return end subroutine jdayi(julday,iyear,imon,iday) c This routine calculates the month (imon) and day (iday) c from the Julian day (julday) and year (iyear). c The appropriate correction is made for leap year. c dimension ndmon(12),nsum(13) c c Specify the number of days in each month ndmon(1) = 31 ndmon(2) = 28 ndmon(3) = 31 ndmon(4) = 30 ndmon(5) = 31 ndmon(6) = 30 ndmon(7) = 31 ndmon(8) = 31 ndmon(9) = 30 ndmon(10) = 31 ndmon(11) = 30 ndmon(12) = 31 c c Correct for leap year if (mod(iyear,4) .eq. 0) ndmon(2)=29 c c Check for illegal input if (mod(iyear,4) .eq. 0) then mxjul = 366 else mxjul = 365 endif c if (julday .lt. 1 .or. julday .gt. mxjul) then imon = -1 iday = -1 return endif c c Calculate the month and day nsum(1) = 0 do 10 i=1,12 nsum(i+1) = nsum(i) + ndmon(i) 10 continue c do 20 i=2,13 if (julday .le. nsum(i)) then imon = i-1 go to 1000 endif 20 continue 1000 continue c iday = julday - nsum(imon) c return end subroutine fcopy(fxy,fxyc) c This routine copies the 2-D array fxy to fxyc c C.DUP lsdiag,primd parameter(ixmx=121,iymx=91,ipmx=21) parameter(irmx=300,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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 99 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 99 continue c return end 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=121,iymx=91,ipmx=21) parameter(irmx=300,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 first guess lat/ c 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 ++ Start iteration do 99 n=1,niter c do jj=-nsearch,nsearch do ii=-nsearch,nsearch c c Calculate x,y coordinates of center-relative r,theta grid poi 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) 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 99 continue 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 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 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 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 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 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) go to 1000 count = count + 1.0 fbar = fbar + fr(i) enddo c 1000 continue c if (count .ge. 1.0) then fbar = fbar/count else fbar = 0.0 endif c return end subroutine vremove(uxy,vxy,rcf,thetacf,xcen,ycen, + mrcf,mtcf,nrcf,ntcf,uraa,vtaa,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 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=121,iymx=91,ipmx=21) parameter(irmx=300,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) dimension uraa(0:mrcf),vtaa(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 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 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 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 Azimuthally average ur,vt call azavgcf(urrt,urr,mrcf,mtcf,nrcf,ntcf) call azavgcf(vtrt,vtr,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 kstart = ifix(rminvr/dr) kend = nrcf c do k=kstart+1,nrcf if (vtr(k) .lt. 0.0) then c write(6,*) 'k,vtr(k)',k,vtr(k) kend = k-1 go to 1000 endif enddo 1000 continue 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 subtracting 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) go to 2000 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-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 2000 continue c enddo enddo c return end 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 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 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 c ** cfindxy code marker (end) ** subroutine splcal(plev,u,v,cx,cy,alpha,ml, + dw,w,ubard,vbard,pbard,ubar,vbar,pbar) c This routine calculates the weights w to determine the vertically c averaged horizontal wind that is as close as possible to c the storm motion cx,cy. The weights are used to determine the c pressure of the steering level. c c Input: c plev - 1-D array containing the pressure levels (mb) c u,v - The enviromental wind at the levels plev c cx,cy - The components of the storm motion vector c alpha - The coefficient for contraining the steering c weights so that they are not "too far" from c the mass weighted average weights c ml - The number of pressure levels c c Output: c dw - The weights for a mass-weighted average c w - The weights for the optimal steering c ubard,vbard - The mass-weighted average horizontal wind c ubar,vbar - The optimally weighted horizontal mean c pbard - The mass-weighted pressure c pbar - The optimally weighted pressure c dimension plev(ml),u(ml),v(ml),dw(ml),w(ml) c c Local variables parameter (np=12,mp=2) c dimension a(np,np),b(np,mp) c c Set ipen=1 for penalty in terms of (W-M) or c ipen=2 for penalty in terms of (W/M-1) ipen=2 c c Make sure np is large enough if (ml .gt. np) then write(6,100) 100 format(/,' np too small in routine wcal') stop endif c n=ml-1 c c Calculate deep-layer mean weights call dlmw(plev,dw,ml) c c Calculate column matrix on the RHS of the linear system for w uk = u(ml) vk = v(ml) dk= dw(ml) c if (ipen .eq. 1) then do 10 k=1,n b(k,1) = (cx-uk)*(u(k)-uk) + (cy-vk)*(v(k)-vk) + + alpha*(1.0 + dw(k) - dk) 10 continue else do 11 k=1,n b(k,1) = (cx-uk)*(u(k)-uk) + (cy-vk)*(v(k)-vk) + + alpha*(1.0/(dk*dk) + 1.0/dw(k) - 1.0/dk) 11 continue endif c c Calculate w coefficient matrix if (ipen .eq. 1) then do 15 j=1,n do 15 i=1,n if (i .eq. j) then ac = 2.0 else ac = 1.0 endif c a(i,j) = (u(i)-uk)*(u(j)-uk) + + (v(i)-vk)*(v(j)-vk) + ac*alpha 15 continue else do 16 j=1,n do 16 i=1,n if (i .eq. j) then ac = (1.0/dk)**2 + (1.0/dw(j))**2 else ac = (1.0/dk)**2 endif c a(i,j) = (u(i)-uk)*(u(j)-uk) + + (v(i)-vk)*(v(j)-vk) + ac*alpha 16 continue endif c c Calculate optimal weights call gaussj(a,n,np,b,1,mp) c do 30 i=1,n w(i) = b(i,1) 30 continue c w(ml) = 1.0 do 40 i=1,n w(ml) = w(ml) - w(i) 40 continue c c Calculate vertically weighted variables ubard = 0.0 vbard = 0.0 ubar = 0.0 vbar = 0.0 pbard = 0.0 pbar = 0.0 do 50 k=1,ml ubard = ubard + dw(k)*u(k) vbard = vbard + dw(k)*v(k) pbard = pbard + dw(k)*plev(k) ubar = ubar + w(k)*u(k) vbar = vbar + w(k)*v(k) pbar = pbar + w(k)*plev(k) 50 continue c return end subroutine dlmw(plev,dw,ml) dimension plev(ml),dw(ml) c if (ml .eq. 1) then dw(1) = 1.0 return endif c pdeep = plev(ml) - plev(1) c dw( 1) = 0.5*(plev(2)-plev(1))/pdeep dw(ml) = 0.5*(plev(ml)-plev(ml-1))/pdeep c if (ml .eq. 2) return c do 10 k=2,ml-1 dw(k) = 0.5*(plev(k+1)-plev(k-1))/pdeep 10 continue c return end SUBROUTINE GAUSSJ(A,N,NP,B,M,MP) PARAMETER (NMAX=50) DIMENSION A(NP,NP),B(NP,MP),IPIV(NMAX),INDXR(NMAX),INDXC(NMAX) DO 11 J=1,N IPIV(J)=0 11 CONTINUE DO 22 I=1,N BIG=0. DO 13 J=1,N IF(IPIV(J).NE.1)THEN DO 12 K=1,N IF (IPIV(K).EQ.0) THEN IF (ABS(A(J,K)).GE.BIG)THEN BIG=ABS(A(J,K)) IROW=J ICOL=K ENDIF ELSE IF (IPIV(K).GT.1) THEN PAUSE 'Singular matrix' ENDIF 12 CONTINUE ENDIF 13 CONTINUE IPIV(ICOL)=IPIV(ICOL)+1 IF (IROW.NE.ICOL) THEN DO 14 L=1,N DUM=A(IROW,L) A(IROW,L)=A(ICOL,L) A(ICOL,L)=DUM 14 CONTINUE DO 15 L=1,M DUM=B(IROW,L) B(IROW,L)=B(ICOL,L) B(ICOL,L)=DUM 15 CONTINUE ENDIF INDXR(I)=IROW INDXC(I)=ICOL IF (A(ICOL,ICOL).EQ.0.) PAUSE 'Singular matrix.' PIVINV=1./A(ICOL,ICOL) A(ICOL,ICOL)=1. DO 16 L=1,N A(ICOL,L)=A(ICOL,L)*PIVINV 16 CONTINUE DO 17 L=1,M B(ICOL,L)=B(ICOL,L)*PIVINV 17 CONTINUE DO 21 LL=1,N IF(LL.NE.ICOL)THEN DUM=A(LL,ICOL) A(LL,ICOL)=0. DO 18 L=1,N A(LL,L)=A(LL,L)-A(ICOL,L)*DUM 18 CONTINUE DO 19 L=1,M B(LL,L)=B(LL,L)-B(ICOL,L)*DUM 19 CONTINUE ENDIF 21 CONTINUE 22 CONTINUE DO 24 L=N,1,-1 IF(INDXR(L).NE.INDXC(L))THEN DO 23 K=1,N DUM=A(K,INDXR(L)) A(K,INDXR(L))=A(K,INDXC(L)) A(K,INDXC(L))=DUM 23 CONTINUE ENDIF 24 CONTINUE RETURN END SUBROUTINE PCMIN(SST,PSL,P,T,R,NA,N,PMIN,VMAX,IFL) C C Revised on 9/24/2005 to fix convergence problems at high pressure C C *** This subroutine calculates the maximum wind speed *** C *** and mimimum central pressure *** C *** achievable in tropical cyclones, given a sounding *** C *** and a sea surface temperature. *** C C INPUT: SST: Sea surface temperature in C C C PSL: Sea level pressure (mb) C C P,T,R: One-dimensional arrays of dimension NA C containing pressure (mb), temperature (C), C and mixing ratio (g/kg). The arrays MUST be C arranged so that the lowest index corresponds C to the lowest model level, with increasing index C corresponding to decreasing pressure. The temperature C sounding should extend to at least the tropopause and C preferably to the lower stratosphere, however the C mixing ratios are not important above the boundary C layer. Missing mixing ratios can be replaced by zeros. C C NA: The dimension of P,T and R C C N: The actual number of points in the sounding C (N is less than or equal to NA) C C OUTPUT: PMIN is the minimum central pressure, in mb C C VMAX is the maximum surface wind speed, in m/s C (reduced to reflect surface drag) C C IFL is a flag: A value of 1 means OK; a value of 0 C indicates no convergence (hypercane); a value of 2 C means that the CAPE routine failed. C C----------------------------------------------------------------------------- REAL T(NA), P(NA), R(NA) C C *** Adjustable constant: Ratio of C_k to C_D *** C CKCD=0.9 C C *** Adjustable constant for buoyancy of displaced parcels: *** C *** 0=Reversible ascent; 1=Pseudo-adiabatic ascent *** C SIG=0.0 C C *** Adjustable switch: if IDISS = 0, no dissipative heating is *** C *** allowed; otherwise, it is *** C IDISS=1 C C *** Exponent, b, in assumed profile of azimuthal velocity in eye, *** C *** V=V_m(r/r_m)^b. Used only in calculation of central pressure *** C b=2.0 C C *** Set level from which parcels lifted *** C NK=1 C C *** Factor to reduce gradient wind to 10 m wind C VREDUC=0.8 C C *** Normalize certain quantities *** C SSTK=SST+273.15 ES0=6.112*EXP(17.67*SST/(243.5+SST)) DO 40 I=1,N R(I)=R(I)*0.001 T(I)=T(I)+273.15 40 CONTINUE C C *** Default values *** C VMAX=0.0 PMIN=PSL IFL=1 C NP=0 PM=950.0 C C *** Find environmental CAPE *** C TP=T(NK) RP=R(NK) PP=P(NK) CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEA,TOA,IFLAG) IF(IFLAG.NE.1)IFL=2 C C *** Begin iteration to find mimimum pressure *** C 100 CONTINUE C C *** Find CAPE at radius of maximum winds *** C TP=T(NK) PP=MIN(PM,1000.0) RP=0.622*R(NK)*PSL/(PP*(0.622+R(NK))-R(NK)*PSL) CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEM,TOM,IFLAG) IF(IFLAG.NE.1)IFL=2 RAT=SSTK/TOM IF(IDISS.EQ.0)RAT=1.0 C C *** Find saturation CAPE at radius of maximum winds *** C TP=SSTK PP=MIN(PM,1000.0) RP=0.622*ES0/(PP-ES0) CALL CAPE(TP,RP,PP,T,R,P,NA,N,SIG,CAPEMS,TOMS,IFLAG) IF(IFLAG.NE.1)IFL=2 C C *** Initial estimate of minimum pressure *** C RS0=RP TV1=T(1)*(1.+R(1)/0.622)/(1.+R(1)) TVAV=0.5*(TV1+SSTK*(1.+RS0/0.622)/(1.+RS0)) C CAT=0.5*CKCD*RAT*(CAPEMS-CAPEM) CAT=CAPEM-CAPEA+0.5*CKCD*RAT*(CAPEMS-CAPEM) CAT=MAX(CAT,0.0) PNEW=PSL*EXP(-CAT/(287.04*TVAV)) C C *** Test for convergence *** C IF(ABS(PNEW-PM).GT.0.2)THEN PM=PNEW NP=NP+1 IF(NP.GT.1000.OR.PM.LT.400.0)THEN PMIN=PSL IFL=0 GOTO 900 END IF GOTO 100 ELSE CATFAC=0.5*(1.+1./b) C CAT=CKCD*RAT*CATFAC*(CAPEMS-CAPEM) CAT=CAPEM-CAPEA+CKCD*RAT*CATFAC*(CAPEMS-CAPEM) CAT=MAX(CAT,0.0) PMIN=PSL*EXP(-CAT/(287.04*TVAV)) END IF 900 CONTINUE FAC=MAX(0.0,(CAPEMS-CAPEM)) VMAX=VREDUC*SQRT(CKCD*RAT*FAC) C C *** Renormalize sounding arrays *** C DO 910 I=1,N R(I)=R(I)*1000.0 T(I)=T(I)-273.15 910 CONTINUE C RETURN END C SUBROUTINE CAPE(TP,RP,PP,T,R,P,ND,N,SIG,CAPED,TOB,IFLAG) C C This subroutine calculates the CAPE of a parcel with pressure PP (mb), C temperature TP (K) and mixing ratio RP (gm/gm) given a sounding C of temperature (T in K) and mixing ratio (R in gm/gm) as a function C of pressure (P in mb). ND is the dimension of the arrays T,R and P, C while N is the actual number of points in the sounding. CAPED is C the calculated value of CAPE and TOB is the temperature at the C level of neutral buoyancy. IFLAG is a flag C integer. If IFLAG = 1, routine is successful; if it is 0, routine did C not run owing to improper sounding (e.g.no water vapor at parcel level). C IFLAG=2 indicates that routine did not converge. C REAL T(ND),R(ND),P(ND),TVRDIF(100) REAL NA C C *** Default values *** C CAPED=0.0 TOB=T(1) IFLAG=1 C C *** Check that sounding is suitable *** C IF(RP.LT.1.0E-6.OR.TP.LT.200.0)THEN IFLAG=0 RETURN END IF C C *** Assign values of thermodynamic constants *** C CPD=1005.7 CPV=1870.0 C CL=4190.0 CL=2500.0 CPVMCL=CPV-CL RV=461.5 RD=287.04 EPS=RD/RV ALV0=2.501E6 C C *** Define various parcel quantities, including reversible *** C *** entropy, S. *** C TPC=TP-273.15 ESP=6.112*EXP(17.67*TPC/(243.5+TPC)) EVP=RP*PP/(EPS+RP) RH=EVP/ESP RH=MIN(RH,1.0) ALV=ALV0-CPVMCL*TPC S=(CPD+RP*CL)*LOG(TP)-RD*LOG(PP-EVP)+ 1 ALV*RP/TP-RP*RV*LOG(RH) C C *** Find lifted condensation pressure, PLCL *** C CHI=TP/(1669.0-122.0*RH-TP) PLCL=PP*(RH**CHI) C C *** Begin updraft loop *** C NCMAX=0 DO J=1,N TVRDIF(J)=0.0 END DO C JMIN=1E6 DO 200 J=1,N C C *** Don't bother lifting parcel above 60 mb and skip sections of sounding below parcel level *** C IF(P(J).LT.59.0.OR.P(J).GE.PP)GOTO 200 C JMIN=MIN(JMIN,J) C C *** Parcel quantities below lifted condensation level *** C IF(P(J).GE.PLCL)THEN TG=TP*(P(J)/PP)**(RD/CPD) RG=RP C C *** Calculate buoyancy *** C TLVR=TG*(1.+RG/EPS)/(1.+RG) TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J)) ELSE C C *** Parcel quantities above lifted condensation level *** C TG=T(J) TJC=T(J)-273.15 ES=6.112*EXP(17.67*TJC/(243.5+TJC)) RG=EPS*ES/(P(J)-ES) C C *** Iteratively calculate lifted parcel temperature and mixing *** C *** ratio for reversible ascent *** C NC=0 120 CONTINUE NC=NC+1 C C *** Calculate estimates of the rates of change of the entropy *** C *** with temperature at constant pressure *** C ALV=ALV0-CPVMCL*(TG-273.15) SL=(CPD+RP*CL+ALV*ALV*RG/(RV*TG*TG))/TG EM=RG*P(J)/(EPS+RG) SG=(CPD+RP*CL)*LOG(TG)-RD*LOG(P(J)-EM)+ 1 ALV*RG/TG IF(NC.LT.3)THEN AP=0.3 ELSE AP=1.0 END IF TGNEW=TG+AP*(S-SG)/SL C C *** Test for convergence *** C IF(ABS(TGNEW-TG).GT.0.001)THEN TG=TGNEW TC=TG-273.15 ENEW=6.112*EXP(17.67*TC/(243.5+TC)) C C *** Bail out if things get out of hand *** C IF(NC.GT.500.OR.ENEW.GT.(P(J)-1.0))THEN IFLAG=2 RETURN END IF RG=EPS*ENEW/(P(J)-ENEW) GOTO 120 END IF NCMAX=MAX(NC,NCMAX) C C *** Calculate buoyancy *** C RMEAN=SIG*RG+(1.-SIG)*RP TLVR=TG*(1.+RG/EPS)/(1.+RMEAN) TVRDIF(J)=TLVR-T(J)*(1.+R(J)/EPS)/(1.+R(J)) END IF 200 CONTINUE C C *** Begin loop to find NA, PA, and CAPE from reversible ascent *** C NA=0.0 PA=0.0 C C *** Find maximum level of positive buoyancy, INB *** C INB=1 DO 550 J=N,JMIN,-1 IF(TVRDIF(J).GT.0.0)INB=MAX(INB,J) 550 CONTINUE IF(INB.EQ.1)RETURN C C *** Find positive and negative areas and CAPE *** C IF(INB.GT.1)THEN DO 600 J=JMIN+1,INB PFAC=RD*(TVRDIF(J)+TVRDIF(J-1))*(P(J-1)-P(J))/(P(J)+P(J-1)) PA=PA+MAX(PFAC,0.0) NA=NA-MIN(PFAC,0.0) 600 CONTINUE C C *** Find area between parcel pressure and first level above it *** C PMA=(PP+P(JMIN)) PFAC=RD*(PP-P(JMIN))/PMA PA=PA+PFAC*MAX(TVRDIF(JMIN),0.0) NA=NA-PFAC*MIN(TVRDIF(JMIN),0.0) C C *** Find residual positive area above INB and TO *** C PAT=0.0 TOB=T(INB) IF(INB.LT.N)THEN PINB=(P(INB+1)*TVRDIF(INB)-P(INB)*TVRDIF(INB+1))/ 1 (TVRDIF(INB)-TVRDIF(INB+1)) PAT=RD*TVRDIF(INB)*(P(INB)-PINB)/(P(INB)+PINB) TOB=(T(INB)*(PINB-P(INB+1))+T(INB+1)*(P(INB)-PINB))/ 1 (P(INB)-P(INB+1)) END IF C C *** Find CAPE *** C CAPED=PA+PAT-NA CAPED=MAX(CAPED,0.0) END IF C RETURN END subroutine virt(t,p,rh,tv) c This routine calculates the virtual temperature c given the temperature, pressure and relative humidity c input: t = temperature (C) c p = pressure (mb) c rh = relative humidity (%) c output: tv = virtual temperature (C) c c Check input if (t .lt. -273.15 .or. p .gt. 1050. .or. p .lt. 10.0 + .or. rh .gt. 100.0 .or. rh .lt. 0.) then tv = 0.0 return endif c c Specify constants tcon = 273.15 epsil = 0.622 c c Convert input to mks units tk = t + tcon pp = 100.0*p rhr = rh/100.0 c c Calculate saturation vapor pressure call teten(tk,es) c c Calculate vapor pressure ws = es*epsil/(pp-es) w = rhr*ws e = w*pp/(w+epsil) c c Calculate virtual temperature and convert back to Celsius tv = tk/( 1.0 -(e/pp)*(1.0-epsil) ) tv = tv - tcon c return end subroutine virti(tv,p,rh,t) c This routine calculates the temperature c given the virtual temperature, pressure and relative humidity. c A simple iteration method is used for the calculation. c input: tv = virtual temperature (C) c p = pressure (mb) c rh = relative humidity (%) c output: t = virtual temperature (C) c c Check input if (tv .lt. -273.15 .or. p .gt. 1050. .or. p .lt. 10.0 + .or. rh .gt. 100.0 .or. rh .lt. 0.) then tv = 0.0 return endif c c Specify constants tcon = 273.15 epsil = 0.622 c c Convert input to mks units tvk = tv + tcon pp = 100.0*p rhr = rh/100.0 c c Make first guess for the temperature tk = tvk c c Iterate to fine the actual temperature nit = 5 do 10 i=1,nit c Calculate saturation vapor pressure call teten(tk,es) c c Calculate vapor pressure ws = es*epsil/(pp-es) w = rhr*ws e = w*pp/(w+epsil) c c Calculate temperature tk = tvk*( 1.0 -(e/pp)*(1.0-epsil) ) c write(6,*) ' ',i,tk 10 continue c c Convert temperature to Celsius t = tk - tcon c return end subroutine teten(tk,es) c This routine calculates the saturation vapor pressure using c Teten's formula c input: tk = temperature (K) c output: es = saturation vapor pressure (Pa) c c Check for valid temperature if (tk .lt. 40.0) then es = 0.0 return endif c c1 = 610.78 c2 = 17.269 c3 = 35.86 c es = c1*exp( c2*(tk-273.16)/(tk-35.86) ) c return end subroutine parsub(pb,tb,rhb,dp,pmin,iadj,mxd,psat,pp,tp,tvp,npp) c This program calculates temperature/pressure profiles c for parcels lifted under various assumptions. The thermodynamic c formulation described by Ooyama (1990) is used for c the calculations. The ice phase is not included in the calculations. c c Input: c pb = boundary layer pressure (Pa) c tb = boundary layer temperature (K) c rhb = boundary layer relative humidity (%) c dp = Pressure increment for profile (Pa) c pmin = Minimum pressure of profile (Pa) c iadj = 0 to carry liquid water c = 1 to remove liquid water as it condenses c mxd = Maximum dimension of the pp,tp,tvp arrays c c Output: c psat = Pressure (Pa) where parcel first becomes saturated c pp = Pressure profile (Pa) c tp = Temperature profile (K) c tvp = Virtual temperature profile c npp = No. of points in pressure profile c dimension pp(mxd),tp(mxd),tvp(mxd) c common /gc/ ra,rv,cpa,cpv,cva,cvv,cck common /rs/ t0,p0,z0,e0,etas0,cl0 common /sc/ a,d0,d1,d2 c c Specify number of iterations for the pressure loop nit = 10 c c Specify gas constants and conversion from C to K (mks units) ra = 287.0 rv = 462.5 cpa = 1005.7 cpv = 1875.0 cva = cpa - ra cvv = cpv - rv cck = 273.15 eps = ra/rv c c Specify constants for saturation vapor pressure equation a = 611.2 d0 = 17.67 d1 = 273.15 d2 = 29.65 c c Specify reference state temperature and dry air pressure t0 = 273.15 p0 = 1000.0e+2 c c Calculate dry air density of the reference state z0 = p0/(ra*t0) c c Calculate saturation vapor pressure and density of reference state call esat(t0,e0,de0,d2e0) etas0 = e0/(rv*t0) c c Calculate lambda of the reference state cl0 = rv*t0*de0/e0 c c Set initial value of condensate to zero etac = 0.0 c c Calculate saturation vapor pressure and saturation c vapor density of parcel call esat(tb,eb,deb,d2eb) etasb = eb/(rv*tb) c c Calculate actual vapor density and vapor pressure of parcel etavb = rhb*etasb/100.0 pvb = etavb*rv*tb c c Calculate dry air pressure and density of the parcel pab = pb - pvb zb = pab/(ra*tb) c c Calculate initial mixing ratio and saturation mixing ratio wb = etavb/zb wsb = etasb/zb c c Calculate initial condensate mixing ratio wc = etac/zb c c Calculate dry air entropy call sacal(tb,zb,sab) c c Calculate moisture entropy (saturated and unsaturated) call smcal(tb,etavb,sm1b,sm2b) c c Calculate entropy density sg1b = zb*sab + etavb*sm1b sg2b = zb*sab + etavb*sm2b c sb = sg1b/zb c c Start the pressure loop c c Initialize variables with those from the boundary layer c before starting the pressure loop pv = pvb t = tb s = sb w = wb c pp(1) = pb tp(1) = tb tvp(1) = tb*(1.0 + wb/eps)/(1.0 + wc + wb) c psat = -100.0 i = 2 c c Start pressure loop with pressure to the nearest 10 mb. p = 1000.0*float(ifix( (pb-1.0)/1000.0 ) ) c 2000 continue c pp(i) = p c c Start iteration n = 0 c 1000 continue pa = p - pv z = pa/(ra*t) eta = w*z c call tcal(t,z,eta,s,t1,t2) t = t1 if (t2 .gt. t1) t=t2 c call esat(t,e,des,d2es) etas = e/(rv*t) etac = eta - etas if (etac .lt. 0.0) etac=0.0 etav = eta - etac c ws = etas/z wc = etac/z c wv = (w - wc) c pv = etav*rv*t rh = etav/etas c n = n+1 if (n .lt. nit) go to 1000 c c Check for saturation pressure if (psat .lt. 0.0 .and. rh .gt. .999) psat = p c if (iadj .eq. 1) then c Remove condensate as it forms c and adjust mixing ratio and entropy density accordingly if (etac .gt. 0.0) then etac = 0.0 eta = etas + etac w = eta/z c call sacal(t,z,sa) call smcal(t,eta,sm1,sm2) c c Calculate entropy density sg = z*sa + eta*sm2 s = sg/z endif endif c c Calculate virtual temperature tv = t*(1.0 + wv/eps)/(1.0 + wc + wv) c tp(i) = t tvp(i) = tv c i = i+1 p = p - dp if (p .ge. pmin .and. i .le. mxd) go to 2000 c npp = i-1 c return end subroutine esat(t,e,de,d2e) c This routine calculates the saturation vapor pressure c (e) in Pa and its first and second derivatives, given c the temperature (T) in K using eqn. (10) from Bolton c (1980), except that T is input in K rather than C. c common /sc/ a,d0,d1,d2 c if (t .le. d2) then e = 0.0 de = 0.0 d2e = 0.0 else e = a*exp(d0*(t-d1)/(t-d2)) de = e*d0*(d1-d2)/( (t-d2)*(t-d2) ) d2e = de*de/e -2.0*de/(t-d2) endif c return end subroutine sacal(t,z,sa) c This routine calculates the dry air entropy (sa) given c the temperature (K) and dry air density z (kg/m3). c common /gc/ ra,rv,cpa,cpv,cva,cvv,cck common /rs/ t0,p0,z0,e0,etas0,cl0 c sa = cva*alog(t/t0) - ra*alog(z/z0) c return end subroutine smcal(t,eta,sm1,sm2) c This routine calculates the entropy of the water substance c (sm1 and sm2) given the temperature (K) and moisture density c (kg/m3). sm1 is valid for unsaturated air and sm2 is valid c for saturated air. c common /gc/ ra,rv,cpa,cpv,cva,cvv,cck common /rs/ t0,p0,z0,e0,etas0,cl0 c sm1 = cvv*alog(t/t0) - rv*alog(eta/etas0) + cl0 c call esat(t,e,de,d2e) etas = e/(rv*t) c sm2 = cvv*alog(t/t0) - rv*t*de/e + cl0 - rv*alog(etas/etas0) + + de/eta c return end subroutine tcal(t,z,eta,s,t1,t2) c This routine inverts the definition of the entropy density c using Newton-Raphson iteration. c c The necessary input is as follows: c t = first guess temperature c z = dry air density c eta = mass density of total moisture c s = entropy density c c There are two temperature outputs. The greater of the two c is the proper temperature. If t2 > t1 then the air is c saturated. If t1 > t2 then the air is not saturated. c common /gc/ ra,rv,cpa,cpv,cva,cvv,cck c c Specify the number of iterations nit = 10 c c Initialize t1 and t2 and counter t1 = t t2 = t n = 0 c c Start iteration 1000 continue c c Calculations for t1 call sacal(t1,z,sa) call smcal(t1,eta,sm1,sm2) c f = sa + eta*sm1/z - s c dsadt1 = cva/t1 dsmdt1 = cvv/t1 c dfdt1 = dsadt1 + eta*dsmdt1/z c t1 = t1 - f/dfdt1 c c Calculations for t2 call sacal(t2,z,sa) call smcal(t2,eta,sm1,sm2) c f = sa + eta*sm2/z - s c call esat(t2,e,de,d2e) c dsadt2 = cva/t2 dsmdt2 = cpv/t2 - (rv*de/e)*(2.0 - t2*de/e) + + d2e*(1.0/eta - rv*t2/e) c dfdt2 = dsadt2 + eta*dsmdt2/z c t2 = t2 - f/dfdt2 c n = n + 1 c c write(6,*) ' n,t1,t2',n,t1,t2 if (n .lt. nit) go to 1000 c return end subroutine pstcal(z1,z2,t1,t2,p1,ps,ts) c This routine calculates the surface pressure (psfc) c from thermodynamic variables near the surface. c c level 1 = level closest to the surface (usually 1000 mb) c level 2 = next level up (usually 850 or 925 mb) c c input: p1 = pressure (Pa) c zi = geopotential height (m) (i=1,2) c ti = temperature (K) (i=1,2) c c output: ps = surface pressure (Pa) c ts = surface temperature (K) c c Specify physical constants (mks units) rd = 287.0 g = 9.81 c c Calculate the lapse rate gamma = -(t2-t1)/(z2-z1) c if (gamma .gt. 0.0) then c Assume constant lapse rate atmosphere ps = p1*( (t1/(t1-gamma*z1))**(g/(rd*gamma)) ) ts = t1 + gamma*z1 else c Assume isothermal atmosphere ps = p1*exp(g*z1/(rd*t1)) ts = t1 endif c return end subroutine stndz(p,z,t,theta) c This routine calculates the standard height z (m) from the c pressure p (mb). The temperature t (K) and potential temperature c theta (K) at p are also calculated. C g = 9.80665 r = 287.05 cp = 1004.0 b = 0.0065 p0 = 1013.25 t0 = 288.15 p00 = 1000.0 p1 = 226.32 t1 = 216.65 z1 = 11000.0 cap = r/cp a = r*b/g C z2 = 20000.0 b2 = -0.0010 p2 = 54.75 t2 = t1 a2 = r*b2/g C if (p .ge. p1) then z = (t0/b)*(1.0 - (p/p0)**a) t = t0 - b*z elseif (p .lt. p1 .and. p .ge. p2) then z = z1 + (r*t1/g)*alog(p1/p) t = t1 else z = z2 + (t2/b2)*(1.0 - (p/p2)**a2) t = t2 - b2*(z-z2) endif C theta = t*( (p00/p)**cap ) c return end subroutine thetae(tk,p,rh,pl,tl,w,te) c This routine calculates the equivalent potential temperature c using the formula described in Bolton (1980) MWR c c Input: tk - Temperature (K) c p - Total pressure (hPa) c rh - Relative humidity (%) c c Output: tl - Temperature (K) of the LCL c pl - Pressure (hPa) of the LCL c te - Thetae (K) c w - Mixing ratio (g/Kg) c c Note: If any input values is outside the normal c atmospheric range, all output variables are set to -999. c c Check input if (tk .lt. 100. .or. tk .gt. 350.) go to 900 if (p .le. 0 .or. p .gt. 1200.) go to 900 if (rh .le. 0. .or. rh .gt. 100.) go to 900 c ctk = 273.15 t = tk - ctk c es = 6.112*exp( 17.67*t/(t+243.5) ) ws = 622.0*es/(p-es) w = ws*rh/100.0 c tl = 55.0 + 1.0/( 1.0/(tk-55.0) - alog(rh/100.0)/2840. ) pl = p*(tl/tk)**(1.0/0.2854) c aa = .2854*(1.0 - .00028*w) bb = ( (3.376/tl) - .00254 )*w*(1.+.00081*w) te = (tk*(1000.0/p)**aa)*exp(bb) c return c 900 continue pl = -999. tl = -999. te = -999. c return end subroutine ctorh(x,y,r,theta) c This routine converts from Cartesion coordinates c to radial coordinates, where theta is in c degrees measured clockwise from c the +y-axis (standard meteorological heading). 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 c Convert theta to heading theta = 90.0 - theta if (theta .lt. 0.0) theta = theta + 360.0 c return end subroutine rhtoc(r,thetah,x,y) c This routine converts from radial coordinates c to Cartesian coordinates, where theta is in c degrees measured lockwise from c the +y-axis (standard meteorological heading). c c Convert theta from heading to standard definition theta = 90.0 - thetah if (theta .lt. 0.0) theta = theta + 360.0 c rtd = 57.296 x = r*cos(theta/rtd) y = r*sin(theta/rtd) c return end