- Timestamp:
- 2020-04-01T09:44:12+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/UKMO/NEMO_4.0.1_meltpond_improvements/src/ICE/icethd_pnd.F90
r11715 r12636 38 38 INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) 39 39 40 REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp 41 40 42 !! * Substitutions 41 43 # include "vectopt_loop_substitute.h90" … … 89 91 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 90 92 a_ip_frac_1d(ji) = rn_apnd 93 a_ip_eff_1d(ji) = rn_apnd 91 94 h_ip_1d(ji) = rn_hpnd 92 95 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 93 96 ELSE 94 97 a_ip_frac_1d(ji) = 0._wp 98 a_ip_eff_1d(ji) = 0._wp 95 99 h_ip_1d(ji) = 0._wp 96 100 a_ip_1d(ji) = 0._wp … … 129 133 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 130 134 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 131 ! 135 REAL(wp), PARAMETER :: max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1) 136 REAL(wp), PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 137 REAL(wp), PARAMETER :: pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 138 ! 139 REAL(wp) :: tot_mlt ! Total ice and snow surface melt (some goes into ponds, some into the ocean) 132 140 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 133 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 141 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding (equivalent volume change) 134 142 REAL(wp) :: z1_Tp ! inverse reference temperature 135 143 REAL(wp) :: z1_rhow ! inverse freshwater density 136 144 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 145 REAL(wp) :: z1_rhoi ! inverse ice density 137 146 REAL(wp) :: zfac, zdum 138 ! 139 INTEGER :: ji ! loop indices 147 REAL(wp) :: t_grad ! Temperature deficit for refreezing 148 REAL(wp) :: omega_dt ! Time independent accumulated variables used for freezing 149 REAL(wp) :: lh_ip_end ! Lid thickness at end of timestep (temporary variable) 150 REAL(wp) :: zdh_frz ! Amount of melt pond that freezes (m) 151 REAL(wp) :: v_ip_old ! Pond volume before leaking back to the ocean 152 REAL(wp) :: dh_ip_over ! Pond thickness change due to overflow or leaking 153 REAL(wp) :: dh_i_pndleak ! Grid box mean change in water depth due to leaking back to the ocean 154 REAL(wp) :: h_gravity_head ! Height above sea level of the top of the melt pond 155 REAL(wp) :: h_percolation ! Distance between the base of the melt pond and the base of the sea ice 156 REAL(wp) :: Sbr ! Brine salinity 157 REAL(wp), DIMENSION(nlay_i) :: phi ! liquid fraction 158 REAL(wp) :: perm ! Permeability of the sea ice 159 REAL(wp) :: za_ip ! Temporary pond fraction 160 REAL(wp) :: max_h_diff_ts ! Maximum meltpond depth change due to leaking or overflow (m per ts) 161 REAL(wp) :: lfrac_pnd ! The fraction of the meltpond exposed (not inder a frozen lid) 162 163 ! 164 INTEGER :: ji, jk ! loop indices 140 165 !!------------------------------------------------------------------- 141 166 z1_rhow = 1._wp / rhow 142 167 z1_zpnd_aspect = 1._wp / zpnd_aspect 143 168 z1_Tp = 1._wp / zTp 169 z1_rhoi = 1._wp / rhoi 170 max_h_diff_ts = max_h_diff_s * rdt_ice 171 172 ! Define time-independent field for use in refreezing 173 omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 144 174 145 175 DO ji = 1, npti 146 ! !--------------------------------! 147 IF( h_i_1d(ji) < rn_himin) THEN ! Case ice thickness < rn_himin ! 148 ! !--------------------------------! 149 !--- Remove ponds on thin ice 176 177 ! !----------------------------------------------------! 178 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 179 ! !----------------------------------------------------! 180 !--- Remove ponds on thin ice or tiny ice fractions 150 181 a_ip_1d(ji) = 0._wp 151 182 a_ip_frac_1d(ji) = 0._wp 152 183 h_ip_1d(ji) = 0._wp 184 lh_ip_1d(ji) = 0._wp 153 185 ! !--------------------------------! 154 186 ELSE ! Case ice thickness >= rn_himin ! 155 187 ! !--------------------------------! 156 188 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step 157 ! 158 ! available meltwater for melt ponding [m, >0] and fraction 159 zdv_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 160 zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) ! from CICE doc 161 !zfr_mlt = zrmin + zrmax * a_i_1d(ji) ! from Holland paper 189 190 ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 191 za_ip = a_ip_1d(ji) 192 IF ( za_ip < epsi06 ) za_ip = epsi06 193 ! 194 ! available meltwater for melt ponding [m, >0] and fraction 195 tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow 196 IF ( ln_pnd_totfrac ) THEN 197 zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction 198 ELSE 199 zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji) ! Use category ice fraction 200 ENDIF 201 zdv_mlt = zfr_mlt * tot_mlt 162 202 ! 163 203 !--- Pond gowth ---! 164 ! v_ip should never be negative, otherwise code crashes 165 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) 204 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 205 ! 206 !--- Lid shrinking. ---! 207 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip 166 208 ! 167 209 ! melt pond mass flux (<0) 168 IF( zdv_mlt > 0._wp ) THEN169 zfac = z fr_mlt * zdv_mlt * rhow * r1_rdtice210 IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN 211 zfac = zdv_mlt * rhow * r1_rdtice 170 212 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 171 213 ! … … 177 219 ! 178 220 !--- Pond contraction (due to refreezing) ---! 179 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 221 IF ( ln_pnd_lids ) THEN 222 223 ! Code to use if using melt pond lids 224 IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN 225 t_grad = (zTp+rt0) - t_su_1d(ji) 226 227 ! The following equation is a rearranged form of: 228 ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) 229 ! where: lid_thickness_start = lh_ip_1d(ji) 230 ! lid_thickness_end = lh_ip_end 231 ! omega_dt is a bunch of terms in the equation that do not change 232 ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier 233 ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density). 234 235 lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) 236 zdh_frz = lh_ip_end - lh_ip_1d(ji) 237 238 ! Pond shrinking 239 v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) 240 241 ! Lid growing 242 IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz 243 ELSE 244 zdh_frz = 0._wp 245 END IF 246 247 ELSE 248 ! Code to use if not using melt pond lids 249 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) 250 ENDIF 251 252 ! 253 ! Make sure pond volume or lid thickness has not gone negative 254 IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp 255 IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp 180 256 ! 181 257 ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac … … 184 260 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 185 261 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 262 263 !--- Pond overflow and vertical flushing ---! 264 IF ( ln_pnd_overflow ) THEN 265 266 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 267 IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN 268 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 269 dh_ip_over = zpnd_aspect * zfr_mlt - h_ip_1d(ji) ! This will be a negative number 270 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 271 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 272 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 273 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 274 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 275 ENDIF 276 277 ! If pond depth exceeds half the ice thickness then reduce the pond volume 278 IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN 279 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 280 dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji) ! This will be a negative number 281 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 282 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 283 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 284 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 285 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 286 ENDIF 287 288 !-- Vertical flushing of pond water --! 289 ! 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. 290 ! This assumes the pond is sitting on top of the ice. 291 h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow 292 293 ! The depth through which water percolates is the distance under the melt pond 294 h_percolation = h_i_1d(ji) - h_ip_1d(ji) 295 296 ! Calculate the permeability of the ice (Assur 1958) 297 DO jk = 1, nlay_i 298 Sbr = - 1.2_wp & 299 - 21.8_wp * (t_i_1d(ji,jk)-rt0) & 300 - 0.919_wp * (t_i_1d(ji,jk)-rt0)**2 & 301 - 0.01878_wp * (t_i_1d(ji,jk)-rt0)**3 302 phi(jk) = sz_i_1d(ji,jk)/Sbr 303 END DO 304 perm = 3.0e-08_wp * (minval(phi))**3 305 306 ! Do the drainage using Darcy's law 307 IF ( perm > 0._wp ) THEN 308 dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number 309 dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative 310 311 v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use 312 dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit 313 h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) 314 a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect 315 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 316 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 317 ENDIF 318 ENDIF 319 320 ! If lid thickness is ten times greater than pond thickness then remove pond 321 IF ( ln_pnd_lids ) THEN 322 IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 323 a_ip_1d(ji) = 0._wp 324 a_ip_frac_1d(ji) = 0._wp 325 h_ip_1d(ji) = 0._wp 326 lh_ip_1d(ji) = 0._wp 327 v_ip_1d(ji) = 0._wp 328 ENDIF 329 ENDIF 330 331 ! If any of the previous changes has removed all the ice thickness then remove ice area. 332 IF ( h_i_1d(ji) == 0._wp ) THEN 333 a_i_1d(ji) = 0._wp 334 h_s_1d(ji) = 0._wp 335 ENDIF 336 186 337 ! 187 338 ENDIF 339 340 ! Calculate how much meltpond is exposed to the atmosphere (not under a frozen lid) 341 IF ( lh_ip_1d(ji) < pnd_lid_min ) THEN ! Pond lid is very thin. Expose all the pond. 342 lfrac_pnd = 1.0 343 ELSE 344 IF ( lh_ip_1d(ji) > pnd_lid_max ) THEN ! Pond lid is very thick. Cover all the pond up with ice and nosw. 345 lfrac_pnd = 0.0 346 ELSE ! Pond lid is of intermediate thickness. Expose part of the melt pond. 347 lfrac_pnd = ( lh_ip_1d(ji) - pnd_lid_min ) / (pnd_lid_max - pnd_lid_min) 348 ENDIF 349 ENDIF 350 351 ! Calculate the melt pond effective area 352 a_ip_eff_1d(ji) = a_ip_frac_1d(ji) * lfrac_pnd 353 188 354 END DO 189 355 ! … … 205 371 INTEGER :: ios, ioptio ! Local integer 206 372 !! 207 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 373 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb, & 374 rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac, & 375 ln_use_pndmass 208 376 !!------------------------------------------------------------------- 209 377 ! … … 223 391 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 224 392 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 393 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_pnd_min = ', rn_pnd_min 394 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_pnd_max = ', rn_pnd_max 395 WRITE(numout,*) ' Use total ice fraction instead of category ice fraction ln_pnd_totfrac = ',ln_pnd_totfrac 396 WRITE(numout,*) ' Allow ponds to overflow and have vertical flushing ln_pnd_overflow = ',ln_pnd_overflow 397 WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ',ln_pnd_lids 398 WRITE(numout,*) ' Use melt pond mass flux diagnostic, passing to ocean ln_use_pndmass = ',ln_use_pndmass 225 399 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 226 400 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
Note: See TracChangeset
for help on using the changeset viewer.
