Changeset 3100

Timestamp:

Oct 25, 2023, 10:34:54 AM (13 months ago)

Author:

bhatnags

Message:

Generic-PCM

This commit updates the slab ocean module to a parallelisable dynamic slab ocean module. This is particularly relevant if you want to be able to use oceanic heat transport in parallel mode.

It has the following features:

(a) Computes sea ice creation and evolution.
(b) Snow has thermodynamic properties.
(c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
(d) Can be used in parallel mode.

Required callphys.def flags:
The slab ocean and its dependencies can be activated with the following flags (already added to deftank):
## Ocean options
## ~
# Model slab-ocean (Main flag for slab ocean)
ok_slab_ocean = .true.
# The following flags can only be set to true if ok_slab_ocean is true
# Ekman transport
slab_ekman = .true.
# Ekman zonal advection
slab_ekman_zonadv = .true.
# Horizontal diffusion (default coef_hdiff=25000., can be changed)
slab_hdiff = .true.
# Slab-ocean timestep (in physics timesteps)
cpl_pas = 1
# Gent-McWilliams? Scheme (can only be true if slab_ekman is true)
slab_gm = .true.

Notes:
In the current state, the model crashes if moistadjustment = .true. Unsure whether this is due to the slab or is an inherent issue with moistadj (under investigation).

SB and EM

Location:

trunk/LMDZ.GENERIC

Files:

• deftank/callphys.earthslab (moved) (moved from trunk/LMDZ.GENERIC/deftank/callphys.earth) (10 diffs)
• deftank/run.def.earthslab.64x48x30 (added)
• libf/aeronostd/photolysis_online.F (modified) (1 diff)
• libf/dynphy_lonlat/phystd/ini_archive.F (modified) (3 diffs)
• libf/dynphy_lonlat/phystd/iniphysiq_mod.F90 (modified) (2 diffs)
• libf/dynphy_lonlat/phystd/lect_start_archive.F (modified) (6 diffs)
• libf/dynphy_lonlat/phystd/newstart.F (modified) (1 diff)
• libf/dynphy_lonlat/phystd/start2archive.F (modified) (5 diffs)
• libf/dynphy_lonlat/phystd/write_archive.F (modified) (2 diffs)
• libf/phystd/calcul_fluxs_mod.F90 (added)
• libf/phystd/hydrol.F90 (modified) (10 diffs)
• libf/phystd/inifis_mod.F90 (modified) (1 diff)
• libf/phystd/iostart.F90 (modified) (6 diffs)
• libf/phystd/ocean_slab_mod.F90 (modified) (2 diffs)
• libf/phystd/phyetat0_mod.F90 (modified) (5 diffs)
• libf/phystd/phyredem.F90 (modified) (3 diffs)
• libf/phystd/phys_state_var_mod.F90 (modified) (4 diffs)
• libf/phystd/physiq_mod.F90 (modified) (34 diffs)
• libf/phystd/slab_heat_transp_mod.F90 (added)
• libf/phystd/slab_ice_h.F90 (deleted)
• libf/phystd/surf_heat_transp_mod.F90 (deleted)

trunk/LMDZ.GENERIC/deftank/callphys.earthslab

@@ -4 +4 @@
 tracer    = .true.
 # Diurnal cycle ?  if diurnal=false, diurnally averaged solar heating
-diurnal   = .false.
+diurnal   = .true.
 # Seasonal cycle ? if season=false, Ls stays constant, to value set in "start"
 season    = .true. 
@@ -17 +17 @@
 # Test energy conservation of model physics ?
 enertest  = .false.
-# Check to see if cpp, mugaz values used match gas mixture defined in gases.def (recommended) ?
-check_cpp_match=.true.
+# check if cpp and mugaz from start.nc are consistent with values computed by comp_cpp_mugaz with gases.def
+check_cpp_match = .false.
+# Check if physics inputs and outputs are ok
+check_physics_inputs = .true.
+check_physics_outputs = .true.
+
+## Directory where external input files are
+## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+datadir     = datadir
 
 ## Radiative transfer options
@@ -26 +33 @@
 # the rad. transfer is computed every "iradia" physical timestep
 iradia     = 6
+# Use blackbody for stellar spectrum ?
+stelbbody    = .false.
+# Stellar blackbody temperature ?
+stelTbb      = 2000.000
 # call multilayer correlated-k radiative transfer ?
 corrk      = .true.
-# Include continuum absorption in radiative transfer (note CO2 is treated separately) ?
-continuum  = .true.
 # folder in which correlated-k data is stored ?
-corrkdir   = N2CO2poor_H2Ovar
+corrkdir = Earth_JL13_extend
+# corrkdir    = Earth_110-710K
+# corrkdir   = N2CO2poor_H2Ovar
+# corrkdir   = megaCO2
+# corrkdir   = null
 # call visible gaseous absorption in radiative transfer ?
 callgasvis = .true.
+# call continuum in radiative transfer ?
+Continuum = .true.
 # Include Rayleigh scattering in the visible ?
 rayleigh   = .true.
+# Uniform absorption coefficient in radiative transfer?
+graybody     = .false.
+# Constant absorption coefficient in visible 
+#      (in m^2/kg; only if graybody=true):
+#      tau_surf= kappa*P/g
+kappa_VI = 1.e-4
+# Constant absorption coefficient in IR
+#      (in m^2/kg; only if graybody=true):
+kappa_IR = 5.e-1 
+# Use Newtonian cooling in place of radiative transfer ?
+newtonian    = .false.
+# Radiative timescale for Newtonian cooling ? [only if newtonian = T]
+tau_relax    = 30.00000
+# Test physics timescale in 1D ?
+testradtimes = .false.
 # Characteristic planetary equilibrium (black body) temperature
 # This is used only in the aerosol radiative transfer setup. (see aerave.F)
 tplanet    = 215.
-# Output spectral OLR in 3D?
+# Output spectral OLR in 1D/3D?
 specOLR    = .false.
 # Output global radiative balance in file 'rad_bal.out' - slow for 1D!!
-meanOLR    = .true.
+meanOLR    = .false.
 # Variable gas species: Radiatively active ?
 varactive  = .true.
 # Variable gas species: Fixed vertical distribution ?
+#   (not to be used in time integration mode)
 varfixed   = .false.
 # Variable gas species: Saturation percentage value at ground ?
-satval     = 0.0
+satval     = .8
+# Use fixed vertical profile, 1 step, no iteration ?
+kastprof     = .false.
+# Remove lower boundary (e.g. for gas giant sims)
+noradsurf    = .false. 
 
 ## Star type
@@ -60 +95 @@
 #       startype = 3            GJ644
 #       startype = 4            HD128167
-#       startype = 9            TRAPPIST-1
-#       startype = 10           Proxima Centauri
 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Stellar flux at 1 AU. Examples:
@@ -67 +100 @@
 # 1024.5 W m-2          Sol today x 0.75 = weak early Sun
 # 18.462 W m-2          The feeble Gl581
+# 19.960 W m-2          Gl581 with e=0.38 orbital average
 Fat1AU = 1366.0
+
 
 ## Tracer and aerosol options 
 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~
-# Gravitational sedimentation of tracers (KEEP FALSE FOR NOW) ?
-sedimentation = .false.
+# Gravitational sedimentation of tracers (just H2O ice for now) ?
+sedimentation = .true.
 
 ## Other physics options
@@ -92 +127 @@
 ## ~~~~~~~~~~~~~~~~~~~~~~~~~~
 # Number of radiatively active aerosols
-naerkind=1
+naerkind = 1
+# Varying H2O cloud fraction?
+CLFvarying    = .true.
+# H2O cloud fraction?
+CLFfixval     = 1.
+# number mixing ratio of CO2 ice particles
+Nmix_co2      = 1.e5
+# basic dust opacity
+dusttau       = 0.0
+# water cloud pressure level (norm. by psurf)
+cloudlvl      = 0.0
+# atm mass update due to tracer evaporation/condensation?
+mass_redistrib = .true.
 # Radiatively active CO2 aerosol?
 aeroco2       = .false.
@@ -101 +148 @@
 # Fixed water aerosol distribution?
 aerofixh2o  = .false.
-# basic dust opacity
-dusttau       = 0.0
-# Varying H2O cloud fraction?
-CLFvarying    = .true.
-# H2O cloud fraction if fixed?
-CLFfixval     = 0.5
-# fixed radii for cloud particles?
-radfixed=.false.
-# number mixing ratio of CO2 ice particles
-Nmix_co2      = 100000.
+# Radiatively active sulfur aersol?
+aeroh2so4     = .false.
+# fixed radii for h2o cloud particles?
+radfixed=.true.
 # number mixing ratio of water particles (for rafixed=.false.)
-Nmix_h2o      = 1.e7
+Nmix_h2o      = 4e6
 # number mixing ratio of water ice particles (for rafixed=.false.)
-Nmix_h2o_ice      = 5.e5
-# radius of H2O water particles (for rafixed=.true.):
-rad_h2o=10.e-6
+Nmix_h2o_ice      = 2.e4
+# radius of H2O water particles (for rafixed=.true.): (CHANGED FROM 10 TO 12 AFTER BENJAMIN)
+rad_h2o=12.e-6
 # radius of H2O ice particles (for rafixed=.true.):
 rad_h2o_ice=35.e-6
-# atm mass update due to tracer evaporation/condensation?
-mass_redistrib = .false.
 
 ## Water options 
@@ -130 +169 @@
 # Model water precipitation (including coagulation etc.)
 waterrain     = .true.
+# Moist adjustment
+moistadjustment = .true.
 # Use simple precipitation scheme?
 precip_scheme=4
 # multiplicative constant in Boucher 95 precip scheme
-Cboucher=1.
-# WATER: hydrology ?
+Cboucher=0.6
+# Include hydrology ?
 hydrology     = .true.
-# Spectral Dependant Albedo ?
-albedo_spectral_mode=.false.
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-# If albedo_spectral_mode=.true., albedosnow becomes the 0.5 micron snow albedo.
-#
-# albedosnow = 0.95  (0.73 Sun-integrated) for fresh snow.
-#            = 0.50  (0.39 Sun-integrated) for dirty snow.
-#            = 0.78  (0.60 Sun-integrated) for 'realistic' snow.
-# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+# active runoff ?
+activerunoff  = .true.
 # H2O snow (and ice) albedo ?
-albedosnow    = 0.60
+albedosnow    = 0.65
 # Maximum sea ice thickness ?
 maxicethick   = 10.
@@ -153 +187 @@
 # Evolve surface water sources ?
 sourceevol    = .false.
+# compute lightning rate ?
+compute_lightning     = .true.
 
 ## CO2 options 
 ## ~~~~~~~~~~~
-# Co2 ice albedo ?
-albedoco2ice   = 0.5
-# gas is non-ideal CO2 ?
-nonideal      = .false.
 # call CO2 condensation ?
 co2cond       = .false.
@@ -165 +197 @@
 nearco2cond   = .false.
 
+## Subsurface options
+## ~~~~~~~~~~~~~~~~~~
+# Number of subsurface layers (For Earth 14 is OK; for Mars 18)
+nsoilmx=14
+# Thickness of topmost soil layer (m) (For Earth 0.03 is OK; for Mars 0.0002)
+lay1_soil=0.03
+#Coefficient for soil layer thickness distribution (default is 2)
+alpha_soil=2
+
+## Ocean options
+## ~~~~~~~~~~~~~
+# Model slab-ocean (Main flag for slab ocean)
+ok_slab_ocean = .true.
+# The following flags can only be set to true if ok_slab_ocean is true
+# Ekman transport
+slab_ekman = .true.
+# Ekman zonal advection
+slab_ekman_zonadv = .true.
+# Horizontal diffusion (default coef_hdiff=25000., can be changed)
+slab_hdiff = .true.
+# Slab-ocean timestep (in physics timesteps)
+cpl_pas = 1
+# sea-ice
+#ok_slab_sic   = .true.
+# Gent-McWilliams Scheme (can only be true if slab_ekman is true)
+slab_gm   = .true.
+# Slab convective adjustment? 0 - no, 1 - yes
+slab_cadj = 1
+# H2O snow (and ice) albedo for sea ?
+#albedosnowsea    = 0.65
+

trunk/LMDZ.GENERIC/libf/aeronostd/photolysis_online.F

@@ -1172 +1172 @@
 
       use chimiedata_h,  only: albedo_snow_chim, albedo_co2_ice_chim
-      use slab_ice_h,    only: h_alb_ice, 
-     &                         alb_ice_min, alb_ice_max
-      use tracer_h,      only: igcm_h2o_ice, igcm_co2_ice
+!      use slab_ice_h,    only: h_alb_ice, alb_ice_min, alb_ice_max
+      use ocean_slab_mod, only: h_alb_ice
+      use ocean_slab_mod, only: alb_ice_min
+      use ocean_slab_mod, only: alb_ice_max
       use callkeys_mod,  only: ok_slab_ocean, co2cond, alb_ocean
       use phys_state_var_mod, only: albedo_bareground, 

trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/ini_archive.F

@@ -36 +36 @@
 
       USE comsoil_h
-      USE slab_ice_h, only: noceanmx
+!      USE slab_ice_h, only: noceanmx
+      USE ocean_slab_mod, ONLY: nslay
 !      use control_mod
       USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
@@ -81 +82 @@
       INTEGER idim_tim
       INTEGER idim_nsoilmx ! "subsurface_layers" dimension ID #
-      INTEGER idim_noceanmx ! "ocean_layers" dimension ID #
+      INTEGER idim_nslay ! "ocean_layers" dimension ID #
       INTEGER nid,nvarid
       real sig_s(llm),s(llm)
@@ -164 +165 @@
       ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
       ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
-      ierr = NF_DEF_DIM (nid,"ocean_layers",noceanmx,idim_noceanmx)
+      ierr = NF_DEF_DIM (nid,"ocean_layers",nslay,idim_nslay)
 
       ierr = NF_DEF_DIM (nid,"index", length, idim_index)

trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/iniphysiq_mod.F90

@@ -11 +11 @@
 
 use control_mod, only: nday 
-use surf_heat_transp_mod, only: ini_surf_heat_transp
+!use surf_heat_transp_mod, only: ini_surf_heat_transp
+USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
 use infotrac, only : nqtot ! number of advected tracers
 use planete_mod, only: ini_planete_mod
@@ -83 +84 @@
 call getin_p("ok_slab_ocean",ok_slab_ocean)
 if (ok_slab_ocean) then
-  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
-                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
-                            aireu,airev)
+!  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
+!                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
+!                            aireu,airev)
+   CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,&
+                                  cu,cuvsurcv,cv,cvusurcu, &
+                                  aire,apoln,apols, &
+                                  aireu,airev,rlatvdyn)
 endif
 

trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/lect_start_archive.F

@@ -9 +9 @@
       USE tracer_h, ONLY: igcm_co2_ice
       USE infotrac, ONLY: tname, nqtot
-      USE slab_ice_h, ONLY: noceanmx
+!      USE slab_ice_h, ONLY: noceanmx
+      USE ocean_slab_mod, ONLY: nslay
       USE callkeys_mod, ONLY: ok_slab_ocean
       USE comvert_mod, ONLY: ap,bp,aps,bps,preff
@@ -68 +69 @@
       REAL,INTENT(OUT) :: emis(ngrid)
       REAL,INTENT(OUT) :: q2(ngrid,llm+1),qsurf(ngrid,nqtot)
-      REAL,INTENT(OUT) :: tslab(ngrid,noceanmx)
+      REAL,INTENT(OUT) :: tslab(ngrid,nslay)
       REAL ,INTENT(OUT) ::rnat(ngrid),pctsrf_sic(ngrid)
       REAL,INTENT(OUT) :: tsea_ice(ngrid),sea_ice(ngrid)
@@ -98 +99 @@
       real emisS(iip1,jjp1)
       REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
-      real tslabS(iip1,jjp1,noceanmx)
+      real tslabS(iip1,jjp1,nslay)
       real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
       real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
@@ -287 +288 @@
       allocate(mlayerold(nsoilold))
       allocate(qsurfold(imold+1,jmold+1,nqtot))
-      allocate(tslabold(imold+1,jmold+1,noceanmx))
+      allocate(tslabold(imold+1,jmold+1,nslay))
       allocate(rnatold(imold+1,jmold+1))
       allocate(pctsrf_sicold(imold+1,jmold+1))
@@ -684 +685 @@
       if(ok_slab_ocean) then
       start=(/1,1,1,memo/)
-      count=(/imold+1,jmold+1,noceanmx,1/)
+      count=(/imold+1,jmold+1,nslay,1/)
 
        ierr=NF_INQ_VARID(nid,"tslab",nvarid)
@@ -1298 +1299 @@
 c 6.3   Slab Ocean :
 c-----------------------------------------------------------------------
-      call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,noceanmx,
-     &                   rlonuold,rlatvold,rlonu,rlatv)
-      call gr_dyn_fi (noceanmx,iim+1,jjm+1,ngrid,tslabs,tslab)
+      call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,nslay,
+     &                   rlonuold,rlatvold,rlonu,rlatv)
+      call gr_dyn_fi (nslay,iim+1,jjm+1,ngrid,tslabs,tslab)
 
       call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,

trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/newstart.F

@@ -30 +30 @@
       use phyredem, only: physdem0, physdem1
       use iostart, only: open_startphy
-      use slab_ice_h, only:noceanmx
+!      use slab_ice_h, only:noceanmx
+      USE ocean_slab_mod, ONLY: nslay
       use filtreg_mod, only: inifilr
       USE mod_const_mpi, ONLY: COMM_LMDZ

trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/start2archive.F

@@ -22 +22 @@
       USE comsoil_h
       
-      use slab_ice_h, only: noceanmx
+!      use slab_ice_h, only: noceanmx
+      USE ocean_slab_mod, ONLY: nslay
       USE ioipsl_getincom, only: getin
       USE planete_mod, only: year_day
@@ -84 +85 @@
 !     added by BC for slab ocean
       REAL rnat(ngridmx),pctsrf_sic(ngridmx),sea_ice(ngridmx)
-      REAL tslab(ngridmx,noceanmx),tsea_ice(ngridmx)
+      REAL, ALLOCATABLE :: tslab(:,:)
+      REAL tsea_ice(ngridmx)
 
 
@@ -104 +106 @@
 !     added by BC for slab ocean
       REAL rnatS(ip1jmp1),pctsrf_sicS(ip1jmp1),sea_iceS(ip1jmp1)
-      REAL tslabS(ip1jmp1,noceanmx),tsea_iceS(ip1jmp1)
+      REAL, ALLOCATABLE :: tslabS(:,:)
+      REAL tsea_iceS(ip1jmp1)
 
 !     For non-orographic GW
@@ -156 +159 @@
       allocate(qsurf(ngridmx,nqtot))
       allocate(qsurfS(ip1jmp1,nqtot))
+      allocate(tslab(ngridmx,nslay)) !Added by SB for slab ocean
+      allocate(tslabS(ip1jmp1,nslay)) !Added by SB for slab ocean
 ! other array allocations:
 !      call ini_comsoil_h(ngridmx) ! done via iniphysiq
@@ -356 +361 @@
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsea_ice,tsea_iceS)
       call gr_fi_dyn(1,ngridmx,iip1,jjp1,sea_ice,sea_iceS)
-      call gr_fi_dyn(noceanmx,ngridmx,iip1,jjp1,tslab,tslabS)
+      call gr_fi_dyn(nslay,ngridmx,iip1,jjp1,tslab,tslabS)
 
       call gr_fi_dyn(llm,ngridmx,iip1,jjp1,du_nonoro_gwd,du_nonoro_gwdS)

trunk/LMDZ.GENERIC/libf/dynphy_lonlat/phystd/write_archive.F

@@ -33 +33 @@
 
       use comsoil_h, only: nsoilmx
-      use slab_ice_h, only: noceanmx
+!      use slab_ice_h, only: noceanmx
+      USE ocean_slab_mod, ONLY: nslay
 
       implicit none
@@ -195 +196 @@
         edges(1)=iip1
         edges(2)=jjp1
-        edges(3)=noceanmx
+        edges(3)=nslay
         edges(4)=1
 #ifdef NC_DOUBLE

trunk/LMDZ.GENERIC/libf/phystd/hydrol.F90

@@ -7 +7 @@
 
   use ioipsl_getin_p_mod, only: getin_p
+#ifndef MESOSCALE
+  use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
+  use mod_phys_lmdz_para, only : is_master, gather, scatter
+#endif
   use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
   USE surfdat_h
@@ -12 +16 @@
   USE geometry_mod, only: cell_area
   USE tracer_h
-  use slab_ice_h
+!  use slab_ice_h
+  USE ocean_slab_mod, ONLY: alb_ice_max
+  USE ocean_slab_mod, ONLY: alb_ice_min
+  USE ocean_slab_mod, ONLY: h_alb_ice
   use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
   use radinc_h, only : L_NSPECTV
@@ -100 +107 @@
 !$OMP THREADPRIVATE(ivap,iliq,iice)
 
-      logical, save :: firstcall
+      logical, save :: firstcall=.true.
 !$OMP THREADPRIVATE(firstcall)
 
-      data firstcall /.true./
+     real :: runoffamount(ngrid)
+!#ifdef CPP_PARA
+      real :: runoffamount_glo(klon_glo)
+      real :: zqsurf_iliq_glo(klon_glo)
+      real :: rnat_glo(klon_glo)
+      real :: oceanarea_glo
+      real :: cell_area_glo(klon_glo)
+!#else
+!      real :: runoffamount_glo(ngrid)
+!      real :: zqsurf_iliq_glo(ngrid)
+!#endif
 
 
@@ -137 +154 @@
 !                      Soon to be extended to the entire water cycle...
 
-!     Total ocean surface area
+!     LOCAL ocean surface area
          oceanarea=0.
          do ig=1,ngrid
@@ -160 +177 @@
       endif
 
+!       call writediagfi(ngrid,"hydrol_rnat"," "," ",2,rnat)
+!!       write (*,*) "oceanarea", oceanarea
+
 !     add physical tendencies already calculated
 !     ------------------------------------------
@@ -170 +190 @@
          enddo    
       enddo
+
+!      call writediagfi(ngrid,"hydrol_zqsurf_iliq_1"," "," ",2,zqsurf(1:ngrid,iliq))
  
       do ig=1,ngrid
@@ -294 +316 @@
 !     deal with runoff
             if(activerunoff)then
-
                runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
                if(ngrid.gt.1)then ! runoff only exists in 3D
@@ -331 +352 @@
       end do ! ig=1,ngrid
 
+!      call writediagfi(ngrid,"hydrol_zqsurf_iliq_2"," "," ",2,zqsurf(1:ngrid,iliq))
+
 !     perform crude bulk averaging of temperature in ocean
 !     ----------------------------------------------------
@@ -363 +386 @@
       if(activerunoff)then
 
-         totalrunoff=0.
+!         totalrunoff=0.
          do ig=1,ngrid
-            if (nint(rnat(ig)).eq.1) then
-               totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
-            endif
+            runoffamount(ig) = cell_area(ig)*runoff(ig)
+!            if (nint(rnat(ig)).eq.1) then
+!               totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
+!            endif
          enddo
+
+!    collect on the full grid
+        call gather(runoffamount,runoffamount_glo)
+        call gather(zqsurf(1:ngrid,iliq),zqsurf_iliq_glo)
+        call gather(rnat,rnat_glo)
+        call gather(cell_area,cell_area_glo)
+
+        if (is_master) then
+                totalrunoff=0.
+                oceanarea_glo=0.
+                do ig=1,klon_glo
+                   if (nint(rnat_glo(ig)).eq.1) then
+                        totalrunoff = totalrunoff + runoffamount_glo(ig)
+                   endif
+                   if (nint(rnat_glo(ig)).eq.0) then
+                        oceanarea_glo = oceanarea_glo + cell_area_glo(ig)
+                   endif
+                enddo
+
+                do ig=1,klon_glo
+                   if (nint(rnat_glo(ig)).eq.0) then
+                        zqsurf_iliq_glo(ig) = zqsurf_iliq_glo(ig) + &
+                    totalrunoff/oceanarea_glo
+                   endif
+                enddo
+
+        endif! is_master
+
+!    scatter the field back on all processes
+        call scatter(zqsurf_iliq_glo,zqsurf(1:ngrid,iliq))
+
+        
          
-         do ig=1,ngrid
-            if (nint(rnat(ig)).eq.0) then
-               zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
-                    totalrunoff/oceanarea
-            endif
-         enddo
-
-      endif         
-
+!         do ig=1,ngrid
+!            if (nint(rnat(ig)).eq.0) then
+!               zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
+!                    totalrunoff/oceanarea
+!            endif
+!         enddo
+
+      endif !activerunoff         
+
+!      call writediagfi(ngrid,"hydrol_zqsurf_iliq_3"," "," ",2,zqsurf(1:ngrid,iliq))
 
 !     Re-add the albedo effects of CO2 ice if necessary
@@ -400 +457 @@
       enddo
 
+!      call writediagfi(ngrid,"hydrol_dqs_hyd_iliq"," "," ",2,dqs_hyd(1:ngrid,iliq))  
+
       if (activerunoff) then
         call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)

trunk/LMDZ.GENERIC/libf/phystd/inifis_mod.F90

@@ -446 +446 @@
      if (is_master) write(*,*) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean
      ! Sanity check: for now slab oncean only works in serial mode
-     if (ok_slab_ocean.and.is_parallel) then
-       call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
-     endif
-
-     if (is_master) write(*,*) trim(rname)//": Use slab-sea-ice ?"
-     ok_slab_sic=.true.         ! default value
-     call getin_p("ok_slab_sic",ok_slab_sic)
-     if (is_master) write(*,*) trim(rname)//": ok_slab_sic = ",ok_slab_sic
-
-     if (is_master) write(*,*) trim(rname)//&
-       ": Use heat transport for the ocean ?"
-     ok_slab_heat_transp=.true.   ! default value
-     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
-     if (is_master) write(*,*) trim(rname)//&
-       ": ok_slab_heat_transp = ",ok_slab_heat_transp
+!     if (ok_slab_ocean.and.is_parallel) then
+!       call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
+!     endif
+
+!     if (is_master) write(*,*) trim(rname)//": Use slab-sea-ice ?"
+!     ok_slab_sic=.true.         ! default value
+!     call getin_p("ok_slab_sic",ok_slab_sic)
+!     if (is_master) write(*,*) trim(rname)//": ok_slab_sic = ",ok_slab_sic
+
+!     if (is_master) write(*,*) trim(rname)//&
+!       ": Use heat transport for the ocean ?"
+!     ok_slab_heat_transp=.true.   ! default value
+!     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
+!     if (is_master) write(*,*) trim(rname)//&
+!       ": ok_slab_heat_transp = ",ok_slab_heat_transp
 
 ! Photochemistry and chemistry in the thermosphere

trunk/LMDZ.GENERIC/libf/phystd/iostart.F90

@@ -469 +469 @@
   USE tracer_h, only: nqtot
   USE comsoil_h, only: nsoilmx
-  USE slab_ice_h, only: noceanmx
+!  USE slab_ice_h, only: noceanmx
+  USE ocean_slab_mod, ONLY: nslay
 
   IMPLICIT NONE
@@ -558 +559 @@
       ENDIF
 
-      ierr=NF90_DEF_DIM(nid_restart,"ocean_layers",noceanmx,idim8)
+      ierr=NF90_DEF_DIM(nid_restart,"ocean_layers",nslay,idim8)
       IF (ierr/=NF90_NOERR) THEN
         write(*,*)'open_restartphy: problem defining oceanic layer dimension '
@@ -646 +647 @@
   USE mod_grid_phy_lmdz
   USE mod_phys_lmdz_para
-  USE slab_ice_h, only: noceanmx
+!  USE slab_ice_h, only: noceanmx
+  USE ocean_slab_mod, ONLY: nslay
+
   IMPLICIT NONE
   CHARACTER(LEN=*),INTENT(IN)    :: field_name
@@ -832 +835 @@
         endif ! of if (.not.present(time))
 
-      ELSE IF (field_size==noceanmx) THEN
+      ELSE IF (field_size==nslay) THEN
         ! input is a 2D "oceanic field" array
         if (.not.present(time)) then ! for a time-independent field
@@ -948 +951 @@
   USE comsoil_h, only: nsoilmx
   USE mod_phys_lmdz_para, only: is_master
-  USE slab_ice_h, only: noceanmx
+!  USE slab_ice_h, only: noceanmx
+  USE ocean_slab_mod, ONLY: nslay
   IMPLICIT NONE
      CHARACTER(LEN=*),INTENT(IN) :: var_name
@@ -1000 +1004 @@
         ! We know it is an  "mlayer" kind of 1D array
         idim1d=idim3
-      ELSEIF (var_size==noceanmx) THEN
+      ELSEIF (var_size==nslay) THEN
         ! We know it is an  "mlayer" kind of 1D array
         idim1d=idim8

trunk/LMDZ.GENERIC/libf/phystd/ocean_slab_mod.F90

@@ -1 @@
-!
-! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ocean_slab_mod.F90,v 1.3 2008-02-04 16:24:28 fairhead Exp $
-!
+!Completed
 MODULE ocean_slab_mod
 !
-! This module is used for both surface ocean and sea-ice when using the slab ocean,
-! "ocean=slab".
-!
-
-
-      use slab_ice_h
-      use watercommon_h, only: T_h2O_ice_liq
-      use surf_heat_transp_mod
-      implicit none
-
-
-
-
-
-  
-  !LOGICAL, PRIVATE, SAVE  :: ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS
-  !!$OMP THREADPRIVATE(ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS)
-  !INTEGER, PRIVATE, SAVE                           :: slab_ekman, slab_cadj
-  !!$OMP THREADPRIVATE(slab_ekman,slab_cadj)
-  INTEGER, PRIVATE, SAVE                           :: lmt_pas, julien, idayvrai
-  !$OMP THREADPRIVATE(lmt_pas,julien,idayvrai)
+!==================================================================
+!
+!     Purpose
+!     -------
+!     The dynamical slab ocean model of the Generic-PCM. It has the following features:
+!         (a) Computes sea ice creation and evolution.
+!         (b) Snow has thermodynamic properties.
+!         (c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
+!         (d) Can be used in parallel mode.
+!
+!     Authors
+!     -------
+!     S. Bhatnagar and E. Millour (2023)
+!     Adapted from the ocean modules of LMDZ Earth (F. Codron) and the Generic-PCM (B. Charnay, 2013).
+!
+!     Notes
+!     -----
+!     Compared to the old model, the new model has the following changes (non-exhaustive):
+!         (a) More realistic description of sea ice creation and evolution - simultaneous
+!             surface, side and bottom melting / freezing depending on fluxes.
+!         (b) Snow has an effective heat capacity.
+!         (c) Snow has "weight"; it can sink an ice block if there is too much of it.
+!         (d) Snow can be blown off by wind.
+!         (e) The two-layer ocean allows for convective adjustment.
+!         (f) Diffusion can follow the Gent-McWilliams scheme + Eddy diffusivity.
+!         (g) Can be used in parallel mode.
+!
+!==================================================================
+
+  USE dimphy, ONLY: klon
+  USE mod_grid_phy_lmdz, ONLY: klon_glo
+  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
+
+  IMPLICIT NONE
+  PRIVATE
+  PUBLIC :: ocean_slab_init, ocean_slab_ice, ocean_slab_noice, &
+            ocean_slab_frac, ocean_slab_get_vars, ocean_slab_final
+
+!***********************************************************************************
+! Global saved variables
+!***********************************************************************************
+  ! number of slab vertical layers
+  INTEGER, PUBLIC, SAVE :: nslay=2
+  !$OMP THREADPRIVATE(nslay)
+  ! number of oceanic grid points
+  INTEGER, PUBLIC, SAVE :: knon
+  !$OMP THREADPRIVATE(knon)
+  ! timestep for coupling (update slab temperature) in timesteps
   INTEGER, PRIVATE, SAVE                           :: cpl_pas
   !$OMP THREADPRIVATE(cpl_pas)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_seaice
-  !$OMP THREADPRIVATE(tmp_seaice)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: tmp_tslab_loc
-  !$OMP THREADPRIVATE(tmp_tslab_loc)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slab_bils
-  !$OMP THREADPRIVATE(slab_bils)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE , SAVE  :: lmt_bils
-  !$OMP THREADPRIVATE(lmt_bils)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tmp_pctsrf_slab
-  !$OMP THREADPRIVATE(tmp_pctsrf_slab)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: tmp_tslab
-  !$OMP THREADPRIVATE(tmp_tslab)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_radsol
-  !$OMP THREADPRIVATE(tmp_radsol)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: dt_hdiff
-  !$OMP THREADPRIVATE(dt_hdiff)
-  REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: dt_ekman
-  !$OMP THREADPRIVATE(dt_ekman)
-  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_flux_o, tmp_flux_g
-  !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
+  ! cyang = 1/heat capacity of top layer (rho.c.H)
+  REAL, PRIVATE, SAVE                              :: cyang
+  !$OMP THREADPRIVATE(cyang)
+  ! depth of slab layers (1st or 2nd layer)
   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
   !$OMP THREADPRIVATE(slabh)
-  LOGICAL, PRIVATE, SAVE                           :: check = .FALSE.
-  !$OMP THREADPRIVATE(check)
+  ! slab temperature
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
+  !$OMP THREADPRIVATE(tslab)
+  ! heat flux convergence due to Ekman
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
+  !$OMP THREADPRIVATE(dt_ekman)
+  ! heat flux convergence due to horiz diffusion
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
+  !$OMP THREADPRIVATE(dt_hdiff)
+  ! heat flux convergence due to GM eddy advection
+  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_gm
+  !$OMP THREADPRIVATE(dt_gm)
+  ! fraction of ocean covered by sea ice (sic / (oce+sic))
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
+  !$OMP THREADPRIVATE(fsic)
+  ! temperature of the sea ice
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
+  !$OMP THREADPRIVATE(tice)
+  ! sea ice thickness, in kg/m2
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
+  !$OMP THREADPRIVATE(seaice)
+  ! net surface heat flux, weighted by open ocean fraction
+  ! slab_bils accumulated over cpl_pas timesteps
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
+  !$OMP THREADPRIVATE(bils_cum)
+  ! net heat flux into the ocean below the ice : conduction + solar radiation
+  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
+  !$OMP THREADPRIVATE(slab_bilg)
+  ! slab_bilg cululated over cpl_pas timesteps
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
+  !$OMP THREADPRIVATE(bilg_cum)
+  ! wind stress saved over cpl_pas timesteps
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
+  !$OMP THREADPRIVATE(taux_cum)
+  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
+  !$OMP THREADPRIVATE(tauy_cum)
+
+!***********************************************************************************
+! Parameters (could be read in def file: move to slab_init)
+!***********************************************************************************
+! snow and ice physical characteristics:
+    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp [in K]
+    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp [in K]
+    REAL, PARAMETER :: sno_den=300. !mean snow density [in kg/m3]
+    REAL, PARAMETER :: ice_den=917. ! ice density [in kg/m3]
+    REAL, PARAMETER :: sea_den=1026. ! sea water density [in kg/m3]
+    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice [in W/(m.K) or (W.kg)/(K.m4)]
+    REAL, PRIVATE, SAVE :: sno_cond ! conductivity of snow [in W/(m.K) or (W.kg)/(K.m4)]
+    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice [in J/(kg.K)]
+    REAL, PARAMETER :: sea_cap=3994.   ! specific heat capacity, seawater [in J/(kg.K)]
+    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice [in J/kg]
+    REAL, PARAMETER :: ice_sub=2834000. ! latent heat of sublimation for snow and ice [in J/kg]
+
+! control of snow and ice cover & freeze / melt (heights in m converted to kg/m2)
+    REAL, PARAMETER :: snow_min=0.05*sno_den ! critical snow height [in kg/m2]
+    REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean [in kg/m2]
+    REAL, PARAMETER :: ice_frac_min=0.001
+    REAL, PRIVATE, SAVE :: ice_frac_max ! Max ice fraction (leads)
+    REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness [in kg/m2]
+    REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness [in kg/m2]
+    ! below ice_thin, priority is to melt lateral / grow height
+    ! ice_thin is also height of new ice
+    REAL, PRIVATE, SAVE :: h_ice_thick ! thin ice thickness
+    ! above ice_thick, priority is melt height / grow lateral
+    REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice [in kg/m2]
+    REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height [in kg/m2]
+
+    REAL, PARAMETER :: epsfra=1.0E-05 ! minimial grid fraction size below which there is no ice
+
+    REAL, PARAMETER, PUBLIC :: capcalocean=50.*4.228e+06 ! surface heat capacity [J.K-1.m-2] (assuming 50 m slab ocean)
+    REAL, PARAMETER, PUBLIC :: capcalseaice=5.1444e+06*0.15
+    REAL, PARAMETER, PUBLIC :: capcalsno=2.3867e+06*0.15
+
+    REAL, PARAMETER, PUBLIC :: h_alb_ice=0.5*ice_den ! height for full ice albedo
+    REAL, PARAMETER, PUBLIC :: h_sno_alb=0.02*sno_den ! height for control of snow fraction
+
+    REAL, PARAMETER, PUBLIC :: alb_ice_min=0.2
+    REAL, PARAMETER, PUBLIC :: alb_ice_max=0.65
+
+! Horizontal transport parameters
+   ! flag for horizontal diffusion
+   LOGICAL, PUBLIC, SAVE :: slab_hdiff
+   !$OMP THREADPRIVATE(slab_hdiff)
+   ! flag for GM eddy diffusivity
+   LOGICAL, PUBLIC, SAVE :: slab_gm
+   !$OMP THREADPRIVATE(slab_gm)
+   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
+   !$OMP THREADPRIVATE(coef_hdiff)
+   ! flags for Ekman, conv adjustment
+   LOGICAL, PUBLIC, SAVE :: slab_ekman
+   !$OMP THREADPRIVATE(slab_ekman)
+   INTEGER, PUBLIC, SAVE :: slab_cadj
+   !$OMP THREADPRIVATE(slab_cadj)
+
+!***********************************************************************************
 
 CONTAINS
 !
-!****************************************************************************************
-!
-  SUBROUTINE ocean_slab_init(ngrid,dtime, tslab_rst, seaice_rst, pctsrf_rst)
-
-      use slab_ice_h
-
-
-
-
-! Input variables
-!****************************************************************************************
-        
-    integer,intent(in) :: ngrid ! number of atmospherci columns
+!***********************************************************************************
+!
+  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst, tslab_rst, tice_rst, seaice_rst, zmasq)
+
+! This routine
+! (1) allocates variables initialised from restart fields
+! (2) allocates some other variables internal to the ocean module
+
+    USE ioipsl_getin_p_mod, ONLY : getin_p
+    USE mod_phys_lmdz_transfert_para, ONLY : gather
+    USE slab_heat_transp_mod, ONLY : ini_slab_transp
+
+    ! Input variables
+!***********************************************************************************
     REAL, INTENT(IN)                         :: dtime
 ! Variables read from restart file
-    REAL, DIMENSION(ngrid,noceanmx), INTENT(IN)  :: tslab_rst         
-    REAL, DIMENSION(ngrid), INTENT(IN)        :: seaice_rst
-    REAL, DIMENSION(ngrid), INTENT(IN) :: pctsrf_rst
-
+    REAL, DIMENSION(klon), INTENT(IN) :: pctsrf_rst
+    REAL, DIMENSION(klon,nslay), INTENT(IN) :: tslab_rst
+    REAL, DIMENSION(klon), INTENT(IN) :: tice_rst
+    REAL, DIMENSION(klon), INTENT(IN) :: seaice_rst
+    REAL, DIMENSION(klon), INTENT(IN) :: zmasq
 
 ! Local variables
-!****************************************************************************************
+!************************************************************************************
     INTEGER                :: error
+    REAL, DIMENSION(klon_glo) :: zmasq_glo
     CHARACTER (len = 80)   :: abort_message
-    CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
-
-
-    print*,'************************'
-    print*,'SLAB OCEAN est actif, prenez precautions !'
-    print*,'************************'    
-
-! Lecture des parametres:
-!    CALL getpar('ok_slab_sic',.true.,ok_slab_sic,'glace de mer dans slab')
-!    CALL getpar('slab_ekman',0,slab_ekman,'type transport Ekman pour slab (0=rien)')
-!    CALL getpar('slab_cadj',1,slab_cadj,'type ajustement conv slab 2 couches')
-
-! Allocate variables initialize from restart fields
-      ALLOCATE(tmp_tslab(ngrid,noceanmx), stat = error)
-      IF (error /= 0) THEN
-         abort_message='Pb allocation tmp_tslab'
-         CALL abort_physic(modname,abort_message,1)
-      ENDIF
-      tmp_tslab(:,:) = tslab_rst(:,:)
-      ALLOCATE(tmp_tslab_loc(ngrid,noceanmx), stat = error)
-      IF (error /= 0) THEN
-         abort_message='Pb allocation tmp_tslab_loc'
-         CALL abort_physic(modname,abort_message,1)
-      ENDIF
-      tmp_tslab_loc(:,:) = tslab_rst(:,:)
-
-    ALLOCATE(tmp_seaice(ngrid), stat = error)
+    CHARACTER (len = 20)   :: modname = 'ocean_slab_init'
+
+!***********************************************************************************
+! Define some parameters
+!***********************************************************************************
+!
+! cpl_pas  coupling period (update of tslab and ice fraction)
+! for a calculation at each physical timestep, cpl_pas=1
+    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
+    CALL getin_p('cpl_pas',cpl_pas)
+    print *,'cpl_pas',cpl_pas
+! Number of slab layers
+!    nslay=2
+!    CALL getin_p('slab_layers',nslay)
+    print *,'number of slab layers : ',nslay
+! Layer thickness
+    ALLOCATE(slabh(nslay), stat = error)
     IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_seaice'
+       abort_message='Pb allocation slabh'
        CALL abort_physic(modname,abort_message,1)
     ENDIF
-    tmp_seaice(:) = seaice_rst(:)
-
-    ALLOCATE(tmp_pctsrf_slab(ngrid), stat = error)
+    slabh(1)=50. ! Height of first ocean layer (wind-mixed layer)
+    CALL getin_p('slab_depth',slabh(1))
+    IF (nslay.GT.1) THEN
+        slabh(2)=150. ! Height of second ocean layer (deep ocean layer)
+    END IF
+! cyang = 1/heat capacity of top layer (rho.c.H)
+    cyang=1/(slabh(1)*sea_den*sea_cap)
+
+! ********** Sea Ice parameters ***********
+    ice_frac_max = 0.99 ! frac = 1 may lead to some problems.
+    CALL getin_p('ice_frac_max',ice_frac_max)
+    h_ice_thick = 1.5
+    CALL getin_p('h_ice_thick',h_ice_thick)
+    h_ice_thick = h_ice_thick * ice_den
+    sno_cond = 0.31
+    CALL getin_p('sno_cond',sno_cond)
+    sno_cond = sno_cond * sno_den
+
+! ********** Heat Transport parameters ****
+! Ekman transport
+!    slab_ekman=0
+    slab_ekman=.FALSE.
+    CALL getin_p('slab_ekman',slab_ekman)
+! GM eddy advection (2-layers only)
+    slab_gm=.FALSE.
+    CALL getin_p('slab_gm',slab_gm)
+!    IF (slab_ekman.LT.2) THEN
+    IF (.NOT.slab_ekman) THEN
+       slab_gm=.FALSE.
+    ENDIF
+! Horizontal diffusion
+    slab_hdiff=.FALSE.
+    CALL getin_p('slab_hdiff',slab_hdiff)
+    IF (slab_gm) THEN
+       coef_hdiff=8000. ! non-dimensional; coef_hdiff should be 25000 if GM is off
+    ELSE
+       coef_hdiff=25000.
+    ENDIF
+    CALL getin_p('coef_hdiff',coef_hdiff)
+
+! Convective adjustment
+!    IF (nslay.EQ.1) THEN
+!        slab_cadj=0
+!    ELSE
+        slab_cadj=1
+!    END IF
+    CALL getin_p('slab_cadj',slab_cadj)
+
+!************************************************************************************
+! Allocate surface fraction read from restart file
+!************************************************************************************
+    ALLOCATE(fsic(klon), stat = error)
     IF (error /= 0) THEN
        abort_message='Pb allocation tmp_pctsrf_slab'
        CALL abort_physic(modname,abort_message,1)
     ENDIF
-    tmp_pctsrf_slab(:) = pctsrf_rst(:)
-    
-! Allocate some other variables internal in module mod_oceanslab
-    ALLOCATE(tmp_radsol(ngrid), stat = error)
+    fsic(:)=0.
+    !zmasq = continent fraction
+    WHERE (1.-zmasq(:)>EPSFRA)
+!        fsic(:) = MIN(pctsrf_rst(:,is_sic)/(1.-zmasq(:)),ice_frac_max)
+        fsic(:) = MIN(pctsrf_rst(:)/(1.-zmasq(:)),ice_frac_max)
+    END WHERE
+
+!************************************************************************************
+! Allocate saved fields
+!************************************************************************************
+    ALLOCATE(tslab(klon,nslay), stat=error)
+       IF (error /= 0) CALL abort_physic &
+         (modname,'pb allocation tslab', 1)
+       tslab(:,:) = tslab_rst(:,:)
+
+    ALLOCATE(bils_cum(klon), stat = error)
     IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_radsol'
+       abort_message='Pb allocation slab_bils_cum'
        CALL abort_physic(modname,abort_message,1)
     ENDIF
-
-    ALLOCATE(tmp_flux_o(ngrid), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_o'
-       CALL abort_physic(modname,abort_message,1)
+    bils_cum(:) = 0.0
+
+!    IF (version_ocean=='sicINT') THEN ! interactive sea ice
+        ALLOCATE(slab_bilg(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_bilg'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        slab_bilg(:) = 0.0
+        ALLOCATE(bilg_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_bilg_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        bilg_cum(:) = 0.0
+        ALLOCATE(tice(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_tice'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        tice(:) = tice_rst(:)
+        ALLOCATE(seaice(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation slab_seaice'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        seaice(:) = seaice_rst(:)
+!    END IF
+
+    IF (slab_hdiff) THEN !horizontal diffusion
+        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_hdiff'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_hdiff(:,:) = 0.0
     ENDIF
-    
-    ALLOCATE(tmp_flux_g(ngrid), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation tmp_flux_g'
-       CALL abort_physic(modname,abort_message,1)
+
+    IF (slab_gm) THEN !GM advection
+        ALLOCATE(dt_gm(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_gm'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_gm(:,:) = 0.0
     ENDIF
 
-! a mettre un slab_bils aussi en force !!!
-    ALLOCATE(slab_bils(ngrid), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation slab_bils'
-       CALL abort_physic(modname,abort_message,1)
+!    IF (slab_ekman.GT.0) THEN ! ekman transport
+    IF (slab_ekman) THEN ! ekman transport
+        ALLOCATE(dt_ekman(klon,nslay), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation dt_ekman'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        dt_ekman(:,:) = 0.0
+        ALLOCATE(taux_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation taux_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        taux_cum(:) = 0.0
+        ALLOCATE(tauy_cum(klon), stat = error)
+        IF (error /= 0) THEN
+           abort_message='Pb allocation tauy_cum'
+           CALL abort_physic(modname,abort_message,1)
+        ENDIF
+        tauy_cum(:) = 0.0
     ENDIF
-    slab_bils(:) = 0.0   
-
-    ALLOCATE(dt_hdiff(ngrid,noceanmx), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation dt_hdiff'
-       CALL abort_physic(modname,abort_message,1)
+
+! Initialize transport
+    IF (slab_hdiff.OR.(slab_ekman)) THEN
+      CALL gather(zmasq,zmasq_glo)
+! Master thread/process only
+!$OMP MASTER
+      IF (is_mpi_root) THEN
+          CALL ini_slab_transp(zmasq_glo)
+      END IF
+!$OMP END MASTER
+    END IF
+
+ END SUBROUTINE ocean_slab_init
+!
+!***********************************************************************************
+!
+  SUBROUTINE ocean_slab_frac(pctsrf_chg, zmasq)
+
+! This routine sends back the sea ice and ocean fraction to the main physics routine.
+! Called only with interactive sea ice.
+
+! Arguments
+!************************************************************************************
+    REAL, DIMENSION(klon), INTENT(IN)                        :: zmasq   ! proxy of rnat
+    REAL, DIMENSION(klon), INTENT(OUT) :: pctsrf_chg  ! sea-ice fraction
+
+       pctsrf_chg(:)=fsic(:)*(1.-zmasq(:))
+
+  END SUBROUTINE ocean_slab_frac
+!
+!************************************************************************************
+!
+  SUBROUTINE ocean_slab_noice(itime, dtime, knon, knindex, &
+       precip_snow, tsurf_in, &
+       radsol, snow, fluxsens, &
+       tsurf_new, flux_u1, flux_v1, zmasq)
+
+    USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman2,slab_gmdiff
+    USE mod_phys_lmdz_para
+
+! This routine
+! (1) computes surface turbulent fluxes over points with some open ocean
+! (2) reads additional Q-flux (everywhere)
+! (3) computes horizontal transport (diffusion & Ekman)
+! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
+
+! Note :
+! klon total number of points
+! knon number of points with open ocean (varies with time)
+! knindex gives position of the knon points within klon.
+! In general, local saved variables have klon values
+! variables exchanged with PBL module have knon.
+
+! Input arguments
+!***********************************************************************************
+    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
+    INTEGER, INTENT(IN)                  :: knon  ! number of points
+    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
+    REAL, INTENT(IN) :: dtime  ! timestep (s)
+    REAL, DIMENSION(klon), INTENT(IN)    :: precip_snow !, precip_rain
+
+    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
+    REAL, DIMENSION(klon), INTENT(IN) :: radsol ! net surface (radiative) flux
+    REAL, DIMENSION(klon), INTENT(IN)   :: flux_u1, flux_v1 ! Comes from turbdiff
+    REAL, DIMENSION(klon), INTENT(IN)   :: fluxsens !, sensible heat flux
+    REAL, DIMENSION(klon), INTENT(IN)   :: zmasq   ! proxy of rnat
+
+! In/Output arguments
+!************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
+
+! Output arguments
+!************************************************************************************
+    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
+
+! Local variables
+!************************************************************************************
+    INTEGER               :: i,ki,k
+    REAL                  :: t_cadj
+
+    ! for new ice creation
+    REAL                  :: e_freeze
+    REAL, DIMENSION(klon) :: slab_bils ! weighted surface heat flux
+    ! horizontal diffusion and Ekman local vars
+    ! dimension = global domain (klon_glo) instead of // subdomains
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo
+    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
+    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
+    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
+    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
+
+!****************************************************************************************
+! 1) Surface fluxes calculation
+!    Points with some open ocean only
+!****************************************************************************************
+! save total cumulated heat fluxes locally
+! radiative + turbulent + melt of falling snow
+    slab_bils(:)=0.
+    DO i=1,knon
+        ki=knindex(i)
+        slab_bils(ki)=(1.-fsic(ki))*(fluxsens(ki)+radsol(ki) &
+                      -precip_snow(ki)*ice_lat*(1.+snow_wfact*fsic(ki)))
+        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
+    END DO
+
+!  Compute surface wind stress
+!    CALL calcul_flux_wind(knon, dtime, &
+!         u0, v0, u1, v1, gustiness, cdragm, &
+!         AcoefU, AcoefV, BcoefU, BcoefV, &
+!         p1lay, temp_air, &
+!         flux_u1, flux_v1)
+
+! save cumulated wind stress
+    IF (slab_ekman) THEN
+      DO i=1,knon
+          ki=knindex(i)
+          taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*(1.-fsic(ki))/cpl_pas
+          tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*(1.-fsic(ki))/cpl_pas
+      END DO
     ENDIF
-    dt_hdiff = 0.0   
-
-    ALLOCATE(dt_ekman(ngrid,noceanmx), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation dt_hdiff'
-       CALL abort_physic(modname,abort_message,1)
+
+!****************************************************************************************
+! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file
+! limit_slab.nc
+!
+!****************************************************************************************
+!    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
+    ! lmt_bils and diff_sst,siv saved by limit_slab
+    ! qflux = total QFlux correction (in W/m2)
+!    IF (qflux_bnd.EQ.2) THEN
+!        dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400.
+!        dt_qflux_sic(:) = -diff_siv(:)*ice_den*ice_lat/86400.
+!    ELSE
+!        dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400.  &
+!                  - diff_siv(:)*ice_den*ice_lat/86400.
+!    END IF
+!    IF (nslay.GT.1) THEN
+!       dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay)
+!    END IF
+
+!****************************************************************************************
+! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
+!    Bring to freezing temp and make sea ice if necessary
+!
+!***********************************************o*****************************************
+    tsurf_new=tsurf_in
+    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
+! ***********************************
+!  Horizontal transport
+! ***********************************
+      IF (slab_ekman) THEN
+          ! copy wind stress to global var
+          CALL gather(taux_cum,taux_glo)
+          CALL gather(tauy_cum,tauy_glo)
+      END IF
+
+      IF (slab_hdiff.OR.(slab_ekman)) THEN
+        CALL gather(tslab,tslab_glo)
+      ! Compute horiz transport on one process only
+        IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus
+          IF (slab_hdiff) THEN
+              dt_hdiff_glo(:,:)=0.
+          END IF
+          IF (slab_ekman) THEN
+              dt_ekman_glo(:,:)=0.
+          END IF
+          IF (slab_gm) THEN
+              dt_gm_glo(:,:)=0.
+          END IF
+          DO i=1,cpl_pas ! time splitting for numerical stability
+!            IF (slab_ekman.GT.0) THEN
+!              SELECT CASE (slab_ekman)
+!                CASE (1) ! 1.5 layer scheme
+!                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
+!                CASE (2) ! 2-layers
+!                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
+!                CASE DEFAULT
+!                  dt_ekman_tmp(:,:)=0.
+!              END SELECT
+            IF (slab_ekman) THEN
+              CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
+
+              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
+              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1
+              DO k=1,nslay
+                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
+              ENDDO
+              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
+              IF (slab_gm) THEN ! Gent-McWilliams eddy advection
+                dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
+                ! convert dt from K.s-1.(kg.m-2) to K.s-1
+                DO k=1,nslay
+                  dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den)
+                END DO
+                tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
+              END IF
+            ENDIF
+            IF (slab_hdiff) THEN ! horizontal diffusion
+              ! laplacian of slab T
+              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
+              ! multiply by diff coef and normalize to 50m slab equivalent
+              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
+              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
+              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
+            END IF
+          END DO ! time splitting
+          IF (slab_hdiff) THEN
+            !dt_hdiff_glo saved in W/m2
+            DO k=1,nslay
+              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
+            END DO
+          END IF
+          IF (slab_gm) THEN
+            !dt_hdiff_glo saved in W/m2
+            dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas
+          END IF
+          IF (slab_ekman) THEN
+            ! dt_ekman_glo saved in W/m2
+            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
+          END IF
+        END IF ! master process
+!$OMP BARRIER
+        ! Send new fields back to all processes
+        CALL Scatter(tslab_glo,tslab)
+        IF (slab_hdiff) THEN
+          CALL Scatter(dt_hdiff_glo,dt_hdiff)
+        END IF
+        IF (slab_gm) THEN
+          CALL Scatter(dt_gm_glo,dt_gm)
+        END IF
+        IF (slab_ekman) THEN
+          CALL Scatter(dt_ekman_glo,dt_ekman)
+          ! clear wind stress
+          taux_cum(:)=0.
+          tauy_cum(:)=0.
+        END IF
+      ENDIF ! transport
+
+! ***********************************
+! Other heat fluxes
+! ***********************************
+      ! Add cumulated ocean surface fluxes
+      tslab(:,1) = tslab(:,1) + bils_cum(:)*cyang*dtime
+      ! Convective adjustment if 2 layers
+      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
+        DO i=1,klon
+          IF (tslab(i,2).GT.tslab(i,1)) THEN
+            ! mean (mass-weighted) temperature
+            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
+            tslab(i,1)=t_cadj
+            tslab(i,2)=t_cadj
+          END IF
+        END DO
+      END IF
+      ! Add read QFlux
+!      IF (qflux_bnd.EQ.1) THEN
+!          ! QFlux from ocean circ. cannot cool tslab below freezing.
+!          dq_freeze = (t_freeze - tslab(:,1)) / (cyang*dtime*cpl_pas)
+!          dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.))
+!      ELSE IF (qflux_bnd.EQ.2) THEN
+!          dq_freeze = (t_freeze - tslab(:,1)) / (cyang*dtime*cpl_pas)  &
+!                     - dt_qflux_sic(:)
+!          dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.))  &
+!                       + dt_qflux_sic(:)
+!      END IF
+!      DO k=1,nslay
+!          tslab(:,k) = tslab(:,k) + dt_qflux(:,k)*cyang*dtime*cpl_pas &
+!                     * slabh(1)/slabh(k)
+!      END DO
+
+! ***********************************
+! Update surface temperature and ice
+! ***********************************
+!      SELECT CASE(version_ocean)
+!      CASE('sicNO') ! no sea ice even below freezing !
+!          DO i=1,knon
+!              ki=knindex(i)
+!              tsurf_new(i)=tslab(ki,1)
+!          END DO
+!      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
+!        ! tslab cannot be below freezing
+!          DO i=1,knon
+!              ki=knindex(i)
+!              IF (tslab(ki,1).LT.t_freeze) THEN
+!                  tslab(ki,1)=t_freeze
+!              END IF
+!              tsurf_new(i)=tslab(ki,1)
+!          END DO
+!      CASE('sicINT') ! interactive sea ice
+          DO i=1,knon
+              ki=knindex(i)
+              ! Check if new slab temperature is below freezing
+              IF (tslab(ki,1).LT.t_freeze) THEN
+                  ! We need to form ice now over ice-free points
+                  ! Else points not seen by slab_ice
+                  IF (fsic(ki)*(1.-zmasq(ki)).LT.epsfra) THEN
+                      ! No ice present yet.
+                      ! quantity of new ice formed
+                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat  &
+                               +fsic(ki)*seaice(ki)
+                      ! new ice forms at h_ice_thin
+                      tsurf_new(ki)=t_freeze
+                      tice(ki)=t_freeze
+                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
+                      IF (fsic(ki).GT.ice_frac_min) THEN
+                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
+                          tslab(ki,1)=t_freeze
+                      ELSE
+                          fsic(ki)=0.
+                      END IF
+                  END IF ! sea ice present
+                  ! if ice already present, new ice formed in slab_ice routine.
+!              ELSE ! temperature above freezing
+!                  tsurf_new(i)=tslab(ki,1)
+              END IF
+          END DO
+!      END SELECT
+      bils_cum(:)=0.0! clear cumulated fluxes
+    END IF ! coupling time
+  END SUBROUTINE ocean_slab_noice
+!
+!*****************************************************************************
+
+!  SUBROUTINE ocean_slab_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, &
+!       radsol, snow, qsurf, qsol, agesno, &
+!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
+!       tsurf_new, dflux_s, dflux_l, swnet)
+
+  SUBROUTINE ocean_slab_ice(itime, dtime, knon, knindex, &
+       precip_snow, tsurf_in, &
+       radsol, snow, fluxsens, &
+       tsurf_new, evap, flux_u1, flux_v1, zmasq)
+
+!   USE calcul_fluxs_mod
+
+!   INCLUDE "YOMCST.h"
+!   INCLUDE "clesphys.h"
+
+! This routine
+! (1) computes surface turbulent fluxes over points with some sea ice
+! (2) computes condutive fluxes in the snow and ice, and ice-ocean flux
+! (3) computes snow/ice albedo
+! (4) updates snow/ice temperature, surface melt if needed.
+! (5) lateral growth if tslab < freezing
+! (6) bottom & side melt / growth depending on bottom fluxes
+! (7) updates slab temperature every cpl_pas
+
+! Note :
+! klon total number of points
+! knon number of points with open ocean (varies with time)
+! knindex gives position of the knon points within klon.
+! In general, local saved variables have klon values
+! variables exchanged with PBL module have knon.
+
+! Input arguments
+!****************************************************************************************
+    INTEGER, INTENT(IN)                  :: itime, knon !, jour
+    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_snow !, precip_rain
+    REAL, DIMENSION(klon), INTENT(IN)   :: evap, fluxsens!,fluxlat
+    REAL, DIMENSION(klon), INTENT(IN)   :: flux_u1, flux_v1
+    REAL, DIMENSION(klon), INTENT(IN)   :: zmasq   ! proxy of rnat
+!    REAL, DIMENSION(klon), INTENT(IN)    :: spechum, temp_air
+!    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)    :: swnet
+
+! In/Output arguments
+!****************************************************************************************
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow!, qsol
+!    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
+    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
+
+! 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
+
+! Local variables
+!****************************************************************************************
+    INTEGER               :: i,ki
+!    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
+!    REAL, DIMENSION(klon) :: u0, v0
+!    REAL, DIMENSION(klon) :: u1_lay, v1_lay
+    REAL, DIMENSION(klon) :: f_bot ! bottom ocean - ice flux
+
+    ! intermediate heat fluxes:
+    ! (conduction snow / ice, shortwave penetration, bottom turbulent fluxes)
+    REAL                  :: f_cond_s, f_cond_i, f_turb
+    ! friction velocity, 1/C and 1/tau conduction for ice
+    REAL                  :: ustar
+!    REAL                  :: 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, e_freeze, snow_evap, h_test, h_new
+    ! dhsic, dfsic change in ice mass, fraction.
+    ! frac_mf ratio of lateral / thickness growth / melt
+    REAL                  :: dhsic, dfsic, frac_mf
+
+!*******************************************************************************
+! 1) Update surface , ice and slab temperature
+!*******************************************************************************
+! Wind stress
+! flux_u1, flux_v1 from physics
+! save cumulated wind stress
+! Use ocean-ice stress = wind stress * (1.-fsic)
+!    IF (slab_ekman.GT.0) THEN
+    IF (slab_ekman) THEN
+      DO i=1,knon
+          ki=knindex(i)
+          taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*fsic(ki)*(1.-fsic(ki))/cpl_pas
+          tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*fsic(ki)*(1.-fsic(ki))/cpl_pas
+      END DO
     ENDIF
-    dt_ekman = 0.0   
-
-
-    ALLOCATE(lmt_bils(ngrid), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation lmt_bils'
-       CALL abort_physic(modname,abort_message,1)
+
+! Initialize ice-ocean flux
+    slab_bilg(:)=0.
+
+    ! Old, explicit scheme for snow & ice conductive fluxes
+    ! radsol is total surface fluxes (input) : radiative + turbulent
+        DO i=1,knon
+        ki=knindex(i) ! For PCM : you can probably replace ki by i
+            ! ocean-ice turbulent heat flux
+            ! turbulent velocity = (tau/rho)^1/2
+            ustar = MAX(5e-4, SQRT((1.-fsic(ki))   &
+                   * SQRT(flux_u1(ki)**2 + flux_v1(ki)**2) / sea_den))
+            f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze)
+            ! f_turb >0 and cannot bring tslab below zero
+            f_turb = MAX(0., MIN(f_turb, &
+                        (tslab(ki,1)-t_freeze) / (cyang*dtime*cpl_pas)))
+
+            ! ice conductive flux (pos up)
+            IF (seaice(ki).GT.0) THEN
+                f_cond_i = ice_cond*(t_freeze-tice(ki))/seaice(ki)
+            ELSE
+                f_cond_i = 0
+            END IF
+
+            ! If snow layer present, tsurf = tsnow
+            ! Problem here, as tsurf_in # tsnow ?
+            IF (snow(ki).GT.snow_min) THEN
+                ! snow conductive flux (pos up)
+                f_cond_s=sno_cond*(tice(ki)-tsurf_in(ki))/snow(ki)
+                ! update ice temperature
+                tice(ki)=tice(ki) + 2./ice_cap/(snow(ki)+seaice(ki)) &
+                *(f_cond_i-f_cond_s)*dtime
+                ! update snow temperature
+                tsurf_new(ki) = tsurf_in(ki) + 2./ice_cap/snow(ki) &
+                    *(fluxsens(ki)+radsol(ki)+f_cond_s)*dtime
+            ELSE IF (seaice(ki).GT.0) THEN ! bare ice. tsurf = tice
+                ! update ice temperature
+                        tice(ki) = tice(ki) + 2./ice_cap/seaice(ki) &
+                                *(fluxsens(ki)+radsol(ki)+f_cond_i)*dtime
+                        tsurf_new(ki) = tice(ki)
+            END IF
+            ! bottom flux (used to grow ice from below)
+            f_bot(ki) = f_turb - f_cond_i
+            slab_bilg(ki) = -f_turb
+        END DO
+!
+!! Surface turbulent fluxes (sens, lat etc) and update surface temp.
+!    dif_grnd(:)=0.
+!    beta(:) = 1.
+!    CALL calcul_fluxs(knon, is_sic, dtime, &
+!        tsurf_in, 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)
+!    DO i=1,knon
+!        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
+!    END DO
+
+! Surface, snow-ice and ice-ocean fluxes.
+! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd)
+
+! Surface turbulent fluxes (sens, lat etc) and update surface temp.
+!    beta(:) = 1.
+!    CALL calcul_fluxs(knon, is_sic, dtime, &
+!        tsurf_in, 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)
+
+!! Update remaining temperature and fluxes
+!    DO i=1,knon
+!    ki=knindex(i)
+!        ! ocean-ice turbulent heat flux
+!        ! turbulent velocity = (tau/rho)^1/2 for low ice cover
+!        ! min = 5e-4 for 1cm/s current
+!        ustar = MAX(5e-4, SQRT((1.-fsic(ki))   &
+!               * SQRT(flux_u1(i)**2 + flux_v1(i)**2) / sea_den))
+!        f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze)
+!        ! ice_ocean fluxes cannot bring tslab below freezing
+!        f_turb = MAX(0., MIN(f_turb, slab_bilg(ki) + (tslab(ki,1)-t_freeze) &
+!                          / (fsic(ki)*cyang*dtime*cpl_pas)))
+!        IF (snow(i).GT.snow_min) THEN
+!            ! snow conductive flux after calcul_fluxs
+!            f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i)
+!            ! 1 / heat capacity and conductive timescale
+!            uscap = 2. / ice_cap / (snow(i)+seaice(ki))
+!            ustau = uscap * ice_cond / seaice(ki)
+!            ! 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)) / seaice(ki)
+!        f_bot(i) = f_turb - f_cond_i
+!        slab_bilg(ki) = slab_bilg(ki)-f_turb
+!    END DO
+
+    ! weight fluxes to ocean by sea ice fraction
+    slab_bilg(:)=slab_bilg(:)*fsic(:)
+
+!****************************************************************************************
+! 2) Update snow and ice surface : thickness and fraction
+!****************************************************************************************
+    DO i=1,knon
+        ki=knindex(i)
+! snow precip (could be before fluxes above ?)
+        IF (precip_snow(ki) > 0.) THEN
+            snow(ki) = snow(ki)+precip_snow(ki)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
+        END IF
+! snow and ice sublimation
+        IF (evap(ki) > 0.) THEN
+           snow_evap = MIN (snow(ki) / dtime, evap(ki))
+           snow(ki) = snow(ki) - snow_evap * dtime
+           snow(ki) = MAX(0.0, snow(ki))
+           seaice(ki) = MAX(0.0,seaice(ki)-(evap(ki)-snow_evap)*dtime)
+        ENDIF
+! Melt / Freeze snow from above if Tsurf>0
+        IF (tsurf_new(ki).GT.t_melt) THEN
+            ! energy available for melting snow (in kg of melted snow /m2)
+            e_melt = MIN(MAX(snow(ki)*(tsurf_new(ki)-t_melt)*ice_cap/2. &
+               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(ki))
+            ! remove snow
+            IF (snow(ki).GT.e_melt) THEN
+                snow(ki)=snow(ki)-e_melt
+                tsurf_new(ki)=t_melt
+            ELSE ! all snow is melted
+                ! add remaining heat flux to ice
+                e_melt=e_melt-snow(ki)
+                snow(ki)=0.0
+                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
+                tsurf_new(ki)=tice(ki)
+            END IF
+        END IF
+! Bottom melt / grow
+! bottom freeze if bottom flux (cond + oce-ice) <0
+        IF (f_bot(ki).LT.0) THEN
+           ! larger fraction of bottom growth
+           frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thick)   &
+                  / (h_ice_max-h_ice_thick)))
+           ! quantity of new ice (formed at mean ice temp)
+           e_melt= -f_bot(ki) * 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-seaice(ki),e_melt/fsic(ki)))
+           seaice(ki)=dhsic+seaice(ki)
+           e_melt=e_melt-fsic(ki)*dhsic
+           IF (e_melt.GT.0.) THEN
+           ! frac_mf fraction used for lateral increase
+           dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
+           fsic(ki)=fsic(ki)+dfsic
+           e_melt=e_melt-dfsic*seaice(ki)
+           ! rest used to increase height
+           seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
+           END IF
+       ELSE
+! melt from below if bottom flux >0
+           ! larger fraction of lateral melt from warm ocean
+           frac_mf=MIN(1.,MAX(0.,(seaice(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(ki) * dtime * fsic(ki) &
+                   / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
+           ! first decrease height to h_thick
+           IF (fsic(ki).GT.0) THEN
+             dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
+           ELSE
+             dhsic=0
+           ENDIF
+           seaice(ki)=seaice(ki)-dhsic
+           e_melt=e_melt-fsic(ki)*dhsic
+           IF (e_melt.GT.0) THEN
+           ! frac_mf fraction used for height decrease
+           dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
+           seaice(ki)=seaice(ki)-dhsic
+           e_melt=e_melt-fsic(ki)*dhsic
+           ! rest used to decrease fraction (up to 0!)
+           dfsic=MIN(fsic(ki),e_melt/seaice(ki))
+           ! keep remaining in ocean if everything melted
+           e_melt=e_melt-dfsic*seaice(ki)
+           slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
+           ELSE
+           dfsic=0
+           END IF
+           fsic(ki)=fsic(ki)-dfsic
+       END IF
+! melt ice from above if Tice>0
+        IF (tice(ki).GT.t_melt) THEN
+            ! quantity of ice melted (kg/m2)
+            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
+             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
+            ! melt from above, height only
+            dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
+            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*fsic(ki))
+              ! if all melted add remaining heat to ocean
+              e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
+              slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
+              ! update height and fraction
+              fsic(ki)=fsic(ki)-dfsic
+            END IF
+            seaice(ki)=seaice(ki)-dhsic
+            ! surface temperature at melting point
+            tice(ki)=t_melt
+            tsurf_new(ki)=t_melt
+        END IF
+! convert snow to ice if below floating line
+        h_test=(seaice(ki)+snow(ki))*ice_den-seaice(ki)*sea_den
+        IF ((h_test.GT.0.).AND.(seaice(ki).GT.h_ice_min)) THEN !snow under water
+            ! extra snow converted to ice (with added frozen sea water)
+            dhsic=h_test/(sea_den-ice_den+sno_den)
+            seaice(ki)=seaice(ki)+dhsic
+            snow(ki)=snow(ki)-dhsic*sno_den/ice_den
+            ! available energy (freeze sea water + bring to tice)
+            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
+                   ice_cap/2.*(t_freeze-tice(ki)))
+            ! update ice temperature
+            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(ki)+seaice(ki))
+        END IF
+    END DO
+
+!*******************************************************************************
+! 3) cumulate ice-ocean fluxes, update tslab, lateral grow
+!***********************************************o*******************************
+    !cumul fluxes.
+    bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
+    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab
+        ! Add cumulated surface fluxes
+        tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
+        bilg_cum(:)=0.
+! If slab temperature below freezing, new lateral growth
+        DO i=1,knon
+            ki=knindex(i)
+            IF (tslab(ki,1).LT.t_freeze) THEN ! create more ice
+                ! quantity of new ice formed over open ocean
+                ! (formed at mean ice temperature)
+                e_freeze=(t_freeze-tslab(ki,1))/cyang &
+                         /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
+                ! new ice height and fraction
+                h_new=MAX(h_ice_thin,MIN(h_ice_new,seaice(ki))) ! new height
+!                h_new=MIN(h_ice_new,seaice(ki))
+                dfsic=MIN(ice_frac_max-fsic(ki)      &
+                          ,MAX(ice_frac_min,e_freeze/h_new))
+                ! update average sea ice height
+                seaice(ki)=(seaice(ki)*fsic(ki)+e_freeze) &
+                           / (dfsic+fsic(ki))
+                ! update snow thickness
+                snow(ki) = snow(ki) * fsic(ki) / (dfsic+fsic(ki))
+                ! update tslab to freezing
+                tslab(ki,1)=t_freeze
+                ! update sea ice fraction
+                fsic(ki)=fsic(ki)+dfsic
+            END IF ! tslab below freezing
+        END DO
+    END IF ! coupling time
+
+!****************************************************************************************
+! 4) Compute sea-ice and snow albedo
+!****************************************************************************************
+! Removed all albedo stuff as it is computed through hydrol in the Generic model
+! ------ End Albedo ----------
+
+    !tests remaining ice fraction
+    WHERE ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min))
+        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
+        tice=t_melt
+        fsic=0.
+        seaice=0.
+    END WHERE
+
+  END SUBROUTINE ocean_slab_ice
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_final
+
+!****************************************************************************************
+! Deallocate module variables
+!****************************************************************************************
+    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
+    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
+    IF (ALLOCATED(tice)) DEALLOCATE(tice)
+    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
+    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
+    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
+    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
+    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
+    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
+    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
+    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
+    IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm)
+!    IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux)
+!    IF (ALLOCATED(dt_qflux_sic)) DEALLOCATE(dt_qflux_sic)
+
+  END SUBROUTINE ocean_slab_final
+!
+!****************************************************************************************
+!
+  SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_g_loc, &
+       dt_hdiff_loc,dt_ekman_loc)
+
+! "Get some variables from module ocean_slab_mod"
+
+    INTEGER, INTENT(IN)                     :: ngrid
+    REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: tslab_loc
+    REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
+    REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
+    REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_hdiff_loc ! [in W/m2]
+    REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_ekman_loc ! [in W/m2]
+    INTEGER :: i
+
+
+! Set the output variables
+    tslab_loc(:,:) = 0.
+    dt_hdiff_loc(:,:)=0.
+    dt_ekman_loc(:,:)=0.
+    tslab_loc(:,:) = tslab(:,:)
+    seaice_loc(:) = seaice(:)
+    flux_g_loc(:) = slab_bilg(:)
+
+!!      dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2
+!!      dt_ekman_loc(:,i) = dt_ekman(:,i)*slabh(i)*1000.*4228.
+
+    IF (slab_hdiff) THEN
+        DO i=1,nslay
+           dt_hdiff_loc(:,i) = dt_hdiff(:,i)
+        ENDDO
     ENDIF
-    lmt_bils(:) = 0.0
-
-    ALLOCATE(slabh(noceanmx), stat = error)
-    IF (error /= 0) THEN
-       abort_message='Pb allocation slabh'
-       CALL abort_physic(modname,abort_message,1)
+    IF (slab_ekman) THEN
+        DO i=1,nslay
+           dt_ekman_loc(:,i) = dt_ekman(:,i)
+        ENDDO
     ENDIF
-    slabh(1)=50.
-    IF (noceanmx.GE.2) slabh(2)=150.
-    IF (noceanmx.GE.3) slabh(3)=2800.
-
-
-!       CALL init_masquv
-
-
-
-
-  END SUBROUTINE ocean_slab_init
-!
-!****************************************************************************************
-!
- 
-!
-!****************************************************************************************
-!
-  SUBROUTINE ocean_slab_ice(dtime, &
-       ngrid, knindex, tsurf, radsol, &
-       precip_snow, snow, evap, &
-       fluxsens,tsurf_new, pctsrf_sic, &
-       taux_slab,tauy_slab,icount)
-
-      use slab_ice_h
-      use callkeys_mod, only: ok_slab_sic
-
-
-! Input arguments  
-!****************************************************************************************
-    INTEGER, INTENT(IN)                     :: ngrid
-    INTEGER, DIMENSION(ngrid), INTENT(IN) :: knindex
-    REAL, INTENT(IN)                        :: dtime
-    REAL, DIMENSION(ngrid), INTENT(IN)    :: tsurf
-    REAL, DIMENSION(ngrid), INTENT(IN)    :: taux_slab
-    REAL, DIMENSION(ngrid), INTENT(IN)    :: tauy_slab
-    REAL, DIMENSION(ngrid), INTENT(IN)    :: evap, fluxsens
-    REAL, DIMENSION(ngrid), INTENT(IN)    :: precip_snow
-    REAL, DIMENSION(ngrid), INTENT(IN)    :: radsol
-    INTEGER, INTENT(IN)                     :: icount
-
-
-!In/Output arguments
-!****************************************************************************************l
-    REAL, DIMENSION(ngrid), INTENT(INOUT)          :: snow
-
-!    REAL, DIMENSION(ngridmx), INTENT(INOUT)          :: agesno
-
-
-! Output arguments
-!****************************************************************************************
-!    REAL, DIMENSION(ngridmx), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
-!    REAL, DIMENSION(ngridmx), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
-    REAL, DIMENSION(ngrid), INTENT(OUT)            :: tsurf_new     
-    REAL, DIMENSION(ngrid), INTENT(OUT)            :: pctsrf_sic
-
-! Local variables
-!****************************************************************************************
-    INTEGER                                 :: i
-    REAL, DIMENSION(ngrid)                :: cal, beta, dif_grnd, capsol
-    REAL, DIMENSION(ngrid)                :: alb_neig, alb_ice!, tsurf_temp
-    REAL, DIMENSION(ngrid)                :: soilcap, soilflux
-    REAL, DIMENSION(ngrid)                :: zfra
-    REAL                                    :: snow_evap, fonte 
-    REAL, DIMENSION(ngrid,noceanmx)        :: tslab
-    REAL, DIMENSION(ngrid)                :: seaice,tsea_ice ! glace de mer (kg/m2)
-    REAL, DIMENSION(ngrid)                :: pctsrf_new
-
-
-!*****************************************************************************
-
-
-! Initialization of output variables
-!    alb1_new(:) = 0.0
-
-!******************************************************************************
-! F. Codron: 3 cas 
-! -Glace interactive, quantité seaice : T glace suit modèle simple
-! -Pas de glace: T oce peut descendre en dessous de 0°C sans geler
-!*****************************************************************************
-
-       pctsrf_new(:)=tmp_pctsrf_slab(:)
-       tmp_radsol(:)=radsol(:)
-       tmp_flux_o(:)=fluxsens(:)
-
-    DO i = 1, ngrid
-       tsurf_new(i) = tsurf(knindex(i))
-       seaice(i)=tmp_seaice(knindex(i))
-
-
-
-       IF (pctsrf_new(knindex(i)) < EPSFRA) THEN
-          snow(i) = 0.0
-          tsurf_new(i) = T_h2O_ice_liq - 1.8
-          !IF (soil_model) tsoil(i,:) = T_h2O_ice_liq -1.8
-       ENDIF
-    ENDDO
-    tmp_flux_g(:) = 0.0
-    tsea_ice(:)=T_h2O_ice_liq-1.8
-! Calculs flux glace/mer et glace/air
-    IF (ok_slab_sic) THEN   
-!*********************************
-! Calcul de beta, heat capacity
-!********************************* 
-!    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
-    
-!    IF ((soil_model)) THEN 
-
-!    ELSE 
-!       dif_grnd = 1.0 / tau_gl
-!       cal = RCPD * calice
-!       WHERE (snow > 0.0) cal = RCPD * calsno 
-!    ENDIF
-!    tsurf_temp = tsurf_new
-!    beta = 1.0
-
-
-! **********************************************
-! Evolution avec glace interactive:
-! *************************************
-!      tsurf_new=tsurf_temp
-      DO i = 1, ngrid
-        IF (pctsrf_new(knindex(i)) .GT. epsfra) THEN
-! *************************************
-! + Calcul Flux entre glace et océan 
-! *************************************
-          tmp_flux_g(knindex(i))=(tsurf_new(i)-(T_h2O_ice_liq-1.8)) &
-            * ice_cond*ice_den/seaice(i)
-
-! ****************************************
-! + Calcul nouvelle température de la glace 
-! ****************************************
-          tsurf_new(i)=tsurf_new(i)+2*(radsol(i)+fluxsens(i) &
-             -tmp_flux_g(knindex(i))) &
-             /(ice_cap*seaice(i))*dtime
-
-! ***************************************
-! + Precip and evaporation of snow and ice
-! ***************************************
-! Add precip
-          IF (precip_snow(i).GT.0.) THEN
-            snow(i)=snow(i)+precip_snow(i)*dtime
-          ENDIF
-! Evaporation 
-          snow_evap=0.
-          IF (evap(i).GT.0.) THEN
-            snow_evap = MIN (snow(i) / dtime, evap(i))
-            snow(i) = snow(i) - snow_evap*dtime
-            snow(i) = MAX(0.0, snow(i))
-            seaice(i)=MAX(0.0,seaice(i)-(evap(i)-snow_evap)*dtime)
-          ENDIF
-        
-! *****************************************
-! + Fonte neige & glace par le haut
-! *****************************************
-! Snow melt
-          IF ((snow(i) > epsfra) .AND. (tsurf_new(i) >= T_h2O_ice_liq)) THEN
-            fonte=MIN(MAX((tsurf_new(i)-T_h2O_ice_liq)*seaice(i) &
-                  /2.*ice_cap/ice_lat,0.0),snow(i))
-            snow(i) = MAX(0.,(snow(i)-fonte))
-            tsurf_new(i)=tsurf_new(i)-fonte*2*ice_lat/(seaice(i)*ice_cap)
-          ENDIF
-          snow(i)=MIN(snow(i),3000.)
-! Ice melt
-          IF ((seaice(i) > epsfra) .AND. (tsurf_new(i) >= T_h2O_ice_liq)) THEN       
-             fonte=seaice(i)*ice_cap*(tsurf_new(i)-T_h2O_ice_liq) &
-                  /(2*ice_lat+ice_cap*1.8)
-             CALL icemelt(fonte,pctsrf_new(knindex(i)),seaice(i))
-             tmp_flux_g(knindex(i))=tmp_flux_g(knindex(i)) &
-                                 +fonte*ice_lat/dtime
-             tsurf_new(i)=T_h2O_ice_liq
-          ENDIF  
-          tmp_seaice(knindex(i))=seaice(i)
-        ENDIF!(pctsrf_new(knindex(i)) .GT. epsfra)
-          tsea_ice(knindex(i))=tsurf_new(i)
-
-      ENDDO
-    ENDIF
-! ******************************************
-! CALL interfoce:
-! cumul/calcul nouvelle T_h2O_ice_liqce
-! fonte/formation de glace en dessous
-! ******************************************   
-    tslab = tmp_tslab
-    
-    CALL interfoce_slab(ngrid, dtime, &
-         tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
-         tslab, tsea_ice, pctsrf_new,taux_slab,tauy_slab,icount)
-    
-    tmp_pctsrf_slab(:)=pctsrf_new(:)
-!    tmp_pctsrf_slab(:,is_oce)=1.0-tmp_pctsrf_slab(:) &
-!            -tmp_pctsrf_slab(:,is_lic)-tmp_pctsrf_slab(:,is_ter)
-
-    DO i=1, ngrid
-       tmp_tslab(knindex(i),:)=tslab(knindex(i),:)
-       seaice(i)=tmp_seaice(knindex(i))
-       tsurf_new(i)=tsea_ice(knindex(i))
-
-    ENDDO
-
-! ****************************
-! calcul new albedo
-! ****************************
-! Albedo neige
-!    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:))  
-!    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
-! Fraction neige (hauteur critique 45kg/m2~15cm)
-!    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/45.0))
-! Albedo glace
-!    IF (ok_slab_sicOBS) THEN
-!       alb_ice=0.6
-!    ELSE
-!       alb_ice(1:knon)=alb_ice_max-(alb_ice_max-alb_ice_min) &
-!         *exp(-seaice(1:knon)/h_alb_ice)
-!    ENDIF
-!Albedo final
-!    alb1_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + & 
-!         alb_ice * (1.0-zfra(1:knon))
-!    alb2_new(:) = alb1_new(:)
-
-        
-! ********************************
-! Return the fraction of sea-ice
-! ********************************
-    pctsrf_sic(:) =  tmp_pctsrf_slab(:)
-
-
-  END SUBROUTINE ocean_slab_ice
-
-
-! 
-!***************************************************************************
-!
-  SUBROUTINE interfoce_slab(ngrid, dtime, &
-       radsol, fluxo, fluxg, pctsrf, &
-       tslab, tsea_ice, pctsrf_slab, &
-       taux_slab, tauy_slab,icount)
-!
-! Cette routine calcule la temperature d'un slab ocean, la glace de mer 
-! et les pourcentages de la maille couverte par l'ocean libre et/ou 
-! la glace de mer pour un "slab" ocean de 50m
-!
-! Conception: Laurent Li
-! Re-ecriture + adaptation LMDZ4: I. Musat
-! Transport, nouveau modèle glace, 2 couches: F.Codron
-!
-! input:
-!   klon         nombre T_h2O_ice_liqtal de points de grille
-!   itap         numero du pas de temps
-!   dtime        pas de temps de la physique (en s)
-!   ijour        jour dans l'annee en cours
-!   radsol       rayonnement net au sol (LW + SW)
-!   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques 
-!   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
-!   pctsrf       tableau des pourcentages de surface de chaque maille
-! output: 
-!   tslab        temperature de l'ocean libre
-!   tsea_ice         temperature de la glace (surface)
-!   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
-
-    use slab_ice_h 
-    use callkeys_mod, only: ok_slab_sic,ok_slab_heat_transp
-
-! Input arguments
-!****************************************************************************************
-    INTEGER, INTENT(IN)                       :: ngrid,icount
-!    INTEGER, INTENT(IN)                       :: itap
-    REAL, INTENT(IN)                          :: dtime       ! not used
-!    INTEGER, INTENT(IN)                       :: ijour
-    REAL, DIMENSION(ngrid), INTENT(IN)         :: radsol
-    REAL, DIMENSION(ngrid), INTENT(IN)         :: fluxo
-    REAL, DIMENSION(ngrid), INTENT(IN)         :: fluxg
-    REAL, DIMENSION(ngrid), INTENT(IN)  :: pctsrf
-    REAL, DIMENSION(ngrid), INTENT(IN)  :: taux_slab
-    REAL, DIMENSION(ngrid), INTENT(IN)  :: tauy_slab
-
-! Output arguments
-!****************************************************************************************
-    REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: tslab
-    REAL, DIMENSION(ngrid), INTENT(INOUT)       :: pctsrf_slab
-    REAL, DIMENSION(ngrid), INTENT(INOUT)       :: tsea_ice
-
-! Local variables
-!****************************************************************************************
-    REAL                    :: fonte,t_cadj
-    INTEGER                 :: i,k,kt,kb
-    REAL                    :: zz, za, zb
-    REAL                    :: cyang ! capacite calorifique slab J/(m2 K)
-    REAL, PARAMETER         :: tau_conv=432000. ! temps ajust conv (5 jours)
-    REAL, DIMENSION(ngrid,noceanmx) :: dtdiff_loc,dtekman_loc
-
-
-
-!***************************************************
-! Capacite calorifique de la couche de surface
-    cyang=capcalocean!slabh(1)*4.228e+06
-    cpl_pas=1!4*iradia
-
-!
-! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
-!!    julien = MOD(ijour,360)
-!!    IF (ok_slab_heaT_h2O_ice_liqBS .and. (MOD(itap,lmt_pas) .EQ. 1)) THEN  
-!!       ! 1er pas de temps de la journee
-!!       idayvrai = ijour
-!!       CALL condsurf(julien,idayvrai, lmt_bils)
-!!    ENDIF
-
-! ----------------------------------------------
-! A chaque pas de temps: cumul flux de chaleur
-! ----------------------------------------------
-! bilan du flux de chaleur dans l'ocean :
-!
-    DO i = 1, ngrid
-       za = radsol(i) + fluxo(i)
-       zb = fluxg(i)
-!       IF(((1-pctsrf(i)).GT.epsfra).OR. &
-!            (pctsrf(i).GT.epsfra)) THEN
-          slab_bils(i)=slab_bils(i)+(za*(1-pctsrf(i)) &
-               +zb*pctsrf(i))/ cpl_pas
-!       ENDIF
-
-    ENDDO !klon
-
-! ---------------------------------------------
-! T_h2O_ice_liqus les cpl_pas: update de tslab et seaice
-! ---------------------------------------------
-
-    IF (mod(icount-1,cpl_pas).eq.0) THEN !fin de journee
-!
-! ---------------------------------------------
-! Transport de chaleur par circulation
-! Decoupage par pas de temps pour stabilite numérique 
-! (diffusion, schema centre pour advection)
-! ---------------------------------------------
-      dt_hdiff(:,:)=0.
-      dt_ekman(:,:)=0.
-
-      IF (ok_slab_heat_transp) THEN
-!       DO i=1,cpl_pas
-!  Transport diffusif
-!         IF (ok_soil_hdiff) THEN
-             CALL divgrad_phy(ngrid,noceanmx,tmp_tslab_loc,dtdiff_loc)
-             dtdiff_loc=dtdiff_loc*soil_hdiff*50./SUM(slabh)!*100.
- !           dtdiff_loc(:,1)=dtdiff_loc(:,1)*soil_hdiff*50./SUM(slabh)*0.8
- !           dtdiff_loc(:,2)=dtdiff_loc(:,2)*soil_hdiff*50./SUM(slabh)*0.2
-             dt_hdiff=dt_hdiff+dtdiff_loc
-             tmp_tslab_loc=tmp_tslab_loc+dtdiff_loc*dtime
-!         END IF
-
-! Calcul  de transport par Ekman
-
-          CALL slab_ekman2(ngrid,taux_slab,tauy_slab,tslab,dtekman_loc)
-
-
-!        END SELECT
-!        IF (slab_ekman.GT.0) THEN
-          do k=1,noceanmx
-            dtekman_loc(:,k)=dtekman_loc(:,k)/(slabh(k)*1000.)!*0.
-          enddo
-          dt_ekman(:,:)=dt_ekman(:,:)+dtekman_loc(:,:)
-          tmp_tslab_loc=tmp_tslab_loc+dtekman_loc*dtime
-!        ENDIF
-!      ENDDO ! time splitting 1,cpl_pas     
-!      IF (slab_ekman.GT.0) THEN
-!         taux_slab(:)=0.
-!         tauy_slab(:)=0.
-!      ENDIF
-
-!       print*, 'slab_bils=',slab_bils(1)
-
-
-      ENDIF!(ok_slab_heat_transp)
-
-      DO i = 1, ngrid
-!      IF (((1-pctsrf(i)).GT.epsfra).OR. &
-!             (pctsrf(i).GT.epsfra)) THEN
-! ---------------------------------------------
-! Ajout des flux de chaleur dans l'océan
-! ---------------------------------------
-
-!print*, 'T_h2O_ice_liqcean1=',tmp_tslab_loc(i,1)
-        tmp_tslab_loc(i,1) = tmp_tslab_loc(i,1) + &
-             slab_bils(i)/cyang*dtime*cpl_pas
-
-!print*, 'capcalocean=',capcalocean
-!print*, 'cyang=',cyang
-!print*, 'dT_h2O_ice_liqcean=',slab_bils(i)/cyang*dtime
-!print*, 'T_h2O_ice_liqcean2=',tmp_tslab_loc(i,1)
-
-
-! on remet l'accumulation a 0
-        slab_bils(i) = 0.
-! ---------------------------------------------
-! Glace interactive 
-! ---------------------------------------------
-        IF (ok_slab_sic) THEN
-! Fondre la glace si température > 0. 
-! -----------------------------------
-          IF ((tmp_tslab_loc(i,1).GT.T_h2O_ice_liq-1.8) .AND. (tmp_seaice(i).GT.0.0)) THEN
-            fonte=(tmp_tslab_loc(i,1)-T_h2O_ice_liq+1.8)*cyang &
-                /(ice_lat+ice_cap/2*(T_h2O_ice_liq-1.8-tsea_ice(i)))
-            CALL icemelt(fonte,pctsrf_slab(i),tmp_seaice(i))
-                 tmp_tslab_loc(i,1)=T_h2O_ice_liq-1.8+fonte*ice_lat/cyang
-          ENDIF 
-! fabriquer de la glace si congelation atteinte:
-! ----------------------------------------------
-          IF (tmp_tslab_loc(i,1).LT.(T_h2O_ice_liq-1.8)) THEN
-
-            IF (tmp_seaice(i).LT.h_ice_min) THEN
-! Creation glace nouvelle
-!             IF (pctsrf_slab(i).LT.(epsfra)) THEN
-                 fonte=(T_h2O_ice_liq-1.8-tmp_tslab_loc(i,1))*cyang/ice_lat
-              IF (fonte.GT.h_ice_min*ice_frac_min) THEN
-                tmp_seaice(i)=MIN(h_ice_thin,fonte/ice_frac_min)
-                tmp_seaice(i)=MAX(tmp_seaice(i),fonte/ice_frac_max)
-!             IF (fonte.GT.0) THEN
-!               tmp_seaice(i)=fonte
-                tsea_ice(i)=T_h2O_ice_liq-1.8
-                pctsrf_slab(i)=(1-pctsrf_slab(i))*fonte/tmp_seaice(i)
-!                pctsrf_slab(i)=1.
-                tmp_tslab_loc(i,1)=T_h2O_ice_liq-1.8
-              ENDIF
-            ELSE  
-! glace déjà présente 
-! Augmenter epaisseur
-              fonte=(T_h2O_ice_liq-1.8-tmp_tslab_loc(i,1))*cyang &
-                   /(ice_lat+ice_cap/2*(T_h2O_ice_liq-1.8-tsea_ice(i)))           
-              zz=tmp_seaice(i)
-              tmp_seaice(i)=MAX(zz,MIN(h_ice_thick,fonte+zz))
-! Augmenter couverture (oce libre et h>h_thick) 
-              za=1-pctsrf_slab(i)
-              zb=pctsrf_slab(i)
-              fonte=fonte*za+MAX(0.0,fonte+zz-tmp_seaice(i))*zb
-              pctsrf_slab(i)=MIN(zb+fonte/tmp_seaice(i), &
-                                    (za+zb)*ice_frac_max)
-              fonte=MAX(0.0,fonte-(pctsrf_slab(i)-zb)*tmp_seaice(i))
-! Augmenter epaisseur si couverture complete
-              tmp_seaice(i)=tmp_seaice(i)+fonte/pctsrf_slab(i)
-              tmp_tslab_loc(i,1) = T_h2O_ice_liq-1.8
-            ENDIF ! presence glace
-          ENDIF !congelation
-! vérifier limites de hauteur de glace
-          IF(tmp_seaice(i).GT.h_ice_min) THEN
-            tmp_seaice(i) = MIN(tmp_seaice(i),h_ice_max)
-          ELSE
-            tmp_seaice(i) = 0. 
-            pctsrf_slab(i)=0.
-          ENDIF! (tmp_seaice(i).GT.h_ice_min)
-       
-        ENDIF !(ok_slab_sic) Glace interactive
-! ----------------------------------
-! Ajustement convectif si > 1 layers
-! ----------------------------------
-
-        IF ((noceanmx.GT.1)) THEN
-          DO kt=1,noceanmx-1
-            kb=kt           
-            DO k=kt+1,noceanmx !look for instability
-              IF (tmp_tslab_loc(i,k).GT.tmp_tslab_loc(i,kt)) THEN
-                kb=k
-              ENDIF
-            ENDDO
-            IF (kb.GT.kt) THEN !ajust conv
-!               IF (slab_cadj.EQ.1) THEN
-!             :   t_cadj=SUM(tmp_tslab_loc(i,kt:kb)*slabh(kt:kb))/SUM(slabh(kt:kb))
-!                DO k=kt,kb
-!                  tmp_tslab_loc(i,k)=t_cadj
-!                ENDDO
-!              ELSEIF (slab_cadj.EQ.2) THEN
-                t_cadj=(tmp_tslab_loc(i,kb)-tmp_tslab_loc(i,kt))*dtime/tau_conv*cpl_pas
-                tmp_tslab_loc(i,kt)=tmp_tslab_loc(i,kt)+t_cadj
-                tmp_tslab_loc(i,kb)=tmp_tslab_loc(i,kb)-t_cadj*slabh(kt)/slabh(kb)
-              !ENDIF 
-            ENDIF
-          ENDDO 
-        ENDIF !ajust 2 layers
-!      ENDIF !pctsrf
-      ENDDO !klon
-
-! On met a jour la temperature et la glace
-    tslab  = tmp_tslab_loc
-
-
-    ENDIF !(mod(icount-1,cpl_pas).eq.0)
-
-  END SUBROUTINE interfoce_slab
-!
-!******************************************************************************  
-  SUBROUTINE icemelt(fonte,pctsrf,seaice)
-
-      use slab_ice_h
-
-
-
-     REAL, INTENT(INOUT)  :: pctsrf
-     REAL , INTENT(INOUT)   ::fonte,seaice !kg/m2
-     REAL :: hh !auxilliary
-     REAL :: ff !auxilliary
-
-
-! ice gt h_ice_thick: decrease thickness up T_h2O_ice_liq h_ice_thick
-     IF (seaice.GT.h_ice_thick) THEN
-        hh=seaice
-!        ff=fonte
-        seaice=max(h_ice_thick,hh-fonte)    
-        fonte=max(0.0,fonte+h_ice_thick-hh)
-
-!        seaice=max(0.,hh-fonte)    
-!        fonte=max(0.0,fonte-(seaice-hh))
-!     IF (seaice.LT.epsfra) THEN
-!        pctsrf=0.
-!        seaice=0.
-!        fonte=ff-hh
-!     ENDIF
-
-     ENDIF
-! ice gt h_ice_thin: partially decrease thickness
-     IF ((seaice.GE.h_ice_thin).AND.(fonte.GT.0.0)) THEN
-       hh=seaice
-       seaice=MAX(hh-0.6*fonte,h_ice_thin)
-       fonte=MAX(0.0,fonte-hh+seaice)
-     ENDIF
-! use rest T_h2O_ice_liq decrease area
-       hh=pctsrf
-       pctsrf=MIN(hh,MAX(0.0,hh-fonte/seaice))
-       fonte=MAX(0.0,fonte-(hh-pctsrf)*seaice)
-
-    END SUBROUTINE icemelt
-
-!****************************************************************************************
-!
-  SUBROUTINE ocean_slab_final!(tslab_rst, seaice_rst)
-
-! This subroutine will send T_h2O_ice_liq phyredem the variables concerning the slab 
-! ocean that should be written T_h2O_ice_liq restart file.
-
-!****************************************************************************************
-
-!    REAL, DIMENSION(ngridmx,noceanmx), INTENT(OUT) :: tslab_rst
-!    REAL, DIMENSION(ngridmx), INTENT(OUT) :: seaice_rst
-
-!****************************************************************************************
-! Set the output variables
-!    tslab_rst(:,:)  = tmp_tslab(:,:)
-!    tslab_rst(:)  = tmp_tslab_loc(:)
-!    seaice_rst(:) = tmp_seaice(:)
-
-! Deallocation of all variables in module
-    DEALLOCATE(tmp_tslab, tmp_tslab_loc, tmp_pctsrf_slab)
-    DEALLOCATE(tmp_seaice, tmp_radsol, tmp_flux_o, tmp_flux_g)
-    DEALLOCATE(slab_bils, lmt_bils)
-    DEALLOCATE(dt_hdiff)
-
-  END SUBROUTINE ocean_slab_final
-!
-!****************************************************************************************
-!
-  SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_o_loc, flux_g_loc, &
-       dt_hdiff_loc,dt_ekman_loc)
-  
-! "Get some variables from module ocean_slab_mod"
-! This subroutine prints variables T_h2O_ice_liq a external routine
-
-    INTEGER, INTENT(IN)                     :: ngrid
-    REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: tslab_loc
-    REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
-    REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_o_loc
-    REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
-    REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: dt_hdiff_loc
-    REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: dt_ekman_loc
-    INTEGER :: i
-
-
-! Set the output variables
-    tslab_loc=0.
-    dt_hdiff_loc=0.
-    dt_ekman_loc=0.
-    tslab_loc  = tmp_tslab
-    seaice_loc(:) = tmp_seaice(:)
-    flux_o_loc(:) = tmp_flux_o(:)
-    flux_g_loc(:) = tmp_flux_g(:)
-    DO i=1,noceanmx
-      dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2
-      dt_ekman_loc(:,i) = dt_ekman(:,i)*slabh(i)*1000.*4228. 
-    ENDDO
-  
-  
+
+
 
   END SUBROUTINE ocean_slab_get_vars
@@ -764 +1127 @@
 !****************************************************************************************
 !
+
 END MODULE ocean_slab_mod

trunk/LMDZ.GENERIC/libf/phystd/phyetat0_mod.F90

@@ -13 +13 @@
 ! to use  'getin_p'
       use ioipsl_getin_p_mod, only: getin_p
-
+!!
+  use write_field_phy, only: Writefield_phy
+!!
   use tabfi_mod, only: tabfi
-  USE tracer_h, ONLY: noms
+  USE tracer_h, ONLY: noms, igcm_h2o_vap
   USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
   use iostart, only: nid_start, open_startphy, close_startphy, &
                      get_field, get_var, inquire_field, &
                      inquire_dimension, inquire_dimension_length
-  use slab_ice_h, only: noceanmx
+!  use slab_ice_h, only: noceanmx
+  USE ocean_slab_mod, ONLY: nslay
   use callkeys_mod, only: CLFvarying,surfalbedo,surfemis, callsoil
   use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
@@ -50 +53 @@
   real,intent(out) :: cloudfrac(ngrid,nlayer)
   real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
-  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,noceanmx)  
+  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,nslay)  
   real,intent(out) :: tsea_ice(ngrid),sea_ice(ngrid)
   real,intent(out) :: rnat(ngrid)
+!  real,intent(out) :: tice(ngrid)
 
 !======================================================================
@@ -314 +318 @@
   if (.not.found) then
     write(*,*) "phyetat0: Failed loading <tslab>"
-    do iq=1,noceanmx
+    do iq=1,nslay
       tslab(1:ngrid,iq)=tsurf(1:ngrid)
     enddo
   endif
 else
-  do iq=1,noceanmx
+  do iq=1,nslay
     tslab(1:ngrid,iq)=tsurf(1:ngrid)
   enddo
@@ -338 +342 @@
 write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
              minval(tsea_ice), maxval(tsea_ice)
+
+!if (startphy_file) then
+  ! Oceanic ice temperature
+!  call get_field("tice",tice,found,indextime)
+!  if (.not.found) then
+!    write(*,*) "phyetat0: Failed loading <tice>"
+!    tice(1:ngrid)=273.15-1.8
+!  endif
+!else
+!  tice(1:ngrid)=273.15-1.8
+!endif ! of if (startphy_file)
+!write(*,*) "phyetat0: Oceanic ice temperature <tice> range:", &
+!             minval(tice), maxval(tice)
 
 if (startphy_file) then
@@ -384 +401 @@
 endif ! of if (nq.ge.1)
 
+!!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
 
 if (startphy_file) then

trunk/LMDZ.GENERIC/libf/phystd/phyredem.F90

@@ -141 +141 @@
                       put_var, put_field
   use tracer_h, only: noms
-  use slab_ice_h, only: noceanmx
+!  use slab_ice_h, only: noceanmx
+  USE ocean_slab_mod, ONLY: nslay
   use callkeys_mod, only: ok_slab_ocean, calllott_nonoro
   use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
@@ -165 +166 @@
   real,intent(in) :: rnat(ngrid)
   real,intent(in) :: pctsrf_sic(ngrid)
-  real,intent(in) :: tslab(ngrid,noceanmx)
+  real,intent(in) :: tslab(ngrid,nslay)
   real,intent(in) :: tsea_ice(ngrid)
   real,intent(in) :: sea_ice(ngrid)
+!  real,intent(in) :: tice(ngrid)
 
   integer :: iq
@@ -200 +202 @@
    call put_field("pctsrf_sic","Pourcentage sea ice",pctsrf_sic)
    call put_field("tslab","Temperature slab ocean",tslab)
-   call put_field("tsea_ice","Temperature sea ice",tsea_ice)
+   call put_field("tsea_ice","Temperature of top layer (seaice or snow)",tsea_ice)
+!   call put_field("tice","Temperature of seaice if snow on top",tice)
    call put_field("sea_ice","Sea ice mass",sea_ice)
  endif!(ok_slab_ocean)

trunk/LMDZ.GENERIC/libf/phystd/phys_state_var_mod.F90

@@ -13 +13 @@
       use comsaison_h, only: mu0, fract
       use radcommon_h, only: glat
-      use slab_ice_h, only : noceanmx
+!      use slab_ice_h, only : noceanmx
+      USE ocean_slab_mod, ONLY: nslay
       USE radinc_h, only : L_NSPECTI, L_NSPECTV,naerkind
       use surfdat_h, only: phisfi, albedodat,  &
@@ -94 +95 @@
       real,dimension(:),allocatable,save ::  pctsrf_sic
       real,dimension(:,:),allocatable,save :: tslab        ! Slab_ocean temperature (K)
-      real,dimension(:),allocatable,save ::  tsea_ice      ! Sea ice temperature (K)
+      real,dimension(:),allocatable,save ::  tsea_ice      ! Temperature of the top layer (K), either snow or ice
+!      real,dimension(:),allocatable,save :: tice           ! Sea ice temperature (K) if there is snow on top of it
       real,dimension(:),allocatable,save :: sea_ice
       real,dimension(:),allocatable,save :: zmasq
@@ -191 +193 @@
         ALLOCATE(tau_col(klon))
         ALLOCATE(pctsrf_sic(klon))
-        ALLOCATE(tslab(klon,noceanmx))
+        ALLOCATE(tslab(klon,nslay))
         ALLOCATE(tsea_ice(klon))
+!        ALLOCATE(tice(klon))
         ALLOCATE(sea_ice(klon))
         ALLOCATE(zmasq(klon))
@@ -285 +288 @@
         DEALLOCATE(tslab)
         DEALLOCATE(tsea_ice)
+ !       DEALLOCATE(tice)
         DEALLOCATE(sea_ice)
         DEALLOCATE(zmasq)

trunk/LMDZ.GENERIC/libf/phystd/physiq_mod.F90

@@ -13 +13 @@
                   pdu,pdv,pdt,pdq,pdpsrf)
 
+!!
+      use write_field_phy, only: Writefield_phy  
+!!
       use ioipsl_getin_p_mod, only: getin_p 
       use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
@@ -38 +41 @@
       use wstats_mod, only: callstats, wstats, mkstats
       use phyredem, only: physdem0, physdem1
-      use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
-                            noceanmx
-      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
-                                ini_surf_heat_transp_mod, &
-                                ocean_slab_get_vars,ocean_slab_final
-      use surf_heat_transp_mod,only :init_masquv
+!!      use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
+!!                            noceanmx
+!!      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
+!!                                capcalocean, capcalseaice,capcalsno, nslay, &
+!!                                ocean_slab_final
+!!      use surf_heat_transp_mod,only :init_masquv
+      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_noice, ocean_slab_ice, &
+                                ocean_slab_frac, ocean_slab_get_vars, &
+                                capcalocean, capcalseaice,capcalsno, nslay, knon, &
+                                ocean_slab_final
       use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
       use mod_phys_lmdz_para, only : is_master
@@ -310 +317 @@
       real zdtrain_generic(ngrid,nlayer)                      ! Rain_generic routine.
       real dtmoist(ngrid,nlayer)                              ! Moistadj routine.
-      real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine.
+      real dt_ekman(ngrid,nslay), dt_hdiff(ngrid,nslay) ! Slab_ocean routine.
       real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
       real zdtchim(ngrid,nlayer)                              ! Calchim routine.
@@ -459 +466 @@
 
       real :: tsurf2(ngrid)
-      real :: flux_o(ngrid),flux_g(ngrid)
+!!      real :: flux_o(ngrid),flux_g(ngrid)
+      real :: flux_g(ngrid)
       real :: flux_sens_lat(ngrid)
       real :: qsurfint(ngrid,nq)
@@ -530 +538 @@
 !        Read 'startfi.nc' file.
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
-#ifndef MESOSCALE
+#ifndef MESOSCALE 
          call phyetat0(startphy_file,                                 &
                        ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
@@ -536 +544 @@
                        cloudfrac,totcloudfrac,hice,                   &
                        rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
-#else
+
+!!         call WriteField_phy("post_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+#else   
+
          day_ini = pday
 #endif
@@ -570 +581 @@
          albedo_co2_ice_SPECTV(:)=0.0
          call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV) 
+
+!!        call WriteField_phy("post_surfini1_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!        call WriteField_phy("post_surfini2_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
          
 !        Initialize orbital calculation. 
@@ -586 +600 @@
          if (ok_slab_ocean)then
 
-           call ocean_slab_init(ngrid,ptimestep, tslab, &
-                                sea_ice, pctsrf_sic)
-
-           call ini_surf_heat_transp_mod()
+!!           call ocean_slab_init(ngrid,ptimestep, tslab, &
+!!                                sea_ice, pctsrf_sic)
+
+!!           call ini_surf_heat_transp_mod()
            
            knindex(:) = 0
+           knon = 0
 
            do ig=1,ngrid
               zmasq(ig)=1
-              knindex(ig) = ig
-              if (nint(rnat(ig)).eq.0) then
+!              knindex(ig) = ig
+              if (nint(rnat(ig)).eq.0) then ! rnat = 0 (ocean), 1 (continent)
                  zmasq(ig)=0
+                 knon=knon+1
+                 knindex(knon)=ig
               endif
            enddo
 
-           CALL init_masquv(ngrid,zmasq)
+           call ocean_slab_init(ptimestep, pctsrf_sic, tslab, tsea_ice, sea_ice, zmasq)
+
+!!           CALL init_masquv(ngrid,zmasq)
 
          endif ! end of 'ok_slab_ocean'.
@@ -635 +654 @@
             
          endif ! end of 'callsoil'.
+
+!!         call WriteField_phy("post_callsoil_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
          
          icount=1
@@ -687 +708 @@
          qsurf_hist(:,:)=qsurf(:,:)
 
+!!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+
 !        Initialize variable for dynamical heating and zonal wind tendency diagnostic
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -724 +747 @@
          endif
 
+!!         call WriteField_phy("post_ice_update_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+
 !        Set metallicity for GCS
 !        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
@@ -734 +759 @@
             CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV)
          endif
-         
+
          call su_watercycle ! even if we don't have a water cycle, we might
                             ! need epsi for the wvp definitions in callcorrk.F
                             ! or RETV, RLvCp for the thermal plume model
+
+!!         call WriteField_phy("post_watercycle_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
 #ifndef MESOSCALE
          if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
@@ -744 +771 @@
                          albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
          endif
+
+!!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
 #endif         
          if (corrk) then 
@@ -761 +790 @@
             call suaer_corrk       ! Set up aerosol optical properties.
          endif
+
+!!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
          ! XIOS outputs
 #ifdef CPP_XIOS
@@ -768 +799 @@
                                      year_day,presnivs,pseudoalt,WNOI,WNOV)
 #endif
+         
+!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+
          write(*,*) "physiq: end of firstcall"
       endif ! end of 'firstcall'
+
+!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
+!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))  
 
       if (check_physics_inputs) then
@@ -775 +812 @@
          call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
       endif
+
+!      call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
+!      call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
 
 ! ------------------------------------------------------
@@ -938 +978 @@
             if(rings_shadow) then
                 call call_rings(ngrid, ptime, pday, diurnal)
-            endif    
+            endif   
+
+!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 
 
 
@@ -1142 +1185 @@
       endif ! of if (callrad)
 
+!!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
 
@@ -1194 +1239 @@
          endif
 
+!!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
 
          if (tracer) then 
@@ -1200 +1246 @@
          end if ! of if (tracer)
 
+!!         call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
          ! test energy conservation
@@ -1412 +1460 @@
          dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid)
 
+!!         call writediagfi(ngrid,"condense_co2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"condense_co2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+
          ! test energy conservation
          if(enertest)then
@@ -1445 +1496 @@
                   ! Moist Convective Adjustment.
                   ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
+!!!!
+!                  call writediagfi(ngrid,"pre_moistadj_pt"," "," ",3,pt)
+!                  call writediagfi(ngrid,"pre_moistadj_pq_h2o_ice"," "," ",3,pq(:,:,igcm_h2o_ice))
+!                  call writediagfi(ngrid,"pre_moistadj_pdq_h2o_ice"," "," ",3,pdq(:,:,igcm_h2o_ice))
+!                  call writediagfi(ngrid,"pre_moistadj_pq_h2o_vap"," "," ",3,pq(:,:,igcm_h2o_vap))
+!                  call writediagfi(ngrid,"pre_moistadj_pdq_h2o_vap"," "," ",3,pdq(:,:,igcm_h2o_vap))
+!!!!
                   call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
                
@@ -1536 +1594 @@
                dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
                dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid) 
-
+               
+!!               call writediagfi(ngrid,"rain_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!               call writediagfi(ngrid,"rain_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+                
                ! Test energy conservation.
                if(enertest)then
@@ -1696 +1757 @@
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
 
+!!            call writediagfi(ngrid,"rain_generic_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+
             ! Test energy conservation.
             if(enertest)then
@@ -1779 +1842 @@
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
 
+!!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+
             ! Test water conservation
             if(watertest)then
@@ -1792 +1857 @@
          end if ! end of 'sedimentation'
 
+!!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
   ! ---------------
@@ -1829 +1896 @@
          endif
 
+!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+
+!!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
+
   ! ------------------
   !   VI.5. Slab Ocean
@@ -1839 +1911 @@
             enddo
 
-            call ocean_slab_ice(ptimestep,                          &
-                                ngrid, knindex, tsea_ice, fluxrad,  &
-                                zdqssnow, qsurf(:,igcm_h2o_ice),    &
-                              - zdqsdif(:,igcm_h2o_vap),            &
-                                flux_sens_lat,tsea_ice, pctsrf_sic, &
-                                taux,tauy,icount)
-
-
-            call ocean_slab_get_vars(ngrid,tslab,      &
-                                     sea_ice, flux_o,  &
-                                     flux_g, dt_hdiff, &
-                                     dt_ekman)
-   
+!!            call ocean_slab_ice(ptimestep,                          &
+!!                                ngrid, knindex, tsea_ice, fluxrad,  &
+!!                               zdqssnow, qsurf(:,igcm_h2o_ice),    &
+!!                              - zdqsdif(:,igcm_h2o_vap),            &
+!!                                flux_sens_lat,tsea_ice, pctsrf_sic, &
+!!                               taux,tauy,icount)
+
+
+!!            call ocean_slab_get_vars(ngrid,tslab,      &
+!!                                     sea_ice, flux_o,  &
+!!                                     flux_g, dt_hdiff, &
+!!                                     dt_ekman)
+            
+            call ocean_slab_noice(icount, ptimestep, knon, knindex, &
+                      zdqssnow, tsea_ice, &
+                      fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
+                      tsea_ice, taux, tauy, zmasq)
+
+            call ocean_slab_ice(icount, ptimestep, knon, knindex, &
+                      zdqssnow, tsea_ice, &
+                      fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
+                      tsea_ice, -zdqsdif(:,igcm_h2o_vap), taux, tauy, zmasq)
+
+            call ocean_slab_frac(pctsrf_sic, zmasq)
+
+            call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
+                      dt_hdiff, dt_ekman)
+  
             do ig=1,ngrid
                if (nint(rnat(ig)).eq.1)then
@@ -1857 +1944 @@
                   tslab(ig,2) = 0.
                   tsea_ice(ig) = 0.
+!                  tice(ig) = 0.
                   sea_ice(ig) = 0.
                   pctsrf_sic(ig) = 0.
@@ -1864 +1952 @@
 
          endif ! end of 'ok_slab_ocean'
+
+!!         call writediagfi(ngrid,"hydrol_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!         call writediagfi(ngrid,"hydrol_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
   ! -----------------------------
@@ -1876 +1967 @@
                         mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
                         sea_ice)
+                
+!!            call writediagfi(ngrid,"hydrol_post0_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
             zdtsurf(1:ngrid)     = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
             dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq)
+
             
             qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
+
+!!            call writediagfi(ngrid,"hydrol_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
+!!            call writediagfi(ngrid,"hydrol_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
 
             ! Test energy conservation
@@ -1951 +2048 @@
                capcal(ig)=capcalocean
                fluxgrd(ig)=0.
-               fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
+!!               fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
+               ! Dividing by cell area to have flux in W/m2
+               fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))/cell_area(ig)
                do iq=1,nsoilmx
                   tsoil(ig,iq)=tsurf(ig)
@@ -2345 +2444 @@
             call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
             call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
-            call wstats(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
+            call wstats(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
+!            call wstats(ngrid,"tice","ice temp if snow on top","K",2,tice)
             call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
             call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
@@ -2387 +2487 @@
       if(ok_slab_ocean) then
          call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
-         call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
+         call writediagfi(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
+!         call writediagfi(ngrid,"tice","ice temp if snow on top","K",2,tice)
          call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
          call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))