Changeset 13553


Ignore:
Timestamp:
2020-10-01T13:33:30+02:00 (4 months ago)
Author:
hadcv
Message:

Merge in trunk up to [13550]

Location:
NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling
Files:
179 edited

Legend:

Unmodified
Added
Removed
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/1_context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/1_namelist_cfg

    r13208 r13553  
    175175!!                                                                    !! 
    176176!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    177 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    178 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     177!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     178!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    179179!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    180180!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    353353&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T) 
    354354!----------------------------------------------------------------------- 
    355       rn_eice     =   0       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4    
    356355/ 
    357356!!====================================================================== 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/2_context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/3_context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AGRIF_DEMO/EXPREF/namelist_cfg

    r13286 r13553  
    180180!!                                                                    !! 
    181181!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    182 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    183 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     182!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     183!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    184184!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    185185!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    354354&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T) 
    355355!----------------------------------------------------------------------- 
    356       rn_eice     =   0       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4    
    357356/ 
    358357!!====================================================================== 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AMM12/EXPREF/context_nemo.xml

    r12377 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/AMM12/EXPREF/namelist_cfg

    r12489 r13553  
    212212!!                                                                    !! 
    213213!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    214 !!   namdrg_top    top    friction                                      (ln_OFF =F & ln_isfcav=T) 
    215 !!   namdrg_bot    bottom friction                                      (ln_OFF =F) 
     214!!   namdrg_top    top    friction                                      (ln_drg_OFF =F & ln_isfcav=T) 
     215!!   namdrg_bot    bottom friction                                      (ln_drg_OFF =F) 
    216216!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    217217!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/C1D_PAPA/EXPREF/namelist_cfg

    r12933 r13553  
    258258!!                                                                    !! 
    259259!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    260 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    261 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     260!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     261!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    262262!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    263263!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    270270/ 
    271271!----------------------------------------------------------------------- 
    272 &namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
    273 !----------------------------------------------------------------------- 
    274 / 
    275 !----------------------------------------------------------------------- 
    276 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     272&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T) 
     273!----------------------------------------------------------------------- 
     274/ 
     275!----------------------------------------------------------------------- 
     276&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    277277!----------------------------------------------------------------------- 
    278278/ 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/GYRE_BFM/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/GYRE_BFM/EXPREF/namelist_cfg

    r12489 r13553  
    101101!!                                                                    !! 
    102102!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    103 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    104 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     103!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     104!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    105105!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    106106!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/GYRE_PISCES/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/GYRE_PISCES/EXPREF/namelist_cfg

    r12489 r13553  
    9999!!                                                                    !! 
    100100!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    101 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    102 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     101!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     102!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    103103!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    104104!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_ABL/EXPREF/namelist_cfg

    r13208 r13553  
    217217!!                                                                    !! 
    218218!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    219 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    220 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     219!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     220!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    221221!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    222222!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_PISCES/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg

    r13286 r13553  
    207207!!                                                                    !! 
    208208!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    209 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    210 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     209!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     210!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    211211!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    212212!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    378378                               !        = 2 add a tke source just at the base of the ML 
    379379                               !        = 3 as = 1 applied on HF part of the stress           (ln_cpl=T) 
    380       rn_eice     =   0       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4    
    381380/ 
    382381!----------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_OFF_PISCES/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_OFF_PISCES/EXPREF/namelist_cfg

    r12489 r13553  
    190190!!                                                                    !! 
    191191!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    192 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    193 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     192!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     193!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    194194!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    195195!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    201201/ 
    202202!----------------------------------------------------------------------- 
    203 &namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
    204 !----------------------------------------------------------------------- 
    205 / 
    206 !----------------------------------------------------------------------- 
    207 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     203&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T) 
     204!----------------------------------------------------------------------- 
     205/ 
     206!----------------------------------------------------------------------- 
     207&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    208208!----------------------------------------------------------------------- 
    209209/ 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_OFF_TRC/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_OFF_TRC/EXPREF/namelist_cfg

    r12489 r13553  
    188188!!                                                                    !! 
    189189!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    190 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    191 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     190!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     191!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    192192!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    193193!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    199199/ 
    200200!----------------------------------------------------------------------- 
    201 &namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
    202 !----------------------------------------------------------------------- 
    203 / 
    204 !----------------------------------------------------------------------- 
    205 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     201&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T) 
     202!----------------------------------------------------------------------- 
     203/ 
     204!----------------------------------------------------------------------- 
     205&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    206206!----------------------------------------------------------------------- 
    207207/ 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_SAS_ICE/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg

    r13286 r13553  
    122122!!                                                                    !! 
    123123!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    124 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    125 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     124!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     125!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    126126!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    127127!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SHARED/field_def_nemo-ice.xml

    r12377 r13553  
    4949          <field id="icehpnd"      long_name="melt pond depth"                                         standard_name="sea_ice_meltpond_depth"                    unit="m" />  
    5050          <field id="icevpnd"      long_name="melt pond volume"                                        standard_name="sea_ice_meltpond_volume"                   unit="m" />  
     51          <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" />  
     52          <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" />  
    5153      
    5254     <!-- heat --> 
     
    8183          <field id="icediv"       long_name="Divergence of the sea-ice velocity field"                standard_name="divergence_of_sea_ice_velocity"            unit="s-1"  /> 
    8284          <field id="iceshe"       long_name="Maximum shear of sea-ice velocity field"                 standard_name="maximum_shear_of_sea_ice_velocity"         unit="s-1"  /> 
    83       
     85          <field id="beta_evp"     long_name="Relaxation parameter of ice rheology (beta)"             standard_name="relaxation_parameter_of_ice_rheology"      unit=""  />    
     86  
    8487     <!-- surface heat fluxes --> 
    8588          <field id="qt_ice"       long_name="total heat flux at ice surface"                          standard_name="surface_downward_heat_flux_in_air"         unit="W/m2" /> 
     
    173176          <field id="frq_m"    unit="-"    /> 
    174177 
     178          <!-- rheology convergence tests --> 
     179          <field id="uice_cvg"   long_name="sea ice velocity convergence"      standard_name="sea_ice_velocity_convergence"      unit="m/s" /> 
     180 
    175181     <!-- ================= --> 
    176182          <!-- Add-ons for SIMIP --> 
     
    211217          <field id="dmisum"       long_name="sea-ice mass change through surface melting"             standard_name="tendency_of_sea_ice_amount_due_to_surface_melting"                       unit="kg/m2/s" /> 
    212218          <field id="dmibom"       long_name="sea-ice mass change through bottom melting"              standard_name="tendency_of_sea_ice_amount_due_to_basal_melting"                         unit="kg/m2/s" /> 
     219          <field id="dmilam"       long_name="sea-ice mass change through lateral melting"             standard_name="tendency_of_sea_ice_amount_due_to_lateral_melting"                       unit="kg/m2/s" /> 
    213220          <field id="dmsspr"       long_name="snow mass change through snow fall"                      standard_name="snowfall_flux"                                                           unit="kg/m2/s" /> 
    214221          <field id="dmsmel"       long_name="snow mass change through melt"                           standard_name="surface_snow_melt_flux"                                                  unit="kg/m2/s" /> 
     
    289296          <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        />  
    290297          <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" />  
     298          <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" />  
    291299          <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        />  
     300          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        />  
    292301          <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        /> 
    293302          <field id="iceage_cat"   long_name="Ice age per category"                              unit="days"    detect_missing_value="true" /> 
     
    300309          <field id="snwthic_cat_cmip"     long_name="Snow thickness in thickness categories"          standard_name="snow_thickness_over_categories"        detect_missing_value="true" unit="m"  > snwthic_cat      * icemask_cat + $missval * (1.-icemask_cat) </field> 
    301310          <field id="iceconc_cat_pct_cmip" long_name="Sea-ice area fractions in thickness categories"  standard_name="sea_ice_area_fraction_over_categories" detect_missing_value="true" unit="%"  > iceconc_cat*100. * icemask_cat + $missval * (1.-icemask_cat) </field> 
     311 
     312          <!-- heat diffusion convergence tests --> 
     313          <field id="tice_cvgerr" long_name="sea ice temperature convergence error"      standard_name="sea_ice_temperature_convergence_err" unit="K" /> 
     314          <field id="tice_cvgstp" long_name="sea ice temperature convergence iterations" standard_name="sea_ice_temperature_convergence_stp" unit=""  /> 
    302315 
    303316   </field_group> <!-- SBC_3D --> 
     
    560573          <field field_ref="dmisum"           name="sidmassmelttop"   /> 
    561574          <field field_ref="dmibom"           name="sidmassmeltbot"   /> 
     575          <field field_ref="dmilam"           name="sidmassmeltlat"   /> 
    562576          <field field_ref="dmsspr"           name="sndmasssnf"       /> 
    563577          <field field_ref="dmsmel"           name="sndmassmelt"      /> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SHARED/field_def_nemo-oce.xml

    r13214 r13553  
    129129        <!-- AGRIF sponge --> 
    130130        <field id="agrif_spt"         long_name=" AGRIF t-sponge coefficient"   unit=" " /> 
     131    
     132   <!-- additions to diawri.F90 --> 
     133        <field id="socegrad"    long_name="module of salinity gradient"              unit="psu/m"   grid_ref="grid_T_3D"/> 
     134        <field id="socegrad2"   long_name="square of module of salinity gradient"    unit="psu2/m2" grid_ref="grid_T_3D"/> 
     135        <field id="ke"          long_name="kinetic energy"          standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_3D" /> 
     136        <field id="ke_int"      long_name="vertical integration of kinetic energy"   unit="m3/s2"   /> 
     137        <field id="relvor"      long_name="relative vorticity"                       unit="s-1"    grid_ref="grid_T_3D"/> 
     138        <field id="absvor"      long_name="absolute vorticity"                       unit="s-1"    grid_ref="grid_T_3D"/> 
     139        <field id="potvor"      long_name="potential vorticity"                      unit="s-1"    grid_ref="grid_T_3D"/> 
     140        <field id="salt2c"      long_name="Salt content vertically integrated"       unit="1e-3*kg/m2" /> 
    131141 
    132142        <!-- t-eddy viscosity coefficients (ldfdyn) --> 
     
    177187        <field id="alpha"        long_name="thermal expansion"                                                         unit="degC-1" grid_ref="grid_T_3D" /> 
    178188        <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" /> 
    179         <field id="bn2"          long_name="squared Brunt-Vaisala frequency"                                           unit="s-1"    grid_ref="grid_T_3D" /> 
    180189        <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" /> 
    181190 
    182191        <!-- Energy - horizontal divergence --> 
    183         <field id="eken"         long_name="kinetic energy"          standard_name="specific_kinetic_energy_of_sea_water"   unit="m2/s2"  grid_ref="grid_T_3D" /> 
    184192        <field id="hdiv"         long_name="horizontal divergence"                                                          unit="s-1"    grid_ref="grid_T_3D" /> 
    185193 
     
    499507        <field id="uocetr_vsum_op"    long_name="ocean current along i-axis * e3u * e2u summed on the vertical"  read_access="true"  freq_op="1mo"    field_ref="e2u"       unit="m3/s"> @uocetr_vsum </field> 
    500508        <field id="uocetr_vsum_cumul" long_name="ocean current along i-axis * e3u * e2u cumulated from southwest point" freq_offset="_reset_" operation="instant" freq_op="1mo"  unit="m3/s" /> 
    501         <field id="msftbarot"         long_name="ocean_barotropic_mass_streamfunction"   unit="kg s-1" > uocetr_vsum_cumul * $rau0 </field> 
     509        <field id="msftbarot"         long_name="ocean_barotropic_mass_streamfunction"   unit="kg s-1" > uocetr_vsum_cumul * $rho0 </field> 
    502510 
    503511 
     
    655663        <field id="w_masstr2"    long_name="square of vertical mass transport"              standard_name="square_of_upward_ocean_mass_transport"   unit="kg2/s2" /> 
    656664 
     665        <!-- EOS --> 
     666        <field id="bn2"          long_name="squared Brunt-Vaisala frequency"                unit="s-2" /> 
     667 
    657668      </field_group> 
    658669 
     
    700711         <field id="uocetr_vsum_section"  long_name="Total 2D transport in i-direction"               field_ref="uoce_e3u_ave_vsum"    grid_ref="grid_U_scalar"  detect_missing_value="true"> this * e2u </field> 
    701712         <field id="uocetr_strait"        long_name="Total transport across lines in i-direction"     field_ref="uocetr_vsum_section"  grid_ref="grid_U_4strait" /> 
    702          <field id="u_masstr_strait"      long_name="Sea water transport across line in i-direction"  field_ref="uocetr_strait"        grid_ref="grid_U_4strait_hsum" unit="kg/s"> this * maskMFO_u * $rau0 </field> 
     713         <field id="u_masstr_strait"      long_name="Sea water transport across line in i-direction"  field_ref="uocetr_strait"        grid_ref="grid_U_4strait_hsum" unit="kg/s"> this * maskMFO_u * $rho0 </field> 
    703714 
    704715         <field id="voce_e3v_ave"         long_name="Monthly average of v*e3v"                        field_ref="voce_e3v"                    freq_op="1mo"   freq_offset="_reset_" > @voce_e3v </field> 
     
    706717         <field id="vocetr_vsum_section"  long_name="Total 2D transport of in j-direction"            field_ref="voce_e3v_ave_vsum"    grid_ref="grid_V_scalar"  detect_missing_value="true"> this * e1v </field> 
    707718         <field id="vocetr_strait"        long_name="Total transport across lines in j-direction"     field_ref="vocetr_vsum_section"  grid_ref="grid_V_4strait"  /> 
    708          <field id="v_masstr_strait"      long_name="Sea water transport across line in j-direction"  field_ref="vocetr_strait"        grid_ref="grid_V_4strait_hsum" unit="kg/s"> this * maskMFO_v * $rau0 </field> 
     719         <field id="v_masstr_strait"      long_name="Sea water transport across line in j-direction"  field_ref="vocetr_strait"        grid_ref="grid_V_4strait_hsum" unit="kg/s"> this * maskMFO_v * $rho0 </field> 
    709720 
    710721         <field id="masstr_strait"        long_name="Sea water transport across line"                                                  grid_ref="grid_4strait"  > u_masstr_strait + v_masstr_strait </field> 
    711722      </field_group> 
    712  
    713723 
    714724      <!-- variables available with ln_floats --> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SHARED/namelist_ice_ref

    r12377 r13553  
    4343   ln_cat_usr       = .false.         !  ice categories are defined by rn_catbnd below (m) 
    4444      rn_catbnd     =   0.,0.45,1.1,2.1,3.7,6.0   
    45    rn_himin         =   0.1           !  minimum ice thickness (m) used in remapping 
     45   rn_himin         =   0.1           !  minimum ice thickness (m) allowed 
     46   rn_himax         =  99.0           !  maximum ice thickness (m) allowed 
    4647/ 
    4748!------------------------------------------------------------------------------ 
     
    5657   rn_ishlat        =   2.            !  lbc : free slip (0) ; partial slip (0-2) ; no slip (2) ; strong slip (>2) 
    5758   ln_landfast_L16  = .false.         !  landfast: parameterization from Lemieux 2016 
    58       rn_depfra     =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
     59      rn_lf_depfra  =   0.125         !        fraction of ocean depth that ice must reach to initiate landfast 
    5960                                      !          recommended range: [0.1 ; 0.25] 
    60       rn_icebfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
    61       rn_lfrelax    =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
    62       rn_tensile    =   0.05          !        isotropic tensile strength [0-0.5??] 
     61      rn_lf_bfr     =  15.            !        maximum bottom stress per unit volume [N/m3] 
     62      rn_lf_relax   =   1.e-5         !        relaxation time scale to reach static friction [s-1] 
     63      rn_lf_tensile =   0.05          !        isotropic tensile strength [0-0.5??] 
    6364/ 
    6465!------------------------------------------------------------------------------ 
     
    9192!------------------------------------------------------------------------------ 
    9293   ln_rhg_EVP       = .true.          !  EVP rheology 
    93       ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
     94      ln_aEVP       = .true.          !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
    9495      rn_creepl     =   2.0e-9        !     creep limit [1/s] 
    9596      rn_ecc        =   2.0           !     eccentricity of the elliptical yield curve           
    96       nn_nevp       = 120             !     number of EVP subcycles                              
     97      nn_nevp       = 100             !     number of EVP subcycles                              
    9798      rn_relast     =   0.333         !     ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast  
    98                                       !        advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 
     99                                      !        advised value: 1/3 (nn_nevp=100) or 1/9 (nn_nevp=300) 
     100   nn_rhg_chkcvg    =   0             !  check convergence of rheology 
     101                                      !     = 0  no check 
     102                                      !     = 1  check at the main time step (output xml: uice_cvg) 
     103                                      !     = 2  check at both main and rheology time steps (additional output: ice_cvg.nc) 
     104                                      !          this option 2 asks a lot of communications between cpu 
    99105/ 
    100106!------------------------------------------------------------------------------ 
    101107&namdyn_adv     !   Ice advection 
    102108!------------------------------------------------------------------------------ 
    103    ln_adv_Pra       = .true.         !  Advection scheme (Prather) 
    104    ln_adv_UMx       = .false.          !  Advection scheme (Ultimate-Macho) 
     109   ln_adv_Pra       = .true.          !  Advection scheme (Prather) 
     110   ln_adv_UMx       = .false.         !  Advection scheme (Ultimate-Macho) 
    105111      nn_UMx        =   5             !     order of the scheme for UMx (1-5 ; 20=centered 2nd order) 
    106112/ 
     
    109115!------------------------------------------------------------------------------ 
    110116   rn_cio           =   5.0e-03       !  ice-ocean drag coefficient (-) 
    111    rn_blow_s        =   0.66          !  mesure of snow blowing into the leads 
     117   nn_snwfra        =   2             !  calculate the fraction of ice covered by snow (for zdf and albedo) 
     118                                      !     = 0  fraction = 1 (if snow) or 0 (if no snow) 
     119                                      !     = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation] 
     120                                      !     = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation] 
     121   rn_snwblow       =   0.66          !  mesure of snow blowing into the leads 
    112122                                      !     = 1 => no snow blowing, < 1 => some snow blowing 
    113123   nn_flxdist       =  -1             !  Redistribute heat flux over ice categories 
     
    118128   ln_cndflx        = .false.         !  Use conduction flux as surface boundary conditions (i.e. for Jules coupling) 
    119129      ln_cndemulate = .false.         !     emulate conduction flux (if not provided in the inputs) 
     130   nn_qtrice        =   1             !  Solar flux transmitted thru the surface scattering layer: 
     131                                      !     = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow)  
     132                                      !     = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities) 
    120133/ 
    121134!------------------------------------------------------------------------------ 
     
    126139   ln_icedO         = .true.          !  activate ice growth in open-water (T) or not (F) 
    127140   ln_icedS         = .true.          !  activate brine drainage (T) or not (F) 
     141   ! 
     142   ln_leadhfx       = .true.          !  heat in the leads is used to melt sea-ice before warming the ocean 
    128143/ 
    129144!------------------------------------------------------------------------------ 
     
    135150   rn_cnd_s         =   0.31          !  thermal conductivity of the snow (0.31 W/m/K, Maykut and Untersteiner, 1971) 
    136151                                      !     Obs: 0.1-0.5 (Lecomte et al, JAMES 2013) 
    137    rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice [1/m] 
     152   rn_kappa_i       =   1.0           !  radiation attenuation coefficient in sea ice                     [1/m] 
     153   rn_kappa_s       =  10.0           !  nn_qtrice = 0: radiation attenuation coefficient in snow         [1/m] 
     154   rn_kappa_smlt    =   7.0           !  nn_qtrice = 1: radiation attenuation coefficient in melting snow [1/m] 
     155   rn_kappa_sdry    =  10.0           !                 radiation attenuation coefficient in dry snow     [1/m] 
     156   ln_zdf_chkcvg    = .false.         !  check convergence of heat diffusion scheme (outputs: tice_cvgerr, tice_cvgstp) 
    138157/ 
    139158!------------------------------------------------------------------------------ 
     
    175194&namthd_pnd     !   Melt ponds 
    176195!------------------------------------------------------------------------------ 
    177    ln_pnd           = .false.         !  activate melt ponds or not 
    178      ln_pnd_H12     = .false.         !  activate evolutive melt ponds (from Holland et al 2012) 
    179      ln_pnd_CST     = .false.         !  activate constant  melt ponds 
    180        rn_apnd      =   0.2           !     prescribed pond fraction, at Tsu=0 degC 
    181        rn_hpnd      =   0.05          !     prescribed pond depth,    at Tsu=0 degC 
    182      ln_pnd_alb     = .false.         !  melt ponds affect albedo or not 
     196   ln_pnd            = .true.         !  activate melt ponds or not 
     197      ln_pnd_LEV     = .true.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 
     198         rn_apnd_min =   0.15         !     minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 
     199         rn_apnd_max =   0.85         !     maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 
     200      ln_pnd_CST     = .false.        !  constant  melt ponds 
     201         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC 
     202         rn_hpnd     =   0.05         !     prescribed pond depth,    at Tsu=0 degC 
     203      ln_pnd_lids    = .true.         !  frozen lids on top of the ponds (only for ln_pnd_LEV) 
     204      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo 
    183205/ 
    184206!------------------------------------------------------------------------------ 
     
    186208!------------------------------------------------------------------------------ 
    187209   ln_iceini        = .true.          !  activate ice initialization (T) or not (F) 
    188    ln_iceini_file   = .false.         !  netcdf file provided for initialization (T) or not (F) 
     210   nn_iceini_file   =   0             !     0 = Initialise sea ice based on SSTs 
     211                                      !     1 = Initialise sea ice from single category netcdf file 
     212                                      !     2 = Initialise sea ice from multi category restart file 
    189213   rn_thres_sst     =   2.0           !  max temp. above Tfreeze with initial ice = (sst - tfreeze) 
    190214   rn_hti_ini_n     =   3.0           !  initial ice thickness       (m), North 
     
    206230   rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North 
    207231   rn_hpd_ini_s     =   0.05          !        "            "             South 
    208    ! -- for ln_iceini_file = T 
    209    sn_hti = 'Ice_initialization'    , -12 ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    210    sn_hts = 'Ice_initialization'    , -12 ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    211    sn_ati = 'Ice_initialization'    , -12 ,'ati'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    212    sn_smi = 'Ice_initialization'    , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    213    sn_tmi = 'Ice_initialization'    , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    214    sn_tsu = 'Ice_initialization'    , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    215    sn_tms = 'NOT USED'              , -12 ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     232   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North 
     233   rn_hld_ini_s     =   0.0           !        "            "             South 
     234   ! -- for nn_iceini_file = 1 
     235   sn_hti = 'Ice_initialization'    , -12. ,'hti'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     236   sn_hts = 'Ice_initialization'    , -12. ,'hts'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     237   sn_ati = 'Ice_initialization'    , -12. ,'ati'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     238   sn_smi = 'Ice_initialization'    , -12. ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     239   sn_tmi = 'Ice_initialization'    , -12. ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     240   sn_tsu = 'Ice_initialization'    , -12. ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     241   sn_tms = 'NOT USED'              , -12. ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    216242   !      melt ponds (be careful, sn_apd is the pond concentration (not fraction), so it differs from rn_apd) 
    217    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    218    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     243   sn_apd = 'NOT USED'              , -12. ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     244   sn_hpd = 'NOT USED'              , -12. ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     245   sn_hld = 'NOT USED'              , -12. ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    219246   cn_dir='./' 
    220247/ 
     
    238265   ln_icediahsb     = .false.         !  output the heat, mass & salt budgets (T) or not (F) 
    239266   ln_icectl        = .false.         !  ice points output for debug (T or F) 
    240    iiceprt          =  10             !  i-index for debug 
    241    jiceprt          =  10             !  j-index for debug 
    242 / 
     267      iiceprt       =  10             !     i-index for debug 
     268      jiceprt       =  10             !     j-index for debug 
     269/ 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SHARED/namelist_ref

    r13514 r13553  
    303303   sn_uoatm    = 'NOT USED'                   ,    6.        , 'UOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Uoceatm', '' 
    304304   sn_voatm    = 'NOT USED'                   ,    6.        , 'VOATM'   ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , 'Voceatm', '' 
     305   sn_cc       = 'NOT USED'                   ,   24.        , 'CC'      ,   .false.   , .true. , 'yearly'  , 'weights_core_orca2_bilinear_noc.nc' , ''       , '' 
    305306   sn_hpgi     = 'NOT USED'                   ,   24.        , 'uhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'UG'     , '' 
    306307   sn_hpgj     = 'NOT USED'                   ,   24.        , 'vhpg'    ,   .false.   , .false., 'monthly' , 'weights_ERAI3D_F128_2_ORCA2_bicubic', 'VG'     , '' 
     
    342343&namsbc_cpl    !   coupled ocean/atmosphere model                       ("key_oasis3") 
    343344!----------------------------------------------------------------------- 
    344    nn_cplmodel   =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data 
    345    ln_usecplmask = .false. !  use a coupling mask file to merge data received from several models 
    346    !                       !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
    347    nn_cats_cpl   =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1) 
     345   nn_cplmodel       =     1   !  Maximum number of models to/from which NEMO is potentially sending/receiving data 
     346   ln_usecplmask     = .false. !  use a coupling mask file to merge data received from several models 
     347   !                           !   -> file cplmask.nc with the float variable called cplmask (jpi,jpj,nn_cplmodel) 
     348   ln_scale_ice_flux = .false. !  use ice fluxes that are already "ice weighted" ( i.e. multiplied ice concentration) 
     349   nn_cats_cpl       =     5   !  Number of sea ice categories over which coupling is to be carried out (if not 1) 
    348350   !_____________!__________________________!____________!_____________!______________________!________! 
    349351   !             !        description       !  multiple  !    vector   !       vector         ! vector ! 
     
    551553         !           !  file name  ! frequency (hours) ! variable  ! time interp.!  clim   ! 'yearly'/ ! weights  ! rotation ! land/sea mask ! 
    552554         !           !             !  (if <0  months)  !   name    !  (logical)  !  (T/F)  ! 'monthly' ! filename ! pairing  ! filename      ! 
    553          sn_isfpar_zmax = 'isfmlt_par',       0        ,'sozisfmax',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    554          sn_isfpar_zmin = 'isfmlt_par',       0        ,'sozisfmin',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     555         sn_isfpar_zmax = 'isfmlt_par',       0.       ,'sozisfmax',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
     556         sn_isfpar_zmin = 'isfmlt_par',       0.       ,'sozisfmin',  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    555557         !* 'spe' and 'oasis' case 
    556          sn_isfpar_fwf = 'isfmlt_par' ,      -12.      ,'sofwfisf' ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    '' 
     558         sn_isfpar_fwf = 'isfmlt_par' ,      -12.      ,'sofwfisf' ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    557559         !* 'bg03' case 
    558          sn_isfpar_Leff = 'isfmlt_par',       0.       ,'Leff'     ,  .false.    , .true.  , 'yearly'   ,    ''    ,   ''     ,    '' 
     560         sn_isfpar_Leff = 'isfmlt_par',       0.       ,'Leff'     ,  .false.    , .true.  , 'yearly'  ,    ''    ,   ''     ,    '' 
    559561      ! 
    560562      ! ---------------- ice sheet coupling ------------------------------- 
     
    739741   bn_aip      = 'NOT USED'              ,         24.       , 'siapnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    740742   bn_hip      = 'NOT USED'              ,         24.       , 'sihpnd'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
     743   bn_hil      = 'NOT USED'              ,         24.       , 'sihlid'  ,    .true.   , .false.,  'daily'  ,    ''            ,   ''     ,     '' 
    741744   ! if bn_t_i etc are "not used", then define arbitrary temperatures and salinity and ponds 
    742745   rn_ice_tem  = 270.         !  arbitrary temperature               of incoming sea ice 
     
    745748   rn_ice_apnd = 0.2          !       --   pond fraction = a_ip/a_i            -- 
    746749   rn_ice_hpnd = 0.05         !       --   pond depth                          -- 
     750   rn_ice_hlid = 0.0          !       --   pond lid depth                      -- 
    747751/ 
    748752!----------------------------------------------------------------------- 
     
    757761!!                                                                    !! 
    758762!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    759 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    760 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     763!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     764!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    761765!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    762766!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    766770&namdrg        !   top/bottom drag coefficient                          (default: NO selection) 
    767771!----------------------------------------------------------------------- 
    768    ln_OFF      = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
     772   ln_drg_OFF  = .false.   !  free-slip       : Cd = 0                  (F => fill namdrg_bot 
    769773   ln_lin      = .false.   !      linear  drag: Cd = Cd0 Uc0                   &   namdrg_top) 
    770774   ln_non_lin  = .false.   !  non-linear  drag: Cd = Cd0 |U| 
     
    772776   ! 
    773777   ln_drgimp   = .true.    !  implicit top/bottom friction flag 
    774 / 
    775 !----------------------------------------------------------------------- 
    776 &namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
     778      ln_drgice_imp = .true. ! implicit ice-ocean drag 
     779/ 
     780!----------------------------------------------------------------------- 
     781&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T) 
    777782!----------------------------------------------------------------------- 
    778783   rn_Cd0      =  1.e-3    !  drag coefficient [-] 
     
    785790/ 
    786791!----------------------------------------------------------------------- 
    787 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     792&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    788793!----------------------------------------------------------------------- 
    789794   rn_Cd0      =  1.e-3    !  drag coefficient [-] 
     
    838843                                 ! 
    839844   !                     ! S-EOS coefficients (ln_seos=T): 
    840    !                             !  rd(T,S,Z)*rau0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 
     845   !                             !  rd(T,S,Z)*rho0 = -a0*(1+.5*lambda*dT+mu*Z+nu*dS)*dT+b0*dS 
    841846   rn_a0       =  1.6550e-1      !  thermal expension coefficient 
    842847   rn_b0       =  7.6554e-1      !  saline  expension coefficient 
     
    11421147   rn_bshear   =   1.e-20  ! background shear (>0) currently a numerical threshold (do not change it) 
    11431148   nn_pdl      =   1       !  Prandtl number function of richarson number (=1, avt=pdl(Ri)*avm) or not (=0, avt=avm) 
    1144    nn_mxl      =   2       !  mixing length: = 0 bounded by the distance to surface and bottom 
     1149   nn_mxl      =   3       !  mixing length: = 0 bounded by the distance to surface and bottom 
    11451150   !                       !                 = 1 bounded by the local vertical scale factor 
    11461151   !                       !                 = 2 first vertical derivative of mixing length bounded by 1 
    11471152   !                       !                 = 3 as =2 with distinct dissipative an mixing length scale 
    11481153   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F) 
    1149       nn_mxlice    = 0        ! type of scaling under sea-ice 
     1154      nn_mxlice    = 2        ! type of scaling under sea-ice 
    11501155                              !    = 0 no scaling under sea-ice 
    11511156                              !    = 1 scaling with constant sea-ice thickness 
    1152                               !    = 2  scaling with mean sea-ice thickness ( only with SI3 sea-ice model ) 
    1153                               !    = 3  scaling with maximum sea-ice thickness 
     1157                              !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model ) 
     1158                              !    = 3 scaling with maximum sea-ice thickness 
    11541159      rn_mxlice   = 10.       ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1) 
    11551160   rn_mxl0     =   0.04    !  surface  buoyancy lenght scale minimum value 
    1156    ln_drg      = .false.   !  top/bottom friction added as boundary condition of TKE 
    11571161   ln_lc       = .true.    !  Langmuir cell parameterisation (Axell 2002) 
    11581162      rn_lc       =   0.15    !  coef. associated to Langmuir cells 
     
    11651169                              !        = 0  constant 10 m length scale 
    11661170                              !        = 1  0.5m at the equator to 30m poleward of 40 degrees 
    1167       rn_eice     =   4       !  below sea ice: =0 ON ; =4 OFF when ice fraction > 1/4 
     1171   nn_eice     =   1       !  attenutaion of langmuir & surface wave breaking under ice 
     1172   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking 
     1173   !                       !           = 1 weigthed by 1-TANH(10*fr_i) 
     1174   !                       !           = 2 weighted by 1-fr_i 
     1175   !                       !           = 3 weighted by 1-MIN(1,4*fr_i)    
    11681176/ 
    11691177!----------------------------------------------------------------------- 
     
    11781186   rn_charn      = 70000.  !  Charnock constant for wb induced roughness length 
    11791187   rn_hsro       =  0.02   !  Minimum surface roughness 
     1188   rn_hsri       =  0.03   !  Ice-ocean roughness 
    11801189   rn_frac_hs    =   1.3   !  Fraction of wave height as roughness (if nn_z0_met>1) 
    11811190   nn_z0_met     =     2   !  Method for surface roughness computation (0/1/2/3) 
    1182    !                             ! =3 requires ln_wave=T 
     1191   !                       !     = 3 requires ln_wave=T 
     1192   nn_z0_ice     =   1     !  attenutaion of surface wave breaking under ice 
     1193   !                       !           = 0 no impact of ice cover 
     1194   !                       !           = 1 roughness uses rn_hsri and is weigthed by 1-TANH(10*fr_i) 
     1195   !                       !           = 2 roughness uses rn_hsri and is weighted by 1-fr_i 
     1196   !                       !           = 3 roughness uses rn_hsri and is weighted by 1-MIN(1,4*fr_i) 
    11831197   nn_bc_surf    =     1   !  surface condition (0/1=Dir/Neum) 
    11841198   nn_bc_bot     =     1   !  bottom condition (0/1=Dir/Neum) 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SPITZ12/EXPREF/context_nemo.xml

    r12276 r13553  
    1111       <variable id="ref_month" type="int"> 01 </variable> 
    1212       <variable id="ref_day"   type="int"> 01 </variable> 
    13        <variable id="rau0"      type="float" > 1026.0 </variable> 
     13       <variable id="rho0"      type="float" > 1026.0 </variable> 
    1414       <variable id="cpocean"   type="float" > 3991.86795711963 </variable> 
    1515       <variable id="convSpsu"  type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SPITZ12/EXPREF/namelist_cfg

    r12489 r13553  
    205205!!                                                                    !! 
    206206!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    207 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    208 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     207!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     208!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    209209!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    210210!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    216216   ln_loglayer = .true.   !  logarithmic drag: Cd = vkarmn/log(z/z0) |U| 
    217217   ln_drgimp   = .true.   !  implicit top/bottom friction flag 
    218 / 
    219 !----------------------------------------------------------------------- 
    220 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     218      ln_drgice_imp = .true. ! implicit ice-ocean drag 
     219/ 
     220!----------------------------------------------------------------------- 
     221&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    221222!----------------------------------------------------------------------- 
    222223   rn_Cd0      =  2.5e-3   !  drag coefficient [-] 
     
    339340   nn_havtb    =    1         !  horizontal shape for avtb (=1) or not (=0) 
    340341/ 
     342!----------------------------------------------------------------------- 
     343&namzdf_tke    !   turbulent eddy kinetic dependent vertical diffusion  (ln_zdftke =T) 
     344!----------------------------------------------------------------------- 
     345   ln_mxl0     = .true.    !  surface mixing length scale = F(wind stress) (T) or not (F) 
     346      nn_mxlice    = 0        ! type of scaling under sea-ice 
     347                              !    = 0 no scaling under sea-ice 
     348                              !    = 1 scaling with constant sea-ice thickness 
     349                              !    = 2 scaling with mean sea-ice thickness ( only with SI3 sea-ice model ) 
     350                              !    = 3 scaling with maximum sea-ice thickness 
     351   nn_eice     =   0       !  attenutaion of langmuir & surface wave breaking under ice 
     352   !                       !           = 0 no impact of ice cover on langmuir & surface wave breaking 
     353   !                       !           = 1 weigthed by 1-TANH(10*fr_i) 
     354   !                       !           = 2 weighted by 1-fr_i 
     355   !                       !           = 3 weighted by 1-MIN(1,4*fr_i) 
     356/ 
    341357!!====================================================================== 
    342358!!                  ***  Diagnostics namelists  ***                   !! 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/SPITZ12/EXPREF/namelist_ice_cfg

    r11731 r13553  
    5555&namsbc         !   Ice surface boundary conditions 
    5656!------------------------------------------------------------------------------ 
     57   nn_snwfra        =   0             !  calculate the fraction of ice covered by snow (for zdf and albedo) 
     58                                      !     = 0  fraction = 1 (if snow) or 0 (if no snow) 
     59                                      !     = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation] 
     60                                      !     = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation] 
     61   nn_qtrice        =   0             !  Solar flux transmitted thru the surface scattering layer: 
     62                                      !     = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow) 
     63                                      !     = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities) 
    5764/ 
    5865!------------------------------------------------------------------------------ 
     
    8188&namthd_pnd     !   Melt ponds 
    8289!------------------------------------------------------------------------------ 
    83    ln_pnd           = .true.          !  activate melt ponds or not 
    84      ln_pnd_H12     = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
    85      ln_pnd_alb     = .true.          !  melt ponds affect albedo or not 
     90   ln_pnd           = .false.          !  activate melt ponds or not 
     91     ln_pnd_LEV     = .false.          !  activate level ice melt ponds 
    8692/ 
    8793 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/WED025/EXPREF/context_nemo.xml

    r11487 r13553  
    99    <!-- Year of time origin for NetCDF files; defaults to 1800 --> 
    1010       <variable id="ref_year" type="int"   > 1800 </variable> 
    11        <variable id="rau0"     type="float" > 1026.0 </variable> 
     11       <variable id="rho0"     type="float" > 1026.0 </variable> 
    1212       <variable id="cpocean"  type="float" > 3991.86795711963 </variable> 
    1313       <variable id="convSpsu" type="float" > 0.99530670233846  </variable> 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/WED025/EXPREF/namelist_cfg

    r13208 r13553  
    362362!!                                                                    !! 
    363363!!   namdrg        top/bottom drag coefficient                          (default: NO selection) 
    364 !!   namdrg_top    top    friction                                      (ln_OFF=F & ln_isfcav=T) 
    365 !!   namdrg_bot    bottom friction                                      (ln_OFF=F) 
     364!!   namdrg_top    top    friction                                      (ln_drg_OFF=F & ln_isfcav=T) 
     365!!   namdrg_bot    bottom friction                                      (ln_drg_OFF=F) 
    366366!!   nambbc        bottom temperature boundary condition                (default: OFF) 
    367367!!   nambbl        bottom boundary layer scheme                         (default: OFF) 
     
    374374/ 
    375375!----------------------------------------------------------------------- 
    376 &namdrg_top    !   TOP friction                                         (ln_OFF =F & ln_isfcav=T) 
     376&namdrg_top    !   TOP friction                                         (ln_drg_OFF =F & ln_isfcav=T) 
    377377!----------------------------------------------------------------------- 
    378378   rn_Cd0      =  2.5e-3    !  drag coefficient [-] 
    379379/ 
    380380!----------------------------------------------------------------------- 
    381 &namdrg_bot    !   BOTTOM friction                                      (ln_OFF =F) 
     381&namdrg_bot    !   BOTTOM friction                                      (ln_drg_OFF =F) 
    382382!----------------------------------------------------------------------- 
    383383   rn_Cd0      =  2.5e-3    !  drag coefficient [-] 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/cfgs/WED025/EXPREF/namelist_ice_cfg

    r12905 r13553  
    4242&namdyn_rhg     !   Ice rheology 
    4343!------------------------------------------------------------------------------ 
     44   ln_rhg_EVP       = .true.          !  EVP rheology 
     45      ln_aEVP       = .false.         !     adaptive rheology (Kimmritz et al. 2016 & 2017) 
    4446/ 
    4547!------------------------------------------------------------------------------ 
     
    5355&namsbc         !   Ice surface boundary conditions 
    5456!------------------------------------------------------------------------------ 
     57   nn_snwfra        =   0             !  calculate the fraction of ice covered by snow (for zdf and albedo) 
     58                                      !     = 0  fraction = 1 (if snow) or 0 (if no snow) 
     59                                      !     = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation] 
     60                                      !     = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation] 
     61   nn_qtrice        =   0             !  Solar flux transmitted thru the surface scattering layer: 
     62                                      !     = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow) 
     63                                      !     = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities) 
    5564/ 
    5665!------------------------------------------------------------------------------ 
     
    7988&namthd_pnd     !   Melt ponds 
    8089!------------------------------------------------------------------------------ 
    81    ln_pnd           = .true.          !  activate melt ponds or not 
    82      ln_pnd_H12     = .true.          !  activate evolutive melt ponds (from Holland et al 2012) 
    83      ln_pnd_alb     = .true.          !  melt ponds affect albedo or not 
     90   ln_pnd            = .false.         !  activate melt ponds or not 
     91      ln_pnd_LEV     = .false.         !  level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 
    8492/ 
    85  
    8693!------------------------------------------------------------------------------ 
    8794&namini         !   Ice initialization 
    8895!------------------------------------------------------------------------------ 
    8996   ln_iceini        = .true.          !  activate ice initialization (T) or not (F) 
    90    ln_iceini_file   = .true.          !  netcdf file provided for initialization (T) or not (F) 
     97   nn_iceini_file   =   1             !     0 = Initialise sea ice based on SSTs 
     98                                      !     1 = Initialise sea ice from single category netcdf file 
     99                                      !     2 = Initialise sea ice from multi category restart file 
    91100   ! -- for ln_iceini_file = T 
    92    sn_hti = 'WED025_init_JRA_200001.nc', -12 ,'icethic_cea',  .false.  , .true., 'yearly'  , '' , '', '' 
    93    sn_hts = 'WED025_init_JRA_200001.nc', -12 ,'icesnow_cea',  .false.  , .true., 'yearly'  , '' , '', '' 
    94    sn_ati = 'WED025_init_JRA_200001.nc', -12 ,'ice_cover'  ,  .false.  , .true., 'yearly'  , '' , '', '' 
    95    sn_smi = 'NOT USED'              , -12 ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    96    sn_tmi = 'NOT USED'              , -12 ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    97    sn_tsu = 'NOT USED'              , -12 ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    98    sn_tms = 'NOT USED'              , -12 ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     101   sn_hti = 'WED025_init_JRA_200001.nc', -12. ,'icethic_cea',  .false.  , .true., 'yearly'  , '' , '', '' 
     102   sn_hts = 'WED025_init_JRA_200001.nc', -12. ,'icesnow_cea',  .false.  , .true., 'yearly'  , '' , '', '' 
     103   sn_ati = 'WED025_init_JRA_200001.nc', -12. ,'ice_cover'  ,  .false.  , .true., 'yearly'  , '' , '', '' 
     104   sn_smi = 'NOT USED'              , -12. ,'smi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     105   sn_tmi = 'NOT USED'              , -12. ,'tmi'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     106   sn_tsu = 'NOT USED'              , -12. ,'tsu'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     107   sn_tms = 'NOT USED'              , -12. ,'tms'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    99108   !      melt ponds (be careful, sn_apd is the pond concentration (not fraction), so it differs from rn_apd) 
    100    sn_apd = 'NOT USED'              , -12 ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    101    sn_hpd = 'NOT USED'              , -12 ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     109   sn_apd = 'NOT USED'              , -12. ,'apd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     110   sn_hpd = 'NOT USED'              , -12. ,'hpd'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
     111   sn_hld = 'NOT USED'              , -12. ,'hld'   ,  .false.  , .true., 'yearly'  , '' , '', '' 
    102112   cn_dir='./' 
    103113/ 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/ice.F90

    r12489 r13553  
    7070   !! a_ip        |      -      |    Ice pond concentration       |       | 
    7171   !! v_ip        |      -      |    Ice pond volume per unit area| m     | 
     72   !! v_il        |    v_il_1d  |    Ice pond lid volume per area | m     | 
    7273   !!                                                                     | 
    7374   !!-------------|-------------|---------------------------------|-------| 
     
    8586   !! t_su        ! t_su_1d     |    Sea ice surface temperature  ! K     | 
    8687   !! h_ip        | h_ip_1d     |    Ice pond thickness           | m     | 
     88   !! h_il        | h_il_1d     |    Ice pond lid thickness       | m     | 
    8789   !!                                                                     | 
    8890   !! notes: the ice model only sees a bulk (i.e., vertically averaged)   | 
     
    112114   !! hm_ip       |      -      |    Mean ice pond depth          | m     | 
    113115   !! vt_ip       |      -      |    Total ice pond vol. per unit area| m | 
     116   !! hm_il       |      -      |    Mean ice pond lid depth      | m     | 
     117   !! vt_il       |      -      |    Total ice pond lid vol. per area | m | 
    114118   !!===================================================================== 
    115119 
     
    137141   REAL(wp), PUBLIC ::   rn_ishlat        !: lateral boundary condition for sea-ice 
    138142   LOGICAL , PUBLIC ::   ln_landfast_L16  !: landfast ice parameterizationfrom lemieux2016  
    139    REAL(wp), PUBLIC ::   rn_depfra        !:    fraction of ocean depth that ice must reach to initiate landfast ice 
    140    REAL(wp), PUBLIC ::   rn_icebfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
    141    REAL(wp), PUBLIC ::   rn_lfrelax       !:    relaxation time scale (s-1) to reach static friction 
    142    REAL(wp), PUBLIC ::   rn_tensile       !:    isotropic tensile strength 
     143   REAL(wp), PUBLIC ::   rn_lf_depfra     !:    fraction of ocean depth that ice must reach to initiate landfast ice 
     144   REAL(wp), PUBLIC ::   rn_lf_bfr        !:    maximum bottom stress per unit area of contact (lemieux2016) or per unit volume (home)  
     145   REAL(wp), PUBLIC ::   rn_lf_relax      !:    relaxation time scale (s-1) to reach static friction 
     146   REAL(wp), PUBLIC ::   rn_lf_tensile    !:    isotropic tensile strength 
    143147   ! 
    144148   !                                     !!** ice-ridging/rafting namelist (namdyn_rdgrft) ** 
     
    151155   INTEGER , PUBLIC ::   nn_nevp          !: number of iterations for subcycling 
    152156   REAL(wp), PUBLIC ::   rn_relast        !: ratio => telast/rDt_ice (1/3 or 1/9 depending on nb of subcycling nevp)  
     157   INTEGER , PUBLIC ::   nn_rhg_chkcvg    !: check ice rheology convergence  
    153158   ! 
    154159   !                                     !!** ice-advection namelist (namdyn_adv) ** 
     
    158163   !                                     !!** ice-surface boundary conditions namelist (namsbc) ** 
    159164                                          ! -- icethd_dh -- ! 
    160    REAL(wp), PUBLIC ::   rn_blow_s        !: coef. for partitioning of snowfall between leads and sea ice 
     165   REAL(wp), PUBLIC ::   rn_snwblow       !: coef. for partitioning of snowfall between leads and sea ice 
     166                                          ! -- icethd_zdf and icealb -- ! 
     167   INTEGER , PUBLIC ::   nn_snwfra        !: calculate the fraction of ice covered by snow 
     168   !                                      !   = 0  fraction = 1 (if snow) or 0 (if no snow) 
     169   !                                      !   = 1  fraction = 1-exp(-0.2*rhos*hsnw) [MetO formulation] 
     170   !                                      !   = 2  fraction = hsnw / (hsnw+0.02)    [CICE formulation] 
    161171                                          ! -- icethd -- ! 
    162172   REAL(wp), PUBLIC ::   rn_cio           !: drag coefficient for oceanic stress 
     
    166176   !                                      !   = 1  Average N(cat) fluxes then redistribute over the N(cat) ice using T-ice and albedo sensitivity 
    167177   !                                      !   = 2  Redistribute a single flux over categories 
     178                                          ! -- icethd_zdf -- ! 
    168179   LOGICAL , PUBLIC ::   ln_cndflx        !: use conduction flux as surface boundary condition (instead of qsr and qns)  
    169180   LOGICAL , PUBLIC ::   ln_cndemulate    !: emulate conduction flux (if not provided)  
     
    172183   INTEGER, PUBLIC, PARAMETER ::   np_cnd_ON  = 1  !: forcing from conduction flux (SM0L) (compute qcn and qsr_tr via sbcblk.F90 or sbccpl.F90) 
    173184   INTEGER, PUBLIC, PARAMETER ::   np_cnd_EMU = 2  !: emulate conduction flux via icethd_zdf.F90 (BL99) (1st round compute qcn and qsr_tr, 2nd round use it) 
    174  
     185   INTEGER, PUBLIC ::   nn_qtrice         !: Solar flux transmitted thru the surface scattering layer: 
     186   !                                      !   = 0  Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow)  
     187   !                                      !   = 1  Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities) 
     188   ! 
    175189   !                                     !!** ice-vertical diffusion namelist (namthd_zdf) ** 
    176190   LOGICAL , PUBLIC ::   ln_cndi_U64      !: thermal conductivity: Untersteiner (1964) 
    177191   LOGICAL , PUBLIC ::   ln_cndi_P07      !: thermal conductivity: Pringle et al (2007) 
    178    REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 
    179192   REAL(wp), PUBLIC ::   rn_cnd_s         !: thermal conductivity of the snow [W/m/K]    
     193   REAL(wp), PUBLIC ::   rn_kappa_i       !: coef. for the extinction of radiation in sea ice, Grenfell et al. (2006) [1/m] 
     194   REAL(wp), PUBLIC ::   rn_kappa_s       !: coef. for the extinction of radiation in snw (nn_qtrice=0) [1/m] 
     195   REAL(wp), PUBLIC ::   rn_kappa_smlt    !: coef. for the extinction of radiation in melt snw (nn_qtrice=1) [1/m] 
     196   REAL(wp), PUBLIC ::   rn_kappa_sdry    !: coef. for the extinction of radiation in dry  snw (nn_qtrice=1) [1/m] 
     197   LOGICAL , PUBLIC ::   ln_zdf_chkcvg    !: check convergence of heat diffusion scheme 
    180198 
    181199   !                                     !!** ice-salinity namelist (namthd_sal) ** 
     
    190208   !                                     !!** ice-ponds namelist (namthd_pnd) 
    191209   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F) 
    192    LOGICAL , PUBLIC ::   ln_pnd_H12       !: Melt ponds scheme from Holland et al 2012 
     210   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010) 
     211   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum ice fraction that contributes to melt ponds 
     212   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum ice fraction that contributes to melt ponds 
    193213   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
    194214   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1) 
    195215   REAL(wp), PUBLIC ::   rn_hpnd          !: prescribed pond depth    (0<rn_hpnd<1) 
     216   LOGICAL,  PUBLIC ::   ln_pnd_lids      !: Allow ponds to have frozen lids 
    196217   LOGICAL , PUBLIC ::   ln_pnd_alb       !: melt ponds affect albedo 
    197218 
     
    218239 
    219240   !                                     !!** define arrays 
    220    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_oce,v_oce !: surface ocean velocity used in ice dynamics 
    221    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ht_i_new    !: ice collection thickness accreted in leads 
    222    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   strength    !: ice strength 
    223    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element 
    224    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   delta_i     !: ice rheology elta factor (Flato & Hibler 95) [s-1] 
    225    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   divu_i      !: Divergence of the velocity field             [s-1] 
    226    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   shear_i     !: Shear of the velocity field                  [s-1] 
    227    ! 
    228    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]      
    229    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qlead       !: heat balance of the lead (or of the open ocean) 
    230    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qsb_ice_bot !: net downward heat flux from the ice to the ocean 
    231    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   fhld        !: heat flux from the lead used for bottom melting 
    232  
    233    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw     !: mass flux from snow-ocean mass exchange             [kg.m-2.s-1] 
    234    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sni !: mass flux from snow ice growth component of wfx_snw [kg.m-2.s-1] 
    235    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sum !: mass flux from surface melt component of wfx_snw    [kg.m-2.s-1] 
    236    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_pnd     !: mass flux from melt pond-ocean mass exchange        [kg.m-2.s-1] 
    237    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr     !: mass flux from snow precipitation on ice            [kg.m-2.s-1] 
    238    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sub     !: mass flux from sublimation of snow/ice              [kg.m-2.s-1] 
    239    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_sub !: mass flux from snow sublimation                     [kg.m-2.s-1] 
    240    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice_sub !: mass flux from ice sublimation                      [kg.m-2.s-1] 
    241  
    242    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_snw_dyn !: mass flux from dynamical component of wfx_snw       [kg.m-2.s-1] 
    243  
    244    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_ice     !: mass flux from ice-ocean mass exchange                   [kg.m-2.s-1] 
    245    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sni     !: mass flux from snow ice growth component of wfx_ice      [kg.m-2.s-1] 
    246    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_opw     !: mass flux from lateral ice growth component of wfx_ice   [kg.m-2.s-1] 
    247    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bog     !: mass flux from bottom ice growth component of wfx_ice    [kg.m-2.s-1] 
    248    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_dyn     !: mass flux from dynamical ice growth component of wfx_ice [kg.m-2.s-1] 
    249    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_bom     !: mass flux from bottom melt component of wfx_ice          [kg.m-2.s-1] 
    250    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum     !: mass flux from surface melt component of wfx_ice         [kg.m-2.s-1] 
    251    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_lam     !: mass flux from lateral melt component of wfx_ice         [kg.m-2.s-1] 
    252    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res     !: mass flux from residual component of wfx_ice             [kg.m-2.s-1] 
    253    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_err_sub !: mass flux error after sublimation                        [kg.m-2.s-1] 
    254  
    255    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice bottom growth                   [pss.kg.m-2.s-1 => g.m-2.s-1] 
    256    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bom     !: salt flux due to ice bottom melt                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
    257    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_lam     !: salt flux due to ice lateral melt                    [pss.kg.m-2.s-1 => g.m-2.s-1] 
    258    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sum     !: salt flux due to ice surface melt                    [pss.kg.m-2.s-1 => g.m-2.s-1] 
    259    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sni     !: salt flux due to snow-ice growth                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
    260    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_opw     !: salt flux due to growth in open water                [pss.kg.m-2.s-1 => g.m-2.s-1] 
    261    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bri     !: salt flux due to brine rejection                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
    262    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_dyn     !: salt flux due to porous ridged ice formation         [pss.kg.m-2.s-1 => g.m-2.s-1] 
    263    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: salt flux due to correction on ice thick. (residual) [pss.kg.m-2.s-1 => g.m-2.s-1] 
    264    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_sub     !: salt flux due to ice sublimation                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
    265  
    266    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth           [W.m-2] 
    267    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt             [W.m-2] 
    268    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt            [W.m-2] 
    269    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation    [W.m-2] 
    270    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice      [W.m-2] 
    271    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt                             [W.m-2] 
    272    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_dif !: heat flux remaining due to change in non-solar flux [W.m-2] 
    273    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping => must be 0   [W.m-2] 
    274    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_atm_oi   !: heat flux at the interface atm-[oce+ice]            [W.m-2] 
    275    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   qt_oce_ai   !: heat flux at the interface oce-[atm+ice]            [W.m-2] 
     241   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   u_oce,v_oce     !: surface ocean velocity used in ice dynamics 
     242   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   ht_i_new        !: ice collection thickness accreted in leads 
     243   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   strength        !: ice strength 
     244   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   stress1_i, stress2_i, stress12_i   !: 1st, 2nd & diagonal stress tensor element 
     245   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   delta_i         !: ice rheology elta factor (Flato & Hibler 95) [s-1] 
     246   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   divu_i          !: Divergence of the velocity field             [s-1] 
     247   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   shear_i         !: Shear of the velocity field                  [s-1] 
     248   ! 
     249   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   t_bo            !: Sea-Ice bottom temperature [Kelvin]      
     250   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qlead           !: heat balance of the lead (or of the open ocean) 
     251   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qsb_ice_bot     !: net downward heat flux from the ice to the ocean 
     252   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   fhld            !: heat flux from the lead used for bottom melting 
     253 
     254   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_snw         !: mass flux from snow-ocean mass exchange             [kg.m-2.s-1] 
     255   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_snw_sni     !: mass flux from snow ice growth component of wfx_snw [kg.m-2.s-1] 
     256   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_snw_sum     !: mass flux from surface melt component of wfx_snw    [kg.m-2.s-1] 
     257   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_pnd         !: mass flux from melt pond-ocean mass exchange        [kg.m-2.s-1] 
     258   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_spr         !: mass flux from snow precipitation on ice            [kg.m-2.s-1] 
     259   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_sub         !: mass flux from sublimation of snow/ice              [kg.m-2.s-1] 
     260   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_snw_sub     !: mass flux from snow sublimation                     [kg.m-2.s-1] 
     261   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_ice_sub     !: mass flux from ice sublimation                      [kg.m-2.s-1] 
     262 
     263   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_snw_dyn     !: mass flux from dynamical component of wfx_snw       [kg.m-2.s-1] 
     264 
     265   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_ice         !: mass flux from ice-ocean mass exchange                   [kg.m-2.s-1] 
     266   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_sni         !: mass flux from snow ice growth component of wfx_ice      [kg.m-2.s-1] 
     267   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_opw         !: mass flux from lateral ice growth component of wfx_ice   [kg.m-2.s-1] 
     268   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_bog         !: mass flux from bottom ice growth component of wfx_ice    [kg.m-2.s-1] 
     269   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_dyn         !: mass flux from dynamical ice growth component of wfx_ice [kg.m-2.s-1] 
     270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_bom         !: mass flux from bottom melt component of wfx_ice          [kg.m-2.s-1] 
     271   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_sum         !: mass flux from surface melt component of wfx_ice         [kg.m-2.s-1] 
     272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_lam         !: mass flux from lateral melt component of wfx_ice         [kg.m-2.s-1] 
     273   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_res         !: mass flux from residual component of wfx_ice             [kg.m-2.s-1] 
     274   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   wfx_err_sub     !: mass flux error after sublimation                        [kg.m-2.s-1] 
     275 
     276   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_bog         !: salt flux due to ice bottom growth                   [pss.kg.m-2.s-1 => g.m-2.s-1] 
     277   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_bom         !: salt flux due to ice bottom melt                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     278   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_lam         !: salt flux due to ice lateral melt                    [pss.kg.m-2.s-1 => g.m-2.s-1] 
     279   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_sum         !: salt flux due to ice surface melt                    [pss.kg.m-2.s-1 => g.m-2.s-1] 
     280   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_sni         !: salt flux due to snow-ice growth                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     281   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_opw         !: salt flux due to growth in open water                [pss.kg.m-2.s-1 => g.m-2.s-1] 
     282   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_bri         !: salt flux due to brine rejection                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     283   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_dyn         !: salt flux due to porous ridged ice formation         [pss.kg.m-2.s-1 => g.m-2.s-1] 
     284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_res         !: salt flux due to correction on ice thick. (residual) [pss.kg.m-2.s-1 => g.m-2.s-1] 
     285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   sfx_sub         !: salt flux due to ice sublimation                     [pss.kg.m-2.s-1 => g.m-2.s-1] 
     286 
     287   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_bog         !: total heat flux causing bottom ice growth           [W.m-2] 
     288   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_bom         !: total heat flux causing bottom ice melt             [W.m-2] 
     289   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_sum         !: total heat flux causing surface ice melt            [W.m-2] 
     290   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_opw         !: total heat flux causing open water ice formation    [W.m-2] 
     291   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_dif         !: total heat flux causing Temp change in the ice      [W.m-2] 
     292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_snw         !: heat flux for snow melt                             [W.m-2] 
     293   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_err_dif     !: heat flux remaining due to change in non-solar flux [W.m-2] 
     294   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qt_atm_oi       !: heat flux at the interface atm-[oce+ice]            [W.m-2] 
     295   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   qt_oce_ai       !: heat flux at the interface oce-[atm+ice]            [W.m-2] 
    276296    
    277297   ! heat flux associated with ice-atmosphere mass exchange 
    278    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation            [W.m-2] 
    279    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  [W.m-2] 
     298   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_sub         !: heat flux for sublimation            [W.m-2] 
     299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_spr         !: heat flux of the snow precipitation  [W.m-2] 
    280300 
    281301   ! heat flux associated with ice-ocean mass exchange 
    282    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (icethd_dh) [W.m-2] 
    283    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from ridging                      [W.m-2] 
    284    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: heat flux due to correction on ice thick. (residual)  [W.m-2] 
    285  
    286    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d     !: maximum ice concentration 2d array 
    287    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_ice_bot    !: transmitted solar radiation under ice 
    288    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice         !: temperature of the first layer                (ln_cndflx=T) [K] 
    289    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice        !: effective conductivity at the top of ice/snow (ln_cndflx=T) [W.m-2.K-1] 
     302   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_thd         !: ice-ocean heat flux from thermo processes (icethd_dh) [W.m-2] 
     303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_dyn         !: ice-ocean heat flux from ridging                      [W.m-2] 
     304   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hfx_res         !: heat flux due to correction on ice thick. (residual)  [W.m-2] 
     305 
     306   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   rn_amax_2d      !: maximum ice concentration 2d array 
     307   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qtr_ice_bot     !: transmitted solar radiation under ice 
     308   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice          !: temperature of the first layer          (ln_cndflx=T) [K] 
     309   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 
    290310 
    291311   !!---------------------------------------------------------------------- 
     
    293313   !!---------------------------------------------------------------------- 
    294314   !! Variables defined for each ice category 
    295    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_i       !: Ice thickness                           (m) 
    296    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_i       !: Ice fractional areas (concentration) 
    297    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_i       !: Ice volume per unit area                (m) 
    298    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_s       !: Snow volume per unit area               (m) 
    299    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_s       !: Snow thickness                          (m) 
    300    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_su      !: Sea-Ice Surface Temperature             (K) 
    301    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   s_i       !: Sea-Ice Bulk salinity                   (pss) 
    302    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   sv_i      !: Sea-Ice Bulk salinity * volume per area (pss.m) 
    303    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i       !: Sea-Ice Age                             (s) 
    304    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i      !: Sea-Ice Age times ice area              (s) 
    305    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i      !: brine volume 
     315   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   h_i           !: Ice thickness                           (m) 
     316   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i           !: Ice fractional areas (concentration) 
     317   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_i           !: Ice volume per unit area                (m) 
     318   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s           !: Snow volume per unit area               (m) 
     319   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   h_s           !: Snow thickness                          (m) 
     320   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   t_su          !: Sea-Ice Surface Temperature             (K) 
     321   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   s_i           !: Sea-Ice Bulk salinity                   (pss) 
     322   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sv_i          !: Sea-Ice Bulk salinity * volume per area (pss.m) 
     323   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   o_i           !: Sea-Ice Age                             (s) 
     324   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   oa_i          !: Sea-Ice Age times ice area              (s) 
     325   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   bv_i          !: brine volume 
    306326 
    307327   !! Variables summed over all categories, or associated to all the ice in a single grid cell 
    308    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice !: components of the ice velocity                          (m/s) 
    309    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s  !: ice and snow total volume per unit area                 (m) 
    310    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   st_i         !: Total ice salinity content                              (pss.m) 
    311    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i         !: ice total fractional area (ice concentration) 
    312    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i        !: =1-at_i ; total open water fractional area 
    313    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i , et_s  !: ice and snow total heat content                         (J/m2) 
    314    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i         !: mean ice temperature over all categories                (K) 
    315    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_s         !: mean snw temperature over all categories                (K) 
    316    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bvm_i        !: brine volume averaged over all categories 
    317    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sm_i         !: mean sea ice salinity averaged over all categories      (pss) 
    318    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_su        !: mean surface temperature over all categories            (K) 
    319    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hm_i         !: mean ice  thickness over all categories                 (m) 
    320    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hm_s         !: mean snow thickness over all categories                 (m) 
    321    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories                        (s) 
    322    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tau_icebfr   !: ice friction on ocean bottom (landfast param activated) 
    323  
    324    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s      !: Snow temperatures     [K] 
    325    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s      !: Snow enthalpy         [J/m2] 
    326    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i      !: ice temperatures      [K] 
    327    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i      !: ice enthalpy          [J/m2] 
    328    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sz_i     !: ice salinity          [PSS] 
    329  
    330    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip       !: melt pond concentration 
    331    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   v_ip       !: melt pond volume per grid cell area      [m] 
    332    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   a_ip_frac  !: melt pond fraction (a_ip/a_i) 
    333    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   h_ip       !: melt pond depth                          [m] 
    334  
    335    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   at_ip      !: total melt pond concentration 
    336    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   hm_ip      !: mean melt pond depth                     [m] 
    337    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   vt_ip      !: total melt pond volume per gridcell area [m] 
    338  
    339    !!---------------------------------------------------------------------- 
    340    !! * Old values of global variables 
    341    !!---------------------------------------------------------------------- 
    342    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b, h_ip_b    !: snow and ice volumes/thickness 
    343    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b, oa_i_b                 !: 
    344    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                                 !: snow heat content 
    345    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                                 !: ice temperatures 
    346    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b                      !: ice velocity 
    347    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                                !: ice concentration (total) 
     328   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice, v_ice  !: components of the ice velocity                          (m/s) 
     329   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   vt_i , vt_s   !: ice and snow total volume per unit area                 (m) 
     330   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   st_i          !: Total ice salinity content                              (pss.m) 
     331   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i          !: ice total fractional area (ice concentration) 
     332   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   ato_i         !: =1-at_i ; total open water fractional area 
     333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   et_i , et_s   !: ice and snow total heat content                         (J/m2) 
     334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_i          !: mean ice temperature over all categories                (K) 
     335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_s          !: mean snw temperature over all categories                (K) 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   bvm_i         !: brine volume averaged over all categories 
     337   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   sm_i          !: mean sea ice salinity averaged over all categories      (pss) 
     338   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tm_su         !: mean surface temperature over all categories            (K) 
     339   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hm_i          !: mean ice  thickness over all categories                 (m) 
     340   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hm_s          !: mean snow thickness over all categories                 (m) 
     341   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   om_i          !: mean ice age over all categories                        (s) 
     342   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   tau_icebfr    !: ice friction on ocean bottom (landfast param activated) 
     343 
     344   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s           !: Snow temperatures     [K] 
     345   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s           !: Snow enthalpy         [J/m2] 
     346   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_i           !: ice temperatures      [K] 
     347   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i           !: ice enthalpy          [J/m2] 
     348   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   sz_i          !: ice salinity          [PSS] 
     349 
     350   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_ip          !: melt pond concentration 
     351   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_ip          !: melt pond volume per grid cell area      [m] 
     352   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_ip_frac     !: melt pond fraction (a_ip/a_i) 
     353   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_ip_eff      !: melt pond effective fraction (not covered up by lid) (a_ip/a_i) 
     354   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   h_ip          !: melt pond depth                          [m] 
     355   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_il          !: melt pond lid volume                     [m] 
     356   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   h_il          !: melt pond lid thickness                  [m] 
     357 
     358   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_ip         !: total melt pond concentration 
     359   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hm_ip         !: mean melt pond depth                     [m] 
     360   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   vt_ip         !: total melt pond volume per gridcell area [m] 
     361   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   hm_il         !: mean melt pond lid depth                     [m] 
     362   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   vt_il         !: total melt pond lid volume per gridcell area [m] 
     363 
     364   !!---------------------------------------------------------------------- 
     365   !! * Global variables at before time step 
     366   !!---------------------------------------------------------------------- 
     367   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 
     368   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   a_i_b, sv_i_b              !: 
     369   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_s_b                      !: snow heat content 
     370   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   e_i_b                      !: ice temperatures 
     371   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   u_ice_b, v_ice_b           !: ice velocity 
     372   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)     ::   at_i_b                     !: ice concentration (total) 
    348373             
    349374   !!---------------------------------------------------------------------- 
    350375   !! * Ice thickness distribution variables 
    351376   !!---------------------------------------------------------------------- 
    352    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hi_max         !: Boundary of ice thickness categories in thickness space 
    353    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hi_mean        !: Mean ice thickness in catgories  
     377   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_max            !: Boundary of ice thickness categories in thickness space 
     378   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:)   ::   hi_mean           !: Mean ice thickness in catgories  
    354379   ! 
    355380   !!---------------------------------------------------------------------- 
    356381   !! * Ice diagnostics 
    357382   !!---------------------------------------------------------------------- 
    358    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi   !: transport of ice volume 
    359    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs   !: transport of snw volume 
    360    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_ei   !: transport of ice enthalpy [W/m2] 
    361    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_es   !: transport of snw enthalpy [W/m2] 
    362    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_sv   !: transport of salt content 
    363    ! 
    364    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_heat     !: snw/ice heat content variation   [W/m2]  
    365    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_sice     !: ice salt content variation   []  
    366    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice     !: ice volume variation   [m/s]  
    367    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw     !: snw volume variation   [m/s]  
    368  
     383   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vi       !: transport of ice volume 
     384   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_vs       !: transport of snw volume 
     385   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_ei       !: transport of ice enthalpy [W/m2] 
     386   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_es       !: transport of snw enthalpy [W/m2] 
     387   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_trp_sv       !: transport of salt content 
     388   ! 
     389   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_heat         !: snw/ice heat content variation   [W/m2]  
     390   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_sice         !: ice salt content variation   []  
     391   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vice         !: ice volume variation   [m/s]  
     392   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_vsnw         !: snw volume variation   [m/s]  
     393   ! 
    369394   !!---------------------------------------------------------------------- 
    370395   !! * Ice conservation 
    371396   !!---------------------------------------------------------------------- 
    372    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_v        !: conservation of ice volume 
    373    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_s        !: conservation of ice salt 
    374    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_t        !: conservation of ice heat 
    375    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fv       !: conservation of ice volume 
    376    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fs       !: conservation of ice salt 
    377    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_ft       !: conservation of ice heat 
     397   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_v            !: conservation of ice volume 
     398   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_s            !: conservation of ice salt 
     399   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_t            !: conservation of ice heat 
     400   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fv           !: conservation of ice volume 
     401   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_fs           !: conservation of ice salt 
     402   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   diag_ft           !: conservation of ice heat 
    378403   ! 
    379404   !!---------------------------------------------------------------------- 
     
    381406   !!---------------------------------------------------------------------- 
    382407   ! Extra sea ice diagnostics to address the data request 
    383    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_si          !: Temperature at Snow-ice interface (K)  
    384    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tm_si         !: mean temperature at the snow-ice interface (K)  
    385    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_bot   !: Bottom  conduction flux (W/m2) 
    386    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_top   !: Surface conduction flux (W/m2) 
    387  
     408   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t_si            !: Temperature at Snow-ice interface (K)  
     409   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   tm_si           !: mean temperature at the snow-ice interface (K)  
     410   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_bot     !: Bottom  conduction flux (W/m2) 
     411   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   qcn_ice_top     !: Surface conduction flux (W/m2) 
    388412   ! 
    389413   !!---------------------------------------------------------------------- 
     
    424448         &      hfx_sum    (jpi,jpj) , hfx_bom   (jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) ,     & 
    425449         &      hfx_opw    (jpi,jpj) , hfx_thd   (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,     & 
    426          &      hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj)             , STAT=ierr(ii) ) 
     450         &      hfx_err_dif(jpi,jpj) , wfx_err_sub(jpi,jpj)                   , STAT=ierr(ii) ) 
    427451 
    428452      ! * Ice global state variables 
     
    448472 
    449473      ii = ii + 1 
    450       ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl) , STAT = ierr(ii) ) 
    451  
    452       ii = ii + 1 
    453       ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , STAT = ierr(ii) ) 
     474      ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl),  & 
     475         &      v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) ) 
     476 
     477      ii = ii + 1 
     478      ALLOCATE( at_ip(jpi,jpj) , hm_ip(jpi,jpj) , vt_ip(jpi,jpj) , hm_il(jpi,jpj) , vt_il(jpi,jpj) , STAT = ierr(ii) ) 
    454479 
    455480      ! * Old values of global variables 
    456481      ii = ii + 1 
    457       ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl), h_ip_b(jpi,jpj,jpl),  & 
    458          &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) ,               & 
    459          &      oa_i_b(jpi,jpj,jpl)                                                   , STAT=ierr(ii) ) 
     482      ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl)        , h_i_b(jpi,jpj,jpl),         & 
     483         &      a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 
     484         &      STAT=ierr(ii) ) 
    460485 
    461486      ii = ii + 1 
     
    484509      IF( ice_alloc /= 0 )   CALL ctl_stop( 'STOP', 'ice_alloc: failed to allocate arrays.' ) 
    485510      ! 
     511 
    486512   END FUNCTION ice_alloc 
    487513 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/ice1d.F90

    r10786 r13553  
    5151   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d 
    5252   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dyn_1d 
    53    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d 
    5453   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_dif_1d 
    5554   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   qt_oce_ai_1d 
     
    124123   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   oa_i_1d       !: 
    125124   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   o_i_1d        !: 
    126    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !: 
     125   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_1d       !: ice ponds 
    127126   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_ip_1d       !: 
    128127   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_ip_1d       !: 
    129    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   a_ip_frac_1d  !: 
     128   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   v_il_1d       !: Ice pond lid 
     129   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   h_il_1d       !: 
    130130 
    131131   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_s_1d      !: corresponding to the 2D var  t_s 
     
    146146   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sss_1d 
    147147 
     148   ! convergence check 
     149   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tice_cvgerr_1d   !: convergence of ice/snow temp (dT)          [K] 
     150   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   tice_cvgstp_1d   !: convergence of ice/snow temp (subtimestep) [-] 
    148151   !  
    149152   !!---------------------- 
     
    157160   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   a_ip_2d 
    158161   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_ip_2d  
     162   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   v_il_2d  
    159163   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_su_2d  
    160164   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   h_i_2d 
     
    175179      !!---------------------------------------------------------------------! 
    176180      INTEGER ::   ice1D_alloc   ! return value 
    177       INTEGER ::   ierr(7), ii 
     181      INTEGER ::   ierr(8), ii 
    178182      !!---------------------------------------------------------------------! 
    179183      ierr(:) = 0 
     
    189193         &      hfx_thd_1d(jpij) , hfx_spr_1d    (jpij) ,                      & 
    190194         &      hfx_snw_1d(jpij) , hfx_sub_1d    (jpij) ,                      & 
    191          &      hfx_res_1d(jpij) , hfx_err_rem_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) ) 
     195         &      hfx_res_1d(jpij) , hfx_err_dif_1d(jpij) , qt_oce_ai_1d(jpij), STAT=ierr(ii) ) 
    192196      ! 
    193197      ii = ii + 1 
     
    208212         &      dh_s_tot(jpij) , dh_i_sum(jpij) , dh_i_itm  (jpij) , dh_i_bom(jpij) , dh_i_bog(jpij) ,  &     
    209213         &      dh_i_sub(jpij) , dh_s_mlt(jpij) , dh_snowice(jpij) , s_i_1d  (jpij) , s_i_new (jpij) ,  & 
    210          &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) ,                   & 
    211          &      h_ip_1d (jpij) , a_ip_frac_1d(jpij) ,                                                   & 
     214         &      a_ip_1d (jpij) , v_ip_1d (jpij) , v_i_1d    (jpij) , v_s_1d  (jpij) , v_il_1d (jpij) ,  & 
     215         &      h_il_1d (jpij) , h_ip_1d (jpij) ,                                                       & 
    212216         &      sv_i_1d (jpij) , oa_i_1d (jpij) , o_i_1d    (jpij) , STAT=ierr(ii) ) 
    213217      ! 
     
    224228      ! 
    225229      ii = ii + 1 
     230      ALLOCATE( tice_cvgerr_1d(jpij) , tice_cvgstp_1d(jpij) , STAT=ierr(ii) ) 
     231      ! 
     232      ii = ii + 1 
    226233      ALLOCATE( a_i_2d (jpij,jpl) , a_ib_2d(jpij,jpl) , h_i_2d (jpij,jpl) , h_ib_2d(jpij,jpl) ,  & 
    227234         &      v_i_2d (jpij,jpl) , v_s_2d (jpij,jpl) , oa_i_2d(jpij,jpl) , sv_i_2d(jpij,jpl) ,  & 
    228          &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) ,                      & 
     235         &      a_ip_2d(jpij,jpl) , v_ip_2d(jpij,jpl) , t_su_2d(jpij,jpl) , v_il_2d(jpij,jpl) ,  & 
    229236         &      STAT=ierr(ii) ) 
    230237 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icealb.F90

    r13295 r13553  
    1414   !!   ice_alb_init   : initialisation of albedo computation 
    1515   !!---------------------------------------------------------------------- 
    16    USE ice, ONLY: jpl ! sea-ice: number of categories 
    1716   USE phycst         ! physical constants 
    1817   USE dom_oce        ! domain: ocean 
     18   USE ice, ONLY: jpl ! sea-ice: number of categories 
     19   USE icevar         ! sea-ice: operations 
    1920   ! 
    2021   USE in_out_manager ! I/O manager 
     
    4748CONTAINS 
    4849 
    49    SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, palb_cs, palb_os ) 
     50   SUBROUTINE ice_alb( pt_su, ph_ice, ph_snw, ld_pnd_alb, pafrac_pnd, ph_pnd, pcloud_fra, palb_ice ) 
    5051      !!---------------------------------------------------------------------- 
    5152      !!               ***  ROUTINE ice_alb  *** 
     
    99100      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   pafrac_pnd   !  melt pond relative fraction (per unit ice area) 
    100101      REAL(wp), INTENT(in   ), DIMENSION(:,:,:) ::   ph_pnd       !  melt pond depth 
    101       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_cs      !  albedo of ice under clear    sky 
    102       REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_os      !  albedo of ice under overcast sky 
    103       ! 
     102      REAL(wp), INTENT(in   ), DIMENSION(:,:)   ::   pcloud_fra   !  cloud fraction 
     103      REAL(wp), INTENT(  out), DIMENSION(:,:,:) ::   palb_ice     !  albedo of ice 
     104      ! 
     105      REAL(wp), DIMENSION(jpi,jpj,jpl) :: za_s_fra   ! ice fraction covered by snow 
    104106      INTEGER  ::   ji, jj, jl                ! dummy loop indices 
    105107      REAL(wp) ::   z1_c1, z1_c2,z1_c3, z1_c4 ! local scalar 
     
    108110      REAL(wp) ::   zalb_ice, zafrac_ice      ! bare sea ice albedo & relative ice fraction 
    109111      REAL(wp) ::   zalb_snw, zafrac_snw      ! snow-covered sea ice albedo & relative snow fraction 
     112      REAL(wp) ::   zalb_cs, zalb_os          ! albedo of ice under clear/overcast sky 
    110113      !!--------------------------------------------------------------------- 
    111114      ! 
     
    118121      z1_c4 = 1. / 0.03 
    119122      ! 
     123      CALL ice_var_snwfra( ph_snw, za_s_fra )   ! calculate ice fraction covered by snow 
     124      ! 
    120125      DO jl = 1, jpl 
    121126         DO_2D( 1, 1, 1, 1 ) 
    122             !                       !--- Specific snow, ice and pond fractions (for now, we prevent melt ponds and snow at the same time) 
    123             IF( ph_snw(ji,jj,jl) == 0._wp ) THEN 
    124                zafrac_snw = 0._wp 
    125                IF( ld_pnd_alb ) THEN 
    126                   zafrac_pnd = pafrac_pnd(ji,jj,jl) 
    127                ELSE 
    128                   zafrac_pnd = 0._wp 
    129                ENDIF 
    130                zafrac_ice = 1._wp - zafrac_pnd 
     127            ! 
     128            !---------------------------------------------! 
     129            !--- Specific snow, ice and pond fractions ---! 
     130            !---------------------------------------------!                
     131            zafrac_snw = za_s_fra(ji,jj,jl) 
     132            IF( ld_pnd_alb ) THEN 
     133               zafrac_pnd = MIN( pafrac_pnd(ji,jj,jl), 1._wp - zafrac_snw ) ! make sure (a_ip_eff + a_s_fra) <= 1 
    131134            ELSE 
    132                zafrac_snw = 1._wp      ! Snow fully "shades" melt ponds and ice 
    133135               zafrac_pnd = 0._wp 
    134                zafrac_ice = 0._wp 
    135             ENDIF 
    136             ! 
     136            ENDIF 
     137            zafrac_ice = MAX( 0._wp, 1._wp - zafrac_pnd - zafrac_snw ) ! max for roundoff errors 
     138            ! 
     139            !---------------! 
     140            !--- Albedos ---! 
     141            !---------------!                
    137142            !                       !--- Bare ice albedo (for hi > 150cm) 
    138143            IF( ld_pnd_alb ) THEN 
    139144               zalb_ice = rn_alb_idry 
    140145            ELSE 
    141                IF( ph_snw(ji,jj,jl) == 0._wp .AND. pt_su(ji,jj,jl) >= rt0 ) THEN  ;   zalb_ice = rn_alb_imlt 
    142                ELSE                                                               ;   zalb_ice = rn_alb_idry   ;   ENDIF 
     146               IF( ph_snw(ji,jj,jl) == 0._wp .AND. pt_su(ji,jj,jl) >= rt0 ) THEN   ;   zalb_ice = rn_alb_imlt 
     147               ELSE                                                                ;   zalb_ice = rn_alb_idry   ;   ENDIF 
    143148            ENDIF 
    144149            !                       !--- Bare ice albedo (for hi < 150cm) 
     
    156161            ENDIF 
    157162            !                       !--- Ponded ice albedo 
    158             IF( ld_pnd_alb ) THEN 
    159                zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd )  
    160             ELSE 
    161                zalb_pnd = rn_alb_dpnd 
    162             ENDIF 
     163            zalb_pnd = rn_alb_dpnd - ( rn_alb_dpnd - zalb_ice ) * EXP( - ph_pnd(ji,jj,jl) * z1_href_pnd )  
     164            ! 
    163165            !                       !--- Surface albedo is weighted mean of snow, ponds and bare ice contributions 
    164             palb_os(ji,jj,jl) = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1) 
    165             ! 
    166             palb_cs(ji,jj,jl) = palb_os(ji,jj,jl)  & 
    167                &                - ( - 0.1010 * palb_os(ji,jj,jl) * palb_os(ji,jj,jl)  & 
    168                &                    + 0.1933 * palb_os(ji,jj,jl) - 0.0148 ) * tmask(ji,jj,1) 
    169             ! 
     166            zalb_os = ( zafrac_snw * zalb_snw + zafrac_pnd * zalb_pnd + zafrac_ice * zalb_ice ) * tmask(ji,jj,1) 
     167            ! 
     168            zalb_cs = zalb_os - ( - 0.1010 * zalb_os * zalb_os  & 
     169               &                  + 0.1933 * zalb_os - 0.0148 ) * tmask(ji,jj,1) 
     170            ! 
     171            ! albedo depends on cloud fraction because of non-linear spectral effects 
     172            palb_ice(ji,jj,jl) = ( 1._wp - pcloud_fra(ji,jj) ) * zalb_cs + pcloud_fra(ji,jj) * zalb_os 
     173 
    170174         END_2D 
    171175      END DO 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icecor.F90

    r13295 r13553  
    8181      DO jl = 1, jpl 
    8282         WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 
    83       END DO 
    84      
     83      END DO     
     84      !                             !----------------------------------------------------- 
     85      !                             !  Rebin categories with thickness out of bounds     ! 
     86      !                             !----------------------------------------------------- 
     87      IF ( jpl > 1 )   CALL ice_itd_reb( kt ) 
     88      ! 
    8589      !                             !----------------------------------------------------- 
    8690      IF ( nn_icesal == 2 ) THEN    !  salinity must stay in bounds [Simin,Simax]        ! 
     
    96100      ENDIF 
    97101      !                             !----------------------------------------------------- 
    98       !                             !  Rebin categories with thickness out of bounds     ! 
    99       !                             !----------------------------------------------------- 
    100       IF ( jpl > 1 )   CALL ice_itd_reb( kt ) 
    101  
    102       !                             !----------------------------------------------------- 
    103102      CALL ice_var_zapsmall         !  Zap small values                                  ! 
    104103      !                             !----------------------------------------------------- 
     
    106105      !                             !----------------------------------------------------- 
    107106      IF( kn == 2 ) THEN            !  Ice drift case: Corrections to avoid wrong values ! 
    108          DO_2D( 0, 0, 0, 0 ) 
     107         DO_2D( 0, 0, 0, 0 )        !----------------------------------------------------- 
    109108            IF ( at_i(ji,jj) == 0._wp ) THEN    ! what to do if there is no ice 
    110109               IF ( at_i(ji+1,jj) == 0._wp )   u_ice(ji  ,jj) = 0._wp   ! right side 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icectl.F90

    r13295 r13553  
    350350      !!                   ***  ROUTINE ice_ctl ***  
    351351      !!                  
    352       !! ** Purpose :   Alerts in case of model crash 
     352      !! ** Purpose :   control checks 
    353353      !!------------------------------------------------------------------- 
    354354      INTEGER, INTENT(in) ::   kt      ! ocean time step 
    355       INTEGER  ::   ji, jj, jk,  jl   ! dummy loop indices 
    356       INTEGER  ::   inb_altests       ! number of alert tests (max 20) 
    357       INTEGER  ::   ialert_id         ! number of the current alert 
    358       REAL(wp) ::   ztmelts           ! ice layer melting point 
     355      INTEGER  ::   ja, ji, jj, jk, jl ! dummy loop indices 
     356      INTEGER  ::   ialert_id          ! number of the current alert 
     357      REAL(wp) ::   ztmelts            ! ice layer melting point 
    359358      CHARACTER (len=30), DIMENSION(20) ::   cl_alname   ! name of alert 
    360359      INTEGER           , DIMENSION(20) ::   inb_alp     ! number of alerts positive 
    361360      !!------------------------------------------------------------------- 
    362  
    363       inb_altests = 10 
    364       inb_alp(:)  =  0 
    365  
    366       ! Alert if incompatible volume and concentration 
    367       ialert_id = 2 ! reference number of this alert 
    368       cl_alname(ialert_id) = ' Incompat vol and con         '    ! name of the alert 
     361      inb_alp(:) = 0 
     362      ialert_id = 0 
     363       
     364      ! Alert if very high salinity 
     365      ialert_id = ialert_id + 1 ! reference number of this alert 
     366      cl_alname(ialert_id) = ' Very high salinity ' ! name of the alert 
    369367      DO jl = 1, jpl 
    370368         DO_2D( 1, 1, 1, 1 ) 
    371             IF(  v_i(ji,jj,jl) /= 0._wp   .AND.   a_i(ji,jj,jl) == 0._wp   ) THEN 
    372                WRITE(numout,*) ' ALERTE 2 :   Incompatible volume and concentration ' 
    373                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     369            IF( v_i(ji,jj,jl) > epsi10  ) THEN 
     370               IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) > rn_simax ) THEN 
     371                  WRITE(numout,*) ' ALERTE :   Very high salinity ',sv_i(ji,jj,jl)/v_i(ji,jj,jl) 
     372                  WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 
     373                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     374               ENDIF 
    374375            ENDIF 
    375376         END_2D 
    376377      END DO 
    377378 
    378       ! Alerte if very thick ice 
    379       ialert_id = 3 ! reference number of this alert 
    380       cl_alname(ialert_id) = ' Very thick ice               ' ! name of the alert 
    381       jl = jpl  
    382       DO_2D( 1, 1, 1, 1 ) 
    383          IF(   h_i(ji,jj,jl)  >  50._wp   ) THEN 
    384             WRITE(numout,*) ' ALERTE 3 :   Very thick ice' 
    385             !CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 :   Very thick ice ' ) 
    386             inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    387          ENDIF 
    388       END_2D 
    389  
    390       ! Alert if very fast ice 
    391       ialert_id = 4 ! reference number of this alert 
    392       cl_alname(ialert_id) = ' Very fast ice               ' ! name of the alert 
    393       DO_2D( 1, 1, 1, 1 ) 
    394          IF(   MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2.  .AND.  & 
    395             &  at_i(ji,jj) > 0._wp   ) THEN 
    396             WRITE(numout,*) ' ALERTE 4 :   Very fast ice' 
    397             !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 :   Very fast ice ' ) 
    398             inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    399          ENDIF 
    400       END_2D 
    401  
    402       ! Alert on salt flux 
    403       ialert_id = 5 ! reference number of this alert 
    404       cl_alname(ialert_id) = ' High salt flux               ' ! name of the alert 
    405       DO_2D( 1, 1, 1, 1 ) 
    406          IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN  ! = 1 psu/day for 1m ocean depth 
    407             WRITE(numout,*) ' ALERTE 5 :   High salt flux' 
    408             !CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 :   High salt flux ' ) 
    409             inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    410          ENDIF 
    411       END_2D 
    412  
    413       ! Alert if there is ice on continents 
    414       ialert_id = 6 ! reference number of this alert 
    415       cl_alname(ialert_id) = ' Ice on continents           ' ! name of the alert 
    416       DO_2D( 1, 1, 1, 1 ) 
    417          IF(   tmask(ji,jj,1) <= 0._wp   .AND.   at_i(ji,jj) > 0._wp   ) THEN  
    418             WRITE(numout,*) ' ALERTE 6 :   Ice on continents' 
    419             !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 :   Ice on continents ' ) 
    420             inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    421          ENDIF 
    422       END_2D 
    423  
    424 ! 
    425 !     ! Alert if very fresh ice 
    426       ialert_id = 7 ! reference number of this alert 
    427       cl_alname(ialert_id) = ' Very fresh ice               ' ! name of the alert 
     379      ! Alert if very low salinity 
     380      ialert_id = ialert_id + 1 ! reference number of this alert 
     381      cl_alname(ialert_id) = ' Very low salinity ' ! name of the alert 
    428382      DO jl = 1, jpl 
    429383         DO_2D( 1, 1, 1, 1 ) 
    430             IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN 
    431                WRITE(numout,*) ' ALERTE 7 :   Very fresh ice' 
    432 !                 CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 :   Very fresh ice ' ) 
    433                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     384            IF( v_i(ji,jj,jl) > epsi10  ) THEN 
     385               IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) < rn_simin ) THEN 
     386                  WRITE(numout,*) ' ALERTE :   Very low salinity ',sv_i(ji,jj,jl),v_i(ji,jj,jl) 
     387                  WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 
     388                  inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     389               ENDIF 
    434390            ENDIF 
    435391         END_2D 
    436392      END DO 
    437 ! 
    438       ! Alert if qns very big 
    439       ialert_id = 8 ! reference number of this alert 
    440       cl_alname(ialert_id) = ' fnsolar very big             ' ! name of the alert 
    441       DO_2D( 1, 1, 1, 1 ) 
    442          IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN 
    443             ! 
    444             WRITE(numout,*) ' ALERTE 8 :   Very high non-solar heat flux' 
    445             !CALL ice_prt( kt, ji, jj, 2, '   ') 
    446             inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    447             ! 
    448          ENDIF 
    449       END_2D 
    450       !+++++ 
    451  
    452 !     ! Alert if too old ice 
    453       ialert_id = 9 ! reference number of this alert 
    454       cl_alname(ialert_id) = ' Very old   ice               ' ! name of the alert 
    455       DO jl = 1, jpl 
    456          DO_2D( 1, 1, 1, 1 ) 
    457             IF ( ( ( ABS( o_i(ji,jj,jl) ) > rDt_ice ) .OR. & 
    458                    ( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. & 
    459                           ( a_i(ji,jj,jl) > 0._wp ) ) THEN 
    460                WRITE(numout,*) ' ALERTE 9 :   Wrong ice age' 
    461                !CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 :   Wrong ice age ') 
    462                inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    463             ENDIF 
    464          END_2D 
    465       END DO 
    466    
    467       ! Alert if very warm ice 
    468       ialert_id = 10 ! reference number of this alert 
    469       cl_alname(ialert_id) = ' Very warm ice                ' ! name of the alert 
    470       inb_alp(ialert_id) = 0 
     393 
     394      ! Alert if very cold ice 
     395      ialert_id = ialert_id + 1 ! reference number of this alert 
     396      cl_alname(ialert_id) = ' Very cold ice ' ! name of the alert 
    471397      DO jl = 1, jpl 
    472398         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
    473399            ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0 
    474             IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > 1.e-10   & 
    475                &                            .AND.  a_i(ji,jj,jl) > 0._wp   ) THEN 
    476                WRITE(numout,*) ' ALERTE 10 :   Very warm ice' 
     400            IF( t_i(ji,jj,jk,jl) < -50.+rt0  .AND.  v_i(ji,jj,jl) > epsi10 ) THEN 
     401               WRITE(numout,*) ' ALERTE :   Very cold ice ',(t_i(ji,jj,jk,jl)-rt0) 
     402               WRITE(numout,*) ' at i,j,k,l = ',ji,jj,jk,jl 
    477403              inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
    478404            ENDIF 
    479405         END_3D 
    480406      END DO 
     407   
     408      ! Alert if very warm ice 
     409      ialert_id = ialert_id + 1 ! reference number of this alert 
     410      cl_alname(ialert_id) = ' Very warm ice ' ! name of the alert 
     411      DO jl = 1, jpl 
     412         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     413            ztmelts    =  -rTmlt * sz_i(ji,jj,jk,jl) + rt0 
     414            IF( t_i(ji,jj,jk,jl) > ztmelts  .AND.  v_i(ji,jj,jl) > epsi10 ) THEN 
     415               WRITE(numout,*) ' ALERTE :   Very warm ice',(t_i(ji,jj,jk,jl)-rt0) 
     416               WRITE(numout,*) ' at i,j,k,l = ',ji,jj,jk,jl 
     417              inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     418            ENDIF 
     419         END_3D 
     420      END DO 
     421       
     422      ! Alerte if very thick ice 
     423      ialert_id = ialert_id + 1 ! reference number of this alert 
     424      cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert 
     425      jl = jpl  
     426      DO_2D( 1, 1, 1, 1 ) 
     427         IF( h_i(ji,jj,jl) > 50._wp ) THEN 
     428            WRITE(numout,*) ' ALERTE :   Very thick ice ',h_i(ji,jj,jl) 
     429            WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 
     430            inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     431         ENDIF 
     432      END_2D 
     433 
     434      ! Alerte if very thin ice 
     435      ialert_id = ialert_id + 1 ! reference number of this alert 
     436      cl_alname(ialert_id) = ' Very thin ice ' ! name of the alert 
     437      jl = 1  
     438      DO_2D( 1, 1, 1, 1 ) 
     439         IF( h_i(ji,jj,jl) < rn_himin ) THEN 
     440            WRITE(numout,*) ' ALERTE :   Very thin ice ',h_i(ji,jj,jl) 
     441            WRITE(numout,*) ' at i,j,l = ',ji,jj,jl 
     442            inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     443         ENDIF 
     444      END_2D 
     445 
     446      ! Alert if very fast ice 
     447      ialert_id = ialert_id + 1 ! reference number of this alert 
     448      cl_alname(ialert_id) = ' Very fast ice ' ! name of the alert 
     449      DO_2D( 1, 1, 1, 1 ) 
     450         IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2. ) THEN 
     451            WRITE(numout,*) ' ALERTE :   Very fast ice ',MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) 
     452            WRITE(numout,*) ' at i,j = ',ji,jj 
     453            inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     454         ENDIF 
     455      END_2D 
     456 
     457      ! Alert if there is ice on continents 
     458      ialert_id = ialert_id + 1 ! reference number of this alert 
     459      cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert 
     460      DO_2D( 1, 1, 1, 1 ) 
     461         IF( tmask(ji,jj,1) == 0._wp .AND. ( at_i(ji,jj) > 0._wp .OR. vt_i(ji,jj) > 0._wp ) ) THEN  
     462            WRITE(numout,*) ' ALERTE :   Ice on continents ',at_i(ji,jj),vt_i(ji,jj) 
     463            WRITE(numout,*) ' at i,j = ',ji,jj 
     464            inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     465         ENDIF 
     466      END_2D 
     467 
     468      ! Alert if incompatible ice concentration and volume 
     469      ialert_id = ialert_id + 1 ! reference number of this alert 
     470      cl_alname(ialert_id) = ' Incompatible ice conc and vol ' ! name of the alert 
     471      DO_2D( 1, 1, 1, 1 ) 
     472         IF(  ( vt_i(ji,jj) == 0._wp .AND. at_i(ji,jj) >  0._wp ) .OR. & 
     473            & ( vt_i(ji,jj) >  0._wp .AND. at_i(ji,jj) == 0._wp ) ) THEN  
     474            WRITE(numout,*) ' ALERTE :   Incompatible ice conc and vol ',at_i(ji,jj),vt_i(ji,jj) 
     475            WRITE(numout,*) ' at i,j = ',ji,jj 
     476            inb_alp(ialert_id) = inb_alp(ialert_id) + 1 
     477         ENDIF 
     478      END_2D 
    481479 
    482480      ! sum of the alerts on all processors 
    483481      IF( lk_mpp ) THEN 
    484          DO ialert_id = 1, inb_altests 
    485             CALL mpp_sum('icectl', inb_alp(ialert_id)) 
     482         DO ja = 1, ialert_id 
     483            CALL mpp_sum('icectl', inb_alp(ja)) 
    486484         END DO 
    487485      ENDIF 
     
    489487      ! print alerts 
    490488      IF( lwp ) THEN 
    491          ialert_id = 1                                 ! reference number of this alert 
    492          cl_alname(ialert_id) = ' NO alerte 1      '   ! name of the alert 
    493489         WRITE(numout,*) ' time step ',kt 
    494490         WRITE(numout,*) ' All alerts at the end of ice model ' 
    495          DO ialert_id = 1, inb_altests 
    496             WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! ' 
     491         DO ja = 1, ialert_id 
     492            WRITE(numout,*) ja, cl_alname(ja)//' : ', inb_alp(ja), ' times ! ' 
    497493         END DO 
    498494      ENDIF 
     
    543539               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    544540               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    545                WRITE(numout,*) 
    546541               WRITE(numout,*) ' - Cell values ' 
    547542               WRITE(numout,*) '   ~~~~~~~~~~~ ' 
     
    552547               DO jl = 1, jpl 
    553548                  WRITE(numout,*) ' - Category (', jl,')' 
     549                  WRITE(numout,*) '   ~~~~~~~~~~~ ' 
    554550                  WRITE(numout,*) ' a_i           : ', a_i(ji,jj,jl) 
    555551                  WRITE(numout,*) ' h_i           : ', h_i(ji,jj,jl) 
     
    588584               WRITE(numout,*) ' v_ice(i  ,j)  : ', v_ice(ji,jj) 
    589585               WRITE(numout,*) ' strength      : ', strength(ji,jj) 
    590                WRITE(numout,*) ' u_ice_b       : ', u_ice_b(ji,jj)    , ' v_ice_b       : ', v_ice_b(ji,jj)   
    591586               WRITE(numout,*) 
    592587                
     
    605600                  WRITE(numout,*) ' e_snow     : ', e_s(ji,jj,1,jl)            , ' e_snow_b   : ', e_s_b(ji,jj,1,jl)  
    606601                  WRITE(numout,*) ' sv_i       : ', sv_i(ji,jj,jl)             , ' sv_i_b     : ', sv_i_b(ji,jj,jl)    
    607                   WRITE(numout,*) ' oa_i       : ', oa_i(ji,jj,jl)             , ' oa_i_b     : ', oa_i_b(ji,jj,jl) 
    608602               END DO !jl 
    609603                
     
    713707         CALL prt_ctl(tab2d_1=v_i        (:,:,jl)        , clinfo1= ' v_i         : ') 
    714708         CALL prt_ctl(tab2d_1=v_s        (:,:,jl)        , clinfo1= ' v_s         : ') 
    715          CALL prt_ctl(tab2d_1=e_i        (:,:,1,jl)      , clinfo1= ' e_i1        : ') 
    716709         CALL prt_ctl(tab2d_1=e_s        (:,:,1,jl)      , clinfo1= ' e_snow      : ') 
    717710         CALL prt_ctl(tab2d_1=sv_i       (:,:,jl)        , clinfo1= ' sv_i        : ') 
     
    721714            CALL prt_ctl_info(' - Layer : ', ivar=jk) 
    722715            CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i       : ') 
     716            CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' e_i       : ') 
    723717         END DO 
    724718      END DO 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn.F90

    r13295 r13553  
    100100      WHERE( a_ip(:,:,:) >= epsi20 ) 
    101101         h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
     102         h_il(:,:,:) = v_il(:,:,:) / a_ip(:,:,:) 
    102103      ELSEWHERE 
    103104         h_ip(:,:,:) = 0._wp 
     105         h_il(:,:,:) = 0._wp 
    104106      END WHERE 
    105107      ! 
     
    127129         ! Then for dx = 2m and dt = 1s => rn_uice = u (1/6th) = 1m/s  
    128130         DO_2D( 1, 1, 1, 1 ) 
    129             zcoefu = ( REAL(jpiglo+1)*0.5 - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5 - 1. ) 
    130             zcoefv = ( REAL(jpjglo+1)*0.5 - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5 - 1. ) 
    131             u_ice(ji,jj) = rn_uice * 1.5 * SIGN( 1.0_wp, zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
    132             v_ice(ji,jj) = rn_vice * 1.5 * SIGN( 1.0_wp, zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
     131            zcoefu = ( REAL(jpiglo+1)*0.5_wp - REAL(ji+nimpp-1) ) / ( REAL(jpiglo+1)*0.5_wp - 1._wp ) 
     132            zcoefv = ( REAL(jpjglo+1)*0.5_wp - REAL(jj+njmpp-1) ) / ( REAL(jpjglo+1)*0.5_wp - 1._wp ) 
     133            u_ice(ji,jj) = rn_uice * 1.5_wp * SIGN( 1.0_wp, zcoefu ) * ABS( zcoefu ) * umask(ji,jj,1) 
     134            v_ice(ji,jj) = rn_vice * 1.5_wp * SIGN( 1.0_wp, zcoefv ) * ABS( zcoefv ) * vmask(ji,jj,1) 
    133135         END_2D 
    134136         ! --- 
     
    218220      NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice,  & 
    219221         &             rn_ishlat ,                                                           & 
    220          &             ln_landfast_L16, rn_depfra, rn_icebfr, rn_lfrelax, rn_tensile 
     222         &             ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile 
    221223      !!------------------------------------------------------------------- 
    222224      ! 
     
    239241         WRITE(numout,*) '      lateral boundary condition for sea ice dynamics        rn_ishlat       = ', rn_ishlat 
    240242         WRITE(numout,*) '      Landfast: param from Lemieux 2016                      ln_landfast_L16 = ', ln_landfast_L16 
    241          WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_depfra       = ', rn_depfra 
    242          WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_icebfr       = ', rn_icebfr 
    243          WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lfrelax      = ', rn_lfrelax 
    244          WRITE(numout,*) '         isotropic tensile strength                          rn_tensile      = ', rn_tensile 
     243         WRITE(numout,*) '         fraction of ocean depth that ice must reach         rn_lf_depfra    = ', rn_lf_depfra 
     244         WRITE(numout,*) '         maximum bottom stress per unit area of contact      rn_lf_bfr       = ', rn_lf_bfr 
     245         WRITE(numout,*) '         relax time scale (s-1) to reach static friction     rn_lf_relax     = ', rn_lf_relax 
     246         WRITE(numout,*) '         isotropic tensile strength                          rn_lf_tensile   = ', rn_lf_tensile 
    245247         WRITE(numout,*) 
    246248      ENDIF 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv.F90

    r12489 r13553  
    8282         !                             !-----------------------! 
    8383         CALL ice_dyn_adv_umx( nn_UMx, kt, u_ice, v_ice, h_i, h_s, h_ip, & 
    84             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     84            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i ) 
    8585         !                             !-----------------------! 
    8686      CASE( np_advPRA )                ! PRATHER scheme        ! 
    8787         !                             !-----------------------! 
    8888         CALL ice_dyn_adv_pra(         kt, u_ice, v_ice, h_i, h_s, h_ip, & 
    89             &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, e_s, e_i ) 
     89            &                          ato_i, v_i, v_s, sv_i, oa_i, a_i, a_ip, v_ip, v_il, e_s, e_i ) 
    9090      END SELECT 
    9191 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv_pra.F90

    r13295 r13553  
    4444   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxap , syap , sxxap , syyap , sxyap    ! melt pond fraction 
    4545   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvp , syvp , sxxvp , syyvp , sxyvp    ! melt pond volume 
     46   REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:,:)   ::   sxvl , syvl , sxxvl , syyvl , sxyvl    ! melt pond lid volume 
    4647 
    4748   !! * Substitutions 
     
    5556 
    5657   SUBROUTINE ice_dyn_adv_pra(         kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    57       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     58      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    5859      !!---------------------------------------------------------------------- 
    5960      !!                **  routine ice_dyn_adv_pra  ** 
     
    8182      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond fraction 
    8283      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     84      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid thickness 
    8385      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8486      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
    8587      ! 
    86       INTEGER  ::   ji,jj, jk, jl, jt       ! dummy loop indices 
     88      INTEGER  ::   ji, jj, jk, jl, jt      ! dummy loop indices 
    8789      INTEGER  ::   icycle                  ! number of sub-timestep for the advection 
    8890      REAL(wp) ::   zdt                     !   -      - 
     
    9092      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
    9193      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx 
    92       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max 
     94      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     95      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     96      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    9397      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zarea 
    9498      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ice, z0snw, z0ai, z0smi, z0oi 
    95       REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z0ap , z0vp, z0vl 
    96100      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   z0es 
    97101      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   z0ei 
     
    100104      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_pra: Prather advection scheme' 
    101105      ! 
    102       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     106      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     107      ! thickness and salinity 
     108      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     109      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     110      END WHERE 
    103111      DO jl = 1, jpl 
    104112         DO_2D( 0, 0, 0, 0 ) 
     
    115123               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    116124               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
     125            zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
     126               &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
     127               &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
     128               &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
    117129         END_2D 
    118130      END DO 
    119       CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 
     131      CALL lbc_lnk_multi( 'icedyn_adv_pra', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
     132      ! 
     133      ! enthalpies 
     134      DO jk = 1, nlay_i 
     135         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     136         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     137         END WHERE 
     138      END DO 
     139      DO jk = 1, nlay_s 
     140         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     141         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     142         END WHERE 
     143      END DO 
     144      DO jl = 1, jpl 
     145         DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
     146            zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
     147               &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
     148               &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
     149               &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
     150         END_3D 
     151      END DO 
     152      DO jl = 1, jpl 
     153         DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
     154            zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
     155               &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
     156               &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
     157               &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
     158         END_3D 
     159      END DO 
     160      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
     161      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
     162      ! 
    120163      ! 
    121164      ! --- If ice drift is too fast, use  subtime steps for advection (CFL test for stability) --- ! 
     
    156199               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    157200            END DO 
    158             IF ( ln_pnd_H12 ) THEN 
    159                z0ap(:,:,jl)  = pa_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond fraction 
    160                z0vp(:,:,jl)  = pv_ip(:,:,jl) * e1e2t(:,:)     ! Melt pond volume 
     201            IF ( ln_pnd_LEV ) THEN 
     202               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
     203               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
     204               IF ( ln_pnd_lids ) THEN 
     205                  z0vl(:,:,jl) = pv_il(:,:,jl) * e1e2t(:,:)   ! Melt pond lid volume 
     206               ENDIF 
    161207            ENDIF 
    162208         END DO 
     
    189235            END DO 
    190236            ! 
    191             IF ( ln_pnd_H12 ) THEN 
     237            IF ( ln_pnd_LEV ) THEN 
    192238               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    193239               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    194240               CALL adv_x( zdt , zudy , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    195241               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )  
     242               IF ( ln_pnd_lids ) THEN 
     243                  CALL adv_x( zdt , zudy , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     244                  CALL adv_y( zdt , zvdx , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     245               ENDIF 
    196246            ENDIF 
    197247            !                                                               !--------------------------------------------! 
     
    220270                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    221271            END DO 
    222             IF ( ln_pnd_H12 ) THEN 
     272            IF ( ln_pnd_LEV ) THEN 
    223273               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    224274               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    225275               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp )    !--- melt pond volume 
    226276               CALL adv_x( zdt , zudy , 0._wp , zarea , z0vp , sxvp , sxxvp , syvp , syyvp , sxyvp ) 
    227             ENDIF 
     277               IF ( ln_pnd_lids ) THEN 
     278                  CALL adv_y( zdt , zvdx , 1._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl ) !--- melt pond lid volume 
     279                  CALL adv_x( zdt , zudy , 0._wp , zarea , z0vl , sxvl , sxxvl , syvl , syyvl , sxyvl )  
     280               ENDIF 
     281           ENDIF 
    228282            ! 
    229283         ENDIF 
     
    242296               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    243297            END DO 
    244             IF ( ln_pnd_H12 ) THEN 
     298            IF ( ln_pnd_LEV ) THEN 
    245299               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    246300               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     301               IF ( ln_pnd_lids ) THEN 
     302                  pv_il(:,:,jl) = z0vl(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
     303               ENDIF 
    247304            ENDIF 
    248305         END DO 
     
    259316         !     Remove negative values (conservation is ensured) 
    260317         !     (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    261          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     318         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    262319         ! 
    263320         ! --- Make sure ice thickness is not too big --- ! 
    264321         !     (because ice thickness can be too large where ice concentration is very small) 
    265          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     322         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     323            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    266324         ! 
    267325         ! --- Ensure snow load is not too big --- ! 
     
    325383 
    326384         !  Calculate fluxes and moments between boxes i<-->i+1               
    327          DO_2D( 0, 0, 1, 1 ) 
     385         DO_2D( 0, 0, 1, 1 )                   !  Flux from i to i+1 WHEN u GT 0 
    328386            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, put(ji,jj) ) ) 
    329387            zalf         =  MAX( 0._wp, put(ji,jj) ) * pdt / psm(ji,jj,jl) 
     
    350408         END_2D 
    351409 
    352          DO_2D( 0, 0, 1, 0 ) 
     410         DO_2D( 0, 0, 1, 0 )                   !  Flux from i+1 to i when u LT 0. 
    353411            zalf          = MAX( 0._wp, -put(ji,jj) ) * pdt / psm(ji+1,jj,jl)  
    354412            zalg  (ji,jj) = zalf 
     
    369427         END_2D 
    370428 
    371          DO_2D( 0, 0, 0, 0 ) 
     429         DO_2D( 0, 0, 0, 0 )                   !  Readjust moments remaining in the box. 
    372430            zbt  =       zbet(ji-1,jj) 
    373431            zbt1 = 1.0 - zbet(ji-1,jj) 
     
    383441 
    384442         !   Put the temporary moments into appropriate neighboring boxes.     
    385          DO_2D( 0, 0, 0, 0 ) 
     443         DO_2D( 0, 0, 0, 0 )                   !   Flux from i to i+1 IF u GT 0. 
    386444            zbt  =       zbet(ji-1,jj) 
    387445            zbt1 = 1.0 - zbet(ji-1,jj) 
     
    403461         END_2D 
    404462 
    405          DO_2D( 0, 0, 0, 0 ) 
     463         DO_2D( 0, 0, 0, 0 )                   !  Flux from i+1 to i IF u LT 0. 
    406464            zbt  =       zbet(ji,jj) 
    407465            zbt1 = 1.0 - zbet(ji,jj) 
     
    482540  
    483541         !  Calculate fluxes and moments between boxes j<-->j+1               
    484          DO_2D( 1, 1, 0, 0 ) 
     542         DO_2D( 1, 1, 0, 0 )                   !  Flux from j to j+1 WHEN v GT 0 
    485543            zbet(ji,jj)  =  MAX( 0._wp, SIGN( 1._wp, pvt(ji,jj) ) ) 
    486544            zalf         =  MAX( 0._wp, pvt(ji,jj) ) * pdt / psm(ji,jj,jl) 
     
    507565         END_2D 
    508566         ! 
    509          DO_2D( 1, 0, 0, 0 ) 
     567         DO_2D( 1, 0, 0, 0 )                   !  Flux from j+1 to j when v LT 0. 
    510568            zalf          = MAX( 0._wp, -pvt(ji,jj) ) * pdt / psm(ji,jj+1,jl)  
    511569            zalg  (ji,jj) = zalf 
     
    541599 
    542600         !   Put the temporary moments into appropriate neighboring boxes.     
    543          DO_2D( 0, 0, 0, 0 ) 
     601         DO_2D( 0, 0, 0, 0 )                   !  Flux from j to j+1 IF v GT 0. 
    544602            zbt  =       zbet(ji,jj-1) 
    545603            zbt1 = 1.0 - zbet(ji,jj-1) 
     
    562620         END_2D 
    563621 
    564          DO_2D( 0, 0, 0, 0 ) 
     622         DO_2D( 0, 0, 0, 0 )                   !  Flux from j+1 to j IF v LT 0. 
    565623            zbt  =       zbet(ji,jj) 
    566624            zbt1 = 1.0 - zbet(ji,jj) 
     
    591649 
    592650 
    593    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     651   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     652      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    594653      !!------------------------------------------------------------------- 
    595654      !!                  ***  ROUTINE Hbig  *** 
     
    605664      !! ** input   : Max thickness of the surrounding 9-points 
    606665      !!------------------------------------------------------------------- 
    607       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    608       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    609       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     666      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     667      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     668      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     669      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     670      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    610671      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    611       ! 
    612       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    613       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     672      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     673      ! 
     674      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     675      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    614676      !!------------------------------------------------------------------- 
    615677      ! 
     
    617679      ! 
    618680      DO jl = 1, jpl 
    619  
    620681         DO_2D( 1, 1, 1, 1 ) 
    621682            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     
    623684               !                               ! -- check h_ip -- ! 
    624685               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    625                IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     686               IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    626687                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    627688                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     
    650711               ENDIF            
    651712               !                   
     713               !                               ! -- check s_i -- ! 
     714               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     715               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     716               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     717                  zfra = psi_max(ji,jj,jl) / zsi 
     718                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     719                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     720               ENDIF 
     721               ! 
    652722            ENDIF 
    653723         END_2D 
    654724      END DO  
     725      ! 
     726      !                                           ! -- check e_i/v_i -- ! 
     727      DO jl = 1, jpl 
     728         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     729            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     730               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     731               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     732               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     733                  zfra = pei_max(ji,jj,jk,jl) / zei 
     734                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     735                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     736               ENDIF 
     737            ENDIF 
     738         END_3D 
     739      END DO 
     740      !                                           ! -- check e_s/v_s -- ! 
     741      DO jl = 1, jpl 
     742         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
     743            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     744               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     745               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     746               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     747                  zfra = pes_max(ji,jj,jk,jl) / zes 
     748                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     749                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     750               ENDIF 
     751            ENDIF 
     752         END_3D 
     753      END DO 
    655754      ! 
    656755   END SUBROUTINE Hbig 
     
    724823         &      sxsal(jpi,jpj,jpl) , sysal(jpi,jpj,jpl) , sxxsal(jpi,jpj,jpl) , syysal(jpi,jpj,jpl) , sxysal(jpi,jpj,jpl) ,   & 
    725824         &      sxage(jpi,jpj,jpl) , syage(jpi,jpj,jpl) , sxxage(jpi,jpj,jpl) , syyage(jpi,jpj,jpl) , sxyage(jpi,jpj,jpl) ,   & 
    726          &      sxap(jpi,jpj,jpl)  , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
    727          &      sxvp(jpi,jpj,jpl)  , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     825         &      sxap (jpi,jpj,jpl) , syap (jpi,jpj,jpl) , sxxap (jpi,jpj,jpl) , syyap (jpi,jpj,jpl) , sxyap (jpi,jpj,jpl) ,   & 
     826         &      sxvp (jpi,jpj,jpl) , syvp (jpi,jpj,jpl) , sxxvp (jpi,jpj,jpl) , syyvp (jpi,jpj,jpl) , sxyvp (jpi,jpj,jpl) ,   & 
     827         &      sxvl (jpi,jpj,jpl) , syvl (jpi,jpj,jpl) , sxxvl (jpi,jpj,jpl) , syyvl (jpi,jpj,jpl) , sxyvl (jpi,jpj,jpl) ,   & 
    728828         ! 
    729829         &      sxc0 (jpi,jpj,nlay_s,jpl) , syc0 (jpi,jpj,nlay_s,jpl) , sxxc0(jpi,jpj,nlay_s,jpl) , & 
     
    820920            END DO 
    821921            ! 
    822             IF( ln_pnd_H12 ) THEN                                    ! melt pond fraction 
    823                CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap  ) 
    824                CALL iom_get( numrir, jpdom_auto, 'syap' , syap  ) 
    825                CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
    826                CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
    827                CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
    828                !                                                     ! melt pond volume 
    829                CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp  ) 
    830                CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp  ) 
    831                CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
    832                CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
    833                CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 
     922            IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
     923               IF( iom_varid( numror, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
     924                  CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap  ) 
     925                  CALL iom_get( numrir, jpdom_auto, 'syap' , syap  ) 
     926                  CALL iom_get( numrir, jpdom_auto, 'sxxap', sxxap ) 
     927                  CALL iom_get( numrir, jpdom_auto, 'syyap', syyap ) 
     928                  CALL iom_get( numrir, jpdom_auto, 'sxyap', sxyap ) 
     929                  !                                                     ! melt pond volume 
     930                  CALL iom_get( numrir, jpdom_auto, 'sxvp' , sxvp  ) 
     931                  CALL iom_get( numrir, jpdom_auto, 'syvp' , syvp  ) 
     932                  CALL iom_get( numrir, jpdom_auto, 'sxxvp', sxxvp ) 
     933                  CALL iom_get( numrir, jpdom_auto, 'syyvp', syyvp ) 
     934                  CALL iom_get( numrir, jpdom_auto, 'sxyvp', sxyvp ) 
     935               ELSE 
     936                  sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp   ! melt pond fraction 
     937                  sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp   ! melt pond volume 
     938               ENDIF 
     939                  ! 
     940               IF ( ln_pnd_lids ) THEN                               ! melt pond lid volume 
     941                  IF( iom_varid( numror, 'sxvl', ldstop = .FALSE. ) > 0 ) THEN 
     942                     CALL iom_get( numrir, jpdom_auto, 'sxvl' , sxvl  ) 
     943                     CALL iom_get( numrir, jpdom_auto, 'syvl' , syvl  ) 
     944                     CALL iom_get( numrir, jpdom_auto, 'sxxvl', sxxvl ) 
     945                     CALL iom_get( numrir, jpdom_auto, 'syyvl', syyvl ) 
     946                     CALL iom_get( numrir, jpdom_auto, 'sxyvl', sxyvl ) 
     947                  ELSE 
     948                     sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp   ! melt pond lid volume 
     949                  ENDIF 
     950               ENDIF 
    834951            ENDIF 
    835952            ! 
     
    845962            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    846963            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    847             IF( ln_pnd_H12 ) THEN 
    848                sxap  = 0._wp   ;   syap  = 0._wp   ;   sxxap  = 0._wp   ;   syyap  = 0._wp   ;   sxyap  = 0._wp   ! melt pond fraction 
    849                sxvp  = 0._wp   ;   syvp  = 0._wp   ;   sxxvp  = 0._wp   ;   syyvp  = 0._wp   ;   sxyvp  = 0._wp   ! melt pond volume 
     964            IF( ln_pnd_LEV ) THEN 
     965               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
     966               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
     967               IF ( ln_pnd_lids ) THEN 
     968                  sxvl = 0._wp; syvl = 0._wp    ;   sxxvl = 0._wp    ;   syyvl = 0._wp    ;   sxyvl = 0._wp       ! melt pond lid volume 
     969               ENDIF 
    850970            ENDIF 
    851971         ENDIF 
     
    9101030         END DO 
    9111031         ! 
    912          IF( ln_pnd_H12 ) THEN                                       ! melt pond fraction 
     1032         IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction 
    9131033            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    9141034            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 
     
    9221042            CALL iom_rstput( iter, nitrst, numriw, 'syyvp', syyvp ) 
    9231043            CALL iom_rstput( iter, nitrst, numriw, 'sxyvp', sxyvp ) 
     1044            ! 
     1045            IF ( ln_pnd_lids ) THEN                                  ! melt pond lid volume 
     1046               CALL iom_rstput( iter, nitrst, numriw, 'sxvl' , sxvl  ) 
     1047               CALL iom_rstput( iter, nitrst, numriw, 'syvl' , syvl  ) 
     1048               CALL iom_rstput( iter, nitrst, numriw, 'sxxvl', sxxvl ) 
     1049               CALL iom_rstput( iter, nitrst, numriw, 'syyvl', syyvl ) 
     1050               CALL iom_rstput( iter, nitrst, numriw, 'sxyvl', sxyvl ) 
     1051            ENDIF 
    9241052         ENDIF 
    9251053         ! 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_adv_umx.F90

    r13295 r13553  
    6060 
    6161   SUBROUTINE ice_dyn_adv_umx( kn_umx, kt, pu_ice, pv_ice, ph_i, ph_s, ph_ip,  & 
    62       &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     62      &                        pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    6363      !!---------------------------------------------------------------------- 
    6464      !!                  ***  ROUTINE ice_dyn_adv_umx  *** 
     
    8585      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pa_ip      ! melt pond concentration 
    8686      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_ip      ! melt pond volume 
     87      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_il      ! melt pond lid volume 
    8788      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s       ! snw heat content 
    8889      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i       ! ice heat content 
     
    9293      REAL(wp) ::   zamsk                   ! 1 if advection of concentration, 0 if advection of other tracers 
    9394      REAL(wp) ::   zdt, zvi_cen 
    94       REAL(wp), DIMENSION(1)           ::   zcflprv, zcflnow   ! for global communication 
    95       REAL(wp), DIMENSION(jpi,jpj)     ::   zudy, zvdx, zcu_box, zcv_box 
    96       REAL(wp), DIMENSION(jpi,jpj)     ::   zati1, zati2 
    97       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zu_cat, zv_cat 
    98       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zua_ho, zva_ho, zua_ups, zva_ups 
    99       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_ai , z1_aip, zhvar 
    100       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zhi_max, zhs_max, zhip_max 
     95      REAL(wp), DIMENSION(1)                  ::   zcflprv, zcflnow   ! for global communication 
     96      REAL(wp), DIMENSION(jpi,jpj)            ::   zudy, zvdx, zcu_box, zcv_box 
     97      REAL(wp), DIMENSION(jpi,jpj)            ::   zati1, zati2 
     98      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zu_cat, zv_cat 
     99      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zua_ho, zva_ho, zua_ups, zva_ups 
     100      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   z1_ai , z1_aip, zhvar 
     101      REAL(wp), DIMENSION(jpi,jpj,jpl)        ::   zhi_max, zhs_max, zhip_max, zs_i, zsi_max 
     102      REAL(wp), DIMENSION(jpi,jpj,nlay_i,jpl) ::   ze_i, zei_max 
     103      REAL(wp), DIMENSION(jpi,jpj,nlay_s,jpl) ::   ze_s, zes_max 
    101104      ! 
    102105      REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) ::   zuv_ho, zvv_ho, zuv_ups, zvv_ups, z1_vi, z1_vs  
     
    105108      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_adv_umx: Ultimate-Macho advection scheme' 
    106109      ! 
    107       ! --- Record max of the surrounding 9-pts ice thick. (for call Hbig) --- ! 
     110      ! --- Record max of the surrounding 9-pts (for call Hbig) --- ! 
     111      ! thickness and salinity 
     112      WHERE( pv_i(:,:,:) >= epsi10 ) ; zs_i(:,:,:) = psv_i(:,:,:) / pv_i(:,:,:) 
     113      ELSEWHERE                      ; zs_i(:,:,:) = 0._wp 
     114      END WHERE 
    108115      DO jl = 1, jpl 
    109116         DO_2D( 0, 0, 0, 0 ) 
     
    120127               &                                               ph_s (ji+1,jj+1,jl), ph_s (ji-1,jj-1,jl), & 
    121128               &                                               ph_s (ji+1,jj-1,jl), ph_s (ji-1,jj+1,jl) ) 
    122          END_2D 
    123       END DO 
    124       CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1.0_wp, zhs_max, 'T', 1.0_wp, zhip_max, 'T', 1.0_wp ) 
     129            zsi_max (ji,jj,jl) = MAX( epsi20, zs_i (ji,jj,jl), zs_i (ji+1,jj  ,jl), zs_i (ji  ,jj+1,jl), & 
     130               &                                               zs_i (ji-1,jj  ,jl), zs_i (ji  ,jj-1,jl), & 
     131               &                                               zs_i (ji+1,jj+1,jl), zs_i (ji-1,jj-1,jl), & 
     132               &                                               zs_i (ji+1,jj-1,jl), zs_i (ji-1,jj+1,jl) ) 
     133         END_2D 
     134      END DO 
     135      CALL lbc_lnk_multi( 'icedyn_adv_umx', zhi_max, 'T', 1._wp, zhs_max, 'T', 1._wp, zhip_max, 'T', 1._wp, zsi_max, 'T', 1._wp ) 
     136      ! 
     137      ! enthalpies 
     138      DO jk = 1, nlay_i 
     139         WHERE( pv_i(:,:,:) >= epsi10 ) ; ze_i(:,:,jk,:) = pe_i(:,:,jk,:) / pv_i(:,:,:) 
     140         ELSEWHERE                      ; ze_i(:,:,jk,:) = 0._wp 
     141         END WHERE 
     142      END DO 
     143      DO jk = 1, nlay_s 
     144         WHERE( pv_s(:,:,:) >= epsi10 ) ; ze_s(:,:,jk,:) = pe_s(:,:,jk,:) / pv_s(:,:,:) 
     145         ELSEWHERE                      ; ze_s(:,:,jk,:) = 0._wp 
     146         END WHERE 
     147      END DO 
     148      DO jl = 1, jpl 
     149         DO_3D( 0, 0, 0, 0, 1, nlay_i ) 
     150            zei_max(ji,jj,jk,jl) = MAX( epsi20, ze_i(ji,jj,jk,jl), ze_i(ji+1,jj  ,jk,jl), ze_i(ji  ,jj+1,jk,jl), & 
     151               &                                                   ze_i(ji-1,jj  ,jk,jl), ze_i(ji  ,jj-1,jk,jl), & 
     152               &                                                   ze_i(ji+1,jj+1,jk,jl), ze_i(ji-1,jj-1,jk,jl), & 
     153               &                                                   ze_i(ji+1,jj-1,jk,jl), ze_i(ji-1,jj+1,jk,jl) ) 
     154         END_3D 
     155      END DO 
     156      DO jl = 1, jpl 
     157         DO_3D( 0, 0, 0, 0, 1, nlay_s ) 
     158            zes_max(ji,jj,jk,jl) = MAX( epsi20, ze_s(ji,jj,jk,jl), ze_s(ji+1,jj  ,jk,jl), ze_s(ji  ,jj+1,jk,jl), & 
     159               &                                                   ze_s(ji-1,jj  ,jk,jl), ze_s(ji  ,jj-1,jk,jl), & 
     160               &                                                   ze_s(ji+1,jj+1,jk,jl), ze_s(ji-1,jj-1,jk,jl), & 
     161               &                                                   ze_s(ji+1,jj-1,jk,jl), ze_s(ji-1,jj+1,jk,jl) ) 
     162         END_3D 
     163      END DO 
     164      CALL lbc_lnk( 'icedyn_adv_pra', zei_max, 'T', 1. ) 
     165      CALL lbc_lnk( 'icedyn_adv_pra', zes_max, 'T', 1. ) 
    125166      ! 
    126167      ! 
     
    318359         ! 
    319360         !== melt ponds ==! 
    320          IF ( ln_pnd_H12 ) THEN 
     361         IF ( ln_pnd_LEV ) THEN 
    321362            ! concentration 
    322363            zamsk = 1._wp 
     
    328369            CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
    329370               &                                      zhvar, pv_ip, zua_ups, zva_ups ) 
     371            ! lid 
     372            IF ( ln_pnd_lids ) THEN 
     373               zamsk = 0._wp 
     374               zhvar(:,:,:) = pv_il(:,:,:) * z1_aip(:,:,:) 
     375               CALL adv_umx( zamsk, kn_umx, jt, kt, zdt, zudy , zvdx , zua_ho , zva_ho , zcu_box, zcv_box, & 
     376                  &                                      zhvar, pv_il, zua_ups, zva_ups ) 
     377            ENDIF 
    330378         ENDIF 
    331379         ! 
     
    342390         ! Remove negative values (conservation is ensured) 
    343391         !    (because advected fields are not perfectly bounded and tiny negative values can occur, e.g. -1.e-20) 
    344          CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pe_s, pe_i ) 
     392         CALL ice_var_zapneg( zdt, pato_i, pv_i, pv_s, psv_i, poa_i, pa_i, pa_ip, pv_ip, pv_il, pe_s, pe_i ) 
    345393         ! 
    346394         ! --- Make sure ice thickness is not too big --- ! 
    347395         !     (because ice thickness can be too large where ice concentration is very small) 
    348          CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     396         CALL Hbig( zdt, zhi_max, zhs_max, zhip_max, zsi_max, zes_max, zei_max, & 
     397            &            pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    349398         ! 
    350399         ! --- Ensure snow load is not too big --- ! 
     
    9571006      !                                                     !--  Laplacian in j-direction  --! 
    9581007      DO jl = 1, jpl 
    959          DO_2D( 1, 0, 0, 0 ) 
     1008         DO_2D( 1, 0, 0, 0 )         ! First derivative (gradient) 
    9601009            ztv1(ji,jj,jl) = ( pt(ji,jj+1,jl) - pt(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    9611010         END_2D 
    962          DO_2D( 0, 0, 0, 0 ) 
     1011         DO_2D( 0, 0, 0, 0 )         ! Second derivative (Laplacian) 
    9631012            ztv2(ji,jj,jl) = ( ztv1(ji,jj,jl) - ztv1(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    9641013         END_2D 
     
    9681017      !                                                     !--  BiLaplacian in j-direction  --! 
    9691018      DO jl = 1, jpl 
    970          DO_2D( 1, 0, 0, 0 ) 
     1019         DO_2D( 1, 0, 0, 0 )         ! First derivative 
    9711020            ztv3(ji,jj,jl) = ( ztv2(ji,jj+1,jl) - ztv2(ji,jj,jl) ) * r1_e2v(ji,jj) * vmask(ji,jj,1) 
    9721021         END_2D 
    973          DO_2D( 0, 0, 0, 0 ) 
     1022         DO_2D( 0, 0, 0, 0 )         ! Second derivative 
    9741023            ztv4(ji,jj,jl) = ( ztv3(ji,jj,jl) - ztv3(ji,jj-1,jl) ) * r1_e2t(ji,jj) 
    9751024         END_2D 
     
    14091458 
    14101459 
    1411    SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, pv_i, pv_s, pa_i, pa_ip, pv_ip, pe_s ) 
     1460   SUBROUTINE Hbig( pdt, phi_max, phs_max, phip_max, psi_max, pes_max, pei_max, & 
     1461      &                  pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i, pe_s, pe_i ) 
    14121462      !!------------------------------------------------------------------- 
    14131463      !!                  ***  ROUTINE Hbig  *** 
     
    14231473      !! ** input   : Max thickness of the surrounding 9-points 
    14241474      !!------------------------------------------------------------------- 
    1425       REAL(wp)                    , INTENT(in   ) ::   pdt                          ! tracer time-step 
    1426       REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max   ! max ice thick from surrounding 9-pts 
    1427       REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip 
     1475      REAL(wp)                    , INTENT(in   ) ::   pdt                                   ! tracer time-step 
     1476      REAL(wp), DIMENSION(:,:,:)  , INTENT(in   ) ::   phi_max, phs_max, phip_max, psi_max   ! max ice thick from surrounding 9-pts 
     1477      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pes_max 
     1478      REAL(wp), DIMENSION(:,:,:,:), INTENT(in   ) ::   pei_max 
     1479      REAL(wp), DIMENSION(:,:,:)  , INTENT(inout) ::   pv_i, pv_s, pa_i, pa_ip, pv_ip, psv_i 
    14281480      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_s 
    1429       ! 
    1430       INTEGER  ::   ji, jj, jl         ! dummy loop indices 
    1431       REAL(wp) ::   z1_dt, zhip, zhi, zhs, zfra 
     1481      REAL(wp), DIMENSION(:,:,:,:), INTENT(inout) ::   pe_i 
     1482      ! 
     1483      INTEGER  ::   ji, jj, jk, jl         ! dummy loop indices 
     1484      REAL(wp) ::   z1_dt, zhip, zhi, zhs, zsi, zes, zei, zfra 
    14321485      !!------------------------------------------------------------------- 
    14331486      ! 
     
    14351488      ! 
    14361489      DO jl = 1, jpl 
    1437  
    14381490         DO_2D( 1, 1, 1, 1 ) 
    14391491            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     
    14411493               !                               ! -- check h_ip -- ! 
    14421494               ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1443                IF( ln_pnd_H12 .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1495               IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    14441496                  zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    14451497                  IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
     
    14681520               ENDIF            
    14691521               !                   
     1522               !                               ! -- check s_i -- ! 
     1523               ! if s_i is larger than the surrounding 9 pts => put salt excess in the ocean 
     1524               zsi = psv_i(ji,jj,jl) / pv_i(ji,jj,jl) 
     1525               IF( zsi > psi_max(ji,jj,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1526                  zfra = psi_max(ji,jj,jl) / zsi 
     1527                  sfx_res(ji,jj) = sfx_res(ji,jj) + psv_i(ji,jj,jl) * ( 1._wp - zfra ) * rhoi * z1_dt 
     1528                  psv_i(ji,jj,jl) = psv_i(ji,jj,jl) * zfra 
     1529               ENDIF 
     1530               ! 
    14701531            ENDIF 
    14711532         END_2D 
    14721533      END DO  
     1534      ! 
     1535      !                                           ! -- check e_i/v_i -- ! 
     1536      DO jl = 1, jpl 
     1537         DO_3D( 1, 1, 1, 1, 1, nlay_i ) 
     1538            IF ( pv_i(ji,jj,jl) > 0._wp ) THEN 
     1539               ! if e_i/v_i is larger than the surrounding 9 pts => put the heat excess in the ocean 
     1540               zei = pe_i(ji,jj,jk,jl) / pv_i(ji,jj,jl) 
     1541               IF( zei > pei_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1542                  zfra = pei_max(ji,jj,jk,jl) / zei 
     1543                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_i(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1544                  pe_i(ji,jj,jk,jl) = pe_i(ji,jj,jk,jl) * zfra 
     1545               ENDIF 
     1546            ENDIF 
     1547         END_3D 
     1548      END DO 
     1549      !                                           ! -- check e_s/v_s -- ! 
     1550      DO jl = 1, jpl 
     1551         DO_3D( 1, 1, 1, 1, 1, nlay_s ) 
     1552            IF ( pv_s(ji,jj,jl) > 0._wp ) THEN 
     1553               ! if e_s/v_s is larger than the surrounding 9 pts => put the heat excess in the ocean 
     1554               zes = pe_s(ji,jj,jk,jl) / pv_s(ji,jj,jl) 
     1555               IF( zes > pes_max(ji,jj,jk,jl) .AND. pa_i(ji,jj,jl) < 0.15 ) THEN 
     1556                  zfra = pes_max(ji,jj,jk,jl) / zes 
     1557                  hfx_res(ji,jj) = hfx_res(ji,jj) - pe_s(ji,jj,jk,jl) * ( 1._wp - zfra ) * z1_dt ! W.m-2 <0 
     1558                  pe_s(ji,jj,jk,jl) = pe_s(ji,jj,jk,jl) * zfra 
     1559               ENDIF 
     1560            ENDIF 
     1561         END_3D 
     1562      END DO 
    14731563      ! 
    14741564   END SUBROUTINE Hbig 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_rdgrft.F90

    r13295 r13553  
    502502      REAL(wp)                  ::   airdg1, oirdg1, aprdg1, virdg1, sirdg1 
    503503      REAL(wp)                  ::   airft1, oirft1, aprft1 
    504       REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg  ! area etc of new ridges 
    505       REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft  ! area etc of rafted ice 
     504      REAL(wp), DIMENSION(jpij) ::   airdg2, oirdg2, aprdg2, virdg2, sirdg2, vsrdg, vprdg, vlrdg  ! area etc of new ridges 
     505      REAL(wp), DIMENSION(jpij) ::   airft2, oirft2, aprft2, virft , sirft , vsrft, vprft, vlrft  ! area etc of rafted ice 
    506506      ! 
    507507      REAL(wp), DIMENSION(jpij) ::   ersw             ! enth of water trapped into ridges 
     
    530530      DO jl1 = 1, jpl 
    531531 
    532          CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) ) 
     532         IF( nn_icesal /= 2 )  THEN       
     533            CALL tab_2d_1d( npti, nptidx(1:npti), s_i_1d(1:npti), s_i(:,:,jl1) ) 
     534         ENDIF 
    533535 
    534536         DO ji = 1, npti 
     
    573575               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft  
    574576 
    575                IF ( ln_pnd_H12 ) THEN 
     577               IF ( ln_pnd_LEV ) THEN 
    576578                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
    577579                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
     
    580582                  aprft2(ji) = a_ip_2d(ji,jl1) * afrft * hi_hrft 
    581583                  vprft (ji) = v_ip_2d(ji,jl1) * afrft 
     584                  IF ( ln_pnd_lids ) THEN 
     585                     vlrdg (ji) = v_il_2d(ji,jl1) * afrdg 
     586                     vlrft (ji) = v_il_2d(ji,jl1) * afrft 
     587                  ENDIF 
    582588               ENDIF 
    583589 
     
    606612               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji) 
    607613               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1 
    608                IF ( ln_pnd_H12 ) THEN 
     614               IF ( ln_pnd_LEV ) THEN 
    609615                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    610616                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
     617                  IF ( ln_pnd_lids ) THEN 
     618                     v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - vlrdg(ji) - vlrft(ji) 
     619                  ENDIF 
    611620               ENDIF 
    612621            ENDIF 
     
    700709                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  & 
    701710                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
    702                   IF ( ln_pnd_H12 ) THEN 
     711                  IF ( ln_pnd_LEV ) THEN 
    703712                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    704713                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 
    705714                     a_ip_2d (ji,jl2) = a_ip_2d(ji,jl2) + (   aprdg2(ji) * rn_fpndrdg * farea         &  
    706715                        &                                   + aprft2(ji) * rn_fpndrft * zswitch(ji)   ) 
     716                     IF ( ln_pnd_lids ) THEN 
     717                        v_il_2d (ji,jl2) = v_il_2d(ji,jl2) + (   vlrdg(ji) * rn_fpndrdg * fvol   (ji) & 
     718                           &                                   + vlrft(ji) * rn_fpndrft * zswitch(ji) ) 
     719                     ENDIF 
    707720                  ENDIF 
    708721                   
     
    735748      !---------------- 
    736749      ! In case ridging/rafting lead to very small negative values (sometimes it happens) 
    737       CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, ze_s_2d, ze_i_2d ) 
     750      CALL ice_var_roundoff( a_i_2d, v_i_2d, v_s_2d, sv_i_2d, oa_i_2d, a_ip_2d, v_ip_2d, v_il_2d, ze_s_2d, ze_i_2d ) 
    738751      ! 
    739752   END SUBROUTINE rdgrft_shift 
     
    841854         CALL tab_3d_2d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    842855         CALL tab_3d_2d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     856         CALL tab_3d_2d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) ) 
    843857         DO jl = 1, jpl 
    844858            DO jk = 1, nlay_s 
     
    867881         CALL tab_2d_3d( npti, nptidx(1:npti), a_ip_2d(1:npti,1:jpl), a_ip(:,:,:) ) 
    868882         CALL tab_2d_3d( npti, nptidx(1:npti), v_ip_2d(1:npti,1:jpl), v_ip(:,:,:) ) 
     883         CALL tab_2d_3d( npti, nptidx(1:npti), v_il_2d(1:npti,1:jpl), v_il(:,:,:) ) 
    869884         DO jl = 1, jpl 
    870885            DO jk = 1, nlay_s 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_rhg.F90

    r12377 r13553  
    108108      INTEGER ::   ios, ioptio   ! Local integer output status for namelist read 
    109109      !! 
    110       NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast 
     110      NAMELIST/namdyn_rhg/  ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg 
    111111      !!------------------------------------------------------------------- 
    112112      ! 
     
    122122         WRITE(numout,*) '~~~~~~~~~~~~~~~' 
    123123         WRITE(numout,*) '   Namelist : namdyn_rhg:' 
    124          WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP = ', ln_rhg_EVP 
    125          WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP    = ', ln_aEVP 
    126          WRITE(numout,*) '         creep limit                                       rn_creepl  = ', rn_creepl 
    127          WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc     = ', rn_ecc 
    128          WRITE(numout,*) '         number of iterations for subcycling               nn_nevp    = ', nn_nevp 
    129          WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast  = ', rn_relast 
     124         WRITE(numout,*) '      rheology EVP (icedyn_rhg_evp)                        ln_rhg_EVP    = ', ln_rhg_EVP 
     125         WRITE(numout,*) '         use adaptive EVP (aEVP)                           ln_aEVP       = ', ln_aEVP 
     126         WRITE(numout,*) '         creep limit                                       rn_creepl     = ', rn_creepl 
     127         WRITE(numout,*) '         eccentricity of the elliptical yield curve        rn_ecc        = ', rn_ecc 
     128         WRITE(numout,*) '         number of iterations for subcycling               nn_nevp       = ', nn_nevp 
     129         WRITE(numout,*) '         ratio of elastic timescale over ice time step     rn_relast     = ', rn_relast 
     130         WRITE(numout,*) '      check convergence of rheology                        nn_rhg_chkcvg = ', nn_rhg_chkcvg 
     131         IF    ( nn_rhg_chkcvg == 0 ) THEN   ;   WRITE(numout,*) '         no check' 
     132         ELSEIF( nn_rhg_chkcvg == 1 ) THEN   ;   WRITE(numout,*) '         check cvg at the main time step' 
     133         ELSEIF( nn_rhg_chkcvg == 2 ) THEN   ;   WRITE(numout,*) '         check cvg at both main and rheology time steps' 
     134         ENDIF 
    130135      ENDIF 
    131136      ! 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/icedyn_rhg_evp.F90

    r13295 r13553  
    4141   USE prtctl         ! Print control 
    4242 
     43   USE netcdf         ! NetCDF library for convergence test 
    4344   IMPLICIT NONE 
    4445   PRIVATE 
     
    5051#  include "do_loop_substitute.h90" 
    5152#  include "domzgr_substitute.h90" 
     53 
     54   !! for convergence tests 
     55   INTEGER ::   ncvgid   ! netcdf file id 
     56   INTEGER ::   nvarid   ! netcdf variable id 
     57   REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zmsk00, zmsk15 
    5258   !!---------------------------------------------------------------------- 
    5359   !! NEMO/ICE 4.0 , NEMO Consortium (2018) 
     
    121127      REAL(wp) ::   ecc2, z1_ecc2                                       ! square of yield ellipse eccenticity 
    122128      REAL(wp) ::   zalph1, z1_alph1, zalph2, z1_alph2                  ! alpha coef from Bouillon 2009 or Kimmritz 2017 
     129      REAl(wp) ::   zbetau, zbetav 
    123130      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV, zvU, zvV             ! ice/snow mass and volume 
    124       REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2       ! temporary scalars 
     131      REAL(wp) ::   zp_delf, zds2, zdt, zdt2, zdiv, zdiv2               ! temporary scalars 
    125132      REAL(wp) ::   zTauO, zTauB, zRHS, zvel                            ! temporary scalars 
    126133      REAL(wp) ::   zkt                                                 ! isotropic tensile strength for landfast ice 
    127134      REAL(wp) ::   zvCr                                                ! critical ice volume above which ice is landfast 
    128135      ! 
    129       REAL(wp) ::   zresm                                               ! Maximal error on ice velocity 
    130136      REAL(wp) ::   zintb, zintn                                        ! dummy argument 
    131137      REAL(wp) ::   zfac_x, zfac_y 
    132138      REAL(wp) ::   zshear, zdum1, zdum2 
    133139      ! 
    134       REAL(wp), DIMENSION(jpi,jpj) ::   zp_delt                         ! P/delta at T points 
     140      REAL(wp), DIMENSION(jpi,jpj) ::   zdelta, zp_delt                 ! delta and P/delta at T points 
    135141      REAL(wp), DIMENSION(jpi,jpj) ::   zbeta                           ! beta coef from Kimmritz 2017 
    136142      ! 
     
    139145      REAL(wp), DIMENSION(jpi,jpj) ::   zmU_t, zmV_t                    ! (ice-snow_mass / dt) on U/V points 
    140146      REAL(wp), DIMENSION(jpi,jpj) ::   zmf                             ! coriolis parameter at T points 
    141       REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     147      REAL(wp), DIMENSION(jpi,jpj) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points 
    142148      ! 
    143149      REAL(wp), DIMENSION(jpi,jpj) ::   zds                             ! shear 
    144150      REAL(wp), DIMENSION(jpi,jpj) ::   zs1, zs2, zs12                  ! stress tensor components 
    145 !!$      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice, zresr           ! check convergence 
    146151      REAL(wp), DIMENSION(jpi,jpj) ::   zsshdyn                         ! array used for the calculation of ice surface slope: 
    147152      !                                                                 !    ocean surface (ssh_m) if ice is not embedded 
     
    157162      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk01x, zmsk01y                ! dummy arrays 
    158163      REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00x, zmsk00y                ! mask for ice presence 
    159       REAL(wp), DIMENSION(jpi,jpj) ::   zfmask, zwf                     ! mask at F points for the ice 
     164      REAL(wp), DIMENSION(jpi,jpj) ::   zfmask                          ! mask at F points for the ice 
    160165 
    161166      REAL(wp), PARAMETER          ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
    162167      REAL(wp), PARAMETER          ::   zmmin  = 1._wp                  ! ice mass (kg/m2)  below which ice velocity becomes very small 
    163168      REAL(wp), PARAMETER          ::   zamin  = 0.001_wp               ! ice concentration below which ice velocity becomes very small 
     169      !! --- check convergence 
     170      REAL(wp), DIMENSION(jpi,jpj) ::   zu_ice, zv_ice 
    164171      !! --- diags 
    165       REAL(wp), DIMENSION(jpi,jpj) ::   zmsk00 
    166172      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zsig1, zsig2, zsig3 
    167173      !! --- SIMIP diags 
     
    176182      IF( kt == nit000 .AND. lwp )   WRITE(numout,*) '-- ice_dyn_rhg_evp: EVP sea-ice rheology' 
    177183      ! 
    178 !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
     184      ! for diagnostics and convergence tests 
     185      ALLOCATE( zmsk00(jpi,jpj), zmsk15(jpi,jpj) ) 
     186      DO_2D( 1, 1, 1, 1 ) 
     187         zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06  ) ) ! 1 if ice    , 0 if no ice 
     188         zmsk15(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.15_wp ) ) ! 1 if 15% ice, 0 if less 
     189      END_2D 
     190      ! 
     191      !!gm for Clem:  OPTIMIZATION:  I think zfmask can be computed one for all at the initialization.... 
    179192      !------------------------------------------------------------------------------! 
    180193      ! 0) mask at F points for the ice 
     
    187200 
    188201      ! Lateral boundary conditions on velocity (modify zfmask) 
    189       zwf(:,:) = zfmask(:,:) 
    190202      DO_2D( 0, 0, 0, 0 ) 
    191203         IF( zfmask(ji,jj) == 0._wp ) THEN 
    192             zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
     204            zfmask(ji,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(ji,jj,1), umask(ji,jj+1,1), & 
     205               &                                          vmask(ji,jj,1), vmask(ji+1,jj,1) ) ) 
    193206         ENDIF 
    194207      END_2D 
    195208      DO jj = 2, jpjm1 
    196209         IF( zfmask(1,jj) == 0._wp ) THEN 
    197             zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     210            zfmask(1  ,jj) = rn_ishlat * MIN( 1._wp , MAX( vmask(2,jj,1), umask(1,jj+1,1), umask(1,jj,1) ) ) 
    198211         ENDIF 
    199212         IF( zfmask(jpi,jj) == 0._wp ) THEN 
    200             zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
    201          ENDIF 
     213            zfmask(jpi,jj) = rn_ishlat * MIN( 1._wp , MAX( umask(jpi,jj+1,1), vmask(jpim1,jj,1), umask(jpi,jj-1,1) ) ) 
     214        ENDIF 
    202215      END DO 
    203216      DO ji = 2, jpim1 
    204217         IF( zfmask(ji,1) == 0._wp ) THEN 
    205             zfmask(ji,1  ) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     218            zfmask(ji, 1 ) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,1,1), umask(ji,2,1), vmask(ji,1,1) ) ) 
    206219         ENDIF 
    207220         IF( zfmask(ji,jpj) == 0._wp ) THEN 
    208             zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     221            zfmask(ji,jpj) = rn_ishlat * MIN( 1._wp , MAX( vmask(ji+1,jpj,1), vmask(ji-1,jpj,1), umask(ji,jpjm1,1) ) ) 
    209222         ENDIF 
    210223      END DO 
     
    220233      z1_ecc2 = 1._wp / ecc2 
    221234 
    222       ! Time step for subcycling 
    223       zdtevp   = rDt_ice / REAL( nn_nevp ) 
    224       z1_dtevp = 1._wp / zdtevp 
    225  
    226235      ! alpha parameters (Bouillon 2009) 
    227236      IF( .NOT. ln_aEVP ) THEN 
    228          zalph1 = ( 2._wp * rn_relast * rDt_ice ) * z1_dtevp 
     237         zdtevp   = rDt_ice / REAL( nn_nevp ) 
     238         zalph1 =   2._wp * rn_relast * REAL( nn_nevp ) 
    229239         zalph2 = zalph1 * z1_ecc2 
    230240 
    231241         z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
    232242         z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     243      ELSE 
     244         zdtevp   = rdt_ice 
     245         ! zalpha parameters set later on adaptatively 
    233246      ENDIF 
     247      z1_dtevp = 1._wp / zdtevp 
    234248          
    235249      ! Initialise stress tensor  
     
    242256 
    243257      ! landfast param from Lemieux(2016): add isotropic tensile strength (following Konig Beatty and Holland, 2010) 
    244       IF( ln_landfast_L16 ) THEN   ;   zkt = rn_tensile 
     258      IF( ln_landfast_L16 ) THEN   ;   zkt = rn_lf_tensile 
    245259      ELSE                         ;   zkt = 0._wp 
    246260      ENDIF 
     
    310324            zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 
    311325            ! ice-bottom stress at U points 
    312             zvCr = zaU(ji,jj) * rn_depfra * hu(ji,jj,Kmm) 
    313             ztaux_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
     326            zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 
     327            ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 
    314328            ! ice-bottom stress at V points 
    315             zvCr = zaV(ji,jj) * rn_depfra * hv(ji,jj,Kmm) 
    316             ztauy_base(ji,jj) = - rn_icebfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
     329            zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 
     330            ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 
    317331            ! ice_bottom stress at T points 
    318             zvCr = at_i(ji,jj) * rn_depfra * ht(ji,jj) 
    319             tau_icebfr(ji,jj) = - rn_icebfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
     332            zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 
     333            tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 
    320334         END_2D 
    321335         CALL lbc_lnk( 'icedyn_rhg_evp', tau_icebfr(:,:), 'T', 1.0_wp ) 
     
    337351         l_full_nf_update = jter == nn_nevp   ! false: disable full North fold update (performances) for iter = 1 to nn_nevp-1 
    338352         ! 
    339 !!$         IF(sn_cfctl%l_prtctl) THEN   ! Convergence test 
    340 !!$            DO jj = 1, jpjm1 
    341 !!$               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
    342 !!$               zv_ice(:,jj) = v_ice(:,jj) 
    343 !!$            END DO 
    344 !!$         ENDIF 
     353         ! convergence test 
     354         IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2  ) THEN 
     355            DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) 
     356               zu_ice(ji,jj) = u_ice(ji,jj) * umask(ji,jj,1) ! velocity at previous time step 
     357               zv_ice(ji,jj) = v_ice(ji,jj) * vmask(ji,jj,1) 
     358            END_2D 
     359         ENDIF 
    345360 
    346361         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     
    353368 
    354369         END_2D 
    355          CALL lbc_lnk( 'icedyn_rhg_evp', zds, 'F', 1.0_wp ) 
    356  
    357          DO_2D( 0, 1, 0, 1 ) 
     370 
     371         DO_2D( 0, 0, 0, 0 ) 
    358372 
    359373            ! shear**2 at T points (doc eq. A16) 
     
    375389             
    376390            ! delta at T points 
    377             zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
    378  
    379             ! P/delta at T points 
    380             zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 
    381  
    382             ! alpha & beta for aEVP 
     391            zdelta(ji,jj) = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 )   
     392 
     393         END_2D 
     394         CALL lbc_lnk( 'icedyn_rhg_evp', zdelta, 'T', 1.0_wp ) 
     395 
     396         ! P/delta at T points 
     397         DO_2D( 1, 1, 1, 1 ) 
     398            zp_delt(ji,jj) = strength(ji,jj) / ( zdelta(ji,jj) + rn_creepl ) 
     399         END_2D 
     400 
     401         DO_2D( 0, 1, 0, 1 )   ! loop ends at jpi,jpj so that no lbc_lnk are needed for zs1 and zs2 
     402 
     403            ! divergence at T points (duplication to avoid communications) 
     404            zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     405               &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     406               &    ) * r1_e1e2t(ji,jj) 
     407             
     408            ! tension at T points (duplication to avoid communications) 
     409            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     410               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     411               &   ) * r1_e1e2t(ji,jj) 
     412             
     413            ! alpha for aEVP 
    383414            !   gamma = 0.5*P/(delta+creepl) * (c*pi)**2/Area * dt/m 
    384415            !   alpha = beta = sqrt(4*gamma) 
     
    388419               zalph2   = zalph1 
    389420               z1_alph2 = z1_alph1 
     421               ! explicit: 
     422               ! z1_alph1 = 1._wp / zalph1 
     423               ! z1_alph2 = 1._wp / zalph1 
     424               ! zalph1 = zalph1 - 1._wp 
     425               ! zalph2 = zalph1 
    390426            ENDIF 
    391427             
    392428            ! stress at T points (zkt/=0 if landfast) 
    393             zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv * (1._wp + zkt) - zdelta * (1._wp - zkt) ) ) * z1_alph1 
    394             zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
     429            zs1(ji,jj) = ( zs1(ji,jj)*zalph1 + zp_delt(ji,jj) * ( zdiv*(1._wp + zkt) - zdelta(ji,jj)*(1._wp - zkt) ) ) * z1_alph1 
     430            zs2(ji,jj) = ( zs2(ji,jj)*zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 * (1._wp + zkt) ) ) * z1_alph2 
    395431           
    396432         END_2D 
    397          CALL lbc_lnk( 'icedyn_rhg_evp', zp_delt, 'T', 1.0_wp ) 
    398  
     433 
     434         ! Save beta at T-points for further computations 
     435         IF( ln_aEVP ) THEN 
     436            DO_2D( 1, 1, 1, 1 ) 
     437               zbeta(ji,jj) = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     438            END_2D 
     439         ENDIF 
     440          
    399441         DO_2D( 1, 0, 1, 0 ) 
    400442 
    401             ! alpha & beta for aEVP 
     443            ! alpha for aEVP 
    402444            IF( ln_aEVP ) THEN 
    403                zalph2   = MAX( 50._wp, rpi * SQRT( 0.5_wp * zp_delt(ji,jj) * r1_e1e2t(ji,jj) * zdt_m(ji,jj) ) ) 
     445               zalph2   = MAX( zbeta(ji,jj), zbeta(ji+1,jj), zbeta(ji,jj+1), zbeta(ji+1,jj+1) ) 
    404446               z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
    405                zbeta(ji,jj) = zalph2 
     447               ! explicit: 
     448               ! z1_alph2 = 1._wp / zalph2 
     449               ! zalph2 = zalph2 - 1._wp 
    406450            ENDIF 
    407451             
     
    469513               ! 
    470514               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    471                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    472                      &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    473                      &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    474                      &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    475                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     515                  zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 
     516                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     517                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     518                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     519                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   &  
     520                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     521                     &                                    ) / ( zbetav + 1._wp )                                              & 
     522                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    476523                     &           )   * zmsk00y(ji,jj) 
    477524               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    478                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    479                      &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    480                      &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    481                      &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    482                      &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    483                      &            )   * zmsk00y(ji,jj) 
     525                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
     526                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     527                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     528                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     529                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     530                     &            )  * zmsk00y(ji,jj) 
    484531               ENDIF 
    485532            END_2D 
     
    518565               ! 
    519566               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    520                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    521                      &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    522                      &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    523                      &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    524                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     567                  zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 
     568                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     569                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     570                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     571                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   & 
     572                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     573                     &                                    ) / ( zbetau + 1._wp )                                              & 
     574                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    525575                     &           )   * zmsk00x(ji,jj) 
    526576               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    527                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
    528                      &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    529                      &                                     / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    530                      &                + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    531                      &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    532                      &            )   * zmsk00x(ji,jj) 
     577                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
     578                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     579                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     580                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     581                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     582                     &           )   * zmsk00x(ji,jj) 
    533583               ENDIF 
    534584            END_2D 
     
    569619               ! 
    570620               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    571                   u_ice(ji,jj) = ( (          rswitch * ( zmU_t(ji,jj) * ( zbeta(ji,jj) * u_ice(ji,jj) + u_ice_b(ji,jj) )       & ! previous velocity 
    572                      &                                  + zRHS + zTauO * u_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    573                      &                                  / MAX( zepsi, zmU_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    574                      &               + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    575                      &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
     621                  zbetau = MAX( zbeta(ji,jj), zbeta(ji+1,jj) ) 
     622                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * ( zbetau * u_ice(ji,jj) + u_ice_b(ji,jj) )         & ! previous velocity 
     623                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     624                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) * ( zbetau + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     625                     &            + ( 1._wp - rswitch ) * (  u_ice_b(ji,jj)                                                   & 
     626                     &                                     + u_ice  (ji,jj) * MAX( 0._wp, zbetau - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     627                     &                                    ) / ( zbetau + 1._wp )                                              & 
     628                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin  
    576629                     &           )   * zmsk00x(ji,jj) 
    577630               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    578                   u_ice(ji,jj) = ( (           rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
    579                      &                                     + zRHS + zTauO * u_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    580                      &                                     / MAX( zepsi, zmU_t(ji,jj) + 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                      &              ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    583                      &            )   * zmsk00x(ji,jj) 
     631                  u_ice(ji,jj) = ( (          rswitch   * ( zmU_t(ji,jj) * u_ice(ji,jj)                                       & ! previous velocity 
     632                     &                                    + zRHS + zTauO * u_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     633                     &                                    ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     634                     &            + ( 1._wp - rswitch ) *   u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     635                     &             ) * zmsk01x(ji,jj) + u_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01x(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     636                     &           )   * zmsk00x(ji,jj) 
    584637               ENDIF 
    585638            END_2D 
     
    618671               ! 
    619672               IF( ln_aEVP ) THEN !--- ice velocity using aEVP (Kimmritz et al 2016 & 2017) 
    620                   v_ice(ji,jj) = ( (          rswitch * ( zmV_t(ji,jj) * ( zbeta(ji,jj) * v_ice(ji,jj) + v_ice_b(ji,jj) )       & ! previous velocity 
    621                      &                                  + zRHS + zTauO * v_ice(ji,jj) )                                         & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    622                      &                                  / MAX( zepsi, zmV_t(ji,jj) * ( zbeta(ji,jj) + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
    623                      &               + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )           & ! static friction => slow decrease to v=0 
    624                      &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                     & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     673                  zbetav = MAX( zbeta(ji,jj), zbeta(ji,jj+1) ) 
     674                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * ( zbetav * v_ice(ji,jj) + v_ice_b(ji,jj) )         & ! previous velocity 
     675                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     676                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) * ( zbetav + 1._wp ) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 
     677                     &            + ( 1._wp - rswitch ) * (  v_ice_b(ji,jj)                                                   & 
     678                     &                                     + v_ice  (ji,jj) * MAX( 0._wp, zbetav - zdtevp * rn_lf_relax )     & ! static friction => slow decrease to v=0 
     679                     &                                    ) / ( zbetav + 1._wp )                                              &  
     680                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    625681                     &           )   * zmsk00y(ji,jj) 
    626682               ELSE               !--- ice velocity using EVP implicit formulation (cf Madec doc & Bouillon 2009) 
    627                   v_ice(ji,jj) = ( (           rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
    628                      &                                     + zRHS + zTauO * v_ice(ji,jj) )                                      & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
    629                      &                                     / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                         & ! m/dt + tau_io(only ice part) + landfast 
    630                      &                + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax )          & ! static friction => slow decrease to v=0 
    631                      &              ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                    & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
    632                      &            )   * zmsk00y(ji,jj) 
     683                  v_ice(ji,jj) = ( (          rswitch   * ( zmV_t(ji,jj) * v_ice(ji,jj)                                       & ! previous velocity 
     684                     &                                    + zRHS + zTauO * v_ice(ji,jj)                                       & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     685                     &                                    ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB )                      & ! m/dt + tau_io(only ice part) + landfast 
     686                     &            + ( 1._wp - rswitch ) *   v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lf_relax )         & ! static friction => slow decrease to v=0 
     687                     &             ) * zmsk01y(ji,jj) + v_oce(ji,jj) * 0.01_wp * ( 1._wp - zmsk01y(ji,jj) )                   & ! v_ice = v_oce/100 if mass < zmmin & conc < zamin 
     688                     &           )   * zmsk00y(ji,jj) 
    633689               ENDIF 
    634690            END_2D 
     
    643699         ENDIF 
    644700 
    645 !!$         IF(sn_cfctl%l_prtctl) THEN   ! Convergence test 
    646 !!$            DO jj = 2 , jpjm1 
    647 !!$               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    648 !!$            END DO 
    649 !!$            zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 
    650 !!$            CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
    651 !!$         ENDIF 
     701         ! convergence test 
     702         IF( nn_rhg_chkcvg == 2 )   CALL rhg_cvg( kt, jter, nn_nevp, u_ice, v_ice, zu_ice, zv_ice ) 
    652703         ! 
    653704         !                                                ! ==================== ! 
    654705      END DO                                              !  end loop over jter  ! 
    655706      !                                                   ! ==================== ! 
     707      IF( ln_aEVP )   CALL iom_put( 'beta_evp' , zbeta ) 
    656708      ! 
    657709      !------------------------------------------------------------------------------! 
     
    667719      END_2D 
    668720       
    669       DO_2D( 0, 0, 0, 0 ) 
     721      DO_2D( 0, 0, 0, 0 )   ! no vector loop 
    670722          
    671723         ! tension**2 at T points 
     
    689741          
    690742         ! delta at T points 
    691          zdelta         = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
    692          rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
    693          pdelta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     743         zdelta(ji,jj)   = SQRT( pdivu_i(ji,jj) * pdivu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 )   
     744         rswitch         = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta(ji,jj) ) ) ! 0 if delta=0 
     745         pdelta_i(ji,jj) = zdelta(ji,jj) + rn_creepl * rswitch 
    694746 
    695747      END_2D 
     
    706758      ! 5) diagnostics 
    707759      !------------------------------------------------------------------------------! 
    708       DO_2D( 1, 1, 1, 1 ) 
    709          zmsk00(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) ! 1 if ice, 0 if no ice 
    710       END_2D 
    711  
    712760      ! --- ice-ocean, ice-atm. & ice-oceanbottom(landfast) stresses --- ! 
    713761      IF(  iom_use('utau_oi') .OR. iom_use('vtau_oi') .OR. iom_use('utau_ai') .OR. iom_use('vtau_ai') .OR. & 
     
    764812         DEALLOCATE( zsig1 , zsig2 , zsig3 ) 
    765813      ENDIF 
    766        
     814 
    767815      ! --- SIMIP --- ! 
    768816      IF(  iom_use('dssh_dx') .OR. iom_use('dssh_dy') .OR. & 
     
    818866      ENDIF 
    819867      ! 
     868      ! --- convergence tests --- ! 
     869      IF( nn_rhg_chkcvg == 1 .OR. nn_rhg_chkcvg == 2 ) THEN 
     870         IF( iom_use('uice_cvg') ) THEN 
     871            IF( ln_aEVP ) THEN   ! output: beta * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
     872               CALL iom_put( 'uice_cvg', MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * zbeta(:,:) * umask(:,:,1) , & 
     873                  &                           ABS( v_ice(:,:) - zv_ice(:,:) ) * zbeta(:,:) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     874            ELSE                 ! output: nn_nevp * ( u(t=nn_nevp) - u(t=nn_nevp-1) ) 
     875               CALL iom_put( 'uice_cvg', REAL( nn_nevp ) * MAX( ABS( u_ice(:,:) - zu_ice(:,:) ) * umask(:,:,1) , & 
     876                  &                                             ABS( v_ice(:,:) - zv_ice(:,:) ) * vmask(:,:,1) ) * zmsk15(:,:) ) 
     877            ENDIF 
     878         ENDIF 
     879      ENDIF       
     880      ! 
     881      DEALLOCATE( zmsk00, zmsk15 ) 
     882      ! 
    820883   END SUBROUTINE ice_dyn_rhg_evp 
     884 
     885 
     886   SUBROUTINE rhg_cvg( kt, kiter, kitermax, pu, pv, pub, pvb ) 
     887      !!---------------------------------------------------------------------- 
     888      !!                    ***  ROUTINE rhg_cvg  *** 
     889      !!                      
     890      !! ** Purpose :   check convergence of oce rheology 
     891      !! 
     892      !! ** Method  :   create a file ice_cvg.nc containing the convergence of ice velocity 
     893      !!                during the sub timestepping of rheology so as: 
     894      !!                  uice_cvg = MAX( u(t+1) - u(t) , v(t+1) - v(t) ) 
     895      !!                This routine is called every sub-iteration, so it is cpu expensive 
     896      !! 
     897      !! ** Note    :   for the first sub-iteration, uice_cvg is set to 0 (too large otherwise)    
     898      !!---------------------------------------------------------------------- 
     899      INTEGER ,                 INTENT(in) ::   kt, kiter, kitermax       ! ocean time-step index 
     900      REAL(wp), DIMENSION(:,:), INTENT(in) ::   pu, pv, pub, pvb          ! now and before velocities 
     901      !! 
     902      INTEGER           ::   it, idtime, istatus 
     903      INTEGER           ::   ji, jj          ! dummy loop indices 
     904      REAL(wp)          ::   zresm           ! local real  
     905      CHARACTER(len=20) ::   clname 
     906      REAL(wp), DIMENSION(jpi,jpj) ::   zres           ! check convergence 
     907      !!---------------------------------------------------------------------- 
     908 
     909      ! create file 
     910      IF( kt == nit000 .AND. kiter == 1 ) THEN 
     911         ! 
     912         IF( lwp ) THEN 
     913            WRITE(numout,*) 
     914            WRITE(numout,*) 'rhg_cvg : ice rheology convergence control' 
     915            WRITE(numout,*) '~~~~~~~' 
     916         ENDIF 
     917         ! 
     918         IF( lwm ) THEN 
     919            clname = 'ice_cvg.nc' 
     920            IF( .NOT. Agrif_Root() )   clname = TRIM(Agrif_CFixed())//"_"//TRIM(clname) 
     921            istatus = NF90_CREATE( TRIM(clname), NF90_CLOBBER, ncvgid ) 
     922            istatus = NF90_DEF_DIM( ncvgid, 'time'  , NF90_UNLIMITED, idtime ) 
     923            istatus = NF90_DEF_VAR( ncvgid, 'uice_cvg', NF90_DOUBLE   , (/ idtime /), nvarid ) 
     924            istatus = NF90_ENDDEF(ncvgid) 
     925         ENDIF 
     926         ! 
     927      ENDIF 
     928 
     929      ! time 
     930      it = ( kt - 1 ) * kitermax + kiter 
     931       
     932      ! convergence 
     933      IF( kiter == 1 ) THEN ! remove the first iteration for calculations of convergence (always very large) 
     934         zresm = 0._wp 
     935      ELSE 
     936         DO_2D( 1, 1, 1, 1 ) 
     937            zres(ji,jj) = MAX( ABS( pu(ji,jj) - pub(ji,jj) ) * umask(ji,jj,1), & 
     938               &               ABS( pv(ji,jj) - pvb(ji,jj) ) * vmask(ji,jj,1) ) * zmsk15(ji,jj) 
     939         END_2D 
     940         zresm = MAXVAL( zres ) 
     941         CALL mpp_max( 'icedyn_rhg_evp', zresm )   ! max over the global domain 
     942      ENDIF 
     943 
     944      IF( lwm ) THEN 
     945         ! write variables 
     946         istatus = NF90_PUT_VAR( ncvgid, nvarid, (/zresm/), (/it/), (/1/) ) 
     947         ! close file 
     948         IF( kt == nitend )   istatus = NF90_CLOSE(ncvgid) 
     949      ENDIF 
     950       
     951   END SUBROUTINE rhg_cvg 
    821952 
    822953 
     
    8761007   END SUBROUTINE rhg_evp_rst 
    8771008 
     1009    
    8781010#else 
    8791011   !!---------------------------------------------------------------------- 
  • NEMO/branches/2020/dev_r13383_HPC-02_Daley_Tiling/src/ICE/iceistate.F90

    r13295 r13553  
    4747   !                             !! ** namelist (namini) ** 
    4848   LOGICAL, PUBLIC  ::   ln_iceini        !: Ice initialization or not 
    49    LOGICAL, PUBLIC  ::   ln_iceini_file   !: Ice initialization from 2D netcdf file 
     49   INTEGER, PUBLIC  ::   nn_iceini_file   !: Ice initialization: 
     50                                  !        0 = Initialise sea ice based on SSTs 
     51                                  !        1 = Initialise sea ice from single category netcdf file 
     52                                  !        2 = Initialise sea ice from multi category restart file 
    5053   REAL(wp) ::   rn_thres_sst 
    5154   REAL(wp) ::   rn_hti_ini_n, rn_hts_ini_n, rn_ati_ini_n, rn_smi_ini_n, rn_tmi_ini_n, rn_tsu_ini_n, rn_tms_ini_n 
    5255   REAL(wp) ::   rn_hti_ini_s, rn_hts_ini_s, rn_ati_ini_s, rn_smi_ini_s, rn_tmi_ini_s, rn_tsu_ini_s, rn_tms_ini_s 
    53    REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n 
    54    REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s 
     56   REAL(wp) ::   rn_apd_ini_n, rn_hpd_ini_n, rn_hld_ini_n 
     57   REAL(wp) ::   rn_apd_ini_s, rn_hpd_ini_s, rn_hld_ini_s 
    5558   ! 
    56    !                              ! if ln_iceini_file = T 
    57    INTEGER , PARAMETER ::   jpfldi = 9           ! maximum number of files to read 
     59   !                              ! if nn_iceini_file = 1 
     60   INTEGER , PARAMETER ::   jpfldi = 10          ! maximum number of files to read 
    5861   INTEGER , PARAMETER ::   jp_hti = 1           ! index of ice thickness    (m) 
    5962   INTEGER , PARAMETER ::   jp_hts = 2           ! index of snw thickness    (m) 
     
    6568   INTEGER , PARAMETER ::   jp_apd = 8           ! index of pnd fraction     (-) 
    6669   INTEGER , PARAMETER ::   jp_hpd = 9           ! index of pnd depth        (m) 
     70   INTEGER , PARAMETER ::   jp_hld = 10          ! index of pnd lid depth    (m) 
    6771   TYPE(FLD), ALLOCATABLE, DIMENSION(:) ::   si  ! structure of input fields (file informations, fields read) 
    6872 
     
    8993      !! ** Steps   :   1) Set initial surface and basal temperatures 
    9094      !!                2) Recompute or read sea ice state variables 
    91       !!                3) Fill in the ice thickness distribution using gaussian 
    92       !!                4) Fill in space-dependent arrays for state variables 
    93       !!                5) snow-ice mass computation 
    94       !!                6) store before fields 
     95      !!                3) Fill in space-dependent arrays for state variables 
     96      !!                4) snow-ice mass computation 
    9597      !! 
    9698      !! ** Notes   : o_i, t_su, t_s, t_i, sz_i must be filled everywhere, even 
     
    107109      REAL(wp), DIMENSION(jpi,jpj)     ::   zht_i_ini, zat_i_ini, ztm_s_ini            !data from namelist or nc file 
    108110      REAL(wp), DIMENSION(jpi,jpj)     ::   zt_su_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 
    109       REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini                       !data from namelist or nc file 
    110       REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !locak arrays 
    111       !! 
    112       REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d 
     111      REAL(wp), DIMENSION(jpi,jpj)     ::   zapnd_ini, zhpnd_ini, zhlid_ini            !data from namelist or nc file 
     112      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   zti_3d , zts_3d                            !temporary arrays 
     113      !! 
     114      REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   zhi_2d, zhs_2d, zai_2d, zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d 
    113115      !-------------------------------------------------------------------- 
    114116 
     
    164166      a_ip     (:,:,:) = 0._wp 
    165167      v_ip     (:,:,:) = 0._wp 
    166       a_ip_frac(:,:,:) = 0._wp 
     168      v_il     (:,:,:) = 0._wp 
     169      a_ip_eff (:,:,:) = 0._wp 
    167170      h_ip     (:,:,:) = 0._wp 
     171      h_il     (:,:,:) = 0._wp 
    168172      ! 
    169173      ! ice velocities 
     
    174178      ! 2) overwrite some of the fields with namelist parameters or netcdf file 
    175179      !------------------------------------------------------------------------ 
    176  
    177  
    178180      IF( ln_iceini ) THEN 
    179          !                             !---------------! 
    180           
     181         ! 
    181182         IF( Agrif_Root() ) THEN 
    182  
    183             IF( ln_iceini_file )THEN      ! Read a file   ! 
     183            !                             !---------------! 
     184            IF( nn_iceini_file == 1 )THEN ! Read a file   ! 
    184185               !                          !---------------! 
    185186               WHERE( ff_t(:,:) >= 0._wp )   ;   zswitch(:,:) = 1._wp 
     
    195196 
    196197               ! -- optional fields -- ! 
    197                !    if fields do not exist then set them to the values present in the namelist (except for snow and surface temperature) 
     198               !    if fields do not exist then set them to the values present in the namelist (except for temperatures) 
    198199               ! 
    199200               ! ice salinity 
     
    207208                  si(jp_tsu)%fnow(:,:,1) = ( rn_tsu_ini_n * zswitch + rn_tsu_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    208209                  si(jp_tms)%fnow(:,:,1) = ( rn_tms_ini_n * zswitch + rn_tms_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    209                ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
    210                   si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 
    211                ELSEIF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
    212                   si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 
    213                ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) THEN ! if T_s is read and not T_su, set T_su = T_s 
    214                   si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 
    215                ELSEIF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_su, set T_su = T_i 
    216                   si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    217                ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) THEN ! if T_su is read and not T_s, set T_s = T_su 
    218                   si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 
    219                ELSEIF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) THEN ! if T_i is read and not T_s, set T_s = T_i 
    220                   si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    221210               ENDIF 
     211               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_i, set T_i = (T_s + T_freeze)/2 
     212                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tms)%fnow(:,:,1) + 271.15 ) 
     213               IF( TRIM(si(jp_tmi)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_i, set T_i = (T_su + T_freeze)/2 
     214                  &     si(jp_tmi)%fnow(:,:,1) = 0.5_wp * ( si(jp_tsu)%fnow(:,:,1) + 271.15 ) 
     215               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tms)%clrootname) /= 'NOT USED' ) & ! if T_s is read and not T_su, set T_su = T_s 
     216                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tms)%fnow(:,:,1) 
     217               IF( TRIM(si(jp_tsu)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_su, set T_su = T_i 
     218                  &     si(jp_tsu)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
     219               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tsu)%clrootname) /= 'NOT USED' ) & ! if T_su is read and not T_s, set T_s = T_su 
     220                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tsu)%fnow(:,:,1) 
     221               IF( TRIM(si(jp_tms)%clrootname) == 'NOT USED' .AND. TRIM(si(jp_tmi)%clrootname) /= 'NOT USED' ) & ! if T_i is read and not T_s, set T_s = T_i 
     222                  &     si(jp_tms)%fnow(:,:,1) = si(jp_tmi)%fnow(:,:,1) 
    222223               ! 
    223224               ! pond concentration 
     
    229230               IF( TRIM(si(jp_hpd)%clrootname) == 'NOT USED' ) & 
    230231                  &     si(jp_hpd)%fnow(:,:,1) = ( rn_hpd_ini_n * zswitch + rn_hpd_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
     232               ! 
     233               ! pond lid depth 
     234               IF( TRIM(si(jp_hld)%clrootname) == 'NOT USED' ) & 
     235                  &     si(jp_hld)%fnow(:,:,1) = ( rn_hld_ini_n * zswitch + rn_hld_ini_s * (1._wp - zswitch) ) * tmask(:,:,1) 
    231236               ! 
    232237               zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) * tmask(:,:,1) 
     
    236241               zapnd_ini(:,:) = si(jp_apd)%fnow(:,:,1) * tmask(:,:,1) 
    237242               zhpnd_ini(:,:) = si(jp_hpd)%fnow(:,:,1) * tmask(:,:,1) 
     243               zhlid_ini(:,:) = si(jp_hld)%fnow(:,:,1) * tmask(:,:,1) 
    238244               ! 
    239245               ! change the switch for the following 
     
    261267                  zapnd_ini(:,:) = rn_apd_ini_n * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc.  
    262268                  zhpnd_ini(:,:) = rn_hpd_ini_n * zswitch(:,:) 
     269                  zhlid_ini(:,:) = rn_hld_ini_n * zswitch(:,:) 
    263270               ELSEWHERE 
    264271                  zht_i_ini(:,:) = rn_hti_ini_s * zswitch(:,:) 
     
    271278                  zapnd_ini(:,:) = rn_apd_ini_s * zswitch(:,:) * zat_i_ini(:,:) ! rn_apd = pond fraction => rn_apd * a_i = pond conc. 
    272279                  zhpnd_ini(:,:) = rn_hpd_ini_s * zswitch(:,:) 
     280                  zhlid_ini(:,:) = rn_hld_ini_s * zswitch(:,:) 
    273281               END WHERE 
    274282               ! 
     
    281289               zapnd_ini(:,:) = 0._wp 
    282290               zhpnd_ini(:,:) = 0._wp 
     291               zhlid_ini(:,:) = 0._wp 
    283292            ENDIF 
    284293             
    285             !-------------! 
    286             ! fill fields ! 
    287             !-------------! 
     294            IF ( .NOT.ln_pnd_lids ) THEN 
     295               zhlid_ini(:,:) = 0._wp 
     296            ENDIF 
     297             
     298            !----------------! 
     299            ! 3) fill fields ! 
     300            !----------------! 
    288301            ! select ice covered grid points 
    289302            npti = 0 ; nptidx(:) = 0 
     
    305318            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d(1:npti)  , zapnd_ini ) 
    306319            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d(1:npti)  , zhpnd_ini ) 
    307  
     320            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d(1:npti)  , zhlid_ini ) 
     321             
    308322            ! allocate temporary arrays 
    309             ALLOCATE( zhi_2d(npti,jpl), zhs_2d(npti,jpl), zai_2d (npti,jpl), & 
    310                &      zti_2d(npti,jpl), zts_2d(npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), zaip_2d(npti,jpl), zhip_2d(npti,jpl) ) 
    311              
     323            ALLOCATE( zhi_2d (npti,jpl), zhs_2d (npti,jpl), zai_2d (npti,jpl), & 
     324               &      zti_2d (npti,jpl), zts_2d (npti,jpl), ztsu_2d(npti,jpl), zsi_2d(npti,jpl), & 
     325               &      zaip_2d(npti,jpl), zhip_2d(npti,jpl), zhil_2d(npti,jpl) ) 
     326 
    312327            ! distribute 1-cat into jpl-cat: (jpi*jpj) -> (jpi*jpj,jpl) 
    313             CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                                                   & 
    314                &              zhi_2d          , zhs_2d          , zai_2d         ,                                                   & 
    315                &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti), s_i_1d(1:npti), a_ip_1d(1:npti), h_ip_1d(1:npti), & 
    316                &              zti_2d          , zts_2d          , ztsu_2d        , zsi_2d        , zaip_2d        , zhip_2d ) 
     328            CALL ice_var_itd( h_i_1d(1:npti)  , h_s_1d(1:npti)  , at_i_1d(1:npti),                  & 
     329               &              zhi_2d          , zhs_2d          , zai_2d         ,                  & 
     330               &              t_i_1d(1:npti,1), t_s_1d(1:npti,1), t_su_1d(1:npti),                  & 
     331               &              s_i_1d(1:npti)  , a_ip_1d(1:npti) , h_ip_1d(1:npti), h_il_1d(1:npti), & 
     332               &              zti_2d          , zts_2d          , ztsu_2d        ,                  & 
     333               &              zsi_2d          , zaip_2d         , zhip_2d        , zhil_2d ) 
    317334 
    318335            ! move to 3D arrays: (jpi*jpj,jpl) -> (jpi,jpj,jpl) 
     
    330347            CALL tab_2d_3d( npti, nptidx(1:npti), zaip_2d  , a_ip   ) 
    331348            CALL tab_2d_3d( npti, nptidx(1:npti), zhip_2d  , h_ip   ) 
     349            CALL tab_2d_3d( npti, nptidx(1:npti), zhil_2d  , h_il   ) 
    332350 
    333351            ! deallocate temporary arrays 
    334352            DEALLOCATE( zhi_2d, zhs_2d, zai_2d , & 
    335                &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d ) 
     353               &        zti_2d, zts_2d, ztsu_2d, zsi_2d, zaip_2d, zhip_2d, zhil_2d ) 
    336354 
    337355            ! calculate extensive and intensive variables 
     
    363381               END_3D 
    364382            END DO 
    365  
    366             ! Melt ponds 
    367             WHERE( a_i > epsi10 ) 
    368                a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
    369             ELSEWHERE 
    370                a_ip_frac(:,:,:) = 0._wp 
    371             END WHERE 
    372             v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
    373               
    374             ! specific temperatures for coupled runs 
    375             tn_ice(:,:,:) = t_su(:,:,:) 
    376             t1_ice(:,:,:) = t_i (:,:,1,:) 
    377             ! 
    378           
     383             
    379384#if  defined key_agrif 
    380385         ELSE 
     
    391396            Agrif_UseSpecialValue = .FALSE. 
    392397        ! lbc ????  
    393    ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, t_su, e_s, e_i 
     398   ! Here we know : a_i, v_i, v_s, sv_i, oa_i, a_ip, v_ip, v_il, t_su, e_s, e_i 
    394399            CALL ice_var_glo2eqv 
    395400            CALL ice_var_zapsmall 
    396401            CALL ice_var_agg(2) 
    397  
    398             ! Melt ponds 
    399             WHERE( a_i > epsi10 ) 
    400                a_ip_frac(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
    401             ELSEWHERE 
    402                a_ip_frac(:,:,:) = 0._wp 
    403             END WHERE 
    404             WHERE( a_ip > 0._wp )       ! ???????     
    405                h_ip(:,:,:) = v_ip(:,:,:) / a_ip(:,:,:) 
    406             ELSEWHERE 
    407                h_ip(:,:,:) = 0._wp 
    408             END WHERE    
    409  
    410             tn_ice(:,:,:) = t_su(:,:,:) 
    411             t1_ice(:,:,:) = t_i (:,:,1,:) 
    412402#endif 
    413           ENDIF ! Agrif_Root 
     403         ENDIF ! Agrif_Root 
     404         ! 
     405         ! Melt ponds 
     406         WHERE( a_i > epsi10 )   ;   a_ip_eff(:,:,:) = a_ip(:,:,:) / a_i(:,:,:) 
     407         ELSEWHERE               ;   a_ip_eff(:,:,:) = 0._wp 
     408         END WHERE 
     409         v_ip(:,:,:) = h_ip(:,:,:) * a_ip(:,:,:) 
     410         v_il(:,:,:) = h_il(:,:,:) * a_ip(:,:,:) 
     411          
     412         ! specific temperatures for coupled runs 
     413         tn_ice(:,:,:) = t_su(:,:,:) 
     414         t1_ice(:,:,:) = t_i (:,:,1,:) 
     415         ! 
     416         ! ice concentration should not exceed amax 
     417         at_i(:,:) = SUM( a_i, dim=3 ) 
     418         DO jl = 1, jpl 
     419            WHERE( at_i(:,:) > rn_amax_2d(:,:) )   a_i(:,:,jl) = a_i(:,:,jl) * rn_amax_2d(:,:) / at_i(:,:) 
     420         END DO 
     421         at_i(:,:) = SUM( a_i, dim=3 ) 
     422         ! 
    414423      ENDIF ! ln_iceini 
    415424      ! 
    416       at_i(:,:) = SUM( a_i, dim=3 ) 
    417       ! 
    418425      !---------------------------------------------- 
    419       ! 3) Snow-ice mass (case ice is fully embedded) 
     426      ! 4) Snow-ice mass (case ice is fully embedded) 
    420427      !---------------------------------------------- 
    421428      snwice_mass  (:,:) = tmask(:,:,1) * SUM( rhos * v_s(:,:,:) + rhoi * v_i(:,:,:), dim=3  )   ! snow+ice mass 
     
    469476!          ENDIF 
    470477      ENDIF 
    471        
    472       !------------------------------------ 
    473       ! 4) store fields at before time-step 
    474       !------------------------------------ 
    475       ! it is only necessary for the 1st interpolation by Agrif 
    476       a_i_b  (:,:,:)   = a_i  (:,:,:) 
    477       e_i_b  (:,:,:,:) = e_i  (:,:,:,:) 
    478       v_i_b  (:,:,:)   = v_i  (:,:,:) 
    479       v_s_b  (:,:,:)   = v_s  (:,:,:) 
    480       e_s_b  (:,:,:,:) = e_s  (:,:,:,:) 
    481       sv_i_b (:,:,:)   = sv_i (:,:,:) 
    482       oa_i_b (:,:,:)   = oa_i (:,:,:) 
    483       u_ice_b(:,:)     = u_ice(:,:) 
    484       v_ice_b(:,:)     = v_ice(:,:) 
    485       ! total concentration is needed for Lupkes parameterizations 
    486       at_i_b (:,:)     = at_i (:,:)  
    487  
    488 !!clem: output of initial state should be written here but it is impossible because 
    489 !!      the ocean and ice are in the same file 
    490 !!      CALL dia_wri_state( Kmm, 'output.init' ) 
     478 
     479      !!clem: output of initial state should be written here but it is impossible because 
     480      !!      the ocean and ice are in the same file 
     481      !!      CALL dia_wri_state( 'output.init' ) 
    491482      ! 
    492483   END SUBROUTINE ice_istate 
     
    505496      !! 
    506497      !!----------------------------------------------------------------------------- 
    507       INTEGER ::   ios, ifpr, ierror   ! Local integers 
    508  
     498      INTEGER ::   ios   ! Local integer output status for namelist read 
     499      INTEGER ::   ifpr, ierror 
    509500      ! 
    510501      CHARACTER(len=256) ::  cn_dir          ! Root directory for location of ice files 
    511       TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd 
     502      TYPE(FLD_N)                    ::   sn_hti, sn_hts, sn_ati, sn_smi, sn_tmi, sn_tsu, sn_tms, sn_apd, sn_hpd, sn_hld 
    512503      TYPE(FLD_N), DIMENSION(jpfldi) ::   slf_i                 ! array of namelist informations on the fields to read 
    513504      ! 
    514       NAMELIST/namini/ ln_iceini, ln_iceini_file, rn_thres_sst, & 
     505      NAMELIST/namini/ ln_iceini, nn_iceini_file, rn_thres_sst, & 
    515506         &             rn_hti_ini_n, rn_hti_ini_s, rn_hts_ini_n, rn_hts_ini_s, & 
    516507         &             rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, rn_smi_ini_s, & 
    517508         &             rn_tmi_ini_n, rn_tmi_ini_s, rn_tsu_ini_n, rn_tsu_ini_s, rn_tms_ini_n, rn_tms_ini_s, & 
    518          &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, & 
    519          &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, cn_dir 
     509         &             rn_apd_ini_n, rn_apd_ini_s, rn_hpd_ini_n, rn_hpd_ini_s, rn_hld_ini_n, rn_hld_ini_s, & 
     510         &             sn_hti, sn_hts, sn_ati, sn_tsu, sn_tmi, sn_smi, sn_tms, sn_apd, sn_hpd, sn_hld, cn_dir 
    520511      !!----------------------------------------------------------------------------- 
    521512      ! 
     
    529520      slf_i(jp_ati) = sn_ati  ;  slf_i(jp_smi) = sn_smi 
    530521      slf_i(jp_tmi) = sn_tmi  ;  slf_i(jp_tsu) = sn_tsu   ;   slf_i(jp_tms) = sn_tms 
    531       slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd 
     522      slf_i(jp_apd) = sn_apd  ;  slf_i(jp_hpd) = sn_hpd   ;   slf_i(jp_hld) = sn_hld 
    532523      ! 
    533524      IF(lwp) THEN                          ! control print 
     
    537528         WRITE(numout,*) '   Namelist namini:' 
    538529         WRITE(numout,*) '      ice initialization (T) or not (F)                ln_iceini      = ', ln_iceini 
    539          WRITE(numout,*) '      ice initialization from a netcdf file            ln_iceini_file = ', ln_iceini_file 
     530         WRITE(numout,*) '      ice initialization from a netcdf file            nn_iceini_file = ', nn_iceini_file 
    540531         WRITE(numout,*) '      max ocean temp. above Tfreeze with initial ice   rn_thres_sst   = ', rn_thres_sst 
    541          IF( ln_iceini .AND. .NOT.ln_iceini_file ) THEN 
     532         IF( ln_iceini .AND. nn_iceini_file == 0 ) THEN 
    542533            WRITE(numout,*) '      initial snw thickness in the north-south         rn_hts_ini     = ', rn_hts_ini_n,rn_hts_ini_s  
    543534            WRITE(numout,*) '      initial ice thickness in the north-south         rn_hti_ini     = ', rn_hti_ini_n,rn_hti_ini_s 
     
    549540            WRITE(numout,*) '      initial pnd fraction  in the north-south         rn_apd_ini     = ', rn_apd_ini_n,rn_apd_ini_s 
    550541            WRITE(numout,*) '      initial pnd depth     in the north-south         rn_hpd_ini     = ', rn_hpd_ini_n,rn_hpd_ini_s 
     542            WRITE(numout,*) '      initial pnd lid depth in the north-south         rn_hld_ini     = ', rn_hld_ini_n,rn_hld_ini_s 
    551543         ENDIF 
    552544      ENDIF 
    553545      ! 
    554       IF( ln_iceini_file ) THEN                      ! Ice initialization using input file 
     546      IF( nn_iceini_file == 1 ) THEN                      ! Ice initialization using input file 
    555547         ! 
    556548         ! set si structure 
     
    573565         rn_apd_ini_n = 0. ; rn_apd_ini_s = 0. 
    574566         rn_hpd_ini_n = 0. ; rn_hpd_ini_s = 0. 
    575          CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 when no ponds' ) 
     567         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0. 
     568         CALL ctl_warn( 'rn_apd_ini & rn_hpd_ini = 0 & rn_hld_ini = 0 when no ponds' ) 
     569      ENDIF 
     570      ! 
     571      IF( .NOT.ln_pnd_lids ) THEN 
     572         rn_hld_ini_n = 0. ; rn_hld_ini_s = 0. 
    576573      ENDIF 
    577574      ! 
    </