Changeset 4649


Ignore:
Timestamp:
2014-05-27T11:28:12+02:00 (6 years ago)
Author:
clem
Message:

finalizing LIM3 heat budget conservation + multiple minor bugs corrections

Location:
branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM
Files:
26 edited

Legend:

Unmodified
Added
Removed
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/field_def.xml

    r4634 r4649  
    152152    
    153153         <!-- *_oce variables available with ln_blk_clio or ln_blk_core --> 
    154          <field id="qns_oce"      long_name="Non solar Downward Heat Flux over open ocean"                 unit="W/m2"     /> 
    155154         <field id="qlw_oce"      long_name="Longwave Downward Heat Flux over open ocean"                  unit="W/m2"     /> 
    156155         <field id="qsb_oce"      long_name="Sensible Downward Heat Flux over open ocean"                  unit="W/m2"     /> 
     
    222221         <field id="qsr_oce"      long_name="solar heat flux at ocean surface"                             unit="W/m2"     /> 
    223222         <field id="qns_oce"      long_name="non-solar heat flux at ocean surface"                         unit="W/m2"     /> 
     223         <field id="qt_ice"       long_name="total flux at ice surface"                                    unit="W/m2"     /> 
     224         <field id="qsr_ice"      long_name="solar heat flux at ice surface"                               unit="W/m2"     /> 
     225         <field id="qns_ice"      long_name="non-solar heat flux at ice surface"                           unit="W/m2"     /> 
     226         <field id="qtr_ice"      long_name="solar heat flux transmitted thru the ice"                     unit="W/m2"     /> 
    224227         <field id="utau_ice"     long_name="Wind stress along i-axis over the ice at i-point"             unit="N/m2"     /> 
    225228         <field id="vtau_ice"     long_name="Wind stress along j-axis over the ice at i-point"             unit="N/m2"     /> 
    226     <field id="qsr_io"       long_name="Ice-Oce downward solar heat flux"                             unit="W/m2"     /> 
    227     <field id="qns_io"       long_name="Ice-Oce downward non-solar heat flux"                         unit="W/m2"     /> 
    228229         <field id="micesalt"     long_name="Mean ice salinity"                                            unit="psu"      /> 
    229230         <field id="miceage"      long_name="Mean ice age"                                                 unit="years"    /> 
     
    237238 
    238239         <field id="micet"        long_name="Mean ice temperature"                                         unit="degC"     /> 
    239          <field id="icehc"        long_name="ice total heat content"                                       unit="10^9 J"   />  
     240         <field id="icehc"        long_name="ice total heat content"                                       unit="10^9J"   />  
    240241         <field id="isnowhc"      long_name="snow total heat content"                                      unit="10^9J"    /> 
    241242         <field id="icest"        long_name="ice surface temperature"                                      unit="degC"     /> 
     
    251252         <field id="icetrp"       long_name="ice volume transport"                                         unit="m/day"   /> 
    252253         <field id="snwtrp"       long_name="snw volume transport"                                         unit="m/day"   /> 
    253          <field id="deitrp"       long_name="advected ice enhalpy"                                         unit="W/2"   /> 
    254          <field id="destrp"       long_name="advected snw enhalpy"                                         unit="W/2"   /> 
     254         <field id="deitrp"       long_name="advected ice enhalpy"                                         unit="W/m2"   /> 
     255         <field id="destrp"       long_name="advected snw enhalpy"                                         unit="W/m2"   /> 
    255256 
    256257         <field id="sfxbri"       long_name="brine salt flux"                                              unit="psu*kg/m2/day" /> 
     
    274275         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   /> 
    275276         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   /> 
    276  
    277          <field id="hfxdhc1"   long_name="Heat content variation in snow and ice"   unit="W/m2" /> 
    278          <field id="hfxspr"    long_name="Heat content of snow precip"              unit="W/m2" /> 
    279          <field id="hfxqsr"    long_name="solar fluxes given to ocean"             unit="W/m2"  /> 
    280          <field id="hfxqns"    long_name="non solar fluxes given to ocean"         unit="W/m2"  /> 
    281  
    282          <field id="hfxthd"   long_name="heat fluxes from ice-ocean exchange during thermo"              unit="W/m2"  /> 
    283          <field id="hfxdyn"   long_name="heat fluxes from ice-ocean exchange during dynamic"             unit="W/m2"  /> 
    284          <field id="hfxres"   long_name="heat fluxes from ice-ocean exchange during resultant"           unit="W/m2"  /> 
    285          <field id="hfxsnw"   long_name="heat fluxes from snow-ocean exchange"                           unit="W/m2"  /> 
    286          <field id="hfxsub"   long_name="heat fluxes from sublimation"                                   unit="W/m2"  /> 
    287          <field id="hfxerr"   long_name="heat fluxes error after heat diffusion"                         unit="W/m2"  /> 
    288          <field id="hfxerr_rem" long_name="heat fluxes error after remapping"                         unit="W/m2"  /> 
    289          <field id="hfxtot"   long_name="heat fluxes total used by ice"                                  unit="W/m2"  /> 
    290          <field id="hfxout"   long_name="non solar heat fluxes received by the ocean"             unit="W/m2"  /> 
    291          <field id="hfxin"    long_name="total     heat fluxes at the ice surface"             unit="W/m2"  /> 
     277         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   /> 
     278 
     279         <field id="hfxsum"   long_name="heat fluxes causing surface ice melt"            unit="W/m2"  /> 
     280         <field id="hfxbom"   long_name="heat fluxes causing bottom ice melt"             unit="W/m2"  /> 
     281         <field id="hfxbog"   long_name="heat fluxes causing bottom ice growth"           unit="W/m2"  /> 
     282         <field id="hfxdif"   long_name="heat fluxes causing ice temperature change"      unit="W/m2"  /> 
     283         <field id="hfxopw"   long_name="heat fluxes causing open water ice formation"    unit="W/m2"  /> 
     284         <field id="hfxsnw"   long_name="heat fluxes causing snow melt"                   unit="W/m2"  /> 
     285         <field id="hfxerr"   long_name="heat fluxes error after heat diffusion"          unit="W/m2"  /> 
     286         <field id="hfxerr_rem" long_name="heat fluxes error after remapping"             unit="W/m2"  /> 
     287         <field id="hfxout"   long_name="non solar heat fluxes received by the ocean"     unit="W/m2"  /> 
     288         <field id="hfxin"    long_name="total     heat fluxes at the ice surface"        unit="W/m2"  /> 
     289 
     290    <!-- heat flux associated with mass exchange --> 
     291         <field id="hfxthd"   long_name="heat fluxes from ice-ocean mass exchange during thermo"              unit="W/m2"  /> 
     292         <field id="hfxdyn"   long_name="heat fluxes from ice-ocean mass exchange during dynamic"             unit="W/m2"  /> 
     293         <field id="hfxres"   long_name="heat fluxes from ice-ocean mass exchange during resultant"           unit="W/m2"  /> 
     294         <field id="hfxsub"   long_name="heat fluxes from ice-atm. mass exchange during sublimation"          unit="W/m2"  /> 
     295         <field id="hfxspr"   long_name="heat fluxes from ice-atm. mass exchange during snow precip"          unit="W/m2" /> 
     296 
     297    <!-- diags --> 
     298         <field id="hfxdhc"    long_name="Heat content variation in snow and ice"   unit="W/m2" /> 
     299         <field id="hfxtur"    long_name="turbulent heat flux at the ice base"      unit="W/m2"  /> 
    292300 
    293301      </field_group> 
     
    403411    <field id="ibgvfxsum"    long_name="global mean volume flux (surface melt)"      unit="m/day"   /> 
    404412    <field id="ibgvfxres"    long_name="global mean volume flux (resultant)"         unit="m/day"   /> 
     413    <field id="ibgvfxspr"    long_name="global mean volume flux (snow precip)"       unit="m/day"   /> 
     414    <field id="ibgvfxsnw"    long_name="global mean volume flux (snow melt)"         unit="m/day"   /> 
     415    <field id="ibgvfxsub"    long_name="global mean volume flux (snow sublimation)"  unit="m/day"   /> 
    405416 
    406417    <field id="ibgsfx"       long_name="global mean salt flux (total)"            unit="psu*m/day"   /> 
     
    415426 
    416427 
    417         <field id="ibghfxdhc1"   long_name="Heat content variation in snow and ice"   unit="W" /> 
     428        <field id="ibghfxdhc"    long_name="Heat content variation in snow and ice"   unit="W" /> 
    418429        <field id="ibghfxspr"    long_name="Heat content of snow precip"              unit="W" /> 
    419         <field id="ibghfxqsr"    long_name="solar fluxes given to ocean"             unit="W"  /> 
    420         <field id="ibghfxqns"    long_name="non solar fluxes given to ocean"         unit="W"  /> 
    421430 
    422431        <field id="ibghfxthd"   long_name="heat fluxes from ice-ocean exchange during thermo"              unit="W"  /> 
     432        <field id="ibghfxsum"   long_name="heat fluxes causing surface ice melt"              unit="W"  /> 
     433        <field id="ibghfxbom"   long_name="heat fluxes causing bottom ice melt"              unit="W"  /> 
     434        <field id="ibghfxbog"   long_name="heat fluxes causing bottom ice growth"              unit="W"  /> 
     435        <field id="ibghfxdif"   long_name="heat fluxes causing ice temperature change"              unit="W"  /> 
     436        <field id="ibghfxopw"   long_name="heat fluxes causing open water ice formation"              unit="W"  /> 
    423437        <field id="ibghfxdyn"   long_name="heat fluxes from ice-ocean exchange during dynamic"             unit="W"  /> 
    424438        <field id="ibghfxres"   long_name="heat fluxes from ice-ocean exchange during resultant"           unit="W"  /> 
    425439        <field id="ibghfxsub"   long_name="heat fluxes from sublimation"                                   unit="W"  /> 
    426440        <field id="ibghfxsnw"   long_name="heat fluxes from snow-ocean exchange"                           unit="W"  /> 
    427         <field id="ibghfxtot"   long_name="heat fluxes total used by ice"                                  unit="W"  /> 
    428441        <field id="ibghfxout"   long_name="non solar heat fluxes received by the ocean"                    unit="W"  /> 
    429442        <field id="ibghfxin"    long_name="total heat fluxes at the ice surface"                           unit="W"  /> 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml

    r4634 r4649  
    5353   <file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc --> 
    5454     <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     /> 
    55      <field field_ref="qsr"          name="hfxqsr"   long_name="surface_net_downward_shortwave_flux"           /> 
    56      <field field_ref="qt"           name="hfxtot"   long_name="surface_net_downward_total_heat_flux"          /> 
     55     <field field_ref="qsr_oce"      name="qsr_oce"  long_name="downward shortwave flux at ocean surface"           /> 
     56     <field field_ref="qns_oce"      name="qns_oce"  long_name="downward non solar flux at ocean surface"           /> 
     57     <field field_ref="qt_oce"       name="qt_oce"   long_name="downward total flux at ocean surface"           /> 
     58     <field field_ref="qsr_ice"      name="qsr_ice"  long_name="downward shortwave flux at ice surface"           /> 
     59     <field field_ref="qns_ice"      name="qns_ice"  long_name="downward non solar flux at ice surface"           /> 
     60     <field field_ref="qtr_ice"      name="qtr_ice"  long_name="shortwave flux transmitted thru the ice"           /> 
     61     <field field_ref="qt_ice"       name="qt_ice"   long_name="downward total flux at ice surface"           /> 
    5762     <field field_ref="saltflx"      name="sfx"  /> 
    5863     <field field_ref="taum"         name="taum" /> 
     
    6368          <field field_ref="utau_ice"         name="utau_ice" /> 
    6469          <field field_ref="vtau_ice"         name="vtau_ice" /> 
    65      <!-- clem 
    66           <field field_ref="qsr_io"           name="iicesflx" /> 
    67           <field field_ref="qns_io"           name="iicenflx" /> 
    68      --> 
    6970 
    7071   </file> 
     
    9798 
    9899   <file id="file6" name_suffix="_icemod" description="ice variables" enabled=".true." > 
    99      <field field_ref="snowthic_cea"     name="sndept"     long_name="surface_snow_thickness"   /> 
     100     <field field_ref="snowthic_cea"     name="snthic"     long_name="surface_snow_thickness"   /> 
    100101     <field field_ref="icethic_cea"      name="sithic"     long_name="sea_ice_thickness"        /> 
    101102          <field field_ref="icevolu"          name="sivolu" /> 
     
    113114          <field field_ref="vfxsnw"          name="vfxsnw" /> 
    114115          <field field_ref="vfxsub"          name="vfxsub" /> 
     116          <field field_ref="vfxspr"          name="vfxspr" /> 
    115117 
    116118          <field field_ref="icetrp"          name="sivtrp" /> 
     
    129131          <field field_ref="sfx"              name="sfx" /> 
    130132 
    131           <field field_ref="hfxdhc1"          name="hfxdhc1"    /> 
    132           <field field_ref="hfxspr"           name="hfxspr"    /> 
    133           <field field_ref="hfxqsr"           name="hfxqsr"    /> 
    134           <field field_ref="hfxqns"           name="hfxqns"    /> 
    135  
    136           <field field_ref="hfxthd"          name="hfxthd"    /> 
     133          <field field_ref="hfxsum"          name="hfxsum"    /> 
     134          <field field_ref="hfxbom"          name="hfxbom"    /> 
     135          <field field_ref="hfxbog"          name="hfxbog"    /> 
     136          <field field_ref="hfxdif"          name="hfxdif"    /> 
     137          <field field_ref="hfxopw"          name="hfxopw"    /> 
     138          <field field_ref="hfxout"          name="hfxout"    /> 
     139          <field field_ref="hfxin"           name="hfxin"    /> 
     140          <field field_ref="hfxsnw"          name="hfxsnw"    /> 
     141          <field field_ref="hfxerr"          name="hfxerr"    /> 
     142          <field field_ref="hfxerr_rem"      name="hfxerr_rem"    /> 
     143 
     144     <!-- ice-ocean heat flux from mass exchange --> 
    137145          <field field_ref="hfxdyn"          name="hfxdyn"    /> 
    138146          <field field_ref="hfxres"          name="hfxres"    /> 
    139           <field field_ref="hfxout"          name="hfxout"    /> 
    140           <field field_ref="hfxin"           name="hfxin"    /> 
    141           <field field_ref="hfxtot"          name="hfxtot"    /> 
    142           <field field_ref="hfxsnw"          name="hfxsnw"    /> 
     147          <field field_ref="hfxthd"          name="hfxthd"    /> 
     148     <!-- ice-atm. heat flux from mass exchange --> 
    143149          <field field_ref="hfxsub"          name="hfxsub"    /> 
    144           <field field_ref="hfxerr"          name="hfxerr"    /> 
    145           <field field_ref="hfxerr_rem"      name="hfxerr_rem"    /> 
     150          <field field_ref="hfxspr"          name="hfxspr"    /> 
     151 
     152     <!-- diags --> 
     153          <field field_ref="hfxdhc"          name="hfxdhc"    /> 
     154          <field field_ref="hfxtur"          name="hfxtur"    /> 
    146155 
    147156          <field field_ref="isst"             name="sst" /> 
     
    167176          <field field_ref="iceconc_cat"      name="siconcat"/> 
    168177          <field field_ref="icethic_cat"      name="sithicat"/> 
    169           <field field_ref="snowthic_cat"     name="sndeptcat"/> 
     178          <field field_ref="snowthic_cat"     name="snthicat"/> 
    170179          <field field_ref="salinity_cat"     name="salincat"/> 
    171           <field field_ref="brinevol_cat"     name="sibrincat"/> 
     180          <field field_ref="brinevol_cat"     name="sibricat"/> 
    172181 
    173182   </file> 
     
    191200          <field field_ref="bgfrctem"     name="bgfrctem"    /> 
    192201          <field field_ref="bgfrcsal"     name="bgfrcsal"    /> 
    193      <!-- 
    194           <field field_ref="bgmistem"     name="bgmistem"    /> 
    195           <field field_ref="bgmissal"     name="bgmissal"    /> 
    196      --> 
     202 
    197203        </file> 
    198204 
     
    215221          <field field_ref="ibgvfxsum"    name="ibgvfxsum"      /> 
    216222          <field field_ref="ibgvfxres"    name="ibgvfxres"      /> 
     223          <field field_ref="ibgvfxspr"    name="ibgvfxspr"      /> 
     224          <field field_ref="ibgvfxsnw"    name="ibgvfxsnw"      /> 
     225          <field field_ref="ibgvfxsub"    name="ibgvfxsub"      /> 
    217226 
    218227          <field field_ref="ibgsfx"       name="ibgsfx"     /> 
     
    226235          <field field_ref="ibgsfxsum"    name="ibgsfxsum"      /> 
    227236 
    228           <field field_ref="ibghfxdhc1"    name="ibghfxdhc1"    /> 
     237          <field field_ref="ibghfxdhc"    name="ibghfxdhc"    /> 
    229238          <field field_ref="ibghfxspr"    name="ibghfxspr"    /> 
    230           <field field_ref="ibghfxqsr"    name="ibghfxqsr"    /> 
    231           <field field_ref="ibghfxqns"    name="ibghfxqns"    /> 
    232239 
    233240          <field field_ref="ibghfxres"    name="ibghfxres"    /> 
     
    235242          <field field_ref="ibghfxdyn"    name="ibghfxdyn"    /> 
    236243          <field field_ref="ibghfxthd"    name="ibghfxthd"    /> 
     244          <field field_ref="ibghfxsum"    name="ibghfxsum"    /> 
     245          <field field_ref="ibghfxbom"    name="ibghfxbom"    /> 
     246          <field field_ref="ibghfxbog"    name="ibghfxbog"    /> 
     247          <field field_ref="ibghfxdif"    name="ibghfxdif"    /> 
     248          <field field_ref="ibghfxopw"    name="ibghfxopw"    /> 
    237249          <field field_ref="ibghfxout"    name="ibghfxout"    /> 
    238250          <field field_ref="ibghfxin"    name="ibghfxin"    /> 
    239           <field field_ref="ibghfxtot"    name="ibghfxtot"    /> 
    240251          <field field_ref="ibghfxsnw"    name="ibghfxsnw"    /> 
    241252 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef_default.xml

    r4634 r4649  
    5353   <file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc --> 
    5454     <field field_ref="empmr"        name="wfo"      long_name="water_flux_into_sea_water"                     /> 
    55      <field field_ref="qsr"          name="hfxqsr"   long_name="surface_net_downward_shortwave_flux"           /> 
    56      <field field_ref="qt"           name="hfxtot"   long_name="surface_net_downward_total_heat_flux"          /> 
     55     <field field_ref="qsr_oce"      name="qsr_oce"  long_name="downward shortwave flux at ocean surface"           /> 
     56     <field field_ref="qns_oce"      name="qns_oce"  long_name="downward non solar flux at ocean surface"           /> 
     57     <field field_ref="qt_oce"       name="qt_oce"   long_name="downward total flux at ocean surface"           /> 
     58     <field field_ref="qsr_ice"      name="qsr_ice"  long_name="downward shortwave flux at ice surface"           /> 
     59     <field field_ref="qns_ice"      name="qns_ice"  long_name="downward non solar flux at ice surface"           /> 
     60     <field field_ref="qtr_ice"      name="qtr_ice"  long_name="shortwave flux transmitted thru the ice"           /> 
     61     <field field_ref="qt_ice"       name="qt_ice"   long_name="downward total flux at ice surface"           /> 
    5762     <field field_ref="saltflx"      name="sfx"  /> 
    5863     <field field_ref="taum"         name="taum" /> 
     
    6368          <field field_ref="utau_ice"         name="utau_ice" /> 
    6469          <field field_ref="vtau_ice"         name="vtau_ice" /> 
    65      <!-- clem 
    66           <field field_ref="qsr_io"           name="iicesflx" /> 
    67           <field field_ref="qns_io"           name="iicenflx" /> 
    68      --> 
    6970 
    7071   </file> 
     
    9798 
    9899   <file id="file6" name_suffix="_icemod" description="ice variables" enabled=".true." > 
    99      <field field_ref="snowthic_cea"     name="sndept"     long_name="surface_snow_thickness"   /> 
     100     <field field_ref="snowthic_cea"     name="snthic"     long_name="surface_snow_thickness"   /> 
    100101     <field field_ref="icethic_cea"      name="sithic"     long_name="sea_ice_thickness"        /> 
    101102          <field field_ref="icevolu"          name="sivolu" /> 
     
    113114          <field field_ref="vfxsnw"          name="vfxsnw" /> 
    114115          <field field_ref="vfxsub"          name="vfxsub" /> 
     116          <field field_ref="vfxspr"          name="vfxspr" /> 
    115117 
    116118          <field field_ref="icetrp"          name="sivtrp" /> 
     
    129131          <field field_ref="sfx"              name="sfx" /> 
    130132 
    131           <field field_ref="hfxdhc1"          name="hfxdhc1"    /> 
    132           <field field_ref="hfxspr"           name="hfxspr"    /> 
    133           <field field_ref="hfxqsr"           name="hfxqsr"    /> 
    134           <field field_ref="hfxqns"           name="hfxqns"    /> 
    135  
    136           <field field_ref="hfxthd"          name="hfxthd"    /> 
     133          <field field_ref="hfxsum"          name="hfxsum"    /> 
     134          <field field_ref="hfxbom"          name="hfxbom"    /> 
     135          <field field_ref="hfxbog"          name="hfxbog"    /> 
     136          <field field_ref="hfxdif"          name="hfxdif"    /> 
     137          <field field_ref="hfxopw"          name="hfxopw"    /> 
     138          <field field_ref="hfxout"          name="hfxout"    /> 
     139          <field field_ref="hfxin"           name="hfxin"    /> 
     140          <field field_ref="hfxsnw"          name="hfxsnw"    /> 
     141          <field field_ref="hfxerr"          name="hfxerr"    /> 
     142          <field field_ref="hfxerr_rem"      name="hfxerr_rem"    /> 
     143 
     144     <!-- ice-ocean heat flux from mass exchange --> 
    137145          <field field_ref="hfxdyn"          name="hfxdyn"    /> 
    138146          <field field_ref="hfxres"          name="hfxres"    /> 
    139           <field field_ref="hfxout"          name="hfxout"    /> 
    140           <field field_ref="hfxin"           name="hfxin"    /> 
    141           <field field_ref="hfxtot"          name="hfxtot"    /> 
    142           <field field_ref="hfxsnw"          name="hfxsnw"    /> 
     147          <field field_ref="hfxthd"          name="hfxthd"    /> 
     148     <!-- ice-atm. heat flux from mass exchange --> 
    143149          <field field_ref="hfxsub"          name="hfxsub"    /> 
    144           <field field_ref="hfxerr"          name="hfxerr"    /> 
    145           <field field_ref="hfxerr_rem"      name="hfxerr_rem"    /> 
     150          <field field_ref="hfxspr"          name="hfxspr"    /> 
     151 
     152     <!-- diags --> 
     153          <field field_ref="hfxdhc"          name="hfxdhc"    /> 
     154          <field field_ref="hfxtur"          name="hfxtur"    /> 
    146155 
    147156          <field field_ref="isst"             name="sst" /> 
     
    167176          <field field_ref="iceconc_cat"      name="siconcat"/> 
    168177          <field field_ref="icethic_cat"      name="sithicat"/> 
    169           <field field_ref="snowthic_cat"     name="sndeptcat"/> 
     178          <field field_ref="snowthic_cat"     name="snthicat"/> 
    170179          <field field_ref="salinity_cat"     name="salincat"/> 
    171           <field field_ref="brinevol_cat"     name="sibrincat"/> 
     180          <field field_ref="brinevol_cat"     name="sibricat"/> 
    172181 
    173182   </file> 
     
    191200          <field field_ref="bgfrctem"     name="bgfrctem"    /> 
    192201          <field field_ref="bgfrcsal"     name="bgfrcsal"    /> 
    193      <!-- 
    194           <field field_ref="bgmistem"     name="bgmistem"    /> 
    195           <field field_ref="bgmissal"     name="bgmissal"    /> 
    196      --> 
     202 
    197203        </file> 
    198204 
     
    215221          <field field_ref="ibgvfxsum"    name="ibgvfxsum"      /> 
    216222          <field field_ref="ibgvfxres"    name="ibgvfxres"      /> 
     223          <field field_ref="ibgvfxspr"    name="ibgvfxspr"      /> 
     224          <field field_ref="ibgvfxsnw"    name="ibgvfxsnw"      /> 
     225          <field field_ref="ibgvfxsub"    name="ibgvfxsub"      /> 
    217226 
    218227          <field field_ref="ibgsfx"       name="ibgsfx"     /> 
     
    226235          <field field_ref="ibgsfxsum"    name="ibgsfxsum"      /> 
    227236 
    228           <field field_ref="ibghfxdhc1"    name="ibghfxdhc1"    /> 
     237          <field field_ref="ibghfxdhc"    name="ibghfxdhc"    /> 
    229238          <field field_ref="ibghfxspr"    name="ibghfxspr"    /> 
    230           <field field_ref="ibghfxqsr"    name="ibghfxqsr"    /> 
    231           <field field_ref="ibghfxqns"    name="ibghfxqns"    /> 
    232239 
    233240          <field field_ref="ibghfxres"    name="ibghfxres"    /> 
     
    235242          <field field_ref="ibghfxdyn"    name="ibghfxdyn"    /> 
    236243          <field field_ref="ibghfxthd"    name="ibghfxthd"    /> 
     244          <field field_ref="ibghfxsum"    name="ibghfxsum"    /> 
     245          <field field_ref="ibghfxbom"    name="ibghfxbom"    /> 
     246          <field field_ref="ibghfxbog"    name="ibghfxbog"    /> 
     247          <field field_ref="ibghfxdif"    name="ibghfxdif"    /> 
     248          <field field_ref="ibghfxopw"    name="ibghfxopw"    /> 
    237249          <field field_ref="ibghfxout"    name="ibghfxout"    /> 
    238250          <field field_ref="ibghfxin"    name="ibghfxin"    /> 
    239           <field field_ref="ibghfxtot"    name="ibghfxtot"    /> 
    240251          <field field_ref="ibghfxsnw"    name="ibghfxsnw"    /> 
    241252 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice

    r4634 r4649  
    1515!----------------------------------------------------------------------- 
    1616   cn_icerst_in  = "restart_ice"   !  suffix of ice restart name (input) 
    17    cn_icerst_out = "restart_ice"      !  suffix of ice restart name (output) 
    18    ln_limdyn   = .true.    !  ice dynamics (T) or thermodynamics only (F) 
    19    amax        = 0.999      !  maximum ice concentration 
    20    cai         =  1.40e-3  !  atmospheric drag over sea ice 
    21    cao         =  1.00e-3  !  atmospheric drag over ocean 
    22    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  = .true.    !  output the heat and salt budgets (T) or not (F) 
     17   cn_icerst_out = "restart_ice"   !  suffix of ice restart name (output) 
     18   ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F) 
     19   amax          = 0.999           !  maximum ice concentration 
     20   cai           = 1.40e-3         !  atmospheric drag over sea ice (clio) 
     21   cao           = 1.00e-3         !  atmospheric drag over ocean   (clio) 
     22   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  = .true.          !  output the heat and salt budgets (T) or not (F) 
    2525/ 
    2626!----------------------------------------------------------------------- 
     
    147147!&namicehsb       !  Heat and salt budgets  
    148148!!----------------------------------------------------------------------- 
    149 ! 
    150149!/ 
    151 !----------------------------------------------------------------------- 
    152 &namiceout     !   parameters for outputs 
    153 !----------------------------------------------------------------------- 
    154 !   noumef      =   43      !  number of fields 
    155 !SF   add_diag_swi=    1      !  1 -> diagnose distribution in thickness space 
    156 !SF                           !  0 -> only simple diagnostics 
    157 ! 
    158 !           !         title of the field           !  name     !   units   !  save  ! multipl. ! additive ! 
    159 !           !                                      !           !           ! or not !  factor  !  factor  ! 
    160 /  
    161150 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_lim3

    r4634 r4649  
    1515!----------------------------------------------------------------------- 
    1616   cn_icerst_in  = "restart_ice"   !  suffix of ice restart name (input) 
    17    cn_icerst_out = "restart_ice"      !  suffix of ice restart name (output) 
    18    ln_limdyn   = .true.    !  ice dynamics (T) or thermodynamics only (F) 
    19    amax        = 0.999      !  maximum ice concentration 
    20    cai         =  1.40e-3  !  atmospheric drag over sea ice 
    21    cao         =  1.00e-3  !  atmospheric drag over ocean 
    22    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  = .true.    !  output the heat and salt budgets (T) or not (F) 
     17   cn_icerst_out = "restart_ice"   !  suffix of ice restart name (output) 
     18   ln_limdyn     = .true.          !  ice dynamics (T) or thermodynamics only (F) 
     19   amax          = 0.999           !  maximum ice concentration 
     20   cai           = 1.40e-3         !  atmospheric drag over sea ice (clio) 
     21   cao           = 1.00e-3         !  atmospheric drag over ocean   (clio) 
     22   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  = .true.          !  output the heat and salt budgets (T) or not (F) 
    2525/ 
    2626!----------------------------------------------------------------------- 
     
    147147!&namicehsb       !  Heat and salt budgets  
    148148!!----------------------------------------------------------------------- 
    149 ! 
    150149!/ 
    151 !----------------------------------------------------------------------- 
    152 &namiceout     !   parameters for outputs 
    153 !----------------------------------------------------------------------- 
    154 !   noumef      =   43      !  number of fields 
    155 !SF   add_diag_swi=    1      !  1 -> diagnose distribution in thickness space 
    156 !SF                           !  0 -> only simple diagnostics 
    157 ! 
    158 !           !         title of the field           !  name     !   units   !  save  ! multipl. ! additive ! 
    159 !           !                                      !           !           ! or not !  factor  !  factor  ! 
    160 /  
    161150 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r4634 r4649  
    270270   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_sum   ! vertical surface melt 
    271271   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_res   ! production (growth+melt) due to limupdate 
     272   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   wfx_spr   ! snow precipitation on ice 
    272273 
    273274   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_bog     !: salt flux due to ice growth/melt                      [PSU/m2/s] 
     
    280281   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sfx_res     !: residual salt flux due to correction of ice thickness [PSU/m2/s] 
    281282 
    282    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh)  
    283    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me)  
    284    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_tot     !: total heat flux lost/gained by ice  
    285    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  
    286    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness 
     283   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bog     !: total heat flux causing bottom ice growth  
     284   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_bom     !: total heat flux causing bottom ice melt  
     285   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sum     !: total heat flux causing surface ice melt  
     286   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_opw     !: total heat flux causing open water ice formation 
     287   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dif     !: total heat flux causing Temp change in the ice  
    287288   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_snw     !: heat flux for snow melt  
    288    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  
    289289   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err     !: heat flux error after heat diffusion  
    290290   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_err_rem !: heat flux error after heat remapping  
    291291   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_in      !: heat flux available for thermo transformations  
    292292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_out     !: heat flux remaining at the end of thermo transformations  
     293 
     294   ! heat flux associated with ice-atmosphere mass exchange 
     295   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_sub     !: heat flux for sublimation  
     296   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_spr     !: heat flux of the snow precipitation  
     297 
     298   ! heat flux associated with ice-ocean mass exchange 
     299   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_thd     !: ice-ocean heat flux from thermo processes (limthd_dh)  
     300   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_dyn     !: ice-ocean heat flux from mecanical processes (limitd_me)  
     301   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   hfx_res     !: residual heat flux due to correction of ice thickness 
    293302 
    294303   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ftr_ice     !: transmitted solar radiation under ice 
     
    416425   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_trp_es   ! transport of snw enthalpy (W/m2) 
    417426   ! 
    418    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_heat_dhc1    ! snw/ice heat content variation   [W/m2]  
     427   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:)   ::   diag_heat_dhc    ! snw/ice heat content variation   [W/m2]  
    419428   ! 
    420429   INTEGER , PUBLIC ::   jiindx, jjindx        !: indexes of the debugging point 
     
    450459 
    451460      ii = ii + 1 
    452       ALLOCATE( sist     (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,      & 
    453          &      frld     (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,      & 
    454          &      wfx_snw  (jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,    & 
    455          &      wfx_bog(jpi,jpj)  , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     & 
    456          &      wfx_res(jpi,jpj)  , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) ,  qlead  (jpi,jpj) ,     & 
    457          &      fhtur    (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      & 
    458          &      sfx_res  (jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      & 
    459          &      sfx_bog  (jpi,jpj) , sfx_bom  (jpi,jpj) , sfx_sum  (jpi,jpj) ,  sfx_sni  (jpi,jpj) ,  sfx_opw  (jpi,jpj) ,   & 
    460          &      hfx_res  (jpi,jpj) , hfx_snw  (jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj), hfx_err_rem(jpi,jpj), & 
    461          &      hfx_in   (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,  & 
    462          &      hfx_tot  (jpi,jpj) , hfx_thd  (jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj),  STAT=ierr(ii) ) 
     461      ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,      & 
     462         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,      & 
     463         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,    & 
     464         &      wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) ,     & 
     465         &      wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) ,  qlead  (jpi,jpj) ,     & 
     466         &      fhtur  (jpi,jpj) , ftr_ice(jpi,jpj,jpl) ,      & 
     467         &      sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) ,      & 
     468         &      sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) ,   & 
     469         &      hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , hfx_err_rem(jpi,jpj), & 
     470         &      hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,  & 
     471         &      hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) , & 
     472         &      hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) ,  STAT=ierr(ii) ) 
    463473 
    464474      ii = ii + 1 
     
    526536      ALLOCATE( dv_dt_thd(jpi,jpj,jpl) ,     & 
    527537         &      izero    (jpi,jpj,jpl)  , diag_trp_vi(jpi,jpj) , diag_trp_vs(jpi,jpj), diag_trp_ei(jpi,jpj), diag_trp_es(jpi,jpj),     &  
    528          &      diag_heat_dhc1(jpi,jpj) ,  STAT=ierr(ii) ) 
     538         &      diag_heat_dhc(jpi,jpj) ,  STAT=ierr(ii) ) 
    529539 
    530540      ice_alloc = MAXVAL( ierr(:) ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r4045 r4649  
    77   !!            3.0  ! 2004-06  (M. Vancoppenolle)   Energy Conservation  
    88   !!            4.0  ! 2011-02  (G. Madec)  add mpp considerations 
     9   !!             -   ! 2014-05  (C. Rousset) add lim_cons_hsm 
    910   !!---------------------------------------------------------------------- 
    1011#if defined key_lim3 
     
    1415   !!    lim_cons     :   checks whether energy, mass and salt are conserved  
    1516   !!---------------------------------------------------------------------- 
     17   USE phycst         ! physical constants 
    1618   USE par_ice        ! LIM-3 parameter 
    1719   USE ice            ! LIM-3 variables 
     
    2830   PUBLIC   lim_column_sum_energy 
    2931   PUBLIC   lim_cons_check 
     32   PUBLIC   lim_cons_hsm 
    3033 
    3134   !!---------------------------------------------------------------------- 
     
    151154   END SUBROUTINE lim_cons_check 
    152155 
     156 
     157   SUBROUTINE lim_cons_hsm( icount, cd_routine, zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b ) 
     158      !!------------------------------------------------------------------- 
     159      !!               ***  ROUTINE lim_cons_hsm *** 
     160      !! 
     161      !! ** Purpose : Test the conservation of heat, salt and mass for each routine 
     162      !! 
     163      !! ** Method  : 
     164      !!--------------------------------------------------------------------- 
     165      INTEGER         , INTENT(in)    :: icount      ! determine wether this is the beggining of the routine (0) or the end (1) 
     166      CHARACTER(len=*), INTENT(in)    :: cd_routine  ! name of the routine 
     167      REAL(wp)        , INTENT(inout) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
     168      REAL(wp)                        :: zvi,   zsmv,   zei,   zfs,   zfw,   zft 
     169      REAL(wp)                        :: zvmin, zamin, zamax  
     170 
     171      IF( icount == 0 ) THEN 
     172 
     173         zvi_b  = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
     174         zsmv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
     175         zei_b  = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
     176         zfw_b  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
     177            &                   wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) 
     178         zfs_b  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     179            &                   sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
     180         zft_b  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     181            &                 - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
     182 
     183      ELSEIF( icount == 1 ) THEN 
     184 
     185         zfs  = glob_sum(   ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) +  & 
     186            &                 sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zfs_b 
     187         zfw  = glob_sum( - ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) +  & 
     188            &                 wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) ) * area(:,:) * tms(:,:) ) - zfw_b 
     189         zft  = glob_sum(   ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:)  &  
     190            &               - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zft_b 
     191  
     192         zvi  = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zvi_b ) * r1_rdtice - zfw  
     193         zsmv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zsmv_b ) * r1_rdtice + ( zfs / rhoic ) 
     194         zei  =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zei_b * r1_rdtice + zft 
     195 
     196         zvmin = glob_min(v_i) 
     197         zamax = glob_max(SUM(a_i,dim=3)) 
     198         zamin = glob_min(a_i) 
     199        
     200         IF(lwp) THEN 
     201            IF ( ABS( zvi    ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (',cd_routine,') = ',(zvi * rday) 
     202            IF ( ABS( zsmv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (',cd_routine,') = ',(zsmv * rday) 
     203            IF ( ABS( zei    ) >  1.    ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (',cd_routine,') = ',(zei) 
     204            IF ( zvmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [m]          (',cd_routine,') = ',(zvmin) 
     205            IF( cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' .AND. zamax > amax+1.e-10 ) THEN 
     206                                          WRITE(numout,*) 'violation a_i>amax            (',cd_routine,') = ',zamax 
     207            ENDIF 
     208            IF ( zamin <  0.            ) WRITE(numout,*) 'violation a_i<0               (',cd_routine,') = ',zamin 
     209         ENDIF 
     210 
     211      ENDIF 
     212 
     213   END SUBROUTINE lim_cons_hsm 
     214 
    153215#else 
    154216   !!---------------------------------------------------------------------- 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r4634 r4649  
    3535   !!PUBLIC   lim_diahsb_rst   ! routine called by ice_init.F90 
    3636 
    37    REAL(wp) ::   frc_sal, frc_vol   ! global forcing trends 
    38    REAL(wp) ::   bg_grme            ! global ice growth+melt trends 
     37   REAL(dp) ::   frc_sal, frc_vol   ! global forcing trends 
     38   REAL(dp) ::   bg_grme            ! global ice growth+melt trends 
    3939   REAL(wp) ::   epsi06 = 1.e-6_wp  ! small number 
    4040 
     
    4747   !! Software governed by the CeCILL licence     (NEMOGCM/NEMO_CeCILL.txt) 
    4848   !!---------------------------------------------------------------------- 
     49 
    4950CONTAINS 
    5051 
     
    5758      !!--------------------------------------------------------------------------- 
    5859      !! 
    59       REAL(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 
    60       REAL(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni, zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn  
    61       REAL(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn, zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res  
    62       REAL(wp)   ::   zbg_hfx_dhc1, zbg_hfx_spr, zbg_hfx_qsr, zbg_hfx_qns 
    63       REAL(wp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_tot, zbg_hfx_out, zbg_hfx_in    
    64       REAL(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme  
    65       REAL(wp)   ::   z1_area, zcoef                     
    66       REAL(wp)   ::   zinda, zindb 
     60      REAL(dp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 
     61      REAL(dp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni, zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn  
     62      REAL(dp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn 
     63      REAL(dp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub   
     64      REAL(dp)   ::   zbg_hfx_dhc, zbg_hfx_spr 
     65      REAL(dp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in    
     66      REAL(dp)   ::   zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw 
     67      REAL(dp)   ::   z_frc_vol, z_frc_sal, z_bg_grme  
     68      REAL(dp)   ::   z1_area                     !    -     - 
     69      REAL(dp)   ::   zinda, zindb 
    6770      !!--------------------------------------------------------------------------- 
    6871      IF( nn_timing == 1 )   CALL timing_start('lim_diahsb') 
     
    7073      IF( numit == nstart ) CALL lim_diahsb_init  
    7174 
    72       z1_area = 1._wp / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 ) 
    73  
    74       zinda = MAX( 0._wp , SIGN( 1._wp , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) ) 
     75      ! 1/area 
     76      z1_area = 1.d0 / MAX( glob_sum( area(:,:) * tms(:,:) ), epsi06 ) 
     77 
     78      zinda = MAX( 0.d0 , SIGN( 1.d0 , glob_sum( area(:,:) * tms(:,:) ) - epsi06 ) ) 
    7579      ! ----------------------- ! 
    7680      ! 1 -  Content variations ! 
     
    7983      zbg_svo = glob_sum( vt_s(:,:) * area(:,:) * tms(:,:) ) ! volume snow 
    8084      zbg_are = glob_sum( at_i(:,:) * area(:,:) * tms(:,:) ) ! area 
    81       zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 )      * area(:,:) * tms(:,:) )  ! mean salt content 
     85      zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) )       ! mean salt content 
    8286      zbg_tem = glob_sum( ( tm_i(:,:) - rtt ) * vt_i(:,:) * area(:,:) * tms(:,:) )  ! mean temp content 
    8387 
    84       zcoef = zinda * z1_area * r1_rau0 * rday 
     88      !zbg_ihc = glob_sum( et_i(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 
     89      !zbg_shc = glob_sum( et_s(:,:) * area(:,:) * tms(:,:) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 
     90 
    8591      ! Volume 
    86       zbg_vfx     = glob_sum(     emp(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    87       zbg_vfx_bog = glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    88       zbg_vfx_opw = glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    89       zbg_vfx_sni = glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    90       zbg_vfx_dyn = glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    91       zbg_vfx_bom = glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    92       zbg_vfx_sum = glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    93       zbg_vfx_res = glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) ) * zcoef 
     92      zbg_vfx     = zinda * glob_sum(      emp(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     93      zbg_vfx_bog = zinda * glob_sum( wfx_bog(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     94      zbg_vfx_opw = zinda * glob_sum( wfx_opw(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     95      zbg_vfx_sni = zinda * glob_sum( wfx_sni(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     96      zbg_vfx_dyn = zinda * glob_sum( wfx_dyn(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     97      zbg_vfx_bom = zinda * glob_sum( wfx_bom(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     98      zbg_vfx_sum = zinda * glob_sum( wfx_sum(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     99      zbg_vfx_res = zinda * glob_sum( wfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     100      zbg_vfx_spr = zinda * glob_sum( wfx_spr(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     101      zbg_vfx_snw = zinda * glob_sum( wfx_snw(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     102      zbg_vfx_sub = zinda * glob_sum( wfx_sub(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
    94103 
    95104      ! Salt 
    96       zbg_sfx     = glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    97       zbg_sfx_bri = glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    98       zbg_sfx_res = glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    99       zbg_sfx_dyn = glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    100  
    101       zbg_sfx_bog = glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    102       zbg_sfx_opw = glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    103       zbg_sfx_sni = glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    104       zbg_sfx_bom = glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) ) * zcoef 
    105       zbg_sfx_sum = glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) ) * zcoef 
     105      zbg_sfx     = zinda * glob_sum(     sfx(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     106      zbg_sfx_bri = zinda * glob_sum( sfx_bri(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     107      zbg_sfx_res = zinda * glob_sum( sfx_res(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     108      zbg_sfx_dyn = zinda * glob_sum( sfx_dyn(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     109 
     110      zbg_sfx_bog = zinda * glob_sum( sfx_bog(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     111      zbg_sfx_opw = zinda * glob_sum( sfx_opw(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     112      zbg_sfx_sni = zinda * glob_sum( sfx_sni(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     113      zbg_sfx_bom = zinda * glob_sum( sfx_bom(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
     114      zbg_sfx_sum = zinda * glob_sum( sfx_sum(:,:) * area(:,:) * tms(:,:) ) * z1_area * r1_rau0 * rday 
    106115 
    107116      ! Heat budget 
    108117      zbg_ihc      = glob_sum( et_i(:,:) * 1.e-20 ) ! ice heat content  [1.e-20 J] 
    109118      zbg_shc      = glob_sum( et_s(:,:) * 1.e-20 ) ! snow heat content [1.e-20 J] 
    110       zbg_hfx_dhc1 = glob_sum( diag_heat_dhc1(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     119      zbg_hfx_dhc  = glob_sum( diag_heat_dhc(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    111120      zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    112       zbg_hfx_qsr  = glob_sum( qsr(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    113       zbg_hfx_qns  = glob_sum( qns(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    114121 
    115122      zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     
    118125      zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    119126      zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    120       zbg_hfx_tot  = glob_sum( hfx_tot(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     127      zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     128      zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     129      zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     130      zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     131      zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    121132      zbg_hfx_out  = glob_sum( hfx_out(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
    122133      zbg_hfx_in   = glob_sum(  hfx_in(:,:) * area(:,:) * tms(:,:) ) ! [in W] 
     
    127138      z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * area(:,:) * tms(:,:) ) ! volume fluxes 
    128139      z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * area(:,:) * tms(:,:) ) ! salt fluxes 
    129       z_bg_grme = glob_sum( ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 
    130                           &   wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) ) / rhoic * area(:,:) * tms(:,:) ) ! volume fluxes 
    131       ! 
    132       frc_vol  = frc_vol + z_frc_vol * rdt_ice 
    133       frc_sal  = frc_sal + z_frc_sal * rdt_ice 
    134       bg_grme  = bg_grme + z_bg_grme * rdt_ice 
    135              
     140      z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 
     141                          &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + wfx_sub(:,:) ) * area(:,:) * tms(:,:) ) ! volume fluxes 
     142      ! 
     143      frc_vol  = frc_vol  + z_frc_vol  * rdt_ice 
     144      frc_sal  = frc_sal  + z_frc_sal  * rdt_ice 
     145      bg_grme  = bg_grme  + z_bg_grme  * rdt_ice 
     146       
     147      ! difference 
     148      !frc_vol = zbg_ivo - frc_vol 
     149      !frc_sal = zbg_sal - frc_sal 
     150       
    136151      ! ----------------------- ! 
    137152      ! 3 - Diagnostics writing ! 
    138153      ! ----------------------- ! 
    139       zindb = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) ) 
     154      zindb = MAX( 0.d0 , SIGN( 1.d0 , zbg_ivo - epsi06 ) ) 
    140155      ! 
    141156      CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)          
     
    156171      CALL iom_put( 'ibgvfxsum' , zbg_vfx_sum                              )   ! volume flux surface melt      - 
    157172      CALL iom_put( 'ibgvfxres' , zbg_vfx_res                              )   ! volume flux resultant         - 
     173      CALL iom_put( 'ibgvfxspr' , zbg_vfx_spr                              )   ! volume flux from snow precip         - 
     174      CALL iom_put( 'ibgvfxsnw' , zbg_vfx_snw                              )   ! volume flux from snow melt         - 
     175      CALL iom_put( 'ibgvfxsub' , zbg_vfx_sub                              )   ! volume flux from sublimation         - 
    158176           
    159177      CALL iom_put( 'ibgsfx'    , zbg_sfx                                  )   ! salt flux         -(psu*m/day equivalent liquid)        
     
    167185      CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      - 
    168186 
    169       CALL iom_put( 'ibghfxdhc1', zbg_hfx_dhc1                            )   ! Heat content variation in snow and ice [W] 
     187      CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W] 
    170188      CALL iom_put( 'ibghfxspr' , zbg_hfx_spr                              )   ! Heat content of snow precip [W] 
    171       CALL iom_put( 'ibghfxqsr' , zbg_hfx_qsr                              )   !     solar fluxes used by snw/ice [W] 
    172       CALL iom_put( 'ibghfxqns' , zbg_hfx_qns                              )   ! non solar fluxes used by snw/ice [W] 
    173189 
    174190      CALL iom_put( 'ibghfxres' , zbg_hfx_res                              )   !  
     
    177193      CALL iom_put( 'ibghfxthd' , zbg_hfx_thd                              )   !  
    178194      CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw                              )   !  
    179       CALL iom_put( 'ibghfxtot' , zbg_hfx_tot                              )   !  
     195      CALL iom_put( 'ibghfxsum' , zbg_hfx_sum                              )   !  
     196      CALL iom_put( 'ibghfxbom' , zbg_hfx_bom                              )   !  
     197      CALL iom_put( 'ibghfxbog' , zbg_hfx_bog                              )   !  
     198      CALL iom_put( 'ibghfxdif' , zbg_hfx_dif                              )   !  
     199      CALL iom_put( 'ibghfxopw' , zbg_hfx_opw                              )   !  
    180200      CALL iom_put( 'ibghfxout' , zbg_hfx_out                              )   !  
    181       CALL iom_put( 'ibghfxin'  , zbg_hfx_in                              )   !  
     201      CALL iom_put( 'ibghfxin'  , zbg_hfx_in                               )   !  
    182202 
    183203      CALL iom_put( 'ibgfrcvol' , frc_vol * 1.e-9                          )   ! vol - forcing     (km3 equivalent liquid)  
    184204      CALL iom_put( 'ibgfrcsfx' , frc_sal * 1.e-9                          )   ! sal - forcing     (psu*km3 equivalent liquid)    
    185       CALL iom_put( 'ibgvolgrm' , bg_grme * rhoic * r1_rau0 * 1.e-9        )   ! vol growth + melt (km3 equivalent liquid)          
     205      CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9                )   ! vol growth + melt (km3 equivalent liquid)          
    186206 
    187207      ! 
     
    189209      ! 
    190210      IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb') 
    191       ! 
     211! 
    192212   END SUBROUTINE lim_diahsb 
    193213 
     
    223243      ! 2 - initial conservation variables ! 
    224244      ! ---------------------------------- ! 
     245      !frc_vol = 0.d0                                           ! volume       trend due to forcing 
     246      !frc_sal = 0.d0                                           ! salt content   -    -   -    -          
     247      !bg_grme = 0.d0                                           ! ice growth + melt volume trend 
    225248      ! 
    226249      CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files 
     
    256279           IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state ' 
    257280           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    258            frc_vol  = 0._wp                                            
    259            frc_sal  = 0._wp                                                   
    260            bg_grme  = 0._wp                                         
     281           frc_vol  = 0.d0                                            
     282           frc_sal  = 0.d0                                                   
     283           bg_grme  = 0.d0                                         
    261284       ENDIF    
    262285 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90

    r4634 r4649  
    3030   USE lib_fortran      ! glob_sum 
    3131   USE timing          ! Timing 
     32   USE limcons        ! conservation tests 
    3233 
    3334   IMPLICIT NONE 
     
    6667      REAL(wp), POINTER, DIMENSION(:)   ::   zmsk           ! i-averaged of tmask 
    6768      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_io, zv_io   ! ice-ocean velocity 
    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) 
     69      ! 
     70      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    7071     !!--------------------------------------------------------------------- 
    7172 
     
    7576      CALL wrk_alloc( jpj, zind, zmsk ) 
    7677 
    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( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    83          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    84       ENDIF 
    85       !- check conservation (C Rousset) 
    86       ! ------------------------------- 
    87  
    8878      IF( kt == nit000 )   CALL lim_dyn_init   ! Initialization (first time-step only) 
    8979 
    9080      IF( ln_limdyn ) THEN 
    9181         ! 
     82         ! conservation test 
     83         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     84 
    9285         old_u_ice(:,:) = u_ice(:,:) * tmu(:,:) 
    9386         old_v_ice(:,:) = v_ice(:,:) * tmv(:,:) 
     
    171164            END DO 
    172165         END DO 
     166         ! 
     167         ! conservation test 
     168         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    173169         ! 
    174170      ELSE      ! no ice dynamics : transmit directly the atmospheric stress to the ocean 
     
    224220      ENDIF 
    225221      ! 
    226       ! ------------------------------- 
    227       !- check conservation (C Rousset) 
    228       IF (ln_limdiahsb) THEN 
    229          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    230          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    231   
    232          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - ( zchk_fw / rhoic ) 
    233          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( 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  
    250222      CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 
    251223      CALL wrk_dealloc( jpj, zind, zmsk ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90

    r4634 r4649  
    2626   USE dom_ice          ! sea-ice domain 
    2727   USE in_out_manager   ! I/O manager 
    28    USE lbclnk           ! lateral boundary condition - MPP exchanges 
    2928   USE lib_mpp          ! MPP library 
    3029   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     
    390389      tn_ice (:,:,:) = t_su (:,:,:) 
    391390 
    392       !-------------------------------------------------------------------- 
    393       ! 4) Lateral boundary conditions                                    |  
    394       !-------------------------------------------------------------------- 
    395       DO jl = 1, jpl 
    396  
    397          CALL lbc_lnk( a_i(:,:,jl)  , 'T', 1. ) 
    398          CALL lbc_lnk( v_i(:,:,jl)  , 'T', 1. ) 
    399          CALL lbc_lnk( v_s(:,:,jl)  , 'T', 1. ) 
    400          CALL lbc_lnk( smv_i(:,:,jl), 'T', 1. ) 
    401          CALL lbc_lnk( oa_i(:,:,jl) , 'T', 1. ) 
    402  
    403          CALL lbc_lnk( ht_i(:,:,jl) , 'T', 1. ) 
    404          CALL lbc_lnk( ht_s(:,:,jl) , 'T', 1. ) 
    405          CALL lbc_lnk( sm_i(:,:,jl) , 'T', 1. ) 
    406          CALL lbc_lnk( o_i(:,:,jl)  , 'T', 1. ) 
    407          CALL lbc_lnk( t_su(:,:,jl) , 'T', 1. ) 
    408          CALL lbc_lnk( tn_ice(:,:,jl) , 'T', 1. ) 
    409          DO jk = 1, nlay_s 
    410             CALL lbc_lnk(t_s(:,:,jk,jl), 'T', 1. ) 
    411             CALL lbc_lnk(e_s(:,:,jk,jl), 'T', 1. ) 
    412          END DO 
    413          DO jk = 1, nlay_i 
    414             CALL lbc_lnk(t_i(:,:,jk,jl), 'T', 1. ) 
    415             CALL lbc_lnk(e_i(:,:,jk,jl), 'T', 1. ) 
    416          END DO 
    417          ! 
    418 !         a_i(:,:,jl) = tms(:,:) * a_i(:,:,jl) 
    419       END DO 
    420  
    421391      ELSE  
    422392         ! if ln_limini=false 
     
    449419         at_i (:,:) = at_i (:,:) + a_i (:,:,jl) 
    450420      END DO 
    451       CALL lbc_lnk( at_i , 'T', 1. ) 
    452 !      at_i(:,:) = tms(:,:) * at_i(:,:) 
    453421      ! 
    454422      !-------------------------------------------------------------------- 
    455       ! 5) Global ice variables for output diagnostics                    |  
     423      ! 4) Global ice variables for output diagnostics                    |  
    456424      !-------------------------------------------------------------------- 
    457425      u_ice (:,:)     = 0._wp 
     
    462430 
    463431      !-------------------------------------------------------------------- 
    464       ! 6) Moments for advection 
     432      ! 5) Moments for advection 
    465433      !-------------------------------------------------------------------- 
    466434 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r4634 r4649  
    2222   USE limthd_lac       ! LIM 
    2323   USE limvar           ! LIM 
    24    USE limcons          ! LIM 
    2524   USE in_out_manager   ! I/O manager 
    2625   USE lbclnk           ! lateral boundary condition - MPP exchanges 
     
    3332   USE limdiahsb 
    3433   USE timing          ! Timing 
     34   USE limcons        ! conservation tests 
    3535 
    3636   IMPLICIT NONE 
     
    143143      REAL(wp), POINTER, DIMENSION(:,:) ::   esnow_mlt       ! energy needed to melt snow in ocean (J m-2) 
    144144      REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final  !  ice volume summed over categories 
    145       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b  
    146       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     145      ! 
     146      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    147147      !!----------------------------------------------------------------------------- 
    148148      IF( nn_timing == 1 )  CALL timing_start('limitd_me') 
     
    158158 
    159159      IF( ln_limdyn ) THEN          !   Start ridging and rafting   ! 
    160       ! ------------------------------- 
    161       !- check conservation (C Rousset) 
    162       IF (ln_limdiahsb) THEN 
    163          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    164          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    165          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    166          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    167          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    168          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    169       ENDIF 
    170       !- check conservation (C Rousset) 
    171       ! ------------------------------- 
     160 
     161      ! conservation test 
     162      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    172163 
    173164      !-----------------------------------------------------------------------------! 
     
    355346            ! 5) Heat, salt and freshwater fluxes 
    356347            !-----------------------------------------------------------------------------! 
    357             wfx_snw(ji,jj) = wfx_snw(ji,jj) - msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
     348            wfx_snw(ji,jj) = wfx_snw(ji,jj) + msnow_mlt(ji,jj) * r1_rdtice     ! fresh water source for ocean 
    358349            hfx_dyn(ji,jj) = hfx_dyn(ji,jj) + esnow_mlt(ji,jj) * unit_fac / area(ji,jj) * r1_rdtice  ! heat sink for ocean (<0, W.m-2) 
    359350 
     
    425416      ENDIF 
    426417 
    427       ! ------------------------------- 
    428       !- check conservation (C Rousset) 
    429       IF (ln_limdiahsb) THEN 
    430          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    431          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    432          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    433   
    434          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    435          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    436          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    437  
    438          zchk_vmin = glob_min(v_i) 
    439          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    440          zchk_amin = glob_min(a_i) 
    441         
    442          IF(lwp) THEN 
    443             IF ( ABS( zchk_v_i   ) >  1.e-2 ) WRITE(numout,*) 'violation volume [kg/day]     (limitd_me) = ',(zchk_v_i * rday) 
    444             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_me) = ',(zchk_smv * rday) 
    445             IF ( ABS( zchk_e_i   ) >  1.e-4 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limitd_me) = ',(zchk_e_i) 
    446             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_me) = ',(zchk_vmin * 1.e-3) 
    447             IF ( zchk_amax >  kamax+epsi10  ) WRITE(numout,*) 'violation a_i>amax            (limitd_me) = ',zchk_amax 
    448             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_me) = ',zchk_amin 
    449          ENDIF 
    450       ENDIF 
    451       !- check conservation (C Rousset) 
    452       ! ------------------------------- 
     418      ! conservation test 
     419      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    453420 
    454421      ENDIF  ! ln_limdyn=.true. 
     
    11081075             
    11091076            sfx_dyn(ji,jj) = sfx_dyn(ji,jj) - smsw(ji,jj) * rhoic * r1_rdtice 
    1110             wfx_dyn(ji,jj) = wfx_dyn(ji,jj) + vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids              
    1111             ! MV HC 2014 this previous line seems ok, i'm not sure at this moment of the sign convention 
     1077            wfx_dyn(ji,jj) = wfx_dyn(ji,jj) - vsw (ji,jj) * rhoic * r1_rdtice   ! gurvan: increase in ice volume du to seawater frozen in voids              
    11121078 
    11131079            !------------------------------------             
     
    14401406      CALL wrk_alloc( jpi, jpj, zmask ) 
    14411407 
     1408      ! to be sure that at_i is the sum of a_i(jl) 
     1409      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
     1410 
    14421411      DO jl = 1, jpl 
    1443  
    14441412         !----------------------------------------------------------------- 
    14451413         ! Count categories to be zapped. 
    1446          ! Abort model in case of negative area. 
    14471414         !----------------------------------------------------------------- 
    14481415         icells = 0 
     
    14501417         DO jj = 1, jpj 
    14511418            DO ji = 1, jpi 
    1452 !               IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <  0._wp   ) .OR.   & 
    1453 !                  & ( a_i(ji,jj,jl) >  0._wp   .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
    1454 !                  & ( v_i(ji,jj,jl) == 0._wp   .AND. a_i(ji,jj,jl) >  0._wp   ) .OR.   & 
    1455 !                  & ( v_i(ji,jj,jl) >  0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
    1456                IF(  ( a_i(ji,jj,jl) >= -epsi10 .AND. a_i(ji,jj,jl) <= epsi10  ) .OR.   & 
    1457                   & ( v_i(ji,jj,jl) >= 0._wp   .AND. v_i(ji,jj,jl) <= epsi10  ) )   zmask(ji,jj) = 1._wp 
     1419               IF( a_i(ji,jj,jl) <= epsi10 .OR. v_i(ji,jj,jl) <= epsi10 .OR. at_i(ji,jj) <= epsi10 ) THEN 
     1420                  zmask(ji,jj) = 1._wp 
     1421               ENDIF 
    14581422            END DO 
    14591423         END DO 
     
    14701434                  zei  = e_i(ji,jj,jk,jl) 
    14711435                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) 
     1436                  t_i(ji,jj,jk,jl) = t_i(ji,jj,jk,jl) * ( 1._wp - zmask(ji,jj) ) + rtt * zmask(ji,jj) 
    14721437                  ! update exchanges with ocean 
    14731438                  hfx_res(ji,jj)   = hfx_res(ji,jj) + ( e_i(ji,jj,jk,jl) - zei ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     
    14971462               ! zap ice and snow volume, add water and salt to ocean 
    14981463               !----------------------------------------------------------------- 
    1499  
    1500                !           xtmp = (rhoi*vicen(i,j,n) + rhos*vsnon(i,j,n)) / dt 
    1501                !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj)                  )   & 
    1502                !                                            * rhosn * v_s(ji,jj,jl) * r1_rdtice 
    1503                !           sfx_res(ji,jj)  = sfx_res(ji,jj) + ( sss_m(ji,jj) - sm_i(ji,jj,jl) )   &  
    1504                !                                            * rhoic * v_i(ji,jj,jl) * r1_rdtice 
    1505                !           sfx (i,j)      = sfx (i,j)      + xtmp 
    1506  
    1507                ato_i(ji,jj)    = a_i  (ji,jj,jl) *       zmask(ji,jj)   + ato_i(ji,jj) 
     1464               ato_i(ji,jj)    = a_i  (ji,jj,jl) *           zmask(ji,jj)   + ato_i(ji,jj) 
    15081465               a_i  (ji,jj,jl) = a_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
    15091466               v_i  (ji,jj,jl) = v_i  (ji,jj,jl) * ( 1._wp - zmask(ji,jj) ) 
     
    15201477               ! update exchanges with ocean 
    15211478               sfx_res(ji,jj)  = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    1522                wfx_res(ji,jj)  = wfx_res(ji,jj) + ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
    1523                wfx_snw(ji,jj)  = wfx_snw(ji,jj) + ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
     1479               wfx_res(ji,jj)  = wfx_res(ji,jj) - ( v_i(ji,jj,jl)   - zvi  ) * rhoic * r1_rdtice 
     1480               wfx_snw(ji,jj)  = wfx_snw(ji,jj) - ( v_s(ji,jj,jl)   - zvs  ) * rhosn * r1_rdtice 
    15241481               hfx_res(ji,jj)  = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
    15251482            END DO 
    15261483         END DO 
    1527          ! 
    1528       END DO                 ! jl  
     1484      END DO ! jl  
     1485 
     1486      ! to be sure that at_i is the sum of a_i(jl) 
     1487      at_i(:,:) = SUM( a_i(:,:,:), dim=3 ) 
    15291488      ! 
    15301489      CALL wrk_dealloc( jpi, jpj, zmask ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90

    r4634 r4649  
    3535   USE lib_fortran      ! to use key_nosignedzero 
    3636   USE timing          ! Timing 
     37   USE limcons        ! conservation tests 
    3738 
    3839   IMPLICIT NONE 
     
    6667      ! 
    6768      INTEGER ::   ji,jj, jk, jl, ja, jm, jbnd1, jbnd2   ! ice types    dummy loop index          
    68       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 
    69       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     69      ! 
     70      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    7071      !!------------------------------------------------------------------ 
    7172      IF( nn_timing == 1 )  CALL timing_start('limitd_th') 
    7273 
    73       ! ------------------------------- 
    74       !- check conservation (C Rousset) 
    75       IF (ln_limdiahsb) THEN 
    76          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    77          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    78          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    79          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    80          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    81          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    82        ENDIF 
    83       !- check conservation (C Rousset) 
    84       ! ------------------------------- 
     74      ! conservation test 
     75      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    8576 
    8677      IF( kt == nit000 .AND. lwp ) THEN 
     
    10798      !  3) Add frazil ice growing in leads. 
    10899      !------------------------------------------------------------------------------| 
    109  
    110100      CALL lim_thd_lac 
    111101      CALL lim_var_glo2eqv    ! only for info 
     
    143133      ENDIF 
    144134      ! 
    145       ! ------------------------------- 
    146       !- check conservation (C Rousset) 
    147       IF( ln_limdiahsb ) THEN 
    148          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    149          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    150          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    151   
    152          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    153          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    154          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    155  
    156          zchk_vmin = glob_min(v_i) 
    157          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    158          zchk_amin = glob_min(a_i) 
    159  
    160          IF(lwp) THEN 
    161             IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limitd_th) = ',(zchk_v_i * rday) 
    162             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limitd_th) = ',(zchk_smv * rday) 
    163             IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limitd_th) = ',(zchk_e_i) 
    164             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limitd_th) = ',(zchk_vmin * 1.e-3) 
    165             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limitd_th) = ',zchk_amax 
    166             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limitd_th) = ',zchk_amin 
    167          ENDIF 
    168        ENDIF 
    169       !- check conservation (C Rousset) 
    170       ! ------------------------------- 
     135      ! conservation test 
     136      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_th', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    171137      ! 
    172138     IF( nn_timing == 1 )  CALL timing_stop('limitd_th') 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r4634 r4649  
    4242   USE lib_fortran      ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
    4343   USE traqsr           ! clem: add penetration of solar flux into the calculation of heat budget 
     44   USE iom 
    4445 
    4546   IMPLICIT NONE 
     
    104105      INTEGER  ::   ji, jj, jl, jk           ! dummy loop indices 
    105106      REAL(wp) ::   zinda, zemp      ! local scalars 
    106       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
    107       REAL(wp) ::   ztmelts         ! clem 2014: for HC diags 
    108  
    109107      REAL(wp) ::   zf_mass         ! Heat flux associated with mass exchange ice->ocean (W.m-2) 
    110108      REAL(wp) ::   zfcm1           ! New solar flux received by the ocean 
     109      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zalb, zalbp     ! 2D/3D workspace 
    111110      !!--------------------------------------------------------------------- 
    112111       
    113112      IF( lk_cpl )   CALL wrk_alloc( jpi, jpj, jpl, zalb, zalbp ) 
     113 
     114      ! make calls for heat fluxes before it is modified 
     115      CALL iom_put( "qsr_oce" , qsr(:,:) * pfrld(:,:) )   !     solar flux at ocean surface 
     116      CALL iom_put( "qns_oce" , qns(:,:) * pfrld(:,:) )   ! non-solar flux at ocean surface 
     117      CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux at ice surface 
     118      CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  ! non-solar flux at ice surface 
     119      CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * old_a_i(:,:,:), dim=3 ) )  !     solar flux transmitted thru ice 
     120      CALL iom_put( "qt_oce"  , ( qsr(:,:) + qns(:,:) ) * pfrld(:,:) )   
     121      CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) ) * old_a_i(:,:,:), dim=3 ) ) 
    114122 
    115123      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
     
    165173            IF( lk_cpl ) THEN  
    166174               zemp = - emp_tot(ji,jj) + emp_ice(ji,jj) * ( 1. - pfrld(ji,jj) )    &   ! 
    167                   &   - wfx_snw(ji,jj) 
     175                  &   + wfx_snw(ji,jj) 
    168176            ELSE 
    169177               zemp =   emp(ji,jj)     *           pfrld(ji,jj)            &   ! evaporation over oceanic fraction 
     
    176184 
    177185            ! mass flux at the ocean/ice interface 
    178             fmmflx(ji,jj) = wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model 
    179             emp(ji,jj)    = zemp + wfx_ice(ji,jj) + wfx_snw(ji,jj)     ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
     186            fmmflx(ji,jj) = - wfx_ice(ji,jj) * rdt_ice                   ! F/M mass flux save at least for biogeochemical model 
     187            emp(ji,jj)    = zemp - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_sub(ji,jj)   ! mass flux + F/M mass flux (always ice/ocean mass exchange) 
    180188             
    181189         END DO 
     
    213221      ENDIF 
    214222 
    215       ! ------------------------------------------------- 
    216       ! C. Rousset Begin Diagnostics for heat in W/m2 
    217       ! ------------------------------------------------- 
    218       DO jj = 1, jpj 
    219          DO ji = 1, jpi             
    220             diag_heat_dhc1(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  &  
    221                &                      SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj)    
    222          END DO 
    223       END DO 
    224       ! ------------------------------------------------- 
    225       ! C. Rousset End Diagnostics 
    226       ! ------------------------------------------------- 
    227223 
    228224      IF(ln_ctl) THEN 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd.F90

    r4634 r4649  
    4444   USE timing         ! Timing 
    4545   USE cpl_oasis3, ONLY : lk_cpl 
     46   USE limcons        ! conservation tests 
    4647 
    4748   IMPLICIT NONE 
     
    8687      INTEGER  :: nbpb             ! nb of icy pts for thermo. cal. 
    8788      INTEGER  :: ii, ij           ! temporary dummy loop index 
    88       REAL(wp) :: zfric_umin = 5e-03_wp    ! lower bound for the friction velocity 
    89       REAL(wp) :: zfric_umax = 2e-02_wp    ! upper bound for the friction velocity 
    90       REAL(wp) :: zinda, zindb, zfric_u     ! local scalar 
    91       REAL(wp) :: zareamin  !    -         - 
    92       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b  
    93       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
    94       REAL(wp) :: zqld, zqfr 
     89      REAL(wp) :: zfric_umin = 0._wp        ! lower bound for the friction velocity (cice value=5.e-04) 
     90      REAL(wp) :: zch        = 0.0057_wp    ! heat transfer coefficient 
     91      REAL(wp) :: zinda, zindb, zareamin  
     92      REAL(wp) :: zfric_u, zqld, zqfr 
    9593      REAL(wp), POINTER, DIMENSION(:) :: zdq, zq_ini, zhfx, zqfx 
    9694      REAL(wp)                        :: zhfx_err, ztest 
     95      ! 
     96      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    9797      !!------------------------------------------------------------------- 
    9898      IF( nn_timing == 1 )  CALL timing_start('limthd') 
     
    103103      zdq(:) = 0._wp ; zq_ini(:) = 0._wp ; zhfx(:) = 0._wp ; zqfx(:) = 0._wp       
    104104 
    105       ! ------------------------------- 
    106       !- check conservation (C Rousset) 
    107       IF (ln_limdiahsb) THEN 
    108          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    109          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    110          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    111          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    112          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    113          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    114       ENDIF 
    115       !- check conservation (C Rousset) 
    116       ! ------------------------------- 
     105      ! conservation test 
     106      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    117107 
    118108      !------------------------------------------------------------------------------! 
     
    158148!CDIR NOVERRCHK 
    159149         DO ji = 1, jpi 
    160             zinda          = tms(ji,jj) * ( 1.0 - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
     150            zinda          = tms(ji,jj) * ( 1._wp - MAX( 0._wp , SIGN( 1._wp , - at_i(ji,jj) + epsi10 ) ) ) ! 0 if no ice 
    161151            ! 
    162152            !           !  solar irradiance transmission at the mixed layer bottom and used in the lead heat budget 
     
    165155            !           !  net downward heat flux from the ice to the ocean, expressed as a function of ocean  
    166156            !           !  temperature and turbulent mixing (McPhee, 1992) 
    167  
    168             ! friction velocity 
    169             zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin )  
    170  
    171             !-- Energy from the turbulent oceanic heat flux. here the drag will depend on ice thickness and type (0.006) 
    172             fhtur(ji,jj)  = zinda * rau0 * rcp * 0.006  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ! W.m-2  
    173             ! clem: why not the following?  
    174             !fhtur(ji,jj)  = zinda * rau0 * rcp * 0.006  * SQRT( ust2s(ji,jj) ) * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 
    175  
    176             !-- Energy received in the lead, zqld is defined everywhere (J.m-2) 
    177             !   It includes turbulent ocean heat flux (only in the leads, the rest is used for bottom melting) 
    178             zqld =  tms(ji,jj) * rdt_ice *                               & 
    179                &  ( pfrld(ji,jj)        * ( qsr(ji,jj) * oatte(ji,jj)           &   ! solar heat + clem modif 
    180                &                          + qns(ji,jj)                          &   ! non solar heat 
    181                &                          + fhtur(ji,jj) )                      &   ! turbulent ice-ocean heat (0 if no ice) 
     157            ! 
     158            ! --- Energy received in the lead, zqld is defined everywhere (J.m-2) --- ! 
     159            zqld =  tms(ji,jj) * rdt_ice *                                       & 
     160               &  ( pfrld(ji,jj)         * ( qsr(ji,jj) * oatte(ji,jj)           &   ! solar heat + clem modif 
     161               &                           + qns(ji,jj) )                        &   ! non solar heat 
    182162               ! latent heat of precip (note that precip is included in qns but not in qns_ice) 
    183163               &    + ( pfrld(ji,jj)**betas - pfrld(ji,jj) ) * sprecip(ji,jj) * ( cpic * ( MIN( tatm_ice(ji,jj), rt0_snow ) - rtt ) - lfus )  & 
    184164               &    + ( 1._wp - pfrld(ji,jj) ) * ( tprecip(ji,jj) - sprecip(ji,jj) ) * rcp * ( tatm_ice(ji,jj) - rtt ) ) 
    185165 
    186             !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) 
     166            !-- Energy needed to bring ocean surface layer until its freezing (<0, J.m-2) --- ! 
    187167            zqfr = tms(ji,jj) * rau0 * rcp * fse3t_m(ji,jj,1) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) 
    188168 
     
    192172            ! If there is ice and leads are warming, then transfer energy from the lead budget and use it for bottom melting  
    193173            IF( at_i(ji,jj) > epsi10 .AND. zqld > 0._wp ) THEN 
    194                fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by a_i since this is (re)multiplied by a_i in limthd_dh.F90 
     174               fhld (ji,jj) = zqld * r1_rdtice / at_i(ji,jj) ! divided by at_i since this is (re)multiplied by a_i in limthd_dh.F90 
    195175               qlead(ji,jj) = 0._wp 
    196176            ENDIF 
    197177            ! 
    198             IF( qlead(ji,jj) == 0._wp )  zqld = 0._wp ; zqfr = 0._wp 
    199             ! 
     178            !-- Energy from the turbulent oceanic heat flux --- ! 
     179            !clem zfric_u        = MAX ( MIN( SQRT( ust2s(ji,jj) ) , zfric_umax ) , zfric_umin ) 
     180            zfric_u      = MAX( SQRT( ust2s(ji,jj) ), zfric_umin )  
     181            fhtur(ji,jj) = MAX( 0._wp, zinda * rau0 * rcp * zch  * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) ) ! W.m-2  
     182            ! upper bound for fhtur: we do not want SST to drop below Tfreeze.  
     183            ! So we say that the heat retrieved from the ocean (fhtur+fhld) must be < to the heat necessary to reach Tfreeze (zqfr)    
     184            ! This is not a clean budget, so that should be corrected at some point 
     185            fhtur(ji,jj) = zinda * MIN( fhtur(ji,jj), - fhld(ji,jj) - zqfr * r1_rdtice / MAX( at_i(ji,jj), epsi10 ) ) 
     186 
    200187            ! ----------------------------------------- 
    201188            ! Net heat flux on top of ice-ocean [W.m-2] 
     
    317304            CALL tab_2d_1d( nbpb, wfx_sum_1d (1:nbpb), wfx_sum         , jpi, jpj, npb(1:nbpb) ) 
    318305            CALL tab_2d_1d( nbpb, wfx_sni_1d (1:nbpb), wfx_sni         , jpi, jpj, npb(1:nbpb) ) 
     306            CALL tab_2d_1d( nbpb, wfx_res_1d (1:nbpb), wfx_res         , jpi, jpj, npb(1:nbpb) ) 
     307            CALL tab_2d_1d( nbpb, wfx_spr_1d (1:nbpb), wfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    319308 
    320309            CALL tab_2d_1d( nbpb, sfx_bog_1d (1:nbpb), sfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     
    323312            CALL tab_2d_1d( nbpb, sfx_sni_1d (1:nbpb), sfx_sni         , jpi, jpj, npb(1:nbpb) ) 
    324313            CALL tab_2d_1d( nbpb, sfx_bri_1d (1:nbpb), sfx_bri         , jpi, jpj, npb(1:nbpb) ) 
     314            CALL tab_2d_1d( nbpb, sfx_res_1d (1:nbpb), sfx_res         , jpi, jpj, npb(1:nbpb) ) 
    325315 
    326316            CALL tab_2d_1d( nbpb, iatte_1d   (1:nbpb), iatte           , jpi, jpj, npb(1:nbpb) )  
     
    329319            CALL tab_2d_1d( nbpb, hfx_thd_1d (1:nbpb), hfx_thd         , jpi, jpj, npb(1:nbpb) ) 
    330320            CALL tab_2d_1d( nbpb, hfx_spr_1d (1:nbpb), hfx_spr         , jpi, jpj, npb(1:nbpb) ) 
    331             CALL tab_2d_1d( nbpb, hfx_tot_1d (1:nbpb), hfx_tot         , jpi, jpj, npb(1:nbpb) ) 
     321            CALL tab_2d_1d( nbpb, hfx_sum_1d (1:nbpb), hfx_sum         , jpi, jpj, npb(1:nbpb) ) 
     322            CALL tab_2d_1d( nbpb, hfx_bom_1d (1:nbpb), hfx_bom         , jpi, jpj, npb(1:nbpb) ) 
     323            CALL tab_2d_1d( nbpb, hfx_bog_1d (1:nbpb), hfx_bog         , jpi, jpj, npb(1:nbpb) ) 
     324            CALL tab_2d_1d( nbpb, hfx_dif_1d (1:nbpb), hfx_dif         , jpi, jpj, npb(1:nbpb) ) 
     325            CALL tab_2d_1d( nbpb, hfx_opw_1d (1:nbpb), hfx_opw         , jpi, jpj, npb(1:nbpb) ) 
    332326            CALL tab_2d_1d( nbpb, hfx_snw_1d (1:nbpb), hfx_snw         , jpi, jpj, npb(1:nbpb) ) 
    333327            CALL tab_2d_1d( nbpb, hfx_sub_1d (1:nbpb), hfx_sub         , jpi, jpj, npb(1:nbpb) ) 
     
    365359               ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1 
    366360               hfx_in (ii,ij) = hfx_in (ii,ij) + a_i_b(ji) * ( qsr_ice_1d(ji) + qns_ice_1d(ji) ) 
    367  
    368361            END DO 
    369362 
     
    379372            CALL lim_thd_dh( 1, nbpb, jl )     
    380373 
    381             ! --- Ice/Snow enthalpy remapping --- ! 
    382             CALL lim_thd_ent( 1, nbpb, jl )  
     374            ! --- Ice enthalpy remapping --- ! 
     375            CALL lim_thd_ent( 1, nbpb, jl, q_i_b(1:nbpb,:) )  
    383376            !                                 
    384377            ! --- diag error on heat remapping - PART 2 --- ! 
     
    391384 
    392385            !---------------------------------! 
    393             ! Ice salinity                    ! 
     386            ! --- Ice salinity --- ! 
    394387            !---------------------------------! 
    395388            CALL lim_thd_sal( 1, nbpb )     
    396389 
    397             !           CALL lim_thd_enmelt(1,nbpb)   ! computes sea ice energy of melting 
     390            !---------------------------------! 
     391            ! --- temperature update --- ! 
     392            !---------------------------------! 
     393            CALL lim_thd_temp( 1, nbpb ) 
     394 
    398395            !-------------------------------- 
    399396            ! 4.4) Move 1D to 2D vectors 
     
    424421               CALL tab_1d_2d( nbpb, wfx_sum       , npb, wfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    425422               CALL tab_1d_2d( nbpb, wfx_sni       , npb, wfx_sni_1d(1:nbpb)   , jpi, jpj ) 
     423               CALL tab_1d_2d( nbpb, wfx_res       , npb, wfx_res_1d(1:nbpb)   , jpi, jpj ) 
     424               CALL tab_1d_2d( nbpb, wfx_spr       , npb, wfx_spr_1d(1:nbpb)   , jpi, jpj ) 
    426425 
    427426               CALL tab_1d_2d( nbpb, sfx_bog       , npb, sfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     
    429428               CALL tab_1d_2d( nbpb, sfx_sum       , npb, sfx_sum_1d(1:nbpb)   , jpi, jpj ) 
    430429               CALL tab_1d_2d( nbpb, sfx_sni       , npb, sfx_sni_1d(1:nbpb)   , jpi, jpj ) 
     430               CALL tab_1d_2d( nbpb, sfx_res       , npb, sfx_res_1d(1:nbpb)   , jpi, jpj ) 
    431431            ! 
    432432            IF( num_sal == 2 ) THEN 
     
    436436              CALL tab_1d_2d( nbpb, hfx_thd       , npb, hfx_thd_1d(1:nbpb)   , jpi, jpj ) 
    437437              CALL tab_1d_2d( nbpb, hfx_spr       , npb, hfx_spr_1d(1:nbpb)   , jpi, jpj ) 
    438               CALL tab_1d_2d( nbpb, hfx_tot       , npb, hfx_tot_1d(1:nbpb)   , jpi, jpj ) 
     438              CALL tab_1d_2d( nbpb, hfx_sum       , npb, hfx_sum_1d(1:nbpb)   , jpi, jpj ) 
     439              CALL tab_1d_2d( nbpb, hfx_bom       , npb, hfx_bom_1d(1:nbpb)   , jpi, jpj ) 
     440              CALL tab_1d_2d( nbpb, hfx_bog       , npb, hfx_bog_1d(1:nbpb)   , jpi, jpj ) 
     441              CALL tab_1d_2d( nbpb, hfx_dif       , npb, hfx_dif_1d(1:nbpb)   , jpi, jpj ) 
     442              CALL tab_1d_2d( nbpb, hfx_opw       , npb, hfx_opw_1d(1:nbpb)   , jpi, jpj ) 
    439443              CALL tab_1d_2d( nbpb, hfx_snw       , npb, hfx_snw_1d(1:nbpb)   , jpi, jpj ) 
    440444              CALL tab_1d_2d( nbpb, hfx_sub       , npb, hfx_sub_1d(1:nbpb)   , jpi, jpj ) 
     
    521525      ENDIF 
    522526      ! 
    523       ! ------------------------------- 
    524       !- check conservation (C Rousset) 
    525       IF (ln_limdiahsb) THEN 
    526          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    527          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    528          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    529   
    530          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    531          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    532          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    533  
    534          zchk_vmin = glob_min(v_i) 
    535          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    536          zchk_amin = glob_min(a_i) 
    537         
    538          IF(lwp) THEN 
    539             IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limthd) = ',(zchk_v_i * rday) 
    540             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limthd) = ',(zchk_smv * rday) 
    541             IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limthd) = ',(zchk_e_i) 
    542             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limthd) = ',(zchk_vmin * 1.e-3) 
    543             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limthd) = ',zchk_amax 
    544             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limthd) = ',zchk_amin 
    545          ENDIF 
    546       ENDIF 
    547       !- check conservation (C Rousset) 
    548       ! ------------------------------- 
     527      ! conservation test 
     528      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limthd', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    549529      ! 
    550530      CALL wrk_dealloc( jpij, zdq, zq_ini, zhfx, zqfx ) 
     
    558538      !!                   ***  ROUTINE lim_thd_enmelt ***  
    559539      !!                  
    560       !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) 
     540      !! ** Purpose :   Computes sea ice energy of melting q_i (J.m-3) from temperature 
    561541      !! 
    562542      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     
    565545      !! 
    566546      INTEGER  ::   ji, jk   ! dummy loop indices 
    567       REAL(wp) ::   ztmelts  ! local scalar  
     547      REAL(wp) ::   ztmelts, zindb  ! local scalar  
    568548      !!------------------------------------------------------------------- 
    569549      ! 
    570550      DO jk = 1, nlay_i             ! Sea ice energy of melting 
    571551         DO ji = kideb, kiut 
    572             ztmelts      =  - tmut  * s_i_b(ji,jk) + rtt  
    573             q_i_b(ji,jk) =    rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                 & 
    574                &                      + lfus * ( 1.0 - (ztmelts-rtt) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   & 
    575                &                      - rcp  * ( ztmelts-rtt  )  )  
     552            ztmelts      = - tmut  * s_i_b(ji,jk) + rtt  
     553            zindb        = MAX( 0._wp , SIGN( 1._wp , -(t_i_b(ji,jk) - rtt) - epsi10 ) ) 
     554            q_i_b(ji,jk) = rhoic * ( cpic * ( ztmelts - t_i_b(ji,jk) )                                             & 
     555               &                   + lfus * ( 1.0 - zindb * ( ztmelts-rtt ) / MIN( t_i_b(ji,jk)-rtt, -epsi10 ) )   & 
     556               &                   - rcp  *                 ( ztmelts-rtt )  )  
    576557         END DO 
    577558      END DO 
     
    584565   END SUBROUTINE lim_thd_enmelt 
    585566 
     567   SUBROUTINE lim_thd_temp( kideb, kiut ) 
     568      !!----------------------------------------------------------------------- 
     569      !!                   ***  ROUTINE lim_thd_temp ***  
     570      !!                  
     571      !! ** Purpose :   Computes sea ice temperature (Kelvin) from enthalpy 
     572      !! 
     573      !! ** Method  :   Formula (Bitz and Lipscomb, 1999) 
     574      !!------------------------------------------------------------------- 
     575      INTEGER, INTENT(in) ::   kideb, kiut   ! bounds for the spatial loop 
     576      !! 
     577      INTEGER  ::   ji, jk   ! dummy loop indices 
     578      REAL(wp) ::   ztmelts, zswitch, zaaa, zbbb, zccc, zdiscrim  ! local scalar  
     579      !!------------------------------------------------------------------- 
     580      ! Recover ice temperature 
     581      DO jk = 1, nlay_i 
     582         DO ji = kideb, kiut 
     583            ztmelts       =  -tmut * s_i_b(ji,jk) + rtt 
     584            ! Conversion q(S,T) -> T (second order equation) 
     585            zaaa          =  cpic 
     586            zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 
     587            zccc          =  lfus * ( ztmelts - rtt ) 
     588            zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 
     589            t_i_b(ji,jk)  =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 
     590             
     591            ! mask temperature 
     592            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )  
     593            t_i_b(ji,jk) =  zswitch * t_i_b(ji,jk) + ( 1._wp - zswitch ) * rtt 
     594         END DO  
     595      END DO  
     596 
     597   END SUBROUTINE lim_thd_temp 
    586598 
    587599   SUBROUTINE lim_thd_init 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90

    r4634 r4649  
    7272      INTEGER  ::   ji , jk        ! dummy loop indices 
    7373      INTEGER  ::   ii, ij         ! 2D corresponding indices to ji 
    74       INTEGER  ::   i_ice_switch   ! ice thickness above a certain treshold or not 
    7574      INTEGER  ::   iter 
    7675 
     
    114113 
    115114      ! mass and salt flux (clem) 
    116       REAL(wp) :: zdvres, zswitch_sal 
     115      REAL(wp) :: zdvres, zswitch_sal, zswitch 
    117116 
    118117      ! Heat conservation  
     
    180179         zq_su (ji) = MAX( 0._wp, zfdum     * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) ) 
    181180         zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice ) 
    182       END DO ! ji 
     181      END DO 
    183182 
    184183      ! 
     
    192191            hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice 
    193192            ! Contribution to mass flux 
    194             wfx_snw_1d(ji) =  wfx_snw_1d(ji) - rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice 
     193            wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice 
    195194            ! updates 
    196195            ht_s_b(ji)   = 0._wp 
     
    247246         zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn 
    248247         ! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K) 
    249          zqprec    (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
     248         zqprec   (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )    
    250249         IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp 
    251250         ! heat flux from snow precip (>0, W.m-2) 
    252251         hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice 
     252         ! mass flux, <0 
     253         wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_b(ji) * zdh_s_pre(ji) * r1_rdtice 
    253254         ! update thickness 
    254255         ht_s_b    (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) ) 
     
    258259         !--------------------- 
    259260         ! thickness change 
     261         IF( zdh_s_pre(ji) > 0._wp ) THEN 
    260262         zindq          = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) ) 
    261263         zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 ) 
    262264         zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting  
    263          ! Heat flux associated with snow melt of falling snow (W.m-2, <0) 
    264          hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice  
    265265         ! heat used to melt snow (W.m-2, >0) 
    266          hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice 
    267          ! snow melting only = water into the ocean (then without snow precip) 
    268          wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice 
     266         hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice 
     267         ! snow melting only = water into the ocean (then without snow precip), >0 
     268         wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice 
    269269          
    270270         ! updates available heat + thickness 
     
    275275         ! clem debug: variation of enthalpy (J.m-2) 
    276276         dq_s(ji) = dq_s(ji) + ( zdh_s_pre(ji) + zdh_s_mel(ji) ) * zqprec(ji)   
    277  
     277         ENDIF 
    278278      END DO 
    279279 
     
    288288            zdeltah  (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting 
    289289            zdh_s_mel(ji)    = zdh_s_mel(ji) + zdeltah(ji,jk)     
    290             ! heat flux associated with snow melt(W.m-2, <0) 
    291             hfx_snw_1d(ji)   = hfx_snw_1d(ji) + zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice 
    292290            ! heat used to melt snow(W.m-2, >0) 
    293             hfx_tot_1d(ji)   = hfx_tot_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice  
     291            hfx_snw_1d(ji)   = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice  
    294292            ! snow melting only = water into the ocean (then without snow precip) 
    295             wfx_snw_1d(ji)   = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
     293            wfx_snw_1d(ji)   = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
    296294 
    297295            ! updates available heat + thickness 
     
    308306      !---------------------- 
    309307      ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 
     308      ! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean) 
     309      ! clem comment: ice should also sublimate 
    310310      IF( lk_cpl ) THEN 
    311311         ! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice) 
     
    320320               &  ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) )  & 
    321321               &  * a_i_b(ji) * r1_rdtice 
    322             hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff ! diag only (to close heat budget) 
    323             ! heat used for sublimation (>0, W.m-2) 
    324             !!? hfx_tot_1d(ji) = hfx_tot_1d(ji) - zcoeff 
     322            hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff 
    325323            ! Mass flux by sublimation 
    326             wfx_sub_1d(ji) =  wfx_sub_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice ! diag only 
    327             wfx_snw_1d(ji) =  wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice 
     324            wfx_sub_1d(ji) =  wfx_sub_1d(ji) - rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice 
    328325            ! new snow thickness 
    329326            ht_s_b(ji)     =  MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) ) 
     
    388385            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice 
    389386 
    390             ! Total heat flux used in this process [W.m-2], < 
    391             hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 
     387            ! Total heat flux used in this process [W.m-2], > 
     388            hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 
    392389 
    393390            ! Contribution to mass flux 
    394             wfx_sum_1d(ji) =  wfx_sum_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
     391            wfx_sum_1d(ji) =  wfx_sum_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
    395392            
    396393            ! record which layers have disappeared (for bottom melting)  
     
    424421      ! Basal growth is driven by heat imbalance at the ice-ocean interface, 
    425422      ! between the inner conductive flux  (fc_bo_i), from the open water heat flux  
    426       ! (fhldb) and the turbulent ocean flux (fhtur).  
     423      ! (fhld) and the turbulent ocean flux (fhtur).  
    427424      ! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice  
    428425 
     
    490487            ! New ice growth 
    491488                                     
    492             zfmdt          = - rhoic * dh_i_bott(ji)                       ! Mass flux x time step (kg/m2, < 0) 
     489            zfmdt          = - rhoic * dh_i_bott(ji)             ! Mass flux x time step (kg/m2, < 0) 
     490 
     491            ztmelts        = - tmut * s_i_new(ji) + rtt          ! New ice melting point (K) 
     492             
     493            zt_i_new       = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i) 
     494             
     495            zEi            = cpic * ( zt_i_new - ztmelts ) &     ! Specific enthalpy of forming ice (J/kg, <0)       
     496               &               - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) )   & 
     497               &               + rcp  * ( ztmelts-rtt )           
     498             
     499            zEw            = rcp  * ( t_bo_b(ji) - rt0 )         ! Specific enthalpy of seawater (J/kg, < 0) 
     500             
     501            zdE            = zEi - zEw                           ! Specific enthalpy difference (J/kg, <0) 
    493502             
    494503            ! Contribution to heat flux to the ocean [W.m-2], >0   
    495504            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice 
    496             ! Total heat flux used in this process [W.m-2]   
    497             hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 
     505 
     506            ! Total heat flux used in this process [W.m-2], <0   
     507            hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 
    498508             
    499             ! Contribution to salt flux  () 
     509            ! Contribution to salt flux, <0 
    500510            sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice 
    501511 
    502             ! Contribution to mass flux 
    503             wfx_bog_1d(ji) =  wfx_bog_1d(ji) + rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice 
     512            ! Contribution to mass flux, <0 
     513            wfx_bog_1d(ji) =  wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice 
    504514 
    505515            ! clem debug: variation of enthalpy (J.m-2) 
     
    542552                  hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice 
    543553 
     554                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok) 
     555                  sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
     556                                     
     557                  ! Contribution to mass flux 
     558                  wfx_res_1d(ji) =  wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
     559 
    544560                  ! clem debug: variation of enthalpy (J.m-2) 
    545561                  dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)   
     
    573589                  ! Contribution to heat flux to the ocean [W.m-2], <0   
    574590                  hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice 
     591 
     592                  ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok) 
     593                  sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
     594                   
     595                  ! Total heat flux used in this process [W.m-2], >0   
     596                  hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 
     597                   
     598                  ! Contribution to mass flux 
     599                  wfx_bom_1d(ji) =  wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
    575600 
    576601                  ! clem debug: variation of enthalpy (J.m-2) 
     
    581606                  h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 
    582607               ENDIF 
    583  
    584                ! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok) 
    585                sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice 
    586  
    587                ! Total heat flux used in this process [W.m-2]   
    588                hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice 
    589  
    590                ! Contribution to mass flux 
    591                wfx_bom_1d(ji) =  wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice 
    592608            
    593609            ENDIF 
     
    602618      ! clem bug: I think this should be included above, so we would not have to  
    603619      !           track heat/salt/mass fluxes backwards 
    604       IF( jpl == 1 ) THEN 
    605          DO ji = kideb, kiut 
    606             IF(  zf_tt(ji)  >=  0._wp  ) THEN 
    607                zdh            = MAX( hmelt , dh_i_bott(ji) ) 
    608                zdvres         = zdh - dh_i_bott(ji) ! >=0 
    609                dh_i_bott(ji)  = zdh 
    610  
    611                ! excessive energy is sent to lateral ablation 
    612                zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) ) 
    613                zq_1cat(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0 
    614  
    615                ! correct salt and mass fluxes 
    616                sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 
    617                wfx_bom_1d(ji) = wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdvres * r1_rdtice 
    618             ENDIF 
    619          END DO 
    620       ENDIF 
     620!      IF( jpl == 1 ) THEN 
     621!         DO ji = kideb, kiut 
     622!            IF(  zf_tt(ji)  >=  0._wp  ) THEN 
     623!               zdh            = MAX( hmelt , dh_i_bott(ji) ) 
     624!               zdvres         = zdh - dh_i_bott(ji) ! >=0 
     625!               dh_i_bott(ji)  = zdh 
     626! 
     627!               ! excessive energy is sent to lateral ablation 
     628!               zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) ) 
     629!               zq_1cat(ji) =  zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0 
     630! 
     631!               ! correct salt and mass fluxes 
     632!               sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation 
     633!               wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdvres * r1_rdtice 
     634!            ENDIF 
     635!         END DO 
     636!      ENDIF 
    621637 
    622638      !------------------------------------------- 
     
    644660!         
    645661!         zq_rema(ji)     = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji)                ! update available heat (J.m-2) 
    646 !         ! Heat flux associated with snow melt 
    647 !         hfx_snw_1d(ji)  = hfx_snw_1d(ji) + zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (<0) 
    648662!         ! heat used to melt snow 
    649 !         hfx_tot_1d(ji)  = hfx_tot_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 
     663!         hfx_snw_1d(ji)  = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0) 
    650664!         ! Contribution to mass flux 
    651 !         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice 
     665!         wfx_snw_1d(ji)  =  wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice 
    652666!         ! clem debug: variation of enthalpy (J.m-2) 
    653667!         dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)   
     
    680694         ! new salinity difference stored (to be used in limthd_ent.F90) 
    681695         IF (  num_sal == 2  ) THEN 
    682             i_ice_switch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) ) 
     696            zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) ) 
    683697            ! salinity dif due to snow-ice formation 
    684             dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch      
     698            dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch      
    685699            ! salinity dif due to bottom growth  
    686700            IF (  zf_tt(ji)  < 0._wp ) THEN 
    687                dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch 
     701               dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch 
    688702            ENDIF 
    689703         ENDIF 
     
    704718         ! Contribution to mass flux 
    705719         ! All snow is thrown in the ocean, and seawater is taken to replace the volume 
    706          wfx_sni_1d(ji) = wfx_sni_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
    707          wfx_snw_1d(ji) = wfx_snw_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
     720         wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice 
     721         wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice 
    708722 
    709723         ! clem debug: variation of enthalpy (J.m-2) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dif.F90

    r4634 r4649  
    3232   PUBLIC   lim_thd_dif   ! called by lim_thd 
    3333 
    34    REAL(wp) ::   epsi10      = 1.e-10_wp    ! 
     34   REAL(wp) ::   epsi10 = 1.e-10_wp    ! 
    3535   !!---------------------------------------------------------------------- 
    3636   !! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2011) 
     
    9292      !!           (04-2007) Energy conservation tested by M. Vancoppenolle 
    9393      !!------------------------------------------------------------------ 
    94       INTEGER , INTENT (in) ::   kideb   ! Start point on which the  the computation is applied 
    95       INTEGER , INTENT (in) ::   kiut    ! End point on which the  the computation is applied 
    96       INTEGER , INTENT (in) ::   jl      ! Category number 
     94      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
     95      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
    9796 
    9897      !! * Local variables 
     
    103102      INTEGER ::   nconv       ! number of iterations in iterative procedure 
    104103      INTEGER ::   minnumeqmin, maxnumeqmax 
    105       INTEGER, DIMENSION(kiut) ::   numeqmin   ! reference number of top equation 
    106       INTEGER, DIMENSION(kiut) ::   numeqmax   ! reference number of bottom equation 
    107       INTEGER, DIMENSION(kiut) ::   isnow      ! switch for presence (1) or absence (0) of snow 
     104      INTEGER, POINTER, DIMENSION(:) ::   numeqmin   ! reference number of top equation 
     105      INTEGER, POINTER, DIMENSION(:) ::   numeqmax   ! reference number of bottom equation 
     106      INTEGER, POINTER, DIMENSION(:) ::   isnow      ! switch for presence (1) or absence (0) of snow 
    108107      REAL(wp) ::   zg1s      =  2._wp        ! for the tridiagonal system 
    109108      REAL(wp) ::   zg1       =  2._wp        ! 
     
    115114      REAL(wp) ::   ztmelt_i    ! ice melting temperature 
    116115      REAL(wp) ::   zerritmax   ! current maximal error on temperature  
    117       REAL(wp), DIMENSION(kiut) ::   ztfs        ! ice melting point 
    118       REAL(wp), DIMENSION(kiut) ::   ztsuold     ! old surface temperature (before the iterative procedure ) 
    119       REAL(wp), DIMENSION(kiut) ::   ztsuoldit   ! surface temperature at previous iteration 
    120       REAL(wp), DIMENSION(kiut) ::   zh_i        ! ice layer thickness 
    121       REAL(wp), DIMENSION(kiut) ::   zh_s        ! snow layer thickness 
    122       REAL(wp), DIMENSION(kiut) ::   zfsw        ! solar radiation absorbed at the surface 
    123       REAL(wp), DIMENSION(kiut) ::   zf          ! surface flux function 
    124       REAL(wp), DIMENSION(kiut) ::   dzf         ! derivative of the surface flux function 
    125       REAL(wp), DIMENSION(kiut) ::   zerrit      ! current error on temperature 
    126       REAL(wp), DIMENSION(kiut) ::   zdifcase    ! case of the equation resolution (1->4) 
    127       REAL(wp), DIMENSION(kiut) ::   zftrice     ! solar radiation transmitted through the ice 
    128       REAL(wp), DIMENSION(kiut) ::   zihic, zhsu 
    129       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztcond_i    ! Ice thermal conductivity 
    130       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradtr_i    ! Radiation transmitted through the ice 
    131       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zradab_i    ! Radiation absorbed in the ice 
    132       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zkappa_i    ! Kappa factor in the ice 
    133       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztiold      ! Old temperature in the ice 
    134       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zeta_i      ! Eta factor in the ice 
    135       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
    136       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   zspeche_i   ! Ice specific heat 
    137       REAL(wp), DIMENSION(kiut,0:nlay_i) ::   z_i         ! Vertical cotes of the layers in the ice 
    138       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradtr_s    ! Radiation transmited through the snow 
    139       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zradab_s    ! Radiation absorbed in the snow 
    140       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zkappa_s    ! Kappa factor in the snow 
    141       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   zeta_s       ! Eta factor in the snow 
    142       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztstemp      ! Temporary temperature in the snow to check the convergence 
    143       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   ztsold       ! Temporary temperature in the snow 
    144       REAL(wp), DIMENSION(kiut,0:nlay_s) ::   z_s          ! Vertical cotes of the layers in the snow 
    145       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindterm   ! Independent term 
    146       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zindtbis   ! temporary independent term 
    147       REAL(wp), DIMENSION(kiut,jkmax+2) ::   zdiagbis 
    148       REAL(wp), DIMENSION(kiut,jkmax+2,3) ::   ztrid   ! tridiagonal system terms 
    149       REAL(wp) ::   ztemp   ! local scalar  
     116      REAL(wp), POINTER, DIMENSION(:) ::   ztfs        ! ice melting point 
     117      REAL(wp), POINTER, DIMENSION(:) ::   ztsuold     ! old surface temperature (before the iterative procedure ) 
     118      REAL(wp), POINTER, DIMENSION(:) ::   ztsuoldit   ! surface temperature at previous iteration 
     119      REAL(wp), POINTER, DIMENSION(:) ::   zh_i        ! ice layer thickness 
     120      REAL(wp), POINTER, DIMENSION(:) ::   zh_s        ! snow layer thickness 
     121      REAL(wp), POINTER, DIMENSION(:) ::   zfsw        ! solar radiation absorbed at the surface 
     122      REAL(wp), POINTER, DIMENSION(:) ::   zf          ! surface flux function 
     123      REAL(wp), POINTER, DIMENSION(:) ::   dzf         ! derivative of the surface flux function 
     124      REAL(wp), POINTER, DIMENSION(:) ::   zerrit      ! current error on temperature 
     125      REAL(wp), POINTER, DIMENSION(:) ::   zdifcase    ! case of the equation resolution (1->4) 
     126      REAL(wp), POINTER, DIMENSION(:) ::   zftrice     ! solar radiation transmitted through the ice 
     127      REAL(wp), POINTER, DIMENSION(:) ::   zihic, zhsu 
     128      REAL(wp), POINTER, DIMENSION(:,:) ::   ztcond_i    ! Ice thermal conductivity 
     129      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_i    ! Radiation transmitted through the ice 
     130      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_i    ! Radiation absorbed in the ice 
     131      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_i    ! Kappa factor in the ice 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   ztiold      ! Old temperature in the ice 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_i      ! Eta factor in the ice 
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   ztitemp     ! Temporary temperature in the ice to check the convergence 
     135      REAL(wp), POINTER, DIMENSION(:,:) ::   zspeche_i   ! Ice specific heat 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   z_i         ! Vertical cotes of the layers in the ice 
     137      REAL(wp), POINTER, DIMENSION(:,:) ::   zradtr_s    ! Radiation transmited through the snow 
     138      REAL(wp), POINTER, DIMENSION(:,:) ::   zradab_s    ! Radiation absorbed in the snow 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   zkappa_s    ! Kappa factor in the snow 
     140      REAL(wp), POINTER, DIMENSION(:,:) ::   zeta_s       ! Eta factor in the snow 
     141      REAL(wp), POINTER, DIMENSION(:,:) ::   ztstemp      ! Temporary temperature in the snow to check the convergence 
     142      REAL(wp), POINTER, DIMENSION(:,:) ::   ztsold       ! Temporary temperature in the snow 
     143      REAL(wp), POINTER, DIMENSION(:,:) ::   z_s          ! Vertical cotes of the layers in the snow 
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   zindterm   ! Independent term 
     145      REAL(wp), POINTER, DIMENSION(:,:) ::   zindtbis   ! temporary independent term 
     146      REAL(wp), POINTER, DIMENSION(:,:) ::   zdiagbis 
     147      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ztrid   ! tridiagonal system terms 
    150148      !!------------------------------------------------------------------      
    151149      !  
     150      CALL wrk_alloc( jpij, numeqmin, numeqmax, isnow ) 
     151      CALL wrk_alloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 
     152      CALL wrk_alloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 
     153      CALL wrk_alloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart=0) 
     154      CALL wrk_alloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart=0) 
     155      CALL wrk_alloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis  ) 
     156      CALL wrk_alloc( jpij, jkmax+2, 3, ztrid ) 
     157 
    152158      !------------------------------------------------------------------------------! 
    153159      ! 1) Initialization                                                            ! 
     
    318324            DO layer = 1, nlay_i-1 
    319325               DO ji = kideb , kiut 
    320                   ztemp = t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2._wp * rtt 
    321                   ztcond_i(ji,layer) = rcdic + 0.0900_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
    322                      &                                   / MIN( -2.0_wp * epsi10, ztemp )   & 
    323                      &                       - 0.0055_wp * ztemp 
     326                  ztcond_i(ji,layer) = rcdic + 0.090_wp * ( s_i_b(ji,layer) + s_i_b(ji,layer+1) )   & 
     327                     &                                  / MIN(-2.0_wp * epsi10, t_i_b(ji,layer)+t_i_b(ji,layer+1) - 2.0_wp * rtt)   & 
     328                     &                       - 0.0055_wp* ( t_i_b(ji,layer) + t_i_b(ji,layer+1) - 2.0*rtt )   
    324329                  ztcond_i(ji,layer) = MAX( ztcond_i(ji,layer), zkimin ) 
    325330               END DO 
    326331            END DO 
    327332            DO ji = kideb , kiut 
    328                ztemp = t_bo_b(ji) - rtt 
    329                ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN( -epsi10, ztemp )   & 
    330                   &                        - 0.011_wp * ztemp   
     333               ztcond_i(ji,nlay_i) = rcdic + 0.090_wp * s_i_b(ji,nlay_i) / MIN(-epsi10,t_bo_b(ji)-rtt)   & 
     334                  &                        - 0.011_wp * ( t_bo_b(ji) - rtt )   
    331335               ztcond_i(ji,nlay_i) = MAX( ztcond_i(ji,nlay_i), zkimin ) 
    332336            END DO 
     
    395399         ! 
    396400         DO ji = kideb , kiut 
    397  
    398401            ! update of the non solar flux according to the update in T_su 
    399402            qns_ice_1d(ji) = qns_ice_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) 
     
    403406               + qns_ice_1d(ji)                  ! non solar total flux  
    404407            ! (LWup, LWdw, SH, LH) 
    405  
    406             ! heat flux used to change surface temperature 
    407             !hfx_tot_1d(ji) = hfx_tot_1d(ji) + dqns_ice_1d(ji) * ( t_su_b(ji) - ztsuoldit(ji) ) * a_i_b(ji) 
    408408         END DO 
    409409 
     
    721721      DO ji = kideb, kiut 
    722722         IF( t_su_b(ji) < rtt ) THEN  ! case T_su < 0degC 
    723             hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
     723            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( qns_ice_1d(ji) + qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
    724724         ELSE                                    ! case T_su = 0degC 
    725             hfx_tot_1d(ji) = hfx_tot_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
     725            hfx_dif_1d(ji) = hfx_dif_1d(ji) + ( fc_su(ji) + i0(ji) * qsr_ice_1d(ji) - zradtr_i(ji,nlay_i) - fc_bo_i(ji) ) * a_i_b(ji) 
    726726         ENDIF 
    727727      END DO 
    728728 
    729729      ! 
     730      CALL wrk_dealloc( jpij, numeqmin, numeqmax, isnow ) 
     731      CALL wrk_dealloc( jpij, ztfs, ztsuold, ztsuoldit, zh_i, zh_s, zfsw ) 
     732      CALL wrk_dealloc( jpij, zf, dzf, zerrit, zdifcase, zftrice, zihic, zhsu ) 
     733      CALL wrk_dealloc( jpij, nlay_i+1, ztcond_i, zradtr_i, zradab_i, zkappa_i, ztiold, zeta_i, ztitemp, z_i, zspeche_i, kjstart = 0 ) 
     734      CALL wrk_dealloc( jpij, nlay_s+1,           zradtr_s, zradab_s, zkappa_s, ztsold, zeta_s, ztstemp, z_s, kjstart = 0 ) 
     735      CALL wrk_dealloc( jpij, jkmax+2, zindterm, zindtbis, zdiagbis ) 
     736      CALL wrk_dealloc( jpij, jkmax+2, 3, ztrid ) 
     737 
    730738   END SUBROUTINE lim_thd_dif 
    731739 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_ent.F90

    r4634 r4649  
    4848CONTAINS 
    4949  
    50    SUBROUTINE lim_thd_ent( kideb, kiut, jl ) 
     50   SUBROUTINE lim_thd_ent( kideb, kiut, jl, qnew ) 
    5151      !!------------------------------------------------------------------- 
    5252      !!               ***   ROUTINE lim_thd_ent  *** 
     
    6262      !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 
    6363      !!            2) linear remapping on the new layers 
    64       !!            3) Ice salinity update + recover temperature from enthalpies 
     64      !! 
     65      !! ------------ cum0(0)                        ------------- cum1(0) 
     66      !!                                    NEW      ------------- 
     67      !! ------------ cum0(1)               ==>      ------------- 
     68      !!     ...                                     ------------- 
     69      !! ------------                                ------------- 
     70      !! ------------ cum0(nlay_i+2)                 ------------- cum1(nlay_i) 
     71      !! 
    6572      !! 
    6673      !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 
     
    6875      INTEGER , INTENT(in) ::   kideb, kiut   ! Start/End point on which the  the computation is applied 
    6976      INTEGER , INTENT(in) ::   jl            ! Thickness cateogry number 
     77 
     78      REAL(wp), INTENT(inout), DIMENSION(:,:) :: qnew          ! new enthlapies (remapped) 
    7079 
    7180      INTEGER  :: ji,ii,ij   !  dummy loop indices 
     
    7685      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum0, zh_cum0   ! old cumulative enthlapies and layers interfaces 
    7786      REAL(wp), POINTER, DIMENSION(:,:) :: zqh_cum1, zh_cum1   ! new cumulative enthlapies and layers interfaces 
     87      REAL(wp), POINTER, DIMENSION(:)   :: zhnew               ! new layers thicknesses 
    7888      !!------------------------------------------------------------------- 
    7989 
    8090      CALL wrk_alloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 
    8191      CALL wrk_alloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 ) 
     92      CALL wrk_alloc( jpij, zhnew ) 
    8293 
    8394      !-------------------------------------------------------------------------- 
     
    96107      !  2) Interpolation on the new layers 
    97108      !------------------------------------ 
     109      ! new layer thickesses 
     110      DO ji = kideb, kiut 
     111         zhnew(ji) = SUM( h_i_old(ji,0:nlay_i+1) ) / REAL( nlay_i )   
     112      ENDDO 
     113 
    98114      ! new layers interfaces 
    99115      zh_cum1(:,0:nlay_i) = 0._wp 
    100116      DO jk1 = 1, nlay_i 
    101117         DO ji = kideb, kiut 
    102             zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + ht_i_b(ji) / REAL( nlay_i ) 
     118            zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 
    103119         ENDDO 
    104120      ENDDO 
     
    123139      DO jk1 = 1, nlay_i 
    124140         DO ji = kideb, kiut 
    125             zswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) + epsi20 ) )  
    126             q_i_b(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) * REAL( nlay_i ) / MAX( ht_i_b(ji), epsi20 ) 
     141            zswitch      =  1._wp - MAX( 0._wp , SIGN( 1._wp , - zhnew(ji) + epsi10 ) )  
     142            qnew(ji,jk1) = zswitch * ( zqh_cum1(ji,jk1) - zqh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi10 ) 
    127143         ENDDO 
    128144      ENDDO 
    129  
    130       !--------------------------------------------------------- 
    131       !  3) Update ice salinity and recover ice temperature 
    132       !--------------------------------------------------------- 
    133       ! Update salinity (basal entrapment, snow ice formation) 
    134       DO ji = kideb, kiut 
    135          sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 
    136       END DO !ji 
    137  
    138       ! Recover ice temperature 
    139       DO jk1 = 1, nlay_i 
    140          DO ji = kideb, kiut 
    141             ztmelts       =  -tmut * s_i_b(ji,jk1) + rtt 
    142             ! Conversion q(S,T) -> T (second order equation) 
    143             zaaa          =  cpic 
    144             zbbb          =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk1) / rhoic - lfus 
    145             zccc          =  lfus * ( ztmelts - rtt ) 
    146             zdiscrim      =  SQRT( MAX( zbbb * zbbb - 4._wp * zaaa * zccc, 0._wp ) ) 
    147             t_i_b(ji,jk1) =  rtt - ( zbbb + zdiscrim ) / ( 2._wp * zaaa ) 
    148              
    149             ! mask temperature 
    150             zswitch       =  1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )  
    151             t_i_b(ji,jk1) =  zswitch * t_i_b(ji,jk1) + ( 1._wp - zswitch ) * rtt 
    152          END DO  
    153       END DO  
    154  
    155145      ! 
    156146      CALL wrk_dealloc( jpij, nlay_i+3, zqh_cum0, zh_cum0, kjstart = 0 ) 
    157147      CALL wrk_dealloc( jpij, nlay_i+1, zqh_cum1, zh_cum1, kjstart = 0 ) 
     148      CALL wrk_dealloc( jpij, zhnew ) 
    158149      ! 
    159150   END SUBROUTINE lim_thd_ent 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_lac.F90

    r4634 r4649  
    3030   USE wrk_nemo       ! work arrays 
    3131   USE lib_fortran    ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)   
     32   USE limthd_ent 
    3233 
    3334   IMPLICIT NONE 
     
    8586      REAL(wp) ::   zEw          ! seawater specific enthalpy (J/kg) 
    8687      REAL(wp) ::   zfmdt        ! mass flux x time step (kg/m2, >0 towards ocean) 
    87         
    88       INTEGER , POINTER, DIMENSION(:) ::   zcatac      ! indexes of categories where new ice grows 
     88      
     89      REAL(wp) ::   zv_newfra 
     90   
     91      INTEGER , POINTER, DIMENSION(:) ::   jcat      ! indexes of categories where new ice grows 
    8992      REAL(wp), POINTER, DIMENSION(:) ::   zswinew     ! switch for new ice or not 
    9093 
     
    97100      REAL(wp), POINTER, DIMENSION(:) ::   zdv_res     ! residual volume in case of excessive heat budget 
    98101      REAL(wp), POINTER, DIMENSION(:) ::   zda_res     ! residual area in case of excessive heat budget 
    99       REAL(wp), POINTER, DIMENSION(:) ::   zat_i_ac    ! total ice fraction     
     102      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_1d    ! total ice fraction     
    100103      REAL(wp), POINTER, DIMENSION(:) ::   zat_i_lev   ! total ice fraction for level ice only (type 1)    
    101104      REAL(wp), POINTER, DIMENSION(:) ::   zv_frazb   ! accretion of frazil ice at the ice bottom 
    102       REAL(wp), POINTER, DIMENSION(:) ::   zvrel_ac    ! relative ice / frazil velocity (1D vector) 
     105      REAL(wp), POINTER, DIMENSION(:) ::   zvrel_1d    ! relative ice / frazil velocity (1D vector) 
    103106 
    104107      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_old      ! old volume of ice in category jl 
    105108      REAL(wp), POINTER, DIMENSION(:,:) ::   za_old      ! old area of ice in category jl 
    106       REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_ac     ! 1-D version of a_i 
    107       REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_ac     ! 1-D version of v_i 
    108       REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_ac    ! 1-D version of oa_i 
    109       REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_ac   ! 1-D version of smv_i 
    110  
    111       REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_ac   !: 1-D version of e_i 
    112  
    113       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zqm0      ! old layer-system heat content 
    114       REAL(wp), POINTER, DIMENSION(:,:,:) ::   zthick0   ! old ice thickness 
    115  
    116       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_i_init, vt_i_final   ! ice volume summed over categories 
    117       REAL(wp), POINTER, DIMENSION(:,:) ::   vt_s_init, vt_s_final   !  snow volume summed over categories 
    118       REAL(wp), POINTER, DIMENSION(:,:) ::   et_i_init, et_i_final   !  ice energy summed over categories 
    119       REAL(wp), POINTER, DIMENSION(:,:) ::   et_s_init               !  snow energy summed over categories 
     109      REAL(wp), POINTER, DIMENSION(:,:) ::   za_i_1d     ! 1-D version of a_i 
     110      REAL(wp), POINTER, DIMENSION(:,:) ::   zv_i_1d     ! 1-D version of v_i 
     111      REAL(wp), POINTER, DIMENSION(:,:) ::   zoa_i_1d    ! 1-D version of oa_i 
     112      REAL(wp), POINTER, DIMENSION(:,:) ::   zsmv_i_1d   ! 1-D version of smv_i 
     113 
     114      REAL(wp), POINTER, DIMENSION(:,:,:) ::   ze_i_1d   !: 1-D version of e_i 
     115 
    120116      REAL(wp), POINTER, DIMENSION(:,:) ::   zvrel                   ! relative ice / frazil velocity 
    121117      !!-----------------------------------------------------------------------! 
    122118 
    123       CALL wrk_alloc( jpij, zcatac )   ! integer 
     119      CALL wrk_alloc( jpij, jcat )   ! integer 
    124120      CALL wrk_alloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    125       CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zv_frazb, zvrel_ac ) 
    126       CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 
    127       CALL wrk_alloc( jpij,jkmax,jpl, ze_i_ac ) 
    128       CALL wrk_alloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 
    129       CALL wrk_alloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 
    130  
    131       et_i_init(:,:) = 0._wp 
    132       et_s_init(:,:) = 0._wp 
    133       vt_i_init(:,:) = 0._wp 
    134       vt_s_init(:,:) = 0._wp 
    135  
    136       !------------------------------------------------------------------------------! 
    137       ! 1) Conservation check and changes in each ice category 
    138       !------------------------------------------------------------------------------! 
    139       IF( con_i ) THEN 
    140          CALL lim_column_sum        ( jpl, v_i          , vt_i_init) 
    141          CALL lim_column_sum        ( jpl, v_s          , vt_s_init) 
    142          CALL lim_column_sum_energy ( jpl, nlay_i , e_i , et_i_init) 
    143          CALL lim_column_sum        ( jpl, e_s(:,:,1,:) , et_s_init) 
    144       ENDIF 
     121      CALL wrk_alloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 
     122      CALL wrk_alloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
     123      CALL wrk_alloc( jpij,jkmax,jpl, ze_i_1d ) 
     124      CALL wrk_alloc( jpi,jpj, zvrel ) 
    145125 
    146126      !------------------------------------------------------------------------------| 
     
    204184                     &          +   vtau_ice(ji  ,jj  ) * tmv(ji  ,jj  ) ) * 0.5_wp 
    205185                  ! Square root of wind stress 
    206                   ztenagm       =  SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 
     186                  ztenagm       =  SQRT( SQRT( ztaux**2 + ztauy**2 ) ) 
    207187 
    208188                  !--------------------- 
    209189                  ! Frazil ice velocity 
    210190                  !--------------------- 
    211                   zvfrx         = zgamafr * zsqcd * ztaux / MAX(ztenagm,epsi10) 
    212                   zvfry         = zgamafr * zsqcd * ztauy / MAX(ztenagm,epsi10) 
     191                  zindb = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 
     192                  zvfrx = zindb * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 
     193                  zvfry = zindb * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 
    213194 
    214195                  !------------------- 
     
    305286      IF ( nbpac > 0 ) THEN 
    306287 
    307          CALL tab_2d_1d( nbpac, zat_i_ac  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) ) 
     288         CALL tab_2d_1d( nbpac, zat_i_1d  (1:nbpac)     , at_i         , jpi, jpj, npac(1:nbpac) ) 
    308289         DO jl = 1, jpl 
    309             CALL tab_2d_1d( nbpac, za_i_ac  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    310             CALL tab_2d_1d( nbpac, zv_i_ac  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    311             CALL tab_2d_1d( nbpac, zoa_i_ac (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    312             CALL tab_2d_1d( nbpac, zsmv_i_ac(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     290            CALL tab_2d_1d( nbpac, za_i_1d  (1:nbpac,jl), a_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     291            CALL tab_2d_1d( nbpac, zv_i_1d  (1:nbpac,jl), v_i  (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     292            CALL tab_2d_1d( nbpac, zoa_i_1d (1:nbpac,jl), oa_i (:,:,jl), jpi, jpj, npac(1:nbpac) ) 
     293            CALL tab_2d_1d( nbpac, zsmv_i_1d(1:nbpac,jl), smv_i(:,:,jl), jpi, jpj, npac(1:nbpac) ) 
    313294            DO jk = 1, nlay_i 
    314                CALL tab_2d_1d( nbpac, ze_i_ac(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
     295               CALL tab_2d_1d( nbpac, ze_i_1d(1:nbpac,jk,jl), e_i(:,:,jk,jl) , jpi, jpj, npac(1:nbpac) ) 
    315296            END DO ! jk 
    316297         END DO ! jl 
     
    322303         CALL tab_2d_1d( nbpac, wfx_opw_1d(1:nbpac)     , wfx_opw, jpi, jpj, npac(1:nbpac) ) 
    323304         CALL tab_2d_1d( nbpac, hicol_b   (1:nbpac)     , hicol  , jpi, jpj, npac(1:nbpac) ) 
    324          CALL tab_2d_1d( nbpac, zvrel_ac  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
     305         CALL tab_2d_1d( nbpac, zvrel_1d  (1:nbpac)     , zvrel  , jpi, jpj, npac(1:nbpac) ) 
    325306 
    326307         CALL tab_2d_1d( nbpac, hfx_thd_1d(1:nbpac)     , hfx_thd, jpi, jpj, npac(1:nbpac) ) 
    327          CALL tab_2d_1d( nbpac, hfx_tot_1d(1:nbpac)     , hfx_tot, jpi, jpj, npac(1:nbpac) ) 
     308         CALL tab_2d_1d( nbpac, hfx_opw_1d(1:nbpac)     , hfx_opw, jpi, jpj, npac(1:nbpac) ) 
    328309 
    329310         !------------------------------------------------------------------------------! 
     
    334315         ! Keep old ice areas and volume in memory 
    335316         !----------------------------------------- 
    336          zv_old(:,:) = zv_i_ac(:,:)  
    337          za_old(:,:) = za_i_ac(:,:) 
     317         zv_old(:,:) = zv_i_1d(:,:)  
     318         za_old(:,:) = za_i_1d(:,:) 
    338319 
    339320         !---------------------- 
     
    343324            zh_newice(ji) = hiccrit(1) 
    344325         END DO 
    345          IF( fraz_swi == 1.0 )   zh_newice(:) = hicol_b(:) 
     326         IF( fraz_swi == 1.0 ) zh_newice(:) = hicol_b(:) 
    346327 
    347328         !---------------------- 
     
    370351            ztmelts       = - tmut * zs_newice(ji) + rtt                  ! Melting point (K) 
    371352            ze_newice(ji) =   rhoic * (  cpic * ( ztmelts - t_bo_b(ji) )                             & 
    372                &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / ( t_bo_b(ji) - rtt ) )   & 
     353               &                       + lfus * ( 1.0 - ( ztmelts - rtt ) / MIN( t_bo_b(ji) - rtt, -epsi10 ) )   & 
    373354               &                       - rcp  *         ( ztmelts - rtt )  ) 
    374             ! MV HC 2014 comment I dont see why this line below is here... ? 
    375             ! This implies that ze_newice gets to rhoic*Lfus if it was negative, but this should never happen 
    376             ze_newice(ji) =   MAX( ze_newice(ji) , 0._wp )    & 
    377                &          +   MAX(  0.0 , SIGN( 1.0 , - ze_newice(ji) )  ) * rhoic * lfus 
    378355         END DO ! ji 
    379356         !---------------- 
     
    405382            hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * zEw * r1_rdtice 
    406383            ! Total heat flux used in this process [W.m-2]   
    407             hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * zdE * r1_rdtice 
     384            hfx_opw_1d(ji) = hfx_opw_1d(ji) - zfmdt * zdE * r1_rdtice 
    408385            ! mass flux 
    409             wfx_opw_1d(ji) = wfx_opw_1d(ji) + zv_newice(ji) * rhoic * r1_rdtice 
     386            wfx_opw_1d(ji) = wfx_opw_1d(ji) - zv_newice(ji) * rhoic * r1_rdtice 
    410387            ! salt flux 
    411388            sfx_opw_1d(ji) = sfx_opw_1d(ji) - zv_newice(ji) * rhoic * zs_newice(ji) * r1_rdtice 
    412389 
    413390            ! A fraction zfrazb of frazil ice is accreted at the ice bottom 
    414             zfrazb        = ( TANH ( Cfrazb * ( zvrel_ac(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
     391            zinda         = 1._wp - MAX( 0._wp, SIGN( 1._wp , - zat_i_1d(ji) ) ) 
     392            zfrazb        = zinda * ( TANH ( Cfrazb * ( zvrel_1d(ji) - vfrazb ) ) + 1.0 ) * 0.5 * maxfrazb 
    415393            zv_frazb(ji)  =         zfrazb   * zv_newice(ji) 
    416394            zv_newice(ji) = ( 1.0 - zfrazb ) * zv_newice(ji) 
    417395         END DO 
    418396 
    419          !------------------------------------ 
    420          ! Diags for energy conservation test 
    421          !------------------------------------ 
    422          DO ji = 1, nbpac 
    423             ii = MOD( npac(ji) - 1 , jpi ) + 1 
    424             ij =    ( npac(ji) - 1 ) / jpi + 1 
    425             ! 
    426             zde = ze_newice(ji) / unit_fac * area(ii,ij) * zv_newice(ji) 
    427             !zde = ze_newice(ji) * area(ii,ij) * zv_newice(ji) 
    428             ! 
    429             ! clem: change that? 
    430             vt_i_init(ii,ij) = vt_i_init(ii,ij) + zv_newice(ji)             ! volume 
    431             et_i_init(ii,ij) = et_i_init(ii,ij) + zde                       ! Energy 
    432  
    433          END DO 
    434397 
    435398         !----------------- 
     
    438401         DO ji = 1, nbpac 
    439402            za_newice(ji) = zv_newice(ji) / zh_newice(ji) 
    440          END DO !ji 
     403         END DO 
    441404 
    442405         !------------------------------------------------------------------------------! 
     
    450413         ! we keep the excessive volume in memory and attribute it later to bottom accretion 
    451414         DO ji = 1, nbpac 
    452             IF ( za_newice(ji) >  ( amax - zat_i_ac(ji) ) ) THEN 
    453                zda_res(ji)   = za_newice(ji) - ( amax - zat_i_ac(ji) ) 
     415            IF ( za_newice(ji) >  ( amax - zat_i_1d(ji) ) ) THEN 
     416               zda_res(ji)   = za_newice(ji) - ( amax - zat_i_1d(ji) ) 
    454417               zdv_res(ji)   = zda_res  (ji) * zh_newice(ji)  
    455418               za_newice(ji) = za_newice(ji) - zda_res  (ji) 
     
    464427         ! Laterally redistribute new ice volume and area 
    465428         !------------------------------------------------ 
    466          zat_i_ac(:) = 0._wp 
     429         zat_i_1d(:) = 0._wp 
    467430         DO jl = 1, jpl 
    468431            DO ji = 1, nbpac 
    469                IF(  hi_max   (jl-1)  <   zh_newice(ji)   .AND.   & 
    470                   & zh_newice(ji)    <=  hi_max   (jl)         ) THEN 
    471                   za_i_ac (ji,jl) = za_i_ac (ji,jl) + za_newice(ji) 
    472                   zv_i_ac (ji,jl) = zv_i_ac (ji,jl) + zv_newice(ji) 
    473                   zcatac  (ji)    = jl 
     432               IF( zh_newice(ji) > hi_max(jl-1) .AND. zh_newice(ji) <= hi_max(jl) ) THEN 
     433                  za_i_1d (ji,jl) = za_i_1d (ji,jl) + za_newice(ji) 
     434                  zv_i_1d (ji,jl) = zv_i_1d (ji,jl) + zv_newice(ji) 
     435                  jcat  (ji)    = jl 
    474436               ENDIF 
    475                zat_i_ac(ji)    = zat_i_ac(ji)    + za_i_ac  (ji,jl) 
     437               zat_i_1d(ji) = zat_i_1d(ji) + za_i_1d  (ji,jl) 
    476438            END DO 
    477439         END DO 
     
    481443         !---------------------------------- 
    482444         DO ji = 1, nbpac 
    483             jl = zcatac(ji)                                                           ! categroy in which new ice is put 
    484             zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) + epsi10 ) )   ! 0 if old ice 
     445            jl = jcat(ji)                                                  ! categroy in which new ice is put 
     446            zswinew  (ji) = MAX( 0._wp , SIGN( 1._wp , - za_old(ji,jl) ) )   ! 0 if old ice 
    485447         END DO 
    486448 
    487449         DO jk = 1, nlay_i 
    488450            DO ji = 1, nbpac 
    489                jl = zcatac(ji) 
    490                ze_i_ac(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji)                & 
    491                   &      + ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_ac(ji,jk,jl) * zv_old(ji,jl) ) / zv_i_ac(ji,jl) 
     451               jl = jcat(ji) 
     452               zinda = MAX( 0._wp, SIGN( 1._wp , zv_i_1d(ji,jl) - epsi20 ) ) 
     453               ze_i_1d(ji,jk,jl) = zswinew(ji)   *   ze_newice(ji) +                                                      & 
     454                  &        ( 1.0 - zswinew(ji) ) * ( ze_newice(ji) * zv_newice(ji) + ze_i_1d(ji,jk,jl) * zv_old(ji,jl) )  & 
     455                  &        * zinda / MAX( zv_i_1d(ji,jl), epsi20 ) 
    492456            END DO 
    493457         END DO 
     
    496460         ! Add excessive volume of new ice at the bottom 
    497461         !----------------------------------------------- 
    498          jm = 1 
    499  
    500          ! ---  Redistributing energy on the new grid (energy is equally distributed in every layer) --- ! 
    501 !         DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    502 !            DO jk = 1, nlay_i 
    503 !               DO ji = 1, nbpac 
    504 !                  ze_i_ac(ji,jk,jl) = ( ze_i_ac(ji,jk,jl) * zv_i_ac(ji,jl) + ze_newice(ji) * ( zdv_res(ji) + zv_frazb(ji) ) ) / & 
    505 !                     &                ( zv_i_ac(ji,jl) + ( zdv_res(ji) + zv_frazb(ji) ) )  
    506 !               END DO 
    507 !            END DO 
    508 !         END DO 
    509  
    510          ! --- Redistributing energy on the new grid (energy is sent to the bottom) PART 1 --- ! 
    511          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
     462         DO jl = 1, jpl 
     463            h_i_old (1:nbpac,0:nlay_i+1) = 0._wp 
     464            qh_i_old(1:nbpac,0:nlay_i+1) = 0._wp 
     465 
    512466            DO jk = 1, nlay_i 
    513467               DO ji = 1, nbpac 
    514                   zthick0(ji,jk,jl) =  zv_i_ac(ji,jl) / REAL( nlay_i ) 
    515                   zqm0   (ji,jk,jl) =  ze_i_ac(ji,jk,jl) * zthick0(ji,jk,jl) 
     468                  h_i_old (ji,jk) = zv_i_1d(ji,jl) / REAL( nlay_i ) 
     469                  qh_i_old(ji,jk) = ze_i_1d(ji,jk,jl) * h_i_old(ji,jk) 
    516470               END DO 
    517471            END DO 
    518          END DO 
    519          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
     472 
    520473            DO ji = 1, nbpac 
    521                zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_ac(ji) - epsi10 ) ) 
    522                zthick0(ji,nlay_i+1,jl) =  zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_ac(ji,jl) / MAX( zat_i_ac(ji) , epsi10 ) 
    523                zqm0   (ji,nlay_i+1,jl) =  ze_newice(ji) * zthick0(ji,nlay_i+1,jl) 
    524             END DO ! ji 
    525          END DO ! jl 
    526  
    527          ze_i_ac(:,:,:) = 0._wp 
    528          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    529             DO jk = 1, nlay_i 
    530                DO layer = 1, nlay_i + 1 
    531                   DO ji = 1, nbpac 
    532                       zindb   =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zthick0(ji,layer,jl) + epsi10 ) )  
    533                       zweight = zindb * MAX( 0._wp,                                                                                              & 
    534                          &      MIN( zv_i_ac(ji,jl) * REAL( layer     ), ( zv_i_ac(ji,jl) +  zthick0(ji,nlay_i+1,jl) ) * REAL( jk     ) )     & 
    535                          &    - MAX( zv_i_ac(ji,jl) * REAL( layer - 1 ), ( zv_i_ac(ji,jl) +  zthick0(ji,nlay_i+1,jl) ) * REAL( jk - 1 ) ) )   & 
    536                          &  / ( REAL( nlay_i ) * MAX( zthick0(ji,layer,jl), epsi10 ) ) 
    537                       ze_i_ac(ji,jk,jl) = ze_i_ac(ji,jk,jl) + zweight * zqm0(ji,layer,jl)   
    538                   END DO  
    539                END DO  
    540             END DO  
    541          END DO 
    542  
    543          ! --- new volumes and layer thickness --- 
    544          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    545             DO ji = 1, nbpac 
    546                zinda = MAX( 0._wp, SIGN( 1._wp , zat_i_ac(ji) - epsi10 ) ) 
    547                zv_i_ac(ji,jl) = zv_i_ac(ji,jl) + zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_ac(ji,jl) / MAX( zat_i_ac(ji) , epsi10 ) 
    548             END DO 
    549          END DO 
    550  
    551          ! --- Redistributing energy on the new grid (energy is sent to the bottom) PART 2 --- ! 
    552          DO jl = ice_cat_bounds(jm,1), ice_cat_bounds(jm,2) 
    553             DO jk = 1, nlay_i 
    554                DO ji = 1, nbpac 
    555                   zindb =  1._wp -  MAX( 0._wp , SIGN( 1._wp , - zv_i_ac(ji,jl) + epsi10 ) )  
    556                   ze_i_ac(ji,jk,jl) = zindb * ze_i_ac(ji,jk,jl) / MAX( zv_i_ac(ji,jl), epsi10 ) * REAL( nlay_i ) 
    557                END DO 
    558             END DO 
    559          END DO 
    560  
     474               zinda          = MAX( 0._wp, SIGN( 1._wp , zat_i_1d(ji) - epsi20 ) ) 
     475               zv_newfra      = zinda * ( zdv_res(ji) + zv_frazb(ji) ) * za_i_1d(ji,jl) / MAX( zat_i_1d(ji) , epsi20 ) 
     476               za_i_1d(ji,jl) = zinda * za_i_1d(ji,jl)                
     477 
     478               h_i_old (ji,nlay_i+1) = zv_newfra 
     479               qh_i_old(ji,nlay_i+1) = ze_newice(ji) * zv_newfra 
     480 
     481               zv_i_1d(ji,jl) = zv_i_1d(ji,jl) + zv_newfra 
     482            ENDDO 
     483 
     484            ! --- Ice enthalpy remapping --- ! 
     485            CALL lim_thd_ent( 1, nbpac, jl, ze_i_1d(1:nbpac,:,jl) )  
     486 
     487         ENDDO 
    561488 
    562489         !------------ 
     
    565492         DO jl = 1, jpl 
    566493            DO ji = 1, nbpac 
    567                zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_ac(ji,jl) + epsi10 ) )  ! 0 if no ice and 1 if yes 
    568                zoa_i_ac(ji,jl)  = za_old(ji,jl) * zoa_i_ac(ji,jl) / MAX( za_i_ac(ji,jl) , epsi10 ) * zindb    
     494               zindb = 1._wp - MAX( 0._wp , SIGN( 1._wp , - za_i_1d(ji,jl) + epsi20 ) )  ! 0 if no ice and 1 if yes 
     495               zoa_i_1d(ji,jl)  = za_old(ji,jl) * zoa_i_1d(ji,jl) / MAX( za_i_1d(ji,jl) , epsi20 ) * zindb    
    569496            END DO  
    570497         END DO    
     
    576503            DO jl = 1, jpl 
    577504               DO ji = 1, nbpac 
    578                   zdv   = zv_i_ac(ji,jl) - zv_old(ji,jl) 
    579                   zsmv_i_ac(ji,jl) = zsmv_i_ac(ji,jl) + zdv * zs_newice(ji) 
     505                  zdv   = zv_i_1d(ji,jl) - zv_old(ji,jl) 
     506                  zsmv_i_1d(ji,jl) = zsmv_i_1d(ji,jl) + zdv * zs_newice(ji) 
    580507               END DO 
    581508            END DO    
     
    586513         !------------------------------------------------------------------------------! 
    587514         DO jl = 1, jpl 
    588             CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_ac (1:nbpac,jl), jpi, jpj ) 
    589             CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_ac (1:nbpac,jl), jpi, jpj ) 
    590             CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_ac(1:nbpac,jl), jpi, jpj ) 
    591             CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_ac(1:nbpac,jl) , jpi, jpj ) 
     515            CALL tab_1d_2d( nbpac, a_i (:,:,jl), npac(1:nbpac), za_i_1d (1:nbpac,jl), jpi, jpj ) 
     516            CALL tab_1d_2d( nbpac, v_i (:,:,jl), npac(1:nbpac), zv_i_1d (1:nbpac,jl), jpi, jpj ) 
     517            CALL tab_1d_2d( nbpac, oa_i(:,:,jl), npac(1:nbpac), zoa_i_1d(1:nbpac,jl), jpi, jpj ) 
     518            CALL tab_1d_2d( nbpac, smv_i (:,:,jl), npac(1:nbpac), zsmv_i_1d(1:nbpac,jl) , jpi, jpj ) 
    592519            DO jk = 1, nlay_i 
    593                CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_ac(1:nbpac,jk,jl), jpi, jpj ) 
     520               CALL tab_1d_2d( nbpac, e_i(:,:,jk,jl), npac(1:nbpac), ze_i_1d(1:nbpac,jk,jl), jpi, jpj ) 
    594521            END DO 
    595522         END DO 
     
    599526 
    600527         CALL tab_1d_2d( nbpac, hfx_thd, npac(1:nbpac), hfx_thd_1d(1:nbpac), jpi, jpj ) 
    601          CALL tab_1d_2d( nbpac, hfx_tot, npac(1:nbpac), hfx_tot_1d(1:nbpac), jpi, jpj ) 
     528         CALL tab_1d_2d( nbpac, hfx_opw, npac(1:nbpac), hfx_opw_1d(1:nbpac), jpi, jpj ) 
    602529         ! 
    603530      ENDIF ! nbpac > 0 
     
    612539                  ! heat content in Joules 
    613540                  e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) * unit_fac )  
    614                   !e_i(ji,jj,jk,jl) = e_i(ji,jj,jk,jl) * area(ji,jj) * v_i(ji,jj,jl) / ( REAL( nlay_i ) )  
    615541               END DO 
    616542            END DO 
     
    618544      END DO 
    619545 
    620       !------------------------------------------------------------------------------| 
    621       ! 10) Conservation check and changes in each ice category 
    622       !------------------------------------------------------------------------------| 
    623       IF( con_i ) THEN  
    624          CALL lim_column_sum (jpl,   v_i, vt_i_final) 
    625          fieldid = 'v_i, limthd_lac' 
    626          CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)  
    627          ! 
    628          CALL lim_column_sum_energy(jpl, nlay_i, e_i, et_i_final) 
    629          fieldid = 'e_i, limthd_lac' 
    630          CALL lim_cons_check (et_i_final, et_i_final, 1.0e-3, fieldid)  
    631          ! 
    632          CALL lim_column_sum (jpl,   v_s, vt_s_final) 
    633          fieldid = 'v_s, limthd_lac' 
    634          CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)  
    635          ! 
    636          !     CALL lim_column_sum (jpl,   e_s(:,:,1,:) , et_s_init) 
    637          !     fieldid = 'e_s, limthd_lac' 
    638          !     CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)  
    639          IF( ln_nicep ) THEN 
    640             DO ji = mi0(jiindx), mi1(jiindx) 
    641                DO jj = mj0(jjindx), mj1(jjindx) 
    642                   WRITE(numout,*) ' vt_i_init : ', vt_i_init (ji,jj) 
    643                   WRITE(numout,*) ' vt_i_final: ', vt_i_final(ji,jj) 
    644                   WRITE(numout,*) ' et_i_init : ', et_i_init (ji,jj) 
    645                   WRITE(numout,*) ' et_i_final: ', et_i_final(ji,jj) 
    646                END DO 
    647             END DO 
    648          ENDIF 
    649          ! 
    650       ENDIF 
    651546      ! 
    652       CALL wrk_dealloc( jpij, zcatac )   ! integer 
     547      CALL wrk_dealloc( jpij, jcat )   ! integer 
    653548      CALL wrk_dealloc( jpij, zswinew, zv_newice, za_newice, zh_newice, ze_newice, zs_newice, zo_newice ) 
    654       CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_ac, zat_i_lev, zv_frazb, zvrel_ac ) 
    655       CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_ac, zv_i_ac, zoa_i_ac, zsmv_i_ac ) 
    656       CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_ac ) 
    657       CALL wrk_dealloc( jpij,jkmax+1,jpl, zqm0, zthick0 ) 
    658       CALL wrk_dealloc( jpi,jpj, vt_i_init, vt_i_final, vt_s_init, vt_s_final, et_i_init, et_i_final, et_s_init, zvrel ) 
     549      CALL wrk_dealloc( jpij, zdv_res, zda_res, zat_i_1d, zat_i_lev, zv_frazb, zvrel_1d ) 
     550      CALL wrk_dealloc( jpij,jpl, zv_old, za_old, za_i_1d, zv_i_1d, zoa_i_1d, zsmv_i_1d ) 
     551      CALL wrk_dealloc( jpij,jkmax,jpl, ze_i_1d ) 
     552      CALL wrk_dealloc( jpi,jpj, zvrel ) 
    659553      ! 
    660554   END SUBROUTINE lim_thd_lac 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_sal.F90

    r4634 r4649  
    5353      ! 
    5454      INTEGER  ::   ji, jk     ! dummy loop indices  
    55       REAL(wp) ::   iflush, igravdr, ztmelts   ! local scalars 
    56       REAL(wp) ::   zaaa, zbbb, zccc, zdiscrim   ! local scalars 
     55      REAL(wp) ::   iflush, igravdr   ! local scalars 
    5756      !!--------------------------------------------------------------------- 
    5857 
     58      !--------------------------------------------------------- 
     59      !  0) Update ice salinity from snow-ice and bottom growth 
     60      !--------------------------------------------------------- 
     61      DO ji = kideb, kiut 
     62         sm_i_b(ji) = sm_i_b(ji) + dsm_i_se_1d(ji) + dsm_i_si_1d(ji) 
     63      END DO 
     64  
    5965      !------------------------------------------------------------------------------| 
    6066      ! 1) Constant salinity, constant in time                                       | 
     
    7177      !  Module 2 : Constant salinity varying in time                                | 
    7278      !------------------------------------------------------------------------------| 
    73  
    7479      IF(  num_sal == 2  ) THEN 
    7580 
     
    7883            ! Switches  
    7984            !---------- 
    80             iflush       =        MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt )        )    ! =1 if summer  
    81             igravdr      =        MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo 
     85            iflush  = MAX( 0._wp , SIGN( 1._wp , t_su_b(ji) - rtt )        )    ! =1 if summer  
     86            igravdr = MAX( 0._wp , SIGN( 1._wp , t_bo_b(ji) - t_su_b(ji) ) )    ! =1 if t_su < t_bo 
    8287 
    8388            !--------------------- 
    8489            ! Salinity tendencies 
    8590            !--------------------- 
    86             !                                
    87             dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  ! drainage by gravity 
    88             dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice  ! drainage by flushing  
    89             ! 
     91            ! drainage by gravity drainage 
     92            dsm_i_gd_1d(ji) = - igravdr * MAX( sm_i_b(ji) - sal_G , 0._wp ) / time_G * rdt_ice  
     93            ! drainage by flushing   
     94            dsm_i_fl_1d(ji) = - iflush  * MAX( sm_i_b(ji) - sal_F , 0._wp ) / time_F * rdt_ice 
     95 
    9096            !----------------- 
    9197            ! Update salinity    
     
    104110         ! Salinity profile 
    105111         CALL lim_var_salprof1d( kideb, kiut ) 
    106  
    107  
    108          ! Only necessary for conservation check since salinity is modified 
    109          !-------------------- 
    110          ! Temperature update 
    111          !-------------------- 
    112          DO jk = 1, nlay_i 
    113             DO ji = kideb, kiut 
    114                ztmelts    =  -tmut*s_i_b(ji,jk) + rtt 
    115                !Conversion q(S,T) -> T (second order equation) 
    116                zaaa         =  cpic 
    117                zbbb         =  ( rcp - cpic ) * ( ztmelts - rtt ) + q_i_b(ji,jk) / rhoic - lfus 
    118                zccc         =  lfus * ( ztmelts - rtt ) 
    119                zdiscrim     =  SQRT(  MAX( zbbb*zbbb - 4.0*zaaa*zccc, 0._wp )  ) 
    120                t_i_b(ji,jk) =  rtt - ( zbbb + zdiscrim ) / ( 2.0 *zaaa ) 
    121             END DO 
    122          END DO 
    123112         ! 
    124113      ENDIF  
     
    127116      !  Module 3 : Profile of salinity, constant in time                            | 
    128117      !------------------------------------------------------------------------------| 
    129  
    130118      IF(  num_sal == 3  )   CALL lim_var_salprof1d( kideb, kiut ) 
    131119 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limtrp.F90

    r4634 r4649  
    3030   USE limvar          ! clem for ice thickness correction 
    3131   USE timing          ! Timing 
     32   USE limcons        ! conservation tests 
    3233 
    3334   IMPLICIT NONE 
     
    7273      REAL(wp), POINTER, DIMENSION(:,:,:)    ::   zs0ice, zs0sn, zs0a, zs0c0 , zs0sm , zs0oi 
    7374      REAL(wp), POINTER, DIMENSION(:,:,:,:)  ::   zs0e 
    74       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 
    75       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax, zchk_umax ! Check errors (C Rousset) 
    7675      ! mass and salt flux (clem) 
    7776      REAL(wp), POINTER, DIMENSION(:,:,:) ::   zviold, zvsold   ! old ice volume... 
     
    7978      REAL(wp), POINTER, DIMENSION(:,:)   ::   zeiold, zesold   ! old enthalpies 
    8079      REAL(wp) :: zdv, zda, zvi, zvs, zsmv, zes, zei 
     80      ! 
     81      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    8182      !!--------------------------------------------------------------------- 
    8283      IF( nn_timing == 1 )  CALL timing_start('limtrp') 
     
    8788 
    8889      CALL wrk_alloc( jpi, jpj, jpl, zaiold, zhimax, zviold, zvsold )   ! clem 
    89  
    90       ! ------------------------------- 
    91       !- check conservation (C Rousset) 
    92       IF( ln_limdiahsb ) THEN 
    93          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    94          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    95          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    96          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    97          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    98          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    99       ENDIF 
    100       !- check conservation (C Rousset) 
    101       ! ------------------------------- 
    10290 
    10391      IF( numit == nstart .AND. lwp ) THEN 
     
    114102      IF( ln_limdyn ) THEN          !   Advection of sea ice properties   ! 
    115103         !                          !-------------------------------------! 
     104 
     105         ! conservation test 
     106         IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     107 
    116108         ! mass and salt flux init (clem) 
    117109         zviold(:,:,:)  = v_i(:,:,:) 
     
    368360 
    369361                 ! Update fluxes 
    370                   wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice  
    371                   wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
     362                  wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice  
     363                  wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    372364                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    373365                  hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     
    425417 
    426418                     ! Update mass fluxes 
    427                      wfx_res(ji,jj) = wfx_res(ji,jj) + ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
    428                      wfx_snw(ji,jj) = wfx_snw(ji,jj) + ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
     419                     wfx_res(ji,jj) = wfx_res(ji,jj) - ( v_i(ji,jj,jl) - zvi ) * rhoic * r1_rdtice 
     420                     wfx_snw(ji,jj) = wfx_snw(ji,jj) - ( v_s(ji,jj,jl) - zvs ) * rhosn * r1_rdtice 
    429421                     sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsmv ) * rhoic * r1_rdtice  
    430422                     hfx_res(ji,jj) = hfx_res(ji,jj) + ( e_s(ji,jj,1,jl) - zes ) * unit_fac / area(ji,jj) * r1_rdtice ! W.m-2 <0 
     
    474466            END DO 
    475467         END DO       
     468 
     469         ! conservation test 
     470         IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limtrp', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    476471 
    477472      ENDIF 
     
    508503         END DO 
    509504      ENDIF 
    510       ! ------------------------------- 
    511       !- check conservation (C Rousset) 
    512       IF( ln_limdiahsb ) THEN 
    513          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    514          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    515          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    516   
    517          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    518          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    519          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    520  
    521          zchk_vmin = glob_min(v_i) 
    522          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    523          zchk_amin = glob_min(a_i) 
    524          zchk_umax = glob_max(SQRT(u_ice**2 + v_ice**2)) 
    525  
    526          IF(lwp) THEN 
    527             IF ( ABS( zchk_v_i   ) >  1.e-4 ) THEN 
    528                WRITE(numout,*) 'violation volume [kg/day]     (limtrp) = ',(zchk_v_i * rday) 
    529                WRITE(numout,*) 'u_ice max [m/s]               (limtrp) = ',zchk_umax 
    530                WRITE(numout,*) 'number of time steps          (limtrp) =',kt 
    531             ENDIF 
    532             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limtrp) = ',(zchk_smv * rday) 
    533             IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limtrp) = ',(zchk_e_i) 
    534             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limtrp) = ',(zchk_vmin * 1.e-3) 
    535             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limtrp) = ',zchk_amin 
    536          ENDIF 
    537       ENDIF 
    538       !- check conservation (C Rousset) 
    539       ! ------------------------------- 
    540505      ! 
    541506      CALL wrk_dealloc( jpi, jpj, zui_u, zvi_v, zsm, zs0at, zs0ow, zeiold, zesold ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate1.F90

    r4634 r4649  
    3838   USE wrk_nemo         ! work arrays 
    3939   USE lib_fortran     ! glob_sum 
    40    ! Check budget (Rousset) 
    4140   USE in_out_manager   ! I/O manager 
    4241   USE iom              ! I/O manager 
    4342   USE lib_mpp          ! MPP library 
    4443   USE timing          ! Timing 
     44   USE limcons        ! conservation tests 
    4545 
    4646   IMPLICIT NONE 
     
    8282      REAL(wp) ::   zh, zdvres, zsal 
    8383      REAL(wp) ::   zat_i_old, ztmelts 
    84  
    85       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 
    86       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     84      ! 
     85      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    8786      !!------------------------------------------------------------------- 
    8887      IF( nn_timing == 1 )  CALL timing_start('limupdate1') 
    8988 
    90       !------------------------------------------------------------------------------ 
    91       ! 1. Update of Global variables                                               | 
    92       !------------------------------------------------------------------------------ 
    93  
    94       !----------------- 
    95       !  Trend terms 
    96       !----------------- 
    97  
    98       ! ------------------------------- 
    99       !- check conservation (C Rousset) 
    100       IF (ln_limdiahsb) THEN 
    101          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    102          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    103          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    104          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    105          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    106          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    107       ENDIF 
    108       !- check conservation (C Rousset) 
    109       ! ------------------------------- 
    110  
     89      IF( ln_limdyn ) THEN  
     90 
     91      ! conservation test 
     92      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     93 
     94      !----------------- 
    11195      ! zap small values 
    11296      !----------------- 
     
    11599      CALL lim_var_glo2eqv 
    116100      
    117       at_i(:,:) = 0._wp 
    118       DO jl = 1, jpl 
    119          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    120       END DO 
    121  
    122101      !---------------------------------------------------- 
    123       ! 2.2) Rebin categories with thickness out of bounds 
     102      ! Rebin categories with thickness out of bounds 
    124103      !---------------------------------------------------- 
    125104      DO jm = 1, jpm 
     
    134113      END DO 
    135114 
    136       !--- 2.13 ice concentration should not exceed amax  
     115      !---------------------------------------------------- 
     116      ! ice concentration should not exceed amax  
    137117      !----------------------------------------------------- 
    138118      DO jl  = 1, jpl 
     
    147127      END DO 
    148128 
    149       at_i(:,:) = 0.0 
     129      at_i(:,:) = 0._wp 
    150130      DO jl = 1, jpl 
    151131         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    152132      END DO 
    153    
    154    
     133     
     134      ! -------------------------------------- 
    155135      ! Final thickness distribution rebinning 
    156136      ! -------------------------------------- 
     
    163143      END DO 
    164144 
    165  
     145      !----------------- 
    166146      ! zap small values 
    167147      !----------------- 
     
    169149 
    170150      !--------------------- 
    171       ! 2.11) Ice salinity 
     151      ! Ice salinity bounds 
    172152      !--------------------- 
    173153      IF (  num_sal == 2  ) THEN  
     
    179159                  ! salinity stays in bounds 
    180160                  i_ice_switch    = 1._wp - MAX( 0._wp, SIGN( 1._wp, - v_i(ji,jj,jl) ) ) 
    181                   smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) !+ s_i_min * ( 1._wp - i_ice_switch ) * v_i(ji,jj,jl) 
     161                  smv_i(ji,jj,jl) = i_ice_switch * MAX( MIN( s_i_max * v_i(ji,jj,jl), smv_i(ji,jj,jl) ), s_i_min * v_i(ji,jj,jl) ) 
    182162                  ! associated salt flux 
    183163                  sfx_res(ji,jj) = sfx_res(ji,jj) - ( smv_i(ji,jj,jl) - zsal ) * rhoic * r1_rdtice 
    184                END DO ! ji 
    185             END DO ! jj 
    186          END DO !jl 
     164               END DO 
     165            END DO 
     166         END DO 
    187167      ENDIF 
    188  
    189       ! ------------------- 
    190       at_i(:,:) = a_i(:,:,1) 
    191       DO jl = 2, jpl 
    192          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    193       END DO 
    194    
    195168 
    196169      ! ------------------------------------------------- 
     
    206179      d_oa_i_trp (:,:,:)   = oa_i (:,:,:)   - old_oa_i (:,:,:) 
    207180      d_smv_i_trp(:,:,:)   = 0._wp 
    208       IF(  num_sal == 2  )   d_smv_i_trp(:,:,:)  = smv_i(:,:,:) - old_smv_i(:,:,:) 
    209  
    210       ! ------------------------------- 
    211       !- check conservation (C Rousset) 
    212       IF( ln_limdiahsb ) THEN 
    213          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    214          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    215          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    216   
    217          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    218          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    219          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    220  
    221          zchk_vmin = glob_min(v_i) 
    222          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    223          zchk_amin = glob_min(a_i) 
    224  
    225          IF(lwp) THEN 
    226             IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate1) = ',(zchk_v_i * rday) 
    227             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate1) = ',(zchk_smv * rday) 
    228             IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate1) = ',(zchk_e_i) 
    229             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate1) = ',(zchk_vmin * 1.e-3) 
    230             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate1) = ',zchk_amax 
    231             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate1) = ',zchk_amin 
    232          ENDIF 
    233        ENDIF 
    234       !- check conservation (C Rousset) 
    235       ! ------------------------------- 
     181      IF(  num_sal == 2  ) d_smv_i_trp(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:) 
     182 
     183      ! conservation test 
     184      IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limupdate1', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    236185 
    237186      IF(ln_ctl) THEN   ! Control print 
     
    302251      ENDIF 
    303252    
     253      ENDIF ! ln_limdyn 
     254 
    304255      IF( nn_timing == 1 )  CALL timing_stop('limupdate1') 
    305256   END SUBROUTINE lim_update1 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limupdate2.F90

    r4634 r4649  
    3939   USE lib_fortran     ! glob_sum 
    4040   USE timing          ! Timing 
     41   USE limcons        ! conservation tests 
    4142 
    4243   IMPLICIT NONE 
     
    8485      REAL(wp) ::   zdE          ! specific enthalpy difference (J/kg) 
    8586      REAL(wp) ::   zfmdt        ! exchange mass flux x time step (J/m2), >0 towards the ocean 
    86  
    87       REAL(wp) :: zchk_v_i, zchk_smv, zchk_e_i, zchk_fs, zchk_fw, zchk_ft, zchk_v_i_b, zchk_smv_b, zchk_e_i_b, zchk_fs_b, zchk_fw_b, zchk_ft_b ! Check conservation (C Rousset) 
    88       REAL(wp) :: zchk_vmin, zchk_amin, zchk_amax ! Check errors (C Rousset) 
     87      ! 
     88      REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b  
    8989      !!------------------------------------------------------------------- 
    9090      IF( nn_timing == 1 )  CALL timing_start('limupdate2') 
    9191 
    92       !---------------------------------------------------------------------------------------- 
    93       ! 1. Computation of trend terms       
    94       !---------------------------------------------------------------------------------------- 
    95  
    96       ! ------------------------------- 
    97       !- check conservation (C Rousset) 
    98       IF (ln_limdiahsb) THEN 
    99          zchk_v_i_b = glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) 
    100          zchk_smv_b = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) 
    101          zchk_e_i_b = glob_sum( SUM(   e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) 
    102          zchk_fw_b  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) 
    103          zchk_fs_b  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) 
    104          zchk_ft_b  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) 
    105       ENDIF 
    106       !- check conservation (C Rousset) 
    107       ! ------------------------------- 
    108  
     92      ! conservation test 
     93      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
     94 
     95      !----------------- 
    10996      ! zap small values 
    11097      !----------------- 
     
    113100      CALL lim_var_glo2eqv 
    114101 
    115       !-------------------------------------- 
    116       ! 2. Review of all pathological cases 
    117       !-------------------------------------- 
    118       at_i(:,:) = 0._wp 
    119       DO jl = 1, jpl 
    120          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    121       END DO 
    122  
    123102      !---------------------------------------------------- 
    124       ! 2.2) Rebin categories with thickness out of bounds 
     103      ! Rebin categories with thickness out of bounds 
    125104      !---------------------------------------------------- 
    126105      DO jm = 1, jpm 
     
    130109      END DO 
    131110 
    132  
    133 !clem debug: it is done in limthd_dh now 
    134 !      ! Melt of snow 
    135 !      !-------------- 
    136 !      DO jl = 1, jpl 
    137 !         DO jj = 1, jpj  
    138 !            DO ji = 1, jpi 
    139 !               IF( v_s(ji,jj,jl) >= epsi20 ) THEN 
    140 !                  ! If snow energy of melting smaller then Lf 
    141 !                  ! Then all snow melts and heat go to the ocean 
    142 !                  !IF ( zEs <= lfus ) THEN  
    143 !                  IF( t_s(ji,jj,1,jl) >= rtt ) THEN 
    144 !                     zdvres = - v_s(ji,jj,jl) 
    145 !                     zEs    = - e_s(ji,jj,1,jl) * unit_fac / ( area(ji,jj) * MAX( v_s(ji,jj,jl), epsi20 ) )  ! snow energy of melting (J.m-3) 
    146 !                     ! Contribution to heat flux to the ocean [W.m-2], < 0   
    147 !                     hfx_res(ji,jj) = hfx_res(ji,jj) - zEs * zdvres * r1_rdtice 
    148 !                     ! Contribution to mass flux 
    149 !                     wfx_snw(ji,jj) =  wfx_snw(ji,jj) + rhosn * zdvres * r1_rdtice 
    150 !                     ! updates 
    151 !                     v_s (ji,jj,jl)   = 0._wp 
    152 !                     ht_s(ji,jj,jl)   = 0._wp 
    153 !                     e_s (ji,jj,1,jl) = 0._wp 
    154 !                     t_s (ji,jj,1,jl) = rtt 
    155 !                  ENDIF 
    156 !               ENDIF 
    157 !            END DO 
    158 !         END DO 
    159 !      END DO 
    160 !clem debug 
    161  
    162       !--- 2.12 Constrain the thickness of the smallest category above 10 cm 
     111      !---------------------------------------------------------------------- 
     112      ! Constrain the thickness of the smallest category above hiclim 
    163113      !---------------------------------------------------------------------- 
    164114      DO jm = 1, jpm 
     
    176126      END DO !jm 
    177127       
    178       !--- 2.13 ice concentration should not exceed amax  
    179128      !----------------------------------------------------- 
    180       at_i(:,:) = 0.0 
     129      ! ice concentration should not exceed amax  
     130      !----------------------------------------------------- 
     131      at_i(:,:) = 0._wp 
    181132      DO jl = 1, jpl 
    182133         at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
     
    199150      END DO 
    200151 
     152      ! -------------------------------------- 
    201153      ! Final thickness distribution rebinning 
    202154      ! -------------------------------------- 
     
    209161      END DO 
    210162 
     163      !----------------- 
    211164      ! zap small values 
    212165      !----------------- 
     
    232185      ENDIF 
    233186 
    234  
    235       ! ------------------- 
    236       at_i(:,:) = a_i(:,:,1) 
    237       DO jl = 2, jpl 
    238          at_i(:,:) = a_i(:,:,jl) + at_i(:,:) 
    239       END DO 
    240  
    241187      !------------------------------------------------------------------------------ 
    242188      ! 2) Corrections to avoid wrong values                                        | 
     
    246192      DO jj = 2, jpjm1 
    247193         DO ji = 2, jpim1 
    248             IF ( at_i(ji,jj) .EQ. 0.0 ) THEN ! what to do if there is no ice 
    249                IF ( at_i(ji+1,jj) .EQ. 0.0 ) u_ice(ji,jj)   = 0.0 ! right side 
    250                IF ( at_i(ji-1,jj) .EQ. 0.0 ) u_ice(ji-1,jj) = 0.0 ! left side 
    251                IF ( at_i(ji,jj+1) .EQ. 0.0 ) v_ice(ji,jj)   = 0.0 ! upper side 
    252                IF ( at_i(ji,jj-1) .EQ. 0.0 ) v_ice(ji,jj-1) = 0.0 ! bottom side 
     194            IF ( at_i(ji,jj) == 0._wp ) THEN ! what to do if there is no ice 
     195               IF ( at_i(ji+1,jj) == 0._wp ) u_ice(ji,jj)   = 0._wp ! right side 
     196               IF ( at_i(ji-1,jj) == 0._wp ) u_ice(ji-1,jj) = 0._wp ! left side 
     197               IF ( at_i(ji,jj+1) == 0._wp ) v_ice(ji,jj)   = 0._wp ! upper side 
     198               IF ( at_i(ji,jj-1) == 0._wp ) v_ice(ji,jj-1) = 0._wp ! bottom side 
    253199            ENDIF 
    254200         END DO 
     
    261207      v_ice(:,:) = v_ice(:,:) * tmv(:,:) 
    262208  
    263  
    264209      ! ------------------------------------------------- 
    265210      ! Diagnostics 
     
    276221      dv_dt_thd(:,:,:) = d_v_i_thd(:,:,:) * r1_rdtice * rday 
    277222 
    278       ! ------------------------------- 
    279       !- check conservation (C Rousset) 
    280       IF( ln_limdiahsb ) THEN 
    281          zchk_fs  = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fs_b 
    282          zchk_fw  = glob_sum( ( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) ) * area(:,:) * tms(:,:) ) - zchk_fw_b 
    283          zchk_ft  = glob_sum( ( hfx_tot(:,:) - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) ) * area(:,:) / unit_fac * tms(:,:) ) - zchk_ft_b 
    284   
    285          zchk_v_i = ( glob_sum( SUM(   v_i(:,:,:)*rhoic + v_s(:,:,:)*rhosn, dim=3 ) * area(:,:) * tms(:,:) ) - zchk_v_i_b ) * r1_rdtice - zchk_fw  
    286          zchk_smv = ( glob_sum( SUM( smv_i(:,:,:), dim=3 ) * area(:,:) * tms(:,:) ) - zchk_smv_b ) * r1_rdtice + ( zchk_fs / rhoic ) 
    287          zchk_e_i =   glob_sum( SUM( e_i(:,:,1:nlay_i,:), dim=3 ) + SUM( e_s(:,:,1:nlay_s,:), dim=3 ) ) * r1_rdtice - zchk_e_i_b * r1_rdtice + zchk_ft 
    288  
    289          zchk_vmin = glob_min(v_i) 
    290          zchk_amax = glob_max(SUM(a_i,dim=3)) 
    291          zchk_amin = glob_min(a_i) 
    292  
    293          IF(lwp) THEN 
    294             IF ( ABS( zchk_v_i   ) >  1.e-4 ) WRITE(numout,*) 'violation volume [kg/day]     (limupdate2) = ',(zchk_v_i * rday) 
    295             IF ( ABS( zchk_smv   ) >  1.e-4 ) WRITE(numout,*) 'violation saline [psu*m3/day] (limupdate2) = ',(zchk_smv * rday) 
    296             IF ( ABS( zchk_e_i   ) >  1.e-2 ) WRITE(numout,*) 'violation enthalpy [1e9 J]    (limupdate2) = ',(zchk_e_i) 
    297             IF ( zchk_vmin <  0.            ) WRITE(numout,*) 'violation v_i<0  [mm]         (limupdate2) = ',(zchk_vmin * 1.e-3) 
    298             IF ( zchk_amax >  amax+epsi10   ) WRITE(numout,*) 'violation a_i>amax            (limupdate2) = ',zchk_amax 
    299             IF ( zchk_amin <  0.            ) WRITE(numout,*) 'violation a_i<0               (limupdate2) = ',zchk_amin 
    300          ENDIF 
    301        ENDIF 
    302       !- check conservation (C Rousset) 
    303       ! ------------------------------- 
     223      ! heat content variation (W.m-2) 
     224      DO jj = 1, jpj 
     225         DO ji = 1, jpi             
     226            diag_heat_dhc(ji,jj) = ( SUM( d_e_i_trp(ji,jj,1:nlay_i,:) + d_e_i_thd(ji,jj,1:nlay_i,:) ) +  &  
     227               &                     SUM( d_e_s_trp(ji,jj,1:nlay_s,:) + d_e_s_thd(ji,jj,1:nlay_s,:) ) ) * unit_fac * r1_rdtice / area(ji,jj)    
     228         END DO 
     229      END DO 
     230 
     231      ! conservation test 
     232      IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limupdate2', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 
    304233 
    305234      IF(ln_ctl) THEN   ! Control print 
     
    370299    
    371300      IF( nn_timing == 1 )  CALL timing_stop('limupdate2') 
     301 
    372302   END SUBROUTINE lim_update2 
    373303#else 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r4635 r4649  
    102102      CALL iom_put( "isst"         , sst_m               )       ! sea surface temperature 
    103103      CALL iom_put( "isss"         , sss_m               )       ! sea surface salinity 
    104       CALL iom_put( "qt_oce"       , qns + qsr           )       ! total flux at ocean surface 
    105104      ! 
    106105      DO jj = 2 , jpjm1 
     
    120119      CALL iom_put( "vice_ipa"     , z2db                )       ! ice velocity v component 
    121120      CALL iom_put( "icevel"       , z2d                 )       ! ice velocity module 
    122 !!SF BE CAREFUL : qsr_oce qnd qns_oce are after penetration over ice 
    123       CALL iom_put( "qsr_oce"      , qsr                 )       ! solar flux at ocean surface 
    124       CALL iom_put( "qns_oce"      , qns                 )       ! non-solar flux at ocean surface 
    125 !!SF end be careful 
    126121      CALL iom_put( "utau_ice"     , utau_ice            )       ! wind stress over ice along i-axis at I-point 
    127122      CALL iom_put( "vtau_ice"     , vtau_ice            )       ! wind stress over ice along j-axis at I-point 
    128 !!SF commented because this computation is not ok 
    129  !SF because qsr is not qsr_ocean but it contains already qsr_ice 
    130       !SF 
    131       !SF DO jj = 1 , jpj 
    132       !SF    DO ji = 1 , jpi 
    133       !SF       z2d(ji,jj)  = ( 1._wp - at_i(ji,jj) ) * qsr(ji,jj) 
    134       !SF   END DO 
    135       !SF END DO 
    136       !SF CALL iom_put( "qsr_io"       , z2d                 )        ! solar flux at ice/ocean surface 
    137       !SF DO jj = 1 , jpj 
    138       !SF    DO ji = 1 , jpi 
    139       !SF       z2d(ji,jj)  = ( 1._wp - at_i(ji,jj) ) * qns(ji,jj) 
    140       !SF   END DO 
    141       !SF END DO 
    142       !SF CALL iom_put( "qns_io"       , z2d                 )        ! non-solar flux at ice/ocean surface 
    143123      CALL iom_put( "snowpre"      , sprecip             )        ! snow precipitation  
    144124      CALL iom_put( "micesalt"     , smt_i               )        ! mean ice salinity 
     
    211191      CALL iom_put( "vfxsnw"     , wfx_snw * rday / rhoic  )        ! total snw growth/melt  
    212192      CALL iom_put( "vfxsub"     , wfx_sub * rday / rhoic  )        ! sublimation (snow)  
    213  
    214        CALL iom_put ('hfxdhc1', diag_heat_dhc1(:,:) )          ! Heat content variation in snow and ice  
    215        CALL iom_put ('hfxspr', hfx_spr(:,:) )          ! Heat content of snow precip  
    216        CALL iom_put ('hfxqsr', qsr(:,:) )          !  solar fluxes used by snw/ice 
    217        CALL iom_put ('hfxqns', qns(:,:) )          !  non solar fluxes used by snw/ice 
     193      CALL iom_put( "vfxspr"     , wfx_spr * rday / rhoic  )        ! precip (snow)  
    218194 
    219195       CALL iom_put ('hfxthd', hfx_thd(:,:) )   !   
     
    222198       CALL iom_put ('hfxout', hfx_out(:,:) )   !   
    223199       CALL iom_put ('hfxin' , hfx_in(:,:) )   !   
    224        CALL iom_put ('hfxtot', hfx_tot(:,:) )   !   
    225200       CALL iom_put ('hfxsnw', hfx_snw(:,:) )   !   
    226201       CALL iom_put ('hfxsub', hfx_sub(:,:) )   !   
    227202       CALL iom_put ('hfxerr', hfx_err(:,:) )   !   
    228203       CALL iom_put ('hfxerr_rem', hfx_err_rem(:,:) )   !   
     204 
     205       CALL iom_put ('hfxsum', hfx_sum(:,:) )   !   
     206       CALL iom_put ('hfxbom', hfx_bom(:,:) )   !   
     207       CALL iom_put ('hfxbog', hfx_bog(:,:) )   !   
     208       CALL iom_put ('hfxdif', hfx_dif(:,:) )   !   
     209       CALL iom_put ('hfxopw', hfx_opw(:,:) )   !   
     210       CALL iom_put ('hfxtur', fhtur(:,:) * at_i(:,:) )   ! turbulent heat flux at ice base  
     211       CALL iom_put ('hfxdhc', diag_heat_dhc(:,:) )          ! Heat content variation in snow and ice  
     212       CALL iom_put ('hfxspr', hfx_spr(:,:) )          ! Heat content of snow precip  
    229213 
    230214      !-------------------------------- 
     
    327311      CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    328312      CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    329       CALL histdef( kid, "iisfxthd", "Salt flux from thermo"   , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    330313      CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    331314      CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""      , jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     
    354337      CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) ) 
    355338      CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) ) 
    356       !CALL histwrite( kid, "iisfxthd", kt, sfx_thd        , jpi*jpj, (/1/) ) 
    357339      CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) ) 
    358340      CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/thd_ice.F90

    r4634 r4649  
    5757   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   t_bo_b        !: <==> the 2D  t_bo 
    5858 
    59    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_tot_1d 
    60    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_thd_1d 
    61    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_spr_1d 
     59   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sum_1d 
     60   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_bom_1d 
     61   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_bog_1d 
     62   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_dif_1d 
     63   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_opw_1d 
    6264   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_snw_1d 
    63    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sub_1d 
    64    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_res_1d 
    6565   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_1d 
    6666   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_err_rem_1d 
     67 
     68   ! heat flux associated with ice-atmosphere mass exchange 
     69   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_sub_1d 
     70   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_spr_1d 
     71 
     72   ! heat flux associated with ice-ocean mass exchange 
     73   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_thd_1d 
     74   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   hfx_res_1d 
    6775 
    6876   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_ice_1d    !: <==> the 2D  wfx_ice 
     
    7583   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_sni_1d    !: <==> the 2D  wfx_ice 
    7684   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_opw_1d    !: <==> the 2D  wfx_ice 
     85   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_res_1d    !: <==> the 2D  wfx_ice 
     86   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   wfx_spr_1d    !: <==> the 2D  wfx_ice 
    7787 
    7888   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_bri_1d    !: <==> the 2D sfx_bri 
     
    8292   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_sni_1d    !:  
    8393   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_opw_1d    !: 
     94   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sfx_res_1d    !: 
    8495 
    8596   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   sprecip_1d    !: <==> the 2D  sprecip 
    8697   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   frld_1d       !: <==> the 2D  frld 
    87    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  frld 
     98   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   at_i_b        !: <==> the 2D  at_i 
    8899   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhtur_1d       !: <==> the 2D  fhtur 
    89100   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) ::   fhld_1d       !: <==> the 2D  fhld 
     
    150161         &      qsr_ice_1d (jpij) ,     & 
    151162         &      fr1_i0_1d(jpij) , fr2_i0_1d(jpij) , qns_ice_1d(jpij) ,     & 
    152          &      t_bo_b   (jpij) , iatte_1d   (jpij) ,     & 
    153          &      oatte_1d (jpij) , hfx_tot_1d(jpij), hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , & 
    154          &      hfx_snw_1d(jpij), hfx_sub_1d(jpij), hfx_err_1d(jpij) , hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) ) 
     163         &      t_bo_b   (jpij) , iatte_1d  (jpij) , oatte_1d (jpij) ,     & 
     164         &      hfx_sum_1d(jpij) , hfx_bom_1d(jpij) ,hfx_bog_1d(jpij) ,hfx_dif_1d(jpij) ,hfx_opw_1d(jpij) , & 
     165         &      hfx_thd_1d(jpij) , hfx_spr_1d(jpij) , & 
     166         &      hfx_snw_1d(jpij) , hfx_sub_1d(jpij) , hfx_err_1d(jpij) , hfx_res_1d(jpij) , hfx_err_rem_1d(jpij),       STAT=ierr(1) ) 
    155167      ! 
    156168      ALLOCATE( sprecip_1d (jpij) , frld_1d    (jpij) , at_i_b     (jpij) ,     & 
    157          &      fhtur_1d   (jpij) , wfx_ice_1d (jpij) , wfx_snw_1d (jpij) ,     & 
    158          &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  & 
     169         &      fhtur_1d   (jpij) , wfx_ice_1d (jpij) , wfx_snw_1d (jpij) , wfx_spr_1d (jpij) ,     & 
     170         &      fhld_1d    (jpij) , wfx_sub_1d (jpij) , wfx_bog_1d(jpij) , wfx_bom_1d(jpij) , wfx_sum_1d(jpij) , wfx_sni_1d (jpij) , wfx_opw_1d (jpij) ,  wfx_res_1d (jpij) ,  & 
    159171         &      dqns_ice_1d(jpij) , qla_ice_1d (jpij) , dqla_ice_1d(jpij) ,     & 
    160172         &      tatm_ice_1d(jpij) ,      &    
    161173         &      i0         (jpij) ,     &   
    162          &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , & 
     174         &      sfx_bri_1d (jpij) , sfx_bog_1d (jpij) , sfx_bom_1d (jpij) ,sfx_sum_1d (jpij) ,sfx_sni_1d (jpij) , sfx_opw_1d (jpij) , sfx_res_1d (jpij) , & 
    163175         &      dsm_i_fl_1d(jpij) , dsm_i_gd_1d(jpij) , dsm_i_se_1d(jpij) ,     &      
    164176         &      dsm_i_si_1d(jpij) , hicol_b    (jpij)                     , STAT=ierr(2) ) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r4634 r4649  
    6969 
    7070   PUBLIC sbc_ice_lim  ! routine called by sbcmod.F90 
     71   PUBLIC lim_prt_state 
    7172    
    7273   !! * Substitutions 
     
    170171         ! 
    171172         IF( ln_nicep ) THEN      ! control print at a given point 
    172             jiindx = 3    ;   jjindx =  49 
     173            jiindx = 15    ;   jjindx =  44 
    173174            IF(lwp) WRITE(numout,*) ' The debugging point is : jiindx : ',jiindx, ' jjindx : ',jjindx 
    174175         ENDIF 
     
    183184         u_oce(:,:) = ssu_m(:,:)                     ! mean surface ocean current at ice velocity point 
    184185         v_oce(:,:) = ssv_m(:,:)                     ! (C-grid dynamics :  U- & V-points as the ocean) 
    185          ! 
    186          t_bo(:,:) = tfreez( sss_m ) +  rt0          ! masked sea surface freezing temperature [Kelvin] 
    187          !                                           ! (set to rt0 over land) 
     186 
     187         ! masked sea surface freezing temperature [Kelvin] 
     188         t_bo(:,:) = ( tfreez( sss_m ) +  rt0 ) * tmask(:,:,1) + rt0 * ( 1. - tmask(:,:,1) ) 
     189 
    188190         CALL albedo_ice( t_su, ht_i, ht_s, zalb_ice_cs, zalb_ice_os )  ! ... ice albedo 
    189191 
     
    310312         ! salt, heat and mass fluxes 
    311313         sfx    (:,:) = 0._wp   ; 
    312          sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp   ;   sfx_res(:,:) = 0._wp 
     314         sfx_bri(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp  
    313315         sfx_sni(:,:) = 0._wp   ;   sfx_opw(:,:) = 0._wp 
    314316         sfx_bog(:,:) = 0._wp   ;   sfx_dyn(:,:) = 0._wp 
    315317         sfx_bom(:,:) = 0._wp   ;   sfx_sum(:,:) = 0._wp 
    316  
    317          hfx_thd(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp   ;   hfx_snw(:,:) = 0._wp 
    318          hfx_tot(:,:) = 0._wp   ;   hfx_spr(:,:) = 0._wp   ;   hfx_res(:,:) = 0._wp 
    319          hfx_sub(:,:) = 0._wp   ;   hfx_err(:,:) = 0._wp   ;   hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
    320          hfx_err_rem(:,:) = 0._wp 
    321  
    322          wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     318         sfx_res(:,:) = 0._wp 
     319 
     320         wfx_snw(:,:) = 0._wp   ;   wfx_ice(:,:) = 0._wp 
    323321         wfx_sni(:,:) = 0._wp   ;   wfx_opw(:,:) = 0._wp 
    324322         wfx_bog(:,:) = 0._wp   ;   wfx_dyn(:,:) = 0._wp 
    325323         wfx_bom(:,:) = 0._wp   ;   wfx_sum(:,:) = 0._wp 
    326          wfx_res(:,:) = 0._wp   ;    
    327          ! 
    328          fhld  (:,:) = 0._wp  
    329          fmmflx(:,:) = 0._wp      
    330          ftr_ice(:,:,:) = 0._wp   ! part of solar radiation transmitted through the ice 
     324         wfx_res(:,:) = 0._wp   ;   wfx_sub(:,:) = 0._wp 
     325         wfx_spr(:,:) = 0._wp   ;    
     326 
     327         hfx_in (:,:) = 0._wp   ;   hfx_out(:,:) = 0._wp 
     328         hfx_thd(:,:) = 0._wp   ;    
     329         hfx_snw(:,:) = 0._wp   ;   hfx_opw(:,:) = 0._wp 
     330         hfx_bog(:,:) = 0._wp   ;   hfx_dyn(:,:) = 0._wp 
     331         hfx_bom(:,:) = 0._wp   ;   hfx_sum(:,:) = 0._wp 
     332         hfx_res(:,:) = 0._wp   ;   hfx_sub(:,:) = 0._wp 
     333         hfx_spr(:,:) = 0._wp   ;   hfx_dif(:,:) = 0._wp  
     334         hfx_err(:,:) = 0._wp   ;   hfx_err_rem(:,:) = 0._wp 
     335 
     336         ! 
     337         fhld  (:,:)    = 0._wp  
     338         fmmflx(:,:)    = 0._wp      
     339         ! part of solar radiation transmitted through the ice 
     340         ftr_ice(:,:,:) = 0._wp 
    331341 
    332342         ! diags 
    333          diag_trp_vi(:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp  ;  
    334          diag_heat_dhc1(:,:) = 0._wp   ;    
     343         diag_trp_vi  (:,:) = 0._wp  ; diag_trp_vs(:,:) = 0._wp  ;  diag_trp_ei(:,:) = 0._wp  ;  diag_trp_es(:,:) = 0._wp 
     344         diag_heat_dhc(:,:) = 0._wp   
    335345 
    336346         ! dynamical invariants 
     
    820830               WRITE(numout,*) ' hfx_in       : ', hfx_in(ji,jj) 
    821831               WRITE(numout,*) ' hfx_out      : ', hfx_out(ji,jj) 
    822                WRITE(numout,*) ' hfx_tot      : ', hfx_tot(ji,jj) 
    823                WRITE(numout,*) ' dhc          : ', diag_heat_dhc1(ji,jj)               
     832               WRITE(numout,*) ' dhc          : ', diag_heat_dhc(ji,jj)               
    824833               WRITE(numout,*) 
    825834               WRITE(numout,*) ' hfx_dyn      : ', hfx_dyn(ji,jj) 
  • branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/OPA_SRC/SBC/sbcmod.F90

    r4634 r4649  
    391391         CALL iom_put( "qt"    , qns  + qsr )                   ! total heat flux  
    392392         CALL iom_put( "qns"   , qns        )                   ! solar heat flux 
    393          CALL iom_put( "qsr"   ,       qsr )                   ! solar heat flux 
     393         CALL iom_put( "qsr"   ,        qsr )                   ! solar heat flux 
    394394         IF( nn_ice > 0 )   CALL iom_put( "ice_cover", fr_i )   ! ice fraction  
    395395      ENDIF 
Note: See TracChangeset for help on using the changeset viewer.