Double gyre setup

config/namelist.config.toy_dbgyre

@@ -0,0 +1,79 @@
+! This is the namelist file for model general configuration
+
+&modelname
+runid='fesom'
+/
+
+&timestep
+step_per_day=45 !96 !96 !72 !72 !45 !72 !96
+run_length=1 !62 !62 !62 !28
+run_length_unit='d' ! y, m, d, s
+/
+
+&clockinit              ! the model starts at
+timenew=0.0
+daynew=1
+yearnew=1900
+/
+
+&paths
+MeshPath=''
+ResultPath='../results_tmp/'
+/
+
+&restart_log
+restart_length=1           ! --> do netcdf restart ( only required for d,h,s cases, y, m take 1)
+restart_length_unit='y'     !output period: y, d, h, s, off
+raw_restart_length=1        ! --> do core dump restart
+raw_restart_length_unit='off' ! e.g. y, d, h, s, off
+bin_restart_length=1        ! --> do derived type binary restart 
+bin_restart_length_unit='off' ! e.g. y, d, h, s, off
+logfile_outfreq=72     !in logfile info. output frequency, # steps
+/
+
+&ale_def
+which_ALE='linfs'       ! 'linfs','zlevel', 'zstar'
+use_partial_cell=.false.
+/
+
+&geometry
+cartesian=.false.
+fplane=.false.
+cyclic_length=90       ![degree]
+rotated_grid=.false.     !option only valid for coupled model case now
+force_rotation=.false.
+alphaEuler=0          ![degree] Euler angles, convention:
+betaEuler=0           ![degree] first around z, then around new x,
+gammaEuler=0         ![degree] then around new z.
+/
+
+&calendar
+include_fleapyear=.false.
+/
+
+&run_config
+use_ice=.false.                  ! ocean+ice
+use_cavity=.false.              !
+use_cavity_partial_cell=.false. 
+use_floatice = .false.
+use_sw_pene=.true.
+flag_debug=.false.
+flag_warn_cflz=.false.
+toy_ocean=.true.
+which_toy="dbgyre"
+flag_debug=.false.
+/
+
+&machine
+n_levels=1
+n_part=24
+/
+
+&icebergs
+use_icesheet_coupling=.false.
+ib_num=1
+use_icebergs=.false.
+steps_per_ib_step=8
+ib_async_mode=0
+/
+

config/namelist.dyn.toy_dbgyre

@@ -0,0 +1,25 @@
+&dynamics_visc
+visc_gamma0  = 0.005 ! [m/s],   backgroung viscosity= gamma0*len, it should be as small a s possible (keep it < 0.01 m/s).
+visc_gamma1  = 0.3   ! [nodim], for computation of the flow aware viscosity
+visc_gamma2  = 0.20 ! [s/m],   is only used in easy backscatter option
+visc_easybsreturn= 0.0
+
+opt_visc     = 5    
+check_opt_visc=.false.  ! check if optvisc=5 is valid based on ratio resol/rossbyR
+! 5=Kinematic (easy) Backscatter
+! 6=Biharmonic flow aware (viscosity depends on velocity Laplacian)
+! 7=Biharmonic flow aware (viscosity depends on velocity differences)
+! 8=Dynamic Backscatter
+
+use_ivertvisc= .true.  
+/
+
+&dynamics_general
+momadv_opt   = 2       ! option for momentum advection in moment only =2
+use_freeslip = .false. ! Switch on free slip
+use_wsplit   = .false. ! Switch for implicite/explicte splitting of vert. velocity
+wsplit_maxcfl= 1.0     ! maximum allowed CFL criteria in vertical (0.5 < w_max_cfl < 1.) 
+                       ! in older FESOM it used to be w_exp_max=1.e-3
+ldiag_KE=.false.       ! activates energy diagnostics
+AB_order=2
+/

config/namelist.oce.toy_dbgyre

@@ -0,0 +1,28 @@
+! The namelist file for the finite-volume ocean model
+
+&oce_dyn
+state_equation=0       ! 1 - full equation of state, 0 - linear equation of state
+C_d=0.001             ! Bottom drag, nondimensional
+A_ver= 1.e-4	       ! Vertical viscosity, m^2/s
+scale_area=5.8e9       ! Visc. and diffus. are for an element with scale_area
+SPP=.false.                 ! Salt Plume Parameterization
+Fer_GM=.false.               ! to swith on/off GM after Ferrari et al. 2010
+K_GM_max     = 2000.0       ! max. GM thickness diffusivity (m2/s)
+K_GM_min     = 2.0          ! max. GM thickness diffusivity (m2/s)
+K_GM_bvref   = 2            ! def of bvref in ferreira scaling 0=srf,1=bot mld,2=mean over mld,3=weighted mean over mld
+K_GM_rampmax = -1.0         ! Resol >K_GM_rampmax[km] GM on
+K_GM_rampmin = -1.0         ! Resol <K_GM_rampmin[km] GM off, in between linear scaled down
+K_GM_resscalorder = 1
+
+scaling_Ferreira   =.false.  ! GM vertical scaling after Ferreira et al.(2005) (as also implemented by Qiang in FESOM 1.4)
+scaling_Rossby     =.false. ! GM is smoothly switched off according to Rossby radius (from 1. in coarse areas to 0. where resolution reaches 2 points/Rossby radius)
+scaling_resolution =.true.  ! GM is spatially scaled with resolution; A value of K_GM corresponds then to a resolution of 100km
+scaling_FESOM14    =.false.  ! special treatment of GM in the NH (as also implemented by Qiang in FESOM 1.4; it is zero within the boundary layer)
+
+Redi  =.false.
+visc_sh_limit=5.0e-3       ! for KPP, max visc due to shear instability
+mix_scheme='PP'           ! vertical mixing scheme: KPP, PP 
+Ricr   = 0.3               ! critical bulk Richardson Number
+concv  = 1.6               ! constant for pure convection (eqn. 23) (Large 1.5-1.6; MOM default 1.8)
+/
+

config/namelist.tra.toy_dbgyre

@@ -0,0 +1,57 @@
+&tracer_listsize
+num_tracers=100 !number of tracers to allocate. shallbe large or equal to the number of streams in &nml_list
+/
+
+&tracer_list
+nml_tracer_list =  
+1  , 'MFCT',       'QR4C',       'FCT ',        1.,          1., 
+2  , 'MFCT',       'QR4C',       'FCT ',        1.,          1., 
+!101, 'UPW1',       'UPW1',       'NON ',        0.,          0.
+/
+
+&tracer_init3d                            ! initial conditions for tracers
+n_ic3d   = 2                              ! number of tracers to initialize
+idlist   = 2, 1                           ! their IDs (0 is temperature, 1 is salinity, etc.). The reading order is defined here!
+filelist = 'phc3.0_winter.nc', 'phc3.0_winter.nc' ! list of files in ClimateDataPath to read (one file per tracer), same order as idlist
+varlist  = 'salt', 'temp'                 ! variables to read from specified files
+t_insitu = .true.                         ! if T is insitu it will be converted to potential after reading it
+/
+
+&tracer_init2d                                      ! initial conditions for 2D tracers (sea ice)
+n_ic2d   = 3                                        ! number of tracers to initialize
+idlist   = 1, 2, 3                                  ! their IDs (0 is a_ice, 1 is m_ice, 3 m_snow). The reading order is defined here!
+filelist = 'a_ice.nc', 'm_ice.nc', 'm_snow.nc' ! list of files in ClimateDataPath to read (one file per tracer), same order as idlist
+varlist  = 'a_ice', 'm_ice', 'm_snow'                 ! variables to read from specified files
+ini_ice_from_file=.false.
+/
+
+&tracer_general
+! bharmonic diffusion for tracers. We recommend to use this option in very high resolution runs (Redi is generally off there).
+smooth_bh_tra =.false.     ! use biharmonic diffusion (filter implementation) for tracers
+gamma0_tra    = 0.0005     ! gammaX_tra are analogous to those in the dynamical part
+gamma1_tra    = 0.0125
+gamma2_tra    = 0.
+i_vert_diff   =.true.
+/
+
+&tracer_phys
+use_momix     = .false.     ! switch on/off  !Monin-Obukhov -> TB04 mixing
+momix_lat     = -50.0      ! latitidinal treshhold for TB04, =90 --> global
+momix_kv      = 0.01       ! PP/KPP, mixing coefficient within MO length
+use_instabmix = .true.     ! enhance convection in case of instable stratification
+instabmix_kv  = 0.1
+use_windmix   = .false.    ! enhance mixing trough wind only for PP mixing (for stability)
+windmix_kv    = 1.e-3
+windmix_nl    = 2
+diff_sh_limit=5.0e-3       ! for KPP, max diff due to shear instability
+Kv0_const=.false.
+double_diffusion=.false.   ! for KPP,dd switch
+K_ver=1.0e-5
+K_hor=10.
+surf_relax_T=0.0
+surf_relax_S=1.929e-06   ! 50m/300days 6.43e-07! m/s 10./(180.*86400.)
+balance_salt_water =.false. ! balance virtual-salt or freshwater flux or not
+clim_relax=0.0	           ! 1/s, geometrical information has to be supplied
+ref_sss_local=.true.
+ref_sss=34.
+/

env.sh

@@ -89,6 +89,8 @@ elif [[ $LOGINHOST =~ \.bullx$ ]]; then
    STRATEGY="atosecmwf"
 elif [[ $LOGINHOST =~ uan[0-9][0-9] ]]; then
    STRATEGY="lumi"
+elif [[ $LOGINHOST =~ nesh-login[1-3] ]]; then
+   STRATEGY="nesh"
 elif [[ -d $DIR/env/$LOGINHOST ]]; then # check if directory with LOGINHOST exists in env
 STRATEGY=$LOGINHOST
 else

env/nesh/shell

@@ -0,0 +1,15 @@
+# Inspired by environment for albedo
+module load oneapi2023-env/2023.2.0
+
+module load cmake/3.27.4
+module load oneapi/2023.2.0
+module load oneapi-mkl/2023.1.0
+module load oneapi-mpi/2021.10.0
+module load netcdf-c/4.9.2-with-oneapi-mpi-2021.10.0
+module load netcdf-fortran/4.6.0-with-oneapi-mpi-2021.10.0
+export FC=mpiifort CC=mpiicc CXX=mpiicpc
+
+
+# haven't set any environment variables. hopefully fine by now
+export I_MPI_PMI=pmi2
+export I_MPI_PMI_LIBRARY=/usr/lib64/libpmi2.so

src/gen_forcing_init.F90

@@ -10,6 +10,19 @@ subroutine forcing_array_setup(partit, mesh)
   end interface
 end module
 
+
+module forcing_array_setup_dbgyre_interfaces
+  interface
+    subroutine forcing_array_setup_dbgyre(partit, mesh)
+      use mod_mesh
+      USE MOD_PARTIT
+      USE MOD_PARSUP
+      type(t_mesh),   intent(in),    target :: mesh
+      type(t_partit), intent(inout), target :: partit
+    end subroutine
+  end interface
+end module
+
 ! Adapted from FESOM code by Q. Wang. 
 ! Added the driving routine forcing_setup.
 ! S.D 05.04.12
@@ -21,6 +34,7 @@ subroutine forcing_setup(partit, mesh)
 USE MOD_PARTIT
 USE MOD_PARSUP
 use forcing_array_setup_interfaces
+use forcing_array_setup_dbgyre_interfaces
 implicit none
 type(t_mesh),   intent(in),    target :: mesh
 type(t_partit), intent(inout), target :: partit
@@ -32,7 +46,34 @@ subroutine forcing_setup(partit, mesh)
      call sbc_ini(partit, mesh)         ! initialize forcing fields
 #endif
   endif 
+  if ((toy_ocean) .AND. TRIM(which_toy)=="dbgyre" .AND. (use_sw_pene)) then
+     call forcing_array_setup_dbgyre(partit, mesh)
+  endif
 end subroutine forcing_setup
+
+! ==========================================================
+subroutine forcing_array_setup_dbgyre(partit, mesh)
+  use g_forcing_arrays
+  use mod_mesh
+  USE MOD_PARTIT
+  implicit none
+  type(t_mesh),   intent(in),    target :: mesh
+  type(t_partit), intent(inout), target :: partit
+  integer    :: n2
+
+#include "associate_part_def.h"
+#include "associate_mesh_def.h"
+#include "associate_part_ass.h"
+#include "associate_mesh_ass.h"
+
+  n2=myDim_nod2D+eDim_nod2D 
+
+  allocate(chl(n2))
+  allocate(sw_3d(nl,n2))
+  chl=0.1_WP
+
+end subroutine forcing_array_setup_dbgyre
+
 ! ==========================================================
 subroutine forcing_array_setup(partit, mesh)
   !inializing forcing fields 

src/gen_modules_backscatter.F90

@@ -361,7 +361,7 @@ subroutine uke_update(dynamics, partit, mesh)
 
                 !Taking out that one place where it is always weird (Pacific Southern Ocean)
                 !Should not really be used later on, once we fix the issue with the 1/4 degree grid
-                if(.not. (TRIM(which_toy)=="soufflet")) then
+                if(.not. (TRIM(which_toy)=="soufflet") .AND. .not. (TRIM(which_toy)=="dbgyre") ) then
                 call elem_center(ed, ex, ey)
                 !a1=-104.*rad
                 !a2=-49.*rad

src/io_meandata.F90

@@ -210,6 +210,11 @@ subroutine ini_mean_io(ice, dynamics, tracers, partit, mesh)
     call def_stream(nod2D, myDim_nod2D, 'ssh',      'sea surface elevation',          'm',      dynamics%eta_n,                          io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
 CASE ('vve_5     ')
     call def_stream(nod2D, myDim_nod2D, 'vve_5',    'vertical velocity at 5th level', 'm/s',    dynamics%w(5,:),                         io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+CASE ('t_star        ')
+    call def_stream(nod2D, myDim_nod2D,'t_star'        , 'air temperature'      , 'C'  , t_star(:)       , io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+CASE ('qsr        ')
+    call def_stream(nod2D, myDim_nod2D,'qsr'        , 'solar radiation'      , 'W/s^2'  , qsr_c(:)       , io_list(i)%freq, io_list(i)%unit, io_list(i)%precision, partit, mesh)
+
 
 ! 2d ssh diagnostic variables
 CASE ('ssh_rhs    ')

src/oce_ale.F90

@@ -3114,7 +3114,7 @@ subroutine impl_vert_visc_ale(dynamics, partit, mesh)
     USE MOD_PARTIT
     USE MOD_PARSUP
     USE MOD_DYN
-    USE g_CONFIG,only: dt
+    USE g_CONFIG !,only: dt
     IMPLICIT NONE
     type(t_dyn)   , intent(inout), target :: dynamics
     type(t_partit), intent(inout), target :: partit
@@ -3239,8 +3239,16 @@ subroutine impl_vert_visc_ale(dynamics, partit, mesh)
         zinv=1.0_WP*dt/(zbar_n(nzmax-1)-zbar_n(nzmax))
         !!PS friction=-C_d*sqrt(UV(1,nlevels(elem)-1,elem)**2+ &
         !!PS             UV(2,nlevels(elem)-1,elem)**2)
-        friction=-C_d*sqrt(UV(1,nzmax-1,elem)**2+ &
-                        UV(2,nzmax-1,elem)**2)
+
+        if ((toy_ocean) .AND. (TRIM(which_toy)=="soufflet")) then
+           friction=-0.005
+        elseif ((toy_ocean) .AND. (TRIM(which_toy)=="dbgyre")) then
+           friction=-0.001
+        else
+           friction=-C_d*sqrt(UV(1,nzmax-1,elem)**2+ &
+                           UV(2,nzmax-1,elem)**2)
+        end if
+
         ur(nzmax-1)=ur(nzmax-1)+zinv*friction*UV(1,nzmax-1,elem)
         vr(nzmax-1)=vr(nzmax-1)+zinv*friction*UV(2,nzmax-1,elem)
 

src/oce_ale_pressure_bv.F90

@@ -3068,6 +3068,8 @@ SUBROUTINE density_linear(t, s, bulk_0, bulk_pz, bulk_pz2, rho_out)
 
   IF((toy_ocean) .AND. (TRIM(which_toy)=="soufflet")) THEN
       rho_out  = density_0 - 0.00025_WP*(t - 10.0_WP)*density_0
+  ELSE IF((toy_ocean) .AND. (TRIM(which_toy)=="dbgyre")) THEN
+      rho_out  = density_0 - density_0*0.0002052_WP*(t - 10.0_WP) + density_0*0.00079_WP*(s - 35.0_WP)
   ELSE
       rho_out  = density_0 + 0.8_WP*(s - 34.0_WP) - 0.2*(t - 20.0_WP)
   END IF

src/oce_ale_tracer.F90

@@ -158,6 +158,9 @@ subroutine solve_tracers_ale(ice, dynamics, tracers, partit, mesh)
     use g_comm_auto
     use o_tracers
     use Toy_Channel_Soufflet
+    use Toy_Channel_Dbgyre
+    use o_ARRAYS, only: heat_flux
+    use g_forcing_arrays, only: sw_3d
     use diff_tracers_ale_interface
     use oce_adv_tra_driver_interfaces
 #if defined(__recom)
@@ -981,12 +984,20 @@ subroutine diff_ver_part_impl_ale(tr_num, dynamics, tracers, ice, partit, mesh)
 
         !_______________________________________________________________________
         ! case of activated shortwave penetration into the ocean, ad 3d contribution
-        if (use_sw_pene .and. tracers%data(tr_num)%ID==1) then
+        if (use_sw_pene .and. tracers%data(tr_num)%ID==1 .and. .not. toy_ocean) then
             do nz=nzmin, nzmax-1
                 zinv=1.0_WP*dt  !/(zbar(nz)-zbar(nz+1)) ale!
                 !!PS tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * ( area(nz+1,n)/areasvol(nz,n)) ) * zinv
                 tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n) * area(nz+1,n)/areasvol(nz,n)) * zinv
             end do
+        elseif (use_sw_pene .and. (tracers%data(tr_num)%ID==1) .and. toy_ocean .and. TRIM(which_toy)=="dbgyre") then
+
+         call cal_shortwave_rad_dbgyre(ice, tracers, partit, mesh)
+            do nz=nzmin, nzmax-1
+                zinv=1.0_WP*dt  !/(zbar(nz)-zbar(nz+1)) ale!
+                tr(nz)=tr(nz)+(sw_3d(nz, n)-sw_3d(nz+1, n)*area(nz+1,n)/area(nz,n))*zinv
+            end do
+
         end if
 
         !_______________________________________________________________________

src/oce_modules.F90

@@ -199,6 +199,7 @@ MODULE o_ARRAYS
 real(kind=WP), allocatable    :: Ssurf_ib(:) ! kh 15.03.21 additional array for asynchronous iceberg computations
 real(kind=WP), allocatable    :: virtual_salt(:), relax_salt(:)
 real(kind=WP), allocatable    :: Tclim(:,:), Sclim(:,:)
+real(kind=WP), allocatable    :: t_star(:), qsr_c(:)
 
 !--------------
 ! LA: add iceberg tracer arrays 2023-02-08