program lsdiag c This program performs diagnostic calculations using c large-scale analyses in packed form. c c Last modified April 2005 for 6 hour data increment 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 (1996 version) c call svcal c c Calculate SHIPS variables (1997-2000 version) c call svncal c c Calculate SHIPS variables (2001-2004 version with 5-day input files c call sv5cal c c Calculate SHIPS variables (2005 version with 5-day input files and 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 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,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 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) dimension fxy(ixmx,iymx) c character *4 sname character *130 hline,tline,iline(60) c c Variables for vertical profiles parameter (nvp=9) c dimension pvp(nvp),uvp(nvp),vvp(nvp),tvp(nvp),rvp(nvp) dimension wvp(nvp),wt1(nvp),wt2(nvp),tevp(nvp),tevps(nvp) dimension zvp(nvp) 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 radinn = 200.0 radout = 800.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 momentum fluxes radefc = 600.0 c c Specify minimum points for area average calculations minpts=5 c c Specify vertical pressure levels (mb) for analysis pvp( 1) = 150.0 pvp( 2) = 200.0 pvp( 3) = 250.0 pvp( 4) = 300.0 pvp( 5) = 400.0 pvp( 6) = 500.0 pvp( 7) = 700.0 pvp( 8) = 850.0 pvp( 9) = 1000.0 c c Specify indices of pressures in pvp array i150 = 1 i200 = 2 i250 = 3 i300 = 4 i400 = 5 i500 = 6 i700 = 7 i850 = 8 i000 = 9 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 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 102 format(1x,a4,1x,3(i2),1x,i2) 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 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. '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 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 irhlo(kt) = 9999 irhhi(kt) = 9999 irhmd(kt) = 9999 ipslv(kt) = 9999 20 continue c c Start time loop do 99 kt=0,itmax 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) rawdfn 300 format(' First data file must exist, skip case ',a40) nskip=nskip+1 go to 2001 endif c if (ierr .eq. 1) then c if (ierr .eq. 1 .or. kt .ge. 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 aavg(fxy,radinn,radout,minpts,1,ierr,uvp(n)) c call varget('V',ptem,fxy,1,iexist) call aavg(fxy,radinn,radout,minpts,1,ierr,vvp(n)) c call varget('T',ptem,fxy,1,iexist) call aavg(fxy,radinn,radout,minpts,1,ierr,tvp(n)) c call varget('Z',ptem,fxy,1,iexist) 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 aavg(fxy,radinn,radout,minpts,1,ierr,rvp(n)) else rvp(n) = 999.9 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)) 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 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 write(lulog,750) pvp(ll),uvp(ll),vvp(ll), + wt1(ll),wt2(ll),wt2(ll)/wt1(ll) + ,tvp(ll),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,radefc,refc) irefc(kt) = ifix(refc) 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 Write out synoptic variables for this case write(luout,210) (iu200(k),k=0,itmax) 210 format(10x,21(1x,i4),1x,'U200') 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,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') 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 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') write(luout,235) (irefc(k),k=0,itmax) 235 format(10x,21(1x,i4),1x,'REFC') c write(luout,236) (it000(k),k=0,itmax) 236 format(10x,21(1x,i4),1x,'T000') write(luout,237) (ir000(k),k=0,itmax) 237 format(10x,21(1x,i4),1x,'R000') write(luout,238) (iz000(k),k=0,itmax) 238 format(10x,21(1x,i4),1x,'Z000') 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 sv5cal c This routine calculates the variables needed for c the SHIPS intensity forecast model (2003 version) c for each case in the input file. c c This version uses the input file with 5-day forecasts. parameter (itmax=10) 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,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 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 irhlo(0:itmax),irhhi(0:itmax),irhmd(0:itmax) dimension ipslv(0:itmax) c dimension ie000(0:itmax) dimension iepos(0:itmax),ieneg(0:itmax) dimension iepss(0:itmax),ienss(0:itmax) dimension fxy(ixmx,iymx) c character *4 sname character *80 hline,tline,iline(60) c c Variables for vertical profiles parameter (nvp=9) c dimension pvp(nvp),uvp(nvp),vvp(nvp),tvp(nvp),rvp(nvp) dimension wvp(nvp),wt1(nvp),wt2(nvp),tevp(nvp),tevps(nvp) c c ******************** c c Set iatype=0 to use analyses to calculate predictors, c or iatype=1 to use forecasts for predictors iatype=0 c c If iatype=0, set last year (yyyy) to use re-analysis files iray = 2000 c c Specify radii (km) for vertical profile variables radinn = 200.0 radout = 800.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 momentum fluxes radefc = 600.0 c c Specify minimum points for area average calculations minpts=5 c c Specify vertical pressure levels (mb) for analysis pvp( 1) = 150.0 pvp( 2) = 200.0 pvp( 3) = 250.0 pvp( 4) = 300.0 pvp( 5) = 400.0 pvp( 6) = 500.0 pvp( 7) = 700.0 pvp( 8) = 850.0 pvp( 9) = 1000.0 c c Specify indices of pressures in pvp array i150 = 1 i200 = 2 i250 = 3 i300 = 4 i400 = 5 i500 = 6 i700 = 7 i850 = 8 i000 = 9 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 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(a80) c c Extract storm name and date information from header line read(hline,102) sname,imyear,immon,imday,imtime 102 format(1x,a4,1x,3(i2),1x,i2) 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 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(67:70) .eq. 'LAT ') then read(tline,105) (islat(k),k=-2,itmax) 105 format(13(1x,i4)) endif C if (tline(67:70) .eq. 'LON ') then read(tline,105) (islon(k),k=-2,itmax) endif c if (tline(67:70) .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 15 continue c c Initialize synoptic variables to missing do 20 kt=0,itmax iu200(kt) = 9999 it200(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 irhlo(kt) = 9999 irhhi(kt) = 9999 irhmd(kt) = 9999 ipslv(kt) = 9999 20 continue c c Start time loop do 99 kt=0,itmax 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 if (imyear4 .le. iray) then rawdfn = 'R00055_X0000_PACK.DAT' else rawdfn = 'A00055_X0000_PACK.DAT' endif c if (imtime .eq. 6 .or. imtime .eq. 18) then rawdfn(8:8) = 'Y' endif c ktime = kt*12 c c Calculate date information for next analysis file call tadd(imyear,immon,imday,imtime,ktime, + iayear,iamon,iaday,iatime) 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*12 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-1) .lt. 900.0 .and. + slat(kt ) .lt. 900.0) then t2lat = slat(kt) t2lon = slon(kt) t1lat = slat(kt-1) t1lon = slon(kt-1) icok = 1 else if (slat(kt+1) .lt. 900.0 .and. + slat(kt ) .lt. 900.0) then t2lat = slat(kt+1) t2lon = slon(kt+1) t1lat = slat(kt) t1lon = slon(kt) icok = 1 else t2lat = slat(kt) t2lon = slon(kt) t1lat = slat(kt) t1lon = slon(kt) 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)/(12.0*3600.0) cy = 1000.0*(temy2-temy1)/(12.0*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 aavg(fxy,radinn,radout,minpts,1,ierr,uvp(n)) c call varget('V',ptem,fxy,1,iexist) call aavg(fxy,radinn,radout,minpts,1,ierr,vvp(n)) c call varget('T',ptem,fxy,1,iexist) call aavg(fxy,radinn,radout,minpts,1,ierr,tvp(n)) c if (ptem .ge. 300.0) then call varget('R',ptem,fxy,1,iexist) call aavg(fxy,radinn,radout,minpts,1,ierr,rvp(n)) else rvp(n) = 999.9 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)) 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 200 mb (in deg C) it200(kt) = ifix( 10.0*(tvp(i200)-273.15) ) 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 write(lulog,750) pvp(ll),uvp(ll),vvp(ll), + wt1(ll),wt2(ll),wt2(ll)/wt1(ll) + ,tvp(ll),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,radefc,refc) irefc(kt) = ifix(refc) 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 Write out synoptic variables for this case write(luout,210) (iu200(k),k=0,itmax) 210 format(10x,11(1x,i4),1x,'U200') write(luout,211) (it200(k),k=0,itmax) 211 format(10x,11(1x,i4),1x,'T200') write(luout,212) (ie000(k),k=0,itmax) 212 format(10x,11(1x,i4),1x,'E000') write(luout,213) (iepos(k),k=0,itmax) 213 format(10x,11(1x,i4),1x,'EPOS') write(luout,214) (ieneg(k),k=0,itmax) 214 format(10x,11(1x,i4),1x,'ENEG') write(luout,226) (iepss(k),k=0,itmax) 226 format(10x,11(1x,i4),1x,'EPSS') write(luout,238) (ienss(k),k=0,itmax) 238 format(10x,11(1x,i4),1x,'ENSS') write(luout,216) (irhlo(k),k=0,itmax) 216 format(10x,11(1x,i4),1x,'RHLO') write(luout,217) (irhmd(k),k=0,itmax) 217 format(10x,11(1x,i4),1x,'RHMD') write(luout,218) (irhhi(k),k=0,itmax) 218 format(10x,11(1x,i4),1x,'RHHI') write(luout,220) (ishrd(k),k=0,itmax) 220 format(10x,11(1x,i4),1x,'SHRD') write(luout,222) (ishtd(k),k=0,itmax) 222 format(10x,11(1x,i4),1x,'SHTD') write(luout,225) (ishrs(k),k=0,itmax) 225 format(10x,11(1x,i4),1x,'SHRS') write(luout,227) (ishts(k),k=0,itmax) 227 format(10x,11(1x,i4),1x,'SHTS') write(luout,228) (ishrg(k),k=0,itmax) 228 format(10x,11(1x,i4),1x,'SHRG') c if (imxshr .eq. 1) then write(luout,229) (ishxu(k),k=0,itmax) 229 format(10x,11(1x,i4),1x,'SHXU') write(luout,230) (ishxl(k),k=0,itmax) 230 format(10x,11(1x,i4),1x,'SHXL') endif c write(luout,231) (ipslv(k),k=0,itmax) 231 format(10x,11(1x,i4),1x,'PSLV') write(luout,232) (ivorb(k),k=0,itmax) 232 format(10x,11(1x,i4),1x,'Z850') write(luout,233) (idivb(k),k=0,itmax) 233 format(10x,11(1x,i4),1x,'D200') write(luout,235) (irefc(k),k=0,itmax) 235 format(10x,11(1x,i4),1x,'REFC') 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 svncal c This routine calculates the variables needed for c the SHIPS intensity forecast model (1998 version) c 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,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,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:7),islon(-2:7) dimension slat(-2:7), slon(-2:7) c dimension iu200(0:6),it200(0:6) dimension ishrd(0:6),ishrs(0:6) dimension ishtd(0:6),ishts(0:6) dimension ivorb(0:6),idivb(0:6),irefc(0:6) dimension if850(0:6),if500(0:6) c dimension fxy(ixmx,iymx) c character *4 sname character *72 hline,tline,iline(50) c c ******************** c c Specify radius (km) for shear calculations radksh=600.0 c c Specify radius (km) for area average of u200,t200 radkut=1000.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 momentum fluxes radefc = 600.0 c c Specify radius (km) for horizontal deformation calculation raddef = 600.0 c c Specify minimum points for area average calculations minpts=5 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(a72) c c Extract storm name and date information from header line read(hline,102) sname,imyear,immon,imday,imtime 102 format(1x,a4,1x,3(i2),1x,i2) c go to 5100 751 write(lulog,*) ' Error reading input file (luinp)' stop 5100 continue C iline(1) = hline 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(52:54) .eq. 'LAT') then read(tline,105) (islat(k),k=-2,6) 105 format(9(1x,i4)) endif C if (tline(52:54) .eq. 'LON') then read(tline,105) (islon(k),k=-2,6) endif c if (tline(52:55) .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:',/,A72) C do 15 k=-2,6 slat(k) = float(islat(k))/10.0 slon(k) = -1.0*float(islon(k))/10.0 15 continue c c Initialize synoptic variables to missing do 20 kt=0,6 iu200(kt) = 9999 it200(kt) = 9999 ishrd(kt) = 9999 ishtd(kt) = 9999 ishrs(kt) = 9999 ishts(kt) = 9999 ivorb(kt) = 9999 idivb(kt) = 9999 irefc(kt) = 9999 if850(kt) = 9999 if500(kt) = 9999 20 continue c c Start time loop do 99 kt=0,4 c c Check for missing track positions if (slat(kt) .gt. 900.0) go to 99 c c Construct rawd file rawdfn = 'S00055_X0000_PACK.DAT' c if (imtime .eq. 6 .or. imtime .eq. 18) then rawdfn(1:1) = 'T' endif c if (sname(1:1) .ne. ' ') then rawdfn(8:8) = sname(1:1) endif c ktime = kt*12 write(rawdfn(2:4),110) ktime 110 format(i3.3) c write(rawdfn(5:6),115) imyear 115 format(i2.2) c write(rawdfn(9:10),115) immon c if (imtime .ge. 12) then imdayt = imday + 50 else imdayt = imday endif write(rawdfn(11:12),115) imdayt 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 area averaged u at 200 mb ptem=200.0 call varget('U',ptem,fxy,0,iexist) if (iexist .eq. 1) then call havg(fxy,radkut,minpts,1,ierr,ubar) ubar = 1.944*ubar if (ierr .eq. 0) then iu200(kt) = ifix(10.0*ubar) else iu200(kt) = 9999 endif else iu200(kt) = 9999 endif c c Calculate area averaged T at 200 mb ptem=200.0 call varget('T',ptem,fxy,0,iexist) if (iexist .eq. 1) then call havg(fxy,radkut,minpts,1,ierr,tbar) tbar = tbar - 273.15 if (ierr .eq. 0) then it200(kt) = ifix(10.0*tbar) else it200(kt) = 9999 endif else it200(kt) = 9999 endif c c Calculate 850 mb horizontal deformation ptem = 850.0 call defcal(ptem,fxy,ierr) if (ierr .eq. 0) then call havg(fxy,raddef,minpts,1,ierr,defbar) if (ierr .eq. 0) then if850(kt) = ifix((1.0e+7)*defbar) else if850(kt) = 9999 endif else if850(kt) = 9999 endif c c Calculate 500 mb horizontal deformation ptem = 500.0 call defcal(ptem,fxy,ierr) if (ierr .eq. 0) then call havg(fxy,raddef,minpts,1,ierr,defbar) if (ierr .eq. 0) then if500(kt) = ifix((1.0e+7)*defbar) else if500(kt) = 9999 endif else if500(kt) = 9999 endif c c Calculate the 850-200 mb vertical shear ptem1 = 850.0 ptem2 = 200.0 call vshear(ptem1,ptem2,slon(kt),slat(kt),radksh, + minpts,shrd,shtd) ishrd(kt) = ifix(10.0*shrd) ishtd(kt) = ifix(shtd) c c Calculate the 850-500 mb vertical shear ptem1 = 850.0 ptem2 = 500.0 call vshear(ptem1,ptem2,slon(kt),slat(kt),radksh, + minpts,shrs,shts) ishrs(kt) = ifix(10.0*shrs) ishts(kt) = ifix(shts) 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 x,y components of storm motion if (slat(KT-1) .lt. 900.0 .and. + slat(KT ) .lt. 900.0) then t2lat = slat(kt) t2lon = slon(kt) t1lat = slat(kt-1) t1lon = slon(kt-1) else t2lat = slat(kt+1) t2lon = slon(kt+1) t1lat = slat(kt) t1lon = slon(kt) 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)/(12.0*3600.0) cy = 1000.0*(temy2-temy1)/(12.0*3600.0) c c Calculate relative eddy flux convergence of angular c momentum ptem = 200.0 call refcal(ptem,slon(kt),slat(kt),cx,cy,radefc,refc) irefc(kt) = ifix(refc) 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 Write out synoptic variables for this case write(luout,210) (iu200(k),k=0,6) 210 format(10x,7(1x,i4),6x,'U200') write(luout,215) (it200(k),k=0,6) 215 format(10x,7(1x,i4),6x,'T200') write(luout,220) (ishrd(k),k=0,6) 220 format(10x,7(1x,i4),6x,'SHRD') write(luout,222) (ishtd(k),k=0,6) 222 format(10x,7(1x,i4),6x,'SHTD') write(luout,225) (ishrs(k),k=0,6) 225 format(10x,7(1x,i4),6x,'SHRS') write(luout,227) (ishts(k),k=0,6) 227 format(10x,7(1x,i4),6x,'SHTS') write(luout,230) (ivorb(k),k=0,6) 230 format(10x,7(1x,i4),6x,'Z850') write(luout,232) (idivb(k),k=0,6) 232 format(10x,7(1x,i4),6x,'D200') write(luout,235) (irefc(k),k=0,6) 235 format(10x,7(1x,i4),6x,'REFC') write(luout,236) (if850(k),k=0,6) 236 format(10x,7(1x,i4),6x,'F850') write(luout,237) (if500(k),k=0,6) 237 format(10x,7(1x,i4),6x,'F500') 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 svcal c This routine calculates the variables needed for c the SHIPS intensity forecast model (1996 version) c 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,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 character *72 hline,tline,iline(15) c dimension islat(-2:7),islon(-2:7) dimension slat(-2:7), slon(-2:7) dimension shear(0:6) dimension fb(irmx),pefc(irmx),refc(irmx) c dimension itemp(nplevh) c parameter (nflux=5) dimension ipefc(nflux),irefc(nflux) c c ******************** c c Generic rawd file name for the input rawd files rawdfn = 'A00000_X0000_PACK.DAT' c c Start the case by case loop nrawd=0 nskip=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,10 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 do 15 k=-2,6 slat(k) = float(islat(k))/10.0 slon(k) = -1.0*float(islon(k))/10.0 15 continue c c Create rawd file name rawdfn(5:6) = HLINE(7:8) c rawdfn(8:8) = HLINE(2:2) rawdfn(9:10) = HLINE(9:10) c read(hline(11:12),400) iday read(hline(14:15),400) itime if (itime .ge. 12) iday = iday+50 write(rawdfn(11:12),405) iday 400 format(i2) 405 format(i2.2) c c Unpack the rawd file call aunpack(ierr) if (ierr .ne. 0) then write(6,205) write(lulog,205) 205 format(1x,' Data file error, this case skipped') nskip = nskip + 1 go to 2001 endif c c Calculate storm speed, storm size, shear, and eddy fluxes call ssfcal(slon,slat,ssx,ssy,ssize,shear,refc,pefc) c c Calculate vertical soundings sradk=1000.0 call vsound(slon(0),slat(0),sradk) c c Modify input lines and write to output file c c ** Write storm size on line 1** isize = ifix(sign(0.5,ssize)+10.0*ssize) tline = iline(1) write(tline(27:30),410) isize 410 format(i4) c c ** Write momentum fluxes on line 1** do 20 i=1,nflux ravg=0.0 pavg=0.0 ks = 3*(i-1) + 2 ke = ks + 2 rnk = float(1+ke-ks) do 25 k=ks,ke ravg = ravg + refc(k) pavg = pavg + pefc(k) 25 continue ravg = ravg/rnk pavg = pavg/rnk irefc(i) = ifix(sign(0.50,ravg) + ravg) ipefc(i) = ifix(sign(0.50,pavg) + pavg) 20 continue c do 30 k=1,nflux ks = 32 + 4*(k-1) ke = ks + 3 ls = 32 + 4*nflux + 4*(k-1) le = ls + 3 write(tline(ks:ke),410) irefc(k) write(tline(ls:le),410) ipefc(k) 30 continue iline(1) = tline c c ** Write vertical shear on line 5 ** tline = iline(5) do 35 k=0,6 ishear = ifix(0.5+10.0*shear(k)) ks = 12 + k*5 ke = ks + 3 write(tline(ks:ke),410) ishear 35 continue iline(5) = tline c c ** Write first 10 lines of output file ** do 40 i=1,10 write(luout,100) iline(i) 40 continue c c ** Write vertical sounding variables on output file ** do 45 k=1,nplevh itemp(k) = ifix(0.5+10.0*up(k)) 45 continue write(luout,420) (itemp(k),k=1,nplevh) 420 format(16(1x,i4)) c do 50 k=1,nplevh itemp(k) = ifix(0.5+10.0*vp(k)) 50 continue write(luout,420) (itemp(k),k=1,nplevh) c do 55 k=1,nplevh itemp(k) = ifix(0.5+zp(k)) if (itemp(k) .gt. 9990 .or. itemp(k) .lt. -999) + itemp(k) = 9999 55 continue write(luout,420) (itemp(k),k=1,nplevh) c do 60 k=1,nplevh itemp(k) = ifix(0.5+10.0*tp(k)) 60 continue write(luout,420) (itemp(k),k=1,nplevh) c do 65 k=1,nplevh itemp(k) = ifix(0.5+10.0*rp(k)) 65 continue write(luout,420) (itemp(k),k=1,nplevh) c do 70 k=1,nplevh itemp(k) = ifix(plevh(k)) 70 continue write(luout,420) (itemp(k),k=1,nplevh) 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 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 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 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 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 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 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