Changeset 7355
- Timestamp:
- 2016-11-28T18:21:42+01:00 (8 years ago)
- Location:
- branches/2016/dev_CNRS_2016/NEMOGCM
- Files:
-
- 20 added
- 1 deleted
- 59 edited
- 5 copied
Legend:
- Unmodified
- Added
- Removed
-
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/ISOMIP/MY_SRC/usrdef_sbc.F90
r7164 r7355 26 26 PRIVATE 27 27 28 PUBLIC usr_def_sbc ! routine called by sbcmod.F90 28 PUBLIC usrdef_sbc_oce ! routine called by sbcmod.F90 29 PUBLIC usrdef_sbc_ice_tau ! routine called by sbcmod.F90 30 PUBLIC usrdef_sbc_ice_flx ! routine called by sbcmod.F90 29 31 30 32 !! * Substitutions … … 37 39 CONTAINS 38 40 39 SUBROUTINE usr _def_sbc( kt )41 SUBROUTINE usrdef_sbc_oce( kt ) 40 42 !!--------------------------------------------------------------------- 41 43 !! *** ROUTINE usr_def_sbc *** … … 71 73 ENDIF 72 74 ! 73 END SUBROUTINE usr_def_sbc 75 END SUBROUTINE usrdef_sbc_oce 76 77 SUBROUTINE usrdef_sbc_ice_tau( kt ) 78 INTEGER, INTENT(in) :: kt ! ocean time step 79 END SUBROUTINE usrdef_sbc_ice_tau 80 81 SUBROUTINE usrdef_sbc_ice_flx( kt ) 82 INTEGER, INTENT(in) :: kt ! ocean time step 83 END SUBROUTINE usrdef_sbc_ice_flx 74 84 75 85 !!====================================================================== -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/LOCK_EXCHANGE/MY_SRC/usrdef_sbc.F90
r7164 r7355 26 26 PRIVATE 27 27 28 PUBLIC usr_def_sbc ! routine called by sbcmod.F90 28 PUBLIC usrdef_sbc_oce ! routine called in sbcmod module 29 PUBLIC usrdef_sbc_ice_tau ! routine called by sbcice_lim.F90 for ice dynamics 30 PUBLIC usrdef_sbc_ice_flx ! routine called by sbcice_lim.F90 for ice thermo 29 31 30 32 !! * Substitutions … … 37 39 CONTAINS 38 40 39 SUBROUTINE usr _def_sbc( kt )41 SUBROUTINE usrdef_sbc_oce( kt ) 40 42 !!--------------------------------------------------------------------- 41 43 !! *** ROUTINE usr_def_sbc *** … … 71 73 ENDIF 72 74 ! 73 END SUBROUTINE usr_def_sbc 75 END SUBROUTINE usrdef_sbc_oce 76 77 SUBROUTINE usrdef_sbc_ice_tau( kt ) 78 INTEGER, INTENT(in) :: kt ! ocean time step 79 END SUBROUTINE usrdef_sbc_ice_tau 80 81 SUBROUTINE usrdef_sbc_ice_flx( kt ) 82 INTEGER, INTENT(in) :: kt ! ocean time step 83 END SUBROUTINE usrdef_sbc_ice_flx 74 84 75 85 !!====================================================================== -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/iodef.xml
r5517 r7355 61 61 <file id="file2" name_suffix="_SBC" description="surface fluxes variables" > <!-- time step automaticaly defined based on nn_fsbc --> 62 62 <field field_ref="empmr" name="wfo" /> 63 <field field_ref="emp_oce" name="emp_oce" long_name="Evap minus Precip over ocean" /> 64 <field field_ref="emp_ice" name="emp_ice" long_name="Evap minus Precip over ice" /> 63 65 <field field_ref="qsr_oce" name="qsr_oce" /> 64 66 <field field_ref="qns_oce" name="qns_oce" /> -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/ORCA2_LIM3/EXP00/namelist_ice_cfg
r4690 r7355 1 1 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2 !! NEMO/LIM-2 : Ice configuration namelist. Overwrites SHARED/namelist_ice_lim2_ref 2 !! LIM3 configuration namelist: Overwrites SHARED/namelist_ice_lim3_ref 3 !! 1 - Generic parameters (namicerun) 4 !! 2 - Diagnostics (namicediag) 5 !! 3 - Ice initialization (namiceini) 6 !! 4 - Ice discretization (namiceitd) 7 !! 5 - Ice dynamics and transport (namicedyn) 8 !! 6 - Ice diffusion (namicehdf) 9 !! 7 - Ice thermodynamics (namicethd) 10 !! 8 - Ice salinity (namicesal) 11 !! 9 - Ice mechanical redistribution (namiceitdme) 3 12 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 4 5 !----------------------------------------------------------------------- 6 &namicerun ! Share parameters for dynamics/advection/thermo 7 !----------------------------------------------------------------------- 13 !------------------------------------------------------------------------------ 14 &namicerun ! Generic parameters 15 !------------------------------------------------------------------------------ 8 16 / 9 !----------------------------------------------------------------------- 10 &namice ini ! ice initialisation11 !----------------------------------------------------------------------- 17 !------------------------------------------------------------------------------ 18 &namicediag ! Diagnostics 19 !------------------------------------------------------------------------------ 12 20 / 13 !----------------------------------------------------------------------- 14 &namice dyn ! ice dynamic15 !----------------------------------------------------------------------- 21 !------------------------------------------------------------------------------ 22 &namiceini ! Ice initialization 23 !------------------------------------------------------------------------------ 16 24 / 17 !----------------------------------------------------------------------- 18 &namice thd ! ice thermodynamic19 !----------------------------------------------------------------------- 25 !------------------------------------------------------------------------------ 26 &namiceitd ! Ice discretization 27 !------------------------------------------------------------------------------ 20 28 / 21 !----------------------------------------------------------------------- 22 &namice sal ! ice salinity23 !----------------------------------------------------------------------- 29 !------------------------------------------------------------------------------ 30 &namicedyn ! Ice dynamics and transport 31 !------------------------------------------------------------------------------ 24 32 / 25 !----------------------------------------------------------------------- 26 &namice itdme ! parameters for mechanical redistribution of ice27 !----------------------------------------------------------------------- 33 !------------------------------------------------------------------------------ 34 &namicehdf ! Ice horizontal diffusion 35 !------------------------------------------------------------------------------ 28 36 / 29 !----------------------------------------------------------------------- 30 &namice dia ! ice diagnostics31 !----------------------------------------------------------------------- 37 !------------------------------------------------------------------------------ 38 &namicethd ! Ice thermodynamics 39 !------------------------------------------------------------------------------ 32 40 / 41 !------------------------------------------------------------------------------ 42 &namicesal ! Ice salinity 43 !------------------------------------------------------------------------------ 44 / 45 !------------------------------------------------------------------------------ 46 &namiceitdme ! Ice mechanical redistribution (ridging and rafting) 47 !------------------------------------------------------------------------------ 48 / -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/OVERFLOW/MY_SRC/usrdef_sbc.F90
r7164 r7355 26 26 PRIVATE 27 27 28 PUBLIC usr_def_sbc ! routine called by sbcmod.F90 28 PUBLIC usrdef_sbc_oce ! routine called in sbcmod module 29 PUBLIC usrdef_sbc_ice_tau ! routine called by sbcice_lim.F90 for ice dynamics 30 PUBLIC usrdef_sbc_ice_flx ! routine called by sbcice_lim.F90 for ice thermo 29 31 30 32 !! * Substitutions … … 37 39 CONTAINS 38 40 39 SUBROUTINE usr _def_sbc( kt )41 SUBROUTINE usrdef_sbc_oce( kt ) 40 42 !!--------------------------------------------------------------------- 41 43 !! *** ROUTINE usr_def_sbc *** … … 71 73 ENDIF 72 74 ! 73 END SUBROUTINE usr_def_sbc 75 END SUBROUTINE usrdef_sbc_oce 76 77 SUBROUTINE usrdef_sbc_ice_tau( kt ) 78 INTEGER, INTENT(in) :: kt ! ocean time step 79 END SUBROUTINE usrdef_sbc_ice_tau 80 81 SUBROUTINE usrdef_sbc_ice_flx( kt ) 82 INTEGER, INTENT(in) :: kt ! ocean time step 83 END SUBROUTINE usrdef_sbc_ice_flx 74 84 75 85 !!====================================================================== -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/SHARED/field_def.xml
r7278 r7355 239 239 <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" /> 240 240 <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" /> 241 <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" /> 241 <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" /> 242 <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" /> 243 <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" /> 242 244 <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" /> 243 245 <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" /> … … 317 319 <field id="icevolu" long_name="ice volume" unit="m" /> 318 320 <field id="snowvol" long_name="snow volume" unit="m" /> 321 <field id="tau_icebfr" long_name="ice friction on ocean bottom for landfast ice" unit="N/m2" /> 319 322 320 323 <field id="icetrp" long_name="ice volume transport" unit="m/day" /> … … 330 333 <field id="sfxbom" long_name="salt flux from bot melt" unit="1e-3*kg/m2/day" /> 331 334 <field id="sfxsum" long_name="salt flux from surf melt" unit="1e-3*kg/m2/day" /> 335 <field id="sfxlam" long_name="salt flux from lateral melt" unit="1e-3*kg/m2/day" /> 332 336 <field id="sfxsni" long_name="salt flux from snow-ice formation" unit="1e-3*kg/m2/day" /> 333 337 <field id="sfxopw" long_name="salt flux from open water ice formation" unit="1e-3*kg/m2/day" /> … … 340 344 <field id="vfxsni" long_name="daily snowice ice prod." unit="m/day" /> 341 345 <field id="vfxsum" long_name="surface melt" unit="m/day" /> 346 <field id="vfxlam" long_name="lateral melt" unit="m/day" /> 342 347 <field id="vfxbom" long_name="bottom melt" unit="m/day" /> 343 348 <field id="vfxres" long_name="daily resultant ice prod./melting from limupdate" unit="m/day" /> … … 345 350 <field id="vfxsnw" long_name="snw melt/growth" unit="m/day" /> 346 351 <field id="vfxsub" long_name="snw sublimation" unit="m/day" /> 352 <field id="vfxsub_err" long_name="excess of snw sublimation sent to ocean" unit="m/day" /> 347 353 <field id="vfxspr" long_name="snw precipitation on ice" unit="m/day" /> 348 <field id="vfxthin" long_name="daily thermo ice prod. for thin ice( <20cm) + open water"unit="m/day" />354 <field id="vfxthin" long_name="daily thermo ice prod. for thin ice(20cm) + open water" unit="m/day" /> 349 355 350 356 <field id="afxtot" long_name="area tendency (total)" unit="day-1" /> … … 521 527 522 528 <!-- available with ln_diahsb --> 523 <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" />529 <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" /> 524 530 <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" /> 525 <field id="bgheatco" long_name="drift in global mean heat content wrt timestep 1" unit="10^9J" /> 526 <field id="bgsaltco" long_name="drift in global mean salt content wrt timestep 1" unit="1e-3*m3" /> 531 <field id="bgheatco" long_name="drift in global mean heat content wrt timestep 1" unit="1.e20J" /> 532 <field id="bgheatfx" long_name="drift in global mean heat flux wrt timestep 1" unit="W/m2" /> 533 <field id="bgsaltco" long_name="drift in global mean salt content wrt timestep 1" unit="1e-3*km3" /> 527 534 <field id="bgvolssh" long_name="drift in global mean ssh volume wrt timestep 1" unit="km3" /> 528 535 <field id="bgvole3t" long_name="drift in global mean volume variation (e3t) wrt timestep 1" unit="km3" /> 529 <field id="bgvoltot" long_name="drift in global mean volume wrt timestep 1" unit="km3" /> 530 <!-- NOTE: No matching iom_put call --> 531 <field id="bgsshtot" long_name="drift in global mean ssh wrt timestep 1" standard_name="global_average_sea_level_change" unit="m" /> 532 <field id="bgfrcvol" long_name="drift in global mean volume from forcing wrt timestep 1" unit="km3" /> 533 <field id="bgfrctem" long_name="drift in global mean heat content from forcing wrt timestep 1" unit="10^9J" /> 534 <field id="bgfrcsal" long_name="drift in global mean salt content from forcing wrt timestep 1" unit="1e-3*km3" /> 535 <field id="bgmistem" long_name="global mean temperature error due to free surface" unit="degC" /> 536 <field id="bgmissal" long_name="global mean salinity error due to free surface" unit="1e-3" /> 537 </field_group> 536 <field id="bgfrcvol" long_name="global mean volume from forcing" unit="km3" /> 537 <field id="bgfrctem" long_name="global mean heat content from forcing" unit="1.e20J" /> 538 <field id="bgfrchfx" long_name="global mean heat flux from forcing" unit="W/m2" /> 539 <field id="bgfrcsal" long_name="global mean salt content from forcing" unit="1e-3*km3" /> 540 <field id="bgmistem" long_name="global mean temperature error due to free surface (no vvl)" unit="degC" /> 541 <field id="bgmissal" long_name="global mean salinity error due to free surface (no vvl)" unit="1e-3" /> 542 </field_group> 538 543 539 544 <!-- LIM3 scalar variables --> … … 541 546 <field_group id="SBC_scalar" domain_ref="1point" > 542 547 <!-- available with ln_limdiaout --> 543 <field id="ibgvoltot" long_name="global mean ice volume" unit="km3" /> 544 <field id="sbgvoltot" long_name="global mean snow volume" unit="km3" /> 545 <field id="ibgarea" long_name="global mean ice area" unit="km2" /> 546 <field id="ibgsaline" long_name="global mean ice salinity" unit="1e-3" /> 547 <field id="ibgtemper" long_name="global mean ice temperature" unit="degC" /> 548 <field id="ibgheatco" long_name="global mean ice heat content" unit="10^20J" /> 549 <field id="sbgheatco" long_name="global mean snow heat content" unit="10^20J" /> 550 <field id="ibgsaltco" long_name="global mean ice salt content" unit="1e-3*km3" /> 551 552 <field id="ibgvfx" long_name="global mean volume flux (emp)" unit="m/day" /> 553 <field id="ibgvfxbog" long_name="global mean volume flux (bottom growth)" unit="m/day" /> 554 <field id="ibgvfxopw" long_name="global mean volume flux (open water growth)" unit="m/day" /> 555 <field id="ibgvfxsni" long_name="global mean volume flux (snow-ice growth)" unit="m/day" /> 556 <field id="ibgvfxdyn" long_name="global mean volume flux (dynamic growth)" unit="m/day" /> 557 <field id="ibgvfxbom" long_name="global mean volume flux (bottom melt)" unit="m/day" /> 558 <field id="ibgvfxsum" long_name="global mean volume flux (surface melt)" unit="m/day" /> 559 <field id="ibgvfxres" long_name="global mean volume flux (resultant)" unit="m/day" /> 560 <field id="ibgvfxspr" long_name="global mean volume flux (snow precip)" unit="m/day" /> 561 <field id="ibgvfxsnw" long_name="global mean volume flux (snow melt)" unit="m/day" /> 562 <field id="ibgvfxsub" long_name="global mean volume flux (snow sublimation)" unit="m/day" /> 563 564 <field id="ibgsfx" long_name="global mean salt flux (total)" unit="1e-3*m/day" /> 565 <field id="ibgsfxbri" long_name="global mean salt flux (brines)" unit="1e-3*m/day" /> 566 <field id="ibgsfxdyn" long_name="global mean salt flux (dynamic)" unit="1e-3*m/day" /> 567 <field id="ibgsfxres" long_name="global mean salt flux (resultant)" unit="1e-3*m/day" /> 568 <field id="ibgsfxbog" long_name="global mean salt flux (thermo)" unit="1e-3*m/day" /> 569 <field id="ibgsfxopw" long_name="global mean salt flux (thermo)" unit="1e-3*m/day" /> 570 <field id="ibgsfxsni" long_name="global mean salt flux (thermo)" unit="1e-3*m/day" /> 571 <field id="ibgsfxbom" long_name="global mean salt flux (thermo)" unit="1e-3*m/day" /> 572 <field id="ibgsfxsum" long_name="global mean salt flux (thermo)" unit="1e-3*m/day" /> 573 <field id="ibgsfxsub" long_name="global mean salt flux (thermo)" unit="1e-3*m/day" /> 574 575 <field id="ibghfxdhc" long_name="Heat content variation in snow and ice" unit="W" /> 576 <field id="ibghfxspr" long_name="Heat content of snow precip" unit="W" /> 577 578 <field id="ibghfxthd" long_name="heat fluxes from ice-ocean exchange during thermo" unit="W" /> 579 <field id="ibghfxsum" long_name="heat fluxes causing surface ice melt" unit="W" /> 580 <field id="ibghfxbom" long_name="heat fluxes causing bottom ice melt" unit="W" /> 581 <field id="ibghfxbog" long_name="heat fluxes causing bottom ice growth" unit="W" /> 582 <field id="ibghfxdif" long_name="heat fluxes causing ice temperature change" unit="W" /> 583 <field id="ibghfxopw" long_name="heat fluxes causing open water ice formation" unit="W" /> 584 <field id="ibghfxdyn" long_name="heat fluxes from ice-ocean exchange during dynamic" unit="W" /> 585 <field id="ibghfxres" long_name="heat fluxes from ice-ocean exchange during resultant" unit="W" /> 586 <field id="ibghfxsub" long_name="heat fluxes from sublimation" unit="W" /> 587 <field id="ibghfxsnw" long_name="heat fluxes from snow-ocean exchange" unit="W" /> 588 <field id="ibghfxout" long_name="non solar heat fluxes received by the ocean" unit="W" /> 589 <field id="ibghfxin" long_name="total heat fluxes at the ice surface" unit="W" /> 590 591 <field id="ibgfrcvol" long_name="global mean forcing volume (emp)" unit="km3" /> 592 <field id="ibgfrcsfx" long_name="global mean forcing salt (sfx)" unit="1e-3*km3" /> 593 <field id="ibgvolgrm" long_name="global mean ice growth+melt volume" unit="km3" /> 548 <field id="ibgfrcvoltop" long_name="global mean ice/snow forcing at interface ice/snow-atm (volume equivalent ocean volume)" unit="km3" /> 549 <field id="ibgfrcvolbot" long_name="global mean ice/snow forcing at interface ice/snow-ocean (volume equivalent ocean volume)" unit="km3" /> 550 <field id="ibgfrctemtop" long_name="global mean heat on top of ice/snw/ocean-atm " unit="1e20J" /> 551 <field id="ibgfrctembot" long_name="global mean heat below ice (on top of ocean) " unit="1e20J" /> 552 <field id="ibgfrcsal" long_name="global mean ice/snow forcing (salt equivalent ocean volume)" unit="pss*km3" /> 553 <field id="ibgfrchfxtop" long_name="global mean heat flux on top of ice/snw/ocean-atm " unit="W/m2" /> 554 <field id="ibgfrchfxbot" long_name="global mean heat flux below ice (on top of ocean) " unit="W/m2" /> 555 556 <field id="ibgvolume" long_name="drift in ice/snow volume (equivalent ocean volume)" unit="km3" /> 557 <field id="ibgsaltco" long_name="drift in ice salt content (equivalent ocean volume)" unit="pss*km3" /> 558 <field id="ibgheatco" long_name="drift in ice/snow heat content" unit="1e20J" /> 559 <field id="ibgheatfx" long_name="drift in ice/snow heat flux" unit="W/m2" /> 560 561 <field id="ibgvol_tot" long_name="global mean ice volume" unit="km3" /> 562 <field id="sbgvol_tot" long_name="global mean snow volume" unit="km3" /> 563 <field id="ibgarea_tot" long_name="global mean ice area" unit="km2" /> 564 <field id="ibgsalt_tot" long_name="global mean ice salt content" unit="1e-3*km3" /> 565 <field id="ibgheat_tot" long_name="global mean ice heat content" unit="1e20J" /> 566 <field id="sbgheat_tot" long_name="global mean snow heat content" unit="1e20J" /> 594 567 </field_group> 595 568 -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/SHARED/namelist_ice_lim3_ref
r7278 r7355 1 1 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 2 !! LIM3 namelist 2 !! LIM3 namelist: 3 3 !! 1 - Generic parameters (namicerun) 4 !! 2 - Ice initialization (namiceini) 5 !! 3 - Ice discretization (namiceitd) 6 !! 4 - Ice dynamics and transport (namicedyn) 7 !! 5 - Ice thermodynamics (namicethd) 8 !! 6 - Ice salinity (namicesal) 9 !! 7 - Ice mechanical redistribution (namiceitdme) 10 !>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 4 !! 2 - Diagnostics (namicediag) 5 !! 3 - Ice initialization (namiceini) 6 !! 4 - Ice discretization (namiceitd) 7 !! 5 - Ice dynamics and transport (namicedyn) 8 !! 6 - Ice diffusion (namicehdf) 9 !! 7 - Ice thermodynamics (namicethd) 10 !! 8 - Ice salinity (namicesal) 11 !! 9 - Ice mechanical redistribution (namiceitdme) 12 !!>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> 11 13 ! 12 14 !------------------------------------------------------------------------------ 13 15 &namicerun ! Generic parameters 14 16 !------------------------------------------------------------------------------ 15 jpl = 5 ! number of ice categories 16 nlay_i = 2 ! number of ice layers 17 nlay_s = 1 ! number of snow layers (only 1 is working) 18 cn_icerst_in = "restart_ice" ! suffix of ice restart name (input) 19 cn_icerst_indir = "." ! directory from which to read input ice restarts 20 cn_icerst_out = "restart_ice" ! suffix of ice restart name (output) 21 cn_icerst_outdir = "." ! directory in which to write output ice restarts 22 ln_limdyn = .true. ! ice dynamics (T) or thermodynamics only (F) 23 rn_amax_n = 0.999 ! maximum tolerated ice concentration NH 24 rn_amax_s = 0.999 ! maximum tolerated ice concentration SH 25 ln_limdiahsb = .false. ! check the heat and salt budgets (T) or not (F) 26 ln_limdiaout = .true. ! output the heat and salt budgets (T) or not (F) 27 ln_icectl = .false. ! ice points output for debug (T or F) 28 iiceprt = 10 ! i-index for debug 29 jiceprt = 10 ! j-index for debug 17 jpl = 5 ! number of ice categories 18 nlay_i = 2 ! number of ice layers 19 nlay_s = 1 ! number of snow layers (only 1 is working) 20 rn_amax_n = 0.997 ! maximum tolerated ice concentration NH 21 rn_amax_s = 0.997 ! maximum tolerated ice concentration SH 22 cn_icerst_in = "restart_ice" ! suffix of ice restart name (input) 23 cn_icerst_out = "restart_ice" ! suffix of ice restart name (output) 24 cn_icerst_indir = "." ! directory to read input ice restarts 25 cn_icerst_outdir = "." ! directory to write output ice restarts 26 ln_limthd = .true. ! ice thermo (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 27 ln_limdyn = .true. ! ice dynamics (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 28 nn_limdyn = 2 ! (ln_limdyn=T) switch for ice dynamics 29 ! 2: total 30 ! 1: advection only (no diffusion, no ridging/rafting) 31 ! 0: advection only (as 1 but with prescribed velocity, bypass rheology) 32 rn_uice = 0.00001 ! (nn_limdyn=0) ice u-velocity 33 rn_vice = -0.00001 ! (nn_limdyn=0) ice v-velocity 34 / 35 !------------------------------------------------------------------------------ 36 &namicediag ! Diagnostics 37 !------------------------------------------------------------------------------ 38 ln_limdiachk = .false. ! check online the heat, mass & salt budgets (T) or not (F) 39 ln_limdiahsb = .false. ! output the heat, mass & salt budgets (T) or not (F) 40 ln_limctl = .false. ! ice points output for debug (T or F) 41 iiceprt = 10 ! i-index for debug 42 jiceprt = 10 ! j-index for debug 30 43 / 31 44 !------------------------------------------------------------------------------ 32 45 &namiceini ! Ice initialization 33 46 !------------------------------------------------------------------------------ 34 ln_iceini = .true. ! activate ice initialization (T) or not (F) 35 rn_thres_sst = 2.0 ! maximum water temperature with initial ice (degC) 36 rn_hts_ini_n = 0.3 ! initial real snow thickness (m), North 37 rn_hts_ini_s = 0.3 ! " " South 38 rn_hti_ini_n = 3.0 ! initial real ice thickness (m), North 39 rn_hti_ini_s = 1.0 ! " " South 40 rn_ati_ini_n = 0.9 ! initial ice concentration (-), North 41 rn_ati_ini_s = 0.9 ! " " South 42 rn_smi_ini_n = 6.3 ! initial ice salinity (g/kg), North 43 rn_smi_ini_s = 6.3 ! " " South 44 rn_tmi_ini_n = 270. ! initial ice/snw temperature (K), North 45 rn_tmi_ini_s = 270. ! " " South 47 ! -- limistate -- ! 48 ln_limini = .true. ! activate ice initialization (T) or not (F) 49 ln_limini_file = .false. ! netcdf file provided for initialization (T) or not (F) 50 rn_thres_sst = 2.0 ! maximum water temperature with initial ice (degC) 51 rn_hts_ini_n = 0.3 ! initial real snow thickness (m), North 52 rn_hts_ini_s = 0.3 ! " " South 53 rn_hti_ini_n = 3.0 ! initial real ice thickness (m), North 54 rn_hti_ini_s = 1.0 ! " " South 55 rn_ati_ini_n = 0.9 ! initial ice concentration (-), North 56 rn_ati_ini_s = 0.9 ! " " South 57 rn_smi_ini_n = 6.3 ! initial ice salinity (g/kg), North 58 rn_smi_ini_s = 6.3 ! " " South 59 rn_tmi_ini_n = 270. ! initial ice/snw temperature (K), North 60 rn_tmi_ini_s = 270. ! " " South 46 61 / 47 62 !------------------------------------------------------------------------------ … … 56 71 &namicedyn ! Ice dynamics and transport 57 72 !------------------------------------------------------------------------------ 58 nn_icestr = 0 ! ice strength parameteriztaion 59 ! 0: Hibler_79 P = pstar*<h>*exp(-c_rhg*A) 60 ! 1: Rothrock_75 P = Cf*coeff*integral(wr.h^2) 61 ln_icestr_bvf = .false. ! ice strength function brine volume (T) or not (F) 62 rn_pe_rdg = 17.0 ! ridging work divided by pot. energy change in ridging, if nn_icestr = 1 63 rn_pstar = 2.0e+04 ! ice strength thickness parameter (N/m2), nn_icestr = 0 64 rn_crhg = 20.0 ! ice strength conc. parameter (-), nn_icestr = 0 65 rn_cio = 5.0e-03 ! ice-ocean drag coefficient (-) 66 rn_creepl = 1.0e-12 ! creep limit (s-1) 67 rn_ecc = 2.0 ! eccentricity of the elliptical yield curve 68 nn_nevp = 120 ! number of EVP subcycles 69 rn_relast = 0.333 ! ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 70 ! advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 71 nn_ahi0 = 2 ! horizontal diffusivity computation 72 ! 0: use rn_ahi0_ref 73 ! 1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 74 ! 2: use rn_ahi0_ref x grid cell length / ( 2deg mean grid cell length ) 75 rn_ahi0_ref = 350.0 ! horizontal sea ice diffusivity (m2/s) 76 ! if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 73 ! -- limtrp & limadv -- ! 74 nn_limadv = 0 ! choose the advection scheme (-1=Prather ; 0=Ultimate-Macho) 75 nn_limadv_ord = 5 ! choose the order of the advection scheme (if nn_limadv=0) 76 ! -- limitd_me -- ! 77 nn_icestr = 0 ! ice strength parameteriztaion 78 ! 0: Hibler_79 P = pstar*<h>*exp(-c_rhg*A) 79 ! 1: Rothrock_75 P = Cf*coeff*integral(wr.h^2) 80 rn_pe_rdg = 17.0 ! (nn_icestr=1) ridging work divided by pot. energy change in ridging 81 rn_pstar = 2.0e+04 ! (nn_icestr=0) ice strength thickness parameter (N/m2) 82 rn_crhg = 20.0 ! (nn_icestr=0) ice strength conc. parameter (-) 83 ln_icestr_bvf = .false. ! ice strength function brine volume (T) or not (F) 84 ! 85 ! -- limdyn & limrhg -- ! 86 rn_cio = 5.0e-03 ! ice-ocean drag coefficient (-) 87 rn_creepl = 1.0e-12 ! creep limit (s-1) 88 rn_ecc = 2.0 ! eccentricity of the elliptical yield curve 89 nn_nevp = 120 ! number of EVP subcycles 90 rn_relast = 0.333 ! ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 91 ! advised value: 1/3 (rn_nevp=120) or 1/9 (rn_nevp=300) 92 ln_landfast = .false. ! landfast ice parameterization (T or F) 93 rn_gamma = 0.15 ! (ln_landfast=T) fraction of ocean depth that ice must reach to initiate landfast 94 ! recommended range: [0.1 ; 0.25] 95 rn_icebfr = 10. ! (ln_landfast=T) maximum bottom stress per unit area of contact (N/m2) 96 ! a very large value ensures ice velocity=0 even with a small contact area 97 ! recommended range: ?? (should be greater than atm-ice stress => >0.1 N/m2) 98 rn_lfrelax = 1.e-5 ! (ln_landfast=T) relaxation time scale to reach static friction (s-1) 77 99 / 78 100 !------------------------------------------------------------------------------ 79 101 &namicehdf ! Ice horizontal diffusion 80 102 !------------------------------------------------------------------------------ 81 nn_convfrq = 5 ! convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 103 ! -- limhdf -- ! 104 nn_ahi0 = -1 ! horizontal diffusivity computation 105 ! -1: no diffusion (bypass limhdf) 106 ! 0: use rn_ahi0_ref 107 ! 1: use rn_ahi0_ref x mean grid cell length / ( 2deg mean grid cell length ) 108 ! 2: use rn_ahi0_ref x grid cell length / ( 2deg mean grid cell length ) 109 rn_ahi0_ref = 350.0 ! horizontal sea ice diffusivity (m2/s) 110 ! if nn_ahi0 > 0, rn_ahi0_ref is the reference value at a nominal 2 deg resolution 82 111 / 83 112 !------------------------------------------------------------------------------ 84 113 &namicethd ! Ice thermodynamics 85 114 !------------------------------------------------------------------------------ 86 rn_hnewice = 0.1 ! thickness for new ice formation in open water (m) 87 ln_frazil = .false. ! use frazil ice collection thickness as a function of wind (T) or not (F) 88 rn_maxfrazb = 1.0 ! maximum fraction of frazil ice collecting at the ice base 89 rn_vfrazb = 0.417 ! thresold drift speed for frazil ice collecting at the ice bottom (m/s) 90 rn_Cfrazb = 5.0 ! squeezing coefficient for frazil ice collecting at the ice bottom 91 rn_himin = 0.10 ! minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice 92 rn_betas = 0.66 ! exponent in lead-ice repratition of snow precipitation 93 ! betas = 1 -> equipartition, betas < 1 -> more on leads 94 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice (m-1) 95 nn_conv_dif = 50 ! maximal number of iterations for heat diffusion computation 96 rn_terr_dif = 0.0001 ! maximum temperature after heat diffusion (degC) 97 nn_ice_thcon= 1 ! sea ice thermal conductivity 98 ! 0: k = k0 + beta.S/T (Untersteiner, 1964) 99 ! 1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 100 nn_monocat = 0 ! virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0) 101 ! 2: simple piling instead of ridging --- temporary option 102 ! 3: activate G(he) only --- temporary option 103 ! 4: activate lateral melting only --- temporary option 104 ln_it_qnsice = .true. ! iterate the surface non-solar flux with surface temperature (T) or not (F) 115 ! -- limthd_dif -- ! 116 rn_kappa_i = 1.0 ! radiation attenuation coefficient in sea ice (m-1) 117 nn_conv_dif = 50 ! maximal number of iterations for heat diffusion computation 118 rn_terr_dif = 1.0e-04 ! maximum temperature after heat diffusion (degC) 119 nn_ice_thcon = 1 ! sea ice thermal conductivity 120 ! 0: k = k0 + beta.S/T (Untersteiner, 1964) 121 ! 1: k = k0 + beta1.S/T - beta2.T (Pringle et al., 2007) 122 ln_it_qnsice = .true. ! iterate the surface non-solar flux with surface temperature (T) or not (F) 123 nn_monocat = 0 ! virtual ITD mono-category parameterizations (1, jpl = 1 only) or not (0) 124 ! 2: simple piling instead of ridging --- temporary option 125 ! 3: activate G(he) only --- temporary option 126 ! 4: activate extra lateral melting only --- temporary option 127 ! -- limthd_dh -- ! 128 ln_limdH = .true. ! activate ice thickness change from growing/melting (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 129 rn_betas = 0.66 ! exponent in lead-ice repratition of snow precipitation 130 ! betas = 1 -> equipartition, betas < 1 -> more on leads 131 ! -- limthd_da -- ! 132 ln_limdA = .true. ! activate lateral melting param. (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 133 rn_beta = 1.0 ! (ln_latmelt=T) coef. beta for lateral melting param. Recommended range=[0.8-1.2] 134 ! => decrease = more melt and melt peaks toward higher concentration (A~0.5 for beta=1 ; A~0.8 for beta=0.2) 135 ! 0.3 = best fit for western Fram Strait and Antarctica 136 ! 1.4 = best fit for eastern Fram Strait 137 rn_dmin = 8. ! (ln_latmelt=T) minimum floe diameter for lateral melting param. Recommended range=[6-10] 138 ! => 6 vs 8m = +40% melting at the peak (A~0.5) 139 ! 10 vs 8m = -20% melting 140 ! -- limthd_lac -- ! 141 ln_limdO = .true. ! activate ice growth in open-water (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 142 rn_hnewice = 0.1 ! thickness for new ice formation in open water (m) 143 ln_frazil = .false. ! Frazil ice parameterization (ice collection as a function of wind) 144 rn_maxfrazb = 1.0 ! (ln_frazil=T) maximum fraction of frazil ice collecting at the ice base 145 rn_vfrazb = 0.417 ! (ln_frazil=T) thresold drift speed for frazil ice collecting at the ice bottom (m/s) 146 rn_Cfrazb = 5.0 ! (ln_frazil=T) squeezing coefficient for frazil ice collecting at the ice bottom 147 ! -- limitd_th -- ! 148 rn_himin = 0.1 ! minimum ice thickness (m) used in remapping, must be smaller than rn_hnewice 105 149 / 106 150 !------------------------------------------------------------------------------ 107 151 &namicesal ! Ice salinity 108 152 !------------------------------------------------------------------------------ 109 nn_icesal = 2 ! ice salinity option 110 ! 1: constant ice salinity (S=rn_icesal) 111 ! 2: varying salinity parameterization S(z,t) 112 ! 3: prescribed salinity profile S(z), Schwarzacher, 1959 113 rn_icesal = 4. ! ice salinity (g/kg, nn_icesal = 1 only) 114 rn_sal_gd = 5. ! restoring ice salinity, gravity drainage (g/kg) 115 rn_time_gd = 1.73e+6 ! restoring time scale, gravity drainage (s) 116 rn_sal_fl = 2. ! restoring ice salinity, flushing (g/kg) 117 rn_time_fl = 8.64e+5 ! restoring time scale, flushing (s) 118 rn_simax = 20. ! maximum tolerated ice salinity (g/kg) 119 rn_simin = 0.1 ! minimum tolerated ice salinity (g/kg) 153 ! -- limthd_sal -- ! 154 ln_limdS = .true. ! activate gravity drainage and flushing (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 155 nn_icesal = 2 ! ice salinity option 156 ! 1: constant ice salinity (S=rn_icesal) 157 ! 2: varying salinity parameterization S(z,t) 158 ! 3: prescribed salinity profile S(z), Schwarzacher, 1959 159 rn_icesal = 4. ! (nn_icesal=1) ice salinity (g/kg) 160 rn_sal_gd = 5. ! restoring ice salinity, gravity drainage (g/kg) 161 rn_time_gd = 1.73e+6 ! restoring time scale, gravity drainage (s) 162 rn_sal_fl = 2. ! restoring ice salinity, flushing (g/kg) 163 rn_time_fl = 8.64e+5 ! restoring time scale, flushing (s) 164 rn_simax = 20. ! maximum tolerated ice salinity (g/kg) 165 rn_simin = 0.1 ! minimum tolerated ice salinity (g/kg) 120 166 / 121 167 !------------------------------------------------------------------------------ 122 168 &namiceitdme ! Ice mechanical redistribution (ridging and rafting) 123 169 !------------------------------------------------------------------------------ 124 rn_Cs = 0.5 ! fraction of shearing energy contributing to ridging 125 rn_fsnowrdg = 0.5 ! snow volume fraction that survives in ridging 126 rn_fsnowrft = 0.5 ! snow volume fraction that survives in rafting 127 nn_partfun = 1 ! type of ridging participation function 128 ! 0: linear (Thorndike et al, 1975) 129 ! 1: exponential (Lipscomb, 2007 130 rn_gstar = 0.15 ! fractional area of thin ice being ridged (nn_partfun = 0) 131 rn_astar = 0.05 ! exponential measure of ridging ice fraction (nn_partfun = 1) 132 rn_hstar = 100.0 ! determines the maximum thickness of ridged ice (m) (Hibler, 1980) 133 ln_rafting = .true. ! rafting activated (T) or not (F) 134 rn_hraft = 0.75 ! threshold thickness for rafting (m) 135 rn_craft = 5.0 ! squeezing coefficient used in the rafting function 136 rn_por_rdg = 0.3 ! porosity of newly ridged ice (Lepparanta et al., 1995) 170 ! -- limitd_me -- ! 171 rn_cs = 0.5 ! fraction of shearing energy contributing to ridging 172 nn_partfun = 1 ! type of ridging participation function 173 ! 0: linear (Thorndike et al, 1975) 174 ! 1: exponential (Lipscomb, 2007) 175 rn_gstar = 0.15 ! (nn_partfun = 0) fractional area of thin ice being ridged 176 rn_astar = 0.05 ! (nn_partfun = 1) exponential measure of ridging ice fraction 177 ln_ridging = .true. ! ridging activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 178 rn_hstar = 100.0 ! (ln_ridging = T) determines the maximum thickness of ridged ice (m) (Hibler, 1980) 179 rn_por_rdg = 0.3 ! (ln_ridging = T) porosity of newly ridged ice (Lepparanta et al., 1995) 180 rn_fsnowrdg = 0.5 ! (ln_ridging = T) snow volume fraction that survives in ridging 181 ln_rafting = .true. ! rafting activated (T) or not (F) => DO NOT TOUCH UNLESS U KNOW WHAT U DO 182 rn_hraft = 0.75 ! (ln_rafting = T) threshold thickness for rafting (m) 183 rn_craft = 5.0 ! (ln_rafting = T) squeezing coefficient used in the rafting function 184 rn_fsnowrft = 0.5 ! (ln_rafting = T) snow volume fraction that survives in rafting 137 185 / -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/SHARED/namelist_ref
r7280 r7355 11 11 !! 6 - Tracer (nameos, namtra_adv, namtra_ldf, namtra_ldfeiv, namtra_dmp) 12 12 !! 7 - dynamics (namdyn_adv, namdyn_vor, namdyn_hpg, namdyn_spg, namdyn_ldf) 13 !! 8 - Verical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_ddm, namzdf_tmx )13 !! 8 - Verical physics (namzdf, namzdf_ric, namzdf_tke, namzdf_ddm, namzdf_tmx, namzdf_tmx_new) 14 14 !! 9 - diagnostics (namnc4, namtrd, namspr, namflo, namhsb, namsto) 15 15 !! 10 - miscellaneous (nammpp, namctl) … … 219 219 / 220 220 !----------------------------------------------------------------------- 221 &namsbc_flx ! surface boundary condition : flux formulation (ln_ana = T)221 &namsbc_flx ! surface boundary condition : flux formulation 222 222 !----------------------------------------------------------------------- 223 223 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! … … 260 260 rn_vfac = 0. ! multiplicative factor for ocean/ice velocity 261 261 ! in the calculation of the wind stress (0.=absolute winds or 1.=relative winds) 262 ln_Cd_L12 = .false. ! Modify the drag ice-atm and oce-atm depending on ice concentration 263 ! This parameterization is from Lupkes et al. (JGR 2012) 262 264 / 263 265 !----------------------------------------------------------------------- … … 291 293 &namsbc_sas ! Stand Alone Surface boundary condition 292 294 !----------------------------------------------------------------------- 293 ! ! file name ! frequency (hours) ! variable ! time interp. ! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 294 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 295 sn_usp = 'sas_grid_U' , 120 , 'vozocrtx' , .true. , .true. , 'yearly' , '' , '' , '' 296 sn_vsp = 'sas_grid_V' , 120 , 'vomecrty' , .true. , .true. , 'yearly' , '' , '' , '' 297 sn_tem = 'sas_grid_T' , 120 , 'sosstsst' , .true. , .true. , 'yearly' , '' , '' , '' 298 sn_sal = 'sas_grid_T' , 120 , 'sosaline' , .true. , .true. , 'yearly' , '' , '' , '' 299 sn_ssh = 'sas_grid_T' , 120 , 'sossheig' , .true. , .true. , 'yearly' , '' , '' , '' 300 sn_e3t = 'sas_grid_T' , 120 , 'e3t_m' , .true. , .true. , 'yearly' , '' , '' , '' 301 sn_frq = 'sas_grid_T' , 120 , 'frq_m' , .true. , .true. , 'yearly' , '' , '' , '' 295 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 296 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 297 l_sasread = .TRUE. ! Read fields in a file if .TRUE. , or initialize to 0. in sbcssm.F90 if .FALSE. 298 sn_usp = 'sas_grid_U', 120 , 'vozocrtx', .true. , .true. , 'yearly' , '' , '' , '' 299 sn_vsp = 'sas_grid_V', 120 , 'vomecrty', .true. , .true. , 'yearly' , '' , '' , '' 300 sn_tem = 'sas_grid_T', 120 , 'sosstsst', .true. , .true. , 'yearly' , '' , '' , '' 301 sn_sal = 'sas_grid_T', 120 , 'sosaline', .true. , .true. , 'yearly' , '' , '' , '' 302 sn_ssh = 'sas_grid_T', 120 , 'sossheig', .true. , .true. , 'yearly' , '' , '' , '' 303 sn_e3t = 'sas_grid_T', 120 , 'e3t_m' , .true. , .true. , 'yearly' , '' , '' , '' 304 sn_frq = 'sas_grid_T', 120 , 'frq_m' , .true. , .true. , 'yearly' , '' , '' , '' 302 305 303 306 ln_3d_uve = .true. ! specify whether we are supplying a 3D u,v and e3 field -
branches/2016/dev_CNRS_2016/NEMOGCM/CONFIG/cfg.txt
r7297 r7355 11 11 GYRE OPA_SRC 12 12 ISOMIP OPA_SRC 13 ORCA2_LIM3_PISCES OPA_SRC LIM_SRC_3 TOP_SRC NST_SRC 13 14 ORCA2_LIM3 OPA_SRC LIM_SRC_3 NST_SRC 14 ORCA2_LIM3_PISCES OPA_SRC LIM_SRC_3 TOP_SRC NST_SRC 15 SAS_BIPER OPA_SRC SAS_SRC LIM_SRC_3 NST_SRC 16 GYRE_LONG OPA_SRC 17 GYRE_32 OPA_SRC 18 ORCA2LIM3_LONG OPA_SRC LIM_SRC_3 NST_SRC -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_2/iceini_2.F90
r6140 r7355 83 83 CALL ice_run_2 ! read in namelist some run parameters 84 84 ! 85 rdt_ice = nn_fsbc * rdt 85 rdt_ice = nn_fsbc * rdt ! sea-ice time step 86 86 numit = nit000 - 1 87 87 ! -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/ice.F90
r7278 r7355 146 146 !! smt_i | - | Mean sea ice salinity | ppt | 147 147 !! tm_i | - | Mean sea ice temperature | K | 148 !! ot_i ! - ! Sea ice areal age content | day |149 148 !! et_i ! - ! Total ice enthalpy | J/m2 | 150 149 !! et_s ! - ! Total snow enthalpy | J/m2 | 151 !! bv_i ! - ! Mean relative brine volume| ??? |150 !! bv_i ! - ! relative brine volume | ??? | 152 151 !!===================================================================== 153 152 … … 157 156 !! * Share Module variables 158 157 !!-------------------------------------------------------------------------- 158 ! !!** ice-generic parameters namelist (namicerun) ** 159 INTEGER , PUBLIC :: jpl !: number of ice categories 160 INTEGER , PUBLIC :: nlay_i !: number of ice layers 161 INTEGER , PUBLIC :: nlay_s !: number of snow layers 162 REAL(wp) , PUBLIC :: rn_amax_n !: maximum ice concentration Northern hemisphere 163 REAL(wp) , PUBLIC :: rn_amax_s !: maximum ice concentration Southern hemisphere 164 CHARACTER(len=32) , PUBLIC :: cn_icerst_in !: suffix of ice restart name (input) 165 CHARACTER(len=32) , PUBLIC :: cn_icerst_out !: suffix of ice restart name (output) 166 CHARACTER(len=256), PUBLIC :: cn_icerst_indir !: ice restart input directory 167 CHARACTER(len=256), PUBLIC :: cn_icerst_outdir!: ice restart output directory 168 LOGICAL , PUBLIC :: ln_limthd !: flag for ice thermo (T) or not (F) 169 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F) 170 INTEGER , PUBLIC :: nn_limdyn !: flag for ice dynamics 171 REAL(wp) , PUBLIC :: rn_uice !: prescribed u-vel (case nn_limdyn=0) 172 REAL(wp) , PUBLIC :: rn_vice !: prescribed v-vel (case nn_limdyn=0) 173 174 ! !!** ice-diagnostics namelist (namicediag) ** 175 LOGICAL , PUBLIC :: ln_limdiachk !: flag for ice diag (T) or not (F) 176 LOGICAL , PUBLIC :: ln_limdiahsb !: flag for ice diag (T) or not (F) 177 LOGICAL , PUBLIC :: ln_limctl !: flag for sea-ice points output (T) or not (F) 178 INTEGER , PUBLIC :: iiceprt !: debug i-point 179 INTEGER , PUBLIC :: jiceprt !: debug j-point 180 181 ! !!** ice-init namelist (namiceini) ** 182 ! -- limistate -- ! 183 LOGICAL , PUBLIC :: ln_limini ! initialization or not 184 LOGICAL , PUBLIC :: ln_limini_file ! Ice initialization state from 2D netcdf file 185 REAL(wp), PUBLIC :: rn_thres_sst ! threshold water temperature for initial sea ice 186 REAL(wp), PUBLIC :: rn_hts_ini_n ! initial snow thickness in the north 187 REAL(wp), PUBLIC :: rn_hts_ini_s ! initial snow thickness in the south 188 REAL(wp), PUBLIC :: rn_hti_ini_n ! initial ice thickness in the north 189 REAL(wp), PUBLIC :: rn_hti_ini_s ! initial ice thickness in the south 190 REAL(wp), PUBLIC :: rn_ati_ini_n ! initial leads area in the north 191 REAL(wp), PUBLIC :: rn_ati_ini_s ! initial leads area in the south 192 REAL(wp), PUBLIC :: rn_smi_ini_n ! initial salinity 193 REAL(wp), PUBLIC :: rn_smi_ini_s ! initial salinity 194 REAL(wp), PUBLIC :: rn_tmi_ini_n ! initial temperature 195 REAL(wp), PUBLIC :: rn_tmi_ini_s ! initial temperature 196 197 ! !!** ice-thickness distribution namelist (namiceitd) ** 198 INTEGER , PUBLIC :: nn_catbnd !: categories distribution following: tanh function (1), or h^(-alpha) function (2) 199 REAL(wp), PUBLIC :: rn_himean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only) 200 201 ! !!** ice-dynamics namelist (namicedyn) ** 202 ! -- limtrp & limadv -- ! 203 INTEGER , PUBLIC :: nn_limadv !: choose the advection scheme (-1=Prather ; 0=Ultimate-Macho) 204 INTEGER , PUBLIC :: nn_limadv_ord !: choose the order of the advection scheme (if Ultimate-Macho) 205 ! -- limitd_me -- ! 206 INTEGER , PUBLIC :: nn_icestr !: ice strength parameterization (0=Hibler79 1=Rothrock75) 207 REAL(wp), PUBLIC :: rn_pe_rdg !: ridging work divided by pot. energy change in ridging, nn_icestr = 1 208 REAL(wp), PUBLIC :: rn_pstar !: determines ice strength, Hibler JPO79 209 REAL(wp), PUBLIC :: rn_crhg !: determines changes in ice strength 210 LOGICAL , PUBLIC :: ln_icestr_bvf !: use brine volume to diminish ice strength 211 ! -- limdyn & limrhg -- ! 212 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress 213 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9 214 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve 215 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 216 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 217 LOGICAL , PUBLIC :: ln_landfast !: landfast ice parameterization (T or F) 218 REAL(wp), PUBLIC :: rn_gamma !: fraction of ocean depth that ice must reach to initiate landfast ice 219 REAL(wp), PUBLIC :: rn_icebfr !: maximum bottom stress per unit area of contact (landfast ice) 220 REAL(wp), PUBLIC :: rn_lfrelax !: relaxation time scale (s-1) to reach static friction (landfast ice) 221 222 ! !!** ice-diffusion namelist (namicehdf) ** 223 INTEGER , PUBLIC :: nn_ahi0 !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation) 224 REAL(wp), PUBLIC :: rn_ahi0_ref !: sea-ice hor. eddy diffusivity coeff. (m2/s) 225 226 ! !!** ice-thermodynamics namelist (namicethd) ** 227 ! -- limthd_dif -- ! 228 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m] 229 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion 230 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion 231 INTEGER , PUBLIC :: nn_ice_thcon !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007) 232 LOGICAL , PUBLIC :: ln_it_qnsice !: iterate surface flux with changing surface temperature or not (F) 233 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1) or not (0) 234 ! -- limthd_dh -- ! 235 LOGICAL , PUBLIC :: ln_limdH !: activate ice thickness change from growing/melting (T) or not (F) 236 REAL(wp), PUBLIC :: rn_betas !: coef. for partitioning of snowfall between leads and sea ice 237 ! -- limthd_da -- ! 238 LOGICAL , PUBLIC :: ln_limdA !: activate lateral melting param. (T) or not (F) 239 REAL(wp), PUBLIC :: rn_beta !: coef. beta for lateral melting param. 240 REAL(wp), PUBLIC :: rn_dmin !: minimum floe diameter for lateral melting param. 241 ! -- limthd_lac -- ! 242 LOGICAL , PUBLIC :: ln_limdO !: activate ice growth in open-water (T) or not (F) 243 REAL(wp), PUBLIC :: rn_hnewice !: thickness for new ice formation (m) 244 LOGICAL , PUBLIC :: ln_frazil !: use of frazil ice collection as function of wind (T) or not (F) 245 REAL(wp), PUBLIC :: rn_maxfrazb !: maximum portion of frazil ice collecting at the ice bottom 246 REAL(wp), PUBLIC :: rn_vfrazb !: threshold drift speed for collection of bottom frazil ice 247 REAL(wp), PUBLIC :: rn_Cfrazb !: squeezing coefficient for collection of bottom frazil ice 248 ! -- limitd_th -- ! 249 REAL(wp), PUBLIC :: rn_himin !: minimum ice thickness 250 251 ! !!** ice-salinity namelist (namicesal) ** 252 LOGICAL , PUBLIC :: ln_limdS !: activate gravity drainage and flushing (T) or not (F) 253 INTEGER , PUBLIC :: nn_icesal !: salinity configuration used in the model 254 ! ! 1 - constant salinity in both space and time 255 ! ! 2 - prognostic salinity (s(z,t)) 256 ! ! 3 - salinity profile, constant in time 257 REAL(wp), PUBLIC :: rn_icesal !: bulk salinity (ppt) in case of constant salinity 258 REAL(wp), PUBLIC :: rn_sal_gd !: restoring salinity for gravity drainage [PSU] 259 REAL(wp), PUBLIC :: rn_time_gd !: restoring time constant for gravity drainage (= 20 days) [s] 260 REAL(wp), PUBLIC :: rn_sal_fl !: restoring salinity for flushing [PSU] 261 REAL(wp), PUBLIC :: rn_time_fl !: restoring time constant for gravity drainage (= 10 days) [s] 262 REAL(wp), PUBLIC :: rn_simax !: maximum ice salinity [PSU] 263 REAL(wp), PUBLIC :: rn_simin !: minimum ice salinity [PSU] 264 265 ! !!** ice-mechanical redistribution namelist (namiceitdme) 266 REAL(wp), PUBLIC :: rn_cs !: fraction of shearing energy contributing to ridging 267 INTEGER , PUBLIC :: nn_partfun !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007) 268 REAL(wp), PUBLIC :: rn_gstar !: fractional area of young ice contributing to ridging 269 REAL(wp), PUBLIC :: rn_astar !: equivalent of G* for an exponential participation function 270 LOGICAL , PUBLIC :: ln_ridging !: ridging of ice or not 271 REAL(wp), PUBLIC :: rn_hstar !: thickness that determines the maximal thickness of ridged ice 272 REAL(wp), PUBLIC :: rn_por_rdg !: initial porosity of ridges (0.3 regular value) 273 REAL(wp), PUBLIC :: rn_fsnowrdg !: fractional snow loss to the ocean during ridging 274 LOGICAL , PUBLIC :: ln_rafting !: rafting of ice or not 275 REAL(wp), PUBLIC :: rn_hraft !: threshold thickness (m) for rafting / ridging 276 REAL(wp), PUBLIC :: rn_craft !: coefficient for smoothness of the hyperbolic tangent in rafting 277 REAL(wp), PUBLIC :: rn_fsnowrft !: fractional snow loss to the ocean during ridging 278 279 ! !!** some other parameters 159 280 INTEGER , PUBLIC :: nstart !: iteration number of the begining of the run 160 281 INTEGER , PUBLIC :: nlast !: iteration number of the end of the run … … 163 284 REAL(wp), PUBLIC :: rdt_ice !: ice time step 164 285 REAL(wp), PUBLIC :: r1_rdtice !: = 1. / rdt_ice 165 166 ! !!** ice-thickness distribution namelist (namiceitd) **167 INTEGER , PUBLIC :: nn_catbnd !: categories distribution following: tanh function (1), or h^(-alpha) function (2)168 REAL(wp), PUBLIC :: rn_himean !: mean thickness of the domain (used to compute the distribution, nn_itdshp = 2 only)169 170 ! !!** ice-dynamics namelist (namicedyn) **171 LOGICAL , PUBLIC :: ln_icestr_bvf !: use brine volume to diminish ice strength172 INTEGER , PUBLIC :: nn_icestr !: ice strength parameterization (0=Hibler79 1=Rothrock75)173 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling174 INTEGER , PUBLIC :: nn_ahi0 !: sea-ice hor. eddy diffusivity coeff. (3 ways of calculation)175 REAL(wp), PUBLIC :: rn_pe_rdg !: ridging work divided by pot. energy change in ridging, nn_icestr = 1176 REAL(wp), PUBLIC :: rn_cio !: drag coefficient for oceanic stress177 REAL(wp), PUBLIC :: rn_pstar !: determines ice strength (N/M), Hibler JPO79178 REAL(wp), PUBLIC :: rn_crhg !: determines changes in ice strength179 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9180 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve181 REAL(wp), PUBLIC :: rn_ahi0_ref !: sea-ice hor. eddy diffusivity coeff. (m2/s)182 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rdt_ice (1/3 or 1/9 depending on nb of subcycling nevp)183 184 ! !!** ice-salinity namelist (namicesal) **185 REAL(wp), PUBLIC :: rn_simax !: maximum ice salinity [PSU]186 REAL(wp), PUBLIC :: rn_simin !: minimum ice salinity [PSU]187 REAL(wp), PUBLIC :: rn_sal_gd !: restoring salinity for gravity drainage [PSU]188 REAL(wp), PUBLIC :: rn_sal_fl !: restoring salinity for flushing [PSU]189 REAL(wp), PUBLIC :: rn_time_gd !: restoring time constant for gravity drainage (= 20 days) [s]190 REAL(wp), PUBLIC :: rn_time_fl !: restoring time constant for gravity drainage (= 10 days) [s]191 REAL(wp), PUBLIC :: rn_icesal !: bulk salinity (ppt) in case of constant salinity192 193 ! !!** ice-salinity namelist (namicesal) **194 INTEGER , PUBLIC :: nn_icesal !: salinity configuration used in the model195 ! ! 1 - constant salinity in both space and time196 ! ! 2 - prognostic salinity (s(z,t))197 ! ! 3 - salinity profile, constant in time198 INTEGER , PUBLIC :: nn_ice_thcon !: thermal conductivity: =0 Untersteiner (1964) ; =1 Pringle et al (2007)199 INTEGER , PUBLIC :: nn_monocat !: virtual ITD mono-category parameterizations (1) or not (0)200 LOGICAL , PUBLIC :: ln_it_qnsice !: iterate surface flux with changing surface temperature or not (F)201 202 ! !!** ice-mechanical redistribution namelist (namiceitdme)203 REAL(wp), PUBLIC :: rn_cs !: fraction of shearing energy contributing to ridging204 REAL(wp), PUBLIC :: rn_fsnowrdg !: fractional snow loss to the ocean during ridging205 REAL(wp), PUBLIC :: rn_fsnowrft !: fractional snow loss to the ocean during ridging206 REAL(wp), PUBLIC :: rn_gstar !: fractional area of young ice contributing to ridging207 REAL(wp), PUBLIC :: rn_astar !: equivalent of G* for an exponential participation function208 REAL(wp), PUBLIC :: rn_hstar !: thickness that determines the maximal thickness of ridged ice209 REAL(wp), PUBLIC :: rn_hraft !: threshold thickness (m) for rafting / ridging210 REAL(wp), PUBLIC :: rn_craft !: coefficient for smoothness of the hyperbolic tangent in rafting211 REAL(wp), PUBLIC :: rn_por_rdg !: initial porosity of ridges (0.3 regular value)212 REAL(wp), PUBLIC :: rn_betas !: coef. for partitioning of snowfall between leads and sea ice213 REAL(wp), PUBLIC :: rn_kappa_i !: coef. for the extinction of radiation Grenfell et al. (2006) [1/m]214 REAL(wp), PUBLIC :: nn_conv_dif !: maximal number of iterations for heat diffusion215 REAL(wp), PUBLIC :: rn_terr_dif !: maximal tolerated error (C) for heat diffusion216 217 ! !!** ice-mechanical redistribution namelist (namiceitdme)218 LOGICAL , PUBLIC :: ln_rafting !: rafting of ice or not219 INTEGER , PUBLIC :: nn_partfun !: participation function: =0 Thorndike et al. (1975), =1 Lipscomb et al. (2007)220 221 REAL(wp), PUBLIC :: usecc2 !: = 1.0 / ( rn_ecc * rn_ecc )222 REAL(wp), PUBLIC :: rhoco !: = rau0 * cio223 286 REAL(wp), PUBLIC :: r1_nlay_i !: 1 / nlay_i 224 287 REAL(wp), PUBLIC :: r1_nlay_s !: 1 / nlay_s 225 ! 226 ! !!** switch for presence of ice or not 227 REAL(wp), PUBLIC :: rswitch 228 ! 229 ! !!** define some parameters 288 REAL(wp), PUBLIC :: rswitch !: switch for the presence of ice (1) or not (0) 230 289 REAL(wp), PUBLIC, PARAMETER :: epsi06 = 1.e-06_wp !: small number 231 290 REAL(wp), PUBLIC, PARAMETER :: epsi10 = 1.e-10_wp !: small number 232 291 REAL(wp), PUBLIC, PARAMETER :: epsi20 = 1.e-20_wp !: small number 233 292 234 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_oce, v_oce !: surface ocean velocity used in ice dynamics235 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahiu , ahiv !: hor. diffusivity coeff. at U- and V-points [m2/s]236 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ust2s, hicol !: friction velocity, ice collection thickness accreted in leads237 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: strp1, strp2 !: strength at previous time steps238 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: strength 293 ! !!** define arrays 294 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_oce, v_oce !: surface ocean velocity used in ice dynamics 295 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ahiu , ahiv !: hor. diffusivity coeff. at U- and V-points [m2/s] 296 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hicol !: ice collection thickness accreted in leads 297 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: strength !: ice strength 239 298 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: stress1_i, stress2_i, stress12_i !: 1st, 2nd & diagonal stress tensor element 240 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_i 241 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: divu_i 242 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: shear_i 299 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: delta_i !: ice rheology elta factor (Flato & Hibler 95) [s-1] 300 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: divu_i !: Divergence of the velocity field [s-1] 301 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: shear_i !: Shear of the velocity field [s-1] 243 302 ! 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 303 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_bo !: Sea-Ice bottom temperature [Kelvin] 247 304 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: frld !: Leads fraction = 1 - ice fraction … … 252 309 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: fhld !: heat flux from the lead used for bottom melting 253 310 254 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw !: snow-ocean mass exchange [kg.m-2.s-1] 255 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_spr !: snow precipitation on ice [kg.m-2.s-1] 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sub !: snow/ice sublimation [kg.m-2.s-1] 257 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ice !: ice-ocean mass exchange [kg.m-2.s-1] 259 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sni !: snow ice growth component of wfx_ice [kg.m-2.s-1] 260 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_opw !: lateral ice growth component of wfx_ice [kg.m-2.s-1] 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_bog !: bottom ice growth component of wfx_ice [kg.m-2.s-1] 262 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_dyn !: dynamical ice growth component of wfx_ice [kg.m-2.s-1] 263 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_bom !: bottom melt component of wfx_ice [kg.m-2.s-1] 264 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sum !: surface melt component of wfx_ice [kg.m-2.s-1] 265 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_res !: residual component of wfx_ice [kg.m-2.s-1] 311 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_snw !: snow-ocean mass exchange [kg.m-2.s-1] 312 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_spr !: snow precipitation on ice [kg.m-2.s-1] 313 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sub !: snow/ice sublimation [kg.m-2.s-1] 314 315 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_ice !: ice-ocean mass exchange [kg.m-2.s-1] 316 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sni !: snow ice growth component of wfx_ice [kg.m-2.s-1] 317 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_opw !: lateral ice growth component of wfx_ice [kg.m-2.s-1] 318 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_bog !: bottom ice growth component of wfx_ice [kg.m-2.s-1] 319 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_dyn !: dynamical ice growth component of wfx_ice [kg.m-2.s-1] 320 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_bom !: bottom melt component of wfx_ice [kg.m-2.s-1] 321 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_sum !: surface melt component of wfx_ice [kg.m-2.s-1] 322 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_lam !: lateral melt component of wfx_ice [kg.m-2.s-1] 323 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: wfx_res !: residual component of wfx_ice [kg.m-2.s-1] 266 324 267 325 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: afx_tot !: ice concentration tendency (total) [s-1] … … 271 329 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bog !: salt flux due to ice growth/melt [PSU/m2/s] 272 330 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_bom !: salt flux due to ice growth/melt [PSU/m2/s] 331 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_lam !: salt flux due to ice growth/melt [PSU/m2/s] 273 332 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_sum !: salt flux due to ice growth/melt [PSU/m2/s] 274 333 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: sfx_sni !: salt flux due to ice growth/melt [PSU/m2/s] … … 302 361 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: hfx_res !: residual heat flux due to correction of ice thickness [W.m-2] 303 362 304 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ,:) :: ftr_ice !: transmitted solar radiation under ice305 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: pahu3D , pahv3D306 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,: ) :: rn_amax_2d !: maximum ice concentration 2d array363 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rn_amax_2d !: maximum ice concentration 2d array 364 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ftr_ice !: transmitted solar radiation under ice 365 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: pahu3D, pahv3D !: ice hor. eddy diffusivity coef. at U- and V-points 307 366 308 367 !!-------------------------------------------------------------------------- … … 310 369 !!-------------------------------------------------------------------------- 311 370 !! Variables defined for each ice category 312 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ht_i !: Ice thickness (m)313 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i !: Ice fractional areas (concentration)314 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_i !: Ice volume per unit area (m)315 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s !: Snow volume per unit area(m)316 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ht_s !: Snow thickness (m)317 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: t_su !: Sea-Ice Surface Temperature (K)318 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sm_i !: Sea-Ice Bulk salinity (ppt)319 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: smv_i !: Sea-Ice Bulk salinity times volume per area (ppt.m)320 ! ! this is an extensive variable that has to be transported321 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: o_i !: Sea-Ice Age (days)322 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: o v_i !: Sea-Ice Age times volume per area (days.m)323 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: oa_i !: Sea-Ice Age times ice area (days)371 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ht_i !: Ice thickness (m) 372 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i !: Ice fractional areas (concentration) 373 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_i !: Ice volume per unit area (m) 374 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s !: Snow volume per unit area(m) 375 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: ht_s !: Snow thickness (m) 376 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: t_su !: Sea-Ice Surface Temperature (K) 377 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: sm_i !: Sea-Ice Bulk salinity (ppt) 378 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: smv_i !: Sea-Ice Bulk salinity times volume per area (ppt.m) 379 ! ! this is an extensive variable that has to be transported 380 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: o_i !: Sea-Ice Age (days) 381 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: oa_i !: Sea-Ice Age times ice area (days) 382 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: bv_i !: brine volume 324 383 325 384 !! Variables summed over all categories, or associated to all the ice in a single grid cell 326 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 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_i , vt_s !: ice and snow total volume per unit area (m) 329 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i !: ice total fractional area (ice concentration) 330 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ato_i !: =1-at_i ; total open water fractional area 331 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] 336 337 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_s !: Snow temperatures [K] 338 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s !: Snow ... 385 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice, v_ice !: components of the ice velocity (m/s) 386 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_i , vt_s !: ice and snow total volume per unit area (m) 387 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i !: ice total fractional area (ice concentration) 388 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: ato_i !: =1-at_i ; total open water fractional area 389 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: et_i , et_s !: ice and snow total heat content 390 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tm_i !: mean ice temperature over all categories 391 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: bvm_i !: brine volume averaged over all categories 392 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: smt_i !: mean sea ice salinity averaged over all categories [PSU] 393 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tm_su !: mean surface temperature over all categories 394 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htm_i !: mean ice thickness over all categories 395 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: htm_s !: mean snow thickness over all categories 396 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: om_i !: mean ice age over all categories 397 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tau_icebfr !: ice friction with bathy (landfast param activated) 398 399 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_s !: Snow temperatures [K] 400 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s !: Snow ... 339 401 340 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_i 341 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i 342 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: s_i 402 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_i !: ice temperatures [K] 403 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i !: ice thermal contents [J/m2] 404 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: s_i !: ice salinities [PSU] 343 405 344 406 !!-------------------------------------------------------------------------- … … 362 424 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_i_b !: ice temperatures 363 425 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: u_ice_b, v_ice_b !: ice velocity 426 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: at_i_b !: ice concentration (total) 364 427 365 428 !!-------------------------------------------------------------------------- … … 368 431 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_max !: Boundary of ice thickness categories in thickness space 369 432 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:) :: hi_mean !: Mean ice thickness in catgories 370 371 !!--------------------------------------------------------------------------372 !! * Ice Run373 !!--------------------------------------------------------------------------374 ! !!: ** Namelist namicerun read in sbc_lim_init **375 INTEGER , PUBLIC :: jpl !: number of ice categories376 INTEGER , PUBLIC :: nlay_i !: number of ice layers377 INTEGER , PUBLIC :: nlay_s !: number of snow layers378 CHARACTER(len=32) , PUBLIC :: cn_icerst_in !: suffix of ice restart name (input)379 CHARACTER(len=256), PUBLIC :: cn_icerst_indir !: ice restart input directory380 CHARACTER(len=32) , PUBLIC :: cn_icerst_out !: suffix of ice restart name (output)381 CHARACTER(len=256), PUBLIC :: cn_icerst_outdir!: ice restart output directory382 LOGICAL , PUBLIC :: ln_limdyn !: flag for ice dynamics (T) or not (F)383 LOGICAL , PUBLIC :: ln_icectl !: flag for sea-ice points output (T) or not (F)384 REAL(wp) , PUBLIC :: rn_amax_n !: maximum ice concentration Northern hemisphere385 REAL(wp) , PUBLIC :: rn_amax_s !: maximum ice concentration Southern hemisphere386 INTEGER , PUBLIC :: iiceprt !: debug i-point387 INTEGER , PUBLIC :: jiceprt !: debug j-point388 433 ! 389 434 !!-------------------------------------------------------------------------- 390 435 !! * Ice diagnostics 391 436 !!-------------------------------------------------------------------------- 392 ! Increment of global variables393 437 ! thd refers to changes induced by thermodynamics 394 438 ! trp '' '' '' advection (transport of ice) 395 LOGICAL , PUBLIC :: ln_limdiahsb !: flag for ice diag (T) or not (F) 396 LOGICAL , PUBLIC :: ln_limdiaout !: flag for ice diag (T) or not (F) 439 ! 397 440 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vi !: transport of ice volume 398 441 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_trp_vs !: transport of snw volume … … 419 462 INTEGER :: ice_alloc 420 463 ! 421 INTEGER :: ierr(1 7), ii464 INTEGER :: ierr(15), ii 422 465 !!----------------------------------------------------------------- 423 466 … … 427 470 ! stay within Fortran's max-line length limit. 428 471 ii = 1 429 ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , & 430 & ahiu (jpi,jpj) , ahiv (jpi,jpj) , & 431 & ust2s (jpi,jpj) , hicol (jpi,jpj) , & 432 & strp1 (jpi,jpj) , strp2 (jpi,jpj) , strength (jpi,jpj) , & 433 & stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , & 434 & delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , STAT=ierr(ii) ) 435 436 ii = ii + 1 437 ALLOCATE( sist (jpi,jpj) , icethi (jpi,jpj) , t_bo (jpi,jpj) , & 438 & frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) , & 439 & wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , & 472 ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , & 473 & ahiu (jpi,jpj) , ahiv (jpi,jpj) , hicol (jpi,jpj) , & 474 & strength(jpi,jpj) , stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , & 475 & delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , STAT=ierr(ii) ) 476 477 ii = ii + 1 478 ALLOCATE( t_bo (jpi,jpj) , frld (jpi,jpj) , pfrld (jpi,jpj) , phicif (jpi,jpj) , & 479 & wfx_snw(jpi,jpj) , wfx_ice(jpi,jpj) , wfx_sub(jpi,jpj) , wfx_lam(jpi,jpj) , & 440 480 & wfx_bog(jpi,jpj) , wfx_dyn(jpi,jpj) , wfx_bom(jpi,jpj) , wfx_sum(jpi,jpj) , & 441 481 & wfx_res(jpi,jpj) , wfx_sni(jpi,jpj) , wfx_opw(jpi,jpj) , wfx_spr(jpi,jpj) , & 442 & afx_tot(jpi,jpj) , afx_thd(jpi,jpj), afx_dyn(jpi,jpj) , & 443 & fhtur (jpi,jpj) , ftr_ice(jpi,jpj,jpl), pahu3D(jpi,jpj,jpl+1), pahv3D(jpi,jpj,jpl+1), & 444 & rn_amax_2d (jpi,jpj) , qlead (jpi,jpj) , & 445 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj), & 446 & sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , & 482 & afx_tot(jpi,jpj) , afx_thd(jpi,jpj), afx_dyn(jpi,jpj) , rn_amax_2d(jpi,jpj), & 483 & fhtur (jpi,jpj) , qlead (jpi,jpj) , & 484 & sfx_res(jpi,jpj) , sfx_bri(jpi,jpj) , sfx_dyn(jpi,jpj) , sfx_sub(jpi,jpj) , sfx_lam(jpi,jpj) , & 485 & sfx_bog(jpi,jpj) , sfx_bom(jpi,jpj) , sfx_sum(jpi,jpj) , sfx_sni(jpi,jpj) , sfx_opw(jpi,jpj) , & 447 486 & hfx_res(jpi,jpj) , hfx_snw(jpi,jpj) , hfx_sub(jpi,jpj) , hfx_err(jpi,jpj) , & 448 & hfx_ err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj) ,&449 & hfx_ in (jpi,jpj) , hfx_out(jpi,jpj) , fhld(jpi,jpj) ,&450 & hfx_ sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , hfx_opw(jpi,jpj) ,&451 & hfx_ thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj), STAT=ierr(ii) )487 & hfx_in (jpi,jpj) , hfx_out(jpi,jpj) , fhld (jpi,jpj) , & 488 & hfx_sum(jpi,jpj) , hfx_bom(jpi,jpj) , hfx_bog(jpi,jpj) , hfx_dif(jpi,jpj) , & 489 & hfx_opw(jpi,jpj) , hfx_thd(jpi,jpj) , hfx_dyn(jpi,jpj) , hfx_spr(jpi,jpj) , & 490 & hfx_err_dif(jpi,jpj) , hfx_err_rem(jpi,jpj) , wfx_err_sub(jpi,jpj) , STAT=ierr(ii) ) 452 491 453 492 ! * Ice global state variables 454 493 ii = ii + 1 455 ALLOCATE( ht_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , & 456 & v_s (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , & 457 & sm_i (jpi,jpj,jpl) , smv_i(jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , & 458 & ov_i (jpi,jpj,jpl) , oa_i (jpi,jpj,jpl) , STAT=ierr(ii) ) 459 ii = ii + 1 460 ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , tio_u(jpi,jpj) , tio_v(jpi,jpj) , & 494 ALLOCATE( ftr_ice(jpi,jpj,jpl) , pahu3D(jpi,jpj,jpl+1) , pahv3D(jpi,jpj,jpl+1) , & 495 & ht_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , & 496 & v_s (jpi,jpj,jpl) , ht_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , & 497 & sm_i (jpi,jpj,jpl) , smv_i (jpi,jpj,jpl) , o_i (jpi,jpj,jpl) , & 498 & oa_i (jpi,jpj,jpl) , bv_i (jpi,jpj,jpl) , STAT=ierr(ii) ) 499 ii = ii + 1 500 ALLOCATE( u_ice(jpi,jpj) , v_ice(jpi,jpj) , & 461 501 & 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) ) 502 & et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i (jpi,jpj) , bvm_i(jpi,jpj) , & 503 & smt_i(jpi,jpj) , tm_su(jpi,jpj) , htm_i(jpi,jpj) , htm_s(jpi,jpj) , & 504 & om_i (jpi,jpj) , tau_icebfr(jpi,jpj) , STAT=ierr(ii) ) 464 505 ii = ii + 1 465 506 ALLOCATE( t_s(jpi,jpj,nlay_s,jpl) , e_s(jpi,jpj,nlay_s,jpl) , STAT=ierr(ii) ) … … 488 529 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 489 530 & a_i_b (jpi,jpj,jpl) , smv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , & 490 & oa_i_b (jpi,jpj,jpl) , u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , STAT=ierr(ii) ) 531 & oa_i_b (jpi,jpj,jpl) , STAT=ierr(ii) ) 532 ii = ii + 1 533 ALLOCATE( u_ice_b(jpi,jpj) , v_ice_b(jpi,jpj) , at_i_b(jpi,jpj) , STAT=ierr(ii) ) 491 534 492 535 ! * Ice thickness distribution variables … … 496 539 ! * Ice diagnostics 497 540 ii = ii + 1 498 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj), diag_trp_ei(jpi,jpj), &499 & diag_trp_es(jpi,jpj) , diag_trp_smv(jpi,jpj), diag_heat (jpi,jpj), &500 & diag_smvi (jpi,jpj) , diag_vice (jpi,jpj), diag_vsnw (jpi,jpj), STAT=ierr(ii) )541 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 542 & diag_trp_es(jpi,jpj) , diag_trp_smv(jpi,jpj) , diag_heat (jpi,jpj), & 543 & diag_smvi (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), STAT=ierr(ii) ) 501 544 502 545 ice_alloc = MAXVAL( ierr(:) ) 503 IF( ice_alloc /= 0 ) CALL ctl_warn('ice_alloc _2: failed to allocate arrays.')546 IF( ice_alloc /= 0 ) CALL ctl_warn('ice_alloc: failed to allocate arrays.') 504 547 ! 505 548 END FUNCTION ice_alloc … … 513 556 !!====================================================================== 514 557 END MODULE ice 515 -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limcons.F90
r7278 r7355 18 18 USE phycst ! physical constants 19 19 USE ice ! LIM-3 variables 20 USE dom_ice ! LIM-3 domain21 20 USE dom_oce ! ocean domain 22 21 USE in_out_manager ! I/O manager … … 165 164 !! + test if ice concentration and volume are > 0 166 165 !! 167 !! ** Method : This is an online diagnostics which can be activated with ln_limdia hsb=true166 !! ** Method : This is an online diagnostics which can be activated with ln_limdiachk=true 168 167 !! It prints in ocean.output if there is a violation of conservation at each time-step 169 168 !! The thresholds (zv_sill, zs_sill, zh_sill) which determine violations are set to … … 185 184 ! salt flux 186 185 zfs_b = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 187 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) 186 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 188 187 & ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) 189 188 190 189 ! water flux 191 zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + &192 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) 190 zfw_b = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 191 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) + wfx_lam(:,:) & 193 192 & ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) 194 193 … … 210 209 ! salt flux 211 210 zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + & 212 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) 211 & sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) & 213 212 & ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfs_b 214 213 215 214 ! water flux 216 zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + &217 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) &215 zfw = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + wfx_opw(:,:) + & 216 & wfx_res(:,:) + wfx_dyn(:,:) + wfx_snw(:,:) + wfx_sub(:,:) + wfx_spr(:,:) + wfx_lam(:,:) & 218 217 & ) * e1e2t(:,:) * tmask(:,:,1) * zconv ) - zfw_b 219 218 … … 260 259 & cd_routine /= 'limtrp' .AND. cd_routine /= 'limitd_me' ) THEN 261 260 WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax 261 IF ( zamax > 1._wp ) WRITE(numout,*) 'violation a_i>1 (',cd_routine,') = ',zamax 262 262 ENDIF 263 263 IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin … … 274 274 !! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step 275 275 !! 276 !! ** Method : This is an online diagnostics which can be activated with ln_limdia hsb=true276 !! ** Method : This is an online diagnostics which can be activated with ln_limdiachk=true 277 277 !! It prints in ocean.output if there is a violation of conservation at each time-step 278 278 !! The thresholds (zv_sill, zs_sill, zh_sill) which determine the violation are set to … … 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 & * e1e2t * 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 & ) * e1e2t * tmask(:,:,1) * zconv ) 292 293 ! salt flux 293 294 zsfx = glob_sum( ( sfx + diag_smvi ) * e1e2t * tmask(:,:,1) * zconv ) * rday -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limctl.F90
r5836 r7355 5 5 !!====================================================================== 6 6 !! History : 3.5 ! 2015-01 (M. Vancoppenolle) Original code 7 !! 3.7 ! 2016-10 (C. Rousset) Add routine lim_prt3D 7 8 !!---------------------------------------------------------------------- 8 9 #if defined key_lim3 … … 12 13 !! lim_ctl : control prints in case of crash 13 14 !! lim_prt : ice control print at a given grid point 15 !! lim_prt3D : control prints of ice arrays 14 16 !!---------------------------------------------------------------------- 15 17 USE oce ! ocean dynamics and tracers … … 17 19 USE ice ! LIM-3: ice variables 18 20 USE thd_ice ! LIM-3: thermodynamical variables 19 USE dom_ice ! LIM-3: ice domain20 21 USE sbc_oce ! Surface boundary condition: ocean fields 21 22 USE sbc_ice ! Surface boundary condition: ice fields … … 35 36 PUBLIC lim_ctl 36 37 PUBLIC lim_prt 38 PUBLIC lim_prt3D 37 39 38 40 !! * Substitutions … … 445 447 END SUBROUTINE lim_prt 446 448 449 SUBROUTINE lim_prt3D( cd_routine ) 450 !!--------------------------------------------------------------------------------------------------------- 451 !! *** ROUTINE lim_prt3D *** 452 !! 453 !! ** Purpose : CTL prints of ice arrays in case ln_ctl is activated 454 !! 455 !!--------------------------------------------------------------------------------------------------------- 456 CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine 457 INTEGER :: jk, jl ! dummy loop indices 458 459 CALL prt_ctl_info(' ========== ') 460 CALL prt_ctl_info( cd_routine ) 461 CALL prt_ctl_info(' ========== ') 462 CALL prt_ctl_info(' - Cell values : ') 463 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 464 CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' cell area :') 465 CALL prt_ctl(tab2d_1=at_i , clinfo1=' at_i :') 466 CALL prt_ctl(tab2d_1=ato_i , clinfo1=' ato_i :') 467 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' vt_i :') 468 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' vt_s :') 469 CALL prt_ctl(tab2d_1=divu_i , clinfo1=' divu_i :') 470 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :') 471 CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' stress1_i :') 472 CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' stress2_i :') 473 CALL prt_ctl(tab2d_1=stress12_i , clinfo1=' stress12_i :') 474 CALL prt_ctl(tab2d_1=strength , clinfo1=' strength :') 475 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :') 476 CALL prt_ctl(tab2d_1=u_ice , clinfo1=' u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :') 477 478 DO jl = 1, jpl 479 CALL prt_ctl_info(' ') 480 CALL prt_ctl_info(' - Category : ', ivar1=jl) 481 CALL prt_ctl_info(' ~~~~~~~~~~') 482 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' ht_i : ') 483 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' ht_s : ') 484 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' t_su : ') 485 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' t_snow : ') 486 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' sm_i : ') 487 CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' o_i : ') 488 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' a_i : ') 489 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' v_i : ') 490 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' v_s : ') 491 CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' e_i1 : ') 492 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' e_snow : ') 493 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' smv_i : ') 494 CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' oa_i : ') 495 496 DO jk = 1, nlay_i 497 CALL prt_ctl_info(' - Layer : ', ivar1=jk) 498 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i : ') 499 END DO 500 END DO 501 502 CALL prt_ctl_info(' ') 503 CALL prt_ctl_info(' - Heat / FW fluxes : ') 504 CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ') 505 CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' sst : ', tab2d_2=sss_m , clinfo2= ' sss : ') 506 CALL prt_ctl(tab2d_1=qsr , clinfo1= ' qsr : ', tab2d_2=qns , clinfo2= ' qns : ') 507 CALL prt_ctl(tab2d_1=emp , clinfo1= ' emp : ', tab2d_2=sfx , clinfo2= ' sfx : ') 508 509 CALL prt_ctl_info(' ') 510 CALL prt_ctl_info(' - Stresses : ') 511 CALL prt_ctl_info(' ~~~~~~~~~~ ') 512 CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', tab2d_2=vtau , clinfo2= ' vtau : ') 513 CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ') 514 CALL prt_ctl(tab2d_1=u_oce , clinfo1= ' u_oce : ', tab2d_2=v_oce , clinfo2= ' v_oce : ') 515 516 END SUBROUTINE lim_prt3D 517 447 518 #else 448 519 !!-------------------------------------------------------------------------- … … 454 525 SUBROUTINE lim_prt ! Empty routine 455 526 END SUBROUTINE lim_prt 527 SUBROUTINE lim_prt3D ! Empty routine 528 END SUBROUTINE lim_prt3D 456 529 #endif 457 530 !!====================================================================== -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limdiahsb.F90
r7278 r7355 14 14 !!---------------------------------------------------------------------- 15 15 USE ice ! LIM-3: sea-ice variable 16 USE dom_ice ! LIM-3: sea-ice domain17 16 USE dom_oce ! ocean domain 18 17 USE sbc_oce ! surface boundary condition: ocean fields … … 31 30 32 31 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 32 PUBLIC lim_diahsb_init ! routine called in sbcice_lim.F90 33 34 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: vol_loc_ini, sal_loc_ini, tem_loc_ini ! initial volume, salt and heat contents 35 REAL(wp) :: frc_sal, frc_voltop, frc_volbot, frc_temtop, frc_tembot ! global forcing trends 36 37 37 !! * Substitutions 38 38 # include "vectopt_loop_substitute.h90" … … 46 46 CONTAINS 47 47 48 SUBROUTINE lim_diahsb 48 SUBROUTINE lim_diahsb( kt ) 49 49 !!--------------------------------------------------------------------------- 50 50 !! *** ROUTINE lim_diahsb *** … … 53 53 !! 54 54 !!--------------------------------------------------------------------------- 55 INTEGER, INTENT(in) :: kt ! number of iteration 55 56 !! 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 57 real(wp) :: zbg_ivol, zbg_svol, zbg_area, zbg_isal, zbg_item ,zbg_stem 58 REAL(wp) :: z_frc_voltop, z_frc_volbot, z_frc_sal, z_frc_temtop, z_frc_tembot 59 REAL(wp) :: zdiff_vol, zdiff_sal, zdiff_tem 67 60 !!--------------------------------------------------------------------------- 68 61 IF( nn_timing == 1 ) CALL timing_start('lim_diahsb') 69 62 70 IF( numit == nstart ) CALL lim_diahsb_init 71 72 ! 1/area 73 z1_area = 1._wp / MAX( glob_sum( e1e2t(:,:) * tmask(:,:,1) ), epsi06 ) 74 75 rswitch = MAX( 0._wp , SIGN( 1._wp , glob_sum( e1e2t(:,:) * tmask(:,:,1) ) - epsi06 ) ) 76 ! ----------------------- ! 77 ! 1 - Content variations ! 78 ! ----------------------- ! 79 zbg_ivo = glob_sum( vt_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! volume ice 80 zbg_svo = glob_sum( vt_s(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! volume snow 81 zbg_are = glob_sum( at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! area 82 zbg_sal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) * tmask(:,:,1) ) ! mean salt content 83 zbg_tem = glob_sum( ( tm_i(:,:) - rt0 ) * vt_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! mean temp content 84 85 !zbg_ihc = glob_sum( et_i(:,:) * e1e2t(:,:) * tmask(:,:,1) ) / MAX( zbg_ivo,epsi06 ) ! ice heat content 86 !zbg_shc = glob_sum( et_s(:,:) * e1e2t(:,:) * 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(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 91 zbg_vfx_bog = ztmp * glob_sum( wfx_bog(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 92 zbg_vfx_opw = ztmp * glob_sum( wfx_opw(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 93 zbg_vfx_sni = ztmp * glob_sum( wfx_sni(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 94 zbg_vfx_dyn = ztmp * glob_sum( wfx_dyn(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 95 zbg_vfx_bom = ztmp * glob_sum( wfx_bom(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 96 zbg_vfx_sum = ztmp * glob_sum( wfx_sum(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 97 zbg_vfx_res = ztmp * glob_sum( wfx_res(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 98 zbg_vfx_spr = ztmp * glob_sum( wfx_spr(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 99 zbg_vfx_snw = ztmp * glob_sum( wfx_snw(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 100 zbg_vfx_sub = ztmp * glob_sum( wfx_sub(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 101 102 ! Salt 103 zbg_sfx = ztmp * glob_sum( sfx(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 104 zbg_sfx_bri = ztmp * glob_sum( sfx_bri(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 105 zbg_sfx_res = ztmp * glob_sum( sfx_res(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 106 zbg_sfx_dyn = ztmp * glob_sum( sfx_dyn(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 107 108 zbg_sfx_bog = ztmp * glob_sum( sfx_bog(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 109 zbg_sfx_opw = ztmp * glob_sum( sfx_opw(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 110 zbg_sfx_sni = ztmp * glob_sum( sfx_sni(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 111 zbg_sfx_bom = ztmp * glob_sum( sfx_bom(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 112 zbg_sfx_sum = ztmp * glob_sum( sfx_sum(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 113 zbg_sfx_sub = ztmp * glob_sum( sfx_sub(:,:) * e1e2t(:,:) * tmask(:,:,1) ) 114 115 ! Heat budget 116 zbg_ihc = glob_sum( et_i(:,:) * e1e2t(:,:) * 1.e-20 ) ! ice heat content [1.e20 J] 117 zbg_shc = glob_sum( et_s(:,:) * e1e2t(:,:) * 1.e-20 ) ! snow heat content [1.e20 J] 118 zbg_hfx_dhc = glob_sum( diag_heat(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 119 zbg_hfx_spr = glob_sum( hfx_spr(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 120 121 zbg_hfx_thd = glob_sum( hfx_thd(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 122 zbg_hfx_dyn = glob_sum( hfx_dyn(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 123 zbg_hfx_res = glob_sum( hfx_res(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 124 zbg_hfx_sub = glob_sum( hfx_sub(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 125 zbg_hfx_snw = glob_sum( hfx_snw(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 126 zbg_hfx_sum = glob_sum( hfx_sum(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 127 zbg_hfx_bom = glob_sum( hfx_bom(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 128 zbg_hfx_bog = glob_sum( hfx_bog(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 129 zbg_hfx_dif = glob_sum( hfx_dif(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 130 zbg_hfx_opw = glob_sum( hfx_opw(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 131 zbg_hfx_out = glob_sum( hfx_out(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! [in W] 132 zbg_hfx_in = glob_sum( hfx_in(:,:) * e1e2t(:,:) * 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(:,:) * e1e2t(:,:) * tmask(:,:,1) ) ! volume fluxes 138 z_frc_sal = r1_rau0 * glob_sum( sfx(:,:) * e1e2t(:,:) * 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(:,:) ) * e1e2t(:,:) * 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 63 ! ----------------------- ! 64 ! 1 - Contents ! 65 ! ----------------------- ! 66 zbg_ivol = glob_sum( vt_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! ice volume (km3) 67 zbg_svol = glob_sum( vt_s(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! snow volume (km3) 68 zbg_area = glob_sum( at_i(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-6 ) ! area (km2) 69 zbg_isal = glob_sum( SUM( smv_i(:,:,:), dim=3 ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt content (pss*km3) 70 zbg_item = glob_sum( et_i * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat content (1.e20 J) 71 zbg_stem = glob_sum( et_s * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat content (1.e20 J) 146 72 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 73 ! ---------------------------! 74 ! 2 - Trends due to forcing ! 75 ! ---------------------------! 76 z_frc_volbot = r1_rau0 * glob_sum( - ( wfx_ice(:,:) + wfx_snw(:,:) + wfx_err_sub(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! freshwater flux ice/snow-ocean 77 z_frc_voltop = r1_rau0 * glob_sum( - ( wfx_sub(:,:) + wfx_spr(:,:) ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! freshwater flux ice/snow-atm 78 z_frc_sal = r1_rau0 * glob_sum( - sfx(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) ! salt fluxes ice/snow-ocean 79 z_frc_tembot = glob_sum( hfx_out(:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat on top of ocean (and below ice) 80 z_frc_temtop = glob_sum( hfx_in (:,:) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 ) ! heat on top of ice-coean 81 ! 82 frc_voltop = frc_voltop + z_frc_voltop * rdt_ice ! km3 83 frc_volbot = frc_volbot + z_frc_volbot * rdt_ice ! km3 84 frc_sal = frc_sal + z_frc_sal * rdt_ice ! km3*pss 85 frc_temtop = frc_temtop + z_frc_temtop * rdt_ice ! 1.e20 J 86 frc_tembot = frc_tembot + z_frc_tembot * rdt_ice ! 1.e20 J 87 88 ! ----------------------- ! 89 ! 3 - Content variations ! 90 ! ----------------------- ! 91 zdiff_vol = r1_rau0 * glob_sum( ( rhoic * vt_i(:,:) + rhosn * vt_s(:,:) - vol_loc_ini(:,:) & ! freshwater trend (km3) 92 & ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) 93 zdiff_sal = r1_rau0 * glob_sum( ( rhoic * SUM( smv_i(:,:,:), dim=3 ) - sal_loc_ini(:,:) & ! salt content trend (km3*pss) 94 & ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-9 ) 95 zdiff_tem = glob_sum( ( et_i(:,:) + et_s(:,:) - tem_loc_ini(:,:) & ! heat content trend (1.e20 J) 96 ! & + SUM( qevap_ice * a_i_b, dim=3 ) & !! clem: I think this line should be commented (but needs a check) 97 & ) * e1e2t(:,:) * tmask(:,:,1) * 1.e-20 ) 98 99 ! ----------------------- ! 100 ! 4 - Drifts ! 101 ! ----------------------- ! 102 zdiff_vol = zdiff_vol - ( frc_voltop + frc_volbot ) 103 zdiff_sal = zdiff_sal - frc_sal 104 zdiff_tem = zdiff_tem - ( frc_tembot - frc_temtop ) 105 106 ! ----------------------- ! 107 ! 5 - Diagnostics writing ! 108 ! ----------------------- ! 109 ! 110 IF( iom_use('ibgvolume') ) CALL iom_put( 'ibgvolume' , zdiff_vol ) ! ice/snow volume drift (km3 equivalent ocean water) 111 IF( iom_use('ibgsaltco') ) CALL iom_put( 'ibgsaltco' , zdiff_sal ) ! ice salt content drift (psu*km3 equivalent ocean water) 112 IF( iom_use('ibgheatco') ) CALL iom_put( 'ibgheatco' , zdiff_tem ) ! ice/snow heat content drift (1.e20 J) 113 IF( iom_use('ibgheatfx') ) CALL iom_put( 'ibgheatfx' , zdiff_tem / & ! ice/snow heat flux drift (W/m2) 114 & glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 115 116 IF( iom_use('ibgfrcvoltop') ) CALL iom_put( 'ibgfrcvoltop' , frc_voltop ) ! vol forcing ice/snw-atm (km3 equivalent ocean water) 117 IF( iom_use('ibgfrcvolbot') ) CALL iom_put( 'ibgfrcvolbot' , frc_volbot ) ! vol forcing ice/snw-ocean (km3 equivalent ocean water) 118 IF( iom_use('ibgfrcsal') ) CALL iom_put( 'ibgfrcsal' , frc_sal ) ! sal - forcing (psu*km3 equivalent ocean water) 119 IF( iom_use('ibgfrctemtop') ) CALL iom_put( 'ibgfrctemtop' , frc_temtop ) ! heat on top of ice/snw/ocean (1.e20 J) 120 IF( iom_use('ibgfrctembot') ) CALL iom_put( 'ibgfrctembot' , frc_tembot ) ! heat on top of ocean(below ice) (1.e20 J) 121 IF( iom_use('ibgfrchfxtop') ) CALL iom_put( 'ibgfrchfxtop' , frc_temtop / & ! heat on top of ice/snw/ocean (W/m2) 122 & glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 123 IF( iom_use('ibgfrchfxbot') ) CALL iom_put( 'ibgfrchfxbot' , frc_tembot / & ! heat on top of ocean(below ice) (W/m2) 124 & glob_sum( e1e2t(:,:) * tmask(:,:,1) * 1.e-20 * kt*rdt ) ) 125 126 IF( iom_use('ibgvol_tot' ) ) CALL iom_put( 'ibgvol_tot' , zbg_ivol ) ! ice volume (km3) 127 IF( iom_use('sbgvol_tot' ) ) CALL iom_put( 'sbgvol_tot' , zbg_svol ) ! snow volume (km3) 128 IF( iom_use('ibgarea_tot') ) CALL iom_put( 'ibgarea_tot' , zbg_area ) ! ice area (km2) 129 IF( iom_use('ibgsalt_tot') ) CALL iom_put( 'ibgsalt_tot' , zbg_isal ) ! ice salinity content (pss*km3) 130 IF( iom_use('ibgheat_tot') ) CALL iom_put( 'ibgheat_tot' , zbg_item ) ! ice heat content (1.e20 J) 131 IF( iom_use('sbgheat_tot') ) CALL iom_put( 'sbgheat_tot' , zbg_stem ) ! snow heat content (1.e20 J) 215 132 ! 216 133 IF( lrst_ice ) CALL lim_diahsb_rst( numit, 'WRITE' ) 217 134 ! 218 135 IF( nn_timing == 1 ) CALL timing_stop('lim_diahsb') 219 !136 ! 220 137 END SUBROUTINE lim_diahsb 221 138 … … 233 150 !! - Compute coefficients for conversion 234 151 !!--------------------------------------------------------------------------- 235 INTEGER :: jk ! dummy loop indice236 152 INTEGER :: ierror ! local integer 237 153 !! … … 247 163 WRITE(numout,*) '~~~~~~~~~~~~' 248 164 ENDIF 249 ! 165 ! 166 ALLOCATE( vol_loc_ini(jpi,jpj), sal_loc_ini(jpi,jpj), tem_loc_ini(jpi,jpj), STAT=ierror ) 167 IF( ierror > 0 ) THEN 168 CALL ctl_stop( 'lim_diahsb: unable to allocate vol_loc_ini' ) 169 RETURN 170 ENDIF 171 250 172 CALL lim_diahsb_rst( nstart, 'READ' ) !* read or initialize all required files 251 173 ! … … 263 185 CHARACTER(len=*), INTENT(in) :: cdrw ! "READ"/"WRITE" flag 264 186 ! 265 INTEGER :: id1, id2, id3 ! local integers266 187 !!---------------------------------------------------------------------- 267 188 ! 268 189 IF( TRIM(cdrw) == 'READ' ) THEN ! Read/initialise 269 190 IF( ln_rstart ) THEN !* Read the restart file 270 !id1 = iom_varid( numrir, 'frc_vol' , ldstop = .TRUE. )271 191 ! 272 192 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 ) 193 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst read at it= ', kt,' date= ', ndastp 194 IF(lwp) WRITE(numout,*) '~~~~~~~' 195 CALL iom_get( numrir, 'frc_voltop' , frc_voltop ) 196 CALL iom_get( numrir, 'frc_volbot' , frc_volbot ) 197 CALL iom_get( numrir, 'frc_temtop' , frc_temtop ) 198 CALL iom_get( numrir, 'frc_tembot' , frc_tembot ) 199 CALL iom_get( numrir, 'frc_sal' , frc_sal ) 200 CALL iom_get( numrir, jpdom_autoglo, 'vol_loc_ini', vol_loc_ini ) 201 CALL iom_get( numrir, jpdom_autoglo, 'tem_loc_ini', tem_loc_ini ) 202 CALL iom_get( numrir, jpdom_autoglo, 'sal_loc_ini', sal_loc_ini ) 278 203 ELSE 279 204 IF(lwp) WRITE(numout,*) '~~~~~~~' 280 205 IF(lwp) WRITE(numout,*) ' lim_diahsb at initial state ' 281 206 IF(lwp) WRITE(numout,*) '~~~~~~~' 282 frc_vol = 0._wp 283 frc_sal = 0._wp 284 bg_grme = 0._wp 207 ! set trends to 0 208 frc_voltop = 0._wp 209 frc_volbot = 0._wp 210 frc_temtop = 0._wp 211 frc_tembot = 0._wp 212 frc_sal = 0._wp 213 ! record initial ice volume, salt and temp 214 vol_loc_ini(:,:) = rhoic * vt_i(:,:) + rhosn * vt_s(:,:) ! ice/snow volume (kg/m2) 215 tem_loc_ini(:,:) = et_i(:,:) + et_s(:,:) ! ice/snow heat content (J) 216 sal_loc_ini(:,:) = rhoic * SUM( smv_i(:,:,:), dim=3 ) ! ice salt content (pss*kg/m2) 217 285 218 ENDIF 286 219 … … 288 221 ! ! ------------------- 289 222 IF(lwp) WRITE(numout,*) '~~~~~~~' 290 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst at it= ', kt,' date= ', ndastp223 IF(lwp) WRITE(numout,*) ' lim_diahsb_rst write at it= ', kt,' date= ', ndastp 291 224 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 ) 225 CALL iom_rstput( kt, nitrst, numriw, 'frc_voltop' , frc_voltop ) 226 CALL iom_rstput( kt, nitrst, numriw, 'frc_volbot' , frc_volbot ) 227 CALL iom_rstput( kt, nitrst, numriw, 'frc_temtop' , frc_temtop ) 228 CALL iom_rstput( kt, nitrst, numriw, 'frc_tembot' , frc_tembot ) 229 CALL iom_rstput( kt, nitrst, numriw, 'frc_sal' , frc_sal ) 230 CALL iom_rstput( kt, nitrst, numriw, 'vol_loc_ini', vol_loc_ini ) 231 CALL iom_rstput( kt, nitrst, numriw, 'tem_loc_ini', tem_loc_ini ) 232 CALL iom_rstput( kt, nitrst, numriw, 'sal_loc_ini', sal_loc_ini ) 295 233 ! 296 234 ENDIF -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limdyn.F90
r5836 r7355 17 17 USE phycst ! physical constants 18 18 USE dom_oce ! ocean space and time domain 19 USE sbc_oce ! Surface boundary condition: ocean fields20 19 USE sbc_ice ! Surface boundary condition: ice fields 21 20 USE ice ! LIM-3 variables 22 USE dom_ice ! LIM-3 domain23 21 USE limrhg ! LIM-3 rheology 24 22 USE lbclnk ! lateral boundary conditions - MPP exchanges … … 26 24 USE wrk_nemo ! work arrays 27 25 USE in_out_manager ! I/O manager 28 USE prtctl ! Print control29 26 USE lib_fortran ! glob_sum 30 USE timing ! Timing 31 USE limcons ! conservation tests 27 USE timing ! Timing 28 USE limcons ! conservation tests 29 USE limctl ! control prints 32 30 USE limvar 33 31 … … 35 33 PRIVATE 36 34 37 PUBLIC lim_dyn ! routine called by ice_step 35 PUBLIC lim_dyn ! routine called by sbcice_lim.F90 36 PUBLIC lim_dyn_init ! routine called by sbcice_lim.F90 38 37 39 38 !! * Substitutions … … 50 49 !! *** ROUTINE lim_dyn *** 51 50 !! 52 !! ** Purpose : compute ice velocity and ocean-ice stress51 !! ** Purpose : compute ice velocity 53 52 !! 54 53 !! ** Method : … … 56 55 !! ** Action : - Initialisation 57 56 !! - Call of the dynamic routine for each hemisphere 58 !! - computation of the stress at the ocean surface59 !! - treatment of the case if no ice dynamic60 57 !!------------------------------------------------------------------------------------ 61 58 INTEGER, INTENT(in) :: kt ! number of iteration 62 59 !! 63 INTEGER :: ji, jj, jl, ja ! dummy loop indices 64 INTEGER :: i_j1, i_jpj ! Starting/ending j-indices for rheology 65 REAL(wp) :: zcoef ! local scalar 66 REAL(wp), POINTER, DIMENSION(:) :: zswitch ! i-averaged indicator of sea-ice 67 REAL(wp), POINTER, DIMENSION(:) :: zmsk ! i-averaged of tmask 68 REAL(wp), POINTER, DIMENSION(:,:) :: zu_io, zv_io ! ice-ocean velocity 69 ! 60 INTEGER :: jl, jk ! dummy loop indices 70 61 REAL(wp) :: zvi_b, zsmv_b, zei_b, zfs_b, zfw_b, zft_b 71 62 !!--------------------------------------------------------------------- … … 73 64 IF( nn_timing == 1 ) CALL timing_start('limdyn') 74 65 75 CALL wrk_alloc( jpi, jpj, zu_io, zv_io ) 76 CALL wrk_alloc( jpj, zswitch, zmsk ) 77 78 CALL lim_var_agg(1) ! aggregate ice categories 79 80 IF( kt == nit000 ) CALL lim_dyn_init ! Initialization (first time-step only) 81 82 IF( ln_limdyn ) THEN 83 ! 84 ! conservation test 85 IF( ln_limdiahsb ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 86 87 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 88 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 89 90 ! Rheology (ice dynamics) 91 ! ======== 92 93 ! Define the j-limits where ice rheology is computed 94 ! --------------------------------------------------- 95 96 IF( lk_mpp .OR. lk_mpp_rep ) THEN ! mpp: compute over the whole domain 97 i_j1 = 1 98 i_jpj = jpj 99 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 100 CALL lim_rhg( i_j1, i_jpj ) 101 ELSE ! optimization of the computational area 102 ! 103 DO jj = 1, jpj 104 zswitch(jj) = SUM( 1.0 - at_i(:,jj) ) ! = REAL(jpj) if ocean everywhere on a j-line 105 zmsk (jj) = SUM( tmask(:,jj,1) ) ! = 0 if land everywhere on a j-line 106 END DO 107 108 IF( l_jeq ) THEN ! local domain include both hemisphere 109 ! ! Rheology is computed in each hemisphere 110 ! ! only over the ice cover latitude strip 111 ! Northern hemisphere 112 i_j1 = njeq 113 i_jpj = jpj 114 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 115 i_j1 = i_j1 + 1 116 END DO 117 i_j1 = MAX( 1, i_j1-2 ) 118 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : NH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 119 CALL lim_rhg( i_j1, i_jpj ) 120 ! 121 ! Southern hemisphere 122 i_j1 = 1 123 i_jpj = njeq 124 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 125 i_jpj = i_jpj - 1 126 END DO 127 i_jpj = MIN( jpj, i_jpj+1 ) 128 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : SH i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 129 ! 130 CALL lim_rhg( i_j1, i_jpj ) 131 ! 132 ELSE ! local domain extends over one hemisphere only 133 ! ! Rheology is computed only over the ice cover 134 ! ! latitude strip 135 i_j1 = 1 136 DO WHILE ( i_j1 <= jpj .AND. zswitch(i_j1) == FLOAT(jpi) .AND. zmsk(i_j1) /=0 ) 137 i_j1 = i_j1 + 1 138 END DO 139 i_j1 = MAX( 1, i_j1-2 ) 140 141 i_jpj = jpj 142 DO WHILE ( i_jpj >= 1 .AND. zswitch(i_jpj) == FLOAT(jpi) .AND. zmsk(i_jpj) /=0 ) 143 i_jpj = i_jpj - 1 144 END DO 145 i_jpj = MIN( jpj, i_jpj+1) 146 ! 147 IF(ln_ctl) CALL prt_ctl_info( 'lim_dyn : one hemisphere: i_j1 = ', ivar1=i_j1, clinfo2=' ij_jpj = ', ivar2=i_jpj ) 148 ! 149 CALL lim_rhg( i_j1, i_jpj ) 150 ! 151 ENDIF 152 ! 153 ENDIF 154 155 ! computation of friction velocity 156 ! -------------------------------- 157 ! ice-ocean velocity at U & V-points (u_ice v_ice at U- & V-points ; ssu_m, ssv_m at U- & V-points) 158 zu_io(:,:) = u_ice(:,:) - ssu_m(:,:) 159 zv_io(:,:) = v_ice(:,:) - ssv_m(:,:) 160 ! frictional velocity at T-point 161 zcoef = 0.5_wp * rn_cio 162 DO jj = 2, jpjm1 163 DO ji = fs_2, fs_jpim1 ! vector opt. 164 ust2s(ji,jj) = zcoef * ( zu_io(ji,jj) * zu_io(ji,jj) + zu_io(ji-1,jj) * zu_io(ji-1,jj) & 165 & + zv_io(ji,jj) * zv_io(ji,jj) + zv_io(ji,jj-1) * zv_io(ji,jj-1) ) * tmask(ji,jj,1) 166 END DO 167 END DO 168 ! 169 ! conservation test 170 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 171 ! 172 ELSE ! no ice dynamics : transmit directly the atmospheric stress to the ocean 173 ! 174 zcoef = SQRT( 0.5_wp ) * r1_rau0 175 DO jj = 2, jpjm1 176 DO ji = fs_2, fs_jpim1 ! vector opt. 177 ust2s(ji,jj) = zcoef * SQRT( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 178 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) * tmask(ji,jj,1) 179 END DO 180 END DO 181 ! 182 ENDIF 183 CALL lbc_lnk( ust2s, 'T', 1. ) ! T-point 184 185 IF(ln_ctl) THEN ! Control print 186 CALL prt_ctl_info(' ') 187 CALL prt_ctl_info(' - Cell values : ') 188 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ') 189 CALL prt_ctl(tab2d_1=ust2s , clinfo1=' lim_dyn : ust2s :') 190 CALL prt_ctl(tab2d_1=divu_i , clinfo1=' lim_dyn : divu_i :') 191 CALL prt_ctl(tab2d_1=delta_i , clinfo1=' lim_dyn : delta_i :') 192 CALL prt_ctl(tab2d_1=strength , clinfo1=' lim_dyn : strength :') 193 CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' lim_dyn : cell area :') 194 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_dyn : at_i :') 195 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_dyn : vt_i :') 196 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_dyn : vt_s :') 197 CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' lim_dyn : stress1_i :') 198 CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' lim_dyn : stress2_i :') 199 CALL prt_ctl(tab2d_1=stress12_i, clinfo1=' lim_dyn : stress12_i:') 66 CALL lim_var_agg(1) ! aggregate ice categories 67 ! 68 ! conservation test 69 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 70 71 ! ice velocities before rheology 72 u_ice_b(:,:) = u_ice(:,:) * umask(:,:,1) 73 v_ice_b(:,:) = v_ice(:,:) * vmask(:,:,1) 74 75 ! Landfast ice parameterization: define max bottom friction 76 tau_icebfr(:,:) = 0._wp 77 IF( ln_landfast ) THEN 200 78 DO jl = 1, jpl 201 CALL prt_ctl_info(' ') 202 CALL prt_ctl_info(' - Category : ', ivar1=jl) 203 CALL prt_ctl_info(' ~~~~~~~~~~') 204 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_dyn : a_i : ') 205 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_dyn : ht_i : ') 206 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_dyn : ht_s : ') 207 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_dyn : v_i : ') 208 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_dyn : v_s : ') 209 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_dyn : e_s : ') 210 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_dyn : t_su : ') 211 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_dyn : t_snow : ') 212 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_dyn : sm_i : ') 213 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_dyn : smv_i : ') 214 DO ja = 1, nlay_i 215 CALL prt_ctl_info(' ') 216 CALL prt_ctl_info(' - Layer : ', ivar1=ja) 217 CALL prt_ctl_info(' ~~~~~~~') 218 CALL prt_ctl(tab2d_1=t_i(:,:,ja,jl) , clinfo1= ' lim_dyn : t_i : ') 219 CALL prt_ctl(tab2d_1=e_i(:,:,ja,jl) , clinfo1= ' lim_dyn : e_i : ') 220 END DO 79 WHERE( ht_i(:,:,jl) > ht_n(:,:) * rn_gamma ) tau_icebfr(:,:) = tau_icebfr(:,:) + a_i(:,:,jl) * rn_icebfr 221 80 END DO 222 81 ENDIF 82 83 ! Rheology (ice dynamics) 84 ! ======== 85 CALL lim_rhg 223 86 ! 224 CALL wrk_dealloc( jpi, jpj, zu_io, zv_io ) 225 CALL wrk_dealloc( jpj, zswitch, zmsk ) 87 ! conservation test 88 IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limdyn', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 89 90 ! Control prints 91 IF( ln_ctl ) CALL lim_prt3D( 'limdyn' ) 226 92 ! 227 93 IF( nn_timing == 1 ) CALL timing_stop('limdyn') … … 243 109 !!------------------------------------------------------------------- 244 110 INTEGER :: ios ! Local integer output status for namelist read 245 NAMELIST/namicedyn/ nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 246 & nn_nevp, rn_relast, nn_ahi0, rn_ahi0_ref 247 INTEGER :: ji, jj 248 REAL(wp) :: za00, zd_max 111 NAMELIST/namicedyn/ nn_limadv, nn_limadv_ord, & 112 & nn_icestr, ln_icestr_bvf, rn_pe_rdg, rn_pstar, rn_crhg, rn_cio, rn_creepl, rn_ecc, & 113 & nn_nevp, rn_relast, ln_landfast, rn_gamma, rn_icebfr, rn_lfrelax 249 114 !!------------------------------------------------------------------- 250 115 … … 262 127 WRITE(numout,*) 'lim_dyn_init : ice parameters for ice dynamics ' 263 128 WRITE(numout,*) '~~~~~~~~~~~~' 264 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 265 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 266 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 267 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 268 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 269 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 270 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 271 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 272 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 273 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 274 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 275 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 129 ! limtrp 130 WRITE(numout,*)' choose the advection scheme (-1=Prather, 0=Ulimate-Macho) nn_limadv = ', nn_limadv 131 WRITE(numout,*)' choose the order of the scheme (if ultimate) nn_limadv_ord = ', nn_limadv_ord 132 ! limrhg 133 WRITE(numout,*)' ice strength parameterization (0=Hibler 1=Rothrock) nn_icestr = ', nn_icestr 134 WRITE(numout,*)' Including brine volume in ice strength comp. ln_icestr_bvf = ', ln_icestr_bvf 135 WRITE(numout,*)' Ratio of ridging work to PotEner change in ridging rn_pe_rdg = ', rn_pe_rdg 136 WRITE(numout,*) ' drag coefficient for oceanic stress rn_cio = ', rn_cio 137 WRITE(numout,*) ' first bulk-rheology parameter rn_pstar = ', rn_pstar 138 WRITE(numout,*) ' second bulk-rhelogy parameter rn_crhg = ', rn_crhg 139 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 140 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 141 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 142 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 143 WRITE(numout,*) ' Landfast: param (T or F) ln_landfast = ', ln_landfast 144 WRITE(numout,*) ' Landfast: fraction of ocean depth that ice must reach rn_gamma = ', rn_gamma 145 WRITE(numout,*) ' Landfast: maximum bottom stress per unit area of contact rn_icebfr = ', rn_icebfr 146 WRITE(numout,*) ' Landfast: relax time scale (s-1) to reach static friction rn_lfrelax = ', rn_lfrelax 276 147 ENDIF 277 148 ! 278 usecc2 = 1._wp / ( rn_ecc * rn_ecc )279 rhoco = rau0 * rn_cio280 !281 ! Diffusion coefficients282 SELECT CASE( nn_ahi0 )283 284 CASE( 0 )285 ahiu(:,:) = rn_ahi0_ref286 ahiv(:,:) = rn_ahi0_ref287 288 IF(lwp) WRITE(numout,*) ''289 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref'290 291 CASE( 1 )292 293 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )294 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain295 296 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2297 ! (60° = min latitude for ice cover)298 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp299 300 IF(lwp) WRITE(numout,*) ''301 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')'302 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp303 304 CASE( 2 )305 306 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) )307 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain308 309 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60° latitude in orca2310 ! (60° = min latitude for ice cover)311 DO jj = 1, jpj312 DO ji = 1, jpi313 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1)314 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1)315 END DO316 END DO317 !318 IF(lwp) WRITE(numout,*) ''319 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1'320 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max321 322 END SELECT323 324 149 END SUBROUTINE lim_dyn_init 325 150 -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limhdf.F90
r7278 r7355 7 7 !! - ! 2001-05 (G. Madec, R. Hordoir) opa norm 8 8 !! 1.0 ! 2002-08 (C. Ethe) F90, free form 9 !! 3. 0! 2015-08 (O. Tintó and M. Castrillo) added lim_hdf (multiple)9 !! 3.6 ! 2015-08 (O. Tintó and M. Castrillo) added lim_hdf (multiple) 10 10 !!---------------------------------------------------------------------- 11 11 #if defined key_lim3 … … 28 28 PRIVATE 29 29 30 PUBLIC lim_hdf ! called by lim_trp30 PUBLIC lim_hdf ! called by lim_trp 31 31 PUBLIC lim_hdf_init ! called by sbc_lim_init 32 32 33 33 LOGICAL :: linit = .TRUE. ! initialization flag (set to flase after the 1st call) 34 INTEGER :: nn_convfrq !: convergence check frequency of the Crant-Nicholson scheme35 34 REAL(wp), ALLOCATABLE, SAVE, DIMENSION(:,:) :: efact ! metric coefficient 36 35 … … 44 43 CONTAINS 45 44 46 SUBROUTINE lim_hdf( ptab , ihdf_vars , jpl , nlay_i)45 SUBROUTINE lim_hdf( ptab, ihdf_vars ) 47 46 !!------------------------------------------------------------------- 48 47 !! *** ROUTINE lim_hdf *** … … 55 54 !! ** Action : update ptab with the diffusive contribution 56 55 !!------------------------------------------------------------------- 57 INTEGER :: jpl, nlay_i, isize, ihdf_vars 58 REAL(wp), DIMENSION(:,:,:), INTENT( inout ),TARGET :: ptab ! Field on which the diffusion is applied 59 ! 60 INTEGER :: ji, jj, jk, jl , jm ! dummy loop indices 61 INTEGER :: iter, ierr ! local integers 62 REAL(wp) :: zrlxint ! local scalars 63 REAL(wp), POINTER , DIMENSION ( : ) :: zconv ! local scalars 64 REAL(wp), POINTER , DIMENSION(:,:,:) :: zrlx,zdiv0, ztab0 65 REAL(wp), POINTER , DIMENSION(:,:) :: zflu, zflv, zdiv 66 CHARACTER(lc) :: charout ! local character 67 REAL(wp), PARAMETER :: zrelax = 0.5_wp ! relaxation constant for iterative procedure 68 REAL(wp), PARAMETER :: zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 69 INTEGER , PARAMETER :: its = 100 ! Maximum number of iteration 56 INTEGER, INTENT( in ) :: ihdf_vars ! number of fields to diffuse 57 REAL(wp), DIMENSION(:,:,:), INTENT( inout ), TARGET :: ptab ! Field on which the diffusion is applied 58 ! 59 INTEGER :: ji, jj, jk, jl, jm ! dummy loop indices 60 INTEGER :: iter, ierr, isize ! local integers 61 REAL(wp) :: zrlxint 62 CHARACTER(lc) :: charout ! local character 63 REAL(wp), PARAMETER :: zrelax = 0.5_wp ! relaxation constant for iterative procedure 64 REAL(wp), PARAMETER :: zalfa = 0.5_wp ! =1.0/0.5/0.0 = implicit/Cranck-Nicholson/explicit 65 INTEGER , PARAMETER :: num_iter_max = 100 ! Maximum number of iteration 66 INTEGER , PARAMETER :: num_convfrq = 5 ! convergence check frequency of the Crant-Nicholson scheme (perf. optimization) 67 REAL(wp), POINTER, DIMENSION(:) :: zconv 68 REAL(wp), POINTER, DIMENSION(:,:,:) :: zrlx, zdiv0, ztab0 69 REAL(wp), POINTER, DIMENSION(:,:) :: zflu, zflv, zdiv 70 70 !!------------------------------------------------------------------- 71 71 TYPE(arrayptr) , ALLOCATABLE, DIMENSION(:) :: pt2d_array, zrlx_array 72 CHARACTER(len=1) , ALLOCATABLE, DIMENSION(:) :: type_array ! define the nature of ptab array grid-points 73 ! ! = T , U , V , F , W and I points 74 REAL(wp) , ALLOCATABLE, DIMENSION(:) :: psgn_array ! =-1 the sign change across the north fold boundary 75 76 !!--------------------------------------------------------------------- 77 72 CHARACTER(len=1) , ALLOCATABLE, DIMENSION(:) :: type_array ! define the nature of ptab array grid-points 73 ! ! = T , U , V , F , W and I points 74 REAL(wp) , ALLOCATABLE, DIMENSION(:) :: psgn_array ! =-1 the sign change across the north fold boundary 75 !!------------------------------------------------------------------- 76 78 77 ! !== Initialisation ==! 79 78 ! +1 open water diffusion 80 isize = jpl *(ihdf_vars+nlay_i)+179 isize = jpl * ( ihdf_vars + nlay_i ) + 1 81 80 ALLOCATE( zconv (isize) ) 82 81 ALLOCATE( pt2d_array(isize) , zrlx_array(isize) ) 83 82 ALLOCATE( type_array(isize) ) 84 83 ALLOCATE( psgn_array(isize) ) 84 85 CALL wrk_alloc( jpi,jpj, zflu, zflv, zdiv ) 86 CALL wrk_alloc( jpi,jpj,isize, zrlx, zdiv0, ztab0 ) 85 87 86 CALL wrk_alloc( jpi, jpj, isize, zrlx, zdiv0, ztab0 ) 87 CALL wrk_alloc( jpi, jpj, zflu, zflv, zdiv ) 88 89 DO jk= 1 , isize 90 pt2d_array(jk)%pt2d=>ptab(:,:,jk) 91 zrlx_array(jk)%pt2d=>zrlx(:,:,jk) 92 type_array(jk)='T' 93 psgn_array(jk)=1. 88 DO jk= 1, isize 89 pt2d_array(jk)%pt2d => ptab(:,:,jk) 90 zrlx_array(jk)%pt2d => zrlx(:,:,jk) 91 type_array(jk) = 'T' 92 psgn_array(jk) = 1. 94 93 END DO 95 94 … … 99 98 IF( lk_mpp ) CALL mpp_sum( ierr ) 100 99 IF( ierr /= 0 ) CALL ctl_stop( 'STOP', 'lim_hdf : unable to allocate arrays' ) 101 DO jj = 2, jpjm1 100 DO jj = 2, jpjm1 102 101 DO ji = fs_2 , fs_jpim1 ! vector opt. 103 102 efact(ji,jj) = ( e2u(ji,jj) + e2u(ji-1,jj) + e1v(ji,jj) + e1v(ji,jj-1) ) * r1_e1e2t(ji,jj) … … 106 105 linit = .FALSE. 107 106 ENDIF 108 ! ! Time integration parameters 109 ! 110 zflu (jpi,: ) = 0._wp 111 zflv (jpi,: ) = 0._wp 112 107 ! 108 ! Arrays initialization 109 zflu(jpi,:) = 0._wp 110 zflv(jpi,:) = 0._wp 113 111 DO jk=1 , isize 114 ztab0(:, : , jk ) = ptab(:,:,jk) ! Arrays initialization112 ztab0(:, : , jk ) = ptab(:,:,jk) 115 113 zdiv0(:, 1 , jk ) = 0._wp 116 114 zdiv0(:,jpj, jk ) = 0._wp … … 119 117 END DO 120 118 121 zconv = 1._wp !== horizontal diffusion using a Crant-Nicholson scheme ==! 122 iter = 0 123 ! 124 DO WHILE( MAXVAL(zconv(:)) > ( 2._wp * 1.e-04 ) .AND. iter <= its ) ! Sub-time step loop 119 ! !== horizontal diffusion using a Crant-Nicholson scheme ==! 120 zconv(:) = 1._wp 121 iter = 0 122 ! 123 DO WHILE( MAXVAL( zconv(:) ) > ( 2._wp * 1.e-04 ) .AND. iter <= num_iter_max ) ! Sub-time step loop 125 124 ! 126 125 iter = iter + 1 ! incrementation of the sub-time step number 127 126 ! 128 127 DO jk = 1 , isize 129 jl = ( jk-1) /( ihdf_vars+nlay_i)+1130 IF ( zconv(jk) > ( 2._wp * 1.e-04 )) THEN128 jl = ( jk - 1 ) / ( ihdf_vars + nlay_i ) + 1 129 IF ( zconv(jk) > ( 2._wp * 1.e-04 ) ) THEN 131 130 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 132 131 DO ji = 1 , fs_jpim1 ! vector opt. … … 159 158 CALL lbc_lnk_multi( zrlx_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables 160 159 ! 161 IF ( MOD( iter-1 , nn_convfrq ) == 0 ) THEN !Convergence test every nn_convfrq iterations (perf. optimization ) 162 DO jk=1,isize 160 161 IF ( MOD( iter-1 , num_convfrq ) == 0 ) THEN ! Convergence test every num_convfrq iterations (perf. optimization ) 162 DO jk = 1, isize 163 163 zconv(jk) = 0._wp ! convergence test 164 164 DO jj = 2, jpjm1 … … 175 175 END DO 176 176 ! 177 END DO ! end of sub-time step loop 178 179 ! ----------------------- 180 !!! final step (clem) !!! 177 END DO ! end of sub-time step loop 178 179 ! --- final step --- ! 181 180 DO jk = 1, isize 182 jl = ( jk-1) /( ihdf_vars+nlay_i)+1181 jl = ( jk - 1 ) / ( ihdf_vars + nlay_i ) + 1 183 182 DO jj = 1, jpjm1 ! diffusive fluxes in U- and V- direction 184 183 DO ji = 1 , fs_jpim1 ! vector opt. … … 198 197 CALL lbc_lnk_multi( pt2d_array, type_array , psgn_array , isize ) ! Multiple interchange of all the variables 199 198 200 !!! final step (clem) !!! 201 ! ----------------------- 202 199 ! 203 200 IF(ln_ctl) THEN 204 201 DO jk = 1 , isize … … 209 206 ENDIF 210 207 ! 211 CALL wrk_dealloc( jpi, jpj, isize, zrlx, zdiv0, ztab0)212 CALL wrk_dealloc( jpi, jpj, zflu, zflv, zdiv)213 208 CALL wrk_dealloc( jpi,jpj, zflu, zflv, zdiv ) 209 CALL wrk_dealloc( jpi,jpj,isize, zrlx, zdiv0, ztab0 ) 210 ! 214 211 DEALLOCATE( zconv ) 215 212 DEALLOCATE( pt2d_array , zrlx_array ) … … 219 216 END SUBROUTINE lim_hdf 220 217 221 222 218 223 219 SUBROUTINE lim_hdf_init … … 232 228 !!------------------------------------------------------------------- 233 229 INTEGER :: ios ! Local integer output status for namelist read 234 NAMELIST/namicehdf/ nn_convfrq 235 !!------------------------------------------------------------------- 236 ! 237 IF(lwp) THEN 238 WRITE(numout,*) 239 WRITE(numout,*) 'lim_hdf : Ice horizontal diffusion' 240 WRITE(numout,*) '~~~~~~~' 241 ENDIF 230 NAMELIST/namicehdf/ nn_ahi0, rn_ahi0_ref 231 INTEGER :: ji, jj 232 REAL(wp) :: za00, zd_max 233 !!------------------------------------------------------------------- 242 234 ! 243 235 REWIND( numnam_ice_ref ) ! Namelist namicehdf in reference namelist : Ice horizontal diffusion … … 252 244 IF(lwp) THEN ! control print 253 245 WRITE(numout,*) 254 WRITE(numout,*)' Namelist of ice parameters for ice horizontal diffusion computation ' 255 WRITE(numout,*)' convergence check frequency of the Crant-Nicholson scheme nn_convfrq = ', nn_convfrq 246 WRITE(numout,*) 'lim_hdf_init : Ice horizontal diffusion' 247 WRITE(numout,*) '~~~~~~~~~~~' 248 WRITE(numout,*) ' horizontal diffusivity calculation nn_ahi0 = ', nn_ahi0 249 WRITE(numout,*) ' horizontal diffusivity coeff. (orca2 grid) rn_ahi0_ref = ', rn_ahi0_ref 256 250 ENDIF 251 ! 252 ! Diffusion coefficients 253 SELECT CASE( nn_ahi0 ) 254 255 CASE( 0 ) 256 ahiu(:,:) = rn_ahi0_ref 257 ahiv(:,:) = rn_ahi0_ref 258 259 IF(lwp) WRITE(numout,*) '' 260 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim constant = rn_ahi0_ref' 261 262 CASE( 1 ) 263 264 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 265 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 266 267 ahiu(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 268 ! (60deg = min latitude for ice cover) 269 ahiv(:,:) = rn_ahi0_ref * zd_max * 1.e-05_wp 270 271 IF(lwp) WRITE(numout,*) '' 272 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to max of e1 e2 over the domain (', zd_max, ')' 273 IF(lwp) WRITE(numout,*) ' value for ahim = ', rn_ahi0_ref * zd_max * 1.e-05_wp 274 275 CASE( 2 ) 276 277 zd_max = MAX( MAXVAL( e1t(:,:) ), MAXVAL( e2t(:,:) ) ) 278 IF( lk_mpp ) CALL mpp_max( zd_max ) ! max over the global domain 279 280 za00 = rn_ahi0_ref * 1.e-05_wp ! 1.e05 = 100km = max grid space at 60deg latitude in orca2 281 ! (60deg = min latitude for ice cover) 282 DO jj = 1, jpj 283 DO ji = 1, jpi 284 ahiu(ji,jj) = za00 * MAX( e1t(ji,jj), e2t(ji,jj) ) * umask(ji,jj,1) 285 ahiv(ji,jj) = za00 * MAX( e1f(ji,jj), e2f(ji,jj) ) * vmask(ji,jj,1) 286 END DO 287 END DO 288 ! 289 IF(lwp) WRITE(numout,*) '' 290 IF(lwp) WRITE(numout,*) ' laplacian operator: ahim proportional to e1' 291 IF(lwp) WRITE(numout,*) ' maximum grid-spacing = ', zd_max, ' maximum value for ahim = ', za00*zd_max 292 293 END SELECT 257 294 ! 258 295 END SUBROUTINE lim_hdf_init … … 265 302 !!====================================================================== 266 303 END MODULE limhdf 267 -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limistate.F90
r7280 r7355 23 23 USE ice ! sea-ice variables 24 24 USE par_oce ! ocean parameters 25 USE dom_ice ! sea-ice domain26 25 USE limvar ! lim_var_salprof 27 26 ! … … 38 37 PUBLIC lim_istate ! routine called by lim_init.F90 39 38 40 ! !!** init namelist (namiceini) **41 REAL(wp) :: rn_thres_sst ! threshold water temperature for initial sea ice42 REAL(wp) :: rn_hts_ini_n ! initial snow thickness in the north43 REAL(wp) :: rn_hts_ini_s ! initial snow thickness in the south44 REAL(wp) :: rn_hti_ini_n ! initial ice thickness in the north45 REAL(wp) :: rn_hti_ini_s ! initial ice thickness in the south46 REAL(wp) :: rn_ati_ini_n ! initial leads area in the north47 REAL(wp) :: rn_ati_ini_s ! initial leads area in the south48 REAL(wp) :: rn_smi_ini_n ! initial salinity49 REAL(wp) :: rn_smi_ini_s ! initial salinity50 REAL(wp) :: rn_tmi_ini_n ! initial temperature51 REAL(wp) :: rn_tmi_ini_s ! initial temperature52 53 39 INTEGER , PARAMETER :: jpfldi = 6 ! maximum number of files to read 54 40 INTEGER , PARAMETER :: jp_hti = 1 ! index of ice thickness (m) at T-point … … 59 45 INTEGER , PARAMETER :: jp_smi = 6 ! index of ice sali at T-point 60 46 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: si ! structure of input fields (file informations, fields read) 61 62 LOGICAL :: ln_iceini ! initialization or not63 LOGICAL :: ln_iceini_file ! Ice initialization state from 2D netcdf file64 47 !!---------------------------------------------------------------------- 65 48 !! LIM 3.0, UCL-LOCEAN-IPSL (2008) … … 98 81 REAL(wp) :: ztmelts, zdh 99 82 INTEGER :: i_hemis, i_fill, jl0 100 REAL(wp) :: z test_1, ztest_2, ztest_3, ztest_4, ztests, zsigma, zarg, zA, zV, zA_cons, zV_cons, zconv83 REAL(wp) :: zarg, zV, zconv, zdv 101 84 REAL(wp), POINTER, DIMENSION(:,:) :: zswitch ! ice indicator 102 85 REAL(wp), POINTER, DIMENSION(:,:) :: zht_i_ini, zat_i_ini, zvt_i_ini !data from namelist or nc file 103 86 REAL(wp), POINTER, DIMENSION(:,:) :: zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini !data from namelist or nc file 104 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini, zv_i_ini !data by cattegories to fill 105 !-------------------------------------------------------------------- 106 107 CALL wrk_alloc( jpi,jpj,jpl, zh_i_ini, za_i_ini, zv_i_ini ) 108 CALL wrk_alloc( jpi,jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 109 CALL wrk_alloc( jpi,jpj, zswitch ) 87 REAL(wp), POINTER, DIMENSION(:,:,:) :: zh_i_ini, za_i_ini !data by cattegories to fill 88 INTEGER , POINTER, DIMENSION(:) :: itest 89 !-------------------------------------------------------------------- 90 91 CALL wrk_alloc( jpi, jpj, jpl, zh_i_ini, za_i_ini ) 92 CALL wrk_alloc( jpi, jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 93 CALL wrk_alloc( jpi, jpj, zswitch ) 94 Call wrk_alloc( 4, itest ) 110 95 111 96 IF(lwp) WRITE(numout,*) … … 117 102 !-------------------------------------------------------------------- 118 103 ! 119 CALL lim_istate_init ! reading the initials parameters of the ice120 121 ! surface temperature122 DO jl = 1, jpl ! loop over categories104 CALL lim_istate_init 105 106 ! init surface temperature 107 DO jl = 1, jpl 123 108 t_su (:,:,jl) = rt0 * tmask(:,:,1) 124 109 tn_ice(:,:,jl) = rt0 * tmask(:,:,1) 125 110 END DO 126 111 127 ! basal temperature (considered at freezing point)112 ! init basal temperature (considered at freezing point) 128 113 CALL eos_fzp( sss_m(:,:), t_bo(:,:) ) 129 114 t_bo(:,:) = ( t_bo(:,:) + rt0 ) * tmask(:,:,1) 130 115 131 116 132 IF( ln_iceini ) THEN 117 !-------------------------------------------------------------------- 118 ! 2) Initialization of sea ice state variables 119 !-------------------------------------------------------------------- 120 IF( ln_limini ) THEN 133 121 ! 134 !-------------------------------------------------------------------- 135 ! 2) Basal temperature, ice mask and hemispheric index 136 !-------------------------------------------------------------------- 122 IF( ln_limini_file )THEN 137 123 ! 138 DO jj = 1, jpj ! ice if sst <= t-freez + ttest139 DO ji = 1, jpi140 IF( ( sst_m(ji,jj) - ( t_bo(ji,jj) - rt0 ) ) * tmask(ji,jj,1) >= rn_thres_sst ) THEN141 zswitch(ji,jj) = 0._wp * tmask(ji,jj,1) ! no ice142 ELSE143 zswitch(ji,jj) = 1._wp * tmask(ji,jj,1) ! ice144 ENDIF145 END DO146 END DO147 148 !--------------------------------------------------------------------149 ! 3) Initialization of sea ice state variables150 !--------------------------------------------------------------------151 IF( ln_iceini_file )THEN152 !153 124 zht_i_ini(:,:) = si(jp_hti)%fnow(:,:,1) 154 125 zht_s_ini(:,:) = si(jp_hts)%fnow(:,:,1) … … 158 129 zsm_i_ini(:,:) = si(jp_smi)%fnow(:,:,1) 159 130 ! 160 ELSE ! ln_iceini_file = F 131 WHERE( zat_i_ini(:,:) > 0._wp ) ; zswitch(:,:) = tmask(:,:,1) 132 ELSEWHERE ; zswitch(:,:) = 0._wp 133 END WHERE 161 134 ! 135 ELSE ! ln_limini_file = F 136 137 !-------------------------------------------------------------------- 138 ! 3) Basal temperature, ice mask 139 !-------------------------------------------------------------------- 140 ! no ice if sst <= t-freez + ttest 141 WHERE( ( sst_m(:,:) - (t_bo(:,:) - rt0) ) * tmask(:,:,1) >= rn_thres_sst ) ; zswitch(:,:) = 0._wp 142 ELSEWHERE ; zswitch(:,:) = tmask(:,:,1) 143 END WHERE 144 162 145 !----------------------------- 163 146 ! 3.1) Hemisphere-dependent arrays … … 167 150 DO ji = 1, jpi 168 151 IF( ff_t(ji,jj) >= 0._wp ) THEN 169 zht_i_ini(ji,jj) = rn_hti_ini_n 170 zht_s_ini(ji,jj) = rn_hts_ini_n 171 zat_i_ini(ji,jj) = rn_ati_ini_n 172 zts_u_ini(ji,jj) = rn_tmi_ini_n 173 zsm_i_ini(ji,jj) = rn_smi_ini_n 174 ztm_i_ini(ji,jj) = rn_tmi_ini_n 152 zht_i_ini(ji,jj) = rn_hti_ini_n * zswitch(ji,jj) 153 zht_s_ini(ji,jj) = rn_hts_ini_n * zswitch(ji,jj) 154 zat_i_ini(ji,jj) = rn_ati_ini_n * zswitch(ji,jj) 155 zts_u_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 156 zsm_i_ini(ji,jj) = rn_smi_ini_n * zswitch(ji,jj) 157 ztm_i_ini(ji,jj) = rn_tmi_ini_n * zswitch(ji,jj) 175 158 ELSE 176 zht_i_ini(ji,jj) = rn_hti_ini_s 177 zht_s_ini(ji,jj) = rn_hts_ini_s 178 zat_i_ini(ji,jj) = rn_ati_ini_s 179 zts_u_ini(ji,jj) = rn_tmi_ini_s 180 zsm_i_ini(ji,jj) = rn_smi_ini_s 181 ztm_i_ini(ji,jj) = rn_tmi_ini_s 159 zht_i_ini(ji,jj) = rn_hti_ini_s * zswitch(ji,jj) 160 zht_s_ini(ji,jj) = rn_hts_ini_s * zswitch(ji,jj) 161 zat_i_ini(ji,jj) = rn_ati_ini_s * zswitch(ji,jj) 162 zts_u_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 163 zsm_i_ini(ji,jj) = rn_smi_ini_s * zswitch(ji,jj) 164 ztm_i_ini(ji,jj) = rn_tmi_ini_s * zswitch(ji,jj) 182 165 ENDIF 183 166 END DO 184 167 END DO 185 168 ! 186 ENDIF ! ln_ iceini_file187 169 ENDIF ! ln_limini_file 170 188 171 zvt_i_ini(:,:) = zht_i_ini(:,:) * zat_i_ini(:,:) ! ice volume 189 190 172 !--------------------------------------------------------------------- 191 173 ! 3.2) Distribute ice concentration and thickness into the categories … … 196 178 zh_i_ini(:,:,:) = 0._wp 197 179 za_i_ini(:,:,:) = 0._wp 198 zv_i_ini(:,:,:) = 0._wp199 180 ! 200 181 DO jj = 1, jpj … … 202 183 ! 203 184 IF( zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp )THEN 185 186 !--- jl0: most likely index where cc will be maximum 187 jl0 = jpl 188 DO jl = 1, jpl 189 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 190 jl0 = jl 191 CYCLE 192 ENDIF 193 END DO 204 194 ! 205 ztest_1 = 0 ; ztest_2 = 0 ; ztest_3 = 0 ; ztest_4 = 0 206 ! ztests = 0 207 ! 208 DO i_fill = jpl, 1, -1 195 ! initialisation of tests 196 itest(:) = 0 197 198 i_fill = jpl + 1 !==================================== 199 DO WHILE ( ( SUM( itest(:) ) /= 4 ) .AND. ( i_fill >= 2 ) ) ! iterative loop on i_fill categories 200 ! iteration !==================================== 201 i_fill = i_fill - 1 202 203 ! initialisation of ice variables for each try 204 zh_i_ini(ji,jj,:) = 0._wp 205 za_i_ini(ji,jj,:) = 0._wp 206 itest(:) = 0 209 207 ! 210 ! IF( ztests /= 4 ) THEN 211 IF ( ( ztest_1 + ztest_2 + ztest_3 + ztest_4 ) /= 4) THEN212 !----------------------------213 ! fill the i_fill categories214 !---------------------------- 215 ! *** 1 category to fill216 IF ( i_fill .EQ. 1 ) THEN217 zh_i_ini(ji,jj, 1) = zht_i_ini(ji,jj) 218 za_i_ini(ji,jj, 1) = zat_i_ini(ji,jj)219 zh_i_ini(ji,jj,2:jpl) = 0._wp220 z a_i_ini(ji,jj,2:jpl) = 0._wp221 E LSE208 ! *** case very thin ice: fill only category 1 209 IF ( i_fill == 1 ) THEN 210 zh_i_ini(ji,jj,1) = zht_i_ini(ji,jj) 211 za_i_ini(ji,jj,1) = zat_i_ini(ji,jj) 212 213 ! *** case ice is thicker: fill categories >1 214 ELSE 215 216 ! Fill ice thicknesses in the (i_fill-1) cat by hmean 217 DO jl = 1, i_fill-1 218 zh_i_ini(ji,jj,jl) = hi_mean(jl) 219 END DO 222 220 ! 223 ! *** >1 categores to fill 224 !--- Ice thicknesses in the i_fill - 1 first categories 225 DO jl = 1, i_fill - 1 226 zh_i_ini(ji,jj,jl) = hi_mean(jl) 227 END DO 221 !--- Concentrations 222 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 223 DO jl = 1, i_fill - 1 224 IF( jl /= jl0 )THEN 225 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / ( 0.5_wp * zht_i_ini(ji,jj) ) 226 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 227 ENDIF 228 END DO 228 229 ! 229 !--- jl0: most likely index where cc will be maximum 230 DO jl = 1, jpl 231 IF ( ( zht_i_ini(ji,jj) > hi_max(jl-1) ) .AND. & 232 & ( zht_i_ini(ji,jj) <= hi_max(jl) ) ) THEN 233 jl0 = jl 234 ENDIF 235 END DO 236 jl0 = MIN(jl0, i_fill) 237 ! 238 !--- Concentrations 239 za_i_ini(ji,jj,jl0) = zat_i_ini(ji,jj) / SQRT(REAL(jpl)) 240 DO jl = 1, i_fill - 1 241 IF( jl /= jl0 )THEN 242 zsigma = 0.5 * zht_i_ini(ji,jj) 243 zarg = ( zh_i_ini(ji,jj,jl) - zht_i_ini(ji,jj) ) / zsigma 244 za_i_ini(ji,jj,jl) = za_i_ini(ji,jj,jl0) * EXP(-zarg**2) 245 ENDIF 246 END DO 247 ! 248 zA = 0. ! sum of the areas in the jpl categories 249 DO jl = 1, i_fill - 1 250 zA = zA + za_i_ini(ji,jj,jl) 251 END DO 252 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - zA ! ice conc in the last category 253 IF ( i_fill < jpl ) za_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 254 ! 255 !--- Ice thickness in the last category 256 zV = 0. ! sum of the volumes of the N-1 categories 257 DO jl = 1, i_fill - 1 258 zV = zV + za_i_ini(ji,jj,jl)*zh_i_ini(ji,jj,jl) 259 END DO 260 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / za_i_ini(ji,jj,i_fill) 261 IF ( i_fill < jpl ) zh_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 262 ! 263 !--- volumes 264 zv_i_ini(ji,jj,:) = za_i_ini(ji,jj,:) * zh_i_ini(ji,jj,:) 265 IF ( i_fill < jpl ) zv_i_ini(ji,jj,i_fill+1:jpl) = 0._wp 266 ! 267 ENDIF ! i_fill 268 ! 269 !--------------------- 270 ! Compatibility tests 271 !--------------------- 272 ! Test 1: area conservation 273 zA_cons = SUM(za_i_ini(ji,jj,:)) ; zconv = ABS(zat_i_ini(ji,jj) - zA_cons ) 274 IF( zconv < 1.0e-6 ) THEN ; ztest_1 = 1. 275 ELSE ; ztest_1 = 0. 230 ! Concentration in the last (i_fill) category 231 za_i_ini(ji,jj,i_fill) = zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:i_fill-1) ) 232 233 ! Ice thickness in the last (i_fill) category 234 zV = SUM( za_i_ini(ji,jj,1:i_fill-1) * zh_i_ini(ji,jj,1:i_fill-1) ) 235 zh_i_ini(ji,jj,i_fill) = ( zvt_i_ini(ji,jj) - zV ) / MAX( za_i_ini(ji,jj,i_fill), epsi10 ) 236 237 ! clem: correction if concentration of upper cat is greater than lower cat 238 ! (it should be a gaussian around jl0 but sometimes it is not) 239 IF ( jl0 /= jpl ) THEN 240 DO jl = jpl, jl0+1, -1 241 IF ( za_i_ini(ji,jj,jl) > za_i_ini(ji,jj,jl-1) ) THEN 242 zdv = zh_i_ini(ji,jj,jl) * za_i_ini(ji,jj,jl) 243 zh_i_ini(ji,jj,jl ) = 0._wp 244 za_i_ini(ji,jj,jl ) = 0._wp 245 za_i_ini(ji,jj,1:jl-1) = za_i_ini(ji,jj,1:jl-1) & 246 & + zdv / MAX( REAL(jl-1) * zht_i_ini(ji,jj), epsi10 ) 247 END IF 248 ENDDO 276 249 ENDIF 277 250 ! 278 ! Test 2: volume conservation 279 zV_cons = SUM(zv_i_ini(ji,jj,:)) 280 zconv = ABS(zvt_i_ini(ji,jj) - zV_cons) 281 IF( zconv < 1.0e-6 ) THEN ; ztest_2 = 1. 282 ELSE ; ztest_2 = 0. 283 ENDIF 251 ENDIF ! case ice is thick or thin 252 253 !--------------------- 254 ! Compatibility tests 255 !--------------------- 256 ! Test 1: area conservation 257 zconv = ABS( zat_i_ini(ji,jj) - SUM( za_i_ini(ji,jj,1:jpl) ) ) 258 IF ( zconv < epsi06 ) itest(1) = 1 259 260 ! Test 2: volume conservation 261 zconv = ABS( zat_i_ini(ji,jj) * zht_i_ini(ji,jj) & 262 & - SUM( za_i_ini (ji,jj,1:jpl) * zh_i_ini (ji,jj,1:jpl) ) ) 263 IF ( zconv < epsi06 ) itest(2) = 1 264 265 ! Test 3: thickness of the last category is in-bounds ? 266 IF ( zh_i_ini(ji,jj,i_fill) >= hi_max(i_fill-1) ) itest(3) = 1 267 268 ! Test 4: positivity of ice concentrations 269 itest(4) = 1 270 DO jl = 1, i_fill 271 IF ( za_i_ini(ji,jj,jl) < 0._wp ) itest(4) = 0 272 END DO 273 ! !============================ 274 END DO ! end iteration on categories 275 ! !============================ 284 276 ! 285 ! Test 3: thickness of the last category is in-bounds ? 286 IF( zh_i_ini(ji,jj,i_fill) > hi_max(i_fill-1) ) THEN ; ztest_3 = 1. 287 ELSE ; ztest_3 = 0. 288 ENDIF 289 ! 290 ! Test 4: positivity of ice concentrations 291 ztest_4 = 1 292 DO jl = 1, jpl 293 IF( za_i_ini(ji,jj,jl) < 0._wp ) ztest_4 = 0 294 END DO 295 ! 296 ENDIF ! ztest_1 + ztest_2 + ztest_3 + ztest_4 297 ! 298 ztests = ztest_1 + ztest_2 + ztest_3 + ztest_4 299 ! 300 END DO ! i_fill 301 ! 302 IF(lwp) THEN 303 IF( ztests /= 4 ) THEN 304 WRITE(numout,*) 305 WRITE(numout,*) ' !!!! ALERT !!! ' 306 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 307 WRITE(numout,*) 308 WRITE(numout,*) ' *** ztests is not equal to 4 : ztests : ', ztests 309 WRITE(numout,*) ' *** ztest_i (i=1,4) = ', ztest_1, ztest_2, ztest_3, ztest_4 310 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 311 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 312 ENDIF ! ztests /= 4 277 IF( lwp .AND. SUM(itest) /= 4 ) THEN 278 WRITE(numout,*) 279 WRITE(numout,*) ' !!!! ALERT itest is not equal to 4 !!! ' 280 WRITE(numout,*) ' !!!! Something is wrong in the LIM3 initialization procedure ' 281 WRITE(numout,*) 282 WRITE(numout,*) ' *** itest_i (i=1,4) = ', itest(:) 283 WRITE(numout,*) ' zat_i_ini : ', zat_i_ini(ji,jj) 284 WRITE(numout,*) ' zht_i_ini : ', zht_i_ini(ji,jj) 313 285 ENDIF 314 !315 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zh m_i_ini(ji,jj) > 0._wp286 287 ENDIF ! zat_i_ini(ji,jj) > 0._wp .AND. zht_i_ini(ji,jj) > 0._wp 316 288 ! 317 289 END DO … … 360 332 smv_i = sm_i * v_i 361 333 ENDIF 362 334 363 335 ! Snow temperature and heat content 364 336 DO jk = 1, nlay_s … … 400 372 tn_ice (:,:,:) = t_su (:,:,:) 401 373 402 ELSE ! if ln_ iceini=false374 ELSE ! if ln_limini=false 403 375 a_i (:,:,:) = 0._wp 404 376 v_i (:,:,:) = 0._wp … … 423 395 END DO 424 396 425 ENDIF ! ln_ iceini397 ENDIF ! ln_limini 426 398 427 399 at_i (:,:) = 0.0_wp … … 473 445 sxyage (:,:,:) = 0._wp 474 446 475 476 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini, zv_i_ini ) 447 !!!clem 448 !! ! Output the initial state and forcings 449 !! CALL dia_wri_state( 'output.init', nit000 ) 450 !!! 451 452 CALL wrk_dealloc( jpi, jpj, jpl, zh_i_ini, za_i_ini ) 477 453 CALL wrk_dealloc( jpi, jpj, zht_i_ini, zat_i_ini, zvt_i_ini, zts_u_ini, zht_s_ini, zsm_i_ini, ztm_i_ini ) 478 454 CALL wrk_dealloc( jpi, jpj, zswitch ) 455 Call wrk_dealloc( 4, itest ) 479 456 480 457 END SUBROUTINE lim_istate … … 505 482 TYPE(FLD_N), DIMENSION(jpfldi) :: slf_i ! array of namelist informations on the fields to read 506 483 ! 507 NAMELIST/namiceini/ ln_ iceini, ln_iceini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, &484 NAMELIST/namiceini/ ln_limini, ln_limini_file, rn_thres_sst, rn_hts_ini_n, rn_hts_ini_s, & 508 485 & rn_hti_ini_n, rn_hti_ini_s, rn_ati_ini_n, rn_ati_ini_s, rn_smi_ini_n, & 509 486 & rn_smi_ini_s, rn_tmi_ini_n, rn_tmi_ini_s, & … … 531 508 WRITE(numout,*) 'lim_istate_init : ice parameters inititialisation ' 532 509 WRITE(numout,*) '~~~~~~~~~~~~~~~' 533 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_ iceini = ', ln_iceini534 WRITE(numout,*) ' ice initialization from a netcdf file ln_ iceini_file = ', ln_iceini_file510 WRITE(numout,*) ' initialization with ice (T) or not (F) ln_limini = ', ln_limini 511 WRITE(numout,*) ' ice initialization from a netcdf file ln_limini_file = ', ln_limini_file 535 512 WRITE(numout,*) ' threshold water temp. for initial sea-ice rn_thres_sst = ', rn_thres_sst 536 513 WRITE(numout,*) ' initial snow thickness in the north rn_hts_ini_n = ', rn_hts_ini_n … … 546 523 ENDIF 547 524 548 IF( ln_ iceini_file ) THEN ! Ice initialization using input file525 IF( ln_limini_file ) THEN ! Ice initialization using input file 549 526 ! 550 527 ! set si structure -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limitd_me.F90
r7280 r7355 18 18 USE thd_ice ! LIM thermodynamics 19 19 USE ice ! LIM variables 20 USE dom_ice ! LIM domain21 20 USE limvar ! LIM 22 21 USE lbclnk ! lateral boundary condition - MPP exchanges 23 22 USE lib_mpp ! MPP library 24 23 USE wrk_nemo ! work arrays 25 USE prtctl ! Print control26 24 27 25 USE in_out_manager ! I/O manager 28 26 USE iom ! I/O manager 29 27 USE lib_fortran ! glob_sum 30 USE limdiahsb31 28 USE timing ! Timing 32 29 USE limcons ! conservation tests 30 USE limctl ! control prints 33 31 34 32 IMPLICIT NONE … … 70 68 !! *** ROUTINE lim_itd_me_alloc *** 71 69 !!---------------------------------------------------------------------! 72 ALLOCATE( &70 ALLOCATE( & 73 71 !* Variables shared among ridging subroutines 74 & asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , & 75 & aksum(jpi,jpj) , & 76 & hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 77 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc ) 72 & asum (jpi,jpj) , athorn(jpi,jpj,0:jpl) , aksum (jpi,jpj) , & 73 & hrmin(jpi,jpj,jpl) , hraft(jpi,jpj,jpl) , aridge(jpi,jpj,jpl) , & 74 & hrmax(jpi,jpj,jpl) , krdg (jpi,jpj,jpl) , araft (jpi,jpj,jpl) , STAT=lim_itd_me_alloc ) 78 75 ! 79 76 IF( lim_itd_me_alloc /= 0 ) CALL ctl_warn( 'lim_itd_me_alloc: failed to allocate arrays' ) … … 127 124 CALL wrk_alloc( jpi,jpj, closing_net, divu_adv, opning, closing_gross ) 128 125 129 IF(ln_ctl) THEN130 CALL prt_ctl(tab2d_1=ato_i , clinfo1=' lim_itd_me: ato_i : ', tab2d_2=at_i , clinfo2=' at_i : ')131 CALL prt_ctl(tab2d_1=divu_i, clinfo1=' lim_itd_me: divu_i : ', tab2d_2=delta_i, clinfo2=' delta_i : ')132 ENDIF133 134 IF( ln_limdyn ) THEN ! Start ridging and rafting !135 136 126 ! conservation test 137 IF( ln_limdia hsb) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b)127 IF( ln_limdiachk ) CALL lim_cons_hsm(0, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 138 128 139 129 !-----------------------------------------------------------------------------! … … 211 201 DO ji = 1, jpi 212 202 za = ( opning(ji,jj) - athorn(ji,jj,0) * closing_gross(ji,jj) ) * rdt_ice 213 IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN! would lead to negative ato_i214 zfac = - ato_i(ji,jj) / za203 IF ( za < 0._wp .AND. za > - ato_i(ji,jj) ) THEN ! would lead to negative ato_i 204 zfac = - ato_i(ji,jj) / za 215 205 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) - ato_i(ji,jj) * r1_rdtice 216 206 ELSEIF( za > 0._wp .AND. za > ( asum(ji,jj) - ato_i(ji,jj) ) ) THEN ! would lead to ato_i > asum 217 zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za207 zfac = ( asum(ji,jj) - ato_i(ji,jj) ) / za 218 208 opning(ji,jj) = athorn(ji,jj,0) * closing_gross(ji,jj) + ( asum(ji,jj) - ato_i(ji,jj) ) * r1_rdtice 219 209 ENDIF … … 259 249 closing_net(ji,jj) = 0._wp 260 250 opning (ji,jj) = 0._wp 251 ato_i (ji,jj) = MAX( 0._wp, 1._wp - SUM( a_i(ji,jj,:) ) ) 261 252 ELSE 262 253 iterate_ridging = 1 … … 292 283 ! control prints 293 284 !-----------------------------------------------------------------------------! 294 IF(ln_ctl) THEN295 CALL lim_var_glo2eqv296 297 CALL prt_ctl_info(' ')298 CALL prt_ctl_info(' - Cell values : ')299 CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')300 CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' lim_itd_me : cell area :')301 CALL prt_ctl(tab2d_1=at_i , clinfo1=' lim_itd_me : at_i :')302 CALL prt_ctl(tab2d_1=vt_i , clinfo1=' lim_itd_me : vt_i :')303 CALL prt_ctl(tab2d_1=vt_s , clinfo1=' lim_itd_me : vt_s :')304 DO jl = 1, jpl305 CALL prt_ctl_info(' ')306 CALL prt_ctl_info(' - Category : ', ivar1=jl)307 CALL prt_ctl_info(' ~~~~~~~~~~')308 CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' lim_itd_me : a_i : ')309 CALL prt_ctl(tab2d_1=ht_i (:,:,jl) , clinfo1= ' lim_itd_me : ht_i : ')310 CALL prt_ctl(tab2d_1=ht_s (:,:,jl) , clinfo1= ' lim_itd_me : ht_s : ')311 CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' lim_itd_me : v_i : ')312 CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' lim_itd_me : v_s : ')313 CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' lim_itd_me : e_s : ')314 CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' lim_itd_me : t_su : ')315 CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' lim_itd_me : t_snow : ')316 CALL prt_ctl(tab2d_1=sm_i (:,:,jl) , clinfo1= ' lim_itd_me : sm_i : ')317 CALL prt_ctl(tab2d_1=smv_i (:,:,jl) , clinfo1= ' lim_itd_me : smv_i : ')318 DO jk = 1, nlay_i319 CALL prt_ctl_info(' ')320 CALL prt_ctl_info(' - Layer : ', ivar1=jk)321 CALL prt_ctl_info(' ~~~~~~~')322 CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : t_i : ')323 CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' lim_itd_me : e_i : ')324 END DO325 END DO326 ENDIF327 328 285 ! conservation test 329 IF( ln_limdiahsb ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 330 331 ENDIF ! ln_limdyn=.true. 332 ! 286 IF( ln_limdiachk ) CALL lim_cons_hsm(1, 'limitd_me', zvi_b, zsmv_b, zei_b, zfw_b, zfs_b, zft_b) 287 288 ! control prints 289 IF( ln_ctl ) CALL lim_prt3D( 'limitd_me' ) 290 333 291 CALL wrk_dealloc( jpi, jpj, closing_net, divu_adv, opning, closing_gross ) 334 292 ! … … 368 326 rswitch = MAX( 0._wp , SIGN( 1._wp, a_i(ji,jj,jl) - epsi20 ) ) 369 327 ht_i(ji,jj,jl) = v_i (ji,jj,jl) / MAX( a_i(ji,jj,jl) , epsi20 ) * rswitch 370 328 END DO 371 329 END DO 372 330 END DO … … 438 396 ENDIF 439 397 440 IF( ln_rafting ) THEN ! Ridging and rafting ice participation functions 398 ! --- Ridging and rafting participation concentrations --- ! 399 IF( ln_rafting .AND. ln_ridging ) THEN 441 400 ! 442 401 DO jl = 1, jpl … … 445 404 zdummy = TANH ( rn_craft * ( ht_i(ji,jj,jl) - rn_hraft ) ) 446 405 aridge(ji,jj,jl) = ( 1._wp + zdummy ) * 0.5_wp * athorn(ji,jj,jl) 447 araft (ji,jj,jl) = ( 1._wp - zdummy ) * 0.5_wp * athorn(ji,jj,jl)406 araft (ji,jj,jl) = athorn(ji,jj,jl) - aridge(ji,jj,jl) 448 407 END DO 449 408 END DO 450 409 END DO 451 452 ELSE 410 ! 411 ELSEIF( ln_ridging .AND. .NOT. ln_rafting ) THEN 453 412 ! 454 413 DO jl = 1, jpl 455 414 aridge(:,:,jl) = athorn(:,:,jl) 415 END DO 416 ! 417 ELSEIF( ln_rafting .AND. .NOT. ln_ridging ) THEN 418 ! 419 DO jl = 1, jpl 420 araft(:,:,jl) = athorn(:,:,jl) 456 421 END DO 457 422 ! … … 657 622 & - sm_i(ji,jj,jl1) * vsw(ij) * rhoic * r1_rdtice ! and get sm_i from the ocean 658 623 ENDIF 659 624 660 625 !------------------------------------------ 661 626 ! 3.7 Put the snow somewhere in the ocean … … 795 760 INTEGER :: numts_rm ! number of time steps for the P smoothing 796 761 REAL(wp) :: zp, z1_3 ! local scalars 797 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 762 REAL(wp), POINTER, DIMENSION(:,:) :: zworka ! temporary array used here 763 REAL(wp), POINTER, DIMENSION(:,:) :: zstrp1, zstrp2 ! strength at previous time steps 798 764 !!---------------------------------------------------------------------- 799 765 800 CALL wrk_alloc( jpi, jpj, zworka)766 CALL wrk_alloc( jpi,jpj, zworka, zstrp1, zstrp2 ) 801 767 802 768 !------------------------------------------------------------------------------! … … 844 810 END DO 845 811 846 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) 812 strength(:,:) = rn_pe_rdg * Cp * strength(:,:) / aksum(:,:) * tmask(:,:,1) 847 813 ! where Cp = (g/2)*(rhow-rhoi)*(rhoi/rhow) and rn_pe_rdg accounts for frictional dissipation 848 814 ksmooth = 1 849 815 850 851 852 816 !------------------------------------------------------------------------------! 817 ! 4) Hibler (1979)' method 818 !------------------------------------------------------------------------------! 853 819 ELSE ! kstrngth ne 1: Hibler (1979) form 854 820 ! 855 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) 821 strength(:,:) = rn_pstar * vt_i(:,:) * EXP( - rn_crhg * ( 1._wp - at_i(:,:) ) ) * tmask(:,:,1) 856 822 ! 857 823 ksmooth = 1 … … 866 832 DO jj = 1, jpj 867 833 DO ji = 1, jpi 868 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bv _i(ji,jj),0.0)))834 strength(ji,jj) = strength(ji,jj) * exp(-5.88*SQRT(MAX(bvm_i(ji,jj),0.0))) 869 835 END DO 870 836 END DO … … 880 846 IF ( ksmooth == 1 ) THEN 881 847 882 CALL lbc_lnk( strength, 'T', 1. )883 884 848 DO jj = 2, jpjm1 885 849 DO ji = 2, jpim1 886 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN850 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 887 851 zworka(ji,jj) = ( 4.0 * strength(ji,jj) & 888 852 & + strength(ji-1,jj) * tmask(ji-1,jj,1) + strength(ji+1,jj) * tmask(ji+1,jj,1) & … … 907 871 ! Temporal smoothing 908 872 !-------------------- 909 IF ( numit == nit000 + nn_fsbc - 1 ) THEN910 strp1(:,:) = 0.0911 strp2(:,:) = 0.0912 ENDIF913 914 873 IF ( ksmooth == 2 ) THEN 915 874 916 CALL lbc_lnk( strength, 'T', 1. ) 917 918 DO jj = 1, jpj - 1 919 DO ji = 1, jpi - 1 920 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp) THEN 875 IF ( numit == nit000 + nn_fsbc - 1 ) THEN 876 zstrp1(:,:) = 0._wp 877 zstrp2(:,:) = 0._wp 878 ENDIF 879 880 DO jj = 2, jpjm1 881 DO ji = 2, jpim1 882 IF ( ( asum(ji,jj) - ato_i(ji,jj) ) > 0._wp ) THEN 921 883 numts_rm = 1 ! number of time steps for the running mean 922 IF ( strp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1923 IF ( strp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1924 zp = ( strength(ji,jj) + strp1(ji,jj) +strp2(ji,jj) ) / numts_rm925 strp2(ji,jj) =strp1(ji,jj)926 strp1(ji,jj) = strength(ji,jj)884 IF ( zstrp1(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 885 IF ( zstrp2(ji,jj) > 0._wp ) numts_rm = numts_rm + 1 886 zp = ( strength(ji,jj) + zstrp1(ji,jj) + zstrp2(ji,jj) ) / numts_rm 887 zstrp2(ji,jj) = zstrp1(ji,jj) 888 zstrp1(ji,jj) = strength(ji,jj) 927 889 strength(ji,jj) = zp 928 929 890 ENDIF 930 891 END DO 931 892 END DO 932 893 894 CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions 895 933 896 ENDIF ! ksmooth 934 897 935 CALL lbc_lnk( strength, 'T', 1. ) ! Boundary conditions 936 937 CALL wrk_dealloc( jpi, jpj, zworka ) 898 CALL wrk_dealloc( jpi,jpj, zworka, zstrp1, zstrp2 ) 938 899 ! 939 900 END SUBROUTINE lim_itd_me_icestrength … … 953 914 !!------------------------------------------------------------------- 954 915 INTEGER :: ios ! Local integer output status for namelist read 955 NAMELIST/namiceitdme/ rn_cs, rn_fsnowrdg, rn_fsnowrft, & 956 & rn_gstar, rn_astar, rn_hstar, ln_rafting, rn_hraft, rn_craft, rn_por_rdg, & 957 & nn_partfun 916 NAMELIST/namiceitdme/ rn_cs, nn_partfun, rn_gstar, rn_astar, & 917 & ln_ridging, rn_hstar, rn_por_rdg, rn_fsnowrdg, ln_rafting, rn_hraft, rn_craft, rn_fsnowrft 958 918 !!------------------------------------------------------------------- 959 919 ! … … 969 929 IF (lwp) THEN ! control print 970 930 WRITE(numout,*) 971 WRITE(numout,*) 'lim_itd_me_init : ice parameters for mechanical ice redistribution '972 WRITE(numout,*) '~~~~~~~~~~~~~~~'973 WRITE(numout,*) ' Namelist namiceitdme :'974 WRITE(numout,*) ' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs975 WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg976 WRITE(numout,*) ' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft977 WRITE(numout,*) ' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar978 WRITE(numout,*) ' Equivalent to G* for an exponential part function rn_astar = ', rn_astar979 WRITE(numout,*) ' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar980 WRITE(numout,*) ' Rafting of ice sheets or not ln_rafting = ', ln_rafting981 WRITE(numout,*) ' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft982 WRITE(numout,*) ' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft983 WRITE(numout,*) ' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg984 WRITE(numout,*) ' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun931 WRITE(numout,*)'lim_itd_me_init : ice parameters for mechanical ice redistribution ' 932 WRITE(numout,*)'~~~~~~~~~~~~~~~' 933 WRITE(numout,*)' Fraction of shear energy contributing to ridging rn_cs = ', rn_cs 934 WRITE(numout,*)' Switch for part. function (0) linear (1) exponential nn_partfun = ', nn_partfun 935 WRITE(numout,*)' Fraction of total ice coverage contributing to ridging rn_gstar = ', rn_gstar 936 WRITE(numout,*)' Equivalent to G* for an exponential part function rn_astar = ', rn_astar 937 WRITE(numout,*)' Ridging of ice sheets or not ln_ridging = ', ln_ridging 938 WRITE(numout,*)' Quantity playing a role in max ridged ice thickness rn_hstar = ', rn_hstar 939 WRITE(numout,*)' Initial porosity of ridges rn_por_rdg = ', rn_por_rdg 940 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrdg = ', rn_fsnowrdg 941 WRITE(numout,*)' Rafting of ice sheets or not ln_rafting = ', ln_rafting 942 WRITE(numout,*)' Parmeter thickness (threshold between ridge-raft) rn_hraft = ', rn_hraft 943 WRITE(numout,*)' Rafting hyperbolic tangent coefficient rn_craft = ', rn_craft 944 WRITE(numout,*)' Fraction of snow volume conserved during ridging rn_fsnowrft = ', rn_fsnowrft 985 945 ENDIF 986 946 ! -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limitd_th.F90
r5407 r7355 18 18 !! lim_itd_shiftice : 19 19 !!---------------------------------------------------------------------- 20 USE dom_ice ! LIM-3 domain21 20 USE par_oce ! ocean parameters 22 21 USE dom_oce ! ocean domain -
branches/2016/dev_CNRS_2016/NEMOGCM/NEMO/LIM_SRC_3/limrhg.F90
r7278 r7355 9 9 !! 3.3 ! 2009-05 (G.Garric) addition of the lim2_evp cas 10 10 !! 3.4 ! 2011-01 (A. Porter) dynamical allocation 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 11 !! 3.5 ! 2012-08 (R. Benshila) AGRIF 12 !! 3.6 ! 2016-06 (C. Rousset) Rewriting + landfast ice + possibility to use mEVP (Bouillon 2013) 12 13 !!---------------------------------------------------------------------- 13 #if defined key_lim3 || ( defined key_lim2 && ! defined key_lim2_vp )14 #if defined key_lim3 14 15 !!---------------------------------------------------------------------- 15 !! 'key_lim3' OR LIM-3 sea-ice model 16 !! 'key_lim2' AND NOT 'key_lim2_vp' EVP LIM-2 sea-ice model 16 !! 'key_lim3' LIM-3 sea-ice model 17 17 !!---------------------------------------------------------------------- 18 18 !! lim_rhg : computes ice velocities … … 24 24 USE sbc_oce ! Surface boundary condition: ocean fields 25 25 USE sbc_ice ! Surface boundary condition: ice fields 26 #if defined key_lim3 27 USE ice ! LIM-3: ice variables 28 USE dom_ice ! LIM-3: ice domain 29 USE limitd_me ! LIM-3: 30 #else 31 USE ice_2 ! LIM-2: ice variables 32 USE dom_ice_2 ! LIM-2: ice domain 33 #endif 26 USE ice ! ice variables 27 USE limitd_me ! ice strength 34 28 USE lbclnk ! Lateral Boundary Condition / MPP link 35 29 USE lib_mpp ! MPP library … … 38 32 USE prtctl ! Print control 39 33 USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined) 40 #if defined key_agrif && defined key_lim241 USE agrif_lim 2_interp34 #if defined key_agrif 35 USE agrif_lim3_interp 42 36 #endif 43 37 #if defined key_bdy … … 48 42 PRIVATE 49 43 50 PUBLIC lim_rhg ! routine called by lim_dyn (or lim_dyn_2)44 PUBLIC lim_rhg ! routine called by lim_dyn 51 45 52 46 !! * Substitutions … … 59 53 CONTAINS 60 54 61 SUBROUTINE lim_rhg ( k_j1, k_jpj )55 SUBROUTINE lim_rhg 62 56 !!------------------------------------------------------------------- 63 57 !! *** SUBROUTINE lim_rhg *** … … 106 100 !! e.g. in the Canadian Archipelago 107 101 !! 102 !! ** Notes : There is the possibility to use mEVP from Bouillon 2013 103 !! (by uncommenting some lines in part 3 and changing alpha and beta parameters) 104 !! but this solution appears very unstable (see Kimmritz et al 2016) 105 !! 108 106 !! References : Hunke and Dukowicz, JPO97 109 107 !! Bouillon et al., Ocean Modelling 2009 108 !! Bouillon et al., Ocean Modelling 2013 110 109 !!------------------------------------------------------------------- 111 INTEGER, INTENT(in) :: k_j1 ! southern j-index for ice computation 112 INTEGER, INTENT(in) :: k_jpj ! northern j-index for ice computation 113 !! 114 INTEGER :: ji, jj ! dummy loop indices 115 INTEGER :: jter ! local integers 110 INTEGER :: ji, jj ! dummy loop indices 111 INTEGER :: jter ! local integers 116 112 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 113 114 REAL(wp) :: zrhoco ! rau0 * rn_cio 115 REAL(wp) :: zdtevp, z1_dtevp ! time step for subcycling 116 REAL(wp) :: ecc2, z1_ecc2 ! square of yield ellipse eccenticity 117 REAL(wp) :: zbeta, zalph1, z1_alph1, zalph2, z1_alph2 ! alpha and beta from Bouillon 2009 and 2013 118 REAL(wp) :: zm1, zm2, zm3, zmassU, zmassV ! ice/snow mass 119 REAL(wp) :: zdelta, zp_delf, zds2, zdt, zdt2, zdiv, zdiv2 ! temporary scalars 120 REAL(wp) :: zTauO, zTauB, zTauE, zCor, zvel ! temporary scalars 121 122 REAL(wp) :: zsig1, zsig2 ! internal ice stress 123 REAL(wp) :: zresm ! Maximal error on ice velocity 124 REAL(wp) :: zintb, zintn ! dummy argument 144 125 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 126 REAL(wp), POINTER, DIMENSION(:,:) :: z1_e1t0, z1_e2t0 ! scale factors 127 REAL(wp), POINTER, DIMENSION(:,:) :: zp_delt ! P/delta at T points 128 ! 129 REAL(wp), POINTER, DIMENSION(:,:) :: zaU , zaV ! ice fraction on U/V points 130 REAL(wp), POINTER, DIMENSION(:,:) :: zmU_t, zmV_t ! ice/snow mass/dt on U/V points 131 REAL(wp), POINTER, DIMENSION(:,:) :: zmf ! coriolis parameter at T points 132 REAL(wp), POINTER, DIMENSION(:,:) :: zTauU_ia , ztauV_ia ! ice-atm. stress at U-V points 133 REAL(wp), POINTER, DIMENSION(:,:) :: zspgU , zspgV ! surface pressure gradient at U/V points 134 REAL(wp), POINTER, DIMENSION(:,:) :: v_oceU, u_oceV, v_iceU, u_iceV ! ocean/ice u/v component on V/U points 135 REAL(wp), POINTER, DIMENSION(:,:) :: zfU , zfV ! internal stresses 136 137 REAL(wp), POINTER, DIMENSION(:,:) :: zds ! shear 138 REAL(wp), POINTER, DIMENSION(:,:) :: zs1, zs2, zs12 ! stress tensor components 139 REAL(wp), POINTER, DIMENSION(:,:) :: zu_ice, zv_ice, zresr ! check convergence 140 REAL(wp), POINTER, DIMENSION(:,:) :: zpice ! array used for the calculation of ice surface slope: 141 ! ocean surface (ssh_m) if ice is not embedded 142 ! ice top surface if ice is embedded 143 REAL(wp), POINTER, DIMENSION(:,:) :: zswitchU, zswitchV ! dummy arrays 144 REAL(wp), POINTER, DIMENSION(:,:) :: zmaskU, zmaskV ! mask for ice presence 145 REAL(wp), POINTER, DIMENSION(:,:) :: zfmask, zwf ! mask at F points for the ice 146 147 REAL(wp), PARAMETER :: zepsi = 1.0e-20_wp ! tolerance parameter 148 REAL(wp), PARAMETER :: zmmin = 1._wp ! ice mass (kg/m2) below which ice velocity equals ocean velocity 149 REAL(wp), PARAMETER :: zshlat = 2._wp ! boundary condition for sea-ice velocity (2=no slip ; 0=free slip) 156 150 !!------------------------------------------------------------------- 157 151 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 ) 162 163 #if defined key_lim2 && ! defined key_lim2_vp 164 # if defined key_agrif 165 USE ice_2, vt_s => hsnm 166 USE ice_2, vt_i => hicm 167 # else 168 vt_s => hsnm 169 vt_i => hicm 170 # endif 171 at_i(:,:) = 1. - frld(:,:) 172 #endif 173 #if defined key_agrif && defined key_lim2 174 CALL agrif_rhg_lim2_load ! First interpolation of coarse values 152 CALL wrk_alloc( jpi,jpj, z1_e1t0, z1_e2t0, zp_delt ) 153 CALL wrk_alloc( jpi,jpj, zaU, zaV, zmU_t, zmV_t, zmf, zTauU_ia, ztauV_ia ) 154 CALL wrk_alloc( jpi,jpj, zspgU, zspgV, v_oceU, u_oceV, v_iceU, u_iceV, zfU, zfV ) 155 CALL wrk_alloc( jpi,jpj, zds, zs1, zs2, zs12, zu_ice, zv_ice, zresr, zpice ) 156 CALL wrk_alloc( jpi,jpj, zswitchU, zswitchV, zmaskU, zmaskV, zfmask, zwf ) 157 158 #if defined key_agrif 159 CALL agrif_interp_lim3( 'U', 0, nn_nevp ) ! First interpolation of coarse values 160 CALL agrif_interp_lim3( 'V', 0, nn_nevp ) 175 161 #endif 176 162 ! 177 163 !------------------------------------------------------------------------------! 178 ! 1) Ice strength (zpresh) ! 179 !------------------------------------------------------------------------------! 164 ! 0) mask at F points for the ice 165 !------------------------------------------------------------------------------! 166 ! ocean/land mask 167 DO jj = 1, jpjm1 168 DO ji = 1, jpim1 ! NO vector opt. 169 zfmask(ji,jj) = tmask(ji,jj,1) * tmask(ji+1,jj,1) * tmask(ji,jj+1,1) * tmask(ji+1,jj+1,1) 170 END DO 171 END DO 172 CALL lbc_lnk( zfmask, 'F', 1._wp ) 173 174 ! Lateral boundary conditions on velocity (modify zfmask) 175 zwf(:,:) = zfmask(:,:) 176 DO jj = 2, jpjm1 177 DO ji = fs_2, fs_jpim1 ! vector opt. 178 IF( zfmask(ji,jj) == 0._wp ) THEN 179 zfmask(ji,jj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jj), zwf(ji,jj+1), zwf(ji-1,jj), zwf(ji,jj-1) ) ) 180 ENDIF 181 END DO 182 END DO 183 DO jj = 2, jpjm1 184 IF( zfmask(1,jj) == 0._wp ) THEN 185 zfmask(1 ,jj) = zshlat * MIN( 1._wp , MAX( zwf(2,jj), zwf(1,jj+1), zwf(1,jj-1) ) ) 186 ENDIF 187 IF( zfmask(jpi,jj) == 0._wp ) THEN 188 zfmask(jpi,jj) = zshlat * MIN( 1._wp , MAX( zwf(jpi,jj+1), zwf(jpim1,jj), zwf(jpi,jj-1) ) ) 189 ENDIF 190 END DO 191 DO ji = 2, jpim1 192 IF( zfmask(ji,1) == 0._wp ) THEN 193 zfmask(ji,1 ) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,1), zwf(ji,2), zwf(ji-1,1) ) ) 194 ENDIF 195 IF( zfmask(ji,jpj) == 0._wp ) THEN 196 zfmask(ji,jpj) = zshlat * MIN( 1._wp , MAX( zwf(ji+1,jpj), zwf(ji-1,jpj), zwf(ji,jpjm1) ) ) 197 ENDIF 198 END DO 199 CALL lbc_lnk( zfmask, 'F', 1._wp ) 200 201 !------------------------------------------------------------------------------! 202 ! 1) define some variables and initialize arrays 203 !------------------------------------------------------------------------------! 204 zrhoco = rau0 * rn_cio 205 206 ! ecc2: square of yield ellipse eccenticrity 207 ecc2 = rn_ecc * rn_ecc 208 z1_ecc2 = 1._wp / ecc2 209 210 ! Time step for subcycling 211 zdtevp = rdt_ice / REAL( nn_nevp ) 212 z1_dtevp = 1._wp / zdtevp 213 214 ! alpha parameters (Bouillon 2009) 215 zalph1 = ( 2._wp * rn_relast * rdt_ice ) * z1_dtevp 216 zalph2 = zalph1 * z1_ecc2 217 218 ! alpha and beta parameters (Bouillon 2013) 219 !!zalph1 = 40. 220 !!zalph2 = 40. 221 !!zbeta = 3000. 222 !!zbeta = REAL( nn_nevp ) ! close to classical EVP of Hunke (2001) 223 224 z1_alph1 = 1._wp / ( zalph1 + 1._wp ) 225 z1_alph2 = 1._wp / ( zalph2 + 1._wp ) 226 227 ! Initialise stress tensor 228 zs1 (:,:) = stress1_i (:,:) 229 zs2 (:,:) = stress2_i (:,:) 230 zs12(:,:) = stress12_i(:,:) 231 232 ! Ice strength 233 CALL lim_itd_me_icestrength( nn_icestr ) 234 235 ! scale factors 236 DO jj = 2, jpjm1 237 DO ji = fs_2, fs_jpim1 238 z1_e1t0(ji,jj) = 1._wp / ( e1t(ji+1,jj ) + e1t(ji,jj ) ) 239 z1_e2t0(ji,jj) = 1._wp / ( e2t(ji ,jj+1) + e2t(ji,jj ) ) 240 END DO 241 END DO 242 180 243 ! 181 ! Put every vector to 0182 delta_i(:,:) = 0._wp ;183 zpresh (:,:) = 0._wp ;184 zpreshc(:,:) = 0._wp185 u_ice2 (:,:) = 0._wp ; v_ice1(:,:) = 0._wp186 divu_i (:,:) = 0._wp ; zdt (:,:) = 0._wp ; zds(:,:) = 0._wp187 shear_i(:,:) = 0._wp188 189 #if defined key_lim3190 CALL lim_itd_me_icestrength( nn_icestr ) ! LIM-3: Ice strength on T-points191 #endif192 193 DO jj = k_j1 , k_jpj ! Ice mass and temp variables194 DO ji = 1 , jpi195 #if defined key_lim3196 zpresh(ji,jj) = tmask(ji,jj,1) * strength(ji,jj)197 #endif198 #if defined key_lim2199 zpresh(ji,jj) = tmask(ji,jj,1) * pstar * vt_i(ji,jj) * EXP( -c_rhg * (1. - at_i(ji,jj) ) )200 #endif201 ! zmask = 1 where there is ice or on land202 zmask(ji,jj) = 1._wp - ( 1._wp - MAX( 0._wp , SIGN ( 1._wp , vt_i(ji,jj) - zepsi ) ) ) * tmask(ji,jj,1)203 END DO204 END DO205 206 ! Ice strength on grid cell corners (zpreshc)207 ! needed for calculation of shear stress208 DO jj = k_j1+1, k_jpj-1209 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 DO216 END DO217 CALL lbc_lnk( zpreshc(:,:), 'F', 1. )218 !219 244 !------------------------------------------------------------------------------! 220 245 ! 2) Wind / ocean stress, mass terms, coriolis terms 221 246 !------------------------------------------------------------------------------! 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 247 234 248 IF( nn_ice_embd == 2 ) THEN !== embedded sea ice: compute representative ice top surface ==! … … 242 256 zintb = REAL( nn_fsbc + 1 ) / REAL( nn_fsbc ) * 0.5_wp 243 257 ! 244 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:)) * r1_rau0258 zpice(:,:) = ssh_m(:,:) + ( zintn * snwice_mass(:,:) + zintb * snwice_mass_b(:,:) ) * r1_rau0 245 259 ! 246 260 ELSE !== non-embedded sea ice: use ocean surface for slope calculation ==! … … 248 262 ENDIF 249 263 250 DO jj = k_j1+1, k_jpj-1264 DO jj = 2, jpjm1 251 265 DO ji = fs_2, fs_jpim1 252 266 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) * ff_t(ji,jj) + e1t(ji,jj) * ff_t(ji+1,jj) ) & 270 & / ( e1t(ji,jj) + e1t(ji+1,jj) + zepsi ) 271 zcorl2(ji,jj) = zmass2(ji,jj) * ( e2t(ji,jj+1) * ff_t(ji,jj) + e2t(ji,jj) * ff_t(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 267 ! ice fraction at U-V points 268 zaU(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji+1,jj) * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 269 zaV(ji,jj) = 0.5_wp * ( at_i(ji,jj) * e1e2t(ji,jj) + at_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 270 271 ! Ice/snow mass at U-V points 272 zm1 = ( rhosn * vt_s(ji ,jj ) + rhoic * vt_i(ji ,jj ) ) 273 zm2 = ( rhosn * vt_s(ji+1,jj ) + rhoic * vt_i(ji+1,jj ) ) 274 zm3 = ( rhosn * vt_s(ji ,jj+1) + rhoic * vt_i(ji ,jj+1) ) 275 zmassU = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm2 * e1e2t(ji+1,jj) ) * r1_e1e2u(ji,jj) * umask(ji,jj,1) 276 zmassV = 0.5_wp * ( zm1 * e1e2t(ji,jj) + zm3 * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 277 278 ! Ocean currents at U-V points 279 v_oceU(ji,jj) = 0.5_wp * ( ( v_oce(ji ,jj) + v_oce(ji ,jj-1) ) * e1t(ji+1,jj) & 280 & + ( v_oce(ji+1,jj) + v_oce(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 281 282 u_oceV(ji,jj) = 0.5_wp * ( ( u_oce(ji,jj ) + u_oce(ji-1,jj ) ) * e2t(ji,jj+1) & 283 & + ( u_oce(ji,jj+1) + u_oce(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 284 285 ! Coriolis at T points (m*f) 286 zmf(ji,jj) = zm1 * ff_t(ji,jj) 287 288 ! m/dt 289 zmU_t(ji,jj) = zmassU * z1_dtevp 290 zmV_t(ji,jj) = zmassV * z1_dtevp 291 292 ! Drag ice-atm. 293 zTauU_ia(ji,jj) = zaU(ji,jj) * utau_ice(ji,jj) 294 zTauV_ia(ji,jj) = zaV(ji,jj) * vtau_ice(ji,jj) 295 296 ! Surface pressure gradient (- m*g*GRAD(ssh)) at U-V points 297 zspgU(ji,jj) = - zmassU * grav * ( zpice(ji+1,jj) - zpice(ji,jj) ) * r1_e1u(ji,jj) 298 zspgV(ji,jj) = - zmassV * grav * ( zpice(ji,jj+1) - zpice(ji,jj) ) * r1_e2v(ji,jj) 299 300 ! masks 301 zmaskU(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassU ) ) ! 0 if no ice 302 zmaskV(ji,jj) = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zmassV ) ) ! 0 if no ice 303 304 ! switches 305 zswitchU(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassU - zmmin ) ) ! 0 if ice mass < zmmin 306 zswitchV(ji,jj) = MAX( 0._wp, SIGN( 1._wp, zmassV - zmmin ) ) ! 0 if ice mass < zmmin 298 307 299 308 END DO 300 309 END DO 301 310 CALL lbc_lnk( zmf, 'T', 1. ) 302 311 ! 303 312 !------------------------------------------------------------------------------! … … 305 314 !------------------------------------------------------------------------------! 306 315 ! 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 316 ! !----------------------! 326 317 DO jter = 1 , nn_nevp ! loop over jter ! 327 318 ! !----------------------! 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_e1e2t(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_e1e2t(ji,jj) 362 363 ! 319 IF(ln_ctl) THEN ! Convergence test 320 DO jj = 1, jpjm1 321 zu_ice(:,jj) = u_ice(:,jj) ! velocity at previous time step 322 zv_ice(:,jj) = v_ice(:,jj) 323 END DO 324 ENDIF 325 326 ! --- divergence, tension & shear (Appendix B of Hunke & Dukowicz, 2002) --- ! 327 DO jj = 1, jpjm1 ! loops start at 1 since there is no boundary condition (lbc_lnk) at i=1 and j=1 for F points 328 DO ji = 1, jpim1 329 330 ! shear at F points 364 331 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 332 & + ( 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_e1e2f(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 382 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_e1e2t(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_e1e2f(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_e1e2f(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. ) 333 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 334 335 END DO 336 END DO 337 CALL lbc_lnk( zds, 'F', 1. ) 338 339 DO jj = 2, jpjm1 340 DO ji = 2, jpim1 ! no vector loop 341 342 ! shear**2 at T points (doc eq. A16) 343 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 344 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 345 & ) * 0.25_wp * r1_e1e2t(ji,jj) 346 347 ! divergence at T points 348 zdiv = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 349 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 350 & ) * r1_e1e2t(ji,jj) 351 zdiv2 = zdiv * zdiv 352 353 ! tension at T points 354 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) & 355 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 356 & ) * r1_e1e2t(ji,jj) 357 zdt2 = zdt * zdt 358 359 ! delta at T points 360 zdelta = SQRT( zdiv2 + ( zdt2 + zds2 ) * z1_ecc2 ) 361 362 ! P/delta at T points 363 zp_delt(ji,jj) = strength(ji,jj) / ( zdelta + rn_creepl ) 364 365 ! stress at T points 366 zs1(ji,jj) = ( zs1(ji,jj) * zalph1 + zp_delt(ji,jj) * ( zdiv - zdelta ) ) * z1_alph1 367 zs2(ji,jj) = ( zs2(ji,jj) * zalph2 + zp_delt(ji,jj) * ( zdt * z1_ecc2 ) ) * z1_alph2 368 369 END DO 370 END DO 371 CALL lbc_lnk( zp_delt, 'T', 1. ) 372 373 DO jj = 1, jpjm1 374 DO ji = 1, jpim1 375 376 ! P/delta at F points 377 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) ) 378 379 ! stress at F points 380 zs12(ji,jj)= ( zs12(ji,jj) * zalph2 + zp_delf * ( zds(ji,jj) * z1_ecc2 ) * 0.5_wp ) * z1_alph2 381 382 END DO 383 END DO 384 CALL lbc_lnk_multi( zs1, 'T', 1., zs2, 'T', 1., zs12, 'F', 1. ) 417 385 418 ! Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) 419 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_e1e2u(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_e1e2v(ji,jj) 386 387 ! --- Ice internal stresses (Appendix C of Hunke and Dukowicz, 2002) --- ! 388 DO jj = 2, jpjm1 389 DO ji = fs_2, fs_jpim1 390 391 ! U points 392 zfU(ji,jj) = 0.5_wp * ( ( zs1(ji+1,jj) - zs1(ji,jj) ) * e2u(ji,jj) & 393 & + ( zs2(ji+1,jj) * e2t(ji+1,jj) * e2t(ji+1,jj) - zs2(ji,jj) * e2t(ji,jj) * e2t(ji,jj) & 394 & ) * r1_e2u(ji,jj) & 395 & + ( zs12(ji,jj) * e1f(ji,jj) * e1f(ji,jj) - zs12(ji,jj-1) * e1f(ji,jj-1) * e1f(ji,jj-1) & 396 & ) * 2._wp * r1_e1u(ji,jj) & 397 & ) * r1_e1e2u(ji,jj) 398 399 ! V points 400 zfV(ji,jj) = 0.5_wp * ( ( zs1(ji,jj+1) - zs1(ji,jj) ) * e1v(ji,jj) & 401 & - ( zs2(ji,jj+1) * e1t(ji,jj+1) * e1t(ji,jj+1) - zs2(ji,jj) * e1t(ji,jj) * e1t(ji,jj) & 402 & ) * r1_e1v(ji,jj) & 403 & + ( zs12(ji,jj) * e2f(ji,jj) * e2f(ji,jj) - zs12(ji-1,jj) * e2f(ji-1,jj) * e2f(ji-1,jj) & 404 & ) * 2._wp * r1_e2v(ji,jj) & 405 & ) * r1_e1e2v(ji,jj) 406 407 ! u_ice at V point 408 u_iceV(ji,jj) = 0.5_wp * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj+1) & 409 & + ( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj ) ) * z1_e2t0(ji,jj) * vmask(ji,jj,1) 410 411 ! v_ice at U point 412 v_iceU(ji,jj) = 0.5_wp * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji+1,jj) & 413 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji ,jj) ) * z1_e1t0(ji,jj) * umask(ji,jj,1) 414 431 415 END DO 432 416 END DO 433 417 ! 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 440 DO jj = k_j1+1, k_jpj-1 418 ! --- Computation of ice velocity --- ! 419 ! Bouillon et al. 2013 (eq 47-48) => unstable unless alpha, beta are chosen wisely and large nn_nevp 420 ! Bouillon et al. 2009 (eq 34-35) => stable 421 IF( MOD(jter,2) .EQ. 0 ) THEN ! even iterations 422 423 DO jj = 2, jpjm1 441 424 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 425 426 ! tau_io/(v_oce - v_ice) 427 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 428 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 429 430 ! tau_bottom/v_ice 431 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 432 zTauB = - tau_icebfr(ji,jj) / zvel 433 434 ! Coriolis at V-points (energy conserving formulation) 435 zCor = - 0.25_wp * r1_e2v(ji,jj) * & 436 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 437 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 438 439 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 440 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 441 442 ! landfast switch => 0 = static friction ; 1 = sliding friction 443 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 444 445 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 446 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 447 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 448 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 449 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 450 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 451 & ) * zmaskV(ji,jj) 452 ! Bouillon 2013 453 !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 454 !! & + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 455 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 456 455 457 END DO 456 458 END DO 457 458 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 459 #if defined key_agrif && defined key_lim2 460 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 459 CALL lbc_lnk( v_ice, 'V', -1. ) 460 461 #if defined key_agrif 462 !! CALL agrif_interp_lim3( 'V', jter, nn_nevp ) 463 CALL agrif_interp_lim3( 'V' ) 461 464 #endif 462 465 #if defined key_bdy 463 CALL bdy_ice_lim_dyn( 'U' )466 CALL bdy_ice_lim_dyn( 'V' ) 464 467 #endif 465 468 466 DO jj = k_j1+1, k_jpj-1469 DO jj = 2, jpjm1 467 470 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 471 472 ! tau_io/(u_oce - u_ice) 473 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 474 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 475 476 ! tau_bottom/u_ice 477 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 478 zTauB = - tau_icebfr(ji,jj) / zvel 479 480 ! Coriolis at U-points (energy conserving formulation) 481 zCor = 0.25_wp * r1_e1u(ji,jj) * & 482 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 483 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 484 485 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 486 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 487 488 ! landfast switch => 0 = static friction ; 1 = sliding friction 489 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 490 491 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 492 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 493 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 494 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 495 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 496 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 497 & ) * zmaskU(ji,jj) 498 ! Bouillon 2013 499 !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) & 500 !! & + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 501 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 481 502 END DO 482 503 END DO 483 484 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 485 #if defined key_agrif && defined key_lim2 486 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 504 CALL lbc_lnk( u_ice, 'U', -1. ) 505 506 #if defined key_agrif 507 !! CALL agrif_interp_lim3( 'U', jter, nn_nevp ) 508 CALL agrif_interp_lim3( 'U' ) 487 509 #endif 488 510 #if defined key_bdy 489 CALL bdy_ice_lim_dyn( 'V' )511 CALL bdy_ice_lim_dyn( 'U' ) 490 512 #endif 491 513 492 ELSE 493 DO jj = k_j1+1, k_jpj-1 514 ELSE ! odd iterations 515 516 DO jj = 2, jpjm1 494 517 DO ji = fs_2, fs_jpim1 495 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass2(ji,jj) ) ) ) * vmask(ji,jj,1) 496 z0 = zmass2(ji,jj) * z1_dtevp 497 ! SB modif because ocean has no slip boundary condition 498 zu_ice2 = 0.5 * ( ( u_ice(ji,jj ) + u_ice(ji-1,jj ) ) * e2t(ji,jj) & 499 & +( u_ice(ji,jj+1) + u_ice(ji-1,jj+1) ) * e2t(ji,jj+1) ) & 500 & / ( e2t(ji,jj+1) + e2t(ji,jj) ) * vmask(ji,jj,1) 501 502 za = rhoco * SQRT( ( zu_ice2 - u_oce2(ji,jj) )**2 + & 503 & ( v_ice(ji,jj) - v_oce(ji,jj) )**2 ) * ( 1.0 - zfrld2(ji,jj) ) 504 zr = z0 * v_ice(ji,jj) + zf2(ji,jj) + za2ct(ji,jj) + za * v_oce(ji,jj) 505 zcca = z0 + za 506 zccb = zcorl2(ji,jj) 507 v_ice(ji,jj) = ( zr - zccb * zu_ice2 ) / ( zcca + zepsi ) * rswitch 518 519 ! tau_io/(u_oce - u_ice) 520 zTauO = zaU(ji,jj) * zrhoco * SQRT( ( u_ice (ji,jj) - u_oce (ji,jj) ) * ( u_ice (ji,jj) - u_oce (ji,jj) ) & 521 & + ( v_iceU(ji,jj) - v_oceU(ji,jj) ) * ( v_iceU(ji,jj) - v_oceU(ji,jj) ) ) 522 523 ! tau_bottom/u_ice 524 zvel = MAX( zepsi, SQRT( v_iceU(ji,jj) * v_iceU(ji,jj) + u_ice(ji,jj) * u_ice(ji,jj) ) ) 525 zTauB = - tau_icebfr(ji,jj) / zvel 526 527 ! Coriolis at U-points (energy conserving formulation) 528 zCor = 0.25_wp * r1_e1u(ji,jj) * & 529 & ( zmf(ji ,jj) * ( e1v(ji ,jj) * v_ice(ji ,jj) + e1v(ji ,jj-1) * v_ice(ji ,jj-1) ) & 530 & + zmf(ji+1,jj) * ( e1v(ji+1,jj) * v_ice(ji+1,jj) + e1v(ji+1,jj-1) * v_ice(ji+1,jj-1) ) ) 531 532 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 533 zTauE = zfU(ji,jj) + zTauU_ia(ji,jj) + zCor + zspgU(ji,jj) + zTauO * ( u_oce(ji,jj) - u_ice(ji,jj) ) 534 535 ! landfast switch => 0 = static friction ; 1 = sliding friction 536 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, ztauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 537 538 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 539 u_ice(ji,jj) = ( ( rswitch * ( zmU_t(ji,jj) * u_ice(ji,jj) & ! previous velocity 540 & + zTauE + zTauO * u_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 541 & ) / MAX( zepsi, zmU_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 542 & + ( 1._wp - rswitch ) * u_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 543 & ) * zswitchU(ji,jj) + u_oce(ji,jj) * ( 1._wp - zswitchU(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 544 & ) * zmaskU(ji,jj) 545 ! Bouillon 2013 546 !!u_ice(ji,jj) = ( zmU_t(ji,jj) * ( zbeta * u_ice(ji,jj) + u_ice_b(ji,jj) ) & 547 !! & + zfU(ji,jj) + zCor + zTauU_ia(ji,jj) + zTauO * u_oce(ji,jj) + zspgU(ji,jj) & 548 !! & ) / MAX( zmU_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchU(ji,jj) 508 549 END DO 509 550 END DO 510 511 CALL lbc_lnk( v_ice(:,:), 'V', -1. ) 512 #if defined key_agrif && defined key_lim2 513 CALL agrif_rhg_lim2( jter, nn_nevp, 'V' ) 551 CALL lbc_lnk( u_ice, 'U', -1. ) 552 553 #if defined key_agrif 554 !! CALL agrif_interp_lim3( 'U', jter, nn_nevp ) 555 CALL agrif_interp_lim3( 'U' ) 514 556 #endif 515 557 #if defined key_bdy 516 CALL bdy_ice_lim_dyn( 'V' )558 CALL bdy_ice_lim_dyn( 'U' ) 517 559 #endif 518 560 519 DO jj = k_j1+1, k_jpj-1561 DO jj = 2, jpjm1 520 562 DO ji = fs_2, fs_jpim1 521 rswitch = ( 1.0 - MAX( 0._wp, SIGN( 1._wp, -zmass1(ji,jj) ) ) ) * umask(ji,jj,1) 522 z0 = zmass1(ji,jj) * z1_dtevp 523 zv_ice1 = 0.5 * ( ( v_ice(ji ,jj) + v_ice(ji ,jj-1) ) * e1t(ji,jj) & 524 & + ( v_ice(ji+1,jj) + v_ice(ji+1,jj-1) ) * e1t(ji+1,jj) ) & 525 & / ( e1t(ji+1,jj) + e1t(ji,jj) ) * umask(ji,jj,1) 526 527 za = rhoco * SQRT( ( u_ice(ji,jj) - u_oce(ji,jj) )**2 + & 528 & ( zv_ice1 - v_oce1(ji,jj) )**2 ) * ( 1.0 - zfrld1(ji,jj) ) 529 zr = z0 * u_ice(ji,jj) + zf1(ji,jj) + za1ct(ji,jj) + za * u_oce(ji,jj) 530 zcca = z0 + za 531 zccb = zcorl1(ji,jj) 532 u_ice(ji,jj) = ( zr + zccb * zv_ice1 ) / ( zcca + zepsi ) * rswitch 563 564 ! tau_io/(v_oce - v_ice) 565 zTauO = zaV(ji,jj) * zrhoco * SQRT( ( v_ice (ji,jj) - v_oce (ji,jj) ) * ( v_ice (ji,jj) - v_oce (ji,jj) ) & 566 & + ( u_iceV(ji,jj) - u_oceV(ji,jj) ) * ( u_iceV(ji,jj) - u_oceV(ji,jj) ) ) 567 568 ! tau_bottom/v_ice 569 zvel = MAX( zepsi, SQRT( v_ice(ji,jj) * v_ice(ji,jj) + u_iceV(ji,jj) * u_iceV(ji,jj) ) ) 570 ztauB = - tau_icebfr(ji,jj) / zvel 571 572 ! Coriolis at V-points (energy conserving formulation) 573 zCor = - 0.25_wp * r1_e2v(ji,jj) * & 574 & ( zmf(ji,jj ) * ( e2u(ji,jj ) * u_ice(ji,jj ) + e2u(ji-1,jj ) * u_ice(ji-1,jj ) ) & 575 & + zmf(ji,jj+1) * ( e2u(ji,jj+1) * u_ice(ji,jj+1) + e2u(ji-1,jj+1) * u_ice(ji-1,jj+1) ) ) 576 577 ! Sum of external forces (explicit solution) = F + tau_ia + Coriolis + spg + tau_io 578 zTauE = zfV(ji,jj) + zTauV_ia(ji,jj) + zCor + zspgV(ji,jj) + zTauO * ( v_oce(ji,jj) - v_ice(ji,jj) ) 579 580 ! landfast switch => 0 = static friction (tau_icebfr > zTauE); 1 = sliding friction 581 rswitch = 1._wp - MIN( 1._wp, ABS( SIGN( 1._wp, zTauE - tau_icebfr(ji,jj) ) - SIGN( 1._wp, zTauE ) ) ) 582 583 ! ice velocity using implicit formulation (cf Madec doc & Bouillon 2009) 584 v_ice(ji,jj) = ( ( rswitch * ( zmV_t(ji,jj) * v_ice(ji,jj) & ! previous velocity 585 & + zTauE + zTauO * v_ice(ji,jj) & ! F + tau_ia + Coriolis + spg + tau_io(only ocean part) 586 & ) / MAX( zepsi, zmV_t(ji,jj) + zTauO - zTauB ) & ! m/dt + tau_io(only ice part) + landfast 587 & + ( 1._wp - rswitch ) * v_ice(ji,jj) * MAX( 0._wp, 1._wp - zdtevp * rn_lfrelax ) & ! static friction => slow decrease to v=0 588 & ) * zswitchV(ji,jj) + v_oce(ji,jj) * ( 1._wp - zswitchV(ji,jj) ) & ! v_ice = v_oce if mass < zmmin 589 & ) * zmaskV(ji,jj) 590 ! Bouillon 2013 591 !!v_ice(ji,jj) = ( zmV_t(ji,jj) * ( zbeta * v_ice(ji,jj) + v_ice_b(ji,jj) ) & 592 !! & + zfV(ji,jj) + zCor + zTauV_ia(ji,jj) + zTauO * v_oce(ji,jj) + zspgV(ji,jj) & 593 !! & ) / MAX( zmV_t(ji,jj) * ( zbeta + 1._wp ) + zTauO - zTauB, zepsi ) * zswitchV(ji,jj) 533 594 END DO 534 595 END DO 535 536 CALL lbc_lnk( u_ice(:,:), 'U', -1. ) 537 #if defined key_agrif && defined key_lim2 538 CALL agrif_rhg_lim2( jter, nn_nevp, 'U' ) 596 CALL lbc_lnk( v_ice, 'V', -1. ) 597 598 #if defined key_agrif 599 !! CALL agrif_interp_lim3( 'V', jter, nn_nevp ) 600 CALL agrif_interp_lim3( 'V' ) 539 601 #endif 540 602 #if defined key_bdy 541 CALL bdy_ice_lim_dyn( 'U' )603 CALL bdy_ice_lim_dyn( 'V' ) 542 604 #endif 543 605 544 606 ENDIF 545 607 546 IF(ln_ctl) THEN 547 !--- Convergence test. 548 DO jj = k_j1+1 , k_jpj-1 608 IF(ln_ctl) THEN ! Convergence test 609 DO jj = 2 , jpjm1 549 610 zresr(:,jj) = MAX( ABS( u_ice(:,jj) - zu_ice(:,jj) ), ABS( v_ice(:,jj) - zv_ice(:,jj) ) ) 550 611 END DO 551 zresm = MAXVAL( zresr( 1:jpi, k_j1+1:k_jpj-1 ) )612 zresm = MAXVAL( zresr( 1:jpi, 2:jpjm1 ) ) 552 613 IF( lk_mpp ) CALL mpp_max( zresm ) ! max over the global domain 553 614 ENDIF 554 615 ! 555 616 ! ! ==================== ! 556 617 END DO ! end loop over jter ! … … 558 619 ! 559 620 !------------------------------------------------------------------------------! 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 621 ! 4) Recompute delta, shear and div (inputs for mechanical redistribution) 622 !------------------------------------------------------------------------------! 623 DO jj = 1, jpjm1 624 DO ji = 1, jpim1 625 626 ! shear at F points 627 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) & 628 & + ( v_ice(ji+1,jj) * r1_e2v(ji+1,jj) - v_ice(ji,jj) * r1_e2v(ji,jj) ) * e2f(ji,jj) * e2f(ji,jj) & 629 & ) * r1_e1e2f(ji,jj) * zfmask(ji,jj) 630 631 END DO 632 END DO 633 CALL lbc_lnk( zds, 'F', 1. ) 634 635 DO jj = 2, jpjm1 636 DO ji = 2, jpim1 ! no vector loop 637 638 ! tension**2 at T points 639 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) & 640 & - ( v_ice(ji,jj) * r1_e1v(ji,jj) - v_ice(ji,jj-1) * r1_e1v(ji,jj-1) ) * e1t(ji,jj) * e1t(ji,jj) & 641 & ) * r1_e1e2t(ji,jj) 642 zdt2 = zdt * zdt 643 644 ! shear**2 at T points (doc eq. A16) 645 zds2 = ( zds(ji,jj ) * zds(ji,jj ) * e1e2f(ji,jj ) + zds(ji-1,jj ) * zds(ji-1,jj ) * e1e2f(ji-1,jj ) & 646 & + zds(ji,jj-1) * zds(ji,jj-1) * e1e2f(ji,jj-1) + zds(ji-1,jj-1) * zds(ji-1,jj-1) * e1e2f(ji-1,jj-1) & 647 & ) * 0.25_wp * r1_e1e2t(ji,jj) 648 649 ! shear at T points 650 shear_i(ji,jj) = SQRT( zdt2 + zds2 ) 651 652 ! divergence at T points 653 divu_i(ji,jj) = ( e2u(ji,jj) * u_ice(ji,jj) - e2u(ji-1,jj) * u_ice(ji-1,jj) & 654 & + e1v(ji,jj) * v_ice(ji,jj) - e1v(ji,jj-1) * v_ice(ji,jj-1) & 655 & ) * r1_e1e2t(ji,jj) 656 657 ! delta at T points 658 zdelta = SQRT( divu_i(ji,jj) * divu_i(ji,jj) + ( zdt2 + zds2 ) * z1_ecc2 ) 659 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, -zdelta ) ) ! 0 if delta=0 660 delta_i(ji,jj) = zdelta + rn_creepl * rswitch 661 570 662 END DO 571 663 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_e1e2t(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_e1e2t(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_e1e2f(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_e1e2t(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_e1e2t(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 664 CALL lbc_lnk_multi( shear_i, 'T', 1., divu_i, 'T', 1., delta_i, 'T', 1. ) 665 666 ! --- Store the stress tensor for the next time step --- ! 648 667 stress1_i (:,:) = zs1 (:,:) 649 668 stress2_i (:,:) = zs2 (:,