| 0,0 → 1,970 |
|---|
| PROGRAM modify_anomaly |
| |
| c ****************************************************************************** |
| c * Read the modified and unmodified PV from the R file and * |
| c * perform some non-standard modifications: (a) Change of amplitude, * |
| c * (b) Stretching and shrinking along an axis, (c) rotation, or * |
| c * (d) change in position. |
| c * |
| c # modify_anomaly.sh R_20060116_18 shifty=-10,rot=10,stretchy=1.5,rot=-10,shifty=10 |
| c # modify_anomaly.sh R_20060116_18 fac=2 |
| c # modify_anomaly.sh R_20060116_18 cex=-30,cey=-30,rot=20 |
| c * * |
| c * Michael Sprenger / Spring, Summer 2007 * |
| c ****************************************************************************** |
| |
| implicit none |
| |
| c ----------------------------------------------------------------------------- |
| c Declaration of parameters and variables |
| c ----------------------------------------------------------------------------- |
| |
| c Input/output file and command string |
| character*80 pvsrcfile |
| character*80 commandstr |
| |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real mdv |
| |
| c 3d fields for calculation of qgPV and Ertel's PV |
| real,allocatable,dimension (:,:,:) :: pv1,pv2,ano |
| |
| c Numerical epsilon |
| real eps |
| parameter (eps=0.01) |
| |
| c Parameters for the transformations |
| real cex,cey ! Centre point for rotation |
| real angle ! Rotation angle |
| |
| c Auxiliary variables |
| real zpos |
| integer i,j,k |
| integer stat |
| character*80 varname |
| integer n |
| real par |
| character*80 com |
| |
| c ----------------------------------------------------------------------------- |
| c Preparations |
| c ----------------------------------------------------------------------------- |
| |
| print*,'********************************************************' |
| print*,'* MODIFY_ANOMALY *' |
| print*,'********************************************************' |
| |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) pvsrcfile |
| read(10,*) commandstr |
| close(10) |
| print* |
| print*,'Input file : ',trim(pvsrcfile) |
| print*,'Command : ',trim(commandstr) |
| 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* |
| |
| 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 modified and non-modified PV |
| allocate(pv1 (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv1' |
| allocate(pv2 (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv2' |
| allocate(ano (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating ano' |
| |
| c Read data from file |
| varname='PV' |
| call read_inp (pv1,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='PV_FILT' |
| call read_inp (pv2,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| |
| c Define the anomaly |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if ( (abs(pv1(i,j,k)-mdv).gt.eps).and. |
| > (abs(pv2(i,j,k)-mdv).gt.eps) ) then |
| ano(i,j,k)=pv1(i,j,k)-pv2(i,j,k) |
| else |
| ano(i,j,k)=0. |
| endif |
| enddo |
| enddo |
| enddo |
| |
| |
| c ------------------------------------------------------------------------------- |
| c Modifications |
| c ------------------------------------------------------------------------------- |
| |
| c Set the default values for parameters |
| cex = 0. |
| cey = 0. |
| angle = 0. |
| |
| c Set the counter for the command string |
| n=1 |
| |
| c Loop over all commands |
| 100 continue |
| |
| c Extract new command/parameter pair; exit if no new command |
| call next_command(commandstr,n,com,par) |
| if (com.eq.'nil') goto 200 |
| print*,trim(com),par |
| |
| c Multiply the anomaly by a constant factor |
| if (com.eq.'fac') then |
| call mod_factor (ano,par,nx,ny,nz,mdv) |
| endif |
| |
| c Set the centre point (needed for rotations) |
| if (com.eq.'cex') then |
| cex=par |
| endif |
| if (com.eq.'cey') then |
| cey=par |
| endif |
| |
| c Rotation around the centrre point |
| if (com.eq.'rot') then |
| call mod_rotation (ano,par,cex,cey, |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| endif |
| |
| c Shift in x or in y direction |
| if (com.eq.'shiftx') then |
| call mod_shift (ano,par,0., |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| endif |
| if (com.eq.'shifty') then |
| call mod_shift (ano,0.,par, |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| endif |
| |
| c Stretch/shrink in x or in y direction |
| if (com.eq.'stretchx') then |
| call mod_stretch (ano,par,1., |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| endif |
| if (com.eq.'stretchy') then |
| call mod_stretch (ano,1.,par, |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| endif |
| |
| |
| c Goto next command/parameter pair |
| goto 100 |
| |
| c All commands handled |
| 200 continue |
| |
| c -------------------------------------------------------------------------------- |
| c Write modified fields |
| c -------------------------------------------------------------------------------- |
| |
| c Build the modified PV distribution from the anomaly |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if ( (abs(pv1(i,j,k)-mdv).gt.eps).and. |
| > (abs(ano(i,j,k)-mdv).gt.eps) ) then |
| pv2(i,j,k)=pv1(i,j,k)-ano(i,j,k) |
| else |
| pv2(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| |
| c Write result to netcdf file |
| varname='PV_FILT' |
| call write_inp (pv2,varname,pvsrcfile,nx,ny,nz) |
| |
| c Write the modified anomaly also to the file |
| varname='PV_ANOM' |
| call write_inp (ano,varname,pvsrcfile,nx,ny,nz) |
| |
| end |
| |
| c ******************************************************************************** |
| c * Command handling and transformations * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Extract next command from command string |
| c -------------------------------------------------------------------------------- |
| |
| subroutine next_command (commandstr,n,com,par) |
| |
| c Given the command string <commandstr>, extract the next command/parameter pair |
| c <com,par> starting at position <n> of the command string. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| character*80 commandstr |
| character*80 com |
| real par |
| integer n |
| |
| c Auxiliary variables |
| integer i,j,k |
| |
| c Check whether end of command line reached |
| if (n.ge.80) then |
| com='nil' |
| par=0. |
| goto 120 |
| endif |
| |
| c Set indices to next next command/parameter pair |
| i=n |
| j=n |
| 100 if ( (commandstr(j:j).ne.'=').and.(j.lt.80) ) then |
| j=j+1 |
| goto 100 |
| endif |
| k=j+1 |
| 110 if ( (commandstr(k:k).ne.',' ).and. |
| > (commandstr(k:k).ne.';' ).and. |
| > (commandstr(k:k).ne.' ' ).and. |
| > (k .lt.80 ) ) then |
| k=k+1 |
| goto 110 |
| endif |
| |
| c Check whether this is a valid command/parameter pair |
| if ( ((j-1).lt.i).or.((k-1).lt.(j+1)) ) then |
| com='nil' |
| par=0. |
| goto 120 |
| endif |
| |
| c Extract tzhe command and the parameter |
| com=commandstr(i:j-1) |
| read(commandstr(j+1:k-1),*) par |
| |
| c Set the counter to the next command/parameter pair |
| n=k+1 |
| |
| c Exit point |
| 120 continue |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Multiply the anomaly by a factor |
| c -------------------------------------------------------------------------------- |
| |
| subroutine mod_factor (field,factor,nx,ny,nz,mdv) |
| |
| c Multiply the anomaly <field> by a constant factor <factor>. The grid and the |
| c missing data value are given by <nx,ny,nz,mdv>. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| real mdv |
| real factor |
| |
| c Parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer i,j,k |
| |
| c Do the transformation |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if ( (abs(field(i,j,k)-mdv).gt.eps) ) then |
| field(i,j,k)=factor*field(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Rotate the anomaly |
| c -------------------------------------------------------------------------------- |
| |
| subroutine mod_rotation (field,angle,cex,cey, |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| |
| c Rotate the anomaly <field> in the horizontal by the angle <angle>. The centre |
| c for the rotation is given by <cex,cey>, expressed in rotated longitude, latitude. |
| c The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>. |
| |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| real xmin,ymin,dx,dy |
| real mdv |
| real angle |
| real cex,cey |
| |
| c Parameters |
| real eps |
| parameter (eps=0.01) |
| real pi180 |
| parameter (pi180=3.14159/180.) |
| |
| c Auxiliary variables |
| integer i,j,k |
| real lon,lat,rlon,rlat |
| real xmax,ymax |
| real tmp1(0:nx,0:ny),tmp2(0:nx,0:ny) |
| |
| c Externals |
| real int2d |
| external int2d |
| |
| c Maximal grid extension |
| xmax=xmin+real(nx-1)*dx |
| ymax=ymin+real(ny-1)*dy |
| |
| c Do the Transformation |
| do k=0,nz |
| |
| c Copy level to 2d array |
| do i=0,nx |
| do j=0,ny |
| tmp1(i,j)=field(i,j,k) |
| enddo |
| enddo |
| |
| c Rotate each grid point |
| do i=0,nx |
| do j=0,ny |
| |
| c Lon/lat coordinates of grid point |
| lon=xmin+real(i)*dx-cex |
| lat=ymin+real(j)*dy-cey |
| |
| c Rotation |
| rlon = cex + lon*cos(angle*pi180) + lat*sin(angle*pi180) |
| rlat = cey - lon*sin(angle*pi180) + lat*cos(angle*pi180) |
| |
| c Do the interpolation |
| if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and. |
| > (rlat.gt.ymin).and.(rlat.lt.ymax) ) |
| > then |
| tmp2(i,j)=int2d(tmp1,rlon,rlat, |
| > dx,dy,xmin,ymin,nx,ny,mdv) |
| else |
| tmp2(i,j)=0. |
| endif |
| |
| enddo |
| enddo |
| |
| c Copy 2d array to level |
| do i=0,nx |
| do j=0,ny |
| field(i,j,k)=tmp2(i,j) |
| enddo |
| enddo |
| |
| enddo |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Shift the anomaly in x and y direction |
| c -------------------------------------------------------------------------------- |
| |
| subroutine mod_shift (field,shiftx,shifty, |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| |
| c Shift the anomaly <field> in the horizontal by the vector <shiftx,shifty>. |
| c The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>. |
| |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| real xmin,ymin,dx,dy |
| real mdv |
| real shiftx,shifty |
| |
| c Parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer i,j,k |
| real lon,lat,rlon,rlat |
| real xmax,ymax |
| real tmp1(0:nx,0:ny),tmp2(0:nx,0:ny) |
| |
| c Externals |
| real int2d |
| external int2d |
| |
| c Maximal grid extension |
| xmax=xmin+real(nx-1)*dx |
| ymax=ymin+real(ny-1)*dy |
| |
| c Do the Transformation |
| do k=0,nz |
| |
| c Copy level to 2d array |
| do i=0,nx |
| do j=0,ny |
| tmp1(i,j)=field(i,j,k) |
| enddo |
| enddo |
| |
| c Rotate each grid point |
| do i=0,nx |
| do j=0,ny |
| |
| c Lon/lat coordinates of grid point |
| lon=xmin+real(i)*dx |
| lat=ymin+real(j)*dy |
| |
| c shifted coordinates |
| rlon = lon - shiftx |
| rlat = lat - shifty |
| |
| c Do the interpolation |
| if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and. |
| > (rlat.gt.ymin).and.(rlat.lt.ymax) ) |
| > then |
| tmp2(i,j)=int2d(tmp1,rlon,rlat, |
| > dx,dy,xmin,ymin,nx,ny,mdv) |
| else |
| tmp2(i,j)=0. |
| endif |
| |
| enddo |
| enddo |
| |
| c Copy 2d array to level |
| do i=0,nx |
| do j=0,ny |
| field(i,j,k)=tmp2(i,j) |
| enddo |
| enddo |
| |
| enddo |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Stretch/shrink the anomaly in x and y direction |
| c -------------------------------------------------------------------------------- |
| |
| subroutine mod_stretch (field,stretchx,stretchy, |
| > nx,ny,nz,xmin,ymin,dx,dy,mdv) |
| |
| c Stretch the anomaly <field> in the horizontal by the fatcors <stretchx,stretchy>. |
| c The grid and the missing data value are given by <nx,ny,nz,xmin,ymin,dx,dy,mdv>. |
| |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| real xmin,ymin,dx,dy |
| real mdv |
| real stretchx,stretchy |
| |
| c Parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer i,j,k |
| real lon,lat,rlon,rlat |
| real xmax,ymax |
| real tmp1(0:nx,0:ny),tmp2(0:nx,0:ny) |
| |
| c Externals |
| real int2d |
| external int2d |
| |
| c Maximal grid extension |
| xmax=xmin+real(nx-1)*dx |
| ymax=ymin+real(ny-1)*dy |
| |
| c Do the Transformation |
| do k=0,nz |
| |
| c Copy level to 2d array |
| do i=0,nx |
| do j=0,ny |
| tmp1(i,j)=field(i,j,k) |
| enddo |
| enddo |
| |
| c Rotate each grid point |
| do i=0,nx |
| do j=0,ny |
| |
| c Lon/lat coordinates of grid point |
| lon=xmin+real(i)*dx |
| lat=ymin+real(j)*dy |
| |
| c shifted coordinates |
| rlon = 1./stretchx * lon |
| rlat = 1./stretchy * lat |
| |
| c Do the interpolation |
| if ( (rlon.gt.xmin).and.(rlon.lt.xmax).and. |
| > (rlat.gt.ymin).and.(rlat.lt.ymax) ) |
| > then |
| tmp2(i,j)=int2d(tmp1,rlon,rlat, |
| > dx,dy,xmin,ymin,nx,ny,mdv) |
| else |
| tmp2(i,j)=0. |
| endif |
| |
| enddo |
| enddo |
| |
| c Copy 2d array to level |
| do i=0,nx |
| do j=0,ny |
| field(i,j,k)=tmp2(i,j) |
| enddo |
| enddo |
| |
| enddo |
| |
| end |
| |
| c ------------------------------------------------------------------------ |
| c Two-dimensional interpolation |
| c ------------------------------------------------------------------------ |
| |
| real function int2d(ar,lon,lat, |
| > dx,dy,xmin,ymin,nx,ny,mdv) |
| |
| c Interpolate the field <ar(nx,ny)> to the position <lon,lat>. The |
| c grid is specified by <dx,dy,xmin,ymin,nx,ny>. |
| |
| implicit none |
| |
| c Declaration of subroutine paramters |
| integer nx,ny |
| real ar(0:nx,0:ny) |
| real lon,lat |
| real dx,dy,xmin,ymin |
| real mdv |
| |
| c Parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| real ri,rj |
| integer i,j,ir,ju |
| real frac0i,frac0j,frac1i,frac1j |
| real tmp,nor |
| |
| c Get index |
| ri=(lon-xmin)/dx |
| rj=(lat-ymin)/dy |
| i=int(ri) |
| j=int(rj) |
| if ((i.lt.0).or.(i.gt.nx).or. |
| > (j.lt.0).or.(j.gt.ny)) then |
| print*,'lat/lon interpolation not possible....' |
| stop |
| endif |
| |
| c Get inidices of left and upper neighbours |
| ir=i+1 |
| if (ir.gt.nx) ir=nx |
| ju=j+1 |
| if (ju.gt.ny) ju=ny |
| |
| c Get the weights for the bilinear interpolation |
| frac0i=ri-float(i) |
| frac0j=rj-float(j) |
| frac1i=1.-frac0i |
| frac1j=1.-frac0j |
| |
| c Bilinear interpolation with missing data check |
| tmp=0. |
| nor=0. |
| if ( (abs(ar(i ,j )-mdv).gt.eps) ) then |
| tmp = tmp + ar(i ,j ) * frac1i * frac1j |
| nor = nor + frac1i * frac1j |
| endif |
| if ( (abs(ar(i ,ju)-mdv).gt.eps) ) then |
| tmp = tmp + ar(i ,ju) * frac1i * frac0j |
| nor = nor + frac1i * frac0j |
| endif |
| if ( (abs(ar(ir,j )-mdv).gt.eps) ) then |
| tmp = tmp + ar(ir,j ) * frac0i * frac1j |
| nor = nor + frac0i * frac1j |
| endif |
| if ( (abs(ar(ir,ju)-mdv).gt.eps) ) then |
| tmp = tmp + ar(ir,ju) * frac0i * frac0j |
| nor = nor + frac0i * frac0j |
| endif |
| |
| c Return result |
| int2d=tmp/nor |
| |
| end |
| |
| c ******************************************************************************** |
| c * NETCDF INPUT AND OUTPUT * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Write input field to netcdf |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz) |
| |
| 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. |
| |
| 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 |
| |
| c Auxiliary variables |
| integer cdfid,stat |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| integer ndimin,outid,i,j,k |
| real tmp(0:nx,0:ny,0:nz) |
| integer ntimes |
| real time(200) |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| |
| c Get grid parameters from PV |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname='PV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| 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 |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Save variables (write definition, if necessary) |
| call cdfwopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname=fieldname |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time(1),0,field,stat) |
| print*,'W ',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 |
| |
| return |
| |
| c Exception handling |
| 998 print*,'Write_Inp: 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 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 from theta |
| 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 <THETA> 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='PV' |
| 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 |
| Property changes: |
|---|
| Added: svn:executable |
|---|