- Timestamp:
- 2020-03-02T17:57:51+01:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_pnd.F90
r12497 r12500 38 38 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 39 39 40 REAL(wp), PUBLIC, PARAMETER :: a_pnd_avail = 0.7_wp ! Fraction of sea ice available for melt ponding41 40 REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 42 41 REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation … … 137 136 REAL(wp), PARAMETER :: max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1) 138 137 ! 138 REAL(wp) :: tot_mlt ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 139 139 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 140 140 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding (equivalent volume change) … … 151 151 REAL(wp) :: dh_ip_over ! Pond thickness change due to overflow or leaking 152 152 REAL(wp) :: dh_i_pndleak ! Grid box mean change in water depth due to leaking back to the ocean 153 REAL(wp) :: weighted_lid_snow ! Lid to go on ponds under snow if snow thickness exceeds snow_lid_min154 153 REAL(wp) :: h_gravity_head ! Height above sea level of the top of the melt pond 155 154 REAL(wp) :: h_percolation ! Distance between the base of the melt pond and the base of the sea ice … … 182 181 h_ip_1d(ji) = 0._wp 183 182 lh_ip_1d(ji) = 0._wp 184 185 zdv_mlt = 0._wp186 zdh_frz = 0._wp187 183 ! !--------------------------------! 188 184 ELSE ! Case ice thickness >= rn_himin ! … … 194 190 IF ( za_ip < epsi06 ) za_ip = epsi06 195 191 ! 196 ! available meltwater for melt ponding 197 ! This is the change in ice volume due to melt 198 zdv_mlt = -(dh_i_pnd(ji)*rhoi + dh_s_pnd(ji)*rhos) * z1_rhow * a_i_1d(ji) 192 ! available meltwater for melt ponding [m, >0] and fraction 193 tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 194 IF ( ln_pnd_totfrac ) THEN 195 zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction 196 ELSE 197 zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji) ! Use category ice fraction 198 ENDIF 199 zdv_mlt = zfr_mlt * tot_mlt 199 200 ! 200 201 !--- Pond gowth ---! … … 202 203 ! 203 204 !--- Lid shrinking. ---! 204 lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 205 ! 205 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 206 ! 207 ! melt pond mass flux (<0) 208 IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN 209 zfac = zdv_mlt * rhow * r1_rdtice 210 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 211 ! 212 ! adjust ice/snow melting flux to balance melt pond flux (>0) 213 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 214 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 215 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 216 ENDIF 206 217 ! 207 218 !--- Pond contraction (due to refreezing) ---! 208 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 209 t_grad = (zTp+rt0) - t_su_1d(ji) 210 211 ! The following equation is a rearranged form of: 212 ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 213 ! where: lid_thickness_start = lh_ip_1d(ji) 214 ! lid_thickness_end = lh_ip_end 215 ! omega_dt is a bunch of terms in the equation that do not change 216 ! note the use of rhow instead of rhoi as we are working with volumes and it is easier if the water and ice volumes (for the lid and the pond) are the same 217 ! (have the same density). The lid will eventually remelt anyway so it doesn't matter if we make this simplification. 218 219 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 220 zdh_frz = lh_ip_end - lh_ip_1d(ji) 221 222 ! Pond shrinking 223 v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 224 225 ! Lid growing 226 lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 219 IF ( ln_pnd_lids ) THEN 220 221 ! Code to use if using melt pond lids 222 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 223 t_grad = (zTp+rt0) - t_su_1d(ji) 224 225 ! The following equation is a rearranged form of: 226 ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 227 ! where: lid_thickness_start = lh_ip_1d(ji) 228 ! lid_thickness_end = lh_ip_end 229 ! omega_dt is a bunch of terms in the equation that do not change 230 ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 231 ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density). 232 233 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 234 zdh_frz = lh_ip_end - lh_ip_1d(ji) 235 236 ! Pond shrinking 237 v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 238 239 ! Lid growing 240 F ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 241 ELSE 242 zdh_frz = 0._wp 243 END IF 244 227 245 ELSE 228 zdh_frz = 0._wp 229 END IF 246 ! Code to use if not using melt pond lids 247 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 248 ENDIF 230 249 231 250 ! … … 239 258 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 240 259 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 241 242 ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume 243 IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN 244 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 245 dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji) ! This will be a negative number 246 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 247 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 248 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 249 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 250 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 251 260 261 !--- Pond overflow and vertical flushing ---! 262 IF ( ln_pnd_overflow ) THEN 263 264 ! If pond area exceeds a_pnd_avail_1d(ji) * a_i_1d(ji) then reduce the pond volume 265 IF ( a_ip_1d(ji) > a_pnd_avail_1d(ji) * a_i_1d(ji) ) THEN 266 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 267 dh_ip_over = zpnd_aspect * a_pnd_avail_1d(ji) - h_ip_1d(ji) ! This will be a negative number 268 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 269 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 270 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 271 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 272 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 273 ENDIF 274 275 ! If pond depth exceeds half the ice thickness then reduce the pond volume 276 IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 277 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 278 dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji) ! This will be a negative number 279 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 280 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 281 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 282 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 283 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 284 ENDIF 285 286 !-- Vertical flushing of pond water --! 287 ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth. 288 ! This assumes the pond is sitting on top of the ice. 289 h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 290 291 ! The depth through which water percolates is the distance under the melt pond 292 h_percolation = h_i_1d(ji) - h_ip_1d(ji) 293 294 ! Calculate the permeability of the ice (Assur 1958) 295 DO jk = 1, nlay_i 296 Sbr = - 1.2_wp & 297 - 21.8_wp * (t_i_1d(ji,jk)-rt0) & 298 - 0.919_wp * (t_i_1d(ji,jk)-rt0)**2 & 299 - 0.01878_wp * (t_i_1d(ji,jk)-rt0)**3 300 phi(jk) = sz_i_1d(ji,jk)/Sbr 301 END DO 302 perm = 3.0e-08_wp * (minval(phi))**3 303 304 ! Do the drainage using Darcy's law 305 IF ( perm > 0._wp ) THEN 306 dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number 307 dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative 308 309 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 310 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 311 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 312 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 313 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 314 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 315 ENDIF 252 316 ENDIF 253 317 254 ! If pond depth exceeds half the ice thickness then reduce the pond volume 255 IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 256 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 257 dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji) ! This will be a negative number 258 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 259 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 260 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 261 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 262 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 263 318 ! If lid thickness is ten times greater than pond thickness then remove pond 319 IF ( ln_pnd_lids ) THEN 320 IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 321 a_ip_1d(ji) = 0._wp 322 a_ip_frac_1d(ji) = 0._wp 323 h_ip_1d(ji) = 0._wp 324 lh_ip_1d(ji) = 0._wp 325 v_ip_1d(ji) = 0._wp 326 ENDIF 264 327 ENDIF 265 266 !-- Vertical flushing of pond water --!267 ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth.268 ! This assumes the pond is sitting on top of the ice.269 h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow270 271 ! The depth through which water percolates is the distance under the melt pond272 h_percolation = h_i_1d(ji) - h_ip_1d(ji)273 274 ! Calculate the permeability of the ice (Assur 1958)275 DO jk = 1, nlay_i276 Sbr = - 1.2_wp &277 - 21.8_wp * (t_i_1d(ji,jk)-rt0) &278 - 0.919_wp * (t_i_1d(ji,jk)-rt0)**2 &279 - 0.01878_wp * (t_i_1d(ji,jk)-rt0)**3280 phi(jk) = sz_i_1d(ji,jk)/Sbr281 END DO282 perm = 3.0e-08_wp * (minval(phi))**3283 284 ! Do the drainage using Darcy's law285 IF ( perm > 0._wp ) THEN286 dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number287 dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative288 289 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use290 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit291 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over)292 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect293 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji)294 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji)295 296 ENDIF297 298 ! If lid thickness is ten times greater than pond thickness then remove pond299 IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN300 a_ip_1d(ji) = 0._wp301 a_ip_frac_1d(ji) = 0._wp302 h_ip_1d(ji) = 0._wp303 lh_ip_1d(ji) = 0._wp304 v_ip_1d(ji) = 0._wp305 END IF306 328 307 329 ! If any of the previous changes has removed all the ice thickness then remove ice area. … … 333 355 INTEGER :: ios, ioptio ! Local integer 334 356 !! 335 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 357 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb, & 358 rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac, & 359 ln_use_pndmass 336 360 !!------------------------------------------------------------------- 337 361 ! … … 351 375 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 352 376 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 377 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_pnd_min = ', rn_pnd_min 378 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_pnd_max = ', rn_pnd_max 379 WRITE(numout,*) ' Use total ice fraction instead of category ice fraction ln_pnd_totfrac = ',ln_pnd_totfrac 380 WRITE(numout,*) ' Allow ponds to overflow and have vertical flushing ln_pnd_overflow = ',ln_pnd_overflow 381 WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ',ln_pnd_lids 382 WRITE(numout,*) ' Use melt pond mass flux diagnostic, passing to ocean ln_use_pndmass = ',ln_use_pndmass 353 383 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 354 384 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
Note: See TracChangeset
for help on using the changeset viewer.
