Commit cbd85c2a authored by kbk's avatar kbk
Browse files

Poor Man's Z-coordinates

parent d7e26b9f
!$Id: coordinates.F90,v 1.4 2003-09-02 14:45:46 kbk Exp $
!$Id: coordinates.F90,v 1.5 2004-01-05 13:23:27 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -6,21 +6,23 @@
! !ROUTINE: coordinates() - defines the vertical coordinate system.
!
! !INTERFACE:
subroutine coordinates(cord_type,cord_relax)
subroutine coordinates(cord_type,cord_relax,maxdepth)
!
! !DESCRIPTION:
!
! !USES:
use halo_zones, only: update_3d_halo,wait_halo,H_TAG,HU_TAG,HV_TAG
use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HU,HV,az,au,av,min_depth
use domain, only: ga,ddu,ddl,d_gamma
use domain, only: iimin,iimax,jjmin,jjmax,kmax,H,HU,HV,az,au,av,min_depth
use domain, only: ga,ddu,ddl,d_gamma,gamma_surf
use variables_3d, only: dt,kmin,kumin,kvmin,ho,hn,huo,hun,hvo,hvn
use variables_3d, only: kmin_pmz,kumin_pmz,kvmin_pmz
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: cord_type
REALTYPE, intent(in) :: cord_relax
REALTYPE, intent(in) :: maxdepth
!
! !INPUT/OUTPUT PARAMETERS:
!
......@@ -30,7 +32,10 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: coordinates.F90,v $
! Revision 1.4 2003-09-02 14:45:46 kbk
! Revision 1.5 2004-01-05 13:23:27 kbk
! Poor Man's Z-coordinates
!
! Revision 1.4 2003/09/02 14:45:46 kbk
! calculate in HALO-zones instead of using update_3d_halo()
!
! Revision 1.3 2003/04/23 12:16:27 kbk
......@@ -92,12 +97,12 @@
! initial import into CVS
!
! !LOCAL VARIABLES:
integer :: i,j,k,rc,kk
REALTYPE:: tmp,kmaxm1,alpha,maxdepth
LOGICAL :: gamma_surf=.true.
REALTYPE:: HH,zz,r
logical, save :: first=.true.,equiv_sigma=.false.
REALTYPE, save, dimension(:), allocatable :: dga,be,sig
integer :: i,j,k,rc,kk
REALTYPE :: tmp,kmaxm1,alpha
REALTYPE :: HH,zz,r,hnmin=0.01
logical, save :: first=.true.,equiv_sigma=.false.
logical :: kminset
REALTYPE, save, dimension(:), allocatable :: dga,be,sig,zlev
REALTYPE, save, dimension(:,:,:), allocatable :: gga
!EOP
!-----------------------------------------------------------------------
......@@ -110,7 +115,7 @@
if (first) then
first = .false.
allocate(ga(0:kmax),stat=rc)
if (.not. allocated(ga)) allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
select case (cord_type)
case (1) ! sigma coordinates
......@@ -124,6 +129,8 @@
else
! Non-equidistant sigma coordinates
! This zooming routine is from Antoine Garapon, ICCH, DK
if (ddu .lt. _ZERO_) ddu=_ZERO_
if (ddl .lt. _ZERO_) ddl=_ZERO_
allocate(dga(0:kmax),stat=rc)
if (rc /= 0) STOP 'coordinates: Error allocating (dga)'
ga(0)= -_ONE_
......@@ -137,9 +144,78 @@
kmin=1
kumin=1
kvmin=1
kmin_pmz=1
kumin_pmz=1
kvmin_pmz=1
case (2) ! z-level
STDERR 'z-level not coded yet'
stop 'coordinates'
allocate(zlev(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (zlev)'
allocate(gga(I3DFIELD),stat=rc) ! dimensionless gamma-coordinate
if (rc /= 0) stop 'coordinates: Error allocating memory (gga)'
! Calculate profile of zlevels
zlev(0)=-maxdepth
do k=1,kmax
zlev(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
zlev(k)=zlev(k)/(tanh(ddl)+tanh(ddu)) - _ONE_
zlev(k)=zlev(k)*maxdepth
end do
kmin_pmz=1
do j=jjmin,jjmax
do i=iimin,iimax
HH=max(sseo(i,j)+H(i,j),min_depth)
kminset=.true.
if (az(i,j).eq.0) then
do k=1,kmax
hn(i,j,k)=HH/float(kmax)
gga(i,j,k)=hn(i,j,k)/HH
end do
kmin_pmz(i,j)=1
else
k=kmax
111 if (zlev(k-1)-(k-1)*hnmin.ge.-H(i,j)) then
hn(i,j,k)=zlev(k)-zlev(k-1)
else
hn(i,j,k)=max(zlev(k)-(-H(i,j)+(k-1)*hnmin),hnmin)
if (kminset.and.hn(i,j,k).eq.hnmin) then
kmin_pmz(i,j)=k+1
kminset=.false.
end if
end if
gga(i,j,k)=hn(i,j,k)/HH
k=k-1
if (k.ge.1) goto 111
end if
end do
end do
do i=iimin-1,iimax
do j=jjmin,jjmax
do k=1,kmax
huo(i,j,k)=min(hn(i,j,k),hn(i+1,j,k))
hun(i,j,k)=huo(i,j,k)
end do
kumin_pmz(i,j)=max(kmin_pmz(i,j),kmin_pmz(i+1,j))
end do
end do
do i=iimin,iimax
do j=jjmin-1,jjmax
do k=1,kmax
hvo(i,j,k)=min(hn(i,j,k),hn(i+1,j,k))
hvn(i,j,k)=hvo(i,j,k)
end do
kvmin_pmz(i,j)=max(kmin_pmz(i,j),kmin_pmz(i,j+1))
end do
end do
do k=1,kmax
do j=jjmin-1,jjmax
do i=iimin,iimax
hvo(i,j,k)=min(hn(i,j,k),hn(i+1,j,k))
hvn(i,j,k)=hvo(i,j,k)
end do
end do
end do
kmin=1
kumin=1
kvmin=1
case (3) ! general vertical coordinates
do k=0,kmax
ga(k) = k
......@@ -152,6 +228,8 @@
if (rc /= 0) stop 'coordinates: Error allocating memory (gga)'
be(0)= -_ONE_
sig(0)= -_ONE_
if (ddu .lt. _ZERO_) ddu=_ZERO_
if (ddl .lt. _ZERO_) ddl=_ZERO_
do k=1,kmax
be(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
be(k)=be(k)/(tanh(ddl)+tanh(ddu))-_ONE_
......@@ -164,7 +242,7 @@
end if
do j=jjmin-HALO,jjmax+HALO
do i=iimin-HALO,iimax+HALO
HH=max(H(i,j),min_depth)
HH=max(sseo(i,j)+H(i,j),min_depth)
alpha=min(&
((be(kk)-be(kk-1))-D_gamma/HH&
*(sig(kk)-sig(kk-1)))&
......@@ -187,6 +265,9 @@
kmin=1
kumin=1
kvmin=1
kmin_pmz=1
kumin_pmz=1
kvmin_pmz=1
! Here, the initial layer distribution is calculated.
do k=1,kmax
......@@ -219,7 +300,9 @@
end do
end do
end do
case default
end select
end if ! first
......@@ -272,42 +355,96 @@
end do
end do
end if
case (2) ! z-level
! Maybe use mask and not loop over k, e.g. ho(i,j,:)=hn(i,j,:)
!kbk#define USE_ARRAY_SYNTAX
do j=jjmin-HALO,jjmax+HALO
do i=iimin-HALO,iimax+HALO
HH=ssen(i,j)+H(i,j)
#ifdef USE_ARRAY_SYNTAX
ho(i,j,:)=hn(i,j,:)
hn(i,j,:)=gga(i,j,:)*HH
#else
do k=1,kmax
ho(i,j,k)=hn(i,j,k)
hn(i,j,k)=gga(i,j,k)*HH
end do
#endif
end do
end do
do j=jjmin-HALO,jjmax+HALO
do i=iimin-HALO,iimax+HALO-1
HH=ssun(i,j)+HU(i,j)
#ifdef USE_ARRAY_SYNTAX
huo(i,j,:)=hun(i,j,:)
hun(i,j,:)=0.5*(gga(i,j,:)+gga(i+1,j,:))*HH
#else
do k=1,kmax
huo(i,j,k)=hun(i,j,k)
hun(i,j,k)=0.5*(gga(i,j,k)+gga(i+1,j,k))*HH
end do
#endif
end do
end do
do j=jjmin-HALO,jjmax+HALO-1
do i=iimin-HALO,iimax+HALO
HH=ssvn(i,j)+HV(i,j)
#ifdef USE_ARRAY_SYNTAX
hvo(i,j,:)=hvn(i,j,:)
hvn(i,j,:)=0.5*(gga(i,j,:)+gga(i,j+1,:))*HH
#else
do k=1,kmax
hvo(i,j,k)=hvn(i,j,k)
hvn(i,j,k)=0.5*(gga(i,j,k)+gga(i,j+1,k))*HH
end do
#endif
#undef USE_ARRAY_SYNTAX
end do
end do
case (3) ! general vertical coordinates
! The general vertical coordinates can be relaxed towards the new layer
! thicknesses by the following relaxation time scale r. This should
! later be generalised also for sigma coordinates.
maxdepth=600. ! needs to be calculated later ... HB
do j=jjmin-HALO,jjmax+HALO
do i=iimin-HALO,iimax+HALO
r = cord_relax/dt*H(i,j)/maxdepth
zz = -H(i,j)
do k=1,kmax-1
ho(i,j,k) = hn(i,j,k)
HH=max(ssen(i,j)+H(i,j),min_depth)
hn(i,j,k)=(ho(i,j,k)*r+HH*(gga(i,j,k)-gga(i,j,k-1)))/(r+1.)
zz = zz + hn(i,j,k)
end do
ho(i,j,kmax) = hn(i,j,kmax)
hn(i,j,kmax)=ssen(i,j)-zz
r=cord_relax/dt*H(i,j)/maxdepth
HH=ssen(i,j)+H(i,j)
if (HH .lt. D_gamma) then
do k=1,kmax
ho(i,j,k)=hn(i,j,k)
hn(i,j,k)=HH/float(kmax)
end do
else
zz=-H(i,j)
do k=1,kmax-1
ho(i,j,k)=hn(i,j,k)
hn(i,j,k)=(ho(i,j,k)*r+HH*(gga(i,j,k)-gga(i,j,k-1)))/(r+1.)
zz=zz+hn(i,j,k)
end do
ho(i,j,kmax)=hn(i,j,kmax)
hn(i,j,kmax)=ssen(i,j)-zz
end if
end do
end do
do j=jjmin-HALO,jjmax+HALO
do i=iimin-HALO,iimax+HALO-1
!KBK if (au(i,j) .gt. 0) then
r = cord_relax/dt*HU(i,j)/maxdepth
zz = -HU(i,j)
r=cord_relax/dt*HU(i,j)/maxdepth
zz=-HU(i,j)
HH=ssun(i,j)+HU(i,j)
do k=1,kmax-1
huo(i,j,k) = hun(i,j,k)
HH=max(ssun(i,j)+HU(i,j),min_depth)
huo(i,j,k)=hun(i,j,k)
hun(i,j,k)=(huo(i,j,k)*r+HH*0.5*(gga(i,j,k)-gga(i,j,k-1) &
+gga(i+1,j,k)-gga(i+1,j,k-1)))/(r+1.)
zz = zz + hun(i,j,k)
zz=zz+hun(i,j,k)
end do
huo(i,j,kmax) = hun(i,j,kmax)
huo(i,j,kmax)=hun(i,j,kmax)
hun(i,j,kmax)=ssun(i,j)-zz
!KBK end if
end do
......@@ -316,16 +453,16 @@
do j=jjmin-HALO,jjmax+HALO-1
do i=iimin-HALO,iimax+HALO
!KBK if (av(i,j).gt.0) then
r = cord_relax/dt*HV(i,j)/maxdepth
zz = -HV(i,j)
r=cord_relax/dt*HV(i,j)/maxdepth
zz=-HV(i,j)
HH=ssvn(i,j)+HV(i,j)
do k=1,kmax-1
hvo(i,j,k) = hvn(i,j,k)
HH=max(ssvn(i,j)+HV(i,j),min_depth)
hvo(i,j,k)=hvn(i,j,k)
hvn(i,j,k)=(hvo(i,j,k)*r+HH*0.5*(gga(i,j,k)-gga(i,j,k-1) &
+gga(i,j+1,k)-gga(i,j+1,k-1)))/(r+1.)
zz=zz+hvn(i,j,k)
end do
hvo(i,j,kmax) = hvn(i,j,kmax)
hvo(i,j,kmax)=hvn(i,j,kmax)
hvn(i,j,kmax)=ssvn(i,j)-zz
!KBK end if
end do
......
!$Id: m3d.F90,v 1.10 2004-01-02 13:54:24 kbk Exp $
!$Id: m3d.F90,v 1.11 2004-01-05 13:23:27 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -19,7 +19,7 @@
!
! !USES:
use parameters, only: avmmol
use domain, only: vert_cord
use domain, only: maxdepth,vert_cord
use m2d, only: Am
use variables_2d, only: D,z,UEx,VEx
#ifndef NO_BAROCLINIC
......@@ -50,7 +50,10 @@
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: m3d.F90,v $
! Revision 1.10 2004-01-02 13:54:24 kbk
! Revision 1.11 2004-01-05 13:23:27 kbk
! Poor Man's Z-coordinates
!
! Revision 1.10 2004/01/02 13:54:24 kbk
! read equation of state info from namelist - Ruiz
!
! Revision 1.9 2003/12/16 17:02:44 kbk
......@@ -244,7 +247,7 @@
! Needed for interpolation of temperature and salinity
if (.not. hotstart) then
call start_macro()
call coordinates(vert_cord,cord_relax)
call coordinates(vert_cord,cord_relax,maxdepth)
end if
#ifndef NO_BAROCLINIC
......@@ -312,7 +315,7 @@
if (bdy3d) call do_bdy_3d(0,T)
#endif
call coordinates(vert_cord,cord_relax)
call coordinates(vert_cord,cord_relax,maxdepth)
#ifndef NO_BOTTFRIC
if (kmax .gt. 1) then
call bottom_friction_3d()
......
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment