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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diawri.F90 in NEMO/trunk/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/trunk/src/OCE/DIA/diawri.F90 @ 15136

Last change on this file since 15136 was 15136, checked in by smasson, 3 years ago
trunk: SWE passes sette tests in debug
Property svn:keywords set to `Id`
File size: 64.6 KB

Rev	Line
[3]	1	MODULE diawri
	2	!!======================================================================
	3	!! * MODULE diawri *
	4	!! Ocean diagnostics : write ocean output files
	5	!!=====================================================================
[2528]	6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
	7	!! 4.0 ! 1991-11 (G. Madec)
	8	!! ! 1992-06 (M. Imbard) correction restart file
	9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
	10	!! ! 1993-03 (M. Imbard) suppress writibm
	11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
	12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
	13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
	14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
	15	!! - ! 2002-09 (G. Madec) F90: Free form and module
	16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
	17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
	18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
[5836]	19	!! 3.7 ! 2014-01 (G. Madec) remove eddy induced velocity from no-IOM output
	20	!! ! change name of output variables in dia_wri_state
[14053]	21	!! 4.0 ! 2020-10 (A. Nasser, S. Techene) add diagnostic for SWE
[2528]	22	!!----------------------------------------------------------------------
[3]	23
	24	!!----------------------------------------------------------------------
[2528]	25	!! dia_wri : create the standart output files
	26	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
	27	!!----------------------------------------------------------------------
[9019]	28	USE oce ! ocean dynamics and tracers
[12377]	29	USE isf_oce
	30	USE isfcpl
	31	USE abl ! abl variables in case ln_abl = .true.
[9019]	32	USE dom_oce ! ocean space and time domain
	33	USE phycst ! physical constants
	34	USE dianam ! build name of file (routine)
	35	USE diahth ! thermocline diagnostics
	36	USE dynadv , ONLY: ln_dynadv_vec
	37	USE icb_oce ! Icebergs
	38	USE icbdia ! Iceberg budgets
	39	USE ldftra ! lateral physics: eddy diffusivity coef.
	40	USE ldfdyn ! lateral physics: eddy viscosity coef.
	41	USE sbc_oce ! Surface boundary condition: ocean fields
	42	USE sbc_ice ! Surface boundary condition: ice fields
	43	USE sbcssr ! restoring term toward SST/SSS climatology
	44	USE sbcwave ! wave parameters
	45	USE wet_dry ! wetting and drying
	46	USE zdf_oce ! ocean vertical physics
	47	USE zdfdrg ! ocean vertical physics: top/bottom friction
	48	USE zdfmxl ! mixed layer
[14045]	49	USE zdfosm ! mixed layer
[6140]	50	!
[9019]	51	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	52	USE in_out_manager ! I/O manager
	53	USE dia25h ! 25h Mean output
	54	USE iom !
	55	USE ioipsl !
[5463]	56
[9570]	57	#if defined key_si3
[10425]	58	USE ice
[9019]	59	USE icewri
[1482]	60	#endif
[2715]	61	USE lib_mpp ! MPP library
[3294]	62	USE timing ! preformance summary
[12377]	63	USE diu_bulk ! diurnal warm layer
	64	USE diu_coolskin ! Cool skin
[2528]	65
[3]	66	IMPLICIT NONE
	67	PRIVATE
	68
[2528]	69	PUBLIC dia_wri ! routines called by step.F90
	70	PUBLIC dia_wri_state
[2715]	71	PUBLIC dia_wri_alloc ! Called by nemogcm module
[14239]	72	#if ! defined key_xios
[12377]	73	PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.)
	74	#endif
[2528]	75	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	76	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	77	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	78	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	79	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
[12377]	80	INTEGER :: nid_A, nz_A, nh_A, ndim_A, ndim_hA ! grid_ABL file
[2528]	81	INTEGER :: ndex(1) ! ???
[2715]	82	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
[12377]	83	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
[2715]	84	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	85	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	86
	87	!! * Substitutions
[12377]	88	# include "do_loop_substitute.h90"
[13237]	89	# include "domzgr_substitute.h90"
[3]	90	!!----------------------------------------------------------------------
[9598]	91	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	92	!! $Id$
[10068]	93	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	94	!!----------------------------------------------------------------------
	95	CONTAINS
	96
[14239]	97	#if defined key_xios
[3]	98	!!----------------------------------------------------------------------
[14239]	99	!! 'key_xios' use IOM library
[2528]	100	!!----------------------------------------------------------------------
[9652]	101	INTEGER FUNCTION dia_wri_alloc()
	102	!
	103	dia_wri_alloc = 0
	104	!
	105	END FUNCTION dia_wri_alloc
[2715]	106
[9652]	107
[12377]	108	SUBROUTINE dia_wri( kt, Kmm )
[1482]	109	!!---------------------------------------------------------------------
	110	!! * ROUTINE dia_wri *
	111	!!
	112	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	113	!! NETCDF format is used by default
	114	!!
	115	!! ** Method : use iom_put
	116	!!----------------------------------------------------------------------
	117	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	118	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[1756]	119	!!
[9019]	120	INTEGER :: ji, jj, jk ! dummy loop indices
	121	INTEGER :: ikbot ! local integer
	122	REAL(wp):: zztmp , zztmpx ! local scalar
	123	REAL(wp):: zztmp2, zztmpy ! - -
[14053]	124	REAL(wp):: ze3
[15017]	125	REAL(wp), DIMENSION(A2D( 0)) :: z2d ! 2D workspace
	126	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
[1482]	127	!!----------------------------------------------------------------------
	128	!
[9124]	129	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	130	!
[1482]	131	! Output the initial state and forcings
	132	IF( ninist == 1 ) THEN
[12377]	133	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	134	ninist = 0
	135	ENDIF
[3]	136
[15087]	137	! initialize arrays
	138	z2d(:,:) = 0._wp
	139	z3d(:,:,:) = 0._wp
	140
[6351]	141	! Output of initial vertical scale factor
	142	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	143	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	144	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[14053]	145	CALL iom_put("e3f_0", e3f_0(:,:,:) )
[6351]	146	!
[14086]	147	IF ( iom_use("tpt_dep") ) THEN
[15087]	148	DO_3D( 0, 0, 0, 0, 1, jpk )
	149	z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
	150	END_3D
[15017]	151	CALL iom_put( "tpt_dep", z3d )
[14086]	152	ENDIF
	153
[15087]	154	! --- vertical scale factors --- !
[13237]	155	IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t
[15087]	156	DO_3D( 0, 0, 0, 0, 1, jpk )
	157	z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
	158	END_3D
[15017]	159	CALL iom_put( "e3t", z3d )
	160	IF ( iom_use("e3tdef") ) THEN
[15087]	161	DO_3D( 0, 0, 0, 0, 1, jpk )
	162	z3d(ji,jj,jk) = ( ( z3d(ji,jj,jk) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
	163	END_3D
[15017]	164	CALL iom_put( "e3tdef", z3d )
	165	ENDIF
[13237]	166	ENDIF
	167	IF ( iom_use("e3u") ) THEN ! time-varying e3u
[15087]	168	DO_3D( 0, 0, 0, 0, 1, jpk )
	169	z3d(ji,jj,jk) = e3u(ji,jj,jk,Kmm)
	170	END_3D
[15017]	171	CALL iom_put( "e3u" , z3d )
[13237]	172	ENDIF
	173	IF ( iom_use("e3v") ) THEN ! time-varying e3v
[15087]	174	DO_3D( 0, 0, 0, 0, 1, jpk )
	175	z3d(ji,jj,jk) = e3v(ji,jj,jk,Kmm)
	176	END_3D
[15017]	177	CALL iom_put( "e3v" , z3d )
[13237]	178	ENDIF
	179	IF ( iom_use("e3w") ) THEN ! time-varying e3w
[15087]	180	DO_3D( 0, 0, 0, 0, 1, jpk )
	181	z3d(ji,jj,jk) = e3w(ji,jj,jk,Kmm)
	182	END_3D
[15017]	183	CALL iom_put( "e3w" , z3d )
[13237]	184	ENDIF
[14053]	185	IF ( iom_use("e3f") ) THEN ! time-varying e3f caution here at Kaa
[15087]	186	DO_3D( 0, 0, 0, 0, 1, jpk )
	187	z3d(ji,jj,jk) = e3f(ji,jj,jk)
	188	END_3D
[15017]	189	CALL iom_put( "e3f" , z3d )
[14053]	190	ENDIF
[5461]	191
[15017]	192	IF ( iom_use("ssh") ) THEN
	193	IF( ll_wd ) THEN ! sea surface height (brought back to the reference used for wetting and drying)
	194	CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) )
	195	ELSE
	196	CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height
	197	ENDIF
[9023]	198	ENDIF
	199
[14053]	200	IF( iom_use("wetdep") ) CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) ! wet depth
	201
	202	#if defined key_qco
	203	IF( iom_use("ht") ) CALL iom_put( "ht" , ht(:,:) ) ! water column at t-point
	204	IF( iom_use("hu") ) CALL iom_put( "hu" , hu(:,:,Kmm) ) ! water column at u-point
	205	IF( iom_use("hv") ) CALL iom_put( "hv" , hv(:,:,Kmm) ) ! water column at v-point
	206	IF( iom_use("hf") ) CALL iom_put( "hf" , hf_0(:,:)*( 1._wp + r3f(:,:) ) ) ! water column at f-point (caution here at Naa)
	207	#endif
[15087]	208
	209	! --- tracers T&S --- !
[12377]	210	CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature
	211	CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature
[15087]	212
[5107]	213	IF ( iom_use("sbt") ) THEN
[13472]	214	DO_2D( 0, 0, 0, 0 )
[12377]	215	ikbot = mbkt(ji,jj)
	216	z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
	217	END_2D
[5107]	218	CALL iom_put( "sbt", z2d ) ! bottom temperature
	219	ENDIF
	220
[12377]	221	CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
	222	CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity
[5107]	223	IF ( iom_use("sbs") ) THEN
[13472]	224	DO_2D( 0, 0, 0, 0 )
[12377]	225	ikbot = mbkt(ji,jj)
	226	z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
	227	END_2D
[5107]	228	CALL iom_put( "sbs", z2d ) ! bottom salinity
	229	ENDIF
[5463]	230
[14179]	231	IF( .NOT.lk_SWE ) CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
[13089]	232
[15087]	233	! --- momentum --- !
[5463]	234	IF ( iom_use("taubot") ) THEN ! bottom stress
[15017]	235	zztmp = rho0 * 0.25_wp
[7753]	236	z2d(:,:) = 0._wp
[13295]	237	DO_2D( 0, 0, 0, 0 )
[12377]	238	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 &
	239	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 &
	240	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 &
	241	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2
	242	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
	243	!
	244	END_2D
[5463]	245	CALL iom_put( "taubot", z2d )
	246	ENDIF
[5107]	247
[12377]	248	CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current
	249	CALL iom_put( "ssu", uu(:,:,1,Kmm) ) ! surface i-current
[5107]	250	IF ( iom_use("sbu") ) THEN
[13472]	251	DO_2D( 0, 0, 0, 0 )
[12377]	252	ikbot = mbku(ji,jj)
	253	z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
	254	END_2D
[5107]	255	CALL iom_put( "sbu", z2d ) ! bottom i-current
	256	ENDIF
	257
[12377]	258	CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current
	259	CALL iom_put( "ssv", vv(:,:,1,Kmm) ) ! surface j-current
[5107]	260	IF ( iom_use("sbv") ) THEN
[13472]	261	DO_2D( 0, 0, 0, 0 )
[12377]	262	ikbot = mbkv(ji,jj)
	263	z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
	264	END_2D
[5107]	265	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	266	ENDIF
[1482]	267
[14476]	268	! ! vertical velocity
[15087]	269	IF( ln_zad_Aimp ) THEN
	270	IF( iom_use('woce') ) THEN
	271	DO_3D( 0, 0, 0, 0, 1, jpk )
	272	z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
	273	END_3D
	274	CALL iom_put( "woce", z3d ) ! explicit plus implicit parts
	275	ENDIF
	276	ELSE
	277	CALL iom_put( "woce", ww )
[14476]	278	ENDIF
[13237]	279
[5461]	280	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
[14476]	281	! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[15087]	282	IF( ln_zad_Aimp ) THEN
	283	DO_3D( 0, 0, 0, 0, 1, jpk )
	284	z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
	285	END_3D
	286	ELSE
	287	DO_3D( 0, 0, 0, 0, 1, jpk )
	288	z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
	289	END_3D
	290	ENDIF
[5461]	291	CALL iom_put( "w_masstr" , z3d )
[14476]	292	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d * z3d )
[5461]	293	ENDIF
	294
[9019]	295	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	296	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	297	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
[5107]	298
[9019]	299	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	300	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
[6351]	301
[15017]	302	IF ( iom_use("sssgrad") .OR. iom_use("sssgrad2") ) THEN
	303	DO_2D( 0, 0, 0, 0 ) ! sss gradient
	304	zztmp = ts(ji,jj,1,jp_sal,Kmm)
	305	zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
	306	zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
	307	z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	308	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
	309	END_2D
	310	CALL iom_put( "sssgrad2", z2d ) ! square of module of sss gradient
	311	IF ( iom_use("sssgrad") ) THEN
	312	DO_2D( 0, 0, 0, 0 )
	313	z2d(ji,jj) = SQRT( z2d(ji,jj) )
	314	END_2D
	315	CALL iom_put( "sssgrad", z2d ) ! module of sss gradient
	316	ENDIF
[13472]	317	ENDIF
	318
[5107]	319	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[15017]	320	DO_2D( 0, 0, 0, 0 ) ! sst gradient
[12377]	321	zztmp = ts(ji,jj,1,jp_tem,Kmm)
	322	zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
	323	zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
[15017]	324	z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	325	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
[12377]	326	END_2D
[9019]	327	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
[15017]	328	IF ( iom_use("sstgrad") ) THEN
	329	DO_2D( 0, 0, 0, 0 )
	330	z2d(ji,jj) = SQRT( z2d(ji,jj) )
	331	END_2D
	332	CALL iom_put( "sstgrad", z2d ) ! module of sst gradient
	333	ENDIF
[4990]	334	ENDIF
	335
[9019]	336	! heat and salt contents
[4990]	337	IF( iom_use("heatc") ) THEN
[7753]	338	z2d(:,:) = 0._wp
[13472]	339	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	340	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
	341	END_3D
[12489]	342	CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
[4990]	343	ENDIF
	344
	345	IF( iom_use("saltc") ) THEN
[7753]	346	z2d(:,:) = 0._wp
[13472]	347	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	348	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
	349	END_3D
[15017]	350	CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[4990]	351	ENDIF
[4840]	352	!
[13472]	353	IF( iom_use("salt2c") ) THEN
	354	z2d(:,:) = 0._wp
	355	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	356	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
	357	END_3D
[15017]	358	CALL iom_put( "salt2c", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[13472]	359	ENDIF
	360	!
	361	IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
[15087]	362	DO_3D( 0, 0, 0, 0, 1, jpk )
[15017]	363	zztmpx = uu(ji-1,jj ,jk,Kmm) + uu(ji,jj,jk,Kmm)
	364	zztmpy = vv(ji ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
	365	z3d(ji,jj,jk) = 0.25_wp * ( zztmpxzztmpx + zztmpyzztmpy )
[12377]	366	END_3D
[13472]	367	CALL iom_put( "ke", z3d ) ! kinetic energy
	368
	369	z2d(:,:) = 0._wp
	370	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	371	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
	372	END_3D
[15017]	373	CALL iom_put( "ke_int", z2d ) ! vertically integrated kinetic energy
[4990]	374	ENDIF
[6351]	375	!
[15017]	376	IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
[14053]	377	z2d(:,:) = 0._wp
	378	DO_2D( 0, 0, 0, 0 )
	379	z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
	380	& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
	381	& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
	382	& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
	383	& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
	384	END_2D
	385	IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d )
	386	ENDIF
	387	!
[15017]	388	IF ( iom_use("ssKEf") ) THEN ! surface kinetic energy at F point
	389	z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
[14053]	390	DO_2D( 0, 0, 0, 0 )
	391	z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
	392	& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
	393	& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
	394	& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
	395	& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
	396	END_2D
[14143]	397	CALL iom_put( "ssKEf", z2d )
[14053]	398	ENDIF
	399	!
[15017]	400	CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
[6351]	401	!
[7646]	402	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[15017]	403
[15087]	404	DO_3D( 0, 0, 0, 0, 1, jpk )
	405	z3d(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
	406	END_3D
[15017]	407	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	408
	409	IF( iom_use("u_masstr_vint") ) THEN
	410	z2d(:,:) = 0._wp
	411	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[15087]	412	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk)
[15017]	413	END_3D
	414	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
	415	ENDIF
	416	IF( iom_use("u_heattr") ) THEN
	417	z2d(:,:) = 0._wp
	418	zztmp = 0.5_wp * rcp
	419	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	420	z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
	421	END_3D
	422	CALL iom_put( "u_heattr", z2d ) ! heat transport in i-direction
	423	ENDIF
	424	IF( iom_use("u_salttr") ) THEN
	425	z2d(:,:) = 0._wp
	426	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	427	z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
	428	END_3D
	429	CALL iom_put( "u_salttr", z2d ) ! heat transport in i-direction
	430	ENDIF
	431
[4990]	432	ENDIF
	433
	434	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[15017]	435
[15087]	436	DO_3D( 0, 0, 0, 0, 1, jpk )
	437	z3d(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
	438	END_3D
[15017]	439	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
	440
	441	IF( iom_use("v_heattr") ) THEN
	442	z2d(:,:) = 0._wp
	443	zztmp = 0.5_wp * rcp
	444	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	445	z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
	446	END_3D
	447	CALL iom_put( "v_heattr", z2d ) ! heat transport in j-direction
	448	ENDIF
	449	IF( iom_use("v_salttr") ) THEN
	450	z2d(:,:) = 0._wp
	451	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	452	z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
	453	END_3D
	454	CALL iom_put( "v_salttr", z2d ) ! heat transport in j-direction
	455	ENDIF
[4761]	456
[1756]	457	ENDIF
[7646]	458
	459	IF( iom_use("tosmint") ) THEN
[9019]	460	z2d(:,:) = 0._wp
[13295]	461	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[15017]	462	z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
[12377]	463	END_3D
[15017]	464	CALL iom_put( "tosmint", z2d ) ! Vertical integral of temperature
[7646]	465	ENDIF
	466	IF( iom_use("somint") ) THEN
[15017]	467	z2d(:,:) = 0._wp
[13295]	468	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[15017]	469	z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
[12377]	470	END_3D
[15017]	471	CALL iom_put( "somint", z2d ) ! Vertical integral of salinity
[7646]	472	ENDIF
	473
[15017]	474	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
[12377]	475
[15017]	476	IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging
[14053]	477
[14143]	478	! Output of surface vorticity terms
[15017]	479	!
	480	CALL iom_put( "ssplavor", ff_f ) ! planetary vorticity ( f )
	481	!
	482	IF ( iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
	483	& iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
[14053]	484	!
	485	z2d(:,:) = 0._wp
[15017]	486	DO_2D( 0, 0, 0, 0 )
[14053]	487	z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
	488	& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
	489	END_2D
[15017]	490	CALL iom_put( "ssrelvor", z2d ) ! relative vorticity ( zeta )
[14053]	491	!
[15017]	492	IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
	493	DO_2D( 0, 0, 0, 0 )
	494	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	495	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	496	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	497	ELSE ; ze3 = 0._wp
	498	ENDIF
	499	z2d(ji,jj) = ze3 * z2d(ji,jj)
	500	END_2D
	501	CALL iom_put( "ssrelpotvor", z2d ) ! relative potential vorticity (zeta/h)
	502	!
	503	IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
	504	DO_2D( 0, 0, 0, 0 )
	505	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	506	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	507	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	508	ELSE ; ze3 = 0._wp
	509	ENDIF
	510	z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)
	511	END_2D
	512	CALL iom_put( "ssabspotvor", z2d ) ! absolute potential vorticity ( q )
	513	!
	514	IF ( iom_use("ssEns") ) THEN
	515	DO_2D( 0, 0, 0, 0 )
	516	z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj)
	517	END_2D
	518	CALL iom_put( "ssEns", z2d ) ! potential enstrophy ( 1/2*q2 )
	519	ENDIF
[14053]	520	ENDIF
[15017]	521	ENDIF
[14053]	522	ENDIF
[6140]	523
[9124]	524	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	525	!
[1482]	526	END SUBROUTINE dia_wri
	527
	528	#else
[2528]	529	!!----------------------------------------------------------------------
	530	!! Default option use IOIPSL library
	531	!!----------------------------------------------------------------------
	532
[9652]	533	INTEGER FUNCTION dia_wri_alloc()
	534	!!----------------------------------------------------------------------
	535	INTEGER, DIMENSION(2) :: ierr
	536	!!----------------------------------------------------------------------
[12493]	537	IF( nn_write == -1 ) THEN
	538	dia_wri_alloc = 0
	539	ELSE
	540	ierr = 0
	541	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	542	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	543	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
[9652]	544	!
[12493]	545	dia_wri_alloc = MAXVAL(ierr)
	546	CALL mpp_sum( 'diawri', dia_wri_alloc )
	547	!
	548	ENDIF
[9652]	549	!
	550	END FUNCTION dia_wri_alloc
[12377]	551
	552	INTEGER FUNCTION dia_wri_alloc_abl()
	553	!!----------------------------------------------------------------------
	554	ALLOCATE( ndex_hA(jpijpj), ndex_A (jpijpj*jpkam1), STAT=dia_wri_alloc_abl)
	555	CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
	556	!
	557	END FUNCTION dia_wri_alloc_abl
[9652]	558
	559
[12377]	560	SUBROUTINE dia_wri( kt, Kmm )
[3]	561	!!---------------------------------------------------------------------
	562	!! * ROUTINE dia_wri *
	563	!!
	564	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	565	!! NETCDF format is used by default
	566	!!
	567	!! ** Method : At the beginning of the first time step (nit000),
	568	!! define all the NETCDF files and fields
	569	!! At each time step call histdef to compute the mean if ncessary
[11536]	570	!! Each nn_write time step, output the instantaneous or mean fields
[3]	571	!!----------------------------------------------------------------------
[5836]	572	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	573	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[5836]	574	!
[2528]	575	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	576	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	577	INTEGER :: inum = 11 ! temporary logical unit
[3294]	578	INTEGER :: ji, jj, jk ! dummy loop indices
	579	INTEGER :: ierr ! error code return from allocation
[2528]	580	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[12377]	581	INTEGER :: ipka ! ABL
[3609]	582	INTEGER :: jn, ierror ! local integers
[6140]	583	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	584	!
[9019]	585	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
[13237]	586	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept ! 3D workspace
[12377]	587	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
[3]	588	!!----------------------------------------------------------------------
[1482]	589	!
[9019]	590	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[12377]	591	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	592	ninist = 0
	593	ENDIF
	594	!
[11536]	595	IF( nn_write == -1 ) RETURN ! we will never do any output
	596	!
	597	IF( ln_timing ) CALL timing_start('dia_wri')
	598	!
[3]	599	! 0. Initialisation
	600	! -----------------
[632]	601
[9019]	602	ll_print = .FALSE. ! local variable for debugging
[3]	603	ll_print = ll_print .AND. lwp
	604
	605	! Define frequency of output and means
[5566]	606	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	607	#if defined key_diainstant
[12489]	608	zsto = nn_write * rn_Dt
[1312]	609	clop = "inst("//TRIM(clop)//")"
[3]	610	#else
[12489]	611	zsto=rn_Dt
[1312]	612	clop = "ave("//TRIM(clop)//")"
[3]	613	#endif
[12489]	614	zout = nn_write * rn_Dt
	615	zmax = ( nitend - nit000 + 1 ) * rn_Dt
[3]	616
	617	! Define indices of the horizontal output zoom and vertical limit storage
[13286]	618	iimi = Nis0 ; iima = Nie0
	619	ijmi = Njs0 ; ijma = Nje0
[3]	620	ipk = jpk
[12377]	621	IF(ln_abl) ipka = jpkam1
[3]	622
	623	! define time axis
[1334]	624	it = kt
	625	itmod = kt - nit000 + 1
[3]	626
[13237]	627	! store e3t for subsitute
	628	DO jk = 1, jpk
	629	ze3t (:,:,jk) = e3t (:,:,jk,Kmm)
	630	zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
	631	END DO
[3]	632
[13237]	633
[3]	634	! 1. Define NETCDF files and fields at beginning of first time step
	635	! -----------------------------------------------------------------
	636
	637	IF( kt == nit000 ) THEN
	638
	639	! Define the NETCDF files (one per grid)
[632]	640
[3]	641	! Compute julian date from starting date of the run
[12489]	642	CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
[1309]	643	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	644	IF(lwp)WRITE(numout,*)
	645	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	646	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	647	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	648	' limit storage in depth = ', ipk
	649
	650	! WRITE root name in date.file for use by postpro
[1581]	651	IF(lwp) THEN
[11536]	652	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	653	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	654	WRITE(inum,*) clhstnam
	655	CLOSE(inum)
	656	ENDIF
[632]	657
[3]	658	! Define the T grid FILE ( nid_T )
[632]	659
[11536]	660	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	661	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	662	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	663	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	664	& nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	665	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	666	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	667	! ! Index of ocean points
	668	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	669	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	670	!
	671	IF( ln_icebergs ) THEN
	672	!
	673	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	674	!! that routine is called from nemogcm, so do it here immediately before its needed
	675	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	676	CALL mpp_sum( 'diawri', ierror )
[3609]	677	IF( ierror /= 0 ) THEN
	678	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	679	RETURN
	680	ENDIF
	681	!
	682	!! iceberg vertical coordinate is class number
	683	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	684	& "number", nclasses, class_num, nb_T )
	685	!
	686	!! each class just needs the surface index pattern
	687	ndim_bT = 3
	688	DO jn = 1,nclasses
	689	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	690	ENDDO
	691	!
	692	ENDIF
