! ! $Id: ocean_forced_mod.F90 6070 2026-02-06 09:58:17Z fhourdin $ ! MODULE ocean_forced_mod ! ! This module is used for both the sub-surfaces ocean and sea-ice for the case of a ! forced ocean, "ocean=force". ! IMPLICIT NONE CONTAINS ! !**************************************************************************************** ! SUBROUTINE ocean_forced_noice( & itime, dtime, jour, knon, knindex, & p1lay, cdragh, cdragq, cdragm, precip_rain, precip_snow, & temp_air, spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ps, u1, v1, gustiness, tsurf_in, & radsol, snow, agesno, & qsurf, evap, fluxsens, fluxlat, flux_u1, flux_v1, & tsurf_new, dflux_s, dflux_l, sens_prec_liq, rhoa, & dthetadz300, pctsrf, Ampl & #ifdef ISO , xtprecip_rain, xtprecip_snow, xtspechum, Roce, rlat, & xtsnow, xtevap, h1 & #endif ) !$gpum horizontal knon klon !=========================================================================================== ! This subroutine treats the "open ocean" tile in forced mode that is, when LMDZ is not ! coupled to the slab or to NEMO. ! The routine receives data from boundary-file limit.nc and does some calculations at the ! surface. !=========================================================================================== USE dimphy USE calcul_fluxs_mod USE limit_read_mod USE mod_grid_phy_lmdz USE indice_sol_mod USE surface_data, ONLY: iflag_leads USE phys_output_var_mod, ONLY: sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o use config_ocean_skin_m, only: activate_ocean_skin USE calbeta_mod, ONLY: calbeta #ifdef ISO USE infotrac_phy, ONLY: ntiso, niso USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall #ifdef ISOVERIF USE isotopes_mod, ONLY: iso_eau, ridicule !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix USE isotopes_verif_mod #endif #endif USE flux_arp_mod_h USE clesphys_mod_h USE yomcst_mod_h ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime, jour, knon INTEGER, DIMENSION(knon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime REAL, DIMENSION(knon), INTENT(IN) :: p1lay REAL, DIMENSION(knon), INTENT(IN) :: cdragh, cdragq, cdragm REAL, DIMENSION(knon), INTENT(IN) :: precip_rain, precip_snow REAL, DIMENSION(knon), INTENT(IN) :: temp_air, spechum REAL, DIMENSION(knon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ REAL, DIMENSION(knon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV REAL, DIMENSION(knon), INTENT(IN) :: ps REAL, DIMENSION(knon), INTENT(IN) :: u1, v1, gustiness REAL, DIMENSION(knon), INTENT(IN) :: tsurf_in REAL, INTENT(IN) :: rhoa(knon) ! (knon) density of moist air (kg / m3) REAL, DIMENSION(knon), INTENT(IN) :: dthetadz300 REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf #ifdef ISO REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtspechum REAL, DIMENSION(klon), INTENT(IN) :: rlat #endif ! In/Output arguments !**************************************************************************************** REAL, DIMENSION(knon), INTENT(INOUT) :: radsol REAL, DIMENSION(knon), INTENT(INOUT) :: snow REAL, DIMENSION(knon), INTENT(INOUT) :: agesno !? put to 0 in ocean #ifdef ISO REAL, DIMENSION(niso, knon), INTENT(IN) :: xtsnow REAL, DIMENSION(niso, knon), INTENT(INOUT):: Roce #endif ! Output arguments !**************************************************************************************** REAL, DIMENSION(knon), INTENT(OUT) :: qsurf REAL, DIMENSION(knon), INTENT(OUT) :: evap, fluxsens, fluxlat REAL, DIMENSION(knon), INTENT(OUT) :: flux_u1, flux_v1 REAL, DIMENSION(knon), INTENT(OUT) :: tsurf_new REAL, DIMENSION(knon), INTENT(OUT) :: dflux_s, dflux_l REAL, INTENT(out):: sens_prec_liq(:) ! (knon) REAL, DIMENSION(knon), INTENT(OUT) :: Ampl #ifdef ISO REAL, DIMENSION(ntiso, knon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux REAL, DIMENSION(knon), INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation #endif ! Local variables !**************************************************************************************** INTEGER :: i, j REAL, DIMENSION(knon) :: cal, beta, dif_grnd REAL, DIMENSION(knon) :: alb_neig, tsurf_lim, zx_sl REAL, DIMENSION(knon) :: u0, v0 REAL, DIMENSION(knon) :: u1_lay, v1_lay LOGICAL :: check = .FALSE. REAL, DIMENSION(knon) :: sens_prec_sol REAL, DIMENSION(knon) :: lat_prec_liq, lat_prec_sol REAL, DIMENSION(knon) :: l_CBL, sicfra #ifdef ISO REAL, PARAMETER :: t_coup = 273.15 #endif !**************************************************************************************** ! Start calculation !**************************************************************************************** IF (check) WRITE (*, *) ' Entering ocean_forced_noice' #ifdef ISO #ifdef ISOVERIF DO i = 1, knon IF (iso_eau > 0) THEN CALL iso_verif_egalite_choix(xtspechum(iso_eau, i), & & spechum(i), 'ocean_forced_mod 111', & & errmax, errmaxrel) CALL iso_verif_egalite_choix(snow(i), & & xtsnow(iso_eau, i), 'ocean_forced_mod 117', & & errmax, errmaxrel) END IF !IF (iso_eau > 0) THEN END DO !DO i=1,knon #endif #endif !**************************************************************************************** ! 1) ! Read sea-surface temperature from file limit.nc ! !**************************************************************************************** IF (klon_glo .eq. 1) THEN ! single-column model DO i = 1, knon tsurf_lim(i) = tg END DO ELSE ! GCM CALL limit_read_sst(knon, knindex, tsurf_lim & #ifdef ISO , Roce, rlat & #endif ) END IF ! knon !**************************************************************************************** ! 2) ! Flux calculation ! !**************************************************************************************** ! Set some variables for calcul_fluxs ! Use calbeta to calculate beta ! Need to initialize qsurf for calbeta but it is not modified by this routine qsurf(:) = 0. CALL calbeta(dtime, is_oce, knon, snow, qsurf, beta, cal, dif_grnd) alb_neig(:) = 0. agesno(:) = 0. lat_prec_liq = 0.; lat_prec_sol = 0. ! Suppose zero surface speed u0(:) = 0.0 v0(:) = 0.0 u1_lay(:) = u1(:) - u0(:) v1_lay(:) = v1(:) - v0(:) ! Computation of tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf ! EV: does that make sense to have activate_ocean_skin in ocean_forced? ! in anycase; cal=dif=0 for ocean surfaces, so tsurf_new=tsurf_in. ! the two options below should give the same result in principle. IF (activate_ocean_skin == 2) THEN CALL calcul_fluxs(knon, is_oce, dtime, & tsurf_in, p1lay, cal, & beta, cdragh, cdragq, ps, & precip_rain, precip_snow, snow, qsurf, & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) tsurf_new = tsurf_lim ELSE CALL calcul_fluxs(knon, is_oce, dtime, & tsurf_lim, p1lay, cal, & beta, cdragh, cdragq, ps, & precip_rain, precip_snow, snow, qsurf, & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) END IF do j = 1, knon i = knindex(j) sens_prec_liq_o(i, 1) = sens_prec_liq(j) sens_prec_sol_o(i, 1) = sens_prec_sol(j) lat_prec_liq_o(i, 1) = lat_prec_liq(j) lat_prec_sol_o(i, 1) = lat_prec_sol(j) end do IF (iflag_leads .GE. 1) THEN !ym hyper faux, melange de champ compresse et non compresse, je rétabli... !ym l_CBL = -52381.*dthetadz300 + 3008.1 !ym Ampl = 6.012e-08*l_CBL**2 - 4.036e-04*l_CBL + 1.4979 !ym WHERE(Ampl(:)>1.2) Ampl(:)=1.2 !ym sicfra(:)=pctsrf(:,is_sic)/(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter)) !ym WHERE(pctsrf(:,is_sic)+pctsrf(:,is_oce)0.7).and.(sicfra<0.9)) Ampl=((sicfra-0.7)/0.2)*Ampl+((0.9-sicfra)/0.2) !ym fluxsens=Ampl*fluxsens !ym dflux_s=Ampl*dflux_s ! Amplification of surface fluxes over leads regions by G. Gastineau ! Ref: Tian Tian et al. 2025, The Cryosphere ! Isau I. 2007, JGR ! Davy and Gao 2019, zenodo ! Only sensible heat flux is accounted for if iflag_leads=1, evaporative and latent ! heat fluxes are also considered if >1. DO j = 1, knon i = knindex(j) l_CBL(j) = -52381.*dthetadz300(j) + 3008.1 Ampl(j) = 6.012e-08*l_CBL(j)**2 - 4.036e-04*l_CBL(j) + 1.4979 IF (Ampl(j) > 1.2) Ampl(j) = 1.2 sicfra(j) = pctsrf(i, is_sic)/(1.-pctsrf(i, is_lic) - pctsrf(i, is_ter)) IF (pctsrf(i, is_sic) + pctsrf(i, is_oce) < EPSFRA) sicfra(j) = 0. IF (sicfra(j) < 0.7) Ampl(j) = 1. IF (sicfra(j) > 0.7 .and. sicfra(j) < 0.9) Ampl(j) = ((sicfra(j) - 0.7)/0.2)*Ampl(j) + ((0.9 - sicfra(j))/0.2) fluxsens(j) = Ampl(j)*fluxsens(j) dflux_s(j) = Ampl(j)*dflux_s(j) END DO IF (iflag_leads .GE. 2) THEN fluxlat(j) = Ampl(j)*fluxlat(j) dflux_l(j) = Ampl(j)*dflux_l(j) evap(j) = Ampl(j)*evap(j) END IF END IF ! - Flux calculation at first model level for U and V CALL calcul_flux_wind(knon, dtime, & u0, v0, u1, v1, gustiness, cdragm, & AcoefU, AcoefV, BcoefU, BcoefV, & p1lay, temp_air, & flux_u1, flux_v1) #ifdef ISO CALL calcul_iso_surf_oce_vectall(knon, knon, t_coup, & ps, tsurf_new, spechum, u1_lay, v1_lay, xtspechum, & evap, Roce, xtevap, h1 & #ifdef ISOTRAC , knindex & #endif ) #endif #ifdef ISO #ifdef ISOVERIF IF (iso_eau > 0) THEN DO i = 1, knon CALL iso_verif_egalite_choix(snow(i), & xtsnow(iso_eau, i), 'ocean_forced_mod 180', & errmax, errmaxrel) END DO ! DO j=1,knon END IF !IF (iso_eau > 0) THEN #endif #endif END SUBROUTINE ocean_forced_noice ! !*************************************************************************************** ! SUBROUTINE ocean_forced_ice( & itime, dtime, jour, knon, knindex, & tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, & AcoefH, AcoefQ, BcoefH, BcoefQ, & AcoefU, AcoefV, BcoefU, BcoefV, & ps, u1, v1, gustiness, pctsrf, & radsol, snow, qsol, agesno, tsoil, & qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, & icesub, icemelt, & tsurf_new, dflux_s, dflux_l, rhoa, swnet, hice, tice, bilg_cumul, & fcds, fcdi, dh_basal_growth, dh_basal_melt, dh_top_melt, dh_snow2sic, & dtice_melt, dtice_snow2sic & #ifdef ISO , xtprecip_rain, xtprecip_snow, xtspechum, Roce, & xtsnow, xtsol, xtevap, Rland_ice & #endif ) !$gpum horizontal knon klon ! !=========================================================================================== ! This subroutine treats the sea ice tile in when LMDZ is not coupled to the slab or NEMO ! The routine reads data from boundary-condition file limit.nc and does flux calculations ! at the surface. !=========================================================================================== USE dimphy USE geometry_mod, ONLY: longitude, latitude USE calcul_fluxs_mod USE surface_data, ONLY: calice, calsno, iflag_seaice, iflag_seaice_alb, & sice_cond, sisno_cond, sisno_den, sisno_min, sithick_min, sisno_wfact, & amax_s, amax_n, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, & si_pen_frac, si_pen_ext, fseaN, fseaS, iflag_leads, sic_hice_fixed USE geometry_mod, ONLY: longitude, latitude, latitude_deg USE limit_read_mod USE simplehydrol_mod, ONLY: simplehydrol USE indice_sol_mod USE albsno_mod, ONLY: albsno USE soil_mod, ONLY: soil USE calbeta_mod, ONLY: calbeta USE phys_output_var_mod, ONLY: sens_prec_liq_o, sens_prec_sol_o, lat_prec_liq_o, lat_prec_sol_o #ifdef ISO USE infotrac_phy, ONLY: niso, ntiso USE isotopes_routines_mod, ONLY: calcul_iso_surf_oce_vectall, calcul_iso_surf_sic_vectall #ifdef ISOVERIF USE isotopes_mod, ONLY: iso_eau, ridicule !USE isotopes_verif_mod, ONLY: errmax,errmaxrel,iso_verif_egalite_choix USE isotopes_verif_mod #endif #endif USE flux_arp_mod_h USE clesphys_mod_h USE yomcst_mod_h USE dimsoil_mod_h, ONLY: nsoilmx ! Input arguments !**************************************************************************************** INTEGER, INTENT(IN) :: itime, jour, knon INTEGER, DIMENSION(knon), INTENT(IN) :: knindex REAL, INTENT(IN) :: dtime REAL, DIMENSION(knon), INTENT(IN) :: tsurf_in REAL, DIMENSION(knon), INTENT(IN) :: p1lay REAL, DIMENSION(knon), INTENT(IN) :: cdragh, cdragm REAL, DIMENSION(knon), INTENT(IN) :: precip_rain, precip_snow REAL, DIMENSION(knon), INTENT(IN) :: temp_air, spechum REAL, DIMENSION(knon), INTENT(IN) :: AcoefH, AcoefQ, BcoefH, BcoefQ REAL, DIMENSION(knon), INTENT(IN) :: AcoefU, AcoefV, BcoefU, BcoefV REAL, DIMENSION(knon), INTENT(IN) :: ps REAL, DIMENSION(knon), INTENT(IN) :: u1, v1, gustiness real, intent(in):: rhoa(knon) ! (knon) density of moist air (kg / m3) REAL, DIMENSION(knon), INTENT(IN) :: swnet REAL, DIMENSION(klon, nbsrf), INTENT(IN) :: pctsrf #ifdef ISO REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtprecip_rain, xtprecip_snow REAL, DIMENSION(ntiso, knon), INTENT(IN) :: xtspechum REAL, DIMENSION(niso, knon), INTENT(IN) :: Roce REAL, DIMENSION(niso, knon), INTENT(IN) :: Rland_ice #endif ! In/Output arguments !**************************************************************************************** REAL, DIMENSION(knon), INTENT(INOUT) :: radsol REAL, DIMENSION(knon), INTENT(INOUT) :: snow, qsol REAL, DIMENSION(knon), INTENT(INOUT) :: agesno REAL, DIMENSION(knon, nsoilmx), INTENT(INOUT) :: tsoil REAL, DIMENSION(klon), INTENT(INOUT) :: hice !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: tice !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: fcds !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: fcdi !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_growth !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: dh_basal_melt !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: dh_top_melt !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: dh_snow2sic !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_melt !ym WARNING uncompressed REAL, DIMENSION(klon), INTENT(INOUT) :: dtice_snow2sic !ym WARNING uncompressed #ifdef ISO REAL, DIMENSION(niso, knon), INTENT(INOUT) :: xtsnow REAL, DIMENSION(niso, knon), INTENT(IN) :: xtsol #endif ! Output arguments !**************************************************************************************** REAL, DIMENSION(knon), INTENT(OUT) :: qsurf REAL, DIMENSION(knon), INTENT(OUT) :: alb1_new ! new albedo in visible SW interval REAL, DIMENSION(knon), INTENT(OUT) :: alb2_new ! new albedo in near IR interval REAL, DIMENSION(knon), INTENT(OUT) :: evap, fluxsens, fluxlat REAL, DIMENSION(knon), INTENT(OUT) :: flux_u1, flux_v1 REAL, DIMENSION(knon), INTENT(OUT) :: tsurf_new REAL, DIMENSION(knon), INTENT(OUT) :: dflux_s, dflux_l REAL, DIMENSION(knon), INTENT(OUT) :: icesub, icemelt #ifdef ISO REAL, DIMENSION(ntiso, knon), INTENT(OUT) :: xtevap #endif ! Local variables !**************************************************************************************** LOGICAL, PARAMETER :: check = .FALSE. INTEGER :: i, j REAL :: zfra REAL, PARAMETER :: t_grnd = 271.35 REAL, DIMENSION(knon) :: cal, beta, dif_grnd, capsol REAL, DIMENSION(knon) :: alb_neig, tsurf_tmp REAL, DIMENSION(knon) :: soilcap, soilflux REAL, DIMENSION(knon) :: u0, v0 REAL, DIMENSION(knon) :: u1_lay, v1_lay REAL, DIMENSION(knon) :: sens_prec_liq, sens_prec_sol REAL, DIMENSION(knon) :: lat_prec_liq, lat_prec_sol REAL, DIMENSION(knon) :: yhice ! compressed version of hice INTEGER :: ki INTEGER :: cpl_pas REAL, DIMENSION(knon) :: bilg REAL, DIMENSION(knon) :: fsic REAL, DIMENSION(knon) :: f_bot REAL, PARAMETER :: latent_ice = 334.0e3 REAL, PARAMETER :: rau_ice = 917.0 REAL, PARAMETER :: kice = 2.2 REAL :: f_cond, f_swpen, f_cond_s, f_cond_i REAL :: ustar, uscap, ustau ! for snow/ice albedo: REAL :: alb_snow, alb_ice, alb_pond REAL :: frac_snow, frac_ice, frac_pond REAL :: z1_i, z2_i, z1_s, zlog ! height parameters ! for ice melt / freeze REAL :: e_melt, snow_evap, h_test ! dhsic, dfsic change in ice mass, fraction. REAL :: dhsic, dfsic, frac_mf REAL :: fsea, amax REAL :: hice_i, tice_i, fsic_new ! snow and ice physical characteristics: REAL, PARAMETER :: t_freeze = 271.35 ! freezing sea water temp REAL, PARAMETER :: t_melt = 273.15 ! melting ice temp REAL :: sno_den!=sisno_den !mean snow density, kg/m3 REAL, PARAMETER :: ice_den = 917. ! ice density REAL, PARAMETER :: sea_den = 1025. ! sea water density REAL :: ice_cond!=sice_cond*ice_den !conductivity of ice REAL :: sno_cond!=sisno_cond*sno_den ! conductivity of snow REAL, PARAMETER :: ice_cap = 2067. ! specific heat capacity, snow and ice REAL, PARAMETER :: sea_cap = 3995. ! specific heat capacity, water REAL, PARAMETER :: ice_lat = 334000. ! freeze /melt latent heat snow and ice ! control of snow and ice cover & freeze / melt (heights converted to kg/m2) REAL :: snow_min!=sisno_min*sno_den !critical snow height 5 cm REAL :: snow_wfact!=sisno_wfact ! max fraction of falling snow blown into ocean REAL, PARAMETER :: ice_frac_min = 0.005 REAL :: h_ice_min!=sithick_min ! min ice thickness ! below ice_thin, priority is melt lateral / grow height ! ice_thin is also height of new ice REAL, PARAMETER :: h_ice_max = 7 ! max ice height ! Ice thickness parameter for lateral growth REAL, PARAMETER :: h_ice_thick = 1.5 REAL, PARAMETER :: h_ice_thin = 0.15 ! albedo and radiation parameters INTEGER :: iflag_sic_albedo ! albedo old or NEMO REAL :: alb_sno_dry!=rn_alb_sdry !dry snow albedo REAL :: alb_sno_wet!=rn_alb_smlt !wet snow albedo REAL :: alb_ice_dry!=rn_alb_idry !dry thick ice REAL :: alb_ice_wet!=rn_alb_imlt !melting thick ice ! Toyoda 2020' albedo ! Values for snow / ice, dry / melting, visible / near IR REAL, PARAMETER :: alb_sdry_vis = 0.98 REAL, PARAMETER :: alb_smlt_vis = 0.88 REAL, PARAMETER :: alb_sdry_nir = 0.7 REAL, PARAMETER :: alb_smlt_nir = 0.55 REAL, PARAMETER :: alb_idry_vis = 0.78 REAL, PARAMETER :: alb_imlt_vis = 0.705 REAL, PARAMETER :: alb_idry_nir = 0.36 REAL, PARAMETER :: alb_imlt_nir = 0.285 REAL, PARAMETER :: h_ice_alb = 0.5*ice_den ! height for full ice albedo REAL, PARAMETER :: h_sno_alb = 0.02*300 ! height for control of snow fraction REAL :: pen_frac !=si_pen_frac !fraction of shortwave penetrating into the ! ice (not snow). Should be visible only, not NIR REAL :: pen_ext !=si_pen_ext !extinction length of penetrating shortwave (m-1) REAl :: lon(knon), lat(knon) ! for indexation #ifdef ISO REAL, PARAMETER :: t_coup = 273.15 REAL, DIMENSION(knon) :: fq_fonte_diag REAL, DIMENSION(knon) :: fqfonte_diag REAL, DIMENSION(knon) :: snow_evap_diag REAL, DIMENSION(knon) :: fqcalving_diag REAL, DIMENSION(knon) :: run_off_lic_diag REAL :: coeff_rel_diag REAL :: max_eau_sol_diag REAL, DIMENSION(knon) :: runoff_diag INTEGER IXT REAL, DIMENSION(niso, knon) :: xtsnow_prec, xtsol_prec REAL, DIMENSION(knon) :: snow_prec, qsol_prec #endif !**************************************************************************************** ! Initialisation !**************************************************************************************** IF (check) WRITE (*, *) 'Entering surface_seaice, knon=', knon iflag_sic_albedo = iflag_seaice_alb IF (iflag_seaice .GE. 1) THEN sno_den = sisno_den !mean snow density, kg/m3 ice_cond = sice_cond*ice_den !conductivity of ice sno_cond = sisno_cond*sno_den ! conductivity of snow snow_min = sisno_min*sno_den !critical snow height 5 cm snow_wfact = sisno_wfact ! max fraction of falling snow blown into ocean h_ice_min = sithick_min ! min ice thickness alb_sno_dry = rn_alb_sdry !dry snow albedo alb_sno_wet = rn_alb_smlt !wet snow albedo alb_ice_dry = rn_alb_idry !dry thick ice alb_ice_wet = rn_alb_imlt !melting thick ice pen_frac = si_pen_frac !fraction of shortwave penetrating into the pen_ext = si_pen_ext !extinction length of penetrating shortwave (m-1) cpl_pas = NINT(86400./dtime*1.0) ! once a day END IF !*************************************************************************************************** ! 1) Heat diffusion in the snow/sea-ice, surface turbulent flux and surface temperature calculation ! if iflag_seaice=0, we use the old-style method considering a constant inertia for ! sea ice (or uperlying snow) and eventually using the soil routine to compute heat diffusion ! in sea ice ! if iflag_seaice>1, we compute the sea-ice thermal properties using the scheme of the slab ! from Liu, Risi, Codron et al. 2021, nature comm. ! if iflag_seaice=1, we read the sea-ice thickness from limit.nc ! if iflag_seaice=2, sea-ice thickness is computed using a prognostic evolution model. ! see also PhD thesis by N. Michalezyk (LOCEAN) !*************************************************************************************************** tsurf_tmp(:) = tsurf_in(:) IF (iflag_seaice == 0) THEN ! Old LMDZ sea ice surface ! calculate the parameters cal, beta, capsol and dif_grnd and then recalculate cal CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd) IF (soil_model) THEN ! update tsoil and calculate soilcap and soilflux lon(1:knon) = longitude(knindex(1:knon)) lat(1:knon) = latitude(knindex(1:knon)) CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, & & lon, lat, tsoil, soilcap, soilflux) cal(1:knon) = RCPD/soilcap(1:knon) radsol(1:knon) = radsol(1:knon) + soilflux(1:knon) dif_grnd = 1.0/tau_gl ELSE dif_grnd = 1.0/tau_gl cal = RCPD*calice WHERE (snow > 0.0) cal = RCPD*calsno END IF ELSEIF (iflag_seaice .GE. 1) THEN ! if iflag_seaice>=1, advanced computation of snow and sea-ice thermal properties. ! They become dependent on sea-ice thickness hice. bilg(:) = 0. dif_grnd(:) = 0. beta(:) = 1. cpl_pas = NINT(86400./dtime*1.0) ! once a day DO i = 1, knon ki = knindex(i) fsic(i) = pctsrf(ki, is_sic) IF (snow(i) .GT. snow_min) THEN ! 1 / snow-layer heat capacity cal(i) = 2.*RCPD/(snow(i)*ice_cap) ! adjustment time-scale of conductive flux dif_grnd(i) = cal(i)*sno_cond/snow(i)/RCPD ! for conductive flux f_cond_s = sno_cond*(tice(ki) - t_freeze)/snow(i) radsol(i) = radsol(i) + f_cond_s ! all shortwave flux absorbed f_swpen = 0. ELSE ! bare ice. f_cond_s = 0. ! 1 / ice-layer heat capacity cal(i) = 2.*RCPD/(hice(ki)*ice_den*ice_cap) ! adjustment time-scale of conductive flux dif_grnd(i) = cal(i)*ice_cond/(hice(ki)*ice_den)/RCPD ! penetrative shortwave flux... f_swpen = swnet(i)*pen_frac*exp(-pen_ext*hice(ki)) radsol(i) = radsol(i) - f_swpen ! GG no conductive flux in this case? END IF bilg(i) = f_swpen - f_cond_s END DO END IF lat_prec_liq(:) = 0.0 lat_prec_sol(:) = 0.0 ! Suppose zero surface speed u0(:) = 0.0 v0(:) = 0.0 u1_lay(:) = u1(:) - u0(:) v1_lay(:) = v1(:) - v0(:) CALL calcul_fluxs(knon, is_sic, dtime, & tsurf_tmp, p1lay, cal, beta, cdragh, cdragh, ps, & precip_rain, precip_snow, snow, qsurf, & radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, & f_qsat_oce, AcoefH, AcoefQ, BcoefH, BcoefQ, & tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l, & sens_prec_liq, sens_prec_sol, lat_prec_liq, lat_prec_sol, rhoa) DO j = 1, knon i = knindex(j) sens_prec_liq_o(i, 2) = sens_prec_liq(j) sens_prec_sol_o(i, 2) = sens_prec_sol(j) lat_prec_liq_o(i, 2) = lat_prec_liq(j) lat_prec_sol_o(i, 2) = lat_prec_sol(j) END DO ! - Flux calculation at first modele level for U and V CALL calcul_flux_wind(knon, dtime, & u0, v0, u1, v1, gustiness, cdragm, & AcoefU, AcoefV, BcoefU, BcoefV, & p1lay, temp_air, & flux_u1, flux_v1) ! - If iflag_seaice >=1, we also compute ice temperature and heat flux in the ice IF (iflag_seaice .GE. 1) THEN DO i = 1, knon ki = knindex(i) ! ocean-ice heat flux fsea = fseaS amax = amax_s if (latitude(ki) > 0) THEN fsea = fseaN amax = amax_n END IF IF (snow(i) .GT. snow_min) THEN ! snow conductive flux after calcul_fluxs (pos up) f_cond_s = sno_cond*(tice(ki) - tsurf_new(i))/snow(i) ! 1 / heat capacity and conductive timescale uscap = 2./ice_cap/(snow(i) + hice(ki)*ice_den) ustau = uscap*ice_cond/(hice(ki)*ice_den) ! update ice temp tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) & /(1 + dtime*ustau) ELSE ! bare ice tice(ki) = tsurf_new(i) END IF ! ice conductive flux (pos up) f_cond_i = ice_cond*(t_freeze - tice(ki))/(hice(ki)*ice_den) f_bot(i) = fsea - f_cond_i fcdi(ki) = f_cond_i - fsea fcds(i) = f_cond_s !bilg(i) = bilg(i)+f_cond_i END DO END IF !**************************************************************************************** ! 2) Treatment of snow and ice hydrology and computation of sea ice thickness ! two different treatments if iflag_seaice = 0 (old-style) or >0 (slab style) !**************************************************************************************** IF (iflag_seaice .EQ. 0) THEN ! if iflag_seaice = 0, we fix the sea-ice thickness to a constant value ! this only serves for albedo calculation if iflag_sic_albedo >=1 DO i = 1, knon ki = knindex(i) hice(ki) = sic_hice_fixed END DO #ifdef ISO #ifdef ISOVERIF DO i = 1, knon IF (iso_eau > 0) THEN IF (snow(i) > ridicule) THEN CALL iso_verif_egalite_choix(xtsnow(iso_eau, i), snow(i), & & 'interfsurf 964', errmax, errmaxrel) END IF !IF ((snow(i) > ridicule)) THEN END IF !IF (iso_eau > 0) THEN END DO !DO i=1,knon #endif DO i = 1, knon snow_prec(i) = snow(i) DO ixt = 1, niso xtsnow_prec(ixt, i) = xtsnow(ixt, i) END DO !DO ixt=1,niso ! initialisation: fq_fonte_diag(i) = 0.0 fqfonte_diag(i) = 0.0 snow_evap_diag(i) = 0.0 END DO !DO i=1,knon #endif CALL simplehydrol(knon, is_sic, knindex, dtime, & tsurf_tmp, precip_rain, precip_snow, & snow, qsol, tsurf_new, evap, icesub, icemelt & #ifdef ISO , fq_fonte_diag, fqfonte_diag, snow_evap_diag, fqcalving_diag & , max_eau_sol_diag, runoff_diag, run_off_lic_diag, coeff_rel_diag & #endif ) #ifdef ISO CALL calcul_iso_surf_sic_vectall(knon, knon, & evap, snow_evap_diag, Tsurf_new, Roce, snow, & fq_fonte_diag, fqfonte_diag, dtime, t_coup, & precip_snow, xtprecip_snow, xtprecip_rain, snow_prec, xtsnow_prec, & xtspechum, spechum, ps, & xtevap, xtsnow, fqcalving_diag, & knindex, is_sic, run_off_lic_diag, coeff_rel_diag, Rland_ice) #ifdef ISOVERIF IF (iso_eau > 0) THEN DO i = 1, knon CALL iso_verif_egalite_choix(snow(i), & xtsnow(iso_eau, i), 'ocean_forced_mod 396', & errmax, errmaxrel) END DO ! DO j=1,knon END IF !IF (iso_eau > 0) then #endif !#ifdef ISOVERIF #endif !#ifdef ISO ELSE IF (iflag_seaice .EQ. 1) THEN ! if iflag_seaice=1 we read the thickness of the sea ice in the limit.nc file !ym totally wrong since "hice" is an uncompressed field (klon) and limit_read_hice return a compressed field (knon) !ym CALL limit_read_hice(knon,knindex,hice) CALL limit_read_hice(knon, knindex, yhice) DO i = 1, knon ki = knindex(i) hice(ki) = yhice(i) END DO END IF !NB another idea could be testing the formula by Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min)) DO i = 1, knon ki = knindex(i) !ym WARNING : fsic uninitilazed, take initialisation similar to iflag_seaice==2 fsic(i) = pctsrf(ki, is_sic) IF (precip_snow(i) > 0.) THEN snow(i) = snow(i) + precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(i))) END IF ! snow and ice sublimation IF (evap(i) > 0.) THEN snow_evap = MIN(snow(i)/dtime, evap(i)) snow(i) = snow(i) - snow_evap*dtime snow(i) = MAX(0.0, snow(i)) IF (iflag_seaice == 2) THEN hice(ki) = MAX(0.0, hice(ki) - (evap(i) - snow_evap)*dtime/ice_den) END IF END IF ! Melt / Freeze snow from above if Tsurf>0 IF (tsurf_new(i) .GT. t_melt) THEN ! energy available for melting snow (in kg of melted snow /m2) e_melt = MIN(MAX(snow(i)*(tsurf_new(i) - t_melt)*ice_cap/2. & /(ice_lat + ice_cap/2.*(t_melt - tice(ki))), 0.0), snow(i)) ! remove snow tice_i = tice(ki) IF (snow(i) .GT. e_melt) THEN snow(i) = snow(i) - e_melt tsurf_new(i) = t_melt ELSE ! all snow is melted ! add remaining heat flux to ice e_melt = e_melt - snow(i) tice(ki) = tice(ki) + e_melt*ice_lat*2./(ice_cap*hice(ki)*ice_den) tsurf_new(i) = tice(ki) END IF dtice_melt(ki) = (tice(ki) - tice_i)/dtime END IF IF (iflag_seaice .EQ. 2) THEN ! Interactive calculation of hice if iflag_seaice ==2 ! Bottom melt / grow ! bottom freeze if bottom flux (cond + oce-ice) <0 IF (f_bot(i) .LT. 0) THEN ! larger fraction of bottom growth frac_mf = MIN(1., MAX(0., (hice(ki) - h_ice_thick) & /(h_ice_max - h_ice_thick))) ! quantity of new ice (formed at mean ice temp) e_melt = -f_bot(i)*dtime*fsic(i) & /(ice_lat + ice_cap/2.*(t_freeze - tice(ki))) ! first increase height to h_thick dhsic = MAX(0., MIN(h_ice_thick - hice(ki), e_melt/(fsic(i)*ice_den))) hice_i = hice(ki) hice(ki) = dhsic + hice(ki) e_melt = e_melt - fsic(i)*dhsic IF (e_melt .GT. 0.) THEN ! frac_mf fraction used for lateral increase dfsic = MIN(amax - fsic(i), e_melt*frac_mf/(hice(ki)*ice_den)) ! No lateral growth -> forced ocean !fsic(ki)=fsic(ki)+dfsic e_melt = e_melt - dfsic*hice(ki)*ice_den ! rest used to increase height hice(ki) = MIN(h_ice_max, hice(ki) + e_melt/(fsic(i)*ice_den)) END IF dh_basal_growth(ki) = (hice(ki) - hice_i)/dtime ! melt from below if bottom flux >0 ELSE ! larger fraction of lateral melt from warm ocean frac_mf = MIN(1., MAX(0., (hice(ki) - h_ice_thin) & /(h_ice_thick - h_ice_thin))) ! bring ice to freezing and melt from below ! quantity of melted ice e_melt = f_bot(i)*dtime*fsic(i) & /(ice_lat + ice_cap/2.*(tice(ki) - t_freeze)) ! first decrease height to h_thick hice_i = hice(ki) dhsic = MAX(0., MIN(hice(ki) - h_ice_thick, e_melt/(fsic(i)*ice_den))) hice(ki) = hice(ki) - dhsic e_melt = e_melt - fsic(i)*dhsic*ice_den IF (e_melt .GT. 0) THEN ! frac_mf fraction used for height decrease dhsic = MAX(0., MIN(hice(ki) - h_ice_min, e_melt/ice_den*frac_mf/fsic(i))) hice(ki) = hice(ki) - dhsic e_melt = e_melt - fsic(i)*dhsic*ice_den ! rest used to decrease fraction (up to 0!) dfsic = MIN(fsic(i), e_melt/(hice(ki)*ice_den)) ! Remaining heat not used if everything melted e_melt = e_melt - dfsic*hice(ki)*ice_den ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime END IF dh_basal_melt(ki) = (hice(ki) - hice_i)/dtime END IF END IF ! melt ice from above if Tice>0 tice_i = tice(ki) IF (tice(ki) .GT. t_melt) THEN IF (iflag_seaice .EQ. 2) THEN ! quantity of ice melted (kg/m2) e_melt = MAX(hice(ki)*ice_den*(tice(ki) - t_melt)*ice_cap/2. & /(ice_lat + ice_cap/2.*(t_melt - t_freeze)), 0.0) ! melt from above, height only hice_i = hice(ki) dhsic = MIN(hice(ki) - h_ice_min, e_melt/ice_den) dh_top_melt(i) = dhsic e_melt = e_melt - dhsic IF (e_melt .GT. 0) THEN ! lateral melt if ice too thin dfsic = MAX(fsic(i) - ice_frac_min, e_melt/(h_ice_min*ice_den)*fsic(i)) ! if all melted do nothing with remaining heat e_melt = MAX(0., e_melt*fsic(i) - dfsic*h_ice_min*ice_den) ! slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime END IF hice(ki) = hice(ki) - dhsic dh_top_melt(ki) = (hice(ki) - hice_i)/dtime ! surface temperature at melting point END IF tice(ki) = t_melt tsurf_new(i) = t_melt END IF dtice_melt(ki) = dtice_melt(ki) + tice(ki) - tice_i ! convert snow to ice if below floating line h_test = (hice(ki)*ice_den + snow(i)) - hice(ki)*sea_den IF ((h_test .GT. 0.) .AND. (hice(ki) .GT. h_ice_min)) THEN !snow under water ! extra snow converted to ice (with added frozen sea water) IF (iflag_seaice .EQ. 2) THEN dhsic = h_test/(sea_den - ice_den + sno_den) hice(ki) = hice(ki) + dhsic END IF snow(i) = snow(i) - dhsic*sno_den ! available energy (freeze sea water + bring to tice) e_melt = dhsic*ice_den*(1.-sno_den/ice_den)*(ice_lat + & ice_cap/2.*(t_freeze - tice(ki))) ! update ice temperature tice_i = tice(ki) tice(ki) = tice(ki) + 2.*e_melt/ice_cap/(snow(i) + hice(ki)*ice_den) IF (iflag_seaice .EQ. 2) THEN dh_snow2sic(ki) = dhsic/dtime END IF dtice_snow2sic(ki) = (tice(ki) - tice_i)/dtime END IF END DO ! cumulated ice-ocean fluxes, update tslab, lateral grow DO j = i, knon ki = knindex(i) bilg_cumul(ki) = bilg_cumul(ki) + bilg(j)/float(cpl_pas) END DO !YM Pb for GPU port, done on an compressed field inside a compressed Kernel !YM ==> move outside of subroutine !YM CALL ocean_forced_ice_reset_bilg_cumul() !YM IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab !YM bilg_cumul(1:klon)=0. !YM END IF ! coupling time !tests ice fraction !YM Pb for GPU port : done on uncompressed field inside a compressed kernel !YM update fsic only on compressed index !YM WHERE (fsic.LT.ice_frac_min) !YM tice=t_melt !YM hice=h_ice_min !YM END WHERE DO j = i, knon ki = knindex(i) IF (fsic(i) .LT. ice_frac_min) THEN tice(ki) = t_melt hice(ki) = h_ice_min END IF END DO END IF ! if iflag_seaice .eq. 0 !**************************************************************************************** ! 3) Compute sea-ice / snow albedo ! there are currently 4 options !**************************************************************************************** SELECT CASE (iflag_sic_albedo) CASE (0) ! default option. Calculation of albedo at snow (alb_neig) and update the age of snow (agesno) ! Note that where the is no snow, the albedo of bare sea ice is taken as the ! value of aged snow (computation with agesno=0). ! when looking at albsno, it is in fact a typical value of a desert (0.55)! CALL albsno(knon, knon, dtime, agesno(:), alb_neig(:), precip_snow(:)) WHERE (snow(1:knon) .LT. 0.0001) agesno(1:knon) = 0. alb1_new(:) = 0.0 DO i = 1, knon zfra = MAX(0.0, MIN(1.0, snow(i)/(snow(i) + 10.0))) alb1_new(i) = alb_neig(i)*zfra + 0.6*(1.0 - zfra) END DO alb2_new(:) = alb1_new(:) CASE (1) ! New visible and IR albedos, dry / melting snow ! based on Toyoda et al, 2020, from the slab model DO i = 1, knon ki = knindex(i) ! snow fraction frac_snow = snow(i)/(snow(i) + h_sno_alb) ! dependence of ice albedo with ice thickness frac_ice = MIN(1., ATAN(4.*hice(ki)*ice_den)/ATAN(4.*h_ice_alb)) ! Total (for ice, min = 0.066 = alb_ocean) IF (tice(ki) .GT. t_melt) THEN alb_ice = 0.066 + (alb_imlt_vis - 0.066)*frac_ice alb1_new(i) = alb_smlt_vis*frac_snow + alb_ice*(1.-frac_snow) alb_ice = 0.066 + (alb_imlt_nir - 0.066)*frac_ice alb2_new(i) = alb_smlt_nir*frac_snow + alb_ice*(1.-frac_snow) ELSEIF (tice(ki) .GT. t_melt - 1.) THEN frac_pond = tice(ki) - t_freeze alb_snow = alb_smlt_vis*frac_pond + alb_sdry_vis*(1.-frac_pond) alb_ice = alb_imlt_vis*frac_pond + alb_idry_vis*(1.-frac_pond) alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice alb1_new(i) = alb_snow*frac_snow + alb_ice*(1.-frac_snow) alb_snow = alb_smlt_nir*frac_pond + alb_sdry_nir*(1.-frac_pond) alb_ice = alb_imlt_nir*frac_pond + alb_idry_nir*(1.-frac_pond) alb_ice = 0.066 + (alb_ice - 0.66)*frac_ice alb2_new(i) = alb_snow*frac_snow + alb_ice*(1.-frac_snow) ELSE alb_ice = 0.066 + (alb_idry_vis - 0.066)*frac_ice alb1_new(i) = alb_sdry_vis*frac_snow + alb_ice*(1.-frac_snow) alb_ice = 0.066 + (alb_idry_nir - 0.066)*frac_ice alb2_new(i) = alb_sdry_nir*frac_snow + alb_ice*(1.-frac_snow) END IF END DO CASE (2) ! LIM3 scheme. Uses clear sky / overcast value, assuming 50% clear sky ! note that this is what is implemented in the coupled model (CMIP5, 6, 7) z1_i = 1.5*ice_den z2_i = 0.05*ice_den zlog = 1./(LOG(1.5) - LOG(0.05)) z1_s = 1./(0.025*sno_den) DO i = 1, knon ki = knindex(i) ! temperature above / below 0 IF (tice(ki) .GE. t_melt) THEN alb_ice = alb_ice_wet alb_snow = alb_sno_wet ELSE alb_ice = alb_ice_dry alb_snow = alb_sno_dry END IF ! ice thickness IF (hice(ki)*ice_den .LT. z2_i) THEN alb_ice = 0.066 + 0.114*hice(ki)*ice_den/z2_i ELSEIF (hice(ki)*ice_den .LT. z1_i) THEN alb_ice = alb_ice + (0.18 - alb_ice) & *(LOG(z1_i) - LOG((hice(ki)+0.001)*ice_den))*zlog END IF ! ice or snow depending on snow thickness alb_snow = alb_snow - (alb_snow - alb_ice)*EXP(-snow(i)*z1_s) ! Effect of clouds (polynomial fit with 50% clouds) alb1_new(i) = alb_snow - 0.5*(-0.1010*alb_snow*alb_snow & + 0.1933*alb_snow - 0.0148) alb2_new(i) = alb1_new(i) END DO CASE (3) ! old slab albedo : single value. age of snow, melt ponds. DO i = 1, knon ki = knindex(i) ! snow albedo: update snow age IF (snow(i) .GT. 0.0001) THEN agesno(i) = (agesno(i) + (1.-agesno(i)/50.)*dtime/86400.) & *EXP(-1.*MAX(0.0, precip_snow(i))*dtime/5.) ELSE agesno(i) = 0.0 END IF ! snow albedo alb_snow = alb_sno_wet + (alb_sno_dry - alb_sno_wet)*EXP(-agesno(i)/50.) ! ice albedo (varies with ice thickness and temp) alb_ice = MAX(0.0, 0.13*LOG(100.*(hice(ki)+0.001)) + 0.1) !alb_ice=MAX(0.0,0.13*LOG(100.*seaice(ki)/ice_den)+0.1) IF (tice(ki) .GT. t_freeze - 0.01) THEN alb_ice = MIN(alb_ice, alb_ice_wet) ELSE alb_ice = MIN(alb_ice, alb_ice_dry) END IF ! pond albedo alb_pond = 0.36 - 0.1*(2.0 + MIN(0.0, MAX(tice(ki) - t_melt, -2.0))) ! pond fraction frac_pond = 0.2*(2.0 + MIN(0.0, MAX(tice(ki) - t_melt, -2.0))) ! snow fraction frac_snow = MAX(0.0, MIN(1.0 - frac_pond, snow(i)/snow_min)) ! ice fraction frac_ice = MAX(0.0, 1.-frac_pond - frac_snow) ! total albedo alb1_new(i) = alb_snow*frac_snow + alb_ice*frac_ice + alb_pond*frac_pond END DO alb2_new(:) = alb1_new(:) END SELECT ! ------ End Albedo ---------- END SUBROUTINE ocean_forced_ice !************************************************************************************* SUBROUTINE ocean_forced_ice_reset_bilg_cumul(itime, dtime, bilg_cumul) USE dimphy, ONLY: klon USE surface_data, ONLY: iflag_seaice, type_ocean, version_ocean IMPLICIT NONE INTEGER, INTENT(IN) :: itime REAL, INTENT(IN) :: dtime REAL, DIMENSION(klon), INTENT(INOUT) :: bilg_cumul INTEGER :: cpl_pas IF (iflag_seaice /= 0) THEN cpl_pas = NINT(86400./dtime*1.0) ! une fois par jour IF (MOD(itime, cpl_pas) .EQ. 0) THEN ! time to update tslab IF (.NOT. (type_ocean == 'couple' .OR. (type_ocean == 'slab' .AND. version_ocean == 'sicINT'))) THEN bilg_cumul(1:klon) = 0. END IF END IF ! coupling time END IF END SUBROUTINE ocean_forced_ice_reset_bilg_cumul END MODULE ocean_forced_mod