MODULE icethd_pnd !!====================================================================== !! *** MODULE icethd_pnd *** !! sea-ice: Melt ponds on top of sea ice !!====================================================================== !! history : ! 2012 (O. Lecomte) Adaptation from Flocco and Turner !! ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation !! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube] !!---------------------------------------------------------------------- #if defined key_si3 !!---------------------------------------------------------------------- !! 'key_si3' : SI3 sea-ice model !!---------------------------------------------------------------------- !! ice_thd_pnd_init : some initialization and namelist read !! ice_thd_pnd : main calling routine !!---------------------------------------------------------------------- USE phycst ! physical constants USE dom_oce ! ocean space and time domain USE ice ! sea-ice: variables USE ice1D ! sea-ice: thermodynamics variables USE icetab ! sea-ice: 1D <==> 2D transformation ! USE in_out_manager ! I/O manager USE lib_mpp ! MPP library USE lib_fortran ! fortran utilities (glob_sum + no signed zero) USE timing ! Timing IMPLICIT NONE PRIVATE PUBLIC ice_thd_pnd_init ! routine called by icestp.F90 PUBLIC ice_thd_pnd ! routine called by icestp.F90 INTEGER :: nice_pnd ! choice of the type of pond scheme ! ! associated indices: INTEGER, PARAMETER :: np_pndNO = 0 ! No pond scheme INTEGER, PARAMETER :: np_pndCST = 1 ! Constant ice pond scheme INTEGER, PARAMETER :: np_pndLEV = 2 ! Level ice pond scheme INTEGER, PARAMETER :: np_pndTOPO = 3 ! Level ice pond scheme !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2018) !! $Id$ !! Software governed by the CeCILL license (see ./LICENSE) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_thd_pnd !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_pnd *** !! !! ** Purpose : change melt pond fraction and thickness !! !! Note: Melt ponds affect currently affect radiative transfer !! They carry no heat, and the melt water they carry is not !! exchanged with the ocean !! !! This means freshwater is directly released after surface & bottom melt in ice_thd_dh !! !! A wfx_pnd has been coded for diagnostic purposes !! Each time wfx_pnd is updated, wfx_sum / wfx_snw_sum must be updated !! !!! The current diagnostic lacks a contribution from drainage !! !!------------------------------------------------------------------- ! SELECT CASE ( nice_pnd ) ! CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! ! CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==! ! CASE (np_pndTOPO) ; CALL pnd_TOPO & !== Topographic melt ponds ==! ( at_i, a_i, & & vt_i, v_i, & & v_s, & & t_i, s_i, & & a_ip_frac, h_ip, & & t_su ) ! END SELECT ! END SUBROUTINE ice_thd_pnd SUBROUTINE pnd_CST !!------------------------------------------------------------------- !! *** ROUTINE pnd_CST *** !! !! ** Purpose : Compute melt pond evolution !! !! ** Method : Melt pond fraction and thickness are prescribed !! to non-zero values when t_su = 0C !! !! ** Tunable parameters : pond fraction (rn_apnd), pond depth (rn_hpnd) !! !! ** Note : Coupling with such melt ponds is only radiative !! Advection, ridging, rafting... are bypassed !! !! ** References : Bush, G.W., and Trump, D.J. (2017) !!------------------------------------------------------------------- INTEGER :: ji ! loop indices !!------------------------------------------------------------------- DO ji = 1, npti ! IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN h_ip_1d(ji) = rn_hpnd a_ip_1d(ji) = rn_apnd * a_i_1d(ji) h_il_1d(ji) = 0._wp ! no pond lids whatsoever ELSE h_ip_1d(ji) = 0._wp a_ip_1d(ji) = 0._wp h_il_1d(ji) = 0._wp ENDIF ! END DO ! END SUBROUTINE pnd_CST SUBROUTINE pnd_LEV !!------------------------------------------------------------------- !! *** ROUTINE pnd_LEV *** !! !! ** Purpose : Compute melt pond evolution !! !! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing !! We work with volumes and then redistribute changes into thickness and concentration !! assuming linear relationship between the two. !! !! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 --- !! dVmelt = (1-r)/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i !! dh_i = meltwater from ice surface melt !! dh_s = meltwater from snow melt !! (1-r) = fraction of melt water that is not flushed !! !! - limtations: a_ip must not exceed (1-r)*a_i !! h_ip must not exceed 0.5*h_i !! !! - pond shrinking: !! if lids: Vp = Vp -dH * a_ip !! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 --- !! !! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H !! H = lid thickness !! Lf = latent heat of fusion !! Tp = -2C !! !! And solved implicitely as: !! H(t+dt)**2 -H(t) * H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0 !! !! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 --- !! !! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi --- from Flocco et al 2007 --- !! perm = permability of sea-ice !! visc = water viscosity !! Hp = height of top of the pond above sea-level !! Hi = ice thickness thru which there is flushing !! !! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness !! !! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip: !! a_ip/a_i = a_ip_frac = h_ip / zaspect !! !! ** Tunable parameters : ln_pnd_lids, rn_apnd_max, rn_apnd_min !! !! ** Note : mostly stolen from CICE !! !! ** References : Flocco and Feltham (JGR, 2007) !! Flocco et al (JGR, 2010) !! Holland et al (J. Clim, 2012) !!------------------------------------------------------------------- REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array !! REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity !! REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh REAL(wp) :: zv_ip_max ! max pond volume allowed REAL(wp) :: zdT ! zTp-t_su REAL(wp) :: zsbr ! Brine salinity REAL(wp) :: zperm ! permeability of sea ice REAL(wp) :: zfac, zdum ! temporary arrays REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse !! INTEGER :: ji, jk ! loop indices !!------------------------------------------------------------------- z1_rhow = 1._wp / rhow z1_aspect = 1._wp / zaspect z1_Tp = 1._wp / zTp DO ji = 1, npti ! !----------------------------------------------------! IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction ! ! !----------------------------------------------------! !--- Remove ponds on thin ice or tiny ice fractions a_ip_1d(ji) = 0._wp h_ip_1d(ji) = 0._wp h_il_1d(ji) = 0._wp ! !--------------------------------! ELSE ! Case ice thickness >= rn_himin ! ! !--------------------------------! v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji) ! !------------------! ! case ice melting ! !------------------! ! !--- available meltwater for melt ponding ---! zdum = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow * a_i_1d(ji) 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 zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! >0, max for roundoff errors? ! !--- overflow ---! ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume ! a_ip_max = zfr_mlt * a_i ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: zv_ip_max = zfr_mlt**2 * a_i_1d(ji) * zaspect zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) ! If pond depth exceeds half the ice thickness then reduce the pond volume ! h_ip_max = 0.5 * h_i ! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as: 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 zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) ) !--- Pond growing ---! v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt ! !--- Lid melting ---! IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0 ! !--- mass flux ---! IF( zdv_mlt > 0._wp ) THEN ! MV add comment on what that mass flux means zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1] wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac ! ! MV ! why surface melt and snow fluxes must be adjusted is not clear ! sounds like things are counted twice ! zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum) wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum) ENDIF !-------------------! ! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0) !-------------------! ! zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) ! !--- Pond contraction (due to refreezing) ---! IF( ln_pnd_lids ) THEN ! !--- Lid growing and subsequent pond shrinking ---! zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0 & SQRT( v_il_1d(ji)**2 + a_ip_1d(ji)**2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors ! Lid growing v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz ) ! Pond shrinking v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz ) ELSE ! Pond shrinking v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6) ENDIF ! !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac ! v_ip = h_ip * a_ip ! 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) 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 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) !---------------! ! Pond flushing ! !---------------! ! height of top of the pond above sea-level zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0 ! Calculate the permeability of the ice (Assur 1958, see Flocco 2010) DO jk = 1, nlay_i zsbr = - 1.2_wp & & - 21.8_wp * ( t_i_1d(ji,jk) - rt0 ) & & - 0.919_wp * ( t_i_1d(ji,jk) - rt0 )**2 & & - 0.0178_wp * ( t_i_1d(ji,jk) - rt0 )**3 ! zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt ztmp(jk) = sz_i_1d(ji,jk) / zsbr END DO zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 ) ! Do the drainage using Darcy's law zdv_flush = -zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) ) v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush ! MV --- why pond drainage does not give back water into freshwater flux ? !--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac 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 h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji) !--- Corrections and lid thickness ---! IF( ln_pnd_lids ) THEN !--- retrieve lid thickness from volume ---! IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji) ELSE ; h_il_1d(ji) = 0._wp ENDIF !--- remove ponds if lids are much larger than ponds ---! IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN a_ip_1d(ji) = 0._wp h_ip_1d(ji) = 0._wp h_il_1d(ji) = 0._wp ENDIF ENDIF ! ENDIF END DO ! END SUBROUTINE pnd_LEV SUBROUTINE pnd_TOPO (aice, aicen, & vice, vicen, & vsnon, & ticen, salin, & a_ip_frac, h_ip, & Tsfc ) !!------------------------------------------------------------------- !! *** ROUTINE pnd_TOPO *** !! !! ** Purpose : Compute melt pond evolution !! !! ** Purpose : Compute melt pond evolution based on the ice !! topography as inferred from the ice thickness !! distribution. !! !! ** Method : This code is initially based on Flocco and Feltham !! (2007) and Flocco et al. (2010). More to come... !! !! ** Tunable parameters : !! !! ** Note : !! !! ** References !! Flocco, D. and D. L. Feltham, 2007. A continuum model of melt pond !! evolution on Arctic sea ice. J. Geophys. Res. 112, C08016, doi: !! 10.1029/2006JC003836. !! Flocco, D., D. L. Feltham and A. K. Turner, 2010. Incorporation of !! a physically based melt pond scheme into the sea ice component of a !! climate model. J. Geophys. Res. 115, C08012, !! doi: 10.1029/2009JC005568. !! !!------------------------------------------------------------------- !js 190423: the lid on melt ponds appears only in the analog subroutine of CICE 5.1.2 REAL (wp), DIMENSION (jpi,jpj), & INTENT(IN) :: & aice, & ! total ice area fraction vice ! total ice volume (m) REAL (wp), DIMENSION (jpi,jpj,jpl), & INTENT(IN) :: & aicen, & ! ice area fraction, per category vsnon, & ! snow volume, per category (m) vicen ! ice volume, per category (m) REAL (wp), DIMENSION (jpi,jpj,nlay_i,jpl), & INTENT(IN) :: & ticen, & ! ice temperature per category (K) salin REAL (wp), DIMENSION (jpi,jpj,jpl), & INTENT(INOUT) :: & a_ip_frac , & ! pond area fraction of ice, per ice category h_ip ! pond depth, per ice category (m) REAL (wp), DIMENSION (jpi,jpj,jpl), & INTENT(IN) :: & Tsfc ! snow/sea ice surface temperature ! local variables REAL (wp), DIMENSION (jpi,jpj,jpl) :: & zTsfcn, & ! ice/snow surface temperature (C) zvolpn ! pond volume per unit area, per category (m) REAL (wp), DIMENSION (jpi,jpj,jpl) :: & zrfrac, & ! fraction of available meltwater retained for melt ponding zapondn,& ! pond area fraction, per category zhpondn ! pond depth, per category (m) REAL (wp), DIMENSION (jpi,jpj) :: & zvolp, & ! total volume of pond, per unit area of pond (m) zwfx_tmp ! temporary array for melt water REAL (wp) :: & zhi, & ! ice thickness (m) zTavg, & ! mean surface temperature across categories (C) z1_rhow, & ! inverse freshwater density zdTs, & ! temperature difference for freeze-up (C) zvpold, & ! dummy pond volume zdvn ! change in melt pond volume for fresh water budget INTEGER, DIMENSION (jpi*jpj) :: & indxi, indxj ! compressed indices for cells with ice melting INTEGER :: ij,icells,ji,jj,jl ! loop indices REAL (wp), PARAMETER :: & ! MV ouate de phoque!!! constants hard coded ???? ! 917 = rhoi ! 0.334 =Lfus!!!! zr1_rlfus = 1._wp / 0.334e+6 / 917._wp , & ! (J/m^3) zTp = -0.15_wp, & ! pond freezing temperature (C) zmin_volp = 1.e-4_wp, & ! minimum pond volume (m) zrexp = 0.01_wp, & ! constant melt pond freeze-up rate z01 = 0.01_wp, & z25 = 0.25_wp, & z5 = 0.5_wp z1_rhow = 1. / rhow !--------------------------------------------------------------- ! Initialization !--------------------------------------------------------------- zhpondn (:,:,:) = 0._wp zapondn (:,:,:) = 0._wp zvolpn (:,:,:) = 0._wp zTsfcn(:,:,:) = Tsfc(:,:,:) - rt0 ! Convert in Celsius ! IF ( ln_pnd_fw ) THEN ! v_ip_b(:,:,:) = v_ip(:,:,:) ! ELSE ! v_ip_b(:,:,:) = 0._wp ! ENDIF !------------------------------------------------------------------ ! Available melt water for melt ponding and corresponding fraction !------------------------------------------------------------------ !js 03/05/19 unset restriction on sign of wfx_pnd_in; mask values close to zero for future division !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:), 0._wp ) ! available meltwater for melt ponding !wfx_pnd_in(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) * MAX( 0._wp, SIGN( 1._wp, wfx_sum(:,:) + wfx_snw_sum(:,:) - epsi10 ) ) !wfx_pnd_in(:,:) = (wfx_sum(:,:) + wfx_snw_sum(:,:)) * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) ) ! MV ! NB: wfx_pnd_in can be slightly negative for very small values (why?) ! This can in some occasions give negative ! v_ip in the first category, which then gives crazy pond ! fractions and crashes the code as soon as the melt-pond ! radiative coupling is activated ! if we understand and remove why wfx_sum or wfx_snw could be ! negative, then, we can remove the MAX ! NB: I now changed to wfx_snw_sum, this may fix the problem. ! We should check ! OLI 07/2017: when we (MV & OL) first started the inclusion of melt ! ponds in the model, we removed the Holland et al. (2012, see ! CESM scheme above) parameterization. I put it back here, ! because I think it is needed. In summary, the sinks of FW for ! ponds are: 1/ Runoff through cracks/leads => depends on the ! total ice area only ! 2/ overflow, including Lüthje et al. (2006) limitation ! (max a_ip fraction function of h_i). This is in ! fact an other form of runoff that depends on the ! ITD ! 3/ Flushing - losses by permeability ! 4/ Refreezing ! 5/ Removal of ponds on thin ice ! I think 1 is needed because it is different from 2. However, ! test runs should/could be done, to check the sensitivity and ! the real usefulness of that stuff. ! Note : the Holland et al. param was wrongly wired in NEMO3.1 (using ! a_i instead of at_i), which might well explain why I had a too ! weak melt pond cover in my simulations (compared to MODIS, in ! situ obs. and CICE simulations. !js 23/04/19: rewired back to a fraction with a_i !!! zrfrac(:,:,:) = rn_pnd_fracmin + ( rn_pnd_fracmax - rn_pnd_fracmin ) * aicen(:,:,:) ! MV2020 !!!jzrfrac(:,:,:) = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i_1d(:) ! = ( 1 - r ) = fraction of melt water that is not flushed !!!zwfx_tmp(:,:) = 0._wp !--- Add retained melt water to melt ponds ! v_ip should never be negative, otherwise code crashes ! Rq MV: as far as I saw, UM5 can create very small negative v_ip values ! hence I added the max, which was not required with Prather (1 yr run) ! OLI: Here I use vt_ip, so I don't know if the max is ! required... !zvolp(:,:) = MAX(vt_ip(:,:),0._wp) + zrfrac(:,:) * wfx_pnd_in(:,:) * z1_rhow * rdt_ice ! Total available melt water, to be distributed as melt ponds ! OLI: 07/2017 Bugfix above, removed " * aice(:,:)" !js 19/04/18: change zrfrac to use aicen zvolp(:,:) = vt_ip(:,:) DO jl = 1, jpl ! Melt water, to be distributed as melt ponds ! zvolp(:,:) = zvolp(:,:) - zrfrac(:,:,jl) & ! * ( dh_i_pnd(:,:,jl)*rhoi + dh_s_pnd(:,:,jl)*rhosn ) & ! * z1_rhow * a_i(:,:,jl) ! MV ---> use expression from level ice melt ponds (dv_mlt) END DO !!! MV ---> rewrite this accounting for level ice melt ponds !!! 1) wfx_pnd_in is obsolete, ln_pnd_fw as well !!! 2) !!! wfx_sum !!! !!! !js 03/05/19: we truncate negative values after calculating zvolp, in a ! similar manner to the subroutine ice_thd_pnd_cesm. Variation dh_i_pnd and ! dh_s_pnd are negative, indicating a loss of ice or snow. But we can expect them ! to be negative for some reasons. We keep this behaviour as it is, for ! fluxes conservation reasons. If some dh are positive, then we remove water ! indirectly from the ponds. zvolp(:,:) = MAX( zvolp(:,:) , 0._wp ) ! ! Fresh water flux going into the ponds ! wfx_pnd_in(:,:) = wfx_pnd_in(:,:) + rhow * ( zvolp(:,:) - vt_ip(:,:) ) * r1_rdtice !--- Remove retained meltwater from surface fluxes ! IF ( ln_pnd_fw ) THEN !wfx_snw_sum(:,:) = wfx_snw_sum(:,:) * ( 1. - zrfrac(:,:) ) !wfx_sum(:,:) = wfx_sum(:,:) * ( 1. - zrfrac(:,:) ) !js 190419: we change the code to use a_i in zrfrac. To be tested, but !it should be conservative. zwfx_tmp is the flux accumulated in the !ponds. wfx_pnd_in is the total surface melt fluxes. zwfx_tmp(:,:) = MAX( wfx_sum(:,:) + wfx_snw_sum(:,:) , epsi10 ) & & * MAX( 0._wp, SIGN( 1._wp, ABS(wfx_sum(:,:) + wfx_snw_sum(:,:)) - epsi10 ) ) ! WHERE ( ABS(zwfx_tmp(:,:)) > epsi10 ) ! zwfx_tmp(:,:) = wfx_pnd_in(:,:) / zwfx_tmp(:,:) ! ELSEWHERE ! zwfx_tmp(:,:) = 0._wp ! ENDWHERE wfx_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_sum(:,:) wfx_snw_sum(:,:) = ( 1._wp - zwfx_tmp(:,:) ) * wfx_snw_sum(:,:) ! ENDIF !----------------------------------------------------------------- ! Identify grid cells with ponds !----------------------------------------------------------------- icells = 0 DO jj = 1, jpj DO ji = 1, jpi IF ( aice(ji,jj) > epsi10 ) THEN zhi = vice(ji,jj) / aice(ji,jj) ELSE zhi = 0._wp END IF IF ( aice(ji,jj) > z01 .and. zhi > rn_himin .and. & zvolp(ji,jj) > zmin_volp*aice(ji,jj)) THEN icells = icells + 1 indxi(icells) = ji indxj(icells) = jj ELSE ! remove ponds on thin ice, or too small ponds zvolpn(ji,jj,:) = 0._wp zvolp (ji,jj) = 0._wp a_ip(ji,jj,:) = 0._wp v_ip(ji,jj,:) = 0._wp a_ip_frac(ji,jj,:) = 0._wp h_ip(ji,jj,:) = 0._wp vt_ip(ji,jj) = 0._wp at_ip(ji,jj) = 0._wp ! IF ( ln_pnd_fw ) & !--- Give freshwater to the ocean ! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zvolp(ji,jj) * rhow * r1_rdtice END IF END DO ! ji END DO ! jj !!! MV sounds like this should be replaced by lid DO ij = 1, icells ji = indxi(ij) jj = indxj(ij) !-------------------------------------------------------------- ! Shrink pond due to refreezing !-------------------------------------------------------------- ! OLI 07/2017: Done like for empirical melt pond scheme (CESM). ! Therefore, I chose to put this part of the code before the main ! routines ice_thd_pnd_area/depth (contrary to the original code), ! seeing the freeze-up as a global sink of ! freshwater for melt ponds in the whole grid cell. If this was done ! after, I would need to make an additional assumption on the shape of ! melt ponds, which I don't want to do (for the CESM scheme, this ! assumption was on the aspect ratio). So I remove some water due to ! refreezing first (using zTavg instead of zTsfcn in each category) and ! then let the FF07 routines do their job for the fractional areas and ! depths of melt ponds. ! The whole ice lid related stuff from FF07 was thus removed and replaced ! by this. As mentionned below, this should be improved, but is much ! easier to conserve heat and freshwater this way. ! Average surface temperature is needed to compute freeze-up at the cell ! scale zTavg = 0._wp DO jl = 1, jpl zTavg = zTavg + zTsfcn(ji,jj,jl)*aicen(ji,jj,jl) END DO zTavg = zTavg / aice(ji,jj) ! The freezing temperature for meltponds is assumed slightly below 0C, ! as if meltponds had a little salt in them (hence the use of zTp). ! The salt budget is not ! altered for meltponds, but if it were then an actual pond freezing ! temperature could be computed. zdTs = MAX ( zTp - zTavg, 0. ) zvpold = zvolp(ji,jj) zvolp(ji,jj) = zvolp(ji,jj) * EXP( zrexp * zdTs / zTp ) !--- Dump meltwater due to refreezing ( of course this is wrong !--- but this parameterization is too simple ) ! IF ( ln_pnd_fw ) & ! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + rhow * ( zvpold - zvolp(ji,jj) ) * r1_rdtice ! ! OLI 07/2017 : Bugfix above, zvpold - zvolp instead of the ! ! opposite, otherwise negative contribution !-------------------------------------------------------------- ! calculate pond area and depth !-------------------------------------------------------------- zdvn = 0._wp CALL ice_thd_pnd_area(aice(ji,jj),vice(ji,jj), & aicen(ji,jj,:), vicen(ji,jj,:), vsnon(ji,jj,:), & ticen(ji,jj,:,:), salin(ji,jj,:,:), & zvolpn(ji,jj,:), zvolp(ji,jj), & zapondn(ji,jj,:),zhpondn(ji,jj,:), zdvn) ! outputs are ! - zdvn ! - zvolpn ! - zvolp ! - zapondn ! - zhpondn ! IF ( ln_pnd_fw ) & ! wfx_pnd_out(ji,jj) = wfx_pnd_out(ji,jj) + zdvn * rhow * r1_rdtice ! update flux from ponds to ocean !--------------------------------------------------------------- ! Update pond volume and fraction !--------------------------------------------------------------- DO jl = 1, jpl a_ip(ji,jj,jl) = zapondn(ji,jj,jl) v_ip(ji,jj,jl) = zvolpn(ji,jj,jl) a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / MAX(aicen(ji,jj,jl), epsi10) & * MAX(0._wp, SIGN(1._wp, aicen(ji,jj,jl) - epsi10)) h_ip (ji,jj,jl) = zhpondn(ji,jj,jl) END DO END DO ! ij ! IF ( ln_pnd_fw ) THEN !js 15/05/19: water going out of the ponds give a positive freshwater ! flux. ! wfx_pnd_out(:,:) = SUM(MAX(0._wp, v_ip_b(:,:,:) - v_ip(:,:,:)), DIM=3) * rhow * r1_rdtice ! ELSE ! wfx_pnd_out(:,:) = 0._wp ! ENDIF END SUBROUTINE pnd_TOPO SUBROUTINE ice_thd_pnd_area(aice, vice, & aicen, vicen, vsnon, ticen, & salin, zvolpn, zvolp, & zapondn,zhpondn,dvolp) !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_pnd_area *** !! !! ** Purpose : Given the total volume of meltwater, update !! pond fraction (a_ip) and depth (should be volume) !! !! ** !! !!------------------------------------------------------------------ REAL (wp), INTENT(IN) :: & aice, vice REAL (wp), DIMENSION(jpl), INTENT(IN) :: & aicen, vicen, vsnon REAL (wp), DIMENSION(nlay_i,jpl), INTENT(IN) :: & ticen, salin REAL (wp), DIMENSION(jpl), INTENT(INOUT) :: & zvolpn REAL (wp), INTENT(INOUT) :: & zvolp, dvolp REAL (wp), DIMENSION(jpl), INTENT(OUT) :: & zapondn, zhpondn INTEGER :: & n, ns, & m_index, & permflag REAL (wp), DIMENSION(jpl) :: & hicen, & hsnon, & asnon, & alfan, & betan, & cum_max_vol, & reduced_aicen REAL (wp), DIMENSION(0:jpl) :: & cum_max_vol_tmp REAL (wp) :: & hpond, & drain, & floe_weight, & pressure_head, & hsl_rel, & deltah, & perm, & msno REAL (wp), parameter :: & viscosity = 1.79e-3_wp, & ! kinematic water viscosity in kg/m/s z0 = 0.0_wp , & c1 = 1.0_wp , & p4 = 0.4_wp , & p6 = 0.6_wp !-----------| ! | ! |-----------| !___________|___________|______________________________________sea-level ! | | ! | |---^--------| ! | | | | ! | | | |-----------| |------- ! | | |alfan(n)| | | ! | | | | |--------------| ! | | | | | | !---------------------------v------------------------------------------- ! | | ^ | | | ! | | | | |--------------| ! | | |betan(n)| | | ! | | | |-----------| |------- ! | | | | ! | |---v------- | ! | | ! |-----------| ! | !-----------| !------------------------------------------------------------------- ! initialize !------------------------------------------------------------------- DO n = 1, jpl zapondn(n) = z0 zhpondn(n) = z0 !---------------------------------------- ! X) compute the effective snow fraction !---------------------------------------- IF (aicen(n) < epsi10) THEN hicen(n) = z0 hsnon(n) = z0 reduced_aicen(n) = z0 asnon(n) = z0 !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables ELSE hicen(n) = vicen(n) / aicen(n) hsnon(n) = vsnon(n) / aicen(n) reduced_aicen(n) = c1 ! n=jpl !js: initial code in NEMO_DEV !IF (n < jpl) reduced_aicen(n) = aicen(n) & ! * (-0.024_wp*hicen(n) + 0.832_wp) !js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large IF (n < jpl) reduced_aicen(n) = aicen(n) & * max(0.2_wp,(-0.024_wp*hicen(n) + 0.832_wp)) asnon(n) = reduced_aicen(n) ! effective snow fraction (empirical) ! MV should check whether this makes sense to have the same effective snow fraction in here ! OLI: it probably doesn't END IF ! This choice for alfa and beta ignores hydrostatic equilibium of categories. ! Hydrostatic equilibium of the entire ITD is accounted for below, assuming ! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all ! categories. alfa and beta partition the ITD - they are areas not thicknesses! ! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area. ! Here, alfa = 60% of the ice area (and since hice is constant in a category, ! alfan = 60% of the ice volume) in each category lies above the reference line, ! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required. ! MV: ! Note that this choice is not in the original FF07 paper and has been adopted in CICE ! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe ! Where does that choice come from ? => OLI : Coz' Chuck Norris said so... alfan(n) = 0.6 * hicen(n) betan(n) = 0.4 * hicen(n) cum_max_vol(n) = z0 cum_max_vol_tmp(n) = z0 END DO ! jpl cum_max_vol_tmp(0) = z0 drain = z0 dvolp = z0 !---------------------------------------------------------- ! x) Drain overflow water, update pond fraction and volume !---------------------------------------------------------- !-------------------------------------------------------------------------- ! the maximum amount of water that can be contained up to each ice category !-------------------------------------------------------------------------- ! MV ! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW) ! Then the excess volume cum_max_vol(jl) drains out of the system ! It should be added to wfx_pnd_out ! END MV !js 18/04/19: XXX do something about this flux thing DO n = 1, jpl-1 ! last category can not hold any volume IF (alfan(n+1) >= alfan(n) .and. alfan(n+1) > z0) THEN ! total volume in level including snow cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) + & (alfan(n+1) - alfan(n)) * sum(reduced_aicen(1:n)) ! subtract snow solid volumes from lower categories in current level DO ns = 1, n cum_max_vol_tmp(n) = cum_max_vol_tmp(n) & - rhos/rhow * & ! free air fraction that can be filled by water asnon(ns) * & ! effective areal fraction of snow in that category max(min(hsnon(ns)+alfan(ns)-alfan(n), alfan(n+1)-alfan(n)), z0) END DO ELSE ! assume higher categories unoccupied cum_max_vol_tmp(n) = cum_max_vol_tmp(n-1) END IF !IF (cum_max_vol_tmp(n) < z0) THEN ! CALL abort_ice('negative melt pond volume') !END IF END DO cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl) !---------------------------------------------------------------- ! is there more meltwater than can be held in the floe? !---------------------------------------------------------------- IF (zvolp >= cum_max_vol(jpl)) THEN drain = zvolp - cum_max_vol(jpl) + epsi10 zvolp = zvolp - drain ! update meltwater volume available dvolp = drain ! this is the drained water IF (zvolp < epsi10) THEN dvolp = dvolp + zvolp zvolp = z0 END IF END IF ! height and area corresponding to the remaining volume CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) DO n=1, m_index !zhpondn(n) = hpond - alfan(n) + alfan(1) ! here oui choulde update ! ! volume instead, no ? zhpondn(n) = max((hpond - alfan(n) + alfan(1)), z0) !js: from CICE 5.1.2 zapondn(n) = reduced_aicen(n) ! in practise, pond fraction depends on the empirical snow fraction ! so in turn on ice thickness END DO !zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d !------------------------------------------------------------------------ ! Drainage through brine network (permeability) !------------------------------------------------------------------------ !!! drainage due to ice permeability - Darcy's law ! sea water level msno = z0 DO n=1,jpl msno = msno + vsnon(n) * rhos END DO floe_weight = (msno + rhoi*vice + rau0*zvolp) / aice hsl_rel = floe_weight / rau0 & - ((sum(betan(:)*aicen(:))/aice) + alfan(1)) deltah = hpond - hsl_rel pressure_head = grav * rau0 * max(deltah, z0) ! drain if ice is permeable permflag = 0 IF (pressure_head > z0) THEN DO n = 1, jpl-1 IF (hicen(n) /= z0) THEN !IF (hicen(n) > z0) THEN !js: from CICE 5.1.2 perm = 0._wp ! MV ugly dummy patch CALL ice_thd_pnd_perm(ticen(:,n), salin(:,n), perm) IF (perm > z0) permflag = 1 drain = perm*zapondn(n)*pressure_head*rdt_ice / & (viscosity*hicen(n)) dvolp = dvolp + min(drain, zvolp) zvolp = max(zvolp - drain, z0) IF (zvolp < epsi10) THEN dvolp = dvolp + zvolp zvolp = z0 END IF END IF END DO ! adjust melt pond dimensions IF (permflag > 0) THEN ! recompute pond depth CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) DO n=1, m_index zhpondn(n) = hpond - alfan(n) + alfan(1) zapondn(n) = reduced_aicen(n) END DO !zapond = sum(zapondn(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d END IF END IF ! pressure_head !------------------------------- ! X) remove water from the snow !------------------------------- !------------------------------------------------------------------------ ! total melt pond volume in category does not include snow volume ! snow in melt ponds is not melted !------------------------------------------------------------------------ ! Calculate pond volume for lower categories DO n=1,m_index-1 zvolpn(n) = zapondn(n) * zhpondn(n) & ! what is not in the snow - (rhos/rhow) * asnon(n) * min(hsnon(n), zhpondn(n)) END DO ! Calculate pond volume for highest category = remaining pond volume ! The following is completely unclear to Martin at least ! Could we redefine properly and recode in a more readable way ? ! m_index = last category with melt pond IF (m_index == 1) zvolpn(m_index) = zvolp ! volume of mw in 1st category is the total volume of melt water IF (m_index > 1) THEN IF (zvolp > sum(zvolpn(1:m_index-1))) THEN zvolpn(m_index) = zvolp - sum(zvolpn(1:m_index-1)) ELSE zvolpn(m_index) = z0 zhpondn(m_index) = z0 zapondn(m_index) = z0 ! If remaining pond volume is negative reduce pond volume of ! lower category IF (zvolp+epsi10 < sum(zvolpn(1:m_index-1))) & zvolpn(m_index-1) = zvolpn(m_index-1) - sum(zvolpn(1:m_index-1)) + zvolp END IF END IF DO n=1,m_index IF (zapondn(n) > epsi10) THEN zhpondn(n) = zvolpn(n) / zapondn(n) ELSE dvolp = dvolp + zvolpn(n) zhpondn(n) = z0 zvolpn(n) = z0 zapondn(n) = z0 END IF END DO DO n = m_index+1, jpl zhpondn(n) = z0 zapondn(n) = z0 zvolpn (n) = z0 END DO END SUBROUTINE ice_thd_pnd_area SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index) !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_pnd_depth *** !! !! ** Purpose : Compute melt pond depth !!------------------------------------------------------------------- REAL (wp), DIMENSION(jpl), INTENT(IN) :: & aicen, & asnon, & hsnon, & alfan, & cum_max_vol REAL (wp), INTENT(IN) :: & zvolp REAL (wp), INTENT(OUT) :: & hpond INTEGER, INTENT(OUT) :: & m_index INTEGER :: n, ns REAL (wp), DIMENSION(0:jpl+1) :: & hitl, & aicetl REAL (wp) :: & rem_vol, & area, & vol, & tmp, & z0 = 0.0_wp !---------------------------------------------------------------- ! hpond is zero if zvolp is zero - have we fully drained? !---------------------------------------------------------------- IF (zvolp < epsi10) THEN hpond = z0 m_index = 0 ELSE !---------------------------------------------------------------- ! Calculate the category where water fills up to !---------------------------------------------------------------- !----------| ! | ! | ! |----------| -- -- !__________|__________|_________________________________________ ^ ! | | rem_vol ^ | Semi-filled ! | |----------|-- -- -- - ---|-- ---- -- -- --v layer ! | | | | ! | | | |hpond ! | | |----------| | |------- ! | | | | | | ! | | | |---v-----| ! | | m_index | | | !------------------------------------------------------------- m_index = 0 ! 1:m_index categories have water in them DO n = 1, jpl IF (zvolp <= cum_max_vol(n)) THEN m_index = n IF (n == 1) THEN rem_vol = zvolp ELSE rem_vol = zvolp - cum_max_vol(n-1) END IF exit ! to break out of the loop END IF END DO m_index = min(jpl-1, m_index) !---------------------------------------------------------------- ! semi-filled layer may have m_index different snow in it !---------------------------------------------------------------- !----------------------------------------------------------- ^ ! | alfan(m_index+1) ! | !hitl(3)--> |----------| | !hitl(2)--> |------------| * * * * *| | !hitl(1)--> |----------|* * * * * * |* * * * * | | !hitl(0)-->------------------------------------------------- | ^ ! various snow from lower categories | |alfa(m_index) ! hitl - heights of the snow layers from thinner and current categories ! aicetl - area of each snow depth in this layer hitl(:) = z0 aicetl(:) = z0 DO n = 1, m_index hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), & alfan(m_index+1) - alfan(m_index)), z0) aicetl(n) = asnon(n) aicetl(0) = aicetl(0) + (aicen(n) - asnon(n)) END DO hitl(m_index+1) = alfan(m_index+1) - alfan(m_index) aicetl(m_index+1) = z0 !---------------------------------------------------------------- ! reorder array according to hitl ! snow heights not necessarily in height order !---------------------------------------------------------------- DO ns = 1, m_index+1 DO n = 0, m_index - ns + 1 IF (hitl(n) > hitl(n+1)) THEN ! swap order tmp = hitl(n) hitl(n) = hitl(n+1) hitl(n+1) = tmp tmp = aicetl(n) aicetl(n) = aicetl(n+1) aicetl(n+1) = tmp END IF END DO END DO !---------------------------------------------------------------- ! divide semi-filled layer into set of sublayers each vertically homogenous !---------------------------------------------------------------- !hitl(3)---------------------------------------------------------------- ! | * * * * * * * * ! |* * * * * * * * * !hitl(2)---------------------------------------------------------------- ! | * * * * * * * * | * * * * * * * * ! |* * * * * * * * * |* * * * * * * * * !hitl(1)---------------------------------------------------------------- ! | * * * * * * * * | * * * * * * * * | * * * * * * * * ! |* * * * * * * * * |* * * * * * * * * |* * * * * * * * * !hitl(0)---------------------------------------------------------------- ! aicetl(0) aicetl(1) aicetl(2) aicetl(3) ! move up over layers incrementing volume DO n = 1, m_index+1 area = sum(aicetl(:)) - & ! total area of sub-layer (rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1) exit ELSE ! still in sub-layer below the sub-layer with the depth rem_vol = rem_vol - vol END IF END DO END IF END SUBROUTINE ice_thd_pnd_depth SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm) !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_pnd_perm *** !! !! ** Purpose : Determine the liquid fraction of brine in the ice !! and its permeability !!------------------------------------------------------------------- REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: & ticen, & ! internal ice temperature (K) salin ! salinity (ppt) !js: ppt according to cice REAL (wp), INTENT(OUT) :: & perm ! permeability REAL (wp) :: & Sbr ! brine salinity REAL (wp), DIMENSION(nlay_i) :: & Tin, & ! ice temperature phi ! liquid fraction INTEGER :: k !----------------------------------------------------------------- ! Compute ice temperatures from enthalpies using quadratic formula !----------------------------------------------------------------- DO k = 1,nlay_i Tin(k) = ticen(k) - rt0 !js: from K to degC END DO !----------------------------------------------------------------- ! brine salinity and liquid fraction !----------------------------------------------------------------- DO k = 1, nlay_i Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus) ! Best expression to date is that one ! Sbr = - 18.7 * Tin(k) − 0.519 * Tin(k)**2 − 0.00535 * Tin(k) ***3 phi(k) = salin(k) / Sbr END DO !----------------------------------------------------------------- ! permeability !----------------------------------------------------------------- perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007) END SUBROUTINE ice_thd_pnd_perm !---------------------------------------------------------------------------------------------------------------------- SUBROUTINE ice_thd_pnd_init !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_pnd_init *** !! !! ** Purpose : Physical constants and parameters linked to melt ponds !! over sea ice !! !! ** Method : Read the namthd_pnd namelist and check the melt pond !! parameter values called at the first timestep (nit000) !! !! ** input : Namelist namthd_pnd !!------------------------------------------------------------------- INTEGER :: ios, ioptio ! Local integer !! NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, & & ln_pnd_CST , rn_apnd, rn_hpnd, & & ln_pnd_TOPO , & & ln_pnd_lids, ln_pnd_alb !!------------------------------------------------------------------- ! REWIND( numnam_ice_ref ) ! Namelist namthd_pnd in reference namelist : Melt Ponds READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901) 901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in reference namelist' ) REWIND( numnam_ice_cfg ) ! Namelist namthd_pnd in configuration namelist : Melt Ponds READ ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 ) 902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' ) IF(lwm) WRITE ( numoni, namthd_pnd ) ! IF(lwp) THEN ! control print WRITE(numout,*) WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds' WRITE(numout,*) '~~~~~~~~~~~~~~~~' WRITE(numout,*) ' Namelist namicethd_pnd:' WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb ENDIF ! ! !== set the choice of ice pond scheme ==! ioptio = 0 IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF IF( ioptio /= 1 ) & & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_LEV or ln_pnd_CST)' ) ! SELECT CASE( nice_pnd ) CASE( np_pndNO ) IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF CASE( np_pndCST ) IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF END SELECT ! END SUBROUTINE ice_thd_pnd_init #else !!---------------------------------------------------------------------- !! Default option Empty module NO SI3 sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icethd_pnd