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 branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2016/dev_v3_6_STABLE_r6506_AGRIF_LIM3/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8142

Last change on this file since 8142 was 7158, checked in by clem, 8 years ago
debug branch
Property svn:keywords set to `Id`
File size: 57.1 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
	19	!!----------------------------------------------------------------------
[3]	20
	21	!!----------------------------------------------------------------------
[2528]	22	!! dia_wri : create the standart output files
	23	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
	24	!!----------------------------------------------------------------------
[3]	25	USE oce ! ocean dynamics and tracers
	26	USE dom_oce ! ocean space and time domain
[4292]	27	USE dynadv, ONLY: ln_dynadv_vec
[3]	28	USE zdf_oce ! ocean vertical physics
	29	USE ldftra_oce ! ocean active tracers: lateral physics
	30	USE ldfdyn_oce ! ocean dynamics: lateral physics
[3294]	31	USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
[3]	32	USE sol_oce ! solver variables
[888]	33	USE sbc_oce ! Surface boundary condition: ocean fields
	34	USE sbc_ice ! Surface boundary condition: ice fields
[3609]	35	USE icb_oce ! Icebergs
	36	USE icbdia ! Iceberg budgets
[888]	37	USE sbcssr ! restoring term toward SST/SSS climatology
[3]	38	USE phycst ! physical constants
	39	USE zdfmxl ! mixed layer
	40	USE dianam ! build name of file (routine)
	41	USE zdfddm ! vertical physics: double diffusion
	42	USE diahth ! thermocline diagnostics
	43	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
	44	USE in_out_manager ! I/O manager
[216]	45	USE diadimg ! dimg direct access file format output
[1359]	46	USE iom
[460]	47	USE ioipsl
[5463]	48	USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities
	49
[1482]	50	#if defined key_lim2
	51	USE limwri_2
[4161]	52	#elif defined key_lim3
	53	USE limwri
[1482]	54	#endif
[2715]	55	USE lib_mpp ! MPP library
[3294]	56	USE timing ! preformance summary
	57	USE wrk_nemo ! working array
[2528]	58
[3]	59	IMPLICIT NONE
	60	PRIVATE
	61
[2528]	62	PUBLIC dia_wri ! routines called by step.F90
	63	PUBLIC dia_wri_state
[2715]	64	PUBLIC dia_wri_alloc ! Called by nemogcm module
[3]	65
[2528]	66	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	67	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	68	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	69	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	70	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
	71	INTEGER :: ndex(1) ! ???
[2715]	72	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
	73	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	74	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	75
	76	!! * Substitutions
	77	# include "zdfddm_substitute.h90"
[1756]	78	# include "domzgr_substitute.h90"
	79	# include "vectopt_loop_substitute.h90"
[3]	80	!!----------------------------------------------------------------------
[2528]	81	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
[5217]	82	!! $Id$
[2528]	83	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
[3]	84	!!----------------------------------------------------------------------
	85	CONTAINS
	86
[2715]	87	INTEGER FUNCTION dia_wri_alloc()
	88	!!----------------------------------------------------------------------
	89	INTEGER, DIMENSION(2) :: ierr
	90	!!----------------------------------------------------------------------
	91	ierr = 0
	92	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	93	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	94	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
	95	!
	96	dia_wri_alloc = MAXVAL(ierr)
	97	IF( lk_mpp ) CALL mpp_sum( dia_wri_alloc )
	98	!
	99	END FUNCTION dia_wri_alloc
	100
[216]	101	#if defined key_dimgout
[3]	102	!!----------------------------------------------------------------------
[2528]	103	!! 'key_dimgout' DIMG output file
[3]	104	!!----------------------------------------------------------------------
	105	# include "diawri_dimg.h90"
	106
	107	#else
	108	!!----------------------------------------------------------------------
	109	!! Default option NetCDF output file
	110	!!----------------------------------------------------------------------
[2528]	111	# if defined key_iomput
[3]	112	!!----------------------------------------------------------------------
[2528]	113	!! 'key_iomput' use IOM library
	114	!!----------------------------------------------------------------------
[2715]	115
[1561]	116	SUBROUTINE dia_wri( kt )
[1482]	117	!!---------------------------------------------------------------------
	118	!! * ROUTINE dia_wri *
	119	!!
	120	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	121	!! NETCDF format is used by default
	122	!!
	123	!! ** Method : use iom_put
	124	!!----------------------------------------------------------------------
[1756]	125	!!
[1482]	126	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[1756]	127	!!
	128	INTEGER :: ji, jj, jk ! dummy loop indices
[5460]	129	INTEGER :: jkbot !
[1756]	130	REAL(wp) :: zztmp, zztmpx, zztmpy !
[3294]	131	!!
[4761]	132	REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace
[3294]	133	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace
[1482]	134	!!----------------------------------------------------------------------
	135	!
[3294]	136	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
	137	!
[4990]	138	CALL wrk_alloc( jpi , jpj , z2d )
[3294]	139	CALL wrk_alloc( jpi , jpj, jpk , z3d )
[2715]	140	!
[1482]	141	! Output the initial state and forcings
	142	IF( ninist == 1 ) THEN
	143	CALL dia_wri_state( 'output.init', kt )
	144	ninist = 0
	145	ENDIF
[3]	146
[6315]	147	! Output of initial vertical scale factor
	148	CALL iom_put("e3t_0", e3t_0(:,:,:) )
	149	CALL iom_put("e3u_0", e3t_0(:,:,:) )
	150	CALL iom_put("e3v_0", e3t_0(:,:,:) )
	151	!
[6348]	152	CALL iom_put( "e3t" , fse3t_n(:,:,:) )
	153	CALL iom_put( "e3u" , fse3u_n(:,:,:) )
	154	CALL iom_put( "e3v" , fse3v_n(:,:,:) )
	155	CALL iom_put( "e3w" , fse3w_n(:,:,:) )
	156	IF( iom_use("e3tdef") ) &
	157	CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
[5461]	158
[6348]	159
[5461]	160	CALL iom_put( "ssh" , sshn ) ! sea surface height
[5107]	161
	162	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
	163	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
	164	IF ( iom_use("sbt") ) THEN
[4990]	165	DO jj = 1, jpj
	166	DO ji = 1, jpi
[5460]	167	jkbot = mbkt(ji,jj)
	168	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
[4990]	169	END DO
[5107]	170	END DO
	171	CALL iom_put( "sbt", z2d ) ! bottom temperature
	172	ENDIF
	173
	174	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
	175	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
	176	IF ( iom_use("sbs") ) THEN
[4990]	177	DO jj = 1, jpj
	178	DO ji = 1, jpi
[5460]	179	jkbot = mbkt(ji,jj)
	180	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
[4990]	181	END DO
[5107]	182	END DO
	183	CALL iom_put( "sbs", z2d ) ! bottom salinity
	184	ENDIF
[5463]	185
	186	IF ( iom_use("taubot") ) THEN ! bottom stress
	187	z2d(:,:) = 0._wp
	188	DO jj = 2, jpjm1
	189	DO ji = fs_2, fs_jpim1 ! vector opt.
	190	zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) &
	191	& + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) )
	192	zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) &
	193	& + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) )
	194	z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)
	195	!
	196	ENDDO
	197	ENDDO
	198	CALL lbc_lnk( z2d, 'T', 1. )
	199	CALL iom_put( "taubot", z2d )
	200	ENDIF
[5107]	201
	202	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
	203	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
	204	IF ( iom_use("sbu") ) THEN
[4990]	205	DO jj = 1, jpj
	206	DO ji = 1, jpi
[5460]	207	jkbot = mbku(ji,jj)
	208	z2d(ji,jj) = un(ji,jj,jkbot)
[4990]	209	END DO
[5107]	210	END DO
	211	CALL iom_put( "sbu", z2d ) ! bottom i-current
	212	ENDIF
[5463]	213	#if defined key_dynspg_ts
	214	CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current
	215	#else
	216	CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current
	217	#endif
[5107]	218
	219	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
	220	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
	221	IF ( iom_use("sbv") ) THEN
[4990]	222	DO jj = 1, jpj
	223	DO ji = 1, jpi
[5460]	224	jkbot = mbkv(ji,jj)
	225	z2d(ji,jj) = vn(ji,jj,jkbot)
[4990]	226	END DO
[5107]	227	END DO
	228	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	229	ENDIF
[5463]	230	#if defined key_dynspg_ts
	231	CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current
	232	#else
	233	CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current
	234	#endif
[1482]	235
[5461]	236	CALL iom_put( "woce", wn ) ! vertical velocity
	237	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
	238	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
	239	z2d(:,:) = rau0 * e12t(:,:)
	240	DO jk = 1, jpk
	241	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
	242	END DO
	243	CALL iom_put( "w_masstr" , z3d )
	244	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
	245	ENDIF
	246
[5107]	247	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	248	CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef.
	249	CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef. (useful only with key_zdfddm)
[6348]	250	! Log of eddy diff coef
	251	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt (:,:,:) ) ) )
	252	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) )
[5107]	253
	254	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[4990]	255	DO jj = 2, jpjm1 ! sst gradient
	256	DO ji = fs_2, fs_jpim1 ! vector opt.
	257	zztmp = tsn(ji,jj,1,jp_tem)
	258	zztmpx = ( tsn(ji+1,jj ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) / e1u(ji-1,jj )
	259	zztmpy = ( tsn(ji ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) / e2v(ji ,jj-1)
	260	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	261	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	262	END DO
[1756]	263	END DO
[4990]	264	CALL lbc_lnk( z2d, 'T', 1. )
	265	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
	266	z2d(:,:) = SQRT( z2d(:,:) )
	267	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
	268	ENDIF
	269
[4761]	270	! clem: heat and salt content
[4990]	271	IF( iom_use("heatc") ) THEN
	272	z2d(:,:) = 0._wp
	273	DO jk = 1, jpkm1
[5107]	274	DO jj = 1, jpj
	275	DO ji = 1, jpi
[4990]	276	z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
	277	END DO
[4761]	278	END DO
	279	END DO
[4990]	280	CALL iom_put( "heatc", (rau0 * rcp) * z2d ) ! vertically integrated heat content (J/m2)
	281	ENDIF
	282
	283	IF( iom_use("saltc") ) THEN
	284	z2d(:,:) = 0._wp
	285	DO jk = 1, jpkm1
[5107]	286	DO jj = 1, jpj
	287	DO ji = 1, jpi
[4990]	288	z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
	289	END DO
	290	END DO
	291	END DO
	292	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
	293	ENDIF
[4840]	294	!
[4990]	295	IF ( iom_use("eken") ) THEN
	296	rke(:,:,jk) = 0._wp ! kinetic energy
	297	DO jk = 1, jpkm1
	298	DO jj = 2, jpjm1
	299	DO ji = fs_2, fs_jpim1 ! vector opt.
	300	zztmp = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
	301	zztmpx = 0.5 * ( un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) &
	302	& + un(ji ,jj,jk) * un(ji ,jj,jk) * e2u(ji ,jj) * fse3u(ji ,jj,jk) ) &
	303	& * zztmp
	304	!
	305	zztmpy = 0.5 * ( vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) &
	306	& + vn(ji,jj ,jk) * vn(ji,jj ,jk) * e1v(ji,jj ) * fse3v(ji,jj ,jk) ) &
	307	& * zztmp
	308	!
	309	rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
	310	!
	311	ENDDO
[4840]	312	ENDDO
	313	ENDDO
[4990]	314	CALL lbc_lnk( rke, 'T', 1. )
	315	CALL iom_put( "eken", rke )
	316	ENDIF
[6348]	317	!
	318	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
	319	!
[4990]	320	IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[1756]	321	z3d(:,:,jpk) = 0.e0
	322	DO jk = 1, jpkm1
[4761]	323	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
[1756]	324	END DO
	325	CALL iom_put( "u_masstr", z3d ) ! mass transport in i-direction
[4990]	326	ENDIF
	327
	328	IF( iom_use("u_heattr") ) THEN
	329	z2d(:,:) = 0.e0
	330	DO jk = 1, jpkm1
	331	DO jj = 2, jpjm1
	332	DO ji = fs_2, fs_jpim1 ! vector opt.
	333	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
	334	END DO
	335	END DO
	336	END DO
	337	CALL lbc_lnk( z2d, 'U', -1. )
	338	CALL iom_put( "u_heattr", (0.5 * rcp) * z2d ) ! heat transport in i-direction
	339	ENDIF
[4761]	340
[4990]	341	IF( iom_use("u_salttr") ) THEN
[1756]	342	z2d(:,:) = 0.e0
	343	DO jk = 1, jpkm1
	344	DO jj = 2, jpjm1
	345	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	346	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
[1756]	347	END DO
	348	END DO
	349	END DO
	350	CALL lbc_lnk( z2d, 'U', -1. )
[4990]	351	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
	352	ENDIF
[4761]	353
[4990]	354
	355	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[4761]	356	z3d(:,:,jpk) = 0.e0
[1756]	357	DO jk = 1, jpkm1
[4761]	358	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
[1756]	359	END DO
	360	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
[4990]	361	ENDIF
	362
	363	IF( iom_use("v_heattr") ) THEN
	364	z2d(:,:) = 0.e0
	365	DO jk = 1, jpkm1
	366	DO jj = 2, jpjm1
	367	DO ji = fs_2, fs_jpim1 ! vector opt.
	368	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
	369	END DO
	370	END DO
	371	END DO
	372	CALL lbc_lnk( z2d, 'V', -1. )
	373	CALL iom_put( "v_heattr", (0.5 * rcp) * z2d ) ! heat transport in j-direction
	374	ENDIF
[4761]	375
[4990]	376	IF( iom_use("v_salttr") ) THEN
[1756]	377	z2d(:,:) = 0.e0
	378	DO jk = 1, jpkm1
	379	DO jj = 2, jpjm1
	380	DO ji = fs_2, fs_jpim1 ! vector opt.
[4990]	381	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
[1756]	382	END DO
	383	END DO
	384	END DO
	385	CALL lbc_lnk( z2d, 'V', -1. )
[4990]	386	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
[1756]	387	ENDIF
[2528]	388	!
[4990]	389	CALL wrk_dealloc( jpi , jpj , z2d )
[3294]	390	CALL wrk_dealloc( jpi , jpj, jpk , z3d )
[2715]	391	!
[3294]	392	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
	393	!
[1482]	394	END SUBROUTINE dia_wri
	395
	396	#else
[2528]	397	!!----------------------------------------------------------------------
	398	!! Default option use IOIPSL library
	399	!!----------------------------------------------------------------------
	400
[1561]	401	SUBROUTINE dia_wri( kt )
[3]	402	!!---------------------------------------------------------------------
	403	!! * ROUTINE dia_wri *
	404	!!
	405	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	406	!! NETCDF format is used by default
	407	!!
	408	!! ** Method : At the beginning of the first time step (nit000),
	409	!! define all the NETCDF files and fields
	410	!! At each time step call histdef to compute the mean if ncessary
	411	!! Each nwrite time step, output the instantaneous or mean fields
	412	!!----------------------------------------------------------------------
[2715]	413	!!
[3]	414	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[2528]	415	!!
	416	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	417	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	418	INTEGER :: inum = 11 ! temporary logical unit
[3294]	419	INTEGER :: ji, jj, jk ! dummy loop indices
	420	INTEGER :: ierr ! error code return from allocation
[2528]	421	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[3609]	422	INTEGER :: jn, ierror ! local integers
[2528]	423	REAL(wp) :: zsto, zout, zmax, zjulian, zdt ! local scalars
[3294]	424	!!
	425	REAL(wp), POINTER, DIMENSION(:,:) :: zw2d ! 2D workspace
	426	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d ! 3D workspace
[3]	427	!!----------------------------------------------------------------------
[3294]	428	!
	429	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
[1482]	430	!
[3294]	431	CALL wrk_alloc( jpi , jpj , zw2d )
[4292]	432	IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_alloc( jpi , jpj , jpk , zw3d )
[2715]	433	!
[1482]	434	! Output the initial state and forcings
	435	IF( ninist == 1 ) THEN
	436	CALL dia_wri_state( 'output.init', kt )
	437	ninist = 0
	438	ENDIF
	439	!
[3]	440	! 0. Initialisation
	441	! -----------------
[632]	442
[3]	443	! local variable for debugging
	444	ll_print = .FALSE.
	445	ll_print = ll_print .AND. lwp
	446
	447	! Define frequency of output and means
	448	zdt = rdt
[5565]	449	clop = "x" ! no use of the mask value (require less cpu time, and otherwise the model crashes)
[3]	450	#if defined key_diainstant
	451	zsto = nwrite * zdt
[1312]	452	clop = "inst("//TRIM(clop)//")"
[3]	453	#else
	454	zsto=zdt
[1312]	455	clop = "ave("//TRIM(clop)//")"
[3]	456	#endif
	457	zout = nwrite * zdt
	458	zmax = ( nitend - nit000 + 1 ) * zdt
	459
	460	! Define indices of the horizontal output zoom and vertical limit storage
	461	iimi = 1 ; iima = jpi
	462	ijmi = 1 ; ijma = jpj
	463	ipk = jpk
	464
	465	! define time axis
[1334]	466	it = kt
	467	itmod = kt - nit000 + 1
[3]	468
	469
	470	! 1. Define NETCDF files and fields at beginning of first time step
	471	! -----------------------------------------------------------------
	472
	473	IF( kt == nit000 ) THEN
	474
	475	! Define the NETCDF files (one per grid)
[632]	476
[3]	477	! Compute julian date from starting date of the run
[1309]	478	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	479	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	480	IF(lwp)WRITE(numout,*)
	481	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	482	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	483	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	484	' limit storage in depth = ', ipk
	485
	486	! WRITE root name in date.file for use by postpro
[1581]	487	IF(lwp) THEN
[895]	488	CALL dia_nam( clhstnam, nwrite,' ' )
[1581]	489	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	490	WRITE(inum,*) clhstnam
	491	CLOSE(inum)
	492	ENDIF
[632]	493
[3]	494	! Define the T grid FILE ( nid_T )
[632]	495
[3]	496	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
	497	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	498	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	499	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	500	& nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	501	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	502	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	503	! ! Index of ocean points
	504	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	505	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	506	!
	507	IF( ln_icebergs ) THEN
	508	!
	509	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	510	!! that routine is called from nemogcm, so do it here immediately before its needed
	511	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
	512	IF( lk_mpp ) CALL mpp_sum( ierror )
	513	IF( ierror /= 0 ) THEN
	514	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	515	RETURN
	516	ENDIF
	517	!
	518	!! iceberg vertical coordinate is class number
	519	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	520	& "number", nclasses, class_num, nb_T )
	521	!
	522	!! each class just needs the surface index pattern
	523	ndim_bT = 3
	524	DO jn = 1,nclasses
	525	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	526	ENDDO
	527	!
	528	ENDIF
[3]	529
	530	! Define the U grid FILE ( nid_U )
	531
	532	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
	533	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	534	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	535	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	536	& nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	537	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	538	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	539	! ! Index of ocean points
	540	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	541	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	542
	543	! Define the V grid FILE ( nid_V )
	544
	545	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
	546	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	547	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	548	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	549	& nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	550	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	551	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	552	! ! Index of ocean points
	553	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	554	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	555
	556	! Define the W grid FILE ( nid_W )
	557
	558	CALL dia_nam( clhstnam, nwrite, 'grid_W' ) ! filename
	559	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	560	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	561	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	562	& nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	563	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	564	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	565
[632]	566
[3]	567	! Declare all the output fields as NETCDF variables
	568
	569	! !!! nid_T : 3D
	570	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	571	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	572	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	573	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[4292]	574	IF( lk_vvl ) THEN
	575	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	576	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	577	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	578	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	579	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	580	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	581	ENDIF
[3]	582	! !!! nid_T : 2D
	583	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	584	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	585	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	586	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	587	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	588	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	589	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	590	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	591	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	592	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	593	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	594	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4292]	595	IF( .NOT. lk_vvl ) THEN
	596	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	597	& , "KgC/m2/s", & ! sosst_cd
[4292]	598	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	599	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	600	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	601	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	602	ENDIF
[888]	603	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	604	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	605	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	606	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	607	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	608	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	609	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	610	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	611	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	612	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	613	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	614	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	615	!
	616	IF( ln_icebergs ) THEN
	617	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	618	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	619	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	620	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	621	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	622	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	623	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	624	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	625	IF( ln_bergdia ) THEN
	626	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	627	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	628	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	629	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	630	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	631	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	632	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	633	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	634	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	635	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	636	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	637	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	638	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	639	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	640	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	641	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	642	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	643	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	644	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	645	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	646	ENDIF
	647	ENDIF
	648
[5407]	649	IF( .NOT. ln_cpl ) THEN
[4990]	650	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	651	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	652	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	653	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	654	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	655	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	656	ENDIF
[3]	657
[5407]	658	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	659	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	660	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	661	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	662	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	663	CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping" , "Kg/m2/s", & ! erp * sn
	664	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	665	ENDIF
	666
[3]	667	clmx ="l_max(only(x))" ! max index on a period
	668	CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	669	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
	670	#if defined key_diahth
	671	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	672	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	673	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	674	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	675	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	676	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	677	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "W" , & ! htc3
	678	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	679	#endif
	680
[5407]	681	IF( ln_cpl .AND. nn_ice == 2 ) THEN
[4990]	682	CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature" , "K" , & ! tn_ice
	683	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	684	CALL histdef( nid_T,"soicealb" , "Ice Albedo" , "[0,1]" , & ! alb_ice
	685	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	686	ENDIF
[3]	687
[2528]	688	CALL histend( nid_T, snc4chunks=snc4set )
[3]	689
	690	! !!! nid_U : 3D
	691	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	692	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[3294]	693	IF( ln_traldf_gdia ) THEN
	694	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
	695	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	696	ELSE
[3]	697	#if defined key_diaeiv
[3294]	698	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
[3]	699	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	700	#endif
[3294]	701	END IF
[3]	702	! !!! nid_U : 2D
[888]	703	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	704	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	705
[2528]	706	CALL histend( nid_U, snc4chunks=snc4set )
[3]	707
	708	! !!! nid_V : 3D
	709	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	710	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[3294]	711	IF( ln_traldf_gdia ) THEN
	712	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
	713	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	714	ELSE
[3]	715	#if defined key_diaeiv
[3294]	716	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
[3]	717	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	718	#endif
[3294]	719	END IF
[3]	720	! !!! nid_V : 2D
[888]	721	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	722	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	723
[2528]	724	CALL histend( nid_V, snc4chunks=snc4set )
[3]	725
	726	! !!! nid_W : 3D
	727	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	728	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[3294]	729	IF( ln_traldf_gdia ) THEN
	730	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
	731	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	732	ELSE
[3]	733	#if defined key_diaeiv
[3294]	734	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
	735	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[3]	736	#endif
[3294]	737	END IF
[3]	738	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	739	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[255]	740	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avmu
	741	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	742
[3]	743	IF( lk_zdfddm ) THEN
	744	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	745	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	746	ENDIF
	747	! !!! nid_W : 2D
	748	#if defined key_traldf_c2d
	749	CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity" , "m2/s" , & ! ahtw
	750	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
[23]	751	# if defined key_traldf_eiv
[3]	752	CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s", & ! aeiw
	753	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
[23]	754	# endif
[3]	755	#endif
	756
[2528]	757	CALL histend( nid_W, snc4chunks=snc4set )
[3]	758
	759	IF(lwp) WRITE(numout,*)
	760	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	761	IF(ll_print) CALL FLUSH(numout )
	762
	763	ENDIF
	764
	765	! 2. Start writing data
	766	! ---------------------
	767
[4292]	768	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	769	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	770	! donne le nombre d'elements, et ndex la liste des indices a sortir
	771
[1334]	772	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
[3]	773	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	774	WRITE(numout,*) '~~~~~~ '
	775	ENDIF
	776
[4292]	777	IF( lk_vvl ) THEN
[4492]	778	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	779	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	780	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	781	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	782	ELSE
	783	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	784	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	785	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	786	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	787	ENDIF
	788	IF( lk_vvl ) THEN
	789	zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	790	CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	791	CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
	792	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	793	ENDIF
[3]	794	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	795	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	796	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	797	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	798	! (includes virtual salt flux beneath ice
	799	! in linear free surface case)
[4292]	800	IF( .NOT. lk_vvl ) THEN
	801	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
	802	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
	803	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
	804	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	805	ENDIF
[888]	806	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	807	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	808	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	809	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	810	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	811	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	812	!
	813	IF( ln_icebergs ) THEN
	814	!
	815	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	816	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	817	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	818	!
	819	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	820	!
	821	IF( ln_bergdia ) THEN
	822	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	823	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	824	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	825	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	826	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	827	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	828	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	829	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	830	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	831	!
	832	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	833	ENDIF
	834	ENDIF
	835
[5407]	836	IF( .NOT. ln_cpl ) THEN
[4990]	837	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	838	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[3294]	839	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	840	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	841	ENDIF
[5407]	842	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	843	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	844	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
	845	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
	846	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	847	ENDIF
	848	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	849	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	850
	851	#if defined key_diahth
	852	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	853	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	854	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	855	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	856	#endif
[888]	857
[5407]	858	IF( ln_cpl .AND. nn_ice == 2 ) THEN
[4990]	859	CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT ) ! surf. ice temperature
	860	CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT ) ! ice albedo
	861	ENDIF
	862
[3]	863	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[3294]	864	IF( ln_traldf_gdia ) THEN
	865	IF (.not. ALLOCATED(psix_eiv))THEN
	866	ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
	867	IF( lk_mpp ) CALL mpp_sum ( ierr )
	868	IF( ierr > 0 ) CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
	869	psix_eiv(:,:,:) = 0.0_wp
	870	psiy_eiv(:,:,:) = 0.0_wp
	871	ENDIF
	872	DO jk=1,jpkm1
	873	zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz
	874	END DO
	875	zw3d(:,:,jpk) = 0._wp
	876	CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U ) ! i-eiv current
	877	ELSE
[3]	878	#if defined key_diaeiv
[3294]	879	CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U ) ! i-eiv current
[3]	880	#endif
[3294]	881	ENDIF
[888]	882	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	883
	884	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[3294]	885	IF( ln_traldf_gdia ) THEN
	886	DO jk=1,jpk-1
	887	zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz
	888	END DO
	889	zw3d(:,:,jpk) = 0._wp
	890	CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V ) ! j-eiv current
	891	ELSE
[3]	892	#if defined key_diaeiv
[3294]	893	CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V ) ! j-eiv current
[3]	894	#endif
[3294]	895	ENDIF
[888]	896	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	897
	898	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
[3294]	899	IF( ln_traldf_gdia ) THEN
	900	DO jk=1,jpk-1
	901	DO jj = 2, jpjm1
	902	DO ji = fs_2, fs_jpim1 ! vector opt.
	903	zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
	904	& (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
	905	END DO
	906	END DO
	907	END DO
	908	zw3d(:,:,jpk) = 0._wp
	909	CALL histwrite( nid_W, "voveeivw", it, zw3d , ndim_T, ndex_T ) ! vert. eiv current
	910	ELSE
[3]	911	# if defined key_diaeiv
[3294]	912	CALL histwrite( nid_W, "voveeivw", it, w_eiv , ndim_T, ndex_T ) ! vert. eiv current
[3]	913	# endif
[3294]	914	ENDIF
[3]	915	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[255]	916	CALL histwrite( nid_W, "votkeavm", it, avmu , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
[3]	917	IF( lk_zdfddm ) THEN
	918	CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T ) ! S vert. eddy diff. coef.
	919	ENDIF
	920	#if defined key_traldf_c2d
	921	CALL histwrite( nid_W, "soleahtw", it, ahtw , ndim_hT, ndex_hT ) ! lateral eddy diff. coef.
[23]	922	# if defined key_traldf_eiv
[3]	923	CALL histwrite( nid_W, "soleaeiw", it, aeiw , ndim_hT, ndex_hT ) ! EIV coefficient at w-point
[23]	924	# endif
[3]	925	#endif
	926
[1318]	927	! 3. Close all files
	928	! ---------------------------------------
[1561]	929	IF( kt == nitend ) THEN
[3]	930	CALL histclo( nid_T )
	931	CALL histclo( nid_U )
	932	CALL histclo( nid_V )
	933	CALL histclo( nid_W )
	934	ENDIF
[2528]	935	!
[3294]	936	CALL wrk_dealloc( jpi , jpj , zw2d )
[4292]	937	IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_dealloc( jpi , jpj , jpk , zw3d )
[2715]	938	!
[3294]	939	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
	940	!
[3]	941	END SUBROUTINE dia_wri
[1482]	942	# endif
[3]	943
[1567]	944	#endif
	945
[1334]	946	SUBROUTINE dia_wri_state( cdfile_name, kt )
[3]	947	!!---------------------------------------------------------------------
	948	!! * ROUTINE dia_wri_state *
	949	!!
	950	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	951	!! the instantaneous ocean state and forcing fields.
	952	!! Used to find errors in the initial state or save the last
	953	!! ocean state in case of abnormal end of a simulation
	954	!!
	955	!! ** Method : NetCDF files using ioipsl
	956	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	957	!! File 'output.abort.nc' is created in case of abnormal job end
	958	!!----------------------------------------------------------------------
[1334]	959	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
	960	INTEGER , INTENT( in ) :: kt ! ocean time-step index
[2528]	961	!!
[648]	962	CHARACTER (len=32) :: clname
[3]	963	CHARACTER (len=40) :: clop
[2528]	964	INTEGER :: id_i , nz_i, nh_i
	965	INTEGER, DIMENSION(1) :: idex ! local workspace
	966	REAL(wp) :: zsto, zout, zmax, zjulian, zdt
[3]	967	!!----------------------------------------------------------------------
[3294]	968	!
[3625]	969	! IF( nn_timing == 1 ) CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
[3]	970
	971	! 0. Initialisation
	972	! -----------------
[632]	973
[648]	974	! Define name, frequency of output and means
	975	clname = cdfile_name
[1792]	976	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
[3]	977	zdt = rdt
	978	zsto = rdt
	979	clop = "inst(x)" ! no use of the mask value (require less cpu time)
	980	zout = rdt
	981	zmax = ( nitend - nit000 + 1 ) * zdt
	982
[648]	983	IF(lwp) WRITE(numout,*)
	984	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	985	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
	986	IF(lwp) WRITE(numout,*) ' and named :', clname, '.nc'
	987
	988
[3]	989	! 1. Define NETCDF files and fields at beginning of first time step
	990	! -----------------------------------------------------------------
	991
	992	! Compute julian date from starting date of the run
[1310]	993	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis
	994	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[648]	995	CALL histbeg( clname, jpi, glamt, jpj, gphit, &
[2528]	996	1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
[3]	997	CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept
[4292]	998	"m", jpk, gdept_1d, nz_i, "down")
[3]	999
	1000	! Declare all the output fields as NetCDF variables
	1001
	1002	CALL histdef( id_i, "vosaline", "Salinity" , "PSU" , & ! salinity
	1003	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1004	CALL histdef( id_i, "votemper", "Temperature" , "C" , & ! temperature
	1005	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[359]	1006	CALL histdef( id_i, "sossheig", "Sea Surface Height" , "m" , & ! ssh
	1007	& jpi, jpj, nh_i, 1 , 1, 1 , nz_i, 32, clop, zsto, zout )
[3]	1008	CALL histdef( id_i, "vozocrtx", "Zonal Current" , "m/s" , & ! zonal current
	1009	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1010	CALL histdef( id_i, "vomecrty", "Meridional Current" , "m/s" , & ! meridonal current
	1011	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1012	CALL histdef( id_i, "vovecrtz", "Vertical Velocity" , "m/s" , & ! vertical current
	1013	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1014	CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S", & ! net freshwater
	1015	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1016	CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2" , & ! net heat flux
	1017	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1018	CALL histdef( id_i, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! solar flux
	1019	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	1020	CALL histdef( id_i, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	1021	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1022	CALL histdef( id_i, "sozotaux", "Zonal Wind Stress" , "N/m2" , & ! i-wind stress
	1023	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1024	CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress
	1025	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4292]	1026	IF( lk_vvl ) THEN
	1027	CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth
	1028	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[6204]	1029	CALL histdef( id_i, "vovvle3t", "T point thickness" , "m" , & ! t-point depth
	1030	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[4292]	1031	END IF
[3]	1032
[1482]	1033	#if defined key_lim2
	1034	CALL lim_wri_state_2( kt, id_i, nh_i )
[4161]	1035	#elif defined key_lim3
	1036	CALL lim_wri_state( kt, id_i, nh_i )
[1482]	1037	#else
[2528]	1038	CALL histend( id_i, snc4chunks=snc4set )
[1482]	1039	#endif
[3]	1040
	1041	! 2. Start writing data
	1042	! ---------------------
	1043	! idex(1) est utilise ssi l'avant dernier argument est diffferent de
	1044	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	1045	! donne le nombre d'elements, et idex la liste des indices a sortir
	1046	idex(1) = 1 ! init to avoid compil warning
[632]	1047
[3]	1048	! Write all fields on T grid
[3294]	1049	CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpijpjjpk, idex ) ! now temperature
	1050	CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpijpjjpk, idex ) ! now salinity
	1051	CALL histwrite( id_i, "sossheig", kt, sshn , jpi*jpj , idex ) ! sea surface height
	1052	CALL histwrite( id_i, "vozocrtx", kt, un , jpijpjjpk, idex ) ! now i-velocity
	1053	CALL histwrite( id_i, "vomecrty", kt, vn , jpijpjjpk, idex ) ! now j-velocity
	1054	CALL histwrite( id_i, "vovecrtz", kt, wn , jpijpjjpk, idex ) ! now k-velocity
	1055	CALL histwrite( id_i, "sowaflup", kt, (emp-rnf ) , jpi*jpj , idex ) ! freshwater budget
	1056	CALL histwrite( id_i, "sohefldo", kt, qsr + qns , jpi*jpj , idex ) ! total heat flux
	1057	CALL histwrite( id_i, "soshfldo", kt, qsr , jpi*jpj , idex ) ! solar heat flux
	1058	CALL histwrite( id_i, "soicecov", kt, fr_i , jpi*jpj , idex ) ! ice fraction
	1059	CALL histwrite( id_i, "sozotaux", kt, utau , jpi*jpj , idex ) ! i-wind stress
	1060	CALL histwrite( id_i, "sometauy", kt, vtau , jpi*jpj , idex ) ! j-wind stress
[6204]	1061	IF( lk_vvl ) THEN
	1062	CALL histwrite( id_i, "vovvldep", kt, fsdept_n(:,:,:), jpijpjjpk, idex )! T-cell depth
	1063	CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpijpjjpk, idex )! T-cell thickness
	1064	END IF
[3]	1065
	1066	! 3. Close the file
	1067	! -----------------
	1068	CALL histclo( id_i )
[1567]	1069	#if ! defined key_iomput && ! defined key_dimgout
[1561]	1070	IF( ninist /= 1 ) THEN
	1071	CALL histclo( nid_T )
	1072	CALL histclo( nid_U )
	1073	CALL histclo( nid_V )
	1074	CALL histclo( nid_W )
	1075	ENDIF
	1076	#endif
[3294]	1077
[3625]	1078	! IF( nn_timing == 1 ) CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
[3294]	1079	!
[3]	1080
	1081	END SUBROUTINE dia_wri_state
	1082	!!======================================================================
	1083	END MODULE diawri

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

Download in other formats: