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 4161 – NEMO

Changeset 4161


Ignore:
Timestamp:
2013-11-07T11:01:27+01:00 (11 years ago)
Author:
cetlod
Message:

dev_LOCEAN_2013 : merge in the 3rd dev branch dev_r4028_CNRS_LIM3, see ticket #1169

Location:
branches/2013/dev_LOCEAN_2013/NEMOGCM
Files:
9 added
8 deleted
49 edited
6 copied

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/ORCA2_LIM/EXP00/iodef_default.xml

    r4153 r4161  
    2626 
    2727      <file_group id="1h" output_freq="1h"  output_level="10" enabled=".TRUE."/> <!-- 1h files --> 
     28 
    2829      <file_group id="2h" output_freq="2h"  output_level="10" enabled=".TRUE."/> <!-- 2h files --> 
     30 
    2931      <file_group id="3h" output_freq="3h"  output_level="10" enabled=".TRUE."/> <!-- 3h files -->      
     32 
    3033      <file_group id="4h" output_freq="4h"  output_level="10" enabled=".TRUE."/> <!-- 4h files --> 
    31       <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files --> 
    32       
    33       <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE."> <!-- 1d files --> 
     34 
     35      <file_group id="6h" output_freq="6h"  output_level="10" enabled=".TRUE."/> <!-- 6h files -->      
     36 
     37      <file_group id="1d" output_freq="1d"  output_level="10" enabled=".TRUE." > <!-- 1d files --> 
    3438 
    3539   <file id="file1" name_suffix="_grid_T" description="ocean T grid variables" > 
     
    3741     <field field_ref="sss"          name="sos"      long_name="sea_surface_salinity"                          /> 
    3842     <field field_ref="ssh"          name="zos"      long_name="sea_surface_height_above_geoid"                /> 
    39    </file> 
    40  
    41    <file id="file2" name_suffix="_grid_U" description="ocean U grid variables" > 
    42      <field field_ref="suoce"        name="uos"     long_name="sea_surface_x_velocity"    /> 
    43    </file> 
    44     
    45    <file id="file3" name_suffix="_grid_V" description="ocean V grid variables" > 
    46      <field field_ref="svoce"        name="vos"     long_name="sea_surface_y_velocity"    /> 
    47    </file> 
    48     
    49       </file_group> 
    50  
    51       <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->     
    52  
    53       <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE.">  <!-- 5d files -->    
    54  
    55    <file id="file4" name_suffix="_grid_T" description="ocean T grid variables" > 
    5643     <field field_ref="toce"         name="thetao"   long_name="sea_water_potential_temperature"               /> 
    5744     <field field_ref="soce"         name="so"       long_name="sea_water_salinity"                            /> 
    58      <field field_ref="sst"          name="tos"      long_name="sea_surface_temperature"                       /> 
    5945     <field field_ref="sst2"         name="tossq"    long_name="square_of_sea_surface_temperature"             /> 
    60      <field field_ref="sss"          name="sos"      long_name="sea_surface_salinity"                          /> 
    61      <field field_ref="ssh"          name="zos"      long_name="sea_surface_height_above_geoid"                /> 
    6246     <field field_ref="ssh2"         name="zossq"    long_name="square_of_sea_surface_height_above_geoid"      /> 
     47     <field field_ref="mldkz5"       /> 
     48     <field field_ref="mldr10_1"     /> 
     49   </file> 
     50 
     51   <file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc --> 
    6352     <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     /> 
    6453     <field field_ref="qsr"          name="rsntds"   long_name="surface_net_downward_shortwave_flux"           /> 
    6554     <field field_ref="qt"           name="tohfls"   long_name="surface_net_downward_total_heat_flux"          /> 
    66      <field field_ref="taum"         /> 
    67      <field field_ref="mldkz5"       /> 
    68      <field field_ref="mldr10_1"     /> 
    69    </file> 
    70     
    71    <file id="file5" name_suffix="_grid_U" description="ocean U grid variables" > 
     55     <field field_ref="saltflx"      name="sosflxdo"  /> 
     56     <field field_ref="taum"         name="taum" /> 
     57     <field field_ref="wspd"         name="sowindsp"  /> 
     58          <field field_ref="precip"       name="soprecip" /> 
     59     <!-- ice and snow --> 
     60     <field field_ref="snowpre"      /> 
     61          <field field_ref="qsr_io"           name="iicesflx" /> 
     62          <field field_ref="qns_io"           name="iicenflx" /> 
     63          <field field_ref="utau_ice"         name="iicestru" /> 
     64          <field field_ref="vtau_ice"         name="iicestrv" /> 
     65 
     66   </file> 
     67 
     68   <file id="file3" name_suffix="_grid_U" description="ocean U grid variables" > 
     69     <field field_ref="suoce"        name="uos"     long_name="sea_surface_x_velocity"    /> 
    7270     <field field_ref="uoce"         name="uo"      long_name="sea_water_x_velocity"      /> 
    73      <field field_ref="suoce"        name="uos"     long_name="sea_surface_x_velocity"    /> 
    7471     <field field_ref="utau"         name="tauuo"   long_name="surface_downward_x_stress" /> 
    75    </file> 
    76     
    77    <file id="file6" name_suffix="_grid_V" description="ocean V grid variables" > 
     72          <!-- available with key_diaar5 --> 
     73     <field field_ref="u_masstr"     name="vozomatr"  /> 
     74     <field field_ref="u_heattr"     name="sozohetr"  /> 
     75   </file> 
     76    
     77   <file id="file4" name_suffix="_grid_V" description="ocean V grid variables" > 
     78     <field field_ref="svoce"        name="vos"     long_name="sea_surface_y_velocity"    /> 
    7879     <field field_ref="voce"         name="vo"      long_name="sea_water_y_velocity"      /> 
    79      <field field_ref="svoce"        name="vos"     long_name="sea_surface_y_velocity"    /> 
    8080     <field field_ref="vtau"         name="tauvo"   long_name="surface_downward_y_stress" /> 
    81    </file> 
    82     
    83    <file id="file7" name_suffix="_grid_W" description="ocean W grid variables" > 
     81          <!-- available with key_diaar5 --> 
     82     <field field_ref="v_masstr"     name="vomematr"  /> 
     83     <field field_ref="v_heattr"     name="somehetr"  /> 
     84   </file> 
     85    
     86   <file id="file5" name_suffix="_grid_W" description="ocean W grid variables" > 
    8487     <field field_ref="woce"         name="wo"      long_name="ocean vertical velocity"         /> 
    8588     <field field_ref="avt"          name="difvho"  long_name="ocean_vertical_heat_diffusivity" /> 
    86    </file> 
    87     
    88    <file id="file8" name_suffix="_icemod" description="ice variables" > 
    89      <field field_ref="ice_pres"                     /> 
     89     <field field_ref="w_masstr"     name="vovematr"  /> 
     90   </file> 
     91 
     92   <file id="file6" name_suffix="_icemod" description="ice variables" > 
    9093     <field field_ref="snowthic_cea" name="snd"     long_name="surface_snow_thickness"   /> 
    9194     <field field_ref="icethic_cea"  name="sit"     long_name="sea_ice_thickness"        /> 
    92      <field field_ref="iceprod_cea"  name="sip"     long_name="sea_ice_thickness"        /> 
    93      <field field_ref="ist_ipa"      /> 
    94      <field field_ref="ioceflxb"     /> 
    95      <field field_ref="uice_ipa"     /> 
    96      <field field_ref="vice_ipa"     /> 
    97      <field field_ref="utau_ice"     /> 
    98      <field field_ref="vtau_ice"     /> 
    99      <field field_ref="qsr_io_cea"   /> 
    100      <field field_ref="qns_io_cea"   /> 
    101      <field field_ref="snowpre"      /> 
    102    </file> 
    103     
     95          <field field_ref="icevolu"      name="iicevolu" /> 
     96          <field field_ref="snowvol"      name="isnowvol" /> 
     97          <field field_ref="iceconc"      name="iiceconc"  /> 
     98 
     99          <field field_ref="icebopr"          name="iicebopr" /> 
     100          <field field_ref="icedypr"          name="iicedypr" /> 
     101          <field field_ref="ioceflxb"         name="ioceflxb" /> 
     102          <field field_ref="uice_ipa"         name="iicevelu" /> 
     103          <field field_ref="vice_ipa"         name="iicevelv" /> 
     104          <field field_ref="isst"             name="isstempe" /> 
     105          <field field_ref="isss"             name="isssalin" /> 
     106          <field field_ref="micesalt"         name="iicesali" /> 
     107          <field field_ref="miceage"          name="iiceages" /> 
     108          <field field_ref="icelapr"          name="iicelapr" /> 
     109          <field field_ref="icesipr"          name="iicesipr" /> 
     110          <field field_ref="micet"            name="iicetemp" /> 
     111          <field field_ref="icehc"            name="iiceheco" /> 
     112          <field field_ref="isnowhc"          name="isnoheco" /> 
     113          <field field_ref="icest"            name="iicesurt" /> 
     114          <field field_ref="sfxbri"            name="iicefsbr" /> 
     115          <field field_ref="sfxthd"            name="iicefseq" /> 
     116          <field field_ref="ibrinv"           name="ibrinvol" /> 
     117          <field field_ref="icecolf"          name="iicecolf" /> 
     118          <field field_ref="icestr"           name="iicestre" /> 
     119          <field field_ref="icevel"           name="iicevelo" /> 
     120          <field field_ref="isume"            name="iicesume" /> 
     121          <field field_ref="ibome"            name="iicebome" /> 
     122          <field field_ref="idive"            name="iicedive" /> 
     123          <field field_ref="ishear"           name="iiceshea" /> 
     124          <field field_ref="icerepr"          name="iicerepr" /> 
     125          <field field_ref="sfxmec"            name="iicefsrp" /> 
     126          <field field_ref="sfxres"            name="iicefsre" /> 
     127          <field field_ref="icetrp"           name="iicevtrp" /> 
     128   </file> 
     129 
     130 
     131   <file id="file7" name_suffix="_scalar" description="scalar variables" > 
     132          <field field_ref="voltot"       name="scvoltot"  /> 
     133          <field field_ref="sshtot"       name="scsshtot"  /> 
     134          <field field_ref="sshsteric"    name="scsshste"  /> 
     135          <field field_ref="sshthster"    name="scsshtst"  /> 
     136          <field field_ref="masstot"      name="scmastot"  /> 
     137          <field field_ref="temptot"      name="sctemtot"  /> 
     138          <field field_ref="saltot"       name="scsaltot"  /> 
     139    
     140          <!-- available with ln_diahsb --> 
     141     <field field_ref="bgtemper"     name="bgtemper"    /> 
     142     <field field_ref="bgsaline"     name="bgsaline"    /> 
     143     <field field_ref="bgheatco"     name="bgheatco"    /> 
     144     <field field_ref="bgsaltco"     name="bgsaltco"    /> 
     145     <field field_ref="bgvolssh"     name="bgvolssh"    /> 
     146     <field field_ref="bgvoltot"     name="bgvoltot"    /> 
     147     <field field_ref="bgsshtot"     name="bgsshtot"    /> 
     148     <field field_ref="bgfrcvol"     name="bgfrcvol"    /> 
     149     <field field_ref="bgfrctem"     name="bgfrctem"    /> 
     150     <field field_ref="bgfrcsal"     name="bgfrcsal"    /> 
     151     <!-- available with ln_limdiahsb --> 
     152     <field field_ref="ibgvoltot"    name="ibgvoltot"   /> 
     153     <field field_ref="sbgvoltot"    name="sbgvoltot"   /> 
     154     <field field_ref="ibgarea"      name="ibgarea"     /> 
     155     <field field_ref="ibgsaline"    name="ibgsaline"   /> 
     156     <field field_ref="ibgtemper"    name="ibgtemper"   /> 
     157     <field field_ref="ibgheatco"    name="ibgheatco"   /> 
     158     <field field_ref="ibgsaltco"    name="ibgsaltco"   /> 
     159     <field field_ref="sbgheatco"    name="sbgheatco"   /> 
     160     <field field_ref="ibgfrcsfx"    name="ibgfrcsfx"  /> 
     161     <field field_ref="ibgfrcemp"    name="ibgfrcemp"   /> 
     162     <field field_ref="ibgsfx"       name="ibgsfx"     /> 
     163     <field field_ref="ibgemp"       name="ibgemp"      /> 
     164     <field field_ref="ibgsfxbri"    name="ibgsfxbri"    /> 
     165     <field field_ref="ibgsfxthd"    name="ibgsfxthd"    /> 
     166     <field field_ref="ibgsfxres"    name="ibgsfxres" /> 
     167     <field field_ref="ibgsfxmec"    name="ibgsfxmec" /> 
     168     <field field_ref="ibggrpme"     name="ibggrpme"    /> 
     169 
     170   </file> 
     171 
     172   <!-- 
     173   <file id="file8" name_suffix="_Tides" description="tidal harmonics" > 
     174     <field field_ref="M2x"          name="M2x"      long_name="M2 Elevation harmonic real part"                       /> 
     175     <field field_ref="M2y"          name="M2y"      long_name="M2 Elevation harmonic imaginary part"                  /> 
     176     <field field_ref="M2x_u"        name="M2x_u"    long_name="M2 current barotrope along i-axis harmonic real part "       /> 
     177     <field field_ref="M2y_u"        name="M2y_u"    long_name="M2 current barotrope along i-axis harmonic imaginary part "  /> 
     178     <field field_ref="M2x_v"        name="M2x_v"    long_name="M2 current barotrope along j-axis harmonic real part "       /> 
     179     <field field_ref="M2y_v"        name="M2y_v"    long_name="M2 current barotrope along j-axis harmonic imaginary part "  /> 
     180   </file> 
     181   --> 
     182 
    104183      </file_group> 
    105184 
     185      <file_group id="3d" output_freq="3d"  output_level="10" enabled=".TRUE."/> <!-- 3d files -->     
     186      <file_group id="5d" output_freq="5d"  output_level="10" enabled=".TRUE."/>  <!-- 5d files -->    
     187 
    106188      <file_group id="1m" output_freq="1mo" output_level="10" enabled=".TRUE."/> <!-- real monthly files --> 
     189 
     190 
    107191      <file_group id="2m" output_freq="2mo" output_level="10" enabled=".TRUE."/> <!-- real 2m files --> 
    108192      <file_group id="3m" output_freq="3mo" output_level="10" enabled=".TRUE."/> <!-- real 3m files --> 
     
    154238        We must have buffer_size > jpi*jpj*jpk*8 (with jpi and jpj the subdomain size) 
    155239--> 
    156      <variable id="buffer_size"               type="integer">25000000</variable> 
     240     <variable id="buffer_size"               type="integer">5000000</variable> 
    157241     <variable id="buffer_server_factor_size" type="integer">2</variable> 
    158242     <variable id="info_level"                type="integer">0</variable> 
    159      <variable id="using_server"              type="boolean">false</variable> 
     243     <variable id="using_server"              type="boolean">true</variable> 
    160244     <variable id="using_oasis"               type="boolean">false</variable> 
    161245     <variable id="oasis_codes_id"            type="string" >oceanx</variable> 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/ORCA2_LIM3/cpp_ORCA2_LIM3.fcm

    r4099 r4161  
    1  bld::tool::fppkeys key_trabbl key_orca_r2 key_lim3 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi 
     1 bld::tool::fppkeys key_trabbl key_lim3 key_dynspg_flt key_diaeiv key_ldfslp key_traldf_c2d key_traldf_eiv key_dynldf_c3d key_zdftke key_zdfddm key_zdftmx key_iomput key_mpp_mpi 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/1_namelist_ref

    r4153 r4161  
    281281   ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data 
    282282   rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
     283   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
     284   rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
     285                           !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 
    283286/ 
    284287!----------------------------------------------------------------------- 
     
    350353   rn_si0      =   0.35    !  RGB & 2 bands: shortess depth of extinction 
    351354   rn_si1      =   23.0    !  2 bands: longest depth of extinction 
     355   ln_qsr_ice  = .true.    !  light penetration for ice-model LIM3 
    352356/ 
    353357!----------------------------------------------------------------------- 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/field_def.xml

    r4153 r4161  
    1515      <!-- T grid --> 
    1616       
    17      <field_group id="grid_T" grid_ref="grid_T_2D"> 
     17     <field_group id="grid_T" grid_ref="grid_T_2D" > 
    1818         <field id="toce"         long_name="temperature"                               unit="degC" grid_ref="grid_T_3D"/> 
    1919         <field id="soce"         long_name="salinity"                                  unit="psu"  grid_ref="grid_T_3D"/> 
     
    5050     </field_group> 
    5151 
     52     <field_group id="Tides_T" grid_ref="grid_T_2D" operation="once" > 
     53    <!-- tidal composante --> 
     54    <field id="M2x"          long_name="M2 Elevation harmonic real part "                             unit="m"        /> 
     55    <field id="M2y"          long_name="M2 Elevation harmonic imaginary part"                         unit="m"        /> 
     56    <field id="S2x"          long_name="S2 Elevation harmonic real part "                             unit="m"        /> 
     57    <field id="S2y"          long_name="S2 Elevation harmonic imaginary part"                         unit="m"        /> 
     58    <field id="N2x"          long_name="N2 Elevation harmonic real part "                             unit="m"        /> 
     59    <field id="N2y"          long_name="N2 Elevation harmonic imaginary part"                         unit="m"        /> 
     60    <field id="K1x"          long_name="K1 Elevation harmonic real part "                             unit="m"        /> 
     61    <field id="K1y"          long_name="K1 Elevation harmonic imaginary part"                         unit="m"        /> 
     62    <field id="O1x"          long_name="O1 Elevation harmonic real part "                             unit="m"        /> 
     63    <field id="O1y"          long_name="O1 Elevation harmonic imaginary part"                         unit="m"        /> 
     64    <field id="Q1x"          long_name="Q1 Elevation harmonic real part "                             unit="m"        /> 
     65    <field id="Q1y"          long_name="Q1 Elevation harmonic imaginary part"                         unit="m"        /> 
     66    <field id="M4x"          long_name="M4 Elevation harmonic real part "                             unit="m"        /> 
     67    <field id="M4y"          long_name="M4 Elevation harmonic imaginary part"                         unit="m"        /> 
     68    <field id="K2x"          long_name="K2 Elevation harmonic real part "                             unit="m"        /> 
     69    <field id="K2y"          long_name="K2 Elevation harmonic imaginary part"                         unit="m"        /> 
     70    <field id="P1x"          long_name="P1 Elevation harmonic real part "                             unit="m"        /> 
     71    <field id="P1y"          long_name="P1 Elevation harmonic imaginary part"                         unit="m"        /> 
     72    <field id="Mfx"          long_name="Mf Elevation harmonic real part "                             unit="m"        /> 
     73    <field id="Mfy"          long_name="Mf Elevation harmonic imaginary part"                         unit="m"        /> 
     74    <field id="Mmx"          long_name="Mm Elevation harmonic real part "                             unit="m"        /> 
     75    <field id="Mmy"          long_name="Mm Elevation harmonic imaginary part"                         unit="m"        /> 
     76     </field_group> 
     77     
     78     <field_group id="Tides_U" grid_ref="grid_U_2D" operation="once" > 
     79    <field id="M2x_u"          long_name="M2 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     80    <field id="M2y_u"          long_name="M2 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     81    <field id="S2x_u"          long_name="S2 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     82    <field id="S2y_u"          long_name="S2 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     83    <field id="N2x_u"          long_name="N2 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     84    <field id="N2y_u"          long_name="N2 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     85    <field id="K1x_u"          long_name="K1 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     86    <field id="K1y_u"          long_name="K1 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     87    <field id="O1x_u"          long_name="O1 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     88    <field id="O1y_u"          long_name="O1 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     89    <field id="Q1x_u"          long_name="Q1 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     90    <field id="Q1y_u"          long_name="Q1 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     91    <field id="M4x_u"          long_name="M4 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     92    <field id="M4y_u"          long_name="M4 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     93    <field id="K2x_u"          long_name="K2 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     94    <field id="K2y_u"          long_name="K2 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     95    <field id="P1x_u"          long_name="P1 current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     96    <field id="P1y_u"          long_name="P1 current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     97    <field id="Mfx_u"          long_name="Mf current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     98    <field id="Mfy_u"          long_name="Mf current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     99    <field id="Mmx_u"          long_name="Mm current barotrope along i-axis harmonic real part "      unit="m/s"        /> 
     100    <field id="Mmy_u"          long_name="Mm current barotrope along i-axis harmonic imaginary part " unit="m/s"        /> 
     101     </field_group> 
     102     
     103     <field_group id="Tides_V" grid_ref="grid_V_2D" operation="once" > 
     104    <field id="M2x_v"          long_name="M2 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     105    <field id="M2y_v"          long_name="M2 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     106    <field id="S2x_v"          long_name="S2 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     107    <field id="S2y_v"          long_name="S2 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     108    <field id="N2x_v"          long_name="N2 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     109    <field id="N2y_v"          long_name="N2 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     110    <field id="K1x_v"          long_name="K1 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     111    <field id="K1y_v"          long_name="K1 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     112    <field id="O1x_v"          long_name="O1 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     113    <field id="O1y_v"          long_name="O1 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     114    <field id="Q1x_v"          long_name="Q1 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     115    <field id="Q1y_v"          long_name="Q1 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     116    <field id="M4x_v"          long_name="M4 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     117    <field id="M4y_v"          long_name="M4 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     118    <field id="K2x_v"          long_name="K2 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     119    <field id="K2y_v"          long_name="K2 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     120    <field id="P1x_v"          long_name="P1 current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     121    <field id="P1y_v"          long_name="P1 current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     122    <field id="Mfx_v"          long_name="Mf current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     123    <field id="Mfy_v"          long_name="Mf current barotrope along j-axis harmonic imaginary part " unit="m/s"        /> 
     124    <field id="Mmx_v"          long_name="Mm current barotrope along j-axis harmonic real part "      unit="m/s"        /> 
     125    <field id="Mmy_v"          long_name="Mm current barotrope along j-axis harmonic imaginary part " unit="m/s"        />   
     126     </field_group> 
     127 
     128 
    52129      <!-- SBC --> 
    53130       
     
    59136         <field id="snowpre"      long_name="Snow precipitation"                                           unit="kg/m2/s"  /> 
    60137         <field id="runoffs"      long_name="River Runoffs"                                                unit="Kg/m2/s"  /> 
     138         <field id="precip"       long_name="Total precipitation"                                           unit="kg/m2/s"  /> 
     139 
    61140 
    62141         <field id="qt"           long_name="Net Downward Heat Flux"                                       unit="W/m2"     /> 
     
    79158         <field id="qhc_oce"      long_name="Downward Heat Content of E-P over open ocean"                 unit="W/m2"     /> 
    80159         <field id="taum_oce"     long_name="wind stress module over open ocean"                           unit="N/m2"     /> 
    81  
    82          <field id="ice_cover"    long_name="Ice fraction"                                                 unit="1"        /> 
    83  
    84          <field id="ioceflxb"     long_name="Oceanic flux at the ice base"                                 unit="W/m2"     /> 
    85          <field id="qsr_ai_cea"   long_name="Air-Ice downward solar heat flux (cell average)"              unit="W/m2"     /> 
    86          <field id="qns_ai_cea"   long_name="Air-Ice downward non-solar heat flux (cell average)"          unit="W/m2"     /> 
    87          <field id="qla_ai_cea"   long_name="Air-Ice downward Latent heat flux (cell average)"             unit="W/m2"     /> 
    88           
    89          <field id="qsr_io_cea"   long_name="Ice-Oce downward solar heat flux (cell average)"              unit="W/m2"     /> 
    90          <field id="qns_io_cea"   long_name="Ice-Oce downward non-solar heat flux (cell average)"          unit="W/m2"     /> 
    91           
    92          <field id="snowthic_cea" long_name="Snow thickness (cell average)"                                unit="m"        /> 
    93          <field id="icethic_cea"  long_name="Ice thickness (cell average)"                                 unit="m"        /> 
    94          <field id="iceprod_cea"  long_name="Ice production (cell average)"                                unit="m/s"      /> 
    95           
    96          <field id="ice_pres"     long_name="Ice presence"                                                 unit="-"        /> 
    97          <field id="ist_cea"      long_name="Ice surface temperature (cell average)"                       unit="degC"     /> 
    98          <field id="ist_ipa"      long_name="Ice surface temperature (ice presence average)"               unit="degC"     />       
    99          <field id="uice_ipa"     long_name="Ice velocity along i-axis at I-point (ice presence average)"  unit="m/s"      />       
    100          <field id="vice_ipa"     long_name="Ice velocity along j-axis at I-point (ice presence average)"  unit="m/s"      />       
    101           
    102          <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             unit="N/m2"     /> 
    103          <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             unit="N/m2"     /> 
    104           
    105          <field id="u_imasstr"    long_name="Sea-ice mass transport along i-axis"                          unit="kg/s"     /> 
    106          <field id="v_imasstr"    long_name="Sea-ice mass transport along j-axis"                          unit="kg/s"     /> 
    107          <field id="emp_x_sst"    long_name="Concentration/Dilution term on SST"                           unit="kgC/m2/s" /> 
    108          <field id="emp_x_sss"    long_name="Concentration/Dilution term on SSS"                         unit="kgPSU/m2/s" /> 
    109160 
    110161         <!-- available key_coupled --> 
     
    132183         <field id="sntoice_cea"   long_name="Snow-Ice Formation Rate (cell average)"                      unit="kg/m2/s"  /> 
    133184         <field id="ticemel_cea"      long_name="Rate of Melt at Upper Surface of Sea Ice (cell average)"     unit="kg/m2/s"  /> 
     185 
     186    <!-- ice fields --> 
     187 
     188         <field id="ice_cover"    long_name="Ice fraction"                                                 unit="1"        /> 
     189 
     190         <field id="ioceflxb"     long_name="Oceanic flux at the ice base"                                 unit="W/m2"     /> 
     191         <field id="qsr_ai_cea"   long_name="Air-Ice downward solar heat flux (cell average)"              unit="W/m2"     /> 
     192         <field id="qns_ai_cea"   long_name="Air-Ice downward non-solar heat flux (cell average)"          unit="W/m2"     /> 
     193         <field id="qla_ai_cea"   long_name="Air-Ice downward Latent heat flux (cell average)"             unit="W/m2"     /> 
     194          
     195         <field id="qsr_io_cea"   long_name="Ice-Oce downward solar heat flux (cell average)"              unit="W/m2"     /> 
     196         <field id="qns_io_cea"   long_name="Ice-Oce downward non-solar heat flux (cell average)"          unit="W/m2"     /> 
     197          
     198         <field id="snowthic_cea" long_name="Snow thickness (cell average)"                                unit="m"        /> 
     199         <field id="icethic_cea"  long_name="Ice thickness (cell average)"                                 unit="m"        /> 
     200         <field id="iceprod_cea"  long_name="Ice production (cell average)"                                unit="m/s"      /> 
     201          
     202         <field id="ice_pres"     long_name="Ice presence"                                                 unit="-"        /> 
     203         <field id="ist_cea"      long_name="Ice surface temperature (cell average)"                       unit="degC"     /> 
     204         <field id="ist_ipa"      long_name="Ice surface temperature (ice presence average)"               unit="degC"     />       
     205         <field id="uice_ipa"     long_name="Ice velocity along i-axis at I-point (ice presence average)"  unit="m/s"      />       
     206         <field id="vice_ipa"     long_name="Ice velocity along j-axis at I-point (ice presence average)"  unit="m/s"      />       
     207          
     208         <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             unit="N/m2"     /> 
     209         <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             unit="N/m2"     /> 
     210          
     211         <field id="u_imasstr"    long_name="Sea-ice mass transport along i-axis"                          unit="kg/s"     /> 
     212         <field id="v_imasstr"    long_name="Sea-ice mass transport along j-axis"                          unit="kg/s"     /> 
     213         <field id="emp_x_sst"    long_name="Concentration/Dilution term on SST"                           unit="kgC/m2/s" /> 
     214         <field id="emp_x_sss"    long_name="Concentration/Dilution term on SSS"                         unit="kgPSU/m2/s" />         
     215          
     216          
     217         <field id="iceconc"      long_name="ice concentration"                                            unit="%"        /> 
     218         <field id="icebopr"      long_name="daily bottom thermo ice prod."                                unit="km3/day"   /> 
     219         <field id="icedypr"      long_name="daily  dynamic ice prod."                                     unit="km3/day"   /> 
     220    <field id="ioceflxb"     long_name="Oceanic flux at the ice base"                                 unit="W/m2"     /> 
     221         <field id="uice_ipa"     long_name="Ice velocity along i-axis at I-point (ice presence average)"  unit="m/s"      /> 
     222         <field id="vice_ipa"     long_name="Ice velocity along j-axis at I-point (ice presence average)"  unit="m/s"      /> 
     223          <field id="isst"         long_name="sea surface temperature"                                      unit="degC"     /> 
     224         <field id="isss"         long_name="sea surface salinity"                                         unit="psu"      />  
     225         <field id="qt_oce"       long_name="total flux at ocean surface"                                  unit="W/m2"     /> 
     226         <field id="qsr_oce"      long_name="solar heat flux at ocean surface"                             unit="W/m2"     /> 
     227         <field id="qns_oce"      long_name="non-solar heat flux at ocean surface"                         unit="W/m2"     /> 
     228         <field id="hfbri"        long_name="heat flux due to brine release"                               unit="W/m2"     /> 
     229         <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             unit="N/m2"     /> 
     230         <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             unit="N/m2"     /> 
     231    <field id="qsr_io"       long_name="Ice-Oce downward solar heat flux"                             unit="W/m2"     /> 
     232    <field id="qns_io"       long_name="Ice-Oce downward non-solar heat flux"                         unit="W/m2"     /> 
     233         <field id="micesalt"     long_name="Mean ice salinity"                                            unit="psu"      /> 
     234         <field id="miceage"      long_name="Mean ice age"                                                 unit="years"    /> 
     235         <field id="icelapr"      long_name="daily lateral thermo ice prod."                               unit="km3/day"   /> 
     236         <field id="icesipr"      long_name="daily snowice ice prod."                                      unit="km3/day"   /> 
     237         <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     /> 
     238         <field id="icehc"        long_name="ice total heat content"                                       unit="10^9 J"   />  
     239         <field id="isnowhc"      long_name="snow total heat content"                                      unit="10^9J"    /> 
     240         <field id="icest"        long_name="ice surface temperature"                                      unit="degC"     /> 
     241         <field id="sfxbri"       long_name="brine salt flux"                                              unit="psu*kg/m2/day" /> 
     242         <field id="sfxthd"       long_name="equivalent FW salt flux"                                      unit="psu*kg/m2/day" /> 
     243         <field id="ibrinv"       long_name="brine volume"                                                 unit="%"        /> 
     244         <field id="icecolf"      long_name="frazil ice collection thickness"                              unit="m"        /> 
     245         <field id="icestr"       long_name="ice strength"                                                 unit="N/m"      /> 
     246         <field id="icevel"       long_name="ice velocity"                                                 unit="m/s"      /> 
     247         <field id="isume"        long_name="surface melt"                                                 unit="km3/day"   /> 
     248         <field id="ibome"        long_name="bottom melt"                                                  unit="km3/day"   /> 
     249         <field id="idive"        long_name="divergence"                                                   unit="10-8s-1"  /> 
     250         <field id="ishear"       long_name="shear"                                                        unit="10-8s-1"  /> 
     251         <field id="icerepr"      long_name="daily resultant ice prod./melting from limupdate"             unit="km3/day"   /> 
     252         <field id="icevolu"      long_name="ice volume"                                                   unit="km3"      /> 
     253         <field id="snowvol"      long_name="snow volume"                                                  unit="km3"      /> 
     254         <field id="sfxmec"       long_name="salt flux from ridging rafting"                               unit="psu*kg/m2/day" /> 
     255         <field id="sfxres"       long_name="salt flux from lipupdate (resultant)"                         unit="psu*kg/m2/day" /> 
     256         <field id="icetrp"       long_name="ice volume transport"                                         unit="km3/day"   /> 
     257 
    134258 
    135259      </field_group> 
     
    214338         <field id="saltot"     long_name="global mean salinity"                       unit="psu"  /> 
    215339         <field id="fram_trans" long_name="Sea Ice Mass Transport Through Fram Strait" unit="kg/s" /> 
     340       <!-- available with ln_diahsb --> 
     341    <field id="bgtemper"     long_name="global mean temperature"                  unit="degC"   /> 
     342    <field id="bgsaline"     long_name="global mean salinity"                     unit="psu"    /> 
     343    <field id="bgheatco"     long_name="global mean heat content"                 unit="10^9J"  /> 
     344    <field id="bgsaltco"     long_name="global mean salt content"                 unit="psu*m3" /> 
     345    <field id="bgvolssh"     long_name="global mean ssh volume"                   unit="km3"     /> 
     346    <field id="bgvoltot"     long_name="global mean volume"                       unit="km3"     /> 
     347    <field id="bgsshtot"     long_name="global mean ssh"                          unit="m"      /> 
     348    <field id="bgfrcvol"     long_name="global mean volume from forcing"          unit="km3"     /> 
     349    <field id="bgfrctem"     long_name="global mean heat content from forcing"    unit="10^9J"  /> 
     350    <field id="bgfrcsal"     long_name="global mean salt content from forcing"    unit="psu*km3" /> 
     351      </field_group> 
     352 
     353      <field_group id="SBC_scalar"  domain_ref="1point" > 
     354         <!-- available with ln_limdiahsb --> 
     355    <field id="ibgvoltot"    long_name="global mean ice volume"                   unit="km3"   /> 
     356    <field id="sbgvoltot"    long_name="global mean snow volume"                  unit="km3"   /> 
     357    <field id="ibgarea"      long_name="global mean ice area"                     unit="km2"   /> 
     358    <field id="ibgsaline"    long_name="global mean ice salinity"                 unit="psu"   /> 
     359    <field id="ibgtemper"    long_name="global mean ice temperature"              unit="degC"   /> 
     360    <field id="ibgheatco"    long_name="global mean ice heat content"             unit="10^9J"   /> 
     361    <field id="ibgsaltco"    long_name="global mean ice salt content"             unit="psu*km3"   /> 
     362    <field id="sbgheatco"    long_name="global mean snow heat content"            unit="10^9J"   /> 
     363    <field id="ibgfrcsfx"    long_name="global mean salt content from sfx"        unit="psu*km3"   /> 
     364    <field id="ibgfrcemp"    long_name="global mean volume from emp"              unit="km3"      /> 
     365    <field id="ibgsfx"       long_name="global mean emps"                         unit="psu*kg/m2/day"   /> 
     366    <field id="ibgemp"       long_name="global mean emp"                          unit="kg/m2/day"   /> 
     367    <field id="ibgsfxbri"    long_name="global mean ice sfx_bri"                  unit="psu*kg/m2/day"   /> 
     368    <field id="ibgsfxthd"    long_name="global mean ice sfx_thd"                  unit="psu*kg/m2/day"   /> 
     369    <field id="ibgsfxres"    long_name="global mean ice sfx_res"                  unit="psu*kg/m2/day"   /> 
     370    <field id="ibgsfxmec"    long_name="global mean ice fsalt_rpo"                unit="psu*kg/m2/day"   /> 
     371    <field id="ibggrpme"     long_name="global mean ice growth+melt volume"       unit="km3"      /> 
    216372      </field_group> 
    217373   
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim2_ref

    r4148 r4161  
    5454   telast      =   9600    !  timescale for EVP elastic waves 
    5555   alphaevp    =   1.0     !  coefficient for the solution of EVP int. stresses 
     56   hminrhg     =   0.05     !  ice thickness (m) below which ice velocity equal ocean velocity 
    5657/ 
    5758!----------------------------------------------------------------------- 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref

    r4147 r4161  
    1414&namicerun     !   Share parameters for dynamics/advection/thermo 
    1515!----------------------------------------------------------------------- 
    16    cn_icerst_in  = "restart_ice_in"   !  suffix of ice restart name (input) 
     16   cn_icerst_in  = "restart_ice"   !  suffix of ice restart name (input) 
    1717   cn_icerst_out = "restart_ice"      !  suffix of ice restart name (output) 
    1818   ln_limdyn   = .true.    !  ice dynamics (T) or thermodynamics only (F) 
    19    acrit       =  1.0e-02 , 1.0e-02  !  minimum fraction for leads in the Northern (Southern) Hemisphere 
    20    hsndif      =  0.0      !  computation of temperature in snow (=0.0) or not 
    21    hicdif      =  0.0      !  computation of temperature in ice  (=0.0) or not (=9999.0) 
     19   amax        = 0.999      !  maximum ice concentration 
    2220   cai         =  1.40e-3  !  atmospheric drag over sea ice 
    2321   cao         =  1.00e-3  !  atmospheric drag over ocean 
    2422   ln_nicep    = .false.   !  Ice points output for debug (yes or no) 
     23   ln_limdiahsb  = .false.    !  check the heat and salt budgets (T) or not (F) 
     24   ln_limdiaout  = .false.    !  output the heat and salt budgets (T) or not (F) 
    2525/ 
    2626!----------------------------------------------------------------------- 
     
    3030   hninn       =  0.3      !  initial snow thickness in the north 
    3131   hnins       =  0.1      !        "            "          south 
    32    hginn_u     =  3.50     !  initial undeformed ice thickness in the north 
    33    hgins_u     =  1.0      !        "            "              "     south 
    34    aginn_u     =  0.95     !  initial undeformed ice concentration in the north 
    35    agins_u     =  0.9      !        "            "              "         south 
    36    hginn_d     =  0.0      !  initial  deformed  ice thickness     in the north 
    37    hgins_d     =  0.0      !   
    38    aginn_d     =  0.00     !  initial  deformed  ice concentration in the north 
    39    agins_d     =  0.00     !        "            "              "         south 
     32   hginn       =  3.50     !  initial undeformed ice thickness in the north 
     33   hgins       =  1.0      !        "            "              "     south 
     34   aginn       =  0.95     !  initial undeformed ice concentration in the north 
     35   agins       =  0.9      !        "            "              "         south 
    4036   sinn        =  6.301    !  initial salinity in the north 
    4137   sins        =  6.301    !        "            "    south 
     
    6258   telast      =9600.0     !  timescale for elastic waves, SB, 720.0 
    6359   alphaevp    =   1.0     !  coefficient for the solution of internal ice stresses 
     60   hminrhg     =   0.05     !  ice thickness (m) below which ice velocity equal ocean velocity 
    6461/ 
    6562!----------------------------------------------------------------------- 
     
    8077   hicmin      = 0.2       !  ice thickness corr. to max. energy stored in brine pocket 
    8178   hiclim      = 0.10      !  minimum ice thickness 
    82    amax        = 0.999     !  maximum lead fraction 
    8379   sbeta       = 1.        !  numerical caracteritic of the scheme for diffusion in ice 
    8480                           !        Cranck-Nicholson (=0.5), implicit (=1), explicit (=0) 
     
    145141   ninfo       = 1         !  frequency of ouputs on file ice_evolu in case of averaging 
    146142/ 
     143!!----------------------------------------------------------------------- 
     144!&namicehsb       !  Heat and salt budgets  
     145!!----------------------------------------------------------------------- 
     146! 
     147!/ 
    147148!----------------------------------------------------------------------- 
    148149&namiceout     !   parameters for outputs 
    149150!----------------------------------------------------------------------- 
    150    noumef      =   37      !  number of fields 
     151   noumef      =   43      !  number of fields 
    151152   add_diag_swi=    1      !  1 -> diagnose distribution in thickness space 
    152153                           !  0 -> only simple diagnostics 
     
    157158   field_2  = 'Ice thickness                      ', 'iicethic', 'm       ',    1   ,  1.0     ,    0.0 
    158159   field_3  = 'Snow thickness                     ', 'isnowthi', 'm       ',    1   ,  1.0     ,    0.0 
    159    field_4  = 'Daily bottom thermo ice production ', 'iicebopr', 'cm/day  ',    1   , 100.     ,    0.0 
    160    field_5  = 'Daily dynamic ice production       ', 'iicedypr', 'cm/day  ',    1   , 100.     ,    0.0 
     160   field_4  = 'Daily bottom thermo ice production ', 'iicebopr', 'km3/day ',    1   , 1.0e-9   ,    0.0 
     161   field_5  = 'Daily dynamic ice production       ', 'iicedypr', 'km3/day ',    1   , 1.0e-9   ,    0.0 
    161162   field_6  = 'Oceanic flux at the ice base       ', 'ioceflxb', 'w/m2    ',    1   ,  1.0     ,    0.0 
    162163   field_7  = 'Ice velocity u                     ', 'iicevelu', 'm/s     ',    1   ,  1.0     ,    0.0 
     
    172173   field_17 = 'Solar flux at ice/ocean surface    ', 'iicesflx', 'w/m2    ',    1   ,  1.0     ,    0.0 
    173174   field_18 = 'Non-solar flux at ice/ocean surface', 'iicenflx', 'w/m2    ',    1   ,  1.0     ,    0.0 
    174    field_19 = 'Snow precipitation                 ', 'isnowpre', 'kg/day ',    1   ,  1.0     ,    0.0 
     175   field_19 = 'Snow precipitation                 ', 'isnowpre', 'kg/m2/d ',    1   ,  1.0     ,    0.0 
    175176   field_20 = 'Mean ice salinity                  ', 'iicesali', 'psu     ',    1   ,  1.0     ,    0.0 
    176177   field_21 = 'Mean ice age                       ', 'iiceages', 'years   ',    1   ,  0.002739,    0.0 
    177    field_22 = 'Daily lateral thermo ice prod.     ', 'iicelapr', 'cm/day  ',    1   ,100.      ,    0.0 
    178    field_23 = 'Daily snowice ice production       ', 'iicesipr', 'cm/day  ',    1   ,100.      ,    0.0 
     178   field_22 = 'Daily lateral thermo ice prod.     ', 'iicelapr', 'km3/day ',    1   ,1.0e-9    ,    0.0 
     179   field_23 = 'Daily snowice ice production       ', 'iicesipr', 'km3/day ',    1   ,1.0e-9    ,    0.0 
    179180   field_24 = 'Mean ice temperature               ', 'iicetemp', 'C       ',    1   ,  1.0     , -273.15 
    180181   field_25 = 'Ice total heat content             ', 'iiceheco', '10^9 J  ',    1   ,  1.0     ,    0.0 
    181182   field_26 = 'Ice surface temperature            ', 'iicesurt', 'C       ',    1   ,  1.0     , -273.15 
    182183   field_27 = 'Snow temperature                   ', 'isnotem2', 'C       ',    1   ,  1.0     , -273.15 
    183    field_28 = 'Fsbri - brine salt flux            ', 'iicfsbri', 'kg/m2/s ',    1   ,  1.0     ,    0.0 
    184    field_29 = 'Fseqv - equivalent FW salt flux    ', 'iicfseqv', 'kg/m2/s ',    1   ,  1.0     ,    0.0 
     184   field_28 = 'Fsbri - brine salt flux            ', 'iicefsbr', 'kg/m2/d ',    1   ,  1.0     ,    0.0 
     185   field_29 = 'Fseqv - equivalent FW salt flux    ', 'iicefseq', 'kg/m2/d ',    1   ,  1.0     ,    0.0 
    185186   field_30 = 'Brine volume                       ', 'ibrinvol', '%       ',    1   ,  100.0   ,    0.0 
    186187   field_31 = 'Frazil ice collection thickness    ', 'iicecolf', 'm       ',    1   ,  1.0     ,    0.0 
    187188   field_32 = 'Ice strength                       ', 'iicestre', 'N/m     ',    1   ,  0.001   ,    0.0 
    188189   field_33 = 'Ice velocity                       ', 'iicevelo', 'm/s     ',    1   ,  1.0     ,    0.0 
    189    field_34 = 'Surface melt                       ', 'iicesume', 'cm/day  ',    1   ,100.      ,    0.0 
    190    field_35 = 'Bottom melt                        ', 'iicebome', 'cm/day  ',    1   ,100.      ,    0.0 
     190   field_34 = 'Surface melt                       ', 'iicesume', 'km3/day ',    1   ,1.0e-9    ,    0.0 
     191   field_35 = 'Bottom melt                        ', 'iicebome', 'km3/day ',    1   ,1.0e-9    ,    0.0 
    191192   field_36 = 'Divergence                         ', 'iicedive', '10-8s-1 ',    1   ,  1.0e8   ,    0.0 
    192193   field_37 = 'Shear                              ', 'iiceshea', '10-8s-1 ',    1   ,  1.0e8   ,    0.0 
    193 /       
     194   field_38 = 'Daily resultant ice prod/melt      ', 'iicerepr', 'km3/day ',    1   ,  1.0e-9  ,    0.0 
     195   field_39 = 'Ice volume                         ', 'iicevolu', 'km3     ',    1   ,  1.0e-9  ,    0.0 
     196   field_40 = 'Snow volume                        ', 'isnowvol', 'km3     ',    1   ,  1.0e-9  ,    0.0 
     197   field_41 = 'Fsrpo - salt flux from ridg/raft   ', 'iicefsrp', 'kg/m2/d ',    1   ,  1.0     ,    0.0 
     198   field_42 = 'Fsres - salt flux from limupdate   ', 'iicefsre', 'kg/m2/d ',    1   ,  1.0     ,    0.0 
     199   field_43 = 'Ice volume transport               ', 'iicevtrp', 'km3/day ',    1   ,1.0e-9    ,    0.0 
     200/  
     201 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/SHARED/namelist_ref

    r4152 r4161  
    283283   ln_taudif   = .false.   !  HF tau contribution: use "mean of stress module - module of the mean stress" data 
    284284   rn_pfac     = 1.        !  multiplicative factor for precipitation (total & snow) 
     285   rn_vfac     = 0.        !  multiplicative factor for ocean/ice velocity 
     286   rn_efac     = 1.        !  multiplicative factor for evaporation (0. or 1.) 
     287                           !  in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 
    285288/ 
    286289!----------------------------------------------------------------------- 
     
    352355   rn_si0      =   0.35    !  RGB & 2 bands: shortess depth of extinction 
    353356   rn_si1      =   23.0    !  2 bands: longest depth of extinction 
     357   ln_qsr_ice  = .true.    !  light penetration for ice-model LIM3 
    354358/ 
    355359!----------------------------------------------------------------------- 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/CONFIG/cfg.txt

    r4159 r4161  
    1 ORCA2_LIM3 OPA_SRC LIM_SRC_3 
     1GYRE OPA_SRC 
    22GYRE_BFM OPA_SRC TOP_SRC 
    3 GYRE OPA_SRC 
     3GYRE_PISCES OPA_SRC TOP_SRC 
    44AMM12 OPA_SRC 
     5ORCA2_LIM OPA_SRC LIM_SRC_2 NST_SRC 
     6ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 
     7ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
     8ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    59ORCA2_SAS_LIM OPA_SRC SAS_SRC LIM_SRC_2 NST_SRC 
    6 ORCA2_LIM_CFC_C14b OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    7 GYRE_PISCES OPA_SRC TOP_SRC 
    8 ORCA2_LIM_PISCES OPA_SRC LIM_SRC_2 NST_SRC TOP_SRC 
    910ORCA2_OFF_PISCES OPA_SRC OFF_SRC TOP_SRC 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_2/ice_2.F90

    r4147 r4161  
    5151   REAL(wp), PUBLIC ::   ahi0        !: sea-ice hor. eddy diffusivity coeff. (m2/s) 
    5252   REAL(wp), PUBLIC ::   alphaevp    !: coefficient for the solution of EVP int. stresses 
     53   REAL(wp), PUBLIC ::   hminrhg = 0.001_wp    !: clem : ice volume (a*h in m) below which ice velocity is set to ocean velocity 
    5354 
    5455   REAL(wp), PUBLIC ::   usecc2                !:  = 1.0 / ( ecc * ecc ) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_2/limdyn_2.F90

    r4147 r4161  
    228228         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    229229         &                c_rhg, etamn, creepl, ecc, ahi0,                  & 
    230          &                nevp, telast,alphaevp 
     230         &                nevp, telast, alphaevp, hminrhg 
    231231      !!------------------------------------------------------------------- 
    232232                     
     
    262262         WRITE(numout,*) '       timescale for elastic waves telast = ', telast 
    263263         WRITE(numout,*) '       coefficient for the solution of int. stresses alphaevp = ', alphaevp 
     264         WRITE(numout,*) '       min ice thickness for rheology calculations     hminrhg = ', hminrhg 
    264265      ENDIF 
    265266      ! 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/dom_ice.F90

    r3625 r4161  
    3131 
    3232   !!---------------------------------------------------------------------- 
    33    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     33   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    3434   !! $Id$ 
    3535   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4147 r4161  
    188188   REAL(wp), PUBLIC ::   alphaevp         !: coeficient of the internal stresses !SB 
    189189   REAL(wp), PUBLIC ::   unit_fac = 1.e+09_wp  !: conversion factor for ice / snow enthalpy 
     190   REAL(wp), PUBLIC ::   hminrhg = 0.001_wp    !: clem : ice volume (a*h, in m) below which ice velocity is set to ocean velocity 
    190191 
    191192   !                                     !!** ice-salinity namelist (namicesal) ** 
     
    405406   !!-------------------------------------------------------------------------- 
    406407   !! Check if everything down here is necessary 
     408   LOGICAL , PUBLIC                                      ::   ln_limdiahsb  !: flag for ice diag (T) or not (F) 
     409   LOGICAL , PUBLIC                                      ::   ln_limdiaout  !: flag for ice diag (T) or not (F) 
    407410   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   v_newice   !: volume of ice formed in the leads 
    408411   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   dv_dt_thd  !: thermodynamic growth rates  
     
    414417   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_bot_me   ! vertical bottom melt  
    415418   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_sur_me   ! vertical surface melt 
     419   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_res_pr   ! production (growth+melt) due to limupdate 
     420   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_vi   ! transport of ice volume 
    416421   INTEGER , PUBLIC ::   jiindx, jjindx        !: indexes of the debugging point 
    417422 
    418423   !!---------------------------------------------------------------------- 
    419    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
     424   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    420425   !! $Id$ 
    421426   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    527532         &      izero    (jpi,jpj,jpl) , diag_bot_gr(jpi,jpj) , diag_dyn_gr(jpi,jpj) ,     & 
    528533         &      fstroc   (jpi,jpj,jpl) , diag_bot_me(jpi,jpj) , diag_sur_me(jpi,jpj) ,     & 
    529          &      fhbricat (jpi,jpj,jpl) , v_newice   (jpi,jpj)                        , STAT=ierr(ii) ) 
     534         &      fhbricat (jpi,jpj,jpl) , diag_res_pr(jpi,jpj) , diag_trp_vi(jpi,jpj) , v_newice(jpi,jpj) , STAT=ierr(ii) ) 
    530535 
    531536      ice_alloc = MAXVAL( ierr(:) ) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/iceini.F90

    r4147 r4161  
    4040 
    4141   !!---------------------------------------------------------------------- 
    42    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     42   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4343   !! $Id$ 
    4444   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    128128      !! ** input   :   Namelist namicerun 
    129129      !!------------------------------------------------------------------- 
    130       NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, acrit, hsndif, hicdif, cai, cao, ln_nicep 
     130      NAMELIST/namicerun/ cn_icerst_in, cn_icerst_out, ln_limdyn, amax, cai, cao, ln_nicep, ln_limdiahsb, ln_limdiaout 
    131131      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    132132      !!------------------------------------------------------------------- 
     
    144144         ln_nicep = .FALSE. 
    145145         CALL ctl_warn( 'ice_run : specific control print for LIM3 desactivated with MPI' ) 
    146       ENDIF        
     146      ENDIF 
    147147      ! 
    148148      IF(lwp) THEN                        ! control print 
     
    151151         WRITE(numout,*) ' ~~~~~~' 
    152152         WRITE(numout,*) '   switch for ice dynamics (1) or not (0)      ln_limdyn   = ', ln_limdyn 
    153          WRITE(numout,*) '   minimum fraction for leads in the NH (SH)  acrit(1/2)   = ', acrit(:) 
    154          WRITE(numout,*) '   computation of temp. in snow (=0) or not (=9999) hsndif = ', hsndif 
    155          WRITE(numout,*) '   computation of temp. in ice  (=0) or not (=9999) hicdif = ', hicdif 
     153         WRITE(numout,*) '   maximum ice concentration                               = ', amax  
    156154         WRITE(numout,*) '   atmospheric drag over sea ice                           = ', cai 
    157155         WRITE(numout,*) '   atmospheric drag over ocean                             = ', cao 
    158156         WRITE(numout,*) '   Several ice points in the ice or not in ocean.output    = ', ln_nicep 
     157         WRITE(numout,*) '   Diagnose heat/salt budget or not          ln_limdiahsb  = ', ln_limdiahsb 
     158         WRITE(numout,*) '   Output   heat/salt budget or not          ln_limdiaout  = ', ln_limdiaout 
    159159      ENDIF 
    160160      ! 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limadv.F90

    r3764 r4161  
    1515   !!   lim_adv_y  : advection of sea ice on y axis 
    1616   !!---------------------------------------------------------------------- 
    17    USE dom_oce        ! ocean domain 
    18    USE ice            ! LIM-3 variables 
    19    USE dom_ice        ! LIM-3 domain 
    20    USE lbclnk         ! lateral boundary condition - MPP exchanges 
    21    USE in_out_manager ! I/O manager 
    22    USE prtctl         ! Print control 
    23    USE lib_mpp        ! MPP library 
    24    USE wrk_nemo       ! work arrays 
    25    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    26    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     17   USE dom_oce          ! ocean domain 
     18   USE dom_ice          ! LIM-3 domain 
     19   USE ice              ! LIM-3 variables 
     20   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     21   USE in_out_manager   ! I/O manager 
     22   USE prtctl           ! Print control 
     23   USE lib_mpp          ! MPP library 
     24   USE wrk_nemo         ! work arrays 
     25   USE lib_fortran      ! to use key_nosignedzero 
    2726 
    2827   IMPLICIT NONE 
     
    3938#  include "vectopt_loop_substitute.h90" 
    4039   !!---------------------------------------------------------------------- 
    41    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4241   !! $Id$ 
    4342   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r3625 r4161  
    3030 
    3131   !!---------------------------------------------------------------------- 
    32    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     32   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    3333   !! $Id$ 
    3434   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r4147 r4161  
    1515   !!    lim_dyn_init : initialization and namelist read 
    1616   !!---------------------------------------------------------------------- 
    17    USE phycst         ! physical constants 
    18    USE dom_oce        ! ocean space and time domain 
    19    USE sbc_oce        ! Surface boundary condition: ocean fields 
    20    USE sbc_ice        ! Surface boundary condition: ice   fields 
    21    USE ice            ! LIM-3 variables 
    22    USE par_ice        ! LIM-3 parameters 
    23    USE dom_ice        ! LIM-3 domain 
    24    USE limrhg         ! LIM-3 rheology 
    25    USE lbclnk         ! lateral boundary conditions - MPP exchanges 
    26    USE lib_mpp        ! MPP library 
    27    USE wrk_nemo       ! work arrays 
    28    USE in_out_manager ! I/O manager 
    29    USE prtctl         ! Print control 
    30    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     17   USE phycst           ! physical constants 
     18   USE dom_oce          ! ocean space and time domain 
     19   USE sbc_oce          ! Surface boundary condition: ocean fields 
     20   USE sbc_ice          ! Surface boundary condition: ice   fields 
     21   USE ice              ! LIM-3 variables 
     22   USE par_ice          ! LIM-3 parameters 
     23   USE dom_ice          ! LIM-3 domain 
     24   USE limrhg           ! LIM-3 rheology 
     25   USE lbclnk           ! lateral boundary conditions - MPP exchanges 
     26   USE lib_mpp          ! MPP library 
     27   USE wrk_nemo         ! work arrays 
     28   USE in_out_manager   ! I/O manager 
     29   USE prtctl           ! Print control 
     30   USE lib_fortran      ! glob_sum 
     31   USE timing          ! Timing 
    3132 
    3233   IMPLICIT NONE 
     
    3839#  include "vectopt_loop_substitute.h90" 
    3940   !!---------------------------------------------------------------------- 
    40    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4142   !! $Id$ 
    4243   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    6566      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask 
    6667      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity 
    67       !!--------------------------------------------------------------------- 
     68      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     69      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     70     !!--------------------------------------------------------------------- 
     71 
     72      IF( nn_timing == 1 )  CALL timing_start('limdyn') 
    6873 
    6974      CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 
    7075      CALL wrk_alloc( jpj, zind, zmsk ) 
     76 
     77      ! ------------------------------- 
     78      !- check conservation (C Rousset) 
     79      IF (ln_limdiahsb) THEN 
     80         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     81         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     82         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
     83         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     84      ENDIF 
     85      !- check conservation (C Rousset) 
     86      ! ------------------------------- 
    7187 
    7288      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
     
    208224      ENDIF 
    209225      ! 
     226      ! ------------------------------- 
     227      !- check conservation (C Rousset) 
     228      IF (ln_limdiahsb) THEN 
     229         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     230         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     231  
     232         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) / rdt_ice 
     233         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) / rdt_ice + ( zchk_fs / rhoic ) 
     234 
     235         zchk_vmin = glob_min(v_i) 
     236         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     237         zchk_amin = glob_min(a_i) 
     238 
     239         IF(lwp) THEN 
     240            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limdyn) = ',(zchk_v_i * rday) 
     241            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limdyn) = ',(zchk_smv * rday) 
     242            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limdyn) = ',(zchk_vmin * 1.e-3) 
     243            !IF ( zchk_amax >  amax+1.e-10   ) WRITE(numout,*) 'violation a_i>amax            (limdyn) = ',zchk_amax 
     244            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limdyn) = ',zchk_amin 
     245         ENDIF 
     246      ENDIF 
     247      !- check conservation (C Rousset) 
     248      ! ------------------------------- 
     249 
    210250      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 
    211251      CALL wrk_dealloc( jpj, zind, zmsk ) 
    212252      ! 
     253      IF( nn_timing == 1 )  CALL timing_stop('limdyn') 
     254 
    213255   END SUBROUTINE lim_dyn 
    214256 
     
    230272         &                dm, nbiter, nbitdr, om, resl, cw, angvg, pstar,   & 
    231273         &                c_rhg, etamn, creepl, ecc, ahi0, & 
    232          &                nevp, telast, alphaevp 
     274         &                nevp, telast, alphaevp, hminrhg 
    233275      !!------------------------------------------------------------------- 
    234276 
     
    264306         WRITE(numout,*) '   timescale for elastic waves                      telast = ', telast 
    265307         WRITE(numout,*) '   coefficient for the solution of int. stresses  alphaevp = ', alphaevp 
     308         WRITE(numout,*) '   min ice thickness for rheology calculations     hminrhg = ', hminrhg 
    266309      ENDIF 
    267310      ! 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90

    r3625 r4161  
    3535#  include "vectopt_loop_substitute.h90" 
    3636   !!---------------------------------------------------------------------- 
    37    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
     37   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    3838   !! $Id$ 
    3939   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    137137      END DO                                       ! end of sub-time step loop 
    138138 
     139      ! ----------------------- 
     140      !!! final step (clem) !!! 
     141      DO jj = 1, jpjm1                                ! diffusive fluxes in U- and V- direction 
     142         DO ji = 1 , fs_jpim1   ! vector opt. 
     143            zflu(ji,jj) = pahu(ji,jj) * e2u(ji,jj) / e1u(ji,jj) * ( ptab(ji+1,jj) - ptab(ji,jj) ) 
     144            zflv(ji,jj) = pahv(ji,jj) * e1v(ji,jj) / e2v(ji,jj) * ( ptab(ji,jj+1) - ptab(ji,jj) ) 
     145         END DO 
     146      END DO 
     147      ! 
     148      DO jj= 2, jpjm1                                 ! diffusive trend : divergence of the fluxes 
     149         DO ji = fs_2 , fs_jpim1   ! vector opt.  
     150            zdiv (ji,jj) = (  zflu(ji,jj) - zflu(ji-1,jj  )   & 
     151                 &            + zflv(ji,jj) - zflv(ji  ,jj-1)  ) / ( e1t (ji,jj) * e2t (ji,jj) ) 
     152            ptab(ji,jj) = ztab0(ji,jj) + 0.5 * ( zdiv(ji,jj) + zdiv0(ji,jj) ) 
     153         END DO 
     154      END DO 
     155      CALL lbc_lnk( ptab, 'T', 1. )                   ! lateral boundary condition 
     156      !!! final step (clem) !!! 
     157      ! ----------------------- 
     158 
    139159      IF(ln_ctl)   THEN 
    140160         zrlx(:,:) = ptab(:,:) - ztab0(:,:) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4147 r4161  
    55   !!====================================================================== 
    66   !! History :  2.0  ! 2004-01 (C. Ethe, G. Madec)  Original code 
    7    !!            4.0  ! 2011-02  (G. Madec) dynamical allocation 
     7   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     8   !!             -   ! 2012    (C. Rousset) add par_oce (for jp_sal)...bug? 
    89   !!---------------------------------------------------------------------- 
    910#if defined key_lim3 
     
    1819   USE dom_oce          ! ocean domain 
    1920   USE sbc_oce          ! Surface boundary condition: ocean fields 
     21   USE sbc_ice          ! Surface boundary condition: ice fields 
    2022   USE eosbn2           ! equation of state 
    2123   USE ice              ! sea-ice variables 
    2224   USE par_ice          ! ice parameters 
     25   USE par_oce          ! ocean parameters 
    2326   USE dom_ice          ! sea-ice domain 
    2427   USE in_out_manager   ! I/O manager 
    2528   USE lbclnk           ! lateral boundary condition - MPP exchanges 
    2629   USE lib_mpp          ! MPP library 
     30   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2731   USE wrk_nemo         ! work arrays 
    28    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    2932 
    3033   IMPLICIT NONE 
     
    4952 
    5053   !!---------------------------------------------------------------------- 
    51    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     54   !!   LIM 3.0,  UCL-LOCEAN-IPSL (2008) 
    5255   !! $Id$ 
    53    !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    54    !!---------------------------------------------------------------------- 
     56   !! Software governed by the CeCILL licence (modipsl/doc/NEMO_CeCILL.txt) 
     57   !!---------------------------------------------------------------------- 
     58 
    5559CONTAINS 
    5660 
     
    6165      !! ** Purpose :   defined the sea-ice initial state 
    6266      !! 
    63       !! ** Method  :   restart from a state defined in a binary file 
    64       !!                or from arbitrary sea-ice conditions 
    65       !!------------------------------------------------------------------- 
    66       INTEGER  ::   ji, jj, jk, jl             ! dummy loop indices 
    67       REAL(wp) ::   zeps6, zeps, ztmelts, epsi06   ! local scalars 
    68       REAL(wp) ::   zvol, zare, zh, zh1, zh2, zh3, zan, zbn, zas, zbs  
    69       REAL(wp), POINTER, DIMENSION(:)   ::   zgfactorn, zhin  
    70       REAL(wp), POINTER, DIMENSION(:)   ::   zgfactors, zhis 
    71       REAL(wp), POINTER, DIMENSION(:,:) ::   zidto      ! ice indicator 
    72       !-------------------------------------------------------------------- 
    73  
    74       CALL wrk_alloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 
     67      !! ** Method  :    
     68      !!                This routine will put some ice where ocean 
     69      !!                is at the freezing point, then fill in ice  
     70      !!                state variables using prescribed initial  
     71      !!                values in the namelist             
     72      !! 
     73      !! ** Steps   :    
     74      !!                1) Read namelist 
     75      !!                2) Basal temperature; ice and hemisphere masks 
     76      !!                3) Fill in the ice thickness distribution using gaussian 
     77      !!                4) Fill in space-dependent arrays for state variables 
     78      !!                5) Diagnostic arrays 
     79      !!                6) Lateral boundary conditions 
     80      !! 
     81      !! History : 
     82      !!   2.0  !  01-04  (C. Ethe, G. Madec)  Original code 
     83      !!   3.0  !  2007   (M. Vancoppenolle)   Rewrite for ice cats 
     84      !!   4.0  !  09-11  (M. Vancoppenolle)   Enhanced version for ice cats 
     85      !!-------------------------------------------------------------------- 
     86 
     87      !! * Local variables 
     88      INTEGER    :: ji, jj, jk, jl             ! dummy loop indices 
     89      REAL(wp)   :: epsi06, epsi20, ztmelts 
     90      INTEGER    :: i_hemis, i_fill, jl0   
     91      REAL(wp)   :: ztest_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv 
     92      REAL(wp), POINTER, DIMENSION(:)     :: zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini 
     93      REAL(wp), POINTER, DIMENSION(:,:)   :: zht_i_ini, za_i_ini, zv_i_ini 
     94      REAL(wp), POINTER, DIMENSION(:,:)   :: zidto    ! ice indicator 
     95      INTEGER,  POINTER, DIMENSION(:,:)   :: zhemis   ! hemispheric index 
     96      !-------------------------------------------------------------------- 
     97 
    7598      CALL wrk_alloc( jpi, jpj, zidto ) 
    76  
    77       !-------------------------------------------------------------------- 
    78       ! 1) Preliminary things  
    79       !-------------------------------------------------------------------- 
    80       epsi06 = 1.e-6_wp 
     99      CALL wrk_alloc( jpi, jpj, zhemis ) 
     100      CALL wrk_alloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini ) 
     101      CALL wrk_alloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
     102 
     103      epsi06   = 1.0e-6 
     104      epsi20   = 1.0e-20 
     105      IF(lwp) WRITE(numout,*) 
     106      IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 
     107      IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 
     108 
     109      !-------------------------------------------------------------------- 
     110      ! 1) Read namelist 
     111      !-------------------------------------------------------------------- 
    81112 
    82113      CALL lim_istate_init     !  reading the initials parameters of the ice 
     
    87118 
    88119      !-------------------------------------------------------------------- 
    89       ! 2) Ice initialization (hi,hs,frld,t_su,sm_i,t_i,t_s)              |  
    90       !-------------------------------------------------------------------- 
    91  
    92       IF(lwp) WRITE(numout,*) 
    93       IF(lwp) WRITE(numout,*) 'lim_istate : Ice initialization ' 
    94       IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ' 
    95  
     120      ! 2) Basal temperature, ice mask and hemispheric index 
     121      !-------------------------------------------------------------------- 
     122 
     123      ! Basal temperature is set to the freezing point of seawater in Celsius 
    96124      t_bo(:,:) = tfreez( tsn(:,:,1,jp_sal) ) * tmask(:,:,1)       ! freezing/melting point of sea water [Celcius] 
    97125 
    98126      DO jj = 1, jpj                                       ! ice if sst <= t-freez + ttest 
    99127         DO ji = 1, jpi 
    100             IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0.e0      ! no ice 
    101             ELSE                                                     ;   zidto(ji,jj) = 1.e0      !    ice 
     128            IF( tsn(ji,jj,1,jp_tem)  - t_bo(ji,jj) >= ttest ) THEN   ;   zidto(ji,jj) = 0._wp      ! no ice 
     129            ELSE                                                     ;   zidto(ji,jj) = 1._wp      !    ice 
    102130            ENDIF 
    103131         END DO 
    104132      END DO 
    105133 
    106       t_bo(:,:) = t_bo(:,:) + rt0                          ! t_bo converted from Celsius to Kelvin (rt0 over land) 
    107  
    108       ! constants for heat contents 
    109       zeps   = 1.e-20_wp 
    110       zeps6  = 1.e-06_wp 
    111  
    112       ! zgfactor for initial ice distribution 
    113       zgfactorn(:) = 0._wp 
    114       zgfactors(:) = 0._wp 
    115  
    116       ! first ice type 
    117       DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) 
    118          zhin (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 
    119          zgfactorn(1) = zgfactorn(1) + exp(-(zhin(1)-hginn_u)*(zhin(1)-hginn_u) * 0.5_wp ) 
    120          zhis (1)     = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 
    121          zgfactors(1) = zgfactors(1) + exp(-(zhis(1)-hgins_u)*(zhis(1)-hgins_u) * 0.5_wp ) 
    122       END DO ! jl 
    123       zgfactorn(1) = aginn_u / zgfactorn(1) 
    124       zgfactors(1) = agins_u / zgfactors(1) 
    125  
    126       ! ------------- 
    127       ! new distribution, polynom of second order, conserving area and volume 
    128       zh1 = 0._wp 
    129       zh2 = 0._wp 
    130       zh3 = 0._wp 
    131       DO jl = 1, jpl 
    132          zh = ( hi_max(jl-1) + hi_max(jl) ) * 0.5_wp 
    133          zh1 = zh1 + zh 
    134          zh2 = zh2 + zh * zh 
    135          zh3 = zh3 + zh * zh * zh 
    136       END DO 
    137       IF(lwp) WRITE(numout,*) ' zh1 : ', zh1 
    138       IF(lwp) WRITE(numout,*) ' zh2 : ', zh2 
    139       IF(lwp) WRITE(numout,*) ' zh3 : ', zh3 
    140  
    141       zvol = aginn_u * hginn_u 
    142       zare = aginn_u 
    143       IF( jpl >= 2 ) THEN 
    144          zbn = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 
    145          zan = ( zare - zbn*zh1 ) / zh2 
    146       ENDIF 
    147  
    148       IF(lwp) WRITE(numout,*) ' zvol: ', zvol 
    149       IF(lwp) WRITE(numout,*) ' zare: ', zare 
    150       IF(lwp) WRITE(numout,*) ' zbn : ', zbn  
    151       IF(lwp) WRITE(numout,*) ' zan : ', zan  
    152  
    153       zvol = agins_u * hgins_u 
    154       zare = agins_u 
    155       IF( jpl >= 2 ) THEN 
    156          zbs = ( zvol*zh2 - zare*zh3 ) / ( zh2*zh2 - zh1*zh3) 
    157          zas = ( zare - zbs*zh1 ) / zh2 
    158       ENDIF 
    159  
    160       IF(lwp) WRITE(numout,*) ' zvol: ', zvol 
    161       IF(lwp) WRITE(numout,*) ' zare: ', zare 
    162       IF(lwp) WRITE(numout,*) ' zbn : ', zbn  
    163       IF(lwp) WRITE(numout,*) ' zan : ', zan  
    164  
    165       !end of new lines 
    166       ! ------------- 
    167 !!! 
    168       ! retour a LIMA_MEC 
    169       !     ! second ice type 
    170       !     zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    171       !     hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    172  
    173       !     ! here to change !!!! 
    174       !     jm = 2 
    175       !     DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    176       !        zhin (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    177       !        zhin (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
    178       !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
    179       !        zgfactorn(2) = zgfactorn(2) + exp(-(zhin(2)-hginn_d)*(zhin(2)-hginn_d)/2.0) 
    180       !        zhis (2)     = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    181       !        zhis (2)     = ( hi_max_typ(jl-ice_cat_bounds(2,1),jm    ) + & 
    182       !                         hi_max_typ(jl-ice_cat_bounds(2,1) + 1,jm)   ) / 2.0 
    183       !        zgfactors(2) = zgfactors(2) + exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0) 
    184       !     END DO ! jl 
    185       !     zgfactorn(2) = aginn_d / zgfactorn(2) 
    186       !     zgfactors(2) = agins_d / zgfactors(2) 
    187       !     hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    188       ! END retour a LIMA_MEC 
    189 !!! 
    190  
    191 !!gm  optimisation :  loop over the ice categories inside the ji, jj loop !!! 
    192  
     134      t_bo(:,:) = t_bo(:,:) + rt0                          ! conversion to Kelvin 
     135 
     136      ! Hemispheric index 
     137      ! MV 2011 new initialization 
    193138      DO jj = 1, jpj 
    194139         DO ji = 1, jpi 
    195  
    196             !--- Northern hemisphere 
    197             !---------------------------------------------------------------- 
    198140            IF( fcor(ji,jj) >= 0._wp ) THEN     
    199  
    200                !----------------------- 
    201                ! Ice area / thickness 
    202                !----------------------- 
    203  
    204                IF ( jpl .EQ. 1) THEN ! one category 
    205  
    206                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 
    207                      a_i(ji,jj,jl)    = zidto(ji,jj) * aginn_u 
    208                      ht_i(ji,jj,jl)   = zidto(ji,jj) * hginn_u 
    209                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    210                   END DO 
    211  
    212                ELSE ! several categories 
    213  
    214                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 
    215                      zhin(1)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    216                      a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(1) * exp(-(zhin(1)-hginn_u)* &  
    217                         (zhin(1)-hginn_u)/2.0) , epsi06) 
    218                      ! new line 
    219                      a_i(ji,jj,jl)    = zidto(ji,jj) * ( zan * zhin(1) * zhin(1) + zbn * zhin(1) ) 
    220                      ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(1)  
    221                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    222                   END DO 
    223  
    224                ENDIF 
    225  
    226  
    227 !!! 
    228                ! retour a LIMA_MEC 
    229                !              !ridged ice 
    230                !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    231                !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    232                !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) ! loop over ice thickness categories 
    233                !                 zhin(2)          = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    234                !                 a_i(ji,jj,jl)    = zidto(ji,jj) * MAX( zgfactorn(2) * exp(-(zhin(2)-hginn_d)* & 
    235                !                                    (zhin(2)-hginn_d)/2.0) , epsi06) 
    236                !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhin(2)  
    237                !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    238                !              END DO 
    239                !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    240  
    241                !              !rafted ice 
    242                !              jl = 6 
    243                !              a_i(ji,jj,jl)       = 0.0 
    244                !              ht_i(ji,jj,jl)      = 0.0 
    245                !              v_i(ji,jj,jl)       = 0.0 
    246                ! END retour a LIMA_MEC 
    247 !!! 
    248  
    249                DO jl = 1, jpl 
    250  
    251                   !------------- 
    252                   ! Snow depth 
    253                   !------------- 
    254                   ht_s(ji,jj,jl)   = zidto(ji,jj) * hninn 
    255                   v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 
    256  
    257                   !--------------- 
    258                   ! Ice salinity 
    259                   !--------------- 
    260                   sm_i(ji,jj,jl)   = zidto(ji,jj) * sinn  + ( 1.0 - zidto(ji,jj) ) * 0.1 
    261                   smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    262  
    263                   !---------- 
    264                   ! Ice age 
    265                   !---------- 
    266                   o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) ) 
    267                   oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    268  
    269                   !------------------------------ 
    270                   ! Sea ice surface temperature 
    271                   !------------------------------ 
    272  
    273                   t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 
    274  
    275                   !------------------------------------ 
    276                   ! Snow temperature and heat content 
    277                   !------------------------------------ 
    278  
    279                   DO jk = 1, nlay_s 
    280                      t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
    281                      ! Snow energy of melting 
    282                      e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    283                      ! Change dimensions 
    284                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    285                      ! Multiply by volume, so that heat content in 10^9 Joules 
    286                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    287                         v_s(ji,jj,jl)  / nlay_s 
    288                   END DO !jk 
    289  
    290                   !----------------------------------------------- 
    291                   ! Ice salinities, temperature and heat content  
    292                   !----------------------------------------------- 
    293  
    294                   DO jk = 1, nlay_i 
    295                      t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt  
    296                      s_i(ji,jj,jk,jl) = zidto(ji,jj) * sinn + ( 1.0 - zidto(ji,jj) ) * 0.1 
    297                      ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    298  
    299                      ! heat content per unit volume 
    300                      e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 
    301                         (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    302                         +   lfus    * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
    303                         - rcp      * ( ztmelts - rtt ) & 
    304                         ) 
    305  
    306                      ! Correct dimensions to avoid big values 
    307                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    308  
    309                      ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    310                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    311                         area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
    312                         nlay_i 
    313                   END DO ! jk 
    314  
    315                END DO ! jl  
    316  
    317             ELSE ! on fcor  
    318  
    319                !--- Southern hemisphere 
    320                !---------------------------------------------------------------- 
    321  
    322                !----------------------- 
    323                ! Ice area / thickness 
    324                !----------------------- 
    325  
    326                IF ( jpl .EQ. 1) THEN ! one category 
    327  
    328                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) ! loop over ice thickness categories 
    329                      a_i(ji,jj,jl)    = zidto(ji,jj) * agins_u 
    330                      ht_i(ji,jj,jl)   = zidto(ji,jj) * hgins_u 
    331                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    332                   END DO 
    333  
    334                ELSE ! several categories 
    335  
    336                   !level ice 
    337                   DO jl = ice_cat_bounds(1,1), ice_cat_bounds(1,2) !over thickness categories 
    338  
    339                      zhis(1)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    340                      a_i(ji,jj,jl) = zidto(ji,jj) * MAX( zgfactors(1) * exp(-(zhis(1)-hgins_u) * &  
    341                         (zhis(1)-hgins_u)/2.0) , epsi06 ) 
    342                      ! new line square distribution volume conserving 
    343                      a_i(ji,jj,jl)    = zidto(ji,jj) * ( zas * zhis(1) * zhis(1) + zbs * zhis(1) ) 
    344                      ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(1)  
    345                      v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    346  
    347                   END DO ! jl 
    348  
    349                ENDIF 
    350  
    351 !!! 
    352                ! retour a LIMA_MEC 
    353                !              !ridged ice 
    354                !              zdummy  = hi_max(ice_cat_bounds(2,1)-1) 
    355                !              hi_max(ice_cat_bounds(2,1)-1) = 0.0 
    356                !              DO jl = ice_cat_bounds(2,1), ice_cat_bounds(2,2) !over thickness categories 
    357                !                 zhis(2)       = ( hi_max(jl-1) + hi_max(jl) ) / 2.0 
    358                !                 a_i(ji,jj,jl) = zidto(ji,jj)*MAX( zgfactors(2)   & 
    359                !                    &          * exp(-(zhis(2)-hgins_d)*(zhis(2)-hgins_d)/2.0), epsi06 ) 
    360                !                 ht_i(ji,jj,jl)   = zidto(ji,jj) * zhis(2)  
    361                !                 v_i(ji,jj,jl)    = ht_i(ji,jj,jl)*a_i(ji,jj,jl) 
    362                !              END DO 
    363                !              hi_max(ice_cat_bounds(2,1)-1) = zdummy 
    364  
    365                !              !rafted ice 
    366                !              jl = 6 
    367                !              a_i(ji,jj,jl)       = 0.0 
    368                !              ht_i(ji,jj,jl)      = 0.0 
    369                !              v_i(ji,jj,jl)       = 0.0 
    370                ! END retour a LIMA_MEC 
    371 !!! 
    372  
    373                DO jl = 1, jpl !over thickness categories 
    374  
    375                   !--------------- 
    376                   ! Snow depth 
    377                   !--------------- 
    378  
    379                   ht_s(ji,jj,jl)   = zidto(ji,jj) * hnins 
    380                   v_s(ji,jj,jl)    = ht_s(ji,jj,jl)*a_i(ji,jj,jl) 
    381  
    382                   !--------------- 
    383                   ! Ice salinity 
    384                   !--------------- 
    385  
    386                   sm_i(ji,jj,jl)   = zidto(ji,jj) * sins  + ( 1.0 - zidto(ji,jj) ) * 0.1 
    387                   smv_i(ji,jj,jl)  = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) 
    388  
    389                   !---------- 
    390                   ! Ice age 
    391                   !---------- 
    392  
    393                   o_i(ji,jj,jl)    = zidto(ji,jj) * 1.0   + ( 1.0 - zidto(ji,jj) ) 
    394                   oa_i(ji,jj,jl)   = o_i(ji,jj,jl) * a_i(ji,jj,jl) 
    395  
    396                   !------------------------------ 
    397                   ! Sea ice surface temperature 
    398                   !------------------------------ 
    399  
    400                   t_su(ji,jj,jl)   = zidto(ji,jj) * 270.0 + ( 1.0 - zidto(ji,jj) ) * t_bo(ji,jj) 
    401  
    402                   !---------------------------------- 
    403                   ! Snow temperature / heat content 
    404                   !---------------------------------- 
    405  
    406                   DO jk = 1, nlay_s 
    407                      t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1.0 - zidto(ji,jj) ) * rtt 
    408                      ! Snow energy of melting 
    409                      e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
    410                      ! Change dimensions 
    411                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
    412                      ! Multiply by volume, so that heat content in 10^9 Joules 
    413                      e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * & 
    414                         v_s(ji,jj,jl)  / nlay_s 
    415                   END DO 
    416  
    417                   !--------------------------------------------- 
    418                   ! Ice temperature, salinity and heat content 
    419                   !--------------------------------------------- 
    420  
    421                   DO jk = 1, nlay_i 
    422                      t_i(ji,jj,jk,jl) = zidto(ji,jj)*270.00 + ( 1.0 - zidto(ji,jj) ) * rtt  
    423                      s_i(ji,jj,jk,jl) = zidto(ji,jj) * sins + ( 1.0 - zidto(ji,jj) ) * 0.1 
    424                      ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
    425  
    426                      ! heat content per unit volume 
    427                      e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * & 
    428                         (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
    429                         +   lfus  * ( 1.0 - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-zeps) ) & 
    430                         - rcp      * ( ztmelts - rtt ) & 
    431                         ) 
    432  
    433                      ! Correct dimensions to avoid big values 
    434                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
    435  
    436                      ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    437                      e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * &  
    438                         area(ji,jj) * a_i(ji,jj,jl) * ht_i(ji,jj,jl) / & 
    439                         nlay_i 
    440                   END DO !jk 
    441  
    442                END DO ! jl  
    443  
    444             ENDIF ! on fcor 
    445  
     141               zhemis(ji,jj) = 1 ! Northern hemisphere 
     142            ELSE 
     143               zhemis(ji,jj) = 2 ! Southern hemisphere 
     144            ENDIF 
    446145         END DO 
    447146      END DO 
    448  
    449       !-------------------------------------------------------------------- 
    450       ! 3) Global ice variables for output diagnostics                    |  
    451       !-------------------------------------------------------------------- 
    452  
    453       fsbbq (:,:)     = 0.e0 
    454       u_ice (:,:)     = 0.e0 
    455       v_ice (:,:)     = 0.e0 
    456       stress1_i(:,:)  = 0.0 
    457       stress2_i(:,:)  = 0.0 
    458       stress12_i(:,:) = 0.0 
    459  
    460       !-------------------------------------------------------------------- 
    461       ! 4) Moments for advection 
    462       !-------------------------------------------------------------------- 
    463  
    464       sxopw (:,:) = 0.e0  
    465       syopw (:,:) = 0.e0  
    466       sxxopw(:,:) = 0.e0  
    467       syyopw(:,:) = 0.e0  
    468       sxyopw(:,:) = 0.e0 
    469  
    470       sxice (:,:,:)  = 0.e0   ;   sxsn (:,:,:)  = 0.e0   ;   sxa  (:,:,:)  = 0.e0 
    471       syice (:,:,:)  = 0.e0   ;   sysn (:,:,:)  = 0.e0   ;   sya  (:,:,:)  = 0.e0 
    472       sxxice(:,:,:)  = 0.e0   ;   sxxsn(:,:,:)  = 0.e0   ;   sxxa (:,:,:)  = 0.e0 
    473       syyice(:,:,:)  = 0.e0   ;   syysn(:,:,:)  = 0.e0   ;   syya (:,:,:)  = 0.e0 
    474       sxyice(:,:,:)  = 0.e0   ;   sxysn(:,:,:)  = 0.e0   ;   sxya (:,:,:)  = 0.e0 
    475  
    476       sxc0  (:,:,:)  = 0.e0   ;   sxe  (:,:,:,:)= 0.e0    
    477       syc0  (:,:,:)  = 0.e0   ;   sye  (:,:,:,:)= 0.e0    
    478       sxxc0 (:,:,:)  = 0.e0   ;   sxxe (:,:,:,:)= 0.e0    
    479       syyc0 (:,:,:)  = 0.e0   ;   syye (:,:,:,:)= 0.e0    
    480       sxyc0 (:,:,:)  = 0.e0   ;   sxye (:,:,:,:)= 0.e0    
    481  
    482       sxsal  (:,:,:)  = 0.e0 
    483       sysal  (:,:,:)  = 0.e0 
    484       sxxsal (:,:,:)  = 0.e0 
    485       syysal (:,:,:)  = 0.e0 
    486       sxysal (:,:,:)  = 0.e0 
    487  
    488       !-------------------------------------------------------------------- 
    489       ! 5) Lateral boundary conditions                                    |  
     147      ! END MV 2011 new initialization 
     148 
     149      !-------------------------------------------------------------------- 
     150      ! 3) Initialization of sea ice state variables 
     151      !-------------------------------------------------------------------- 
     152 
     153      !----------------------------- 
     154      ! 3.1) Hemisphere-dependent arrays 
     155      !----------------------------- 
     156      ! assign initial thickness, concentration, snow depth and salinity to 
     157      ! an hemisphere-dependent array 
     158      zhm_i_ini(1) = hginn ; zhm_i_ini(2) = hgins  ! ice thickness 
     159      zat_i_ini(1) = aginn ; zat_i_ini(2) = agins  ! ice concentration 
     160      zvt_i_ini(:) = zhm_i_ini(:) * zat_i_ini(:)   ! ice volume 
     161      zhm_s_ini(1) = hninn ; zhm_s_ini(2) = hnins  ! snow depth 
     162      zsm_i_ini(1) = sinn  ; zsm_i_ini(2) = sins   ! bulk ice salinity 
     163 
     164      !--------------------------------------------------------------------- 
     165      ! 3.2) Distribute ice concentration and thickness into the categories 
     166      !--------------------------------------------------------------------- 
     167      ! a gaussian distribution for ice concentration is used 
     168      ! then we check whether the distribution fullfills 
     169      ! volume and area conservation, positivity and ice categories bounds 
     170      DO i_hemis = 1, 2 
     171 
     172      ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 
     173 
     174      ! note for the great nemo engineers:  
     175      ! only very few of the WRITE statements are necessary for the reference version 
     176      ! they were one day useful, but now i personally doubt of their 
     177      ! potential for bringing anything useful 
     178 
     179      DO i_fill = jpl, 1, -1 
     180         IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) .NE. 4 ) THEN 
     181            !---------------------------- 
     182            ! fill the i_fill categories 
     183            !---------------------------- 
     184            ! *** 1 category to fill 
     185            IF ( i_fill .EQ. 1 ) THEN 
     186               zht_i_ini(1,i_hemis)       = zhm_i_ini(i_hemis) 
     187               za_i_ini(1,i_hemis)        = zat_i_ini(i_hemis) 
     188               zht_i_ini(2:jpl,i_hemis)   = 0._wp 
     189               za_i_ini(2:jpl,i_hemis)    = 0._wp 
     190            ELSE 
     191 
     192            ! *** >1 categores to fill 
     193            !--- Ice thicknesses in the i_fill - 1 first categories 
     194               DO jl = 1, i_fill - 1 
     195                  zht_i_ini(jl,i_hemis)    = 0.5 * ( hi_max(jl) + hi_max(jl-1) ) 
     196               END DO 
     197 
     198            !--- jl0: most likely index where cc will be maximum 
     199               DO jl = 1, jpl 
     200                  IF ( ( zhm_i_ini(i_hemis) .GT. hi_max(jl-1) ) .AND. & 
     201                       ( zhm_i_ini(i_hemis) .LE. hi_max(jl)   ) ) THEN 
     202                     jl0 = jl 
     203                  ENDIF 
     204               END DO 
     205               jl0 = MIN(jl0, i_fill) 
     206 
     207            !--- Concentrations 
     208               za_i_ini(jl0,i_hemis)      = zat_i_ini(i_hemis) / SQRT(REAL(jpl)) 
     209               DO jl = 1, i_fill - 1 
     210                  IF ( jl .NE. jl0 ) THEN 
     211                     zsigma               = 0.5 * zhm_i_ini(i_hemis) 
     212                     zarg                 = ( zht_i_ini(jl,i_hemis) - zhm_i_ini(i_hemis) ) / zsigma 
     213                     za_i_ini(jl,i_hemis) = za_i_ini(jl0,i_hemis) * EXP(-zarg**2) 
     214                  ENDIF 
     215               END DO  
     216 
     217               zA = 0. ! sum of the areas in the jpl categories  
     218               DO jl = 1, i_fill - 1 
     219                 zA = zA + za_i_ini(jl,i_hemis) 
     220               END DO 
     221               za_i_ini(i_fill,i_hemis)   = zat_i_ini(i_hemis) - zA ! ice conc in the last category 
     222               IF ( i_fill .LT. jpl ) za_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     223          
     224            !--- Ice thickness in the last category 
     225               zV = 0. ! sum of the volumes of the N-1 categories 
     226               DO jl = 1, i_fill - 1 
     227                  zV = zV + za_i_ini(jl,i_hemis)*zht_i_ini(jl,i_hemis) 
     228               END DO 
     229               zht_i_ini(i_fill,i_hemis) = ( zvt_i_ini(i_hemis) - zV ) / za_i_ini(i_fill,i_hemis)  
     230               IF ( i_fill .LT. jpl ) zht_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     231 
     232            !--- volumes 
     233               zv_i_ini(:,i_hemis) = za_i_ini(:,i_hemis) * zht_i_ini(:,i_hemis) 
     234               IF ( i_fill .LT. jpl ) zv_i_ini(i_fill+1:jpl, i_hemis) = 0._wp 
     235 
     236            ENDIF ! i_fill 
     237 
     238            !--------------------- 
     239            ! Compatibility tests 
     240            !--------------------- 
     241            ! Test 1: area conservation 
     242            zA_cons = SUM(za_i_ini(:,i_hemis)) ; zconv = ABS(zat_i_ini(i_hemis) - zA_cons ) 
     243            IF ( zconv .LT. 1.0e-6 ) THEN 
     244               ztest_1 = 1 
     245            ELSE  
     246              ! this write is useful 
     247              IF(lwp)  WRITE(numout,*) ' * TEST1 AREA NOT CONSERVED *** zA_cons = ', zA_cons,' zat_i_ini = ',zat_i_ini(i_hemis)  
     248               ztest_1 = 0 
     249            ENDIF 
     250 
     251            ! Test 2: volume conservation 
     252            zV_cons = SUM(zv_i_ini(:,i_hemis)) 
     253            zconv = ABS(zvt_i_ini(i_hemis) - zV_cons) 
     254 
     255            IF ( zconv .LT. 1.0e-6 ) THEN 
     256               ztest_2 = 1 
     257            ELSE 
     258              ! this write is useful 
     259              IF(lwp)  WRITE(numout,*) ' * TEST2 VOLUME NOT CONSERVED *** zV_cons = ', zV_cons, & 
     260                            ' zvt_i_ini = ', zvt_i_ini(i_hemis) 
     261               ztest_2 = 0 
     262            ENDIF 
     263 
     264            ! Test 3: thickness of the last category is in-bounds ? 
     265            IF ( zht_i_ini(i_fill, i_hemis) .GT. hi_max(i_fill-1) ) THEN 
     266               ztest_3 = 1 
     267            ELSE 
     268               ! this write is useful 
     269               IF(lwp) WRITE(numout,*) ' * TEST 3 THICKNESS OF THE LAST CATEGORY OUT OF BOUNDS *** zht_i_ini(i_fill,i_hemis) = ', & 
     270               zht_i_ini(i_fill,i_hemis), ' hi_max(jpl-1) = ', hi_max(i_fill-1) 
     271               ztest_3 = 0 
     272            ENDIF 
     273 
     274            ! Test 4: positivity of ice concentrations 
     275            ztest_4 = 1 
     276            DO jl = 1, jpl 
     277               IF ( za_i_ini(jl,i_hemis) .LT. 0._wp ) THEN  
     278                  ! this write is useful 
     279                  IF(lwp) WRITE(numout,*) ' * TEST 4 POSITIVITY NOT OK FOR CAT ', jl, ' WITH A = ', za_i_ini(jl,i_hemis) 
     280                  ztest_4 = 0 
     281               ENDIF 
     282            END DO 
     283 
     284         ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 
     285  
     286         ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 
     287 
     288      END DO ! i_fill 
     289 
     290      IF(lwp) THEN  
     291         WRITE(numout,*), ' ztests : ', ztests 
     292         IF ( ztests .NE. 4 ) THEN 
     293            WRITE(numout,*) 
     294            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     295            WRITE(numout,*), ' !!!! RED ALERT                  !!! ' 
     296            WRITE(numout,*), ' !!!! BIIIIP BIIIP BIIIIP BIIIIP !!!' 
     297            WRITE(numout,*), ' !!!! Something is wrong in the LIM3 initialization procedure ' 
     298            WRITE(numout,*), ' !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! ' 
     299            WRITE(numout,*) 
     300            WRITE(numout,*), ' *** ztests is not equal to 4 ' 
     301            WRITE(numout,*), ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 
     302            WRITE(numout,*), ' zat_i_ini : ', zat_i_ini(i_hemis) 
     303            WRITE(numout,*), ' zhm_i_ini : ', zhm_i_ini(i_hemis) 
     304         ENDIF ! ztests .NE. 4 
     305      ENDIF 
     306       
     307      END DO ! i_hemis 
     308 
     309      !--------------------------------------------------------------------- 
     310      ! 3.3) Space-dependent arrays for ice state variables 
     311      !--------------------------------------------------------------------- 
     312 
     313      ! Ice concentration, thickness and volume, snow depth, ice 
     314      ! salinity, ice age, surface temperature 
     315      DO jl = 1, jpl ! loop over categories 
     316         DO jj = 1, jpj 
     317            DO ji = 1, jpi 
     318               a_i(ji,jj,jl)   = zidto(ji,jj) * za_i_ini (jl,zhemis(ji,jj))  ! concentration 
     319               ht_i(ji,jj,jl)  = zidto(ji,jj) * zht_i_ini(jl,zhemis(ji,jj))  ! ice thickness 
     320               ht_s(ji,jj,jl)  = zidto(ji,jj) * zhm_s_ini(zhemis(ji,jj))     ! snow depth 
     321               sm_i(ji,jj,jl)  = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj))     ! salinity 
     322               o_i(ji,jj,jl)   = zidto(ji,jj) * 1._wp + ( 1._wp - zidto(ji,jj) ) ! age 
     323               t_su(ji,jj,jl)  = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * t_bo(ji,jj) ! surf temp 
     324  
     325               ! ice volume, snow volume, salt content, age content 
     326               v_i(ji,jj,jl)   = ht_i(ji,jj,jl) * a_i(ji,jj,jl)              ! ice volume 
     327               v_s(ji,jj,jl)   = ht_s(ji,jj,jl) * a_i(ji,jj,jl)              ! snow volume 
     328               smv_i(ji,jj,jl) = MIN( sm_i(ji,jj,jl) , sss_m(ji,jj) ) * v_i(ji,jj,jl) ! salt content 
     329               oa_i(ji,jj,jl)  = o_i(ji,jj,jl) * a_i(ji,jj,jl)               ! age content 
     330            END DO ! ji 
     331         END DO ! jj 
     332      END DO ! jl 
     333 
     334      ! Snow temperature and heat content 
     335      DO jk = 1, nlay_s 
     336         DO jl = 1, jpl ! loop over categories 
     337            DO jj = 1, jpj 
     338               DO ji = 1, jpi 
     339                   t_s(ji,jj,jk,jl) = zidto(ji,jj) * 270.0 + ( 1._wp - zidto(ji,jj) ) * rtt 
     340                   ! Snow energy of melting 
     341                   e_s(ji,jj,jk,jl) = zidto(ji,jj) * rhosn * ( cpic * ( rtt - t_s(ji,jj,jk,jl) ) + lfus ) 
     342                   ! Change dimensions 
     343                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / unit_fac 
     344                   ! Multiply by volume, so that heat content in 10^9 Joules 
     345                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) * area(ji,jj) * v_s(ji,jj,jl) / nlay_s 
     346               END DO ! ji 
     347            END DO ! jj 
     348         END DO ! jl 
     349      END DO ! jk 
     350 
     351      ! Ice salinity, temperature and heat content 
     352      DO jk = 1, nlay_i 
     353         DO jl = 1, jpl ! loop over categories 
     354            DO jj = 1, jpj 
     355               DO ji = 1, jpi 
     356                   t_i(ji,jj,jk,jl) = zidto(ji,jj) * 270.00 + ( 1._wp - zidto(ji,jj) ) * rtt  
     357                   s_i(ji,jj,jk,jl) = zidto(ji,jj) * zsm_i_ini(zhemis(ji,jj)) + ( 1._wp - zidto(ji,jj) ) * s_i_min 
     358                   ztmelts          = - tmut * s_i(ji,jj,jk,jl) + rtt !Melting temperature in K 
     359 
     360                   ! heat content per unit volume 
     361                   e_i(ji,jj,jk,jl) = zidto(ji,jj) * rhoic * (   cpic    * ( ztmelts - t_i(ji,jj,jk,jl) ) & 
     362                      +   lfus    * ( 1._wp - (ztmelts-rtt) / MIN((t_i(ji,jj,jk,jl)-rtt),-epsi20) ) & 
     363                      -   rcp     * ( ztmelts - rtt ) ) 
     364 
     365                   ! Correct dimensions to avoid big values 
     366                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / unit_fac  
     367 
     368                   ! Mutliply by ice volume, and divide by number of layers  
     369                   ! to get heat content in 10^9 J 
     370                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / nlay_i 
     371               END DO ! ji 
     372            END DO ! jj 
     373         END DO ! jl 
     374      END DO ! jk 
     375 
     376      !-------------------------------------------------------------------- 
     377      ! 4) Global ice variables for output diagnostics                    |  
     378      !-------------------------------------------------------------------- 
     379      fsbbq (:,:)     = 0._wp 
     380      u_ice (:,:)     = 0._wp 
     381      v_ice (:,:)     = 0._wp 
     382      stress1_i(:,:)  = 0._wp 
     383      stress2_i(:,:)  = 0._wp 
     384      stress12_i(:,:) = 0._wp 
     385 
     386# if defined key_coupled 
     387      albege(:,:)   = 0.8 * tms(:,:) 
     388# endif 
     389 
     390      !-------------------------------------------------------------------- 
     391      ! 5) Moments for advection 
     392      !-------------------------------------------------------------------- 
     393 
     394      sxopw (:,:) = 0._wp  
     395      syopw (:,:) = 0._wp  
     396      sxxopw(:,:) = 0._wp  
     397      syyopw(:,:) = 0._wp  
     398      sxyopw(:,:) = 0._wp 
     399 
     400      sxice (:,:,:)  = 0._wp   ;   sxsn (:,:,:)  = 0._wp   ;   sxa  (:,:,:)  = 0._wp 
     401      syice (:,:,:)  = 0._wp   ;   sysn (:,:,:)  = 0._wp   ;   sya  (:,:,:)  = 0._wp 
     402      sxxice(:,:,:)  = 0._wp   ;   sxxsn(:,:,:)  = 0._wp   ;   sxxa (:,:,:)  = 0._wp 
     403      syyice(:,:,:)  = 0._wp   ;   syysn(:,:,:)  = 0._wp   ;   syya (:,:,:)  = 0._wp 
     404      sxyice(:,:,:)  = 0._wp   ;   sxysn(:,:,:)  = 0._wp   ;   sxya (:,:,:)  = 0._wp 
     405 
     406      sxc0  (:,:,:)  = 0._wp   ;   sxe  (:,:,:,:)= 0._wp    
     407      syc0  (:,:,:)  = 0._wp   ;   sye  (:,:,:,:)= 0._wp    
     408      sxxc0 (:,:,:)  = 0._wp   ;   sxxe (:,:,:,:)= 0._wp    
     409      syyc0 (:,:,:)  = 0._wp   ;   syye (:,:,:,:)= 0._wp    
     410      sxyc0 (:,:,:)  = 0._wp   ;   sxye (:,:,:,:)= 0._wp    
     411 
     412      sxsal  (:,:,:)  = 0._wp 
     413      sysal  (:,:,:)  = 0._wp 
     414      sxxsal (:,:,:)  = 0._wp 
     415      syysal (:,:,:)  = 0._wp 
     416      sxysal (:,:,:)  = 0._wp 
     417 
     418      !-------------------------------------------------------------------- 
     419      ! 6) Lateral boundary conditions                                    |  
    490420      !-------------------------------------------------------------------- 
    491421 
    492422      DO jl = 1, jpl 
     423 
    493424         CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. ) 
    494425         CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. ) 
     
    496427         CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) 
    497428         CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) 
    498          ! 
     429 
    499430         CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) 
    500431         CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) 
     
    513444         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) 
    514445      END DO 
     446       
     447      at_i (:,:) = 0.0_wp 
     448      DO jl = 1, jpl 
     449         at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 
     450      END DO 
    515451 
    516452      CALL lbc_lnk( at_i , 'T', 1. ) 
     
    519455      CALL lbc_lnk( fsbbq  , 'T', 1. ) 
    520456      ! 
    521       CALL wrk_dealloc( jpm, zgfactorn, zgfactors, zhin, zhis ) 
     457      !-------------------------------------------------------------------- 
     458      ! 6) ????                                                           |  
     459      !-------------------------------------------------------------------- 
     460      tn_ice (:,:,:) = t_su (:,:,:) 
     461 
    522462      CALL wrk_dealloc( jpi, jpj, zidto ) 
    523       ! 
     463      CALL wrk_dealloc( jpi, jpj, zhemis ) 
     464      CALL wrk_dealloc( jpl,   2, zht_i_ini,  za_i_ini,  zv_i_ini ) 
     465      CALL wrk_dealloc(   2,      zhm_i_ini, zat_i_ini, zvt_i_ini, zhm_s_ini, zsm_i_ini ) 
     466 
    524467   END SUBROUTINE lim_istate 
    525  
    526468 
    527469   SUBROUTINE lim_istate_init 
     
    531473      !! ** Purpose : Definition of initial state of the ice  
    532474      !! 
    533       !! ** Method :   Read the namiceini namelist and check the parameter  
    534       !!             values called at the first timestep (nit000) 
    535       !! 
    536       !! ** input  :   namelist namiceini 
     475      !! ** Method : Read the namiceini namelist and check the parameter  
     476      !!       values called at the first timestep (nit000) 
     477      !! 
     478      !! ** input :  
     479      !!        Namelist namiceini 
     480      !! 
     481      !! history : 
     482      !!  8.5  ! 03-08 (C. Ethe) original code  
     483      !!  8.5  ! 07-11 (M. Vancoppenolle) rewritten initialization 
    537484      !!----------------------------------------------------------------------------- 
     485<<<<<<< .courant 
    538486      INTEGER :: ios                 ! Local integer output status for namelist read 
    539487      NAMELIST/namiceini/ ttest, hninn, hginn_u, aginn_u, hginn_d, aginn_d, hnins,   & 
    540488         &                hgins_u, agins_u, hgins_d, agins_d, sinn, sins 
     489======= 
     490      NAMELIST/namiceini/ ttest, hninn, hnins, hginn, hgins, aginn, agins, sinn, sins 
     491>>>>>>> .fusion-droit.r4160 
    541492      !!----------------------------------------------------------------------------- 
     493<<<<<<< .courant 
    542494      ! 
    543495      REWIND( numnam_ice_ref )              ! Namelist namiceini in reference namelist : Ice initial state 
     
    549501902   IF( ios /= 0 ) CALL ctl_nam ( ios , 'namiceini in configuration namelist', lwp ) 
    550502      WRITE ( numoni, namiceini ) 
    551       ! 
    552       IF(lwp) THEN                        ! control print 
     503======= 
     504 
     505      ! Define the initial parameters 
     506      ! ------------------------- 
     507 
     508      ! Read Namelist namiceini  
     509      REWIND ( numnam_ice ) 
     510      READ   ( numnam_ice , namiceini ) 
     511>>>>>>> .fusion-droit.r4160 
     512      IF(lwp) THEN 
    553513         WRITE(numout,*) 
    554514         WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 
     
    556516         WRITE(numout,*) '   threshold water temp. for initial sea-ice    ttest      = ', ttest 
    557517         WRITE(numout,*) '   initial snow thickness in the north          hninn      = ', hninn 
    558          WRITE(numout,*) '   initial undef ice thickness in the north     hginn_u    = ', hginn_u 
    559          WRITE(numout,*) '   initial undef ice concentr. in the north     aginn_u    = ', aginn_u           
    560          WRITE(numout,*) '   initial  def  ice thickness in the north     hginn_d    = ', hginn_d 
    561          WRITE(numout,*) '   initial  def  ice concentr. in the north     aginn_d    = ', aginn_d           
    562518         WRITE(numout,*) '   initial snow thickness in the south          hnins      = ', hnins  
    563          WRITE(numout,*) '   initial undef ice thickness in the north     hgins_u    = ', hgins_u 
    564          WRITE(numout,*) '   initial undef ice concentr. in the north     agins_u    = ', agins_u           
    565          WRITE(numout,*) '   initial  def  ice thickness in the north     hgins_d    = ', hgins_d 
    566          WRITE(numout,*) '   initial  def  ice concentr. in the north     agins_d    = ', agins_d           
    567          WRITE(numout,*) '   initial  ice salinity       in the north     sinn       = ', sinn 
    568          WRITE(numout,*) '   initial  ice salinity       in the south     sins       = ', sins 
     519         WRITE(numout,*) '   initial ice thickness  in the north          hginn      = ', hginn 
     520         WRITE(numout,*) '   initial ice thickness  in the south          hgins      = ', hgins 
     521         WRITE(numout,*) '   initial ice concentr.  in the north          aginn      = ', aginn 
     522         WRITE(numout,*) '   initial ice concentr.  in the north          agins      = ', agins 
     523         WRITE(numout,*) '   initial  ice salinity  in the north          sinn       = ', sinn 
     524         WRITE(numout,*) '   initial  ice salinity  in the south          sins       = ', sins 
    569525      ENDIF 
    570       ! 
     526 
    571527   END SUBROUTINE lim_istate_init 
    572528 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4147 r4161  
    55   !!====================================================================== 
    66   !! History :  LIM  ! 2006-02  (M. Vancoppenolle) Original code  
    7    !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & fsalt_rpo 
     7   !!            3.2  ! 2009-07  (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in smsw & sfx_mec 
    88   !!            4.0  ! 2011-02  (G. Madec) dynamical allocation 
    99   !!---------------------------------------------------------------------- 
     
    1212   !!   'key_lim3'                                      LIM-3 sea-ice model 
    1313   !!---------------------------------------------------------------------- 
    14    USE par_oce        ! ocean parameters 
    15    USE dom_oce        ! ocean domain 
    16    USE phycst         ! physical constants (ocean directory)  
    17    USE sbc_oce        ! surface boundary condition: ocean fields 
    18    USE thd_ice        ! LIM thermodynamics 
    19    USE ice            ! LIM variables 
    20    USE par_ice        ! LIM parameters 
    21    USE dom_ice        ! LIM domain 
    22    USE limthd_lac     ! LIM 
    23    USE limvar         ! LIM 
    24    USE limcons        ! LIM 
    25    USE in_out_manager ! I/O manager 
    26    USE lbclnk         ! lateral boundary condition - MPP exchanges 
    27    USE lib_mpp        ! MPP library 
    28    USE wrk_nemo       ! work arrays 
    29    USE prtctl         ! Print control 
    30    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    31    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     14   USE par_oce          ! ocean parameters 
     15   USE dom_oce          ! ocean domain 
     16   USE phycst           ! physical constants (ocean directory)  
     17   USE sbc_oce          ! surface boundary condition: ocean fields 
     18   USE thd_ice          ! LIM thermodynamics 
     19   USE ice              ! LIM variables 
     20   USE par_ice          ! LIM parameters 
     21   USE dom_ice          ! LIM domain 
     22   USE limthd_lac       ! LIM 
     23   USE limvar           ! LIM 
     24   USE limcons          ! LIM 
     25   USE in_out_manager   ! I/O manager 
     26   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     27   USE lib_mpp          ! MPP library 
     28   USE wrk_nemo         ! work arrays 
     29   USE prtctl           ! Print control 
     30  ! Check budget (Rousset) 
     31   USE iom              ! I/O manager 
     32   USE lib_fortran     ! glob_sum 
     33   USE limdiahsb 
     34   USE timing          ! Timing 
    3235 
    3336   IMPLICIT NONE 
     
    6265   REAL(wp), PARAMETER ::   krdgmin = 1.1_wp    ! min ridge thickness multiplier 
    6366   REAL(wp), PARAMETER ::   kraft   = 2.0_wp    ! rafting multipliyer 
     67   REAL(wp), PARAMETER ::   kamax   = 1.0 
    6468 
    6569   REAL(wp) ::   Cp                             !  
     
    141145      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    142146      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
     147      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     148      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     149      ! mass and salt flux (clem) 
     150      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold, zsmvold   ! old ice volume... 
    143151      !!----------------------------------------------------------------------------- 
     152      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
    144153 
    145154      CALL wrk_alloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
     155 
     156      CALL wrk_alloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
    146157 
    147158      IF( numit == nstart  )   CALL lim_itd_me_init   ! Initialization (first time-step only) 
     
    151162         CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ') 
    152163      ENDIF 
     164 
     165      IF( ln_limdyn ) THEN          !   Start ridging and rafting   ! 
     166      ! ------------------------------- 
     167      !- check conservation (C Rousset) 
     168      IF (ln_limdiahsb) THEN 
     169         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     170         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     171         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
     172         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     173      ENDIF 
     174      !- check conservation (C Rousset) 
     175      ! ------------------------------- 
     176 
     177      ! mass and salt flux init (clem) 
     178      zviold(:,:,:) = v_i(:,:,:) 
     179      zvsold(:,:,:) = v_s(:,:,:) 
     180      zsmvold(:,:,:) = smv_i(:,:,:) 
    153181 
    154182      !-----------------------------------------------------------------------------! 
     
    204232            ! to give asum = 1.0 after ridging. 
    205233 
    206             divu_adv(ji,jj) = ( 1._wp - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
     234            divu_adv(ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice  ! asum found in ridgeprep 
    207235 
    208236            IF( divu_adv(ji,jj) < 0._wp )   closing_net(ji,jj) = MAX( closing_net(ji,jj), -divu_adv(ji,jj) ) 
     
    288316         DO jj = 1, jpj 
    289317            DO ji = 1, jpi 
    290                IF (ABS(asum(ji,jj) - 1.0) .LT. epsi11) THEN 
     318               IF (ABS(asum(ji,jj) - kamax ) .LT. epsi11) THEN 
    291319                  closing_net(ji,jj) = 0._wp 
    292320                  opning     (ji,jj) = 0._wp 
    293321               ELSE 
    294322                  iterate_ridging    = 1 
    295                   divu_adv   (ji,jj) = (1._wp - asum(ji,jj)) * r1_rdtice 
     323                  divu_adv   (ji,jj) = ( kamax - asum(ji,jj) ) * r1_rdtice 
    296324                  closing_net(ji,jj) = MAX( 0._wp, -divu_adv(ji,jj) ) 
    297325                  opning     (ji,jj) = MAX( 0._wp,  divu_adv(ji,jj) ) 
     
    330358         DO ji = 1, jpi 
    331359 
    332             IF( ABS( asum(ji,jj) - 1.0 ) > epsi11 )  asum_error = .true. 
     360            IF(ABS(asum(ji,jj) - kamax) > epsi11 ) asum_error = .true. 
    333361 
    334362            dardg1dt(ji,jj) = dardg1dt(ji,jj) * r1_rdtice 
     
    349377      DO jj = 1, jpj 
    350378         DO ji = 1, jpi 
    351             IF( ABS( asum(ji,jj) - 1._wp )  >  epsi11 ) THEN   ! there is a bug 
     379            IF( ABS( asum(ji,jj) - kamax)  >  epsi11 ) THEN   ! there is a bug 
    352380               WRITE(numout,*) ' ' 
    353381               WRITE(numout,*) ' ALERTE : Ridging error: total area = ', asum(ji,jj) 
     
    377405      CALL lim_var_glo2eqv 
    378406      CALL lim_itd_me_zapsmall 
     407 
     408      !-------------------------------- 
     409      ! Update mass/salt fluxes (clem) 
     410      !-------------------------------- 
     411      DO jl = 1, jpl 
     412         DO jj = 1, jpj  
     413            DO ji = 1, jpi 
     414               diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * r1_rdtice 
     415               rdm_ice(ji,jj) = rdm_ice(ji,jj) + ( v_i(ji,jj,jl) - zviold(ji,jj,jl) ) * rhoic  
     416               rdm_snw(ji,jj) = rdm_snw(ji,jj) + ( v_s(ji,jj,jl) - zvsold(ji,jj,jl) ) * rhosn  
     417               sfx_mec(ji,jj) = sfx_mec(ji,jj) - ( smv_i(ji,jj,jl) - zsmvold(ji,jj,jl) ) * rhoic * r1_rdtice  
     418            END DO 
     419         END DO 
     420      END DO 
    379421 
    380422      !----------------- 
     
    425467      ENDIF 
    426468 
     469      ! ------------------------------- 
     470      !- check conservation (C Rousset) 
     471      IF (ln_limdiahsb) THEN 
     472         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     473         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     474  
     475         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
     476         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
     477 
     478         zchk_vmin = glob_min(v_i) 
     479         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     480         zchk_amin = glob_min(a_i) 
     481        
     482         IF(lwp) THEN 
     483            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_me) = ',(zchk_v_i * rday) 
     484            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 
     485            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3) 
     486            IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax 
     487            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin 
     488         ENDIF 
     489      ENDIF 
     490      !- check conservation (C Rousset) 
     491      ! ------------------------------- 
     492 
    427493      !-------------------------! 
    428494      ! Back to initial values 
     
    448514 
    449515      ! heat content has to be corrected before ice volume 
    450       DO jl = 1, jpl 
    451          DO jk = 1, nlay_i 
    452             DO jj = 1, jpj 
    453                DO ji = 1, jpi 
    454                   IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
    455                      ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
    456                      old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
    457                      d_e_i_trp(ji,jj,jk,jl) = 0._wp 
    458                   ENDIF 
    459                END DO 
    460             END DO 
    461          END DO 
    462       END DO 
    463  
    464       DO jl = 1, jpl 
    465          DO jj = 1, jpj 
    466             DO ji = 1, jpi 
    467                IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. & 
    468                     d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
    469                   old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
    470                   d_v_i_trp (ji,jj,jl)   = 0._wp 
    471                   old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
    472                   d_a_i_trp (ji,jj,jl)   = 0._wp 
    473                   old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
    474                   d_v_s_trp (ji,jj,jl)   = 0._wp 
    475                   old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
    476                   d_e_s_trp (ji,jj,1,jl) = 0._wp 
    477                   old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
    478                   d_oa_i_trp(ji,jj,jl)   = 0._wp 
    479                   IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
    480                   d_smv_i_trp(ji,jj,jl)  = 0._wp 
    481                ENDIF 
    482             END DO 
    483          END DO 
    484       END DO 
     516!clem@order 
     517!      DO jl = 1, jpl 
     518!         DO jk = 1, nlay_i 
     519!            DO jj = 1, jpj 
     520!               DO ji = 1, jpi 
     521!                  IF ( ( old_v_i(ji,jj,jl) < epsi06 ) .AND. & 
     522!                     ( d_v_i_trp(ji,jj,jl) > epsi06 ) ) THEN 
     523!                     old_e_i(ji,jj,jk,jl)   = d_e_i_trp(ji,jj,jk,jl) 
     524!                     d_e_i_trp(ji,jj,jk,jl) = 0._wp 
     525!                  ENDIF 
     526!               END DO 
     527!            END DO 
     528!         END DO 
     529!      END DO 
     530! 
     531!      DO jl = 1, jpl 
     532!         DO jj = 1, jpj 
     533!            DO ji = 1, jpi 
     534!               IF(  old_v_i  (ji,jj,jl) < epsi06  .AND. & 
     535!                    d_v_i_trp(ji,jj,jl) > epsi06    ) THEN 
     536!                  old_v_i   (ji,jj,jl)   = d_v_i_trp(ji,jj,jl) 
     537!                  d_v_i_trp (ji,jj,jl)   = 0._wp 
     538!                  old_a_i   (ji,jj,jl)   = d_a_i_trp(ji,jj,jl) 
     539!                  d_a_i_trp (ji,jj,jl)   = 0._wp 
     540!                  old_v_s   (ji,jj,jl)   = d_v_s_trp(ji,jj,jl) 
     541!                  d_v_s_trp (ji,jj,jl)   = 0._wp 
     542!                  old_e_s   (ji,jj,1,jl) = d_e_s_trp(ji,jj,1,jl) 
     543!                  d_e_s_trp (ji,jj,1,jl) = 0._wp 
     544!                  old_oa_i  (ji,jj,jl)   = d_oa_i_trp(ji,jj,jl) 
     545!                  d_oa_i_trp(ji,jj,jl)   = 0._wp 
     546!                  IF(  num_sal == 2  )   old_smv_i(ji,jj,jl) = d_smv_i_trp(ji,jj,jl) 
     547!                  d_smv_i_trp(ji,jj,jl)  = 0._wp 
     548!               ENDIF 
     549!            END DO 
     550!         END DO 
     551!      END DO 
     552!clem@order 
     553      ENDIF  ! ln_limdyn=.true. 
    485554      ! 
    486555      CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross, msnow_mlt, esnow_mlt, vt_i_init, vt_i_final ) 
    487556      ! 
     557      CALL wrk_dealloc( jpi, jpj, jpl, zviold, zvsold, zsmvold )   ! clem 
     558      ! 
     559      IF( nn_timing == 1 )  CALL timing_stop('limitd_me') 
    488560   END SUBROUTINE lim_itd_me 
    489561 
     
    10861158            afrft(ji,jj) = arft1(ji,jj) / aicen_init(ji,jj,jl1) !rafting 
    10871159 
    1088             IF (afrac(ji,jj) > 1.0 + epsi11) THEN  !riging 
     1160            IF (afrac(ji,jj) > kamax + epsi11) THEN  !riging 
    10891161               large_afrac = .true. 
    1090             ELSEIF (afrac(ji,jj) > 1.0) THEN  ! roundoff error 
    1091                afrac(ji,jj) = 1.0 
     1162            ELSEIF (afrac(ji,jj) > kamax) THEN  ! roundoff error 
     1163               afrac(ji,jj) = kamax 
    10921164            ENDIF 
    1093             IF (afrft(ji,jj) > 1.0 + epsi11) THEN !rafting 
     1165            IF (afrft(ji,jj) > kamax + epsi11) THEN !rafting 
    10941166               large_afrft = .true. 
    1095             ELSEIF (afrft(ji,jj) > 1.0) THEN  ! roundoff error 
    1096                afrft(ji,jj) = 1.0 
     1167            ELSEIF (afrft(ji,jj) > kamax) THEN  ! roundoff error 
     1168               afrft(ji,jj) = kamax 
    10971169            ENDIF 
    10981170 
     
    11371209             
    11381210            !                                                             ! excess of salt is flushed into the ocean 
    1139             sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
    1140  
    1141             rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic / rau0   ! increase in ice volume du to seawater frozen in voids 
    1142              
     1211            !sfx_mec(ji,jj) = sfx_mec(ji,jj) + ( zsrdg2 - srdg2(ji,jj) ) * rhoic * r1_rdtice 
     1212 
     1213            !rdm_ice(ji,jj) = rdm_ice(ji,jj) + vsw(ji,jj) * rhoic    ! gurvan: increase in ice volume du to seawater frozen in voids              
     1214 
    11431215            !------------------------------------             
    11441216            ! 3.6 Increment ridging diagnostics 
     
    11501222            dardg1dt   (ji,jj) = dardg1dt(ji,jj) + ardg1(ji,jj) + arft1(ji,jj) 
    11511223            dardg2dt   (ji,jj) = dardg2dt(ji,jj) + ardg2(ji,jj) + arft2(ji,jj) 
    1152             diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
     1224            !clem diag_dyn_gr(ji,jj) = diag_dyn_gr(ji,jj) + ( vrdg2(ji,jj) + virft(ji,jj) ) * r1_rdtice 
    11531225            opening    (ji,jj) = opening (ji,jj) + opning(ji,jj) * rdt_ice 
    11541226 
     
    12171289 
    12181290               ! Mutliply by ice volume, and divide by number of layers to get heat content in 10^9 J 
    1219                ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / nlay_i 
     1291               ersw (ji,jj,jk)  = ersw(ji,jj,jk) * area(ji,jj) * vsw(ji,jj) / REAL( nlay_i ) 
    12201292 
    12211293               erdg2(ji,jj,jk)  = erdg1(ji,jj,jk) + ersw(ji,jj,jk) 
     
    12401312               ji = indxi(ij) 
    12411313               jj = indxj(ij) 
    1242                IF( afrac(ji,jj) > 1.0 + epsi11 ) THEN  
     1314               IF( afrac(ji,jj) > kamax + epsi11 ) THEN  
    12431315                  WRITE(numout,*) '' 
    12441316                  WRITE(numout,*) ' ardg > a_i' 
     
    12521324               ji = indxi(ij) 
    12531325               jj = indxj(ij) 
    1254                IF( afrft(ji,jj) > 1.0 + epsi11 ) THEN  
     1326               IF( afrft(ji,jj) > kamax + epsi11 ) THEN  
    12551327                  WRITE(numout,*) '' 
    12561328                  WRITE(numout,*) ' arft > a_i' 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r3764 r4161  
    1919   !!   lim_itd_shiftice : 
    2020   !!---------------------------------------------------------------------- 
    21    USE par_oce        ! ocean parameters 
    22    USE dom_oce        ! ocean domain 
    23    USE phycst         ! physical constants (ocean directory)  
    24    USE ice            ! LIM-3 variables 
    25    USE par_ice        ! LIM-3 parameters 
    26    USE dom_ice        ! LIM-3 domain 
    27    USE thd_ice        ! LIM-3 thermodynamic variables 
    28    USE limthd_lac     ! LIM-3 lateral accretion 
    29    USE limvar         ! LIM-3 variables 
    30    USE limcons        ! LIM-3 conservation 
    31    USE prtctl         ! Print control 
    32    USE in_out_manager ! I/O manager 
    33    USE lib_mpp        ! MPP library 
    34    USE wrk_nemo       ! work arrays 
    35    USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    36    USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 
     21   USE dom_ice          ! LIM-3 domain 
     22   USE par_oce          ! ocean parameters 
     23   USE dom_oce          ! ocean domain 
     24   USE phycst           ! physical constants (ocean directory)  
     25   USE thd_ice          ! LIM-3 thermodynamic variables 
     26   USE ice              ! LIM-3 variables 
     27   USE par_ice          ! LIM-3 parameters 
     28   USE limthd_lac       ! LIM-3 lateral accretion 
     29   USE limvar           ! LIM-3 variables 
     30   USE limcons          ! LIM-3 conservation 
     31   USE prtctl           ! Print control 
     32   USE in_out_manager   ! I/O manager 
     33   USE lib_mpp          ! MPP library 
     34   USE wrk_nemo         ! work arrays 
     35   USE lib_fortran      ! to use key_nosignedzero 
     36   USE timing          ! Timing 
    3737 
    3838   IMPLICIT NONE 
     
    4545   PUBLIC   lim_itd_shiftice 
    4646 
    47    REAL(wp) ::   epsi20 = 1e-20_wp   ! constant values 
    48    REAL(wp) ::   epsi13 = 1e-13_wp   ! 
    49    REAL(wp) ::   epsi10 = 1e-10_wp   ! 
     47   REAL(wp) ::   epsi20 = 1.e-20_wp   ! constant values 
     48   REAL(wp) ::   epsi13 = 1.e-13_wp   ! 
     49   REAL(wp) ::   epsi10 = 1.e-10_wp   ! 
    5050 
    5151   !!---------------------------------------------------------------------- 
    52    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
     52   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    5353   !! $Id$ 
    5454   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    6767      ! 
    6868      INTEGER ::   jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index          
    69  
    70       !!------------------------------------------------------------------ 
     69      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     70      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     71      !!------------------------------------------------------------------ 
     72      IF( nn_timing == 1 )  CALL timing_start('limitd_th') 
     73 
     74      ! ------------------------------- 
     75      !- check conservation (C Rousset) 
     76      IF (ln_limdiahsb) THEN 
     77         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     78         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     79         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
     80         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     81       ENDIF 
     82      !- check conservation (C Rousset) 
     83      ! ------------------------------- 
    7184 
    7285      IF( kt == nit000 .AND. lwp ) THEN 
     
    107120      d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)  
    108121      d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:) 
    109  
     122      !?? d_oa_i_thd(:,:,:)  = oa_i (:,:,:) - old_oa_i (:,:,:) 
    110123      d_smv_i_thd(:,:,:) = 0._wp 
    111       IF( num_sal == 2  )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     124      IF( num_sal == 2 )   d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     125 
     126      ! diag only (clem) 
     127      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    112128 
    113129      IF(ln_ctl) THEN   ! Control print 
     
    142158         END DO 
    143159      ENDIF 
    144  
     160      ! 
     161      ! ------------------------------- 
     162      !- check conservation (C Rousset) 
     163      IF( ln_limdiahsb ) THEN 
     164         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     165         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     166  
     167         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
     168         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
     169 
     170         zchk_vmin = glob_min(v_i) 
     171         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     172         zchk_amin = glob_min(a_i) 
     173 
     174         IF(lwp) THEN 
     175            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limitd_th) = ',(zchk_v_i * rday) 
     176            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 
     177            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(zchk_vmin * 1.e-3) 
     178            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_th) = ',zchk_amax 
     179            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_th) = ',zchk_amin 
     180         ENDIF 
     181       ENDIF 
     182      !- check conservation (C Rousset) 
     183      ! ------------------------------- 
     184      ! 
    145185      !- Recover Old values 
    146       a_i(:,:,:)   = old_a_i(:,:,:) 
    147       v_s(:,:,:)   = old_v_s(:,:,:) 
    148       v_i(:,:,:)   = old_v_i(:,:,:) 
    149       e_s(:,:,:,:) = old_e_s(:,:,:,:) 
    150       e_i(:,:,:,:) = old_e_i(:,:,:,:) 
    151       ! 
     186      a_i(:,:,:)   = old_a_i (:,:,:) 
     187      v_s(:,:,:)   = old_v_s (:,:,:) 
     188      v_i(:,:,:)   = old_v_i (:,:,:) 
     189      e_s(:,:,:,:) = old_e_s (:,:,:,:) 
     190      e_i(:,:,:,:) = old_e_i (:,:,:,:) 
     191      !?? oa_i(:,:,:)  = old_oa_i(:,:,:) 
    152192      IF( num_sal == 2 )   smv_i(:,:,:) = old_smv_i(:,:,:) 
    153193      ! 
     194      IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
    154195   END SUBROUTINE lim_itd_th 
    155196   ! 
     
    172213      ! 
    173214      INTEGER  ::   ji, jj, jl     ! dummy loop index 
    174       INTEGER  ::   zji, zjj, nd   ! local integer 
     215      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
     216      INTEGER  ::   nd             ! local integer 
    175217      REAL(wp) ::   zx1, zwk1, zdh0, zetamin, zdamax   ! local scalars 
    176       REAL(wp) ::   zx2, zwk2, zda0, zetamax, zhimin   !   -      - 
     218      REAL(wp) ::   zx2, zwk2, zda0, zetamax           !   -      - 
    177219      REAL(wp) ::   zx3,             zareamin, zindb   !   -      - 
    178220      CHARACTER (len = 15) :: fieldid 
     
    210252      CALL wrk_alloc( jpi,jpj, zhb0,zhb1,vt_i_init,vt_i_final,vt_s_init,vt_s_final,et_i_init,et_i_final,et_s_init,et_s_final ) 
    211253 
    212       zhimin   = 0.1      !minimum ice thickness tolerated by the model 
    213254      zareamin = epsi10   !minimum area in thickness categories tolerated by the conceptors of the model 
    214255 
     
    240281         DO jj = 1, jpj 
    241282            DO ji = 1, jpi 
    242                zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)))     !0 if no ice and 1 if yes 
     283               zindb             = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl)+epsi10))     !0 if no ice and 1 if yes 
    243284               ht_i(ji,jj,jl)    = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),epsi10) * zindb 
    244                zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes 
     285               zindb             = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
    245286               zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),epsi10) * zindb 
    246287               IF( a_i(ji,jj,jl) > 1e-6 )   zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)  
     
    285326      DO jl = klbnd, kubnd - 1 
    286327         DO ji = 1, nbrem 
    287             zji = nind_i(ji) 
    288             zjj = nind_j(ji) 
     328            ii = nind_i(ji) 
     329            ij = nind_j(ji) 
    289330            ! 
    290             IF ( ( zht_i_o(zji,zjj,jl)  .GT.epsi10 ) .AND. &  
    291                ( zht_i_o(zji,zjj,jl+1).GT.epsi10 ) ) THEN 
     331            IF ( ( zht_i_o(ii,ij,jl)  .GT.epsi10 ) .AND. &  
     332               ( zht_i_o(ii,ij,jl+1).GT.epsi10 ) ) THEN 
    292333               !interpolate between adjacent category growth rates 
    293                zslope = ( zdhice(zji,zjj,jl+1)     - zdhice(zji,zjj,jl) ) / & 
    294                   ( zht_i_o   (zji,zjj,jl+1) - zht_i_o   (zji,zjj,jl) ) 
    295                zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + & 
    296                   zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) ) 
    297             ELSEIF (zht_i_o(zji,zjj,jl).gt.epsi10) THEN 
    298                zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) 
    299             ELSEIF (zht_i_o(zji,zjj,jl+1).gt.epsi10) THEN 
    300                zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1) 
     334               zslope = ( zdhice(ii,ij,jl+1)     - zdhice(ii,ij,jl) ) / & 
     335                  ( zht_i_o   (ii,ij,jl+1) - zht_i_o   (ii,ij,jl) ) 
     336               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) + & 
     337                  zslope * ( hi_max(jl) - zht_i_o(ii,ij,jl) ) 
     338            ELSEIF (zht_i_o(ii,ij,jl).gt.epsi10) THEN 
     339               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl) 
     340            ELSEIF (zht_i_o(ii,ij,jl+1).gt.epsi10) THEN 
     341               zhbnew(ii,ij,jl) = hi_max(jl) + zdhice(ii,ij,jl+1) 
    301342            ELSE 
    302                zhbnew(zji,zjj,jl) = hi_max(jl) 
     343               zhbnew(ii,ij,jl) = hi_max(jl) 
    303344            ENDIF 
    304345         END DO 
     
    307348         DO ji = 1, nbrem 
    308349            ! jl, ji 
    309             zji = nind_i(ji) 
    310             zjj = nind_j(ji) 
     350            ii = nind_i(ji) 
     351            ij = nind_j(ji) 
    311352            ! jl, ji 
    312             IF ( ( a_i(zji,zjj,jl) .GT.epsi10) .AND. &  
    313                ( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) & 
     353            IF ( ( a_i(ii,ij,jl) .GT.epsi10) .AND. &  
     354               ( ht_i(ii,ij,jl).GE. zhbnew(ii,ij,jl) ) & 
    314355               ) THEN 
    315                zremap_flag(zji,zjj) = 0 
    316             ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. epsi10 ) .AND. & 
    317                ( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) & 
     356               zremap_flag(ii,ij) = 0 
     357            ELSEIF ( ( a_i(ii,ij,jl+1) .GT. epsi10 ) .AND. & 
     358               ( ht_i(ii,ij,jl+1).LE. zhbnew(ii,ij,jl) ) & 
    318359               ) THEN 
    319                zremap_flag(zji,zjj) = 0 
     360               zremap_flag(ii,ij) = 0 
    320361            ENDIF 
    321362 
    322363            !- 4.3 Check that each zhbnew does not exceed maximal values hi_max   
    323364            ! jl, ji 
    324             IF (zhbnew(zji,zjj,jl).gt.hi_max(jl+1)) THEN 
    325                zremap_flag(zji,zjj) = 0 
     365            IF (zhbnew(ii,ij,jl).gt.hi_max(jl+1)) THEN 
     366               zremap_flag(ii,ij) = 0 
    326367            ENDIF 
    327368            ! jl, ji 
    328             IF (zhbnew(zji,zjj,jl).lt.hi_max(jl-1)) THEN 
    329                zremap_flag(zji,zjj) = 0 
     369            IF (zhbnew(ii,ij,jl).lt.hi_max(jl-1)) THEN 
     370               zremap_flag(ii,ij) = 0 
    330371            ENDIF 
    331372            ! jl, ji 
     
    379420      !- 7.2 Area lost due to melting of thin ice (first category,  klbnd) 
    380421      DO ji = 1, nbrem 
    381          zji = nind_i(ji)  
    382          zjj = nind_j(ji)  
     422         ii = nind_i(ji)  
     423         ij = nind_j(ji)  
    383424 
    384425         !ji 
    385          IF (a_i(zji,zjj,klbnd) .gt. epsi10) THEN 
    386             zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category 
     426         IF (a_i(ii,ij,klbnd) .gt. epsi10) THEN 
     427            zdh0 = zdhice(ii,ij,klbnd) !decrease of ice thickness in the lower category 
    387428            ! ji, a_i > epsi10 
    388429            IF (zdh0 .lt. 0.0) THEN !remove area from category 1 
     
    391432 
    392433               !Integrate g(1) from 0 to dh0 to estimate area melted 
    393                zetamax = MIN(zdh0,hR(zji,zjj,klbnd)) - hL(zji,zjj,klbnd) 
     434               zetamax = MIN(zdh0,hR(ii,ij,klbnd)) - hL(ii,ij,klbnd) 
    394435               IF (zetamax.gt.0.0) THEN 
    395436                  zx1  = zetamax 
    396437                  zx2  = 0.5 * zetamax*zetamax  
    397                   zda0 = g1(zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed 
     438                  zda0 = g1(ii,ij,klbnd) * zx2 + g0(ii,ij,klbnd) * zx1 !ice area removed 
    398439                  ! Constrain new thickness <= ht_i 
    399                   zdamax = a_i(zji,zjj,klbnd) * &  
    400                      (1.0 - ht_i(zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0 
     440                  zdamax = a_i(ii,ij,klbnd) * &  
     441                     (1.0 - ht_i(ii,ij,klbnd)/zht_i_o(ii,ij,klbnd)) ! zdamax > 0 
    401442                  !ice area lost due to melting of thin ice 
    402443                  zda0   = MIN(zda0, zdamax) 
    403444 
    404445                  ! Remove area, conserving volume 
    405                   ht_i(zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &  
    406                      * a_i(zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 ) 
    407                   a_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd) - zda0 
    408                   v_i(zji,zjj,klbnd)  = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd) 
     446                  ht_i(ii,ij,klbnd) = ht_i(ii,ij,klbnd) &  
     447                     * a_i(ii,ij,klbnd) / ( a_i(ii,ij,klbnd) - zda0 ) 
     448                  a_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd) - zda0 
     449                  v_i(ii,ij,klbnd)  = a_i(ii,ij,klbnd)*ht_i(ii,ij,klbnd) ! clem@useless ? 
    409450               ENDIF     ! zetamax > 0 
    410451               ! ji, a_i > epsi10 
     
    412453            ELSE ! if ice accretion 
    413454               ! ji, a_i > epsi10; zdh0 > 0 
    414                IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
     455               IF ( ntyp .EQ. 1 ) zhbnew(ii,ij,klbnd-1) = MIN(zdh0,hi_max(klbnd))  
    415456               ! zhbnew was 0, and is shifted to the right to account for thin ice 
    416457               ! growth in openwater (F0 = f1) 
    417                IF ( ntyp .NE. 1 ) zhbnew(zji,zjj,0) = 0  
     458               IF ( ntyp .NE. 1 ) zhbnew(ii,ij,0) = 0  
    418459               ! in other types there is 
    419460               ! no open water growth (F0 = 0) 
     
    446487 
    447488         DO ji = 1, nbrem 
    448             zji = nind_i(ji) 
    449             zjj = nind_j(ji) 
    450  
    451             IF (zhbnew(zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 
     489            ii = nind_i(ji) 
     490            ij = nind_j(ji) 
     491 
     492            IF (zhbnew(ii,ij,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1 
    452493 
    453494               ! left and right integration limits in eta space 
    454                zvetamin(ji) = MAX(hi_max(jl), hL(zji,zjj,jl)) - hL(zji,zjj,jl) 
    455                zvetamax(ji) = MIN(zhbnew(zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl) 
    456                zdonor(zji,zjj,jl) = jl 
     495               zvetamin(ji) = MAX(hi_max(jl), hL(ii,ij,jl)) - hL(ii,ij,jl) 
     496               zvetamax(ji) = MIN(zhbnew(ii,ij,jl), hR(ii,ij,jl)) - hL(ii,ij,jl) 
     497               zdonor(ii,ij,jl) = jl 
    457498 
    458499            ELSE  ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl 
     
    460501               ! left and right integration limits in eta space 
    461502               zvetamin(ji) = 0.0 
    462                zvetamax(ji) = MIN(hi_max(jl), hR(zji,zjj,jl+1)) - hL(zji,zjj,jl+1) 
    463                zdonor(zji,zjj,jl) = jl + 1 
     503               zvetamax(ji) = MIN(hi_max(jl), hR(ii,ij,jl+1)) - hL(ii,ij,jl+1) 
     504               zdonor(ii,ij,jl) = jl + 1 
    464505 
    465506            ENDIF  ! zhbnew(jl) > hi_max(jl) 
     
    475516            zwk2 = zwk2 * zetamax 
    476517            zx3  = 1.0/3.0 * (zwk2 - zwk1) 
    477             nd   = zdonor(zji,zjj,jl) 
    478             zdaice(zji,zjj,jl) = g1(zji,zjj,nd)*zx2 + g0(zji,zjj,nd)*zx1 
    479             zdvice(zji,zjj,jl) = g1(zji,zjj,nd)*zx3 + g0(zji,zjj,nd)*zx2 + & 
    480                zdaice(zji,zjj,jl)*hL(zji,zjj,nd) 
     518            nd   = zdonor(ii,ij,jl) 
     519            zdaice(ii,ij,jl) = g1(ii,ij,nd)*zx2 + g0(ii,ij,nd)*zx1 
     520            zdvice(ii,ij,jl) = g1(ii,ij,nd)*zx3 + g0(ii,ij,nd)*zx2 + & 
     521               zdaice(ii,ij,jl)*hL(ii,ij,nd) 
    481522 
    482523         END DO ! ji 
     
    493534 
    494535      DO ji = 1, nbrem 
    495          zji = nind_i(ji) 
    496          zjj = nind_j(ji) 
    497          IF ( ( zhimin .GT. 0.0 ) .AND. &  
    498             ( ( a_i(zji,zjj,1) .GT. epsi10 ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) & 
    499             ) THEN 
    500             a_i(zji,zjj,1)  = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin  
    501             ht_i(zji,zjj,1) = zhimin 
    502             v_i(zji,zjj,1)  = a_i(zji,zjj,1)*ht_i(zji,zjj,1) 
     536         ii = nind_i(ji) 
     537         ij = nind_j(ji) 
     538         IF ( ( a_i(ii,ij,1) > epsi10 ) .AND. ( ht_i(ii,ij,1) < hiclim ) ) THEN 
     539            a_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) / hiclim  
     540            ht_i(ii,ij,1) = hiclim 
     541            v_i(ii,ij,1)  = a_i(ii,ij,1) * ht_i(ii,ij,1) !clem@useless 
    503542         ENDIF 
    504543      END DO !ji 
     
    625664 
    626665      INTEGER ::   ji, jj, jl, jl2, jl1, jk   ! dummy loop indices 
    627       INTEGER ::   zji, zjj          ! indices when changing from 2D-1D is done 
     666      INTEGER ::   ii, ij          ! indices when changing from 2D-1D is done 
    628667 
    629668      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zaTsfn 
     
    759798 
    760799         DO ji = 1, nbrem  
    761             zji = nind_i(ji) 
    762             zjj = nind_j(ji) 
    763  
    764             jl1 = zdonor(zji,zjj,jl) 
    765             zindb             = MAX( 0.0 , SIGN( 1.0 , v_i(zji,zjj,jl1) - epsi10 ) ) 
    766             zworka(zji,zjj)   = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),epsi10) * zindb 
     800            ii = nind_i(ji) 
     801            ij = nind_j(ji) 
     802 
     803            jl1 = zdonor(ii,ij,jl) 
     804            zindb             = MAX( 0.0 , SIGN( 1.0 , v_i(ii,ij,jl1) - epsi10 ) ) 
     805            zworka(ii,ij)   = zdvice(ii,ij,jl) / MAX(v_i(ii,ij,jl1),epsi10) * zindb 
    767806            IF( jl1 == jl) THEN   ;   jl2 = jl1+1 
    768807            ELSE                    ;   jl2 = jl  
     
    773812            !-------------- 
    774813 
    775             a_i(zji,zjj,jl1) = a_i(zji,zjj,jl1) - zdaice(zji,zjj,jl) 
    776             a_i(zji,zjj,jl2) = a_i(zji,zjj,jl2) + zdaice(zji,zjj,jl) 
     814            a_i(ii,ij,jl1) = a_i(ii,ij,jl1) - zdaice(ii,ij,jl) 
     815            a_i(ii,ij,jl2) = a_i(ii,ij,jl2) + zdaice(ii,ij,jl) 
    777816 
    778817            !-------------- 
     
    780819            !-------------- 
    781820 
    782             v_i(zji,zjj,jl1) = v_i(zji,zjj,jl1) - zdvice(zji,zjj,jl)  
    783             v_i(zji,zjj,jl2) = v_i(zji,zjj,jl2) + zdvice(zji,zjj,jl) 
     821            v_i(ii,ij,jl1) = v_i(ii,ij,jl1) - zdvice(ii,ij,jl)  
     822            v_i(ii,ij,jl2) = v_i(ii,ij,jl2) + zdvice(ii,ij,jl) 
    784823 
    785824            !-------------- 
     
    787826            !-------------- 
    788827 
    789             zdvsnow          = v_s(zji,zjj,jl1) * zworka(zji,zjj) 
    790             v_s(zji,zjj,jl1) = v_s(zji,zjj,jl1) - zdvsnow 
    791             v_s(zji,zjj,jl2) = v_s(zji,zjj,jl2) + zdvsnow  
     828            zdvsnow          = v_s(ii,ij,jl1) * zworka(ii,ij) 
     829            v_s(ii,ij,jl1) = v_s(ii,ij,jl1) - zdvsnow 
     830            v_s(ii,ij,jl2) = v_s(ii,ij,jl2) + zdvsnow  
    792831 
    793832            !-------------------- 
     
    795834            !-------------------- 
    796835 
    797             zdesnow              = e_s(zji,zjj,1,jl1) * zworka(zji,zjj) 
    798             e_s(zji,zjj,1,jl1)   = e_s(zji,zjj,1,jl1) - zdesnow 
    799             e_s(zji,zjj,1,jl2)   = e_s(zji,zjj,1,jl2) + zdesnow 
     836            zdesnow              = e_s(ii,ij,1,jl1) * zworka(ii,ij) 
     837            e_s(ii,ij,1,jl1)   = e_s(ii,ij,1,jl1) - zdesnow 
     838            e_s(ii,ij,1,jl2)   = e_s(ii,ij,1,jl2) + zdesnow 
    800839 
    801840            !-------------- 
     
    803842            !-------------- 
    804843 
    805             zdo_aice             = oa_i(zji,zjj,jl1) * zdaice(zji,zjj,jl) 
    806             oa_i(zji,zjj,jl1)    = oa_i(zji,zjj,jl1) - zdo_aice 
    807             oa_i(zji,zjj,jl2)    = oa_i(zji,zjj,jl2) + zdo_aice 
     844            zdo_aice             = oa_i(ii,ij,jl1) * zdaice(ii,ij,jl) 
     845            oa_i(ii,ij,jl1)    = oa_i(ii,ij,jl1) - zdo_aice 
     846            oa_i(ii,ij,jl2)    = oa_i(ii,ij,jl2) + zdo_aice 
    808847 
    809848            !-------------- 
     
    811850            !-------------- 
    812851 
    813             zdsm_vice            = smv_i(zji,zjj,jl1) * zworka(zji,zjj) 
    814             smv_i(zji,zjj,jl1)   = smv_i(zji,zjj,jl1) - zdsm_vice 
    815             smv_i(zji,zjj,jl2)   = smv_i(zji,zjj,jl2) + zdsm_vice 
     852            zdsm_vice            = smv_i(ii,ij,jl1) * zworka(ii,ij) 
     853            smv_i(ii,ij,jl1)   = smv_i(ii,ij,jl1) - zdsm_vice 
     854            smv_i(ii,ij,jl2)   = smv_i(ii,ij,jl2) + zdsm_vice 
    816855 
    817856            !--------------------- 
     
    819858            !--------------------- 
    820859 
    821             zdaTsf               = t_su(zji,zjj,jl1) * zdaice(zji,zjj,jl) 
    822             zaTsfn(zji,zjj,jl1)  = zaTsfn(zji,zjj,jl1) - zdaTsf 
    823             zaTsfn(zji,zjj,jl2)  = zaTsfn(zji,zjj,jl2) + zdaTsf  
     860            zdaTsf               = t_su(ii,ij,jl1) * zdaice(ii,ij,jl) 
     861            zaTsfn(ii,ij,jl1)  = zaTsfn(ii,ij,jl1) - zdaTsf 
     862            zaTsfn(ii,ij,jl2)  = zaTsfn(ii,ij,jl2) + zdaTsf  
    824863 
    825864         END DO                 ! ji 
     
    832871!CDIR NODEP 
    833872            DO ji = 1, nbrem 
    834                zji = nind_i(ji) 
    835                zjj = nind_j(ji) 
    836  
    837                jl1 = zdonor(zji,zjj,jl) 
     873               ii = nind_i(ji) 
     874               ij = nind_j(ji) 
     875 
     876               jl1 = zdonor(ii,ij,jl) 
    838877               IF (jl1 .EQ. jl) THEN 
    839878                  jl2 = jl+1 
     
    842881               ENDIF 
    843882 
    844                zdeice = e_i(zji,zjj,jk,jl1) * zworka(zji,zjj) 
    845                e_i(zji,zjj,jk,jl1) =  e_i(zji,zjj,jk,jl1) - zdeice 
    846                e_i(zji,zjj,jk,jl2) =  e_i(zji,zjj,jk,jl2) + zdeice  
     883               zdeice = e_i(ii,ij,jk,jl1) * zworka(ii,ij) 
     884               e_i(ii,ij,jk,jl1) =  e_i(ii,ij,jk,jl1) - zdeice 
     885               e_i(ii,ij,jk,jl2) =  e_i(ii,ij,jk,jl2) + zdeice  
    847886            END DO              ! ji 
    848887         END DO                 ! jk 
     
    860899                  ht_i(ji,jj,jl)  =  v_i   (ji,jj,jl) / a_i(ji,jj,jl)  
    861900                  t_su(ji,jj,jl)  =  zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)  
    862                   zindsn          =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes 
     901                  zindsn          =  1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl)+epsi10)) !0 if no ice and 1 if yes 
    863902               ELSE 
    864903                  ht_i(ji,jj,jl)  = 0._wp 
     
    9671006                  zshiftflag        = 1 
    9681007                  zdonor(ji,jj,jl)  = jl  
    969                   zdaice(ji,jj,jl)  = a_i(ji,jj,jl) 
    970                   zdvice(ji,jj,jl)  = v_i(ji,jj,jl) 
     1008                  ! begin TECLIM change 
     1009                  !zdaice(ji,jj,jl)  = a_i(ji,jj,jl) 
     1010                  !zdvice(ji,jj,jl)  = v_i(ji,jj,jl) 
     1011                  zdaice(ji,jj,jl)  = a_i(ji,jj,jl)/2 
     1012                  zdvice(ji,jj,jl)  = v_i(ji,jj,jl)-zdaice(ji,jj,jl)*(hi_max(jl)+hi_max(jl-1))/2 
     1013                  ! end TECLIM change  
    9711014               ENDIF 
    9721015            END DO                 ! ji 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limmsh.F90

    r3625 r4161  
    2626 
    2727   !!---------------------------------------------------------------------- 
    28    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     28   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    2929   !! $Id$ 
    3030   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r3791 r4161  
    4141   USE agrif_lim2_interp 
    4242#endif 
     43#if defined key_bdy 
     44   USE bdyice_lim 
     45#endif 
    4346 
    4447   IMPLICIT NONE 
     
    5356#  include "vectopt_loop_substitute.h90" 
    5457   !!---------------------------------------------------------------------- 
    55    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     58   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    5659   !! $Id$ 
    5760   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    413416 
    414417               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zzdst*zzdst ) * usecc2 )   
    415                deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 
    416 !!gm faster to replace the line above with simply: 
    417 !!                deltat(ji,jj) = MAX( delta, creepl ) 
    418 !!gm end  
    419  
     418               ! MV rewriting 
     419               ! deltat(ji,jj) = MAX( SQRT(zdd(ji,jj)**2 + (zdt(ji,jj)**2 + zzdst**2)*usecc2), creepl ) 
     420               !!gm faster to replace the line above with simply: 
     421               !!                deltat(ji,jj) = MAX( delta, creepl ) 
     422               !!gm end   
     423               deltat(ji,jj) = delta + creepl 
     424               ! END MV 
    420425               !-Calculate stress tensor components zs1 and zs2  
    421426               !-at centre of grid cells (see section 3.5 of CICE user's guide). 
     
    472477 
    473478         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
     479 
     480!#if defined key_bdy 
     481!         ! clem: change zs1, zs2, zs12 at the boundary for each iteration 
     482!         CALL bdy_ice_lim_dyn( 2, zs1, zs2, zs12 ) 
     483!         CALL lbc_lnk( zs1 (:,:), 'T', 1. ) 
     484!         CALL lbc_lnk( zs2 (:,:), 'T', 1. ) 
     485!         CALL lbc_lnk( zs12(:,:), 'F', 1. ) 
     486!#endif          
    474487 
    475488         ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
     
    520533 
    521534            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    522 #if defined key_agrif 
     535#if defined key_agrif && defined key_lim2 
    523536            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    524537#endif 
     
    548561 
    549562            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    550 #if defined key_agrif 
     563#if defined key_agrif && defined key_lim2 
    551564            CALL agrif_rhg_lim2( jter, nevp, 'V' ) 
    552565#endif 
     
    577590 
    578591            CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    579 #if defined key_agrif 
     592#if defined key_agrif && defined key_lim2 
    580593            CALL agrif_rhg_lim2( jter, nevp , 'V' ) 
    581594#endif 
     
    608621 
    609622            CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    610 #if defined key_agrif 
     623#if defined key_agrif && defined key_lim2 
    611624            CALL agrif_rhg_lim2( jter, nevp, 'U' ) 
    612625#endif 
    613626 
    614627         ENDIF 
     628          
     629!#if defined key_bdy 
     630!         ! clem: change u_ice and v_ice at the boundary for each iteration 
     631!         CALL bdy_ice_lim_dyn( 1 ) 
     632!#endif          
    615633 
    616634         IF(ln_ctl) THEN 
     
    624642         ENDIF 
    625643 
    626          !                                                   ! ==================== ! 
     644         !                                                ! ==================== ! 
    627645      END DO                                              !  end loop over jter  ! 
    628646      !                                                   ! ==================== ! 
    629  
    630647      ! 
    631648      !------------------------------------------------------------------------------! 
    632649      ! 4) Prevent ice velocities when the ice is thin 
    633650      !------------------------------------------------------------------------------! 
    634       ! 
    635       ! If the ice thickness is below 1cm then ice velocity should equal the 
     651      !clem : add hminrhg in the namelist 
     652      ! 
     653      ! If the ice thickness is below hminrhg (5cm) then ice velocity should equal the 
    636654      ! ocean velocity,  
    637655      ! This prevents high velocity when ice is thin 
     
    641659         DO ji = fs_2, fs_jpim1 
    642660            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    643             zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
    644             IF ( zdummy .LE. 5.0e-2 ) THEN 
     661            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     662            zdummy = vt_i(ji,jj) 
     663            IF ( zdummy .LE. hminrhg ) THEN 
    645664               u_ice(ji,jj) = u_oce(ji,jj) 
    646665               v_ice(ji,jj) = v_oce(ji,jj) 
     
    651670      CALL lbc_lnk( u_ice(:,:), 'U', -1. )  
    652671      CALL lbc_lnk( v_ice(:,:), 'V', -1. )  
    653 #if defined key_agrif 
     672#if defined key_agrif && defined key_lim2 
    654673      CALL agrif_rhg_lim2( nevp , nevp, 'U' ) 
    655674      CALL agrif_rhg_lim2( nevp , nevp, 'V' ) 
    656675#endif 
     676#if defined key_bdy 
     677      ! clem: change u_ice and v_ice at the boundary 
     678      CALL bdy_ice_lim_dyn( 1 ) 
     679#endif          
    657680 
    658681      DO jj = k_j1+1, k_jpj-1  
    659682         DO ji = fs_2, fs_jpim1 
    660683            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    661             zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
    662             IF ( zdummy .LE. 5.0e-2 ) THEN 
     684            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     685            zdummy = vt_i(ji,jj) 
     686            IF ( zdummy .LE. hminrhg ) THEN 
    663687               v_ice1(ji,jj)  = 0.5*( (v_ice(ji,jj)+v_ice(ji,jj-1))*e1t(ji+1,jj)   & 
    664688                  &                 +(v_ice(ji+1,jj)+v_ice(ji+1,jj-1))*e1t(ji,jj)) & 
     
    683707            !- zds(:,:): shear on northeast corner of grid cells 
    684708            zindb  = MAX( 0.0, SIGN( 1.0, at_i(ji,jj) - 1.0e-6 ) )  
    685             zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
    686  
    687             IF ( zdummy .LE. 5.0e-2 ) THEN 
     709            !zdummy = zindb * vt_i(ji,jj) / MAX(at_i(ji,jj) , 1.0e-06 ) 
     710            zdummy = vt_i(ji,jj) 
     711            IF ( zdummy .LE. hminrhg ) THEN 
    688712 
    689713               zdd(ji,jj) = ( e2u(ji,jj)*u_ice(ji,jj)                      & 
     
    719743                  &           - e1v( ji  , jj-1 ) * u_ice2(ji  ,jj-1)  ) / area(ji,jj) 
    720744 
    721                deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   &  
    722                   &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 &  
    723                   &                          ) + creepl 
    724  
     745!              deltat(ji,jj) = SQRT(    zdd(ji,jj)*zdd(ji,jj)   &  
     746!                  &                 + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 &  
     747!                  &                          ) + creepl 
     748               ! MV rewriting 
     749               delta = SQRT( zdd(ji,jj)*zdd(ji,jj) + ( zdt(ji,jj)*zdt(ji,jj) + zdst(ji,jj)*zdst(ji,jj) ) * usecc2 )   
     750               deltat(ji,jj) = delta + creepl 
     751               ! END MV 
     752             
    725753            ENDIF ! zdummy 
    726754 
     
    738766            divu_i (ji,jj) = zdd   (ji,jj) 
    739767            delta_i(ji,jj) = deltat(ji,jj) 
     768            ! begin TECLIM change  
     769            zdst(ji,jj)= (  e2u( ji  , jj   ) * v_ice1(ji,jj)           &    
     770               &          - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)         &    
     771               &          + e1v( ji  , jj   ) * u_ice2(ji,jj)           &    
     772               &          - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) / area(ji,jj)  
    740773            shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst(ji,jj) * zdst(ji,jj) ) 
     774            ! end TECLIM change 
    741775         END DO 
    742776      END DO 
    743       CALL lbc_lnk( divu_i (:,:), 'T', 1. )      ! Lateral boundary condition 
     777 
     778      ! Lateral boundary condition 
     779      CALL lbc_lnk( divu_i (:,:), 'T', 1. ) 
    744780      CALL lbc_lnk( delta_i(:,:), 'T', 1. ) 
     781      ! CALL lbc_lnk( shear_i(:,:), 'F', 1. ) 
    745782      CALL lbc_lnk( shear_i(:,:), 'T', 1. ) 
    746783 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limrst.F90

    r3625 r4161  
    3838 
    3939   !!---------------------------------------------------------------------- 
    40    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     40   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4141   !! $Id$ 
    4242   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    162162      CALL iom_rstput( iter, nitrst, numriw, 'v_ice'     , v_ice      ) 
    163163      CALL iom_rstput( iter, nitrst, numriw, 'fsbbq'     , fsbbq      ) 
     164      CALL iom_rstput( iter, nitrst, numriw, 'iatte'     , iatte      ) ! clem modif 
     165      CALL iom_rstput( iter, nitrst, numriw, 'oatte'     , oatte      ) ! clem modif 
    164166      CALL iom_rstput( iter, nitrst, numriw, 'stress1_i' , stress1_i  ) 
    165167      CALL iom_rstput( iter, nitrst, numriw, 'stress2_i' , stress2_i  ) 
     
    340342      !Control of date 
    341343 
    342       IF( ( nit000 - INT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   & 
     344      IF( ( nit000 - NINT(ziter) ) /= 1 .AND. ABS( nrstdt ) == 1 )   & 
    343345         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nit000 in ice restart',  & 
    344346         &                   '   verify the file or rerun with the value 0 for the',        & 
    345347         &                   '   control of time parameter  nrstdt' ) 
    346       IF( INT(zfice) /= nn_fsbc          .AND. ABS( nrstdt ) == 1 )   & 
     348      IF( NINT(zfice) /= nn_fsbc          .AND. ABS( nrstdt ) == 1 )   & 
    347349         &     CALL ctl_stop( 'lim_rst_read ===>>>> : problem with nn_fsbc in ice restart',  & 
    348350         &                   '   verify the file or rerun with the value 0 for the',         & 
     
    369371         CALL iom_get( numrir, jpdom_autoglo, znam , z2d ) 
    370372         t_su(:,:,jl) = z2d(:,:) 
     373         tn_ice (:,:,:) = t_su (:,:,:) 
    371374      END DO 
    372375 
     
    437440      CALL iom_get( numrir, jpdom_autoglo, 'v_ice'     , v_ice      ) 
    438441      CALL iom_get( numrir, jpdom_autoglo, 'fsbbq'     , fsbbq      ) 
     442      CALL iom_get( numrir, jpdom_autoglo, 'iatte'     , iatte      ) ! clem modif 
     443      CALL iom_get( numrir, jpdom_autoglo, 'oatte'     , oatte      ) ! clem modif 
    439444      CALL iom_get( numrir, jpdom_autoglo, 'stress1_i' , stress1_i  ) 
    440445      CALL iom_get( numrir, jpdom_autoglo, 'stress2_i' , stress2_i  ) 
     
    563568      END DO 
    564569      ! 
    565       CALL iom_close( numrir ) 
     570      !clem CALL iom_close( numrir ) 
    566571      ! 
    567572      CALL wrk_dealloc( nlay_i, zs_zero ) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4148 r4161  
    1010   !!                 !                  + simplification of the ice-ocean stress calculation 
    1111   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
     12   !!             -   ! 2012    (D. Iovino) salt flux change 
     13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux 
    1214   !!            3.5  ! 2012-10 (A. Coward, G. Madec) salt fluxes ; ice+snow mass 
    1315   !!---------------------------------------------------------------------- 
     
    3537   USE prtctl           ! Print control 
    3638   USE cpl_oasis3, ONLY : lk_cpl 
     39   USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget 
    3740   USE oce,        ONLY : sshn, sshb, snwice_mass, snwice_mass_b, snwice_fmass, sshu_b, sshv_b, sshu_n, sshv_n, sshf_n 
    3841   USE dom_ice,    ONLY : tms 
     
    5760#  include "vectopt_loop_substitute.h90" 
    5861   !!---------------------------------------------------------------------- 
    59    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     62   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    6063   !! $Id$ 
    6164   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    99102      INTEGER, INTENT(in) ::   kt    ! number of iteration 
    100103      ! 
    101       INTEGER  ::   ji, jj           ! dummy loop indices 
     104      INTEGER  ::   ji, jj, jl           ! dummy loop indices 
    102105      INTEGER  ::   ierr, ifvt, i1mfr, idfr           ! local integer 
    103106      INTEGER  ::   iflt, ial , iadv , ifral, ifrdv   !   -      - 
     
    106109      REAL(wp) ::   zfcm1 , zfcm2                     !   -      - 
    107110      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
     111      REAL(wp) ::   zzfcm1, zfscmbq ! clem: for light penetration 
    108112      !!--------------------------------------------------------------------- 
    109113       
     
    119123         DO ji = 1, jpi 
    120124            zinda   = 1.0 - MAX( rzero , SIGN( rone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
    121             ifvt    = zinda  *  MAX( rzero , SIGN( rone, -phicif  (ji,jj) ) )  !subscripts are bad here 
    122             i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - ( at_i(ji,jj)       ) ) ) 
     125            ifvt    = zinda  *  MAX( rzero , SIGN( rone, - phicif(ji,jj) ) )  !subscripts are bad here 
     126            i1mfr   = 1.0 - MAX( rzero , SIGN( rone ,  - at_i(ji,jj) ) ) 
    123127            idfr    = 1.0 - MAX( rzero , SIGN( rone , ( 1.0 - at_i(ji,jj) ) - pfrld(ji,jj) ) ) 
    124128            iflt    = zinda  * (1 - i1mfr) * (1 - ifvt ) 
     
    141145 
    142146            !   computation the solar flux at ocean surface 
    143             zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) 
     147            IF (lk_cpl) THEN ! be carfeful: not being tested yet 
     148               ! original line 
     149               !zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) 
     150               ! new line to include solar penetration (not tested) 
     151               zfcm1 = qsr_tot(ji,jj) + fstric(ji,jj) * at_i(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 
     152               DO jl = 1, jpl 
     153                  zfcm1 = zfcm1 - qsr_ice(ji,jj,jl) * a_i(ji,jj,jl) 
     154               END DO 
     155            ELSE 
     156               zfcm1   = pfrld(ji,jj) * qsr(ji,jj)  + & 
     157                    &    ( 1._wp - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 
     158            ENDIF 
    144159            ! fstric     Solar flux transmitted trough the ice 
    145160            ! qsr        Net short wave heat flux on free ocean 
    146161            ! new line 
    147             fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) 
     162            fscmbq(ji,jj) = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) / ( 1.0 - zinda + zinda * iatte(ji,jj) ) 
     163 
     164            ! solar flux and fscmbq with light penetration (clem) 
     165            zzfcm1  = pfrld(ji,jj) * qsr(ji,jj) * oatte(ji,jj) + ( 1. - pfrld(ji,jj) ) * fstric(ji,jj) 
     166            zfscmbq = ( 1.0 - pfrld(ji,jj) ) * fstric(ji,jj) 
    148167 
    149168            !  computation the non solar heat flux at ocean surface 
    150             zfcm2 = - zfcm1                                                                     & ! ??? 
    151                &    + iflt    * fscmbq(ji,jj)                                                   & ! total ablation: heat given to the ocean 
     169            zfcm2 = - zzfcm1                                                                    & ! 
     170               &    + iflt    * zfscmbq                                                         & ! total ablation: heat given to the ocean 
    152171               &    + ifral   * ( ial * qcmif(ji,jj) + (1 - ial) * qldif(ji,jj) ) * r1_rdtice   & 
    153172               &    + ifrdv   * (       qfvbq(ji,jj) +             qdtcn(ji,jj) ) * r1_rdtice   & 
     
    170189            !                           ! fdtcn : turbulent oceanic heat flux 
    171190 
    172 !!gm   this IF prevents the vertorisation of the whole loop 
    173             IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
    174                WRITE(numout,*) ' lim_sbc : heat fluxes ' 
    175                WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx) 
    176                WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx) 
    177                WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx) 
    178                WRITE(numout,*) 
    179                WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx) 
    180                WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    181                WRITE(numout,*) ' ifral     : ', ifral 
    182                WRITE(numout,*) ' ial       : ', ial   
    183                WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx) 
    184                WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx) 
    185                WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 
    186                WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 
    187                WRITE(numout,*) ' ifrdv     : ', ifrdv 
    188                WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx) 
    189                WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx) 
    190                WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 
    191                WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 
    192                WRITE(numout,*) ' ' 
    193                WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
    194                WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx) 
    195                WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx) 
    196                WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx) 
    197                WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
    198             ENDIF 
    199 !!gm   end 
     191            !!gm   this IF prevents the vertorisation of the whole loop 
     192          !  IF ( ( ji == jiindx ) .AND. ( jj == jjindx) ) THEN 
     193          !     WRITE(numout,*) ' lim_sbc : heat fluxes ' 
     194          !     WRITE(numout,*) ' qsr       : ', qsr(jiindx,jjindx) 
     195          !     WRITE(numout,*) ' pfrld     : ', pfrld(jiindx,jjindx) 
     196          !     WRITE(numout,*) ' fstric    : ', fstric (jiindx,jjindx) 
     197          !     WRITE(numout,*) 
     198          !     WRITE(numout,*) ' qns       : ', qns(jiindx,jjindx) 
     199          !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
     200          !     WRITE(numout,*) ' ifral     : ', ifral 
     201          !     WRITE(numout,*) ' ial       : ', ial   
     202          !     WRITE(numout,*) ' qcmif     : ', qcmif(jiindx,jjindx) 
     203          !     WRITE(numout,*) ' qldif     : ', qldif(jiindx,jjindx) 
     204          !     !WRITE(numout,*) ' qcmif / dt: ', qcmif(jiindx,jjindx) * r1_rdtice 
     205          !     !WRITE(numout,*) ' qldif / dt: ', qldif(jiindx,jjindx) * r1_rdtice 
     206          !     WRITE(numout,*) ' ifrdv     : ', ifrdv 
     207          !     WRITE(numout,*) ' qfvbq     : ', qfvbq(jiindx,jjindx) 
     208          !     WRITE(numout,*) ' qdtcn     : ', qdtcn(jiindx,jjindx) 
     209          !     !WRITE(numout,*) ' qfvbq / dt: ', qfvbq(jiindx,jjindx) * r1_rdtice 
     210          !     !WRITE(numout,*) ' qdtcn / dt: ', qdtcn(jiindx,jjindx) * r1_rdtice 
     211          !     WRITE(numout,*) ' ' 
     212          !     WRITE(numout,*) ' fdtcn     : ', fdtcn(jiindx,jjindx) 
     213          !     WRITE(numout,*) ' fhmec     : ', fhmec(jiindx,jjindx) 
     214          !     WRITE(numout,*) ' fheat_mec : ', fheat_mec(jiindx,jjindx) 
     215          !     WRITE(numout,*) ' fhbri     : ', fhbri(jiindx,jjindx) 
     216          !     WRITE(numout,*) ' fheat_res : ', fheat_res(jiindx,jjindx) 
     217          !  ENDIF 
     218            !!gm   end 
    200219         END DO 
    201220      END DO 
     
    218237 
    219238            !  computing freshwater exchanges at the ice/ocean interface 
    220             zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
    221                &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
    222                &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
    223                &   - fmmec(ji,jj)                                         ! snow falling when ridging 
     239            IF (lk_cpl) THEN  
     240               zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
     241                  &   - rdm_snw(ji,jj) / rdt_ice 
     242            ELSE 
     243               zemp =   emp(ji,jj)     * ( 1.0 - at_i(ji,jj)          )  &   ! evaporation over oceanic fraction 
     244                  &   - tprecip(ji,jj) *         at_i(ji,jj)             &   ! all precipitation reach the ocean 
     245                  &   + sprecip(ji,jj) * ( 1. - (pfrld(ji,jj)**betas) )  &   ! except solid precip intercepted by sea-ice 
     246                  &   - fmmec(ji,jj)                                         ! snow falling when ridging 
     247            ENDIF 
    224248 
    225249            ! mass flux at the ocean/ice interface (sea ice fraction) 
     
    370394      !! ** input   : Namelist namicedia 
    371395      !!------------------------------------------------------------------- 
     396      REAL(wp) :: zsum, zarea 
    372397      ! 
    373398      INTEGER  ::   ji, jj                          ! dummy loop indices 
     
    390415         END WHERE 
    391416      ENDIF 
     417      ! clem modif 
     418      iatte(:,:) = 1._wp 
     419      oatte(:,:) = 1._wp 
     420      ! 
    392421      !                                      ! embedded sea ice 
    393422      IF( nn_ice_embd /= 0 ) THEN            ! mass exchanges between ice and ocean (case 1 or 2) set the snow+ice mass 
     
    435464      ENDIF 
    436465      ! 
     466!!?      IF( .NOT. ln_rstart ) THEN           ! delete the initial ssh below sea-ice area 
     467!!?         ! 
     468!!?         zarea     = glob_sum( e1e2t(:,:) )           ! interior global domain surface 
     469!!?         zsum      = glob_sum( e1e2t(:,:) * ( snwice_mass(:,:) ) ) / zarea * r1_rau0 
     470!!?         sshn(:,:) = sshn(:,:) - zsum  
     471!!?         sshb(:,:) = sshb(:,:) - zsum 
     472!!?      ENDIF 
     473      ! 
     474 
    437475   END SUBROUTINE lim_sbc_init 
    438476 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limtab.F90

    r3625 r4161  
    2020 
    2121   !!---------------------------------------------------------------------- 
    22    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2010) 
     22   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    2323   !! $Id$ 
    2424   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r4147 r4161  
    1111   !!            3.3  ! 2010-11 (G. Madec) corrected snow melting heat (due to factor betas) 
    1212   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     13   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux 
    1314   !!---------------------------------------------------------------------- 
    1415#if defined key_lim3 
     
    4041   USE prtctl         ! Print control 
    4142   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     43   USE timing         ! Timing 
    4244 
    4345   IMPLICIT NONE 
     
    9294      REAL(wp) ::   zfntlat, zpareff, zareamin, zcoef   !    -         - 
    9395      REAL(wp), POINTER, DIMENSION(:,:) ::   zqlbsbq   ! link with lead energy budget qldif 
     96      REAL(wp) :: zchk_v_i, zchk_smv, zchk_fs, zchk_fw, zchk_v_i_b, zchk_smv_b, zchk_fs_b, zchk_fw_b ! Check conservation (C Rousset) 
     97      REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    9498      !!------------------------------------------------------------------- 
     99      IF( nn_timing == 1 )  CALL timing_start('limthd') 
    95100 
    96101      CALL wrk_alloc( jpi, jpj, zqlbsbq ) 
    97102    
     103      ! ------------------------------- 
     104      !- check conservation (C Rousset) 
     105      IF (ln_limdiahsb) THEN 
     106         zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     107         zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     108         zchk_fw_b  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) 
     109         zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) 
     110      ENDIF 
     111      !- check conservation (C Rousset) 
     112      ! ------------------------------- 
     113 
    98114      !------------------------------------------------------------------------------! 
    99115      ! 1) Initialization of diagnostic variables                                    ! 
     
    109125               DO ji = 1, jpi 
    110126                  !Energy of melting q(S,T) [J.m-3] 
    111                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * nlay_i 
     127                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_i(ji,jj,jl) , epsi06 ) ) * REAL( nlay_i ) 
    112128                  !0 if no ice and 1 if yes 
    113129                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_i(ji,jj,jl) )  )  
     
    121137               DO ji = 1, jpi 
    122138                  !Energy of melting q(S,T) [J.m-3] 
    123                   e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * nlay_s 
     139                  e_s(ji,jj,jk,jl) = e_s(ji,jj,jk,jl) / ( area(ji,jj) * MAX( v_s(ji,jj,jl) , epsi06 ) ) * REAL( nlay_s ) 
    124140                  !0 if no ice and 1 if yes 
    125141                  zindb = 1.0 - MAX(  0.0 , SIGN( 1.0 , - ht_s(ji,jj,jl) )  )  
     
    134150      ! 1.3) Set some dummies to 0 
    135151      !----------------------------- 
    136       rdvosif(:,:) = 0.e0   ! variation of ice volume at surface 
    137       rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom 
    138       fdvolif(:,:) = 0.e0   ! total variation of ice volume 
    139       rdvonif(:,:) = 0.e0   ! lateral variation of ice volume 
    140       fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice 
    141       ffltbif(:,:) = 0.e0   ! linked with fstric 
    142       qfvbq  (:,:) = 0.e0   ! linked with fstric 
    143       rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area 
    144       rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area 
    145       hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.  
    146       sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean 
    147       fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean 
    148       sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay 
     152      !clem rdvosif(:,:) = 0.e0   ! variation of ice volume at surface 
     153      !clem rdvobif(:,:) = 0.e0   ! variation of ice volume at bottom 
     154      !clem fdvolif(:,:) = 0.e0   ! total variation of ice volume 
     155      !clem rdvonif(:,:) = 0.e0   ! lateral variation of ice volume 
     156      !clem fstric (:,:) = 0.e0   ! part of solar radiation transmitted through the ice 
     157      !clem ffltbif(:,:) = 0.e0   ! linked with fstric 
     158      !clem qfvbq  (:,:) = 0.e0   ! linked with fstric 
     159      !clem rdm_snw(:,:) = 0.e0   ! variation of snow mass per unit area 
     160      !clem rdm_ice(:,:) = 0.e0   ! variation of ice mass per unit area 
     161      !clem hicifp (:,:) = 0.e0   ! daily thermodynamic ice production.  
     162      !clem sfx_bri(:,:) = 0.e0   ! brine flux contribution to salt flux to the ocean 
     163      !clem fhbri  (:,:) = 0.e0   ! brine flux contribution to heat flux to the ocean 
     164      !clem sfx_thd(:,:) = 0.e0   ! equivalent salt flux to the ocean due to ice/growth decay 
    149165 
    150166      !----------------------------------- 
     
    165181!CDIR NOVERRCHK 
    166182         DO ji = 1, jpi 
    167             zthsnice       = SUM( ht_s(ji,jj,1:jpl) ) + SUM( ht_i(ji,jj,1:jpl) ) 
    168             zindb          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - zthsnice ) ) )  
    169183            phicif(ji,jj)  = vt_i(ji,jj) 
    170184            pfrld(ji,jj)   = 1.0 - at_i(ji,jj) 
    171             zinda          = 1.0 - MAX( zzero , SIGN( zone , - ( 1.0 - pfrld(ji,jj) ) ) ) 
     185            zinda          = tms(ji,jj) * ( 1.0 - MAX( zzero , SIGN( zone , - at_i(ji,jj) ) ) ) 
    172186            ! 
    173187            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
     
    180194 
    181195            ! here the drag will depend on ice thickness and type (0.006) 
    182             fdtcn(ji,jj)  = zindb * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )  
     196            fdtcn(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( (sst_m(ji,jj) + rt0) - t_bo(ji,jj) )  
    183197            ! also category dependent 
    184198            !           !-- Energy from the turbulent oceanic heat flux heat flux coming in the lead  
    185             qdtcn(ji,jj)  = zindb * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 
     199            qdtcn(ji,jj)  = zinda * fdtcn(ji,jj) * (1.0 - at_i(ji,jj)) * rdt_ice 
    186200            !                        
    187201            !           !-- Lead heat budget, qldif (part 1, next one is in limthd_dh)  
    188202            !           !   caution: exponent betas used as more snow can fallinto leads 
    189203            qldif(ji,jj) =  tms(ji,jj) * rdt_ice  * (                             & 
    190                &   pfrld(ji,jj)        * (  qsr(ji,jj)                            &   ! solar heat 
     204               &   pfrld(ji,jj)        * (  qsr(ji,jj) * oatte(ji,jj)             &   ! solar heat + clem modif 
    191205               &                            + qns(ji,jj)                          &   ! non solar heat 
    192206               &                            + fdtcn(ji,jj)                        &   ! turbulent ice-ocean heat 
    193                &                            + fsbbq(ji,jj) * ( 1.0 - zindb )  )   &   ! residual heat from previous step 
     207               &                            + fsbbq(ji,jj) * ( 1.0 - zinda )  )   &   ! residual heat from previous step 
    194208               & - pfrld(ji,jj)**betas * sprecip(ji,jj) * lfus                    )   ! latent heat of sprecip melting 
    195209            ! 
     
    206220            ! 
    207221            ! Energy needed to bring ocean surface layer until its freezing (qcmif, limflx) 
    208             qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) * ( 1. - zinda ) 
     222            qcmif  (ji,jj) =  rau0 * rcp * fse3t(ji,jj,1) * ( t_bo(ji,jj) - (sst_m(ji,jj) + rt0) ) 
    209223            ! 
    210224            ! oceanic heat flux (limthd_dh) 
    211             fbif   (ji,jj) = zindb * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 
     225            fbif   (ji,jj) = zinda * (  fsbbq(ji,jj) / MAX( at_i(ji,jj) , epsi20 ) + fdtcn(ji,jj) ) 
    212226            ! 
    213227         END DO 
     
    294308            CALL tab_2d_1d( nbpb, qfvbq_1d   (1:nbpb), qfvbq           , jpi, jpj, npb(1:nbpb) ) 
    295309 
     310            CALL tab_2d_1d( nbpb, iatte_1d   (1:nbpb), iatte           , jpi, jpj, npb(1:nbpb) ) ! clem modif 
     311            CALL tab_2d_1d( nbpb, oatte_1d   (1:nbpb), oatte           , jpi, jpj, npb(1:nbpb) ) ! clem modif 
    296312            !-------------------------------- 
    297313            ! 4.3) Thermodynamic processes 
     
    411427      ! 5.4) Diagnostic thermodynamic growth rates 
    412428      !-------------------------------------------- 
    413       d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes  
    414       dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
     429!clem@useless      d_v_i_thd(:,:,:) = v_i      (:,:,:) - old_v_i(:,:,:)    ! ice volumes  
     430!clem@mv-to-itd    dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    415431 
    416432      IF( con_i )   fbif(:,:) = fbif(:,:) + zqlbsbq(:,:) 
     
    448464      ENDIF 
    449465      ! 
     466      ! ------------------------------- 
     467      !- check conservation (C Rousset) 
     468      IF (ln_limdiahsb) THEN 
     469         zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_thd(:,:) + sfx_res(:,:) + sfx_mec(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
     470         zchk_fw  = glob_sum( rdm_ice(:,:) * area(:,:) * tms(:,:) ) - zchk_fw_b 
     471  
     472         zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b - ( zchk_fw / rhoic ) ) * r1_rdtice 
     473         zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
     474 
     475         zchk_vmin = glob_min(v_i) 
     476         zchk_amax = glob_max(SUM(a_i,dim=3)) 
     477         zchk_amin = glob_min(a_i) 
     478        
     479         IF(lwp) THEN 
     480            IF ( ABS( zchk_v_i   ) >  1.e-5 ) WRITE(numout,*) 'violation volume [m3/day]     (limthd) = ',(zchk_v_i * rday) 
     481            IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 
     482            IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(zchk_vmin * 1.e-3) 
     483            IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limthd) = ',zchk_amax 
     484            IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limthd) = ',zchk_amin 
     485         ENDIF 
     486      ENDIF 
     487      !- check conservation (C Rousset) 
     488      ! ------------------------------- 
     489      ! 
    450490      CALL wrk_dealloc( jpi, jpj, zqlbsbq ) 
    451491      ! 
     492      IF( nn_timing == 1 )  CALL timing_stop('limthd') 
    452493   END SUBROUTINE lim_thd 
    453494 
     
    472513      DO jk = 1, nlay_i                ! total q over all layers, ice [J.m-2] 
    473514         DO ji = kideb, kiut 
    474             etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 
     515            etilayer(ji,jk) = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 
    475516            eti     (ji,jl) = eti(ji,jl) + etilayer(ji,jk)  
    476517         END DO 
    477518      END DO 
    478519      DO ji = kideb, kiut              ! total q over all layers, snow [J.m-2] 
    479          ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / nlay_s 
     520         ets(ji,jl) = ets(ji,jl) + q_s_b(ji,1) * ht_s_b(ji) / REAL( nlay_s ) 
    480521      END DO 
    481522      ! 
     
    498539 
    499540      INTEGER  ::   ji, jk         ! loop indices 
    500       INTEGER  ::   zji, zjj 
     541      INTEGER  ::   ii, ij 
    501542      INTEGER  ::   numce          ! number of points for which conservation is violated 
    502543      REAL(wp) ::   meance         ! mean conservation error 
     
    521562      !---------------------------------------- 
    522563      DO ji = kideb, kiut 
    523          zji = MOD( npb(ji) - 1 , jpi ) + 1 
    524          zjj =    ( npb(ji) - 1 ) / jpi + 1 
     564         ii = MOD( npb(ji) - 1 , jpi ) + 1 
     565         ij =    ( npb(ji) - 1 ) / jpi + 1 
    525566         fatm     (ji,jl) = qnsr_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) 
    526          sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(zji,zjj,jl) 
     567         sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji) + qsr_ice_1d(ji) * i0(ji) - fstroc(ii,ij,jl) 
    527568      END DO 
    528569 
     
    579620         IF ( ( ( t_su_b(ji) .LT. rtt ) .AND. ( surf_error(ji,jl) .GT. max_surf_err ) ) .OR. & 
    580621            ( cons_error(ji,jl) .GT. max_cons_err  ) ) THEN 
    581             zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    582             zjj                 = ( npb(ji) - 1 ) / jpi + 1 
     622            ii                 = MOD( npb(ji) - 1, jpi ) + 1 
     623            ij                 = ( npb(ji) - 1 ) / jpi + 1 
    583624            ! 
    584625            WRITE(numout,*) ' alerte 1     ' 
     
    586627            WRITE(numout,*) ' heat diffusion in the ice ' 
    587628            WRITE(numout,*) ' Category   : ', jl 
    588             WRITE(numout,*) ' zji , zjj  : ', zji, zjj 
    589             WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj) 
     629            WRITE(numout,*) ' ii , ij  : ', ii, ij 
     630            WRITE(numout,*) ' lat, lon   : ', gphit(ii,ij), glamt(ii,ij) 
    590631            WRITE(numout,*) ' cons_error : ', cons_error(ji,jl) 
    591632            WRITE(numout,*) ' surf_error : ', surf_error(ji,jl) 
     
    615656            WRITE(numout,*) ' fc_bo      : ', - fc_bo_i  (ji) 
    616657            WRITE(numout,*) ' foc        : ', fbif_1d(ji) 
    617             WRITE(numout,*) ' fstroc     : ', fstroc   (zji,zjj,jl) 
     658            WRITE(numout,*) ' fstroc     : ', fstroc   (ii,ij,jl) 
    618659            WRITE(numout,*) ' i0         : ', i0(ji) 
    619660            WRITE(numout,*) ' qsr_ice    : ', (1.0-i0(ji))*qsr_ice_1d(ji) 
     
    651692      ! 
    652693      INTEGER  ::   ji                ! loop indices 
    653       INTEGER  ::   zji, zjj, numce         ! local integers 
     694      INTEGER  ::   ii, ij, numce         ! local integers 
    654695      REAL(wp) ::   meance, max_cons_err    !local scalar 
    655696      !!--------------------------------------------------------------------- 
     
    669710      !---------------------------------------- 
    670711      DO ji = kideb, kiut 
    671          zji = MOD( npb(ji) - 1 , jpi ) + 1 
    672          zjj =    ( npb(ji) - 1 ) / jpi + 1 
     712         ii = MOD( npb(ji) - 1 , jpi ) + 1 
     713         ij =    ( npb(ji) - 1 ) / jpi + 1 
    673714 
    674715         fatm      (ji,jl) = qnsr_ice_1d(ji) + qsr_ice_1d(ji)                       ! total heat flux 
    675          sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(zji,zjj,jl)  
     716         sum_fluxq (ji,jl) = fatm(ji,jl) + fbif_1d(ji) - ftotal_fin(ji) - fstroc(ii,ij,jl)  
    676717         cons_error(ji,jl) = ABS( dq_i(ji,jl) * r1_rdtice + sum_fluxq(ji,jl) ) 
    677718      END DO 
     
    706747      DO ji = kideb, kiut 
    707748         IF ( cons_error(ji,jl) .GT. max_cons_err  ) THEN 
    708             zji = MOD( npb(ji) - 1, jpi ) + 1 
    709             zjj =    ( npb(ji) - 1 ) / jpi + 1 
     749            ii = MOD( npb(ji) - 1, jpi ) + 1 
     750            ij =    ( npb(ji) - 1 ) / jpi + 1 
    710751            ! 
    711752            WRITE(numout,*) ' alerte 1 - category : ', jl 
    712753            WRITE(numout,*) ' Untolerated conservation error after limthd_ent ' 
    713             WRITE(numout,*) ' zji , zjj  : ', zji, zjj 
    714             WRITE(numout,*) ' lat, lon   : ', gphit(zji,zjj), glamt(zji,zjj) 
     754            WRITE(numout,*) ' ii , ij  : ', ii, ij 
     755            WRITE(numout,*) ' lat, lon   : ', gphit(ii,ij), glamt(ii,ij) 
    715756            WRITE(numout,*) ' * ' 
    716757            WRITE(numout,*) ' Ftotal     : ', sum_fluxq(ji,jl) 
     
    724765            WRITE(numout,*) ' foce       : ', fbif_1d(ji) 
    725766            WRITE(numout,*) ' fres       : ', ftotal_fin(ji) 
    726             WRITE(numout,*) ' fhbri      : ', fhbricat(zji,zjj,jl) 
     767            WRITE(numout,*) ' fhbri      : ', fhbricat(ii,ij,jl) 
    727768            WRITE(numout,*) ' * ' 
    728769            WRITE(numout,*) ' Heat contents --- : ' 
     
    793834      INTEGER  ::   ios                 ! Local integer output status for namelist read 
    794835      NAMELIST/namicethd/ hmelt , hiccrit, fraz_swi, maxfrazb, vfrazb, Cfrazb,   & 
    795          &                hicmin, hiclim, amax  ,                                & 
     836         &                hicmin, hiclim,                                        & 
    796837         &                sbeta  , parlat, hakspl, hibspl, exld,                 & 
    797838         &                hakdif, hnzst  , thth  , parsub, alphs, betas,         &  
     
    825866         WRITE(numout,*)'      ice thick. corr. to max. energy stored in brine pocket  hicmin       = ', hicmin   
    826867         WRITE(numout,*)'      minimum ice thickness                                   hiclim       = ', hiclim  
    827          WRITE(numout,*)'      maximum lead fraction                                   amax         = ', amax  
    828868         WRITE(numout,*)'      numerical carac. of the scheme for diffusion in ice ' 
    829869         WRITE(numout,*)'      Cranck-Nicholson (=0.5), implicit (=1), explicit (=0)   sbeta        = ', sbeta 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r3808 r4161  
    66   !! History :  LIM  ! 2003-05 (M. Vancoppenolle) Original code in 1D 
    77   !!                 ! 2005-06 (M. Vancoppenolle) 3D version  
    8    !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdmsnif & rdmicif 
     8   !!            3.2  ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in rdm_snw & rdm_ice 
    99   !!            3.4  ! 2011-02 (G. Madec) dynamical allocation 
    1010   !!            3.5  ! 2012-10 (G. Madec & co) salt flux + bug fixes  
     
    3939 
    4040   !!---------------------------------------------------------------------- 
    41    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     41   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010) 
    4242   !! $Id$ 
    4343   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7373      !!  
    7474      INTEGER  ::   ji , jk        ! dummy loop indices 
    75       INTEGER  ::   ii, ij       ! 2D corresponding indices to ji 
     75      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
    7676      INTEGER  ::   isnow          ! switch for presence (1) or absence (0) of snow 
    7777      INTEGER  ::   isnowic        ! snow ice formation not 
     
    8484      REAL(wp) ::   zdhnm, zhnnew, zhisn, zihic, zzc       ! 
    8585      REAL(wp) ::   zfracs       ! fractionation coefficient for bottom salt entrapment 
    86       REAL(wp) ::   zds          ! increment of bottom ice salinity 
    8786      REAL(wp) ::   zcoeff       ! dummy argument for snowfall partitioning over ice and leads 
    8887      REAL(wp) ::   zsm_snowice  ! snow-ice salinity 
     
    107106      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_pre   ! snow precipitation  
    108107      REAL(wp), POINTER, DIMENSION(:) ::   zdh_s_sub   ! snow sublimation 
    109       REAL(wp), POINTER, DIMENSION(:) ::   zsfx_melt   ! salt flux due to ice melt 
    110108 
    111109      REAL(wp), POINTER, DIMENSION(:,:) ::   zdeltah 
     
    120118      REAL(wp), POINTER, DIMENSION(:,:) ::   zqt_i_lay   ! total ice heat content 
    121119 
     120      ! mass and salt flux (clem) 
     121      REAL(wp) :: zdvres, zdvsur, zdvbot 
     122      REAL(wp), POINTER, DIMENSION(:) ::   zviold, zvsold   ! old ice volume... 
     123 
    122124      ! Heat conservation  
    123125      INTEGER  ::   num_iter_max, numce_dh 
     
    128130 
    129131      CALL wrk_alloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    130       CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zsfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
     132      CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    131133      CALL wrk_alloc( jpij, zinnermelt, zfbase, zdq_i ) 
    132134      CALL wrk_alloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
    133135 
    134       zsfx_melt (:) = 0._wp 
     136      CALL wrk_alloc( jpij, zviold, zvsold ) ! clem 
     137       
    135138      ftotal_fin(:) = 0._wp 
    136139      zfdt_init (:) = 0._wp 
    137140      zfdt_final(:) = 0._wp 
    138141 
     142      dh_i_surf (:) = 0._wp 
     143      dh_i_bott (:) = 0._wp 
     144      dh_snowice(:) = 0._wp 
     145 
    139146      DO ji = kideb, kiut 
    140147         old_ht_i_b(ji) = ht_i_b(ji) 
    141148         old_ht_s_b(ji) = ht_s_b(ji) 
     149         zviold(ji) = a_i_b(ji) * ht_i_b(ji) ! clem 
     150         zvsold(ji) = a_i_b(ji) * ht_s_b(ji) ! clem 
    142151      END DO 
    143152      ! 
     
    164173      ! 
    165174      DO ji = kideb, kiut     ! Layer thickness 
    166          zh_i(ji) = ht_i_b(ji) / nlay_i 
    167          zh_s(ji) = ht_s_b(ji) / nlay_s 
     175         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
     176         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 
    168177      END DO 
    169178      ! 
     
    171180      DO jk = 1, nlay_s 
    172181         DO ji = kideb, kiut 
    173             zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s 
     182            zqt_s(ji) =  zqt_s(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s ) 
    174183         END DO 
    175184      END DO 
     
    178187      DO jk = 1, nlay_i 
    179188         DO ji = kideb, kiut 
    180             zzc = q_i_b(ji,jk) * ht_i_b(ji) / nlay_i 
     189            zzc = q_i_b(ji,jk) * ht_i_b(ji) / REAL( nlay_i ) 
    181190            zqt_i(ji)        =  zqt_i(ji) + zzc 
    182191            zqt_i_lay(ji,jk) =              zzc 
     
    244253         zhn            =  1.0 - MAX(  zzero , SIGN( zone , - zhsnew )  ) 
    245254         ht_s_b(ji)     =  MAX( zzero , zhsnew ) 
     255         ! we recompute dh_s_tot (clem)  
     256         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    246257         ! Volume and mass variations of snow 
    247258         dvsbq_1d  (ji) =  a_i_b(ji) * ( ht_s_b(ji) - zhsold(ji) - zdh_s_pre(ji) ) 
    248259         dvsbq_1d  (ji) =  MIN( zzero, dvsbq_1d(ji) ) 
    249          rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
     260         !clem rdm_snw_1d(ji) =  rdm_snw_1d(ji) + rhosn * dvsbq_1d(ji) 
    250261      END DO ! ji 
    251262 
     
    254265      !-------------------------- 
    255266      DO ji = kideb, kiut  
    256          dh_i_surf(ji) =  0._wp 
    257267         z_f_surf (ji) =  zqfont_su(ji) * r1_rdtice   ! heat conservation test 
    258268         zdq_i    (ji) =  0._wp 
     
    272282            zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    273283            ! 
    274             !                                                    ! contribution to ice-ocean salt flux  
    275             zsfx_melt(ji)  = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
     284            ! clem 
     285            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
     286               &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic / rdt_ice 
    276287         END DO 
    277288      END DO 
     
    334345         DO ji = kideb,kiut 
    335346            q_s_b    (ji,jk) = rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) 
    336             zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / nlay_s            ! heat conservation 
     347            zqt_dummy(ji)    =  zqt_dummy(ji) + q_s_b(ji,jk) * ht_s_b(ji) / REAL( nlay_s )            ! heat conservation 
    337348         END DO 
    338349      END DO 
     
    375386               ! Basal growth rate = - F*dt / q 
    376387               dh_i_bott(ji)       =  - rdt_ice * ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) / q_i_b(ji,nlay_i+1)  
     388               sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    377389            ENDIF 
    378390         END DO 
     
    416428                  zfracs = zswi1  * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) )   & 
    417429                     &                   + zswi2  * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) )  
    418                   zds         = zfracs * sss_m(ii,ij) - s_i_new(ji) 
     430                  zfracs = MIN( 0.5 , zfracs ) 
    419431                  s_i_new(ji) = zfracs * sss_m(ii,ij) 
    420432               ENDIF ! fc_bo_i 
     
    425437         DO ji = kideb, kiut 
    426438            IF( ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) ) .LT. 0.0  ) THEN 
    427                ! New ice salinity must not exceed 15 psu 
     439               ! New ice salinity must not exceed 20 psu 
    428440               s_i_new(ji) = MIN( s_i_new(ji), s_i_max ) 
    429441               ! Metling point in K 
     
    437449               ! Salinity update 
    438450               ! entrapment during bottom growth 
    439                dsm_i_se_1d(ji) = ( s_i_new(ji) * dh_i_bott(ji) + sm_i_b(ji) * ht_i_b(ji) )    & 
    440                   &            / MAX( ht_i_b(ji) + dh_i_bott(ji) ,epsi13 ) - sm_i_b(ji) 
     451               sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * dh_i_bott(ji) * rhoic * r1_rdtice 
    441452            ENDIF ! heat budget 
    442453         END DO 
     
    476487                  zdq_i    (ji   ) = zdq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk) * r1_rdtice 
    477488               ENDIF 
    478                ! contribution to salt flux 
    479                zsfx_melt(ji) = zsfx_melt(ji) - sm_i_b(ji) * a_i_b(ji) * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice  
     489               ! clem: contribution to salt flux 
     490               sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji)    & 
     491                    &                              * MIN( zdeltah(ji,jk) , 0._wp ) * rhoic * r1_rdtice 
    480492            ENDIF 
    481493         END DO ! ji 
     
    528540         ELSE                  ;   zdhbf =              dh_i_bott(ji)  
    529541         ENDIF 
     542         zdvres        = zdhbf - dh_i_bott(ji) 
     543         dh_i_bott(ji) = zdhbf 
     544         sfx_thd_1d(ji)  = sfx_thd_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice 
    530545         !                     ! excessive energy is sent to lateral ablation 
    531          fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 )   & 
    532             &                          * ( zdhbf - dh_i_bott(ji) ) * r1_rdtice 
    533          dh_i_bott(ji)  = zdhbf 
    534          !                     !since ice volume is only used for outputs, we keep it global for all categories 
    535          dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
    536          !                     !new ice thickness 
    537          zhgnew   (ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
    538          !                     ! diagnostic ( bottom ice growth ) 
    539          ii = MOD( npb(ji) - 1, jpi ) + 1 
    540          ij = ( npb(ji) - 1 ) / jpi + 1 
    541          diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    542          diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
    543          diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
     546         fsup     (ji) =  rhoic * lfus * at_i_b(ji) / MAX( 1.0 - at_i_b(ji) , epsi13 ) * zdvres * r1_rdtice 
    544547      END DO 
    545548 
     
    552555         ! Adapt the remaining energy if too much ice melts 
    553556         !-------------------------------------------------- 
    554          zihgnew =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    555          ! 0 if no more ice 
    556          zhgnew    (ji) =         zihgnew   * zhgnew(ji)      ! ice thickness is put to 0 
    557          ! remaining heat 
     557         zdvres     = MAX( 0._wp, - ht_i_b(ji) - dh_i_surf(ji) - dh_i_bott(ji) ) 
     558         zdvsur     = MIN( 0._wp, dh_i_surf(ji) + zdvres ) - dh_i_surf(ji) ! fill the surface first 
     559         zdvbot     = MAX( 0._wp, zdvres - zdvsur ) ! then the bottom 
     560         dh_i_surf (ji) = dh_i_surf(ji) + zdvsur ! clem 
     561         dh_i_bott (ji) = dh_i_bott(ji) + zdvbot ! clem 
     562 
     563         ! new ice thickness (clem) 
     564         zhgnew(ji) = ht_i_b(ji) + dh_i_surf(ji) + dh_i_bott(ji) 
     565         zihgnew    = 1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) ) !1 if ice 
     566         zhgnew(ji) = zihgnew * zhgnew(ji)      ! ice thickness is put to 0 
     567  
     568         !                     !since ice volume is only used for outputs, we keep it global for all categories 
     569         dvbbq_1d (ji) = a_i_b(ji) * dh_i_bott(ji) 
     570 
     571        ! remaining heat 
    558572         zfdt_final(ji) = ( 1.0 - zihgnew ) * ( zqfont_su(ji) +  zqfont_bo(ji) )  
    559573 
     
    569583         ht_s_b(ji)     =  MAX( zzero , zhnfi ) 
    570584         zqt_s(ji)      =  zqt_s(ji) * ht_s_b(ji) 
     585         ! we recompute dh_s_tot (clem) 
     586         dh_s_tot (ji)  =  ht_s_b(ji) - zhsold(ji) 
    571587 
    572588         ! Mass variations of ice and snow 
     
    579595         ! 
    580596         !                                              ! mass variation cumulated over category 
    581          rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
    582          rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
     597         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + zzfmass_s                     ! snow  
     598         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + zzfmass_i                     ! ice  
    583599 
    584600         ! Remaining heat to the ocean  
     
    586602         focea(ji)  = - zfdt_final(ji) * r1_rdtice         ! focea is in W.m-2 * dt 
    587603 
     604         ! residual salt flux (clem) 
     605         !-------------------------- 
     606         ! surface 
     607         sfx_thd_1d(ji)    = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvsur * rhoic * r1_rdtice 
     608         ! bottom 
     609         IF ( fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji) >= 0._wp ) THEN ! melting 
     610            sfx_thd_1d(ji) = sfx_thd_1d(ji) - sm_i_b(ji)  * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 
     611         ELSE                                                          ! growth 
     612            sfx_thd_1d(ji) = sfx_thd_1d(ji) - s_i_new(ji) * a_i_b(ji) * zdvbot * rhoic * r1_rdtice 
     613         ENDIF 
     614         ! 
     615         ! diagnostic ( bottom ice growth ) 
     616         ii = MOD( npb(ji) - 1, jpi ) + 1 
     617         ij = ( npb(ji) - 1 ) / jpi + 1 
     618         diag_bot_gr(ii,ij) = diag_bot_gr(ii,ij) + MAX(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
     619         diag_sur_me(ii,ij) = diag_sur_me(ii,ij) + MIN(dh_i_surf(ji),0.0)*a_i_b(ji) * r1_rdtice 
     620         diag_bot_me(ii,ij) = diag_bot_me(ii,ij) + MIN(dh_i_bott(ji),0.0)*a_i_b(ji) * r1_rdtice 
    588621      END DO 
    589622 
     
    591624 
    592625      !--------------------------- 
    593       ! Salt flux and heat fluxes                     
     626      ! heat fluxes                     
    594627      !--------------------------- 
    595628      DO ji = kideb, kiut 
    596629         zihgnew    =  1.0 - MAX( zzero , SIGN( zone , - zhgnew(ji) ) )   ! =1 if ice 
    597          ! 
    598          ! Salt flux 
    599          sfx_thd_1d(ji) = sfx_thd_1d(ji) +        zihgnew  * zsfx_melt(ji)               & 
    600             &                            - (1.0 - zihgnew) * zfmass_i (ji) * sm_i_b(ji)  * r1_rdtice 
    601630         ! 
    602631         ! Heat flux 
     
    646675         dmgwi_1d  (ji) = dmgwi_1d(ji) + a_i_b(ji) * ( ht_s_b(ji) - zhnnew ) * rhosn 
    647676 
    648          ! All snow is thrown in the ocean, and seawater is taken to replace the volume 
    649          rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic * ( 1. - rhosn / rhoic ) 
    650          rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn 
     677         !clem rdm_ice_1d(ji) = rdm_ice_1d(ji) + a_i_b(ji) * ( zhgnew(ji) - ht_i_b(ji) ) * rhoic  
     678         !clem rdm_snw_1d(ji) = rdm_snw_1d(ji) + a_i_b(ji) * ( zhnnew     - ht_s_b(ji) ) * rhosn  
    651679 
    652680         !        Equivalent salt flux (1) Snow-ice formation component 
     
    658686         ELSE                      ;   zsm_snowice = sm_i_b(ji)    
    659687         ENDIF 
    660          sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    661          ! 
    662688         ! entrapment during snow ice formation 
    663          i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - ht_i_b(ji) + 1.0e-6 ) ) 
    664          isnowic      = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - dh_snowice(ji)      ) ) * i_ice_switch 
    665          IF(  num_sal == 2  )   & 
    666             dsm_i_si_1d(ji) = (  zsm_snowice * dh_snowice(ji)    & 
    667             &                  + sm_i_b(ji) * ht_i_b(ji) / MAX( ht_i_b(ji) + dh_snowice(ji), epsi13 )   & 
    668             &                  - sm_i_b(ji)  ) * isnowic      
     689         ! clem: new salinity difference stored (to be used in limthd_ent.F90) 
     690         IF (  num_sal == 2  ) THEN 
     691            i_ice_switch = 1.0 - MAX( 0.e0 , SIGN( 1.0 , - zhgnew(ji) + epsi13 ) ) 
     692            ! salinity dif due to snow-ice formation 
     693            dsm_i_si_1d(ji) = ( zsm_snowice - sm_i_b(ji) ) * dh_snowice(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch      
     694            ! salinity dif due to bottom growth  
     695            IF (  fc_bo_i(ji) + fbif_1d(ji) + qlbbq_1d(ji)  < 0._wp ) THEN 
     696               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( zhgnew(ji), epsi13 ) * i_ice_switch 
     697            ENDIF 
     698         ENDIF 
    669699 
    670700         !  Actualize new snow and ice thickness. 
     
    680710         diag_sni_gr(ii,ij)  = diag_sni_gr(ii,ij) + dh_snowice(ji)*a_i_b(ji) * r1_rdtice 
    681711         ! 
     712         ! salt flux 
     713         sfx_thd_1d(ji) = sfx_thd_1d(ji) - zsm_snowice * a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
     714         !-------------------------------- 
     715         ! Update mass fluxes (clem) 
     716         !-------------------------------- 
     717         rdm_ice_1d(ji) = rdm_ice_1d(ji) + ( a_i_b(ji) * ht_i_b(ji) - zviold(ji) ) * rhoic  
     718         rdm_snw_1d(ji) = rdm_snw_1d(ji) + ( a_i_b(ji) * ht_s_b(ji) - zvsold(ji) ) * rhosn  
     719 
    682720      END DO !ji 
    683721      ! 
    684722      CALL wrk_dealloc( jpij, zh_i, zh_s, ztfs, zhsold, zqprec, zqfont_su, zqfont_bo, z_f_surf, zhgnew, zfmass_i ) 
    685       CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zsfx_melt, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
     723      CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zfdt_init, zfdt_final, zqt_i, zqt_s, zqt_dummy ) 
    686724      CALL wrk_dealloc( jpij, zinnermelt, zfbase, zdq_i ) 
    687725      CALL wrk_dealloc( jpij, jkmax, zdeltah, zqt_i_lay ) 
     726      ! 
     727      CALL wrk_dealloc( jpij, zviold, zvsold ) ! clem 
    688728      ! 
    689729   END SUBROUTINE lim_thd_dh 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r3808 r4161  
    1010   !!                 ! 04-2007 (M. Vancoppenolle) Energy conservation 
    1111   !!            4.0  ! 2011-02 (G. Madec) dynamical allocation 
     12   !!             -   ! 2012-05 (C. Rousset) add penetration solar flux 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_lim3 
     
    3435 
    3536   !!---------------------------------------------------------------------- 
    36    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     37   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    3738   !! $Id$ 
    3839   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    156157      DO ji = kideb , kiut 
    157158         ! is there snow or not 
    158          isnow(ji)= INT(  1._wp - MAX(  0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
     159         isnow(ji)= NINT(  1._wp - MAX( 0._wp , SIGN(1._wp, - ht_s_b(ji) ) )  ) 
    159160         ! surface temperature of fusion 
    160161!!gm ???  ztfs(ji) = rtt !!!???? 
    161          ztfs(ji) = isnow(ji) * rtt + (1.0-isnow(ji)) * rtt 
     162         ztfs(ji) = REAL( isnow(ji) ) * rtt + REAL( 1 - isnow(ji) ) * rtt 
    162163         ! layer thickness 
    163          zh_i(ji) = ht_i_b(ji) / nlay_i 
    164          zh_s(ji) = ht_s_b(ji) / nlay_s 
     164         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
     165         zh_s(ji) = ht_s_b(ji) / REAL( nlay_s ) 
    165166      END DO 
    166167 
     
    174175      DO layer = 1, nlay_s            ! vert. coord of the up. lim. of the layer-th snow layer 
    175176         DO ji = kideb , kiut 
    176             z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / nlay_s 
     177            z_s(ji,layer) = z_s(ji,layer-1) + ht_s_b(ji) / REAL( nlay_s ) 
    177178         END DO 
    178179      END DO 
     
    180181      DO layer = 1, nlay_i            ! vert. coord of the up. lim. of the layer-th ice layer 
    181182         DO ji = kideb , kiut 
    182             z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / nlay_i 
     183            z_i(ji,layer) = z_i(ji,layer-1) + ht_i_b(ji) / REAL( nlay_i ) 
    183184         END DO 
    184185      END DO 
     
    201202      DO ji = kideb , kiut 
    202203         ! switches 
    203          isnow(ji) = INT(  1._wp - MAX(  0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
     204         isnow(ji) = NINT(  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )  )  
    204205         ! hs > 0, isnow = 1 
    205206         zhsu (ji) = hnzst  ! threshold for the computation of i0 
    206207         zihic(ji) = MAX( 0._wp , 1._wp - ( ht_i_b(ji) / zhsu(ji) ) )      
    207208 
    208          i0(ji)    = ( 1._wp - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
     209         i0(ji)    = REAL( 1 - isnow(ji) ) * ( fr1_i0_1d(ji) + zihic(ji) * fr2_i0_1d(ji) ) 
    209210         !fr1_i0_1d = i0 for a thin ice surface 
    210211         !fr1_i0_2d = i0 for a thick ice surface 
     
    243244 
    244245      DO ji = kideb, kiut           ! ice initialization 
    245          zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * isnow(ji) + zftrice(ji) * ( 1._wp - isnow(ji) ) 
     246         zradtr_i(ji,0) = zradtr_s(ji,nlay_s) * REAL( isnow(ji) ) + zftrice(ji) * REAL( 1 - isnow(ji) ) 
    246247      END DO 
    247248 
     
    256257 
    257258      DO ji = kideb, kiut           ! Radiation transmitted below the ice 
    258          fstbif_1d(ji) = fstbif_1d(ji) + zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) 
     259         fstbif_1d(ji) = fstbif_1d(ji) + iatte_1d(ji) * zradtr_i(ji,nlay_i) * a_i_b(ji) / at_i_b(ji) ! clem modif 
    259260      END DO 
    260261 
     
    264265         ii = MOD( npb(ji) - 1 , jpi ) + 1 
    265266         ij =    ( npb(ji) - 1 ) / jpi + 1 
    266          fstroc(ii,ij,jl) = zradtr_i(ji,nlay_i) 
     267         fstroc(ii,ij,jl) = iatte_1d(ji) * zradtr_i(ji,nlay_i) ! clem modif 
    267268      END DO 
    268269      ! +++++ 
     
    376377            zkappa_s(ji,nlay_s)   = 2.0*rcdsn*ztcond_i(ji,0)/MAX(zeps, & 
    377378               (ztcond_i(ji,0)*zh_s(ji) + rcdsn*zh_i(ji))) 
    378             zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*isnow(ji) & 
    379                + zkappa_i(ji,0)*(1.0-isnow(ji)) 
     379            zkappa_i(ji,0)        = zkappa_s(ji,nlay_s)*REAL( isnow(ji) ) & 
     380               + zkappa_i(ji,0)*REAL( 1 - isnow(ji) ) 
    380381         END DO 
    381382         ! 
     
    658659               t_s_b(ji,nlay_s)     =  (zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) & 
    659660               *  t_i_b(ji,1))/zdiagbis(ji,nlay_s+1) & 
    660                *        MAX(0.0,SIGN(1.0,ht_s_b(ji)-zeps))  
     661               *        MAX(0.0,SIGN(1.0,ht_s_b(ji)))  
    661662 
    662663            ! surface temperature 
    663             isnow(ji)     = INT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  ) 
     664            isnow(ji)     = NINT(  1.0 - MAX( 0.0 , SIGN( 1.0 , -ht_s_b(ji) )  )  ) 
    664665            ztsuoldit(ji) = t_su_b(ji) 
    665             IF( t_su_b(ji) < ztfs(ji) )   & 
    666                t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( isnow(ji)*t_s_b(ji,1)   & 
    667                &          + (1.0-isnow(ji))*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
     666            IF( t_su_b(ji) < ztfs(ji) ) & 
     667               t_su_b(ji) = ( zindtbis(ji,numeqmin(ji)) - ztrid(ji,numeqmin(ji),3)* ( REAL( isnow(ji) )*t_s_b(ji,1)   & 
     668               &          + REAL( 1 - isnow(ji) )*t_i_b(ji,1) ) ) / zdiagbis(ji,numeqmin(ji))   
    668669         END DO 
    669670         ! 
     
    721722#endif 
    722723         !                                ! surface ice conduction flux 
    723          isnow(ji)       = INT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  ) 
    724          fc_su(ji)       =  -           isnow(ji)  * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   & 
    725             &               - ( 1._wp - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji)) 
     724         isnow(ji)       = NINT(  1._wp - MAX( 0._wp, SIGN( 1._wp, -ht_s_b(ji) ) )  ) 
     725         fc_su(ji)       =  -     REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * (t_s_b(ji,1) - t_su_b(ji))   & 
     726            &               - REAL( 1 - isnow(ji) ) * zkappa_i(ji,0) * zg1  * (t_i_b(ji,1) - t_su_b(ji)) 
    726727         !                                ! bottom ice conduction flux 
    727728         fc_bo_i(ji)     =  - zkappa_i(ji,nlay_i) * ( zg1*(t_bo_b(ji) - t_i_b(ji,nlay_i)) ) 
     
    734735         DO ji = kideb, kiut 
    735736            ! Upper snow value 
    736             fc_s(ji,0) = - isnow(ji) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )  
     737            fc_s(ji,0) = - REAL( isnow(ji) ) * zkappa_s(ji,0) * zg1s * ( t_s_b(ji,1) - t_su_b(ji) )  
    737738            ! Bott. snow value 
    738             fc_s(ji,1) = - isnow(ji)* zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )  
     739            fc_s(ji,1) = - REAL( isnow(ji) ) * zkappa_s(ji,1) * ( t_i_b(ji,1) - t_s_b(ji,1) )  
    739740         END DO 
    740741         DO ji = kideb, kiut         ! Upper ice layer 
    741             fc_i(ji,0) = - isnow(ji) * &  ! interface flux if there is snow 
     742            fc_i(ji,0) = - REAL( isnow(ji) ) * &  ! interface flux if there is snow 
    742743               ( zkappa_i(ji,0)  * ( t_i_b(ji,1) - t_s_b(ji,nlay_s ) ) ) & 
    743                - ( 1.0 - isnow(ji) ) * ( zkappa_i(ji,0) * &  
     744               - REAL( 1 - isnow(ji) ) * ( zkappa_i(ji,0) * &  
    744745               zg1 * ( t_i_b(ji,1) - t_su_b(ji) ) ) ! upper flux if not 
    745746         END DO 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r3625 r4161  
    4444 
    4545   !!---------------------------------------------------------------------- 
    46    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     46   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4747   !! $Id$ 
    4848   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7575 
    7676      INTEGER ::   ji,jk   !  dummy loop indices 
    77       INTEGER ::   zji, zjj       ,   &  !  dummy indices 
     77      INTEGER ::   ii, ij       ,   &  !  dummy indices 
    7878         ntop0          ,   &  !  old layer top index 
    7979         nbot1          ,   &  !  new layer bottom index 
     
    145145 
    146146      DO ji = kideb, kiut 
    147          zh_i(ji) = old_ht_i_b(ji) / nlay_i  
    148          zh_s(ji) = old_ht_s_b(ji) / nlay_s 
     147         zh_i(ji) = old_ht_i_b(ji) / REAL( nlay_i )  
     148         zh_s(ji) = old_ht_s_b(ji) / REAL( nlay_s ) 
    149149      END DO 
    150150 
     
    166166      DO jk = 1, nlays0 
    167167         DO ji = kideb, kiut 
    168             snind(ji)  = jk        *      INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20))) & 
    169                + snind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)-epsi20)))) 
     168            snind(ji)  = jk        *      NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji)))) & 
     169               + snind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_s_tot(ji)-zdeltah(ji))))) 
    170170            zdeltah(ji)= zdeltah(ji) + zh_s(ji) 
    171171         END DO ! ji 
     
    175175      !              0 if not 
    176176      DO ji = kideb, kiut 
    177          snswi(ji)     = MAX(0,INT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 
     177         snswi(ji)     = MAX(0,NINT(-dh_s_tot(ji)/MAX(epsi20,ABS(dh_s_tot(ji))))) 
    178178      END DO ! ji 
    179179 
     
    190190      DO jk = 1, nlayi0 
    191191         DO ji = kideb, kiut 
    192             icsuind(ji) = jk          *      INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20))) & 
    193                + icsuind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)-epsi20)))) 
     192            icsuind(ji) = jk          *      NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji)))) & 
     193               + icsuind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_surf(ji)-zdeltah(ji))))) 
    194194            zdeltah(ji) = zdeltah(ji) + zh_i(ji) 
    195195         END DO ! ji 
     
    200200      !     0 if not 
    201201      DO ji = kideb, kiut 
    202          icsuswi(ji)  = MAX(0,INT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 
     202         icsuswi(ji)  = MAX(0,NINT(-dh_i_surf(ji)/MAX(epsi20 , ABS(dh_i_surf(ji)) ) ) ) 
    203203      ENDDO 
    204204 
     
    216216      DO jk = nlayi0, 1, -1 
    217217         DO ji = kideb, kiut 
    218             icboind(ji) = (nlayi0+1-jk) *      INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))) & 
    219                &          + icboind(ji) * (1 - INT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)-epsi20))))  
     218            icboind(ji) = (nlayi0+1-jk) *      NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))) & 
     219               &          + icboind(ji) * (1 - NINT(MAX(0.0,SIGN(1.0,-dh_i_bott(ji)-zdeltah(ji)))))  
    220220            zdeltah(ji) = zdeltah(ji) + zh_i(ji) 
    221221         END DO 
     
    232232      !     0 if ablation is on the way 
    233233      DO ji = kideb, kiut  
    234          icboswi(ji) = MAX(0,INT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 
     234         icboswi(ji) = MAX(0,NINT(dh_i_bott(ji) / MAX(epsi20,ABS(dh_i_bott(ji))))) 
    235235      END DO 
    236236 
     
    248248         DO ji = kideb, kiut 
    249249            snicind(ji) = (nlays0+1-jk) & 
    250                *      INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20))) + snicind(ji)   & 
    251                * (1 - INT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)-epsi20)))) 
     250               *      NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji)))) + snicind(ji)   & 
     251               * (1 - NINT(MAX(0.0,SIGN(1.0,dh_snowice(ji)-zdeltah(ji))))) 
    252252            zdeltah(ji) = zdeltah(ji) + zh_s(ji) 
    253253         END DO 
     
    258258      !     0 if not 
    259259      DO ji = kideb, kiut 
    260          snicswi(ji)   = MAX(0,INT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 
     260         snicswi(ji)   = MAX(0,NINT(dh_snowice(ji)/MAX(epsi20,ABS(dh_snowice(ji))))) 
    261261      ENDDO 
    262262 
     
    279279 
    280280      DO ji = kideb, kiut 
    281          nbot0(ji) =  nlays0  + 1 - snind(ji) + ( 1. - snicind(ji) ) * snicswi(ji) 
     281         nbot0(ji) =  nlays0  + 1 - snind(ji) + ( 1 - snicind(ji) ) * snicswi(ji) 
    282282         ! cotes of the top of the layers 
    283283         zm0(ji,0) =  0._wp 
     
    291291            limsum = ( 1 - snswi(ji) ) * ( jk - 1 ) + snswi(ji) * ( jk + snind(ji) - 1 ) 
    292292            limsum = MIN( limsum , nlay_s ) 
    293             zm0(ji,jk) =  dh_s_tot(ji) + zh_s(ji) * limsum 
    294          END DO 
    295       END DO 
    296  
    297       DO ji = kideb, kiut 
    298          zm0(ji,nbot0(ji)) =  dh_s_tot(ji) - snicswi(ji) * dh_snowice(ji) + zh_s(ji) * nlays0 
    299          zm0(ji,1)         =  dh_s_tot(ji) * (1 -snswi(ji) ) + snswi(ji) * zm0(ji,1) 
     293            zm0(ji,jk) =  dh_s_tot(ji) + zh_s(ji) * REAL( limsum ) 
     294         END DO 
     295      END DO 
     296 
     297      DO ji = kideb, kiut 
     298         zm0(ji,nbot0(ji)) =  dh_s_tot(ji) - REAL( snicswi(ji) ) * dh_snowice(ji) + zh_s(ji) * REAL( nlays0 ) 
     299         zm0(ji,1)         =  dh_s_tot(ji) * REAL( 1 - snswi(ji) ) + REAL( snswi(ji) ) * zm0(ji,1) 
    300300      END DO 
    301301 
     
    309309 
    310310      DO ji = kideb, kiut         ! layer heat content 
    311          qm0    (ji,1) =  rhosn * (  cpic * ( rtt - ( 1. - snswi(ji) ) * tatm_ice_1d(ji)        & 
    312             &                                            - snswi(ji)  * t_s_b      (ji,1)  )   & 
     311         qm0    (ji,1) =  rhosn * (  cpic * ( rtt - REAL( 1 - snswi(ji) ) * tatm_ice_1d(ji)        & 
     312            &                                         - REAL( snswi(ji) ) * t_s_b      (ji,1)  )   & 
    313313            &                      + lfus  ) * zthick0(ji,1) 
    314314         zqts_in(ji)   =  zqts_in(ji) + qm0(ji,1)  
     
    320320            limsum      = MIN( limsum , nlay_s ) 
    321321            qm0(ji,jk)  = rhosn * ( cpic * ( rtt - t_s_b(ji,limsum) ) + lfus ) * zthick0(ji,jk) 
    322             zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, epsi20 - ht_s_b(ji) ) ) 
    323             zqts_in(ji) = zqts_in(ji) + ( 1. - snswi(ji) ) * qm0(ji,jk) * zswitch 
     322            zswitch = 1.0 - MAX (0.0, SIGN ( 1.0, - ht_s_b(ji) ) ) 
     323            zqts_in(ji) = zqts_in(ji) + REAL( 1 - snswi(ji) ) * qm0(ji,jk) * zswitch 
    324324         END DO ! jk 
    325325      END DO ! ji 
     
    360360      !------------------- 
    361361      DO ji = kideb, kiut 
    362          zh_s(ji)  = ht_s_b(ji) / nlay_s 
     362         zh_s(ji)  = ht_s_b(ji) / REAL( nlay_s ) 
    363363         z_s(ji,0) =  0._wp 
    364364      ENDDO 
     
    366366      DO jk = 1, nlay_s 
    367367         DO ji = kideb, kiut 
    368             z_s(ji,jk) =  zh_s(ji) * jk 
     368            z_s(ji,jk) =  zh_s(ji) * REAL( jk ) 
    369369         END DO 
    370370      END DO 
     
    394394                  &                 - MAX(zm0(ji,layer0-1), z_s(ji,layer1-1))) / MAX(zhl0(ji,layer0),epsi10))  
    395395               q_s_b(ji,layer1) = q_s_b(ji,layer1) + zrl01(layer1,layer0)*qm0(ji,layer0)   & 
    396                   &                                * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+epsi20)) 
     396                  &                                * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 
    397397            END DO 
    398398         END DO 
     
    410410         DO ji = kideb, kiut 
    411411            IF ( ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    412                zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    413                zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    414                WRITE(numout,*) ' violation of heat conservation : ',             & 
    415                   ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 
    416                WRITE(numout,*) ' ji, jj   : ', zji, zjj 
     412               ii                 = MOD( npb(ji) - 1, jpi ) + 1 
     413               ij                 = ( npb(ji) - 1 ) / jpi + 1 
     414               WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqts_in(ji) - zqts_fin(ji) ) * r1_rdtice 
     415               WRITE(numout,*) ' ji, jj   : ', ii, ij 
    417416               WRITE(numout,*) ' ht_s_b   : ', ht_s_b(ji) 
    418417               WRITE(numout,*) ' zqts_in  : ', zqts_in (ji) * r1_rdtice 
     
    441440      DO jk = 1, nlay_s 
    442441         DO ji = kideb, kiut 
    443             zswitch = MAX ( 0.0 , SIGN ( 1.0, epsi20 - ht_s_b(ji) ) ) 
     442            zswitch = MAX ( 0.0 , SIGN ( 1.0, - ht_s_b(ji) ) ) 
    444443            t_s_b(ji,jk) = rtt + ( 1.0 - zswitch ) * ( - zfac1 * q_s_b(ji,jk) + zfac2 ) 
    445444         END DO 
     
    480479            limsum    =  ( (icsuswi(ji)*(icsuind(ji)+jk-1) + &  
    481480               (1-icsuswi(ji))*jk))*(1-snicswi(ji)) + (jk-1)*snicswi(ji) 
    482             zm0(ji,jk)=  icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) & 
    483                +  limsum * zh_i(ji) 
    484          END DO 
    485       END DO 
    486  
    487       DO ji = kideb, kiut 
    488          zm0(ji,nbot0(ji)) =  icsuswi(ji)*dh_i_surf(ji) + snicswi(ji)*dh_snowice(ji) + dh_i_bott(ji) & 
    489             +  zh_i(ji) * nlayi0 
    490          zm0(ji,1)         =  snicswi(ji)*dh_snowice(ji) + (1-snicswi(ji))*zm0(ji,1) 
     481            zm0(ji,jk)=  REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) & 
     482               +  REAL(limsum) * zh_i(ji) 
     483         END DO 
     484      END DO 
     485 
     486      DO ji = kideb, kiut 
     487         zm0(ji,nbot0(ji)) =  REAL(icsuswi(ji))*dh_i_surf(ji) + REAL(snicswi(ji))*dh_snowice(ji) + dh_i_bott(ji) & 
     488            +  zh_i(ji) * REAL(nlayi0) 
     489         zm0(ji,1)         =  REAL(snicswi(ji))*dh_snowice(ji) + REAL(1-snicswi(ji))*zm0(ji,1) 
    491490      END DO 
    492491 
     
    521520      !---------------------------- 
    522521      DO ji = kideb, kiut         
    523          ztmelts    = ( 1.0 - icboswi(ji) ) * (-tmut * s_i_b  (ji,nlayi0) )   &   ! case of melting ice 
    524             &       +         icboswi(ji)  * (-tmut * s_i_new(ji)        )   &   ! case of forming ice 
     522         ztmelts    = REAL( 1 - icboswi(ji) ) * (-tmut * s_i_b  (ji,nlayi0) )   &   ! case of melting ice 
     523            &       +     REAL( icboswi(ji) ) * (-tmut * s_i_new(ji)        )   &   ! case of forming ice 
    525524            &       + rtt                                                         ! in Kelvin 
    526525 
     
    528527         ztform = t_i_b(ji,nlay_i) 
    529528         IF(  num_sal == 2  )   ztform = t_bo_b(ji) 
    530          qm0(ji,nbot0(ji)) = ( 1.0 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice 
    531             &              + icboswi(ji) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice 
     529         qm0(ji,nbot0(ji)) = REAL( 1 - icboswi(ji) )*qm0(ji,nbot0(ji))             &   ! case of melting ice 
     530            &              + REAL( icboswi(ji) ) * rhoic * ( cpic*(ztmelts-ztform)       &   ! case of forming ice 
    532531            + lfus *( 1.0-(ztmelts-rtt) / MIN ( (ztform-rtt) , - epsi10 ) )      &  
    533532            - rcp*(ztmelts-rtt) ) * zthick0(ji,nbot0(ji)  ) 
     
    540539         ! energy of the flooding seawater 
    541540         zqsnic = rau0 * rcp * ( rtt - t_bo_b(ji) ) * dh_snowice(ji) * & 
    542             (rhoic - rhosn) / rhoic * snicswi(ji) ! generally positive 
     541            (rhoic - rhosn) / rhoic * REAL(snicswi(ji)) ! generally positive 
    543542         ! Heat conservation diagnostic 
    544543         qt_i_in(ji,jl) = qt_i_in(ji,jl) + zqsnic  
     
    549548         ! = enthalpy of snow + enthalpy of frozen water 
    550549         zqsnic         =  zqsnow(ji) + zqsnic 
    551          qm0(ji,1)      =  snicswi(ji) * zqsnic + ( 1 - snicswi(ji) ) * qm0(ji,1) 
     550         qm0(ji,1)      =  REAL(snicswi(ji)) * zqsnic + REAL( 1 - snicswi(ji) ) * qm0(ji,1) 
    552551 
    553552      END DO ! ji 
     
    556555         DO ji = kideb, kiut 
    557556            ! Heat conservation 
    558             zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06+epsi20) ) & 
    559                &                                   * MAX( 0.0 , SIGN( 1. , nbot0(ji) - jk + epsi20 ) ) 
     557            zqti_in(ji) = zqti_in(ji) + qm0(ji,jk) * MAX( 0.0 , SIGN(1.0,ht_i_b(ji)-epsi06) ) & 
     558               &                                   * MAX( 0.0 , SIGN( 1. , REAL(nbot0(ji) - jk) ) ) 
    560559         END DO 
    561560      END DO 
     
    575574      !------------------ 
    576575      DO ji = kideb, kiut 
    577          zh_i(ji) = ht_i_b(ji) / nlay_i 
     576         zh_i(ji) = ht_i_b(ji) / REAL( nlay_i ) 
    578577      ENDDO 
    579578 
     
    606605               q_i_b(ji,layer1) = q_i_b(ji,layer1) &  
    607606                  + zrl01(layer1,layer0)*qm0(ji,layer0) & 
    608                   * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06+epsi20)) & 
    609                   * MAX(0.0,SIGN(1.0,nbot0(ji)-layer0+epsi20)) 
     607                  * MAX(0.0,SIGN(1.0,ht_i_b(ji)-epsi06)) & 
     608                  * MAX(0.0,SIGN(1.0,REAL(nbot0(ji)-layer0))) 
    610609            END DO 
    611610         END DO 
     
    622621      END DO 
    623622      ! 
    624       DO ji = kideb, kiut 
    625          IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
    626             zji                 = MOD( npb(ji) - 1, jpi ) + 1 
    627             zjj                 = ( npb(ji) - 1 ) / jpi + 1 
    628             WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 
    629             WRITE(numout,*) ' ji, jj   : ', zji, zjj 
    630             WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji) 
    631             WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice 
    632             WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 
    633             WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 
    634             WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 
    635             WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 
    636             WRITE(numout,*) ' icsuswi  : ', icsuswi(ji) 
    637             WRITE(numout,*) ' icboswi  : ', icboswi(ji) 
    638             WRITE(numout,*) ' snicswi  : ', snicswi(ji) 
    639          ENDIF 
    640       END DO 
     623      IF ( con_i ) THEN 
     624         DO ji = kideb, kiut 
     625            IF ( ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice  >  1.0e-6 ) THEN 
     626               ii                 = MOD( npb(ji) - 1, jpi ) + 1 
     627               ij                 = ( npb(ji) - 1 ) / jpi + 1 
     628               WRITE(numout,*) ' violation of heat conservation : ', ABS ( zqti_in(ji) - zqti_fin(ji) ) * r1_rdtice 
     629               WRITE(numout,*) ' ji, jj   : ', ii, ij 
     630               WRITE(numout,*) ' ht_i_b   : ', ht_i_b(ji) 
     631               WRITE(numout,*) ' zqti_in  : ', zqti_in (ji) * r1_rdtice 
     632               WRITE(numout,*) ' zqti_fin : ', zqti_fin(ji) * r1_rdtice 
     633               WRITE(numout,*) ' dh_i_bott: ', dh_i_bott(ji) 
     634               WRITE(numout,*) ' dh_i_surf: ', dh_i_surf(ji) 
     635               WRITE(numout,*) ' dh_snowice:', dh_snowice(ji) 
     636               WRITE(numout,*) ' icsuswi  : ', icsuswi(ji) 
     637               WRITE(numout,*) ' icboswi  : ', icboswi(ji) 
     638               WRITE(numout,*) ' snicswi  : ', snicswi(ji) 
     639            ENDIF 
     640         END DO 
     641      ENDIF 
    641642 
    642643      !---------------------- 
  • branches/2013/dev_LOCEAN_2013/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r3625 r4161  
    4646 
    4747   !!---------------------------------------------------------------------- 
    48    !! NEMO/LIM3 3.4 , UCL - NEMO Consortium (2011) 
     48   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
    4949   !! $Id$ 
    5050   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
     
    7878      !!               update ht_s_b, ht_i_b and tbif_1d(:,:)       
    7979      !!------------------------------------------------------------------------ 
    80       INTEGER  ::   ji,jj,jk,jl,jm   ! dummy loop indices 
    81       INTEGER  ::   layer, nbpac     ! local integers  
    82       INTEGER  ::   zji, zjj, iter   !   -       - 
    83       REAL(wp) ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zde   ! local scalars 
     80      INTEGER ::   ji,jj,jk,jl,jm   ! dummy loop indices 
     81      INTEGER ::   layer, nbpac     ! local integers  
     82      INTEGER ::   ii, ij, iter   !   -       - 
     83      REAL(wp)  ::   ztmelts, zdv, zqold, zfrazb, zweight, zalphai, zindb, zinda, zde  ! local scalars 
    8484      REAL(wp) ::   zgamafr, zvfrx, zvgx, ztaux, ztwogp, zf , zhicol_new        !   -      - 
    8585      REAL(wp) ::   ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, zsqcd , zhicrit   !   -      - 
    86       REAL(wp) ::   zcoef                                                       !   -      - 
    8786      LOGICAL  ::   iterate_frazil   ! iterate frazil ice collection thickness 
    8887      CHARACTER (len = 15) :: fieldid 
     
    160159               DO ji = 1, jpi 
    161160                  !Energy of melting q(S,T) [J.m-3] 
    162                   e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * nlay_i 
     161                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) / MAX( area(ji,jj) * v_i(ji,jj,jl) ,  epsi10 ) * REAL( nlay_i ) 
    163162                  zindb = 1._wp - MAX(  0._wp , SIGN( 1._wp , -v_i(ji,jj,jl) )  )   !0 if no ice and 1 if yes 
    164163                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * unit_fac * zindb 
     
    342341         CASE ( 2 )                    ! Sice = F(z,t) [Vancoppenolle et al (2005)] 
    343342            DO ji = 1, nbpac 
    344                zji =   MOD( npac(ji) - 1 , jpi ) + 1 
    345                zjj =      ( npac(ji) - 1 ) / jpi + 1 
    346                zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(zji,zjj)  ) 
     343               ii =   MOD( npac(ji) - 1 , jpi ) + 1 
     344               ij =      ( npac(ji) - 1 ) / jpi + 1 
     345               zs_newice(ji) = MIN(  4.606 + 0.91 / zh_newice(ji) , s_i_max , 0.5 * sss_m(ii,ij)  ) 
    347346            END DO 
    348347         CASE ( 3 )                    ! Sice = F(z) [multiyear ice] 
     
    389388         END DO 
    390389 
    391          !--------------------------------- 
    392          ! Salt flux due to new ice growth 
    393          !--------------------------------- 
    394          ! note that for constant salinity zs_newice() = bulk_sal (see top of the subroutine) 
    395          DO ji = 1, nbpac 
    396             sfx_thd_1d(ji) = sfx_thd_1d(ji) - zs_newice(ji) * rhoic * zv_newice(ji) * r1_rdtice 
    397             rdm_ice_1d(ji) = rdm_ice_1d(ji) +                 rhoic * zv_newice(ji) 
    398          END DO ! ji 
    399  
    400390         !------------------------------------ 
    401391         ! Diags for energy conservation test 
    402392         !------------------------------------ 
    403393         DO ji = 1, nbpac 
    404             zji = MOD( npac(ji) - 1 , jpi ) + 1 
    405             zjj =    ( npac(ji) - 1 ) / jpi + 1 
     394            ii = MOD( npac(ji) - 1 , jpi ) + 1 
     395            ij =    ( npac(ji) - 1 ) / jpi + 1 
    406396            ! 
    407             zde = ze_newice(ji) / unit_fac * area(zji,zjj) * zv_newice(ji) 
     397            zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 
    408398            ! 
    409