- Timestamp:
- 2022-01-19T19:34:39+01:00 (2 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.4_GO8_paquage_branch/src/ICE/icethd_pnd.F90
r14075 r15658 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 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 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 44 zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m) 45 zTd = 273.0_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 53 ! 54 ! dV/dt = mlt + drn + lid + rnf 55 ! mlt = input from surface melting 56 ! drn = drainage through brine network 57 ! lid = lid growth & melt 58 ! rnf = runoff (water directly removed out of surface melting + overflow) 59 ! 60 ! In topo mode, the pond water lost because it is in the snow is not included in the budget 61 ! 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] 75 ! 76 !-------------------------------------------------------------------------- 39 77 40 78 !! * Substitutions … … 52 90 !! 53 91 !! ** Purpose : change melt pond fraction and thickness 54 !! 92 !! 93 !! Note: Melt ponds affect only radiative transfer for now 94 !! 95 !! No heat, no salt. 96 !! The melt water they carry is collected but 97 !! not removed from fw budget or released to the ocean 98 !! 99 !! A wfx_pnd has been coded for diagnostic purposes 100 !! It is not fully consistent yet. 101 !! 102 !! The current diagnostic lacks a contribution from drainage 103 !! 55 104 !!------------------------------------------------------------------- 56 ! 105 !! 106 107 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 108 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 109 110 diag_dvpn_mlt(:,:) = 0._wp ; diag_dvpn_drn(:,:) = 0._wp 111 diag_dvpn_lid(:,:) = 0._wp ; diag_dvpn_rnf(:,:) = 0._wp 112 diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp 113 diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp 114 57 115 SELECT CASE ( nice_pnd ) 58 116 ! 59 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==!117 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 60 118 ! 61 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 119 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 120 ! 121 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 62 122 ! 63 123 END SELECT 64 124 ! 125 126 IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting 127 IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid 128 IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage 129 IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow 130 131 DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf ) 132 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 133 65 134 END SUBROUTINE ice_thd_pnd 66 135 … … 75 144 !! to non-zero values when t_su = 0C 76 145 !! 77 !! ** Tunable parameters : pond fraction (rn_apnd), ponddepth (rn_hpnd)146 !! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd) 78 147 !! 79 148 !! ** Note : Coupling with such melt ponds is only radiative … … 145 214 !! a_ip/a_i = a_ip_frac = h_ip / zaspect 146 215 !! 147 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min216 !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 148 217 !! 149 !! ** Note : mostly stolen from CICE 150 !! 151 !! ** References : Flocco and Feltham (JGR, 2007) 152 !! Flocco et al (JGR, 2010) 153 !! Holland et al (J. Clim, 2012) 218 !! ** Note : Mostly stolen from CICE but not only 219 !! 220 !! ** References : Holland et al (J. Clim, 2012), Hunke et al (OM 2012) 221 !! 154 222 !!------------------------------------------------------------------- 155 223 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array … … 168 236 REAL(wp) :: zfac, zdum ! temporary arrays 169 237 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 170 !! 171 INTEGER :: ji, jk ! loop indices 238 REAL(wp) :: zvold ! 239 !! 240 INTEGER :: ji, jj, jk, jl ! loop indices 241 172 242 !!------------------------------------------------------------------- 173 243 z1_rhow = 1._wp / rhow 174 244 z1_aspect = 1._wp / zaspect 175 245 z1_Tp = 1._wp / zTp 176 177 DO ji = 1, npti 178 ! !----------------------------------------------------! 179 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 180 ! !----------------------------------------------------! 181 !--- Remove ponds on thin ice or tiny ice fractions 182 a_ip_1d(ji) = 0._wp 183 h_ip_1d(ji) = 0._wp 184 h_il_1d(ji) = 0._wp 185 ! !--------------------------------! 186 ELSE ! Case ice thickness >= rn_himin ! 187 ! !--------------------------------! 188 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 189 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 190 ! 191 !------------------! 192 ! case ice melting ! 193 !------------------! 194 ! 195 !--- available meltwater for melt ponding ---! 196 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 197 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 198 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 199 ! 200 !--- overflow ---! 201 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 202 ! a_ip_max = zfr_mlt * a_i 203 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 204 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 205 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 206 207 ! If pond depth exceeds half the ice thickness then reduce the pond volume 208 ! h_ip_max = 0.5 * h_i 209 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 210 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 211 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 246 247 !----------------------------------------------------------------------------------------------- 248 ! Identify grid cells with ice 249 !----------------------------------------------------------------------------------------------- 250 at_i(:,:) = SUM( a_i, dim=3 ) 251 ! 252 npti = 0 ; nptidx(:) = 0 253 DO jj = 1, jpj 254 DO ji = 1, jpi 255 IF ( at_i(ji,jj) > epsi10 ) THEN 256 npti = npti + 1 257 nptidx( npti ) = (jj - 1) * jpi + ji 258 ENDIF 259 END DO 260 END DO 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), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 273 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 274 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 275 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 276 277 DO jl = 1, jpl 278 279 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 280 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 281 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 282 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 283 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 284 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 285 286 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum (1:npti), dh_i_sum_2d (:,:,jl) ) 287 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt (1:npti), dh_s_mlt_2d (:,:,jl) ) 212 288 213 !--- Pond growing ---! 214 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 215 ! 216 !--- Lid melting ---! 217 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 218 ! 219 !--- mass flux ---! 220 IF( zdv_mlt > 0._wp ) THEN 221 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 222 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 223 ! 224 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 225 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 226 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 227 ENDIF 228 229 !-------------------! 230 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 231 !-------------------! 232 ! 233 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 234 ! 235 !--- Pond contraction (due to refreezing) ---! 236 IF( ln_pnd_lids ) THEN 237 ! 238 !--- Lid growing and subsequent pond shrinking ---! 239 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 240 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 241 242 ! Lid growing 243 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 244 245 ! Pond shrinking 246 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 247 248 ELSE 249 ! Pond shrinking 250 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 251 ENDIF 252 ! 253 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 254 ! v_ip = h_ip * a_ip 255 ! 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) 256 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 257 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 258 259 !---------------! 260 ! Pond flushing ! 261 !---------------! 262 ! height of top of the pond above sea-level 263 zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 289 DO jk = 1, nlay_i 290 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 291 END DO 264 292 265 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 266 DO jk = 1, nlay_i 267 zsbr = - 1.2_wp & 268 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 269 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 270 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 271 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 272 END DO 273 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 274 275 ! Do the drainage using Darcy's law 276 zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 277 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 278 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 279 280 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 281 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 282 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 283 284 !--- Corrections and lid thickness ---! 285 IF( ln_pnd_lids ) THEN 286 !--- retrieve lid thickness from volume ---! 287 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 288 ELSE ; h_il_1d(ji) = 0._wp 289 ENDIF 290 !--- remove ponds if lids are much larger than ponds ---! 291 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 293 !----------------------------------------------------------------------------------------------- 294 ! Go for ponds 295 !----------------------------------------------------------------------------------------------- 296 297 DO ji = 1, npti 298 ! !----------------------------------------------------! 299 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 300 ! !----------------------------------------------------! 301 !--- Remove ponds on thin ice or tiny ice fractions 292 302 a_ip_1d(ji) = 0._wp 293 303 h_ip_1d(ji) = 0._wp 294 304 h_il_1d(ji) = 0._wp 305 ! !--------------------------------! 306 ELSE ! Case ice thickness >= rn_himin ! 307 ! !--------------------------------! 308 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 309 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 310 ! 311 !------------------! 312 ! case ice melting ! 313 !------------------! 314 ! 315 !--- available meltwater for melt ponding ---! 316 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 317 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 318 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 319 320 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_rdtice ! surface melt input diag 321 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_rdtice ! runoff diag 322 ! 323 !--- overflow ---! 324 ! 325 ! 1) area driven overflow 326 ! 327 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume 328 ! a_ip_max = zfr_mlt * a_i 329 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 330 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 331 zvold = zdv_mlt 332 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 333 334 ! 335 ! 2) depth driven overflow 336 ! 337 ! If pond depth exceeds half the ice thickness then reduce the pond volume 338 ! h_ip_max = 0.5 * h_i 339 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 340 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 341 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 342 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_rdtice ! runoff diag - overflow contribution 343 344 !--- Pond growing ---! 345 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 346 ! 347 !--- Lid melting ---! 348 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 349 ! 350 !--- mass flux ---! 351 ! MV I would recommend to remove that 352 ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 353 354 IF( zdv_mlt > 0._wp ) THEN 355 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 356 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 357 ! 358 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 359 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 360 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 361 ENDIF 362 363 !-------------------! 364 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 365 !-------------------! 366 ! 367 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 368 ! 369 !--- Pond contraction (due to refreezing) ---! 370 zvold = v_ip_1d(ji) ! for diag 371 372 IF( ln_pnd_lids ) THEN 373 ! 374 !--- Lid growing and subsequent pond shrinking ---! 375 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 376 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 377 378 ! Lid growing 379 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 380 381 ! Pond shrinking 382 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 383 384 ELSE 385 386 ! Pond shrinking 387 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 388 389 ENDIF 390 391 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_rdtice ! shrinking counted as lid diagnostic 392 393 ! 394 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 395 ! v_ip = h_ip * a_ip 396 ! 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) 397 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 398 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 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) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 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 ! MV linear expression more consistent & simpler zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 414 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 415 END DO 416 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 417 418 ! Do the drainage using Darcy's law 419 zdv_flush = - zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush 420 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 421 ! zdv_flush = 0._wp ! MV remove pond drainage for now 422 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 423 424 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_rdtice ! shrinking counted as lid diagnostic 425 426 ! MV --- why pond drainage does not give back water into freshwater flux ? 427 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 428 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 429 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 430 431 !--- Corrections and lid thickness ---! 432 IF( ln_pnd_lids ) THEN 433 !--- retrieve lid thickness from volume ---! 434 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 435 ELSE ; h_il_1d(ji) = 0._wp 436 ENDIF 437 !--- remove ponds if lids are much larger than ponds ---! 438 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 439 a_ip_1d(ji) = 0._wp 440 h_ip_1d(ji) = 0._wp 441 h_il_1d(ji) = 0._wp 442 ENDIF 443 ENDIF 444 ! 295 445 ENDIF 446 447 END DO ! ji 448 449 !----------------------------------------------------------------------------------------------- 450 ! Retrieve 2D arrays 451 !----------------------------------------------------------------------------------------------- 452 453 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 454 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 455 456 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 457 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 458 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 459 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 460 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 461 DO jk = 1, nlay_i 462 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 463 END DO 464 465 END DO ! ji 466 467 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 468 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 469 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 470 471 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 472 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 473 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 474 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 475 476 ! 477 ENDIF 478 479 END SUBROUTINE pnd_LEV 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 513 ! local variables 514 REAL(wp) :: & 515 zdHui, & ! change in thickness of ice lid (m) 516 zomega, & ! conduction 517 zdTice, & ! temperature difference across ice lid (C) 518 zdvice, & ! change in ice volume (m) 519 zTavg, & ! mean surface temperature across categories (C) 520 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 521 zTp, & ! pond freezing temperature (C) 522 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 523 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 524 z1_rhow, & ! inverse water density 525 zv_pnd , & ! volume of meltwater contributing to ponds 526 zv_mlt ! total amount of meltwater produced 527 528 REAL(wp), DIMENSION(jpi,jpj) :: zvolp_ini, & !! total melt pond water available before redistribution and drainage 529 zvolp , & !! total melt pond water volume 530 zvolp_res !! remaining melt pond water available after drainage 531 532 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 533 534 INTEGER :: ji, jj, jk, jl ! loop indices 535 536 INTEGER :: i_test 537 538 ! Note 539 ! equivalent for CICE translation 540 ! a_ip -> apond 541 ! a_ip_frac -> apnd 542 543 !--------------------------------------------------------------- 544 ! Initialise 545 !--------------------------------------------------------------- 546 547 ! Parameters & constants (move to parameters) 548 zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) 549 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 550 z1_rhow = 1._wp / rhow 551 552 ! Set required ice variables (hard-coded here for now) 553 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 554 555 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 556 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 557 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 558 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 559 560 ! thickness 561 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 562 ELSEWHERE ; z1_a_i(:,:,:) = 0._wp 563 END WHERE 564 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 565 566 !--------------------------------------------------------------- 567 ! Change 2D to 1D 568 !--------------------------------------------------------------- 569 ! MV 570 ! a less computing-intensive version would have 2D-1D passage here 571 ! use what we have in iceitd.F90 (incremental remapping) 572 573 !-------------------------------------------------------------- 574 ! Collect total available pond water volume 575 !-------------------------------------------------------------- 576 ! Assuming that meltwater (+rain in principle) runsoff the surface 577 ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 578 ! I cite her words, they are very talkative 579 ! "grid cells with very little ice cover (and hence more open water area) 580 ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 581 ! "This results in the same runoff fraction r for each ice category within a grid cell" 582 583 zvolp(:,:) = 0._wp 584 zvolp_res(:,:) = 0._wp 585 586 DO jl = 1, jpl 587 DO jj = 1, jpj 588 DO ji = 1, jpi 589 590 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 591 592 !--- Available and contributing meltwater for melt ponding ---! 593 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 594 & * z1_rhow * a_i(ji,jj,jl) 595 ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 596 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds 597 zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl 598 599 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice ! diags 600 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_rdtice 601 602 !--- Create possible new ponds 603 ! if pond does not exist, create new pond over full ice area 604 !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 605 IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 606 a_ip(ji,jj,jl) = a_i(ji,jj,jl) 607 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 608 ENDIF 609 610 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 611 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 612 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 613 614 !--- Total available pond water volume (pre-existing + newly produced)j 615 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 616 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 617 618 ENDIF ! a_i 619 620 END DO! jl 621 END DO ! jj 622 END DO ! ji 623 624 zvolp_ini(:,:) = zvolp(:,:) 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 jj = 1, jpj 638 DO ji = 1, jpi 639 640 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp_ini(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 641 642 !-------------------------- 643 ! Pond lid growth and melt 644 !-------------------------- 645 ! Mean surface temperature 646 zTavg = 0._wp 647 DO jl = 1, jpl 648 zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 649 END DO 650 zTavg = zTavg / at_i(ji,jj) !!! could get a division by zero here 651 652 DO jl = 1, jpl-1 653 654 IF ( v_il(ji,jj,jl) > epsi10 ) THEN 655 656 !---------------------------------------------------------------- 657 ! Lid melting: floating upper ice layer melts in whole or part 658 !---------------------------------------------------------------- 659 ! Use Tsfc for each category 660 IF ( t_su(ji,jj,jl) > zTp ) THEN 661 662 zdvice = MIN( -dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 663 664 IF ( zdvice > epsi10 ) THEN 665 666 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 667 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 668 ! as it is already counted in surface melt 669 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 670 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 671 672 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 673 ! ice lid melted and category is pond covered 674 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 675 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 676 v_il(ji,jj,jl) = 0._wp 677 ENDIF 678 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 679 680 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 681 682 ENDIF 683 684 !---------------------------------------------------------------- 685 ! Freeze pre-existing lid 686 !---------------------------------------------------------------- 687 688 ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 689 690 ! differential growth of base of surface floating ice layer 691 ! zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 ! line valid for temp in C 692 zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0 693 zomega = rcnd_i * zdTice / zrhoi_L 694 zdHui = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 695 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 696 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 697 698 IF ( zdvice > epsi10 ) THEN 699 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 700 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 701 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 702 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 703 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 704 705 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 706 707 ENDIF 708 709 ENDIF ! Tsfcn(i,j,n) 710 711 !---------------------------------------------------------------- 712 ! Freeze new lids 713 !---------------------------------------------------------------- 714 ! upper ice layer begins to form 715 ! note: albedo does not change 716 717 ELSE ! v_il < epsi10 718 719 ! thickness of newly formed ice 720 ! the surface temperature of a meltpond is the same as that 721 ! of the ice underneath (0C), and the thermodynamic surface 722 ! flux is the same 723 724 !!! we need net surface energy flux, excluding conduction 725 !!! fsurf is summed over categories in CICE 726 !!! we have the category-dependent flux, let us use it ? 727 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 728 zdHui = MAX ( -zfsurf * rdt_ice / zrhoi_L , 0._wp ) 729 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 730 731 IF ( zdvice > epsi10 ) THEN 732 733 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 734 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 735 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 736 737 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 738 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 739 740 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 741 ENDIF 742 743 ENDIF ! v_il 744 745 !---------------------------------------------------------------- 746 ! Remove ice lid if there is no liquid pond 747 !---------------------------------------------------------------- 748 749 IF (v_ip(ji,jj,jl) < epsi10) THEN 750 v_il(ji,jj,jl) = 0._wp 751 ENDIF 752 753 END DO ! jl 754 755 ELSE ! remove ponds on thin ice 756 757 v_ip(ji,jj,:) = 0._wp 758 v_il(ji,jj,:) = 0._wp 759 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 760 ! zvolp(ji,jj) = 0._wp 761 296 762 ENDIF 297 ! 298 ENDIF 299 300 END DO 301 ! 302 END SUBROUTINE pnd_LEV 303 763 764 END DO ! jj 765 END DO ! ji 766 767 ENDIF ! ln_pnd_lids 768 769 !--------------------------------------------------------------- 770 ! Clean-up variables (probably duplicates what icevar would do) 771 !--------------------------------------------------------------- 772 ! MV comment 773 ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 774 ! icevar should recompute all other variables (if needed at all) 775 776 DO jl = 1, jpl 777 778 DO jj = 1, jpj 779 DO ji = 1, jpi 780 781 ! ! zap lids on small ponds 782 ! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 783 ! .AND. v_il(ji,jj,jl) > epsi10) THEN 784 ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 785 ! ENDIF 786 787 ! recalculate equivalent pond variables 788 IF ( a_ip(ji,jj,jl) > epsi10) THEN 789 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 790 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 791 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 792 ENDIF 793 ! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 794 ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 795 ! ENDIF 796 797 END DO ! ji 798 END DO ! jj 799 800 END DO ! jl 801 802 END SUBROUTINE pnd_TOPO 803 804 805 SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 806 807 !!------------------------------------------------------------------- 808 !! *** ROUTINE ice_thd_pnd_area *** 809 !! 810 !! ** Purpose : Given the total volume of available pond water, 811 !! redistribute and drain water 812 !! 813 !! ** Method 814 !! 815 !-----------| 816 ! | 817 ! |-----------| 818 !___________|___________|______________________________________sea-level 819 ! | | 820 ! | |---^--------| 821 ! | | | | 822 ! | | | |-----------| |------- 823 ! | | | alfan | | | 824 ! | | | | |--------------| 825 ! | | | | | | 826 !---------------------------v------------------------------------------- 827 ! | | ^ | | | 828 ! | | | | |--------------| 829 ! | | | betan | | | 830 ! | | | |-----------| |------- 831 ! | | | | 832 ! | |---v------- | 833 ! | | 834 ! |-----------| 835 ! | 836 !-----------| 837 ! 838 !! 839 !!------------------------------------------------------------------ 840 841 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 842 zvolp, & ! total available pond water 843 zdvolp ! remaining meltwater after redistribution 844 845 INTEGER :: & 846 ns, & 847 m_index, & 848 permflag 849 850 REAL (wp), DIMENSION(jpl) :: & 851 hicen, & 852 hsnon, & 853 asnon, & 854 alfan, & 855 betan, & 856 cum_max_vol, & 857 reduced_aicen 858 859 REAL (wp), DIMENSION(0:jpl) :: & 860 cum_max_vol_tmp 861 862 REAL (wp) :: & 863 hpond, & 864 drain, & 865 floe_weight, & 866 pressure_head, & 867 hsl_rel, & 868 deltah, & 869 perm, & 870 msno 871 872 REAL (wp), parameter :: & 873 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 874 875 INTEGER :: ji, jj, jk, jl ! loop indices 876 877 a_ip(:,:,:) = 0._wp 878 h_ip(:,:,:) = 0._wp 879 880 DO jj = 1, jpj 881 DO ji = 1, jpi 882 883 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 884 885 !------------------------------------------------------------------- 886 ! initialize 887 !------------------------------------------------------------------- 888 889 DO jl = 1, jpl 890 891 !---------------------------------------- 892 ! compute the effective snow fraction 893 !---------------------------------------- 894 895 IF (a_i(ji,jj,jl) < epsi10) THEN 896 hicen(jl) = 0._wp 897 hsnon(jl) = 0._wp 898 reduced_aicen(jl) = 0._wp 899 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 900 ELSE 901 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 902 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 903 reduced_aicen(jl) = 1._wp ! n=jpl 904 905 !js: initial code in NEMO_DEV 906 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 907 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 908 909 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 910 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 911 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 912 913 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 914 ! MV should check whether this makes sense to have the same effective snow fraction in here 915 ! OLI: it probably doesn't 916 END IF 917 918 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 919 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 920 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 921 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 922 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 923 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 924 ! alfan = 60% of the ice volume) in each category lies above the reference line, 925 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 926 927 ! MV: 928 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 929 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 930 931 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 932 933 alfan(jl) = 0.6 * hicen(jl) 934 betan(jl) = 0.4 * hicen(jl) 935 936 cum_max_vol(jl) = 0._wp 937 cum_max_vol_tmp(jl) = 0._wp 938 939 END DO ! jpl 940 941 cum_max_vol_tmp(0) = 0._wp 942 drain = 0._wp 943 zdvolp(ji,jj) = 0._wp 944 945 !---------------------------------------------------------- 946 ! Drain overflow water, update pond fraction and volume 947 !---------------------------------------------------------- 948 949 !-------------------------------------------------------------------------- 950 ! the maximum amount of water that can be contained up to each ice category 951 !-------------------------------------------------------------------------- 952 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 953 ! Then the excess volume cum_max_vol(jl) drains out of the system 954 ! It should be added to wfx_pnd_out 955 956 DO jl = 1, jpl-1 ! last category can not hold any volume 957 958 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 959 960 ! total volume in level including snow 961 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 962 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 963 964 ! subtract snow solid volumes from lower categories in current level 965 DO ns = 1, jl 966 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 967 - rhos/rhow * & ! free air fraction that can be filled by water 968 asnon(ns) * & ! effective areal fraction of snow in that category 969 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 970 END DO 971 972 ELSE ! assume higher categories unoccupied 973 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 974 END IF 975 !IF (cum_max_vol_tmp(jl) < z0) THEN 976 ! CALL abort_ice('negative melt pond volume') 977 !END IF 978 END DO 979 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 980 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 981 982 !---------------------------------------------------------------- 983 ! is there more meltwater than can be held in the floe? 984 !---------------------------------------------------------------- 985 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 986 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 987 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 988 989 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 990 991 zdvolp(ji,jj) = drain ! this is the drained water 992 IF (zvolp(ji,jj) < epsi10) THEN 993 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 994 zvolp(ji,jj) = 0._wp 995 END IF 996 END IF 997 998 ! height and area corresponding to the remaining volume 999 ! routine leaves zvolp unchanged 1000 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1001 1002 DO jl = 1, m_index 1003 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 1004 ! ! volume instead, no ? 1005 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 1006 a_ip(ji,jj,jl) = reduced_aicen(jl) 1007 ! in practise, pond fraction depends on the empirical snow fraction 1008 ! so in turn on ice thickness 1009 END DO 1010 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1011 1012 !------------------------------------------------------------------------ 1013 ! Drainage through brine network (permeability) 1014 !------------------------------------------------------------------------ 1015 !!! drainage due to ice permeability - Darcy's law 1016 1017 ! sea water level 1018 msno = 0._wp 1019 DO jl = 1 , jpl 1020 msno = msno + v_s(ji,jj,jl) * rhos 1021 END DO 1022 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rau0*zvolp(ji,jj) ) / at_i(ji,jj) 1023 hsl_rel = floe_weight / rau0 & 1024 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1025 1026 deltah = hpond - hsl_rel 1027 pressure_head = grav * rau0 * max(deltah, 0._wp) 1028 1029 ! drain if ice is permeable 1030 permflag = 0 1031 1032 IF (pressure_head > 0._wp) THEN 1033 DO jl = 1, jpl-1 1034 IF ( hicen(jl) /= 0._wp ) THEN 1035 1036 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1037 1038 perm = 0._wp ! MV ugly dummy patch 1039 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1040 IF (perm > 0._wp) permflag = 1 1041 1042 drain = perm*a_ip(ji,jj,jl)*pressure_head*rdt_ice / & 1043 (viscosity*hicen(jl)) 1044 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1045 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1046 1047 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1048 1049 IF (zvolp(ji,jj) < epsi10) THEN 1050 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 1051 zvolp(ji,jj) = 0._wp 1052 END IF 1053 END IF 1054 END DO 1055 1056 ! adjust melt pond dimensions 1057 IF (permflag > 0) THEN 1058 ! recompute pond depth 1059 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1060 DO jl = 1, m_index 1061 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 1062 a_ip(ji,jj,jl) = reduced_aicen(jl) 1063 END DO 1064 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1065 END IF 1066 END IF ! pressure_head 1067 1068 !------------------------------- 1069 ! remove water from the snow 1070 !------------------------------- 1071 !------------------------------------------------------------------------ 1072 ! total melt pond volume in category does not include snow volume 1073 ! snow in melt ponds is not melted 1074 !------------------------------------------------------------------------ 1075 1076 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1077 ! how much, so I did not diagnose it 1078 ! so if there is a problem here, nobody is going to see it... 1079 1080 1081 ! Calculate pond volume for lower categories 1082 DO jl = 1,m_index-1 1083 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 1084 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1085 END DO 1086 1087 ! Calculate pond volume for highest category = remaining pond volume 1088 1089 ! The following is completely unclear to Martin at least 1090 ! Could we redefine properly and recode in a more readable way ? 1091 1092 ! m_index = last category with melt pond 1093 1094 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 1095 1096 IF (m_index > 1) THEN 1097 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1098 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 1099 ELSE 1100 v_ip(ji,jj,m_index) = 0._wp 1101 h_ip(ji,jj,m_index) = 0._wp 1102 a_ip(ji,jj,m_index) = 0._wp 1103 ! If remaining pond volume is negative reduce pond volume of 1104 ! lower category 1105 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 1106 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) 1107 END IF 1108 END IF 1109 1110 DO jl = 1,m_index 1111 IF (a_ip(ji,jj,jl) > epsi10) THEN 1112 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 1113 ELSE 1114 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1115 h_ip(ji,jj,jl) = 0._wp 1116 v_ip(ji,jj,jl) = 0._wp 1117 a_ip(ji,jj,jl) = 0._wp 1118 END IF 1119 END DO 1120 DO jl = m_index+1, jpl 1121 h_ip(ji,jj,jl) = 0._wp 1122 a_ip(ji,jj,jl) = 0._wp 1123 v_ip(ji,jj,jl) = 0._wp 1124 END DO 1125 1126 ENDIF 1127 END DO ! ji 1128 END DO ! jj 1129 1130 END SUBROUTINE ice_thd_pnd_area 1131 1132 1133 SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 1134 !!------------------------------------------------------------------- 1135 !! *** ROUTINE ice_thd_pnd_depth *** 1136 !! 1137 !! ** Purpose : Compute melt pond depth 1138 !!------------------------------------------------------------------- 1139 1140 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1141 aicen, & 1142 asnon, & 1143 hsnon, & 1144 alfan, & 1145 cum_max_vol 1146 1147 REAL (wp), INTENT(IN) :: & 1148 zvolp 1149 1150 REAL (wp), INTENT(OUT) :: & 1151 hpond 1152 1153 INTEGER, INTENT(OUT) :: & 1154 m_index 1155 1156 INTEGER :: n, ns 1157 1158 REAL (wp), DIMENSION(0:jpl+1) :: & 1159 hitl, & 1160 aicetl 1161 1162 REAL (wp) :: & 1163 rem_vol, & 1164 area, & 1165 vol, & 1166 tmp, & 1167 z0 = 0.0_wp 1168 1169 !---------------------------------------------------------------- 1170 ! hpond is zero if zvolp is zero - have we fully drained? 1171 !---------------------------------------------------------------- 1172 1173 IF (zvolp < epsi10) THEN 1174 hpond = z0 1175 m_index = 0 1176 ELSE 1177 1178 !---------------------------------------------------------------- 1179 ! Calculate the category where water fills up to 1180 !---------------------------------------------------------------- 1181 1182 !----------| 1183 ! | 1184 ! | 1185 ! |----------| -- -- 1186 !__________|__________|_________________________________________ ^ 1187 ! | | rem_vol ^ | Semi-filled 1188 ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1189 ! | | | | 1190 ! | | | |hpond 1191 ! | | |----------| | |------- 1192 ! | | | | | | 1193 ! | | | |---v-----| 1194 ! | | m_index | | | 1195 !------------------------------------------------------------- 1196 1197 m_index = 0 ! 1:m_index categories have water in them 1198 DO n = 1, jpl 1199 IF (zvolp <= cum_max_vol(n)) THEN 1200 m_index = n 1201 IF (n == 1) THEN 1202 rem_vol = zvolp 1203 ELSE 1204 rem_vol = zvolp - cum_max_vol(n-1) 1205 END IF 1206 exit ! to break out of the loop 1207 END IF 1208 END DO 1209 m_index = min(jpl-1, m_index) 1210 1211 !---------------------------------------------------------------- 1212 ! semi-filled layer may have m_index different snow in it 1213 !---------------------------------------------------------------- 1214 1215 !----------------------------------------------------------- ^ 1216 ! | alfan(m_index+1) 1217 ! | 1218 !hitl(3)--> |----------| | 1219 !hitl(2)--> |------------| * * * * *| | 1220 !hitl(1)--> |----------|* * * * * * |* * * * * | | 1221 !hitl(0)-->------------------------------------------------- | ^ 1222 ! various snow from lower categories | |alfa(m_index) 1223 1224 ! hitl - heights of the snow layers from thinner and current categories 1225 ! aicetl - area of each snow depth in this layer 1226 1227 hitl(:) = z0 1228 aicetl(:) = z0 1229 DO n = 1, m_index 1230 hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1231 alfan(m_index+1) - alfan(m_index)), z0) 1232 aicetl(n) = asnon(n) 1233 1234 aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1235 END DO 1236 1237 hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1238 aicetl(m_index+1) = z0 1239 1240 !---------------------------------------------------------------- 1241 ! reorder array according to hitl 1242 ! snow heights not necessarily in height order 1243 !---------------------------------------------------------------- 1244 1245 DO ns = 1, m_index+1 1246 DO n = 0, m_index - ns + 1 1247 IF (hitl(n) > hitl(n+1)) THEN ! swap order 1248 tmp = hitl(n) 1249 hitl(n) = hitl(n+1) 1250 hitl(n+1) = tmp 1251 tmp = aicetl(n) 1252 aicetl(n) = aicetl(n+1) 1253 aicetl(n+1) = tmp 1254 END IF 1255 END DO 1256 END DO 1257 1258 !---------------------------------------------------------------- 1259 ! divide semi-filled layer into set of sublayers each vertically homogenous 1260 !---------------------------------------------------------------- 1261 1262 !hitl(3)---------------------------------------------------------------- 1263 ! | * * * * * * * * 1264 ! |* * * * * * * * * 1265 !hitl(2)---------------------------------------------------------------- 1266 ! | * * * * * * * * | * * * * * * * * 1267 ! |* * * * * * * * * |* * * * * * * * * 1268 !hitl(1)---------------------------------------------------------------- 1269 ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1270 ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1271 !hitl(0)---------------------------------------------------------------- 1272 ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1273 1274 ! move up over layers incrementing volume 1275 DO n = 1, m_index+1 1276 1277 area = sum(aicetl(:)) - & ! total area of sub-layer 1278 (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1279 1280 vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1281 1282 IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1283 hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 1284 1285 exit 1286 ELSE ! still in sub-layer below the sub-layer with the depth 1287 rem_vol = rem_vol - vol 1288 END IF 1289 1290 END DO 1291 1292 END IF 1293 1294 END SUBROUTINE ice_thd_pnd_depth 1295 1296 1297 SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 1298 !!------------------------------------------------------------------- 1299 !! *** ROUTINE ice_thd_pnd_perm *** 1300 !! 1301 !! ** Purpose : Determine the liquid fraction of brine in the ice 1302 !! and its permeability 1303 !!------------------------------------------------------------------- 1304 1305 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1306 ticen, & ! internal ice temperature (K) 1307 salin ! salinity (ppt) !js: ppt according to cice 1308 1309 REAL (wp), INTENT(OUT) :: & 1310 perm ! permeability 1311 1312 REAL (wp) :: & 1313 Sbr ! brine salinity 1314 1315 REAL (wp), DIMENSION(nlay_i) :: & 1316 Tin, & ! ice temperature 1317 phi ! liquid fraction 1318 1319 INTEGER :: k 1320 1321 !----------------------------------------------------------------- 1322 ! Compute ice temperatures from enthalpies using quadratic formula 1323 !----------------------------------------------------------------- 1324 1325 DO k = 1,nlay_i 1326 Tin(k) = ticen(k) - rt0 !js: from K to degC 1327 END DO 1328 1329 !----------------------------------------------------------------- 1330 ! brine salinity and liquid fraction 1331 !----------------------------------------------------------------- 1332 1333 DO k = 1, nlay_i 1334 1335 ! Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 1336 ! Best expression to date is that one 1337 Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 1338 phi(k) = salin(k) / Sbr 1339 1340 END DO 1341 1342 !----------------------------------------------------------------- 1343 ! permeability 1344 !----------------------------------------------------------------- 1345 1346 perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 1347 1348 END SUBROUTINE ice_thd_pnd_perm 1349 1350 1351 !---------------------------------------------------------------------------------------------------------------------- 304 1352 305 1353 SUBROUTINE ice_thd_pnd_init … … 317 1365 INTEGER :: ios, ioptio ! Local integer 318 1366 !! 319 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 320 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1367 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 1368 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1369 & ln_pnd_TOPO , & 321 1370 & ln_pnd_lids, ln_pnd_alb 322 1371 !!------------------------------------------------------------------- … … 336 1385 WRITE(numout,*) ' Namelist namicethd_pnd:' 337 1386 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 1387 WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO 338 1388 WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV 339 1389 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 340 1390 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 1391 WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush 341 1392 WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST 342 1393 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd … … 351 1402 IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF 352 1403 IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF 1404 IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF 353 1405 IF( ioptio /= 1 ) & 354 & 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)' )1406 & 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)' ) 355 1407 ! 356 1408 SELECT CASE( nice_pnd )
Note: See TracChangeset
for help on using the changeset viewer.
