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