Changeset 15388
- Timestamp:
- 2021-10-17T13:33:47+02:00 (3 years ago)
- Location:
- NEMO/trunk/src/ICE
- Files:
-
- 5 edited
-
ice.F90 (modified) (1 diff)
-
icesbc.F90 (modified) (6 diffs)
-
icethd.F90 (modified) (4 diffs)
-
icethd_do.F90 (modified) (3 diffs)
-
icewri.F90 (modified) (4 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/trunk/src/ICE/ice.F90
r15334 r15388 196 196 ! ! = 0 Grenfell and Maykut 1977 (depends on cloudiness and is 0 when there is snow) 197 197 ! ! = 1 Lebrun 2019 (equals 0.3 anytime with different melting/dry snw conductivities) 198 199 ! !!** namelist (namthd) ** 200 LOGICAL , PUBLIC :: ln_icedH ! activate ice thickness change from growing/melting (T) or not (F) 201 LOGICAL , PUBLIC :: ln_icedA ! activate lateral melting param. (T) or not (F) 202 LOGICAL , PUBLIC :: ln_icedO ! activate ice growth in open-water (T) or not (F) 203 LOGICAL , PUBLIC :: ln_icedS ! activate gravity drainage and flushing (T) or not (F) 204 LOGICAL , PUBLIC :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean 205 ! 206 ! !!** namelist (namthd_do) ** 207 REAL(wp), PUBLIC :: rn_hinew ! thickness for new ice formation (m) 208 LOGICAL , PUBLIC :: ln_frazil ! use of frazil ice collection as function of wind (T) or not (F) 209 REAL(wp), PUBLIC :: rn_maxfraz ! maximum portion of frazil ice collecting at the ice bottom 210 REAL(wp), PUBLIC :: rn_vfraz ! threshold drift speed for collection of bottom frazil ice 211 REAL(wp), PUBLIC :: rn_Cfraz ! squeezing coefficient for collection of bottom frazil ice 198 212 ! 199 213 ! !!** ice-vertical diffusion namelist (namthd_zdf) ** -
NEMO/trunk/src/ICE/icesbc.F90
r14595 r15388 109 109 !! dqns_ice = non solar heat sensistivity [W/m2] 110 110 !! qemp_oce, qemp_ice, qprec_ice, qevap_ice = sensible heat (associated with evap & precip) [W/m2] 111 !! + these fields 112 !! qsb_ice_bot = sensible heat at the ice bottom [W/m2] 113 !! fhld, qlead = heat budget in the leads [W/m2] 111 114 !! + some fields that are not used outside this module: 112 115 !! qla_ice = latent heat flux over ice [W/m2] … … 117 120 INTEGER, INTENT(in) :: kt ! ocean time step 118 121 INTEGER, INTENT(in) :: ksbc ! flux formulation (user defined, bulk or Pure Coupled) 119 !120 INTEGER :: ji, jj, jl ! dummy loop index121 REAL(wp) :: zmiss_val ! missing value retrieved from xios122 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zalb, zmsk00 ! 2D workspace123 122 !!-------------------------------------------------------------------- 124 123 ! … … 130 129 WRITE(numout,*)'~~~~~~~~~~~~~~~' 131 130 ENDIF 132 133 ! get missing value from xml 134 CALL iom_miss_val( "icetemp", zmiss_val ) 135 136 ! --- ice albedo --- ! 131 ! !== ice albedo ==! 137 132 CALL ice_alb( t_su, h_i, h_s, ln_pnd_alb, a_ip_eff, h_ip, cloud_fra, alb_ice ) 138 139 133 ! 140 134 SELECT CASE( ksbc ) !== fluxes over sea ice ==! … … 142 136 CASE( jp_usr ) !--- user defined formulation 143 137 CALL usrdef_sbc_ice_flx( kt, h_s, h_i ) 144 CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation138 CASE( jp_blk, jp_abl ) !--- bulk formulation & ABL formulation 145 139 CALL blk_ice_2 ( t_su, h_s, h_i, alb_ice, & 146 140 & theta_air_zt(:,:), q_air_zt(:,:), & ! #LB: known from "sbc_oce" module... … … 156 150 IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) 157 151 END SELECT 158 159 !--- output ice albedo and surface albedo ---! 160 IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN 161 162 ALLOCATE( zalb(jpi,jpj), zmsk00(jpi,jpj) ) 163 164 WHERE( at_i_b < 1.e-03 ) 165 zmsk00(:,:) = 0._wp 166 zalb (:,:) = rn_alb_oce 167 ELSEWHERE 168 zmsk00(:,:) = 1._wp 169 zalb (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 170 END WHERE 171 ! ice albedo 172 CALL iom_put( 'icealb' , zalb * zmsk00 + zmiss_val * ( 1._wp - zmsk00 ) ) 173 ! ice+ocean albedo 174 zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 175 CALL iom_put( 'albedo' , zalb ) 176 177 DEALLOCATE( zalb, zmsk00 ) 178 179 ENDIF 152 ! !== some fluxes at the ice-ocean interface and in the leads 153 CALL ice_flx_other 180 154 ! 181 155 IF( ln_timing ) CALL timing_stop('icesbc') … … 270 244 271 245 246 SUBROUTINE ice_flx_other 247 !!----------------------------------------------------------------------- 248 !! *** ROUTINE ice_flx_other *** 249 !! 250 !! ** Purpose : prepare necessary fields for thermo calculations 251 !! 252 !! ** Inputs : u_ice, v_ice, ssu_m, ssv_m, utau, vtau 253 !! frq_m, qsr_oce, qns_oce, qemp_oce, e3t_m, sst_m 254 !! ** Outputs : qsb_ice_bot, fhld, qlead 255 !!----------------------------------------------------------------------- 256 INTEGER :: ji, jj ! dummy loop indices 257 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1 258 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04) 259 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient 260 REAL(wp), DIMENSION(jpi,jpj) :: zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2) 261 !!----------------------------------------------------------------------- 262 ! 263 ! computation of friction velocity at T points 264 IF( ln_icedyn ) THEN 265 DO_2D( 0, 0, 0, 0 ) 266 zu_io = u_ice(ji ,jj ) - ssu_m(ji ,jj ) 267 zu_iom1 = u_ice(ji-1,jj ) - ssu_m(ji-1,jj ) 268 zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj ) 269 zv_iom1 = v_ice(ji ,jj-1) - ssv_m(ji ,jj-1) 270 ! 271 zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * tmask(ji,jj,1) 272 zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) + & 273 & ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) ) 274 END_2D 275 ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean 276 DO_2D( 0, 0, 0, 0 ) 277 zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * & 278 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) & 279 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1) 280 zvel(ji,jj) = 0._wp 281 END_2D 282 ENDIF 283 CALL lbc_lnk( 'icesbc', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp ) 284 ! 285 !--------------------------------------------------------------------! 286 ! Partial computation of forcing for the thermodynamic sea ice model 287 !--------------------------------------------------------------------! 288 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! needed for qlead 289 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice 290 ! 291 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- ! 292 zqld = tmask(ji,jj,1) * rDt_ice * & 293 & ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + & 294 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) ) 295 296 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- ! 297 ! (mostly<0 but >0 if supercooling) 298 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst) 299 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0 300 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0 301 302 ! --- Sensible ocean-to-ice heat flux (W/m2) --- ! 303 ! (mostly>0 but <0 if supercooling) 304 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin ) 305 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) ) 306 307 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach 308 ! the freezing point, so that we do not have SST < T_freeze 309 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg 310 ! The following formulation is ok for both normal conditions and supercooling 311 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) ) 312 313 ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously 314 ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary) 315 IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN 316 zqfr = 0._wp 317 zqfr_pos = 0._wp 318 qsb_ice_bot(ji,jj) = 0._wp 319 ENDIF 320 ! 321 ! --- Energy Budget of the leads (qlead, J.m-2) --- ! 322 ! qlead is the energy received from the atm. in the leads. 323 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2) 324 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2) 325 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN 326 ! upper bound for fhld: fhld should be equal to zqld 327 ! but we have to make sure that this heat will not make the sst drop below the freezing point 328 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos 329 ! The following formulation is ok for both normal conditions and supercooling 330 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90 331 & - qsb_ice_bot(ji,jj) ) 332 qlead(ji,jj) = 0._wp 333 ELSE 334 fhld (ji,jj) = 0._wp 335 ! upper bound for qlead: qlead should be equal to zqld 336 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point. 337 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb) 338 ! and freezing point is reached if zqfr = zqld - qsb*a/dt 339 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr 340 ! The following formulation is ok for both normal conditions and supercooling 341 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr ) 342 ENDIF 343 ! 344 ! If ice is landfast and ice concentration reaches its max 345 ! => stop ice formation in open water 346 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp 347 ! 348 ! If the grid cell is almost fully covered by ice (no leads) 349 ! => stop ice formation in open water 350 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp 351 ! 352 ! If ln_leadhfx is false 353 ! => do not use energy of the leads to melt sea-ice 354 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp 355 ! 356 END_2D 357 358 ! In case we bypass open-water ice formation 359 IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp 360 ! In case we bypass growing/melting from top and bottom 361 IF( .NOT. ln_icedH ) THEN 362 qsb_ice_bot(:,:) = 0._wp 363 fhld (:,:) = 0._wp 364 ENDIF 365 366 END SUBROUTINE ice_flx_other 367 368 272 369 SUBROUTINE ice_sbc_init 273 370 !!------------------------------------------------------------------- -
NEMO/trunk/src/ICE/icethd.F90
r15385 r15388 20 20 USE sbc_oce , ONLY : sss_m, sst_m, e3t_m, utau, vtau, ssu_m, ssv_m, frq_m, sprecip, ln_cpl 21 21 USE sbc_ice , ONLY : qsr_oce, qns_oce, qemp_oce, qsr_ice, qns_ice, dqns_ice, evap_ice, qprec_ice, qevap_ice, & 22 & qml_ice, qcn_ice, qtr_ice_top , utau_ice, vtau_ice22 & qml_ice, qcn_ice, qtr_ice_top 23 23 USE ice1D ! sea-ice: thermodynamics variables 24 24 USE icethd_zdf ! sea-ice: vertical heat diffusion … … 48 48 PUBLIC ice_thd_init ! called by ice_init 49 49 50 !!** namelist (namthd) **51 LOGICAL :: ln_icedH ! activate ice thickness change from growing/melting (T) or not (F)52 LOGICAL :: ln_icedA ! activate lateral melting param. (T) or not (F)53 LOGICAL :: ln_icedO ! activate ice growth in open-water (T) or not (F)54 LOGICAL :: ln_icedS ! activate gravity drainage and flushing (T) or not (F)55 LOGICAL :: ln_leadhfx ! heat in the leads is used to melt sea-ice before warming the ocean56 57 50 !! for convergence tests 58 51 REAL(wp), ALLOCATABLE, DIMENSION(:,:,:) :: ztice_cvgerr, ztice_cvgstp … … 109 102 ztice_cvgerr = 0._wp ; ztice_cvgstp = 0._wp 110 103 ENDIF 111 112 ! --- calculate (some) heat fluxes and frazil ice --- ! 113 ! out => sensible heat (qsb_ice_bot) & heat budget in the leads (fhld & qlead) 114 ! out => ice collection thickness (ht_i_new) and fraction of frazil (fraz_frac) 115 CALL thd_prep 116 104 ! 105 CALL ice_thd_frazil !--- frazil ice: collection thickness (ht_i_new) & fraction of frazil (fraz_frac) 106 ! 117 107 !-------------------------------------------------------------------------------------------! 118 108 ! Thermodynamic computation (only on grid points covered by ice) => loop over ice categories … … 260 250 ! 261 251 END SUBROUTINE ice_thd_mono 262 263 SUBROUTINE thd_prep264 !!-----------------------------------------------------------------------265 !! *** ROUTINE thd_prep ***266 !!267 !! ** Purpose : prepare necessary fields for thermo calculations268 !!269 !! For the fluxes270 !! ** Inputs : u_ice, v_ice, ssu_m, ssv_m, utau, vtau271 !! frq_m, qsr_oce, qns_oce, qemp_oce, e3t_m, sst_m272 !! ** Outputs : qsb_ice_bot, fhld, qlead273 !!274 !! For the collection thickness (frazil)275 !! ** Inputs : u_ice, v_ice, utau_ice, vtau_ice276 !! ** Ouputs : ht_i_new, fraz_frac277 !!-----------------------------------------------------------------------278 INTEGER :: ji, jj ! dummy loop indices279 REAL(wp) :: zfric_u, zqld, zqfr, zqfr_neg, zqfr_pos, zu_io, zv_io, zu_iom1, zv_iom1280 REAL(wp), PARAMETER :: zfric_umin = 0._wp ! lower bound for the friction velocity (cice value=5.e-04)281 REAL(wp), PARAMETER :: zch = 0.0057_wp ! heat transfer coefficient282 REAL(wp), DIMENSION(jpi,jpj) :: zfric, zvel ! ice-ocean velocity (m/s) and frictional velocity (m2/s2)283 !284 ! for frazil ice285 INTEGER :: iter286 REAL(wp) :: zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp287 REAL(wp), PARAMETER :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used)288 REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness289 REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag)290 REAL(wp), PARAMETER :: zgamafr = 0.03_wp291 !!-----------------------------------------------------------------------292 !293 ! computation of friction velocity at T points294 IF( ln_icedyn ) THEN295 DO_2D( 0, 0, 0, 0 )296 zu_io = u_ice(ji ,jj ) - ssu_m(ji ,jj )297 zu_iom1 = u_ice(ji-1,jj ) - ssu_m(ji-1,jj )298 zv_io = v_ice(ji ,jj ) - ssv_m(ji ,jj )299 zv_iom1 = v_ice(ji ,jj-1) - ssv_m(ji ,jj-1)300 !301 zfric(ji,jj) = rn_cio * ( 0.5_wp * ( zu_io*zu_io + zu_iom1*zu_iom1 + zv_io*zv_io + zv_iom1*zv_iom1 ) ) * tmask(ji,jj,1)302 zvel (ji,jj) = 0.5_wp * SQRT( ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) * ( u_ice(ji-1,jj ) + u_ice(ji,jj) ) + &303 & ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) * ( v_ice(ji ,jj-1) + v_ice(ji,jj) ) )304 END_2D305 ELSE ! if no ice dynamics => transfer directly the atmospheric stress to the ocean306 DO_2D( 0, 0, 0, 0 )307 zfric(ji,jj) = r1_rho0 * SQRT( 0.5_wp * &308 & ( utau(ji,jj) * utau(ji,jj) + utau(ji-1,jj) * utau(ji-1,jj) &309 & + vtau(ji,jj) * vtau(ji,jj) + vtau(ji,jj-1) * vtau(ji,jj-1) ) ) * tmask(ji,jj,1)310 zvel(ji,jj) = 0._wp311 END_2D312 ENDIF313 CALL lbc_lnk( 'icethd', zfric, 'T', 1.0_wp, zvel, 'T', 1.0_wp )314 !315 !--------------------------------------------------------------------!316 ! Partial computation of forcing for the thermodynamic sea ice model317 !--------------------------------------------------------------------!318 DO_2D( nn_hls, nn_hls, nn_hls, nn_hls ) ! needed for qlead319 rswitch = tmask(ji,jj,1) * MAX( 0._wp , SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) ! 0 if no ice320 !321 ! --- Energy received in the lead from atm-oce exchanges, zqld is defined everywhere (J.m-2) --- !322 zqld = tmask(ji,jj,1) * rDt_ice * &323 & ( ( 1._wp - at_i_b(ji,jj) ) * qsr_oce(ji,jj) * frq_m(ji,jj) + &324 & ( 1._wp - at_i_b(ji,jj) ) * qns_oce(ji,jj) + qemp_oce(ji,jj) )325 326 ! --- Energy needed to bring ocean surface layer until its freezing, zqfr is defined everywhere (J.m-2) --- !327 ! (mostly<0 but >0 if supercooling)328 zqfr = rho0 * rcp * e3t_m(ji,jj) * ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) * tmask(ji,jj,1) ! both < 0 (t_bo < sst) and > 0 (t_bo > sst)329 zqfr_neg = MIN( zqfr , 0._wp ) ! only < 0330 zqfr_pos = MAX( zqfr , 0._wp ) ! only > 0331 332 ! --- Sensible ocean-to-ice heat flux (W/m2) --- !333 ! (mostly>0 but <0 if supercooling)334 zfric_u = MAX( SQRT( zfric(ji,jj) ), zfric_umin )335 qsb_ice_bot(ji,jj) = rswitch * rho0 * rcp * zch * zfric_u * ( ( sst_m(ji,jj) + rt0 ) - t_bo(ji,jj) )336 337 ! upper bound for qsb_ice_bot: the heat retrieved from the ocean must be smaller than the heat necessary to reach338 ! the freezing point, so that we do not have SST < T_freeze339 ! This implies: qsb_ice_bot(ji,jj) * at_i(ji,jj) * rtdice <= - zqfr_neg340 ! The following formulation is ok for both normal conditions and supercooling341 qsb_ice_bot(ji,jj) = rswitch * MIN( qsb_ice_bot(ji,jj), - zqfr_neg * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) )342 343 ! If conditions are always supercooled (such as at the mouth of ice-shelves), then ice grows continuously344 ! ==> stop ice formation by artificially setting up the turbulent fluxes to 0 when volume > 20m (arbitrary)345 IF( ( t_bo(ji,jj) - ( sst_m(ji,jj) + rt0 ) ) > 0._wp .AND. vt_i(ji,jj) >= 20._wp ) THEN346 zqfr = 0._wp347 zqfr_pos = 0._wp348 qsb_ice_bot(ji,jj) = 0._wp349 ENDIF350 !351 ! --- Energy Budget of the leads (qlead, J.m-2) --- !352 ! qlead is the energy received from the atm. in the leads.353 ! If warming (zqld >= 0), then the energy in the leads is used to melt ice (bottom melting) => fhld (W/m2)354 ! If cooling (zqld < 0), then the energy in the leads is used to grow ice in open water => qlead (J.m-2)355 IF( zqld >= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN356 ! upper bound for fhld: fhld should be equal to zqld357 ! but we have to make sure that this heat will not make the sst drop below the freezing point358 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr_pos359 ! The following formulation is ok for both normal conditions and supercooling360 fhld (ji,jj) = rswitch * MAX( 0._wp, ( zqld - zqfr_pos ) * r1_Dt_ice / MAX( at_i(ji,jj), epsi10 ) & ! divided by at_i since this is (re)multiplied by a_i in icethd_dh.F90361 & - qsb_ice_bot(ji,jj) )362 qlead(ji,jj) = 0._wp363 ELSE364 fhld (ji,jj) = 0._wp365 ! upper bound for qlead: qlead should be equal to zqld366 ! but before using this heat for ice formation, we suppose that the ocean cools down till the freezing point.367 ! The energy for this cooling down is zqfr. Also some heat will be removed from the ocean from turbulent fluxes (qsb)368 ! and freezing point is reached if zqfr = zqld - qsb*a/dt369 ! so the max heat that can be pulled out of the ocean is zqld - qsb - zqfr370 ! The following formulation is ok for both normal conditions and supercooling371 qlead(ji,jj) = MIN( 0._wp , zqld - ( qsb_ice_bot(ji,jj) * at_i(ji,jj) * rDt_ice ) - zqfr )372 ENDIF373 !374 ! If ice is landfast and ice concentration reaches its max375 ! => stop ice formation in open water376 IF( zvel(ji,jj) <= 5.e-04_wp .AND. at_i(ji,jj) >= rn_amax_2d(ji,jj)-epsi06 ) qlead(ji,jj) = 0._wp377 !378 ! If the grid cell is almost fully covered by ice (no leads)379 ! => stop ice formation in open water380 IF( at_i(ji,jj) >= (1._wp - epsi10) ) qlead(ji,jj) = 0._wp381 !382 ! If ln_leadhfx is false383 ! => do not use energy of the leads to melt sea-ice384 IF( .NOT.ln_leadhfx ) fhld(ji,jj) = 0._wp385 !386 END_2D387 388 ! In case we bypass open-water ice formation389 IF( .NOT. ln_icedO ) qlead(:,:) = 0._wp390 ! In case we bypass growing/melting from top and bottom391 IF( .NOT. ln_icedH ) THEN392 qsb_ice_bot(:,:) = 0._wp393 fhld (:,:) = 0._wp394 ENDIF395 396 !---------------------------------------------------------!397 ! Collection thickness of ice formed in leads and polynyas398 !---------------------------------------------------------!399 ! ht_i_new is the thickness of new ice formed in open water400 ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T)401 ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge402 ! where it forms aggregates of a specific thickness called collection thickness.403 !404 fraz_frac(:,:) = 0._wp405 !406 ! Default new ice thickness407 WHERE( qlead(:,:) < 0._wp ) ! cooling408 ht_i_new(:,:) = rn_hinew409 ELSEWHERE410 ht_i_new(:,:) = 0._wp411 END WHERE412 413 IF( ln_frazil ) THEN414 ztwogp = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) ) ! reduced grav415 !416 DO_2D( 0, 0, 0, 0 )417 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling418 ! -- Wind stress -- !419 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp420 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp421 ! Square root of wind stress422 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) )423 424 ! -- Frazil ice velocity -- !425 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) )426 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 )427 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 )428 429 ! -- Pack ice velocity -- !430 zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp431 zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp432 433 ! -- Relative frazil/pack ice velocity -- !434 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) )435 zvrel2 = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch436 437 ! -- fraction of frazil ice -- !438 fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz439 440 ! -- new ice thickness (iterative loop) -- !441 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1_wp ) &442 & / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) - zhicrit * zhicrit ) * ztwogp * zvrel2443 iter = 1444 DO WHILE ( iter < 20 )445 zf = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) - &446 & ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2447 zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2448 449 ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 )450 iter = iter + 1451 END DO452 !453 ! bound ht_i_new (though I don't see why it should be necessary)454 ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) )455 !456 ELSE457 ht_i_new(ji,jj) = 0._wp458 ENDIF459 !460 END_2D461 !462 CALL lbc_lnk( 'icethd', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp )463 464 ENDIF465 466 END SUBROUTINE thd_prep467 252 468 253 SUBROUTINE ice_thd_1d2d( kl, kn ) -
NEMO/trunk/src/ICE/icethd_do.F90
r15334 r15388 17 17 USE phycst ! physical constants 18 18 USE sbc_oce , ONLY : sss_m 19 USE sbc_ice , ONLY : utau_ice, vtau_ice 19 20 USE ice1D ! sea-ice: thermodynamics variables 20 21 USE ice ! sea-ice: variables … … 34 35 35 36 PUBLIC ice_thd_do ! called by ice_thd 37 PUBLIC ice_thd_frazil ! called by ice_thd 36 38 PUBLIC ice_thd_do_init ! called by ice_stp 37 38 ! !!** namelist (namthd_do) **39 REAL(wp), PUBLIC :: rn_hinew ! thickness for new ice formation (m)40 LOGICAL , PUBLIC :: ln_frazil ! use of frazil ice collection as function of wind (T) or not (F)41 REAL(wp), PUBLIC :: rn_maxfraz ! maximum portion of frazil ice collecting at the ice bottom42 REAL(wp), PUBLIC :: rn_vfraz ! threshold drift speed for collection of bottom frazil ice43 REAL(wp), PUBLIC :: rn_Cfraz ! squeezing coefficient for collection of bottom frazil ice44 39 45 40 !! * Substitutions … … 337 332 338 333 334 SUBROUTINE ice_thd_frazil 335 !!----------------------------------------------------------------------- 336 !! *** ROUTINE ice_thd_frazil *** 337 !! 338 !! ** Purpose : frazil ice collection thickness and fraction 339 !! 340 !! ** Inputs : u_ice, v_ice, utau_ice, vtau_ice 341 !! ** Ouputs : ht_i_new, fraz_frac 342 !!----------------------------------------------------------------------- 343 INTEGER :: ji, jj ! dummy loop indices 344 INTEGER :: iter 345 REAL(wp) :: zvfrx, zvgx, ztaux, zf, ztenagm, zvfry, zvgy, ztauy, zvrel2, zfp, ztwogp 346 REAL(wp), PARAMETER :: zcai = 1.4e-3_wp ! ice-air drag (clem: should be dependent on coupling/forcing used) 347 REAL(wp), PARAMETER :: zhicrit = 0.04_wp ! frazil ice thickness 348 REAL(wp), PARAMETER :: zsqcd = 1.0_wp / SQRT( 1.3_wp * zcai ) ! 1/SQRT(airdensity*drag) 349 REAL(wp), PARAMETER :: zgamafr = 0.03_wp 350 !!----------------------------------------------------------------------- 351 ! 352 !---------------------------------------------------------! 353 ! Collection thickness of ice formed in leads and polynyas 354 !---------------------------------------------------------! 355 ! ht_i_new is the thickness of new ice formed in open water 356 ! ht_i_new can be either prescribed (ln_frazil=F) or computed (ln_frazil=T) 357 ! Frazil ice forms in open water, is transported by wind, accumulates at the edge of the consolidated ice edge 358 ! where it forms aggregates of a specific thickness called collection thickness. 359 ! 360 fraz_frac(:,:) = 0._wp 361 ! 362 ! Default new ice thickness 363 WHERE( qlead(:,:) < 0._wp ) ! cooling 364 ht_i_new(:,:) = rn_hinew 365 ELSEWHERE 366 ht_i_new(:,:) = 0._wp 367 END WHERE 368 369 IF( ln_frazil ) THEN 370 ztwogp = 2._wp * rho0 / ( grav * 0.3_wp * ( rho0 - rhoi ) ) ! reduced grav 371 ! 372 DO_2D( 0, 0, 0, 0 ) 373 IF ( qlead(ji,jj) < 0._wp ) THEN ! cooling 374 ! -- Wind stress -- ! 375 ztaux = ( utau_ice(ji-1,jj ) * umask(ji-1,jj ,1) + utau_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 376 ztauy = ( vtau_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + vtau_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 377 ! Square root of wind stress 378 ztenagm = SQRT( SQRT( ztaux * ztaux + ztauy * ztauy ) ) 379 380 ! -- Frazil ice velocity -- ! 381 rswitch = MAX( 0._wp, SIGN( 1._wp , ztenagm - epsi10 ) ) 382 zvfrx = rswitch * zgamafr * zsqcd * ztaux / MAX( ztenagm, epsi10 ) 383 zvfry = rswitch * zgamafr * zsqcd * ztauy / MAX( ztenagm, epsi10 ) 384 385 ! -- Pack ice velocity -- ! 386 zvgx = ( u_ice(ji-1,jj ) * umask(ji-1,jj ,1) + u_ice(ji,jj) * umask(ji,jj,1) ) * 0.5_wp 387 zvgy = ( v_ice(ji ,jj-1) * vmask(ji ,jj-1,1) + v_ice(ji,jj) * vmask(ji,jj,1) ) * 0.5_wp 388 389 ! -- Relative frazil/pack ice velocity -- ! 390 rswitch = MAX( 0._wp, SIGN( 1._wp , at_i(ji,jj) - epsi10 ) ) 391 zvrel2 = MAX( (zvfrx - zvgx)*(zvfrx - zvgx) + (zvfry - zvgy)*(zvfry - zvgy), 0.15_wp*0.15_wp ) * rswitch 392 393 ! -- fraction of frazil ice -- ! 394 fraz_frac(ji,jj) = rswitch * ( TANH( rn_Cfraz * ( SQRT(zvrel2) - rn_vfraz ) ) + 1._wp ) * 0.5_wp * rn_maxfraz 395 396 ! -- new ice thickness (iterative loop) -- ! 397 ht_i_new(ji,jj) = zhicrit + ( zhicrit + 0.1_wp ) & 398 & / ( ( zhicrit + 0.1_wp ) * ( zhicrit + 0.1_wp ) - zhicrit * zhicrit ) * ztwogp * zvrel2 399 iter = 1 400 DO WHILE ( iter < 20 ) 401 zf = ( ht_i_new(ji,jj) - zhicrit ) * ( ht_i_new(ji,jj) * ht_i_new(ji,jj) - zhicrit * zhicrit ) - & 402 & ht_i_new(ji,jj) * zhicrit * ztwogp * zvrel2 403 zfp = ( ht_i_new(ji,jj) - zhicrit ) * ( 3.0_wp * ht_i_new(ji,jj) + zhicrit ) - zhicrit * ztwogp * zvrel2 404 405 ht_i_new(ji,jj) = ht_i_new(ji,jj) - zf / MAX( zfp, epsi20 ) 406 iter = iter + 1 407 END DO 408 ! 409 ! bound ht_i_new (though I don't see why it should be necessary) 410 ht_i_new(ji,jj) = MAX( 0.01_wp, MIN( ht_i_new(ji,jj), hi_max(jpl) ) ) 411 ! 412 ELSE 413 ht_i_new(ji,jj) = 0._wp 414 ENDIF 415 ! 416 END_2D 417 ! 418 CALL lbc_lnk( 'icethd_frazil', fraz_frac, 'T', 1.0_wp, ht_i_new, 'T', 1.0_wp ) 419 420 ENDIF 421 END SUBROUTINE ice_thd_frazil 422 423 339 424 SUBROUTINE ice_thd_do_init 340 425 !!----------------------------------------------------------------------- -
NEMO/trunk/src/ICE/icewri.F90
r14997 r15388 20 20 USE ice ! sea-ice: variables 21 21 USE icevar ! sea-ice: operations 22 USE icealb , ONLY : rn_alb_oce 22 23 ! 23 24 USE ioipsl ! … … 53 54 REAL(wp) :: z2da, z2db, zrho1, zrho2 54 55 REAL(wp) :: zmiss_val ! missing value retrieved from xios 55 REAL(wp), DIMENSION(jpi,jpj) :: z2d , zfast! 2D workspace56 REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace 56 57 REAL(wp), DIMENSION(jpi,jpj) :: zmsk00, zmsk05, zmsk15, zmsksn ! O%, 5% and 15% concentration mask and snow mask 57 58 REAL(wp), DIMENSION(jpi,jpj,jpl) :: zmsk00l, zmsksnl ! cat masks 59 REAL(wp), DIMENSION(:,:), ALLOCATABLE :: zfast, zalb, zmskalb ! 2D workspace 58 60 ! 59 61 ! Global ice diagnostics (SIMIP) … … 131 133 IF( iom_use('vice' ) ) CALL iom_put( 'vice' , v_ice ) ! ice velocity v 132 134 ! 133 IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN ! module of ice velocity 135 IF( iom_use('icevel') .OR. iom_use('fasticepres') ) THEN ! module of ice velocity & fast ice 136 ALLOCATE( zfast(jpi,jpj) ) 134 137 DO_2D( 0, 0, 0, 0 ) 135 138 z2da = u_ice(ji,jj) + u_ice(ji-1,jj) … … 144 147 END WHERE 145 148 CALL iom_put( 'fasticepres', zfast ) 146 ENDIF 147 149 DEALLOCATE( zfast ) 150 ENDIF 151 ! 152 IF( iom_use('icealb') .OR. iom_use('albedo') ) THEN ! ice albedo and surface albedo 153 ALLOCATE( zalb(jpi,jpj), zmskalb(jpi,jpj) ) 154 ! ice albedo 155 WHERE( at_i_b < 1.e-03 ) 156 zmskalb(:,:) = 0._wp 157 zalb (:,:) = rn_alb_oce 158 ELSEWHERE 159 zmskalb(:,:) = 1._wp 160 zalb (:,:) = SUM( alb_ice * a_i_b, dim=3 ) / at_i_b 161 END WHERE 162 CALL iom_put( 'icealb' , zalb * zmskalb + zmiss_val * ( 1._wp - zmskalb ) ) 163 ! ice+ocean albedo 164 zalb(:,:) = SUM( alb_ice * a_i_b, dim=3 ) + rn_alb_oce * ( 1._wp - at_i_b ) 165 CALL iom_put( 'albedo' , zalb ) 166 DEALLOCATE( zalb, zmskalb ) 167 ENDIF 168 ! 148 169 ! --- category-dependent fields --- ! 149 170 IF( iom_use('icemask_cat' ) ) CALL iom_put( 'icemask_cat' , zmsk00l ) ! ice mask 0%
Note: See TracChangeset
for help on using the changeset viewer.
