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

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

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

Context Navigation

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

Last change on this file since 14210 was 14179, checked in by techene, 4 years ago
#2385 restart with qco bug fix for euler 1st and RK3
Property svn:keywords set to `Id`
File size: 64.7 KB

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

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

Download in other formats: