| /trunk/diag/calc_qgpv.f |
|---|
| 0,0 → 1,819 |
| PROGRAM main_calc_qgpv |
| c ******************************************************************************** |
| c * CALCULATE QUASI-GEOSTROPHIC PV ACCORDING TO ITS DEFINITION * |
| c * Michael Sprenger / Summer, Autumn 2006 * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Declaration of variables, parameters, externals and common blocks |
| c -------------------------------------------------------------------------------- |
| implicit none |
| c Input and output file |
| character*80 diagfile |
| character*80 referfile |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real deltax,deltay,deltaz |
| real mdv |
| real, allocatable,dimension (:,:) :: coriol |
| c Reference state |
| real, allocatable,dimension (:) :: nsqref |
| real, allocatable,dimension (:) :: thetaref |
| real, allocatable,dimension (:) :: rhoref |
| real, allocatable,dimension (:) :: pressref |
| real, allocatable,dimension (:) :: zref |
| c 3d fields for calculation of qgPV |
| real,allocatable,dimension (:,:,:) :: qgpv,th,uu,vv |
| c Auxiliary variables |
| integer i,j,k |
| integer stat |
| character*80 varname |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| print*,'********************************************************' |
| print*,'* CALC_QGPV *' |
| print*,'********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) diagfile |
| read(10,*) referfile |
| close(10) |
| print* |
| print*,trim(diagfile) |
| print*,trim(referfile) |
| print* |
| c Get lat/lon gid parameters from input file |
| call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > diagfile) |
| 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 Coriolis parameters |
| allocate(coriol(0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating coriol' |
| 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 3d fields |
| allocate(qgpv(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating qgpv' |
| allocate(uu(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating uu' |
| allocate(vv(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating vv' |
| allocate(th(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating th' |
| c -------------------------------------------------------------------------------- |
| c Calculate the qgPV from definition and put it onto file |
| c -------------------------------------------------------------------------------- |
| c Read data from file |
| varname='U' |
| call read_inp (uu,varname,diagfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='V' |
| call read_inp (vv,varname,diagfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='TH' |
| call read_inp (th,varname,diagfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| c Read reference profile and grid parameters |
| call read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol, |
| > referfile) |
| deltax=1000.*deltax |
| deltay=1000.*deltay |
| print*,'Deltax,deltay,deltaz =',deltax,deltay,deltaz |
| c Calculate qgPV |
| print*,'C qgPV' |
| call calc_qgpv (qgpv,uu,vv,th, |
| > rhoref,pressref,nsqref,thetaref,coriol, |
| > nx,ny,nz,deltax,deltay,deltaz,mdv) |
| c Write result to netcdf file |
| varname='QGPV_DIA' |
| call write_inp (qgpv,varname,diagfile,nx,ny,nz) |
| end |
| c ******************************************************************************** |
| c * NETCDF 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 max_th |
| 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 THETA |
| 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='TH' |
| 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 |
| 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 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 Read refernece profile from netcdf |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol, |
| > pvsrcfile) |
| c Read the reference profile from 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 nx,nny,nz : Grid dimension in x,y,z direction |
| c deltax,deltay,deltaz : Grid spacings used for calculations (m) |
| c coriol : Coriolis parameter (s^-1) |
| c pvsrcfile : Input 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) |
| real deltax,deltay,deltaz |
| real coriol (0:nx,0:ny) |
| 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 |
| real x(0:nx,0:ny),y(0:nx,0:ny) |
| real mean,count |
| 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 cdfopn(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 |
| c Read data from netcdf file |
| isok=0 |
| varname='NSQREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,nsqref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='RHOREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,rhoref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='THETAREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,thetaref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='PREREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,pressref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='ZREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,zref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='CORIOL' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,coriol,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='X' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,x,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='Y' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,y,stat) |
| if (stat.ne.0) goto 997 |
| c Close netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| c Determine the grid spacings <deltax, deltay, deltaz> |
| mean=0. |
| count=0. |
| do i=1,nx |
| do j=0,ny |
| mean=mean+abs(x(i,j)-x(i-1,j)) |
| count=count+1. |
| enddo |
| enddo |
| deltax=mean/count |
| mean=0. |
| count=0. |
| do j=1,ny |
| do i=0,nx |
| mean=mean+abs(y(i,j)-y(i,j-1)) |
| count=count+1. |
| enddo |
| enddo |
| deltay=mean/count |
| mean=0. |
| count=0. |
| do k=1,nz-1 |
| mean=mean+abs(zref(k+1)-zref(k-1)) |
| count=count+1. |
| enddo |
| deltaz=mean/count |
| return |
| c Exception handling |
| 997 print*,'Read_Ref: 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='TH' |
| 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 Calculate qgPV from wind and theta |
| c -------------------------------------------------------------------------------- |
| subroutine calc_qgpv (qgpv,uu,vv,th, |
| > rhoref,pressref,nsqref,thetaref,coriol, |
| > nx,ny,nz,deltax,deltay,deltaz,mdv) |
| c Calculate qgPV from wind and potential temperature according to |
| c equation 2.9 p 16 Thesis Rene Fehlmann. Note a cartesian grid with |
| c equidistant grid spacings <deltax,deltay,deltaz> is assumend. No |
| c 'correction' is made for spherical geoemtry. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real qgpv(0:nx,0:ny,0:nz) |
| real uu(0:nx,0:ny,0:nz) |
| real vv(0:nx,0:ny,0:nz) |
| real th(0:nx,0:ny,0:nz) |
| real rhoref(0:2*nz) |
| real nsqref(0:2*nz) |
| real thetaref(0:2*nz) |
| real pressref(0:2*nz) |
| real deltax,deltay,deltaz |
| real mdv |
| real coriol(0:nx,0:ny) |
| c Numerical epsilon and physical constants |
| real g |
| parameter (g=9.81) |
| real eps |
| parameter (eps=0.01) |
| real scale |
| parameter (scale=1e6) |
| c Auxiliary variables |
| integer i,j,k |
| integer kk,jj |
| real dvdx(0:nx,0:ny,0:nz) |
| real dudy(0:nx,0:ny,0:nz) |
| real dtdz(0:nx,0:ny,0:nz) |
| real t1,t2 |
| c Calculate horizontal derivatives dudy and dvdx of velocity |
| do k=0,nz |
| do j=0,ny |
| do i=0,nx |
| c Calculate dudy |
| if (j.eq.0) then |
| if ( (abs(uu(i,1,k)-mdv).gt.eps).and. |
| > (abs(uu(i,0,k)-mdv).gt.eps)) then |
| dudy(i,0,k)=(uu(i,1,k)-uu(i,0,k))/deltay |
| else |
| dudy(i,0,k)=mdv |
| endif |
| elseif (j.eq.ny) then |
| if ( (abs(uu(i,ny, k)-mdv).gt.eps).and. |
| > (abs(uu(i,ny-1,k)-mdv).gt.eps)) then |
| dudy(i,ny,k)=(uu(i,ny,k)-uu(i,ny-1,k))/deltay |
| else |
| dudy(i,ny,k)=mdv |
| endif |
| else |
| if ( (abs(uu(i,j+1,k)-mdv).gt.eps).and. |
| > (abs(uu(i,j-1,k)-mdv).gt.eps)) then |
| dudy(i,j,k)=(uu(i,j+1,k)-uu(i,j-1,k))/(2.*deltay) |
| else |
| dudy(i,j,k)=mdv |
| endif |
| endif |
| c Calculate dvdx |
| if (i.eq.0) then |
| if ((abs(vv(1,j,k)-mdv).gt.eps).and. |
| > (abs(vv(0,j,k)-mdv).gt.eps)) then |
| dvdx(0,j,k)=(vv(1,j,k)-vv(0,j,k))/deltax |
| else |
| dvdx(0,j,k)=mdv |
| endif |
| elseif (i.eq.nx) then |
| if ((abs(vv(nx, j,k)-mdv).gt.eps).and. |
| > (abs(vv(nx-1,j,k)-mdv).gt.eps)) then |
| dvdx(nx,j,k)=(vv(nx,j,k)-vv(nx-1,j,k))/deltax |
| else |
| dvdx(nx,j,k)=mdv |
| endif |
| else |
| if ((abs(vv(i+1,j,k)-mdv).gt.eps).and. |
| > (abs(vv(i-1,j,k)-mdv).gt.eps)) then |
| dvdx(i,j,k)=(vv(i+1,j,k)-vv(i-1,j,k))/(2.*deltax) |
| else |
| dvdx(i,j,k)=mdv |
| endif |
| endif |
| enddo |
| enddo |
| enddo |
| c Calculate vertical derivative of potential temperature |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if (k.eq.0) then |
| if ((abs(th(i,j,2)-mdv).gt.eps).and. |
| > (abs(th(i,j,1)-mdv).gt.eps)) then |
| t1=rhoref(2)*th(i,j,1)/thetaref(2)/nsqref(2)*g |
| t2=rhoref(0)*th(i,j,0)/thetaref(0)/nsqref(0)*g |
| dtdz(i,j,0)=(t1-t2)/deltaz |
| else |
| dtdz(i,j,0)=mdv |
| endif |
| else if (k.eq.nz) then |
| if ((abs(th(i,j,nz )-mdv).gt.eps).and. |
| > (abs(th(i,j,nz-1)-mdv).gt.eps)) then |
| kk=2*nz |
| t1=rhoref(kk )*th(i,j,nz )/ |
| > thetaref(kk)/nsqref(kk)*g |
| t2=rhoref(kk-2)*th(i,j,nz-1)/ |
| > thetaref(kk-2)/nsqref(kk-2)*9.8 |
| dtdz(i,j,nz)=(t1-t2)/deltaz |
| else |
| dtdz(i,j,nz)=mdv |
| endif |
| else |
| if ((abs(th(i,j,k+1)-mdv).gt.eps).and. |
| > (abs(th(i,j,k )-mdv).gt.eps).and. |
| > (abs(th(i,j,k-1)-mdv).gt.eps)) then |
| kk=2*k |
| t1=rhoref(kk+1)*(th(i,j,k)+th(i,j,k+1))/2. |
| > /thetaref(kk+1)/nsqref(kk+1)*g |
| t2=rhoref(kk-1)*(th(i,j,k)+th(i,j,k-1))/2. |
| > /thetaref(kk-1)/nsqref(kk-1)*g |
| dtdz(i,j,k)=(t1-t2)/deltaz |
| else |
| dtdz(i,j,k)=mdv |
| endif |
| endif |
| enddo |
| enddo |
| enddo |
| c Calculate qgPV |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| kk=2*k |
| if ((abs(dudy(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdz(i,j,k)-mdv).gt.eps)) then |
| qgpv(i,j,k)=-dudy(i,j,k)+dvdx(i,j,k)+ |
| > coriol(i,j)*dtdz(i,j,k)/rhoref(kk) |
| else |
| qgpv(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| Property changes: |
| Added: svn:executable |
| /trunk/diag/check_boundcon.f |
|---|
| 0,0 → 1,871 |
| PROGRAM set_boundcon |
| c ************************************************************************ |
| c * Set boundary conditions for inversion; lower and upper boundary * |
| c * conditions for potential temperature; lateral boundary conditions * |
| c * for zonal and meridional wind; in particular, missing data values * |
| c * are removed. * |
| c * * |
| c * Michael Sprenger / Summer 2006 * |
| c ************************************************************************ |
| c -------------------------------------------------------------------------------- |
| c Declaration of variables, parameters, externals and common blocks |
| c -------------------------------------------------------------------------------- |
| implicit none |
| c Input and output file |
| character*80 anomafile |
| character*80 referfile |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real mdv |
| real deltax,deltay,deltaz |
| real, allocatable,dimension (:,:) :: coriol |
| c Reference state |
| real, allocatable,dimension (:) :: nsqref |
| real, allocatable,dimension (:) :: thetaref |
| real, allocatable,dimension (:) :: rhoref |
| real, allocatable,dimension (:) :: pressref |
| real, allocatable,dimension (:) :: zref |
| C Boundary conditions |
| real, allocatable,dimension (:,:) :: thetatop |
| real, allocatable,dimension (:,:) :: thetabot |
| real, allocatable,dimension (:,:) :: ua |
| real, allocatable,dimension (:,:) :: ub |
| real, allocatable,dimension (:,:) :: va |
| real, allocatable,dimension (:,:) :: vb |
| c 3d arrays |
| real,allocatable,dimension (:,:,:) :: th_anom,pv_anom |
| real,allocatable,dimension (:,:,:) :: uu_anom,vv_anom |
| c Auxiliary variables |
| integer i,j,k,kk |
| integer stat |
| character*80 varname |
| integer n1,n2 |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| print*,'********************************************************' |
| print*,'* CHECK_BOUNDCON *' |
| print*,'********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) anomafile |
| read(10,*) referfile |
| close(10) |
| print* |
| print*,trim(anomafile) |
| print*,trim(referfile) |
| print* |
| c Get lat/lon gid parameters from input file |
| call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > anomafile) |
| 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 3d arrays |
| allocate(pv_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv_anom' |
| allocate(th_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating th_anom' |
| allocate(uu_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating uu_anom' |
| allocate(vv_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating vv_anom' |
| 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 Coriolis parameter |
| allocate(coriol(0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating f' |
| c Allocate memory for boundary conditions |
| allocate(thetatop(0:nx,0:ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array thetatop ***' |
| allocate(thetabot(0:nx,0:ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array thetabot ***' |
| allocate(ua(0:nx,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ua ***' |
| allocate(ub(0:nx,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ub ***' |
| allocate(va(0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array va ***' |
| allocate(vb(0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array vb ***' |
| c Read reference profile and ngrid parameters |
| call read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol, |
| > referfile) |
| deltax=1000.*deltax |
| deltay=1000.*deltay |
| print*,'Deltax,deltay,deltaz =',deltax,deltay,deltaz |
| c Read data from MOD file |
| varname='QGPV' |
| call read_inp (pv_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='TH' |
| call read_inp (th_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='U' |
| call read_inp (uu_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='V' |
| call read_inp (vv_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| c -------------------------------------------------------------------------------- |
| c Consistency check for boundary conditions and adaptions |
| c -------------------------------------------------------------------------------- |
| c Copy 3d to boundary conditions |
| do i=0,nx |
| do j=0,ny |
| thetatop(i,j)=th_anom(i,j,nz) |
| thetabot(i,j)=th_anom(i,j,0) |
| enddo |
| enddo |
| do i=0,nx |
| do k=0,nz |
| ua(i,k)=uu_anom(i, 0,k) |
| ub(i,k)=uu_anom(i,ny,k) |
| enddo |
| enddo |
| do j=0,ny |
| do k=0,nz |
| va(j,k)=vv_anom( 0,j,k) |
| vb(j,k)=vv_anom(nx,j,k) |
| enddo |
| enddo |
| c Check the lower and upper boundary condition for consistency check |
| print* |
| call combouncon(pv_anom,nsqref,rhoref,thetatop, |
| > thetabot,thetaref,coriol,ua,ub,va,vb, |
| > nx,ny,nz,deltax,deltay,deltaz) |
| print* |
| end |
| c ******************************************************************************** |
| c * NETCDF INPUT AND OUTPUT * |
| c ******************************************************************************** |
| 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 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='QGPV' |
| 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 Read refernece profile from netcdf |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol, |
| > pvsrcfile) |
| c Read the reference profile from 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 nx,nny,nz : Grid dimension in x,y,z direction |
| c deltax,deltay,deltaz : Grid spacings used for calculations (m) |
| c coriol : Coriolis parameter (s^-1) |
| c pvsrcfile : Input 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) |
| real deltax,deltay,deltaz |
| real coriol (0:nx,0:ny) |
| 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 |
| real x(0:nx,0:ny),y(0:nx,0:ny) |
| real mean,count |
| 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 cdfopn(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 |
| c Read data from netcdf file |
| isok=0 |
| varname='NSQREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,nsqref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='RHOREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,rhoref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='THETAREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,thetaref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='PREREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,pressref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='ZREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,zref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='CORIOL' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,coriol,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='X' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,x,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='Y' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,y,stat) |
| if (stat.ne.0) goto 997 |
| c Close netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| c Determine the grid spacings <deltax, deltay, deltaz> |
| mean=0. |
| count=0. |
| do i=1,nx |
| do j=0,ny |
| mean=mean+abs(x(i,j)-x(i-1,j)) |
| count=count+1. |
| enddo |
| enddo |
| deltax=mean/count |
| mean=0. |
| count=0. |
| do j=1,ny |
| do i=0,nx |
| mean=mean+abs(y(i,j)-y(i,j-1)) |
| count=count+1. |
| enddo |
| enddo |
| deltay=mean/count |
| mean=0. |
| count=0. |
| do k=1,nz-1 |
| mean=mean+abs(zref(k+1)-zref(k-1)) |
| count=count+1. |
| enddo |
| deltaz=mean/count |
| return |
| c Exception handling |
| 997 print*,'Read_Ref: Problem with input netcdf file... Stop' |
| stop |
| end |
| c ******************************************************************************** |
| c * BOUNDARY CONDITIONS - CONSISTENCY CHECK AND ADAPTIONS * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Boundary condition |
| c -------------------------------------------------------------------------------- |
| subroutine combouncon(pv,nsq,rhoref,thetatop, |
| > thetabot,thetaref,coriol, |
| > ua,ub,va,vb,nx,ny,nz,dx,dy,dz) |
| c Evaluate the consistency integrals A.7 from Rene's dissertation. This inegral |
| c is a necessary condition that the von Neumann problem has a unique solution. |
| c Adjust the upper and lower boundary conditions on <thetabot> and <thetatop>, so |
| c that the consitency check is ok. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real dx,dy,dz |
| real pv(0:nx,0:ny,0:nz) |
| real nsq(0:2*nz) |
| real rhoref(0:2*nz) |
| real thetatop(0:nx,0:ny) |
| real thetabot(0:nx,0:ny) |
| real thetaref(0:2*nz) |
| real coriol(0:nx,0:ny) |
| real ua(0:nx,0:nz) |
| real ub(0:nx,0:nz) |
| real va(0:ny,0:nz) |
| real vb(0:ny,0:nz) |
| c Numerical and physical parameters |
| real g |
| parameter (g=9.81) |
| c Auxiliary variables |
| integer i,j,k |
| real dxy,dxz,dyz,dxyz |
| real integr,denombot,denomtop,denom |
| real shifttop,shiftbot |
| c Set the area and volume infinitesimals for integration |
| dxy =dx*dy |
| dxz =dx*dz |
| dyz =dy*dz |
| dxyz=dx*dy*dz |
| c Init integration variables |
| integr=0. |
| c Inner |
| do k=1,nz-1 |
| do i=1,nx-1 |
| do j=1,ny-1 |
| integr=integr+dxyz*rhoref(2*k)*pv(i,j,k) |
| enddo |
| enddo |
| enddo |
| c ZY plane |
| do k=1,nz-1 |
| do j=1,ny-1 |
| integr=integr+dyz* |
| > rhoref(2*k)*(dx*pv(0, j,k)+va(j,k)) |
| c |
| integr=integr+dyz* |
| > rhoref(2*k)*(dx*pv(nx,j,k)-vb(j,k)) |
| enddo |
| enddo |
| c ZX plane |
| do k=1,nz-1 |
| do i=1,nx-1 |
| integr=integr+dxz* |
| > rhoref(2*k)*(dy*pv(i,0,k)-ua(i,k)) |
| c |
| integr=integr+dxz* |
| > rhoref(2*k)*(dy*pv(i,ny,k)+ub(i,k)) |
| enddo |
| enddo |
| c XY plane |
| do i=1,nx-1 |
| do j=1,ny-1 |
| integr=integr+dxy*rhoref(0)*( |
| > dz*pv(i,j,0)+coriol(i,j)*g*thetabot(i,j)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dxy*rhoref(2*nz)*( |
| > dz*pv(i,j,nz)-coriol(i,j)*g*thetatop(i,j)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| enddo |
| enddo |
| c X edges |
| do i=1,nx-1 |
| integr=integr+dx* |
| > rhoref(0)*(dyz*pv(i,0,0)- |
| > dz*ua(i,0)+dy*coriol(i,0)*g*thetabot(i,0)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dx* |
| > rhoref(0)*(dyz*pv(i,ny,0)+ |
| > dz*ub(i,0)+dy*coriol(i,ny)*g*thetabot(i,ny)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dx* |
| > rhoref(2*nz)*(dyz*pv(i,0,nz)- |
| > dz*ua(i,nz)-dy*coriol(i,0)*g*thetatop(i,0)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| integr=integr+dx* |
| > rhoref(2*nz)*(dyz*pv(i,ny,nz)+ |
| > dz*ub(i,nz)-dy*coriol(i,ny)*g*thetatop(i,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| enddo |
| c Y edges |
| do j=1,ny-1 |
| integr=integr+dy* |
| > rhoref(0)*(dxz*pv(0,j,0)+ |
| > dz*va(j,0)+dx*coriol(0,j)*g*thetabot(0,j)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dy* |
| > rhoref(0)*(dxz*pv(nx,j,0)- |
| > dz*vb(j,0)+dx*coriol(nx,j)*g*thetabot(nx,j)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dy* |
| > rhoref(2*nz)*(dxz*pv(0,j,nz)+ |
| > dz*va(j,nz)-dx*coriol(0,j)*g*thetatop(0,j)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+dy* |
| > rhoref(2*nz)*(dxz*pv(nx,j,nz)- |
| > dz*vb(j,nz)-dx*coriol(nx,j)*g*thetatop(nx,j)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| enddo |
| c Z edges |
| do k=1,nz-1 |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(0,0,k)+ |
| > dy*va(0,k)-dx*ua(0,k)) |
| c |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(nx,0,k)- |
| > dy*vb(0,k)-dx*ua(nx,k)) |
| c |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(0,ny,k)+ |
| > dy*va(ny,k)+dx*ub(0,k)) |
| c |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(nx,ny,k)- |
| > dy*vb(ny,k)+dx*ub(nx,k)) |
| enddo |
| c Points |
| integr=integr+rhoref(0)*(dxyz*pv(0,0,0)+ |
| > dyz*va(0,0)-dxz*ua(0,0)+ |
| > dxy*coriol(0,0)*g*thetabot(0,0)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(0)*(dxyz*pv(nx,0,0)- |
| > dyz*vb(0,0)-dxz*ua(nx,0)+ |
| > dxy*coriol(nx,0)*g*thetabot(nx,0)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(0)*(dxyz*pv(0,ny,0)+ |
| > dyz*va(ny,0)+dxz*ub(0,0)+ |
| > dxy*coriol(0,ny)*g*thetabot(0,ny)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(0)*(dxyz*pv(nx,ny,0)- |
| > dyz*vb(ny,0)+dxz*ub(nx,0)+ |
| > dxy*coriol(nx,ny)*g*thetabot(nx,ny)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(0,0,nz)+ |
| > dyz*va(0,nz)-dxz*ua(0,nz)- |
| > dxy*coriol(0,0)*g*thetatop(0,0)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(nx,0,nz)- |
| > dyz*vb(0,nz)-dxz*ua(nx,nz)- |
| > dxy*coriol(nx,0)*g*thetatop(nx,0)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(0,ny,nz)+ |
| > dyz*va(ny,nz)+dxz*ub(0,nz)- |
| > dxy*coriol(0,ny)*g*thetatop(0,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(nx,ny,nz)- |
| > dyz*vb(ny,nz)+dxz*ub(nx,nz)- |
| > dxy*coriol(nx,ny)*g*thetatop(nx,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| c Get the integrals from the reference state at bottom and top |
| denombot=0. |
| denomtop=0. |
| do i=0,nx |
| do j=0,ny |
| denombot=denombot+dxy* |
| > rhoref(0)*coriol(i,j)*g/ |
| > (nsq(0)*thetaref(0)) |
| c |
| denomtop=denomtop+dxy* |
| > rhoref(2*nz)*coriol(i,j)*g/ |
| > (nsq(2*nz)*thetaref(2*nz)) |
| enddo |
| enddo |
| denom=denomtop-denombot |
| c Determine the deviation of potential temperature from reference profile |
| shiftbot=0. |
| shifttop=0. |
| do i=0,nx |
| do j=0,ny |
| shifttop=shifttop+thetatop(i,j) |
| shiftbot=shiftbot+thetabot(i,j) |
| enddo |
| enddo |
| shifttop=shifttop/real((nx+1)*(ny+1)) |
| shiftbot=shiftbot/real((nx+1)*(ny+1)) |
| c Write some information about the consitency integrals |
| print*,'Consistency Check for boundary' |
| print*,' integ = ', integr |
| print*,' denombot = ', denombot |
| print*,' denomtop = ', denomtop |
| print*,' denom = ', denom |
| print*,' theta adjustment = ', integr/denom |
| print*,' theta shift @ top = ', shifttop, |
| > thetaref(2*nz) |
| print*,' theta shift @ bot = ', shiftbot, |
| > thetaref(0) |
| end |
| Property changes: |
| Added: svn:executable |
| /trunk/diag/check_boundcon.m |
|---|
| 0,0 → 1,267 |
| % ------------------------------------------------------------------------- |
| % Load files |
| % ------------------------------------------------------------------------- |
| % Base directory and filename |
| base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/'; |
| folder = base; |
| filename = 'Z1_20060115_18'; |
| disp([folder filename]) |
| % First image (otherwise first image is not correctly written) |
| figname = [base '/test.eps']; |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % Load variables from file (on model levels) |
| z_th = cdf_loadV(folder,filename,'TH'); |
| z_uu = cdf_loadV(folder,filename,'U'); |
| z_vv = cdf_loadV(folder,filename,'V'); |
| % ------------------------------------------------------------------------- |
| % Lower boundary condition for potential temperature |
| % ------------------------------------------------------------------------- |
| ilev=1; |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set projection |
| load coast |
| h=axesm('MapProjection','eqdcylin'); |
| setm(gca,'FLatLimit',[z_th.latmin z_th.latmax],'FLonLimit',[z_th.lonmin z_th.lonmax]); |
| h=plotm(lat,long); |
| gridm; |
| % Scale the plotting field for color map |
| fld=squeeze(z_th.var(ilev,:,:)); |
| c_map = scale_col(230:5:310,fld); |
| % Plot TH |
| lat=z_th.ymin + (0:z_th.ny-1) * z_th.dy; |
| lon=z_th.xmin + (0:z_th.nx-1) * z_th.dx; |
| [C,h]=contourfm(lat,lon,c_map.data,c_map.xtick); |
| for icnt = 1: length(h) |
| set( h(icnt), 'EdgeColor', 'none' ) |
| end |
| % Add color bar |
| colormap('default'); |
| ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,2); |
| caxis(c_map.caxis); |
| q=colorbar('vert'); |
| set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| title([ 'Lower Boundary Condition : Potential Temperature @' num2str(z_th.aklay(ilev)) ' m ASL' ]); |
| % Save figure |
| figname = [ base '/bound_lower_th_' filename '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % ------------------------------------------------------------------------- |
| % Upper boundary condition for potential temperature |
| % ------------------------------------------------------------------------- |
| ilev=z_th.nz; |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set projection |
| load coast |
| h=axesm('MapProjection','eqdcylin'); |
| setm(gca,'FLatLimit',[z_th.latmin z_th.latmax],'FLonLimit',[z_th.lonmin z_th.lonmax]); |
| h=plotm(lat,long); |
| gridm; |
| % Scale the plotting field for color map |
| fld=squeeze(z_th.var(ilev,:,:)); |
| c_map = scale_col(450:10:550,fld); |
| % Plot TH |
| lat=z_th.ymin + (0:z_th.ny-1) * z_th.dy; |
| lon=z_th.xmin + (0:z_th.nx-1) * z_th.dx; |
| [C,h]=contourfm(lat,lon,c_map.data,c_map.xtick); |
| for icnt = 1: length(h) |
| set( h(icnt), 'EdgeColor', 'none' ) |
| end |
| % Add color bar |
| colormap('default'); |
| ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,2); |
| caxis(c_map.caxis); |
| q=colorbar('vert'); |
| set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| title([ 'Upper Boundary Condition : Potential Temperature @ ' num2str(z_th.aklay(ilev)) ' m ASL' ]); |
| % Save figure |
| figname = [ base '/bound_upper_th_' filename '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % ------------------------------------------------------------------------- |
| % Southern lateral boundary condition for zonal wind |
| % ------------------------------------------------------------------------- |
| ilev=1; |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Scale the plotting field for color map |
| fld=squeeze(z_uu.var(:,ilev,:)); |
| c_map = scale_col(-50:10:50,fld); |
| % Plot |
| lev=z_uu.aklev/10000; |
| lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.dx; |
| [C,h]=contourf(lon,lev,c_map.data,c_map.xtick); |
| for icnt = 1: length(h) |
| set( h(icnt), 'EdgeColor', 'none' ) |
| end |
| % Add color bar |
| colormap('default'); |
| ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,2); |
| caxis(c_map.caxis); |
| q=colorbar('vert'); |
| set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| title([ 'Southern lateral Boundary Condition : Zonal Wind @ ' num2str(z_uu.ymin) '^\circ N' ]); |
| xlabel('Longitude'); |
| ylabel('Height [km]'); |
| % Save figure |
| figname = [ base '/bound_south_uu_' filename '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % ------------------------------------------------------------------------- |
| % Northern lateral boundary condition for zonal wind |
| % ------------------------------------------------------------------------- |
| ilev=z_uu.ny; |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Scale the plotting field for color map |
| fld=squeeze(z_uu.var(:,ilev,:)); |
| c_map = scale_col(-50:10:50,fld); |
| % Plot |
| lev=z_uu.aklev/10000; |
| lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.dx; |
| [C,h]=contourf(lon,lev,c_map.data,c_map.xtick); |
| for icnt = 1: length(h) |
| set( h(icnt), 'EdgeColor', 'none' ) |
| end |
| % Add color bar |
| colormap('default'); |
| ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,2); |
| caxis(c_map.caxis); |
| q=colorbar('vert'); |
| set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| title([ 'Northern lateral Boundary Condition : Zonal Wind @ ' num2str(z_uu.ymax) '^\circ N' ]); |
| xlabel('Longitude'); |
| ylabel('Height [km]'); |
| % Save figure |
| figname = [ base '/bound_north_uu_' filename '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % ------------------------------------------------------------------------- |
| % Western lateral boundary condition for meridional wind |
| % ------------------------------------------------------------------------- |
| ilev=1; |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Scale the plotting field for color map |
| fld=squeeze(z_vv.var(:,:,ilev)); |
| c_map = scale_col(-50:10:50,fld); |
| % Plot |
| lev=z_vv.aklev/10000; |
| lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy; |
| [C,h]=contourf(lat,lev,c_map.data,c_map.xtick); |
| for icnt = 1: length(h) |
| set( h(icnt), 'EdgeColor', 'none' ) |
| end |
| % Add color bar |
| colormap('default'); |
| ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,2); |
| caxis(c_map.caxis); |
| q=colorbar('vert'); |
| set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| title([ 'Western lateral Boundary Condition : Meridional Wind @ ' num2str(z_vv.xmin) '^\circ E' ]); |
| xlabel('Latitude'); |
| ylabel('Height [km]'); |
| % Save figure |
| figname = [ base '/bound_west_vv_' filename '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % ------------------------------------------------------------------------- |
| % Eastern lateral boundary condition for meridional wind |
| % ------------------------------------------------------------------------- |
| ilev=z_vv.nx; |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Scale the plotting field for color map |
| fld=squeeze(z_vv.var(:,:,ilev)); |
| c_map = scale_col(-50:10:50,fld); |
| % Plot |
| lev=z_vv.aklev/10000; |
| lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy; |
| [C,h]=contourf(lat,lev,c_map.data,c_map.xtick); |
| for icnt = 1: length(h) |
| set( h(icnt), 'EdgeColor', 'none' ) |
| end |
| % Add color bar |
| colormap('default'); |
| ctb(['/home/sprenger/ECFC_STE_Forecast/ctb_isen'],c_map.xtick,2); |
| caxis(c_map.caxis); |
| q=colorbar('vert'); |
| set(q,'ytick',c_map.xtick,'YTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| title([ 'Eastern lateral Boundary Condition : Meridional Wind @ ' num2str(z_vv.xmax) '^\circ E' ]); |
| xlabel('Latitude'); |
| ylabel('Height [km]'); |
| % Save figure |
| figname = [ base '/bound_east_vv_' filename '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| Property changes: |
| Added: svn:executable |
| /trunk/diag/difference.f |
|---|
| 0,0 → 1,241 |
| PROGRAM difference |
| c ************************************************************************** |
| c * Get the difference between two runs * |
| c * Michael Sprenger / Autumn 2006 * |
| c ************************************************************************** |
| implicit none |
| c -------------------------------------------------------------------------- |
| c Declaration of subroutine parameters |
| c -------------------------------------------------------------------------- |
| c Physical and numerical parameters |
| real eps |
| parameter (eps=0.01) |
| integer nmax |
| parameter (nmax=300*300*200) |
| c Variables for input file 1 |
| character*80 i1_filename |
| real i1_varmin(4),i1_varmax(4),i1_stag(4) |
| integer i1_vardim(4) |
| real i1_mdv |
| integer i1_ndim |
| integer i1_nx,i1_ny,i1_nz |
| real i1_time |
| integer i1_nvars |
| character*80 i1_vnam(100) |
| real i1_field(nmax) |
| c Variables for input file 2 |
| character*80 i2_filename |
| real i2_varmin(4),i2_varmax(4),i2_stag(4) |
| integer i2_vardim(4) |
| real i2_mdv |
| integer i2_ndim |
| integer i2_nx,i2_ny,i2_nz |
| real i2_time |
| integer i2_nvars |
| character*80 i2_vnam(100) |
| real i2_field(nmax) |
| c Variables for output file |
| character*80 o_filename |
| real o_varmin(4),o_varmax(4),o_stag(4) |
| integer o_vardim(4) |
| real o_mdv |
| integer o_ndim |
| integer o_nx,o_ny,o_nz |
| real o_time |
| real o_field(nmax) |
| c Auxiliary variables |
| integer i,j,k |
| integer cdfid,ierr |
| integer isok |
| character*80 cfn,varname |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| print*,'********************************************************' |
| print*,'* DIFFERENCE *' |
| print*,'********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) i1_filename |
| read(10,*) i2_filename |
| read(10,*) o_filename |
| close(10) |
| print* |
| print*,trim(i1_filename) |
| print*,trim(i2_filename) |
| print*,trim(o_filename) |
| print* |
| c Get a list of all variables on the two input files |
| call cdfopn(i1_filename,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,i1_nvars,i1_vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,i1_vnam(i1_nvars),i1_ndim,i1_mdv,i1_vardim, |
| > i1_varmin,i1_varmax,i1_stag,ierr) |
| call clscdf(cdfid,ierr) |
| print*,(trim(i1_vnam(i))//' ',i=1,i1_nvars) |
| call cdfopn(i2_filename,cdfid,ierr) |
| if (ierr.ne.0) goto 997 |
| call getvars(cdfid,i2_nvars,i2_vnam,ierr) |
| if (ierr.ne.0) goto 997 |
| call clscdf(cdfid,ierr) |
| print*,(trim(i2_vnam(i))//' ',i=1,i2_nvars) |
| print* |
| c Create the output file |
| o_ndim=i1_ndim |
| do i=1,i1_ndim |
| o_varmin(i)=i1_varmin(i) |
| o_varmax(i)=i1_varmax(i) |
| enddo |
| cfn=trim(o_filename)//'_cst' |
| call crecdf(trim(o_filename),cdfid,o_varmin,o_varmax, |
| > o_ndim,trim(cfn),ierr) |
| if (ierr.ne.0) goto 996 |
| call clscdf(cdfid,ierr) |
| c -------------------------------------------------------------------------------- |
| c Loop through all variables |
| c -------------------------------------------------------------------------------- |
| do i=1,i1_nvars |
| c Check wether the variable is available on both files |
| varname=i1_vnam(i) |
| if (varname.eq.'time') goto 100 |
| isok=0 |
| call check_varok (isok,varname,i2_vnam,i2_nvars) |
| if (isok.eq.0) goto 100 |
| c Read first input file |
| call cdfopn(i1_filename,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname,i1_ndim,i1_mdv,i1_vardim, |
| > i1_varmin,i1_varmax,i1_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,0.,0,i1_field,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| c Read second input file |
| call cdfopn(i2_filename,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname,i2_ndim,i2_mdv,i2_vardim, |
| > i2_varmin,i2_varmax,i2_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,0.,0,i2_field,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| c Consistency check |
| if (i1_ndim.ne.i2_ndim) then |
| print*,'Inconsistent input files... Stop',i1_ndim,i2_ndim |
| stop |
| endif |
| do j=1,3 |
| if ( (i1_vardim(j).ne.i2_vardim(j)).or. |
| > (abs(i1_varmin(j)-i2_varmin(j)).gt.eps).or. |
| > (abs(i1_varmax(j)-i2_varmax(j)).gt.eps).or. |
| > (abs(i1_stag(j)-i2_stag(j)).gt.eps)) then |
| print*,'Inconsistent input files... Stop' |
| print*,j,i1_varmin(j),i2_varmin(j) |
| print*,j,i1_varmax(j),i2_varmax(j) |
| print*,j,i1_vardim(j),i2_vardim(j) |
| print*,j,i1_stag(j), i2_stag(j) |
| stop |
| endif |
| enddo |
| c Get the difference |
| do j=1,i1_vardim(1)*i1_vardim(2)*i1_vardim(3) |
| if ( (abs(i1_field(j)-i1_mdv).gt.eps).and. |
| > (abs(i2_field(j)-i2_mdv).gt.eps) ) then |
| o_field(j)=i1_field(j)-i2_field(j) |
| else |
| o_field(j)=i1_mdv |
| endif |
| enddo |
| c Write to output file |
| o_ndim=i1_ndim |
| o_mdv =i1_mdv |
| do j=1,i1_ndim |
| o_vardim(j)=i1_vardim(j) |
| o_varmin(j)=i1_varmin(j) |
| o_varmax(j)=i1_varmax(j) |
| o_stag(j) =i1_stag(j) |
| enddo |
| call cdfwopn(o_filename,cdfid,ierr) |
| if (ierr.ne.0) goto 996 |
| call putdef(cdfid,varname,o_ndim,o_mdv,o_vardim, |
| > o_varmin,o_varmax,o_stag,ierr) |
| if (ierr.ne.0) goto 996 |
| call putdat(cdfid,varname,0.,0,o_field,ierr) |
| if (ierr.ne.0) goto 996 |
| print*,'W ',trim(varname),' ',trim(o_filename) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 996 |
| c Next |
| 100 continue |
| enddo |
| c ----------------------------------------------------------------- |
| c Exception handling and format specs |
| c ----------------------------------------------------------------- |
| stop |
| 998 print*,'Problems with input file 1 ',trim(i1_filename) |
| stop |
| 997 print*,'Problems with input file 2 ',trim(i2_filename) |
| stop |
| 996 print*,'Problems with output file ',trim(o_filename) |
| 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 |
| Property changes: |
| Added: svn:executable |
| /trunk/diag/hydrostatic.f |
|---|
| 0,0 → 1,542 |
| PROGRAM hydrostatic |
| c Calculate the geopotential and add it to the P file |
| c Michael Sprenger / Spring 2006 |
| implicit none |
| c --------------------------------------------------------------- |
| c Declaration of variables |
| c --------------------------------------------------------------- |
| c Variables for input P file : model level |
| character*80 ml_pfn |
| real ml_varmin(4),ml_varmax(4),ml_stag(4) |
| integer ml_vardim(4) |
| real ml_mdv |
| integer ml_ndim |
| integer ml_nx,ml_ny,ml_nz |
| real ml_xmin,ml_xmax,ml_ymin,ml_ymax,ml_dx,ml_dy |
| integer ml_ntimes |
| real ml_aklev(500),ml_bklev(500) |
| real ml_aklay(500),ml_bklay(500) |
| real ml_time |
| real ml_pollon,ml_pollat |
| integer ml_nvars |
| character*80 ml_vnam(100) |
| integer ml_idate(5) |
| real,allocatable, dimension (:,:) :: ml_ps,ml_zb |
| real,allocatable, dimension (:,:,:) :: ml_t3,ml_q3,ml_p3,ml_tv3 |
| real,allocatable, dimension (:,:,:) :: ml_zlay3 |
| c Variables for input Z file : pressure level |
| character*80 pl_zfn |
| real pl_varmin(4),pl_varmax(4),pl_stag(4) |
| integer pl_vardim(4) |
| real pl_mdv |
| integer pl_ndim |
| integer pl_nx,pl_ny,pl_nz |
| real pl_xmin,pl_xmax,pl_ymin,pl_ymax,pl_dx,pl_dy |
| integer pl_ntimes |
| real pl_aklev(500),pl_bklev(500) |
| real pl_aklay(500),pl_bklay(500) |
| real pl_time |
| real pl_pollon,pl_pollat |
| integer pl_nvars |
| character*80 pl_vnam(100) |
| integer pl_idate(5) |
| real,allocatable, dimension (:,:,:) :: pl_z3,pl_p3 |
| c Physical and numerical parameters |
| real g |
| parameter (g=9.80665) |
| real eps |
| parameter (eps=0.01) |
| real tzero |
| parameter (tzero=273.15) |
| real kappa |
| parameter (kappa=0.6078) |
| real zerodiv |
| parameter (zerodiv=0.0000000001) |
| real dpmin |
| parameter (dpmin=10.) |
| real rdg |
| parameter (rdg=29.271) |
| c Auxiliary variables |
| integer ierr |
| integer cdfid,cstid |
| character*80 cfn |
| integer stat |
| real time |
| real tv1(1000),z1(1000),p1(1000),f1(1000) |
| real spline_tv1(1000),spline_f1(1000),spline_z1(1000) |
| real pu,po,zu,zo,p,z,dp,p0,tvu,tvo,ff |
| integer i,j,k,l |
| integer lmin,n |
| character*80 varname,cdfname |
| integer isok |
| integer mode |
| c ----------------------------------------------------------------- |
| c Read input fields |
| c ----------------------------------------------------------------- |
| print*,'*********************************************************' |
| print*,'* hydrostatic *' |
| print*,'*********************************************************' |
| c Read in the parameter file |
| open(10,file='fort.10') |
| read(10,*) ml_pfn |
| read(10,*) pl_zfn |
| close(10) |
| c Decide which mode is used (1: reference from Z, 2: reference from ORO/PS) |
| if ( trim(ml_pfn).ne.trim(pl_zfn) ) then |
| mode=1 |
| print*,'Taking reference from Z ',trim(pl_zfn) |
| else |
| mode=2 |
| print*,'Taking reference from ORO/PS ',trim(pl_zfn) |
| endif |
| print*,trim(ml_pfn) |
| print*,trim(pl_zfn) |
| c Get grid description for P file : model level |
| call cdfopn(ml_pfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,ml_nvars,ml_vnam,ierr) |
| varname='T' |
| isok=0 |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.ne.1) goto 998 |
| call getdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| ml_nx =ml_vardim(1) |
| ml_ny =ml_vardim(2) |
| ml_nz =ml_vardim(3) |
| ml_xmin=ml_varmin(1) |
| ml_ymin=ml_varmin(2) |
| call getlevs(cstid,ml_nz,ml_aklev,ml_bklev,ml_aklay,ml_bklay,ierr) |
| call getgrid(cstid,ml_dx,ml_dy,ierr) |
| ml_xmax=ml_xmin+real(ml_nx-1)*ml_dx |
| ml_ymax=ml_ymin+real(ml_ny-1)*ml_dy |
| call gettimes(cdfid,ml_time,ml_ntimes,ierr) |
| call getstart(cstid,ml_idate,ierr) |
| call getpole(cstid,ml_pollon,ml_pollat,ierr) |
| call clscdf(cstid,ierr) |
| call clscdf(cdfid,ierr) |
| c Get grid description reference: either ORO(P) or Z(Z) |
| call cdfopn(pl_zfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,pl_nvars,pl_vnam,ierr) |
| if (mode.eq.1) then |
| varname='Z' |
| isok=0 |
| call check_varok (isok,varname,pl_vnam,pl_nvars) |
| if (isok.ne.1) goto 998 |
| call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim, |
| > pl_varmin,pl_varmax,pl_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| call getlevs(cstid,pl_nz,pl_aklev,pl_bklev, |
| > pl_aklay,pl_bklay,ierr) |
| call getgrid(cstid,pl_dx,pl_dy,ierr) |
| call gettimes(cdfid,pl_time,pl_ntimes,ierr) |
| call getstart(cstid,pl_idate,ierr) |
| call getpole(cstid,pl_pollon,pl_pollat,ierr) |
| else if (mode.eq.2) then |
| varname='ORO' |
| isok=0 |
| call check_varok (isok,varname,pl_vnam,pl_nvars) |
| if (isok.ne.1) goto 998 |
| call getdef(cdfid,varname,pl_ndim,pl_mdv,pl_vardim, |
| > pl_varmin,pl_varmax,pl_stag,ierr) |
| if (ierr.ne.0) goto 998 |
| call getgrid(cstid,pl_dx,pl_dy,ierr) |
| call gettimes(cdfid,pl_time,pl_ntimes,ierr) |
| call getstart(cstid,pl_idate,ierr) |
| call getpole(cstid,pl_pollon,pl_pollat,ierr) |
| endif |
| pl_nx =pl_vardim(1) |
| pl_ny =pl_vardim(2) |
| pl_nz =pl_vardim(3) |
| pl_xmin=pl_varmin(1) |
| pl_ymin=pl_varmin(2) |
| pl_xmax=pl_xmin+real(pl_nx-1)*pl_dx |
| pl_ymax=pl_ymin+real(pl_ny-1)*pl_dy |
| call clscdf(cstid,ierr) |
| call clscdf(cdfid,ierr) |
| c Consitency check for the grids |
| if ( (ml_nx.ne.pl_nx).or. |
| > (ml_ny.ne.pl_ny).or. |
| > (abs(ml_xmin-pl_xmin ).gt.eps).or. |
| > (abs(ml_ymin-pl_ymin ).gt.eps).or. |
| > (abs(ml_xmax-pl_xmax ).gt.eps).or. |
| > (abs(ml_ymax-pl_ymax ).gt.eps).or. |
| > (abs(ml_dx -pl_dx ).gt.eps).or. |
| > (abs(ml_dy -pl_dy ).gt.eps).or. |
| > (abs(ml_time-pl_time ).gt.eps).or. |
| > (abs(ml_pollon-pl_pollon).gt.eps).or. |
| > (abs(ml_pollat-pl_pollat).gt.eps)) then |
| print*,'Input P and Z grids are not consistent... Stop' |
| print*,'Xmin ',ml_xmin ,pl_xmin |
| print*,'Ymin ',ml_ymin ,pl_ymin |
| print*,'Xmax ',ml_xmax ,pl_xmax |
| print*,'Ymax ',ml_ymax ,pl_ymax |
| print*,'Dx ',ml_dx ,pl_dx |
| print*,'Dy ',ml_dy ,pl_dy |
| print*,'Pollon ',ml_pollon ,pl_pollon |
| print*,'Pollat ',ml_pollat ,pl_pollat |
| print*,'Time ',ml_time ,pl_time |
| stop |
| endif |
| c Allocate memory for all fields |
| allocate(ml_ps(ml_nx,ml_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_ps ***' |
| allocate(ml_zb(ml_nx,ml_ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_zb ***' |
| allocate(ml_p3(ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_p3 ***' |
| allocate(ml_t3(ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_t3 ***' |
| allocate(ml_q3(ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_q3 ***' |
| allocate(ml_tv3(ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_tv3 ***' |
| allocate(ml_zlay3(ml_nx,ml_ny,ml_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ml_zlay3 ***' |
| allocate(pl_z3(pl_nx,pl_ny,pl_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array pl_z3 ***' |
| allocate(pl_p3(pl_nx,pl_ny,pl_nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array pl_p3 ***' |
| c Read T, Q, PS from P file |
| call cdfopn(ml_pfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname='T' |
| call check_varok (isok,varname, ml_vnam,ml_nvars) |
| varname='Q' |
| call check_varok (isok,varname, ml_vnam,ml_nvars) |
| varname='PS' |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.ne.3) goto 998 |
| print*,'R T ',trim(ml_pfn) |
| call getdat(cdfid,'T',ml_time,0,ml_t3,ierr) |
| print*,'R Q ',trim(ml_pfn) |
| call getdat(cdfid,'Q',ml_time,0,ml_q3,ierr) |
| print*,'R PS ',trim(ml_pfn) |
| call getdat(cdfid,'PS',ml_time,1,ml_ps,ierr) |
| call clscdf(cdfid,ierr) |
| c Read ORO from P or Z from Z file |
| call cdfopn(pl_zfn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| if (mode.eq.1) then |
| isok=0 |
| varname='Z' |
| call check_varok (isok,varname,pl_vnam,pl_nvars) |
| if (isok.ne.1) goto 998 |
| print*,'R Z ',trim(pl_zfn) |
| call getdat(cdfid,varname,pl_time,0,pl_z3,ierr) |
| else if (mode.eq.2) then |
| isok=0 |
| varname='ORO' |
| call check_varok (isok,varname,pl_vnam,pl_nvars) |
| if (isok.ne.1) goto 998 |
| print*,'R ORO ',trim(pl_zfn) |
| call getdat(cdfid,varname,pl_time,1,pl_z3,ierr) |
| endif |
| call clscdf(cdfid,ierr) |
| c Set the values for the pressure on the pressure levels |
| do i=1,pl_nx |
| do j=1,pl_ny |
| do k=1,pl_nz |
| if (mode.eq.1) then |
| pl_p3(i,j,k)=pl_aklay(k) |
| else if (mode.eq.2) then |
| pl_p3(i,j,k)=ml_ps(i,j) |
| endif |
| enddo |
| enddo |
| enddo |
| c Calculate 3d pressure field |
| print*,'C P' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| ml_p3(i,j,k)=ml_aklay(k)+ml_bklay(k)*ml_ps(i,j) |
| enddo |
| enddo |
| enddo |
| c Calculate 3d virtual temperature |
| print*,'C TV' |
| do k=1,ml_nz |
| do i=1,ml_nx |
| do j=1,ml_ny |
| ml_tv3(i,j,k) = (ml_t3(i,j,k)+tzero)* |
| > (1.+kappa*ml_q3(i,j,k)) |
| enddo |
| enddo |
| enddo |
| c ----------------------------------------------------------------- |
| c Calculate geopotential |
| c ----------------------------------------------------------------- |
| c Integrate hydrostatic equation towards layers |
| print*,'C HYDROSTATIC EQUATION (LAYERS)' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| c Make the virtual temperature profile available |
| do k=1,ml_nz |
| p1 (ml_nz-k+1)=ml_p3 (i,j,k) |
| tv1(ml_nz-k+1)=ml_tv3(i,j,k) |
| enddo |
| call spline (p1,tv1,ml_nz,1.e30,1.e30,spline_tv1) |
| c Loop over all model levels |
| do k=1,ml_nz |
| c Get pressure at the grid point |
| p = ml_aklay(k)+ml_bklay(k)*ml_ps(i,j) |
| c Find nearest pressure level which is above topography |
| if (mode.eq.1) then |
| lmin=pl_nz |
| do l=1,pl_nz |
| if ((abs(p-pl_p3(i,j,l))).lt. |
| > (abs(p-pl_p3(i,j,lmin))) |
| > .and. |
| > (pl_p3(i,j,l).lt.ml_ps(i,j)) ) then |
| lmin=l |
| endif |
| enddo |
| else if (mode.eq.2) then |
| lmin=1 |
| endif |
| c Integrate hydrostatic equation from this level to the grid point |
| p0 = pl_p3(i,j,lmin) |
| n = nint(abs(p-p0)/dpmin) |
| if (n.lt.1) n=1 |
| dp = (p-p0)/real(n) |
| pu = p0 |
| z = pl_z3(i,j,lmin) |
| call splint(p1,tv1,spline_tv1,ml_nz,pu,tvu) |
| do l=1,n |
| po = pu+dp |
| call splint(p1,tv1,spline_tv1,ml_nz,po,tvo) |
| z = z + rdg*0.5*(tvu+tvo)*alog(pu/po) |
| tvu = tvo |
| pu = po |
| enddo |
| c Set the geopotential at the grid point |
| ml_zlay3(i,j,k) = z |
| enddo |
| enddo |
| enddo |
| c ----------------------------------------------------------------- |
| c Calculate height of topography |
| c ----------------------------------------------------------------- |
| if (mode.eq.1) then |
| print*,'C TOPOGRAPHY' |
| do i=1,ml_nx |
| do j=1,ml_ny |
| c Make the z/p profile available |
| do k=1,ml_nz |
| p1(ml_nz-k+1)=ml_p3(i,j,k) |
| z1(ml_nz-k+1)=ml_zlay3(i,j,k) |
| enddo |
| c Cubic spline interpolation |
| call spline (p1,z1,ml_nz,1.e30,1.e30,spline_z1) |
| call splint (p1,z1,spline_z1,ml_nz,ml_ps(i,j),ml_zb(i,j)) |
| enddo |
| enddo |
| endif |
| c ----------------------------------------------------------------- |
| c Write the topography and geopotential to P file |
| c ----------------------------------------------------------------- |
| c Open P file |
| call cdfwopn(trim(ml_pfn),cdfid,ierr) |
| c Write orography (levels) |
| if (mode.eq.1) then |
| varname='ORO' |
| print*,'W ',trim(varname),' ',trim(ml_pfn) |
| isok=0 |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| ml_vardim(3)=1 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| ml_vardim(3)=ml_nz |
| if (ierr.ne.0) goto 997 |
| endif |
| call putdat(cdfid,varname,ml_time,1,ml_zb,ierr) |
| if (ierr.ne.0) goto 997 |
| endif |
| c Write geopotential on layers |
| varname='Z_DIA' |
| print*,'W ',trim(varname),' ',trim(ml_pfn) |
| isok=0 |
| call check_varok (isok,varname,ml_vnam,ml_nvars) |
| if (isok.eq.0) then |
| ml_stag(3)=-0.5 |
| call putdef(cdfid,varname,ml_ndim,ml_mdv,ml_vardim, |
| > ml_varmin,ml_varmax,ml_stag,ierr) |
| if (ierr.ne.0) goto 997 |
| endif |
| call putdat(cdfid,varname,ml_time,0,ml_zlay3,ierr) |
| if (ierr.ne.0) goto 997 |
| c Close P file |
| call clscdf(cdfid,ierr) |
| c ----------------------------------------------------------------- |
| c Exception handling and format specs |
| c ----------------------------------------------------------------- |
| stop |
| 998 print*,'Z: Problems with input from m level' |
| stop |
| 997 print*,'Z: Problems with output on m level' |
| stop |
| 996 print*,'F: Problems with input from m level' |
| stop |
| 995 print*,'F: Problems with output on z level' |
| stop |
| end |
| c ------------------------------------------------------------- |
| c Natural cubic spline |
| c ------------------------------------------------------------- |
| SUBROUTINE spline(x,y,n,yp1,ypn,y2) |
| INTEGER n,NMAX |
| REAL yp1,ypn,x(n),y(n),y2(n) |
| PARAMETER (NMAX=500) |
| INTEGER i,k |
| REAL p,qn,sig,un,u(NMAX) |
| if (yp1.gt..99e30) then |
| y2(1)=0. |
| u(1)=0. |
| else |
| y2(1)=-0.5 |
| u(1)=(3./(x(2)-x(1)))*((y(2)-y(1))/(x(2)-x(1))-yp1) |
| endif |
| do 11 i=2,n-1 |
| sig=(x(i)-x(i-1))/(x(i+1)-x(i-1)) |
| p=sig*y2(i-1)+2. |
| y2(i)=(sig-1.)/p |
| u(i)=(6.*((y(i+1)-y(i))/(x(i+ |
| *1)-x(i))-(y(i)-y(i-1))/(x(i)-x(i-1)))/(x(i+1)-x(i-1))-sig* |
| *u(i-1))/p |
| 11 continue |
| if (ypn.gt..99e30) then |
| qn=0. |
| un=0. |
| else |
| qn=0.5 |
| un=(3./(x(n)-x(n-1)))*(ypn-(y(n)-y(n-1))/(x(n)-x(n-1))) |
| endif |
| y2(n)=(un-qn*u(n-1))/(qn*y2(n-1)+1.) |
| do 12 k=n-1,1,-1 |
| y2(k)=y2(k)*y2(k+1)+u(k) |
| 12 continue |
| return |
| END |
| SUBROUTINE splint(xa,ya,y2a,n,x,y) |
| INTEGER n |
| REAL x,y,xa(n),y2a(n),ya(n) |
| INTEGER k,khi,klo |
| REAL a,b,h |
| klo=1 |
| khi=n |
| 1 if (khi-klo.gt.1) then |
| k=(khi+klo)/2 |
| if(xa(k).gt.x)then |
| khi=k |
| else |
| klo=k |
| endif |
| goto 1 |
| endif |
| h=xa(khi)-xa(klo) |
| if (h.eq.0.) pause 'bad xa input in splint' |
| a=(xa(khi)-x)/h |
| b=(x-xa(klo))/h |
| y=a*ya(klo)+b*ya(khi)+((a**3-a)*y2a(klo)+(b**3-b)*y2a(khi))*(h** |
| *2)/6. |
| return |
| 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)>. |
| c If this is the case, <isok> is incremented by 1. Otherwise <isok> |
| c 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 |
| Property changes: |
| Added: svn:executable |
| /trunk/diag/qvec_analysis.f |
|---|
| 0,0 → 1,1468 |
| PROGRAM z2s |
| c Calculate secondary fields on z levels |
| c Michael Sprenger / Summer 2006 |
| implicit none |
| c --------------------------------------------------------------- |
| c Declaration of parameters |
| c --------------------------------------------------------------- |
| real tzero |
| parameter (tzero=273.15) |
| real kappa |
| parameter (kappa=0.6078) |
| real rd |
| parameter (rd=287.) |
| c --------------------------------------------------------------- |
| c Declaration of variables |
| c --------------------------------------------------------------- |
| c Variables for output Z file : height level |
| character*80 cfn |
| real varmin(4),varmax(4),stag(4) |
| integer vardim(4) |
| real mdv |
| integer ndim |
| integer nx,ny,nz |
| real xmin,xmax,ymin,ymax,dx,dy |
| integer ntimes |
| real aklev(1000),bklev(1000) |
| real aklay(1000),bklay(1000) |
| real time |
| real pollon,pollat |
| integer idate(5) |
| integer nfields |
| character*80 field_nam(100) |
| real,allocatable, dimension (:,:,:,:) :: field_dat |
| real,allocatable, dimension (:,:,:) :: z3 |
| real,allocatable, dimension (:,:) :: x2,y2,f2,oro |
| real,allocatable, dimension (:,:,:) :: out1,out2 |
| real,allocatable, dimension (:,:,:) :: inp |
| real,allocatable,dimension (:) :: nsqref |
| real,allocatable,dimension (:) :: thetaref |
| real,allocatable,dimension (:) :: tref |
| real,allocatable,dimension (:) :: rhoref |
| real,allocatable,dimension (:) :: pressref |
| real,allocatable,dimension (:) :: zref |
| real deltax,deltay,deltaz |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| integer cdfid,cstid |
| character*3 mode |
| c Parameter file |
| character*80 fieldname |
| integer nfilt |
| character*80 ofn,gri |
| c Auxiliary variables |
| integer ierr,stat |
| integer i,j,k,n |
| real,allocatable, dimension (:,:) :: tmp |
| character*80 vnam2(100) |
| integer nvars2 |
| c --------------------------------------------------------------- |
| c Preparations |
| c --------------------------------------------------------------- |
| print*,'*********************************************************' |
| print*,'* qvec_analysis *' |
| print*,'*********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) fieldname |
| read(10,*) ofn |
| read(10,*) gri |
| read(10,*) nfilt |
| close(10) |
| c Decide whether to enter ANO or MOD/ORG mode |
| mode = ofn(1:3) |
| if ( (mode.ne.'ANO').and. |
| > (mode.ne.'MOD').and. |
| > (mode.ne.'ORG') ) |
| >then |
| print*,'Unsupported mode ',trim(mode) |
| stop |
| endif |
| c Get grid description for Z file : height level |
| call cdfopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getcfn(cdfid,cfn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(cfn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,'T',ndim,mdv,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) |
| 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=xmin+real(nx-1)*dx |
| ymax=ymin+real(ny-1)*dy |
| call gettimes(cdfid,time,ntimes,ierr) |
| if (ierr.ne.0) goto 998 |
| call getstart(cstid,idate,ierr) |
| if (ierr.ne.0) goto 998 |
| call getpole(cstid,pollon,pollat,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c Get a list of all variables on the GRID file |
| if (trim(gri).ne.trim(ofn)) then |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars2,vnam2,ierr) |
| if (ierr.ne.0) goto 998 |
| do i=1,nvars2 |
| vnam(nvars+i)=vnam2(i) |
| enddo |
| nvars=nvars+nvars2 |
| endif |
| c Check whether calculation of <fieldname> is supported |
| if ( (fieldname.ne.'GEO' ).and. |
| > (fieldname.ne.'AGEO' ).and. |
| > (fieldname.ne.'DIV_UV' ).and. |
| > (fieldname.ne.'QVEC' ).and. |
| > (fieldname.ne.'DIV_Q' ) ) then |
| print*,'Variable not supported ',trim(fieldname) |
| stop |
| endif |
| c Set dependencies |
| if (fieldname.eq.'GEO') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| elseif (fieldname.eq.'AGEO') then |
| nfields=4 |
| field_nam(1)='U' |
| field_nam(2)='UG' |
| field_nam(3)='V' |
| field_nam(4)='VG' |
| else if (fieldname.eq.'DIV_UV') then |
| nfields=2 |
| field_nam(1)='U' |
| field_nam(2)='V' |
| else if (fieldname.eq.'QVEC') then |
| nfields=3 |
| field_nam(1)='UG' |
| field_nam(2)='VG' |
| field_nam(3)='TH' |
| else if (fieldname.eq.'DIV_Q') then |
| nfields=2 |
| field_nam(1)='QX' |
| field_nam(2)='QY' |
| endif |
| c Allocate memory |
| allocate(field_dat(nfields,nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating field_dat' |
| allocate(out1(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating out1' |
| allocate(out2(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating out2' |
| allocate(inp(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating inp' |
| allocate(z3(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating z3' |
| allocate(x2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating x2' |
| allocate(y2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating y2' |
| allocate(f2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating f2' |
| allocate(oro(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating oro' |
| 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' |
| allocate(tref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating tref' |
| c Allocate auxiliary fields |
| allocate(tmp(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating tmp' |
| c Read reference profile and grid parameters |
| call read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,f2,oro,gri) |
| c Read X grid |
| varname='X' |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| print*,'Variable ',trim(varname),' missing... Stop' |
| stop |
| endif |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,x2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(gri) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c Read Y grid |
| varname='Y' |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| print*,'Variable ',trim(varname),' missing... Stop' |
| stop |
| endif |
| call cdfopn(gri,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,y2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(gri) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c Calculate reference profile of temperature |
| print*,'C T_ref = TH_ref * ( P_ref/1000)^kappa' |
| do i=0,2*nz |
| tref(i) = thetaref(i) * ( pressref(i)/1000. ) ** kappa |
| enddo |
| c Init height levels |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| z3(i,j,k)=zref(2*k-1) |
| enddo |
| enddo |
| enddo |
| c Load needed variables |
| do n=1,nfields |
| c Check whether variable is available on file |
| varname=field_nam(n) |
| isok=0 |
| call check_varok (isok,varname,vnam,nvars) |
| c Variable is available on file |
| if (isok.eq.1) then |
| call cdfopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdat(cdfid,varname,time,0,inp,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(ofn) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k)=inp(i,j,k) |
| enddo |
| enddo |
| enddo |
| else |
| print*,'Variable missing : ',trim(varname) |
| stop |
| endif |
| enddo |
| c Add reference profile if ANO file is provided |
| if ( mode.eq.'ANO' ) then |
| print*,'C T -> T_ano + T_ref' |
| n=0 |
| do i=1,nfields |
| if (field_nam(i).eq.'T') n=i |
| enddo |
| if (n.ne.0) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k) = field_dat(n,i,j,k) + tref(2*k-1) |
| enddo |
| enddo |
| enddo |
| endif |
| print*,'C TH -> TH_ano + TH_ref' |
| n=0 |
| do i=1,nfields |
| if (field_nam(i).eq.'TH') n=i |
| enddo |
| if (n.ne.0) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k) = field_dat(n,i,j,k)+thetaref(2*k-1) |
| enddo |
| enddo |
| enddo |
| endif |
| print*,'C P -> P_ano + P_ref' |
| n=0 |
| do i=1,nfields |
| if (field_nam(i).eq.'P') n=i |
| enddo |
| if (n.ne.0) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k) = field_dat(n,i,j,k)+pressref(2*k-1) |
| enddo |
| enddo |
| enddo |
| endif |
| endif |
| c Change units: P (hPa -> Pa), T(K -> degC) |
| n=0 |
| do i=1,nfields |
| if (field_nam(i).eq.'P') n=i |
| enddo |
| if (n.ne.0) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| field_dat(n,i,j,k)=100.*field_dat(n,i,j,k) |
| enddo |
| enddo |
| enddo |
| endif |
| n=0 |
| do i=1,nfields |
| if (field_nam(i).eq.'T') n=i |
| enddo |
| if (n.ne.0) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ( field_dat(n,i,j,1).gt.100. ) then |
| field_dat(n,i,j,k)=field_dat(n,i,j,k) - tzero |
| endif |
| enddo |
| enddo |
| enddo |
| endif |
| c ---------------------------------------------------------------- |
| c Calculations |
| c ---------------------------------------------------------------- |
| c Geostrophic wind (GEO) |
| if (fieldname.eq.'GEO') then |
| print*,'C ',trim(fieldname) |
| field_nam(1)='RHO' |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ( mode.eq.'ANO' ) then |
| field_dat(1,i,j,k)=rhoref(2*k-1) |
| else |
| field_dat(1,i,j,k)= |
| > 1./rd*field_dat(2,i,j,k)/(field_dat(1,i,j,k)+tzero) |
| endif |
| enddo |
| enddo |
| enddo |
| call calc_geo (out1,out2, ! (Ug,Vg) |
| > field_dat(1,:,:,:), ! RHO |
| > field_dat(2,:,:,:), ! P |
| > z3,x2,y2,f2, ! Z,X,Y,CORIOL |
| > nx,ny,nz,mdv) |
| c Ageostrophic wind (AGEO) |
| elseif (fieldname.eq.'AGEO') then |
| print*,'C ',trim(fieldname) |
| call calc_ageo (out1,out2, ! (Ua,Va) |
| > field_dat(1,:,:,:), ! U |
| > field_dat(2,:,:,:), ! UG |
| > field_dat(3,:,:,:), ! V |
| > field_dat(4,:,:,:), ! VG |
| > nx,ny,nz,mdv) |
| c Divergence of wind (DIV_UV) |
| else if (fieldname.eq.'DIV_UV') then |
| print*,'C ',trim(fieldname) |
| call calc_div_uv (out1, ! Div(U,V)) |
| > field_dat(1,:,:,:), ! U |
| > field_dat(2,:,:,:), ! V |
| > z3,x2,y2,f2, ! Z,X,Y,CORIOL |
| > nx,ny,nz,mdv) |
| c Q vector (QVEC) |
| else if (fieldname.eq.'QVEC') then |
| print*,'C ',trim(fieldname) |
| call calc_qvec (out1,out2, ! (Qx,Qy) |
| > field_dat(1,:,:,:), ! UG |
| > field_dat(2,:,:,:), ! VG |
| > field_dat(3,:,:,:), ! TH |
| > z3,x2,y2,f2, ! Z,X,Y,CORIOL |
| > nx,ny,nz,mdv) |
| c Divergence of Q vector (DIV_Q) |
| else if (fieldname.eq.'DIV_Q') then |
| print*,'C ',trim(fieldname) |
| call calc_div_q (out1, ! Div(U,V)) |
| > field_dat(1,:,:,:), ! QX |
| > field_dat(2,:,:,:), ! QY |
| > z3,x2,y2,f2, ! Z,X,Y,CORIOL |
| > nx,ny,nz,mdv) |
| endif |
| c Horizontal filtering of the resulting fields |
| print*,'F ',trim(fieldname) |
| do k=1,nz |
| do i=1,nx |
| do j=1,ny |
| tmp(i,j)=out1(i,j,k) |
| enddo |
| enddo |
| do n=1,nfilt |
| call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0) |
| enddo |
| do i=1,nx |
| do j=1,ny |
| out1(i,j,k)=tmp(i,j) |
| enddo |
| enddo |
| do i=1,nx |
| do j=1,ny |
| tmp(i,j)=out2(i,j,k) |
| enddo |
| enddo |
| do n=1,nfilt |
| call filt2d (tmp,tmp,nx,ny,1.,mdv,0,0,0,0) |
| enddo |
| do i=1,nx |
| do j=1,ny |
| out2(i,j,k)=tmp(i,j) |
| enddo |
| enddo |
| enddo |
| c ---------------------------------------------------------------- |
| c Save result onto netcdf file |
| c ---------------------------------------------------------------- |
| c Open output file |
| call cdfwopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c Save geostrophic wind |
| if (fieldname.eq.'GEO') then |
| isok=0 |
| varname='UG' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out1,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| isok=0 |
| varname='VG' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| c Save ageostrophic wind |
| elseif (fieldname.eq.'AGEO') then |
| isok=0 |
| varname='UA' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out1,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| isok=0 |
| varname='VA' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| c Save divergence of wind field |
| else if (fieldname.eq.'DIV_UV') then |
| isok=0 |
| varname='DIV_UV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out1,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| c Save components of Q vector |
| else if (fieldname.eq.'QVEC') then |
| isok=0 |
| varname='QX' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out1,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| isok=0 |
| varname='QY' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| c Save divergence of wind field |
| else if (fieldname.eq.'DIV_Q') then |
| isok=0 |
| varname='DIV_Q' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndim,mdv,vardim, |
| > varmin,varmax,stag,ierr) |
| if (ierr.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time,0,out1,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| endif |
| c Close output file |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c ---------------------------------------------------------------- |
| c Exception handling |
| c ---------------------------------------------------------------- |
| stop |
| 998 print*,'Problem with netcdf file' |
| stop |
| end |
| c **************************************************************** |
| c * SUBROUTINE SECTION: AUXILIARY ROUTINES * |
| c **************************************************************** |
| 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)>. |
| c If this is the case, <isok> is incremented by 1. Otherwise <isok> |
| c 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 * SUBROUTINE SECTION: CALCULATE SECONDARY FIELDS * |
| c **************************************************************** |
| c ---------------------------------------------------------------- |
| c Calculate geostrophic wind components |
| c ---------------------------------------------------------------- |
| subroutine calc_geo (ug,vg,rho,p, |
| > z3,x2,y2,f2,nx,ny,nz,mdv) |
| c Calculate the geostrophic wind components (ug,vg) if the temperature |
| c (t) and the pressure (p) are given. The grid and the missing data |
| c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real ug (nx,ny,nz) |
| real vg (nx,ny,nz) |
| real rho(nx,ny,nz) |
| real p (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real f2 (nx,ny) |
| real mdv |
| c Physical parameters and numerical constants |
| real eps |
| parameter (eps=0.01) |
| real g |
| parameter (g=9.80616) |
| c Auxiliray variables |
| integer i,j,k |
| real dpdx(nx,ny,nz) |
| real dpdy(nx,ny,nz) |
| c Calculate horizontal derivatives of pressure |
| call deriv(dpdx,p,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv(dpdy,p,'y',z3,x2,y2,nx,ny,nz,mdv) |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(rho (i,j,k)-mdv).gt.eps).and. |
| > (abs(dpdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dpdy(i,j,k)-mdv).gt.eps)) then |
| ug(i,j,k)=-1./(rho(i,j,k)*f2(i,j))*dpdy(i,j,k) |
| vg(i,j,k)= 1./(rho(i,j,k)*f2(i,j))*dpdx(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ---------------------------------------------------------------- |
| c Calculate ageostrophic wind components |
| c ---------------------------------------------------------------- |
| subroutine calc_ageo (ua,va,u,ug,v,vg,nx,ny,nz,mdv) |
| c Calculate the geostrophic wind components (ug,vg) if the temperature |
| c (t) and the pressure (p) are given. The grid and the missing data |
| c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real ua (nx,ny,nz) |
| real va (nx,ny,nz) |
| real u (nx,ny,nz) |
| real ug (nx,ny,nz) |
| real v (nx,ny,nz) |
| real vg (nx,ny,nz) |
| real mdv |
| c Parameters |
| real eps |
| parameter (eps=0.01) |
| c Auxiliray variables |
| integer i,j,k |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(u (i,j,k)-mdv).gt.eps).and. |
| > (abs(ug(i,j,k)-mdv).gt.eps).and. |
| > (abs(v (i,j,k)-mdv).gt.eps).and. |
| > (abs(vg(i,j,k)-mdv).gt.eps)) then |
| ua(i,j,k)=u(i,j,k) - ug(i,j,k) |
| va(i,j,k)=v(i,j,k) - vg(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ---------------------------------------------------------------- |
| c Calculate divergence of wind field |
| c ---------------------------------------------------------------- |
| subroutine calc_div_uv (div_uv,u,v, |
| > z3,x2,y2,f2,nx,ny,nz,mdv) |
| c Calculate the divergence (div_uv) of the horizontal wind field if+ |
| c the wind components (u,v) are given. The grid and the missing data |
| c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real div_uv(nx,ny,nz) |
| real u (nx,ny,nz) |
| real v (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real f2 (nx,ny) |
| real mdv |
| c Physical parameters and numerical constants |
| real eps |
| parameter (eps=0.01) |
| c Auxiliray variables |
| integer i,j,k |
| real rho |
| real dudx(nx,ny,nz) |
| real dvdy(nx,ny,nz) |
| c Calculate horizontal derivatives of pressure |
| call deriv(dudx,u,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv(dvdy,v,'y',z3,x2,y2,nx,ny,nz,mdv) |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(dudx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdy(i,j,k)-mdv).gt.eps)) then |
| div_uv(i,j,k)=dudx(i,j,k)+dvdy(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ---------------------------------------------------------------- |
| c Calculate Q vector |
| c ---------------------------------------------------------------- |
| subroutine calc_qvec (qx3,qy3, |
| > th3,u3,v3, |
| > z3,x2,y2,f2,nx,ny,nz,mdv) |
| c Calculate teh Q vector components <qx3> and <qy3>, as well as the divergence |
| c of the Q vector <divq3> on the model grid. The grid is specified in the horizontal |
| c by <xmin,ymin,dx,dy,nx,ny>. The number of vertical levels is <nz>. The input field |
| c are: potential temperature <th3>, horizontal wind <u3> and <v3>, pressure <p3>. |
| c The calculation follows the one described in "Weather Analysis, Dusan Djuric" |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real qx3(nx,ny,nz) |
| real qy3(nx,ny,nz) |
| real th3(nx,ny,nz) |
| real u3 (nx,ny,nz) |
| real v3 (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real f2 (nx,ny) |
| real mdv |
| c Physical and numerical parameters |
| real scale1,scale2 |
| parameter (scale1=1.E10,scale2=1.E14) |
| real eps |
| parameter (eps=0.01) |
| real g |
| parameter (g=9.80616) |
| real tzero |
| parameter (tzero=273.16) |
| c Auxiliary variables |
| real dudx(nx,ny,nz) |
| real dudy(nx,ny,nz) |
| real dvdx(nx,ny,nz) |
| real dvdy(nx,ny,nz) |
| real dtdx(nx,ny,nz) |
| real dtdy(nx,ny,nz) |
| integer i,j,k |
| c Needed derivatives |
| call deriv (dudx, u3,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dvdy, v3,'y',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dtdx,th3,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dtdy,th3,'y',z3,x2,y2,nx,ny,nz,mdv) |
| c Calculate vector components of the Q vector |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| c Evaluate Q vector formula with missing data check |
| if ((abs(dudx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dudy(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdy(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdy(i,j,k)-mdv).gt.eps)) then |
| qx3(i,j,k) = -g/tzero * (dudx(i,j,k)*dtdx(i,j,k) |
| > +dvdx(i,j,k)*dtdy(i,j,k)) |
| qy3(i,j,k) = -g/tzero * (dudy(i,j,k)*dtdx(i,j,k) |
| > +dvdy(i,j,k)*dtdy(i,j,k)) |
| else |
| qx3(i,j,k)=mdv |
| qy3(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| c Scale the output |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if (abs(qx3(i,j,k)-mdv).gt.eps) then |
| qx3(i,j,k)=scale1*qx3(i,j,k) |
| endif |
| if (abs(qy3(i,j,k)-mdv).gt.eps) then |
| qy3(i,j,k)=scale1*qy3(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ---------------------------------------------------------------- |
| c Calculate divergence of wind field |
| c ---------------------------------------------------------------- |
| subroutine calc_div_q (div_q,qx,qy, |
| > z3,x2,y2,f2,nx,ny,nz,mdv) |
| c Calculate the divergence (div_q) of the Q vector field if |
| c the components (qx,qy) are given. The grid and the missing data |
| c value are specified by (z3,x2,y2,f2,nx,ny,nz,mdv) |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real div_q(nx,ny,nz) |
| real qx (nx,ny,nz) |
| real qy (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real f2 (nx,ny) |
| real mdv |
| c Physical parameters and numerical constants |
| real eps |
| parameter (eps=0.01) |
| c Auxiliray variables |
| integer i,j,k |
| real rho |
| real dqxdx(nx,ny,nz) |
| real dqydy(nx,ny,nz) |
| c Calculate horizontal derivatives of pressure |
| call deriv(dqxdx,qx,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv(dqydy,qy,'y',z3,x2,y2,nx,ny,nz,mdv) |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(dqxdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dqydy(i,j,k)-mdv).gt.eps)) then |
| div_q(i,j,k)=dqxdx(i,j,k)+dqydy(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c **************************************************************** |
| c * SUBROUTINE SECTION: GRID HANDLING * |
| c **************************************************************** |
| c ----------------------------------------------------------------- |
| c Horizontal and vertical derivatives for 3d fields |
| c ----------------------------------------------------------------- |
| subroutine deriv (df,f,direction, |
| > z3,x2,y2,nx,ny,nz,mdv) |
| c Calculate horizontal and vertical derivatives of the 3d field <f>. |
| c The direction of the derivative is specified in <direction> |
| c 'x','y' : Horizontal derivative in x and y direction |
| c 'p','z','t','m' : Vertical derivative (pressure, height, theta, model) |
| c The 3d field <z3> specifies the isosurfaces along which the horizontal |
| c derivatives are calculated or the levels for the vertical derivatives. |
| implicit none |
| c Input and output parameters |
| integer nx,ny,nz |
| real df (nx,ny,nz) |
| real f (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| character direction |
| real mdv |
| c Numerical and physical parameters |
| real pi180 |
| parameter (pi180=3.141592654/180.) |
| real deltay |
| parameter (deltay=111.1775E3) |
| real zerodiv |
| parameter (zerodiv=0.00000001) |
| real eps |
| parameter (eps=0.01) |
| c Auxiliary variables |
| integer i,j,k |
| real vmin,vmax |
| real scale,lat |
| real vu,vl,vuvl,vlvu |
| integer o,u,w,e,n,s |
| c Vertical derivative |
| if ((direction.eq.'z').or. |
| > (direction.eq.'th').or. |
| > (direction.eq.'p').or. |
| > (direction.eq.'m').and. |
| > (nz.gt.1)) then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| o=k+1 |
| if (o.gt.nz) o=nz |
| u=k-1 |
| if (u.lt.1) u=1 |
| if ((abs(f(i,j,o)-mdv).gt.eps).and. |
| > (abs(f(i,j,u)-mdv).gt.eps).and. |
| > (abs(f(i,j,k)-mdv).gt.eps)) then |
| vu = z3(i,j,k)-z3(i,j,o) |
| vl = z3(i,j,u)-z3(i,j,k) |
| vuvl = vu/(vl+zerodiv) |
| vlvu = 1./(vuvl+zerodiv) |
| df(i,j,k) = 1./(vu+vl) |
| > * (vuvl*(f(i,j,u)-f(i,j,k)) |
| > + vlvu*(f(i,j,k)-f(i,j,o))) |
| else |
| df(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| c Horizontal derivative in the y direction: 3d |
| elseif (direction.eq.'y') then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| s=j-1 |
| if (s.lt.1) s=1 |
| n=j+1 |
| if (n.gt.ny) n=ny |
| if ((abs(f(i,n,k)-mdv).gt.eps).and. |
| > (abs(f(i,j,k)-mdv).gt.eps).and. |
| > (abs(f(i,s,k)-mdv).gt.eps)) then |
| vu = 1000.*(y2(i,j)-y2(i,n)) |
| vl = 1000.*(y2(i,s)-y2(i,j)) |
| vuvl = vu/(vl+zerodiv) |
| vlvu = 1./(vuvl+zerodiv) |
| df(i,j,k) = 1./(vu+vl) |
| > * (vuvl*(f(i,s,k)-f(i,j,k)) |
| > + vlvu*(f(i,j,k)-f(i,n,k))) |
| else |
| df(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| c Horizontal derivative in the x direction: 3d |
| elseif (direction.eq.'x') then |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| w=i-1 |
| if (w.lt.1) w=1 |
| e=i+1 |
| if (e.gt.nx) e=nx |
| if ((abs(f(w,j,k)-mdv).gt.eps).and. |
| > (abs(f(i,j,k)-mdv).gt.eps).and. |
| > (abs(f(e,j,k)-mdv).gt.eps)) then |
| vu = 1000.*(x2(i,j)-x2(e,j)) |
| vl = 1000.*(x2(w,j)-x2(i,j)) |
| vuvl = vu/(vl+zerodiv) |
| vlvu = 1./(vuvl+zerodiv) |
| df(i,j,k) = 1./(vu+vl) |
| > * (vuvl*(f(w,j,k)-f(i,j,k)) |
| > + vlvu*(f(i,j,k)-f(e,j,k))) |
| else |
| df(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| c Undefined direction for derivative |
| else |
| print*,'Invalid direction of derivative... Stop' |
| stop |
| endif |
| end |
| c ----------------------------------------------------------------- |
| c Horizontal filter |
| c ----------------------------------------------------------------- |
| subroutine filt2d (a,af,nx,ny,fil,misdat, |
| & iperx,ipery,ispol,inpol) |
| c Apply a conservative diffusion operator onto the 2d field a, |
| c with full missing data checking. |
| c |
| c a real inp array to be filtered, dimensioned (nx,ny) |
| c af real out filtered array, dimensioned (nx,ny), can be |
| c equivalenced with array a in the calling routine |
| c f1 real workarray, dimensioned (nx+1,ny) |
| c f2 real workarray, dimensioned (nx,ny+1) |
| c fil real inp filter-coeff., 0<afil<=1. Maximum filtering with afil=1 |
| c corresponds to one application of the linear filter. |
| c misdat real inp missing-data value, a(i,j)=misdat indicates that |
| c the corresponding value is not available. The |
| c misdat-checking can be switched off with with misdat=0. |
| c iperx int inp periodic boundaries in the x-direction (1=yes,0=no) |
| c ipery int inp periodic boundaries in the y-direction (1=yes,0=no) |
| c inpol int inp northpole at j=ny (1=yes,0=no) |
| c ispol int inp southpole at j=1 (1=yes,0=no) |
| c |
| c Christoph Schaer, 1993 |
| c argument declaration |
| integer nx,ny |
| real a(nx,ny),af(nx,ny),fil,misdat |
| integer iperx,ipery,inpol,ispol |
| c local variable declaration |
| integer i,j,is |
| real fh |
| real f1(nx+1,ny),f2(nx,ny+1) |
| c compute constant fh |
| fh=0.125*fil |
| c compute fluxes in x-direction |
| if (misdat.eq.0.) then |
| do j=1,ny |
| do i=2,nx |
| f1(i,j)=a(i-1,j)-a(i,j) |
| enddo |
| enddo |
| else |
| do j=1,ny |
| do i=2,nx |
| if ((a(i,j).eq.misdat).or.(a(i-1,j).eq.misdat)) then |
| f1(i,j)=0. |
| else |
| f1(i,j)=a(i-1,j)-a(i,j) |
| endif |
| enddo |
| enddo |
| endif |
| if (iperx.eq.1) then |
| c do periodic boundaries in the x-direction |
| do j=1,ny |
| f1(1,j)=f1(nx,j) |
| f1(nx+1,j)=f1(2,j) |
| enddo |
| else |
| c set boundary-fluxes to zero |
| do j=1,ny |
| f1(1,j)=0. |
| f1(nx+1,j)=0. |
| enddo |
| endif |
| c compute fluxes in y-direction |
| if (misdat.eq.0.) then |
| do j=2,ny |
| do i=1,nx |
| f2(i,j)=a(i,j-1)-a(i,j) |
| enddo |
| enddo |
| else |
| do j=2,ny |
| do i=1,nx |
| if ((a(i,j).eq.misdat).or.(a(i,j-1).eq.misdat)) then |
| f2(i,j)=0. |
| else |
| f2(i,j)=a(i,j-1)-a(i,j) |
| endif |
| enddo |
| enddo |
| endif |
| c set boundary-fluxes to zero |
| do i=1,nx |
| f2(i,1)=0. |
| f2(i,ny+1)=0. |
| enddo |
| if (ipery.eq.1) then |
| c do periodic boundaries in the x-direction |
| do i=1,nx |
| f2(i,1)=f2(i,ny) |
| f2(i,ny+1)=f2(i,2) |
| enddo |
| endif |
| if (iperx.eq.1) then |
| if (ispol.eq.1) then |
| c do south-pole |
| is=(nx-1)/2 |
| do i=1,nx |
| f2(i,1)=-f2(mod(i-1+is,nx)+1,2) |
| enddo |
| endif |
| if (inpol.eq.1) then |
| c do north-pole |
| is=(nx-1)/2 |
| do i=1,nx |
| f2(i,ny+1)=-f2(mod(i-1+is,nx)+1,ny) |
| enddo |
| endif |
| endif |
| c compute flux-convergence -> filter |
| if (misdat.eq.0.) then |
| do j=1,ny |
| do i=1,nx |
| af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1)) |
| enddo |
| enddo |
| else |
| do j=1,ny |
| do i=1,nx |
| if (a(i,j).eq.misdat) then |
| af(i,j)=misdat |
| else |
| af(i,j)=a(i,j)+fh*(f1(i,j)-f1(i+1,j)+f2(i,j)-f2(i,j+1)) |
| endif |
| enddo |
| enddo |
| endif |
| end |
| c -------------------------------------------------------------------------------- |
| c Read refernece profile from netcdf |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol,oro, |
| > pvsrcfile) |
| c Read the reference profile from 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 nx,nny,nz : Grid dimension in x,y,z direction |
| c deltax,deltay,deltaz : Grid spacings used for calculations (m) |
| c coriol : Coriolis parameter (s^-1) |
| c oro : Height of orography (m) |
| c pvsrcfile : Input 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) |
| real deltax,deltay,deltaz |
| real coriol (0:nx,0:ny) |
| real oro (0:nx,0:ny) |
| 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 |
| real x(0:nx,0:ny),y(0:nx,0:ny) |
| real mean,count |
| 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 cdfopn(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 |
| c Read data from netcdf file |
| isok=0 |
| varname='NSQREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,nsqref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='RHOREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,rhoref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='THETAREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,thetaref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='PREREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,pressref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='ZREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,zref,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='CORIOL' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,coriol,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='ORO' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,oro,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='X' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,x,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='Y' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,y,stat) |
| if (stat.ne.0) goto 997 |
| c Close netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| c Determine the grid spacings <deltax, deltay, deltaz> |
| mean=0. |
| count=0. |
| do i=1,nx |
| do j=0,ny |
| mean=mean+abs(x(i)-x(i-1)) |
| count=count+1. |
| enddo |
| enddo |
| deltax=mean/count |
| mean=0. |
| count=0. |
| do j=1,ny |
| do i=0,nx |
| mean=mean+abs(y(j)-y(j-1)) |
| count=count+1. |
| enddo |
| enddo |
| deltay=mean/count |
| mean=0. |
| count=0. |
| do k=1,nz-1 |
| mean=mean+abs(zref(k+1)-zref(k-1)) |
| count=count+1. |
| enddo |
| deltaz=mean/count |
| return |
| c Exception handling |
| 997 print*,'Read_Ref: Problem with input netcdf file... Stop' |
| stop |
| end |
| Property changes: |
| Added: svn:executable |