Changeset 7067


Ignore:
Timestamp:
2016-10-21T17:35:49+02:00 (4 years ago)
Author:
cetlod
Message:

Offline vvl phased with the head of 3.6 stable

Location:
branches/2015/dev_r6946_Offline_vvl/NEMOGCM
Files:
19 edited

Legend:

Unmodified
Added
Removed
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml

    r6316 r7067  
    8888     <field field_ref="u_masstr"     name="vozomatr" /> 
    8989     <field field_ref="u_heattr"     name="sozohetr" /> 
    90       <field field_ref="u_salttr"     name="sozosatr" /> 
     90     <field field_ref="u_salttr"     name="sozosatr" /> 
    9191   </file> 
    9292    
     
    9999     <field field_ref="v_masstr"     name="vomematr" /> 
    100100     <field field_ref="v_heattr"     name="somehetr" /> 
    101       <field field_ref="v_salttr"     name="somesatr" /> 
     101     <field field_ref="v_salttr"     name="somesatr" /> 
    102102   </file> 
    103103    
     
    116116          <field field_ref="iceconc"         name="siconc" /> 
    117117 
    118           <field field_ref="vfxbog"          name="vfxbog" /> 
    119           <field field_ref="vfxdyn"          name="vfxdyn" /> 
    120           <field field_ref="vfxopw"          name="vfxopw" /> 
    121           <field field_ref="vfxsni"          name="vfxsni" /> 
    122           <field field_ref="vfxsum"          name="vfxsum" /> 
    123           <field field_ref="vfxbom"          name="vfxbom" /> 
    124           <field field_ref="vfxres"          name="vfxres" /> 
    125118          <field field_ref="vfxice"          name="vfxice" /> 
    126119          <field field_ref="vfxsnw"          name="vfxsnw" /> 
    127120          <field field_ref="vfxsub"          name="vfxsub" /> 
     121          <field field_ref="vfxsub_err"      name="vfxsub_err" /> 
    128122          <field field_ref="vfxspr"          name="vfxspr" /> 
    129123 
     
    134128          <field field_ref="destrp"          name="destrp" /> 
    135129 
    136           <field field_ref="sfxbri"          name="sfxbri" /> 
    137           <field field_ref="sfxdyn"          name="sfxdyn" /> 
    138           <field field_ref="sfxres"          name="sfxres" /> 
    139           <field field_ref="sfxbog"          name="sfxbog" /> 
    140           <field field_ref="sfxbom"          name="sfxbom" /> 
    141           <field field_ref="sfxsum"          name="sfxsum" /> 
    142           <field field_ref="sfxsni"          name="sfxsni" /> 
    143           <field field_ref="sfxopw"          name="sfxopw" /> 
    144130          <field field_ref="sfx"             name="sfx"    /> 
    145131 
    146           <field field_ref="hfxsum"          name="hfxsum"     /> 
    147           <field field_ref="hfxbom"          name="hfxbom"     /> 
    148           <field field_ref="hfxbog"          name="hfxbog"     /> 
    149           <field field_ref="hfxdif"          name="hfxdif"     /> 
    150           <field field_ref="hfxopw"          name="hfxopw"     /> 
    151132          <field field_ref="hfxout"          name="hfxout"     /> 
    152133          <field field_ref="hfxin"           name="hfxin"      /> 
    153           <field field_ref="hfxsnw"          name="hfxsnw"     /> 
    154           <field field_ref="hfxerr"          name="hfxerr"     /> 
    155           <field field_ref="hfxerr_rem"      name="hfxerr_rem" /> 
    156  
    157      <!-- ice-ocean heat flux from mass exchange --> 
    158           <field field_ref="hfxdyn"          name="hfxdyn" /> 
    159           <field field_ref="hfxres"          name="hfxres" /> 
    160           <field field_ref="hfxthd"          name="hfxthd" /> 
    161      <!-- ice-atm. heat flux from mass exchange --> 
    162           <field field_ref="hfxsub"          name="hfxsub" /> 
    163           <field field_ref="hfxspr"          name="hfxspr" /> 
     134 
    164135 
    165136     <!-- diags --> 
    166           <field field_ref="hfxdhc"          name="hfxdhc" /> 
    167           <field field_ref="hfxtur"          name="hfxtur" /> 
    168  
    169137          <field field_ref="isst"            name="sst"    /> 
    170138          <field field_ref="isss"            name="sss"    /> 
     
    209177          <field field_ref="bgsaline"     name="bgsaline"   /> 
    210178          <field field_ref="bgheatco"     name="bgheatco"   /> 
     179          <field field_ref="bgheatfx"     name="bgheatfx"   /> 
    211180          <field field_ref="bgsaltco"     name="bgsaltco"   /> 
    212181          <field field_ref="bgvolssh"     name="bgvolssh"   />  
     
    214183          <field field_ref="bgfrcvol"     name="bgfrcvol"   /> 
    215184          <field field_ref="bgfrctem"     name="bgfrctem"   /> 
     185          <field field_ref="bgfrchfx"     name="bgfrchfx"   /> 
    216186          <field field_ref="bgfrcsal"     name="bgfrcsal"   /> 
    217187 
    218           <field field_ref="ibgvoltot"    name="ibgvoltot"  /> 
    219           <field field_ref="sbgvoltot"    name="sbgvoltot"  /> 
    220           <field field_ref="ibgarea"      name="ibgarea"    /> 
    221           <field field_ref="ibgsaline"    name="ibgsaline"  /> 
    222           <field field_ref="ibgtemper"    name="ibgtemper"  /> 
    223           <field field_ref="ibgheatco"    name="ibgheatco"  /> 
    224           <field field_ref="sbgheatco"    name="sbgheatco"  /> 
    225           <field field_ref="ibgsaltco"    name="ibgsaltco"  /> 
    226  
    227           <field field_ref="ibgvfx"       name="ibgvfx"     /> 
    228           <field field_ref="ibgvfxbog"    name="ibgvfxbog"  /> 
    229           <field field_ref="ibgvfxopw"    name="ibgvfxopw"  /> 
    230           <field field_ref="ibgvfxsni"    name="ibgvfxsni"  /> 
    231           <field field_ref="ibgvfxdyn"    name="ibgvfxdyn"  /> 
    232           <field field_ref="ibgvfxbom"    name="ibgvfxbom"  /> 
    233           <field field_ref="ibgvfxsum"    name="ibgvfxsum"  /> 
    234           <field field_ref="ibgvfxres"    name="ibgvfxres"  /> 
    235           <field field_ref="ibgvfxspr"    name="ibgvfxspr"  /> 
    236           <field field_ref="ibgvfxsnw"    name="ibgvfxsnw"  /> 
    237           <field field_ref="ibgvfxsub"    name="ibgvfxsub"  /> 
    238  
    239           <field field_ref="ibgsfx"       name="ibgsfx"     /> 
    240           <field field_ref="ibgsfxbri"    name="ibgsfxbri"  /> 
    241           <field field_ref="ibgsfxdyn"    name="ibgsfxdyn"  /> 
    242           <field field_ref="ibgsfxres"    name="ibgsfxres"  /> 
    243           <field field_ref="ibgsfxbog"    name="ibgsfxbog"  /> 
    244           <field field_ref="ibgsfxopw"    name="ibgsfxopw"  /> 
    245           <field field_ref="ibgsfxsni"    name="ibgsfxsni"  /> 
    246           <field field_ref="ibgsfxbom"    name="ibgsfxbom"  /> 
    247           <field field_ref="ibgsfxsum"    name="ibgsfxsum"  /> 
    248  
    249           <field field_ref="ibghfxdhc"    name="ibghfxdhc"  /> 
    250           <field field_ref="ibghfxspr"    name="ibghfxspr"  /> 
    251  
    252           <field field_ref="ibghfxres"    name="ibghfxres"  /> 
    253           <field field_ref="ibghfxsub"    name="ibghfxsub"  /> 
    254           <field field_ref="ibghfxdyn"    name="ibghfxdyn"  /> 
    255           <field field_ref="ibghfxthd"    name="ibghfxthd"  /> 
    256           <field field_ref="ibghfxsum"    name="ibghfxsum"  /> 
    257           <field field_ref="ibghfxbom"    name="ibghfxbom"  /> 
    258           <field field_ref="ibghfxbog"    name="ibghfxbog"  /> 
    259           <field field_ref="ibghfxdif"    name="ibghfxdif"  /> 
    260           <field field_ref="ibghfxopw"    name="ibghfxopw"  /> 
    261           <field field_ref="ibghfxout"    name="ibghfxout"  /> 
    262           <field field_ref="ibghfxin"     name="ibghfxin"   /> 
    263           <field field_ref="ibghfxsnw"    name="ibghfxsnw"  /> 
    264  
    265           <field field_ref="ibgfrcvol"    name="ibgfrcvol"  /> 
    266           <field field_ref="ibgfrcsfx"    name="ibgfrcsfx"  /> 
    267           <field field_ref="ibgvolgrm"    name="ibgvolgrm"  /> 
     188     <field field_ref="ibgvol_tot"   name="ibgvol_tot"  /> 
     189     <field field_ref="sbgvol_tot"   name="sbgvol_tot"  /> 
     190     <field field_ref="ibgarea_tot"  name="ibgarea_tot" /> 
     191     <field field_ref="ibgsalt_tot"  name="ibgsalt_tot" /> 
     192     <field field_ref="ibgheat_tot"  name="ibgheat_tot" /> 
     193     <field field_ref="sbgheat_tot"  name="sbgheat_tot" /> 
     194      
     195     <field field_ref="ibgvolume"    name="ibgvolume"   /> 
     196     <field field_ref="ibgsaltco"    name="ibgsaltco"   /> 
     197     <field field_ref="ibgheatco"    name="ibgheatco"   /> 
     198          <field field_ref="ibgheatfx"    name="ibgheatfx"   /> 
     199      
     200     <field field_ref="ibgfrcvoltop" name="ibgfrcvoltop" /> 
     201     <field field_ref="ibgfrcvolbot" name="ibgfrcvolbot" /> 
     202     <field field_ref="ibgfrctemtop" name="ibgfrctemtop" /> 
     203     <field field_ref="ibgfrctembot" name="ibgfrctembot" /> 
     204     <field field_ref="ibgfrcsal"    name="ibgfrcsal"    /> 
     205          <field field_ref="ibgfrchfxtop" name="ibgfrchfxtop" /> 
     206          <field field_ref="ibgfrchfxbot" name="ibgfrchfxbot" /> 
    268207 
    269208        </file> 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/CONFIG/SHARED/field_def.xml

    r6471 r7067  
    225225         <field id="hflx_rain_cea" long_name="heat flux due to rainfall"                                standard_name="temperature_flux_due_to_rainfall_expressed_as_heat_flux_into_sea_water"        unit="W/m2"     /> 
    226226         <field id="hflx_evap_cea" long_name="heat flux due to evaporation"                             standard_name="temperature_flux_due_to_evaporation_expressed_as_heat_flux_out_of_sea_water"   unit="W/m2"     /> 
    227          <field id="hflx_snow_cea" long_name="heat flux due to snow falling over ice-free ocean"        standard_name="heat_flux_into_sea_water_due_to_snow_thermodynamics"                           unit="W/m2"     /> 
     227         <field id="hflx_snow_cea" long_name="heat flux due to snow falling"                            standard_name="heat_flux_onto_ocean_and_ice_due_to_snow_thermodynamics"                       unit="W/m2"     /> 
     228         <field id="hflx_snow_ai_cea" long_name="heat flux due to snow falling over ice"                standard_name="heat_flux_onto_ice_due_to_snow_thermodynamics"                                 unit="W/m2"     /> 
     229         <field id="hflx_snow_ao_cea" long_name="heat flux due to snow falling over ice-free ocean"     standard_name="heat_flux_onto_sea_water_due_to_snow_thermodynamics"                           unit="W/m2"     /> 
    228230         <field id="hflx_ice_cea"  long_name="heat flux due to ice thermodynamics"                      standard_name="heat_flux_into_sea_water_due_to_sea_ice_thermodynamics"                        unit="W/m2"     /> 
    229231         <field id="hflx_rnf_cea"  long_name="heat flux due to runoffs"                                 standard_name="temperature_flux_due_to_runoff_expressed_as_heat_flux_into_sea_water"          unit="W/m2"     /> 
     
    331333         <field id="vfxsnw"       long_name="snw melt/growth"                                              unit="m/day"   /> 
    332334         <field id="vfxsub"       long_name="snw sublimation"                                              unit="m/day"   /> 
     335         <field id="vfxsub_err"   long_name="excess of snw sublimation sent to ocean"                      unit="m/day"   /> 
    333336         <field id="vfxspr"       long_name="snw precipitation on ice"                                     unit="m/day"   /> 
    334337         <field id="vfxthin"      long_name="daily thermo ice prod. for thin ice(<20cm) + open water"      unit="m/day"   /> 
     
    499502       <field id="bgtemper"     long_name="drift in global mean temperature wrt timestep 1"                 standard_name="change_over_time_in_sea_water_potential_temperature"   unit="degC"     /> 
    500503       <field id="bgsaline"     long_name="drift in global mean salinity wrt timestep 1"                    standard_name="change_over_time_in_sea_water_practical_salinity"      unit="1e-3"     /> 
    501        <field id="bgheatco"     long_name="drift in global mean heat content wrt timestep 1"                                                                                      unit="10^9J"    /> 
    502        <field id="bgsaltco"     long_name="drift in global mean salt content wrt timestep 1"                                                                                      unit="1e-3*m3"  /> 
     504       <field id="bgheatco"     long_name="drift in global mean heat content wrt timestep 1"                                                                                      unit="1.e20J"   /> 
     505       <field id="bgheatfx"     long_name="drift in global mean heat flux    wrt timestep 1"                                                                                      unit="W/m2"     /> 
     506       <field id="bgsaltco"     long_name="drift in global mean salt content wrt timestep 1"                                                                                      unit="1e-3*km3" /> 
    503507       <field id="bgvolssh"     long_name="drift in global mean ssh volume wrt timestep 1"                                                                                        unit="km3"      /> 
    504508         <field id="bgvole3t"     long_name="drift in global mean volume variation (e3t) wrt timestep 1"                                                                            unit="km3"      /> 
    505        <field id="bgvoltot"     long_name="drift in global mean volume wrt timestep 1"                                                                                            unit="km3"      /> 
    506          <!-- NOTE: No matching iom_put call --> 
    507        <field id="bgsshtot"     long_name="drift in global mean ssh wrt timestep 1"                         standard_name="global_average_sea_level_change"                       unit="m"        /> 
    508        <field id="bgfrcvol"     long_name="drift in global mean volume from forcing wrt timestep 1"                                                                               unit="km3"      /> 
    509        <field id="bgfrctem"     long_name="drift in global mean heat content from forcing wrt timestep 1"                                                                         unit="10^9J"    /> 
    510        <field id="bgfrcsal"     long_name="drift in global mean salt content from forcing wrt timestep 1"                                                                         unit="1e-3*km3" /> 
    511        <field id="bgmistem"     long_name="global mean temperature error due to free surface"                                                                                     unit="degC"     /> 
    512        <field id="bgmissal"     long_name="global mean salinity error due to free surface"                                                                                        unit="1e-3"     /> 
     509       <field id="bgfrcvol"     long_name="global mean volume from forcing"                                                                                                       unit="km3"      /> 
     510       <field id="bgfrctem"     long_name="global mean heat content from forcing"                                                                                                 unit="1.e20J"   /> 
     511         <field id="bgfrchfx"     long_name="global mean heat flux from forcing"                                                                                                    unit="W/m2"     /> 
     512       <field id="bgfrcsal"     long_name="global mean salt content from forcing"                                                                                                 unit="1e-3*km3" /> 
     513       <field id="bgmistem"     long_name="global mean temperature error due to free surface (no vvl)"                                                                            unit="degC"     /> 
     514       <field id="bgmissal"     long_name="global mean salinity error due to free surface (no vvl)"                                                                               unit="1e-3"     /> 
    513515      </field_group> 
    514516 
     
    517519      <field_group id="SBC_scalar"  domain_ref="1point" > 
    518520         <!-- available with ln_limdiaout --> 
    519          <field id="ibgvoltot"    long_name="global mean ice volume"                                 unit="km3"        /> 
    520          <field id="sbgvoltot"    long_name="global mean snow volume"                                unit="km3"        /> 
    521          <field id="ibgarea"      long_name="global mean ice area"                                   unit="km2"        /> 
    522          <field id="ibgsaline"    long_name="global mean ice salinity"                               unit="1e-3"       /> 
    523          <field id="ibgtemper"    long_name="global mean ice temperature"                            unit="degC"       /> 
    524          <field id="ibgheatco"    long_name="global mean ice heat content"                           unit="10^20J"     /> 
    525          <field id="sbgheatco"    long_name="global mean snow heat content"                          unit="10^20J"     /> 
    526          <field id="ibgsaltco"    long_name="global mean ice salt content"                           unit="1e-3*km3"   /> 
    527  
    528          <field id="ibgvfx"       long_name="global mean volume flux (emp)"                          unit="m/day"      /> 
    529          <field id="ibgvfxbog"    long_name="global mean volume flux (bottom growth)"                unit="m/day"      /> 
    530          <field id="ibgvfxopw"    long_name="global mean volume flux (open water growth)"            unit="m/day"      /> 
    531          <field id="ibgvfxsni"    long_name="global mean volume flux (snow-ice growth)"              unit="m/day"      /> 
    532          <field id="ibgvfxdyn"    long_name="global mean volume flux (dynamic growth)"               unit="m/day"      /> 
    533          <field id="ibgvfxbom"    long_name="global mean volume flux (bottom melt)"                  unit="m/day"      /> 
    534          <field id="ibgvfxsum"    long_name="global mean volume flux (surface melt)"                 unit="m/day"      /> 
    535          <field id="ibgvfxres"    long_name="global mean volume flux (resultant)"                    unit="m/day"      /> 
    536          <field id="ibgvfxspr"    long_name="global mean volume flux (snow precip)"                  unit="m/day"      /> 
    537          <field id="ibgvfxsnw"    long_name="global mean volume flux (snow melt)"                    unit="m/day"      /> 
    538          <field id="ibgvfxsub"    long_name="global mean volume flux (snow sublimation)"             unit="m/day"      /> 
    539  
    540          <field id="ibgsfx"       long_name="global mean salt flux (total)"                          unit="1e-3*m/day" /> 
    541          <field id="ibgsfxbri"    long_name="global mean salt flux (brines)"                         unit="1e-3*m/day" /> 
    542          <field id="ibgsfxdyn"    long_name="global mean salt flux (dynamic)"                        unit="1e-3*m/day" /> 
    543          <field id="ibgsfxres"    long_name="global mean salt flux (resultant)"                      unit="1e-3*m/day" /> 
    544          <field id="ibgsfxbog"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" /> 
    545          <field id="ibgsfxopw"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" /> 
    546          <field id="ibgsfxsni"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" /> 
    547          <field id="ibgsfxbom"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" /> 
    548          <field id="ibgsfxsum"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" /> 
    549          <field id="ibgsfxsub"    long_name="global mean salt flux (thermo)"                         unit="1e-3*m/day" /> 
    550  
    551          <field id="ibghfxdhc"    long_name="Heat content variation in snow and ice"                 unit="W"          /> 
    552          <field id="ibghfxspr"    long_name="Heat content of snow precip"                            unit="W"          /> 
    553  
    554          <field id="ibghfxthd"    long_name="heat fluxes from ice-ocean exchange during thermo"      unit="W"          /> 
    555          <field id="ibghfxsum"    long_name="heat fluxes causing surface ice melt"                   unit="W"          /> 
    556          <field id="ibghfxbom"    long_name="heat fluxes causing bottom ice melt"                    unit="W"          /> 
    557          <field id="ibghfxbog"    long_name="heat fluxes causing bottom ice growth"                  unit="W"          /> 
    558          <field id="ibghfxdif"    long_name="heat fluxes causing ice temperature change"             unit="W"          /> 
    559          <field id="ibghfxopw"    long_name="heat fluxes causing open water ice formation"           unit="W"          /> 
    560          <field id="ibghfxdyn"    long_name="heat fluxes from ice-ocean exchange during dynamic"     unit="W"          /> 
    561          <field id="ibghfxres"    long_name="heat fluxes from ice-ocean exchange during resultant"   unit="W"          /> 
    562          <field id="ibghfxsub"    long_name="heat fluxes from sublimation"                           unit="W"          /> 
    563          <field id="ibghfxsnw"    long_name="heat fluxes from snow-ocean exchange"                   unit="W"          /> 
    564          <field id="ibghfxout"    long_name="non solar heat fluxes received by the ocean"            unit="W"          /> 
    565          <field id="ibghfxin"     long_name="total heat fluxes at the ice surface"                   unit="W"          /> 
    566  
    567          <field id="ibgfrcvol"    long_name="global mean forcing volume (emp)"                       unit="km3"        /> 
    568          <field id="ibgfrcsfx"    long_name="global mean forcing salt   (sfx)"                       unit="1e-3*km3"   /> 
    569          <field id="ibgvolgrm"    long_name="global mean ice growth+melt volume"                     unit="km3"        /> 
     521         <field id="ibgfrcvoltop"    long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)"   unit="km3"       /> 
     522         <field id="ibgfrcvolbot"    long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3"       /> 
     523         <field id="ibgfrctemtop"    long_name="global mean heat on top of ice/snw/ocean-atm "                                             unit="1e20J"     /> 
     524         <field id="ibgfrctembot"    long_name="global mean heat below ice (on top of ocean) "                                             unit="1e20J"     /> 
     525         <field id="ibgfrcsal"       long_name="global mean ice/snow forcing (salt equivalent ocean volume)"                               unit="pss*km3"   /> 
     526         <field id="ibgfrchfxtop"    long_name="global mean heat flux on top of ice/snw/ocean-atm "                                        unit="W/m2"      /> 
     527         <field id="ibgfrchfxbot"    long_name="global mean heat flux below ice (on top of ocean) "                                        unit="W/m2"      /> 
     528  
     529         <field id="ibgvolume"    long_name="drift in ice/snow volume (equivalent ocean volume)"            unit="km3"        /> 
     530         <field id="ibgsaltco"    long_name="drift in ice salt content (equivalent ocean volume)"           unit="pss*km3"    /> 
     531         <field id="ibgheatco"    long_name="drift in ice/snow heat content"                                unit="1e20J"      /> 
     532         <field id="ibgheatfx"    long_name="drift in ice/snow heat flux"                                   unit="W/m2"       /> 
     533 
     534         <field id="ibgvol_tot"    long_name="global mean ice volume"                                       unit="km3"        /> 
     535         <field id="sbgvol_tot"    long_name="global mean snow volume"                                      unit="km3"        /> 
     536         <field id="ibgarea_tot"   long_name="global mean ice area"                                         unit="km2"        /> 
     537         <field id="ibgsalt_tot"   long_name="global mean ice salt content"                                 unit="1e-3*km3"   /> 
     538         <field id="ibgheat_tot"   long_name="global mean ice heat content"                                 unit="1e20J"      /> 
     539         <field id="sbgheat_tot"   long_name="global mean snow heat content"                                unit="1e20J"      /> 
    570540      </field_group> 
    571541   
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/ice.F90

    r6477 r7067  
    243243   ! 
    244244   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   sist        !: Average Sea-Ice Surface Temperature [Kelvin] 
    245    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   icethi      !: total ice thickness (for all categories) (diag only) 
    246245   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   t_bo        !: Sea-Ice bottom temperature [Kelvin]      
    247246   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   frld        !: Leads fraction = 1 - ice fraction 
     
    320319   !                                                                  !  this is an extensive variable that has to be transported 
    321320   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   o_i     !: Sea-Ice Age (days) 
    322    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   ov_i    !: Sea-Ice Age times volume per area (days.m) 
    323321   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   oa_i    !: Sea-Ice Age times ice area (days) 
     322   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) ::   bv_i    !: brine volume 
    324323 
    325324   !! Variables summed over all categories, or associated to all the ice in a single grid cell 
    326325   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   u_ice, v_ice   !: components of the ice velocity (m/s) 
    327    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tio_u, tio_v   !: components of the ice-ocean stress (N/m2) 
    328326   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   vt_i , vt_s    !: ice and snow total volume per unit area (m) 
    329327   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   at_i           !: ice total fractional area (ice concentration) 
    330328   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ato_i          !: =1-at_i ; total open water fractional area 
    331329   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   et_i , et_s    !: ice and snow total heat content 
    332    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   ot_i           !: mean age over all categories 
    333    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i           !: mean ice temperature over all categories 
    334    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bv_i           !: brine volume averaged over all categories 
    335    REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i          !: mean sea ice salinity averaged over all categories [PSU] 
     330   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_i         !: mean ice temperature over all categories 
     331   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   bvm_i        !: brine volume averaged over all categories 
     332   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   smt_i        !: mean sea ice salinity averaged over all categories [PSU] 
     333   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   tm_su        !: mean surface temperature over all categories 
     334   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_i        !: mean ice  thickness over all categories 
     335   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   htm_s        !: mean snow thickness over all categories 
     336   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) ::   om_i         !: mean ice age over all categories 
    336337 
    337338   REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) ::   t_s        !: Snow temperatures [K] 
     
    435436 
    436437      ii = ii + 1 
    437       ALLOCATE( sist   (jpi,jpj) , icethi (jpi,jpj) , t_bo   (jpi,jpj) ,                        & 
     438      ALLOCATE( sist   (jpi,jpj) , t_bo   (jpi,jpj) ,                        & 
    438439         &      frld   (jpi,jpj) , pfrld  (jpi,jpj) , phicif (jpi,jpj) ,                        & 
    439440         &      wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) ,                        & 
     
    456457         &      v_s  (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) ,     & 
    457458         &      sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i  (jpi,jpj,jpl) ,     & 
    458          &      ov_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl)                      , STAT=ierr(ii) ) 
    459       ii = ii + 1 
    460       ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) ,     & 
     459         &      oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) ) 
     460      ii = ii + 1 
     461      ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) ,      & 
    461462         &      vt_i (jpi,jpj) , vt_s (jpi,jpj) , at_i (jpi,jpj) , ato_i(jpi,jpj) ,     & 
    462          &      et_i (jpi,jpj) , et_s (jpi,jpj) , ot_i (jpi,jpj) , tm_i (jpi,jpj) ,     & 
    463          &      bv_i (jpi,jpj) , smt_i(jpi,jpj)                                   , STAT=ierr(ii) ) 
     463         &      et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) ,     & 
     464         &      smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) ,     & 
     465         &      om_i (jpi,jpj)                              , STAT=ierr(ii) ) 
    464466      ii = ii + 1 
    465467      ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90

    r6399 r7067  
    288288#if ! defined key_bdy 
    289289      ! heat flux 
    290       zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es - SUM( qevap_ice * a_i_b, dim=3 ) )  & 
    291          &              * e12t * tmask(:,:,1) * zconv )  
     290      zhfx  = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es  & 
     291      !  &              - SUM( qevap_ice * a_i_b, dim=3 )                           & !!clem: I think this line must be commented (but need check) 
     292         &              ) * e12t * tmask(:,:,1) * zconv )  
    292293      ! salt flux 
    293294      zsfx  = glob_sum( ( sfx + diag_smvi ) * e12t * tmask(:,:,1) * zconv ) * rday 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90

    r6417 r7067  
    3131 
    3232   PUBLIC   lim_diahsb        ! routine called by ice_step.F90 
    33  
    34    real(wp) ::   frc_sal, frc_vol   ! global forcing trends 
    35    real(wp) ::   bg_grme            ! global ice growth+melt trends 
    36  
     33   PUBLIC   lim_diahsb_init   ! routine called in sbcice_lim.F90 
     34 
     35   REAL(wp), DIMENSION(:,:), ALLOCATABLE ::   vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents 
     36   REAL(wp)                              ::   frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot  ! global forcing trends 
     37    
    3738   !! * Substitutions 
    3839#  include "vectopt_loop_substitute.h90" 
     
    4647CONTAINS 
    4748 
    48    SUBROUTINE lim_diahsb 
     49   SUBROUTINE lim_diahsb( kt ) 
    4950      !!--------------------------------------------------------------------------- 
    5051      !!                  ***  ROUTINE lim_diahsb  *** 
     
    5354      !!  
    5455      !!--------------------------------------------------------------------------- 
     56      INTEGER, INTENT(in) :: kt    ! number of iteration 
    5557      !! 
    56       real(wp)   ::   zbg_ivo, zbg_svo, zbg_are, zbg_sal ,zbg_tem ,zbg_ihc ,zbg_shc 
    57       real(wp)   ::   zbg_sfx, zbg_sfx_bri, zbg_sfx_bog, zbg_sfx_bom, zbg_sfx_sum, zbg_sfx_sni,   & 
    58       &               zbg_sfx_opw, zbg_sfx_res, zbg_sfx_dyn, zbg_sfx_sub  
    59       real(wp)   ::   zbg_vfx, zbg_vfx_bog, zbg_vfx_opw, zbg_vfx_sni, zbg_vfx_dyn 
    60       real(wp)   ::   zbg_vfx_bom, zbg_vfx_sum, zbg_vfx_res, zbg_vfx_spr, zbg_vfx_snw, zbg_vfx_sub   
    61       real(wp)   ::   zbg_hfx_dhc, zbg_hfx_spr 
    62       real(wp)   ::   zbg_hfx_res, zbg_hfx_sub, zbg_hfx_dyn, zbg_hfx_thd, zbg_hfx_snw, zbg_hfx_out, zbg_hfx_in    
    63       real(wp)   ::   zbg_hfx_sum, zbg_hfx_bom, zbg_hfx_bog, zbg_hfx_dif, zbg_hfx_opw 
    64       real(wp)   ::   z_frc_vol, z_frc_sal, z_bg_grme  
    65       real(wp)   ::   z1_area                     !    -     - 
    66       REAL(wp)   ::   ztmp 
     58      real(wp)   ::   zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem 
     59      REAL(wp)   ::   z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot   
     60      REAL(wp)   ::   zdiff_vol, zdiff_sal, zdiff_tem   
    6761      !!--------------------------------------------------------------------------- 
    6862      IF( nn_timing == 1 )   CALL timing_start('lim_diahsb') 
    6963 
    70       IF( numit == nstart ) CALL lim_diahsb_init  
    71  
    72       ! 1/area 
    73       z1_area = 1._wp / MAX( glob_sum( e12t(:,:) * tmask(:,:,1) ), epsi06 ) 
    74  
    75       rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e12t(:,:) * tmask(:,:,1) ) - epsi06 ) ) 
    76       ! ----------------------- ! 
    77       ! 1 -  Content variations ! 
    78       ! ----------------------- ! 
    79       zbg_ivo = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume ice  
    80       zbg_svo = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume snow 
    81       zbg_are = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) ) ! area 
    82       zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) )       ! mean salt content 
    83       zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e12t(:,:) * tmask(:,:,1) )  ! mean temp content 
    84  
    85       !zbg_ihc = glob_sum( et_i(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 
    86       !zbg_shc = glob_sum( et_s(:,:) * e12t(:,:) * tmask(:,:,1) ) / MAX( zbg_svo,epsi06 ) ! snow heat content 
    87  
    88       ! Volume 
    89       ztmp = rswitch * z1_area * r1_rau0 * rday 
    90       zbg_vfx     = ztmp * glob_sum(     emp(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    91       zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    92       zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    93       zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    94       zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    95       zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    96       zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    97       zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    98       zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    99       zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    100       zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    101  
    102       ! Salt 
    103       zbg_sfx     = ztmp * glob_sum(     sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    104       zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    105       zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    106       zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    107  
    108       zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    109       zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    110       zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    111       zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    112       zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    113       zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) 
    114  
    115       ! Heat budget 
    116       zbg_ihc      = glob_sum( et_i(:,:) * e12t(:,:) * 1.e-20 ) ! ice heat content  [1.e20 J] 
    117       zbg_shc      = glob_sum( et_s(:,:) * e12t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 
    118       zbg_hfx_dhc  = glob_sum( diag_heat(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    119       zbg_hfx_spr  = glob_sum( hfx_spr(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    120  
    121       zbg_hfx_thd  = glob_sum( hfx_thd(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    122       zbg_hfx_dyn  = glob_sum( hfx_dyn(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    123       zbg_hfx_res  = glob_sum( hfx_res(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    124       zbg_hfx_sub  = glob_sum( hfx_sub(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    125       zbg_hfx_snw  = glob_sum( hfx_snw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    126       zbg_hfx_sum  = glob_sum( hfx_sum(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    127       zbg_hfx_bom  = glob_sum( hfx_bom(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    128       zbg_hfx_bog  = glob_sum( hfx_bog(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    129       zbg_hfx_dif  = glob_sum( hfx_dif(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    130       zbg_hfx_opw  = glob_sum( hfx_opw(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    131       zbg_hfx_out  = glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    132       zbg_hfx_in   = glob_sum(  hfx_in(:,:) * e12t(:,:) * tmask(:,:,1) ) ! [in W] 
    133      
    134       ! --------------------------------------------- ! 
    135       ! 2 - Trends due to forcing and ice growth/melt ! 
    136       ! --------------------------------------------- ! 
    137       z_frc_vol = r1_rau0 * glob_sum( - emp(:,:) * e12t(:,:) * tmask(:,:,1) ) ! volume fluxes 
    138       z_frc_sal = r1_rau0 * glob_sum(   sfx(:,:) * e12t(:,:) * tmask(:,:,1) ) ! salt fluxes 
    139       z_bg_grme = glob_sum( - ( wfx_bog(:,:) + wfx_opw(:,:) + wfx_sni(:,:) + wfx_dyn(:,:) + & 
    140                           &     wfx_bom(:,:) + wfx_sum(:,:) + wfx_res(:,:) + wfx_snw(:,:) + & 
    141                           &     wfx_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) ) ! 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 
     64      ! ----------------------- ! 
     65      ! 1 -  Contents ! 
     66      ! ----------------------- ! 
     67      zbg_ivol = glob_sum( vt_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! ice volume (km3) 
     68      zbg_svol = glob_sum( vt_s(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                  ! snow volume (km3) 
     69      zbg_area = glob_sum( at_i(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-6 )                  ! area (km2) 
     70      zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3) 
     71      zbg_item = glob_sum( et_i * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J) 
     72      zbg_stem = glob_sum( et_s * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                      ! heat content (1.e20 J) 
    14673       
    147       ! difference 
    148       !frc_vol = zbg_ivo - frc_vol 
    149       !frc_sal = zbg_sal - frc_sal 
    150        
    151       ! ----------------------- ! 
    152       ! 3 - Diagnostics writing ! 
    153       ! ----------------------- ! 
    154       rswitch = MAX( 0._wp , SIGN( 1._wp , zbg_ivo - epsi06 ) ) 
    155       ! 
    156       IF( iom_use('ibgvoltot') )   & 
    157       CALL iom_put( 'ibgvoltot' , zbg_ivo * rhoic * r1_rau0 * 1.e-9        )   ! ice volume (km3 equivalent liquid)          
    158       IF( iom_use('sbgvoltot') )   & 
    159       CALL iom_put( 'sbgvoltot' , zbg_svo * rhosn * r1_rau0 * 1.e-9        )   ! snw volume (km3 equivalent liquid)        
    160       IF( iom_use('ibgarea') )   & 
    161       CALL iom_put( 'ibgarea'   , zbg_are * 1.e-6                          )   ! ice area   (km2) 
    162       IF( iom_use('ibgsaline') )   & 
    163       CALL iom_put( 'ibgsaline' , rswitch * zbg_sal / MAX( zbg_ivo, epsi06 ) )   ! ice saline (psu) 
    164       IF( iom_use('ibgtemper') )   & 
    165       CALL iom_put( 'ibgtemper' , rswitch * zbg_tem / MAX( zbg_ivo, epsi06 ) )   ! ice temper (C) 
    166       CALL iom_put( 'ibgheatco' , zbg_ihc                                  )   ! ice heat content (1.e20 J)         
    167       CALL iom_put( 'sbgheatco' , zbg_shc                                  )   ! snw heat content (1.e20 J) 
    168       IF( iom_use('ibgsaltco') )   & 
    169       CALL iom_put( 'ibgsaltco' , zbg_sal * rhoic * r1_rau0 * 1.e-9        )   ! ice salt content (psu*km3 equivalent liquid)         
    170  
    171       CALL iom_put( 'ibgvfx'    , zbg_vfx                                  )   ! volume flux emp (m/day liquid) 
    172       CALL iom_put( 'ibgvfxbog' , zbg_vfx_bog                              )   ! volume flux bottom growth     -(m/day equivalent liquid) 
    173       CALL iom_put( 'ibgvfxopw' , zbg_vfx_opw                              )   ! volume flux open water growth - 
    174       CALL iom_put( 'ibgvfxsni' , zbg_vfx_sni                              )   ! volume flux snow ice growth   - 
    175       CALL iom_put( 'ibgvfxdyn' , zbg_vfx_dyn                              )   ! volume flux dynamic growth    - 
    176       CALL iom_put( 'ibgvfxbom' , zbg_vfx_bom                              )   ! volume flux bottom melt       - 
    177       CALL iom_put( 'ibgvfxsum' , zbg_vfx_sum                              )   ! volume flux surface melt      - 
    178       CALL iom_put( 'ibgvfxres' , zbg_vfx_res                              )   ! volume flux resultant         - 
    179       CALL iom_put( 'ibgvfxspr' , zbg_vfx_spr                              )   ! volume flux from snow precip         - 
    180       CALL iom_put( 'ibgvfxsnw' , zbg_vfx_snw                              )   ! volume flux from snow melt         - 
    181       CALL iom_put( 'ibgvfxsub' , zbg_vfx_sub                              )   ! volume flux from sublimation         - 
    182            
    183       CALL iom_put( 'ibgsfx'    , zbg_sfx                                  )   ! salt flux         -(psu*m/day equivalent liquid)        
    184       CALL iom_put( 'ibgsfxbri' , zbg_sfx_bri                              )   ! salt flux brines  -       
    185       CALL iom_put( 'ibgsfxdyn' , zbg_sfx_dyn                              )   ! salt flux dynamic -     
    186       CALL iom_put( 'ibgsfxres' , zbg_sfx_res                              )   ! salt flux result  -     
    187       CALL iom_put( 'ibgsfxbog' , zbg_sfx_bog                              )   ! salt flux bottom growth    
    188       CALL iom_put( 'ibgsfxopw' , zbg_sfx_opw                              )   ! salt flux open water growth - 
    189       CALL iom_put( 'ibgsfxsni' , zbg_sfx_sni                              )   ! salt flux snow ice growth   - 
    190       CALL iom_put( 'ibgsfxbom' , zbg_sfx_bom                              )   ! salt flux bottom melt       - 
    191       CALL iom_put( 'ibgsfxsum' , zbg_sfx_sum                              )   ! salt flux surface melt      - 
    192       CALL iom_put( 'ibgsfxsub' , zbg_sfx_sub                              )   ! salt flux sublimation      - 
    193  
    194       CALL iom_put( 'ibghfxdhc' , zbg_hfx_dhc                              )   ! Heat content variation in snow and ice [W] 
    195       CALL iom_put( 'ibghfxspr' , zbg_hfx_spr                              )   ! Heat content of snow precip [W] 
    196  
    197       CALL iom_put( 'ibghfxres' , zbg_hfx_res                              )   !  
    198       CALL iom_put( 'ibghfxsub' , zbg_hfx_sub                              )   !  
    199       CALL iom_put( 'ibghfxdyn' , zbg_hfx_dyn                              )   !  
    200       CALL iom_put( 'ibghfxthd' , zbg_hfx_thd                              )   !  
    201       CALL iom_put( 'ibghfxsnw' , zbg_hfx_snw                              )   !  
    202       CALL iom_put( 'ibghfxsum' , zbg_hfx_sum                              )   !  
    203       CALL iom_put( 'ibghfxbom' , zbg_hfx_bom                              )   !  
    204       CALL iom_put( 'ibghfxbog' , zbg_hfx_bog                              )   !  
    205       CALL iom_put( 'ibghfxdif' , zbg_hfx_dif                              )   !  
    206       CALL iom_put( 'ibghfxopw' , zbg_hfx_opw                              )   !  
    207       CALL iom_put( 'ibghfxout' , zbg_hfx_out                              )   !  
    208       CALL iom_put( 'ibghfxin'  , zbg_hfx_in                               )   !  
    209  
    210       CALL iom_put( 'ibgfrcvol' , frc_vol * 1.e-9                          )   ! vol - forcing     (km3 equivalent liquid)  
    211       CALL iom_put( 'ibgfrcsfx' , frc_sal * 1.e-9                          )   ! sal - forcing     (psu*km3 equivalent liquid)    
    212       IF( iom_use('ibgvolgrm') )   & 
    213       CALL iom_put( 'ibgvolgrm' , bg_grme * r1_rau0 * 1.e-9                )   ! vol growth + melt (km3 equivalent liquid)          
    214  
     74      ! ---------------------------! 
     75      ! 2 - Trends due to forcing  ! 
     76      ! ---------------------------! 
     77      z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )  ! freshwater flux ice/snow-ocean  
     78      z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                     ! freshwater flux ice/snow-atm 
     79      z_frc_sal    = r1_rau0 * glob_sum( - sfx(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )                                            ! salt fluxes ice/snow-ocean 
     80      z_frc_tembot =           glob_sum( hfx_out(:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ocean (and below ice) 
     81      z_frc_temtop =           glob_sum( hfx_in (:,:) * e12t(:,:) * tmask(:,:,1) * 1.e-20 )                                         ! heat on top of ice-coean 
     82      ! 
     83      frc_voltop  = frc_voltop  + z_frc_voltop  * rdt_ice ! km3 
     84      frc_volbot  = frc_volbot  + z_frc_volbot  * rdt_ice ! km3 
     85      frc_sal     = frc_sal     + z_frc_sal     * rdt_ice ! km3*pss 
     86      frc_temtop  = frc_temtop  + z_frc_temtop  * rdt_ice ! 1.e20 J 
     87      frc_tembot  = frc_tembot  + z_frc_tembot  * rdt_ice ! 1.e20 J 
     88             
     89      ! ----------------------- ! 
     90      ! 3 -  Content variations ! 
     91      ! ----------------------- ! 
     92      zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:)  &  ! freshwater trend (km3)  
     93         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 )  
     94      zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:)     &  ! salt content trend (km3*pss) 
     95         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-9 ) 
     96      zdiff_tem =           glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:)                  &  ! heat content trend (1.e20 J) 
     97      !  &                            + SUM( qevap_ice * a_i_b, dim=3 ) &     !! clem: I think this line should be commented (but needs a check) 
     98         &                            ) * e12t(:,:) * tmask(:,:,1) * 1.e-20 ) 
     99 
     100      ! ----------------------- ! 
     101      ! 4 -  Drifts             ! 
     102      ! ----------------------- ! 
     103      zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot ) 
     104      zdiff_sal = zdiff_sal - frc_sal 
     105      zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop ) 
     106 
     107      ! ----------------------- ! 
     108      ! 5 - Diagnostics writing ! 
     109      ! ----------------------- ! 
     110      ! 
     111      IF( iom_use('ibgvolume') )  CALL iom_put( 'ibgvolume' , zdiff_vol        )   ! ice/snow volume  drift            (km3 equivalent ocean water)          
     112      IF( iom_use('ibgsaltco') )  CALL iom_put( 'ibgsaltco' , zdiff_sal        )   ! ice salt content drift            (psu*km3 equivalent ocean water) 
     113      IF( iom_use('ibgheatco') )  CALL iom_put( 'ibgheatco' , zdiff_tem        )   ! ice/snow heat content drift       (1.e20 J) 
     114      IF( iom_use('ibgheatfx') )  CALL iom_put( 'ibgheatfx' , zdiff_tem /      &   ! ice/snow heat flux drift          (W/m2) 
     115         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
     116 
     117      IF( iom_use('ibgfrcvoltop') )  CALL iom_put( 'ibgfrcvoltop' , frc_voltop )   ! vol  forcing ice/snw-atm          (km3 equivalent ocean water)  
     118      IF( iom_use('ibgfrcvolbot') )  CALL iom_put( 'ibgfrcvolbot' , frc_volbot )   ! vol  forcing ice/snw-ocean        (km3 equivalent ocean water)  
     119      IF( iom_use('ibgfrcsal') )     CALL iom_put( 'ibgfrcsal'    , frc_sal    )   ! sal - forcing                     (psu*km3 equivalent ocean water)    
     120      IF( iom_use('ibgfrctemtop') )  CALL iom_put( 'ibgfrctemtop' , frc_temtop )   ! heat on top of ice/snw/ocean      (1.e20 J)    
     121      IF( iom_use('ibgfrctembot') )  CALL iom_put( 'ibgfrctembot' , frc_tembot )   ! heat on top of ocean(below ice)   (1.e20 J)    
     122      IF( iom_use('ibgfrchfxtop') )  CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean      (W/m2)  
     123         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
     124      IF( iom_use('ibgfrchfxbot') )  CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice)   (W/m2)  
     125         &                                                    glob_sum( e12t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 
     126 
     127      IF( iom_use('ibgvol_tot' ) )  CALL iom_put( 'ibgvol_tot'  , zbg_ivol     )   ! ice volume                        (km3) 
     128      IF( iom_use('sbgvol_tot' ) )  CALL iom_put( 'sbgvol_tot'  , zbg_svol     )   ! snow volume                       (km3) 
     129      IF( iom_use('ibgarea_tot') )  CALL iom_put( 'ibgarea_tot' , zbg_area     )   ! ice area                          (km2) 
     130      IF( iom_use('ibgsalt_tot') )  CALL iom_put( 'ibgsalt_tot' , zbg_isal     )   ! ice salinity content              (pss*km3) 
     131      IF( iom_use('ibgheat_tot') )  CALL iom_put( 'ibgheat_tot' , zbg_item     )   ! ice heat content                  (1.e20 J) 
     132      IF( iom_use('sbgheat_tot') )  CALL iom_put( 'sbgheat_tot' , zbg_stem     )   ! snow heat content                 (1.e20 J) 
    215133      ! 
    216134      IF( lrst_ice )   CALL lim_diahsb_rst( numit, 'WRITE' ) 
    217135      ! 
    218136      IF( nn_timing == 1 )   CALL timing_stop('lim_diahsb') 
    219 ! 
     137      ! 
    220138   END SUBROUTINE lim_diahsb 
    221139 
     
    233151      !!             - Compute coefficients for conversion 
    234152      !!--------------------------------------------------------------------------- 
    235       INTEGER            ::   jk       ! dummy loop indice 
    236153      INTEGER            ::   ierror   ! local integer 
    237154      !! 
     
    247164         WRITE(numout,*) '~~~~~~~~~~~~' 
    248165      ENDIF 
    249       ! 
     166      !       
     167      ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror ) 
     168      IF( ierror > 0 )  THEN 
     169         CALL ctl_stop( 'lim_diahsb: unable to allocate vol_loc_ini' ) 
     170         RETURN 
     171      ENDIF 
     172 
    250173      CALL lim_diahsb_rst( nstart, 'READ' )  !* read or initialize all required files 
    251174      ! 
     
    263186     CHARACTER(len=*), INTENT(in) ::   cdrw   ! "READ"/"WRITE" flag 
    264187     ! 
    265      INTEGER ::   id1, id2, id3   ! local integers 
    266188     !!---------------------------------------------------------------------- 
    267189     ! 
    268190     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    269191        IF( ln_rstart ) THEN                   !* Read the restart file 
    270            !id1 = iom_varid( numrir, 'frc_vol'  , ldstop = .TRUE. ) 
    271192           ! 
    272193           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    273            IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp 
    274            IF(lwp) WRITE(numout,*) '~~~~~~~' 
    275            CALL iom_get( numrir, 'frc_vol', frc_vol ) 
    276            CALL iom_get( numrir, 'frc_sal', frc_sal ) 
    277            CALL iom_get( numrir, 'bg_grme', bg_grme ) 
     194           IF(lwp) WRITE(numout,*) ' lim_diahsb_rst read at it= ', kt,' date= ', ndastp 
     195           IF(lwp) WRITE(numout,*) '~~~~~~~' 
     196           CALL iom_get( numrir, 'frc_voltop' , frc_voltop  ) 
     197           CALL iom_get( numrir, 'frc_volbot' , frc_volbot  ) 
     198           CALL iom_get( numrir, 'frc_temtop' , frc_temtop  ) 
     199           CALL iom_get( numrir, 'frc_tembot' , frc_tembot  ) 
     200           CALL iom_get( numrir, 'frc_sal'    , frc_sal     ) 
     201           CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini ) 
     202           CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini ) 
     203           CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini ) 
    278204        ELSE 
    279205           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    280206           IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state ' 
    281207           IF(lwp) WRITE(numout,*) '~~~~~~~' 
    282            frc_vol  = 0._wp                                           
    283            frc_sal  = 0._wp                                                  
    284            bg_grme  = 0._wp                                        
     208           ! set trends to 0 
     209           frc_voltop  = 0._wp                                           
     210           frc_volbot  = 0._wp                                           
     211           frc_temtop  = 0._wp                                                  
     212           frc_tembot  = 0._wp                                                  
     213           frc_sal     = 0._wp                                                  
     214           ! record initial ice volume, salt and temp 
     215           vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:)  ! ice/snow volume (kg/m2) 
     216           tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:)                  ! ice/snow heat content (J) 
     217           sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 )     ! ice salt content (pss*kg/m2) 
     218            
    285219       ENDIF 
    286220 
     
    288222        !                                   ! ------------------- 
    289223        IF(lwp) WRITE(numout,*) '~~~~~~~' 
    290         IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp 
     224        IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write at it= ', kt,' date= ', ndastp 
    291225        IF(lwp) WRITE(numout,*) '~~~~~~~' 
    292         CALL iom_rstput( kt, nitrst, numriw, 'frc_vol'   , frc_vol     ) 
    293         CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'   , frc_sal     ) 
    294         CALL iom_rstput( kt, nitrst, numriw, 'bg_grme'   , bg_grme     ) 
     226        CALL iom_rstput( kt, nitrst, numriw, 'frc_voltop' , frc_voltop  ) 
     227        CALL iom_rstput( kt, nitrst, numriw, 'frc_volbot' , frc_volbot  ) 
     228        CALL iom_rstput( kt, nitrst, numriw, 'frc_temtop' , frc_temtop  ) 
     229        CALL iom_rstput( kt, nitrst, numriw, 'frc_tembot' , frc_tembot  ) 
     230        CALL iom_rstput( kt, nitrst, numriw, 'frc_sal'    , frc_sal     ) 
     231        CALL iom_rstput( kt, nitrst, numriw, 'vol_loc_ini', vol_loc_ini ) 
     232        CALL iom_rstput( kt, nitrst, numriw, 'tem_loc_ini', tem_loc_ini ) 
     233        CALL iom_rstput( kt, nitrst, numriw, 'sal_loc_ini', sal_loc_ini ) 
    295234        ! 
    296235     ENDIF 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90

    r6469 r7067  
    866866         DO jj = 1, jpj 
    867867            DO ji = 1, jpi 
    868                strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv_i(ji,jj),0.0))) 
     868               strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0))) 
    869869            END DO 
    870870         END DO 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90

    r5888 r7067  
    1010   !!            3.4  !  2011-01  (A. Porter)  dynamical allocation  
    1111   !!            3.5  !  2012-08  (R. Benshila)  AGRIF  
     12   !!            3.6  !  2016-06  (C. Rousset) Rewriting (conserves energy) 
    1213   !!---------------------------------------------------------------------- 
    1314#if defined key_lim3 || (  defined key_lim2 && ! defined key_lim2_vp ) 
     
    9596      !!                 coriolis terms of the momentum equation 
    9697      !!              3) Solve the momentum equation (iterative procedure) 
    97       !!              4) Prevent high velocities if the ice is thin 
    98       !!              5) Recompute invariants of the strain rate tensor 
     98      !!              4) Recompute invariants of the strain rate tensor 
    9999      !!                 which are inputs of the ITD, store stress 
    100100      !!                 for the next time step 
    101       !!              6) Control prints of residual (convergence) 
     101      !!              5) Control prints of residual (convergence) 
    102102      !!                 and charge ellipse. 
    103103      !!                 The user should make sure that the parameters 
     
    106106      !!                 e.g. in the Canadian Archipelago 
    107107      !! 
     108      !! ** Notes   : Boundary condition for ice is chosen no-slip  
     109      !!              but can be adjusted with param rn_shlat 
     110      !! 
    108111      !! References : Hunke and Dukowicz, JPO97 
    109112      !!              Bouillon et al., Ocean Modelling 2009 
     
    115118      INTEGER ::   jter     ! local integers 
    116119      CHARACTER (len=50) ::   charout 
    117       REAL(wp) ::   zt11, zt12, zt21, zt22, ztagnx, ztagny, delta                         ! 
    118       REAL(wp) ::   za, zstms          ! local scalars 
    119       REAL(wp) ::   zc1, zc2, zc3      ! ice mass 
    120  
    121       REAL(wp) ::   dtevp , z1_dtevp              ! time step for subcycling 
    122       REAL(wp) ::   dtotel, z1_dtotel, ecc2, ecci ! square of yield ellipse eccenticity 
    123       REAL(wp) ::   z0, zr, zcca, zccb            ! temporary scalars 
    124       REAL(wp) ::   zu_ice2, zv_ice1              ! 
    125       REAL(wp) ::   zddc, zdtc                    ! delta on corners and on centre 
    126       REAL(wp) ::   zdst                          ! shear at the center of the grid point 
    127       REAL(wp) ::   zdsshx, zdsshy                ! term for the gradient of ocean surface 
    128       REAL(wp) ::   sigma1, sigma2                ! internal ice stress 
    129  
    130       REAL(wp) ::   zresm         ! Maximal error on ice velocity 
    131       REAL(wp) ::   zintb, zintn  ! dummy argument 
    132  
    133       REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh           ! temporary array for ice strength 
    134       REAL(wp), POINTER, DIMENSION(:,:) ::   zpreshc          ! Ice strength on grid cell corners (zpreshc) 
    135       REAL(wp), POINTER, DIMENSION(:,:) ::   zfrld1, zfrld2   ! lead fraction on U/V points 
    136       REAL(wp), POINTER, DIMENSION(:,:) ::   zmass1, zmass2   ! ice/snow mass on U/V points 
    137       REAL(wp), POINTER, DIMENSION(:,:) ::   zcorl1, zcorl2   ! coriolis parameter on U/V points 
    138       REAL(wp), POINTER, DIMENSION(:,:) ::   za1ct , za2ct    ! temporary arrays 
    139       REAL(wp), POINTER, DIMENSION(:,:) ::   v_oce1           ! ocean u/v component on U points                            
    140       REAL(wp), POINTER, DIMENSION(:,:) ::   u_oce2           ! ocean u/v component on V points 
    141       REAL(wp), POINTER, DIMENSION(:,:) ::   u_ice2, v_ice1   ! ice u/v component on V/U point 
    142       REAL(wp), POINTER, DIMENSION(:,:) ::   zf1   , zf2      ! arrays for internal stresses 
    143       REAL(wp), POINTER, DIMENSION(:,:) ::   zmask            ! mask ocean grid points 
     120 
     121      REAL(wp) ::   zdtevp, z1_dtevp                                         ! time step for subcycling 
     122      REAL(wp) ::   ecc2, z1_ecc2                                            ! square of yield ellipse eccenticity 
     123      REAL(wp) ::   zbeta, zalph1, z1_alph1, zalph2, z1_alph2                ! alpha and beta from Bouillon 2009 and 2013 
     124      REAL(wp) ::   zm1, zm2, zm3, zmassU, zmassV                            ! ice/snow mass 
     125      REAL(wp) ::   zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2            ! temporary scalars 
     126      REAL(wp) ::   zTauO, zTauE, zCor                                       ! temporary scalars 
     127 
     128      REAL(wp) ::   zsig1, zsig2                                             ! internal ice stress 
     129      REAL(wp) ::   zresm                                                    ! Maximal error on ice velocity 
     130      REAL(wp) ::   zintb, zintn                                             ! dummy argument 
    144131       
    145       REAL(wp), POINTER, DIMENSION(:,:) ::   zdt              ! tension at centre of grid cells 
    146       REAL(wp), POINTER, DIMENSION(:,:) ::   zds              ! Shear on northeast corner of grid cells 
    147       REAL(wp), POINTER, DIMENSION(:,:) ::   zs1   , zs2      ! Diagonal stress tensor components zs1 and zs2  
    148       REAL(wp), POINTER, DIMENSION(:,:) ::   zs12             ! Non-diagonal stress tensor component zs12 
    149       REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr   ! Local error on velocity 
    150       REAL(wp), POINTER, DIMENSION(:,:) ::   zpice            ! array used for the calculation of ice surface slope: 
    151                                                               !   ocean surface (ssh_m) if ice is not embedded 
    152                                                               !   ice top surface if ice is embedded    
    153  
    154       REAL(wp), PARAMETER               ::   zepsi = 1.0e-20_wp ! tolerance parameter 
    155       REAL(wp), PARAMETER               ::   zvmin = 1.0e-03_wp ! ice volume below which ice velocity equals ocean velocity 
     132      REAL(wp), POINTER, DIMENSION(:,:) ::   zpresh                          ! temporary array for ice strength 
     133      REAL(wp), POINTER, DIMENSION(:,:) ::   z1_e1t0, z1_e2t0                ! scale factors 
     134      REAL(wp), POINTER, DIMENSION(:,:) ::   zp_delt                         ! P/delta at T points 
     135      ! 
     136      REAL(wp), POINTER, DIMENSION(:,:) ::   zaU   , zaV                     ! ice fraction on U/V points 
     137      REAL(wp), POINTER, DIMENSION(:,:) ::   zmU_t, zmV_t                    ! ice/snow mass/dt on U/V points 
     138      REAL(wp), POINTER, DIMENSION(:,:) ::   zmf                             ! coriolis parameter at T points 
     139      REAL(wp), POINTER, DIMENSION(:,:) ::   zTauU_ia , ztauV_ia             ! ice-atm. stress at U-V points 
     140      REAL(wp), POINTER, DIMENSION(:,:) ::   zspgU , zspgV                   ! surface pressure gradient at U/V points 
     141      REAL(wp), POINTER, DIMENSION(:,:) ::   v_oceU, u_oceV, v_iceU, u_iceV  ! ocean/ice u/v component on V/U points                            
     142      REAL(wp), POINTER, DIMENSION(:,:) ::   zfU   , zfV                     ! internal stresses 
     143       
     144      REAL(wp), POINTER, DIMENSION(:,:) ::   zds                             ! shear 
     145      REAL(wp), POINTER, DIMENSION(:,:) ::   zs1, zs2, zs12                  ! stress tensor components 
     146      REAL(wp), POINTER, DIMENSION(:,:) ::   zu_ice, zv_ice, zresr           ! check convergence 
     147      REAL(wp), POINTER, DIMENSION(:,:) ::   zpice                           ! array used for the calculation of ice surface slope: 
     148                                                                             !   ocean surface (ssh_m) if ice is not embedded 
     149                                                                             !   ice top surface if ice is embedded    
     150      REAL(wp), POINTER, DIMENSION(:,:) ::   zswitchU, zswitchV              ! dummy arrays 
     151      REAL(wp), POINTER, DIMENSION(:,:) ::   zmaskU, zmaskV                  ! mask for ice presence 
     152      REAL(wp), POINTER, DIMENSION(:,:) ::   zfmask, zwf                     ! mask at F points for the ice 
     153 
     154      REAL(wp), PARAMETER               ::   zepsi  = 1.0e-20_wp             ! tolerance parameter 
     155      REAL(wp), PARAMETER               ::   zmmin  = 1._wp                  ! ice mass (kg/m2) below which ice velocity equals ocean velocity 
     156      REAL(wp), PARAMETER               ::   zshlat = 2._wp                  ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip) 
    156157      !!------------------------------------------------------------------- 
    157158 
    158       CALL wrk_alloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    159       CALL wrk_alloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    160       CALL wrk_alloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    161       CALL wrk_alloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
     159      CALL wrk_alloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
     160      CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     161      CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
     162      CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     163      CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    162164 
    163165#if  defined key_lim2 && ! defined key_lim2_vp 
     
    176178      ! 
    177179      !------------------------------------------------------------------------------! 
    178       ! 1) Ice strength (zpresh)                                ! 
    179       !------------------------------------------------------------------------------! 
    180       ! 
    181       ! Put every vector to 0 
    182       delta_i(:,:) = 0._wp   ; 
    183       zpresh (:,:) = 0._wp   ;   
    184       zpreshc(:,:) = 0._wp 
    185       u_ice2 (:,:) = 0._wp   ;   v_ice1(:,:) = 0._wp 
    186       divu_i (:,:) = 0._wp   ;   zdt   (:,:) = 0._wp   ;   zds(:,:) = 0._wp 
    187       shear_i(:,:) = 0._wp 
    188  
     180      ! 0) mask at F points for the ice (on the whole domain, not only k_j1,k_jpj)  
     181      !------------------------------------------------------------------------------! 
     182      ! ocean/land mask 
     183      DO jj = 1, jpjm1 
     184         DO ji = 1, jpim1      ! NO vector opt. 
     185            zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 
     186         END DO 
     187      END DO 
     188      CALL lbc_lnk( zfmask, 'F', 1._wp ) 
     189 
     190      ! Lateral boundary conditions on velocity (modify zfmask) 
     191      zwf(:,:) = zfmask(:,:) 
     192      DO jj = 2, jpjm1 
     193         DO ji = fs_2, fs_jpim1   ! vector opt. 
     194            IF( zfmask(ji,jj) == 0._wp ) THEN 
     195               zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 
     196            ENDIF 
     197         END DO 
     198      END DO 
     199      DO jj = 2, jpjm1 
     200         IF( zfmask(1,jj) == 0._wp ) THEN 
     201            zfmask(1  ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 
     202         ENDIF 
     203         IF( zfmask(jpi,jj) == 0._wp ) THEN 
     204            zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 
     205         ENDIF 
     206      END DO 
     207      DO ji = 2, jpim1 
     208         IF( zfmask(ji,1) == 0._wp ) THEN 
     209            zfmask(ji,1  ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 
     210         ENDIF 
     211         IF( zfmask(ji,jpj) == 0._wp ) THEN 
     212            zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 
     213         ENDIF 
     214      END DO 
     215      CALL lbc_lnk( zfmask, 'F', 1._wp ) 
     216 
     217      !------------------------------------------------------------------------------! 
     218      ! 1) define some variables and initialize arrays 
     219      !------------------------------------------------------------------------------! 
     220      ! ecc2: square of yield ellipse eccenticrity 
     221      ecc2    = rn_ecc * rn_ecc 
     222      z1_ecc2 = 1._wp / ecc2 
     223 
     224      ! Time step for subcycling 
     225      zdtevp   = rdt_ice / REAL( nn_nevp ) 
     226      z1_dtevp = 1._wp / zdtevp 
     227 
     228      ! alpha parameters (Bouillon 2009) 
    189229#if defined key_lim3 
    190       CALL lim_itd_me_icestrength( nn_icestr )      ! LIM-3: Ice strength on T-points 
    191 #endif 
    192  
    193       DO jj = k_j1 , k_jpj       ! Ice mass and temp variables 
    194          DO ji = 1 , jpi 
     230      zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 
     231#else 
     232      zalph1 = ( 2._wp * telast ) * z1_dtevp 
     233#endif 
     234      zalph2 = zalph1 * z1_ecc2 
     235 
     236      z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 
     237      z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 
     238 
     239      ! Initialise stress tensor  
     240      zs1 (:,:) = stress1_i (:,:)  
     241      zs2 (:,:) = stress2_i (:,:) 
     242      zs12(:,:) = stress12_i(:,:) 
     243 
     244      ! Ice strength 
    195245#if defined key_lim3 
    196             zpresh(ji,jj) = tmask(ji,jj,1) *  strength(ji,jj) 
    197 #endif 
    198 #if defined key_lim2 
    199             zpresh(ji,jj) = tmask(ji,jj,1) *  pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) ) 
    200 #endif 
    201             ! zmask = 1 where there is ice or on land 
    202             zmask(ji,jj)    = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1) 
     246      CALL lim_itd_me_icestrength( nn_icestr ) 
     247      zpresh(:,:) = tmask(:,:,1) *  strength(:,:) 
     248#else 
     249      zpresh(:,:) = tmask(:,:,1) *  pstar * vt_i(:,:) * EXP( -c_rhg * (1. - at_i(:,:) ) ) 
     250#endif 
     251 
     252      ! scale factors 
     253      DO jj = k_j1+1, k_jpj-1 
     254         DO ji = fs_2, fs_jpim1 
     255            z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj  ) + e1t(ji,jj  ) ) 
     256            z1_e2t0(ji,jj) = 1._wp / ( e2t(ji  ,jj+1) + e2t(ji,jj  ) ) 
    203257         END DO 
    204258      END DO 
    205  
    206       ! Ice strength on grid cell corners (zpreshc) 
    207       ! needed for calculation of shear stress  
    208       DO jj = k_j1+1, k_jpj-1 
    209          DO ji = 2, jpim1 !RB caution no fs_ (ji+1,jj+1) 
    210             zstms          =  tmask(ji+1,jj+1,1) * wght(ji+1,jj+1,2,2) + tmask(ji,jj+1,1) * wght(ji+1,jj+1,1,2) +   & 
    211                &              tmask(ji+1,jj,1)   * wght(ji+1,jj+1,2,1) + tmask(ji,jj,1)   * wght(ji+1,jj+1,1,1) 
    212             zpreshc(ji,jj) = ( zpresh(ji+1,jj+1) * wght(ji+1,jj+1,2,2) + zpresh(ji,jj+1) * wght(ji+1,jj+1,1,2) +   & 
    213                &               zpresh(ji+1,jj)   * wght(ji+1,jj+1,2,1) + zpresh(ji,jj)   * wght(ji+1,jj+1,1,1)     & 
    214                &             ) / MAX( zstms, zepsi ) 
    215          END DO 
    216       END DO 
    217       CALL lbc_lnk( zpreshc(:,:), 'F', 1. ) 
     259             
    218260      ! 
    219261      !------------------------------------------------------------------------------! 
    220262      ! 2) Wind / ocean stress, mass terms, coriolis terms 
    221263      !------------------------------------------------------------------------------! 
    222       ! 
    223       !  Wind stress, coriolis and mass terms on the sides of the squares         
    224       !  zfrld1: lead fraction on U-points                                       
    225       !  zfrld2: lead fraction on V-points                                      
    226       !  zmass1: ice/snow mass on U-points                                     
    227       !  zmass2: ice/snow mass on V-points                                    
    228       !  zcorl1: Coriolis parameter on U-points                              
    229       !  zcorl2: Coriolis parameter on V-points                             
    230       !  (ztagnx,ztagny): wind stress on U/V points                        
    231       !  v_oce1: ocean v component on u points                           
    232       !  u_oce2: ocean u component on v points                          
    233264 
    234265      IF( nn_ice_embd == 2 ) THEN             !== embedded sea ice: compute representative ice top surface ==! 
     
    242273         zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 
    243274         ! 
    244          zpice(:,:) = ssh_m(:,:) + (  zintn * snwice_mass(:,:) +  zintb * snwice_mass_b(:,:) ) * r1_rau0 
     275         zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 
    245276         ! 
    246277      ELSE                                    !== non-embedded sea ice: use ocean surface for slope calculation ==! 
     
    251282         DO ji = fs_2, fs_jpim1 
    252283 
    253             zc1 = tmask(ji  ,jj  ,1) * ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
    254             zc2 = tmask(ji+1,jj  ,1) * ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
    255             zc3 = tmask(ji  ,jj+1,1) * ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
    256  
    257             zt11 = tmask(ji  ,jj,1) * e1t(ji  ,jj) 
    258             zt12 = tmask(ji+1,jj,1) * e1t(ji+1,jj) 
    259             zt21 = tmask(ji,jj  ,1) * e2t(ji,jj  ) 
    260             zt22 = tmask(ji,jj+1,1) * e2t(ji,jj+1) 
    261  
    262             ! Leads area. 
    263             zfrld1(ji,jj) = ( zt12 * ( 1.0 - at_i(ji,jj) ) + zt11 * ( 1.0 - at_i(ji+1,jj) ) ) / ( zt11 + zt12 + zepsi ) 
    264             zfrld2(ji,jj) = ( zt22 * ( 1.0 - at_i(ji,jj) ) + zt21 * ( 1.0 - at_i(ji,jj+1) ) ) / ( zt21 + zt22 + zepsi ) 
    265  
    266             ! Mass, coriolis coeff. and currents 
    267             zmass1(ji,jj) = ( zt12 * zc1 + zt11 * zc2 ) / ( zt11 + zt12 + zepsi ) 
    268             zmass2(ji,jj) = ( zt22 * zc1 + zt21 * zc3 ) / ( zt21 + zt22 + zepsi ) 
    269             zcorl1(ji,jj) = zmass1(ji,jj) * ( e1t(ji+1,jj) * fcor(ji,jj) + e1t(ji,jj) * fcor(ji+1,jj) )   & 
    270                &                          / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 
    271             zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * fcor(ji,jj) + e2t(ji,jj) * fcor(ji,jj+1) )   & 
    272                &                          / ( e2t(ji,jj+1) + e2t(ji,jj) + zepsi ) 
    273             ! 
    274             ! Ocean has no slip boundary condition 
    275             v_oce1(ji,jj)  = 0.5 * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji,jj)      & 
    276                &                   + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji+1,jj) )  & 
    277                &                   / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)   
    278  
    279             u_oce2(ji,jj)  = 0.5 * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj)      & 
    280                &                   + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj+1) )  & 
    281                &                   / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    282  
    283             ! Wind stress at U,V-point 
    284             ztagnx = ( 1. - zfrld1(ji,jj) ) * utau_ice(ji,jj) 
    285             ztagny = ( 1. - zfrld2(ji,jj) ) * vtau_ice(ji,jj) 
    286  
    287             ! Computation of the velocity field taking into account the ice internal interaction. 
    288             ! Terms that are independent of the velocity field. 
    289  
    290             ! SB On utilise maintenant le gradient de la pente de l'ocean 
    291             ! include it later 
    292  
    293             zdsshx =  ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
    294             zdsshy =  ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
    295  
    296             za1ct(ji,jj) = ztagnx - zmass1(ji,jj) * grav * zdsshx 
    297             za2ct(ji,jj) = ztagny - zmass2(ji,jj) * grav * zdsshy 
     284            ! ice fraction at U-V points 
     285            zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji+1,jj) * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1) 
     286            zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e12t(ji,jj) + at_i(ji,jj+1) * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1) 
     287 
     288            ! Ice/snow mass at U-V points 
     289            zm1 = ( rhosn * vt_s(ji  ,jj  ) + rhoic * vt_i(ji  ,jj  ) ) 
     290            zm2 = ( rhosn * vt_s(ji+1,jj  ) + rhoic * vt_i(ji+1,jj  ) ) 
     291            zm3 = ( rhosn * vt_s(ji  ,jj+1) + rhoic * vt_i(ji  ,jj+1) ) 
     292            zmassU = 0.5_wp * ( zm1 * e12t(ji,jj) + zm2 * e12t(ji+1,jj) ) * r1_e12u(ji,jj) * umask(ji,jj,1) 
     293            zmassV = 0.5_wp * ( zm1 * e12t(ji,jj) + zm3 * e12t(ji,jj+1) ) * r1_e12v(ji,jj) * vmask(ji,jj,1) 
     294 
     295            ! Ocean currents at U-V points 
     296            v_oceU(ji,jj)   = 0.5_wp * ( ( v_oce(ji  ,jj) + v_oce(ji  ,jj-1) ) * e1t(ji+1,jj)    & 
     297               &                       + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     298             
     299            u_oceV(ji,jj)   = 0.5_wp * ( ( u_oce(ji,jj  ) + u_oce(ji-1,jj  ) ) * e2t(ji,jj+1)    & 
     300               &                       + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     301 
     302            ! Coriolis at T points (m*f) 
     303            zmf(ji,jj)      = zm1 * fcor(ji,jj) 
     304 
     305            ! m/dt 
     306            zmU_t(ji,jj)    = zmassU * z1_dtevp 
     307            zmV_t(ji,jj)    = zmassV * z1_dtevp 
     308 
     309            ! Drag ice-atm. 
     310            zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 
     311            zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 
     312 
     313            ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 
     314            zspgU(ji,jj)    = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 
     315            zspgV(ji,jj)    = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 
     316 
     317            ! masks 
     318            zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) )  ! 0 if no ice 
     319            zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) )  ! 0 if no ice 
     320 
     321            ! switches 
     322            zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 
     323            zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 
    298324 
    299325         END DO 
    300326      END DO 
    301  
     327      CALL lbc_lnk( zmf, 'T', 1. ) 
    302328      ! 
    303329      !------------------------------------------------------------------------------! 
     
    305331      !------------------------------------------------------------------------------! 
    306332      ! 
    307       ! Time step for subcycling 
    308       dtevp  = rdt_ice / nn_nevp 
    309 #if defined key_lim3 
    310       dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice ) 
    311 #else 
    312       dtotel = dtevp / ( 2._wp * telast ) 
    313 #endif 
    314       z1_dtotel = 1._wp / ( 1._wp + dtotel ) 
    315       z1_dtevp  = 1._wp / dtevp 
    316       !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter) 
    317       ecc2 = rn_ecc * rn_ecc 
    318       ecci = 1. / ecc2 
    319  
    320       !-Initialise stress tensor  
    321       zs1 (:,:) = stress1_i (:,:)  
    322       zs2 (:,:) = stress2_i (:,:) 
    323       zs12(:,:) = stress12_i(:,:) 
    324  
    325333      !                                               !----------------------! 
    326334      DO jter = 1 , nn_nevp                           !    loop over jter    ! 
    327335         !                                            !----------------------!         
    328          DO jj = k_j1, k_jpj-1 
    329             zu_ice(:,jj) = u_ice(:,jj)    ! velocity at previous time step 
    330             zv_ice(:,jj) = v_ice(:,jj) 
    331          END DO 
    332  
    333          DO jj = k_j1+1, k_jpj-1 
    334             DO ji = fs_2, fs_jpim1   !RB bug no vect opt due to zmask 
    335  
    336                !   
    337                !- Divergence, tension and shear (Section a. Appendix B of Hunke & Dukowicz, 2002) 
    338                !- divu_i(:,:), zdt(:,:): divergence and tension at centre of grid cells 
    339                !- zds(:,:): shear on northeast corner of grid cells 
    340                ! 
    341                !- IMPORTANT REMINDER: Dear Gurvan, note that, the way these terms are coded,  
    342                !                      there are many repeated calculations.  
    343                !                      Speed could be improved by regrouping terms. For 
    344                !                      the moment, however, the stress is on clarity of coding to avoid 
    345                !                      bugs (Martin, for Miguel). 
    346                ! 
    347                !- ALSO: arrays zdt, zds and delta could  
    348                !  be removed in the future to minimise memory demand. 
    349                ! 
    350                !- MORE NOTES: Note that we are calculating deformation rates and stresses on the corners of 
    351                !              grid cells, exactly as in the B grid case. For simplicity, the indexation on 
    352                !              the corners is the same as in the B grid. 
    353                ! 
    354                ! 
    355                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
    356                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
    357                   &            ) * r1_e12t(ji,jj) 
    358  
    359                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
    360                   &         - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
    361                   &         ) * r1_e12t(ji,jj) 
    362  
    363                ! 
     336         IF(ln_ctl) THEN   ! Convergence test 
     337            DO jj = k_j1, k_jpj-1 
     338               zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 
     339               zv_ice(:,jj) = v_ice(:,jj) 
     340            END DO 
     341         ENDIF 
     342 
     343         ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 
     344         DO jj = k_j1, k_jpj-1         ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 
     345            DO ji = 1, jpim1 
     346 
     347               ! shear at F points 
    364348               zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
    365349                  &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
    366                   &         ) * r1_e12f(ji,jj) * ( 2._wp - fmask(ji,jj,1) )   & 
    367                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    368  
    369  
    370                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
    371                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    372                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1)  
    373  
    374                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    375                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    376                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    377             END DO 
    378          END DO 
    379  
    380          CALL lbc_lnk_multi( v_ice1, 'U', -1., u_ice2, 'V', -1. )      ! lateral boundary cond. 
    381           
     350                  &         ) * r1_e12f(ji,jj) * zfmask(ji,jj) 
     351 
     352            END DO 
     353         END DO 
     354         CALL lbc_lnk( zds, 'F', 1. ) 
     355 
    382356         DO jj = k_j1+1, k_jpj-1 
    383             DO ji = fs_2, fs_jpim1 
    384  
    385                !- Calculate Delta at centre of grid cells 
    386                zdst          = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )   & 
    387                   &            + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1)   & 
    388                   &            ) * r1_e12t(ji,jj) 
    389  
    390                delta          = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    391                delta_i(ji,jj) = delta + rn_creepl 
    392  
    393                !- Calculate Delta on corners 
    394                zddc  = (  ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    395                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    396                   &    ) * r1_e12f(ji,jj) 
    397  
    398                zdtc  = (- ( v_ice1(ji,jj+1) * r1_e1u(ji,jj+1) - v_ice1(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    399                   &     + ( u_ice2(ji+1,jj) * r1_e2v(ji+1,jj) - u_ice2(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    400                   &    ) * r1_e12f(ji,jj) 
    401  
    402                zddc = SQRT( zddc**2 + ( zdtc**2 + zds(ji,jj)**2 ) * usecc2 ) + rn_creepl 
    403  
    404                !-Calculate stress tensor components zs1 and zs2 at centre of grid cells (see section 3.5 of CICE user's guide). 
    405                zs1(ji,jj)  = ( zs1 (ji,jj) + dtotel * ( divu_i(ji,jj) - delta ) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    406                   &          ) * z1_dtotel 
    407                zs2(ji,jj)  = ( zs2 (ji,jj) + dtotel *         ecci * zdt(ji,jj) / delta_i(ji,jj)   * zpresh(ji,jj)   & 
    408                   &          ) * z1_dtotel 
    409                !-Calculate stress tensor component zs12 at corners 
    410                zs12(ji,jj) = ( zs12(ji,jj) + dtotel *         ecci * zds(ji,jj) / ( 2._wp * zddc ) * zpreshc(ji,jj)  & 
    411                   &          ) * z1_dtotel  
    412  
    413             END DO 
    414          END DO 
    415  
    416          CALL lbc_lnk_multi( zs1 , 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
     357            DO ji = 2, jpim1 ! no vector loop 
     358 
     359               ! shear**2 at T points (doc eq. A16) 
     360               zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     361                  &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     362                  &   ) * 0.25_wp * r1_e12t(ji,jj) 
     363               
     364               ! divergence at T points 
     365               zdiv  = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     366                  &    + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     367                  &    ) * r1_e12t(ji,jj) 
     368               zdiv2 = zdiv * zdiv 
     369                
     370               ! tension at T points 
     371               zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     372                  &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     373                  &   ) * r1_e12t(ji,jj) 
     374               zdt2 = zdt * zdt 
     375                
     376               ! delta at T points 
     377               zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * usecc2 )   
     378 
     379               ! P/delta at T points 
     380               zp_delt(ji,jj) = zpresh(ji,jj) / ( zdelta + rn_creepl ) 
     381                
     382               ! stress at T points 
     383               zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 
     384               zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 
     385              
     386            END DO 
     387         END DO 
     388         CALL lbc_lnk( zp_delt, 'T', 1. ) 
     389 
     390         DO jj = k_j1, k_jpj-1 
     391            DO ji = 1, jpim1 
     392 
     393               ! P/delta at F points 
     394               zp_delf = 0.25_wp * ( zp_delt(ji,jj) + zp_delt(ji+1,jj) + zp_delt(ji,jj+1) + zp_delt(ji+1,jj+1) ) 
     395                
     396               ! stress at F points 
     397               zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 
     398 
     399            END DO 
     400         END DO 
     401         CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 
    417402  
    418          ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 
     403         ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 
    419404         DO jj = k_j1+1, k_jpj-1 
    420             DO ji = fs_2, fs_jpim1 
    421                !- contribution of zs1, zs2 and zs12 to zf1 
    422                zf1(ji,jj) = 0.5 * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)  & 
    423                   &             + ( zs2(ji+1,jj) * e2t(ji+1,jj)**2 - zs2(ji,jj) * e2t(ji,jj)**2 ) * r1_e2u(ji,jj)          & 
    424                   &             + 2.0 * ( zs12(ji,jj) * e1f(ji,jj)**2 - zs12(ji,jj-1) * e1f(ji,jj-1)**2 ) * r1_e1u(ji,jj)  & 
    425                   &                ) * r1_e12u(ji,jj) 
    426                ! contribution of zs1, zs2 and zs12 to zf2 
    427                zf2(ji,jj) = 0.5 * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)  & 
    428                   &             - ( zs2(ji,jj+1) * e1t(ji,jj+1)**2 - zs2(ji,jj) * e1t(ji,jj)**2 ) * r1_e1v(ji,jj)          & 
    429                   &             + 2.0 * ( zs12(ji,jj) * e2f(ji,jj)**2 - zs12(ji-1,jj) * e2f(ji-1,jj)**2 ) * r1_e2v(ji,jj)  & 
    430                   &               )  * r1_e12v(ji,jj) 
     405            DO ji = fs_2, fs_jpim1                
     406 
     407               ! U points 
     408               zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj)                                             & 
     409                  &                  + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj)    & 
     410                  &                    ) * r1_e2u(ji,jj)                                                                      & 
     411                  &                  + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1)  & 
     412                  &                    ) * 2._wp * r1_e1u(ji,jj)                                                              & 
     413                  &                  ) * r1_e12u(ji,jj) 
     414 
     415               ! V points 
     416               zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj)                                             & 
     417                  &                  - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj)    & 
     418                  &                    ) * r1_e1v(ji,jj)                                                                      & 
     419                  &                  + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj)  & 
     420                  &                    ) * 2._wp * r1_e2v(ji,jj)                                                              & 
     421                  &                  ) * r1_e12v(ji,jj) 
     422 
     423               ! u_ice at V point 
     424               u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
     425                  &                     + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 
     426                
     427               ! v_ice at U point 
     428               v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji+1,jj)     & 
     429                  &                     + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 
     430 
    431431            END DO 
    432432         END DO 
    433433         ! 
    434          ! Computation of ice velocity 
    435          ! 
    436          ! Both the Coriolis term and the ice-ocean drag are solved semi-implicitly. 
    437          ! 
    438          IF (MOD(jter,2).eq.0) THEN  
    439  
     434         ! --- Computation of ice velocity --- ! 
     435         !  Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 
     436         !  Bouillon et al. 2009 (eq 34-35) => stable 
     437         IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 
     438             
    440439            DO jj = k_j1+1, k_jpj-1 
    441440               DO ji = fs_2, fs_jpim1 
    442                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    443                   z0           = zmass1(ji,jj) * z1_dtevp 
    444  
    445                   ! SB modif because ocean has no slip boundary condition 
    446                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji  ,jj)     & 
    447                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    448                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    449                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    450                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    451                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    452                   zcca         = z0 + za 
    453                   zccb         = zcorl1(ji,jj) 
    454                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
     441 
     442                  ! tau_io/(v_oce - v_ice) 
     443                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     444                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     445 
     446                  ! Coriolis at V-points (energy conserving formulation) 
     447                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     448                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     449                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     450 
     451                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     452                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     453                   
     454                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     455                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     456                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     457                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     458                     &           ) * zmaskV(ji,jj) 
    455459               END DO 
    456460            END DO 
    457  
    458             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
     461            CALL lbc_lnk( v_ice, 'V', -1. ) 
     462             
     463#if defined key_agrif && defined key_lim2 
     464            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
     465#endif 
     466#if defined key_bdy 
     467            CALL bdy_ice_lim_dyn( 'V' ) 
     468#endif          
     469 
     470            DO jj = k_j1+1, k_jpj-1 
     471               DO ji = fs_2, fs_jpim1 
     472                                
     473                  ! tau_io/(u_oce - u_ice) 
     474                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     475                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     476 
     477                  ! Coriolis at U-points (energy conserving formulation) 
     478                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     479                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     480                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     481                   
     482                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     483                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     484 
     485                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     486                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     487                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     488                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     489                     &           ) * zmaskU(ji,jj) 
     490               END DO 
     491            END DO 
     492            CALL lbc_lnk( u_ice, 'U', -1. ) 
     493             
    459494#if defined key_agrif && defined key_lim2 
    460495            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    461496#endif 
    462497#if defined key_bdy 
    463          CALL bdy_ice_lim_dyn( 'U' ) 
     498            CALL bdy_ice_lim_dyn( 'U' ) 
    464499#endif          
     500 
     501         ELSE ! odd iterations 
    465502 
    466503            DO jj = k_j1+1, k_jpj-1 
    467504               DO ji = fs_2, fs_jpim1 
    468  
    469                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    470                   z0           = zmass2(ji,jj) * z1_dtevp 
    471                   ! SB modif because ocean has no slip boundary condition 
    472                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    473                      &                 + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    474                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    475                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  &  
    476                      &                         ( v_ice(ji,jj) - v_oce(ji,jj))**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    477                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    478                   zcca         = z0 + za 
    479                   zccb         = zcorl2(ji,jj) 
    480                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
     505                                
     506                  ! tau_io/(u_oce - u_ice) 
     507                  zTauO = zaU(ji,jj) * rhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) )  & 
     508                     &                             + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 
     509 
     510                  ! Coriolis at U-points (energy conserving formulation) 
     511                  zCor  =   0.25_wp * r1_e1u(ji,jj) *  & 
     512                     &    ( zmf(ji  ,jj) * ( e1v(ji  ,jj) * v_ice(ji  ,jj) + e1v(ji  ,jj-1) * v_ice(ji  ,jj-1) )  & 
     513                     &    + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 
     514                   
     515                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     516                  zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 
     517 
     518                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     519                  u_ice(ji,jj) = ( ( zmU_t(ji,jj) * u_ice(ji,jj) + zTauE + zTauO * u_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     520                     &             ) / MAX( zepsi, zmU_t(ji,jj) + zTauO ) * zswitchU(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     521                     &             + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin  
     522                     &           ) * zmaskU(ji,jj) 
    481523               END DO 
    482524            END DO 
    483  
    484             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
     525            CALL lbc_lnk( u_ice, 'U', -1. ) 
     526             
     527#if defined key_agrif && defined key_lim2 
     528            CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
     529#endif 
     530#if defined key_bdy 
     531            CALL bdy_ice_lim_dyn( 'U' ) 
     532#endif          
     533 
     534           DO jj = k_j1+1, k_jpj-1 
     535               DO ji = fs_2, fs_jpim1 
     536 
     537                  ! tau_io/(v_oce - v_ice) 
     538                  zTauO = zaV(ji,jj) * rhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) )  & 
     539                     &                             + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 
     540 
     541                  ! Coriolis at V-points (energy conserving formulation) 
     542                  zCor  = - 0.25_wp * r1_e2v(ji,jj) *  & 
     543                     &    ( zmf(ji,jj  ) * ( e2u(ji,jj  ) * u_ice(ji,jj  ) + e2u(ji-1,jj  ) * u_ice(ji-1,jj  ) )  & 
     544                     &    + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 
     545 
     546                  ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 
     547                  zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 
     548                   
     549                  ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 
     550                  v_ice(ji,jj) = ( ( zmV_t(ji,jj) * v_ice(ji,jj) + zTauE + zTauO * v_ice(ji,jj)  &  ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 
     551                     &             ) / MAX( zepsi, zmV_t(ji,jj) + zTauO ) * zswitchV(ji,jj)      &  ! m/dt + tau_io(only ice part) 
     552                     &             + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) )                  &  ! v_ice = v_oce if mass < zmmin 
     553                     &           ) * zmaskV(ji,jj) 
     554               END DO 
     555            END DO 
     556            CALL lbc_lnk( v_ice, 'V', -1. ) 
     557             
    485558#if defined key_agrif && defined key_lim2 
    486559            CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    487560#endif 
    488561#if defined key_bdy 
    489          CALL bdy_ice_lim_dyn( 'V' ) 
     562            CALL bdy_ice_lim_dyn( 'V' ) 
    490563#endif          
    491564 
    492          ELSE  
     565         ENDIF 
     566          
     567         IF(ln_ctl) THEN   ! Convergence test 
    493568            DO jj = k_j1+1, k_jpj-1 
    494                DO ji = fs_2, fs_jpim1 
    495                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 
    496                   z0           = zmass2(ji,jj) * z1_dtevp 
    497                   ! SB modif because ocean has no slip boundary condition 
    498                   zu_ice2      = 0.5 * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj)       & 
    499                      &                  +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) )   & 
    500                      &                 / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1)    
    501  
    502                   za           = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 +  & 
    503                      &                         ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 
    504                   zr           = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 
    505                   zcca         = z0 + za 
    506                   zccb         = zcorl2(ji,jj) 
    507                   v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 
    508                END DO 
    509             END DO 
    510  
    511             CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 
    512 #if defined key_agrif && defined key_lim2 
    513             CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 
    514 #endif 
    515 #if defined key_bdy 
    516          CALL bdy_ice_lim_dyn( 'V' ) 
    517 #endif          
    518  
    519             DO jj = k_j1+1, k_jpj-1 
    520                DO ji = fs_2, fs_jpim1 
    521                   rswitch      = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 
    522                   z0           = zmass1(ji,jj) * z1_dtevp 
    523                   zv_ice1      = 0.5 * ( ( v_ice(ji  ,jj) + v_ice(ji  ,jj-1) ) * e1t(ji,jj)       & 
    524                      &                 + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) )   & 
    525                      &                 / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    526  
    527                   za           = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 +  & 
    528                      &                         ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 
    529                   zr           = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 
    530                   zcca         = z0 + za 
    531                   zccb         = zcorl1(ji,jj) 
    532                   u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch  
    533                END DO 
    534             END DO 
    535  
    536             CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 
    537 #if defined key_agrif && defined key_lim2 
    538             CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 
    539 #endif 
    540 #if defined key_bdy 
    541          CALL bdy_ice_lim_dyn( 'U' ) 
    542 #endif          
    543  
    544          ENDIF 
    545           
    546          IF(ln_ctl) THEN 
    547             !---  Convergence test. 
    548             DO jj = k_j1+1 , k_jpj-1 
    549569               zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 
    550570            END DO 
     
    552572            IF( lk_mpp )   CALL mpp_max( zresm )   ! max over the global domain 
    553573         ENDIF 
    554  
     574         ! 
    555575         !                                                ! ==================== ! 
    556576      END DO                                              !  end loop over jter  ! 
     
    558578      ! 
    559579      !------------------------------------------------------------------------------! 
    560       ! 4) Prevent ice velocities when the ice is thin 
    561       !------------------------------------------------------------------------------! 
    562       ! If the ice volume is below zvmin then ice velocity should equal the 
    563       ! ocean velocity. This prevents high velocity when ice is thin 
    564       DO jj = k_j1+1, k_jpj-1 
    565          DO ji = fs_2, fs_jpim1 
    566             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    567                u_ice(ji,jj) = u_oce(ji,jj) 
    568                v_ice(ji,jj) = v_oce(ji,jj) 
    569             ENDIF 
     580      ! 4) Recompute delta, shear and div (inputs for mechanical redistribution)  
     581      !------------------------------------------------------------------------------! 
     582      DO jj = k_j1, k_jpj-1  
     583         DO ji = 1, jpim1 
     584 
     585            ! shear at F points 
     586            zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)   & 
     587               &         + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)   & 
     588               &         ) * r1_e12f(ji,jj) * zfmask(ji,jj) 
     589 
     590         END DO 
     591      END DO            
     592      CALL lbc_lnk( zds, 'F', 1. ) 
     593       
     594      DO jj = k_j1+1, k_jpj-1  
     595         DO ji = 2, jpim1 ! no vector loop 
     596             
     597            ! tension**2 at T points 
     598            zdt  = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)   & 
     599               &   - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)   & 
     600               &   ) * r1_e12t(ji,jj) 
     601            zdt2 = zdt * zdt 
     602             
     603            ! shear**2 at T points (doc eq. A16) 
     604            zds2 = ( zds(ji,jj  ) * zds(ji,jj  ) * e12f(ji,jj  ) + zds(ji-1,jj  ) * zds(ji-1,jj  ) * e12f(ji-1,jj  )  & 
     605               &   + zds(ji,jj-1) * zds(ji,jj-1) * e12f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e12f(ji-1,jj-1)  & 
     606               &   ) * 0.25_wp * r1_e12t(ji,jj) 
     607             
     608            ! shear at T points 
     609            shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 
     610 
     611            ! divergence at T points 
     612            divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj)   & 
     613               &            + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1)   & 
     614               &            ) * r1_e12t(ji,jj) 
     615             
     616            ! delta at T points 
     617            zdelta         = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * usecc2 )   
     618            rswitch        = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 
     619            delta_i(ji,jj) = zdelta + rn_creepl * rswitch 
     620 
    570621         END DO 
    571622      END DO 
    572  
    573       CALL lbc_lnk_multi( u_ice(:,:), 'U', -1., v_ice(:,:), 'V', -1. ) 
    574  
    575 #if defined key_agrif && defined key_lim2 
    576       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'U' ) 
    577       CALL agrif_rhg_lim2( nn_nevp , nn_nevp, 'V' ) 
    578 #endif 
    579 #if defined key_bdy 
    580       CALL bdy_ice_lim_dyn( 'U' ) 
    581       CALL bdy_ice_lim_dyn( 'V' ) 
    582 #endif          
    583  
    584       DO jj = k_j1+1, k_jpj-1  
    585          DO ji = fs_2, fs_jpim1 
    586             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    587                v_ice1(ji,jj)  = 0.5_wp * ( ( v_ice(ji  ,jj) + v_ice(ji,  jj-1) ) * e1t(ji+1,jj)     & 
    588                   &                      + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji  ,jj) )   & 
    589                   &                      / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 
    590  
    591                u_ice2(ji,jj)  = 0.5_wp * ( ( u_ice(ji,jj  ) + u_ice(ji-1,jj  ) ) * e2t(ji,jj+1)     & 
    592                   &                      + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj  ) )   & 
    593                   &                      / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 
    594             ENDIF  
    595          END DO 
    596       END DO 
    597  
    598       CALL lbc_lnk_multi( u_ice2(:,:), 'V', -1., v_ice1(:,:), 'U', -1. ) 
    599  
    600       ! Recompute delta, shear and div, inputs for mechanical redistribution  
    601       DO jj = k_j1+1, k_jpj-1 
    602          DO ji = fs_2, jpim1   !RB bug no vect opt due to zmask 
    603             !- divu_i(:,:), zdt(:,:): divergence and tension at centre  
    604             !- zds(:,:): shear on northeast corner of grid cells 
    605             IF ( vt_i(ji,jj) <= zvmin ) THEN 
    606  
    607                divu_i(ji,jj) = (  e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj  ) * u_ice(ji-1,jj  )   & 
    608                   &             + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji  ,jj-1) * v_ice(ji  ,jj-1)   & 
    609                   &            ) * r1_e12t(ji,jj) 
    610  
    611                zdt(ji,jj) = ( ( u_ice(ji,jj) * r1_e2u(ji,jj) - u_ice(ji-1,jj) * r1_e2u(ji-1,jj) ) * e2t(ji,jj) * e2t(ji,jj)  & 
    612                   &          -( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj)  & 
    613                   &         ) * r1_e12t(ji,jj) 
    614                ! 
    615                ! SB modif because ocean has no slip boundary condition  
    616                zds(ji,jj) = ( ( u_ice(ji,jj+1) * r1_e1u(ji,jj+1) - u_ice(ji,jj) * r1_e1u(ji,jj) ) * e1f(ji,jj) * e1f(ji,jj)  & 
    617                   &          +( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj)  & 
    618                   &         ) * r1_e12f(ji,jj) * ( 2.0 - fmask(ji,jj,1) )                                     & 
    619                   &         * zmask(ji,jj) * zmask(ji,jj+1) * zmask(ji+1,jj) * zmask(ji+1,jj+1) 
    620  
    621                zdst = ( e2u(ji,jj) * v_ice1(ji,jj) - e2u(ji-1,jj  ) * v_ice1(ji-1,jj  )    & 
    622                   &   + e1v(ji,jj) * u_ice2(ji,jj) - e1v(ji  ,jj-1) * u_ice2(ji  ,jj-1) ) * r1_e12t(ji,jj) 
    623  
    624                delta = SQRT( divu_i(ji,jj)**2 + ( zdt(ji,jj)**2 + zdst**2 ) * usecc2 )   
    625                delta_i(ji,jj) = delta + rn_creepl 
    626              
    627             ENDIF 
    628          END DO 
    629       END DO 
    630       ! 
    631       !------------------------------------------------------------------------------! 
    632       ! 5) Store stress tensor and its invariants 
    633       !------------------------------------------------------------------------------! 
    634       ! * Invariants of the stress tensor are required for limitd_me 
    635       !   (accelerates convergence and improves stability) 
    636       DO jj = k_j1+1, k_jpj-1 
    637          DO ji = fs_2, fs_jpim1 
    638             zdst           = (  e2u(ji,jj) * v_ice1(ji,jj) - e2u( ji-1, jj   ) * v_ice1(ji-1,jj)  &    
    639                &              + e1v(ji,jj) * u_ice2(ji,jj) - e1v( ji  , jj-1 ) * u_ice2(ji,jj-1) ) * r1_e12t(ji,jj)  
    640             shear_i(ji,jj) = SQRT( zdt(ji,jj) * zdt(ji,jj) + zdst * zdst ) 
    641          END DO 
    642       END DO 
    643  
    644       ! Lateral boundary condition 
    645       CALL lbc_lnk_multi( divu_i (:,:), 'T', 1., delta_i(:,:), 'T', 1.,  shear_i(:,:), 'T', 1. ) 
    646  
    647       ! * Store the stress tensor for the next time step 
     623      CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 
     624       
     625      ! --- Store the stress tensor for the next time step --- ! 
    648626      stress1_i (:,:) = zs1 (:,:) 
    649627      stress2_i (:,:) = zs2 (:,:) 
     
    652630      ! 
    653631      !------------------------------------------------------------------------------! 
    654       ! 6) Control prints of residual and charge ellipse 
     632      ! 5) Control prints of residual and charge ellipse 
    655633      !------------------------------------------------------------------------------! 
    656634      ! 
     
    675653               DO ji = 2, jpim1 
    676654                  IF (zpresh(ji,jj) > 1.0) THEN 
    677                      sigma1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
    678                      sigma2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
     655                     zsig1 = ( zs1(ji,jj) + (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) )  
     656                     zsig2 = ( zs1(ji,jj) - (zs2(ji,jj)**2 + 4*zs12(ji,jj)**2 )**0.5 ) / ( 2*zpresh(ji,jj) ) 
    679657                     WRITE(charout,FMT="('lim_rhg  :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 
    680658                     CALL prt_ctl_info(charout) 
     
    687665      ENDIF 
    688666      ! 
    689       CALL wrk_dealloc( jpi,jpj, zpresh, zfrld1, zmass1, zcorl1, za1ct , zpreshc, zfrld2, zmass2, zcorl2, za2ct ) 
    690       CALL wrk_dealloc( jpi,jpj, u_oce2, u_ice2, v_oce1 , v_ice1 , zmask               ) 
    691       CALL wrk_dealloc( jpi,jpj, zf1   , zu_ice, zf2   , zv_ice , zdt    , zds  ) 
    692       CALL wrk_dealloc( jpi,jpj, zs1   , zs2   , zs12   , zresr , zpice                 ) 
     667      CALL wrk_dealloc( jpi,jpj, zpresh, z1_e1t0, z1_e2t0, zp_delt ) 
     668      CALL wrk_dealloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 
     669      CALL wrk_dealloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 
     670      CALL wrk_dealloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 
     671      CALL wrk_dealloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 
    693672 
    694673   END SUBROUTINE lim_rhg 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90

    r6399 r7067  
    110110      !!--------------------------------------------------------------------- 
    111111 
    112       ! make calls for heat fluxes before it is modified 
    113       ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
    114       IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
    115       IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
    116       IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
    117       IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
    118       IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
    119       IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )   
    120       IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
    121          &                                                      * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) 
    122       IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
    123       IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
    124       IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   ! emp over ocean (taking into account the snow blown away from the ice) 
    125       IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   ! emp over ice   (taking into account the snow blown away from the ice) 
    126  
    127       ! albedo output 
     112      ! make call for albedo output before it is modified 
    128113      CALL wrk_alloc( jpi,jpj, zalb )     
    129114 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90

    r6469 r7067  
    5454   PUBLIC   lim_var_eqv2glo       
    5555   PUBLIC   lim_var_salprof       
    56    PUBLIC   lim_var_icetm         
    5756   PUBLIC   lim_var_bv            
    5857   PUBLIC   lim_var_salprof1d     
     
    8988      ! Compute variables 
    9089      !-------------------- 
    91       vt_i (:,:) = 0._wp 
    92       vt_s (:,:) = 0._wp 
    93       at_i (:,:) = 0._wp 
    94       ato_i(:,:) = 1._wp 
    95       ! 
    96       DO jl = 1, jpl 
    97          DO jj = 1, jpj 
    98             DO ji = 1, jpi 
    99                ! 
    100                vt_i(ji,jj) = vt_i(ji,jj) + v_i(ji,jj,jl) ! ice volume 
    101                vt_s(ji,jj) = vt_s(ji,jj) + v_s(ji,jj,jl) ! snow volume 
    102                at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl) ! ice concentration 
    103                ! 
    104                rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )  
    105                icethi(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch  ! ice thickness 
    106             END DO 
    107          END DO 
    108       END DO 
    109  
     90      ! integrated values 
     91      vt_i (:,:) = SUM( v_i, dim=3 ) 
     92      vt_s (:,:) = SUM( v_s, dim=3 ) 
     93      at_i (:,:) = SUM( a_i, dim=3 ) 
     94      et_s(:,:)  = SUM( SUM( e_s(:,:,:,:), dim=4 ), dim=3 ) 
     95      et_i(:,:)  = SUM( SUM( e_i(:,:,:,:), dim=4 ), dim=3 ) 
     96      ! 
    11097      DO jj = 1, jpj 
    11198         DO ji = 1, jpi 
     
    115102 
    116103      IF( kn > 1 ) THEN 
    117          et_s (:,:) = 0._wp 
    118          ot_i (:,:) = 0._wp 
     104         ! 
     105         ! mean ice/snow thickness 
     106         DO jj = 1, jpj 
     107            DO ji = 1, jpi 
     108               rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )  
     109               htm_i(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 
     110               htm_s(ji,jj) = vt_s(ji,jj) / MAX( at_i(ji,jj) , epsi10 ) * rswitch 
     111            ENDDO 
     112         ENDDO 
     113 
     114         ! mean temperature (K), salinity and age 
    119115         smt_i(:,:) = 0._wp 
    120          et_i (:,:) = 0._wp 
    121          ! 
     116         tm_i(:,:)  = 0._wp 
     117         tm_su(:,:) = 0._wp 
     118         om_i (:,:) = 0._wp 
    122119         DO jl = 1, jpl 
     120             
    123121            DO jj = 1, jpj 
    124122               DO ji = 1, jpi 
    125                   et_s(ji,jj)  = et_s(ji,jj)  + e_s(ji,jj,1,jl)                                           ! snow heat content 
    126                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi20 ) )  
    127                   smt_i(ji,jj) = smt_i(ji,jj) + smv_i(ji,jj,jl) / MAX( vt_i(ji,jj) , epsi20 ) * rswitch   ! ice salinity 
    128                   rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi20 ) )  
    129                   ot_i(ji,jj)  = ot_i(ji,jj)  + oa_i(ji,jj,jl)  / MAX( at_i(ji,jj) , epsi20 ) * rswitch   ! ice age 
    130                END DO 
    131             END DO 
    132          END DO 
    133          ! 
    134          DO jl = 1, jpl 
     123                  rswitch      = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 
     124                  tm_su(ji,jj) = tm_su(ji,jj) + rswitch * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi10 ) 
     125                  om_i (ji,jj) = om_i (ji,jj) + rswitch *   oa_i(ji,jj,jl)                         / MAX( at_i(ji,jj) , epsi10 ) 
     126               END DO 
     127            END DO 
     128             
    135129            DO jk = 1, nlay_i 
    136                et_i(:,:) = et_i(:,:) + e_i(:,:,jk,jl)       ! ice heat content 
    137             END DO 
    138          END DO 
     130               DO jj = 1, jpj 
     131                  DO ji = 1, jpi 
     132                     rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
     133                     tm_i(ji,jj)  = tm_i(ji,jj)  + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
     134                        &            / MAX( vt_i(ji,jj) , epsi10 ) 
     135                     smt_i(ji,jj) = smt_i(ji,jj) + r1_nlay_i * rswitch * s_i(ji,jj,jk,jl) * v_i(ji,jj,jl)  & 
     136                        &            / MAX( vt_i(ji,jj) , epsi10 ) 
     137                  END DO 
     138               END DO 
     139            END DO 
     140         END DO 
     141         tm_i  = tm_i  + rt0 
     142         tm_su = tm_su + rt0 
    139143         ! 
    140144      ENDIF 
     
    246250      ! Mean temperature 
    247251      !------------------- 
    248       vt_i (:,:) = 0._wp 
    249       DO jl = 1, jpl 
    250          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    251       END DO 
     252      ! integrated values 
     253      vt_i (:,:) = SUM( v_i, dim=3 ) 
     254      vt_s (:,:) = SUM( v_s, dim=3 ) 
     255      at_i (:,:) = SUM( a_i, dim=3 ) 
    252256 
    253257      tm_i(:,:) = 0._wp 
     
    397401   END SUBROUTINE lim_var_salprof 
    398402 
    399  
    400    SUBROUTINE lim_var_icetm 
    401       !!------------------------------------------------------------------ 
    402       !!                ***  ROUTINE lim_var_icetm *** 
    403       !! 
    404       !! ** Purpose :   computes mean sea ice temperature 
     403   SUBROUTINE lim_var_bv 
     404      !!------------------------------------------------------------------ 
     405      !!                ***  ROUTINE lim_var_bv *** 
     406      !! 
     407      !! ** Purpose :   computes mean brine volume (%) in sea ice 
     408      !! 
     409      !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
     410      !! 
     411      !! References : Vancoppenolle et al., JGR, 2007 
    405412      !!------------------------------------------------------------------ 
    406413      INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    407414      !!------------------------------------------------------------------ 
    408  
    409       ! Mean sea ice temperature 
    410       vt_i (:,:) = 0._wp 
    411       DO jl = 1, jpl 
    412          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    413       END DO 
    414  
    415       tm_i(:,:) = 0._wp 
     415      ! 
     416      bvm_i(:,:)   = 0._wp 
     417      bv_i (:,:,:) = 0._wp 
    416418      DO jl = 1, jpl 
    417419         DO jk = 1, nlay_i 
    418420            DO jj = 1, jpj 
    419421               DO ji = 1, jpi 
    420                   rswitch = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
    421                   tm_i(ji,jj) = tm_i(ji,jj) + r1_nlay_i * rswitch * ( t_i(ji,jj,jk,jl) - rt0 ) * v_i(ji,jj,jl)  & 
    422                      &            / MAX( vt_i(ji,jj) , epsi10 ) 
    423                END DO 
    424             END DO 
    425          END DO 
    426       END DO 
    427       tm_i = tm_i + rt0 
    428  
    429    END SUBROUTINE lim_var_icetm 
    430  
    431  
    432    SUBROUTINE lim_var_bv 
    433       !!------------------------------------------------------------------ 
    434       !!                ***  ROUTINE lim_var_bv *** 
    435       !! 
    436       !! ** Purpose :   computes mean brine volume (%) in sea ice 
    437       !! 
    438       !! ** Method  : e = - 0.054 * S (ppt) / T (C) 
    439       !! 
    440       !! References : Vancoppenolle et al., JGR, 2007 
    441       !!------------------------------------------------------------------ 
    442       INTEGER  ::   ji, jj, jk, jl   ! dummy loop indices 
    443       REAL(wp) ::   zbvi             ! local scalars 
    444       !!------------------------------------------------------------------ 
    445       ! 
    446       vt_i (:,:) = 0._wp 
    447       DO jl = 1, jpl 
    448          vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl) 
    449       END DO 
    450  
    451       bv_i(:,:) = 0._wp 
    452       DO jl = 1, jpl 
    453          DO jk = 1, nlay_i 
    454             DO jj = 1, jpj 
    455                DO ji = 1, jpi 
    456                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  ) 
    457                   zbvi  = - rswitch * tmut * s_i(ji,jj,jk,jl) / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 )   & 
    458                      &                   * v_i(ji,jj,jl) * r1_nlay_i 
    459                   rswitch = (  1._wp - MAX( 0._wp , SIGN( 1._wp , - vt_i(ji,jj) + epsi20 ) )  ) 
    460                   bv_i(ji,jj) = bv_i(ji,jj) + rswitch * zbvi  / MAX( vt_i(ji,jj) , epsi20 ) 
    461                END DO 
     422                  rswitch        = ( 1._wp - MAX( 0._wp , SIGN( 1._wp , (t_i(ji,jj,jk,jl) - rt0) + epsi10 ) )  ) 
     423                  bv_i(ji,jj,jl) = bv_i(ji,jj,jl) - rswitch * tmut * s_i(ji,jj,jk,jl) * r1_nlay_i  & 
     424                     &                            / MIN( t_i(ji,jj,jk,jl) - rt0, - epsi10 ) 
     425               END DO 
     426            END DO 
     427         END DO 
     428          
     429         DO jj = 1, jpj 
     430            DO ji = 1, jpi 
     431               rswitch      = MAX( 0._wp , SIGN( 1._wp , vt_i(ji,jj) - epsi10 ) ) 
     432               bvm_i(ji,jj) = bvm_i(ji,jj) + rswitch * bv_i(ji,jj,jl) * v_i(ji,jj,jl) / MAX( vt_i(ji,jj), epsi10 ) 
    462433            END DO 
    463434         END DO 
     
    715686            zht_i(ji,1:jpl) = 0._wp 
    716687            za_i (ji,1:jpl) = 0._wp 
    717              
     688            itest(:)        = 0       
     689       
    718690            ! *** case very thin ice: fill only category 1 
    719691            IF ( i_fill == 1 ) THEN 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90

    r6417 r7067  
    1717   USE sbc_oce         ! Surface boundary condition: ocean fields 
    1818   USE sbc_ice         ! Surface boundary condition: ice fields 
    19    USE dom_ice 
    2019   USE ice 
    2120   USE limvar 
     
    4039   !!---------------------------------------------------------------------- 
    4140CONTAINS 
    42  
    43 #if defined key_dimgout 
    44 # include "limwri_dimg.h90" 
    45 #else 
    4641 
    4742   SUBROUTINE lim_wri( kindic ) 
     
    5954      INTEGER  ::  ji, jj, jk, jl  ! dummy loop indices 
    6055      REAL(wp) ::  z1_365 
    61       REAL(wp) ::  ztmp 
    62       REAL(wp), POINTER, DIMENSION(:,:,:) ::  zoi, zei, zt_i, zt_s 
    63       REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, z2da, z2db, zswi    ! 2D workspace 
     56      REAL(wp) ::  z2da, z2db, ztmp 
     57      REAL(wp), POINTER, DIMENSION(:,:,:) ::  zswi2 
     58      REAL(wp), POINTER, DIMENSION(:,:)   ::  z2d, zswi    ! 2D workspace 
    6459      !!------------------------------------------------------------------- 
    6560 
    6661      IF( nn_timing == 1 )  CALL timing_start('limwri') 
    6762 
    68       CALL wrk_alloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s ) 
    69       CALL wrk_alloc( jpi, jpj     , z2d, z2da, z2db, zswi ) 
     63      CALL wrk_alloc( jpi, jpj, jpl, zswi2 ) 
     64      CALL wrk_alloc( jpi, jpj     , z2d, zswi ) 
    7065 
    7166      !----------------------------- 
     
    7469      z1_365 = 1._wp / 365._wp 
    7570 
    76       CALL lim_var_icetm      ! mean sea ice temperature 
    77  
    78       CALL lim_var_bv         ! brine volume 
    79  
    80       DO jj = 1, jpj          ! presence indicator of ice 
     71      ! brine volume 
     72      CALL lim_var_bv  
     73 
     74      ! tresholds for outputs 
     75      DO jj = 1, jpj 
    8176         DO ji = 1, jpi 
    8277            zswi(ji,jj)  = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) 
    8378         END DO 
    8479      END DO 
    85       ! 
    86       ! 
    87       !                                              
    88       IF ( iom_use( "icethic_cea" ) ) THEN                       ! mean ice thickness 
    89          DO jj = 1, jpj  
     80      DO jl = 1, jpl 
     81         DO jj = 1, jpj 
    9082            DO ji = 1, jpi 
    91                z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 
     83               zswi2(ji,jj,jl)  = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    9284            END DO 
    9385         END DO 
    94          CALL iom_put( "icethic_cea"  , z2d              ) 
    95       ENDIF 
    96  
    97       IF ( iom_use( "snowthic_cea" ) ) THEN                      ! snow thickness = mean snow thickness over the cell  
    98          DO jj = 1, jpj                                             
    99             DO ji = 1, jpi 
    100                z2d(ji,jj)  = vt_s(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) 
    101             END DO 
    102          END DO 
    103          CALL iom_put( "snowthic_cea" , z2d              )        
    104       ENDIF 
     86      END DO 
    10587      ! 
     88      ! fluxes  
     89      ! pfrld is the lead fraction at the previous time step (actually between TRP and THD) 
     90      IF( iom_use('qsr_oce') )   CALL iom_put( "qsr_oce" , qsr_oce(:,:) * pfrld(:,:) )                                   !     solar flux at ocean surface 
     91      IF( iom_use('qns_oce') )   CALL iom_put( "qns_oce" , qns_oce(:,:) * pfrld(:,:) + qemp_oce(:,:) )                   ! non-solar flux at ocean surface 
     92      IF( iom_use('qsr_ice') )   CALL iom_put( "qsr_ice" , SUM( qsr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux at ice surface 
     93      IF( iom_use('qns_ice') )   CALL iom_put( "qns_ice" , SUM( qns_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) + qemp_ice(:,:) ) ! non-solar flux at ice surface 
     94      IF( iom_use('qtr_ice') )   CALL iom_put( "qtr_ice" , SUM( ftr_ice(:,:,:) * a_i_b(:,:,:), dim=3 ) )                 !     solar flux transmitted thru ice 
     95      IF( iom_use('qt_oce' ) )   CALL iom_put( "qt_oce"  , ( qsr_oce(:,:) + qns_oce(:,:) ) * pfrld(:,:) + qemp_oce(:,:) )   
     96      IF( iom_use('qt_ice' ) )   CALL iom_put( "qt_ice"  , SUM( ( qns_ice(:,:,:) + qsr_ice(:,:,:) )   & 
     97         &                                                      * a_i_b(:,:,:),dim=3 ) + qemp_ice(:,:) ) 
     98      IF( iom_use('qemp_oce') )  CALL iom_put( "qemp_oce" , qemp_oce(:,:) )   
     99      IF( iom_use('qemp_ice') )  CALL iom_put( "qemp_ice" , qemp_ice(:,:) )   
     100      IF( iom_use('emp_oce' ) )  CALL iom_put( "emp_oce"  , emp_oce(:,:) )   !emp over ocean (taking into account the snow blown away from the ice) 
     101      IF( iom_use('emp_ice' ) )  CALL iom_put( "emp_ice"  , emp_ice(:,:) )   !emp over ice   (taking into account the snow blown away from the ice) 
     102 
     103      ! velocity 
    106104      IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN  
    107105         DO jj = 2 , jpjm1 
    108106            DO ji = 2 , jpim1 
    109                z2da(ji,jj)  = (  u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 
    110                z2db(ji,jj)  = (  v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 
     107               z2da  = ( u_ice(ji,jj) * umask(ji,jj,1) + u_ice(ji-1,jj) * umask(ji-1,jj,1) ) * 0.5_wp 
     108               z2db  = ( v_ice(ji,jj) * vmask(ji,jj,1) + v_ice(ji,jj-1) * vmask(ji,jj-1,1) ) * 0.5_wp 
     109               z2d(ji,jj) = SQRT( z2da * z2da + z2db * z2db ) 
    111110           END DO 
    112111         END DO 
    113          CALL lbc_lnk( z2da, 'T', -1. ) 
    114          CALL lbc_lnk( z2db, 'T', -1. ) 
    115          CALL iom_put( "uice_ipa"     , z2da             )       ! ice velocity u component 
    116          CALL iom_put( "vice_ipa"     , z2db             )       ! ice velocity v component 
    117          DO jj = 1, jpj                                  
    118             DO ji = 1, jpi 
    119                z2d(ji,jj)  = SQRT( z2da(ji,jj) * z2da(ji,jj) + z2db(ji,jj) * z2db(ji,jj) )  
    120             END DO 
    121          END DO 
    122          CALL iom_put( "icevel"       , z2d              )       ! ice velocity module 
     112         CALL lbc_lnk( z2d, 'T', 1. ) 
     113         CALL iom_put( "uice_ipa"     , u_ice      )       ! ice velocity u component 
     114         CALL iom_put( "vice_ipa"     , v_ice      )       ! ice velocity v component 
     115         CALL iom_put( "icevel"       , z2d        )       ! ice velocity module 
    123116      ENDIF 
    124117      ! 
    125       IF ( iom_use( "miceage" ) ) THEN  
    126          z2d(:,:) = 0.e0 
    127          DO jl = 1, jpl 
    128             DO jj = 1, jpj 
    129                DO ji = 1, jpi 
    130                   rswitch    = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
    131                   z2d(ji,jj) = z2d(ji,jj) + rswitch * oa_i(ji,jj,jl) / MAX( at_i(ji,jj), 0.1 ) 
    132                END DO 
    133             END DO 
    134          END DO 
    135          CALL iom_put( "miceage"     , z2d * z1_365      )        ! mean ice age 
    136       ENDIF 
    137  
    138       IF ( iom_use( "micet" ) ) THEN  
    139          DO jj = 1, jpj 
    140             DO ji = 1, jpi 
    141                z2d(ji,jj) = ( tm_i(ji,jj) - rt0 ) * zswi(ji,jj) 
    142             END DO 
    143          END DO 
    144          CALL iom_put( "micet"       , z2d               )        ! mean ice temperature 
    145       ENDIF 
     118      IF ( iom_use( "miceage" ) )       CALL iom_put( "miceage"     , om_i * zswi * z1_365   )  ! mean ice age 
     119      IF ( iom_use( "icethic_cea" ) )   CALL iom_put( "icethic_cea" , htm_i * zswi           )  ! ice thickness mean 
     120      IF ( iom_use( "snowthic_cea" ) )  CALL iom_put( "snowthic_cea", htm_s * zswi           )  ! snow thickness mean 
     121      IF ( iom_use( "micet" ) )         CALL iom_put( "micet"       , ( tm_i  - rt0 ) * zswi )  ! ice mean    temperature 
     122      IF ( iom_use( "icest" ) )         CALL iom_put( "icest"       , ( tm_su - rt0 ) * zswi )  ! ice surface temperature 
     123      IF ( iom_use( "icecolf" ) )       CALL iom_put( "icecolf"     , hicol                  )  ! frazil ice collection thickness 
    146124      ! 
    147       IF ( iom_use( "icest" ) ) THEN  
    148          z2d(:,:) = 0.e0 
    149          DO jl = 1, jpl 
    150             DO jj = 1, jpj 
    151                DO ji = 1, jpi 
    152                   z2d(ji,jj) = z2d(ji,jj) + zswi(ji,jj) * ( t_su(ji,jj,jl) - rt0 ) * a_i(ji,jj,jl) / MAX( at_i(ji,jj) , epsi06 ) 
    153                END DO 
    154             END DO 
    155          END DO 
    156          CALL iom_put( "icest"       , z2d              )        ! ice surface temperature 
    157       ENDIF 
    158  
    159       IF ( iom_use( "icecolf" ) )   CALL iom_put( "icecolf", hicol )  ! frazil ice collection thickness 
    160  
    161125      CALL iom_put( "isst"        , sst_m               )        ! sea surface temperature 
    162126      CALL iom_put( "isss"        , sss_m               )        ! sea surface salinity 
    163       CALL iom_put( "iceconc"     , at_i                )        ! ice concentration 
    164       CALL iom_put( "icevolu"     , vt_i                )        ! ice volume = mean ice thickness over the cell 
    165       CALL iom_put( "icehc"       , et_i                )        ! ice total heat content 
    166       CALL iom_put( "isnowhc"     , et_s                )        ! snow total heat content 
    167       CALL iom_put( "ibrinv"      , bv_i * 100._wp      )        ! brine volume 
     127      CALL iom_put( "iceconc"     , at_i  * zswi        )        ! ice concentration 
     128      CALL iom_put( "icevolu"     , vt_i  * zswi        )        ! ice volume = mean ice thickness over the cell 
     129      CALL iom_put( "icehc"       , et_i  * zswi        )        ! ice total heat content 
     130      CALL iom_put( "isnowhc"     , et_s  * zswi        )        ! snow total heat content 
     131      CALL iom_put( "ibrinv"      , bvm_i * zswi * 100. )        ! brine volume 
    168132      CALL iom_put( "utau_ice"    , utau_ice            )        ! wind stress over ice along i-axis at I-point 
    169133      CALL iom_put( "vtau_ice"    , vtau_ice            )        ! wind stress over ice along j-axis at I-point 
    170134      CALL iom_put( "snowpre"     , sprecip * 86400.    )        ! snow precipitation  
    171       CALL iom_put( "micesalt"    , smt_i               )        ! mean ice salinity 
    172  
    173       CALL iom_put( "icestr"      , strength * 0.001    )        ! ice strength 
    174       CALL iom_put( "idive"       , divu_i * 1.0e8      )        ! divergence 
    175       CALL iom_put( "ishear"      , shear_i * 1.0e8     )        ! shear 
    176       CALL iom_put( "snowvol"     , vt_s                )        ! snow volume 
     135      CALL iom_put( "micesalt"    , smt_i   * zswi      )        ! mean ice salinity 
     136 
     137      CALL iom_put( "icestr"      , strength * zswi )    ! ice strength 
     138      CALL iom_put( "idive"       , divu_i * 1.0e8      )    ! divergence 
     139      CALL iom_put( "ishear"      , shear_i * 1.0e8     )    ! shear 
     140      CALL iom_put( "snowvol"     , vt_s   * zswi       )        ! snow volume 
    177141       
    178142      CALL iom_put( "icetrp"      , diag_trp_vi * rday  )        ! ice volume transport 
     
    183147 
    184148      CALL iom_put( "sfxbog"      , sfx_bog * rday      )        ! salt flux from bottom growth 
    185       CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melt 
    186       CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melt 
     149      CALL iom_put( "sfxbom"      , sfx_bom * rday      )        ! salt flux from bottom melting 
     150      CALL iom_put( "sfxsum"      , sfx_sum * rday      )        ! salt flux from surface melting 
    187151      CALL iom_put( "sfxsni"      , sfx_sni * rday      )        ! salt flux from snow ice formation 
    188152      CALL iom_put( "sfxopw"      , sfx_opw * rday      )        ! salt flux from open water formation 
    189153      CALL iom_put( "sfxdyn"      , sfx_dyn * rday      )        ! salt flux from ridging rafting 
    190       CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from residual 
     154      CALL iom_put( "sfxres"      , sfx_res * rday      )        ! salt flux from limupdate (resultant) 
    191155      CALL iom_put( "sfxbri"      , sfx_bri * rday      )        ! salt flux from brines 
    192156      CALL iom_put( "sfxsub"      , sfx_sub * rday      )        ! salt flux from sublimation 
     
    202166      CALL iom_put( "vfxbom"     , wfx_bom * ztmp       )        ! bottom melt  
    203167      CALL iom_put( "vfxice"     , wfx_ice * ztmp       )        ! total ice growth/melt  
     168 
     169      IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
     170         WHERE( htm_i(:,:) < 0.2 .AND. htm_i(:,:) > 0. ) ; z2d = wfx_bog 
     171         ELSEWHERE                                       ; z2d = 0._wp 
     172         END WHERE 
     173         CALL iom_put( "vfxthin", ( wfx_opw + z2d ) * ztmp ) 
     174      ENDIF 
     175 
     176      ztmp = rday / rhosn 
     177      CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    204178      CALL iom_put( "vfxsnw"     , wfx_snw * ztmp       )        ! total snw growth/melt  
    205       CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow)  
    206       CALL iom_put( "vfxspr"     , wfx_spr * ztmp       )        ! precip (snow) 
    207        
     179      CALL iom_put( "vfxsub"     , wfx_sub * ztmp       )        ! sublimation (snow/ice)  
     180      CALL iom_put( "vfxsub_err" , wfx_err_sub * ztmp   )        ! "excess" of sublimation sent to ocean       
     181  
    208182      CALL iom_put( "afxtot"     , afx_tot * rday       )        ! concentration tendency (total) 
    209183      CALL iom_put( "afxdyn"     , afx_dyn * rday       )        ! concentration tendency (dynamics) 
     
    225199      CALL iom_put ('hfxdif'     , hfx_dif(:,:)         )   !   
    226200      CALL iom_put ('hfxopw'     , hfx_opw(:,:)         )   !   
    227       CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM(a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base  
     201      CALL iom_put ('hfxtur'     , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3 ) ) ! turbulent heat flux at ice base  
    228202      CALL iom_put ('hfxdhc'     , diag_heat(:,:)       )   ! Heat content variation in snow and ice  
    229203      CALL iom_put ('hfxspr'     , hfx_spr(:,:)         )   ! Heat content of snow precip  
    230204 
    231  
    232       IF ( iom_use( "vfxthin" ) ) THEN   ! ice production for open water + thin ice (<20cm) => comparable to observations   
    233          DO jj = 1, jpj  
    234             DO ji = 1, jpi 
    235                z2d(ji,jj)  = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness 
    236             END DO 
    237          END DO 
    238          WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog 
    239          ELSEWHERE                                   ; z2da = 0._wp 
    240          END WHERE 
    241          CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp ) 
    242       ENDIF 
    243205       
    244206      !-------------------------------- 
    245207      ! Output values for each category 
    246208      !-------------------------------- 
    247       CALL iom_put( "iceconc_cat"      , a_i         )        ! area for categories 
    248       CALL iom_put( "icethic_cat"      , ht_i        )        ! thickness for categories 
    249       CALL iom_put( "snowthic_cat"     , ht_s        )        ! snow depth for categories 
    250       CALL iom_put( "salinity_cat"     , sm_i        )        ! salinity for categories 
    251  
     209      IF ( iom_use( "iceconc_cat"  ) )  CALL iom_put( "iceconc_cat"      , a_i   * zswi2   )        ! area for categories 
     210      IF ( iom_use( "icethic_cat"  ) )  CALL iom_put( "icethic_cat"      , ht_i  * zswi2   )        ! thickness for categories 
     211      IF ( iom_use( "snowthic_cat" ) )  CALL iom_put( "snowthic_cat"     , ht_s  * zswi2   )        ! snow depth for categories 
     212      IF ( iom_use( "salinity_cat" ) )  CALL iom_put( "salinity_cat"     , sm_i  * zswi2   )        ! salinity for categories 
    252213      ! ice temperature 
    253       IF ( iom_use( "icetemp_cat" ) ) THEN  
    254          zt_i(:,:,:) = SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i 
    255          CALL iom_put( "icetemp_cat"   , zt_i - rt0  ) 
    256       ENDIF 
    257        
     214      IF ( iom_use( "icetemp_cat"  ) )  CALL iom_put( "icetemp_cat", ( SUM( t_i(:,:,:,:), dim=3 ) * r1_nlay_i - rt0 ) * zswi2 ) 
    258215      ! snow temperature 
    259       IF ( iom_use( "snwtemp_cat" ) ) THEN  
    260          zt_s(:,:,:) = SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s 
    261          CALL iom_put( "snwtemp_cat"   , zt_s - rt0  ) 
    262       ENDIF 
    263  
    264       ! Compute ice age 
    265       IF ( iom_use( "iceage_cat" ) ) THEN  
    266          DO jl = 1, jpl  
    267             DO jj = 1, jpj 
    268                DO ji = 1, jpi 
    269                   rswitch = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - 0.1 ) ) 
    270                   rswitch = rswitch * MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - 0.1 ) ) 
    271                   zoi(ji,jj,jl) = oa_i(ji,jj,jl)  / MAX( a_i(ji,jj,jl) , 0.1 ) * rswitch 
    272                END DO 
    273             END DO 
    274          END DO 
    275          CALL iom_put( "iceage_cat"   , zoi * z1_365 )        ! ice age for categories 
    276       ENDIF 
    277  
    278       ! Compute brine volume 
    279       IF ( iom_use( "brinevol_cat" ) ) THEN  
    280          zei(:,:,:) = 0._wp 
    281          DO jl = 1, jpl  
    282             DO jk = 1, nlay_i 
    283                DO jj = 1, jpj 
    284                   DO ji = 1, jpi 
    285                      rswitch = MAX( 0._wp , SIGN( 1._wp , a_i(ji,jj,jl) - epsi06 ) ) 
    286                      zei(ji,jj,jl) = zei(ji,jj,jl) + 100.0 *  & 
    287                         ( - tmut * s_i(ji,jj,jk,jl) / MIN( ( t_i(ji,jj,jk,jl) - rt0 ), - epsi06 ) ) * & 
    288                         rswitch * r1_nlay_i 
    289                   END DO 
    290                END DO 
    291             END DO 
    292          END DO 
    293          CALL iom_put( "brinevol_cat"     , zei      )        ! brine volume for categories 
    294       ENDIF 
     216      IF ( iom_use( "snwtemp_cat"  ) )  CALL iom_put( "snwtemp_cat", ( SUM( t_s(:,:,:,:), dim=3 ) * r1_nlay_s - rt0 ) * zswi2 ) 
     217      ! ice age 
     218      IF ( iom_use( "iceage_cat"   ) )  CALL iom_put( "iceage_cat" , o_i * zswi2 * z1_365 ) 
     219      ! brine volume 
     220      IF ( iom_use( "brinevol_cat" ) )  CALL iom_put( "brinevol_cat", bv_i * 100. * zswi2 ) 
    295221 
    296222      !     !  Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s 
     
    298224      !     not yet implemented 
    299225       
    300       CALL wrk_dealloc( jpi, jpj, jpl, zoi, zei, zt_i, zt_s ) 
    301       CALL wrk_dealloc( jpi, jpj     , z2d, zswi, z2da, z2db ) 
     226      CALL wrk_dealloc( jpi, jpj, jpl, zswi2 ) 
     227      CALL wrk_dealloc( jpi, jpj     , z2d, zswi ) 
    302228 
    303229      IF( nn_timing == 1 )  CALL timing_stop('limwri') 
    304230       
    305231   END SUBROUTINE lim_wri 
    306 #endif 
    307232 
    308233  
     
    319244      !!   4.0  !  2013-06  (C. Rousset) 
    320245      !!---------------------------------------------------------------------- 
    321       INTEGER, INTENT( in ) ::   kt               ! ocean time-step index) 
    322       INTEGER, INTENT( in ) ::   kid , kh_i        
     246      INTEGER, INTENT( in )   ::   kt               ! ocean time-step index) 
     247      INTEGER, INTENT( in )   ::   kid , kh_i 
     248      INTEGER                 ::   nz_i, jl 
     249      REAL(wp), DIMENSION(jpl) :: jcat 
    323250      !!---------------------------------------------------------------------- 
    324  
    325       CALL histdef( kid, "iicethic", "Ice thickness"           , "m"      ,   & 
    326       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    327       CALL histdef( kid, "iiceconc", "Ice concentration"       , "%"      ,   & 
    328       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    329       CALL histdef( kid, "iicetemp", "Ice temperature"         , "C"      ,   & 
    330       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    331       CALL histdef( kid, "iicevelu", "i-Ice speed (I-point)"   , "m/s"    ,   & 
    332       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    333       CALL histdef( kid, "iicevelv", "j-Ice speed (I-point)"   , "m/s"    ,   & 
    334       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    335       CALL histdef( kid, "iicestru", "i-Wind stress over ice (I-pt)", "Pa",   & 
    336       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    337       CALL histdef( kid, "iicestrv", "j-Wind stress over ice (I-pt)", "Pa",   & 
    338       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    339       CALL histdef( kid, "iicesflx", "Solar flux over ocean"     , "w/m2" ,   & 
    340       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    341       CALL histdef( kid, "iicenflx", "Non-solar flux over ocean" , "w/m2" ,   & 
     251      DO jl = 1, jpl 
     252         jcat(jl) = REAL(jl) 
     253      ENDDO 
     254       
     255      CALL histvert( kid, "ncatice", "Ice Categories","", jpl, jcat, nz_i, "up") 
     256 
     257      CALL histdef( kid, "sithic", "Ice thickness"           , "m"      ,   & 
     258      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     259      CALL histdef( kid, "siconc", "Ice concentration"       , "%"      ,   & 
     260      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     261      CALL histdef( kid, "sitemp", "Ice temperature"         , "C"      ,   & 
     262      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     263      CALL histdef( kid, "sivelu", "i-Ice speed "            , "m/s"    ,   & 
     264      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     265      CALL histdef( kid, "sivelv", "j-Ice speed "            , "m/s"    ,   & 
     266      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     267      CALL histdef( kid, "sistru", "i-Wind stress over ice " , "Pa"     ,   & 
     268      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     269      CALL histdef( kid, "sistrv", "j-Wind stress over ice " , "Pa"     ,   & 
     270      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     271      CALL histdef( kid, "sisflx", "Solar flux over ocean"     , "w/m2" ,   & 
     272      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     273      CALL histdef( kid, "sinflx", "Non-solar flux over ocean" , "w/m2" ,   & 
    342274      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    343275      CALL histdef( kid, "isnowpre", "Snow precipitation"      , "kg/m2/s",   & 
    344276      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    345       CALL histdef( kid, "iicesali", "Ice salinity"            , "PSU"    ,   & 
    346       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    347       CALL histdef( kid, "iicevolu", "Ice volume"              , "m"      ,   & 
    348       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    349       CALL histdef( kid, "iicedive", "Ice divergence"          , "10-8s-1",   & 
    350       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
    351       CALL histdef( kid, "iicebopr", "Ice bottom production"   , "m/s"    ,   & 
    352       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    353       CALL histdef( kid, "iicedypr", "Ice dynamic production"  , "m/s"    ,   & 
    354       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    355       CALL histdef( kid, "iicelapr", "Ice open water prod"     , "m/s"    ,   & 
     277      CALL histdef( kid, "sisali", "Ice salinity"            , "PSU"    ,   & 
     278      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     279      CALL histdef( kid, "sivolu", "Ice volume"              , "m"      ,   & 
     280      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     281      CALL histdef( kid, "sidive", "Ice divergence"          , "10-8s-1",   & 
     282      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt )  
     283 
     284      CALL histdef( kid, "vfxbog", "Ice bottom production"   , "m/s"    ,   & 
     285      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     286      CALL histdef( kid, "vfxdyn", "Ice dynamic production"  , "m/s"    ,   & 
     287      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     288      CALL histdef( kid, "vfxopw", "Ice open water prod"     , "m/s"    ,   & 
    356289      &       jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    357       CALL histdef( kid, "iicesipr", "Snow ice production "    , "m/s"    ,   & 
    358       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    359       CALL histdef( kid, "iicerepr", "Ice prod from limupdate" , "m/s"    ,   & 
    360       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    361       CALL histdef( kid, "iicebome", "Ice bottom melt"         , "m/s"    ,   & 
    362       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    363       CALL histdef( kid, "iicesume", "Ice surface melt"        , "m/s"    ,   & 
    364       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    365       CALL histdef( kid, "iisfxdyn", "Salt flux from dynmics"  , ""       ,   & 
    366       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
    367       CALL histdef( kid, "iisfxres", "Salt flux from limupdate", ""       ,   & 
    368       &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     290      CALL histdef( kid, "vfxsni", "Snow ice production "    , "m/s"    ,   & 
     291      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     292      CALL histdef( kid, "vfxres", "Ice prod from limupdate" , "m/s"    ,   & 
     293      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     294      CALL histdef( kid, "vfxbom", "Ice bottom melt"         , "m/s"    ,   & 
     295      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     296      CALL histdef( kid, "vfxsum", "Ice surface melt"        , "m/s"    ,   & 
     297      &      jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 
     298 
     299      CALL histdef( kid, "sithicat", "Ice thickness"         , "m"      ,   & 
     300      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     301      CALL histdef( kid, "siconcat", "Ice concentration"     , "%"      ,   & 
     302      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     303      CALL histdef( kid, "sisalcat", "Ice salinity"           , ""      ,   & 
     304      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     305      CALL histdef( kid, "sitemcat", "Ice temperature"       , "C"      ,   & 
     306      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     307      CALL histdef( kid, "snthicat", "Snw thickness"         , "m"      ,   & 
     308      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
     309      CALL histdef( kid, "sntemcat", "Snw temperature"       , "C"      ,   & 
     310      &      jpi, jpj, kh_i, jpl, 1, jpl, nz_i, 32, "inst(x)", rdt, rdt ) 
    369311 
    370312      CALL histend( kid, snc4set )   ! end of the file definition 
    371313 
    372       CALL histwrite( kid, "iicethic", kt, icethi        , jpi*jpj, (/1/) )     
    373       CALL histwrite( kid, "iiceconc", kt, at_i          , jpi*jpj, (/1/) ) 
    374       CALL histwrite( kid, "iicetemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) 
    375       CALL histwrite( kid, "iicevelu", kt, u_ice          , jpi*jpj, (/1/) ) 
    376       CALL histwrite( kid, "iicevelv", kt, v_ice          , jpi*jpj, (/1/) ) 
    377       CALL histwrite( kid, "iicestru", kt, utau_ice       , jpi*jpj, (/1/) ) 
    378       CALL histwrite( kid, "iicestrv", kt, vtau_ice       , jpi*jpj, (/1/) ) 
    379       CALL histwrite( kid, "iicesflx", kt, qsr , jpi*jpj, (/1/) ) 
    380       CALL histwrite( kid, "iicenflx", kt, qns , jpi*jpj, (/1/) ) 
     314      CALL histwrite( kid, "sithic", kt, htm_i         , jpi*jpj, (/1/) )     
     315      CALL histwrite( kid, "siconc", kt, at_i          , jpi*jpj, (/1/) ) 
     316      CALL histwrite( kid, "sitemp", kt, tm_i - rt0    , jpi*jpj, (/1/) ) 
     317      CALL histwrite( kid, "sivelu", kt, u_ice          , jpi*jpj, (/1/) ) 
     318      CALL histwrite( kid, "sivelv", kt, v_ice          , jpi*jpj, (/1/) ) 
     319      CALL histwrite( kid, "sistru", kt, utau_ice       , jpi*jpj, (/1/) ) 
     320      CALL histwrite( kid, "sistrv", kt, vtau_ice       , jpi*jpj, (/1/) ) 
     321      CALL histwrite( kid, "sisflx", kt, qsr , jpi*jpj, (/1/) ) 
     322      CALL histwrite( kid, "sinflx", kt, qns , jpi*jpj, (/1/) ) 
    381323      CALL histwrite( kid, "isnowpre", kt, sprecip        , jpi*jpj, (/1/) ) 
    382       CALL histwrite( kid, "iicesali", kt, smt_i          , jpi*jpj, (/1/) ) 
    383       CALL histwrite( kid, "iicevolu", kt, vt_i           , jpi*jpj, (/1/) ) 
    384       CALL histwrite( kid, "iicedive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
    385  
    386       CALL histwrite( kid, "iicebopr", kt, wfx_bog        , jpi*jpj, (/1/) ) 
    387       CALL histwrite( kid, "iicedypr", kt, wfx_dyn        , jpi*jpj, (/1/) ) 
    388       CALL histwrite( kid, "iicelapr", kt, wfx_opw        , jpi*jpj, (/1/) ) 
    389       CALL histwrite( kid, "iicesipr", kt, wfx_sni        , jpi*jpj, (/1/) ) 
    390       CALL histwrite( kid, "iicerepr", kt, wfx_res        , jpi*jpj, (/1/) ) 
    391       CALL histwrite( kid, "iicebome", kt, wfx_bom        , jpi*jpj, (/1/) ) 
    392       CALL histwrite( kid, "iicesume", kt, wfx_sum        , jpi*jpj, (/1/) ) 
    393       CALL histwrite( kid, "iisfxdyn", kt, sfx_dyn        , jpi*jpj, (/1/) ) 
    394       CALL histwrite( kid, "iisfxres", kt, sfx_res        , jpi*jpj, (/1/) ) 
     324      CALL histwrite( kid, "sisali", kt, smt_i          , jpi*jpj, (/1/) ) 
     325      CALL histwrite( kid, "sivolu", kt, vt_i           , jpi*jpj, (/1/) ) 
     326      CALL histwrite( kid, "sidive", kt, divu_i*1.0e8   , jpi*jpj, (/1/) ) 
     327 
     328      CALL histwrite( kid, "vfxbog", kt, wfx_bog        , jpi*jpj, (/1/) ) 
     329      CALL histwrite( kid, "vfxdyn", kt, wfx_dyn        , jpi*jpj, (/1/) ) 
     330      CALL histwrite( kid, "vfxopw", kt, wfx_opw        , jpi*jpj, (/1/) ) 
     331      CALL histwrite( kid, "vfxsni", kt, wfx_sni        , jpi*jpj, (/1/) ) 
     332      CALL histwrite( kid, "vfxres", kt, wfx_res        , jpi*jpj, (/1/) ) 
     333      CALL histwrite( kid, "vfxbom", kt, wfx_bom        , jpi*jpj, (/1/) ) 
     334      CALL histwrite( kid, "vfxsum", kt, wfx_sum        , jpi*jpj, (/1/) ) 
     335 
     336      CALL histwrite( kid, "sithicat", kt, ht_i        , jpi*jpj*jpl, (/1/) )     
     337      CALL histwrite( kid, "siconcat", kt, a_i         , jpi*jpj*jpl, (/1/) )     
     338      CALL histwrite( kid, "sisalcat", kt, sm_i        , jpi*jpj*jpl, (/1/) )     
     339      CALL histwrite( kid, "sitemcat", kt, tm_i - rt0  , jpi*jpj*jpl, (/1/) )     
     340      CALL histwrite( kid, "snthicat", kt, ht_s        , jpi*jpj*jpl, (/1/) )     
     341      CALL histwrite( kid, "sntemcat", kt, tm_su - rt0 , jpi*jpj*jpl, (/1/) )     
    395342 
    396343      ! Close the file 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90

    r5628 r7067  
    3838   PUBLIC   dia_hsb        ! routine called by step.F90 
    3939   PUBLIC   dia_hsb_init   ! routine called by nemogcm.F90 
    40    PUBLIC   dia_hsb_rst    ! routine called by step.F90 
    4140 
    4241   LOGICAL, PUBLIC ::   ln_diahsb   !: check the heat and salt budgets 
     
    8685      !!--------------------------------------------------------------------------- 
    8786      IF( nn_timing == 1 )   CALL timing_start('dia_hsb')       
     87      ! 
    8888      CALL wrk_alloc( jpi,jpj,   z2d0, z2d1 ) 
    8989      ! 
     
    174174      ENDDO 
    175175 
    176       ! Substract forcing from heat content, salt content and volume variations 
     176      ! ------------------------ ! 
     177      ! 3 -  Drifts              ! 
     178      ! ------------------------ ! 
    177179      zdiff_v1 = zdiff_v1 - frc_v 
    178180      IF( lk_vvl )   zdiff_v2 = zdiff_v2 - frc_v 
     
    187189 
    188190      ! ----------------------- ! 
    189       ! 3 - Diagnostics writing ! 
     191      ! 4 - Diagnostics writing ! 
    190192      ! ----------------------- ! 
    191193      zvol_tot = 0._wp                    ! total ocean volume (calculated with scale factors) 
     
    200202!!gm end 
    201203 
     204      CALL iom_put(   'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
     205      CALL iom_put(   'bgfrctem' , frc_t    * rau0 * rcp * 1.e-20 )   ! hc  - surface forcing (1.e20 J)  
     206      CALL iom_put(   'bgfrchfx' , frc_t    * rau0 * rcp /  &         ! hc  - surface forcing (W/m2)  
     207         &                       ( surf_tot * kt * rdt )        ) 
     208      CALL iom_put(   'bgfrcsal' , frc_s    * 1.e-9    )              ! sc  - surface forcing (psu*km3)  
     209 
    202210      IF( lk_vvl ) THEN 
    203         CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature variation (C)  
    204         CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    variation (psu) 
    205         CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content variation (1.e20 J)  
    206         CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content variation (psu*km3) 
    207         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)   
    208         CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t variation (km3)   
    209         CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
    210         CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)  
    211         CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu)  
     211        CALL iom_put( 'bgtemper' , zdiff_hc / zvol_tot )              ! Temperature drift     (C)  
     212        CALL iom_put( 'bgsaline' , zdiff_sc / zvol_tot )              ! Salinity    drift     (pss) 
     213        CALL iom_put( 'bgheatco' , zdiff_hc * 1.e-20 * rau0 * rcp )   ! Heat content drift    (1.e20 J)  
     214        CALL iom_put( 'bgheatfx' , zdiff_hc * rau0 * rcp /  &         ! Heat flux drift       (W/m2)  
     215           &                       ( surf_tot * kt * rdt )        ) 
     216        CALL iom_put( 'bgsaltco' , zdiff_sc * 1.e-9    )              ! Salt content drift    (psu*km3) 
     217        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)   
     218        CALL iom_put( 'bgvole3t' , zdiff_v2 * 1.e-9    )              ! volume e3t drift      (km3)   
    212219      ELSE 
    213         CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content variation (C)  
    214         CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content variation (psu) 
    215         CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content variation (1.e20 J)  
    216         CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content variation (psu*km3) 
    217         CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh variation (km3)   
    218         CALL iom_put( 'bgfrcvol' , frc_v    * 1.e-9    )              ! vol - surface forcing (km3)  
    219         CALL iom_put( 'bgfrctem' , frc_t / zvol_tot    )              ! hc  - surface forcing (C)  
    220         CALL iom_put( 'bgfrcsal' , frc_s / zvol_tot    )              ! sc  - surface forcing (psu)  
     220        CALL iom_put( 'bgtemper' , zdiff_hc1 / zvol_tot)              ! Heat content drift    (C)  
     221        CALL iom_put( 'bgsaline' , zdiff_sc1 / zvol_tot)              ! Salt content drift    (pss) 
     222        CALL iom_put( 'bgheatco' , zdiff_hc1 * 1.e-20 * rau0 * rcp )  ! Heat content drift    (1.e20 J)  
     223        CALL iom_put( 'bgheatfx' , zdiff_hc1 * rau0 * rcp /  &        ! Heat flux drift       (W/m2)  
     224           &                       ( surf_tot * kt * rdt )         ) 
     225        CALL iom_put( 'bgsaltco' , zdiff_sc1 * 1.e-9    )             ! Salt content drift    (psu*km3) 
     226        CALL iom_put( 'bgvolssh' , zdiff_v1 * 1.e-9    )              ! volume ssh drift      (km3)   
    221227        CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot )              ! hc  - error due to free surface (C) 
    222228        CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot )              ! sc  - error due to free surface (psu) 
     
    244250     ! 
    245251     INTEGER ::   ji, jj, jk   ! dummy loop indices 
    246      INTEGER ::   id1          ! local integers 
    247252     !!---------------------------------------------------------------------- 
    248253     ! 
    249254     IF( TRIM(cdrw) == 'READ' ) THEN        ! Read/initialise  
    250255        IF( ln_rstart ) THEN                   !* Read the restart file 
    251            !id1 = iom_varid( numror, 'frc_vol'  , ldstop = .FALSE. ) 
    252256           ! 
    253257           IF(lwp) WRITE(numout,*) '~~~~~~~' 
     
    261265              CALL iom_get( numror, 'frc_wn_s', frc_wn_s ) 
    262266           ENDIF 
    263            CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini ) 
    264            CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini ) 
    265            CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini ) 
    266            CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini ) 
     267           CALL iom_get( numror, jpdom_autoglo, 'ssh_ini', ssh_ini(:,:) ) 
     268           CALL iom_get( numror, jpdom_autoglo, 'e3t_ini', e3t_ini(:,:,:) ) 
     269           CALL iom_get( numror, jpdom_autoglo, 'hc_loc_ini', hc_loc_ini(:,:,:) ) 
     270           CALL iom_get( numror, jpdom_autoglo, 'sc_loc_ini', sc_loc_ini(:,:,:) ) 
    267271           IF( .NOT. lk_vvl ) THEN 
    268               CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini ) 
    269               CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini ) 
     272              CALL iom_get( numror, jpdom_autoglo, 'ssh_hc_loc_ini', ssh_hc_loc_ini(:,:) ) 
     273              CALL iom_get( numror, jpdom_autoglo, 'ssh_sc_loc_ini', ssh_sc_loc_ini(:,:) ) 
    270274           ENDIF 
    271275       ELSE 
     
    312316           CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s ) 
    313317        ENDIF 
    314         CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini ) 
    315         CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini ) 
    316         CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini ) 
    317         CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini ) 
     318        CALL iom_rstput( kt, nitrst, numrow, 'ssh_ini', ssh_ini(:,:) ) 
     319        CALL iom_rstput( kt, nitrst, numrow, 'e3t_ini', e3t_ini(:,:,:) ) 
     320        CALL iom_rstput( kt, nitrst, numrow, 'hc_loc_ini', hc_loc_ini(:,:,:) ) 
     321        CALL iom_rstput( kt, nitrst, numrow, 'sc_loc_ini', sc_loc_ini(:,:,:) ) 
    318322        IF( .NOT. lk_vvl ) THEN 
    319            CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini ) 
    320            CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini ) 
     323           CALL iom_rstput( kt, nitrst, numrow, 'ssh_hc_loc_ini', ssh_hc_loc_ini(:,:) ) 
     324           CALL iom_rstput( kt, nitrst, numrow, 'ssh_sc_loc_ini', ssh_sc_loc_ini(:,:) ) 
    321325        ENDIF 
     326 
    322327        ! 
    323328     ENDIF 
     
    338343      !!             - Compute coefficients for conversion 
    339344      !!--------------------------------------------------------------------------- 
    340       INTEGER ::   jk       ! dummy loop indice 
    341345      INTEGER ::   ierror   ! local integer 
    342346      INTEGER ::   ios 
     
    344348      NAMELIST/namhsb/ ln_diahsb 
    345349      !!---------------------------------------------------------------------- 
    346  
    347       IF(lwp) THEN 
    348          WRITE(numout,*) 
    349          WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
    350          WRITE(numout,*) '~~~~~~~~ ' 
    351       ENDIF 
    352350 
    353351      REWIND( numnam_ref )              ! Namelist namhsb in reference namelist 
     
    360358      IF(lwm) WRITE ( numond, namhsb ) 
    361359 
    362       ! 
    363       IF(lwp) THEN                   ! Control print 
     360      IF(lwp) THEN 
    364361         WRITE(numout,*) 
    365          WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets' 
    366          WRITE(numout,*) '~~~~~~~~~~~~' 
    367          WRITE(numout,*) '   Namelist namhsb : set hsb parameters' 
    368          WRITE(numout,*) '      Switch for hsb diagnostic (T) or not (F)  ln_diahsb  = ', ln_diahsb 
    369          WRITE(numout,*) 
    370       ENDIF 
    371  
     362         WRITE(numout,*) 'dia_hsb_init' 
     363         WRITE(numout,*) '~~~~~~~~ ' 
     364         WRITE(numout,*) '  check the heat and salt budgets (T) or not (F)       ln_diahsb = ', ln_diahsb 
     365      ENDIF 
     366      ! 
    372367      IF( .NOT. ln_diahsb )   RETURN 
    373368         !      IF( .NOT. lk_mpp_rep ) & 
     
    382377         &      e3t_ini(jpi,jpj,jpk), surf(jpi,jpj),  ssh_ini(jpi,jpj), STAT=ierror ) 
    383378      IF( ierror > 0 ) THEN 
    384          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
    385       ENDIF 
    386  
    387       IF(.NOT. lk_vvl ) ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj),STAT=ierror ) 
    388       IF( ierror > 0 ) THEN 
    389          CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' )   ;   RETURN 
     379         CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) 
     380         RETURN 
     381      ENDIF 
     382 
     383      IF( .NOT. lk_vvl ) THEN 
     384         ALLOCATE( ssh_hc_loc_ini(jpi,jpj), ssh_sc_loc_ini(jpi,jpj), STAT=ierror ) 
     385         IF( ierror > 0 )   THEN 
     386            CALL ctl_stop( 'dia_hsb: unable to allocate hc_loc_ini' ) 
     387            RETURN 
     388         ENDIF 
    390389      ENDIF 
    391390 
     
    393392      ! 2 - Time independant variables and file opening ! 
    394393      ! ----------------------------------------------- ! 
    395       IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated" 
    396       IF(lwp) WRITE(numout,*) '~~~~~~~' 
    397394      surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:)      ! masked surface grid cell area 
    398       surf_tot  = glob_sum( surf(:,:) )                                       ! total ocean surface area 
     395      surf_tot  = glob_sum( surf(:,:) )                   ! total ocean surface area 
    399396 
    400397      IF( lk_bdy ) CALL ctl_warn( 'dia_hsb does not take open boundary fluxes into account' )          
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/SBC/sbccpl.F90

    r6721 r7067  
    16071607      ENDIF 
    16081608 
    1609       !! clem: we should output qemp_oce and qemp_ice (at least) 
    1610       IF( iom_use('hflx_snow_cea') )   CALL iom_put( 'hflx_snow_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) )   ! heat flux from snow (cell average) 
    1611       !! these diags are not outputed yet 
    1612 !!      IF( iom_use('hflx_rain_cea') )   CALL iom_put( 'hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )   ! heat flux from rain (cell average) 
    1613 !!      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put( 'hflx_snow_ao_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average) 
    1614 !!      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put( 'hflx_snow_ai_cea', sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) ) ! heat flux from snow (cell average) 
     1609      ! some more outputs 
     1610      IF( iom_use('hflx_snow_cea') )    CALL iom_put('hflx_snow_cea',   sprecip(:,:) * ( zcptn(:,:) - Lfus ) )                       ! heat flux from snow (cell average) 
     1611      IF( iom_use('hflx_rain_cea') )    CALL iom_put('hflx_rain_cea', ( tprecip(:,:) - sprecip(:,:) ) * zcptn(:,:) )                 ! heat flux from rain (cell average) 
     1612      IF( iom_use('hflx_snow_ao_cea') ) CALL iom_put('hflx_snow_ao_cea',sprecip(:,:) * ( zcptn(:,:) - Lfus ) * (1._wp - zsnw(:,:)) ) ! heat flux from snow (cell average) 
     1613      IF( iom_use('hflx_snow_ai_cea') ) CALL iom_put('hflx_snow_ai_cea',sprecip(:,:) * ( zcptn(:,:) - Lfus ) * zsnw(:,:) )           ! heat flux from snow (cell average) 
    16151614 
    16161615#else 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90

    r6399 r7067  
    229229         CALL lim_sbc_flx( kt )                     ! Update surface ocean mass, heat and salt fluxes 
    230230         ! 
    231          IF(ln_limdiaout) CALL lim_diahsb           ! Diagnostics and outputs  
     231         IF(ln_limdiaout) CALL lim_diahsb( kt )     ! Diagnostics and outputs  
    232232         ! 
    233233         CALL lim_wri( 1 )                          ! Ice outputs  
     
    310310         numit = nit000 - 1 
    311311      ENDIF 
    312       CALL lim_var_agg(1) 
     312      CALL lim_var_agg(2) 
    313313      CALL lim_var_glo2eqv 
    314314      ! 
    315315      CALL lim_sbc_init                 ! ice surface boundary condition    
     316      ! 
     317      IF( ln_limdiaout) CALL lim_diahsb_init  ! initialization for diags 
    316318      ! 
    317319      fr_i(:,:)     = at_i(:,:)         ! initialisation of sea-ice fraction 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90

    r4624 r7067  
    162162                  &                               + avmv(ji,jj,jk) + avmv(ji,jj-1,jk)  )   & 
    163163                  &          + avtb(jk) * tmask(ji,jj,jk) 
    164                !                                            ! Add the background coefficient on eddy viscosity 
     164            END DO 
     165         END DO 
     166         DO jj = 2, jpjm1                                   ! Add the background coefficient on eddy viscosity 
     167            DO ji = 2, jpim1 
    165168               avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk) 
    166169               avmv(ji,jj,jk) = avmv(ji,jj,jk) + avmb(jk) * vmask(ji,jj,jk) 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/step.F90

    r6405 r7067  
    337337      IF( lk_vvl           )   CALL dom_vvl_sf_swp( kstp )  ! swap of vertical scale factors 
    338338      ! 
     339      IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics 
     340 
    339341      IF( lrst_oce         )   CALL rst_write( kstp )       ! write output ocean restart file 
    340342      IF( ln_sto_eos       )   CALL sto_rst_write( kstp )   ! write restart file for stochastic parameters 
     
    351353      ENDIF 
    352354#endif 
    353       IF( ln_diahsb        )   CALL dia_hsb( kstp )         ! - ML - global conservation diagnostics 
    354       IF( lk_diaobs  )         CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
     355      IF( lk_diaobs        )   CALL dia_obs( kstp )         ! obs-minus-model (assimilation) diagnostics (call after dynamics update) 
    355356 
    356357      !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90

    r6943 r7067  
    251251      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout)            ::  pe1 , pe2 , pe3   !  PAR ( R-G-B) 
    252252      REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL  ::  pe0  
    253       REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(out)  , OPTIONAL  ::  pqsr100   
     253      REAL(wp), DIMENSION(jpi,jpj)    , INTENT(out)  , OPTIONAL  ::  pqsr100   
    254254      !! * local variables 
    255255      INTEGER    ::   ji, jj, jk     ! dummy loop indices 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90

    r6941 r7067  
    6969                                CALL trc_adv( kstp )            ! horizontal & vertical advection  
    7070         IF( ln_zps ) THEN 
    71            IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kt, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom 
    72            ELSE                 ; CALL zps_hde    ( kt, jptra, trb, gtru, gtrv )                                      !  only bottom 
     71           IF( ln_isfcav ) THEN ; CALL zps_hde_isf( kstp, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi )  ! both top & bottom 
     72           ELSE                 ; CALL zps_hde    ( kstp, jptra, trb, gtru, gtrv )                                      !  only bottom 
    7373           ENDIF 
    7474         ENDIF 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcrst.F90

    r5513 r7067  
    304304         IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift 
    305305      END DO 
    306       WRITE(numout,*)  
     306      IF(lwp) WRITE(numout,*)  
    3073079000  FORMAT(' tracer nb :',i2,'    name :',a10,'    mean :',e18.10,'    min :',e18.10, & 
    308308      &      '    max :',e18.10,'    drift :',e18.10, ' %') 
  • branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcstp.F90

    r6952 r7067  
    3333   REAL(wp) :: rdt_sampl 
    3434   INTEGER  :: nb_rec_per_day 
    35    INTEGER  :: isecfst, iseclast 
     35   REAL(wp) :: rsecfst, rseclast 
    3636   LOGICAL  :: llnew 
    3737 
     
    131131      INTEGER, INTENT(in) ::   kt 
    132132      INTEGER  :: jn 
    133       REAL(wp) :: zsecfst 
     133      REAL(wp) :: zkt 
    134134      CHARACTER(len=1)               ::   cl1                      ! 1 character 
    135135      CHARACTER(len=2)               ::   cl2                      ! 2 characters 
     
    137137      IF( kt == nittrc000 ) THEN 
    138138         IF( ln_cpl )  THEN   
    139             rdt_sampl = 86400. / ncpl_qsr_freq 
     139            rdt_sampl = rday / ncpl_qsr_freq 
    140140            nb_rec_per_day = ncpl_qsr_freq 
    141141         ELSE   
    142             rdt_sampl = MAX( 3600., rdt * nn_dttrc ) 
    143             nb_rec_per_day = INT( 86400 / rdt_sampl ) 
     142            rdt_sampl = MAX( 3600., rdttrc(1) ) 
     143            nb_rec_per_day = INT( rday / rdt_sampl ) 
    144144         ENDIF 
    145145         ! 
     
    155155         IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. & 
    156156                            iom_varid( numrtr, 'qsr_arr_1', ldstop = .FALSE. ) > 0 .AND. & 
    157                             iom_varid( numrtr, 'zsecfst'  , ldstop = .FALSE. ) > 0 ) THEN  
    158             IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file' 
     157                            iom_varid( numrtr, 'ktdcy'    , ldstop = .FALSE. ) > 0 ) THEN  
     158            CALL iom_get( numrtr, 'ktdcy', zkt )   !  A mean of qsr 
     159            rsecfst = INT( zkt ) * rdttrc(1) 
     160            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean read in the restart file at time-step rsecfst =', rsecfst, ' s ' 
    159161            CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean )   !  A mean of qsr 
    160             CALL iom_get( numrtr, 'zsecfst', zsecfst )   !  A mean of qsr 
    161             isecfst = INT( zsecfst ) 
    162162            DO jn = 1, nb_rec_per_day  
    163163             IF( jn <= 9 )  THEN 
     
    171171         ELSE                                         !* no restart: set from nit000 values 
    172172            IF(lwp) WRITE(numout,*) 'trc_qsr_mean:   qsr_mean set to nit000 values' 
    173             isecfst  = nsec_year + nsec1jan000   !   number of seconds between Jan. 1st 00h of nit000 year and the middle of time step 
     173            rsecfst  = kt * rdttrc(1) 
    174174            ! 
    175175            qsr_mean(:,:) = qsr(:,:) 
     
    181181      ENDIF 
    182182      ! 
    183       iseclast = nsec_year + nsec1jan000 
    184       ! 
    185       llnew   = ( iseclast - isecfst )  > INT( rdt_sampl )   !   new shortwave to store 
     183      rseclast = kt * rdttrc(1) 
     184      ! 
     185      llnew   = ( rseclast - rsecfst ) .ge.  rdt_sampl    !   new shortwave to store 
    186186      IF( llnew ) THEN 
    187187          IF( lwp ) WRITE(numout,*) ' New shortwave to sample for TOP at time kt = ', kt, & 
    188              &                      ' time = ', (iseclast+rdt*nn_dttrc/2.)/3600.,'hours ' 
    189           isecfst = iseclast 
     188             &                      ' time = ', rseclast/3600.,'hours ' 
     189          rsecfst = rseclast 
    190190          DO jn = 1, nb_rec_per_day - 1 
    191191             qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1) 
     
    199199         IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file  kt =', kt 
    200200         IF(lwp) WRITE(numout,*) '~~~~~~~' 
     201         zkt = REAL( kt, wp ) 
     202         CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 
    201203          DO jn = 1, nb_rec_per_day  
    202204             IF( jn <= 9 )  THEN 
     
    209211         ENDDO 
    210212         CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) ) 
    211          zsecfst = REAL( isecfst, wp ) 
    212          CALL iom_rstput( kt, nitrst, numrtw, 'zsecfst', zsecfst ) 
    213213      ENDIF 
    214214      ! 
Note: See TracChangeset for help on using the changeset viewer.