Changeset 13850
- Timestamp:
- 2020-11-23T14:10:46+01:00 (3 years ago)
- Location:
- NEMO/branches/2020/SI3-05_MeltPonds_topo
- Files:
-
- 6 edited
-
cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg (modified) (1 diff)
-
cfgs/SHARED/namelist_ice_ref (modified) (1 diff)
-
src/ICE/ice.F90 (modified) (3 diffs)
-
src/ICE/icestp.F90 (modified) (1 diff)
-
src/ICE/icethd.F90 (modified) (3 diffs)
-
src/ICE/icethd_pnd.F90 (modified) (11 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/ORCA2_SAS_ICE/EXPREF/namelist_cfg
r13278 r13850 23 23 cn_exp = "ORCA2_SAS" ! experience name 24 24 nn_it000 = 1 ! first time step 25 nn_itend = 100 ! last time step (std 5475)25 nn_itend = 5840 ! last time step (std 5840) 26 26 / 27 27 !----------------------------------------------------------------------- -
NEMO/branches/2020/SI3-05_MeltPonds_topo/cfgs/SHARED/namelist_ice_ref
r13346 r13850 195 195 !------------------------------------------------------------------------------ 196 196 ln_pnd = .true. ! activate melt ponds or not 197 ln_pnd_LEV = .true. ! level ice melt ponds (from Flocco et al 2007,2010 & Holland et al 2012) 197 ln_pnd_TOPO = .false. ! topographic melt ponds 198 ln_pnd_LEV = .true. ! level ice melt ponds 198 199 rn_apnd_min = 0.15 ! minimum ice fraction that contributes to melt pond. range: 0.0 -- 0.15 ?? 199 200 rn_apnd_max = 0.85 ! maximum ice fraction that contributes to melt pond. range: 0.7 -- 0.85 ?? -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/ice.F90
r13640 r13850 208 208 ! !!** ice-ponds namelist (namthd_pnd) 209 209 LOGICAL , PUBLIC :: ln_pnd !: Melt ponds (T) or not (F) 210 LOGICAL , PUBLIC :: ln_pnd_LEV !: Melt ponds scheme from Holland et al (2012), Flocco et al (2007, 2010) 210 LOGICAL , PUBLIC :: ln_pnd_TOPO !: Topographic Melt ponds scheme (Flocco et al 2007, 2010) 211 LOGICAL , PUBLIC :: ln_pnd_LEV !: Simple melt pond scheme 211 212 REAL(wp), PUBLIC :: rn_apnd_min !: Minimum ice fraction that contributes to melt ponds 212 213 REAL(wp), PUBLIC :: rn_apnd_max !: Maximum ice fraction that contributes to melt ponds … … 308 309 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: t1_ice !: temperature of the first layer (ln_cndflx=T) [K] 309 310 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: cnd_ice !: effective conductivity of the 1st layer (ln_cndflx=T) [W.m-2.K-1] 311 312 ! meltwater arrays to save for melt ponds 313 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_i_sum_2d !: surface melt (2d arrays for ponds) [m] 314 REAL(wp), PUBLIC, ALLOCATABLE, SAVE, DIMENSION(:,:,:) :: dh_s_mlt_2d !: snow surf melt (2d arrays for ponds) [m] 310 315 311 316 !!---------------------------------------------------------------------- … … 458 463 ii = ii + 1 459 464 ALLOCATE( qtr_ice_bot(jpi,jpj,jpl) , cnd_ice(jpi,jpj,jpl) , t1_ice(jpi,jpj,jpl) , & 465 & dh_i_sum_2d(jpi,jpj,jpl) , dh_s_mlt_2d(jpi,jpj,jpl) , & 460 466 & h_i (jpi,jpj,jpl) , a_i (jpi,jpj,jpl) , v_i (jpi,jpj,jpl) , & 461 467 & v_s (jpi,jpj,jpl) , h_s (jpi,jpj,jpl) , t_su (jpi,jpj,jpl) , & -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icestp.F90
r13720 r13850 414 414 wfx_pnd(ji,jj) = 0._wp 415 415 416 dh_i_sum_2d(:,:,:) = 0._wp 417 dh_s_mlt_2d(:,:,:) = 0._wp 418 416 419 hfx_thd(ji,jj) = 0._wp ; 417 420 hfx_snw(ji,jj) = 0._wp ; hfx_opw(ji,jj) = 0._wp -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd.F90
r13642 r13850 247 247 IF( ln_icedH ) THEN ! --- Growing/Melting --- ! 248 248 CALL ice_thd_dh ! Ice-Snow thickness 249 CALL ice_thd_pnd ! Melt ponds formation250 249 CALL ice_thd_ent( e_i_1d(1:npti,:) ) ! Ice enthalpy remapping 251 250 ENDIF … … 268 267 IF( ln_icediachk ) CALL ice_cons2D (1, 'icethd', diag_v, diag_s, diag_t, diag_fv, diag_fs, diag_ft) 269 268 ! 269 IF ( ln_pnd ) CALL ice_thd_pnd ! --- Melt ponds 270 270 271 IF( jpl > 1 ) CALL ice_itd_rem( kt ) ! --- Transport ice between thickness categories --- ! 271 272 ! … … 536 537 CALL tab_1d_2d( npti, nptidx(1:npti), cnd_ice_1d(1:npti), cnd_ice(:,:,kl) ) 537 538 CALL tab_1d_2d( npti, nptidx(1:npti), t1_ice_1d (1:npti), t1_ice (:,:,kl) ) 539 ! ponds 540 CALL tab_1d_2d( npti, nptidx(1:npti), dh_i_sum (1:npti) , dh_i_sum_2d(:,:,kl) ) 541 CALL tab_1d_2d( npti, nptidx(1:npti), dh_s_mlt (1:npti) , dh_s_mlt_2d(:,:,kl) ) 538 542 ! SIMIP diagnostics 539 543 CALL tab_1d_2d( npti, nptidx(1:npti), t_si_1d (1:npti), t_si (:,:,kl) ) -
NEMO/branches/2020/SI3-05_MeltPonds_topo/src/ICE/icethd_pnd.F90
r13284 r13850 1 1 MODULE icethd_pnd 2 2 !!====================================================================== 3 !! --- closest to CICE version 3 4 !! *** MODULE icethd_pnd *** 4 5 !! sea-ice: Melt ponds on top of sea ice … … 20 21 USE ice1D ! sea-ice: thermodynamics variables 21 22 USE icetab ! sea-ice: 1D <==> 2D transformation 23 USE sbc_ice ! surface energy budget 22 24 ! 23 25 USE in_out_manager ! I/O manager … … 34 36 INTEGER :: nice_pnd ! choice of the type of pond scheme 35 37 ! ! associated indices: 36 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 37 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 38 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 38 INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme 39 INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme 40 INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme 41 INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme 39 42 40 43 !! * Substitutions … … 52 55 !! 53 56 !! ** Purpose : change melt pond fraction and thickness 57 !! 58 !! Note: Melt ponds affect only radiative transfer for now 59 !! 60 !! No heat, no salt. 61 !! The melt water they carry is collected but 62 !! not removed from fw budget or released to the ocean 63 !! 64 !! A wfx_pnd has been coded for diagnostic purposes 65 !! It is not fully consistent yet. 66 !! 67 !! The current diagnostic lacks a contribution from drainage 68 !! 69 !! Each time wfx_pnd is updated, wfx_sum / wfx_snw_sum must be updated 54 70 !! 55 71 !!------------------------------------------------------------------- … … 57 73 SELECT CASE ( nice_pnd ) 58 74 ! 59 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==!75 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 60 76 ! 61 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 77 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 78 ! 79 CASE (np_pndTOPO) ; CALL pnd_TOPO !& !== Topographic melt ponds ==! 62 80 ! 63 81 END SELECT … … 75 93 !! to non-zero values when t_su = 0C 76 94 !! 77 !! ** Tunable parameters : pond fraction (rn_apnd), ponddepth (rn_hpnd)95 !! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd) 78 96 !! 79 97 !! ** Note : Coupling with such melt ponds is only radiative … … 147 165 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min 148 166 !! 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) 167 !! ** Note : mostly stolen from CICE but not only 168 !! 169 !! ** References : Holland et al (J. Clim, 2012) 170 !! 154 171 !!------------------------------------------------------------------- 172 155 173 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array 156 174 !! … … 169 187 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 170 188 !! 171 INTEGER :: ji, jk ! loop indices 189 INTEGER :: ji, jj, jk, jl ! loop indices 190 172 191 !!------------------------------------------------------------------- 192 173 193 z1_rhow = 1._wp / rhow 174 194 z1_aspect = 1._wp / zaspect 175 195 z1_Tp = 1._wp / zTp 176 177 DO ji = 1, npti 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 182 a_ip_1d(ji) = 0._wp 183 h_ip_1d(ji) = 0._wp 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) ) ) 196 197 !----------------------------------------------------------------------------------------------- 198 ! Identify grid cells with ice 199 !----------------------------------------------------------------------------------------------- 200 at_i(:,:) = SUM( a_i, dim=3 ) 201 ! 202 npti = 0 ; nptidx(:) = 0 203 DO jj = 1, jpj 204 DO ji = 1, jpi 205 IF ( at_i(ji,jj) > epsi10 ) THEN 206 npti = npti + 1 207 nptidx( npti ) = (jj - 1) * jpi + ji 208 ENDIF 209 END DO 210 END DO 211 212 !----------------------------------------------------------------------------------------------- 213 ! Prepare 1D arrays 214 !----------------------------------------------------------------------------------------------- 215 216 IF( npti > 0 ) THEN 217 218 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 219 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 220 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 221 222 DO jl = 1, jpl 223 224 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 225 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 226 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 227 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 228 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 229 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 230 231 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum (1:npti), dh_i_sum_2d (:,:,jl) ) 232 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt (1:npti), dh_s_mlt_2d (:,:,jl) ) 212 233 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 ---! 220 IF( zdv_mlt > 0._wp ) THEN 221 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 222 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 223 ! 224 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 225 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 226 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 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 ) 234 ! 235 !--- Pond contraction (due to refreezing) ---! 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 234 DO jk = 1, nlay_i 235 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 236 END DO 264 237 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 238 !----------------------------------------------------------------------------------------------- 239 ! Go for ponds 240 !----------------------------------------------------------------------------------------------- 241 242 243 DO ji = 1, npti 244 ! !----------------------------------------------------! 245 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 246 ! !----------------------------------------------------! 247 !--- Remove ponds on thin ice or tiny ice fractions 292 248 a_ip_1d(ji) = 0._wp 293 249 h_ip_1d(ji) = 0._wp 294 250 h_il_1d(ji) = 0._wp 251 ! !--------------------------------! 252 ELSE ! Case ice thickness >= rn_himin ! 253 ! !--------------------------------! 254 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 255 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 256 ! 257 !------------------! 258 ! case ice melting ! 259 !------------------! 260 ! 261 !--- available meltwater for melt ponding ---! 262 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 263 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 264 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 265 ! 266 !--- overflow ---! 267 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 268 ! a_ip_max = zfr_mlt * a_i 269 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 270 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 271 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 272 273 ! If pond depth exceeds half the ice thickness then reduce the pond volume 274 ! h_ip_max = 0.5 * h_i 275 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 276 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear 277 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 278 279 !--- Pond growing ---! 280 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 281 ! 282 !--- Lid melting ---! 283 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 284 ! 285 !--- mass flux ---! 286 ! MV I would recommend to remove that 287 ! Since melt ponds carry no freshwater there is no point in modifying water fluxes 288 289 IF( zdv_mlt > 0._wp ) THEN 290 zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] 291 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 292 ! 293 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 294 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 295 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 296 ENDIF 297 298 !-------------------! 299 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 300 !-------------------! 301 ! 302 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 303 ! 304 !--- Pond contraction (due to refreezing) ---! 305 IF( ln_pnd_lids ) THEN 306 ! 307 !--- Lid growing and subsequent pond shrinking ---! 308 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 309 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 310 311 ! Lid growing 312 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 313 314 ! Pond shrinking 315 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 316 317 ELSE 318 ! Pond shrinking 319 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 320 ENDIF 321 ! 322 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 323 ! v_ip = h_ip * a_ip 324 ! 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) 325 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 326 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 327 328 !---------------! 329 ! Pond flushing ! 330 !---------------! 331 ! height of top of the pond above sea-level 332 zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 333 334 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 335 DO jk = 1, nlay_i 336 ! MV Assur is inconsistent with SI3 337 zsbr = - 1.2_wp & 338 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 339 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 340 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 341 ! MV linear expression more consistent & simpler zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 342 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 343 END DO 344 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 345 346 ! Do the drainage using Darcy's law 347 zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 348 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 349 zdv_flush = 0._wp ! MV remove pond drainage for now 350 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 351 352 ! MV --- why pond drainage does not give back water into freshwater flux ? 353 354 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 355 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 356 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 357 358 !--- Corrections and lid thickness ---! 359 IF( ln_pnd_lids ) THEN 360 !--- retrieve lid thickness from volume ---! 361 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 362 ELSE ; h_il_1d(ji) = 0._wp 363 ENDIF 364 !--- remove ponds if lids are much larger than ponds ---! 365 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 366 a_ip_1d(ji) = 0._wp 367 h_ip_1d(ji) = 0._wp 368 h_il_1d(ji) = 0._wp 369 ENDIF 370 ENDIF 371 ! 295 372 ENDIF 373 374 END DO ! ji 375 376 !----------------------------------------------------------------------------------------------- 377 ! Retrieve 2D arrays 378 !----------------------------------------------------------------------------------------------- 379 380 v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti) 381 v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti) 382 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 383 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 384 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 385 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 386 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 387 DO jk = 1, nlay_i 388 CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 389 END DO 390 391 END DO ! ji 392 393 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum ) 394 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum ) 395 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 396 397 ! 398 ENDIF 399 400 END SUBROUTINE pnd_LEV 401 402 SUBROUTINE pnd_TOPO 403 404 !!------------------------------------------------------------------- 405 !! *** ROUTINE pnd_TOPO *** 406 !! 407 !! ** Purpose : Compute melt pond evolution 408 !! 409 !! ** Purpose : Compute melt pond evolution based on the ice 410 !! topography as inferred from the ice thickness 411 !! distribution. 412 !! 413 !! ** Method : This code is initially based on Flocco and Feltham 414 !! (2007) and Flocco et al. (2010). More to come... 415 !! 416 !! ** Tunable parameters : 417 !! 418 !! ** Note : 419 !! 420 !! ** References 421 !! 422 !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 423 !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 424 !! 10.1029/2006JC003836. 425 !! 426 !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 427 !! a physically based melt pond scheme into the sea ice component of a 428 !! climate model. J. Geophys. Res. 115, C08012, 429 !! doi: 10.1029/2009JC005568. 430 !! 431 !!------------------------------------------------------------------- 432 433 ! local variables 434 REAL(wp) :: & 435 zdHui, & ! change in thickness of ice lid (m) 436 zomega, & ! conduction 437 zdTice, & ! temperature difference across ice lid (C) 438 zdvice, & ! change in ice volume (m) 439 zTavg, & ! mean surface temperature across categories (C) 440 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 441 zTp, & ! pond freezing temperature (C) 442 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 443 zdvn , & ! change in melt pond volume for fresh water budget (not used for now) 444 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 445 z1_rhow, & ! inverse water density 446 zpond , & ! dummy variable 447 zdum ! dummy variable 448 449 REAL(wp), PARAMETER :: & 450 zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m) 451 zTd = 0.15_wp , & ! temperature difference for freeze-up (C) 452 zvp_min = 1.e-4_wp ! minimum pond volume (m) 453 454 REAL(wp), DIMENSION(jpi,jpj) :: zvolp !! total melt pond water available before redistribution and drainage 455 456 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 457 458 INTEGER :: ji, jj, jk, jl ! loop indices 459 460 INTEGER :: i_test 461 462 ! Note 463 ! equivalent for CICE translation 464 ! a_ip -> apond 465 ! a_ip_frac -> apnd 466 467 !--------------------------------------------------------------- 468 ! Initialize 469 !--------------------------------------------------------------- 470 471 ! Parameters & constants (move to parameters) 472 zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) 473 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 474 z1_rhow = 1._wp / rhow 475 476 ! Set required ice variables (hard-coded here for now) 477 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 478 479 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 480 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 481 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 482 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 483 484 ! thickness 485 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 486 ELSEWHERE ; z1_a_i(:,:,:) = 0._wp 487 END WHERE 488 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 489 490 !----------------------------------------------------------------- 491 ! Start pond loop 492 !----------------------------------------------------------------- 493 494 zvolp(:,:) = 0._wp 495 496 DO jj = 1, jpj 497 DO ji = 1, jpi 498 499 !-------------------------------------------------------------- 500 ! Collect meltwater on ponds 501 !-------------------------------------------------------------- 502 DO jl = 1, jpl 503 504 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 505 506 !--- Available meltwater for melt ponding ---! 507 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * a_i(ji,jj,jl) ! = ( 1 - r ) = fraction of melt water that is not flushed 508 zdum = -( dh_i_sum_2d(ji,jj,jl)*rhoi + dh_s_mlt_2d(ji,jj,jl)*rhos ) * z1_rhow * a_i(ji,jj,jl) 509 zpond = zfr_mlt * zdum ! check 510 511 !--- Create possible new ponds 512 ! if pond does not exist, create new pond over full ice area 513 IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 514 h_ip(ji,jj,jl) = 0._wp 515 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 516 ENDIF 517 518 !--- Deepen existing ponds before redistribution and drainage(later on) 519 h_ip(ji,jj,jl) = ( zpond + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) ) / a_ip_frac(ji,jj,jl) 520 zvolp(ji,jj) = zvolp(ji,jj) + h_ip(ji,jj,jl) * a_ip_frac(ji,jj,jl) 521 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) !useless 522 523 ENDIF ! a_i 524 525 END DO 526 527 h_ip(ji,jj,:) = 0._wp ! pond thickness 528 a_ip(ji,jj,:) = 0._wp ! pond fraction per grid area 529 530 !----------------------------------------------------------------- 531 ! Identify grid cells with ice 532 !----------------------------------------------------------------- 533 534 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. vt_ip(ji,jj) > zvp_min *at_i(ji,jj) ) THEN 535 536 !---------------------------------------------------------------------------- 537 ! Calculate pond area and depth (redistribution of melt water on topography) 538 !---------------------------------------------------------------------------- 539 ! CALL ice_thd_pnd_area( at_i(ji,jj) , vt_i(ji,jj), & 540 ! a_i (ji,jj,:), v_i(ji,jj,:), v_s(ji,jj,:), & 541 ! t_i(ji,jj,:,:), sz_i(ji,jj,:,:), & 542 ! v_ip(ji,jj,:) , zvolp(ji,jj), & 543 ! a_ip(ji,jj,:), h_ip(ji,jj,:), zdvn) 544 545 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvn ! ---> zdvn is drainage (useless for now since we dont change fw budget, no ?) 546 547 !-------------------------- 548 ! Pond lid growth and melt 549 !-------------------------- 550 ! Mean surface temperature 551 zTavg = 0._wp 552 DO jl = 1, jpl 553 zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 554 END DO 555 zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 556 557 DO jl = 1, jpl-1 558 559 IF ( v_il(ji,jj,jl) > epsi10 ) THEN 560 561 !---------------------------------------------------------------- 562 ! Lid melting: floating upper ice layer melts in whole or part 563 !---------------------------------------------------------------- 564 ! Use Tsfc for each category 565 IF ( t_su(ji,jj,jl) > zTp ) THEN 566 567 zdvice = min( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl)) 568 IF ( zdvice > epsi10 ) then 569 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 570 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice 571 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 572 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 573 574 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 575 ! ice lid melted and category is pond covered 576 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 577 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 578 v_il(ji,jj,jl) = 0._wp 579 ENDIF 580 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 581 ENDIF 582 583 !---------------------------------------------------------------- 584 ! Freeze pre-existing lid 585 !---------------------------------------------------------------- 586 587 ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 588 589 ! differential growth of base of surface floating ice layer 590 zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 591 zomega = rcnd_i * zdTice / zrhoi_L 592 zdHui = SQRT( 2._wp * zomega * rdt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 593 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 594 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 595 596 IF ( zdvice > epsi10 ) THEN 597 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 598 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 599 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 600 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 601 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 602 ENDIF 603 604 ENDIF ! Tsfcn(i,j,n) 605 606 !---------------------------------------------------------------- 607 ! Freeze new lids 608 !---------------------------------------------------------------- 609 ! upper ice layer begins to form 610 ! note: albedo does not change 611 612 ELSE ! v_il < epsi10 613 614 ! thickness of newly formed ice 615 ! the surface temperature of a meltpond is the same as that 616 ! of the ice underneath (0C), and the thermodynamic surface 617 ! flux is the same 618 619 !!!!! FSURF !!!!! 620 !!! we need net surface energy flux, excluding conduction 621 !!! fsurf is summed over categories in CICE 622 !!! we have the category-dependent flux, let us use it ? 623 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 624 zdHui = MAX ( -zfsurf * rdt_ice/zrhoi_L , 0._wp ) 625 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 626 IF ( zdvice > epsi10 ) THEN 627 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 628 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 629 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 630 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 631 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! in principle, this is useless as h_ip is computed in icevar 632 ENDIF 633 634 ENDIF ! v_il 635 636 END DO ! jl 637 638 ELSE ! remove ponds on thin ice 639 640 v_ip(ji,jj,:) = 0._wp 641 v_il (ji,jj,:) = 0._wp 642 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 643 ! zvolp(ji,jj) = 0._wp 644 296 645 ENDIF 297 ! 646 647 !--------------------------------------------------------------- 648 ! remove ice lid if there is no liquid pond 649 ! v_il may be nonzero on category ncat due to dynamics 650 !--------------------------------------------------------------- 651 652 DO jl = 1, jpl 653 654 IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 655 .AND. v_il (ji,jj,jl) > epsi10) THEN 656 v_il(ji,jj,jl) = 0._wp 657 ENDIF 658 659 ! reload tracers 660 IF ( a_ip(ji,jj,jl) > epsi10) THEN 661 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 662 ELSE 663 v_il(ji,jj,jl) = 0._wp 664 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as a_ip_frac computed in icevar 665 ENDIF 666 IF ( a_ip(ji,jj,jl) > epsi10 ) THEN 667 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 668 !h_ip(ji,jj,jl) = zhpondn(ji,jj,jl) 669 ELSE 670 a_ip_frac(ji,jj,jl) = 0._wp 671 h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 672 h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 298 673 ENDIF 299 674 300 END DO 301 ! 302 END SUBROUTINE pnd_LEV 303 675 END DO ! jl 676 677 END DO ! jj 678 END DO ! ji 679 680 END SUBROUTINE pnd_TOPO 681 682 683 SUBROUTINE ice_thd_pnd_area(aice, vice, & 684 aicen, vicen, vsnon, ticen, & 685 salin, zvolpn, zvolp, & 686 zapondn,zhpondn,dvolp) 687 688 !!------------------------------------------------------------------- 689 !! *** ROUTINE ice_thd_pnd_area *** 690 !! 691 !! ** Purpose : Given the total volume of meltwater, update 692 !! pond fraction (a_ip) and depth (should be volume) 693 !! 694 !! ** 695 !! 696 !!------------------------------------------------------------------ 697 698 REAL (wp), INTENT(IN) :: & 699 aice, vice 700 701 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 702 aicen, vicen, vsnon 703 704 REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & 705 ticen, salin 706 707 REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & 708 zvolpn 709 710 REAL (wp), INTENT(INOUT) :: & 711 zvolp, dvolp 712 713 REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & 714 zapondn, zhpondn 715 716 INTEGER :: & 717 n, ns, & 718 m_index, & 719 permflag 720 721 REAL (wp), DIMENSION(jpl) :: & 722 hicen, & 723 hsnon, & 724 asnon, & 725 alfan, & 726 betan, & 727 cum_max_vol, & 728 reduced_aicen 729 730 REAL (wp), DIMENSION(0:jpl) :: & 731 cum_max_vol_tmp 732 733 REAL (wp) :: & 734 hpond, & 735 drain, & 736 floe_weight, & 737 pressure_head, & 738 hsl_rel, & 739 deltah, & 740 perm, & 741 msno 742 743 REAL (wp), parameter :: & 744 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 745 746 !-----------| 747 ! | 748 ! |-----------| 749 !___________|___________|______________________________________sea-level 750 ! | | 751 ! | |---^--------| 752 ! | | | | 753 ! | | | |-----------| |------- 754 ! | | |alfan(n)| | | 755 ! | | | | |--------------| 756 ! | | | | | | 757 !---------------------------v------------------------------------------- 758 ! | | ^ | | | 759 ! | | | | |--------------| 760 ! | | |betan(n)| | | 761 ! | | | |-----------| |------- 762 ! | | | | 763 ! | |---v------- | 764 ! | | 765 ! |-----------| 766 ! | 767 !-----------| 768 769 !------------------------------------------------------------------- 770 ! initialize 771 !------------------------------------------------------------------- 772 773 DO n = 1, jpl 774 775 zapondn(n) = 0._wp 776 zhpondn(n) = 0._wp 777 778 !---------------------------------------- 779 ! X) compute the effective snow fraction 780 !---------------------------------------- 781 IF (aicen(n) < epsi10) THEN 782 hicen(n) = 0._wp 783 hsnon(n) = 0._wp 784 reduced_aicen(n) = 0._wp 785 asnon(n) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 786 ELSE 787 hicen(n) = vicen(n) / aicen(n) 788 hsnon(n) = vsnon(n) / aicen(n) 789 reduced_aicen(n) = 1._wp ! n=jpl 790 791 !js: initial code in NEMO_DEV 792 !IF (n < jpl) reduced_aicen(n) = aicen(n) & 793 ! * (-0.024_wp*hicen(n) + 0.832_wp) 794 795 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 796 IF (n < jpl) reduced_aicen(n) = aicen(n) & 797 * max(0.2_wp,(-0.024_wp*hicen(n) + 0.832_wp)) 798 799 asnon(n) = reduced_aicen(n) ! effective snow fraction (empirical) 800 ! MV should check whether this makes sense to have the same effective snow fraction in here 801 ! OLI: it probably doesn't 802 END IF 803 804 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 805 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 806 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 807 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 808 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 809 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 810 ! alfan = 60% of the ice volume) in each category lies above the reference line, 811 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 812 813 ! MV: 814 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 815 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 816 817 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 818 819 alfan(n) = 0.6 * hicen(n) 820 betan(n) = 0.4 * hicen(n) 821 822 cum_max_vol(n) = 0._wp 823 cum_max_vol_tmp(n) = 0._wp 824 825 END DO ! jpl 826 827 cum_max_vol_tmp(0) = 0._wp 828 drain = 0._wp 829 dvolp = 0._wp 830 831 !---------------------------------------------------------- 832 ! x) Drain overflow water, update pond fraction and volume 833 !---------------------------------------------------------- 834 835 !-------------------------------------------------------------------------- 836 ! the maximum amount of water that can be contained up to each ice category 837 !-------------------------------------------------------------------------- 838 839 ! MV 840 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 841 ! Then the excess volume cum_max_vol(jl) drains out of the system 842 ! It should be added to wfx_pnd_out 843 ! END MV 844 !js 18/04/19: XXX do something about this flux thing 845 846 DO n = 1, jpl-1 ! last category can not hold any volume 847 848 IF (alfan(n+1) >= alfan(n) .AND. alfan(n+1) > 0._wp ) THEN 849 850 ! total volume in level including snow 851 cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & 852 (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) 853 854 ! subtract snow solid volumes from lower categories in current level 855 DO ns = 1, n 856 cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & 857 - rhos/rhow * & ! free air fraction that can be filled by water 858 asnon(ns) * & ! effective areal fraction of snow in that category 859 max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)-alfan(n)), 0._wp) 860 END DO 861 862 ELSE ! assume higher categories unoccupied 863 cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) 864 END IF 865 !IF (cum_max_vol_tmp(n) < z0) THEN 866 ! CALL abort_ice('negative melt pond volume') 867 !END IF 868 END DO 869 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 870 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 871 872 !---------------------------------------------------------------- 873 ! is there more meltwater than can be held in the floe? 874 !---------------------------------------------------------------- 875 IF (zvolp >= cum_max_vol(jpl)) THEN 876 drain = zvolp - cum_max_vol(jpl) + epsi10 877 zvolp = zvolp - drain ! update meltwater volume available 878 dvolp = drain ! this is the drained water 879 IF (zvolp < epsi10) THEN 880 dvolp = dvolp + zvolp 881 zvolp = 0._wp 882 END IF 883 END IF 884 885 ! height and area corresponding to the remaining volume 886 887 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 888 889 DO n=1, m_index 890 !zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update 891 ! ! volume instead, no ? 892 zhpondn(n) = max((hpond - alfan(n) + alfan(1)), 0._wp) !js: from CICE 5.1.2 893 zapondn(n) = reduced_aicen(n) 894 ! in practise, pond fraction depends on the empirical snow fraction 895 ! so in turn on ice thickness 896 END DO 897 !zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 898 899 !------------------------------------------------------------------------ 900 ! Drainage through brine network (permeability) 901 !------------------------------------------------------------------------ 902 !!! drainage due to ice permeability - Darcy's law 903 904 ! sea water level 905 msno =0._wp 906 DO n=1,jpl 907 msno = msno + vsnon(n) * rhos 908 END DO 909 floe_weight = (msno + rhoi*vice + rau0*zvolp) / aice 910 hsl_rel = floe_weight / rau0 & 911 - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) 912 913 deltah = hpond - hsl_rel 914 pressure_head = grav * rau0 * max(deltah, 0._wp) 915 916 ! drain if ice is permeable 917 permflag = 0 918 IF (pressure_head > 0._wp) THEN 919 DO n = 1, jpl-1 920 IF (hicen(n) /= 0._wp) THEN 921 !IF (hicen(n) > 0._wp) THEN !js: from CICE 5.1.2 922 perm = 0._wp ! MV ugly dummy patch 923 CALL ice_thd_pnd_perm(ticen(:,n), salin(:,n), perm) 924 IF (perm > 0._wp) permflag = 1 925 926 drain = perm*zapondn(n)*pressure_head*rdt_ice / & 927 (viscosity*hicen(n)) 928 dvolp = dvolp + min(drain, zvolp) 929 zvolp = max(zvolp - drain, 0._wp) 930 IF (zvolp < epsi10) THEN 931 dvolp = dvolp + zvolp 932 zvolp = 0._wp 933 END IF 934 END IF 935 END DO 936 937 ! adjust melt pond dimensions 938 IF (permflag > 0) THEN 939 ! recompute pond depth 940 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 941 DO n=1, m_index 942 zhpondn(n) = hpond - alfan(n) + alfan(1) 943 zapondn(n) = reduced_aicen(n) 944 END DO 945 !zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 946 END IF 947 END IF ! pressure_head 948 949 !------------------------------- 950 ! X) remove water from the snow 951 !------------------------------- 952 !------------------------------------------------------------------------ 953 ! total melt pond volume in category does not include snow volume 954 ! snow in melt ponds is not melted 955 !------------------------------------------------------------------------ 956 957 ! Calculate pond volume for lower categories 958 DO n=1,m_index-1 959 zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow 960 - (rhos/rhow) * asnon(n) * min(hsnon(n), zhpondn(n)) 961 END DO 962 963 ! Calculate pond volume for highest category = remaining pond volume 964 965 ! The following is completely unclear to Martin at least 966 ! Could we redefine properly and recode in a more readable way ? 967 968 ! m_index = last category with melt pond 969 970 IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water 971 972 IF (m_index > 1) THEN 973 IF (zvolp > sum(zvolpn(1:m_index-1))) THEN 974 zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) 975 ELSE 976 zvolpn(m_index) = 0._wp 977 zhpondn(m_index) = 0._wp 978 zapondn(m_index) = 0._wp 979 ! If remaining pond volume is negative reduce pond volume of 980 ! lower category 981 IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & 982 zvolpn(m_index-1) = zvolpn(m_index-1) - sum(zvolpn(1:m_index-1)) + zvolp 983 END IF 984 END IF 985 986 DO n=1,m_index 987 IF (zapondn(n) > epsi10) THEN 988 zhpondn(n) = zvolpn(n) / zapondn(n) 989 ELSE 990 dvolp = dvolp + zvolpn(n) 991 zhpondn(n) = 0._wp 992 zvolpn(n) = 0._wp 993 zapondn(n) = 0._wp 994 END IF 995 END DO 996 DO n = m_index+1, jpl 997 zhpondn(n) = 0._wp 998 zapondn(n) = 0._wp 999 zvolpn (n) = 0._wp 1000 END DO 1001 1002 END SUBROUTINE ice_thd_pnd_area 1003 1004 1005 SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 1006 !!------------------------------------------------------------------- 1007 !! *** ROUTINE ice_thd_pnd_depth *** 1008 !! 1009 !! ** Purpose : Compute melt pond depth 1010 !!------------------------------------------------------------------- 1011 1012 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1013 aicen, & 1014 asnon, & 1015 hsnon, & 1016 alfan, & 1017 cum_max_vol 1018 1019 REAL (wp), INTENT(IN) :: & 1020 zvolp 1021 1022 REAL (wp), INTENT(OUT) :: & 1023 hpond 1024 1025 INTEGER, INTENT(OUT) :: & 1026 m_index 1027 1028 INTEGER :: n, ns 1029 1030 REAL (wp), DIMENSION(0:jpl+1) :: & 1031 hitl, & 1032 aicetl 1033 1034 REAL (wp) :: & 1035 rem_vol, & 1036 area, & 1037 vol, & 1038 tmp, & 1039 z0 = 0.0_wp 1040 1041 !---------------------------------------------------------------- 1042 ! hpond is zero if zvolp is zero - have we fully drained? 1043 !---------------------------------------------------------------- 1044 1045 IF (zvolp < epsi10) THEN 1046 hpond = z0 1047 m_index = 0 1048 ELSE 1049 1050 !---------------------------------------------------------------- 1051 ! Calculate the category where water fills up to 1052 !---------------------------------------------------------------- 1053 1054 !----------| 1055 ! | 1056 ! | 1057 ! |----------| -- -- 1058 !__________|__________|_________________________________________ ^ 1059 ! | | rem_vol ^ | Semi-filled 1060 ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1061 ! | | | | 1062 ! | | | |hpond 1063 ! | | |----------| | |------- 1064 ! | | | | | | 1065 ! | | | |---v-----| 1066 ! | | m_index | | | 1067 !------------------------------------------------------------- 1068 1069 m_index = 0 ! 1:m_index categories have water in them 1070 DO n = 1, jpl 1071 IF (zvolp <= cum_max_vol(n)) THEN 1072 m_index = n 1073 IF (n == 1) THEN 1074 rem_vol = zvolp 1075 ELSE 1076 rem_vol = zvolp - cum_max_vol(n-1) 1077 END IF 1078 exit ! to break out of the loop 1079 END IF 1080 END DO 1081 m_index = min(jpl-1, m_index) 1082 1083 !---------------------------------------------------------------- 1084 ! semi-filled layer may have m_index different snow in it 1085 !---------------------------------------------------------------- 1086 1087 !----------------------------------------------------------- ^ 1088 ! | alfan(m_index+1) 1089 ! | 1090 !hitl(3)--> |----------| | 1091 !hitl(2)--> |------------| * * * * *| | 1092 !hitl(1)--> |----------|* * * * * * |* * * * * | | 1093 !hitl(0)-->------------------------------------------------- | ^ 1094 ! various snow from lower categories | |alfa(m_index) 1095 1096 ! hitl - heights of the snow layers from thinner and current categories 1097 ! aicetl - area of each snow depth in this layer 1098 1099 hitl(:) = z0 1100 aicetl(:) = z0 1101 DO n = 1, m_index 1102 hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1103 alfan(m_index+1) - alfan(m_index)), z0) 1104 aicetl(n) = asnon(n) 1105 1106 aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1107 END DO 1108 1109 hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1110 aicetl(m_index+1) = z0 1111 1112 !---------------------------------------------------------------- 1113 ! reorder array according to hitl 1114 ! snow heights not necessarily in height order 1115 !---------------------------------------------------------------- 1116 1117 DO ns = 1, m_index+1 1118 DO n = 0, m_index - ns + 1 1119 IF (hitl(n) > hitl(n+1)) THEN ! swap order 1120 tmp = hitl(n) 1121 hitl(n) = hitl(n+1) 1122 hitl(n+1) = tmp 1123 tmp = aicetl(n) 1124 aicetl(n) = aicetl(n+1) 1125 aicetl(n+1) = tmp 1126 END IF 1127 END DO 1128 END DO 1129 1130 !---------------------------------------------------------------- 1131 ! divide semi-filled layer into set of sublayers each vertically homogenous 1132 !---------------------------------------------------------------- 1133 1134 !hitl(3)---------------------------------------------------------------- 1135 ! | * * * * * * * * 1136 ! |* * * * * * * * * 1137 !hitl(2)---------------------------------------------------------------- 1138 ! | * * * * * * * * | * * * * * * * * 1139 ! |* * * * * * * * * |* * * * * * * * * 1140 !hitl(1)---------------------------------------------------------------- 1141 ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1142 ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1143 !hitl(0)---------------------------------------------------------------- 1144 ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1145 1146 ! move up over layers incrementing volume 1147 DO n = 1, m_index+1 1148 1149 area = sum(aicetl(:)) - & ! total area of sub-layer 1150 (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1151 1152 vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1153 1154 IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1155 hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 1156 1157 exit 1158 ELSE ! still in sub-layer below the sub-layer with the depth 1159 rem_vol = rem_vol - vol 1160 END IF 1161 1162 END DO 1163 1164 END IF 1165 1166 END SUBROUTINE ice_thd_pnd_depth 1167 1168 1169 SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 1170 !!------------------------------------------------------------------- 1171 !! *** ROUTINE ice_thd_pnd_perm *** 1172 !! 1173 !! ** Purpose : Determine the liquid fraction of brine in the ice 1174 !! and its permeability 1175 !!------------------------------------------------------------------- 1176 1177 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1178 ticen, & ! internal ice temperature (K) 1179 salin ! salinity (ppt) !js: ppt according to cice 1180 1181 REAL (wp), INTENT(OUT) :: & 1182 perm ! permeability 1183 1184 REAL (wp) :: & 1185 Sbr ! brine salinity 1186 1187 REAL (wp), DIMENSION(nlay_i) :: & 1188 Tin, & ! ice temperature 1189 phi ! liquid fraction 1190 1191 INTEGER :: k 1192 1193 !----------------------------------------------------------------- 1194 ! Compute ice temperatures from enthalpies using quadratic formula 1195 !----------------------------------------------------------------- 1196 1197 DO k = 1,nlay_i 1198 Tin(k) = ticen(k) - rt0 !js: from K to degC 1199 END DO 1200 1201 !----------------------------------------------------------------- 1202 ! brine salinity and liquid fraction 1203 !----------------------------------------------------------------- 1204 1205 DO k = 1, nlay_i 1206 1207 Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 1208 ! Best expression to date is that one 1209 Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 1210 phi(k) = salin(k) / Sbr 1211 1212 END DO 1213 1214 !----------------------------------------------------------------- 1215 ! permeability 1216 !----------------------------------------------------------------- 1217 1218 perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 1219 1220 END SUBROUTINE ice_thd_pnd_perm 1221 1222 1223 !---------------------------------------------------------------------------------------------------------------------- 304 1224 305 1225 SUBROUTINE ice_thd_pnd_init … … 319 1239 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & 320 1240 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1241 & ln_pnd_TOPO , & 321 1242 & ln_pnd_lids, ln_pnd_alb 322 1243 !!------------------------------------------------------------------- … … 336 1257 WRITE(numout,*) ' Namelist namicethd_pnd:' 337 1258 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 1259 WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO 338 1260 WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV 339 1261 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min … … 351 1273 IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF 352 1274 IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF 1275 IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF 353 1276 IF( ioptio /= 1 ) & 354 1277 & 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)' )
Note: See TracChangeset
for help on using the changeset viewer.
