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/UKMO/NEMO_4.0.1_fix_cpl/src/ICE – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.1_fix_cpl/src/ICE/icethd_dh.F90 @ 12484

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

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

Download in other formats: