Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

icectl.F90 @ 8879

Last change on this file since 8879 was 8879, checked in by frrh, 6 years ago
Merge in http://fcm3/projects/NEMO.xm/log/branches/UKMO/dev_r8183_ICEMODEL_svn_removed revisions 8738:8847 inclusive.
File size: 34.7 KB

Line
1	MODULE icectl
2	!!======================================================================
3	!! * MODULE icectl *
4	!! sea-ice : controls and prints
5	!!======================================================================
6	!! History : 3.5 ! 2015-01 (M. Vancoppenolle) Original code
7	!! 3.7 ! 2016-10 (C. Rousset) Add routine ice_prt3D
8	!!----------------------------------------------------------------------
9	#if defined key_lim3
10	!!----------------------------------------------------------------------
11	!! 'key_lim3' ESIM sea-ice model
12	!!----------------------------------------------------------------------
13	!! ice_cons_hsm : conservation tests on heat, salt and mass
14	!! ice_cons_final : conservation tests on heat, salt and mass at end of time step
15	!! ice_ctl : control prints in case of crash
16	!! ice_prt : control prints at a given grid point
17	!! ice_prt3D : control prints of ice arrays
18	!!----------------------------------------------------------------------
19	USE phycst ! physical constants
20	USE oce ! ocean dynamics and tracers
21	USE dom_oce ! ocean space and time domain
22	USE ice ! sea-ice: variables
23	USE ice1D ! sea-ice: thermodynamics variables
24	USE sbc_oce ! Surface boundary condition: ocean fields
25	USE sbc_ice ! Surface boundary condition: ice fields
26	!
27	USE in_out_manager ! I/O manager
28	USE lib_mpp ! MPP library
29	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
30	USE timing ! Timing
31	USE prtctl ! Print control
32
33	IMPLICIT NONE
34	PRIVATE
35
36	PUBLIC ice_cons_hsm
37	PUBLIC ice_cons_final
38	PUBLIC ice_ctl
39	PUBLIC ice_prt
40	PUBLIC ice_prt3D
41
42	!! * Substitutions
43	# include "vectopt_loop_substitute.h90"
44	!!----------------------------------------------------------------------
45	!! NEMO/ICE 4.0 , NEMO Consortium (2017)
46	!! $Id: icectl.F90 5043 2015-01-28 16:44:18Z clem $
47	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
48	!!----------------------------------------------------------------------
49	CONTAINS
50
51	SUBROUTINE ice_cons_hsm( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft )
52	!!-------------------------------------------------------------------
53	!! * ROUTINE ice_cons_hsm *
54	!!
55	!! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
56	!! + test if ice concentration and volume are > 0
57	!!
58	!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
59	!! It prints in ocean.output if there is a violation of conservation at each time-step
60	!! The thresholds (zv_sill, zs_sill, zt_sill) which determine violations are set to
61	!! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
62	!! For salt and heat thresholds, ice is considered to have a salinity of 10
63	!! and a heat content of 3e5 J/kg (=latent heat of fusion)
64	!!-------------------------------------------------------------------
65	INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end
66	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
67	REAL(wp) , INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
68	!!
69	REAL(wp) :: zv, zs, zt, zfs, zfv, zft
70	REAL(wp) :: zvmin, zamin, zamax
71	REAL(wp) :: zvtrp, zetrp
72	REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill
73	REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt
74	!!-------------------------------------------------------------------
75	!
76	IF( icount == 0 ) THEN
77	! ! water flux
78	pdiag_fv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + &
79	& wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + &
80	& wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + &
81	& wfx_ice_sub(:,:) + wfx_spr(:,:) &
82	& ) * e1e2t(:,:) ) * zconv
83	!
84	! ! salt flux
85	pdiag_fs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + &
86	& sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) &
87	& ) * e1e2t(:,:) ) * zconv
88	!
89	! ! heat flux
90	pdiag_ft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &
91	& - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &
92	& ) * e1e2t(:,:) ) * zconv
93
94	pdiag_v = glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t * zconv )
95
96	pdiag_s = glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t * zconv )
97
98	pdiag_t = glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) &
99	& + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv
100
101	ELSEIF( icount == 1 ) THEN
102
103	! water flux
104	zfv = glob_sum( -( wfx_bog(:,:) + wfx_bom(:,:) + wfx_sum(:,:) + wfx_sni(:,:) + &
105	& wfx_opw(:,:) + wfx_res(:,:) + wfx_dyn(:,:) + wfx_lam(:,:) + wfx_pnd(:,:) + &
106	& wfx_snw_sni(:,:) + wfx_snw_sum(:,:) + wfx_snw_dyn(:,:) + wfx_snw_sub(:,:) + &
107	& wfx_ice_sub(:,:) + wfx_spr(:,:) &
108	& ) * e1e2t(:,:) ) * zconv - pdiag_fv
109
110	! salt flux
111	zfs = glob_sum( ( sfx_bri(:,:) + sfx_bog(:,:) + sfx_bom(:,:) + sfx_sum(:,:) + sfx_sni(:,:) + &
112	& sfx_opw(:,:) + sfx_res(:,:) + sfx_dyn(:,:) + sfx_sub(:,:) + sfx_lam(:,:) &
113	& ) * e1e2t(:,:) ) * zconv - pdiag_fs
114
115	! heat flux
116	zft = glob_sum( ( hfx_sum(:,:) + hfx_bom(:,:) + hfx_bog(:,:) + hfx_dif(:,:) + hfx_opw(:,:) + hfx_snw(:,:) &
117	& - hfx_thd(:,:) - hfx_dyn(:,:) - hfx_res(:,:) - hfx_sub(:,:) - hfx_spr(:,:) &
118	& ) * e1e2t(:,:) ) * zconv - pdiag_ft
119
120	! outputs
121	zv = ( ( glob_sum( SUM( v_i * rhoic + v_s * rhosn, dim=3 ) * e1e2t ) * zconv &
122	& - pdiag_v ) * r1_rdtice - zfv ) * rday
123
124	zs = ( ( glob_sum( SUM( sv_i * rhoic , dim=3 ) * e1e2t ) * zconv &
125	& - pdiag_s ) * r1_rdtice + zfs ) * rday
126
127	zt = ( glob_sum( ( SUM( SUM( e_i(:,:,1:nlay_i,:), dim=4 ), dim=3 ) &
128	& + SUM( SUM( e_s(:,:,1:nlay_s,:), dim=4 ), dim=3 ) ) * e1e2t ) * zconv &
129	& - pdiag_t ) * r1_rdtice + zft
130
131	! zvtrp and zetrp must be close to 0 if the advection scheme is conservative
132	zvtrp = glob_sum( ( diag_trp_vi * rhoic + diag_trp_vs * rhosn ) * e1e2t ) * zconv * rday
133	zetrp = glob_sum( ( diag_trp_ei + diag_trp_es ) * e1e2t ) * zconv
134
135	zvmin = glob_min( v_i )
136	zamax = glob_max( SUM( a_i, dim=3 ) )
137	zamin = glob_min( a_i )
138
139	! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
140	zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2
141	zv_sill = zarea * 2.5e-5
142	zs_sill = zarea * 25.e-5
143	zt_sill = zarea * 10.e-5
144
145	IF(lwp) THEN
146	IF ( ABS( zv ) > zv_sill ) WRITE(numout,*) 'violation volume [Mt/day] (',cd_routine,') = ',zv
147	IF ( ABS( zs ) > zs_sill ) WRITE(numout,) 'violation saline [psuMt/day] (',cd_routine,') = ',zs
148	IF ( ABS( zt ) > zt_sill ) WRITE(numout,*) 'violation enthalpy [GW] (',cd_routine,') = ',zt
149	IF ( ABS(zvtrp ) > zv_sill .AND. cd_routine == 'icedyn_adv' ) THEN
150	WRITE(numout,*) 'violation vtrp [Mt/day] (',cd_routine,') = ',zvtrp
151	WRITE(numout,*) 'violation etrp [GW] (',cd_routine,') = ',zetrp
152	ENDIF
153	IF ( zvmin < -epsi10 ) WRITE(numout,*) 'violation v_i<0 [m] (',cd_routine,') = ',zvmin
154	IF ( zamax > MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &
155	& WRITE(numout,*) 'violation a_i>amax (',cd_routine,') = ',zamax
156	IF ( zamin < -epsi10 ) WRITE(numout,*) 'violation a_i<0 (',cd_routine,') = ',zamin
157	ENDIF
158	!
159	ENDIF
160
161	END SUBROUTINE ice_cons_hsm
162
163
164	SUBROUTINE ice_cons_final( cd_routine )
165	!!-------------------------------------------------------------------
166	!! * ROUTINE ice_cons_final *
167	!!
168	!! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
169	!!
170	!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
171	!! It prints in ocean.output if there is a violation of conservation at each time-step
172	!! The thresholds (zv_sill, zs_sill, zt_sill) which determine the violation are set to
173	!! a minimum of 1 mm of ice (over the ice area) that is lost/gained spuriously during 100 years.
174	!! For salt and heat thresholds, ice is considered to have a salinity of 10
175	!! and a heat content of 3e5 J/kg (=latent heat of fusion)
176	!!-------------------------------------------------------------------
177	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
178	REAL(wp) :: zhfx, zsfx, zvfx
179	REAL(wp) :: zarea, zv_sill, zs_sill, zt_sill
180	REAL(wp), PARAMETER :: zconv = 1.e-9 ! convert W to GW and kg to Mt
181	!!-------------------------------------------------------------------
182
183	! water flux
184	zvfx = glob_sum( ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + diag_vice + diag_vsnw ) * e1e2t ) * zconv * rday
185
186	! salt flux
187	zsfx = glob_sum( ( sfx + diag_sice ) * e1e2t ) * zconv * rday
188
189	! heat flux
190	zhfx = glob_sum( ( hfx_in - hfx_out - diag_heat - diag_trp_ei - diag_trp_es &
191	! & - SUM( qevap_ice * a_i_b, dim=3 ) & !!clem: I think this line must be commented (but need check)
192	& ) * e1e2t ) * zconv
193
194	! set threshold values and calculate the ice area (+epsi10 to set a threshold > 0 when there is no ice)
195	zarea = glob_sum( SUM( a_i + epsi10, dim=3 ) * e1e2t ) * zconv ! in 1.e9 m2
196	zv_sill = zarea * 2.5e-5
197	zs_sill = zarea * 25.e-5
198	zt_sill = zarea * 10.e-5
199
200	IF(lwp) THEN
201	IF( ABS( zvfx ) > zv_sill ) WRITE(numout,*) 'violation vfx [Mt/day] (',cd_routine,') = ',zvfx
202	IF( ABS( zsfx ) > zs_sill ) WRITE(numout,) 'violation sfx [psuMt/day] (',cd_routine,') = ',zsfx
203	IF( ABS( zhfx ) > zt_sill ) WRITE(numout,*) 'violation hfx [GW] (',cd_routine,') = ',zhfx
204	ENDIF
205	!
206	END SUBROUTINE ice_cons_final
207
208
209	SUBROUTINE ice_ctl( kt )
210	!!-------------------------------------------------------------------
211	!! * ROUTINE ice_ctl *
212	!!
213	!! ** Purpose : Alerts in case of model crash
214	!!-------------------------------------------------------------------
215	INTEGER, INTENT(in) :: kt ! ocean time step
216	INTEGER :: ji, jj, jk, jl ! dummy loop indices
217	INTEGER :: inb_altests ! number of alert tests (max 20)
218	INTEGER :: ialert_id ! number of the current alert
219	REAL(wp) :: ztmelts ! ice layer melting point
220	CHARACTER (len=30), DIMENSION(20) :: cl_alname ! name of alert
221	INTEGER , DIMENSION(20) :: inb_alp ! number of alerts positive
222	!!-------------------------------------------------------------------
223
224	inb_altests = 10
225	inb_alp(:) = 0
226
227	! Alert if incompatible volume and concentration
228	ialert_id = 2 ! reference number of this alert
229	cl_alname(ialert_id) = ' Incompat vol and con ' ! name of the alert
230
231	DO jl = 1, jpl
232	DO jj = 1, jpj
233	DO ji = 1, jpi
234	IF( v_i(ji,jj,jl) /= 0._wp .AND. a_i(ji,jj,jl) == 0._wp ) THEN
235	!WRITE(numout,*) ' ALERTE 2 : Incompatible volume and concentration '
236	!WRITE(numout,*) ' at_i ', at_i(ji,jj)
237	!WRITE(numout,*) ' Point - category', ji, jj, jl
238	!WRITE(numout,) ' a_i ** a_i_b ', a_i (ji,jj,jl), a_i_b (ji,jj,jl)
239	!WRITE(numout,) ' v_i ** v_i_b ', v_i (ji,jj,jl), v_i_b (ji,jj,jl)
240	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
241	ENDIF
242	END DO
243	END DO
244	END DO
245
246	! Alerte if very thick ice
247	ialert_id = 3 ! reference number of this alert
248	cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert
249	jl = jpl
250	DO jj = 1, jpj
251	DO ji = 1, jpi
252	IF( h_i(ji,jj,jl) > 50._wp ) THEN
253	!CALL ice_prt( kt, ji, jj, 2, ' ALERTE 3 : Very thick ice ' )
254	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
255	ENDIF
256	END DO
257	END DO
258
259	! Alert if very fast ice
260	ialert_id = 4 ! reference number of this alert
261	cl_alname(ialert_id) = ' Very fast ice ' ! name of the alert
262	DO jj = 1, jpj
263	DO ji = 1, jpi
264	IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 1.5 .AND. &
265	& at_i(ji,jj) > 0._wp ) THEN
266	!CALL ice_prt( kt, ji, jj, 1, ' ALERTE 4 : Very fast ice ' )
267	!WRITE(numout,*) ' ice strength : ', strength(ji,jj)
268	!WRITE(numout,*) ' oceanic stress utau : ', utau(ji,jj)
269	!WRITE(numout,*) ' oceanic stress vtau : ', vtau(ji,jj)
270	!WRITE(numout,*) ' sea-ice stress utau_ice : ', utau_ice(ji,jj)
271	!WRITE(numout,*) ' sea-ice stress vtau_ice : ', vtau_ice(ji,jj)
272	!WRITE(numout,*) ' sst : ', sst_m(ji,jj)
273	!WRITE(numout,*) ' sss : ', sss_m(ji,jj)
274	!WRITE(numout,*)
275	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
276	ENDIF
277	END DO
278	END DO
279
280	! Alert if there is ice on continents
281	ialert_id = 6 ! reference number of this alert
282	cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert
283	DO jj = 1, jpj
284	DO ji = 1, jpi
285	IF( tmask(ji,jj,1) <= 0._wp .AND. at_i(ji,jj) > 0._wp ) THEN
286	!CALL ice_prt( kt, ji, jj, 1, ' ALERTE 6 : Ice on continents ' )
287	!WRITE(numout,*) ' masks s, u, v : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1)
288	!WRITE(numout,*) ' sst : ', sst_m(ji,jj)
289	!WRITE(numout,*) ' sss : ', sss_m(ji,jj)
290	!WRITE(numout,*) ' at_i(ji,jj) : ', at_i(ji,jj)
291	!WRITE(numout,*) ' v_ice(ji,jj) : ', v_ice(ji,jj)
292	!WRITE(numout,*) ' v_ice(ji,jj-1) : ', v_ice(ji,jj-1)
293	!WRITE(numout,*) ' u_ice(ji-1,jj) : ', u_ice(ji-1,jj)
294	!WRITE(numout,*) ' u_ice(ji,jj) : ', v_ice(ji,jj)
295	!
296	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
297	ENDIF
298	END DO
299	END DO
300
301	!
302	! ! Alert if very fresh ice
303	ialert_id = 7 ! reference number of this alert
304	cl_alname(ialert_id) = ' Very fresh ice ' ! name of the alert
305	DO jl = 1, jpl
306	DO jj = 1, jpj
307	DO ji = 1, jpi
308	IF( s_i(ji,jj,jl) < 0.1 .AND. a_i(ji,jj,jl) > 0._wp ) THEN
309	! CALL ice_prt(kt,ji,jj,1, ' ALERTE 7 : Very fresh ice ' )
310	! WRITE(numout,*) ' sst : ', sst_m(ji,jj)
311	! WRITE(numout,*) ' sss : ', sss_m(ji,jj)
312	! WRITE(numout,*)
313	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
314	ENDIF
315	END DO
316	END DO
317	END DO
318	!
319
320	! ! Alert if too old ice
321	ialert_id = 9 ! reference number of this alert
322	cl_alname(ialert_id) = ' Very old ice ' ! name of the alert
323	DO jl = 1, jpl
324	DO jj = 1, jpj
325	DO ji = 1, jpi
326	IF ( ( ( ABS( o_i(ji,jj,jl) ) > rdt_ice ) .OR. &
327	( ABS( o_i(ji,jj,jl) ) < 0._wp) ) .AND. &
328	( a_i(ji,jj,jl) > 0._wp ) ) THEN
329	!CALL ice_prt( kt, ji, jj, 1, ' ALERTE 9 : Wrong ice age ')
330	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
331	ENDIF
332	END DO
333	END DO
334	END DO
335
336	! Alert on salt flux
337	ialert_id = 5 ! reference number of this alert
338	cl_alname(ialert_id) = ' High salt flux ' ! name of the alert
339	DO jj = 1, jpj
340	DO ji = 1, jpi
341	IF( ABS( sfx (ji,jj) ) > 1.0e-2 ) THEN ! = 1 psu/day for 1m ocean depth
342	!CALL ice_prt( kt, ji, jj, 3, ' ALERTE 5 : High salt flux ' )
343	!DO jl = 1, jpl
344	!WRITE(numout,*) ' Category no: ', jl
345	!WRITE(numout,*) ' a_i : ', a_i (ji,jj,jl) , ' a_i_b : ', a_i_b (ji,jj,jl)
346	!WRITE(numout,*) ' v_i : ', v_i (ji,jj,jl) , ' v_i_b : ', v_i_b (ji,jj,jl)
347	!WRITE(numout,*) ' '
348	!END DO
349	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
350	ENDIF
351	END DO
352	END DO
353
354	! Alert if qns very big
355	ialert_id = 8 ! reference number of this alert
356	cl_alname(ialert_id) = ' fnsolar very big ' ! name of the alert
357	DO jj = 1, jpj
358	DO ji = 1, jpi
359	IF( ABS( qns(ji,jj) ) > 1500._wp .AND. at_i(ji,jj) > 0._wp ) THEN
360	!
361	!WRITE(numout,*) ' ALERTE 8 : Very high non-solar heat flux'
362	!WRITE(numout,*) ' ji, jj : ', ji, jj
363	!WRITE(numout,*) ' qns : ', qns(ji,jj)
364	!WRITE(numout,*) ' sst : ', sst_m(ji,jj)
365	!WRITE(numout,*) ' sss : ', sss_m(ji,jj)
366	!
367	!CALL ice_prt( kt, ji, jj, 2, ' ')
368	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
369	!
370	ENDIF
371	END DO
372	END DO
373	!+++++
374
375	! Alert if very warm ice
376	ialert_id = 10 ! reference number of this alert
377	cl_alname(ialert_id) = ' Very warm ice ' ! name of the alert
378	inb_alp(ialert_id) = 0
379	DO jl = 1, jpl
380	DO jk = 1, nlay_i
381	DO jj = 1, jpj
382	DO ji = 1, jpi
383	ztmelts = -tmut * sz_i(ji,jj,jk,jl) + rt0
384	IF( t_i(ji,jj,jk,jl) >= ztmelts .AND. v_i(ji,jj,jl) > 1.e-10 &
385	& .AND. a_i(ji,jj,jl) > 0._wp ) THEN
386	!WRITE(numout,*) ' ALERTE 10 : Very warm ice'
387	!WRITE(numout,*) ' ji, jj, jk, jl : ', ji, jj, jk, jl
388	!WRITE(numout,*) ' t_i : ', t_i(ji,jj,jk,jl)
389	!WRITE(numout,*) ' e_i : ', e_i(ji,jj,jk,jl)
390	!WRITE(numout,*) ' sz_i: ', sz_i(ji,jj,jk,jl)
391	!WRITE(numout,*) ' ztmelts : ', ztmelts
392	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
393	ENDIF
394	END DO
395	END DO
396	END DO
397	END DO
398
399	! sum of the alerts on all processors
400	IF( lk_mpp ) THEN
401	DO ialert_id = 1, inb_altests
402	CALL mpp_sum(inb_alp(ialert_id))
403	END DO
404	ENDIF
405
406	! print alerts
407	IF( lwp ) THEN
408	ialert_id = 1 ! reference number of this alert
409	cl_alname(ialert_id) = ' NO alerte 1 ' ! name of the alert
410	WRITE(numout,*) ' time step ',kt
411	WRITE(numout,*) ' All alerts at the end of ice model '
412	DO ialert_id = 1, inb_altests
413	WRITE(numout,*) ialert_id, cl_alname(ialert_id)//' : ', inb_alp(ialert_id), ' times ! '
414	END DO
415	ENDIF
416	!
417	END SUBROUTINE ice_ctl
418
419
420	SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 )
421	!!-------------------------------------------------------------------
422	!! * ROUTINE ice_prt *
423	!!
424	!! ** Purpose : Writes global ice state on the (i,j) point
425	!! in ocean.ouput
426	!! 3 possibilities exist
427	!! n = 1/-1 -> simple ice state (plus Mechanical Check if -1)
428	!! n = 2 -> exhaustive state
429	!! n = 3 -> ice/ocean salt fluxes
430	!!
431	!! ** input : point coordinates (i,j)
432	!! n : number of the option
433	!!-------------------------------------------------------------------
434	INTEGER , INTENT(in) :: kt ! ocean time step
435	INTEGER , INTENT(in) :: ki, kj, kn ! ocean gridpoint indices
436	CHARACTER(len=*), INTENT(in) :: cd1 !
437	!!
438	INTEGER :: jl, ji, jj
439	!!-------------------------------------------------------------------
440
441	DO ji = mi0(ki), mi1(ki)
442	DO jj = mj0(kj), mj1(kj)
443
444	WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title
445
446	!----------------
447	! Simple state
448	!----------------
449
450	IF ( kn == 1 .OR. kn == -1 ) THEN
451	WRITE(numout,*) ' ice_prt - Point : ',ji,jj
452	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
453	WRITE(numout,*) ' Simple state '
454	WRITE(numout,*) ' masks s,u,v : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1)
455	WRITE(numout,*) ' lat - long : ', gphit(ji,jj), glamt(ji,jj)
456	WRITE(numout,*) ' - Ice drift '
457	WRITE(numout,*) ' ~~~~~~~~~~~ '
458	WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj)
459	WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj)
460	WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ji,jj-1)
461	WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj)
462	WRITE(numout,*) ' strength : ', strength(ji,jj)
463	WRITE(numout,*)
464	WRITE(numout,*) ' - Cell values '
465	WRITE(numout,*) ' ~~~~~~~~~~~ '
466	WRITE(numout,*) ' cell area : ', e1e2t(ji,jj)
467	WRITE(numout,*) ' at_i : ', at_i(ji,jj)
468	WRITE(numout,*) ' vt_i : ', vt_i(ji,jj)
469	WRITE(numout,*) ' vt_s : ', vt_s(ji,jj)
470	DO jl = 1, jpl
471	WRITE(numout,*) ' - Category (', jl,')'
472	WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl)
473	WRITE(numout,*) ' h_i : ', h_i(ji,jj,jl)
474	WRITE(numout,*) ' h_s : ', h_s(ji,jj,jl)
475	WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl)
476	WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl)
477	WRITE(numout,*) ' e_s : ', e_s(ji,jj,1,jl)
478	WRITE(numout,*) ' e_i : ', e_i(ji,jj,1:nlay_i,jl)
479	WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl)
480	WRITE(numout,*) ' t_snow : ', t_s(ji,jj,1,jl)
481	WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl)
482	WRITE(numout,*) ' s_i : ', s_i(ji,jj,jl)
483	WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl)
484	WRITE(numout,*)
485	END DO
486	ENDIF
487	IF( kn == -1 ) THEN
488	WRITE(numout,) ' Mechanical Check ************* '
489	WRITE(numout,*) ' Check what means ice divergence '
490	WRITE(numout,*) ' Total ice concentration ', at_i (ji,jj)
491	WRITE(numout,*) ' Total lead fraction ', ato_i(ji,jj)
492	WRITE(numout,*) ' Sum of both ', ato_i(ji,jj) + at_i(ji,jj)
493	WRITE(numout,*) ' Sum of both minus 1 ', ato_i(ji,jj) + at_i(ji,jj) - 1.00
494	ENDIF
495
496
497	!--------------------
498	! Exhaustive state
499	!--------------------
500
501	IF ( kn .EQ. 2 ) THEN
502	WRITE(numout,*) ' ice_prt - Point : ',ji,jj
503	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
504	WRITE(numout,*) ' Exhaustive state '
505	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
506	WRITE(numout,*)
507	WRITE(numout,*) ' - Cell values '
508	WRITE(numout,*) ' ~~~~~~~~~~~ '
509	WRITE(numout,*) ' cell area : ', e1e2t(ji,jj)
510	WRITE(numout,*) ' at_i : ', at_i(ji,jj)
511	WRITE(numout,*) ' vt_i : ', vt_i(ji,jj)
512	WRITE(numout,*) ' vt_s : ', vt_s(ji,jj)
513	WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj)
514	WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj)
515	WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ji,jj-1)
516	WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj)
517	WRITE(numout,*) ' strength : ', strength(ji,jj)
518	WRITE(numout,*) ' u_ice_b : ', u_ice_b(ji,jj) , ' v_ice_b : ', v_ice_b(ji,jj)
519	WRITE(numout,*)
520
521	DO jl = 1, jpl
522	WRITE(numout,*) ' - Category (',jl,')'
523	WRITE(numout,*) ' ~~~~~~~~ '
524	WRITE(numout,*) ' h_i : ', h_i(ji,jj,jl) , ' h_s : ', h_s(ji,jj,jl)
525	WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl)
526	WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl) , ' t_s : ', t_s(ji,jj,1,jl)
527	WRITE(numout,*) ' s_i : ', s_i(ji,jj,jl) , ' o_i : ', o_i(ji,jj,jl)
528	WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) , ' a_i_b : ', a_i_b(ji,jj,jl)
529	WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) , ' v_i_b : ', v_i_b(ji,jj,jl)
530	WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) , ' v_s_b : ', v_s_b(ji,jj,jl)
531	WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl) , ' ei1 : ', e_i_b(ji,jj,1,jl)
532	WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl) , ' ei2_b : ', e_i_b(ji,jj,2,jl)
533	WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl)
534	WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl) , ' sv_i_b : ', sv_i_b(ji,jj,jl)
535	WRITE(numout,*) ' oa_i : ', oa_i(ji,jj,jl) , ' oa_i_b : ', oa_i_b(ji,jj,jl)
536	END DO !jl
537
538	WRITE(numout,*)
539	WRITE(numout,*) ' - Heat / FW fluxes '
540	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
541	WRITE(numout,) ' - Heat fluxes in and out the ice **'
542	WRITE(numout,) ' qsr_ini : ', (1._wp-at_i_b(ji,jj)) qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) )
543	WRITE(numout,) ' qns_ini : ', (1._wp-at_i_b(ji,jj)) qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) )
544	WRITE(numout,*)
545	WRITE(numout,*)
546	WRITE(numout,*) ' sst : ', sst_m(ji,jj)
547	WRITE(numout,*) ' sss : ', sss_m(ji,jj)
548	WRITE(numout,*)
549	WRITE(numout,*) ' - Stresses '
550	WRITE(numout,*) ' ~~~~~~~~ '
551	WRITE(numout,*) ' utau_ice : ', utau_ice(ji,jj)
552	WRITE(numout,*) ' vtau_ice : ', vtau_ice(ji,jj)
553	WRITE(numout,*) ' utau : ', utau (ji,jj)
554	WRITE(numout,*) ' vtau : ', vtau (ji,jj)
555	ENDIF
556
557	!---------------------
558	! Salt / heat fluxes
559	!---------------------
560
561	IF ( kn .EQ. 3 ) THEN
562	WRITE(numout,*) ' ice_prt - Point : ',ji,jj
563	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
564	WRITE(numout,*) ' - Salt / Heat Fluxes '
565	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
566	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
567	WRITE(numout,*)
568	WRITE(numout,) ' - Heat fluxes at bottom interface **'
569	WRITE(numout,*) ' qsr : ', qsr(ji,jj)
570	WRITE(numout,*) ' qns : ', qns(ji,jj)
571	WRITE(numout,*)
572	WRITE(numout,*) ' hfx_mass : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj)
573	WRITE(numout,*) ' hfx_in : ', hfx_in(ji,jj)
574	WRITE(numout,*) ' hfx_out : ', hfx_out(ji,jj)
575	WRITE(numout,*) ' dhc : ', diag_heat(ji,jj)
576	WRITE(numout,*)
577	WRITE(numout,*) ' hfx_dyn : ', hfx_dyn(ji,jj)
578	WRITE(numout,*) ' hfx_thd : ', hfx_thd(ji,jj)
579	WRITE(numout,*) ' hfx_res : ', hfx_res(ji,jj)
580	WRITE(numout,*) ' fhtur : ', fhtur(ji,jj)
581	WRITE(numout,) ' qlead : ', qlead(ji,jj) r1_rdtice
582	WRITE(numout,*)
583	WRITE(numout,) ' - Salt fluxes at bottom interface **'
584	WRITE(numout,*) ' emp : ', emp (ji,jj)
585	WRITE(numout,*) ' sfx : ', sfx (ji,jj)
586	WRITE(numout,*) ' sfx_res : ', sfx_res(ji,jj)
587	WRITE(numout,*) ' sfx_bri : ', sfx_bri(ji,jj)
588	WRITE(numout,*) ' sfx_dyn : ', sfx_dyn(ji,jj)
589	WRITE(numout,*)
590	WRITE(numout,*) ' - Momentum fluxes '
591	WRITE(numout,*) ' utau : ', utau(ji,jj)
592	WRITE(numout,*) ' vtau : ', vtau(ji,jj)
593	ENDIF
594	WRITE(numout,*) ' '
595	!
596	END DO
597	END DO
598	!
599	END SUBROUTINE ice_prt
600
601	SUBROUTINE ice_prt3D( cd_routine )
602	!!-------------------------------------------------------------------
603	!! * ROUTINE ice_prt3D *
604	!!
605	!! ** Purpose : CTL prints of ice arrays in case ln_ctl is activated
606	!!
607	!!-------------------------------------------------------------------
608	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
609	INTEGER :: jk, jl ! dummy loop indices
610
611	CALL prt_ctl_info(' ========== ')
612	CALL prt_ctl_info( cd_routine )
613	CALL prt_ctl_info(' ========== ')
614	CALL prt_ctl_info(' - Cell values : ')
615	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
616	CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' cell area :')
617	CALL prt_ctl(tab2d_1=at_i , clinfo1=' at_i :')
618	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' ato_i :')
619	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' vt_i :')
620	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' vt_s :')
621	CALL prt_ctl(tab2d_1=divu_i , clinfo1=' divu_i :')
622	CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :')
623	CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' stress1_i :')
624	CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' stress2_i :')
625	CALL prt_ctl(tab2d_1=stress12_i , clinfo1=' stress12_i :')
626	CALL prt_ctl(tab2d_1=strength , clinfo1=' strength :')
627	CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :')
628	CALL prt_ctl(tab2d_1=u_ice , clinfo1=' u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
629
630	DO jl = 1, jpl
631	CALL prt_ctl_info(' ')
632	CALL prt_ctl_info(' - Category : ', ivar1=jl)
633	CALL prt_ctl_info(' ~~~~~~~~~~')
634	CALL prt_ctl(tab2d_1=h_i (:,:,jl) , clinfo1= ' h_i : ')
635	CALL prt_ctl(tab2d_1=h_s (:,:,jl) , clinfo1= ' h_s : ')
636	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' t_su : ')
637	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' t_snow : ')
638	CALL prt_ctl(tab2d_1=s_i (:,:,jl) , clinfo1= ' s_i : ')
639	CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' o_i : ')
640	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' a_i : ')
641	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' v_i : ')
642	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' v_s : ')
643	CALL prt_ctl(tab2d_1=e_i (:,:,1,jl) , clinfo1= ' e_i1 : ')
644	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' e_snow : ')
645	CALL prt_ctl(tab2d_1=sv_i (:,:,jl) , clinfo1= ' sv_i : ')
646	CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' oa_i : ')
647
648	DO jk = 1, nlay_i
649	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
650	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i : ')
651	END DO
652	END DO
653
654	CALL prt_ctl_info(' ')
655	CALL prt_ctl_info(' - Heat / FW fluxes : ')
656	CALL prt_ctl_info(' ~~~~~~~~~~~~~~~~~~ ')
657	CALL prt_ctl(tab2d_1=sst_m , clinfo1= ' sst : ', tab2d_2=sss_m , clinfo2= ' sss : ')
658	CALL prt_ctl(tab2d_1=qsr , clinfo1= ' qsr : ', tab2d_2=qns , clinfo2= ' qns : ')
659	CALL prt_ctl(tab2d_1=emp , clinfo1= ' emp : ', tab2d_2=sfx , clinfo2= ' sfx : ')
660
661	CALL prt_ctl_info(' ')
662	CALL prt_ctl_info(' - Stresses : ')
663	CALL prt_ctl_info(' ~~~~~~~~~~ ')
664	CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', tab2d_2=vtau , clinfo2= ' vtau : ')
665	CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ')
666
667	END SUBROUTINE ice_prt3D
668
669	#else
670	!!----------------------------------------------------------------------
671	!! Default option Empty Module No ESIM sea-ice model
672	!!----------------------------------------------------------------------
673	#endif
674
675	!!======================================================================
676	END MODULE icectl

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

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

Context Navigation

source: branches/UKMO/dev_r8126_LIM3_couple/NEMOGCM/NEMO/LIM_SRC_3/icectl.F90 @ 8879

Download in other formats: