MODULE ocean_slab_h
!
!==================================================================
!
!     Purpose
!     -------
!     Constants and runtime parameters shared by the slab ocean of the
!     Generic PCM ("ocean_slab_mod", "slab_heat_transp_mod") and by some
!     other subroutines ("hydrol", "physiq_mod", "photolysis_online").
!
!     Authors
!     -------
!     Refactor (JB Clement, 2026)
!     Constants gathered from "ocean_slab_mod" (S. Bhatnagar and E. Millour,
!     2023) and "slab_heat_transp_mod".
!
!     Notes
!     -----
!     The runtime parameters below are read from "callphys.def" either:
!         - in "ocean_slab_init", for the ones which only make sense
!           when the slab ocean is activated;
!         - in "inifis", for the ones which are also used by "hydrol" (lakes)
!           and "photolysis_online".
!
!==================================================================

use watercommon_h, only: T_h2o_ice_liq, RLFTT, rhowaterice

implicit none

!******************************************************************
! Material constants
!******************************************************************
real, parameter :: t_melt = T_h2o_ice_liq  ! Melting ice temp [K]
real, parameter :: ice_den = rhowaterice   ! Ice density [kg/m3]
real, parameter :: ice_lat = RLFTT         ! Freeze / melt latent heat of snow and ice [J/kg]
real, parameter :: sno_den = 300.          ! Mean snow density [kg/m3]
real, parameter :: sea_den = 1026.         ! Sea water density [kg/m3]
real, parameter :: ice_cond = 2.17*ice_den ! Conductivity of ice [W/(m.K) or (W.kg)/(K.m4)]
real, parameter :: ice_cap = 2067.         ! Specific heat capacity, snow and ice [J/(kg.K)]
real, parameter :: sea_cap = 3994.         ! Specific heat capacity, seawater [J/(kg.K)]
real, parameter :: snow_min = 0.05*sno_den ! Critical snow height [kg/m2]
real, parameter :: epsfra = 1.0E-05        ! Minimal grid fraction size below which there is no ice

! Thermally active thickness [m] behind the surface heat capacity of sea ice and the snow lying on it
! The volumetric capacities used to be unrelated hardcoded values:
!     capcalseaice = 5.1444e+06*0.15 = 771660 [J.K-1.m-2]
!     capcalsno    = 2.3867e+06*0.15 = 358005 [J.K-1.m-2]
! These parameters reproduce them exactly because using 'h_ice_thin' would
! divide capcalseaice by 2.7 and capcalsno by 3.85 and change how fast the
! surface temperature of sea ice and of the snow on it responds.
! To be decided to keep or correct!!!
real, parameter :: h_capcal_ice = 0.4057866 ! 771660./(ice_den*ice_cap)
real, parameter :: h_capcal_sno = 0.5773343 ! 358005./(sno_den*ice_cap)

! Mixed-layer depth [m] for which coef_hdiff was calibrated. The horizontal
! diffusion tendency is applied uniformly over the whole slab column, so it is
! rescaled by h_hdiff_ref/sum(slabh) to keep the transported heat, rather than
! the tendency, independent of the slab configuration.
real, parameter :: h_hdiff_ref = 50.

!******************************************************************
! Slab geometry, set in ocean_slab_init
!******************************************************************
! Number of slab vertical layers
integer, save :: nslay = 2
!$OMP THREADPRIVATE(nslay)
! Depth of the slab layers [m] (1st = wind-mixed layer, 2nd = deep layer)
real, allocatable, dimension(:), save :: slabh
!$OMP THREADPRIVATE(slabh)
! cyang = 1/heat capacity of top layer (rho.c.H) = 1/capcalocean
real, save :: cyang
!$OMP THREADPRIVATE(cyang)
! Surface heat capacities seen by the physics [J.K-1.m-2]
real, save :: capcalocean, capcalseaice, capcalsno
!$OMP THREADPRIVATE(capcalocean,capcalseaice,capcalsno)

!******************************************************************
! Sea ice parameters, set in ocean_slab_init
!******************************************************************
real, save :: t_freeze ! Freezing sea water temp [K] = t_melt + Tsaldiff
!$OMP THREADPRIVATE(t_freeze)
real, save :: sno_cond ! Conductivity of snow [W/(m.K) or (W.kg)/(K.m4)]
!$OMP THREADPRIVATE(sno_cond)
real, save :: snow_wfact ! Max fraction of falling snow blown into the ocean
!$OMP THREADPRIVATE(snow_wfact)
real, save :: ice_frac_min ! Min ice fraction below which the ice is discarded
!$OMP THREADPRIVATE(ice_frac_min)
real, save :: ice_frac_max ! Max ice fraction (leads)
!$OMP THREADPRIVATE(ice_frac_max)
real, save :: h_ice_min ! Min ice thickness [kg/m2]
!$OMP THREADPRIVATE(h_ice_min)
! Below h_ice_thin, priority is to melt lateral/grow height.
! h_ice_thin is also the height of new ice.
real, save :: h_ice_thin ! Thin ice thickness [kg/m2]
!$OMP THREADPRIVATE(h_ice_thin)
! Above h_ice_thick, priority is to melt height/grow lateral
real, save :: h_ice_thick ! thick ice thickness [kg/m2]
!$OMP THREADPRIVATE(h_ice_thick)
real, save :: h_ice_new ! Max height of new open ocean ice [kg/m2]
!$OMP THREADPRIVATE(h_ice_new)
real, save :: h_ice_ref ! Reference ice height for the growth partition between lateral extension and vertical thickening [kg/m2]
!$OMP THREADPRIVATE(h_ice_ref)
real, save :: h_ice_warn ! Ice height above which a warning is issued (diagnostic only) [kg/m2]
!$OMP THREADPRIVATE(h_ice_warn)

!******************************************************************
! Sea ice albedo law, set in inifis
!******************************************************************
! Height [kg/m2] used in the calculation of the sea ice albedo vs thickness
! (30 cm rather than 50 cm based on comparisons with Brandt et al. 2005)
! More info on the slab ocean wiki page: https://lmdz-forge.lmd.jussieu.fr/mediawiki/Planets/index.php/Slab_ocean_model
real, save :: h_alb_ice
!$OMP THREADPRIVATE(h_alb_ice)
! Minimum sea ice albedo, reached for a vanishing ice thickness
real, save :: alb_ice_min
!$OMP THREADPRIVATE(alb_ice_min)

END MODULE ocean_slab_h