Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Normal
Revision Log

diawri.F90 @ 11639

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

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

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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package_reanalysis4/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 11639

Download in other formats: