Only necessary changes for Double Gyre setup

Author: alexwilms

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_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)=="nemo") ) 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)=="nemo")) 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)=="nemo")) 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_Nemo
+    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)=="nemo") then
+
+         call cal_shortwave_rad_nemo(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

src/oce_setup_step.F90

@@ -91,6 +91,7 @@ subroutine ocean_setup(dynamics, tracers, partit, mesh)
     use g_cvmix_tidal
     use g_backscatter
     use Toy_Channel_Soufflet
+    use Toy_Channel_Nemo
     use oce_initial_state_interface
     use oce_adv_tra_fct_interfaces
     use init_ale_interface
@@ -226,6 +227,8 @@ subroutine ocean_setup(dynamics, tracers, partit, mesh)
               call compute_zonal_mean_ini(partit, mesh)  
               call compute_zonal_mean(dynamics, tracers, partit, mesh)
            end if
+         CASE("nemo")
+             call initial_state_nemo(dynamics, tracers, partit, mesh)
        END SELECT
     else
        if (flag_debug .and. partit%mype==0)  print *, achar(27)//'[36m'//'     --> call oce_initial_state'//achar(27)//'[0m' 
@@ -748,6 +751,8 @@ SUBROUTINE arrays_init(num_tracers, partit, mesh)
     ! ================     
     !allocate(stress_diag(2, elem_size))!delete me
     !!PS allocate(Visc(nl-1, elem_size))
+    allocate(t_star(node_size))
+    allocate(qsr_c(node_size))
     ! ================
     ! elevation and its rhs
     ! ================

src/oce_shortwave_pene_nemo.F90

@@ -0,0 +1,127 @@
+subroutine cal_shortwave_rad_nemo(ice, tracers, partit, mesh)
+  ! This routine is inherited from FESOM 1.4 and adopted appropreately. It calculates 
+  ! shortwave penetration into the ocean assuming the constant chlorophyll concentration.
+  ! No penetration under the ice is applied. A decent way for ice region is to be discussed.
+  ! This routine should be called after ice2oce coupling done if ice model is used.
+  ! Ref.: Morel and Antoine 1994, Sweeney et al. 2005
+  USE MOD_MESH
+  USE MOD_ICE
+  USE o_PARAM
+  USE o_ARRAYS
+  USE MOD_PARSUP
+  USE MOD_PARTIT
+  USE MOD_TRACER
+  USE g_CONFIG
+  use g_forcing_arrays!, only: chl, sw_3d
+  use g_comm_auto
+ !  use i_param
+ !  use i_arrays
+ ! use i_therm_param
+ ! use Toy_Channel_Nemo!, only: qsr_c, t_star
+  IMPLICIT NONE
+  type(t_ice) ,intent(in), target :: ice
+  type(t_tracer), intent(inout), target :: tracers
+  type(t_partit), intent(inout), target :: partit
+  type(t_mesh), intent(in) , target :: mesh
+
+
+  integer      :: m, n2, n3, k, nzmax, nzmin
+  real(kind=WP):: swsurf, aux, zTstar, ztrp
+  real(kind=WP):: c, c2, c3, c4, c5
+  real(kind=WP):: v1, v2, sc1, sc2
+  real(kind=WP), pointer :: albw
+
+
+#include "associate_part_def.h"
+#include "associate_mesh_def.h"
+#include "associate_part_ass.h"
+#include "associate_mesh_ass.h"
+
+  sw_3d=0.0_WP
+  albw => ice%thermo%albw  
+
+  
+  do n2=1, myDim_nod2D+eDim_nod2D     
+     
+     !calculate heat flux
+     zTstar    = 28.3      ! intensity from 28.3 a -5 deg
+     ztrp = 4.0        ! retroaction term on heat fluxes (W/m2/K)
+     
+     !lat = coord_nod2D(2,n2)
+     
+     !t_star (n2) = zTstar * cos(pi*((lat / rad - 5.0 ) / 107.0))
+     heat_flux(n2) = ztrp * (tracers%data(1)%values(1,n2) - t_star(n2)) 
+     
+     
+     
+     ! shortwave rad.
+     swsurf=(1.0_WP-albw)*qsr_c(n2)
+     ! the visible part (300nm-750nm)
+     swsurf=swsurf*0.54_WP
+     ! subtract visible sw rad. from heat_flux, which is '+' for upward
+     
+     !if (mype==0) write(*,*) "heat_flux_routine", heat_flux(100)
+     !if (mype==0) write(*,*) "radiation_routine", qsr_c(100)
+     
+     heat_flux(n2)=heat_flux(n2)+swsurf
+ 
+     ! attenuation func. for vis. sw rad. according to Morel/Antoine param.
+     ! the four parameters in the func.
+     
+
+
+     ! limit chl from below
+     if (chl(n2) < 0.02_WP) chl(n2)=0.02_WP
+     
+     
+     c=log10(chl(n2))
+     c2=c*c
+     c3=c2*c
+     c4=c3*c
+     c5=c4*c
+     ! --> coefficients come from Sweeney et al. 2005, "Impacts of shortwave
+     ! penetration depthon large scale ocean circulation and heat transport" see 
+     ! Appendix A
+     
+
+     
+     v1=0.008_WP*c+0.132_WP*c2+0.038_WP*c3-0.017_WP*c4-0.007_WP*c5
+     v2=0.679_WP-v1
+     v1=0.321_WP+v1
+     sc1=1.54_WP-0.197_WP*c+0.166_WP*c2-0.252_WP*c3-0.055_WP*c4+0.042_WP*c5
+     sc2=7.925_WP-6.644_WP*c+3.662_WP*c2-1.815_WP*c3-0.218_WP*c4+0.502_WP*c5
+    
+
+    
+     ! convert from heat flux [W/m2] to temperature flux [K m/s]
+     swsurf=swsurf/vcpw
+     ! vis. sw. rad. in the colume
+     nzmax=(mesh%nlevels(n2))
+     nzmin=(mesh%ulevels(n2))
+     sw_3d(nzmin, n2)=swsurf
+     do k=nzmin+1, nzmax
+        aux=(v1*exp(mesh%zbar_3d_n(k,n2)/sc1)+v2*exp(mesh%zbar_3d_n(k,n2)/sc2))
+        sw_3d(k, n2)=swsurf*aux
+        if (aux < 1.e-5_WP .OR. k==nzmax) then 
+           sw_3d(k, n2)=0.0_WP
+           exit
+        end if
+     end do
+     
+     ! sw_3d --> TEMPERATURE FLUX through full depth level interfaces into/out off 
+     ! the tracer volume 
+     ! sum(sw_3d(1:nlevels(n2)-1,n2)-sw_3d(2:nlevels(n2),n2)) = swsurf !!!
+     
+!for testing the subroutine
+!if (mype==30 .and. n2==100) then
+!write(*,*) 'heat_flux=', heat_flux(n2)
+!write(*,*) 'short/longwave=', shortwave(n2), longwave(n2), swsurf*vcpw
+!do k=1, nzmax
+!   write(*,*) sw_3d(k, n2)*vcpw
+!end do
+!end if
+
+  end do
+!call par_ex
+!stop
+end subroutine cal_shortwave_rad_nemo