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