program calccai_nam_10day implicit none integer m,nc,nr,i,j,nfiles,nrecs,np,k,cnc,cnr,kk,ntime,nt,ntleap parameter(nfiles=57,nc=144,nr=37,nrecs=80000,np=20,nt=365) parameter(cnc=144,cnr=73,ntime=672,ntleap=366) integer x(nrecs),y(nrecs),krec,startx,starty,endx,endy integer mm,nn,pp,ii,storm(nrecs),bb,cc,p,jj real dx,dy,slat,slon real data(20),lat(nrecs),lon(nrecs),day(nrecs) integer month(nrecs),jday(nrecs) real slp(nrecs), dslp(nrecs),year(nrecs) real cnt(nc,nr,12),sum(nc,nr,12) real cslp(cnc,cnr,nt),fcslp(nc,nr,nt) real cslpleap(cnc,cnr,ntleap),fcslpleap(nc,nr,ntleap) real tempslp(cnc,cnr),cai(nc,nr),regsum1(ntime,8) real regsum2(ntime,8),regsum3(ntime,8) real msum(nc,nr,12),mcnt(nc,nr,12),mavg(nc,nr,12) real tempmavg(nc,nr),gsum,gcnt,gavg integer mstmp,ystmp,flag,r,s,ten1 character*55 filename character*4 outfile ******************************************************************** * Set starting points of output grid for binning of lat/lon from ASCII * storm track files slat = 0.0 slon = 0.0 dx = 2.5 dy = 2.5 ******************************************************************** ******************************************************************** * Read in climatological daily SLP (non-leap year first, then leap year) * Place in different arrays. Easiest to have two separate arrays. open(30,file="mslp_daily_clim_7995.bin", & status='old',access='direct', & recl=cnc*cnr*4) do m=1,nt read(30,rec=m) tempslp do i=1,cnr do j=1,cnc cslp(j,i,m) = tempslp(j,i) c print*, i,j,m,cslp(j,i,m) enddo enddo enddo close(30) open(31,file="mslp_daily_clim_7995_leap.bin", & status='old',access='direct', & recl=cnc*cnr*4) do m=1,ntleap read(31,rec=m) tempslp do i=1,cnr do j=1,cnc cslpleap(j,i,m) = tempslp(j,i) c print*, i,j,m,cslp(j,i,m) enddo enddo enddo close(31) ********************************************************************** ********************************************************************** * Subset grid for the Northern Hemisphere centered on 0.0 (row 37 of * original grid). Two separate arrays for leap and non-leap years. do m=1,nt do j=1,cnc ii = 0 do i=37,cnr c print*, i,j,cslp(j,i,m)/100.0 ii = ii + 1 fcslp(j,ii,m) = cslp(j,i,m)/100.0 c print*, m,i,j,ii,cslp(j,i,m),fcslp(j,ii,m) enddo enddo enddo do m=1,ntleap do j=1,cnc ii = 0 do i=37,cnr ii = ii + 1 fcslpleap(j,ii,m) = cslpleap(j,i,m)/100.0 c print*, m,i,j,ii,cslp(j,i,m),fcslp(j,ii,m) enddo enddo enddo *********************************************************************** *********************************************************************** * Calculation of CAI * Read input filename list. The files in this list are looped through. open(10,file="./input_jan2007_cai_list",status='old') * Open output file for calculated CAI as a function of grid point open(20,file="./CAI_10daysum_jan_2007_NA",status='unknown', & access='direct',recl=100*21*4) * Now loop through storm track ASCII files for each year do m=1,nfiles * Initialize counter and sums do pp=1,3 ! 3 10 day period of January do mm=1,nc do nn=1,nr cnt(mm,nn,pp) = 0.0 sum(mm,nn,pp) = 0.0 enddo enddo enddo * Read filename from input list file and then yearly storm track file read(10,'(A55)') filename outfile = filename(32:35) if(outfile.eq."1952".or.outfile.eq."1956".or.outfile.eq."1960".or. & outfile.eq."1964".or.outfile.eq."1968".or.outfile.eq."1972".or. & outfile.eq."1976".or.outfile.eq."1980".or.outfile.eq."1984".or. & outfile.eq."1988".or.outfile.eq."1992".or.outfile.eq."1996".or. & outfile.eq."2000".or.outfile.eq."2004".or.outfile.eq."2008")then flag = 1 else flag = 0 endif print*, filename print*, outfile print*, flag open(11,file=filename,status='old') do i=1,nrecs read(11,*,end=999) (data(k),k=1,np) * Storing needed information for CAI calculation year(i)=data(3) month(i)=int(data(4)) day(i)=data(5) jday(i) = int(data(1)) slp(i)=data(10) lat(i) = data(14) lon(i) = data(15) storm(i) = data(20) if (jday(i) .lt. 31.0) then c print*, year(i),month(i),day(i),slp(i),lat(i),lon(i),storm(i) if (jday(i).lt.11.0) then ten1 = 1 elseif (jday(i).ge.11.0.and.jday(i).lt.21.0) then ten1 = 2 else ten1 = 3 endif print*, jday(i),ten1 * Bin individual storm centers into the GR-2. Determine grid point. y(i) = int( (lat(i) - slat) / dy) + 1.0 x(i) = int( (lon(i) - slon) / dx) + 1.0 * Calculate difference between storm SLP and daily climatological SLP if (flag .eq. 0) then dslp(i) = slp(i) - fcslp(x(i),y(i),jday(i)) else dslp(i) = slp(i) - fcslpleap(x(i),y(i),jday(i)) endif * Sum over time period. Here it is over months. cnt(x(i),y(i),ten1) = cnt(x(i),y(i),ten1) + 1.0 sum(x(i),y(i),ten1) = sum(x(i),y(i),ten1) + dslp(i) c print*, "I= ",i,month(i),lat(i),lon(i),y(i),x(i) c print*, "I= ",i,jday(i),fcslp(x(i),y(i),jday(i)),slp(i),dslp(i) c write(*,66) i,y(i),x(i),month(i),slp(i),fcslp(x(i),y(i), c & month(i)), c & dslp(i), c & cnt(x(i),y(i),month(i)),sum(x(i),y(i),month(i)) c 66 format(i6,1x,i3,1x,i3,1x,i2,1x,f8.2,1x,f8.2,1x,f8.2,1x, c & f5.1,1x,f10.2) endif enddo ****************************************************************************** 999 continue ****************************************************************************** * Write grid of monthly values for year ccc krec = ((m*12) - 12) + 1 krec = ((m*3) - 3) + 1 ccc do k=1,12 do k=1,3 c print*, k, krec do j=1,nc do i=1,nr cai(j,i) = -1.0 * sum(j,i,k) enddo enddo **** 45-144 (110 - 357.5 E), 9-29 (20-70 N) write(20,rec=krec) ((cai(r,s),r=45,144),s=9,29) krec = krec + 1 77 format(a4,1x,i3,1x,f7.2) enddo close(11) enddo close(10) close(20) end