| 0,0 → 1,723 |
|---|
| PROGRAM inv_prep |
| |
| c ******************************************************************************** |
| c * CALCULATE REFERENCE PROFILE, CORIOLIS PARAMETER AND GRID PARAMETERS * |
| c * Rene Fehlmann 1994 / Code re-organization: Michael Sprenger, 2006 * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Declaration of variables, parameters, externals and common blocks |
| c -------------------------------------------------------------------------------- |
| |
| implicit none |
| |
| c Input and output file |
| character*80 pvsrcfile,referfile |
| integer mode |
| real radius |
| |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real mdv |
| |
| c Reference state |
| real, allocatable,dimension (:) :: nsqref |
| real, allocatable,dimension (:) :: thetaref |
| real, allocatable,dimension (:) :: rhoref |
| real, allocatable,dimension (:) :: pressref |
| real, allocatable,dimension (:) :: zref |
| real pressn,thetan |
| |
| c 3d fields for calculation of reference profile |
| real,allocatable,dimension (:,:,:) :: th,rho,nsq,p,z |
| |
| c 2d weight for mean |
| real,allocatable,dimension (:,:) :: weight |
| |
| c Auxiliary variables |
| integer i,j,k |
| integer stat |
| integer jj,kk |
| character*80 varname |
| integer istep |
| real sum,max,min |
| integer cnt,nmd |
| integer step |
| integer i0,j0 |
| real lon0,lat0,lon1,lat1,dist |
| integer vardim(4),ndim,cdfid,ierr |
| real varmin(4),varmax(4),stag(4) |
| real mdv |
| character*80 name |
| |
| c Externals |
| real sdis |
| external sdis |
| |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| |
| print*,'********************************************************' |
| print*,'* ref_grid *' |
| print*,'********************************************************' |
| |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) pvsrcfile |
| read(10,*) referfile |
| read(10,*) name,radius |
| if ( trim(name).ne.'REF_R') stop |
| close(10) |
| print* |
| print*,trim(pvsrcfile) |
| print*,trim(referfile) |
| print*,radius |
| print* |
| |
| c Get lat/lon gid parameters from input file |
| call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > pvsrcfile) |
| print*,'Read_Dim: nx,ny,nz = ',nx,ny,nz |
| print*,' dx,dy,dz = ',dx,dy,dz |
| print*,' xmin,ymin,zmin = ',xmin,ymin,zmin |
| print*,' mdv = ',mdv |
| print* |
| xmax = xmin + real(nx-1) * dx |
| ymax = ymin + real(ny-1) * dy |
| |
| c Count from 0, not from 1: consistent with <inv_cart.f>. |
| nx=nx-1 |
| ny=ny-1 |
| nz=nz-1 |
| |
| c Allocate memory for reference profile |
| allocate(rhoref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating rhoref' |
| allocate(pressref(0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pressref' |
| allocate(thetaref(0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating thetaref' |
| allocate(nsqref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating nsqref' |
| allocate(zref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating zref' |
| |
| c Allocate memory for calculatation of reference profile |
| allocate(th (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating th' |
| allocate(rho(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating rho' |
| allocate(p (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating p' |
| allocate(nsq(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating nsq' |
| allocate(z(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating z' |
| |
| c Allocate memory for weight |
| allocate(weight(0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating weight' |
| |
| c -------------------------------------------------------------------------------- |
| c Calculate the reference profile and put it onto netcdf file |
| c -------------------------------------------------------------------------------- |
| |
| c Read data from file |
| varname='TH' |
| call read_inp (th,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='RHO' |
| call read_inp (rho,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='NSQ' |
| call read_inp (nsq,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='P' |
| call read_inp (p,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| |
| c Init the height field (not really necessary, but code becomes more symmetrical) |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| z(i,j,k)=zmin+real(k)*dz |
| enddo |
| enddo |
| enddo |
| |
| c Define the weight |
| if ( radius.lt.0 ) then |
| do i=0,nx |
| do j=0,ny |
| weight(i,j) = 1. |
| enddo |
| enddo |
| else |
| i0 = nx/2 |
| j0 = ny/2 |
| lon0 = xmin + real(i0-1) * dx |
| lat0 = ymin + real(j0-1) * dy |
| weight(i0,j0)=1. |
| do i=0,nx |
| do j=0,ny |
| lon1 = xmin + real(i-1) * dx |
| lat1 = ymin + real(j-1) * dy |
| dist = sdis(lon0,lat0,lon1,lat1) |
| if ( dist.lt.radius ) then |
| weight(i,j) = 1. |
| else |
| weight(i,j) = 0. |
| endif |
| enddo |
| enddo |
| endif |
| |
| c Determine the reference profile (mean over domain, split levels and layers) |
| call mean(zref, z, nx,ny,nz,mdv,weight) |
| call mean(nsqref, nsq,nx,ny,nz,mdv,weight) |
| call mean(rhoref, rho,nx,ny,nz,mdv,weight) |
| call mean(thetaref,th, nx,ny,nz,mdv,weight) |
| call mean(pressref,p, nx,ny,nz,mdv,weight) |
| |
| c Write reference file to netcdf file |
| call write_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nz,referfile) |
| |
| c Write some info |
| print* |
| print*,'Ref: z p rho nsq theta' |
| step=2*nz/10 |
| if (step.lt.1) step=1 |
| do k=0,2*nz,step |
| write(*,'(8x,f10.1,2f10.2,f10.6,f10.2)') |
| > zref(k),pressref(k),rhoref(k),nsqref(k),thetaref(k) |
| enddo |
| print* |
| |
| c Write weighht to REF file |
| call cdfwopn(trim(referfile),cdfid,ierr) |
| varname='WEIGHT' |
| vardim(1)=nx+1 |
| vardim(2)=ny+1 |
| vardim(3)=1 |
| vardim(4)=1 |
| varmin(1)=xmin |
| varmin(2)=ymin |
| varmin(3)=0. |
| varmax(1)=xmax |
| varmax(2)=ymax |
| varmax(3)=0. |
| ndim=3 |
| mdv=-999. |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| call putdat(cdfid,varname,0.,1,weight,ierr) |
| |
| c -------------------------------------------------------------------------------- |
| c Format specifications |
| c -------------------------------------------------------------------------------- |
| |
| 111 format (5f20.9) |
| 106 format (2f20.9) |
| |
| end |
| |
| c -------------------------------------------------------------------------- |
| c Spherical distance between lat/lon points |
| c -------------------------------------------------------------------------- |
| |
| real function sdis(xp,yp,xq,yq) |
| c |
| c calculates spherical distance (in km) between two points given |
| c by their spherical coordinates (xp,yp) and (xq,yq), respectively. |
| c |
| real re |
| parameter (re=6370.) |
| real pi180 |
| parameter (pi180=3.14159/180.) |
| real xp,yp,xq,yq,arg |
| |
| arg=sin(pi180*yp)*sin(pi180*yq)+ |
| > cos(pi180*yp)*cos(pi180*yq)*cos(pi180*(xp-xq)) |
| if (arg.lt.-1.) arg=-1. |
| if (arg.gt.1.) arg=1. |
| |
| sdis=re*acos(arg) |
| |
| end |
| |
| c ******************************************************************************** |
| c * NETCDF INPUT AND OUTPUT * |
| c ******************************************************************************** |
| |
| |
| c -------------------------------------------------------------------------------- |
| c Write reference file to netcdf |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE write_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nz,pvsrcfile) |
| |
| c Write the reference profile to file |
| c |
| c thetaref : Reference potential temperature (K) |
| c pressref : Reference pressure (Pa) |
| c rhoref : Reference density (kg/m^3) |
| c nsqref : Stratification (s^-1) |
| c zref : Reference height (m) |
| c nz : Grid dimension in z direction |
| c pvsrcfile : Output file |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real nsqref (0:2*nz) |
| real thetaref(0:2*nz) |
| real rhoref (0:2*nz) |
| real pressref(0:2*nz) |
| real zref (0:2*nz) |
| character*80 pvsrcfile |
| |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer cdfid,stat |
| integer vardim(4) |
| real misdat |
| integer ndimin |
| real varmin(4),varmax(4),stag(4) |
| integer i,j,k,nf1 |
| integer ntimes |
| real time(200) |
| character*80 vnam(100),varname |
| integer nvars |
| integer isok,ierr |
| |
| c Get grid description from topography |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 997 |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='ORO' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 997 |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Open output netcdf file |
| call cdfwopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Create the variable if necessary |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='NSQREF' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| vardim(1)=1 |
| vardim(2)=1 |
| vardim(3)=2*nz+1 |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| endif |
| isok=0 |
| varname='RHOREF' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| vardim(1)=1 |
| vardim(2)=1 |
| vardim(3)=2*nz+1 |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| endif |
| isok=0 |
| varname='PREREF' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| vardim(1)=1 |
| vardim(2)=1 |
| vardim(3)=2*nz+1 |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| endif |
| isok=0 |
| varname='THETAREF' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| vardim(1)=1 |
| vardim(2)=1 |
| vardim(3)=2*nz+1 |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| endif |
| isok=0 |
| varname='ZREF' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| vardim(1)=1 |
| vardim(2)=1 |
| vardim(3)=2*nz+1 |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| endif |
| |
| c Write data |
| varname='NSQREF' |
| print*,'W NSQREF ',trim(pvsrcfile) |
| call putdat(cdfid,varname,time(1),0,nsqref,stat) |
| if (stat.ne.0) goto 997 |
| varname='RHOREF' |
| print*,'W RHOREF ',trim(pvsrcfile) |
| call putdat(cdfid,varname,time(1),0,rhoref,stat) |
| if (stat.ne.0) goto 997 |
| varname='THETAREF' |
| print*,'W THETAREF ',trim(pvsrcfile) |
| call putdat(cdfid,varname,time(1),0,thetaref,stat) |
| if (stat.ne.0) goto 997 |
| varname='PREREF' |
| print*,'W PREREF ',trim(pvsrcfile) |
| call putdat(cdfid,varname,time(1),0,pressref,stat) |
| if (stat.ne.0) goto 997 |
| varname='ZREF' |
| print*,'W ZREF ',trim(pvsrcfile) |
| call putdat(cdfid,varname,time(1),0,zref,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Close output netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| |
| return |
| |
| c Exception handling |
| 997 print*,'Write_Ref: Problem with input netcdf file... Stop' |
| stop |
| |
| end |
| |
| |
| c -------------------------------------------------------------------------------- |
| c Read input fields for reference profile |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE read_inp (field,fieldname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| |
| c Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified |
| c by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input |
| c files are consitent with this grid. The missing data value is set to <mdv>. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| character*80 fieldname |
| character*80 pvsrcfile |
| real dx,dy,dz,xmin,ymin,zmin |
| real mdv |
| |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer cdfid,stat,cdfid99 |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| integer ndimin,outid,i,j,k |
| real max_th |
| real tmp(nx,ny,nz) |
| integer ntimes |
| real time(200) |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| |
| c Open the input netcdf file |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Check whether needed variables are on file |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname=trim(fieldname) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| |
| c Get the grid parameters |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Check whether grid parameters are consistent |
| if ( (vardim(1).ne.(nx+1)).or. |
| > (vardim(2).ne.(ny+1)).or. |
| > (vardim(3).ne.(nz+1)).or. |
| > (abs(varmin(1)-xmin).gt.eps).or. |
| > (abs(varmin(2)-ymin).gt.eps).or. |
| > (abs(varmin(3)-zmin).gt.eps).or. |
| > (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or. |
| > (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or. |
| > (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) |
| >then |
| print*,'Input grid inconsitency...' |
| print*,' Nx = ',vardim(1),nx+1 |
| print*,' Ny = ',vardim(2),ny+1 |
| print*,' Nz = ',vardim(3),nz+1 |
| print*,' Varminx = ',varmin(1),xmin |
| print*,' Varminy = ',varmin(2),ymin |
| print*,' Varminz = ',varmin(3),zmin |
| print*,' Dx = ',(varmax(1)-varmin(1))/real(nx-1),dx |
| print*,' Dy = ',(varmax(2)-varmin(2))/real(ny-1),dy |
| print*,' Dz = ',(varmax(3)-varmin(3))/real(nz-1),dz |
| goto 998 |
| endif |
| |
| c Load variables |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| call getdat(cdfid,varname,time(1),0,field,stat) |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| |
| c Close input netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Set missing data value to <mdv> |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if (abs(field(i,j,k)-misdat).lt.eps) then |
| field(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| |
| return |
| |
| c Exception handling |
| 998 print*,'Read_Inp: Problem with input netcdf file... Stop' |
| stop |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Check whether variable is found on netcdf file |
| c -------------------------------------------------------------------------------- |
| |
| subroutine check_varok (isok,varname,varlist,nvars) |
| |
| c Check whether the variable <varname> is in the list <varlist(nvars)>. If this is |
| C the case, <isok> is incremented by 1. Otherwise <isok> keeps its value. |
| |
| implicit none |
| |
| c Declaraion of subroutine parameters |
| integer isok |
| integer nvars |
| character*80 varname |
| character*80 varlist(nvars) |
| |
| c Auxiliary variables |
| integer i |
| |
| c Main |
| do i=1,nvars |
| if (trim(varname).eq.trim(varlist(i))) isok=isok+1 |
| enddo |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Get grid parameters |
| c -------------------------------------------------------------------------------- |
| |
| subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > pvsrcfile) |
| |
| c Get the grid parameters from the variable <TH> on the input file <pvsrcfile>. |
| c The grid parameters are |
| c nx,ny,nz : Number of grid points in x, y and z direction |
| c xmin,ymin,zmin : Minimum domain coordinates in x, y and z direction |
| c xmax,ymax,zmax : Maximal domain coordinates in x, y and z direction |
| c dx,dy,dz : Horizontal and vertical resolution |
| c Additionally, it is checked whether the vertical grid is equally spaced. If ok, |
| c the grid paramters are transformed from lon/lat to distance (in meters) |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| character*80 pvsrcfile |
| integer nx,ny,nz |
| real dx,dy,dz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real mdv |
| |
| c Numerical epsilon and other physical/geoemtrical parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer cdfid,cstid |
| integer ierr |
| character*80 vnam(100),varname |
| integer nvars |
| integer isok |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| real aklev(1000),bklev(1000),aklay(1000),bklay(1000) |
| real dh |
| character*80 csn |
| integer ndim |
| integer i |
| |
| c Get all grid parameters |
| call cdfopn(pvsrcfile,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname='T' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| call getcfn(cdfid,csn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(csn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax, |
| > stag,ierr) |
| if (ierr.ne.0) goto 998 |
| nx=vardim(1) |
| ny=vardim(2) |
| nz=vardim(3) |
| xmin=varmin(1) |
| ymin=varmin(2) |
| zmin=varmin(3) |
| call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr) |
| if (ierr.ne.0) goto 998 |
| call getgrid(cstid,dx,dy,ierr) |
| if (ierr.ne.0) goto 998 |
| xmax=varmax(1) |
| ymax=varmax(2) |
| zmax=varmax(3) |
| dz=(zmax-zmin)/real(nz-1) |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Check whether the grid is equallay spaced in the vertical |
| do i=1,nz-1 |
| dh=aklev(i+1)-aklev(i) |
| if (abs(dh-dz).gt.eps) then |
| print*,'Aklev: Vertical grid must be equally spaced... Stop' |
| print*,(aklev(i),i=1,nz) |
| stop |
| endif |
| dh=aklay(i+1)-aklay(i) |
| if (abs(dh-dz).gt.eps) then |
| print*,'Aklay: Vertical grid must be equally spaced... Stop' |
| print*,(aklay(i),i=1,nz) |
| stop |
| endif |
| enddo |
| |
| c Set missing data value |
| mdv=misdat |
| |
| return |
| |
| c Exception handling |
| 998 print*,'Read_Dim: Problem with input netcdf file... Stop' |
| stop |
| |
| end |
| |
| |
| c ******************************************************************************** |
| c * DEFINE REFERENCE PROFILE * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Calculate area mean |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE mean(a,m,nx,ny,nz,mdv,weight) |
| |
| c Calculate the area-mean of <m> and save it on <a>. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| real mdv |
| integer nx,ny,nz |
| real m(0:nx,0:ny,0:nz),a(0:2*nz) |
| real weight(0:nx,0:ny) |
| |
| c Numerical epsilon |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary varaibles |
| real mea(0:nz) |
| integer i,j,k,kk |
| real counter |
| |
| c Determine the mean over all levels (handle missing data) |
| do k=0,nz |
| mea(k)=0. |
| counter=0. |
| do i=0,nx |
| do j=0,ny |
| if (abs(m(i,j,k)-mdv).gt.eps) then |
| mea(k)=mea(k)+m(i,j,k)*weight(i,j) |
| counter=counter+weight(i,j) |
| endif |
| |
| enddo |
| enddo |
| if (counter.gt.0) then |
| mea(k)=mea(k)/counter |
| else |
| mea(k)=mdv |
| endif |
| enddo |
| |
| c Prepare the output array: split layers and levels |
| do k=0,nz-1 |
| kk=2*k |
| a(kk)=mea(k) |
| a(kk+1)=0.5*(mea(k)+mea(k+1)) |
| enddo |
| a(2*nz)=mea(nz) |
| |
| end |
| Property changes: |
|---|
| Added: svn:executable |
|---|