numerical dissipation analysis added

src/3d/dynamic_allocations_3d.h

@@ -97,6 +97,12 @@
   allocate(nummix2d_T(I2DFIELD),stat=rc)
   if (rc /= 0) stop 'init_3d: Error allocating memory (nummix2d_T)'
 
+  allocate(numdis3d(I3DFIELD),stat=rc)
+  if (rc /= 0) stop 'init_3d: Error allocating memory (numdis3d)'
+
+  allocate(numdis2d(I2DFIELD),stat=rc)
+  if (rc /= 0) stop 'init_3d: Error allocating memory (numdis2d)'
+
   allocate(phymix3d_S(I3DFIELD),stat=rc)
   if (rc /= 0) stop 'init_3d: Error allocating memory (phymix3d_S)'
 

src/3d/dynamic_declarations_3d.h

@@ -38,6 +38,8 @@
   REALTYPE, dimension(:,:,:), allocatable   :: phymix3d_S,phymix3d_T
   REALTYPE, dimension(:,:), allocatable     :: nummix2d_S,nummix2d_T
   REALTYPE, dimension(:,:), allocatable     :: phymix2d_S,phymix2d_T
+  REALTYPE, dimension(:,:,:), allocatable   :: numdis3d
+  REALTYPE, dimension(:,:), allocatable     :: numdis2d
 #endif
 
 ! suspended matter

src/3d/static_3d.h

@@ -61,6 +61,8 @@
   REALTYPE                             :: nummix2d_S(I2DFIELD)
   REALTYPE                             :: nummix3d_T(I3DFIELD)
   REALTYPE                             :: nummix2d_T(I2DFIELD)
+  REALTYPE                             :: numdis3d(I3DFIELD)
+  REALTYPE                             :: numdis2d(I2DFIELD)
   REALTYPE                             :: phymix3d_S(I3DFIELD)
   REALTYPE                             :: phymix2d_S(I2DFIELD)
   REALTYPE                             :: phymix3d_T(I3DFIELD)

src/3d/uv_advect_3d.F90

@@ -242,11 +242,19 @@
 ! (\ref{SxA}) and (\ref{SyA}).
 !
 ! With this method, all higher-order directional-split advection schemes
-! are now available for the momentum advection. The advective
+! are available for the momentum advection. The advective
 ! fluxes needed for this have to be averaged from the conservative
 ! advective fluxes resulting from the continuity equation
-! Continuity will
+! (\ref{ContiLayerInt}). Continuity will
 ! still be retained due to the linearity of the continuity equation.
+! 
+! \paragraph{Numerical dissipation.}\label{uvadvect-dissipation}
+!
+! For the directional split method, numerical dissipation is calculated
+! if {\tt do\_mixing\_analysis} is set to {\tt .true.},
+! using the method suggested by \cite{BURCHARD12}.
+!
+!
 !
 ! !
 ! !USES:
@@ -263,6 +271,9 @@
    use advection_3d, only: do_advection_3d
    use halo_zones, only: update_3d_halo,wait_halo,U_TAG,V_TAG
    use getm_timers, only: tic, toc, TIM_UVADV3D, TIM_UVADV3DH
+   use variables_3d, only: do_mixing_analysis
+   use variables_3d, only: numdis3d,numdis2d
+
 !$ use omp_lib
    IMPLICIT NONE
 !
@@ -287,6 +298,9 @@
    REALTYPE                  :: area_inv(I2DFIELD)
    REALTYPE                  :: AH=_ZERO_
    REALTYPE                  :: dti,dxdyi
+   REALTYPE                  :: vel2(I3DFIELD)
+   REALTYPE                  :: vel2o(I3DFIELD)
+   REALTYPE                  :: numdiss(I3DFIELD)
 #endif
 !EOP
 !-----------------------------------------------------------------------
@@ -363,16 +377,53 @@
       end do
    end do
 
-
    call tic(TIM_UVADV3DH)
    call update_3d_halo(uuEx,uuEx,au,imin,jmin,imax,jmax,kmax,U_TAG)
    call wait_halo(U_TAG)
    call toc(TIM_UVADV3DH)
 
+   if (do_mixing_analysis) then
+      do k=1,kmax ! calculate square of u-velocity before advection step 
+         do j=jmin,jmax
+            do i=imin,imax
+               vel2(i,j,k)=uuEx(i,j,k)**2 
+            end do
+         end do
+      end do
+      vel2o=vel2
+   end if   
+
    call do_advection_3d(dt,uuEx,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv,&
                         dxuadv,dxvadv,dyuadv,dyvadv,area_inv,          &
                         azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
 
+   if (do_mixing_analysis) then
+      call do_advection_3d(dt,vel2,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv,&
+                           dxuadv,dxvadv,dyuadv,dyvadv,area_inv,          &
+                           azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
+
+      do k=1,kmax ! calculate kinetic energy dissipaion rate for u-velocity
+         do j=jmin,jmax
+            do i=imin,imax
+               numdiss(i,j,k)=(vel2(i,j,k)-uuEx(i,j,k)**2)/dt
+            end do
+         end do
+      end do
+
+      numdis2d=_ZERO_
+      do k=1,kmax ! calculate kinetic energy dissipaion rate for u-velocity
+         do j=jmin,jmax
+            do i=imin,imax
+               numdis3d(i,j,k)=0.5*(numdiss(i,j,k)+numdiss(i-1,j,k))
+               numdis2d(i,j)=numdis2d(i,j)                             &
+                             +0.5*(numdiss(i,j,k)*hun(i,j,k)           &
+                                  +numdiss(i-1,j,k)*hun(i-1,j,k))
+            end do
+         end do
+      end do
+   end if
+   
+
 !$OMP PARALLEL DEFAULT(SHARED) PRIVATE(i,j,k)
 
 
@@ -442,10 +493,47 @@
    call wait_halo(V_TAG)
    call toc(TIM_UVADV3DH)
 
+   if (do_mixing_analysis) then
+      do k=1,kmax ! calculate square of v-velocity before advection step 
+         do j=jmin,jmax
+            do i=imin,imax
+               vel2(i,j,k)=vvEx(i,j,k)**2 
+            end do
+         end do
+      end do
+      vel2o=vel2
+   end if
+
    call do_advection_3d(dt,vvEx,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv,&
                         dxuadv,dxvadv,dyuadv,dyvadv,area_inv,          &
                         azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
 
+   if (do_mixing_analysis) then
+      call do_advection_3d(dt,vel2,uadv,vadv,wadv,huadv,hvadv,hoadv,hnadv,&
+                           dxuadv,dxvadv,dyuadv,dyvadv,area_inv,          &
+                           azadv,auadv,avadv,hor_adv,ver_adv,adv_split,AH)
+
+      do k=1,kmax ! calculate kinetic energy dissipaion rate for u-velocity
+         do j=jmin,jmax
+            do i=imin,imax
+               numdiss(i,j,k)=(vel2(i,j,k)-vvEx(i,j,k)**2)/dt
+            end do
+         end do
+      end do
+
+      do k=1,kmax ! calculate kinetic energy dissipaion rate for u-velocity
+         do j=jmin,jmax
+            do i=imin,imax
+               numdis3d(i,j,k)=numdis3d(i,j,k)                         &
+                              +0.5*(numdiss(i,j,k)+numdiss(i,j-1,k))
+               numdis2d(i,j)=numdis2d(i,j)                             &
+                             +0.5*(numdiss(i,j,k)*hvn(i,j,k)           &
+                                  +numdiss(i,j-1,k)*hvn(i,j-1,k))
+            end do
+         end do
+      end do
+   end if
+
 ! OMP-NOTE: It might not pay off to thread this loop (due to OMP overhead)
 !   vvEx=-(vvEx*hvn-vv)*dti ! Here, vvEx is the advection term.
 

src/3d/variables_3d.F90

@@ -209,6 +209,7 @@
    idpdy=_ZERO_
    nummix3d_S = _ZERO_ ; nummix2d_S = _ZERO_
    nummix3d_T = _ZERO_ ; nummix2d_T = _ZERO_
+   numdis3d   = _ZERO_ ; numdis2d   = _ZERO_
 #endif
 
    adv_schemes(1) = "3D first-order upstream advection"

src/ncdf/init_3d_ncdf.F90

@@ -318,6 +318,12 @@
       mv = nummix_missing
       vr(1) = -100.0
       vr(2) = 100.0
+      err = nf90_def_var(ncid,'numdis3d',NCDF_FLOAT_PRECISION,f4_dims,nm3d_id)
+      if (err .NE. NF90_NOERR) go to 10
+      call set_attributes(ncid,nm3d_id, &
+          long_name='numerical dissipation', &
+          units='W/kg',&
+          FillValue=fv,missing_value=mv,valid_range=vr)
 
       if (calc_salt) then
          err = nf90_def_var(ncid,'nummix3d_S',NCDF_FLOAT_PRECISION,f4_dims,nm3dS_id)

src/ncdf/init_mean_ncdf.F90

@@ -170,6 +170,19 @@
       mv = nummix_missing
       vr(1) = -100.0
       vr(2) =  100.0
+      err = nf90_def_var(ncid,'numdis3d',NCDF_FLOAT_PRECISION,f4_dims,nm3d_id)
+      if (err .NE. NF90_NOERR) go to 10
+      call set_attributes(ncid,nm3d_id, &
+          long_name='mean numerical dissipation', &
+          units='W/kg',&
+          FillValue=fv,missing_value=mv,valid_range=vr)
+
+      err = nf90_def_var(ncid,'numdis2d',NCDF_FLOAT_PRECISION,f3_dims,nm2d_id)
+      if (err .NE. NF90_NOERR) go to 10
+      call set_attributes(ncid,nm2d_id, &
+          long_name='mean, vert. integrated numerical dissipation', &
+          units='Wm/kg',&
+          FillValue=fv,missing_value=mv,valid_range=vr)
 
       if (calc_salt) then
          err = nf90_def_var(ncid,'nummix3d_S',NCDF_FLOAT_PRECISION,f4_dims,nm3dS_id)

src/ncdf/ncdf_3d.F90

@@ -39,6 +39,7 @@
 #endif
    integer                             :: nm3dS_id,nm3dT_id,nm2dS_id,nm2dT_id
    integer                             :: pm3dS_id,pm3dT_id,pm2dS_id,pm2dT_id
+   integer                             :: nm3d_id,nm2d_id
 
    REALTYPE, dimension(:,:,:), allocatable :: ws
 

src/ncdf/ncdf_mean.F90

@@ -25,6 +25,7 @@
    integer                             :: saltmean_id,tempmean_id,hmean_id
    integer                             :: nm3dS_id,nm3dT_id,nm2dS_id,nm2dT_id
    integer                             :: pm3dS_id,pm3dT_id,pm2dS_id,pm2dT_id
+   integer                             :: nm3d_id,nm2d_id
 #ifdef GETM_BIO
    integer, allocatable                :: biomean_id(:)
 #endif

src/ncdf/save_3d_ncdf.F90

@@ -29,6 +29,7 @@
 #ifndef NO_BAROCLINIC
    use variables_3d, only: S,T,rho,rad,NN
    use variables_3d, only: nummix3d_S,nummix3d_T,phymix3d_S,phymix3d_T
+   use variables_3d, only: numdis3d
 #endif
    use variables_3d, only: tke,num,nuh,eps
 #ifdef SPM
@@ -293,6 +294,11 @@
 
 #ifndef NO_BAROCLINIC
    if (save_mix_analysis) then
+      call cnv_3d(imin,jmin,imax,jmax,kmin,kmax,az,numdis3d,nummix_missing, &
+                  imin,imax,jmin,jmax,0,kmax,ws)
+      err = nf90_put_var(ncid,nm3d_id,ws(_3D_W_),start,edges)
+      if (err .NE. NF90_NOERR) go to 10
+
       if (calc_salt) then
          call cnv_3d(imin,jmin,imax,jmax,kmin,kmax,az,nummix3d_S,nummix_missing, &
                      imin,imax,jmin,jmax,0,kmax,ws)

src/ncdf/save_mean_ncdf.F90

@@ -88,6 +88,7 @@
    edges(3) = zlen
    edges(4) = 1
 
+
 !  layer thickness
    if (hmean_id .gt. 0) then
       call cnv_3d(imin,jmin,imax,jmax,kmin,kmax,az,hmean,h_missing, &
@@ -132,6 +133,11 @@
    end if
 
    if (save_mix_analysis) then
+      call cnv_3d(imin,jmin,imax,jmax,kmin,kmax,az, &
+                  numdis3d_mean,nummix_missing, &
+                  imin,imax,jmin,jmax,0,kmax,ws3d)
+      err = nf90_put_var(ncid, nm3d_id,ws3d(_3D_W_),start,edges)
+      if (err .NE. NF90_NOERR) go to 10
       if (calc_salt) then
          call cnv_3d(imin,jmin,imax,jmax,kmin,kmax,az, &
                      nummix3d_S_mean,nummix_missing, &
@@ -167,6 +173,11 @@
       edges(2) = ylen
       edges(3) = 1
 
+      call cnv_2d(imin,jmin,imax,jmax,az,numdis2d_mean,nummix_missing, &
+                  imin,jmin,imax,jmax,ws2d)
+      err = nf90_put_var(ncid, nm2d_id,ws2d(_2D_W_),start,edges)
+      if (err .NE. NF90_NOERR) go to 10
+
       if (calc_salt) then
          call cnv_2d(imin,jmin,imax,jmax,az,nummix2d_S_mean,nummix_missing, &
                      imin,jmin,imax,jmax,ws2d)

src/output/calc_mean_fields.F90

@@ -20,6 +20,7 @@
    use variables_3d, only: S,T
    use variables_3d, only: nummix3d_S,nummix2d_S,nummix3d_T,nummix2d_T
    use variables_3d, only: phymix3d_S,phymix2d_S,phymix3d_T,phymix2d_T
+   use variables_3d, only: numdis3d,numdis2d
 #endif
 #ifdef GETM_BIO
    use bio, only: bio_calc
@@ -89,6 +90,12 @@
           stop 'calc_mean_fields.F90: Error allocating memory (Smean)'
 
       if (do_mixing_analysis) then
+         allocate(numdis3d_mean(I3DFIELD),stat=rc)
+           if (rc /= 0) &
+              stop 'calc_mean_fields.F90: Error allocating memory (numdis3d_mean)'
+         allocate(numdis2d_mean(I2DFIELD),stat=rc)
+           if (rc /= 0) &
+              stop 'calc_mean_fields.F90: Error allocating memory (numdis2d_mean)'
          if (calc_temp) then
             allocate(nummix3d_T_mean(I3DFIELD),stat=rc)
             if (rc /= 0) &
@@ -133,6 +140,7 @@
 #ifndef NO_BAROCLINIC
       Tmean=_ZERO_; Smean=_ZERO_
       if (do_mixing_analysis) then
+         numdis3d_mean=_ZERO_; numdis2d_mean=_ZERO_
          if (calc_temp) then
             nummix3d_T_mean=_ZERO_; nummix2d_T_mean=_ZERO_
             phymix3d_T_mean=_ZERO_; phymix2d_T_mean=_ZERO_
@@ -180,6 +188,8 @@
       Tmean = Tmean + T
       Smean = Smean + S
       if (do_mixing_analysis) then
+         numdis3d_mean = numdis3d_mean + numdis3d
+         numdis2d_mean = numdis2d_mean + numdis2d
          if (calc_temp) then
             nummix3d_T_mean = nummix3d_T_mean + nummix3d_T
             nummix2d_T_mean = nummix2d_T_mean + nummix2d_T
@@ -216,6 +226,8 @@
          Tmean = Tmean / step
          Smean = Smean / step
          if (do_mixing_analysis) then
+            numdis3d_mean = numdis3d_mean / step
+            numdis2d_mean = numdis2d_mean / step
             if (calc_temp) then
                nummix3d_T_mean = nummix3d_T_mean / step
                nummix2d_T_mean = nummix2d_T_mean / step

src/output/diagnostic_variables.F90

@@ -32,6 +32,8 @@
    REALTYPE,dimension(:,:), allocatable :: nummix2d_S_mean
    REALTYPE,dimension(:,:,:), allocatable :: nummix3d_T_mean
    REALTYPE,dimension(:,:), allocatable :: nummix2d_T_mean
+   REALTYPE,dimension(:,:,:), allocatable :: numdis3d_mean
+   REALTYPE,dimension(:,:), allocatable :: numdis2d_mean
    REALTYPE,dimension(:,:,:), allocatable :: phymix3d_S_mean
    REALTYPE,dimension(:,:), allocatable :: phymix2d_S_mean
    REALTYPE,dimension(:,:,:), allocatable :: phymix3d_T_mean