- Timestamp:
- 2021-11-28T18:59:49+01:00 (3 years ago)
- Location:
- NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icesbc.F90 (modified) (6 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk
- Property svn:externals
-
old new 9 9 10 10 # SETTE 11 ^/utils/CI/sette@14244 sette 11 ^/utils/CI/sette@HEAD sette 12
-
- Property svn:externals
-
NEMO/branches/2021/ticket2632_r14588_theta_sbcblk/src/ICE/icesbc.F90
r14433 r15548 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... 147 141 & sf(jp_slp)%fnow(:,:,1), sf(jp_qlw)%fnow(:,:,1), & 148 142 & sf(jp_prec)%fnow(:,:,1), sf(jp_snow)%fnow(:,:,1) ) 149 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )143 IF( ln_mixcpl ) CALL sbc_cpl_ice_flx( kt, picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 150 144 IF( nn_flxdist /= -1 ) CALL ice_flx_dist ( t_su, alb_ice, qns_ice, qsr_ice, dqns_ice, evap_ice, devap_ice, nn_flxdist ) 151 145 ! ! compute conduction flux and surface temperature (as in Jules surface module) … … 153 147 & CALL blk_ice_qcn ( ln_virtual_itd, t_su, t_bo, h_s, h_i ) 154 148 CASE ( jp_purecpl ) !--- coupled formulation 155 CALL sbc_cpl_ice_flx( picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i )149 CALL sbc_cpl_ice_flx( kt, picefr=at_i_b, palbi=alb_ice, psst=sst_m, pist=t_su, phs=h_s, phi=h_i ) 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 !!-------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.
