| /trunk/pvin/inv_cart.f |
|---|
| 0,0 → 1,1529 |
| PROGRAM inv_cart |
| c ******************************************************************************** |
| c * CALCULATES FOR A PV AND THETA DISTRIBUTION OTHER PROGNOSTIC FIELDS BY * |
| c * MEANS OF A PV INVERSION * |
| c * * |
| c * Rene Fehlmann 1994 / Code re-organization: Michael Sprenger, 2006 * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Declaration of variables, parameters, externals and common blocks |
| c -------------------------------------------------------------------------------- |
| implicit none |
| c 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 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 Potentiual vorticity and stream function |
| real, allocatable,dimension (:,:,:) :: psi |
| real, allocatable,dimension (:,:,:) :: pv |
| c Auxiliary arrays for numerics |
| real, allocatable,dimension (:,:,:) :: a |
| real, allocatable,dimension (:,:,:) :: b |
| real, allocatable,dimension (:,:,:) :: c |
| c Input parameters |
| character*80 pvsrcfile |
| character*80 referfile |
| c Auxiliary variables |
| integer i,j,k |
| integer stat |
| c -------------------------------------------------------------------------------- |
| c Preparations |
| c -------------------------------------------------------------------------------- |
| print*,'********************************************************' |
| print*,'* INV_CART *' |
| print*,'********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) pvsrcfile |
| read(10,*) referfile |
| close(10) |
| print* |
| print*,trim(pvsrcfile) |
| c Get lat/lon gid parameters from input file |
| call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > pvsrcfile) |
| print* |
| 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 |
| nx=nx-1 |
| ny=ny-1 |
| nz=nz-1 |
| 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 Allocate memory for 3d PV and stream function |
| allocate(psi(0:nx,0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array psi ***' |
| allocate(pv(0:nx,0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array pv ***' |
| c Alllocate memory for matrix elements for inversion operator |
| allocate(a(0:nx,0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array a ***' |
| allocate(b(0:nx,0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array b ***' |
| allocate(c(0:nx,0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array c ***' |
| c Allocate memory for reference profile |
| allocate(nsqref(0:2*nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array nsqref ***' |
| allocate(thetaref(0:2*nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array thetaref ***' |
| allocate(rhoref(0:2*nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array rhoref ***' |
| allocate(pressref(0:2*nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array pressref ***' |
| allocate(zref(0:2*nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array zref ***' |
| c Allocate memory for Coriolis parameter |
| allocate(coriol(0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating coriol' |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| 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 |
| c Read input fields from netcdf |
| call read_inp (pv,thetabot,thetatop,ua,ub,va,vb, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin, |
| > pvsrcfile) |
| c -------------------------------------------------------------------------------- |
| c Perform the inversion |
| c -------------------------------------------------------------------------------- |
| C Init matrix elements for inversion |
| call matrixel(a,b,c,coriol,nx,ny,nz,nsqref,rhoref, |
| > deltax,deltay,deltaz) |
| c Inversion |
| call sor(psi,nsqref,thetatop,thetabot,thetaref,rhoref,coriol, |
| > pv,ua,ub,va,vb,a,b,c,nx,ny,nz,deltax,deltay,deltaz) |
| c -------------------------------------------------------------------------------- |
| c Output |
| c -------------------------------------------------------------------------------- |
| c Write output to netcdf |
| print* |
| call write_out (psi,thetabot,thetatop,ua,ub,va,vb, |
| > nx,ny,nz,deltax,deltay,deltaz, |
| > coriol,thetaref,rhoref,pressref,pvsrcfile) |
| END |
| c ******************************************************************************** |
| c * NETCDF AND ASCII INPUT AND OUTPUT * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Output to netcdf file |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE write_out (psi,thetabot,thetatop,ua,ub,va,vb,nx,ny,nz, |
| > dx,dy,dz,coriol,thetaref,rhoref,pref, |
| > pvsrcfile) |
| c Write the result of the inversion to output netcdf |
| c |
| c psi : streamm function as calculated from inversion |
| c thetabot,thetatop : potential temperature at lower and upper boundary |
| c ua,ub : Zonal wind at western and eastern boundary |
| c va,vb : Meridional wind at southern and northern boundary |
| c nx,ny,nz : Grid dimension in x, y and z direction |
| c dx,dy,dz : Grid resolution |
| c coriol : Coriolis parameter |
| c thetaref : Reference profile of potential temperature |
| c rhoref : Reference profile of density |
| c pref : Reference profile of pressure |
| c pvsrcfile : Name of the output netcdf file |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real psi(0:nx,0:ny,0:nz) |
| real thetatop(0:nx,0:ny) |
| real thetabot(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) |
| character*(30) pvsrcfile |
| real dx,dy,dz |
| real thetaref(0:2*nz) |
| real rhoref(0:2*nz) |
| real pref(0:2*nz) |
| real coriol(0:nx,0:ny) |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| real g |
| parameter (g=9.81) |
| real preunit |
| parameter (preunit=0.01) |
| real kappa |
| parameter (kappa=0.286) |
| 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,jj,kk |
| real tmp(0:nx,0:ny,0:nz) |
| real pr (0:nx,0:ny,0:nz) |
| real th (0:nx,0:ny,0:nz) |
| integer ntimes |
| real time(200) |
| character*80 vnam(100),varname |
| integer nvars |
| integer isok,ierr |
| real meanpsi,meancnt |
| c Get grid description |
| 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='QGPV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 997 |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| c Open output netcdf file |
| call cdfwopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 997 |
| c Create the variable if necessary |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='PSI' |
| 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 997 |
| endif |
| isok=0 |
| varname='U' |
| 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 997 |
| endif |
| isok=0 |
| varname='V' |
| 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 997 |
| endif |
| isok=0 |
| varname='TH' |
| 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 997 |
| endif |
| isok=0 |
| varname='T' |
| 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 997 |
| endif |
| isok=0 |
| varname='P' |
| 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 997 |
| endif |
| c Write stream function |
| varname='PSI' |
| call putdat(cdfid,varname,time(1),0,psi,ierr) |
| if (stat.ne.0) goto 997 |
| print*,'W PSI ',trim(pvsrcfile) |
| c Calculate and write velocity U: keep southern and northern boundary |
| do k=0,nz |
| do i=0,nx |
| do j=1,ny-1 |
| tmp(i,j,k)=(psi(i,j-1,k)-psi(i,j+1,k))/ |
| > (2.*dy*coriol(i,j)) |
| enddo |
| tmp(i,0,k) =ua(i,k) |
| tmp(i,ny,k)=ub(i,k) |
| enddo |
| enddo |
| varname='U' |
| call putdat(cdfid,varname,time(1),0,tmp,ierr) |
| if (stat.ne.0) goto 997 |
| print*,'W U ',trim(pvsrcfile) |
| C Calculate and write velocity V: keep western and eastern boundary |
| do k=0,nz |
| do j=0,ny |
| do i=1,nx-1 |
| tmp(i,j,k)=(psi(i+1,j,k)-psi(i-1,j,k))/ |
| > (2.*dx*coriol(i,j)) |
| enddo |
| tmp(0,j,k)=va(j,k) |
| tmp(nx,j,k)=vb(j,k) |
| enddo |
| enddo |
| varname='V' |
| call putdat(cdfid,varname,time(1),0,tmp,ierr) |
| if (stat.ne.0) goto 997 |
| print*,'W V ',trim(pvsrcfile) |
| c Calculate and write potential temperature: keep lower and upper boundary |
| c Potential temperature is needed for calculation of temperature: keep it |
| do i=0,nx |
| do j=0,ny |
| th(i,j,0)=thetabot(i,j) |
| do k=1,nz-1 |
| th(i,j,k)=thetaref(2*k)* |
| > (psi(i,j,k+1)-psi(i,j,k-1))/(2.*dz*g) |
| enddo |
| th(i,j,nz)=thetatop(i,j) |
| enddo |
| enddo |
| varname='TH' |
| call putdat(cdfid,varname,time(1),0,th,ierr) |
| if (stat.ne.0) goto 997 |
| print*,'W TH ',trim(pvsrcfile) |
| c Calculate and write pressure: The pressure is directly proportional to the |
| c streamfunction. But the streamfunction is determined only up to an additive |
| c constant. Shift the streamfunction in such a way that it vanish in the mean |
| c on the boundaries. Pressure is needed for calculation of temperature: keep it |
| meanpsi=0. |
| meancnt=0. |
| do i=0,nx |
| do j=0,ny |
| meanpsi=meanpsi+psi(i,j,0)+psi(i,j,nz) |
| meancnt=meancnt+2 |
| enddo |
| enddo |
| do i=0,nx |
| do k=0,nz |
| meanpsi=meanpsi+psi(i,0,k)+psi(i,ny,k) |
| meancnt=meancnt+2 |
| enddo |
| enddo |
| do j=0,ny |
| do k=0,nz |
| meanpsi=meanpsi+psi(0,j,k)+psi(nx,j,k) |
| meancnt=meancnt+2 |
| enddo |
| enddo |
| meanpsi=meanpsi/meancnt |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| kk=2*k |
| pr(i,j,k)=preunit*rhoref(kk)*(psi(i,j,k)-meanpsi) |
| enddo |
| enddo |
| enddo |
| varname='P' |
| call putdat(cdfid,varname,time(1),0,pr,ierr) |
| if (stat.ne.0) goto 997 |
| print*,'W P ',trim(pvsrcfile) |
| c Calculate and write temperature |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| kk=2*k |
| tmp(i,j,k)=((pref(kk)/1000.)**kappa) * |
| > (th(i,j,k)+kappa*thetaref(kk)*pr(i,j,k)/pref(kk)) |
| enddo |
| enddo |
| enddo |
| varname='T' |
| call putdat(cdfid,varname,time(1),0,tmp,ierr) |
| if (stat.ne.0) goto 997 |
| print*,'W T ',trim(pvsrcfile) |
| c Close output netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| return |
| c Exception handling |
| 997 print*,'Problem with output netcdf file... Stop' |
| stop |
| end |
| c -------------------------------------------------------------------------------- |
| c Read the reference file |
| 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 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 the input netcdf file |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE read_inp (pv,thetabot,thetatop,ua,ub,va,vb, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin, |
| > pvsrcfile) |
| c Read all needed field from netcdf file <pvsrcfile>. The input fields are: |
| c pv : quasigeostrophic potential vorticity |
| c thetabot,thetatop : potential temperature at lower and upper boundary |
| c ua,ub : Zonal wind at western and eastern boundary |
| c va,vb : Meridional wind at southern and northern boundary |
| c The grid is specified by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed |
| c whether the input files are consitent with this grid. The input netcdf file must |
| c contain the variables <QGPV,THETA,U,V>. If the netcdf file also contains the fields |
| c <DQGPV,DTHETA,DU,DV>, these increments are added to <QGPV,THETA,U,V>. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real pv(0:nx,0:ny,0:nz) |
| real thetatop(0:nx,0:ny) |
| real thetabot(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) |
| character*(30) pvsrcfile |
| real dx,dy,dz,xmin,ymin,zmin |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| 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 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='TH' |
| call check_varok(isok,varname,vnam,nvars) |
| varname='U' |
| call check_varok(isok,varname,vnam,nvars) |
| varname='V' |
| call check_varok(isok,varname,vnam,nvars) |
| varname='QGPV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.ne.4) goto 998 |
| c Get the grid parameters |
| varname='QGPV' |
| 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 consitent |
| 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(nx)-dx).gt.eps).or. |
| > (abs((varmax(2)-varmin(2))/real(ny)-dy).gt.eps).or. |
| > (abs((varmax(3)-varmin(3))/real(nz)-dz).gt.eps) ) then |
| print*,'Input grid inconsitency...' |
| print*,'xmin : ',xmin,varmin(1) |
| print*,'ymin : ',ymin,varmin(2) |
| print*,'zmin : ',zmin,varmin(3) |
| print*,'dx : ',dx,(varmax(1)-varmin(1))/real(nx) |
| print*,'dy : ',dy,(varmax(2)-varmin(2))/real(ny) |
| print*,'dz : ',dz,(varmax(3)-varmin(3))/real(nz) |
| print*,'nx : ',nx |
| print*,'ny : ',ny |
| print*,'nz : ',nz |
| goto 998 |
| endif |
| c THETA: Load upper and lower boundary values |
| varname='TH' |
| call getdat(cdfid,varname,time(1),0,tmp,stat) |
| print*,'R TH ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| do i=0,nx |
| do j=0,ny |
| if ( abs(tmp(i,j,0)-misdat).lt.eps ) then |
| thetabot(i,j)=0. |
| else |
| thetabot(i,j)=tmp(i,j,0) |
| endif |
| if ( abs(tmp(i,j,nz)-misdat).lt.eps ) then |
| thetatop(i,j)=0. |
| else |
| thetatop(i,j)=tmp(i,j,nz) |
| endif |
| enddo |
| enddo |
| c U: Load zonal velocity at southern and northern boundary |
| varname='U' |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| call getdat(cdfid,varname,time(1),0,tmp,stat) |
| print*,'R U ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| do i=0,nx |
| do k=0,nz |
| if ( abs(tmp(i,0,k)-misdat).lt.eps ) then |
| ua(i,k)=0. |
| else |
| ua(i,k)=tmp(i,0,k) |
| endif |
| if ( abs(tmp(i,ny,k)-misdat).lt.eps ) then |
| ub(i,k)=0. |
| else |
| ub(i,k)=tmp(i,ny,k) |
| endif |
| enddo |
| enddo |
| c Load meridional velocity at western and eastern boundary |
| varname='V' |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| call getdat(cdfid,varname,time(1),0,tmp,stat) |
| print*,'R V ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| do j=0,ny |
| do k=0,nz |
| if ( abs(tmp(0,j,k)-misdat).lt.eps ) then |
| va(j,k)=0. |
| else |
| va(j,k)=tmp(0,j,k) |
| endif |
| if ( abs(tmp(nx,j,k)-misdat).lt.eps ) then |
| vb(j,k)=0. |
| else |
| vb(j,k)=tmp(nx,j,k) |
| endif |
| enddo |
| enddo |
| c Load qgPV |
| varname='QGPV' |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| call getdat(cdfid,varname,time(1),0,tmp,stat) |
| print*,'R QGPV ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if ( abs(tmp(i,j,k)-misdat).lt.eps ) then |
| pv(i,j,k)=0. |
| else |
| pv(i,j,k)=tmp(i,j,k) |
| endif |
| enddo |
| enddo |
| enddo |
| c Close input netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| return |
| c Exception handling |
| 998 print*,'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 * INVERSION ROUTINES * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c SOR algorithm (successive over relaxation) |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE sor(psi,nsq,thetatop,thetabot,thetaref,rhoref, |
| > coriol,pv,ua,ub,va,vb,a,b,c,nx,ny,nz,dx,dy,dz) |
| c Solve the qgPV equation by succesive over relaxation (SOR). The subroutine |
| c parameters are: |
| c |
| c psi : Streamfunction, i.e. result of the PV inversion |
| c nsq,rhoref,thetaref : Reference profile |
| c thetatop,thetabot : Upper and lower boundary condition |
| c pv : quasigeostrophic potential vorticity (qgPV) |
| c ua,ub,va,vb : lateral boundary condition for wind |
| c a,b,c : Matrices for the inversion operator |
| c nx,ny,nz,dx,dy,dz : Grid specification |
| c coriol : Coriolis parameter |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real dx,dy,dz |
| real psi (0:nx,0:ny,0:nz) |
| real nsq (0:2*nz) |
| real thetatop(0:nx,0:ny) |
| real thetabot(0:nx,0:ny) |
| real thetaref(0:2*nz) |
| real rhoref (0:2*nz) |
| real pv (0:nx,0:ny,0:nz) |
| real ua (0:nx,0:nz) |
| real ub (0:nx,0:nz) |
| real va (0:ny,0:nz) |
| real vb (0:ny,0:nz) |
| real a (0:nx,0:ny,0:nz) |
| real b (0:nx,0:ny,0:nz) |
| real c (0:nx,0:ny,0:nz) |
| real coriol (0:nx,0:ny) |
| c Numerical and physical parameters |
| real maxspec |
| parameter (maxspec=2.0) |
| integer nofiter |
| parameter (nofiter=500) |
| real omega |
| parameter (omega=1.81) |
| c Auxiliary variables |
| integer counter |
| integer i,j,k |
| real deltasq,psigauge |
| real specx,specy,specz |
| real helpx,helpy,helpz |
| c Init the output array |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| psi(i,j,k)=0. |
| enddo |
| enddo |
| enddo |
| c Calculate the spectrum of the matrix |
| i=nx/2 |
| specx=4.*a(i,0,0)/ |
| > (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1)) |
| specy=2.*(b(i,0,0)+b(i,1,0))/ |
| > (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1)) |
| specz=2.*(c(i,0,0)+c(i,0,1))/ |
| > (2.*a(i,0,0)+b(i,0,0)+b(i,1,0)+c(i,0,0)+c(i,0,1)) |
| do k=1,nz-2 |
| do j=1,ny-2 |
| helpx=4.*a(i,j,k)/ |
| > (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1)) |
| if (helpx.gt.specx) specx=helpx |
| helpy=2.*(b(i,j,k)+b(i,j+1,k))/ |
| > (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1)) |
| if (helpy.gt.specy) specy=helpy |
| helpz=2.*(c(i,j,k)+c(i,j,k+1))/ |
| > (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k)+c(i,j,k-1)) |
| if (helpz.gt.specz) specz=helpz |
| enddo |
| enddo |
| c Check whether the dimensions of the grid are sufficient |
| print * |
| print *, 'Spectrum of the matrix in each direction ' |
| print *, 'Spectrum = ', specx, specy, specz |
| print * |
| if ((maxspec*specx.lt.specy).or.(maxspec*specx.lt.specz)) then |
| print*,' Nx too small... Stop' |
| stop |
| endif |
| if ((maxspec*specy.lt.specx).or.(maxspec*specy.lt.specz)) then |
| print*,'Ny too small... Stop' |
| stop |
| endif |
| if ((maxspec*specz.lt.specx).or.(maxspec*specz.lt.specy)) then |
| print*,'Nz too small... Stop' |
| stop |
| endif |
| c Calculate error: control variable for the iteration |
| psigauge=0. |
| deltasq=0. |
| do k=1,nz-1 |
| do i=1,nx-1 |
| do j=1,ny-1 |
| deltasq=deltasq+(-pv(i,j,k)+( |
| > a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)- |
| > 2.*psi(i,j,k)) + |
| > b(i,j,k)*(psi(i,j+1,k)-psi(i,j,k))- |
| > b(i,j-1,k)*(psi(i,j,k)-psi(i,j-1,k))+ |
| > c(i,j,k)*(psi(i,j,k+1)-psi(i,j,k))- |
| > c(i,j,k-1)*(psi(i,j,k)-psi(i,j,k-1)) |
| > )/(dx*dy*dz*rhoref(2*k)))**2. |
| enddo |
| enddo |
| enddo |
| print 102, 'psigauge', psigauge, 'deltasq', |
| > deltasq/(real(nx)*real(ny)*real(nz)) |
| c Iterations |
| do counter=1,nofiter |
| C Perform one iteration step |
| call psiappsor(omega,pv,psi,nsq,rhoref,thetatop, |
| > thetabot,thetaref,coriol,ua,ub,va,vb, |
| > a,b,c,nx,ny,nz,dx,dy,dz) |
| c Adjustment |
| if (mod(counter,100).eq.0) then |
| psigauge=0. |
| do i=0,nx |
| do j=0,ny |
| if (psi(i,j,0).lt.psigauge) then |
| psigauge=psi(i,j,0) |
| endif |
| enddo |
| enddo |
| do k=0,nz |
| do i=0,nx |
| do j=0,ny |
| psi(i,j,k)=psi(i,j,k)-psigauge |
| enddo |
| enddo |
| enddo |
| endif |
| c Calculate error: control variable for the iteration |
| if (mod(counter,nofiter/10).eq.0) then |
| deltasq=0. |
| do k=1,nz-1 |
| do i=1,nx-1 |
| do j=1,ny-1 |
| deltasq=deltasq+(-pv(i,j,k)+( |
| > a(i,j,k)*(psi(i+1,j,k)+psi(i-1,j,k)- |
| > 2.*psi(i,j,k)) + |
| > b(i,j,k)*(psi(i,j+1,k)-psi(i,j,k))- |
| > b(i,j-1,k)*(psi(i,j,k)-psi(i,j-1,k))+ |
| > c(i,j,k)*(psi(i,j,k+1)-psi(i,j,k))- |
| > c(i,j,k-1)*(psi(i,j,k)-psi(i,j,k-1)) |
| > )/(dx*dy*dz*rhoref(2*k)))**2. |
| enddo |
| enddo |
| enddo |
| print 102, 'psigauge', psigauge, 'deltasq', |
| > deltasq/(real(nx)*real(ny)*real(nz)) |
| endif |
| enddo |
| return |
| c Format specifications |
| 102 format (a11, ' = ',e10.3,a11, ' = ',e10.3) |
| end |
| c -------------------------------------------------------------------------------- |
| c SOR algorithm (successive over relaxation) |
| c -------------------------------------------------------------------------------- |
| subroutine psiappsor(omega,pv,psi,nsq,rhoref,thetatop, |
| > thetabot,thetaref,coriol,ua,ub,va,vb, |
| > a,b,c,nx,ny,nz,dx,dy,dz) |
| c Perform one relaxation step |
| c |
| c psi : Streamfunction, i.e. result of the PV inversion |
| c nsq,rhoref,thetaref : Reference profile |
| c thetatop,thetabot : Upper and lower boundary condition |
| c pv : quasigeostrophic potential vorticity (qgPV) |
| c ua,ub,va,vb : lateral boundary condition for wind |
| c a,b,c : Matrices for the inversion operator |
| c nx,ny,nz,dx,dy,dz : Grid specification |
| c nofiter : Number of iterations |
| c omega : Relaxation parameter |
| c coriol : Coriolis parameter |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real pv(0:nx,0:ny,0:nz) |
| real psi(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 ua(0:nx,0:nz) |
| real ub(0:nx,0:nz) |
| real va(0:ny,0:nz) |
| real vb(0:ny,0:nz) |
| real a(0:nx,0:ny,0:nz) |
| real b(0:nx,0:ny,0:nz) |
| real c(0:nx,0:ny,0:nz) |
| real coriol(0:nx,0:ny) |
| real dx,dy,dz |
| real omega |
| c Numerical and physical parameters |
| real g |
| parameter (g=9.81) |
| c Auxiliary variables |
| integer i,j,k |
| real dxy,dxz,dyz,dxyz |
| c Set the area and volume infinitesimals for integration |
| dxy=dx*dy |
| dxz=dx*dz |
| dyz=dy*dz |
| dxyz=dx*dy*dz |
| c Inner |
| do k=1,nz-1 |
| do i=1,nx-1 |
| do j=1,ny-1 |
| psi(i,j,k)=omega*(-dxyz* |
| > rhoref(2*k)*pv(i,j,k)+a(i,j,k)* |
| > (psi(i+1,j,k)+psi(i-1,j,k))+b(i,j,k)* |
| > psi(i,j+1,k)+b(i,j-1,k)*psi(i,j-1,k)+c(i,j,k)* |
| > psi(i,j,k+1)+c(i,j,k-1)*psi(i,j,k-1))/ |
| > (2.*a(i,j,k)+b(i,j,k)+b(i,j-1,k)+c(i,j,k-1)+ |
| > c(i,j,k))+(1.-omega)*psi(i,j,k) |
| enddo |
| enddo |
| enddo |
| c ZY plane |
| do k=1,nz-1 |
| do j=1,ny-1 |
| psi(0,j,k)=omega*(-dyz* |
| > rhoref(2*k)*(dx*pv(0,j,k)+va(j,k))+ |
| > a(0,j,k)*psi(1,j,k)+ |
| > b(0,j,k)*psi(0,j+1,k)+b(0,j-1,k)*psi(0,j-1,k)+ |
| > c(0,j,k)*psi(0,j,k+1)+c(0,j,k-1)*psi(0,j,k-1))/ |
| > (a(0,j,k)+b(0,j,k)+b(0,j-1,k)+c(0,j,k-1)+c(0,j,k)) |
| > +(1.-omega)*psi(0,j,k) |
| c |
| psi(nx,j,k)=omega*(-dyz* |
| > rhoref(2*k)*(dx*pv(nx,j,k)-vb(j,k))+ |
| > a(nx,j,k)*psi(nx-1,j,k)+ |
| > b(nx,j,k)*psi(nx,j+1,k)+b(nx,j-1,k)*psi(nx,j-1,k)+ |
| > c(nx,j,k)*psi(nx,j,k+1)+c(nx,j,k-1)*psi(nx,j,k-1))/ |
| > (a(nx,j,k)+b(nx,j-1,k)+b(nx,j,k)+c(nx,j,k-1)+c(nx,j,k)) |
| > +(1.-omega)*psi(nx,j,k) |
| enddo |
| enddo |
| c ZX plane |
| do k=1,nz-1 |
| do i=1,nx-1 |
| psi(i,0,k)=omega*(-dxz* |
| > rhoref(2*k)*(dy*pv(i,0,k)-ua(i,k))+ |
| > a(i,0,k)*(psi(i+1,0,k)+psi(i-1,0,k))+ |
| > b(i,0,k)*psi(i,1,k)+ |
| > c(i,0,k)*psi(i,0,k+1)+c(i,0,k-1)*psi(i,0,k-1))/ |
| > (2.*a(i,0,k)+b(i,0,k)+c(i,0,k-1)+c(i,0,k)) |
| > +(1.-omega)*psi(i,0,k) |
| c |
| psi(i,ny,k)=omega*(-dxz* |
| > rhoref(2*k)*(dy*pv(i,ny,k)+ub(i,k))+ |
| > a(i,ny-1,k)*(psi(i+1,ny,k)+psi(i-1,ny,k))+ |
| > b(i,ny-1,k)*psi(i,ny-1,k)+ |
| > c(i,ny-1,k)*psi(i,ny,k+1)+c(i,ny-1,k-1)* |
| > psi(i,ny,k-1))/(2.*a(i,ny-1,k)+b(i,ny-1,k)+ |
| > c(i,ny-1,k-1)+c(i,ny-1,k)) |
| > +(1.-omega)*psi(i,ny,k) |
| enddo |
| enddo |
| c XY plane |
| do i=1,nx-1 |
| do j=1,ny-1 |
| psi(i,j,0)=omega*(-dxy*rhoref(0)*( |
| > dz*pv(i,j,0)+g*coriol(i,j)*thetabot(i,j)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(i,j,0)*(psi(i+1,j,0)+psi(i-1,j,0))+ |
| > b(i,j,0)*psi(i,j+1,0)+b(i,j-1,0)*psi(i,j-1,0)+ |
| > c(i,j,0)*psi(i,j,1))/ |
| > (2.*a(i,j,0)+b(i,j-1,0)+b(i,j,0)+c(i,j,0)) |
| > +(1.-omega)*psi(i,j,0) |
| c |
| psi(i,j,nz)=omega*(-dxy*rhoref(2*nz)*( |
| > dz*pv(i,j,nz)-g*coriol(i,j)*thetatop(i,j)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(i,j,nz)*(psi(i+1,j,nz)+psi(i-1,j,nz))+ |
| > b(i,j,nz)*psi(i,j+1,nz)+b(i,j-1,nz)*psi(i,j-1,nz)+ |
| > c(i,j,nz-1)*psi(i,j,nz-1))/ |
| > (2.*a(i,j,nz)+b(i,j-1,nz)+b(i,j,nz)+c(i,j,nz-1)) |
| > +(1.-omega)*psi(i,j,nz) |
| enddo |
| enddo |
| c Y edges |
| do j=1,ny-1 |
| psi(0,j,0)=omega*(-dy*rhoref(0)*(dxz*pv(0,j,0)+ |
| > dz*va(j,0)+dx*g*coriol(0,j)*thetabot(0,j)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(0,j,0)*psi(1,j,0)+ |
| > b(0,j,0)*psi(0,j+1,0)+b(0,j-1,0)*psi(0,j-1,0)+ |
| > c(0,j,0)*psi(0,j,1))/ |
| > (a(0,j,0)+b(0,j-1,0)+b(0,j,0)+c(0,j,0)) |
| > +(1.-omega)*psi(0,j,0) |
| c |
| psi(nx,j,0)=omega*(-dy*rhoref(0)*(dxz*pv(nx,j,0)- |
| > dz*vb(j,0)+dx*g*coriol(nx,j)*thetabot(nx,j)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(nx,j,0)*psi(nx-1,j,0)+ |
| > b(nx,j,0)*psi(nx,j+1,0)+b(nx,j-1,0)*psi(nx,j-1,0)+ |
| > c(nx,j,0)*psi(nx,j,1))/ |
| > (a(nx,j,0)+b(nx,j-1,0)+b(nx,j,0)+c(nx,j,0)) |
| > +(1.-omega)*psi(nx,j,0) |
| c |
| psi(0,j,nz)=omega*(-dy*rhoref(2*nz)*(dxz*pv(0,j,nz)+ |
| > dz*va(j,nz)-dx*g*coriol(0,j)*thetatop(0,j)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(0,j,nz)*psi(1,j,nz)+ |
| > b(0,j,nz)*psi(0,j+1,nz)+b(0,j-1,nz)*psi(0,j-1,nz)+ |
| > c(0,j,nz-1)*psi(0,j,nz-1))/ |
| > (a(0,j,nz)+b(0,j-1,nz)+b(0,j,nz)+c(0,j,nz-1)) |
| > +(1.-omega)*psi(0,j,nz) |
| c |
| psi(nx,j,nz)=omega*(-dy*rhoref(2*nz)*(dxz*pv(nx,j,nz)- |
| > dz*vb(j,nz)-dx*g*coriol(nx,j)*thetatop(nx,j)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(nx,j,nz)*psi(nx-1,j,nz)+ |
| > b(nx,j,nz)*psi(nx,j+1,nz)+b(nx,j-1,nz)*psi(nx,j-1,nz)+ |
| > c(nx,j,nz-1)*psi(nx,j,nz-1))/ |
| > (a(nx,j,nz)+b(nx,j-1,nz)+b(nx,j,nz)+c(nx,j,nz-1)) |
| > +(1.-omega)*psi(nx,j,nz) |
| enddo |
| c X edges |
| do i=1,nx-1 |
| psi(i,0,0)=omega*(-dx*rhoref(0)*(dyz*pv(i,0,0)- |
| > dz*ua(i,0)+dy*g*coriol(i,0)*thetabot(i,0)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(i,0,0)*(psi(i+1,0,0)+psi(i-1,0,0))+ |
| > b(i,0,0)*psi(i,1,0)+ |
| > c(i,0,0)*psi(i,0,1))/ |
| > (2.*a(i,0,0)+b(i,0,0)+c(i,0,0)) |
| > +(1.-omega)*psi(i,0,0) |
| c |
| psi(i,ny,0)=omega*(-dx*rhoref(0)*(dyz*pv(i,ny,0)+ |
| > dz*ub(i,0)+dy*g*coriol(i,ny)*thetabot(i,ny)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(i,ny,0)*(psi(i+1,ny,0)+psi(i-1,ny,0))+ |
| > b(i,ny-1,0)*psi(i,ny-1,0)+ |
| > c(i,ny,0)*psi(i,ny,1))/ |
| > (2.*a(i,ny,0)+b(i,ny-1,0)+c(i,ny,0)) |
| > +(1.-omega)*psi(i,ny,0) |
| c |
| psi(i,0,nz)=omega*(-dx*rhoref(2*nz)*(dyz*pv(i,0,nz)- |
| > dz*ua(i,nz)-dy*g*coriol(i,0)*thetatop(i,0)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(i,0,nz)*(psi(i+1,0,nz)+psi(i-1,0,nz))+ |
| > b(i,0,nz)*psi(i,1,nz)+ |
| > c(i,0,nz-1)*psi(i,0,nz-1))/ |
| > (2.*a(i,0,nz)+b(i,0,nz)+c(i,0,nz-1)) |
| > +(1.-omega)*psi(i,0,nz) |
| c |
| psi(i,ny,nz)=omega*(-dx*rhoref(2*nz)*(dyz*pv(i,ny,nz)+ |
| > dz*ub(i,nz)-dy*g*coriol(i,ny)*thetatop(i,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(i,ny,nz)*(psi(i+1,ny,nz)+psi(i-1,ny,nz))+ |
| > b(i,ny-1,nz)*psi(i,ny-1,nz)+ |
| > c(i,ny,nz-1)*psi(i,ny,nz-1))/ |
| > (2.*a(i,ny,nz)+b(i,ny-1,nz)+c(i,ny,nz-1)) |
| > +(1.-omega)*psi(i,ny,nz) |
| enddo |
| c Z edges |
| do k=1,nz-1 |
| psi(0,0,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(0,0,k)+ |
| > dy*va(0,k)-dx*ua(0,k))+ |
| > a(0,0,k)*psi(1,0,k)+ |
| > b(0,0,k)*psi(0,1,k)+ |
| > c(0,0,k)*psi(0,0,k+1)+c(0,0,k-1)*psi(0,0,k-1))/ |
| > (a(0,0,k)+b(0,0,k)+c(0,0,k-1)+c(0,0,k)) |
| > +(1.-omega)*psi(0,0,k) |
| c |
| psi(nx,0,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(nx,0,k)- |
| > dy*vb(0,k)-dx*ua(nx,k))+ |
| > a(nx,0,k)*psi(nx-1,0,k)+ |
| > b(nx,0,k)*psi(nx,1,k)+ |
| > c(nx,0,k)*psi(nx,0,k+1)+c(nx,0,k-1)*psi(nx,0,k-1))/ |
| > (a(nx,0,k)+b(nx,0,k)+c(nx,0,k-1)+c(nx,0,k)) |
| > +(1.-omega)*psi(nx,0,k) |
| c |
| psi(0,ny,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(0,ny,k)+ |
| > dy*va(ny,k)+dx*ub(0,k))+ |
| > a(0,ny,k)*psi(1,ny,k)+ |
| > b(0,ny-1,k)*psi(0,ny-1,k)+ |
| > c(0,ny,k)*psi(0,ny,k+1)+c(0,ny,k-1)*psi(0,ny,k-1))/ |
| > (a(0,ny,k)+b(0,ny-1,k)+c(0,ny,k-1)+c(0,ny,k)) |
| > +(1.-omega)*psi(0,ny,k) |
| c |
| psi(nx,ny,k)=omega*(-dz*rhoref(2*k)*(dxy*pv(nx,ny,k)- |
| > dy*vb(ny,k)+dx*ub(nx,k))+ |
| > a(nx,ny,k)*psi(nx-1,ny,k)+ |
| > b(nx,ny-1,k)*psi(nx,ny-1,k)+ |
| > c(nx,ny,k)*psi(nx,ny,k+1)+c(nx,ny,k-1)*psi(nx,ny,k-1))/ |
| > (a(nx,ny,k)+b(nx,ny-1,k)+c(nx,ny,k-1)+c(nx,ny,k)) |
| > +(1.-omega)*psi(nx,ny,k) |
| enddo |
| c Points |
| psi(0,0,0)=omega*(-rhoref(0)*(dxyz*pv(0,0,0)+dyz*va(0,0)- |
| > dxz*ua(0,0)+dxy*g*coriol(0,0)*thetabot(0,0)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(0,0,0)*psi(1,0,0)+ |
| > b(0,0,0)*psi(0,1,0)+ |
| > c(0,0,0)*psi(0,0,1))/ |
| > (a(0,0,0)+b(0,0,0)+c(0,0,0))+ |
| > (1.-omega)*psi(0,0,0) |
| c |
| psi(nx,0,0)=omega*(-rhoref(0)*(dxyz*pv(nx,0,0)-dyz*vb(0,0)- |
| > dxz*ua(nx,0)+dxy*g*coriol(nx,0)*thetabot(nx,0)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(nx,0,0)*psi(nx-1,0,0)+ |
| > b(nx,0,0)*psi(nx,1,0)+ |
| > c(nx,0,0)*psi(nx,0,1))/ |
| > (a(nx,0,0)+b(nx,0,0)+c(nx,0,0))+ |
| > (1.-omega)*psi(nx,0,0) |
| c |
| psi(0,ny,0)=omega*(-rhoref(0)*(dxyz*pv(0,ny,0)+dyz*va(ny,0)+ |
| > dxz*ub(0,0)+dxy*g*coriol(0,ny)*thetabot(0,ny)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(0,ny,0)*psi(1,ny,0)+ |
| > b(0,ny-1,0)*psi(0,ny-1,0)+ |
| > c(0,ny,0)*psi(0,ny,1))/ |
| > (a(0,ny,0)+b(0,ny-1,0)+c(0,ny,0))+ |
| > (1.-omega)*psi(0,ny,0) |
| c |
| psi(nx,ny,0)=omega*(-rhoref(0)*(dxyz*pv(nx,ny,0)-dyz*vb(ny,0)+ |
| > dxz*ub(nx,0)+dxy*g*coriol(nx,ny)*thetabot(nx,ny)/ |
| > (nsq(0)*thetaref(0)))+ |
| > a(nx,ny,0)*psi(nx-1,ny,0)+ |
| > b(nx,ny-1,0)*psi(nx,ny-1,0)+ |
| > c(nx,ny,0)*psi(nx,ny,1))/ |
| > (a(nx,ny,0)+b(nx,ny-1,0)+c(nx,ny,0))+ |
| > (1.-omega)*psi(nx,ny,0) |
| c |
| psi(0,0,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(0,0,nz)+dyz*va(0,nz)- |
| > dxz*ua(0,nz)-dxy*g*coriol(0,0)*thetatop(0,0)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(0,0,nz)*psi(1,0,nz)+ |
| > b(0,0,nz)*psi(0,1,nz)+ |
| > c(0,0,nz-1)*psi(0,0,nz-1))/ |
| > (a(0,0,nz)+b(0,0,nz)+c(0,0,nz-1))+ |
| > (1.-omega)*psi(0,0,nz) |
| c |
| psi(nx,0,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(nx,0,nz)-dyz*vb(0,nz) |
| > -dxz*ua(nx,nz)-dxy*g*coriol(nx,0)*thetatop(nx,0)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(nx,0,nz)*psi(nx-1,0,nz)+ |
| > b(nx,0,nz)*psi(nx,1,nz)+ |
| > c(nx,0,nz-1)*psi(nx,0,nz-1))/ |
| > (a(nx,0,nz)+b(nx,0,nz)+c(nx,0,nz-1))+ |
| > (1.-omega)*psi(nx,0,nz) |
| c |
| psi(0,ny,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(0,ny,nz)+ |
| > dyz*va(ny,nz)+dxz*ub(0,nz)- |
| > dxy*g*coriol(0,ny)*thetatop(0,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(0,ny,nz)*psi(1,ny,nz)+ |
| > b(0,ny-1,nz)*psi(0,ny-1,nz)+ |
| > c(0,ny,nz-1)*psi(0,ny,nz-1))/ |
| > (a(0,ny,nz)+b(0,ny-1,nz)+c(0,ny,nz-1))+ |
| > (1.-omega)*psi(0,ny,nz) |
| c |
| psi(nx,ny,nz)=omega*(-rhoref(2*nz)*(dxyz*pv(nx,ny,nz)- |
| > dyz*vb(ny,nz)+dxz*ub(nx,nz)- |
| > dxy*g*coriol(nx,ny)*thetatop(nx,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz)))+ |
| > a(nx,ny,nz)*psi(nx-1,ny,nz)+ |
| > b(nx,ny-1,nz)*psi(nx,ny-1,nz)+ |
| > c(nx,ny,nz-1)*psi(nx,ny,nz-1))/ |
| > (a(nx,ny,nz)+b(nx,ny-1,nz)+c(nx,ny,nz-1))+ |
| > (1.-omega)*psi(nx,ny,nz) |
| c |
| end |
| c -------------------------------------------------------------------------------- |
| c Init matrix elements for the inversion |
| c -------------------------------------------------------------------------------- |
| subroutine matrixel(a,b,c,coriol,nx,ny,nz,nsq,rhoref,dx,dy,dz) |
| c Define the coefficients for the inversion problem (see page 119ff in Rene's |
| c dissertation). |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,nz,ny |
| real a (0:nx,0:ny,0:nz) |
| real b (0:nx,0:ny,0:nz) |
| real c (0:nx,0:ny,0:nz) |
| real coriol(0:nx,0:ny) |
| real nsq (0:2*nz) |
| real rhoref(0:2*nz) |
| real dx,dy,dz |
| c Auxiliary variables |
| integer i,j,k |
| c Calculate coefficients |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| a(i,j,k)=dy*dz*rhoref(2*k)/(dx*coriol(i,j)) |
| if (j.lt.ny) then |
| b(i,j,k)=dx*dz*rhoref(2*k)/ |
| > (dy*0.5*(coriol(i,j)+coriol(i,j+1))) |
| else |
| b(i,j,k)=0. |
| endif |
| if (k.lt.nz) then |
| c(i,j,k)=dx*dy*rhoref(2*k+1)*coriol(i,j)/ |
| > (dz*nsq(2*k+1)) |
| else |
| c(i,j,k)=0. |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| Property changes: |
| Added: svn:executable |
| /trunk/pvin/inv_cart.m |
|---|
| 0,0 → 1,259 |
| % ------------------------------------------------------------------------- |
| % Load files |
| % ------------------------------------------------------------------------- |
| % Base directory and filename |
| base = '/net/rossby/lhome/sprenger/PV_Inversion_Tool/run/'; |
| folder = base; |
| filename = 'ANO_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'); |
| z_qgpv = cdf_loadV(folder,filename,'QGPV'); |
| z_psi = cdf_loadV(folder,filename,'PSI'); |
| % ------------------------------------------------------------------------- |
| % Plot stream function |
| % ------------------------------------------------------------------------- |
| for ilev=1:z_psi.nz |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set the geographical projection |
| m_proj('Equidistant Cylindrical','long',[z_psi.lonmin z_psi.lonmax],'lat',[z_psi.latmin z_psi.latmax]); |
| % Scale the plotting field for color map |
| fld=squeeze(z_psi.var(ilev,:,:)); |
| c_map = scale_col(-4000:500:4000,fld); |
| % Plot Stream Function PSI |
| lat=z_psi.ymin + (0:z_psi.ny-1) * z_psi.dy; |
| lon=z_psi.xmin + (0:z_psi.nx-1) * z_psi.dx; |
| [C,h]=m_contourf(lon,lat,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,1); |
| 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 |
| m_grid; |
| m_coast('linewidth',1,'color','k'); |
| title(num2str(z_psi.aklay(ilev))); |
| % Save figure |
| pre=''; |
| if ( z_psi.aklay(ilev) < 10 ) |
| pre='0'; |
| end |
| if ( z_psi.aklay(ilev) < 100 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_psi.aklay(ilev) < 1000 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_psi.aklay(ilev) < 10000 ) |
| pre=[pre '0' ]; |
| end |
| figname = [ base '/sf_2d_' filename '_' pre num2str(z_psi.aklay(ilev)) '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % End loop over all levels |
| end |
| % ------------------------------------------------------------------------- |
| % Plot anomaly of potential temperature |
| % ------------------------------------------------------------------------- |
| for ilev=1:z_th.nz |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set the geographical projection |
| m_proj('Equidistant Cylindrical','long',[z_th.lonmin z_th.lonmax],'lat',[z_th.latmin z_th.latmax]); |
| % Scale the plotting field for color map |
| fld=squeeze(z_th.var(ilev,:,:)); |
| c_map = scale_col(-20:2.5:20,fld); |
| % Plot potential temperature |
| 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]=m_contourf(lon,lat,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,1); |
| 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 |
| m_grid; |
| m_coast('linewidth',1,'color','k'); |
| title(num2str(z_th.aklay(ilev))); |
| % Save figure |
| pre=''; |
| if ( z_th.aklay(ilev) < 10 ) |
| pre='0'; |
| end |
| if ( z_th.aklay(ilev) < 100 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_th.aklay(ilev) < 1000 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_th.aklay(ilev) < 10000 ) |
| pre=[pre '0' ]; |
| end |
| figname = [ base '/ta_2d_' filename '_' pre num2str(z_th.aklay(ilev)) '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % End loop over all levels |
| end |
| % ------------------------------------------------------------------------- |
| % Plot anomaly of meridional wind |
| % ------------------------------------------------------------------------- |
| for ilev=1:z_vv.nz |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set the geographical projection |
| m_proj('Equidistant Cylindrical','long',[z_vv.lonmin z_vv.lonmax],'lat',[z_vv.latmin z_vv.latmax]); |
| % Scale the plotting field for color map |
| fld=squeeze(z_vv.var(ilev,:,:)); |
| c_map = scale_col(-30:3:30,fld); |
| % Plot potential temperature |
| lat=z_vv.ymin + (0:z_vv.ny-1) * z_vv.dy; |
| lon=z_vv.xmin + (0:z_vv.nx-1) * z_vv.dx; |
| [C,h]=m_contourf(lon,lat,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,1); |
| 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 |
| m_grid; |
| m_coast('linewidth',1,'color','k'); |
| title(num2str(z_vv.aklay(ilev))); |
| % Save figure |
| pre=''; |
| if ( z_vv.aklay(ilev) < 10 ) |
| pre='0'; |
| end |
| if ( z_vv.aklay(ilev) < 100 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_vv.aklay(ilev) < 1000 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_vv.aklay(ilev) < 10000 ) |
| pre=[pre '0' ]; |
| end |
| figname = [ base '/va_2d_' filename '_' pre num2str(z_vv.aklay(ilev)) '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % End loop over all levels |
| end |
| % ------------------------------------------------------------------------- |
| % Plot anomaly of zonal wind |
| % ------------------------------------------------------------------------- |
| for ilev=1:z_uu.nz |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set the geographical projection |
| m_proj('Equidistant Cylindrical','long',[z_uu.lonmin z_uu.lonmax],'lat',[z_uu.latmin z_uu.latmax]); |
| % Scale the plotting field for color map |
| fld=squeeze(z_uu.var(ilev,:,:)); |
| c_map = scale_col(-30:3:30,fld); |
| % Plot potential temperature |
| lat=z_uu.ymin + (0:z_uu.ny-1) * z_uu.dy; |
| lon=z_uu.xmin + (0:z_uu.nx-1) * z_uu.dx; |
| [C,h]=m_contourf(lon,lat,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,1); |
| 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 |
| m_grid; |
| m_coast('linewidth',1,'color','k'); |
| title(num2str(z_uu.aklay(ilev))); |
| % Save figure |
| pre=''; |
| if ( z_uu.aklay(ilev) < 10 ) |
| pre='0'; |
| end |
| if ( z_uu.aklay(ilev) < 100 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_uu.aklay(ilev) < 1000 ) |
| pre=[pre '0' ]; |
| end |
| if ( z_uu.aklay(ilev) < 10000 ) |
| pre=[pre '0' ]; |
| end |
| figname = [ base '/ua_2d_' filename '_' pre num2str(z_uu.aklay(ilev)) '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % End loop over all levels |
| end |
| Property changes: |
| Added: svn:executable |
| /trunk/pvin/prep_iteration.f |
|---|
| 0,0 → 1,513 |
| PROGRAM prep_iteration |
| c ************************************************************************ |
| c * Prepare the next step for the PV inversion * |
| c * Michael Sprenger / Summer, Autumn 2006 * |
| c ************************************************************************ |
| implicit none |
| c ------------------------------------------------------------------------ |
| c Declaration of variables and parameters |
| c ------------------------------------------------------------------------ |
| c Input and output file |
| character*80 anomafile |
| character*80 iterafile |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real mdv |
| c Numerical epsilon and other variables |
| real eps |
| parameter (eps=0.01) |
| real alpha |
| c 3d arrays |
| real,allocatable,dimension (:,:,:) :: v_iter,v_anom |
| real,allocatable,dimension (:,:,:) :: u_iter,u_anom |
| real,allocatable,dimension (:,:,:) :: t_iter,t_anom |
| real,allocatable,dimension (:,:,:) :: p_iter,p_anom |
| c Auciliary variables |
| integer i,j,k |
| integer stat |
| character*80 varname |
| character*80 name |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| print*,'********************************************************' |
| print*,'* PREP_ITERATION *' |
| print*,'********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) iterafile |
| read(10,*) anomafile |
| read(10,*) name,alpha |
| if ( trim(name).ne.'ALPHA' ) stop |
| close(10) |
| print* |
| print*,trim(anomafile) |
| print*,trim(iterafile) |
| 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(u_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating u_anom' |
| allocate(v_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating v_anom' |
| allocate(t_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating t_anom' |
| allocate(p_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating p_anom' |
| allocate(u_iter (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating u_iter' |
| allocate(v_iter (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating v_iter' |
| allocate(t_iter (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating t_iter' |
| allocate(p_iter (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating p_iter' |
| c Read anomaly and iteration fields |
| varname='U' |
| call read_inp (u_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='V' |
| call read_inp (v_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='T' |
| call read_inp (t_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='P' |
| call read_inp (p_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='U' |
| call read_inp (u_iter,varname,iterafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='V' |
| call read_inp (v_iter,varname,iterafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='T' |
| call read_inp (t_iter,varname,iterafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='P' |
| call read_inp (p_iter,varname,iterafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| c -------------------------------------------------------------------------------- |
| c Modify the iteration fields |
| c -------------------------------------------------------------------------------- |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| c Update zonal velocity |
| if ((abs(u_anom(i,j,k)-mdv).gt.eps).and. |
| > (abs(u_iter(i,j,k)-mdv).gt.eps)) then |
| u_iter(i,j,k)=u_iter(i,j,k)-alpha*u_anom(i,j,k) |
| else |
| u_iter(i,j,k)=mdv |
| endif |
| c Update meridional velocity |
| if ((abs(v_anom(i,j,k)-mdv).gt.eps).and. |
| > (abs(v_iter(i,j,k)-mdv).gt.eps)) then |
| v_iter(i,j,k)=v_iter(i,j,k)-alpha*v_anom(i,j,k) |
| else |
| v_iter(i,j,k)=mdv |
| endif |
| c Update temperature |
| if ((abs(t_anom(i,j,k)-mdv).gt.eps).and. |
| > (abs(t_iter(i,j,k)-mdv).gt.eps)) then |
| t_iter(i,j,k)=t_iter(i,j,k)-alpha*t_anom(i,j,k) |
| else |
| t_iter(i,j,k)=mdv |
| endif |
| c Update pressure |
| if ((abs(p_anom(i,j,k)-mdv).gt.eps).and. |
| > (abs(p_iter(i,j,k)-mdv).gt.eps)) then |
| p_iter(i,j,k)=p_iter(i,j,k)-alpha*p_anom(i,j,k) |
| else |
| p_iter(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| c -------------------------------------------------------------------------------- |
| c Write output |
| c -------------------------------------------------------------------------------- |
| varname='U' |
| call write_inp (u_iter,varname,iterafile,nx,ny,nz) |
| varname='V' |
| call write_inp (v_iter,varname,iterafile,nx,ny,nz) |
| varname='T' |
| call write_inp (t_iter,varname,iterafile,nx,ny,nz) |
| varname='P' |
| call write_inp (p_iter,varname,iterafile,nx,ny,nz) |
| end |
| c ******************************************************************************** |
| c * NETCDF INPUT AND OUTPUT * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Write input field to netcdf |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz) |
| c Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified |
| c by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input |
| c files are consitent with this grid. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field (0:nx,0:ny,0:nz) |
| character*80 fieldname |
| character*80 pvsrcfile |
| c Auxiliary variables |
| integer cdfid,stat |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| integer ndimin,outid,i,j,k |
| real 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 |
| 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 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='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 |
| Property changes: |
| Added: svn:executable |
| /trunk/pvin/pv_to_qgpv.f |
|---|
| 0,0 → 1,894 |
| PROGRAM pv_to_qgpv |
| c ******************************************************************************** |
| c * TRANSFORM ERTEL'S PV TO QUASI-GEOSTROPHIC PV * |
| c * Rene Fehlmann 1994 / Code re-organization: Michael Sprenger, 2006 * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Declaration of variables, parameters, externals and common blocks |
| c -------------------------------------------------------------------------------- |
| implicit none |
| c Input and output file |
| character*80 pvsrcfile |
| character*80 referfile |
| character*80 anomafile |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real mdv |
| c Numerical and physical parameters |
| real pi180 ! Pi/180 |
| parameter (pi180=3.141592654/180.) |
| real rerd ! Earth's radius |
| parameter (rerd=6.371229e6) |
| real eps ! Numerical epsilon |
| parameter (eps=0.01) |
| real scale ! Scale for PV unit |
| parameter (scale=1e6) |
| real minagl ! No PV and qgPV below this height AGL |
| parameter (minagl=1000.) |
| c Reference state and grid parameters |
| real, allocatable,dimension (:) :: nsqref |
| real, allocatable,dimension (:) :: thetaref |
| real, allocatable,dimension (:) :: rhoref |
| real, allocatable,dimension (:) :: pressref |
| real, allocatable,dimension (:) :: zref |
| real, allocatable,dimension (:,:) :: coriol |
| real, allocatable,dimension (:,:) :: oro |
| real deltax,deltay,deltaz |
| c 3d fields for calculation of qgPV and Ertel's PV |
| real,allocatable,dimension (:,:,:) :: qgpv,pv1,pv2,pv |
| c Auxiliary variables |
| real zpos |
| integer i,j,k |
| integer stat |
| character*80 varname |
| integer istep |
| real mean,rmsq,min,max |
| integer step |
| real,allocatable,dimension (:,:) :: tmp |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| print*,'********************************************************' |
| print*,'* PV_TO_QGPV *' |
| print*,'********************************************************' |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) pvsrcfile |
| read(10,*) referfile |
| read(10,*) anomafile |
| close(10) |
| print* |
| print*,trim(pvsrcfile) |
| print*,trim(referfile) |
| print*,trim(anomafile) |
| print* |
| c Get lat/lon gid parameters from input file |
| call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > pvsrcfile) |
| print*,'Read_Dim: nx,ny,nz = ',nx,ny,nz |
| print*,' dx,dy,dz = ',dx,dy,dz |
| print*,' xmin,ymin,zmin = ',xmin,ymin,zmin |
| print*,' mdv = ',mdv |
| print* |
| c Count from 0, not from 1: consistent with <inv_cart.f>. |
| nx=nx-1 |
| ny=ny-1 |
| nz=nz-1 |
| c Allocate memory for reference profile and grid parameters |
| 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(coriol (0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating coriol' |
| allocate(oro (0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating oro' |
| c Allocate memory for calculation of qgPV and Ertel's PV |
| allocate(pv1 (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv1' |
| allocate(pv2 (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv2' |
| allocate(pv (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv' |
| allocate(qgpv(0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating qgpv' |
| c Allocate memory for temporary array |
| allocate(tmp(0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating tmp' |
| c -------------------------------------------------------------------------------- |
| c Calculate the qgPV from Ertel's PV and put it onto file |
| c -------------------------------------------------------------------------------- |
| c Read data from file |
| varname='PV' |
| call read_inp (pv1,varname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='PV_AIM' |
| call read_inp (pv2,varname,pvsrcfile, |
| > 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,oro, |
| > referfile) |
| c If the PV is negative, set it to zero |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if (pv1(i,j,k).lt.0.) pv1(i,j,k)=0. |
| if (pv2(i,j,k).lt.0.) pv2(i,j,k)=0. |
| enddo |
| enddo |
| enddo |
| c Get the difference of Ertel's PV and set all missing values to 0 |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| if ((abs(pv1(i,j,k)-mdv).gt.eps).and. |
| > (abs(pv2(i,j,k)-mdv).gt.eps)) then |
| pv(i,j,k)=pv1(i,j,k)-pv2(i,j,k) |
| else |
| pv(i,j,k)=0. |
| endif |
| enddo |
| enddo |
| enddo |
| c Calculate qgPV |
| call epv_to_qgpv (qgpv,pv, |
| > rhoref,pressref,nsqref,thetaref, |
| > nx,ny,nz,mdv) |
| c Set values on the boundaries to zero |
| do i=0,nx |
| do j=0,ny |
| qgpv(i,j, 0)=0. |
| qgpv(i,j,nz)=0. |
| enddo |
| enddo |
| do i=0,nx |
| do k=0,nz |
| qgpv(i, 0,k)=0. |
| qgpv(i,ny,k)=0. |
| enddo |
| enddo |
| do j=0,ny |
| do k=0,nz |
| qgpv( 0,j,k)=0. |
| qgpv(nx,j,k)=0. |
| enddo |
| enddo |
| c Set all values to zero which are too near to the surface |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| zpos=zmin+real(k)*dz |
| if (zpos.lt.(oro(i,j)+minagl)) then |
| pv(i,j,k)=0. |
| qgpv(i,j,k)=0. |
| endif |
| enddo |
| enddo |
| enddo |
| c Write result to netcdf file |
| varname='QGPV' |
| call write_inp (qgpv,varname,anomafile,nx,ny,nz) |
| varname='PV' |
| call write_inp (pv,varname,anomafile,nx,ny,nz) |
| c Write some info |
| print* |
| print*,'PV -> qgPV: k z min max mean rmsq' |
| step=nz/10 |
| if (step.lt.1) step=1 |
| do k=0,nz,step |
| do i=0,nx |
| do j=0,ny |
| tmp(i,j)=qgpv(i,j,k) |
| enddo |
| enddo |
| call calc_error(min,max,mean,rmsq,tmp,nx,ny) |
| write(*,'(8x,i3,f10.1,4f10.2)') |
| > k,zmin+real(k)*dz,scale*min,scale*max, |
| > scale*mean,scale*rmsq |
| enddo |
| print* |
| c -------------------------------------------------------------------------------- |
| c Format specifications |
| c -------------------------------------------------------------------------------- |
| 111 format (5f20.9) |
| 106 format (2f20.9) |
| end |
| c ******************************************************************************** |
| c * NETCDF INPUT AND OUTPUT * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Write input field to netcdf |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE write_inp (field,fieldname,pvsrcfile,nx,ny,nz) |
| c Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified |
| c by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input |
| c files are consitent with this grid. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field (0:nx,0:ny,0:nz) |
| character*80 fieldname |
| character*80 pvsrcfile |
| c Auxiliary variables |
| integer cdfid,stat |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| integer ndimin,outid,i,j,k |
| real tmp(0:nx,0:ny,0:nz) |
| integer ntimes |
| real time(200) |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| c Get grid parameters from PV |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname='PV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 998 |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| c Save variables (write definition, if necessary) |
| call cdfwopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname=fieldname |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) then |
| call putdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| endif |
| call putdat(cdfid,varname,time(1),0,field,stat) |
| print*,'W ',trim(varname),' ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| c Close input netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| return |
| c Exception handling |
| 998 print*,'Write_Inp: Problem with input netcdf file... Stop' |
| stop |
| end |
| c -------------------------------------------------------------------------------- |
| c Read input fields for reference profile |
| c -------------------------------------------------------------------------------- |
| SUBROUTINE read_inp (field,fieldname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| c Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified |
| c by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input |
| c files are consitent with this grid. The missing data value is set to <mdv>. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| character*80 fieldname |
| character*80 pvsrcfile |
| real dx,dy,dz,xmin,ymin,zmin |
| real mdv |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| c Auxiliary variables |
| integer cdfid,stat,cdfid99 |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| integer ndimin,outid,i,j,k |
| real tmp(nx,ny,nz) |
| integer ntimes |
| real time(200) |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| c Open the input netcdf file |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| c Check whether needed variables are on file |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname=trim(fieldname) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| c Get the grid parameters from theta |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 998 |
| c Check whether grid parameters are consistent |
| if ( (vardim(1).ne.(nx+1)).or. |
| > (vardim(2).ne.(ny+1)).or. |
| > (vardim(3).ne.(nz+1)).or. |
| > (abs(varmin(1)-xmin).gt.eps).or. |
| > (abs(varmin(2)-ymin).gt.eps).or. |
| > (abs(varmin(3)-zmin).gt.eps).or. |
| > (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or. |
| > (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or. |
| > (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) |
| >then |
| print*,'Input grid inconsitency...' |
| print*,' Nx = ',vardim(1),nx+1 |
| print*,' Ny = ',vardim(2),ny+1 |
| print*,' Nz = ',vardim(3),nz+1 |
| print*,' Varminx = ',varmin(1),xmin |
| print*,' Varminy = ',varmin(2),ymin |
| print*,' Varminz = ',varmin(3),zmin |
| print*,' Dx = ',(varmax(1)-varmin(1))/real(nx-1),dx |
| print*,' Dy = ',(varmax(2)-varmin(2))/real(ny-1),dy |
| print*,' Dz = ',(varmax(3)-varmin(3))/real(nz-1),dz |
| goto 998 |
| endif |
| c Load variables |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| call getdat(cdfid,varname,time(1),0,field,stat) |
| print*, 'R ',trim(varname),' ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| c Close input netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| c Set missing data value to <mdv> |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if (abs(field(i,j,k)-misdat).lt.eps) then |
| field(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| return |
| c Exception handling |
| 998 print*,'Read_Inp: Problem with input netcdf file... Stop' |
| stop |
| end |
| c -------------------------------------------------------------------------------- |
| c 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 |
| c -------------------------------------------------------------------------------- |
| c Check whether variable is found on netcdf file |
| c -------------------------------------------------------------------------------- |
| subroutine check_varok (isok,varname,varlist,nvars) |
| c Check whether the variable <varname> is in the list <varlist(nvars)>. If this is |
| C the case, <isok> is incremented by 1. Otherwise <isok> keeps its value. |
| implicit none |
| c Declaraion of subroutine parameters |
| integer isok |
| integer nvars |
| character*80 varname |
| character*80 varlist(nvars) |
| c Auxiliary variables |
| integer i |
| c Main |
| do i=1,nvars |
| if (trim(varname).eq.trim(varlist(i))) isok=isok+1 |
| enddo |
| end |
| c -------------------------------------------------------------------------------- |
| c Get grid parameters |
| c -------------------------------------------------------------------------------- |
| subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > pvsrcfile) |
| c Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>. |
| c The grid parameters are |
| c nx,ny,nz : Number of grid points in x, y and z direction |
| c xmin,ymin,zmin : Minimum domain coordinates in x, y and z direction |
| c xmax,ymax,zmax : Maximal domain coordinates in x, y and z direction |
| c dx,dy,dz : Horizontal and vertical resolution |
| c Additionally, it is checked whether the vertical grid is equally spaced. If ok, |
| c the grid paramters are transformed from lon/lat to distance (in meters) |
| implicit none |
| c Declaration of subroutine parameters |
| character*80 pvsrcfile |
| integer nx,ny,nz |
| real dx,dy,dz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real mdv |
| c Numerical epsilon and other physical/geoemtrical parameters |
| real eps |
| parameter (eps=0.01) |
| c Auxiliary variables |
| integer cdfid,cstid |
| integer ierr |
| character*80 vnam(100),varname |
| integer nvars |
| integer isok |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| real aklev(1000),bklev(1000),aklay(1000),bklay(1000) |
| real dh |
| character*80 csn |
| integer ndim |
| integer i |
| c Get all grid parameters |
| call cdfopn(pvsrcfile,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname='PV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| call getcfn(cdfid,csn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(csn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax, |
| > stag,ierr) |
| if (ierr.ne.0) goto 998 |
| nx=vardim(1) |
| ny=vardim(2) |
| nz=vardim(3) |
| xmin=varmin(1) |
| ymin=varmin(2) |
| zmin=varmin(3) |
| call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr) |
| if (ierr.ne.0) goto 998 |
| call getgrid(cstid,dx,dy,ierr) |
| if (ierr.ne.0) goto 998 |
| xmax=varmax(1) |
| ymax=varmax(2) |
| zmax=varmax(3) |
| dz=(zmax-zmin)/real(nz-1) |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c Check whether the grid is equallay spaced in the vertical |
| do i=1,nz-1 |
| dh=aklev(i+1)-aklev(i) |
| if (abs(dh-dz).gt.eps) then |
| print*,'Aklev: Vertical grid must be equally spaced... Stop' |
| print*,(aklev(i),i=1,nz) |
| stop |
| endif |
| dh=aklay(i+1)-aklay(i) |
| if (abs(dh-dz).gt.eps) then |
| print*,'Aklay: Vertical grid must be equally spaced... Stop' |
| print*,(aklay(i),i=1,nz) |
| stop |
| endif |
| enddo |
| c Set missing data value |
| mdv=misdat |
| return |
| c Exception handling |
| 998 print*,'Read_Dim: Problem with input netcdf file... Stop' |
| stop |
| end |
| c ******************************************************************************** |
| c * CALCULATION * |
| c ******************************************************************************** |
| c -------------------------------------------------------------------------------- |
| c Calculate qgPV from Ertels's PV |
| c -------------------------------------------------------------------------------- |
| subroutine epv_to_qgpv (qgpv,pv, |
| > rhoref,pressref,nsqref,thetaref, |
| > nx,ny,nz,mdv) |
| c Calculate the qgPV from Ertel's PV according to equation 2.11 p16, Thesis |
| c from Rene Fehlmann. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real qgpv(0:nx,0:ny,0:nz),pv(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 mdv |
| c Numerical epsilon |
| 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 |
| c Calculation |
| do i=0,nx |
| do j=0,ny |
| do k=0,nz |
| kk=2*k |
| if (( abs(rhoref(kk) -mdv).gt.eps).and. |
| > ( abs(thetaref(kk)-mdv).gt.eps).and. |
| > ( abs(nsqref(kk) -mdv).gt.eps).and. |
| > ( abs(pv(i,j,k) -mdv).gt.eps)) then |
| qgpv(i,j,k)=rhoref(kk)*g*pv(i,j,k)/thetaref(kk)/ |
| > nsqref(kk)/scale |
| else |
| qgpv(i,j,k)=0. |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c -------------------------------------------------------------------------------- |
| c Calculate error statistics |
| c -------------------------------------------------------------------------------- |
| subroutine calc_error (min,max,mean,rmsq,tmp,nx,ny) |
| c Calculate the error statistics for the two-dimensional error field <tmp>. The |
| c following error measures are calculated: the minimum <min>, the maximum <max>, |
| c the mean <mean>, the root-mean square <rmsq> |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny |
| real tmp(0:nx,0:ny) |
| real mean,rmsq |
| real min,max |
| c Auxiliary variables |
| integer i,j |
| real sum |
| integer cnt |
| c Calculate the minimum and maximum |
| min=tmp(0,0) |
| max=tmp(0,0) |
| do i=0,nx |
| do j=0,ny |
| if (tmp(i,j).lt.min) min=tmp(i,j) |
| if (tmp(i,j).gt.max) max=tmp(i,j) |
| enddo |
| enddo |
| c Calculate the mean |
| sum=0. |
| cnt=0 |
| do i=0,nx |
| do j=0,ny |
| cnt=cnt+1 |
| sum=sum+tmp(i,j) |
| enddo |
| enddo |
| if (cnt.ge.1) then |
| mean=sum/real(cnt) |
| else |
| mean=0. |
| endif |
| c Calculate rmsq |
| rmsq=0. |
| cnt=0 |
| do i=0,nx |
| do j=0,ny |
| cnt=cnt+1 |
| rmsq=rmsq+(tmp(i,j)-mean)**2 |
| enddo |
| enddo |
| if (cnt.ge.1) then |
| rmsq=1./real(cnt)*sqrt(rmsq) |
| else |
| rmsq=0. |
| endif |
| end |
| Property changes: |
| Added: svn:executable |
| /trunk/pvin/pv_to_qgpv.m |
|---|
| 0,0 → 1,81 |
| folder ='/lhome/sprenger/PV_Inversion_Tool/'; |
| filename='Z1_20060115_18'; |
| % -------------------------------------------------------------------------- |
| % Plot qgPV anomaly (horizontal sections) |
| % -------------------------------------------------------------------------- |
| % 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) |
| m_pv = cdf_loadV(folder,filename,'QGPV_ANOM'); |
| % Loop over all levels |
| for ilev=1:m_pv.nz |
| % Create a new figure |
| close; |
| fh=figure('Units','pixels','Position',[100 100 900 900]) |
| % Set the geographical projection |
| m_proj('Equidistant Cylindrical','long',[m_pv.lonmin m_pv.lonmax],'lat',[m_pv.latmin m_pv.latmax]); |
| % Scale the plotting field for color map |
| fld=1e6*squeeze(m_pv.var(ilev,:,:)); |
| c_map = scale_col(0:50:600,fld); |
| % Plot PV |
| lat=m_pv.ymin + (0:m_pv.ny-1) * m_pv.dy; |
| lon=m_pv.xmin + (0:m_pv.nx-1) * m_pv.dx; |
| [C,h]=m_contourf(lon,lat,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,1); |
| caxis(c_map.caxis); |
| q=colorbar('hori'); |
| set(q,'xtick',c_map.xtick,'XTickLabel',c_map.label); |
| % Add the grid and the coast lines to the plot |
| m_grid; |
| m_coast('linewidth',1,'color','k'); |
| title(num2str(m_pv.aklay(ilev))); |
| % Save figure |
| pre=''; |
| if ( m_pv.aklay(ilev) < 10 ) |
| pre='0'; |
| end |
| if ( m_pv.aklay(ilev) < 100 ) |
| pre=[pre '0' ]; |
| end |
| if ( m_pv.aklay(ilev) < 1000 ) |
| pre=[pre '0' ]; |
| end |
| if ( m_pv.aklay(ilev) < 10000 ) |
| pre=[pre '0' ]; |
| end |
| figname = [ base '/qg_2d_' filename '_' pre num2str(m_pv.aklay(ilev)) '.eps' ]; |
| set(gcf, 'PaperPosition', [2 1 15 10]); |
| print('-depsc2','-r0',figname); |
| % End loop over all levels |
| end |
| Property changes: |
| Added: svn:executable |
| /trunk/pvin/z2s.f |
|---|
| 0,0 → 1,978 |
| PROGRAM z2s |
| c Calculate secondary fields on z levels |
| c Michael Sprenger / Summer 2006 |
| implicit none |
| 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 |
| real,allocatable, dimension (:,:,:) :: out |
| real,allocatable, dimension (:,:,:) :: inp |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| integer cdfid,cstid |
| 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 (:,:) :: out2 |
| character*80 vnam2(100) |
| integer nvars2 |
| c --------------------------------------------------------------- |
| c Preparations |
| c --------------------------------------------------------------- |
| print*,'*********************************************************' |
| print*,'* z2s *' |
| 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 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.'TH' ).and. |
| > (fieldname.ne.'NSQ') .and. |
| > (fieldname.ne.'PV' ) .and. |
| > (fieldname.ne.'RHO')) then |
| print*,'Variable not supported ',trim(fieldname) |
| stop |
| endif |
| c Set dependencies |
| if (fieldname.eq.'TH') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| else if (fieldname.eq.'RHO') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| else if (fieldname.eq.'NSQ') then |
| nfields=2 |
| field_nam(1)='T' |
| field_nam(2)='Q' |
| else if (fieldname.eq.'PV') then |
| nfields=4 |
| field_nam(1)='T' |
| field_nam(2)='P' |
| field_nam(3)='U' |
| field_nam(4)='V' |
| endif |
| c Allocate memory |
| allocate(field_dat(nfields,nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating field_dat' |
| allocate(out(nx,ny,nz),stat=stat) |
| if (stat.ne.0) print*,'error allocating out' |
| 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' |
| c Allocate auxiliary fields |
| allocate(out2(nx,ny),stat=stat) |
| if (stat.ne.0) print*,'error allocating out2' |
| 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 Read Coriolis parametzer |
| varname='CORIOL' |
| 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,f2,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'R ',trim(varname),' ',trim(gri) |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| c Init height levels |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| z3(i,j,k)=aklay(k) |
| 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 Change unit of pressure from hPa to Pa |
| 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 |
| c ---------------------------------------------------------------- |
| c Calculations |
| c ---------------------------------------------------------------- |
| c Call to the defining routines |
| if (fieldname.eq.'RHO') then |
| print*,'C ',trim(fieldname) |
| call calc_rho (out, ! RHO |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! P |
| > nx,ny,nz,mdv) |
| else if (fieldname.eq.'TH') then |
| print*,'C ',trim(fieldname) |
| call calc_theta (out, ! TH |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! P |
| > nx,ny,nz,mdv) |
| else if (fieldname.eq.'NSQ') then |
| print*,'C ',trim(fieldname) |
| call calc_nsq (out, ! NSQ |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! Q |
| > z3, ! Z |
| > nx,ny,nz,mdv) |
| else if (fieldname.eq.'PV') then |
| print*,'C ',trim(fieldname) |
| call calc_pv (out, ! PV |
| > field_dat(1,:,:,:), ! T |
| > field_dat(2,:,:,:), ! P |
| > field_dat(3,:,:,:), ! U |
| > field_dat(4,:,:,:), ! V |
| > 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 |
| out2(i,j)=out(i,j,k) |
| enddo |
| enddo |
| do n=1,nfilt |
| call filt2d (out2,out2,nx,ny,1.,mdv,0,0,0,0) |
| enddo |
| do i=1,nx |
| do j=1,ny |
| out(i,j,k)=out2(i,j) |
| enddo |
| enddo |
| enddo |
| c ---------------------------------------------------------------- |
| c Save result onto netcdf file |
| c ---------------------------------------------------------------- |
| call cdfwopn(ofn,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname=fieldname |
| 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,out,ierr) |
| if (ierr.ne.0) goto 998 |
| print*,'W ',trim(varname),' ',trim(ofn) |
| 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 density |
| c ----------------------------------------------------------------- |
| subroutine calc_rho (rho3,t3,p3, |
| > nx,ny,nz,mdv) |
| c Calculate the density <rho3> (in kg/m^3) from temperature <t3> |
| c (in deg C) and pressure <p3> (in hPa). |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real mdv |
| real rho3(nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real p3 (nx,ny,nz) |
| c Declaration of physical constants |
| real eps |
| parameter (eps=0.01) |
| real rd |
| parameter (rd=287.) |
| real tzero |
| parameter (tzero=273.15) |
| c Auxiliary variables |
| integer i,j,k |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(t3(i,j,k)-mdv).gt.eps).and. |
| > (abs(p3(i,j,k)-mdv).gt.eps)) then |
| rho3(i,j,k)=1./rd*p3(i,j,k)/(t3(i,j,k)+tzero) |
| else |
| rho3(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ----------------------------------------------------------------- |
| c Calculate potential temperature |
| c ----------------------------------------------------------------- |
| subroutine calc_theta (th3,t3,p3, |
| > nx,ny,nz,mdv) |
| c Calculate potential temperature <th3> (in K) from temperature <t3> |
| c (in deg C) and pressure <p3> (in hPa). |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real mdv |
| real th3(nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real p3 (nx,ny,nz) |
| c Declaration of physical constants |
| real rdcp,tzero,p0 |
| parameter (rdcp=0.286) |
| parameter (tzero=273.15) |
| parameter (p0=100000.) |
| real eps |
| parameter (eps=0.01) |
| c Auxiliary variables |
| integer i,j,k |
| c Calculation |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(t3(i,j,k)-mdv).gt.eps).and. |
| > (abs(p3(i,j,k)-mdv).gt.eps)) then |
| th3(i,j,k)=(t3(i,j,k)+tzero)*( (p0/p3(i,j,k))**rdcp ) |
| else |
| th3(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ----------------------------------------------------------------- |
| c Calculate stratification |
| c ----------------------------------------------------------------- |
| subroutine calc_nsq (nsq3,t3,q3, |
| > z3,nx,ny,nz,mdv) |
| c Calculate stratification <nsq3> on the model grid. The grid is |
| c specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number |
| c of vertical levels is <nz>. The input field are: temperature <t3>, |
| c specific humidity <q3>, height <z3> and horizontal grid <x2>, <y2>. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real nsq3(nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real q3 (nx,ny,nz) |
| real z3 (nx,ny,nz) |
| real x2 (nx,ny) |
| real y2 (nx,ny) |
| real mdv |
| c Physical and numerical parameters |
| real scale |
| parameter (scale=1.) |
| 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 cp |
| parameter (cp=1005.) |
| real zerodiv |
| parameter (zerodiv=0.0000000001) |
| real R |
| parameter (R=287.) |
| c Auxiliary variables |
| real tv (nx,ny,nz) |
| real dtvdz(nx,ny,nz) |
| integer i,j,k |
| real scaledz |
| c Calculate 3d virtual temperature |
| do k=1,nz |
| do i=1,nx |
| do j=1,ny |
| if ((abs(t3(i,j,k)-mdv).gt.eps).and. |
| > (abs(q3(i,j,k)-mdv).gt.eps)) |
| > then |
| tv(i,j,k) = (t3(i,j,k)+tzero)*(1.+kappa*q3(i,j,k)) |
| else |
| tv(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| c Vertical derivative of virtual temperature |
| call deriv (dtvdz,tv,'z',z3,x2,y2,nx,ny,nz,mdv) |
| c Stratification |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if ((abs(dtvdz(i,j,k)-mdv).gt.eps).and. |
| > (abs(tv (i,j,k)-mdv).gt.eps)) |
| > then |
| nsq3(i,j,k) = g/tv(i,j,k) * (dtvdz(i,j,k) + g/cp) |
| else |
| nsq3(i,j,k) = mdv |
| endif |
| enddo |
| enddo |
| enddo |
| end |
| c ----------------------------------------------------------------- |
| c Calculate potential vorticity |
| c ----------------------------------------------------------------- |
| subroutine calc_pv (pv3,t3,p3,u3,v3, |
| > z3,x2,y2,f2,nx,ny,nz,mdv) |
| c Calculate potential vorticity <pv3> on the model grid. The grid is |
| c specified in the horizontal by <xmin,ymin,dx,dy,nx,ny>. The number |
| c of vertical levels is <nz>. The input field are: potential temperature |
| c <th3>, horizontal wind <u3> and <v3>, density <rho3>, height <z3>. |
| c The terms including the vertical velocity are neglected in the calculation |
| c of the PV. |
| implicit none |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real pv3 (nx,ny,nz) |
| real t3 (nx,ny,nz) |
| real p3 (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 scale |
| parameter (scale=1.E6) |
| real omega |
| parameter (omega=7.292E-5) |
| real pi180 |
| parameter (pi180=3.141592654/180.) |
| real eps |
| parameter (eps=0.01) |
| c Auxiliary variables |
| real dtdz(nx,ny,nz) |
| real dtdx(nx,ny,nz) |
| real dtdy(nx,ny,nz) |
| real dudz(nx,ny,nz) |
| real dvdz(nx,ny,nz) |
| real dvdx(nx,ny,nz) |
| real dudy(nx,ny,nz) |
| real rho3(nx,ny,nz) |
| real th3 (nx,ny,nz) |
| integer i,j,k |
| c Calculate density and potential temperature |
| call calc_rho (rho3,t3,p3,nx,ny,nz,mdv) |
| call calc_theta (th3 ,t3,p3,nx,ny,nz,mdv) |
| c Calculate needed derivatives |
| call deriv (dudz, u3,'z',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dvdz, v3,'z',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dtdz,th3,'z',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) |
| call deriv (dvdx, v3,'x',z3,x2,y2,nx,ny,nz,mdv) |
| call deriv (dudy, u3,'y',z3,x2,y2,nx,ny,nz,mdv) |
| c Calculate potential vorticity |
| do j=1,ny |
| do i=1,nx |
| do k=1,nz |
| c Evaluate PV formula with missing data check |
| if ((abs(dtdz(i,j,k)-mdv).gt.eps).and. |
| > (abs(dudz(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdy(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdz(i,j,k)-mdv).gt.eps).and. |
| > (abs(rho3(i,j,k)-mdv).gt.eps).and. |
| > (abs(dtdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dvdx(i,j,k)-mdv).gt.eps).and. |
| > (abs(dudy(i,j,k)-mdv).gt.eps)) then |
| pv3(i,j,k)=scale*1./rho3(i,j,k)* |
| > ((dvdx(i,j,k)-dudy(i,j,k)+f2(i,j))*dtdz(i,j,k) |
| > +dudz(i,j,k)*dtdy(i,j,k) |
| > -dvdz(i,j,k)*dtdx(i,j,k)) |
| else |
| pv3(i,j,k)=mdv |
| 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 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 |
| Property changes: |
| Added: svn:executable |