Commit b4e340bc authored by kbk's avatar kbk
Browse files

added NetCDF IO rewrite + de-stag of velocities - Umlauf

parent 5b50155f
!$Id: vv_momentum_3d.F90,v 1.6 2004-07-28 14:58:18 hb Exp $
!$Id: vv_momentum_3d.F90,v 1.7 2005-04-25 09:32:34 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
......@@ -12,7 +12,7 @@
!
! Three-dimensional velocity equation in northern direction.
! If #MUDFLAT is defined, fitting of profiles is made with
! respect to the new surface elevation, otherwise to the
! respect to the new surface elevation, otherwise to the
! old surface elevation.
!
! !USES:
......@@ -50,7 +50,10 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! $Log: vv_momentum_3d.F90,v $
! Revision 1.6 2004-07-28 14:58:18 hb
! Revision 1.7 2005-04-25 09:32:34 kbk
! added NetCDF IO rewrite + de-stag of velocities - Umlauf
!
! Revision 1.6 2004/07/28 14:58:18 hb
! Changing subroutine calling order via MUDFLAT
!
! Revision 1.5 2004/04/20 16:49:37 hb
......
#$Id: Makefile,v 1.2 2003-04-07 16:45:07 kbk Exp $
#$Id: Makefile,v 1.3 2005-04-25 09:32:34 kbk Exp $
#
# Makefile to build the domain library - libdomain.a
#
......@@ -9,14 +9,16 @@ INCS = static_domain.h dynamic_declarations_domain.h \
dynamic_allocations_domain.h
LIB = $(LIBDIR)/libdomain${buildtype}.a
MODSRC = domain.F90
MODSRC = topo_interface.F90 domain.F90
LIBSRC = part_domain.F90 bdy_spec.F90 mirror_bdy_2d.F90 print_bdy.F90
LIBSRC = part_domain.F90 bdy_spec.F90 mirror_bdy_2d.F90 print_bdy.F90
SRC = $(MODSRC) $(LIBSRC)
MOD = \
${LIB}(domain.o)
${LIB}(topo_interface.o) \
${LIB}(domain.o)
OBJ = \
${LIB}(part_domain.o) \
......@@ -39,7 +41,7 @@ doc:
$(PROTEX) $(SRC) > $(DOCDIR)/domain.tex
clean:
$(RM) $(LIB) $(MODDIR)/domain.{m.mod}
$(RM) $(LIB) $(MODDIR)/domain.{m.mod} $(MODDIR)/topo_interface.{m.mod}
realclean: clean
$(RM) *.o
......
This diff is collapsed.
allocate(H(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (H)'
allocate(lonc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (lonc)'
! mask
allocate(az(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (az)'
allocate(latc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (latc)'
allocate(au(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (au)'
allocate(conv(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (conv)'
allocate(av(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (av)'
allocate(cor(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (cor)'
allocate(ax(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (ax)'
allocate(coru(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (coru)'
allocate(corv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (corv)'
! bathymetry
allocate(H(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (H)'
allocate(HU(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (HU)'
......@@ -39,27 +36,24 @@
if (rc /= 0) stop 'init_domain: Error allocating memory (dry_v)'
dry_v = _ONE_
allocate(az(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (az)'
allocate(au(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (au)'
allocate(av(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (av)'
! coriolis terms
allocate(cor(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (cor)'
allocate(ax(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (ax)'
allocate(coru(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (coru)'
#if ! ( defined(CURVILINEAR) || defined(SPHERICAL) )
allocate(corv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (corv)'
allocate(xc(imin-HALO:imax+HALO),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (xc)'
allocate(yc(jmin-HALO:jmax+HALO),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (yc)'
! lat/lon
allocate(lonc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (lonc)'
#else
allocate(latc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (latc)'
allocate(lonx(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (lonx)'
......@@ -79,15 +73,53 @@
allocate(latv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (latv)'
! grid convergence
allocate(convc(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (convc)'
allocate(convx(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (convx)'
allocate(angle(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (angle)'
! grid points
allocate(xx(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xx)'
allocate(yx(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yx)'
allocate(xc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xc)'
allocate(yc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yc)'
allocate(xu(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xu)'
allocate(yu(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yu)'
allocate(xv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xv)'
allocate(yv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yv)'
! metric parameters
allocate(dxdyc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (dxdyc)'
allocate(dydxc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (dydxc)'
allocate(angle(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (angle)'
allocate(dxc(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (dxc)'
......@@ -121,31 +153,6 @@
allocate(arvd1(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_domain: Error allocating memory (arvd1)'
#if defined(CURVILINEAR)
allocate(xc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xc)'
allocate(yc(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yc)'
allocate(xx(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xx)'
allocate(yx(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yx)'
allocate(xu(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xu)'
allocate(yu(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yu)'
allocate(xv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (xv)'
allocate(yv(E2DFIELD),stat=rc)
if (rc /=0) stop 'init_domain: Error allocating memory (yv)'
#endif
#endif
......@@ -3,26 +3,39 @@
integer :: iimin=-1,iimax=-1,jjmin=-1,jjmax=-1
integer :: kmax=1
! mask
integer, dimension(:,:), allocatable :: az,au,av,ax
! bathymetry
REALTYPE, dimension(:,:), allocatable :: H,HU,HV
REALTYPE, dimension(:,:), allocatable :: lonc,latc,conv
REALTYPE, dimension(:,:), allocatable :: dry_z,dry_u,dry_v
! coriolis terms
REALTYPE, dimension(:,:), allocatable :: cor,coru,corv
#if ! ( defined(CURVILINEAR) || defined(SPHERICAL) )
REALTYPE, dimension(:), allocatable :: xc,yc
REALTYPE :: dx,dy,ard1
#else
! lat/lon
REALTYPE, dimension(:,:), allocatable :: lonc,latc
REALTYPE, dimension(:,:), allocatable :: lonx,latx
REALTYPE, dimension(:,:), allocatable :: lonu,latu
REALTYPE, dimension(:,:), allocatable :: lonv,latv
REALTYPE, dimension(:,:), allocatable :: dxdyc,dydxc,angle
! grid convergence
REALTYPE, dimension(:,:), allocatable :: angle,convc,convx
! grid points
REALTYPE, dimension(:,:), allocatable :: xx,yx
REALTYPE, dimension(:,:), allocatable :: xc,yc
REALTYPE, dimension(:,:), allocatable :: xu,yu
REALTYPE, dimension(:,:), allocatable :: xv,yv
! metric parameters
REALTYPE :: dx,dy,x0,y0,ard1
REALTYPE :: dlon,dlat,lon0,lat0
REALTYPE, dimension(:,:), allocatable :: dxdyc,dydxc
REALTYPE, dimension(:,:), allocatable :: dxc,dxu,dxv,dxx
REALTYPE, dimension(:,:), allocatable :: dyc,dyu,dyv,dyx
REALTYPE, dimension(:,:), allocatable :: arcd1,arud1,arvd1
#if defined(CURVILINEAR)
REALTYPE, dimension(:,:), allocatable:: xc,yc,xu,yu,xv,yv,xx,yx
#endif
#endif
#include "dimensions.h"
! mask
integer :: az(E2DFIELD)
integer :: au(E2DFIELD)
integer :: av(E2DFIELD)
integer :: ax(E2DFIELD)
! bathymetry
REALTYPE :: H(E2DFIELD)
REALTYPE :: HU(E2DFIELD)
REALTYPE :: HV(E2DFIELD)
REALTYPE :: dry_z(E2DFIELD)
REALTYPE :: dry_u(E2DFIELD)
REALTYPE :: dry_v(E2DFIELD)
! coriolis terms
REALTYPE :: cor(E2DFIELD)
REALTYPE :: coru(E2DFIELD)
REALTYPE :: corv(E2DFIELD)
! lat/lon
REALTYPE :: lonc(E2DFIELD)
REALTYPE :: latc(E2DFIELD)
REALTYPE :: conv(E2DFIELD)
#if ( ! defined(SPHERICAL) && ! defined(CURVILINEAR) )
REALTYPE :: xc(imin-HALO:imax+HALO)
REALTYPE :: yc(jmin-HALO:jmax+HALO)
REALTYPE :: zc(0:kmax)
REALTYPE :: dx,dy,ard1
#else
REALTYPE :: lonx(E2DFIELD)
REALTYPE :: latx(E2DFIELD)
REALTYPE :: lonu(E2DFIELD)
REALTYPE :: latu(E2DFIELD)
REALTYPE :: lonv(E2DFIELD)
REALTYPE :: latv(E2DFIELD)
! grid convergence
REALTYPE :: angle(E2DFIELD)
REALTYPE :: convc(E2DFIELD)
REALTYPE :: convx(E2DFIELD)
! grid points
REALTYPE :: xx(E2DFIELD)
REALTYPE :: yx(E2DFIELD)
REALTYPE :: xc(E2DFIELD)
REALTYPE :: yc(E2DFIELD)
REALTYPE :: xu(E2DFIELD)
REALTYPE :: yu(E2DFIELD)
REALTYPE :: xv(E2DFIELD)
REALTYPE :: yv(E2DFIELD)
! metric parameters
REALTYPE :: dx,dy,x0,y0,ard1
REALTYPE :: dlon,dlat,lon0,lat0
REALTYPE :: dxdyc(E2DFIELD)
REALTYPE :: dydxc(E2DFIELD)
REALTYPE :: angle(E2DFIELD)
REALTYPE :: dxc(E2DFIELD)
REALTYPE :: dxu(E2DFIELD)
REALTYPE :: dxv(E2DFIELD)
......@@ -44,14 +60,4 @@
REALTYPE :: arcd1(E2DFIELD)
REALTYPE :: arud1(E2DFIELD)
REALTYPE :: arvd1(E2DFIELD)
#if defined(CURVILINEAR)
REALTYPE :: xx(E2DFIELD)
REALTYPE :: yx(E2DFIELD)
REALTYPE :: xc(E2DFIELD)
REALTYPE :: yc(E2DFIELD)
REALTYPE :: xu(E2DFIELD)
REALTYPE :: yu(E2DFIELD)
REALTYPE :: xv(E2DFIELD)
REALTYPE :: yv(E2DFIELD)
#endif
#endif
!$Id: topo_interface.F90,v 1.1 2005-04-25 09:32:34 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !MODULE: topo_interface --- for type safety of grid-related functions calls
!
! !INTERFACE:
module topo_interface
!
! !DESCRIPTION:
! This module contains explicit interfaces for some external functions related
! to checking and reading the topography and grid. These explicit interfaces
! prevent functions in the module {\tt domain} from calling the routines with
! any incorrect type and number of arguments (type-safety).
!
! !USES:
IMPLICIT NONE
!
! !REVISION HISTORY:
! Original author(s): Lars Umlauf
!
! $Log: topo_interface.F90,v $
! Revision 1.1 2005-04-25 09:32:34 kbk
! added NetCDF IO rewrite + de-stag of velocities - Umlauf
!
!
!EOP
!-----------------------------------------------------------------------
interface
subroutine check_grid(filename,filetype,iextr,jextr)
character(len=*) :: filename
integer, intent(in) :: filetype
#ifdef STATIC
integer, intent(in) :: iextr
integer, intent(in) :: jextr
#else
integer, intent(out) :: iextr
integer, intent(out) :: jextr
#endif
end subroutine check_grid
end interface
interface
subroutine get_grid(filetype,H,Hland, &
iextr,jextr,ioff,joff,imin,imax,jmin,jmax)
integer, intent(in) :: filetype
integer, intent(in) :: iextr,jextr,ioff,joff
integer, intent(in) :: imin,imax,jmin,jmax
REALTYPE, intent(out) :: H(E2DFIELD)
REALTYPE, intent(out) :: Hland
end subroutine get_grid
end interface
!-----------------------------------------------------------------------
end module topo_interface
!-----------------------------------------------------------------------
! Copyright (C) 2005 - Lars Umlauf, Hans Burchard and Karsten Bolding
!-----------------------------------------------------------------------
#$Id: Makefile,v 1.5 2005-04-25 07:55:49 kbk Exp $
#$Id: Makefile,v 1.6 2005-04-25 09:32:34 kbk Exp $
#
# Makefile to build utilities written in Fortran90 - libfutils.a
#
......@@ -11,7 +11,9 @@ MODSRC = exceptions.F90 parallel.F90 parameters.F90 time.F90 grid_interpol.F90
LIBSRC = ver_interpol.F90 kbk_interpol.F90 tridiagonal.F90 pos.F90 \
cnv_2d.F90 cnv_3d.F90 eta_mask.F90 col_interpol.F90 \
to_2d_vel.F90 to_3d_vel.F90 tow.F90
to_2d_vel.F90 to_2d_u.F90 to_2d_v.F90 \
to_3d_vel.F90 to_3d_uu.F90 to_3d_vv.F90 \
tow.F90 c2x.F90
SRC = $(MODSRC) $(LIBSRC)
......@@ -39,8 +41,13 @@ ${LIB}(cnv_3d.o) \
${LIB}(eta_mask.o) \
${LIB}(col_interpol.o) \
${LIB}(to_2d_vel.o) \
${LIB}(to_2d_u.o) \
${LIB}(to_2d_v.o) \
${LIB}(to_3d_vel.o) \
${LIB}(tow.o)
${LIB}(to_3d_uu.o) \
${LIB}(to_3d_vv.o) \
${LIB}(tow.o) \
${LIB}(c2x.o)
all: modules objects
......
!$Id: c2x.F90,v 1.1 2005-04-25 09:32:34 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: Interpolate fields from central T-points to X-points
!
! !INTERFACE:
subroutine c2x(iimin,iimax,jjmin,jjmax,cfield,xfield)
!
! !DESCRIPTION:
! This routine interpolates a variable given on the
! T-points to the X-points. At the edges and corners,
! an extrapolation is performed, since the outermost
! points are X-points in the GETM grid layout.
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: iimin,jjmin,iimax,jjmax
REALTYPE, intent(in) :: cfield(I2DFIELD)
!
! !INPUT/OUTPUT PARAMETERS:
!
! !OUTPUT PARAMETERS:
REALTYPE, intent(out) :: xfield(I2DFIELD)
!
! !REVISION HISTORY:
! Original author(s): Lars Umlauf
!
! $Log: c2x.F90,v $
! Revision 1.1 2005-04-25 09:32:34 kbk
! added NetCDF IO rewrite + de-stag of velocities - Umlauf
!
!
! !LOCAL VARIABLES:
integer :: i,j
REALTYPE :: ufield(I2DFIELD),vfield(I2DFIELD)
!EOP
!-----------------------------------------------------------------------
!BOC
! do the interior X-points
do j=jjmin,jjmax-1
do i=iimin,iimax-1
xfield(i,j) = 0.25*( cfield(i ,j) + cfield(i+1,j+1) &
+ cfield(i+1,j) + cfield(i ,j+1) )
end do
end do
! do the interior U and V-points as an intermediate step
do j=jjmin,jjmax
do i=iimin,iimax-1
ufield(i,j) = 0.5*( cfield(i,j) + cfield(i+1,j) )
end do
end do
do j=jjmin,jjmax-1
do i=iimin,iimax
vfield(i,j) = 0.5*( cfield(i,j) + cfield(i,j+1) )
end do
end do
! do the edges
do i=iimin,iimax-1
xfield(i,jjmin-1) = 2.0*ufield(i,jjmin) - xfield(i,jjmin )
xfield(i,jjmax ) = 2.0*ufield(i,jjmax) - xfield(i,jjmax-1)
end do
do i=jjmin,jjmax-1
xfield(iimin-1,j) = 2.0*vfield(iimin,j) - xfield(iimin,j )
xfield(iimax ,j) = 2.0*vfield(iimax,j) - xfield(iimax-1,j)
end do
! do the exterior corners
xfield(iimin-1,jjmin-1) = 2.0*ufield(iimin-1,jjmin) - xfield(iimin-1,jjmin )
xfield(iimin-1,jjmax ) = 2.0*ufield(iimin-1,jjmax) - xfield(iimin-1,jjmax-1)
xfield(iimax ,jjmin-1) = 2.0*ufield(iimax ,jjmin) - xfield(iimax ,jjmin )
xfield(iimax ,jjmax ) = 2.0*ufield(iimax ,jjmax) - xfield(iimax ,jjmax-1)
return
end subroutine c2x
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2005 - Lars Umlauf, Hans Burchard and Karsten Bolding
!-----------------------------------------------------------------------
!$Id: to_2d_u.F90,v 1.1 2005-04-25 09:32:34 kbk Exp $
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: Mask U-velocity and interpolate to T-points
!
! !INTERFACE:
subroutine to_2d_u(imin,jmin,imax,jmax,az,u,DU,missing, &
il,jl,ih,jh,vel)
!
! !DESCRIPTION:
! This routine interpolates the $U$-velocity component from the $U$-points
! to the $T$-points. If the mask for the $T$-points is zero, the positions
! are filled up with "missing values". The result is written in a single-precision
! vector four netCDF output.
!
! !USES:
IMPLICIT NONE
!
! !INPUT PARAMETERS:
integer, intent(in) :: imin,jmin,imax,jmax
integer, intent(in) :: az(E2DFIELD)
REALTYPE, intent(in) :: u(E2DFIELD)
REALTYPE, intent(in) :: DU(E2DFIELD)
REALTYPE, intent(in) :: missing
integer, intent(in) :: il,jl,ih,jh
! !OUTPUT PARAMETERS:
REAL_4B, intent(out) :: vel(*)
!
! !REVISION HISTORY:
! Original author(s): Lars Umlauf
!
! $Log: to_2d_u.F90,v $
! Revision 1.1 2005-04-25 09:32:34 kbk
! added NetCDF IO rewrite + de-stag of velocities - Umlauf
!
!
! !LOCAL VARIABLES:
integer :: i,j,indx
!EOP
!-----------------------------------------------------------------------
!BOC
indx = 1
do j=jl,jh
do i=il,ih
if (az(i,j) .gt. 0) then