Context Navigation

icectl.F90 @ 15795

Last change on this file since 15795 was 15377, checked in by clem, 3 years ago
reduce drastically the number of global com (from 250 to 50) when the logical ln_icediachk is on
Property svn:keywords set to `Id`
File size: 44.4 KB

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

Download in other formats: