- Timestamp:
- 2020-12-01T23:04:18+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3_martin_ponds/src/ICE/icethd_pnd.F90
r13931 r13957 36 36 INTEGER :: nice_pnd ! choice of the type of pond scheme 37 37 ! ! associated indices: 38 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme39 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme40 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme38 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 41 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 42 42 43 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 44 zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m) 45 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 46 zvp_min = 1.e-4_wp ! minimum pond volume (m) 47 48 !-------------------------------------------------------------------------- 49 ! 50 ! Pond volume per area budget diags 51 ! 52 ! The idea of diags is the volume of ponds per grid cell area is 43 !-------------------------------------------------------------------------- 44 ! Diagnostics for pond volume per area 53 45 ! 54 46 ! dV/dt = mlt + drn + lid + rnf … … 59 51 ! 60 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 61 54 ! 62 ! In level mode, all terms are incorporated 63 64 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics 65 diag_dvpn_mlt, & !! meltwater pond volume input [m/day] 66 diag_dvpn_drn, & !! pond volume lost by drainage [m/day] 67 diag_dvpn_lid, & !! exchange with lid / refreezing [m/day] 68 diag_dvpn_rnf !! meltwater pond lost to runoff [m/day] 69 70 REAL(wp), ALLOCATABLE, DIMENSION(:) :: & ! pond volume budget diagnostics (1d) 71 diag_dvpn_mlt_1d, & !! meltwater pond volume input [m/day] 72 diag_dvpn_drn_1d, & !! pond volume lost by drainage [m/day] 73 diag_dvpn_lid_1d, & !! exchange with lid / refreezing [m/day] 74 diag_dvpn_rnf_1d !! meltwater pond lost to runoff [m/day] 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 [-] 75 63 76 64 !! * Substitutions … … 90 78 !! ** Purpose : change melt pond fraction and thickness 91 79 !! 92 !! Note: Melt ponds affect only radiative transfer for now 93 !! 94 !! No heat, no salt. 95 !! The melt water they carry is collected but 96 !! not removed from fw budget or released to the ocean 97 !! 98 !! A wfx_pnd has been coded for diagnostic purposes 99 !! It is not fully consistent yet. 100 !! 101 !! The current diagnostic lacks a contribution from drainage 102 !! 80 !! ** Note : Melt ponds affect only radiative transfer for now 81 !! No heat, no salt. 82 !! The current diagnostics lacks a contribution from drainage 103 83 !!------------------------------------------------------------------- 104 !! 84 INTEGER :: ji, jj, jl ! loop indices 85 !!------------------------------------------------------------------- 105 86 106 87 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 107 88 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 108 109 diag_dvpn_mlt (:,:) = 0._wp ; diag_dvpn_drn(:,:)= 0._wp110 diag_dvpn_lid (:,:) = 0._wp ; diag_dvpn_rnf(:,:)= 0._wp89 ! 90 diag_dvpn_mlt (:,:) = 0._wp ; diag_dvpn_drn (:,:) = 0._wp 91 diag_dvpn_lid (:,:) = 0._wp ; diag_dvpn_rnf (:,:) = 0._wp 111 92 diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp 112 93 diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp 113 94 114 SELECT CASE ( nice_pnd ) 95 !------------------------------------- 96 ! Remove ponds where ice has vanished 97 !------------------------------------- 98 at_i(:,:) = SUM( a_i, dim=3 ) 115 99 ! 116 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 117 ! 118 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 119 ! 120 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 121 ! 122 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 123 147 ! 124 125 IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 126 IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 127 IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 128 IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 129 130 DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 148 DEALLOCATE( diag_dvpn_mlt , diag_dvpn_lid , diag_dvpn_drn , diag_dvpn_rnf ) 131 149 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 132 150 133 151 END SUBROUTINE ice_thd_pnd 134 135 152 136 153 … … 151 168 !! ** References : Bush, G.W., and Trump, D.J. (2017) 152 169 !!------------------------------------------------------------------- 153 INTEGER :: ji ! loop indices 170 INTEGER :: ji, jl ! loop indices 171 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 154 172 !!------------------------------------------------------------------- 155 DO ji = 1, npti 156 ! 157 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 158 h_ip_1d(ji) = rn_hpnd 159 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 160 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 161 ELSE 162 h_ip_1d(ji) = 0._wp 163 a_ip_1d(ji) = 0._wp 164 h_il_1d(ji) = 0._wp 165 ENDIF 166 ! 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 167 211 END DO 168 212 ! … … 203 247 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 204 248 !! 205 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi 206 !! 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) 207 251 !! visc = water viscosity 208 252 !! Hp = height of top of the pond above sea-level 209 253 !! Hi = ice thickness thru which there is flushing 254 !! flush= correction otherwise flushing is excessive 210 255 !! 211 256 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness … … 218 263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 219 264 !! 220 !! ** References : Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 221 !! 222 !!------------------------------------------------------------------- 223 265 !! ** References : Flocco and Feltham (JGR, 2007) 266 !! Flocco et al (JGR, 2010) 267 !! Holland et al (J. Clim, 2012) 268 !! Hunke et al (OM 2012) 269 !!------------------------------------------------------------------- 224 270 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 225 271 !! … … 228 274 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 229 275 !! 230 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 231 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 232 279 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 233 280 REAL(wp) :: zv_ip_max ! max pond volume allowed 234 281 REAL(wp) :: zdT ! zTp-t_su 235 REAL(wp) :: zsbr 282 REAL(wp) :: zsbr, ztmelts ! Brine salinity 236 283 REAL(wp) :: zperm ! permeability of sea ice 237 284 REAL(wp) :: zfac, zdum ! temporary arrays 238 285 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 239 REAL(wp) :: zvold ! 240 !! 241 INTEGER :: ji, jj, jk, jl ! loop indices 242 286 !! 287 INTEGER :: ji, jk, jl ! loop indices 243 288 !!------------------------------------------------------------------- 244 245 289 z1_rhow = 1._wp / rhow 246 290 z1_aspect = 1._wp / zaspect 247 291 z1_Tp = 1._wp / zTp 248 292 249 !----------------------------------------------------------------------------------------------- 250 ! Identify grid cells with ice 251 !----------------------------------------------------------------------------------------------- 252 at_i(:,:) = SUM( a_i, dim=3 ) 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) 452 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 453 ! 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 456 ! 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 ! 468 END DO 253 469 ! 254 npti = 0 ; nptidx(:) = 0 255 DO_2D( 1, 1, 1, 1 ) 256 IF ( at_i(ji,jj) > epsi10 ) THEN 257 npti = npti + 1 258 nptidx( npti ) = (jj - 1) * jpi + ji 259 ENDIF 260 END_2D 261 262 !----------------------------------------------------------------------------------------------- 263 ! Prepare 1D arrays 264 !----------------------------------------------------------------------------------------------- 265 266 IF( npti > 0 ) THEN 267 268 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 269 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti) , wfx_sum ) 270 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti) , wfx_pnd ) 271 272 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti) , at_i ) 273 274 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 275 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 276 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 277 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 278 279 DO jl = 1, jpl 280 281 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 282 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 283 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 284 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 285 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 286 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 287 288 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum (1:npti), dh_i_sum_2d (:,:,jl) ) 289 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt (1:npti), dh_s_mlt_2d (:,:,jl) ) 290 291 DO jk = 1, nlay_i 292 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 293 END DO 294 295 !----------------------------------------------------------------------------------------------- 296 ! Go for ponds 297 !----------------------------------------------------------------------------------------------- 298 299 DO ji = 1, npti 300 ! !----------------------------------------------------! 301 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 302 ! !----------------------------------------------------! 303 !--- Remove ponds on thin ice or tiny ice fractions 304 a_ip_1d(ji) = 0._wp 305 h_ip_1d(ji) = 0._wp 306 h_il_1d(ji) = 0._wp 307 ! !--------------------------------! 308 ELSE ! Case ice thickness >= rn_himin ! 309 ! !--------------------------------! 310 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 311 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 312 ! 313 !------------------! 314 ! case ice melting ! 315 !------------------! 316 ! 317 !--- available meltwater for melt ponding ---! 318 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 319 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 320 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 321 322 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_Dt_ice ! surface melt input diag 323 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_Dt_ice ! runoff diag 324 ! 325 !--- overflow ---! 326 ! 327 ! 1) area driven overflow 328 ! 329 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 330 ! a_ip_max = zfr_mlt * a_i 331 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 332 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 333 zvold = zdv_mlt 334 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 335 336 ! 337 ! 2) depth driven overflow 338 ! 339 ! If pond depth exceeds half the ice thickness then reduce the pond volume 340 ! h_ip_max = 0.5 * h_i 341 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 342 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 343 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 344 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_Dt_ice ! runoff diag - overflow contribution 345 346 !--- Pond growing ---! 347 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 348 ! 349 !--- Lid melting ---! 350 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 351 ! 352 !--- mass flux ---! 353 ! MV I would recommend to remove that 354 ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 355 356 IF( zdv_mlt > 0._wp ) THEN 357 zfac = zdv_mlt * rhow * r1_Dt_ice ! melt pond mass flux < 0 [kg.m-2.s-1] 358 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 359 ! 360 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 361 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 362 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 363 ENDIF 364 365 !-------------------! 366 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 367 !-------------------! 368 ! 369 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 370 ! 371 !--- Pond contraction (due to refreezing) ---! 372 zvold = v_ip_1d(ji) ! for diag 373 374 IF( ln_pnd_lids ) THEN 375 ! 376 !--- Lid growing and subsequent pond shrinking ---! 377 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 378 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 379 380 ! Lid growing 381 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 382 383 ! Pond shrinking 384 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 385 386 ELSE 387 388 ! Pond shrinking 389 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 390 391 ENDIF 392 393 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_Dt_ice ! shrinking counted as lid diagnostic 394 395 ! 396 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 397 ! v_ip = h_ip * a_ip 398 ! 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) 399 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 400 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 401 402 !------------------------------------------------! 403 ! Pond drainage through brine network (flushing) ! 404 !------------------------------------------------! 405 ! height of top of the pond above sea-level 406 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 407 408 ! Calculate the permeability of the ice 409 ! mv- note here that a linear expression for permeability is fine, 410 ! simpler and more consistent with the rest of SI3 code 411 DO jk = 1, nlay_i 412 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 413 END DO 414 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 415 416 ! Do the drainage using Darcy's law 417 zdv_flush = - zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 418 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 419 ! zdv_flush = 0._wp ! MV remove pond drainage for now 420 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 421 422 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_Dt_ice ! shrinking counted as lid diagnostic 423 424 ! MV --- why pond drainage does not give back water into freshwater flux ? 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) > epsi10 ) 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 ENDIF 444 445 END DO ! ji 446 447 !----------------------------------------------------------------------------------------------- 448 ! Retrieve 2D arrays 449 !----------------------------------------------------------------------------------------------- 450 451 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 452 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 453 454 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 455 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 456 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 457 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 458 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 459 460 END DO ! ji 461 462 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 465 466 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 467 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 468 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 469 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 470 470 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd ) 471 471 ! 472 ENDIF 473 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 ! 474 477 END SUBROUTINE pnd_LEV 475 478 … … 507 510 !! 508 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 509 516 510 517 ! local variables … … 536 543 ! a_ip -> apond 537 544 ! a_ip_frac -> apnd 545 546 CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 538 547 539 548 !--------------------------------------------------------------- … … 628 637 DO_2D( 1, 1, 1, 1 ) 629 638 630 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN639 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 631 640 632 641 !-------------------------- … … 774 783 775 784 END DO ! jl 785 776 786 777 787 END SUBROUTINE pnd_TOPO … … 848 858 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 849 859 850 INTEGER :: ji, jj, jk, jl ! loop indices 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 851 864 852 865 a_ip(:,:,:) = 0._wp … … 855 868 DO_2D( 1, 1, 1, 1 ) 856 869 857 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN870 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 858 871 859 872 !-------------------------------------------------------------------
Note: See TracChangeset
for help on using the changeset viewer.