[3]	693
	694	! Define the U grid FILE ( nid_U )
	695
[11536]	696	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	697	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	698	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	699	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	700	& nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	701	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	702	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	703	! ! Index of ocean points
	704	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	705	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	706
	707	! Define the V grid FILE ( nid_V )
	708
[11536]	709	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	710	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	711	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	712	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	713	& nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	714	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	715	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	716	! ! Index of ocean points
	717	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	718	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	719
	720	! Define the W grid FILE ( nid_W )
	721
[11536]	722	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	723	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	724	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	725	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	726	& nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	727	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	728	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	729
[12377]	730	IF( ln_abl ) THEN
	731	! Define the ABL grid FILE ( nid_A )
	732	CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
	733	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	734	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	735	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	736	& nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
[12377]	737	CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept
	738	& "m", ipka, ght_abl(2:jpka), nz_A, "up" )
	739	! ! Index of ocean points
	740	ALLOCATE( zw3d_abl(jpi,jpj,ipka) )
	741	zw3d_abl(:,:,:) = 1._wp
	742	CALL wheneq( jpijpjipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume
	743	CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface
	744	DEALLOCATE(zw3d_abl)
	745	ENDIF
[13237]	746	!
[632]	747
[3]	748	! Declare all the output fields as NETCDF variables
	749
	750	! !!! nid_T : 3D
	751	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	752	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	753	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	754	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	755	IF( .NOT.ln_linssh ) THEN
[13237]	756	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t n
[4292]	757	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13237]	758	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t n
[4292]	759	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13237]	760	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t n
[4292]	761	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	762	ENDIF
[3]	763	! !!! nid_T : 2D
	764	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	765	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	766	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	767	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	768	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	769	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	770	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	771	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	772	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	773	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	774	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	775	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	776	IF( ln_linssh ) THEN
[12377]	777	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts(:,:,1,jp_tem,Kmm)
[3625]	778	& , "KgC/m2/s", & ! sosst_cd
[4292]	779	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	780	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts(:,:,1,jp_sal,Kmm)
[3625]	781	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	782	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	783	ENDIF
[888]	784	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	785	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	786	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	787	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[15136]	788	IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
	789	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	790	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	791	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	792	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	793	ENDIF
[1037]	794	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	795	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	796	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	797	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	798	!
	799	IF( ln_abl ) THEN
	800	CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl
	801	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	802	CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl
	803	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	804	CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl
	805	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	806	CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl
	807	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	808	CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl
	809	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	810	CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl
	811	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	812	CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl
	813	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	814	CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh
	815	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	816	#if defined key_si3
	817	CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i
	818	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	819	#endif
	820	CALL histend( nid_A, snc4chunks=snc4set )
	821	ENDIF
	822	!
[3609]	823	IF( ln_icebergs ) THEN
	824	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	825	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	826	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	827	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	828	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	829	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	830	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	831	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	832	IF( ln_bergdia ) THEN
	833	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	834	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	835	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	836	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	837	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	838	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	839	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	840	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	841	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	842	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	843	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	844	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	845	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	846	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	847	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	848	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	849	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	850	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	851	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	852	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	853	ENDIF
	854	ENDIF
	855
[11536]	856	IF( ln_ssr ) THEN
[4990]	857	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	858	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	859	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	860	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	861	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	862	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	863	ENDIF
[11536]	864
[3]	865	clmx ="l_max(only(x))" ! max index on a period
[5836]	866	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	867	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	868	#if defined key_diahth
	869	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	870	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	871	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	872	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	873	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	874	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	875	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	876	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	877	#endif
	878
[2528]	879	CALL histend( nid_T, snc4chunks=snc4set )
[3]	880
	881	! !!! nid_U : 3D
[12377]	882	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! uu(:,:,:,Kmm)
[3]	883	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	884	IF( ln_wave .AND. ln_sdw) THEN
	885	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	886	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	887	ENDIF
[3]	888	! !!! nid_U : 2D
[888]	889	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	890	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	891
[2528]	892	CALL histend( nid_U, snc4chunks=snc4set )
[3]	893
	894	! !!! nid_V : 3D
[12377]	895	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vv(:,:,:,Kmm)
[3]	896	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	897	IF( ln_wave .AND. ln_sdw) THEN
	898	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	899	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	900	ENDIF
[3]	901	! !!! nid_V : 2D
[888]	902	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	903	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	904
[2528]	905	CALL histend( nid_V, snc4chunks=snc4set )
[3]	906
	907	! !!! nid_W : 3D
[12377]	908	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! ww
[3]	909	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	910	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	911	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	912	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	913	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	914
[9019]	915	IF( ln_zdfddm ) THEN
[3]	916	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	917	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	918	ENDIF
[7646]	919
	920	IF( ln_wave .AND. ln_sdw) THEN
	921	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	922	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	923	ENDIF
[3]	924	! !!! nid_W : 2D
[2528]	925	CALL histend( nid_W, snc4chunks=snc4set )
[3]	926
	927	IF(lwp) WRITE(numout,*)
	928	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	929	IF(ll_print) CALL FLUSH(numout )
	930
	931	ENDIF
	932
	933	! 2. Start writing data
	934	! ---------------------
	935
[4292]	936	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	937	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	938	! donne le nombre d'elements, et ndex la liste des indices a sortir
	939
[11536]	940	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[3]	941	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	942	WRITE(numout,*) '~~~~~~ '
	943	ENDIF
	944
[6140]	945	IF( .NOT.ln_linssh ) THEN
[13237]	946	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! heat content
	947	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! salt content
	948	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	949	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	950	ELSE
[12377]	951	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
	952	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
	953	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT ) ! sea surface temperature
	954	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity
[4292]	955	ENDIF
[6140]	956	IF( .NOT.ln_linssh ) THEN
[13237]	957	zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	958	CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:) , ndim_T , ndex_T ) ! level thickness
	959	CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T ) ! t-point depth
[4292]	960	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	961	ENDIF
[12377]	962	CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
[2528]	963	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	964	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	965	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	966	! (includes virtual salt flux beneath ice
	967	! in linear free surface case)
[6140]	968	IF( ln_linssh ) THEN
[12377]	969	zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
[4292]	970	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[12377]	971	zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
[4292]	972	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	973	ENDIF
[888]	974	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	975	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[15136]	976	IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
	977	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
	978	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
	979	ENDIF
[1037]	980	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	981	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[12377]	982	!
	983	IF( ln_abl ) THEN
	984	ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
	985	IF( ln_mskland ) THEN
	986	DO jk=1,jpka
	987	zw3d_abl(:,:,jk) = tmask(:,:,1)
	988	END DO
	989	ELSE
	990	zw3d_abl(:,:,:) = 1._wp
	991	ENDIF
	992	CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh
	993	CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl
	994	CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl
	995	CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl
	996	CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl
	997	CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl
	998	CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl
	999	CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl
	1000	#if defined key_si3
	1001	CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i
	1002	#endif
	1003	DEALLOCATE(zw3d_abl)
	1004	ENDIF
	1005	!
[3609]	1006	IF( ln_icebergs ) THEN
	1007	!
	1008	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	1009	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	1010	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	1011	!
	1012	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	1013	!
	1014	IF( ln_bergdia ) THEN
	1015	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	1016	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	1017	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	1018	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	1019	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	1020	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	1021	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	1022	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	1023	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	1024	!
	1025	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	1026	ENDIF
	1027	ENDIF
	1028
[11536]	1029	IF( ln_ssr ) THEN
[4990]	1030	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	1031	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[12377]	1032	zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
[4990]	1033	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	1034	ENDIF
	1035	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	1036	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	1037
	1038	#if defined key_diahth
	1039	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	1040	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	1041	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	1042	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	1043	#endif
[888]	1044
[12377]	1045	CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
[888]	1046	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	1047
[12377]	1048	CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
[888]	1049	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	1050
[11418]	1051	IF( ln_zad_Aimp ) THEN
[12377]	1052	CALL histwrite( nid_W, "vovecrtz", it, ww + wi , ndim_T, ndex_T ) ! vert. current
[11418]	1053	ELSE
[12377]	1054	CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current
[11418]	1055	ENDIF
[3]	1056	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	1057	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	1058	IF( ln_zdfddm ) THEN
	1059	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	1060	ENDIF
	1061
[7646]	1062	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	1063	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	1064	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	1065	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	1066	ENDIF
	1067
[1318]	1068	! 3. Close all files
	1069	! ---------------------------------------
[1561]	1070	IF( kt == nitend ) THEN
[3]	1071	CALL histclo( nid_T )
	1072	CALL histclo( nid_U )
	1073	CALL histclo( nid_V )
	1074	CALL histclo( nid_W )
[12377]	1075	IF(ln_abl) CALL histclo( nid_A )
[3]	1076	ENDIF
[2528]	1077	!
[9124]	1078	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	1079	!
[3]	1080	END SUBROUTINE dia_wri
[1567]	1081	#endif
	1082
[12377]	1083	SUBROUTINE dia_wri_state( Kmm, cdfile_name )
[3]	1084	!!---------------------------------------------------------------------
	1085	!! * ROUTINE dia_wri_state *
	1086	!!
	1087	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	1088	!! the instantaneous ocean state and forcing fields.
	1089	!! Used to find errors in the initial state or save the last
	1090	!! ocean state in case of abnormal end of a simulation
	1091	!!
	1092	!! ** Method : NetCDF files using ioipsl
	1093	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	1094	!! File 'output.abort.nc' is created in case of abnormal job end
	1095	!!----------------------------------------------------------------------
[12377]	1096	INTEGER , INTENT( in ) :: Kmm ! time level index
[1334]	1097	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	1098	!!
[12377]	1099	INTEGER :: inum, jk
[14053]	1100	REAL(wp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept ! 3D workspace for qco substitution
[3]	1101	!!----------------------------------------------------------------------
[3294]	1102	!
[13237]	1103	IF(lwp) THEN
	1104	WRITE(numout,*)
	1105	WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	1106	WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
	1107	WRITE(numout,*) ' and named :', cdfile_name, '...nc'
	1108	ENDIF
[12649]	1109	!
[13237]	1110	DO jk = 1, jpk
	1111	ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
	1112	zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
	1113	END DO
	1114	!
[12649]	1115	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
	1116	!
[12377]	1117	CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature
	1118	CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity
	1119	CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm) ) ! sea surface height
	1120	CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm) ) ! now i-velocity
	1121	CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity
[11418]	1122	IF( ln_zad_Aimp ) THEN
[12377]	1123	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi ) ! now k-velocity
[11418]	1124	ELSE
[12377]	1125	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity
[11418]	1126	ENDIF
[12377]	1127	CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity
[13237]	1128	CALL iom_rstput( 0, 0, inum, 'ht' , ht(:,:) ) ! now water column height
[12649]	1129	!
[12377]	1130	IF ( ln_isf ) THEN
	1131	IF (ln_isfcav_mlt) THEN
	1132	CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav ) ! now k-velocity
	1133	CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity
	1134	CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity
[12649]	1135	CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) ) ! now k-velocity
	1136	CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) ) ! now k-velocity
	1137	CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
[12377]	1138	END IF
	1139	IF (ln_isfpar_mlt) THEN
[12649]	1140	CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) ) ! now k-velocity
[12377]	1141	CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity
	1142	CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity
	1143	CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity
[12649]	1144	CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) ) ! now k-velocity
	1145	CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) ) ! now k-velocity
	1146	CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
[12377]	1147	END IF
	1148	END IF
[12649]	1149	!
[7646]	1150	IF( ALLOCATED(ahtu) ) THEN
[10425]	1151	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	1152	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	1153	ENDIF
	1154	IF( ALLOCATED(ahmt) ) THEN
[10425]	1155	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	1156	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	1157	ENDIF
[10425]	1158	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	1159	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	1160	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	1161	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	1162	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	1163	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	1164	IF( .NOT.ln_linssh ) THEN
[13237]	1165	CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept ) ! T-cell depth
	1166	CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t ) ! T-cell thickness
[10425]	1167	END IF
[7646]	1168	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	1169	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	1170	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	1171	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	1172	ENDIF
[12377]	1173	IF ( ln_abl ) THEN
	1174	CALL iom_rstput ( 0, 0, inum, "uz1_abl", u_abl(:,:,2,nt_a ) ) ! now first level i-wind
	1175	CALL iom_rstput ( 0, 0, inum, "vz1_abl", v_abl(:,:,2,nt_a ) ) ! now first level j-wind
	1176	CALL iom_rstput ( 0, 0, inum, "tz1_abl", tq_abl(:,:,2,nt_a,1) ) ! now first level temperature
	1177	CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity
	1178	ENDIF
[14045]	1179	IF( ln_zdfosm ) THEN
	1180	CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth
	1181	CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth
	1182	CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion
	1183	CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity
	1184	CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing
	1185	CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing
	1186	CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask ) ! non-local u forcing
	1187	CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask ) ! non-local v forcing
	1188	IF( ln_osm_mle ) THEN
	1189	CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth
	1190	END IF
	1191	ENDIF
[12649]	1192	!
	1193	CALL iom_close( inum )
	1194	!
[10425]	1195	#if defined key_si3
	1196	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
[12649]	1197	CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' )
[10425]	1198	CALL ice_wri_state( inum )
[12649]	1199	CALL iom_close( inum )
[1561]	1200	ENDIF
[12933]	1201	!
[1561]	1202	#endif
[3]	1203	END SUBROUTINE dia_wri_state
[6140]	1204
[3]	1205	!!======================================================================
	1206	END MODULE diawri

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

Download in other formats: