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/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/2015/nemo_v3_6_STABLE/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6348

Last change on this file since 6348 was 6348, checked in by cetlod, 8 years ago
bugfix: move the output of scale factor before time swapping and output some variables needs for offline, see ticket #1682
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
	449	IF( nacc == 1 ) zdt = rdtmin
[5565]	450	clop = "x" ! no use of the mask value (require less cpu time, and otherwise the model crashes)
[3]	451	#if defined key_diainstant
	452	zsto = nwrite * zdt
[1312]	453	clop = "inst("//TRIM(clop)//")"
[3]	454	#else
	455	zsto=zdt
[1312]	456	clop = "ave("//TRIM(clop)//")"
[3]	457	#endif
	458	zout = nwrite * zdt
	459	zmax = ( nitend - nit000 + 1 ) * zdt
	460
	461	! Define indices of the horizontal output zoom and vertical limit storage
	462	iimi = 1 ; iima = jpi
	463	ijmi = 1 ; ijma = jpj
	464	ipk = jpk
	465
	466	! define time axis
[1334]	467	it = kt
	468	itmod = kt - nit000 + 1
[3]	469
	470
	471	! 1. Define NETCDF files and fields at beginning of first time step
	472	! -----------------------------------------------------------------
	473
	474	IF( kt == nit000 ) THEN
	475
	476	! Define the NETCDF files (one per grid)
[632]	477
[3]	478	! Compute julian date from starting date of the run
[1309]	479	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
	480	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	481	IF(lwp)WRITE(numout,*)
	482	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	483	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	484	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	485	' limit storage in depth = ', ipk
	486
	487	! WRITE root name in date.file for use by postpro
[1581]	488	IF(lwp) THEN
[895]	489	CALL dia_nam( clhstnam, nwrite,' ' )
[1581]	490	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	491	WRITE(inum,*) clhstnam
	492	CLOSE(inum)
	493	ENDIF
[632]	494
[3]	495	! Define the T grid FILE ( nid_T )
[632]	496
[3]	497	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
	498	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	499	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	500	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	501	& nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	502	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	503	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	504	! ! Index of ocean points
	505	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	506	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	507	!
	508	IF( ln_icebergs ) THEN
	509	!
	510	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	511	!! that routine is called from nemogcm, so do it here immediately before its needed
	512	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
	513	IF( lk_mpp ) CALL mpp_sum( ierror )
	514	IF( ierror /= 0 ) THEN
	515	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	516	RETURN
	517	ENDIF
	518	!
	519	!! iceberg vertical coordinate is class number
	520	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	521	& "number", nclasses, class_num, nb_T )
	522	!
	523	!! each class just needs the surface index pattern
	524	ndim_bT = 3
	525	DO jn = 1,nclasses
	526	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	527	ENDDO
	528	!
	529	ENDIF
[3]	530
	531	! Define the U grid FILE ( nid_U )
	532
	533	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
	534	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	535	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	536	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	537	& nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	538	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	539	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	540	! ! Index of ocean points
	541	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	542	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	543
	544	! Define the V grid FILE ( nid_V )
	545
	546	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
	547	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	548	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	549	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	550	& nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	551	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	552	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	553	! ! Index of ocean points
	554	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	555	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	556
	557	! Define the W grid FILE ( nid_W )
	558
	559	CALL dia_nam( clhstnam, nwrite, 'grid_W' ) ! filename
	560	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	561	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	562	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[2528]	563	& nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	564	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	565	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	566
[632]	567
[3]	568	! Declare all the output fields as NETCDF variables
	569
	570	! !!! nid_T : 3D
	571	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	572	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	573	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	574	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[4292]	575	IF( lk_vvl ) THEN
	576	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
	577	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	578	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
	579	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	580	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
	581	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	582	ENDIF
[3]	583	! !!! nid_T : 2D
	584	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	585	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	586	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	587	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	588	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	589	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	590	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	591	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	592	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	593	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	594	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	595	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4292]	596	IF( .NOT. lk_vvl ) THEN
	597	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
[3625]	598	& , "KgC/m2/s", & ! sosst_cd
[4292]	599	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	600	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
[3625]	601	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	602	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	603	ENDIF
[888]	604	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	605	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	606	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	607	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	608	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	609	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	610	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	611	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	612	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	613	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	614	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	615	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3609]	616	!
	617	IF( ln_icebergs ) THEN
	618	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	619	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	620	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	621	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	622	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	623	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	624	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	625	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	626	IF( ln_bergdia ) THEN
	627	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	628	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	629	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	630	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	631	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	632	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	633	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	634	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	635	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	636	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	637	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	638	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	639	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	640	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	641	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	642	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	643	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	644	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	645	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	646	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	647	ENDIF
	648	ENDIF
	649
[5407]	650	IF( .NOT. ln_cpl ) THEN
[4990]	651	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	652	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	653	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	654	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	655	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	656	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	657	ENDIF
[3]	658
[5407]	659	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	660	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	661	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	662	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	663	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	664	CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping" , "Kg/m2/s", & ! erp * sn
	665	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	666	ENDIF
	667
[3]	668	clmx ="l_max(only(x))" ! max index on a period
	669	CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	670	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
	671	#if defined key_diahth
	672	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	673	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	674	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	675	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	676	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	677	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	678	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "W" , & ! htc3
	679	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	680	#endif
	681
[5407]	682	IF( ln_cpl .AND. nn_ice == 2 ) THEN
[4990]	683	CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature" , "K" , & ! tn_ice
	684	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	685	CALL histdef( nid_T,"soicealb" , "Ice Albedo" , "[0,1]" , & ! alb_ice
	686	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	687	ENDIF
[3]	688
[2528]	689	CALL histend( nid_T, snc4chunks=snc4set )
[3]	690
	691	! !!! nid_U : 3D
	692	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
	693	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[3294]	694	IF( ln_traldf_gdia ) THEN
	695	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
	696	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	697	ELSE
[3]	698	#if defined key_diaeiv
[3294]	699	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
[3]	700	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	701	#endif
[3294]	702	END IF
[3]	703	! !!! nid_U : 2D
[888]	704	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	705	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	706
[2528]	707	CALL histend( nid_U, snc4chunks=snc4set )
[3]	708
	709	! !!! nid_V : 3D
	710	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
	711	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[3294]	712	IF( ln_traldf_gdia ) THEN
	713	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
	714	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	715	ELSE
[3]	716	#if defined key_diaeiv
[3294]	717	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
[3]	718	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	719	#endif
[3294]	720	END IF
[3]	721	! !!! nid_V : 2D
[888]	722	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	723	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	724
[2528]	725	CALL histend( nid_V, snc4chunks=snc4set )
[3]	726
	727	! !!! nid_W : 3D
	728	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
	729	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[3294]	730	IF( ln_traldf_gdia ) THEN
	731	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
	732	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	733	ELSE
[3]	734	#if defined key_diaeiv
[3294]	735	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
	736	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[3]	737	#endif
[3294]	738	END IF
[3]	739	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	740	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[255]	741	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avmu
	742	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	743
[3]	744	IF( lk_zdfddm ) THEN
	745	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	746	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	747	ENDIF
	748	! !!! nid_W : 2D
	749	#if defined key_traldf_c2d
	750	CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity" , "m2/s" , & ! ahtw
	751	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
[23]	752	# if defined key_traldf_eiv
[3]	753	CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s", & ! aeiw
	754	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
[23]	755	# endif
[3]	756	#endif
	757
[2528]	758	CALL histend( nid_W, snc4chunks=snc4set )
[3]	759
	760	IF(lwp) WRITE(numout,*)
	761	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	762	IF(ll_print) CALL FLUSH(numout )
	763
	764	ENDIF
	765
	766	! 2. Start writing data
	767	! ---------------------
	768
[4292]	769	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	770	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	771	! donne le nombre d'elements, et ndex la liste des indices a sortir
	772
[1334]	773	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
[3]	774	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	775	WRITE(numout,*) '~~~~~~ '
	776	ENDIF
	777
[4292]	778	IF( lk_vvl ) THEN
[4492]	779	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
	780	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * fse3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
	781	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	782	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * fse3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	783	ELSE
	784	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
	785	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
	786	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
	787	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
	788	ENDIF
	789	IF( lk_vvl ) THEN
	790	zw3d(:,:,:) = ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	791	CALL histwrite( nid_T, "vovvle3t", it, fse3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
	792	CALL histwrite( nid_T, "vovvldep", it, fsdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
	793	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	794	ENDIF
[3]	795	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
[2528]	796	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	797	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	798	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	799	! (includes virtual salt flux beneath ice
	800	! in linear free surface case)
[4292]	801	IF( .NOT. lk_vvl ) THEN
	802	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
	803	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
	804	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
	805	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	806	ENDIF
[888]	807	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	808	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	809	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	810	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	811	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	812	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[3609]	813	!
	814	IF( ln_icebergs ) THEN
	815	!
	816	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	817	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	818	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	819	!
	820	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	821	!
	822	IF( ln_bergdia ) THEN
	823	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	824	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	825	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	826	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	827	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	828	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	829	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	830	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	831	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	832	!
	833	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	834	ENDIF
	835	ENDIF
	836
[5407]	837	IF( .NOT. ln_cpl ) THEN
[4990]	838	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	839	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[3294]	840	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
[4990]	841	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	842	ENDIF
[5407]	843	IF( ln_cpl .AND. nn_ice <= 1 ) THEN
[4990]	844	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	845	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
	846	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
	847	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	848	ENDIF
	849	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	850	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	851
	852	#if defined key_diahth
	853	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	854	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	855	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	856	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	857	#endif
[888]	858
[5407]	859	IF( ln_cpl .AND. nn_ice == 2 ) THEN
[4990]	860	CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT ) ! surf. ice temperature
	861	CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT ) ! ice albedo
	862	ENDIF
	863
[3]	864	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
[3294]	865	IF( ln_traldf_gdia ) THEN
	866	IF (.not. ALLOCATED(psix_eiv))THEN
	867	ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
	868	IF( lk_mpp ) CALL mpp_sum ( ierr )
	869	IF( ierr > 0 ) CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
	870	psix_eiv(:,:,:) = 0.0_wp
	871	psiy_eiv(:,:,:) = 0.0_wp
	872	ENDIF
	873	DO jk=1,jpkm1
	874	zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz
	875	END DO
	876	zw3d(:,:,jpk) = 0._wp
	877	CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U ) ! i-eiv current
	878	ELSE
[3]	879	#if defined key_diaeiv
[3294]	880	CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U ) ! i-eiv current
[3]	881	#endif
[3294]	882	ENDIF
[888]	883	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	884
	885	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
[3294]	886	IF( ln_traldf_gdia ) THEN
	887	DO jk=1,jpk-1
	888	zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz
	889	END DO
	890	zw3d(:,:,jpk) = 0._wp
	891	CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V ) ! j-eiv current
	892	ELSE
[3]	893	#if defined key_diaeiv
[3294]	894	CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V ) ! j-eiv current
[3]	895	#endif
[3294]	896	ENDIF
[888]	897	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	898
	899	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
[3294]	900	IF( ln_traldf_gdia ) THEN
	901	DO jk=1,jpk-1
	902	DO jj = 2, jpjm1
	903	DO ji = fs_2, fs_jpim1 ! vector opt.
	904	zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
	905	& (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
	906	END DO
	907	END DO
	908	END DO
	909	zw3d(:,:,jpk) = 0._wp
	910	CALL histwrite( nid_W, "voveeivw", it, zw3d , ndim_T, ndex_T ) ! vert. eiv current
	911	ELSE
[3]	912	# if defined key_diaeiv
[3294]	913	CALL histwrite( nid_W, "voveeivw", it, w_eiv , ndim_T, ndex_T ) ! vert. eiv current
[3]	914	# endif
[3294]	915	ENDIF
[3]	916	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[255]	917	CALL histwrite( nid_W, "votkeavm", it, avmu , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
[3]	918	IF( lk_zdfddm ) THEN
	919	CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T ) ! S vert. eddy diff. coef.
	920	ENDIF
	921	#if defined key_traldf_c2d
	922	CALL histwrite( nid_W, "soleahtw", it, ahtw , ndim_hT, ndex_hT ) ! lateral eddy diff. coef.
[23]	923	# if defined key_traldf_eiv
[3]	924	CALL histwrite( nid_W, "soleaeiw", it, aeiw , ndim_hT, ndex_hT ) ! EIV coefficient at w-point
[23]	925	# endif
[3]	926	#endif
	927
[1318]	928	! 3. Close all files
	929	! ---------------------------------------
[1561]	930	IF( kt == nitend ) THEN
[3]	931	CALL histclo( nid_T )
	932	CALL histclo( nid_U )
	933	CALL histclo( nid_V )
	934	CALL histclo( nid_W )
	935	ENDIF
[2528]	936	!
[3294]	937	CALL wrk_dealloc( jpi , jpj , zw2d )
[4292]	938	IF ( ln_traldf_gdia .OR. lk_vvl ) call wrk_dealloc( jpi , jpj , jpk , zw3d )
[2715]	939	!
[3294]	940	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
	941	!
[3]	942	END SUBROUTINE dia_wri
[1482]	943	# endif
[3]	944
[1567]	945	#endif
	946
[1334]	947	SUBROUTINE dia_wri_state( cdfile_name, kt )
[3]	948	!!---------------------------------------------------------------------
	949	!! * ROUTINE dia_wri_state *
	950	!!
	951	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	952	!! the instantaneous ocean state and forcing fields.
	953	!! Used to find errors in the initial state or save the last
	954	!! ocean state in case of abnormal end of a simulation
	955	!!
	956	!! ** Method : NetCDF files using ioipsl
	957	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	958	!! File 'output.abort.nc' is created in case of abnormal job end
	959	!!----------------------------------------------------------------------
[1334]	960	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
	961	INTEGER , INTENT( in ) :: kt ! ocean time-step index
[2528]	962	!!
[648]	963	CHARACTER (len=32) :: clname
[3]	964	CHARACTER (len=40) :: clop
[2528]	965	INTEGER :: id_i , nz_i, nh_i
	966	INTEGER, DIMENSION(1) :: idex ! local workspace
	967	REAL(wp) :: zsto, zout, zmax, zjulian, zdt
[3]	968	!!----------------------------------------------------------------------
[3294]	969	!
[3625]	970	! IF( nn_timing == 1 ) CALL timing_start('dia_wri_state') ! not sure this works for routines not called in first timestep
[3]	971
	972	! 0. Initialisation
	973	! -----------------
[632]	974
[648]	975	! Define name, frequency of output and means
	976	clname = cdfile_name
[1792]	977	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
[3]	978	zdt = rdt
	979	zsto = rdt
	980	clop = "inst(x)" ! no use of the mask value (require less cpu time)
	981	zout = rdt
	982	zmax = ( nitend - nit000 + 1 ) * zdt
	983
[648]	984	IF(lwp) WRITE(numout,*)
	985	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	986	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
	987	IF(lwp) WRITE(numout,*) ' and named :', clname, '.nc'
	988
	989
[3]	990	! 1. Define NETCDF files and fields at beginning of first time step
	991	! -----------------------------------------------------------------
	992
	993	! Compute julian date from starting date of the run
[1310]	994	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis
	995	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[648]	996	CALL histbeg( clname, jpi, glamt, jpj, gphit, &
[2528]	997	1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
[3]	998	CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept
[4292]	999	"m", jpk, gdept_1d, nz_i, "down")
[3]	1000
	1001	! Declare all the output fields as NetCDF variables
	1002
	1003	CALL histdef( id_i, "vosaline", "Salinity" , "PSU" , & ! salinity
	1004	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1005	CALL histdef( id_i, "votemper", "Temperature" , "C" , & ! temperature
	1006	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[359]	1007	CALL histdef( id_i, "sossheig", "Sea Surface Height" , "m" , & ! ssh
	1008	& jpi, jpj, nh_i, 1 , 1, 1 , nz_i, 32, clop, zsto, zout )
[3]	1009	CALL histdef( id_i, "vozocrtx", "Zonal Current" , "m/s" , & ! zonal current
	1010	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1011	CALL histdef( id_i, "vomecrty", "Meridional Current" , "m/s" , & ! meridonal current
	1012	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1013	CALL histdef( id_i, "vovecrtz", "Vertical Velocity" , "m/s" , & ! vertical current
	1014	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
	1015	CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S", & ! net freshwater
	1016	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1017	CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2" , & ! net heat flux
	1018	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1019	CALL histdef( id_i, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! solar flux
	1020	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	1021	CALL histdef( id_i, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	1022	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1023	CALL histdef( id_i, "sozotaux", "Zonal Wind Stress" , "N/m2" , & ! i-wind stress
	1024	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	1025	CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress
	1026	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4292]	1027	IF( lk_vvl ) THEN
	1028	CALL histdef( id_i, "vovvldep", "T point depth" , "m" , & ! t-point depth
	1029	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[6204]	1030	CALL histdef( id_i, "vovvle3t", "T point thickness" , "m" , & ! t-point depth
	1031	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
[4292]	1032	END IF
[3]	1033
[1482]	1034	#if defined key_lim2
	1035	CALL lim_wri_state_2( kt, id_i, nh_i )
[4161]	1036	#elif defined key_lim3
	1037	CALL lim_wri_state( kt, id_i, nh_i )
[1482]	1038	#else
[2528]	1039	CALL histend( id_i, snc4chunks=snc4set )
[1482]	1040	#endif
[3]	1041
	1042	! 2. Start writing data
	1043	! ---------------------
	1044	! idex(1) est utilise ssi l'avant dernier argument est diffferent de
	1045	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	1046	! donne le nombre d'elements, et idex la liste des indices a sortir
	1047	idex(1) = 1 ! init to avoid compil warning
[632]	1048
[3]	1049	! Write all fields on T grid
[3294]	1050	CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpijpjjpk, idex ) ! now temperature
	1051	CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpijpjjpk, idex ) ! now salinity
	1052	CALL histwrite( id_i, "sossheig", kt, sshn , jpi*jpj , idex ) ! sea surface height
	1053	CALL histwrite( id_i, "vozocrtx", kt, un , jpijpjjpk, idex ) ! now i-velocity
	1054	CALL histwrite( id_i, "vomecrty", kt, vn , jpijpjjpk, idex ) ! now j-velocity
	1055	CALL histwrite( id_i, "vovecrtz", kt, wn , jpijpjjpk, idex ) ! now k-velocity
	1056	CALL histwrite( id_i, "sowaflup", kt, (emp-rnf ) , jpi*jpj , idex ) ! freshwater budget
	1057	CALL histwrite( id_i, "sohefldo", kt, qsr + qns , jpi*jpj , idex ) ! total heat flux
	1058	CALL histwrite( id_i, "soshfldo", kt, qsr , jpi*jpj , idex ) ! solar heat flux
	1059	CALL histwrite( id_i, "soicecov", kt, fr_i , jpi*jpj , idex ) ! ice fraction
	1060	CALL histwrite( id_i, "sozotaux", kt, utau , jpi*jpj , idex ) ! i-wind stress
	1061	CALL histwrite( id_i, "sometauy", kt, vtau , jpi*jpj , idex ) ! j-wind stress
[6204]	1062	IF( lk_vvl ) THEN
	1063	CALL histwrite( id_i, "vovvldep", kt, fsdept_n(:,:,:), jpijpjjpk, idex )! T-cell depth
	1064	CALL histwrite( id_i, "vovvle3t", kt, fse3t_n (:,:,:), jpijpjjpk, idex )! T-cell thickness
	1065	END IF
[3]	1066
	1067	! 3. Close the file
	1068	! -----------------
	1069	CALL histclo( id_i )
[1567]	1070	#if ! defined key_iomput && ! defined key_dimgout
[1561]	1071	IF( ninist /= 1 ) THEN
	1072	CALL histclo( nid_T )
	1073	CALL histclo( nid_U )
	1074	CALL histclo( nid_V )
	1075	CALL histclo( nid_W )
	1076	ENDIF
	1077	#endif
[3294]	1078
[3625]	1079	! IF( nn_timing == 1 ) CALL timing_stop('dia_wri_state') ! not sure this works for routines not called in first timestep
[3294]	1080	!
[3]	1081
	1082	END SUBROUTINE dia_wri_state
	1083	!!======================================================================
	1084	END MODULE diawri

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

Download in other formats: