New URL for NEMO forge! http://forge.nemo-ocean.eu

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

icethd_pnd.F90 in NEMO/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icethd_pnd.F90 @ 14686

Last change on this file since 14686 was 14252, checked in by vancop, 3 years ago
bug fixes in topo ponds for trunk
Property svn:keywords set to `Id`
File size: 65.7 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats: