Changeset 14037 for NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icethd_pnd.F90
- Timestamp:
- 2020-12-03T12:20:38+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd_pnd.F90 (modified) (7 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG
- Property svn:externals
-
old new 8 8 9 9 # SETTE 10 ^/utils/CI/sette @13292sette10 ^/utils/CI/sette_wave@13990 sette
-
- Property svn:externals
-
NEMO/branches/2020/dev_r13333_KERNEL-08_techene_gm_HPG_SPG/src/ICE/icethd_pnd.F90
r12489 r14037 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 pond scheme 38 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 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 !! ** Purpose : change melt pond fraction 52 !! 53 !! ** Method : brut force 78 !! ** Purpose : change melt pond fraction and thickness 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 54 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) ) 55 89 ! 56 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 ) 57 99 ! 58 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 59 ! 60 CASE (np_pndH12) ; CALL pnd_H12 !== Holland et al 2012 melt ponds ==! 61 ! 62 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 63 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 64 151 END SUBROUTINE ice_thd_pnd 65 152 … … 81 168 !! ** References : Bush, G.W., and Trump, D.J. (2017) 82 169 !!------------------------------------------------------------------- 83 INTEGER :: ji ! loop indices 170 INTEGER :: ji, jl ! loop indices 171 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 84 172 !!------------------------------------------------------------------- 85 DO ji = 1, npti 86 ! 87 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 88 a_ip_frac_1d(ji) = rn_apnd 89 h_ip_1d(ji) = rn_hpnd 90 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 91 ELSE 92 a_ip_frac_1d(ji) = 0._wp 93 h_ip_1d(ji) = 0._wp 94 a_ip_1d(ji) = 0._wp 95 ENDIF 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 211 END DO 212 ! 213 END SUBROUTINE pnd_CST 214 215 216 SUBROUTINE pnd_LEV 217 !!------------------------------------------------------------------- 218 !! *** ROUTINE pnd_LEV *** 219 !! 220 !! ** Purpose : Compute melt pond evolution 221 !! 222 !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing 223 !! We work with volumes and then redistribute changes into thickness and concentration 224 !! assuming linear relationship between the two. 225 !! 226 !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- 227 !! dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 228 !! dh_i = meltwater from ice surface melt 229 !! dh_s = meltwater from snow melt 230 !! (1-r) = fraction of melt water that is not flushed 231 !! 232 !! - limtations: a_ip must not exceed (1-r)*a_i 233 !! h_ip must not exceed 0.5*h_i 234 !! 235 !! - pond shrinking: 236 !! if lids: Vp = Vp -dH * a_ip 237 !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- 238 !! 239 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 240 !! H = lid thickness 241 !! Lf = latent heat of fusion 242 !! Tp = -2C 243 !! 244 !! And solved implicitely as: 245 !! H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0 246 !! 247 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 248 !! 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) 251 !! visc = water viscosity 252 !! Hp = height of top of the pond above sea-level 253 !! Hi = ice thickness thru which there is flushing 254 !! flush= correction otherwise flushing is excessive 255 !! 256 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness 257 !! 258 !! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip: 259 !! a_ip/a_i = a_ip_frac = h_ip / zaspect 260 !! 261 !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 262 !! 263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 264 !! 265 !! ** References : Flocco and Feltham (JGR, 2007) 266 !! Flocco et al (JGR, 2010) 267 !! Holland et al (J. Clim, 2012) 268 !! Hunke et al (OM 2012) 269 !!------------------------------------------------------------------- 270 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 271 !! 272 REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio 273 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 274 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 275 !! 276 REAL(wp) :: zfr_mlt, zdv_mlt, zdv_avail ! fraction and volume of available meltwater retained for melt ponding 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 279 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 280 REAL(wp) :: zv_ip_max ! max pond volume allowed 281 REAL(wp) :: zdT ! zTp-t_su 282 REAL(wp) :: zsbr, ztmelts ! Brine salinity 283 REAL(wp) :: zperm ! permeability of sea ice 284 REAL(wp) :: zfac, zdum ! temporary arrays 285 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 286 !! 287 INTEGER :: ji, jk, jl ! loop indices 288 !!------------------------------------------------------------------- 289 z1_rhow = 1._wp / rhow 290 z1_aspect = 1._wp / zaspect 291 z1_Tp = 1._wp / zTp 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) 452 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 453 ! 454 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 455 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 456 ! 457 END DO 458 459 !-------------------------------------------------------------------- 460 ! Retrieve 2D arrays 461 !-------------------------------------------------------------------- 462 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) ) 466 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) ) 96 467 ! 97 468 END DO 98 469 ! 99 END SUBROUTINE pnd_CST 100 101 102 SUBROUTINE pnd_H12 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 ! 477 END SUBROUTINE pnd_LEV 478 479 480 481 SUBROUTINE pnd_TOPO 482 103 483 !!------------------------------------------------------------------- 104 !! *** ROUTINE pnd_H12 *** 105 !! 106 !! ** Purpose : Compute melt pond evolution 107 !! 108 !! ** Method : Empirical method. A fraction of meltwater is accumulated in ponds 109 !! and sent to ocean when surface is freezing 110 !! 111 !! pond growth: Vp = Vp + dVmelt 112 !! with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 113 !! pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) 114 !! with Tp = -2degC 115 !! 116 !! ** Tunable parameters : (no real expertise yet, ideas?) 117 !! 118 !! ** Note : Stolen from CICE for quick test of the melt pond 119 !! radiation and freshwater interfaces 120 !! Coupling can be radiative AND freshwater 121 !! Advection, ridging, rafting are called 122 !! 123 !! ** References : Holland, M. M. et al (J Clim 2012) 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 !! 124 511 !!------------------------------------------------------------------- 125 REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding 126 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum - - - - - 127 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 128 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 129 ! 130 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 131 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 132 REAL(wp) :: z1_Tp ! inverse reference temperature 133 REAL(wp) :: z1_rhow ! inverse freshwater density 134 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 135 REAL(wp) :: zfac, zdum 136 ! 137 INTEGER :: ji ! loop indices 138 !!------------------------------------------------------------------- 139 z1_rhow = 1._wp / rhow 140 z1_zpnd_aspect = 1._wp / zpnd_aspect 141 z1_Tp = 1._wp / zTp 142 143 DO ji = 1, npti 144 ! !--------------------------------! 145 IF( h_i_1d(ji) < rn_himin) THEN ! Case ice thickness < rn_himin ! 146 ! !--------------------------------! 147 !--- Remove ponds on thin ice 148 a_ip_1d(ji) = 0._wp 149 a_ip_frac_1d(ji) = 0._wp 150 h_ip_1d(ji) = 0._wp 151 ! !--------------------------------! 152 ELSE ! Case ice thickness >= rn_himin ! 153 ! !--------------------------------! 154 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step 155 ! 156 ! available meltwater for melt ponding [m, >0] and fraction 157 zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 158 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc 159 !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper 160 ! 161 !--- Pond gowth ---! 162 ! v_ip should never be negative, otherwise code crashes 163 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 164 ! 165 ! melt pond mass flux (<0) 166 IF( zdv_mlt > 0._wp ) THEN 167 zfac = zfr_mlt * zdv_mlt * rhow * r1_Dt_ice 168 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 169 ! 170 ! adjust ice/snow melting flux to balance melt pond flux (>0) 171 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 172 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 173 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 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 174 749 ENDIF 175 ! 176 !--- Pond contraction (due to refreezing) ---! 177 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 178 ! 179 ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 180 ! h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i 181 a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) ) 182 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 183 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 184 ! 185 ENDIF 186 END DO 187 ! 188 END SUBROUTINE pnd_H12 189 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 190 1336 191 1337 SUBROUTINE ice_thd_pnd_init … … 203 1349 INTEGER :: ios, ioptio ! Local integer 204 1350 !! 205 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 1351 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 1352 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1353 & ln_pnd_TOPO, & 1354 & ln_pnd_lids, ln_pnd_alb 206 1355 !!------------------------------------------------------------------- 207 1356 ! … … 217 1366 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 218 1367 WRITE(numout,*) ' Namelist namicethd_pnd:' 219 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 220 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 221 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 222 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 223 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 224 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 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 1370 WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV 1371 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 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 1374 WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST 1375 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 1376 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 1377 WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids 1378 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 225 1379 ENDIF 226 1380 ! … … 229 1383 IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF 230 1384 IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF 231 IF( ln_pnd_H12 ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndH12 ; ENDIF 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 232 1387 IF( ioptio /= 1 ) & 233 & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_ H12 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)' ) 234 1389 ! 235 1390 SELECT CASE( nice_pnd ) 236 1391 CASE( np_pndNO ) 237 IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF 1392 IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF 1393 IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF 1394 CASE( np_pndCST ) 1395 IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF 238 1396 END SELECT 239 1397 !
Note: See TracChangeset
for help on using the changeset viewer.
