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