Skip to content

GitLab

G
getm
Commit c2406931 authored by Knut's avatar Knut
Browse files
Options

new numdis output

parent 41d8170b
Hide whitespace changes
Inline Side-by-side
Showing with 296 additions and 181 deletions
src/2d/m2d.F90
View file @ c2406931
......@@ -43,11 +43,12 @@
REALTYPE,dimension(E2DFIELD),intent(out) :: DU,DV
end subroutine depth_update
subroutine uv_advect(U,V,Dold,Dnew,DU,DV)
subroutine uv_advect(U,V,Dold,Dnew,DU,DV,numdis)
use domain, only: imin,imax,jmin,jmax
IMPLICIT NONE
REALTYPE,dimension(E2DFIELD),intent(in) :: U,V
REALTYPE,dimension(E2DFIELD),target,intent(in) :: Dold,Dnew,DU,DV
REALTYPE,dimension(:,:),pointer,intent(out),optional :: numdis
end subroutine uv_advect
subroutine uv_diffusion(An_method,U,V,D,DU,DV)
......@@ -371,16 +372,13 @@
allocate(phydis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (phydis_2d)'
phydis_2d = _ZERO_
allocate(numdis_u_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_u_2d)'
numdis_u_2d = _ZERO_
allocate(numdis_v_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_v_2d)'
numdis_v_2d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_2d(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d)'
numdis_2d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_2d_old(E2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_2d: Error allocating memory (numdis_2d_old)'
numdis_2d_old = _ZERO_
#endif
end if
......@@ -492,7 +490,7 @@
end if
call tic(TIM_INTEGR2D)
call uv_advect(U,V,Dlast,D,DU,DV)
call uv_advect(U,V,Dlast,D,DU,DV,numdis_2d)
call uv_diffusion(An_method,U,V,D,DU,DV) ! Has to be called after uv_advect.
if (do_numerical_analyses_2d) call physical_dissipation(U,V,DU,DV,Am,phydis_2d)
call toc(TIM_INTEGR2D)
......
src/2d/physical_dissipation.F90
View file @ c2406931
......@@ -136,7 +136,7 @@
end do
do j=jmin-HALO,jmax+HALO-1
#endif
do i=imax+HALO-1,imin-HALO+1 ! loop order MUST NOT be changed!!!
do i=imax+HALO-1,imin-HALO+1,-1 ! loop order MUST NOT be changed!!!
! Note (KK): slip condition dudyV(av=0)=0
! prolonged outflow condition dvdxV(av=3)=0
! shearV(av=3) would require shearX outside open boundary
......@@ -150,7 +150,7 @@
end do
#ifndef SLICE_MODEL
end do
do j=jmax+HALO-1,jmin-HALO+1 ! loop order MUST NOT be changed!!!
do j=jmax+HALO-1,jmin-HALO+1,-1 ! loop order MUST NOT be changed!!!
do i=imin-HALO+1,imax+HALO-1
if (az(i,j) .eq. 1) then
! calculate shearC
......
src/2d/uv_advect.F90
View file @ c2406931
......@@ -5,7 +5,7 @@
! !ROUTINE: uv_advect - 2D advection of momentum \label{sec-uv-advect}
!
! !INTERFACE:
subroutine uv_advect(U,V,Dold,Dnew,DU,DV)
subroutine uv_advect(U,V,Dold,Dnew,DU,DV,numdis)
! Note (KK): keep in sync with interface in m2d.F90
!
......@@ -18,14 +18,14 @@
! !USES:
use domain, only: imin,imax,jmin,jmax,az,au,av,ax
#if defined(SPHERICAL) || defined(CURVILINEAR)
use domain, only: dxv,dyu
use domain, only: dxv,dyu,arcd1,arud1,arvd1
#else
use domain, only: dx,dy
use domain, only: dx,dy,ard1
#endif
use m2d, only: dtm,vel2d_adv_split,vel2d_adv_hor
use variables_2d, only: UEx,VEx
use variables_2d, only: do_numerical_analyses_2d
use variables_2d, only: numdis_u_2d,numdis_v_2d,numdis_2d
use variables_2d, only: numdis_2d_old
use advection, only: NOADV,UPSTREAM,J7,do_advection
use halo_zones, only: update_2d_halo,wait_halo,U_TAG,V_TAG
use getm_timers, only: tic,toc,TIM_UVADV,TIM_UVADVH
......@@ -36,15 +36,21 @@
REALTYPE,dimension(E2DFIELD),intent(in) :: U,V
REALTYPE,dimension(E2DFIELD),target,intent(in) :: Dold,Dnew,DU,DV
!
! !OUTPUT PARAMETERS:
REALTYPE,dimension(:,:),pointer,intent(out),optional :: numdis
!
! !REVISION HISTORY:
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
logical :: calc_numdis
integer :: i,j
REALTYPE,dimension(E2DFIELD) :: fadv,Uadv,Vadv,DUadv,DVadv
REALTYPE,dimension(E2DFIELD),target :: Dadv
REALTYPE,dimension(:,:),pointer :: pDadv
REALTYPE,dimension(E2DFIELD) :: fadv,Uadv,Vadv,DUadv,DVadv,Dires
REALTYPE,dimension(E2DFIELD),target :: Dadv,nvd
REALTYPE,dimension(:,:),pointer :: pDadv,p_nvd
#ifdef _NUMERICAL_ANALYSES_OLD_
REALTYPE,dimension(E2DFIELD) :: work2d
#endif
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -58,6 +64,19 @@
#endif
call tic(TIM_UVADV)
if (present(numdis)) then
calc_numdis = associated(numdis)
else
calc_numdis = .false.
end if
if (calc_numdis) then
p_nvd => nvd
else
p_nvd => null()
end if
if (vel2d_adv_hor .ne. NOADV) then
!$OMP PARALLEL DEFAULT(SHARED) &
......@@ -201,9 +220,43 @@
!$OMP SINGLE
call do_advection(dtm,fadv,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,U_TAG, &
advres=UEx,nvd=numdis_u_2d)
Dires=Dires,advres=UEx,nvd=p_nvd)
!$OMP END SINGLE
if (calc_numdis) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
nvd(i,j) = nvd(i,j)*Dires(i,j)/ARUD1
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
!$OMP SINGLE
call update_2d_halo(nvd,nvd,au,imin,jmin,imax,jmax,U_TAG)
call wait_halo(U_TAG)
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis(i,j) = _HALF_*( nvd(i-1,j) + nvd(i,j) )/Dnew(i,j)*ARCD1
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end if
#ifdef _NUMERICAL_ANALYSES_OLD_
if (do_numerical_analyses_2d) then
......@@ -217,7 +270,7 @@
do j=jmin,jmax ! calculate kinetic energy dissipaion rate for u-velocity
#endif
do i=imin,imax
work2d(i,j) = ( work2d(i,j) - fadv(i,j)**2 ) / dtm
work2d(i,j) = ( work2d(i,j) - fadv(i,j)**2 )/dtm*Dires(i,j)/ARUD1
end do
#ifndef SLICE_MODEL
end do
......@@ -235,7 +288,7 @@
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis_2d(i,j) = _HALF_*( work2d(i-1,j) + work2d(i,j) )
numdis_2d_old(i,j) = _HALF_*( work2d(i-1,j) + work2d(i,j) )/Dnew(i,j)*ARCD1
end if
end do
#ifndef SLICE_MODEL
......@@ -369,9 +422,44 @@
!$OMP SINGLE
call do_advection(dtm,fadv,Uadv,Vadv,DUadv,DVadv,pDadv,pDadv, &
vel2d_adv_split,vel2d_adv_hor,_ZERO_,V_TAG, &
advres=VEx,nvd=numdis_v_2d)
Dires=Dires,advres=VEx,nvd=p_nvd)
!$OMP END SINGLE
if (calc_numdis) then
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
nvd(i,j) = nvd(i,j)*Dires(i,j)/ARVD1
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
!$OMP SINGLE
call update_2d_halo(nvd,nvd,av,imin,jmin,imax,jmax,V_TAG)
call wait_halo(V_TAG)
!$OMP END SINGLE
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis(i,j) = numdis(i,j) &
+_HALF_*( nvd(i,j-1) + nvd(i,j) )/Dnew(i,j)*ARCD1
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end if
#ifdef _NUMERICAL_ANALYSES_OLD_
if (do_numerical_analyses_2d) then
......@@ -385,7 +473,7 @@
do j=jmin,jmax ! calculate kinetic energy dissipaion rate for v-velocity
#endif
do i=imin,imax
work2d(i,j) = ( work2d(i,j) - fadv(i,j)**2 ) / dtm
work2d(i,j) = ( work2d(i,j) - fadv(i,j)**2 )/dtm*Dires(i,j)/ARVD1
end do
#ifndef SLICE_MODEL
end do
......@@ -407,8 +495,8 @@
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis_2d(i,j) = numdis_2d(i,j) &
+_HALF_*( work2d(i,j-1) + work2d(i,j) )
numdis_2d_old(i,j) = numdis_2d_old(i,j) &
+_HALF_*( work2d(i,j-1) + work2d(i,j) )/Dnew(i,j)*ARCD1
end if
end do
#ifndef SLICE_MODEL
......
src/2d/variables_2d.F90
View file @ c2406931
......@@ -30,9 +30,8 @@
#include "dynamic_declarations_2d.h"
#endif
REALTYPE,dimension(:,:),pointer :: numdis_u_2d=>null()
REALTYPE,dimension(:,:),pointer :: numdis_v_2d=>null()
REALTYPE,dimension(:,:),allocatable :: numdis_2d,phydis_2d
REALTYPE,dimension(:,:),pointer :: numdis_2d=>null()
REALTYPE,dimension(:,:),allocatable :: numdis_2d_old,phydis_2d
integer :: size2d_field
integer :: mem2d
......
src/3d/m3d.F90
View file @ c2406931
......@@ -313,16 +313,13 @@
allocate(phydis_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (phydis_int)'
phydis_int = _ZERO_
allocate(numdis_u_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_u_3d)'
numdis_u_3d = _ZERO_
allocate(numdis_v_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_v_3d)'
numdis_v_3d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_3d(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_3d)'
numdis_3d = _ZERO_
#ifdef _NUMERICAL_ANALYSES_OLD_
allocate(numdis_3d_old(I3DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_3d_old)'
numdis_3d_old = _ZERO_
allocate(numdis_int(I2DFIELD),stat=rc)
if (rc /= 0) stop 'postinit_3d: Error allocating memory (numdis_int)'
numdis_int = _ZERO_
......
src/3d/uv_advect_3d.F90
View file @ c2406931
......@@ -21,9 +21,9 @@
! !USES:
use domain, only: imin,imax,jmin,jmax,kmax,az,au,av,ax
#if defined(SPHERICAL) || defined(CURVILINEAR)
use domain, only: dxv,dyu
use domain, only: dxv,dyu,arcd1,arud1,arvd1
#else
use domain, only: dx,dy
use domain, only: dx,dy,ard1
#endif
use m3d, only: vel3d_adv_split,vel3d_adv_hor,vel3d_adv_ver
use variables_3d, only: dt,uu,vv,ww,ho,hn,hun,hvn,uuEx,vvEx
......@@ -31,7 +31,7 @@
use advection_3d, only: do_advection_3d
use halo_zones, only: update_3d_halo,wait_halo,U_TAG,V_TAG
use variables_3d, only: do_numerical_analyses_3d
use variables_3d, only: numdis_u_3d,numdis_v_3d,numdis_3d,numdis_int
use variables_3d, only: numdis_3d,numdis_3d_old,numdis_int
#ifdef _MOMENTUM_TERMS_
use domain, only: dry_u,dry_v
use variables_3d, only: adv_u,adv_v
......@@ -44,11 +44,14 @@
! Original author(s): Hans Burchard & Karsten Bolding
!
! !LOCAL VARIABLES:
logical :: calc_numdis
integer :: i,j,k
REALTYPE,dimension(I3DFIELD) :: fadv3d,uuadv,vvadv,wwadv,huadv,hvadv
REALTYPE,dimension(I3DFIELD),target :: hnadv
REALTYPE,dimension(:,:,:),pointer :: phadv
REALTYPE,dimension(I3DFIELD) :: work3d,hires
REALTYPE,dimension(I3DFIELD) :: fadv3d,uuadv,vvadv,wwadv,huadv,hvadv,hires
REALTYPE,dimension(I3DFIELD),target :: hnadv,nvd
REALTYPE,dimension(:,:,:),pointer :: phadv,p_nvd
#ifdef _NUMERICAL_ANALYSES_OLD_
REALTYPE,dimension(I3DFIELD) :: work3d
#endif
!EOP
!-----------------------------------------------------------------------
!BOC
......@@ -62,6 +65,13 @@
#endif
call tic(TIM_UVADV3D)
calc_numdis = associated(numdis_3d)
if (calc_numdis) then
p_nvd => nvd
else
p_nvd => null()
end if
if (vel3d_adv_hor.ne.NOADV .or. vel3d_adv_ver.ne.NOADV) then
! Note (KK): wwadv(kmax) will be overwritten by ww(kmax) anyway
......@@ -214,32 +224,53 @@
!$OMP SINGLE
call do_advection_3d(dt,fadv3d,uuadv,vvadv,wwadv,huadv,hvadv,phadv,phadv, &
vel3d_adv_split,vel3d_adv_hor,vel3d_adv_ver,_ZERO_,U_TAG, &
advres=uuEx,nvd=numdis_u_3d)
hires=hires,advres=uuEx,nvd=p_nvd)
!$OMP END SINGLE
#ifdef _MOMENTUM_TERMS_
do k=1,kmax
if (calc_numdis) then
do k=1,kmax
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
do j=jmin,jmax
#endif
do i=imin,imax
adv_u(i,j,k) = dry_u(i,j) * uuEx(i,j,k)
end do
do i=imin,imax
nvd(i,j,k) = nvd(i,j,k)*hires(i,j,k)/ARUD1
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end do
!$OMP SINGLE
call update_3d_halo(nvd,nvd,au,imin,jmin,imax,jmax,kmax,U_TAG)
call wait_halo(U_TAG)
!$OMP END SINGLE
do k=1,kmax
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
!$OMP END DO NOWAIT
end do
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis_3d(i,j,k) = _HALF_*( nvd(i-1,j,k) + nvd(i,j,k) )/hn(i,j,k)*ARCD1
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end do
end if
#ifdef _NUMERICAL_ANALYSES_OLD_
if (do_numerical_analyses_3d) then
!$OMP SINGLE
call do_advection_3d(dt,work3d,uuadv,vvadv,wwadv,huadv,hvadv,phadv,phadv, &
vel3d_adv_split,vel3d_adv_hor,vel3d_adv_ver,_ZERO_,U_TAG, &
hires=hires)
vel3d_adv_split,vel3d_adv_hor,vel3d_adv_ver,_ZERO_,U_TAG)
numdis_int = _ZERO_
!$OMP END SINGLE
......@@ -250,7 +281,7 @@
do j=jmin,jmax
#endif
do i=imin,imax
work3d(i,j,k) = ( work3d(i,j,k) - fadv3d(i,j,k)**2 ) / dt
work3d(i,j,k) = ( work3d(i,j,k) - fadv3d(i,j,k)**2 )/dt*hires(i,j,k)/ARUD1
end do
#ifndef SLICE_MODEL
end do
......@@ -271,10 +302,9 @@
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis_3d(i,j,k) = _HALF_*( work3d(i-1,j,k) + work3d(i,j,k) )
numdis_int(i,j) = numdis_int(i,j) &
+_HALF_*( work3d(i-1,j,k)*hires(i-1,j,k) &
+work3d(i ,j,k)*hires(i ,j,k) )
numdis_3d_old(i,j,k) = _HALF_*( work3d(i-1,j,k) + work3d(i,j,k) )/hn(i,j,k)*ARCD1
numdis_int(i,j) = numdis_int(i,j) &
+_HALF_*( work3d(i-1,j,k) + work3d(i,j,k) )*ARCD1
end if
end do
#ifndef SLICE_MODEL
......@@ -286,6 +316,21 @@
end if
#endif
#ifdef _MOMENTUM_TERMS_
do k=1,kmax
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
adv_u(i,j,k) = dry_u(i,j) * uuEx(i,j,k)
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO NOWAIT
end do
#endif
! Here begins dimensional split advection for v-velocity
......@@ -416,24 +461,47 @@
!$OMP SINGLE
call do_advection_3d(dt,fadv3d,uuadv,vvadv,wwadv,huadv,hvadv,phadv,phadv, &
vel3d_adv_split,vel3d_adv_hor,vel3d_adv_ver,_ZERO_,V_TAG, &
advres=vvEx,nvd=numdis_v_3d)
hires=hires,advres=vvEx,nvd=p_nvd)
!$OMP END SINGLE
#ifdef _MOMENTUM_TERMS_
do k=1,kmax
if (calc_numdis) then
do k=1,kmax
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
do j=jmin,jmax
#endif
do i=imin,imax
adv_v(i,j,k) = dry_v(i,j) * vvEx(i,j,k)
end do
do i=imin,imax
nvd(i,j,k) = nvd(i,j,k)*hires(i,j,k)/ARVD1
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end do
!$OMP SINGLE
call update_3d_halo(nvd,nvd,av,imin,jmin,imax,jmax,kmax,V_TAG)
call wait_halo(V_TAG)
!$OMP END SINGLE
do k=1,kmax
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
!$OMP END DO NOWAIT
end do
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis_3d(i,j,k) = numdis_3d(i,j,k) &
+_HALF_*( nvd(i,j-1,k) + nvd(i,j,k) )/hn(i,j,k)*ARCD1
end if
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO
end do
end if
#ifdef _NUMERICAL_ANALYSES_OLD_
if (do_numerical_analyses_3d) then
......@@ -450,7 +518,7 @@
do j=jmin,jmax
#endif
do i=imin,imax
work3d(i,j,k) = ( work3d(i,j,k) - fadv3d(i,j,k)**2 ) / dt
work3d(i,j,k) = ( work3d(i,j,k) - fadv3d(i,j,k)**2 )/dt*hires(i,j,k)/ARVD1
end do
#ifndef SLICE_MODEL
end do
......@@ -475,11 +543,10 @@
#endif
do i=imin,imax
if (az(i,j) .eq. 1) then
numdis_3d(i,j,k) = numdis_3d(i,j,k) &
+_HALF_*( work3d(i,j-1,k) + work3d(i,j,k) )
numdis_int(i,j) = numdis_int(i,j) &
+_HALF_*( work3d(i,j-1,k)*hires(i,j-1,k) &
+work3d(i,j ,k)*hires(i,j ,k) )
numdis_3d_old(i,j,k) = numdis_3d_old(i,j,k) &
+_HALF_*( work3d(i,j-1,k) + work3d(i,j,k) )/hn(i,j,k)*ARCD1
numdis_int(i,j) = numdis_int(i,j) &
+_HALF_*( work3d(i,j-1,k) + work3d(i,j,k) )*ARCD1
end if
end do
#ifndef SLICE_MODEL
......@@ -491,6 +558,22 @@
end if
#endif
#ifdef _MOMENTUM_TERMS_
do k=1,kmax
!$OMP DO SCHEDULE(RUNTIME)
#ifndef SLICE_MODEL
do j=jmin,jmax
#endif
do i=imin,imax
adv_v(i,j,k) = dry_v(i,j) * vvEx(i,j,k)
end do
#ifndef SLICE_MODEL
end do
#endif
!$OMP END DO NOWAIT
end do
#endif
!$OMP END PARALLEL
end if
......
src/3d/variables_3d.F90
View file @ c2406931
......@@ -126,9 +126,8 @@
#include "dynamic_declarations_3d.h"
#endif
REALTYPE,dimension(:,:,:),pointer :: numdis_u_3d=>null()
REALTYPE,dimension(:,:,:),pointer :: numdis_v_3d=>null()
REALTYPE, dimension(:,:,:), allocatable :: numdis_3d,phydis_3d
REALTYPE,dimension(:,:,:),pointer :: numdis_3d=>null()
REALTYPE, dimension(:,:,:), allocatable :: numdis_3d_old,phydis_3d
REALTYPE, dimension(:,:), allocatable :: numdis_int,phydis_int
REALTYPE,dimension(:,:,:),pointer :: nummix_S=>null()
REALTYPE,dimension(:,:,:),pointer :: nummix_T=>null()
......
src/ncdf/init_2d_ncdf.F90
View file @ c2406931
......@@ -120,23 +120,17 @@
mv = nummix_missing
vr(1) = -100.0
vr(2) = 100.0
err = nf90_def_var(ncid,'numdis_u_2d',NCDF_FLOAT_PRECISION,f3_dims,ndu2d_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,ndu2d_id, &
long_name='numerical dissipation (u)', &
units='W/kg',&
FillValue=fv,missing_value=mv,valid_range=vr)
err = nf90_def_var(ncid,'numdis_v_2d',NCDF_FLOAT_PRECISION,f3_dims,ndv2d_id)
err = nf90_def_var(ncid,'numdis_2d',NCDF_FLOAT_PRECISION,f3_dims,nd2d_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,ndv2d_id, &
long_name='numerical dissipation (v)', &
call set_attributes(ncid,nd2d_id, &
long_name='numerical dissipation', &
units='W/kg',&
FillValue=fv,missing_value=mv,valid_range=vr)
#ifdef _NUMERICAL_ANALYSES_OLD_
err = nf90_def_var(ncid,'numdis_2d',NCDF_FLOAT_PRECISION,f3_dims,nd2d_id)
err = nf90_def_var(ncid,'numdis_2d_old',NCDF_FLOAT_PRECISION,f3_dims,nd2do_id)
if (err .NE. NF90_NOERR) go to 10
call set_attributes(ncid,nd2d_id, &
long_name='numerical dissipation', &
call set_attributes(ncid,nd2do_id, &
long_name='numerical dissipation (old)', &
units='W/kg',&
FillValue=fv,missing_value=mv,valid_range=vr)