date: Mon Mar 4 13:58:40 2002 from: Ian Harris subject: getgales.for to: r.teck@uea.ac.uk program getgales c c Temporarily-non-interactive version. Hardcoded in the data definitions c (lats, lons, ilats). Lats must be multiples of 2.5; lons are in between c gridcells so are multiples of 2.5 + 1.25. ilats are the row numbers for c the coefficients arrays (ax,bx,cx,dx). An ilat of 1 would give the c coefficients for 75N, wheras a value of 5 corresponds to 65N. c c Interactive program to produce gale indices for a given c location, using NCEP Reanalysis Data and routines from c Phil Jones. Preset for the five Holsmeer areas: c iceland 23W 65N (nearest = -23.75, 65) c gbight 9.5E 55N (nearest = 8.75, 55) c ireland 8.4W 54.75N (nearest = -8.75, 55) c portugal 9W 38.5N (nearest = -8.75, 37.5) c w england 5W 55N (nearest = -3.75, 55) - for Phil J integer i,j,k,m,n,sy,ey,nd,nl,dimo,err,cyr,cmo,cdy,nyr,nmo,ndy, * ilats(5) real lats(5),lons(5),grid(21,53),w,s,zw,zs,f,z,g,p,wd, * ax(17),bx(17),cx(17),dx(17) character filroot*47,fil*26,outf*20 logical dbg data dbg/.true./ data filroot/'/cru/ncep/reanalysis/window1/daily/surface/slp/'/ data cyr/0/,cmo/0/,cdy/0/,nyr/0/,nmo/0/,ndy/0/ c *** change the next two lines to alter the locations *** data lats/65,55,55,37.5,55/,lons/-23.75,8.75,-8.75,-8.75,-3.75/ data ilats/5,9,9,16,9/ c *** change the next line to limit the temporal extent (daily files) *** data sy/1958/,ey/2001/ call getconstants(ax,bx,cx,dx) c open(10,status='unknown',file='constants.dat') c do k = 1,17 c write(10,'(f4.1,4f8.3)')(75-(k-1)*2.5),ax(k),bx(k),cx(k),dx(k) c enddo c close(10) c stop do k = 1,5 if (lons(k).lt.0) then write(outf,'("gales.",i4.4,"N",i4.4,"W.dat")') * int(lats(k)*100),-1*int(lons(k)*100) else write(outf,'("gales.",i4.4,"N",i4.4,"E.dat")') * int(lats(k)*100),int(lons(k)*100) endif open(10+k,status='unknown',file=outf) enddo write(*,'(/"GETGALES - Extract Gale Indices from NCEP Data")') c write(*,'(/"Please enter a latitude between c * 33N and 77N (33 and 77): ",$)') c read(*,*)inlat c write(*,'(/"Please enter a longitude between c * 52W and 62E (-52 and 62): ",$)') c read(*,*)inlon c write(*,'(/"Please enter a start year between c * 1958 and 1998: ",$)') c read(*,'(i4)')sy c write(*,'(/"Please enter and end year between ",i4, c *" and 1998: ",$)')sy c read(*,'(i4)')ey c write(*,'(/"Please enter an output filename: ",$)') c read(*,'(a)')outf c open(11,status='unknown',file=outf) cyr = sy cmo = 1 cdy = 0 nyr = 0 nmo = 0 ndy = 0 do i=sy,ey write(fil,'("slp.",i4,".d.w1.53x21.dat.gz")')i write(*,'(/"Processing ",a)')fil call system('gunzip -c '//filroot//fil//' > gg.tmp') call system('wc -l gg.tmp > gg.wcl') open(10,status='old',file='gg.wcl') read(10,*)nl close(10) if (mod(nl,21).ne.0) then write(*,'("ERROR: Unexpected line count (",i10.10, * ") in ",a)')nl,filroot//fil stop endif nd = nl/21 if (nd.ne.(337+dimo(i,2))) then write(*,'("ERROR: Unexpected day count (",i3.3, * ") in ",a)')nd,filroot//fil stop endif open(10,status='old',file='gg.tmp') do j = 1,nd call tnext(cyr,cmo,cdy,-1,-1,-1,nyr,nmo,ndy,k,m,n) if (nyr.ne.i) then write(*,'("ERROR: year mismatch from tnext: ",3i6)')nyr,i stop endif cyr = nyr cmo = nmo cdy = ndy do k = 1,21 read(10,'(53f12.0)')(grid(k,m),m=1,53) enddo do k = 1,5 c call gale2(grid,lats(k),lons(k),w,s,zw,zs,f,z,g,p,wd,err) m = ilats(k) ! get constants for this latitude if (dbg) write(*,'("A1: latlon: ",2f8.2)')lats(k),lons(k) if (dbg) write(*,'("A2: abcd: ",4f8.2)')ax(m),bx(m),cx(m),dx(m) call gale2(grid,lats(k),lons(k),ax(m),bx(m),cx(m),dx(m), * g,dbg,err) if (err.eq.0) then write(10+k,'(3i4,f10.4)')i,cmo,cdy,g else write(10+k,'(3i4,f10.4)')i,cmo,cdy,mv endif enddo enddo close(10) enddo do k = 1,5 close(10+k) enddo call system('rm gg.tmp gg.wcl') stop end subroutine gale2(grid,rlat,rlon,a,b,c,d,g,dbg,err) c *** Passed a standard grid (ie from NCEP Reanalysis Data), and a c *** lat/lon pair, this routine returns the gale index. c *** Grid must be 21r x 53c, ie 'window 1' from the CRU NCEP Data c *** pages. First row is 80N. Lats and lons must give c *** [lat = multiple of 2.5] and [lon = multiple of 2.5 + 1.25]. c *** Lats must run from -60 (west) to +70 (east). integer i,j,k,ilat,ilon,er real rlat,rlon,lat,lon,grid(21,53),sgrid(5,8),w,s,zw,zs,f, * z,g,ml,mv logical dbg data mv/-9999.00/ lat = rlat ! preserve anonymity of 'global' lon = rlon ! lat and lon fields if (dbg) write(*,'(" Fa")') err = 0 ilat = nint(33-(lat/2.5)) ! convert lat to array row if (dbg) write(*,'(" Fb: ilat: ",i3)')ilat lon = lon + 60 ! false lon to start from 0 not -60 ilon = nint((lon-1.25)/2.5 - 2) ! convert lon to westernmost array column if (dbg) write(*,'(" Fc: ilon: ",i3)')ilon do i = ilat-2,ilat+2 ! extract & scale appropriate sgrid do j = ilon,ilon+7 if (grid(i,j).eq.mv) then err = 1 return endif sgrid(i-ilat+3,j-ilon+1) = grid(i,j)/100 enddo enddo if (dbg) then write(*,'(" Fd: sgrid:")') do i = 1,5 write(*,'(" ",8i8)')(nint(sgrid(i,j)),j=1,8) enddo endif c next bit lifted from clarejcindices.for (Clare Goodess) w=0.25*(sgrid(4,3)+sgrid(4,4)+sgrid(4,5)+sgrid(4,6)) - * 0.25*(sgrid(2,3)+sgrid(2,4)+sgrid(2,5)+sgrid(2,6)) if (dbg) write(*,'(" Fe: w: ",f8.2)')w s=a*(0.125*(sgrid(2,6)+2*sgrid(3,6)+sgrid(4,6)+sgrid(2,5)+ * 2*sgrid(3,5)+sgrid(4,5))-0.125*(sgrid(2,3)+2*sgrid(3,3)+ * sgrid(4,3)+sgrid(2,4)+2*sgrid(3,4)+sgrid(4,4))) f = (s**2 + w**2)**0.5 if (dbg) write(*,'(" Ff: s,f: "2f8.2)')s,f zw=b*(0.25*(sgrid(5,3)+sgrid(5,4)+sgrid(5,5)+sgrid(5,6)) - * 0.25*(sgrid(3,3)+sgrid(3,4)+sgrid(3,5)+sgrid(3,6))) - * c*(0.25*(sgrid(3,3)+sgrid(3,4)+sgrid(3,5)+sgrid(3,6)) - * 0.25*(sgrid(1,3)+sgrid(1,4)+sgrid(1,5)+sgrid(1,6))) if (dbg) write(*,'(" Fg: zw: "f8.2)')zw zs=d*( 0.125*(sgrid(2,8)+2*sgrid(3,8)+sgrid(4,8)+ * sgrid(2,7)+2*sgrid(3,7)+sgrid(4,7)) - * 0.125*(sgrid(2,6)+2*sgrid(3,6)+sgrid(4,6)+ * sgrid(2,5)+2*sgrid(3,5)+sgrid(4,5)) ) - * d*( 0.125*(sgrid(2,4)+2*sgrid(3,4)+sgrid(4,4)+ * sgrid(2,3)+2*sgrid(3,3)+sgrid(4,3)) - * 0.125*(sgrid(2,2)+2*sgrid(3,2)+sgrid(4,2)+ * sgrid(2,1)+2*sgrid(3,1)+sgrid(4,1)) ) z = zw+zs g = (f**2+0.25*z**2)**0.5 ! Gale Index - the bit we actually *want* if (dbg) write(*,'(" Fh: zs,z,g: ",3f8.2)')zs,z,g return end subroutine getconstants(ax,bx,cx,dx) c Clare G's subroutine to calculate the scaling factors for different c latitudes. Could probably be improved. real centrey,ax(17),bx(17),cx(17),dx(17),xlat,xrat centrey = 75 xlat = 2.5 xrat = 1.75 do i = 1,17 ax(i) = 1/cosd(centrey) bx(i) = sind(centrey)/sind(centrey-xlat) cx(i) = sind(centrey)/sind(centrey+xlat) dx(i) = xrat/(4*cosd(centrey)**2) centrey = centrey - xlat enddo return end subroutine tnext(cyr,cmo,cdy,chr,cmi,cse,nyr,nmo,ndy,nhr,nmi,nse) c Annoying subroutine that, passed a date and time, returns the next c date and time in the sequence. Now you'd think that would be a bit c of a lifesaver, rather than 'annoying' - especially as you can give c it any temporal resolution by filling in 0..n positions at start c and/or end with '-1' - but instead it's annoying simply because I c already wrote it once, last year, and now I can't find the damn thing. c c Example usage: c call tnext(1998,11,30,14,59,-1,a,b,c,d,e,f) c Returns: c a 1998 c b 11 c c 30 c d 15 c e 0 c f -1 integer cyr,cmo,cdy,chr,cmi,cse,nyr,nmo,ndy,nhr,nmi,nse, * roll,dimo roll = 1 ! flag to indicate roll-up still required if (cse.eq.-1) then nse = -1 elseif (cse.eq.59) then nse = 0 else nse = cse + 1 roll = 0 endif if (cmi.eq.-1) then nmi = -1 elseif (roll.eq.0) then nmi = cmi else if (cmi.eq.59) then nmi = 0 else nmi = cmi + 1 roll = 0 endif endif if (chr.eq.-1) then nhr = -1 elseif (roll.eq.0) then nhr = chr else if (chr.eq.23) then nhr = 0 else nhr = chr + 1 roll = 0 endif endif if (cdy.eq.-1) then ndy = -1 elseif (roll.eq.0) then ndy = cdy else if (cdy.eq.dimo(cyr,cmo)) then ndy = 1 else ndy = cdy + 1 roll = 0 endif endif if (cmo.eq.-1) then nmo = -1 elseif (roll.eq.0) then nmo = cmo else if (cmo.eq.12) then nmo = 1 else nmo = cmo + 1 roll = 0 endif endif if (cyr.eq.-1) then nyr = -1 elseif (roll.eq.0) then nyr = cyr else nyr = cyr + 1 endif return end integer function dimo(yr,mo) integer yr,mo if (mo.eq.4.or.mo.eq.6.or.mo.eq.9.or.mo.eq.11) then dimo = 30 elseif (mo.ne.2) then dimo = 31 else dimo = 28 if ((mod(yr,4).eq.0.and.mod(yr,400).ne.0).or.(yr.eq.2000)) * dimo = 29 endif return end