Merge pull request #253 from matsbn/feature-diffusivity_anisot3D

New mesoscale eddy diffusivity options

cime_config/buildnml

@@ -178,6 +178,7 @@ set DKTZL                  = 1
 set EITMTH   = "'gm'"
 set EDRITP   = "'large scale'"
 set EDWMTH   = "'smooth'"
+set EDDF2D   = .false.
 set EDSPRS   = .true.
 set EGC      = 0.85
 set EGGAM    = 200.
@@ -191,6 +192,11 @@ else
   set EGMXDF   = 1500.
 endif
 set EGIDFQ   = 1.
+set TBFILE   = "'unset'"
+set RHISCF   = 0.
+set EDANIS   = .false.
+set REDI3D   = .false.
+set RHSCTP   = .false.
 set RI0      = 1.2
 set BDMTYP   = 2
 if ($BLOM_UNIT == cgs) then
@@ -1121,6 +1127,7 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
 !            horizontal grid spacing. Valid methods: 'smooth', 'step' (a)
 ! MLRTTP   : Type of mixed layer restratification time scale. Valid
 !            types: 'variable', 'constant', 'limited' (a)
+! EDDF2D   : If true, eddy diffusivity has a 2d structure (l)
 ! EDSPRS   : Apply eddy mixing suppression away from steering level (l)
 ! EGC      : Parameter c in Eden and Greatbatch (2008) parameterization (f)
 ! EGGAM    : Parameter gamma in E. & G. (2008) param. (f)
@@ -1130,6 +1137,13 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
 ! EGIDFQ   : Factor relating the isopycnal diffusivity to the layer
 !            interface diffusivity in the Eden and Greatbatch (2008)
 !            parameterization. egidfq=difint/difiso () (f)
+! TBFILE   : Name of file containing topographic beta parameter (a)
+! RHISCF   : Linear scaling parameter for topographic rhines scale () (f)
+! EDANIS   : If true, apply anisotropy correction to eddy diffusivity (l)
+! REDI3D   : If true, then isopycnal/neutral diffusion will have 3D
+!            structure based in the 3D structure of anisotropy (l)
+! RHSCTP   : If true, use the minimum of planetary and topographic beta
+!            to define the Rhines scale (l)
 ! RI0      : Critical gradient richardson number for shear driven
 !            vertical mixing () (f)
 ! BDMTYP   : Type of background diapycnal mixing. If bdmtyp=1 the
@@ -1149,13 +1163,19 @@ cat >>! $RUNDIR/ocn_in$inststr << EOF
   EITMTH   = $EITMTH
   EDRITP   = $EDRITP
   EDWMTH   = $EDWMTH
+  EDDF2D   = $EDDF2D
   EDSPRS   = $EDSPRS
   EGC      = $EGC
   EGGAM    = $EGGAM
   EGLSMN   = $EGLSMN
   EGMNDF   = $EGMNDF
   EGMXDF   = $EGMXDF
   EGIDFQ   = $EGIDFQ
+  TBFILE   = $TBFILE
+  RHISCF   = $RHISCF
+  EDANIS   = $EDANIS
+  REDI3D   = $REDI3D
+  RHSCTP   = $RHSCTP
   RI0      = $RI0
   BDMTYP   = $BDMTYP
   BDMC1    = $BDMC1

phy/geoenv_file.F

@@ -1,5 +1,6 @@
 ! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2022 Mats Bentsen, Ping-Gin Chiu, Mehmet Ilicak
+! Copyright (C) 2015-2023 Mats Bentsen, Ping-Gin Chiu, Mehmet Ilicak,
+!                         Aleksi Nummelin
 !
 ! This file is part of BLOM.
 !
@@ -27,12 +28,13 @@ subroutine geoenv_file
       use mod_config, only: inst_suffix
       use mod_constants, only: rearth, pi, radian, L_mks2cgs
       use mod_xc
+      use mod_diffusion, only: rhsctp, tbfile
       use mod_grid, only: grfile, qclon, qclat, pclon, pclat, uclon,
      .                    uclat, vclon, vclat, scqx, scqy, scpx, scpy,
      .                    scux, scuy, scvx, scvy, scq2, scp2, scu2,
      .                    scv2, qlon, qlat, plon, plat, ulon, ulat,
      .                    vlon, vlat, depths, corioq, coriop, betafp,
-     .                    angle, cosang, sinang, nwp
+     .                    betatp, angle, cosang, sinang, hangle, nwp
       use netcdf
 c
       implicit none
@@ -49,6 +51,9 @@ subroutine geoenv_file
       integer, dimension(3) :: start,count
       integer i,j,k,status,ncid,dimid,varid,ios,ncwm,l
       logical fexist
+c
+      real, parameter ::
+     .  iL_mks2cgs = 1./L_mks2cgs
 c
       namelist /cwmod/ cwmtag,cwmedg,cwmi,cwmj,cwmwth
 c
@@ -700,6 +705,69 @@ subroutine geoenv_file
       endif
 c
 c --- ------------------------------------------------------------------
+c --- read topographic beta if needed
+c --- ------------------------------------------------------------------
+c
+      if (rhsctp) then
+c
+        if (mnproc.eq.1) then
+          write (lp,'(2a)') ' reading topographic beta from ',
+     .                      trim(tbfile)
+          call flush(lp)
+          status=nf90_open(tbfile,nf90_nowrite,ncid)
+          if (status.ne.nf90_noerr) then
+            write(lp,'(4a)') ' nf90_open: ',trim(tbfile),': ',
+     .                       nf90_strerror(status)
+            call xchalt('(geoenv_file)')
+                   stop '(geoenv_file)'
+          endif
+          status=nf90_inq_varid(ncid,'topo_beta',varid)
+          if (status.ne.nf90_noerr) then
+            write(lp,'(2a)') ' nf90_inq_varid: topo_beta: ',
+     .                       nf90_strerror(status)
+            call xchalt('(geoenv_file)')
+                   stop '(geoenv_file)'
+          endif
+          status=nf90_get_var(ncid,varid,tmpg)
+          if (status.ne.nf90_noerr) then
+            write(lp,'(2a)') ' nf90_get_var: topo_beta: ',
+     .                       nf90_strerror(status)
+            call xchalt('(geoenv_file)')
+                   stop '(geoenv_file)'
+          endif
+        endif
+        call xcaput(tmpg,betatp,1)
+c
+        if (mnproc.eq.1) then
+          status=nf90_inq_varid(ncid,'hangle',varid)
+          if (status.ne.nf90_noerr) then
+            write(lp,'(2a)') ' nf90_inq_varid: hangle: ',
+     .                       nf90_strerror(status)
+            call xchalt('(geoenv_file)')
+                   stop '(geoenv_file)'
+          endif
+          status=nf90_get_var(ncid,varid,tmpg)
+          if (status.ne.nf90_noerr) then
+            write(lp,'(2a)') ' nf90_get_var: hangle: ',
+     .                       nf90_strerror(status)
+            call xchalt('(geoenv_file)')
+                   stop '(geoenv_file)'
+          endif
+        endif
+        call xcaput(tmpg,hangle,1)
+c
+        if (mnproc.eq.1) then
+          status=nf90_close(ncid)
+          if (status.ne.nf90_noerr) then
+            write(lp,'(4a)') ' nf90_close: ',trim(tbfile),': ',
+     .                       nf90_strerror(status)
+            call xchalt('(geoenv_file)')
+                   stop '(geoenv_file)'
+          endif
+        endif
+c
+      endif
+c --- ------------------------------------------------------------------
 c --- Apply channel width modifications if specified in namelist
 c --- ------------------------------------------------------------------
 c
@@ -788,9 +856,9 @@ subroutine geoenv_file
       endif
 c
 c --- ------------------------------------------------------------------
-c --- Get correct units of scale factors, precompute cosine and sine of
-c --- local angle of i-direction and with eastward direction, and
-c --- compute Coriolis and beta plane parameter
+c --- Get correct units of scale factors and topographic beta,
+c --- precompute cosine and sine of local angle of i-direction and with
+c --- eastward direction, and compute Coriolis and beta plane parameter
 c --- ------------------------------------------------------------------
 c
 c$OMP PARALLEL DO PRIVATE(i)
@@ -820,6 +888,16 @@ subroutine geoenv_file
         enddo
       enddo
 c$OMP END PARALLEL DO
+c
+      if (rhsctp) then
+c$OMP PARALLEL DO PRIVATE(i)
+        do j=1,jj
+          do i=1,ii
+            betatp(i,j)=betatp(i,j)*iL_mks2cgs
+          enddo
+        enddo
+c$OMP END PARALLEL DO
+      endif
 c
       return
       end

phy/geoenv_test.F

@@ -1,5 +1,5 @@
 ! ------------------------------------------------------------------------------
-! Copyright (C) 2015-2020 Mats Bentsen
+! Copyright (C) 2015-2023 Mats Bentsen, Aleksi Nummelin
 !
 ! This file is part of BLOM.
 !
@@ -28,8 +28,8 @@ subroutine geoenv_test
      .                    vclon, vclat, scqx, scqy, scpx, scpy, scux,
      .                    scuy, scvx, scvy, scq2, scp2, scu2, scv2,
      .                    qlon, qlat, plon, plat, ulon, ulat, vlon,
-     .                    vlat, depths, corioq, coriop, betafp, angle,
-     .                    cosang, sinang, nwp
+     .                    vlat, depths, corioq, coriop, betafp, betatp,
+     .                    angle, cosang, sinang, hangle, nwp
 c
       implicit none
 c
@@ -67,8 +67,10 @@ subroutine geoenv_test
       corioq=0.
       coriop=0.
       betafp=0.
+      betatp=0.
       cosang=0.
       sinang=0.
+      hangle=0.
 c
       return
       end

phy/mod_advect.F

@@ -93,10 +93,6 @@ subroutine advect(m,n,mm,nn,k1m,k1n)
       enddo
 c$OMP END PARALLEL DO
 c
-      call xctilr(ubflxs_p, 1,2, 2,2, halo_uv)
-      call xctilr(vbflxs_p, 1,2, 2,2, halo_vv)
-      call xctilr(pbu, 1,2, 2,2, halo_us)
-      call xctilr(pbv, 1,2, 2,2, halo_vs)
 #ifdef TRC
       do nt=1,ntr
 #  if defined(TKE) && !defined(TKEADV)

phy/mod_dia.F

@@ -1096,7 +1096,7 @@ subroutine diaacc(m,n,mm,nn,k1m,k1n)
 c$OMP END PARALLEL DO
       endif
 c
-      if (sum(ACC_MLD(1:nphy)).ne.0) then
+      if (sum(ACC_MLD(1:nphy)+ACC_MAXMLD(1:nphy)).ne.0) then
         select case (vcoord_type_tag)
           case (isopyc_bulkml)
 c$OMP PARALLEL DO PRIVATE(l,i)