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 pond scheme INTEGER, PARAMETER :: np_pndH12 = 2 ! Evolutive pond scheme (Holland et al. 2012) REAL(wp), PUBLIC, PARAMETER :: pnd_lid_max = 0.015_wp ! pond lid thickness above which the ponds disappear from the albedo calculation REAL(wp), PUBLIC, PARAMETER :: pnd_lid_min = 0.005_wp ! pond lid thickness below which the full pond area is used in the albedo calculation REAL(wp), PARAMETER :: viscosity_dyn = 1.79e-3_wp !! * 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 !! !! ** Method : brut force !!------------------------------------------------------------------- ! SELECT CASE ( nice_pnd ) ! CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==! ! CASE (np_pndH12) ; CALL pnd_H12 !== Holland et al 2012 melt ponds ==! ! 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 a_ip_frac_1d(ji) = rn_apnd h_ip_1d(ji) = rn_hpnd a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) ELSE a_ip_frac_1d(ji) = 0._wp h_ip_1d(ji) = 0._wp a_ip_1d(ji) = 0._wp ENDIF ! END DO ! END SUBROUTINE pnd_CST SUBROUTINE pnd_H12 !!------------------------------------------------------------------- !! *** ROUTINE pnd_H12 *** !! !! ** Purpose : Compute melt pond evolution !! !! ** Method : Empirical method. A fraction of meltwater is accumulated in ponds !! and sent to ocean when surface is freezing !! !! pond growth: Vp = Vp + dVmelt !! with dVmelt = R/rhow * ( rhoi*dh_i + rhos*dh_s ) * a_i !! pond contraction: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) !! with Tp = -2degC !! !! ** Tunable parameters : (no real expertise yet, ideas?) !! !! ** Note : Stolen from CICE for quick test of the melt pond !! radiation and freshwater interfaces !! Coupling can be radiative AND freshwater !! Advection, ridging, rafting are called !! !! ** References : Holland, M. M. et al (J Clim 2012) !!------------------------------------------------------------------- REAL(wp), PARAMETER :: zrmin = 0.15_wp ! minimum fraction of available meltwater retained for melt ponding REAL(wp), PARAMETER :: zrmax = 0.70_wp ! maximum - - - - - REAL(wp), PARAMETER :: zpnd_aspect = 0.8_wp ! pond aspect ratio REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature REAL(wp), PARAMETER :: max_h_diff_s = -1.0E-6 ! Maximum meltpond depth change due to leaking or overflow (m s-1) ! REAL(wp) :: tot_mlt ! Total ice and snow surface melt (some goes into ponds, some into the ocean) REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding REAL(wp) :: zdv_mlt ! available meltwater for melt ponding (equivalent volume change) REAL(wp) :: z1_Tp ! inverse reference temperature REAL(wp) :: z1_rhow ! inverse freshwater density REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio REAL(wp) :: z1_rhoi ! inverse ice density REAL(wp) :: zfac, zdum REAL(wp) :: t_grad ! Temperature deficit for refreezing REAL(wp) :: omega_dt ! Time independent accumulated variables used for freezing REAL(wp) :: lh_ip_end ! Lid thickness at end of timestep (temporary variable) REAL(wp) :: zdh_frz ! Amount of melt pond that freezes (m) REAL(wp) :: v_ip_old ! Pond volume before leaking back to the ocean REAL(wp) :: dh_ip_over ! Pond thickness change due to overflow or leaking REAL(wp) :: dh_i_pndleak ! Grid box mean change in water depth due to leaking back to the ocean REAL(wp) :: h_gravity_head ! Height above sea level of the top of the melt pond REAL(wp) :: h_percolation ! Distance between the base of the melt pond and the base of the sea ice REAL(wp) :: Sbr ! Brine salinity REAL(wp), DIMENSION(nlay_i) :: phi ! liquid fraction REAL(wp) :: perm ! Permeability of the sea ice REAL(wp) :: za_ip ! Temporary pond fraction REAL(wp) :: max_h_diff_ts ! Maximum meltpond depth change due to leaking or overflow (m per ts) ! INTEGER :: ji, jk ! loop indices !!------------------------------------------------------------------- z1_rhow = 1._wp / rhow z1_zpnd_aspect = 1._wp / zpnd_aspect z1_Tp = 1._wp / zTp z1_rhoi = 1._wp / rhoi max_h_diff_ts = max_h_diff_s * rdt_ice ! Define time-independent field for use in refreezing omega_dt = 2.0_wp * rcnd_i * rdt_ice / (rLfus * rhow) 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 a_ip_frac_1d(ji) = 0._wp h_ip_1d(ji) = 0._wp lh_ip_1d(ji) = 0._wp ! !--------------------------------! ELSE ! Case ice thickness >= rn_himin ! ! !--------------------------------! v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! record pond volume at previous time step ! To avoid divide by zero errors in some calculations we will use a temporary pond fraction variable that has a minimum value of epsi06 za_ip = a_ip_1d(ji) IF ( za_ip < epsi06 ) za_ip = epsi06 ! ! available meltwater for melt ponding [m, >0] and fraction tot_mlt = -( dh_i_sum(ji)*rhoi + dh_s_mlt(ji)*rhos ) * z1_rhow IF ( ln_pnd_totfrac ) THEN zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * at_i_1d(ji) ! Use total ice fraction ELSE zfr_mlt = rn_pnd_min + ( rn_pnd_max - rn_pnd_min ) * a_i_1d(ji) ! Use category ice fraction ENDIF zdv_mlt = zfr_mlt * tot_mlt ! !--- Pond gowth ---! v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt ! !--- Lid shrinking. ---! IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) - zdv_mlt / za_ip ! ! melt pond mass flux (<0) IF( ln_use_pndmass .AND. zdv_mlt > 0._wp ) THEN zfac = zdv_mlt * rhow * r1_rdtice wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac ! ! adjust ice/snow melting flux to balance melt pond flux (>0) zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) 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 ! !--- Pond contraction (due to refreezing) ---! IF ( ln_pnd_lids ) THEN ! Code to use if using melt pond lids IF ( t_su_1d(ji) < (zTp+rt0) .AND. v_ip_1d(ji) > 0._wp ) THEN t_grad = (zTp+rt0) - t_su_1d(ji) ! The following equation is a rearranged form of: ! lid_thickness_end - lid_thickness_start = rcnd_i * t_grad * rdt_ice / (0.5*(lid_thickness_end + lid_thickness_start) * rLfus * rhow) ! where: lid_thickness_start = lh_ip_1d(ji) ! lid_thickness_end = lh_ip_end ! omega_dt is a bunch of terms in the equation that do not change ! Note the use of rhow instead of rhoi as we are working with volumes and it is mathematically easier ! if the water and ice specific volumes (for the lid and the pond) are the same (have the same density). lh_ip_end = SQRT(omega_dt * t_grad + lh_ip_1d(ji)**2) zdh_frz = lh_ip_end - lh_ip_1d(ji) ! Pond shrinking v_ip_1d(ji) = v_ip_1d(ji) - zdh_frz * a_ip_1d(ji) ! Lid growing IF ( ln_pnd_lids ) lh_ip_1d(ji) = lh_ip_1d(ji) + zdh_frz ELSE zdh_frz = 0._wp END IF ELSE ! Code to use if not using melt pond lids v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) ENDIF ! ! Make sure pond volume or lid thickness has not gone negative IF ( v_ip_1d(ji) < 0._wp ) v_ip_1d(ji) = 0._wp IF ( lh_ip_1d(ji) < 0._wp ) lh_ip_1d(ji) = 0._wp ! ! Set new pond area and depth assuming linear relation between h_ip and a_ip_frac ! h_ip = zpnd_aspect * a_ip_frac = zpnd_aspect * a_ip/a_i a_ip_1d(ji) = SQRT( v_ip_1d(ji) * z1_zpnd_aspect * a_i_1d(ji) ) a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) !--- Pond overflow and vertical flushing ---! IF ( ln_pnd_overflow ) THEN ! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond volume IF ( a_ip_1d(ji) > zfr_mlt * a_i_1d(ji) ) THEN v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use dh_ip_over = zpnd_aspect * zfr_mlt - h_ip_1d(ji) ! This will be a negative number dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ENDIF ! If pond depth exceeds half the ice thickness then reduce the pond volume IF ( h_ip_1d(ji) > 0.5_wp * h_i_1d(ji) ) THEN v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use dh_ip_over = 0.5_wp * h_i_1d(ji) - h_ip_1d(ji) ! This will be a negative number dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ENDIF !-- Vertical flushing of pond water --! ! The height above sea level of the top of the melt pond is the ratios of density of ice and water times the ice depth. ! This assumes the pond is sitting on top of the ice. h_gravity_head = h_i_1d(ji) * (rhow - rhoi) * z1_rhow ! The depth through which water percolates is the distance under the melt pond h_percolation = h_i_1d(ji) - h_ip_1d(ji) ! Calculate the permeability of the ice (Assur 1958) DO jk = 1, nlay_i Sbr = - 1.2_wp & - 21.8_wp * (t_i_1d(ji,jk)-rt0) & - 0.919_wp * (t_i_1d(ji,jk)-rt0)**2 & - 0.01878_wp * (t_i_1d(ji,jk)-rt0)**3 phi(jk) = sz_i_1d(ji,jk)/Sbr END DO perm = 3.0e-08_wp * (minval(phi))**3 ! Do the drainage using Darcy's law IF ( perm > 0._wp ) THEN dh_ip_over = -1.0_wp*perm*rhow*grav*h_gravity_head*rdt_ice / (viscosity_dyn*h_percolation) ! This should be a negative number dh_ip_over = MIN(dh_ip_over, 0._wp) ! Make sure it is negative v_ip_old = v_ip_1d(ji) ! Save original volume before leak for future use dh_ip_over = MAX(dh_ip_over,max_h_diff_ts) ! Apply a limit h_ip_1d(ji) = MAX(0._wp, h_ip_1d(ji) + dh_ip_over) a_ip_frac_1d(ji) = h_ip_1d(ji) / zpnd_aspect a_ip_1d(ji) = a_ip_frac_1d(ji) * a_i_1d(ji) v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ENDIF ENDIF ! If lid thickness is ten times greater than pond thickness then remove pond IF ( ln_pnd_lids ) THEN IF ( lh_ip_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN a_ip_1d(ji) = 0._wp a_ip_frac_1d(ji) = 0._wp h_ip_1d(ji) = 0._wp lh_ip_1d(ji) = 0._wp v_ip_1d(ji) = 0._wp ENDIF ENDIF ! If any of the previous changes has removed all the ice thickness then remove ice area. IF ( h_i_1d(ji) == 0._wp ) THEN a_i_1d(ji) = 0._wp h_s_1d(ji) = 0._wp ENDIF ! ENDIF END DO ! END SUBROUTINE pnd_H12 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_H12, ln_pnd_CST, rn_apnd, rn_hpnd, ln_pnd_alb, & rn_pnd_min, rn_pnd_max, ln_pnd_overflow, ln_pnd_lids, ln_pnd_totfrac, & ln_use_pndmass !!------------------------------------------------------------------- ! 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,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_pnd_min = ', rn_pnd_min WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_pnd_max = ', rn_pnd_max WRITE(numout,*) ' Use total ice fraction instead of category ice fraction ln_pnd_totfrac = ',ln_pnd_totfrac WRITE(numout,*) ' Allow ponds to overflow and have vertical flushing ln_pnd_overflow = ',ln_pnd_overflow WRITE(numout,*) ' Melt ponds can have frozen lids ln_pnd_lids = ',ln_pnd_lids WRITE(numout,*) ' Use melt pond mass flux diagnostic, passing to ocean ln_use_pndmass = ',ln_use_pndmass WRITE(numout,*) ' Prescribed melt pond fraction and depth 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,*) ' 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_H12 ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndH12 ; ENDIF IF( ioptio /= 1 ) & & CALL ctl_stop( 'ice_thd_pnd_init: choose either none (ln_pnd=F) or only one pond scheme (ln_pnd_H12 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 END SELECT ! END SUBROUTINE ice_thd_pnd_init #else !!---------------------------------------------------------------------- !! Default option Empty module NO SI3 sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icethd_pnd