Context Navigation

icethd_pnd.F90 @ 15696

Last change on this file since 15696 was 14669, checked in by dancopsey, 3 years ago
Add fixes to melt pond lids
File size: 65.1 KB

Rev	Line
[8637]	1	MODULE icethd_pnd
	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
[14659]	22	USE sbc_ice ! surface energy budget
[8637]	23	!
	24	USE in_out_manager ! I/O manager
[14659]	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:
[14659]	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
[14659]	43	REAL(wp), PARAMETER :: & ! shared parameters for topographic melt ponds
	44	zhi_min = 0.1_wp , & ! minimum ice thickness with ponds (m)
	45	zTd = 273.0_wp , & ! temperature difference for freeze-up (C)
	46	zvp_min = 1.e-4_wp ! minimum pond volume (m)
	47
	48	!--------------------------------------------------------------------------
	49	!
	50	! Pond volume per area budget diags
	51	!
	52	! The idea of diags is the volume of ponds per grid cell area is
	53	!
	54	! dV/dt = mlt + drn + lid + rnf
	55	! mlt = input from surface melting
	56	! drn = drainage through brine network
	57	! lid = lid growth & melt
	58	! rnf = runoff (water directly removed out of surface melting + overflow)
	59	!
	60	! In topo mode, the pond water lost because it is in the snow is not included in the budget
	61	!
	62	! In level mode, all terms are incorporated
	63
	64	REAL(wp), ALLOCATABLE, DIMENSION(:,:) :: & ! pond volume budget diagnostics
	65	diag_dvpn_mlt, & !! meltwater pond volume input [m/day]
	66	diag_dvpn_drn, & !! pond volume lost by drainage [m/day]
	67	diag_dvpn_lid, & !! exchange with lid / refreezing [m/day]
	68	diag_dvpn_rnf !! meltwater pond lost to runoff [m/day]
	69
	70	REAL(wp), ALLOCATABLE, DIMENSION(:) :: & ! pond volume budget diagnostics (1d)
	71	diag_dvpn_mlt_1d, & !! meltwater pond volume input [m/day]
	72	diag_dvpn_drn_1d, & !! pond volume lost by drainage [m/day]
	73	diag_dvpn_lid_1d, & !! exchange with lid / refreezing [m/day]
	74	diag_dvpn_rnf_1d !! meltwater pond lost to runoff [m/day]
	75	!
	76	!--------------------------------------------------------------------------
	77
[8637]	78	!! * Substitutions
	79	# include "vectopt_loop_substitute.h90"
	80	!!----------------------------------------------------------------------
[9598]	81	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
[10069]	82	!! $Id$
[10068]	83	!! Software governed by the CeCILL license (see ./LICENSE)
[8637]	84	!!----------------------------------------------------------------------
	85	CONTAINS
	86
	87	SUBROUTINE ice_thd_pnd
	88	!!-------------------------------------------------------------------
	89	!! * ROUTINE ice_thd_pnd *
	90	!!
[13284]	91	!! ** Purpose : change melt pond fraction and thickness
[14659]	92	!!
	93	!! Note: Melt ponds affect only radiative transfer for now
	94	!!
	95	!! No heat, no salt.
	96	!! The melt water they carry is collected but
	97	!! not removed from fw budget or released to the ocean
	98	!!
	99	!! A wfx_pnd has been coded for diagnostic purposes
	100	!! It is not fully consistent yet.
	101	!!
	102	!! The current diagnostic lacks a contribution from drainage
	103	!!
[8637]	104	!!-------------------------------------------------------------------
[14659]	105	!!
	106
	107	ALLOCATE( diag_dvpn_mlt(jpi,jpj), diag_dvpn_lid(jpi,jpj), diag_dvpn_drn(jpi,jpj), diag_dvpn_rnf(jpi,jpj) )
	108	ALLOCATE( diag_dvpn_mlt_1d(jpij), diag_dvpn_lid_1d(jpij), diag_dvpn_drn_1d(jpij), diag_dvpn_rnf_1d(jpij) )
	109
	110	diag_dvpn_mlt(:,:) = 0._wp ; diag_dvpn_drn(:,:) = 0._wp
	111	diag_dvpn_lid(:,:) = 0._wp ; diag_dvpn_rnf(:,:) = 0._wp
	112	diag_dvpn_mlt_1d(:) = 0._wp ; diag_dvpn_drn_1d(:) = 0._wp
	113	diag_dvpn_lid_1d(:) = 0._wp ; diag_dvpn_rnf_1d(:) = 0._wp
	114
[8637]	115	SELECT CASE ( nice_pnd )
[9169]	116	!
[14659]	117	CASE (np_pndCST) ; CALL pnd_CST !== Constant melt ponds ==!
[9169]	118	!
[14659]	119	CASE (np_pndLEV) ; CALL pnd_LEV !== Level ice melt ponds ==!
[9169]	120	!
[14659]	121	CASE (np_pndTOPO) ; CALL pnd_TOPO !== Topographic melt ponds ==!
	122	!
[8637]	123	END SELECT
[9169]	124	!
[14659]	125
	126	IF( iom_use('dvpn_mlt' ) ) CALL iom_put( 'dvpn_mlt', diag_dvpn_mlt * 100._wp * 86400._wp ) ! input from melting
	127	IF( iom_use('dvpn_lid' ) ) CALL iom_put( 'dvpn_lid', diag_dvpn_lid * 100._wp * 86400._wp ) ! exchanges with lid
	128	IF( iom_use('dvpn_drn' ) ) CALL iom_put( 'dvpn_drn', diag_dvpn_drn * 100._wp * 86400._wp ) ! vertical drainage
	129	IF( iom_use('dvpn_rnf' ) ) CALL iom_put( 'dvpn_rnf', diag_dvpn_rnf * 100._wp * 86400._wp ) ! runoff + overflow
	130
	131	DEALLOCATE( diag_dvpn_mlt, diag_dvpn_lid, diag_dvpn_drn, diag_dvpn_rnf )
	132	DEALLOCATE( diag_dvpn_mlt_1d, diag_dvpn_lid_1d, diag_dvpn_drn_1d, diag_dvpn_rnf_1d )
	133
[8637]	134	END SUBROUTINE ice_thd_pnd
	135
[9169]	136
[8637]	137	SUBROUTINE pnd_CST
	138	!!-------------------------------------------------------------------
	139	!! * ROUTINE pnd_CST *
	140	!!
[9169]	141	!! ** Purpose : Compute melt pond evolution
[8637]	142	!!
[9169]	143	!! ** Method : Melt pond fraction and thickness are prescribed
[9604]	144	!! to non-zero values when t_su = 0C
[8637]	145	!!
[14659]	146	!! ** Tunable parameters : Pond fraction (rn_apnd) & depth (rn_hpnd)
[8637]	147	!!
[9169]	148	!! ** Note : Coupling with such melt ponds is only radiative
[9604]	149	!! Advection, ridging, rafting... are bypassed
[8637]	150	!!
	151	!! ** References : Bush, G.W., and Trump, D.J. (2017)
	152	!!-------------------------------------------------------------------
	153	INTEGER :: ji ! loop indices
	154	!!-------------------------------------------------------------------
	155	DO ji = 1, npti
[9169]	156	!
[8637]	157	IF( a_i_1d(ji) > 0._wp .AND. t_su_1d(ji) >= rt0 ) THEN
	158	h_ip_1d(ji) = rn_hpnd
[13284]	159	a_ip_1d(ji) = rn_apnd * a_i_1d(ji)
	160	h_il_1d(ji) = 0._wp ! no pond lids whatsoever
[8637]	161	ELSE
	162	h_ip_1d(ji) = 0._wp
	163	a_ip_1d(ji) = 0._wp
[13284]	164	h_il_1d(ji) = 0._wp
[8637]	165	ENDIF
[9169]	166	!
[8637]	167	END DO
[9169]	168	!
[8637]	169	END SUBROUTINE pnd_CST
	170
[9169]	171
[13284]	172	SUBROUTINE pnd_LEV
[8637]	173	!!-------------------------------------------------------------------
[13284]	174	!! * ROUTINE pnd_LEV *
[8637]	175	!!
[13284]	176	!! ** Purpose : Compute melt pond evolution
[8637]	177	!!
[13284]	178	!! ** Method : A fraction of meltwater is accumulated in ponds and sent to ocean when surface is freezing
	179	!! We work with volumes and then redistribute changes into thickness and concentration
	180	!! assuming linear relationship between the two.
[8637]	181	!!
[13284]	182	!! ** Action : - pond growth: Vp = Vp + dVmelt --- from Holland et al 2012 ---
	183	!! dVmelt = (1-r)/rhow * ( rhoidh_i + rhosdh_s ) * a_i
	184	!! dh_i = meltwater from ice surface melt
	185	!! dh_s = meltwater from snow melt
	186	!! (1-r) = fraction of melt water that is not flushed
	187	!!
	188	!! - limtations: a_ip must not exceed (1-r)*a_i
	189	!! h_ip must not exceed 0.5*h_i
	190	!!
	191	!! - pond shrinking:
	192	!! if lids: Vp = Vp -dH * a_ip
	193	!! dH = lid thickness change. Retrieved from this eq.: --- from Flocco et al 2010 ---
	194	!!
	195	!! rhoi * Lf * dH/dt = ki * MAX(Tp-Tsu,0) / H
	196	!! H = lid thickness
	197	!! Lf = latent heat of fusion
	198	!! Tp = -2C
	199	!!
	200	!! And solved implicitely as:
	201	!! H(t+dt)*2 -H(t) H(t+dt) -ki * (Tp-Tsu) * dt / (rhoi*Lf) = 0
	202	!!
	203	!! if no lids: Vp = Vp * exp(0.01*MAX(Tp-Tsu,0)/Tp) --- from Holland et al 2012 ---
	204	!!
	205	!! - Flushing: w = -perm/visc * rho_oce * grav * Hp / Hi --- from Flocco et al 2007 ---
	206	!! perm = permability of sea-ice
	207	!! visc = water viscosity
	208	!! Hp = height of top of the pond above sea-level
	209	!! Hi = ice thickness thru which there is flushing
	210	!!
	211	!! - Corrections: remove melt ponds when lid thickness is 10 times the pond thickness
	212	!!
	213	!! - pond thickness and area is retrieved from pond volume assuming a linear relationship between h_ip and a_ip:
	214	!! a_ip/a_i = a_ip_frac = h_ip / zaspect
	215	!!
[14659]	216	!! ** Tunable parameters : rn_apnd_max, rn_apnd_min, rn_pnd_flush
[8637]	217	!!
[14659]	218	!! ** Note : Mostly stolen from CICE but not only
[8637]	219	!!
[14659]	220	!! ** References : Holland et al (J. Clim, 2012), Hunke et al (OM 2012)
	221	!!
[8637]	222	!!-------------------------------------------------------------------
[13284]	223	REAL(wp), DIMENSION(nlay_i) :: ztmp ! temporary array
	224	!!
	225	REAL(wp), PARAMETER :: zaspect = 0.8_wp ! pond aspect ratio
	226	REAL(wp), PARAMETER :: zTp = -2._wp ! reference temperature
	227	REAL(wp), PARAMETER :: zvisc = 1.79e-3_wp ! water viscosity
	228	!!
	229	REAL(wp) :: zfr_mlt, zdv_mlt ! fraction and volume of available meltwater retained for melt ponding
	230	REAL(wp) :: zdv_frz, zdv_flush ! Amount of melt pond that freezes, flushes
	231	REAL(wp) :: zhp ! heigh of top of pond lid wrt ssh
	232	REAL(wp) :: zv_ip_max ! max pond volume allowed
	233	REAL(wp) :: zdT ! zTp-t_su
	234	REAL(wp) :: zsbr ! Brine salinity
	235	REAL(wp) :: zperm ! permeability of sea ice
	236	REAL(wp) :: zfac, zdum ! temporary arrays
	237	REAL(wp) :: z1_rhow, z1_aspect, z1_Tp ! inverse
[14659]	238	REAL(wp) :: zvold !
[13284]	239	!!
[14659]	240	INTEGER :: ji, jj, jk, jl ! loop indices
	241
[8637]	242	!!-------------------------------------------------------------------
[13284]	243	z1_rhow = 1._wp / rhow
	244	z1_aspect = 1._wp / zaspect
	245	z1_Tp = 1._wp / zTp
[14659]	246
	247	!-----------------------------------------------------------------------------------------------
	248	! Identify grid cells with ice
	249	!-----------------------------------------------------------------------------------------------
	250	at_i(:,:) = SUM( a_i, dim=3 )
	251	!
	252	npti = 0 ; nptidx(:) = 0
	253	DO jj = 1, jpj
	254	DO ji = 1, jpi
	255	IF ( at_i(ji,jj) > epsi10 ) THEN
	256	npti = npti + 1
	257	nptidx( npti ) = (jj - 1) * jpi + ji
	258	ENDIF
	259	END DO
	260	END DO
	261
	262	!-----------------------------------------------------------------------------------------------
	263	! Prepare 1D arrays
	264	!-----------------------------------------------------------------------------------------------
[8637]	265
[14659]	266	IF( npti > 0 ) THEN
	267
	268	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
	269	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_sum_1d (1:npti) , wfx_sum )
	270	CALL tab_2d_1d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti) , wfx_pnd )
	271
	272	CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt )
	273	CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn )
	274	CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid )
	275	CALL tab_2d_1d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf )
	276
	277	DO jl = 1, jpl
	278
	279	CALL tab_2d_1d( npti, nptidx(1:npti), a_i_1d (1:npti), a_i (:,:,jl) )
	280	CALL tab_2d_1d( npti, nptidx(1:npti), h_i_1d (1:npti), h_i (:,:,jl) )
	281	CALL tab_2d_1d( npti, nptidx(1:npti), t_su_1d (1:npti), t_su (:,:,jl) )
	282	CALL tab_2d_1d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) )
	283	CALL tab_2d_1d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) )
	284	CALL tab_2d_1d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) )
[13284]	285
[14659]	286	CALL tab_2d_1d( npti, nptidx(1:npti), dh_i_sum (1:npti), dh_i_sum_2d (:,:,jl) )
	287	CALL tab_2d_1d( npti, nptidx(1:npti), dh_s_mlt (1:npti), dh_s_mlt_2d (:,:,jl) )
[13284]	288
[14659]	289	DO jk = 1, nlay_i
	290	CALL tab_2d_1d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) )
	291	END DO
	292
	293	!-----------------------------------------------------------------------------------------------
	294	! Go for ponds
	295	!-----------------------------------------------------------------------------------------------
[13284]	296
[14659]	297	DO ji = 1, npti
	298	! !----------------------------------------------------!
	299	IF( h_i_1d(ji) < rn_himin .OR. a_i_1d(ji) < epsi10 ) THEN ! Case ice thickness < rn_himin or tiny ice fraction !
	300	! !----------------------------------------------------!
	301	!--- Remove ponds on thin ice or tiny ice fractions
	302	a_ip_1d(ji) = 0._wp
	303	h_ip_1d(ji) = 0._wp
	304	h_il_1d(ji) = 0._wp
	305	! !--------------------------------!
	306	ELSE ! Case ice thickness >= rn_himin !
	307	! !--------------------------------!
	308	v_ip_1d(ji) = h_ip_1d(ji) * a_ip_1d(ji) ! retrieve volume from thickness
	309	v_il_1d(ji) = h_il_1d(ji) * a_ip_1d(ji)
	310	!
	311	!------------------!
	312	! case ice melting !
	313	!------------------!
	314	!
	315	!--- available meltwater for melt ponding ---!
	316	zdum = -( dh_i_sum(ji)rhoi + dh_s_mlt(ji)rhos ) * z1_rhow * a_i_1d(ji)
	317	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
	318	zdv_mlt = MAX( 0._wp, zfr_mlt * zdum ) ! max for roundoff errors?
	319
	320	diag_dvpn_mlt_1d(ji) = diag_dvpn_mlt_1d(ji) + zdum * r1_rdtice ! surface melt input diag
	321	diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( 1. - zfr_mlt ) * zdum * r1_rdtice ! runoff diag
	322	!
	323	!--- overflow ---!
	324	!
	325	! 1) area driven overflow
	326	!
	327	! If pond area exceeds zfr_mlt * a_i_1d(ji) then reduce the pond water volume
	328	! a_ip_max = zfr_mlt * a_i
	329	! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
	330	zv_ip_max = zfr_mlt*2 a_i_1d(ji) * zaspect
	331	zvold = zdv_mlt
	332	zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
	333
	334	!
	335	! 2) depth driven overflow
	336	!
	337	! If pond depth exceeds half the ice thickness then reduce the pond volume
	338	! h_ip_max = 0.5 * h_i
	339	! => from zaspect = h_ip / (a_ip / a_i), set v_ip_max as:
	340	zv_ip_max = z1_aspect * a_i_1d(ji) * 0.25 * h_i_1d(ji) * h_i_1d(ji) ! MV dimensions are wrong here or comment is unclear
	341	zdv_mlt = MAX( 0._wp, MIN( zdv_mlt, zv_ip_max - v_ip_1d(ji) ) )
	342	diag_dvpn_rnf_1d(ji) = diag_dvpn_rnf_1d(ji) + ( zdv_mlt - zvold ) * r1_rdtice ! runoff diag - overflow contribution
	343
	344	!--- Pond growing ---!
	345	v_ip_1d(ji) = v_ip_1d(ji) + zdv_mlt
	346	!
	347	!--- Lid melting ---!
	348	IF( ln_pnd_lids ) v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) - zdv_mlt ) ! must be bounded by 0
	349	!
	350	!--- mass flux ---!
	351	! MV I would recommend to remove that
	352	! Since melt ponds carry no freshwater there is no point in modifying water fluxes
	353
	354	IF( zdv_mlt > 0._wp ) THEN
	355	zfac = zdv_mlt * rhow * r1_rdtice ! melt pond mass flux < 0 [kg.m-2.s-1]
	356	wfx_pnd_1d(ji) = wfx_pnd_1d(ji) - zfac
	357	!
	358	zdum = zfac / ( wfx_snw_sum_1d(ji) + wfx_sum_1d(ji) ) ! adjust ice/snow melting flux > 0 to balance melt pond flux
	359	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) * (1._wp + zdum)
	360	wfx_sum_1d(ji) = wfx_sum_1d(ji) * (1._wp + zdum)
	361	ENDIF
	362
	363	!-------------------!
	364	! case ice freezing ! i.e. t_su_1d(ji) < (zTp+rt0)
	365	!-------------------!
	366	!
	367	zdT = MAX( zTp+rt0 - t_su_1d(ji), 0._wp )
	368	!
	369	!--- Pond contraction (due to refreezing) ---!
	370	zvold = v_ip_1d(ji) ! for diag
	371
	372	IF( ln_pnd_lids ) THEN
	373	!
	374	!--- Lid growing and subsequent pond shrinking ---!
	375	zdv_frz = 0.5_wp * MAX( 0._wp, -v_il_1d(ji) + & ! Flocco 2010 (eq. 5) solved implicitly as aH**2 + bH + c = 0
	376	& SQRT( v_il_1d(ji)2 + a_ip_1d(ji)2 * 4._wp * rcnd_i * zdT * rdt_ice / (rLfus * rhow) ) ) ! max for roundoff errors
[13284]	377
[14659]	378	! Lid growing
	379	v_il_1d(ji) = MAX( 0._wp, v_il_1d(ji) + zdv_frz )
[13284]	380
[14659]	381	! Pond shrinking
	382	v_ip_1d(ji) = MAX( 0._wp, v_ip_1d(ji) - zdv_frz )
	383
	384	ELSE
	385
	386	! Pond shrinking
	387	v_ip_1d(ji) = v_ip_1d(ji) * EXP( 0.01_wp * zdT * z1_Tp ) ! Holland 2012 (eq. 6)
	388
	389	ENDIF
[13284]	390
[14659]	391	diag_dvpn_lid_1d(ji) = diag_dvpn_lid_1d(ji) + ( v_ip_1d(ji) - zvold ) * r1_rdtice ! shrinking counted as lid diagnostic
[13284]	392
[14659]	393	!
	394	!--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
	395	! v_ip = h_ip * a_ip
	396	! 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)
	397	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
	398	h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
	399
	400	!------------------------------------------------!
	401	! Pond drainage through brine network (flushing) !
	402	!------------------------------------------------!
	403	! height of top of the pond above sea-level
	404	zhp = ( h_i_1d(ji) * ( rau0 - rhoi ) + h_ip_1d(ji) * ( rau0 - rhow * a_ip_1d(ji) / a_i_1d(ji) ) ) * r1_rau0
[13284]	405
[14659]	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	! MV linear expression more consistent & simpler zsbr = - ( t_i_1d(ji,jk) - rt0 ) / rTmlt
	414	ztmp(jk) = sz_i_1d(ji,jk) / zsbr
	415	END DO
	416	zperm = MAX( 0._wp, 3.e-08_wp * MINVAL(ztmp)**3 )
[13284]	417
[14659]	418	! Do the drainage using Darcy's law
	419	zdv_flush = - zperm * rau0 * grav * zhp * rdt_ice / (zvisc * h_i_1d(ji)) * a_ip_1d(ji) * rn_pnd_flush
	420	zdv_flush = MAX( zdv_flush, -v_ip_1d(ji) )
	421	! zdv_flush = 0._wp ! MV remove pond drainage for now
	422	v_ip_1d(ji) = v_ip_1d(ji) + zdv_flush
	423
	424	diag_dvpn_drn_1d(ji) = diag_dvpn_drn_1d(ji) + zdv_flush * r1_rdtice ! shrinking counted as lid diagnostic
[13284]	425
[14659]	426	! MV --- why pond drainage does not give back water into freshwater flux ?
	427	!--- Set new pond area and depth ---! assuming linear relation between h_ip and a_ip_frac
	428	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
	429	h_ip_1d(ji) = zaspect * a_ip_1d(ji) / a_i_1d(ji)
[13284]	430
[14659]	431	!--- Corrections and lid thickness ---!
	432	IF( ln_pnd_lids ) THEN
	433	!--- retrieve lid thickness from volume ---!
	434	IF( a_ip_1d(ji) > epsi10 ) THEN ; h_il_1d(ji) = v_il_1d(ji) / a_ip_1d(ji)
	435	ELSE ; h_il_1d(ji) = 0._wp
	436	ENDIF
	437	!--- remove ponds if lids are much larger than ponds ---!
	438	IF ( h_il_1d(ji) > h_ip_1d(ji) * 10._wp ) THEN
	439	a_ip_1d(ji) = 0._wp
	440	h_ip_1d(ji) = 0._wp
	441	h_il_1d(ji) = 0._wp
	442	ENDIF
	443	ENDIF
	444	!
[13284]	445	ENDIF
[14659]	446
	447	END DO ! ji
	448
	449	!-----------------------------------------------------------------------------------------------
	450	! Retrieve 2D arrays
	451	!-----------------------------------------------------------------------------------------------
	452
	453	v_ip_1d(1:npti) = h_ip_1d(1:npti) * a_ip_1d(1:npti)
	454	v_il_1d(1:npti) = h_il_1d(1:npti) * a_ip_1d(1:npti)
	455
	456	CALL tab_1d_2d( npti, nptidx(1:npti), a_ip_1d (1:npti), a_ip (:,:,jl) )
	457	CALL tab_1d_2d( npti, nptidx(1:npti), h_ip_1d (1:npti), h_ip (:,:,jl) )
	458	CALL tab_1d_2d( npti, nptidx(1:npti), h_il_1d (1:npti), h_il (:,:,jl) )
	459	CALL tab_1d_2d( npti, nptidx(1:npti), v_ip_1d (1:npti), v_ip (:,:,jl) )
	460	CALL tab_1d_2d( npti, nptidx(1:npti), v_il_1d (1:npti), v_il (:,:,jl) )
	461	DO jk = 1, nlay_i
	462	CALL tab_1d_2d( npti, nptidx(1:npti), sz_i_1d(1:npti,jk), sz_i(:,:,jk,jl) )
	463	END DO
	464
	465	END DO ! ji
	466
	467	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_snw_sum_1d(1:npti), wfx_snw_sum )
	468	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_sum_1d (1:npti), wfx_sum )
	469	CALL tab_1d_2d( npti, nptidx(1:npti), wfx_pnd_1d (1:npti), wfx_pnd )
[13284]	470
[14659]	471	CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_mlt_1d (1:npti), diag_dvpn_mlt )
	472	CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_drn_1d (1:npti), diag_dvpn_drn )
	473	CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_lid_1d (1:npti), diag_dvpn_lid )
	474	CALL tab_1d_2d( npti, nptidx(1:npti), diag_dvpn_rnf_1d (1:npti), diag_dvpn_rnf )
	475
[9169]	476	!
[14659]	477	ENDIF
	478
[13284]	479	END SUBROUTINE pnd_LEV
[14659]	480
	481	SUBROUTINE pnd_TOPO
	482
	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
	490	!! (2007) and Flocco et al. (2010).
	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	!!-------------------------------------------------------------------
[8637]	512
[14659]	513	! local variables
	514	REAL(wp) :: &
	515	zdHui, & ! change in thickness of ice lid (m)
	516	zomega, & ! conduction
	517	zdTice, & ! temperature difference across ice lid (C)
	518	zdvice, & ! change in ice volume (m)
	519	zTavg, & ! mean surface temperature across categories (C)
	520	zfsurf, & ! net heat flux, excluding conduction and transmitted radiation (W/m2)
	521	zTp, & ! pond freezing temperature (C)
[14669]	522	lid_Ttop,& ! lid top temperature (K)
[14659]	523	zrhoi_L, & ! volumetric latent heat of sea ice (J/m^3)
	524	zfr_mlt, & ! fraction and volume of available meltwater retained for melt ponding
	525	z1_rhow, & ! inverse water density
	526	zv_mlt ! total amount of meltwater produced
[9169]	527
[14659]	528	REAL(wp), DIMENSION(jpi,jpj) :: zvolp_ini, & !! total melt pond water available before redistribution and drainage
	529	zvolp , & !! total melt pond water volume
	530	zvolp_res !! remaining melt pond water available after drainage
	531
[14669]	532	REAL(wp), DIMENSION(jpi,jpj,jpl) :: z1_a_i, &
	533	zv_pnd !! volume of meltwater contributing to ponds
[14659]	534
	535	INTEGER :: ji, jj, jk, jl ! loop indices
	536
	537	INTEGER :: i_test
	538
	539	! Note
	540	! equivalent for CICE translation
	541	! a_ip -> apond
	542	! a_ip_frac -> apnd
	543
	544	!---------------------------------------------------------------
	545	! Initialise
	546	!---------------------------------------------------------------
	547
	548	! Parameters & constants (move to parameters)
	549	zrhoi_L = rhoi * rLfus ! volumetric latent heat (J/m^3)
	550	zTp = rt0 - 0.15_wp ! pond freezing point, slightly below 0C (ponds are bid saline)
	551	z1_rhow = 1._wp / rhow
	552
	553	! Set required ice variables (hard-coded here for now)
	554	! zfpond(:,:) = 0._wp ! contributing freshwater flux (?)
	555
	556	at_i (:,:) = SUM( a_i (:,:,:), dim=3 ) ! ice fraction
	557	vt_i (:,:) = SUM( v_i (:,:,:), dim=3 ) ! volume per grid area
	558	vt_ip(:,:) = SUM( v_ip(:,:,:), dim=3 ) ! pond volume per grid area
	559	vt_il(:,:) = SUM( v_il(:,:,:), dim=3 ) ! lid volume per grid area
	560
	561	! thickness
	562	WHERE( a_i(:,:,:) > epsi20 ) ; z1_a_i(:,:,:) = 1._wp / a_i(:,:,:)
	563	ELSEWHERE ; z1_a_i(:,:,:) = 0._wp
	564	END WHERE
	565	h_i(:,:,:) = v_i (:,:,:) * z1_a_i(:,:,:)
	566
	567	!---------------------------------------------------------------
	568	! Change 2D to 1D
	569	!---------------------------------------------------------------
	570	! MV
	571	! a less computing-intensive version would have 2D-1D passage here
	572	! use what we have in iceitd.F90 (incremental remapping)
	573
	574	!--------------------------------------------------------------
	575	! Collect total available pond water volume
	576	!--------------------------------------------------------------
	577	! Assuming that meltwater (+rain in principle) runsoff the surface
	578	! Holland et al (2012) suggest that the fraction of runoff decreases with total ice fraction
	579	! I cite her words, they are very talkative
	580	! "grid cells with very little ice cover (and hence more open water area)
	581	! have a higher runoff fraction to rep- resent the greater proximity of ice to open water."
	582	! "This results in the same runoff fraction r for each ice category within a grid cell"
	583
	584	zvolp(:,:) = 0._wp
	585
	586	DO jl = 1, jpl
	587	DO jj = 1, jpj
	588	DO ji = 1, jpi
	589
	590	IF ( a_i(ji,jj,jl) > epsi10 ) THEN
	591
	592	!--- Available and contributing meltwater for melt ponding ---!
	593	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
	594	& * z1_rhow * a_i(ji,jj,jl)
	595	! MV -> could move this directly in ice_thd_dh and get an array (ji,jj,jl) for surface melt water volume per grid area
	596	zfr_mlt = rn_apnd_min + ( rn_apnd_max - rn_apnd_min ) * at_i(ji,jj) ! fraction of surface meltwater going to ponds
[14669]	597	zv_pnd(ji,jj,jl) = zfr_mlt * zv_mlt ! contributing meltwater volume for category jl
[14659]	598
	599	diag_dvpn_mlt(ji,jj) = diag_dvpn_mlt(ji,jj) + zv_mlt * r1_rdtice ! diags
	600	diag_dvpn_rnf(ji,jj) = diag_dvpn_rnf(ji,jj) + ( 1. - zfr_mlt ) * zv_mlt * r1_rdtice
	601
	602	!--- Create possible new ponds
	603	! if pond does not exist, create new pond over full ice area
	604	!IF ( a_ip_frac(ji,jj,jl) < epsi10 ) THEN
	605	IF ( a_ip(ji,jj,jl) < epsi10 ) THEN
	606	a_ip(ji,jj,jl) = a_i(ji,jj,jl)
	607	a_ip_frac(ji,jj,jl) = 1.0_wp ! pond fraction of sea ice (apnd for CICE)
	608	ENDIF
	609
	610	!--- Deepen existing ponds with no change in pond fraction, before redistribution and drainage
[14669]	611	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) + zv_pnd(ji,jj,jl) ! use pond water to increase thickness
[14659]	612	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
	613
	614	!--- Total available pond water volume (pre-existing + newly produced)j
	615	zvolp(ji,jj) = zvolp(ji,jj) + v_ip(ji,jj,jl)
	616	! zfpond(ji,jj) = zfpond(ji,jj) + zpond * a_ip_frac(ji,jj,jl) ! useless for now
	617
	618	ENDIF ! a_i
	619
	620	END DO! jl
	621	END DO ! jj
	622	END DO ! ji
	623
	624	zvolp_ini(:,:) = zvolp(:,:)
	625
	626	!--------------------------------------------------------------
	627	! Redistribute and drain water from ponds
	628	!--------------------------------------------------------------
	629	CALL ice_thd_pnd_area( zvolp, zvolp_res )
	630
	631	!--------------------------------------------------------------
	632	! Melt pond lid growth and melt
	633	!--------------------------------------------------------------
	634
	635	IF( ln_pnd_lids ) THEN
	636
	637	DO jj = 1, jpj
	638	DO ji = 1, jpi
	639
	640	IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp_ini(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
	641
	642	!--------------------------
	643	! Pond lid growth and melt
	644	!--------------------------
	645	! Mean surface temperature
	646	zTavg = 0._wp
	647	DO jl = 1, jpl
	648	zTavg = zTavg + t_su(ji,jj,jl)*a_i(ji,jj,jl)
	649	END DO
	650	!! EWB: bug here - should be divided by at_i(ji,jj) but zTavg is not used anywhere
	651	zTavg = zTavg / a_i(ji,jj,jl) !!! could get a division by zero here
	652
	653	DO jl = 1, jpl-1
	654
	655
[14669]	656	!----------------------------------------------------------------
	657	! Lid melting: floating upper ice layer melts in whole or part
	658	!----------------------------------------------------------------
[14659]	659
[14669]	660	zdvice = MIN( zv_pnd(ji,jj,jl), v_il(ji,jj,jl) )
[14659]	661
[14669]	662	v_il (ji,jj,jl) = v_il (ji,jj,jl) - zdvice
[14659]	663
[14669]	664	IF ( v_il(ji,jj,jl) < epsi10 ) THEN
	665	! ice lid melted
	666	v_il(ji,jj,jl) = 0._wp
	667	ENDIF
[14659]	668
[14669]	669	diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) + zdvice ! diag
	670
[14659]	671	!----------------------------------------------------------------
[14669]	672	! Freeze lid
[14659]	673	!----------------------------------------------------------------
	674
[14669]	675	! The temperature of the lid's top surface is warmer than the ice top temperature
	676	! by about 2 degrees (warmed by the melt pond under the lid)
	677	lid_Ttop = t_su(ji,jj,jl) + 2.0_wp
[14659]	678
[14669]	679	! differential growth of base of surface floating ice layer
	680	zdTice = MAX( - ( lid_Ttop - zTp ) , 0._wp ) ! > 0
	681	zomega = rcnd_i * zdTice / zrhoi_L
	682	zdvice = SQRT( 2._wp * zomega * rdt_ice * a_ip(ji,jj,jl)2 + v_il(ji,jj,jl)2 ) &
	683	- v_il(ji,jj,jl)
	684	zdvice = min( zdvice , v_ip(ji,jj,jl) )
[14659]	685
[14669]	686	IF ( zdvice > epsi10 ) THEN
	687	v_il (ji,jj,jl) = v_il(ji,jj,jl) + zdvice
	688	v_ip(ji,jj,jl) = v_ip(ji,jj,jl) - zdvice
[14659]	689
[14669]	690	! Possibly this check might not be needed to avoid divide by zero errors as at
	691	! the end of ice_thd_pnd_area all ponds with an area less than epsi10 were zerod anyway
	692	IF ( a_ip(ji,jj,jl) > epsi10 ) THEN
	693	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
	694	ENDIF
[14659]	695
[14669]	696	diag_dvpn_lid(ji,jj) = diag_dvpn_lid(ji,jj) - zdvice ! diag
	697
	698	ENDIF
	699
	700	!----------------------------------------------------------------
	701	! Remove ponds if lids are much larger than ponds
	702	!----------------------------------------------------------------
	703
	704	IF ( a_ip(ji,jj,jl) > epsi10) THEN
	705	h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl)
	706	IF ( h_il(ji,jj,jl) > h_ip(ji,jj,jl) * 10._wp ) THEN
	707	a_ip(ji,jj,jl) = 0._wp
	708	h_ip(ji,jj,jl) = 0._wp
	709	v_ip(ji,jj,jl) = 0._wp
	710	h_il(ji,jj,jl) = 0._wp
	711	v_il(ji,jj,jl) = 0._wp
[14659]	712	ENDIF
[14669]	713	ENDIF
[14659]	714
	715	END DO ! jl
	716
	717	ELSE ! remove ponds on thin ice
	718
	719	v_ip(ji,jj,:) = 0._wp
	720	v_il(ji,jj,:) = 0._wp
	721	! zfpond(ji,jj) = zfpond(ji,jj)- zvolp(ji,jj)
	722	! zvolp(ji,jj) = 0._wp
	723
	724	ENDIF
	725
	726	END DO ! jj
	727	END DO ! ji
	728
	729	ENDIF ! ln_pnd_lids
	730
	731	!---------------------------------------------------------------
	732	! Clean-up variables (probably duplicates what icevar would do)
	733	!---------------------------------------------------------------
	734	! MV comment
	735	! In the ideal world, the lines above should update only v_ip, a_ip, v_il
	736	! icevar should recompute all other variables (if needed at all)
	737
	738	DO jl = 1, jpl
	739
	740	DO jj = 1, jpj
	741	DO ji = 1, jpi
	742
	743	! ! zap lids on small ponds
	744	! IF ( a_i(ji,jj,jl) > epsi10 .AND. v_ip(ji,jj,jl) < epsi10 &
	745	! .AND. v_il(ji,jj,jl) > epsi10) THEN
	746	! v_il(ji,jj,jl) = 0._wp ! probably uselesss now since we get zap_small
	747	! ENDIF
	748
	749	! recalculate equivalent pond variables
	750	IF ( a_ip(ji,jj,jl) > epsi10) THEN
	751	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
	752	a_ip_frac(ji,jj,jl) = a_ip(ji,jj,jl) / a_i(ji,jj,jl) ! MV in principle, useless as computed in icevar
	753	h_il(ji,jj,jl) = v_il(ji,jj,jl) / a_ip(ji,jj,jl) ! MV in principle, useless as computed in icevar
	754	ENDIF
	755	! h_ip(ji,jj,jl) = 0._wp ! MV in principle, useless as computed in icevar
	756	! h_il(ji,jj,jl) = 0._wp ! MV in principle, useless as omputed in icevar
	757	! ENDIF
	758
	759	END DO ! ji
	760	END DO ! jj
	761
	762	END DO ! jl
	763
	764	END SUBROUTINE pnd_TOPO
	765
	766
	767	SUBROUTINE ice_thd_pnd_area( zvolp , zdvolp )
	768
	769	!!-------------------------------------------------------------------
	770	!! * ROUTINE ice_thd_pnd_area *
	771	!!
	772	!! ** Purpose : Given the total volume of available pond water,
	773	!! redistribute and drain water
	774	!!
	775	!! ** Method
	776	!!
	777	!-----------\|
	778	! \|
	779	! \|-----------\|
	780	!___________\|___________\|______________________________________sea-level
	781	! \| \|
	782	! \| \|---^--------\|
	783	! \| \| \| \|
	784	! \| \| \| \|-----------\| \|-------
	785	! \| \| \| alfan \| \| \|
	786	! \| \| \| \| \|--------------\|
	787	! \| \| \| \| \| \|
	788	!---------------------------v-------------------------------------------
	789	! \| \| ^ \| \| \|
	790	! \| \| \| \| \|--------------\|
	791	! \| \| \| betan \| \| \|
	792	! \| \| \| \|-----------\| \|-------
	793	! \| \| \| \|
	794	! \| \|---v------- \|
	795	! \| \|
	796	! \|-----------\|
	797	! \|
	798	!-----------\|
	799	!
	800	!!
	801	!!------------------------------------------------------------------
	802
	803	REAL (wp), DIMENSION(jpi,jpj), INTENT(INOUT) :: &
	804	zvolp, & ! total available pond water
	805	zdvolp ! remaining meltwater after redistribution
	806
	807	INTEGER :: &
	808	ns, &
	809	m_index, &
	810	permflag
	811
	812	REAL (wp), DIMENSION(jpl) :: &
	813	hicen, &
	814	hsnon, &
	815	asnon, &
	816	alfan, &
	817	betan, &
	818	cum_max_vol, &
	819	reduced_aicen
	820
	821	REAL (wp), DIMENSION(0:jpl) :: &
	822	cum_max_vol_tmp
	823
	824	REAL (wp) :: &
	825	hpond, &
	826	drain, &
	827	floe_weight, &
	828	pressure_head, &
	829	hsl_rel, &
	830	deltah, &
	831	perm, &
	832	msno
	833
	834	REAL (wp), parameter :: &
	835	viscosity = 1.79e-3_wp ! kinematic water viscosity in kg/m/s
	836
	837	INTEGER :: ji, jj, jk, jl ! loop indices
	838
	839	a_ip(:,:,:) = 0._wp
	840	h_ip(:,:,:) = 0._wp
	841
	842	DO jj = 1, jpj
	843	DO ji = 1, jpi
	844
	845	IF ( at_i(ji,jj) > 0.01 .AND. hm_i(ji,jj) > zhi_min .AND. zvolp(ji,jj) > zvp_min * at_i(ji,jj) ) THEN
	846
	847	!-------------------------------------------------------------------
	848	! initialize
	849	!-------------------------------------------------------------------
	850
	851	DO jl = 1, jpl
	852
	853	!----------------------------------------
	854	! compute the effective snow fraction
	855	!----------------------------------------
	856
	857	IF (a_i(ji,jj,jl) < epsi10) THEN
	858	hicen(jl) = 0._wp
	859	hsnon(jl) = 0._wp
	860	reduced_aicen(jl) = 0._wp
	861	asnon(jl) = 0._wp !js: in CICE 5.1.2: make sense as the compiler may not initiate the variables
	862	ELSE
	863	hicen(jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
	864	hsnon(jl) = v_s(ji,jj,jl) / a_i(ji,jj,jl)
	865	reduced_aicen(jl) = 1._wp ! n=jpl
	866
	867	!js: initial code in NEMO_DEV
	868	!IF (n < jpl) reduced_aicen(jl) = aicen(jl) &
	869	! * (-0.024_wp*hicen(jl) + 0.832_wp)
	870
	871	!js: from CICE 5.1.2: this limit reduced_aicen to 0.2 when hicen is too large
	872	IF (jl < jpl) reduced_aicen(jl) = a_i(ji,jj,jl) &
	873	* max(0.2_wp,(-0.024_wp*hicen(jl) + 0.832_wp))
	874
	875	asnon(jl) = reduced_aicen(jl) ! effective snow fraction (empirical)
	876	! MV should check whether this makes sense to have the same effective snow fraction in here
	877	! OLI: it probably doesn't
	878	END IF
	879
	880	! This choice for alfa and beta ignores hydrostatic equilibium of categories.
	881	! Hydrostatic equilibium of the entire ITD is accounted for below, assuming
	882	! a surface topography implied by alfa=0.6 and beta=0.4, and rigidity across all
	883	! categories. alfa and beta partition the ITD - they are areas not thicknesses!
	884	! Multiplying by hicen, alfan and betan (below) are thus volumes per unit area.
	885	! Here, alfa = 60% of the ice area (and since hice is constant in a category,
	886	! alfan = 60% of the ice volume) in each category lies above the reference line,
	887	! and 40% below. Note: p6 is an arbitrary choice, but alfa+beta=1 is required.
	888
	889	! MV:
	890	! Note that this choice is not in the original FF07 paper and has been adopted in CICE
	891	! No reason why is explained in the doc, but I guess there is a reason. I'll try to investigate, maybe
	892
	893	! Where does that choice come from ? => OLI : Coz' Chuck Norris said so...
	894
	895	alfan(jl) = 0.6 * hicen(jl)
	896	betan(jl) = 0.4 * hicen(jl)
	897
	898	cum_max_vol(jl) = 0._wp
	899	cum_max_vol_tmp(jl) = 0._wp
	900
	901	END DO ! jpl
	902
	903	cum_max_vol_tmp(0) = 0._wp
	904	drain = 0._wp
	905	zdvolp(ji,jj) = 0._wp
	906
	907	!----------------------------------------------------------
	908	! Drain overflow water, update pond fraction and volume
	909	!----------------------------------------------------------
	910
	911	!--------------------------------------------------------------------------
	912	! the maximum amount of water that can be contained up to each ice category
	913	!--------------------------------------------------------------------------
	914	! If melt ponds are too deep to be sustainable given the ITD (OVERFLOW)
	915	! Then the excess volume cum_max_vol(jl) drains out of the system
	916	! It should be added to wfx_pnd_out
	917
	918	DO jl = 1, jpl-1 ! last category can not hold any volume
	919
	920	IF (alfan(jl+1) >= alfan(jl) .AND. alfan(jl+1) > 0._wp ) THEN
	921
	922	! total volume in level including snow
	923	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1) + &
	924	(alfan(jl+1) - alfan(jl)) * sum(reduced_aicen(1:jl))
	925
	926	! subtract snow solid volumes from lower categories in current level
	927	DO ns = 1, jl
	928	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl) &
	929	- rhos/rhow * & ! free air fraction that can be filled by water
	930	asnon(ns) * & ! effective areal fraction of snow in that category
	931	max(min(hsnon(ns)+alfan(ns)-alfan(jl), alfan(jl+1)-alfan(jl)), 0._wp)
	932	END DO
	933
	934	ELSE ! assume higher categories unoccupied
	935	cum_max_vol_tmp(jl) = cum_max_vol_tmp(jl-1)
	936	END IF
	937	!IF (cum_max_vol_tmp(jl) < z0) THEN
	938	! CALL abort_ice('negative melt pond volume')
	939	!END IF
	940	END DO
	941	cum_max_vol_tmp(jpl) = cum_max_vol_tmp(jpl-1) ! last category holds no volume
	942	cum_max_vol (1:jpl) = cum_max_vol_tmp(1:jpl)
	943
	944	!----------------------------------------------------------------
	945	! is there more meltwater than can be held in the floe?
	946	!----------------------------------------------------------------
	947	IF (zvolp(ji,jj) >= cum_max_vol(jpl)) THEN
	948	drain = zvolp(ji,jj) - cum_max_vol(jpl) + epsi10
	949	zvolp(ji,jj) = zvolp(ji,jj) - drain ! update meltwater volume available
	950
	951	diag_dvpn_rnf(ji,jj) = - drain ! diag - overflow counted in the runoff part (arbitrary choice)
	952
	953	zdvolp(ji,jj) = drain ! this is the drained water
	954	IF (zvolp(ji,jj) < epsi10) THEN
	955	zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
	956	zvolp(ji,jj) = 0._wp
	957	END IF
	958	END IF
	959
	960	! height and area corresponding to the remaining volume
	961	! routine leaves zvolp unchanged
	962	CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
	963
	964	DO jl = 1, m_index
	965	!h_ip(jl) = hpond - alfan(jl) + alfan(1) ! here oui choulde update
	966	! ! volume instead, no ?
	967	h_ip(ji,jj,jl) = max((hpond - alfan(jl) + alfan(1)), 0._wp) !js: from CICE 5.1.2
	968	a_ip(ji,jj,jl) = reduced_aicen(jl)
	969	! in practise, pond fraction depends on the empirical snow fraction
	970	! so in turn on ice thickness
	971	END DO
	972	!zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
	973
	974	!------------------------------------------------------------------------
	975	! Drainage through brine network (permeability)
	976	!------------------------------------------------------------------------
	977	!!! drainage due to ice permeability - Darcy's law
	978
	979	! sea water level
	980	msno = 0._wp
	981	DO jl = 1 , jpl
	982	msno = msno + v_s(ji,jj,jl) * rhos
	983	END DO
	984	floe_weight = ( msno + rhoivt_i(ji,jj) + rau0zvolp(ji,jj) ) / at_i(ji,jj)
	985	hsl_rel = floe_weight / rau0 &
	986	- ( ( sum(betan(:)*a_i(ji,jj,:)) / at_i(ji,jj) ) + alfan(1) )
	987
	988	deltah = hpond - hsl_rel
	989	pressure_head = grav * rau0 * max(deltah, 0._wp)
	990
	991	! drain if ice is permeable
	992	permflag = 0
	993
	994	IF (pressure_head > 0._wp) THEN
	995	DO jl = 1, jpl-1
	996	IF ( hicen(jl) /= 0._wp ) THEN
	997
	998	!IF (hicen(jl) > 0._wp) THEN !js: from CICE 5.1.2
	999
	1000	perm = 0._wp ! MV ugly dummy patch
	1001	CALL ice_thd_pnd_perm(t_i(ji,jj,:,jl), sz_i(ji,jj,:,jl), perm) ! bof
	1002	IF (perm > 0._wp) permflag = 1
	1003
	1004	drain = perma_ip(ji,jj,jl)pressure_head*rdt_ice / &
	1005	(viscosity*hicen(jl))
	1006	zdvolp(ji,jj) = zdvolp(ji,jj) + min(drain, zvolp(ji,jj))
	1007	zvolp(ji,jj) = max(zvolp(ji,jj) - drain, 0._wp)
	1008
	1009	diag_dvpn_drn(ji,jj) = - drain ! diag (could be better coded)
	1010
	1011	IF (zvolp(ji,jj) < epsi10) THEN
	1012	zdvolp(ji,jj) = zdvolp(ji,jj) + zvolp(ji,jj)
	1013	zvolp(ji,jj) = 0._wp
	1014	END IF
	1015	END IF
	1016	END DO
	1017
	1018	! adjust melt pond dimensions
	1019	IF (permflag > 0) THEN
	1020	! recompute pond depth
	1021	CALL ice_thd_pnd_depth(reduced_aicen, asnon, hsnon, alfan, zvolp(ji,jj), cum_max_vol, hpond, m_index)
	1022	DO jl = 1, m_index
	1023	h_ip(ji,jj,jl) = hpond - alfan(jl) + alfan(1)
	1024	a_ip(ji,jj,jl) = reduced_aicen(jl)
	1025	END DO
	1026	!zapond = sum(a_ip(1:m_index)) !js: from CICE 5.1.2; not in Icepack1.1.0-6-gac6195d
	1027	END IF
	1028	END IF ! pressure_head
	1029
	1030	!-------------------------------
	1031	! remove water from the snow
	1032	!-------------------------------
	1033	!------------------------------------------------------------------------
	1034	! total melt pond volume in category does not include snow volume
	1035	! snow in melt ponds is not melted
	1036	!------------------------------------------------------------------------
	1037
	1038	! MV here, it seems that we remove some meltwater from the ponds, but I can't really tell
	1039	! how much, so I did not diagnose it
	1040	! so if there is a problem here, nobody is going to see it...
	1041
	1042
	1043	! Calculate pond volume for lower categories
	1044	DO jl = 1,m_index-1
	1045	v_ip(ji,jj,jl) = a_ip(ji,jj,jl) * h_ip(ji,jj,jl) & ! what is not in the snow
	1046	- (rhos/rhow) * asnon(jl) * min(hsnon(jl), h_ip(ji,jj,jl))
	1047	END DO
	1048
	1049	! Calculate pond volume for highest category = remaining pond volume
	1050
	1051	! The following is completely unclear to Martin at least
	1052	! Could we redefine properly and recode in a more readable way ?
	1053
	1054	! m_index = last category with melt pond
	1055
	1056	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
	1057
	1058	IF (m_index > 1) THEN
	1059	IF (zvolp(ji,jj) > sum( v_ip(ji,jj,1:m_index-1))) THEN
	1060	v_ip(ji,jj,m_index) = zvolp(ji,jj) - sum(v_ip(ji,jj,1:m_index-1))
	1061	ELSE
	1062	v_ip(ji,jj,m_index) = 0._wp
	1063	h_ip(ji,jj,m_index) = 0._wp
	1064	a_ip(ji,jj,m_index) = 0._wp
	1065	! If remaining pond volume is negative reduce pond volume of
	1066	! lower category
	1067	IF ( zvolp(ji,jj) + epsi10 < SUM(v_ip(ji,jj,1:m_index-1))) &
	1068	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)
	1069	END IF
	1070	END IF
	1071
	1072	DO jl = 1,m_index
	1073	IF (a_ip(ji,jj,jl) > epsi10) THEN
	1074	h_ip(ji,jj,jl) = v_ip(ji,jj,jl) / a_ip(ji,jj,jl)
	1075	ELSE
	1076	zdvolp(ji,jj) = zdvolp(ji,jj) + v_ip(ji,jj,jl)
	1077	h_ip(ji,jj,jl) = 0._wp
	1078	v_ip(ji,jj,jl) = 0._wp
	1079	a_ip(ji,jj,jl) = 0._wp
	1080	END IF
	1081	END DO
	1082	DO jl = m_index+1, jpl
	1083	h_ip(ji,jj,jl) = 0._wp
	1084	a_ip(ji,jj,jl) = 0._wp
	1085	v_ip(ji,jj,jl) = 0._wp
	1086	END DO
	1087
	1088	ENDIF
	1089	END DO ! ji
	1090	END DO ! jj
	1091
	1092	END SUBROUTINE ice_thd_pnd_area
	1093
	1094
	1095	SUBROUTINE ice_thd_pnd_depth(aicen, asnon, hsnon, alfan, zvolp, cum_max_vol, hpond, m_index)
	1096	!!-------------------------------------------------------------------
	1097	!! * ROUTINE ice_thd_pnd_depth *
	1098	!!
	1099	!! ** Purpose : Compute melt pond depth
	1100	!!-------------------------------------------------------------------
	1101
	1102	REAL (wp), DIMENSION(jpl), INTENT(IN) :: &
	1103	aicen, &
	1104	asnon, &
	1105	hsnon, &
	1106	alfan, &
	1107	cum_max_vol
	1108
	1109	REAL (wp), INTENT(IN) :: &
	1110	zvolp
	1111
	1112	REAL (wp), INTENT(OUT) :: &
	1113	hpond
	1114
	1115	INTEGER, INTENT(OUT) :: &
	1116	m_index
	1117
	1118	INTEGER :: n, ns
	1119
	1120	REAL (wp), DIMENSION(0:jpl+1) :: &
	1121	hitl, &
	1122	aicetl
	1123
	1124	REAL (wp) :: &
	1125	rem_vol, &
	1126	area, &
	1127	vol, &
	1128	tmp, &
	1129	z0 = 0.0_wp
	1130
	1131	!----------------------------------------------------------------
	1132	! hpond is zero if zvolp is zero - have we fully drained?
	1133	!----------------------------------------------------------------
	1134
	1135	IF (zvolp < epsi10) THEN
	1136	hpond = z0
	1137	m_index = 0
	1138	ELSE
	1139
	1140	!----------------------------------------------------------------
	1141	! Calculate the category where water fills up to
	1142	!----------------------------------------------------------------
	1143
	1144	!----------\|
	1145	! \|
	1146	! \|
	1147	! \|----------\| -- --
	1148	!__________\|__________\|_________________________________________ ^
	1149	! \| \| rem_vol ^ \| Semi-filled
	1150	! \| \|----------\|-- -- -- - ---\|-- ---- -- -- --v layer
	1151	! \| \| \| \|
	1152	! \| \| \| \|hpond
	1153	! \| \| \|----------\| \| \|-------
	1154	! \| \| \| \| \| \|
	1155	! \| \| \| \|---v-----\|
	1156	! \| \| m_index \| \| \|
	1157	!-------------------------------------------------------------
	1158
	1159	m_index = 0 ! 1:m_index categories have water in them
	1160	DO n = 1, jpl
	1161	IF (zvolp <= cum_max_vol(n)) THEN
	1162	m_index = n
	1163	IF (n == 1) THEN
	1164	rem_vol = zvolp
	1165	ELSE
	1166	rem_vol = zvolp - cum_max_vol(n-1)
	1167	END IF
	1168	exit ! to break out of the loop
	1169	END IF
	1170	END DO
	1171	m_index = min(jpl-1, m_index)
	1172
	1173	!----------------------------------------------------------------
	1174	! semi-filled layer may have m_index different snow in it
	1175	!----------------------------------------------------------------
	1176
	1177	!----------------------------------------------------------- ^
	1178	! \| alfan(m_index+1)
	1179	! \|
	1180	!hitl(3)--> \|----------\| \|
	1181	!hitl(2)--> \|------------\| * * * * *\| \|
	1182	!hitl(1)--> \|----------\|* * * * * * \|* * * * * \| \|
	1183	!hitl(0)-->------------------------------------------------- \| ^
	1184	! various snow from lower categories \| \|alfa(m_index)
	1185
	1186	! hitl - heights of the snow layers from thinner and current categories
	1187	! aicetl - area of each snow depth in this layer
	1188
	1189	hitl(:) = z0
	1190	aicetl(:) = z0
	1191	DO n = 1, m_index
	1192	hitl(n) = max(min(hsnon(n) + alfan(n) - alfan(m_index), &
	1193	alfan(m_index+1) - alfan(m_index)), z0)
	1194	aicetl(n) = asnon(n)
	1195
	1196	aicetl(0) = aicetl(0) + (aicen(n) - asnon(n))
	1197	END DO
	1198
	1199	hitl(m_index+1) = alfan(m_index+1) - alfan(m_index)
	1200	aicetl(m_index+1) = z0
	1201
	1202	!----------------------------------------------------------------
	1203	! reorder array according to hitl
	1204	! snow heights not necessarily in height order
	1205	!----------------------------------------------------------------
	1206
	1207	DO ns = 1, m_index+1
	1208	DO n = 0, m_index - ns + 1
	1209	IF (hitl(n) > hitl(n+1)) THEN ! swap order
	1210	tmp = hitl(n)
	1211	hitl(n) = hitl(n+1)
	1212	hitl(n+1) = tmp
	1213	tmp = aicetl(n)
	1214	aicetl(n) = aicetl(n+1)
	1215	aicetl(n+1) = tmp
	1216	END IF
	1217	END DO
	1218	END DO
	1219
	1220	!----------------------------------------------------------------
	1221	! divide semi-filled layer into set of sublayers each vertically homogenous
	1222	!----------------------------------------------------------------
	1223
	1224	!hitl(3)----------------------------------------------------------------
	1225	! \| * * * * * * * *
	1226	! \|* * * * * * * * *
	1227	!hitl(2)----------------------------------------------------------------
	1228	! \| * * * * * * * * \| * * * * * * * *
	1229	! \|* * * * * * * * * \|* * * * * * * * *
	1230	!hitl(1)----------------------------------------------------------------
	1231	! \| * * * * * * * * \| * * * * * * * * \| * * * * * * * *
	1232	! \|* * * * * * * * * \|* * * * * * * * * \|* * * * * * * * *
	1233	!hitl(0)----------------------------------------------------------------
	1234	! aicetl(0) aicetl(1) aicetl(2) aicetl(3)
	1235
	1236	! move up over layers incrementing volume
	1237	DO n = 1, m_index+1
	1238
	1239	area = sum(aicetl(:)) - & ! total area of sub-layer
	1240	(rhos/rau0) * sum(aicetl(n:jpl+1)) ! area of sub-layer occupied by snow
	1241
	1242	vol = (hitl(n) - hitl(n-1)) * area ! thickness of sub-layer times area
	1243
	1244	IF (vol >= rem_vol) THEN ! have reached the sub-layer with the depth within
	1245	hpond = rem_vol / area + hitl(n-1) + alfan(m_index) - alfan(1)
	1246
	1247	exit
	1248	ELSE ! still in sub-layer below the sub-layer with the depth
	1249	rem_vol = rem_vol - vol
	1250	END IF
	1251
	1252	END DO
	1253
	1254	END IF
	1255
	1256	END SUBROUTINE ice_thd_pnd_depth
	1257
	1258
	1259	SUBROUTINE ice_thd_pnd_perm(ticen, salin, perm)
	1260	!!-------------------------------------------------------------------
	1261	!! * ROUTINE ice_thd_pnd_perm *
	1262	!!
	1263	!! ** Purpose : Determine the liquid fraction of brine in the ice
	1264	!! and its permeability
	1265	!!-------------------------------------------------------------------
	1266
	1267	REAL (wp), DIMENSION(nlay_i), INTENT(IN) :: &
	1268	ticen, & ! internal ice temperature (K)
	1269	salin ! salinity (ppt) !js: ppt according to cice
	1270
	1271	REAL (wp), INTENT(OUT) :: &
	1272	perm ! permeability
	1273
	1274	REAL (wp) :: &
	1275	Sbr ! brine salinity
	1276
	1277	REAL (wp), DIMENSION(nlay_i) :: &
	1278	Tin, & ! ice temperature
	1279	phi ! liquid fraction
	1280
	1281	INTEGER :: k
	1282
	1283	!-----------------------------------------------------------------
	1284	! Compute ice temperatures from enthalpies using quadratic formula
	1285	!-----------------------------------------------------------------
	1286
	1287	DO k = 1,nlay_i
	1288	Tin(k) = ticen(k) - rt0 !js: from K to degC
	1289	END DO
	1290
	1291	!-----------------------------------------------------------------
	1292	! brine salinity and liquid fraction
	1293	!-----------------------------------------------------------------
	1294
	1295	DO k = 1, nlay_i
	1296
	1297	! Sbr = - Tin(k) / rTmlt ! Consistent expression with SI3 (linear liquidus)
	1298	! Best expression to date is that one
	1299	Sbr = - 18.7 * Tin(k) - 0.519 * Tin(k)*2 - 0.00535 Tin(k) **3
	1300	phi(k) = salin(k) / Sbr
	1301
	1302	END DO
	1303
	1304	!-----------------------------------------------------------------
	1305	! permeability
	1306	!-----------------------------------------------------------------
	1307
	1308	perm = 3.0e-08_wp * (minval(phi))**3 ! Golden et al. (2007)
	1309
	1310	END SUBROUTINE ice_thd_pnd_perm
	1311
	1312
	1313	!----------------------------------------------------------------------------------------------------------------------
	1314
[8637]	1315	SUBROUTINE ice_thd_pnd_init
	1316	!!-------------------------------------------------------------------
	1317	!! * ROUTINE ice_thd_pnd_init *
	1318	!!
	1319	!! ** Purpose : Physical constants and parameters linked to melt ponds
	1320	!! over sea ice
	1321	!!
	1322	!! ** Method : Read the namthd_pnd namelist and check the melt pond
	1323	!! parameter values called at the first timestep (nit000)
	1324	!!
	1325	!! ** input : Namelist namthd_pnd
	1326	!!-------------------------------------------------------------------
[9169]	1327	INTEGER :: ios, ioptio ! Local integer
	1328	!!
[14659]	1329	NAMELIST/namthd_pnd/ ln_pnd, ln_pnd_LEV , rn_apnd_min, rn_apnd_max, rn_pnd_flush, &
	1330	& ln_pnd_CST , rn_apnd, rn_hpnd, &
	1331	& ln_pnd_TOPO , &
[13284]	1332	& ln_pnd_lids, ln_pnd_alb
[8637]	1333	!!-------------------------------------------------------------------
[9169]	1334	!
[8637]	1335	REWIND( numnam_ice_ref ) ! Namelist namthd_pnd in reference namelist : Melt Ponds
	1336	READ ( numnam_ice_ref, namthd_pnd, IOSTAT = ios, ERR = 901)
[11536]	1337	901 IF( ios /= 0 ) CALL ctl_nam ( ios , 'namthd_pnd in reference namelist' )
[8637]	1338	REWIND( numnam_ice_cfg ) ! Namelist namthd_pnd in configuration namelist : Melt Ponds
	1339	READ ( numnam_ice_cfg, namthd_pnd, IOSTAT = ios, ERR = 902 )
[11536]	1340	902 IF( ios > 0 ) CALL ctl_nam ( ios , 'namthd_pnd in configuration namelist' )
[8637]	1341	IF(lwm) WRITE ( numoni, namthd_pnd )
[9169]	1342	!
[8637]	1343	IF(lwp) THEN ! control print
	1344	WRITE(numout,*)
	1345	WRITE(numout,*) 'ice_thd_pnd_init: ice parameters for melt ponds'
	1346	WRITE(numout,*) '~~~~~~~~~~~~~~~~'
	1347	WRITE(numout,*) ' Namelist namicethd_pnd:'
[13284]	1348	WRITE(numout,*) ' Melt ponds activated or not ln_pnd = ', ln_pnd
[14659]	1349	WRITE(numout,*) ' Topographic melt pond scheme ln_pnd_TOPO = ', ln_pnd_TOPO
[13284]	1350	WRITE(numout,*) ' Level ice melt pond scheme ln_pnd_LEV = ', ln_pnd_LEV
	1351	WRITE(numout,*) ' Minimum ice fraction that contributes to melt ponds rn_apnd_min = ', rn_apnd_min
	1352	WRITE(numout,*) ' Maximum ice fraction that contributes to melt ponds rn_apnd_max = ', rn_apnd_max
[14659]	1353	WRITE(numout,*) ' Pond flushing efficiency rn_pnd_flush = ', rn_pnd_flush
[13284]	1354	WRITE(numout,*) ' Constant ice melt pond scheme ln_pnd_CST = ', ln_pnd_CST
	1355	WRITE(numout,*) ' Prescribed pond fraction rn_apnd = ', rn_apnd
	1356	WRITE(numout,*) ' Prescribed pond depth rn_hpnd = ', rn_hpnd
	1357	WRITE(numout,*) ' Frozen lids on top of melt ponds ln_pnd_lids = ', ln_pnd_lids
	1358	WRITE(numout,*) ' Melt ponds affect albedo or not ln_pnd_alb = ', ln_pnd_alb
[8637]	1359	ENDIF
	1360	!
	1361	! !== set the choice of ice pond scheme ==!
	1362	ioptio = 0
[11536]	1363	IF( .NOT.ln_pnd ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndNO ; ENDIF
	1364	IF( ln_pnd_CST ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndCST ; ENDIF
[13284]	1365	IF( ln_pnd_LEV ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndLEV ; ENDIF
[14659]	1366	IF( ln_pnd_TOPO ) THEN ; ioptio = ioptio + 1 ; nice_pnd = np_pndTOPO ; ENDIF
[11536]	1367	IF( ioptio /= 1 ) &
[14659]	1368	& 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]	1369	!
[8637]	1370	SELECT CASE( nice_pnd )
	1371	CASE( np_pndNO )
[13284]	1372	IF( ln_pnd_alb ) THEN ; ln_pnd_alb = .FALSE. ; CALL ctl_warn( 'ln_pnd_alb=false when no ponds' ) ; ENDIF
	1373	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when no ponds' ) ; ENDIF
	1374	CASE( np_pndCST )
	1375	IF( ln_pnd_lids ) THEN ; ln_pnd_lids = .FALSE. ; CALL ctl_warn( 'ln_pnd_lids=false when constant ponds' ) ; ENDIF
[8637]	1376	END SELECT
	1377	!
	1378	END SUBROUTINE ice_thd_pnd_init
	1379
	1380	#else
	1381	!!----------------------------------------------------------------------
[9570]	1382	!! Default option Empty module NO SI3 sea-ice model
[8637]	1383	!!----------------------------------------------------------------------
	1384	#endif
	1385
	1386	!!======================================================================
	1387	END MODULE icethd_pnd

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_fix_topo_ponds/src/ICE/icethd_pnd.F90 @ 15696

Download in other formats: