Context Navigation

diawri.F90 @ 15199

Last change on this file since 15199 was 15141, checked in by smasson, 3 years ago
trunk: avoid implicit loops in diawri
Property svn:keywords set to `Id`
File size: 65.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
[14239]	72	#if ! defined key_xios
[12377]	73	PUBLIC dia_wri_alloc_abl ! Called by sbcabl module (if ln_abl = .true.)
	74	#endif
[2528]	75	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
[3609]	76	INTEGER :: nb_T , ndim_bT ! grid_T file
[2528]	77	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
	78	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
	79	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
[12377]	80	INTEGER :: nid_A, nz_A, nh_A, ndim_A, ndim_hA ! grid_ABL file
[2528]	81	INTEGER :: ndex(1) ! ???
[2715]	82	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
[12377]	83	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hA, ndex_A ! ABL
[2715]	84	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
[3609]	85	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
[3]	86
	87	!! * Substitutions
[12377]	88	# include "do_loop_substitute.h90"
[13237]	89	# include "domzgr_substitute.h90"
[3]	90	!!----------------------------------------------------------------------
[9598]	91	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
[5217]	92	!! $Id$
[10068]	93	!! Software governed by the CeCILL license (see ./LICENSE)
[3]	94	!!----------------------------------------------------------------------
	95	CONTAINS
	96
[14239]	97	#if defined key_xios
[3]	98	!!----------------------------------------------------------------------
[14239]	99	!! 'key_xios' use IOM library
[2528]	100	!!----------------------------------------------------------------------
[9652]	101	INTEGER FUNCTION dia_wri_alloc()
	102	!
	103	dia_wri_alloc = 0
	104	!
	105	END FUNCTION dia_wri_alloc
[2715]	106
[9652]	107
[12377]	108	SUBROUTINE dia_wri( kt, Kmm )
[1482]	109	!!---------------------------------------------------------------------
	110	!! * ROUTINE dia_wri *
	111	!!
	112	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	113	!! NETCDF format is used by default
	114	!!
	115	!! ** Method : use iom_put
	116	!!----------------------------------------------------------------------
	117	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	118	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[1756]	119	!!
[9019]	120	INTEGER :: ji, jj, jk ! dummy loop indices
	121	INTEGER :: ikbot ! local integer
	122	REAL(wp):: zztmp , zztmpx ! local scalar
	123	REAL(wp):: zztmp2, zztmpy ! - -
[14053]	124	REAL(wp):: ze3
[15017]	125	REAL(wp), DIMENSION(A2D( 0)) :: z2d ! 2D workspace
	126	REAL(wp), DIMENSION(A2D(nn_hls),jpk) :: z3d ! 3D workspace
[1482]	127	!!----------------------------------------------------------------------
	128	!
[9124]	129	IF( ln_timing ) CALL timing_start('dia_wri')
[3294]	130	!
[1482]	131	! Output the initial state and forcings
	132	IF( ninist == 1 ) THEN
[12377]	133	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	134	ninist = 0
	135	ENDIF
[3]	136
[15087]	137	! initialize arrays
	138	z2d(:,:) = 0._wp
	139	z3d(:,:,:) = 0._wp
	140
[6351]	141	! Output of initial vertical scale factor
	142	CALL iom_put("e3t_0", e3t_0(:,:,:) )
[10114]	143	CALL iom_put("e3u_0", e3u_0(:,:,:) )
	144	CALL iom_put("e3v_0", e3v_0(:,:,:) )
[14053]	145	CALL iom_put("e3f_0", e3f_0(:,:,:) )
[6351]	146	!
[14086]	147	IF ( iom_use("tpt_dep") ) THEN
[15087]	148	DO_3D( 0, 0, 0, 0, 1, jpk )
	149	z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm)
	150	END_3D
[15017]	151	CALL iom_put( "tpt_dep", z3d )
[14086]	152	ENDIF
	153
[15087]	154	! --- vertical scale factors --- !
[13237]	155	IF ( iom_use("e3t") .OR. iom_use("e3tdef") ) THEN ! time-varying e3t
[15087]	156	DO_3D( 0, 0, 0, 0, 1, jpk )
	157	z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm)
	158	END_3D
[15017]	159	CALL iom_put( "e3t", z3d )
	160	IF ( iom_use("e3tdef") ) THEN
[15087]	161	DO_3D( 0, 0, 0, 0, 1, jpk )
	162	z3d(ji,jj,jk) = ( ( z3d(ji,jj,jk) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
	163	END_3D
[15017]	164	CALL iom_put( "e3tdef", z3d )
	165	ENDIF
[13237]	166	ENDIF
	167	IF ( iom_use("e3u") ) THEN ! time-varying e3u
[15087]	168	DO_3D( 0, 0, 0, 0, 1, jpk )
	169	z3d(ji,jj,jk) = e3u(ji,jj,jk,Kmm)
	170	END_3D
[15017]	171	CALL iom_put( "e3u" , z3d )
[13237]	172	ENDIF
	173	IF ( iom_use("e3v") ) THEN ! time-varying e3v
[15087]	174	DO_3D( 0, 0, 0, 0, 1, jpk )
	175	z3d(ji,jj,jk) = e3v(ji,jj,jk,Kmm)
	176	END_3D
[15017]	177	CALL iom_put( "e3v" , z3d )
[13237]	178	ENDIF
	179	IF ( iom_use("e3w") ) THEN ! time-varying e3w
[15087]	180	DO_3D( 0, 0, 0, 0, 1, jpk )
	181	z3d(ji,jj,jk) = e3w(ji,jj,jk,Kmm)
	182	END_3D
[15017]	183	CALL iom_put( "e3w" , z3d )
[13237]	184	ENDIF
[14053]	185	IF ( iom_use("e3f") ) THEN ! time-varying e3f caution here at Kaa
[15087]	186	DO_3D( 0, 0, 0, 0, 1, jpk )
	187	z3d(ji,jj,jk) = e3f(ji,jj,jk)
	188	END_3D
[15017]	189	CALL iom_put( "e3f" , z3d )
[14053]	190	ENDIF
[5461]	191
[15017]	192	IF ( iom_use("ssh") ) THEN
	193	IF( ll_wd ) THEN ! sea surface height (brought back to the reference used for wetting and drying)
	194	CALL iom_put( "ssh" , (ssh(:,:,Kmm)+ssh_ref)*ssmask(:,:) )
	195	ELSE
	196	CALL iom_put( "ssh" , ssh(:,:,Kmm) ) ! sea surface height
	197	ENDIF
[9023]	198	ENDIF
	199
[14053]	200	IF( iom_use("wetdep") ) CALL iom_put( "wetdep" , ht_0(:,:) + ssh(:,:,Kmm) ) ! wet depth
	201
	202	#if defined key_qco
	203	IF( iom_use("ht") ) CALL iom_put( "ht" , ht(:,:) ) ! water column at t-point
	204	IF( iom_use("hu") ) CALL iom_put( "hu" , hu(:,:,Kmm) ) ! water column at u-point
	205	IF( iom_use("hv") ) CALL iom_put( "hv" , hv(:,:,Kmm) ) ! water column at v-point
	206	IF( iom_use("hf") ) CALL iom_put( "hf" , hf_0(:,:)*( 1._wp + r3f(:,:) ) ) ! water column at f-point (caution here at Naa)
	207	#endif
[15087]	208
	209	! --- tracers T&S --- !
[12377]	210	CALL iom_put( "toce", ts(:,:,:,jp_tem,Kmm) ) ! 3D temperature
	211	CALL iom_put( "sst", ts(:,:,1,jp_tem,Kmm) ) ! surface temperature
[15087]	212
[5107]	213	IF ( iom_use("sbt") ) THEN
[13472]	214	DO_2D( 0, 0, 0, 0 )
[12377]	215	ikbot = mbkt(ji,jj)
	216	z2d(ji,jj) = ts(ji,jj,ikbot,jp_tem,Kmm)
	217	END_2D
[5107]	218	CALL iom_put( "sbt", z2d ) ! bottom temperature
	219	ENDIF
	220
[12377]	221	CALL iom_put( "soce", ts(:,:,:,jp_sal,Kmm) ) ! 3D salinity
	222	CALL iom_put( "sss", ts(:,:,1,jp_sal,Kmm) ) ! surface salinity
[5107]	223	IF ( iom_use("sbs") ) THEN
[13472]	224	DO_2D( 0, 0, 0, 0 )
[12377]	225	ikbot = mbkt(ji,jj)
	226	z2d(ji,jj) = ts(ji,jj,ikbot,jp_sal,Kmm)
	227	END_2D
[5107]	228	CALL iom_put( "sbs", z2d ) ! bottom salinity
	229	ENDIF
[5463]	230
[14179]	231	IF( .NOT.lk_SWE ) CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
[13089]	232
[15087]	233	! --- momentum --- !
[5463]	234	IF ( iom_use("taubot") ) THEN ! bottom stress
[15017]	235	zztmp = rho0 * 0.25_wp
[7753]	236	z2d(:,:) = 0._wp
[13295]	237	DO_2D( 0, 0, 0, 0 )
[12377]	238	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * uu(ji ,jj,mbku(ji ,jj),Kmm) )**2 &
	239	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * uu(ji-1,jj,mbku(ji-1,jj),Kmm) )**2 &
	240	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vv(ji,jj ,mbkv(ji,jj ),Kmm) )**2 &
	241	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vv(ji,jj-1,mbkv(ji,jj-1),Kmm) )**2
	242	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
	243	!
	244	END_2D
[5463]	245	CALL iom_put( "taubot", z2d )
	246	ENDIF
[5107]	247
[12377]	248	CALL iom_put( "uoce", uu(:,:,:,Kmm) ) ! 3D i-current
	249	CALL iom_put( "ssu", uu(:,:,1,Kmm) ) ! surface i-current
[5107]	250	IF ( iom_use("sbu") ) THEN
[13472]	251	DO_2D( 0, 0, 0, 0 )
[12377]	252	ikbot = mbku(ji,jj)
	253	z2d(ji,jj) = uu(ji,jj,ikbot,Kmm)
	254	END_2D
[5107]	255	CALL iom_put( "sbu", z2d ) ! bottom i-current
	256	ENDIF
	257
[12377]	258	CALL iom_put( "voce", vv(:,:,:,Kmm) ) ! 3D j-current
	259	CALL iom_put( "ssv", vv(:,:,1,Kmm) ) ! surface j-current
[5107]	260	IF ( iom_use("sbv") ) THEN
[13472]	261	DO_2D( 0, 0, 0, 0 )
[12377]	262	ikbot = mbkv(ji,jj)
	263	z2d(ji,jj) = vv(ji,jj,ikbot,Kmm)
	264	END_2D
[5107]	265	CALL iom_put( "sbv", z2d ) ! bottom j-current
[4990]	266	ENDIF
[1482]	267
[14476]	268	! ! vertical velocity
[15087]	269	IF( ln_zad_Aimp ) THEN
	270	IF( iom_use('woce') ) THEN
	271	DO_3D( 0, 0, 0, 0, 1, jpk )
	272	z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
	273	END_3D
	274	CALL iom_put( "woce", z3d ) ! explicit plus implicit parts
	275	ENDIF
	276	ELSE
	277	CALL iom_put( "woce", ww )
[14476]	278	ENDIF
[13237]	279
[5461]	280	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
[14476]	281	! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[15087]	282	IF( ln_zad_Aimp ) THEN
	283	DO_3D( 0, 0, 0, 0, 1, jpk )
	284	z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ( ww(ji,jj,jk) + wi(ji,jj,jk) )
	285	END_3D
	286	ELSE
	287	DO_3D( 0, 0, 0, 0, 1, jpk )
	288	z3d(ji,jj,jk) = rho0 * e1e2t(ji,jj) * ww(ji,jj,jk)
	289	END_3D
	290	ENDIF
[5461]	291	CALL iom_put( "w_masstr" , z3d )
[14476]	292	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d * z3d )
[5461]	293	ENDIF
	294
[9019]	295	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	296	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	297	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
[5107]	298
[9019]	299	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	300	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
[6351]	301
[15017]	302	IF ( iom_use("sssgrad") .OR. iom_use("sssgrad2") ) THEN
	303	DO_2D( 0, 0, 0, 0 ) ! sss gradient
	304	zztmp = ts(ji,jj,1,jp_sal,Kmm)
	305	zztmpx = (ts(ji+1,jj,1,jp_sal,Kmm) - zztmp) * r1_e1u(ji,jj) + (zztmp - ts(ji-1,jj ,1,jp_sal,Kmm)) * r1_e1u(ji-1,jj)
	306	zztmpy = (ts(ji,jj+1,1,jp_sal,Kmm) - zztmp) * r1_e2v(ji,jj) + (zztmp - ts(ji ,jj-1,1,jp_sal,Kmm)) * r1_e2v(ji,jj-1)
	307	z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	308	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
	309	END_2D
	310	CALL iom_put( "sssgrad2", z2d ) ! square of module of sss gradient
	311	IF ( iom_use("sssgrad") ) THEN
	312	DO_2D( 0, 0, 0, 0 )
	313	z2d(ji,jj) = SQRT( z2d(ji,jj) )
	314	END_2D
	315	CALL iom_put( "sssgrad", z2d ) ! module of sss gradient
	316	ENDIF
[13472]	317	ENDIF
	318
[5107]	319	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[15017]	320	DO_2D( 0, 0, 0, 0 ) ! sst gradient
[12377]	321	zztmp = ts(ji,jj,1,jp_tem,Kmm)
	322	zztmpx = ( ts(ji+1,jj,1,jp_tem,Kmm) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - ts(ji-1,jj ,1,jp_tem,Kmm) ) * r1_e1u(ji-1,jj)
	323	zztmpy = ( ts(ji,jj+1,1,jp_tem,Kmm) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - ts(ji ,jj-1,1,jp_tem,Kmm) ) * r1_e2v(ji,jj-1)
[15017]	324	z2d(ji,jj) = 0.25_wp * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	325	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * vmask(ji,jj-1,1)
[12377]	326	END_2D
[9019]	327	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
[15017]	328	IF ( iom_use("sstgrad") ) THEN
	329	DO_2D( 0, 0, 0, 0 )
	330	z2d(ji,jj) = SQRT( z2d(ji,jj) )
	331	END_2D
	332	CALL iom_put( "sstgrad", z2d ) ! module of sst gradient
	333	ENDIF
[4990]	334	ENDIF
	335
[9019]	336	! heat and salt contents
[4990]	337	IF( iom_use("heatc") ) THEN
[7753]	338	z2d(:,:) = 0._wp
[13472]	339	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	340	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
	341	END_3D
[12489]	342	CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
[4990]	343	ENDIF
	344
	345	IF( iom_use("saltc") ) THEN
[7753]	346	z2d(:,:) = 0._wp
[13472]	347	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	348	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
	349	END_3D
[15017]	350	CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[4990]	351	ENDIF
[4840]	352	!
[13472]	353	IF( iom_use("salt2c") ) THEN
	354	z2d(:,:) = 0._wp
	355	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	356	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
	357	END_3D
[15017]	358	CALL iom_put( "salt2c", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[13472]	359	ENDIF
	360	!
	361	IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
[15087]	362	DO_3D( 0, 0, 0, 0, 1, jpk )
[15017]	363	zztmpx = uu(ji-1,jj ,jk,Kmm) + uu(ji,jj,jk,Kmm)
	364	zztmpy = vv(ji ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm)
	365	z3d(ji,jj,jk) = 0.25_wp * ( zztmpxzztmpx + zztmpyzztmpy )
[12377]	366	END_3D
[13472]	367	CALL iom_put( "ke", z3d ) ! kinetic energy
	368
	369	z2d(:,:) = 0._wp
	370	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	371	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
	372	END_3D
[15017]	373	CALL iom_put( "ke_int", z2d ) ! vertically integrated kinetic energy
[4990]	374	ENDIF
[6351]	375	!
[15017]	376	IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
[14053]	377	z2d(:,:) = 0._wp
	378	DO_2D( 0, 0, 0, 0 )
	379	z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
	380	& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
	381	& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
	382	& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
	383	& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
	384	END_2D
	385	IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d )
	386	ENDIF
	387	!
