- Timestamp:
- 2020-04-08T18:54:44+02:00 (4 years ago)
- File:
-
- 1 edited
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/r4.0-HEAD_r12713_clem_dan_fixcpl/src/ICE/icethd_pnd.F90
r11536 r12720 89 89 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 90 90 a_ip_frac_1d(ji) = rn_apnd 91 a_ip_eff_1d(ji) = rn_apnd 91 92 h_ip_1d(ji) = rn_hpnd 92 93 a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) 94 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 93 95 ELSE 94 96 a_ip_frac_1d(ji) = 0._wp 97 a_ip_eff_1d(ji) = 0._wp 95 98 h_ip_1d(ji) = 0._wp 96 99 a_ip_1d(ji) = 0._wp 100 h_il_1d(ji) = 0._wp 97 101 ENDIF 98 102 ! … … 106 110 !! *** ROUTINE pnd_H12 *** 107 111 !! 108 !! ** Purpose : Compute melt pond evolution 109 !! 110 !! ** Method : Empirical method. A fraction of meltwater is accumulated in ponds 111 !! and sent to ocean when surface is freezing 112 !! 113 !! pond growth: Vp = Vp + dVmelt 114 !! with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 115 !! pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) 116 !! with Tp = -2degC 117 !! 118 !! ** Tunable parameters : (no real expertise yet, ideas?) 112 !! ** Purpose : Compute melt pond evolution 113 !! 114 !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing 115 !! We work with volumes and then redistribute changes into thickness and concentration 116 !! assuming linear relationship between the two. 117 !! 118 !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- 119 !! dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i 120 !! dh_i = meltwater from ice surface melt 121 !! dh_s = meltwater from snow melt 122 !! (1-r) = fraction of melt water that is not flushed 123 !! 124 !! - limtations: a_ip must not exceed (1-r)*a_i 125 !! h_ip must not exceed 0.5*h_i 126 !! 127 !! - pond shrinking: 128 !! if lids: Vp = Vp -dH * a_ip 129 !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- 130 !! 131 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 132 !! H = lid thickness 133 !! Lf = latent heat of fusion 134 !! Tp = -2C 135 !! 136 !! And solved implicitely as: 137 !! H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0 138 !! 139 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 140 !! 141 !! - Overflow: w = -perm/visc * rho_oce * grav * Hp / Hi --- from Flocco et al 2007 --- 142 !! perm = permability of sea-ice 143 !! visc = water viscosity 144 !! Hp = height of top of the pond above sea-level 145 !! Hi = ice thickness thru which there is flushing 146 !! 147 !! 148 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness 149 !! 150 !! - effective pond area: to be used for albedo 151 !! 152 !! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip: 153 !! a_ip/a_i = a_ip_frac = h_ip / zaspect 154 !! 155 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 119 156 !! 120 !! ** Note : Stolen from CICE for quick test of the melt pond 121 !! radiation and freshwater interfaces 122 !! Coupling can be radiative AND freshwater 123 !! Advection, ridging, rafting are called 124 !! 125 !! ** References : Holland, M. M. et al (J Clim 2012) 126 !!------------------------------------------------------------------- 127 REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding 128 REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum - - - - - 129 REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio 130 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 131 ! 132 REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding 133 REAL(wp) :: zdv_mlt ! available meltwater for melt ponding 134 REAL(wp) :: z1_Tp ! inverse reference temperature 135 REAL(wp) :: z1_rhow ! inverse freshwater density 136 REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio 137 REAL(wp) :: zfac, zdum 138 ! 139 INTEGER :: ji ! loop indices 140 !!------------------------------------------------------------------- 141 z1_rhow = 1._wp / rhow 142 z1_zpnd_aspect = 1._wp / zpnd_aspect 143 z1_Tp = 1._wp / zTp 157 !! ** Note : mostly stolen from CICE 158 !! 159 !! ** References : Flocco and Feltham (JGR, 2007) 160 !! Flocco et al (JGR, 2010) 161 !! Holland et al (J. Clim, 2012) 162 !!------------------------------------------------------------------- 163 REAL(wp), DIMENSION(nlay_i) :: zperm ! Permeability of sea ice 164 !! 165 REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio 166 REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature 167 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 168 REAL(wp), PARAMETER :: zhl_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation 169 REAL(wp), PARAMETER :: zhl_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation 170 !! 171 REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding 172 REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes 173 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 174 REAL(wp) :: v_ip_max ! max pond volume allowed 175 REAL(wp) :: zdT ! zTp-t_su 176 REAL(wp) :: zsbr ! Brine salinity 177 REAL(wp) :: zfac, zdum ! temporary arrays 178 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 179 !! 180 INTEGER :: ji, jk ! loop indices 181 !!------------------------------------------------------------------- 182 z1_rhow = 1._wp / rhow 183 z1_aspect = 1._wp / zaspect 184 z1_Tp = 1._wp / zTp 144 185 145 186 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 187 ! !----------------------------------------------------! 188 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 189 ! !----------------------------------------------------! 190 !--- Remove ponds on thin ice or tiny ice fractions 150 191 a_ip_1d(ji) = 0._wp 151 192 a_ip_frac_1d(ji) = 0._wp 152 193 h_ip_1d(ji) = 0._wp 153 ! !--------------------------------! 154 ELSE ! Case ice thickness >= rn_himin ! 155 ! !--------------------------------! 156 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 162 ! 163 !--- 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 ) 166 ! 167 ! melt pond mass flux (<0) 194 h_il_1d(ji) = 0._wp 195 ! 196 ! clem: problem with conservation or not ? 197 ! !--------------------------------! 198 ELSE ! Case ice thickness >= rn_himin ! 199 ! !--------------------------------! 200 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 201 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 202 ! 203 !------------------! 204 ! case ice melting ! 205 !------------------! 206 ! 207 !--- available meltwater for melt ponding ---! 208 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 209 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(ji) ! = ( 1 - r ) in H12 = fraction of melt water that is not flushed 210 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 211 ! 212 !--- overflow ---! 213 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 214 ! a_ip_max = zfr_mlt * a_i 215 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 216 v_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 217 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) ) 218 219 ! If pond depth exceeds half the ice thickness then reduce the pond volume 220 ! h_ip_max = 0.5 * h_i 221 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 222 v_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 223 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, v_ip_max - v_ip_1d(ji) ) ) 224 225 !--- Pond growing ---! 226 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 227 ! 228 !--- Lid melting ---! 229 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 230 ! 231 !--- mass flux ---! 168 232 IF( zdv_mlt > 0._wp ) THEN 169 zfac = z fr_mlt * zdv_mlt * rhow * r1_rdtice233 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 170 234 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 171 235 ! 172 ! adjust ice/snow melting flux to balance melt pond flux (>0) 173 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 236 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 174 237 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 175 238 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 176 239 ENDIF 240 241 !-------------------! 242 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 243 !-------------------! 244 ! 245 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 177 246 ! 178 247 !--- 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 ) 180 ! 181 ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac 182 ! h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i 183 a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) ) 248 IF( ln_pnd_lids ) THEN 249 ! 250 !--- Lid growing and subsequent pond shrinking ---! 251 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 252 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 253 254 ! Lid growing 255 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 256 257 ! Pond shrinking 258 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 259 260 ELSE 261 ! Pond shrinking 262 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 263 ENDIF 264 ! 265 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 266 ! v_ip = h_ip * a_ip 267 ! 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) 268 a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) 184 269 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 185 h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) 270 h_ip_1d(ji) = zaspect * a_ip_frac_1d(ji) 271 272 !---------------! 273 ! Pond flushing ! 274 !---------------! 275 IF( ln_pnd_flush ) THEN 276 ! height of top of the pond above sea-level 277 zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_frac_1d(ji) ) ) * r1_rau0 278 279 ! Calculate the permeability of the ice (Assur 1958) 280 DO jk = 1, nlay_i 281 zsbr = - 1.2_wp & 282 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 283 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 284 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 ! clem: error here the factor was 0.01878 instead of 0.0178 (cf Flocco 2010) 285 zperm(jk) = MAX( 0._wp, 3.e-08_wp * (sz_i_1d(ji,jk) / zsbr)**3 ) 286 END DO 287 288 ! Do the drainage using Darcy's law 289 zdv_flush = -MINVAL(zperm(:)) * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 290 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 291 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 292 293 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 294 a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_aspect * a_i_1d(ji) ) 295 a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) 296 h_ip_1d(ji) = zaspect * a_ip_frac_1d(ji) 297 298 ENDIF 299 300 !--- Corrections and lid thickness ---! 301 IF( ln_pnd_lids ) THEN 302 !--- remove ponds if lids are much larger than ponds ---! 303 IF ( v_il_1d(ji) > v_ip_1d(ji) * 10._wp ) THEN 304 a_ip_1d(ji) = 0._wp 305 a_ip_frac_1d(ji) = 0._wp 306 h_ip_1d(ji) = 0._wp 307 v_il_1d(ji) = 0._wp 308 ENDIF 309 !--- retrieve lid thickness from volume ---! 310 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 311 ELSE ; h_il_1d(ji) = 0._wp ; 312 ENDIF 313 ENDIF 186 314 ! 187 315 ENDIF 316 188 317 END DO 318 319 !-------------------------------------------------! 320 ! How much melt pond is exposed to the atmosphere ! 321 !-------------------------------------------------! 322 ! Calculate the melt pond effective area (used for albedo) 323 WHERE ( h_il_1d(1:npti) <= zhl_min ) ; a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti) ! lid is very thin. Expose all the pond 324 ELSEWHERE( h_il_1d(1:npti) >= zhl_max ) ; a_ip_eff_1d(1:npti) = 0._wp ! lid is very thick. Cover all the pond up with ice and snow 325 ELSEWHERE ; a_ip_eff_1d(1:npti) = a_ip_frac_1d(1:npti) * & ! lid is in between. Expose part of the pond 326 & ( h_il_1d(1:npti) - zhl_min ) / ( zhl_max - zhl_min ) 327 END WHERE 189 328 ! 190 329 END SUBROUTINE pnd_H12 … … 205 344 INTEGER :: ios, ioptio ! Local integer 206 345 !! 207 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb 346 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_H12, ln_pnd_lids, ln_pnd_flush, rn_apnd_min, rn_apnd_max, & 347 & ln_pnd_CST, rn_apnd, rn_hpnd, & 348 & ln_pnd_alb 208 349 !!------------------------------------------------------------------- 209 350 ! … … 221 362 WRITE(numout,*) '~~~~~~~~~~~~~~~~' 222 363 WRITE(numout,*) ' Namelist namicethd_pnd:' 223 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 224 WRITE(numout,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 225 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 226 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 227 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 228 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 364 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 365 WRITE(numout,*) ' Evolutive melt pond fraction and depth ln_pnd_H12 = ', ln_pnd_H12 366 WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ', ln_pnd_lids 367 WRITE(numout,*) ' Allow ponds to flush thru the ice ln_pnd_flush = ', ln_pnd_flush 368 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 369 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 370 WRITE(numout,*) ' Prescribed melt pond fraction and depth ln_pnd_CST = ', ln_pnd_CST 371 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd 372 WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd 373 WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb 229 374 ENDIF 230 375 !
Note: See TracChangeset
for help on using the changeset viewer.
