Changeset 3100 for trunk/LMDZ.GENERIC


Ignore:
Timestamp:
Oct 25, 2023, 10:34:54 AM (14 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:
3 added
2 deleted
15 edited
1 moved

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/deftank/callphys.earthslab

    r3099 r3100  
    44tracer    = .true.
    55# Diurnal cycle ?  if diurnal=false, diurnally averaged solar heating
    6 diurnal   = .false.
     6diurnal   = .true.
    77# Seasonal cycle ? if season=false, Ls stays constant, to value set in "start"
    88season    = .true.
    … …  
    1717# Test energy conservation of model physics ?
    1818enertest  = .false.
    19 # Check to see if cpp, mugaz values used match gas mixture defined in gases.def (recommended) ?
    20 check_cpp_match=.true.
     19# check if cpp and mugaz from start.nc are consistent with values computed by comp_cpp_mugaz with gases.def
     20check_cpp_match = .false.
     21# Check if physics inputs and outputs are ok
     22check_physics_inputs = .true.
     23check_physics_outputs = .true.
     24
     25## Directory where external input files are
     26## ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     27datadir     = datadir
    2128
    2229## Radiative transfer options
    … …  
    2633# the rad. transfer is computed every "iradia" physical timestep
    2734iradia     = 6
     35# Use blackbody for stellar spectrum ?
     36stelbbody    = .false.
     37# Stellar blackbody temperature ?
     38stelTbb      = 2000.000
    2839# call multilayer correlated-k radiative transfer ?
    2940corrk      = .true.
    30 # Include continuum absorption in radiative transfer (note CO2 is treated separately) ?
    31 continuum  = .true.
    3241# folder in which correlated-k data is stored ?
    33 corrkdir   = N2CO2poor_H2Ovar
     42corrkdir = Earth_JL13_extend
     43# corrkdir    = Earth_110-710K
     44# corrkdir   = N2CO2poor_H2Ovar
     45# corrkdir   = megaCO2
     46# corrkdir   = null
    3447# call visible gaseous absorption in radiative transfer ?
    3548callgasvis = .true.
     49# call continuum in radiative transfer ?
     50Continuum = .true.
    3651# Include Rayleigh scattering in the visible ?
    3752rayleigh   = .true.
     53# Uniform absorption coefficient in radiative transfer?
     54graybody     = .false.
     55# Constant absorption coefficient in visible
     56#      (in m^2/kg; only if graybody=true):
     57#      tau_surf= kappa*P/g
     58kappa_VI = 1.e-4
     59# Constant absorption coefficient in IR
     60#      (in m^2/kg; only if graybody=true):
     61kappa_IR = 5.e-1
     62# Use Newtonian cooling in place of radiative transfer ?
     63newtonian    = .false.
     64# Radiative timescale for Newtonian cooling ? [only if newtonian = T]
     65tau_relax    = 30.00000
     66# Test physics timescale in 1D ?
     67testradtimes = .false.
    3868# Characteristic planetary equilibrium (black body) temperature
    3969# This is used only in the aerosol radiative transfer setup. (see aerave.F)
    4070tplanet    = 215.
    41 # Output spectral OLR in 3D?
     71# Output spectral OLR in 1D/3D?
    4272specOLR    = .false.
    4373# Output global radiative balance in file 'rad_bal.out' - slow for 1D!!
    44 meanOLR    = .true.
     74meanOLR    = .false.
    4575# Variable gas species: Radiatively active ?
    4676varactive  = .true.
    4777# Variable gas species: Fixed vertical distribution ?
     78#   (not to be used in time integration mode)
    4879varfixed   = .false.
    4980# Variable gas species: Saturation percentage value at ground ?
    50 satval     = 0.0
     81satval     = .8
     82# Use fixed vertical profile, 1 step, no iteration ?
     83kastprof     = .false.
     84# Remove lower boundary (e.g. for gas giant sims)
     85noradsurf    = .false.
    5186
    5287## Star type
    … …  
    6095#       startype = 3            GJ644
    6196#       startype = 4            HD128167
    62 #       startype = 9            TRAPPIST-1
    63 #       startype = 10           Proxima Centauri
    6497# ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    6598# Stellar flux at 1 AU. Examples:
    … …  
    67100# 1024.5 W m-2          Sol today x 0.75 = weak early Sun
    68101# 18.462 W m-2          The feeble Gl581
     102# 19.960 W m-2          Gl581 with e=0.38 orbital average
    69103Fat1AU = 1366.0
     104
    70105
    71106## Tracer and aerosol options
    72107## ~~~~~~~~~~~~~~~~~~~~~~~~~~
    73 # Gravitational sedimentation of tracers (KEEP FALSE FOR NOW) ?
    74 sedimentation = .false.
     108# Gravitational sedimentation of tracers (just H2O ice for now) ?
     109sedimentation = .true.
    75110
    76111## Other physics options
    … …  
    92127## ~~~~~~~~~~~~~~~~~~~~~~~~~~
    93128# Number of radiatively active aerosols
    94 naerkind=1
     129naerkind = 1
     130# Varying H2O cloud fraction?
     131CLFvarying    = .true.
     132# H2O cloud fraction?
     133CLFfixval     = 1.
     134# number mixing ratio of CO2 ice particles
     135Nmix_co2      = 1.e5
     136# basic dust opacity
     137dusttau       = 0.0
     138# water cloud pressure level (norm. by psurf)
     139cloudlvl      = 0.0
     140# atm mass update due to tracer evaporation/condensation?
     141mass_redistrib = .true.
    95142# Radiatively active CO2 aerosol?
    96143aeroco2       = .false.
    … …  
    101148# Fixed water aerosol distribution?
    102149aerofixh2o  = .false.
    103 # basic dust opacity
    104 dusttau       = 0.0
    105 # Varying H2O cloud fraction?
    106 CLFvarying    = .true.
    107 # H2O cloud fraction if fixed?
    108 CLFfixval     = 0.5
    109 # fixed radii for cloud particles?
    110 radfixed=.false.
    111 # number mixing ratio of CO2 ice particles
    112 Nmix_co2      = 100000.
     150# Radiatively active sulfur aersol?
     151aeroh2so4     = .false.
     152# fixed radii for h2o cloud particles?
     153radfixed=.true.
    113154# number mixing ratio of water particles (for rafixed=.false.)
    114 Nmix_h2o      = 1.e7
     155Nmix_h2o      = 4e6
    115156# number mixing ratio of water ice particles (for rafixed=.false.)
    116 Nmix_h2o_ice      = 5.e5
    117 # radius of H2O water particles (for rafixed=.true.):
    118 rad_h2o=10.e-6
     157Nmix_h2o_ice      = 2.e4
     158# radius of H2O water particles (for rafixed=.true.): (CHANGED FROM 10 TO 12 AFTER BENJAMIN)
     159rad_h2o=12.e-6
    119160# radius of H2O ice particles (for rafixed=.true.):
    120161rad_h2o_ice=35.e-6
    121 # atm mass update due to tracer evaporation/condensation?
    122 mass_redistrib = .false.
    123162
    124163## Water options
    … …  
    130169# Model water precipitation (including coagulation etc.)
    131170waterrain     = .true.
     171# Moist adjustment
     172moistadjustment = .true.
    132173# Use simple precipitation scheme?
    133174precip_scheme=4
    134175# multiplicative constant in Boucher 95 precip scheme
    135 Cboucher=1.
    136 # WATER: hydrology ?
     176Cboucher=0.6
     177# Include hydrology ?
    137178hydrology     = .true.
    138 # Spectral Dependant Albedo ?
    139 albedo_spectral_mode=.false.
    140 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    141 # If albedo_spectral_mode=.true., albedosnow becomes the 0.5 micron snow albedo.
    142 #
    143 # albedosnow = 0.95  (0.73 Sun-integrated) for fresh snow.
    144 #            = 0.50  (0.39 Sun-integrated) for dirty snow.
    145 #            = 0.78  (0.60 Sun-integrated) for 'realistic' snow.
    146 # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     179# active runoff ?
     180activerunoff  = .true.
    147181# H2O snow (and ice) albedo ?
    148 albedosnow    = 0.60
     182albedosnow    = 0.65
    149183# Maximum sea ice thickness ?
    150184maxicethick   = 10.
    … …  
    153187# Evolve surface water sources ?
    154188sourceevol    = .false.
     189# compute lightning rate ?
     190compute_lightning     = .true.
    155191
    156192## CO2 options
    157193## ~~~~~~~~~~~
    158 # Co2 ice albedo ?
    159 albedoco2ice   = 0.5
    160 # gas is non-ideal CO2 ?
    161 nonideal      = .false.
    162194# call CO2 condensation ?
    163195co2cond       = .false.
    … …  
    165197nearco2cond   = .false.
    166198
     199## Subsurface options
     200## ~~~~~~~~~~~~~~~~~~
     201# Number of subsurface layers (For Earth 14 is OK; for Mars 18)
     202nsoilmx=14
     203# Thickness of topmost soil layer (m) (For Earth 0.03 is OK; for Mars 0.0002)
     204lay1_soil=0.03
     205#Coefficient for soil layer thickness distribution (default is 2)
     206alpha_soil=2
     207
     208## Ocean options
     209## ~~~~~~~~~~~~~
     210# Model slab-ocean (Main flag for slab ocean)
     211ok_slab_ocean = .true.
     212# The following flags can only be set to true if ok_slab_ocean is true
     213# Ekman transport
     214slab_ekman = .true.
     215# Ekman zonal advection
     216slab_ekman_zonadv = .true.
     217# Horizontal diffusion (default coef_hdiff=25000., can be changed)
     218slab_hdiff = .true.
     219# Slab-ocean timestep (in physics timesteps)
     220cpl_pas = 1
     221# sea-ice
     222#ok_slab_sic   = .true.
     223# Gent-McWilliams Scheme (can only be true if slab_ekman is true)
     224slab_gm   = .true.
     225# Slab convective adjustment? 0 - no, 1 - yes
     226slab_cadj = 1
     227# H2O snow (and ice) albedo for sea ?
     228#albedosnowsea    = 0.65
     229

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

    r3014 r3100  
    11721172
    11731173      use chimiedata_h,  only: albedo_snow_chim, albedo_co2_ice_chim
    1174       use slab_ice_h,    only: h_alb_ice,
    1175      &                         alb_ice_min, alb_ice_max
    1176       use tracer_h,      only: igcm_h2o_ice, igcm_co2_ice
     1174!      use slab_ice_h,    only: h_alb_ice, alb_ice_min, alb_ice_max
     1175      use ocean_slab_mod, only: h_alb_ice
     1176      use ocean_slab_mod, only: alb_ice_min
     1177      use ocean_slab_mod, only: alb_ice_max
    11771178      use callkeys_mod,  only: ok_slab_ocean, co2cond, alb_ocean
    11781179      use phys_state_var_mod, only: albedo_bareground,

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

    r1478 r3100  
    3636
    3737      USE comsoil_h
    38       USE slab_ice_h, only: noceanmx
     38!      USE slab_ice_h, only: noceanmx
     39      USE ocean_slab_mod, ONLY: nslay
    3940!      use control_mod
    4041      USE comvert_mod, ONLY: ap,bp,aps,bps,presnivs,pseudoalt
    … …  
    8182      INTEGER idim_tim
    8283      INTEGER idim_nsoilmx ! "subsurface_layers" dimension ID #
    83       INTEGER idim_noceanmx ! "ocean_layers" dimension ID #
     84      INTEGER idim_nslay ! "ocean_layers" dimension ID #
    8485      INTEGER nid,nvarid
    8586      real sig_s(llm),s(llm)
    … …  
    164165      ierr = NF_DEF_DIM (nid, "altitude", llm, idim_llm)
    165166      ierr = NF_DEF_DIM (nid,"subsurface_layers",nsoilmx,idim_nsoilmx)
    166       ierr = NF_DEF_DIM (nid,"ocean_layers",noceanmx,idim_noceanmx)
     167      ierr = NF_DEF_DIM (nid,"ocean_layers",nslay,idim_nslay)
    167168
    168169      ierr = NF_DEF_DIM (nid,"index", length, idim_index)

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

    r1682 r3100  
    1111
    1212use control_mod, only: nday
    13 use surf_heat_transp_mod, only: ini_surf_heat_transp
     13!use surf_heat_transp_mod, only: ini_surf_heat_transp
     14USE slab_heat_transp_mod, ONLY: ini_slab_transp_geom
    1415use infotrac, only : nqtot ! number of advected tracers
    1516use planete_mod, only: ini_planete_mod
    … …  
    8384call getin_p("ok_slab_ocean",ok_slab_ocean)
    8485if (ok_slab_ocean) then
    85   call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
    86                             cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
    87                             aireu,airev)
     86!  call ini_surf_heat_transp(ip1jm,ip1jmp1,unsairez,fext,unsaire, &
     87!                            cu,cuvsurcv,cv,cvusurcu,aire,apoln,apols, &
     88!                            aireu,airev)
     89   CALL ini_slab_transp_geom(ip1jm,ip1jmp1,unsairez,fext,unsaire,&
     90                                  cu,cuvsurcv,cv,cvusurcu, &
     91                                  aire,apoln,apols, &
     92                                  aireu,airev,rlatvdyn)
    8893endif
    8994

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

    r2484 r3100  
    99      USE tracer_h, ONLY: igcm_co2_ice
    1010      USE infotrac, ONLY: tname, nqtot
    11       USE slab_ice_h, ONLY: noceanmx
     11!      USE slab_ice_h, ONLY: noceanmx
     12      USE ocean_slab_mod, ONLY: nslay
    1213      USE callkeys_mod, ONLY: ok_slab_ocean
    1314      USE comvert_mod, ONLY: ap,bp,aps,bps,preff
    … …  
    6869      REAL,INTENT(OUT) :: emis(ngrid)
    6970      REAL,INTENT(OUT) :: q2(ngrid,llm+1),qsurf(ngrid,nqtot)
    70       REAL,INTENT(OUT) :: tslab(ngrid,noceanmx)
     71      REAL,INTENT(OUT) :: tslab(ngrid,nslay)
    7172      REAL ,INTENT(OUT) ::rnat(ngrid),pctsrf_sic(ngrid)
    7273      REAL,INTENT(OUT) :: tsea_ice(ngrid),sea_ice(ngrid)
    … …  
    9899      real emisS(iip1,jjp1)
    99100      REAL q2S(iip1,jjp1,llm+1),qsurfS(iip1,jjp1,nqtot)
    100       real tslabS(iip1,jjp1,noceanmx)
     101      real tslabS(iip1,jjp1,nslay)
    101102      real pctsrf_sicS(iip1,jjp1),tsea_iceS(iip1,jjp1)
    102103      real rnatS(iip1,jjp1), sea_iceS(iip1,jjp1)
    … …  
    287288      allocate(mlayerold(nsoilold))
    288289      allocate(qsurfold(imold+1,jmold+1,nqtot))
    289       allocate(tslabold(imold+1,jmold+1,noceanmx))
     290      allocate(tslabold(imold+1,jmold+1,nslay))
    290291      allocate(rnatold(imold+1,jmold+1))
    291292      allocate(pctsrf_sicold(imold+1,jmold+1))
    … …  
    684685      if(ok_slab_ocean) then
    685686      start=(/1,1,1,memo/)
    686       count=(/imold+1,jmold+1,noceanmx,1/)
     687      count=(/imold+1,jmold+1,nslay,1/)
    687688
    688689       ierr=NF_INQ_VARID(nid,"tslab",nvarid)
    … …  
    12981299c 6.3   Slab Ocean :
    12991300c-----------------------------------------------------------------------
    1300       call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,noceanmx,
    1301      &                   rlonuold,rlatvold,rlonu,rlatv)
    1302       call gr_dyn_fi (noceanmx,iim+1,jjm+1,ngrid,tslabs,tslab)
     1301      call interp_horiz (tslabold,tslabs,imold,jmold,iim,jjm,nslay,
     1302     &                   rlonuold,rlatvold,rlonu,rlatv)
     1303      call gr_dyn_fi (nslay,iim+1,jjm+1,ngrid,tslabs,tslab)
    13031304
    13041305      call interp_horiz (rnatold,rnats,imold,jmold,iim,jjm,1,

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

    r2785 r3100  
    3030      use phyredem, only: physdem0, physdem1
    3131      use iostart, only: open_startphy
    32       use slab_ice_h, only:noceanmx
     32!      use slab_ice_h, only:noceanmx
     33      USE ocean_slab_mod, ONLY: nslay
    3334      use filtreg_mod, only: inifilr
    3435      USE mod_const_mpi, ONLY: COMM_LMDZ

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

    r2785 r3100  
    2222      USE comsoil_h
    2323     
    24       use slab_ice_h, only: noceanmx
     24!      use slab_ice_h, only: noceanmx
     25      USE ocean_slab_mod, ONLY: nslay
    2526      USE ioipsl_getincom, only: getin
    2627      USE planete_mod, only: year_day
    … …  
    8485!     added by BC for slab ocean
    8586      REAL rnat(ngridmx),pctsrf_sic(ngridmx),sea_ice(ngridmx)
    86       REAL tslab(ngridmx,noceanmx),tsea_ice(ngridmx)
     87      REAL, ALLOCATABLE :: tslab(:,:)
     88      REAL tsea_ice(ngridmx)
    8789
    8890
    … …  
    104106!     added by BC for slab ocean
    105107      REAL rnatS(ip1jmp1),pctsrf_sicS(ip1jmp1),sea_iceS(ip1jmp1)
    106       REAL tslabS(ip1jmp1,noceanmx),tsea_iceS(ip1jmp1)
     108      REAL, ALLOCATABLE :: tslabS(:,:)
     109      REAL tsea_iceS(ip1jmp1)
    107110
    108111!     For non-orographic GW
    … …  
    156159      allocate(qsurf(ngridmx,nqtot))
    157160      allocate(qsurfS(ip1jmp1,nqtot))
     161      allocate(tslab(ngridmx,nslay)) !Added by SB for slab ocean
     162      allocate(tslabS(ip1jmp1,nslay)) !Added by SB for slab ocean
    158163! other array allocations:
    159164!      call ini_comsoil_h(ngridmx) ! done via iniphysiq
    … …  
    356361      call gr_fi_dyn(1,ngridmx,iip1,jjp1,tsea_ice,tsea_iceS)
    357362      call gr_fi_dyn(1,ngridmx,iip1,jjp1,sea_ice,sea_iceS)
    358       call gr_fi_dyn(noceanmx,ngridmx,iip1,jjp1,tslab,tslabS)
     363      call gr_fi_dyn(nslay,ngridmx,iip1,jjp1,tslab,tslabS)
    359364
    360365      call gr_fi_dyn(llm,ngridmx,iip1,jjp1,du_nonoro_gwd,du_nonoro_gwdS)

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

    r2336 r3100  
    3333
    3434      use comsoil_h, only: nsoilmx
    35       use slab_ice_h, only: noceanmx
     35!      use slab_ice_h, only: noceanmx
     36      USE ocean_slab_mod, ONLY: nslay
    3637
    3738      implicit none
    … …  
    195196        edges(1)=iip1
    196197        edges(2)=jjp1
    197         edges(3)=noceanmx
     198        edges(3)=nslay
    198199        edges(4)=1
    199200#ifdef NC_DOUBLE

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

    r2482 r3100  
    77
    88  use ioipsl_getin_p_mod, only: getin_p
     9#ifndef MESOSCALE
     10  use mod_grid_phy_lmdz, only : klon_glo ! # of physics point on full grid
     11  use mod_phys_lmdz_para, only : is_master, gather, scatter
     12#endif
    913  use watercommon_h, only: T_h2O_ice_liq, RLFTT, rhowater, mx_eau_sol
    1014  USE surfdat_h
    … …  
    1216  USE geometry_mod, only: cell_area
    1317  USE tracer_h
    14   use slab_ice_h
     18!  use slab_ice_h
     19  USE ocean_slab_mod, ONLY: alb_ice_max
     20  USE ocean_slab_mod, ONLY: alb_ice_min
     21  USE ocean_slab_mod, ONLY: h_alb_ice
    1522  use callkeys_mod, only: albedosnow,alb_ocean,albedoco2ice,ok_slab_ocean,Tsaldiff,maxicethick,co2cond
    1623  use radinc_h, only : L_NSPECTV
    … …  
    100107!$OMP THREADPRIVATE(ivap,iliq,iice)
    101108
    102       logical, save :: firstcall
     109      logical, save :: firstcall=.true.
    103110!$OMP THREADPRIVATE(firstcall)
    104111
    105       data firstcall /.true./
     112     real :: runoffamount(ngrid)
     113!#ifdef CPP_PARA
     114      real :: runoffamount_glo(klon_glo)
     115      real :: zqsurf_iliq_glo(klon_glo)
     116      real :: rnat_glo(klon_glo)
     117      real :: oceanarea_glo
     118      real :: cell_area_glo(klon_glo)
     119!#else
     120!      real :: runoffamount_glo(ngrid)
     121!      real :: zqsurf_iliq_glo(ngrid)
     122!#endif
    106123
    107124
    … …  
    137154!                      Soon to be extended to the entire water cycle...
    138155
    139 !     Total ocean surface area
     156!     LOCAL ocean surface area
    140157         oceanarea=0.
    141158         do ig=1,ngrid
    … …  
    160177      endif
    161178
     179!       call writediagfi(ngrid,"hydrol_rnat"," "," ",2,rnat)
     180!!       write (*,*) "oceanarea", oceanarea
     181
    162182!     add physical tendencies already calculated
    163183!     ------------------------------------------
    … …  
    170190         enddo   
    171191      enddo
     192
     193!      call writediagfi(ngrid,"hydrol_zqsurf_iliq_1"," "," ",2,zqsurf(1:ngrid,iliq))
    172194 
    173195      do ig=1,ngrid
    … …  
    294316!     deal with runoff
    295317            if(activerunoff)then
    296 
    297318               runoff(ig) = max(zqsurf(ig,iliq) - mx_eau_sol, 0.0)
    298319               if(ngrid.gt.1)then ! runoff only exists in 3D
    … …  
    331352      end do ! ig=1,ngrid
    332353
     354!      call writediagfi(ngrid,"hydrol_zqsurf_iliq_2"," "," ",2,zqsurf(1:ngrid,iliq))
     355
    333356!     perform crude bulk averaging of temperature in ocean
    334357!     ----------------------------------------------------
    … …  
    363386      if(activerunoff)then
    364387
    365          totalrunoff=0.
     388!         totalrunoff=0.
    366389         do ig=1,ngrid
    367             if (nint(rnat(ig)).eq.1) then
    368                totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
    369             endif
     390            runoffamount(ig) = cell_area(ig)*runoff(ig)
     391!            if (nint(rnat(ig)).eq.1) then
     392!               totalrunoff = totalrunoff + cell_area(ig)*runoff(ig)
     393!            endif
    370394         enddo
     395
     396!    collect on the full grid
     397        call gather(runoffamount,runoffamount_glo)
     398        call gather(zqsurf(1:ngrid,iliq),zqsurf_iliq_glo)
     399        call gather(rnat,rnat_glo)
     400        call gather(cell_area,cell_area_glo)
     401
     402        if (is_master) then
     403                totalrunoff=0.
     404                oceanarea_glo=0.
     405                do ig=1,klon_glo
     406                   if (nint(rnat_glo(ig)).eq.1) then
     407                        totalrunoff = totalrunoff + runoffamount_glo(ig)
     408                   endif
     409                   if (nint(rnat_glo(ig)).eq.0) then
     410                        oceanarea_glo = oceanarea_glo + cell_area_glo(ig)
     411                   endif
     412                enddo
     413
     414                do ig=1,klon_glo
     415                   if (nint(rnat_glo(ig)).eq.0) then
     416                        zqsurf_iliq_glo(ig) = zqsurf_iliq_glo(ig) + &
     417                    totalrunoff/oceanarea_glo
     418                   endif
     419                enddo
     420
     421        endif! is_master
     422
     423!    scatter the field back on all processes
     424        call scatter(zqsurf_iliq_glo,zqsurf(1:ngrid,iliq))
     425
     426       
    371427         
    372          do ig=1,ngrid
    373             if (nint(rnat(ig)).eq.0) then
    374                zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
    375                     totalrunoff/oceanarea
    376             endif
    377          enddo
    378 
    379       endif         
    380 
     428!         do ig=1,ngrid
     429!            if (nint(rnat(ig)).eq.0) then
     430!               zqsurf(ig,iliq) = zqsurf(ig,iliq) + &
     431!                    totalrunoff/oceanarea
     432!            endif
     433!         enddo
     434
     435      endif !activerunoff         
     436
     437!      call writediagfi(ngrid,"hydrol_zqsurf_iliq_3"," "," ",2,zqsurf(1:ngrid,iliq))
    381438
    382439!     Re-add the albedo effects of CO2 ice if necessary
    … …  
    400457      enddo
    401458
     459!      call writediagfi(ngrid,"hydrol_dqs_hyd_iliq"," "," ",2,dqs_hyd(1:ngrid,iliq)) 
     460
    402461      if (activerunoff) then
    403462        call writediagfi(ngrid,'runoff','Runoff amount',' ',2,runoff)

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

    r3099 r3100  
    446446     if (is_master) write(*,*) trim(rname)//": ok_slab_ocean = ",ok_slab_ocean
    447447     ! Sanity check: for now slab oncean only works in serial mode
    448      if (ok_slab_ocean.and.is_parallel) then
    449        call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
    450      endif
    451 
    452      if (is_master) write(*,*) trim(rname)//": Use slab-sea-ice ?"
    453      ok_slab_sic=.true.         ! default value
    454      call getin_p("ok_slab_sic",ok_slab_sic)
    455      if (is_master) write(*,*) trim(rname)//": ok_slab_sic = ",ok_slab_sic
    456 
    457      if (is_master) write(*,*) trim(rname)//&
    458        ": Use heat transport for the ocean ?"
    459      ok_slab_heat_transp=.true.   ! default value
    460      call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
    461      if (is_master) write(*,*) trim(rname)//&
    462        ": ok_slab_heat_transp = ",ok_slab_heat_transp
     448!     if (ok_slab_ocean.and.is_parallel) then
     449!       call abort_physic(rname," Error: slab ocean should only be used in serial mode!",1)
     450!     endif
     451
     452!     if (is_master) write(*,*) trim(rname)//": Use slab-sea-ice ?"
     453!     ok_slab_sic=.true.         ! default value
     454!     call getin_p("ok_slab_sic",ok_slab_sic)
     455!     if (is_master) write(*,*) trim(rname)//": ok_slab_sic = ",ok_slab_sic
     456
     457!     if (is_master) write(*,*) trim(rname)//&
     458!       ": Use heat transport for the ocean ?"
     459!     ok_slab_heat_transp=.true.   ! default value
     460!     call getin_p("ok_slab_heat_transp",ok_slab_heat_transp)
     461!     if (is_master) write(*,*) trim(rname)//&
     462!       ": ok_slab_heat_transp = ",ok_slab_heat_transp
    463463
    464464! Photochemistry and chemistry in the thermosphere

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

    r1621 r3100  
    469469  USE tracer_h, only: nqtot
    470470  USE comsoil_h, only: nsoilmx
    471   USE slab_ice_h, only: noceanmx
     471!  USE slab_ice_h, only: noceanmx
     472  USE ocean_slab_mod, ONLY: nslay
    472473
    473474  IMPLICIT NONE
    … …  
    558559      ENDIF
    559560
    560       ierr=NF90_DEF_DIM(nid_restart,"ocean_layers",noceanmx,idim8)
     561      ierr=NF90_DEF_DIM(nid_restart,"ocean_layers",nslay,idim8)
    561562      IF (ierr/=NF90_NOERR) THEN
    562563        write(*,*)'open_restartphy: problem defining oceanic layer dimension '
    … …  
    646647  USE mod_grid_phy_lmdz
    647648  USE mod_phys_lmdz_para
    648   USE slab_ice_h, only: noceanmx
     649!  USE slab_ice_h, only: noceanmx
     650  USE ocean_slab_mod, ONLY: nslay
     651
    649652  IMPLICIT NONE
    650653  CHARACTER(LEN=*),INTENT(IN)    :: field_name
    … …  
    832835        endif ! of if (.not.present(time))
    833836
    834       ELSE IF (field_size==noceanmx) THEN
     837      ELSE IF (field_size==nslay) THEN
    835838        ! input is a 2D "oceanic field" array
    836839        if (.not.present(time)) then ! for a time-independent field
    … …  
    948951  USE comsoil_h, only: nsoilmx
    949952  USE mod_phys_lmdz_para, only: is_master
    950   USE slab_ice_h, only: noceanmx
     953!  USE slab_ice_h, only: noceanmx
     954  USE ocean_slab_mod, ONLY: nslay
    951955  IMPLICIT NONE
    952956     CHARACTER(LEN=*),INTENT(IN) :: var_name
    … …  
    10001004        ! We know it is an  "mlayer" kind of 1D array
    10011005        idim1d=idim3
    1002       ELSEIF (var_size==noceanmx) THEN
     1006      ELSEIF (var_size==nslay) THEN
    10031007        ! We know it is an  "mlayer" kind of 1D array
    10041008        idim1d=idim8

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

    r1682 r3100  
    1 !
    2 ! $Header: /home/cvsroot/LMDZ4/libf/phylmd/ocean_slab_mod.F90,v 1.3 2008-02-04 16:24:28 fairhead Exp $
    3 !
     1!Completed
    42MODULE ocean_slab_mod
    53!
    6 ! This module is used for both surface ocean and sea-ice when using the slab ocean,
    7 ! "ocean=slab".
    8 !
    9 
    10 
    11       use slab_ice_h
    12       use watercommon_h, only: T_h2O_ice_liq
    13       use surf_heat_transp_mod
    14       implicit none
    15 
    16 
    17 
    18 
    19 
    20  
    21   !LOGICAL, PRIVATE, SAVE  :: ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS
    22   !!$OMP THREADPRIVATE(ok_slab_sic,ok_slab_heaT_h2O_ice_liqBS)
    23   !INTEGER, PRIVATE, SAVE                           :: slab_ekman, slab_cadj
    24   !!$OMP THREADPRIVATE(slab_ekman,slab_cadj)
    25   INTEGER, PRIVATE, SAVE                           :: lmt_pas, julien, idayvrai
    26   !$OMP THREADPRIVATE(lmt_pas,julien,idayvrai)
     4!==================================================================
     5!
     6!     Purpose
     7!     -------
     8!     The dynamical slab ocean model of the Generic-PCM. It has the following features:
     9!         (a) Computes sea ice creation and evolution.
     10!         (b) Snow has thermodynamic properties.
     11!         (c) Computes oceanic horizontal transport (diffusion & surface-wind driven Ekman transport).
     12!         (d) Can be used in parallel mode.
     13!
     14!     Authors
     15!     -------
     16!     S. Bhatnagar and E. Millour (2023)
     17!     Adapted from the ocean modules of LMDZ Earth (F. Codron) and the Generic-PCM (B. Charnay, 2013).
     18!
     19!     Notes
     20!     -----
     21!     Compared to the old model, the new model has the following changes (non-exhaustive):
     22!         (a) More realistic description of sea ice creation and evolution - simultaneous
     23!             surface, side and bottom melting / freezing depending on fluxes.
     24!         (b) Snow has an effective heat capacity.
     25!         (c) Snow has "weight"; it can sink an ice block if there is too much of it.
     26!         (d) Snow can be blown off by wind.
     27!         (e) The two-layer ocean allows for convective adjustment.
     28!         (f) Diffusion can follow the Gent-McWilliams scheme + Eddy diffusivity.
     29!         (g) Can be used in parallel mode.
     30!
     31!==================================================================
     32
     33  USE dimphy, ONLY: klon
     34  USE mod_grid_phy_lmdz, ONLY: klon_glo
     35  USE mod_phys_lmdz_mpi_data, ONLY: is_mpi_root
     36
     37  IMPLICIT NONE
     38  PRIVATE
     39  PUBLIC :: ocean_slab_init, ocean_slab_ice, ocean_slab_noice, &
     40            ocean_slab_frac, ocean_slab_get_vars, ocean_slab_final
     41
     42!***********************************************************************************
     43! Global saved variables
     44!***********************************************************************************
     45  ! number of slab vertical layers
     46  INTEGER, PUBLIC, SAVE :: nslay=2
     47  !$OMP THREADPRIVATE(nslay)
     48  ! number of oceanic grid points
     49  INTEGER, PUBLIC, SAVE :: knon
     50  !$OMP THREADPRIVATE(knon)
     51  ! timestep for coupling (update slab temperature) in timesteps
    2752  INTEGER, PRIVATE, SAVE                           :: cpl_pas
    2853  !$OMP THREADPRIVATE(cpl_pas)
    29   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_seaice
    30   !$OMP THREADPRIVATE(tmp_seaice)
    31   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: tmp_tslab_loc
    32   !$OMP THREADPRIVATE(tmp_tslab_loc)
    33   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slab_bils
    34   !$OMP THREADPRIVATE(slab_bils)
    35   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE , SAVE  :: lmt_bils
    36   !$OMP THREADPRIVATE(lmt_bils)
    37   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE :: tmp_pctsrf_slab
    38   !$OMP THREADPRIVATE(tmp_pctsrf_slab)
    39   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: tmp_tslab
    40   !$OMP THREADPRIVATE(tmp_tslab)
    41   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_radsol
    42   !$OMP THREADPRIVATE(tmp_radsol)
    43   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: dt_hdiff
    44   !$OMP THREADPRIVATE(dt_hdiff)
    45   REAL, ALLOCATABLE, DIMENSION(:,:), PRIVATE, SAVE   :: dt_ekman
    46   !$OMP THREADPRIVATE(dt_ekman)
    47   REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: tmp_flux_o, tmp_flux_g
    48   !$OMP THREADPRIVATE(tmp_flux_o,tmp_flux_g)
     54  ! cyang = 1/heat capacity of top layer (rho.c.H)
     55  REAL, PRIVATE, SAVE                              :: cyang
     56  !$OMP THREADPRIVATE(cyang)
     57  ! depth of slab layers (1st or 2nd layer)
    4958  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE   :: slabh
    5059  !$OMP THREADPRIVATE(slabh)
    51   LOGICAL, PRIVATE, SAVE                           :: check = .FALSE.
    52   !$OMP THREADPRIVATE(check)
     60  ! slab temperature
     61  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: tslab
     62  !$OMP THREADPRIVATE(tslab)
     63  ! heat flux convergence due to Ekman
     64  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_ekman
     65  !$OMP THREADPRIVATE(dt_ekman)
     66  ! heat flux convergence due to horiz diffusion
     67  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_hdiff
     68  !$OMP THREADPRIVATE(dt_hdiff)
     69  ! heat flux convergence due to GM eddy advection
     70  REAL, ALLOCATABLE, DIMENSION(:,:), PUBLIC, SAVE   :: dt_gm
     71  !$OMP THREADPRIVATE(dt_gm)
     72  ! fraction of ocean covered by sea ice (sic / (oce+sic))
     73  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: fsic
     74  !$OMP THREADPRIVATE(fsic)
     75  ! temperature of the sea ice
     76  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: tice
     77  !$OMP THREADPRIVATE(tice)
     78  ! sea ice thickness, in kg/m2
     79  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: seaice
     80  !$OMP THREADPRIVATE(seaice)
     81  ! net surface heat flux, weighted by open ocean fraction
     82  ! slab_bils accumulated over cpl_pas timesteps
     83  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bils_cum
     84  !$OMP THREADPRIVATE(bils_cum)
     85  ! net heat flux into the ocean below the ice : conduction + solar radiation
     86  REAL, ALLOCATABLE, DIMENSION(:), PUBLIC, SAVE  :: slab_bilg
     87  !$OMP THREADPRIVATE(slab_bilg)
     88  ! slab_bilg cululated over cpl_pas timesteps
     89  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: bilg_cum
     90  !$OMP THREADPRIVATE(bilg_cum)
     91  ! wind stress saved over cpl_pas timesteps
     92  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: taux_cum
     93  !$OMP THREADPRIVATE(taux_cum)
     94  REAL, ALLOCATABLE, DIMENSION(:), PRIVATE, SAVE  :: tauy_cum
     95  !$OMP THREADPRIVATE(tauy_cum)
     96
     97!***********************************************************************************
     98! Parameters (could be read in def file: move to slab_init)
     99!***********************************************************************************
     100! snow and ice physical characteristics:
     101    REAL, PARAMETER :: t_freeze=271.35 ! freezing sea water temp [in K]
     102    REAL, PARAMETER :: t_melt=273.15   ! melting ice temp [in K]
     103    REAL, PARAMETER :: sno_den=300. !mean snow density [in kg/m3]
     104    REAL, PARAMETER :: ice_den=917. ! ice density [in kg/m3]
     105    REAL, PARAMETER :: sea_den=1026. ! sea water density [in kg/m3]
     106    REAL, PARAMETER :: ice_cond=2.17*ice_den !conductivity of ice [in W/(m.K) or (W.kg)/(K.m4)]
     107    REAL, PRIVATE, SAVE :: sno_cond ! conductivity of snow [in W/(m.K) or (W.kg)/(K.m4)]
     108    REAL, PARAMETER :: ice_cap=2067.   ! specific heat capacity, snow and ice [in J/(kg.K)]
     109    REAL, PARAMETER :: sea_cap=3994.   ! specific heat capacity, seawater [in J/(kg.K)]
     110    REAL, PARAMETER :: ice_lat=334000. ! freeze /melt latent heat snow and ice [in J/kg]
     111    REAL, PARAMETER :: ice_sub=2834000. ! latent heat of sublimation for snow and ice [in J/kg]
     112
     113! control of snow and ice cover & freeze / melt (heights in m converted to kg/m2)
     114    REAL, PARAMETER :: snow_min=0.05*sno_den ! critical snow height [in kg/m2]
     115    REAL, PARAMETER :: snow_wfact=0.4 ! max fraction of falling snow blown into ocean [in kg/m2]
     116    REAL, PARAMETER :: ice_frac_min=0.001
     117    REAL, PRIVATE, SAVE :: ice_frac_max ! Max ice fraction (leads)
     118    REAL, PARAMETER :: h_ice_min=0.01*ice_den ! min ice thickness [in kg/m2]
     119    REAL, PARAMETER :: h_ice_thin=0.15*ice_den ! thin ice thickness [in kg/m2]
     120    ! below ice_thin, priority is to melt lateral / grow height
     121    ! ice_thin is also height of new ice
     122    REAL, PRIVATE, SAVE :: h_ice_thick ! thin ice thickness
     123    ! above ice_thick, priority is melt height / grow lateral
     124    REAL, PARAMETER :: h_ice_new=1.*ice_den ! max height of new open ocean ice [in kg/m2]
     125    REAL, PARAMETER :: h_ice_max=10.*ice_den ! max ice height [in kg/m2]
     126
     127    REAL, PARAMETER :: epsfra=1.0E-05 ! minimial grid fraction size below which there is no ice
     128
     129    REAL, PARAMETER, PUBLIC :: capcalocean=50.*4.228e+06 ! surface heat capacity [J.K-1.m-2] (assuming 50 m slab ocean)
     130    REAL, PARAMETER, PUBLIC :: capcalseaice=5.1444e+06*0.15
     131    REAL, PARAMETER, PUBLIC :: capcalsno=2.3867e+06*0.15
     132
     133    REAL, PARAMETER, PUBLIC :: h_alb_ice=0.5*ice_den ! height for full ice albedo
     134    REAL, PARAMETER, PUBLIC :: h_sno_alb=0.02*sno_den ! height for control of snow fraction
     135
     136    REAL, PARAMETER, PUBLIC :: alb_ice_min=0.2
     137    REAL, PARAMETER, PUBLIC :: alb_ice_max=0.65
     138
     139! Horizontal transport parameters
     140   ! flag for horizontal diffusion
     141   LOGICAL, PUBLIC, SAVE :: slab_hdiff
     142   !$OMP THREADPRIVATE(slab_hdiff)
     143   ! flag for GM eddy diffusivity
     144   LOGICAL, PUBLIC, SAVE :: slab_gm
     145   !$OMP THREADPRIVATE(slab_gm)
     146   REAL, PRIVATE, SAVE    :: coef_hdiff ! coefficient for horizontal diffusion
     147   !$OMP THREADPRIVATE(coef_hdiff)
     148   ! flags for Ekman, conv adjustment
     149   LOGICAL, PUBLIC, SAVE :: slab_ekman
     150   !$OMP THREADPRIVATE(slab_ekman)
     151   INTEGER, PUBLIC, SAVE :: slab_cadj
     152   !$OMP THREADPRIVATE(slab_cadj)
     153
     154!***********************************************************************************
    53155
    54156CONTAINS
    55157!
    56 !****************************************************************************************
    57 !
    58   SUBROUTINE ocean_slab_init(ngrid,dtime, tslab_rst, seaice_rst, pctsrf_rst)
    59 
    60       use slab_ice_h
    61 
    62 
    63 
    64 
    65 ! Input variables
    66 !****************************************************************************************
    67        
    68     integer,intent(in) :: ngrid ! number of atmospherci columns
     158!***********************************************************************************
     159!
     160  SUBROUTINE ocean_slab_init(dtime, pctsrf_rst, tslab_rst, tice_rst, seaice_rst, zmasq)
     161
     162! This routine
     163! (1) allocates variables initialised from restart fields
     164! (2) allocates some other variables internal to the ocean module
     165
     166    USE ioipsl_getin_p_mod, ONLY : getin_p
     167    USE mod_phys_lmdz_transfert_para, ONLY : gather
     168    USE slab_heat_transp_mod, ONLY : ini_slab_transp
     169
     170    ! Input variables
     171!***********************************************************************************
    69172    REAL, INTENT(IN)                         :: dtime
    70173! Variables read from restart file
    71     REAL, DIMENSION(ngrid,noceanmx), INTENT(IN)  :: tslab_rst         
    72     REAL, DIMENSION(ngrid), INTENT(IN)        :: seaice_rst
    73     REAL, DIMENSION(ngrid), INTENT(IN) :: pctsrf_rst
    74 
     174    REAL, DIMENSION(klon), INTENT(IN) :: pctsrf_rst
     175    REAL, DIMENSION(klon,nslay), INTENT(IN) :: tslab_rst
     176    REAL, DIMENSION(klon), INTENT(IN) :: tice_rst
     177    REAL, DIMENSION(klon), INTENT(IN) :: seaice_rst
     178    REAL, DIMENSION(klon), INTENT(IN) :: zmasq
    75179
    76180! Local variables
    77 !****************************************************************************************
     181!************************************************************************************
    78182    INTEGER                :: error
     183    REAL, DIMENSION(klon_glo) :: zmasq_glo
    79184    CHARACTER (len = 80)   :: abort_message
    80     CHARACTER (len = 20)   :: modname = 'ocean_slab_intit'
    81 
    82 
    83     print*,'************************'
    84     print*,'SLAB OCEAN est actif, prenez precautions !'
    85     print*,'************************'   
    86 
    87 ! Lecture des parametres:
    88 !    CALL getpar('ok_slab_sic',.true.,ok_slab_sic,'glace de mer dans slab')
    89 !    CALL getpar('slab_ekman',0,slab_ekman,'type transport Ekman pour slab (0=rien)')
    90 !    CALL getpar('slab_cadj',1,slab_cadj,'type ajustement conv slab 2 couches')
    91 
    92 ! Allocate variables initialize from restart fields
    93       ALLOCATE(tmp_tslab(ngrid,noceanmx), stat = error)
    94       IF (error /= 0) THEN
    95          abort_message='Pb allocation tmp_tslab'
    96          CALL abort_physic(modname,abort_message,1)
    97       ENDIF
    98       tmp_tslab(:,:) = tslab_rst(:,:)
    99       ALLOCATE(tmp_tslab_loc(ngrid,noceanmx), stat = error)
    100       IF (error /= 0) THEN
    101          abort_message='Pb allocation tmp_tslab_loc'
    102          CALL abort_physic(modname,abort_message,1)
    103       ENDIF
    104       tmp_tslab_loc(:,:) = tslab_rst(:,:)
    105 
    106     ALLOCATE(tmp_seaice(ngrid), stat = error)
     185    CHARACTER (len = 20)   :: modname = 'ocean_slab_init'
     186
     187!***********************************************************************************
     188! Define some parameters
     189!***********************************************************************************
     190!
     191! cpl_pas  coupling period (update of tslab and ice fraction)
     192! for a calculation at each physical timestep, cpl_pas=1
     193    cpl_pas = NINT(86400./dtime * 1.0) ! une fois par jour
     194    CALL getin_p('cpl_pas',cpl_pas)
     195    print *,'cpl_pas',cpl_pas
     196! Number of slab layers
     197!    nslay=2
     198!    CALL getin_p('slab_layers',nslay)
     199    print *,'number of slab layers : ',nslay
     200! Layer thickness
     201    ALLOCATE(slabh(nslay), stat = error)
    107202    IF (error /= 0) THEN
    108        abort_message='Pb allocation tmp_seaice'
     203       abort_message='Pb allocation slabh'
    109204       CALL abort_physic(modname,abort_message,1)
    110205    ENDIF
    111     tmp_seaice(:) = seaice_rst(:)
    112 
    113     ALLOCATE(tmp_pctsrf_slab(ngrid), stat = error)
     206    slabh(1)=50. ! Height of first ocean layer (wind-mixed layer)
     207    CALL getin_p('slab_depth',slabh(1))
     208    IF (nslay.GT.1) THEN
     209        slabh(2)=150. ! Height of second ocean layer (deep ocean layer)
     210    END IF
     211! cyang = 1/heat capacity of top layer (rho.c.H)
     212    cyang=1/(slabh(1)*sea_den*sea_cap)
     213
     214! ********** Sea Ice parameters ***********
     215    ice_frac_max = 0.99 ! frac = 1 may lead to some problems.
     216    CALL getin_p('ice_frac_max',ice_frac_max)
     217    h_ice_thick = 1.5
     218    CALL getin_p('h_ice_thick',h_ice_thick)
     219    h_ice_thick = h_ice_thick * ice_den
     220    sno_cond = 0.31
     221    CALL getin_p('sno_cond',sno_cond)
     222    sno_cond = sno_cond * sno_den
     223
     224! ********** Heat Transport parameters ****
     225! Ekman transport
     226!    slab_ekman=0
     227    slab_ekman=.FALSE.
     228    CALL getin_p('slab_ekman',slab_ekman)
     229! GM eddy advection (2-layers only)
     230    slab_gm=.FALSE.
     231    CALL getin_p('slab_gm',slab_gm)
     232!    IF (slab_ekman.LT.2) THEN
     233    IF (.NOT.slab_ekman) THEN
     234       slab_gm=.FALSE.
     235    ENDIF
     236! Horizontal diffusion
     237    slab_hdiff=.FALSE.
     238    CALL getin_p('slab_hdiff',slab_hdiff)
     239    IF (slab_gm) THEN
     240       coef_hdiff=8000. ! non-dimensional; coef_hdiff should be 25000 if GM is off
     241    ELSE
     242       coef_hdiff=25000.
     243    ENDIF
     244    CALL getin_p('coef_hdiff',coef_hdiff)
     245
     246! Convective adjustment
     247!    IF (nslay.EQ.1) THEN
     248!        slab_cadj=0
     249!    ELSE
     250        slab_cadj=1
     251!    END IF
     252    CALL getin_p('slab_cadj',slab_cadj)
     253
     254!************************************************************************************
     255! Allocate surface fraction read from restart file
     256!************************************************************************************
     257    ALLOCATE(fsic(klon), stat = error)
    114258    IF (error /= 0) THEN
    115259       abort_message='Pb allocation tmp_pctsrf_slab'
    116260       CALL abort_physic(modname,abort_message,1)
    117261    ENDIF
    118     tmp_pctsrf_slab(:) = pctsrf_rst(:)
    119    
    120 ! Allocate some other variables internal in module mod_oceanslab
    121     ALLOCATE(tmp_radsol(ngrid), stat = error)
     262    fsic(:)=0.
     263    !zmasq = continent fraction
     264    WHERE (1.-zmasq(:)>EPSFRA)
     265!        fsic(:) = MIN(pctsrf_rst(:,is_sic)/(1.-zmasq(:)),ice_frac_max)
     266        fsic(:) = MIN(pctsrf_rst(:)/(1.-zmasq(:)),ice_frac_max)
     267    END WHERE
     268
     269!************************************************************************************
     270! Allocate saved fields
     271!************************************************************************************
     272    ALLOCATE(tslab(klon,nslay), stat=error)
     273       IF (error /= 0) CALL abort_physic &
     274         (modname,'pb allocation tslab', 1)
     275       tslab(:,:) = tslab_rst(:,:)
     276
     277    ALLOCATE(bils_cum(klon), stat = error)
    122278    IF (error /= 0) THEN
    123        abort_message='Pb allocation tmp_radsol'
     279       abort_message='Pb allocation slab_bils_cum'
    124280       CALL abort_physic(modname,abort_message,1)
    125281    ENDIF
    126 
    127     ALLOCATE(tmp_flux_o(ngrid), stat = error)
    128     IF (error /= 0) THEN
    129        abort_message='Pb allocation tmp_flux_o'
    130        CALL abort_physic(modname,abort_message,1)
     282    bils_cum(:) = 0.0
     283
     284!    IF (version_ocean=='sicINT') THEN ! interactive sea ice
     285        ALLOCATE(slab_bilg(klon), stat = error)
     286        IF (error /= 0) THEN
     287           abort_message='Pb allocation slab_bilg'
     288           CALL abort_physic(modname,abort_message,1)
     289        ENDIF
     290        slab_bilg(:) = 0.0
     291        ALLOCATE(bilg_cum(klon), stat = error)
     292        IF (error /= 0) THEN
     293           abort_message='Pb allocation slab_bilg_cum'
     294           CALL abort_physic(modname,abort_message,1)
     295        ENDIF
     296        bilg_cum(:) = 0.0
     297        ALLOCATE(tice(klon), stat = error)
     298        IF (error /= 0) THEN
     299           abort_message='Pb allocation slab_tice'
     300           CALL abort_physic(modname,abort_message,1)
     301        ENDIF
     302        tice(:) = tice_rst(:)
     303        ALLOCATE(seaice(klon), stat = error)
     304        IF (error /= 0) THEN
     305           abort_message='Pb allocation slab_seaice'
     306           CALL abort_physic(modname,abort_message,1)
     307        ENDIF
     308        seaice(:) = seaice_rst(:)
     309!    END IF
     310
     311    IF (slab_hdiff) THEN !horizontal diffusion
     312        ALLOCATE(dt_hdiff(klon,nslay), stat = error)
     313        IF (error /= 0) THEN
     314           abort_message='Pb allocation dt_hdiff'
     315           CALL abort_physic(modname,abort_message,1)
     316        ENDIF
     317        dt_hdiff(:,:) = 0.0
    131318    ENDIF
    132    
    133     ALLOCATE(tmp_flux_g(ngrid), stat = error)
    134     IF (error /= 0) THEN
    135        abort_message='Pb allocation tmp_flux_g'
    136        CALL abort_physic(modname,abort_message,1)
     319
     320    IF (slab_gm) THEN !GM advection
     321        ALLOCATE(dt_gm(klon,nslay), stat = error)
     322        IF (error /= 0) THEN
     323           abort_message='Pb allocation dt_gm'
     324           CALL abort_physic(modname,abort_message,1)
     325        ENDIF
     326        dt_gm(:,:) = 0.0
    137327    ENDIF
    138328
    139 ! a mettre un slab_bils aussi en force !!!
    140     ALLOCATE(slab_bils(ngrid), stat = error)
    141     IF (error /= 0) THEN
    142        abort_message='Pb allocation slab_bils'
    143        CALL abort_physic(modname,abort_message,1)
     329!    IF (slab_ekman.GT.0) THEN ! ekman transport
     330    IF (slab_ekman) THEN ! ekman transport
     331        ALLOCATE(dt_ekman(klon,nslay), stat = error)
     332        IF (error /= 0) THEN
     333           abort_message='Pb allocation dt_ekman'
     334           CALL abort_physic(modname,abort_message,1)
     335        ENDIF
     336        dt_ekman(:,:) = 0.0
     337        ALLOCATE(taux_cum(klon), stat = error)
     338        IF (error /= 0) THEN
     339           abort_message='Pb allocation taux_cum'
     340           CALL abort_physic(modname,abort_message,1)
     341        ENDIF
     342        taux_cum(:) = 0.0
     343        ALLOCATE(tauy_cum(klon), stat = error)
     344        IF (error /= 0) THEN
     345           abort_message='Pb allocation tauy_cum'
     346           CALL abort_physic(modname,abort_message,1)
     347        ENDIF
     348        tauy_cum(:) = 0.0
    144349    ENDIF
    145     slab_bils(:) = 0.0   
    146 
    147     ALLOCATE(dt_hdiff(ngrid,noceanmx), stat = error)
    148     IF (error /= 0) THEN
    149        abort_message='Pb allocation dt_hdiff'
    150        CALL abort_physic(modname,abort_message,1)
     350
     351! Initialize transport
     352    IF (slab_hdiff.OR.(slab_ekman)) THEN
     353      CALL gather(zmasq,zmasq_glo)
     354! Master thread/process only
     355!$OMP MASTER
     356      IF (is_mpi_root) THEN
     357          CALL ini_slab_transp(zmasq_glo)
     358      END IF
     359!$OMP END MASTER
     360    END IF
     361
     362 END SUBROUTINE ocean_slab_init
     363!
     364!***********************************************************************************
     365!
     366  SUBROUTINE ocean_slab_frac(pctsrf_chg, zmasq)
     367
     368! This routine sends back the sea ice and ocean fraction to the main physics routine.
     369! Called only with interactive sea ice.
     370
     371! Arguments
     372!************************************************************************************
     373    REAL, DIMENSION(klon), INTENT(IN)                        :: zmasq   ! proxy of rnat
     374    REAL, DIMENSION(klon), INTENT(OUT) :: pctsrf_chg  ! sea-ice fraction
     375
     376       pctsrf_chg(:)=fsic(:)*(1.-zmasq(:))
     377
     378  END SUBROUTINE ocean_slab_frac
     379!
     380!************************************************************************************
     381!
     382  SUBROUTINE ocean_slab_noice(itime, dtime, knon, knindex, &
     383       precip_snow, tsurf_in, &
     384       radsol, snow, fluxsens, &
     385       tsurf_new, flux_u1, flux_v1, zmasq)
     386
     387    USE slab_heat_transp_mod, ONLY: divgrad_phy,slab_ekman2,slab_gmdiff
     388    USE mod_phys_lmdz_para
     389
     390! This routine
     391! (1) computes surface turbulent fluxes over points with some open ocean
     392! (2) reads additional Q-flux (everywhere)
     393! (3) computes horizontal transport (diffusion & Ekman)
     394! (4) updates slab temperature every cpl_pas ; creates new ice if needed.
     395
     396! Note :
     397! klon total number of points
     398! knon number of points with open ocean (varies with time)
     399! knindex gives position of the knon points within klon.
     400! In general, local saved variables have klon values
     401! variables exchanged with PBL module have knon.
     402
     403! Input arguments
     404!***********************************************************************************
     405    INTEGER, INTENT(IN)                  :: itime ! current timestep INTEGER,
     406    INTEGER, INTENT(IN)                  :: knon  ! number of points
     407    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
     408    REAL, INTENT(IN) :: dtime  ! timestep (s)
     409    REAL, DIMENSION(klon), INTENT(IN)    :: precip_snow !, precip_rain
     410
     411    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in ! surface temperature
     412    REAL, DIMENSION(klon), INTENT(IN) :: radsol ! net surface (radiative) flux
     413    REAL, DIMENSION(klon), INTENT(IN)   :: flux_u1, flux_v1 ! Comes from turbdiff
     414    REAL, DIMENSION(klon), INTENT(IN)   :: fluxsens !, sensible heat flux
     415    REAL, DIMENSION(klon), INTENT(IN)   :: zmasq   ! proxy of rnat
     416
     417! In/Output arguments
     418!************************************************************************************
     419    REAL, DIMENSION(klon), INTENT(INOUT) :: snow ! in kg/m2
     420
     421! Output arguments
     422!************************************************************************************
     423    REAL, DIMENSION(klon), INTENT(OUT)   :: tsurf_new ! new surface tempearture
     424
     425! Local variables
     426!************************************************************************************
     427    INTEGER               :: i,ki,k
     428    REAL                  :: t_cadj
     429
     430    ! for new ice creation
     431    REAL                  :: e_freeze
     432    REAL, DIMENSION(klon) :: slab_bils ! weighted surface heat flux
     433    ! horizontal diffusion and Ekman local vars
     434    ! dimension = global domain (klon_glo) instead of // subdomains
     435    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_glo,dt_ekman_glo,dt_gm_glo
     436    ! dt_ekman_glo saved for diagnostic, dt_ekman_tmp used for time loop
     437    REAL, DIMENSION(klon_glo,nslay) :: dt_hdiff_tmp, dt_ekman_tmp
     438    REAL, DIMENSION(klon_glo,nslay) :: tslab_glo
     439    REAL, DIMENSION(klon_glo) :: taux_glo,tauy_glo
     440
     441!****************************************************************************************
     442! 1) Surface fluxes calculation
     443!    Points with some open ocean only
     444!****************************************************************************************
     445! save total cumulated heat fluxes locally
     446! radiative + turbulent + melt of falling snow
     447    slab_bils(:)=0.
     448    DO i=1,knon
     449        ki=knindex(i)
     450        slab_bils(ki)=(1.-fsic(ki))*(fluxsens(ki)+radsol(ki) &
     451                      -precip_snow(ki)*ice_lat*(1.+snow_wfact*fsic(ki)))
     452        bils_cum(ki)=bils_cum(ki)+slab_bils(ki)
     453    END DO
     454
     455!  Compute surface wind stress
     456!    CALL calcul_flux_wind(knon, dtime, &
     457!         u0, v0, u1, v1, gustiness, cdragm, &
     458!         AcoefU, AcoefV, BcoefU, BcoefV, &
     459!         p1lay, temp_air, &
     460!         flux_u1, flux_v1)
     461
     462! save cumulated wind stress
     463    IF (slab_ekman) THEN
     464      DO i=1,knon
     465          ki=knindex(i)
     466          taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*(1.-fsic(ki))/cpl_pas
     467          tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*(1.-fsic(ki))/cpl_pas
     468      END DO
    151469    ENDIF
    152     dt_hdiff = 0.0   
    153 
    154     ALLOCATE(dt_ekman(ngrid,noceanmx), stat = error)
    155     IF (error /= 0) THEN
    156        abort_message='Pb allocation dt_hdiff'
    157        CALL abort_physic(modname,abort_message,1)
     470
     471!****************************************************************************************
     472! 2) Q-Flux : get global variables lmt_bils, diff_sst and diff_siv from file
     473! limit_slab.nc
     474!
     475!****************************************************************************************
     476!    CALL limit_slab(itime, dtime, jour, lmt_bils, diff_sst, diff_siv)
     477    ! lmt_bils and diff_sst,siv saved by limit_slab
     478    ! qflux = total QFlux correction (in W/m2)
     479!    IF (qflux_bnd.EQ.2) THEN
     480!        dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400.
     481!        dt_qflux_sic(:) = -diff_siv(:)*ice_den*ice_lat/86400.
     482!    ELSE
     483!        dt_qflux(:,1) = lmt_bils(:,1)+diff_sst(:)/cyang/86400.  &
     484!                  - diff_siv(:)*ice_den*ice_lat/86400.
     485!    END IF
     486!    IF (nslay.GT.1) THEN
     487!       dt_qflux(:,2:nslay)=lmt_bils(:,2:nslay)
     488!    END IF
     489
     490!****************************************************************************************
     491! 3) Recalculate new temperature (add Surf fluxes, Q-Flux, Ocean transport)
     492!    Bring to freezing temp and make sea ice if necessary
     493!
     494!***********************************************o*****************************************
     495    tsurf_new=tsurf_in
     496    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab & fraction
     497! ***********************************
     498!  Horizontal transport
     499! ***********************************
     500      IF (slab_ekman) THEN
     501          ! copy wind stress to global var
     502          CALL gather(taux_cum,taux_glo)
     503          CALL gather(tauy_cum,tauy_glo)
     504      END IF
     505
     506      IF (slab_hdiff.OR.(slab_ekman)) THEN
     507        CALL gather(tslab,tslab_glo)
     508      ! Compute horiz transport on one process only
     509        IF (is_mpi_root .AND. is_omp_root) THEN ! Only master processus
     510          IF (slab_hdiff) THEN
     511              dt_hdiff_glo(:,:)=0.
     512          END IF
     513          IF (slab_ekman) THEN
     514              dt_ekman_glo(:,:)=0.
     515          END IF
     516          IF (slab_gm) THEN
     517              dt_gm_glo(:,:)=0.
     518          END IF
     519          DO i=1,cpl_pas ! time splitting for numerical stability
     520!            IF (slab_ekman.GT.0) THEN
     521!              SELECT CASE (slab_ekman)
     522!                CASE (1) ! 1.5 layer scheme
     523!                  CALL slab_ekman1(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp)
     524!                CASE (2) ! 2-layers
     525!                  CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
     526!                CASE DEFAULT
     527!                  dt_ekman_tmp(:,:)=0.
     528!              END SELECT
     529            IF (slab_ekman) THEN
     530              CALL slab_ekman2(taux_glo,tauy_glo,tslab_glo,dt_ekman_tmp,dt_hdiff_tmp,slab_gm)
     531
     532              dt_ekman_glo(:,:)=dt_ekman_glo(:,:)+dt_ekman_tmp(:,:)
     533              ! convert dt_ekman from K.s-1.(kg.m-2) to K.s-1
     534              DO k=1,nslay
     535                dt_ekman_tmp(:,k)=dt_ekman_tmp(:,k)/(slabh(k)*sea_den)
     536              ENDDO
     537              tslab_glo=tslab_glo+dt_ekman_tmp*dtime
     538              IF (slab_gm) THEN ! Gent-McWilliams eddy advection
     539                dt_gm_glo(:,:)=dt_gm_glo(:,:)+ dt_hdiff_tmp(:,:)
     540                ! convert dt from K.s-1.(kg.m-2) to K.s-1
     541                DO k=1,nslay
     542                  dt_hdiff_tmp(:,k)=dt_hdiff_tmp(:,k)/(slabh(k)*sea_den)
     543                END DO
     544                tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
     545              END IF
     546            ENDIF
     547            IF (slab_hdiff) THEN ! horizontal diffusion
     548              ! laplacian of slab T
     549              CALL divgrad_phy(nslay,tslab_glo,dt_hdiff_tmp)
     550              ! multiply by diff coef and normalize to 50m slab equivalent
     551              dt_hdiff_tmp=dt_hdiff_tmp*coef_hdiff*50./SUM(slabh)
     552              dt_hdiff_glo(:,:)=dt_hdiff_glo(:,:)+ dt_hdiff_tmp(:,:)
     553              tslab_glo=tslab_glo+dt_hdiff_tmp*dtime
     554            END IF
     555          END DO ! time splitting
     556          IF (slab_hdiff) THEN
     557            !dt_hdiff_glo saved in W/m2
     558            DO k=1,nslay
     559              dt_hdiff_glo(:,k)=dt_hdiff_glo(:,k)*slabh(k)*sea_den*sea_cap/cpl_pas
     560            END DO
     561          END IF
     562          IF (slab_gm) THEN
     563            !dt_hdiff_glo saved in W/m2
     564            dt_gm_glo(:,:)=dt_gm_glo(:,:)*sea_cap/cpl_pas
     565          END IF
     566          IF (slab_ekman) THEN
     567            ! dt_ekman_glo saved in W/m2
     568            dt_ekman_glo(:,:)=dt_ekman_glo(:,:)*sea_cap/cpl_pas
     569          END IF
     570        END IF ! master process
     571!$OMP BARRIER
     572        ! Send new fields back to all processes
     573        CALL Scatter(tslab_glo,tslab)
     574        IF (slab_hdiff) THEN
     575          CALL Scatter(dt_hdiff_glo,dt_hdiff)
     576        END IF
     577        IF (slab_gm) THEN
     578          CALL Scatter(dt_gm_glo,dt_gm)
     579        END IF
     580        IF (slab_ekman) THEN
     581          CALL Scatter(dt_ekman_glo,dt_ekman)
     582          ! clear wind stress
     583          taux_cum(:)=0.
     584          tauy_cum(:)=0.
     585        END IF
     586      ENDIF ! transport
     587
     588! ***********************************
     589! Other heat fluxes
     590! ***********************************
     591      ! Add cumulated ocean surface fluxes
     592      tslab(:,1) = tslab(:,1) + bils_cum(:)*cyang*dtime
     593      ! Convective adjustment if 2 layers
     594      IF ((nslay.GT.1).AND.(slab_cadj.GT.0)) THEN
     595        DO i=1,klon
     596          IF (tslab(i,2).GT.tslab(i,1)) THEN
     597            ! mean (mass-weighted) temperature
     598            t_cadj=SUM(tslab(i,:)*slabh(:))/SUM(slabh(:))
     599            tslab(i,1)=t_cadj
     600            tslab(i,2)=t_cadj
     601          END IF
     602        END DO
     603      END IF
     604      ! Add read QFlux
     605!      IF (qflux_bnd.EQ.1) THEN
     606!          ! QFlux from ocean circ. cannot cool tslab below freezing.
     607!          dq_freeze = (t_freeze - tslab(:,1)) / (cyang*dtime*cpl_pas)
     608!          dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.))
     609!      ELSE IF (qflux_bnd.EQ.2) THEN
     610!          dq_freeze = (t_freeze - tslab(:,1)) / (cyang*dtime*cpl_pas)  &
     611!                     - dt_qflux_sic(:)
     612!          dt_qflux(:,1) = MAX(dt_qflux(:,1), MIN(dq_freeze,0.))  &
     613!                       + dt_qflux_sic(:)
     614!      END IF
     615!      DO k=1,nslay
     616!          tslab(:,k) = tslab(:,k) + dt_qflux(:,k)*cyang*dtime*cpl_pas &
     617!                     * slabh(1)/slabh(k)
     618!      END DO
     619
     620! ***********************************
     621! Update surface temperature and ice
     622! ***********************************
     623!      SELECT CASE(version_ocean)
     624!      CASE('sicNO') ! no sea ice even below freezing !
     625!          DO i=1,knon
     626!              ki=knindex(i)
     627!              tsurf_new(i)=tslab(ki,1)
     628!          END DO
     629!      CASE('sicOBS') ! "realistic" case, for prescribed sea ice
     630!        ! tslab cannot be below freezing
     631!          DO i=1,knon
     632!              ki=knindex(i)
     633!              IF (tslab(ki,1).LT.t_freeze) THEN
     634!                  tslab(ki,1)=t_freeze
     635!              END IF
     636!              tsurf_new(i)=tslab(ki,1)
     637!          END DO
     638!      CASE('sicINT') ! interactive sea ice
     639          DO i=1,knon
     640              ki=knindex(i)
     641              ! Check if new slab temperature is below freezing
     642              IF (tslab(ki,1).LT.t_freeze) THEN
     643                  ! We need to form ice now over ice-free points
     644                  ! Else points not seen by slab_ice
     645                  IF (fsic(ki)*(1.-zmasq(ki)).LT.epsfra) THEN
     646                      ! No ice present yet.
     647                      ! quantity of new ice formed
     648                      e_freeze=(t_freeze-tslab(ki,1))/cyang/ice_lat  &
     649                               +fsic(ki)*seaice(ki)
     650                      ! new ice forms at h_ice_thin
     651                      tsurf_new(ki)=t_freeze
     652                      tice(ki)=t_freeze
     653                      fsic(ki)=MIN(ice_frac_max,e_freeze/h_ice_thin)
     654                      IF (fsic(ki).GT.ice_frac_min) THEN
     655                          seaice(ki)=MIN(e_freeze/fsic(ki),h_ice_max)
     656                          tslab(ki,1)=t_freeze
     657                      ELSE
     658                          fsic(ki)=0.
     659                      END IF
     660                  END IF ! sea ice present
     661                  ! if ice already present, new ice formed in slab_ice routine.
     662!              ELSE ! temperature above freezing
     663!                  tsurf_new(i)=tslab(ki,1)
     664              END IF
     665          END DO
     666!      END SELECT
     667      bils_cum(:)=0.0! clear cumulated fluxes
     668    END IF ! coupling time
     669  END SUBROUTINE ocean_slab_noice
     670!
     671!*****************************************************************************
     672
     673!  SUBROUTINE ocean_slab_ice(   &
     674!       itime, dtime, jour, knon, knindex, &
     675!       tsurf_in, p1lay, cdragh, cdragm, precip_rain, precip_snow, temp_air, spechum, &
     676!       AcoefH, AcoefQ, BcoefH, BcoefQ, &
     677!       AcoefU, AcoefV, BcoefU, BcoefV, &
     678!       ps, u1, v1, gustiness, &
     679!       radsol, snow, qsurf, qsol, agesno, &
     680!       alb1_new, alb2_new, evap, fluxsens, fluxlat, flux_u1, flux_v1, &
     681!       tsurf_new, dflux_s, dflux_l, swnet)
     682
     683  SUBROUTINE ocean_slab_ice(itime, dtime, knon, knindex, &
     684       precip_snow, tsurf_in, &
     685       radsol, snow, fluxsens, &
     686       tsurf_new, evap, flux_u1, flux_v1, zmasq)
     687
     688!   USE calcul_fluxs_mod
     689
     690!   INCLUDE "YOMCST.h"
     691!   INCLUDE "clesphys.h"
     692
     693! This routine
     694! (1) computes surface turbulent fluxes over points with some sea ice
     695! (2) computes condutive fluxes in the snow and ice, and ice-ocean flux
     696! (3) computes snow/ice albedo
     697! (4) updates snow/ice temperature, surface melt if needed.
     698! (5) lateral growth if tslab < freezing
     699! (6) bottom & side melt / growth depending on bottom fluxes
     700! (7) updates slab temperature every cpl_pas
     701
     702! Note :
     703! klon total number of points
     704! knon number of points with open ocean (varies with time)
     705! knindex gives position of the knon points within klon.
     706! In general, local saved variables have klon values
     707! variables exchanged with PBL module have knon.
     708
     709! Input arguments
     710!****************************************************************************************
     711    INTEGER, INTENT(IN)                  :: itime, knon !, jour
     712    INTEGER, DIMENSION(klon), INTENT(IN) :: knindex
     713    REAL, INTENT(IN)                     :: dtime
     714    REAL, DIMENSION(klon), INTENT(IN)    :: tsurf_in
     715!    REAL, DIMENSION(klon), INTENT(IN)    :: p1lay
     716!    REAL, DIMENSION(klon), INTENT(IN)    :: cdragh, cdragm
     717    REAL, DIMENSION(klon), INTENT(IN)    :: precip_snow !, precip_rain
     718    REAL, DIMENSION(klon), INTENT(IN)   :: evap, fluxsens!,fluxlat
     719    REAL, DIMENSION(klon), INTENT(IN)   :: flux_u1, flux_v1
     720    REAL, DIMENSION(klon), INTENT(IN)   :: zmasq   ! proxy of rnat
     721!    REAL, DIMENSION(klon), INTENT(IN)    :: spechum, temp_air
     722!    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefH, AcoefQ, BcoefH, BcoefQ
     723!    REAL, DIMENSION(klon), INTENT(IN)    :: AcoefU, AcoefV, BcoefU, BcoefV
     724!    REAL, DIMENSION(klon), INTENT(IN)    :: ps
     725!    REAL, DIMENSION(klon), INTENT(IN)    :: u1, v1, gustiness
     726!    REAL, DIMENSION(klon), INTENT(IN)    :: swnet
     727
     728! In/Output arguments
     729!****************************************************************************************
     730    REAL, DIMENSION(klon), INTENT(INOUT)          :: snow!, qsol
     731!    REAL, DIMENSION(klon), INTENT(INOUT)          :: agesno
     732    REAL, DIMENSION(klon), INTENT(INOUT)          :: radsol
     733
     734! Output arguments
     735!****************************************************************************************
     736!    REAL, DIMENSION(klon), INTENT(OUT)            :: qsurf
     737!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
     738!    REAL, DIMENSION(klon), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
     739!    REAL, DIMENSION(klon), INTENT(OUT)            :: evap, fluxsens!, fluxlat
     740!    REAL, DIMENSION(klon), INTENT(OUT)            :: flux_u1, flux_v1
     741    REAL, DIMENSION(klon), INTENT(OUT)            :: tsurf_new
     742!    REAL, DIMENSION(klon), INTENT(OUT)            :: dflux_s, dflux_l
     743
     744! Local variables
     745!****************************************************************************************
     746    INTEGER               :: i,ki
     747!    REAL, DIMENSION(klon) :: cal, beta, dif_grnd
     748!    REAL, DIMENSION(klon) :: u0, v0
     749!    REAL, DIMENSION(klon) :: u1_lay, v1_lay
     750    REAL, DIMENSION(klon) :: f_bot ! bottom ocean - ice flux
     751
     752    ! intermediate heat fluxes:
     753    ! (conduction snow / ice, shortwave penetration, bottom turbulent fluxes)
     754    REAL                  :: f_cond_s, f_cond_i, f_turb
     755    ! friction velocity, 1/C and 1/tau conduction for ice
     756    REAL                  :: ustar
     757!    REAL                  :: uscap, ustau
     758    ! for snow/ice albedo:
     759!    REAL                  :: alb_snow, alb_ice, alb_pond
     760!    REAL                  :: frac_snow, frac_ice, frac_pond
     761!    REAL                  :: z1_i, z2_i, z1_s, zlog ! height parameters
     762    ! for ice melt / freeze
     763    REAL                  :: e_melt, e_freeze, snow_evap, h_test, h_new
     764    ! dhsic, dfsic change in ice mass, fraction.
     765    ! frac_mf ratio of lateral / thickness growth / melt
     766    REAL                  :: dhsic, dfsic, frac_mf
     767
     768!*******************************************************************************
     769! 1) Update surface , ice and slab temperature
     770!*******************************************************************************
     771! Wind stress
     772! flux_u1, flux_v1 from physics
     773! save cumulated wind stress
     774! Use ocean-ice stress = wind stress * (1.-fsic)
     775!    IF (slab_ekman.GT.0) THEN
     776    IF (slab_ekman) THEN
     777      DO i=1,knon
     778          ki=knindex(i)
     779          taux_cum(ki)=taux_cum(ki)+flux_u1(ki)*fsic(ki)*(1.-fsic(ki))/cpl_pas
     780          tauy_cum(ki)=tauy_cum(ki)+flux_v1(ki)*fsic(ki)*(1.-fsic(ki))/cpl_pas
     781      END DO
    158782    ENDIF
    159     dt_ekman = 0.0   
    160 
    161 
    162     ALLOCATE(lmt_bils(ngrid), stat = error)
    163     IF (error /= 0) THEN
    164        abort_message='Pb allocation lmt_bils'
    165        CALL abort_physic(modname,abort_message,1)
     783
     784! Initialize ice-ocean flux
     785    slab_bilg(:)=0.
     786
     787    ! Old, explicit scheme for snow & ice conductive fluxes
     788    ! radsol is total surface fluxes (input) : radiative + turbulent
     789        DO i=1,knon
     790        ki=knindex(i) ! For PCM : you can probably replace ki by i
     791            ! ocean-ice turbulent heat flux
     792            ! turbulent velocity = (tau/rho)^1/2
     793            ustar = MAX(5e-4, SQRT((1.-fsic(ki))   &
     794                   * SQRT(flux_u1(ki)**2 + flux_v1(ki)**2) / sea_den))
     795            f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze)
     796            ! f_turb >0 and cannot bring tslab below zero
     797            f_turb = MAX(0., MIN(f_turb, &
     798                        (tslab(ki,1)-t_freeze) / (cyang*dtime*cpl_pas)))
     799
     800            ! ice conductive flux (pos up)
     801            IF (seaice(ki).GT.0) THEN
     802                f_cond_i = ice_cond*(t_freeze-tice(ki))/seaice(ki)
     803            ELSE
     804                f_cond_i = 0
     805            END IF
     806
     807            ! If snow layer present, tsurf = tsnow
     808            ! Problem here, as tsurf_in # tsnow ?
     809            IF (snow(ki).GT.snow_min) THEN
     810                ! snow conductive flux (pos up)
     811                f_cond_s=sno_cond*(tice(ki)-tsurf_in(ki))/snow(ki)
     812                ! update ice temperature
     813                tice(ki)=tice(ki) + 2./ice_cap/(snow(ki)+seaice(ki)) &
     814                *(f_cond_i-f_cond_s)*dtime
     815                ! update snow temperature
     816                tsurf_new(ki) = tsurf_in(ki) + 2./ice_cap/snow(ki) &
     817                    *(fluxsens(ki)+radsol(ki)+f_cond_s)*dtime
     818            ELSE IF (seaice(ki).GT.0) THEN ! bare ice. tsurf = tice
     819                ! update ice temperature
     820                        tice(ki) = tice(ki) + 2./ice_cap/seaice(ki) &
     821                                *(fluxsens(ki)+radsol(ki)+f_cond_i)*dtime
     822                        tsurf_new(ki) = tice(ki)
     823            END IF
     824            ! bottom flux (used to grow ice from below)
     825            f_bot(ki) = f_turb - f_cond_i
     826            slab_bilg(ki) = -f_turb
     827        END DO
     828!
     829!! Surface turbulent fluxes (sens, lat etc) and update surface temp.
     830!    dif_grnd(:)=0.
     831!    beta(:) = 1.
     832!    CALL calcul_fluxs(knon, is_sic, dtime, &
     833!        tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
     834!        precip_rain, precip_snow, snow, qsurf,  &
     835!        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
     836!        f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
     837!        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     838!    DO i=1,knon
     839!        IF (snow(i).LT.snow_min) tice(knindex(i))=tsurf_new(i)
     840!    END DO
     841
     842! Surface, snow-ice and ice-ocean fluxes.
     843! Prepare call to calcul_fluxs (cal, beta, radsol, dif_grnd)
     844
     845! Surface turbulent fluxes (sens, lat etc) and update surface temp.
     846!    beta(:) = 1.
     847!    CALL calcul_fluxs(knon, is_sic, dtime, &
     848!        tsurf_in, p1lay, cal, beta, cdragh, cdragh, ps, &
     849!        precip_rain, precip_snow, snow, qsurf,  &
     850!        radsol, dif_grnd, temp_air, spechum, u1_lay, v1_lay, gustiness, &
     851!        f_qsat_oce,AcoefH, AcoefQ, BcoefH, BcoefQ, &
     852!        tsurf_new, evap, fluxlat, fluxsens, dflux_s, dflux_l)
     853
     854!! Update remaining temperature and fluxes
     855!    DO i=1,knon
     856!    ki=knindex(i)
     857!        ! ocean-ice turbulent heat flux
     858!        ! turbulent velocity = (tau/rho)^1/2 for low ice cover
     859!        ! min = 5e-4 for 1cm/s current
     860!        ustar = MAX(5e-4, SQRT((1.-fsic(ki))   &
     861!               * SQRT(flux_u1(i)**2 + flux_v1(i)**2) / sea_den))
     862!        f_turb = 0.0057 * sea_den * sea_cap * ustar * (tslab(ki,1) - t_freeze)
     863!        ! ice_ocean fluxes cannot bring tslab below freezing
     864!        f_turb = MAX(0., MIN(f_turb, slab_bilg(ki) + (tslab(ki,1)-t_freeze) &
     865!                          / (fsic(ki)*cyang*dtime*cpl_pas)))
     866!        IF (snow(i).GT.snow_min) THEN
     867!            ! snow conductive flux after calcul_fluxs
     868!            f_cond_s = sno_cond * (tice(ki)-tsurf_new(i)) / snow(i)
     869!            ! 1 / heat capacity and conductive timescale
     870!            uscap = 2. / ice_cap / (snow(i)+seaice(ki))
     871!            ustau = uscap * ice_cond / seaice(ki)
     872!            ! update ice temp
     873!            tice(ki) = (tice(ki) + dtime*(ustau*t_freeze - uscap*f_cond_s)) &
     874!                     / (1 + dtime*ustau)
     875!        ELSE ! bare ice
     876!            tice(ki)=tsurf_new(i)
     877!        END IF
     878!        ! ice conductive flux (pos up)
     879!        f_cond_i = ice_cond * (t_freeze-tice(ki)) / seaice(ki)
     880!        f_bot(i) = f_turb - f_cond_i
     881!        slab_bilg(ki) = slab_bilg(ki)-f_turb
     882!    END DO
     883
     884    ! weight fluxes to ocean by sea ice fraction
     885    slab_bilg(:)=slab_bilg(:)*fsic(:)
     886
     887!****************************************************************************************
     888! 2) Update snow and ice surface : thickness and fraction
     889!****************************************************************************************
     890    DO i=1,knon
     891        ki=knindex(i)
     892! snow precip (could be before fluxes above ?)
     893        IF (precip_snow(ki) > 0.) THEN
     894            snow(ki) = snow(ki)+precip_snow(ki)*dtime*(1.-snow_wfact*(1.-fsic(ki)))
     895        END IF
     896! snow and ice sublimation
     897        IF (evap(ki) > 0.) THEN
     898           snow_evap = MIN (snow(ki) / dtime, evap(ki))
     899           snow(ki) = snow(ki) - snow_evap * dtime
     900           snow(ki) = MAX(0.0, snow(ki))
     901           seaice(ki) = MAX(0.0,seaice(ki)-(evap(ki)-snow_evap)*dtime)
     902        ENDIF
     903! Melt / Freeze snow from above if Tsurf>0
     904        IF (tsurf_new(ki).GT.t_melt) THEN
     905            ! energy available for melting snow (in kg of melted snow /m2)
     906            e_melt = MIN(MAX(snow(ki)*(tsurf_new(ki)-t_melt)*ice_cap/2. &
     907               /(ice_lat+ice_cap/2.*(t_melt-tice(ki))),0.0),snow(ki))
     908            ! remove snow
     909            IF (snow(ki).GT.e_melt) THEN
     910                snow(ki)=snow(ki)-e_melt
     911                tsurf_new(ki)=t_melt
     912            ELSE ! all snow is melted
     913                ! add remaining heat flux to ice
     914                e_melt=e_melt-snow(ki)
     915                snow(ki)=0.0
     916                tice(ki)=tice(ki)+e_melt*ice_lat*2./(ice_cap*seaice(ki))
     917                tsurf_new(ki)=tice(ki)
     918            END IF
     919        END IF
     920! Bottom melt / grow
     921! bottom freeze if bottom flux (cond + oce-ice) <0
     922        IF (f_bot(ki).LT.0) THEN
     923           ! larger fraction of bottom growth
     924           frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thick)   &
     925                  / (h_ice_max-h_ice_thick)))
     926           ! quantity of new ice (formed at mean ice temp)
     927           e_melt= -f_bot(ki) * dtime * fsic(ki) &
     928                   / (ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     929           ! first increase height to h_thick
     930           dhsic=MAX(0.,MIN(h_ice_thick-seaice(ki),e_melt/fsic(ki)))
     931           seaice(ki)=dhsic+seaice(ki)
     932           e_melt=e_melt-fsic(ki)*dhsic
     933           IF (e_melt.GT.0.) THEN
     934           ! frac_mf fraction used for lateral increase
     935           dfsic=MIN(ice_frac_max-fsic(ki),e_melt*frac_mf/seaice(ki))
     936           fsic(ki)=fsic(ki)+dfsic
     937           e_melt=e_melt-dfsic*seaice(ki)
     938           ! rest used to increase height
     939           seaice(ki)=MIN(h_ice_max,seaice(ki)+e_melt/fsic(ki))
     940           END IF
     941       ELSE
     942! melt from below if bottom flux >0
     943           ! larger fraction of lateral melt from warm ocean
     944           frac_mf=MIN(1.,MAX(0.,(seaice(ki)-h_ice_thin)   &
     945                  / (h_ice_thick-h_ice_thin)))
     946           ! bring ice to freezing and melt from below
     947           ! quantity of melted ice
     948           e_melt= f_bot(ki) * dtime * fsic(ki) &
     949                   / (ice_lat+ice_cap/2.*(tice(ki)-t_freeze))
     950           ! first decrease height to h_thick
     951           IF (fsic(ki).GT.0) THEN
     952             dhsic=MAX(0.,MIN(seaice(ki)-h_ice_thick,e_melt/fsic(ki)))
     953           ELSE
     954             dhsic=0
     955           ENDIF
     956           seaice(ki)=seaice(ki)-dhsic
     957           e_melt=e_melt-fsic(ki)*dhsic
     958           IF (e_melt.GT.0) THEN
     959           ! frac_mf fraction used for height decrease
     960           dhsic=MAX(0.,MIN(seaice(ki)-h_ice_min,e_melt*frac_mf/fsic(ki)))
     961           seaice(ki)=seaice(ki)-dhsic
     962           e_melt=e_melt-fsic(ki)*dhsic
     963           ! rest used to decrease fraction (up to 0!)
     964           dfsic=MIN(fsic(ki),e_melt/seaice(ki))
     965           ! keep remaining in ocean if everything melted
     966           e_melt=e_melt-dfsic*seaice(ki)
     967           slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
     968           ELSE
     969           dfsic=0
     970           END IF
     971           fsic(ki)=fsic(ki)-dfsic
     972       END IF
     973! melt ice from above if Tice>0
     974        IF (tice(ki).GT.t_melt) THEN
     975            ! quantity of ice melted (kg/m2)
     976            e_melt=MAX(seaice(ki)*(tice(ki)-t_melt)*ice_cap/2. &
     977             /(ice_lat+ice_cap/2.*(t_melt-t_freeze)),0.0)
     978            ! melt from above, height only
     979            dhsic=MIN(seaice(ki)-h_ice_min,e_melt)
     980            e_melt=e_melt-dhsic
     981            IF (e_melt.GT.0) THEN
     982              ! lateral melt if ice too thin
     983              dfsic=MAX(fsic(ki)-ice_frac_min,e_melt/h_ice_min*fsic(ki))
     984              ! if all melted add remaining heat to ocean
     985              e_melt=MAX(0.,e_melt*fsic(ki)-dfsic*h_ice_min)
     986              slab_bilg(ki) = slab_bilg(ki) + e_melt*ice_lat/dtime
     987              ! update height and fraction
     988              fsic(ki)=fsic(ki)-dfsic
     989            END IF
     990            seaice(ki)=seaice(ki)-dhsic
     991            ! surface temperature at melting point
     992            tice(ki)=t_melt
     993            tsurf_new(ki)=t_melt
     994        END IF
     995! convert snow to ice if below floating line
     996        h_test=(seaice(ki)+snow(ki))*ice_den-seaice(ki)*sea_den
     997        IF ((h_test.GT.0.).AND.(seaice(ki).GT.h_ice_min)) THEN !snow under water
     998            ! extra snow converted to ice (with added frozen sea water)
     999            dhsic=h_test/(sea_den-ice_den+sno_den)
     1000            seaice(ki)=seaice(ki)+dhsic
     1001            snow(ki)=snow(ki)-dhsic*sno_den/ice_den
     1002            ! available energy (freeze sea water + bring to tice)
     1003            e_melt=dhsic*(1.-sno_den/ice_den)*(ice_lat+ &
     1004                   ice_cap/2.*(t_freeze-tice(ki)))
     1005            ! update ice temperature
     1006            tice(ki)=tice(ki)+2.*e_melt/ice_cap/(snow(ki)+seaice(ki))
     1007        END IF
     1008    END DO
     1009
     1010!*******************************************************************************
     1011! 3) cumulate ice-ocean fluxes, update tslab, lateral grow
     1012!***********************************************o*******************************
     1013    !cumul fluxes.
     1014    bilg_cum(:)=bilg_cum(:)+slab_bilg(:)
     1015    IF (MOD(itime,cpl_pas).EQ.0) THEN ! time to update tslab
     1016        ! Add cumulated surface fluxes
     1017        tslab(:,1)=tslab(:,1)+bilg_cum(:)*cyang*dtime
     1018        bilg_cum(:)=0.
     1019! If slab temperature below freezing, new lateral growth
     1020        DO i=1,knon
     1021            ki=knindex(i)
     1022            IF (tslab(ki,1).LT.t_freeze) THEN ! create more ice
     1023                ! quantity of new ice formed over open ocean
     1024                ! (formed at mean ice temperature)
     1025                e_freeze=(t_freeze-tslab(ki,1))/cyang &
     1026                         /(ice_lat+ice_cap/2.*(t_freeze-tice(ki)))
     1027                ! new ice height and fraction
     1028                h_new=MAX(h_ice_thin,MIN(h_ice_new,seaice(ki))) ! new height
     1029!                h_new=MIN(h_ice_new,seaice(ki))
     1030                dfsic=MIN(ice_frac_max-fsic(ki)      &
     1031                          ,MAX(ice_frac_min,e_freeze/h_new))
     1032                ! update average sea ice height
     1033                seaice(ki)=(seaice(ki)*fsic(ki)+e_freeze) &
     1034                           / (dfsic+fsic(ki))
     1035                ! update snow thickness
     1036                snow(ki) = snow(ki) * fsic(ki) / (dfsic+fsic(ki))
     1037                ! update tslab to freezing
     1038                tslab(ki,1)=t_freeze
     1039                ! update sea ice fraction
     1040                fsic(ki)=fsic(ki)+dfsic
     1041            END IF ! tslab below freezing
     1042        END DO
     1043    END IF ! coupling time
     1044
     1045!****************************************************************************************
     1046! 4) Compute sea-ice and snow albedo
     1047!****************************************************************************************
     1048! Removed all albedo stuff as it is computed through hydrol in the Generic model
     1049! ------ End Albedo ----------
     1050
     1051    !tests remaining ice fraction
     1052    WHERE ((fsic.LT.ice_frac_min).OR.(seaice.LT.h_ice_min))
     1053        tslab(:,1)=tslab(:,1)-fsic*seaice*ice_lat*cyang
     1054        tice=t_melt
     1055        fsic=0.
     1056        seaice=0.
     1057    END WHERE
     1058
     1059  END SUBROUTINE ocean_slab_ice
     1060!
     1061!****************************************************************************************
     1062!
     1063  SUBROUTINE ocean_slab_final
     1064
     1065!****************************************************************************************
     1066! Deallocate module variables
     1067!****************************************************************************************
     1068    IF (ALLOCATED(tslab)) DEALLOCATE(tslab)
     1069    IF (ALLOCATED(fsic)) DEALLOCATE(fsic)
     1070    IF (ALLOCATED(tice)) DEALLOCATE(tice)
     1071    IF (ALLOCATED(seaice)) DEALLOCATE(seaice)
     1072    IF (ALLOCATED(slab_bilg)) DEALLOCATE(slab_bilg)
     1073    IF (ALLOCATED(bilg_cum)) DEALLOCATE(bilg_cum)
     1074    IF (ALLOCATED(bils_cum)) DEALLOCATE(bils_cum)
     1075    IF (ALLOCATED(taux_cum)) DEALLOCATE(taux_cum)
     1076    IF (ALLOCATED(tauy_cum)) DEALLOCATE(tauy_cum)
     1077    IF (ALLOCATED(dt_ekman)) DEALLOCATE(dt_ekman)
     1078    IF (ALLOCATED(dt_hdiff)) DEALLOCATE(dt_hdiff)
     1079    IF (ALLOCATED(dt_gm)) DEALLOCATE(dt_gm)
     1080!    IF (ALLOCATED(dt_qflux)) DEALLOCATE(dt_qflux)
     1081!    IF (ALLOCATED(dt_qflux_sic)) DEALLOCATE(dt_qflux_sic)
     1082
     1083  END SUBROUTINE ocean_slab_final
     1084!
     1085!****************************************************************************************
     1086!
     1087  SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_g_loc, &
     1088       dt_hdiff_loc,dt_ekman_loc)
     1089
     1090! "Get some variables from module ocean_slab_mod"
     1091
     1092    INTEGER, INTENT(IN)                     :: ngrid
     1093    REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: tslab_loc
     1094    REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
     1095    REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
     1096    REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_hdiff_loc ! [in W/m2]
     1097    REAL, DIMENSION(ngrid,nslay), INTENT(OUT) :: dt_ekman_loc ! [in W/m2]
     1098    INTEGER :: i
     1099
     1100
     1101! Set the output variables
     1102    tslab_loc(:,:) = 0.
     1103    dt_hdiff_loc(:,:)=0.
     1104    dt_ekman_loc(:,:)=0.
     1105    tslab_loc(:,:) = tslab(:,:)
     1106    seaice_loc(:) = seaice(:)
     1107    flux_g_loc(:) = slab_bilg(:)
     1108
     1109!!      dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2
     1110!!      dt_ekman_loc(:,i) = dt_ekman(:,i)*slabh(i)*1000.*4228.
     1111
     1112    IF (slab_hdiff) THEN
     1113        DO i=1,nslay
     1114           dt_hdiff_loc(:,i) = dt_hdiff(:,i)
     1115        ENDDO
    1661116    ENDIF
    167     lmt_bils(:) = 0.0
    168 
    169     ALLOCATE(slabh(noceanmx), stat = error)
    170     IF (error /= 0) THEN
    171        abort_message='Pb allocation slabh'
    172        CALL abort_physic(modname,abort_message,1)
     1117    IF (slab_ekman) THEN
     1118        DO i=1,nslay
     1119           dt_ekman_loc(:,i) = dt_ekman(:,i)
     1120        ENDDO
    1731121    ENDIF
    174     slabh(1)=50.
    175     IF (noceanmx.GE.2) slabh(2)=150.
    176     IF (noceanmx.GE.3) slabh(3)=2800.
    177 
    178 
    179 !       CALL init_masquv
    180 
    181 
    182 
    183 
    184   END SUBROUTINE ocean_slab_init
    185 !
    186 !****************************************************************************************
    187 !
    188  
    189 !
    190 !****************************************************************************************
    191 !
    192   SUBROUTINE ocean_slab_ice(dtime, &
    193        ngrid, knindex, tsurf, radsol, &
    194        precip_snow, snow, evap, &
    195        fluxsens,tsurf_new, pctsrf_sic, &
    196        taux_slab,tauy_slab,icount)
    197 
    198       use slab_ice_h
    199       use callkeys_mod, only: ok_slab_sic
    200 
    201 
    202 ! Input arguments 
    203 !****************************************************************************************
    204     INTEGER, INTENT(IN)                     :: ngrid
    205     INTEGER, DIMENSION(ngrid), INTENT(IN) :: knindex
    206     REAL, INTENT(IN)                        :: dtime
    207     REAL, DIMENSION(ngrid), INTENT(IN)    :: tsurf
    208     REAL, DIMENSION(ngrid), INTENT(IN)    :: taux_slab
    209     REAL, DIMENSION(ngrid), INTENT(IN)    :: tauy_slab
    210     REAL, DIMENSION(ngrid), INTENT(IN)    :: evap, fluxsens
    211     REAL, DIMENSION(ngrid), INTENT(IN)    :: precip_snow
    212     REAL, DIMENSION(ngrid), INTENT(IN)    :: radsol
    213     INTEGER, INTENT(IN)                     :: icount
    214 
    215 
    216 !In/Output arguments
    217 !****************************************************************************************l
    218     REAL, DIMENSION(ngrid), INTENT(INOUT)          :: snow
    219 
    220 !    REAL, DIMENSION(ngridmx), INTENT(INOUT)          :: agesno
    221 
    222 
    223 ! Output arguments
    224 !****************************************************************************************
    225 !    REAL, DIMENSION(ngridmx), INTENT(OUT)            :: alb1_new  ! new albedo in visible SW interval
    226 !    REAL, DIMENSION(ngridmx), INTENT(OUT)            :: alb2_new  ! new albedo in near IR interval
    227     REAL, DIMENSION(ngrid), INTENT(OUT)            :: tsurf_new     
    228     REAL, DIMENSION(ngrid), INTENT(OUT)            :: pctsrf_sic
    229 
    230 ! Local variables
    231 !****************************************************************************************
    232     INTEGER                                 :: i
    233     REAL, DIMENSION(ngrid)                :: cal, beta, dif_grnd, capsol
    234     REAL, DIMENSION(ngrid)                :: alb_neig, alb_ice!, tsurf_temp
    235     REAL, DIMENSION(ngrid)                :: soilcap, soilflux
    236     REAL, DIMENSION(ngrid)                :: zfra
    237     REAL                                    :: snow_evap, fonte
    238     REAL, DIMENSION(ngrid,noceanmx)        :: tslab
    239     REAL, DIMENSION(ngrid)                :: seaice,tsea_ice ! glace de mer (kg/m2)
    240     REAL, DIMENSION(ngrid)                :: pctsrf_new
    241 
    242 
    243 !*****************************************************************************
    244 
    245 
    246 ! Initialization of output variables
    247 !    alb1_new(:) = 0.0
    248 
    249 !******************************************************************************
    250 ! F. Codron: 3 cas
    251 ! -Glace interactive, quantité seaice : T glace suit modèle simple
    252 ! -Pas de glace: T oce peut descendre en dessous de 0°C sans geler
    253 !*****************************************************************************
    254 
    255        pctsrf_new(:)=tmp_pctsrf_slab(:)
    256        tmp_radsol(:)=radsol(:)
    257        tmp_flux_o(:)=fluxsens(:)
    258 
    259     DO i = 1, ngrid
    260        tsurf_new(i) = tsurf(knindex(i))
    261        seaice(i)=tmp_seaice(knindex(i))
    262 
    263 
    264 
    265        IF (pctsrf_new(knindex(i)) < EPSFRA) THEN
    266           snow(i) = 0.0
    267           tsurf_new(i) = T_h2O_ice_liq - 1.8
    268           !IF (soil_model) tsoil(i,:) = T_h2O_ice_liq -1.8
    269        ENDIF
    270     ENDDO
    271     tmp_flux_g(:) = 0.0
    272     tsea_ice(:)=T_h2O_ice_liq-1.8
    273 ! Calculs flux glace/mer et glace/air
    274     IF (ok_slab_sic) THEN   
    275 !*********************************
    276 ! Calcul de beta, heat capacity
    277 !*********************************
    278 !    CALL calbeta(dtime, is_sic, knon, snow, qsol, beta, capsol, dif_grnd)
    279    
    280 !    IF ((soil_model)) THEN
    281 
    282 !    ELSE
    283 !       dif_grnd = 1.0 / tau_gl
    284 !       cal = RCPD * calice
    285 !       WHERE (snow > 0.0) cal = RCPD * calsno
    286 !    ENDIF
    287 !    tsurf_temp = tsurf_new
    288 !    beta = 1.0
    289 
    290 
    291 ! **********************************************
    292 ! Evolution avec glace interactive:
    293 ! *************************************
    294 !      tsurf_new=tsurf_temp
    295       DO i = 1, ngrid
    296         IF (pctsrf_new(knindex(i)) .GT. epsfra) THEN
    297 ! *************************************
    298 ! + Calcul Flux entre glace et océan
    299 ! *************************************
    300           tmp_flux_g(knindex(i))=(tsurf_new(i)-(T_h2O_ice_liq-1.8)) &
    301             * ice_cond*ice_den/seaice(i)
    302 
    303 ! ****************************************
    304 ! + Calcul nouvelle température de la glace
    305 ! ****************************************
    306           tsurf_new(i)=tsurf_new(i)+2*(radsol(i)+fluxsens(i) &
    307              -tmp_flux_g(knindex(i))) &
    308              /(ice_cap*seaice(i))*dtime
    309 
    310 ! ***************************************
    311 ! + Precip and evaporation of snow and ice
    312 ! ***************************************
    313 ! Add precip
    314           IF (precip_snow(i).GT.0.) THEN
    315             snow(i)=snow(i)+precip_snow(i)*dtime
    316           ENDIF
    317 ! Evaporation
    318           snow_evap=0.
    319           IF (evap(i).GT.0.) THEN
    320             snow_evap = MIN (snow(i) / dtime, evap(i))
    321             snow(i) = snow(i) - snow_evap*dtime
    322             snow(i) = MAX(0.0, snow(i))
    323             seaice(i)=MAX(0.0,seaice(i)-(evap(i)-snow_evap)*dtime)
    324           ENDIF
    325        
    326 ! *****************************************
    327 ! + Fonte neige & glace par le haut
    328 ! *****************************************
    329 ! Snow melt
    330           IF ((snow(i) > epsfra) .AND. (tsurf_new(i) >= T_h2O_ice_liq)) THEN
    331             fonte=MIN(MAX((tsurf_new(i)-T_h2O_ice_liq)*seaice(i) &
    332                   /2.*ice_cap/ice_lat,0.0),snow(i))
    333             snow(i) = MAX(0.,(snow(i)-fonte))
    334             tsurf_new(i)=tsurf_new(i)-fonte*2*ice_lat/(seaice(i)*ice_cap)
    335           ENDIF
    336           snow(i)=MIN(snow(i),3000.)
    337 ! Ice melt
    338           IF ((seaice(i) > epsfra) .AND. (tsurf_new(i) >= T_h2O_ice_liq)) THEN       
    339              fonte=seaice(i)*ice_cap*(tsurf_new(i)-T_h2O_ice_liq) &
    340                   /(2*ice_lat+ice_cap*1.8)
    341              CALL icemelt(fonte,pctsrf_new(knindex(i)),seaice(i))
    342              tmp_flux_g(knindex(i))=tmp_flux_g(knindex(i)) &
    343                                  +fonte*ice_lat/dtime
    344              tsurf_new(i)=T_h2O_ice_liq
    345           ENDIF 
    346           tmp_seaice(knindex(i))=seaice(i)
    347         ENDIF!(pctsrf_new(knindex(i)) .GT. epsfra)
    348           tsea_ice(knindex(i))=tsurf_new(i)
    349 
    350       ENDDO
    351     ENDIF
    352 ! ******************************************
    353 ! CALL interfoce:
    354 ! cumul/calcul nouvelle T_h2O_ice_liqce
    355 ! fonte/formation de glace en dessous
    356 ! ******************************************   
    357     tslab = tmp_tslab
    358    
    359     CALL interfoce_slab(ngrid, dtime, &
    360          tmp_radsol, tmp_flux_o, tmp_flux_g, tmp_pctsrf_slab, &
    361          tslab, tsea_ice, pctsrf_new,taux_slab,tauy_slab,icount)
    362    
    363     tmp_pctsrf_slab(:)=pctsrf_new(:)
    364 !    tmp_pctsrf_slab(:,is_oce)=1.0-tmp_pctsrf_slab(:) &
    365 !            -tmp_pctsrf_slab(:,is_lic)-tmp_pctsrf_slab(:,is_ter)
    366 
    367     DO i=1, ngrid
    368        tmp_tslab(knindex(i),:)=tslab(knindex(i),:)
    369        seaice(i)=tmp_seaice(knindex(i))
    370        tsurf_new(i)=tsea_ice(knindex(i))
    371 
    372     ENDDO
    373 
    374 ! ****************************
    375 ! calcul new albedo
    376 ! ****************************
    377 ! Albedo neige
    378 !    CALL albsno(klon,knon,dtime,agesno(:),alb_neig(:), precip_snow(:)) 
    379 !    WHERE (snow(1 : knon) .LT. 0.0001) agesno(1 : knon) = 0.
    380 ! Fraction neige (hauteur critique 45kg/m2~15cm)
    381 !    zfra(1:knon) = MAX(0.0,MIN(1.0,snow(1:knon)/45.0))
    382 ! Albedo glace
    383 !    IF (ok_slab_sicOBS) THEN
    384 !       alb_ice=0.6
    385 !    ELSE
    386 !       alb_ice(1:knon)=alb_ice_max-(alb_ice_max-alb_ice_min) &
    387 !         *exp(-seaice(1:knon)/h_alb_ice)
    388 !    ENDIF
    389 !Albedo final
    390 !    alb1_new(1 : knon) = alb_neig(1 : knon) *zfra(1:knon) + &
    391 !         alb_ice * (1.0-zfra(1:knon))
    392 !    alb2_new(:) = alb1_new(:)
    393 
    394        
    395 ! ********************************
    396 ! Return the fraction of sea-ice
    397 ! ********************************
    398     pctsrf_sic(:) =  tmp_pctsrf_slab(:)
    399 
    400 
    401   END SUBROUTINE ocean_slab_ice
    402 
    403 
    404 !
    405 !***************************************************************************
    406 !
    407   SUBROUTINE interfoce_slab(ngrid, dtime, &
    408        radsol, fluxo, fluxg, pctsrf, &
    409        tslab, tsea_ice, pctsrf_slab, &
    410        taux_slab, tauy_slab,icount)
    411 !
    412 ! Cette routine calcule la temperature d'un slab ocean, la glace de mer
    413 ! et les pourcentages de la maille couverte par l'ocean libre et/ou
    414 ! la glace de mer pour un "slab" ocean de 50m
    415 !
    416 ! Conception: Laurent Li
    417 ! Re-ecriture + adaptation LMDZ4: I. Musat
    418 ! Transport, nouveau modèle glace, 2 couches: F.Codron
    419 !
    420 ! input:
    421 !   klon         nombre T_h2O_ice_liqtal de points de grille
    422 !   itap         numero du pas de temps
    423 !   dtime        pas de temps de la physique (en s)
    424 !   ijour        jour dans l'annee en cours
    425 !   radsol       rayonnement net au sol (LW + SW)
    426 !   fluxo        flux turbulent (sensible + latent) sur les mailles oceaniques
    427 !   fluxg        flux de conduction entre la surface de la glace de mer et l'ocean
    428 !   pctsrf       tableau des pourcentages de surface de chaque maille
    429 ! output:
    430 !   tslab        temperature de l'ocean libre
    431 !   tsea_ice         temperature de la glace (surface)
    432 !   pctsrf_slab  "pourcentages" (valeurs entre 0. et 1.) surfaces issus du slab
    433 
    434     use slab_ice_h
    435     use callkeys_mod, only: ok_slab_sic,ok_slab_heat_transp
    436 
    437 ! Input arguments
    438 !****************************************************************************************
    439     INTEGER, INTENT(IN)                       :: ngrid,icount
    440 !    INTEGER, INTENT(IN)                       :: itap
    441     REAL, INTENT(IN)                          :: dtime       ! not used
    442 !    INTEGER, INTENT(IN)                       :: ijour
    443     REAL, DIMENSION(ngrid), INTENT(IN)         :: radsol
    444     REAL, DIMENSION(ngrid), INTENT(IN)         :: fluxo
    445     REAL, DIMENSION(ngrid), INTENT(IN)         :: fluxg
    446     REAL, DIMENSION(ngrid), INTENT(IN)  :: pctsrf
    447     REAL, DIMENSION(ngrid), INTENT(IN)  :: taux_slab
    448     REAL, DIMENSION(ngrid), INTENT(IN)  :: tauy_slab
    449 
    450 ! Output arguments
    451 !****************************************************************************************
    452     REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: tslab
    453     REAL, DIMENSION(ngrid), INTENT(INOUT)       :: pctsrf_slab
    454     REAL, DIMENSION(ngrid), INTENT(INOUT)       :: tsea_ice
    455 
    456 ! Local variables
    457 !****************************************************************************************
    458     REAL                    :: fonte,t_cadj
    459     INTEGER                 :: i,k,kt,kb
    460     REAL                    :: zz, za, zb
    461     REAL                    :: cyang ! capacite calorifique slab J/(m2 K)
    462     REAL, PARAMETER         :: tau_conv=432000. ! temps ajust conv (5 jours)
    463     REAL, DIMENSION(ngrid,noceanmx) :: dtdiff_loc,dtekman_loc
    464 
    465 
    466 
    467 !***************************************************
    468 ! Capacite calorifique de la couche de surface
    469     cyang=capcalocean!slabh(1)*4.228e+06
    470     cpl_pas=1!4*iradia
    471 
    472 !
    473 ! lecture du bilan au sol lmt_bils issu d'une simulation forcee en debut de journee
    474 !!    julien = MOD(ijour,360)
    475 !!    IF (ok_slab_heaT_h2O_ice_liqBS .and. (MOD(itap,lmt_pas) .EQ. 1)) THEN 
    476 !!       ! 1er pas de temps de la journee
    477 !!       idayvrai = ijour
    478 !!       CALL condsurf(julien,idayvrai, lmt_bils)
    479 !!    ENDIF
    480 
    481 ! ----------------------------------------------
    482 ! A chaque pas de temps: cumul flux de chaleur
    483 ! ----------------------------------------------
    484 ! bilan du flux de chaleur dans l'ocean :
    485 !
    486     DO i = 1, ngrid
    487        za = radsol(i) + fluxo(i)
    488        zb = fluxg(i)
    489 !       IF(((1-pctsrf(i)).GT.epsfra).OR. &
    490 !            (pctsrf(i).GT.epsfra)) THEN
    491           slab_bils(i)=slab_bils(i)+(za*(1-pctsrf(i)) &
    492                +zb*pctsrf(i))/ cpl_pas
    493 !       ENDIF
    494 
    495     ENDDO !klon
    496 
    497 ! ---------------------------------------------
    498 ! T_h2O_ice_liqus les cpl_pas: update de tslab et seaice
    499 ! ---------------------------------------------
    500 
    501     IF (mod(icount-1,cpl_pas).eq.0) THEN !fin de journee
    502 !
    503 ! ---------------------------------------------
    504 ! Transport de chaleur par circulation
    505 ! Decoupage par pas de temps pour stabilite numérique
    506 ! (diffusion, schema centre pour advection)
    507 ! ---------------------------------------------
    508       dt_hdiff(:,:)=0.
    509       dt_ekman(:,:)=0.
    510 
    511       IF (ok_slab_heat_transp) THEN
    512 !       DO i=1,cpl_pas
    513 !  Transport diffusif
    514 !         IF (ok_soil_hdiff) THEN
    515              CALL divgrad_phy(ngrid,noceanmx,tmp_tslab_loc,dtdiff_loc)
    516              dtdiff_loc=dtdiff_loc*soil_hdiff*50./SUM(slabh)!*100.
    517  !           dtdiff_loc(:,1)=dtdiff_loc(:,1)*soil_hdiff*50./SUM(slabh)*0.8
    518  !           dtdiff_loc(:,2)=dtdiff_loc(:,2)*soil_hdiff*50./SUM(slabh)*0.2
    519              dt_hdiff=dt_hdiff+dtdiff_loc
    520              tmp_tslab_loc=tmp_tslab_loc+dtdiff_loc*dtime
    521 !         END IF
    522 
    523 ! Calcul  de transport par Ekman
    524 
    525           CALL slab_ekman2(ngrid,taux_slab,tauy_slab,tslab,dtekman_loc)
    526 
    527 
    528 !        END SELECT
    529 !        IF (slab_ekman.GT.0) THEN
    530           do k=1,noceanmx
    531             dtekman_loc(:,k)=dtekman_loc(:,k)/(slabh(k)*1000.)!*0.
    532           enddo
    533           dt_ekman(:,:)=dt_ekman(:,:)+dtekman_loc(:,:)
    534           tmp_tslab_loc=tmp_tslab_loc+dtekman_loc*dtime
    535 !        ENDIF
    536 !      ENDDO ! time splitting 1,cpl_pas     
    537 !      IF (slab_ekman.GT.0) THEN
    538 !         taux_slab(:)=0.
    539 !         tauy_slab(:)=0.
    540 !      ENDIF
    541 
    542 !       print*, 'slab_bils=',slab_bils(1)
    543 
    544 
    545       ENDIF!(ok_slab_heat_transp)
    546 
    547       DO i = 1, ngrid
    548 !      IF (((1-pctsrf(i)).GT.epsfra).OR. &
    549 !             (pctsrf(i).GT.epsfra)) THEN
    550 ! ---------------------------------------------
    551 ! Ajout des flux de chaleur dans l'océan
    552 ! ---------------------------------------
    553 
    554 !print*, 'T_h2O_ice_liqcean1=',tmp_tslab_loc(i,1)
    555         tmp_tslab_loc(i,1) = tmp_tslab_loc(i,1) + &
    556              slab_bils(i)/cyang*dtime*cpl_pas
    557 
    558 !print*, 'capcalocean=',capcalocean
    559 !print*, 'cyang=',cyang
    560 !print*, 'dT_h2O_ice_liqcean=',slab_bils(i)/cyang*dtime
    561 !print*, 'T_h2O_ice_liqcean2=',tmp_tslab_loc(i,1)
    562 
    563 
    564 ! on remet l'accumulation a 0
    565         slab_bils(i) = 0.
    566 ! ---------------------------------------------
    567 ! Glace interactive
    568 ! ---------------------------------------------
    569         IF (ok_slab_sic) THEN
    570 ! Fondre la glace si température > 0.
    571 ! -----------------------------------
    572           IF ((tmp_tslab_loc(i,1).GT.T_h2O_ice_liq-1.8) .AND. (tmp_seaice(i).GT.0.0)) THEN
    573             fonte=(tmp_tslab_loc(i,1)-T_h2O_ice_liq+1.8)*cyang &
    574                 /(ice_lat+ice_cap/2*(T_h2O_ice_liq-1.8-tsea_ice(i)))
    575             CALL icemelt(fonte,pctsrf_slab(i),tmp_seaice(i))
    576                  tmp_tslab_loc(i,1)=T_h2O_ice_liq-1.8+fonte*ice_lat/cyang
    577           ENDIF
    578 ! fabriquer de la glace si congelation atteinte:
    579 ! ----------------------------------------------
    580           IF (tmp_tslab_loc(i,1).LT.(T_h2O_ice_liq-1.8)) THEN
    581 
    582             IF (tmp_seaice(i).LT.h_ice_min) THEN
    583 ! Creation glace nouvelle
    584 !             IF (pctsrf_slab(i).LT.(epsfra)) THEN
    585                  fonte=(T_h2O_ice_liq-1.8-tmp_tslab_loc(i,1))*cyang/ice_lat
    586               IF (fonte.GT.h_ice_min*ice_frac_min) THEN
    587                 tmp_seaice(i)=MIN(h_ice_thin,fonte/ice_frac_min)
    588                 tmp_seaice(i)=MAX(tmp_seaice(i),fonte/ice_frac_max)
    589 !             IF (fonte.GT.0) THEN
    590 !               tmp_seaice(i)=fonte
    591                 tsea_ice(i)=T_h2O_ice_liq-1.8
    592                 pctsrf_slab(i)=(1-pctsrf_slab(i))*fonte/tmp_seaice(i)
    593 !                pctsrf_slab(i)=1.
    594                 tmp_tslab_loc(i,1)=T_h2O_ice_liq-1.8
    595               ENDIF
    596             ELSE 
    597 ! glace déjà présente
    598 ! Augmenter epaisseur
    599               fonte=(T_h2O_ice_liq-1.8-tmp_tslab_loc(i,1))*cyang &
    600                    /(ice_lat+ice_cap/2*(T_h2O_ice_liq-1.8-tsea_ice(i)))           
    601               zz=tmp_seaice(i)
    602               tmp_seaice(i)=MAX(zz,MIN(h_ice_thick,fonte+zz))
    603 ! Augmenter couverture (oce libre et h>h_thick)
    604               za=1-pctsrf_slab(i)
    605               zb=pctsrf_slab(i)
    606               fonte=fonte*za+MAX(0.0,fonte+zz-tmp_seaice(i))*zb
    607               pctsrf_slab(i)=MIN(zb+fonte/tmp_seaice(i), &
    608                                     (za+zb)*ice_frac_max)
    609               fonte=MAX(0.0,fonte-(pctsrf_slab(i)-zb)*tmp_seaice(i))
    610 ! Augmenter epaisseur si couverture complete
    611               tmp_seaice(i)=tmp_seaice(i)+fonte/pctsrf_slab(i)
    612               tmp_tslab_loc(i,1) = T_h2O_ice_liq-1.8
    613             ENDIF ! presence glace
    614           ENDIF !congelation
    615 ! vérifier limites de hauteur de glace
    616           IF(tmp_seaice(i).GT.h_ice_min) THEN
    617             tmp_seaice(i) = MIN(tmp_seaice(i),h_ice_max)
    618           ELSE
    619             tmp_seaice(i) = 0.
    620             pctsrf_slab(i)=0.
    621           ENDIF! (tmp_seaice(i).GT.h_ice_min)
    622        
    623         ENDIF !(ok_slab_sic) Glace interactive
    624 ! ----------------------------------
    625 ! Ajustement convectif si > 1 layers
    626 ! ----------------------------------
    627 
    628         IF ((noceanmx.GT.1)) THEN
    629           DO kt=1,noceanmx-1
    630             kb=kt           
    631             DO k=kt+1,noceanmx !look for instability
    632               IF (tmp_tslab_loc(i,k).GT.tmp_tslab_loc(i,kt)) THEN
    633                 kb=k
    634               ENDIF
    635             ENDDO
    636             IF (kb.GT.kt) THEN !ajust conv
    637 !               IF (slab_cadj.EQ.1) THEN
    638 !             :   t_cadj=SUM(tmp_tslab_loc(i,kt:kb)*slabh(kt:kb))/SUM(slabh(kt:kb))
    639 !                DO k=kt,kb
    640 !                  tmp_tslab_loc(i,k)=t_cadj
    641 !                ENDDO
    642 !              ELSEIF (slab_cadj.EQ.2) THEN
    643                 t_cadj=(tmp_tslab_loc(i,kb)-tmp_tslab_loc(i,kt))*dtime/tau_conv*cpl_pas
    644                 tmp_tslab_loc(i,kt)=tmp_tslab_loc(i,kt)+t_cadj
    645                 tmp_tslab_loc(i,kb)=tmp_tslab_loc(i,kb)-t_cadj*slabh(kt)/slabh(kb)
    646               !ENDIF
    647             ENDIF
    648           ENDDO
    649         ENDIF !ajust 2 layers
    650 !      ENDIF !pctsrf
    651       ENDDO !klon
    652 
    653 ! On met a jour la temperature et la glace
    654     tslab  = tmp_tslab_loc
    655 
    656 
    657     ENDIF !(mod(icount-1,cpl_pas).eq.0)
    658 
    659   END SUBROUTINE interfoce_slab
    660 !
    661 !****************************************************************************** 
    662   SUBROUTINE icemelt(fonte,pctsrf,seaice)
    663 
    664       use slab_ice_h
    665 
    666 
    667 
    668      REAL, INTENT(INOUT)  :: pctsrf
    669      REAL , INTENT(INOUT)   ::fonte,seaice !kg/m2
    670      REAL :: hh !auxilliary
    671      REAL :: ff !auxilliary
    672 
    673 
    674 ! ice gt h_ice_thick: decrease thickness up T_h2O_ice_liq h_ice_thick
    675      IF (seaice.GT.h_ice_thick) THEN
    676         hh=seaice
    677 !        ff=fonte
    678         seaice=max(h_ice_thick,hh-fonte)   
    679         fonte=max(0.0,fonte+h_ice_thick-hh)
    680 
    681 !        seaice=max(0.,hh-fonte)   
    682 !        fonte=max(0.0,fonte-(seaice-hh))
    683 !     IF (seaice.LT.epsfra) THEN
    684 !        pctsrf=0.
    685 !        seaice=0.
    686 !        fonte=ff-hh
    687 !     ENDIF
    688 
    689      ENDIF
    690 ! ice gt h_ice_thin: partially decrease thickness
    691      IF ((seaice.GE.h_ice_thin).AND.(fonte.GT.0.0)) THEN
    692        hh=seaice
    693        seaice=MAX(hh-0.6*fonte,h_ice_thin)
    694        fonte=MAX(0.0,fonte-hh+seaice)
    695      ENDIF
    696 ! use rest T_h2O_ice_liq decrease area
    697        hh=pctsrf
    698        pctsrf=MIN(hh,MAX(0.0,hh-fonte/seaice))
    699        fonte=MAX(0.0,fonte-(hh-pctsrf)*seaice)
    700 
    701     END SUBROUTINE icemelt
    702 
    703 !****************************************************************************************
    704 !
    705   SUBROUTINE ocean_slab_final!(tslab_rst, seaice_rst)
    706 
    707 ! This subroutine will send T_h2O_ice_liq phyredem the variables concerning the slab
    708 ! ocean that should be written T_h2O_ice_liq restart file.
    709 
    710 !****************************************************************************************
    711 
    712 !    REAL, DIMENSION(ngridmx,noceanmx), INTENT(OUT) :: tslab_rst
    713 !    REAL, DIMENSION(ngridmx), INTENT(OUT) :: seaice_rst
    714 
    715 !****************************************************************************************
    716 ! Set the output variables
    717 !    tslab_rst(:,:)  = tmp_tslab(:,:)
    718 !    tslab_rst(:)  = tmp_tslab_loc(:)
    719 !    seaice_rst(:) = tmp_seaice(:)
    720 
    721 ! Deallocation of all variables in module
    722     DEALLOCATE(tmp_tslab, tmp_tslab_loc, tmp_pctsrf_slab)
    723     DEALLOCATE(tmp_seaice, tmp_radsol, tmp_flux_o, tmp_flux_g)
    724     DEALLOCATE(slab_bils, lmt_bils)
    725     DEALLOCATE(dt_hdiff)
    726 
    727   END SUBROUTINE ocean_slab_final
    728 !
    729 !****************************************************************************************
    730 !
    731   SUBROUTINE ocean_slab_get_vars(ngrid,tslab_loc, seaice_loc, flux_o_loc, flux_g_loc, &
    732        dt_hdiff_loc,dt_ekman_loc)
    733  
    734 ! "Get some variables from module ocean_slab_mod"
    735 ! This subroutine prints variables T_h2O_ice_liq a external routine
    736 
    737     INTEGER, INTENT(IN)                     :: ngrid
    738     REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: tslab_loc
    739     REAL, DIMENSION(ngrid), INTENT(OUT) :: seaice_loc
    740     REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_o_loc
    741     REAL, DIMENSION(ngrid), INTENT(OUT) :: flux_g_loc
    742     REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: dt_hdiff_loc
    743     REAL, DIMENSION(ngrid,noceanmx), INTENT(OUT) :: dt_ekman_loc
    744     INTEGER :: i
    745 
    746 
    747 ! Set the output variables
    748     tslab_loc=0.
    749     dt_hdiff_loc=0.
    750     dt_ekman_loc=0.
    751     tslab_loc  = tmp_tslab
    752     seaice_loc(:) = tmp_seaice(:)
    753     flux_o_loc(:) = tmp_flux_o(:)
    754     flux_g_loc(:) = tmp_flux_g(:)
    755     DO i=1,noceanmx
    756       dt_hdiff_loc(:,i) = dt_hdiff(:,i)*slabh(i)*1000.*4228. !Convert en W/m2
    757       dt_ekman_loc(:,i) = dt_ekman(:,i)*slabh(i)*1000.*4228.
    758     ENDDO
    759  
    760  
     1122
     1123
    7611124
    7621125  END SUBROUTINE ocean_slab_get_vars
    … …  
    7641127!****************************************************************************************
    7651128!
     1129
    7661130END MODULE ocean_slab_mod

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

    r2557 r3100  
    1313! to use  'getin_p'
    1414      use ioipsl_getin_p_mod, only: getin_p
    15 
     15!!
     16  use write_field_phy, only: Writefield_phy
     17!!
    1618  use tabfi_mod, only: tabfi
    17   USE tracer_h, ONLY: noms
     19  USE tracer_h, ONLY: noms, igcm_h2o_vap
    1820  USE surfdat_h, only: phisfi, albedodat, zmea, zstd, zsig, zgam, zthe
    1921  use iostart, only: nid_start, open_startphy, close_startphy, &
    2022                     get_field, get_var, inquire_field, &
    2123                     inquire_dimension, inquire_dimension_length
    22   use slab_ice_h, only: noceanmx
     24!  use slab_ice_h, only: noceanmx
     25  USE ocean_slab_mod, ONLY: nslay
    2326  use callkeys_mod, only: CLFvarying,surfalbedo,surfemis, callsoil
    2427  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
    … …  
    5053  real,intent(out) :: cloudfrac(ngrid,nlayer)
    5154  real,intent(out) :: hice(ngrid), totcloudfrac(ngrid)
    52   real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,noceanmx) 
     55  real,intent(out) :: pctsrf_sic(ngrid),tslab(ngrid,nslay) 
    5356  real,intent(out) :: tsea_ice(ngrid),sea_ice(ngrid)
    5457  real,intent(out) :: rnat(ngrid)
     58!  real,intent(out) :: tice(ngrid)
    5559
    5660!======================================================================
    … …  
    314318  if (.not.found) then
    315319    write(*,*) "phyetat0: Failed loading <tslab>"
    316     do iq=1,noceanmx
     320    do iq=1,nslay
    317321      tslab(1:ngrid,iq)=tsurf(1:ngrid)
    318322    enddo
    319323  endif
    320324else
    321   do iq=1,noceanmx
     325  do iq=1,nslay
    322326    tslab(1:ngrid,iq)=tsurf(1:ngrid)
    323327  enddo
    … …  
    338342write(*,*) "phyetat0: Oceanic ice temperature <tsea_ice> range:", &
    339343             minval(tsea_ice), maxval(tsea_ice)
     344
     345!if (startphy_file) then
     346  ! Oceanic ice temperature
     347!  call get_field("tice",tice,found,indextime)
     348!  if (.not.found) then
     349!    write(*,*) "phyetat0: Failed loading <tice>"
     350!    tice(1:ngrid)=273.15-1.8
     351!  endif
     352!else
     353!  tice(1:ngrid)=273.15-1.8
     354!endif ! of if (startphy_file)
     355!write(*,*) "phyetat0: Oceanic ice temperature <tice> range:", &
     356!             minval(tice), maxval(tice)
    340357
    341358if (startphy_file) then
    … …  
    384401endif ! of if (nq.ge.1)
    385402
     403!!call WriteField_phy("in_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    386404
    387405if (startphy_file) then

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

    r2299 r3100  
    141141                      put_var, put_field
    142142  use tracer_h, only: noms
    143   use slab_ice_h, only: noceanmx
     143!  use slab_ice_h, only: noceanmx
     144  USE ocean_slab_mod, ONLY: nslay
    144145  use callkeys_mod, only: ok_slab_ocean, calllott_nonoro
    145146  use nonoro_gwd_ran_mod, only: du_nonoro_gwd, dv_nonoro_gwd, &
    … …  
    165166  real,intent(in) :: rnat(ngrid)
    166167  real,intent(in) :: pctsrf_sic(ngrid)
    167   real,intent(in) :: tslab(ngrid,noceanmx)
     168  real,intent(in) :: tslab(ngrid,nslay)
    168169  real,intent(in) :: tsea_ice(ngrid)
    169170  real,intent(in) :: sea_ice(ngrid)
     171!  real,intent(in) :: tice(ngrid)
    170172
    171173  integer :: iq
    … …  
    200202   call put_field("pctsrf_sic","Pourcentage sea ice",pctsrf_sic)
    201203   call put_field("tslab","Temperature slab ocean",tslab)
    202    call put_field("tsea_ice","Temperature sea ice",tsea_ice)
     204   call put_field("tsea_ice","Temperature of top layer (seaice or snow)",tsea_ice)
     205!   call put_field("tice","Temperature of seaice if snow on top",tice)
    203206   call put_field("sea_ice","Sea ice mass",sea_ice)
    204207 endif!(ok_slab_ocean)

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

    r2714 r3100  
    1313      use comsaison_h, only: mu0, fract
    1414      use radcommon_h, only: glat
    15       use slab_ice_h, only : noceanmx
     15!      use slab_ice_h, only : noceanmx
     16      USE ocean_slab_mod, ONLY: nslay
    1617      USE radinc_h, only : L_NSPECTI, L_NSPECTV,naerkind
    1718      use surfdat_h, only: phisfi, albedodat,  &
    … …  
    9495      real,dimension(:),allocatable,save ::  pctsrf_sic
    9596      real,dimension(:,:),allocatable,save :: tslab        ! Slab_ocean temperature (K)
    96       real,dimension(:),allocatable,save ::  tsea_ice      ! Sea ice temperature (K)
     97      real,dimension(:),allocatable,save ::  tsea_ice      ! Temperature of the top layer (K), either snow or ice
     98!      real,dimension(:),allocatable,save :: tice           ! Sea ice temperature (K) if there is snow on top of it
    9799      real,dimension(:),allocatable,save :: sea_ice
    98100      real,dimension(:),allocatable,save :: zmasq
    … …  
    191193        ALLOCATE(tau_col(klon))
    192194        ALLOCATE(pctsrf_sic(klon))
    193         ALLOCATE(tslab(klon,noceanmx))
     195        ALLOCATE(tslab(klon,nslay))
    194196        ALLOCATE(tsea_ice(klon))
     197!        ALLOCATE(tice(klon))
    195198        ALLOCATE(sea_ice(klon))
    196199        ALLOCATE(zmasq(klon))
    … …  
    285288        DEALLOCATE(tslab)
    286289        DEALLOCATE(tsea_ice)
     290 !       DEALLOCATE(tice)
    287291        DEALLOCATE(sea_ice)
    288292        DEALLOCATE(zmasq)

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

    r3081 r3100  
    1313                  pdu,pdv,pdt,pdq,pdpsrf)
    1414
     15!!
     16      use write_field_phy, only: Writefield_phy 
     17!!
    1518      use ioipsl_getin_p_mod, only: getin_p
    1619      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind, corrkdir, banddir
    … …  
    3841      use wstats_mod, only: callstats, wstats, mkstats
    3942      use phyredem, only: physdem0, physdem1
    40       use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
    41                             noceanmx
    42       use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
    43                                 ini_surf_heat_transp_mod, &
    44                                 ocean_slab_get_vars,ocean_slab_final
    45       use surf_heat_transp_mod,only :init_masquv
     43!!      use slab_ice_h, only: capcalocean, capcalseaice,capcalsno, &
     44!!                            noceanmx
     45!!      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_ice, &
     46!!                                capcalocean, capcalseaice,capcalsno, nslay, &
     47!!                                ocean_slab_final
     48!!      use surf_heat_transp_mod,only :init_masquv
     49      use ocean_slab_mod, only :ocean_slab_init, ocean_slab_noice, ocean_slab_ice, &
     50                                ocean_slab_frac, ocean_slab_get_vars, &
     51                                capcalocean, capcalseaice,capcalsno, nslay, knon, &
     52                                ocean_slab_final
    4653      use planetwide_mod, only: planetwide_minval,planetwide_maxval,planetwide_sumval
    4754      use mod_phys_lmdz_para, only : is_master
    … …  
    310317      real zdtrain_generic(ngrid,nlayer)                      ! Rain_generic routine.
    311318      real dtmoist(ngrid,nlayer)                              ! Moistadj routine.
    312       real dt_ekman(ngrid,noceanmx), dt_hdiff(ngrid,noceanmx) ! Slab_ocean routine.
     319      real dt_ekman(ngrid,nslay), dt_hdiff(ngrid,nslay) ! Slab_ocean routine.
    313320      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
    314321      real zdtchim(ngrid,nlayer)                              ! Calchim routine.
    … …  
    459466
    460467      real :: tsurf2(ngrid)
    461       real :: flux_o(ngrid),flux_g(ngrid)
     468!!      real :: flux_o(ngrid),flux_g(ngrid)
     469      real :: flux_g(ngrid)
    462470      real :: flux_sens_lat(ngrid)
    463471      real :: qsurfint(ngrid,nq)
    … …  
    530538!        Read 'startfi.nc' file.
    531539!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    532 #ifndef MESOSCALE
     540#ifndef MESOSCALE 
    533541         call phyetat0(startphy_file,                                 &
    534542                       ngrid,nlayer,"startfi.nc",0,0,nsoilmx,nq,      &
    … …  
    536544                       cloudfrac,totcloudfrac,hice,                   &
    537545                       rnat,pctsrf_sic,tslab, tsea_ice,sea_ice)
    538 #else
     546
     547!!         call WriteField_phy("post_phyetat0_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     548#else   
     549
    539550         day_ini = pday
    540551#endif
    … …  
    570581         albedo_co2_ice_SPECTV(:)=0.0
    571582         call surfini(ngrid,nq,qsurf,albedo,albedo_bareground,albedo_snow_SPECTV,albedo_co2_ice_SPECTV)
     583
     584!!        call WriteField_phy("post_surfini1_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     585!!        call WriteField_phy("post_surfini2_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    572586         
    573587!        Initialize orbital calculation.
    … …  
    586600         if (ok_slab_ocean)then
    587601
    588            call ocean_slab_init(ngrid,ptimestep, tslab, &
    589                                 sea_ice, pctsrf_sic)
    590 
    591            call ini_surf_heat_transp_mod()
     602!!           call ocean_slab_init(ngrid,ptimestep, tslab, &
     603!!                                sea_ice, pctsrf_sic)
     604
     605!!           call ini_surf_heat_transp_mod()
    592606           
    593607           knindex(:) = 0
     608           knon = 0
    594609
    595610           do ig=1,ngrid
    596611              zmasq(ig)=1
    597               knindex(ig) = ig
    598               if (nint(rnat(ig)).eq.0) then
     612!              knindex(ig) = ig
     613              if (nint(rnat(ig)).eq.0) then ! rnat = 0 (ocean), 1 (continent)
    599614                 zmasq(ig)=0
     615                 knon=knon+1
     616                 knindex(knon)=ig
    600617              endif
    601618           enddo
    602619
    603            CALL init_masquv(ngrid,zmasq)
     620           call ocean_slab_init(ptimestep, pctsrf_sic, tslab, tsea_ice, sea_ice, zmasq)
     621
     622!!           CALL init_masquv(ngrid,zmasq)
    604623
    605624         endif ! end of 'ok_slab_ocean'.
    … …  
    635654           
    636655         endif ! end of 'callsoil'.
     656
     657!!         call WriteField_phy("post_callsoil_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    637658         
    638659         icount=1
    … …  
    687708         qsurf_hist(:,:)=qsurf(:,:)
    688709
     710!!         call WriteField_phy("post_qsurf_hist_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     711
    689712!        Initialize variable for dynamical heating and zonal wind tendency diagnostic
    690713!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    724747         endif
    725748
     749!!         call WriteField_phy("post_ice_update_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     750
    726751!        Set metallicity for GCS
    727752!        ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    … …  
    734759            CALL init_thermcell_mod(g, rcp, r, pi, T_h2o_ice_liq, RV)
    735760         endif
    736          
     761
    737762         call su_watercycle ! even if we don't have a water cycle, we might
    738763                            ! need epsi for the wvp definitions in callcorrk.F
    739764                            ! or RETV, RLvCp for the thermal plume model
     765
     766!!         call WriteField_phy("post_watercycle_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    740767#ifndef MESOSCALE
    741768         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
    … …  
    744771                         albedo_bareground,inertiedat,zmea,zstd,zsig,zgam,zthe)
    745772         endif
     773
     774!!         call WriteField_phy("post_physdem_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    746775#endif         
    747776         if (corrk) then
    … …  
    761790            call suaer_corrk       ! Set up aerosol optical properties.
    762791         endif
     792
     793!!         call WriteField_phy("post_corrk_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
    763794         ! XIOS outputs
    764795#ifdef CPP_XIOS
    … …  
    768799                                     year_day,presnivs,pseudoalt,WNOI,WNOV)
    769800#endif
     801         
     802!!         call WriteField_phy("post_xios_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     803
    770804         write(*,*) "physiq: end of firstcall"
    771805      endif ! end of 'firstcall'
     806
     807!!      call WriteField_phy("post_firstcall_qsurf",qsurf(1:ngrid,igcm_h2o_vap),1)
     808!!      call writediagfi(ngrid,"firstcall_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap)) 
    772809
    773810      if (check_physics_inputs) then
    … …  
    775812         call check_physics_fields("begin physiq:", pt, pu, pv, pplev, pq)
    776813      endif
     814
     815!      call writediagfi(ngrid,"pre_physical_rnat"," "," ",2,rnat)
     816!      call writediagfi(ngrid,"pre_physical_capcal"," "," ",2,capcal)
    777817
    778818! ------------------------------------------------------
    … …  
    938978            if(rings_shadow) then
    939979                call call_rings(ngrid, ptime, pday, diurnal)
    940             endif   
     980            endif   
     981
     982!!            call writediagfi(ngrid,"corrk_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     983!!            call writediagfi(ngrid,"corrk_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    941984
    942985
    … …  
    11421185      endif ! of if (callrad)
    11431186
     1187!!      call writediagfi(ngrid,"vdifc_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1188!!      call writediagfi(ngrid,"vdifc_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    11441189
    11451190
    … …  
    11941239         endif
    11951240
     1241!!         call writediagfi(ngrid,"vdifc_post_zdqsdif"," "," ",2,zdqsdif(1:ngrid,igcm_h2o_vap))
    11961242
    11971243         if (tracer) then
    … …  
    12001246         end if ! of if (tracer)
    12011247
     1248!!         call writediagfi(ngrid,"vdifc_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1249!!         call writediagfi(ngrid,"vdifc_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    12021250
    12031251         ! test energy conservation
    … …  
    14121460         dqsurf(1:ngrid,igcm_co2_ice) = dqsurf(1:ngrid,igcm_co2_ice) + zdqsurfc(1:ngrid)
    14131461
     1462!!         call writediagfi(ngrid,"condense_co2_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1463!!         call writediagfi(ngrid,"condense_co2_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1464
    14141465         ! test energy conservation
    14151466         if(enertest)then
    … …  
    14451496                  ! Moist Convective Adjustment.
    14461497                  ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     1498!!!!
     1499!                  call writediagfi(ngrid,"pre_moistadj_pt"," "," ",3,pt)
     1500!                  call writediagfi(ngrid,"pre_moistadj_pq_h2o_ice"," "," ",3,pq(:,:,igcm_h2o_ice))
     1501!                  call writediagfi(ngrid,"pre_moistadj_pdq_h2o_ice"," "," ",3,pdq(:,:,igcm_h2o_ice))
     1502!                  call writediagfi(ngrid,"pre_moistadj_pq_h2o_vap"," "," ",3,pq(:,:,igcm_h2o_vap))
     1503!                  call writediagfi(ngrid,"pre_moistadj_pdq_h2o_vap"," "," ",3,pdq(:,:,igcm_h2o_vap))
     1504!!!!
    14471505                  call moistadj(ngrid,nlayer,nq,pt,pq,pdq,pplev,pplay,dtmoist,dqmoist,ptimestep,rneb_man)
    14481506               
    … …  
    15361594               dqsurf(1:ngrid,igcm_h2o_vap) = dqsurf(1:ngrid,igcm_h2o_vap)+zdqsrain(1:ngrid)
    15371595               dqsurf(1:ngrid,igcm_h2o_ice) = dqsurf(1:ngrid,igcm_h2o_ice)+zdqssnow(1:ngrid)
    1538 
     1596               
     1597!!               call writediagfi(ngrid,"rain_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1598!!               call writediagfi(ngrid,"rain_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1599               
    15391600               ! Test energy conservation.
    15401601               if(enertest)then
    … …  
    16961757            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq)+zdqsrain_generic(1:ngrid,1:nq)
    16971758
     1759!!            call writediagfi(ngrid,"rain_generic_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1760
    16981761            ! Test energy conservation.
    16991762            if(enertest)then
    … …  
    17791842            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + zdqssed(1:ngrid,1:nq)
    17801843
     1844!!            call writediagfi(ngrid,"callsedim_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1845
    17811846            ! Test water conservation
    17821847            if(watertest)then
    … …  
    17921857         end if ! end of 'sedimentation'
    17931858
     1859!!         call writediagfi(ngrid,"mass_redist_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1860!!         call writediagfi(ngrid,"mass_redist_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    17941861
    17951862  ! ---------------
    … …  
    18291896         endif
    18301897
     1898!         call writediagfi(ngrid,"mass_redistribution_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1899
     1900!!         call writediagfi(ngrid,"slab_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1901!!         call writediagfi(ngrid,"slab_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
     1902
    18311903  ! ------------------
    18321904  !   VI.5. Slab Ocean
    … …  
    18391911            enddo
    18401912
    1841             call ocean_slab_ice(ptimestep,                          &
    1842                                 ngrid, knindex, tsea_ice, fluxrad,  &
    1843                                 zdqssnow, qsurf(:,igcm_h2o_ice),    &
    1844                               - zdqsdif(:,igcm_h2o_vap),            &
    1845                                 flux_sens_lat,tsea_ice, pctsrf_sic, &
    1846                                 taux,tauy,icount)
    1847 
    1848 
    1849             call ocean_slab_get_vars(ngrid,tslab,      &
    1850                                      sea_ice, flux_o,  &
    1851                                      flux_g, dt_hdiff, &
    1852                                      dt_ekman)
    1853    
     1913!!            call ocean_slab_ice(ptimestep,                          &
     1914!!                                ngrid, knindex, tsea_ice, fluxrad,  &
     1915!!                               zdqssnow, qsurf(:,igcm_h2o_ice),    &
     1916!!                              - zdqsdif(:,igcm_h2o_vap),            &
     1917!!                                flux_sens_lat,tsea_ice, pctsrf_sic, &
     1918!!                               taux,tauy,icount)
     1919
     1920
     1921!!            call ocean_slab_get_vars(ngrid,tslab,      &
     1922!!                                     sea_ice, flux_o,  &
     1923!!                                     flux_g, dt_hdiff, &
     1924!!                                     dt_ekman)
     1925           
     1926            call ocean_slab_noice(icount, ptimestep, knon, knindex, &
     1927                      zdqssnow, tsea_ice, &
     1928                      fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
     1929                      tsea_ice, taux, tauy, zmasq)
     1930
     1931            call ocean_slab_ice(icount, ptimestep, knon, knindex, &
     1932                      zdqssnow, tsea_ice, &
     1933                      fluxrad, qsurf(:,igcm_h2o_ice), flux_sens_lat, &
     1934                      tsea_ice, -zdqsdif(:,igcm_h2o_vap), taux, tauy, zmasq)
     1935
     1936            call ocean_slab_frac(pctsrf_sic, zmasq)
     1937
     1938            call ocean_slab_get_vars(ngrid, tslab, sea_ice, flux_g, &
     1939                      dt_hdiff, dt_ekman)
     1940 
    18541941            do ig=1,ngrid
    18551942               if (nint(rnat(ig)).eq.1)then
    … …  
    18571944                  tslab(ig,2) = 0.
    18581945                  tsea_ice(ig) = 0.
     1946!                  tice(ig) = 0.
    18591947                  sea_ice(ig) = 0.
    18601948                  pctsrf_sic(ig) = 0.
    … …  
    18641952
    18651953         endif ! end of 'ok_slab_ocean'
     1954
     1955!!         call writediagfi(ngrid,"hydrol_pre_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1956!!         call writediagfi(ngrid,"hydrol_pre_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    18661957
    18671958  ! -----------------------------
    … …  
    18761967                        mu0,zdtsurf,zdtsurf_hyd,hice,pctsrf_sic,            &
    18771968                        sea_ice)
     1969               
     1970!!            call writediagfi(ngrid,"hydrol_post0_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    18781971
    18791972            zdtsurf(1:ngrid)     = zdtsurf(1:ngrid) + zdtsurf_hyd(1:ngrid)
    18801973            dqsurf(1:ngrid,1:nq) = dqsurf(1:ngrid,1:nq) + dqs_hyd(1:ngrid,1:nq)
     1974
    18811975           
    18821976            qsurf(1:ngrid,1:nq) = qsurf(1:ngrid,1:nq) + ptimestep*dqsurf(1:ngrid,1:nq)
     1977
     1978!!            call writediagfi(ngrid,"hydrol_post_dqsurf"," "," ",2,dqsurf(1:ngrid,igcm_h2o_vap))
     1979!!            call writediagfi(ngrid,"hydrol_post_qsurf"," "," ",2,qsurf(1:ngrid,igcm_h2o_vap))
    18831980
    18841981            ! Test energy conservation
    … …  
    19512048               capcal(ig)=capcalocean
    19522049               fluxgrd(ig)=0.
    1953                fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
     2050!!               fluxgrdocean(ig)=pctsrf_sic(ig)*flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))
     2051               ! Dividing by cell area to have flux in W/m2
     2052               fluxgrdocean(ig)=flux_g(ig)+(1-pctsrf_sic(ig))*(dt_hdiff(ig,1)+dt_ekman(ig,1))/cell_area(ig)
    19542053               do iq=1,nsoilmx
    19552054                  tsoil(ig,iq)=tsurf(ig)
    … …  
    23452444            call wstats(ngrid,"tslab2","tslab2","K",2,tslab(:,2))
    23462445            call wstats(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
    2347             call wstats(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
     2446            call wstats(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
     2447!            call wstats(ngrid,"tice","ice temp if snow on top","K",2,tice)
    23482448            call wstats(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
    23492449            call wstats(ngrid,"rnat","nature of the surface","",2,rnat)
    … …  
    23872487      if(ok_slab_ocean) then
    23882488         call writediagfi(ngrid,"pctsrf_sic","pct ice/sea","",2,pctsrf_sic)
    2389          call writediagfi(ngrid,"tsea_ice","tsea_ice","K",2,tsea_ice)
     2489         call writediagfi(ngrid,"tsea_ice","top layer temp, snow/ice","K",2,tsea_ice)
     2490!         call writediagfi(ngrid,"tice","ice temp if snow on top","K",2,tice)
    23902491         call writediagfi(ngrid,"sea_ice","sea ice","kg/m2",2,sea_ice)
    23912492         call writediagfi(ngrid,"tslab1","tslab1","K",2,tslab(:,1))
Note: See TracChangeset for help on using the changeset viewer.