diff --git a/config/namelist.config.toy_dbgyre b/config/namelist.config.toy_dbgyre new file mode 100644 index 000000000..3b8966e4f --- /dev/null +++ b/config/namelist.config.toy_dbgyre @@ -0,0 +1,79 @@ +! This is the namelist file for model general configuration + +&modelname +runid='fesom' +/ + +×tep +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 +/ + diff --git a/config/namelist.dyn.toy_dbgyre b/config/namelist.dyn.toy_dbgyre new file mode 100644 index 000000000..07c9a7eeb --- /dev/null +++ b/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 +/ diff --git a/config/namelist.oce.toy_dbgyre b/config/namelist.oce.toy_dbgyre new file mode 100644 index 000000000..bb6d84856 --- /dev/null +++ b/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, use C_d=0.001 for double gyre setup based on Bagaeva et al.(2024) +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 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. +/ diff --git a/env.sh b/env.sh index 6ab2d3799..e6cf2ce6b 100755 --- a/env.sh +++ b/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 diff --git a/env/nesh/shell b/env/nesh/shell new file mode 100644 index 000000000..132d85719 --- /dev/null +++ b/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 diff --git a/src/gen_forcing_init.F90 b/src/gen_forcing_init.F90 old mode 100755 new mode 100644 index 97dcdbd42..8d92be366 --- a/src/gen_forcing_init.F90 +++ b/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 diff --git a/src/gen_modules_backscatter.F90 b/src/gen_modules_backscatter.F90 index 1a3c620cd..55a6aa985 100644 --- a/src/gen_modules_backscatter.F90 +++ b/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 diff --git a/src/io_meandata.F90 b/src/io_meandata.F90 index c5f958320..9de8c783e 100644 --- a/src/io_meandata.F90 +++ b/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 ') diff --git a/src/oce_ale.F90 b/src/oce_ale.F90 index 681bfa51e..8dfd42b9c 100644 --- a/src/oce_ale.F90 +++ b/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,14 @@ 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)=="dbgyre")) then + friction=-C_d + 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) diff --git a/src/oce_ale_pressure_bv.F90 b/src/oce_ale_pressure_bv.F90 index 0a7b3c564..95ddc6e5a 100644 --- a/src/oce_ale_pressure_bv.F90 +++ b/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 diff --git a/src/oce_ale_tracer.F90 b/src/oce_ale_tracer.F90 index bb40ff335..244e629e0 100644 --- a/src/oce_ale_tracer.F90 +++ b/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 !_______________________________________________________________________ diff --git a/src/oce_modules.F90 b/src/oce_modules.F90 index 57f647e25..f2a8e22cf 100755 --- a/src/oce_modules.F90 +++ b/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 diff --git a/src/oce_setup_step.F90 b/src/oce_setup_step.F90 index 4a1e96d9d..66cdb046d 100755 --- a/src/oce_setup_step.F90 +++ b/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_Dbgyre 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("dbgyre") + call initial_state_dbgyre(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 ! ================ diff --git a/src/oce_shortwave_pene_dbgyre.F90 b/src/oce_shortwave_pene_dbgyre.F90 new file mode 100644 index 000000000..59a130642 --- /dev/null +++ b/src/oce_shortwave_pene_dbgyre.F90 @@ -0,0 +1,124 @@ +subroutine cal_shortwave_rad_dbgyre(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 + use g_comm_auto + use Toy_Channel_Dbgyre + 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_dbgyre diff --git a/src/toy_channel_dbgyre.F90 b/src/toy_channel_dbgyre.F90 new file mode 100644 index 000000000..f03a287ad --- /dev/null +++ b/src/toy_channel_dbgyre.F90 @@ -0,0 +1,261 @@ +MODULE Toy_Channel_Dbgyre + use mod_mesh + USE o_ARRAYS + USE o_PARAM + USE MOD_PARSUP + USE MOD_PARTIT + USE MOD_TRACER + USE MOD_DYN + USE g_config + + implicit none + SAVE + private + public :: initial_state_dbgyre + real(kind=WP) :: ysize =3180000.0 ! the meridional lenght of the channel [m] + real(kind=WP) :: xsize =2120000.0 ! the zonal lenght of the channel [m] !4.5*pi*r_earth=90018410.49779853 + real(kind=WP) :: zsize =4000.0 ! m The depth + real(kind=WP) :: lat0 =17.5 + +!----------------------------------------------------------------------------------------- +! + contains +! +!-------------------------------------------------------------------------------------------- + +subroutine initial_state_dbgyre(dynamics, tracers, partit, mesh) + ! Profiles NEMO 2010 + implicit none + type(t_mesh), intent(inout) , target :: mesh + type(t_partit), intent(inout), target :: partit + type(t_tracer), intent(inout), target :: tracers + type(t_dyn), intent(inout), target :: dynamics + + integer :: elem, n, nz, elnodes(3),m,k,nzmax + real(kind=8) :: dst!, yn, Fy, Lx + real(kind=WP) :: lon, lat + real(kind=WP) :: ztau, zTstar, ztrp ! wind intensity, atmospheric intensity, heat flux intensity + REAL(kind=WP) :: zrhoa ! Air density kg/m3 + REAL(kind=WP) :: a_constant, b_constant ! constants + REAL(kind=WP) :: t_star_top, t_star_bottom + integer :: wind, tpert, tprofile, sprofile, tflux, elevation + real(kind=WP):: c, c2, c3, c4, c5 + real(kind=WP):: v1, v2, sc1, sc2 + real(kind=WP):: swsurf, aux + +#include "associate_part_def.h" +#include "associate_mesh_def.h" +#include "associate_part_ass.h" +#include "associate_mesh_ass.h" + + ! default values + wind=1 + tpert=0 + tprofile=1 + sprofile=0 + tflux=1 + elevation=0 + stress_surf = 0.0_8 !make sure how it works (y-component is 0) + heat_flux = 0.0_8 + + t_star = 0.0_8 + qsr_c = 0.0_8 + relax2clim = 0.0_8 + tracers%data(2)%values(:,:) = 35.0_8 !tr_arr(:,:,2) = 35.0_8 + !tr_arr(:,:,1) = 20.0_8 + + Tsurf = tracers%data(1)%values(1,:) !tr_arr(1,:,1) + Ssurf = tracers%data(2)%values(1,:) !tr_arr(1,:,2) + water_flux = 0.0_8 + !clim_relax = 0.0_8 + + + +! ======== +! wind forcing (momentum flux): +! ======== + ! take the mean intensity as "spring" (and we don't need rotation coefficient as we don't have i,j coordinates) + ztau = -0.12 !0.105 + zrhoa = 1.22 + + + if (wind==1) then + DO elem=1, myDim_elem2D + elnodes=elem2d_nodes(:,elem) + lat=sum(coord_nod2D(2,elnodes))/3.0_WP + lon=sum(coord_nod2D(1,elnodes))/3.0_WP + + stress_surf(1,elem)= ztau * SIN( pi * (lat/rad - 15.0) / (29.0-15.0) ) + + ! write to the file (in the first time step both x and y components) + + END DO + else + stress_surf = 0.0 + endif + + +! ======== +! elevation +! ======== + if (elevation==1) then + dynamics%eta_n=0.01*(coord_nod2D(2,:)-30.0*rad)/(15.0*rad) + end if + +! ======== +! Initial temperature profile (use Philander analytic profile) +! ======== + !IF (tprofile==1) THEN + ! DO nz=1, nlevels_nod2D(n)-1 + ! tr_arr(nz,:,1) = (1+Z(nz)/4000)*(coord_nod2D(2,:)/rad) + ! END DO + !ENDIF + if (tprofile==1) then + DO n=1, myDim_nod2D+eDim_nod2D + !lat = coord_nod2D(2,n)/rad + DO nz=1, nlevels_nod2D(n)-1 + + !tr_arr(nz,n,1) = (1+Z(nz)/4000)*(coord_nod2D(2,n)/rad) + !tr_arr(nz,n,1)=tr_arr(nz,n,1)-0.95_WP*20*tanh(abs(Z(nz))/300)-abs(Z(nz))/2400.0_WP + tracers%data(1)%values(nz,n) = 7.5 * ( 1.-TANH((ABS(Z(nz))-80.)/30.) ) & + & + 10. * ( 5000. - ABS(Z(nz)) ) /5000. + !tr_arr(nz,n,1)=25.0 + .05*tanh((Z(nz) - 300.0)/5.0) - Z(nz)*0.008_WP + + END DO + END DO + ENDIF + !write(*,*) mype, 'Temperature', maxval(tr_arr(:,:,1)), minval(tr_arr(:,:,1)) + + + Tsurf=tracers%data(1)%values(1,:) !tr_arr(1,:,1) + + + if (sprofile==1) then + DO n=1, myDim_nod2D+eDim_nod2D + + DO nz=1, nlevels_nod2D(n)-1 + + tracers%data(2)%values(nz,n) = 36. - ABS(Z(nz)) /4000. + + + END DO + END DO + ENDIF + !write(*,*) mype, 'Salinity', maxval(tr_arr(:,:,2)), minval(tr_arr(:,:,2)) + + !Ssurf=tr_arr(1,:,2) + +! ======== +! temperature flux from atmosphere +! ======== + zTstar = 28.3 ! intensity from 28.3 a -5 deg + ztrp = 4.0 ! retroaction term on heat fluxes (W/m2/K) + + IF (tflux==1) then + DO n=1, myDim_nod2D+eDim_nod2D + lat = coord_nod2D(2,n) + !the function varies from 27C at 50N to 7C at 15N, again the analogy with spring + t_star(n) = zTstar * cos(pi*((lat / rad - 5.0 ) / 107.0)) + !solar radiation component + qsr_c(n) = 230 * COS( pi * ((lat / rad) / 162.0 )) + + heat_flux(n) = ztrp * (tracers%data(1)%values(1,n) - t_star(n)) !+ qsr_c(n) + + !write(*,*) mype, 't_star', t_star + + ! shortwave rad. + !swsurf=(1.0_WP-albw)*qsr_c(n) + ! the visible part (300nm-750nm) + !swsurf=swsurf*0.54_WP + ! subtract visible sw rad. from heat_flux, which is '+' for upward + !heat_flux(n)=heat_flux(n)+swsurf + + ! limit chl from below + !if (chl(n) < 0.02_WP) chl(n)=0.02_WP + + + !c=log10(chl(n)) + !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=(nlevels(n)) + !sw_3d(1, n)=swsurf + !do k=2, nzmax + ! aux=(v1*exp(zbar_3d_n(k,n)/sc1)+v2*exp(zbar_3d_n(k,n)/sc2)) + ! sw_3d(k, n)=swsurf*aux + ! if (aux < 1.e-5_WP .OR. k==nzmax) then + ! sw_3d(k, n)=0.0_WP + ! exit + ! end if + !end do + + + END DO + ENDIF + + + + !if (mype==0) write(*,*) "heat_flux_channel", heat_flux(100) + !if (mype==0) write(*,*) "radiation_channel", qsr_c(100) + + + ! ========= + ! Coriolis calculation + ! ========= + + DO n=1,myDim_elem2D + elnodes=elem2D_nodes(:,n) + !dst - change in latitude + dst=(sum(coord_nod2D(2, elnodes))/3.0-lat0*rad)*r_earth-ysize/2 + mesh%coriolis(n)=1.0e-4+dst*1.8e-11 + END DO + write(*,*) mype, 'COR', maxval(mesh%coriolis*10000.0), minval(mesh%coriolis*10000.0) + + DO n=1,myDim_nod2D+eDim_nod2D + dst=(coord_nod2D(2, n)-lat0*rad)*r_earth-ysize/2 + mesh%coriolis_node(n)=1.0e-4+dst*1.8e-11 + END DO + write(*,*) mype, 'COR_n', maxval(mesh%coriolis_node*10000.0), minval(mesh%coriolis_node*10000.0) + + + + ! ========= + ! Temperature perturbation to check + ! ========= + +if (tpert==1) then + do n=1, myDim_nod2D+eDim_nod2D + lat=coord_nod2D(2,n)-15.0*rad + lon=coord_nod2D(1,n) + if (lon>cyclic_length/2) lon=lon-cyclic_length + if (lon<-cyclic_length/2) lon=lon+cyclic_length + dst=sqrt((lat)**2+(lon)**2) + if (dst>1.5*rad) cycle + do nz=1, nlevels_nod2D(n)-1 + !if(abs(Z(nz)+500)<300) then + tracers%data(1)%values(nz,n)=tracers%data(1)%values(nz,n)+1.0*cos(pi*dst/2.0/1.5/rad) !exp(-(dst/(1.5*rad))**2) + !end if + end do + end do +end if +!write(*,*) mype, 'dst', dst +!write(*,*) mype, 'rad', rad + +end subroutine initial_state_dbgyre +END MODULE Toy_Channel_Dbgyre diff --git a/work/job_nesh b/work/job_nesh new file mode 100755 index 000000000..ca718f3e0 --- /dev/null +++ b/work/job_nesh @@ -0,0 +1,29 @@ +#!/bin/bash +#SBATCH --job-name=nemo_dg_test +#SBATCH --partition=base +#SBATCH --time=02:45:00 +#SBATCH --nodes=1 +#SBATCH --tasks-per-node=20 +#SBATCH --cpus-per-task=1 +#SBATCH --mem=60G + +#SBATCH --output %x_%j.out +#SBATCH --error %x_%j.err + +# disable hyperthreading +#SBATCH --hint=nomultithread + +module purge +source ../env/nesh/shell +export OMP_NUM_THREADS=1 + +ln -sf ../bin/fesom.x . + +#___DETERMINE SLURM JOBID+OUTPUTFILE____________________________________________ +jobid=$(echo $SLURM_JOB_ID | cut -d"." -f1) +fname="${SLURM_JOB_NAME}_${jobid}.out" + +#___PUT JOB IN QUEUE____________________________________________________________ +date +srun --mpi=pmi2 ./fesom.x >> ${fname} +date