Changeset 14789 for NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icethd_pnd.F90
- Timestamp:
- 2021-05-05T13:18:04+02:00 (3 years ago)
- Location:
- NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Files:
-
- 2 edited
-
. (modified) (1 prop)
-
src/ICE/icethd_pnd.F90 (modified) (14 diffs)
Legend:
- Unmodified
- Added
- Removed
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU
- Property svn:externals
-
old new 3 3 ^/utils/build/mk@HEAD mk 4 4 ^/utils/tools@HEAD tools 5 ^/vendors/AGRIF/dev _r12970_AGRIF_CMEMSext/AGRIF5 ^/vendors/AGRIF/dev@HEAD ext/AGRIF 6 6 ^/vendors/FCM@HEAD ext/FCM 7 7 ^/vendors/IOIPSL@HEAD ext/IOIPSL 8 ^/vendors/PPR@HEAD ext/PPR 8 9 9 10 # SETTE 10 ^/utils/CI/sette@1 3559sette11 ^/utils/CI/sette@14244 sette
-
- Property svn:externals
-
NEMO/branches/2021/dev_r13747_HPC-11_mcastril_HPDAonline_DiagGPU/src/ICE/icethd_pnd.F90
r13472 r14789 1 MODULE icethd_pnd 1 MODULE icethd_pnd 2 2 !!====================================================================== 3 3 !! *** MODULE icethd_pnd *** … … 20 20 USE ice1D ! sea-ice: thermodynamics variables 21 21 USE icetab ! sea-ice: 1D <==> 2D transformation 22 USE sbc_ice ! surface energy budget 22 23 ! 23 24 USE in_out_manager ! I/O manager 25 USE iom ! I/O manager library 24 26 USE lib_mpp ! MPP library 25 27 USE lib_fortran ! fortran utilities (glob_sum + no signed zero) … … 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 39 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 42 43 !-------------------------------------------------------------------------- 44 ! Diagnostics for pond volume per area 45 ! 46 ! dV/dt = mlt + drn + lid + rnf 47 ! mlt = input from surface melting 48 ! drn = drainage through brine network 49 ! lid = lid growth & melt 50 ! rnf = runoff (water directly removed out of surface melting + overflow) 51 ! 52 ! In topo mode, the pond water lost because it is in the snow is not included in the budget 53 ! In level mode, all terms are incorporated 54 ! 55 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_mlt ! meltwater pond volume input [kg/m2/s] 56 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_drn ! pond volume lost by drainage [-] 57 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_lid ! exchange with lid / refreezing [-] 58 REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_rnf ! meltwater pond lost to runoff [-] 59 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_mlt_1d ! meltwater pond volume input [-] 60 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_drn_1d ! pond volume lost by drainage [-] 61 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_lid_1d ! exchange with lid / refreezing [-] 62 REAL(wp), ALLOCATABLE, DIMENSION(:) :: diag_dvpn_rnf_1d ! meltwater pond lost to runoff [-] 63 64 !! * Substitutions 65 # include "do_loop_substitute.h90" 40 66 !!---------------------------------------------------------------------- 41 67 !! NEMO/ICE 4.0 , NEMO Consortium (2018) … … 46 72 47 73 SUBROUTINE ice_thd_pnd 74 48 75 !!------------------------------------------------------------------- 49 76 !! *** ROUTINE ice_thd_pnd *** 50 !! 77 !! 51 78 !! ** Purpose : change melt pond fraction and thickness 52 !! 79 !! 80 !! ** Note : Melt ponds affect only radiative transfer for now 81 !! No heat, no salt. 82 !! The current diagnostics lacks a contribution from drainage 53 83 !!------------------------------------------------------------------- 84 INTEGER :: ji, jj, jl ! loop indices 85 !!------------------------------------------------------------------- 86 87 ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) ) 88 ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) ) 54 89 ! 55 SELECT CASE ( nice_pnd ) 90 diag_dvpn_mlt (:,:) = 0._wp ; diag_dvpn_drn (:,:) = 0._wp 91 diag_dvpn_lid (:,:) = 0._wp ; diag_dvpn_rnf (:,:) = 0._wp 92 diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp 93 diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp 94 95 !------------------------------------- 96 ! Remove ponds where ice has vanished 97 !------------------------------------- 98 at_i(:,:) = SUM( a_i, dim=3 ) 56 99 ! 57 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 58 ! 59 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 60 ! 61 END SELECT 100 DO jl = 1, jpl 101 DO_2D( 1, 1, 1, 1 ) 102 IF( v_i(ji,jj,jl) < epsi10 .OR. at_i(ji,jj) < epsi10 ) THEN 103 wfx_pnd (ji,jj) = wfx_pnd(ji,jj) + ( v_ip(ji,jj,jl) + v_il(ji,jj,jl) ) * rhow * r1_Dt_ice 104 a_ip (ji,jj,jl) = 0._wp 105 v_ip (ji,jj,jl) = 0._wp 106 v_il (ji,jj,jl) = 0._wp 107 h_ip (ji,jj,jl) = 0._wp 108 h_il (ji,jj,jl) = 0._wp 109 a_ip_frac(ji,jj,jl) = 0._wp 110 ENDIF 111 END_2D 112 END DO 113 114 !------------------------------ 115 ! Identify grid cells with ice 116 !------------------------------ 117 npti = 0 ; nptidx(:) = 0 118 DO_2D( 1, 1, 1, 1 ) 119 IF( at_i(ji,jj) >= epsi10 ) THEN 120 npti = npti + 1 121 nptidx( npti ) = (jj - 1) * jpi + ji 122 ENDIF 123 END_2D 124 125 !------------------------------------ 126 ! Select melt pond scheme to be used 127 !------------------------------------ 128 IF( npti > 0 ) THEN 129 SELECT CASE ( nice_pnd ) 130 ! 131 CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! 132 ! 133 CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! 134 ! 135 CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==! 136 ! 137 END SELECT 138 ENDIF 139 140 !------------------------------------ 141 ! Diagnostics 142 !------------------------------------ 143 CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt ) ! input from melting 144 CALL iom_put( 'dvpn_lid', diag_dvpn_lid ) ! exchanges with lid 145 CALL iom_put( 'dvpn_drn', diag_dvpn_drn ) ! vertical drainage 146 CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf ) ! runoff + overflow 62 147 ! 63 END SUBROUTINE ice_thd_pnd 64 65 66 SUBROUTINE pnd_CST 148 DEALLOCATE( diag_dvpn_mlt , diag_dvpn_lid , diag_dvpn_drn , diag_dvpn_rnf ) 149 DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d ) 150 151 END SUBROUTINE ice_thd_pnd 152 153 154 SUBROUTINE pnd_CST 67 155 !!------------------------------------------------------------------- 68 156 !! *** ROUTINE pnd_CST *** … … 70 158 !! ** Purpose : Compute melt pond evolution 71 159 !! 72 !! ** Method : Melt pond fraction and thickness are prescribed 160 !! ** Method : Melt pond fraction and thickness are prescribed 73 161 !! to non-zero values when t_su = 0C 74 162 !! 75 163 !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd) 76 !! 164 !! 77 165 !! ** Note : Coupling with such melt ponds is only radiative 78 166 !! Advection, ridging, rafting... are bypassed … … 80 168 !! ** References : Bush, G.W., and Trump, D.J. (2017) 81 169 !!------------------------------------------------------------------- 82 INTEGER :: ji ! loop indices 170 INTEGER :: ji, jl ! loop indices 171 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 83 172 !!------------------------------------------------------------------- 84 DO ji = 1, npti 85 ! 86 IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN 87 h_ip_1d(ji) = rn_hpnd 88 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 89 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 90 ELSE 91 h_ip_1d(ji) = 0._wp 92 a_ip_1d(ji) = 0._wp 93 h_il_1d(ji) = 0._wp 94 ENDIF 95 ! 173 DO jl = 1, jpl 174 175 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 176 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) ) 177 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 178 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 179 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 180 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 181 182 DO ji = 1, npti 183 ! 184 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 185 ! 186 IF( a_i_1d(ji) >= 0.01_wp .AND. t_su_1d(ji) >= rt0 ) THEN 187 h_ip_1d(ji) = rn_hpnd 188 a_ip_1d(ji) = rn_apnd * a_i_1d(ji) 189 h_il_1d(ji) = 0._wp ! no pond lids whatsoever 190 ELSE 191 h_ip_1d(ji) = 0._wp 192 a_ip_1d(ji) = 0._wp 193 h_il_1d(ji) = 0._wp 194 ENDIF 195 ! 196 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 197 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 198 ! 199 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 200 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 201 ! 202 END DO 203 204 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) ) 205 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) ) 206 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) ) 207 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) ) 208 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) ) 209 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd(:,:) ) 210 96 211 END DO 97 212 ! … … 107 222 !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing 108 223 !! We work with volumes and then redistribute changes into thickness and concentration 109 !! assuming linear relationship between the two. 224 !! assuming linear relationship between the two. 110 225 !! 111 226 !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- … … 122 237 !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- 123 238 !! 124 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 239 !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H 125 240 !! H = lid thickness 126 241 !! Lf = latent heat of fusion … … 132 247 !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- 133 248 !! 134 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi 135 !! perm = permability of sea-ice 249 !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi * flush --- from Flocco et al 2007 --- 250 !! perm = permability of sea-ice + correction from Hunke et al 2012 (flush) 136 251 !! visc = water viscosity 137 252 !! Hp = height of top of the pond above sea-level 138 253 !! Hi = ice thickness thru which there is flushing 254 !! flush= correction otherwise flushing is excessive 139 255 !! 140 256 !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness … … 143 259 !! a_ip/a_i = a_ip_frac = h_ip / zaspect 144 260 !! 145 !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min146 !! 147 !! ** Note : mostly stolen from CICE261 !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush 262 !! 263 !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. 148 264 !! 149 265 !! ** References : Flocco and Feltham (JGR, 2007) 150 266 !! Flocco et al (JGR, 2010) 151 267 !! Holland et al (J. Clim, 2012) 268 !! Hunke et al (OM 2012) 152 269 !!------------------------------------------------------------------- 153 270 REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array … … 157 274 REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity 158 275 !! 159 REAL(wp) :: zfr_mlt, zdv_mlt 276 REAL(wp) :: zfr_mlt, zdv_mlt, zdv_avail ! fraction and volume of available meltwater retained for melt ponding 160 277 REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes 278 REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids 161 279 REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh 162 280 REAL(wp) :: zv_ip_max ! max pond volume allowed 163 281 REAL(wp) :: zdT ! zTp-t_su 164 REAL(wp) :: zsbr 282 REAL(wp) :: zsbr, ztmelts ! Brine salinity 165 283 REAL(wp) :: zperm ! permeability of sea ice 166 284 REAL(wp) :: zfac, zdum ! temporary arrays 167 285 REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse 168 286 !! 169 INTEGER :: ji, jk 287 INTEGER :: ji, jk, jl ! loop indices 170 288 !!------------------------------------------------------------------- 171 z1_rhow = 1._wp / rhow 289 z1_rhow = 1._wp / rhow 172 290 z1_aspect = 1._wp / zaspect 173 z1_Tp = 1._wp / zTp 174 175 DO ji = 1, npti 176 ! !----------------------------------------------------! 177 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 178 ! !----------------------------------------------------! 179 !--- Remove ponds on thin ice or tiny ice fractions 180 a_ip_1d(ji) = 0._wp 181 h_ip_1d(ji) = 0._wp 182 h_il_1d(ji) = 0._wp 183 ! !--------------------------------! 184 ELSE ! Case ice thickness >= rn_himin ! 185 ! !--------------------------------! 186 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 291 z1_Tp = 1._wp / zTp 292 293 CALL tab_2d_1d( npti, nptidx(1:npti), at_i_1d (1:npti), at_i ) 294 CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd ) 295 296 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 297 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 298 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 299 CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 300 301 DO jl = 1, jpl 302 303 CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) ) 304 CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) ) 305 CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su(:,:,jl) ) 306 CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip(:,:,jl) ) 307 CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip(:,:,jl) ) 308 CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il(:,:,jl) ) 309 310 CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum(1:npti), dh_i_sum_2d(:,:,jl) ) 311 CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt(1:npti), dh_s_mlt_2d(:,:,jl) ) 312 313 DO jk = 1, nlay_i 314 CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) ) 315 CALL tab_2d_1d( npti, nptidx(1:npti), t_i_1d (1:npti,jk), t_i (:,:,jk,jl) ) 316 END DO 317 318 !----------------------- 319 ! Melt pond calculations 320 !----------------------- 321 DO ji = 1, npti 322 ! 323 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) 324 ! !----------------------------------------------------! 325 IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < 0.01_wp ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! 326 ! !----------------------------------------------------! 327 !--- Remove ponds on thin ice or tiny ice fractions 328 a_ip_1d(ji) = 0._wp 329 h_ip_1d(ji) = 0._wp 330 h_il_1d(ji) = 0._wp 331 ! !--------------------------------! 332 ELSE ! Case ice thickness >= rn_himin ! 333 ! !--------------------------------! 334 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness 335 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 336 ! 337 !------------------! 338 ! case ice melting ! 339 !------------------! 340 ! 341 !--- available meltwater for melt ponding (zdv_avail) ---! 342 zdv_avail = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) ! > 0 343 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 344 zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? 345 ! 346 !--- overflow ---! 347 ! 348 ! area driven overflow 349 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 350 ! a_ip_max = zfr_mlt * a_i 351 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 352 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 353 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 354 355 ! depth driven overflow 356 ! If pond depth exceeds half the ice thickness then reduce the pond volume 357 ! h_ip_max = 0.5 * h_i 358 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 359 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 360 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 361 362 !--- Pond growing ---! 363 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 364 ! 365 !--- Lid melting ---! 366 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 367 ! 368 !-------------------! 369 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 370 !-------------------! 371 ! 372 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 373 ! 374 !--- Pond contraction (due to refreezing) ---! 375 IF( ln_pnd_lids ) THEN 376 ! 377 !--- Lid growing and subsequent pond shrinking ---! 378 zdv_frz = - 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 379 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rDt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 380 381 ! Lid growing 382 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) 383 384 ! Pond shrinking 385 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 386 387 ELSE 388 zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) 389 ! Pond shrinking 390 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) 391 ENDIF 392 ! 393 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 394 ! v_ip = h_ip * a_ip 395 ! 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) 396 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 397 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 398 ! 399 400 !------------------------------------------------! 401 ! Pond drainage through brine network (flushing) ! 402 !------------------------------------------------! 403 ! height of top of the pond above sea-level 404 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 405 406 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 407 DO jk = 1, nlay_i 408 ! MV Assur is inconsistent with SI3 409 !!zsbr = - 1.2_wp & 410 !! & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 411 !! & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 412 !! & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 413 !!ztmp(jk) = sz_i_1d(ji,jk) / zsbr 414 ! MV linear expression more consistent & simpler: zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt 415 ztmelts = -rTmlt * sz_i_1d(ji,jk) 416 ztmp(jk) = ztmelts / MIN( ztmelts, t_i_1d(ji,jk) - rt0 ) 417 END DO 418 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 419 420 ! Do the drainage using Darcy's law 421 zdv_flush = -zperm * rho0 * grav * zhp * rDt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush ! zflush comes from Hunke et al. (2012) 422 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 423 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 424 425 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 426 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 427 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 428 429 !--- Corrections and lid thickness ---! 430 IF( ln_pnd_lids ) THEN 431 !--- retrieve lid thickness from volume ---! 432 IF( a_ip_1d(ji) > 0.01_wp ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 433 ELSE ; h_il_1d(ji) = 0._wp 434 ENDIF 435 !--- remove ponds if lids are much larger than ponds ---! 436 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 437 a_ip_1d(ji) = 0._wp 438 h_ip_1d(ji) = 0._wp 439 h_il_1d(ji) = 0._wp 440 ENDIF 441 ENDIF 442 443 ! diagnostics: dvpnd = mlt+rnf+lid+drn 444 diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + rhow * zdv_avail * r1_Dt_ice ! > 0, surface melt input 445 diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + rhow * ( zdv_mlt - zdv_avail ) * r1_Dt_ice ! < 0, runoff 446 diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + rhow * zdv_frz * r1_Dt_ice ! < 0, shrinking 447 diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + rhow * zdv_flush * r1_Dt_ice ! < 0, drainage 448 ! 449 ENDIF 450 ! 451 v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) 187 452 v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) 188 453 ! 189 !------------------! 190 ! case ice melting ! 191 !------------------! 454 zdv_pnd = ( h_ip_1d(ji) + h_il_1d(ji) ) * a_ip_1d(ji) - zdv_pnd 455 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zdv_pnd * rhow * r1_Dt_ice 192 456 ! 193 !--- available meltwater for melt ponding ---! 194 zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 195 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 196 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors? 197 ! 198 !--- overflow ---! 199 ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume 200 ! a_ip_max = zfr_mlt * a_i 201 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 202 zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect 203 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 204 205 ! If pond depth exceeds half the ice thickness then reduce the pond volume 206 ! h_ip_max = 0.5 * h_i 207 ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 208 zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) 209 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) 210 211 !--- Pond growing ---! 212 v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt 213 ! 214 !--- Lid melting ---! 215 IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 216 ! 217 !--- mass flux ---! 218 IF( zdv_mlt > 0._wp ) THEN 219 zfac = zdv_mlt * rhow * r1_Dt_ice ! melt pond mass flux < 0 [kg.m-2.s-1] 220 wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac 221 ! 222 zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux 223 wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) 224 wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) 225 ENDIF 226 227 !-------------------! 228 ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) 229 !-------------------! 230 ! 231 zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) 232 ! 233 !--- Pond contraction (due to refreezing) ---! 234 IF( ln_pnd_lids ) THEN 235 ! 236 !--- Lid growing and subsequent pond shrinking ---! 237 zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 238 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors 239 240 ! Lid growing 241 v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) 242 243 ! Pond shrinking 244 v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) 245 246 ELSE 247 ! Pond shrinking 248 v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) 249 ENDIF 250 ! 251 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 252 ! v_ip = h_ip * a_ip 253 ! 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) 254 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 255 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 256 257 !---------------! 258 ! Pond flushing ! 259 !---------------! 260 ! height of top of the pond above sea-level 261 zhp = ( h_i_1d(ji) * ( rho0 - rhoi ) + h_ip_1d(ji) * ( rho0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rho0 262 263 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) 264 DO jk = 1, nlay_i 265 zsbr = - 1.2_wp & 266 & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & 267 & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & 268 & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 269 ztmp(jk) = sz_i_1d(ji,jk) / zsbr 270 END DO 271 zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) 272 273 ! Do the drainage using Darcy's law 274 zdv_flush = -zperm * rho0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) 275 zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) 276 v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush 277 278 !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 279 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 280 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) 281 282 !--- Corrections and lid thickness ---! 283 IF( ln_pnd_lids ) THEN 284 !--- retrieve lid thickness from volume ---! 285 IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) 286 ELSE ; h_il_1d(ji) = 0._wp 287 ENDIF 288 !--- remove ponds if lids are much larger than ponds ---! 289 IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN 290 a_ip_1d(ji) = 0._wp 291 h_ip_1d(ji) = 0._wp 292 h_il_1d(ji) = 0._wp 293 ENDIF 294 ENDIF 295 ! 296 ENDIF 297 457 END DO 458 459 !-------------------------------------------------------------------- 460 ! Retrieve 2D arrays 461 !-------------------------------------------------------------------- 462 CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d(1:npti), a_ip(:,:,jl) ) 463 CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d(1:npti), h_ip(:,:,jl) ) 464 CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d(1:npti), h_il(:,:,jl) ) 465 CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d(1:npti), v_ip(:,:,jl) ) 466 CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d(1:npti), v_il(:,:,jl) ) 467 ! 298 468 END DO 299 469 ! 470 CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d(1:npti), wfx_pnd ) 471 ! 472 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt ) 473 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn ) 474 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid ) 475 CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf ) 476 ! 300 477 END SUBROUTINE pnd_LEV 301 478 302 479 303 SUBROUTINE ice_thd_pnd_init 480 481 SUBROUTINE pnd_TOPO 482 483 !!------------------------------------------------------------------- 484 !! *** ROUTINE pnd_TOPO *** 485 !! 486 !! ** Purpose : Compute melt pond evolution based on the ice 487 !! topography inferred from the ice thickness distribution 488 !! 489 !! ** Method : This code is initially based on Flocco and Feltham 490 !! (2007) and Flocco et al. (2010). 491 !! 492 !! - Calculate available pond water base on surface meltwater 493 !! - Redistribute water as a function of topography, drain water 494 !! - Exchange water with the lid 495 !! 496 !! ** Tunable parameters : 497 !! 498 !! ** Note : 499 !! 500 !! ** References 501 !! 502 !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond 503 !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: 504 !! 10.1029/2006JC003836. 505 !! 506 !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of 507 !! a physically based melt pond scheme into the sea ice component of a 508 !! climate model. J. Geophys. Res. 115, C08012, 509 !! doi: 10.1029/2009JC005568. 510 !! 511 !!------------------------------------------------------------------- 512 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 513 zTd = 273._wp , & ! temperature difference for freeze-up (K) 514 zvp_min = 1.e-4_wp ! minimum pond volume (m) 515 516 517 ! local variables 518 REAL(wp) :: & 519 zdHui, & ! change in thickness of ice lid (m) 520 zomega, & ! conduction 521 zdTice, & ! temperature difference across ice lid (C) 522 zdvice, & ! change in ice volume (m) 523 zTavg, & ! mean surface temperature across categories (C) 524 zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2) 525 zTp, & ! pond freezing temperature (C) 526 zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3) 527 zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding 528 z1_rhow, & ! inverse water density 529 zv_pnd , & ! volume of meltwater contributing to ponds 530 zv_mlt ! total amount of meltwater produced 531 532 REAL(wp), DIMENSION(jpi,jpj) :: zvolp_ini , & !! total melt pond water available before redistribution and drainage 533 zvolp , & !! total melt pond water volume 534 zvolp_res !! remaining melt pond water available after drainage 535 536 REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i 537 538 INTEGER :: ji, jj, jk, jl ! loop indices 539 540 INTEGER :: i_test 541 542 ! Note 543 ! equivalent for CICE translation 544 ! a_ip -> apond 545 ! a_ip_frac -> apnd 546 547 CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) 548 549 !--------------------------------------------------------------- 550 ! Initialise 551 !--------------------------------------------------------------- 552 553 ! Parameters & constants (move to parameters) 554 zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) 555 zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) 556 z1_rhow = 1._wp / rhow 557 558 ! Set required ice variables (hard-coded here for now) 559 ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) 560 561 at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction 562 vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area 563 vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area 564 vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area 565 566 ! thickness 567 WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) 568 ELSEWHERE ; z1_a_i(:,:,:) = 0._wp 569 END WHERE 570 h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) 571 572 !--------------------------------------------------------------- 573 ! Change 2D to 1D 574 !--------------------------------------------------------------- 575 ! MV 576 ! a less computing-intensive version would have 2D-1D passage here 577 ! use what we have in iceitd.F90 (incremental remapping) 578 579 !-------------------------------------------------------------- 580 ! Collect total available pond water volume 581 !-------------------------------------------------------------- 582 ! Assuming that meltwater (+rain in principle) runsoff the surface 583 ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction 584 ! I cite her words, they are very talkative 585 ! "grid cells with very little ice cover (and hence more open water area) 586 ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." 587 ! "This results in the same runoff fraction r for each ice category within a grid cell" 588 589 zvolp(:,:) = 0._wp 590 591 DO jl = 1, jpl 592 DO_2D( 1, 1, 1, 1 ) 593 594 IF ( a_i(ji,jj,jl) > epsi10 ) THEN 595 596 !--- Available and contributing meltwater for melt ponding ---! 597 zv_mlt = - ( dh_i_sum_2d(ji,jj,jl) * rhoi + dh_s_mlt_2d(ji,jj,jl) * rhos ) & ! available volume of surface melt water per grid area 598 & * z1_rhow * a_i(ji,jj,jl) 599 ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area 600 zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds 601 zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl 602 603 diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice ! diags 604 diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice 605 606 !--- Create possible new ponds 607 ! if pond does not exist, create new pond over full ice area 608 !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN 609 IF ( a_ip(ji,jj,jl) < epsi10 ) THEN 610 a_ip(ji,jj,jl) = a_i(ji,jj,jl) 611 a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) 612 ENDIF 613 614 !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage 615 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness 616 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 617 618 !--- Total available pond water volume (pre-existing + newly produced)j 619 zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) 620 ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now 621 622 ENDIF ! a_i 623 624 END_2D 625 END DO ! ji 626 627 zvolp_ini(:,:) = zvolp(:,:) 628 629 !-------------------------------------------------------------- 630 ! Redistribute and drain water from ponds 631 !-------------------------------------------------------------- 632 CALL ice_thd_pnd_area( zvolp, zvolp_res ) 633 634 !-------------------------------------------------------------- 635 ! Melt pond lid growth and melt 636 !-------------------------------------------------------------- 637 638 IF( ln_pnd_lids ) THEN 639 640 DO_2D( 1, 1, 1, 1 ) 641 642 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp_ini(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 643 644 !-------------------------- 645 ! Pond lid growth and melt 646 !-------------------------- 647 ! Mean surface temperature 648 zTavg = 0._wp 649 DO jl = 1, jpl 650 zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl) 651 END DO 652 zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here 653 654 DO jl = 1, jpl-1 655 656 IF ( v_il(ji,jj,jl) > epsi10 ) THEN 657 658 !---------------------------------------------------------------- 659 ! Lid melting: floating upper ice layer melts in whole or part 660 !---------------------------------------------------------------- 661 ! Use Tsfc for each category 662 IF ( t_su(ji,jj,jl) > zTp ) THEN 663 664 zdvice = MIN( dh_i_sum_2d(ji,jj,jl)*a_ip(ji,jj,jl), v_il(ji,jj,jl) ) 665 666 IF ( zdvice > epsi10 ) THEN 667 668 v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice 669 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - 670 ! as it is already counted in surface melt 671 ! zvolp(ji,jj) = zvolp(ji,jj) + zdvice ! pointless to calculate total volume (done in icevar) 672 ! zfpond(ji,jj) = fpond(ji,jj) + zdvice ! pointless to follow fw budget (ponds have no fw) 673 674 IF ( v_il(ji,jj,jl) < epsi10 .AND. v_ip(ji,jj,jl) > epsi10) THEN 675 ! ice lid melted and category is pond covered 676 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + v_il(ji,jj,jl) 677 ! zfpond(ji,jj) = zfpond (ji,jj) + v_il(ji,jj,jl) 678 v_il(ji,jj,jl) = 0._wp 679 ENDIF 680 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) !!! could get a division by zero here 681 682 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag 683 684 ENDIF 685 686 !---------------------------------------------------------------- 687 ! Freeze pre-existing lid 688 !---------------------------------------------------------------- 689 690 ELSE IF ( v_ip(ji,jj,jl) > epsi10 ) THEN ! Tsfcn(i,j,n) <= Tp 691 692 ! differential growth of base of surface floating ice layer 693 zdTice = MAX( - ( t_su(ji,jj,jl) - zTd ) , 0._wp ) ! > 0 694 zomega = rcnd_i * zdTice / zrhoi_L 695 zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & 696 - v_il(ji,jj,jl) / a_i(ji,jj,jl) 697 zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 698 699 IF ( zdvice > epsi10 ) THEN 700 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 701 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 702 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 703 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 704 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 705 706 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 707 708 ENDIF 709 710 ENDIF ! Tsfcn(i,j,n) 711 712 !---------------------------------------------------------------- 713 ! Freeze new lids 714 !---------------------------------------------------------------- 715 ! upper ice layer begins to form 716 ! note: albedo does not change 717 718 ELSE ! v_il < epsi10 719 720 ! thickness of newly formed ice 721 ! the surface temperature of a meltpond is the same as that 722 ! of the ice underneath (0C), and the thermodynamic surface 723 ! flux is the same 724 725 !!! we need net surface energy flux, excluding conduction 726 !!! fsurf is summed over categories in CICE 727 !!! we have the category-dependent flux, let us use it ? 728 zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) 729 zdHui = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp ) 730 zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) 731 IF ( zdvice > epsi10 ) THEN 732 v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice 733 v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice 734 735 diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag 736 ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice 737 ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice 738 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) ! MV - in principle, this is useless as h_ip is computed in icevar 739 ENDIF 740 741 ENDIF ! v_il 742 743 END DO ! jl 744 745 ELSE ! remove ponds on thin ice 746 747 v_ip(ji,jj,:) = 0._wp 748 v_il(ji,jj,:) = 0._wp 749 ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) 750 ! zvolp(ji,jj) = 0._wp 751 752 ENDIF 753 754 END_2D 755 756 ENDIF ! ln_pnd_lids 757 758 !--------------------------------------------------------------- 759 ! Clean-up variables (probably duplicates what icevar would do) 760 !--------------------------------------------------------------- 761 ! MV comment 762 ! In the ideal world, the lines above should update only v_ip, a_ip, v_il 763 ! icevar should recompute all other variables (if needed at all) 764 765 DO jl = 1, jpl 766 767 DO_2D( 1, 1, 1, 1 ) 768 769 ! ! zap lids on small ponds 770 ! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & 771 ! .AND. v_il(ji,jj,jl) > epsi10) THEN 772 ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small 773 ! ENDIF 774 775 ! recalculate equivalent pond variables 776 IF ( a_ip(ji,jj,jl) > epsi10) THEN 777 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_i(ji,jj,jl) 778 a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar 779 h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar 780 ENDIF 781 ! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar 782 ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar 783 ! ENDIF 784 785 END_2D 786 787 END DO ! jl 788 789 790 END SUBROUTINE pnd_TOPO 791 792 793 SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) 794 795 !!------------------------------------------------------------------- 796 !! *** ROUTINE ice_thd_pnd_area *** 797 !! 798 !! ** Purpose : Given the total volume of available pond water, 799 !! redistribute and drain water 800 !! 801 !! ** Method 802 !! 803 !-----------| 804 ! | 805 ! |-----------| 806 !___________|___________|______________________________________sea-level 807 ! | | 808 ! | |---^--------| 809 ! | | | | 810 ! | | | |-----------| |------- 811 ! | | | alfan | | | 812 ! | | | | |--------------| 813 ! | | | | | | 814 !---------------------------v------------------------------------------- 815 ! | | ^ | | | 816 ! | | | | |--------------| 817 ! | | | betan | | | 818 ! | | | |-----------| |------- 819 ! | | | | 820 ! | |---v------- | 821 ! | | 822 ! |-----------| 823 ! | 824 !-----------| 825 ! 826 !! 827 !!------------------------------------------------------------------ 828 829 REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & 830 zvolp, & ! total available pond water 831 zdvolp ! remaining meltwater after redistribution 832 833 INTEGER :: & 834 ns, & 835 m_index, & 836 permflag 837 838 REAL (wp), DIMENSION(jpl) :: & 839 hicen, & 840 hsnon, & 841 asnon, & 842 alfan, & 843 betan, & 844 cum_max_vol, & 845 reduced_aicen 846 847 REAL (wp), DIMENSION(0:jpl) :: & 848 cum_max_vol_tmp 849 850 REAL (wp) :: & 851 hpond, & 852 drain, & 853 floe_weight, & 854 pressure_head, & 855 hsl_rel, & 856 deltah, & 857 perm, & 858 msno 859 860 REAL (wp), parameter :: & 861 viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s 862 863 REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds 864 zvp_min = 1.e-4_wp ! minimum pond volume (m) 865 866 INTEGER :: ji, jj, jk, jl ! loop indices 867 868 a_ip(:,:,:) = 0._wp 869 h_ip(:,:,:) = 0._wp 870 871 DO_2D( 1, 1, 1, 1 ) 872 873 IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN 874 875 !------------------------------------------------------------------- 876 ! initialize 877 !------------------------------------------------------------------- 878 879 DO jl = 1, jpl 880 881 !---------------------------------------- 882 ! compute the effective snow fraction 883 !---------------------------------------- 884 885 IF (a_i(ji,jj,jl) < epsi10) THEN 886 hicen(jl) = 0._wp 887 hsnon(jl) = 0._wp 888 reduced_aicen(jl) = 0._wp 889 asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables 890 ELSE 891 hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) 892 hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) 893 reduced_aicen(jl) = 1._wp ! n=jpl 894 895 !js: initial code in NEMO_DEV 896 !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & 897 ! * (-0.024_wp*hicen(jl) + 0.832_wp) 898 899 !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large 900 IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & 901 * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) 902 903 asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) 904 ! MV should check whether this makes sense to have the same effective snow fraction in here 905 ! OLI: it probably doesn't 906 END IF 907 908 ! This choice for alfa and beta ignores hydrostatic equilibium of categories. 909 ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming 910 ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all 911 ! categories. alfa and beta partition the ITD - they are areas not thicknesses! 912 ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. 913 ! Here, alfa = 60% of the ice area (and since hice is constant in a category, 914 ! alfan = 60% of the ice volume) in each category lies above the reference line, 915 ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. 916 917 ! MV: 918 ! Note that this choice is not in the original FF07 paper and has been adopted in CICE 919 ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe 920 921 ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... 922 923 alfan(jl) = 0.6 * hicen(jl) 924 betan(jl) = 0.4 * hicen(jl) 925 926 cum_max_vol(jl) = 0._wp 927 cum_max_vol_tmp(jl) = 0._wp 928 929 END DO ! jpl 930 931 cum_max_vol_tmp(0) = 0._wp 932 drain = 0._wp 933 zdvolp(ji,jj) = 0._wp 934 935 !---------------------------------------------------------- 936 ! Drain overflow water, update pond fraction and volume 937 !---------------------------------------------------------- 938 939 !-------------------------------------------------------------------------- 940 ! the maximum amount of water that can be contained up to each ice category 941 !-------------------------------------------------------------------------- 942 ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) 943 ! Then the excess volume cum_max_vol(jl) drains out of the system 944 ! It should be added to wfx_pnd_out 945 946 DO jl = 1, jpl-1 ! last category can not hold any volume 947 948 IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN 949 950 ! total volume in level including snow 951 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & 952 (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) 953 954 ! subtract snow solid volumes from lower categories in current level 955 DO ns = 1, jl 956 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & 957 - rhos/rhow * & ! free air fraction that can be filled by water 958 asnon(ns) * & ! effective areal fraction of snow in that category 959 max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) 960 END DO 961 962 ELSE ! assume higher categories unoccupied 963 cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) 964 END IF 965 !IF (cum_max_vol_tmp(jl) < z0) THEN 966 ! CALL abort_ice('negative melt pond volume') 967 !END IF 968 END DO 969 cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume 970 cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) 971 972 !---------------------------------------------------------------- 973 ! is there more meltwater than can be held in the floe? 974 !---------------------------------------------------------------- 975 IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN 976 drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 977 zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available 978 979 diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) 980 981 zdvolp(ji,jj) = drain ! this is the drained water 982 IF (zvolp(ji,jj) < epsi10) THEN 983 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 984 zvolp(ji,jj) = 0._wp 985 END IF 986 END IF 987 988 ! height and area corresponding to the remaining volume 989 ! routine leaves zvolp unchanged 990 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 991 992 DO jl = 1, m_index 993 !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update 994 ! ! volume instead, no ? 995 h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 996 a_ip(ji,jj,jl) = reduced_aicen(jl) 997 ! in practise, pond fraction depends on the empirical snow fraction 998 ! so in turn on ice thickness 999 END DO 1000 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1001 1002 !------------------------------------------------------------------------ 1003 ! Drainage through brine network (permeability) 1004 !------------------------------------------------------------------------ 1005 !!! drainage due to ice permeability - Darcy's law 1006 1007 ! sea water level 1008 msno = 0._wp 1009 DO jl = 1 , jpl 1010 msno = msno + v_s(ji,jj,jl) * rhos 1011 END DO 1012 floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj) 1013 hsl_rel = floe_weight / rho0 & 1014 - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) 1015 1016 deltah = hpond - hsl_rel 1017 pressure_head = grav * rho0 * max(deltah, 0._wp) 1018 1019 ! drain if ice is permeable 1020 permflag = 0 1021 1022 IF (pressure_head > 0._wp) THEN 1023 DO jl = 1, jpl-1 1024 IF ( hicen(jl) /= 0._wp ) THEN 1025 1026 !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 1027 1028 perm = 0._wp ! MV ugly dummy patch 1029 CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof 1030 IF (perm > 0._wp) permflag = 1 1031 1032 drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / & 1033 (viscosity*hicen(jl)) 1034 zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) 1035 zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) 1036 1037 diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) 1038 1039 IF (zvolp(ji,jj) < epsi10) THEN 1040 zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) 1041 zvolp(ji,jj) = 0._wp 1042 END IF 1043 END IF 1044 END DO 1045 1046 ! adjust melt pond dimensions 1047 IF (permflag > 0) THEN 1048 ! recompute pond depth 1049 CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) 1050 DO jl = 1, m_index 1051 h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) 1052 a_ip(ji,jj,jl) = reduced_aicen(jl) 1053 END DO 1054 !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d 1055 END IF 1056 END IF ! pressure_head 1057 1058 !------------------------------- 1059 ! remove water from the snow 1060 !------------------------------- 1061 !------------------------------------------------------------------------ 1062 ! total melt pond volume in category does not include snow volume 1063 ! snow in melt ponds is not melted 1064 !------------------------------------------------------------------------ 1065 1066 ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell 1067 ! how much, so I did not diagnose it 1068 ! so if there is a problem here, nobody is going to see it... 1069 1070 1071 ! Calculate pond volume for lower categories 1072 DO jl = 1,m_index-1 1073 v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow 1074 - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) 1075 END DO 1076 1077 ! Calculate pond volume for highest category = remaining pond volume 1078 1079 ! The following is completely unclear to Martin at least 1080 ! Could we redefine properly and recode in a more readable way ? 1081 1082 ! m_index = last category with melt pond 1083 1084 IF (m_index == 1) v_ip(ji,jj,m_index) = zvolp(ji,jj) ! volume of mw in 1st category is the total volume of melt water 1085 1086 IF (m_index > 1) THEN 1087 IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN 1088 v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) 1089 ELSE 1090 v_ip(ji,jj,m_index) = 0._wp 1091 h_ip(ji,jj,m_index) = 0._wp 1092 a_ip(ji,jj,m_index) = 0._wp 1093 ! If remaining pond volume is negative reduce pond volume of 1094 ! lower category 1095 IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & 1096 v_ip(ji,jj,m_index-1) = v_ip(ji,jj,m_index-1) - sum(v_ip(ji,jj,1:m_index-1)) + zvolp(ji,jj) 1097 END IF 1098 END IF 1099 1100 DO jl = 1,m_index 1101 IF (a_ip(ji,jj,jl) > epsi10) THEN 1102 h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) 1103 ELSE 1104 zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) 1105 h_ip(ji,jj,jl) = 0._wp 1106 v_ip(ji,jj,jl) = 0._wp 1107 a_ip(ji,jj,jl) = 0._wp 1108 END IF 1109 END DO 1110 DO jl = m_index+1, jpl 1111 h_ip(ji,jj,jl) = 0._wp 1112 a_ip(ji,jj,jl) = 0._wp 1113 v_ip(ji,jj,jl) = 0._wp 1114 END DO 1115 1116 ENDIF 1117 1118 END_2D 1119 1120 END SUBROUTINE ice_thd_pnd_area 1121 1122 1123 SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) 1124 !!------------------------------------------------------------------- 1125 !! *** ROUTINE ice_thd_pnd_depth *** 1126 !! 1127 !! ** Purpose : Compute melt pond depth 1128 !!------------------------------------------------------------------- 1129 1130 REAL (wp), DIMENSION(jpl), INTENT(IN) :: & 1131 aicen, & 1132 asnon, & 1133 hsnon, & 1134 alfan, & 1135 cum_max_vol 1136 1137 REAL (wp), INTENT(IN) :: & 1138 zvolp 1139 1140 REAL (wp), INTENT(OUT) :: & 1141 hpond 1142 1143 INTEGER, INTENT(OUT) :: & 1144 m_index 1145 1146 INTEGER :: n, ns 1147 1148 REAL (wp), DIMENSION(0:jpl+1) :: & 1149 hitl, & 1150 aicetl 1151 1152 REAL (wp) :: & 1153 rem_vol, & 1154 area, & 1155 vol, & 1156 tmp, & 1157 z0 = 0.0_wp 1158 1159 !---------------------------------------------------------------- 1160 ! hpond is zero if zvolp is zero - have we fully drained? 1161 !---------------------------------------------------------------- 1162 1163 IF (zvolp < epsi10) THEN 1164 hpond = z0 1165 m_index = 0 1166 ELSE 1167 1168 !---------------------------------------------------------------- 1169 ! Calculate the category where water fills up to 1170 !---------------------------------------------------------------- 1171 1172 !----------| 1173 ! | 1174 ! | 1175 ! |----------| -- -- 1176 !__________|__________|_________________________________________ ^ 1177 ! | | rem_vol ^ | Semi-filled 1178 ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer 1179 ! | | | | 1180 ! | | | |hpond 1181 ! | | |----------| | |------- 1182 ! | | | | | | 1183 ! | | | |---v-----| 1184 ! | | m_index | | | 1185 !------------------------------------------------------------- 1186 1187 m_index = 0 ! 1:m_index categories have water in them 1188 DO n = 1, jpl 1189 IF (zvolp <= cum_max_vol(n)) THEN 1190 m_index = n 1191 IF (n == 1) THEN 1192 rem_vol = zvolp 1193 ELSE 1194 rem_vol = zvolp - cum_max_vol(n-1) 1195 END IF 1196 exit ! to break out of the loop 1197 END IF 1198 END DO 1199 m_index = min(jpl-1, m_index) 1200 1201 !---------------------------------------------------------------- 1202 ! semi-filled layer may have m_index different snow in it 1203 !---------------------------------------------------------------- 1204 1205 !----------------------------------------------------------- ^ 1206 ! | alfan(m_index+1) 1207 ! | 1208 !hitl(3)--> |----------| | 1209 !hitl(2)--> |------------| * * * * *| | 1210 !hitl(1)--> |----------|* * * * * * |* * * * * | | 1211 !hitl(0)-->------------------------------------------------- | ^ 1212 ! various snow from lower categories | |alfa(m_index) 1213 1214 ! hitl - heights of the snow layers from thinner and current categories 1215 ! aicetl - area of each snow depth in this layer 1216 1217 hitl(:) = z0 1218 aicetl(:) = z0 1219 DO n = 1, m_index 1220 hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & 1221 alfan(m_index+1) - alfan(m_index)), z0) 1222 aicetl(n) = asnon(n) 1223 1224 aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) 1225 END DO 1226 1227 hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) 1228 aicetl(m_index+1) = z0 1229 1230 !---------------------------------------------------------------- 1231 ! reorder array according to hitl 1232 ! snow heights not necessarily in height order 1233 !---------------------------------------------------------------- 1234 1235 DO ns = 1, m_index+1 1236 DO n = 0, m_index - ns + 1 1237 IF (hitl(n) > hitl(n+1)) THEN ! swap order 1238 tmp = hitl(n) 1239 hitl(n) = hitl(n+1) 1240 hitl(n+1) = tmp 1241 tmp = aicetl(n) 1242 aicetl(n) = aicetl(n+1) 1243 aicetl(n+1) = tmp 1244 END IF 1245 END DO 1246 END DO 1247 1248 !---------------------------------------------------------------- 1249 ! divide semi-filled layer into set of sublayers each vertically homogenous 1250 !---------------------------------------------------------------- 1251 1252 !hitl(3)---------------------------------------------------------------- 1253 ! | * * * * * * * * 1254 ! |* * * * * * * * * 1255 !hitl(2)---------------------------------------------------------------- 1256 ! | * * * * * * * * | * * * * * * * * 1257 ! |* * * * * * * * * |* * * * * * * * * 1258 !hitl(1)---------------------------------------------------------------- 1259 ! | * * * * * * * * | * * * * * * * * | * * * * * * * * 1260 ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * 1261 !hitl(0)---------------------------------------------------------------- 1262 ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) 1263 1264 ! move up over layers incrementing volume 1265 DO n = 1, m_index+1 1266 1267 area = sum(aicetl(:)) - & ! total area of sub-layer 1268 (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow 1269 1270 vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area 1271 1272 IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within 1273 hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) 1274 1275 exit 1276 ELSE ! still in sub-layer below the sub-layer with the depth 1277 rem_vol = rem_vol - vol 1278 END IF 1279 1280 END DO 1281 1282 END IF 1283 1284 END SUBROUTINE ice_thd_pnd_depth 1285 1286 1287 SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) 1288 !!------------------------------------------------------------------- 1289 !! *** ROUTINE ice_thd_pnd_perm *** 1290 !! 1291 !! ** Purpose : Determine the liquid fraction of brine in the ice 1292 !! and its permeability 1293 !!------------------------------------------------------------------- 1294 1295 REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & 1296 ticen, & ! internal ice temperature (K) 1297 salin ! salinity (ppt) !js: ppt according to cice 1298 1299 REAL (wp), INTENT(OUT) :: & 1300 perm ! permeability 1301 1302 REAL (wp) :: & 1303 Sbr ! brine salinity 1304 1305 REAL (wp), DIMENSION(nlay_i) :: & 1306 Tin, & ! ice temperature 1307 phi ! liquid fraction 1308 1309 INTEGER :: k 1310 1311 !----------------------------------------------------------------- 1312 ! Compute ice temperatures from enthalpies using quadratic formula 1313 !----------------------------------------------------------------- 1314 1315 DO k = 1,nlay_i 1316 Tin(k) = ticen(k) - rt0 !js: from K to degC 1317 END DO 1318 1319 !----------------------------------------------------------------- 1320 ! brine salinity and liquid fraction 1321 !----------------------------------------------------------------- 1322 1323 DO k = 1, nlay_i 1324 1325 Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) 1326 ! Best expression to date is that one (Vancoppenolle et al JGR 2019) 1327 ! Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 1328 phi(k) = salin(k) / Sbr 1329 1330 END DO 1331 1332 !----------------------------------------------------------------- 1333 ! permeability 1334 !----------------------------------------------------------------- 1335 1336 perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) 1337 1338 END SUBROUTINE ice_thd_pnd_perm 1339 1340 SUBROUTINE ice_thd_pnd_init 304 1341 !!------------------------------------------------------------------- 305 1342 !! *** ROUTINE ice_thd_pnd_init *** … … 308 1345 !! over sea ice 309 1346 !! 310 !! ** Method : Read the namthd_pnd namelist and check the melt pond 1347 !! ** Method : Read the namthd_pnd namelist and check the melt pond 311 1348 !! parameter values called at the first timestep (nit000) 312 1349 !! 313 !! ** input : Namelist namthd_pnd 1350 !! ** input : Namelist namthd_pnd 314 1351 !!------------------------------------------------------------------- 315 1352 INTEGER :: ios, ioptio ! Local integer 316 1353 !! 317 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, &1354 NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & 318 1355 & ln_pnd_CST , rn_apnd, rn_hpnd, & 1356 & ln_pnd_TOPO, & 319 1357 & ln_pnd_lids, ln_pnd_alb 320 1358 !!------------------------------------------------------------------- … … 332 1370 WRITE(numout,*) ' Namelist namicethd_pnd:' 333 1371 WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd 1372 WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO 334 1373 WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV 335 1374 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min 336 1375 WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max 1376 WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush 337 1377 WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST 338 1378 WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd … … 347 1387 IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF 348 1388 IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF 1389 IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF 349 1390 IF( ioptio /= 1 ) & 350 & 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)' )1391 & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV, ln_pnd_CST or ln_pnd_TOPO)' ) 351 1392 ! 352 1393 SELECT CASE( nice_pnd ) 353 CASE( np_pndNO ) 1394 CASE( np_pndNO ) 354 1395 IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF 355 1396 IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF 356 CASE( np_pndCST ) 1397 CASE( np_pndCST ) 357 1398 IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF 358 1399 END SELECT 359 1400 ! 360 1401 END SUBROUTINE ice_thd_pnd_init 361 1402 362 1403 #else 363 1404 !!---------------------------------------------------------------------- 364 1405 !! Default option Empty module NO SI3 sea-ice model 365 1406 !!---------------------------------------------------------------------- 366 #endif 1407 #endif 367 1408 368 1409 !!====================================================================== 369 END MODULE icethd_pnd 1410 END MODULE icethd_pnd
Note: See TracChangeset
for help on using the changeset viewer.
