Changeset 7067 for branches/2015
- Timestamp:
- 2016-10-21T17:35:49+02:00 (8 years ago)
- Location:
- branches/2015/dev_r6946_Offline_vvl/NEMOGCM
- Files:
-
- 19 edited
-
CONFIG/ORCA2_LIM3/EXP00/iodef.xml (modified) (6 diffs)
-
CONFIG/SHARED/field_def.xml (modified) (4 diffs)
-
NEMO/LIM_SRC_3/ice.F90 (modified) (4 diffs)
-
NEMO/LIM_SRC_3/limcons.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limdiahsb.F90 (modified) (7 diffs)
-
NEMO/LIM_SRC_3/limitd_me.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limrhg.F90 (modified) (13 diffs)
-
NEMO/LIM_SRC_3/limsbc.F90 (modified) (1 diff)
-
NEMO/LIM_SRC_3/limvar.F90 (modified) (6 diffs)
-
NEMO/LIM_SRC_3/limwri.F90 (modified) (9 diffs)
-
NEMO/OPA_SRC/DIA/diahsb.F90 (modified) (13 diffs)
-
NEMO/OPA_SRC/SBC/sbccpl.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/SBC/sbcice_lim.F90 (modified) (2 diffs)
-
NEMO/OPA_SRC/ZDF/zdfric.F90 (modified) (1 diff)
-
NEMO/OPA_SRC/step.F90 (modified) (2 diffs)
-
NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90 (modified) (1 diff)
-
NEMO/TOP_SRC/TRP/trctrp.F90 (modified) (1 diff)
-
NEMO/TOP_SRC/trcrst.F90 (modified) (1 diff)
-
NEMO/TOP_SRC/trcstp.F90 (modified) (8 diffs)
Legend:
- Unmodified
- Added
- Removed
-
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml
r6316 r7067 88 88 <field field_ref="u_masstr" name="vozomatr" /> 89 89 <field field_ref="u_heattr" name="sozohetr" /> 90 90 <field field_ref="u_salttr" name="sozosatr" /> 91 91 </file> 92 92 … … 99 99 <field field_ref="v_masstr" name="vomematr" /> 100 100 <field field_ref="v_heattr" name="somehetr" /> 101 101 <field field_ref="v_salttr" name="somesatr" /> 102 102 </file> 103 103 … … 116 116 <field field_ref="iceconc" name="siconc" /> 117 117 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" />125 118 <field field_ref="vfxice" name="vfxice" /> 126 119 <field field_ref="vfxsnw" name="vfxsnw" /> 127 120 <field field_ref="vfxsub" name="vfxsub" /> 121 <field field_ref="vfxsub_err" name="vfxsub_err" /> 128 122 <field field_ref="vfxspr" name="vfxspr" /> 129 123 … … 134 128 <field field_ref="destrp" name="destrp" /> 135 129 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" />144 130 <field field_ref="sfx" name="sfx" /> 145 131 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" />151 132 <field field_ref="hfxout" name="hfxout" /> 152 133 <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 164 135 165 136 <!-- diags --> 166 <field field_ref="hfxdhc" name="hfxdhc" />167 <field field_ref="hfxtur" name="hfxtur" />168 169 137 <field field_ref="isst" name="sst" /> 170 138 <field field_ref="isss" name="sss" /> … … 209 177 <field field_ref="bgsaline" name="bgsaline" /> 210 178 <field field_ref="bgheatco" name="bgheatco" /> 179 <field field_ref="bgheatfx" name="bgheatfx" /> 211 180 <field field_ref="bgsaltco" name="bgsaltco" /> 212 181 <field field_ref="bgvolssh" name="bgvolssh" /> … … 214 183 <field field_ref="bgfrcvol" name="bgfrcvol" /> 215 184 <field field_ref="bgfrctem" name="bgfrctem" /> 185 <field field_ref="bgfrchfx" name="bgfrchfx" /> 216 186 <field field_ref="bgfrcsal" name="bgfrcsal" /> 217 187 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" /> 268 207 269 208 </file> -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/CONFIG/SHARED/field_def.xml
r6471 r7067 225 225 <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" /> 226 226 <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" /> 228 230 <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" /> 229 231 <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" /> … … 331 333 <field id="vfxsnw" long_name="snw melt/growth" unit="m/day" /> 332 334 <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" /> 333 336 <field id="vfxspr" long_name="snw precipitation on ice" unit="m/day" /> 334 337 <field id="vfxthin" long_name="daily thermo ice prod. for thin ice(<20cm) + open water" unit="m/day" /> … … 499 502 <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" /> 500 503 <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" /> 503 507 <field id="bgvolssh" long_name="drift in global mean ssh volume wrt timestep 1" unit="km3" /> 504 508 <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" /> 513 515 </field_group> 514 516 … … 517 519 <field_group id="SBC_scalar" domain_ref="1point" > 518 520 <!-- 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" /> 570 540 </field_group> 571 541 -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r6477 r7067 243 243 ! 244 244 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)246 245 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_bo !: Sea-Ice bottom temperature [Kelvin] 247 246 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frld !: Leads fraction = 1 - ice fraction … … 320 319 ! ! this is an extensive variable that has to be transported 321 320 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)323 321 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 324 323 325 324 !! Variables summed over all categories, or associated to all the ice in a single grid cell 326 325 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)328 326 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_i , vt_s !: ice and snow total volume per unit area (m) 329 327 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i !: ice total fractional area (ice concentration) 330 328 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ato_i !: =1-at_i ; total open water fractional area 331 329 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 336 337 337 338 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_s !: Snow temperatures [K] … … 435 436 436 437 ii = ii + 1 437 ALLOCATE( sist (jpi,jpj) , icethi (jpi,jpj) ,t_bo (jpi,jpj) , &438 ALLOCATE( sist (jpi,jpj) , t_bo (jpi,jpj) , & 438 439 & frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) , & 439 440 & wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , & … … 456 457 & v_s (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , & 457 458 & sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , & 458 & o v_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) , & 461 462 & 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) ) 464 466 ii = ii + 1 465 467 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 288 288 #if ! defined key_bdy 289 289 ! 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 ) 292 293 ! salt flux 293 294 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 31 31 32 32 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 37 38 !! * Substitutions 38 39 # include "vectopt_loop_substitute.h90" … … 46 47 CONTAINS 47 48 48 SUBROUTINE lim_diahsb 49 SUBROUTINE lim_diahsb( kt ) 49 50 !!--------------------------------------------------------------------------- 50 51 !! *** ROUTINE lim_diahsb *** … … 53 54 !! 54 55 !!--------------------------------------------------------------------------- 56 INTEGER, INTENT(in) :: kt ! number of iteration 55 57 !! 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 67 61 !!--------------------------------------------------------------------------- 68 62 IF( nn_timing == 1 ) CALL timing_start('lim_diahsb') 69 63 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) 146 73 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) 215 133 ! 216 134 IF( lrst_ice ) CALL lim_diahsb_rst( numit, 'WRITE' ) 217 135 ! 218 136 IF( nn_timing == 1 ) CALL timing_stop('lim_diahsb') 219 !137 ! 220 138 END SUBROUTINE lim_diahsb 221 139 … … 233 151 !! - Compute coefficients for conversion 234 152 !!--------------------------------------------------------------------------- 235 INTEGER :: jk ! dummy loop indice236 153 INTEGER :: ierror ! local integer 237 154 !! … … 247 164 WRITE(numout,*) '~~~~~~~~~~~~' 248 165 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 250 173 CALL lim_diahsb_rst( nstart, 'READ' ) !* read or initialize all required files 251 174 ! … … 263 186 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 264 187 ! 265 INTEGER :: id1, id2, id3 ! local integers266 188 !!---------------------------------------------------------------------- 267 189 ! 268 190 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 269 191 IF( ln_rstart ) THEN !* Read the restart file 270 !id1 = iom_varid( numrir, 'frc_vol' , ldstop = .TRUE. )271 192 ! 272 193 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 ) 278 204 ELSE 279 205 IF(lwp) WRITE(numout,*) '~~~~~~~' 280 206 IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state ' 281 207 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 285 219 ENDIF 286 220 … … 288 222 ! ! ------------------- 289 223 IF(lwp) WRITE(numout,*) '~~~~~~~' 290 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp224 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write at it= ', kt,' date= ', ndastp 291 225 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 ) 295 234 ! 296 235 ENDIF -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r6469 r7067 866 866 DO jj = 1, jpj 867 867 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))) 869 869 END DO 870 870 END DO -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r5888 r7067 10 10 !! 3.4 ! 2011-01 (A. Porter) dynamical allocation 11 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 !! 3.6 ! 2016-06 (C. Rousset) Rewriting (conserves energy) 12 13 !!---------------------------------------------------------------------- 13 14 #if defined key_lim3 || ( defined key_lim2 && ! defined key_lim2_vp ) … … 95 96 !! coriolis terms of the momentum equation 96 97 !! 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 99 99 !! which are inputs of the ITD, store stress 100 100 !! for the next time step 101 !! 6) Control prints of residual (convergence)101 !! 5) Control prints of residual (convergence) 102 102 !! and charge ellipse. 103 103 !! The user should make sure that the parameters … … 106 106 !! e.g. in the Canadian Archipelago 107 107 !! 108 !! ** Notes : Boundary condition for ice is chosen no-slip 109 !! but can be adjusted with param rn_shlat 110 !! 108 111 !! References : Hunke and Dukowicz, JPO97 109 112 !! Bouillon et al., Ocean Modelling 2009 … … 115 118 INTEGER :: jter ! local integers 116 119 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 144 131 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) 156 157 !!------------------------------------------------------------------- 157 158 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 ) 162 164 163 165 #if defined key_lim2 && ! defined key_lim2_vp … … 176 178 ! 177 179 !------------------------------------------------------------------------------! 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) 189 229 #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 195 245 #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 ) ) 203 257 END DO 204 258 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 218 260 ! 219 261 !------------------------------------------------------------------------------! 220 262 ! 2) Wind / ocean stress, mass terms, coriolis terms 221 263 !------------------------------------------------------------------------------! 222 !223 ! Wind stress, coriolis and mass terms on the sides of the squares224 ! zfrld1: lead fraction on U-points225 ! zfrld2: lead fraction on V-points226 ! zmass1: ice/snow mass on U-points227 ! zmass2: ice/snow mass on V-points228 ! zcorl1: Coriolis parameter on U-points229 ! zcorl2: Coriolis parameter on V-points230 ! (ztagnx,ztagny): wind stress on U/V points231 ! v_oce1: ocean v component on u points232 ! u_oce2: ocean u component on v points233 264 234 265 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! … … 242 273 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 243 274 ! 244 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:)) * r1_rau0275 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 245 276 ! 246 277 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! … … 251 282 DO ji = fs_2, fs_jpim1 252 283 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 298 324 299 325 END DO 300 326 END DO 301 327 CALL lbc_lnk( zmf, 'T', 1. ) 302 328 ! 303 329 !------------------------------------------------------------------------------! … … 305 331 !------------------------------------------------------------------------------! 306 332 ! 307 ! Time step for subcycling308 dtevp = rdt_ice / nn_nevp309 #if defined key_lim3310 dtotel = dtevp / ( 2._wp * rn_relast * rdt_ice )311 #else312 dtotel = dtevp / ( 2._wp * telast )313 #endif314 z1_dtotel = 1._wp / ( 1._wp + dtotel )315 z1_dtevp = 1._wp / dtevp316 !-ecc2: square of yield ellipse eccenticrity (reminder: must become a namelist parameter)317 ecc2 = rn_ecc * rn_ecc318 ecci = 1. / ecc2319 320 !-Initialise stress tensor321 zs1 (:,:) = stress1_i (:,:)322 zs2 (:,:) = stress2_i (:,:)323 zs12(:,:) = stress12_i(:,:)324 325 333 ! !----------------------! 326 334 DO jter = 1 , nn_nevp ! loop over jter ! 327 335 ! !----------------------! 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 364 348 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) & 365 349 & + ( 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 382 356 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. ) 417 402 418 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002)403 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 419 404 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 431 431 END DO 432 432 END DO 433 433 ! 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 440 439 DO jj = k_j1+1, k_jpj-1 441 440 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) 455 459 END DO 456 460 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 459 494 #if defined key_agrif && defined key_lim2 460 495 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 461 496 #endif 462 497 #if defined key_bdy 463 CALL bdy_ice_lim_dyn( 'U' )498 CALL bdy_ice_lim_dyn( 'U' ) 464 499 #endif 500 501 ELSE ! odd iterations 465 502 466 503 DO jj = k_j1+1, k_jpj-1 467 504 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) 481 523 END DO 482 524 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 485 558 #if defined key_agrif && defined key_lim2 486 559 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 487 560 #endif 488 561 #if defined key_bdy 489 CALL bdy_ice_lim_dyn( 'V' )562 CALL bdy_ice_lim_dyn( 'V' ) 490 563 #endif 491 564 492 ELSE 565 ENDIF 566 567 IF(ln_ctl) THEN ! Convergence test 493 568 DO jj = k_j1+1, k_jpj-1 494 DO ji = fs_2, fs_jpim1495 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1)496 z0 = zmass2(ji,jj) * z1_dtevp497 ! SB modif because ocean has no slip boundary condition498 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 + za506 zccb = zcorl2(ji,jj)507 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch508 END DO509 END DO510 511 CALL lbc_lnk( v_ice(:,:), 'V', -1. )512 #if defined key_agrif && defined key_lim2513 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' )514 #endif515 #if defined key_bdy516 CALL bdy_ice_lim_dyn( 'V' )517 #endif518 519 DO jj = k_j1+1, k_jpj-1520 DO ji = fs_2, fs_jpim1521 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1)522 z0 = zmass1(ji,jj) * z1_dtevp523 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 + za531 zccb = zcorl1(ji,jj)532 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch533 END DO534 END DO535 536 CALL lbc_lnk( u_ice(:,:), 'U', -1. )537 #if defined key_agrif && defined key_lim2538 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' )539 #endif540 #if defined key_bdy541 CALL bdy_ice_lim_dyn( 'U' )542 #endif543 544 ENDIF545 546 IF(ln_ctl) THEN547 !--- Convergence test.548 DO jj = k_j1+1 , k_jpj-1549 569 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 570 END DO … … 552 572 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 553 573 ENDIF 554 574 ! 555 575 ! ! ==================== ! 556 576 END DO ! end loop over jter ! … … 558 578 ! 559 579 !------------------------------------------------------------------------------! 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 570 621 END DO 571 622 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 --- ! 648 626 stress1_i (:,:) = zs1 (:,:) 649 627 stress2_i (:,:) = zs2 (:,:) … … 652 630 ! 653 631 !------------------------------------------------------------------------------! 654 ! 6) Control prints of residual and charge ellipse632 ! 5) Control prints of residual and charge ellipse 655 633 !------------------------------------------------------------------------------! 656 634 ! … … 675 653 DO ji = 2, jpim1 676 654 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) ) 679 657 WRITE(charout,FMT="('lim_rhg :', I4, I4, D23.16, D23.16, D23.16, D23.16, A10)") 680 658 CALL prt_ctl_info(charout) … … 687 665 ENDIF 688 666 ! 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 ) 693 672 694 673 END SUBROUTINE lim_rhg -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limsbc.F90
r6399 r7067 110 110 !!--------------------------------------------------------------------- 111 111 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 128 113 CALL wrk_alloc( jpi,jpj, zalb ) 129 114 -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limvar.F90
r6469 r7067 54 54 PUBLIC lim_var_eqv2glo 55 55 PUBLIC lim_var_salprof 56 PUBLIC lim_var_icetm57 56 PUBLIC lim_var_bv 58 57 PUBLIC lim_var_salprof1d … … 89 88 ! Compute variables 90 89 !-------------------- 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 ! 110 97 DO jj = 1, jpj 111 98 DO ji = 1, jpi … … 115 102 116 103 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 119 115 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 122 119 DO jl = 1, jpl 120 123 121 DO jj = 1, jpj 124 122 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 135 129 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 139 143 ! 140 144 ENDIF … … 246 250 ! Mean temperature 247 251 !------------------- 248 vt_i (:,:) = 0._wp249 DO jl = 1, jpl250 vt_i(:,:) = vt_i(:,:) + v_i(:,:,jl)251 END DO252 ! 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 ) 252 256 253 257 tm_i(:,:) = 0._wp … … 397 401 END SUBROUTINE lim_var_salprof 398 402 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 405 412 !!------------------------------------------------------------------ 406 413 INTEGER :: ji, jj, jk, jl ! dummy loop indices 407 414 !!------------------------------------------------------------------ 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 416 418 DO jl = 1, jpl 417 419 DO jk = 1, nlay_i 418 420 DO jj = 1, jpj 419 421 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 ) 462 433 END DO 463 434 END DO … … 715 686 zht_i(ji,1:jpl) = 0._wp 716 687 za_i (ji,1:jpl) = 0._wp 717 688 itest(:) = 0 689 718 690 ! *** case very thin ice: fill only category 1 719 691 IF ( i_fill == 1 ) THEN -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/LIM_SRC_3/limwri.F90
r6417 r7067 17 17 USE sbc_oce ! Surface boundary condition: ocean fields 18 18 USE sbc_ice ! Surface boundary condition: ice fields 19 USE dom_ice20 19 USE ice 21 20 USE limvar … … 40 39 !!---------------------------------------------------------------------- 41 40 CONTAINS 42 43 #if defined key_dimgout44 # include "limwri_dimg.h90"45 #else46 41 47 42 SUBROUTINE lim_wri( kindic ) … … 59 54 INTEGER :: ji, jj, jk, jl ! dummy loop indices 60 55 REAL(wp) :: z1_365 61 REAL(wp) :: z tmp62 REAL(wp), POINTER, DIMENSION(:,:,:) :: z oi, zei, zt_i, zt_s63 REAL(wp), POINTER, DIMENSION(:,:) :: z2d, z 2da, z2db, zswi ! 2D workspace56 REAL(wp) :: z2da, z2db, ztmp 57 REAL(wp), POINTER, DIMENSION(:,:,:) :: zswi2 58 REAL(wp), POINTER, DIMENSION(:,:) :: z2d, zswi ! 2D workspace 64 59 !!------------------------------------------------------------------- 65 60 66 61 IF( nn_timing == 1 ) CALL timing_start('limwri') 67 62 68 CALL wrk_alloc( jpi, jpj, jpl, z oi, zei, zt_i, zt_s)69 CALL wrk_alloc( jpi, jpj , z2d, z 2da, z2db, zswi )63 CALL wrk_alloc( jpi, jpj, jpl, zswi2 ) 64 CALL wrk_alloc( jpi, jpj , z2d, zswi ) 70 65 71 66 !----------------------------- … … 74 69 z1_365 = 1._wp / 365._wp 75 70 76 CALL lim_var_icetm ! mean sea ice temperature77 78 CALL lim_var_bv ! brine volume 79 80 DO jj = 1, jpj ! presence indicator of ice71 ! brine volume 72 CALL lim_var_bv 73 74 ! tresholds for outputs 75 DO jj = 1, jpj 81 76 DO ji = 1, jpi 82 77 zswi(ji,jj) = MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi06 ) ) 83 78 END DO 84 79 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 90 82 DO ji = 1, jpi 91 z 2d(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 ) ) 92 84 END DO 93 85 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 105 87 ! 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 106 104 IF ( iom_use( "uice_ipa" ) .OR. iom_use( "vice_ipa" ) .OR. iom_use( "icevel" ) ) THEN 107 105 DO jj = 2 , jpjm1 108 106 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 ) 111 110 END DO 112 111 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 123 116 ENDIF 124 117 ! 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 146 124 ! 147 IF ( iom_use( "icest" ) ) THEN148 z2d(:,:) = 0.e0149 DO jl = 1, jpl150 DO jj = 1, jpj151 DO ji = 1, jpi152 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 DO154 END DO155 END DO156 CALL iom_put( "icest" , z2d ) ! ice surface temperature157 ENDIF158 159 IF ( iom_use( "icecolf" ) ) CALL iom_put( "icecolf", hicol ) ! frazil ice collection thickness160 161 125 CALL iom_put( "isst" , sst_m ) ! sea surface temperature 162 126 CALL iom_put( "isss" , sss_m ) ! sea surface salinity 163 CALL iom_put( "iceconc" , at_i 164 CALL iom_put( "icevolu" , vt_i 165 CALL iom_put( "icehc" , et_i 166 CALL iom_put( "isnowhc" , et_s 167 CALL iom_put( "ibrinv" , bv _i * 100._wp) ! brine volume127 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 168 132 CALL iom_put( "utau_ice" , utau_ice ) ! wind stress over ice along i-axis at I-point 169 133 CALL iom_put( "vtau_ice" , vtau_ice ) ! wind stress over ice along j-axis at I-point 170 134 CALL iom_put( "snowpre" , sprecip * 86400. ) ! snow precipitation 171 CALL iom_put( "micesalt" , smt_i 172 173 CALL iom_put( "icestr" , strength * 0.001 )! ice strength174 CALL iom_put( "idive" , divu_i * 1.0e8 ) 175 CALL iom_put( "ishear" , shear_i * 1.0e8 ) 176 CALL iom_put( "snowvol" , vt_s 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 177 141 178 142 CALL iom_put( "icetrp" , diag_trp_vi * rday ) ! ice volume transport … … 183 147 184 148 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 187 151 CALL iom_put( "sfxsni" , sfx_sni * rday ) ! salt flux from snow ice formation 188 152 CALL iom_put( "sfxopw" , sfx_opw * rday ) ! salt flux from open water formation 189 153 CALL iom_put( "sfxdyn" , sfx_dyn * rday ) ! salt flux from ridging rafting 190 CALL iom_put( "sfxres" , sfx_res * rday ) ! salt flux from residual154 CALL iom_put( "sfxres" , sfx_res * rday ) ! salt flux from limupdate (resultant) 191 155 CALL iom_put( "sfxbri" , sfx_bri * rday ) ! salt flux from brines 192 156 CALL iom_put( "sfxsub" , sfx_sub * rday ) ! salt flux from sublimation … … 202 166 CALL iom_put( "vfxbom" , wfx_bom * ztmp ) ! bottom melt 203 167 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) 204 178 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( "vfxs pr" , 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 208 182 CALL iom_put( "afxtot" , afx_tot * rday ) ! concentration tendency (total) 209 183 CALL iom_put( "afxdyn" , afx_dyn * rday ) ! concentration tendency (dynamics) … … 225 199 CALL iom_put ('hfxdif' , hfx_dif(:,:) ) ! 226 200 CALL iom_put ('hfxopw' , hfx_opw(:,:) ) ! 227 CALL iom_put ('hfxtur' , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3) ) ! turbulent heat flux at ice base201 CALL iom_put ('hfxtur' , fhtur(:,:) * SUM( a_i_b(:,:,:), dim=3 ) ) ! turbulent heat flux at ice base 228 202 CALL iom_put ('hfxdhc' , diag_heat(:,:) ) ! Heat content variation in snow and ice 229 203 CALL iom_put ('hfxspr' , hfx_spr(:,:) ) ! Heat content of snow precip 230 204 231 232 IF ( iom_use( "vfxthin" ) ) THEN ! ice production for open water + thin ice (<20cm) => comparable to observations233 DO jj = 1, jpj234 DO ji = 1, jpi235 z2d(ji,jj) = vt_i(ji,jj) / MAX( at_i(ji,jj), epsi06 ) * zswi(ji,jj) ! mean ice thickness236 END DO237 END DO238 WHERE( z2d(:,:) < 0.2 .AND. z2d(:,:) > 0. ) ; z2da = wfx_bog239 ELSEWHERE ; z2da = 0._wp240 END WHERE241 CALL iom_put( "vfxthin", ( wfx_opw + z2da ) * ztmp )242 ENDIF243 205 244 206 !-------------------------------- 245 207 ! Output values for each category 246 208 !-------------------------------- 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 252 213 ! 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 ) 258 215 ! 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 ) 295 221 296 222 ! ! Create an output files (output.lim.abort.nc) if S < 0 or u > 20 m/s … … 298 224 ! not yet implemented 299 225 300 CALL wrk_dealloc( jpi, jpj, jpl, z oi, 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 ) 302 228 303 229 IF( nn_timing == 1 ) CALL timing_stop('limwri') 304 230 305 231 END SUBROUTINE lim_wri 306 #endif307 232 308 233 … … 319 244 !! 4.0 ! 2013-06 (C. Rousset) 320 245 !!---------------------------------------------------------------------- 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 323 250 !!---------------------------------------------------------------------- 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" , & 342 274 & jpi, jpj, kh_i, 1, 1, 1, -99, 32, "inst(x)", rdt, rdt ) 343 275 CALL histdef( kid, "isnowpre", "Snow precipitation" , "kg/m2/s", & 344 276 & 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" , & 356 289 & 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 ) 369 311 370 312 CALL histend( kid, snc4set ) ! end of the file definition 371 313 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/) ) 381 323 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/) ) 395 342 396 343 ! Close the file -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/DIA/diahsb.F90
r5628 r7067 38 38 PUBLIC dia_hsb ! routine called by step.F90 39 39 PUBLIC dia_hsb_init ! routine called by nemogcm.F90 40 PUBLIC dia_hsb_rst ! routine called by step.F9041 40 42 41 LOGICAL, PUBLIC :: ln_diahsb !: check the heat and salt budgets … … 86 85 !!--------------------------------------------------------------------------- 87 86 IF( nn_timing == 1 ) CALL timing_start('dia_hsb') 87 ! 88 88 CALL wrk_alloc( jpi,jpj, z2d0, z2d1 ) 89 89 ! … … 174 174 ENDDO 175 175 176 ! Substract forcing from heat content, salt content and volume variations 176 ! ------------------------ ! 177 ! 3 - Drifts ! 178 ! ------------------------ ! 177 179 zdiff_v1 = zdiff_v1 - frc_v 178 180 IF( lk_vvl ) zdiff_v2 = zdiff_v2 - frc_v … … 187 189 188 190 ! ----------------------- ! 189 ! 3- Diagnostics writing !191 ! 4 - Diagnostics writing ! 190 192 ! ----------------------- ! 191 193 zvol_tot = 0._wp ! total ocean volume (calculated with scale factors) … … 200 202 !!gm end 201 203 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 202 210 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) 212 219 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) 221 227 CALL iom_put( 'bgmistem' , zerr_hc1 / zvol_tot ) ! hc - error due to free surface (C) 222 228 CALL iom_put( 'bgmissal' , zerr_sc1 / zvol_tot ) ! sc - error due to free surface (psu) … … 244 250 ! 245 251 INTEGER :: ji, jj, jk ! dummy loop indices 246 INTEGER :: id1 ! local integers247 252 !!---------------------------------------------------------------------- 248 253 ! 249 254 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 250 255 IF( ln_rstart ) THEN !* Read the restart file 251 !id1 = iom_varid( numror, 'frc_vol' , ldstop = .FALSE. )252 256 ! 253 257 IF(lwp) WRITE(numout,*) '~~~~~~~' … … 261 265 CALL iom_get( numror, 'frc_wn_s', frc_wn_s ) 262 266 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(:,:,:) ) 267 271 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(:,:) ) 270 274 ENDIF 271 275 ELSE … … 312 316 CALL iom_rstput( kt, nitrst, numrow, 'frc_wn_s', frc_wn_s ) 313 317 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(:,:,:) ) 318 322 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(:,:) ) 321 325 ENDIF 326 322 327 ! 323 328 ENDIF … … 338 343 !! - Compute coefficients for conversion 339 344 !!--------------------------------------------------------------------------- 340 INTEGER :: jk ! dummy loop indice341 345 INTEGER :: ierror ! local integer 342 346 INTEGER :: ios … … 344 348 NAMELIST/namhsb/ ln_diahsb 345 349 !!---------------------------------------------------------------------- 346 347 IF(lwp) THEN348 WRITE(numout,*)349 WRITE(numout,*) 'dia_hsb_init : check the heat and salt budgets'350 WRITE(numout,*) '~~~~~~~~ '351 ENDIF352 350 353 351 REWIND( numnam_ref ) ! Namelist namhsb in reference namelist … … 360 358 IF(lwm) WRITE ( numond, namhsb ) 361 359 362 ! 363 IF(lwp) THEN ! Control print 360 IF(lwp) THEN 364 361 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 ! 372 367 IF( .NOT. ln_diahsb ) RETURN 373 368 ! IF( .NOT. lk_mpp_rep ) & … … 382 377 & e3t_ini(jpi,jpj,jpk), surf(jpi,jpj), ssh_ini(jpi,jpj), STAT=ierror ) 383 378 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 390 389 ENDIF 391 390 … … 393 392 ! 2 - Time independant variables and file opening ! 394 393 ! ----------------------------------------------- ! 395 IF(lwp) WRITE(numout,*) "dia_hsb: heat salt volume budgets activated"396 IF(lwp) WRITE(numout,*) '~~~~~~~'397 394 surf(:,:) = e1t(:,:) * e2t(:,:) * tmask_i(:,:) ! masked surface grid cell area 398 surf_tot = glob_sum( surf(:,:) ) 395 surf_tot = glob_sum( surf(:,:) ) ! total ocean surface area 399 396 400 397 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 1607 1607 ENDIF 1608 1608 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) 1615 1614 1616 1615 #else -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/SBC/sbcice_lim.F90
r6399 r7067 229 229 CALL lim_sbc_flx( kt ) ! Update surface ocean mass, heat and salt fluxes 230 230 ! 231 IF(ln_limdiaout) CALL lim_diahsb 231 IF(ln_limdiaout) CALL lim_diahsb( kt ) ! Diagnostics and outputs 232 232 ! 233 233 CALL lim_wri( 1 ) ! Ice outputs … … 310 310 numit = nit000 - 1 311 311 ENDIF 312 CALL lim_var_agg( 1)312 CALL lim_var_agg(2) 313 313 CALL lim_var_glo2eqv 314 314 ! 315 315 CALL lim_sbc_init ! ice surface boundary condition 316 ! 317 IF( ln_limdiaout) CALL lim_diahsb_init ! initialization for diags 316 318 ! 317 319 fr_i(:,:) = at_i(:,:) ! initialisation of sea-ice fraction -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/OPA_SRC/ZDF/zdfric.F90
r4624 r7067 162 162 & + avmv(ji,jj,jk) + avmv(ji,jj-1,jk) ) & 163 163 & + 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 165 168 avmu(ji,jj,jk) = avmu(ji,jj,jk) + avmb(jk) * umask(ji,jj,jk) 166 169 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 337 337 IF( lk_vvl ) CALL dom_vvl_sf_swp( kstp ) ! swap of vertical scale factors 338 338 ! 339 IF( ln_diahsb ) CALL dia_hsb( kstp ) ! - ML - global conservation diagnostics 340 339 341 IF( lrst_oce ) CALL rst_write( kstp ) ! write output ocean restart file 340 342 IF( ln_sto_eos ) CALL sto_rst_write( kstp ) ! write restart file for stochastic parameters … … 351 353 ENDIF 352 354 #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) 355 356 356 357 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/PISCES/P4Z/p4zopt.F90
r6943 r7067 251 251 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout) :: pe1 , pe2 , pe3 ! PAR ( R-G-B) 252 252 REAL(wp), DIMENSION(jpi,jpj,jpk), INTENT(inout), OPTIONAL :: pe0 253 REAL(wp), DIMENSION(jpi,jpj ,jpk), INTENT(out) , OPTIONAL :: pqsr100253 REAL(wp), DIMENSION(jpi,jpj) , INTENT(out) , OPTIONAL :: pqsr100 254 254 !! * local variables 255 255 INTEGER :: ji, jj, jk ! dummy loop indices -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/TRP/trctrp.F90
r6941 r7067 69 69 CALL trc_adv( kstp ) ! horizontal & vertical advection 70 70 IF( ln_zps ) THEN 71 IF( ln_isfcav ) THEN ; CALL zps_hde_isf( k t, jptra, trb, pgtu=gtru, pgtv=gtrv, pgtui=gtrui, pgtvi=gtrvi ) ! both top & bottom72 ELSE ; CALL zps_hde ( k t, jptra, trb, gtru, gtrv ) ! only bottom71 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 73 73 ENDIF 74 74 ENDIF -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcrst.F90
r5513 r7067 304 304 IF(lwp) WRITE(numout,9000) jn, TRIM( ctrcnm(jn) ), zmean, zmin, zmax, zdrift 305 305 END DO 306 WRITE(numout,*)306 IF(lwp) WRITE(numout,*) 307 307 9000 FORMAT(' tracer nb :',i2,' name :',a10,' mean :',e18.10,' min :',e18.10, & 308 308 & ' max :',e18.10,' drift :',e18.10, ' %') -
branches/2015/dev_r6946_Offline_vvl/NEMOGCM/NEMO/TOP_SRC/trcstp.F90
r6952 r7067 33 33 REAL(wp) :: rdt_sampl 34 34 INTEGER :: nb_rec_per_day 35 INTEGER :: isecfst, iseclast35 REAL(wp) :: rsecfst, rseclast 36 36 LOGICAL :: llnew 37 37 … … 131 131 INTEGER, INTENT(in) :: kt 132 132 INTEGER :: jn 133 REAL(wp) :: z secfst133 REAL(wp) :: zkt 134 134 CHARACTER(len=1) :: cl1 ! 1 character 135 135 CHARACTER(len=2) :: cl2 ! 2 characters … … 137 137 IF( kt == nittrc000 ) THEN 138 138 IF( ln_cpl ) THEN 139 rdt_sampl = 86400./ ncpl_qsr_freq139 rdt_sampl = rday / ncpl_qsr_freq 140 140 nb_rec_per_day = ncpl_qsr_freq 141 141 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 ) 144 144 ENDIF 145 145 ! … … 155 155 IF( ln_rsttr .AND. iom_varid( numrtr, 'qsr_mean' , ldstop = .FALSE. ) > 0 .AND. & 156 156 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 ' 159 161 CALL iom_get( numrtr, jpdom_autoglo, 'qsr_mean', qsr_mean ) ! A mean of qsr 160 CALL iom_get( numrtr, 'zsecfst', zsecfst ) ! A mean of qsr161 isecfst = INT( zsecfst )162 162 DO jn = 1, nb_rec_per_day 163 163 IF( jn <= 9 ) THEN … … 171 171 ELSE !* no restart: set from nit000 values 172 172 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 step173 rsecfst = kt * rdttrc(1) 174 174 ! 175 175 qsr_mean(:,:) = qsr(:,:) … … 181 181 ENDIF 182 182 ! 183 iseclast = nsec_year + nsec1jan000184 ! 185 llnew = ( iseclast - isecfst ) > INT( rdt_sampl )! new shortwave to store183 rseclast = kt * rdttrc(1) 184 ! 185 llnew = ( rseclast - rsecfst ) .ge. rdt_sampl ! new shortwave to store 186 186 IF( llnew ) THEN 187 187 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 = iseclast188 & ' time = ', rseclast/3600.,'hours ' 189 rsecfst = rseclast 190 190 DO jn = 1, nb_rec_per_day - 1 191 191 qsr_arr(:,:,jn) = qsr_arr(:,:,jn+1) … … 199 199 IF(lwp) WRITE(numout,*) 'trc_mean_qsr : write qsr_mean in restart file kt =', kt 200 200 IF(lwp) WRITE(numout,*) '~~~~~~~' 201 zkt = REAL( kt, wp ) 202 CALL iom_rstput( kt, nitrst, numrtw, 'ktdcy', zkt ) 201 203 DO jn = 1, nb_rec_per_day 202 204 IF( jn <= 9 ) THEN … … 209 211 ENDDO 210 212 CALL iom_rstput( kt, nitrst, numrtw, 'qsr_mean', qsr_mean(:,:) ) 211 zsecfst = REAL( isecfst, wp )212 CALL iom_rstput( kt, nitrst, numrtw, 'zsecfst', zsecfst )213 213 ENDIF 214 214 !
Note: See TracChangeset
for help on using the changeset viewer.
