Commit 310937b7 authored by Knut's avatar Knut

[U|V]into=>[U|V]adv + incorporate slow_terms into stop_macro

parent ffa3cfc1
......@@ -89,7 +89,6 @@
<File RelativePath="..\..\src\3d\salinity.F90"/>
<File RelativePath="..\..\src\3d\shear_frequency.F90"/>
<File RelativePath="..\..\src\3d\sigma_coordinates.F90"/>
<File RelativePath="..\..\src\3d\slow_terms.F90"/>
<File RelativePath="..\..\src\3d\start_macro.F90"/>
<File RelativePath="..\..\src\3d\stop_macro.F90"/>
<File RelativePath="..\..\src\3d\stresses_3d.F90"/>
......
......@@ -79,7 +79,6 @@
<File RelativePath="..\..\src\3d\salinity.F90"/>
<File RelativePath="..\..\src\3d\shear_frequency.F90"/>
<File RelativePath="..\..\src\3d\sigma_coordinates.F90"/>
<File RelativePath="..\..\src\3d\slow_terms.F90"/>
<File RelativePath="..\..\src\3d\start_macro.F90"/>
<File RelativePath="..\..\src\3d\stop_macro.F90"/>
<File RelativePath="..\..\src\3d\stresses_3d.F90"/>
......
......@@ -13,7 +13,7 @@ TEXSRC = m3d.F90 variables_3d.F90 vertical_coordinates.F90 \
temperature.F90 salinity.F90 eqstate.F90 \
internal_pressure.F90 ip_blumberg_mellor.F90 ip_blumberg_mellor_lin.F90 ip_z_interpol.F90 ip_song_wright.F90 ip_chu_fan.F90 ip_shchepetkin_mcwilliams.F90 ip_stelling_vankester.F90 \
bdy_3d.F90 rivers.F90 spm.F90 \
start_macro.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 slow_terms.F90 stop_macro.F90 \
start_macro.F90 uu_momentum_3d.F90 vv_momentum_3d.F90 ww_momentum_3d.F90 uv_advect_3d.F90 uv_diffusion_3d.F90 stop_macro.F90 \
shear_frequency.F90 buoyancy_frequency.F90 stresses_3d.F90 gotm.F90 tke_eps_advect_3d.F90 physical_dissipation_3d.F90 numerical_mixing.F90 physical_mixing.F90 structure_friction_3d.F90
MOD = \
......@@ -77,7 +77,6 @@ $(LIB)(physical_mixing.o)
endif
OBJ += \
$(LIB)(eddyviscosity.o) \
$(LIB)(slow_terms.o) \
$(LIB)(adv_split_w.o) \
$(LIB)(stresses_3d.o) \
$(LIB)(shear_frequency.o) \
......
......@@ -177,6 +177,12 @@
allocate(Dvn(I2DFIELD),stat=rc) ! depth after macro time step (v-column)
if (rc /= 0) stop 'init_3d: Error allocating memory (Dvn)'
allocate(Uadv(I2DFIELD),stat=rc)
if (rc /= 0) stop 'init_3d: Error allocating memory (Uadv)'
allocate(Vadv(E2DFIELD),stat=rc)
if (rc /= 0) stop 'init_3d: Error allocating memory (Vadv)'
allocate(rru(I2DFIELD),stat=rc) ! Bottom drag term in u-vel. points (3D)
if (rc /= 0) stop 'init_3d: Error allocating memory (rru)'
......
......@@ -52,6 +52,7 @@
REALTYPE, dimension(:,:), allocatable :: ssuo,ssun
REALTYPE, dimension(:,:), allocatable :: ssvo,ssvn
REALTYPE,dimension(:,:),allocatable,target :: t_Dn,t_Dold,Dun,Dvn
REALTYPE,dimension(:,:),allocatable :: Uadv,Vadv
! 3D friction in 3D domain
REALTYPE, dimension(:,:), allocatable :: rru,rrv,zub,zvb
......
......@@ -359,6 +359,7 @@
do i=imin-HALO,imax+HALO
if (au(i,j) .eq. 0) then
uu(i,j,:) = _ZERO_
Uadv(i,j) = _ZERO_
end if
end do
end do
......@@ -366,6 +367,7 @@
do i=imin-HALO,imax+HALO
if (av(i,j) .eq. 0) then
vv(i,j,:) = _ZERO_
Vadv(i,j) = _ZERO_
end if
end do
end do
......@@ -587,12 +589,8 @@
#endif
#ifndef NO_BAROTROPIC
call slow_terms()
#endif
call tic(TIM_INTEGR3D)
call stop_macro()
call toc(TIM_INTEGR3D)
#endif
#ifdef DEBUG
write(debug,*) 'Leaving integrate_3d()'
......
#include "cppdefs.h"
!-----------------------------------------------------------------------
!BOP
!
! !ROUTINE: slow_terms - calculation of slow terms \label{sec-slow-terms}
!
! !INTERFACE:
subroutine slow_terms
!
! !DESCRIPTION:
!
! Here, the calculation of the so-called slow terms (which are the
! interaction terms between the barotropic and the baroclinic mode) is
! completed. The mathematical form of these slow terms is given by
! equations (\ref{Slowfirst}) - (\ref{Slowlast}), see section
! \ref{SectionVerticalIntegrated}.
! These calculations have been prepared in the routines
! {\tt slow\_bottom\_friction}, {\tt slow\_advection} and
! {\tt slow\_diffusion}.
!
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,au,av
use variables_2d, only: Uint,Vint,UEx,VEx,Slru,Slrv,SlUx,SlVx,ru,rv
use variables_3d, only: kumin,kvmin,uu,vv,hun,hvn,Dn,Dun,Dvn
use variables_3d, only: uuEx,vvEx,rru,rrv
#ifndef NO_BAROCLINIC
use variables_3d, only: idpdx,idpdy
#endif
#ifdef STRUCTURE_FRICTION
use variables_3d, only: sf
#endif
use m3d, only: ip_fac
use m2d, only: uv_advect,uv_diffusion,bottom_friction
use getm_timers, only: tic, toc, TIM_SLOWTERMS
!$ use omp_lib
IMPLICIT NONE
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: vertsum
!EOP
!-----------------------------------------------------------------------
!BOC
#ifdef DEBUG
integer, save :: Ncall = 0
Ncall = Ncall+1
write(debug,*) 'slow_terms() # ',Ncall
#endif
call tic(TIM_SLOWTERMS)
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,vertsum)
if (kmax .gt. 1) then
call bottom_friction(Uint,Vint,Dun,Dvn,ru,rv)
call uv_advect(Uint,Vint,Dun,Dvn)
call uv_diffusion(0,Uint,Vint,Dn,Dun,Dvn) ! Has to be called after uv_advect.
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (au(i,j) .ge. 1) then
vertsum=-UEx(i,j)
do k=kumin(i,j),kmax
#ifdef NO_BAROCLINIC
vertsum=vertsum+uuEx(i,j,k)
#else
vertsum=vertsum+uuEx(i,j,k)-ip_fac*idpdx(i,j,k)
#endif
end do
SlUx(i,j)=vertsum
#ifdef NO_SLR
STDERR 'NO_SLR U'
Slru(i,j)= _ZERO_
#else
k=kumin(i,j)
Slru(i,j) = rru(i,j)*uu(i,j,k)/hun(i,j,k) &
- ru(i,j)*Uint(i,j)/Dun(i,j)
#endif
#ifdef STRUCTURE_FRICTION
do k=kumin(i,j),kmax
Slru(i,j)=Slru(i,j)+uu(i,j,k)*_HALF_*(sf(i,j,k)+sf(i+1,j,k))
end do
#endif
end if
if (av(i,j) .ge. 1) then
vertsum=-VEx(i,j)
do k=kvmin(i,j),kmax
#ifdef NO_BAROCLINIC
vertsum=vertsum+vvEx(i,j,k)
#else
vertsum=vertsum+vvEx(i,j,k)-ip_fac*idpdy(i,j,k)
#endif
end do
SlVx(i,j)=vertsum
#ifdef NO_SLR
STDERR 'NO_SLR V'
Slrv(i,j)= _ZERO_
#else
k=kvmin(i,j)
Slrv(i,j) = rrv(i,j)*vv(i,j,k)/hvn(i,j,k) &
- rv(i,j)*Vint(i,j)/Dvn(i,j)
#endif
#ifdef STRUCTURE_FRICTION
do k=kvmin(i,j),kmax
Slrv(i,j)=Slrv(i,j)+vv(i,j,k)*_HALF_*(sf(i,j,k)+sf(i,j+1,k))
end do
#endif
end if
end do
end do
!$OMP END DO
#ifndef NO_BAROCLINIC
else
!
! Here kmax=1, so the loops degenerate and there is no need
! to test for k .ge. kumin(i,j).
k=1
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (au(i,j) .ge. 1) then
SlUx(i,j)=-ip_fac*idpdx(i,j,k)
end if
if (av(i,j) .ge. 1) then
SlVx(i,j)=-ip_fac*idpdy(i,j,k)
end if
end do
end do
!$OMP END DO
#endif
end if
!$OMP END PARALLEL
call toc(TIM_SLOWTERMS)
#ifdef DEBUG
write(debug,*) 'Leaving slow_terms()'
write(debug,*)
#endif
return
end subroutine slow_terms
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
......@@ -31,7 +31,7 @@
use domain, only: imin,imax,jmin,jmax,H,HU,HV,min_depth
use m2d, only: z,Uint,Vint
use m3d, only: M
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn,Dn,Dold,Dun,Dvn
use variables_3d, only: sseo,ssen,ssuo,ssun,ssvo,ssvn,Dn,Dold,Dun,Dvn,Uadv,Vadv
use getm_timers, only: tic, toc, TIM_STARTMCR
IMPLICIT NONE
!
......@@ -84,9 +84,9 @@
! Defining vertically integrated, conservative
! u- and v-transport for macro time step
split = _ONE_/float(M)
Uint = split*Uint
Vint = split*Vint
split = _ONE_/M
Uadv = split*Uint
Vadv = split*Vint
call toc(TIM_STARTMCR)
#ifdef DEBUG
......
......@@ -82,6 +82,7 @@
REALTYPE :: ssvo(I2DFIELD)
REALTYPE :: ssvn(I2DFIELD)
REALTYPE,dimension(I2DFIELD),target :: t_Dn,t_Dold,Dun,Dvn
REALTYPE,dimension(I2DFIELD) :: Uadv,Vadv
! 3D friction in 3D domain
REALTYPE :: rru(I2DFIELD)
......
......@@ -10,18 +10,36 @@
! !DESCRIPTION:
!
! This routine should be called from {\tt m3d} at the end of each macro
! time step in order to reinitialise the transports {\tt Uint} and
! {\tt Vint} to zero.
! time step in order to calculate the so-called slow terms and to
! reinitialise the transports {\tt Uint} and {\tt Vint} to zero.
! The mathematical form of the interaction terms between the barotropic
! and the baroclinic mode is given by
! equations (\ref{Slowfirst}) - (\ref{Slowlast}), see section
! \ref{SectionVerticalIntegrated}.
!
! !USES:
use variables_2d, only: Uint,Vint
use domain, only: imin,imax,jmin,jmax,kmax,au,av
use variables_2d, only: Uint,Vint,UEx,VEx,Slru,Slrv,SlUx,SlVx,ru,rv
use variables_3d, only: kumin,kvmin,uu,vv,hun,hvn,Dn,Dun,Dvn
use variables_3d, only: Uadv,Vadv,uuEx,vvEx,rru,rrv
#ifndef NO_BAROCLINIC
use variables_3d, only: idpdx,idpdy
#endif
#ifdef STRUCTURE_FRICTION
use variables_3d, only: sf
#endif
use m3d, only: ip_fac
use m2d, only: uv_advect,uv_diffusion,bottom_friction
use getm_timers, only: tic, toc, TIM_STOPMCR
!$ use omp_lib
IMPLICIT NONE
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
integer :: i,j,k
REALTYPE :: vertsum
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -30,10 +48,112 @@
Ncall = Ncall+1
write(debug,*) 'stop_macro() # ',Ncall
#endif
call tic(TIM_STOPMCR)
!$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k,vertsum)
if (kmax .gt. 1) then
call bottom_friction(Uadv,Vadv,Dun,Dvn,ru,rv)
call uv_advect(Uadv,Vadv,Dun,Dvn)
call uv_diffusion(0,Uadv,Vadv,Dn,Dun,Dvn) ! Has to be called after uv_advect.
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (au(i,j) .ge. 1) then
vertsum=-UEx(i,j)
do k=kumin(i,j),kmax
#ifdef NO_BAROCLINIC
vertsum=vertsum+uuEx(i,j,k)
#else
vertsum=vertsum+uuEx(i,j,k)-ip_fac*idpdx(i,j,k)
#endif
end do
SlUx(i,j)=vertsum
#ifdef NO_SLR
STDERR 'NO_SLR U'
Slru(i,j)= _ZERO_
#else
k=kumin(i,j)
Slru(i,j) = rru(i,j)*uu(i,j,k)/hun(i,j,k) &
- ru(i,j)*Uadv(i,j)/Dun(i,j)
#endif
#ifdef STRUCTURE_FRICTION
do k=kumin(i,j),kmax
Slru(i,j)=Slru(i,j)+uu(i,j,k)*_HALF_*(sf(i,j,k)+sf(i+1,j,k))
end do
#endif
end if
Uint= _ZERO_
Vint= _ZERO_
if (av(i,j) .ge. 1) then
vertsum=-VEx(i,j)
do k=kvmin(i,j),kmax
#ifdef NO_BAROCLINIC
vertsum=vertsum+vvEx(i,j,k)
#else
vertsum=vertsum+vvEx(i,j,k)-ip_fac*idpdy(i,j,k)
#endif
end do
SlVx(i,j)=vertsum
#ifdef NO_SLR
STDERR 'NO_SLR V'
Slrv(i,j)= _ZERO_
#else
k=kvmin(i,j)
Slrv(i,j) = rrv(i,j)*vv(i,j,k)/hvn(i,j,k) &
- rv(i,j)*Vadv(i,j)/Dvn(i,j)
#endif
#ifdef STRUCTURE_FRICTION
do k=kvmin(i,j),kmax
Slrv(i,j)=Slrv(i,j)+vv(i,j,k)*_HALF_*(sf(i,j,k)+sf(i,j+1,k))
end do
#endif
end if
end do
end do
!$OMP END DO
#ifndef NO_BAROCLINIC
else
!
! Here kmax=1, so the loops degenerate and there is no need
! to test for k .ge. kumin(i,j).
k=1
!$OMP DO SCHEDULE(RUNTIME)
do j=jmin,jmax
do i=imin,imax
if (au(i,j) .ge. 1) then
SlUx(i,j)=-ip_fac*idpdx(i,j,k)
end if
if (av(i,j) .ge. 1) then
SlVx(i,j)=-ip_fac*idpdy(i,j,k)
end if
end do
end do
!$OMP END DO
#endif
end if
!$OMP END PARALLEL
Uint= _ZERO_
Vint= _ZERO_
call toc(TIM_STOPMCR)
#ifdef DEBUG
write(debug,*) 'Leaving stop_macro()'
write(debug,*)
......@@ -41,7 +161,6 @@
return
end subroutine stop_macro
!EOC
!-----------------------------------------------------------------------
! Copyright (C) 2001 - Hans Burchard and Karsten Bolding !
!-----------------------------------------------------------------------
......@@ -52,10 +52,10 @@
#else
use domain, only: dx,dy
#endif
use variables_2d, only: Uint,D
use bdy_3d, only: do_bdy_3d
use variables_3d, only: dt,cnpar,kumin,uu,vv,huo,hun,hvo,uuEx,ww,hvn
use variables_3d, only: num,nuh,sseo,Dun,rru
use variables_3d, only: Uadv,Dn
#ifdef _MOMENTUM_TERMS_
use variables_3d, only: tdv_u,cor_u,ipg_u,epg_u,vsd_u,hsd_u
#endif
......@@ -190,8 +190,8 @@
! Barotropic pressure gradient
#ifndef NO_BAROTROPIC
zp=max(sseo(i+1,j),-H(i ,j)+min(min_depth,D(i+1,j)))
zm=max(sseo(i ,j),-H(i+1,j)+min(min_depth,D(i ,j)))
zp=max(sseo(i+1,j),-H(i ,j)+min(min_depth,Dn(i+1,j)))
zm=max(sseo(i ,j),-H(i+1,j)+min(min_depth,Dn(i ,j)))
zx=(zp-zm+(airp(i+1,j)-airp(i,j))*gammai)/DXU
#else
zx=_ZERO_
......@@ -236,13 +236,13 @@
call getm_tridiagonal(kmax,kumin(i,j),kmax,a1,a2,a3,a4,Res)
! Transport correction: the integral of the new velocities has to
! be the same than the transport calculated by the external mode, Uint.
! be the same than the transport calculated by the external mode, Uadv.
ResInt= _ZERO_
do k=kumin(i,j),kmax
ResInt=ResInt+Res(k)
end do
Diff=(Uint(i,j)-ResInt)/Dun(i,j)
Diff=(Uadv(i,j)-ResInt)/Dun(i,j)
do k=kumin(i,j),kmax
#ifdef _MOMENTUM_TERMS_
......@@ -296,7 +296,7 @@
#endif
end do
else ! if (kmax .eq. kumin(i,j))
uu(i,j,kmax)=Uint(i,j)
uu(i,j,kmax)=Uadv(i,j)
end if
end if
end do
......
......@@ -204,6 +204,7 @@
#endif
ssen = _ZERO_ ; ssun = _ZERO_ ; ssvn = _ZERO_
Dn = _ZERO_ ; Dold = _ZERO_ ; Dun = _ZERO_ ; Dvn = _ZERO_
Uadv = _ZERO_ ; Vadv = _ZERO_
zub = -9999.0 ; zvb = -9999.0 ! must be initialised for gotm
if (bottfric_method .eq. 1) then
......
......@@ -57,8 +57,8 @@
#else
use domain, only: dx,dy
#endif
use variables_2d, only: Vint,D
use bdy_3d, only: do_bdy_3d
use variables_3d, only: Vadv,Dn
use variables_3d, only: dt,cnpar,kvmin,uu,vv,huo,hvo,hvn,vvEx,ww,hun
use variables_3d, only: num,nuh,sseo,Dvn,rrv
#ifdef _MOMENTUM_TERMS_
......@@ -210,8 +210,8 @@
! Barotropic pressure gradient
#ifndef NO_BAROTROPIC
zp=max(sseo(i,j+1),-H(i,j )+min(min_depth,D(i,j+1)))
zm=max(sseo(i,j ),-H(i,j+1)+min(min_depth,D(i,j )))
zp=max(sseo(i,j+1),-H(i,j )+min(min_depth,Dn(i,j+1)))
zm=max(sseo(i,j ),-H(i,j+1)+min(min_depth,Dn(i,j )))
zy=(zp-zm+(airp(i,j+1)-airp(i,j))/gamma)/DYV
#else
zy=_ZERO_
......@@ -256,13 +256,13 @@
call getm_tridiagonal(kmax,kvmin(i,j),kmax,a1,a2,a3,a4,Res)
! Transport correction: the integral of the new velocities has to
! be the same than the transport calculated by the external mode, Vint.
! be the same than the transport calculated by the external mode, Vadv.
ResInt= _ZERO_
do k=kvmin(i,j),kmax
ResInt=ResInt+Res(k)
end do
Diff=(Vint(i,j)-ResInt)/Dvn(i,j)
Diff=(Vadv(i,j)-ResInt)/Dvn(i,j)
do k=kvmin(i,j),kmax
#ifdef _MOMENTUM_TERMS_
......@@ -317,7 +317,7 @@
#endif
end do
else ! (kmax .eq. kvmin(i,j))
vv(i,j,kmax)=Vint(i,j)
vv(i,j,kmax)=Vadv(i,j)
end if
end if
end do
......
......@@ -50,7 +50,6 @@
integer, parameter :: TIM_SS = 40 ! 3d shear_frequency
integer, parameter :: TIM_NN = 41 ! 3d buoyancy_frequency
integer, parameter :: TIM_STRESSES3D = 44 ! 3d stresses_3d
integer, parameter :: TIM_SLOWTERMS = 46 ! 3d slow_terms
integer, parameter :: TIM_TEMP = 52 ! 3d do_temperature
integer, parameter :: TIM_TEMPH = 53 ! 3d temperature halo (presently in m3d/do_integrate_3d)
integer, parameter :: TIM_SALT = 54 ! 3d do_salinity
......@@ -198,12 +197,12 @@
timernames(TIM_SS) = 'shear_frequency'
timernames(TIM_NN) = 'buoyancy_frequency'
timernames(TIM_STRESSES3D) = 'stresses_3d'
timernames(TIM_SLOWTERMS) = 'slow_terms'
timernames(TIM_TEMP) = 'do_temperature'
timernames(TIM_SALT) = 'do_salinity'
timernames(TIM_COORDS) = 'coordinates'
timernames(TIM_INTPRESS) = 'do_internal_pressure'
timernames(TIM_STARTMCR) = 'start_macro'
timernames(TIM_STOPMCR) = 'start_macro'
timernames(TIM_EQSTATE) = 'eq_state'
timernames(TIM_CALCMEANF) = 'calc_mean_fields'
......
......@@ -24,8 +24,7 @@
#else
use domain, only: dx,dy,ard1
#endif
use variables_2d, only: Uinto,Vinto
use variables_3d, only: ssen,Dn,Dold,Dun,Dvn
use variables_3d, only: ssen,Dn,Dold,Dun,Dvn,Uadv,Vadv
use variables_3d, only: dt,kmin,ho,hn,uu,hun,vv,hvn,ww,hcc,SS
use variables_3d, only: taubx,tauby
#ifdef _MOMENTUM_TERMS_
......@@ -119,13 +118,13 @@
! transports
if (u_adv_id .ne. -1) then
call cnv_2d(imin,jmin,imax,jmax,az,Uinto,vel_missing, &
call cnv_2d(imin,jmin,imax,jmax,az,Uadv,vel_missing, &
imin,jmin,imax,jmax,ws2d)
err = nf90_put_var(ncid,u_adv_id,ws2d(_2D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
end if
if (v_adv_id .ne. -1) then
call cnv_2d(imin,jmin,imax,jmax,az,Vinto,vel_missing, &
call cnv_2d(imin,jmin,imax,jmax,az,Vadv,vel_missing, &
imin,jmin,imax,jmax,ws2d)
err = nf90_put_var(ncid,v_adv_id,ws2d(_2D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
......@@ -141,7 +140,7 @@
#else
dx,dy,ard1, &
#endif
xc,xu,xv,Dn,Dold,Uinto,Dun,Vinto,Dvn,wrk2d,wrk2d,vel_missing,ws2d)
xc,xu,xv,Dn,Dold,Uadv,Dun,Vadv,Dvn,wrk2d,wrk2d,vel_missing,ws2d)
err = nf90_put_var(ncid,velx_adv_id,ws2d(_2D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
end if
......@@ -154,7 +153,7 @@
#else
dx,dy,ard1, &
#endif
yc,yu,yv,Dn,Dold,Uinto,Dun,Vinto,Dvn,wrk2d,wrk2d,vel_missing,ws2d)
yc,yu,yv,Dn,Dold,Uadv,Dun,Vadv,Dvn,wrk2d,wrk2d,vel_missing,ws2d)
err = nf90_put_var(ncid,vely_adv_id,ws2d(_2D_W_),start,edges)
if (err .NE. NF90_NOERR) go to 10
end if
......