Commit 2491fa10 authored by kbk's avatar kbk
Browse files

should be fixed and used later

parent ea351e1f
!$Id: coordinates.F90.devel,v 1.1 2003-04-07 17:07:06 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: coordinates() - defines the vertical coordinate system.
!
! !INTERFACE:
subroutine coordinates(cord_type,cord_relax)
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,kmax
use domain, only: H,HU,HV,az,au,av,min_depth
use domain, only: ga,ddu,ddl,d_gamma
use variables_3d, only: dt,kmin,kumin,kvmin,ho,hn,huo,hun,hvo,hvn
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn
use halo_zones, only: update_3d_halo,wait_halo,H_TAG,HU_TAG,HV_TAG
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: cord_type
REALTYPE, intent(in) :: cord_relax
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: coordinates.F90.devel,v $
! Revision 1.1 2003-04-07 17:07:06 kbk
! should be fixed and used later
!
!
! !LOCAL VARIABLES:
integer :: i,j,k,rc,kk
REALTYPE :: alpha,maxdepth
LOGICAL :: gamma_surf=.true.
REALTYPE :: HH,zz,r
logical, save :: first=.true.,equiv_sigma=.false.
REALTYPE, save :: kmaxm1
REALTYPE, save, dimension(:), allocatable :: dga,be,sig
REALTYPE, save, dimension(:,:,:), allocatable :: gga
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'coordinates() # ',Ncall
#endif
if (first) then
first = .false.
allocate(ga(0:kmax),stat=rc)
if (rc /= 0) stop 'coordinates: Error allocating (ga)'
allocate(dga(0:kmax),stat=rc)
if (rc /= 0) STOP 'coordinates: Error allocating (dga)'
ga = _ZERO_; dga = _ZERO_
select case (cord_type)
case (1) ! sigma coordinates
if (ddu .le. _ZERO_ .and. ddl .le. _ZERO_) then
equiv_sigma=.true.
kmaxm1= _ONE_/float(kmax)
ga(0) = -_ONE_
do k=1,kmax
ga(k) = ga(k-1) + _ONE_/kmax
dga(k)=ga(k)-ga(k-1)
end do
ga(kmax) = _ZERO_
else
! Non-equidistant sigma coordinates
! This zooming routine is from Antoine Garapon, ICCH, DK
ga(0)= -_ONE_
dga(0)= _ZERO_
do k=1,kmax
ga(k)=tanh((ddl+ddu)*k/float(kmax)-ddl)+tanh(ddl)
ga(k)=ga(k)/(tanh(ddl)+tanh(ddu)) - _ONE_
dga(k)=ga(k)-ga(k-1)
end do
end if
kmin=1
kumin=1
kvmin=1
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hn(i,j,1:kmax)=max(sseo(i,j)+H(i,j),min_depth)*dga(1:kmax)
hun(i,j,1:kmax)=max(ssuo(i,j)+HU(i,j),min_depth)*dga(1:kmax)
hvn(i,j,1:kmax)=max(ssvo(i,j)+HV(i,j),min_depth)*dga(1:kmax)
end do
end do
case (2) ! z-level
STDERR 'z-level not coded yet'
stop 'coordinates'
case (3) ! general vertical coordinates
do k=0,kmax
ga(k) = k
end do
allocate(sig(0:kmax),stat=rc) ! dimensionless sigma-coordinate
if (rc /= 0) STOP 'coordinates: Error allocating (sig)'
allocate(be(0:kmax),stat=rc) ! dimensionless beta-coordinate
if (rc /= 0) STOP 'coordinates: Error allocating (be)'
allocate(gga(I3DFIELD),stat=rc) ! dimensionless gamma-coordinate
if (rc /= 0) stop 'coordinates: Error allocating memory (gga)'
sig = _ZERO_ ; be = _ZERO_ ; gga = _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_
sig(k)=k/float(kmax)-_ONE_
end do
if (gamma_surf) then
kk=kmax
else
kk=1
end if
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
HH=max(H(i,j),min_depth)
alpha=min( &
((be(kk)-be(kk-1))-D_gamma/HH &
*(sig(kk)-sig(kk-1))) &
/((be(kk)-be(kk-1))-(sig(kk)-sig(kk-1))),_ONE_)
gga(i,j,0)=-_ONE_
do k=1,kmax
gga(i,j,k)=alpha*sig(k)+(1.-alpha)*be(k)
if (gga(i,j,k) .lt. gga(i,j,k-1)) then
STDERR kk,(be(kk)-be(kk-1)),(sig(kk)-sig(kk-1))
STDERR D_gamma,HH
STDERR alpha
STDERR k-1,gga(i,j,k-1),be(k-1),sig(k-1)
STDERR k,gga(i,j,k),be(k),sig(k)
stop 'coordinates'
end if
end do
end do
end do
! call update_3d_halo(gga,gga,az,iimin,jjmin,iimax,jjmax,kmax,H_TAG)
kmin=1
kumin=1
kvmin=1
! call wait_halo(H_TAG)
! Initial layer distribution - uses new as temporary var
do k=1,kmax
do j=jjmin,jjmax
do i=iimin-1,iimax
HH=max(ssuo(i,j)+HU(i,j),min_depth)
hun(i,j,k)=HH*0.5* &
(gga(i,j,k)-gga(i,j,k-1)+gga(i+1,j,k)-gga(i+1,j,k-1))
end do
end do
end do
call update_3d_halo(hun,hun,au,iimin,jjmin,iimax,jjmax,kmax,HU_TAG)
do k=1,kmax
do j=jjmin-1,jjmax
do i=iimin,iimax
HH=max(ssvo(i,j)+HV(i,j),min_depth)
hvn(i,j,k)=HH*0.5* &
(gga(i,j,k)-gga(i,j,k-1)+gga(i,j+1,k)-gga(i,j+1,k-1))
end do
end do
end do
call wait_halo(HU_TAG)
call update_3d_halo(hvn,hvn,av,iimin,jjmin,iimax,jjmax,kmax,HV_TAG)
do k=1,kmax
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
HH=max(sseo(i,j)+H(i,j),min_depth)
hn(i,j,k)=HH*(gga(i,j,k)-gga(i,j,k-1))
end do
end do
end do
call wait_halo(HV_TAG)
case default
end select
end if ! first
ho = hn
huo = hun
hvo = hvn
#if 0
STDERR H(iimin,jjmin),H(iimax,jjmin),H(iimax,jjmax),H(iimin,jjmax)
STDERR H(iimin-1,jjmin-1),H(iimax+1,jjmin-1),H(iimax+1,jjmax+1),H(iimin-1,jjmax+1)
STDERR hn(iimin-1,jjmin-1,1:kmax:2)
STDERR hn(iimax+1,jjmin-1,1:kmax:2)
STDERR hn(iimax+1,jjmax+1,1:kmax:2)
STDERR hn(iimin-1,jjmax+1,1:kmax:2)
#endif
select case (cord_type)
case (1) ! sigma coordinates
if (equiv_sigma) then
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hn(i,j,:)=max(ssen(i,j)+H(i,j),min_depth)*kmaxm1
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hun(i,j,:)=max(ssun(i,j)+HU(i,j),min_depth)*kmaxm1
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hvn(i,j,:)=max(ssvn(i,j)+HV(i,j),min_depth)*kmaxm1
end do
end do
else ! non-equivdistant
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hn(i,j,1:kmax)=max(ssen(i,j)+H(i,j),min_depth)*dga(1:kmax)
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hun(i,j,1:kmax)=max(ssun(i,j)+HU(i,j),min_depth)*dga(1:kmax)
end do
end do
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
hvn(i,j,1:kmax)=max(ssvn(i,j)+HV(i,j),min_depth)*dga(1:kmax)
end do
end do
end if
case (2) ! z-level
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.
! hun and hvn needs to communicate with neighbors therefore they are
! done befoe hn
maxdepth=600. ! needs to be calculated later ... HB
do j=jjmin,jjmax
do i=iimin-1,iimax
if (au(i,j).gt.0) then
r = cord_relax/dt*HU(i,j)/maxdepth
zz = -HU(i,j)
do k=1,kmax-1
HH=max(ssun(i,j)+HU(i,j),min_depth)
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)
end do
hun(i,j,kmax)=ssun(i,j)-zz
end if
end do
end do
call update_3d_halo(hun,hun,au,iimin,jjmin,iimax,jjmax,kmax,HU_TAG)
do j=jjmin-1,jjmax
do i=iimin,iimax
if (av(i,j).gt.0) then
r = cord_relax/dt*HV(i,j)/maxdepth
zz = -HV(i,j)
do k=1,kmax-1
HH=max(ssvn(i,j)+HV(i,j),min_depth)
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
hvn(i,j,kmax)=ssvn(i,j)-zz
end if
end do
end do
call wait_halo(HU_TAG)
call update_3d_halo(hvn,hvn,av,iimin,jjmin,iimax,jjmax,kmax,HV_TAG)
do j=jjmin-1,jjmax+1
do i=iimin-1,iimax+1
if(az(i,j) .gt. 0) then
r = cord_relax/dt*H(i,j)/maxdepth
zz = -H(i,j)
do k=1,kmax-1
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
hn(i,j,kmax)=ssen(i,j)-zz
end if
end do
end do
call wait_halo(HV_TAG)
case default
end select
#ifdef DEBUG
write(debug,*) 'Leaving Coordinates()'
write(debug,*)
#endif
return
end subroutine coordinates
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: start_macro.F90.devel,v 1.1 2003-04-07 17:07:06 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: start_macro() - initialize the macro loop.
!
! !INTERFACE:
subroutine start_macro()
!
! !DESCRIPTION:
!
! !USES:
use domain, only: iimin,iimax,jjmin,jjmax,H,HU,HV,min_depth
use variables_2d, only: z,Uint,Vint
use m3d, only: M,dt
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn
use halo_zones, only: update_2d_halo,z_TAG,U_TAG,V_TAG
IMPLICIT NONE
!
! !INPUT PARAMETERS:
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: start_macro.F90.devel,v $
! Revision 1.1 2003-04-07 17:07:06 kbk
! should be fixed and used later
!
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: split
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'start_macro() # ',Ncall
#endif
do j=jjmin-1,jjmax+1 ! Defining 'old' and 'new' sea surface
do i=iimin-1,iimax+1 ! elevation for macro time step
sseo(i,j)=ssen(i,j)
ssen(i,j)=z(i,j)
#if 1
! This does not work for Sylt - pointed out by Manuel - why?
ssen(i,j)=max(ssen(i,j),-H(i,j)+min_depth)
#endif
end do
end do
do j=jjmin-1,jjmax+1 ! Same for U-points
do i=iimin-1,iimax
ssuo(i,j)=ssun(i,j)
ssun(i,j)=0.25*(sseo(i,j)+sseo(i+1,j)+ssen(i,j)+ssen(i+1,j))
ssun(i,j)=max(ssun(i,j),-HU(i,j)+min_depth)
end do
end do
ssuo(iimax+1,:)=ssun(iimax+1,:)
ssun(iimax+1,:)=ssun(iimax,:)
do j=jjmin-1,jjmax ! Same for V-points
do i=iimin-1,iimax+1
ssvo(i,j)=ssvn(i,j)
ssvn(i,j)=0.25*(sseo(i,j)+sseo(i,j+1)+ssen(i,j)+ssen(i,j+1))
ssvn(i,j)=max(ssvn(i,j),-HV(i,j)+min_depth)
end do
end do
ssvo(:,jjmax+1)=ssvn(:,jjmax+1)
ssvn(:,jjmax+1)=ssvn(:,jjmax)
! Defining vertically integrated, conservative
! u- and v-transport for macro time step
split = _ONE_/float(M)
Uint = split*Uint
Vint = split*Vint
#ifdef DEBUG
write(debug,*) 'Leaving start_macro()'
write(debug,*)
#endif
return
end subroutine start_macro
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
!$Id: ncdf_3d_bdy.F90.devel,v 1.1 2003-04-07 17:09:47 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: ncdf_3d_bdy - input in NetCDF format
!
! !INTERFACE:
module ncdf_3d_bdy
!
! !DESCRIPTION:
!
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax
use domain, only: nsbv,NWB,NNB,NEB,NSB
use domain, only: wi,wfj,wlj,nj,nfi,nli,ei,efj,elj,sj,sfi,sli
use domain, only: H
use m2d, only: dtm
use variables_3d, only: hn
use bdy_3d, only: T_bdy,S_bdy
use time, only: string_to_julsecs,TimeDiff,julianday,secondsofday
IMPLICIT NONE
!
private
!
public :: init_3d_bdy_ncdf,do_3d_bdy_ncdf
!
! !PRIVATE DATA MEMBERS:
integer :: ncid
integer :: time_id,temp_id,salt_id
integer :: start(4),edges(4)
integer :: z_dim,time_dim,time_len,bdy_len
logical :: climatology=.false.,from_3d_fields=.false.
REALTYPE :: offset
REAL_4B, allocatable :: bdy_times(:),wrk(:)
REALTYPE, allocatable, dimension(:,:) :: T_old, T_new
REALTYPE, allocatable, dimension(:,:) :: S_old, S_new
REALTYPE, allocatable, dimension(:,:,:) :: T_bdy_clim,S_bdy_clim
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
!
! $Log: ncdf_3d_bdy.F90.devel,v $
! Revision 1.1 2003-04-07 17:09:47 kbk
! should be fixed and used later
!
! Revision 1.2 2003/03/17 15:42:51 gotm
! Support for special 3D boundary data file
!
! Revision 1.1.1.1 2002/05/02 14:01:49 gotm
! recovering after CVS crash
!
! Revision 1.1 2001/10/17 13:28:27 bbh
! Initial import
!
!EOP
!-----------------------------------------------------------------------
contains
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: init_3d_bdy_ncdf -
!
! !INTERFACE:
subroutine init_3d_bdy_ncdf(fname)
!
! !DESCRIPTION:
! kurt,kurt
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
character(len=*), intent(in) :: fname
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
!
! !REVISION HISTORY:
! Original author(s): Karsten Bolding & Hans Burchard
!
! See log for module
!
! !LOCAL VARIABLES:
character(len=256) :: units
integer :: j1,s1
integer :: ndims
integer, allocatable, dimension(:) :: dim_ids,dim_len
character(len=16), allocatable :: dim_name(:)
REAL_4B, allocatable, dimension(:,:,:) :: sdum,tdum
REAL_4B, allocatable, dimension(:) :: zlev
integer :: rc,err
integer :: i,j,k,l,n,id
!EOP
!-------------------------------------------------------------------------
!BOC
include "netcdf.inc"
#ifdef DEBUG
write(debug,*) 'ncdf_init_3d_bdy (NetCDF)'
write(debug,*) 'Reading from: ',trim(fname)
#endif
LEVEL3 'init_3d_bdy_ncdf'
err = nf_open(fname,NCNOWRIT,ncid)
if (err .NE. NF_NOERR) go to 10
err = nf_inq_ndims(ncid,ndims)
if (err .NE. NF_NOERR) go to 10
allocate(dim_ids(ndims),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (dim_ids)'
allocate(dim_len(ndims),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (dim_len)'
allocate(dim_name(ndims),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (dim_name)'
do n=1,ndims
err = nf_inq_dim(ncid, n, dim_name(n), dim_len(n))
if (err .NE. NF_NOERR) go to 10
STDERR n,dim_name(n), dim_len(n)
end do
if(ndims .eq. 4) then
! We are reading boundary values from a full 3D field
! We assume COARDS conventions
! 1 -> lon,x-axis
! 2 -> lat,y-axis
! 3 -> zax,levels
! 4 -> time
LEVEL4 'boundary data from 3D fields'
from_3d_fields=.true.
z_dim = 3
time_dim = 4
allocate(tdum(imin:imax,jmin:jmax,dim_len(3)),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (tdum)'
allocate(sdum(imin:imax,jmin:jmax,dim_len(3)),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (sdum)'
else
LEVEL4 'special boundary data file'
! We are reading boundary values from a special boundary data file
! 1 -> zax,levels
! 2 -> bdy_points
! 3 -> time
z_dim = 1
time_dim = 3
bdy_len = dim_len(2)
end if
STDERR 'z_dim ',z_dim,' dim_len(z_dim) ',dim_len(z_dim)
allocate(zlev(dim_len(z_dim)),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (zlev)'
err = nf_inq_varid(ncid, dim_name(z_dim), id)
if (err .ne. NF_NOERR) go to 10
err = nf_get_var_real(ncid,id,zlev)
if (err .ne. NF_NOERR) go to 10
zlev = -_ONE_*zlev
time_len = dim_len(time_dim)
if( time_len .le. 12) then
climatology=.true.
LEVEL4 'Assuming climatolgical 3D boundary conditions'
LEVEL4 '# of times = ',time_len
end if
err = nf_inq_varid(ncid,'temp',temp_id)
if (err .NE. NF_NOERR) go to 10
err = nf_inq_varid(ncid,'salt',salt_id)
if (err .NE. NF_NOERR) go to 10
allocate(S_new(nsbv,0:kmax),stat=err)
if (err /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (S_new)'
allocate(T_new(nsbv,0:kmax),stat=err)
if (err /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (T_new)'
if (climatology) then
allocate(T_bdy_clim(time_len,nsbv,0:kmax),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (T_bdy_clim)'
allocate(S_bdy_clim(time_len,nsbv,0:kmax),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (S_bdy_clim)'
allocate(wrk(dim_len(z_dim)),stat=rc)
if (rc /= 0) stop 'init_3d_bdy_ncdf: Error allocating memory (wrk)'
if (from_3d_fields) then
start(1) = imin; edges(1) = imax-imin+1;
start(2) = jmin; edges(2) = jmax-jmin+1;
start(3) = 1; edges(3) = dim_len(3);
edges(4) = 1
do l=1,time_len
start(4) = l
err = nf_get_vara_real(ncid,temp_id,start,edges,tdum)
if (err .ne. NF_NOERR) go to 10
err = nf_get_vara_real(ncid,salt_id,start,edges,sdum)
if (err .ne. NF_NOERR) go to 10
k = 0
do n=1,NWB
i = wi(n)
do j=wfj(1),wlj(1)
k = k+1
wrk(:) = sdum(i,j,:)
call interpol(zlev,wrk,H(i,j),kmax,hn(i,j,:),S_bdy_clim(l,k,:))
wrk(:) = tdum(i,j,:)
call interpol(zlev,wrk,H(i,j),kmax,hn(i,j,:),T_bdy_clim(l,k,:))
end do
end do
do n = 1,NNB
j = nj(n)
do i = nfi(n),nli(n)
k = k+1
wrk(:) = sdum(i,j,:)
call interpol(zlev,wrk,H(i,j),kmax,hn(i,j,:),S_bdy_clim(l,k,:))
wrk(:) = tdum(i,j,:)
call interpol(zlev,wrk,H(i,j),kmax,hn(i,j,:),T_bdy_clim(l,k,:))
end do
end do