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/trunk/src/ICE – NEMO

Context Navigation

source: NEMO/trunk/src/ICE/icethd_dh.F90 @ 9598

Last change on this file since 9598 was 9598, checked in by nicolasmartin, 6 years ago
Reorganisation plan for NEMO repository: changes to make compilation succeed with new structure Juste one issue left with AGRIF_NORDIC with AGRIF preprocessing Standardisation of routines header with version 4.0 and year 2018 Fix for some broken symlinks
File size: 35.6 KB

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

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

Download in other formats: