subroutine obs_pos_interpol !*********************************************************************** ! ! Gives 4D (space and time) interpolated output at position and time of ! measurements (flight tracks, soundings, etc.) ! The routine assumes a regular model grid in both lon and lat direction. ! For other grids (e.g. Gaussian) modifications will be necessary ! to that part of the routine which determines the grid cell containing ! the measurement point and the 7 closest neighbour cells. ! ! The routine is called at the end of each model time step. It has to ! be initialized at the beginning of every new month (Here it is ! assumed that each model month is simulated separately. Therefore ! the routine is only initialized in its first call). ! ! Note: separate interpolations for H2O and PO3 are not implemented here. ! ! The initialization does: ! 1) Open the observation position file. ! 2) Read until first observation AFTER the present model time ! 3) Save tracers fields of the selected tracers for time interpolation. ! ! The program is then called after each model time step until the ! end of the month to do the time-interpolation between the present ! and the previous model fields and to do the spatial interpolation ! to all observation positions available between the present and the ! previous model time (ie. observations in the interval [t-timestep,t]) ! ! Interpolation steps: ! 1) the cell containing the measurement is determined ! 2) the seven closest neighbour cells are determined (together ! they make up a cube with 8 edge points) ! 3) linear interpolation in time (between 8 points in previous ! and present model field) and in space (trilinear interpolation ! between the 8 points). Note: vertical interpolation is done ! linearly in log(pressure). ! ! Expected input is the TRADEOFF observation file 'obs_input.dat'. ! The obs_input.dat file has to be overwritten with the appropriate file ! at the beginning of each month simulated. ! ! NOTE1: ! All variables which should be saved until the next ! call are defined with the SAVE attribute. ! NOTE2: Time interpolation could probably be omitted. In this case ! however, all observations in the time interval ! [t-1/2*timestep,itau+0.5*timestep] should be loaded. ! ! 31 July '00, Dominik Brunner, ETH Zurich, Switzerland ! based on subroutine falcon by Ernst Meijer, KNMI, The Netherlands ! !*********************************************************************** use inits ! defines grid geometry (nx,ny,nz,lonmin,latmin,deltalon ! deltalat,phlb), the model time variables (idate,itau), ! the Mole masses xm and xmair, and the mass mixing ratio ! field mmr use dateconv ! date/time conversion routines idate2itau and itau2idate ! which convert between date (YYYY,MM,DD,hh,mm,ss) and ! model time in seconds. implicit none integer, parameter :: no=10 ! number of output species integer, parameter, dimension(no) :: is= & ! indices of output parameters (/io3,ico,ioh,iho2,ih2o2,ino,ino2,ihno3,ipan,iradon/) real, dimension(no) :: fout=-9999. ! output array (interpolated values (ppbv)) integer, save, dimension(6) :: tidate=0 ! observation date/time YYYY,MM,DD,hh,mm real, save, dimension(3) :: txyz ! obs. position (lon,lat,press) integer, save :: titau ! obs. date/time in seconds integer :: tid ! obs. data source ID integer, dimension(3) :: tijl ! obs. position in terms of model indices (i,j,l) real :: tlon,tlat,tplog ! longitude, latitude and log(pressure) ! of cell tijl containing the obs. pos. real, dimension(2) :: vxtmp,vytmp,vztmp ! temporary values for the sequential inter- ! polation in time, vert, lat and lon direction real, dimension(2) :: cx,cy,cz ! lon, lat and log press of cube of neighbour cells integer, dimension(2) :: ci,cj,cl ! indices of cube of neighbour cells real :: ftint ! factor for time interpolation real, dimension(2) :: fhint ! factors for horizontal interpolation real :: fvint ! factor for vertical interolation real :: utc ! time in seconds of day (UTC) real, save, dimension(nx,ny,nz,no) :: mmro ! old mass mixing ratios (MMR) ! real, save, dimension(nx,ny,nz,no) :: h2oo ! old H2O mass mixing ratios ! real, save, dimension(nx,ny,nz,no) :: po3o ! old O3 production rates real, save, dimension(2) :: xat ! current and old model time (as real number) integer, parameter :: kfin=43 ! unit number of observation data file character(len=13), parameter :: obsinname ='obs_input.dat' ! input filename character(len=14), parameter :: obsoutname='obs_output.dat' ! output filename integer, parameter :: kfo=44 ! unit number of model interpolation file logical, save :: fex ! flag for existing observation file logical, save :: first ! flag for initialization character(len=49) :: chead ! header line real, save, dimension(no) :: ppbf ! conversion factors MMR -> ppb integer :: i,j,l,n ! counters data fex /.true./ data first /.true./ ! INITIALISATION if (first) then ! init conversion factors do n=1,no ppbf(n)=1.e9*xmair/xm(is(n)) enddo first = .false. inquire (file=obsinname,exist=fex) if (.not.fex) return open(unit=kfin,file=obsinname,status='old') open(unit=kfo,file=obsoutname) ! read header line read(kfin,95)chead write(*,*) 'Observation interpolation - Initialisation' write(*,*)chead 95 format(a43) ! output header: write(kfo,*)'YYYYMMDD secs lon(deg) lat(deg) p(Pa) ID',(' '//names(is(n)),n=1,no) ! read until first position after model time titau = -99. do while(titau.le.itau) ! first obs. point must be after model time read(kfin,96,end=98) tidate(1:5),txyz(1),txyz(2),txyz(3),tid ! lon,lat,press call idate2itau(tidate,titau) ! write out dummy (-9999.) values because no interpolation is made ! at beginning of month before first model time step! if (titau.le.itau) then write(kfo,97)10000*tidate(1)+100*tidate(2)+tidate(3),& 3600*tidate(4)+60*tidate(5)+tidate(6),& txyz(1),txyz(2),txyz(3),tid,(fout(n),n=1,no) endif enddo !96 format(i4,i2,i2,i3,i2,f10.4,f9.4,f7.0,i3) 96 format(i4,i2,i2,i3,i2,f10.4,f9.4,f7.0,i6) ! 97 format(i8,1x,i5,f10.4,f9.4,f8.0,i3,50e12.4) 97 format(i8,1x,i5,f10.4,f9.4,f8.0,i6,50e12.4) goto 99 98 continue ! no more data points in file fex=.false. close(kfin) close(kfo) 99 continue endif xat(2)=float(itau) ! save actual time as real number ! AS LONG AS DATA POINTS ARE BEFORE MODEL TIME; INTERPOLATE do while (titau.le.itau.and.fex) ! DETERMINE INDICES OF CELL CONTAINING MEASUREMENT ! Calculation based on regular grid as defined in inits.f90 ! first horizontal indices tijl(1) = MOD(int((txyz(1)-(lonmin-deltalon/2.))/deltalon)+1,nx) tijl(2) = MOD(int((txyz(2)-(latmin-deltalat/2.))/deltalat)+1,ny) ! second, vertical index tijl(3)=1 do l=1,nz if ( txyz(3).le.phlb(tijl(1),tijl(2),l).and. & txyz(3).gt.phlb(tijl(1),tijl(2),l+1)) tijl(3) = l enddo ! AVOID INTERPOLATION BELOW SURFACE: If pressure txyz(3) exceeds ! pressure at center of lowest model layer than set txyz(3) ! to the pressure of the lowest layer minus 1 Pa. ! (-1 Pa to ensure that no interpolation is made below the surface) if (txyz(3).gt.(phlb(tijl(1),tijl(2),2)+phlb(tijl(1),tijl(2),1))/2.) then tijl(3)=1 ! assign pressure at layer center txyz(3)=0.5*(phlb(tijl(1),tijl(2),2)+phlb(tijl(1),tijl(2),1))-0.1 endif ! AVOID INTERPOLATION ABOVE MODEL TOP: If pressure txyz(3) exceeds ! pressure at center of highest model layer than set txyz(3) ! to the pressure of the highest layer plus 1 Pa. ! (+1 Pa to ensure that no interpolation is made above top) if (txyz(3).lt.(phlb(tijl(1),tijl(2),nz+1)+phlb(tijl(1),tijl(2),nz))/2.) then tijl(3)=nz ! assign pressure at layer center txyz(3)=0.5*(phlb(tijl(1),tijl(2),nz+1)+phlb(tijl(1),tijl(2),nz))+0.1 endif ! DETERMINE INDICES AND LON,LAT,LOG(PRESS) VALUES OF CELL CENTERS AND ! NEIGHBOURING CELL CENTERS ! Calculation based on regular grid as defined in inits.f90 txyz(3)=log(txyz(3)) ! convert to log(pressure) for interpolation tlon = lonmin+float(tijl(1)-1)*deltalon tlat = latmin+float(tijl(2)-1)*deltalat tplog = log(0.5*(phlb(tijl(1),tijl(2),tijl(3))+ & phlb(tijl(1),tijl(2),tijl(3)+1))) if (tlon.le.txyz(1)) then ci(1) = tijl(1) ci(2) = tijl(1)+1 cx(1) = tlon cx(2) = tlon + deltalon else ci(1) = tijl(1)-1 ci(2) = tijl(1) cx(1) = tlon - deltalon cx(2) = tlon endif if (tlat.le.txyz(2)) then cj(1) = tijl(2) cj(2) = tijl(2)+1 cy(1) = tlat cy(2) = tlat + deltalat else cj(1) = tijl(2)-1 cj(2) = tijl(2) cy(1) = tlat - deltalat cy(2) = tlat endif if (tplog.ge.txyz(3)) then cl(1) = tijl(3) cl(2) = tijl(3)+1 cz(1) = tplog cz(2) = log(0.5*(phlb(tijl(1),tijl(2),cl(2))+ & phlb(tijl(1),tijl(2),cl(2)+1))) else cl(1) = tijl(3)-1 cl(2) = tijl(3) cz(1) = log(0.5*(phlb(tijl(1),tijl(2),cl(1))+ & phlb(tijl(1),tijl(2),cl(1)+1))) cz(2) = tplog endif ! CALCULATE FACTORS FOR LINEAR INTERPOLATION IN 4D (TIME AND SPACE) ftint=(float(titau)-xat(1))/(xat(2)-xat(1)) ! time interpolation factor ! horizontal interpolation factors for lon and lat direction fhint(1)=(txyz(1)-cx(1))/(cx(2)-cx(1)) fhint(2)=(txyz(2)-cy(1))/(cy(2)-cy(1)) ! vertical interpolation factor fvint=(txyz(3)-cz(1))/(cz(2)-cz(1)) ! LOOP OVER TRACERS AND INTERPOLATE BETWEEN CELL AND CLOSEST NEIGHBOURS do n=1,no do i=1,2 ! lon do j=1,2 ! lat do l=1,2 ! vertical ! do time interpolation vztmp(l)=mmro(ci(i),cj(j),cl(l),n)+ftint*& ( mmr(ci(i),cj(j),cl(l),is(n))-mmro(ci(i),cj(j),cl(l),n) ) enddo ! do vertical interpolation vytmp(j)=vztmp(1)+fvint*(vztmp(2)-vztmp(1)) enddo ! do lat interpolation vxtmp(i)=vytmp(1)+fhint(2)*(vytmp(2)-vytmp(1)) enddo ! do lon interpolation and multiply with conversion factor fout(n)=ppbf(n)*(vxtmp(1)+fhint(1)*(vxtmp(2)-vxtmp(1))) enddo ! DO INTERPOLATION ALSO FOR OZONE PRODUCTION RATE AND H2O (NOT IMPLEMENTED) ! OUTPUT LABELED WITH DATE, UTC AND POSITION (LON,LAT,PRES) write(kfo,100)10000*tidate(1)+100*tidate(2)+tidate(3),& 3600*tidate(4)+60*tidate(5)+tidate(6),& txyz(1),txyz(2),exp(txyz(3)),tid,(fout(n),n=1,no) 100 format(i8,1x,i5,f10.4,f9.4,f8.0,i3,50e12.4) ! PREPARE OUTPUT AT NEXT POSITION read(kfin,96,end=101) tidate(1:5),txyz(1),txyz(2),txyz(3),tid ! lon,lat,press call idate2itau(tidate,titau) goto 102 101 continue ! no more data points in file fex=.false. close(kfin) close(kfo) 102 continue enddo ! SAVE THE TIME AND MASS MIXING RATIOS FOR NEXT CALL xat(1)=float(itau) do n=1,no mmro(:,:,:,n) = mmr(:,:,:,is(n)) enddo ! SAVE H2O MMR AND OZONE PRODUCTION RATE ! h2oo=h2o ! po3o=po3 return end subroutine obs_pos_interpol