MODULE icethd_pnd !!====================================================================== !! *** MODULE icethd_pnd *** !! Melt ponds !!====================================================================== !! history : ! Original code by Daniela Flocco and Adrian Turner !! 1.0 ! 2012 (O. Lecomte) Adaptation for scientific tests (NEMO3.1) !! 2.0 ! 2017 (M. Vancoppenolle, O. Lecomte, C. Rousset) Implementation in NEMO4 !!---------------------------------------------------------------------- #if defined key_lim3 !!---------------------------------------------------------------------- !! 'key_lim3' : ESIM 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) !! * Substitutions # include "vectopt_loop_substitute.h90" !!---------------------------------------------------------------------- !! NEMO/ICE 4.0 , NEMO Consortium (2017) !! $Id: icethd_pnd.F90 8420 2017-10-05 13:07:10Z clem $ !! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt) !!---------------------------------------------------------------------- CONTAINS SUBROUTINE ice_thd_pnd !!------------------------------------------------------------------- !! *** ROUTINE ice_thd_pnd *** !! !! ** Purpose : change melt pond fraction !! !! ** Method : brut force !! !! ** Action : - !! - !!------------------------------------------------------------------- IF( nn_timing == 1 ) CALL timing_start('ice_thd_pnd') 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 IF( nn_timing == 1 ) CALL timing_stop('ice_thd_pnd') 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 ! REAL(wp) :: zfr_mlt ! fraction of available meltwater retained for melt ponding REAL(wp) :: zdv_mlt ! available meltwater for melt ponding REAL(wp) :: z1_Tp ! reference temperature REAL(wp) :: z1_rhofw ! inverse freshwater density REAL(wp) :: z1_zpnd_aspect ! inverse pond aspect ratio REAL(wp) :: zvpnd ! dummy pond volume REAL(wp) :: zfac, zdum INTEGER :: ji ! loop indices !!------------------------------------------------------------------- z1_rhofw = 1. / rhofw z1_zpnd_aspect = 1. / zpnd_aspect z1_Tp = 1._wp / zTp DO ji = 1, npti ! !--------------------------------! IF( h_i_1d(ji) < rn_himin) THEN ! Case ice thickness < rn_himin ! ! !--------------------------------! !--- Remove ponds on thin ice and send water into the ocean IF( ln_pnd_fwb ) THEN wfx_pnd_1d(ji) = wfx_pnd_1d(ji) + h_ip_1d(ji) * a_ip_1d(ji) * rhofw * r1_rdtice ENDIF a_ip_1d(ji) = 0._wp a_ip_frac_1d(ji) = 0._wp h_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 zvpnd = v_ip_1d(ji) ! available meltwater for melt ponding [m, >0] and fraction zdv_mlt = -( dh_i_surf(ji)*rhoic + dh_s_mlt(ji)*rhosn ) * z1_rhofw * a_i_1d(ji) zfr_mlt = zrmin + ( zrmax - zrmin ) * a_i_1d(ji) !!clem CICE doc !!! zfr_mlt = zrmin + zrmax * a_i_1d(ji) !!clem Holland paper !--- Pond gowth ! v_ip should never be negative, otherwise code crashes ! MV: as far as I saw, UM5 can create very small negative v_ip values (not Prather) v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) + zfr_mlt * zdv_mlt ) IF( ln_pnd_fwb .AND. zdv_mlt > 0._wp ) THEN ! melt ponds mass flux (<0) zfac = zfr_mlt * zdv_mlt * rhofw * r1_rdtice wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac ! adjust ice/snow melting (>0) to take into account ponding 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) v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * MAX( zTp+rt0 - t_su_1d(ji), 0._wp ) * z1_Tp ) ! --- 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) ) !NB: the SQRT has been a recurring source of crash when v_ip or a_i tuns to be even only slightly negative a_ip_frac_1d(ji) = a_ip_1d(ji) / a_i_1d(ji) h_ip_1d(ji) = zpnd_aspect * a_ip_frac_1d(ji) !! clem: there is no clever way to conserve mass here... ! IF( ln_pnd_fwb ) THEN ! ! melt ponds mass flux (>0 when ponds shrink) ! zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw * r1_rdtice ! wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac!! ! ! ! adjust ice/snow ! zfac = ( v_ip_1d(ji) - zvpnd ) * rhofw / ( a_i_1d(ji)*h_s_1d(ji)*rhosn + a_i_1d(ji)*h_i_1d(ji)*rhoic )! ! ! wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + h_s_1d(ji) * zfac * a_i_1d(ji) * rhosn * r1_rdtice ! wfx_sum_1d(ji) = wfx_sum_1d(ji) + h_i_1d(ji) * zfac * a_i_1d(ji) * rhoic * r1_rdtice ! ! !!h_s_1d(ji) = h_s_1d(ji) * ( 1._wp + zfac ) ! !!h_i_1d(ji) = h_i_1d(ji) * ( 1._wp + zfac )! ! ! 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 output status for namelist read NAMELIST/namthd_pnd/ ln_pnd_H12, ln_pnd_fwb, ln_pnd_CST, rn_apnd, rn_hpnd, 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', lwp ) 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', lwp ) 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,*) ' Evolutive melt pond fraction and depth (Holland et al 2012) ln_pnd_H12 = ', ln_pnd_H12 WRITE(numout,*) ' Melt ponds store fresh water or not ln_pnd_fwb = ', ln_pnd_fwb 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 nice_pnd = np_pndNO 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 one and only one pond scheme (ln_pnd_H12 or ln_pnd_CST)' ) SELECT CASE( nice_pnd ) CASE( np_pndNO ) IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when no ponds' ) ; ENDIF IF(ln_pnd_alb) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF CASE( np_pndCST) IF(ln_pnd_fwb) THEN ; ln_pnd_fwb = .FALSE. ; CALL ctl_warn( 'ln_pnd_fwb=false when ln_pnd_CST=true' ) ; ENDIF END SELECT ! END SUBROUTINE ice_thd_pnd_init #else !!---------------------------------------------------------------------- !! Default option Empty module NO ESIM sea-ice model !!---------------------------------------------------------------------- #endif !!====================================================================== END MODULE icethd_pnd