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

Last change on this file since 4649 was 4649, checked in by clem, 10 years ago
finalizing LIM3 heat budget conservation + multiple minor bugs corrections
Property svn:keywords set to `Id`
File size: 37.7 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 :: iter
75
76	REAL(wp) :: ztmelts ! local scalar
77	REAL(wp) :: zdh, zfdum !
78	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
79	REAL(wp) :: zcoeff ! dummy argument for snowfall partitioning over ice and leads
80	REAL(wp) :: zs_snic ! snow-ice salinity
81	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
82	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
83	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
84	REAL(wp) :: zgrr ! bottom growth rate
85	REAL(wp) :: zt_i_new ! bottom formation temperature
86
87	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
88	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
89	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
90	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
91	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
92	REAL(wp) :: zsstK ! SST in Kelvin
93
94	REAL(wp), POINTER, DIMENSION(:) :: zh_s ! snow layer thickness
95	REAL(wp), POINTER, DIMENSION(:) :: zqprec ! energy of fallen snow (J.m-3)
96	REAL(wp), POINTER, DIMENSION(:) :: zq_su ! heat for surface ablation (J.m-2)
97	REAL(wp), POINTER, DIMENSION(:) :: zq_bo ! heat for bottom ablation (J.m-2)
98	REAL(wp), POINTER, DIMENSION(:) :: zq_1cat ! corrected heat in case 1-cat and hmelt>15cm (J.m-2)
99	REAL(wp), POINTER, DIMENSION(:) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
100	REAL(wp), POINTER, DIMENSION(:) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
101	INTEGER , POINTER, DIMENSION(:) :: icount ! number of layers vanished by melting
102
103	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_mel ! snow melt
104	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_pre ! snow precipitation
105	REAL(wp), POINTER, DIMENSION(:) :: zdh_s_sub ! snow sublimation
106
107	REAL(wp), POINTER, DIMENSION(:,:) :: zdeltah
108	REAL(wp), POINTER, DIMENSION(:,:) :: zh_i ! ice layer thickness
109
110	REAL(wp), POINTER, DIMENSION(:) :: zqh_i ! total ice heat content (J.m-2)
111	REAL(wp), POINTER, DIMENSION(:) :: zqh_s ! total snow heat content (J.m-2)
112	REAL(wp), POINTER, DIMENSION(:) :: zq_s ! total snow enthalpy (J.m-3)
113
114	! mass and salt flux (clem)
115	REAL(wp) :: zdvres, zswitch_sal, zswitch
116
117	! Heat conservation
118	INTEGER :: num_iter_max
119	REAL(wp) :: zinda, zindq, zindh
120	REAL(wp), POINTER, DIMENSION(:) :: zintermelt ! debug
121
122	!!------------------------------------------------------------------
123
124	! Discriminate between varying salinity (num_sal=2) and prescribed cases (other values)
125	SELECT CASE( num_sal ) ! varying salinity or not
126	CASE( 1, 3, 4 ) ; zswitch_sal = 0 ! prescribed salinity profile
127	CASE( 2 ) ; zswitch_sal = 1 ! varying salinity profile
128	END SELECT
129
130	CALL wrk_alloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
131	CALL wrk_alloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
132	CALL wrk_alloc( jpij, zintermelt )
133	CALL wrk_alloc( jpij, jkmax, zdeltah, zh_i )
134	CALL wrk_alloc( jpij, icount )
135
136	dh_i_surf (:) = 0._wp ; dh_i_bott (:) = 0._wp ; dh_snowice(:) = 0._wp
137	dsm_i_se_1d(:) = 0._wp ; dsm_i_si_1d(:) = 0._wp
138
139	zqprec (:) = 0._wp ; zq_su (:) = 0._wp ; zq_bo (:) = 0._wp ; zf_tt (:) = 0._wp
140	zq_1cat(:) = 0._wp ; zq_rema(:) = 0._wp
141
142	zh_s (:) = 0._wp
143	zdh_s_pre(:) = 0._wp
144	zdh_s_mel(:) = 0._wp
145	zdh_s_sub(:) = 0._wp
146	zqh_s (:) = 0._wp
147	zqh_i (:) = 0._wp
148
149	zh_i (:,:) = 0._wp
150	zdeltah (:,:) = 0._wp
151	zintermelt(:) = 0._wp
152	icount (:) = 0
153
154	! debug
155	dq_i(:) = 0._wp
156	dq_s(:) = 0._wp
157
158	! initialize layer thicknesses and enthalpies
159	h_i_old (:,0:nlay_i+1) = 0._wp
160	qh_i_old(:,0:nlay_i+1) = 0._wp
161	DO jk = 1, nlay_i
162	DO ji = kideb, kiut
163	h_i_old (ji,jk) = ht_i_b(ji) / REAL( nlay_i )
164	qh_i_old(ji,jk) = q_i_b(ji,jk) * h_i_old(ji,jk)
165	ENDDO
166	ENDDO
167	!
168	!------------------------------------------------------------------------------!
169	! 1) Calculate available heat for surface and bottom ablation !
170	!------------------------------------------------------------------------------!
171	!
172	DO ji = kideb, kiut
173	zinda = 1._wp - MAX( 0._wp , SIGN( 1._wp , - ht_s_b(ji) ) )
174	ztmelts = zinda * rtt + ( 1._wp - zinda ) * rtt
175
176	zfdum = qns_ice_1d(ji) + ( 1._wp - i0(ji) ) * qsr_ice_1d(ji) - fc_su(ji)
177	zf_tt(ji) = fc_bo_i(ji) + fhtur_1d(ji) + fhld_1d(ji)
178
179	zq_su (ji) = MAX( 0._wp, zfdum * rdt_ice ) * MAX( 0._wp , SIGN( 1._wp, t_su_b(ji) - ztmelts ) )
180	zq_bo (ji) = MAX( 0._wp, zf_tt(ji) * rdt_ice )
181	END DO
182
183	!
184	!------------------------------------------------------------------------------!
185	! If snow temperature is above freezing point, then snow melts
186	! (should not happen but sometimes it does)
187	!------------------------------------------------------------------------------!
188	DO ji = kideb, kiut
189	IF( t_s_b(ji,1) > rtt ) THEN !!! Internal melting
190	! Contribution to heat flux to the ocean [W.m-2], < 0
191	hfx_res_1d(ji) = hfx_res_1d(ji) + q_s_b(ji,1) * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
192	! Contribution to mass flux
193	wfx_snw_1d(ji) = wfx_snw_1d(ji) + rhosn * ht_s_b(ji) * a_i_b(ji) * r1_rdtice
194	! updates
195	ht_s_b(ji) = 0._wp
196	q_s_b (ji,1) = 0._wp
197	t_s_b (ji,1) = rtt
198	END IF
199	END DO
200
201	!------------------------------------------------------------!
202	! 2) Computing layer thicknesses and enthalpies. !
203	!------------------------------------------------------------!
204	!
205	DO ji = kideb, kiut
206	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
207	END DO
208	!
209	DO jk = 1, nlay_s
210	DO ji = kideb, kiut
211	zqh_s(ji) = zqh_s(ji) + q_s_b(ji,jk) * zh_s(ji)
212	END DO
213	END DO
214	!
215	DO jk = 1, nlay_i
216	DO ji = kideb, kiut
217	zh_i(ji,jk) = ht_i_b(ji) / REAL( nlay_i )
218	zqh_i(ji) = zqh_i(ji) + q_i_b(ji,jk) * zh_i(ji,jk)
219	END DO
220	END DO
221	!
222	!------------------------------------------------------------------------------\|
223	! 3) Surface ablation and sublimation \|
224	!------------------------------------------------------------------------------\|
225	!
226	!-------------------------
227	! 3.1 Snow precips / melt
228	!-------------------------
229	! Snow accumulation in one thermodynamic time step
230	! snowfall is partitionned between leads and ice
231	! if snow fall was uniform, a fraction (1-at_i) would fall into leads
232	! but because of the winds, more snow falls on leads than on sea ice
233	! and a greater fraction (1-at_i)^beta of the total mass of snow
234	! (beta < 1) falls in leads.
235	! In reality, beta depends on wind speed,
236	! and should decrease with increasing wind speed but here, it is
237	! considered as a constant. an average value is 0.66
238	! Martin Vancoppenolle, December 2006
239
240	DO ji = kideb, kiut
241	!-----------
242	! Snow fall
243	!-----------
244	! thickness change
245	zcoeff = ( 1._wp - ( 1._wp - at_i_b(ji) )**betas ) / at_i_b(ji)
246	zdh_s_pre(ji) = zcoeff * sprecip_1d(ji) * rdt_ice / rhosn
247	! enthalpy of the precip (>0, J.m-3) (tatm_ice is now in K)
248	zqprec (ji) = rhosn * ( cpic * ( rtt - MIN( tatm_ice_1d(ji), rt0_snow) ) + lfus )
249	IF( sprecip_1d(ji) == 0._wp ) zqprec(ji) = 0._wp
250	! heat flux from snow precip (>0, W.m-2)
251	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_b(ji) * zqprec(ji) * r1_rdtice
252	! mass flux, <0
253	wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhosn * a_i_b(ji) * zdh_s_pre(ji) * r1_rdtice
254	! update thickness
255	ht_s_b (ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_pre(ji) )
256
257	!---------------------
258	! Melt of falling snow
259	!---------------------
260	! thickness change
261	IF( zdh_s_pre(ji) > 0._wp ) THEN
262	zindq = 1._wp - MAX( 0._wp , SIGN( 1._wp , - zqprec(ji) + epsi20 ) )
263	zdh_s_mel (ji) = - zindq * zq_su(ji) / MAX( zqprec(ji) , epsi20 )
264	zdh_s_mel (ji) = MAX( - zdh_s_pre(ji), zdh_s_mel(ji) ) ! bound melting
265	! heat used to melt snow (W.m-2, >0)
266	hfx_snw_1d(ji) = hfx_snw_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), >0
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	ENDIF
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 used to melt snow(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	! snow melting only = water into the ocean (then without snow precip)
293	wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
294
295	! updates available heat + thickness
296	zq_su (ji) = MAX( 0._wp , zq_su (ji) + zdeltah(ji,jk) * q_s_b(ji,jk) )
297	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdeltah(ji,jk) )
298
299	! clem debug: variation of enthalpy (J.m-2)
300	dq_s(ji) = dq_s(ji) + zdeltah(ji,jk) * q_s_b(ji,jk)
301	END DO
302	END DO
303
304	!----------------------
305	! 3.2 Snow sublimation
306	!----------------------
307	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
308	! clem comment: not counted in mass exchange in limsbc since this is an exchange with atm. (not ocean)
309	! clem comment: ice should also sublimate
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
323	! Mass flux by sublimation
324	wfx_sub_1d(ji) = wfx_sub_1d(ji) - rhosn * a_i_b(ji) * zdh_s_sub(ji) * r1_rdtice
325	! new snow thickness
326	ht_s_b(ji) = MAX( 0._wp , ht_s_b(ji) + zdh_s_sub(ji) )
327	! clem debug: variation of enthalpy (J.m-2)
328	dq_s(ji) = dq_s(ji) + zdh_s_sub(ji) * q_s_b(ji,1)
329	END DO
330	ENDIF
331
332	! --- Update snow diags --- !
333	DO ji = kideb, kiut
334	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
335	zh_s(ji) = ht_s_b(ji) / REAL( nlay_s )
336	END DO ! ji
337
338	!-------------------------------------------
339	! 3.3 Update temperature, energy
340	!-------------------------------------------
341	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
342	zq_s(:) = 0._wp
343	DO jk = 1, nlay_s
344	DO ji = kideb,kiut
345	zindh = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) + epsi20 ) )
346	q_s_b(ji,jk) = ( 1._wp - zindh ) / MAX( ht_s_b(ji), epsi20 ) * &
347	& ( ( MAX( 0._wp, dh_s_tot(ji) ) ) * zqprec(ji) + &
348	& ( - MAX( 0._wp, dh_s_tot(ji) ) + ht_s_b(ji) ) * rhosn * ( cpic * ( rtt - t_s_b(ji,jk) ) + lfus ) )
349	zq_s(ji) = zq_s(ji) + q_s_b(ji,jk)
350	END DO
351	END DO
352
353	!--------------------------
354	! 3.4 Surface ice ablation
355	!--------------------------
356	zdeltah(:,:) = 0._wp ! important
357	DO jk = 1, nlay_i
358	DO ji = kideb, kiut
359	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of layer k [J/kg, <0]
360
361	ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer k [K]
362
363	zEw = rcp * ( ztmelts - rt0 ) ! Specific enthalpy of resulting meltwater [J/kg, <0]
364
365	zdE = zEi - zEw ! Specific enthalpy difference < 0
366
367	zfmdt = - zq_su(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
368
369	zdeltah(ji,jk) = - zfmdt / rhoic ! Melt of layer jk [m, <0]
370
371	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]
372
373	zq_su(ji) = MAX( 0._wp , zq_su(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat
374
375	dh_i_surf(ji) = dh_i_surf(ji) + zdeltah(ji,jk) ! Cumulate surface melt
376
377	zfmdt = - rhoic * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
378
379	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
380
381	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
382	sfx_sum_1d(ji) = sfx_sum_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
383
384	! Contribution to heat flux [W.m-2], < 0
385	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
386
387	! Total heat flux used in this process [W.m-2], > 0
388	hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
389
390	! Contribution to mass flux
391	wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
392
393	! record which layers have disappeared (for bottom melting)
394	! => icount=0 : no layer has vanished
395	! => icount=5 : 5 layers have vanished
396	zindh = NINT( MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) ) )
397	icount(ji) = icount(ji) + zindh
398	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
399
400	! clem debug: variation of enthalpy (J.m-2)
401	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
402
403	! update heat content (J.m-2) and layer thickness
404	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
405	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
406	END DO
407	END DO
408	! update ice thickness
409	DO ji = kideb, kiut
410	ht_i_b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_surf(ji) )
411	END DO
412
413	!
414	!------------------------------------------------------------------------------!
415	! 4) Basal growth / melt !
416	!------------------------------------------------------------------------------!
417	!
418	!------------------
419	! 4.1 Basal growth
420	!------------------
421	! Basal growth is driven by heat imbalance at the ice-ocean interface,
422	! between the inner conductive flux (fc_bo_i), from the open water heat flux
423	! (fhld) and the turbulent ocean flux (fhtur).
424	! fc_bo_i is positive downwards. fhtur and fhld are positive to the ice
425
426	! If salinity varies in time, an iterative procedure is required, because
427	! the involved quantities are inter-dependent.
428	! Basal growth (dh_i_bott) depends upon new ice specific enthalpy (zEi),
429	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bott)
430	! -> need for an iterative procedure, which converges quickly
431
432	IF ( num_sal == 2 ) THEN
433	num_iter_max = 5
434	ELSE
435	num_iter_max = 1
436	ENDIF
437
438	!clem debug. Just to be sure that enthalpy at nlay_i+1 is null
439	DO ji = kideb, kiut
440	q_i_b(ji,nlay_i+1) = 0._wp
441	END DO
442
443	! Iterative procedure
444	DO iter = 1, num_iter_max
445	DO ji = kideb, kiut
446	IF( zf_tt(ji) < 0._wp ) THEN
447
448	! New bottom ice salinity (Cox & Weeks, JGR88 )
449	!--- zswi1 if dh/dt < 2.0e-8
450	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
451	!--- zswi2 if dh/dt > 3.6e-7
452	zgrr = MIN( 1.0e-3, MAX ( dh_i_bott(ji) * r1_rdtice , epsi10 ) )
453	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
454	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
455	zswi1 = 1. - zswi2 * zswi12
456	zfracs = MIN ( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
457	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
458
459	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
460
461	s_i_new(ji) = zswitch_sal * zfracs * sss_m(ii,ij) & ! New ice salinity
462	+ ( 1. - zswitch_sal ) * sm_i_b(ji)
463	! New ice growth
464	ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K)
465
466	zt_i_new = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
467
468	zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0)
469	& - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) &
470	& + rcp * ( ztmelts-rtt )
471
472	zEw = rcp * ( t_bo_b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
473
474	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
475
476	dh_i_bott(ji) = rdt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoic ) )
477
478	q_i_b(ji,nlay_i+1) = -zEi * rhoic ! New ice energy of melting (J/m3, >0)
479
480	ENDIF ! fc_bo_i
481	END DO ! ji
482	END DO ! iter
483
484	! Contribution to Energy and Salt Fluxes
485	DO ji = kideb, kiut
486	IF( zf_tt(ji) < 0._wp ) THEN
487	! New ice growth
488
489	zfmdt = - rhoic * dh_i_bott(ji) ! Mass flux x time step (kg/m2, < 0)
490
491	ztmelts = - tmut * s_i_new(ji) + rtt ! New ice melting point (K)
492
493	zt_i_new = zswitch_sal * t_bo_b(ji) + ( 1. - zswitch_sal) * t_i_b(ji, nlay_i)
494
495	zEi = cpic * ( zt_i_new - ztmelts ) & ! Specific enthalpy of forming ice (J/kg, <0)
496	& - lfus * ( 1.0 - ( ztmelts - rtt ) / ( zt_i_new - rtt ) ) &
497	& + rcp * ( ztmelts-rtt )
498
499	zEw = rcp * ( t_bo_b(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
500
501	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
502
503	! Contribution to heat flux to the ocean [W.m-2], >0
504	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
505
506	! Total heat flux used in this process [W.m-2], <0
507	hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
508
509	! Contribution to salt flux, <0
510	sfx_bog_1d(ji) = sfx_bog_1d(ji) + s_i_new(ji) * a_i_b(ji) * zfmdt * r1_rdtice
511
512	! Contribution to mass flux, <0
513	wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoic * a_i_b(ji) * dh_i_bott(ji) * r1_rdtice
514
515	! clem debug: variation of enthalpy (J.m-2)
516	dq_i(ji) = dq_i(ji) + dh_i_bott(ji) * q_i_b(ji,nlay_i+1)
517
518	! update heat content (J.m-2) and layer thickness
519	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)
520	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bott(ji)
521	ENDIF
522	END DO
523
524	!----------------
525	! 4.2 Basal melt
526	!----------------
527	zdeltah(:,:) = 0._wp ! important
528	DO jk = nlay_i, 1, -1
529	DO ji = kideb, kiut
530	IF( zf_tt(ji) >= 0._wp .AND. jk > icount(ji) ) THEN ! do not calculate where layer has already disappeared from surface melting
531
532	ztmelts = - tmut * s_i_b(ji,jk) + rtt ! Melting point of layer jk (K)
533
534	IF( t_i_b(ji,jk) >= ztmelts ) THEN !!! Internal melting
535	zintermelt(ji) = 1._wp
536
537	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
538
539	!!zEw = rcp * ( t_i_b(ji,jk) - rtt ) ! Specific enthalpy of meltwater at T = t_i_b (J/kg, <0)
540
541	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
542	! set up at 0 since no energy is needed to melt water...(it is already melted)
543
544	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
545	! this should normally not happen, but sometimes, heat diffusion leads to this
546
547	dh_i_bott (ji) = dh_i_bott(ji) + zdeltah(ji,jk)
548
549	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
550
551	! Contribution to heat flux to the ocean [W.m-2], <0 (ice enthalpy zEi is "sent" to the ocean)
552	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_b(ji) * zEi * r1_rdtice
553
554	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
555	sfx_res_1d(ji) = sfx_res_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
556
557	! Contribution to mass flux
558	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
559
560	! clem debug: variation of enthalpy (J.m-2)
561	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
562
563	! update heat content (J.m-2) and layer thickness
564	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
565	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
566
567	ELSE !!! Basal melting
568
569	zEi = - q_i_b(ji,jk) / rhoic ! Specific enthalpy of melting ice (J/kg, <0)
570
571	zEw = rcp * ( ztmelts - rtt )! Specific enthalpy of meltwater (J/kg, <0)
572
573	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
574
575	zfmdt = - zq_bo(ji) / zdE ! Mass flux x time step (kg/m2, >0)
576
577	zdeltah(ji,jk) = - zfmdt / rhoic ! Gross thickness change
578
579	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
580
581	zq_bo(ji) = MAX( 0._wp , zq_bo(ji) - zdeltah(ji,jk) * rhoic * zdE ) ! update available heat. MAX is necessary for roundup errors
582
583	dh_i_bott(ji) = dh_i_bott(ji) + zdeltah(ji,jk) ! Update basal melt
584
585	zfmdt = - zdeltah(ji,jk) * rhoic ! Mass flux x time step > 0
586
587	zQm = zfmdt * zEw ! Heat exchanged with ocean
588
589	! Contribution to heat flux to the ocean [W.m-2], <0
590	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
591
592	! Contribution to salt flux (clem: using sm_i_b and not s_i_b(jk) is ok)
593	sfx_bom_1d(ji) = sfx_bom_1d(ji) - sm_i_b(ji) * a_i_b(ji) * zdeltah(ji,jk) * rhoic * r1_rdtice
594
595	! Total heat flux used in this process [W.m-2], >0
596	hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_b(ji) * zdE * r1_rdtice
597
598	! Contribution to mass flux
599	wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdeltah(ji,jk) * r1_rdtice
600
601	! clem debug: variation of enthalpy (J.m-2)
602	dq_i(ji) = dq_i(ji) + zdeltah(ji,jk) * q_i_b(ji,jk)
603
604	! update heat content (J.m-2) and layer thickness
605	qh_i_old(ji,jk) = qh_i_old(ji,jk) + zdeltah(ji,jk) * q_i_b(ji,jk)
606	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
607	ENDIF
608
609	ENDIF
610	END DO ! ji
611	END DO ! jk
612
613	!------------------------------------------------------------------------------!
614	! Excessive ablation in a 1-category model
615	! in a 1-category sea ice model, bottom ablation must not exceed hmelt (-0.15)
616	!------------------------------------------------------------------------------!
617	! ??? keep ???
618	! clem bug: I think this should be included above, so we would not have to
619	! track heat/salt/mass fluxes backwards
620	! IF( jpl == 1 ) THEN
621	! DO ji = kideb, kiut
622	! IF( zf_tt(ji) >= 0._wp ) THEN
623	! zdh = MAX( hmelt , dh_i_bott(ji) )
624	! zdvres = zdh - dh_i_bott(ji) ! >=0
625	! dh_i_bott(ji) = zdh
626	!
627	! ! excessive energy is sent to lateral ablation
628	! zinda = MAX( 0._wp, SIGN( 1._wp , 1._wp - at_i_b(ji) - epsi20 ) )
629	! zq_1cat(ji) = zinda * rhoic * lfus * at_i_b(ji) / MAX( 1._wp - at_i_b(ji) , epsi20 ) * zdvres ! J.m-2 >=0
630	!
631	! ! correct salt and mass fluxes
632	! 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
633	! wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoic * a_i_b(ji) * zdvres * r1_rdtice
634	! ENDIF
635	! END DO
636	! ENDIF
637
638	!-------------------------------------------
639	! Update temperature, energy
640	!-------------------------------------------
641	DO ji = kideb, kiut
642	ht_i_b(ji) = MAX( 0._wp , ht_i_b(ji) + dh_i_bott(ji) )
643	END DO
644
645	!-------------------------------------------
646	! 5. What to do with remaining energy
647	!-------------------------------------------
648	! If heat still available for melting and snow remains, then melt more snow
649	!-------------------------------------------
650	zdeltah(:,:) = 0._wp ! important
651	DO ji = kideb, kiut
652	zq_rema(ji) = zq_su(ji) + zq_bo(ji)
653	! zindh = 1._wp - MAX( 0._wp, SIGN( 1._wp, - ht_s_b(ji) ) ) ! =1 if snow
654	! zindq = 1._wp - MAX( 0._wp, SIGN( 1._wp, - zq_s(ji) + epsi20 ) )
655	! zdeltah (ji,1) = - zindh * zindq * zq_rema(ji) / MAX( zq_s(ji), epsi20 )
656	! zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - ht_s_b(ji) ) ) ! bound melting
657	! zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,1)
658	! dh_s_tot (ji) = dh_s_tot(ji) + zdeltah(ji,1)
659	! ht_s_b (ji) = ht_s_b(ji) + zdeltah(ji,1)
660	!
661	! zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * zq_s(ji) ! update available heat (J.m-2)
662	! ! heat used to melt snow
663	! hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_b(ji) * zq_s(ji) * r1_rdtice ! W.m-2 (>0)
664	! ! Contribution to mass flux
665	! wfx_snw_1d(ji) = wfx_snw_1d(ji) - rhosn * a_i_b(ji) * zdeltah(ji,1) * r1_rdtice
666	! ! clem debug: variation of enthalpy (J.m-2)
667	! dq_s(ji) = dq_s(ji) + zdeltah(ji,1) * q_s_b(ji,1)
668	!
669	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
670	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
671	hfx_out(ii,ij) = hfx_out(ii,ij) + ( zq_1cat(ji) + zq_rema(ji) * a_i_b(ji) ) * r1_rdtice
672
673	IF( ln_nicep .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
674	END DO
675
676	!
677	!------------------------------------------------------------------------------\|
678	! 6) Snow-Ice formation \|
679	!------------------------------------------------------------------------------\|
680	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
681	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
682	DO ji = kideb, kiut
683	!
684	dh_snowice(ji) = MAX( 0._wp , ( rhosn * ht_s_b(ji) + (rhoic-rau0) * ht_i_b(ji) ) / ( rhosn+rau0-rhoic ) )
685
686	ht_i_b(ji) = ht_i_b(ji) + dh_snowice(ji)
687	ht_s_b(ji) = ht_s_b(ji) - dh_snowice(ji)
688
689	! Salinity of snow ice
690	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
691	zs_snic = zswitch_sal * sss_m(ii,ij) * ( rhoic - rhosn ) / rhoic + ( 1. - zswitch_sal ) * sm_i_b(ji)
692
693	! entrapment during snow ice formation
694	! new salinity difference stored (to be used in limthd_ent.F90)
695	IF ( num_sal == 2 ) THEN
696	zswitch = MAX( 0._wp , SIGN( 1._wp , ht_i_b(ji) - epsi10 ) )
697	! salinity dif due to snow-ice formation
698	dsm_i_si_1d(ji) = ( zs_snic - sm_i_b(ji) ) * dh_snowice(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch
699	! salinity dif due to bottom growth
700	IF ( zf_tt(ji) < 0._wp ) THEN
701	dsm_i_se_1d(ji) = ( s_i_new(ji) - sm_i_b(ji) ) * dh_i_bott(ji) / MAX( ht_i_b(ji), epsi10 ) * zswitch
702	ENDIF
703	ENDIF
704
705	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
706	ii = MOD( npb(ji) - 1, jpi ) + 1 ; ij = ( npb(ji) - 1 ) / jpi + 1
707	zfmdt = ( rhosn - rhoic ) * MAX( dh_snowice(ji), 0._wp ) ! <0
708	zsstK = sst_m(ii,ij) + rt0
709	zEw = rcp * ( zsstK - rt0 )
710	zQm = zfmdt * zEw
711
712	! Contribution to heat flux
713	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_b(ji) * zEw * r1_rdtice
714
715	! Contribution to salt flux
716	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_m(ii,ij) * a_i_b(ji) * zfmdt * r1_rdtice
717
718	! Contribution to mass flux
719	! All snow is thrown in the ocean, and seawater is taken to replace the volume
720	wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_b(ji) * dh_snowice(ji) * rhoic * r1_rdtice
721	wfx_snw_1d(ji) = wfx_snw_1d(ji) + a_i_b(ji) * dh_snowice(ji) * rhosn * r1_rdtice
722
723	! clem debug: variation of enthalpy (J.m-2)
724	dq_s(ji) = dq_s(ji) - dh_snowice(ji) * q_s_b(ji,1)
725	dq_i(ji) = dq_i(ji) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
726
727	! update heat content (J.m-2) and layer thickness
728	qh_i_old(ji,0) = qh_i_old(ji,0) + dh_snowice(ji) * q_s_b(ji,1) + zfmdt * zEw
729	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
730
731	! Total ablation (to debug)
732	IF( ht_i_b(ji) <= 0._wp ) a_i_b(ji) = 0._wp
733
734	END DO !ji
735
736	!
737	!-------------------------------------------
738	! Update temperature, energy
739	!-------------------------------------------
740	!clem bug: we should take snow into account here
741	DO ji = kideb, kiut
742	zindh = 1.0 - MAX( 0._wp , SIGN( 1._wp , - ht_i_b(ji) ) )
743	t_su_b(ji) = zindh * t_su_b(ji) + ( 1.0 - zindh ) * rtt
744	END DO ! ji
745
746	DO jk = 1, nlay_s
747	DO ji = kideb,kiut
748	! mask enthalpy
749	zinda = MAX( 0._wp , SIGN( 1._wp, - ht_s_b(ji) ) )
750	q_s_b(ji,jk) = ( 1.0 - zinda ) * q_s_b(ji,jk)
751	! recalculate t_s_b from q_s_b
752	t_s_b(ji,jk) = rtt + ( 1._wp - zinda ) * ( - q_s_b(ji,jk) / ( rhosn * cpic ) + lfus / cpic )
753	END DO
754	END DO
755
756	CALL wrk_dealloc( jpij, zh_s, zqprec, zq_su, zq_bo, zf_tt, zq_1cat, zq_rema )
757	CALL wrk_dealloc( jpij, zdh_s_mel, zdh_s_pre, zdh_s_sub, zqh_i, zqh_s, zq_s )
758	CALL wrk_dealloc( jpij, zintermelt )
759	CALL wrk_dealloc( jpij, jkmax, zdeltah, zh_i )
760	CALL wrk_dealloc( jpij, icount )
761	!
762	!
763	END SUBROUTINE lim_thd_dh
764
765	#else
766	!!----------------------------------------------------------------------
767	!! Default option NO LIM3 sea-ice model
768	!!----------------------------------------------------------------------
769	CONTAINS
770	SUBROUTINE lim_thd_dh ! Empty routine
771	END SUBROUTINE lim_thd_dh
772	#endif
773
774	!!======================================================================
775	END MODULE limthd_dh

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

Download in other formats: