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

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

icethd_dh.F90 in NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r12512_HPC-04_mcastril_Mixed_Precision_implementation/src/ICE/icethd_dh.F90 @ 13181

Last change on this file since 13181 was 13181, checked in by orioltp, 4 years ago
Several bugfixes added by Sam Hatfield (ECMWF).
Property svn:keywords set to `Id`
File size: 35.8 KB

Line
1	MODULE icethd_dh
2	!!======================================================================
3	!! * MODULE icethd_dh *
4	!! seaice : thermodynamic growth and melt
5	!!======================================================================
6	!! History : ! 2003-05 (M. Vancoppenolle) Original code in 1D
7	!! ! 2005-06 (M. Vancoppenolle) 3D version
8	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
9	!!----------------------------------------------------------------------
10	#if defined key_si3
11	!!----------------------------------------------------------------------
12	!! 'key_si3' SI3 sea-ice model
13	!!----------------------------------------------------------------------
14	!! ice_thd_dh : vertical sea-ice growth and melt
15	!! ice_thd_snwblow : distribute snow fall between ice and ocean
16	!!----------------------------------------------------------------------
17	USE dom_oce ! ocean space and time domain
18	USE phycst ! physical constants
19	USE ice ! sea-ice: variables
20	USE ice1D ! sea-ice: thermodynamics variables
21	USE icethd_sal ! sea-ice: salinity profiles
22	!
23	USE in_out_manager ! I/O manager
24	USE lib_mpp ! MPP library
25	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
26
27	IMPLICIT NONE
28	PRIVATE
29
30	PUBLIC ice_thd_dh ! called by ice_thd
31	PUBLIC ice_thd_snwblow ! called in sbcblk/sbccpl and here
32
33	INTERFACE ice_thd_snwblow
34	MODULE PROCEDURE ice_thd_snwblow_1d, ice_thd_snwblow_2d
35	END INTERFACE
36
37	!!----------------------------------------------------------------------
38	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
39	!! $Id$
40	!! Software governed by the CeCILL license (see ./LICENSE)
41	!!----------------------------------------------------------------------
42	CONTAINS
43
44	SUBROUTINE ice_thd_dh
45	!!------------------------------------------------------------------
46	!! * ROUTINE ice_thd_dh *
47	!!
48	!! ** Purpose : compute ice and snow thickness changes due to growth/melting
49	!!
50	!! ** Method : Ice/Snow surface melting arises from imbalance in surface fluxes
51	!! Bottom accretion/ablation arises from flux budget
52	!! Snow thickness can increase by precipitation and decrease by sublimation
53	!! If snow load excesses Archmiede limit, snow-ice is formed by
54	!! the flooding of sea-water in the snow
55	!!
56	!! - Compute available flux of heat for surface ablation
57	!! - Compute snow and sea ice enthalpies
58	!! - Surface ablation and sublimation
59	!! - Bottom accretion/ablation
60	!! - Snow ice formation
61	!!
62	!! References : Bitz and Lipscomb, 1999, J. Geophys. Res.
63	!! Fichefet T. and M. Maqueda 1997, J. Geophys. Res., 102(C6), 12609-12646
64	!! Vancoppenolle, Fichefet and Bitz, 2005, Geophys. Res. Let.
65	!! Vancoppenolle et al.,2009, Ocean Modelling
66	!!------------------------------------------------------------------
67	INTEGER :: ji, jk ! dummy loop indices
68	INTEGER :: iter ! local integer
69
70	REAL(wp) :: ztmelts ! local scalar
71	REAL(wp) :: zdum
72	REAL(wp) :: zfracs ! fractionation coefficient for bottom salt entrapment
73	REAL(wp) :: zswi1 ! switch for computation of bottom salinity
74	REAL(wp) :: zswi12 ! switch for computation of bottom salinity
75	REAL(wp) :: zswi2 ! switch for computation of bottom salinity
76	REAL(wp) :: zgrr ! bottom growth rate
77	REAL(wp) :: zt_i_new ! bottom formation temperature
78	REAL(wp) :: z1_rho ! 1/(rhos+rho0-rhoi)
79
80	REAL(wp) :: zQm ! enthalpy exchanged with the ocean (J/m2), >0 towards the ocean
81	REAL(wp) :: zEi ! specific enthalpy of sea ice (J/kg)
82	REAL(wp) :: zEw ! specific enthalpy of exchanged water (J/kg)
83	REAL(wp) :: zdE ! specific enthalpy difference (J/kg)
84	REAL(wp) :: zfmdt ! exchange mass flux x time step (J/m2), >0 towards the ocean
85
86	REAL(wp), DIMENSION(jpij) :: zqprec ! energy of fallen snow (J.m-3)
87	REAL(wp), DIMENSION(jpij) :: zq_top ! heat for surface ablation (J.m-2)
88	REAL(wp), DIMENSION(jpij) :: zq_bot ! heat for bottom ablation (J.m-2)
89	REAL(wp), DIMENSION(jpij) :: zq_rema ! remaining heat at the end of the routine (J.m-2)
90	REAL(wp), DIMENSION(jpij) :: zf_tt ! Heat budget to determine melting or freezing(W.m-2)
91	REAL(wp), DIMENSION(jpij) :: zevap_rema ! remaining mass flux from sublimation (kg.m-2)
92
93	REAL(wp), DIMENSION(jpij) :: zdh_s_mel ! snow melt
94	REAL(wp), DIMENSION(jpij) :: zdh_s_pre ! snow precipitation
95	REAL(wp), DIMENSION(jpij) :: zdh_s_sub ! snow sublimation
96
97	REAL(wp), DIMENSION(jpij,nlay_s) :: zh_s ! snw layer thickness
98	REAL(wp), DIMENSION(jpij,nlay_i) :: zh_i ! ice layer thickness
99	REAL(wp), DIMENSION(jpij,nlay_i) :: zdeltah
100	INTEGER , DIMENSION(jpij,nlay_i) :: icount ! number of layers vanished by melting
101
102	REAL(wp), DIMENSION(jpij) :: zsnw ! distribution of snow after wind blowing
103
104	REAL(wp) :: zswitch_sal
105
106	INTEGER :: num_iter_max ! Heat conservation
107	!!------------------------------------------------------------------
108
109	! Discriminate between time varying salinity and constant
110	SELECT CASE( nn_icesal ) ! varying salinity or not
111	CASE( 1 , 3 ) ; zswitch_sal = 0._wp ! prescribed salinity profile
112	CASE( 2 ) ; zswitch_sal = 1._wp ! varying salinity profile
113	END SELECT
114
115	! initialize layer thicknesses and enthalpies
116	h_i_old (1:npti,0:nlay_i+1) = 0._wp
117	eh_i_old(1:npti,0:nlay_i+1) = 0._wp
118	DO jk = 1, nlay_i
119	DO ji = 1, npti
120	h_i_old (ji,jk) = h_i_1d(ji) * r1_nlay_i
121	eh_i_old(ji,jk) = e_i_1d(ji,jk) * h_i_old(ji,jk)
122	END DO
123	END DO
124	!
125	! ! ============================================== !
126	! ! Available heat for surface and bottom ablation !
127	! ! ============================================== !
128	!
129	IF( ln_cndflx .AND. .NOT.ln_cndemulate ) THEN
130	!
131	DO ji = 1, npti
132	zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
133	END DO
134	!
135	ELSE
136	!
137	DO ji = 1, npti
138	zdum = qns_ice_1d(ji) + qsr_ice_1d(ji) - qtr_ice_top_1d(ji) - qcn_ice_top_1d(ji)
139	qml_ice_1d(ji) = zdum * MAX( 0._wp , SIGN( 1._wp, t_su_1d(ji) - rt0 ) )
140	zq_top(ji) = MAX( 0._wp, qml_ice_1d(ji) * rDt_ice )
141	END DO
142	!
143	ENDIF
144	!
145	DO ji = 1, npti
146	zf_tt(ji) = qcn_ice_bot_1d(ji) + qsb_ice_bot_1d(ji) + fhld_1d(ji)
147	zq_bot(ji) = MAX( 0._wp, zf_tt(ji) * rDt_ice )
148	END DO
149
150	! Ice and snow layer thicknesses
151	!-------------------------------
152	DO jk = 1, nlay_i
153	DO ji = 1, npti
154	zh_i(ji,jk) = h_i_1d(ji) * r1_nlay_i
155	END DO
156	END DO
157
158	DO jk = 1, nlay_s
159	DO ji = 1, npti
160	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
161	END DO
162	END DO
163
164	! ! ============ !
165	! ! Snow !
166	! ! ============ !
167	!
168	! Internal melting
169	! ----------------
170	! IF snow temperature is above freezing point, THEN snow melts (should not happen but sometimes it does)
171	DO jk = 1, nlay_s
172	DO ji = 1, npti
173	IF( t_s_1d(ji,jk) > rt0 ) THEN
174	hfx_res_1d (ji) = hfx_res_1d (ji) + e_s_1d(ji,jk) * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! heat flux to the ocean [W.m-2], < 0
175	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) + rhos * zh_s(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! mass flux
176	! updates
177	dh_s_mlt(ji) = dh_s_mlt(ji) - zh_s(ji,jk)
178	h_s_1d (ji) = h_s_1d(ji) - zh_s(ji,jk)
179	zh_s (ji,jk) = 0._wp
180	e_s_1d (ji,jk) = 0._wp
181	t_s_1d (ji,jk) = rt0
182	END IF
183	END DO
184	END DO
185
186	! Snow precipitation
187	!-------------------
188	CALL ice_thd_snwblow( 1.0_wp - at_i_1d(1:npti), zsnw(1:npti) ) ! snow distribution over ice after wind blowing
189
190	zdeltah(1:npti,:) = 0._wp
191	DO ji = 1, npti
192	IF( sprecip_1d(ji) > 0._wp ) THEN
193	!
194	! --- precipitation ---
195	zdh_s_pre (ji) = zsnw(ji) * sprecip_1d(ji) * rDt_ice * r1_rhos / at_i_1d(ji) ! thickness change
196	zqprec (ji) = - qprec_ice_1d(ji) ! enthalpy of the precip (>0, J.m-3)
197	!
198	hfx_spr_1d(ji) = hfx_spr_1d(ji) + zdh_s_pre(ji) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat flux from snow precip (>0, W.m-2)
199	wfx_spr_1d(ji) = wfx_spr_1d(ji) - rhos * a_i_1d(ji) * zdh_s_pre(ji) * r1_Dt_ice ! mass flux, <0
200
201	! --- melt of falling snow ---
202	rswitch = MAX( 0._wp , SIGN( 1._wp , zqprec(ji) - epsi20 ) )
203	zdeltah (ji,1) = - rswitch * zq_top(ji) / MAX( zqprec(ji) , epsi20 ) ! thickness change
204	zdeltah (ji,1) = MAX( - zdh_s_pre(ji), zdeltah(ji,1) ) ! bound melting
205	hfx_snw_1d (ji) = hfx_snw_1d (ji) - zdeltah(ji,1) * a_i_1d(ji) * zqprec(ji) * r1_Dt_ice ! heat used to melt snow (W.m-2, >0)
206	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip), >0
207
208	! updates available heat + precipitations after melting
209	dh_s_mlt (ji) = dh_s_mlt(ji) + zdeltah(ji,1)
210	zq_top (ji) = MAX( 0._wp , zq_top (ji) + zdeltah(ji,1) * zqprec(ji) )
211	zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
212
213	! update thickness
214	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_pre(ji) )
215	!
216	ELSE
217	!
218	zdh_s_pre(ji) = 0._wp
219	zqprec (ji) = 0._wp
220	!
221	ENDIF
222	END DO
223
224	! recalculate snow layers
225	DO jk = 1, nlay_s
226	DO ji = 1, npti
227	zh_s(ji,jk) = h_s_1d(ji) * r1_nlay_s
228	END DO
229	END DO
230
231	! Snow melting
232	! ------------
233	! If heat still available (zq_top > 0), then melt more snow
234	zdeltah(1:npti,:) = 0._wp
235	zdh_s_mel(1:npti) = 0._wp
236	DO jk = 1, nlay_s
237	DO ji = 1, npti
238	IF( zh_s(ji,jk) > 0._wp .AND. zq_top(ji) > 0._wp ) THEN
239	!
240	rswitch = MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,jk) - epsi20 ) )
241	zdeltah (ji,jk) = - rswitch * zq_top(ji) / MAX( e_s_1d(ji,jk), epsi20 ) ! thickness change
242	zdeltah (ji,jk) = MAX( zdeltah(ji,jk) , - zh_s(ji,jk) ) ! bound melting
243	zdh_s_mel(ji) = zdh_s_mel(ji) + zdeltah(ji,jk)
244
245	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,jk) * a_i_1d(ji) * e_s_1d (ji,jk) * r1_Dt_ice ! heat used to melt snow(W.m-2, >0)
246	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! snow melting only = water into the ocean (then without snow precip)
247
248	! updates available heat + thickness
249	dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,jk)
250	zq_top (ji) = MAX( 0._wp , zq_top(ji) + zdeltah(ji,jk) * e_s_1d(ji,jk) )
251	h_s_1d (ji) = MAX( 0._wp , h_s_1d(ji) + zdeltah(ji,jk) )
252	zh_s (ji,jk) = MAX( 0._wp , zh_s(ji,jk) + zdeltah(ji,jk) )
253	!
254	ENDIF
255	END DO
256	END DO
257
258	! Snow sublimation
259	!-----------------
260	! qla_ice is always >=0 (upwards), heat goes to the atmosphere, therefore snow sublimates
261	! comment: not counted in mass/heat exchange in iceupdate.F90 since this is an exchange with atm. (not ocean)
262	zdeltah(1:npti,:) = 0._wp
263	DO ji = 1, npti
264	IF( evap_ice_1d(ji) > 0._wp ) THEN
265	!
266	zdh_s_sub (ji) = MAX( - h_s_1d(ji) , - evap_ice_1d(ji) * r1_rhos * rDt_ice )
267	zevap_rema(ji) = evap_ice_1d(ji) * rDt_ice + zdh_s_sub(ji) * rhos ! remaining evap in kg.m-2 (used for ice melting later on)
268	zdeltah (ji,1) = MAX( zdh_s_sub(ji), - zdh_s_pre(ji) )
269
270	hfx_sub_1d (ji) = hfx_sub_1d(ji) + & ! Heat flux by sublimation [W.m-2], < 0 (sublimate snow that had fallen, then pre-existing snow)
271	& ( zdeltah(ji,1) * zqprec(ji) + ( zdh_s_sub(ji) - zdeltah(ji,1) ) * e_s_1d(ji,1) ) &
272	& * a_i_1d(ji) * r1_Dt_ice
273	wfx_snw_sub_1d(ji) = wfx_snw_sub_1d(ji) - rhos * a_i_1d(ji) * zdh_s_sub(ji) * r1_Dt_ice ! Mass flux by sublimation
274
275	! new snow thickness
276	h_s_1d(ji) = MAX( 0._wp , h_s_1d(ji) + zdh_s_sub(ji) )
277	! update precipitations after sublimation and correct sublimation
278	zdh_s_pre(ji) = zdh_s_pre(ji) + zdeltah(ji,1)
279	zdh_s_sub(ji) = zdh_s_sub(ji) - zdeltah(ji,1)
280	!
281	ELSE
282	!
283	zdh_s_sub (ji) = 0._wp
284	zevap_rema(ji) = 0._wp
285	!
286	ENDIF
287	END DO
288
289	! --- Update snow diags --- !
290	DO ji = 1, npti
291	dh_s_tot(ji) = zdh_s_mel(ji) + zdh_s_pre(ji) + zdh_s_sub(ji)
292	END DO
293
294	! Update temperature, energy
295	!---------------------------
296	! new temp and enthalpy of the snow (remaining snow precip + remaining pre-existing snow)
297	DO jk = 1, nlay_s
298	DO ji = 1,npti
299	rswitch = MAX( 0._wp , SIGN( 1._wp, h_s_1d(ji) - epsi20 ) )
300	e_s_1d(ji,jk) = rswitch / MAX( h_s_1d(ji), epsi20 ) * &
301	& ( ( zdh_s_pre(ji) ) * zqprec(ji) + &
302	& ( h_s_1d(ji) - zdh_s_pre(ji) ) * rhos * ( rcpi * ( rt0 - t_s_1d(ji,jk) ) + rLfus ) )
303	END DO
304	END DO
305
306	! ! ============ !
307	! ! Ice !
308	! ! ============ !
309
310	! Surface ice melting
311	!--------------------
312	zdeltah(1:npti,:) = 0._wp ! important
313	DO jk = 1, nlay_i
314	DO ji = 1, npti
315	ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer k [C]
316
317	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
318
319	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
320	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
321	! set up at 0 since no energy is needed to melt water...(it is already melted)
322	zdeltah(ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
323	! this should normally not happen, but sometimes, heat diffusion leads to this
324	zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0
325
326	dh_i_itm(ji) = dh_i_itm(ji) + zdeltah(ji,jk) ! Cumulate internal melting
327
328	zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
329
330	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0
331	! ice enthalpy zEi is "sent" to the ocean
332	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux
333	! using s_i_1d and not sz_i_1d(jk) is ok
334	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux
335
336	ELSE !-- Surface melting
337
338	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of layer k [J/kg, <0]
339	zEw = rcp * ztmelts ! Specific enthalpy of resulting meltwater [J/kg, <0]
340	zdE = zEi - zEw ! Specific enthalpy difference < 0
341
342	zfmdt = - zq_top(ji) / zdE ! Mass flux to the ocean [kg/m2, >0]
343
344	zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Melt of layer jk [m, <0]
345
346	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]
347
348	zq_top(ji) = MAX( 0._wp , zq_top(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat
349
350	dh_i_sum(ji) = dh_i_sum(ji) + zdeltah(ji,jk) ! Cumulate surface melt
351
352	zfmdt = - rhoi * zdeltah(ji,jk) ! Recompute mass flux [kg/m2, >0]
353
354	zQm = zfmdt * zEw ! Energy of the melt water sent to the ocean [J/m2, <0]
355
356	sfx_sum_1d(ji) = sfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0
357	! using s_i_1d and not sz_i_1d(jk) is ok)
358	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux [W.m-2], < 0
359	hfx_sum_1d(ji) = hfx_sum_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], > 0
360	!
361	wfx_sum_1d(ji) = wfx_sum_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux
362
363	END IF
364
365	! Ice sublimation
366	! ---------------
367	zdum = MAX( - ( zh_i(ji,jk) + zdeltah(ji,jk) ) , - zevap_rema(ji) * r1_rhoi )
368	zdeltah (ji,jk) = zdeltah (ji,jk) + zdum
369	dh_i_sub(ji) = dh_i_sub(ji) + zdum
370
371	sfx_sub_1d(ji) = sfx_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * s_i_1d(ji) * r1_Dt_ice ! Salt flux >0
372	! clem: flux is sent to the ocean for simplicity
373	! but salt should remain in the ice except
374	! if all ice is melted. => must be corrected
375	hfx_sub_1d(ji) = hfx_sub_1d(ji) + zdum * e_i_1d(ji,jk) * a_i_1d(ji) * r1_Dt_ice ! Heat flux [W.m-2], < 0
376
377	wfx_ice_sub_1d(ji) = wfx_ice_sub_1d(ji) - rhoi * a_i_1d(ji) * zdum * r1_Dt_ice ! Mass flux > 0
378
379	! update remaining mass flux
380	zevap_rema(ji) = zevap_rema(ji) + zdum * rhoi
381
382	! record which layers have disappeared (for bottom melting)
383	! => icount=0 : no layer has vanished
384	! => icount=5 : 5 layers have vanished
385	rswitch = MAX( 0._wp , SIGN( 1._wp , - ( zh_i(ji,jk) + zdeltah(ji,jk) ) ) )
386	icount(ji,jk) = NINT( rswitch )
387	zh_i(ji,jk) = MAX( 0._wp , zh_i(ji,jk) + zdeltah(ji,jk) )
388
389	! update heat content (J.m-2) and layer thickness
390	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
391	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
392	END DO
393	END DO
394
395	! update ice thickness
396	DO ji = 1, npti
397	h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_sum(ji) + dh_i_itm(ji) + dh_i_sub(ji) )
398	END DO
399
400	! remaining "potential" evap is sent to ocean
401	DO ji = 1, npti
402	wfx_err_sub_1d(ji) = wfx_err_sub_1d(ji) - zevap_rema(ji) * a_i_1d(ji) * r1_Dt_ice ! <=0 (net evap for the ocean in kg.m-2.s-1)
403	END DO
404
405
406	! Ice Basal growth
407	!------------------
408	! Basal growth is driven by heat imbalance at the ice-ocean interface,
409	! between the inner conductive flux (qcn_ice_bot), from the open water heat flux
410	! (fhld) and the sensible ice-ocean flux (qsb_ice_bot).
411	! qcn_ice_bot is positive downwards. qsb_ice_bot and fhld are positive to the ice
412
413	! If salinity varies in time, an iterative procedure is required, because
414	! the involved quantities are inter-dependent.
415	! Basal growth (dh_i_bog) depends upon new ice specific enthalpy (zEi),
416	! which depends on forming ice salinity (s_i_new), which depends on dh/dt (dh_i_bog)
417	! -> need for an iterative procedure, which converges quickly
418
419	num_iter_max = 1
420	IF( nn_icesal == 2 ) num_iter_max = 5 ! salinity varying in time
421
422	DO ji = 1, npti
423	IF( zf_tt(ji) < 0._wp ) THEN
424	DO iter = 1, num_iter_max ! iterations
425
426	! New bottom ice salinity (Cox & Weeks, JGR88 )
427	!--- zswi1 if dh/dt < 2.0e-8
428	!--- zswi12 if 2.0e-8 < dh/dt < 3.6e-7
429	!--- zswi2 if dh/dt > 3.6e-7
430	zgrr = MIN( 1.0e-3, MAX ( dh_i_bog(ji) * r1_Dt_ice , epsi10 ) )
431	zswi2 = MAX( 0._wp , SIGN( 1._wp , zgrr - 3.6e-7 ) )
432	zswi12 = MAX( 0._wp , SIGN( 1._wp , zgrr - 2.0e-8 ) ) * ( 1.0 - zswi2 )
433	zswi1 = 1. - zswi2 * zswi12
434	zfracs = MIN( zswi1 * 0.12 + zswi12 * ( 0.8925 + 0.0568 * LOG( 100.0 * zgrr ) ) &
435	& + zswi2 * 0.26 / ( 0.26 + 0.74 * EXP ( - 724300.0 * zgrr ) ) , 0.5 )
436
437	s_i_new(ji) = zswitch_sal * zfracs * sss_1d(ji) + ( 1. - zswitch_sal ) * s_i_1d(ji) ! New ice salinity
438
439	ztmelts = - rTmlt * s_i_new(ji) ! New ice melting point (C)
440
441	zt_i_new = zswitch_sal * t_bo_1d(ji) + ( 1. - zswitch_sal) * t_i_1d(ji, nlay_i)
442
443	zEi = rcpi * ( zt_i_new - (ztmelts+rt0) ) & ! Specific enthalpy of forming ice (J/kg, <0)
444	& - rLfus * ( 1.0 - ztmelts / ( MIN( zt_i_new - rt0, -epsi10 ) ) ) + rcp * ztmelts
445
446	zEw = rcp * ( t_bo_1d(ji) - rt0 ) ! Specific enthalpy of seawater (J/kg, < 0)
447
448	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
449
450	dh_i_bog(ji) = rDt_ice * MAX( 0._wp , zf_tt(ji) / ( zdE * rhoi ) )
451
452	END DO
453	! Contribution to Energy and Salt Fluxes
454	zfmdt = - rhoi * dh_i_bog(ji) ! Mass flux x time step (kg/m2, < 0)
455
456	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], >0
457	hfx_bog_1d(ji) = hfx_bog_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat flux used in this process [W.m-2], <0
458
459	sfx_bog_1d(ji) = sfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * s_i_new(ji) * r1_Dt_ice ! Salt flux, <0
460
461	wfx_bog_1d(ji) = wfx_bog_1d(ji) - rhoi * a_i_1d(ji) * dh_i_bog(ji) * r1_Dt_ice ! Mass flux, <0
462
463	! update heat content (J.m-2) and layer thickness
464	eh_i_old(ji,nlay_i+1) = eh_i_old(ji,nlay_i+1) + dh_i_bog(ji) * (-zEi * rhoi)
465	h_i_old (ji,nlay_i+1) = h_i_old (ji,nlay_i+1) + dh_i_bog(ji)
466
467	ENDIF
468
469	END DO
470
471	! Ice Basal melt
472	!---------------
473	zdeltah(1:npti,:) = 0._wp ! important
474	DO jk = nlay_i, 1, -1
475	DO ji = 1, npti
476	IF( zf_tt(ji) > 0._wp .AND. jk > icount(ji,jk) ) THEN ! do not calculate where layer has already disappeared by surface melting
477
478	ztmelts = - rTmlt * sz_i_1d(ji,jk) ! Melting point of layer jk (C)
479
480	IF( t_i_1d(ji,jk) >= (ztmelts+rt0) ) THEN !-- Internal melting
481
482	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
483	zdE = 0._wp ! Specific enthalpy difference (J/kg, <0)
484	! set up at 0 since no energy is needed to melt water...(it is already melted)
485	zdeltah (ji,jk) = MIN( 0._wp , - zh_i(ji,jk) ) ! internal melting occurs when the internal temperature is above freezing
486	! this should normally not happen, but sometimes, heat diffusion leads to this
487
488	dh_i_itm (ji) = dh_i_itm(ji) + zdeltah(ji,jk)
489
490	zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0
491
492	hfx_res_1d(ji) = hfx_res_1d(ji) + zfmdt * a_i_1d(ji) * zEi * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0
493	! ice enthalpy zEi is "sent" to the ocean
494	sfx_res_1d(ji) = sfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux
495	! using s_i_1d and not sz_i_1d(jk) is ok
496	wfx_res_1d(ji) = wfx_res_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux
497
498	! update heat content (J.m-2) and layer thickness
499	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
500	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
501
502	ELSE !-- Basal melting
503
504	zEi = - e_i_1d(ji,jk) * r1_rhoi ! Specific enthalpy of melting ice (J/kg, <0)
505	zEw = rcp * ztmelts ! Specific enthalpy of meltwater (J/kg, <0)
506	zdE = zEi - zEw ! Specific enthalpy difference (J/kg, <0)
507
508	zfmdt = - zq_bot(ji) / zdE ! Mass flux x time step (kg/m2, >0)
509
510	zdeltah(ji,jk) = - zfmdt * r1_rhoi ! Gross thickness change
511
512	zdeltah(ji,jk) = MIN( 0._wp , MAX( zdeltah(ji,jk), - zh_i(ji,jk) ) ) ! bound thickness change
513
514	zq_bot(ji) = MAX( 0._wp , zq_bot(ji) - zdeltah(ji,jk) * rhoi * zdE ) ! update available heat. MAX is necessary for roundup errors
515
516	dh_i_bom(ji) = dh_i_bom(ji) + zdeltah(ji,jk) ! Update basal melt
517
518	zfmdt = - zdeltah(ji,jk) * rhoi ! Mass flux x time step > 0
519
520	zQm = zfmdt * zEw ! Heat exchanged with ocean
521
522	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux to the ocean [W.m-2], <0
523	hfx_bom_1d(ji) = hfx_bom_1d(ji) - zfmdt * a_i_1d(ji) * zdE * r1_Dt_ice ! Heat used in this process [W.m-2], >0
524
525	sfx_bom_1d(ji) = sfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * s_i_1d(ji) * r1_Dt_ice ! Salt flux
526	! using s_i_1d and not sz_i_1d(jk) is ok
527
528	wfx_bom_1d(ji) = wfx_bom_1d(ji) - rhoi * a_i_1d(ji) * zdeltah(ji,jk) * r1_Dt_ice ! Mass flux
529
530	! update heat content (J.m-2) and layer thickness
531	eh_i_old(ji,jk) = eh_i_old(ji,jk) + zdeltah(ji,jk) * e_i_1d(ji,jk)
532	h_i_old (ji,jk) = h_i_old (ji,jk) + zdeltah(ji,jk)
533	ENDIF
534
535	ENDIF
536	END DO
537	END DO
538
539	! Update temperature, energy
540	! --------------------------
541	DO ji = 1, npti
542	h_i_1d(ji) = MAX( 0._wp , h_i_1d(ji) + dh_i_bog(ji) + dh_i_bom(ji) )
543	END DO
544
545	! If heat still available then melt more snow
546	!-------------------------------------------
547	zdeltah(1:npti,:) = 0._wp ! important
548	DO ji = 1, npti
549	zq_rema (ji) = zq_top(ji) + zq_bot(ji)
550	rswitch = 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ! =1 if snow
551	rswitch = rswitch * MAX( 0._wp, SIGN( 1._wp, e_s_1d(ji,1) - epsi20 ) )
552	zdeltah (ji,1) = - rswitch * zq_rema(ji) / MAX( e_s_1d(ji,1), epsi20 )
553	zdeltah (ji,1) = MIN( 0._wp , MAX( zdeltah(ji,1) , - h_s_1d(ji) ) ) ! bound melting
554	dh_s_tot(ji) = dh_s_tot(ji) + zdeltah(ji,1)
555	h_s_1d (ji) = h_s_1d (ji) + zdeltah(ji,1)
556
557	zq_rema(ji) = zq_rema(ji) + zdeltah(ji,1) * e_s_1d(ji,1) ! update available heat (J.m-2)
558	hfx_snw_1d(ji) = hfx_snw_1d(ji) - zdeltah(ji,1) * a_i_1d(ji) * e_s_1d(ji,1) * r1_Dt_ice ! Heat used to melt snow, W.m-2 (>0)
559	wfx_snw_sum_1d(ji) = wfx_snw_sum_1d(ji) - rhos * a_i_1d(ji) * zdeltah(ji,1) * r1_Dt_ice ! Mass flux
560	dh_s_mlt(ji) = dh_s_mlt(ji) + zdeltah(ji,1)
561	!
562	! Remaining heat flux (W.m-2) is sent to the ocean heat budget
563	qt_oce_ai_1d(ji) = qt_oce_ai_1d(ji) + ( zq_rema(ji) * a_i_1d(ji) ) * r1_Dt_ice
564
565	IF( ln_icectl .AND. zq_rema(ji) < 0. .AND. lwp ) WRITE(numout,*) 'ALERTE zq_rema <0 = ', zq_rema(ji)
566	END DO
567
568	!
569	! Snow-Ice formation
570	! ------------------
571	! When snow load excesses Archimede's limit, snow-ice interface goes down under sea-level,
572	! flooding of seawater transforms snow into ice dh_snowice is positive for the ice
573	z1_rho = 1._wp / ( rhos+rho0-rhoi )
574	DO ji = 1, npti
575	!
576	dh_snowice(ji) = MAX( 0._wp , ( rhos * h_s_1d(ji) + (rhoi-rho0) * h_i_1d(ji) ) * z1_rho )
577
578	h_i_1d(ji) = h_i_1d(ji) + dh_snowice(ji)
579	h_s_1d(ji) = h_s_1d(ji) - dh_snowice(ji)
580
581	! Contribution to energy flux to the ocean [J/m2], >0 (if sst<0)
582	zfmdt = ( rhos - rhoi ) * dh_snowice(ji) ! <0
583	zEw = rcp * sst_1d(ji)
584	zQm = zfmdt * zEw
585
586	hfx_thd_1d(ji) = hfx_thd_1d(ji) + zfmdt * a_i_1d(ji) * zEw * r1_Dt_ice ! Heat flux
587
588	sfx_sni_1d(ji) = sfx_sni_1d(ji) + sss_1d(ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice ! Salt flux
589
590	! Case constant salinity in time: virtual salt flux to keep salinity constant
591	IF( nn_icesal /= 2 ) THEN
592	sfx_bri_1d(ji) = sfx_bri_1d(ji) - sss_1d (ji) * a_i_1d(ji) * zfmdt * r1_Dt_ice & ! put back sss_m into the ocean
593	& - s_i_1d(ji) * a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice ! and get rn_icesal from the ocean
594	ENDIF
595
596	! Mass flux: All snow is thrown in the ocean, and seawater is taken to replace the volume
597	wfx_sni_1d(ji) = wfx_sni_1d(ji) - a_i_1d(ji) * dh_snowice(ji) * rhoi * r1_Dt_ice
598	wfx_snw_sni_1d(ji) = wfx_snw_sni_1d(ji) + a_i_1d(ji) * dh_snowice(ji) * rhos * r1_Dt_ice
599
600	! update heat content (J.m-2) and layer thickness
601	eh_i_old(ji,0) = eh_i_old(ji,0) + dh_snowice(ji) * e_s_1d(ji,1) + zfmdt * zEw
602	h_i_old (ji,0) = h_i_old (ji,0) + dh_snowice(ji)
603
604	END DO
605
606	!
607	! Update temperature, energy
608	! --------------------------
609	DO ji = 1, npti
610	rswitch = 1._wp - MAX( 0._wp , SIGN( 1._wp , - h_i_1d(ji) ) )
611	t_su_1d(ji) = rswitch * t_su_1d(ji) + ( 1._wp - rswitch ) * rt0
612	END DO
613
614	DO jk = 1, nlay_s
615	DO ji = 1,npti
616	! where there is no ice or no snow
617	rswitch = ( 1._wp - MAX( 0._wp, SIGN( 1._wp, - h_s_1d(ji) ) ) ) * ( 1._wp - MAX( 0._wp, SIGN(1._wp, - h_i_1d(ji) ) ) )
618	! mass & energy loss to the ocean
619	hfx_res_1d(ji) = hfx_res_1d(ji) + ( 1._wp - rswitch ) * &
620	& ( e_s_1d(ji,jk) * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! heat flux to the ocean [W.m-2], < 0
621	wfx_res_1d(ji) = wfx_res_1d(ji) + ( 1._wp - rswitch ) * &
622	& ( rhos * h_s_1d(ji) * r1_nlay_s * a_i_1d(ji) * r1_Dt_ice ) ! mass flux
623	! update energy (mass is updated in the next loop)
624	e_s_1d(ji,jk) = rswitch * e_s_1d(ji,jk)
625	! recalculate t_s_1d from e_s_1d
626	t_s_1d(ji,jk) = rt0 + rswitch * ( - e_s_1d(ji,jk) * r1_rhos * r1_rcpi + rLfus * r1_rcpi )
627	END DO
628	END DO
629
630	! --- ensure that a_i = 0 & h_s = 0 where h_i = 0 ---
631	WHERE( h_i_1d(1:npti) == 0._wp )
632	a_i_1d(1:npti) = 0._wp
633	h_s_1d(1:npti) = 0._wp
634	END WHERE
635	!
636	END SUBROUTINE ice_thd_dh
637
638
639	!!--------------------------------------------------------------------------
640	!! INTERFACE ice_thd_snwblow
641	!!
642	!! ** Purpose : Compute distribution of precip over the ice
643	!!
644	!! Snow accumulation in one thermodynamic time step
645	!! snowfall is partitionned between leads and ice.
646	!! If snow fall was uniform, a fraction (1-at_i) would fall into leads
647	!! but because of the winds, more snow falls on leads than on sea ice
648	!! and a greater fraction (1-at_i)^beta of the total mass of snow
649	!! (beta < 1) falls in leads.
650	!! In reality, beta depends on wind speed,
651	!! and should decrease with increasing wind speed but here, it is
652	!! considered as a constant. an average value is 0.66
653	!!--------------------------------------------------------------------------
654	!!gm I think it can be usefull to set this as a FUNCTION, not a SUBROUTINE....
655	SUBROUTINE ice_thd_snwblow_2d( pin, pout )
656	REAL(wp), DIMENSION(:,:), INTENT(in ) :: pin ! previous fraction lead ( 1. - a_i_b )
657	REAL(wp), DIMENSION(:,:), INTENT(inout) :: pout
658	pout = ( 1._wp - ( pin )**rn_blow_s )
659	END SUBROUTINE ice_thd_snwblow_2d
660
661	SUBROUTINE ice_thd_snwblow_1d( pin, pout )
662	REAL(wp), DIMENSION(:), INTENT(in ) :: pin
663	REAL(wp), DIMENSION(:), INTENT(inout) :: pout
664	pout = ( 1._wp - ( pin )**rn_blow_s )
665	END SUBROUTINE ice_thd_snwblow_1d
666
667	#else
668	!!----------------------------------------------------------------------
669	!! Default option NO SI3 sea-ice model
670	!!----------------------------------------------------------------------
671	#endif
672
673	!!======================================================================
674	END MODULE icethd_dh

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

Download in other formats: