| 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 |
|---|