Changeset 14016
- Timestamp:
- 2020-12-02T16:28:39+01:00 (4 years ago)
- Location:
- NEMO/branches/2020/tickets_icb_1900
- Files:
-
- 48 edited
- 5 copied
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/tickets_icb_1900
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette_ MPI3_LoopFusion@13943sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/tickets_icb_1900/cfgs/ORCA2_ICE_PISCES/EXPREF/namelist_cfg
r13899 r14016 90 90 ! ! =2 annual global mean of e-p-r set to zero 91 91 ln_wave = .false. ! Activate coupling with wave (T => fill namsbc_wave) 92 ln_cdgw = .false. ! Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)93 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)94 nn_sdrift = 0 ! Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift95 ! ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]96 ! ! = 1 Phillips: v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]97 ! ! = 2 Phillips as (1) but using the wave frequency from a wave model98 ln_tauwoc = .false. ! Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)99 ln_tauw = .false. ! Activate ocean stress components from wave model100 ln_stcor = .false. ! Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)101 92 / 102 93 !----------------------------------------------------------------------- … … 167 158 &namsbc_wave ! External fields from wave model (ln_wave=T) 168 159 !----------------------------------------------------------------------- 160 ln_sdw = .false. ! get the 2D Surf Stokes Drift & Compute the 3D stokes drift 161 ln_stcor = .false. ! add Stokes Coriolis and tracer advection terms 162 ln_cdgw = .false. ! Neutral drag coefficient read from wave model 163 ln_tauoc = .false. ! ocean stress is modified by wave induced stress 164 ln_wave_test= .false. ! Test case with constant wave fields 165 ! 166 ln_charn = .false. ! Charnock coefficient read from wave model (IFS only) 167 ln_taw = .false. ! ocean stress is modified by wave induced stress (coupled mode) 168 ln_phioc = .false. ! TKE flux from wave model 169 ln_bern_srfc= .false. ! wave induced pressure. Bernoulli head J term 170 ln_breivikFV_2016 = .false. ! breivik 2016 vertical stokes profile 171 ln_vortex_force = .false. 172 ! 173 cn_dir = './' ! root directory for the waves data location 174 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! 175 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights filename ! rotation ! land/sea mask ! 176 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! ! pairing ! filename ! 177 sn_cdg = 'sdw_ecwaves_orca2' , 6. , 'drag_coeff' , .true. , .true. , 'yearly' , '' , '' , '' 178 sn_usd = 'sdw_ecwaves_orca2' , 6. , 'u_sd2d' , .true. , .true. , 'yearly' , '' , '' , '' 179 sn_vsd = 'sdw_ecwaves_orca2' , 6. , 'v_sd2d' , .true. , .true. , 'yearly' , '' , '' , '' 180 sn_hsw = 'sdw_ecwaves_orca2' , 6. , 'hs' , .true. , .true. , 'yearly' , '' , '' , '' 181 sn_wmp = 'sdw_ecwaves_orca2' , 6. , 'wmp' , .true. , .true. , 'yearly' , '' , '' , '' 182 sn_wnum = 'sdw_ecwaves_orca2' , 6. , 'wave_num' , .true. , .true. , 'yearly' , '' , '' , '' 169 183 / 170 184 !----------------------------------------------------------------------- … … 378 392 ! = 2 add a tke source just at the base of the ML 379 393 ! = 3 as = 1 applied on HF part of the stress (ln_cpl=T) 394 ln_mxhsw = .false. ! surface mixing length scale = F(wave height) 380 395 / 381 396 !----------------------------------------------------------------------- -
NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/field_def_nemo-ice.xml
r13899 r14016 51 51 <field id="icehlid" long_name="melt pond lid depth" standard_name="sea_ice_meltpondlid_depth" unit="m" /> 52 52 <field id="icevlid" long_name="melt pond lid volume" standard_name="sea_ice_meltpondlid_volume" unit="m" /> 53 <field id="dvpn_mlt" long_name="pond volume tendency due to surface melt" standard_name="sea_ice_pondvolume_tendency_melt" unit="kg/m2/s" /> 54 <field id="dvpn_lid" long_name="pond volume tendency due to exchanges with lid" standard_name="sea_ice_pondvolume_tendency_lids" unit="kg/m2/s" /> 55 <field id="dvpn_rnf" long_name="pond volume tendency due to runoff" standard_name="sea_ice_pondvolume_tendency_runoff" unit="kg/m2/s" /> 56 <field id="dvpn_drn" long_name="pond volume tendency due to drainage" standard_name="sea_ice_pondvolume_tendency_drainage" unit="kg/m2/s" /> 53 57 54 58 <!-- heat --> … … 77 81 <field id="sig1_pnorm" long_name="P-normalized 1st principal stress component" unit="" /> 78 82 <field id="sig2_pnorm" long_name="P-normalized 2nd principal stress component" unit="" /> 83 <field id="icedlt" long_name="delta" standard_name="delta" unit="" /> 79 84 <field id="normstr" long_name="Average normal stress in sea ice" standard_name="average_normal_stress" unit="N/m" /> 80 85 <field id="sheastr" long_name="Maximum shear stress in sea ice" standard_name="maximum_shear_stress" unit="N/m" /> … … 82 87 <field id="icediv" long_name="Divergence of the sea-ice velocity field" standard_name="divergence_of_sea_ice_velocity" unit="s-1" /> 83 88 <field id="iceshe" long_name="Maximum shear of sea-ice velocity field" standard_name="maximum_shear_of_sea_ice_velocity" unit="s-1" /> 89 <field id="aniso" long_name="anisotropy of sea ice floe orientation (0.5 - 1)" standard_name="anisotropy" unit="" /> 90 <field id="yield11" long_name="yield surface tensor component 11" standard_name="yield11" unit="N/m" /> 91 <field id="yield22" long_name="yield surface tensor component 22" standard_name="yield22" unit="N/m" /> 92 <field id="yield12" long_name="yield surface tensor component 12" standard_name="yield12" unit="N/m" /> 84 93 <field id="beta_evp" long_name="Relaxation parameter of ice rheology (beta)" standard_name="relaxation_parameter_of_ice_rheology" unit="" /> 85 94 … … 297 306 <field id="snwtemp_cat" long_name="Snow temperature per category" unit="degC" detect_missing_value="true" /> 298 307 <field id="icettop_cat" long_name="Ice/snow surface temperature per category" unit="degC" detect_missing_value="true" /> 299 <field id="iceapnd_cat" long_name="Ice melt pond concentration per category" unit="" /> 308 <field id="iceapnd_cat" long_name="Ice melt pond grid fraction per category" unit="" /> 309 <field id="icevpnd_cat" long_name="Ice melt pond volume per grid area per category" unit="m" /> 300 310 <field id="icehpnd_cat" long_name="Ice melt pond thickness per category" unit="m" detect_missing_value="true" /> 301 311 <field id="icehlid_cat" long_name="Ice melt pond lid thickness per category" unit="m" detect_missing_value="true" /> 302 <field id="iceafpnd_cat" long_name="Ice melt pond fraction per category"unit="" />312 <field id="iceafpnd_cat" long_name="Ice melt pond ice fraction per category" unit="" /> 303 313 <field id="iceaepnd_cat" long_name="Ice melt pond effective fraction per category" unit="" /> 304 314 <field id="icemask_cat" long_name="Fraction of time step with sea ice (per category)" unit="" /> … … 405 415 <field field_ref="sig1_pnorm" name="sig1_pnorm"/> 406 416 <field field_ref="sig2_pnorm" name="sig2_pnorm"/> 417 <field field_ref="icedlt" name="sidelta" /> 407 418 408 419 <!-- heat fluxes --> -
NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/field_def_nemo-oce.xml
r13899 r14016 234 234 <field id="cfl_cw" long_name="w-courant number" unit="#" /> 235 235 236 <!-- variables available with ln_zdfmfc=.true. --> 237 <field id="mf_Tp" long_name="plume_temperature" standard_name="plume_temperature" unit="degC" grid_ref="grid_T_3D" /> 238 <field id="mf_Sp" long_name="plume_salinity" standard_name="plume_salinity" unit="1e-3" grid_ref="grid_T_3D" /> 239 <field id="mf_mf" long_name="mass flux" standard_name="mf_mass_flux" unit="m" grid_ref="grid_T_3D" /> 240 236 241 </field_group> <!-- grid_T --> 237 242 … … 649 654 <field id="avm_evd" long_name="convective enhancement of vertical viscosity" standard_name="ocean_vertical_momentum_diffusivity_due_to_convection" unit="m2/s" /> 650 655 656 <!-- mf_app and mf_wp: available with ln_zdfmfc --> 657 <field id="mf_app" long_name="convective area" standard_name="mf_convective_area" unit="%" grid_ref="grid_W_3D" /> 658 <field id="mf_wp" long_name="convective velocity" standard_name="mf_convective_velo" unit="m/s" grid_ref="grid_W_3D" /> 659 660 651 661 <!-- avt_tide: available with ln_zdfiwm=T --> 652 662 <field id="av_ratio" long_name="S over T diffusivity ratio" standard_name="salinity_over_temperature_diffusivity_ratio" unit="1" /> -
NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ice_ref
r13899 r14016 24 24 jpl = 5 ! number of ice categories 25 25 nlay_i = 2 ! number of ice layers 26 nlay_s = 1 ! number of snow layers (only 1 is working)26 nlay_s = 2 ! number of snow layers 27 27 ln_virtual_itd = .false. ! virtual ITD mono-category parameterization (jpl=1 only) 28 28 ! i.e. enhanced thermal conductivity & virtual thin ice melting … … 62 62 rn_lf_relax = 1.e-5 ! relaxation time scale to reach static friction [s-1] 63 63 rn_lf_tensile = 0.05 ! isotropic tensile strength [0-0.5??] 64 65 cn_dir = './' ! root directory for the grounded icebergs mask data location 66 !___________!________________!___________________!___________!_____________!________!___________!__________!__________!_______________! 67 ! ! file name ! frequency (hours) ! variable ! time interp.! clim ! 'yearly'/ ! weights ! rotation ! land/sea mask ! 68 ! ! ! (if <0 months) ! name ! (logical) ! (T/F) ! 'monthly' ! filename ! pairing ! filename ! 69 sn_icbmsk = 'NOT USED' , -12. , 'icb_mask', .false. , .true. , 'yearly' , '' , '' , '' 64 70 / 65 71 !------------------------------------------------------------------------------ … … 92 98 !------------------------------------------------------------------------------ 93 99 ln_rhg_EVP = .true. ! EVP rheology 100 ln_rhg_EAP = .false. ! EAP rheology 94 101 ln_aEVP = .true. ! adaptive rheology (Kimmritz et al. 2016 & 2017) 95 102 rn_creepl = 2.0e-9 ! creep limit [1/s] … … 98 105 rn_relast = 0.333 ! ratio of elastic timescale to ice time step: Telast = dt_ice * rn_relast 99 106 ! advised value: 1/3 (nn_nevp=100) or 1/9 (nn_nevp=300) 100 nn_rhg_chkcvg = 0 !check convergence of rheology107 nn_rhg_chkcvg = 0 ! check convergence of rheology 101 108 ! = 0 no check 102 109 ! = 1 check at the main time step (output xml: uice_cvg) 103 110 ! = 2 check at both main and rheology time steps (additional output: ice_cvg.nc) 104 111 ! this option 2 asks a lot of communications between cpu 112 ln_rhg_VP = .false. ! VP rheology 113 nn_vp_nout = 10 ! number of outer iterations 114 nn_vp_ninn = 1500 ! number of inner iterations 115 nn_vp_chkcvg = 5 ! iteration step for convergence check 105 116 / 106 117 !------------------------------------------------------------------------------ … … 195 206 !------------------------------------------------------------------------------ 196 207 ln_pnd = .true. ! activate melt ponds or not 197 ln_pnd_LEV = .true. ! level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 198 rn_apnd_min = 0.15 ! minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 199 rn_apnd_max = 0.85 ! maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? 208 ln_pnd_TOPO = .false. ! topographic melt ponds 209 ln_pnd_LEV = .true. ! level ice melt ponds 210 rn_apnd_min = 0.15 ! minimum meltwater fraction contributing to pond growth (TOPO and LEV) 211 rn_apnd_max = 0.85 ! maximum meltwater fraction contributing to pond growth (TOPO and LEV) 212 rn_pnd_flush= 0.01 ! pond flushing efficiency (tuning parameter) (LEV) 200 213 ln_pnd_CST = .false. ! constant melt ponds 201 214 rn_apnd = 0.2 ! prescribed pond fraction, at Tsu=0 degC … … 261 274 ln_icediachk = .false. ! check online heat, mass & salt budgets 262 275 ! ! rate of ice spuriously gained/lost at each time step => rn_icechk=1 <=> 1.e-6 m/hour 263 rn_icechk_cel = 1 00. ! check at each gridcell (1.e-4m/h)=> stops the code if violated (and writes a file)264 rn_icechk_glo = 1. ! check over the entire ice cover (1.e-6m/h)=> only prints warnings276 rn_icechk_cel = 1. ! check at each gridcell (1.e-06m/h)=> stops the code if violated (and writes a file) 277 rn_icechk_glo = 1.e-04 ! check over the entire ice cover (1.e-10m/h)=> only prints warnings 265 278 ln_icediahsb = .false. ! output the heat, mass & salt budgets (T) or not (F) 266 279 ln_icectl = .false. ! ice points output for debug (T or F) -
NEMO/branches/2020/tickets_icb_1900/cfgs/SHARED/namelist_ref
r14012 r14016 237 237 ln_apr_dyn = .false. ! Patm gradient added in ocean & ice Eqs. (T => fill namsbc_apr ) 238 238 ln_wave = .false. ! Activate coupling with wave (T => fill namsbc_wave) 239 ln_cdgw = .false. ! Neutral drag coefficient read from wave model (T => ln_wave=.true. & fill namsbc_wave)240 ln_sdw = .false. ! Read 2D Surf Stokes Drift & Computation of 3D stokes drift (T => ln_wave=.true. & fill namsbc_wave)241 nn_sdrift = 0 ! Parameterization for the calculation of 3D-Stokes drift from the surface Stokes drift242 ! ! = 0 Breivik 2015 parameterization: v_z=v_0*[exp(2*k*z)/(1-8*k*z)]243 ! ! = 1 Phillips: v_z=v_o*[exp(2*k*z)-beta*sqrt(-2*k*pi*z)*erfc(sqrt(-2*k*z))]244 ! ! = 2 Phillips as (1) but using the wave frequency from a wave model245 ln_tauwoc = .false. ! Activate ocean stress modified by external wave induced stress (T => ln_wave=.true. & fill namsbc_wave)246 ln_tauw = .false. ! Activate ocean stress components from wave model247 ln_stcor = .false. ! Activate Stokes Coriolis term (T => ln_wave=.true. & ln_sdw=.true. & fill namsbc_wave)248 239 nn_lsm = 0 ! =0 land/sea mask for input fields is not applied (keep empty land/sea mask filename field) , 249 240 ! =1:n number of iterations of land/sea mask application for input fields (fill land/sea mask filename field) … … 376 367 sn_rcv_cal = 'coupled' , 'no' , '' , '' , '' 377 368 sn_rcv_co2 = 'coupled' , 'no' , '' , '' , '' 378 sn_rcv_hsig = 'none' , 'no' , '' , '' , ''379 369 sn_rcv_iceflx = 'none' , 'no' , '' , '' , '' 380 370 sn_rcv_mslp = 'none' , 'no' , '' , '' , '' 381 sn_rcv_phioc = 'none' , 'no' , '' , '' , ''382 sn_rcv_sdrfx = 'none' , 'no' , '' , '' , ''383 sn_rcv_sdrfy = 'none' , 'no' , '' , '' , ''384 sn_rcv_wper = 'none' , 'no' , '' , '' , ''385 sn_rcv_wnum = 'none' , 'no' , '' , '' , ''386 sn_rcv_wfreq = 'none' , 'no' , '' , '' , ''387 sn_rcv_wdrag = 'none' , 'no' , '' , '' , ''388 371 sn_rcv_ts_ice = 'none' , 'no' , '' , '' , '' 389 372 sn_rcv_isf = 'none' , 'no' , '' , '' , '' 390 373 sn_rcv_icb = 'none' , 'no' , '' , '' , '' 391 sn_rcv_tauwoc = 'none' , 'no' , '' , '' , '' 392 sn_rcv_tauw = 'none' , 'no' , '' , '' , '' 393 sn_rcv_wdrag = 'none' , 'no' , '' , '' , '' 374 sn_rcv_hsig = 'none' , 'no' , '' ' '' , 'T' 375 sn_rcv_phioc = 'none' , 'no' , '' , '' , 'T' 376 sn_rcv_sdrfx = 'none' , 'no' , '' , '' , 'T' 377 sn_rcv_sdrfy = 'none' , 'no' , '' ' '' , 'T' 378 sn_rcv_wper = 'none' , 'no' , '' ' '' , 'T' 379 sn_rcv_wnum = 'none' , 'no' , '' ' '' , 'T' 380 sn_rcv_wstrf = 'none' , 'no' , '' ' '' , 'T' 381 sn_rcv_wdrag = 'none' , 'no' , '' ' '' , 'T' 382 sn_rcv_charn = 'none' , 'no' , '' , '' , 'T' 383 sn_rcv_taw = 'none' , 'no' , '' , '' , 'U,V' 384 sn_rcv_bhd = 'none' , 'no' , '' ' '' , 'T' 385 sn_rcv_tusd = 'none' , 'no' , '' ' '' , 'T' 386 sn_rcv_tvsd = 'none' , 'no' , '' ' '' , 'T' 394 387 / 395 388 !----------------------------------------------------------------------- … … 571 564 &namsbc_wave ! External fields from wave model (ln_wave=T) 572 565 !----------------------------------------------------------------------- 566 ln_sdw = .false. ! get the 2D Surf Stokes Drift & Compute the 3D stokes drift 567 ln_stcor = .false. ! add Stokes Coriolis and tracer advection terms 568 ln_cdgw = .false. ! Neutral drag coefficient read from wave model 569 ln_tauoc = .false. ! ocean stress is modified by wave induced stress 570 ln_wave_test= .false. ! Test case with constant wave fields 571 ! 572 ln_charn = .false. ! Charnock coefficient read from wave model (IFS only) 573 ln_taw = .false. ! ocean stress is modified by wave induced stress (coupled mode) 574 ln_phioc = .false. ! TKE flux from wave model 575 ln_bern_srfc= .false. ! wave induced pressure. Bernoulli head J term 576 ln_breivikFV_2016 = .false. ! breivik 2016 vertical stokes profile 577 ln_vortex_force = .false. ! Vortex Force term 578 ln_stshear = .false. ! include stokes shear in EKE computation 579 ! 573 580 cn_dir = './' ! root directory for the waves data location 574 581 !___________!_________________________!___________________!___________!_____________!________!___________!__________________!__________!_______________! … … 580 587 sn_hsw = 'sdw_ecwaves_orca2' , 6. , 'hs' , .true. , .true. , 'yearly' , '' , '' , '' 581 588 sn_wmp = 'sdw_ecwaves_orca2' , 6. , 'wmp' , .true. , .true. , 'yearly' , '' , '' , '' 582 sn_wfr = 'sdw_ecwaves_orca2' , 6. , 'wfr' , .true. , .true. , 'yearly' , '' , '' , ''583 589 sn_wnum = 'sdw_ecwaves_orca2' , 6. , 'wave_num' , .true. , .true. , 'yearly' , '' , '' , '' 584 sn_tauwoc = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , '' 585 sn_tauwx = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , '' 586 sn_tauwy = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , '' 590 sn_tauoc = 'sdw_ecwaves_orca2' , 6. , 'wave_stress', .true. , .true. , 'yearly' , '' , '' , '' 587 591 / 588 592 !----------------------------------------------------------------------- … … 1114 1118 nn_npc = 1 ! frequency of application of npc 1115 1119 nn_npcp = 365 ! npc control print frequency 1120 ln_zdfmfc = .false. ! Mass Flux Convection 1116 1121 ! 1117 1122 ln_zdfddm = .false. ! double diffusive mixing … … 1164 1169 rn_mxlice = 10. ! max constant ice thickness value when scaling under sea-ice ( nn_mxlice=1) 1165 1170 rn_mxl0 = 0.04 ! surface buoyancy lenght scale minimum value 1171 ln_mxhsw = .false. ! surface mixing length scale = F(wave height) 1166 1172 ln_lc = .true. ! Langmuir cell parameterisation (Axell 2002) 1167 1173 rn_lc = 0.15 ! coef. associated to Langmuir cells … … 1179 1185 ! ! = 2 weighted by 1-fr_i 1180 1186 ! ! = 3 weighted by 1-MIN(1,4*fr_i) 1187 nn_bc_surf = 1 ! surface condition (0/1=Dir/Neum) ! Only applicable for wave coupling (ln_cplwave=1) 1188 nn_bc_bot = 1 ! bottom condition (0/1=Dir/Neum) ! Only applicable for wave coupling (ln_cplwave=1) 1181 1189 / 1182 1190 !----------------------------------------------------------------------- … … 1223 1231 ! ! = 1: Pierson Moskowitz wave spectrum 1224 1232 ! ! = 0: Constant La# = 0.3 1233 / 1234 !----------------------------------------------------------------------- 1235 &namzdf_mfc ! Mass Flux Convection 1236 !----------------------------------------------------------------------- 1237 ln_edmfuv = .false. ! Activate on velocity fields (Not available yet) 1238 rn_cemf = 1. ! entrain/detrain coef. (<0 => cte; >0 % depending on dW/dz 1239 rn_cwmf = -0. ! entrain/detrain coef. (<0 => cte; >0 % depending on dW/dz 1240 rn_cent = 2.e-5 ! entrain of convective area 1241 rn_cdet = 3.e-5 ! detrain of convective area 1242 rn_cap = 0.9 ! Coef. for CAP estimation 1243 App_max = 0.1 ! Maximum convection area (% of the cell) 1225 1244 / 1226 1245 !----------------------------------------------------------------------- -
NEMO/branches/2020/tickets_icb_1900/doc/namelists/namdyn_rhg
r13899 r14016 3 3 !------------------------------------------------------------------------------ 4 4 ln_rhg_EVP = .true. ! EVP rheology 5 ln_rhg_EAP = .false. ! EAP rheology 5 6 ln_aEVP = .false. ! adaptive rheology (Kimmritz et al. 2016 & 2017) 6 7 rn_creepl = 2.0e-9 ! creep limit [1/s] -
NEMO/branches/2020/tickets_icb_1900/src/ICE/ice.F90
r13899 r14016 150 150 ! 151 151 ! !!** ice-rheology namelist (namdyn_rhg) ** 152 ! -- evp 153 LOGICAL , PUBLIC :: ln_rhg_EVP ! EVP rheology switch, used for rdgrft and rheology 154 LOGICAL , PUBLIC :: ln_rhg_EAP ! EAP rheology switch, used for rdgrft and rheology 152 155 LOGICAL , PUBLIC :: ln_aEVP !: using adaptive EVP (T or F) 153 REAL(wp), PUBLIC :: rn_creepl !: creep limit : has to be under 1.0e-9156 REAL(wp), PUBLIC :: rn_creepl !: creep limit (has to be low enough, circa 10-9 m/s, depending on rheology) 154 157 REAL(wp), PUBLIC :: rn_ecc !: eccentricity of the elliptical yield curve 155 158 INTEGER , PUBLIC :: nn_nevp !: number of iterations for subcycling 156 159 REAL(wp), PUBLIC :: rn_relast !: ratio => telast/rDt_ice (1/3 or 1/9 depending on nb of subcycling nevp) 157 160 INTEGER , PUBLIC :: nn_rhg_chkcvg !: check ice rheology convergence 161 ! -- vp 162 LOGICAL , PUBLIC :: ln_rhg_VP !: VP rheology 163 INTEGER , PUBLIC :: nn_vp_nout !: Number of outer iterations 164 INTEGER , PUBLIC :: nn_vp_ninn !: Number of inner iterations (linear system solver) 165 INTEGER , PUBLIC :: nn_vp_chkcvg !: Number of iterations every each convergence is checked 158 166 ! 159 167 ! !!** ice-advection namelist (namdyn_adv) ** … … 208 216 ! !!** ice-ponds namelist (namthd_pnd) 209 217 LOGICAL , PUBLIC :: ln_pnd !: Melt ponds (T) or not (F) 210 LOGICAL , PUBLIC :: ln_pnd_LEV !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010) 211 REAL(wp), PUBLIC :: rn_apnd_min !: Minimum ice fraction that contributes to melt ponds 212 REAL(wp), PUBLIC :: rn_apnd_max !: Maximum ice fraction that contributes to melt ponds 218 LOGICAL , PUBLIC :: ln_pnd_TOPO !: Topographic Melt ponds scheme (Flocco et al 2007, 2010) 219 LOGICAL , PUBLIC :: ln_pnd_LEV !: Simple melt pond scheme 220 REAL(wp), PUBLIC :: rn_apnd_min !: Minimum fraction of melt water contributing to ponds 221 REAL(wp), PUBLIC :: rn_apnd_max !: Maximum fraction of melt water contributing to ponds 222 REAL(wp), PUBLIC :: rn_pnd_flush !: Pond flushing efficiency (tuning parameter) 213 223 LOGICAL , PUBLIC :: ln_pnd_CST !: Melt ponds scheme with constant fraction and depth 214 224 REAL(wp), PUBLIC :: rn_apnd !: prescribed pond fraction (0<rn_apnd<1) … … 246 256 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: divu_i !: Divergence of the velocity field [s-1] 247 257 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: shear_i !: Shear of the velocity field [s-1] 258 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: aniso_11, aniso_12 !: structure tensor elements 259 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: rdg_conv 248 260 ! 249 261 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: t_bo !: Sea-Ice bottom temperature [Kelvin] … … 341 353 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: om_i !: mean ice age over all categories (s) 342 354 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: tau_icebfr !: ice friction on ocean bottom (landfast param activated) 355 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: icb_mask !: mask of grounded icebergs if landfast [0-1] 343 356 344 357 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: t_s !: Snow temperatures [K] … … 362 375 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: vt_il !: total melt pond lid volume per gridcell area [m] 363 376 377 ! meltwater arrays to save for melt ponds (mv - could be grouped in a single meltwater volume array) 378 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_i_sum_2d !: surface melt (2d arrays for ponds) [m] 379 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_s_mlt_2d !: snow surf melt (2d arrays for ponds) [m] 380 364 381 !!---------------------------------------------------------------------- 365 382 !! * Global variables at before time step 366 383 !!---------------------------------------------------------------------- 367 384 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_s_b, v_i_b, h_s_b, h_i_b !: snow and ice volumes/thickness 385 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: v_ip_b, v_il_b !: ponds and lids volumes 368 386 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: a_i_b, sv_i_b !: 369 387 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:,:) :: e_s_b !: snow heat content … … 392 410 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vsnw !: snw volume variation [m/s] 393 411 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_aice !: ice conc. variation [s-1] 412 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_vpnd !: pond volume variation [m/s] 394 413 ! 395 414 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:) :: diag_adv_mass !: advection of mass (kg/m2/s) … … 436 455 ALLOCATE( u_oce (jpi,jpj) , v_oce (jpi,jpj) , ht_i_new (jpi,jpj) , strength(jpi,jpj) , & 437 456 & stress1_i(jpi,jpj) , stress2_i(jpi,jpj) , stress12_i(jpi,jpj) , & 438 & delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , STAT=ierr(ii) ) 457 & delta_i (jpi,jpj) , divu_i (jpi,jpj) , shear_i (jpi,jpj) , & 458 & aniso_11 (jpi,jpj) , aniso_12 (jpi,jpj) , rdg_conv (jpi,jpj) , STAT=ierr(ii) ) 439 459 440 460 ii = ii + 1 … … 468 488 & et_i (jpi,jpj) , et_s (jpi,jpj) , tm_i(jpi,jpj) , tm_s(jpi,jpj) , & 469 489 & sm_i (jpi,jpj) , tm_su(jpi,jpj) , hm_i(jpi,jpj) , hm_s(jpi,jpj) , & 470 & om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj) 490 & om_i (jpi,jpj) , bvm_i(jpi,jpj) , tau_icebfr(jpi,jpj), icb_mask(jpi,jpj), STAT=ierr(ii) ) 471 491 472 492 ii = ii + 1 … … 478 498 ii = ii + 1 479 499 ALLOCATE( a_ip(jpi,jpj,jpl) , v_ip(jpi,jpj,jpl) , a_ip_frac(jpi,jpj,jpl) , h_ip(jpi,jpj,jpl), & 480 & v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , STAT = ierr(ii) ) 500 & v_il(jpi,jpj,jpl) , h_il(jpi,jpj,jpl) , a_ip_eff (jpi,jpj,jpl) , & 501 & dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , STAT = ierr(ii) ) 481 502 482 503 ii = ii + 1 … … 486 507 ii = ii + 1 487 508 ALLOCATE( v_s_b (jpi,jpj,jpl) , v_i_b (jpi,jpj,jpl) , h_s_b(jpi,jpj,jpl) , h_i_b(jpi,jpj,jpl), & 509 & v_ip_b(jpi,jpj,jpl) , v_il_b(jpi,jpj,jpl) , & 488 510 & a_i_b (jpi,jpj,jpl) , sv_i_b(jpi,jpj,jpl) , e_i_b(jpi,jpj,nlay_i,jpl) , e_s_b(jpi,jpj,nlay_s,jpl) , & 489 511 & STAT=ierr(ii) ) … … 500 522 ALLOCATE( diag_trp_vi(jpi,jpj) , diag_trp_vs (jpi,jpj) , diag_trp_ei(jpi,jpj), & 501 523 & diag_trp_es(jpi,jpj) , diag_trp_sv (jpi,jpj) , diag_heat (jpi,jpj), & 502 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), &524 & diag_sice (jpi,jpj) , diag_vice (jpi,jpj) , diag_vsnw (jpi,jpj), diag_aice(jpi,jpj), diag_vpnd(jpi,jpj), & 503 525 & diag_adv_mass(jpi,jpj), diag_adv_salt(jpi,jpj), diag_adv_heat(jpi,jpj), STAT=ierr(ii) ) 504 526 -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icectl.F90
r13899 r14016 85 85 !! 86 86 REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, & 87 & zdiag_vmin, zdiag_amin, zdiag_amax, zdiag_eimin, zdiag_esmin, zdiag_smin 87 & zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, & 88 & zdiag_eimin, zdiag_esmin, zdiag_simin 88 89 REAL(wp) :: zvtrp, zetrp 89 90 REAL(wp) :: zarea … … 92 93 IF( icount == 0 ) THEN 93 94 94 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos , dim=3 ) * e1e2t )95 pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) 95 96 pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) 96 97 pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t ) … … 112 113 113 114 ! -- mass diag -- ! 114 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos, dim=3 ) * e1e2t ) - pdiag_v ) * r1_Dt_ice & 115 zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) & 116 & - pdiag_v ) * r1_Dt_ice & 115 117 & + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + & 116 118 & wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + & … … 132 134 133 135 ! -- min/max diag -- ! 134 zdiag_amax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 135 zdiag_vmin = glob_min( 'icectl', v_i ) 136 zdiag_amin = glob_min( 'icectl', a_i ) 137 zdiag_smin = glob_min( 'icectl', sv_i ) 136 zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) ) 137 zdiag_vimin = glob_min( 'icectl', v_i ) 138 zdiag_vsmin = glob_min( 'icectl', v_s ) 139 zdiag_vpmin = glob_min( 'icectl', v_ip ) 140 zdiag_vlmin = glob_min( 'icectl', v_il ) 141 zdiag_aimin = glob_min( 'icectl', a_i ) 142 zdiag_simin = glob_min( 'icectl', sv_i ) 138 143 zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) ) 139 144 zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) ) … … 155 160 & WRITE(numout,*) cd_routine,' : violation heat cons. [J] = ',zdiag_heat * rDt_ice 156 161 ! check negative values 157 IF( zdiag_vmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vmin 158 IF( zdiag_amin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_amin 159 IF( zdiag_smin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_smin 160 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 161 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 IF( zdiag_vimin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vimin 163 IF( zdiag_vsmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_s < 0 = ',zdiag_vsmin 164 IF( zdiag_vpmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_ip < 0 = ',zdiag_vpmin 165 IF( zdiag_vlmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_il < 0 = ',zdiag_vlmin 166 IF( zdiag_aimin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_aimin 167 IF( zdiag_simin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_simin 168 IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin 169 IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin 162 170 ! check maximum ice concentration 163 IF( zdiag_a max >MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &164 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_a max171 IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) & 172 & WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_aimax 165 173 ! check if advection scheme is conservative 166 174 IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 167 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * r dt_ice175 & WRITE(numout,*) cd_routine,' : violation adv scheme [kg] = ',zvtrp * rDt_ice 168 176 IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) & 169 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * r dt_ice177 & WRITE(numout,*) cd_routine,' : violation adv scheme [J] = ',zetrp * rDt_ice 170 178 ENDIF 171 179 ! … … 193 201 ! water flux 194 202 ! -- mass diag -- ! 195 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub&196 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t )203 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd & 204 & + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) 197 205 198 206 ! -- salt diag -- ! … … 200 208 201 209 ! -- heat diag -- ! 202 zdiag_heat 210 zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t ) 203 211 ! equivalent to this: 204 212 !!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw & … … 245 253 IF( icount == 0 ) THEN 246 254 247 pdiag_v = SUM( v_i * rhoi + v_s * rhos , dim=3 )248 pdiag_s = SUM( sv_i * rhoi 255 pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) 256 pdiag_s = SUM( sv_i * rhoi , dim=3 ) 249 257 pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) 250 258 … … 261 269 262 270 ! -- mass diag -- ! 263 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos , dim=3 ) - pdiag_v ) * r1_Dt_ice&271 zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_Dt_ice & 264 272 & + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + & 265 273 & wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) & … … 352 360 CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) ! 353 361 CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) ! 362 ! mean state 363 CALL iom_rstput( 0, 0, inum, 'icecon' , SUM(a_i ,dim=3) , ktype = jp_r8 ) ! 364 CALL iom_rstput( 0, 0, inum, 'icevol' , SUM(v_i ,dim=3) , ktype = jp_r8 ) ! 365 CALL iom_rstput( 0, 0, inum, 'snwvol' , SUM(v_s ,dim=3) , ktype = jp_r8 ) ! 366 CALL iom_rstput( 0, 0, inum, 'pndvol' , SUM(v_ip,dim=3) , ktype = jp_r8 ) ! 367 CALL iom_rstput( 0, 0, inum, 'lidvol' , SUM(v_il,dim=3) , ktype = jp_r8 ) ! 354 368 355 369 CALL iom_close( inum ) … … 776 790 ! -- mass diag -- ! 777 791 zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub & 778 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * r dt_ice792 & + diag_vice + diag_vsnw - diag_adv_mass ) * e1e2t ) * rDt_ice 779 793 zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rDt_ice 780 794 781 795 ! -- salt diag -- ! 782 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * r dt_ice * 1.e-3796 zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rDt_ice * 1.e-3 783 797 zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rDt_ice * 1.e-3 784 798 -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn.F90
r13899 r14016 29 29 USE lbclnk ! lateral boundary conditions (or mpp links) 30 30 USE timing ! Timing 31 USE fldread ! read input fields 31 32 32 33 IMPLICIT NONE … … 51 52 REAL(wp) :: rn_vice ! prescribed v-vel (case np_dynADV1D & np_dynADV2D) 52 53 54 TYPE(FLD), ALLOCATABLE, DIMENSION(:) :: sf_icbmsk ! structure of input grounded icebergs mask (file informations, fields read) 55 53 56 !! * Substitutions 54 57 # include "do_loop_substitute.h90" … … 81 84 ! 82 85 ! controls 83 IF( ln_timing ) CALL timing_start('ice dyn')86 IF( ln_timing ) CALL timing_start('ice_dyn') 84 87 ! 85 88 IF( kt == nit000 .AND. lwp ) THEN … … 106 109 END WHERE 107 110 ! 111 IF( ln_landfast_L16 ) THEN 112 CALL fld_read( kt, 1, sf_icbmsk ) 113 icb_mask(:,:) = sf_icbmsk(1)%fnow(:,:,1) 114 ENDIF 108 115 ! 109 116 SELECT CASE( nice_dyn ) !-- Set which dynamics is running … … 172 179 ! 173 180 ! controls 174 IF( ln_timing ) CALL timing_stop ('ice dyn')181 IF( ln_timing ) CALL timing_stop ('ice_dyn') 175 182 ! 176 183 END SUBROUTINE ice_dyn … … 216 223 !! ** input : Namelist namdyn 217 224 !!------------------------------------------------------------------- 218 INTEGER :: ios, ioptio ! Local integer output status for namelist read 225 INTEGER :: ios, ioptio, ierror ! Local integer output status for namelist read 226 ! 227 CHARACTER(len=256) :: cn_dir ! Root directory for location of ice files 228 TYPE(FLD_N) :: sn_icbmsk ! informations about the grounded icebergs field to be read 219 229 !! 220 230 NAMELIST/namdyn/ ln_dynALL, ln_dynRHGADV, ln_dynADV1D, ln_dynADV2D, rn_uice, rn_vice, & 221 231 & rn_ishlat , & 222 & ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile 232 & ln_landfast_L16, rn_lf_depfra, rn_lf_bfr, rn_lf_relax, rn_lf_tensile, & 233 & sn_icbmsk, cn_dir 223 234 !!------------------------------------------------------------------- 224 235 ! … … 269 280 IF( .NOT.ln_landfast_L16 ) tau_icebfr(:,:) = 0._wp 270 281 ! 282 ! !--- allocate and fill structure for grounded icebergs mask 283 IF( ln_landfast_L16 ) THEN 284 ALLOCATE( sf_icbmsk(1), STAT=ierror ) 285 IF( ierror > 0 ) THEN 286 CALL ctl_stop( 'ice_dyn_init: unable to allocate sf_icbmsk structure' ) ; RETURN 287 ENDIF 288 ! 289 CALL fld_fill( sf_icbmsk, (/ sn_icbmsk /), cn_dir, 'ice_dyn_init', & 290 & 'landfast ice is a function of read grounded icebergs', 'icedyn' ) 291 ! 292 ALLOCATE( sf_icbmsk(1)%fnow(jpi,jpj,1) ) 293 IF( sf_icbmsk(1)%ln_tint ) ALLOCATE( sf_icbmsk(1)%fdta(jpi,jpj,1,2) ) 294 IF( TRIM(sf_icbmsk(1)%clrootname) == 'NOT USED' ) sf_icbmsk(1)%fnow(:,:,1) = 0._wp ! not used field (set to 0) 295 ELSE 296 icb_mask(:,:) = 0._wp 297 ENDIF 298 ! !--- other init 271 299 CALL ice_dyn_rdgrft_init ! set ice ridging/rafting parameters 272 300 CALL ice_dyn_rhg_init ! set ice rheology parameters -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_adv_pra.F90
r14012 r14016 156 156 157 157 ! diagnostics 158 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 158 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 159 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 159 160 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 160 161 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & … … 178 179 z0ei(:,:,jk,jl) = pe_i(:,:,jk,jl) * e1e2t(:,:) ! Ice heat content 179 180 END DO 180 IF ( ln_pnd_LEV ) THEN181 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 181 182 z0ap(:,:,jl) = pa_ip(:,:,jl) * e1e2t(:,:) ! Melt pond fraction 182 183 z0vp(:,:,jl) = pv_ip(:,:,jl) * e1e2t(:,:) ! Melt pond volume … … 214 215 END DO 215 216 ! 216 IF ( ln_pnd_LEV ) THEN217 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 217 218 CALL adv_x( zdt , zudy , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 218 219 CALL adv_y( zdt , zvdx , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) … … 249 250 & sxxe(:,:,jk,:), sye(:,:,jk,:), syye(:,:,jk,:), sxye(:,:,jk,:) ) 250 251 END DO 251 IF ( ln_pnd_LEV ) THEN252 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 252 253 CALL adv_y( zdt , zvdx , 1._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) !--- melt pond fraction 253 254 CALL adv_x( zdt , zudy , 0._wp , zarea , z0ap , sxap , sxxap , syap , syyap , sxyap ) … … 278 279 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ei , 'T', 1._wp, sxe , 'T', -1._wp, sye , 'T', -1._wp & ! ice enthalpy 279 280 & , sxxe , 'T', 1._wp, syye , 'T', 1._wp, sxye , 'T', 1._wp ) 280 IF ( ln_pnd_LEV ) THEN281 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 281 282 CALL lbc_lnk_multi( 'icedyn_adv_pra', z0ap , 'T', 1._wp, sxap , 'T', -1._wp, syap , 'T', -1._wp & ! melt pond fraction 282 283 & , sxxap, 'T', 1._wp, syyap, 'T', 1._wp, sxyap, 'T', 1._wp & … … 302 303 pe_i(:,:,jk,jl) = z0ei(:,:,jk,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 303 304 END DO 304 IF ( ln_pnd_LEV ) THEN305 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 305 306 pa_ip(:,:,jl) = z0ap(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) 306 307 pv_ip(:,:,jl) = z0vp(:,:,jl) * r1_e1e2t(:,:) * tmask(:,:,1) … … 320 321 ! 321 322 ! --- diagnostics --- ! 322 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 323 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 324 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 323 325 & - zdiag_adv_mass(:,:) ) * z1_dt 324 326 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & … … 769 771 ! ! -- check h_ip -- ! 770 772 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 771 IF( ln_pnd_LEV . AND. pv_ip(ji,jj,jl) > 0._wp ) THEN773 IF( ln_pnd_LEV .OR. ln_pnd_TOPO .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 772 774 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 773 775 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN … … 1015 1017 END DO 1016 1018 ! 1017 IF( ln_pnd_LEV ) THEN ! melt pond fraction1019 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN ! melt pond fraction 1018 1020 IF( iom_varid( numrir, 'sxap', ldstop = .FALSE. ) > 0 ) THEN 1019 1021 CALL iom_get( numrir, jpdom_auto, 'sxap' , sxap , psgn = -1._wp ) … … 1057 1059 sxc0 = 0._wp ; syc0 = 0._wp ; sxxc0 = 0._wp ; syyc0 = 0._wp ; sxyc0 = 0._wp ! snow layers heat content 1058 1060 sxe = 0._wp ; sye = 0._wp ; sxxe = 0._wp ; syye = 0._wp ; sxye = 0._wp ! ice layers heat content 1059 IF( ln_pnd_LEV ) THEN1061 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 1060 1062 sxap = 0._wp ; syap = 0._wp ; sxxap = 0._wp ; syyap = 0._wp ; sxyap = 0._wp ! melt pond fraction 1061 1063 sxvp = 0._wp ; syvp = 0._wp ; sxxvp = 0._wp ; syyvp = 0._wp ; sxyvp = 0._wp ! melt pond volume … … 1135 1137 END DO 1136 1138 ! 1137 IF( ln_pnd_LEV ) THEN ! melt pond fraction1139 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN ! melt pond fraction 1138 1140 CALL iom_rstput( iter, nitrst, numriw, 'sxap' , sxap ) 1139 1141 CALL iom_rstput( iter, nitrst, numriw, 'syap' , syap ) -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_adv_umx.F90
r13899 r14016 182 182 183 183 ! diagnostics 184 zdiag_adv_mass(:,:) = SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos 184 zdiag_adv_mass(:,:) = SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 185 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow 185 186 zdiag_adv_salt(:,:) = SUM( psv_i(:,:,:) , dim=3 ) * rhoi 186 187 zdiag_adv_heat(:,:) = - SUM(SUM( pe_i(:,:,1:nlay_i,:) , dim=4 ), dim=3 ) & … … 338 339 ! 339 340 !== melt ponds ==! 340 IF ( ln_pnd_LEV ) THEN341 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 341 342 ! concentration 342 343 zamsk = 1._wp … … 358 359 359 360 ! --- Lateral boundary conditions --- ! 360 IF ( ln_pnd_LEV.AND. ln_pnd_lids ) THEN361 IF ( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. ln_pnd_lids ) THEN 361 362 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 362 363 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp, pv_il,'T',1._wp ) 363 ELSEIF( ln_pnd_LEV.AND. .NOT.ln_pnd_lids ) THEN364 ELSEIF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. .NOT.ln_pnd_lids ) THEN 364 365 CALL lbc_lnk_multi( 'icedyn_adv_umx', pa_i,'T',1._wp, pv_i,'T',1._wp, pv_s,'T',1._wp, psv_i,'T',1._wp, poa_i,'T',1._wp & 365 366 & , pa_ip,'T',1._wp, pv_ip,'T',1._wp ) … … 379 380 ! 380 381 ! --- diagnostics --- ! 381 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i(:,:,:) , dim=3 ) * rhoi + SUM( pv_s(:,:,:) , dim=3 ) * rhos & 382 diag_adv_mass(:,:) = diag_adv_mass(:,:) + ( SUM( pv_i (:,:,:) , dim=3 ) * rhoi + SUM( pv_s (:,:,:) , dim=3 ) * rhos & 383 & + SUM( pv_ip(:,:,:) , dim=3 ) * rhow + SUM( pv_il(:,:,:) , dim=3 ) * rhow & 382 384 & - zdiag_adv_mass(:,:) ) * z1_dt 383 385 diag_adv_salt(:,:) = diag_adv_salt(:,:) + ( SUM( psv_i(:,:,:) , dim=3 ) * rhoi & … … 1497 1499 ! ! -- check h_ip -- ! 1498 1500 ! if h_ip is larger than the surrounding 9 pts => reduce h_ip and increase a_ip 1499 IF( ln_pnd_LEV.AND. pv_ip(ji,jj,jl) > 0._wp ) THEN1501 IF( ( ln_pnd_LEV .OR. ln_pnd_TOPO ) .AND. pv_ip(ji,jj,jl) > 0._wp ) THEN 1500 1502 zhip = pv_ip(ji,jj,jl) / MAX( epsi20, pa_ip(ji,jj,jl) ) 1501 1503 IF( zhip > phip_max(ji,jj,jl) .AND. pa_ip(ji,jj,jl) < 0.15 ) THEN -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_rdgrft.F90
r13899 r14016 140 140 INTEGER , DIMENSION(jpij) :: iptidx ! compute ridge/raft or not 141 141 REAL(wp), DIMENSION(jpij) :: zdivu, zdelt ! 1D divu_i & delta_i 142 REAL(wp), DIMENSION(jpij) :: zconv ! 1D rdg_conv (if EAP rheology) 142 143 ! 143 144 INTEGER, PARAMETER :: jp_itermax = 20 … … 175 176 ! just needed here 176 177 CALL tab_2d_1d( npti, nptidx(1:npti), zdelt (1:npti) , delta_i ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), zconv (1:npti) , rdg_conv ) 177 179 ! needed here and in the iteration loop 178 180 CALL tab_2d_1d( npti, nptidx(1:npti), zdivu (1:npti) , divu_i) ! zdivu is used as a work array here (no change in divu_i) … … 184 186 ! closing_net = rate at which open water area is removed + ice area removed by ridging 185 187 ! - ice area added in new ridges 186 closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 188 IF( ln_rhg_EVP .OR. ln_rhg_VP ) & 189 & closing_net(ji) = rn_csrdg * 0.5_wp * ( zdelt(ji) - ABS( zdivu(ji) ) ) - MIN( zdivu(ji), 0._wp ) 190 IF( ln_rhg_EAP ) closing_net(ji) = zconv(ji) 187 191 ! 188 192 IF( zdivu(ji) < 0._wp ) closing_net(ji) = MAX( closing_net(ji), -zdivu(ji) ) ! make sure the closing rate is large enough … … 575 579 oirft2(ji) = oa_i_2d(ji,jl1) * afrft * hi_hrft 576 580 577 IF ( ln_pnd_LEV ) THEN581 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 578 582 aprdg1 = a_ip_2d(ji,jl1) * afrdg 579 583 aprdg2(ji) = a_ip_2d(ji,jl1) * afrdg * hi_hrdg(ji,jl1) … … 612 616 sv_i_2d(ji,jl1) = sv_i_2d(ji,jl1) - sirdg1 - sirft(ji) 613 617 oa_i_2d(ji,jl1) = oa_i_2d(ji,jl1) - oirdg1 - oirft1 614 IF ( ln_pnd_LEV ) THEN618 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 615 619 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - aprdg1 - aprft1 616 620 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - vprdg(ji) - vprft(ji) … … 709 713 v_s_2d (ji,jl2) = v_s_2d (ji,jl2) + ( vsrdg (ji) * rn_fsnwrdg * fvol(ji) + & 710 714 & vsrft (ji) * rn_fsnwrft * zswitch(ji) ) 711 IF ( ln_pnd_LEV ) THEN715 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 712 716 v_ip_2d (ji,jl2) = v_ip_2d(ji,jl2) + ( vprdg (ji) * rn_fpndrdg * fvol (ji) & 713 717 & + vprft (ji) * rn_fpndrft * zswitch(ji) ) … … 776 780 ! !--------------------------------------------------! 777 781 strength(:,:) = rn_pstar * SUM( v_i(:,:,:), dim=3 ) * EXP( -rn_crhg * ( 1._wp - SUM( a_i(:,:,:), dim=3 ) ) ) 778 ismooth = 1 782 ismooth = 1 ! original code 783 ! ismooth = 0 ! try for EAP stability 779 784 ! !--------------------------------------------------! 780 785 ELSE ! Zero strength ! -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_rhg.F90
r13899 r14016 17 17 USE ice ! sea-ice: variables 18 18 USE icedyn_rhg_evp ! sea-ice: EVP rheology 19 USE icedyn_rhg_eap ! sea-ice: EAP rheology 20 USE icedyn_rhg_vp ! sea-ice: VP rheology 19 21 USE icectl ! sea-ice: control prints 20 22 ! … … 33 35 ! ! associated indices: 34 36 INTEGER, PARAMETER :: np_rhgEVP = 1 ! EVP rheology 35 !! INTEGER, PARAMETER :: np_rhgEAP = 2 ! EAP rheology 37 INTEGER, PARAMETER :: np_rhgEAP = 2 ! EAP rheology 38 INTEGER, PARAMETER :: np_rhgVP = 3 ! VP rheology 36 39 37 ! ** namelist (namrhg) **38 LOGICAL :: ln_rhg_EVP ! EVP rheology39 40 ! 40 41 !!---------------------------------------------------------------------- … … 77 78 ! !------------------------! 78 79 CALL ice_dyn_rhg_evp( kt, Kmm, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i ) 79 ! 80 ! 81 ! !------------------------! 82 CASE( np_rhgVP ) ! Viscous-Plastic ! 83 ! !------------------------! 84 CALL ice_dyn_rhg_vp ( kt, shear_i, divu_i, delta_i ) 85 ! 86 ! !----------------------------! 87 CASE( np_rhgEAP ) ! Elasto-Anisotropic-Plastic ! 88 ! !----------------------------! 89 CALL ice_dyn_rhg_eap( kt, Kmm, stress1_i, stress2_i, stress12_i, shear_i, divu_i, delta_i, aniso_11, aniso_12, rdg_conv ) 80 90 END SELECT 81 91 ! 82 IF( lrst_ice ) THEN !* write EVP fields in the restart file 83 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'WRITE', kt ) 92 IF( lrst_ice ) THEN 93 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'WRITE', kt ) !* write EVP fields in the restart file 94 IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'WRITE', kt ) !* write EAP fields in the restart file 95 ! MV note: no restart needed for VP as there is no time equation for stress tensor 84 96 ENDIF 85 97 ! … … 108 120 INTEGER :: ios, ioptio ! Local integer output status for namelist read 109 121 !! 110 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg 122 NAMELIST/namdyn_rhg/ ln_rhg_EVP, ln_aEVP, ln_rhg_EAP, rn_creepl, rn_ecc , nn_nevp, rn_relast, nn_rhg_chkcvg, & !-- evp 123 & ln_rhg_VP, nn_vp_nout, nn_vp_ninn, nn_vp_chkcvg !-- vp 111 124 !!------------------------------------------------------------------- 112 125 ! … … 124 137 WRITE(numout,*) ' rheology EVP (icedyn_rhg_evp) ln_rhg_EVP = ', ln_rhg_EVP 125 138 WRITE(numout,*) ' use adaptive EVP (aEVP) ln_aEVP = ', ln_aEVP 126 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl 127 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc 139 WRITE(numout,*) ' creep limit rn_creepl = ', rn_creepl ! also used by vp 140 WRITE(numout,*) ' eccentricity of the elliptical yield curve rn_ecc = ', rn_ecc ! also used by vp 128 141 WRITE(numout,*) ' number of iterations for subcycling nn_nevp = ', nn_nevp 129 142 WRITE(numout,*) ' ratio of elastic timescale over ice time step rn_relast = ', rn_relast 130 WRITE(numout,*) ' check convergence of rheology nn_rhg_chkcvg = ', nn_rhg_chkcvg 131 IF ( nn_rhg_chkcvg == 0 ) THEN ; WRITE(numout,*) ' no check' 132 ELSEIF( nn_rhg_chkcvg == 1 ) THEN ; WRITE(numout,*) ' check cvg at the main time step' 133 ELSEIF( nn_rhg_chkcvg == 2 ) THEN ; WRITE(numout,*) ' check cvg at both main and rheology time steps' 143 WRITE(numout,*) ' check convergence of rheology nn_rhg_chkcvg = ', nn_rhg_chkcvg 144 WRITE(numout,*) ' rheology VP (icedyn_rhg_VP) ln_rhg_VP = ', ln_rhg_VP 145 WRITE(numout,*) ' number of outer iterations nn_vp_nout = ', nn_vp_nout 146 WRITE(numout,*) ' number of inner iterations nn_vp_ninn = ', nn_vp_ninn 147 WRITE(numout,*) ' iteration step for convergence check nn_vp_chkcvg = ', nn_vp_chkcvg 148 IF( ln_rhg_EVP ) THEN 149 IF ( nn_rhg_chkcvg == 0 ) THEN ; WRITE(numout,*) ' no check cvg' 150 ELSEIF( nn_rhg_chkcvg == 1 ) THEN ; WRITE(numout,*) ' check cvg at the main time step' 151 ELSEIF( nn_rhg_chkcvg == 2 ) THEN ; WRITE(numout,*) ' check cvg at both main and rheology time steps' 152 ENDIF 134 153 ENDIF 154 WRITE(numout,*) ' rheology EAP (icedyn_rhg_eap) ln_rhg_EAP = ', ln_rhg_EAP 135 155 ENDIF 136 156 ! … … 138 158 ioptio = 0 139 159 IF( ln_rhg_EVP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEVP ; ENDIF 140 !! IF( ln_rhg_EAP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEAP ; ENDIF 160 IF( ln_rhg_EAP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgEAP ; ENDIF 161 IF( ln_rhg_VP ) THEN ; ioptio = ioptio + 1 ; nice_rhg = np_rhgVP ; ENDIF 141 162 IF( ioptio /= 1 ) CALL ctl_stop( 'ice_dyn_rhg_init: choose one and only one ice rheology' ) 142 163 ! 143 164 IF( ln_rhg_EVP ) CALL rhg_evp_rst( 'READ' ) !* read or initialize all required files 165 IF( ln_rhg_EAP ) CALL rhg_eap_rst( 'READ' ) !* read or initialize all required files 166 ! no restart for VP as there is no explicit time dependency in the equation 144 167 ! 145 168 END SUBROUTINE ice_dyn_rhg_init -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icedyn_rhg_evp.F90
r14012 r14016 326 326 zvV = 0.5_wp * ( vt_i(ji,jj) * e1e2t(ji,jj) + vt_i(ji,jj+1) * e1e2t(ji,jj+1) ) * r1_e1e2v(ji,jj) * vmask(ji,jj,1) 327 327 ! ice-bottom stress at U points 328 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) 328 zvCr = zaU(ji,jj) * rn_lf_depfra * hu(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 329 329 ztaux_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvU - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaU(ji,jj) ) ) 330 330 ! ice-bottom stress at V points 331 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) 331 zvCr = zaV(ji,jj) * rn_lf_depfra * hv(ji,jj,Kmm) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 332 332 ztauy_base(ji,jj) = - rn_lf_bfr * MAX( 0._wp, zvV - zvCr ) * EXP( -rn_crhg * ( 1._wp - zaV(ji,jj) ) ) 333 333 ! ice_bottom stress at T points 334 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) 334 zvCr = at_i(ji,jj) * rn_lf_depfra * ht(ji,jj) * ( 1._wp - icb_mask(ji,jj) ) ! if grounded icebergs are read: ocean depth = 0 335 335 tau_icebfr(ji,jj) = - rn_lf_bfr * MAX( 0._wp, vt_i(ji,jj) - zvCr ) * EXP( -rn_crhg * ( 1._wp - at_i(ji,jj) ) ) 336 336 END_2D -
NEMO/branches/2020/tickets_icb_1900/src/ICE/iceistate.F90
r13899 r14016 426 426 ! 4) Snow-ice mass (case ice is fully embedded) 427 427 !---------------------------------------------- 428 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s (:,:,:) + rhoi * v_i(:,:,:), dim=3 ) ! snow+ice mass428 snwice_mass (:,:) = tmask(:,:,1) * SUM( rhos * v_s + rhoi * v_i + rhow * ( v_ip + v_il ), dim=3 ) ! snow+ice mass 429 429 snwice_mass_b(:,:) = snwice_mass(:,:) 430 430 ! -
NEMO/branches/2020/tickets_icb_1900/src/ICE/iceitd.F90
r13899 r14016 29 29 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) 30 30 USE prtctl ! Print control 31 USE timing ! Timing 31 32 32 33 IMPLICIT NONE … … 87 88 REAL(wp), DIMENSION(jpij,0:jpl) :: zhbnew ! new boundaries of ice categories 88 89 !!------------------------------------------------------------------ 90 IF( ln_timing ) CALL timing_start('iceitd_rem') 89 91 90 92 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_rem: remapping ice thickness distribution' … … 315 317 IF ( a_i_1d(ji) > epsi10 .AND. h_i_1d(ji) < rn_himin ) THEN 316 318 a_i_1d(ji) = a_i_1d(ji) * h_i_1d(ji) / rn_himin 317 IF( ln_pnd_LEV ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin319 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) a_ip_1d(ji) = a_ip_1d(ji) * h_i_1d(ji) / rn_himin 318 320 h_i_1d(ji) = rn_himin 319 321 ENDIF … … 328 330 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_rem', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 329 331 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_rem', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 332 IF( ln_timing ) CALL timing_stop ('iceitd_rem') 330 333 ! 331 334 END SUBROUTINE ice_itd_rem … … 486 489 zaTsfn(ji,jl2) = zaTsfn(ji,jl2) + ztrans 487 490 ! 488 IF ( ln_pnd_LEV ) THEN491 IF ( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 489 492 ztrans = a_ip_2d(ji,jl1) * zworka(ji) ! Pond fraction 490 493 a_ip_2d(ji,jl1) = a_ip_2d(ji,jl1) - ztrans 491 494 a_ip_2d(ji,jl2) = a_ip_2d(ji,jl2) + ztrans 492 495 ! 493 ztrans = v_ip_2d(ji,jl1) * zwork a(ji) ! Pond volume (also proportional to da/a)496 ztrans = v_ip_2d(ji,jl1) * zworkv(ji) ! Pond volume 494 497 v_ip_2d(ji,jl1) = v_ip_2d(ji,jl1) - ztrans 495 498 v_ip_2d(ji,jl2) = v_ip_2d(ji,jl2) + ztrans 496 499 ! 497 500 IF ( ln_pnd_lids ) THEN ! Pond lid volume 498 ztrans = v_il_2d(ji,jl1) * zwork a(ji)501 ztrans = v_il_2d(ji,jl1) * zworkv(ji) 499 502 v_il_2d(ji,jl1) = v_il_2d(ji,jl1) - ztrans 500 503 v_il_2d(ji,jl2) = v_il_2d(ji,jl2) + ztrans … … 606 609 REAL(wp), DIMENSION(jpij,jpl-1) :: zdaice, zdvice ! ice area and volume transferred 607 610 !!------------------------------------------------------------------ 611 IF( ln_timing ) CALL timing_start('iceitd_reb') 608 612 ! 609 613 IF( kt == nit000 .AND. lwp ) WRITE(numout,*) '-- ice_itd_reb: rebining ice thickness distribution' … … 635 639 jdonor(ji,jl) = jl 636 640 ! how much of a_i you send in cat sup is somewhat arbitrary 637 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 638 !! zdaice(ji,jl) = a_i_1d(ji) * ( h_i_1d(ji) - hi_max(jl) + epsi10 ) / h_i_1d(ji) 639 !! zdvice(ji,jl) = v_i_1d(ji) - ( a_i_1d(ji) - zdaice(ji,jl) ) * ( hi_max(jl) - epsi10 ) 640 !!clem: these do not work properly after a restart (I do not know why) => not sure it is still true 641 !! zdaice(ji,jl) = a_i_1d(ji) 642 !! zdvice(ji,jl) = v_i_1d(ji) 643 !!clem: these are from UCL and work ok 644 zdaice(ji,jl) = a_i_1d(ji) * 0.5_wp 645 zdvice(ji,jl) = v_i_1d(ji) - zdaice(ji,jl) * ( hi_max(jl) + hi_max(jl-1) ) * 0.5_wp 641 ! these are from CICE => transfer everything 642 !!zdaice(ji,jl) = a_i_1d(ji) 643 !!zdvice(ji,jl) = v_i_1d(ji) 644 ! these are from LLN => transfer only half of the category 645 zdaice(ji,jl) = 0.5_wp * a_i_1d(ji) 646 zdvice(ji,jl) = v_i_1d(ji) - (1._wp - 0.5_wp) * a_i_1d(ji) * hi_mean(jl) 646 647 END DO 647 648 ! … … 686 687 IF( ln_icediachk ) CALL ice_cons_hsm(1, 'iceitd_reb', rdiag_v, rdiag_s, rdiag_t, rdiag_fv, rdiag_fs, rdiag_ft) 687 688 IF( ln_icediachk ) CALL ice_cons2D (1, 'iceitd_reb', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 689 IF( ln_timing ) CALL timing_stop ('iceitd_reb') 688 690 ! 689 691 END SUBROUTINE ice_itd_reb -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icesbc.F90
r13899 r14016 62 62 !!------------------------------------------------------------------- 63 63 ! 64 IF( ln_timing ) CALL timing_start('ice _sbc')64 IF( ln_timing ) CALL timing_start('icesbc') 65 65 ! 66 66 IF( kt == nit000 .AND. lwp ) THEN … … 89 89 ENDIF 90 90 ! 91 IF( ln_timing ) CALL timing_stop('ice _sbc')91 IF( ln_timing ) CALL timing_stop('icesbc') 92 92 ! 93 93 END SUBROUTINE ice_sbc_tau … … 122 122 !!-------------------------------------------------------------------- 123 123 ! 124 IF( ln_timing ) CALL timing_start('ice _sbc_flx')124 IF( ln_timing ) CALL timing_start('icesbc') 125 125 126 126 IF( kt == nit000 .AND. lwp ) THEN … … 176 176 ENDIF 177 177 ! 178 IF( ln_timing ) CALL timing_stop('ice _sbc_flx')178 IF( ln_timing ) CALL timing_stop('icesbc') 179 179 ! 180 180 END SUBROUTINE ice_sbc_flx -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icestp.F90
r14012 r14016 121 121 !!---------------------------------------------------------------------- 122 122 ! 123 IF( ln_timing ) CALL timing_start('ice _stp')123 IF( ln_timing ) CALL timing_start('icestp') 124 124 ! 125 125 ! !-----------------------! … … 215 215 !!gm remark, the ocean-ice stress is not saved in ice diag call above ..... find a solution!!! 216 216 ! 217 IF( ln_timing ) CALL timing_stop('ice _stp')217 IF( ln_timing ) CALL timing_stop('icestp') 218 218 ! 219 219 END SUBROUTINE ice_stp … … 373 373 v_i_b (:,:,:) = v_i (:,:,:) ! ice volume 374 374 v_s_b (:,:,:) = v_s (:,:,:) ! snow volume 375 v_ip_b(:,:,:) = v_ip(:,:,:) ! pond volume 376 v_il_b(:,:,:) = v_il(:,:,:) ! pond lid volume 375 377 sv_i_b(:,:,:) = sv_i(:,:,:) ! salt content 376 378 e_s_b (:,:,:,:) = e_s (:,:,:,:) ! snow thermal energy … … 432 434 diag_heat(ji,jj) = 0._wp ; diag_sice(ji,jj) = 0._wp 433 435 diag_vice(ji,jj) = 0._wp ; diag_vsnw(ji,jj) = 0._wp 436 diag_aice(ji,jj) = 0._wp ; diag_vpnd(ji,jj) = 0._wp 434 437 435 438 tau_icebfr (ji,jj) = 0._wp ! landfast ice param only (clem: important to keep the init here) … … 457 460 qcn_ice (ji,jj,jl) = 0._wp ! conductive flux (ln_cndflx=T & ln_cndemule=T) 458 461 qtr_ice_bot(ji,jj,jl) = 0._wp ! part of solar radiation transmitted through the ice needed at least for outputs 462 ! Melt pond surface melt diagnostics (mv - more efficient: grouped into one water volume flux) 463 dh_i_sum_2d(ji,jj,jl) = 0._wp 464 dh_s_mlt_2d(ji,jj,jl) = 0._wp 459 465 END_2D 460 466 ENDDO … … 485 491 diag_vsnw(:,:) = diag_vsnw(:,:) & 486 492 & + SUM( v_s (:,:,:) - v_s_b (:,:,:) , dim=3 ) * r1_Dt_ice * rhos 493 diag_vpnd(:,:) = diag_vpnd(:,:) & 494 & + SUM( v_ip + v_il - v_ip_b - v_il_b , dim=3 ) * r1_Dt_ice * rhow 487 495 ! 488 496 IF( kn == 2 ) CALL iom_put ( 'hfxdhc' , diag_heat ) ! output of heat trend -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icethd.F90
r13899 r14016 166 166 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 167 167 168 ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 169 ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 170 IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 171 zqfr = 0._wp 172 zqfr_pos = 0._wp 173 qsb_ice_bot(ji,jj) = 0._wp 174 ENDIF 175 ! 168 176 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 169 177 ! qlead is the energy received from the atm. in the leads. … … 239 247 IF( ln_icedH ) THEN ! --- Growing/Melting --- ! 240 248 CALL ice_thd_dh ! Ice-Snow thickness 241 CALL ice_thd_pnd ! Melt ponds formation242 249 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 243 250 ENDIF … … 260 267 IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 261 268 ! 269 IF ( ln_pnd .AND. ln_icedH ) & 270 & CALL ice_thd_pnd ! --- Melt ponds 271 ! 262 272 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! 263 273 ! … … 266 276 CALL ice_cor( kt , 2 ) ! --- Corrections --- ! 267 277 ! 268 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * r dt_ice ! ice natural aging incrementation278 oa_i(:,:,:) = oa_i(:,:,:) + a_i(:,:,:) * rDt_ice ! ice natural aging incrementation 269 279 ! 270 280 ! convergence tests … … 377 387 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 378 388 END DO 379 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )380 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )381 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) )382 389 ! 383 390 CALL tab_2d_1d( npti, nptidx(1:npti), qprec_ice_1d (1:npti), qprec_ice ) … … 409 416 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 410 417 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 411 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )412 418 ! 413 419 CALL tab_2d_1d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) … … 464 470 v_s_1d (1:npti) = h_s_1d (1:npti) * a_i_1d (1:npti) 465 471 sv_i_1d(1:npti) = s_i_1d (1:npti) * v_i_1d (1:npti) 466 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)467 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)468 472 oa_i_1d(1:npti) = o_i_1d (1:npti) * a_i_1d (1:npti) 469 473 … … 483 487 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,kl) ) 484 488 END DO 485 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,kl) )486 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,kl) )487 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,kl) )488 489 ! 489 490 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sni_1d(1:npti), wfx_snw_sni ) … … 501 502 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_spr_1d (1:npti), wfx_spr ) 502 503 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_lam_1d (1:npti), wfx_lam ) 503 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )504 504 ! 505 505 CALL tab_1d_2d( npti, nptidx(1:npti), sfx_bog_1d (1:npti), sfx_bog ) … … 529 529 CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) ) 530 530 CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) ) 531 ! Melt ponds 532 CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum (1:npti) , dh_i_sum_2d(:,:,kl) ) 533 CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt (1:npti) , dh_s_mlt_2d(:,:,kl) ) 531 534 ! SIMIP diagnostics 532 535 CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d (1:npti), t_si (:,:,kl) ) … … 537 540 CALL tab_1d_2d( npti, nptidx(1:npti), v_s_1d (1:npti), v_s (:,:,kl) ) 538 541 CALL tab_1d_2d( npti, nptidx(1:npti), sv_i_1d(1:npti), sv_i(:,:,kl) ) 539 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,kl) )540 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,kl) )541 542 CALL tab_1d_2d( npti, nptidx(1:npti), oa_i_1d(1:npti), oa_i(:,:,kl) ) 542 543 ! check convergence of heat diffusion scheme -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icethd_dh.F90
r13899 r14016 55 55 !! - Snow ice formation 56 56 !! 57 !! ** Note : h=max(0,h+dh) are often used to ensure positivity of h. 58 !! very small negative values can occur otherwise (e.g. -1.e-20) 59 !! 57 60 !! References : Bitz and Lipscomb, 1999, J. Geophys. Res. 58 61 !! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646 … … 79 82 REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean 80 83 81 REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)82 84 REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2) 83 85 REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2) … … 85 87 REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2) 86 88 REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2) 87 88 REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt 89 REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation 90 REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation 91 92 REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness 93 REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness 94 REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah 95 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting 96 89 REAL(wp), DIMENSION(jpij) :: zdeltah 97 90 REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing 91 92 INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanishing by melting 93 REAL(wp), DIMENSION(jpij,0:nlay_i+1) :: zh_i ! ice layer thickness (m) 94 REAL(wp), DIMENSION(jpij,0:nlay_s ) :: zh_s ! snw layer thickness (m) 95 REAL(wp), DIMENSION(jpij,0:nlay_s ) :: ze_s ! snw layer enthalpy (J.m-3) 98 96 99 97 REAL(wp) :: zswitch_sal … … 108 106 END SELECT 109 107 110 ! initialize layer thicknesses and enthalpies 108 ! initialize ice layer thicknesses and enthalpies 109 eh_i_old(1:npti,0:nlay_i+1) = 0._wp 111 110 h_i_old (1:npti,0:nlay_i+1) = 0._wp 112 eh_i_old(1:npti,0:nlay_i+1) = 0._wp111 zh_i (1:npti,0:nlay_i+1) = 0._wp 113 112 DO jk = 1, nlay_i 114 113 DO ji = 1, npti 114 eh_i_old(ji,jk) = h_i_1d(ji) * r1_nlay_i * e_i_1d(ji,jk) 115 115 h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i 116 eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk) 116 zh_i (ji,jk) = h_i_1d(ji) * r1_nlay_i 117 END DO 118 END DO 119 ! 120 ! initialize snw layer thicknesses and enthalpies 121 zh_s(1:npti,0) = 0._wp 122 ze_s(1:npti,0) = 0._wp 123 DO jk = 1, nlay_s 124 DO ji = 1, npti 125 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s 126 ze_s(ji,jk) = e_s_1d(ji,jk) 117 127 END DO 118 128 END DO … … 141 151 zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji) + qtr_ice_bot_1d(ji) * frq_m_1d(ji) 142 152 zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice ) 143 END DO144 145 ! Ice and snow layer thicknesses146 !-------------------------------147 DO jk = 1, nlay_i148 DO ji = 1, npti149 zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i150 END DO151 END DO152 153 DO jk = 1, nlay_s154 DO ji = 1, npti155 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s156 END DO157 153 END DO 158 154 … … 167 163 DO ji = 1, npti 168 164 IF( t_s_1d(ji,jk) > rt0 ) THEN 169 hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0170 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos 165 hfx_res_1d (ji) = hfx_res_1d (ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 166 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux 171 167 ! updates 172 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)173 h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk)168 dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk) 169 h_s_1d (ji) = MAX( 0._wp, h_s_1d (ji) - zh_s(ji,jk) ) 174 170 zh_s (ji,jk) = 0._wp 175 e_s_1d (ji,jk) = 0._wp 176 t_s_1d (ji,jk) = rt0 171 ze_s (ji,jk) = 0._wp 177 172 END IF 178 173 END DO … … 181 176 ! Snow precipitation 182 177 !------------------- 183 CALL ice_var_snwblow( 1.0_wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 184 185 zdeltah(1:npti,:) = 0._wp 178 CALL ice_var_snwblow( 1._wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing 179 186 180 DO ji = 1, npti 187 181 IF( sprecip_1d(ji) > 0._wp ) THEN 182 zh_s(ji,0) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness of precip 183 ze_s(ji,0) = MAX( 0._wp, - qprec_ice_1d(ji) ) ! enthalpy of the precip (>0, J.m-3) 188 184 ! 189 ! --- precipitation --- 190 zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness change 191 zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3) 185 hfx_spr_1d(ji) = hfx_spr_1d(ji) + ze_s(ji,0) * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2) 186 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * zh_s(ji,0) * a_i_1d(ji) * r1_Dt_ice ! mass flux, <0 192 187 ! 193 hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2)194 wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice ! mass flux, <0195 196 ! --- melt of falling snow ---197 rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )198 zdeltah (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change199 zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting200 hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat used to melt snow (W.m-2, >0)201 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip), >0202 203 ! updates available heat + precipitations after melting204 dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)205 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )206 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)207 208 188 ! update thickness 209 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) ) 210 ! 211 ELSE 212 ! 213 zdh_s_pre(ji) = 0._wp 214 zqprec (ji) = 0._wp 215 ! 189 h_s_1d(ji) = h_s_1d(ji) + zh_s(ji,0) 216 190 ENDIF 217 END DO218 219 ! recalculate snow layers220 DO jk = 1, nlay_s221 DO ji = 1, npti222 zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s223 END DO224 191 END DO 225 192 226 193 ! Snow melting 227 194 ! ------------ 228 ! If heat still available (zq_top > 0), then melt more snow 229 zdeltah(1:npti,:) = 0._wp 230 zdh_s_mel(1:npti) = 0._wp 231 DO jk = 1, nlay_s 195 ! If heat still available (zq_top > 0) 196 ! then all snw precip has been melted and we need to melt more snow 197 DO jk = 0, nlay_s 232 198 DO ji = 1, npti 233 199 IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN 234 200 ! 235 rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) ) 236 zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change 237 zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting 238 zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk) 239 240 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0) 241 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip) 201 rswitch = MAX( 0._wp , SIGN( 1._wp , ze_s(ji,jk) - epsi20 ) ) 202 zdum = - rswitch * zq_top(ji) / MAX( ze_s(ji,jk), epsi20 ) ! thickness change 203 zdum = MAX( zdum , - zh_s(ji,jk) ) ! bound melting 204 205 hfx_snw_1d (ji) = hfx_snw_1d (ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0) 206 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! snow melting only = water into the ocean 242 207 243 208 ! updates available heat + thickness 244 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk) 245 zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) ) 246 h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) ) 247 zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) ) 209 dh_s_mlt(ji) = dh_s_mlt(ji) + zdum 210 zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdum * ze_s(ji,jk) ) 211 h_s_1d (ji) = MAX( 0._wp , h_s_1d (ji) + zdum ) 212 zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) 213 !!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp 248 214 ! 249 215 ENDIF … … 255 221 ! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates 256 222 ! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean) 257 zdeltah(1:npti,:) = 0._wp 223 zdeltah (1:npti) = 0._wp ! total snow thickness that sublimates, < 0 224 zevap_rema(1:npti) = 0._wp 258 225 DO ji = 1, npti 259 226 IF( evap_ice_1d(ji) > 0._wp ) THEN 227 zdeltah (ji) = MAX( - evap_ice_1d(ji) * r1_rhos * rDt_ice, - h_s_1d(ji) ) ! amount of snw that sublimates, < 0 228 zevap_rema(ji) = MAX( 0._wp, evap_ice_1d(ji) * rDt_ice + zdeltah(ji) * rhos ) ! remaining evap in kg.m-2 (used for ice sublimation later on) 229 ENDIF 230 END DO 231 232 DO jk = 0, nlay_s 233 DO ji = 1, npti 234 zdum = MAX( -zh_s(ji,jk), zdeltah(ji) ) ! snow layer thickness that sublimates, < 0 260 235 ! 261 zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rDt_ice ) 262 zevap_rema(ji) = evap_ice_1d(ji) * rDt_ice + zdh_s_sub(ji) * rhos ! remaining evap in kg.m-2 (used for ice melting later on) 263 zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) ) 264 265 hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow) 266 & ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) & 267 & * a_i_1d(ji) * r1_Dt_ice 268 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice ! Mass flux by sublimation 269 270 ! new snow thickness 271 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) ) 272 ! update precipitations after sublimation and correct sublimation 273 zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1) 274 zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1) 275 ! 276 ELSE 277 ! 278 zdh_s_sub (ji) = 0._wp 279 zevap_rema(ji) = 0._wp 280 ! 281 ENDIF 282 END DO 283 284 ! --- Update snow diags --- ! 285 DO ji = 1, npti 286 dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji) 287 END DO 288 289 ! Update temperature, energy 290 !--------------------------- 291 ! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow) 292 DO jk = 1, nlay_s 293 DO ji = 1,npti 294 rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) ) 295 e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * & 296 & ( ( zdh_s_pre(ji) ) * zqprec(ji) + & 297 & ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) ) 298 END DO 299 END DO 300 236 hfx_sub_1d (ji) = hfx_sub_1d (ji) + ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux of snw that sublimates [W.m-2], < 0 237 wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux by sublimation 238 239 ! update thickness 240 h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdum ) 241 zh_s (ji,jk) = MAX( 0._wp , zh_s (ji,jk) + zdum ) 242 !!$ IF( zh_s(ji,jk) == 0._wp ) ze_s(ji,jk) = 0._wp 243 244 ! update sublimation left 245 zdeltah(ji) = MIN( zdeltah(ji) - zdum, 0._wp ) 246 END DO 247 END DO 248 249 ! 301 250 ! ! ============ ! 302 251 ! ! Ice ! … … 305 254 ! Surface ice melting 306 255 !-------------------- 307 zdeltah(1:npti,:) = 0._wp ! important308 256 DO jk = 1, nlay_i 309 257 DO ji = 1, npti … … 313 261 314 262 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0] 315 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 316 ! set up at 0 since no energy is needed to melt water...(it is already melted) 317 zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 318 ! this should normally not happen, but sometimes, heat diffusion leads to this 319 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 320 321 dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting 322 323 zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0] 324 325 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 326 ! ice enthalpy zEi is "sent" to the ocean 327 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 328 ! using s_i_1d and not sz_i_1d(jk) is ok 329 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 330 263 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 264 ! set up at 0 since no energy is needed to melt water...(it is already melted) 265 zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 266 ! this should normally not happen, but sometimes, heat diffusion leads to this 267 zfmdt = - zdum * rhoi ! Recompute mass flux [kg/m2, >0] 268 ! 269 dh_i_itm(ji) = dh_i_itm(ji) + zdum ! Cumulate internal melting 270 ! 271 hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 272 ! ice enthalpy zEi is "sent" to the ocean 273 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 274 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux 275 ! using s_i_1d and not sz_i_1d(jk) is ok 331 276 ELSE !-- Surface melting 332 277 … … 337 282 zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0] 338 283 339 zd eltah(ji,jk)= - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]340 341 zd eltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]342 343 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk)* rhoi * zdE ) ! update available heat344 345 dh_i_sum(ji) = dh_i_sum(ji) + zd eltah(ji,jk)! Cumulate surface melt346 347 zfmdt = - rhoi * zd eltah(ji,jk)! Recompute mass flux [kg/m2, >0]284 zdum = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0] 285 286 zdum = MIN( 0._wp , MAX( zdum , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0] 287 288 zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdum * rhoi * zdE ) ! update available heat 289 290 dh_i_sum(ji) = dh_i_sum(ji) + zdum ! Cumulate surface melt 291 292 zfmdt = - rhoi * zdum ! Recompute mass flux [kg/m2, >0] 348 293 349 294 zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0] 350 295 351 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 352 ! using s_i_1d and not sz_i_1d(jk) is ok) 353 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux [W.m-2], < 0 354 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0 355 ! 356 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 357 296 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 297 hfx_sum_1d(ji) = hfx_sum_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0 298 wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 299 sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0 300 ! using s_i_1d and not sz_i_1d(jk) is ok) 358 301 END IF 359 302 ! update thickness 303 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 304 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 305 ! 306 ! update heat content (J.m-2) and layer thickness 307 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 308 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 309 ! 310 ! 360 311 ! Ice sublimation 361 312 ! --------------- 362 zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi ) 363 zdeltah (ji,jk) = zdeltah (ji,jk) + zdum 364 dh_i_sub(ji) = dh_i_sub(ji) + zdum 365 366 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0 367 ! clem: flux is sent to the ocean for simplicity 368 ! but salt should remain in the ice except 369 ! if all ice is melted. => must be corrected 370 hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 371 372 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice ! Mass flux > 0 373 374 ! update remaining mass flux 375 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 376 313 zdum = MAX( - zh_i(ji,jk) , - zevap_rema(ji) * r1_rhoi ) 314 ! 315 hfx_sub_1d(ji) = hfx_sub_1d(ji) + e_i_1d(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0 316 wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux > 0 317 sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux >0 318 ! clem: flux is sent to the ocean for simplicity 319 ! but salt should remain in the ice except 320 ! if all ice is melted. => must be corrected 321 ! update remaining mass flux and thickness 322 zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi 323 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 324 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 325 dh_i_sub(ji) = dh_i_sub(ji) + zdum 326 327 ! update heat content (J.m-2) and layer thickness 328 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 329 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 330 377 331 ! record which layers have disappeared (for bottom melting) 378 332 ! => icount=0 : no layer has vanished 379 333 ! => icount=5 : 5 layers have vanished 380 rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk)) ) )334 rswitch = MAX( 0._wp , SIGN( 1._wp , - zh_i(ji,jk) ) ) 381 335 icount(ji,jk) = NINT( rswitch ) 382 zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )383 336 384 ! update heat content (J.m-2) and layer thickness 385 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 386 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 387 END DO 388 END DO 389 390 ! update ice thickness 391 DO ji = 1, npti 392 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) ) 393 END DO 394 337 END DO 338 END DO 339 395 340 ! remaining "potential" evap is sent to ocean 396 341 DO ji = 1, npti … … 430 375 & + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 ) 431 376 432 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity433 434 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)435 436 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)437 438 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)439 & - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp* ztmelts440 441 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)442 443 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)444 445 dh_i_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )377 s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity 378 379 ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C) 380 381 zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i) 382 383 zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0) 384 & - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts 385 386 zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0) 387 388 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 389 390 dh_i_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) ) 446 391 447 392 END DO 448 393 ! Contribution to Energy and Salt Fluxes 449 zfmdt = - rhoi * dh_i_bog(ji)! Mass flux x time step (kg/m2, < 0)394 zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0) 450 395 451 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0 452 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0 453 454 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice ! Salt flux, <0 455 456 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice ! Mass flux, <0 396 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0 397 hfx_bog_1d(ji) = hfx_bog_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0 398 wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * a_i_1d(ji) * r1_Dt_ice ! Mass flux, <0 399 sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * dh_i_bog(ji) * s_i_new(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux, <0 400 401 ! update thickness 402 zh_i(ji,nlay_i+1) = zh_i(ji,nlay_i+1) + dh_i_bog(ji) 403 h_i_1d(ji) = h_i_1d(ji) + dh_i_bog(ji) 457 404 458 405 ! update heat content (J.m-2) and layer thickness … … 466 413 ! Ice Basal melt 467 414 !--------------- 468 zdeltah(1:npti,:) = 0._wp ! important469 415 DO jk = nlay_i, 1, -1 470 416 DO ji = 1, npti … … 475 421 IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting 476 422 477 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 478 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 479 ! set up at 0 since no energy is needed to melt water...(it is already melted) 480 zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 481 ! this should normally not happen, but sometimes, heat diffusion leads to this 482 483 dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk) 484 485 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 486 487 hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 488 ! ice enthalpy zEi is "sent" to the ocean 489 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 490 ! using s_i_1d and not sz_i_1d(jk) is ok 491 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 492 493 ! update heat content (J.m-2) and layer thickness 494 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 495 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 496 423 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 424 zdE = 0._wp ! Specific enthalpy difference (J/kg, <0) 425 ! set up at 0 since no energy is needed to melt water...(it is already melted) 426 zdum = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing 427 ! this should normally not happen, but sometimes, heat diffusion leads to this 428 dh_i_itm (ji) = dh_i_itm(ji) + zdum 429 ! 430 zfmdt = - zdum * rhoi ! Mass flux x time step > 0 431 ! 432 hfx_res_1d(ji) = hfx_res_1d(ji) + zEi * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 433 ! ice enthalpy zEi is "sent" to the ocean 434 wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 435 sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux 436 ! using s_i_1d and not sz_i_1d(jk) is ok 497 437 ELSE !-- Basal melting 498 438 499 zEi = - e_i_1d(ji,jk) * r1_rhoi! Specific enthalpy of melting ice (J/kg, <0)500 zEw = rcp * ztmelts! Specific enthalpy of meltwater (J/kg, <0)501 zdE = zEi - zEw! Specific enthalpy difference (J/kg, <0)502 503 zfmdt = - zq_bot(ji) / zdE! Mass flux x time step (kg/m2, >0)504 505 zd eltah(ji,jk) = - zfmdt * r1_rhoi! Gross thickness change506 507 zd eltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change439 zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0) 440 zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0) 441 zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0) 442 443 zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0) 444 445 zdum = - zfmdt * r1_rhoi ! Gross thickness change 446 447 zdum = MIN( 0._wp , MAX( zdum, - zh_i(ji,jk) ) ) ! bound thickness change 508 448 509 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 510 511 dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt 512 513 zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0 514 515 zQm = zfmdt * zEw ! Heat exchanged with ocean 516 517 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 518 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat used in this process [W.m-2], >0 519 520 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux 521 ! using s_i_1d and not sz_i_1d(jk) is ok 522 523 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux 524 525 ! update heat content (J.m-2) and layer thickness 526 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk) 527 h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk) 449 zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdum * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors 450 451 dh_i_bom(ji) = dh_i_bom(ji) + zdum ! Update basal melt 452 453 zfmdt = - zdum * rhoi ! Mass flux x time step > 0 454 455 zQm = zfmdt * zEw ! Heat exchanged with ocean 456 457 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0 458 hfx_bom_1d(ji) = hfx_bom_1d(ji) - zdE * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat used in this process [W.m-2], >0 459 wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * zdum * a_i_1d(ji) * r1_Dt_ice ! Mass flux 460 sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * zdum * s_i_1d(ji) * a_i_1d(ji) * r1_Dt_ice ! Salt flux 461 ! using s_i_1d and not sz_i_1d(jk) is ok 528 462 ENDIF 529 463 ! update thickness 464 zh_i(ji,jk) = MAX( 0._wp, zh_i(ji,jk) + zdum ) 465 h_i_1d(ji) = MAX( 0._wp, h_i_1d(ji) + zdum ) 466 ! 467 ! update heat content (J.m-2) and layer thickness 468 eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdum * e_i_1d(ji,jk) 469 h_i_old (ji,jk) = h_i_old (ji,jk) + zdum 530 470 ENDIF 531 471 END DO 532 472 END DO 533 473 534 ! Update temperature, energy 535 ! -------------------------- 536 DO ji = 1, npti 537 h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) ) 538 END DO 539 540 ! If heat still available then melt more snow 541 !------------------------------------------- 542 zdeltah(1:npti,:) = 0._wp ! important 543 DO ji = 1, npti 544 zq_rema (ji) = zq_top(ji) + zq_bot(ji) 545 rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow 546 rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) ) 547 zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 ) 548 zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting 549 dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1) 550 h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1) 551 552 zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2) 553 hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice ! Heat used to melt snow, W.m-2 (>0) 554 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! Mass flux 555 dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1) 556 ! 557 ! Remaining heat flux (W.m-2) is sent to the ocean heat budget 558 !!hfx_res_1d(ji) = hfx_res_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice 559 560 IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji) 561 END DO 562 563 ! 474 ! Remove snow if ice has melted entirely 475 ! -------------------------------------- 476 DO jk = 0, nlay_s 477 DO ji = 1,npti 478 IF( h_i_1d(ji) == 0._wp ) THEN 479 ! mass & energy loss to the ocean 480 hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 481 wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux 482 483 ! update thickness and energy 484 h_s_1d(ji) = 0._wp 485 ze_s (ji,jk) = 0._wp 486 zh_s (ji,jk) = 0._wp 487 ENDIF 488 END DO 489 END DO 490 491 ! Snow load on ice 492 ! ----------------- 493 ! When snow load exceeds Archimede's limit and sst is positive, 494 ! snow-ice formation (next bloc) can lead to negative ice enthalpy. 495 ! Therefore we consider here that this excess of snow falls into the ocean 496 zdeltah(1:npti) = h_s_1d(1:npti) + h_i_1d(1:npti) * (rhoi-rho0) * r1_rhos 497 DO jk = 0, nlay_s 498 DO ji = 1, npti 499 IF( zdeltah(ji) > 0._wp .AND. sst_1d(ji) > 0._wp ) THEN 500 ! snow layer thickness that falls into the ocean 501 zdum = MIN( zdeltah(ji) , zh_s(ji,jk) ) 502 ! mass & energy loss to the ocean 503 hfx_res_1d(ji) = hfx_res_1d(ji) - ze_s(ji,jk) * zdum * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0 504 wfx_res_1d(ji) = wfx_res_1d(ji) + rhos * zdum * a_i_1d(ji) * r1_Dt_ice ! mass flux 505 ! update thickness and energy 506 h_s_1d(ji) = MAX( 0._wp, h_s_1d(ji) - zdum ) 507 zh_s (ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) 508 ! update snow thickness that still has to fall 509 zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum ) 510 ENDIF 511 END DO 512 END DO 513 564 514 ! Snow-Ice formation 565 515 ! ------------------ 566 ! When snow load exce sses Archimede's limit, snow-ice interface goes down under sea-level,567 ! flooding of seawater transforms snow into ice dh_snowice is positive for the ice516 ! When snow load exceeds Archimede's limit, snow-ice interface goes down under sea-level, 517 ! flooding of seawater transforms snow into ice. Thickness that is transformed is dh_snowice (positive for the ice) 568 518 z1_rho = 1._wp / ( rhos+rho0-rhoi ) 519 zdeltah(1:npti) = 0._wp 569 520 DO ji = 1, npti 570 521 ! 571 dh_snowice(ji) = MAX( 522 dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho ) 572 523 573 524 h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji) … … 579 530 zQm = zfmdt * zEw 580 531 581 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux 582 583 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux 532 hfx_thd_1d(ji) = hfx_thd_1d(ji) + zEw * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Heat flux 533 sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice ! Salt flux 584 534 585 535 ! Case constant salinity in time: virtual salt flux to keep salinity constant 586 536 IF( nn_icesal /= 2 ) THEN 587 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice &! put back sss_m into the ocean588 & - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi* r1_Dt_ice ! and get rn_icesal from the ocean537 sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d(ji) * zfmdt * a_i_1d(ji) * r1_Dt_ice & ! put back sss_m into the ocean 538 & - s_i_1d(ji) * dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice ! and get rn_icesal from the ocean 589 539 ENDIF 590 540 591 541 ! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume 592 wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice 593 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice 542 wfx_sni_1d (ji) = wfx_sni_1d (ji) - dh_snowice(ji) * rhoi * a_i_1d(ji) * r1_Dt_ice 543 wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + dh_snowice(ji) * rhos * a_i_1d(ji) * r1_Dt_ice 544 545 ! update thickness 546 zh_i(ji,0) = zh_i(ji,0) + dh_snowice(ji) 547 zdeltah(ji) = dh_snowice(ji) 594 548 595 549 ! update heat content (J.m-2) and layer thickness 596 eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw597 550 h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji) 598 599 END DO 600 601 ! 602 ! Update temperature, energy 603 ! -------------------------- 604 DO ji = 1, npti 605 rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) ) 606 t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0 607 END DO 608 551 eh_i_old(ji,0) = eh_i_old(ji,0) + zfmdt * zEw ! 1st part (sea water enthalpy) 552 553 END DO 554 ! 555 DO jk = nlay_s, 0, -1 ! flooding of snow starts from the base 556 DO ji = 1, npti 557 zdum = MIN( zdeltah(ji), zh_s(ji,jk) ) ! amount of snw that floods, > 0 558 zh_s(ji,jk) = MAX( 0._wp, zh_s(ji,jk) - zdum ) ! remove some snow thickness 559 eh_i_old(ji,0) = eh_i_old(ji,0) + zdum * ze_s(ji,jk) ! 2nd part (snow enthalpy) 560 ! update dh_snowice 561 zdeltah(ji) = MAX( 0._wp, zdeltah(ji) - zdum ) 562 END DO 563 END DO 564 ! 565 ! 566 !!$ ! --- Update snow diags --- ! 567 !!$ !!clem: this is wrong. dh_s_tot is not used anyway 568 !!$ DO ji = 1, npti 569 !!$ dh_s_tot(ji) = dh_s_tot(ji) + dh_s_mlt(ji) + zdeltah(ji) + zdh_s_sub(ji) - dh_snowice(ji) 570 !!$ END DO 571 ! 572 ! 573 ! Remapping of snw enthalpy on a regular grid 574 !-------------------------------------------- 575 CALL snw_ent( zh_s, ze_s, e_s_1d ) 576 577 ! recalculate t_s_1d from e_s_1d 609 578 DO jk = 1, nlay_s 610 579 DO ji = 1,npti 611 ! where there is no ice or no snow 612 rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) ) 613 ! mass & energy loss to the ocean 614 hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * & 615 & ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! heat flux to the ocean [W.m-2], < 0 616 wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * & 617 & ( rhos * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! mass flux 618 ! update energy (mass is updated in the next loop) 619 e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk) 620 ! recalculate t_s_1d from e_s_1d 621 t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 622 END DO 623 END DO 580 IF( h_s_1d(ji) > 0._wp ) THEN 581 t_s_1d(ji,jk) = rt0 + ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi ) 582 ELSE 583 t_s_1d(ji,jk) = rt0 584 ENDIF 585 END DO 586 END DO 587 588 ! Note: remapping of ice enthalpy is done in icethd.F90 624 589 625 590 ! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 --- 626 591 WHERE( h_i_1d(1:npti) == 0._wp ) 627 a_i_1d(1:npti) = 0._wp 628 h_s_1d(1:npti) = 0._wp 592 a_i_1d (1:npti) = 0._wp 593 h_s_1d (1:npti) = 0._wp 594 t_su_1d(1:npti) = rt0 629 595 END WHERE 630 !596 631 597 END SUBROUTINE ice_thd_dh 632 598 599 SUBROUTINE snw_ent( ph_old, pe_old, pe_new ) 600 !!------------------------------------------------------------------- 601 !! *** ROUTINE snw_ent *** 602 !! 603 !! ** Purpose : 604 !! This routine computes new vertical grids in the snow, 605 !! and consistently redistributes temperatures. 606 !! Redistribution is made so as to ensure to energy conservation 607 !! 608 !! 609 !! ** Method : linear conservative remapping 610 !! 611 !! ** Steps : 1) cumulative integrals of old enthalpies/thicknesses 612 !! 2) linear remapping on the new layers 613 !! 614 !! ------------ cum0(0) ------------- cum1(0) 615 !! NEW ------------- 616 !! ------------ cum0(1) ==> ------------- 617 !! ... ------------- 618 !! ------------ ------------- 619 !! ------------ cum0(nlay_s+1) ------------- cum1(nlay_s) 620 !! 621 !! 622 !! References : Bitz & Lipscomb, JGR 99; Vancoppenolle et al., GRL, 2005 623 !!------------------------------------------------------------------- 624 REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: ph_old ! old thicknesses (m) 625 REAL(wp), DIMENSION(jpij,0:nlay_s), INTENT(in ) :: pe_old ! old enthlapies (J.m-3) 626 REAL(wp), DIMENSION(jpij,1:nlay_s), INTENT(inout) :: pe_new ! new enthlapies (J.m-3, remapped) 627 ! 628 INTEGER :: ji ! dummy loop indices 629 INTEGER :: jk0, jk1 ! old/new layer indices 630 ! 631 REAL(wp), DIMENSION(jpij,0:nlay_s+1) :: zeh_cum0, zh_cum0 ! old cumulative enthlapies and layers interfaces 632 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeh_cum1, zh_cum1 ! new cumulative enthlapies and layers interfaces 633 REAL(wp), DIMENSION(jpij) :: zhnew ! new layers thicknesses 634 !!------------------------------------------------------------------- 635 636 !-------------------------------------------------------------------------- 637 ! 1) Cumulative integral of old enthalpy * thickness and layers interfaces 638 !-------------------------------------------------------------------------- 639 zeh_cum0(1:npti,0) = 0._wp 640 zh_cum0 (1:npti,0) = 0._wp 641 DO jk0 = 1, nlay_s+1 642 DO ji = 1, npti 643 zeh_cum0(ji,jk0) = zeh_cum0(ji,jk0-1) + pe_old(ji,jk0-1) * ph_old(ji,jk0-1) 644 zh_cum0 (ji,jk0) = zh_cum0 (ji,jk0-1) + ph_old(ji,jk0-1) 645 END DO 646 END DO 647 648 !------------------------------------ 649 ! 2) Interpolation on the new layers 650 !------------------------------------ 651 ! new layer thickesses 652 DO ji = 1, npti 653 zhnew(ji) = SUM( ph_old(ji,0:nlay_s) ) * r1_nlay_s 654 END DO 655 656 ! new layers interfaces 657 zh_cum1(1:npti,0) = 0._wp 658 DO jk1 = 1, nlay_s 659 DO ji = 1, npti 660 zh_cum1(ji,jk1) = zh_cum1(ji,jk1-1) + zhnew(ji) 661 END DO 662 END DO 663 664 zeh_cum1(1:npti,0:nlay_s) = 0._wp 665 ! new cumulative q*h => linear interpolation 666 DO jk0 = 1, nlay_s+1 667 DO jk1 = 1, nlay_s-1 668 DO ji = 1, npti 669 IF( zh_cum1(ji,jk1) <= zh_cum0(ji,jk0) .AND. zh_cum1(ji,jk1) > zh_cum0(ji,jk0-1) ) THEN 670 zeh_cum1(ji,jk1) = ( zeh_cum0(ji,jk0-1) * ( zh_cum0(ji,jk0) - zh_cum1(ji,jk1 ) ) + & 671 & zeh_cum0(ji,jk0 ) * ( zh_cum1(ji,jk1) - zh_cum0(ji,jk0-1) ) ) & 672 & / ( zh_cum0(ji,jk0) - zh_cum0(ji,jk0-1) ) 673 ENDIF 674 END DO 675 END DO 676 END DO 677 ! to ensure that total heat content is strictly conserved, set: 678 zeh_cum1(1:npti,nlay_s) = zeh_cum0(1:npti,nlay_s+1) 679 680 ! new enthalpies 681 DO jk1 = 1, nlay_s 682 DO ji = 1, npti 683 rswitch = MAX( 0._wp , SIGN( 1._wp , zhnew(ji) - epsi20 ) ) 684 pe_new(ji,jk1) = rswitch * ( zeh_cum1(ji,jk1) - zeh_cum1(ji,jk1-1) ) / MAX( zhnew(ji), epsi20 ) 685 END DO 686 END DO 687 688 END SUBROUTINE snw_ent 689 690 633 691 #else 634 692 !!---------------------------------------------------------------------- -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icethd_pnd.F90
r13899 r14016 20 20 USE ice1D ! sea-ice: thermodynamics variables 21 21 USE icetab ! sea-ice: 1D <==> 2D transformation 22 USE sbc_ice ! surface energy budget 22 23 ! 23 24 USE in_out_manager ! I/O manager 25 USE iom ! I/O manager library 24 26 USE lib_mpp ! MPP library 25 27 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) … … 34 36 INTEGER :: nice_pnd ! choice of the type of pond scheme 35 37 ! ! associated indices: 36 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 37 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 38 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 39 38 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 39 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 40 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 41 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 42 43 !-------------------------------------------------------------------------- 44 ! Diagnostics for pond volume per area 45 ! 46 ! dV/dt = mlt + drn + lid + rnf 47 ! mlt = input from surface melting 48 ! drn = drainage through brine network 49 ! lid = lid growth & melt 50 ! rnf = runoff (water directly removed out of surface melting + overflow) 51 ! 52 ! In topo mode, the pond water lost because it is in the snow is not included in the budget 53 ! In level mode, all terms are incorporated 54 ! 55 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_mlt ! meltwater pond volume input [kg/m2/s] 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_drn ! pond volume lost by drainage [-] 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_lid ! exchange with lid / refreezing [-] 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_rnf ! meltwater pond lost to runoff [-] 59 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_mlt_1d ! meltwater pond volume input [-] 60 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_drn_1d ! pond volume lost by drainage [-] 61 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_lid_1d ! exchange with lid / refreezing [-] 62 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_rnf_1d ! meltwater pond lost to runoff [-] 63 64 !! * Substitutions 65 # include "do_loop_substitute.h90" 40 66 !!---------------------------------------------------------------------- 41 67 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 46 72 47 73 SUBROUTINE ice_thd_pnd 74 48 75 !!------------------------------------------------------------------- 49 76 !! *** ROUTINE ice_thd_pnd *** 50 77 !! 51 78 !! ** Purpose : change melt pond fraction and thickness 52 !! 79 !! 80 !! ** Note : Melt ponds affect only radiative transfer for now 81 !! No heat, no salt. 82 !! The current diagnostics lacks a contribution from drainage 53 83 !!------------------------------------------------------------------- 84 INTEGER :: ji, jj, jl ! loop indices 85 !!------------------------------------------------------------------- 86 87 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 88 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 54 89 ! 55 SELECT CASE ( nice_pnd ) 90 diag_dvpn_mlt (:,:) = 0._wp ; diag_dvpn_drn (:,:) = 0._wp 91 diag_dvpn_lid (:,:) = 0._wp ; diag_dvpn_rnf (:,:) = 0._wp 92 diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp 93 diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp 94 95 !------------------------------------- 96 ! Remove ponds where ice has vanished 97 !------------------------------------- 98 at_i(:,:) = SUM( a_i, dim=3 ) 56 99 ! 57 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 58 ! 59 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 60 ! 61 END SELECT 100 DO jl = 1, jpl 101 DO_2D( 1, 1, 1, 1 ) 102 IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN 103 wfx_pnd (ji,jj) = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 104 a_ip (ji,jj,jl) = 0._wp 105 v_ip (ji,jj,jl) = 0._wp 106 v_il (ji,jj,jl) = 0._wp 107 h_ip (ji,jj,jl) = 0._wp 108 h_il (ji,jj,jl) = 0._wp 109 a_ip_frac(ji,jj,jl) = 0._wp 110 ENDIF 111 END_2D 112 END DO 113 114 !------------------------------ 115 ! Identify grid cells with ice 116 !------------------------------ 117 npti = 0 ; nptidx(:) = 0 118 DO_2D( 1, 1, 1, 1 ) 119 IF( at_i(ji,jj) >= epsi10 ) THEN 120 npti = npti + 1 121 nptidx( npti ) = (jj - 1) * jpi + ji 122 ENDIF 123 END_2D 124 125 !------------------------------------ 126 ! Select melt pond scheme to be used 127 !------------------------------------ 128 IF( npti > 0 ) THEN 129 SELECT CASE ( nice_pnd ) 130 ! 131 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 132 ! 133 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 134 ! 135 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 136 ! 137 END SELECT 138 ENDIF 139 140 !------------------------------------ 141 ! Diagnostics 142 !------------------------------------ 143 CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting 144 CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid 145 CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage 146 CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow 62 147 ! 148 DEALLOCATE( diag_dvpn_mlt , diag_dvpn_lid , diag_dvpn_drn , diag_dvpn_rnf ) 149 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 150 63 151 END SUBROUTINE ice_thd_pnd 64 152 … … 80 168 !! ** References : Bush, G.W., and Trump, D.J. (2017) 81 169 !!------------------------------------------------------------------- 82 INTEGER :: ji ! loop indices 170 INTEGER :: ji, jl ! loop indices 171 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 83 172 !!------------------------------------------------------------------- 84 DO ji = 1, npti 85 ! 86 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 87 h_ip_1d(ji) = rn_hpnd 88 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 89 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 90 ELSE 91 h_ip_1d(ji) = 0._wp 92 a_ip_1d(ji) = 0._wp 93 h_il_1d(ji) = 0._wp 94 ENDIF 95 ! 173 DO jl = 1, jpl 174 175 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 176 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 177 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 179 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 180 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 181 182 DO ji = 1, npti 183 ! 184 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 185 ! 186 IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 187 h_ip_1d(ji) = rn_hpnd 188 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 189 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 190 ELSE 191 h_ip_1d(ji) = 0._wp 192 a_ip_1d(ji) = 0._wp 193 h_il_1d(ji) = 0._wp 194 ENDIF 195 ! 196 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 197 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 198 ! 199 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 200 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 201 ! 202 END DO 203 204 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 205 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 206 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 207 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 208 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 209 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 210 96 211 END DO 97 212 ! … … 132 247 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 133 248 !! 134 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi 135 !! perm = permability of sea-ice 249 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi * flush --- from Flocco et al 2007 --- 250 !! perm = permability of sea-ice + correction from Hunke et al 2012 (flush) 136 251 !! visc = water viscosity 137 252 !! Hp = height of top of the pond above sea-level 138 253 !! Hi = ice thickness thru which there is flushing 254 !! flush= correction otherwise flushing is excessive 139 255 !! 140 256 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness … … 143 259 !! a_ip/a_i = a_ip_frac = h_ip / zaspect 144 260 !! 145 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min261 !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 146 262 !! 147 !! ** Note : mostly stolen from CICE263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 148 264 !! 149 265 !! ** References : Flocco and Feltham (JGR, 2007) 150 266 !! Flocco et al (JGR, 2010) 151 267 !! Holland et al (J. Clim, 2012) 152 !!------------------------------------------------------------------- 268 !! Hunke et al (OM 2012) 269 !!------------------------------------------------------------------- 153 270 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 154 271 !! … … 157 274 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 158 275 !! 159 REAL(wp) :: zfr_mlt, zdv_mlt 276 REAL(wp) :: zfr_mlt, zdv_mlt, zdv_avail ! fraction and volume of available meltwater retained for melt ponding 160 277 REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes 278 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 161 279 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 162 280 REAL(wp) :: zv_ip_max ! max pond volume allowed 163 281 REAL(wp) :: zdT ! zTp-t_su 164 REAL(wp) :: zsbr 282 REAL(wp) :: zsbr, ztmelts ! Brine salinity 165 283 REAL(wp) :: zperm ! permeability of sea ice 166 284 REAL(wp) :: zfac, zdum ! temporary arrays 167 285 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 168 286 !! 169 INTEGER :: ji, jk 287 INTEGER :: ji, jk, jl ! loop indices 170 288 !!------------------------------------------------------------------- 171 289 z1_rhow = 1._wp / rhow 172 290 z1_aspect = 1._wp / zaspect 173 291 z1_Tp = 1._wp / zTp 174 175 DO ji = 1, npti 176 ! !----------------------------------------------------! 177 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 178 ! !----------------------------------------------------! 179 !--- Remove ponds on thin ice or tiny ice fractions 180 a_ip_1d(ji) = 0._wp 181 h_ip_1d(ji) = 0._wp 182 h_il_1d(ji) = 0._wp 183 ! !--------------------------------! 184 ELSE ! Case ice thickness >= rn_himin ! 185 ! !--------------------------------! 186 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 292 293 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti), at_i ) 294 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 295 296 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 297 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 298 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 299 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 300 301 DO jl = 1, jpl 302 303 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 304 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 305 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su(:,:,jl) ) 306 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,jl) ) 307 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip(:,:,jl) ) 308 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il(:,:,jl) ) 309 310 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum(1:npti), dh_i_sum_2d(:,:,jl) ) 311 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt(1:npti), dh_s_mlt_2d(:,:,jl) ) 312 313 DO jk = 1, nlay_i 314 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 315 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) ) 316 END DO 317 318 !----------------------- 319 ! Melt pond calculations 320 !----------------------- 321 DO ji = 1, npti 322 ! 323 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 324 ! !----------------------------------------------------! 325 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 326 ! !----------------------------------------------------! 327 !--- Remove ponds on thin ice or tiny ice fractions 328 a_ip_1d(ji) = 0._wp 329 h_ip_1d(ji) = 0._wp 330 h_il_1d(ji) = 0._wp 331 ! !--------------------------------! 332 ELSE ! Case ice thickness >= rn_himin ! 333 ! !--------------------------------! 334 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 335 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 336 ! 337 !------------------! 338 ! case ice melting ! 339 !------------------! 340 ! 341 !--- available meltwater for melt ponding (zdv_avail) ---! 342 zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 343 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed 344 zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? 345 ! 346 !--- overflow ---! 347 ! 348 ! area driven overflow 349 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 350 ! a_ip_max = zfr_mlt * a_i 351 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 352 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 353 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 354 355 ! depth driven overflow 356 ! If pond depth exceeds half the ice thickness then reduce the pond volume 357 ! h_ip_max = 0.5 * h_i 358 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 359 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 360 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 361 362 !--- Pond growing ---! 363 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 364 ! 365 !--- Lid melting ---! 366 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 367 ! 368 !-------------------! 369 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 370 !-------------------! 371 ! 372 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 373 ! 374 !--- Pond contraction (due to refreezing) ---! 375 IF( ln_pnd_lids ) THEN 376 ! 377 !--- Lid growing and subsequent pond shrinking ---! 378 zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 379 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 380 381 ! Lid growing 382 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 383 384 ! Pond shrinking 385 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 386 387 ELSE 388 zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) 389 ! Pond shrinking 390 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 391 ENDIF 392 ! 393 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 394 ! v_ip = h_ip * a_ip 395 ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 396 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 397 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 398 ! 399 400 !------------------------------------------------! 401 ! Pond drainage through brine network (flushing) ! 402 !------------------------------------------------! 403 ! height of top of the pond above sea-level 404 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 405 406 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 407 DO jk = 1, nlay_i 408 ! MV Assur is inconsistent with SI3 409 !!zsbr = - 1.2_wp & 410 !! & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 411 !! & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 412 !! & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 413 !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 414 ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 415 ztmelts = -rTmlt * sz_i_1d(ji,jk) 416 ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 417 END DO 418 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 419 420 ! Do the drainage using Darcy's law 421 zdv_flush = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012) 422 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 423 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 424 425 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 426 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 427 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 428 429 !--- Corrections and lid thickness ---! 430 IF( ln_pnd_lids ) THEN 431 !--- retrieve lid thickness from volume ---! 432 IF( a_ip_1d(ji) > 0.01_wp ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 433 ELSE ; h_il_1d(ji) = 0._wp 434 ENDIF 435 !--- remove ponds if lids are much larger than ponds ---! 436 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 437 a_ip_1d(ji) = 0._wp 438 h_ip_1d(ji) = 0._wp 439 h_il_1d(ji) = 0._wp 440 ENDIF 441 ENDIF 442 443 ! diagnostics: dvpnd = mlt+rnf+lid+drn 444 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow * zdv_avail * r1_Dt_ice ! > 0, surface melt input 445 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_Dt_ice ! < 0, runoff 446 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + rhow * zdv_frz * r1_Dt_ice ! < 0, shrinking 447 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow * zdv_flush * r1_Dt_ice ! < 0, drainage 448 ! 449 ENDIF 450 ! 451 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 187 452 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 188 453 ! 189 !------------------! 190 ! case ice melting ! 191 !------------------! 454 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 455 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 192 456 ! 193 !--- available meltwater for melt ponding ---! 194 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 195 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) = fraction of melt water that is not flushed 196 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 197 ! 198 !--- overflow ---! 199 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 200 ! a_ip_max = zfr_mlt * a_i 201 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 202 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 203 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 204 205 ! If pond depth exceeds half the ice thickness then reduce the pond volume 206 ! h_ip_max = 0.5 * h_i 207 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 208 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 209 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 210 211 !--- Pond growing ---! 212 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 213 ! 214 !--- Lid melting ---! 215 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 216 ! 217 !--- mass flux ---! 218 IF( zdv_mlt > 0._wp ) THEN 219 zfac = zdv_mlt * rhow * r1_Dt_ice ! melt pond mass flux < 0 [kg.m-2.s-1] 220 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 221 ! 222 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 223 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 224 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 225 ENDIF 226 227 !-------------------! 228 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 229 !-------------------! 230 ! 231 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 232 ! 233 !--- Pond contraction (due to refreezing) ---! 234 IF( ln_pnd_lids ) THEN 235 ! 236 !--- Lid growing and subsequent pond shrinking ---! 237 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 238 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 239 240 ! Lid growing 241 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 242 243 ! Pond shrinking 244 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 245 246 ELSE 247 ! Pond shrinking 248 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 249 ENDIF 250 ! 251 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 252 ! v_ip = h_ip * a_ip 253 ! a_ip/a_i = a_ip_frac = h_ip / zaspect (cf Holland 2012, fitting SHEBA so that knowing v_ip we can distribute it to a_ip and h_ip) 254 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 255 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 256 257 !---------------! 258 ! Pond flushing ! 259 !---------------! 260 ! height of top of the pond above sea-level 261 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 262 263 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 264 DO jk = 1, nlay_i 265 zsbr = - 1.2_wp & 266 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 267 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 268 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 269 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 270 END DO 271 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 272 273 ! Do the drainage using Darcy's law 274 zdv_flush = -zperm * rho0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 275 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 276 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 277 278 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 279 a_ip_1d(ji) = MIN( a_i_1d(ji), SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) ) ! make sure a_ip < a_i 280 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 281 282 !--- Corrections and lid thickness ---! 283 IF( ln_pnd_lids ) THEN 284 !--- retrieve lid thickness from volume ---! 285 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 286 ELSE ; h_il_1d(ji) = 0._wp 287 ENDIF 288 !--- remove ponds if lids are much larger than ponds ---! 289 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 290 a_ip_1d(ji) = 0._wp 291 h_ip_1d(ji) = 0._wp 292 h_il_1d(ji) = 0._wp 293 ENDIF 294 ENDIF 295 ! 296 ENDIF 297 457 END DO 458 459 !-------------------------------------------------------------------- 460 ! Retrieve 2D arrays 461 !-------------------------------------------------------------------- 462 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) ) 466 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) ) 467 ! 298 468 END DO 299 469 ! 470 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd ) 471 ! 472 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 473 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 474 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 475 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 476 ! 300 477 END SUBROUTINE pnd_LEV 301 478 479 480 481 SUBROUTINE pnd_TOPO 482 483 !!------------------------------------------------------------------- 484 !! *** ROUTINE pnd_TOPO *** 485 !! 486 !! ** Purpose : Compute melt pond evolution based on the ice 487 !! topography inferred from the ice thickness distribution 488 !! 489 !! ** Method : This code is initially based on Flocco and Feltham 490 !! (2007) and Flocco et al. (2010). 491 !! 492 !! - Calculate available pond water base on surface meltwater 493 !! - Redistribute water as a function of topography, drain water 494 !! - Exchange water with the lid 495 !! 496 !! ** Tunable parameters : 497 !! 498 !! ** Note : 499 !! 500 !! ** References 501 !! 502 !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 503 !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 504 !! 10.1029/2006JC003836. 505 !! 506 !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 507 !! a physically based melt pond scheme into the sea ice component of a 508 !! climate model. J. Geophys. Res. 115, C08012, 509 !! doi: 10.1029/2009JC005568. 510 !! 511 !!------------------------------------------------------------------- 512 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 513 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 514 zvp_min = 1.e-4_wp ! minimum pond volume (m) 515 516 517 ! local variables 518 REAL(wp) :: & 519 zdHui, & ! change in thickness of ice lid (m) 520 zomega, & ! conduction 521 zdTice, & ! temperature difference across ice lid (C) 522 zdvice, & ! change in ice volume (m) 523 zTavg, & ! mean surface temperature across categories (C) 524 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 525 zTp, & ! pond freezing temperature (C) 526 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 527 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 528 z1_rhow, & ! inverse water density 529 zv_pnd , & ! volume of meltwater contributing to ponds 530 zv_mlt ! total amount of meltwater produced 531 532 REAL(wp), DIMENSION(jpi,jpj) :: zvolp, & !! total melt pond water available before redistribution and drainage 533 zvolp_res !! remaining melt pond water available after drainage 534 535 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 536 537 INTEGER :: ji, jj, jk, jl ! loop indices 538 539 INTEGER :: i_test 540 541 ! Note 542 ! equivalent for CICE translation 543 ! a_ip -> apond 544 ! a_ip_frac -> apnd 545 546 CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 547 548 !--------------------------------------------------------------- 549 ! Initialise 550 !--------------------------------------------------------------- 551 552 ! Parameters & constants (move to parameters) 553 zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) 554 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 555 z1_rhow = 1._wp / rhow 556 557 ! Set required ice variables (hard-coded here for now) 558 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 559 560 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 561 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 562 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 563 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 564 565 ! thickness 566 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 567 ELSEWHERE ; z1_a_i(:,:,:) = 0._wp 568 END WHERE 569 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 570 571 !--------------------------------------------------------------- 572 ! Change 2D to 1D 573 !--------------------------------------------------------------- 574 ! MV 575 ! a less computing-intensive version would have 2D-1D passage here 576 ! use what we have in iceitd.F90 (incremental remapping) 577 578 !-------------------------------------------------------------- 579 ! Collect total available pond water volume 580 !-------------------------------------------------------------- 581 ! Assuming that meltwater (+rain in principle) runsoff the surface 582 ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 583 ! I cite her words, they are very talkative 584 ! "grid cells with very little ice cover (and hence more open water area) 585 ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 586 ! "This results in the same runoff fraction r for each ice category within a grid cell" 587 588 zvolp(:,:) = 0._wp 589 590 DO jl = 1, jpl 591 DO_2D( 1, 1, 1, 1 ) 592 593 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 594 595 !--- Available and contributing meltwater for melt ponding ---! 596 zv_mlt = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) & ! available volume of surface melt water per grid area 597 & * z1_rhow * a_i(ji,jj,jl) 598 ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 599 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds 600 zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl 601 602 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice ! diags 603 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice 604 605 !--- Create possible new ponds 606 ! if pond does not exist, create new pond over full ice area 607 !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 608 IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 609 a_ip(ji,jj,jl) = a_i(ji,jj,jl) 610 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 611 ENDIF 612 613 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 614 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 615 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 616 617 !--- Total available pond water volume (pre-existing + newly produced)j 618 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 619 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 620 621 ENDIF ! a_i 622 623 END_2D 624 END DO ! ji 625 626 !-------------------------------------------------------------- 627 ! Redistribute and drain water from ponds 628 !-------------------------------------------------------------- 629 CALL ice_thd_pnd_area( zvolp, zvolp_res ) 630 631 !-------------------------------------------------------------- 632 ! Melt pond lid growth and melt 633 !-------------------------------------------------------------- 634 635 IF( ln_pnd_lids ) THEN 636 637 DO_2D( 1, 1, 1, 1 ) 638 639 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 640 641 !-------------------------- 642 ! Pond lid growth and melt 643 !-------------------------- 644 ! Mean surface temperature 645 zTavg = 0._wp 646 DO jl = 1, jpl 647 zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 648 END DO 649 zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 650 651 DO jl = 1, jpl-1 652 653 IF ( v_il(ji,jj,jl) > epsi10 ) THEN 654 655 !---------------------------------------------------------------- 656 ! Lid melting: floating upper ice layer melts in whole or part 657 !---------------------------------------------------------------- 658 ! Use Tsfc for each category 659 IF ( t_su(ji,jj,jl) > zTp ) THEN 660 661 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 662 663 IF ( zdvice > epsi10 ) THEN 664 665 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 666 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 667 ! as it is already counted in surface melt 668 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 669 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 670 671 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 672 ! ice lid melted and category is pond covered 673 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 674 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 675 v_il(ji,jj,jl) = 0._wp 676 ENDIF 677 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 678 679 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 680 681 ENDIF 682 683 !---------------------------------------------------------------- 684 ! Freeze pre-existing lid 685 !---------------------------------------------------------------- 686 687 ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 688 689 ! differential growth of base of surface floating ice layer 690 zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 691 zomega = rcnd_i * zdTice / zrhoi_L 692 zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 693 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 694 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 695 696 IF ( zdvice > epsi10 ) THEN 697 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 698 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 699 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 700 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 701 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 702 703 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 704 705 ENDIF 706 707 ENDIF ! Tsfcn(i,j,n) 708 709 !---------------------------------------------------------------- 710 ! Freeze new lids 711 !---------------------------------------------------------------- 712 ! upper ice layer begins to form 713 ! note: albedo does not change 714 715 ELSE ! v_il < epsi10 716 717 ! thickness of newly formed ice 718 ! the surface temperature of a meltpond is the same as that 719 ! of the ice underneath (0C), and the thermodynamic surface 720 ! flux is the same 721 722 !!! we need net surface energy flux, excluding conduction 723 !!! fsurf is summed over categories in CICE 724 !!! we have the category-dependent flux, let us use it ? 725 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 726 zdHui = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp ) 727 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 728 IF ( zdvice > epsi10 ) THEN 729 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 730 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 731 732 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 733 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 734 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 735 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar 736 ENDIF 737 738 ENDIF ! v_il 739 740 END DO ! jl 741 742 ELSE ! remove ponds on thin ice 743 744 v_ip(ji,jj,:) = 0._wp 745 v_il(ji,jj,:) = 0._wp 746 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 747 ! zvolp(ji,jj) = 0._wp 748 749 ENDIF 750 751 END_2D 752 753 ENDIF ! ln_pnd_lids 754 755 !--------------------------------------------------------------- 756 ! Clean-up variables (probably duplicates what icevar would do) 757 !--------------------------------------------------------------- 758 ! MV comment 759 ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 760 ! icevar should recompute all other variables (if needed at all) 761 762 DO jl = 1, jpl 763 764 DO_2D( 1, 1, 1, 1 ) 765 766 ! ! zap lids on small ponds 767 ! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 768 ! .AND. v_il(ji,jj,jl) > epsi10) THEN 769 ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 770 ! ENDIF 771 772 ! recalculate equivalent pond variables 773 IF ( a_ip(ji,jj,jl) > epsi10) THEN 774 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_i(ji,jj,jl) 775 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 776 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 777 ENDIF 778 ! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 779 ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 780 ! ENDIF 781 782 END_2D 783 784 END DO ! jl 785 786 787 END SUBROUTINE pnd_TOPO 788 789 790 SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 791 792 !!------------------------------------------------------------------- 793 !! *** ROUTINE ice_thd_pnd_area *** 794 !! 795 !! ** Purpose : Given the total volume of available pond water, 796 !! redistribute and drain water 797 !! 798 !! ** Method 799 !! 800 !-----------| 801 ! | 802 ! |-----------| 803 !___________|___________|______________________________________sea-level 804 ! | | 805 ! | |---^--------| 806 ! | | | | 807 ! | | | |-----------| |------- 808 ! | | | alfan | | | 809 ! | | | | |--------------| 810 ! | | | | | | 811 !---------------------------v------------------------------------------- 812 ! | | ^ | | | 813 ! | | | | |--------------| 814 ! | | | betan | | | 815 ! | | | |-----------| |------- 816 ! | | | | 817 ! | |---v------- | 818 ! | | 819 ! |-----------| 820 ! | 821 !-----------| 822 ! 823 !! 824 !!------------------------------------------------------------------ 825 826 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 827 zvolp, & ! total available pond water 828 zdvolp ! remaining meltwater after redistribution 829 830 INTEGER :: & 831 ns, & 832 m_index, & 833 permflag 834 835 REAL (wp), DIMENSION(jpl) :: & 836 hicen, & 837 hsnon, & 838 asnon, & 839 alfan, & 840 betan, & 841 cum_max_vol, & 842 reduced_aicen 843 844 REAL (wp), DIMENSION(0:jpl) :: & 845 cum_max_vol_tmp 846 847 REAL (wp) :: & 848 hpond, & 849 drain, & 850 floe_weight, & 851 pressure_head, & 852 hsl_rel, & 853 deltah, & 854 perm, & 855 msno 856 857 REAL (wp), parameter :: & 858 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 859 860 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 861 zvp_min = 1.e-4_wp ! minimum pond volume (m) 862 863 INTEGER :: ji, jj, jk, jl ! loop indices 864 865 a_ip(:,:,:) = 0._wp 866 h_ip(:,:,:) = 0._wp 867 868 DO_2D( 1, 1, 1, 1 ) 869 870 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 871 872 !------------------------------------------------------------------- 873 ! initialize 874 !------------------------------------------------------------------- 875 876 DO jl = 1, jpl 877 878 !---------------------------------------- 879 ! compute the effective snow fraction 880 !---------------------------------------- 881 882 IF (a_i(ji,jj,jl) < epsi10) THEN 883 hicen(jl) = 0._wp 884 hsnon(jl) = 0._wp 885 reduced_aicen(jl) = 0._wp 886 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 887 ELSE 888 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 889 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 890 reduced_aicen(jl) = 1._wp ! n=jpl 891 892 !js: initial code in NEMO_DEV 893 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 894 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 895 896 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 897 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 898 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 899 900 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 901 ! MV should check whether this makes sense to have the same effective snow fraction in here 902 ! OLI: it probably doesn't 903 END IF 904 905 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 906 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 907 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 908 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 909 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 910 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 911 ! alfan = 60% of the ice volume) in each category lies above the reference line, 912 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 913 914 ! MV: 915 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 916 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 917 918 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 919 920 alfan(jl) = 0.6 * hicen(jl) 921 betan(jl) = 0.4 * hicen(jl) 922 923 cum_max_vol(jl) = 0._wp 924 cum_max_vol_tmp(jl) = 0._wp 925 926 END DO ! jpl 927 928 cum_max_vol_tmp(0) = 0._wp 929 drain = 0._wp 930 zdvolp(ji,jj) = 0._wp 931 932 !---------------------------------------------------------- 933 ! Drain overflow water, update pond fraction and volume 934 !---------------------------------------------------------- 935 936 !-------------------------------------------------------------------------- 937 ! the maximum amount of water that can be contained up to each ice category 938 !-------------------------------------------------------------------------- 939 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 940 ! Then the excess volume cum_max_vol(jl) drains out of the system 941 ! It should be added to wfx_pnd_out 942 943 DO jl = 1, jpl-1 ! last category can not hold any volume 944 945 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 946 947 ! total volume in level including snow 948 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 949 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 950 951 ! subtract snow solid volumes from lower categories in current level 952 DO ns = 1, jl 953 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 954 - rhos/rhow * & ! free air fraction that can be filled by water 955 asnon(ns) * & ! effective areal fraction of snow in that category 956 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 957 END DO 958 959 ELSE ! assume higher categories unoccupied 960 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 961 END IF 962 !IF (cum_max_vol_tmp(jl) < z0) THEN 963 ! CALL abort_ice('negative melt pond volume') 964 !END IF 965 END DO 966 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 967 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 968 969 !---------------------------------------------------------------- 970 ! is there more meltwater than can be held in the floe? 971 !---------------------------------------------------------------- 972 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 973 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 974 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 975 976 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 977 978 zdvolp(ji,jj) = drain ! this is the drained water 979 IF (zvolp(ji,jj) < epsi10) THEN 980 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 981 zvolp(ji,jj) = 0._wp 982 END IF 983 END IF 984 985 ! height and area corresponding to the remaining volume 986 ! routine leaves zvolp unchanged 987 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 988 989 DO jl = 1, m_index 990 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 991 ! ! volume instead, no ? 992 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 993 a_ip(ji,jj,jl) = reduced_aicen(jl) 994 ! in practise, pond fraction depends on the empirical snow fraction 995 ! so in turn on ice thickness 996 END DO 997 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 998 999 !------------------------------------------------------------------------ 1000 ! Drainage through brine network (permeability) 1001 !------------------------------------------------------------------------ 1002 !!! drainage due to ice permeability - Darcy's law 1003 1004 ! sea water level 1005 msno = 0._wp 1006 DO jl = 1 , jpl 1007 msno = msno + v_s(ji,jj,jl) * rhos 1008 END DO 1009 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj) 1010 hsl_rel = floe_weight / rho0 & 1011 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1012 1013 deltah = hpond - hsl_rel 1014 pressure_head = grav * rho0 * max(deltah, 0._wp) 1015 1016 ! drain if ice is permeable 1017 permflag = 0 1018 1019 IF (pressure_head > 0._wp) THEN 1020 DO jl = 1, jpl-1 1021 IF ( hicen(jl) /= 0._wp ) THEN 1022 1023 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1024 1025 perm = 0._wp ! MV ugly dummy patch 1026 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1027 IF (perm > 0._wp) permflag = 1 1028 1029 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / & 1030 (viscosity*hicen(jl)) 1031 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1032 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1033 1034 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1035 1036 IF (zvolp(ji,jj) < epsi10) THEN 1037 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 1038 zvolp(ji,jj) = 0._wp 1039 END IF 1040 END IF 1041 END DO 1042 1043 ! adjust melt pond dimensions 1044 IF (permflag > 0) THEN 1045 ! recompute pond depth 1046 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1047 DO jl = 1, m_index 1048 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 1049 a_ip(ji,jj,jl) = reduced_aicen(jl) 1050 END DO 1051 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1052 END IF 1053 END IF ! pressure_head 1054 1055 !------------------------------- 1056 ! remove water from the snow 1057 !------------------------------- 1058 !------------------------------------------------------------------------ 1059 ! total melt pond volume in category does not include snow volume 1060 ! snow in melt ponds is not melted 1061 !------------------------------------------------------------------------ 1062 1063 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1064 ! how much, so I did not diagnose it 1065 ! so if there is a problem here, nobody is going to see it... 1066 1067 1068 ! Calculate pond volume for lower categories 1069 DO jl = 1,m_index-1 1070 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 1071 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1072 END DO 1073 1074 ! Calculate pond volume for highest category = remaining pond volume 1075 1076 ! The following is completely unclear to Martin at least 1077 ! Could we redefine properly and recode in a more readable way ? 1078 1079 ! m_index = last category with melt pond 1080 1081 IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 1082 1083 IF (m_index > 1) THEN 1084 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1085 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 1086 ELSE 1087 v_ip(ji,jj,m_index) = 0._wp 1088 h_ip(ji,jj,m_index) = 0._wp 1089 a_ip(ji,jj,m_index) = 0._wp 1090 ! If remaining pond volume is negative reduce pond volume of 1091 ! lower category 1092 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 1093 v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 1094 END IF 1095 END IF 1096 1097 DO jl = 1,m_index 1098 IF (a_ip(ji,jj,jl) > epsi10) THEN 1099 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 1100 ELSE 1101 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1102 h_ip(ji,jj,jl) = 0._wp 1103 v_ip(ji,jj,jl) = 0._wp 1104 a_ip(ji,jj,jl) = 0._wp 1105 END IF 1106 END DO 1107 DO jl = m_index+1, jpl 1108 h_ip(ji,jj,jl) = 0._wp 1109 a_ip(ji,jj,jl) = 0._wp 1110 v_ip(ji,jj,jl) = 0._wp 1111 END DO 1112 1113 ENDIF 1114 1115 END_2D 1116 1117 END SUBROUTINE ice_thd_pnd_area 1118 1119 1120 SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 1121 !!------------------------------------------------------------------- 1122 !! *** ROUTINE ice_thd_pnd_depth *** 1123 !! 1124 !! ** Purpose : Compute melt pond depth 1125 !!------------------------------------------------------------------- 1126 1127 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1128 aicen, & 1129 asnon, & 1130 hsnon, & 1131 alfan, & 1132 cum_max_vol 1133 1134 REAL (wp), INTENT(IN) :: & 1135 zvolp 1136 1137 REAL (wp), INTENT(OUT) :: & 1138 hpond 1139 1140 INTEGER, INTENT(OUT) :: & 1141 m_index 1142 1143 INTEGER :: n, ns 1144 1145 REAL (wp), DIMENSION(0:jpl+1) :: & 1146 hitl, & 1147 aicetl 1148 1149 REAL (wp) :: & 1150 rem_vol, & 1151 area, & 1152 vol, & 1153 tmp, & 1154 z0 = 0.0_wp 1155 1156 !---------------------------------------------------------------- 1157 ! hpond is zero if zvolp is zero - have we fully drained? 1158 !---------------------------------------------------------------- 1159 1160 IF (zvolp < epsi10) THEN 1161 hpond = z0 1162 m_index = 0 1163 ELSE 1164 1165 !---------------------------------------------------------------- 1166 ! Calculate the category where water fills up to 1167 !---------------------------------------------------------------- 1168 1169 !----------| 1170 ! | 1171 ! | 1172 ! |----------| -- -- 1173 !__________|__________|_________________________________________ ^ 1174 ! | | rem_vol ^ | Semi-filled 1175 ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1176 ! | | | | 1177 ! | | | |hpond 1178 ! | | |----------| | |------- 1179 ! | | | | | | 1180 ! | | | |---v-----| 1181 ! | | m_index | | | 1182 !------------------------------------------------------------- 1183 1184 m_index = 0 ! 1:m_index categories have water in them 1185 DO n = 1, jpl 1186 IF (zvolp <= cum_max_vol(n)) THEN 1187 m_index = n 1188 IF (n == 1) THEN 1189 rem_vol = zvolp 1190 ELSE 1191 rem_vol = zvolp - cum_max_vol(n-1) 1192 END IF 1193 exit ! to break out of the loop 1194 END IF 1195 END DO 1196 m_index = min(jpl-1, m_index) 1197 1198 !---------------------------------------------------------------- 1199 ! semi-filled layer may have m_index different snow in it 1200 !---------------------------------------------------------------- 1201 1202 !----------------------------------------------------------- ^ 1203 ! | alfan(m_index+1) 1204 ! | 1205 !hitl(3)--> |----------| | 1206 !hitl(2)--> |------------| * * * * *| | 1207 !hitl(1)--> |----------|* * * * * * |* * * * * | | 1208 !hitl(0)-->------------------------------------------------- | ^ 1209 ! various snow from lower categories | |alfa(m_index) 1210 1211 ! hitl - heights of the snow layers from thinner and current categories 1212 ! aicetl - area of each snow depth in this layer 1213 1214 hitl(:) = z0 1215 aicetl(:) = z0 1216 DO n = 1, m_index 1217 hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1218 alfan(m_index+1) - alfan(m_index)), z0) 1219 aicetl(n) = asnon(n) 1220 1221 aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1222 END DO 1223 1224 hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1225 aicetl(m_index+1) = z0 1226 1227 !---------------------------------------------------------------- 1228 ! reorder array according to hitl 1229 ! snow heights not necessarily in height order 1230 !---------------------------------------------------------------- 1231 1232 DO ns = 1, m_index+1 1233 DO n = 0, m_index - ns + 1 1234 IF (hitl(n) > hitl(n+1)) THEN ! swap order 1235 tmp = hitl(n) 1236 hitl(n) = hitl(n+1) 1237 hitl(n+1) = tmp 1238 tmp = aicetl(n) 1239 aicetl(n) = aicetl(n+1) 1240 aicetl(n+1) = tmp 1241 END IF 1242 END DO 1243 END DO 1244 1245 !---------------------------------------------------------------- 1246 ! divide semi-filled layer into set of sublayers each vertically homogenous 1247 !---------------------------------------------------------------- 1248 1249 !hitl(3)---------------------------------------------------------------- 1250 ! | * * * * * * * * 1251 ! |* * * * * * * * * 1252 !hitl(2)---------------------------------------------------------------- 1253 ! | * * * * * * * * | * * * * * * * * 1254 ! |* * * * * * * * * |* * * * * * * * * 1255 !hitl(1)---------------------------------------------------------------- 1256 ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1257 ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1258 !hitl(0)---------------------------------------------------------------- 1259 ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1260 1261 ! move up over layers incrementing volume 1262 DO n = 1, m_index+1 1263 1264 area = sum(aicetl(:)) - & ! total area of sub-layer 1265 (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1266 1267 vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1268 1269 IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1270 hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 1271 1272 exit 1273 ELSE ! still in sub-layer below the sub-layer with the depth 1274 rem_vol = rem_vol - vol 1275 END IF 1276 1277 END DO 1278 1279 END IF 1280 1281 END SUBROUTINE ice_thd_pnd_depth 1282 1283 1284 SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 1285 !!------------------------------------------------------------------- 1286 !! *** ROUTINE ice_thd_pnd_perm *** 1287 !! 1288 !! ** Purpose : Determine the liquid fraction of brine in the ice 1289 !! and its permeability 1290 !!------------------------------------------------------------------- 1291 1292 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1293 ticen, & ! internal ice temperature (K) 1294 salin ! salinity (ppt) !js: ppt according to cice 1295 1296 REAL (wp), INTENT(OUT) :: & 1297 perm ! permeability 1298 1299 REAL (wp) :: & 1300 Sbr ! brine salinity 1301 1302 REAL (wp), DIMENSION(nlay_i) :: & 1303 Tin, & ! ice temperature 1304 phi ! liquid fraction 1305 1306 INTEGER :: k 1307 1308 !----------------------------------------------------------------- 1309 ! Compute ice temperatures from enthalpies using quadratic formula 1310 !----------------------------------------------------------------- 1311 1312 DO k = 1,nlay_i 1313 Tin(k) = ticen(k) - rt0 !js: from K to degC 1314 END DO 1315 1316 !----------------------------------------------------------------- 1317 ! brine salinity and liquid fraction 1318 !----------------------------------------------------------------- 1319 1320 DO k = 1, nlay_i 1321 1322 Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 1323 ! Best expression to date is that one (Vancoppenolle et al JGR 2019) 1324 ! Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 1325 phi(k) = salin(k) / Sbr 1326 1327 END DO 1328 1329 !----------------------------------------------------------------- 1330 ! permeability 1331 !----------------------------------------------------------------- 1332 1333 perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 1334 1335 END SUBROUTINE ice_thd_pnd_perm 302 1336 303 1337 SUBROUTINE ice_thd_pnd_init … … 315 1349 INTEGER :: ios, ioptio ! Local integer 316 1350 !! 317 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &1351 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 318 1352 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1353 & ln_pnd_TOPO, & 319 1354 & ln_pnd_lids, ln_pnd_alb 320 1355 !!------------------------------------------------------------------- … … 332 1367 WRITE(numout,*) ' Namelist namicethd_pnd:' 333 1368 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 1369 WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO 334 1370 WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV 335 1371 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 336 1372 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 1373 WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush 337 1374 WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST 338 1375 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd … … 347 1384 IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF 348 1385 IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF 1386 IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF 349 1387 IF( ioptio /= 1 ) & 350 & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' )1388 & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV, ln_pnd_CST or ln_pnd_TOPO)' ) 351 1389 ! 352 1390 SELECT CASE( nice_pnd ) -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icethd_zdf_bl99.F90
r13899 r14016 109 109 REAL(wp), DIMENSION(jpij) :: zdqns_ice_b ! derivative of the surface flux function 110 110 ! 111 REAL(wp), DIMENSION(jpij ) 112 REAL(wp), DIMENSION(jpij,nlay_i) 113 REAL(wp), DIMENSION(jpij,nlay_s) 114 REAL(wp), DIMENSION(jpij,nlay_i) 115 REAL(wp), DIMENSION(jpij,nlay_s) 116 REAL(wp), DIMENSION(jpij,0:nlay_i) 117 REAL(wp), DIMENSION(jpij,0:nlay_i) 118 REAL(wp), DIMENSION(jpij,0:nlay_i) 119 REAL(wp), DIMENSION(jpij,0:nlay_i) 120 REAL(wp), DIMENSION(jpij,0:nlay_i) 121 REAL(wp), DIMENSION(jpij,0:nlay_i) 122 REAL(wp), DIMENSION(jpij,0:nlay_s) 123 REAL(wp), DIMENSION(jpij,0:nlay_s) 124 REAL(wp), DIMENSION(jpij,0:nlay_s) 125 REAL(wp), DIMENSION(jpij,0:nlay_s) 126 REAL(wp), DIMENSION(jpij) 127 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindterm ! 'Ind'ependent term128 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zindtbis ! Temporary 'ind'ependent term129 REAL(wp), DIMENSION(jpij ,nlay_i+3) :: zdiagbis ! Temporary 'dia'gonal term130 REAL(wp), DIMENSION(jpij ,nlay_i+3,3) :: ztrid ! Tridiagonal system terms131 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat132 REAL(wp), DIMENSION(jpij ) :: zghe ! G(he), th. conduct enhancement factor, mono-cat133 REAL(wp), DIMENSION(jpij ) :: za_s_fra ! ice fraction covered by snow134 REAL(wp), DIMENSION(jpij ) :: isnow ! snow presence (1) or not (0)135 REAL(wp), DIMENSION(jpij ) :: isnow_comb ! snow presence for met-office111 REAL(wp), DIMENSION(jpij ) :: ztsuold ! Old surface temperature in the ice 112 REAL(wp), DIMENSION(jpij,nlay_i) :: ztiold ! Old temperature in the ice 113 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsold ! Old temperature in the snow 114 REAL(wp), DIMENSION(jpij,nlay_i) :: ztib ! Temporary temperature in the ice to check the convergence 115 REAL(wp), DIMENSION(jpij,nlay_s) :: ztsb ! Temporary temperature in the snow to check the convergence 116 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i ! Ice thermal conductivity 117 REAL(wp), DIMENSION(jpij,0:nlay_i) :: ztcond_i_cp ! copy 118 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradtr_i ! Radiation transmitted through the ice 119 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zradab_i ! Radiation absorbed in the ice 120 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zkappa_i ! Kappa factor in the ice 121 REAL(wp), DIMENSION(jpij,0:nlay_i) :: zeta_i ! Eta factor in the ice 122 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradtr_s ! Radiation transmited through the snow 123 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zradab_s ! Radiation absorbed in the snow 124 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zkappa_s ! Kappa factor in the snow 125 REAL(wp), DIMENSION(jpij,0:nlay_s) :: zeta_s ! Eta factor in the snow 126 REAL(wp), DIMENSION(jpij) :: zkappa_comb ! Combined snow and ice surface conductivity 127 REAL(wp), DIMENSION(jpij) :: zq_ini ! diag errors on heat 128 REAL(wp), DIMENSION(jpij) :: zghe ! G(he), th. conduct enhancement factor, mono-cat 129 REAL(wp), DIMENSION(jpij) :: za_s_fra ! ice fraction covered by snow 130 REAL(wp), DIMENSION(jpij) :: isnow ! snow presence (1) or not (0) 131 REAL(wp), DIMENSION(jpij) :: isnow_comb ! snow presence for met-office 132 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindterm ! 'Ind'ependent term 133 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zindtbis ! Temporary 'ind'ependent term 134 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1) :: zdiagbis ! Temporary 'dia'gonal term 135 REAL(wp), DIMENSION(jpij,nlay_i+nlay_s+1,3) :: ztrid ! Tridiagonal system terms 136 136 ! 137 137 ! Mono-category … … 533 533 ! Solve the tridiagonal system with Gauss elimination method. 534 534 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 535 jm_maxt = 0 536 jm_mint = nlay_i+5 537 DO ji = 1, npti 538 jm_mint = MIN(jm_min(ji),jm_mint) 539 jm_maxt = MAX(jm_max(ji),jm_maxt) 540 END DO 541 542 DO jk = jm_mint+1, jm_maxt 543 DO ji = 1, npti 544 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 535 !!$ jm_maxt = 0 536 !!$ jm_mint = nlay_i+5 537 !!$ DO ji = 1, npti 538 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 539 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 540 !!$ END DO 541 !!$ !!clem SNWLAY => check why LIM1D does not get this loop. Is nlay_i+5 correct? 542 !!$ 543 !!$ DO jk = jm_mint+1, jm_maxt 544 !!$ DO ji = 1, npti 545 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 546 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 547 !!$ zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 548 !!$ END DO 549 !!$ END DO 550 ! clem: maybe one should find a way to reverse this loop for mpi performance 551 DO ji = 1, npti 552 jm_mint = jm_min(ji) 553 jm_maxt = jm_max(ji) 554 DO jm = jm_mint+1, jm_maxt 545 555 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 546 556 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) … … 564 574 END DO 565 575 576 ! snow temperatures 566 577 DO ji = 1, npti 567 578 ! Variables used after iterations 568 579 ! Value must be frozen after convergence for MPP independance reason 569 IF ( .NOT. l_T_converged(ji) ) THEN 570 ! snow temperatures 571 IF( h_s_1d(ji) > 0._wp ) THEN 572 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 573 ENDIF 574 ! surface temperature 580 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 581 & t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 582 END DO 583 !!clem SNWLAY 584 DO jm = nlay_s, 2, -1 585 DO ji = 1, npti 586 jk = jm - 1 587 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 588 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 589 END DO 590 END DO 591 592 ! surface temperature 593 DO ji = 1, npti 594 IF( .NOT. l_T_converged(ji) ) THEN 575 595 ztsub(ji) = t_su_1d(ji) 576 596 IF( t_su_1d(ji) < rt0 ) THEN 577 t_su_1d(ji) = ( 578 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) *t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji))597 t_su_1d(ji) = ( zindtbis(ji,jm_min(ji)) - ztrid(ji,jm_min(ji),3) * & 598 & ( isnow(ji) * t_s_1d(ji,1) + ( 1._wp - isnow(ji) ) * t_i_1d(ji,1) ) ) / zdiagbis(ji,jm_min(ji)) 579 599 ENDIF 580 600 ENDIF 581 601 END DO 582 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1)583 602 ! 584 603 !-------------------------------------------------------------- … … 727 746 ! Solve the tridiagonal system with Gauss elimination method. 728 747 ! Thomas algorithm, from Computational fluid Dynamics, J.D. ANDERSON, McGraw-Hill 1984 729 jm_maxt = 0 730 jm_mint = nlay_i+5 731 DO ji = 1, npti 732 jm_mint = MIN(jm_min(ji),jm_mint) 733 jm_maxt = MAX(jm_max(ji),jm_maxt) 734 END DO 735 736 DO jk = jm_mint+1, jm_maxt 737 DO ji = 1, npti 738 jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 748 !!$ jm_maxt = 0 749 !!$ jm_mint = nlay_i+5 750 !!$ DO ji = 1, npti 751 !!$ jm_mint = MIN(jm_min(ji),jm_mint) 752 !!$ jm_maxt = MAX(jm_max(ji),jm_maxt) 753 !!$ END DO 754 !!$ 755 !!$ DO jk = jm_mint+1, jm_maxt 756 !!$ DO ji = 1, npti 757 !!$ jm = MIN(MAX(jm_min(ji)+1,jk),jm_max(ji)) 758 !!$ zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 759 !!$ zindtbis(ji,jm) = zindterm(ji,jm) - ztrid(ji,jm,1) * zindtbis(ji,jm-1) / zdiagbis(ji,jm-1) 760 !!$ END DO 761 !!$ END DO 762 ! clem: maybe one should find a way to reverse this loop for mpi performance 763 DO ji = 1, npti 764 jm_mint = jm_min(ji) 765 jm_maxt = jm_max(ji) 766 DO jm = jm_mint+1, jm_maxt 739 767 zdiagbis(ji,jm) = ztrid (ji,jm,2) - ztrid(ji,jm,1) * ztrid (ji,jm-1,3) / zdiagbis(ji,jm-1) 740 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1)/ zdiagbis(ji,jm-1)768 zindtbis(ji,jm) = zindterm(ji,jm ) - ztrid(ji,jm,1) * zindtbis(ji,jm-1 ) / zdiagbis(ji,jm-1) 741 769 END DO 742 770 END DO 743 771 744 772 ! ice temperatures 745 773 DO ji = 1, npti … … 761 789 ! snow temperatures 762 790 DO ji = 1, npti 763 ! Variable used after iterations791 ! Variables used after iterations 764 792 ! Value must be frozen after convergence for MPP independance reason 765 IF ( .NOT. l_T_converged(ji) ) THEN 766 IF( h_s_1d(ji) > 0._wp ) THEN 767 t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 768 ENDIF 769 ENDIF 770 END DO 771 !clem: in order to have several layers of snow, there is a missing loop here for t_s_1d(1:nlay_s-1) 793 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 794 & t_s_1d(ji,nlay_s) = ( zindtbis(ji,nlay_s+1) - ztrid(ji,nlay_s+1,3) * t_i_1d(ji,1) ) / zdiagbis(ji,nlay_s+1) 795 END DO 796 !!clem SNWLAY 797 DO jm = nlay_s, 2, -1 798 DO ji = 1, npti 799 jk = jm - 1 800 IF ( .NOT. l_T_converged(ji) .AND. h_s_1d(ji) > 0._wp ) & 801 & t_s_1d(ji,jk) = ( zindtbis(ji,jm) - ztrid(ji,jm,3) * t_s_1d(ji,jk+1) ) / zdiagbis(ji,jm) 802 END DO 803 END DO 772 804 ! 773 805 !-------------------------------------------------------------- … … 923 955 !--- Snow-ice interfacial temperature (diagnostic SIMIP) 924 956 IF( h_s_1d(ji) >= zhs_ssl ) THEN 925 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji, 1) &926 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) &957 t_si_1d(ji) = ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i * t_s_1d(ji,nlay_s) & 958 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s * t_i_1d(ji,1) ) & 927 959 & / ( rn_cnd_s * h_i_1d(ji) * r1_nlay_i & 928 960 & + ztcond_i(ji,1) * h_s_1d(ji) * r1_nlay_s ) -
NEMO/branches/2020/tickets_icb_1900/src/ICE/iceupdate.F90
r13899 r14016 94 94 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 95 95 !!--------------------------------------------------------------------- 96 IF( ln_timing ) CALL timing_start('ice _update')96 IF( ln_timing ) CALL timing_start('iceupdate') 97 97 98 98 IF( kt == nit000 .AND. lwp ) THEN … … 154 154 ! ice-ocean mass flux 155 155 wfx_ice(ji,jj) = wfx_bog(ji,jj) + wfx_bom(ji,jj) + wfx_sum(ji,jj) + wfx_sni(ji,jj) & 156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) + wfx_pnd(ji,jj)156 & + wfx_opw(ji,jj) + wfx_dyn(ji,jj) + wfx_res(ji,jj) + wfx_lam(ji,jj) 157 157 158 158 ! snw-ocean mass flux … … 160 160 161 161 ! total mass flux at the ocean/ice interface 162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_ err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_ err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux162 fmmflx(ji,jj) = - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! ice-ocean mass flux saved at least for biogeochemical model 163 emp (ji,jj) = emp_oce(ji,jj) - wfx_ice(ji,jj) - wfx_snw(ji,jj) - wfx_pnd(ji,jj) - wfx_err_sub(ji,jj) ! atm-ocean + ice-ocean mass flux 164 164 165 165 ! Salt flux at the ocean surface … … 172 172 snwice_mass_b(ji,jj) = snwice_mass(ji,jj) ! save mass from the previous ice time step 173 173 ! ! new mass per unit area 174 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) )174 snwice_mass (ji,jj) = tmask(ji,jj,1) * ( rhos * vt_s(ji,jj) + rhoi * vt_i(ji,jj) + rhow * (vt_ip(ji,jj) + vt_il(ji,jj)) ) 175 175 ! ! time evolution of snow+ice mass 176 176 snwice_fmass (ji,jj) = ( snwice_mass(ji,jj) - snwice_mass_b(ji,jj) ) * r1_Dt_ice … … 286 286 IF( ln_icectl ) CALL ice_prt (kt, iiceprt, jiceprt, 3, 'Final state ice_update') ! prints 287 287 IF( sn_cfctl%l_prtctl ) CALL ice_prt3D ('iceupdate') ! prints 288 IF( ln_timing ) CALL timing_stop ('ice _update')! timing288 IF( ln_timing ) CALL timing_stop ('iceupdate') ! timing 289 289 ! 290 290 END SUBROUTINE ice_update_flx … … 324 324 REAL(wp) :: zflagi ! - - 325 325 !!--------------------------------------------------------------------- 326 IF( ln_timing ) CALL timing_start('ice_update _tau')326 IF( ln_timing ) CALL timing_start('ice_update') 327 327 328 328 IF( kt == nit000 .AND. lwp ) THEN … … 376 376 CALL lbc_lnk_multi( 'iceupdate', utau, 'U', -1.0_wp, vtau, 'V', -1.0_wp ) ! lateral boundary condition 377 377 ! 378 IF( ln_timing ) CALL timing_stop('ice_update _tau')378 IF( ln_timing ) CALL timing_stop('ice_update') 379 379 ! 380 380 END SUBROUTINE ice_update_tau -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icevar.F90
r13899 r14016 236 236 z1_zhmax = 1._wp / hi_max(jpl) 237 237 WHERE( h_i(:,:,jpl) > zhmax ) ! bound h_i by hi_max (i.e. 99 m) with associated update of ice area 238 h_i (:,:,jpl) = zhmax238 h_i (:,:,jpl) = zhmax 239 239 a_i (:,:,jpl) = v_i(:,:,jpl) * z1_zhmax 240 240 z1_a_i(:,:,jpl) = zhmax * z1_v_i(:,:,jpl) … … 252 252 ELSEWHERE( h_il(:,:,:) >= zhl_max ) ; a_ip_eff(:,:,:) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow 253 253 ELSEWHERE ; a_ip_eff(:,:,:) = a_ip_frac(:,:,:) * & ! lid is in between. Expose part of the pond 254 & ( h_il(:,:,:) - zhl_min) / ( zhl_max - zhl_min )254 & ( zhl_max - h_il(:,:,:) ) / ( zhl_max - zhl_min ) 255 255 END WHERE 256 256 ! … … 534 534 DO_2D( 1, 1, 1, 1 ) 535 535 ! update exchanges with ocean 536 sfx_res(ji,jj) = sfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice 537 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice 538 wfx_res(ji,jj) = wfx_res(ji,jj) + (1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice 536 sfx_res(ji,jj) = sfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * sv_i(ji,jj,jl) * rhoi * r1_Dt_ice 537 wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_i (ji,jj,jl) * rhoi * r1_Dt_ice 538 wfx_res(ji,jj) = wfx_res(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * v_s (ji,jj,jl) * rhos * r1_Dt_ice 539 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + ( 1._wp - zswitch(ji,jj) ) * ( v_ip(ji,jj,jl)+v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 539 540 ! 540 541 a_i (ji,jj,jl) = a_i (ji,jj,jl) * zswitch(ji,jj) … … 551 552 v_ip (ji,jj,jl) = v_ip (ji,jj,jl) * zswitch(ji,jj) 552 553 v_il (ji,jj,jl) = v_il (ji,jj,jl) * zswitch(ji,jj) 554 h_ip (ji,jj,jl) = h_ip (ji,jj,jl) * zswitch(ji,jj) 555 h_il (ji,jj,jl) = h_il (ji,jj,jl) * zswitch(ji,jj) 553 556 ! 554 557 END_2D … … 635 638 psv_i (ji,jj,jl) = 0._wp 636 639 ENDIF 640 IF( pv_ip(ji,jj,jl) < 0._wp .OR. pv_il(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 641 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + pv_il(ji,jj,jl) * rhow * z1_dt 642 pv_il (ji,jj,jl) = 0._wp 643 ENDIF 644 IF( pv_ip(ji,jj,jl) < 0._wp .OR. pa_ip(ji,jj,jl) <= 0._wp ) THEN 645 wfx_pnd(ji,jj) = wfx_pnd(ji,jj) + pv_ip(ji,jj,jl) * rhow * z1_dt 646 pv_ip (ji,jj,jl) = 0._wp 647 ENDIF 637 648 END_2D 638 649 ! … … 643 654 WHERE( pa_i (:,:,:) < 0._wp ) pa_i (:,:,:) = 0._wp 644 655 WHERE( pa_ip (:,:,:) < 0._wp ) pa_ip (:,:,:) = 0._wp 645 WHERE( pv_ip (:,:,:) < 0._wp ) pv_ip (:,:,:) = 0._wp ! in theory one should change wfx_pnd(-) and wfx_sum(+)646 WHERE( pv_il (:,:,:) < 0._wp ) pv_il (:,:,:) = 0._wp ! but it does not change conservation, so keep it this way is ok647 656 ! 648 657 END SUBROUTINE ice_var_zapneg … … 675 684 WHERE( pe_i (1:npti,:,:) < 0._wp ) pe_i (1:npti,:,:) = 0._wp ! e_i must be >= 0 676 685 WHERE( pe_s (1:npti,:,:) < 0._wp ) pe_s (1:npti,:,:) = 0._wp ! e_s must be >= 0 677 IF( ln_pnd_LEV ) THEN686 IF( ln_pnd_LEV .OR. ln_pnd_TOPO ) THEN 678 687 WHERE( pa_ip(1:npti,:) < 0._wp ) pa_ip(1:npti,:) = 0._wp ! a_ip must be >= 0 679 688 WHERE( pv_ip(1:npti,:) < 0._wp ) pv_ip(1:npti,:) = 0._wp ! v_ip must be >= 0 -
NEMO/branches/2020/tickets_icb_1900/src/ICE/icewri.F90
r13899 r14016 160 160 IF( iom_use('icebrv_cat' ) ) CALL iom_put( 'icebrv_cat' , bv_i * 100. * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! brine volume 161 161 IF( iom_use('iceapnd_cat' ) ) CALL iom_put( 'iceapnd_cat' , a_ip * zmsk00l ) ! melt pond frac for categories 162 IF( iom_use('icevpnd_cat' ) ) CALL iom_put( 'icevpnd_cat' , v_ip * zmsk00l ) ! melt pond volume for categories 162 163 IF( iom_use('icehpnd_cat' ) ) CALL iom_put( 'icehpnd_cat' , h_ip * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond thickness for categories 163 164 IF( iom_use('icehlid_cat' ) ) CALL iom_put( 'icehlid_cat' , h_il * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! melt pond lid thickness for categories 164 IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac for categories165 IF( iom_use('iceafpnd_cat') ) CALL iom_put( 'iceafpnd_cat', a_ip_frac * zmsk00l ) ! melt pond frac per ice area for categories 165 166 IF( iom_use('iceaepnd_cat') ) CALL iom_put( 'iceaepnd_cat', a_ip_eff * zmsk00l ) ! melt pond effective frac for categories 166 167 IF( iom_use('icealb_cat' ) ) CALL iom_put( 'icealb_cat' , alb_ice * zmsk00l + zmiss_val * ( 1._wp - zmsk00l ) ) ! ice albedo for categories -
NEMO/branches/2020/tickets_icb_1900/src/OCE/C1D/step_c1d.F90
r14012 r14016 104 104 IF( ln_tradmp ) CALL tra_dmp( kstp, Nbb, Nnn, ts, Nrhs ) ! internal damping trends- tracers 105 105 IF(.NOT.ln_linssh)CALL tra_adv( kstp, Nbb, Nnn, ts, Nrhs ) ! horizontal & vertical advection 106 IF( ln_zdfmfc ) CALL tra_mfc( kstp, Nbb , ts, Nrhs ) ! Mass Flux Convection 106 107 IF( ln_zdfosm ) CALL tra_osm( kstp, Nnn , ts, Nrhs ) ! OSMOSIS non-local tracer fluxes 107 108 CALL tra_zdf( kstp, Nbb, Nnn, Nrhs, ts, Naa ) ! vertical mixing … … 122 123 CALL dyn_atf ( kstp, Nbb, Nnn, Naa , uu, vv, e3t, e3u, e3v ) ! time filtering of "now" fields 123 124 IF(.NOT.ln_linssh)CALL ssh_atf ( kstp, Nbb, Nnn, Naa , ssh ) ! time filtering of "now" sea surface height 124 IF( kstp == nit000 .AND. ln_linssh) THEN 125 ssh(:,:,Naa) = ssh(:,:,Nnn) ! init ssh after in ln_linssh case 125 IF( kstp == nit000 .AND. ln_linssh) THEN 126 ssh(:,:,Naa) = ssh(:,:,Nnn) ! init ssh after in ln_linssh case 126 127 ENDIF 127 128 ! -
NEMO/branches/2020/tickets_icb_1900/src/OCE/DYN/dynspg.F90
r13899 r14016 6 6 !! History : 1.0 ! 2005-12 (C. Talandier, G. Madec, V. Garnier) Original code 7 7 !! 3.2 ! 2009-07 (R. Benshila) Suppression of rigid-lid option 8 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) add Bernoulli Head for 9 !! wave coupling 8 10 !!---------------------------------------------------------------------- 9 11 … … 19 21 USE sbc_ice , ONLY : snwice_mass, snwice_mass_b 20 22 USE sbcapr ! surface boundary condition: atmospheric pressure 23 USE sbcwave, ONLY : bhd_wave 21 24 USE dynspg_exp ! surface pressure gradient (dyn_spg_exp routine) 22 25 USE dynspg_ts ! surface pressure gradient (dyn_spg_ts routine) … … 143 146 ENDIF 144 147 ! 148 IF( ln_wave .and. ln_bern_srfc ) THEN !== Add J terms: depth-independent Bernoulli head 149 DO_2D( 0, 0, 0, 0 ) 150 spgu(ji,jj) = spgu(ji,jj) + ( bhd_wave(ji+1,jj) - bhd_wave(ji,jj) ) / e1u(ji,jj) !++ bhd_wave from wave model in m2/s2 [BHD parameters in WW3] 151 spgv(ji,jj) = spgv(ji,jj) + ( bhd_wave(ji,jj+1) - bhd_wave(ji,jj) ) / e2v(ji,jj) 152 END_2D 153 ENDIF 154 ! 145 155 DO_3D( 0, 0, 0, 0, 1, jpkm1 ) !== Add all terms to the general trend 146 156 puu(ji,jj,jk,Krhs) = puu(ji,jj,jk,Krhs) + spgu(ji,jj) -
NEMO/branches/2020/tickets_icb_1900/src/OCE/DYN/dynvor.F90
r13899 r14016 21 21 !! - ! 2018-03 (G. Madec) add two new schemes (ln_dynvor_enT and ln_dynvor_eet) 22 22 !! - ! 2018-04 (G. Madec) add pre-computed gradient for metric term calculation 23 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) add vortex force trends (ln_vortex_force=T) 23 24 !!---------------------------------------------------------------------- 24 25 … … 37 38 USE trddyn ! trend manager: dynamics 38 39 USE sbcwave ! Surface Waves (add Stokes-Coriolis force) 39 USE sbc_oce , ONLY : ln_stcor! use Stoke-Coriolis force40 USE sbc_oce, ONLY : ln_stcor, ln_vortex_force ! use Stoke-Coriolis force 40 41 ! 41 42 USE lbclnk ! ocean lateral boundary conditions (or mpp link) … … 121 122 ALLOCATE( ztrdu(jpi,jpj,jpk), ztrdv(jpi,jpj,jpk) ) 122 123 ! 123 ztrdu(:,:,:) = puu(:,:,:,Krhs) !* planetary vorticity trend (including Stokes-Coriolis force)124 ztrdu(:,:,:) = puu(:,:,:,Krhs) !* planetary vorticity trend 124 125 ztrdv(:,:,:) = pvv(:,:,:,Krhs) 125 126 SELECT CASE( nvor_scheme ) 126 127 CASE( np_ENS ) ; CALL vor_ens( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! enstrophy conserving scheme 127 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend128 128 CASE( np_ENE, np_MIX ) ; CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme 129 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend130 129 CASE( np_ENT ) ; CALL vor_enT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (T-pts) 131 IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend132 130 CASE( np_EET ) ; CALL vor_eeT( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy conserving scheme (een with e3t) 133 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend134 131 CASE( np_EEN ) ; CALL vor_een( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! energy & enstrophy scheme 135 IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend136 132 END SELECT 137 133 ztrdu(:,:,:) = puu(:,:,:,Krhs) - ztrdu(:,:,:) … … 161 157 CASE( np_ENT ) !* energy conserving scheme (T-pts) 162 158 CALL vor_enT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 163 IF( ln_stcor ) CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 159 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 160 CALL vor_enT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 161 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 162 CALL vor_enT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 163 ENDIF 164 164 CASE( np_EET ) !* energy conserving scheme (een scheme using e3t) 165 165 CALL vor_eeT( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 166 IF( ln_stcor ) CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 166 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 167 CALL vor_eeT( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 168 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 169 CALL vor_eeT( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 170 ENDIF 167 171 CASE( np_ENE ) !* energy conserving scheme 168 172 CALL vor_ene( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 169 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 173 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 174 CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 175 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 176 CALL vor_ene( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 177 ENDIF 170 178 CASE( np_ENS ) !* enstrophy conserving scheme 171 179 CALL vor_ens( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 172 IF( ln_stcor ) CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 180 181 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 182 CALL vor_ens( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 183 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 184 CALL vor_ens( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 185 ENDIF 173 186 CASE( np_MIX ) !* mixed ene-ens scheme 174 187 CALL vor_ens( kt, Kmm, nrvm, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! relative vorticity or metric trend (ens) 175 188 CALL vor_ene( kt, Kmm, ncor, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! planetary vorticity trend (ene) 176 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 189 IF( ln_stcor ) CALL vor_ene( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 190 IF( ln_vortex_force ) CALL vor_ens( kt, Kmm, nrvm, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add vortex force 177 191 CASE( np_EEN ) !* energy and enstrophy conserving scheme 178 192 CALL vor_een( kt, Kmm, ntot, puu(:,:,:,Kmm) , pvv(:,:,:,Kmm) , puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! total vorticity trend 179 IF( ln_stcor ) CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 193 IF( ln_stcor .AND. .NOT. ln_vortex_force ) THEN 194 CALL vor_een( kt, Kmm, ncor, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend 195 ELSE IF( ln_stcor .AND. ln_vortex_force ) THEN 196 CALL vor_een( kt, Kmm, ntot, usd, vsd, puu(:,:,:,Krhs), pvv(:,:,:,Krhs) ) ! add the Stokes-Coriolis trend and vortex force 197 ENDIF 180 198 END SELECT 181 199 ! -
NEMO/branches/2020/tickets_icb_1900/src/OCE/DYN/dynzad.F90
r13899 r14016 16 16 USE trd_oce ! trends: ocean variables 17 17 USE trddyn ! trend manager: dynamics 18 USE sbcwave, ONLY: wsd ! Surface Waves (add vertical Stokes-drift) 18 19 ! 19 20 USE in_out_manager ! I/O manager … … 79 80 DO jk = 2, jpkm1 ! Vertical momentum advection at level w and u- and v- vertical 80 81 DO_2D( 0, 1, 0, 1 ) ! vertical fluxes 82 IF( ln_vortex_force ) THEN 83 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wsd(ji,jj,jk) ) 84 ELSE 81 85 zww(ji,jj) = 0.25_wp * e1e2t(ji,jj) * ww(ji,jj,jk) 86 ENDIF 82 87 END_2D 83 88 DO_2D( 0, 0, 0, 0 ) ! vertical momentum advection at w-point -
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/cpl_oasis3.F90
r13899 r14016 66 66 INTEGER :: nsnd ! total number of fields sent 67 67 INTEGER :: ncplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data 68 INTEGER, PUBLIC, PARAMETER :: nmaxfld=6 0! Maximum number of coupling fields68 INTEGER, PUBLIC, PARAMETER :: nmaxfld=62 ! Maximum number of coupling fields 69 69 INTEGER, PUBLIC, PARAMETER :: nmaxcat=5 ! Maximum number of coupling fields 70 70 INTEGER, PUBLIC, PARAMETER :: nmaxcpl=5 ! Maximum number of coupling fields -
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbc_oce.F90
r13899 r14016 12 12 !! 4.0 ! 2016-06 (L. Brodeau) new unified bulk routine (based on AeroBulk) 13 13 !! 4.0 ! 2019-03 (F. Lemarié, G. Samson) add compatibility with ABL mode 14 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) modified wave parameters in namelist 14 15 !!---------------------------------------------------------------------- 15 16 … … 36 37 LOGICAL , PUBLIC :: ln_blk !: bulk formulation 37 38 LOGICAL , PUBLIC :: ln_abl !: Atmospheric boundary layer model 39 LOGICAL , PUBLIC :: ln_wave !: wave in the system (forced or coupled) 38 40 #if defined key_oasis3 39 41 LOGICAL , PUBLIC :: lk_oasis = .TRUE. !: OASIS used … … 56 58 ! !: = 1 global mean of e-p-r set to zero at each nn_fsbc time step 57 59 ! !: = 2 annual global mean of e-p-r set to zero 58 LOGICAL , PUBLIC :: ln_wave !: true if some coupling with wave model59 LOGICAL , PUBLIC :: ln_cdgw !: true if neutral drag coefficient from wave model60 LOGICAL , PUBLIC :: ln_sdw !: true if 3d stokes drift from wave model61 LOGICAL , PUBLIC :: ln_tauwoc !: true if normalized stress from wave is used62 LOGICAL , PUBLIC :: ln_tauw !: true if ocean stress components from wave is used63 LOGICAL , PUBLIC :: ln_stcor !: true if Stokes-Coriolis term is used64 !65 INTEGER , PUBLIC :: nn_sdrift ! type of parameterization to calculate vertical Stokes drift66 !67 60 LOGICAL , PUBLIC :: ln_icebergs !: Icebergs 68 61 ! … … 71 64 ! !!* namsbc_cpl namelist * 72 65 INTEGER , PUBLIC :: nn_cats_cpl !: Number of sea ice categories over which the coupling is carried out 73 66 ! 67 ! !!* namsbc_wave namelist * 68 LOGICAL , PUBLIC :: ln_sdw !: =T 3d stokes drift from wave model 69 LOGICAL , PUBLIC :: ln_stcor !: =T if Stokes-Coriolis and tracer advection terms are used 70 LOGICAL , PUBLIC :: ln_cdgw !: =T neutral drag coefficient from wave model 71 LOGICAL , PUBLIC :: ln_tauoc !: =T if normalized stress from wave is used 72 LOGICAL , PUBLIC :: ln_wave_test !: =T wave test case (constant Stokes drift) 73 LOGICAL , PUBLIC :: ln_charn !: =T Chranock coefficient from wave model 74 LOGICAL , PUBLIC :: ln_taw !: =T wind stress corrected by wave intake 75 LOGICAL , PUBLIC :: ln_phioc !: =T TKE surface BC from wave model 76 LOGICAL , PUBLIC :: ln_bern_srfc !: Bernoulli head, waves' inuced pressure 77 LOGICAL , PUBLIC :: ln_breivikFV_2016 !: Breivik 2016 profile 78 LOGICAL , PUBLIC :: ln_vortex_force !: vortex force activation 79 LOGICAL , PUBLIC :: ln_stshear !: Stoked Drift shear contribution in zdftke 80 ! 74 81 !!---------------------------------------------------------------------- 75 82 !! switch definition (improve readability) … … 81 88 INTEGER , PUBLIC, PARAMETER :: jp_purecpl = 5 !: Pure ocean-atmosphere Coupled formulation 82 89 INTEGER , PUBLIC, PARAMETER :: jp_none = 6 !: for OPA when doing coupling via SAS module 83 84 !!---------------------------------------------------------------------- 85 !! Stokes drift parametrization definition 86 !!---------------------------------------------------------------------- 87 INTEGER , PUBLIC, PARAMETER :: jp_breivik_2014 = 0 !: Breivik 2014: v_z=v_0*[exp(2*k*z)/(1-8*k*z)] 88 INTEGER , PUBLIC, PARAMETER :: jp_li_2017 = 1 !: Li et al 2017: Stokes drift based on Phillips spectrum (Breivik 2016) 89 ! with depth averaged profile 90 INTEGER , PUBLIC, PARAMETER :: jp_peakfr = 2 !: Li et al 2017: using the peak wave number read from wave model instead 91 ! of the inverse depth scale 92 LOGICAL , PUBLIC :: ll_st_bv2014 = .FALSE. ! logical indicator, .true. if Breivik 2014 parameterisation is active. 93 LOGICAL , PUBLIC :: ll_st_li2017 = .FALSE. ! logical indicator, .true. if Li 2017 parameterisation is active. 94 LOGICAL , PUBLIC :: ll_st_bv_li = .FALSE. ! logical indicator, .true. if either Breivik or Li parameterisation is active. 95 LOGICAL , PUBLIC :: ll_st_peakfr = .FALSE. ! logical indicator, .true. if using Li 2017 with peak wave number 96 90 ! 97 91 !!---------------------------------------------------------------------- 98 92 !! component definition -
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcblk.F90
r13899 r14016 314 314 ENDIF 315 315 END DO 316 !317 IF( ln_wave ) THEN318 !Activated wave module but neither drag nor stokes drift activated319 IF( .NOT.(ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) ) THEN320 CALL ctl_stop( 'STOP', 'Ask for wave coupling but ln_cdgw=F, ln_sdw=F, ln_tauwoc=F, ln_stcor=F' )321 !drag coefficient read from wave model definable only with mfs bulk formulae and core322 ELSEIF(ln_cdgw .AND. .NOT. ln_NCAR ) THEN323 CALL ctl_stop( 'drag coefficient read from wave model definable only with NCAR and CORE bulk formulae')324 ELSEIF(ln_stcor .AND. .NOT. ln_sdw) THEN325 CALL ctl_stop( 'Stokes-Coriolis term calculated only if activated Stokes Drift ln_sdw=T')326 ENDIF327 ELSE328 IF( ln_cdgw .OR. ln_sdw .OR. ln_tauwoc .OR. ln_stcor ) &329 & CALL ctl_stop( 'Not Activated Wave Module (ln_wave=F) but asked coupling ', &330 & 'with drag coefficient (ln_cdgw =T) ' , &331 & 'or Stokes Drift (ln_sdw=T) ' , &332 & 'or ocean stress modification due to waves (ln_tauwoc=T) ', &333 & 'or Stokes-Coriolis term (ln_stcori=T)' )334 ENDIF335 316 ! 336 317 IF( ln_abl ) THEN ! ABL: read 3D fields for wind, temperature, humidity and pressure gradient -
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbcblk_algo_ecmwf.F90
r13899 r14016 17 17 !!---------------------------------------------------------------------- 18 18 !! History : 4.0 ! 2016-02 (L.Brodeau) Original code 19 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) Charnock coeff from wave model 19 20 !!---------------------------------------------------------------------- 20 21 … … 31 32 USE in_out_manager ! I/O manager 32 33 USE prtctl ! Print control 33 USE sbcwave, ONLY : cdn_wave! wave module34 USE sbcwave, ONLY : charn ! wave module 34 35 #if defined key_si3 || defined key_cice 35 36 USE sbc_ice ! Surface boundary condition: ice fields … … 233 234 u_star = 0.035_wp*U_blk*ztmp1/ztmp0 ! (u* = 0.035*Un10) 234 235 235 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 236 IF (ln_charn) THEN ! Charnock value if wave coupling 237 z0 = charn*u_star*u_star/grav + 0.11_wp*znu_a/u_star 238 ELSE 239 z0 = charn0*u_star*u_star/grav + 0.11_wp*znu_a/u_star 240 ENDIF 241 236 242 z0 = MIN( MAX(ABS(z0), 1.E-9) , 1._wp ) ! (prevents FPE from stupid values from masked region later on) 237 243 … … 302 308 ztmp2 = u_star*u_star 303 309 ztmp1 = znu_a/u_star 304 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 310 IF (ln_charn) THEN ! Charnock value if wave coupling 311 z0 = MIN( ABS( alpha_M*ztmp1 + charn*ztmp2/grav ) , 0.001_wp) 312 ELSE 313 z0 = MIN( ABS( alpha_M*ztmp1 + charn0*ztmp2/grav ) , 0.001_wp) 314 ENDIF 305 315 z0t = MIN( ABS( alpha_H*ztmp1 ) , 0.001_wp) ! eq.3.26, Chap.3, p.34, IFS doc - Cy31r1 306 316 z0q = MIN( ABS( alpha_Q*ztmp1 ) , 0.001_wp) -
NEMO/branches/2020/tickets_icb_1900/src/OCE/SBC/sbccpl.F90
r13899 r14016 8 8 !! 3.1 ! 2009_02 (G. Madec, S. Masson, E. Maisonave, A. Caubel) generic coupled interface 9 9 !! 3.4 ! 2011_11 (C. Harris) more flexibility + multi-category fields 10 !! 4.2 ! 2020-12 (G. Madec, E. Clementi) wave coupling updates 10 11 !!---------------------------------------------------------------------- 11 12 … … 106 107 INTEGER, PARAMETER :: jpr_fraqsr = 42 ! fraction of solar net radiation absorbed in the first ocean level 107 108 INTEGER, PARAMETER :: jpr_mslp = 43 ! mean sea level pressure 108 INTEGER, PARAMETER :: jpr_hsig = 44 ! Hsig 109 INTEGER, PARAMETER :: jpr_phioc = 45 ! Wave=>ocean energy flux 110 INTEGER, PARAMETER :: jpr_sdrftx = 46 ! Stokes drift on grid 1 111 INTEGER, PARAMETER :: jpr_sdrfty = 47 ! Stokes drift on grid 2 109 !** surface wave coupling ** 110 INTEGER, PARAMETER :: jpr_hsig = 44 ! Hsig 111 INTEGER, PARAMETER :: jpr_phioc = 45 ! Wave=>ocean energy flux 112 INTEGER, PARAMETER :: jpr_sdrftx = 46 ! Stokes drift on grid 1 113 INTEGER, PARAMETER :: jpr_sdrfty = 47 ! Stokes drift on grid 2 112 114 INTEGER, PARAMETER :: jpr_wper = 48 ! Mean wave period 113 115 INTEGER, PARAMETER :: jpr_wnum = 49 ! Mean wavenumber 114 INTEGER, PARAMETER :: jpr_ tauwoc= 50 ! Stress fraction adsorbed by waves116 INTEGER, PARAMETER :: jpr_wstrf = 50 ! Stress fraction adsorbed by waves 115 117 INTEGER, PARAMETER :: jpr_wdrag = 51 ! Neutral surface drag coefficient 116 INTEGER, PARAMETER :: jpr_isf = 52 117 INTEGER, PARAMETER :: jpr_icb = 53 118 INTEGER, PARAMETER :: jpr_wfreq = 54 ! Wave peak frequency 119 INTEGER, PARAMETER :: jpr_tauwx = 55 ! x component of the ocean stress from waves 120 INTEGER, PARAMETER :: jpr_tauwy = 56 ! y component of the ocean stress from waves 121 INTEGER, PARAMETER :: jpr_ts_ice = 57 ! Sea ice surface temp 122 123 INTEGER, PARAMETER :: jprcv = 57 ! total number of fields received 118 INTEGER, PARAMETER :: jpr_charn = 52 ! Chranock coefficient 119 INTEGER, PARAMETER :: jpr_twox = 53 ! wave to ocean momentum flux 120 INTEGER, PARAMETER :: jpr_twoy = 54 ! wave to ocean momentum flux 121 INTEGER, PARAMETER :: jpr_tawx = 55 ! net wave-supported stress 122 INTEGER, PARAMETER :: jpr_tawy = 56 ! net wave-supported stress 123 INTEGER, PARAMETER :: jpr_bhd = 57 ! Bernoulli head. waves' induced surface pressure 124 INTEGER, PARAMETER :: jpr_tusd = 58 ! zonal stokes transport 125 INTEGER, PARAMETER :: jpr_tvsd = 59 ! meridional stokes tranmport 126 INTEGER, PARAMETER :: jpr_isf = 60 127 INTEGER, PARAMETER :: jpr_icb = 61 128 INTEGER, PARAMETER :: jpr_ts_ice = 62 ! Sea ice surface temp 129 130 INTEGER, PARAMETER :: jprcv = 62 ! total number of fields received 124 131 125 132 INTEGER, PARAMETER :: jps_fice = 1 ! ice fraction sent to the atmosphere … … 184 191 & sn_snd_thick1, sn_snd_cond, sn_snd_mpnd , sn_snd_sstfrz, sn_snd_ttilyr 185 192 ! ! Received from the atmosphere 186 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_ tauw, sn_rcv_dqnsdt, sn_rcv_qsr, &193 TYPE(FLD_C) :: sn_rcv_w10m, sn_rcv_taumod, sn_rcv_tau, sn_rcv_dqnsdt, sn_rcv_qsr, & 187 194 & sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf, sn_rcv_ts_ice 188 195 TYPE(FLD_C) :: sn_rcv_cal, sn_rcv_iceflx, sn_rcv_co2, sn_rcv_mslp, sn_rcv_icb, sn_rcv_isf 189 ! Send to waves196 ! ! Send to waves 190 197 TYPE(FLD_C) :: sn_snd_ifrac, sn_snd_crtw, sn_snd_wlev 191 ! Received from waves192 TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, sn_rcv_tauwoc,&193 sn_rcv_wdrag, sn_rcv_wfreq198 ! ! Received from waves 199 TYPE(FLD_C) :: sn_rcv_hsig, sn_rcv_phioc, sn_rcv_sdrfx, sn_rcv_sdrfy, sn_rcv_wper, sn_rcv_wnum, & 200 & sn_rcv_wstrf, sn_rcv_wdrag, sn_rcv_charn, sn_rcv_taw, sn_rcv_bhd, sn_rcv_tusd, sn_rcv_tvsd 194 201 ! ! Other namelist parameters 195 202 INTEGER :: nn_cplmodel ! Maximum number of models to/from which NEMO is potentialy sending/receiving data … … 274 281 & sn_snd_ifrac , sn_snd_crtw , sn_snd_wlev , sn_rcv_hsig , sn_rcv_phioc , & 275 282 & sn_rcv_w10m , sn_rcv_taumod, sn_rcv_tau , sn_rcv_dqnsdt, sn_rcv_qsr , & 276 & sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_ tauwoc, &277 & sn_rcv_ wdrag , sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal ,&278 & sn_rcv_ iceflx, sn_rcv_co2 , sn_rcv_mslp ,&279 & sn_rcv_ic b , sn_rcv_isf , sn_rcv_wfreq, sn_rcv_tauw , &280 & sn_rcv_ts_ice 283 & sn_rcv_sdrfx , sn_rcv_sdrfy , sn_rcv_wper , sn_rcv_wnum , sn_rcv_wstrf , & 284 & sn_rcv_charn , sn_rcv_taw , sn_rcv_bhd , sn_rcv_tusd , sn_rcv_tvsd, & 285 & sn_rcv_wdrag , sn_rcv_qns , sn_rcv_emp , sn_rcv_rnf , sn_rcv_cal , & 286 & sn_rcv_iceflx, sn_rcv_co2 , sn_rcv_icb , sn_rcv_isf , sn_rcv_ts_ice 287 281 288 !!--------------------------------------------------------------------- 282 289 ! … … 319 326 WRITE(numout,*)' sea ice heat fluxes = ', TRIM(sn_rcv_iceflx%cldes), ' (', TRIM(sn_rcv_iceflx%clcat), ')' 320 327 WRITE(numout,*)' atm co2 = ', TRIM(sn_rcv_co2%cldes ), ' (', TRIM(sn_rcv_co2%clcat ), ')' 328 WRITE(numout,*)' Sea ice surface skin temperature= ', TRIM(sn_rcv_ts_ice%cldes), ' (', TRIM(sn_rcv_ts_ice%clcat), ')' 329 WRITE(numout,*)' surface waves:' 321 330 WRITE(numout,*)' significant wave heigth = ', TRIM(sn_rcv_hsig%cldes ), ' (', TRIM(sn_rcv_hsig%clcat ), ')' 322 331 WRITE(numout,*)' wave to oce energy flux = ', TRIM(sn_rcv_phioc%cldes ), ' (', TRIM(sn_rcv_phioc%clcat ), ')' … … 325 334 WRITE(numout,*)' Mean wave period = ', TRIM(sn_rcv_wper%cldes ), ' (', TRIM(sn_rcv_wper%clcat ), ')' 326 335 WRITE(numout,*)' Mean wave number