Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

icectl.F90

Last change on this file was 14590, checked in by clem, 3 years ago
4.0-HEAD: solve ticket #2627 to allow simulations with conductive fluxes instead of normal fluxes on top of sea ice (MetO requirement). I do not think I tackled all the issues but this is the best I can do without having a proper configuration to test it.
Property svn:keywords set to `Id`
File size: 47.5 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	!! 4.0 ! 2018 (many people) SI3 [aka Sea Ice cube]
9	!!----------------------------------------------------------------------
10	#if defined key_si3
11	!!----------------------------------------------------------------------
12	!! 'key_si3' SI3 sea-ice model
13	!!----------------------------------------------------------------------
14	!! ice_cons_hsm : conservation tests on heat, salt and mass during a time step (global)
15	!! ice_cons_final : conservation tests on heat, salt and mass at end of time step (global)
16	!! ice_cons2D : conservation tests on heat, salt and mass at each gridcell
17	!! ice_ctl : control prints in case of crash
18	!! ice_prt : control prints at a given grid point
19	!! ice_prt3D : control prints of ice arrays
20	!!----------------------------------------------------------------------
21	USE phycst ! physical constants
22	USE oce ! ocean dynamics and tracers
23	USE dom_oce ! ocean space and time domain
24	USE ice ! sea-ice: variables
25	USE ice1D ! sea-ice: thermodynamics variables
26	USE sbc_oce ! Surface boundary condition: ocean fields
27	USE sbc_ice ! Surface boundary condition: ice fields
28	!
29	USE in_out_manager ! I/O manager
30	USE iom ! I/O manager library
31	USE lib_mpp ! MPP library
32	USE lib_fortran ! fortran utilities (glob_sum + no signed zero)
33	USE timing ! Timing
34	USE prtctl ! Print control
35
36	IMPLICIT NONE
37	PRIVATE
38
39	PUBLIC ice_cons_hsm
40	PUBLIC ice_cons_final
41	PUBLIC ice_cons2D
42	PUBLIC ice_ctl
43	PUBLIC ice_prt
44	PUBLIC ice_prt3D
45	PUBLIC ice_drift_wri
46	PUBLIC ice_drift_init
47
48	! thresold rates for conservation
49	! these values are changed by the namelist parameter rn_icechk, so that threshold = zchk * rn_icechk
50	REAL(wp), PARAMETER :: zchk_m = 2.5e-7 ! kg/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost
51	REAL(wp), PARAMETER :: zchk_s = 2.5e-6 ! g/m2/s <=> 1e-6 m of ice per hour spuriously gained/lost (considering s=10g/kg)
52	REAL(wp), PARAMETER :: zchk_t = 7.5e-2 ! W/m2 <=> 1e-6 m of ice per hour spuriously gained/lost (considering Lf=3e5J/kg)
53
54	! for drift outputs
55	CHARACTER(LEN=50) :: clname="icedrift_diagnostics.ascii" ! ascii filename
56	INTEGER :: numicedrift ! outfile unit
57	REAL(wp) :: rdiag_icemass, rdiag_icesalt, rdiag_iceheat
58	REAL(wp) :: rdiag_adv_icemass, rdiag_adv_icesalt, rdiag_adv_iceheat
59
60	!! * Substitutions
61	# include "vectopt_loop_substitute.h90"
62	!!----------------------------------------------------------------------
63	!! NEMO/ICE 4.0 , NEMO Consortium (2018)
64	!! $Id$
65	!! Software governed by the CeCILL license (see ./LICENSE)
66	!!----------------------------------------------------------------------
67	CONTAINS
68
69	SUBROUTINE ice_cons_hsm( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft )
70	!!-------------------------------------------------------------------
71	!! * ROUTINE ice_cons_hsm *
72	!!
73	!! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
74	!! + test if ice concentration and volume are > 0
75	!!
76	!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
77	!! It prints in ocean.output if there is a violation of conservation at each time-step
78	!! The thresholds (zchk_m, zchk_s, zchk_t) determine violations
79	!! For salt and heat thresholds, ice is considered to have a salinity of 10
80	!! and a heat content of 3e5 J/kg (=latent heat of fusion)
81	!!-------------------------------------------------------------------
82	INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end
83	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
84	REAL(wp) , INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
85	!!
86	REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, &
87	& zdiag_vimin, zdiag_vsmin, zdiag_vpmin, zdiag_vlmin, zdiag_aimin, zdiag_aimax, &
88	& zdiag_eimin, zdiag_esmin, zdiag_simin
89	REAL(wp) :: zvtrp, zetrp
90	REAL(wp) :: zarea
91	!!-------------------------------------------------------------------
92	!
93	IF( icount == 0 ) THEN
94
95	pdiag_v = glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t )
96	pdiag_s = glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t )
97	pdiag_t = glob_sum( 'icectl', ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) ) * e1e2t )
98
99	! mass flux
100	pdiag_fv = glob_sum( 'icectl', &
101	& ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
102	& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) * e1e2t )
103	! salt flux
104	pdiag_fs = glob_sum( 'icectl', &
105	& ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + &
106	& sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t )
107	! heat flux
108	pdiag_ft = glob_sum( 'icectl', &
109	& ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
110	& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t )
111
112	ELSEIF( icount == 1 ) THEN
113
114	! -- mass diag -- !
115	zdiag_mass = ( glob_sum( 'icectl', SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) * e1e2t ) &
116	& - pdiag_v ) * r1_rdtice &
117	& + glob_sum( 'icectl', ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + &
118	& wfx_lam + wfx_pnd + wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + &
119	& wfx_ice_sub + wfx_spr ) * e1e2t ) &
120	& - pdiag_fv
121	!
122	! -- salt diag -- !
123	zdiag_salt = ( glob_sum( 'icectl', SUM( sv_i * rhoi , dim=3 ) * e1e2t ) - pdiag_s ) * r1_rdtice &
124	& + glob_sum( 'icectl', ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + &
125	& sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) * e1e2t ) &
126	& - pdiag_fs
127	!
128	! -- heat diag -- !
129	zdiag_heat = ( glob_sum( 'icectl', ( SUM(SUM(e_i, dim=4), dim=3) + SUM(SUM(e_s, dim=4), dim=3) ) * e1e2t ) - pdiag_t &
130	& ) * r1_rdtice &
131	& + glob_sum( 'icectl', ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
132	& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) * e1e2t ) &
133	& - pdiag_ft
134
135	! -- min/max diag -- !
136	zdiag_aimax = glob_max( 'icectl', SUM( a_i, dim=3 ) )
137	zdiag_vimin = glob_min( 'icectl', v_i )
138	zdiag_vsmin = glob_min( 'icectl', v_s )
139	zdiag_vpmin = glob_min( 'icectl', v_ip )
140	zdiag_vlmin = glob_min( 'icectl', v_il )
141	zdiag_aimin = glob_min( 'icectl', a_i )
142	zdiag_simin = glob_min( 'icectl', sv_i )
143	zdiag_eimin = glob_min( 'icectl', SUM( e_i, dim=3 ) )
144	zdiag_esmin = glob_min( 'icectl', SUM( e_s, dim=3 ) )
145
146	! -- advection scheme is conservative? -- !
147	zvtrp = glob_sum( 'icectl', diag_adv_mass * e1e2t )
148	zetrp = glob_sum( 'icectl', diag_adv_heat * e1e2t )
149
150	! ice area (+epsi10 to set a threshold > 0 when there is no ice)
151	zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t )
152
153	IF( lwp ) THEN
154	! check conservation issues
155	IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) &
156	& WRITE(numout,) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass rdt_ice
157	IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) &
158	& WRITE(numout,) cd_routine,' : violation salt cons. [g] = ',zdiag_salt rdt_ice
159	IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) &
160	& WRITE(numout,) cd_routine,' : violation heat cons. [J] = ',zdiag_heat rdt_ice
161	! check negative values
162	IF( zdiag_vimin < 0. ) WRITE(numout,*) cd_routine,' : violation v_i < 0 = ',zdiag_vimin
163	IF( zdiag_vsmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_s < 0 = ',zdiag_vsmin
164	IF( zdiag_vpmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_ip < 0 = ',zdiag_vpmin
165	IF( zdiag_vlmin < 0. ) WRITE(numout,*) cd_routine,' : violation v_il < 0 = ',zdiag_vlmin
166	IF( zdiag_aimin < 0. ) WRITE(numout,*) cd_routine,' : violation a_i < 0 = ',zdiag_aimin
167	IF( zdiag_simin < 0. ) WRITE(numout,*) cd_routine,' : violation s_i < 0 = ',zdiag_simin
168	IF( zdiag_eimin < 0. ) WRITE(numout,*) cd_routine,' : violation e_i < 0 = ',zdiag_eimin
169	IF( zdiag_esmin < 0. ) WRITE(numout,*) cd_routine,' : violation e_s < 0 = ',zdiag_esmin
170	! check maximum ice concentration
171	IF( zdiag_aimax>MAX(rn_amax_n,rn_amax_s)+epsi10 .AND. cd_routine /= 'icedyn_adv' .AND. cd_routine /= 'icedyn_rdgrft' ) &
172	& WRITE(numout,*) cd_routine,' : violation a_i > amax = ',zdiag_aimax
173	! check if advection scheme is conservative
174	IF( ABS(zvtrp) > zchk_m * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
175	& WRITE(numout,) cd_routine,' : violation adv scheme [kg] = ',zvtrp rdt_ice
176	IF( ABS(zetrp) > zchk_t * rn_icechk_glo * zarea .AND. cd_routine == 'icedyn_adv' ) &
177	& WRITE(numout,) cd_routine,' : violation adv scheme [J] = ',zetrp rdt_ice
178	ENDIF
179	!
180	ENDIF
181
182	END SUBROUTINE ice_cons_hsm
183
184	SUBROUTINE ice_cons_final( cd_routine )
185	!!-------------------------------------------------------------------
186	!! * ROUTINE ice_cons_final *
187	!!
188	!! ** Purpose : Test the conservation of heat, salt and mass at the end of each ice time-step
189	!!
190	!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
191	!! It prints in ocean.output if there is a violation of conservation at each time-step
192	!! The thresholds (zchk_m, zchk_s, zchk_t) determine the violations
193	!! For salt and heat thresholds, ice is considered to have a salinity of 10
194	!! and a heat content of 3e5 J/kg (=latent heat of fusion)
195	!!-------------------------------------------------------------------
196	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
197	REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat
198	REAL(wp) :: zarea
199	!!-------------------------------------------------------------------
200
201	! water flux
202	! -- mass diag -- !
203	zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd &
204	& + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t )
205
206	! -- salt diag -- !
207	zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t )
208
209	! -- heat diag -- !
210	zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t )
211	! equivalent to this:
212	!!zdiag_heat = glob_sum( 'icectl', ( -diag_heat + hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
213	!! & - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr &
214	!! & ) * e1e2t )
215
216	! ice area (+epsi10 to set a threshold > 0 when there is no ice)
217	zarea = glob_sum( 'icectl', SUM( a_i + epsi10, dim=3 ) * e1e2t )
218
219	IF( lwp ) THEN
220	IF( ABS(zdiag_mass) > zchk_m * rn_icechk_glo * zarea ) &
221	& WRITE(numout,) cd_routine,' : violation mass cons. [kg] = ',zdiag_mass rdt_ice
222	IF( ABS(zdiag_salt) > zchk_s * rn_icechk_glo * zarea ) &
223	& WRITE(numout,) cd_routine,' : violation salt cons. [g] = ',zdiag_salt rdt_ice
224	IF( ABS(zdiag_heat) > zchk_t * rn_icechk_glo * zarea ) &
225	& WRITE(numout,) cd_routine,' : violation heat cons. [J] = ',zdiag_heat rdt_ice
226	ENDIF
227	!
228	END SUBROUTINE ice_cons_final
229
230	SUBROUTINE ice_cons2D( icount, cd_routine, pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft )
231	!!-------------------------------------------------------------------
232	!! * ROUTINE ice_cons2D *
233	!!
234	!! ** Purpose : Test the conservation of heat, salt and mass for each ice routine
235	!! + test if ice concentration and volume are > 0
236	!!
237	!! ** Method : This is an online diagnostics which can be activated with ln_icediachk=true
238	!! It stops the code if there is a violation of conservation at any gridcell
239	!!-------------------------------------------------------------------
240	INTEGER , INTENT(in) :: icount ! called at: =0 the begining of the routine, =1 the end
241	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
242	REAL(wp) , DIMENSION(jpi,jpj), INTENT(inout) :: pdiag_v, pdiag_s, pdiag_t, pdiag_fv, pdiag_fs, pdiag_ft
243	!!
244	REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass, zdiag_salt, zdiag_heat, &
245	& zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin !!, zdiag_amax
246	INTEGER :: jl, jk
247	LOGICAL :: ll_stop_m = .FALSE.
248	LOGICAL :: ll_stop_s = .FALSE.
249	LOGICAL :: ll_stop_t = .FALSE.
250	CHARACTER(len=120) :: clnam ! filename for the output
251	!!-------------------------------------------------------------------
252	!
253	IF( icount == 0 ) THEN
254
255	pdiag_v = SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 )
256	pdiag_s = SUM( sv_i * rhoi , dim=3 )
257	pdiag_t = SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 )
258
259	! mass flux
260	pdiag_fv = wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
261	& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr
262	! salt flux
263	pdiag_fs = sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam
264	! heat flux
265	pdiag_ft = hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
266	& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr
267
268	ELSEIF( icount == 1 ) THEN
269
270	! -- mass diag -- !
271	zdiag_mass = ( SUM( v_i * rhoi + v_s * rhos + ( v_ip + v_il ) * rhow, dim=3 ) - pdiag_v ) * r1_rdtice &
272	& + ( wfx_bog + wfx_bom + wfx_sum + wfx_sni + wfx_opw + wfx_res + wfx_dyn + wfx_lam + wfx_pnd + &
273	& wfx_snw_sni + wfx_snw_sum + wfx_snw_dyn + wfx_snw_sub + wfx_ice_sub + wfx_spr ) &
274	& - pdiag_fv
275	IF( MAXVAL( ABS(zdiag_mass) ) > zchk_m * rn_icechk_cel ) ll_stop_m = .TRUE.
276	!
277	! -- salt diag -- !
278	zdiag_salt = ( SUM( sv_i * rhoi , dim=3 ) - pdiag_s ) * r1_rdtice &
279	& + ( sfx_bri + sfx_bog + sfx_bom + sfx_sum + sfx_sni + sfx_opw + sfx_res + sfx_dyn + sfx_sub + sfx_lam ) &
280	& - pdiag_fs
281	IF( MAXVAL( ABS(zdiag_salt) ) > zchk_s * rn_icechk_cel ) ll_stop_s = .TRUE.
282	!
283	! -- heat diag -- !
284	zdiag_heat = ( SUM( SUM( e_i, dim=4 ), dim=3 ) + SUM( SUM( e_s, dim=4 ), dim=3 ) - pdiag_t ) * r1_rdtice &
285	& + ( hfx_sum + hfx_bom + hfx_bog + hfx_dif + hfx_opw + hfx_snw &
286	& - hfx_thd - hfx_dyn - hfx_res - hfx_sub - hfx_spr ) &
287	& - pdiag_ft
288	IF( MAXVAL( ABS(zdiag_heat) ) > zchk_t * rn_icechk_cel ) ll_stop_t = .TRUE.
289	!
290	! -- other diags -- !
291	! a_i < 0
292	zdiag_amin(:,:) = 0._wp
293	DO jl = 1, jpl
294	WHERE( a_i(:,:,jl) < 0._wp ) zdiag_amin(:,:) = 1._wp
295	ENDDO
296	! v_i < 0
297	zdiag_vmin(:,:) = 0._wp
298	DO jl = 1, jpl
299	WHERE( v_i(:,:,jl) < 0._wp ) zdiag_vmin(:,:) = 1._wp
300	ENDDO
301	! s_i < 0
302	zdiag_smin(:,:) = 0._wp
303	DO jl = 1, jpl
304	WHERE( s_i(:,:,jl) < 0._wp ) zdiag_smin(:,:) = 1._wp
305	ENDDO
306	! e_i < 0
307	zdiag_emin(:,:) = 0._wp
308	DO jl = 1, jpl
309	DO jk = 1, nlay_i
310	WHERE( e_i(:,:,jk,jl) < 0._wp ) zdiag_emin(:,:) = 1._wp
311	ENDDO
312	ENDDO
313	! a_i > amax
314	!WHERE( SUM( a_i, dim=3 ) > ( MAX( rn_amax_n, rn_amax_s ) + epsi10 ) ; zdiag_amax(:,:) = 1._wp
315	!ELSEWHERE ; zdiag_amax(:,:) = 0._wp
316	!END WHERE
317
318	IF( ll_stop_m .OR. ll_stop_s .OR. ll_stop_t ) THEN
319	clnam = 'diag_ice_conservation_'//cd_routine
320	CALL ice_cons_wri( clnam, zdiag_mass, zdiag_salt, zdiag_heat, zdiag_amin, zdiag_vmin, zdiag_smin, zdiag_emin )
321	ENDIF
322
323	IF( ll_stop_m ) CALL ctl_stop( 'STOP', cd_routine//': ice mass conservation issue' )
324	IF( ll_stop_s ) CALL ctl_stop( 'STOP', cd_routine//': ice salt conservation issue' )
325	IF( ll_stop_t ) CALL ctl_stop( 'STOP', cd_routine//': ice heat conservation issue' )
326
327	ENDIF
328
329	END SUBROUTINE ice_cons2D
330
331	SUBROUTINE ice_cons_wri( cdfile_name, pdiag_mass, pdiag_salt, pdiag_heat, pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin )
332	!!---------------------------------------------------------------------
333	!! * ROUTINE ice_cons_wri *
334	!!
335	!! ** Purpose : create a NetCDF file named cdfile_name which contains
336	!! the instantaneous fields when conservation issue occurs
337	!!
338	!! ** Method : NetCDF files using ioipsl
339	!!----------------------------------------------------------------------
340	CHARACTER(len=*), INTENT( in ) :: cdfile_name ! name of the file created
341	REAL(wp), DIMENSION(:,:), INTENT( in ) :: pdiag_mass, pdiag_salt, pdiag_heat, &
342	& pdiag_amin, pdiag_vmin, pdiag_smin, pdiag_emin !!, pdiag_amax
343	!!
344	INTEGER :: inum
345	!!----------------------------------------------------------------------
346	!
347	IF(lwp) WRITE(numout,*)
348	IF(lwp) WRITE(numout,*) 'ice_cons_wri : single instantaneous ice state'
349	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ named :', cdfile_name, '...nc'
350	IF(lwp) WRITE(numout,*)
351
352	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
353
354	CALL iom_rstput( 0, 0, inum, 'cons_mass', pdiag_mass(:,:) , ktype = jp_r8 ) ! ice mass spurious lost/gain
355	CALL iom_rstput( 0, 0, inum, 'cons_salt', pdiag_salt(:,:) , ktype = jp_r8 ) ! ice salt spurious lost/gain
356	CALL iom_rstput( 0, 0, inum, 'cons_heat', pdiag_heat(:,:) , ktype = jp_r8 ) ! ice heat spurious lost/gain
357	! other diags
358	CALL iom_rstput( 0, 0, inum, 'aneg_count', pdiag_amin(:,:) , ktype = jp_r8 ) !
359	CALL iom_rstput( 0, 0, inum, 'vneg_count', pdiag_vmin(:,:) , ktype = jp_r8 ) !
360	CALL iom_rstput( 0, 0, inum, 'sneg_count', pdiag_smin(:,:) , ktype = jp_r8 ) !
361	CALL iom_rstput( 0, 0, inum, 'eneg_count', pdiag_emin(:,:) , ktype = jp_r8 ) !
362	! mean state
363	CALL iom_rstput( 0, 0, inum, 'icecon' , SUM(a_i ,dim=3) , ktype = jp_r8 ) !
364	CALL iom_rstput( 0, 0, inum, 'icevol' , SUM(v_i ,dim=3) , ktype = jp_r8 ) !
365	CALL iom_rstput( 0, 0, inum, 'snwvol' , SUM(v_s ,dim=3) , ktype = jp_r8 ) !
366	CALL iom_rstput( 0, 0, inum, 'pndvol' , SUM(v_ip,dim=3) , ktype = jp_r8 ) !
367	CALL iom_rstput( 0, 0, inum, 'lidvol' , SUM(v_il,dim=3) , ktype = jp_r8 ) !
368
369	CALL iom_close( inum )
370
371	END SUBROUTINE ice_cons_wri
372
373	SUBROUTINE ice_ctl( kt )
374	!!-------------------------------------------------------------------
375	!! * ROUTINE ice_ctl *
376	!!
377	!! ** Purpose : control checks
378	!!-------------------------------------------------------------------
379	INTEGER, INTENT(in) :: kt ! ocean time step
380	INTEGER :: ja, ji, jj, jk, jl ! dummy loop indices
381	INTEGER :: ialert_id ! number of the current alert
382	REAL(wp) :: ztmelts ! ice layer melting point
383	CHARACTER (len=30), DIMENSION(20) :: cl_alname ! name of alert
384	INTEGER , DIMENSION(20) :: inb_alp ! number of alerts positive
385	!!-------------------------------------------------------------------
386	inb_alp(:) = 0
387	ialert_id = 0
388
389	! Alert if very high salinity
390	ialert_id = ialert_id + 1 ! reference number of this alert
391	cl_alname(ialert_id) = ' Very high salinity ' ! name of the alert
392	DO jl = 1, jpl
393	DO jj = 1, jpj
394	DO ji = 1, jpi
395	IF( v_i(ji,jj,jl) > epsi10 ) THEN
396	IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) > rn_simax ) THEN
397	WRITE(numout,*) ' ALERTE : Very high salinity ',sv_i(ji,jj,jl)/v_i(ji,jj,jl)
398	WRITE(numout,*) ' at i,j,l = ',ji,jj,jl
399	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
400	ENDIF
401	ENDIF
402	END DO
403	END DO
404	END DO
405
406	! Alert if very low salinity
407	ialert_id = ialert_id + 1 ! reference number of this alert
408	cl_alname(ialert_id) = ' Very low salinity ' ! name of the alert
409	DO jl = 1, jpl
410	DO jj = 1, jpj
411	DO ji = 1, jpi
412	IF( v_i(ji,jj,jl) > epsi10 ) THEN
413	IF( sv_i(ji,jj,jl) / v_i(ji,jj,jl) < rn_simin ) THEN
414	WRITE(numout,*) ' ALERTE : Very low salinity ',sv_i(ji,jj,jl),v_i(ji,jj,jl)
415	WRITE(numout,*) ' at i,j,l = ',ji,jj,jl
416	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
417	ENDIF
418	ENDIF
419	END DO
420	END DO
421	END DO
422
423	! Alert if very cold ice
424	ialert_id = ialert_id + 1 ! reference number of this alert
425	cl_alname(ialert_id) = ' Very cold ice ' ! name of the alert
426	DO jl = 1, jpl
427	DO jk = 1, nlay_i
428	DO jj = 1, jpj
429	DO ji = 1, jpi
430	ztmelts = -rTmlt * sz_i(ji,jj,jk,jl) + rt0
431	IF( t_i(ji,jj,jk,jl) < -50.+rt0 .AND. v_i(ji,jj,jl) > epsi10 ) THEN
432	WRITE(numout,*) ' ALERTE : Very cold ice ',(t_i(ji,jj,jk,jl)-rt0)
433	WRITE(numout,*) ' at i,j,k,l = ',ji,jj,jk,jl
434	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
435	ENDIF
436	END DO
437	END DO
438	END DO
439	END DO
440
441	! Alert if very warm ice
442	ialert_id = ialert_id + 1 ! reference number of this alert
443	cl_alname(ialert_id) = ' Very warm ice ' ! name of the alert
444	DO jl = 1, jpl
445	DO jk = 1, nlay_i
446	DO jj = 1, jpj
447	DO ji = 1, jpi
448	ztmelts = -rTmlt * sz_i(ji,jj,jk,jl) + rt0
449	IF( t_i(ji,jj,jk,jl) > ztmelts .AND. v_i(ji,jj,jl) > epsi10 ) THEN
450	WRITE(numout,*) ' ALERTE : Very warm ice',(t_i(ji,jj,jk,jl)-rt0)
451	WRITE(numout,*) ' at i,j,k,l = ',ji,jj,jk,jl
452	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
453	ENDIF
454	END DO
455	END DO
456	END DO
457	END DO
458
459	! Alerte if very thick ice
460	ialert_id = ialert_id + 1 ! reference number of this alert
461	cl_alname(ialert_id) = ' Very thick ice ' ! name of the alert
462	jl = jpl
463	DO jj = 1, jpj
464	DO ji = 1, jpi
465	IF( h_i(ji,jj,jl) > 50._wp ) THEN
466	WRITE(numout,*) ' ALERTE : Very thick ice ',h_i(ji,jj,jl)
467	WRITE(numout,*) ' at i,j,l = ',ji,jj,jl
468	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
469	ENDIF
470	END DO
471	END DO
472
473	! Alerte if very thin ice
474	ialert_id = ialert_id + 1 ! reference number of this alert
475	cl_alname(ialert_id) = ' Very thin ice ' ! name of the alert
476	jl = 1
477	DO jj = 1, jpj
478	DO ji = 1, jpi
479	IF( h_i(ji,jj,jl) < rn_himin ) THEN
480	WRITE(numout,*) ' ALERTE : Very thin ice ',h_i(ji,jj,jl)
481	WRITE(numout,*) ' at i,j,l = ',ji,jj,jl
482	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
483	ENDIF
484	END DO
485	END DO
486
487	! Alert if very fast ice
488	ialert_id = ialert_id + 1 ! reference number of this alert
489	cl_alname(ialert_id) = ' Very fast ice ' ! name of the alert
490	DO jj = 1, jpj
491	DO ji = 1, jpi
492	IF( MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) ) > 2. ) THEN
493	WRITE(numout,*) ' ALERTE : Very fast ice ',MAX( ABS( u_ice(ji,jj) ), ABS( v_ice(ji,jj) ) )
494	WRITE(numout,*) ' at i,j = ',ji,jj
495	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
496	ENDIF
497	END DO
498	END DO
499
500	! Alert if there is ice on continents
501	ialert_id = ialert_id + 1 ! reference number of this alert
502	cl_alname(ialert_id) = ' Ice on continents ' ! name of the alert
503	DO jj = 1, jpj
504	DO ji = 1, jpi
505	IF( tmask(ji,jj,1) == 0._wp .AND. ( at_i(ji,jj) > 0._wp .OR. vt_i(ji,jj) > 0._wp ) ) THEN
506	WRITE(numout,*) ' ALERTE : Ice on continents ',at_i(ji,jj),vt_i(ji,jj)
507	WRITE(numout,*) ' at i,j = ',ji,jj
508	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
509	ENDIF
510	END DO
511	END DO
512
513	! Alert if incompatible ice concentration and volume
514	ialert_id = ialert_id + 1 ! reference number of this alert
515	cl_alname(ialert_id) = ' Incompatible ice conc and vol ' ! name of the alert
516	DO jj = 1, jpj
517	DO ji = 1, jpi
518	IF( ( vt_i(ji,jj) == 0._wp .AND. at_i(ji,jj) > 0._wp ) .OR. &
519	& ( vt_i(ji,jj) > 0._wp .AND. at_i(ji,jj) == 0._wp ) ) THEN
520	WRITE(numout,*) ' ALERTE : Incompatible ice conc and vol ',at_i(ji,jj),vt_i(ji,jj)
521	WRITE(numout,*) ' at i,j = ',ji,jj
522	inb_alp(ialert_id) = inb_alp(ialert_id) + 1
523	ENDIF
524	END DO
525	END DO
526
527	! sum of the alerts on all processors
528	IF( lk_mpp ) THEN
529	DO ja = 1, ialert_id
530	CALL mpp_sum('icectl', inb_alp(ja))
531	END DO
532	ENDIF
533
534	! print alerts
535	IF( lwp ) THEN
536	WRITE(numout,*) ' time step ',kt
537	WRITE(numout,*) ' All alerts at the end of ice model '
538	DO ja = 1, ialert_id
539	WRITE(numout,*) ja, cl_alname(ja)//' : ', inb_alp(ja), ' times ! '
540	END DO
541	ENDIF
542	!
543	END SUBROUTINE ice_ctl
544
545	SUBROUTINE ice_prt( kt, ki, kj, kn, cd1 )
546	!!-------------------------------------------------------------------
547	!! * ROUTINE ice_prt *
548	!!
549	!! ** Purpose : Writes global ice state on the (i,j) point
550	!! in ocean.ouput
551	!! 3 possibilities exist
552	!! n = 1/-1 -> simple ice state
553	!! n = 2 -> exhaustive state
554	!! n = 3 -> ice/ocean salt fluxes
555	!!
556	!! ** input : point coordinates (i,j)
557	!! n : number of the option
558	!!-------------------------------------------------------------------
559	INTEGER , INTENT(in) :: kt ! ocean time step
560	INTEGER , INTENT(in) :: ki, kj, kn ! ocean gridpoint indices
561	CHARACTER(len=*), INTENT(in) :: cd1 !
562	!!
563	INTEGER :: jl, ji, jj
564	!!-------------------------------------------------------------------
565
566	DO ji = mi0(ki), mi1(ki)
567	DO jj = mj0(kj), mj1(kj)
568
569	WRITE(numout,*) ' time step ',kt,' ',cd1 ! print title
570
571	!----------------
572	! Simple state
573	!----------------
574
575	IF ( kn == 1 .OR. kn == -1 ) THEN
576	WRITE(numout,*) ' ice_prt - Point : ',ji,jj
577	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
578	WRITE(numout,*) ' Simple state '
579	WRITE(numout,*) ' masks s,u,v : ', tmask(ji,jj,1), umask(ji,jj,1), vmask(ji,jj,1)
580	WRITE(numout,*) ' lat - long : ', gphit(ji,jj), glamt(ji,jj)
581	WRITE(numout,*) ' - Ice drift '
582	WRITE(numout,*) ' ~~~~~~~~~~~ '
583	WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj)
584	WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj)
585	WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ji,jj-1)
586	WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj)
587	WRITE(numout,*) ' strength : ', strength(ji,jj)
588	WRITE(numout,*) ' - Cell values '
589	WRITE(numout,*) ' ~~~~~~~~~~~ '
590	WRITE(numout,*) ' at_i : ', at_i(ji,jj)
591	WRITE(numout,*) ' ato_i : ', ato_i(ji,jj)
592	WRITE(numout,*) ' vt_i : ', vt_i(ji,jj)
593	WRITE(numout,*) ' vt_s : ', vt_s(ji,jj)
594	DO jl = 1, jpl
595	WRITE(numout,*) ' - Category (', jl,')'
596	WRITE(numout,*) ' ~~~~~~~~~~~ '
597	WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl)
598	WRITE(numout,*) ' h_i : ', h_i(ji,jj,jl)
599	WRITE(numout,*) ' h_s : ', h_s(ji,jj,jl)
600	WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl)
601	WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl)
602	WRITE(numout,*) ' e_s : ', e_s(ji,jj,1:nlay_s,jl)
603	WRITE(numout,*) ' e_i : ', e_i(ji,jj,1:nlay_i,jl)
604	WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl)
605	WRITE(numout,*) ' t_snow : ', t_s(ji,jj,1:nlay_s,jl)
606	WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl)
607	WRITE(numout,*) ' s_i : ', s_i(ji,jj,jl)
608	WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl)
609	WRITE(numout,*)
610	END DO
611	ENDIF
612
613	!--------------------
614	! Exhaustive state
615	!--------------------
616
617	IF ( kn .EQ. 2 ) THEN
618	WRITE(numout,*) ' ice_prt - Point : ',ji,jj
619	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
620	WRITE(numout,*) ' Exhaustive state '
621	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
622	WRITE(numout,*)
623	WRITE(numout,*) ' - Cell values '
624	WRITE(numout,*) ' ~~~~~~~~~~~ '
625	WRITE(numout,*) ' at_i : ', at_i(ji,jj)
626	WRITE(numout,*) ' vt_i : ', vt_i(ji,jj)
627	WRITE(numout,*) ' vt_s : ', vt_s(ji,jj)
628	WRITE(numout,*) ' u_ice(i-1,j) : ', u_ice(ji-1,jj)
629	WRITE(numout,*) ' u_ice(i ,j) : ', u_ice(ji,jj)
630	WRITE(numout,*) ' v_ice(i ,j-1): ', v_ice(ji,jj-1)
631	WRITE(numout,*) ' v_ice(i ,j) : ', v_ice(ji,jj)
632	WRITE(numout,*) ' strength : ', strength(ji,jj)
633	WRITE(numout,*)
634
635	DO jl = 1, jpl
636	WRITE(numout,*) ' - Category (',jl,')'
637	WRITE(numout,*) ' ~~~~~~~~ '
638	WRITE(numout,*) ' h_i : ', h_i(ji,jj,jl) , ' h_s : ', h_s(ji,jj,jl)
639	WRITE(numout,*) ' t_i : ', t_i(ji,jj,1:nlay_i,jl)
640	WRITE(numout,*) ' t_su : ', t_su(ji,jj,jl) , ' t_s : ', t_s(ji,jj,1:nlay_s,jl)
641	WRITE(numout,*) ' s_i : ', s_i(ji,jj,jl) , ' o_i : ', o_i(ji,jj,jl)
642	WRITE(numout,*) ' a_i : ', a_i(ji,jj,jl) , ' a_i_b : ', a_i_b(ji,jj,jl)
643	WRITE(numout,*) ' v_i : ', v_i(ji,jj,jl) , ' v_i_b : ', v_i_b(ji,jj,jl)
644	WRITE(numout,*) ' v_s : ', v_s(ji,jj,jl) , ' v_s_b : ', v_s_b(ji,jj,jl)
645	WRITE(numout,*) ' e_i1 : ', e_i(ji,jj,1,jl) , ' ei1 : ', e_i_b(ji,jj,1,jl)
646	WRITE(numout,*) ' e_i2 : ', e_i(ji,jj,2,jl) , ' ei2_b : ', e_i_b(ji,jj,2,jl)
647	WRITE(numout,*) ' e_snow : ', e_s(ji,jj,1,jl) , ' e_snow_b : ', e_s_b(ji,jj,1,jl)
648	WRITE(numout,*) ' sv_i : ', sv_i(ji,jj,jl) , ' sv_i_b : ', sv_i_b(ji,jj,jl)
649	END DO !jl
650
651	WRITE(numout,*)
652	WRITE(numout,*) ' - Heat / FW fluxes '
653	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
654	WRITE(numout,) ' - Heat fluxes in and out the ice **'
655	WRITE(numout,) ' qsr_ini : ', (1._wp-at_i_b(ji,jj)) qsr(ji,jj) + SUM( a_i_b(ji,jj,:) * qsr_ice(ji,jj,:) )
656	WRITE(numout,) ' qns_ini : ', (1._wp-at_i_b(ji,jj)) qns(ji,jj) + SUM( a_i_b(ji,jj,:) * qns_ice(ji,jj,:) )
657	WRITE(numout,*)
658	WRITE(numout,*)
659	WRITE(numout,*) ' sst : ', sst_m(ji,jj)
660	WRITE(numout,*) ' sss : ', sss_m(ji,jj)
661	WRITE(numout,*)
662	WRITE(numout,*) ' - Stresses '
663	WRITE(numout,*) ' ~~~~~~~~ '
664	WRITE(numout,*) ' utau_ice : ', utau_ice(ji,jj)
665	WRITE(numout,*) ' vtau_ice : ', vtau_ice(ji,jj)
666	WRITE(numout,*) ' utau : ', utau (ji,jj)
667	WRITE(numout,*) ' vtau : ', vtau (ji,jj)
668	ENDIF
669
670	!---------------------
671	! Salt / heat fluxes
672	!---------------------
673
674	IF ( kn .EQ. 3 ) THEN
675	WRITE(numout,*) ' ice_prt - Point : ',ji,jj
676	WRITE(numout,*) ' ~~~~~~~~~~~~~~ '
677	WRITE(numout,*) ' - Salt / Heat Fluxes '
678	WRITE(numout,*) ' ~~~~~~~~~~~~~~~~ '
679	WRITE(numout,*) ' lat - long ', gphit(ji,jj), glamt(ji,jj)
680	WRITE(numout,*)
681	WRITE(numout,) ' - Heat fluxes at bottom interface **'
682	WRITE(numout,*) ' qsr : ', qsr(ji,jj)
683	WRITE(numout,*) ' qns : ', qns(ji,jj)
684	WRITE(numout,*)
685	WRITE(numout,*) ' hfx_mass : ', hfx_thd(ji,jj) + hfx_dyn(ji,jj) + hfx_snw(ji,jj) + hfx_res(ji,jj)
686	WRITE(numout,*) ' qt_atm_oi : ', qt_atm_oi(ji,jj)
687	WRITE(numout,*) ' qt_oce_ai : ', qt_oce_ai(ji,jj)
688	WRITE(numout,*) ' dhc : ', diag_heat(ji,jj)
689	WRITE(numout,*)
690	WRITE(numout,*) ' hfx_dyn : ', hfx_dyn(ji,jj)
691	WRITE(numout,*) ' hfx_thd : ', hfx_thd(ji,jj)
692	WRITE(numout,*) ' hfx_res : ', hfx_res(ji,jj)
693	WRITE(numout,*) ' qsb_ice_bot : ', qsb_ice_bot(ji,jj)
694	WRITE(numout,) ' qlead : ', qlead(ji,jj) r1_rdtice
695	WRITE(numout,*)
696	WRITE(numout,) ' - Salt fluxes at bottom interface **'
697	WRITE(numout,*) ' emp : ', emp (ji,jj)
698	WRITE(numout,*) ' sfx : ', sfx (ji,jj)
699	WRITE(numout,*) ' sfx_res : ', sfx_res(ji,jj)
700	WRITE(numout,*) ' sfx_bri : ', sfx_bri(ji,jj)
701	WRITE(numout,*) ' sfx_dyn : ', sfx_dyn(ji,jj)
702	WRITE(numout,*)
703	WRITE(numout,*) ' - Momentum fluxes '
704	WRITE(numout,*) ' utau : ', utau(ji,jj)
705	WRITE(numout,*) ' vtau : ', vtau(ji,jj)
706	ENDIF
707	WRITE(numout,*) ' '
708	!
709	END DO
710	END DO
711	!
712	END SUBROUTINE ice_prt
713
714	SUBROUTINE ice_prt3D( cd_routine )
715	!!-------------------------------------------------------------------
716	!! * ROUTINE ice_prt3D *
717	!!
718	!! ** Purpose : CTL prints of ice arrays in case ln_ctl is activated
719	!!
720	!!-------------------------------------------------------------------
721	CHARACTER(len=*), INTENT(in) :: cd_routine ! name of the routine
722	INTEGER :: jk, jl ! dummy loop indices
723
724	CALL prt_ctl_info(' ========== ')
725	CALL prt_ctl_info( cd_routine )
726	CALL prt_ctl_info(' ========== ')
727	CALL prt_ctl_info(' - Cell values : ')
728	CALL prt_ctl_info(' ~~~~~~~~~~~~~ ')
729	CALL prt_ctl(tab2d_1=e1e2t , clinfo1=' cell area :')
730	CALL prt_ctl(tab2d_1=at_i , clinfo1=' at_i :')
731	CALL prt_ctl(tab2d_1=ato_i , clinfo1=' ato_i :')
732	CALL prt_ctl(tab2d_1=vt_i , clinfo1=' vt_i :')
733	CALL prt_ctl(tab2d_1=vt_s , clinfo1=' vt_s :')
734	CALL prt_ctl(tab2d_1=divu_i , clinfo1=' divu_i :')
735	CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :')
736	CALL prt_ctl(tab2d_1=stress1_i , clinfo1=' stress1_i :')
737	CALL prt_ctl(tab2d_1=stress2_i , clinfo1=' stress2_i :')
738	CALL prt_ctl(tab2d_1=stress12_i , clinfo1=' stress12_i :')
739	CALL prt_ctl(tab2d_1=strength , clinfo1=' strength :')
740	CALL prt_ctl(tab2d_1=delta_i , clinfo1=' delta_i :')
741	CALL prt_ctl(tab2d_1=u_ice , clinfo1=' u_ice :', tab2d_2=v_ice , clinfo2=' v_ice :')
742
743	DO jl = 1, jpl
744	CALL prt_ctl_info(' ')
745	CALL prt_ctl_info(' - Category : ', ivar1=jl)
746	CALL prt_ctl_info(' ~~~~~~~~~~')
747	CALL prt_ctl(tab2d_1=h_i (:,:,jl) , clinfo1= ' h_i : ')
748	CALL prt_ctl(tab2d_1=h_s (:,:,jl) , clinfo1= ' h_s : ')
749	CALL prt_ctl(tab2d_1=t_su (:,:,jl) , clinfo1= ' t_su : ')
750	CALL prt_ctl(tab2d_1=t_s (:,:,1,jl) , clinfo1= ' t_snow : ')
751	CALL prt_ctl(tab2d_1=s_i (:,:,jl) , clinfo1= ' s_i : ')
752	CALL prt_ctl(tab2d_1=o_i (:,:,jl) , clinfo1= ' o_i : ')
753	CALL prt_ctl(tab2d_1=a_i (:,:,jl) , clinfo1= ' a_i : ')
754	CALL prt_ctl(tab2d_1=v_i (:,:,jl) , clinfo1= ' v_i : ')
755	CALL prt_ctl(tab2d_1=v_s (:,:,jl) , clinfo1= ' v_s : ')
756	CALL prt_ctl(tab2d_1=e_s (:,:,1,jl) , clinfo1= ' e_snow : ')
757	CALL prt_ctl(tab2d_1=sv_i (:,:,jl) , clinfo1= ' sv_i : ')
758	CALL prt_ctl(tab2d_1=oa_i (:,:,jl) , clinfo1= ' oa_i : ')
759
760	DO jk = 1, nlay_i
761	CALL prt_ctl_info(' - Layer : ', ivar1=jk)
762	CALL prt_ctl(tab2d_1=t_i(:,:,jk,jl) , clinfo1= ' t_i : ')
763	CALL prt_ctl(tab2d_1=e_i(:,:,jk,jl) , clinfo1= ' e_i : ')
764	END DO
765	END DO
766
767	CALL prt_ctl_info(' ')
768	CALL prt_ctl_info(' - Stresses : ')
769	CALL prt_ctl_info(' ~~~~~~~~~~ ')
770	CALL prt_ctl(tab2d_1=utau , clinfo1= ' utau : ', tab2d_2=vtau , clinfo2= ' vtau : ')
771	CALL prt_ctl(tab2d_1=utau_ice , clinfo1= ' utau_ice : ', tab2d_2=vtau_ice , clinfo2= ' vtau_ice : ')
772
773	END SUBROUTINE ice_prt3D
774
775
776	SUBROUTINE ice_drift_wri( kt )
777	!!-------------------------------------------------------------------
778	!! * ROUTINE ice_drift_wri *
779	!!
780	!! ** Purpose : conservation of mass, salt and heat
781	!! write the drift in a ascii file at each time step
782	!! and the total run drifts
783	!!-------------------------------------------------------------------
784	INTEGER, INTENT(in) :: kt ! ice time-step index
785	!
786	INTEGER :: ji, jj
787	REAL(wp) :: zdiag_mass, zdiag_salt, zdiag_heat, zdiag_adv_mass, zdiag_adv_salt, zdiag_adv_heat
788	!!REAL(wp), DIMENSION(jpi,jpj) :: zdiag_mass2D, zdiag_salt2D, zdiag_heat2D
789	!!-------------------------------------------------------------------
790	!
791	IF( kt == nit000 .AND. lwp ) THEN
792	WRITE(numout,*)
793	WRITE(numout,*) 'ice_drift_wri: sea-ice drifts'
794	WRITE(numout,*) '~~~~~~~~~~~~~'
795	ENDIF
796	!
797	!clem: the following lines check the ice drift in 2D.
798	! to use this check, uncomment those lines and add the 3 fields in field_def_ice.xml
799	!!! 2D budgets (must be close to 0)
800	!!IF( iom_use('icedrift_mass') .OR. iom_use('icedrift_salt') .OR. iom_use('icedrift_heat') ) THEN
801	!! DO jj = 1, jpj
802	!! DO ji = 1, jpi
803	!! zdiag_mass2D(ji,jj) = wfx_ice(ji,jj) + wfx_snw(ji,jj) + wfx_spr(ji,jj) + wfx_sub(ji,jj) + wfx_pnd(ji,jj) &
804	!! & + diag_vice(ji,jj) + diag_vsnw(ji,jj) + diag_vpnd(ji,jj) - diag_adv_mass(ji,jj)
805	!! zdiag_salt2D(ji,jj) = sfx(ji,jj) + diag_sice(ji,jj) - diag_adv_salt(ji,jj)
806	!! zdiag_heat2D(ji,jj) = qt_oce_ai(ji,jj) - qt_atm_oi(ji,jj) + diag_heat(ji,jj) - diag_adv_heat(ji,jj)
807	!! END DO
808	!! END DO
809	!! !
810	!! ! write outputs
811	!! CALL iom_put( 'icedrift_mass', zdiag_mass2D )
812	!! CALL iom_put( 'icedrift_salt', zdiag_salt2D )
813	!! CALL iom_put( 'icedrift_heat', zdiag_heat2D )
814	!!ENDIF
815
816	! -- mass diag -- !
817	zdiag_mass = glob_sum( 'icectl', ( wfx_ice + wfx_snw + wfx_spr + wfx_sub + wfx_pnd &
818	& + diag_vice + diag_vsnw + diag_vpnd - diag_adv_mass ) * e1e2t ) * rdt_ice
819	zdiag_adv_mass = glob_sum( 'icectl', diag_adv_mass * e1e2t ) * rdt_ice
820
821	! -- salt diag -- !
822	zdiag_salt = glob_sum( 'icectl', ( sfx + diag_sice - diag_adv_salt ) * e1e2t ) * rdt_ice * 1.e-3
823	zdiag_adv_salt = glob_sum( 'icectl', diag_adv_salt * e1e2t ) * rdt_ice * 1.e-3
824
825	! -- heat diag -- !
826	zdiag_heat = glob_sum( 'icectl', ( qt_oce_ai - qt_atm_oi + diag_heat - diag_adv_heat ) * e1e2t )
827	zdiag_adv_heat = glob_sum( 'icectl', diag_adv_heat * e1e2t )
828
829	! ! write out to file
830	IF( lwp ) THEN
831	! check global drift (must be close to 0)
832	WRITE(numicedrift,FMT='(2x,i6,3x,a19,4x,f25.5)') kt, 'mass drift [kg]', zdiag_mass
833	WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift [kg]', zdiag_salt
834	WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift [W] ', zdiag_heat
835	! check drift from advection scheme (can be /=0 with bdy but not sure why)
836	WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'mass drift adv [kg]', zdiag_adv_mass
837	WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'salt drift adv [kg]', zdiag_adv_salt
838	WRITE(numicedrift,FMT='(11x, a19,4x,f25.5)') 'heat drift adv [W] ', zdiag_adv_heat
839	ENDIF
840	! ! drifts
841	rdiag_icemass = rdiag_icemass + zdiag_mass
842	rdiag_icesalt = rdiag_icesalt + zdiag_salt
843	rdiag_iceheat = rdiag_iceheat + zdiag_heat
844	rdiag_adv_icemass = rdiag_adv_icemass + zdiag_adv_mass
845	rdiag_adv_icesalt = rdiag_adv_icesalt + zdiag_adv_salt
846	rdiag_adv_iceheat = rdiag_adv_iceheat + zdiag_adv_heat
847	!
848	! ! output drifts and close ascii file
849	IF( kt == nitend - nn_fsbc + 1 .AND. lwp ) THEN
850	! to ascii file
851	WRITE(numicedrift,) '*****************************************'
852	WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift [kg]', rdiag_icemass
853	WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run mass drift adv [kg]', rdiag_adv_icemass
854	WRITE(numicedrift,) '*****************************************'
855	WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift [kg]', rdiag_icesalt
856	WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run salt drift adv [kg]', rdiag_adv_icesalt
857	WRITE(numicedrift,) '*****************************************'
858	WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift [W] ', rdiag_iceheat
859	WRITE(numicedrift,FMT='(3x,a23,6x,E10.2)') 'Run heat drift adv [W] ', rdiag_adv_iceheat
860	CLOSE( numicedrift )
861	!
862	! to ocean output
863	WRITE(numout,*)
864	WRITE(numout,*) 'ice_drift_wri: ice drifts information for the run '
865	WRITE(numout,*) '~~~~~~~~~~~~~'
866	! check global drift (must be close to 0)
867	WRITE(numout,*) ' sea-ice mass drift [kg] = ', rdiag_icemass
868	WRITE(numout,*) ' sea-ice salt drift [kg] = ', rdiag_icesalt
869	WRITE(numout,*) ' sea-ice heat drift [W] = ', rdiag_iceheat
870	! check drift from advection scheme (can be /=0 with bdy but not sure why)
871	WRITE(numout,*) ' sea-ice mass drift adv [kg] = ', rdiag_adv_icemass
872	WRITE(numout,*) ' sea-ice salt drift adv [kg] = ', rdiag_adv_icesalt
873	WRITE(numout,*) ' sea-ice heat drift adv [W] = ', rdiag_adv_iceheat
874	ENDIF
875	!
876	END SUBROUTINE ice_drift_wri
877
878	SUBROUTINE ice_drift_init
879	!!----------------------------------------------------------------------
880	!! * ROUTINE ice_drift_init *
881	!!
882	!! ** Purpose : create output file, initialise arrays
883	!!----------------------------------------------------------------------
884	!
885	IF( .NOT.ln_icediachk ) RETURN ! exit
886	!
887	IF(lwp) THEN
888	WRITE(numout,*)
889	WRITE(numout,*) 'ice_drift_init: Output ice drifts to ',TRIM(clname), ' file'
890	WRITE(numout,*) '~~~~~~~~~~~~~'
891	WRITE(numout,*)
892	!
893	! create output ascii file
894	CALL ctl_opn( numicedrift, clname, 'UNKNOWN', 'FORMATTED', 'SEQUENTIAL', 1, numout, lwp, narea )
895	WRITE(numicedrift,*) 'Timestep Drifts'
896	WRITE(numicedrift,) '*****************************************'
897	ENDIF
898	!
899	rdiag_icemass = 0._wp
900	rdiag_icesalt = 0._wp
901	rdiag_iceheat = 0._wp
902	rdiag_adv_icemass = 0._wp
903	rdiag_adv_icesalt = 0._wp
904	rdiag_adv_iceheat = 0._wp
905	!
906	END SUBROUTINE ice_drift_init
907
908	#else
909	!!----------------------------------------------------------------------
910	!! Default option Empty Module No SI3 sea-ice model
911	!!----------------------------------------------------------------------
912	#endif
913
914	!!======================================================================
915	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: NEMO/releases/r4.0/r4.0-HEAD/src/ICE/icectl.F90

Download in other formats: