[14072] | 1 | MODULE icethd_pnd |
---|
[8637] | 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 |
---|
[14005] | 22 | USE sbc_ice ! surface energy budget |
---|
[8637] | 23 | ! |
---|
| 24 | USE in_out_manager ! I/O manager |
---|
[14005] | 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: |
---|
[14005] | 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 | |
---|
[14072] | 43 | !-------------------------------------------------------------------------- |
---|
[14005] | 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 [-] |
---|
[14072] | 58 | REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: diag_dvpn_rnf ! meltwater pond lost to runoff [-] |
---|
[14005] | 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" |
---|
[8637] | 66 | !!---------------------------------------------------------------------- |
---|
[9598] | 67 | !! NEMO/ICE 4.0 , NEMO Consortium (2018) |
---|
[10069] | 68 | !! $Id$ |
---|
[10068] | 69 | !! Software governed by the CeCILL license (see ./LICENSE) |
---|
[8637] | 70 | !!---------------------------------------------------------------------- |
---|
| 71 | CONTAINS |
---|
| 72 | |
---|
| 73 | SUBROUTINE ice_thd_pnd |
---|
[14005] | 74 | |
---|
[8637] | 75 | !!------------------------------------------------------------------- |
---|
| 76 | !! *** ROUTINE ice_thd_pnd *** |
---|
[14072] | 77 | !! |
---|
[13472] | 78 | !! ** Purpose : change melt pond fraction and thickness |
---|
[14005] | 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 |
---|
[8637] | 83 | !!------------------------------------------------------------------- |
---|
[14005] | 84 | INTEGER :: ji, jj, jl ! loop indices |
---|
| 85 | !!------------------------------------------------------------------- |
---|
[14072] | 86 | |
---|
[14005] | 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) ) |
---|
[9169] | 89 | ! |
---|
[14005] | 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 ) |
---|
[9169] | 99 | ! |
---|
[14005] | 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 |
---|
[14072] | 113 | |
---|
[14005] | 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 |
---|
[9169] | 147 | ! |
---|
[14005] | 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 ) |
---|
[8637] | 150 | |
---|
[14072] | 151 | END SUBROUTINE ice_thd_pnd |
---|
[9169] | 152 | |
---|
[14072] | 153 | |
---|
| 154 | SUBROUTINE pnd_CST |
---|
[8637] | 155 | !!------------------------------------------------------------------- |
---|
| 156 | !! *** ROUTINE pnd_CST *** |
---|
| 157 | !! |
---|
[9169] | 158 | !! ** Purpose : Compute melt pond evolution |
---|
[8637] | 159 | !! |
---|
[14072] | 160 | !! ** Method : Melt pond fraction and thickness are prescribed |
---|
[9604] | 161 | !! to non-zero values when t_su = 0C |
---|
[8637] | 162 | !! |
---|
| 163 | !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd) |
---|
[14072] | 164 | !! |
---|
[9169] | 165 | !! ** Note : Coupling with such melt ponds is only radiative |
---|
[9604] | 166 | !! Advection, ridging, rafting... are bypassed |
---|
[8637] | 167 | !! |
---|
| 168 | !! ** References : Bush, G.W., and Trump, D.J. (2017) |
---|
| 169 | !!------------------------------------------------------------------- |
---|
[14005] | 170 | INTEGER :: ji, jl ! loop indices |
---|
| 171 | REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids |
---|
[8637] | 172 | !!------------------------------------------------------------------- |
---|
[14005] | 173 | DO jl = 1, jpl |
---|
[14072] | 174 | |
---|
[14005] | 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 |
---|
[14072] | 187 | h_ip_1d(ji) = rn_hpnd |
---|
[14005] | 188 | a_ip_1d(ji) = rn_apnd * a_i_1d(ji) |
---|
| 189 | h_il_1d(ji) = 0._wp ! no pond lids whatsoever |
---|
| 190 | ELSE |
---|
[14072] | 191 | h_ip_1d(ji) = 0._wp |
---|
[14005] | 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 | |
---|
[8637] | 211 | END DO |
---|
[9169] | 212 | ! |
---|
[8637] | 213 | END SUBROUTINE pnd_CST |
---|
| 214 | |
---|
[9169] | 215 | |
---|
[13472] | 216 | SUBROUTINE pnd_LEV |
---|
[8637] | 217 | !!------------------------------------------------------------------- |
---|
[13472] | 218 | !! *** ROUTINE pnd_LEV *** |
---|
[8637] | 219 | !! |
---|
[13472] | 220 | !! ** Purpose : Compute melt pond evolution |
---|
[8637] | 221 | !! |
---|
[13472] | 222 | !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing |
---|
| 223 | !! We work with volumes and then redistribute changes into thickness and concentration |
---|
[14072] | 224 | !! assuming linear relationship between the two. |
---|
[8637] | 225 | !! |
---|
[13472] | 226 | !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- |
---|
| 227 | !! dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i |
---|
| 228 | !! dh_i = meltwater from ice surface melt |
---|
| 229 | !! dh_s = meltwater from snow melt |
---|
| 230 | !! (1-r) = fraction of melt water that is not flushed |
---|
| 231 | !! |
---|
| 232 | !! - limtations: a_ip must not exceed (1-r)*a_i |
---|
| 233 | !! h_ip must not exceed 0.5*h_i |
---|
| 234 | !! |
---|
| 235 | !! - pond shrinking: |
---|
| 236 | !! if lids: Vp = Vp -dH * a_ip |
---|
| 237 | !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- |
---|
| 238 | !! |
---|
[14072] | 239 | !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H |
---|
[13472] | 240 | !! H = lid thickness |
---|
| 241 | !! Lf = latent heat of fusion |
---|
| 242 | !! Tp = -2C |
---|
| 243 | !! |
---|
| 244 | !! And solved implicitely as: |
---|
| 245 | !! H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0 |
---|
| 246 | !! |
---|
| 247 | !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- |
---|
| 248 | !! |
---|
[14005] | 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) |
---|
[13472] | 251 | !! visc = water viscosity |
---|
| 252 | !! Hp = height of top of the pond above sea-level |
---|
| 253 | !! Hi = ice thickness thru which there is flushing |
---|
[14005] | 254 | !! flush= correction otherwise flushing is excessive |
---|
[13472] | 255 | !! |
---|
| 256 | !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness |
---|
| 257 | !! |
---|
| 258 | !! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip: |
---|
| 259 | !! a_ip/a_i = a_ip_frac = h_ip / zaspect |
---|
| 260 | !! |
---|
[14005] | 261 | !! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush |
---|
[8637] | 262 | !! |
---|
[14072] | 263 | !! ** Note : Mostly stolen from CICE but not only. These are between level-ice ponds and CESM ponds. |
---|
| 264 | !! |
---|
[13472] | 265 | !! ** References : Flocco and Feltham (JGR, 2007) |
---|
| 266 | !! Flocco et al (JGR, 2010) |
---|
| 267 | !! Holland et al (J. Clim, 2012) |
---|
[14005] | 268 | !! Hunke et al (OM 2012) |
---|
[14072] | 269 | !!------------------------------------------------------------------- |
---|
[13472] | 270 | REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array |
---|
| 271 | !! |
---|
| 272 | REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio |
---|
| 273 | REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature |
---|
| 274 | REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity |
---|
| 275 | !! |
---|
[14005] | 276 | REAL(wp) :: zfr_mlt, zdv_mlt, zdv_avail ! fraction and volume of available meltwater retained for melt ponding |
---|
[13472] | 277 | REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes |
---|
[14005] | 278 | REAL(wp) :: zdv_pnd ! Amount of water going into the ponds & lids |
---|
[13472] | 279 | REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh |
---|
| 280 | REAL(wp) :: zv_ip_max ! max pond volume allowed |
---|
| 281 | REAL(wp) :: zdT ! zTp-t_su |
---|
[14005] | 282 | REAL(wp) :: zsbr, ztmelts ! Brine salinity |
---|
[13472] | 283 | REAL(wp) :: zperm ! permeability of sea ice |
---|
| 284 | REAL(wp) :: zfac, zdum ! temporary arrays |
---|
| 285 | REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse |
---|
| 286 | !! |
---|
[14005] | 287 | INTEGER :: ji, jk, jl ! loop indices |
---|
[8637] | 288 | !!------------------------------------------------------------------- |
---|
[14072] | 289 | z1_rhow = 1._wp / rhow |
---|
[13472] | 290 | z1_aspect = 1._wp / zaspect |
---|
[14072] | 291 | z1_Tp = 1._wp / zTp |
---|
| 292 | |
---|
[14005] | 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 ) |
---|
[14072] | 295 | |
---|
[14005] | 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 ) |
---|
[8637] | 300 | |
---|
[14005] | 301 | DO jl = 1, jpl |
---|
[13472] | 302 | |
---|
[14005] | 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 |
---|
[14072] | 317 | |
---|
[14005] | 318 | !----------------------- |
---|
| 319 | ! Melt pond calculations |
---|
| 320 | !----------------------- |
---|
| 321 | DO ji = 1, npti |
---|
[13472] | 322 | ! |
---|
[14005] | 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) |
---|
[9169] | 336 | ! |
---|
[14005] | 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 |
---|
[14072] | 344 | zdv_mlt = MAX( 0._wp, zfr_mlt * zdv_avail ) ! max for roundoff errors? |
---|
[14005] | 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 |
---|
[14072] | 351 | ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: |
---|
[14005] | 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) ) ) |
---|
[13472] | 354 | |
---|
[14005] | 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 |
---|
[14072] | 358 | ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: |
---|
[14005] | 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 |
---|
[13472] | 364 | ! |
---|
[14005] | 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 | ! |
---|
[14072] | 377 | !--- Lid growing and subsequent pond shrinking ---! |
---|
[14005] | 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 |
---|
[13472] | 380 | |
---|
[14005] | 381 | ! Lid growing |
---|
| 382 | v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_frz ) |
---|
[13472] | 383 | |
---|
[14005] | 384 | ! Pond shrinking |
---|
| 385 | v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) |
---|
[13472] | 386 | |
---|
[14005] | 387 | ELSE |
---|
[14072] | 388 | zdv_frz = v_ip_1d(ji) * ( EXP( 0.01_wp * zdT * z1_Tp ) - 1._wp ) ! Holland 2012 (eq. 6) |
---|
[14005] | 389 | ! Pond shrinking |
---|
| 390 | v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zdv_frz ) |
---|
[13472] | 391 | ENDIF |
---|
[14005] | 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 | |
---|
[14072] | 400 | !------------------------------------------------! |
---|
[14005] | 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) |
---|
[14072] | 422 | zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) ! < 0 |
---|
[14005] | 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 |
---|
[13472] | 441 | ENDIF |
---|
[14005] | 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 | ! |
---|
[13472] | 449 | ENDIF |
---|
[9169] | 450 | ! |
---|
[14005] | 451 | v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) |
---|
| 452 | v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) |
---|
| 453 | ! |
---|
| 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 |
---|
| 456 | ! |
---|
| 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 | ! |
---|
[8637] | 468 | END DO |
---|
[9169] | 469 | ! |
---|
[14005] | 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 | ! |
---|
[13472] | 477 | END SUBROUTINE pnd_LEV |
---|
[8637] | 478 | |
---|
[9169] | 479 | |
---|
[14005] | 480 | |
---|
[14072] | 481 | SUBROUTINE pnd_TOPO |
---|
| 482 | |
---|
[14005] | 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 |
---|
[14072] | 490 | !! (2007) and Flocco et al. (2010). |
---|
[14005] | 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 = 0.15_wp , & ! temperature difference for freeze-up (C) |
---|
| 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, & !! total melt pond water available before redistribution and drainage |
---|
| 533 | zvolp_res !! remaining melt pond water available after drainage |
---|
[14072] | 534 | |
---|
[14005] | 535 | REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i |
---|
| 536 | |
---|
| 537 | INTEGER :: ji, jj, jk, jl ! loop indices |
---|
| 538 | |
---|
| 539 | INTEGER :: i_test |
---|
| 540 | |
---|
| 541 | ! Note |
---|
| 542 | ! equivalent for CICE translation |
---|
| 543 | ! a_ip -> apond |
---|
| 544 | ! a_ip_frac -> apnd |
---|
| 545 | |
---|
| 546 | CALL ctl_stop( 'STOP', 'icethd_pnd : topographic melt ponds are still an ongoing work' ) |
---|
[14072] | 547 | |
---|
[14005] | 548 | !--------------------------------------------------------------- |
---|
| 549 | ! Initialise |
---|
| 550 | !--------------------------------------------------------------- |
---|
| 551 | |
---|
| 552 | ! Parameters & constants (move to parameters) |
---|
| 553 | zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3) |
---|
| 554 | zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline) |
---|
[14072] | 555 | z1_rhow = 1._wp / rhow |
---|
[14005] | 556 | |
---|
| 557 | ! Set required ice variables (hard-coded here for now) |
---|
[14072] | 558 | ! zfpond(:,:) = 0._wp ! contributing freshwater flux (?) |
---|
| 559 | |
---|
[14005] | 560 | at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction |
---|
| 561 | vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area |
---|
| 562 | vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area |
---|
| 563 | vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area |
---|
[14072] | 564 | |
---|
[14005] | 565 | ! thickness |
---|
| 566 | WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:) |
---|
| 567 | ELSEWHERE ; z1_a_i(:,:,:) = 0._wp |
---|
| 568 | END WHERE |
---|
| 569 | h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:) |
---|
[14072] | 570 | |
---|
[14005] | 571 | !--------------------------------------------------------------- |
---|
| 572 | ! Change 2D to 1D |
---|
| 573 | !--------------------------------------------------------------- |
---|
[14072] | 574 | ! MV |
---|
[14005] | 575 | ! a less computing-intensive version would have 2D-1D passage here |
---|
| 576 | ! use what we have in iceitd.F90 (incremental remapping) |
---|
| 577 | |
---|
| 578 | !-------------------------------------------------------------- |
---|
| 579 | ! Collect total available pond water volume |
---|
| 580 | !-------------------------------------------------------------- |
---|
| 581 | ! Assuming that meltwater (+rain in principle) runsoff the surface |
---|
| 582 | ! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction |
---|
| 583 | ! I cite her words, they are very talkative |
---|
[14072] | 584 | ! "grid cells with very little ice cover (and hence more open water area) |
---|
[14005] | 585 | ! have a higher runoff fraction to rep- resent the greater proximity of ice to open water." |
---|
| 586 | ! "This results in the same runoff fraction r for each ice category within a grid cell" |
---|
[14072] | 587 | |
---|
[14005] | 588 | zvolp(:,:) = 0._wp |
---|
| 589 | |
---|
| 590 | DO jl = 1, jpl |
---|
| 591 | DO_2D( 1, 1, 1, 1 ) |
---|
[14072] | 592 | |
---|
[14005] | 593 | IF ( a_i(ji,jj,jl) > epsi10 ) THEN |
---|
[14072] | 594 | |
---|
[14005] | 595 | !--- Available and contributing meltwater for melt ponding ---! |
---|
| 596 | 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 |
---|
| 597 | & * z1_rhow * a_i(ji,jj,jl) |
---|
| 598 | ! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area |
---|
| 599 | zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds |
---|
| 600 | zv_pnd = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl |
---|
| 601 | |
---|
| 602 | diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_Dt_ice ! diags |
---|
[14072] | 603 | diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_Dt_ice |
---|
[14005] | 604 | |
---|
| 605 | !--- Create possible new ponds |
---|
| 606 | ! if pond does not exist, create new pond over full ice area |
---|
| 607 | !IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN |
---|
| 608 | IF ( a_ip(ji,jj,jl) < epsi10 ) THEN |
---|
| 609 | a_ip(ji,jj,jl) = a_i(ji,jj,jl) |
---|
| 610 | a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE) |
---|
| 611 | ENDIF |
---|
[14072] | 612 | |
---|
[14005] | 613 | !--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage |
---|
| 614 | v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd ! use pond water to increase thickness |
---|
| 615 | h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) |
---|
[14072] | 616 | |
---|
[14005] | 617 | !--- Total available pond water volume (pre-existing + newly produced)j |
---|
[14072] | 618 | zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl) |
---|
[14005] | 619 | ! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now |
---|
[14072] | 620 | |
---|
[14005] | 621 | ENDIF ! a_i |
---|
| 622 | |
---|
| 623 | END_2D |
---|
| 624 | END DO ! ji |
---|
[14072] | 625 | |
---|
[14005] | 626 | !-------------------------------------------------------------- |
---|
| 627 | ! Redistribute and drain water from ponds |
---|
[14072] | 628 | !-------------------------------------------------------------- |
---|
[14005] | 629 | CALL ice_thd_pnd_area( zvolp, zvolp_res ) |
---|
[14072] | 630 | |
---|
[14005] | 631 | !-------------------------------------------------------------- |
---|
| 632 | ! Melt pond lid growth and melt |
---|
[14072] | 633 | !-------------------------------------------------------------- |
---|
| 634 | |
---|
[14005] | 635 | IF( ln_pnd_lids ) THEN |
---|
| 636 | |
---|
| 637 | DO_2D( 1, 1, 1, 1 ) |
---|
| 638 | |
---|
| 639 | IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > rn_himin .AND. vt_ip(ji,jj) > zvp_min * at_i(ji,jj) ) THEN |
---|
[14072] | 640 | |
---|
[14005] | 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 |
---|
[14072] | 650 | |
---|
[14005] | 651 | DO jl = 1, jpl-1 |
---|
[14072] | 652 | |
---|
[14005] | 653 | IF ( v_il(ji,jj,jl) > epsi10 ) THEN |
---|
[14072] | 654 | |
---|
[14005] | 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) ) |
---|
[14072] | 662 | |
---|
[14005] | 663 | IF ( zdvice > epsi10 ) THEN |
---|
[14072] | 664 | |
---|
[14005] | 665 | v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice |
---|
[14072] | 666 | v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zdvice ! MV: not sure i understand dh_i_sum seems counted twice - |
---|
[14005] | 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) |
---|
[14072] | 670 | |
---|
[14005] | 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 |
---|
[14072] | 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) |
---|
[14005] | 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 |
---|
[14072] | 678 | |
---|
[14005] | 679 | diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag |
---|
[14072] | 680 | |
---|
[14005] | 681 | ENDIF |
---|
[14072] | 682 | |
---|
[14005] | 683 | !---------------------------------------------------------------- |
---|
[14072] | 684 | ! Freeze pre-existing lid |
---|
[14005] | 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 |
---|
[14072] | 690 | zdTice = MAX( - t_su(ji,jj,jl) - zTd , 0._wp ) ! > 0 |
---|
[14005] | 691 | zomega = rcnd_i * zdTice / zrhoi_L |
---|
| 692 | zdHui = SQRT( 2._wp * zomega * rDt_ice + ( v_il(ji,jj,jl) / a_i(ji,jj,jl) )**2 ) & |
---|
| 693 | - v_il(ji,jj,jl) / a_i(ji,jj,jl) |
---|
| 694 | zdvice = min( zdHui*a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) |
---|
[14072] | 695 | |
---|
[14005] | 696 | IF ( zdvice > epsi10 ) THEN |
---|
| 697 | v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice |
---|
| 698 | v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice |
---|
| 699 | ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice |
---|
| 700 | ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice |
---|
| 701 | h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) |
---|
[14072] | 702 | |
---|
[14005] | 703 | diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag |
---|
[14072] | 704 | |
---|
[14005] | 705 | ENDIF |
---|
[14072] | 706 | |
---|
[14005] | 707 | ENDIF ! Tsfcn(i,j,n) |
---|
| 708 | |
---|
| 709 | !---------------------------------------------------------------- |
---|
| 710 | ! Freeze new lids |
---|
| 711 | !---------------------------------------------------------------- |
---|
| 712 | ! upper ice layer begins to form |
---|
| 713 | ! note: albedo does not change |
---|
| 714 | |
---|
| 715 | ELSE ! v_il < epsi10 |
---|
[14072] | 716 | |
---|
[14005] | 717 | ! thickness of newly formed ice |
---|
| 718 | ! the surface temperature of a meltpond is the same as that |
---|
[14072] | 719 | ! of the ice underneath (0C), and the thermodynamic surface |
---|
[14005] | 720 | ! flux is the same |
---|
[14072] | 721 | |
---|
[14005] | 722 | !!! we need net surface energy flux, excluding conduction |
---|
| 723 | !!! fsurf is summed over categories in CICE |
---|
| 724 | !!! we have the category-dependent flux, let us use it ? |
---|
[14072] | 725 | zfsurf = qns_ice(ji,jj,jl) + qsr_ice(ji,jj,jl) |
---|
[14005] | 726 | zdHui = MAX ( -zfsurf * rDt_ice/zrhoi_L , 0._wp ) |
---|
| 727 | zdvice = MIN ( zdHui * a_ip(ji,jj,jl) , v_ip(ji,jj,jl) ) |
---|
| 728 | IF ( zdvice > epsi10 ) THEN |
---|
| 729 | v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice |
---|
| 730 | v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice |
---|
[14072] | 731 | |
---|
[14005] | 732 | diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag |
---|
| 733 | ! zvolp(ji,jj) = zvolp(ji,jj) - zdvice |
---|
| 734 | ! zfpond(ji,jj) = zfpond(ji,jj) - zdvice |
---|
| 735 | 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 |
---|
| 736 | ENDIF |
---|
[14072] | 737 | |
---|
[14005] | 738 | ENDIF ! v_il |
---|
[14072] | 739 | |
---|
[14005] | 740 | END DO ! jl |
---|
| 741 | |
---|
| 742 | ELSE ! remove ponds on thin ice |
---|
| 743 | |
---|
| 744 | v_ip(ji,jj,:) = 0._wp |
---|
| 745 | v_il(ji,jj,:) = 0._wp |
---|
| 746 | ! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj) |
---|
[14072] | 747 | ! zvolp(ji,jj) = 0._wp |
---|
[14005] | 748 | |
---|
| 749 | ENDIF |
---|
| 750 | |
---|
| 751 | END_2D |
---|
| 752 | |
---|
| 753 | ENDIF ! ln_pnd_lids |
---|
| 754 | |
---|
| 755 | !--------------------------------------------------------------- |
---|
| 756 | ! Clean-up variables (probably duplicates what icevar would do) |
---|
| 757 | !--------------------------------------------------------------- |
---|
| 758 | ! MV comment |
---|
| 759 | ! In the ideal world, the lines above should update only v_ip, a_ip, v_il |
---|
| 760 | ! icevar should recompute all other variables (if needed at all) |
---|
| 761 | |
---|
| 762 | DO jl = 1, jpl |
---|
| 763 | |
---|
| 764 | DO_2D( 1, 1, 1, 1 ) |
---|
| 765 | |
---|
| 766 | ! ! zap lids on small ponds |
---|
| 767 | ! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 & |
---|
| 768 | ! .AND. v_il(ji,jj,jl) > epsi10) THEN |
---|
| 769 | ! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small |
---|
| 770 | ! ENDIF |
---|
[14072] | 771 | |
---|
[14005] | 772 | ! recalculate equivalent pond variables |
---|
| 773 | IF ( a_ip(ji,jj,jl) > epsi10) THEN |
---|
| 774 | h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_i(ji,jj,jl) |
---|
| 775 | a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar |
---|
| 776 | h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar |
---|
| 777 | ENDIF |
---|
| 778 | ! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar |
---|
| 779 | ! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar |
---|
| 780 | ! ENDIF |
---|
[14072] | 781 | |
---|
[14005] | 782 | END_2D |
---|
| 783 | |
---|
| 784 | END DO ! jl |
---|
| 785 | |
---|
| 786 | |
---|
| 787 | END SUBROUTINE pnd_TOPO |
---|
| 788 | |
---|
[14072] | 789 | |
---|
[14005] | 790 | SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp ) |
---|
| 791 | |
---|
| 792 | !!------------------------------------------------------------------- |
---|
| 793 | !! *** ROUTINE ice_thd_pnd_area *** |
---|
| 794 | !! |
---|
[14072] | 795 | !! ** Purpose : Given the total volume of available pond water, |
---|
[14005] | 796 | !! redistribute and drain water |
---|
| 797 | !! |
---|
| 798 | !! ** Method |
---|
| 799 | !! |
---|
| 800 | !-----------| |
---|
| 801 | ! | |
---|
| 802 | ! |-----------| |
---|
| 803 | !___________|___________|______________________________________sea-level |
---|
| 804 | ! | | |
---|
| 805 | ! | |---^--------| |
---|
| 806 | ! | | | | |
---|
| 807 | ! | | | |-----------| |------- |
---|
| 808 | ! | | | alfan | | | |
---|
| 809 | ! | | | | |--------------| |
---|
| 810 | ! | | | | | | |
---|
| 811 | !---------------------------v------------------------------------------- |
---|
| 812 | ! | | ^ | | | |
---|
| 813 | ! | | | | |--------------| |
---|
| 814 | ! | | | betan | | | |
---|
| 815 | ! | | | |-----------| |------- |
---|
| 816 | ! | | | | |
---|
| 817 | ! | |---v------- | |
---|
| 818 | ! | | |
---|
| 819 | ! |-----------| |
---|
| 820 | ! | |
---|
| 821 | !-----------| |
---|
| 822 | ! |
---|
| 823 | !! |
---|
| 824 | !!------------------------------------------------------------------ |
---|
[14072] | 825 | |
---|
[14005] | 826 | REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: & |
---|
| 827 | zvolp, & ! total available pond water |
---|
| 828 | zdvolp ! remaining meltwater after redistribution |
---|
| 829 | |
---|
| 830 | INTEGER :: & |
---|
| 831 | ns, & |
---|
| 832 | m_index, & |
---|
| 833 | permflag |
---|
| 834 | |
---|
| 835 | REAL (wp), DIMENSION(jpl) :: & |
---|
| 836 | hicen, & |
---|
| 837 | hsnon, & |
---|
| 838 | asnon, & |
---|
| 839 | alfan, & |
---|
| 840 | betan, & |
---|
| 841 | cum_max_vol, & |
---|
| 842 | reduced_aicen |
---|
| 843 | |
---|
| 844 | REAL (wp), DIMENSION(0:jpl) :: & |
---|
| 845 | cum_max_vol_tmp |
---|
| 846 | |
---|
| 847 | REAL (wp) :: & |
---|
| 848 | hpond, & |
---|
| 849 | drain, & |
---|
| 850 | floe_weight, & |
---|
| 851 | pressure_head, & |
---|
| 852 | hsl_rel, & |
---|
| 853 | deltah, & |
---|
| 854 | perm, & |
---|
| 855 | msno |
---|
| 856 | |
---|
| 857 | REAL (wp), parameter :: & |
---|
| 858 | viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s |
---|
| 859 | |
---|
| 860 | REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds |
---|
| 861 | zvp_min = 1.e-4_wp ! minimum pond volume (m) |
---|
| 862 | |
---|
| 863 | INTEGER :: ji, jj, jk, jl ! loop indices |
---|
| 864 | |
---|
| 865 | a_ip(:,:,:) = 0._wp |
---|
| 866 | h_ip(:,:,:) = 0._wp |
---|
[14072] | 867 | |
---|
[14005] | 868 | DO_2D( 1, 1, 1, 1 ) |
---|
[14072] | 869 | |
---|
[14005] | 870 | 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 |
---|
[14072] | 871 | |
---|
[14005] | 872 | !------------------------------------------------------------------- |
---|
| 873 | ! initialize |
---|
| 874 | !------------------------------------------------------------------- |
---|
[14072] | 875 | |
---|
[14005] | 876 | DO jl = 1, jpl |
---|
[14072] | 877 | |
---|
[14005] | 878 | !---------------------------------------- |
---|
| 879 | ! compute the effective snow fraction |
---|
| 880 | !---------------------------------------- |
---|
[14072] | 881 | |
---|
[14005] | 882 | IF (a_i(ji,jj,jl) < epsi10) THEN |
---|
| 883 | hicen(jl) = 0._wp |
---|
| 884 | hsnon(jl) = 0._wp |
---|
| 885 | reduced_aicen(jl) = 0._wp |
---|
| 886 | asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables |
---|
| 887 | ELSE |
---|
| 888 | hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl) |
---|
| 889 | hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl) |
---|
| 890 | reduced_aicen(jl) = 1._wp ! n=jpl |
---|
[14072] | 891 | |
---|
[14005] | 892 | !js: initial code in NEMO_DEV |
---|
| 893 | !IF (n < jpl) reduced_aicen(jl) = aicen(jl) & |
---|
| 894 | ! * (-0.024_wp*hicen(jl) + 0.832_wp) |
---|
[14072] | 895 | |
---|
[14005] | 896 | !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large |
---|
[14072] | 897 | IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) & |
---|
[14005] | 898 | * max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp)) |
---|
[14072] | 899 | |
---|
[14005] | 900 | asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical) |
---|
| 901 | ! MV should check whether this makes sense to have the same effective snow fraction in here |
---|
| 902 | ! OLI: it probably doesn't |
---|
| 903 | END IF |
---|
[14072] | 904 | |
---|
[14005] | 905 | ! This choice for alfa and beta ignores hydrostatic equilibium of categories. |
---|
| 906 | ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming |
---|
| 907 | ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all |
---|
| 908 | ! categories. alfa and beta partition the ITD - they are areas not thicknesses! |
---|
| 909 | ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. |
---|
| 910 | ! Here, alfa = 60% of the ice area (and since hice is constant in a category, |
---|
| 911 | ! alfan = 60% of the ice volume) in each category lies above the reference line, |
---|
| 912 | ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. |
---|
[14072] | 913 | |
---|
[14005] | 914 | ! MV: |
---|
| 915 | ! Note that this choice is not in the original FF07 paper and has been adopted in CICE |
---|
| 916 | ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe |
---|
[14072] | 917 | |
---|
[14005] | 918 | ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... |
---|
[14072] | 919 | |
---|
[14005] | 920 | alfan(jl) = 0.6 * hicen(jl) |
---|
| 921 | betan(jl) = 0.4 * hicen(jl) |
---|
[14072] | 922 | |
---|
[14005] | 923 | cum_max_vol(jl) = 0._wp |
---|
| 924 | cum_max_vol_tmp(jl) = 0._wp |
---|
[14072] | 925 | |
---|
[14005] | 926 | END DO ! jpl |
---|
[14072] | 927 | |
---|
[14005] | 928 | cum_max_vol_tmp(0) = 0._wp |
---|
| 929 | drain = 0._wp |
---|
| 930 | zdvolp(ji,jj) = 0._wp |
---|
[14072] | 931 | |
---|
[14005] | 932 | !---------------------------------------------------------- |
---|
| 933 | ! Drain overflow water, update pond fraction and volume |
---|
| 934 | !---------------------------------------------------------- |
---|
[14072] | 935 | |
---|
[14005] | 936 | !-------------------------------------------------------------------------- |
---|
| 937 | ! the maximum amount of water that can be contained up to each ice category |
---|
| 938 | !-------------------------------------------------------------------------- |
---|
| 939 | ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) |
---|
| 940 | ! Then the excess volume cum_max_vol(jl) drains out of the system |
---|
| 941 | ! It should be added to wfx_pnd_out |
---|
[14072] | 942 | |
---|
[14005] | 943 | DO jl = 1, jpl-1 ! last category can not hold any volume |
---|
[14072] | 944 | |
---|
[14005] | 945 | IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN |
---|
[14072] | 946 | |
---|
[14005] | 947 | ! total volume in level including snow |
---|
| 948 | cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + & |
---|
| 949 | (alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl)) |
---|
[14072] | 950 | |
---|
[14005] | 951 | ! subtract snow solid volumes from lower categories in current level |
---|
| 952 | DO ns = 1, jl |
---|
| 953 | cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) & |
---|
| 954 | - rhos/rhow * & ! free air fraction that can be filled by water |
---|
| 955 | asnon(ns) * & ! effective areal fraction of snow in that category |
---|
| 956 | max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp) |
---|
| 957 | END DO |
---|
[14072] | 958 | |
---|
[14005] | 959 | ELSE ! assume higher categories unoccupied |
---|
| 960 | cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) |
---|
| 961 | END IF |
---|
| 962 | !IF (cum_max_vol_tmp(jl) < z0) THEN |
---|
| 963 | ! CALL abort_ice('negative melt pond volume') |
---|
| 964 | !END IF |
---|
| 965 | END DO |
---|
| 966 | cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume |
---|
| 967 | cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) |
---|
[14072] | 968 | |
---|
[14005] | 969 | !---------------------------------------------------------------- |
---|
| 970 | ! is there more meltwater than can be held in the floe? |
---|
| 971 | !---------------------------------------------------------------- |
---|
| 972 | IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN |
---|
| 973 | drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10 |
---|
| 974 | zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available |
---|
[14072] | 975 | |
---|
[14005] | 976 | diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice) |
---|
[14072] | 977 | |
---|
[14005] | 978 | zdvolp(ji,jj) = drain ! this is the drained water |
---|
| 979 | IF (zvolp(ji,jj) < epsi10) THEN |
---|
| 980 | zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) |
---|
| 981 | zvolp(ji,jj) = 0._wp |
---|
| 982 | END IF |
---|
| 983 | END IF |
---|
[14072] | 984 | |
---|
[14005] | 985 | ! height and area corresponding to the remaining volume |
---|
| 986 | ! routine leaves zvolp unchanged |
---|
| 987 | CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) |
---|
[14072] | 988 | |
---|
[14005] | 989 | DO jl = 1, m_index |
---|
| 990 | !h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update |
---|
| 991 | ! ! volume instead, no ? |
---|
| 992 | h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2 |
---|
| 993 | a_ip(ji,jj,jl) = reduced_aicen(jl) |
---|
| 994 | ! in practise, pond fraction depends on the empirical snow fraction |
---|
| 995 | ! so in turn on ice thickness |
---|
| 996 | END DO |
---|
| 997 | !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d |
---|
[14072] | 998 | |
---|
[14005] | 999 | !------------------------------------------------------------------------ |
---|
| 1000 | ! Drainage through brine network (permeability) |
---|
| 1001 | !------------------------------------------------------------------------ |
---|
| 1002 | !!! drainage due to ice permeability - Darcy's law |
---|
[14072] | 1003 | |
---|
[14005] | 1004 | ! sea water level |
---|
[14072] | 1005 | msno = 0._wp |
---|
[14005] | 1006 | DO jl = 1 , jpl |
---|
| 1007 | msno = msno + v_s(ji,jj,jl) * rhos |
---|
| 1008 | END DO |
---|
| 1009 | floe_weight = ( msno + rhoi*vt_i(ji,jj) + rho0*zvolp(ji,jj) ) / at_i(ji,jj) |
---|
| 1010 | hsl_rel = floe_weight / rho0 & |
---|
| 1011 | - ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) ) |
---|
[14072] | 1012 | |
---|
[14005] | 1013 | deltah = hpond - hsl_rel |
---|
| 1014 | pressure_head = grav * rho0 * max(deltah, 0._wp) |
---|
[14072] | 1015 | |
---|
[14005] | 1016 | ! drain if ice is permeable |
---|
| 1017 | permflag = 0 |
---|
[14072] | 1018 | |
---|
[14005] | 1019 | IF (pressure_head > 0._wp) THEN |
---|
| 1020 | DO jl = 1, jpl-1 |
---|
| 1021 | IF ( hicen(jl) /= 0._wp ) THEN |
---|
[14072] | 1022 | |
---|
[14005] | 1023 | !IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2 |
---|
[14072] | 1024 | |
---|
[14005] | 1025 | perm = 0._wp ! MV ugly dummy patch |
---|
[14072] | 1026 | CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof |
---|
[14005] | 1027 | IF (perm > 0._wp) permflag = 1 |
---|
[14072] | 1028 | |
---|
[14005] | 1029 | drain = perm*a_ip(ji,jj,jl)*pressure_head*rDt_ice / & |
---|
| 1030 | (viscosity*hicen(jl)) |
---|
| 1031 | zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj)) |
---|
| 1032 | zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp) |
---|
[14072] | 1033 | |
---|
[14005] | 1034 | diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded) |
---|
[14072] | 1035 | |
---|
[14005] | 1036 | IF (zvolp(ji,jj) < epsi10) THEN |
---|
| 1037 | zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj) |
---|
| 1038 | zvolp(ji,jj) = 0._wp |
---|
| 1039 | END IF |
---|
| 1040 | END IF |
---|
| 1041 | END DO |
---|
[14072] | 1042 | |
---|
[14005] | 1043 | ! adjust melt pond dimensions |
---|
| 1044 | IF (permflag > 0) THEN |
---|
| 1045 | ! recompute pond depth |
---|
| 1046 | CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index) |
---|
| 1047 | DO jl = 1, m_index |
---|
| 1048 | h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1) |
---|
| 1049 | a_ip(ji,jj,jl) = reduced_aicen(jl) |
---|
| 1050 | END DO |
---|
| 1051 | !zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d |
---|
| 1052 | END IF |
---|
| 1053 | END IF ! pressure_head |
---|
[14072] | 1054 | |
---|
[14005] | 1055 | !------------------------------- |
---|
| 1056 | ! remove water from the snow |
---|
| 1057 | !------------------------------- |
---|
| 1058 | !------------------------------------------------------------------------ |
---|
| 1059 | ! total melt pond volume in category does not include snow volume |
---|
| 1060 | ! snow in melt ponds is not melted |
---|
| 1061 | !------------------------------------------------------------------------ |
---|
[14072] | 1062 | |
---|
[14005] | 1063 | ! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell |
---|
| 1064 | ! how much, so I did not diagnose it |
---|
| 1065 | ! so if there is a problem here, nobody is going to see it... |
---|
[14072] | 1066 | |
---|
| 1067 | |
---|
[14005] | 1068 | ! Calculate pond volume for lower categories |
---|
| 1069 | DO jl = 1,m_index-1 |
---|
| 1070 | v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow |
---|
| 1071 | - (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl)) |
---|
| 1072 | END DO |
---|
[14072] | 1073 | |
---|
[14005] | 1074 | ! Calculate pond volume for highest category = remaining pond volume |
---|
[14072] | 1075 | |
---|
[14005] | 1076 | ! The following is completely unclear to Martin at least |
---|
| 1077 | ! Could we redefine properly and recode in a more readable way ? |
---|
[14072] | 1078 | |
---|
[14005] | 1079 | ! m_index = last category with melt pond |
---|
[14072] | 1080 | |
---|
[14005] | 1081 | 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 |
---|
[14072] | 1082 | |
---|
[14005] | 1083 | IF (m_index > 1) THEN |
---|
| 1084 | IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN |
---|
| 1085 | v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1)) |
---|
| 1086 | ELSE |
---|
[14072] | 1087 | v_ip(ji,jj,m_index) = 0._wp |
---|
[14005] | 1088 | h_ip(ji,jj,m_index) = 0._wp |
---|
| 1089 | a_ip(ji,jj,m_index) = 0._wp |
---|
| 1090 | ! If remaining pond volume is negative reduce pond volume of |
---|
| 1091 | ! lower category |
---|
| 1092 | IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) & |
---|
| 1093 | 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) |
---|
| 1094 | END IF |
---|
| 1095 | END IF |
---|
[14072] | 1096 | |
---|
[14005] | 1097 | DO jl = 1,m_index |
---|
| 1098 | IF (a_ip(ji,jj,jl) > epsi10) THEN |
---|
| 1099 | h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl) |
---|
| 1100 | ELSE |
---|
| 1101 | zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl) |
---|
[14072] | 1102 | h_ip(ji,jj,jl) = 0._wp |
---|
[14005] | 1103 | v_ip(ji,jj,jl) = 0._wp |
---|
| 1104 | a_ip(ji,jj,jl) = 0._wp |
---|
| 1105 | END IF |
---|
| 1106 | END DO |
---|
| 1107 | DO jl = m_index+1, jpl |
---|
[14072] | 1108 | h_ip(ji,jj,jl) = 0._wp |
---|
| 1109 | a_ip(ji,jj,jl) = 0._wp |
---|
| 1110 | v_ip(ji,jj,jl) = 0._wp |
---|
[14005] | 1111 | END DO |
---|
[14072] | 1112 | |
---|
[14005] | 1113 | ENDIF |
---|
| 1114 | |
---|
| 1115 | END_2D |
---|
| 1116 | |
---|
| 1117 | END SUBROUTINE ice_thd_pnd_area |
---|
| 1118 | |
---|
| 1119 | |
---|
| 1120 | SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) |
---|
| 1121 | !!------------------------------------------------------------------- |
---|
| 1122 | !! *** ROUTINE ice_thd_pnd_depth *** |
---|
| 1123 | !! |
---|
| 1124 | !! ** Purpose : Compute melt pond depth |
---|
| 1125 | !!------------------------------------------------------------------- |
---|
| 1126 | |
---|
| 1127 | REAL (wp), DIMENSION(jpl), INTENT(IN) :: & |
---|
| 1128 | aicen, & |
---|
| 1129 | asnon, & |
---|
| 1130 | hsnon, & |
---|
| 1131 | alfan, & |
---|
| 1132 | cum_max_vol |
---|
| 1133 | |
---|
| 1134 | REAL (wp), INTENT(IN) :: & |
---|
| 1135 | zvolp |
---|
| 1136 | |
---|
| 1137 | REAL (wp), INTENT(OUT) :: & |
---|
| 1138 | hpond |
---|
| 1139 | |
---|
| 1140 | INTEGER, INTENT(OUT) :: & |
---|
| 1141 | m_index |
---|
| 1142 | |
---|
| 1143 | INTEGER :: n, ns |
---|
| 1144 | |
---|
| 1145 | REAL (wp), DIMENSION(0:jpl+1) :: & |
---|
| 1146 | hitl, & |
---|
| 1147 | aicetl |
---|
| 1148 | |
---|
| 1149 | REAL (wp) :: & |
---|
| 1150 | rem_vol, & |
---|
| 1151 | area, & |
---|
| 1152 | vol, & |
---|
| 1153 | tmp, & |
---|
| 1154 | z0 = 0.0_wp |
---|
| 1155 | |
---|
| 1156 | !---------------------------------------------------------------- |
---|
| 1157 | ! hpond is zero if zvolp is zero - have we fully drained? |
---|
| 1158 | !---------------------------------------------------------------- |
---|
| 1159 | |
---|
| 1160 | IF (zvolp < epsi10) THEN |
---|
| 1161 | hpond = z0 |
---|
| 1162 | m_index = 0 |
---|
| 1163 | ELSE |
---|
| 1164 | |
---|
| 1165 | !---------------------------------------------------------------- |
---|
| 1166 | ! Calculate the category where water fills up to |
---|
| 1167 | !---------------------------------------------------------------- |
---|
| 1168 | |
---|
| 1169 | !----------| |
---|
| 1170 | ! | |
---|
| 1171 | ! | |
---|
| 1172 | ! |----------| -- -- |
---|
| 1173 | !__________|__________|_________________________________________ ^ |
---|
| 1174 | ! | | rem_vol ^ | Semi-filled |
---|
| 1175 | ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer |
---|
| 1176 | ! | | | | |
---|
| 1177 | ! | | | |hpond |
---|
| 1178 | ! | | |----------| | |------- |
---|
| 1179 | ! | | | | | | |
---|
| 1180 | ! | | | |---v-----| |
---|
| 1181 | ! | | m_index | | | |
---|
| 1182 | !------------------------------------------------------------- |
---|
| 1183 | |
---|
| 1184 | m_index = 0 ! 1:m_index categories have water in them |
---|
| 1185 | DO n = 1, jpl |
---|
| 1186 | IF (zvolp <= cum_max_vol(n)) THEN |
---|
| 1187 | m_index = n |
---|
| 1188 | IF (n == 1) THEN |
---|
| 1189 | rem_vol = zvolp |
---|
| 1190 | ELSE |
---|
| 1191 | rem_vol = zvolp - cum_max_vol(n-1) |
---|
| 1192 | END IF |
---|
| 1193 | exit ! to break out of the loop |
---|
| 1194 | END IF |
---|
| 1195 | END DO |
---|
| 1196 | m_index = min(jpl-1, m_index) |
---|
| 1197 | |
---|
| 1198 | !---------------------------------------------------------------- |
---|
| 1199 | ! semi-filled layer may have m_index different snow in it |
---|
| 1200 | !---------------------------------------------------------------- |
---|
| 1201 | |
---|
| 1202 | !----------------------------------------------------------- ^ |
---|
| 1203 | ! | alfan(m_index+1) |
---|
| 1204 | ! | |
---|
| 1205 | !hitl(3)--> |----------| | |
---|
| 1206 | !hitl(2)--> |------------| * * * * *| | |
---|
| 1207 | !hitl(1)--> |----------|* * * * * * |* * * * * | | |
---|
| 1208 | !hitl(0)-->------------------------------------------------- | ^ |
---|
| 1209 | ! various snow from lower categories | |alfa(m_index) |
---|
| 1210 | |
---|
| 1211 | ! hitl - heights of the snow layers from thinner and current categories |
---|
| 1212 | ! aicetl - area of each snow depth in this layer |
---|
| 1213 | |
---|
| 1214 | hitl(:) = z0 |
---|
| 1215 | aicetl(:) = z0 |
---|
| 1216 | DO n = 1, m_index |
---|
| 1217 | hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & |
---|
| 1218 | alfan(m_index+1) - alfan(m_index)), z0) |
---|
| 1219 | aicetl(n) = asnon(n) |
---|
| 1220 | |
---|
| 1221 | aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) |
---|
| 1222 | END DO |
---|
| 1223 | |
---|
| 1224 | hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) |
---|
| 1225 | aicetl(m_index+1) = z0 |
---|
| 1226 | |
---|
| 1227 | !---------------------------------------------------------------- |
---|
| 1228 | ! reorder array according to hitl |
---|
| 1229 | ! snow heights not necessarily in height order |
---|
| 1230 | !---------------------------------------------------------------- |
---|
| 1231 | |
---|
| 1232 | DO ns = 1, m_index+1 |
---|
| 1233 | DO n = 0, m_index - ns + 1 |
---|
| 1234 | IF (hitl(n) > hitl(n+1)) THEN ! swap order |
---|
| 1235 | tmp = hitl(n) |
---|
| 1236 | hitl(n) = hitl(n+1) |
---|
| 1237 | hitl(n+1) = tmp |
---|
| 1238 | tmp = aicetl(n) |
---|
| 1239 | aicetl(n) = aicetl(n+1) |
---|
| 1240 | aicetl(n+1) = tmp |
---|
| 1241 | END IF |
---|
| 1242 | END DO |
---|
| 1243 | END DO |
---|
| 1244 | |
---|
| 1245 | !---------------------------------------------------------------- |
---|
| 1246 | ! divide semi-filled layer into set of sublayers each vertically homogenous |
---|
| 1247 | !---------------------------------------------------------------- |
---|
| 1248 | |
---|
| 1249 | !hitl(3)---------------------------------------------------------------- |
---|
| 1250 | ! | * * * * * * * * |
---|
| 1251 | ! |* * * * * * * * * |
---|
| 1252 | !hitl(2)---------------------------------------------------------------- |
---|
| 1253 | ! | * * * * * * * * | * * * * * * * * |
---|
| 1254 | ! |* * * * * * * * * |* * * * * * * * * |
---|
| 1255 | !hitl(1)---------------------------------------------------------------- |
---|
| 1256 | ! | * * * * * * * * | * * * * * * * * | * * * * * * * * |
---|
| 1257 | ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * |
---|
| 1258 | !hitl(0)---------------------------------------------------------------- |
---|
| 1259 | ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) |
---|
| 1260 | |
---|
| 1261 | ! move up over layers incrementing volume |
---|
| 1262 | DO n = 1, m_index+1 |
---|
| 1263 | |
---|
| 1264 | area = sum(aicetl(:)) - & ! total area of sub-layer |
---|
| 1265 | (rhos/rho0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow |
---|
| 1266 | |
---|
| 1267 | vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area |
---|
| 1268 | |
---|
| 1269 | IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within |
---|
| 1270 | hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) |
---|
| 1271 | |
---|
| 1272 | exit |
---|
| 1273 | ELSE ! still in sub-layer below the sub-layer with the depth |
---|
| 1274 | rem_vol = rem_vol - vol |
---|
| 1275 | END IF |
---|
| 1276 | |
---|
| 1277 | END DO |
---|
| 1278 | |
---|
| 1279 | END IF |
---|
| 1280 | |
---|
| 1281 | END SUBROUTINE ice_thd_pnd_depth |
---|
| 1282 | |
---|
| 1283 | |
---|
| 1284 | SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) |
---|
| 1285 | !!------------------------------------------------------------------- |
---|
| 1286 | !! *** ROUTINE ice_thd_pnd_perm *** |
---|
| 1287 | !! |
---|
| 1288 | !! ** Purpose : Determine the liquid fraction of brine in the ice |
---|
| 1289 | !! and its permeability |
---|
| 1290 | !!------------------------------------------------------------------- |
---|
| 1291 | |
---|
| 1292 | REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & |
---|
| 1293 | ticen, & ! internal ice temperature (K) |
---|
| 1294 | salin ! salinity (ppt) !js: ppt according to cice |
---|
| 1295 | |
---|
| 1296 | REAL (wp), INTENT(OUT) :: & |
---|
| 1297 | perm ! permeability |
---|
| 1298 | |
---|
| 1299 | REAL (wp) :: & |
---|
| 1300 | Sbr ! brine salinity |
---|
| 1301 | |
---|
| 1302 | REAL (wp), DIMENSION(nlay_i) :: & |
---|
| 1303 | Tin, & ! ice temperature |
---|
| 1304 | phi ! liquid fraction |
---|
| 1305 | |
---|
| 1306 | INTEGER :: k |
---|
| 1307 | |
---|
| 1308 | !----------------------------------------------------------------- |
---|
| 1309 | ! Compute ice temperatures from enthalpies using quadratic formula |
---|
| 1310 | !----------------------------------------------------------------- |
---|
| 1311 | |
---|
| 1312 | DO k = 1,nlay_i |
---|
| 1313 | Tin(k) = ticen(k) - rt0 !js: from K to degC |
---|
| 1314 | END DO |
---|
| 1315 | |
---|
| 1316 | !----------------------------------------------------------------- |
---|
| 1317 | ! brine salinity and liquid fraction |
---|
| 1318 | !----------------------------------------------------------------- |
---|
| 1319 | |
---|
| 1320 | DO k = 1, nlay_i |
---|
[14072] | 1321 | |
---|
[14005] | 1322 | Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) |
---|
| 1323 | ! Best expression to date is that one (Vancoppenolle et al JGR 2019) |
---|
| 1324 | ! Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)**2 - 0.00535 * Tin(k) **3 |
---|
| 1325 | phi(k) = salin(k) / Sbr |
---|
[14072] | 1326 | |
---|
[14005] | 1327 | END DO |
---|
| 1328 | |
---|
| 1329 | !----------------------------------------------------------------- |
---|
| 1330 | ! permeability |
---|
| 1331 | !----------------------------------------------------------------- |
---|
| 1332 | |
---|
| 1333 | perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) |
---|
| 1334 | |
---|
| 1335 | END SUBROUTINE ice_thd_pnd_perm |
---|
| 1336 | |
---|
[14072] | 1337 | SUBROUTINE ice_thd_pnd_init |
---|
[8637] | 1338 | !!------------------------------------------------------------------- |
---|
| 1339 | !! *** ROUTINE ice_thd_pnd_init *** |
---|
| 1340 | !! |
---|
| 1341 | !! ** Purpose : Physical constants and parameters linked to melt ponds |
---|
| 1342 | !! over sea ice |
---|
| 1343 | !! |
---|
[14072] | 1344 | !! ** Method : Read the namthd_pnd namelist and check the melt pond |
---|
[8637] | 1345 | !! parameter values called at the first timestep (nit000) |
---|
| 1346 | !! |
---|
[14072] | 1347 | !! ** input : Namelist namthd_pnd |
---|
[8637] | 1348 | !!------------------------------------------------------------------- |
---|
[9169] | 1349 | INTEGER :: ios, ioptio ! Local integer |
---|
| 1350 | !! |
---|
[14005] | 1351 | NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, & |
---|
[13472] | 1352 | & ln_pnd_CST , rn_apnd, rn_hpnd, & |
---|
[14005] | 1353 | & ln_pnd_TOPO, & |
---|
[13472] | 1354 | & ln_pnd_lids, ln_pnd_alb |
---|
[8637] | 1355 | !!------------------------------------------------------------------- |
---|
[9169] | 1356 | ! |
---|
[8637] | 1357 | READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901) |
---|
[11536] | 1358 | 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in reference namelist' ) |
---|
[8637] | 1359 | READ ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 ) |
---|
[11536] | 1360 | 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' ) |
---|
[8637] | 1361 | IF(lwm) WRITE ( numoni, namthd_pnd ) |
---|
[9169] | 1362 | ! |
---|
[8637] | 1363 | IF(lwp) THEN ! control print |
---|
| 1364 | WRITE(numout,*) |
---|
| 1365 | WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds' |
---|
| 1366 | WRITE(numout,*) '~~~~~~~~~~~~~~~~' |
---|
| 1367 | WRITE(numout,*) ' Namelist namicethd_pnd:' |
---|
[13472] | 1368 | WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd |
---|
[14005] | 1369 | WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO |
---|
[13472] | 1370 | WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV |
---|
| 1371 | WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min |
---|
| 1372 | WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max |
---|
[14005] | 1373 | WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush |
---|
[13472] | 1374 | WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST |
---|
| 1375 | WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd |
---|
| 1376 | WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd |
---|
| 1377 | WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids |
---|
| 1378 | WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb |
---|
[8637] | 1379 | ENDIF |
---|
| 1380 | ! |
---|
| 1381 | ! !== set the choice of ice pond scheme ==! |
---|
| 1382 | ioptio = 0 |
---|
[11536] | 1383 | IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF |
---|
| 1384 | IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF |
---|
[13472] | 1385 | IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF |
---|
[14005] | 1386 | IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF |
---|
[11536] | 1387 | IF( ioptio /= 1 ) & |
---|
[14005] | 1388 | & 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] | 1389 | ! |
---|
[8637] | 1390 | SELECT CASE( nice_pnd ) |
---|
[14072] | 1391 | CASE( np_pndNO ) |
---|
[13472] | 1392 | IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF |
---|
| 1393 | IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF |
---|
[14072] | 1394 | CASE( np_pndCST ) |
---|
[13472] | 1395 | IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF |
---|
[8637] | 1396 | END SELECT |
---|
| 1397 | ! |
---|
| 1398 | END SUBROUTINE ice_thd_pnd_init |
---|
[14072] | 1399 | |
---|
[8637] | 1400 | #else |
---|
| 1401 | !!---------------------------------------------------------------------- |
---|
[9570] | 1402 | !! Default option Empty module NO SI3 sea-ice model |
---|
[8637] | 1403 | !!---------------------------------------------------------------------- |
---|
[14072] | 1404 | #endif |
---|
[8637] | 1405 | |
---|
| 1406 | !!====================================================================== |
---|
[14072] | 1407 | END MODULE icethd_pnd |
---|