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.

limitd_th.F90 in branches/dev_002_LIM/NEMO/LIM_SRC_3 – NEMO

Context Navigation

source: branches/dev_002_LIM/NEMO/LIM_SRC_3/limitd_th.F90 @ 825

Last change on this file since 825 was 825, checked in by ctlod, 16 years ago
dev_002_LIM : add the LIM 3.0 component, see ticketr: #71
File size: 57.3 KB

Line
1	MODULE limitd_th
2	#if defined key_lim3
3	!!======================================================================
4	!! * MODULE limitd_th *
5	!! Thermodynamics of ice thickness distribution
6	!! computation of changes in g(h)
7	!!======================================================================
8
9	!!----------------------------------------------------------------------
10	!! * Modules used
11	USE dom_ice
12	USE par_oce ! ocean parameters
13	USE dom_oce
14	USE phycst ! physical constants (ocean directory)
15	USE ice_oce ! ice variables
16	USE thd_ice
17	USE limicepoints
18	USE limistate
19	USE in_out_manager
20	USE ice
21	USE par_ice
22	USE limthd_lab
23	USE limthd_lac
24	USE limvar
25	USE iceini
26	USE limcons
27
28	IMPLICIT NONE
29	PRIVATE
30
31	!! * Routine accessibility
32	PUBLIC lim_itd_th ! called by ice_stp
33	PUBLIC lim_itd_th_rem
34	PUBLIC lim_itd_th_reb
35	PUBLIC lim_itd_fitline
36	PUBLIC lim_itd_shiftice
37
38	!! * Module variables
39	REAL(wp) :: & ! constant values
40	epsi20 = 1e-20 , &
41	epsi13 = 1e-13 , &
42	zzero = 0.e0 , &
43	zone = 1.e0
44
45	!!----------------------------------------------------------------------
46	!! LIM-@ 4.0, UCL-ASTR (2005)
47	!! (c) UCL-ASTR and Martin Vancoppenolle
48	!!----------------------------------------------------------------------
49
50	!!----------------------------------------------------------------------------------------------
51	!!----------------------------------------------------------------------------------------------
52
53	CONTAINS
54
55	SUBROUTINE lim_itd_th
56	!!------------------------------------------------------------------
57	!! * ROUTINE lim_itd_th *
58	!! ** Purpose :
59	!! This routine computes the thermodynamics of ice thickness
60	!! distribution
61	!! ** Method :
62	!!
63	!! ** Arguments :
64	!! kideb , kiut : Starting and ending points on which the
65	!! the computation is applied
66	!!
67	!! ** Inputs / Ouputs : (global commons)
68	!!
69	!! ** External :
70	!!
71	!! ** References :
72	!!
73	!! ** History :
74	!! (12-2005) Martin Vancoppenolle
75	!! Au moment ou j'ecris ces lignes, je ne me rends pas
76	!! compte du boulot que j'entame... un truc de malate comme
77	!! on dit ici chez les Belches
78	!!
79	!!------------------------------------------------------------------
80	!! * Arguments
81
82	!! * Local variables
83	INTEGER :: ji, & ! spatial dummy loop index
84	jj, & ! spatial dummy loop index
85	jk, & ! vertical layering dummy loop index
86	jl, & ! ice category dummy loop index
87	jm, & ! ice types dummy loop index
88	index, &
89	jbnd1, &
90	jbnd2
91
92	REAL(wp) :: & ! constant values
93	zeps = 1.0e-10, &
94	epsi10 = 1.0e-10
95
96	REAL(wp) :: & ! constant values for ice enthalpy
97	zindb
98
99	!!-- End of declarations
100	!!----------------------------------------------------------------------------------------------
101
102	IF(lwp) THEN
103	WRITE(numout,*)
104	WRITE(numout,*) 'lim_itd_th : Thermodynamics of the ice thickness distribution'
105	WRITE(numout,*) '~~~~~~~~~~~'
106	ENDIF
107
108	!------------------------------------------------------------------------------\|
109	! 1) Transport of ice between thickness categories. \|
110	!------------------------------------------------------------------------------\|
111	! Given thermodynamic growth rates, transport ice between
112	! thickness categories.
113	DO jm = 1, jpm
114	jbnd1 = ice_cat_bounds(jm,1)
115	jbnd2 = ice_cat_bounds(jm,2)
116	IF (ice_ncat_types(jm) .GT. 1 ) CALL lim_itd_th_rem(jbnd1, jbnd2, jm)
117	END DO
118
119	CALL lim_var_glo2eqv ! only for info
120	CALL lim_var_agg(1)
121
122	!+++++
123	WRITE(numout,*) ' From limitd_th : '
124	WRITE(numout,*) ' at_i : ', at_i(jiindex,jjindex)
125	WRITE(numout,*) ' vt_i : ', vt_i(jiindex,jjindex)
126	WRITE(numout,*) ' vt_s : ', vt_s(jiindex,jjindex)
127	DO jl = 1, jpl
128	WRITE(numout,) ' - category number ', jl
129	WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl)
130	WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl)
131	WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,jl)
132	WRITE(numout,*) ' v_s : ', v_s(jiindex,jjindex,jl)
133	WRITE(numout,*) ' ht_s : ', ht_s(jiindex,jjindex,jl)
134	WRITE(numout,*) ' e_s : ', e_s(jiindex,jjindex,1,jl)/1.0e9
135	WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,1:nlay_i,jl)/1.0e9
136	WRITE(numout,*) ' t_su : ', t_su(jiindex,jjindex,jl)
137	WRITE(numout,*) ' t_snow : ', t_s(jiindex,jjindex,1,jl)
138	WRITE(numout,*) ' t_i : ', t_i(jiindex,jjindex,1:nlay_i,jl)
139	WRITE(numout,*) ' smv_i : ', smv_i(jiindex,jjindex,jl)
140	WRITE(numout,*) ' oa_i : ', oa_i(jiindex,jjindex,jl)
141	WRITE(numout,*)
142	END DO
143	!+++++
144
145	!------------------------------------------------------------------------------\|
146	! 2) Melt ice laterally.
147	!------------------------------------------------------------------------------\|
148	! DO jm = 1, jpm
149	! CALL lim_thd_lab(ice_cat_bounds(jm,1),ice_cat_bounds(jm,2))
150	! END DO
151	! CALL lim_thd_lab
152
153	!------------------------------------------------------------------------------\|
154	! 3) Add frazil ice growing in leads.
155	!------------------------------------------------------------------------------\|
156
157	CALL lim_thd_lac
158	CALL lim_var_glo2eqv ! only for info
159
160	!+++++
161	WRITE(numout,) ' limthd_lac, new values **** '
162	DO jl = 1, jpl
163	WRITE(numout,) ' - category number ', jl
164	WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,jl)
165	WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,jl)
166	WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,jl)
167	WRITE(numout,*) ' v_s : ', v_s(jiindex,jjindex,jl)
168	WRITE(numout,*) ' e_i : ', e_i(jiindex,jjindex,1:nlay_i,jl)/1.0e9
169	WRITE(numout,*) ' smv_i : ', smv_i(jiindex,jjindex,jl)
170	WRITE(numout,*) ' t_su : ', t_su(jiindex,jjindex,jl)
171	WRITE(numout,*) ' t_snow : ', t_s(jiindex,jjindex,1,jl)
172	WRITE(numout,*) ' t_i : ', t_i(jiindex,jjindex,1:nlay_i,jl)
173	WRITE(numout,*)
174	END DO
175	!+++++
176	!----------------------------------------------------------------------------------------
177	! 4) Computation of trend terms and get back to old values
178	!----------------------------------------------------------------------------------------
179
180	!- Trend terms
181	d_a_i_thd (:,:,:) = a_i(:,:,:) - old_a_i(:,:,:)
182	d_v_s_thd (:,:,:) = v_s(:,:,:) - old_v_s(:,:,:)
183	d_v_i_thd (:,:,:) = v_i(:,:,:) - old_v_i(:,:,:)
184	d_e_s_thd(:,:,:,:) = e_s(:,:,:,:) - old_e_s(:,:,:,:)
185	d_e_i_thd(:,:,:,:) = e_i(:,:,:,:) - old_e_i(:,:,:,:)
186
187	d_smv_i_thd(:,:,:) = 0.0
188	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
189	d_smv_i_thd(:,:,:) = smv_i(:,:,:) - old_smv_i(:,:,:)
190
191	!- Recover Old values
192	a_i(:,:,:) = old_a_i (:,:,:)
193	v_s(:,:,:) = old_v_s (:,:,:)
194	v_i(:,:,:) = old_v_i (:,:,:)
195	e_s(:,:,:,:) = old_e_s (:,:,:,:)
196	e_i(:,:,:,:) = old_e_i (:,:,:,:)
197
198	IF ( ( num_sal .EQ. 2 ) .OR. ( num_sal .EQ. 4 ) ) &
199	smv_i(:,:,:) = old_smv_i (:,:,:)
200
201	END SUBROUTINE lim_itd_th
202
203	!!----------------------------------------------------------------------------------------------
204	!!----------------------------------------------------------------------------------------------
205
206	SUBROUTINE lim_itd_th_rem(klbnd,kubnd,ntyp)
207	!!------------------------------------------------------------------
208	!! * ROUTINE lim_itd_th_rem *
209	!! ** Purpose :
210	!! This routine computes the redistribution of ice thickness
211	!! after thermodynamic growth of ice thickness
212	!!
213	!! ** Method : Linear remapping
214	!!
215	!! ** Arguments :
216	!! klbnd, kubnd : Starting and ending category index on which the
217	!! the computation is applied
218	!!
219	!! ** Inputs / Ouputs : (global commons)
220	!!
221	!! ** External :
222	!!
223	!! ** References : W.H. Lipscomb, JGR 2001
224	!!
225	!! ** History :
226	!! largely inspired from CICE (c) W. H. Lipscomb and E.C. Hunke
227	!!
228	!! (01-2006) Martin Vancoppenolle, UCL-ASTR, translation from
229	!! CICE
230	!! (06-2006) Adaptation to include salt, age and types
231	!! (04-2007) Mass conservation checked
232	!!
233	!! Je suis d'humeur massacrante aujourd'hui, tout
234	!! le monde m'embete et m'empeche de coder
235	!!
236	!! Muere lentamente
237	!! quien evita una pasion y su remolino de emociones
238	!! justamente estas que regresan el brillo a los ojos
239	!! y restauran los corazones destrozados
240	!!
241	!!------------------------------------------------------------------
242	!! * Arguments
243
244	INTEGER , INTENT (IN) :: &
245	klbnd , & ! Start thickness category index point
246	kubnd , & ! End point on which the the computation is applied
247	ntyp ! Number of the type used
248
249	!! * Local variables
250	INTEGER :: ji, & ! spatial dummy loop index
251	jj, & ! spatial dummy loop index
252	jk, & ! vertical layering dummy loop index
253	jl, & ! ice category dummy loop index
254	index, & ! for ice points
255	zji, zjj, & ! dummy indices used when changing coordinates
256	nd ! used for thickness categories
257
258	INTEGER , DIMENSION(jpi,jpj,jpl-1) :: &
259	zdonor ! donor category index
260
261	REAL(wp) :: & ! constant values
262	zeps = 1.0e-10
263
264	REAL(wp) :: & ! constant values for ice enthalpy
265	zindb , &
266	zareamin , & ! minimum tolerated area in a thickness category
267	zwk1, zwk2, & ! all the following are dummy arguments
268	zx1, zx2, zx3, & !
269	zetamin , & ! minimum value of eta
270	zetamax , & ! maximum value of eta
271	zdh0 , & !
272	zda0 , & !
273	zdamax , & !
274	zhimin
275
276	REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
277	zdhice , & ! ice thickness increment
278	g0 , & ! coefficients for fitting the line of the ITD
279	g1 , & ! coefficients for fitting the line of the ITD
280	hL , & ! left boundary for the ITD for each thickness
281	hR , & ! left boundary for the ITD for each thickness
282	zht_i_o , & ! old ice thickness
283	dummy_es
284
285	REAL(wp), DIMENSION(jpi,jpj,jpl-1) :: &
286	zdaice , & ! local increment of ice area
287	zdvice ! local increment of ice volume
288
289	REAL(wp), DIMENSION(jpi,jpj,0:jpl) :: &
290	zhbnew ! new boundaries of ice categories
291
292	REAL(wp), DIMENSION(jpi,jpj) :: &
293	zhb0, zhb1 ! category boundaries for thinnes categories
294
295	REAL, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
296	zvetamin, zvetamax ! maximum values for etas
297
298	INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
299	nind_i , & ! compressed indices for i/j directions
300	nind_j
301
302	INTEGER :: &
303	nbrem ! number of cells with ice to transfer
304
305	LOGICAL, DIMENSION(jpi,jpj) :: & !:
306	zremap_flag ! compute remapping or not ????
307
308	REAL(wp) :: & ! constant values for ice enthalpy
309	zdummy, zdummy2, &
310	zslope ! used to compute local thermodynamic "speeds"
311
312	REAL (wp), DIMENSION(jpi,jpj) :: & !
313	vt_i_init, vt_i_final, & ! ice volume summed over categories
314	vt_s_init, vt_s_final, & ! snow volume summed over categories
315	et_i_init, et_i_final, & ! ice energy summed over categories
316	et_s_init, et_s_final ! snow energy summed over categories
317
318	CHARACTER (len = 15) :: fieldid
319
320	!!-- End of declarations
321	!!----------------------------------------------------------------------------------------------
322	zhimin = 0.1 !minimum ice thickness tolerated by the model
323	zareamin = zeps !minimum area in thickness categories tolerated by the conceptors of the model
324
325	!!----------------------------------------------------------------------------------------------
326	!! 0) Conservation checkand changes in each ice category
327	!!----------------------------------------------------------------------------------------------
328	IF ( con_i ) THEN
329	CALL lim_column_sum (jpl, v_i, vt_i_init)
330	CALL lim_column_sum (jpl, v_s, vt_s_init)
331	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_init)
332	dummy_es(:,:,:) = e_s(:,:,1,:)
333	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_init)
334	ENDIF
335
336	!!----------------------------------------------------------------------------------------------
337	!! 1) Compute thickness and changes in each ice category
338	!!----------------------------------------------------------------------------------------------
339	IF(lwp) THEN
340	WRITE(numout,*)
341	WRITE(numout,*) 'lim_itd_th_rem : Remapping the ice thickness distribution'
342	WRITE(numout,*) '~~~~~~~~~~~~~~~'
343	WRITE(numout,*) ' klbnd : ', klbnd
344	WRITE(numout,*) ' kubnd : ', kubnd
345	WRITE(numout,*) ' ntyp : ', ntyp
346	ENDIF
347
348	! +++++ [
349	! index = 1
350	! jiindex = arc_sp_grid(index,1)
351	! jjindex = arc_sp_grid(index,2)
352	! WRITE(numout,) '', arc_sp_acro(index), ' ', arc_sp_name(index)
353	! WRITE(numout,*)
354	! WRITE(numout,*) ' a_i : ', a_i(jiindex,jjindex,klbnd:kubnd)
355	! WRITE(numout,*) ' ht_i : ', ht_i(jiindex,jjindex,klbnd:kubnd)
356	! WRITE(numout,*) ' v_i : ', v_i(jiindex,jjindex,klbnd:kubnd)
357	! +++++ ]
358
359	zdhice(:,:,:) = 0.0
360	DO jl = klbnd, kubnd
361	DO jj = 1, jpj
362	DO ji = 1, jpi
363
364	zindb = 1.0-MAX(0.0,SIGN(1.0,-a_i(ji,jj,jl))) !0 if no ice and 1 if yes
365	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / MAX(a_i(ji,jj,jl),zeps) * zindb
366	zindb = 1.0-MAX(0.0,SIGN(1.0,-old_a_i(ji,jj,jl))) !0 if no ice and 1 if yes
367	zht_i_o(ji,jj,jl) = old_v_i(ji,jj,jl) / MAX(old_a_i(ji,jj,jl),zeps) * zindb
368	IF (a_i(ji,jj,jl).gt.1e-6) THEN
369	zdhice(ji,jj,jl) = ht_i(ji,jj,jl) - zht_i_o(ji,jj,jl)
370	ENDIF
371
372	END DO
373	END DO
374	END DO
375	!+++++
376	! WRITE(numout,) ' 1 ** '
377	! WRITE(numout,*) ' klbnd -> kubnd ', klbnd, kubnd
378	! WRITE(numout,*) 'ht_i ', ht_i (jiindex,jjindex,klbnd:kubnd)
379	! WRITE(numout,*) 'zht_i_o', zht_i_o(jiindex,jjindex,klbnd:kubnd)
380	! WRITE(numout,*) 'zdhice ', zdhice (jiindex,jjindex,klbnd:kubnd)
381	!+++++
382
383	!-----------------------------------------------------------------------------------------------
384	! 2) Compute fractional ice area in each grid cell
385	!-----------------------------------------------------------------------------------------------
386	at_i(:,:) = 0.0
387	DO jl = klbnd, kubnd
388	DO jj = 1, jpj
389	DO ji = 1, jpi
390	at_i(ji,jj) = at_i(ji,jj) + a_i(ji,jj,jl)
391	END DO
392	END DO
393	END DO
394
395	!+++++
396	! WRITE(numout,) ' 2 ** '
397	! WRITE(numout,*) ' klbnd -> kubnd ', klbnd, kubnd
398	! WRITE(numout,*) 'a_i ', a_i (jiindex,jjindex,klbnd:kubnd)
399	! WRITE(numout,*) 'at_i ', at_i (jiindex,jjindex)
400	!+++++
401
402	!-----------------------------------------------------------------------------------------------
403	! 3) Identify grid cells with ice
404	!-----------------------------------------------------------------------------------------------
405	nbrem = 0
406	DO jj = 1, jpj
407	DO ji = 1, jpi
408	IF ( at_i(ji,jj) .gt. zareamin ) THEN
409	nbrem = nbrem + 1
410	nind_i(nbrem) = ji
411	nind_j(nbrem) = jj
412	zremap_flag(ji,jj) = .true.
413	ELSE
414	zremap_flag(ji,jj) = .false.
415	ENDIF
416	END DO !ji
417	END DO !jj
418
419	!-----------------------------------------------------------------------------------------------
420	! 4) Compute new category boundaries
421	!-----------------------------------------------------------------------------------------------
422	!- 4.1 Compute category boundaries
423	hi_max(kubnd) = 999.99
424	zhbnew(:,:,:) = 0.0
425
426	DO jl = klbnd, kubnd - 1
427	! jl
428	DO ji = 1, nbrem
429	! jl, ji
430	zji = nind_i(ji)
431	zjj = nind_j(ji)
432	!
433	IF ( ( zht_i_o(zji,zjj,jl) .GT.zeps ) .AND. &
434	( zht_i_o(zji,zjj,jl+1).GT.zeps ) ) THEN
435	!interpolate between adjacent category growth rates
436	zslope = ( zdhice(zji,zjj,jl+1) - zdhice(zji,zjj,jl) ) / &
437	( zht_i_o (zji,zjj,jl+1) - zht_i_o (zji,zjj,jl) )
438	zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl) + &
439	zslope * ( hi_max(jl) - zht_i_o(zji,zjj,jl) )
440	ELSEIF (zht_i_o(zji,zjj,jl).gt.zeps) THEN
441	zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl)
442	ELSEIF (zht_i_o(zji,zjj,jl+1).gt.zeps) THEN
443	zhbnew(zji,zjj,jl) = hi_max(jl) + zdhice(zji,zjj,jl+1)
444	ELSE
445	zhbnew(zji,zjj,jl) = hi_max(jl)
446	ENDIF
447	! jl, ji
448	END DO !ji
449	! jl
450
451	!- 4.2 Check that each zhbnew lies between adjacent values of ice thickness
452	DO ji = 1, nbrem
453	! jl, ji
454	zji = nind_i(ji)
455	zjj = nind_j(ji)
456	! jl, ji
457	IF ( ( a_i(zji,zjj,jl) .GT.zeps) .AND. &
458	( ht_i(zji,zjj,jl).GE. zhbnew(zji,zjj,jl) ) &
459	) THEN
460	zremap_flag(zji,zjj) = .false.
461	ELSEIF ( ( a_i(zji,zjj,jl+1) .GT. zeps ) .AND. &
462	( ht_i(zji,zjj,jl+1).LE. zhbnew(zji,zjj,jl) ) &
463	) THEN
464	zremap_flag(zji,zjj) = .false.
465	ENDIF
466
467	!- 4.3 Check that each zhbnew does not exceed maximal values hi_max
468	! jl, ji
469	IF (zhbnew(zji,zjj,jl).gt.hi_max(jl+1)) THEN
470	zremap_flag(zji,zjj) = .false.
471	ENDIF
472	! jl, ji
473	IF (zhbnew(zji,zjj,jl).lt.hi_max(jl-1)) THEN
474	zremap_flag(zji,zjj) = .false.
475	ENDIF
476	! jl, ji
477	END DO !ji
478	! ji
479	END DO !jl
480	!+++++
481	! WRITE(numout,) ' 4 ** '
482	! WRITE(numout,*) ' klbnd -> kubnd - 1 ', klbnd, kubnd - 1
483	! WRITE(numout,*) ' hi_max ', hi_max(klbnd:kubnd-1)
484	! WRITE(numout,*) ' zhbnew ', zhbnew(jiindex,jjindex,klbnd:kubnd-1)
485	!+++++
486
487	!-----------------------------------------------------------------------------------------------
488	! 5) Identify cells where ITD is to be remapped
489	!-----------------------------------------------------------------------------------------------
490	nbrem = 0
491	DO jj = 1, jpj
492	DO ji = 1, jpi
493	IF ( zremap_flag(ji,jj) ) THEN
494	nbrem = nbrem + 1
495	nind_i(nbrem) = ji
496	nind_j(nbrem) = jj
497	ENDIF
498	END DO !ji
499	END DO !jj
500
501	!-----------------------------------------------------------------------------------------------
502	! 6) Fill arrays with lowermost / uppermost boundaries of 'new' categories
503	!-----------------------------------------------------------------------------------------------
504	DO jj = 1, jpj
505	DO ji = 1, jpi
506	zhb0(ji,jj) = hi_max_typ(0,ntyp) ! 0eme
507	zhb1(ji,jj) = hi_max_typ(1,ntyp) ! 1er
508
509	zhbnew(ji,jj,klbnd-1) = 0.0
510
511	IF ( a_i(ji,jj,kubnd) .GT. zeps ) THEN
512	zhbnew(ji,jj,kubnd) = 3.0ht_i(ji,jj,kubnd) - 2.0zhbnew(ji,jj,kubnd-1)
513	ELSE
514	zhbnew(ji,jj,kubnd) = hi_max(kubnd)
515	ENDIF
516
517	IF ( zhbnew(ji,jj,kubnd) .LT. hi_max(kubnd-1) ) &
518	zhbnew(ji,jj,kubnd) = hi_max(kubnd-1)
519
520	END DO !jj
521	END DO !jj
522
523	!+++++
524	! WRITE(numout,) ' 6 ** '
525	! WRITE(numout,*) ' klbnd -1, kubnd ', klbnd - 1, kubnd
526	! WRITE(numout,*) ' zhb0 ', zhb0(jiindex,jjindex)
527	! WRITE(numout,*) ' zhb1 ', zhb1(jiindex,jjindex)
528	! WRITE(numout,*) ' zhbnew klbnd-1 ', zhbnew(jiindex,jjindex,klbnd-1)
529	! WRITE(numout,*) ' zhbnew kubnd ', zhbnew(jiindex,jjindex,klbnd)
530	!+++++
531
532	!-----------------------------------------------------------------------------------------------
533	! 7) Compute g(h)
534	!-----------------------------------------------------------------------------------------------
535	!- 7.1 g(h) for category 1 at start of time step
536	CALL lim_itd_fitline(klbnd, zhb0, zhb1, zht_i_o(:,:,klbnd), &
537	g0(:,:,klbnd), g1(:,:,klbnd), hL(:,:,klbnd), &
538	hR(:,:,klbnd), zremap_flag)
539
540	!+++++
541	! WRITE(numout,) ' 7a ** klbnd ', klbnd
542	! WRITE(numout,*) ' g0(klbnd) ', g0(jiindex,jjindex,klbnd)
543	! WRITE(numout,*) ' g1(klbnd) ', g1(jiindex,jjindex,klbnd)
544	! WRITE(numout,*) ' hL(klbnd) ', hL(jiindex,jjindex,klbnd)
545	! WRITE(numout,*) ' hR(klbnd) ', hR(jiindex,jjindex,klbnd)
546	!+++++
547
548	!- 7.2 Area lost due to melting of thin ice (first category, klbnd)
549	DO ji = 1, nbrem
550	zji = nind_i(ji)
551	zjj = nind_j(ji)
552
553	!ji
554	IF (a_i(zji,zjj,klbnd) .gt. zeps) THEN
555	zdh0 = zdhice(zji,zjj,klbnd) !decrease of ice thickness in the lower category
556	! ji, a_i > zeps
557	IF (zdh0 .lt. 0.0) THEN !remove area from category 1
558	! ji, a_i > zeps; zdh0 < 0
559	zdh0 = MIN(-zdh0,hi_max(klbnd))
560
561	!Integrate g(1) from 0 to dh0 to estimate area melted
562	zetamax = MIN(zdh0,hR(zji,zjj,klbnd)) - hL(zji,zjj,klbnd)
563	IF (zetamax.gt.0.0) THEN
564	zx1 = zetamax
565	zx2 = 0.5 * zetamax*zetamax
566	zda0 = g1(zji,zjj,klbnd) * zx2 + g0(zji,zjj,klbnd) * zx1 !ice area removed
567	! Constrain new thickness <= ht_i
568	zdamax = a_i(zji,zjj,klbnd) * &
569	(1.0 - ht_i(zji,zjj,klbnd)/zht_i_o(zji,zjj,klbnd)) ! zdamax > 0
570	!ice area lost due to melting of thin ice
571	zda0 = MIN(zda0, zdamax)
572
573	! Remove area, conserving volume
574	ht_i(zji,zjj,klbnd) = ht_i(zji,zjj,klbnd) &
575	* a_i(zji,zjj,klbnd) / ( a_i(zji,zjj,klbnd) - zda0 )
576	a_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd) - zda0
577	v_i(zji,zjj,klbnd) = a_i(zji,zjj,klbnd)*ht_i(zji,zjj,klbnd)
578	ENDIF ! zetamax > 0
579	! ji, a_i > zeps
580
581	ELSE ! if ice accretion
582	! ji, a_i > zeps; zdh0 > 0
583	IF ( ntyp .EQ. 1 ) zhbnew(zji,zjj,klbnd-1) = MIN(zdh0,hi_max(klbnd))
584	! zhbnew was 0, and is shifted to the right to account for thin ice
585	! growth in openwater (F0 = f1)
586	IF ( ntyp .NE. 1 ) zhbnew(zji,zjj,0) = 0
587	! in other types there is
588	! no open water growth (F0 = 0)
589	ENDIF ! zdh0
590
591	! a_i > zeps
592	ENDIF ! a_i > zeps
593
594	END DO ! ji
595
596	!- 7.3 g(h) for each thickness category
597	DO jl = klbnd, kubnd
598	CALL lim_itd_fitline(jl, zhbnew(:,:,jl-1), zhbnew(:,:,jl), ht_i(:,:,jl), &
599	g0(:,:,jl), g1(:,:,jl), hL(:,:,jl), hR(:,:,jl), &
600	zremap_flag)
601	END DO
602
603	!+++++
604	! WRITE(numout,) ' 7b ** klbnd->kubnd ', klbnd, kubnd
605	! WRITE(numout,*) ' g0 ', g0(jiindex,jjindex,klbnd:kubnd)
606	! WRITE(numout,*) ' g1 ', g1(jiindex,jjindex,klbnd:kubnd)
607	! WRITE(numout,*) ' hL ', hL(jiindex,jjindex,klbnd:kubnd)
608	! WRITE(numout,*) ' hR ', hR(jiindex,jjindex,klbnd:kubnd)
609	! WRITE(numout,*)
610	! WRITE(numout,*) ' ht_i ', ht_i(jiindex,jjindex,klbnd:kubnd)
611	! WRITE(numout,*) ' a_i ', a_i (jiindex,jjindex,klbnd:kubnd)
612	! WRITE(numout,*) ' v_i ', v_i (jiindex,jjindex,klbnd:kubnd)
613	!+++++
614
615	!-----------------------------------------------------------------------------------------------
616	! 8) Compute area and volume to be shifted across each boundary
617	!-----------------------------------------------------------------------------------------------
618
619	DO jl = klbnd, kubnd - 1
620	DO jj = 1, jpj
621	DO ji = 1, jpi
622	zdonor(ji,jj,jl) = 0
623	zdaice(ji,jj,jl) = 0.0
624	zdvice(ji,jj,jl) = 0.0
625	END DO
626	END DO
627
628	DO ji = 1, nbrem
629	zji = nind_i(ji)
630	zjj = nind_j(ji)
631
632	IF (zhbnew(zji,zjj,jl) .gt. hi_max(jl)) THEN ! transfer from jl to jl+1
633
634	! left and right integration limits in eta space
635	zvetamin(ji) = MAX(hi_max(jl), hL(zji,zjj,jl)) - hL(zji,zjj,jl)
636	zvetamax(ji) = MIN(zhbnew(zji,zjj,jl), hR(zji,zjj,jl)) - hL(zji,zjj,jl)
637	zdonor(zji,zjj,jl) = jl
638
639	ELSE ! zhbnew(jl) <= hi_max(jl) ; transfer from jl+1 to jl
640
641	! left and right integration limits in eta space
642	zvetamin(ji) = 0.0
643	zvetamax(ji) = MIN(hi_max(jl), hR(zji,zjj,jl+1)) - hL(zji,zjj,jl+1)
644	zdonor(zji,zjj,jl) = jl + 1
645
646	ENDIF ! zhbnew(jl) > hi_max(jl)
647
648	zetamax = MAX(zvetamax(ji), zvetamin(ji)) ! no transfer if etamax < etamin
649	zetamin = zvetamin(ji)
650
651	zx1 = zetamax - zetamin
652	zwk1 = zetamin*zetamin
653	zwk2 = zetamax*zetamax
654	zx2 = 0.5 * (zwk2 - zwk1)
655	zwk1 = zwk1 * zetamin
656	zwk2 = zwk2 * zetamax
657	zx3 = 1.0/3.0 * (zwk2 - zwk1)
658	nd = zdonor(zji,zjj,jl)
659	zdaice(zji,zjj,jl) = g1(zji,zjj,nd)zx2 + g0(zji,zjj,nd)zx1
660	zdvice(zji,zjj,jl) = g1(zji,zjj,nd)zx3 + g0(zji,zjj,nd)zx2 + &
661	zdaice(zji,zjj,jl)*hL(zji,zjj,nd)
662
663	END DO ! ji
664	END DO ! jl klbnd -> kubnd - 1
665
666	!!----------------------------------------------------------------------------------------------
667	!! 9) Shift ice between categories
668	!!----------------------------------------------------------------------------------------------
669	CALL lim_itd_shiftice ( klbnd, kubnd, zdonor, zdaice, zdvice )
670
671	! WRITE(numout,) ' 9 ** '
672	! WRITE(numout,*) ' ht_i ', ht_i(jiindex,jjindex,klbnd:kubnd)
673	! WRITE(numout,*) ' a_i ', a_i (jiindex,jjindex,klbnd:kubnd)
674	! WRITE(numout,*) ' v_i ', v_i (jiindex,jjindex,klbnd:kubnd)
675
676	!!----------------------------------------------------------------------------------------------
677	!! 10) Make sure ht_i >= minimum ice thickness hi_min
678	!!----------------------------------------------------------------------------------------------
679
680	DO ji = 1, nbrem
681	zji = nind_i(ji)
682	zjj = nind_j(ji)
683	IF ( ( zhimin .GT. 0.0 ) .AND. &
684	( ( a_i(zji,zjj,1) .GT. zeps ) .AND. ( ht_i(zji,zjj,1) .LT. zhimin ) ) &
685	) THEN
686	a_i(zji,zjj,1) = a_i(zji,zjj,1) * ht_i(zji,zjj,1) / zhimin
687	ht_i(zji,zjj,1) = zhimin
688	v_i(zji,zjj,1) = a_i(zji,zjj,1)*ht_i(zji,zjj,1)
689	ENDIF
690	END DO !ji
691
692	!!----------------------------------------------------------------------------------------------
693	!! 11) Conservation check
694	!!----------------------------------------------------------------------------------------------
695	IF ( con_i ) THEN
696	CALL lim_column_sum (jpl, v_i, vt_i_final)
697	fieldid = ' v_i : limitd_th '
698	CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
699
700	CALL lim_column_sum_energy (jpl, nlay_i, e_i, et_i_final)
701	fieldid = ' e_i : limitd_th '
702	CALL lim_cons_check (et_i_init, et_i_final, 1.0e-3, fieldid)
703
704	CALL lim_column_sum (jpl, v_s, vt_s_final)
705	fieldid = ' v_s : limitd_th '
706	CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
707
708	dummy_es(:,:,:) = e_s(:,:,1,:)
709	CALL lim_column_sum (jpl, dummy_es(:,:,:) , et_s_final)
710	fieldid = ' e_s : limitd_th '
711	CALL lim_cons_check (et_s_init, et_s_final, 1.0e-3, fieldid)
712	ENDIF
713
714	END SUBROUTINE lim_itd_th_rem
715
716	!!----------------------------------------------------------------------------------------------
717	!!----------------------------------------------------------------------------------------------
718
719	SUBROUTINE lim_itd_fitline(num_cat, HbL, Hbr, hice, g0, g1, hL, hR, zremap_flag )
720
721	!!------------------------------------------------------------------
722	!! * ROUTINE lim_itd_fitline *
723	!! ** Purpose :
724	!! fit g(h) with a line using area, volume constraints
725	!!
726	!! ** Method :
727	!! Fit g(h) with a line, satisfying area and volume constraints.
728	!! To reduce roundoff errors caused by large values of g0 and g1,
729	!! we actually compute g(eta), where eta = h - hL, and hL is the
730	!! left boundary.
731	!!
732	!! ** Arguments :
733	!!
734	!! ** Inputs / Ouputs : (global commons)
735	!!
736	!! ** External :
737	!!
738	!! ** References :
739	!!
740	!! ** History :
741	!! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
742	!! (01-2006) Martin Vancoppenolle
743	!! Au moment ou j'ecris ces lignes, je ne me rends pas
744	!! compte du boulot que j'entame... un truc de malate comme
745	!! on dit ici chez les Belches
746	!! This routine was inspired from CICE (W.H. Lipscomb, E. Hunke, C. M. Bitz)
747	!! the sea ice model of LANL, Los Alamos, USA.
748	!! Thanks to these guys and their team to put their routines online
749	!!
750	!!------------------------------------------------------------------
751	!! * Arguments
752
753	INTEGER, INTENT(in) :: num_cat ! category index
754
755	REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: & !:
756	HbL, HbR ! left and right category boundaries
757
758	REAL(wp), DIMENSION(jpi,jpj), INTENT(IN) :: & !:
759	hice ! ice thickness
760
761	REAL(wp), DIMENSION(jpi,jpj), INTENT(OUT) :: & !:
762	g0, g1 , & ! coefficients in linear equation for g(eta)
763	hL , & ! min value of range over which g(h) > 0
764	hR ! max value of range over which g(h) > 0
765
766	LOGICAL, DIMENSION(jpi,jpj), INTENT(IN) :: & !:
767	zremap_flag
768
769	INTEGER :: &
770	ji,jj ! horizontal indices
771
772	REAL(wp) :: &
773	zh13 , & ! HbL + 1/3 * (HbR - HbL)
774	zh23 , & ! HbL + 2/3 * (HbR - HbL)
775	zdhr , & ! 1 / (hR - hL)
776	zwk1, zwk2 , & ! temporary variables
777	zacrith ! critical minimum concentration in an ice category
778
779	REAL(wp) :: & ! constant values
780	zeps = 1.0e-10
781
782	zacrith = 1.0e-6
783	!!-- End of declarations
784	!!----------------------------------------------------------------------------------------------
785
786	! WRITE(numout,*) ' lim_itd_fitline : linearly fitting the function g(h) '
787	! WRITE(numout,*) ' ~~~~~~~~~~~~~~~ in category number ', num_cat
788
789	DO jj = 1, jpj
790	DO ji = 1, jpi
791
792	IF ( zremap_flag(ji,jj) .AND. a_i(ji,jj,num_cat) .gt. zacrith &
793	.AND. hice(ji,jj) .GT. 0.0 ) THEN
794
795	! Initialize hL and hR
796
797	hL(ji,jj) = HbL(ji,jj)
798	hR(ji,jj) = HbR(ji,jj)
799
800	! Change hL or hR if hice falls outside central third of range
801
802	zh13 = 1.0/3.0 * (2.0*hL(ji,jj) + hR(ji,jj))
803	zh23 = 1.0/3.0 * (hL(ji,jj) + 2.0*hR(ji,jj))
804
805	IF (hice(ji,jj) < zh13) THEN
806	hR(ji,jj) = 3.0hice(ji,jj) - 2.0hL(ji,jj)
807	ELSEIF (hice(ji,jj) > zh23) THEN
808	hL(ji,jj) = 3.0hice(ji,jj) - 2.0hR(ji,jj)
809	ENDIF
810
811	! Compute coefficients of g(eta) = g0 + g1*eta
812
813	zdhr = 1.0 / (hR(ji,jj) - hL(ji,jj))
814	zwk1 = 6.0 * a_i(ji,jj,num_cat) * zdhr
815	zwk2 = (hice(ji,jj) - hL(ji,jj)) * zdhr
816	g0(ji,jj) = zwk1 * (2.0/3.0 - zwk2)
817	g1(ji,jj) = 2.0zdhr zwk1 * (zwk2 - 0.5)
818
819	ELSE ! remap_flag = .false. or a_i < zeps
820
821	hL(ji,jj) = 0.0
822	hR(ji,jj) = 0.0
823	g0(ji,jj) = 0.0
824	g1(ji,jj) = 0.0
825
826	ENDIF ! a_i > zeps
827
828	END DO !ji
829	END DO ! jj
830
831	END SUBROUTINE lim_itd_fitline
832
833	!----------------------------------------------------------------------------------------------
834	!----------------------------------------------------------------------------------------------
835
836	SUBROUTINE lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice)
837	!!------------------------------------------------------------------
838	!! * ROUTINE lim_itd_shiftice *
839	!! ** Purpose : shift ice across category boundaries, conserving everything
840	!! ( area, volume, energy, age*vol, and mass of salt )
841	!!
842	!! ** Method :
843	!!
844	!! ** Arguments :
845	!!
846	!! ** Inputs / Ouputs : (global commons)
847	!!
848	!! ** External :
849	!!
850	!! ** References :
851	!!
852	!! ** History :
853	!! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
854	!! (01-2006) Martin Vancoppenolle
855	!! Et quand fatigues de s'etre souvenus
856	!! Nos souvenirs fatigues ne seront plus que des loques...
857	!! "!!! Il a attrape la ratatinette, l'epouvantable ratatinette"
858	!!
859	!! This routine was largely inspired from CICE
860	!! (W.H. Lipscomb, E. Hunke, C. M. Bitz)
861	!! the sea ice model of LANL, Los Alamos, USA.
862	!! Merci a eux et a leur equipe de mettre leurs routines en ligne
863	!!
864	!! J'ajoute aujourd'hui : le secret de ma vie pour l'instant
865	!! c'est de trouver la balance entre mon ego et l'absence d'ego
866	!! balance entre respect des autres et respect de moi
867	!! trouver comment realiser la compréhension de l'autre
868	!! sans déclencher ma propre aliénation
869	!!
870	!!------------------------------------------------------------------
871	!! * Arguments
872
873	INTEGER , INTENT (IN) :: &
874	klbnd , & ! Start thickness category index point
875	kubnd ! End point on which the the computation is applied
876
877	INTEGER , DIMENSION(jpi,jpj,jpl-1), INTENT(IN) :: &
878	zdonor ! donor category index
879
880	REAL(wp), DIMENSION(jpi,jpj,jpl-1), INTENT(INOUT) :: &
881	zdaice , & ! ice area transferred across boundary
882	zdvice ! ice volume transferred across boundary
883
884	INTEGER :: &
885	ji,jj,jl, & ! horizontal indices, thickness category index
886	jl2, & ! receiver category
887	jl1, & ! donor category
888	jk, & ! ice layer index
889	zji, zjj, & ! indices when changing from 2D-1D is done
890	index ! for ice points referencing
891
892	REAL(wp), DIMENSION(jpi,jpj,jpl) :: &
893	zaTsfn, &
894	zqsnow
895
896	REAL(wp), DIMENSION(jpi,jpj) :: &
897	zworka ! temporary array used here
898
899	REAL(wp) :: &
900	zdvsnow , & ! snow volume transferred
901	zdesnow , & ! snow energy transferred
902	zdeice , & ! ice energy transferred
903	zdsm_vice , & ! ice salinity times volume transferred
904	zsm_v1 , & ! ice salinity times volume
905	zsm_v2 , & ! ice salinity times volume
906	zdo_aice , & ! ice age times volume transferred
907	zo_v1 , & ! ice age times volume
908	zo_v2 , & ! ice age times volume
909	zdaTsf , & ! aicen*Tsfcn transferred
910	zindsn , & ! snow or not
911	zindb ! ice or not
912
913	INTEGER, DIMENSION(1:(jpi+1)*(jpj+1)) :: &
914	nind_i , & ! compressed indices for i/j directions
915	nind_j
916
917	INTEGER :: &
918	nbrem ! number of cells with ice to transfer
919
920	LOGICAL :: &
921	zdaice_negative , & ! true if daice < -puny
922	zdvice_negative , & ! true if dvice < -puny
923	zdaice_greater_aicen , & ! true if daice > aicen
924	zdvice_greater_vicen ! true if dvice > vicen
925
926	REAL(wp) :: & ! constant values
927	zeps = 1.0e-10
928
929	!!-- End of declarations
930	! WRITE(numout,*) ' lim_itd_shiftice : shifting ice between categories '
931	! WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
932	!----------------------------------------------------------------------------------------------
933	! 1) Define a variable equal to a_i*T_su
934	!----------------------------------------------------------------------------------------------
935
936	DO jl = klbnd, kubnd
937	DO jj = 1, jpj
938	DO ji = 1, jpi
939	zaTsfn(ji,jj,jl) = a_i(ji,jj,jl)*t_su(ji,jj,jl)
940	END DO ! ji
941	END DO ! jj
942	END DO ! jl
943
944	!----------------------------------------------------------------------------------------------
945	! 2) Check for daice or dvice out of range, allowing for roundoff error
946	!----------------------------------------------------------------------------------------------
947	! Note: zdaice < 0 or zdvice < 0 usually happens when category jl
948	! has a small area, with h(n) very close to a boundary. Then
949	! the coefficients of g(h) are large, and the computed daice and
950	! dvice can be in error. If this happens, it is best to transfer
951	! either the entire category or nothing at all, depending on which
952	! side of the boundary hice(n) lies.
953	!-----------------------------------------------------------------
954	DO jl = klbnd, kubnd-1
955
956	zdaice_negative = .false.
957	zdvice_negative = .false.
958	zdaice_greater_aicen = .false.
959	zdvice_greater_vicen = .false.
960
961	DO jj = 1, jpj
962	DO ji = 1, jpi
963
964	IF (zdonor(ji,jj,jl) .GT. 0) THEN
965	jl1 = zdonor(ji,jj,jl)
966
967	IF (zdaice(ji,jj,jl) .LT. 0.0) THEN
968	IF (zdaice(ji,jj,jl) .GT. -zeps) THEN
969	IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1) .GT. hi_max(jl) ) &
970	.OR. &
971	( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) &
972	) THEN
973	zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
974	zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
975	ELSE
976	zdaice(ji,jj,jl) = 0.0 ! shift no ice
977	zdvice(ji,jj,jl) = 0.0
978	ENDIF
979	ELSE
980	zdaice_negative = .true.
981	ENDIF
982	ENDIF
983
984	IF (zdvice(ji,jj,jl) .LT. 0.0) THEN
985	IF (zdvice(ji,jj,jl) .GT. -zeps ) THEN
986	IF ( ( jl1.EQ.jl .AND. ht_i(ji,jj,jl1).GT.hi_max(jl) ) &
987	.OR. &
988	( jl1.EQ.jl+1 .AND. ht_i(ji,jj,jl1) .LE. hi_max(jl) ) &
989	) THEN
990	zdaice(ji,jj,jl) = a_i(ji,jj,jl1) ! shift entire category
991	zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
992	ELSE
993	zdaice(ji,jj,jl) = 0.0 ! shift no ice
994	zdvice(ji,jj,jl) = 0.0
995	ENDIF
996	ELSE
997	zdvice_negative = .true.
998	ENDIF
999	ENDIF
1000
1001	! If daice is close to aicen, set daice = aicen.
1002	IF (zdaice(ji,jj,jl) .GT. a_i(ji,jj,jl1) - zeps ) THEN
1003	IF (zdaice(ji,jj,jl) .LT. a_i(ji,jj,jl1)+zeps) THEN
1004	zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
1005	zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
1006	ELSE
1007	zdaice_greater_aicen = .true.
1008	ENDIF
1009	ENDIF
1010
1011	IF (zdvice(ji,jj,jl) .GT. v_i(ji,jj,jl1)-zeps) THEN
1012	IF (zdvice(ji,jj,jl) .LT. v_i(ji,jj,jl1)+zeps) THEN
1013	zdaice(ji,jj,jl) = a_i(ji,jj,jl1)
1014	zdvice(ji,jj,jl) = v_i(ji,jj,jl1)
1015	ELSE
1016	zdvice_greater_vicen = .true.
1017	ENDIF
1018	ENDIF
1019
1020	ENDIF ! donor > 0
1021	END DO ! i
1022	END DO ! j
1023
1024	END DO !jl
1025
1026	!-----------------------------------------------------------------
1027	! error messages
1028	!-----------------------------------------------------------------
1029
1030	! if (daice_negative) then
1031	! do j = jlo,jhi
1032	! do i = ilo,ihi
1033	! if (donor(i,j,n) > 0 .and. daice(i,j,n) <= -puny) then
1034	! write(nu_diag,*) my_task,':',i,j,
1035	! & 'ITD Neg daice =',daice(i,j,n),' boundary',n
1036	! call abort_ice ('ice: ITD Neg daice')
1037	! endif
1038	! enddo
1039	! enddo
1040	! endif
1041
1042	! if (dvice_negative) then
1043	! do j = jlo,jhi
1044	! do i = ilo,ihi
1045	! if (donor(i,j,n) > 0 .and. dvice(i,j,n) <= -puny) then
1046	! write(nu_diag,*) my_task,':',i,j,
1047	! & 'ITD Neg dvice =',dvice(i,j,n),' boundary',n
1048	! call abort_ice ('ice: ITD Neg dvice')
1049	! endif
1050	! enddo
1051	! enddo
1052	! endif
1053
1054	! if (daice_greater_aicen) then
1055	! do j = jlo,jhi
1056	! do i = ilo,ihi
1057	! if (donor(i,j,n) > 0) then
1058	! n1 = donor(i,j,n)
1059	! if (daice(i,j,n) >= aicen(i,j,n1)+puny) then
1060	! write(nu_diag,*) my_task,':',i,j,
1061	! & 'ITD daice > aicen, cat',n1
1062	! write(nu_diag,*) my_task,':',i,j,
1063	! & 'daice =', daice(i,j,n),
1064	! & 'aicen =', aicen(i,j,n1)
1065	! call abort_ice ('ice: ITD daice > aicen')
1066	! endif
1067	! endif
1068	! enddo
1069	! enddo
1070	! endif
1071
1072	! if (dvice_greater_vicen) then
1073	! do j = jlo,jhi
1074	! do i = ilo,ihi
1075	! if (donor(i,j,n) > 0) then
1076	! n1 = donor(i,j,n)
1077	! if (dvice(i,j,n) >= vicen(i,j,n1)+puny) then
1078	! write(nu_diag,*) my_task,':',i,j,
1079	! & 'ITD dvice > vicen, cat',n1
1080	! write(nu_diag,*) my_task,':',i,j,
1081	! & 'dvice =', dvice(i,j,n),
1082	! & 'vicen =', vicen(i,j,n1)
1083	! call abort_ice ('ice: ITD dvice > vicen')
1084	! endif
1085	! endif
1086	! enddo
1087	! enddo
1088	! endif
1089
1090	! enddo ! boundaries 1 to ncat-1
1091
1092	!-------------------------------------------------------------------------------
1093	! 3) Transfer volume and energy between categories
1094	!-------------------------------------------------------------------------------
1095
1096	DO jl = klbnd, kubnd - 1
1097	nbrem = 0
1098	DO jj = 1, jpj
1099	DO ji = 1, jpi
1100	IF (zdaice(ji,jj,jl) .GT. 0.0 ) THEN ! daice(n) can be < puny
1101	nbrem = nbrem + 1
1102	nind_i(nbrem) = ji
1103	nind_j(nbrem) = jj
1104	ENDIF ! tmask
1105	END DO
1106	END DO
1107
1108	DO ji = 1, nbrem
1109	zji = nind_i(ji)
1110	zjj = nind_j(ji)
1111
1112	jl1 = zdonor(zji,zjj,jl)
1113	zindb = MAX( 0.0 , SIGN( 1.0 , v_i(zji,zjj,jl1) - zeps ) )
1114	zworka(zji,zjj) = zdvice(zji,zjj,jl) / MAX(v_i(zji,zjj,jl1),zeps) * zindb
1115	IF (jl1 .eq. jl) THEN
1116	jl2 = jl1+1
1117	ELSE ! n1 = n+1
1118	jl2 = jl
1119	ENDIF
1120
1121	!--------------
1122	! Ice areas
1123	!--------------
1124
1125	a_i(zji,zjj,jl1) = a_i(zji,zjj,jl1) - zdaice(zji,zjj,jl)
1126	a_i(zji,zjj,jl2) = a_i(zji,zjj,jl2) + zdaice(zji,zjj,jl)
1127
1128	!--------------
1129	! Ice volumes
1130	!--------------
1131
1132	v_i(zji,zjj,jl1) = v_i(zji,zjj,jl1) - zdvice(zji,zjj,jl)
1133	v_i(zji,zjj,jl2) = v_i(zji,zjj,jl2) + zdvice(zji,zjj,jl)
1134
1135	!--------------
1136	! Snow volumes
1137	!--------------
1138
1139	zdvsnow = v_s(zji,zjj,jl1) * zworka(zji,zjj)
1140	v_s(zji,zjj,jl1) = v_s(zji,zjj,jl1) - zdvsnow
1141	v_s(zji,zjj,jl2) = v_s(zji,zjj,jl2) + zdvsnow
1142
1143	!--------------------
1144	! Snow heat content
1145	!--------------------
1146
1147	zdesnow = e_s(zji,zjj,1,jl1) * zworka(zji,zjj)
1148	e_s(zji,zjj,1,jl1) = e_s(zji,zjj,1,jl1) - zdesnow
1149	e_s(zji,zjj,1,jl2) = e_s(zji,zjj,1,jl2) + zdesnow
1150
1151	!--------------
1152	! Ice age
1153	!--------------
1154
1155	zdo_aice = oa_i(zji,zjj,jl1) * zdaice(zji,zjj,jl)
1156	oa_i(zji,zjj,jl1) = oa_i(zji,zjj,jl1) - zdo_aice
1157	oa_i(zji,zjj,jl2) = oa_i(zji,zjj,jl2) + zdo_aice
1158
1159	!--------------
1160	! Ice salinity
1161	!--------------
1162
1163	zdsm_vice = smv_i(zji,zjj,jl1) * zworka(zji,zjj)
1164	smv_i(zji,zjj,jl1) = smv_i(zji,zjj,jl1) - zdsm_vice
1165	smv_i(zji,zjj,jl2) = smv_i(zji,zjj,jl2) + zdsm_vice
1166
1167	!---------------------
1168	! Surface temperature
1169	!---------------------
1170
1171	zdaTsf = t_su(zji,zjj,jl1) * zdaice(zji,zjj,jl)
1172	zaTsfn(zji,zjj,jl1) = zaTsfn(zji,zjj,jl1) - zdaTsf
1173	zaTsfn(zji,zjj,jl2) = zaTsfn(zji,zjj,jl2) + zdaTsf
1174
1175	END DO ! ji
1176
1177	!------------------
1178	! Ice heat content
1179	!------------------
1180
1181	DO jk = 1, nlay_i
1182	DO ji = 1, nbrem
1183	zji = nind_i(ji)
1184	zjj = nind_j(ji)
1185
1186	jl1 = zdonor(zji,zjj,jl)
1187	IF (jl1 .EQ. jl) THEN
1188	jl2 = jl+1
1189	ELSE ! n1 = n+1
1190	jl2 = jl
1191	ENDIF
1192
1193	zdeice = e_i(zji,zjj,jk,jl1) * zworka(zji,zjj)
1194	e_i(zji,zjj,jk,jl1) = e_i(zji,zjj,jk,jl1) - zdeice
1195	e_i(zji,zjj,jk,jl2) = e_i(zji,zjj,jk,jl2) + zdeice
1196	END DO ! ji
1197	END DO ! jk
1198
1199	END DO ! boundaries, 1 to ncat-1
1200
1201	!-----------------------------------------------------------------
1202	! Update ice thickness and temperature
1203	!-----------------------------------------------------------------
1204
1205	DO jl = klbnd, kubnd
1206	DO jj = 1, jpj
1207	DO ji = 1, jpi
1208	IF ( a_i(ji,jj,jl) .GT. zeps ) THEN
1209	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
1210	t_su(ji,jj,jl) = zaTsfn(ji,jj,jl) / a_i(ji,jj,jl)
1211	zindsn = 1.0 - MAX(0.0,SIGN(1.0,-v_s(ji,jj,jl))) !0 if no ice and 1 if yes
1212	! t_s(ji,jj,1,jl) = - zQsnow(ji,jj,jl) / MAX(v_s(ji,jj,jl),zeps) * zindsn + rtt
1213	ELSE
1214	ht_i(ji,jj,jl) = 0.0
1215	t_su(ji,jj,jl) = rtt
1216	! t_s(ji,jj,1,jl) = rtt
1217	ENDIF
1218	END DO ! ji
1219	END DO ! jj
1220	END DO ! jl
1221
1222	END SUBROUTINE lim_itd_shiftice
1223
1224	!----------------------------------------------------------------------------------------
1225	!----------------------------------------------------------------------------------------
1226
1227	SUBROUTINE lim_itd_th_reb(klbnd, kubnd, ntyp)
1228	!!------------------------------------------------------------------
1229	!! * ROUTINE lim_itd_th_reb *
1230	!! ** Purpose : rebin - rebins thicknesses into defined categories
1231	!!
1232	!! ** Method :
1233	!!
1234	!! ** Arguments :
1235	!!
1236	!! ** Inputs / Ouputs : (global commons)
1237	!!
1238	!! ** External :
1239	!!
1240	!! ** References :
1241	!!
1242	!! ** History : (2005) Translation from CICE
1243	!! (2006) Adaptation to include salt, age and types
1244	!! (2007) Mass conservation checked
1245	!!
1246	!! authors: William H. Lipscomb, LANL, Elizabeth C. Hunke, LANL
1247	!! (01-2006) Martin Vancoppenolle (adaptation)
1248	!!
1249	!! Quel etre encore que celui-ci! Le Jugement Dernier sera la
1250	!! avant qu'il vous fasse jamais une avance sur votre mois,
1251	!! Seigneur! Tu peux supplier, te mettre en quatre,
1252	!! meme si tu es dans la misere, il ne te donnera rien,
1253	!! le vieux demon! Et quant on pense que, chez lui,
1254	!! sa cuisiniere lui donne des gifles! Je ne vois pas l'interet
1255	!! qu'il y a a travailler dans un ministere. Cela ne rapporte
1256	!! absolument rien.
1257	!!
1258	!! This routine was largely inspired from CICE
1259	!! (W.H. Lipscomb, E. Hunke, C. M. Bitz)
1260	!! the sea ice model of LANL, Los Alamos, USA.
1261	!!
1262	!!------------------------------------------------------------------
1263	!! * Arguments
1264	INTEGER , INTENT (in) :: &
1265	klbnd , & ! Start thickness category index point
1266	kubnd , & ! End point on which the the computation is applied
1267	ntyp ! number of the ice type involved in the rebinning process
1268
1269	INTEGER :: &
1270	ji,jj, & ! horizontal indices
1271	jl ! category index
1272
1273	LOGICAL :: & !:
1274	zshiftflag ! = .true. if ice must be shifted
1275
1276	INTEGER, DIMENSION(jpi,jpj,jpl) :: &
1277	zdonor ! donor category index
1278
1279	REAL(wp), DIMENSION(jpi, jpj, jpl) :: &
1280	zdaice , & ! ice area transferred
1281	zdvice ! ice volume transferred
1282
1283	REAL(wp) :: & ! constant values
1284	zeps = 1.0e-10, &
1285	epsi10 = 1.0e-10, &
1286	zindb
1287
1288	REAL (wp), DIMENSION(jpi,jpj) :: & !
1289	vt_i_init, vt_i_final, & ! ice volume summed over categories
1290	vt_s_init, vt_s_final ! snow volume summed over categories
1291
1292	CHARACTER (len = 15) :: fieldid
1293
1294	!!-- End of declarations
1295	!------------------------------------------------------------------------------
1296	! WRITE(numout,*) ' lim_itd_th_reb '
1297	! WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
1298	! WRITE(numout,*) ' Ice Type no ', ntyp
1299	! WRITE(numout,*) ' bounds of categories ', klbnd, kubnd
1300
1301	! ! conservation check
1302	! CALL lim_column_sum (jpl, v_i, vt_i_init)
1303	! CALL lim_column_sum (jpl, v_s, vt_s_init)
1304
1305	!
1306	!------------------------------------------------------------------------------
1307	! 1) Compute ice thickness.
1308	!------------------------------------------------------------------------------
1309	! nothing to do
1310	DO jl = klbnd, kubnd
1311	DO jj = 1, jpj
1312	DO ji = 1, jpi
1313	IF (a_i(ji,jj,jl) .GT. zeps) THEN
1314	ht_i(ji,jj,jl) = v_i(ji,jj,jl) / a_i(ji,jj,jl)
1315	ELSE
1316	ht_i(ji,jj,jl) = 0.0
1317	ENDIF
1318	END DO ! i
1319	END DO ! j
1320	END DO ! n
1321
1322	!------------------------------------------------------------------------------
1323	! 2) Make sure thickness of cat klbnd is at least hi_max_typ(klbnd)
1324	!------------------------------------------------------------------------------
1325	DO jj = 1, jpj
1326	DO ji = 1, jpi
1327
1328	IF (a_i(ji,jj,klbnd) > zeps) THEN
1329	IF (ht_i(ji,jj,klbnd) .LE. hi_max_typ(0,ntyp) .AND. hi_max_typ(0,ntyp) .GT. 0.0 ) THEN
1330	a_i(ji,jj,klbnd) = v_i(ji,jj,klbnd) / hi_max_typ(0,ntyp)
1331	ht_i(ji,jj,klbnd) = hi_max_typ(0,ntyp)
1332	ENDIF
1333	ENDIF
1334	END DO ! i
1335	END DO ! j
1336
1337	!------------------------------------------------------------------------------
1338	! 3) If a category thickness is not in bounds, shift the
1339	! entire area, volume, and energy to the neighboring category
1340	!------------------------------------------------------------------------------
1341	!-------------------------
1342	! Initialize shift arrays
1343	!-------------------------
1344
1345	DO jl = klbnd, kubnd
1346	DO jj = 1, jpj
1347	DO ji = 1, jpi
1348	zdonor(ji,jj,jl) = 0
1349	zdaice(ji,jj,jl) = 0.0
1350	zdvice(ji,jj,jl) = 0.0
1351	END DO
1352	END DO
1353	END DO
1354
1355	!-------------------------
1356	! Move thin categories up
1357	!-------------------------
1358
1359	DO jl = klbnd, kubnd - 1 ! loop over category boundaries
1360
1361	!---------------------------------------
1362	! identify thicknesses that are too big
1363	!---------------------------------------
1364	zshiftflag = .false.
1365
1366	DO jj = 1, jpj
1367	DO ji = 1, jpi
1368	IF (a_i(ji,jj,jl) .GT. zeps .AND. ht_i(ji,jj,jl) .GT. hi_max(jl) ) THEN
1369	zshiftflag = .true.
1370	zdonor(ji,jj,jl) = jl
1371	zdaice(ji,jj,jl) = a_i(ji,jj,jl)
1372	zdvice(ji,jj,jl) = v_i(ji,jj,jl)
1373	ENDIF
1374	END DO ! ji
1375	END DO ! jj
1376
1377	IF (zshiftflag) THEN
1378
1379	!------------------------------
1380	! Shift ice between categories
1381	!------------------------------
1382	CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice)
1383
1384	!------------------------
1385	! Reset shift parameters
1386	!------------------------
1387	DO jj = 1, jpj
1388	DO ji = 1, jpi
1389	zdonor(ji,jj,jl) = 0
1390	zdaice(ji,jj,jl) = 0.0
1391	zdvice(ji,jj,jl) = 0.0
1392	END DO
1393	END DO
1394
1395	ENDIF ! zshiftflag
1396
1397	END DO ! jl
1398
1399	!----------------------------
1400	! Move thick categories down
1401	!----------------------------
1402
1403	DO jl = kubnd - 1, 1, -1 ! loop over category boundaries
1404
1405	!-----------------------------------------
1406	! Identify thicknesses that are too small
1407	!-----------------------------------------
1408	zshiftflag = .false.
1409
1410	DO jj = 1, jpj
1411	DO ji = 1, jpi
1412	IF (a_i(ji,jj,jl+1) .GT. zeps .AND. &
1413	ht_i(ji,jj,jl+1) .LE. hi_max(jl)) THEN
1414
1415	zshiftflag = .true.
1416	zdonor(ji,jj,jl) = jl + 1
1417	zdaice(ji,jj,jl) = a_i(ji,jj,jl+1)
1418	zdvice(ji,jj,jl) = v_i(ji,jj,jl+1)
1419	ENDIF
1420	END DO ! ji
1421	END DO ! jj
1422
1423	IF (zshiftflag) THEN
1424
1425	!------------------------------
1426	! Shift ice between categories
1427	!------------------------------
1428	CALL lim_itd_shiftice (klbnd, kubnd, zdonor, zdaice, zdvice)
1429
1430	!------------------------
1431	! Reset shift parameters
1432	!------------------------
1433	DO jj = 1, jpj
1434	DO ji = 1, jpi
1435	zdonor(ji,jj,jl) = 0
1436	zdaice(ji,jj,jl) = 0.0
1437	zdvice(ji,jj,jl) = 0.0
1438	END DO
1439	END DO
1440
1441	ENDIF ! zshiftflag
1442
1443	END DO ! jl
1444
1445	!------------------------------------------------------------------------------
1446	! 4) Conservation check
1447	!------------------------------------------------------------------------------
1448
1449	! CALL lim_column_sum (jpl, v_i, vt_i_final)
1450	! fieldid = ' v_i : limitd_reb '
1451	! CALL lim_cons_check (vt_i_init, vt_i_final, 1.0e-6, fieldid)
1452
1453	! CALL lim_column_sum (jpl, v_s, vt_s_final)
1454	! fieldid = ' v_s : limitd_reb '
1455	! CALL lim_cons_check (vt_s_init, vt_s_final, 1.0e-6, fieldid)
1456
1457	END SUBROUTINE lim_itd_th_reb
1458
1459	#else
1460	!!======================================================================
1461	!! * MODULE limitd_th *
1462	!! no sea ice model
1463	!!======================================================================
1464	CONTAINS
1465	SUBROUTINE lim_itd_th ! Empty routines
1466	END SUBROUTINE lim_itd_th
1467	SUBROUTINE lim_itd_th_ini
1468	END SUBROUTINE lim_itd_th_ini
1469	SUBROUTINE lim_itd_th_rem
1470	END SUBROUTINE lim_itd_th_rem
1471	SUBROUTINE lim_itd_fitline
1472	END SUBROUTINE lim_itd_fitline
1473	SUBROUTINE lim_itd_shiftice
1474	END SUBROUTINE lim_itd_shiftice
1475	SUBROUTINE lim_itd_th_reb
1476	END SUBROUTINE lim_itd_th_reb
1477	#endif
1478	END MODULE limitd_th

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

Download in other formats: