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.

limthd_dh.F90 in branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/2013/dev_r4028_CNRS_LIM3/NEMOGCM/NEMO/LIM_SRC_3/limthd_dh.F90 @ 4634

Last change on this file since 4634 was 4634, checked in by clem, 10 years ago
major changes in heat budget
Property svn:keywords set to `Id`
File size: 37.0 KB

Line
1	MODULE limthd_dh
2	!!======================================================================
3	!! * MODULE limthd_dh *
4	!! LIM-3 : thermodynamic growth and decay of the ice
5	!!======================================================================
6	!! History : LIM ! 2003-05 (M. Vancoppenolle) Original code in 1D
7	!! ! 2005-06 (M. Vancoppenolle) 3D version
8	!! 3.2 ! 2009-07 (M. Vancoppenolle, Y. Aksenov, G. Madec) bug correction in wfx_snw & wfx_ice
9	!! 3.4 ! 2011-02 (G. Madec) dynamical allocation
10	!! 3.5 ! 2012-10 (G. Madec & co) salt flux + bug fixes
11	!!----------------------------------------------------------------------
12	#if defined key_lim3
13	!!----------------------------------------------------------------------
14	!! 'key_lim3' LIM3 sea-ice model
15	!!----------------------------------------------------------------------
16	!! lim_thd_dh : vertical accr./abl. and lateral ablation of sea ice
17	!!----------------------------------------------------------------------
18	USE par_oce ! ocean parameters
19	USE phycst ! physical constants (OCE directory)
20	USE sbc_oce ! Surface boundary condition: ocean fields
21	USE ice ! LIM variables
22	USE par_ice ! LIM parameters
23	USE thd_ice ! LIM thermodynamics
24	USE in_out_manager ! I/O manager
25	USE lib_mpp ! MPP library
26	USE wrk_nemo ! work arrays
27	USE lib_fortran ! Fortran utilities (allows no signed zero when 'key_nosignedzero' defined)
28	USE cpl_oasis3, ONLY : lk_cpl
29
30	IMPLICIT NONE
31	PRIVATE
32
33	PUBLIC lim_thd_dh ! called by lim_thd
34
35	REAL(wp) :: epsi20 = 1.e-20 ! constant values
36	REAL(wp) :: epsi10 = 1.e-10 !
37
38	!!----------------------------------------------------------------------
39	!! NEMO/LIM3 4.0 , UCL - NEMO Consortium (2010)
40	!! $Id$
41	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
42	!!----------------------------------------------------------------------
43	CONTAINS
44
45	SUBROUTINE lim_thd_dh( kideb, kiut, jl )
46	!!------------------------------------------------------------------
47	!! * ROUTINE lim_thd_dh *
48	!!
49	!! ** Purpose : determines variations of ice and snow thicknesses.
50	!!
51	!! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes
52	!! Bottom accretion/ablation arises from flux budget
53	!! Snow thickness can increase by precipitation and decrease by sublimation
54	!! If snow load excesses Archmiede limit, snow-ice is formed by
55	!! the flooding of sea-water in the snow
56	!!
57	!! 1) Compute available flux of heat for surface ablation
58	!! 2) Compute snow and sea ice enthalpies
59	!! 3) Surface ablation and sublimation
60	!! 4) Bottom accretion/ablation
61	!! 5) Case of Total ablation
62	!! 6) Snow ice formation
63	!!
64	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
65	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
66	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
67	!! Vancoppenolle et al.,2009, Ocean Modelling
68	!!------------------------------------------------------------------
69	INTEGER , INTENT(in) :: kideb, kiut ! Start/End point on which the the computation is applied
70	INTEGER , INTENT(in) :: jl ! Thickness cateogry number
71	!!
72	INTEGER :: ji , jk ! dummy loop indices
73	INTEGER :: ii, ij ! 2D corresponding indices to ji
74	INTEGER :: i_ice_switch ! ice thickness above a certain treshold or not
75	INTEGER :: iter
76
77	REAL(wp) :: ztmelts ! local scalar
78	REAL(wp) :: zdh, zfdum !
79	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
80	REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads
81	REAL(wp) :: zs_snic ! snow-ice salinity
82	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
83	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
84	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
85	REAL(wp) :: zgrr ! bottom growth rate
86	REAL(wp) :: zt_i_new ! bottom formation temperature
87
88	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
89	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
90	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
91	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
92	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
93	REAL(wp) :: zsstK ! SST in Kelvin
94
95	REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness
96	REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3)
97	REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2)
98	REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2)
99	REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)
100	REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
101	REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
102	INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting
103
104	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt
105	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation
106	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation
107
108	REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah
109	REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness
110
111	REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2)
112	REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2)
113	REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3)
114
115	! mass and salt flux (clem)
116	REAL(wp) :: zdvres, zswitch_sal
117
118	! Heat conservation
119	INTEGER :: num_iter_max
120	REAL(wp) :: zinda, zindq, zindh
121	REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug
122
123	!!------------------------------------------------------------------
124
125	! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
126	SELECT CASE( num_sal ) ! varying salinity or not
127	CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile
128	CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile
129	END SELECT
130
131	CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
132	CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
133	CALL wrk_alloc( jpij, zintermelt )
134	CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )
135	CALL wrk_alloc( jpij, icount )
136
137	dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp
138	dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp
139
140	zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp
141	zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp
142
143	zh_s (:) = 0._wp
144	zdh_s_pre(:) = 0._wp
145	zdh_s_mel(:) = 0._wp
146	zdh_s_sub(:) = 0._wp
147	zqh_s (:) = 0._wp
148	zqh_i (:) = 0._wp
149
150	zh_i (:,:) = 0._wp
151	zdeltah (:,:) = 0._wp
152	zintermelt(:) = 0._wp
153	icount (:) = 0
154
155	! debug
156	dq_i(:) = 0._wp
157	dq_s(:) = 0._wp
158
159	! initialize layer thicknesses and enthalpies
160	h_i_old (:,0:nlay_i+1) = 0._wp
161	qh_i_old(:,0:nlay_i+1) = 0._wp
162	DO jk = 1, nlay_i
163	DO ji = kideb, kiut
164	h_i_old (ji,jk) = ht_i_b(ji) / REAL( nlay_i )
165	qh_i_old(ji,jk) = q_i_b(ji,jk) * h_i_old(ji,jk)
166	ENDDO
167	ENDDO
168	!
169	!------------------------------------------------------------------------------!
170	! 1) Calculate available heat for surface and bottom ablation !
171	!------------------------------------------------------------------------------!
172	!
173	DO ji = kideb, kiut
174	zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
175	ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt
176
177	zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
178	zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)
179
180	zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
181	zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
182	END DO ! ji
183
184	!
185	!------------------------------------------------------------------------------!
186	! If snow temperature is above freezing point, then snow melts
187	! (should not happen but sometimes it does)
188	!------------------------------------------------------------------------------!
189	DO ji = kideb, kiut
190	IF( t_s_b(ji,1) > rtt ) THEN !!! Internal melting
191	! Contribution to heat flux to the ocean [W.m-2], < 0
192	hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
193	! Contribution to mass flux
194	wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
195	! updates
196	ht_s_b(ji) = 0._wp
197	q_s_b (ji,1) = 0._wp
198	t_s_b (ji,1) = rtt
199	END IF
200	END DO
201
202	!------------------------------------------------------------!
203	! 2) Computing layer thicknesses and enthalpies. !
204	!------------------------------------------------------------!
205	!
206	DO ji = kideb, kiut
207	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
208	END DO
209	!
210	DO jk = 1, nlay_s
211	DO ji = kideb, kiut
212	zqh_s(ji) = zqh_s(ji) + q_s_b(ji,jk) * zh_s(ji)
213	END DO
214	END DO
215	!
216	DO jk = 1, nlay_i
217	DO ji = kideb, kiut
218	zh_i(ji,jk) = ht_i_b(ji) / REAL( nlay_i )
219	zqh_i(ji) = zqh_i(ji) + q_i_b(ji,jk) * zh_i(ji,jk)
220	END DO
221	END DO
222	!
223	!------------------------------------------------------------------------------\|
224	! 3) Surface ablation and sublimation \|
225	!------------------------------------------------------------------------------\|
226	!
227	!-------------------------
228	! 3.1 Snow precips / melt
229	!-------------------------
230	! Snow accumulation in one thermodynamic time step
231	! snowfall is partitionned between leads and ice
232	! if snow fall was uniform, a fraction (1-at_i) would fall into leads
233	! but because of the winds, more snow falls on leads than on sea ice
234	! and a greater fraction (1-at_i)^beta of the total mass of snow
235	! (beta < 1) falls in leads.
236	! In reality, beta depends on wind speed,
237	! and should decrease with increasing wind speed but here, it is
238	! considered as a constant. an average value is 0.66
239	! Martin Vancoppenolle, December 2006
240
241	DO ji = kideb, kiut
242	!-----------
243	! Snow fall
244	!-----------
245	! thickness change
246	zcoeff = ( 1._wp - ( 1._wp - at_i_b(ji) )**betas ) / at_i_b(ji)
247	zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
248	! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
249	zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )
250	IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
251	! heat flux from snow precip (>0, W.m-2)
252	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
253	! update thickness
254	ht_s_b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )
255
256	!---------------------
257	! Melt of falling snow
258	!---------------------
259	! thickness change
260	zindq = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
261	zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
262	zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting
263	! Heat flux associated with snow melt of falling snow (W.m-2, <0)
264	hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
265	! heat used to melt snow (W.m-2, >0)
266	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdh_s_mel(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
267	! snow melting only = water into the ocean (then without snow precip)
268	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_mel(ji) * r1_rdtice
269
270	! updates available heat + thickness
271	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdh_s_mel(ji) * zqprec(ji) )
272	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_mel(ji) )
273	zh_s (ji) = ht_s_b(ji) / REAL( nlay_s )
274
275	! clem debug: variation of enthalpy (J.m-2)
276	dq_s(ji) = dq_s(ji) + ( zdh_s_pre(ji) + zdh_s_mel(ji) ) * zqprec(ji)
277
278	END DO
279
280	! If heat still available, then melt more snow
281	zdeltah(:,:) = 0._wp ! important
282	DO jk = 1, nlay_s
283	DO ji = kideb, kiut
284	! thickness change
285	zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) )
286	zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - q_s_b(ji,jk) + epsi20 ) )
287	zdeltah (ji,jk) = - zindh * zindq * zq_su(ji) / MAX( q_s_b(ji,jk), epsi20 )
288	zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji) ) ! bound melting
289	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk)
290	! heat flux associated with snow melt(W.m-2, <0)
291	hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice
292	! heat used to melt snow(W.m-2, >0)
293	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdeltah(ji,jk) * a_i_b(ji) * q_s_b(ji,jk) * r1_rdtice
294	! snow melting only = water into the ocean (then without snow precip)
295	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
296
297	! updates available heat + thickness
298	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_b(ji,jk) )
299	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )
300
301	! clem debug: variation of enthalpy (J.m-2)
302	dq_s(ji) = dq_s(ji) + zdeltah(ji,jk) * q_s_b(ji,jk)
303	END DO
304	END DO
305
306	!----------------------
307	! 3.2 Snow sublimation
308	!----------------------
309	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
310	IF( lk_cpl ) THEN
311	! coupled mode: sublimation already included in emp_ice (to do in limsbc_ice)
312	zdh_s_sub(:) = 0._wp
313	ELSE
314	! forced mode: snow thickness change due to sublimation
315	DO ji = kideb, kiut
316	zdh_s_sub(ji) = MAX( - ht_s_b(ji) , - parsub * qla_ice_1d(ji) / ( rhosn * lsub ) * rdt_ice )
317	! Heat flux by sublimation [W.m-2], < 0
318	! sublimate first snow that had fallen, then pre-existing snow
319	zcoeff = ( MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) * zqprec(ji) + &
320	& ( zdh_s_sub(ji) - MAX( zdh_s_sub(ji), - MAX( 0._wp, zdh_s_pre(ji) + zdh_s_mel(ji) ) ) ) * q_s_b(ji,1) ) &
321	& * a_i_b(ji) * r1_rdtice
322	hfx_sub_1d(ji) = hfx_sub_1d(ji) + zcoeff ! diag only (to close heat budget)
323	! heat used for sublimation (>0, W.m-2)
324	!!? hfx_tot_1d(ji) = hfx_tot_1d(ji) - zcoeff
325	! Mass flux by sublimation
326	wfx_sub_1d(ji) = wfx_sub_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice ! diag only
327	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
328	! new snow thickness
329	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
330	! clem debug: variation of enthalpy (J.m-2)
331	dq_s(ji) = dq_s(ji) + zdh_s_sub(ji) * q_s_b(ji,1)
332	END DO
333	ENDIF
334
335	! --- Update snow diags --- !
336	DO ji = kideb, kiut
337	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
338	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
339	END DO ! ji
340
341	!-------------------------------------------
342	! 3.3 Update temperature, energy
343	!-------------------------------------------
344	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
345	zq_s(:) = 0._wp
346	DO jk = 1, nlay_s
347	DO ji = kideb,kiut
348	zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) + epsi20 ) )
349	q_s_b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) * &
350	& ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + &
351	& ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )
352	zq_s(ji) = zq_s(ji) + q_s_b(ji,jk)
353	END DO
354	END DO
355
356	!--------------------------
357	! 3.4 Surface ice ablation
358	!--------------------------
359	zdeltah(:,:) = 0._wp ! important
360	DO jk = 1, nlay_i
361	DO ji = kideb, kiut
362	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]
363
364	ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer k [K]
365
366	zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0]
367
368	zdE = zEi - zEw ! Specific enthalpy difference < 0
369
370	zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
371
372	zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0]
373
374	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk) , - zh_i(ji,jk) ) ) ! Melt of layer jk cannot exceed the layer thickness [m, <0]
375
376	zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
377
378	dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt
379
380	zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
381
382	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
383
384	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
385	sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
386
387	! Contribution to heat flux [W.m-2], < 0
388	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
389
390	! Total heat flux used in this process [W.m-2], < 0
391	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
392
393	! Contribution to mass flux
394	wfx_sum_1d(ji) = wfx_sum_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
395
396	! record which layers have disappeared (for bottom melting)
397	! => icount=0 : no layer has vanished
398	! => icount=5 : 5 layers have vanished
399	zindh = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) )
400	icount(ji) = icount(ji) + zindh
401	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
402
403	! clem debug: variation of enthalpy (J.m-2)
404	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
405
406	! update heat content (J.m-2) and layer thickness
407	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
408	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
409	END DO
410	END DO
411	! update ice thickness
412	DO ji = kideb, kiut
413	ht_i_b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )
414	END DO
415
416	!
417	!------------------------------------------------------------------------------!
418	! 4) Basal growth / melt !
419	!------------------------------------------------------------------------------!
420	!
421	!------------------
422	! 4.1 Basal growth
423	!------------------
424	! Basal growth is driven by heat imbalance at the ice-ocean interface,
425	! between the inner conductive flux (fc_bo_i), from the open water heat flux
426	! (fhldb) and the turbulent ocean flux (fhtur).
427	! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
428
429	! If salinity varies in time, an iterative procedure is required, because
430	! the involved quantities are inter-dependent.
431	! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
432	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
433	! -> need for an iterative procedure, which converges quickly
434
435	IF ( num_sal == 2 ) THEN
436	num_iter_max = 5
437	ELSE
438	num_iter_max = 1
439	ENDIF
440
441	!clem debug. Just to be sure that enthalpy at nlay_i+1 is null
442	DO ji = kideb, kiut
443	q_i_b(ji,nlay_i+1) = 0._wp
444	END DO
445
446	! Iterative procedure
447	DO iter = 1, num_iter_max
448	DO ji = kideb, kiut
449	IF( zf_tt(ji) < 0._wp ) THEN
450
451	! New bottom ice salinity (Cox & Weeks, JGR88 )
452	!--- zswi1 if dh/dt < 2.0e-8
453	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
454	!--- zswi2 if dh/dt > 3.6e-7
455	zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
456	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
457	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
458	zswi1 = 1. - zswi2 * zswi12
459	zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
460	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
461
462	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
463
464	s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity
465	+ ( 1. - zswitch_sal ) * sm_i_b(ji)
466	! New ice growth
467	ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K)
468
469	zt_i_new = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
470
471	zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0)
472	& - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) &
473	& + rcp * ( ztmelts-rtt )
474
475	zEw = rcp * ( t_bo_b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
476
477	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
478
479	dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
480
481	q_i_b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)
482
483	ENDIF ! fc_bo_i
484	END DO ! ji
485	END DO ! iter
486
487	! Contribution to Energy and Salt Fluxes
488	DO ji = kideb, kiut
489	IF( zf_tt(ji) < 0._wp ) THEN
490	! New ice growth
491
492	zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0)
493
494	! Contribution to heat flux to the ocean [W.m-2], >0
495	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
496	! Total heat flux used in this process [W.m-2]
497	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
498
499	! Contribution to salt flux ()
500	sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
501
502	! Contribution to mass flux
503	wfx_bog_1d(ji) = wfx_bog_1d(ji) + rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
504
505	! clem debug: variation of enthalpy (J.m-2)
506	dq_i(ji) = dq_i(ji) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
507
508	! update heat content (J.m-2) and layer thickness
509	qh_i_old(ji,nlay_i+1) = qh_i_old(ji,nlay_i+1) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
510	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
511	ENDIF
512	END DO
513
514	!----------------
515	! 4.2 Basal melt
516	!----------------
517	zdeltah(:,:) = 0._wp ! important
518	DO jk = nlay_i, 1, -1
519	DO ji = kideb, kiut
520	IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting
521
522	ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer jk (K)
523
524	IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
525	zintermelt(ji) = 1._wp
526
527	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
528
529	!!zEw = rcp * ( t_i_b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b (J/kg, <0)
530
531	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
532	! set up at 0 since no energy is needed to melt water...(it is already melted)
533
534	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
535	! this should normally not happen, but sometimes, heat diffusion leads to this
536
537	dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk)
538
539	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
540
541	! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
542	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice
543
544	! clem debug: variation of enthalpy (J.m-2)
545	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
546
547	! update heat content (J.m-2) and layer thickness
548	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
549	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
550
551	ELSE !!! Basal melting
552
553	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
554
555	zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
556
557	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
558
559	zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)
560
561	zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change
562
563	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
564
565	zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
566
567	dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt
568
569	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
570
571	zQm = zfmdt * zEw ! Heat exchanged with ocean
572
573	! Contribution to heat flux to the ocean [W.m-2], <0
574	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
575
576	! clem debug: variation of enthalpy (J.m-2)
577	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
578
579	! update heat content (J.m-2) and layer thickness
580	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
581	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
582	ENDIF
583
584	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
585	sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
586
587	! Total heat flux used in this process [W.m-2]
588	hfx_tot_1d(ji) = hfx_tot_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
589
590	! Contribution to mass flux
591	wfx_bom_1d(ji) = wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
592
593	ENDIF
594	END DO ! ji
595	END DO ! jk
596
597	!------------------------------------------------------------------------------!
598	! Excessive ablation in a 1-category model
599	! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
600	!------------------------------------------------------------------------------!
601	! ??? keep ???
602	! clem bug: I think this should be included above, so we would not have to
603	! track heat/salt/mass fluxes backwards
604	IF( jpl == 1 ) THEN
605	DO ji = kideb, kiut
606	IF( zf_tt(ji) >= 0._wp ) THEN
607	zdh = MAX( hmelt , dh_i_bott(ji) )
608	zdvres = zdh - dh_i_bott(ji) ! >=0
609	dh_i_bott(ji) = zdh
610
611	! excessive energy is sent to lateral ablation
612	zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
613	zq_1cat(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
614
615	! correct salt and mass fluxes
616	sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdvres * rhoic * r1_rdtice ! this is only a raw approximation
617	wfx_bom_1d(ji) = wfx_bom_1d(ji) + rhoic * a_i_b(ji) * zdvres * r1_rdtice
618	ENDIF
619	END DO
620	ENDIF
621
622	!-------------------------------------------
623	! Update temperature, energy
624	!-------------------------------------------
625	DO ji = kideb, kiut
626	ht_i_b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )
627	END DO
628
629	!-------------------------------------------
630	! 5. What to do with remaining energy
631	!-------------------------------------------
632	! If heat still available for melting and snow remains, then melt more snow
633	!-------------------------------------------
634	zdeltah(:,:) = 0._wp ! important
635	DO ji = kideb, kiut
636	zq_rema(ji) = zq_su(ji) + zq_bo(ji)
637	! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) ) ! =1 if snow
638	! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) )
639	! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
640	! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_b(ji) ) ) ! bound melting
641	! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)
642	! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1)
643	! ht_s_b (ji) = ht_s_b(ji) + zdeltah(ji,1)
644	!
645	! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2)
646	! ! Heat flux associated with snow melt
647	! hfx_snw_1d(ji) = hfx_snw_1d(ji) + zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (<0)
648	! ! heat used to melt snow
649	! hfx_tot_1d(ji) = hfx_tot_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
650	! ! Contribution to mass flux
651	! wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
652	! ! clem debug: variation of enthalpy (J.m-2)
653	! dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)
654	!
655	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
656	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
657	hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_b(ji) ) * r1_rdtice
658
659	IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
660	END DO
661
662	!
663	!------------------------------------------------------------------------------\|
664	! 6) Snow-Ice formation \|
665	!------------------------------------------------------------------------------\|
666	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
667	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
668	DO ji = kideb, kiut
669	!
670	dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )
671
672	ht_i_b(ji) = ht_i_b(ji) + dh_snowice(ji)
673	ht_s_b(ji) = ht_s_b(ji) - dh_snowice(ji)
674
675	! Salinity of snow ice
676	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
677	zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_b(ji)
678
679	! entrapment during snow ice formation
680	! new salinity difference stored (to be used in limthd_ent.F90)
681	IF ( num_sal == 2 ) THEN
682	i_ice_switch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
683	! salinity dif due to snow-ice formation
684	dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch
685	! salinity dif due to bottom growth
686	IF ( zf_tt(ji) < 0._wp ) THEN
687	dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * i_ice_switch
688	ENDIF
689	ENDIF
690
691	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
692	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
693	zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp ) ! <0
694	zsstK = sst_m(ii,ij) + rt0
695	zEw = rcp * ( zsstK - rt0 )
696	zQm = zfmdt * zEw
697
698	! Contribution to heat flux
699	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
700
701	! Contribution to salt flux
702	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice
703
704	! Contribution to mass flux
705	! All snow is thrown in the ocean, and seawater is taken to replace the volume
706	wfx_sni_1d(ji) = wfx_sni_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
707	wfx_snw_1d(ji) = wfx_snw_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
708
709	! clem debug: variation of enthalpy (J.m-2)
710	dq_s(ji) = dq_s(ji) - dh_snowice(ji) * q_s_b(ji,1)
711	dq_i(ji) = dq_i(ji) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
712
713	! update heat content (J.m-2) and layer thickness
714	qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
715	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
716
717	! Total ablation (to debug)
718	IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp
719
720	END DO !ji
721
722	!
723	!-------------------------------------------
724	! Update temperature, energy
725	!-------------------------------------------
726	!clem bug: we should take snow into account here
727	DO ji = kideb, kiut
728	zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )
729	t_su_b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt
730	END DO ! ji
731
732	DO jk = 1, nlay_s
733	DO ji = kideb,kiut
734	! mask enthalpy
735	zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) ) )
736	q_s_b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)
737	! recalculate t_s_b from q_s_b
738	t_s_b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
739	END DO
740	END DO
741
742	CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
743	CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
744	CALL wrk_dealloc( jpij, zintermelt )
745	CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )
746	CALL wrk_dealloc( jpij, icount )
747	!
748	!
749	END SUBROUTINE lim_thd_dh
750
751	#else
752	!!----------------------------------------------------------------------
753	!! Default option NO LIM3 sea-ice model
754	!!----------------------------------------------------------------------
755	CONTAINS
756	SUBROUTINE lim_thd_dh ! Empty routine
757	END SUBROUTINE lim_thd_dh
758	#endif
759
760	!!======================================================================
761	END MODULE limthd_dh

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

Download in other formats: