Context Navigation

source: lim1d_ws/trunk/SOURCES/source_3.20/ice_cons.f @ 2

Last change on this file since 2 was 2, checked in by vancop, 8 years ago
initial import /Users/ioulianikolskaia/Boulot/CODES/LIM1D/ARCHIVE/TMP/LIM1D_v3.20/
File size: 14.5 KB

Line
1	SUBROUTINE ice_th_glohec(eti,ets,etilayer,kideb,kiut,jl,
2	& nlay_s, nlay_i)
3
4	!!-----------------------------------------------------------------------
5	!! * ROUTINE ice_th_glohec *
6	!!
7	!! ** Purpose : Compute total heat content for each category
8	!! Works with 1d vectors only
9	!!
10	!! history :
11	!! 9.9 ! 07-04 (M.Vancoppenolle) original code
12	!!-----------------------------------------------------------------------
13	!! * Local variables
14
15	INCLUDE 'type.com'
16	INCLUDE 'para.com'
17	INCLUDE 'const.com'
18	INCLUDE 'ice.com'
19	INCLUDE 'thermo.com'
20
21	INTEGER, INTENT(in) ::
22	& kideb, kiut, ! bounds for the spatial loop
23	& jl ! category number
24
25	REAL(8), DIMENSION (nbpt,jpl), INTENT(out) ::
26	& eti, ets ! vertically-summed heat content for ice /snow
27
28	REAL(8), DIMENSION (nbpt,maxnlay), INTENT(out) ::
29	& etilayer ! heat content for ice layers
30
31	REAL(8) ::
32	& zdes, ! snow heat content increment (dummy)
33	& zeps ! very small value (1.e-10)
34
35	INTEGER ::
36	& ji,jj,jk ! loop indices
37
38	!!-----------------------------------------------------------------------
39	eti(:,:) = 0.0
40	ets(:,:) = 0.0
41	zeps = 1.0e-10
42
43	! total q over all layers, ice [J.m-2]
44	DO jk = 1, nlay_i
45	DO ji = kideb, kiut
46	etilayer(ji,jk) = q_i_b(ji,jk)
47	& * ht_i_b(ji) / nlay_i
48	eti(ji,jl) = eti(ji,jl) + etilayer(ji,jk)
49	END DO
50	END DO
51
52	! total q over all layers, snow [J.m-2]
53	DO ji = kideb, kiut
54	zdes = q_s_b(ji,1) * ht_s_b(ji) / nlay_s
55	ets(ji,jl) = ets(ji,jl) + zdes
56	END DO
57
58	RETURN
59	!------------------------------------------------------------------------------
60	END SUBROUTINE ice_th_glohec
61
62	!==============================================================================
63
64	SUBROUTINE ice_th_enmelt(kideb, kiut, nlay_s, nlay_i)
65	!!-----------------------------------------------------------------------
66	!! * ROUTINE ice_th_enmelt *
67	!!
68	!! ** Purpose : Computes sea ice energy of melting q_i (J.m-3)
69	!!
70	!! ** Method : Formula (Bitz and Lipscomb, 1999)
71	!!
72	!! history : Martin Vancoppenolle, May 2007
73	!!-------------------------------------------------------------------
74
75	INCLUDE 'type.com'
76	INCLUDE 'para.com'
77	INCLUDE 'const.com'
78	INCLUDE 'ice.com'
79	INCLUDE 'thermo.com'
80
81	REAL(8) :: !: goes to trash
82	& ztmelts , !: sea ice freezing point in K
83	& zeps
84
85	INTEGER ::
86	& ji, !: spatial loop index
87	& jk !: vertical index
88
89	LOGICAL ::
90	& ln_write = .FALSE.
91
92	!!-------------------------------------------------------------------
93	zeps = 1.0e-10
94	IF ( ln_write ) THEN
95	WRITE(numout,) ' * ice_th_enmelt : '
96	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~ '
97	WRITE(numout,*)
98	WRITE(numout,*) ' kideb : ', kideb
99	WRITE(numout,*) ' kiut : ', kiut
100	WRITE(numout,*) ' nlay_s: ', nlay_s
101	WRITE(numout,*) ' nlay_i: ', nlay_i
102	ENDIF
103
104	! Sea ice energy of melting
105	DO jk = 1, nlay_i
106	DO ji = kideb, kiut
107	ztmelts = - tmut * s_i_b(ji,jk) + tpw
108	q_i_b(ji,jk) = rhog * ( cpg * ( ztmelts - t_i_b(ji,jk) )
109	& + lfus*( 1.0 - (ztmelts-tpw) /
110	& MIN((t_i_b(ji,jk)-tpw),-zeps) )
111	& - cpw * ( ztmelts-tpw ) )
112	END DO !ji
113	END DO !jk
114
115	! Snow energy of melting
116	DO jk = 1, nlay_s
117	DO ji = kideb,kiut
118	q_s_b(ji,jk) = rhon * ( cpg * ( tpw - t_s_b(ji,jk) ) +
119	& lfus )
120	END DO !ji
121	END DO !jk
122
123	WRITE(numout,*)
124
125	RETURN
126	!------------------------------------------------------------------------------
127	END SUBROUTINE ice_th_enmelt
128
129	!==============================================================================
130
131	SUBROUTINE ice_th_con_dif(kideb,kiut,nlay_s,nlay_i,jl)
132	!!-----------------------------------------------------------------------
133	!! * ROUTINE ice_th_con_dif *
134	!!
135	!! ** Purpose : Test energy conservation after heat diffusion
136	!!
137	!! history :
138	!! 9.9 ! 07-04 (M.Vancoppenolle) original code
139	!!-------------------------------------------------------------------
140	!! * Local variables
141
142	INCLUDE 'type.com'
143	INCLUDE 'para.com'
144	INCLUDE 'const.com'
145	INCLUDE 'ice.com'
146	INCLUDE 'thermo.com'
147
148	REAL(8) :: !: ! goes to trash
149	& meance, !: mean conservation error
150	& max_cons_err, !: maximum tolerated conservation error
151	& max_surf_err !: maximum tolerated surface error
152
153	INTEGER ::
154	& numce !: number of points for which conservation
155	! is violated
156	INTEGER ::
157	& ji,jj,jk, !: loop indices
158	& zji, zjj
159
160	LOGICAL ::
161	& l_write
162	!!---------------------------------------------------------------------
163
164	max_cons_err = 0.1
165	max_surf_err = 0.001
166	l_write = .TRUE.
167
168	IF ( l_write ) THEN
169	WRITE(numout,) ' * ice_th_con_dif : '
170	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~~ '
171	ENDIF
172
173	!--------------------------
174	! 1) Increment of energy
175	!--------------------------
176	! vertically integrated
177	DO ji = kideb, kiut
178	dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl)
179	& + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
180	END DO
181
182	! layer by layer
183	DO ji = kideb, kiut
184	DO jk = 1, nlay_i
185	dq_i_layer(ji,jk) = q_i_layer_fin(ji,jk) -
186	& q_i_layer_in(ji,jk)
187	END DO
188	END DO
189
190	!-------------------------------------------
191	! 2) Atmospheric heat flux, ice heat budget
192	!-------------------------------------------
193
194	DO ji = kideb, kiut
195	fatm(ji,jl) = ab(ji)*fsolgb(ji) + fratsb(ji)
196	& + fcsb(ji) + fleb(ji)
197	sum_fluxq(ji,jl) = fc_su(ji) - fc_bo_i(ji)
198	& + fsolgb(ji)*(1-ab(ji)) - ftroce
199	END DO
200
201	IF ( l_write ) THEN
202	DO ji = kideb, kiut
203	WRITE(numout,*) ' fc_su : ', fc_su(ji)
204	WRITE(numout,*) ' fc_bo_i : ', fc_bo_i(ji)
205	WRITE(numout,) ' fsolio : ', fsolgb(ji)*(1.-ab(ji))
206	WRITE(numout,*) ' ftroce : ', ftroce
207	END DO
208	ENDIF
209
210	!-----------------------
211	! 3) Conservation error
212	!-----------------------
213	DO ji = kideb, kiut
214	cons_error(ji,jl) = ABS( dq_i(ji,jl) / ddtb +
215	& sum_fluxq(ji,jl) )
216	END DO
217
218	numce = 0
219	meance = 0.0
220	DO ji = kideb, kiut
221	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
222	numce = numce + 1
223	meance = meance + cons_error(ji,jl)
224	ENDIF
225	IF (numce .GT. 0 ) meance = meance / numce
226
227	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
228	WRITE(numout,*) ' Diffusion of heat - large cons. error '
229	WRITE(numout,*) ' Maximum tolerated conservation error : ',
230	& max_cons_err
231	WRITE(numout,*) ' Thickness category : ', jl
232	WRITE(numout,*) ' DIF Mean conservation error ',
233	& meance, numit
234	ENDIF
235
236	ENDDO
237
238	!---------------------------------------------------------
239	! 4) Surface error due to imbalance between Fatm and Fcsu
240	!---------------------------------------------------------
241	numce = 0.0
242	meance = 0.0
243
244	DO ji = kideb, kiut
245	surf_error(ji,jl) = ABS ( fatm(ji,jl) - fc_su(ji) )
246	IF ( ( t_su_b(ji) .LT. tpw ) .AND. ( surf_error(ji,jl) .GT.
247	& max_surf_err ) ) THEN
248	numce = numce + 1
249	meance = meance + surf_error(ji,jl)
250	ENDIF
251	IF (numce .GT. 0 ) meance = meance / numce
252
253	IF ( ( t_su_b(ji) .LT. tpw ) .AND. ( surf_error(ji,jl) .GT.
254	& max_surf_err ) ) THEN
255	WRITE(numout,*) ' Diffusion of heat - large surface error '
256	WRITE(numout,*) ' Surf_error : ', surf_error(1,jl),
257	& ' (W/m2) '
258	WRITE(numout,*) ' Maximum tolerated surface error : ',
259	& max_surf_err
260	WRITE(numout,*) ' Thickness category : ', jl
261	WRITE(numout,*) ' t_su_b: ', t_su_b(ji)
262	WRITE(numout,*) ' t_s_b : ', ( t_s_b(ji,layer),
263	& layer = 1, nlay_s )
264	WRITE(numout,*) ' t_i_b : ', ( t_i_b(ji,layer),
265	& layer = 1, nlay_i )
266	ENDIF
267
268	ENDDO
269
270	!-------------------
271	! 5) Layer by layer
272	!-------------------
273	IF ( l_write ) THEN
274	WRITE(numout,*)
275	WRITE(numout,*) ' Layer by layer ... '
276
277	DO ji = kideb, kiut
278	WRITE(numout,*) ' T_su : ', t_su_b(ji)
279	WRITE(numout,) ' Snow * h_s = ',
280	& ht_s_b(ji)
281
282	IF ( ht_s_b(ji) .GT. 0.0 ) THEN
283	zdqs = ( qt_s_fin(ji,jl) - qt_s_in(ji,jl) ) / ddtb
284	WRITE(numout,*) ' Conservation error (W/m2): ',
285	& ABS( zdqs + fc_s(ji,0) -
286	& fc_s(ji,1) + radab_s(1) )
287	ENDIF
288	END DO
289
290	DO ji = kideb, kiut
291	WRITE(numout,) ' Ice * h_i = ', ht_i_b(ji)
292	DO jk = 1, nlay_i
293	WRITE(numout,*) ' LAYER : ', jk
294	WRITE(numout,*) ' Conservation error (W/m2): ',
295	& ABS( dq_i_layer(ji,jk) / ddtb + fc_i(ji,jk-1) -
296	& fc_i(ji,jk) + radab_i(jk) ) !radab_phy_i(jk) + radab_alg_i(jk) )
297	END DO
298	END DO
299
300	ENDIF
301	!
302	!------------------------------------------------------------------------------
303	END SUBROUTINE ice_th_con_dif
304	!==============================================================================
305
306	SUBROUTINE ice_th_con_dh(kideb,kiut,nlay_s,nlay_i,jl)
307	!!-----------------------------------------------------------------------
308	!! * ROUTINE ice_th_con_dh *
309	!!
310	!! ** Purpose : Test energy conservation after ice growth and melt
311	!!
312	!! history :
313	!! 9.9 ! 07-04 (M.Vancoppenolle) original code
314	!!-------------------------------------------------------------------
315	!! * Local variables
316
317	INCLUDE 'type.com'
318	INCLUDE 'para.com'
319	INCLUDE 'const.com'
320	INCLUDE 'ice.com'
321	INCLUDE 'thermo.com'
322
323	REAL(8) :: !: ! goes to trash
324	& meance, !: mean conservation error
325	& max_cons_err, !: maximum tolerated conservation error
326	& max_surf_err !: maximum tolerated surface error
327
328	INTEGER ::
329	& numce !: number of points for which conservation
330	! is violated
331	INTEGER ::
332	& ji,jj,jk, !: loop indices
333	& zji, zjj
334	LOGICAL ::
335	& l_write
336
337	!!---------------------------------------------------------------------
338
339	l_write = .TRUE.
340	max_cons_err = 0.1
341	max_surf_err = 0.001
342
343	WRITE(numout,) ' * ice_th_con_dh : '
344	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~~~ '
345
346	!------------------------
347	! 1) Increment of energy
348	!------------------------
349	DO ji = kideb, kiut ! vertically integrated
350	dq_i(ji,jl) = qt_i_fin(ji,jl) - qt_i_in(ji,jl)
351	& + qt_s_fin(ji,jl) - qt_s_in(ji,jl)
352	END DO
353
354	DO ji = kideb, kiut ! layer by layer
355	DO jk = 1, nlay_i
356	dq_i_layer(ji,jk) = q_i_layer_fin(ji,jk) -
357	& q_i_layer_in(ji,jk)
358	END DO
359	END DO
360
361	!-------------------------------------------
362	! 2) Atmospheric heat flux, ice heat budget
363	!-------------------------------------------
364	DO ji = kideb, kiut
365	fatm(ji,jl) = fsolgb(ji) + fratsb(ji)
366	& + fcsb(ji) + fleb(ji)
367	sum_fluxq(ji,jl) = fatm(ji,jl) + fbbqb(ji)
368	& - ftroce + fsnic + fprec
369	END DO
370
371	IF ( l_write ) THEN
372	DO ji = kideb, kiut
373	WRITE(numout,*) ' fsolgb : ' , fsolgb(ji)
374	WRITE(numout,*) ' fratsb : ' , fratsb(ji)
375	WRITE(numout,*) ' fcsb : ' , fcsb(ji)
376	WRITE(numout,*) ' fleb : ' , fleb(ji)
377	WRITE(numout,*)
378	WRITE(numout,*) ' fatm : ' , fatm(ji,jl)
379	WRITE(numout,*) ' fbbqb : ' , fbbqb(ji)
380	WRITE(numout,*) ' ftroce : ' , ftroce
381	WRITE(numout,*) ' fprec : ' , fprec
382	WRITE(numout,*) ' fsnic : ' , fsnic
383	WRITE(numout,*)
384	WRITE(numout,*) ' sum_fluxq : ', sum_fluxq(ji,jl)
385	WRITE(numout,*) ' dq_i : ', dq_i(ji,jl) / ddtb
386	END DO
387	ENDIF
388
389	!-----------------------
390	! 3) Conservation error
391	!-----------------------
392	IF ( l_write ) WRITE(numout,*) ' kideb, kiut, jl : ',
393	& kideb, kiut, jl
394
395	DO ji = kideb, kiut
396	cons_error(ji,jl) = ABS( dq_i(ji,jl) / ddtb +
397	& sum_fluxq(ji,jl) )
398	END DO
399
400	numce = 0
401	meance = 0.0
402	DO ji = kideb, kiut
403	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
404	numce = numce + 1
405	meance = meance + cons_error(ji,jl)
406	ENDIF
407	IF (numce .GT. 0 ) meance = meance / numce
408
409	IF ( cons_error(ji,jl) .GT. max_cons_err ) THEN
410	WRITE(numout,*) ' Growth and melt - large cons. error '
411	WRITE(numout,*) ' Maximum tolerated conservation error : ',
412	& max_cons_err
413	WRITE(numout,*) ' DH Mean conservation error ',
414	& meance, numit
415	ENDIF
416	ENDDO
417
418	!------------------------------------------------------------------------------
419	END SUBROUTINE ice_th_con_dh

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

Download in other formats: