c weighted average gridding c created on 13.6.97 c modified for year and box average on 4-8-97 c modified for f90 on 8-10-97 c modified for w/o yr, mn, dt on 28-10-97 c modified for idys problem on 12-10-99 c modified for flexible julian days and for all years of 10 day avrg on 9-9-1 c modified for dimension check on 24-6-2010 c for all years of 10-day (eq) avrg, like for rrr, is rqrd remove cmnt at %%% character *30 if,of,gm*1,mnthly*1,ymd*1,all*1,std*1 character *1 tmpcnv c parameter (m1=362,m2=242,m3=2) parameter (m1=5,m2=4,m3=14) dimension sumw(m1,m2,m3),parm(0:25),idys(0:12),xiy(m3),xim(m3) dimension iyrs(0:64),idysl(0:12),hrzntl(m2),xid(m3) dimension gy(m2),gx(m1) c m1: lat; m2: lon; m3: prd no integer *4 em,ed,cn,xmnlmt,ymnlmt,xmxlmt,ymxlmt,idffy real *8 sumz(m1,m2,m3),stdz(m1,m2,m3) data idys/0,31,59,90,120,151,181,212,243,273,304,334,365/ data idysl/0,31,60,91,121,152,182,213,244,274,305,335,366/ c for 1950 number of days to be added is 0 c for 1951 number of days to be added are 365 (days in 1950) c for 1952 number of days to be added are 365 (days in 1950+51) c for 1953 number of days to be added are 365 (days in 1950+51+52) c so the leap year factor comes for 53 (not for 52) ntd=-365 std='n' !for standard dev if (std.eq.'y') then print * print *,' std calculations!! ' print * end if do i=1950,2014 ichk=mod(i,4) if (ichk.eq.1) then ndd=366 else ndd=365 end if ntd=ntd+ndd iyrs(i-1950)=ntd end do nf=1 ! no of separate files to be read. Rqrd to read two big files separetely all='y' print '(a,$)',' Box/Weighted average (b/w): ' read '(a)',gm if (gm.eq.'w'.and.std.eq.'y') then print * print *,' Standard Deviation not possible with 1 weighted averaging' stop end if mxd=366 mid=366 rmnv=99999 rmxv=-99999 nff=0 ! for checking the file nos ncd=0 nn=0 14 print '(/,a,$)',' input: ' read '(a)',if open (unit=1,file=if,status='old') if (nff.ge.1) go to 12 print '(/,a,$)',' output: ' read '(a)',of open (unit=12,file=of,status='unknown') print '(/,a,$)',' min and max lat: ' read *,ymn,ymx print '(/,a,$)',' min and max lon: ' read *,xmn,xmx print '(/,a,$)',' lat/lon grid sizes: ' read *,dy,dx print '(/,a,$)',' tot no of parm to be read: ' read *,np print '(/,a,$)',' column nos of parms to be gridded: ' read *,cn print '(/,a,$)',' yr/mon/date column nos: ' read *,iy,mon,idt print '(/,a,$)',' value not to be considered: ' read *,thv print '(/,a,$)',' value to be written when not available: ' read *,vna print '(/,a,$)',' lat/lon column nos: ' read *,lac,loc ymd='n' if (iy.ne.0.and.mon.ne.0.and.idt.ne.0) then ymd='y' print '(/,a,$)',' strtng year month and date: ' read *,aisy,aism,aisd isy=aisy ism=aism isd=aisd if (isy.ge.0.and.isy.lt.50) then isy=isy+2000 else if (isy.ge.50.and.isy.le.99) then isy=isy+1900 end if c if (isy.gt.1900) isy=isy-1900 print '(/,a,$)',' endng year month and date: ' read *,aiey,aem,aed iey=aiey em=aem ed=aed if (iey.ge.0.and.iey.lt.50) then iey=iey+2000 else if (iey.ge.50.and.iey.le.99) then iey=iey+1900 end if c if (iey.gt.1900) iyr=iyr-1900 print '(/,a,$)',' no of days to be averaged: ' read *,nad if (nad.eq.30) mnthly='y' !comment if actual 30 day avrg rqrd end if if (gm.eq.'w') then print '(/,a,$)',' wtng pwr: ' read *,ipwr print '(/,a,$)',' srch rds: ' read *,srad end if if (ymd.ne.'n') then idfy=iey-isy stnd=idys(ism-1)+isd irm=mod(isy,4) c irmdr=isy-irm*4 if (irm.eq.0) stnd=idysl(ism-1)+isd endy=idys(em-1)+ed+iyrs(idfy) if (irm.eq.0) endy=idysl(em-1)+ed+iyrs(idfy) ndt=(endy-stnd)/nad end if nx=int((xmx-xmn)/dx)+1. ny=int((ymx-ymn)/dy)+1. c initialise x & y grid mode: do i=1,ny gy(i)=float(i-1)*dy+ymn-dy/2 end do do i=1,nx gx(i)=float(i-1)*dx+xmn-dx/2 end do if (nx.gt.m1.or.ny.gt.m2.or.ndt.gt.m3) then print*, nx, m1 print*, ny, m2 print*, ndt, m3 print *, 'Dimension problem in sumbz' stop end if do i=1,nx do j=1,ny do k=1,ndt sumz(i,j,k)=0. sumw(i,j,k)=0. end do end do end do 12 continue read (1,*) 5 read(1,*,end=90)(parm(i),i=1,np) nn=nn+1 y=parm(lac) x=parm(loc) z=parm(cn) if (ymd.ne.'n') then iyr=parm(iy) c iyr=2002 !%%%if all years of 10-day (for eg) avrg is rqrd if (iyr.ge.0.and.iyr.lt.50) then iyr=iyr+2000 else if (iyr.ge.50.and.iyr.le.99) then iyr=iyr+1900 end if c if (iyr.gt.1900) iyr=iyr-1900 if (iyr.gt.iey) go to 90 imon=parm(mon) idy=parm(idt) idffy=iyr-isy xxx=idys(imon-1)+idy irmm=mod(iyr,4) c irmmdr=iyr-irmm*4 if (irmm.eq.0) xxx=idysl(imon-1)+idy end if if (x.lt.xmn.or.x.gt.xmx.or.y.lt.ymn.or.y.gt.ymx.or.z.eq.thv) 1 go to 5 if (ymd.ne.'n') then !********** if (xxx.gt.endy) go to 5 if (mnthly.eq.'y') then nad=idys(imon)-idys(imon-1) !no change reqd? if (irmm.eq.0)nad=idysl(imon)-idysl(imon-1) !no change reqd? nd=imon-ism+1+idffy*12 nd1=imon else itdys=iyrs(iyr-1950)-iyrs(isy-1950) c print *, itdys,iyrs(iyr-1950),iyrs(isy-1950),iyr,isy nd=(idys(imon-1)+idy+itdys-stnd+nad)/nad c nd=(idys(imon-1)+idy+iyrs(idffy)-stnd+nad)/nad if (irmm.eq.0) nd=(idysl(imon-1)+idy+itdys-stnd+nad)/nad end if c print *,nd if (nd.le.0) go to 5 xiy(nd)=iyr xim(nd)=imon xid(nd)=idy mid=min(mid,nd) mad=max(mad,nd) end if !*********** ncd=ncd+1 if (ymd.eq.'n') then mid=1 mad=1 nd=1 end if if (gm.eq.'w') then xmnlmt=int((x-xmn-srad)/dx)-5 ymnlmt=int((y-ymn-srad)/dy)-5 xmxlmt=int((x-xmn+srad)/dx)+5 ymxlmt=int((y-ymn+srad)/dy)+5 end if if (gm.eq.'b') then xmnlmt=int((x-xmn-dx)/dx)-1 ymnlmt=int((y-ymn-dy)/dy)-1 xmxlmt=int((x-xmn+dx)/dx)+1 ymxlmt=int((y-ymn+dy)/dy)+1 end if if (xmnlmt.lt.0) xmnlmt=1 !0 if (ymnlmt.lt.0) ymnlmt=1 !0 if (xmxlmt.gt.nx) xmxlmt=nx if (ymxlmt.gt.ny) ymxlmt=ny do i=ymnlmt,ymxlmt do j=xmnlmt,xmxlmt dx1=abs(x-gx(j)) dy1=abs(y-gy(i)) dist=sqrt(dx1*dx1+dy1*dy1) cc print '(7f8.2)',dx1,x,gx(j),dy1,y,gy(i),dist if (dist.lt.srad.and.gm.eq.'w') then c weighting functions: c 1: gaussian 2: Cressman 3: inverse square 4: Barnes wf=2 if (wf.eq.1) then wfnc=exp(-1*(dist**ipwr)/(2*srad**ipwr)) else if (wf.eq.2) then wfnc=(-dist**ipwr+srad**ipwr)/(dist**ipwr+srad**ipwr) else if (wf.eq.3) then wfnc=1./(1.+dist**ipwr) else if (wf.eq.4) then wfnc=exp(-4*(dist**ipwr)/srad**ipwr) end if c sumz(i,j,nd)=sumz(i,j,nd)+(z/(dist**ipwr+1.)) c sumw(i,j,nd)=sumw(i,j,nd)+(1./(dist**ipwr+1.)) sumz(i,j,nd)=sumz(i,j,nd)+(z*wfnc) cc std dev cannot be computed with weighted average cc sig x*x is less for weighted avrg cc stdz(i,j,nd)=stdz(i,j,nd)+(z*wfnc)*(z*wfnc) sumw(i,j,nd)=sumw(i,j,nd)+wfnc else if (dx1.lt.dx.and.x.lt.gx(j)+dx/2..and. 1 dy1.lt.dy.and.y.lt.gy(i)+dy/2..and.gm.eq.'b') then sumz(i,j,nd)=sumz(i,j,nd)+z ! (z*wfnc) sumw(i,j,nd)=sumw(i,j,nd)+1 ! (wfnc) stdz(i,j,nd)=stdz(i,j,nd)+z*z end if end do end do go to 5 90 nff=nff+1 close (1) if (nff.lt.nf) go to 14 print '(///,a,$)', ' nx,ny: ' print *,nx,ny c write (12,'('' total/considered values: '',2i10)')nn,ncd print '(//)' print *,' total/considered values: ',nn,ncd tmpcnv='n' if (tmpcnv.eq.'y') then print * print * print *,' Temp conv in active' print * print * print *, ' tmpmn, tmpmx, shmn, shmx' read *,tmpmn,tmpmx,shmn,shmx end if if (mnthly.eq.'y') nad=30 do nd=mid,mad iyy=xiy(nd) mnth=xim(nd) idate=xid(nd) do i=1,ny rmdlt=gy(i) !-dy/2 !remove dy/2 if cntr pnt not rqrd nh=0 do j=1,nx if (sumw(i,j,nd).gt.0) then sumz(i,j,nd)=sumz(i,j,nd)/sumw(i,j,nd) if (tmpcnv.eq.'y') then tmdf=tmpmx-tmpmn shdf=shmx-shmn if (sumz(i,j,nd).gt.shmn.and.sumz(i,j,k).lt.shmx) then tmp=(sumz(i,j,nd)-shmn)*tmdf/shdf else if(sumz(i,j,k).le.shmn) tmp=0 if(sumz(i,j,k).gt.shmx) tmp=tmdf+1 end if sumz(i,j,nd)=tmp+tmpmn !ssh to T conversion end if if (std.eq.'y') then if (sumw(i,j,nd).ge.3) then xx=stdz(i,j,nd)/sumw(i,j,nd) sumz(i,j,nd)=sqrt(xx-sumz(i,j,nd)*sumz(i,j,nd)) else sumz(i,j,nd)=vna end if end if else sumz(i,j,nd)=vna end if if (sumz(i,j,nd).gt.zmx.and.sumw(i,j,nd).gt.0)zmx=sumz(i,j,nd) if (sumz(i,j,nd).lt.zmn.and.sumw(i,j,nd).gt.0)zmn=sumz(i,j,nd) iydf=(iyy-isy)*12 rmdln=gx(j) !-dx/2 yn=ymn-dy/2. xn=xmn-dx/2. if (ymd.ne.'n') then !wth yr,mn,dt c writes only at data available points if (sumz(i,j,nd).ne.vna.and.iyy.ne.0.and.rmdlt x .ge.yn.and.rmdln.ge.xn.and. x rmdln.ge.xmn.and.rmdlt.ge.ymn.and.all.eq.'n') then write(12,'(3i5,3f9.2)') x iyy,mnth,idate,rmdlt,rmdln,sumz(i,j,nd) nh=nh+1 hrzntl(nh)=sumz(i,j,nd) ! for writing in map form (horizontally) end if c writes at all the points c if (iyy.ne.0.and.rmdlt.ge.yn.and.rmdln.ge.xn.and. if (rmdlt.ge.yn.and.rmdln.ge.xn.and. x rmdln.ge.xmn.and.rmdlt.ge.ymn.and.all.eq.'y') then write (12,'(3i5,3f9.2)') x iyy,mnth,idate,rmdlt,rmdln,sumz(i,j,nd) nh=nh+1 hrzntl(nh)=sumz(i,j,nd) ! for writing in map form (horizontally) end if end if if (ymd.eq.'n') then !w/o yr,mn,dt if (sumz(i,j,nd).ne.vna.and.rmdlt x .ge.yn.and.rmdln.ge.xn.and.rmdln.ge.xmn. x and.rmdlt.ge.ymn.and.all.eq.'n') write(12,'(3f9.2)') x rmdlt,rmdln,sumz(i,j,nd) if (rmdlt.ge.yn.and.rmdln.ge.xn.and. x rmdln.ge.xmn.and.rmdlt.ge.ymn.and.all.eq.'y') x write (12,'(3f9.2)') x rmdlt,rmdln,sumz(i,j,nd) end if end do c if (rmdlt.ge.yn.and.rmdln.ge.xn.and. c x rmdln.ge.xmn.and.rmdlt.ge.ymn.and.all.eq.'y') then if (nh.ne.0) write (110,'(20f9.2)')(hrzntl(j),j=1,nh) c end if end do end do end