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

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

Changeset 15658


Ignore:
Timestamp:
2022-01-19T19:34:39+01:00 (2 years ago)
Author:
jpalmier
Message:

Merge with first branch : NEMO_4.0.4_GO8_package

Location:
NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch
Files:
40 edited
2 copied

Legend:

Unmodified
Added
Removed

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-ice.xml

    r13648 r15658  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
    30        <field field_ref="sst_m"            name="sst_m"  /> 
    31        <field field_ref="sss_m"            name="sss_m"  /> 
    32         
     30       <field field_ref="sst_m_pot"        name="sst_m_pot"  /> 
     31       <field field_ref="sss_m_abs"        name="sss_m_abs"  /> 
     32 
    3333       <!-- heat --> 
    3434       <field field_ref="icetemp"          name="sitemp" /> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/ORCA2_ICE_PISCES/EXPREF/file_def_nemo-oce.xml

    r12276 r15658  
    1414        <file id="file11" name_suffix="_grid_T" description="ocean T grid variables" > 
    1515          <field field_ref="e3t"      /> 
    16           <field field_ref="toce"         name="thetao"                                                                      operation="instant" freq_op="5d" > @toce_e3t / @e3t </field> 
    17           <field field_ref="soce"         name="so"                                                                          operation="instant" freq_op="5d" > @soce_e3t / @e3t </field> 
    18           <field field_ref="sst"          name="tos"   /> 
    19           <field field_ref="sss"          name="sos"   /> 
     16          <field field_ref="toce_con"         name="thetao_con"                                                                      operation="instant" freq_op="5d" > @toce_con_e3t / @e3t </field> 
     17          <field field_ref="soce_abs"         name="so_abs"                                                                          operation="instant" freq_op="5d" > @soce_abs_e3t / @e3t </field> 
     18          <field field_ref="sst_con"          name="tos_con"   /> 
     19          <field field_ref="sss_abs"          name="sos_abs"   /> 
    2020          <field field_ref="ssh"          name="zos"   /> 
    21           <field field_ref="sst"          name="tosstd"       long_name="sea surface temperature standard deviation"         operation="instant" freq_op="5d" > sqrt( @sst2 - @sst * @sst ) </field> 
     21          <field field_ref="sst_con"          name="tosstd_con"       long_name="sea surface temperature standard deviation"         operation="instant" freq_op="5d" > sqrt( @sst2_con - @sst_con * @sst_con ) </field> 
    2222          <field field_ref="ssh"          name="zosstd"       long_name="sea surface height above geoid standard deviation"  operation="instant" freq_op="5d" > sqrt( @ssh2 - @ssh * @ssh ) </field> 
    23           <field field_ref="sst"          name="sstdcy"       long_name="amplitude of sst diurnal cycle"                     operation="average" freq_op="1d" > @sstmax - @sstmin </field> 
     23          <field field_ref="sst_con"          name="sstdcy_con"       long_name="amplitude of sst diurnal cycle"                     operation="average" freq_op="1d" > @sstmax_con - @sstmin_con </field> 
    2424          <field field_ref="mldkz5"   /> 
    2525          <field field_ref="mldr10_1" /> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/ORCA2_SAS_ICE/EXPREF/file_def_nemo-ice.xml

    r9896 r15658  
    2828       <field field_ref="iceapnd"          name="siapnd" /> 
    2929       <field field_ref="icevpnd"          name="sivpnd" /> 
     30       <field field_ref="icevlid"          name="sivlid" /> 
    3031       <field field_ref="sst_m"            name="sst_m"  /> 
    3132       <field field_ref="sss_m"            name="sss_m"  /> 
    … …  
    3637       <field field_ref="icetbot"          name="sitbot" /> 
    3738       <field field_ref="icetsni"          name="sitsni" /> 
     39 
     40       <!-- ponds --> 
     41       <field field_ref="dvpn_mlt"         name="dvpn_mlt" /> 
     42       <field field_ref="dvpn_lid"         name="dvpn_lid" /> 
     43       <field field_ref="dvpn_rnf"         name="dvpn_rnf" /> 
     44       <field field_ref="dvpn_drn"         name="dvpn_drn" /> 
    3845        
    3946       <!-- momentum --> 
    … …  
    8592       <field field_ref="icesalt_cat"      name="sisalcat"/> 
    8693       <field field_ref="icetemp_cat"      name="sitemcat"/> 
     94       <field field_ref="iceapnd_cat"      name="siapncat"/> 
     95       <field field_ref="icevpnd_cat"      name="sivpncat"/> 
    8796        
    8897     </file> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/SHARED/domain_def_nemo.xml

    r12276 r15658  
    181181     <domain id="EqW" domain_ref="grid_W" > <zoom_domain id="EqW"/> </domain> 
    182182 
    183               <!--   zonal mean grid   --> 
    184      <domain_group id="gznl"> 
    185         <domain id="gznl" long_name="gznl"/> 
    186         <domain id="ptr" domain_ref="gznl" >  
    187             <zoom_domain id="ptr" ibegin="0000" jbegin="0" ni="1" nj="0000" />  
    188         </domain> 
    189       </domain_group> 
     183      
     184     <!--   zonal mean grid   --> 
     185     <domain id="gznl" long_name="gznl"/> 
     186     <domain id="ptr" domain_ref="gznl" > 
     187    <zoom_domain id="ptr" ibegin="0000" jbegin="0" ni="1" nj="0000" /> 
     188     </domain>   
     189     <domain id="znl_T" domain_ref="gznl" > <zoom_domain id="znl_T"/> </domain> 
     190     <domain id="znl_W" domain_ref="gznl" > <zoom_domain id="znl_W"/> </domain> 
    190191 
    191192      

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/SHARED/field_def_nemo-ice.xml

    r13648 r15658  
    4444          <field id="icefrb"       long_name="Sea-ice freeboard"                                       standard_name="sea_ice_freeboard"                         unit="m"    /> 
    4545          <field id="icealb"       long_name="Sea-ice or snow albedo"                                  standard_name="sea_ice_albedo"                            unit=""    detect_missing_value="true" /> 
    46       
    47      <!-- melt ponds --> 
    48      <field id="iceapnd"      long_name="melt pond concentration"                                 standard_name="sea_ice_meltpond_concentration"            unit=""  />  
     46 
     47          <!-- melt ponds --> 
     48          <field id="iceapnd"      long_name="melt pond concentration"                                 standard_name="sea_ice_meltpond_concentration"            unit=""  />  
    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="iceepnd"      long_name="melt pond effective concentration"                       standard_name="sea_ice_meltpond_effective_concentration"            unit=""  />   
    5152          <field id="icehlid"      long_name="melt pond lid depth"                                     standard_name="sea_ice_meltpondlid_depth"                 unit="m" />  
    5253          <field id="icevlid"      long_name="melt pond lid volume"                                    standard_name="sea_ice_meltpondlid_volume"                unit="m" />  
     54      
     55     <field id="dvpn_mlt"     long_name="pond volume tendency due to surface melt"                standard_name="sea_ice_pondvolume_tendency_melt"          unit="cm/d" />  
     56     <field id="dvpn_lid"     long_name="pond volume tendency due to exchanges with lid"          standard_name="sea_ice_pondvolume_tendency_lids"          unit="cm/d" />  
     57     <field id="dvpn_rnf"     long_name="pond volume tendency due to runoff"                      standard_name="sea_ice_pondvolume_tendency_runoff"        unit="cm/d" />  
     58     <field id="dvpn_drn"     long_name="pond volume tendency due to drainage"                    standard_name="sea_ice_pondvolume_tendency_drainage"      unit="cm/d" />  
    5359      
    5460     <!-- heat --> 
    … …  
    169175 
    170176     <!-- sbcssm variables --> 
    171           <field id="sst_m"    unit="degC" /> 
    172           <field id="sss_m"    unit="psu"  /> 
     177          <field id="sst_m_pot"    unit="degC" /> 
     178     <!-- EOS-80 --> 
     179          <field id="sss_m_pra"    unit="psu"  /> 
     180          <!-- TEOS-10 --> 
     181          <field id="sss_m_abs"    unit="psu"  /> 
     182 
    173183          <field id="ssu_m"    unit="m/s"  /> 
    174184          <field id="ssv_m"    unit="m/s"  /> 
    … …  
    218228          <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" /> 
    219229          <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" /> 
    220           <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" /> 
     230     <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" /> 
    221231          <field id="dmsspr"       long_name="snow mass change through snow fall"                      standard_name="snowfall_flux"                                                           unit="kg/m2/s" /> 
    222232          <field id="dmsmel"       long_name="snow mass change through melt"                           standard_name="surface_snow_melt_flux"                                                  unit="kg/m2/s" /> 
    … …  
    250260         <field id="xmtrpice_ave"     long_name="Monthly average of x-ice mass transport"   field_ref="xmtrpice"          grid_ref="grid_U_2D"        freq_op="1mo" freq_offset="_reset_"  > @xmtrpice </field> 
    251261         <field id="xmtrpice_section"                                                                                     grid_ref="grid_U_scalar"  > xmtrpice_ave </field> 
    252          <field id="xmtrpice_strait"                                                        field_ref="xmtrpice_section"  grid_ref="grid_U_4strait_ice"  /> 
    253          <field id="xstrait_mifl"                                                           field_ref="xmtrpice_strait"   grid_ref="grid_U_4strait_ice_hsum" unit="kg/s" detect_missing_value="true" > this * maskMFO_u_ice </field> 
    254262 
    255263         <field id="ymtrpice_ave"     long_name="Monthly average of y-ice mass transport"   field_ref="ymtrpice"          grid_ref="grid_V_2D"        freq_op="1mo" freq_offset="_reset_"  > @ymtrpice </field> 
    256264         <field id="ymtrpice_section"                                                                                     grid_ref="grid_V_scalar"  > ymtrpice_ave </field> 
    257          <field id="ymtrpice_strait"                                                        field_ref="ymtrpice_section"  grid_ref="grid_V_4strait_ice"  /> 
    258     <field id="ystrait_mifl"                                                           field_ref="ymtrpice_strait"   grid_ref="grid_V_4strait_ice_hsum" unit="kg/s" detect_missing_value="true"  > this * maskMFO_v_ice </field> 
    259265 
    260266         <field id="xmtrpsnw_ave"     long_name="Monthly average of x-snow mass transport"  field_ref="xmtrpsnw"          grid_ref="grid_U_2D"        freq_op="1mo" freq_offset="_reset_"  > @xmtrpsnw </field> 
    261267         <field id="xmtrpsnw_section"                                                                                     grid_ref="grid_U_scalar"  > xmtrpsnw_ave </field> 
    262          <field id="xmtrpsnw_strait"                                                        field_ref="xmtrpsnw_section"  grid_ref="grid_U_4strait_ice"  /> 
    263          <field id="xstrait_msfl"                                                           field_ref="xmtrpsnw_strait"   grid_ref="grid_U_4strait_ice_hsum" unit="kg/s" detect_missing_value="true" > this * maskMFO_u_ice </field> 
    264268 
    265269         <field id="ymtrpsnw_ave"     long_name="Monthly average of y-snow mass transport"  field_ref="ymtrpsnw"          grid_ref="grid_V_2D"        freq_op="1mo" freq_offset="_reset_"  > @ymtrpsnw </field> 
    266270         <field id="ymtrpsnw_section"                                                                                     grid_ref="grid_V_scalar"  > ymtrpsnw_ave </field> 
    267          <field id="ymtrpsnw_strait"                                                        field_ref="ymtrpsnw_section"  grid_ref="grid_V_4strait_ice"  /> 
    268     <field id="ystrait_msfl"                                                           field_ref="ymtrpsnw_strait"   grid_ref="grid_V_4strait_ice_hsum" unit="kg/s" detect_missing_value="true"  > this * maskMFO_v_ice </field> 
    269271 
    270272         <field id="xatrp_ave"        long_name="Monthly average of x-ice area transport"   field_ref="xatrp"             grid_ref="grid_U_2D"        freq_op="1mo" freq_offset="_reset_"  > @xatrp </field> 
    271273         <field id="xatrp_section"                                                                                        grid_ref="grid_U_scalar"  > xatrp_ave </field> 
    272          <field id="xatrp_strait"                                                           field_ref="xatrp_section"     grid_ref="grid_U_4strait_ice"  /> 
    273          <field id="xstrait_arfl"                                                           field_ref="xatrp_strait"      grid_ref="grid_U_4strait_ice_hsum" unit="kg/s" detect_missing_value="true" > this * maskMFO_u_ice </field> 
    274274 
    275275         <field id="yatrp_ave"        long_name="Monthly average of y-ice area transport"   field_ref="yatrp"             grid_ref="grid_V_2D"        freq_op="1mo" freq_offset="_reset_"  > @yatrp </field> 
    276276         <field id="yatrp_section"                                                                                        grid_ref="grid_V_scalar"  > yatrp_ave </field> 
    277          <field id="yatrp_strait"                                                           field_ref="yatrp_section"     grid_ref="grid_V_4strait_ice"  /> 
    278     <field id="ystrait_arfl"                                                           field_ref="yatrp_strait"      grid_ref="grid_V_4strait_ice_hsum" unit="m2/s" detect_missing_value="true"  > this * maskMFO_v_ice </field> 
    279  
    280          <field id="strait_mifl"      long_name="Sea ice mass flux through straits"      standard_name="sea_ice_mass_transport_across_line"   unit="kg/s"  freq_op="1mo"  grid_ref="grid_4strait_ice" > xstrait_mifl + ystrait_mifl </field> 
    281          <field id="strait_msfl"      long_name="Snow mass flux through straits"         standard_name="snow_mass_transport_across_line"      unit="kg/s"  freq_op="1mo"  grid_ref="grid_4strait_ice" > xstrait_msfl + ystrait_msfl </field> 
    282     <field id="strait_arfl"      long_name="Sea ice area flux through straits"      standard_name="sea_area_mass_transport_across_line"  unit="m2/s"  freq_op="1mo"  grid_ref="grid_4strait_ice" > xstrait_arfl + ystrait_arfl </field> 
     277 
    283278 
    284279   </field_group> <!-- SBC_2D --> 
    … …  
    295290          <field id="snwtemp_cat"  long_name="Snow temperature per category"                     unit="degC"    detect_missing_value="true" /> 
    296291          <field id="icettop_cat"  long_name="Ice/snow surface temperature per category"         unit="degC"    detect_missing_value="true" /> 
    297           <field id="iceapnd_cat"  long_name="Ice melt pond concentration per category"          unit=""        />  
     292          <field id="icevpnd_cat"  long_name="Ice melt pond volume per grid area per category"   unit="m"       />  
     293          <field id="iceapnd_cat"  long_name="Ice melt pond grid fraction"                       unit=""        />  
    298294          <field id="icehpnd_cat"  long_name="Ice melt pond thickness per category"              unit="m"       detect_missing_value="true" />  
    299295          <field id="icehlid_cat"  long_name="Ice melt pond lid thickness per category"          unit="m"       detect_missing_value="true" />  
    300           <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"               unit=""        />  
     296          <field id="iceafpnd_cat" long_name="Ice melt pond ice fraction per category"           unit=""        />  
    301297          <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category"     unit=""        />  
    302298          <field id="icemask_cat"  long_name="Fraction of time step with sea ice (per category)" unit=""        /> 
    … …  
    375371     <field field_ref="icehpnd"          name="sihpnd" /> 
    376372     <field field_ref="icevpnd"          name="sivpnd" /> 
     373     <field field_ref="iceepnd"          name="siepnd" /> 
    377374          <field field_ref="iceage"           name="siage"  /> 
    378      <field field_ref="sst_m"            name="sst_m"  /> 
    379      <field field_ref="sss_m"            name="sss_m"  /> 
     375     <field id="sst_m_pot"    unit="degC" /> 
     376     <!-- EOS-80 --> 
     377     <field id="sss_m_pra"    unit="psu"  /> 
     378     <!-- TEOS-10 --> 
     379          <field id="sss_m_abs"    unit="psu"  /> 
    380380      
    381381     <!-- heat --> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/SHARED/field_def_nemo-oce.xml

    r13648 r15658  
    1919       
    2020      <field_group id="grid_T" grid_ref="grid_T_2D" > 
    21         <field id="e3t"          long_name="T-cell thickness"                    standard_name="cell_thickness"        unit="m"   grid_ref="grid_T_3D" /> 
    22         <field id="e3ts"         long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"        unit="m"   grid_ref="grid_T_SFC"/> 
    23         <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness"    unit="m"   grid_ref="grid_T_3D" /> 
    24         <field id="e3tb"         long_name="bottom T-cell thickness"             standard_name="bottom_cell_thickness" unit="m"   grid_ref="grid_T_2D"/>  
     21        <field id="e3t"          long_name="T-cell thickness"                    standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_3D" /> 
     22        <field id="e3ts"     long_name="T-cell thickness"   field_ref="e3t"  standard_name="cell_thickness"     unit="m"   grid_ref="grid_T_SFC"/> 
     23        <field id="e3t_0"        long_name="Initial T-cell thickness"            standard_name="ref_cell_thickness" unit="m"   grid_ref="grid_T_3D" /> 
     24        <field id="e3tb"         long_name="bottom T-cell thickness"             standard_name="bottom_cell_thickness" unit="m"   grid_ref="grid_T_2D"/> 
    2525        <field id="e3t_300"      field_ref="e3t"                grid_ref="grid_T_zoom_300"       detect_missing_value="true" /> 
    26         <field id="e3t_vsum300"  field_ref="e3t_300"            grid_ref="grid_T_vsum"   detect_missing_value="true" /> 
     26   <field id="e3t_vsum300"  field_ref="e3t_300"            grid_ref="grid_T_vsum"   detect_missing_value="true" /> 
    2727   <field id="masscello"    long_name="Sea Water Mass per unit area"   standard_name="sea_water_mass_per_unit_area"   unit="kg/m2"   grid_ref="grid_T_3D"/> 
    28         <field id="volcello"     long_name="Ocean Volume"                   standard_name="ocean_volume"   unit="m3"       grid_ref="grid_T_3D"/>  
    29         <field id="toce"         long_name="temperature"                         standard_name="sea_water_potential_temperature"   unit="degC"     grid_ref="grid_T_3D"/> 
    30         <field id="toce_e3t"     long_name="temperature (thickness weighted)"                                                      unit="degC"     grid_ref="grid_T_3D" > toce * e3t </field > 
    31         <field id="soce"         long_name="salinity"                            standard_name="sea_water_practical_salinity"      unit="1e-3"     grid_ref="grid_T_3D"/> 
    32         <field id="soce_e3t"     long_name="salinity    (thickness weighted)"                                                      unit="1e-3"     grid_ref="grid_T_3D" > soce * e3t </field > 
     28        <field id="volcello"     long_name="Ocean Volume"                   standard_name="ocean_volume"   unit="m3"       grid_ref="grid_T_3D"/> 
     29 
     30        <!-- EOS80 --> 
     31        <field id="toce_pot"         long_name="potential temperature"           standard_name="sea_water_potential_temperature"      unit="degC"     grid_ref="grid_T_3D"/> 
     32        <field id="toce_pot_e3t"     long_name="potential temperature (thickness weighted)"                                           unit="degC"     grid_ref="grid_T_3D" > toce_pot * e3t </field > 
     33        <field id="soce_pra"         long_name="practical salinity"              standard_name="sea_water_practical_salinity"      unit="1e-3"     grid_ref="grid_T_3D"/> 
     34        <field id="soce_pra_e3t"     long_name="practical salinity    (thickness weighted)"                                        unit="1e-3"     grid_ref="grid_T_3D" > soce_pra * e3t </field > 
     35        <!-- TEOS10 --> 
     36        <field id="toce_con"         long_name="conservative temperature"        standard_name="sea_water_conservative_temperature"   unit="degC"     grid_ref="grid_T_3D"/> 
     37        <field id="toce_con_e3t"     long_name="conservative temperature (thickness weighted)"                                        unit="degC"     grid_ref="grid_T_3D" > toce_con * e3t </field >   
     38        <field id="soce_abs"         long_name="absolute salinity"               standard_name="sea_water_absolute_salinity"       unit="1e-3"     grid_ref="grid_T_3D"/> 
     39        <field id="soce_abs_e3t"     long_name="absolute salinity    (thickness weighted)"                                         unit="1e-3"     grid_ref="grid_T_3D" > soce_abs * e3t </field > 
    3340 
    3441        <field id="toce_e3t_300"      field_ref="toce_e3t"          unit="degree_C"     grid_ref="grid_T_zoom_300"      detect_missing_value="true" /> 
    … …  
    4956   <field id="ahmt_3d"      long_name=" 3D      t-eddy viscosity coefficient"   unit="m2/s or m4/s"  grid_ref="grid_T_3D"/> 
    5057 
    51         <field id="sst"          long_name="sea surface temperature"                            standard_name="sea_surface_temperature"             unit="degC"     /> 
    52         <field id="sst2"         long_name="square of sea surface temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst * sst </field > 
    53         <field id="sstmax"       long_name="max of sea surface temperature"   field_ref="sst"   operation="maximum"                                                 /> 
    54         <field id="sstmin"       long_name="min of sea surface temperature"   field_ref="sst"   operation="minimum"                                                 /> 
    55         <field id="sstgrad"      long_name="module of sst gradient"                                                                                 unit="degC/m"   /> 
    56         <field id="sstgrad2"     long_name="square of module of sst gradient"                                                                       unit="degC2/m2" /> 
    57         <field id="sbt"          long_name="sea bottom temperature"                                                                                 unit="degC"     /> 
    58         <field id="tosmint"      long_name="vertical integral of temperature times density"     standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature"  unit="(kg m2) degree_C" /> 
    59         <field id="sst_wl"       long_name="Delta SST of warm layer"                                                                                unit="degC"     /> 
    60         <field id="sst_cs"       long_name="Delta SST of cool skin"                                                                                 unit="degC"     /> 
    61    <field id="temp_3m"      long_name="temperature at 3m"                                                                                      unit="degC"     /> 
    62          
    63         <field id="sss"          long_name="sea surface salinity"                               standard_name="sea_surface_salinity"                unit="1e-3"     /> 
    64         <field id="sss2"         long_name="square of sea surface salinity"                                                                         unit="1e-6"      > sss * sss </field > 
    65         <field id="sssmax"       long_name="max of sea surface salinity"      field_ref="sss"   operation="maximum"                                                 /> 
    66         <field id="sssmin"       long_name="min of sea surface salinity"      field_ref="sss"   operation="minimum"                                                 /> 
    67         <field id="sbs"          long_name="sea bottom salinity"                                                                                    unit="0.001"    /> 
    68         <field id="somint"       long_name="vertical integral of salinity times density"        standard_name="integral_wrt_depth_of_product_of_density_and_salinity"  unit="(kg m2) x (1e-3)" />  
     58        <!-- EOS80 --> 
     59        <field id="sst_pot"          long_name="sea surface potential temperature"                            standard_name="sea_surface_temperature"             unit="degC"     /> 
     60        <field id="sst2_pot"         long_name="square of sea surface potential temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst_pot * sst_pot </field > 
     61        <field id="sstmax_pot"       long_name="max of sea surface potential temperature"   field_ref="sst_pot"   operation="maximum"                                                 /> 
     62        <field id="sstmin_pot"       long_name="min of sea surface potential temperature"   field_ref="sst_pot"   operation="minimum"                                                 /> 
     63        <field id="sstgrad_pot"      long_name="module of potential sst gradient"                                                                                 unit="degC/m"   /> 
     64        <field id="sstgrad2_pot"     long_name="square of module of potential sst gradient"                                                                       unit="degC2/m2" /> 
     65        <field id="sbt_pot"          long_name="sea bottom potential temperature"                                                                                 unit="degC"     /> 
     66        <field id="tosmint_pot"      long_name="vertical integral of potential temperature times density"     standard_name="integral_wrt_depth_of_product_of_density_and_potential_temperature"  unit="(kg m2) degree_C" /> 
     67        <field id="sst_wl_pot"       long_name="Delta potential SST of warm layer"                                                                                unit="degC"     /> 
     68        <field id="sst_cs_pot"       long_name="Delta potential SST of cool skin"                                                                                 unit="degC"     /> 
     69   <field id="temp_3m_pot"      long_name="potential temperature at 3m"                                                                                      unit="degC"     /> 
     70 
     71        <field id="sss_pra"          long_name="sea surface practical salinity"                               standard_name="sea_surface_practical_salinity"                unit="1e-3"     /> 
     72        <field id="sss2_pra"         long_name="square of sea surface practical salinity"                                                                         unit="1e-6"      > sss_pra * sss_pra </field > 
     73        <field id="sssmax_pra"       long_name="max of sea surface practical salinity"      field_ref="sss_pra"   operation="maximum"                                                 /> 
     74        <field id="sssmin_pra"       long_name="min of sea surface practical salinity"      field_ref="sss_pra"   operation="minimum"                                                 /> 
     75        <field id="sbs_pra"          long_name="sea bottom practical salinity"                                                                                    unit="0.001"    /> 
     76        <field id="somint_pra"       long_name="vertical integral of practical salinity times density"        standard_name="integral_wrt_depth_of_product_of_density_and_practical_salinity"  unit="(kg m2) x (1e-3)" />  
     77        <!-- TEOS10 --> 
     78        <field id="sst_con"          long_name="sea surface conservative temperature"                            standard_name="sea_surface_conservative_temperature"             unit="degC"     /> 
     79        <field id="sst2_con"         long_name="square of sea surface conservative temperature"                  standard_name="square_of_sea_surface_temperature"   unit="degC2"     > sst_con * sst_con </field > 
     80        <field id="sstmax_con"       long_name="max of sea surface conservative temperature"   field_ref="sst_con"   operation="maximum"                                                 /> 
     81        <field id="sstmin_con"       long_name="min of sea surface conservative temperature"   field_ref="sst_con"   operation="minimum"                                                 /> 
     82        <field id="sstgrad_con"      long_name="module of conservative sst gradient"                                                                                 unit="degC/m"   /> 
     83        <field id="sstgrad2_con"     long_name="square of module of conservative sst gradient"                                                                       unit="degC2/m2" /> 
     84        <field id="sbt_con"          long_name="sea bottom conservative temperature"                                                                                 unit="degC"     /> 
     85        <field id="tosmint_con"      long_name="vertical integral of conservative temperature times density"     standard_name="integral_wrt_depth_of_product_of_density_and_conservative_temperature"  unit="(kg m2) degree_C" /> 
     86        <field id="sst_wl_con"       long_name="Delta conservative SST of warm layer"                                                                                unit="degC"     /> 
     87        <field id="sst_cs_con"       long_name="Delta conservative SST of cool skin"                                                                                 unit="degC"     /> 
     88   <field id="temp_3m_con"      long_name="conservative temperature at 3m"                                                                                      unit="degC"     /> 
     89 
     90        <field id="sss_abs"          long_name="sea surface absolute salinity"                               standard_name="sea_surface_absolute_salinity"                unit="1e-3"     /> 
     91        <field id="sss2_abs"         long_name="square of sea surface absolute salinity"                                                                         unit="1e-6"      > sss_abs * sss_abs </field > 
     92        <field id="sssmax_abs"       long_name="max of sea surface absolute salinity"      field_ref="sss_abs"   operation="maximum"                                                 /> 
     93        <field id="sssmin_abs"       long_name="min of sea surface absolute salinity"      field_ref="sss_abs"   operation="minimum"                                                 /> 
     94        <field id="sbs_abs"          long_name="sea bottom absolute salinity"                                                                                    unit="0.001"    /> 
     95        <field id="somint_abs"       long_name="vertical integral of absolute salinity times density"        standard_name="integral_wrt_depth_of_product_of_density_and_absolute_salinity"  unit="(kg m2) x (1e-3)" /> 
    6996 
    7097        <field id="taubot"       long_name="bottom stress module"                                                                                   unit="N/m2"     />  
    … …  
    85112        <field id="mldr10_1max"  long_name="Max of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="maximum"                                                                          /> 
    86113        <field id="mldr10_1min"  long_name="Min of Mixed Layer Depth (dsigma = 0.01 wrt 10m)"   field_ref="mldr10_1"   operation="minimum"                                                                          /> 
     114         <field id="mldzint_1"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     115         <field id="mldzint_2"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     116         <field id="mldzint_3"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     117         <field id="mldzint_4"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     118         <field id="mldzint_5"    long_name="Mixed Layer Depth interpolated"                     standard_name="ocean_mixed_layer_thickness"                                                       unit="m"          /> 
     119         <field id="mldhtc_1"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     120         <field id="mldhtc_2"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     121         <field id="mldhtc_3"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     122         <field id="mldhtc_4"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
     123         <field id="mldhtc_5"     long_name="Mixed Layer Depth integrated heat content"          standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    87124        <field id="heatc"        long_name="Heat content vertically integrated"                 standard_name="integral_of_sea_water_potential_temperature_wrt_depth_expressed_as_heat_content"   unit="J/m2"       /> 
    88125        <field id="saltc"        long_name="Salt content vertically integrated"                                                                                                                   unit="1e-3*kg/m2" /> 
    … …  
    92129        <field id="beta"         long_name="haline contraction"                                                        unit="1e3"    grid_ref="grid_T_3D" /> 
    93130        <field id="rhop"         long_name="potential density (sigma0)"        standard_name="sea_water_sigma_theta"   unit="kg/m3"  grid_ref="grid_T_3D" /> 
     131        <field id="N_2d"         long_name="Depth-mean N"                                                              unit="m/s" /> 
     132        <field id="modslp"       long_name="sqrt( slpi^2 + slpj^2 )"                                                   unit="1"      grid_ref="grid_T_3D" /> 
     133        <field id="RossRad"      long_name="Rossby radius"                                                             unit="m" /> 
     134        <field id="RossRadlim"   long_name="Rossby radius (limited)"                                                   unit="m" /> 
     135        <field id="Tclinic_recip" long_name="recip of baroclinic timescale"                                            unit="s-1" /> 
     136        <field id="RR_GS"        long_name="Rossby radius / min(dx,dy)"                                                unit="1" /> 
    94137 
    95138        <!-- Energy - horizontal divergence --> 
    … …  
    360403      
    361404     <!-- sbcssm variables --> 
    362           <field id="sst_m"    unit="degC" /> 
    363           <field id="sss_m"    unit="psu"  /> 
     405          <field id="sst_m_pot"    unit="degC" /> 
     406     <!-- EOS-80 --> 
     407          <field id="sss_m_pra"    unit="psu"  /> 
     408          <!-- TEOS-10 --> 
     409          <field id="sss_m_abs"    unit="psu"  /> 
     410 
    364411          <field id="ssu_m"    unit="m/s"  /> 
    365412          <field id="ssv_m"    unit="m/s"  /> 
    … …  
    393440 
    394441 
     442        <field id="uoce2_e3u"    long_name="ocean current along i-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_U_3D"  > uoce * uoce * e3u </field> 
    395443        <field id="ssu"          long_name="ocean surface current along i-axis"                                                                 unit="m/s"                             /> 
    396444        <field id="sbu"          long_name="ocean bottom current along i-axis"                                                                  unit="m/s"                             /> 
    … …  
    450498        <field id="voce"         long_name="ocean current along j-axis"                             standard_name="sea_water_y_velocity"        unit="m/s"        grid_ref="grid_V_3D" /> 
    451499        <field id="voce_e3v"     long_name="ocean current along j-axis  (thickness weighted)"                                                   unit="m/s"        grid_ref="grid_V_3D"  > voce * e3v </field> 
     500        <field id="voce2_e3v"    long_name="ocean current along j-axis squared (thickness weighted)"                                            unit="m3/s2"      grid_ref="grid_V_3D"  > voce * voce * e3v </field> 
    452501        <field id="ssv"          long_name="ocean surface current along j-axis"                                                                 unit="m/s"                             /> 
    453502        <field id="sbv"          long_name="ocean bottom current along j-axis"                                                                  unit="m/s"                             /> 
    … …  
    551600      <field id="ahmf_2d"      long_name=" surface f-eddy viscosity coefficient"   unit="m2/s or m4/s" /> 
    552601      <field id="ahmf_3d"      long_name=" 3D      f-eddy viscosity coefficient"   unit="m2/s or m4/s"  grid_ref="grid_T_3D"/> 
     602 
     603      <!-- product fields --> 
     604      <field_group id="diaprod"> 
     605   <field id="ut"           long_name="product_of_sea_water_x_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_U_3D"   /> 
     606        <field id="ut_e3u"       long_name="product_of_sea_water_x_velocity_and_potential_temperature * e3u"  unit="degree_C m2/s"   grid_ref="grid_U_3D" > ut * e3u </field > 
     607   <field id="us"           long_name="product_of_sea_water_x_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_U_3D"   /> 
     608        <field id="us_e3u"       long_name="product_of_sea_water_x_velocity_and_salinity * e3u"             unit="PSU m2/s"      grid_ref="grid_U_3D" > us * e3u </field > 
     609   <field id="urhop"        long_name="product_of_sea_water_x_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_U_3D"   /> 
     610        <field id="urhop_e3u"    long_name="product_of_sea_water_x_velocity_and_potential_density * e3u"    unit="(kg/m3).(m2/s)"   grid_ref="grid_U_3D" > urhop * e3u </field > 
     611   <field id="vt"           long_name="product_of_sea_water_y_velocity_and_potential_temperature"      unit="degree_C m/s"      grid_ref="grid_V_3D"   /> 
     612        <field id="vt_e3v"       long_name="product_of_sea_water_y_velocity_and_potential_temperature * e3v"  unit="degree_C m2/s"   grid_ref="grid_V_3D" > vt * e3v </field > 
     613   <field id="vs"           long_name="product_of_sea_water_y_velocity_and_salinity"                   unit="PSU m/s"       grid_ref="grid_V_3D"   /> 
     614        <field id="vs_e3v"       long_name="product_of_sea_water_y_velocity_and_salinity * e3t"             unit="PSU m2/s"      grid_ref="grid_V_3D" > vs * e3v </field > 
     615   <field id="vrhop"        long_name="product_of_sea_water_y_velocity_and_potential_density"          unit="(kg/m3).(m/s)" grid_ref="grid_V_3D"   /> 
     616        <field id="vrhop_e3v"    long_name="product_of_sea_water_y_velocity_and_potential_density * e3t"    unit="(kg/m3).(m2/s)"  grid_ref="grid_V_3D" > vrhop * e3v </field > 
     617   <field id="wt"           long_name="product_of_upward_sea_water_velocity_and_potential_temperature" unit="degree_C m/s"      grid_ref="grid_W_3D"   /> 
     618   <field id="ws"           long_name="product_of_upward_sea_water_velocity_and_salinity"              unit="PSU m/s"       grid_ref="grid_W_3D"   /> 
     619   <field id="wrhop"        long_name="product_of_upward_sea_water_velocity_and_potential_density"     unit="(kg/m3).(m/s)" grid_ref="grid_W_3D"   /> 
     620   <field id="uv"           long_name="product_of_sea_water_x_velocity_and_sea_water_y_velocity"       unit="m2/s2   "      grid_ref="grid_T_3D"   /> 
     621   <field id="uw"           long_name="product_of_upward_sea_water_velocity_and_sea_water_x_velocity"  unit="m2/s2   "      grid_ref="grid_W_3D"   /> 
     622   <field id="vw"           long_name="product_of_upward_sea_water_velocity_and_sea_water_y_velocity"  unit="m2/s2"         grid_ref="grid_W_3D"   /> 
     623      </field_group> 
    553624 
    554625      <field_group id="scalar"  grid_ref="grid_scalar"  > 
    … …  
    587658         <field id="uoce_e3u_ave_vsum"    long_name="Vertical sum of u*e3u"                           field_ref="uoce_e3u_ave"         grid_ref="grid_U_vsum"     /> 
    588659         <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> 
    589          <field id="uocetr_strait"        long_name="Total transport across lines in i-direction"     field_ref="uocetr_vsum_section"  grid_ref="grid_U_4strait" /> 
    590          <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> 
    591660 
    592661         <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> 
    593662         <field id="voce_e3v_ave_vsum"    long_name="Vertical sum of v*e3v"                           field_ref="voce_e3v_ave"         grid_ref="grid_V_vsum"      /> 
    594663         <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> 
    595          <field id="vocetr_strait"        long_name="Total transport across lines in j-direction"     field_ref="vocetr_vsum_section"  grid_ref="grid_V_4strait"  /> 
    596          <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> 
    597  
    598          <field id="masstr_strait"        long_name="Sea water transport across line"                                                  grid_ref="grid_4strait"  > u_masstr_strait + v_masstr_strait </field> 
     664 
    599665      </field_group> 
    600666 
    … …  
    897963     
    898964    <field_group id="mooring" > 
    899       <field field_ref="toce"         name="thetao"   long_name="sea_water_potential_temperature"      /> 
    900       <field field_ref="soce"         name="so"       long_name="sea_water_salinity"                   /> 
     965      <!-- EOS80 --> 
     966      <field field_ref="toce_pot"         name="thetao_pot"   long_name="sea_water_potential_temperature"      /> 
     967      <field field_ref="soce_pra"         name="so_pra"       long_name="sea_water_practical_salinity"                   /> 
     968      <!-- TEOS10 --> 
     969      <field field_ref="toce_con"         name="thetao_con"   long_name="sea_water_conservative_temperature"      /> 
     970      <field field_ref="soce_abs"         name="so_con"       long_name="sea_water_absolute_salinity"                   /> 
     971 
    901972      <field field_ref="uoce"         name="uo"       long_name="sea_water_x_velocity"                 /> 
    902973      <field field_ref="voce"         name="vo"       long_name="sea_water_y_velocity"                 /> 
    … …  
    904975      <field field_ref="avt"          name="difvho"   long_name="ocean_vertical_heat_diffusivity"      /> 
    905976      <field field_ref="avm"          name="difvmo"   long_name="ocean_vertical_momentum_diffusivity"  /> 
    906        
    907       <field field_ref="sst"          name="tos"      long_name="sea_surface_temperature"                       /> 
    908       <field field_ref="sst2"         name="tossq"    long_name="square_of_sea_surface_temperature"             /> 
    909       <field field_ref="sstgrad"      name="tosgrad"  long_name="module_of_sea_surface_temperature_gradient"    /> 
    910       <field field_ref="sss"          name="sos"      long_name="sea_surface_salinity"                          /> 
     977 
     978      <!-- EOS80 --> 
     979      <field field_ref="sst_pot"          name="tos_pot"      long_name="sea_surface_potential_temperature"                       /> 
     980      <field field_ref="sst2_pot"         name="tossq_pot"    long_name="square_of_sea_surface_potential_temperature"             /> 
     981      <field field_ref="sstgrad_pot"      name="tosgrad_pot"  long_name="module_of_sea_surface_potential_temperature_gradient"    /> 
     982      <field field_ref="sss_pra"          name="sos_pra"      long_name="sea_surface_absolute_salinity"                          /> 
     983      <!-- TEOS10 --> 
     984      <field field_ref="sst_con"          name="tos_con"      long_name="sea_surface_conservative_temperature"                       /> 
     985      <field field_ref="sst2_con"         name="tossq_con"    long_name="square_of_sea_surface_conservative_temperature"             /> 
     986      <field field_ref="sstgrad_con"      name="tosgrad_con"  long_name="module_of_sea_surface_conservative_temperature_gradient"    /> 
     987      <field field_ref="sss_abs"          name="sos"      long_name="sea_surface_absolute_salinity"                          /> 
     988 
    911989      <field field_ref="ssh"          name="zos"      long_name="sea_surface_height_above_geoid"                /> 
    912990      <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     /> 
    … …  
    9311009 
    9321010    <field_group id="groupT" > 
    933       <field field_ref="toce"         name="thetao"   long_name="sea_water_potential_temperature"               /> 
    934       <field field_ref="soce"         name="so"       long_name="sea_water_salinity"                            /> 
    935       <field field_ref="sst"          name="tos"      long_name="sea_surface_temperature"                       /> 
    936       <field field_ref="sst2"         name="tossq"    long_name="square_of_sea_surface_temperature"             /> 
    937       <field field_ref="sss"          name="sos"      long_name="sea_surface_salinity"                          /> 
     1011      <!-- EOS80 --> 
     1012      <field field_ref="toce_pot"         name="thetao_pot"   long_name="sea_water_potential_temperature"               /> 
     1013      <field field_ref="soce_pra"         name="so_pra"       long_name="sea_water_practical_salinity"                            /> 
     1014      <field field_ref="sst_pot"          name="tos_pot"      long_name="sea_surface_potential_temperature"                       /> 
     1015      <field field_ref="sst2_pot"         name="tossq_pot"    long_name="square_of_sea_surface_potential_temperature"             /> 
     1016      <field field_ref="sss_pra"          name="sos_pra"      long_name="sea_surface_practical_salinity"                          /> 
     1017      <!-- TEOS10 --> 
     1018      <field field_ref="toce_con"         name="thetao_con"   long_name="sea_water_conservative_temperature"               /> 
     1019      <field field_ref="soce_abs"         name="so_abs"       long_name="sea_water_absolute_salinity"                            /> 
     1020      <field field_ref="sst_con"          name="tos_con"      long_name="sea_surface_conservative_temperature"                       /> 
     1021      <field field_ref="sst2_con"         name="tossq_con"    long_name="square_of_sea_surface_conservative_temperature"             /> 
     1022      <field field_ref="sss_abs"          name="sos_abs"      long_name="sea_surface_absolute_salinity"                          /> 
     1023 
    9381024      <field field_ref="ssh"          name="zos"      long_name="sea_surface_height_above_geoid"                /> 
    9391025      <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     /> 
    … …  
    9671053    </field_group> 
    9681054 
     1055    <!-- TMB diagnostic output --> 
     1056    <field_group  id="1h_grid_T_tmb" grid_ref="grid_T_2D" operation="instant"> 
     1057      <!-- EOS80 --> 
     1058      <field id="top_temp_pot"           name="votemper_top_pot"  unit="degC"  /> 
     1059      <field id="mid_temp_pot"           name="votemper_mid_pot"  unit="degC"  /> 
     1060      <field id="bot_temp_pot"           name="votemper_bot_pot"  unit="degC"  /> 
     1061      <field id="top_sal_pra"            name="vosaline_top_pra"  unit="psu"   /> 
     1062      <field id="mid_sal_pra"            name="vosaline_mid_pra"  unit="psu"   /> 
     1063      <field id="bot_sal_pra"            name="vosaline_bot_pra"  unit="psu"   /> 
     1064      <!-- TEOS10 --> 
     1065      <field id="top_temp_con"           name="votemper_top_con"  unit="degC"  /> 
     1066      <field id="mid_temp_con"           name="votemper_mid_con"  unit="degC"  /> 
     1067      <field id="bot_temp_con"           name="votemper_bot_con"  unit="degC"  /> 
     1068      <field id="top_sal_abs"            name="vosaline_top_abs"  unit="psu"   /> 
     1069      <field id="mid_sal_abs"            name="vosaline_mid_abs"  unit="psu"   /> 
     1070      <field id="bot_sal_abs"            name="vosaline_bot_abs"  unit="psu"   /> 
     1071 
     1072      <field id="sshnmasked"         name="sossheig"      unit="m"     />  
     1073    </field_group> 
     1074 
    9691075    <field_group  id="1h_grid_U_tmb" grid_ref="grid_U_2D" operation="instant"> 
    9701076      <field id="top_u"           name="vozocrtx_top"  unit="m/s"  /> 
    … …  
    9831089    <!-- 25h diagnostic output --> 
    9841090    <field_group id="25h_grid_T" grid_ref="grid_T_3D" operation="instant"> 
    985       <field id="temper25h"         name="potential temperature 25h mean"    unit="degC" /> 
     1091      <!-- EOS80 --> 
     1092      <field id="temper25h_pot"         name="potential temperature 25h mean"    unit="degC" /> 
    9861093      <field id="tempis25h"         name="insitu temperature 25h mean"    unit="degC" /> 
    987       <field id="salin25h"          name="salinity 25h mean"                 unit="psu"  /> 
     1094      <field id="salin25h_pra"          name="practical salinity 25h mean"                 unit="psu"  /> 
     1095      <!-- TEOS10 --> 
     1096      <field id="temper25h_con"         name="conservative temperature 25h mean"    unit="degC" /> 
     1097      <field id="tempis25h"         name="insitu temperature 25h mean"    unit="degC" /> 
     1098      <field id="salin25h_abs"          name="absolute salinity 25h mean"                 unit="psu"  /> 
     1099 
    9881100      <field id="ssh25h"            name="sea surface height 25h mean"  grid_ref="grid_T_2D"      unit="m"    /> 
    9891101    </field_group> 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/SHARED/grid_def_nemo.xml

    r12331 r15658  
    2020         <domain domain_ref="grid_T" /> 
    2121         <axis axis_ref="deptht" /> 
    22        </grid> 
    23         <!--  --> 
    24        <grid id="grid_T_3DS" > 
    25          <domain domain_ref="grid_T" /> 
    26          <axis axis_ref="profsed" /> 
    2722       </grid> 
    2823        <!--  --> 
    … …  
    6762       </grid> 
    6863        <!--  --> 
    69  
    70  
    71        <grid id="grid_znl_T_2D"> 
    72          <domain domain_ref="gznl" /> 
    73          <axis axis_ref="basin" /> 
    74        </grid> 
    75  
    76        <grid id="grid_znl_T_3D"> 
    77          <domain domain_ref="gznl" /> 
    78          <axis axis_ref="deptht"  /> 
    79          <axis axis_ref="basin" /> 
    80        </grid> 
    81  
    82        <grid id="grid_znl_W_3D"> 
    83          <domain domain_ref="gznl" /> 
    84          <axis axis_ref="depthw"  /> 
    85          <axis axis_ref="basin" /> 
    86        </grid> 
    87  
    88       <grid id="grid_ptr_T_2D"> 
    89          <domain domain_ref="ptr" /> 
    90          <axis axis_ref="basin" /> 
    91        </grid> 
    92  
    93        <grid id="grid_ptr_T_3D"> 
    94          <domain  domain_ref="ptr" /> 
    95          <axis axis_ref="deptht"  /> 
    96          <axis axis_ref="basin" /> 
    97        </grid> 
    98  
    99        <grid id="grid_ptr_W_3D"> 
    100          <domain  domain_ref="ptr" /> 
    101          <axis axis_ref="depthw"  /> 
    102          <axis axis_ref="basin" /> 
    103        </grid> 
    104  
    10564       <grid id="grid_ptr_W_GLO"> 
    10665         <domain  domain_ref="ptr" /> 
    … …  
    155114       </grid> 
    156115 
    157        <!-- for ORCA2 grid  --> 
    158        <grid id="cumul_U"> 
    159          <axis axis_ref="cumul_U" n_glo="182" > 
    160            <reduce_domain local="true" operation="sum" direction="jDir" /> 
    161            <reduce_axis operation="sum" /> 
    162          </axis> 
    163          <axis axis_ref="depthu" /> 
    164        </grid> 
    165  
    166        <!-- for eORCA1 grid 
    167  
    168        <grid id="cumul_U"> 
    169          <axis axis_ref="cumul_U" n_glo="362" > 
    170            <reduce_domain local="true" operation="sum" direction="jDir" /> 
    171            <reduce_axis operation="sum" /> 
    172          </axis> 
    173          <axis axis_ref="depthu" /> 
    174        </grid> 
    175  
    176       --> 
    177  
    178116 
    179117       <grid id="grid_T_zoom_300"> 
    … …  
    192130       </grid> 
    193131  
    194        <grid id="grid_U_4strait"> 
    195          <domain domain_ref="grid_U" /> 
    196          <axis axis_ref="section"> 
    197            <duplicate_scalar/> 
    198          </axis> 
    199        </grid> 
    200   
    201        <grid id="grid_V_4strait"> 
    202          <domain domain_ref="grid_V" /> 
    203          <axis axis_ref="section"> 
    204            <duplicate_scalar/> 
    205          </axis> 
    206        </grid> 
    207   
    208        <grid id="grid_U_4strait_hsum"> 
    209          <scalar > 
    210            <reduce_domain operation="sum" local="true"/> 
    211            <reduce_scalar operation="sum" /> 
    212          </scalar> 
    213          <axis axis_ref="section"/> 
    214        </grid> 
    215   
    216        <grid id="grid_V_4strait_hsum"> 
    217          <scalar > 
    218            <reduce_domain operation="sum" local="true"/> 
    219            <reduce_scalar operation="sum" /> 
    220          </scalar> 
    221          <axis axis_ref="section"/> 
    222        </grid> 
    223   
    224        <grid id="grid_4strait"> 
    225          <axis axis_ref="section"/> 
    226        </grid> 
    227  
    228       <grid id="grid_U_4strait_ice"> 
    229         <domain domain_ref="grid_U" /> 
    230         <axis axis_ref="section_ice"> 
    231           <duplicate_scalar/> 
    232         </axis> 
    233       </grid> 
    234  
    235       <grid id="grid_V_4strait_ice"> 
    236         <domain domain_ref="grid_V" /> 
    237         <axis axis_ref="section_ice"> 
    238          <duplicate_scalar/> 
    239         </axis> 
    240       </grid> 
    241  
    242       <grid id="grid_U_4strait_ice_hsum"> 
    243         <scalar > 
    244          <reduce_domain operation="sum" local="true"/> 
    245          <reduce_scalar operation="sum" /> 
    246         </scalar> 
    247         <axis axis_ref="section_ice"/> 
    248       </grid> 
    249  
    250      <grid id="grid_V_4strait_ice_hsum"> 
    251         <scalar > 
    252          <reduce_domain operation="sum" local="true"/> 
    253          <reduce_scalar operation="sum" /> 
    254         </scalar> 
    255         <axis axis_ref="section_ice"/> 
    256      </grid> 
    257  
    258      <grid id="grid_4strait_ice"> 
    259        <axis axis_ref="section_ice"/> 
    260      </grid> 
    261   
    262132      <!-- scalars --> 
    263133      <grid id="grid_scalar" > 
    264134        <scalar/> 
    265135      </grid> 
    266   
    267     </grid_definition>    
    268   
     136        <!--  --> 
     137       <grid id="grid_EqT" > 
     138         <domain id="EqT" /> 
     139       </grid> 
     140        <!--  --> 
     141       <grid id="gznl_T_2D"> 
     142         <domain id="ptr" /> 
     143       </grid> 
     144        <!--  --> 
     145       <grid id="gznl_T_3D"> 
     146         <domain id="ptr" /> 
     147         <axis axis_ref="deptht" /> 
     148       </grid> 
     149        <!--  --> 
     150       <grid id="gznl_W_2D"> 
     151         <domain id="ptr" /> 
     152       </grid> 
     153        <!--  --> 
     154       <grid id="gznl_W_3D"> 
     155         <domain id="ptr" /> 
     156         <axis axis_ref="depthw" /> 
     157       </grid> 
     158       <grid id="vert_sum"> 
     159         <domain id="grid_T"/> 
     160         <scalar> 
     161            <reduce_axis operation="sum" /> 
     162         </scalar> 
     163       </grid> 
     164       <grid id="zoom_300"> 
     165         <domain id="grid_T" /> 
     166         <axis axis_ref="deptht300"/> 
     167       </grid> 
     168       <grid id="zoom_300_sum"> 
     169         <domain id="grid_T" /> 
     170         <scalar> 
     171            <reduce_axis operation="sum" /> 
     172         </scalar> 
     173       </grid> 
     174    </grid_definition> 
     175     

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/SHARED/namelist_ice_ref

    r13346 r15658  
    195195!------------------------------------------------------------------------------ 
    196196   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 ?? 
     197      ln_pnd_TOPO    = .true.         !  topographic melt ponds  
     198      ln_pnd_LEV     = .false.        !  level ice melt ponds  
     199         rn_apnd_min =   0.15         !     minimum meltwater fraction contributing to pond growth (TOPO and LEV) 
     200         rn_apnd_max =   0.85         !     maximum meltwater fraction contributing to pond growth (TOPO and LEV) 
     201         rn_pnd_flush=   0.1          !     pond flushing efficiency (tuning parameter) (LEV) 
    200202      ln_pnd_CST     = .false.        !  constant  melt ponds 
    201203         rn_apnd     =   0.2          !     prescribed pond fraction, at Tsu=0 degC 
    202204         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) 
     205      ln_pnd_lids    = .false.        !  frozen lids on top of the ponds (only for ln_pnd_LEV) 
    204206      ln_pnd_alb     = .true.         !  effect of melt ponds on ice albedo 
    205207/ 
    … …  
    226228   rn_tms_ini_n     = 270.            !  initial snw temperature     (K), North 
    227229   rn_tms_ini_s     = 270.            !        "            "             South 
    228    rn_apd_ini_n     =   0.2           !  initial pond fraction       (-), North 
     230   rn_apd_ini_n     =   0.0           !  initial pond fraction       (-), North 
    229231   rn_apd_ini_s     =   0.2           !        "            "             South 
    230    rn_hpd_ini_n     =   0.05          !  initial pond depth          (m), North 
     232   rn_hpd_ini_n     =   0.00          !  initial pond depth          (m), North 
    231233   rn_hpd_ini_s     =   0.05          !        "            "             South 
    232234   rn_hld_ini_n     =   0.0           !  initial pond lid depth      (m), North 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/cfgs/SHARED/namelist_ref

    r14075 r15658  
    6767   ln_xios_read = .FALSE.  !  use XIOS to read restart file (only for a single file restart) 
    6868   nn_wxios = 0      !  use XIOS to write restart file 0 - no, 1 - single file output, 2 - multiple file output 
     69   ln_rst_eos = .TRUE.     ! check if the equation of state used to produce the restart is consistent with model 
    6970/ 
    7071!----------------------------------------------------------------------- 
    … …  
    745746&nameos        !   ocean Equation Of Seawater                           (default: NO selection) 
    746747!----------------------------------------------------------------------- 
    747    ln_teos10   = .false.         !  = Use TEOS-10 
     748   ln_teos10   = .true.         !  = Use TEOS-10 
    748749   ln_eos80    = .false.         !  = Use EOS80 
    749750   ln_seos     = .false.         !  = Use S-EOS (simplified Eq.) 
    … …  
    837838      !                             !   = 30 F(i,j,k)  =ldf_c2d * ldf_c1d 
    838839      !                        !  time invariant coefficients:  aei0 = 1/2  Ue*Le  
    839       rn_Ue        = 0.02           !  lateral diffusive velocity [m/s] (nn_aht_ijk_t= 0, 10, 20, 30) 
    840       rn_Le        = 200.e+3        !  lateral diffusive length   [m]   (nn_aht_ijk_t= 0, 10) 
     840      rn_Ue        = 0.02           !  lateral diffusive velocity [m/s] (nn_aei_ijk_t= 0, 10, 20, 30) 
     841      rn_Le        = 200.e+3        !  lateral diffusive length   [m]   (nn_aei_ijk_t= 0, 10) 
    841842      ! 
    842843      ln_ldfeiv_dia =.false.   ! diagnose eiv stream function and velocities 
     844      nn_ldfeiv_shape = 0           !  shape of bounding coefficient    (nn_aei_ijk_t= 21 only) 
    843845/ 
    844846!----------------------------------------------------------------------- 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/ice.F90

    r14075 r15658  
    208208   !                                     !!** ice-ponds namelist (namthd_pnd) 
    209209   LOGICAL , PUBLIC ::   ln_pnd           !: Melt ponds (T) or not (F) 
    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 
     210   LOGICAL , PUBLIC ::   ln_pnd_TOPO      !: Topographic Melt ponds scheme (Flocco et al 2007, 2010) 
     211   LOGICAL , PUBLIC ::   ln_pnd_LEV       !: Simple melt pond scheme 
     212   REAL(wp), PUBLIC ::   rn_apnd_min      !: Minimum fraction of melt water contributing to ponds 
     213   REAL(wp), PUBLIC ::   rn_apnd_max      !: Maximum fraction of melt water contributing to ponds 
     214   REAL(wp), PUBLIC ::   rn_pnd_flush     !: Pond flushing efficiency (tuning parameter) 
    213215   LOGICAL , PUBLIC ::   ln_pnd_CST       !: Melt ponds scheme with constant fraction and depth 
    214216   REAL(wp), PUBLIC ::   rn_apnd          !: prescribed pond fraction (0<rn_apnd<1) 
    … …  
    308310   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   t1_ice          !: temperature of the first layer          (ln_cndflx=T) [K] 
    309311   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   cnd_ice         !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 
     312 
     313   ! meltwater arrays to save for melt ponds 
     314   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_i_sum_2d   !: surface melt (2d arrays for ponds)       [m] 
     315   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dh_s_mlt_2d   !: snow surf melt (2d arrays for ponds)     [m] 
    310316 
    311317   !!---------------------------------------------------------------------- 
    … …  
    458464      ii = ii + 1 
    459465      ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) ,  & 
     466         &      dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) ,                    & 
    460467         &      h_i        (jpi,jpj,jpl) , a_i    (jpi,jpj,jpl) , v_i   (jpi,jpj,jpl) ,  & 
    461468         &      v_s        (jpi,jpj,jpl) , h_s    (jpi,jpj,jpl) , t_su  (jpi,jpj,jpl) ,  & 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icedyn_adv_pra.F90

    r14075 r15658  
    178178               z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice  heat content 
    179179            END DO 
    180             IF ( ln_pnd_LEV ) THEN 
     180            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    181181               z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond fraction 
    182182               z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:)      ! Melt pond volume 
    … …  
    214214            END DO 
    215215            ! 
    216             IF ( ln_pnd_LEV ) THEN 
     216            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    217217               CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    218218               CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )  
    … …  
    249249                  &                                 sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 
    250250            END DO 
    251             IF ( ln_pnd_LEV ) THEN 
     251            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    252252               CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap )    !--- melt pond fraction 
    253253               CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) 
    … …  
    278278         CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei  , 'T', 1._wp, sxe   , 'T', -1._wp, sye   , 'T', -1._wp  & ! ice enthalpy 
    279279            &                                , sxxe  , 'T', 1._wp, syye  , 'T',  1._wp, sxye  , 'T',  1._wp  ) 
    280          IF ( ln_pnd_LEV ) THEN 
     280         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    281281            CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp  & ! melt pond fraction 
    282282               &                                , sxxap, 'T', 1._wp, syyap, 'T',  1._wp, sxyap, 'T',  1._wp  & 
    … …  
    302302               pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    303303            END DO 
    304             IF ( ln_pnd_LEV ) THEN 
     304            IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    305305               pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    306306               pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 
    … …  
    780780                  !                               ! -- check h_ip -- ! 
    781781                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    782                   IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     782                  IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    783783                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    784784                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 
    … …  
    10271027            END DO 
    10281028            ! 
    1029             IF( ln_pnd_LEV ) THEN                                    ! melt pond fraction 
     1029            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                    ! melt pond fraction 
    10301030               IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 
    10311031                  CALL iom_get( numrir, jpdom_autoglo, 'sxap' , sxap  ) 
    … …  
    10691069            sxc0  = 0._wp   ;   syc0  = 0._wp   ;   sxxc0  = 0._wp   ;   syyc0  = 0._wp   ;   sxyc0  = 0._wp      ! snow layers heat content 
    10701070            sxe   = 0._wp   ;   sye   = 0._wp   ;   sxxe   = 0._wp   ;   syye   = 0._wp   ;   sxye   = 0._wp      ! ice layers heat content 
    1071             IF( ln_pnd_LEV ) THEN 
     1071            IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    10721072               sxap = 0._wp ;   syap = 0._wp    ;   sxxap = 0._wp    ;   syyap = 0._wp    ;   sxyap = 0._wp       ! melt pond fraction 
    10731073               sxvp = 0._wp ;   syvp = 0._wp    ;   sxxvp = 0._wp    ;   syyvp = 0._wp    ;   sxyvp = 0._wp       ! melt pond volume 
    … …  
    11371137         END DO 
    11381138         ! 
    1139          IF( ln_pnd_LEV ) THEN                                       ! melt pond fraction 
     1139         IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN                                       ! melt pond fraction 
    11401140            CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap  ) 
    11411141            CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap  ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icedyn_adv_umx.F90

    r14075 r15658  
    341341         ! 
    342342         !== melt ponds ==! 
    343          IF ( ln_pnd_LEV ) THEN 
     343         IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    344344            ! concentration 
    345345            zamsk = 1._wp 
    … …  
    361361         ! 
    362362         ! --- Lateral boundary conditions --- ! 
    363          IF    ( ln_pnd_LEV .AND. ln_pnd_lids ) THEN 
     363         IF    ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 
    364364            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
    365365               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 
    366          ELSEIF( ln_pnd_LEV .AND. .NOT.ln_pnd_lids ) THEN 
     366         ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 
    367367            CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 
    368368               &                                , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) 
    … …  
    16111611                  !                               ! -- check h_ip -- ! 
    16121612                  ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 
    1613                   IF( ln_pnd_LEV .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
     1613                  IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 
    16141614                     zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 
    16151615                     IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icedyn_rdgrft.F90

    r14075 r15658  
    567567               oirft2(ji) = oa_i_2d(ji,jl1)   * afrft * hi_hrft  
    568568 
    569                IF ( ln_pnd_LEV ) THEN 
     569               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    570570                  aprdg1     = a_ip_2d(ji,jl1) * afrdg 
    571571                  aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) 
    … …  
    604604               sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1    - sirft(ji) 
    605605               oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1    - oirft1 
    606                IF ( ln_pnd_LEV ) THEN 
     606               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    607607                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1    - aprft1 
    608608                  v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) 
    … …  
    701701                  v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji)  +  & 
    702702                     &                                  vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 
    703                   IF ( ln_pnd_LEV ) THEN 
     703                  IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    704704                     v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + (   vprdg (ji) * rn_fpndrdg * fvol   (ji)   & 
    705705                        &                                   + vprft (ji) * rn_fpndrft * zswitch(ji)   ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/iceitd.F90

    r14075 r15658  
    305305            IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 
    306306               a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin  
    307                IF( ln_pnd_LEV )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
     307               IF( ln_pnd_LEV .OR. ln_pnd_TOPO )   a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 
    308308               h_i_1d(ji) = rn_himin 
    309309            ENDIF 
    … …  
    476476               zaTsfn(ji,jl2)  = zaTsfn(ji,jl2) + ztrans 
    477477               !   
    478                IF ( ln_pnd_LEV ) THEN 
     478               IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    479479                  ztrans          = a_ip_2d(ji,jl1) * zworka(ji)     ! Pond fraction 
    480480                  a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icerst.F90

    r14075 r15658  
    2727   USE in_out_manager ! I/O manager 
    2828   USE iom            ! I/O manager library 
     29   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2930   USE lib_mpp        ! MPP library 
    3031   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    … …  
    5253      INTEGER, INTENT(in) ::   kt       ! number of iteration 
    5354      ! 
     55      INTEGER             ::   iyear, imonth, iday 
     56      REAL (wp)           ::   zsec 
     57      REAL (wp)           ::   zfjulday 
    5458      CHARACTER(len=20)   ::   clkt     ! ocean time-step define as a character 
    5559      CHARACTER(len=50)   ::   clname   ! ice output restart file name 
    … …  
    6771         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN 
    6872            ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    69             IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    70             ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst 
     73            IF ( ln_rstdate ) THEN 
     74               zfjulday = fjulday + (2*nn_fsbc+1)*rdt / rday 
     75               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     76               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     77               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     78            ELSE 
     79               IF( nitrst > 99999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     80               ELSE                           ;   WRITE(clkt, '(i8.8)') nitrst 
     81               ENDIF 
    7182            ENDIF 
    7283            ! create the file 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icestp.F90

    r14075 r15658  
    211211      ! --- Ocean time step --- ! 
    212212      !-------------------------! 
    213       IF( ln_icedyn )                   CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses 
     213      CALL ice_update_tau( kt, ub(:,:,1), vb(:,:,1) )   ! -- update surface ocean stresses 
    214214!!gm   remark, the ocean-ice stress is not saved in ice diag call above .....  find a solution!!! 
    215215      ! 
    … …  
    269269      ENDIF 
    270270      CALL ice_var_glo2eqv 
    271       CALL ice_var_agg(1) 
     271      CALL ice_var_agg(2) 
    272272      ! 
    273273      CALL ice_dyn_init                ! set ice dynamics parameters 
    … …  
    414414            wfx_pnd(ji,jj) = 0._wp 
    415415 
     416            dh_i_sum_2d(:,:,:) = 0._wp 
     417            dh_s_mlt_2d(:,:,:) = 0._wp 
     418 
    416419            hfx_thd(ji,jj) = 0._wp   ; 
    417420            hfx_snw(ji,jj) = 0._wp   ;   hfx_opw(ji,jj) = 0._wp 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icetab.F90

    r14075 r15658  
    4040      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind  ! input index 
    4141      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(in   ) ::   tab2d    ! input 2D field 
    42       REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(  out) ::   tab1d    ! output 1D field 
     42      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(inout) ::   tab1d    ! output 1D field 
    4343      ! 
    4444      INTEGER ::   jl, jn, jid, jjd 
    … …  
    6161      INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind  ! input index 
    6262      REAL(wp), DIMENSION(jpi,jpj), INTENT(in   ) ::   tab2d    ! input 2D field 
    63       REAL(wp), DIMENSION(ndim1d) , INTENT(  out) ::   tab1d    ! output 1D field 
     63      REAL(wp), DIMENSION(ndim1d) , INTENT(inout) ::   tab1d    ! output 1D field 
    6464      ! 
    6565      INTEGER ::   jn , jid, jjd 
    … …  
    8080      INTEGER , DIMENSION(ndim1d)     , INTENT(in   ) ::   tab_ind   ! input index 
    8181      REAL(wp), DIMENSION(ndim1d,jpl) , INTENT(in   ) ::   tab1d     ! input 1D field 
    82       REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(  out) ::   tab2d     ! output 2D field 
     82      REAL(wp), DIMENSION(jpi,jpj,jpl), INTENT(inout) ::   tab2d     ! output 2D field 
    8383      ! 
    8484      INTEGER ::   jl, jn, jid, jjd 
    … …  
    101101      INTEGER , DIMENSION(ndim1d) , INTENT(in   ) ::   tab_ind   ! input index 
    102102      REAL(wp), DIMENSION(ndim1d) , INTENT(in   ) ::   tab1d     ! input 1D field 
    103       REAL(wp), DIMENSION(jpi,jpj), INTENT(  out) ::   tab2d     ! output 2D field 
     103      REAL(wp), DIMENSION(jpi,jpj), INTENT(inout) ::   tab2d     ! output 2D field 
    104104      ! 
    105105      INTEGER ::   jn , jid, jjd 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icethd.F90

    r14075 r15658  
    247247            IF( ln_icedH ) THEN                                     ! --- Growing/Melting --- ! 
    248248                              CALL ice_thd_dh                           ! Ice-Snow thickness    
    249                               CALL ice_thd_pnd                          ! Melt ponds formation 
    250249                              CALL ice_thd_ent( e_i_1d(1:npti,:) )      ! Ice enthalpy remapping 
    251250            ENDIF 
    … …  
    268267      IF( ln_icediachk )   CALL ice_cons2D  (1, 'icethd',  diag_v,  diag_s,  diag_t,  diag_fv,  diag_fs,  diag_ft) 
    269268      !                    
     269      IF ( ln_pnd  )          CALL ice_thd_pnd                      ! --- Melt ponds 
     270 
    270271      IF( jpl > 1  )          CALL ice_itd_rem( kt )                ! --- Transport ice between thickness categories --- ! 
    271272      ! 
    … …  
    536537         CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) ) 
    537538         CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) ) 
     539         ! ponds 
     540         CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum  (1:npti) , dh_i_sum_2d(:,:,kl) ) 
     541         CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt  (1:npti) , dh_s_mlt_2d(:,:,kl) ) 
    538542         ! SIMIP diagnostics          
    539543         CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d       (1:npti), t_si       (:,:,kl) ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icethd_pnd.F90

    r14075 r15658  
    2020   USE ice1D          ! sea-ice: thermodynamics variables 
    2121   USE icetab         ! sea-ice: 1D <==> 2D transformation 
     22   USE sbc_ice        ! surface energy budget 
    2223   ! 
    2324   USE in_out_manager ! I/O manager 
     25   USE iom            ! I/O manager library 
    2426   USE lib_mpp        ! MPP library 
    2527   USE lib_fortran    ! fortran utilities (glob_sum + no signed zero) 
    … …  
    3436   INTEGER ::              nice_pnd    ! choice of the type of pond scheme 
    3537   !                                   ! associated indices: 
    36    INTEGER, PARAMETER ::   np_pndNO  = 0   ! No pond scheme 
    37    INTEGER, PARAMETER ::   np_pndCST = 1   ! Constant ice pond scheme 
    38    INTEGER, PARAMETER ::   np_pndLEV = 2   ! Level ice pond scheme 
     38   INTEGER, PARAMETER ::   np_pndNO   = 0   ! No pond scheme 
     39   INTEGER, PARAMETER ::   np_pndCST  = 1   ! Constant ice pond scheme 
     40   INTEGER, PARAMETER ::   np_pndLEV  = 2   ! Level ice pond scheme 
     41   INTEGER, PARAMETER ::   np_pndTOPO = 3   ! Level ice pond scheme 
     42 
     43   REAL(wp), PARAMETER :: &   ! shared parameters for topographic melt ponds 
     44      zhi_min = 0.1_wp        , & ! minimum ice thickness with ponds (m)  
     45      zTd     = 273.0_wp      , & ! temperature difference for freeze-up (C) 
     46      zvp_min = 1.e-4_wp          ! minimum pond volume (m) 
     47 
     48   !-------------------------------------------------------------------------- 
     49   !  
     50   ! Pond volume per area budget diags 
     51   !   
     52   ! The idea of diags is the volume of ponds per grid cell area is 
     53   ! 
     54   ! dV/dt = mlt + drn + lid + rnf 
     55   ! mlt   = input from surface melting 
     56   ! drn   = drainage through brine network 
     57   ! lid   = lid growth & melt 
     58   ! rnf   = runoff (water directly removed out of surface melting + overflow) 
     59   ! 
     60   ! In topo mode, the pond water lost because it is in the snow is not included in the budget 
     61   ! 
     62   ! In level mode, all terms are incorporated 
     63 
     64   REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics 
     65      diag_dvpn_mlt,  &     !! meltwater pond volume input      [m/day] 
     66      diag_dvpn_drn,  &     !! pond volume lost by drainage     [m/day] 
     67      diag_dvpn_lid,  &     !! exchange with lid / refreezing   [m/day] 
     68      diag_dvpn_rnf         !! meltwater pond lost to runoff    [m/day] 
     69       
     70   REAL(wp), ALLOCATABLE, DIMENSION(:) ::   & ! pond volume budget diagnostics (1d) 
     71      diag_dvpn_mlt_1d,  &  !! meltwater pond volume input      [m/day] 
     72      diag_dvpn_drn_1d,  &  !! pond volume lost by drainage     [m/day] 
     73      diag_dvpn_lid_1d,  &  !! exchange with lid / refreezing   [m/day] 
     74      diag_dvpn_rnf_1d      !! meltwater pond lost to runoff    [m/day] 
     75   ! 
     76   !-------------------------------------------------------------------------- 
    3977 
    4078   !! * Substitutions 
    … …  
    5290      !!                
    5391      !! ** Purpose :   change melt pond fraction and thickness 
    54       !!                 
     92      !! 
     93      !! Note: Melt ponds affect only radiative transfer for now 
     94      !! 
     95      !!       No heat, no salt. 
     96      !!       The melt water they carry is collected but  
     97      !!       not removed from fw budget or released to the ocean 
     98      !! 
     99      !!       A wfx_pnd has been coded for diagnostic purposes 
     100      !!       It is not fully consistent yet. 
     101      !! 
     102      !!       The current diagnostic lacks a contribution from drainage 
     103      !! 
    55104      !!------------------------------------------------------------------- 
    56       ! 
     105      !! 
     106       
     107      ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 
     108      ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 
     109 
     110      diag_dvpn_mlt(:,:) = 0._wp    ;   diag_dvpn_drn(:,:)  = 0._wp 
     111      diag_dvpn_lid(:,:) = 0._wp    ;   diag_dvpn_rnf(:,:)  = 0._wp 
     112      diag_dvpn_mlt_1d(:) = 0._wp   ;   diag_dvpn_drn_1d(:) = 0._wp 
     113      diag_dvpn_lid_1d(:) = 0._wp   ;   diag_dvpn_rnf_1d(:) = 0._wp 
     114 
    57115      SELECT CASE ( nice_pnd ) 
    58116      ! 
    59       CASE (np_pndCST)   ;   CALL pnd_CST    !==  Constant melt ponds  ==! 
     117      CASE (np_pndCST)   ;   CALL pnd_CST                              !==  Constant melt ponds  ==! 
    60118         ! 
    61       CASE (np_pndLEV)   ;   CALL pnd_LEV    !==  Level ice melt ponds  ==! 
     119      CASE (np_pndLEV)   ;   CALL pnd_LEV                              !==  Level ice melt ponds  ==! 
     120         ! 
     121      CASE (np_pndTOPO)  ;   CALL pnd_TOPO                             !==  Topographic melt ponds  ==! 
    62122         ! 
    63123      END SELECT 
    64124      ! 
     125 
     126      IF( iom_use('dvpn_mlt'  ) )   CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 
     127      IF( iom_use('dvpn_lid'  ) )   CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 
     128      IF( iom_use('dvpn_drn'  ) )   CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 
     129      IF( iom_use('dvpn_rnf'  ) )   CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 
     130 
     131      DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 
     132      DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 
     133       
    65134   END SUBROUTINE ice_thd_pnd  
    66135 
    … …  
    75144      !!                to non-zero values when t_su = 0C 
    76145      !! 
    77       !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd) 
     146      !! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd) 
    78147      !!                 
    79148      !! ** Note   : Coupling with such melt ponds is only radiative 
    … …  
    145214      !!                                  a_ip/a_i = a_ip_frac = h_ip / zaspect 
    146215      !! 
    147       !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 
     216      !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 
    148217      !!  
    149       !! ** Note       :   mostly stolen from CICE 
    150       !! 
    151       !! ** References :   Flocco and Feltham (JGR, 2007) 
    152       !!                   Flocco et al       (JGR, 2010) 
    153       !!                   Holland et al      (J. Clim, 2012) 
     218      !! ** Note       :   Mostly stolen from CICE but not only 
     219      !! 
     220      !! ** References :   Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 
     221      !!                    
    154222      !!------------------------------------------------------------------- 
    155223      REAL(wp), DIMENSION(nlay_i) ::   ztmp           ! temporary array 
    … …  
    168236      REAL(wp) ::   zfac, zdum                        ! temporary arrays 
    169237      REAL(wp) ::   z1_rhow, z1_aspect, z1_Tp         ! inverse 
    170       !! 
    171       INTEGER  ::   ji, jk                            ! loop indices 
     238      REAL(wp) ::   zvold                             ! 
     239      !! 
     240      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     241       
    172242      !!------------------------------------------------------------------- 
    173243      z1_rhow   = 1._wp / rhow  
    174244      z1_aspect = 1._wp / zaspect 
    175245      z1_Tp     = 1._wp / zTp  
    176  
    177       DO ji = 1, npti 
    178          !                                                            !----------------------------------------------------! 
    179          IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
    180             !                                                         !----------------------------------------------------! 
    181             !--- Remove ponds on thin ice or tiny ice fractions 
    182             a_ip_1d(ji)      = 0._wp 
    183             h_ip_1d(ji)      = 0._wp 
    184             h_il_1d(ji)      = 0._wp 
    185             !                                                         !--------------------------------! 
    186          ELSE                                                         ! Case ice thickness >= rn_himin ! 
    187             !                                                         !--------------------------------! 
    188             v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
    189             v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
    190             ! 
    191             !------------------! 
    192             ! case ice melting ! 
    193             !------------------! 
    194             ! 
    195             !--- available meltwater for melt ponding ---! 
    196             zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
    197             zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
    198             zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
    199             ! 
    200             !--- overflow ---! 
    201             ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 
    202             !    a_ip_max = zfr_mlt * a_i 
    203             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    204             zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
    205             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
    206  
    207             ! If pond depth exceeds half the ice thickness then reduce the pond volume 
    208             !    h_ip_max = 0.5 * h_i 
    209             !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
    210             zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 
    211             zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     246       
     247      !----------------------------------------------------------------------------------------------- 
     248      !  Identify grid cells with ice 
     249      !----------------------------------------------------------------------------------------------- 
     250      at_i(:,:) = SUM( a_i, dim=3 ) 
     251      ! 
     252      npti = 0   ;   nptidx(:) = 0 
     253      DO jj = 1, jpj 
     254         DO ji = 1, jpi 
     255            IF ( at_i(ji,jj) > epsi10 ) THEN 
     256               npti = npti + 1 
     257               nptidx( npti ) = (jj - 1) * jpi + ji 
     258            ENDIF 
     259         END DO 
     260      END DO 
     261       
     262      !----------------------------------------------------------------------------------------------- 
     263      ! Prepare 1D arrays 
     264      !----------------------------------------------------------------------------------------------- 
     265 
     266      IF( npti > 0 ) THEN 
     267       
     268         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum      ) 
     269         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti)   , wfx_sum          ) 
     270         CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti)   , wfx_pnd          ) 
     271          
     272         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 
     273         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 
     274         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 
     275         CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 
     276       
     277         DO jl = 1, jpl 
     278          
     279            CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d      (1:npti), a_i      (:,:,jl) ) 
     280            CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d      (1:npti), h_i      (:,:,jl) ) 
     281            CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d     (1:npti), t_su     (:,:,jl) ) 
     282            CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
     283            CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     284            CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
     285 
     286            CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum    (1:npti), dh_i_sum_2d     (:,:,jl) ) 
     287            CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt    (1:npti), dh_s_mlt_2d     (:,:,jl) ) 
    212288             
    213             !--- Pond growing ---! 
    214             v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
    215             ! 
    216             !--- Lid melting ---! 
    217             IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
    218             ! 
    219             !--- mass flux ---! 
    220             IF( zdv_mlt > 0._wp ) THEN 
    221                zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
    222                wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
    223                ! 
    224                zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
    225                wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
    226                wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
    227             ENDIF 
    228  
    229             !-------------------! 
    230             ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
    231             !-------------------! 
    232             ! 
    233             zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
    234             ! 
    235             !--- Pond contraction (due to refreezing) ---! 
    236             IF( ln_pnd_lids ) THEN 
    237                ! 
    238                !--- Lid growing and subsequent pond shrinking ---!  
    239                zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
    240                   &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
    241                 
    242                ! Lid growing 
    243                v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
    244                 
    245                ! Pond shrinking 
    246                v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
    247  
    248             ELSE 
    249                ! Pond shrinking 
    250                v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
    251             ENDIF 
    252             ! 
    253             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    254             ! v_ip     = h_ip * a_ip 
    255             ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 
    256             a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    257             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    258  
    259             !---------------!             
    260             ! Pond flushing ! 
    261             !---------------! 
    262             ! height of top of the pond above sea-level 
    263             zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 
     289            DO jk = 1, nlay_i 
     290               CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
     291            END DO 
    264292             
    265             ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
    266             DO jk = 1, nlay_i 
    267                zsbr = - 1.2_wp                                  & 
    268                   &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
    269                   &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
    270                   &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
    271                ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
    272             END DO 
    273             zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
    274              
    275             ! Do the drainage using Darcy's law 
    276             zdv_flush   = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 
    277             zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
    278             v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
    279              
    280             !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
    281             a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
    282             h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
    283  
    284             !--- Corrections and lid thickness ---! 
    285             IF( ln_pnd_lids ) THEN 
    286                !--- retrieve lid thickness from volume ---! 
    287                IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
    288                ELSE                              ;   h_il_1d(ji) = 0._wp 
    289                ENDIF 
    290                !--- remove ponds if lids are much larger than ponds ---! 
    291                IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     293            !----------------------------------------------------------------------------------------------- 
     294            ! Go for ponds 
     295            !----------------------------------------------------------------------------------------------- 
     296 
     297            DO ji = 1, npti 
     298               !                                                            !----------------------------------------------------! 
     299               IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN    ! Case ice thickness < rn_himin or tiny ice fraction ! 
     300                  !                                                         !----------------------------------------------------! 
     301                  !--- Remove ponds on thin ice or tiny ice fractions 
    292302                  a_ip_1d(ji)      = 0._wp 
    293303                  h_ip_1d(ji)      = 0._wp 
    294304                  h_il_1d(ji)      = 0._wp 
     305                  !                                                         !--------------------------------! 
     306               ELSE                                                         ! Case ice thickness >= rn_himin ! 
     307                  !                                                         !--------------------------------! 
     308                  v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)   ! retrieve volume from thickness 
     309                  v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 
     310                  ! 
     311                  !------------------! 
     312                  ! case ice melting ! 
     313                  !------------------! 
     314                  ! 
     315                  !--- available meltwater for melt ponding ---! 
     316                  zdum    = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 
     317                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) !  = ( 1 - r ) = fraction of melt water that is not flushed 
     318                  zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?  
     319                   
     320                  diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_rdtice                      ! surface melt input diag 
     321                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_rdtice   ! runoff diag 
     322                  ! 
     323                  !--- overflow ---! 
     324                  ! 
     325                  ! 1) area driven overflow 
     326                  ! 
     327                  ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 
     328                  !    a_ip_max = zfr_mlt * a_i 
     329                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     330                  zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 
     331                  zvold     = zdv_mlt 
     332                  zdv_mlt   = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )                      
     333                   
     334                  !  
     335                  ! 2) depth driven overflow 
     336                  ! 
     337                  ! If pond depth exceeds half the ice thickness then reduce the pond volume 
     338                  !    h_ip_max = 0.5 * h_i 
     339                  !    => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:  
     340                  zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 
     341                  zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 
     342                  diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_rdtice       ! runoff diag - overflow contribution 
     343             
     344                  !--- Pond growing ---! 
     345                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 
     346                  ! 
     347                  !--- Lid melting ---! 
     348                  IF( ln_pnd_lids )   v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 
     349                  ! 
     350                  !--- mass flux ---! 
     351                  ! MV I would recommend to remove that 
     352                  ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 
     353                   
     354                  IF( zdv_mlt > 0._wp ) THEN 
     355                     zfac = zdv_mlt * rhow * r1_rdtice                        ! melt pond mass flux < 0 [kg.m-2.s-1] 
     356                     wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 
     357                     ! 
     358                     zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) )    ! adjust ice/snow melting flux > 0 to balance melt pond flux 
     359                     wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 
     360                     wfx_sum_1d(ji)     = wfx_sum_1d(ji)     * (1._wp + zdum) 
     361                  ENDIF 
     362 
     363                  !-------------------! 
     364                  ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 
     365                  !-------------------! 
     366                  ! 
     367                  zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 
     368                  ! 
     369                  !--- Pond contraction (due to refreezing) ---! 
     370                  zvold       = v_ip_1d(ji) ! for diag 
     371 
     372                  IF( ln_pnd_lids ) THEN 
     373                     ! 
     374                     !--- Lid growing and subsequent pond shrinking ---!  
     375                     zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 
     376                        &                    SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 
     377                
     378                     ! Lid growing 
     379                     v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 
     380                
     381                     ! Pond shrinking 
     382                     v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 
     383                      
     384                  ELSE 
     385                   
     386                     ! Pond shrinking 
     387                     v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 
     388                      
     389                  ENDIF 
     390 
     391                  diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_rdtice   ! shrinking counted as lid diagnostic 
     392 
     393                  ! 
     394                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     395                  ! v_ip     = h_ip * a_ip 
     396                  ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 
     397                  a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
     398                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     399 
     400                  !------------------------------------------------!             
     401                  ! Pond drainage through brine network (flushing) ! 
     402                  !------------------------------------------------! 
     403                  ! height of top of the pond above sea-level 
     404                  zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 
     405             
     406                  ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 
     407                  DO jk = 1, nlay_i 
     408                     ! MV Assur is inconsistent with SI3 
     409                     zsbr = - 1.2_wp                                  & 
     410                        &   - 21.8_wp    * ( t_i_1d(ji,jk) - rt0 )    & 
     411                        &   - 0.919_wp   * ( t_i_1d(ji,jk) - rt0 )**2 & 
     412                        &   - 0.0178_wp  * ( t_i_1d(ji,jk) - rt0 )**3 
     413                     ! MV linear expression more consistent & simpler              zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 
     414                     ztmp(jk) = sz_i_1d(ji,jk) / zsbr 
     415                  END DO 
     416                  zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 
     417             
     418                  ! Do the drainage using Darcy's law 
     419                  zdv_flush   = - zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 
     420                  zdv_flush   = MAX( zdv_flush, -v_ip_1d(ji) ) 
     421                  ! zdv_flush   = 0._wp ! MV remove pond drainage for now 
     422                  v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 
     423                   
     424                  diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_rdtice   ! shrinking counted as lid diagnostic 
     425             
     426                  ! MV --- why pond drainage does not give back water into freshwater flux ?             
     427                  !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 
     428                  a_ip_1d(ji)      = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 
     429                  h_ip_1d(ji)      = zaspect * a_ip_1d(ji) / a_i_1d(ji) 
     430 
     431                  !--- Corrections and lid thickness ---! 
     432                  IF( ln_pnd_lids ) THEN 
     433                     !--- retrieve lid thickness from volume ---! 
     434                     IF( a_ip_1d(ji) > epsi10 ) THEN   ;   h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 
     435                     ELSE                              ;   h_il_1d(ji) = 0._wp 
     436                     ENDIF 
     437                     !--- remove ponds if lids are much larger than ponds ---! 
     438                     IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 
     439                        a_ip_1d(ji)      = 0._wp 
     440                        h_ip_1d(ji)      = 0._wp 
     441                        h_il_1d(ji)      = 0._wp 
     442                     ENDIF 
     443                  ENDIF 
     444                  ! 
    295445               ENDIF 
     446                
     447            END DO ! ji 
     448 
     449            !----------------------------------------------------------------------------------------------- 
     450            ! Retrieve 2D arrays 
     451            !----------------------------------------------------------------------------------------------- 
     452             
     453            v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 
     454            v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 
     455             
     456            CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d     (1:npti), a_ip     (:,:,jl) ) 
     457            CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d     (1:npti), h_ip     (:,:,jl) ) 
     458            CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d     (1:npti), h_il     (:,:,jl) ) 
     459            CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d     (1:npti), v_ip     (:,:,jl) ) 
     460            CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d     (1:npti), v_il     (:,:,jl) ) 
     461            DO jk = 1, nlay_i 
     462               CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl)  ) 
     463            END DO 
     464    
     465         END DO ! ji 
     466                   
     467         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 
     468         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum        ) 
     469         CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd        ) 
     470          
     471         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt        ) 
     472         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn        ) 
     473         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid        ) 
     474         CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf        )                   
     475                            
     476      ! 
     477      ENDIF 
     478       
     479   END SUBROUTINE pnd_LEV 
     480    
     481   SUBROUTINE pnd_TOPO     
     482                                          
     483      !!------------------------------------------------------------------- 
     484      !!                ***  ROUTINE pnd_TOPO  *** 
     485      !! 
     486      !! ** Purpose :   Compute melt pond evolution based on the ice 
     487      !!                topography inferred from the ice thickness distribution 
     488      !! 
     489      !! ** Method  :   This code is initially based on Flocco and Feltham 
     490      !!                (2007) and Flocco et al. (2010).  
     491      !! 
     492      !!                - Calculate available pond water base on surface meltwater 
     493      !!                - Redistribute water as a function of topography, drain water 
     494      !!                - Exchange water with the lid 
     495      !! 
     496      !! ** Tunable parameters : 
     497      !! 
     498      !! ** Note : 
     499      !! 
     500      !! ** References 
     501      !! 
     502      !!    Flocco, D. and D. L. Feltham, 2007.  A continuum model of melt pond 
     503      !!    evolution on Arctic sea ice.  J. Geophys. Res. 112, C08016, doi: 
     504      !!    10.1029/2006JC003836. 
     505      !! 
     506      !!    Flocco, D., D. L. Feltham and A. K. Turner, 2010.  Incorporation of 
     507      !!    a physically based melt pond scheme into the sea ice component of a 
     508      !!    climate model.  J. Geophys. Res. 115, C08012, 
     509      !!    doi: 10.1029/2009JC005568. 
     510      !! 
     511      !!------------------------------------------------------------------- 
     512 
     513      ! local variables 
     514      REAL(wp) :: & 
     515         zdHui,   &      ! change in thickness of ice lid (m) 
     516         zomega,  &      ! conduction 
     517         zdTice,  &      ! temperature difference across ice lid (C) 
     518         zdvice,  &      ! change in ice volume (m) 
     519         zTavg,   &      ! mean surface temperature across categories (C) 
     520         zfsurf,  &      ! net heat flux, excluding conduction and transmitted radiation (W/m2) 
     521         zTp,     &      ! pond freezing temperature (C) 
     522         zrhoi_L, &      ! volumetric latent heat of sea ice (J/m^3) 
     523         zfr_mlt, &      ! fraction and volume of available meltwater retained for melt ponding 
     524         z1_rhow, &      ! inverse water density 
     525         zv_pnd  , &     ! volume of meltwater contributing to ponds 
     526         zv_mlt          ! total amount of meltwater produced 
     527 
     528      REAL(wp), DIMENSION(jpi,jpj) ::   zvolp_ini, & !! total melt pond water available before redistribution and drainage 
     529                                        zvolp    , & !! total melt pond water volume 
     530                                        zvolp_res    !! remaining melt pond water available after drainage 
     531                                         
     532      REAL(wp), DIMENSION(jpi,jpj,jpl) ::   z1_a_i 
     533 
     534      INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     535 
     536      INTEGER  ::   i_test 
     537 
     538      ! Note 
     539      ! equivalent for CICE translation 
     540      ! a_ip      -> apond 
     541      ! a_ip_frac -> apnd 
     542       
     543      !--------------------------------------------------------------- 
     544      ! Initialise 
     545      !--------------------------------------------------------------- 
     546 
     547      ! Parameters & constants (move to parameters) 
     548      zrhoi_L   = rhoi * rLfus      ! volumetric latent heat (J/m^3) 
     549      zTp       = rt0 - 0.15_wp          ! pond freezing point, slightly below 0C (ponds are bid saline) 
     550      z1_rhow   = 1._wp / rhow  
     551 
     552      ! Set required ice variables (hard-coded here for now) 
     553!      zfpond(:,:) = 0._wp          ! contributing freshwater flux (?)  
     554       
     555      at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 
     556      vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 
     557      vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 
     558      vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 
     559       
     560      ! thickness 
     561      WHERE( a_i(:,:,:) > epsi20 )   ;   z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 
     562      ELSEWHERE                      ;   z1_a_i(:,:,:) = 0._wp 
     563      END WHERE 
     564      h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 
     565       
     566      !--------------------------------------------------------------- 
     567      ! Change 2D to 1D 
     568      !--------------------------------------------------------------- 
     569      ! MV  
     570      ! a less computing-intensive version would have 2D-1D passage here 
     571      ! use what we have in iceitd.F90 (incremental remapping) 
     572 
     573      !-------------------------------------------------------------- 
     574      ! Collect total available pond water volume 
     575      !-------------------------------------------------------------- 
     576      ! Assuming that meltwater (+rain in principle) runsoff the surface 
     577      ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 
     578      ! I cite her words, they are very talkative 
     579      ! "grid cells with very little ice cover (and hence more open water area)  
     580      ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 
     581      ! "This results in the same runoff fraction r for each ice category within a grid cell" 
     582       
     583      zvolp(:,:) = 0._wp 
     584      zvolp_res(:,:) = 0._wp 
     585 
     586      DO jl = 1, jpl 
     587         DO jj = 1, jpj 
     588            DO ji = 1, jpi 
     589                  
     590               IF ( a_i(ji,jj,jl) > epsi10 ) THEN 
     591             
     592                  !--- Available and contributing meltwater for melt ponding ---! 
     593                  zv_mlt  = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) &        ! available volume of surface melt water per grid area 
     594                     &    * z1_rhow * a_i(ji,jj,jl) 
     595                      ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 
     596                  zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj)                  ! fraction of surface meltwater going to ponds 
     597                  zv_pnd  = zfr_mlt * zv_mlt                                                           ! contributing meltwater volume for category jl 
     598 
     599                  diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice                     ! diags 
     600                  diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_rdtice    
     601 
     602                  !--- Create possible new ponds 
     603                  ! if pond does not exist, create new pond over full ice area 
     604                  !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 
     605                  IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 
     606                     a_ip(ji,jj,jl)       = a_i(ji,jj,jl) 
     607                     a_ip_frac(ji,jj,jl)  = 1.0_wp    ! pond fraction of sea ice (apnd for CICE) 
     608                  ENDIF 
     609                   
     610                  !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 
     611                  v_ip(ji,jj,jl) = v_ip(ji,jj,jl) +  zv_pnd                                            ! use pond water to increase thickness 
     612                  h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     613                   
     614                  !--- Total available pond water volume (pre-existing + newly produced)j 
     615                  zvolp(ji,jj)   = zvolp(ji,jj)   + v_ip(ji,jj,jl)  
     616!                 zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 
     617                    
     618               ENDIF ! a_i 
     619 
     620            END DO! jl             
     621         END DO ! jj 
     622      END DO ! ji 
     623 
     624      zvolp_ini(:,:) = zvolp(:,:) 
     625                   
     626      !-------------------------------------------------------------- 
     627      ! Redistribute and drain water from ponds 
     628      !--------------------------------------------------------------    
     629      CALL ice_thd_pnd_area( zvolp, zvolp_res ) 
     630                                    
     631      !-------------------------------------------------------------- 
     632      ! Melt pond lid growth and melt 
     633      !--------------------------------------------------------------    
     634       
     635      IF( ln_pnd_lids ) THEN 
     636 
     637         DO jj = 1, jpj 
     638            DO ji = 1, jpi 
     639 
     640            IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp_ini(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     641                   
     642               !-------------------------- 
     643               ! Pond lid growth and melt 
     644               !-------------------------- 
     645               ! Mean surface temperature 
     646               zTavg = 0._wp 
     647               DO jl = 1, jpl 
     648                  zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 
     649               END DO 
     650               zTavg = zTavg / at_i(ji,jj) !!! could get a division by zero here 
     651          
     652               DO jl = 1, jpl-1 
     653             
     654                  IF ( v_il(ji,jj,jl) > epsi10 ) THEN 
     655                
     656                     !---------------------------------------------------------------- 
     657                     ! Lid melting: floating upper ice layer melts in whole or part 
     658                     !---------------------------------------------------------------- 
     659                     ! Use Tsfc for each category 
     660                     IF ( t_su(ji,jj,jl) > zTp ) THEN 
     661 
     662                        zdvice = MIN( -dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 
     663                         
     664                        IF ( zdvice > epsi10 ) THEN 
     665                         
     666                           v_il (ji,jj,jl) = v_il (ji,jj,jl)   - zdvice 
     667                           v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)    + zdvice ! MV: not sure i understand dh_i_sum seems counted twice -  
     668                                                                        ! as it is already counted in surface melt 
     669!                          zvolp(ji,jj)     = zvolp(ji,jj)     + zdvice ! pointless to calculate total volume (done in icevar) 
     670!                          zfpond(ji,jj)    = fpond(ji,jj)     + zdvice ! pointless to follow fw budget (ponds have no fw) 
     671                      
     672                           IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 
     673                           ! ice lid melted and category is pond covered 
     674                              v_ip(ji,jj,jl)  = v_ip(ji,jj,jl)  + v_il(ji,jj,jl)  
     675!                             zfpond(ji,jj)    = zfpond (ji,jj)    + v_il(ji,jj,jl)  
     676                              v_il(ji,jj,jl)   = 0._wp 
     677                           ENDIF 
     678                           h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 
     679                            
     680                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice   ! diag 
     681                            
     682                        ENDIF 
     683                         
     684                     !---------------------------------------------------------------- 
     685                     ! Freeze pre-existing lid  
     686                     !---------------------------------------------------------------- 
     687 
     688                     ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 
     689 
     690                        ! differential growth of base of surface floating ice layer 
     691                        ! zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0  ! line valid for temp in C  
     692                        zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0    
     693                        zomega = rcnd_i * zdTice / zrhoi_L 
     694                        zdHui  = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 
     695                               - v_il(ji,jj,jl) / a_i(ji,jj,jl) 
     696                        zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     697                   
     698                        IF ( zdvice > epsi10 ) THEN 
     699                           v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     700                           v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     701!                          zvolp(ji,jj)    = zvolp(ji,jj)     - zdvice 
     702!                          zfpond(ji,jj)   = zfpond(ji,jj)    - zdvice 
     703                           h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     704                            
     705                           diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice    ! diag 
     706                            
     707                        ENDIF 
     708                   
     709                     ENDIF ! Tsfcn(i,j,n) 
     710 
     711                  !---------------------------------------------------------------- 
     712                  ! Freeze new lids 
     713                  !---------------------------------------------------------------- 
     714                  !  upper ice layer begins to form 
     715                  ! note: albedo does not change 
     716 
     717                  ELSE ! v_il < epsi10 
     718                     
     719                     ! thickness of newly formed ice 
     720                     ! the surface temperature of a meltpond is the same as that 
     721                     ! of the ice underneath (0C), and the thermodynamic surface  
     722                     ! flux is the same 
     723                      
     724                     !!! we need net surface energy flux, excluding conduction 
     725                     !!! fsurf is summed over categories in CICE 
     726                     !!! we have the category-dependent flux, let us use it ? 
     727                     zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl)                      
     728                     zdHui  = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp ) 
     729                     zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 
     730 
     731                     IF ( zdvice > epsi10 ) THEN 
     732 
     733                        v_il (ji,jj,jl)  = v_il(ji,jj,jl)   + zdvice 
     734                        v_ip(ji,jj,jl)   = v_ip(ji,jj,jl)   - zdvice 
     735                        diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice      ! diag 
     736 
     737!                       zvolp(ji,jj)     = zvolp(ji,jj)     - zdvice 
     738!                       zfpond(ji,jj)    = zfpond(ji,jj)    - zdvice 
     739 
     740                        h_ip(ji,jj,jl)   = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar 
     741                     ENDIF 
     742                
     743                  ENDIF  ! v_il 
     744 
     745                  !---------------------------------------------------------------- 
     746                  ! Remove ice lid if there is no liquid pond 
     747                  !---------------------------------------------------------------- 
     748 
     749                  IF (v_ip(ji,jj,jl) < epsi10) THEN 
     750                     v_il(ji,jj,jl) = 0._wp 
     751                  ENDIF 
     752             
     753               END DO ! jl 
     754 
     755            ELSE  ! remove ponds on thin ice 
     756 
     757               v_ip(ji,jj,:) = 0._wp 
     758               v_il(ji,jj,:) = 0._wp 
     759!              zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 
     760!                 zvolp(ji,jj)    = 0._wp          
     761 
    296762            ENDIF 
    297             ! 
    298          ENDIF 
    299           
    300       END DO 
    301       ! 
    302    END SUBROUTINE pnd_LEV 
    303  
     763 
     764      END DO   ! jj 
     765   END DO   ! ji 
     766 
     767      ENDIF ! ln_pnd_lids 
     768 
     769      !--------------------------------------------------------------- 
     770      ! Clean-up variables (probably duplicates what icevar would do) 
     771      !--------------------------------------------------------------- 
     772      ! MV comment 
     773      ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 
     774      ! icevar should recompute all other variables (if needed at all) 
     775 
     776      DO jl = 1, jpl 
     777 
     778         DO jj = 1, jpj 
     779            DO ji = 1, jpi 
     780 
     781!              ! zap lids on small ponds 
     782!              IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 
     783!                                          .AND. v_il(ji,jj,jl) > epsi10) THEN 
     784!                 v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 
     785!              ENDIF 
     786       
     787               ! recalculate equivalent pond variables 
     788               IF ( a_ip(ji,jj,jl) > epsi10) THEN 
     789                  h_ip(ji,jj,jl)      = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     790                  a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     791                  h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 
     792               ENDIF 
     793!                 h_ip(ji,jj,jl)      = 0._wp ! MV in principle, useless as computed in icevar 
     794!                 h_il(ji,jj,jl)      = 0._wp ! MV in principle, useless as omputed in icevar 
     795!              ENDIF 
     796                
     797            END DO   ! ji 
     798         END DO   ! jj 
     799 
     800      END DO   ! jl 
     801 
     802   END SUBROUTINE pnd_TOPO 
     803 
     804    
     805   SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 
     806 
     807       !!------------------------------------------------------------------- 
     808       !!                ***  ROUTINE ice_thd_pnd_area *** 
     809       !! 
     810       !! ** Purpose : Given the total volume of available pond water,  
     811       !!              redistribute and drain water 
     812       !! 
     813       !! ** Method 
     814       !! 
     815       !-----------| 
     816       !           | 
     817       !           |-----------| 
     818       !___________|___________|______________________________________sea-level 
     819       !           |           | 
     820       !           |           |---^--------| 
     821       !           |           |   |        | 
     822       !           |           |   |        |-----------|              |------- 
     823       !           |           |   | alfan  |           |              | 
     824       !           |           |   |        |           |--------------| 
     825       !           |           |   |        |           |              | 
     826       !---------------------------v------------------------------------------- 
     827       !           |           |   ^        |           |              | 
     828       !           |           |   |        |           |--------------| 
     829       !           |           |   | betan  |           |              | 
     830       !           |           |   |        |-----------|              |------- 
     831       !           |           |   |        | 
     832       !           |           |---v------- | 
     833       !           |           | 
     834       !           |-----------| 
     835       !           | 
     836       !-----------| 
     837       ! 
     838       !! 
     839       !!------------------------------------------------------------------ 
     840        
     841       REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 
     842          zvolp,                                       &  ! total available pond water 
     843          zdvolp                                          ! remaining meltwater after redistribution 
     844 
     845       INTEGER ::  & 
     846          ns,      & 
     847          m_index, & 
     848          permflag 
     849 
     850       REAL (wp), DIMENSION(jpl) :: & 
     851          hicen, & 
     852          hsnon, & 
     853          asnon, & 
     854          alfan, & 
     855          betan, & 
     856          cum_max_vol, & 
     857          reduced_aicen 
     858 
     859       REAL (wp), DIMENSION(0:jpl) :: & 
     860          cum_max_vol_tmp 
     861 
     862       REAL (wp) :: & 
     863          hpond, & 
     864          drain, & 
     865          floe_weight, & 
     866          pressure_head, & 
     867          hsl_rel, & 
     868          deltah, & 
     869          perm, & 
     870          msno 
     871 
     872       REAL (wp), parameter :: & 
     873          viscosity = 1.79e-3_wp     ! kinematic water viscosity in kg/m/s 
     874 
     875       INTEGER  ::   ji, jj, jk, jl                    ! loop indices 
     876 
     877       a_ip(:,:,:) = 0._wp 
     878       h_ip(:,:,:) = 0._wp 
     879  
     880       DO jj = 1, jpj 
     881          DO ji = 1, jpi 
     882  
     883             IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 
     884  
     885        !------------------------------------------------------------------- 
     886        ! initialize 
     887        !------------------------------------------------------------------- 
     888  
     889        DO jl = 1, jpl 
     890  
     891           !---------------------------------------- 
     892           ! compute the effective snow fraction 
     893           !---------------------------------------- 
     894  
     895           IF (a_i(ji,jj,jl) < epsi10)  THEN 
     896              hicen(jl) =  0._wp 
     897              hsnon(jl) =  0._wp 
     898              reduced_aicen(jl) = 0._wp 
     899              asnon(jl) = 0._wp         !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 
     900           ELSE 
     901              hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 
     902              hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 
     903              reduced_aicen(jl) = 1._wp ! n=jpl 
     904  
     905              !js: initial code in NEMO_DEV 
     906              !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 
     907              !                     * (-0.024_wp*hicen(jl) + 0.832_wp) 
     908  
     909              !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 
     910              IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &  
     911                                   * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 
     912  
     913              asnon(jl) = reduced_aicen(jl)  ! effective snow fraction (empirical) 
     914              ! MV should check whether this makes sense to have the same effective snow fraction in here 
     915              ! OLI: it probably doesn't 
     916           END IF 
     917  
     918  ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 
     919  ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 
     920  ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 
     921  ! categories.  alfa and beta partition the ITD - they are areas not thicknesses! 
     922  ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 
     923  ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 
     924  ! alfan = 60% of the ice volume) in each category lies above the reference line, 
     925  ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 
     926  
     927  ! MV: 
     928  ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 
     929  ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 
     930  
     931  ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 
     932  
     933           alfan(jl) = 0.6 * hicen(jl) 
     934           betan(jl) = 0.4 * hicen(jl) 
     935  
     936           cum_max_vol(jl)     = 0._wp 
     937           cum_max_vol_tmp(jl) = 0._wp 
     938  
     939        END DO ! jpl 
     940  
     941        cum_max_vol_tmp(0) = 0._wp 
     942        drain = 0._wp 
     943        zdvolp(ji,jj) = 0._wp 
     944  
     945        !---------------------------------------------------------- 
     946        ! Drain overflow water, update pond fraction and volume 
     947        !---------------------------------------------------------- 
     948  
     949        !-------------------------------------------------------------------------- 
     950        ! the maximum amount of water that can be contained up to each ice category 
     951        !-------------------------------------------------------------------------- 
     952        ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 
     953        ! Then the excess volume cum_max_vol(jl) drains out of the system 
     954        ! It should be added to wfx_pnd_out 
     955  
     956        DO jl = 1, jpl-1 ! last category can not hold any volume 
     957  
     958           IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 
     959  
     960              ! total volume in level including snow 
     961              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 
     962                 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 
     963  
     964              ! subtract snow solid volumes from lower categories in current level 
     965              DO ns = 1, jl 
     966                 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 
     967                    - rhos/rhow * &     ! free air fraction that can be filled by water 
     968                      asnon(ns)  * &    ! effective areal fraction of snow in that category 
     969                      max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 
     970              END DO 
     971  
     972           ELSE ! assume higher categories unoccupied 
     973              cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 
     974           END IF 
     975           !IF (cum_max_vol_tmp(jl) < z0) THEN 
     976           !   CALL abort_ice('negative melt pond volume') 
     977           !END IF 
     978        END DO 
     979        cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1)  ! last category holds no volume 
     980        cum_max_vol  (1:jpl) = cum_max_vol_tmp(1:jpl) 
     981  
     982        !---------------------------------------------------------------- 
     983        ! is there more meltwater than can be held in the floe? 
     984        !---------------------------------------------------------------- 
     985        IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 
     986           drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 
     987           zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 
     988  
     989           diag_dvpn_rnf(ji,jj) = - drain      ! diag - overflow counted in the runoff part (arbitrary choice) 
     990            
     991           zdvolp(ji,jj) = drain         ! this is the drained water 
     992           IF (zvolp(ji,jj) < epsi10) THEN 
     993              zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     994              zvolp(ji,jj) = 0._wp 
     995           END IF 
     996        END IF 
     997  
     998        ! height and area corresponding to the remaining volume 
     999        ! routine leaves zvolp unchanged 
     1000        CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     1001  
     1002        DO jl = 1, m_index 
     1003           !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 
     1004           !                                         !  volume instead, no ? 
     1005           h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp)      !js: from CICE 5.1.2 
     1006           a_ip(ji,jj,jl) = reduced_aicen(jl) 
     1007           ! in practise, pond fraction depends on the empirical snow fraction 
     1008           ! so in turn on ice thickness 
     1009        END DO 
     1010        !zapond = sum(a_ip(1:m_index))     !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1011  
     1012        !------------------------------------------------------------------------ 
     1013        ! Drainage through brine network (permeability) 
     1014        !------------------------------------------------------------------------ 
     1015        !!! drainage due to ice permeability - Darcy's law 
     1016  
     1017        ! sea water level 
     1018        msno = 0._wp  
     1019        DO jl = 1 , jpl 
     1020          msno = msno + v_s(ji,jj,jl) * rhos 
     1021        END DO 
     1022        floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 
     1023        hsl_rel = floe_weight / rau0 & 
     1024                - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 
     1025  
     1026        deltah = hpond - hsl_rel 
     1027        pressure_head = grav * rau0 * max(deltah, 0._wp) 
     1028  
     1029        ! drain if ice is permeable 
     1030        permflag = 0 
     1031  
     1032        IF (pressure_head > 0._wp) THEN 
     1033           DO jl = 1, jpl-1 
     1034              IF ( hicen(jl) /= 0._wp ) THEN 
     1035  
     1036              !IF (hicen(jl) > 0._wp) THEN           !js: from CICE 5.1.2 
     1037  
     1038                 perm = 0._wp ! MV ugly dummy patch 
     1039                 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl),  sz_i(ji,jj,:,jl), perm) ! bof  
     1040                 IF (perm > 0._wp) permflag = 1 
     1041  
     1042                 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 
     1043                                          (viscosity*hicen(jl)) 
     1044                 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 
     1045                 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 
     1046  
     1047                 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 
     1048                  
     1049                 IF (zvolp(ji,jj) < epsi10) THEN 
     1050                    zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 
     1051                    zvolp(ji,jj) = 0._wp 
     1052                 END IF 
     1053             END IF 
     1054          END DO 
     1055  
     1056           ! adjust melt pond dimensions 
     1057           IF (permflag > 0) THEN 
     1058              ! recompute pond depth 
     1059             CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 
     1060              DO jl = 1, m_index 
     1061                 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 
     1062                 a_ip(ji,jj,jl) = reduced_aicen(jl) 
     1063              END DO 
     1064              !zapond = sum(a_ip(1:m_index))       !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 
     1065           END IF 
     1066        END IF ! pressure_head 
     1067  
     1068        !------------------------------- 
     1069        ! remove water from the snow 
     1070        !------------------------------- 
     1071        !------------------------------------------------------------------------ 
     1072        ! total melt pond volume in category does not include snow volume 
     1073        ! snow in melt ponds is not melted 
     1074        !------------------------------------------------------------------------ 
     1075         
     1076        ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 
     1077        ! how much, so I did not diagnose it 
     1078        ! so if there is a problem here, nobody is going to see it... 
     1079         
     1080  
     1081        ! Calculate pond volume for lower categories 
     1082        DO jl = 1,m_index-1 
     1083           v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 
     1084                          - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 
     1085        END DO 
     1086  
     1087        ! Calculate pond volume for highest category = remaining pond volume 
     1088  
     1089        ! The following is completely unclear to Martin at least 
     1090        ! Could we redefine properly and recode in a more readable way ? 
     1091  
     1092        ! m_index = last category with melt pond 
     1093  
     1094        IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 
     1095  
     1096        IF (m_index > 1) THEN 
     1097          IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 
     1098             v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 
     1099          ELSE 
     1100             v_ip(ji,jj,m_index) = 0._wp  
     1101             h_ip(ji,jj,m_index) = 0._wp 
     1102             a_ip(ji,jj,m_index) = 0._wp 
     1103             ! If remaining pond volume is negative reduce pond volume of 
     1104             ! lower category 
     1105             IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 
     1106              v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 
     1107          END IF 
     1108        END IF 
     1109  
     1110        DO jl = 1,m_index 
     1111           IF (a_ip(ji,jj,jl) > epsi10) THEN 
     1112               h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 
     1113           ELSE 
     1114              zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 
     1115              h_ip(ji,jj,jl) = 0._wp  
     1116              v_ip(ji,jj,jl)  = 0._wp 
     1117              a_ip(ji,jj,jl) = 0._wp 
     1118           END IF 
     1119        END DO 
     1120        DO jl = m_index+1, jpl 
     1121           h_ip(ji,jj,jl) = 0._wp  
     1122           a_ip(ji,jj,jl) = 0._wp  
     1123           v_ip(ji,jj,jl) = 0._wp  
     1124        END DO 
     1125         
     1126           ENDIF 
     1127        END DO ! ji 
     1128     END DO ! jj 
     1129 
     1130    END SUBROUTINE ice_thd_pnd_area 
     1131 
     1132 
     1133    SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 
     1134       !!------------------------------------------------------------------- 
     1135       !!                ***  ROUTINE ice_thd_pnd_depth  *** 
     1136       !! 
     1137       !! ** Purpose :   Compute melt pond depth 
     1138       !!------------------------------------------------------------------- 
     1139 
     1140       REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 
     1141          aicen, & 
     1142          asnon, & 
     1143          hsnon, & 
     1144          alfan, & 
     1145          cum_max_vol 
     1146 
     1147       REAL (wp), INTENT(IN) :: & 
     1148          zvolp 
     1149 
     1150       REAL (wp), INTENT(OUT) :: & 
     1151          hpond 
     1152 
     1153       INTEGER, INTENT(OUT) :: & 
     1154          m_index 
     1155 
     1156       INTEGER :: n, ns 
     1157 
     1158       REAL (wp), DIMENSION(0:jpl+1) :: & 
     1159          hitl, & 
     1160          aicetl 
     1161 
     1162       REAL (wp) :: & 
     1163          rem_vol, & 
     1164          area, & 
     1165          vol, & 
     1166          tmp, & 
     1167          z0   = 0.0_wp 
     1168 
     1169       !---------------------------------------------------------------- 
     1170       ! hpond is zero if zvolp is zero - have we fully drained? 
     1171       !---------------------------------------------------------------- 
     1172 
     1173       IF (zvolp < epsi10) THEN 
     1174        hpond = z0 
     1175        m_index = 0 
     1176       ELSE 
     1177 
     1178        !---------------------------------------------------------------- 
     1179        ! Calculate the category where water fills up to 
     1180        !---------------------------------------------------------------- 
     1181 
     1182        !----------| 
     1183        !          | 
     1184        !          | 
     1185        !          |----------|                                     -- -- 
     1186        !__________|__________|_________________________________________ ^ 
     1187        !          |          |             rem_vol     ^                | Semi-filled 
     1188        !          |          |----------|-- -- -- - ---|-- ---- -- -- --v layer 
     1189        !          |          |          |              | 
     1190        !          |          |          |              |hpond 
     1191        !          |          |          |----------|   |     |------- 
     1192        !          |          |          |          |   |     | 
     1193        !          |          |          |          |---v-----| 
     1194        !          |          | m_index  |          |         | 
     1195        !------------------------------------------------------------- 
     1196 
     1197        m_index = 0  ! 1:m_index categories have water in them 
     1198        DO n = 1, jpl 
     1199           IF (zvolp <= cum_max_vol(n)) THEN 
     1200              m_index = n 
     1201              IF (n == 1) THEN 
     1202                 rem_vol = zvolp 
     1203              ELSE 
     1204                 rem_vol = zvolp - cum_max_vol(n-1) 
     1205              END IF 
     1206              exit ! to break out of the loop 
     1207           END IF 
     1208        END DO 
     1209        m_index = min(jpl-1, m_index) 
     1210 
     1211        !---------------------------------------------------------------- 
     1212        ! semi-filled layer may have m_index different snow in it 
     1213        !---------------------------------------------------------------- 
     1214 
     1215        !-----------------------------------------------------------  ^ 
     1216        !                                                             |  alfan(m_index+1) 
     1217        !                                                             | 
     1218        !hitl(3)-->                             |----------|          | 
     1219        !hitl(2)-->                |------------| * * * * *|          | 
     1220        !hitl(1)-->     |----------|* * * * * * |* * * * * |          | 
     1221        !hitl(0)-->-------------------------------------------------  |  ^ 
     1222        !                various snow from lower categories          |  |alfa(m_index) 
     1223 
     1224        ! hitl - heights of the snow layers from thinner and current categories 
     1225        ! aicetl - area of each snow depth in this layer 
     1226 
     1227        hitl(:) = z0 
     1228        aicetl(:) = z0 
     1229        DO n = 1, m_index 
     1230           hitl(n)   = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 
     1231                                  alfan(m_index+1) - alfan(m_index)), z0) 
     1232           aicetl(n) = asnon(n) 
     1233 
     1234           aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 
     1235        END DO 
     1236 
     1237        hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 
     1238        aicetl(m_index+1) = z0 
     1239 
     1240        !---------------------------------------------------------------- 
     1241        ! reorder array according to hitl 
     1242        ! snow heights not necessarily in height order 
     1243        !---------------------------------------------------------------- 
     1244 
     1245        DO ns = 1, m_index+1 
     1246           DO n = 0, m_index - ns + 1 
     1247              IF (hitl(n) > hitl(n+1)) THEN ! swap order 
     1248                 tmp = hitl(n) 
     1249                 hitl(n) = hitl(n+1) 
     1250                 hitl(n+1) = tmp 
     1251                 tmp = aicetl(n) 
     1252                 aicetl(n) = aicetl(n+1) 
     1253                 aicetl(n+1) = tmp 
     1254              END IF 
     1255           END DO 
     1256        END DO 
     1257 
     1258        !---------------------------------------------------------------- 
     1259        ! divide semi-filled layer into set of sublayers each vertically homogenous 
     1260        !---------------------------------------------------------------- 
     1261 
     1262        !hitl(3)---------------------------------------------------------------- 
     1263        !                                                       | * * * * * * * * 
     1264        !                                                       |* * * * * * * * * 
     1265        !hitl(2)---------------------------------------------------------------- 
     1266        !                                    | * * * * * * * *  | * * * * * * * * 
     1267        !                                    |* * * * * * * * * |* * * * * * * * * 
     1268        !hitl(1)---------------------------------------------------------------- 
     1269        !                 | * * * * * * * *  | * * * * * * * *  | * * * * * * * * 
     1270        !                 |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 
     1271        !hitl(0)---------------------------------------------------------------- 
     1272        !    aicetl(0)         aicetl(1)           aicetl(2)          aicetl(3) 
     1273 
     1274        ! move up over layers incrementing volume 
     1275        DO n = 1, m_index+1 
     1276 
     1277           area = sum(aicetl(:)) - &                 ! total area of sub-layer 
     1278                (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 
     1279 
     1280           vol = (hitl(n) - hitl(n-1)) * area      ! thickness of sub-layer times area 
     1281 
     1282           IF (vol >= rem_vol) THEN  ! have reached the sub-layer with the depth within 
     1283              hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 
     1284 
     1285              exit 
     1286           ELSE  ! still in sub-layer below the sub-layer with the depth 
     1287              rem_vol = rem_vol - vol 
     1288           END IF 
     1289 
     1290        END DO 
     1291 
     1292       END IF 
     1293 
     1294    END SUBROUTINE ice_thd_pnd_depth 
     1295 
     1296 
     1297    SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 
     1298       !!------------------------------------------------------------------- 
     1299       !!                ***  ROUTINE ice_thd_pnd_perm *** 
     1300       !! 
     1301       !! ** Purpose :   Determine the liquid fraction of brine in the ice 
     1302       !!                and its permeability 
     1303       !!------------------------------------------------------------------- 
     1304 
     1305       REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 
     1306          ticen,  &     ! internal ice temperature (K) 
     1307          salin         ! salinity (ppt)     !js: ppt according to cice 
     1308 
     1309       REAL (wp), INTENT(OUT) :: & 
     1310          perm      ! permeability 
     1311 
     1312       REAL (wp) ::   & 
     1313          Sbr       ! brine salinity 
     1314 
     1315       REAL (wp), DIMENSION(nlay_i) ::   & 
     1316          Tin, &    ! ice temperature 
     1317          phi       ! liquid fraction 
     1318 
     1319       INTEGER :: k 
     1320 
     1321       !----------------------------------------------------------------- 
     1322       ! Compute ice temperatures from enthalpies using quadratic formula 
     1323       !----------------------------------------------------------------- 
     1324 
     1325       DO k = 1,nlay_i 
     1326          Tin(k) = ticen(k) - rt0   !js: from K to degC 
     1327       END DO 
     1328 
     1329       !----------------------------------------------------------------- 
     1330       ! brine salinity and liquid fraction 
     1331       !----------------------------------------------------------------- 
     1332 
     1333       DO k = 1, nlay_i 
     1334        
     1335          ! Sbr    = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 
     1336          ! Best expression to date is that one 
     1337          Sbr  = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 
     1338          phi(k) = salin(k) / Sbr 
     1339           
     1340       END DO 
     1341 
     1342       !----------------------------------------------------------------- 
     1343       ! permeability 
     1344       !----------------------------------------------------------------- 
     1345 
     1346       perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 
     1347 
     1348   END SUBROUTINE ice_thd_pnd_perm 
     1349    
     1350    
     1351!---------------------------------------------------------------------------------------------------------------------- 
    3041352 
    3051353   SUBROUTINE ice_thd_pnd_init  
    … …  
    3171365      INTEGER  ::   ios, ioptio   ! Local integer 
    3181366      !! 
    319       NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 
    320          &                          ln_pnd_CST , rn_apnd, rn_hpnd,         & 
     1367      NAMELIST/namthd_pnd/  ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 
     1368         &                          ln_pnd_CST , rn_apnd, rn_hpnd,                       & 
     1369         &                          ln_pnd_TOPO ,                                        & 
    3211370         &                          ln_pnd_lids, ln_pnd_alb 
    3221371      !!------------------------------------------------------------------- 
    … …  
    3361385         WRITE(numout,*) '   Namelist namicethd_pnd:' 
    3371386         WRITE(numout,*) '      Melt ponds activated or not                                 ln_pnd       = ', ln_pnd 
     1387         WRITE(numout,*) '         Topographic melt pond scheme                             ln_pnd_TOPO  = ', ln_pnd_TOPO 
    3381388         WRITE(numout,*) '         Level ice melt pond scheme                               ln_pnd_LEV   = ', ln_pnd_LEV 
    3391389         WRITE(numout,*) '            Minimum ice fraction that contributes to melt ponds   rn_apnd_min  = ', rn_apnd_min 
    3401390         WRITE(numout,*) '            Maximum ice fraction that contributes to melt ponds   rn_apnd_max  = ', rn_apnd_max 
     1391         WRITE(numout,*) '            Pond flushing efficiency                              rn_pnd_flush = ', rn_pnd_flush 
    3411392         WRITE(numout,*) '         Constant ice melt pond scheme                            ln_pnd_CST   = ', ln_pnd_CST 
    3421393         WRITE(numout,*) '            Prescribed pond fraction                              rn_apnd      = ', rn_apnd 
    … …  
    3511402      IF( ln_pnd_CST  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndCST    ;   ENDIF 
    3521403      IF( ln_pnd_LEV  ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndLEV    ;   ENDIF 
     1404      IF( ln_pnd_TOPO ) THEN   ;   ioptio = ioptio + 1   ;   nice_pnd = np_pndTOPO   ;   ENDIF 
    3531405      IF( ioptio /= 1 )   & 
    354          & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' ) 
     1406         & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV, ln_pnd_CST or ln_pnd_TOPO)' ) 
    3551407      ! 
    3561408      SELECT CASE( nice_pnd ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/iceupdate.F90

    r14075 r15658  
    146146            ! 
    147147            ! the non-solar is simply derived from the solar flux 
    148             qns(ji,jj) = qt_oce_ai(ji,jj) - zqsr               
     148            qns(ji,jj) = qt_oce_ai(ji,jj) - qsr(ji,jj)               
    149149 
    150150            ! Mass flux at the atm. surface        

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icevar.F90

    r14075 r15658  
    244244      ELSEWHERE( h_il(:,:,:) >= zhl_max )  ;   a_ip_eff(:,:,:) = 0._wp                  ! lid is very thick. Cover all the pond up with ice and snow 
    245245      ELSEWHERE                            ;   a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * &   ! lid is in between. Expose part of the pond 
    246          &                                                       ( h_il(:,:,:) - zhl_min ) / ( zhl_max - zhl_min ) 
     246         &                                                       ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 
    247247      END WHERE 
    248248      ! 
    … …  
    565565            END DO 
    566566         END DO 
     567 
     568         !----------------------------------------------------------------- 
     569         ! zap small ponds 
     570         !----------------------------------------------------------------- 
     571         DO jj = 1, jpj 
     572            DO ji = 1, jpi 
     573               IF ( v_ip(ji,jj,jl) <= epsi10 ) THEN 
     574                  a_ip(ji,jj,jl)      = 0._wp 
     575                  a_ip_frac(ji,jj,jl) = 0._wp 
     576                  v_ip(ji,jj,jl)      = 0._wp 
     577                  h_ip(ji,jj,jl)      = 0._wp 
     578                  v_il(ji,jj,jl)      = 0._wp 
     579                  h_il(ji,jj,jl)      = 0._wp 
     580                ENDIF 
     581            END DO 
     582         END DO 
    567583         ! 
    568584      END DO  
    … …  
    696712      WHERE( pe_i (1:npti,:,:) < 0._wp .AND. pe_i (1:npti,:,:) > -epsi06 )   pe_i (1:npti,:,:) = 0._wp   !  e_i must be >= 0 
    697713      WHERE( pe_s (1:npti,:,:) < 0._wp .AND. pe_s (1:npti,:,:) > -epsi06 )   pe_s (1:npti,:,:) = 0._wp   !  e_s must be >= 0 
    698       IF( ln_pnd_LEV ) THEN 
     714      IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 
    699715         WHERE( pa_ip(1:npti,:) < 0._wp .AND. pa_ip(1:npti,:) > -epsi10 )    pa_ip(1:npti,:)   = 0._wp   ! a_ip must be >= 0 
    700716         WHERE( pv_ip(1:npti,:) < 0._wp .AND. pv_ip(1:npti,:) > -epsi10 )    pv_ip(1:npti,:)   = 0._wp   ! v_ip must be >= 0 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icewri.F90

    r14075 r15658  
    104104      IF( iom_use('snwthic' ) )   CALL iom_put( 'snwthic', hm_s        * zmsk00 )                                           ! snw thickness 
    105105      IF( iom_use('icebrv'  ) )   CALL iom_put( 'icebrv' , bvm_i* 100. * zmsk00 )                                           ! brine volume 
    106       IF( iom_use('iceage'  ) )   CALL iom_put( 'iceage' , om_i / rday * zmsk15 + zmiss_val * ( 1._wp - zmsk15 ) )          ! ice age 
     106      IF( iom_use('iceage'  ) )   CALL iom_put( 'iceage' , om_i / rday * zmsk15 )                                           ! ice age 
    107107      IF( iom_use('icehnew' ) )   CALL iom_put( 'icehnew', ht_i_new             )                                           ! new ice thickness formed in the leads 
    108108      IF( iom_use('snwvolu' ) )   CALL iom_put( 'snwvolu', vt_s        * zmsksn )                                           ! snow volume 
    … …  
    116116      IF( iom_use('icehpnd' ) )   CALL iom_put( 'icehpnd', hm_ip  * zmsk00      )                                           ! melt pond depth 
    117117      IF( iom_use('icevpnd' ) )   CALL iom_put( 'icevpnd', vt_ip  * zmsk00      )                                           ! melt pond total volume per unit area 
     118      IF( iom_use('iceepnd' ) )   CALL iom_put( 'iceepnd', SUM( a_ip_eff * a_i, dim=3 ) * zmsk00  )                         ! melt pond total effective fraction per cell area 
    118119      IF( iom_use('icehlid' ) )   CALL iom_put( 'icehlid', hm_il  * zmsk00      )                                           ! melt pond lid depth 
    119120      IF( iom_use('icevlid' ) )   CALL iom_put( 'icevlid', vt_il  * zmsk00      )                                           ! melt pond lid total volume per unit area 
    120121      ! salt 
    121       IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) ! mean ice salinity 
     122      IF( iom_use('icesalt' ) )   CALL iom_put( 'icesalt', sm_i                 * zmsk00 )                                  ! mean ice salinity 
    122123      IF( iom_use('icesalm' ) )   CALL iom_put( 'icesalm', st_i * rhoi * 1.0e-3 * zmsk00 )                                  ! Mass of salt in sea ice per cell area 
    123124      ! heat 
    124       IF( iom_use('icetemp' ) )   CALL iom_put( 'icetemp', ( tm_i  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! ice mean temperature 
    125       IF( iom_use('snwtemp' ) )   CALL iom_put( 'snwtemp', ( tm_s  - rt0 ) * zmsksn + zmiss_val * ( 1._wp - zmsksn ) )      ! snw mean temperature 
    126       IF( iom_use('icettop' ) )   CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice surface 
    127       IF( iom_use('icetbot' ) )   CALL iom_put( 'icetbot', ( t_bo  - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the ice bottom 
    128       IF( iom_use('icetsni' ) )   CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) )      ! temperature at the snow-ice interface 
     125      IF( iom_use('icetemp' ) )   CALL iom_put( 'icetemp', ( tm_i  - rt0 ) * zmsk00  )                                      ! ice mean temperature 
     126      IF( iom_use('snwtemp' ) )   CALL iom_put( 'snwtemp', ( tm_s  - rt0 ) * zmsksn  )                                      ! snw mean temperature 
     127      IF( iom_use('icettop' ) )   CALL iom_put( 'icettop', ( tm_su - rt0 ) * zmsk00  )                                      ! temperature at the ice surface 
     128      IF( iom_use('icetbot' ) )   CALL iom_put( 'icetbot', ( t_bo  - rt0 ) * zmsk00  )                                      ! temperature at the ice bottom 
     129      IF( iom_use('icetsni' ) )   CALL iom_put( 'icetsni', ( tm_si - rt0 ) * zmsk00  )                                      ! temperature at the snow-ice interface 
    129130      IF( iom_use('icehc'   ) )   CALL iom_put( 'icehc'  ,  -et_i          * zmsk00 )                                       ! ice heat content 
    130131      IF( iom_use('snwhc'   ) )   CALL iom_put( 'snwhc'  ,  -et_s          * zmsksn )                                       ! snow heat content 
    … …  
    153154      IF( iom_use('icemask_cat' ) )   CALL iom_put( 'icemask_cat' ,                  zmsk00l                                   ) ! ice mask 0% 
    154155      IF( iom_use('iceconc_cat' ) )   CALL iom_put( 'iceconc_cat' , a_i            * zmsk00l                                   ) ! area for categories 
    155       IF( iom_use('icethic_cat' ) )   CALL iom_put( 'icethic_cat' , h_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! thickness for categories 
    156       IF( iom_use('snwthic_cat' ) )   CALL iom_put( 'snwthic_cat' , h_s            * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow depth for categories 
    157       IF( iom_use('icesalt_cat' ) )   CALL iom_put( 'icesalt_cat' , s_i            * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! salinity for categories 
    158       IF( iom_use('iceage_cat'  ) )   CALL iom_put( 'iceage_cat'  , o_i / rday     * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice age 
     156      IF( iom_use('icethic_cat' ) )   CALL iom_put( 'icethic_cat' , h_i            * zmsk00l                                   ) ! thickness for categories 
     157      IF( iom_use('snwthic_cat' ) )   CALL iom_put( 'snwthic_cat' , h_s            * zmsksnl                                   ) ! snow depth for categories 
     158      IF( iom_use('icesalt_cat' ) )   CALL iom_put( 'icesalt_cat' , s_i            * zmsk00l                                   ) ! salinity for categories 
     159      IF( iom_use('iceage_cat'  ) )   CALL iom_put( 'iceage_cat'  , o_i / rday     * zmsk00l                                   ) ! ice age 
    159160      IF( iom_use('icetemp_cat' ) )   CALL iom_put( 'icetemp_cat' , ( SUM( t_i, dim=3 ) * r1_nlay_i - rt0 ) & 
    160          &                                                                         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice temperature 
     161         &                                                                         * zmsk00l                                   ) ! ice temperature 
    161162      IF( iom_use('snwtemp_cat' ) )   CALL iom_put( 'snwtemp_cat' , ( SUM( t_s, dim=3 ) * r1_nlay_s - rt0 ) & 
    162          &                                                                         * zmsksnl + zmiss_val * ( 1._wp - zmsksnl ) ) ! snow temperature 
    163       IF( iom_use('icettop_cat' ) )   CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! surface temperature 
    164       IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 
     163         &                                                                         * zmsksnl                                   ) ! snow temperature 
     164      IF( iom_use('icettop_cat' ) )   CALL iom_put( 'icettop_cat' , ( t_su - rt0 ) * zmsk00l                                   ) ! surface temperature 
     165      IF( iom_use('icebrv_cat'  ) )   CALL iom_put( 'icebrv_cat'  ,   bv_i * 100.  * zmsk00l                                   ) ! brine volume 
    165166      IF( iom_use('iceapnd_cat' ) )   CALL iom_put( 'iceapnd_cat' ,   a_ip         * zmsk00l                                   ) ! melt pond frac for categories 
     167      IF( iom_use('icevpnd_cat' ) )   CALL iom_put( 'icevpnd_cat' ,   v_ip         * zmsk00l                                   ) ! melt pond volume for categories 
    166168      IF( iom_use('icehpnd_cat' ) )   CALL iom_put( 'icehpnd_cat' ,   h_ip         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories 
    167169      IF( iom_use('icehlid_cat' ) )   CALL iom_put( 'icehlid_cat' ,   h_il         * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories 
    168       IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac for categories 
     170      IF( iom_use('iceafpnd_cat') )   CALL iom_put( 'iceafpnd_cat',   a_ip_frac    * zmsk00l                                   ) ! melt pond frac per ice area for categories 
    169171      IF( iom_use('iceaepnd_cat') )   CALL iom_put( 'iceaepnd_cat',   a_ip_eff     * zmsk00l                                   ) ! melt pond effective frac for categories 
    170172      IF( iom_use('icealb_cat'  ) )   CALL iom_put( 'icealb_cat'  ,   alb_ice      * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/DIA/diawri.F90

    r14075 r15658  
    4949   USE iom            !  
    5050   USE ioipsl         !  
    51  
     51   USE eosbn2 
    5252#if defined key_si3 
    5353   USE ice  
    … …  
    113113      REAL(wp), DIMENSION(jpi,jpj)     ::   z2d   ! 2D workspace 
    114114      REAL(wp), DIMENSION(jpi,jpj,jpk) ::   z3d   ! 3D workspace 
     115      CHARACTER(len=4),SAVE :: ttype , stype           ! temperature and salinity type 
    115116      !!---------------------------------------------------------------------- 
    116117      !  
     118      IF( kt == nit000 ) THEN 
     119         IF( ln_TEOS10 ) THEN 
     120            IF ( iom_use("toce_pot") .OR. iom_use("soce_pra") .OR. iom_use("sst_pot") .OR. iom_use("sss_pra") & 
     121                  & .OR. iom_use("sbt_pot") .OR. iom_use("sbs_pra") .OR. iom_use("sstgrad_pot") .OR. iom_use("sstgrad2_pot") & 
     122                  & .OR. iom_use("tosmint_pot") .OR. iom_use("somint_pra"))  THEN  
     123               CALL ctl_stop( 'diawri: potential temperature and practical salinity not available with ln_TEOS10' ) 
     124            ELSE 
     125               ttype='con' ; stype='abs'   ! teos-10 using conservative temperature and absolute salinity 
     126            ENDIF  
     127         ELSE IF( ln_EOS80  ) THEN 
     128            IF ( iom_use("toce_con") .OR. iom_use("soce_abs") .OR. iom_use("sst_con") .OR. iom_use("sss_abs") & 
     129                  & .OR. iom_use("sbt_con") .OR. iom_use("sbs_abs") .OR. iom_use("sstgrad_con") .OR. iom_use("sstgrad2_con") & 
     130                  & .OR. iom_use("tosmint_con") .OR. iom_use("somint_abs"))  THEN  
     131               CALL ctl_stop( 'diawri: conservative temperature and absolute salinity not available with ln_EOS80' ) 
     132            ELSE 
     133               ttype='pot' ; stype='pra'   ! eos-80 using potential temperature and practical salinity 
     134            ENDIF 
     135         ELSE IF ( ln_SEOS) THEN 
     136            ttype='seos' ; stype='seos' ! seos using Simplified Equation of state 
     137         ENDIF 
     138      ENDIF 
     139 
    117140      IF( ln_timing )   CALL timing_start('dia_wri') 
    118141      !  
    … …  
    144167         CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) ) 
    145168       
    146       CALL iom_put( "toce", tsn(:,:,:,jp_tem) )    ! 3D temperature 
    147       CALL iom_put(  "sst", tsn(:,:,1,jp_tem) )    ! surface temperature 
    148       IF ( iom_use("sbt") ) THEN 
     169      CALL iom_put( "toce_"//ttype, tsn(:,:,:,jp_tem) )    ! 3D temperature 
     170      CALL iom_put(  "sst_"//ttype, tsn(:,:,1,jp_tem) )    ! surface temperature 
     171      IF ( iom_use("sbt_"//ttype) ) THEN 
    149172         DO jj = 1, jpj 
    150173            DO ji = 1, jpi 
    … …  
    153176            END DO 
    154177         END DO 
    155          CALL iom_put( "sbt", z2d )                ! bottom temperature 
     178         CALL iom_put( "sbt_"//ttype, z2d )                ! bottom temperature 
    156179      ENDIF 
    157180       
    158       CALL iom_put( "soce", tsn(:,:,:,jp_sal) )    ! 3D salinity 
    159       CALL iom_put(  "sss", tsn(:,:,1,jp_sal) )    ! surface salinity 
    160       IF ( iom_use("sbs") ) THEN 
     181      CALL iom_put( "soce_"//stype, tsn(:,:,:,jp_sal) )    ! 3D salinity 
     182      CALL iom_put(  "sss_"//stype, tsn(:,:,1,jp_sal) )    ! surface salinity 
     183      IF ( iom_use("sbs_"//stype) ) THEN 
    161184         DO jj = 1, jpj 
    162185            DO ji = 1, jpi 
    … …  
    165188            END DO 
    166189         END DO 
    167          CALL iom_put( "sbs", z2d )                ! bottom salinity 
     190         CALL iom_put( "sbs_"//stype, z2d )                ! bottom salinity 
    168191      ENDIF 
    169192 
    … …  
    233256      IF( iom_use('logavs') )   CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) ) 
    234257 
    235       IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN 
     258      IF ( iom_use("sstgrad_"//ttype) .OR. iom_use("sstgrad2_"//ttype) ) THEN 
    236259         DO jj = 2, jpjm1                                    ! sst gradient 
    237260            DO ji = fs_2, fs_jpim1   ! vector opt. 
    … …  
    244267         END DO 
    245268         CALL lbc_lnk( 'diawri', z2d, 'T', 1. ) 
    246          CALL iom_put( "sstgrad2",  z2d )          ! square of module of sst gradient 
     269         CALL iom_put( "sstgrad2_"//ttype,  z2d )          ! square of module of sst gradient 
    247270         z2d(:,:) = SQRT( z2d(:,:) ) 
    248          CALL iom_put( "sstgrad" ,  z2d )          ! module of sst gradient 
     271         CALL iom_put( "sstgrad_"//ttype ,  z2d )          ! module of sst gradient 
    249272      ENDIF 
    250273          
    … …  
    365388      ENDIF 
    366389 
    367       IF( iom_use("tosmint") ) THEN 
     390      IF( (.NOT.l_ldfeiv_time) .AND. ( iom_use('RossRad')  .OR. iom_use('RossRadlim') & 
     391            &                     .OR. iom_use('Tclinic_recip') .OR. iom_use('RR_GS')      & 
     392            &                     .OR. iom_use('aeiu_2d')  .OR. iom_use('aeiv_2d') ) ) THEN 
     393         CALL ldf_eiv(kt, 75.0, z2d, z3d(:,:,1)) 
     394         CALL iom_put('aeiu_2d', z2d) 
     395         CALL iom_put('aeiv_2d', z3d(:,:,1)) 
     396      ENDIF 
     397 
     398      IF( iom_use("tosmint_"//ttype) ) THEN 
    368399         z2d(:,:) = 0._wp 
    369400         DO jk = 1, jpkm1 
    … …  
    375406         END DO 
    376407         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    377          CALL iom_put( "tosmint", rau0 * z2d )        ! Vertical integral of temperature 
    378       ENDIF 
    379       IF( iom_use("somint") ) THEN 
     408         CALL iom_put( "tosmint_"//ttype, rau0 * z2d )        ! Vertical integral of temperature 
     409      ENDIF 
     410      IF( iom_use("somint_"//stype) ) THEN 
    380411         z2d(:,:)=0._wp 
    381412         DO jk = 1, jpkm1 
    … …  
    387418         END DO 
    388419         CALL lbc_lnk( 'diawri', z2d, 'T', -1. ) 
    389          CALL iom_put( "somint", rau0 * z2d )         ! Vertical integral of salinity 
     420         CALL iom_put( "somint_"//stype, rau0 * z2d )         ! Vertical integral of salinity 
    390421      ENDIF 
    391422 
    … …  
    930961         CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd            )    ! now StokesDrift k-velocity 
    931962      ENDIF 
    932   
    933963#if defined key_si3 
    934964      IF( nn_ice == 2 ) THEN   ! condition needed in case agrif + ice-model but no-ice in child grid 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/DOM/domain.F90

    r14075 r15658  
    292292         &             nn_it000, nn_itend , nn_date0    , nn_time0     , nn_leapy  , nn_istate ,     & 
    293293         &             nn_stock, nn_write , ln_mskland  , ln_clobber   , nn_chunksz, nn_euler  ,     & 
    294          &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios 
     294         &             ln_cfmeta, ln_iscpl, ln_xios_read, nn_wxios, ln_rstdate, ln_rst_eos 
     295 
    295296      NAMELIST/namdom/ ln_linssh, rn_isfhmin, rn_rdt, rn_atfp, ln_crs, ln_meshmask 
    296297#if defined key_netcdf4 
    … …  
    340341#endif 
    341342         WRITE(numout,*) '      mask land points                ln_mskland      = ', ln_mskland 
     343         WRITE(numout,*) '      date-stamp restart files        ln_rstdate = ', ln_rstdate 
    342344         WRITE(numout,*) '      additional CF standard metadata ln_cfmeta       = ', ln_cfmeta 
    343345         WRITE(numout,*) '      overwrite an existing file      ln_clobber      = ', ln_clobber 
    344346         WRITE(numout,*) '      NetCDF chunksize (bytes)        nn_chunksz      = ', nn_chunksz 
    345347         WRITE(numout,*) '      IS coupling at the restart step ln_iscpl        = ', ln_iscpl 
     348         WRITE(numout,*) '      check restart equation of state ln_rst_eos      = ', ln_rst_eos 
     349 
    346350         IF( TRIM(Agrif_CFixed()) == '0' ) THEN 
    347351            WRITE(numout,*) '      READ restart for a single file using XIOS ln_xios_read =', ln_xios_read 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/DOM/dommsk.F90

    r14075 r15658  
    3232   USE lbclnk         ! ocean lateral boundary conditions (or mpp link) 
    3333   USE lib_mpp        ! Massively Parallel Processing library 
     34   USE iom             ! For shlat2d 
     35   USE fldread         ! for sn_shlat2d 
    3436 
    3537   IMPLICIT NONE 
    … …  
    9395      INTEGER  ::   ios, inum 
    9496      !! 
    95       NAMELIST/namlbc/ rn_shlat, ln_vorlat 
     97      REAL(wp) :: zshlat           !: locally modified shlat for some strait 
     98      REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: zshlat2d 
     99      LOGICAL                         :: ln_shlat2d 
     100      CHARACTER(len = 256)            :: cn_shlat2d_file, cn_shlat2d_var   
     101      !! 
     102      NAMELIST/namlbc/ rn_shlat, ln_vorlat, ln_shlat2d, cn_shlat2d_file, cn_shlat2d_var 
    96103      NAMELIST/nambdy/ ln_bdy ,nb_bdy, ln_coords_file, cn_coords_file,         & 
    97104         &             ln_mask_file, cn_mask_file, cn_dyn2d, nn_dyn2d_dta,     & 
    … …  
    120127      ! 
    121128      IF(lwp) WRITE(numout,*) 
    122       IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  free-slip' 
    123       ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  no-slip' 
    124       ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  partial-slip' 
    125       ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  strong-slip' 
     129 
     130      IF ( ln_shlat2d ) THEN 
     131         IF(lwp) WRITE(numout,*) '         READ shlat as a 2D coefficient in a file ' 
     132         ALLOCATE( zshlat2d(jpi,jpj) ) 
     133         CALL iom_open(TRIM(cn_shlat2d_file), inum) 
     134         CALL iom_get (inum, jpdom_data, TRIM(cn_shlat2d_var), zshlat2d, 1) ! 
     135         CALL iom_close(inum) 
    126136      ELSE 
    127          CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' ) 
     137         IF     (      rn_shlat == 0.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  free-slip' 
     138         ELSEIF (      rn_shlat == 2.               ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  no-slip' 
     139         ELSEIF ( 0. < rn_shlat .AND. rn_shlat < 2. ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  partial-slip' 
     140         ELSEIF ( 2. < rn_shlat                     ) THEN   ;   IF(lwp) WRITE(numout,*) '   ==>>>   ocean lateral  strong-slip' 
     141         ELSE 
     142            CALL ctl_stop( 'dom_msk: wrong value for rn_shlat (i.e. a negalive value). We stop.' ) 
     143         ENDIF 
    128144      ENDIF 
    129145 
    … …  
    240256      ! Lateral boundary conditions on velocity (modify fmask) 
    241257      ! ---------------------------------------   
    242       IF( rn_shlat /= 0 ) THEN      ! Not free-slip lateral boundary condition 
     258      IF( rn_shlat /= 0 .or. ln_shlat2d ) THEN      ! Not free-slip lateral boundary condition everywhere 
    243259         ! 
    244260         DO jk = 1, jpk 
    245             DO jj = 2, jpjm1 
    246                DO ji = fs_2, fs_jpim1   ! vector opt. 
    247                   IF( fmask(ji,jj,jk) == 0._wp ) THEN 
    248                      fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk), & 
    249                         &                                           vmask(ji,jj,jk), vmask(ji+1,jj,jk) ) ) 
    250                   ENDIF 
     261            IF (  ln_shlat2d ) THEN 
     262               DO jj = 2, jpjm1 
     263                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     264                     IF( fmask(ji,jj,jk) == 0._wp ) THEN 
     265                        fmask(ji,jj,jk) = zshlat2d(ji,jj) * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk),   & 
     266                           &                                                  vmask(ji,jj,jk), vmask(ji+1,jj,jk)  )  ) 
     267                     ENDIF 
     268                  END DO 
    251269               END DO 
    252             END DO 
     270            ELSE 
     271               DO jj = 2, jpjm1 
     272                  DO ji = fs_2, fs_jpim1   ! vector opt. 
     273                     IF( fmask(ji,jj,jk) == 0._wp ) THEN 
     274                        fmask(ji,jj,jk) = rn_shlat * MIN( 1._wp , MAX( umask(ji,jj,jk), umask(ji,jj+1,jk),   & 
     275                           &                                           vmask(ji,jj,jk), vmask(ji+1,jj,jk)   )  ) 
     276                     ENDIF 
     277                  END DO 
     278               END DO 
     279            ENDIF 
    253280            DO jj = 2, jpjm1 
    254281               IF( fmask(1,jj,jk) == 0._wp ) THEN 
    … …  
    277304         END DO 
    278305         ! 
     306         IF( ln_shlat2d ) DEALLOCATE( zshlat2d ) 
     307         ! 
    279308         CALL lbc_lnk( 'dommsk', fmask, 'F', 1._wp )      ! Lateral boundary conditions on fmask 
    280309         ! 
    … …  
    284313       
    285314      ! User defined alteration of fmask (use to reduce ocean transport in specified straits) 
     315      ! Only call if we are not using the shlat2d option. 
    286316      ! --------------------------------  
    287317      ! 
    288       CALL usr_def_fmask( cn_cfg, nn_cfg, fmask ) 
     318      IF ( .not. ln_shlat2d ) THEN       
     319         CALL usr_def_fmask( cn_cfg, nn_cfg, fmask ) 
     320      ENDIF 
    289321      ! 
    290322   END SUBROUTINE dom_msk 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/ICB/icbdia.F90

    r14075 r15658  
    8383   REAL(wp)                      ::  heat_to_bergs_net, heat_to_ocean_net, melt_net 
    8484   REAL(wp)                      ::  calving_to_bergs_net 
     85   REAL(wp)                      ::  vel_factor_min 
    8586 
    8687   INTEGER                       ::  nbergs_start, nbergs_end, nbergs_calved 
    8788   INTEGER                       ::  nbergs_melted 
    88    INTEGER                       ::  nspeeding_tickets 
     89   INTEGER , DIMENSION(4)        ::  nspeeding_tickets 
    8990   INTEGER , DIMENSION(nclasses) ::  nbergs_calved_by_class 
    9091 
    … …  
    124125      nbergs_calved             = 0 
    125126      nbergs_calved_by_class(:) = 0 
    126       nspeeding_tickets         = 0 
     127      nspeeding_tickets(:)      = 0 
     128      vel_factor_min            = 1._wp 
    127129      stored_heat_end           = 0._wp 
    128130      floating_heat_end         = 0._wp 
    … …  
    157159      IF( lk_mpp ) THEN 
    158160         ALLOCATE( rsumbuf(23) )          ; rsumbuf(:) = 0._wp 
    159          ALLOCATE( nsumbuf(4+nclasses) )  ; nsumbuf(:) = 0 
     161         ALLOCATE( nsumbuf(7+nclasses) )  ; nsumbuf(:) = 0 
    160162         rsumbuf(1) = floating_mass_start 
    161163         rsumbuf(2) = bergs_mass_start 
    … …  
    265267            nsumbuf(2) = nbergs_calved 
    266268            nsumbuf(3) = nbergs_melted 
    267             nsumbuf(4) = nspeeding_tickets 
     269            nsumbuf(4:7) = nspeeding_tickets(:) 
    268270            DO ik = 1, nclasses 
    269                nsumbuf(4+ik) = nbergs_calved_by_class(ik) 
     271               nsumbuf(7+ik) = nbergs_calved_by_class(ik) 
    270272            END DO 
    271             CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+4), nclasses+4 ) 
     273            CALL mpp_sum( 'icbdia', nsumbuf(1:nclasses+7), nclasses+7 ) 
    272274            ! 
    273275            nbergs_end        = nsumbuf(1) 
    274276            nbergs_calved     = nsumbuf(2) 
    275277            nbergs_melted     = nsumbuf(3) 
    276             nspeeding_tickets = nsumbuf(4) 
     278            nspeeding_tickets(:) = nsumbuf(4:7) 
    277279            DO ik = 1,nclasses 
    278                nbergs_calved_by_class(ik)= nsumbuf(4+ik) 
     280               nbergs_calved_by_class(ik)= nsumbuf(7+ik) 
    279281            END DO 
    280282            ! 
     283            CALL mpp_min( 'icbdia', vel_factor_min, 1 ) 
    281284         ENDIF 
    282285         ! 
    … …  
    329332         IF (nn_verbose_level > 0) THEN 
    330333            WRITE( numicb, '("calved by class = ",i6,20(",",i6))') (nbergs_calved_by_class(ik),ik=1,nclasses) 
    331             IF( nspeeding_tickets > 0 )   WRITE( numicb, '("speeding tickets issued = ",i6)') nspeeding_tickets 
     334            WRITE( numicb, '("n speeding tickets by RK4 stage = ",i6,3(",",i6))') (nspeeding_tickets(ik),ik=1,4) 
     335            IF( SUM(nspeeding_tickets) > 0 ) THEN 
     336               WRITE( numicb, '("min velocity reduction factor = ",f12.8)') vel_factor_min 
     337            ENDIF 
    332338         ENDIF 
    333339         ! 
    … …  
    337343         nbergs_calved             = 0 
    338344         nbergs_calved_by_class(:) = 0 
    339          nspeeding_tickets         = 0 
     345         nspeeding_tickets(:)      = 0 
     346         vel_factor_min            = 1._wp 
    340347         stored_heat_start         = stored_heat_end 
    341348         floating_heat_start       = floating_heat_end 
    … …  
    474481 
    475482 
    476    SUBROUTINE icb_dia_speed() 
    477       !!---------------------------------------------------------------------- 
    478       !!---------------------------------------------------------------------- 
    479       ! 
    480       IF( .NOT.ln_bergdia )   RETURN 
    481       nspeeding_tickets = nspeeding_tickets + 1 
     483   SUBROUTINE icb_dia_speed(pvel_factor, pn_stage) 
     484      !!---------------------------------------------------------------------- 
     485      !!---------------------------------------------------------------------- 
     486      REAL(wp), INTENT(in) ::   pvel_factor   ! factor by which velocity reduced 
     487      INTEGER , INTENT(in) ::   pn_stage  ! which stage of the RK4 calculation are we on 
     488      ! 
     489      IF( (.NOT.ln_bergdia) .OR. pn_stage .lt. 1 .OR. pn_stage .gt. 4 )   RETURN 
     490      nspeeding_tickets(pn_stage) = nspeeding_tickets(pn_stage) + 1 
     491      vel_factor_min = MIN(vel_factor_min,pvel_factor)   ! keep track of minimum reduction factor 
    482492      ! 
    483493   END SUBROUTINE icb_dia_speed 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/ICB/icbdyn.F90

    r14075 r15658  
    8585         !                                         !**   A1 = A(X1,V1) 
    8686         CALL icb_accel( berg , zxi1, ze1, zuvel1, zuvel1, zax1,     & 
    87             &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2 ) 
     87            &                   zyj1, ze2, zvvel1, zvvel1, zay1, zdt_2, 1 ) 
    8888         ! 
    8989         zu1 = zuvel1 / ze1                           !**   V1 in d(i,j)/dt 
    … …  
    102102         !                                         !**   A2 = A(X2,V2) 
    103103         CALL icb_accel( berg , zxi2, ze1, zuvel2, zuvel1, zax2,    & 
    104             &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2 ) 
     104            &                   zyj2, ze2, zvvel2, zvvel1, zay2, zdt_2, 2 ) 
    105105         ! 
    106106         zu2 = zuvel2 / ze1                           !**   V2 in d(i,j)/dt 
    … …  
    118118         !                                         !**   A3 = A(X3,V3) 
    119119         CALL icb_accel( berg , zxi3, ze1, zuvel3, zuvel1, zax3,    & 
    120             &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt ) 
     120            &                   zyj3, ze2, zvvel3, zvvel1, zay3, zdt, 3 ) 
    121121         ! 
    122122         zu3 = zuvel3 / ze1                           !**   V3 in d(i,j)/dt 
    … …  
    134134         !                                         !**   A4 = A(X4,V4) 
    135135         CALL icb_accel( berg , zxi4, ze1, zuvel4, zuvel1, zax4,    & 
    136             &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt ) 
     136            &                   zyj4, ze2, zvvel4, zvvel1, zay4, zdt, 4 ) 
    137137 
    138138         zu4 = zuvel4 / ze1                           !**   V4 in d(i,j)/dt 
    … …  
    235235 
    236236   SUBROUTINE icb_accel( berg , pxi, pe1, puvel, puvel0, pax,                & 
    237       &                         pyj, pe2, pvvel, pvvel0, pay, pdt ) 
     237      &                         pyj, pe2, pvvel, pvvel0, pay, pdt, pstage ) 
    238238      !!---------------------------------------------------------------------- 
    239239      !!                  ***  ROUTINE icb_accel  *** 
    … …  
    250250      REAL(wp)               , INTENT(inout) ::   pax, pay         ! berg acceleration 
    251251      REAL(wp)               , INTENT(in   ) ::   pdt              ! berg time step 
     252      INTEGER                , INTENT(in   ) ::   pstage           ! RK4 stage (for speed limiter) 
    252253      ! 
    253254      REAL(wp), PARAMETER ::   pp_alpha     = 0._wp      ! 
    … …  
    265266      REAL(wp) ::   zampl, zwmod, zCr, zLwavelength, zLcutoff, zLtop 
    266267      REAL(wp) ::   zlambda, zdetA, zA11, zA12, zaxe, zaye, zD_hi 
    267       REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new 
     268      REAL(wp) ::   zuveln, zvveln, zus, zvs, zspeed, zloc_dx, zspeed_new, zvel_factor, zcfl_scale 
    268269      !!---------------------------------------------------------------------- 
    269270 
    … …  
    361362         zspeed = SQRT( zuveln*zuveln + zvveln*zvveln )    ! Speed of berg 
    362363         IF( zspeed > 0._wp ) THEN 
    363             zloc_dx = MIN( pe1, pe2 )                          ! minimum grid spacing 
    364             zspeed_new = zloc_dx / pdt * rn_speed_limit        ! Speed limit as a factor of dx / dt 
     364            zloc_dx = MIN( pe1, pe2 )                                ! minimum grid spacing 
     365            ! cfl scale is function of the RK4 step 
     366            IF( pstage .le. 2 ) THEN 
     367               zcfl_scale=0.5 
     368            ELSE 
     369               zcfl_scale=1.0 
     370            ENDIF 
     371            zspeed_new = zloc_dx / pdt * rn_speed_limit * zcfl_scale ! Speed limit as a factor of dx / dt 
    365372            IF( zspeed_new < zspeed ) THEN 
    366                zuveln = zuveln * ( zspeed_new / zspeed )        ! Scale velocity to reduce speed 
    367                zvveln = zvveln * ( zspeed_new / zspeed )        ! without changing the direction 
    368                CALL icb_dia_speed() 
     373               zvel_factor = zspeed_new / zspeed 
     374               zuveln = zuveln * ( zvel_factor )             ! Scale velocity to reduce speed 
     375               zvveln = zvveln * ( zvel_factor )             ! without changing the direction 
     376               pax = (zuveln - puvel0)/pdt 
     377               pay = (zvveln - pvvel0)/pdt 
     378               WRITE(numicb,*) 'speeding ticket (zspeed_new/zspeed): ',zvel_factor, pdt, zcfl_scale 
     379               CALL FLUSH(numicb) 
     380               CALL icb_dia_speed(zvel_factor,pstage) 
    369381            ENDIF 
    370382         ENDIF 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/ICB/icbrst.F90

    r14075 r15658  
    2525   USE netcdf         ! netcdf routines for IO 
    2626   USE iom 
     27   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2728   USE icb_oce        ! define iceberg arrays 
    2829   USE icbutl         ! iceberg utility routines 
    … …  
    191192      INTEGER ::   idg  ! number of digits 
    192193      INTEGER ::   ix_dim, iy_dim, ik_dim, in_dim 
    193       CHARACTER(len=256)     :: cl_path 
    194       CHARACTER(len=256)     :: cl_filename 
    195       CHARACTER(len=8  )     :: cl_kt 
     194      INTEGER             ::   iyear, imonth, iday 
     195      REAL (wp)           ::   zsec 
     196      REAL (wp)           ::   zfjulday 
     197      CHARACTER(len=256)  :: cl_path 
     198      CHARACTER(len=256)  :: cl_filename 
     199      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    196200      CHARACTER(LEN=12 )     :: clfmt            ! writing format 
    197201      TYPE(iceberg), POINTER :: this 
    … …  
    209213         cl_path = TRIM(cn_ocerst_outdir) 
    210214         IF( cl_path(LEN_TRIM(cl_path):) /= '/' ) cl_path = TRIM(cl_path) // '/' 
    211          WRITE(cl_kt, '(i8.8)') kt 
    212          cl_filename = TRIM(cexper)//"_icebergs_"//cl_kt//"_restart" 
     215         IF ( ln_rstdate ) THEN 
     216            zfjulday = fjulday + rdt / rday 
     217            IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     218            CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     219            WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     220         ELSE 
     221            IF( kt > 999999999 ) THEN   ;   WRITE(clkt, *       ) kt 
     222            ELSE                        ;   WRITE(clkt, '(i8.8)') kt 
     223            ENDIF 
     224         ENDIF 
     225         cl_filename = TRIM(cexper)//"_icebergs_"//TRIM(ADJUSTL(clkt))//"_restart" 
    213226         IF( lk_mpp ) THEN 
    214227            idg = MAX( INT(LOG10(REAL(MAX(1,jpnij-1),wp))) + 1, 4 )          ! how many digits to we need to write? min=4, max=9 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/ICB/icbutl.F90

    r14075 r15658  
    428428      IF( ij == jpj ) THEN ; ij = ij-1 ; ierr = ierr + 1 ; END IF 
    429429      ! 
    430       IF ( ierr > 0 ) CALL ctl_stop('STOP','icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)') 
     430      IF ( ierr > 0 ) CALL ctl_stop('STOP', & 
     431     &                'icb_utl_bilin_e: an icebergs coordinates is out of valid range (out of bound error)', & 
     432     &                'This can be fixed using rn_speed_limit=0.4 in &namberg.') 
    431433      ! 
    432434      IF(    0.0_wp <= zi .AND. zi < 0.5_wp   ) THEN 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/IOM/in_out_manager.F90

    r14075 r15658  
    2828   LOGICAL       ::   ln_rstart        !: start from (F) rest or (T) a restart file 
    2929   LOGICAL       ::   ln_rst_list      !: output restarts at list of times (T) or by frequency (F) 
     30   LOGICAL       ::   ln_rst_eos       !: check equation of state used for the restart is consistent with model 
    3031   INTEGER       ::   nn_rstctl        !: control of the time step (0, 1 or 2) 
    3132   INTEGER       ::   nn_rstssh   = 0  !: hand made initilization of ssh or not (1/0) 
    … …  
    4041   INTEGER, DIMENSION(10) :: nn_stocklist  !: restart dump times 
    4142   LOGICAL       ::   ln_mskland       !: mask land points in NetCDF outputs (costly: + ~15%) 
     43   LOGICAL       ::   ln_rstdate       !: T=> stamp output restart files with date instead of timestep 
    4244   LOGICAL       ::   ln_cfmeta        !: output additional data to netCDF files required for compliance with the CF metadata standard 
    4345   LOGICAL       ::   ln_clobber       !: clobber (overwrite) an existing file 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/IOM/iom.F90

    r14075 r15658  
    368368!from restart.F90 
    369369   CALL iom_set_rstw_var_active("rdt") 
     370   CALL iom_set_rstw_var_active("neos") 
     371 
    370372   IF ( .NOT. ln_diurnal_only ) THEN 
    371373        CALL iom_set_rstw_var_active('ub'  ) 
    … …  
    417419        i = 0 
    418420        i = i + 1; fields(i)%vname="rdt";            fields(i)%grid="grid_scalar" 
     421        i = i + 1; fields(i)%vname="neos";           fields(i)%grid="grid_scalar" 
    419422        i = i + 1; fields(i)%vname="un";             fields(i)%grid="grid_N_3D" 
    420423        i = i + 1; fields(i)%vname="ub";             fields(i)%grid="grid_N_3D" 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/IOM/restart.F90

    r14075 r15658  
    2727   USE in_out_manager  ! I/O manager 
    2828   USE iom             ! I/O module 
     29   USE ioipsl, ONLY : ju2ymds    ! for calendar 
    2930   USE diurnal_bulk 
    3031   USE lib_mpp         ! distribued memory computing library 
    … …  
    5960      INTEGER, INTENT(in) ::   kt     ! ocean time-step 
    6061      !! 
     62      INTEGER             ::   iyear, imonth, iday 
     63      REAL (wp)           ::   zsec 
     64      REAL (wp)           ::   zfjulday 
    6165      CHARACTER(LEN=20)   ::   clkt     ! ocean time-step deine as a character 
    6266      CHARACTER(LEN=50)   ::   clname   ! ocean output restart file name 
    … …  
    9094         IF( nitrst <= nitend .AND. nitrst > 0 ) THEN  
    9195            ! beware of the format used to write kt (default is i8.8, that should be large enough...) 
    92             IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
    93             ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     96            IF ( ln_rstdate ) THEN 
     97               zfjulday = fjulday + rdt / rday 
     98               IF( ABS(zfjulday - REAL(NINT(zfjulday),wp)) < 0.1 / rday )   zfjulday = REAL(NINT(zfjulday),wp)   ! avoid truncation error 
     99               CALL ju2ymds( zfjulday, iyear, imonth, iday, zsec )            
     100               WRITE(clkt, '(i4.4,2i2.2)') iyear, imonth, iday 
     101            ELSE 
     102               IF( nitrst > 999999999 ) THEN   ;   WRITE(clkt, *       ) nitrst 
     103               ELSE                            ;   WRITE(clkt, '(i8.8)') nitrst 
     104               ENDIF 
    94105            ENDIF 
    95106            ! create the file 
    … …  
    173184                  END IF 
    174185      ENDIF 
     186                     CALL iom_rstput( kt, nitrst, numrow, 'neos'    , REAL(neos)      , ldxios = lwxios)   ! equation of state 
     187                     !CALL iom_rstput( kt, nitrst, numrow, 'neos'    , neos      , ktype = jp_i1, ldxios = lwxios)   ! equation of state 
     188 
    175189       
    176190      IF (ln_diurnal) CALL iom_rstput( kt, nitrst, numrow, 'Dsst', x_dsst, ldxios = lwxios )   
    … …  
    249263      !!---------------------------------------------------------------------- 
    250264      REAL(wp) ::   zrdt 
     265      REAL(wp) ::   zeos 
    251266      INTEGER  ::   jk 
    252267      REAL(wp), DIMENSION(jpi, jpj, jpk) :: w3d 
    … …  
    255270      CALL rst_read_open           ! open restart for reading (if not already opened) 
    256271 
     272      IF ( ln_rst_eos ) THEN 
     273         ! Check equation of state used is consistent with the restart 
     274         IF( iom_varid( numror, 'neos') == -1) THEN 
     275            CALL ctl_stop( 'restart, rst_read: variable neos not found. STOP check that the equations of state in the restart file and in the namelist nameos are consistent and use ln_rst_eos=F') 
     276         ELSE 
     277            CALL iom_get( numror, 'neos', zeos, ldxios = lrxios ) 
     278            IF ( INT(zeos) /= neos ) CALL ctl_stop( 'restart, rst_read: equation of state used in restart file differs from namelist nameos') 
     279         ENDIF 
     280      ENDIF 
     281 
    257282      ! Check dynamics and tracer time-step consistency and force Euler restart if changed 
    258       IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 )   THEN 
     283      IF( iom_varid( numror, 'rdt', ldstop = .FALSE. ) > 0 )   THEN  
    259284         CALL iom_get( numror, 'rdt', zrdt, ldxios = lrxios ) 
    260285         IF( zrdt /= rdt )   neuler = 0 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/LDF/ldftra.F90

    r14075 r15658  
    7272   REAL(wp), PUBLIC ::      rn_Ue               !: lateral diffusive velocity  [m/s] 
    7373   REAL(wp), PUBLIC ::      rn_Le               !: lateral diffusive length    [m] 
     74   INTEGER,  PUBLIC ::   nn_ldfeiv_shape     !: shape of bounding coefficient (Treguier et al formulation only) 
    7475    
    7576   !                                  ! Flag to control the type of lateral diffusive operator 
    … …  
    409410      !!---------------------------------------------------------------------- 
    410411      ! 
    411       IF( ln_ldfeiv .AND. nn_aei_ijk_t == 21 ) THEN       ! eddy induced velocity coefficients 
     412      IF( ln_ldfeiv .AND. ( nn_aei_ijk_t == 21 .OR. nn_aei_ijk_t == 22 ) ) THEN        
     413         !                                ! eddy induced velocity coefficients 
    412414         !                                ! =F(growth rate of baroclinic instability) 
    413415         !                                ! max value aeiv_0 ; decreased to 0 within 20N-20S 
    … …  
    502504      !! 
    503505      NAMELIST/namtra_eiv/ ln_ldfeiv   , ln_ldfeiv_dia,   &   ! eddy induced velocity (eiv) 
    504          &                 nn_aei_ijk_t, rn_Ue, rn_Le         ! eiv  coefficient 
     506         &                 nn_aei_ijk_t, rn_Ue, rn_Le,    &   ! eiv  coefficient 
     507         &                 nn_ldfeiv_shape 
    505508      !!---------------------------------------------------------------------- 
    506509      ! 
    … …  
    525528         WRITE(numout,*) '      eiv streamfunction & velocity diag.        ln_ldfeiv_dia = ', ln_ldfeiv_dia 
    526529         WRITE(numout,*) '      coefficients :' 
    527          WRITE(numout,*) '         type of time-space variation            nn_aei_ijk_t  = ', nn_aht_ijk_t 
     530         WRITE(numout,*) '         type of time-space variation            nn_aei_ijk_t  = ', nn_aei_ijk_t 
    528531         WRITE(numout,*) '         lateral diffusive velocity (if cst)     rn_Ue         = ', rn_Ue, ' m/s' 
    529532         WRITE(numout,*) '         lateral diffusive length   (if cst)     rn_Le         = ', rn_Le, ' m' 
    … …  
    595598            IF(lwp) WRITE(numout,*) '                                       = F( growth rate of baroclinic instability )' 
    596599            IF(lwp) WRITE(numout,*) '           maximum allowed value: aei0 = ', aei0, ' m2/s' 
     600            IF(lwp) WRITE(numout,*) '           shape of bounding coefficient : ',nn_ldfeiv_shape 
    597601            ! 
    598602            l_ldfeiv_time = .TRUE.     ! will be calculated by call to ldf_tra routine in step.F90 
    … …  
    638642      !!---------------------------------------------------------------------- 
    639643      INTEGER                         , INTENT(in   ) ::   kt             ! ocean time-step index 
    640       REAL(wp)                        , INTENT(inout) ::   paei0          ! max value            [m2/s] 
     644      REAL(wp)                        , INTENT(in   ) ::   paei0          ! max value            [m2/s] 
    641645      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) ::   paeiu, paeiv   ! eiv coefficient      [m2/s] 
    642646      ! 
    643647      INTEGER  ::   ji, jj, jk    ! dummy loop indices 
    644       REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei    ! local scalars 
    645       REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zaeiw   ! 2D workspace 
     648      REAL(wp) ::   zfw, ze3w, zn2, z1_f20, zaht, zaht_min, zzaei, z2_3    ! local scalars 
     649      REAL(wp), DIMENSION(jpi,jpj) ::   zn, zah, zhw, zRo, zRo_lim, zTclinic_recip, zaeiw, zratio   ! 2D workspace 
     650      REAL(wp), DIMENSION(jpi,jpj,jpk) :: zmodslp ! 3D workspace 
    646651      !!---------------------------------------------------------------------- 
    647652      ! 
    648653      zn (:,:) = 0._wp        ! Local initialization 
     654      zmodslp(:,:,:) = 0._wp 
    649655      zhw(:,:) = 5._wp 
    650656      zah(:,:) = 0._wp 
    651657      zRo(:,:) = 0._wp 
     658      zRo_lim(:,:) = 0._wp 
     659      zTclinic_recip(:,:) = 0._wp 
     660      zratio(:,:) = 0._wp  
     661      zaeiw(:,:) = 0._wp     
    652662      !                       ! Compute lateral diffusive coefficient at T-point 
    653663      IF( ln_traldf_triad ) THEN 
    … …  
    682692                  ! T^-1 = sqrt(m_jpk(N^2*(r1^2+r2^2)*e3w)) 
    683693                  ze3w = e3w_n(ji,jj,jk) * wmask(ji,jj,jk) 
    684                   zah(ji,jj) = zah(ji,jj) + zn2 * ( wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
    685                      &                            + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) ) * ze3w 
     694                  zmodslp(ji,jj,jk) =  wslpi(ji,jj,jk) * wslpi(ji,jj,jk)   & 
     695                     &               + wslpj(ji,jj,jk) * wslpj(ji,jj,jk) 
     696                  zah(ji,jj) = zah(ji,jj) + zn2 * zmodslp(ji,jj,jk) * ze3w 
    686697                  zhw(ji,jj) = zhw(ji,jj) + ze3w 
    687698               END DO 
    … …  
    693704         DO ji = fs_2, fs_jpim1   ! vector opt. 
    694705            zfw = MAX( ABS( 2. * omega * SIN( rad * gphit(ji,jj) ) ) , 1.e-10 ) 
    695             ! Rossby radius at w-point taken betwenn 2 km and  40km 
    696             zRo(ji,jj) = MAX(  2.e3 , MIN( .4 * zn(ji,jj) / zfw, 40.e3 )  ) 
     706            ! Rossby radius at w-point taken between 2 km and  40km 
     707            zRo(ji,jj) = .4 * zn(ji,jj) / zfw 
     708            zRo_lim(ji,jj) = MAX(  2.e3 , MIN( zRo(ji,jj), 40.e3 )  ) 
    697709            ! Compute aeiw by multiplying Ro^2 and T^-1 
    698             zaeiw(ji,jj) = zRo(ji,jj) * zRo(ji,jj) * SQRT( zah(ji,jj) / zhw(ji,jj) ) * tmask(ji,jj,1) 
     710            zTclinic_recip(ji,jj) = SQRT( MAX(zah(ji,jj),0._wp) / zhw(ji,jj) ) * tmask(ji,jj,1) 
     711            zaeiw(ji,jj) = zRo_lim(ji,jj) * zRo_lim(ji,jj) * zTclinic_recip(ji,jj)  
    699712         END DO 
    700713      END DO 
    701  
     714      IF( iom_use('N_2d') ) CALL iom_put('N_2d',zn(:,:)/zhw(:,:)) 
     715      IF( iom_use('modslp') ) CALL iom_put('modslp',SQRT(zmodslp(:,:,:)) ) 
     716      CALL iom_put('RossRad',zRo) 
     717      CALL iom_put('RossRadlim',zRo_lim) 
     718      CALL iom_put('Tclinic_recip',zTclinic_recip) 
    702719      !                                         !==  Bound on eiv coeff.  ==! 
    703720      z1_f20 = 1._wp / (  2._wp * omega * sin( rad * 20._wp )  ) 
    704       DO jj = 2, jpjm1 
    705          DO ji = fs_2, fs_jpim1   ! vector opt. 
    706             zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease 
    707             zaeiw(ji,jj) = MIN( zzaei , paei0 )                                  ! Max value = paei0 
    708          END DO 
    709       END DO 
     721      z2_3 = 2._wp/3._wp 
     722 
     723      SELECT CASE(nn_ldfeiv_shape) 
     724         CASE(0) !! Standard shape applied - decrease in tropics and cap.  
     725            DO jj = 2, jpjm1 
     726               DO ji = fs_2, fs_jpim1   ! vector opt. 
     727                  zzaei = MIN( 1._wp, ABS( ff_t(ji,jj) * z1_f20 ) ) * zaeiw(ji,jj)     ! tropical decrease 
     728                  zaeiw(ji,jj) = MIN( zzaei, paei0 ) 
     729               END DO 
     730            END DO 
     731 
     732         CASE(1) !! Abrupt cut-off on Rossby radius: 
     733! JD : modifications here to introduce scaling by local rossby radius of deformation vs local grid scale 
     734! arbitrary decision that GM is de-activated if local rossy radius larger than 2 times local grid scale 
     735! based on Hallberg (2013) 
     736            DO jj = 2, jpjm1 
     737               DO ji = fs_2, fs_jpim1   ! vector opt. 
     738                  IF ( zRo(ji,jj) >= ( 2._wp * MIN( e1t(ji,jj), e2t(ji,jj) ) ) ) THEN 
     739! TODO : use a version of zRo that integrates over a few time steps ? 
     740                      zaeiw(ji,jj) = 0._wp 
     741                  ELSE 
     742                      zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 ) 
     743                  ENDIF 
     744               END DO 
     745            END DO 
     746 
     747         CASE(2) !! Rossby radius ramp type 1: 
     748            DO jj = 2, jpjm1 
     749               DO ji = fs_2, fs_jpim1   ! vector opt. 
     750                  zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj)) 
     751                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*(2._wp - zratio(ji,jj)) ) ) * paei0 ) 
     752               END DO 
     753            END DO 
     754            CALL iom_put('RR_GS',zratio) 
     755 
     756         CASE(3) !! Rossby radius ramp type 2: 
     757            DO jj = 2, jpjm1 
     758               DO ji = fs_2, fs_jpim1   ! vector opt. 
     759                  zratio(ji,jj) = MIN(e1t(ji,jj),e2t(ji,jj))/zRo(ji,jj) 
     760                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, z2_3*( zratio(ji,jj) - 0.5_wp ) ) ) * paei0 ) 
     761               END DO 
     762            END DO 
     763 
     764         CASE(4) !! bathymetry ramp: 
     765            DO jj = 2, jpjm1 
     766               DO ji = fs_2, fs_jpim1   ! vector opt. 
     767                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), MAX( 0._wp, MIN( 1._wp, 0.001*(ht_0(ji,jj) - 2000._wp) ) ) * paei0 ) 
     768               END DO 
     769            END DO 
     770 
     771         CASE(5) !! Rossby radius ramp type 1 applied to Treguier et al coefficient rather than cap: 
     772                 !! Note the ramp is RR/GS=[2.0,1.0] (not [2.0,0.5] as for cases 2,3) and we ramp up  
     773                 !! to 5% of the Treguier et al coefficient, aiming for peak values of around 100m2/s 
     774                 !! at high latitudes rather than 2000m2/s which is what you get in eORCA025 with an  
     775                 !! uncapped coefficient. 
     776            DO jj = 2, jpjm1 
     777               DO ji = fs_2, fs_jpim1   ! vector opt. 
     778                  zratio(ji,jj) = zRo(ji,jj)/MIN(e1t(ji,jj),e2t(ji,jj)) 
     779                  zaeiw(ji,jj) = MAX( 0._wp, MIN( 1._wp, 2._wp - zratio(ji,jj) ) ) * 0.05 * zaeiw(ji,jj) 
     780                  zaeiw(ji,jj) = MIN( zaeiw(ji,jj), paei0 ) 
     781               END DO 
     782            END DO 
     783            CALL iom_put('RR_GS',zratio) 
     784 
     785         CASE DEFAULT 
     786               CALL ctl_stop('ldf_eiv: Unrecognised option for nn_ldfeiv_shape.')          
     787 
     788      END SELECT 
     789 
    710790      CALL lbc_lnk( 'ldftra', zaeiw(:,:), 'W', 1. )       ! lateral boundary condition 
    711791      !                
    … …  
    722802         paeiv(:,:,jk) = paeiv(:,:,1) * vmask(:,:,jk) 
    723803      END DO 
    724       !   
     804      ! 
    725805   END SUBROUTINE ldf_eiv 
    726806 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/SBC/sbcisf.F90

    r14075 r15658  
    355355                ik = 2 
    356356!!gm potential bug: use gdepw_0 not _n 
    357                 DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_n(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
     357!!DS following Gurvan's suggestion here because it fixes a restartability issue. 
     358                DO WHILE ( ik <= mbkt(ji,jj) .AND. gdepw_0(ji,jj,ik) < rzisf_tbl(ji,jj) ) ;  ik = ik + 1 ;  END DO 
    358359                misfkt(ji,jj) = ik-1 
    359360            END DO 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/SBC/sbcrnf.F90

    r14075 r15658  
    1818   USE dom_oce        ! ocean space and time domain 
    1919   USE phycst         ! physical constants 
     20   USE sbcisf         ! surface boundary condition: ice-shelf 
    2021   USE sbc_oce        ! surface boundary condition variables 
    2122   USE eosbn2         ! Equation Of State 
    … …  
    140141               rnf_tsc(:,:,jp_tem) = sst_m(:,:) * rnf(:,:) * r1_rau0 
    141142            END WHERE 
     143            WHERE( sf_t_rnf(1)%fnow(:,:,1) == -222._wp ) 
     144            ! where fwf comes from melting of ice shelves or iceberg 
     145            rnf_tsc(:,:,jp_tem) = ztfrz(:,:) * rnf(:,:) * r1_rau0 - rnf(:,:) * rLfusisf * r1_rau0_rcp 
     146            END WHERE 
    142147         ELSE                                                        ! use SST as runoffs temperature 
    143148            !CEOD River is fresh water so must at least be 0 unless we consider ice 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/SBC/sbcssm.F90

    r14075 r15658  
    5757      REAL(wp) ::   zcoef, zf_sbc       ! local scalar 
    5858      REAL(wp), DIMENSION(jpi,jpj,jpts) :: zts 
     59      CHARACTER(len=4),SAVE :: stype 
    5960      !!--------------------------------------------------------------------- 
     61      IF( kt == nit000 ) THEN 
     62         IF( ln_TEOS10 ) THEN 
     63            stype='abs'   ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module) 
     64         ELSE IF( ln_EOS80  ) THEN 
     65            stype='pra'   ! eos-80: using practical salinity 
     66         ELSE IF ( ln_SEOS) THEN 
     67            stype='seos' ! seos using Simplified Equation of state (sst is converted to potential temperature for the surface module) 
     68         ENDIF 
     69      ENDIF 
    6070      ! 
    6171      !                                        !* surface T-, U-, V- ocean level variables (T, S, depth, velocity) 
    … …  
    174184         CALL iom_put( 'ssu_m', ssu_m ) 
    175185         CALL iom_put( 'ssv_m', ssv_m ) 
    176          CALL iom_put( 'sst_m', sst_m ) 
    177          CALL iom_put( 'sss_m', sss_m ) 
     186         CALL iom_put( 'sst_m_pot', sst_m ) 
     187         CALL iom_put( 'sss_m_'//stype, sss_m ) 
    178188         CALL iom_put( 'ssh_m', ssh_m ) 
    179189         CALL iom_put( 'e3t_m', e3t_m ) 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/ZDF/zdfmxl.F90

    r14075 r15658  
    1515   USE trc_oce  , ONLY: l_offline         ! ocean space and time domain variables 
    1616   USE zdf_oce        ! ocean vertical physics 
     17   USE eosbn2         ! for zdf_mxl_zint 
    1718   ! 
    1819   USE in_out_manager ! I/O manager 
    … …  
    3132   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlp    !: mixed layer depth  (rho=rho0+zdcrit) [m]   (used by LDF) 
    3233   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hmlpt   !: depth of the last T-point inside the mixed layer [m] (used by LDF) 
     34   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   hmld_zint  !: vertically-interpolated mixed layer depth   [m] 
     35   REAL(wp), PUBLIC, ALLOCATABLE,       DIMENSION(:,:) ::   htc_mld    ! Heat content of hmld_zint 
     36   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)    :: ll_found   ! Is T_b to be found by interpolation ? 
     37   LOGICAL, PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:)  :: ll_belowml ! Flag points below mixed layer when ll_found=F 
    3338 
    3439   REAL(wp), PUBLIC ::   rho_c = 0.01_wp    !: density criterion for mixed layer depth 
    3540   REAL(wp), PUBLIC ::   avt_c = 5.e-4_wp   ! Kz criterion for the turbocline depth 
     41 
     42   TYPE, PUBLIC :: MXL_ZINT   !: Structure for MLD defs 
     43      INTEGER   :: mld_type   ! mixed layer type      
     44      REAL(wp)  :: zref       ! depth of initial T_ref 
     45      REAL(wp)  :: dT_crit    ! Critical temp diff 
     46      REAL(wp)  :: iso_frac   ! Fraction of rn_dT_crit  
     47   END TYPE MXL_ZINT 
    3648 
    3749   !!---------------------------------------------------------------------- 
    … …  
    4860      zdf_mxl_alloc = 0      ! set to zero if no array to be allocated 
    4961      IF( .NOT. ALLOCATED( nmln ) ) THEN 
    50          ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), STAT= zdf_mxl_alloc ) 
     62         ALLOCATE( nmln(jpi,jpj), hmld(jpi,jpj), hmlp(jpi,jpj), hmlpt(jpi,jpj), hmld_zint(jpi,jpj),     & 
     63   &          htc_mld(jpi,jpj), ll_found(jpi,jpj), ll_belowml(jpi,jpj,jpk), STAT= zdf_mxl_alloc )          
    5164         ! 
    5265         CALL mpp_sum ( 'zdfmxl', zdf_mxl_alloc ) 
    … …  
    137150      ENDIF 
    138151      ! 
     152      ! Vertically-interpolated mixed-layer depth diagnostic 
     153      CALL zdf_mxl_zint( kt ) 
     154      ! 
    139155      IF(ln_ctl)   CALL prt_ctl( tab2d_1=REAL(nmln,wp), clinfo1=' nmln : ', tab2d_2=hmlp, clinfo2=' hmlp : ' ) 
    140156      ! 
    141157   END SUBROUTINE zdf_mxl 
     158 
     159   SUBROUTINE zdf_mxl_zint_mld( sf )  
     160      !!----------------------------------------------------------------------------------  
     161      !!                    ***  ROUTINE zdf_mxl_zint_mld  ***  
     162      !                                                                         
     163      !   Calculate vertically-interpolated mixed layer depth diagnostic.  
     164      !             
     165      !   This routine can calculate the mixed layer depth diagnostic suggested by 
     166      !   Kara et al, 2000, JGR, 105, 16803, but is more general and can calculate 
     167      !   vertically-interpolated mixed-layer depth diagnostics with other parameter 
     168      !   settings set in the namzdf_mldzint namelist.   
     169      !  
     170      !   If mld_type=1 the mixed layer depth is calculated as the depth at which the   
     171      !   density has increased by an amount equivalent to a temperature difference of   
     172      !   0.8C at the surface.  
     173      !  
     174      !   For other values of mld_type the mixed layer is calculated as the depth at   
     175      !   which the temperature differs by 0.8C from the surface temperature.   
     176      !                                                                         
     177      !   David Acreman, Daley Calvert                                       
     178      !  
     179      !!-----------------------------------------------------------------------------------  
     180 
     181      TYPE(MXL_ZINT), INTENT(in)  :: sf 
     182 
     183      ! Diagnostic criteria 
     184      INTEGER   :: nn_mld_type   ! mixed layer type      
     185      REAL(wp)  :: rn_zref       ! depth of initial T_ref 
     186      REAL(wp)  :: rn_dT_crit    ! Critical temp diff 
     187      REAL(wp)  :: rn_iso_frac   ! Fraction of rn_dT_crit used 
     188 
     189      ! Local variables 
     190      REAL(wp), PARAMETER :: zepsilon = 1.e-30          ! local small value 
     191      INTEGER, DIMENSION(jpi,jpj) :: ikmt          ! number of active tracer levels  
     192      INTEGER, DIMENSION(jpi,jpj) :: ik_ref        ! index of reference level  
     193      INTEGER, DIMENSION(jpi,jpj) :: ik_iso        ! index of last uniform temp level  
     194      REAL, DIMENSION(jpi,jpj,jpk)  :: zT            ! Temperature or density  
     195      REAL, DIMENSION(jpi,jpj)    :: ppzdep        ! depth for use in calculating d(rho)  
     196      REAL, DIMENSION(jpi,jpj)    :: zT_ref        ! reference temperature  
     197      REAL    :: zT_b                                   ! base temperature  
     198      REAL, DIMENSION(jpi,jpj,jpk)  :: zdTdz         ! gradient of zT  
     199      REAL, DIMENSION(jpi,jpj,jpk)  :: zmoddT        ! Absolute temperature difference  
     200      REAL    :: zdz                                    ! depth difference  
     201      REAL    :: zdT                                    ! temperature difference  
     202      REAL, DIMENSION(jpi,jpj)    :: zdelta_T      ! difference critereon  
     203      REAL, DIMENSION(jpi,jpj)    :: zRHO1, zRHO2  ! Densities  
     204      INTEGER :: ji, jj, jk                             ! loop counter  
     205 
     206      !!-------------------------------------------------------------------------------------  
     207      !   
     208      ! Unpack structure 
     209      nn_mld_type = sf%mld_type 
     210      rn_zref     = sf%zref 
     211      rn_dT_crit  = sf%dT_crit 
     212      rn_iso_frac = sf%iso_frac 
     213 
     214      ! Set the mixed layer depth criterion at each grid point  
     215      IF( nn_mld_type == 0 ) THEN 
     216         zdelta_T(:,:) = rn_dT_crit 
     217         zT(:,:,:) = rhop(:,:,:) 
     218      ELSE IF( nn_mld_type == 1 ) THEN 
     219         ppzdep(:,:)=0.0  
     220         call eos ( tsn(:,:,1,:), ppzdep(:,:), zRHO1(:,:) )  
     221! Use zT temporarily as a copy of tsn with rn_dT_crit added to SST  
     222! [assumes number of tracers less than number of vertical levels]  
     223         zT(:,:,1:jpts)=tsn(:,:,1,1:jpts)  
     224         zT(:,:,jp_tem)=zT(:,:,1)+rn_dT_crit  
     225         CALL eos( zT(:,:,1:jpts), ppzdep(:,:), zRHO2(:,:) )  
     226         zdelta_T(:,:) = abs( zRHO1(:,:) - zRHO2(:,:) ) * rau0  
     227         ! RHO from eos (2d version) doesn't calculate north or east halo:  
     228         CALL lbc_lnk( 'zdfmxl', zdelta_T, 'T', 1. )  
     229         zT(:,:,:) = rhop(:,:,:)  
     230      ELSE  
     231         zdelta_T(:,:) = rn_dT_crit                       
     232         zT(:,:,:) = tsn(:,:,:,jp_tem)                            
     233      END IF  
     234 
     235      ! Calculate the gradient of zT and absolute difference for use later  
     236      DO jk = 1 ,jpk-2  
     237         zdTdz(:,:,jk)  =    ( zT(:,:,jk+1) - zT(:,:,jk) ) / e3w_n(:,:,jk+1)  
     238         zmoddT(:,:,jk) = abs( zT(:,:,jk+1) - zT(:,:,jk) )  
     239      END DO  
     240 
     241      ! Find density/temperature at the reference level (Kara et al use 10m).           
     242      ! ik_ref is the index of the box centre immediately above or at the reference level  
     243      ! Find rn_zref in the array of model level depths and find the ref     
     244      ! density/temperature by linear interpolation.                                    
     245      DO jk = jpkm1, 2, -1  
     246         WHERE ( gdept_n(:,:,jk) > rn_zref )  
     247           ik_ref(:,:) = jk - 1  
     248           zT_ref(:,:) = zT(:,:,jk-1) + zdTdz(:,:,jk-1) * ( rn_zref - gdept_n(:,:,jk-1) )  
     249         END WHERE  
     250      END DO  
     251 
     252      ! If the first grid box centre is below the reference level then use the  
     253      ! top model level to get zT_ref  
     254      WHERE ( gdept_n(:,:,1) > rn_zref )   
     255         zT_ref = zT(:,:,1)  
     256         ik_ref = 1  
     257      END WHERE  
     258 
     259      ! The number of active tracer levels is 1 less than the number of active w levels  
     260      ikmt(:,:) = mbkt(:,:) - 1  
     261 
     262      ! Initialize / reset 
     263      ll_found(:,:) = .false. 
     264 
     265      IF ( rn_iso_frac - zepsilon > 0. ) THEN 
     266         ! Search for a uniform density/temperature region where adjacent levels           
     267         ! differ by less than rn_iso_frac * deltaT.                                       
     268         ! ik_iso is the index of the last level in the uniform layer   
     269         ! ll_found indicates whether the mixed layer depth can be found by interpolation  
     270         ik_iso(:,:)   = ik_ref(:,:)  
     271         DO jj = 1, nlcj  
     272            DO ji = 1, nlci  
     273!CDIR NOVECTOR  
     274               DO jk = ik_ref(ji,jj), ikmt(ji,jj)-1  
     275                  IF ( zmoddT(ji,jj,jk) > ( rn_iso_frac * zdelta_T(ji,jj) ) ) THEN  
     276                     ik_iso(ji,jj)   = jk  
     277                     ll_found(ji,jj) = ( zmoddT(ji,jj,jk) > zdelta_T(ji,jj) )  
     278                     EXIT  
     279                  END IF  
     280               END DO  
     281            END DO  
     282         END DO  
     283 
     284         ! Use linear interpolation to find depth of mixed layer base where possible  
     285         hmld_zint(:,:) = rn_zref  
     286         DO jj = 1, jpj  
     287            DO ji = 1, jpi  
     288               IF (ll_found(ji,jj) .and. tmask(ji,jj,1) == 1.0) THEN  
     289                  zdz =  abs( zdelta_T(ji,jj) / zdTdz(ji,jj,ik_iso(ji,jj)) )  
     290                  hmld_zint(ji,jj) = gdept_n(ji,jj,ik_iso(ji,jj)) + zdz  
     291               END IF  
     292            END DO  
     293         END DO  
     294      END IF 
     295 
     296      ! If ll_found = .false. then calculate MLD using difference of zdelta_T     
     297      ! from the reference density/temperature  
     298  
     299! Prevent this section from working on land points  
     300      WHERE ( tmask(:,:,1) /= 1.0 )  
     301         ll_found = .true.  
     302      END WHERE  
     303  
     304      DO jk=1, jpk  
     305         ll_belowml(:,:,jk) = abs( zT(:,:,jk) - zT_ref(:,:) ) >= zdelta_T(:,:)   
     306      END DO  
     307  
     308! Set default value where interpolation cannot be used (ll_found=false)   
     309      DO jj = 1, jpj  
     310         DO ji = 1, jpi  
     311            IF ( .not. ll_found(ji,jj) )  hmld_zint(ji,jj) = gdept_n(ji,jj,ikmt(ji,jj))  
     312         END DO  
     313      END DO  
     314 
     315      DO jj = 1, jpj  
     316         DO ji = 1, jpi  
     317!CDIR NOVECTOR  
     318            DO jk = ik_ref(ji,jj)+1, ikmt(ji,jj)  
     319               IF ( ll_found(ji,jj) ) EXIT  
     320               IF ( ll_belowml(ji,jj,jk) ) THEN                 
     321                  zT_b = zT_ref(ji,jj) + zdelta_T(ji,jj) * SIGN(1.0, zdTdz(ji,jj,jk-1) )  
     322                  zdT  = zT_b - zT(ji,jj,jk-1)                                       
     323                  zdz  = zdT / zdTdz(ji,jj,jk-1)                                        
     324                  hmld_zint(ji,jj) = gdept_n(ji,jj,jk-1) + zdz  
     325                  EXIT                                                    
     326               END IF  
     327            END DO  
     328         END DO  
     329      END DO  
     330 
     331      hmld_zint(:,:) = hmld_zint(:,:)*tmask(:,:,1)  
     332      !   
     333   END SUBROUTINE zdf_mxl_zint_mld 
     334 
     335   SUBROUTINE zdf_mxl_zint_htc( kt ) 
     336      !!---------------------------------------------------------------------- 
     337      !!                  ***  ROUTINE zdf_mxl_zint_htc  *** 
     338      !!  
     339      !! ** Purpose :    
     340      !! 
     341      !! ** Method  :    
     342      !!---------------------------------------------------------------------- 
     343 
     344      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     345 
     346      INTEGER :: ji, jj, jk 
     347      INTEGER :: ikmax 
     348      REAL(wp) :: zc, zcoef 
     349      ! 
     350      INTEGER,  ALLOCATABLE, DIMENSION(:,:) ::   ilevel 
     351      REAL(wp), ALLOCATABLE, DIMENSION(:,:) ::   zthick_0, zthick 
     352 
     353      !!---------------------------------------------------------------------- 
     354 
     355      IF( .NOT. ALLOCATED(ilevel) ) THEN 
     356         ALLOCATE( ilevel(jpi,jpj), zthick_0(jpi,jpj), & 
     357         &         zthick(jpi,jpj), STAT=ji ) 
     358         IF( lk_mpp  )   CALL mpp_sum( 'zdfmxl', ji ) 
     359         IF( ji /= 0 )   CALL ctl_stop( 'STOP', 'zdf_mxl_zint_htc : unable to allocate arrays' ) 
     360      ENDIF 
     361 
     362      ! Find last whole model T level above the MLD 
     363      ilevel(:,:)   = 0 
     364      zthick_0(:,:) = 0._wp 
     365 
     366      DO jk = 1, jpkm1   
     367         DO jj = 1, jpj 
     368            DO ji = 1, jpi                     
     369               zthick_0(ji,jj) = zthick_0(ji,jj) + e3t_n(ji,jj,jk) 
     370               IF( zthick_0(ji,jj) < hmld_zint(ji,jj) )   ilevel(ji,jj) = jk 
     371            END DO 
     372         END DO 
     373         WRITE(numout,*) 'zthick_0(jk =',jk,') =',zthick_0(2,2) 
     374         WRITE(numout,*) 'gdepw_n(jk+1 =',jk+1,') =',gdepw_n(2,2,jk+1) 
     375      END DO 
     376 
     377      ! Surface boundary condition 
     378      IF( ln_linssh ) THEN  ;   zthick(:,:) = sshn(:,:)   ;   htc_mld(:,:) = tsn(:,:,1,jp_tem) * sshn(:,:) * tmask(:,:,1)    
     379      ELSE                  ;   zthick(:,:) = 0._wp       ;   htc_mld(:,:) = 0._wp                                    
     380      ENDIF 
     381 
     382      ! Deepest whole T level above the MLD 
     383      ikmax = MIN( MAXVAL( ilevel(:,:) ), jpkm1 ) 
     384 
     385      ! Integration down to last whole model T level 
     386      DO jk = 1, ikmax 
     387         DO jj = 1, jpj 
     388            DO ji = 1, jpi 
     389               zc = e3t_n(ji,jj,jk) * REAL( MIN( MAX( 0, ilevel(ji,jj) - jk + 1 ) , 1  )  )    ! 0 below ilevel 
     390               zthick(ji,jj) = zthick(ji,jj) + zc 
     391               htc_mld(ji,jj) = htc_mld(ji,jj) + zc * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk) 
     392            END DO 
     393         END DO 
     394      END DO 
     395 
     396      ! Subsequent partial T level 
     397      zthick(:,:) = hmld_zint(:,:) - zthick(:,:)   !   remaining thickness to reach MLD 
     398 
     399      DO jj = 1, jpj 
     400         DO ji = 1, jpi 
     401            htc_mld(ji,jj) = htc_mld(ji,jj) + tsn(ji,jj,ilevel(ji,jj)+1,jp_tem)  &  
     402      &                      * MIN( e3t_n(ji,jj,ilevel(ji,jj)+1), zthick(ji,jj) ) * tmask(ji,jj,ilevel(ji,jj)+1) 
     403         END DO 
     404      END DO 
     405 
     406      WRITE(numout,*) 'htc_mld(after) =',htc_mld(2,2) 
     407 
     408      ! Convert to heat content 
     409      zcoef = rau0 * rcp 
     410      htc_mld(:,:) = zcoef * htc_mld(:,:) 
     411 
     412   END SUBROUTINE zdf_mxl_zint_htc 
     413 
     414   SUBROUTINE zdf_mxl_zint( kt ) 
     415      !!---------------------------------------------------------------------- 
     416      !!                  ***  ROUTINE zdf_mxl_zint  *** 
     417      !!  
     418      !! ** Purpose :    
     419      !! 
     420      !! ** Method  :    
     421      !!---------------------------------------------------------------------- 
     422 
     423      INTEGER, INTENT(in) ::   kt   ! ocean time-step index 
     424 
     425      INTEGER :: ios 
     426      INTEGER :: jn 
     427 
     428      INTEGER :: nn_mld_diag = 0    ! number of diagnostics 
     429 
     430      CHARACTER(len=1) :: cmld 
     431 
     432      TYPE(MXL_ZINT) :: sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     433      TYPE(MXL_ZINT), SAVE, DIMENSION(5) ::   mld_diags 
     434 
     435      NAMELIST/namzdf_mldzint/ nn_mld_diag, sn_mld1, sn_mld2, sn_mld3, sn_mld4, sn_mld5 
     436 
     437      !!---------------------------------------------------------------------- 
     438       
     439      IF( kt == nit000 ) THEN 
     440         REWIND( numnam_ref )              ! Namelist namzdf_mldzint in reference namelist  
     441         READ  ( numnam_ref, namzdf_mldzint, IOSTAT = ios, ERR = 901) 
     442901      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in reference namelist' ) 
     443 
     444         REWIND( numnam_cfg )              ! Namelist namzdf_mldzint in configuration namelist  
     445         READ  ( numnam_cfg, namzdf_mldzint, IOSTAT = ios, ERR = 902 ) 
     446902      IF( ios /= 0 ) CALL ctl_nam ( ios , 'namzdf_mldzint in configuration namelist' ) 
     447         IF(lwm) WRITE ( numond, namzdf_mldzint ) 
     448 
     449         IF( nn_mld_diag > 5 )   CALL ctl_stop( 'STOP', 'zdf_mxl_ini: Specify no more than 5 MLD definitions' ) 
     450 
     451         mld_diags(1) = sn_mld1 
     452         mld_diags(2) = sn_mld2 
     453         mld_diags(3) = sn_mld3 
     454         mld_diags(4) = sn_mld4 
     455         mld_diags(5) = sn_mld5 
     456 
     457         IF( lwp .AND. (nn_mld_diag > 0) ) THEN 
     458            WRITE(numout,*) '=============== Vertically-interpolated mixed layer ================' 
     459            WRITE(numout,*) '(Diagnostic number, nn_mld_type, rn_zref, rn_dT_crit, rn_iso_frac)' 
     460            DO jn = 1, nn_mld_diag 
     461               WRITE(numout,*) 'MLD criterion',jn,':' 
     462               WRITE(numout,*) '    nn_mld_type =', mld_diags(jn)%mld_type 
     463               WRITE(numout,*) '    rn_zref ='    , mld_diags(jn)%zref 
     464               WRITE(numout,*) '    rn_dT_crit =' , mld_diags(jn)%dT_crit 
     465               WRITE(numout,*) '    rn_iso_frac =', mld_diags(jn)%iso_frac 
     466            END DO 
     467            WRITE(numout,*) '====================================================================' 
     468         ENDIF 
     469      ENDIF 
     470 
     471      IF( nn_mld_diag > 0 ) THEN 
     472         DO jn = 1, nn_mld_diag 
     473            WRITE(cmld,'(I1)') jn 
     474            IF( iom_use( "mldzint_"//cmld ) .OR. iom_use( "mldhtc_"//cmld ) ) THEN 
     475               CALL zdf_mxl_zint_mld( mld_diags(jn) ) 
     476 
     477               IF( iom_use( "mldzint_"//cmld ) ) THEN 
     478                  CALL iom_put( "mldzint_"//cmld, hmld_zint(:,:) ) 
     479               ENDIF 
     480 
     481               IF( iom_use( "mldhtc_"//cmld ) )  THEN 
     482                  CALL zdf_mxl_zint_htc( kt ) 
     483                  CALL iom_put( "mldhtc_"//cmld , htc_mld(:,:)   ) 
     484               ENDIF 
     485            ENDIF 
     486         END DO 
     487      ENDIF 
     488 
     489   END SUBROUTINE zdf_mxl_zint 
    142490 
    143491   !!====================================================================== 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/step.F90

    r14075 r15658  
    207207      IF( ln_diaptr  )   CALL dia_ptr                 ! Poleward adv/ldf TRansports diagnostics 
    208208      IF( ln_diaharm )   CALL dia_harm( kstp )        ! Tidal harmonic analysis 
     209                         CALL dia_prod( kstp )        ! ocean model: product diagnostics 
    209210                         CALL dia_wri ( kstp )        ! ocean model: outputs 
    210211      ! 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/OCE/step_oce.F90

    r14075 r15658  
    8080   USE diahsb          ! heat, salt and volume budgets    (dia_hsb routine) 
    8181   USE diaharm 
     82   USE diaprod 
    8283   USE diacfl 
    8384   USE diaobs          ! Observation operator 

  • NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/SAS/sbcssm.F90

    r14075 r15658  
    7777      REAL(wp) ::   ztinta     ! ratio applied to after  records when doing time interpolation 
    7878      REAL(wp) ::   ztintb     ! ratio applied to before records when doing time interpolation 
    79       !!---------------------------------------------------------------------- 
     79      CHARACTER(len=4),SAVE :: stype 
     80      !!--------------------------------------------------------------------- 
     81      IF( kt == nit000 ) THEN 
     82         IF( ln_TEOS10 ) THEN 
     83            stype='abs'   ! teos-10: using absolute salinity (sst is converted to potential temperature for the surface module) 
     84         ELSE IF( ln_EOS80  ) THEN 
     85            stype='pra'   ! eos-80: using practical salinity 
     86         ELSE IF ( ln_SEOS) THEN 
     87            stype='seos' ! seos using Simplified Equation of state (sst is converted to potential temperature for the surface module) 
     88         ENDIF 
     89      ENDIF 
    8090      ! 
    8191      IF( ln_timing )   CALL timing_start( 'sbc_ssm') 
    … …  
    144154         CALL iom_put( 'ssu_m', ssu_m ) 
    145155         CALL iom_put( 'ssv_m', ssv_m ) 
    146          CALL iom_put( 'sst_m', sst_m ) 
    147          CALL iom_put( 'sss_m', sss_m ) 
     156         CALL iom_put( 'sst_m_pot', sst_m ) 
     157         CALL iom_put( 'sss_m_'//stype, sss_m ) 
    148158         CALL iom_put( 'ssh_m', ssh_m ) 
    149159         IF( .NOT.ln_linssh )   CALL iom_put( 'e3t_m', e3t_m ) 
Note: See TracChangeset for help on using the changeset viewer.