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