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