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 @ 14209

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

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

Download in other formats: