New URL for NEMO forge!   http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.
Changeset 8752 for branches – NEMO

Changeset 8752 for branches


Ignore:
Timestamp:
2017-11-20T13:54:32+01:00 (6 years ago)
Author:
dancopsey
Message:

Merged in main ICEMODEL branch (branches/2017/dev_r8183_ICEMODEL) from revision 8587 to 8726.

Location:
branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM
Files:
3 deleted
29 edited
3 copied

Legend:

Unmodified
Added
Removed

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/ORCA2_LIM3_PISCES/EXP00/namelist_ice_cfg

    r8738 r8752  
    1313!!             11 - Ice growth in open water           (namthd_do) 
    1414!!             12 - Ice salinity                       (namthd_sal) 
    15 !!             13 - Ice melt ponds                     (nammp) 
     15!!             13 - Ice melt ponds                     (namthd_pnd) 
    1616!!             14 - Ice initialization                 (namini) 
    1717!!             15 - Ice/snow albedos                   (namalb) 
    … …  
    2020! 
    2121!------------------------------------------------------------------------------ 
    22 &nampar     !   Generic parameters 
     22&nampar         !   Generic parameters 
    2323!------------------------------------------------------------------------------ 
    2424/ 
    2525!------------------------------------------------------------------------------ 
    26 &namitd     !   Ice discretization 
     26&namitd         !   Ice discretization 
    2727!------------------------------------------------------------------------------ 
    2828/ 
    2929!------------------------------------------------------------------------------ 
    30 &namdyn     !   Ice dynamics 
     30&namdyn         !   Ice dynamics 
    3131!------------------------------------------------------------------------------ 
    3232/ 
    … …  
    4848/ 
    4949!------------------------------------------------------------------------------ 
    50 &namthd     !   Ice thermodynamics 
     50&namthd         !   Ice thermodynamics 
    5151!------------------------------------------------------------------------------ 
    5252/ 
    … …  
    5656/ 
    5757!------------------------------------------------------------------------------ 
    58 &namthd_da     !   Ice lateral melting 
     58&namthd_da      !   Ice lateral melting 
    5959!------------------------------------------------------------------------------ 
    6060/ 
    6161!------------------------------------------------------------------------------ 
    62 &namthd_do     !   Ice growth in open water 
     62&namthd_do      !   Ice growth in open water 
    6363!------------------------------------------------------------------------------ 
    6464/ 
    … …  
    6868/ 
    6969!------------------------------------------------------------------------------ 
    70 &nammp      !   Melt ponds 
     70&namthd_pnd     !   Melt ponds 
    7171!------------------------------------------------------------------------------ 
    7272/ 
    7373!------------------------------------------------------------------------------ 
    74 &namini     !   Ice initialization 
     74&namini         !   Ice initialization 
    7575!------------------------------------------------------------------------------ 
    7676/ 
    7777!------------------------------------------------------------------------------ 
    78 &namalb     !   albedo parameters 
     78&namalb         !   albedo parameters 
    7979!------------------------------------------------------------------------------ 
    8080/ 
    8181!------------------------------------------------------------------------------ 
    82 &namdia     !   Diagnostics 
     82&namdia         !   Diagnostics 
    8383!------------------------------------------------------------------------------ 
    8484/ 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/ORCA2_SAS_LIM3/EXP00/namelist_ice_cfg

    r8738 r8752  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    22!! ESIM configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref 
    3 !!              1 - Generic parameters                 (namice_run) 
    4 !!              2 - Ice thickness discretization       (namice_itd) 
    5 !!              3 - Ice dynamics                       (namice_dyn) 
    6 !!              4 - Ice ridging/rafting                (namice_rdgrft) 
    7 !!              5 - Ice rheology                       (namice_rhg) 
    8 !!              6 - Ice advection                      (namice_adv) 
    9 !!              7 - Ice thermodynamics                 (namice_thd) 
    10 !!              8 - Ice salinity                       (namice_sal) 
    11 !!              9 - Ice melt ponds                     (namice_mp) 
    12 !!             10 - Ice initialization                 (namice_ini) 
    13 !!             11 - Ice/snow albedos                   (namice_alb) 
    14 !!             12 - Ice diagnostics                    (namice_dia) 
     3!!              1 - Generic parameters                 (nampar) 
     4!!              2 - Ice thickness discretization       (namitd) 
     5!!              3 - Ice dynamics                       (namdyn) 
     6!!              4 - Ice ridging/rafting                (namdyn_rdgrft) 
     7!!              5 - Ice rheology                       (namdyn_rhg) 
     8!!              6 - Ice advection                      (namdyn_adv) 
     9!!              7 - Ice surface forcing                (namforcing) 
     10!!              8 - Ice thermodynamics                 (namthd) 
     11!!              9 - Ice heat diffusion                 (namthd_zdf) 
     12!!             10 - Ice lateral melting                (namthd_da) 
     13!!             11 - Ice growth in open water           (namthd_do) 
     14!!             12 - Ice salinity                       (namthd_sal) 
     15!!             13 - Ice melt ponds                     (namthd_pnd) 
     16!!             14 - Ice initialization                 (namini) 
     17!!             15 - Ice/snow albedos                   (namalb) 
     18!!             16 - Ice diagnostics                    (namdia) 
    1519!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1620! 
    1721!------------------------------------------------------------------------------ 
    18 &namice_run     !   Generic parameters 
     22&nampar         !   Generic parameters 
    1923!------------------------------------------------------------------------------ 
    2024/ 
    2125!------------------------------------------------------------------------------ 
    22 &namice_itd     !   Ice discretization 
     26&namitd         !   Ice discretization 
    2327!------------------------------------------------------------------------------ 
    2428/ 
    2529!------------------------------------------------------------------------------ 
    26 &namice_dyn     !   Ice dynamics 
     30&namdyn         !   Ice dynamics 
    2731!------------------------------------------------------------------------------ 
    2832/ 
    2933!------------------------------------------------------------------------------ 
    30 &namice_rdgrft  !   Ice ridging/rafting 
     34&namdyn_rdgrft  !   Ice ridging/rafting 
    3135!------------------------------------------------------------------------------ 
    3236/ 
    3337!------------------------------------------------------------------------------ 
    34 &namice_rhg     !   Ice rheology 
     38&namdyn_rhg     !   Ice rheology 
    3539!------------------------------------------------------------------------------ 
    3640/ 
    3741!------------------------------------------------------------------------------ 
    38 &namice_adv     !   Ice advection 
     42&namdyn_adv     !   Ice advection 
    3943!------------------------------------------------------------------------------ 
    4044/ 
    4145!------------------------------------------------------------------------------ 
    42 &namice_thd     !   Ice thermodynamics 
     46&namforcing     !   Ice surface forcing 
    4347!------------------------------------------------------------------------------ 
    4448/ 
    4549!------------------------------------------------------------------------------ 
    46 &namice_sal     !   Ice salinity 
     50&namthd         !   Ice thermodynamics 
    4751!------------------------------------------------------------------------------ 
    4852/ 
    4953!------------------------------------------------------------------------------ 
    50 &namicemp      !   Melt ponds 
     54&namthd_zdf     !   Ice heat diffusion 
    5155!------------------------------------------------------------------------------ 
    5256/ 
    5357!------------------------------------------------------------------------------ 
    54 &namice_ini     !   Ice initialization 
     58&namthd_da      !   Ice lateral melting 
    5559!------------------------------------------------------------------------------ 
    5660/ 
    5761!------------------------------------------------------------------------------ 
    58 &namice_alb     !   albedo parameters 
     62&namthd_do      !   Ice growth in open water 
    5963!------------------------------------------------------------------------------ 
    6064/ 
    6165!------------------------------------------------------------------------------ 
    62 &namice_dia     !   Diagnostics 
     66&namthd_sal     !   Ice salinity 
    6367!------------------------------------------------------------------------------ 
    6468/ 
     69!------------------------------------------------------------------------------ 
     70&namthd_pnd     !   Melt ponds 
     71!------------------------------------------------------------------------------ 
     72/ 
     73!------------------------------------------------------------------------------ 
     74&namini         !   Ice initialization 
     75!------------------------------------------------------------------------------ 
     76/ 
     77!------------------------------------------------------------------------------ 
     78&namalb         !   albedo parameters 
     79!------------------------------------------------------------------------------ 
     80/ 
     81!------------------------------------------------------------------------------ 
     82&namdia         !   Diagnostics 
     83!------------------------------------------------------------------------------ 
     84/ 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/SHARED/field_def_nemo-lim.xml

    r8738 r8752  
    132132         <field id="icevmp"       long_name="melt pond volume"                                             standard_name="sea_ice_meltpond_volume"                            unit="m"            />  
    133133 
    134          <field id="iceconc_cat"  long_name="Sea-ice area fractions in thickness categories"               unit=""   grid_ref="grid_T_3D_ncatice" /> 
     134         <field id="iceconc_cat"  long_name="Sea-ice concentration in thickness categories"                unit=""   grid_ref="grid_T_3D_ncatice" /> 
    135135         <field id="icethic_cat"  long_name="Sea-ice thickness in thickness categories"                    unit="m"  grid_ref="grid_T_3D_ncatice" /> 
    136136         <field id="snowthic_cat" long_name="Snow thickness in thickness categories"                       unit="m"  grid_ref="grid_T_3D_ncatice" /> 
    … …  
    140140         <field id="icetemp_cat"  long_name="Ice temperature for categories"                               unit="degC"   grid_ref="grid_T_3D_ncatice" /> 
    141141         <field id="snwtemp_cat"  long_name="Snow temperature for categories"                              unit="degC"   grid_ref="grid_T_3D_ncatice" /> 
    142          <field id="iceamp_cat"   long_name="Ice melt pond fraction for categories"                        unit="%"      grid_ref="grid_T_3D_ncatice" />  
     142         <field id="iceamp_cat"   long_name="Ice melt pond concentration for categories"                   unit="%"      grid_ref="grid_T_3D_ncatice" />  
    143143         <field id="icevmp_cat"   long_name="Ice melt pond volume for categories"                          unit="m"      grid_ref="grid_T_3D_ncatice" />  
     144         <field id="icehmp_cat"   long_name="Ice melt pond thickness for categories"                       unit="m"      grid_ref="grid_T_3D_ncatice" />  
     145         <field id="iceafp_cat"   long_name="Ice melt pond fraction for categories"                        unit="m"      grid_ref="grid_T_3D_ncatice" />  
    144146 
    145147         <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     /> 
    … …  
    186188         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="kg/m2/s"   /> 
    187189         <field id="vfxthin"      long_name="thermo ice prod. for thin ice(20cm) + open water"             unit="kg/m2/s"   /> 
     190         <field id="vfxpnd"       long_name="melt pond water flux to ocean"                                unit="kg/m2/s"   /> 
    188191 
    189192         <field id="afxtot"       long_name="area tendency (total)"                                        unit="s-1"   /> 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r8738 r8752  
    1313!!             11 - Ice growth in open water           (namthd_do) 
    1414!!             12 - Ice salinity                       (namthd_sal) 
    15 !!             13 - Ice melt ponds                     (nammp) 
     15!!             13 - Ice melt ponds                     (namthd_pnd) 
    1616!!             14 - Ice initialization                 (namini) 
    1717!!             15 - Ice/snow albedos                   (namalb) 
    … …  
    2020! 
    2121!------------------------------------------------------------------------------ 
    22 &nampar     !   Generic parameters 
     22&nampar         !   Generic parameters 
    2323!------------------------------------------------------------------------------ 
    2424   jpl              =   5             !  number of ice  categories 
    … …  
    3939/ 
    4040!------------------------------------------------------------------------------ 
    41 &namitd     !   Ice discretization 
    42 !------------------------------------------------------------------------------ 
    43    rn_himean        =   2.0           !  expected domain-average ice thickness (m) 
     41&namitd         !   Ice discretization 
     42!------------------------------------------------------------------------------ 
     43   ln_cat_hfn       = .true.          !  ice categories are defined by a function following rn_himean**(-0.05) 
     44      rn_himean     =   2.0           !  expected domain-average ice thickness (m) 
     45   ln_cat_usr       = .false.         !  ice categories are defined by rn_catbnd below (m) 
     46      rn_catbnd     =   0.,0.45,1.1,2.1,3.7,6.0   
    4447   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping 
    4548/ 
    4649!------------------------------------------------------------------------------ 
    47 &namdyn     !   Ice dynamics 
     50&namdyn         !   Ice dynamics 
    4851!------------------------------------------------------------------------------ 
    4952   ln_dynFULL       = .true.          !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
    … …  
    9093!------------------------------------------------------------------------------ 
    9194   ln_rhg_EVP       = .true.          !  EVP rheology 
    92       rn_creepl     =   1.0e-12       !     creep limit [1/s] 
     95      ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
     96      rn_creepl     =   2.0e-9        !     creep limit [1/s] 
    9397      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve           
    9498      nn_nevp       = 120             !     number of EVP subcycles                              
    … …  
    118122                                      !     = 2  Redistribute a single flux over categories 
    119123                                      !          ==> coupled mode only 
    120 / 
    121 !------------------------------------------------------------------------------ 
    122 &namthd     !   Ice thermodynamics 
     124   nice_jules       = 1               !  Jules coupling (0=OFF, 1=EMULATED, 2=ACTIVE) 
     125/ 
     126!------------------------------------------------------------------------------ 
     127&namthd         !   Ice thermodynamics 
    123128!------------------------------------------------------------------------------ 
    124129   ln_icedH         = .true.          !  activate ice thickness change from growing/melting (T) or not (F) 
    125130   ln_icedA         = .true.          !  activate lateral melting param. (T) or not (F) 
    126131   ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F) 
    127    ln_icedS         = .true.          !  activate gravity drainage and flushing (T) or not (F) 
     132   ln_icedS         = .true.          !  activate brine drainage (T) or not (F) 
    128133/ 
    129134!------------------------------------------------------------------------------ 
    … …  
    136141                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    137142   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m] 
    138    ln_dqns_i        = .true.          !  change the surface non-solar flux with surface temperature (T) or not (F) 
    139143/ 
    140144!------------------------------------------------------------------------------ 
    … …  
    150154/ 
    151155!------------------------------------------------------------------------------ 
    152 &namthd_do     !   Ice growth in open water 
     156&namthd_do      !   Ice growth in open water 
    153157!------------------------------------------------------------------------------ 
    154158   rn_hinew         =   0.1           !  thickness for new ice formation in open water (m), must be larger than rn_hnewice 
    … …  
    174178/ 
    175179!------------------------------------------------------------------------------ 
    176 &nammp      !   Melt ponds 
    177 !------------------------------------------------------------------------------ 
    178    ln_pnd           = .false.         !  active melt ponds 
    179    ln_pnd_rad       = .false.         !  active melt ponds radiative coupling 
    180    ln_pnd_fw        = .false.         !  active melt ponds freshwater coupling 
    181    nn_pnd_scheme    =   0             !  type of melt pond scheme  : =0 prescribed ( Tsu=0 ), =1 empirical, =2 topographic 
    182    rn_apnd          =   0.2           !  prescribed pond fraction, at Tsu=0  : (0<rn_apnd<1, nn_pnd_scheme = 0) 
    183    rn_hpnd          =   0.05          !  prescribed pond depth, at Tsu=0     : (0<rn_apnd<1, nn_pnd_scheme = 0) 
    184 / 
    185 !------------------------------------------------------------------------------ 
    186 &namini     !   Ice initialization 
     180&namthd_pnd     !   Melt ponds 
     181!------------------------------------------------------------------------------ 
     182   ln_pnd_H12       = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
     183      ln_pnd_fwb    = .false.         !     melt ponds store freshwater or not 
     184   ln_pnd_CST       = .false.         !  activate constant melt ponds 
     185      rn_apnd       =   0.2           !     prescribed pond fraction, at Tsu=0 
     186      rn_hpnd       =   0.05          !     prescribed pond depth, at Tsu=0 
     187   ln_pnd_alb       = .false.         !  melt ponds affect albedo or not 
     188/ 
     189!------------------------------------------------------------------------------ 
     190&namini         !   Ice initialization 
    187191!------------------------------------------------------------------------------ 
    188192   ln_iceini        = .true.          !  activate ice initialization (T) or not (F) 
    … …  
    209213/ 
    210214!------------------------------------------------------------------------------ 
    211 &namalb     !   albedo parameters 
    212 !------------------------------------------------------------------------------ 
    213    nn_ice_alb       =    1            !  parameterization of ice/snow albedo 
    214                                       !     0: Shine & Henderson-Sellers (JGR 1985), giving clear-sky albedo 
    215                                       !     1: "home made" based on Brandt et al. (JClim 2005) 
    216                                       !                         and Grenfell & Perovich (JGR 2004), giving cloud-sky albedo 
    217                                       !     2: as 1 with melt ponds 
    218    rn_alb_sdry      =   0.85          !  dry snow albedo         : 0.80 (nn_ice_alb = 0); 0.85 (nn_ice_alb = 1); obs 0.85-0.87 (cloud-sky) 
    219    rn_alb_smlt      =   0.75          !  melting snow albedo     : 0.65 ( '' )          ; 0.75 ( '' )          ; obs 0.72-0.82 ( '' ) 
    220    rn_alb_idry      =   0.60          !  dry ice albedo          : 0.72 ( '' )          ; 0.60 ( '' )          ; obs 0.54-0.65 ( '' ) 
    221    rn_alb_imlt      =   0.50          !  bare puddled ice albedo : 0.53 ( '' )          ; 0.50 ( '' )          ; obs 0.49-0.58 ( '' ) 
    222    rn_alb_dpnd      =   0.27          !  ponded ice albedo       : 0.25 ( '' )          ; 0.27 ( '' )          ; obs 0.10-0.30 ( '' ) 
    223 / 
    224 !------------------------------------------------------------------------------ 
    225 &namdia     !   Diagnostics 
     215&namalb         !   albedo parameters 
     216!------------------------------------------------------------------------------ 
     217   !                                  !                          !  obs range (cloud-sky) 
     218   rn_alb_sdry      =   0.85          !  dry snow albedo         :  0.85 -- 0.87 
     219   rn_alb_smlt      =   0.75          !  melting snow albedo     :  0.72 -- 0.82 
     220   rn_alb_idry      =   0.60          !  dry ice albedo          :  0.54 -- 0.65 
     221   rn_alb_imlt      =   0.50          !  bare puddled ice albedo :  0.49 -- 0.58 
     222   rn_alb_dpnd      =   0.27          !  ponded ice albedo       :  0.10 -- 0.30  
     223/ 
     224!------------------------------------------------------------------------------ 
     225&namdia         !   Diagnostics 
    226226!------------------------------------------------------------------------------ 
    227227   ln_icediachk     = .false.         !  check online the heat, mass & salt budgets (T) or not (F) 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/CONFIG/TEST_CASES/SAS_BIPER/EXP00/namelist_ice_cfg

    r8738 r8752  
    11!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    2 !! LIM3 namelist:   
    3 !!              1 - Generic parameters                 (namicerun) 
    4 !!              2 - Diagnostics                        (namicediag) 
    5 !!              3 - Ice initialization                 (namiceini) 
    6 !!              4 - Ice discretization                 (namiceitd) 
    7 !!              5 - Ice dynamics and transport         (namicedyn) 
    8 !!              6 - Ice thermodynamics                 (namicethd) 
    9 !!              7 - Ice salinity                       (namicesal) 
    10 !!              8 - Ice mechanical redistribution      (namiceitdme) 
    11 !!              9 - Ice/snow albedos                   (namicealb) 
     2!! ESIM namelist:   
     3!!              1 - Generic parameters                 (nampar) 
     4!!              2 - Ice thickness discretization       (namitd) 
     5!!              3 - Ice dynamics                       (namdyn) 
     6!!              4 - Ice ridging/rafting                (namdyn_rdgrft) 
     7!!              5 - Ice rheology                       (namdyn_rhg) 
     8!!              6 - Ice advection                      (namdyn_adv) 
     9!!              7 - Ice surface forcing                (namforcing) 
     10!!              8 - Ice thermodynamics                 (namthd) 
     11!!              9 - Ice heat diffusion                 (namthd_zdf) 
     12!!             10 - Ice lateral melting                (namthd_da) 
     13!!             11 - Ice growth in open water           (namthd_do) 
     14!!             12 - Ice salinity                       (namthd_sal) 
     15!!             13 - Ice melt ponds                     (namthd_pnd) 
     16!!             14 - Ice initialization                 (namini) 
     17!!             15 - Ice/snow albedos                   (namalb) 
     18!!             16 - Ice diagnostics                    (namdia) 
    1219!!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
    1320! 
    1421!------------------------------------------------------------------------------ 
    15 &namicerun     !   Generic parameters 
     22&nampar         !   Generic parameters 
    1623!------------------------------------------------------------------------------ 
    17    jpl              =    1          !  number of ice  categories 
    18    nlay_i           =    1          !  number of ice  layers 
    19    ln_limthd        =  .false.       !  ice thermo   (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    20    ln_limdyn        =  .true.       !  ice dynamics (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    21    nn_limdyn        =   0           !     (ln_limdyn=T) switch for ice dynamics 
    22                                     !      2: total 
    23                                     !      1: advection only (no diffusion, no ridging/rafting) 
    24                                     !      0: advection only (as 1 but with prescribed velocity, bypass rheology) 
    25    rn_uice          =   0.5    !     (nn_limdyn=0) ice u-velocity 
    26    rn_vice          =   0.0    !     (nn_limdyn=0) ice v-velocity 
     24   jpl              =   1             !  number of ice  categories 
     25   nlay_i           =   1             !  number of ice  layers 
     26   ln_icedyn        = .true.          !  ice dynamics (T) or not (F) 
     27   ln_icethd        = .false.         !  ice thermo   (T) or not (F) 
    2728/ 
    2829!------------------------------------------------------------------------------ 
    29 &namicediag    !   Diagnostics 
     30&namitd         !   Ice discretization 
     31!------------------------------------------------------------------------------ 
     32   rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping 
     33/ 
     34!------------------------------------------------------------------------------ 
     35&namdyn         !   Ice dynamics 
     36!------------------------------------------------------------------------------ 
     37   ln_dynFULL       = .false.         !  dyn.: full ice dynamics               (rheology + advection + ridging/rafting + correction) 
     38   ln_dynRHGADV     = .false.         !  dyn.: no ridge/raft & no corrections  (rheology + advection) 
     39   ln_dynADV        = .true.          !  dyn.: only advection w prescribed vel.(rn_uvice + advection) 
     40      rn_uice       =   0.5           !        prescribed ice u-velocity 
     41      rn_vice       =   0.            !        prescribed ice v-velocity 
     42/ 
     43!------------------------------------------------------------------------------ 
     44&namdyn_rdgrft  !   Ice ridging/rafting 
    3045!------------------------------------------------------------------------------ 
    3146/ 
    3247!------------------------------------------------------------------------------ 
    33 &namiceini     !   Ice initialization 
    34 !------------------------------------------------------------------------------ 
    35                   ! -- limistate -- ! 
    36    ln_limini      = .true.         !  activate ice initialization (T) or not (F) 
    37    ln_limini_file = .true.         !  netcdf file provided for initialization (T) or not (F) 
    38    cn_dir="./" 
    39    sn_hti = 'initice'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    40    sn_hts = 'initice'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    41    sn_ati = 'initice'    , -12 ,'ati'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    42    sn_tsu = 'initice'    , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    43    sn_tmi = 'initice'    , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    44    sn_smi = 'initice'    , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    45 / 
    46 !------------------------------------------------------------------------------ 
    47 &namiceitd     !   Ice discretization 
     48&namdyn_rhg     !   Ice rheology 
    4849!------------------------------------------------------------------------------ 
    4950/ 
    5051!------------------------------------------------------------------------------ 
    51 &namicedyn     !   Ice dynamics and transport 
     52&namdyn_adv     !   Ice advection 
    5253!------------------------------------------------------------------------------ 
    53                   ! -- limtrp & limadv -- ! 
    54    nn_limadv      =    0            !  choose the advection scheme (-1=Prather ; 0=Ultimate-Macho) 
    55    nn_limadv_ord  =    5            !  choose the order of the advection scheme (if nn_limadv=0) 
    5654/ 
    5755!------------------------------------------------------------------------------ 
    58 &namicethd     !   Ice thermodynamics 
     56&namforcing     !   Ice surface forcing 
    5957!------------------------------------------------------------------------------ 
    60                   ! -- limthd_dh -- ! 
    61    ln_limdH       = .true.          !  activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    62                   ! -- limthd_da -- ! 
    63    ln_limdA       = .true.          !  activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    64                  ! -- limthd_lac -- ! 
    65    ln_limdO       = .true.          !  activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    66    rn_hnewice     = 0.02            !  thickness for new ice formation in open water (m) 
    67                   ! -- limitd_th -- ! 
    68    rn_himin       = 0.01            !  minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice 
    6958/ 
    7059!------------------------------------------------------------------------------ 
    71 &namicesal     !   Ice salinity 
     60&namthd         !   Ice thermodynamics 
    7261!------------------------------------------------------------------------------ 
    73                  ! -- limthd_sal -- ! 
    74    ln_limdS       = .true.          !  activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    7562/ 
    7663!------------------------------------------------------------------------------ 
    77 &namiceitdme   !   Ice mechanical redistribution (ridging and rafting) 
     64&namthd_zdf     !   Ice heat diffusion 
    7865!------------------------------------------------------------------------------ 
    79                   ! -- limitd_me -- ! 
    80    ln_ridging     =   .true.        !  ridging activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    81    ln_rafting     =   .true.        !  rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 
    8266/ 
    83 !----------------------------------------------------------------------- 
    84 &namicealb     !   albedo parameters 
    85 !----------------------------------------------------------------------- 
     67!------------------------------------------------------------------------------ 
     68&namthd_da      !   Ice lateral melting 
     69!------------------------------------------------------------------------------ 
    8670/ 
     71!------------------------------------------------------------------------------ 
     72&namthd_do      !   Ice growth in open water 
     73!------------------------------------------------------------------------------ 
     74/ 
     75!------------------------------------------------------------------------------ 
     76&namthd_sal     !   Ice salinity 
     77!------------------------------------------------------------------------------ 
     78/ 
     79!------------------------------------------------------------------------------ 
     80&namthd_pnd     !   Melt ponds 
     81!------------------------------------------------------------------------------ 
     82/ 
     83!------------------------------------------------------------------------------ 
     84&namini         !   Ice initialization 
     85!------------------------------------------------------------------------------ 
     86/ 
     87!------------------------------------------------------------------------------ 
     88&namalb         !   albedo parameters 
     89!------------------------------------------------------------------------------ 
     90/ 
     91!------------------------------------------------------------------------------ 
     92&namdia         !   Diagnostics 
     93!------------------------------------------------------------------------------ 
     94/ 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r8742 r8752  
    178178   ! 
    179179   !                                     !!** ice-rheology namelist (namrhg) ** 
     180   LOGICAL , PUBLIC ::   ln_aEVP          !: using adaptive EVP (T or F)  
    180181   REAL(wp), PUBLIC ::   rn_creepl        !: creep limit : has to be under 1.0e-9 
    181182   REAL(wp), PUBLIC ::   rn_ecc           !: eccentricity of the elliptical yield curve 
    … …  
    195196   !                                      !   = 2  Redistribute a single flux over categories 
    196197 
     198   INTEGER , PUBLIC            ::   nice_jules            ! Choice of jules coupling  
     199   !                                                      ! Associated indices 
     200   INTEGER , PUBLIC, PARAMETER ::   np_jules_OFF    = 0   ! no Jules coupling (ice thermodynamics forced via qsr and qns) 
     201   INTEGER , PUBLIC, PARAMETER ::   np_jules_EMULE  = 1   ! emulated Jules coupling via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 
     202   INTEGER , PUBLIC, PARAMETER ::   np_jules_ACTIVE = 2   ! active Jules coupling                      (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 
     203 
    197204   !                                     !!** ice-salinity namelist (namthd_sal) ** 
    198205   INTEGER , PUBLIC ::   nn_icesal        !: salinity configuration used in the model 
    … …  
    204211   REAL(wp), PUBLIC ::   rn_simin         !: minimum ice salinity [PSU] 
    205212 
    206    ! MV MP 2016 
    207    !                                     !!** melt pond namelist (nammp) 
    208    LOGICAL , PUBLIC ::   ln_pnd           !: activate ponds or not 
    209    LOGICAL , PUBLIC ::   ln_pnd_rad       !: ponds radiatively active or not 
    210    LOGICAL , PUBLIC ::   ln_pnd_fw        !: ponds active wrt meltwater or not 
    211    INTEGER , PUBLIC ::   nn_pnd_scheme    !: type of melt pond scheme:   =0 prescribed, =1 empirical, =2 topographic 
    212    REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1), only if nn_pnd_scheme = 0 
    213    REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1), only if nn_pnd_scheme = 0 
    214    ! END MV MP 2016 
     213   !                                     !!** namelist namthd_pnd 
     214   LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012 
     215   LOGICAL , PUBLIC ::   ln_pnd_fwb       !: melt ponds store freshwater 
     216   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
     217   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1) 
     218   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1) 
     219   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo 
     220 
    215221   !                                     !!** ice-diagnostics namelist (namdia) ** 
    216222   LOGICAL , PUBLIC ::   ln_icediachk     !: flag for ice diag (T) or not (F) 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icealb.F90

    r8738 r8752  
    3636   ! 
    3737   ! ** albedo namelist (namalb) 
    38    INTEGER  ::   nn_ice_alb       ! type of albedo scheme: 0: Shine & Henderson-Sellers (JGR 1985) 
    39    !                                      !                         1: "home made" based on Brandt et al. (JClim 2005) 
    40    !                                      !                                             and Grenfell & Perovich (JGR 2004) 
    41    !                                      !                         2: Same as 1 but with melt ponds 
    4238   REAL(wp) ::   rn_alb_sdry      ! dry snow albedo 
    4339   REAL(wp) ::   rn_alb_smlt      ! melting snow albedo 
    … …  
    5349CONTAINS 
    5450 
    55    SUBROUTINE ice_alb( pt_ice, ph_ice, ph_snw, pafrac_pnd, ph_pnd, ld_pnd, pa_ice_cs, pa_ice_os ) 
     51   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os ) 
    5652      !!---------------------------------------------------------------------- 
    5753      !!               ***  ROUTINE ice_alb  *** 
    … …  
    5955      !! ** Purpose :   Computation of the albedo of the snow/ice system  
    6056      !!        
    61       !! ** Method  :   Two schemes are available (from namelist parameter nn_ice_alb) 
    62       !!                  0: the scheme is that of Shine & Henderson-Sellers (JGR 1985) for clear-skies 
    63       !!                  1: the scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005) 
    64       !!                                                                           and Grenfell & Perovich (JGR 2004) 
    65       !!                  2: fractional surface-based formulation of scheme 1 (NEW) 
    66       !!                Description of scheme 1: 
     57      !! ** Method  :   The scheme is "home made" (for cloudy skies) and based on Brandt et al. (J. Climate 2005) 
     58      !!                                                                      and Grenfell & Perovich (JGR 2004) 
    6759      !!                  1) Albedo dependency on ice thickness follows the findings from Brandt et al (2005) 
    6860      !!                     which are an update of Allison et al. (JGR 1993) ; Brandt et al. 1999 
    … …  
    7668      !!                  4) The needed 4 parameters are: dry and melting snow, freezing ice and bare puddled ice 
    7769      !! 
    78       !! ** Note    :   The parameterization from Shine & Henderson-Sellers presents several misconstructions: 
     70      !!                     compilation of values from literature (reference overcast sky values) 
     71      !!                        rn_alb_sdry = 0.85      ! dry snow 
     72      !!                        rn_alb_smlt = 0.75      ! melting snow 
     73      !!                        rn_alb_idry = 0.60      ! bare frozen ice 
     74      !!                        rn_alb_imlt = 0.50      ! bare puddled ice albedo 
     75      !!                        rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015) 
     76      !!                                                ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich 
     77      !!                     Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
     78      !!                        rn_alb_sdry = 0.85      ! dry snow 
     79      !!                        rn_alb_smlt = 0.72      ! melting snow 
     80      !!                        rn_alb_idry = 0.65      ! bare frozen ice 
     81      !!                     Brandt et al 2005 (East Antarctica) 
     82      !!                        rn_alb_sdry = 0.87      ! dry snow 
     83      !!                        rn_alb_smlt = 0.82      ! melting snow 
     84      !!                        rn_alb_idry = 0.54      ! bare frozen ice 
     85      !! 
     86      !! ** Note    :   The old parameterization from Shine & Henderson-Sellers (not here anymore) presented several misconstructions 
    7987      !!                  1) ice albedo when ice thick. tends to 0 is different than ocean albedo 
    8088      !!                  2) for small ice thick. covered with some snow (<3cm?), albedo is larger  
    … …  
    8795      !!                Grenfell & Perovich 2004, JGR, vol 109  
    8896      !!---------------------------------------------------------------------- 
    89       REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_ice       !  ice surface temperature (Kelvin) 
     97      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pt_su        !  ice surface temperature (Kelvin) 
    9098      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_ice       !  sea-ice thickness 
    9199      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_snw       !  snow depth 
     100      LOGICAL , INTENT(in   )                   ::   ld_pnd_alb   !  effect of melt ponds on albedo 
    92101      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
    93102      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
    94       LOGICAL , INTENT(in   )                   ::   ld_pnd       !  melt ponds radiatively active or not 
    95       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_cs    !  albedo of ice under clear    sky 
    96       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   pa_ice_os    !  albedo of ice under overcast sky 
     103      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky 
     104      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky 
    97105      ! 
    98106      INTEGER  ::   ji, jj, jl                ! dummy loop indices 
    99       REAL(wp) ::   zswitch, z1_c1, z1_c2     ! local scalar 
    100       REAL(wp) ::   z1_href_pnd               !   -      -                  
    101       REAL(wp) ::   zalb_sm, zalb_sf, zalb_st ! albedo of snow melting, freezing, total 
    102       ! 
    103       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb, zalb_it   ! intermediate variable & albedo of ice (snow free) 
    104 !! MV MP 
    105       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_pnd        ! ponded sea ice albedo 
    106       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_ice        ! bare sea ice albedo 
    107       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zalb_snw        ! snow-covered sea ice albedo 
    108       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_snw      ! relative snow fraction 
    109       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_ice      ! relative ice fraction 
    110       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zafrac_pnd      ! relative ice fraction (effective) 
     107      REAL(wp) ::   z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar 
     108      REAL(wp) ::   z1_href_pnd               ! inverse of the characteristic length scale (Lecomte et al. 2015) 
     109      REAL(wp) ::   zalb_pnd, zafrac_pnd      ! ponded sea ice albedo & relative pound fraction 
     110      REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction 
     111      REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction 
    111112      !!--------------------------------------------------------------------- 
    112113      ! 
    113114      IF( nn_timing == 1 )   CALL timing_start('icealb') 
    114115      ! 
    115       !----------------------------------------------------- 
    116       !  Snow-free albedo (no ice thickness dependence yet) 
    117       !----------------------------------------------------- 
    118       ! 
    119       ! Part common to nn_ice_alb = 0, 1, 2 
    120       ! 
    121       IF ( .NOT. ld_pnd ) THEN   !--- No melt ponds OR radiatively inactive melt ponds 
    122          ! Bare ice albedo is prescribed, with implicit assumption on pond fraction and depth 
    123          WHERE     ( ph_snw == 0._wp .AND. pt_ice >= rt0 )   ;   zalb(:,:,:) = rn_alb_imlt 
    124          ELSEWHERE                                           ;   zalb(:,:,:) = rn_alb_idry 
    125          END WHERE 
    126       ENDIF 
    127  
    128       SELECT CASE ( nn_ice_alb )       ! select a parameterization 
    129       ! 
    130       !              !------------------------------------------ 
    131       CASE( 0 )      !  Shine and Henderson-Sellers (1985)     !   (based on clear sky values) 
    132          !           !------------------------------------------ 
    133          ! 
    134          !                       ! Thickness-dependent bare ice albedo 
    135          WHERE    ( 1.5  < ph_ice                     )   ;   zalb_it = zalb 
    136          ELSEWHERE( 1.0  < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = 0.472  + 2.0 * ( zalb - 0.472 ) * ( ph_ice - 1.0 ) 
    137          ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.0 )   ;   zalb_it = 0.2467 + 0.7049 * ph_ice              & 
    138             &                                                                  - 0.8608 * ph_ice * ph_ice     & 
    139             &                                                                  + 0.3812 * ph_ice * ph_ice * ph_ice 
    140          ELSEWHERE                                        ;   zalb_it = 0.1    + 3.6    * ph_ice 
    141          END WHERE 
    142          ! 
    143          IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo 
    144             z1_href_pnd = 1. / 0.05    ! inverse of the characteristic length scale (Lecomte et al. 2015) 
    145             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd )  
    146             !                          ! Snow-free ice albedo is a function of pond fraction 
    147             WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 )    
    148                zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd  
    149             END WHERE 
    150          ENDIF  
    151          ! 
    152          DO jl = 1, jpl 
    153             DO jj = 1, jpj 
    154                DO ji = 1, jpi 
    155                   !                    ! Freezing snow 
    156                   ! no effect of underlying ice layer IF snow thickness > c1. Albedo does not depend on snow thick if > ppc2 
    157                   zswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ( ph_snw(ji,jj,jl) - ppc1 ) ) ) 
    158                   zalb_sf = ( 1._wp - zswitch ) * (  zalb_it(ji,jj,jl)  & 
    159                      &                                   + ph_snw(ji,jj,jl) * ( rn_alb_sdry - zalb_it(ji,jj,jl) ) * pp1_c1  )   & 
    160                      &    +           zswitch   * rn_alb_sdry   
    161                      ! 
    162                   !                    ! Melting snow 
    163                   ! no effect of underlying ice layer. Albedo does not depend on snow thick IF > ppc2 
    164                   zswitch = MAX( 0._wp , SIGN( 1._wp , ph_snw(ji,jj,jl) - ppc2 ) ) 
    165                   zalb_sm = ( 1._wp - zswitch ) * ( rn_alb_imlt + ph_snw(ji,jj,jl) * ( rn_alb_smlt - rn_alb_imlt ) * pp1_c2 )   & 
    166                      &      +         zswitch   *   rn_alb_smlt  
    167                      ! 
    168                   !                    ! Snow albedo 
    169                   zswitch =  MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
    170                   zalb_st =  zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    171                   ! 
    172                   !                    ! Surface albedo 
    173                   zswitch             = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
    174                   pa_ice_cs(ji,jj,jl) = zswitch * zalb_st + ( 1._wp - zswitch ) * zalb_it(ji,jj,jl) 
    175                   ! 
    176                END DO 
     116      z1_href_pnd = 0.05 
     117      z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
     118      z1_c2 = 1. / 0.05 
     119      z1_c3 = 1. / 0.02 
     120      z1_c4 = 1. / 0.03 
     121      ! 
     122      DO jl = 1, jpl 
     123         DO jj = 1, jpj 
     124            DO ji = 1, jpi 
     125               !                       !--- Specific snow, ice and pond fractions (for now, we prevent melt ponds and snow at the same time) 
     126               IF( ph_snw(ji,jj,jl) == 0._wp ) THEN 
     127                  zafrac_snw = 0._wp 
     128                  IF( ld_pnd_alb ) THEN 
     129                     zafrac_pnd = pafrac_pnd(ji,jj,jl) 
     130                  ELSE 
     131                     zafrac_pnd = 0._wp 
     132                  ENDIF 
     133                  zafrac_ice = 1._wp - zafrac_pnd 
     134               ELSE 
     135                  zafrac_snw = 1._wp      ! Snow fully "shades" melt ponds and ice 
     136                  zafrac_pnd = 0._wp 
     137                  zafrac_ice = 0._wp 
     138               ENDIF 
     139               ! 
     140               !                       !--- Bare ice albedo (for hi > 150cm) 
     141               IF( ld_pnd_alb ) THEN 
     142                  zalb_ice = rn_alb_idry 
     143               ELSE 
     144                  IF( ph_snw(ji,jj,jl) == 0._wp .AND. pt_su(ji,jj,jl) >= rt0 ) THEN  ;   zalb_ice = rn_alb_imlt 
     145                  ELSE                                                               ;   zalb_ice = rn_alb_idry   ;   ENDIF 
     146               ENDIF 
     147               !                       !--- Bare ice albedo (for hi < 150cm) 
     148               IF( 0.05 < ph_ice(ji,jj,jl) .AND. ph_ice(ji,jj,jl) <= 1.5 ) THEN      ! 5cm < hi < 150cm 
     149                  zalb_ice = zalb_ice    + ( 0.18 - zalb_ice   ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice(ji,jj,jl)) ) 
     150               ELSEIF( ph_ice(ji,jj,jl) <= 0.05 ) THEN                               ! 0cm < hi < 5cm 
     151                  zalb_ice = rn_alb_oce  + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice(ji,jj,jl) 
     152               ENDIF 
     153               ! 
     154               !                       !--- Snow-covered ice albedo (freezing, melting cases) 
     155               IF( pt_su(ji,jj,jl) < rt0_snow ) THEN 
     156                  zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c3 ) 
     157               ELSE 
     158                  zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw(ji,jj,jl) * z1_c4 ) 
     159               ENDIF 
     160               !                       !--- Ponded ice albedo 
     161               IF( ld_pnd_alb ) THEN 
     162                  zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd )  
     163               ELSE 
     164                  zalb_pnd = rn_alb_dpnd 
     165               ENDIF 
     166               !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions 
     167               palb_os(ji,jj,jl) = zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice 
     168               ! 
    177169            END DO 
    178170         END DO 
    179          ! 
    180          pa_ice_os(:,:,:) = pa_ice_cs(:,:,:) + ppcloud       ! Oberhuber correction for overcast sky 
    181          ! 
    182          ! 
    183          !           !------------------------------------------ 
    184       CASE( 1 )      !  New parameterization (2016)            !    (based on overcast sky values) 
    185          !           !------------------------------------------ 
    186          ! 
    187          !                    compilation of values from literature (reference overcast sky values) 
    188          !                          rn_alb_sdry = 0.85      ! dry snow 
    189          !                          rn_alb_smlt = 0.75      ! melting snow 
    190          !                          rn_alb_idry = 0.60      ! bare frozen ice 
    191          !                          rn_alb_dpnd = 0.36      ! ponded-ice overcast albedo (Lecomte et al, 2015) 
    192          !                                                  ! early melt pnds 0.27, late melt ponds 0.14 Grenfell & Perovich 
    193          !                    Perovich et al 2002 (Sheba) => the only dataset for which all types of ice/snow were retrieved 
    194          !                          rn_alb_sdry = 0.85      ! dry snow 
    195          !                          rn_alb_smlt = 0.72      ! melting snow 
    196          !                          rn_alb_idry = 0.65      ! bare frozen ice 
    197          !                    Brandt et al 2005 (East Antarctica) 
    198          !                          rn_alb_sdry = 0.87      ! dry snow 
    199          !                          rn_alb_smlt = 0.82      ! melting snow 
    200          !                          rn_alb_idry = 0.54      ! bare frozen ice 
    201          !  
    202          !              !--- Computation of snow-free ice albedo  
    203          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
    204          z1_c2 = 1. / 0.05 
    205  
    206          !              !--- Accounting for the ice-thickness dependency 
    207          WHERE     ( 1.5  < ph_ice                     )   ;   zalb_it = zalb 
    208          ELSE WHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_it = zalb + ( 0.18 - zalb ) * z1_c1 * ( LOG(1.5) - LOG(ph_ice) ) 
    209          ELSE WHERE                                        ;   zalb_it = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
    210          END WHERE 
    211          ! 
    212          IF ( ld_pnd ) THEN      ! Depth-dependent ponded ice albedo 
    213             z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015) 
    214             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_it ) * EXP( - ph_pnd * z1_href_pnd )  
    215             ! 
    216             !                    ! Snow-free ice albedo is weighted mean of ponded ice and bare ice contributions 
    217             WHERE ( ph_snw == 0._wp .AND. pt_ice >= rt0 ) 
    218                zalb_it = zalb_it * ( 1. - pafrac_pnd  ) + zalb_pnd * pafrac_pnd 
    219             END WHERE 
    220          ENDIF  
    221          ! 
    222          !              !--- Overcast sky surface albedo (accounting for snow, ice melt ponds) 
    223          z1_c1 = 1. / 0.02 
    224          z1_c2 = 1. / 0.03 
    225          DO jl = 1, jpl 
    226             DO jj = 1, jpj 
    227                DO ji = 1, jpi 
    228                   ! Snow depth dependence of snow albedo 
    229                   zalb_sf = rn_alb_sdry - ( rn_alb_sdry - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c1 ); 
    230                   zalb_sm = rn_alb_smlt - ( rn_alb_smlt - zalb_it(ji,jj,jl)) * EXP( - ph_snw(ji,jj,jl) * z1_c2 ); 
    231  
    232                   ! Snow albedo (MV I guess we could use rt0 instead of rt0_snow) 
    233                   zswitch = MAX( 0._wp , SIGN( 1._wp , pt_ice(ji,jj,jl) - rt0_snow ) )    
    234                   zalb_st = zswitch * zalb_sm + ( 1._wp - zswitch ) * zalb_sf 
    235  
    236                   ! Surface albedo    
    237                   zswitch             = MAX( 0._wp , SIGN( 1._wp , - ph_snw(ji,jj,jl) ) ) 
    238                   pa_ice_os(ji,jj,jl) = ( 1._wp - zswitch ) * zalb_st + zswitch *  zalb_it(ji,jj,jl) 
    239  
    240               END DO 
    241             END DO 
    242          END DO 
    243          ! 
    244          !              !--- Clear sky surface albedo 
    245          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 );  
    246          ! 
    247          ! 
    248          !           !------------------------------------------ 
    249       CASE( 2 )      !  Fractional surface-based parameterization (2017) 
    250          !           !------------------------------------------ 
    251          !              MV: I propose this formulation that is more elegant, and more easy to expand towards 
    252          !              varying pond and snow fraction. 
    253          !              Formulation 1 has issues to handle ponds and snow together that can't easily be fixed.  
    254          !              This one handles it better, I believe. 
    255          ! 
    256          !                 !==  Snow, bare ice and ponded ice fractions  ==! 
    257          ! 
    258          !                       ! Specific fractions (zafrac) refer to relative area covered by snow within each ice category 
    259          ! 
    260          !                       !--- Effective pond fraction (for now, we prevent melt ponds and snow at the same time) 
    261          zafrac_pnd = 0._wp 
    262          IF ( ld_pnd ) THEN 
    263             WHERE( ph_snw == 0._wp )   zafrac_pnd = pafrac_pnd   ! Snow fully "shades" melt ponds 
    264          ENDIF          
    265          ! 
    266          !                       !--- Specific snow fraction (for now, prescribed) 
    267          WHERE( ph_snw > 0._wp )   ;   zafrac_snw = 1. 
    268          ELSEWHERE                 ;   zafrac_snw = 0. 
    269          END WHERE 
    270          ! 
    271          !                       !--- Specific ice fraction  
    272          zafrac_ice = 1. - zafrac_snw - zafrac_pnd 
    273          ! 
    274          !                 !==  Snow-covered, pond-covered, and bare ice albedos  ==! 
    275          ! 
    276          !                       !--- Bare ice albedo 
    277          z1_c1 = 1. / ( LOG(1.5) - LOG(0.05) )  
    278          z1_c2 = 1. / 0.05 
    279          WHERE    ( 1.5  < ph_ice                     )   ;   zalb_ice = zalb 
    280          ELSEWHERE( 0.05 < ph_ice .AND. ph_ice <= 1.5 )   ;   zalb_ice = zalb       + ( 0.18 - zalb       ) * z1_c1 *  & 
    281             &                                                                       ( LOG(1.5) - LOG(ph_ice) ) 
    282          ELSEWHERE                                        ;   zalb_ice = rn_alb_oce + ( 0.18 - rn_alb_oce ) * z1_c2 * ph_ice 
    283          END WHERE 
    284          ! 
    285          z1_c1 = 1. / 0.02       !--- Snow-covered ice albedo (freezing, melting cases) 
    286          z1_c2 = 1. / 0.03 
    287          ! 
    288          WHERE( pt_ice < rt0_snow )   ;   zalb_snw = rn_alb_sdry - ( rn_alb_sdry - zalb_ice ) * EXP( - ph_snw * z1_c1 ) 
    289          ELSEWHERE                    ;   zalb_snw = rn_alb_smlt - ( rn_alb_smlt - zalb_ice ) * EXP( - ph_snw * z1_c2 ) 
    290          END WHERE 
    291          ! 
    292          IF ( ld_pnd ) THEN      !--- Depth-dependent ponded ice albedo 
    293             z1_href_pnd = 0.05        ! inverse of the characteristic length scale (Lecomte et al. 2015) 
    294             zalb_pnd  = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd * z1_href_pnd )  
    295          ELSE 
    296             zalb_pnd  = rn_alb_dpnd 
    297          ENDIF 
    298          !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions 
    299          pa_ice_os = zafrac_snw * zalb_snw  +  zafrac_pnd * zalb_pnd  +  zafrac_ice * zalb_ice 
    300          pa_ice_cs = pa_ice_os - ( - 0.1010 * pa_ice_os * pa_ice_os + 0.1933 * pa_ice_os - 0.0148 ) 
    301  
    302       END SELECT 
     171      END DO 
     172      ! 
     173      palb_cs(:,:,:) = palb_os(:,:,:) - ( - 0.1010 * palb_os(:,:,:) * palb_os(:,:,:) + 0.1933 * palb_os(:,:,:) - 0.0148 ) 
    303174      ! 
    304175      IF( nn_timing == 1 )   CALL timing_stop('icealb') 
    … …  
    317188      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    318189      !! 
    319       NAMELIST/namalb/ nn_ice_alb, rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd 
     190      NAMELIST/namalb/ rn_alb_sdry, rn_alb_smlt, rn_alb_idry, rn_alb_imlt, rn_alb_dpnd 
    320191      !!---------------------------------------------------------------------- 
    321192      ! 
    … …  
    334205         WRITE(numout,*) '~~~~~~~~~~~~' 
    335206         WRITE(numout,*) '   Namelist namalb:' 
    336          WRITE(numout,*) '      choose the albedo parameterization   nn_ice_alb  = ', nn_ice_alb 
    337207         WRITE(numout,*) '      albedo of dry snow                   rn_alb_sdry = ', rn_alb_sdry 
    338208         WRITE(numout,*) '      albedo of melting snow               rn_alb_smlt = ', rn_alb_smlt 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90

    r8738 r8752  
    7676      IF( icount == 0 ) THEN 
    7777         !                          ! water flux 
    78          pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    79             &                    wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    80             &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)    & 
     78         pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
     79            &                    wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
     80            &                    wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
     81            &                    wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    8182            &                  ) * e1e2t(:,:) ) * zconv 
    8283         ! 
    … …  
    9697 
    9798         pdiag_t = glob_sum( (  SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 )     & 
    98             &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 )   ) * e1e2t ) * zconv 
     99            &                 + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv 
    99100 
    100101      ELSEIF( icount == 1 ) THEN 
    101102 
    102103         ! water flux 
    103          zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:)     + wfx_sum(:,:)     + wfx_sni(:,:)     +                     & 
    104             &                wfx_opw(:,:) + wfx_res(:,:)     + wfx_dyn(:,:)     + wfx_lam(:,:)     + wfx_ice_sub(:,:) +  & 
    105             &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + wfx_spr(:,:)        & 
     104         zfv  = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) +                  & 
     105            &                wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:)  +  & 
     106            &                wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) +  & 
     107            &                wfx_ice_sub(:,:) + wfx_spr(:,:)  & 
    106108            &              ) * e1e2t(:,:) ) * zconv - pdiag_fv 
    107109 
    … …  
    117119  
    118120         ! outputs 
    119          zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t  ) * zconv  & 
     121         zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv  & 
    120122            &     - pdiag_v ) * r1_rdtice - zfv ) * rday 
    121123 
    122          zs = ( ( glob_sum( SUM( sv_i * rhoic            , dim=3 ) * e1e2t ) * zconv  & 
     124         zs = ( ( glob_sum( SUM( sv_i * rhoic             , dim=3 ) * e1e2t ) * zconv  & 
    123125            &     - pdiag_s ) * r1_rdtice + zfs ) * rday 
    124126 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn.F90

    r8738 r8752  
    172172            END DO 
    173173         END DO 
    174       END DO 
    175              
    176       IF ( nn_pnd_scheme > 0 ) THEN               !-- correct pond fraction to avoid a_ip > a_i 
    177          WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    178       ENDIF 
     174      END DO  
     175      !                                           !-- correct pond fraction to avoid a_ip > a_i 
     176      WHERE( a_ip(:,:,:) > a_i(:,:,:) )   a_ip(:,:,:) = a_i(:,:,:) 
    179177      ! 
    180178   END SUBROUTINE Hbig 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_pra.F90

    r8738 r8752  
    120120            z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    121121         END DO 
    122          IF ( nn_pnd_scheme > 0 ) THEN 
     122         IF ( ln_pnd_H12 ) THEN 
    123123            z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    124124            z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
    … …  
    167167                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    168168               END DO 
    169                IF ( nn_pnd_scheme > 0 ) THEN 
     169               IF ( ln_pnd_H12 ) THEN 
    170170                  CALL adv_x( zusnit, pu_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &    !--- melt pond fraction -- 
    171171                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    … …  
    220220                     &                                      syye(:,:,jk,jl), sxye(:,:,jk,jl) ) 
    221221               END DO 
    222                IF ( nn_pnd_scheme > 0 ) THEN 
     222               IF ( ln_pnd_H12 ) THEN 
    223223                  CALL adv_y( zusnit, pv_ice, 1._wp, zarea, z0ap  (:,:,jl), sxap (:,:,jl),   &   !--- melt pond fraction --- 
    224224                     &                                      sxxap (:,:,jl), syap (:,:,jl), syyap (:,:,jl), sxyap (:,:,jl)  ) 
    … …  
    249249         END DO 
    250250         ! MV MP 2016 
    251          IF ( nn_pnd_scheme > 0 ) THEN 
     251         IF ( ln_pnd_H12 ) THEN 
    252252            pa_ip  (:,:,jl)   = z0ap (:,:,jl) * r1_e1e2t(:,:) 
    253253            pv_ip  (:,:,jl)   = z0vp (:,:,jl) * r1_e1e2t(:,:) 
    … …  
    755755                  sxyage(:,:,jl)= z2d(:,:) 
    756756               END DO 
    757                IF ( nn_pnd_scheme > 0 ) THEN 
     757               IF ( ln_pnd_H12 ) THEN 
    758758                  DO jl = 1, jpl  
    759759                     WRITE(zchar,'(I2.2)') jl 
    … …  
    833833               syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp 
    834834               sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp 
    835                IF ( nn_pnd_scheme > 0 ) THEN 
     835               IF ( ln_pnd_H12 ) THEN 
    836836                  sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
    837837                  syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
    … …  
    854854            syyc0 (:,:,:) = 0._wp   ;   syye (:,:,:,:) = 0._wp   ;   syysal (:,:,:) = 0._wp   ;   syyage (:,:,:) = 0._wp 
    855855            sxyc0 (:,:,:) = 0._wp   ;   sxye (:,:,:,:) = 0._wp   ;   sxysal (:,:,:) = 0._wp   ;   sxyage (:,:,:) = 0._wp 
    856             IF ( nn_pnd_scheme > 0 ) THEN 
     856            IF ( ln_pnd_H12 ) THEN 
    857857               sxap  (:,:,:) = 0._wp    ; sxvp  (:,:,:) = 0._wp  
    858858               syap  (:,:,:) = 0._wp    ; syvp  (:,:,:) = 0._wp  
    … …  
    989989            END DO 
    990990         END DO 
    991          IF ( nn_pnd_scheme > 0 ) THEN 
     991         IF ( ln_pnd_H12 ) THEN 
    992992            DO jl = 1, jpl  
    993993               WRITE(zchar,'(I2.2)') jl 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_adv_umx.F90

    r8738 r8752  
    124124            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_s(:,:,jl) )         ! Snow volume 
    125125            CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pe_s(:,:,1,jl) )       ! Snow heat content 
    126             IF ( nn_pnd_scheme > 0 ) THEN 
     126            IF ( ln_pnd_H12 ) THEN 
    127127               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pa_ip(:,:,jl) )     ! Melt pond fraction 
    128128               CALL adv_umx( k_order, kt, zdt, zudy, zvdx, zcu_box, zcv_box, pv_ip(:,:,jl) )     ! Melt pond volume 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rdgrft.F90

    r8738 r8752  
    253253         CALL tab_3d_2d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
    254254         CALL tab_3d_2d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
    255          IF ( nn_pnd_scheme > 0 ) THEN 
    256             CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    257             CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    258          ENDIF 
     255         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
     256         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    259257         DO jl = 1, jpl 
    260258            DO jk = 1, nlay_s 
    … …  
    270268         CALL tab_2d_1d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
    271269         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    272          IF ( nn_pnd_scheme > 0 ) THEN 
    273             CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    274          ENDIF 
     270         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    275271          
    276272         !-----------------------------------------------------------------------------! 
    … …  
    320316         CALL tab_2d_3d( npti, nptidx(1:npti), sv_i_2d(1:npti,1:jpl), sv_i(:,:,:) ) 
    321317         CALL tab_2d_3d( npti, nptidx(1:npti), oa_i_2d(1:npti,1:jpl), oa_i(:,:,:) ) 
    322          IF ( nn_pnd_scheme > 0 ) THEN 
    323             CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    324             CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    325          ENDIF 
     318         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
     319         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
    326320         DO jl = 1, jpl 
    327321            DO jk = 1, nlay_s 
    … …  
    337331         CALL tab_1d_2d( npti, nptidx(1:npti), hfx_dyn_1d    (1:npti), hfx_dyn    (:,:) ) 
    338332         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_dyn_1d(1:npti), wfx_snw_dyn(:,:) ) 
    339          IF ( nn_pnd_scheme > 0 ) THEN 
    340             CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    341          ENDIF 
     333         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 
    342334 
    343335      ENDIF ! npti > 0 
    … …  
    651643 
    652644               !MV MP 2016 
    653                IF ( nn_pnd_scheme > 0 ) THEN 
     645               IF ( ln_pnd_H12 ) THEN 
    654646                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
    655647                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    … …  
    679671               ! Put the melt pond water into the ocean 
    680672               !------------------------------------------             
    681                IF ( ( nn_pnd_scheme > 0 ) .AND. ln_pnd_fw ) THEN 
     673               IF ( ln_pnd_fwb ) THEN 
    682674                  wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + ( rhofw * vprdg(ji) * ( 1._wp - rn_fpndrdg )   &        ! fresh water source for ocean 
    683675                     &                              + rhofw * vprft(ji) * ( 1._wp - rn_fpndrft ) ) * r1_rdtice 
    … …  
    702694               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1 
    703695               ! MV MP 2016 
    704                IF ( nn_pnd_scheme > 0 ) THEN 
     696               IF ( ln_pnd_H12 ) THEN 
    705697                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    706698                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
    … …  
    766758                     &                                    vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
    767759                  ! MV MP 2016 
    768                   IF ( nn_pnd_scheme > 0 ) THEN 
     760                  IF ( ln_pnd_H12 ) THEN 
    769761                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    770762                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg.F90

    r8738 r8752  
    7878      CASE( np_rhgEVP )                ! Elasto-Viscous-Plastic ! 
    7979         !                             !------------------------! 
    80          CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, u_ice, v_ice, shear_i, divu_i, delta_i ) 
     80         CALL ice_dyn_rhg_evp( kt, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 
    8181         !          
    8282      END SELECT 
    … …  
    107107      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    108108      !! 
    109       NAMELIST/namdyn_rhg/  ln_rhg_EVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 
     109      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 
    110110      !!------------------------------------------------------------------- 
    111111      ! 
    … …  
    125125         WRITE(numout,*) '   Namelist namdyn_rhg:' 
    126126         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP 
     127         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP 
    127128         WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl 
    128129         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icedyn_rhg_evp.F90

    r8738 r8752  
    5353CONTAINS 
    5454 
    55    SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, u_ice, v_ice, pshear_i, pdivu_i, pdelta_i ) 
     55   SUBROUTINE ice_dyn_rhg_evp( kt, pstress1_i, pstress2_i, pstress12_i, pshear_i, pdivu_i, pdelta_i ) 
    5656      !!------------------------------------------------------------------- 
    5757      !!                 ***  SUBROUTINE ice_dyn_rhg_evp  *** 
    … …  
    9595      !!              5) Diagnostics including charge ellipse 
    9696      !! 
    97       !! ** Notes   : There is the possibility to use mEVP from Bouillon 2013 
    98       !!              (by uncommenting some lines in part 3 and changing alpha and beta parameters) 
    99       !!              but this solution appears very unstable (see Kimmritz et al 2016) 
     97      !! ** Notes   : There is the possibility to use aEVP from the nice work of Kimmritz et al. (2016 & 2017) 
     98      !!              by setting up ln_aEVP=T (i.e. changing alpha and beta parameters). 
     99      !!              This is an upgraded version of mEVP from Bouillon et al. 2013 
     100      !!              (i.e. more stable and better convergence) 
    100101      !! 
    101102      !! References : Hunke and Dukowicz, JPO97 
    102103      !!              Bouillon et al., Ocean Modelling 2009 
    103104      !!              Bouillon et al., Ocean Modelling 2013 
     105      !!              Kimmritz et al., Ocean Modelling 2016 & 2017 
    104106      !!------------------------------------------------------------------- 
    105107      INTEGER, INTENT(in) ::   kt      ! time step 
    106108      REAL(wp), DIMENSION(:,:), INTENT(inout) ::   pstress1_i, pstress2_i, pstress12_i 
    107       REAL(wp), DIMENSION(:,:), INTENT(inout) ::   u_ice, v_ice  
    108109      REAL(wp), DIMENSION(:,:), INTENT(  out) ::   pshear_i, pdivu_i, pdelta_i  
    109110      !! 
    … …  
    114115      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
    115116      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
    116       REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013 
     117      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                       ! alpha coef from Bouillon 2009 or Kimmritz 2017 
    117118      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass 
    118119      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars 
    … …  
    126127      REAL(wp), DIMENSION(jpi,jpj) ::   z1_e1t0, z1_e2t0                ! scale factors 
    127128      REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
    128       ! 
     129      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
     130      ! 
     131      REAL(wp), DIMENSION(jpi,jpj) ::   zdt_m                           ! (dt / ice-snow_mass) on T points 
    129132      REAL(wp), DIMENSION(jpi,jpj) ::   zaU   , zaV                     ! ice fraction on U/V points 
    130       REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points 
     133      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    131134      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    132135      REAL(wp), DIMENSION(jpi,jpj) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
    … …  
    231234 
    232235      ! alpha parameters (Bouillon 2009) 
    233       zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
    234       zalph2 = zalph1 * z1_ecc2 
    235  
    236       ! alpha and beta parameters (Bouillon 2013) 
    237       !!zalph1 = 40. 
    238       !!zalph2 = 40. 
    239       !!zbeta  = 3000. 
    240       !!zbeta = REAL( nn_nevp )   ! close to classical EVP of Hunke (2001) 
    241  
    242       z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    243       z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    244  
     236      IF( .NOT. ln_aEVP ) THEN 
     237         zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     238         zalph2 = zalph1 * z1_ecc2 
     239 
     240         z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     241         z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     242      ENDIF 
     243          
    245244      ! Initialise stress tensor  
    246245      zs1 (:,:) = pstress1_i (:,:)  
    … …  
    304303            zmf(ji,jj)      = zm1 * ff_t(ji,jj) 
    305304 
     305            ! dt/m at T points (for alpha and beta coefficients) 
     306            zdt_m(ji,jj)    = zdtevp / MAX( zm1, zmmin ) 
     307             
    306308            ! m/dt 
    307309            zmU_t(ji,jj)    = zmassU * z1_dtevp 
    308310            zmV_t(ji,jj)    = zmassV * z1_dtevp 
    309  
     311             
    310312            ! Drag ice-atm. 
    311313            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
    … …  
    326328         END DO 
    327329      END DO 
    328       CALL lbc_lnk( zmf, 'T', 1. ) 
     330      CALL lbc_lnk_multi( zmf, 'T', 1., zdt_m, 'T', 1. ) 
    329331      ! 
    330332      !------------------------------------------------------------------------------! 
    … …  
    380382               ! P/delta at T points 
    381383               zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
     384 
     385               ! alpha & beta for aEVP 
     386               !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
     387               !   alpha = beta = sqrt(4*gamma) 
     388               IF( ln_aEVP ) THEN 
     389                  zalph1   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     390                  z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     391                  zalph2   = zalph1 
     392                  z1_alph2 = z1_alph1 
     393               ENDIF 
    382394                
    383395               ! stress at T points 
    … …  
    392404            DO ji = 1, jpim1 
    393405 
     406               ! alpha & beta for aEVP 
     407               IF( ln_aEVP ) THEN 
     408                  zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     409                  z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     410                  zbeta(ji,jj) = zalph2 
     411               ENDIF 
     412                
    394413               ! P/delta at F points 
    395414               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
    … …  
    411430                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
    412431                  &                  ) * r1_e1e2u(ji,jj) 
    413                   ! 
    414                   !                !--- V points 
     432               ! 
     433               !                !--- V points 
    415434               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
    416435                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
    … …  
    419438                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
    420439                  &                  ) * r1_e1e2v(ji,jj) 
    421                   ! 
    422                   !                !--- u_ice at V point 
     440               ! 
     441               !                !--- u_ice at V point 
    423442               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    424443                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
    425                   ! 
    426                   !                !--- v_ice at U point 
     444               ! 
     445               !                !--- v_ice at U point 
    427446               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    428447                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
    429                ! 
    430448            END DO 
    431449         END DO 
    432450         ! 
    433451         ! --- Computation of ice velocity --- ! 
    434          !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
     452         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta vary as in Kimmritz 2016 & 2017 
    435453         !  Bouillon et al. 2009 (eq 34-35) => stable 
    436454         IF( MOD(jter,2) == 0 ) THEN ! even iterations 
    … …  
    438456            DO jj = 2, jpjm1 
    439457               DO ji = fs_2, fs_jpim1 
    440                      !                 !--- tau_io/(v_oce - v_ice) 
     458                  !                 !--- tau_io/(v_oce - v_ice) 
    441459                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    442460                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    443                      !                 !--- Ocean-to-Ice stress 
     461                  !                 !--- Ocean-to-Ice stress 
    444462                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    445                      ! 
    446                      !                 !--- tau_bottom/v_ice 
     463                  ! 
     464                  !                 !--- tau_bottom/v_ice 
    447465                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    448466                  zTauB = - tau_icebfr(ji,jj) / zvel 
    449                      ! 
    450                      !                 !--- Coriolis at V-points (energy conserving formulation) 
     467                  ! 
     468                  !                 !--- Coriolis at V-points (energy conserving formulation) 
    451469                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    452470                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    453471                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    454                      ! 
    455                      !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     472                  ! 
     473                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    456474                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    457                      ! 
    458                      !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
     475                  ! 
     476                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    459477                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    460                      ! 
    461                      !                 !--- ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     478                  ! 
     479                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     480                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     481                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     482                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     483                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     484                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     485                     &           ) * zmaskV(ji,jj) 
     486                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    462487                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity 
    463488                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    … …  
    466491                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
    467492                     &            ) * zmaskV(ji,jj) 
    468                      ! 
    469                   ! Bouillon 2013 
    470                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    471                   !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    472                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
    473                   ! 
     493                  ENDIF 
    474494               END DO 
    475495            END DO 
    … …  
    483503            ! 
    484504            DO jj = 2, jpjm1 
    485                DO ji = fs_2, fs_jpim1 
    486                                 
    487                   ! tau_io/(u_oce - u_ice) 
     505               DO ji = fs_2, fs_jpim1           
     506                  !                 !--- tau_io/(u_oce - u_ice) 
    488507                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    489508                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    490  
    491                   ! Ocean-to-Ice stress 
     509                  !                 !--- Ocean-to-Ice stress 
    492510                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    493  
    494                   ! tau_bottom/u_ice 
     511                  ! 
     512                  !                 !--- tau_bottom/u_ice 
    495513                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    496514                  zTauB = - tau_icebfr(ji,jj) / zvel 
    497  
    498                   ! Coriolis at U-points (energy conserving formulation) 
     515                  ! 
     516                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    499517                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    500518                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    501519                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    502                    
    503                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     520                  ! 
     521                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    504522                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    505  
    506                   ! landfast switch => 0 = static friction ; 1 = sliding friction 
     523                  ! 
     524                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    507525                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    508  
    509                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     526                  ! 
     527                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     528                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     529                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     530                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     531                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
     532                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     533                     &            ) * zmaskU(ji,jj) 
     534                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    510535                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity 
    511536                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    … …  
    514539                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
    515540                     &            ) * zmaskU(ji,jj) 
    516                   ! Bouillon 2013 
    517                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    518                   !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    519                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
     541                  ENDIF 
    520542               END DO 
    521543            END DO 
    … …  
    532554            DO jj = 2, jpjm1 
    533555               DO ji = fs_2, fs_jpim1 
    534  
    535                   ! tau_io/(u_oce - u_ice) 
     556                  !                 !--- tau_io/(u_oce - u_ice) 
    536557                  zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
    537558                     &                              + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
    538  
    539                   ! Ocean-to-Ice stress 
     559                  !                 !--- Ocean-to-Ice stress 
    540560                  ztaux_oi(ji,jj) = zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
    541  
    542                   ! tau_bottom/u_ice 
     561                  ! 
     562                  !                 !--- tau_bottom/u_ice 
    543563                  zvel  = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 
    544564                  zTauB = - tau_icebfr(ji,jj) / zvel 
    545  
    546                   ! Coriolis at U-points (energy conserving formulation) 
     565                  ! 
     566                  !                 !--- Coriolis at U-points (energy conserving formulation) 
    547567                  zCorx(ji,jj)  =   0.25_wp * r1_e1u(ji,jj) *  & 
    548568                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
    549569                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
    550  
    551                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     570                  ! 
     571                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    552572                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCorx(ji,jj) + zspgU(ji,jj) + ztaux_oi(ji,jj) 
    553  
    554                   ! landfast switch => 0 = static friction ; 1 = sliding friction 
     573                  ! 
     574                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    555575                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    556  
    557                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     576                  ! 
     577                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     578                  u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     579                     &                                     + zTauE + zTauO * u_ice(ji,jj)                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     580                     &                                  ) / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     581                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )              & ! static friction => slow decrease to v=0 
     582                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                              & ! v_ice = v_oce if mass < zmmin  
     583                     &            ) * zmaskU(ji,jj) 
     584                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    558585                  u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj)  * u_ice(ji,jj)                             &  ! previous velocity 
    559586                     &                                     + zTauE + zTauO * u_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    560587                     &                                     ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )             &  ! m/dt + tau_io(only ice part) + landfast 
    561588                     &              + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )  &  ! static friction => slow decrease to v=0 
    562                      &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     589                     &              ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
    563590                     &            ) * zmaskU(ji,jj) 
    564                   ! Bouillon 2013 
    565                   !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) )                  & 
    566                   !!   &           + zfU(ji,jj) + zCorx(ji,jj) + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj)  & 
    567                   !!   &           ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 
     591                  ENDIF 
    568592               END DO 
    569593            END DO 
    … …  
    578602            DO jj = 2, jpjm1 
    579603               DO ji = fs_2, fs_jpim1 
    580           
    581                   ! tau_io/(v_oce - v_ice) 
     604                  !                 !--- tau_io/(v_oce - v_ice) 
    582605                  zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
    583606                     &                              + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
    584  
    585                   ! Ocean-to-Ice stress 
     607                  !                 !--- Ocean-to-Ice stress 
    586608                  ztauy_oi(ji,jj) = zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
    587  
    588                   ! tau_bottom/v_ice 
     609                  ! 
     610                  !                 !--- tau_bottom/v_ice 
    589611                  zvel  = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 
    590612                  ztauB = - tau_icebfr(ji,jj) / zvel 
    591                    
    592                   ! Coriolis at V-points (energy conserving formulation) 
     613                  ! 
     614                  !                 !--- Coriolis at v-points (energy conserving formulation) 
    593615                  zCory(ji,jj)  = - 0.25_wp * r1_e2v(ji,jj) *  & 
    594616                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
    595617                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
    596  
    597                   ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     618                  ! 
     619                  !                 !--- Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
    598620                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCory(ji,jj) + zspgV(ji,jj) + ztauy_oi(ji,jj) 
    599  
    600                   ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction 
     621                  ! 
     622                  !                 !--- landfast switch => 0 = static friction ; 1 = sliding friction 
    601623                  rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 
    602                    
    603                   ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     624                  ! 
     625                  IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
     626                  v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     627                     &                                  + zTauE + zTauO * v_ice(ji,jj)                                            & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     628                     &                                  ) / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     629                     &                 + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
     630                     &             ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                               & ! v_ice = v_oce if mass < zmmin 
     631                     &           ) * zmaskV(ji,jj) 
     632                  ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    604633                  v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj)  * v_ice(ji,jj)                             &  ! previous velocity 
    605634                     &                                     + zTauE + zTauO * v_ice(ji,jj)                             &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    … …  
    608637                     &              ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
    609638                     &            ) * zmaskV(ji,jj) 
    610                   ! Bouillon 2013 
    611                   !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) )                  & 
    612                   !!   &           + zfV(ji,jj) + zCory(ji,jj) + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj)  & 
    613                   !!   &           ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 
     639                  ENDIF 
    614640               END DO 
    615641            END DO 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceforcing.F90

    r8738 r8752  
    9999      !! 
    100100      !! ** Action  : It provides the following fields used in sea ice model: 
    101       !!                fr1_i0  , fr2_i0                         = 1sr & 2nd fraction of qsr penetration in ice  [%] 
    102101      !!                emp_oce , emp_ice                        = E-P over ocean and sea ice                    [Kg/m2/s] 
    103102      !!                sprecip                                  = solid precipitation                           [Kg/m2/s] 
    … …  
    130129 
    131130      ! --- cloud-sky and overcast-sky ice albedos --- ! 
    132       CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) 
     131      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) 
    133132 
    134133      ! albedo depends on cloud fraction because of non-linear spectral effects 
    … …  
    142141         ! 
    143142      CASE( jp_blk )                !--- bulk formulation 
    144                                 CALL blk_ice_flx( t_su, alb_ice )    !  
    145          IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     143                                CALL blk_ice_flx( t_su, h_s, h_i, alb_ice )    !  
     144         IF( ln_mixcpl      )   CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 
    146145         IF( nn_iceflx /= 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx ) 
    147146         ! 
    148147      CASE ( jp_purecpl )           !--- coupled formulation 
    149                                 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su ) 
     148                                CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 
    150149         IF( nn_iceflx == 2 )   CALL ice_flx_dist( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_iceflx ) 
    151150         ! 
    … …  
    268267      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    269268      !! 
    270       NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx 
     269      NAMELIST/namforcing/ rn_cio, rn_blow_s, nn_iceflx, nice_jules 
    271270      !!------------------------------------------------------------------- 
    272271      ! 
    … …  
    285284         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    286285         WRITE(numout,*) '   Namelist namforcing:' 
    287          WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio    = ', rn_cio 
    288          WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s = ', rn_blow_s 
    289          WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx = ', nn_iceflx 
     286         WRITE(numout,*) '      drag coefficient for oceanic stress              rn_cio     = ', rn_cio 
     287         WRITE(numout,*) '      coefficient for ice-lead partition of snowfall   rn_blow_s  = ', rn_blow_s 
     288         WRITE(numout,*) '      Multicategory heat flux formulation              nn_iceflx  = ', nn_iceflx 
     289         WRITE(numout,*) '      Jules coupling (0=no, 1=emulated, 2=active)      nice_jules = ', nice_jules 
    290290      ENDIF 
    291291      ! 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceistate.F90

    r8738 r8752  
    9494      REAL(wp) ::   ztmelts, zdh 
    9595      INTEGER  ::   i_hemis, i_fill, jl0 
    96       REAL(wp) ::   zarg, zV, zconv, zdv 
     96      REAL(wp) ::   zarg, zV, zconv, zdv, zfac 
    9797      INTEGER , DIMENSION(4)           ::   itest 
    9898      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d 
    … …  
    314314 
    315315         ! for constant salinity in time 
    316          IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN 
     316         IF( nn_icesal /= 2 )  THEN 
    317317            CALL ice_var_salprof 
    318318            sv_i = s_i * v_i 
    … …  
    358358         tn_ice (:,:,:) = t_su (:,:,:) 
    359359 
    360          ! MV MP 2016 
    361360         ! Melt pond volume and fraction 
    362          IF ( ln_pnd ) THEN 
    363             DO jl = 1, jpl 
    364                a_ip_frac(:,:,jl) = 0.2 *  zswitch(:,:) 
    365                h_ip     (:,:,jl) = 0.05 * zswitch(:,:) 
    366                a_ip(:,:,jl)      = a_ip_frac(:,:,jl) * a_i (:,:,jl)  
    367                v_ip(:,:,jl)      = h_ip     (:,:,jl) * a_ip(:,:,jl) 
    368             END DO 
    369          ELSE 
    370             a_ip(:,:,:)      = 0._wp 
    371             v_ip(:,:,:)      = 0._wp 
    372             a_ip_frac(:,:,:) = 0._wp 
    373             h_ip     (:,:,:) = 0._wp 
    374          ENDIF 
    375          ! END MV MP 2016 
     361         IF ( ln_pnd_CST .OR. ln_pnd_H12 ) THEN   ;   zfac = 1._wp 
     362         ELSE                                     ;   zfac = 0._wp   ;   ENDIF  
     363         DO jl = 1, jpl 
     364            a_ip_frac(:,:,jl) = rn_apnd * zswitch(:,:) * zfac 
     365            h_ip     (:,:,jl) = rn_hpnd * zswitch(:,:) * zfac 
     366         END DO 
     367         a_ip(:,:,:) = a_ip_frac(:,:,:) * a_i (:,:,:)  
     368         v_ip(:,:,:) = h_ip     (:,:,:) * a_ip(:,:,:) 
    376369 
    377370      ELSE ! if ln_iceini=false 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceitd.F90

    r8738 r8752  
    3636   PUBLIC   ice_itd_reb   ! called in icecor 
    3737 
     38   INTEGER ::              nice_catbnd   ! choice of the type of ice category function 
     39   !                                       ! associated indices: 
     40   INTEGER, PARAMETER ::   np_cathfn = 1   ! categories defined by a function 
     41   INTEGER, PARAMETER ::   np_catusr = 2   ! categories defined by the user 
     42   ! 
    3843   ! ** namelist (namitd) ** 
     44   LOGICAL  ::   ln_cat_hfn  ! ice categories are defined by function like rn_himean**(-0.05) 
    3945   REAL(wp) ::   rn_himean   ! mean thickness of the domain 
    40  
     46   LOGICAL  ::   ln_cat_usr  ! ice categories are defined by rn_catbnd 
     47   REAL(wp), DIMENSION(0:100) ::   rn_catbnd   ! ice categories bounds 
     48   ! 
    4149   !!---------------------------------------------------------------------- 
    4250   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    … …  
    293301            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    294302               a_i_1d (ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    295                ! MV MP 2016 
    296                IF ( nn_pnd_scheme > 0 ) THEN 
    297                   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    298                ENDIF 
    299                ! END MV MP 2016 
     303               IF ( ln_pnd_H12 )    a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    300304               h_i_1d(ji) = rn_himin 
    301305            ENDIF 
    … …  
    457461               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    458462               !   
    459                ! MV MP 2016  
    460                IF ( nn_pnd_scheme > 0 ) THEN 
     463               IF ( ln_pnd_H12 ) THEN 
    461464                  !                                                  ! Pond fraction 
    462465                  ztrans          = a_ip_2d(ji,jl1) * pdaice(ji,jl) !!clem: should be * zworka(ji) but it does not work 
    … …  
    468471                  v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 
    469472               ENDIF 
    470                ! END MV MP 2016 
    471473               ! 
    472474            ENDIF   ! jl1 >0 
    … …  
    643645      !! ** input   :   Namelist namitd 
    644646      !!------------------------------------------------------------------- 
    645       INTEGER  ::   jl    ! dummy loop index 
    646       INTEGER  ::   ios   ! Local integer output status for namelist read 
     647      INTEGER  ::   jl            ! dummy loop index 
     648      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read 
    647649      REAL(wp) ::   zhmax, znum, zden, zalpha   !   -      - 
    648       !! 
    649       NAMELIST/namitd/ rn_himean, rn_himin 
     650      ! 
     651      NAMELIST/namitd/ ln_cat_hfn, rn_himean, ln_cat_usr, rn_catbnd, rn_himin 
    650652      !!------------------------------------------------------------------ 
    651653      ! 
    … …  
    664666         WRITE(numout,*) '~~~~~~~~~~~~' 
    665667         WRITE(numout,*) '   Namelist namitd: ' 
    666          WRITE(numout,*) '      mean ice thickness in the domain               rn_himean = ', rn_himean 
    667          WRITE(numout,*) '      minimum ice thickness                          rn_himin  = ', rn_himin  
     668         WRITE(numout,*) '      Ice categories are defined by a function of rn_himean**(-0.05)    ln_cat_hfn = ', ln_cat_hfn 
     669         WRITE(numout,*) '         mean ice thickness in the domain                               rn_himean  = ', rn_himean 
     670         WRITE(numout,*) '      Ice categories are defined by rn_catbnd                           ln_cat_usr = ', ln_cat_usr 
     671         WRITE(numout,*) '         ice categories boundaries (m)                                  rn_catbnd  = ', rn_catbnd  
     672         WRITE(numout,*) '      minimum ice thickness                                             rn_himin   = ', rn_himin  
    668673      ENDIF 
    669674      ! 
    … …  
    671676      !  Thickness categories boundaries  ! 
    672677      !-----------------------------------! 
    673       ! 
    674       zalpha = 0.05_wp              ! max of each category (from h^(-alpha) function) 
    675       zhmax  = 3._wp * rn_himean 
    676       DO jl = 1, jpl 
    677          znum = jpl * ( zhmax+1 )**zalpha 
    678          zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp ) 
    679          hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
    680       END DO 
     678      !                             !== set the choice of ice categories ==! 
     679      ioptio = 0  
     680      IF( ln_cat_hfn ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_cathfn    ;   ENDIF 
     681      IF( ln_cat_usr ) THEN   ;   ioptio = ioptio + 1   ;   nice_catbnd = np_catusr    ;   ENDIF 
     682      IF( ioptio /= 1 )   CALL ctl_stop( 'ice_itd_init: choose one and only one ice categories boundaries' ) 
     683      ! 
     684      SELECT CASE( nice_catbnd ) 
     685      !                                !------------------------! 
     686      CASE( np_cathfn )                ! h^(-alpha) function 
     687         !                             !------------------------! 
     688         zalpha = 0.05_wp 
     689         zhmax  = 3._wp * rn_himean 
     690         DO jl = 1, jpl 
     691            znum = jpl * ( zhmax+1 )**zalpha 
     692            zden = REAL( jpl-jl , wp ) * ( zhmax + 1._wp )**zalpha + REAL( jl , wp ) 
     693            hi_max(jl) = ( znum / zden )**(1./zalpha) - 1 
     694         END DO 
     695         !                             !------------------------! 
     696      CASE( np_catusr )                ! user defined 
     697         !                             !------------------------! 
     698         DO jl = 0, jpl 
     699            hi_max(jl) = rn_catbnd(jl) 
     700         END DO 
     701         ! 
     702      END SELECT 
    681703      ! 
    682704      DO jl = 1, jpl                ! mean thickness by category 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icerst.F90

    r8738 r8752  
    147147      END DO 
    148148 
    149       ! MV MP 2016 
    150       IF ( nn_pnd_scheme > 0 ) THEN 
    151          DO jl = 1, jpl  
    152             WRITE(zchar,'(I2.2)') jl 
    153             znam = 'a_ip'//'_htc'//zchar 
    154             z2d(:,:) = a_ip(:,:,jl) 
    155             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip 
    156             znam = 'v_ip'//'_htc'//zchar 
    157             z2d(:,:) = v_ip(:,:,jl) 
    158             CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip 
    159          END DO 
    160       ENDIF 
    161       ! END MV MP 2016 
     149      DO jl = 1, jpl  
     150         WRITE(zchar,'(I2.2)') jl 
     151         znam = 'a_ip'//'_htc'//zchar 
     152         z2d(:,:) = a_ip(:,:,jl) 
     153         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! a_ip 
     154         znam = 'v_ip'//'_htc'//zchar 
     155         z2d(:,:) = v_ip(:,:,jl) 
     156         CALL iom_rstput( iter, nitrst, numriw, znam , z2d ) ! v_ip 
     157      END DO 
    162158 
    163159      DO jl = 1, jpl  
    … …  
    193189      !! ** purpose  :   read restart file 
    194190      !!---------------------------------------------------------------------- 
    195       INTEGER  :: jk, jl 
    196       REAL(wp) :: zfice, ziter 
     191      INTEGER  ::   jk, jl 
     192      INTEGER  ::   id1            ! local integer 
     193      REAL(wp) ::   zfice, ziter 
    197194      REAL(wp), DIMENSION(jpi,jpj) ::   z2d 
    198195      CHARACTER(len=25) ::   znam 
    … …  
    248245      END DO 
    249246 
    250       ! MV MP 2016 
    251       IF ( nn_pnd_scheme > 0 ) THEN 
     247      id1 = iom_varid( numrir, 'a_ip_htc01' , ldstop = .FALSE. ) 
     248      IF( id1 > 0 ) THEN                       ! fields exist (melt ponds) 
    252249         DO jl = 1, jpl  
    253250            WRITE(zchar,'(I2.2)') jl 
    … …  
    259256            v_ip(:,:,jl) = z2d(:,:) 
    260257         END DO 
    261       ENDIF 
    262       ! END MV MP 2016 
     258      ELSE                                     ! start from rest 
     259         IF(lwp) WRITE(numout,*) '   ==>>   previous run without melt ponds output then set it' 
     260         a_ip(:,:,:) = 0._wp 
     261         v_ip(:,:,:) = 0._wp 
     262      ENDIF 
    263263 
    264264      DO jl = 1, jpl  

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icestp.F90

    r8738 r8752  
    3535   USE icedyn         ! sea-ice: dynamics 
    3636   USE icethd         ! sea-ice: thermodynamics 
    37    USE limmp          ! sea-ice: melt ponds 
    3837   USE icecor         ! sea-ice: corrections 
    3938   USE iceupdate      ! sea-ice: sea surface boundary condition update 
    … …  
    154153         !------------------------------------------------------! 
    155154         ! It provides the following fields used in sea ice model: 
    156          !    fr1_i0  , fr2_i0     = 1sr & 2nd fraction of qsr penetration in ice  [%] 
    157155         !    emp_oce , emp_ice    = E-P over ocean and sea ice                    [Kg/m2/s] 
    158156         !    sprecip              = solid precipitation                           [Kg/m2/s] 
    … …  
    170168         IF( ln_icethd )                CALL ice_thd( kt )            ! -- Ice thermodynamics       
    171169         ! 
    172          IF ( ln_pnd )                  CALL lim_mp( kt )             ! -- Melt ponds 
    173170         ! 
    174171         IF( ln_icethd )                CALL ice_cor( kt , 2 )        ! -- Corrections 
    … …  
    238235      CALL ice_itd_init                ! ice thickness distribution initialization 
    239236      ! 
    240       CALL lim_mp_init                 ! set melt ponds parameters (clem: important to be located here) 
    241       ! 
     237      IF( ln_icethd ) THEN 
     238         CALL ice_thd_init             ! set ice thermodynics parameters (clem: important to call it first for melt ponds) 
     239      ENDIF    
    242240      !                                ! Initial sea-ice state 
    243241      IF( .NOT. ln_rstart ) THEN              ! start from rest: sea-ice deduced from sst 
    … …  
    255253         CALL ice_dyn_init             ! set ice dynamics parameters 
    256254      ENDIF 
    257       ! 
    258       IF( ln_icethd ) THEN 
    259          CALL ice_thd_init             ! set ice thermodynics parameters 
    260       ENDIF    
    261255      ! 
    262256      CALL ice_update_init             ! ice surface boundary condition 
    … …  
    325319         IF(lwp) WRITE(numout,*) '   nn_monocat forced to 0 as jpl>1, i.e. multi-category case is chosen' 
    326320      ENDIF 
    327       IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 
    328          CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 
    329       ENDIF 
     321!     IF ( jpl == 1 .AND. nn_monocat == 0 ) THEN 
     322!        CALL ctl_stop( 'STOP', 'par_init : if jpl=1 then nn_monocat should be between 1 and 4' ) 
     323!     ENDIF 
    330324      ! 
    331325      IF( ln_bdy .AND. ln_icediachk )   CALL ctl_warn('par_init: online conservation check does not work with BDY') 
    … …  
    398392      wfx_snw_sub(:,:) = 0._wp ; wfx_ice_sub(:,:) = 0._wp 
    399393      wfx_snw_sni(:,:) = 0._wp  
    400       ! MV MP 2016 
    401394      wfx_pnd(:,:) = 0._wp 
    402       ! END MV MP 2016 
    403395 
    404396      hfx_thd(:,:) = 0._wp   ; 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd.F90

    r8738 r8752  
    2626   USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, qns_tot, qsr_tot, sprecip, ln_cpl 
    2727   USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 
    28       &                 fr1_i0, fr2_i0 
     28      &                 qml_ice, qcn_ice, qsr_ice_tr 
    2929   USE ice1D          ! sea-ice: thermodynamics variables 
    3030   USE icethd_zdf     ! sea-ice: vertical heat diffusion 
    … …  
    3434   USE icethd_ent     ! sea-ice: enthalpy redistribution 
    3535   USE icethd_do      ! sea-ice: growth in open water 
     36   USE icethd_pnd     ! sea-ice: melt ponds 
    3637   USE iceitd         ! sea-ice: remapping thickness distribution 
    3738   USE icetab         ! sea-ice: 1D <==> 2D transformation 
    … …  
    8687      !!             - call ice_thd_rem  for remapping thickness distribution 
    8788      !!             - call ice_thd_do   for ice growth in leads 
    88       !!--------------------------------------------------------------------- 
     89      !!------------------------------------------------------------------- 
    8990      INTEGER, INTENT(in) :: kt    ! number of iteration 
    9091      ! 
    … …  
    230231            s_i_new   (1:npti) = 0._wp ; dh_s_tot (1:npti) = 0._wp  ! --- some init --- !  (important to have them here)  
    231232            dh_i_surf (1:npti) = 0._wp ; dh_i_bott(1:npti) = 0._wp 
    232             dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp 
     233            dh_snowice(1:npti) = 0._wp ; dh_i_sub (1:npti) = 0._wp ; dh_s_mlt(1:npti) = 0._wp 
    233234            ! 
    234235            IF( ln_icedH ) THEN                                     ! --- growing/melting --- ! 
    235236                              CALL ice_thd_zdf                             ! Ice/Snow Temperature profile 
    236237                              CALL ice_thd_dh                              ! Ice/Snow thickness    
     238                              CALL ice_thd_pnd                             ! Melt ponds formation 
    237239                              CALL ice_thd_ent( e_i_1d(1:npti,:) )         ! Ice enthalpy remapping 
    238240            ENDIF 
    … …  
    362364         CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             ) 
    363365         CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     ) 
    364          CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     ) 
    365          CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     ) 
     366         CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     ) 
     367         CALL tab_2d_1d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     ) 
    366368         CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     ) 
    367          CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     ) 
     369         CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     ) 
    368370         DO jk = 1, nlay_s 
    369             CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)   ) 
    370             CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)   ) 
     371            CALL tab_2d_1d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    ) 
     372            CALL tab_2d_1d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    ) 
    371373         END DO 
    372374         DO jk = 1, nlay_i 
    373             CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl)   ) 
    374             CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl)   ) 
    375             CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)   ) 
    376          END DO 
     375            CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  ) 
     376            CALL tab_2d_1d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  ) 
     377            CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
     378         END DO 
     379         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
     380         CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
     381         CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
    377382         ! 
    378383         CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d(1:npti), qprec_ice        ) 
    379384         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_1d  (1:npti), qsr_ice (:,:,kl) ) 
    380          CALL tab_2d_1d( npti, nptidx(1:npti), fr1_i0_1d   (1:npti), fr1_i0           ) 
    381          CALL tab_2d_1d( npti, nptidx(1:npti), fr2_i0_1d   (1:npti), fr2_i0           ) 
    382385         CALL tab_2d_1d( npti, nptidx(1:npti), qns_ice_1d  (1:npti), qns_ice (:,:,kl) ) 
    383386         CALL tab_2d_1d( npti, nptidx(1:npti), ftr_ice_1d  (1:npti), ftr_ice (:,:,kl) ) 
    … …  
    388391         CALL tab_2d_1d( npti, nptidx(1:npti), fhtur_1d    (1:npti), fhtur            ) 
    389392         CALL tab_2d_1d( npti, nptidx(1:npti), fhld_1d     (1:npti), fhld             ) 
     393          
     394         CALL tab_2d_1d( npti, nptidx(1:npti), qml_ice_1d   (1:npti), qml_ice      (:,:,kl)  ) 
     395         CALL tab_2d_1d( npti, nptidx(1:npti), qcn_ice_1d   (1:npti), qcn_ice      (:,:,kl) ) 
     396         CALL tab_2d_1d( npti, nptidx(1:npti), qsr_ice_tr_1d(1:npti), qsr_ice_tr   (:,:,kl) ) 
    390397         ! 
    391398         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni   ) 
    … …  
    403410         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr          ) 
    404411         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam          ) 
     412         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd          ) 
    405413         ! 
    406414         CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog          ) 
    … …  
    454462         ! 
    455463         ! Change thickness to volume (replaces routine ice_var_eqv2glo) 
    456          v_i_1d(1:npti)  = h_i_1d(1:npti) * a_i_1d(1:npti) 
    457          v_s_1d(1:npti)  = h_s_1d(1:npti) * a_i_1d(1:npti) 
    458          sv_i_1d(1:npti) = s_i_1d(1:npti) * v_i_1d(1:npti) 
     464         v_i_1d (1:npti) = h_i_1d (1:npti) * a_i_1d (1:npti) 
     465         v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 
     466         sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 
     467         v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
    459468          
    460469         CALL tab_1d_2d( npti, nptidx(1:npti), at_i_1d(1:npti), at_i             ) 
    461470         CALL tab_1d_2d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,kl)     ) 
    462          CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d(1:npti), h_i(:,:,kl)     ) 
    463          CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d(1:npti), h_s(:,:,kl)     ) 
     471         CALL tab_1d_2d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,kl)     ) 
     472         CALL tab_1d_2d( npti, nptidx(1:npti), h_s_1d (1:npti), h_s (:,:,kl)     ) 
    464473         CALL tab_1d_2d( npti, nptidx(1:npti), t_su_1d(1:npti), t_su(:,:,kl)     ) 
    465          CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,kl)     ) 
     474         CALL tab_1d_2d( npti, nptidx(1:npti), s_i_1d (1:npti), s_i (:,:,kl)     ) 
    466475         DO jk = 1, nlay_s 
    467             CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl) ) 
    468             CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl) ) 
     476            CALL tab_1d_2d( npti, nptidx(1:npti), t_s_1d(1:npti,jk), t_s(:,:,jk,kl)    ) 
     477            CALL tab_1d_2d( npti, nptidx(1:npti), e_s_1d(1:npti,jk), e_s(:,:,jk,kl)    ) 
    469478         END DO 
    470479         DO jk = 1, nlay_i 
    471             CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d(1:npti,jk), t_i(:,:,jk,kl) ) 
    472             CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d(1:npti,jk), e_i(:,:,jk,kl) ) 
    473             CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 
    474          END DO 
     480            CALL tab_1d_2d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,kl)  ) 
     481            CALL tab_1d_2d( npti, nptidx(1:npti), e_i_1d (1:npti,jk), e_i (:,:,jk,kl)  ) 
     482            CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl)  ) 
     483         END DO 
     484         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,kl) ) 
     485         CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,kl) ) 
     486         CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_frac_1d(1:npti), a_ip_frac(:,:,kl) ) 
    475487         ! 
    476488         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) 
    … …  
    488500         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr        ) 
    489501         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam        ) 
     502         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
    490503         ! 
    491504         CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog        ) 
    … …  
    523536         CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 
    524537         CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 
     538         CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) ) 
    525539         ! 
    526540      END SELECT 
    … …  
    530544 
    531545   SUBROUTINE ice_thd_init 
    532       !!----------------------------------------------------------------------- 
     546      !!------------------------------------------------------------------- 
    533547      !!                   ***  ROUTINE ice_thd_init ***  
    534548      !!                  
    … …  
    570584      IF( ln_icedO )   CALL ice_thd_do_init    ! set ice growth in open water parameters 
    571585                       CALL ice_thd_sal_init   ! set ice salinity parameters 
    572       ! 
    573       IF( ln_icedS .AND. nn_icesal == 1 ) THEN 
    574          ln_icedS = .FALSE. 
    575          CALL ctl_warn('ln_icedS is set to false since constant ice salinity is chosen (nn_icesal=1)') 
    576       ENDIF 
     586                       CALL ice_thd_pnd_init   ! set melt ponds parameters 
    577587      ! 
    578588   END SUBROUTINE ice_thd_init 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_dh.F90

    r8738 r8752  
    130130      !------------------------------------------------------------------------------! 
    131131      ! 
    132       DO ji = 1, npti 
    133          zdum       = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)  
     132      SELECT CASE ( nice_jules ) 
     133 
     134         CASE ( np_jules_ACTIVE ) 
     135 
     136            DO ji = 1, npti 
     137               zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
     138            END DO 
     139 
     140         CASE ( np_jules_EMULE ) 
     141 
     142            DO ji = 1, npti 
     143               zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) )  
     144               qml_ice_1d(ji)   =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     145               zq_su(ji)        =   MAX( 0._wp, qml_ice_1d(ji) * rdt_ice ) 
     146            END DO 
     147 
     148         CASE ( np_jules_OFF )  
     149 
     150            DO ji = 1, npti 
     151               zdum             =   ( qns_ice_1d(ji) + qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) - fc_su(ji) )  
     152               zdum             =   zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
     153               zq_su(ji)        =   MAX( 0._wp, zdum             * rdt_ice ) 
     154            END DO 
     155 
     156      END SELECT 
     157          
     158      DO ji = 1, npti 
    134159         zf_tt(ji)  = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)  
    135  
    136          zq_su (ji) = MAX( 0._wp, zdum      * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) ) 
    137160         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    138161      END DO 
    … …  
    148171            ! Contribution to mass flux 
    149172            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhosn * h_s_1d(ji) * a_i_1d(ji) * r1_rdtice 
     173            dh_s_mlt(ji)       = dh_s_mlt(ji) - h_s_1d(ji) 
    150174            ! updates 
    151             h_s_1d(ji)   = 0._wp 
     175            h_s_1d(ji)    = 0._wp 
    152176            e_s_1d (ji,1) = 0._wp 
    153177            t_s_1d (ji,1) = rt0 
    … …  
    211235         ! snow melting only = water into the ocean (then without snow precip), >0 
    212236         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     237         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    213238         ! updates available heat + precipitations after melting 
    214239         zq_su     (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,1) * zqprec(ji) )       
    … …  
    233258            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d(ji,jk) * r1_rdtice  
    234259            ! snow melting only = water into the ocean (then without snow precip) 
    235             wfx_snw_sum_1d(ji)   = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     260            wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,jk) * r1_rdtice 
     261            dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,jk) 
    236262            ! updates available heat + thickness 
    237             zq_su (ji)  = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
     263            zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 
    238264            h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 
    239265         END DO 
    … …  
    571597         zdeltah  (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 
    572598         dh_s_tot (ji)   = dh_s_tot(ji)  + zdeltah(ji,1) 
    573          h_s_1d   (ji)  = h_s_1d(ji)   + zdeltah(ji,1) 
     599         h_s_1d   (ji)   = h_s_1d(ji)   + zdeltah(ji,1) 
    574600         
    575601         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1)                ! update available heat (J.m-2) 
    … …  
    577603         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_rdtice ! W.m-2 (>0) 
    578604         ! Contribution to mass flux 
    579          wfx_snw_sum_1d(ji)  =  wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     605         wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhosn * a_i_1d(ji) * zdeltah(ji,1) * r1_rdtice 
     606         dh_s_mlt(ji)       = dh_s_mlt(ji) + zdeltah(ji,1) 
    580607         !     
    581608         ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 
    … …  
    611638 
    612639         ! Case constant salinity in time: virtual salt flux to keep salinity constant 
    613          IF( nn_icesal == 1 .OR. nn_icesal == 3 )  THEN 
     640         IF( nn_icesal /= 2 )  THEN 
    614641            sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt                  * r1_rdtice  & ! put back sss_m     into the ocean 
    615642               &                            - s_i_1d(ji)  * a_i_1d(ji) * dh_snowice(ji) * rhoic * r1_rdtice    ! and get  rn_icesal from the ocean  

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icethd_zdf.F90

    r8738 r8752  
    3434   LOGICAL  ::   ln_cndi_U64      ! thermal conductivity: Untersteiner (1964) 
    3535   LOGICAL  ::   ln_cndi_P07      ! thermal conductivity: Pringle et al (2007) 
    36    REAL(wp) ::   rn_cnd_s         ! thermal conductivity of the snow [W/m/K] 
    3736   REAL(wp) ::   rn_kappa_i       ! coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
    38    LOGICAL  ::   ln_dqns_i        ! change non-solar surface flux with changing surface temperature (T) or not (F) 
    39  
     37   REAL(wp), PUBLIC ::   rn_cnd_s ! thermal conductivity of the snow [W/m/K] 
     38 
     39   INTEGER  ::   nice_zdf         ! Choice of the type of vertical heat diffusion formulation 
     40   !                                           ! associated indices: 
     41   INTEGER, PARAMETER ::   np_BL99    = 1      ! Bitz and Lipscomb (1999) 
     42 
     43   INTEGER , PARAMETER ::   np_zdf_jules_OFF    = 0   !   compute all temperatures from qsr and qns 
     44   INTEGER , PARAMETER ::   np_zdf_jules_SND    = 1   !   compute conductive heat flux and surface temperature from qsr and qns 
     45   INTEGER , PARAMETER ::   np_zdf_jules_RCV    = 2   !   compute snow and inner ice temperatures from qcnd 
     46    
    4047   !!---------------------------------------------------------------------- 
    4148   !! NEMO/ICE 4.0 , NEMO Consortium (2017) 
    … …  
    4552CONTAINS 
    4653 
    47    SUBROUTINE ice_thd_zdf 
     54      SUBROUTINE ice_thd_zdf 
     55       
     56      !! 
    4857      !!------------------------------------------------------------------- 
    4958      !!                ***  ROUTINE ice_thd_zdf  *** 
    5059      !! ** Purpose : 
     60      !!           This chooses between the appropriate routine for the  
     61      !!           computation of vertical diffusion 
     62      !! 
     63      !!------------------------------------------------------------------- 
     64      !! 
     65    
     66      SELECT CASE ( nice_zdf )      ! Choose the vertical heat diffusion solver 
     67       
     68                                       !-------------       
     69         CASE( np_BL99 )               ! BL99 solver 
     70                                       !------------- 
     71          
     72            IF ( nice_jules == np_jules_OFF ) THEN   ! No Jules coupler      
     73 
     74               CALL ice_thd_zdf_BL99 ( np_zdf_jules_OFF ) 
     75                                        
     76            ENDIF 
     77             
     78            IF ( nice_jules == np_jules_EMULE ) THEN   ! Jules coupler is emulated 
     79             
     80               CALL ice_thd_zdf_BL99 ( np_zdf_jules_SND ) 
     81               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 
     82                                                    
     83            ENDIF 
     84 
     85            IF ( nice_jules == np_jules_ACTIVE ) THEN   ! Jules coupler is emulated 
     86             
     87               CALL ice_thd_zdf_BL99 ( np_zdf_jules_RCV ) 
     88                                                    
     89            ENDIF 
     90             
     91      END SELECT 
     92       
     93   END SUBROUTINE ice_thd_zdf 
     94    
     95    
     96    
     97   SUBROUTINE ice_thd_zdf_BL99(k_jules) 
     98      !!------------------------------------------------------------------- 
     99      !!                ***  ROUTINE ice_thd_zdf_BL99  *** 
     100      !! ** Purpose : 
    51101      !!           This routine determines the time evolution of snow and sea-ice  
    52       !!           temperature profiles. 
     102      !!           temperature profiles, using the original Bitz and Lipscomb (1999) algorithm 
     103      !! 
    53104      !! ** Method  : 
    54105      !!           This is done by solving the heat equation diffusion with 
    … …  
    83134      !!           total ice/snow thickness : h_i_1d, h_s_1d 
    84135      !!------------------------------------------------------------------- 
     136      INTEGER, INTENT(in) ::   k_jules     ! Jules coupling (0=OFF, 1=RECEIVE, 2=SEND) 
     137       
    85138      INTEGER ::   ji, jk         ! spatial loop index 
    86139      INTEGER ::   jm             ! current reference number of equation 
    … …  
    101154      REAL(wp) ::   zdti_bnd  =  1.e-4_wp     ! maximal authorized error on temperature  
    102155      REAL(wp) ::   ztmelt_i                  ! ice melting temperature 
    103       REAL(wp) ::   z1_hsu 
    104156      REAL(wp) ::   zdti_max                  ! current maximal error on temperature  
    105157      REAL(wp) ::   zcpi                      ! Ice specific heat 
    … …  
    111163      REAL(wp), DIMENSION(jpij)     ::   zh_i, z1_h_i ! ice layer thickness 
    112164      REAL(wp), DIMENSION(jpij)     ::   zh_s, z1_h_s ! snow layer thickness 
    113       REAL(wp), DIMENSION(jpij)     ::   zfsw        ! solar radiation absorbed at the surface 
    114165      REAL(wp), DIMENSION(jpij)     ::   zqns_ice_b  ! solar radiation absorbed at the surface 
    115166      REAL(wp), DIMENSION(jpij)     ::   zfnet       ! surface flux function 
    116167      REAL(wp), DIMENSION(jpij)     ::   zdqns_ice_b ! derivative of the surface flux function 
    117       REAL(wp), DIMENSION(jpij)     ::   zftrice     ! solar radiation transmitted through the ice 
    118        
     168       
     169      REAL(wp), DIMENSION(jpij       )     ::   ztsuold     ! Old surface temperature in the ice 
    119170      REAL(wp), DIMENSION(jpij,nlay_i)     ::   ztiold      ! Old temperature in the ice 
    120171      REAL(wp), DIMENSION(jpij,nlay_s)     ::   ztsold      ! Old temperature in the snow 
    … …  
    136187      REAL(wp), DIMENSION(jpij)            ::   zq_ini      ! diag errors on heat 
    137188      REAL(wp), DIMENSION(jpij)            ::   zghe        ! G(he), th. conduct enhancement factor, mono-cat 
     189 
     190      REAL(wp) ::   zfr1, zfr2, zfrqsr_tr_i   ! dummy factor 
    138191       
    139192      ! Mono-category 
    … …  
    166219      ELSEWHERE                         ;   z1_h_s(1:npti) = 0._wp 
    167220      END WHERE 
    168       ! 
    169       ! temperatures 
    170       ztsub  (1:npti)   = t_su_1d(1:npti)  ! temperature at the previous iteration 
     221 
     222      ! Store initial temperatures and non solar heat fluxes 
     223      IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
     224 
     225         ztsub  (1:npti)       =   t_su_1d(1:npti)                          ! surface temperature at iteration n-1 
     226         ztsuold(1:npti)       =   t_su_1d(1:npti)                          ! surface temperature initial value 
     227         zdqns_ice_b(1:npti)   =   dqns_ice_1d(1:npti)                      ! derivative of incoming nonsolar flux  
     228         zqns_ice_b (1:npti)   =   qns_ice_1d(1:npti)                       ! store previous qns_ice_1d value 
     229 
     230         t_su_1d(1:npti)       =   MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! required to leave the choice between melting or not 
     231 
     232      ENDIF    
     233 
    171234      ztsold (1:npti,:) = t_s_1d(1:npti,:)   ! Old snow temperature 
    172235      ztiold (1:npti,:) = t_i_1d(1:npti,:)   ! Old ice temperature 
    173       t_su_1d(1:npti)   = MIN( t_su_1d(1:npti), rt0 - ztsu_err )   ! necessary 
    174       ! 
     236 
    175237      !------------- 
    176238      ! 2) Radiation 
    177239      !------------- 
    178       z1_hsu = 1._wp / 0.1_wp ! threshold for the computation of i0 
    179       DO ji = 1, npti 
    180          ! --- Computation of i0 --- ! 
    181          ! i0 describes the fraction of solar radiation which does not contribute 
    182          ! to the surface energy budget but rather penetrates inside the ice. 
    183          ! We assume that no radiation is transmitted through the snow 
    184          ! If there is no no snow 
    185          ! zfsw    = (1-i0).qsr_ice   is absorbed at the surface  
    186          ! zftrice = io.qsr_ice       is below the surface  
    187          ! ftr_ice = io.qsr_ice.exp(-k(h_i)) transmitted below the ice  
    188          ! fr1_i0_1d = i0 for a thin ice cover, fr1_i0_2d = i0 for a thick ice cover 
    189          zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * z1_hsu ) )      
    190          i0(ji) = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zfac * fr2_i0_1d(ji) ) 
    191  
    192          ! --- Solar radiation absorbed / transmitted at the surface --- ! 
    193          !     Derivative of the non solar flux 
    194          zfsw   (ji)     =  qsr_ice_1d(ji) * ( 1 - i0(ji) )   ! Shortwave radiation absorbed at surface 
    195          zftrice(ji)     =  qsr_ice_1d(ji) *       i0(ji)     ! Solar radiation transmitted below the surface layer 
    196          zdqns_ice_b(ji) = dqns_ice_1d(ji)                    ! derivative of incoming nonsolar flux  
    197          zqns_ice_b (ji) =  qns_ice_1d(ji)                    ! store previous qns_ice_1d value 
    198       END DO 
    199  
    200240      ! --- Transmission/absorption of solar radiation in the ice --- ! 
    201       zradtr_s(1:npti,0) = zftrice(1:npti) 
     241!     zfr1 = ( 0.18 * ( 1.0 - 0.81 ) + 0.35 * 0.81 )            !   standard value 
     242!     zfr2 = ( 0.82 * ( 1.0 - 0.81 ) + 0.65 * 0.81 )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     243 
     244!     DO ji = 1, npti 
     245 
     246!        zfac = MAX( 0._wp , 1._wp - ( h_i_1d(ji) * 10._wp ) )      
     247 
     248!              zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     249!              IF ( h_s_1d(ji) >= 0.0_wp )   zfrqsr_tr_i = 0._wp     !   snow fully opaque 
     250 
     251!              qsr_ice_tr_1d(ji) = zfrqsr_tr_i * qsr_ice_1d(ji)   !   transmitted solar radiation  
     252 
     253!              zfsw(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) 
     254!              zftrice(ji) = qsr_ice_tr_1d(ji) 
     255!              i0(ji) = zfrqsr_tr_i 
     256 
     257!     END DO 
     258 
     259      zradtr_s(1:npti,0) = qsr_ice_tr_1d(1:npti) 
    202260      DO jk = 1, nlay_s 
    203261         DO ji = 1, npti 
    … …  
    209267      END DO 
    210268          
    211       zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + zftrice(1:npti) * ( 1._wp - isnow(1:npti) ) 
     269      zradtr_i(1:npti,0) = zradtr_s(1:npti,nlay_s) * isnow(1:npti) + qsr_ice_tr_1d(1:npti) * ( 1._wp - isnow(1:npti) ) 
    212270      DO jk = 1, nlay_i  
    213271         DO ji = 1, npti 
    … …  
    330388            END DO 
    331389         END DO 
    332          ! 
    333          !---------------------------- 
    334          ! 6) surface flux computation 
    335          !---------------------------- 
    336          IF ( ln_dqns_i ) THEN  
    337             DO ji = 1, npti 
    338                ! update of the non solar flux according to the update in T_su 
     390          
     391         ! 
     392         !----------------------------------------! 
     393         !                                        ! 
     394         !       JULES IF (OFF or SND MODE)       ! 
     395         !                                        ! 
     396         !----------------------------------------! 
     397         ! 
     398          
     399         IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN ! OFF or SND mode 
     400          
     401            ! 
     402            ! In OFF mode the original BL99 temperature computation is used 
     403            ! (with qsr_ice, qns_ice and dqns_ice as inputs) 
     404            ! 
     405            ! In SND mode, the computation is required to compute the conduction fluxes 
     406            ! 
     407          
     408            !---------------------------- 
     409            ! 6) surface flux computation 
     410            !---------------------------- 
     411 
     412            DO ji = 1, npti 
     413            ! update of the non solar flux according to the update in T_su 
    339414               qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_1d(ji) - ztsub(ji) ) 
    340415            END DO 
    341          ENDIF 
    342  
    343          DO ji = 1, npti 
    344             zfnet(ji) = zfsw(ji) + qns_ice_1d(ji) ! incoming = net absorbed solar radiation + non solar total flux (LWup, LWdw, SH, LH) 
    345          END DO 
    346          ! 
    347          !---------------------------- 
    348          ! 7) tridiagonal system terms 
    349          !---------------------------- 
    350          !!layer denotes the number of the layer in the snow or in the ice 
    351          !!jm denotes the reference number of the equation in the tridiagonal 
    352          !!system, terms of tridiagonal system are indexed as following : 
    353          !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
    354  
    355          !!ice interior terms (top equation has the same form as the others) 
    356          ztrid   (1:npti,:,:) = 0._wp 
    357          zindterm(1:npti,:)   = 0._wp 
    358          zindtbis(1:npti,:)   = 0._wp 
    359          zdiagbis(1:npti,:)   = 0._wp 
    360  
    361          DO jm = nlay_s + 2, nlay_s + nlay_i  
     416 
     417            DO ji = 1, npti 
     418               zfnet(ji) = qsr_ice_1d(ji) - qsr_ice_tr_1d(ji) + qns_ice_1d(ji) ! net heat flux = net solar - transmitted solar + non solar 
     419            END DO 
     420            ! 
     421            !---------------------------- 
     422            ! 7) tridiagonal system terms 
     423            !---------------------------- 
     424            !!layer denotes the number of the layer in the snow or in the ice 
     425            !!jm denotes the reference number of the equation in the tridiagonal 
     426            !!system, terms of tridiagonal system are indexed as following : 
     427            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
     428 
     429            !!ice interior terms (top equation has the same form as the others) 
     430            ztrid   (1:npti,:,:) = 0._wp 
     431            zindterm(1:npti,:)   = 0._wp 
     432            zindtbis(1:npti,:)   = 0._wp 
     433            zdiagbis(1:npti,:)   = 0._wp 
     434 
     435            DO jm = nlay_s + 2, nlay_s + nlay_i  
    362436            DO ji = 1, npti 
    363437               jk = jm - nlay_s - 1 
    … …  
    367441               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
    368442            END DO 
    369          ENDDO 
    370  
    371          jm =  nlay_s + nlay_i + 1 
    372          DO ji = 1, npti 
     443            ENDDO 
     444 
     445            jm =  nlay_s + nlay_i + 1 
     446            DO ji = 1, npti 
    373447            !!ice bottom term 
    374448            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
    … …  
    377451            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
    378452               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
    379          ENDDO 
    380  
    381  
    382          DO ji = 1, npti 
     453            ENDDO 
     454 
     455 
     456            DO ji = 1, npti 
    383457            !                               !---------------------! 
    384             IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
     458            IF ( h_s_1d(ji) > 0.0 ) THEN    !  snow-covered cells ! 
    385459               !                            !---------------------! 
    386460               ! snow interior terms (bottom equation has the same form as the others) 
    387461               DO jm = 3, nlay_s + 1 
    388                   jk = jm - 1 
    389                   ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
    390                   ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
    391                   ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
    392                   zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     462               jk = jm - 1 
     463               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     464               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
     465               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
     466               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
    393467               END DO 
    394468 
    395469               ! case of only one layer in the ice (ice equation is altered) 
    396470               IF ( nlay_i == 1 ) THEN 
    397                   ztrid(ji,nlay_s+2,3)  = 0.0 
    398                   zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     471               ztrid(ji,nlay_s+2,3)  = 0.0 
     472               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
    399473               ENDIF 
    400474 
    401475               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    402476 
    403                   jm_min(ji) = 1 
    404                   jm_max(ji) = nlay_i + nlay_s + 1 
    405  
    406                   ! surface equation 
    407                   ztrid(ji,1,1)  = 0.0 
    408                   ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
    409                   ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
    410                   zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
    411  
    412                   ! first layer of snow equation 
    413                   ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
    414                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    415                   ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1) 
    416                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
     477               jm_min(ji) = 1 
     478               jm_max(ji) = nlay_i + nlay_s + 1 
     479 
     480               ! surface equation 
     481               ztrid(ji,1,1)  = 0.0 
     482               ztrid(ji,1,2)  = zdqns_ice_b(ji) - zg1s * zkappa_s(ji,0) 
     483               ztrid(ji,1,3)  = zg1s * zkappa_s(ji,0) 
     484               zindterm(ji,1) = zdqns_ice_b(ji) * t_su_1d(ji) - zfnet(ji) 
     485 
     486               ! first layer of snow equation 
     487               ztrid(ji,2,1)  = - zkappa_s(ji,0) * zg1s * zeta_s(ji,1) 
     488               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
     489               ztrid(ji,2,3)  = - zeta_s(ji,1)* zkappa_s(ji,1) 
     490               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) * zradab_s(ji,1) 
    417491 
    418492               ELSE                            !--  case 2 : surface is melting 
    419                   ! 
    420                   jm_min(ji) = 2 
    421                   jm_max(ji) = nlay_i + nlay_s + 1 
    422  
    423                   ! first layer of snow equation 
    424                   ztrid(ji,2,1)  = 0.0 
    425                   ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
    426                   ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
    427                   zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
    428                      &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
     493               ! 
     494               jm_min(ji) = 2 
     495               jm_max(ji) = nlay_i + nlay_s + 1 
     496 
     497               ! first layer of snow equation 
     498               ztrid(ji,2,1)  = 0.0 
     499               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * ( zkappa_s(ji,1) + zkappa_s(ji,0) * zg1s ) 
     500               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
     501               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
     502               &           ( zradab_s(ji,1) + zkappa_s(ji,0) * zg1s * t_su_1d(ji) )  
    429503               ENDIF 
    430504            !                               !---------------------! 
    … …  
    433507               ! 
    434508               IF ( t_su_1d(ji) < rt0 ) THEN   !--  case 1 : no surface melting 
    435                   ! 
    436                   jm_min(ji) = nlay_s + 1 
    437                   jm_max(ji) = nlay_i + nlay_s + 1 
    438  
    439                   ! surface equation    
    440                   ztrid(ji,jm_min(ji),1)  = 0.0 
    441                   ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
    442                   ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1 
    443                   zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 
    444  
    445                   ! first layer of ice equation 
    446                   ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
    447                   ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
    448                   ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)   
    449                   zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
    450  
    451                   ! case of only one layer in the ice (surface & ice equations are altered) 
    452                   IF ( nlay_i == 1 ) THEN 
    453                      ztrid(ji,jm_min(ji),1)    = 0.0 
    454                      ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
    455                      ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0 
    456                      ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
    457                      ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    458                      ztrid(ji,jm_min(ji)+1,3)  = 0.0 
    459                      zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  & 
    460                         &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
    461                   ENDIF 
     509               ! 
     510               jm_min(ji) = nlay_s + 1 
     511               jm_max(ji) = nlay_i + nlay_s + 1 
     512 
     513               ! surface equation    
     514               ztrid(ji,jm_min(ji),1)  = 0.0 
     515               ztrid(ji,jm_min(ji),2)  = zdqns_ice_b(ji) - zkappa_i(ji,0)*zg1     
     516               ztrid(ji,jm_min(ji),3)  = zkappa_i(ji,0)*zg1 
     517               zindterm(ji,jm_min(ji)) = zdqns_ice_b(ji)*t_su_1d(ji) - zfnet(ji) 
     518 
     519               ! first layer of ice equation 
     520               ztrid(ji,jm_min(ji)+1,1)  = - zkappa_i(ji,0) * zg1 * zeta_i(ji,1) 
     521               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 ) 
     522               ztrid(ji,jm_min(ji)+1,3)  = - zeta_i(ji,1) * zkappa_i(ji,1)   
     523               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) * zradab_i(ji,1)   
     524 
     525               ! case of only one layer in the ice (surface & ice equations are altered) 
     526               IF ( nlay_i == 1 ) THEN 
     527               ztrid(ji,jm_min(ji),1)    = 0.0 
     528               ztrid(ji,jm_min(ji),2)    = zdqns_ice_b(ji) - zkappa_i(ji,0) * 2.0 
     529               ztrid(ji,jm_min(ji),3)    = zkappa_i(ji,0) * 2.0 
     530               ztrid(ji,jm_min(ji)+1,1)  = -zkappa_i(ji,0) * 2.0 * zeta_i(ji,1) 
     531               ztrid(ji,jm_min(ji)+1,2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     532               ztrid(ji,jm_min(ji)+1,3)  = 0.0 
     533               zindterm(ji,jm_min(ji)+1) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     534               &                      ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) ) 
     535               ENDIF 
    462536 
    463537               ELSE                            !--  case 2 : surface is melting 
    464538 
    465                   jm_min(ji)    =  nlay_s + 2 
    466                   jm_max(ji)    =  nlay_i + nlay_s + 1 
    467  
    468                   ! first layer of ice equation 
    469                   ztrid(ji,jm_min(ji),1)  = 0.0 
    470                   ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
    471                   ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
    472                   zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
    473                      &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
    474  
    475                   ! case of only one layer in the ice (surface & ice equations are altered) 
    476                   IF ( nlay_i == 1 ) THEN 
    477                      ztrid(ji,jm_min(ji),1)  = 0.0 
    478                      ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
    479                      ztrid(ji,jm_min(ji),3)  = 0.0 
    480                      zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
    481                         &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
    482                   ENDIF 
     539               jm_min(ji)    =  nlay_s + 2 
     540               jm_max(ji)    =  nlay_i + nlay_s + 1 
     541 
     542               ! first layer of ice equation 
     543               ztrid(ji,jm_min(ji),1)  = 0.0 
     544               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,1) + zkappa_i(ji,0) * zg1 )   
     545               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
     546               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     547               &                    ( zradab_i(ji,1) + zkappa_i(ji,0) * zg1 * t_su_1d(ji) )  
     548 
     549               ! case of only one layer in the ice (surface & ice equations are altered) 
     550               IF ( nlay_i == 1 ) THEN 
     551               ztrid(ji,jm_min(ji),1)  = 0.0 
     552               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * ( zkappa_i(ji,0) * 2.0 + zkappa_i(ji,1) ) 
     553               ztrid(ji,jm_min(ji),3)  = 0.0 
     554               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji) )  & 
     555               &                    + t_su_1d(ji) * zeta_i(ji,1) * zkappa_i(ji,0) * 2.0 
     556               ENDIF 
    483557 
    484558               ENDIF 
    … …  
    488562            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
    489563            ! 
    490          END DO 
    491          ! 
    492          !------------------------------ 
    493          ! 8) tridiagonal system solving 
    494          !------------------------------ 
    495          ! Solve the tridiagonal system with Gauss elimination method. 
    496          ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
    497          jm_maxt = 0 
    498          jm_mint = nlay_i+5 
    499          DO ji = 1, npti 
     564            END DO 
     565            ! 
     566            !------------------------------ 
     567            ! 8) tridiagonal system solving 
     568            !------------------------------ 
     569            ! Solve the tridiagonal system with Gauss elimination method. 
     570            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     571            jm_maxt = 0 
     572            jm_mint = nlay_i+5 
     573            DO ji = 1, npti 
    500574            jm_mint = MIN(jm_min(ji),jm_mint) 
    501575            jm_maxt = MAX(jm_max(ji),jm_maxt) 
    502          END DO 
    503  
    504          DO jk = jm_mint+1, jm_maxt 
     576            END DO 
     577 
     578            DO jk = jm_mint+1, jm_maxt 
    505579            DO ji = 1, npti 
    506580               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
    … …  
    508582               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
    509583            END DO 
    510          END DO 
    511  
    512          DO ji = 1, npti 
     584            END DO 
     585 
     586            DO ji = 1, npti 
    513587            ! ice temperatures 
    514588            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
    515          END DO 
    516  
    517          DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
     589            END DO 
     590 
     591            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
    518592            DO ji = 1, npti 
    519593               jk    =  jm - nlay_s - 1 
    520594               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
    521595            END DO 
    522          END DO 
    523  
    524          DO ji = 1, npti 
     596            END DO 
     597 
     598            DO ji = 1, npti 
    525599            ! snow temperatures       
    526600            IF( h_s_1d(ji) > 0._wp ) THEN 
    527601               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
    528                   &                / zdiagbis(ji,nlay_s+1) 
     602               &                / zdiagbis(ji,nlay_s+1) 
    529603            ENDIF 
    530604            ! surface temperature 
    … …  
    532606            IF( t_su_1d(ji) < rt0 ) THEN 
    533607               t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) *  & 
    534                   &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
    535             ENDIF 
    536          END DO 
    537          ! 
    538          !-------------------------------------------------------------- 
    539          ! 9) Has the scheme converged ?, end of the iterative procedure 
    540          !-------------------------------------------------------------- 
    541          ! check that nowhere it has started to melt 
    542          ! zdti_max is a measure of error, it has to be under zdti_bnd 
    543          zdti_max = 0._wp 
    544          DO ji = 1, npti 
     608               &          ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 
     609            ENDIF 
     610            END DO 
     611            ! 
     612            !-------------------------------------------------------------- 
     613            ! 9) Has the scheme converged ?, end of the iterative procedure 
     614            !-------------------------------------------------------------- 
     615            ! check that nowhere it has started to melt 
     616            ! zdti_max is a measure of error, it has to be under zdti_bnd 
     617            zdti_max = 0._wp 
     618            DO ji = 1, npti 
    545619            t_su_1d(ji) = MAX( MIN( t_su_1d(ji) , rt0 ) , rt0 - 100._wp ) 
    546620            zdti_max = MAX( zdti_max, ABS( t_su_1d(ji) - ztsub(ji) ) )      
    547          END DO 
    548  
    549          DO jk  =  1, nlay_s 
     621            END DO 
     622 
     623            DO jk  =  1, nlay_s 
    550624            DO ji = 1, npti 
    551625               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
    552626               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
    553627            END DO 
    554          END DO 
    555  
    556          DO jk  =  1, nlay_i 
     628            END DO 
     629 
     630            DO jk  =  1, nlay_i 
    557631            DO ji = 1, npti 
    558632               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
    … …  
    560634               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
    561635            END DO 
    562          END DO 
    563  
    564          ! Compute spatial maximum over all errors 
    565          ! note that this could be optimized substantially by iterating only the non-converging points 
    566          IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
    567  
     636            END DO 
     637 
     638            ! Compute spatial maximum over all errors 
     639            ! note that this could be optimized substantially by iterating only the non-converging points 
     640            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     641         ! 
     642         !----------------------------------------! 
     643         !                                        ! 
     644         !       JULES IF (RCV MODE)              ! 
     645         !                                        ! 
     646         !----------------------------------------! 
     647         ! 
     648 
     649         ELSE IF ( k_jules == np_zdf_jules_RCV  ) THEN ! RCV mode 
     650          
     651            ! 
     652            ! In RCV mode, we use a modified BL99 solver  
     653            ! with conduction flux (qcn_ice) as forcing term 
     654            ! 
     655            !---------------------------- 
     656            ! 7) tridiagonal system terms 
     657            !---------------------------- 
     658            !!layer denotes the number of the layer in the snow or in the ice 
     659            !!jm denotes the reference number of the equation in the tridiagonal 
     660            !!system, terms of tridiagonal system are indexed as following : 
     661            !!1 is subdiagonal term, 2 is diagonal and 3 is superdiagonal one 
     662 
     663            !!ice interior terms (top equation has the same form as the others) 
     664            ztrid   (1:npti,:,:) = 0._wp 
     665            zindterm(1:npti,:)   = 0._wp 
     666            zindtbis(1:npti,:)   = 0._wp 
     667            zdiagbis(1:npti,:)   = 0._wp 
     668 
     669            DO jm = nlay_s + 2, nlay_s + nlay_i  
     670            DO ji = 1, npti 
     671               jk = jm - nlay_s - 1 
     672               ztrid(ji,jm,1)  = - zeta_i(ji,jk) * zkappa_i(ji,jk-1) 
     673               ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,jk) * ( zkappa_i(ji,jk-1) + zkappa_i(ji,jk) ) 
     674               ztrid(ji,jm,3)  = - zeta_i(ji,jk) * zkappa_i(ji,jk) 
     675               zindterm(ji,jm) = ztiold(ji,jk) + zeta_i(ji,jk) * zradab_i(ji,jk) 
     676            END DO 
     677            ENDDO 
     678 
     679            jm =  nlay_s + nlay_i + 1 
     680            DO ji = 1, npti 
     681            !!ice bottom term 
     682            ztrid(ji,jm,1)  = - zeta_i(ji,nlay_i)*zkappa_i(ji,nlay_i-1)    
     683            ztrid(ji,jm,2)  = 1.0 + zeta_i(ji,nlay_i) * ( zkappa_i(ji,nlay_i) * zg1 + zkappa_i(ji,nlay_i-1) ) 
     684            ztrid(ji,jm,3)  = 0.0 
     685            zindterm(ji,jm) = ztiold(ji,nlay_i) + zeta_i(ji,nlay_i) *  & 
     686               &            ( zradab_i(ji,nlay_i) + zkappa_i(ji,nlay_i) * zg1 *  t_bo_1d(ji) )  
     687            ENDDO 
     688 
     689 
     690            DO ji = 1, npti 
     691            !                              !---------------------! 
     692            IF ( h_s_1d(ji) > 0.0 ) THEN   !  snow-covered cells ! 
     693               !                           !---------------------! 
     694               ! snow interior terms (bottom equation has the same form as the others) 
     695               DO jm = 3, nlay_s + 1 
     696               jk = jm - 1 
     697               ztrid(ji,jm,1)  = - zeta_s(ji,jk) * zkappa_s(ji,jk-1) 
     698               ztrid(ji,jm,2)  = 1.0 + zeta_s(ji,jk) * ( zkappa_s(ji,jk-1) + zkappa_s(ji,jk) ) 
     699               ztrid(ji,jm,3)  = - zeta_s(ji,jk)*zkappa_s(ji,jk) 
     700               zindterm(ji,jm) = ztsold(ji,jk) + zeta_s(ji,jk) * zradab_s(ji,jk) 
     701               END DO 
     702 
     703               ! case of only one layer in the ice (ice equation is altered) 
     704               IF ( nlay_i == 1 ) THEN 
     705               ztrid(ji,nlay_s+2,3)  = 0.0 
     706               zindterm(ji,nlay_s+2) = zindterm(ji,nlay_s+2) + zkappa_i(ji,1) * t_bo_1d(ji)  
     707               ENDIF 
     708 
     709               jm_min(ji) = 2 
     710               jm_max(ji) = nlay_i + nlay_s + 1 
     711 
     712               ! first layer of snow equation 
     713               ztrid(ji,2,1)  = 0.0 
     714               ztrid(ji,2,2)  = 1.0 + zeta_s(ji,1) * zkappa_s(ji,1) 
     715               ztrid(ji,2,3)  = - zeta_s(ji,1)*zkappa_s(ji,1)  
     716               zindterm(ji,2) = ztsold(ji,1) + zeta_s(ji,1) *  & 
     717               &           ( zradab_s(ji,1) + qcn_ice_1d(ji) )  
     718 
     719            !                               !---------------------! 
     720            ELSE                            ! cells without snow  ! 
     721            !                               !---------------------! 
     722 
     723               jm_min(ji)    =  nlay_s + 2 
     724               jm_max(ji)    =  nlay_i + nlay_s + 1 
     725 
     726               ! first layer of ice equation 
     727               ztrid(ji,jm_min(ji),1)  = 0.0 
     728               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1)   
     729               ztrid(ji,jm_min(ji),3)  = - zeta_i(ji,1) * zkappa_i(ji,1) 
     730               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) *  & 
     731               &                    ( zradab_i(ji,1) + qcn_ice_1d(ji) ) 
     732 
     733               ! case of only one layer in the ice (surface & ice equations are altered) 
     734               IF ( nlay_i == 1 ) THEN 
     735               ztrid(ji,jm_min(ji),1)  = 0.0 
     736               ztrid(ji,jm_min(ji),2)  = 1.0 + zeta_i(ji,1) * zkappa_i(ji,1) 
     737               ztrid(ji,jm_min(ji),3)  = 0.0 
     738               zindterm(ji,jm_min(ji)) = ztiold(ji,1) + zeta_i(ji,1) * ( zradab_i(ji,1) + zkappa_i(ji,1) * t_bo_1d(ji)    & 
     739               &                    + qcn_ice_1d(ji) ) 
     740 
     741               ENDIF 
     742 
     743            ENDIF 
     744            ! 
     745            zindtbis(ji,jm_min(ji)) = zindterm(ji,jm_min(ji)) 
     746            zdiagbis(ji,jm_min(ji)) = ztrid(ji,jm_min(ji),2) 
     747            ! 
     748            END DO 
     749            ! 
     750            !------------------------------ 
     751            ! 8) tridiagonal system solving 
     752            !------------------------------ 
     753            ! Solve the tridiagonal system with Gauss elimination method. 
     754            ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 
     755            jm_maxt = 0 
     756            jm_mint = nlay_i+5 
     757            DO ji = 1, npti 
     758            jm_mint = MIN(jm_min(ji),jm_mint) 
     759            jm_maxt = MAX(jm_max(ji),jm_maxt) 
     760            END DO 
     761 
     762            DO jk = jm_mint+1, jm_maxt 
     763            DO ji = 1, npti 
     764               jm = min(max(jm_min(ji)+1,jk),jm_max(ji)) 
     765               zdiagbis(ji,jm) = ztrid(ji,jm,2)  - ztrid(ji,jm,1) * ztrid(ji,jm-1,3)  / zdiagbis(ji,jm-1) 
     766               zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 
     767            END DO 
     768            END DO 
     769 
     770            DO ji = 1, npti 
     771            ! ice temperatures 
     772            t_i_1d(ji,nlay_i) = zindtbis(ji,jm_max(ji)) / zdiagbis(ji,jm_max(ji)) 
     773            END DO 
     774 
     775            DO jm = nlay_i + nlay_s, nlay_s + 2, -1 
     776            DO ji = 1, npti 
     777               jk    =  jm - nlay_s - 1 
     778               t_i_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_i_1d(ji,jk+1) ) / zdiagbis(ji,jm) 
     779            END DO 
     780            END DO 
     781 
     782            DO ji = 1, npti 
     783            ! snow temperatures       
     784            IF( h_s_1d(ji) > 0._wp ) THEN 
     785               t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) )  & 
     786               &                / zdiagbis(ji,nlay_s+1) 
     787            ENDIF 
     788            END DO 
     789            ! 
     790            !-------------------------------------------------------------- 
     791            ! 9) Has the scheme converged ?, end of the iterative procedure 
     792            !-------------------------------------------------------------- 
     793            ! check that nowhere it has started to melt 
     794            ! zdti_max is a measure of error, it has to be under zdti_bnd 
     795            zdti_max = 0._wp 
     796 
     797            DO jk  =  1, nlay_s 
     798            DO ji = 1, npti 
     799               t_s_1d(ji,jk) = MAX( MIN( t_s_1d(ji,jk), rt0 ), rt0 - 100._wp ) 
     800               zdti_max = MAX( zdti_max, ABS( t_s_1d(ji,jk) - ztsb(ji,jk) ) ) 
     801            END DO 
     802            END DO 
     803 
     804            DO jk  =  1, nlay_i 
     805            DO ji = 1, npti 
     806               ztmelt_i      = -tmut * sz_i_1d(ji,jk) + rt0  
     807               t_i_1d(ji,jk) =  MAX( MIN( t_i_1d(ji,jk), ztmelt_i ), rt0 - 100._wp ) 
     808               zdti_max      =  MAX( zdti_max, ABS( t_i_1d(ji,jk) - ztib(ji,jk) ) ) 
     809            END DO 
     810            END DO 
     811 
     812            ! Compute spatial maximum over all errors 
     813            ! note that this could be optimized substantially by iterating only the non-converging points 
     814            IF( lk_mpp ) CALL mpp_max( zdti_max, kcom=ncomm_ice ) 
     815 
     816         ENDIF ! k_jules 
     817          
    568818      END DO  ! End of the do while iterative procedure 
    569  
     819       
    570820      IF( ln_icectl .AND. lwp ) THEN 
    571821         WRITE(numout,*) ' zdti_max : ', zdti_max 
    572822         WRITE(numout,*) ' iconv    : ', iconv 
    573823      ENDIF 
     824       
    574825      ! 
    575826      !----------------------------- 
    576827      ! 10) Fluxes at the interfaces 
    577828      !----------------------------- 
     829      ! 
     830      ! --- update conduction fluxes 
     831      ! 
    578832      DO ji = 1, npti 
    579          !                                ! surface ice conduction flux 
     833      !                                ! surface ice conduction flux 
    580834         fc_su(ji)   =  -           isnow(ji)   * zkappa_s(ji,0) * zg1s * (t_s_1d(ji,1) - t_su_1d(ji))   & 
    581835            &           - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_1d(ji,1) - t_su_1d(ji)) 
    582          !                                ! bottom ice conduction flux 
     836      !                                ! bottom ice conduction flux 
    583837         fc_bo_i(ji) =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_1d(ji) - t_i_1d(ji,nlay_i)) ) 
    584838      END DO 
    585  
    586       ! --- computes sea ice energy of melting compulsory for icethd_dh --- ! 
    587       CALL ice_thd_enmelt 
    588  
    589       ! --- diagnose the change in non-solar flux due to surface temperature change --- ! 
    590       IF ( ln_dqns_i ) THEN 
     839       
     840      ! 
     841      ! --- Diagnose the heat loss due to changing non-solar / conduction flux --- ! 
     842      ! 
     843      DO ji = 1, npti 
     844            IF ( k_jules == np_zdf_jules_OFF .OR. k_jules == np_zdf_jules_SND  ) THEN  
     845                   ! OFF or SND mode 
     846               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji)  
     847            ELSE   ! RCV mode 
     848               hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( fc_su(ji) - qcn_ice_1d(ji) ) * a_i_1d(ji)  
     849            ENDIF 
     850      END DO 
     851       
     852      ! 
     853      ! --- Diagnose the heat loss due to non-fully converged temperature solution (should not be above 10-4 W-m2) --- ! 
     854      ! 
     855 
     856      IF ( ( k_jules == np_zdf_jules_OFF ) .OR. ( k_jules == np_zdf_jules_RCV ) ) THEN ! OFF 
     857       
     858         CALL ice_thd_enmelt        
     859 
     860         !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    591861         DO ji = 1, npti 
    592             hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) - ( qns_ice_1d(ji)  - zqns_ice_b(ji) ) * a_i_1d(ji)  
     862            zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
     863               &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
     864                
     865            IF ( ( k_jules == np_zdf_jules_OFF ) ) THEN 
     866                
     867                  IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
     868                     zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji)    - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)  
     869                  ELSE                          ! case T_su = 0degC 
     870                     zhfx_err = ( fc_su(ji)      + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
     871                  ENDIF 
     872             
     873            ELSE ! RCV CASE 
     874             
     875               zhfx_err = ( fc_su(ji) + qsr_ice_tr_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
     876             
     877            ENDIF 
     878             
     879            ! total heat sink to be sent to the ocean 
     880            hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
     881 
     882            ! hfx_dif = Heat flux diagnostic of sensible heat used to warm/cool ice in W.m-2    
     883            hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
     884    
     885         END DO  
     886            
     887         ! 
     888         ! --- SIMIP diagnostics 
     889         ! 
     890 
     891         DO ji = 1, npti 
     892            !--- Conduction fluxes (positive downwards) 
     893            diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
     894            diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
     895    
     896            !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
     897            zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
     898            IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
     899               t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
     900                  &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
     901            ELSE 
     902               t_si_1d(ji) = t_su_1d(ji) 
     903            ENDIF 
    593904         END DO 
    594       END IF 
    595  
    596       ! --- diag conservation imbalance on heat diffusion - PART 2 --- ! 
    597       !     hfx_dif = Heat flux used to warm/cool ice in W.m-2 
    598       !     zhfx_err = correction on the diagnosed heat flux due to non-convergence of the algorithm used to solve heat equation 
    599       DO ji = 1, npti 
    600          zdq = - zq_ini(ji) + ( SUM( e_i_1d(ji,1:nlay_i) ) * h_i_1d(ji) * r1_nlay_i +  & 
    601             &                   SUM( e_s_1d(ji,1:nlay_s) ) * h_s_1d(ji) * r1_nlay_s ) 
    602  
    603          IF( t_su_1d(ji) < rt0 ) THEN  ! case T_su < 0degC 
    604             zhfx_err = ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji)  
    605          ELSE                          ! case T_su = 0degC 
    606             zhfx_err = ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) + zdq * r1_rdtice ) * a_i_1d(ji) 
    607          ENDIF 
    608          hfx_dif_1d(ji) = hfx_dif_1d(ji) - zdq * r1_rdtice * a_i_1d(ji) 
    609  
    610          ! total heat that is sent to the ocean (i.e. not used in the heat diffusion equation) 
    611          hfx_err_dif_1d(ji) = hfx_err_dif_1d(ji) + zhfx_err 
    612  
    613       END DO    
    614  
    615       ! --- Diagnostics SIMIP --- ! 
    616       DO ji = 1, npti 
    617          !--- Conduction fluxes (positive downwards) 
    618          diag_fc_bo_1d(ji) = diag_fc_bo_1d(ji) + fc_bo_i(ji) * a_i_1d(ji) / at_i_1d(ji) 
    619          diag_fc_su_1d(ji) = diag_fc_su_1d(ji) + fc_su(ji)   * a_i_1d(ji) / at_i_1d(ji) 
    620  
    621          !--- Snow-ice interfacial temperature (diagnostic SIMIP) 
    622          zfac = rn_cnd_s * zh_i(ji) + ztcond_i(ji,1) * zh_s(ji) 
    623          IF( zh_s(ji) >= 1.e-3 .AND. zfac > epsi10 ) THEN 
    624             t_si_1d(ji) = ( rn_cnd_s       * zh_i(ji) * t_s_1d(ji,1) +   & 
    625                &            ztcond_i(ji,1) * zh_s(ji) * t_i_1d(ji,1) ) / zfac 
    626          ELSE 
    627             t_si_1d(ji) = t_su_1d(ji) 
    628          ENDIF 
    629       END DO 
    630       ! 
    631    END SUBROUTINE ice_thd_zdf 
     905      
     906      ENDIF 
     907      ! 
     908      !--------------------------------------------------------------------------------------- 
     909      ! 11) Jules coupling: reset inner snow and ice temperatures, update conduction fluxes 
     910      !--------------------------------------------------------------------------------------- 
     911      ! 
     912      IF ( k_jules == np_zdf_jules_SND ) THEN   ! --- Jules coupling in "SND" mode 
     913       
     914         ! Restore temperatures to their initial values 
     915         t_s_1d(1:npti,:)        = ztsold (1:npti,:) 
     916         t_i_1d(1:npti,:)        = ztiold (1:npti,:) 
     917         qcn_ice_1d(1:npti)      = fc_su(1:npti) 
     918          
     919      ENDIF 
     920       
     921   END SUBROUTINE ice_thd_zdf_BL99 
     922 
    632923 
    633924 
    … …  
    663954 
    664955 
     956 
    665957   SUBROUTINE ice_thd_zdf_init 
    666958      !!----------------------------------------------------------------------- 
    … …  
    675967      !! ** input   :   Namelist namthd_zdf 
    676968      !!------------------------------------------------------------------- 
    677       INTEGER  ::   ios   ! Local integer output status for namelist read 
     969      INTEGER  ::   ios, ioptio   ! Local integer output status for namelist read 
    678970      !! 
    679       NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i, ln_dqns_i 
     971      NAMELIST/namthd_zdf/ ln_zdf_BL99, ln_cndi_U64, ln_cndi_P07, rn_cnd_s, rn_kappa_i 
    680972      !!------------------------------------------------------------------- 
    681973      ! 
    … …  
    694986         WRITE(numout,*) '~~~~~~~~~~~~~~~~' 
    695987         WRITE(numout,*) '   Namelist namthd_zdf:' 
    696          WRITE(numout,*) '      Diffusion follows a Bitz and Lipscomb (1999)            ln_zdf_BL99  = ', ln_zdf_BL99 
     988         WRITE(numout,*) '      Bitz and Lipscomb (1999) formulation                    ln_zdf_BL99  = ', ln_zdf_BL99 
    697989         WRITE(numout,*) '      thermal conductivity in the ice (Untersteiner 1964)     ln_cndi_U64  = ', ln_cndi_U64 
    698990         WRITE(numout,*) '      thermal conductivity in the ice (Pringle et al 2007)    ln_cndi_P07  = ', ln_cndi_P07 
    699991         WRITE(numout,*) '      thermal conductivity in the snow                        rn_cnd_s     = ', rn_cnd_s 
    700992         WRITE(numout,*) '      extinction radiation parameter in sea ice               rn_kappa_i   = ', rn_kappa_i 
    701          WRITE(numout,*) '      change the surface non-solar flux with Tsu or not       ln_dqns_i    = ', ln_dqns_i 
    702993     ENDIF 
     994 
    703995      ! 
    704996      IF ( ( ln_cndi_U64 .AND. ln_cndi_P07 ) .OR. ( .NOT.ln_cndi_U64 .AND. .NOT.ln_cndi_P07 ) ) THEN 
    705997         CALL ctl_stop( 'ice_thd_zdf_init: choose one and only one formulation for thermal conduction (ln_cndi_U64 or ln_cndi_P07)' ) 
    706998      ENDIF 
     999 
     1000      !                             !== set the choice of ice vertical thermodynamic formulation ==! 
     1001      ioptio = 0  
     1002      !      !--- BL99 thermo dynamics                               (linear liquidus + constant thermal properties) 
     1003      IF( ln_zdf_BL99 ) THEN   ;   ioptio = ioptio + 1   ;   nice_zdf = np_BL99   ;   ENDIF 
     1004      IF( ioptio /= 1 )    CALL ctl_stop( 'ice_thd_init: one and only one ice thermo option has to be defined ' ) 
    7071005      ! 
    7081006   END SUBROUTINE ice_thd_zdf_init 
    … …  
    7111009   !!---------------------------------------------------------------------- 
    7121010   !!   Default option       Dummy Module             No ESIM sea-ice model 
    713    !!---------------------------------------------------------------------- 
     1011   !!--------------------------------------------------------------------- 
    7141012#endif 
    7151013 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/iceupdate.F90

    r8738 r8752  
    164164            ! mass flux from ice/ocean 
    165165            wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj)   & 
    166                &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj)  
    167  
    168             IF ( ln_pnd_fw )   wfx_ice(ji,jj) = wfx_ice(ji,jj) + wfx_pnd(ji,jj) 
     166               &           + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj) 
    169167 
    170168            ! add the snow melt water to snow mass flux to the ocean 
    … …  
    199197      ! Snow/ice albedo (only if sent to coupler, useless in forced mode) 
    200198      !------------------------------------------------------------------ 
    201       CALL ice_alb( t_su, h_i, h_s, a_ip_frac, h_ip, ln_pnd_rad, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
     199      CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_frac, h_ip, zalb_cs, zalb_os ) ! cloud-sky and overcast-sky ice albedos 
    202200      ! 
    203201      alb_ice(:,:,:) = ( 1._wp - cldf_ice ) * zalb_cs(:,:,:) + cldf_ice * zalb_os(:,:,:) 
    … …  
    246244      CALL iom_put( "vfxlam"     , wfx_lam              )        ! lateral melt  
    247245      CALL iom_put( "vfxice"     , wfx_ice              )        ! total ice growth/melt  
    248       IF ( ln_pnd )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux 
     246      IF ( ln_pnd_fwb )   CALL iom_put( "vfxpnd", wfx_pnd   )        ! melt pond water flux 
    249247 
    250248      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icevar.F90

    r8738 r8752  
    9797      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
    9898 
    99       ! MV MP 2016 
    100       IF ( ln_pnd ) THEN                     ! Melt pond 
    101          at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) 
    102          vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 
    103       ENDIF 
    104       ! END MP 2016 
     99      at_ip(:,:) = SUM( a_ip(:,:,:), dim=3 ) ! melt ponds 
     100      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) 
    105101 
    106102      ato_i(:,:) = 1._wp - at_i(:,:)    ! open water fraction   
    … …  
    167163!!                a criteria for icy area (i.e. a_i > epsi20 and v_i > epsi20 ) 
    168164 
    169       !------------------------------------------------------- 
    170       ! Ice thickness, snow thickness, ice salinity, ice age 
    171       !------------------------------------------------------- 
     165      !--------------------------------------------------------------- 
     166      ! Ice thickness, snow thickness, ice salinity, ice age and ponds 
     167      !--------------------------------------------------------------- 
    172168      !                                            !--- inverse of the ice area 
    173169      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 
    … …  
    178174      ELSEWHERE                      ;   z1_v_i(:,:,:) = 0._wp 
    179175      END WHERE 
    180       ! 
    181       h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)    !--- ice thickness 
     176      !                                           !--- ice thickness 
     177      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 
    182178 
    183179      zhmax    =          hi_max(jpl) 
    184180      z1_zhmax =  1._wp / hi_max(jpl)                
    185       WHERE( h_i(:,:,jpl) > zhmax )               !--- bound h_i by hi_max (i.e. 99 m) with associated update of ice area 
     181      WHERE( h_i(:,:,jpl) > zhmax )   ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 
    186182         h_i  (:,:,jpl) = zhmax 
    187183         a_i   (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax  
    188          z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl)          ! NB: v_i always /=0 as h_i > hi_max 
     184         z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) 
    189185      END WHERE 
    190  
    191       h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:)    !--- snow thickness 
    192        
    193       o_i(:,:,:)  = oa_i(:,:,:) * z1_a_i(:,:,:)    !--- ice age 
    194  
    195       IF( nn_icesal == 2 ) THEN                    !--- salinity (with a minimum value imposed everywhere) 
     186      !                                           !--- snow thickness 
     187      h_s(:,:,:) = v_s (:,:,:) * z1_a_i(:,:,:) 
     188      !                                           !--- ice age       
     189      o_i(:,:,:) = oa_i(:,:,:) * z1_a_i(:,:,:) 
     190      !                                           !--- pond fraction and thickness       
     191      a_ip_frac(:,:,:) = a_ip(:,:,:) * z1_a_i(:,:,:) 
     192      WHERE( a_ip_frac(:,:,:) > epsi20 )   ;   h_ip(:,:,:) = v_ip(:,:,:) * z1_a_i(:,:,:) / a_ip_frac(:,:,:) 
     193      ELSEWHERE                            ;   h_ip(:,:,:) = 0._wp 
     194      END WHERE 
     195      ! 
     196      !                                           !---  salinity (with a minimum value imposed everywhere)      
     197      IF( nn_icesal == 2 ) THEN 
    196198         WHERE( v_i(:,:,:) > epsi20 )   ;   s_i(:,:,:) = MAX( rn_simin , MIN( rn_simax, sv_i(:,:,:) * z1_v_i(:,:,:) ) ) 
    197199         ELSEWHERE                      ;   s_i(:,:,:) = rn_simin 
    198200         END WHERE 
    199201      ENDIF 
    200  
    201       CALL ice_var_salprof      ! salinity profile 
     202      CALL ice_var_salprof   ! salinity profile 
    202203 
    203204      !------------------- 
    … …  
    242243      vt_s (:,:) = SUM( v_s, dim=3 ) 
    243244      at_i (:,:) = SUM( a_i, dim=3 ) 
    244  
    245 ! MV MP 2016 
    246 ! probably should resum for melt ponds ??? 
    247  
    248245      ! 
    249246   END SUBROUTINE ice_var_glo2eqv 
    … …  
    258255      !!------------------------------------------------------------------- 
    259256      ! 
    260       v_i (:,:,:) = h_i(:,:,:) * a_i(:,:,:) 
    261       v_s (:,:,:) = h_s(:,:,:) * a_i(:,:,:) 
    262       sv_i(:,:,:) = s_i(:,:,:) * v_i(:,:,:) 
     257      v_i (:,:,:) = h_i (:,:,:) * a_i (:,:,:) 
     258      v_s (:,:,:) = h_s (:,:,:) * a_i (:,:,:) 
     259      sv_i(:,:,:) = s_i (:,:,:) * v_i (:,:,:) 
     260      v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    263261      ! 
    264262   END SUBROUTINE ice_var_eqv2glo 
    … …  
    478476               wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl)   * rhosn * r1_rdtice 
    479477               hfx_res(ji,jj)  = hfx_res(ji,jj) - (1._wp - zswitch(ji,jj) ) * e_s (ji,jj,1,jl) * r1_rdtice ! W.m-2 <0 
     478               IF( ln_pnd_fwb ) THEN  
     479                  wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 
     480               ENDIF 
    480481               !----------------------------------------------------------------- 
    481482               ! Zap snow energy  
    … …  
    498499               h_s (ji,jj,jl) = h_s (ji,jj,jl) * zswitch(ji,jj) 
    499500 
    500             END DO 
    501          END DO 
    502  
    503          IF( ln_pnd ) THEN  
    504             DO jj = 1 , jpj 
    505                DO ji = 1 , jpi 
    506                   IF( ln_pnd_fw )   & 
    507                      &   wfx_res(ji,jj)  = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_ip(ji,jj,jl) * rhofw * r1_rdtice 
    508                   a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
    509                   v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
    510                END DO 
    511             END DO 
    512          ENDIF 
    513           
     501               a_ip (ji,jj,jl) = a_ip (ji,jj,jl) * zswitch(ji,jj) 
     502               v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 
     503 
     504            END DO 
     505         END DO 
     506         
    514507      END DO  
    515508 
    … …  
    520513      ! open water = 1 if at_i=0 
    521514      WHERE( at_i(:,:) == 0._wp )   ato_i(:,:) = 1._wp 
     515 
    522516      ! 
    523517   END SUBROUTINE ice_var_zapsmall 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/LIM_SRC_3/icewri.F90

    r8738 r8752  
    128128      CALL iom_put( "snowvol"     , vt_s  * zswi      )        ! snow volume 
    129129       
    130       IF ( ln_pnd ) THEN 
    131          CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction 
    132          CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area 
    133       ENDIF 
     130      CALL iom_put( "iceamp"  , at_ip  * zswi        )   ! melt pond total fraction 
     131      CALL iom_put( "icevmp"  , vt_ip  * zswi        )   ! melt pond total volume per unit area 
    134132 
    135133      !---------------------------------- 
    … …  
    145143      IF ( iom_use('brinevol_cat') )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 )          ! brine volume 
    146144 
    147       IF ( ln_pnd ) THEN 
    148          IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories 
    149          IF ( iom_use('icevmp_cat') )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories 
    150          IF ( iom_use('icehmp_cat') )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories 
    151          IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories 
    152       ENDIF 
     145      IF ( iom_use('iceamp_cat') )  CALL iom_put( "iceamp_cat"     , a_ip       * zswi2   )       ! melt pond frac for categories 
     146      IF ( iom_use('icevmp_cat') )  CALL iom_put( "icevmp_cat"     , v_ip       * zswi2   )       ! melt pond frac for categories 
     147      IF ( iom_use('icehmp_cat') )  CALL iom_put( "icehmp_cat"     , h_ip       * zswi2   )       ! melt pond frac for categories 
     148      IF ( iom_use('iceafp_cat') )  CALL iom_put( "iceafp_cat"     , a_ip_frac  * zswi2   )       ! melt pond frac for categories 
    153149 
    154150      !-------------------------------- 
    … …  
    316312      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   & 
    317313      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    318  
    319       ! MV MP 2016 
    320       IF ( ln_pnd ) THEN 
    321          CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   & 
    322       &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    323          CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   & 
    324       &         jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    325       ENDIF 
    326       ! END MV MP 2016 
    327  
     314      CALL histdef( kid, "si_amp", "Melt pond fraction"      , "%"      ,   & 
     315      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     316      CALL histdef( kid, "si_vmp", "Melt pond volume"        ,  "m"     ,   & 
     317      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    328318      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   & 
    329319      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    … …  
    370360      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    371361 
    372       ! MV MP 2016 
    373       IF ( ln_pnd ) THEN 
    374          CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) ) 
    375          CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) ) 
    376       ENDIF 
    377       ! END MV MP 2016 
     362      CALL histwrite( kid, "si_amp", kt, at_ip         , jpi*jpj, (/1/) ) 
     363      CALL histwrite( kid, "si_vmp", kt, vt_ip         , jpi*jpj, (/1/) ) 
    378364 
    379365      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) ) 
    … …  
    384370      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) ) 
    385371      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) ) 
    386       IF ( ln_pnd ) & 
    387          CALL histwrite( kid, "vfxpnd", kt, wfx_pnd     , jpi*jpj, (/1/) ) 
     372      CALL histwrite( kid, "vfxpnd", kt, wfx_pnd        , jpi*jpj, (/1/) ) 
    388373 
    389374      CALL histwrite( kid, "sithicat", kt, h_i         , jpi*jpj*jpl, (/1/) )     

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/BDY/bdyice.F90

    r8738 r8752  
    1919   USE eosbn2          ! equation of state 
    2020   USE oce             ! ocean dynamics and tracers variables 
    21    USE ice             ! LIM_3 ice variables 
    22    USE icevar 
    23    USE icectl 
     21   USE ice             ! sea-ice: variables 
     22   USE icevar          ! sea-ice: operations 
     23   USE iceitd          ! sea-ice: rebining 
     24   USE icectl          ! sea-ice: control prints 
    2425   USE par_oce         ! ocean parameters 
    2526   USE dom_oce         ! ocean space and time domain variables  
    … …  
    7172      END DO 
    7273      ! 
    73                         CALL ice_var_zapsmall 
    74                         CALL ice_var_agg(1) 
     74      CALL ice_var_zapsmall 
     75      CALL ice_var_agg(1) 
    7576      IF( ln_icectl )   CALL ice_prt( kt, iiceprt, jiceprt, 1, ' - ice thermo bdy - ' ) 
    7677      ! 
    … …  
    248249      END DO !jl 
    249250      ! 
     251      ! --- In case categories are out of bounds, do a remapping --- ! 
     252      !     i.e. inputs have not the same ice thickness distribution  
     253      !          (set by rn_himean) than the regional simulation 
     254      IF( jpl > 1 )   CALL ice_itd_reb( kt ) 
    250255      !       
    251256      IF( nn_timing == 1 )   CALL timing_stop('bdy_ice_frs') 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbc_ice.F90

    r8738 r8752  
    4848   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   alb_ice        !: ice albedo                                       [-] 
    4949 
     50   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qml_ice        !: heat available for snow / ice surface melting     [W/m2]  
     51   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice        !: heat conduction flux in the layer below surface   [W/m2]  
     52   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qsr_ice_tr     !: solar flux transmitted below the ice surface      [W/m2] 
     53 
    5054   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   utau_ice       !: atmos-ice u-stress. VP: I-pt ; EVP: U,V-pts   [N/m2] 
    5155   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vtau_ice       !: atmos-ice v-stress. VP: I-pt ; EVP: U,V-pts   [N/m2] 
    52    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr1_i0         !: Solar surface transmission parameter, thick ice  [-] 
    53    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fr2_i0         !: Solar surface transmission parameter, thin ice   [-] 
    5456   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   emp_ice        !: sublimation - precip over sea ice          [kg/m2/s] 
    5557 
    … …  
    123125      ALLOCATE( qns_ice (jpi,jpj,jpl) , qsr_ice (jpi,jpj,jpl) ,     & 
    124126         &      qla_ice (jpi,jpj,jpl) , dqla_ice(jpi,jpj,jpl) ,     & 
    125          &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl) , alb_ice (jpi,jpj,jpl) ,   & 
     127         &      dqns_ice(jpi,jpj,jpl) , tn_ice  (jpi,jpj,jpl)  , alb_ice (jpi,jpj,jpl) ,   & 
     128         &      qml_ice(jpi,jpj,jpl)  , qcn_ice(jpi,jpj,jpl)   , qsr_ice_tr(jpi,jpj,jpl),  & 
    126129         &      utau_ice(jpi,jpj)     , vtau_ice(jpi,jpj)     , wndm_ice(jpi,jpj)     ,   & 
    127          &      fr1_i0  (jpi,jpj)     , fr2_i0  (jpi,jpj)     ,     & 
    128130         &      evap_ice(jpi,jpj,jpl) , devap_ice(jpi,jpj,jpl) , qprec_ice(jpi,jpj) ,   & 
    129131         &      qemp_ice(jpi,jpj)     , qevap_ice(jpi,jpj,jpl) , qemp_oce (jpi,jpj) ,   & 
    … …  
    139141                a_i(jpi,jpj,ncat)     , topmelt(jpi,jpj,ncat) , botmelt(jpi,jpj,ncat) , & 
    140142                STAT= ierr(2) ) 
    141       IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , fr1_i0(jpi,jpj)       , tn_ice (jpi,jpj,1)    , & 
    142          &                     v_ice(jpi,jpj)        , fr2_i0(jpi,jpj)       , alb_ice(jpi,jpj,1)    , & 
     143      IF( ln_cpl )   ALLOCATE( u_ice(jpi,jpj)        , tn_ice (jpi,jpj,1)    , & 
     144         &                     v_ice(jpi,jpj)        , alb_ice(jpi,jpj,1)    , & 
    143145         &                     emp_ice(jpi,jpj)      , qns_ice(jpi,jpj,1)    , dqns_ice(jpi,jpj,1)   , & 
    144146         &                     STAT= ierr(3) )       
    … …  
    168170   REAL            , PUBLIC, PARAMETER ::   cldf_ice = 0.81       !: cloud fraction over sea ice, summer CLIO value   [-] 
    169171   INTEGER         , PUBLIC, PARAMETER ::   jpl = 1  
    170    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice,fr1_i0,fr2_i0          ! jpi, jpj 
     172   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_ice, v_ice                        ! jpi, jpj 
    171173   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   tn_ice, alb_ice, qns_ice, dqns_ice  ! (jpi,jpj,jpl) 
    172174   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbcblk.F90

    r8738 r8752  
    1717   !!                                          ==> based on AeroBulk (http://aerobulk.sourceforge.net/) 
    1818   !!            4.0  !  2016-10  (G. Madec)  introduce a sbc_blk_init routine 
     19   !!            4.0  !  2016-10  (M. Vancoppenolle)  Introduce Jules emulator (M. Vancoppenolle)  
    1920   !!---------------------------------------------------------------------- 
    2021 
    … …  
    4243   USE ice     , ONLY :   u_ice, v_ice, jpl, a_i_b, at_i_b, tm_su 
    4344   USE icethd_dh      ! for CALL ice_thd_snwblow 
     45   USE icethd_zdf, ONLY:   rn_cnd_s  ! for blk_ice_qcn 
    4446#endif 
    4547   USE sbcblk_algo_ncar     ! => turb_ncar     : NCAR - CORE (Large & Yeager, 2009)  
    … …  
    6466   PUBLIC   blk_ice_tau   ! routine called in icestp module 
    6567   PUBLIC   blk_ice_flx   ! routine called in icestp module 
    66 #endif 
     68   PUBLIC   blk_ice_qcn   ! routine called in icestp module 
     69#endif  
    6770 
    6871!!Lolo: should ultimately be moved in the module with all physical constants ? 
    … …  
    697700 
    698701 
    699    SUBROUTINE blk_ice_flx( ptsu, palb ) 
     702   SUBROUTINE blk_ice_flx( ptsu, phs, phi, palb ) 
    700703      !!--------------------------------------------------------------------- 
    701704      !!                     ***  ROUTINE blk_ice_flx  *** 
    … …  
    710713      !!--------------------------------------------------------------------- 
    711714      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   ptsu   ! sea ice surface temperature 
     715      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phs    ! snow thickness 
     716      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   phi    ! ice thickness 
    712717      REAL(wp), DIMENSION(:,:,:), INTENT(in)  ::   palb   ! ice albedo (all skies) 
    713718      !! 
    … …  
    716721      REAL(wp) ::   zcoef_dqlw, zcoef_dqla   !   -      - 
    717722      REAL(wp) ::   zztmp, z1_lsub           !   -      - 
     723      REAL(wp) ::   zfrqsr_tr_i              ! sea ice shortwave fraction transmitted below through the ice 
     724      REAL(wp) ::   zfr1, zfr2, zfac         ! local variables 
    718725      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qlw         ! long wave heat flux over ice 
    719726      REAL(wp), DIMENSION(:,:,:), POINTER ::   z_qsb         ! sensible  heat flux over ice 
    … …  
    722729      REAL(wp), DIMENSION(:,:)  , POINTER ::   zevap, zsnw   ! evaporation and snw distribution after wind blowing (LIM3) 
    723730      REAL(wp), DIMENSION(:,:)  , POINTER ::   zrhoa 
     731 
    724732      !!--------------------------------------------------------------------- 
    725733      ! 
    … …  
    830838      CALL wrk_dealloc( jpi,jpj,   zevap, zsnw ) 
    831839 
    832       !-------------------------------------------------------------------- 
    833       ! FRACTIONs of net shortwave radiation which is not absorbed in the 
    834       ! thin surface layer and penetrates inside the ice cover 
    835       ! ( Maykut and Untersteiner, 1971 ; Ebert and Curry, 1993 ) 
    836       ! 
    837       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    838       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
    839       ! 
    840       ! 
     840      ! --- absorbed and transmitted shortwave radiation (W/m2) --- ! 
     841      ! 
     842      ! former coding was 
     843      !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     844      !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     845 
     846      ! --- surface transmission parameter (i0, Grenfell Maykut 77) --- ! 
     847      zfr1 = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice )            !   standard value 
     848      zfr2 = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice )            !   zfr2 such that zfr1 + zfr2 to equal 1 
     849 
     850      qsr_ice_tr(:,:,:) = 0._wp 
     851 
     852      DO jl = 1, jpl 
     853         DO jj = 1, jpj 
     854            DO ji = 1, jpi 
     855 
     856               zfac = MAX( 0._wp , 1._wp - phi(ji,jj,jl) * 10._wp )     !   linear weighting factor: =0 for phi=0, =1 for phi = 0.1 
     857               zfrqsr_tr_i = zfr1 + zfac * zfr2                         !   below 10 cm, linearly increase zfrqsr_tr_i until 1 at zero thickness 
     858 
     859               IF ( phs(ji,jj,jl) <= 0.0_wp ) THEN  
     860                  zfrqsr_tr_i  = zfr1            + zfac * zfr2 
     861               ELSE  
     862                  zfrqsr_tr_i  = 0._wp                                  !   snow fully opaque 
     863               ENDIF 
     864 
     865               qsr_ice_tr(ji,jj,jl) = zfrqsr_tr_i * qsr_ice(ji,jj,jl)   !   transmitted solar radiation  
     866 
     867            END DO 
     868         END DO 
     869      END DO 
     870                                                                                       
     871 
    841872      IF(ln_ctl) THEN 
    842873         CALL prt_ctl(tab3d_1=qla_ice , clinfo1=' blk_ice: qla_ice  : ', tab3d_2=z_qsb   , clinfo2=' z_qsb    : ', kdim=jpl) 
    … …  
    854885 
    855886   END SUBROUTINE blk_ice_flx 
     887    
     888    
     889 
     890   SUBROUTINE blk_ice_qcn( k_monocat, ptsu, ptb, phs, phi ) 
     891    
     892      !!--------------------------------------------------------------------- 
     893      !!                     ***  ROUTINE blk_ice_qcn  *** 
     894      !! 
     895      !! ** Purpose :   Compute surface temperature and snow/ice conduction flux 
     896      !!                to force sea ice / snow thermodynamics 
     897      !!                in the case JULES coupler is emulated 
     898      !!                 
     899      !! ** Method  :   compute surface energy balance assuming neglecting heat storage 
     900      !!                following the 0-layer Semtner (1976) approach 
     901      !! 
     902      !! ** Outputs : - ptsu    : sea-ice / snow surface temperature (K) 
     903      !!              - qcn_ice : surface inner conduction flux (W/m2) 
     904      !! 
     905      !!--------------------------------------------------------------------- 
     906      !! 
     907      INTEGER, INTENT(in)                        ::   k_monocat                 ! single-category option 
     908 
     909      REAL(wp), DIMENSION(:,:,:), INTENT(inout)  ::   ptsu                      ! sea ice / snow surface temperature 
     910 
     911      REAL(wp), DIMENSION(:,:)  , INTENT(in)     ::   ptb                       ! sea ice base temperature 
     912      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phs                       ! snow thickness 
     913      REAL(wp), DIMENSION(:,:,:), INTENT(in)     ::   phi                       ! sea ice thickness 
     914       
     915      INTEGER                                    ::   ji, jj, jl                ! dummy loop indices 
     916      INTEGER                                    ::   iter                      ! 
     917      REAL(wp)                                   ::   zfac, zfac2, zfac3        ! dummy factors 
     918      REAL(wp)                                   ::   zkeff_h, ztsu             ! 
     919      REAL(wp)                                   ::   zqc, zqnet                ! 
     920      REAL(wp)                                   ::   zhe, zqa0                 ! 
     921            
     922      INTEGER , PARAMETER                        ::   nit = 10                  ! number of iterations 
     923      REAL(wp), PARAMETER                        ::   zepsilon = 0.1_wp         ! characteristic thickness for enhanced conduction 
     924 
     925      REAL(wp), DIMENSION(:,:,:), POINTER        ::   zgfac                    ! enhanced conduction factor 
     926       
     927      !!--------------------------------------------------------------------- 
     928 
     929      IF( nn_timing == 1 )  CALL timing_start('blk_ice_qcn') 
     930      ! 
     931      CALL wrk_alloc( jpi,jpj,jpl,   zgfac ) 
     932 
     933      ! -------------------------------------! 
     934      !      I   Enhanced conduction factor  ! 
     935      ! -------------------------------------! 
     936      ! 
     937      ! Emulates the enhancement of conduction by unresolved thin ice (k_monocat = 1/3) 
     938      ! Fichefet and Morales Maqueda, JGR 1997 
     939      ! 
     940      zgfac(:,:,:) = 1._wp 
     941       
     942      SELECT CASE ( k_monocat )  
     943       
     944         CASE ( 1 , 3 ) 
     945       
     946            zfac    =   1._wp /  ( rn_cnd_s + rcdic ) 
     947            zfac2   =   EXP(1._wp) * 0.5_wp * zepsilon 
     948            zfac3   =   2._wp / zepsilon 
     949       
     950            DO jl = 1, jpl                 
     951               DO jj = 1 , jpj 
     952                  DO ji = 1, jpi 
     953                     !                                ! Effective thickness 
     954                     zhe               =   ( rn_cnd_s * phi(ji,jj,jl) + rcdic * phs(ji,jj,jl) ) * zfac 
     955                     !                                ! Enhanced conduction factor 
     956                     IF( zhe >=  zfac2 )  & 
     957                        zgfac(ji,jj,jl)   =   MIN( 2._wp, ( 0.5_wp + 0.5 * LOG( zhe * zfac3 ) ) ) 
     958                  END DO 
     959               END DO 
     960            END DO 
     961       
     962      END SELECT 
     963       
     964      ! -------------------------------------------------------------! 
     965      !      II   Surface temperature and conduction flux            ! 
     966      ! -------------------------------------------------------------! 
     967 
     968      zfac   =   rcdic * rn_cnd_s  
     969                                      ! ========================== ! 
     970      DO jl = 1, jpl                  !  Loop over ice categories  ! 
     971         !                            ! ========================== ! 
     972         DO jj = 1 , jpj 
     973            DO ji = 1, jpi 
     974               !                      ! Effective conductivity of the snow-ice system divided by thickness 
     975               zkeff_h   =   zfac * zgfac(ji,jj,jl) / ( rcdic * phs(ji,jj,jl) + rn_cnd_s * phi(ji,jj,jl) ) 
     976                                      ! Store initial surface temperature 
     977               ztsu      =   ptsu(ji,jj,jl) 
     978                                      ! Net initial atmospheric heat flux 
     979               zqa0      =   qsr_ice(ji,jj,jl) - qsr_ice_tr(ji,jj,jl) + qns_ice(ji,jj,jl) 
     980                
     981               DO iter = 1, nit       ! --- Iteration loop 
     982                
     983                   !                  ! Conduction heat flux through snow-ice system (>0 downwards) 
     984                   zqc   =   zkeff_h * ( ztsu - ptb(ji,jj) ) 
     985                   !                  ! Surface energy budget 
     986                   zqnet =   zqa0 + dqns_ice(ji,jj,jl) * ( ztsu - ptsu(ji,jj,jl) ) - zqc 
     987                   !                  ! Temperature update 
     988                   ztsu  =   ztsu - zqnet / ( dqns_ice(ji,jj,jl) - zkeff_h ) 
     989                
     990               END DO 
     991                
     992               ptsu(ji,jj,jl)    = MIN( rt0, ztsu ) 
     993                
     994               qcn_ice(ji,jj,jl) = zkeff_h * ( ptsu(ji,jj,jl) - ptb(ji,jj) ) 
     995                
     996            END DO 
     997         END DO 
     998          
     999      END DO  
     1000       
     1001      CALL wrk_dealloc( jpi,jpj,jpl,   zgfac ) 
     1002       
     1003      IF( nn_timing == 1 )  CALL timing_stop('blk_ice_qcn') 
     1004 
     1005   END SUBROUTINE blk_ice_qcn 
    8561006    
    8571007#endif 
    … …  
    10931243      zqi_sat(:,:) = 0.98_wp * q_sat( tm_su(:,:), sf(jp_slp)%fnow(:,:,1) )  ! saturation humidity over ice   [kg/kg] 
    10941244      ! 
    1095 !!      DO jj = 2, jpjm1 
    1096 !!         DO ji = fs_2, fs_jpim1 
    1097       DO jj = 1, jpj 
    1098          DO ji = 1, jpi 
     1245      DO jj = 2, jpjm1           ! reduced loop is necessary for reproducibility 
     1246         DO ji = fs_2, fs_jpim1 
    10991247            ! Virtual potential temperature [K] 
    11001248            zthetav_os = zst(ji,jj)   * ( 1._wp + rctv0 * zqo_sat(ji,jj) )   ! over ocean 
    … …  
    11401288         END DO 
    11411289      END DO 
    1142 !!      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
     1290      CALL lbc_lnk_multi( Cd, 'T',  1., Ch, 'T', 1. ) 
    11431291      ! 
    11441292   END SUBROUTINE Cdn10_Lupkes2015 

  • branches/UKMO/dev_r8183_ICEMODEL_svn_removed/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r8751 r8752  
    3232   USE geo2ocean      !  
    3333   USE oce   , ONLY : tsn, un, vn, sshn, ub, vb, sshb, fraqsr_1lev 
    34    USE albedooce      !  
     34   USE ocealb         !  
    3535   USE eosbn2         !  
    3636   USE sbcrnf, ONLY : l_rnfcpl 
    … …  
    178178   TYPE( DYNARR ), SAVE, DIMENSION(jprcv) ::   frcv                     ! all fields recieved from the atmosphere 
    179179 
    180    REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   albedo_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
     180   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   alb_oce_mix    ! ocean albedo sent to atmosphere (mix clear/overcast sky) 
    181181 
    182182   REAL(wp) ::   rpref = 101000._wp   ! reference atmospheric pressure[N/m2]  
    … …  
    202202      ierr(:) = 0 
    203203      ! 
    204       ALLOCATE( albedo_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) ) 
     204      ALLOCATE( alb_oce_mix(jpi,jpj), nrcvinfo(jprcv),  STAT=ierr(1) ) 
    205205       
    206206#if ! defined key_lim3 && ! defined key_cice 
    … …  
    737737      !     2. receiving mixed oce-ice solar radiation  
    738738      IF ( TRIM ( sn_snd_alb%cldes ) == 'mixed oce-ice' .OR. TRIM ( sn_rcv_qsr%cldes ) == 'mixed oce-ice' ) THEN 
    739          CALL albedo_oce( zaos, zacs ) 
     739         CALL oce_alb( zaos, zacs ) 
    740740         ! Due to lack of information on nebulosity : mean clear/overcast sky 
    741          albedo_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5 
     741         alb_oce_mix(:,:) = ( zacs(:,:) + zaos(:,:) ) * 0.5 
    742742      ENDIF 
    743743 
    … …  
    15301530    
    15311531 
    1532    SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist ) 
     1532   SUBROUTINE sbc_cpl_ice_flx( picefr, palbi, psst, pist, phs, phi ) 
    15331533      !!---------------------------------------------------------------------- 
    15341534      !!             ***  ROUTINE sbc_cpl_ice_flx  *** 
    … …  
    15851585      REAL(wp), INTENT(in), DIMENSION(:,:  ), OPTIONAL ::   psst       ! sea surface temperature     [Celsius] 
    15861586      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   pist       ! ice surface temperature     [Kelvin] 
    1587       ! 
    1588       INTEGER ::   jl         ! dummy loop index 
     1587      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phs        ! snow depth                  [m] 
     1588      REAL(wp), INTENT(in), DIMENSION(:,:,:), OPTIONAL ::   phi        ! ice thickness               [m] 
     1589      ! 
     1590      INTEGER ::   ji,jj,jl         ! dummy loop index 
    15891591      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw 
    15901592      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice 
    15911593      REAL(wp), POINTER, DIMENSION(:,:  ) ::   zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice 
    1592       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice 
     1594      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i 
    15931595      !!---------------------------------------------------------------------- 
    15941596      ! 
    … …  
    15981600      CALL wrk_alloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    15991601      CALL wrk_alloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    1600       CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
     1602      CALL wrk_alloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i ) 
    16011603 
    16021604      IF( ln_mixcpl )   zmsk(:,:) = 1. - xcplmask(:,:,0) 
    … …  
    18901892!       ( see OASIS3 user guide, 5th edition, p39 ) 
    18911893         zqsr_ice(:,:,1) = frcv(jpr_qsrmix)%z3(:,:,1) * ( 1.- palbi(:,:,1) )   & 
    1892             &            / (  1.- ( albedo_oce_mix(:,:  ) * ziceld(:,:)       & 
     1894            &            / (  1.- ( alb_oce_mix(:,:  ) * ziceld(:,:)       & 
    18931895            &                     + palbi         (:,:,1) * picefr(:,:) ) ) 
    18941896      END SELECT 
    … …  
    19511953      END SELECT 
    19521954 
    1953       ! Surface transimission parameter io (Maykut Untersteiner , 1971 ; Ebert and Curry, 1993 ) 
    1954       ! Used for LIM3 
    1955       ! Coupled case: since cloud cover is not received from atmosphere  
    1956       !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
    1957       fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
    1958       fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     1955      ! --- Transmitted shortwave radiation (W/m2) --- ! 
     1956       
     1957      IF ( nice_jules == 0 ) THEN 
     1958                
     1959         zfrqsr_tr_i(:,:,:) = 0._wp   !   surface transmission parameter 
     1960      
     1961         ! former coding was     
     1962         ! Coupled case: since cloud cover is not received from atmosphere  
     1963         !               ===> used prescribed cloud fraction representative for polar oceans in summer (0.81) 
     1964         !     fr1_i0(:,:) = ( 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice ) 
     1965         !     fr2_i0(:,:) = ( 0.82 * ( 1.0 - cldf_ice ) + 0.65 * cldf_ice ) 
     1966          
     1967         ! to retrieve that coding, we needed to access h_i & h_s from here 
     1968         ! we could even retrieve cloud fraction from the coupler 
     1969                
     1970         DO jl = 1, jpl 
     1971            DO jj = 1 , jpj 
     1972               DO ji = 1, jpi 
     1973             
     1974                  !--- surface transmission parameter (Grenfell Maykut 77) --- ! 
     1975                  zfrqsr_tr_i(ji,jj,jl) = 0.18 * ( 1.0 - cldf_ice ) + 0.35 * cldf_ice  
     1976                
     1977                  ! --- influence of snow and thin ice --- ! 
     1978                  IF ( phs(ji,jj,jl) >= 0.0_wp ) zfrqsr_tr_i(ji,jj,jl) = 0._wp   !   snow fully opaque 
     1979                  IF ( phi(ji,jj,jl) <= 0.1_wp ) zfrqsr_tr_i(ji,jj,jl) = 1._wp   !   thin ice transmits all solar radiation 
     1980               END DO 
     1981            END DO 
     1982         END DO 
     1983          
     1984         qsr_ice_tr(:,:,:) =   zfrqsr_tr_i(:,:,:) * qsr_ice(:,:,:)               !   transmitted solar radiation  
     1985                
     1986      ENDIF 
     1987       
     1988      IF ( nice_jules == 2 ) THEN 
     1989       
     1990         ! here we must receive the qsr_ice_tr array from the coupler 
     1991         ! for now just assume zero 
     1992          
     1993         qsr_ice_tr(:,:,:) = 0.0_wp 
     1994       
     1995      ENDIF 
     1996 
     1997 
    19591998 
    19601999      CALL wrk_dealloc( jpi,jpj,     zcptn, zcptrain, zcptsnw, ziceld, zmsk, zsnw ) 
    19612000      CALL wrk_dealloc( jpi,jpj,     zemp_tot, zemp_ice, zemp_oce, ztprecip, zsprecip, zevap_oce, zevap_ice, zdevap_ice ) 
    19622001      CALL wrk_dealloc( jpi,jpj,     zqns_tot, zqns_oce, zqsr_tot, zqsr_oce, zqprec_ice, zqemp_oce, zqemp_ice ) 
    1963       CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice ) 
     2002      CALL wrk_dealloc( jpi,jpj,jpl, zqns_ice, zqsr_ice, zdqns_ice, zqevap_ice, zfrqsr_tr_i ) 
    19642003      ! 
    19652004      IF( nn_timing == 1 )   CALL timing_stop('sbc_cpl_ice_flx') 
    … …  
    20572096                   ztmp1(:,:) = SUM( alb_ice (:,:,1:jpl) * a_i(:,:,1:jpl), dim=3 ) / SUM( a_i(:,:,1:jpl), dim=3 ) 
    20582097                ELSEWHERE 
    2059                    ztmp1(:,:) = albedo_oce_mix(:,:) 
     2098                   ztmp1(:,:) = alb_oce_mix(:,:) 
    20602099                END WHERE 
    20612100             CASE default   ;   CALL ctl_stop( 'sbc_cpl_snd: wrong definition of sn_snd_alb%clcat' ) 
    … …  
    20852124 
    20862125      IF( ssnd(jps_albmix)%laction ) THEN                         ! mixed ice-ocean 
    2087          ztmp1(:,:) = albedo_oce_mix(:,:) * zfr_l(:,:) 
     2126         ztmp1(:,:) = alb_oce_mix(:,:) * zfr_l(:,:) 
    20882127         DO jl=1,jpl 
    20892128            ztmp1(:,:) = ztmp1(:,:) + alb_ice(:,:,jl) * a_i(:,:,jl) 
Note: See TracChangeset for help on using the changeset viewer.