!
! $Id: ocean_forced_mod.F90 5662 2025-05-20 14:24:41Z fairhead $
!
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            
       )
!
! This subroutine treats the "open ocean", all grid points that are not entierly covered
! by ice.
! The routine receives data from climatologie 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
#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(klon), INTENT(IN)     :: knindex
    REAL, INTENT(IN)                         :: dtime
    REAL, DIMENSION(klon), INTENT(IN)        :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)        :: cdragh, cdragq, cdragm
    REAL, DIMENSION(klon), INTENT(IN)        :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)        :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefH, AcoefQ, BcoefH, BcoefQ
    REAL, DIMENSION(klon), INTENT(IN)        :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)        :: ps
    REAL, DIMENSION(klon), INTENT(IN)        :: u1, v1, gustiness
    REAL, DIMENSION(klon), INTENT(IN)        :: tsurf_in
    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
!GG
     REAL, DIMENSION(klon), INTENT(IN)        :: dthetadz300
     REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
!

#ifdef ISO
    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtprecip_rain, xtprecip_snow
    REAL, DIMENSION(ntiso,klon), INTENT(IN)  :: xtspechum
    REAL, DIMENSION(klon),       INTENT(IN)  :: rlat
#endif

! In/Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)     :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)     :: snow
    REAL, DIMENSION(klon), INTENT(INOUT)     :: agesno !? put to 0 in ocean
#ifdef ISO     
    REAL, DIMENSION(niso,klon), INTENT(IN)   :: xtsnow
    REAL, DIMENSION(niso,klon), INTENT(INOUT):: Roce
#endif

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)       :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)       :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)       :: flux_u1, flux_v1
    REAL, DIMENSION(klon), INTENT(OUT)       :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)       :: dflux_s, dflux_l
    REAL, INTENT(out):: sens_prec_liq(:) ! (knon)
!GG
     REAL, DIMENSION(klon), INTENT(OUT)       :: Ampl
!

#ifdef ISO     
    REAL, DIMENSION(ntiso,klon), INTENT(OUT) :: xtevap ! isotopes in evaporation flux
    REAL, DIMENSION(klon),       INTENT(OUT) :: h1 ! just a diagnostic, not useful for the simulation
#endif

! Local variables
!****************************************************************************************
    INTEGER                     :: i, j
    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd
    REAL, DIMENSION(klon)       :: alb_neig, tsurf_lim, zx_sl
    REAL, DIMENSION(klon)       :: u0, v0
    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
    LOGICAL                     :: check=.FALSE.
    REAL, DIMENSION(knon)       :: sens_prec_sol
    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
! GG
    REAL, DIMENSION(klon)       :: 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)
      ENDIF !IF (iso_eau > 0) THEN
    ENDDO !DO i=1,knon
#endif      
#endif 

!****************************************************************************************
! 1)    
! Read sea-surface temperature from file limit.nc
!
!****************************************************************************************
!--sb:
!!jyg    if (knon.eq.1) then ! single-column model
    if (klon_glo.eq.1) then ! single-column model
      ! EV: now surface Tin flux_arp.h
      !CALL read_tsurf1d(knon,tsurf_lim) ! new
       DO i = 1, knon
        tsurf_lim(i) = tg
       ENDDO

    else ! GCM
      CALL limit_read_sst(knon,knindex,tsurf_lim &
#ifdef ISO
     &     ,Roce,rlat &
#endif     
     &     )
    endif ! knon
!sb--

!****************************************************************************************
! 2)
! Flux calculation
!
!****************************************************************************************
! Set some variables for calcul_fluxs
    !cal = 0.
    !beta = 1.
    !dif_grnd = 0.
    
    
    ! EV: 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(:)

! Calcul de tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l and qsurf
    CALL calcul_fluxs(knon, is_oce, dtime, &
         merge(tsurf_in, tsurf_lim, activate_ocean_skin == 2), 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)
    if (activate_ocean_skin == 2) tsurf_new = tsurf_lim

    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)
    enddo

!GG
    if (iflag_leads == 1) then
      l_CBL = -52381.*dthetadz300 + 3008.1
      Ampl = 6.012e-08*l_CBL**2 - 4.036e-04*l_CBL + 1.4979
      WHERE(Ampl(:)>1.2) Ampl(:)=1.2
      sicfra(:)=pctsrf(:,is_sic)/(1.-pctsrf(:,is_lic)-pctsrf(:,is_ter))
      WHERE(pctsrf(:,is_sic)+pctsrf(:,is_oce)<EPSFRA) sicfra(:)=0.
      WHERE(sicfra<0.7) Ampl(:)=1.
      WHERE((sicfra>0.7).and.(sicfra<0.9)) Ampl=((sicfra-0.7)/0.2)*Ampl+((0.9-sicfra)/0.2)
      fluxsens=Ampl*fluxsens
      dflux_s=Ampl*dflux_s
    endif


! - 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)  

#ifdef ISO     
    CALL calcul_iso_surf_oce_vectall(klon, 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
!          write(*,*) 'ocean_forced_mod 176: sortie de ocean_forced_noice'
    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)
      ENDDO ! DO j=1,knon
    ENDIF !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, &
!GG       ps, u1, v1, gustiness, &
       ps, u1, v1, gustiness, pctsrf, &
!GG
       radsol, snow, qsol, agesno, tsoil, &
       qsurf, alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
!GG       tsurf_new, dflux_s, dflux_l, rhoa)
       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 &
!GG
#ifdef ISO
       ,xtprecip_rain, xtprecip_snow, xtspechum,Roce, &
       xtsnow, xtsol,xtevap,Rland_ice &  
#endif            
       )
!
! This subroutine treats the ocean where there is ice. 
! The routine reads data from climatologie file and does flux calculations at the 
! surface.
!
    USE dimphy
    USE geometry_mod, ONLY: longitude,latitude
    USE calcul_fluxs_mod
!GG    USE surface_data,     ONLY : calice, calsno
    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

    USE geometry_mod, ONLY: longitude,latitude,latitude_deg
!GG
    USE limit_read_mod
    USE fonte_neige_mod,  ONLY : fonte_neige
    USE indice_sol_mod
    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

!   INCLUDE "indicesol.h"


! Input arguments
!****************************************************************************************
    INTEGER, INTENT(IN)                  :: itime, jour, knon
    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
    REAL, INTENT(IN)                     :: dtime
    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
    REAL, DIMENSION(klon), INTENT(IN)    :: precip_rain, precip_snow
    REAL, DIMENSION(klon), INTENT(IN)    :: temp_air, spechum
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
    REAL, DIMENSION(klon), INTENT(IN)    :: ps
    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
    real, intent(in):: rhoa(:) ! (knon) density of moist air  (kg / m3)
!GG
    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
    REAL, DIMENSION(klon,nbsrf), INTENT(IN)  :: pctsrf
!GG
#ifdef ISO
    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtprecip_rain, xtprecip_snow
    REAL, DIMENSION(ntiso,klon), INTENT(IN) :: xtspechum
    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Roce 
    REAL, DIMENSION(niso,klon),  INTENT(IN) :: Rland_ice
#endif

! In/Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow, qsol
    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
    REAL, DIMENSION(klon, nsoilmx), INTENT(INOUT) :: tsoil
!GG
    REAL, DIMENSION(klon), INTENT(INOUT)          :: hice
    REAL, DIMENSION(klon), INTENT(INOUT)          :: tice
    REAL, DIMENSION(klon), INTENT(INOUT)          :: bilg_cumul
    REAL, DIMENSION(klon), INTENT(INOUT)          :: fcds
    REAL, DIMENSION(klon), INTENT(INOUT)          :: fcdi
    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_basal_growth
    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_basal_melt
    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_top_melt
    REAL, DIMENSION(klon), INTENT(INOUT)          :: dh_snow2sic
    REAL, DIMENSION(klon), INTENT(INOUT)          :: dtice_melt
    REAL, DIMENSION(klon), INTENT(INOUT)          :: dtice_snow2sic
!GG
#ifdef ISO     
    REAL, DIMENSION(niso,klon), INTENT(INOUT)     :: xtsnow
    REAL, DIMENSION(niso,klon), INTENT(IN)        :: xtsol
#endif

! Output arguments
!****************************************************************************************
    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens, fluxlat
    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l      
#ifdef ISO     
    REAL, DIMENSION(ntiso,klon), INTENT(OUT)      :: xtevap
#endif      

! Local variables
!****************************************************************************************
    LOGICAL                     :: check=.FALSE.
    INTEGER                     :: i, j
    REAL                        :: zfra
    REAL, PARAMETER             :: t_grnd=271.35
    REAL, DIMENSION(klon)       :: cal, beta, dif_grnd, capsol, icesub
    REAL, DIMENSION(klon)       :: alb_neig, tsurf_tmp
    REAL, DIMENSION(klon)       :: soilcap, soilflux
    REAL, DIMENSION(klon)       :: u0, v0
    REAL, DIMENSION(klon)       :: u1_lay, v1_lay
    REAL, DIMENSION(knon)       :: sens_prec_liq, sens_prec_sol
    REAL, DIMENSION(klon)       :: lat_prec_liq, lat_prec_sol    
!GG
    INTEGER                     :: ki
    INTEGER                     :: cpl_pas
    REAL, DIMENSION(klon)       :: bilg, fsic, 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, SAVE :: 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
! new (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)

! HF from ocean below ice
!    REAL, PARAMETER :: fseaN=2.0 ! NH
!    REAL, PARAMETER :: fseaS=4.0 ! SH    
!GG

#ifdef ISO
    REAL, PARAMETER :: t_coup = 273.15
    REAL, DIMENSION(klon) :: fq_fonte_diag
    REAL, DIMENSION(klon) :: fqfonte_diag
    REAL, DIMENSION(klon) :: snow_evap_diag 
    REAL, DIMENSION(klon) :: fqcalving_diag 
    REAL, DIMENSION(klon) :: run_off_lic_diag
    REAL :: coeff_rel_diag
    REAL :: max_eau_sol_diag  
    REAL, DIMENSION(klon) :: runoff_diag    
    INTEGER IXT
    REAL, DIMENSION(niso,klon) :: xtsnow_prec, xtsol_prec
    REAL, DIMENSION(klon) :: snow_prec, qsol_prec  
#endif

!****************************************************************************************
! Start calculation
!****************************************************************************************
    IF (check) WRITE(*,*)'Entering surface_seaice, knon=',knon 

!****************************************************************************************
! 1) 
! Flux calculation : tsurf_new, evap, fluxlat, fluxsens, flux_u1, flux_v1
!                    dflux_s, dflux_l and qsurf
!****************************************************************************************

    tsurf_tmp(:) = tsurf_in(:)

!GG
    IF (iflag_seaice==0) THEN ! Old LMDZ sea ice surface
!GG
! 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
       CALL soil(dtime, is_sic, knon, snow, tsurf_tmp, qsol, &
        & longitude(knindex(1:knon)), latitude(knindex(1:knon)), 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 
    ENDIF

!GG
    ELSEIF (iflag_seaice==2) 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)

    bilg(:)=0.
    dif_grnd(:)=0.
    beta(:) = 1.
    fsic(:) = pctsrf(:,is_sic)
    cpl_pas =  NINT(86400./dtime * 1.0) ! une fois par jour

    ! Surface, snow-ice and ice-ocean fluxes.
! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd)

    !write(*,*) 'radsol 1',radsol(1:100)
    DO i=1,knon
        ki=knindex(i)
        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(ki)=f_swpen-f_cond_s
    END DO

    endif

!    beta = 1.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(:)
    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)
    enddo

! - 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)  

!****************************************************************************************
! 2)
! Calculations due to snow and runoff
!
!****************************************************************************************
!GG
    if (iflag_seaice==0) then
!GG
#ifdef ISO
   ! verif
#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)
        ENDIF !IF ((snow(i) > ridicule)) THEN
      ENDIF !IF (iso_eau > 0) THEN      
    ENDDO !DO i=1,knon  
#endif
   ! end verif

    DO i = 1, knon
      snow_prec(i) = snow(i)
      DO ixt = 1, niso
      xtsnow_prec(ixt,i) = xtsnow(ixt,i)
      ENDDO !DO ixt=1,niso
      ! initialisation:
      fq_fonte_diag(i) = 0.0
      fqfonte_diag(i)  = 0.0
      snow_evap_diag(i)= 0.0
    ENDDO !DO i=1,knon
#endif


    CALL fonte_neige( knon, is_sic, knindex, dtime, &
         tsurf_tmp, precip_rain, precip_snow, &
         snow, qsol, tsurf_new, evap, icesub &
#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
! isotopes: tout est externalisé
!#ifdef ISOVERIF
!        write(*,*) 'ocean_forced_mod 377: call calcul_iso_surf_sic_vectall'
!        write(*,*) 'klon,knon=',klon,knon
!#endif
    CALL calcul_iso_surf_sic_vectall(klon,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
        !write(*,*) 'ocean_forced_mod 391: sortie calcul_iso_surf_sic_vectall'
    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)
      ENDDO ! DO j=1,knon
    ENDIF !IF (iso_eau > 0) then
#endif 
!#ifdef ISOVERIF
#endif   
!#ifdef ISO
    
! Calculation of albedo at snow (alb_neig) and update the age of snow (agesno)
! 
    CALL albsno(klon, 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)
    ENDDO

    alb2_new(:) = alb1_new(:)

!GG
  else 

        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
           ENDIF

           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)
           ENDIF
           ! 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(ki) = bilg(ki)+f_cond_i
        END DO

!****************************************************************************************
! 2) Update snow and ice surface : thickness 
!****************************************************************************************

    IF (iflag_seaice==1) THEN
!   Read from limit
    CALL limit_read_hice(knon,knindex,hice)
    ENDIF
!   Formula Krinner et al. 1997 : h = (0.2 + 3.8(f_min**2))(1 + 2(f- f_min)) 



    DO i=1,knon
        ki=knindex(i)
        IF (precip_snow(i) > 0.) THEN
            snow(i) = snow(i)+precip_snow(i)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
        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)
           ENDIF
        ENDIF
! 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
! Bottom melt / grow
! bottom freeze if bottom flux (cond + oce-ice) <0
        IF (iflag_seaice==2) THEN
         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(ki) &
                   / (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(ki)*ice_den)))
           hice_i=hice(ki)
           hice(ki)=dhsic+hice(ki)
           e_melt=e_melt-fsic(ki)*dhsic
           IF (e_melt.GT.0.) THEN
           ! frac_mf fraction used for lateral increase
           dfsic=MIN(amax-fsic(ki),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(ki) * 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(ki) &
                   / (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(ki)*ice_den)))
           hice(ki)=hice(ki)-dhsic
           e_melt=e_melt-fsic(ki)*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(ki)))
           hice(ki)=hice(ki)-dhsic
           e_melt=e_melt-fsic(ki)*dhsic*ice_den
           ! rest used to decrease fraction (up to 0!)
           dfsic=MIN(fsic(ki),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==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(ki)-ice_frac_min,e_melt/(h_ice_min*ice_den)*fsic(ki))
              ! if all melted do nothing with remaining heat 
              e_melt=MAX(0.,e_melt*fsic(ki)-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==2) THEN
             dhsic=h_test/(sea_den-ice_den+sno_den)
             hice(ki)=hice(ki)+dhsic
            ENDIF
            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==2) THEN
              dh_snow2sic(ki)=dhsic/dtime
            END IF
            dtice_snow2sic(ki)=(tice(ki)-tice_i)/dtime
        END IF
    END DO

        !write(*,*) 'hice 2',hice(1:100)
        !write(*,*) 'tice 2',tice(1:100)

        iflag_sic_albedo=iflag_seaice_alb

!*******************************************************************************
! 3) cumulate ice-ocean fluxes, update tslab, lateral grow
!***********************************************o*******************************
    !cumul fluxes.
    bilg_cumul(:)=bilg_cumul(:)+bilg(:)/float(cpl_pas)
    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab
        bilg_cumul(:)=0.
    END IF ! coupling time

!   write(*,*) 'hice 3',hice(1:100)
!   write(*,*) 'tice 3',tice(1:100)
    !tests ice fraction 
    WHERE (fsic.LT.ice_frac_min)
        tice=t_melt
        hice=h_ice_min
    END WHERE

    !write(*,*) 'hice 4',hice(1:100)
    !write(*,*) 'tice 4',tice(1:100)

    endif

!****************************************************************************************
! 4) Compute sea-ice and snow albedo
!****************************************************************************************
    IF (iflag_seaice > 0) THEN
    SELECT CASE (iflag_sic_albedo)
      CASE(0)
! 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 tkickness and temp)
          alb_ice=MAX(0.0,0.13*LOG(100.*hice(ki))+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(:)

      CASE(1)
! New visible and IR albedos, dry / melting snow
! based on Toyoda et al, 2020
      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)
          ENDIF
      END DO

      CASE(2)
! LIM3 scheme. Uses clear sky / overcast value, with 50% clear sky
      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
            ENDIF
            ! 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)*ice_den)) * zlog
            ENDIF
            ! 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)
      CALL albsno(klon, 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)
      ENDDO
      alb2_new(:) = alb1_new(:)

      print*,'alb_neig=',alb_neig
      print*,'zfra=',zfra
      print*,'snow=',snow
      print*,'alb1_new=',alb1_new
      print*,'alb2_new=',alb2_new
    END SELECT
    END IF
! ------ End Albedo ----------

!GG
  END SUBROUTINE ocean_forced_ice

!************************************************************************
! 1D case
!************************************************************************
!  SUBROUTINE read_tsurf1d(knon,sst_out)
!
! This subroutine specifies the surface temperature to be used in 1D simulations
!
!      USE dimphy, ONLY : klon
!
!      INTEGER, INTENT(IN)                  :: knon     ! nomber of points on compressed grid
!      REAL, DIMENSION(klon), INTENT(OUT)   :: sst_out  ! tsurf used to force the single-column model
!
!       INTEGER :: i
! COMMON defined in lmdz1d.F:
!       real ts_cur
!       common /sst_forcing/ts_cur
!
!       DO i = 1, knon
!        sst_out(i) = ts_cur
!       ENDDO
!
!      END SUBROUTINE read_tsurf1d
!
!
!************************************************************************
END MODULE ocean_forced_mod