[15017]	388	IF ( iom_use("ssKEf") ) THEN ! surface kinetic energy at F point
	389	z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
[14053]	390	DO_2D( 0, 0, 0, 0 )
	391	z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
	392	& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
	393	& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
	394	& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
	395	& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
	396	END_2D
[14143]	397	CALL iom_put( "ssKEf", z2d )
[14053]	398	ENDIF
	399	!
[15017]	400	CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
[6351]	401	!
[7646]	402	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[15017]	403
[15087]	404	DO_3D( 0, 0, 0, 0, 1, jpk )
	405	z3d(ji,jj,jk) = rho0 * uu(ji,jj,jk,Kmm) * e2u(ji,jj) * e3u(ji,jj,jk,Kmm) * umask(ji,jj,jk)
	406	END_3D
[15017]	407	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	408
	409	IF( iom_use("u_masstr_vint") ) THEN
	410	z2d(:,:) = 0._wp
	411	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[15087]	412	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk)
[15017]	413	END_3D
	414	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
	415	ENDIF
	416	IF( iom_use("u_heattr") ) THEN
	417	z2d(:,:) = 0._wp
	418	zztmp = 0.5_wp * rcp
	419	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	420	z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji+1,jj,jk,jp_tem,Kmm) )
	421	END_3D
	422	CALL iom_put( "u_heattr", z2d ) ! heat transport in i-direction
	423	ENDIF
	424	IF( iom_use("u_salttr") ) THEN
	425	z2d(:,:) = 0._wp
	426	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	427	z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji+1,jj,jk,jp_sal,Kmm) )
	428	END_3D
	429	CALL iom_put( "u_salttr", z2d ) ! heat transport in i-direction
	430	ENDIF
	431
[4990]	432	ENDIF
	433
	434	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[15017]	435
[15087]	436	DO_3D( 0, 0, 0, 0, 1, jpk )
	437	z3d(ji,jj,jk) = rho0 * vv(ji,jj,jk,Kmm) * e1v(ji,jj) * e3v(ji,jj,jk,Kmm) * vmask(ji,jj,jk)
	438	END_3D
[15017]	439	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
	440
	441	IF( iom_use("v_heattr") ) THEN
	442	z2d(:,:) = 0._wp
	443	zztmp = 0.5_wp * rcp
	444	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	445	z2d(ji,jj) = z2d(ji,jj) + zztmp * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_tem,Kmm) + ts(ji,jj+1,jk,jp_tem,Kmm) )
	446	END_3D
	447	CALL iom_put( "v_heattr", z2d ) ! heat transport in j-direction
	448	ENDIF
	449	IF( iom_use("v_salttr") ) THEN
	450	z2d(:,:) = 0._wp
	451	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	452	z2d(ji,jj) = z2d(ji,jj) + 0.5 * z3d(ji,jj,jk) * ( ts(ji,jj,jk,jp_sal,Kmm) + ts(ji,jj+1,jk,jp_sal,Kmm) )
	453	END_3D
	454	CALL iom_put( "v_salttr", z2d ) ! heat transport in j-direction
	455	ENDIF
[4761]	456
[1756]	457	ENDIF
[7646]	458
	459	IF( iom_use("tosmint") ) THEN
[9019]	460	z2d(:,:) = 0._wp
[13295]	461	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[15017]	462	z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
[12377]	463	END_3D
[15017]	464	CALL iom_put( "tosmint", z2d ) ! Vertical integral of temperature
[7646]	465	ENDIF
	466	IF( iom_use("somint") ) THEN
[15017]	467	z2d(:,:) = 0._wp
[13295]	468	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[15017]	469	z2d(ji,jj) = z2d(ji,jj) + rho0 * e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
[12377]	470	END_3D
[15017]	471	CALL iom_put( "somint", z2d ) ! Vertical integral of salinity
[7646]	472	ENDIF
	473
[15017]	474	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
[12377]	475
[15017]	476	IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging
[14053]	477
[14143]	478	! Output of surface vorticity terms
[15017]	479	!
	480	CALL iom_put( "ssplavor", ff_f ) ! planetary vorticity ( f )
	481	!
	482	IF ( iom_use("ssrelvor") .OR. iom_use("ssEns") .OR. &
	483	& iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
[14053]	484	!
	485	z2d(:,:) = 0._wp
[15017]	486	DO_2D( 0, 0, 0, 0 )
[14053]	487	z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
	488	& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
	489	END_2D
[15017]	490	CALL iom_put( "ssrelvor", z2d ) ! relative vorticity ( zeta )
[14053]	491	!
[15017]	492	IF ( iom_use("ssEns") .OR. iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") ) THEN
	493	DO_2D( 0, 0, 0, 0 )
	494	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	495	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	496	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	497	ELSE ; ze3 = 0._wp
	498	ENDIF
	499	z2d(ji,jj) = ze3 * z2d(ji,jj)
	500	END_2D
	501	CALL iom_put( "ssrelpotvor", z2d ) ! relative potential vorticity (zeta/h)
	502	!
	503	IF ( iom_use("ssEns") .OR. iom_use("ssabspotvor") ) THEN
	504	DO_2D( 0, 0, 0, 0 )
	505	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	506	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	507	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	508	ELSE ; ze3 = 0._wp
	509	ENDIF
	510	z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)
	511	END_2D
	512	CALL iom_put( "ssabspotvor", z2d ) ! absolute potential vorticity ( q )
	513	!
	514	IF ( iom_use("ssEns") ) THEN
	515	DO_2D( 0, 0, 0, 0 )
	516	z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj)
	517	END_2D
	518	CALL iom_put( "ssEns", z2d ) ! potential enstrophy ( 1/2*q2 )
	519	ENDIF
[14053]	520	ENDIF
[15017]	521	ENDIF
[14053]	522	ENDIF
[6140]	523
[9124]	524	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	525	!
[1482]	526	END SUBROUTINE dia_wri
	527
	528	#else
[2528]	529	!!----------------------------------------------------------------------
	530	!! Default option use IOIPSL library
	531	!!----------------------------------------------------------------------
	532
[9652]	533	INTEGER FUNCTION dia_wri_alloc()
	534	!!----------------------------------------------------------------------
	535	INTEGER, DIMENSION(2) :: ierr
	536	!!----------------------------------------------------------------------
[12493]	537	IF( nn_write == -1 ) THEN
	538	dia_wri_alloc = 0
	539	ELSE
	540	ierr = 0
	541	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	542	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	543	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
[9652]	544	!
[12493]	545	dia_wri_alloc = MAXVAL(ierr)
	546	CALL mpp_sum( 'diawri', dia_wri_alloc )
	547	!
	548	ENDIF
[9652]	549	!
	550	END FUNCTION dia_wri_alloc
[12377]	551
	552	INTEGER FUNCTION dia_wri_alloc_abl()
	553	!!----------------------------------------------------------------------
	554	ALLOCATE( ndex_hA(jpijpj), ndex_A (jpijpj*jpkam1), STAT=dia_wri_alloc_abl)
	555	CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
	556	!
	557	END FUNCTION dia_wri_alloc_abl
[9652]	558
	559
[12377]	560	SUBROUTINE dia_wri( kt, Kmm )
[3]	561	!!---------------------------------------------------------------------
	562	!! * ROUTINE dia_wri *
	563	!!
	564	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	565	!! NETCDF format is used by default
	566	!!
	567	!! ** Method : At the beginning of the first time step (nit000),
	568	!! define all the NETCDF files and fields
	569	!! At each time step call histdef to compute the mean if ncessary
[11536]	570	!! Each nn_write time step, output the instantaneous or mean fields
[3]	571	!!----------------------------------------------------------------------
[5836]	572	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	573	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[5836]	574	!
[2528]	575	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	576	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	577	INTEGER :: inum = 11 ! temporary logical unit
[3294]	578	INTEGER :: ji, jj, jk ! dummy loop indices
	579	INTEGER :: ierr ! error code return from allocation
[2528]	580	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[12377]	581	INTEGER :: ipka ! ABL
[3609]	582	INTEGER :: jn, ierror ! local integers
[6140]	583	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	584	!
[15141]	585	REAL(wp), DIMENSION(jpi,jpj ) :: z2d ! 2D workspace
	586	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
[12377]	587	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
[3]	588	!!----------------------------------------------------------------------
[1482]	589	!
[9019]	590	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[12377]	591	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	592	ninist = 0
	593	ENDIF
	594	!
[11536]	595	IF( nn_write == -1 ) RETURN ! we will never do any output
	596	!
	597	IF( ln_timing ) CALL timing_start('dia_wri')
	598	!
[3]	599	! 0. Initialisation
	600	! -----------------
[632]	601
[9019]	602	ll_print = .FALSE. ! local variable for debugging
[3]	603	ll_print = ll_print .AND. lwp
	604
	605	! Define frequency of output and means
[5566]	606	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	607	#if defined key_diainstant
[12489]	608	zsto = nn_write * rn_Dt
[1312]	609	clop = "inst("//TRIM(clop)//")"
[3]	610	#else
[12489]	611	zsto=rn_Dt
[1312]	612	clop = "ave("//TRIM(clop)//")"
[3]	613	#endif
[12489]	614	zout = nn_write * rn_Dt
	615	zmax = ( nitend - nit000 + 1 ) * rn_Dt
[3]	616
	617	! Define indices of the horizontal output zoom and vertical limit storage
[13286]	618	iimi = Nis0 ; iima = Nie0
	619	ijmi = Njs0 ; ijma = Nje0
[3]	620	ipk = jpk
[12377]	621	IF(ln_abl) ipka = jpkam1
[3]	622
	623	! define time axis
[1334]	624	it = kt
	625	itmod = kt - nit000 + 1
[3]	626
	627	! 1. Define NETCDF files and fields at beginning of first time step
	628	! -----------------------------------------------------------------
	629
	630	IF( kt == nit000 ) THEN
	631
	632	! Define the NETCDF files (one per grid)
[632]	633
[3]	634	! Compute julian date from starting date of the run
[12489]	635	CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
[1309]	636	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	637	IF(lwp)WRITE(numout,*)
	638	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	639	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	640	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	641	' limit storage in depth = ', ipk
	642
	643	! WRITE root name in date.file for use by postpro
[1581]	644	IF(lwp) THEN
[11536]	645	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	646	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	647	WRITE(inum,*) clhstnam
	648	CLOSE(inum)
	649	ENDIF
[632]	650
[3]	651	! Define the T grid FILE ( nid_T )
[632]	652
[11536]	653	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	654	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	655	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	656	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	657	& nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	658	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	659	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	660	! ! Index of ocean points
	661	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	662	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	663	!
	664	IF( ln_icebergs ) THEN
	665	!
	666	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	667	!! that routine is called from nemogcm, so do it here immediately before its needed
	668	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	669	CALL mpp_sum( 'diawri', ierror )
[3609]	670	IF( ierror /= 0 ) THEN
	671	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	672	RETURN
	673	ENDIF
	674	!
	675	!! iceberg vertical coordinate is class number
	676	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	677	& "number", nclasses, class_num, nb_T )
	678	!
	679	!! each class just needs the surface index pattern
	680	ndim_bT = 3
	681	DO jn = 1,nclasses
	682	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	683	ENDDO
	684	!
	685	ENDIF
[3]	686
	687	! Define the U grid FILE ( nid_U )
	688
[11536]	689	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	690	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	691	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	692	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	693	& nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	694	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	695	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	696	! ! Index of ocean points
	697	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	698	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	699
	700	! Define the V grid FILE ( nid_V )
	701
[11536]	702	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	703	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	704	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	705	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	706	& nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	707	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	708	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	709	! ! Index of ocean points
	710	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	711	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	712
	713	! Define the W grid FILE ( nid_W )
	714
[11536]	715	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	716	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	717	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	718	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	719	& nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	720	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	721	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	722
[12377]	723	IF( ln_abl ) THEN
	724	! Define the ABL grid FILE ( nid_A )
	725	CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
	726	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	727	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	728	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	729	& nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
[12377]	730	CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept
	731	& "m", ipka, ght_abl(2:jpka), nz_A, "up" )
	732	! ! Index of ocean points
	733	ALLOCATE( zw3d_abl(jpi,jpj,ipka) )
	734	zw3d_abl(:,:,:) = 1._wp
	735	CALL wheneq( jpijpjipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume
	736	CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface
	737	DEALLOCATE(zw3d_abl)
	738	ENDIF
[13237]	739	!
[632]	740
[3]	741	! Declare all the output fields as NETCDF variables
	742
	743	! !!! nid_T : 3D
	744	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	745	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	746	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	747	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	748	IF( .NOT.ln_linssh ) THEN
[13237]	749	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t n
[4292]	750	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13237]	751	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t n
[4292]	752	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13237]	753	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t n
[4292]	754	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	755	ENDIF
[3]	756	! !!! nid_T : 2D
	757	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	758	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	759	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	760	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	761	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	762	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	763	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	764	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	765	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	766	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	767	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	768	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	769	IF( ln_linssh ) THEN
[12377]	770	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts(:,:,1,jp_tem,Kmm)
[3625]	771	& , "KgC/m2/s", & ! sosst_cd
[4292]	772	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	773	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts(:,:,1,jp_sal,Kmm)
[3625]	774	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	775	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	776	ENDIF
[888]	777	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	778	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	779	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	780	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[15136]	781	IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
	782	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	783	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	784	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	785	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	786	ENDIF
[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
[15141]	939	DO_3D( 0, 0, 0, 0, 1, jpk )
	940	z3d(ji,jj,jk) = ts(ji,jj,jk,jp_tem,Kmm) * e3t(ji,jj,jk,Kmm)
	941	END_3D
	942	CALL histwrite( nid_T, "votemper", it, z3d, ndim_T , ndex_T ) ! heat content
	943	DO_3D( 0, 0, 0, 0, 1, jpk )
	944	z3d(ji,jj,jk) = ts(ji,jj,jk,jp_sal,Kmm) * e3t(ji,jj,jk,Kmm)
	945	END_3D
	946	CALL histwrite( nid_T, "vosaline", it, z3d, ndim_T , ndex_T ) ! salt content
	947	DO_2D( 0, 0, 0, 0 )
	948	z2d(ji,jj ) = ts(ji,jj, 1,jp_tem,Kmm) * e3t(ji,jj, 1,Kmm)
	949	END_2D
	950	CALL histwrite( nid_T, "sosstsst", it, z2d, ndim_hT, ndex_hT ) ! sea surface heat content
	951	DO_2D( 0, 0, 0, 0 )
	952	z2d(ji,jj ) = ts(ji,jj, 1,jp_sal,Kmm) * e3t(ji,jj, 1,Kmm)
	953	END_2D
	954	CALL histwrite( nid_T, "sosaline", it, z2d, ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	955	ELSE
[12377]	956	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
	957	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
	958	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT ) ! sea surface temperature
	959	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity
[4292]	960	ENDIF
[6140]	961	IF( .NOT.ln_linssh ) THEN
[15141]	962	DO_3D( 0, 0, 0, 0, 1, jpk )
	963	z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
	964	END_3D
	965	CALL histwrite( nid_T, "vovvle3t", it, z3d , ndim_T , ndex_T ) ! level thickness
	966	DO_3D( 0, 0, 0, 0, 1, jpk )
	967	z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
	968	END_3D
	969	CALL histwrite( nid_T, "vovvldep", it, z3d , ndim_T , ndex_T ) ! t-point depth
	970	DO_3D( 0, 0, 0, 0, 1, jpk )
	971	z3d(ji,jj,jk) = ( ( e3t(ji,jj,jk,Kmm) - e3t_0(ji,jj,jk) ) / e3t_0(ji,jj,jk) * 100._wp * tmask(ji,jj,jk) ) ** 2
	972	END_3D
	973	CALL histwrite( nid_T, "vovvldef", it, z3d , ndim_T , ndex_T ) ! level thickness deformation
[4292]	974	ENDIF
[15141]	975	CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
	976	DO_2D( 0, 0, 0, 0 )
	977	z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
	978	END_2D
	979	CALL histwrite( nid_T, "sowaflup", it, z2d , ndim_hT, ndex_hT ) ! upward water flux
[4570]	980	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	981	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	982	! (includes virtual salt flux beneath ice
	983	! in linear free surface case)
[6140]	984	IF( ln_linssh ) THEN
[15141]	985	DO_2D( 0, 0, 0, 0 )
	986	z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_tem,Kmm)
	987	END_2D
	988	CALL histwrite( nid_T, "sosst_cd", it, z2d, ndim_hT, ndex_hT ) ! c/d term on sst
	989	DO_2D( 0, 0, 0, 0 )
	990	z2d(ji,jj) = emp (ji,jj) * ts(ji,jj,1,jp_sal,Kmm)
	991	END_2D
	992	CALL histwrite( nid_T, "sosss_cd", it, z2d, ndim_hT, ndex_hT ) ! c/d term on sss
[4292]	993	ENDIF
[15141]	994	DO_2D( 0, 0, 0, 0 )
	995	z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
	996	END_2D
	997	CALL histwrite( nid_T, "sohefldo", it, z2d , ndim_hT, ndex_hT ) ! total heat flux
[3]	998	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[15136]	999	IF( ALLOCATED(hmld) ) THEN ! zdf_mxl not called by SWE
	1000	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
	1001	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
	1002	ENDIF
[1037]	1003	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	1004	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[12377]	1005	!
	1006	IF( ln_abl ) THEN
	1007	ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
	1008	IF( ln_mskland ) THEN
	1009	DO jk=1,jpka
	1010	zw3d_abl(:,:,jk) = tmask(:,:,1)
	1011	END DO
	1012	ELSE
	1013	zw3d_abl(:,:,:) = 1._wp
	1014	ENDIF
	1015	CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh
	1016	CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl
	1017	CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl
	1018	CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl
	1019	CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl
	1020	CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl
	1021	CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl
	1022	CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl
	1023	#if defined key_si3
	1024	CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i
	1025	#endif
	1026	DEALLOCATE(zw3d_abl)
	1027	ENDIF
	1028	!
[3609]	1029	IF( ln_icebergs ) THEN
	1030	!
	1031	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	1032	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	1033	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	1034	!
	1035	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	1036	!
	1037	IF( ln_bergdia ) THEN
	1038	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	1039	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	1040	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	1041	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	1042	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	1043	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	1044	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	1045	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	1046	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	1047	!
	1048	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	1049	ENDIF
	1050	ENDIF
	1051
[11536]	1052	IF( ln_ssr ) THEN
[4990]	1053	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	1054	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[15141]	1055	DO_2D( 0, 0, 0, 0 )
	1056	z2d(ji,jj) = erp(ji,jj) * ts(ji,jj,1,jp_sal,Kmm) * tmask(ji,jj,1)
	1057	END_2D
	1058	CALL histwrite( nid_T, "sosafldp", it, z2d , ndim_hT, ndex_hT ) ! salt flux damping
[4990]	1059	ENDIF
	1060	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	1061	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	1062
	1063	#if defined key_diahth
	1064	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	1065	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	1066	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	1067	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	1068	#endif
[888]	1069
[15141]	1070	CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
[888]	1071	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	1072
[15141]	1073	CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
[888]	1074	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	1075
[11418]	1076	IF( ln_zad_Aimp ) THEN
[15141]	1077	DO_3D( 0, 0, 0, 0, 1, jpk )
	1078	z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
	1079	END_3D
	1080	CALL histwrite( nid_W, "vovecrtz", it, z3d , ndim_T, ndex_T ) ! vert. current
[11418]	1081	ELSE
[12377]	1082	CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current
[11418]	1083	ENDIF
[3]	1084	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	1085	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	1086	IF( ln_zdfddm ) THEN
	1087	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	1088	ENDIF
	1089
[7646]	1090	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	1091	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	1092	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	1093	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	1094	ENDIF
	1095
[1318]	1096	! 3. Close all files
	1097	! ---------------------------------------
[1561]	1098	IF( kt == nitend ) THEN
[3]	1099	CALL histclo( nid_T )
	1100	CALL histclo( nid_U )
	1101	CALL histclo( nid_V )
	1102	CALL histclo( nid_W )
[12377]	1103	IF(ln_abl) CALL histclo( nid_A )
[3]	1104	ENDIF
[2528]	1105	!
[9124]	1106	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	1107	!
[3]	1108	END SUBROUTINE dia_wri
[1567]	1109	#endif
	1110
[12377]	1111	SUBROUTINE dia_wri_state( Kmm, cdfile_name )
[3]	1112	!!---------------------------------------------------------------------
	1113	!! * ROUTINE dia_wri_state *
	1114	!!
	1115	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	1116	!! the instantaneous ocean state and forcing fields.
	1117	!! Used to find errors in the initial state or save the last
	1118	!! ocean state in case of abnormal end of a simulation
	1119	!!
	1120	!! ** Method : NetCDF files using ioipsl
	1121	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	1122	!! File 'output.abort.nc' is created in case of abnormal job end
	1123	!!----------------------------------------------------------------------
[12377]	1124	INTEGER , INTENT( in ) :: Kmm ! time level index
[1334]	1125	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	1126	!!
[15141]	1127	INTEGER :: ji, jj, jk ! dummy loop indices
	1128	INTEGER :: inum
	1129	REAL(wp), DIMENSION(jpi,jpj) :: z2d
	1130	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d
[3]	1131	!!----------------------------------------------------------------------
[3294]	1132	!
[13237]	1133	IF(lwp) THEN
	1134	WRITE(numout,*)
	1135	WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	1136	WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
	1137	WRITE(numout,*) ' and named :', cdfile_name, '...nc'
	1138	ENDIF
[12649]	1139	!
	1140	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
	1141	!
[12377]	1142	CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature
	1143	CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity
[15141]	1144	CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm) ) ! sea surface height
	1145	CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm) ) ! now i-velocity
	1146	CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity
[11418]	1147	IF( ln_zad_Aimp ) THEN
[15141]	1148	DO_3D( 0, 0, 0, 0, 1, jpk )
	1149	z3d(ji,jj,jk) = ww(ji,jj,jk) + wi(ji,jj,jk)
	1150	END_3D
	1151	CALL iom_rstput( 0, 0, inum, 'vovecrtz', z3d ) ! now k-velocity
[11418]	1152	ELSE
[12377]	1153	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity
[11418]	1154	ENDIF
[15141]	1155	CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep )
[13237]	1156	CALL iom_rstput( 0, 0, inum, 'ht' , ht(:,:) ) ! now water column height
[12649]	1157	!
[12377]	1158	IF ( ln_isf ) THEN
	1159	IF (ln_isfcav_mlt) THEN
[15141]	1160	CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav )
	1161	CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav )
	1162	CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav )
	1163	CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,wp) )
	1164	CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,wp) )
[12649]	1165	CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,wp), ktype = jp_i1 )
[12377]	1166	END IF
	1167	IF (ln_isfpar_mlt) THEN
[15141]	1168	CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,wp) )
	1169	CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par )
	1170	CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par )
	1171	CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par )
	1172	CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,wp) )
	1173	CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,wp) )
[12649]	1174	CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,wp), ktype = jp_i1 )
[12377]	1175	END IF
	1176	END IF
[12649]	1177	!
[7646]	1178	IF( ALLOCATED(ahtu) ) THEN
[10425]	1179	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	1180	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	1181	ENDIF
	1182	IF( ALLOCATED(ahmt) ) THEN
[10425]	1183	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	1184	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	1185	ENDIF
[15141]	1186	DO_2D( 0, 0, 0, 0 )
	1187	z2d(ji,jj) = emp(ji,jj) - rnf(ji,jj)
	1188	END_2D
	1189	CALL iom_rstput( 0, 0, inum, 'sowaflup', z2d ) ! freshwater budget
	1190	DO_2D( 0, 0, 0, 0 )
	1191	z2d(ji,jj) = qsr(ji,jj) + qns(ji,jj)
	1192	END_2D
	1193	CALL iom_rstput( 0, 0, inum, 'sohefldo', z2d ) ! total heat flux
[10425]	1194	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	1195	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	1196	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	1197	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[15141]	1198	IF( .NOT.ln_linssh ) THEN
	1199	DO_3D( 0, 0, 0, 0, 1, jpk )
	1200	z3d(ji,jj,jk) = gdept(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
	1201	END_3D
	1202	CALL iom_rstput( 0, 0, inum, 'vovvldep', z3d ) ! T-cell depth
	1203	DO_3D( 0, 0, 0, 0, 1, jpk )
	1204	z3d(ji,jj,jk) = e3t(ji,jj,jk,Kmm) ! 3D workspace for qco substitution
	1205	END_3D
	1206	CALL iom_rstput( 0, 0, inum, 'vovvle3t', z3d ) ! T-cell thickness
[10425]	1207	END IF
[7646]	1208	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	1209	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	1210	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	1211	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	1212	ENDIF
[12377]	1213	IF ( ln_abl ) THEN
	1214	CALL iom_rstput ( 0, 0, inum, "uz1_abl", u_abl(:,:,2,nt_a ) ) ! now first level i-wind
	1215	CALL iom_rstput ( 0, 0, inum, "vz1_abl", v_abl(:,:,2,nt_a ) ) ! now first level j-wind
	1216	CALL iom_rstput ( 0, 0, inum, "tz1_abl", tq_abl(:,:,2,nt_a,1) ) ! now first level temperature
	1217	CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity
	1218	ENDIF
[14045]	1219	IF( ln_zdfosm ) THEN
	1220	CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth
	1221	CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth
	1222	CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion
	1223	CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity
	1224	CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing
	1225	CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing
	1226	CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask ) ! non-local u forcing
	1227	CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask ) ! non-local v forcing
	1228	IF( ln_osm_mle ) THEN
	1229	CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth
	1230	END IF
	1231	ENDIF
[12649]	1232	!
	1233	CALL iom_close( inum )
	1234	!
[10425]	1235	#if defined key_si3
	1236	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
[12649]	1237	CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' )
[10425]	1238	CALL ice_wri_state( inum )
[12649]	1239	CALL iom_close( inum )
[1561]	1240	ENDIF
[12933]	1241	!
[1561]	1242	#endif
[3]	1243	END SUBROUTINE dia_wri_state
[6140]	1244
[3]	1245	!!======================================================================
	1246	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: NEMO/trunk/src/OCE/DIA/diawri.F90 @ 15199

Download in other formats: