first running version without correct restart

CMakeLists.txt

@@ -12,7 +12,7 @@ set(TOPLEVEL_DIR ${CMAKE_CURRENT_LIST_DIR})
 set(FESOM_COUPLED OFF CACHE BOOL "compile fesom standalone or with oasis support (i.e. coupled)")
 set(OIFS_COUPLED OFF CACHE BOOL "compile fesom coupled to OpenIFS. (Also needs FESOM_COUPLED to work)")
 set(CRAY OFF CACHE BOOL "compile with cray ftn")
-set(USE_ICEPACK OFF CACHE BOOL "compile fesom with the Iceapck modules for sea ice column physics.")
+set(USE_ICEPACK ON CACHE BOOL "compile fesom with the Icepack modules for sea ice column physics.")
 set(OPENMP_REPRODUCIBLE OFF CACHE BOOL "serialize OpenMP loops that are critical for reproducible results")
 
 #set(VERBOSE OFF CACHE BOOL "toggle debug output")

config/namelist.config

@@ -5,31 +5,31 @@ runid='fesom'
 /
 
 &timestep
-step_per_day=32 !96 !96 !72 !72 !45 !72 !96
-run_length=62 !62 !62 !62 !28
+step_per_day=96 !96 !96 !72 !72 !45 !72 !96
+run_length=20 !62 !62 !62 !28
 run_length_unit='y'		! y, m, d, s
 /
 
 &clockinit              ! the model starts at
 timenew=0.0
 daynew=1
-yearnew=1958
+yearnew=1960
 /
 
 &paths
-MeshPath='/work/ollie/projects/clidyn/FESOM2/meshes/core2/'
-ClimateDataPath='/work/ollie/projects/clidyn/FESOM2/hydrography/phc3.0/'
-ResultPath='../result_tmp/'
+MeshPath='/albedo/work/user/pscholz/mesh_fesom2.0/core2_srt_dep@node/'
+ClimateDataPath='/albedo/work/projects/p_fesom/FROM-OLLIE/FESOM2/hydrography/phc3.0/'
+ResultPath='/albedo/scratch/user/fkauker/fesom2/core2_head-ip/'
 /
 
 &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
+restart_length_unit='m'     !output period: y, d, h, s, off
 raw_restart_length=1        ! --> do core dump restart
-raw_restart_length_unit='y' ! e.g. y, d, h, s, off
+raw_restart_length_unit='m' ! e.g. y, d, h, s, off
 bin_restart_length=1        ! --> do derived type binary restart 
-bin_restart_length_unit='y' ! e.g. y, d, h, s, off
-logfile_outfreq=960     !in logfile info. output frequency, # steps
+bin_restart_length_unit='m' ! e.g. y, d, h, s, off
+logfile_outfreq=96          !in logfile info. output frequency, # steps
 /
 
 &ale_def
@@ -62,6 +62,7 @@ flag_debug=.false.
 /
 
 &machine
-n_levels=2
-n_part= 12, 36          ! 432 number of partitions on each hierarchy level
+n_levels=1
+n_part= 508          ! 508 number of partitions on each hierarchy level
 /
+

config/namelist.cvmix

@@ -26,7 +26,9 @@ idemix_mu0         = 0.33333333 ! dissipation parameter (dimensionless)
 idemix_sforcusage  = 0.2
 idemix_n_hor_iwe_prop_iter = 5  ! iterations for contribution from horiz. wave propagation
 idemix_surforc_file = '/work/ollie/clidyn/forcing/IDEMIX/fourier_smooth_2005_cfsr_inert_rgrid.nc'
+idemix_surforc_vname= 'var706'
 idemix_botforc_file = '/work/ollie/clidyn/forcing/IDEMIX/tidal_energy_gx1v6_20090205_rgrid.nc'
+idemix_botforc_vname= 'wave_dissipation'
 /
 
 ! namelist for PP

config/namelist.dyn

@@ -4,7 +4,8 @@ visc_gamma1  = 0.1   ! [nodim], for computation of the flow aware viscosity
 visc_gamma2  = 0.285 ! [s/m],   is only used in easy backscatter option
 visc_easybsreturn= 1.5
 
-opt_visc     = 5          
+opt_visc     = 5     
+check_opt_visc=.true.  ! 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)
@@ -20,5 +21,6 @@ use_wsplit   = .false. ! Switch for implicite/explicte splitting of vert. veloci
 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.icepack

@@ -110,12 +110,16 @@
 !!!!!!!!!!                Icepack output namelist                    !!!!!!!!!!
 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
 
+&nml_listsize
+    io_listsize	    = 100
+/
+
 &nml_list_icepack
 io_list_icepack =  'aicen     ',1, 'm', 4, ! Sea ice concentration
                    'vicen     ',1, 'm', 4, ! Volume per unit area of ice
                    'vsnon     ',1, 'm', 4, ! Volume per unit area of snow
-                   !'Tsfcn     ',1, 'm', 4, ! Sea ice surf. temperature
-                   !'iagen     ',1, 'm', 4, ! Sea ice age
+                   'Tsfcn     ',1, 'm', 4, ! Sea ice surf. temperature
+                   'iagen     ',1, 'm', 4, ! Sea ice age
                    !'FYn       ',1, 'm', 4, ! First year ice
                    !'lvln      ',1, 'm', 4, ! Ridged sea ice area and volume
                    !'pond_cesmn',1, 'm', 4, ! Melt ponds area, volume and refrozen lid thickness
@@ -129,10 +133,10 @@ io_list_icepack =  'aicen     ',1, 'm', 4, ! Sea ice concentration
                    !'aice      ',1, 'm', 4, ! Sea ice concentration
                    !'vice      ',1, 'm', 4, ! Volume per unit area of ice
                    !'vsno      ',1, 'm', 4, ! Volume per unit area of snow
-                   !'Tsfc      ',1, 'm', 4, ! Sea ice surf. temperature
-                   !'iage      ',1, 'm', 4, ! Sea ice age
-                   !'FY        ',1, 'm', 4, ! First year ice
-                   !'lvl       ',1, 'm', 4, ! Ridged sea ice area and volume
+                   'Tsfc      ',1, 'm', 4, ! Sea ice surf. temperature
+                   'iage      ',1, 'm', 4, ! Sea ice age
+                   'FY        ',1, 'm', 4, ! First year ice
+                   'lvl       ',1, 'm', 4, ! Ridged sea ice area and volume
                    !'pond_cesm ',1, 'm', 4, ! Melt ponds area, volume and refrozen lid thickness
                    !'pond_topo ',1, 'm', 4, ! Melt ponds area, volume and refrozen lid thickness
                    !'pond_lvl  ',1, 'm', 4, ! Melt ponds area, volume and refrozen lid thickness

config/namelist.io

@@ -12,7 +12,7 @@ ldiag_extflds    =.false.
 /
 
 &nml_general
-io_listsize    =100 !number of streams to allocate. shallbe large or equal to the number of streams in &nml_list
+io_listsize    =120 !number of streams to allocate. shallbe large or equal to the number of streams in &nml_list
 vec_autorotate =.false.
 /
 
@@ -21,22 +21,25 @@ vec_autorotate =.false.
 ! for 'fer_C', 'bolus_u', 'bolus_v', 'bolus_w', 'fer_K' to work Fer_GM must be .true. otherwise no output
 ! 'otracers' - all other tracers if applicable
 ! for 'dMOC' to work ldiag_dMOC must be .true. otherwise no output
+! for 'utemp', 'vtemp', 'usalt', 'vsalt' output, set ldiag_trflx=.true.
 &nml_list
 io_list =  'sst       ',1, 'm', 4,
            'sss       ',1, 'm', 4,
     	   'ssh       ',1, 'm', 4,
-           'uice      ',1, 'd', 4,
-           'vice      ',1, 'd', 4,
+           'uice      ',1, 'm', 4,
+           'vice      ',1, 'm', 4,
            'a_ice     ',1, 'm', 4,
            'm_ice     ',1, 'm', 4,
            'm_snow    ',1, 'm', 4,
-           'MLD1      ',1, 'm', 4,
+	   'flice     ',1, 'm', 4,
+	   'MLD1      ',1, 'm', 4,
            'MLD2      ',1, 'm', 4,
            'MLD3      ',1, 'm', 4,
            'tx_sur    ',1, 'm', 4,
            'ty_sur    ',1, 'm', 4,
-           'temp      ',1, 'y', 4,
-           'salt      ',1, 'y', 8,
+           'temp      ',1, 'm', 4,
+           'salt      ',1, 'm', 4,
+	   'otracers  ',1, 'm', 4,
            'N2        ',1, 'y', 4,
            'Kv        ',1, 'y', 4,
            'u         ',1, 'y', 4,
@@ -48,4 +51,6 @@ io_list =  'sst       ',1, 'm', 4,
            'bolus_u   ',1, 'y', 4,
            'bolus_v   ',1, 'y', 4,
            'bolus_w   ',1, 'y', 4,
+           'fw        ',1, 'm', 4,
+           'fh        ',1, 'm', 4,
 /

config/namelist.tra

@@ -4,9 +4,8 @@ num_tracers=100 !number of tracers to allocate. shallbe large or equal to the nu
 
 &tracer_list
 nml_tracer_list =  
-1  , 'MFCT',       'QR4C',       'FCT ',        1.,          1., 
-2  , 'MFCT',       'QR4C',       'FCT ',        1.,          1., 
-!101, 'UPW1',       'UPW1',       'NON ',        0.,          0.
+1  , 'MFCT',       'QR4C',       'FCT ',        0.,          1., 
+2  , 'MFCT',       'QR4C',       'FCT ',        0.,          1., 
 /
 
 &tracer_init3d                            ! initial conditions for tracers

src/CMakeLists.txt

@@ -152,6 +152,7 @@ if(${CMAKE_Fortran_COMPILER_ID} STREQUAL  Intel )
       target_compile_options(${PROJECT_NAME} PRIVATE -march=core-avx2 -mtune=core-avx2)
    elseif(${FESOM_PLATFORM_STRATEGY} STREQUAL  albedo)
       target_compile_options(${PROJECT_NAME} PRIVATE -march=core-avx2 -O3 -ip -fPIC -qopt-malloc-options=2 -qopt-prefetch=5 -unroll-aggressive) #NEC mpi option
+      # target_compile_options(${PROJECT_NAME} PRIVATE -march=core-avx2 -O0 -g -traceback -check bounds ) # Franks debug options
    else()
       target_compile_options(${PROJECT_NAME} PRIVATE -xHost)
    endif()

src/ice_EVP.F90

@@ -428,11 +428,11 @@ subroutine EVPdynamics(ice, partit, mesh)
 
                 !_______________________________________________________________
                 ! Hunke and Dukowicz c*h*p*
-#if defined (__icepack)
-                ice_strength(el) = ice%pstar*msum*exp(-ice%c_pressure*(1.0_WP-asum))
-#else
+! #if defined (__icepack) # [email protected]
                 ice_strength(el) = ice%pstar*msum*exp(-ice%c_pressure*(1.0_WP-asum))
-#endif
+! #else
+!                ice_strength(el) = ice%pstar*msum*exp(-ice%c_pressure*(1.0_WP-asum))
+! #endif # [email protected]
                 ice_strength(el) = 0.5_WP*ice_strength(el)
 
                 !_______________________________________________________________
@@ -481,11 +481,11 @@ subroutine EVPdynamics(ice, partit, mesh)
                 asum = sum(a_ice(elnodes))/3.0_WP
 
                 ! ===== Hunke and Dukowicz c*h*p*
-#if defined (__icepack)
-                ice_strength(el) = ice%pstar*msum*exp(-ice%c_pressure*(1.0_WP-asum))
-#else
+! #if defined (__icepack) # [email protected]
                 ice_strength(el) = ice%pstar*msum*exp(-ice%c_pressure*(1.0_WP-asum))
-#endif
+! #else
+!                ice_strength(el) = ice%pstar*msum*exp(-ice%c_pressure*(1.0_WP-asum))
+! #endif # [email protected]
                 ice_strength(el) = 0.5_WP*ice_strength(el)
 
                 ! use rhs_m and rhs_a for storing the contribution from elevation:

src/ice_maEVP.F90

@@ -159,7 +159,7 @@ subroutine stress_tensor_m(ice, partit, mesh)
         eps2=eps11(elem)-eps22(elem)   
 
         ! ====== moduli:
-        delta=eps1**2+vale*(eps2**2+4.0_WP*eps12(elem)**2)
+        delta=eps1**2 + vale * ( eps2**2 + 4.0_WP*eps12(elem)**2 )
         delta=sqrt(delta)
 
 #if defined (__icepack)

src/ice_oce_coupling.F90

@@ -322,18 +322,35 @@ subroutine oce_fluxes(ice, dynamics, tracers, partit, mesh)
                            vsno_out=m_snow,                &
                            fhocn_tot_out=net_heat_flux,    &
                            fresh_tot_out=fresh_wa_flux,    &
+                           dhs_dt_out=thdgrsn,             &
+                           dhi_dt_out=thdgr,               &
                            fsalt_out=real_salt_flux,       &
-                           dhi_dt_out=thdgrsn,             &
-                           dhs_dt_out=thdgr,               &
                            evap_ocn_out=evaporation        )
 
-    heat_flux(:)   = - net_heat_flux(:)
-    water_flux(:)  = - (fresh_wa_flux(:)/1000.0_WP) - runoff(:)
+    ! write(*,*) "************************"
+    ! write(*,*) "min/max a_ice = ", minval(a_ice), maxval(a_ice)
+    ! write(*,*) "min/max net_heat_flux = ", minval(net_heat_flux), maxval(net_heat_flux)
+    ! write(*,*) "min/max fresh_wa_flux = ", minval(fresh_wa_flux), maxval(fresh_wa_flux)
+    ! write(*,*) "min/max real_salt_flux = ", minval(real_salt_flux), maxval(real_salt_flux)
+    ! write(*,*) "min/max evaporation = ", minval(evaporation), maxval(evaporation)
+    ! write(*,*) "-------------------------"
+    ! [email protected]
 
-    ! Evaporation
-    evaporation(:) = - evaporation(:) / 1000.0_WP
-    ice_sublimation(:) = 0.0_WP
+!$OMP PARALLEL DO
+    do n=1, myDim_nod2d+eDim_nod2d  
+       heat_flux(n)   = - net_heat_flux(n)
+       water_flux(n)  = - (fresh_wa_flux(n)/1000.0_WP) - runoff(n)
+       ! Evaporation
+       evaporation(n) = - evaporation(n) / 1000.0_WP
+       ice_sublimation(n) = 0.0_WP
+    end do
+!$OMP END PARALLEL DO
 
+    ! write(*,*) "min/max heat_flux = ", minval(heat_flux), maxval(heat_flux)
+    ! write(*,*) "min/max water_flux = ", minval(water_flux), maxval(water_flux)
+    ! write(*,*) "min/max evaporation = ", minval(evaporation), maxval(evaporation)
+    ! write(*,*) "************************"
+
     call init_flux_atm_ocn()
 
 #else

src/icepack_drivers/associate_mesh.h

@@ -7,8 +7,6 @@ real(kind=WP)   , pointer :: ocean_areawithcav
 integer         , pointer :: nl
 integer         , pointer :: nn_size
 
-
-
 real(kind=WP), dimension(:,:), pointer :: coord_nod2D, geo_coord_nod2D
 integer, dimension(:,:)      , pointer :: elem2D_nodes
 integer, dimension(:,:)      , pointer :: edges       
@@ -62,7 +60,7 @@ ocean_areawithcav  => mesh%ocean_areawithcav
 nl                 => mesh%nl  
 nn_size            => mesh%nn_size
 
-
+! needed to specify dimensions
 myDim_nod2D  => p_partit%myDim_nod2D
 eDim_nod2D   => p_partit%eDim_nod2D
 myDim_elem2D => p_partit%myDim_elem2D
@@ -71,7 +69,6 @@ eXDim_elem2D => p_partit%eXDim_elem2D
 myDim_edge2D => p_partit%myDim_edge2D
 eDim_edge2D  => p_partit%eDim_edge2D
 
-
 coord_nod2D(1:2,1:myDim_nod2D+eDim_nod2D)                  => mesh%coord_nod2D        
 geo_coord_nod2D(1:2,1:myDim_nod2D+eDim_nod2D)              => mesh%geo_coord_nod2D    
 elem2D_nodes(1:3, 1:myDim_elem2D+eDim_elem2D+eXDim_elem2D) => mesh%elem2D_nodes
@@ -84,17 +81,17 @@ edge_cross_dxdy(1:4,1:myDim_edge2D+eDim_edge2D)            => mesh%edge_cross_dx
 elem_cos(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D)          => mesh%elem_cos           
 metric_factor(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D)     => mesh%metric_factor      
 elem_neighbors(1:3,1:myDim_elem2D)                         => mesh%elem_neighbors     
-nod_in_elem2D      => mesh%nod_in_elem2D   ! (maxval(rmax),myDim_nod2D+eDim_nod2D)    
-x_corners          => mesh%x_corners   ! (myDim_nod2D, maxval(rmax)) 
-y_corners          => mesh%y_corners   ! (myDim_nod2D, maxval(rmax))       
+nod_in_elem2D                                              => mesh%nod_in_elem2D   ! (maxval(rmax),myDim_nod2D+eDim_nod2D)    
+x_corners                                                  => mesh%x_corners   ! (myDim_nod2D, maxval(rmax)) 
+y_corners                                                  => mesh%y_corners   ! (myDim_nod2D, maxval(rmax))       
 nod_in_elem2D_num(1:myDim_nod2D+eDim_nod2D)                => mesh%nod_in_elem2D_num  
 depth(1:myDim_nod2D+eDim_nod2D)                            => mesh%depth              
 gradient_vec(1:6,1:myDim_elem2D)                           => mesh%gradient_vec       
 gradient_sca(1:6,1:myDim_elem2D)                           => mesh%gradient_sca       
 bc_index_nod2D(1:myDim_nod2D+eDim_nod2D)                   => mesh%bc_index_nod2D     
 zbar(1:mesh%nl)                                            => mesh%zbar               
 Z(1:mesh%nl-1)                                             => mesh%Z
-elem_depth         => mesh%elem_depth      ! never used, not even allocated
+elem_depth                                                 => mesh%elem_depth      ! never used, not even allocated
 nlevels(1:myDim_elem2D+eDim_elem2D+eXDim_elem2D)           => mesh%nlevels            
 nlevels_nod2D(1:myDim_nod2D+eDim_nod2D)                    => mesh%nlevels_nod2D
 nlevels_nod2D_min(1:myDim_nod2D+eDim_nod2D)                => mesh%nlevels_nod2D_min
@@ -130,5 +127,3 @@ zbar_n_bot(1:myDim_nod2D+eDim_nod2D)                       => mesh%zbar_n_bot
 zbar_e_bot(1:myDim_elem2D+eDim_elem2D)                     => mesh%zbar_e_bot
 zbar_n_srf(1:myDim_nod2D+eDim_nod2D)                       => mesh%zbar_n_srf
 zbar_e_srf(1:myDim_elem2D+eDim_elem2D)                     => mesh%zbar_e_srf
-
-