| 0,0 → 1,871 |
|---|
| |
| PROGRAM set_boundcon |
| |
| c ************************************************************************ |
| c * Set boundary conditions for inversion; lower and upper boundary * |
| c * conditions for potential temperature; lateral boundary conditions * |
| c * for zonal and meridional wind; in particular, missing data values * |
| c * are removed. * |
| c * * |
| c * Michael Sprenger / Summer 2006 * |
| c ************************************************************************ |
| |
| c -------------------------------------------------------------------------------- |
| c Declaration of variables, parameters, externals and common blocks |
| c -------------------------------------------------------------------------------- |
| |
| implicit none |
| |
| c Input and output file |
| character*80 anomafile |
| character*80 referfile |
| |
| c Grid parameters |
| integer nx,ny,nz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real dx,dy,dz |
| real mdv |
| real deltax,deltay,deltaz |
| real, allocatable,dimension (:,:) :: coriol |
| |
| c Reference state |
| real, allocatable,dimension (:) :: nsqref |
| real, allocatable,dimension (:) :: thetaref |
| real, allocatable,dimension (:) :: rhoref |
| real, allocatable,dimension (:) :: pressref |
| real, allocatable,dimension (:) :: zref |
| |
| C Boundary conditions |
| real, allocatable,dimension (:,:) :: thetatop |
| real, allocatable,dimension (:,:) :: thetabot |
| real, allocatable,dimension (:,:) :: ua |
| real, allocatable,dimension (:,:) :: ub |
| real, allocatable,dimension (:,:) :: va |
| real, allocatable,dimension (:,:) :: vb |
| |
| c 3d arrays |
| real,allocatable,dimension (:,:,:) :: th_anom,pv_anom |
| real,allocatable,dimension (:,:,:) :: uu_anom,vv_anom |
| |
| c Auxiliary variables |
| integer i,j,k,kk |
| integer stat |
| character*80 varname |
| integer n1,n2 |
| |
| c -------------------------------------------------------------------------------- |
| c Input |
| c -------------------------------------------------------------------------------- |
| |
| print*,'********************************************************' |
| print*,'* CHECK_BOUNDCON *' |
| print*,'********************************************************' |
| |
| c Read parameter file |
| open(10,file='fort.10') |
| read(10,*) anomafile |
| read(10,*) referfile |
| close(10) |
| print* |
| print*,trim(anomafile) |
| print*,trim(referfile) |
| print* |
| |
| c Get lat/lon gid parameters from input file |
| call read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > anomafile) |
| print*,'Read_Dim: nx,ny,nz = ',nx,ny,nz |
| print*,' dx,dy,dz = ',dx,dy,dz |
| print*,' xmin,ymin,zmin = ',xmin,ymin,zmin |
| print*,' mdv = ',mdv |
| print* |
| |
| c Count from 0, not from 1: consistent with <inv_cart.f>. |
| nx=nx-1 |
| ny=ny-1 |
| nz=nz-1 |
| |
| c Allocate memory for 3d arrays |
| allocate(pv_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pv_anom' |
| allocate(th_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating th_anom' |
| allocate(uu_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating uu_anom' |
| allocate(vv_anom (0:nx,0:ny,0:nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating vv_anom' |
| |
| c Allocate memory for reference profile |
| allocate(rhoref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating rhoref' |
| allocate(pressref(0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating pressref' |
| allocate(thetaref(0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating thetaref' |
| allocate(nsqref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating nsqref' |
| allocate(zref (0:2*nz),STAT=stat) |
| if (stat.ne.0) print*,'error allocating zref' |
| |
| c Allocate memory for Coriolis parameter |
| allocate(coriol(0:nx,0:ny),STAT=stat) |
| if (stat.ne.0) print*,'error allocating f' |
| |
| c Allocate memory for boundary conditions |
| allocate(thetatop(0:nx,0:ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array thetatop ***' |
| allocate(thetabot(0:nx,0:ny),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array thetabot ***' |
| allocate(ua(0:nx,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ua ***' |
| allocate(ub(0:nx,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array ub ***' |
| allocate(va(0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array va ***' |
| allocate(vb(0:ny,0:nz),stat=stat) |
| if (stat.ne.0) print*,'*** error allocating array vb ***' |
| |
| c Read reference profile and ngrid parameters |
| call read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol, |
| > referfile) |
| deltax=1000.*deltax |
| deltay=1000.*deltay |
| print*,'Deltax,deltay,deltaz =',deltax,deltay,deltaz |
| |
| c Read data from MOD file |
| varname='QGPV' |
| call read_inp (pv_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='TH' |
| call read_inp (th_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='U' |
| call read_inp (uu_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| varname='V' |
| call read_inp (vv_anom,varname,anomafile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| |
| |
| c -------------------------------------------------------------------------------- |
| c Consistency check for boundary conditions and adaptions |
| c -------------------------------------------------------------------------------- |
| |
| c Copy 3d to boundary conditions |
| do i=0,nx |
| do j=0,ny |
| thetatop(i,j)=th_anom(i,j,nz) |
| thetabot(i,j)=th_anom(i,j,0) |
| enddo |
| enddo |
| do i=0,nx |
| do k=0,nz |
| ua(i,k)=uu_anom(i, 0,k) |
| ub(i,k)=uu_anom(i,ny,k) |
| enddo |
| enddo |
| do j=0,ny |
| do k=0,nz |
| va(j,k)=vv_anom( 0,j,k) |
| vb(j,k)=vv_anom(nx,j,k) |
| enddo |
| enddo |
| |
| c Check the lower and upper boundary condition for consistency check |
| print* |
| call combouncon(pv_anom,nsqref,rhoref,thetatop, |
| > thetabot,thetaref,coriol,ua,ub,va,vb, |
| > nx,ny,nz,deltax,deltay,deltaz) |
| print* |
| |
| |
| end |
| |
| |
| c ******************************************************************************** |
| c * NETCDF INPUT AND OUTPUT * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Read input fields for reference profile |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE read_inp (field,fieldname,pvsrcfile, |
| > nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv) |
| |
| c Read <fieldname> from netcdf file <pvsrcfile> into <field>. The grid is specified |
| c by <nx,ny,nz,dx,dy,dz,xmin,ymin,zmin>. A check is performed whether the input |
| c files are consitent with this grid. The missing data value is set to <mdv>. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real field(0:nx,0:ny,0:nz) |
| character*80 fieldname |
| character*80 pvsrcfile |
| real dx,dy,dz,xmin,ymin,zmin |
| real mdv |
| |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer cdfid,stat,cdfid99 |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| integer ndimin,outid,i,j,k |
| real max_th |
| real tmp(nx,ny,nz) |
| integer ntimes |
| real time(200) |
| integer nvars |
| character*80 vnam(100),varname |
| integer isok |
| |
| c Open the input netcdf file |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Check whether needed variables are on file |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 998 |
| isok=0 |
| varname=trim(fieldname) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| |
| c Get the grid parameters from theta |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Check whether grid parameters are consistent |
| if ( (vardim(1).ne.(nx+1)).or. |
| > (vardim(2).ne.(ny+1)).or. |
| > (vardim(3).ne.(nz+1)).or. |
| > (abs(varmin(1)-xmin).gt.eps).or. |
| > (abs(varmin(2)-ymin).gt.eps).or. |
| > (abs(varmin(3)-zmin).gt.eps).or. |
| > (abs((varmax(1)-varmin(1))/real(vardim(1)-1)-dx).gt.eps).or. |
| > (abs((varmax(2)-varmin(2))/real(vardim(2)-1)-dy).gt.eps).or. |
| > (abs((varmax(3)-varmin(3))/real(vardim(3)-1)-dz).gt.eps) ) |
| >then |
| print*,'Input grid inconsitency...' |
| print*,' Nx = ',vardim(1),nx+1 |
| print*,' Ny = ',vardim(2),ny+1 |
| print*,' Nz = ',vardim(3),nz+1 |
| print*,' Varminx = ',varmin(1),xmin |
| print*,' Varminy = ',varmin(2),ymin |
| print*,' Varminz = ',varmin(3),zmin |
| print*,' Dx = ',(varmax(1)-varmin(1))/real(nx-1),dx |
| print*,' Dy = ',(varmax(2)-varmin(2))/real(ny-1),dy |
| print*,' Dz = ',(varmax(3)-varmin(3))/real(nz-1),dz |
| goto 998 |
| endif |
| |
| c Load variables |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 998 |
| call getdat(cdfid,varname,time(1),0,field,stat) |
| print*, 'R ',trim(varname),' ',trim(pvsrcfile) |
| if (stat.ne.0) goto 998 |
| |
| c Close input netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 998 |
| |
| c Set missing data value to <mdv> |
| do i=1,nx |
| do j=1,ny |
| do k=1,nz |
| if (abs(field(i,j,k)-misdat).lt.eps) then |
| field(i,j,k)=mdv |
| endif |
| enddo |
| enddo |
| enddo |
| |
| return |
| |
| c Exception handling |
| 998 print*,'Read_Inp: Problem with input netcdf file... Stop' |
| stop |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Check whether variable is found on netcdf file |
| c -------------------------------------------------------------------------------- |
| |
| subroutine check_varok (isok,varname,varlist,nvars) |
| |
| c Check whether the variable <varname> is in the list <varlist(nvars)>. If this is |
| C the case, <isok> is incremented by 1. Otherwise <isok> keeps its value. |
| |
| implicit none |
| |
| c Declaraion of subroutine parameters |
| integer isok |
| integer nvars |
| character*80 varname |
| character*80 varlist(nvars) |
| |
| c Auxiliary variables |
| integer i |
| |
| c Main |
| do i=1,nvars |
| if (trim(varname).eq.trim(varlist(i))) isok=isok+1 |
| enddo |
| |
| end |
| |
| c -------------------------------------------------------------------------------- |
| c Get grid parameters |
| c -------------------------------------------------------------------------------- |
| |
| subroutine read_dim (nx,ny,nz,dx,dy,dz,xmin,ymin,zmin,mdv, |
| > pvsrcfile) |
| |
| c Get the grid parameters from the variable <THETA> on the input file <pvsrcfile>. |
| c The grid parameters are |
| c nx,ny,nz : Number of grid points in x, y and z direction |
| c xmin,ymin,zmin : Minimum domain coordinates in x, y and z direction |
| c xmax,ymax,zmax : Maximal domain coordinates in x, y and z direction |
| c dx,dy,dz : Horizontal and vertical resolution |
| c Additionally, it is checked whether the vertical grid is equally spaced. If ok, |
| c the grid paramters are transformed from lon/lat to distance (in meters) |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| character*80 pvsrcfile |
| integer nx,ny,nz |
| real dx,dy,dz |
| real xmin,ymin,zmin,xmax,ymax,zmax |
| real mdv |
| |
| c Numerical epsilon and other physical/geoemtrical parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer cdfid,cstid |
| integer ierr |
| character*80 vnam(100),varname |
| integer nvars |
| integer isok |
| integer vardim(4) |
| real misdat |
| real varmin(4),varmax(4),stag(4) |
| real aklev(1000),bklev(1000),aklay(1000),bklay(1000) |
| real dh |
| character*80 csn |
| integer ndim |
| integer i |
| |
| c Get all grid parameters |
| call cdfopn(pvsrcfile,cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getvars(cdfid,nvars,vnam,ierr) |
| if (ierr.ne.0) goto 998 |
| isok=0 |
| varname='QGPV' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 998 |
| call getcfn(cdfid,csn,ierr) |
| if (ierr.ne.0) goto 998 |
| call cdfopn(csn,cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call getdef(cdfid,varname,ndim,misdat,vardim,varmin,varmax, |
| > stag,ierr) |
| if (ierr.ne.0) goto 998 |
| nx=vardim(1) |
| ny=vardim(2) |
| nz=vardim(3) |
| xmin=varmin(1) |
| ymin=varmin(2) |
| zmin=varmin(3) |
| call getlevs(cstid,nz,aklev,bklev,aklay,bklay,ierr) |
| if (ierr.ne.0) goto 998 |
| call getgrid(cstid,dx,dy,ierr) |
| if (ierr.ne.0) goto 998 |
| xmax=varmax(1) |
| ymax=varmax(2) |
| zmax=varmax(3) |
| dz=(zmax-zmin)/real(nz-1) |
| call clscdf(cstid,ierr) |
| if (ierr.ne.0) goto 998 |
| call clscdf(cdfid,ierr) |
| if (ierr.ne.0) goto 998 |
| |
| c Check whether the grid is equallay spaced in the vertical |
| do i=1,nz-1 |
| dh=aklev(i+1)-aklev(i) |
| if (abs(dh-dz).gt.eps) then |
| print*,'Aklev: Vertical grid must be equally spaced... Stop' |
| print*,(aklev(i),i=1,nz) |
| stop |
| endif |
| dh=aklay(i+1)-aklay(i) |
| if (abs(dh-dz).gt.eps) then |
| print*,'Aklay: Vertical grid must be equally spaced... Stop' |
| print*,(aklay(i),i=1,nz) |
| stop |
| endif |
| enddo |
| |
| c Set missing data value |
| mdv=misdat |
| |
| return |
| |
| c Exception handling |
| 998 print*,'Read_Dim: Problem with input netcdf file... Stop' |
| stop |
| |
| end |
| |
| |
| c -------------------------------------------------------------------------------- |
| c Read refernece profile from netcdf |
| c -------------------------------------------------------------------------------- |
| |
| SUBROUTINE read_ref (nsqref,rhoref,thetaref,pressref,zref, |
| > nx,ny,nz,deltax,deltay,deltaz,coriol, |
| > pvsrcfile) |
| |
| c Read the reference profile from file |
| c |
| c thetaref : Reference potential temperature (K) |
| c pressref : Reference pressure (Pa) |
| c rhoref : Reference density (kg/m^3) |
| c nsqref : Stratification (s^-1) |
| c zref : Reference height (m) |
| c nx,nny,nz : Grid dimension in x,y,z direction |
| c deltax,deltay,deltaz : Grid spacings used for calculations (m) |
| c coriol : Coriolis parameter (s^-1) |
| c pvsrcfile : Input file |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real nsqref (0:2*nz) |
| real thetaref(0:2*nz) |
| real rhoref (0:2*nz) |
| real pressref(0:2*nz) |
| real zref (0:2*nz) |
| real deltax,deltay,deltaz |
| real coriol (0:nx,0:ny) |
| character*80 pvsrcfile |
| |
| c Numerical and physical parameters |
| real eps |
| parameter (eps=0.01) |
| |
| c Auxiliary variables |
| integer cdfid,stat |
| integer vardim(4) |
| real misdat |
| integer ndimin |
| real varmin(4),varmax(4),stag(4) |
| integer i,j,k,nf1 |
| integer ntimes |
| real time(200) |
| character*80 vnam(100),varname |
| integer nvars |
| integer isok,ierr |
| real x(0:nx,0:ny),y(0:nx,0:ny) |
| real mean,count |
| |
| c Get grid description from topography |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 997 |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 997 |
| isok=0 |
| varname='ORO' |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdef(cdfid,varname,ndimin,misdat,vardim, |
| > varmin,varmax,stag,stat) |
| if (stat.ne.0) goto 997 |
| time(1)=0. |
| call gettimes(cdfid,time,ntimes,stat) |
| if (stat.ne.0) goto 997 |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Open output netcdf file |
| call cdfopn(pvsrcfile,cdfid,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Create the variable if necessary |
| call getvars(cdfid,nvars,vnam,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Read data from netcdf file |
| isok=0 |
| varname='NSQREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,nsqref,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='RHOREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,rhoref,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='THETAREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,thetaref,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='PREREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,pressref,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='ZREF' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,zref,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='CORIOL' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,coriol,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='X' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,x,stat) |
| if (stat.ne.0) goto 997 |
| |
| isok=0 |
| varname='Y' |
| print*,'R ',trim(varname),' ',trim(pvsrcfile) |
| call check_varok(isok,varname,vnam,nvars) |
| if (isok.eq.0) goto 997 |
| call getdat(cdfid,varname,time(1),0,y,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Close netcdf file |
| call clscdf(cdfid,stat) |
| if (stat.ne.0) goto 997 |
| |
| c Determine the grid spacings <deltax, deltay, deltaz> |
| mean=0. |
| count=0. |
| do i=1,nx |
| do j=0,ny |
| mean=mean+abs(x(i,j)-x(i-1,j)) |
| count=count+1. |
| enddo |
| enddo |
| deltax=mean/count |
| |
| mean=0. |
| count=0. |
| do j=1,ny |
| do i=0,nx |
| mean=mean+abs(y(i,j)-y(i,j-1)) |
| count=count+1. |
| enddo |
| enddo |
| deltay=mean/count |
| |
| mean=0. |
| count=0. |
| do k=1,nz-1 |
| mean=mean+abs(zref(k+1)-zref(k-1)) |
| count=count+1. |
| enddo |
| deltaz=mean/count |
| |
| return |
| |
| c Exception handling |
| 997 print*,'Read_Ref: Problem with input netcdf file... Stop' |
| stop |
| |
| end |
| |
| |
| c ******************************************************************************** |
| c * BOUNDARY CONDITIONS - CONSISTENCY CHECK AND ADAPTIONS * |
| c ******************************************************************************** |
| |
| c -------------------------------------------------------------------------------- |
| c Boundary condition |
| c -------------------------------------------------------------------------------- |
| |
| subroutine combouncon(pv,nsq,rhoref,thetatop, |
| > thetabot,thetaref,coriol, |
| > ua,ub,va,vb,nx,ny,nz,dx,dy,dz) |
| |
| c Evaluate the consistency integrals A.7 from Rene's dissertation. This inegral |
| c is a necessary condition that the von Neumann problem has a unique solution. |
| c Adjust the upper and lower boundary conditions on <thetabot> and <thetatop>, so |
| c that the consitency check is ok. |
| |
| implicit none |
| |
| c Declaration of subroutine parameters |
| integer nx,ny,nz |
| real dx,dy,dz |
| real pv(0:nx,0:ny,0:nz) |
| real nsq(0:2*nz) |
| real rhoref(0:2*nz) |
| real thetatop(0:nx,0:ny) |
| real thetabot(0:nx,0:ny) |
| real thetaref(0:2*nz) |
| real coriol(0:nx,0:ny) |
| real ua(0:nx,0:nz) |
| real ub(0:nx,0:nz) |
| real va(0:ny,0:nz) |
| real vb(0:ny,0:nz) |
| |
| c Numerical and physical parameters |
| real g |
| parameter (g=9.81) |
| |
| c Auxiliary variables |
| integer i,j,k |
| real dxy,dxz,dyz,dxyz |
| real integr,denombot,denomtop,denom |
| real shifttop,shiftbot |
| |
| c Set the area and volume infinitesimals for integration |
| dxy =dx*dy |
| dxz =dx*dz |
| dyz =dy*dz |
| dxyz=dx*dy*dz |
| |
| c Init integration variables |
| integr=0. |
| |
| c Inner |
| do k=1,nz-1 |
| do i=1,nx-1 |
| do j=1,ny-1 |
| integr=integr+dxyz*rhoref(2*k)*pv(i,j,k) |
| enddo |
| enddo |
| enddo |
| |
| c ZY plane |
| do k=1,nz-1 |
| do j=1,ny-1 |
| integr=integr+dyz* |
| > rhoref(2*k)*(dx*pv(0, j,k)+va(j,k)) |
| c |
| integr=integr+dyz* |
| > rhoref(2*k)*(dx*pv(nx,j,k)-vb(j,k)) |
| enddo |
| enddo |
| |
| c ZX plane |
| do k=1,nz-1 |
| do i=1,nx-1 |
| integr=integr+dxz* |
| > rhoref(2*k)*(dy*pv(i,0,k)-ua(i,k)) |
| c |
| integr=integr+dxz* |
| > rhoref(2*k)*(dy*pv(i,ny,k)+ub(i,k)) |
| enddo |
| enddo |
| |
| c XY plane |
| do i=1,nx-1 |
| do j=1,ny-1 |
| integr=integr+dxy*rhoref(0)*( |
| > dz*pv(i,j,0)+coriol(i,j)*g*thetabot(i,j)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dxy*rhoref(2*nz)*( |
| > dz*pv(i,j,nz)-coriol(i,j)*g*thetatop(i,j)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| enddo |
| enddo |
| |
| c X edges |
| do i=1,nx-1 |
| integr=integr+dx* |
| > rhoref(0)*(dyz*pv(i,0,0)- |
| > dz*ua(i,0)+dy*coriol(i,0)*g*thetabot(i,0)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dx* |
| > rhoref(0)*(dyz*pv(i,ny,0)+ |
| > dz*ub(i,0)+dy*coriol(i,ny)*g*thetabot(i,ny)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dx* |
| > rhoref(2*nz)*(dyz*pv(i,0,nz)- |
| > dz*ua(i,nz)-dy*coriol(i,0)*g*thetatop(i,0)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| |
| integr=integr+dx* |
| > rhoref(2*nz)*(dyz*pv(i,ny,nz)+ |
| > dz*ub(i,nz)-dy*coriol(i,ny)*g*thetatop(i,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| enddo |
| |
| c Y edges |
| do j=1,ny-1 |
| integr=integr+dy* |
| > rhoref(0)*(dxz*pv(0,j,0)+ |
| > dz*va(j,0)+dx*coriol(0,j)*g*thetabot(0,j)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dy* |
| > rhoref(0)*(dxz*pv(nx,j,0)- |
| > dz*vb(j,0)+dx*coriol(nx,j)*g*thetabot(nx,j)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+dy* |
| > rhoref(2*nz)*(dxz*pv(0,j,nz)+ |
| > dz*va(j,nz)-dx*coriol(0,j)*g*thetatop(0,j)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+dy* |
| > rhoref(2*nz)*(dxz*pv(nx,j,nz)- |
| > dz*vb(j,nz)-dx*coriol(nx,j)*g*thetatop(nx,j)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| enddo |
| |
| c Z edges |
| do k=1,nz-1 |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(0,0,k)+ |
| > dy*va(0,k)-dx*ua(0,k)) |
| c |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(nx,0,k)- |
| > dy*vb(0,k)-dx*ua(nx,k)) |
| c |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(0,ny,k)+ |
| > dy*va(ny,k)+dx*ub(0,k)) |
| c |
| integr=integr+dz* |
| > rhoref(2*k)*(dxy*pv(nx,ny,k)- |
| > dy*vb(ny,k)+dx*ub(nx,k)) |
| enddo |
| |
| c Points |
| integr=integr+rhoref(0)*(dxyz*pv(0,0,0)+ |
| > dyz*va(0,0)-dxz*ua(0,0)+ |
| > dxy*coriol(0,0)*g*thetabot(0,0)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(0)*(dxyz*pv(nx,0,0)- |
| > dyz*vb(0,0)-dxz*ua(nx,0)+ |
| > dxy*coriol(nx,0)*g*thetabot(nx,0)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(0)*(dxyz*pv(0,ny,0)+ |
| > dyz*va(ny,0)+dxz*ub(0,0)+ |
| > dxy*coriol(0,ny)*g*thetabot(0,ny)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(0)*(dxyz*pv(nx,ny,0)- |
| > dyz*vb(ny,0)+dxz*ub(nx,0)+ |
| > dxy*coriol(nx,ny)*g*thetabot(nx,ny)/ |
| > (nsq(0)*thetaref(0))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(0,0,nz)+ |
| > dyz*va(0,nz)-dxz*ua(0,nz)- |
| > dxy*coriol(0,0)*g*thetatop(0,0)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(nx,0,nz)- |
| > dyz*vb(0,nz)-dxz*ua(nx,nz)- |
| > dxy*coriol(nx,0)*g*thetatop(nx,0)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(0,ny,nz)+ |
| > dyz*va(ny,nz)+dxz*ub(0,nz)- |
| > dxy*coriol(0,ny)*g*thetatop(0,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| integr=integr+rhoref(2*nz)*(dxyz*pv(nx,ny,nz)- |
| > dyz*vb(ny,nz)+dxz*ub(nx,nz)- |
| > dxy*coriol(nx,ny)*g*thetatop(nx,ny)/ |
| > (nsq(2*nz)*thetaref(2*nz))) |
| c |
| |
| c Get the integrals from the reference state at bottom and top |
| denombot=0. |
| denomtop=0. |
| do i=0,nx |
| do j=0,ny |
| denombot=denombot+dxy* |
| > rhoref(0)*coriol(i,j)*g/ |
| > (nsq(0)*thetaref(0)) |
| c |
| denomtop=denomtop+dxy* |
| > rhoref(2*nz)*coriol(i,j)*g/ |
| > (nsq(2*nz)*thetaref(2*nz)) |
| enddo |
| enddo |
| denom=denomtop-denombot |
| |
| c Determine the deviation of potential temperature from reference profile |
| shiftbot=0. |
| shifttop=0. |
| do i=0,nx |
| do j=0,ny |
| shifttop=shifttop+thetatop(i,j) |
| shiftbot=shiftbot+thetabot(i,j) |
| enddo |
| enddo |
| shifttop=shifttop/real((nx+1)*(ny+1)) |
| shiftbot=shiftbot/real((nx+1)*(ny+1)) |
| |
| c Write some information about the consitency integrals |
| print*,'Consistency Check for boundary' |
| print*,' integ = ', integr |
| print*,' denombot = ', denombot |
| print*,' denomtop = ', denomtop |
| print*,' denom = ', denom |
| print*,' theta adjustment = ', integr/denom |
| print*,' theta shift @ top = ', shifttop, |
| > thetaref(2*nz) |
| print*,' theta shift @ bot = ', shiftbot, |
| > thetaref(0) |
| |
| end |
| |
| Property changes: |
|---|
| Added: svn:executable |
|---|