Context Navigation

← Previous Revision
Next Revision →
Normal
Revision Log

diawri.F90

Last change on this file was 15540, checked in by sparonuz, 3 years ago
Mixed precision version, tested up to 30 years on ORCA2.
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
[14644]	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
[14986]	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
[2528]	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
[14644]	97	#if defined key_xios
[3]	98	!!----------------------------------------------------------------------
[14644]	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
[15540]	125	REAL(dp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
[9019]	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
[14200]	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
[14644]	253	! ! vertical velocity
	254	IF( ln_zad_Aimp ) THEN ; CALL iom_put( "woce", ww + wi ) ! explicit plus implicit parts
	255	ELSE ; CALL iom_put( "woce", ww )
	256	ENDIF
[13237]	257
[5461]	258	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
[14644]	259	! ! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
[5461]	260	DO jk = 1, jpk
[14644]	261	IF( ln_zad_Aimp ) THEN
	262	z3d(:,:,jk) = rho0 * e1e2t(:,:) * ( ww(:,:,jk) + wi(:,:,jk) )
	263	ELSE
	264	z3d(:,:,jk) = rho0 * e1e2t(:,:) * ww(:,:,jk)
	265	ENDIF
[5461]	266	END DO
	267	CALL iom_put( "w_masstr" , z3d )
[14644]	268	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d * z3d )
[5461]	269	ENDIF
	270
[9019]	271	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
	272	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
	273	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
[5107]	274
[9019]	275	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
	276	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
[6351]	277
[13472]	278	IF ( iom_use("socegrad") .OR. iom_use("socegrad2") ) THEN
	279	z3d(:,:,jpk) = 0.
	280	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	281	zztmp = ts(ji,jj,jk,jp_sal,Kmm)
	282	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)
	283	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)
	284	z3d(ji,jj,jk) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	285	& * umask(ji,jj,jk) * umask(ji-1,jj,jk) * vmask(ji,jj,jk) * umask(ji,jj-1,jk)
	286	END_3D
	287	CALL iom_put( "socegrad2", z3d ) ! square of module of sal gradient
	288	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	289	z3d(ji,jj,jk) = SQRT( z3d(ji,jj,jk) )
	290	END_3D
	291	CALL iom_put( "socegrad" , z3d ) ! module of sal gradient
	292	ENDIF
	293
[5107]	294	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
[13497]	295	DO_2D( 0, 0, 0, 0 ) ! sst gradient
[12377]	296	zztmp = ts(ji,jj,1,jp_tem,Kmm)
	297	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)
	298	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)
	299	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
	300	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
	301	END_2D
[9019]	302	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
[13472]	303	DO_2D( 0, 0, 0, 0 )
	304	z2d(ji,jj) = SQRT( z2d(ji,jj) )
	305	END_2D
[9019]	306	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
[4990]	307	ENDIF
	308
[9019]	309	! heat and salt contents
[4990]	310	IF( iom_use("heatc") ) THEN
[7753]	311	z2d(:,:) = 0._wp
[13472]	312	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	313	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm) * tmask(ji,jj,jk)
	314	END_3D
[12489]	315	CALL iom_put( "heatc", rho0_rcp * z2d ) ! vertically integrated heat content (J/m2)
[4990]	316	ENDIF
	317
	318	IF( iom_use("saltc") ) THEN
[7753]	319	z2d(:,:) = 0._wp
[13472]	320	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	321	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm) * tmask(ji,jj,jk)
	322	END_3D
[12489]	323	CALL iom_put( "saltc", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
[4990]	324	ENDIF
[4840]	325	!
[13472]	326	IF( iom_use("salt2c") ) THEN
	327	z2d(:,:) = 0._wp
	328	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	329	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)
	330	END_3D
	331	CALL iom_put( "salt2c", rho0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
	332	ENDIF
	333	!
	334	IF ( iom_use("ke") .OR. iom_use("ke_int") ) THEN
[9399]	335	z3d(:,:,jpk) = 0._wp
[13295]	336	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[13472]	337	zztmpx = 0.5 * ( uu(ji-1,jj ,jk,Kmm) + uu(ji,jj,jk,Kmm) )
	338	zztmpy = 0.5 * ( vv(ji ,jj-1,jk,Kmm) + vv(ji,jj,jk,Kmm) )
	339	z3d(ji,jj,jk) = 0.5 * ( zztmpxzztmpx + zztmpyzztmpy )
[12377]	340	END_3D
[13472]	341	CALL iom_put( "ke", z3d ) ! kinetic energy
	342
	343	z2d(:,:) = 0._wp
	344	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	345	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * z3d(ji,jj,jk) * e1e2t(ji,jj) * tmask(ji,jj,jk)
	346	END_3D
	347	CALL iom_put( "ke_int", z2d ) ! vertically integrated kinetic energy
[4990]	348	ENDIF
[6351]	349	!
[14053]	350	IF ( iom_use("sKE") ) THEN ! surface kinetic energy at T point
	351	z2d(:,:) = 0._wp
	352	DO_2D( 0, 0, 0, 0 )
	353	z2d(ji,jj) = 0.25_wp * ( uu(ji ,jj,1,Kmm) * uu(ji ,jj,1,Kmm) * e1e2u(ji ,jj) * e3u(ji ,jj,1,Kmm) &
	354	& + uu(ji-1,jj,1,Kmm) * uu(ji-1,jj,1,Kmm) * e1e2u(ji-1,jj) * e3u(ji-1,jj,1,Kmm) &
	355	& + vv(ji,jj ,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji,jj ) * e3v(ji,jj ,1,Kmm) &
	356	& + vv(ji,jj-1,1,Kmm) * vv(ji,jj-1,1,Kmm) * e1e2v(ji,jj-1) * e3v(ji,jj-1,1,Kmm) ) &
	357	& * r1_e1e2t(ji,jj) / e3t(ji,jj,1,Kmm) * ssmask(ji,jj)
	358	END_2D
[14286]	359	CALL lbc_lnk( 'diawri', z2d, 'T', 1._wp )
[14053]	360	IF ( iom_use("sKE" ) ) CALL iom_put( "sKE" , z2d )
	361	ENDIF
	362	!
[14200]	363	IF ( iom_use("ssKEf") ) THEN ! surface kinetic energy at F point
[14053]	364	z2d(:,:) = 0._wp ! CAUTION : only valid in SWE, not with bathymetry
	365	DO_2D( 0, 0, 0, 0 )
	366	z2d(ji,jj) = 0.25_wp * ( uu(ji,jj ,1,Kmm) * uu(ji,jj ,1,Kmm) * e1e2u(ji,jj ) * e3u(ji,jj ,1,Kmm) &
	367	& + uu(ji,jj+1,1,Kmm) * uu(ji,jj+1,1,Kmm) * e1e2u(ji,jj+1) * e3u(ji,jj+1,1,Kmm) &
	368	& + vv(ji ,jj,1,Kmm) * vv(ji,jj ,1,Kmm) * e1e2v(ji ,jj) * e3v(ji ,jj,1,Kmm) &
	369	& + vv(ji+1,jj,1,Kmm) * vv(ji+1,jj,1,Kmm) * e1e2v(ji+1,jj) * e3v(ji+1,jj,1,Kmm) ) &
	370	& * r1_e1e2f(ji,jj) / e3f(ji,jj,1) * ssfmask(ji,jj)
	371	END_2D
[14286]	372	CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp )
[14200]	373	CALL iom_put( "ssKEf", z2d )
[14053]	374	ENDIF
	375	!
[12377]	376	CALL iom_put( "hdiv", hdiv ) ! Horizontal divergence
[13472]	377
	378	IF ( iom_use("relvor") .OR. iom_use("absvor") .OR. iom_use("potvor") ) THEN
	379
	380	z3d(:,:,jpk) = 0._wp
	381	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	382	z3d(ji,jj,jk) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,jk,Kmm) - e2v(ji,jj) * vv(ji,jj,jk,Kmm) &
	383	& - e1u(ji ,jj+1) * uu(ji ,jj+1,jk,Kmm) + e1u(ji,jj) * uu(ji,jj,jk,Kmm) ) * r1_e1e2f(ji,jj)
	384	END_3D
	385	CALL iom_put( "relvor", z3d ) ! relative vorticity
	386
	387	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	388	z3d(ji,jj,jk) = ff_f(ji,jj) + z3d(ji,jj,jk)
	389	END_3D
	390	CALL iom_put( "absvor", z3d ) ! absolute vorticity
	391
	392	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
	393	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) &
	394	& + e3t(ji,jj ,jk,Kmm)tmask(ji,jj ,jk) + e3t(ji+1,jj ,jk,Kmm)tmask(ji+1,jj ,jk) )
	395	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	396	ELSE ; ze3 = 0._wp
	397	ENDIF
	398	z3d(ji,jj,jk) = ze3 * z3d(ji,jj,jk)
	399	END_3D
	400	CALL iom_put( "potvor", z3d ) ! potential vorticity
	401
	402	ENDIF
[6351]	403	!
[7646]	404	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
[7753]	405	z3d(:,:,jpk) = 0.e0
	406	z2d(:,:) = 0.e0
[1756]	407	DO jk = 1, jpkm1
[12489]	408	z3d(:,:,jk) = rho0 * uu(:,:,jk,Kmm) * e2u(:,:) * e3u(:,:,jk,Kmm) * umask(:,:,jk)
[7753]	409	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
[1756]	410	END DO
[9019]	411	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
	412	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
[4990]	413	ENDIF
	414
	415	IF( iom_use("u_heattr") ) THEN
[9019]	416	z2d(:,:) = 0._wp
[13295]	417	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	418	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) )
	419	END_3D
[9019]	420	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
[4990]	421	ENDIF
[4761]	422
[4990]	423	IF( iom_use("u_salttr") ) THEN
[7753]	424	z2d(:,:) = 0.e0
[13295]	425	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	426	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) )
	427	END_3D
[9019]	428	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
[4990]	429	ENDIF
[4761]	430
[4990]	431
	432	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
[7753]	433	z3d(:,:,jpk) = 0.e0
[1756]	434	DO jk = 1, jpkm1
[12489]	435	z3d(:,:,jk) = rho0 * vv(:,:,jk,Kmm) * e1v(:,:) * e3v(:,:,jk,Kmm) * vmask(:,:,jk)
[1756]	436	END DO
[9019]	437	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
[4990]	438	ENDIF
	439
	440	IF( iom_use("v_heattr") ) THEN
[7753]	441	z2d(:,:) = 0.e0
[13295]	442	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	443	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) )
	444	END_3D
[9019]	445	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
[4990]	446	ENDIF
[4761]	447
[4990]	448	IF( iom_use("v_salttr") ) THEN
[9019]	449	z2d(:,:) = 0._wp
[13295]	450	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	451	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) )
	452	END_3D
[9019]	453	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
[1756]	454	ENDIF
[7646]	455
	456	IF( iom_use("tosmint") ) THEN
[9019]	457	z2d(:,:) = 0._wp
[13295]	458	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	459	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_tem,Kmm)
	460	END_3D
[12489]	461	CALL iom_put( "tosmint", rho0 * z2d ) ! Vertical integral of temperature
[7646]	462	ENDIF
	463	IF( iom_use("somint") ) THEN
[7753]	464	z2d(:,:)=0._wp
[13295]	465	DO_3D( 0, 0, 0, 0, 1, jpkm1 )
[12377]	466	z2d(ji,jj) = z2d(ji,jj) + e3t(ji,jj,jk,Kmm) * ts(ji,jj,jk,jp_sal,Kmm)
	467	END_3D
[12489]	468	CALL iom_put( "somint", rho0 * z2d ) ! Vertical integral of salinity
[7646]	469	ENDIF
	470
[9019]	471	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
[2528]	472	!
[12377]	473
	474	IF (ln_dia25h) CALL dia_25h( kt, Kmm ) ! 25h averaging
[14053]	475
[14200]	476	! Output of surface vorticity terms
	477	IF ( iom_use("ssrelvor") .OR. iom_use("ssplavor") .OR. &
	478	& iom_use("ssrelpotvor") .OR. iom_use("ssabspotvor") .OR. &
	479	& iom_use("ssEns") ) THEN
[14053]	480	!
	481	z2d(:,:) = 0._wp
	482	ze3 = 0._wp
	483	DO_2D( 1, 0, 1, 0 )
	484	z2d(ji,jj) = ( e2v(ji+1,jj ) * vv(ji+1,jj ,1,Kmm) - e2v(ji,jj) * vv(ji,jj,1,Kmm) &
	485	& - e1u(ji ,jj+1) * uu(ji ,jj+1,1,Kmm) + e1u(ji,jj) * uu(ji,jj,1,Kmm) ) * r1_e1e2f(ji,jj)
	486	END_2D
[14286]	487	CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp )
[14200]	488	CALL iom_put( "ssrelvor", z2d ) ! relative vorticity ( zeta )
[14053]	489	!
[14200]	490	CALL iom_put( "ssplavor", ff_f ) ! planetary vorticity ( f )
[14053]	491	!
	492	DO_2D( 1, 0, 1, 0 )
	493	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	494	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	495	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	496	ELSE ; ze3 = 0._wp
	497	ENDIF
	498	z2d(ji,jj) = ze3 * z2d(ji,jj)
	499	END_2D
[14286]	500	CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp )
[14200]	501	CALL iom_put( "ssrelpotvor", z2d ) ! relative potential vorticity (zeta/h)
[14053]	502	!
	503	DO_2D( 1, 0, 1, 0 )
	504	ze3 = ( e3t(ji,jj+1,1,Kmm) * e1e2t(ji,jj+1) + e3t(ji+1,jj+1,1,Kmm) * e1e2t(ji+1,jj+1) &
	505	& + e3t(ji,jj ,1,Kmm) * e1e2t(ji,jj ) + e3t(ji+1,jj ,1,Kmm) * e1e2t(ji+1,jj ) ) * r1_e1e2f(ji,jj)
	506	IF( ze3 /= 0._wp ) THEN ; ze3 = 4._wp / ze3
	507	ELSE ; ze3 = 0._wp
	508	ENDIF
	509	z2d(ji,jj) = ze3 * ff_f(ji,jj) + z2d(ji,jj)
	510	END_2D
[14286]	511	CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp )
[14200]	512	CALL iom_put( "ssabspotvor", z2d ) ! absolute potential vorticity ( q )
[14053]	513	!
	514	DO_2D( 1, 0, 1, 0 )
	515	z2d(ji,jj) = 0.5_wp * z2d(ji,jj) * z2d(ji,jj)
	516	END_2D
[14286]	517	CALL lbc_lnk( 'diawri', z2d, 'F', 1._wp )
[14200]	518	CALL iom_put( "ssEns", z2d ) ! potential enstrophy ( 1/2*q2 )
[14053]	519	!
	520	ENDIF
[6140]	521
[9124]	522	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	523	!
[1482]	524	END SUBROUTINE dia_wri
	525
	526	#else
[2528]	527	!!----------------------------------------------------------------------
	528	!! Default option use IOIPSL library
	529	!!----------------------------------------------------------------------
	530
[9652]	531	INTEGER FUNCTION dia_wri_alloc()
	532	!!----------------------------------------------------------------------
	533	INTEGER, DIMENSION(2) :: ierr
	534	!!----------------------------------------------------------------------
[12493]	535	IF( nn_write == -1 ) THEN
	536	dia_wri_alloc = 0
	537	ELSE
	538	ierr = 0
	539	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
	540	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
	541	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
[9652]	542	!
[12493]	543	dia_wri_alloc = MAXVAL(ierr)
	544	CALL mpp_sum( 'diawri', dia_wri_alloc )
	545	!
	546	ENDIF
[9652]	547	!
	548	END FUNCTION dia_wri_alloc
[12377]	549
	550	INTEGER FUNCTION dia_wri_alloc_abl()
	551	!!----------------------------------------------------------------------
	552	ALLOCATE( ndex_hA(jpijpj), ndex_A (jpijpj*jpkam1), STAT=dia_wri_alloc_abl)
	553	CALL mpp_sum( 'diawri', dia_wri_alloc_abl )
	554	!
	555	END FUNCTION dia_wri_alloc_abl
[9652]	556
	557
[12377]	558	SUBROUTINE dia_wri( kt, Kmm )
[3]	559	!!---------------------------------------------------------------------
	560	!! * ROUTINE dia_wri *
	561	!!
	562	!! ** Purpose : Standard output of opa: dynamics and tracer fields
	563	!! NETCDF format is used by default
	564	!!
	565	!! ** Method : At the beginning of the first time step (nit000),
	566	!! define all the NETCDF files and fields
	567	!! At each time step call histdef to compute the mean if ncessary
[11536]	568	!! Each nn_write time step, output the instantaneous or mean fields
[3]	569	!!----------------------------------------------------------------------
[5836]	570	INTEGER, INTENT( in ) :: kt ! ocean time-step index
[12377]	571	INTEGER, INTENT( in ) :: Kmm ! ocean time level index
[5836]	572	!
[2528]	573	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
	574	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
	575	INTEGER :: inum = 11 ! temporary logical unit
[3294]	576	INTEGER :: ji, jj, jk ! dummy loop indices
	577	INTEGER :: ierr ! error code return from allocation
[2528]	578	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
[12377]	579	INTEGER :: ipka ! ABL
[3609]	580	INTEGER :: jn, ierror ! local integers
[6140]	581	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
[5836]	582	!
[9019]	583	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
[13237]	584	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d, ze3t, zgdept ! 3D workspace
[12377]	585	REAL(wp), DIMENSION(:,:,:), ALLOCATABLE :: zw3d_abl ! ABL 3D workspace
[3]	586	!!----------------------------------------------------------------------
[1482]	587	!
[9019]	588	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
[12377]	589	CALL dia_wri_state( Kmm, 'output.init' )
[1482]	590	ninist = 0
	591	ENDIF
	592	!
[11536]	593	IF( nn_write == -1 ) RETURN ! we will never do any output
	594	!
	595	IF( ln_timing ) CALL timing_start('dia_wri')
	596	!
[3]	597	! 0. Initialisation
	598	! -----------------
[632]	599
[9019]	600	ll_print = .FALSE. ! local variable for debugging
[3]	601	ll_print = ll_print .AND. lwp
	602
	603	! Define frequency of output and means
[5566]	604	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
[3]	605	#if defined key_diainstant
[12489]	606	zsto = nn_write * rn_Dt
[1312]	607	clop = "inst("//TRIM(clop)//")"
[3]	608	#else
[12489]	609	zsto=rn_Dt
[1312]	610	clop = "ave("//TRIM(clop)//")"
[3]	611	#endif
[12489]	612	zout = nn_write * rn_Dt
	613	zmax = ( nitend - nit000 + 1 ) * rn_Dt
[3]	614
	615	! Define indices of the horizontal output zoom and vertical limit storage
[13286]	616	iimi = Nis0 ; iima = Nie0
	617	ijmi = Njs0 ; ijma = Nje0
[3]	618	ipk = jpk
[12377]	619	IF(ln_abl) ipka = jpkam1
[3]	620
	621	! define time axis
[1334]	622	it = kt
	623	itmod = kt - nit000 + 1
[3]	624
[13237]	625	! store e3t for subsitute
	626	DO jk = 1, jpk
	627	ze3t (:,:,jk) = e3t (:,:,jk,Kmm)
	628	zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
	629	END DO
[3]	630
[13237]	631
[3]	632	! 1. Define NETCDF files and fields at beginning of first time step
	633	! -----------------------------------------------------------------
	634
	635	IF( kt == nit000 ) THEN
	636
	637	! Define the NETCDF files (one per grid)
[632]	638
[3]	639	! Compute julian date from starting date of the run
[12489]	640	CALL ymds2ju( nyear, nmonth, nday, rn_Dt, zjulian )
[1309]	641	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
[3]	642	IF(lwp)WRITE(numout,*)
	643	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
	644	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
	645	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
	646	' limit storage in depth = ', ipk
	647
	648	! WRITE root name in date.file for use by postpro
[1581]	649	IF(lwp) THEN
[11536]	650	CALL dia_nam( clhstnam, nn_write,' ' )
[1581]	651	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
[895]	652	WRITE(inum,*) clhstnam
	653	CLOSE(inum)
	654	ENDIF
[632]	655
[3]	656	! Define the T grid FILE ( nid_T )
[632]	657
[11536]	658	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
[3]	659	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	660	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	661	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	662	& nit000-1, zjulian, rn_Dt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
[3]	663	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
[4292]	664	& "m", ipk, gdept_1d, nz_T, "down" )
[3]	665	! ! Index of ocean points
	666	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
	667	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
[3609]	668	!
	669	IF( ln_icebergs ) THEN
	670	!
	671	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
	672	!! that routine is called from nemogcm, so do it here immediately before its needed
	673	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
[10425]	674	CALL mpp_sum( 'diawri', ierror )
[3609]	675	IF( ierror /= 0 ) THEN
	676	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
	677	RETURN
	678	ENDIF
	679	!
	680	!! iceberg vertical coordinate is class number
	681	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
	682	& "number", nclasses, class_num, nb_T )
	683	!
	684	!! each class just needs the surface index pattern
	685	ndim_bT = 3
	686	DO jn = 1,nclasses
	687	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
	688	ENDDO
	689	!
	690	ENDIF
[3]	691
	692	! Define the U grid FILE ( nid_U )
	693
[11536]	694	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
[3]	695	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	696	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
	697	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	698	& nit000-1, zjulian, rn_Dt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
[3]	699	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
[4292]	700	& "m", ipk, gdept_1d, nz_U, "down" )
[3]	701	! ! Index of ocean points
	702	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
	703	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
	704
	705	! Define the V grid FILE ( nid_V )
	706
[11536]	707	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
[3]	708	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	709	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
	710	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	711	& nit000-1, zjulian, rn_Dt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
[3]	712	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
[4292]	713	& "m", ipk, gdept_1d, nz_V, "down" )
[3]	714	! ! Index of ocean points
	715	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
	716	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
	717
	718	! Define the W grid FILE ( nid_W )
	719
[11536]	720	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
[3]	721	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
	722	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	723	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	724	& nit000-1, zjulian, rn_Dt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
[3]	725	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
[4292]	726	& "m", ipk, gdepw_1d, nz_W, "down" )
[3]	727
[12377]	728	IF( ln_abl ) THEN
	729	! Define the ABL grid FILE ( nid_A )
	730	CALL dia_nam( clhstnam, nn_write, 'grid_ABL' )
	731	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
	732	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
	733	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
[12489]	734	& nit000-1, zjulian, rn_Dt, nh_A, nid_A, domain_id=nidom, snc4chunks=snc4set )
[12377]	735	CALL histvert( nid_A, "ght_abl", "Vertical T levels", & ! Vertical grid: gdept
	736	& "m", ipka, ght_abl(2:jpka), nz_A, "up" )
	737	! ! Index of ocean points
	738	ALLOCATE( zw3d_abl(jpi,jpj,ipka) )
	739	zw3d_abl(:,:,:) = 1._wp
	740	CALL wheneq( jpijpjipka, zw3d_abl, 1, 1., ndex_A , ndim_A ) ! volume
	741	CALL wheneq( jpi*jpj , zw3d_abl, 1, 1., ndex_hA, ndim_hA ) ! surface
	742	DEALLOCATE(zw3d_abl)
	743	ENDIF
[13237]	744	!
[632]	745
[3]	746	! Declare all the output fields as NETCDF variables
	747
	748	! !!! nid_T : 3D
	749	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
	750	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	751	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
	752	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[6140]	753	IF( .NOT.ln_linssh ) THEN
[13237]	754	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t n
[4292]	755	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13237]	756	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t n
[4292]	757	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
[13237]	758	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t n
[4292]	759	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
	760	ENDIF
[3]	761	! !!! nid_T : 2D
	762	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
	763	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	764	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
	765	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[359]	766	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
[3]	767	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[2528]	768	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
[3]	769	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[4570]	770	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
	771	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3625]	772	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
[3]	773	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[6140]	774	IF( ln_linssh ) THEN
[12377]	775	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * ts(:,:,1,jp_tem,Kmm)
[3625]	776	& , "KgC/m2/s", & ! sosst_cd
[4292]	777	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	778	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * ts(:,:,1,jp_sal,Kmm)
[3625]	779	& , "KgPSU/m2/s",& ! sosss_cd
[4292]	780	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	781	ENDIF
[888]	782	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
[3]	783	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	784	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
	785	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1585]	786	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
	787	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[3]	788	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
	789	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1037]	790	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
[3]	791	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[1649]	792	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
	793	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[12377]	794	!
	795	IF( ln_abl ) THEN
	796	CALL histdef( nid_A, "t_abl", "Potential Temperature" , "K" , & ! t_abl
	797	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	798	CALL histdef( nid_A, "q_abl", "Humidity" , "kg/kg" , & ! q_abl
	799	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	800	CALL histdef( nid_A, "u_abl", "Atmospheric U-wind " , "m/s" , & ! u_abl
	801	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	802	CALL histdef( nid_A, "v_abl", "Atmospheric V-wind " , "m/s" , & ! v_abl
	803	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	804	CALL histdef( nid_A, "tke_abl", "Atmospheric TKE " , "m2/s2" , & ! tke_abl
	805	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	806	CALL histdef( nid_A, "avm_abl", "Atmospheric turbulent viscosity", "m2/s" , & ! avm_abl
	807	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	808	CALL histdef( nid_A, "avt_abl", "Atmospheric turbulent diffusivity", "m2/s2", & ! avt_abl
	809	& jpi, jpj, nh_A, ipka, 1, ipka, nz_A, 32, clop, zsto, zout )
	810	CALL histdef( nid_A, "pblh", "Atmospheric boundary layer height " , "m", & ! pblh
	811	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	812	#if defined key_si3
	813	CALL histdef( nid_A, "oce_frac", "Fraction of open ocean" , " ", & ! ato_i
	814	& jpi, jpj, nh_A, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	815	#endif
	816	CALL histend( nid_A, snc4chunks=snc4set )
	817	ENDIF
	818	!
[3609]	819	IF( ln_icebergs ) THEN
	820	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
	821	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	822	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
	823	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	824	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
	825	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	826	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
	827	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	828	IF( ln_bergdia ) THEN
	829	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
	830	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	831	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
	832	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	833	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
	834	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	835	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
	836	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	837	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
	838	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	839	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
	840	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	841	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
	842	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	843	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
	844	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	845	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
	846	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	847	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
	848	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
	849	ENDIF
	850	ENDIF
	851
[11536]	852	IF( ln_ssr ) THEN
[4990]	853	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
	854	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	855	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
	856	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	857	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
	858	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	859	ENDIF
[11536]	860
[3]	861	clmx ="l_max(only(x))" ! max index on a period
[5836]	862	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
	863	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
[3]	864	#if defined key_diahth
	865	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
	866	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	867	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
	868	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	869	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
	870	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
[7646]	871	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
[3]	872	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
	873	#endif
	874
[2528]	875	CALL histend( nid_T, snc4chunks=snc4set )
[3]	876
	877	! !!! nid_U : 3D
[12377]	878	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! uu(:,:,:,Kmm)
[3]	879	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
[7646]	880	IF( ln_wave .AND. ln_sdw) THEN
	881	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
	882	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
	883	ENDIF
[3]	884	! !!! nid_U : 2D
[888]	885	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
[3]	886	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	887
[2528]	888	CALL histend( nid_U, snc4chunks=snc4set )
[3]	889
	890	! !!! nid_V : 3D
[12377]	891	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vv(:,:,:,Kmm)
[3]	892	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
[7646]	893	IF( ln_wave .AND. ln_sdw) THEN
	894	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
	895	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
	896	ENDIF
[3]	897	! !!! nid_V : 2D
[888]	898	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
[3]	899	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
	900
[2528]	901	CALL histend( nid_V, snc4chunks=snc4set )
[3]	902
	903	! !!! nid_W : 3D
[12377]	904	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! ww
[3]	905	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	906	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
	907	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
[9019]	908	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
[255]	909	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	910
[9019]	911	IF( ln_zdfddm ) THEN
[3]	912	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
	913	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	914	ENDIF
[7646]	915
	916	IF( ln_wave .AND. ln_sdw) THEN
	917	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
	918	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
	919	ENDIF
[3]	920	! !!! nid_W : 2D
[2528]	921	CALL histend( nid_W, snc4chunks=snc4set )
[3]	922
	923	IF(lwp) WRITE(numout,*)
	924	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
	925	IF(ll_print) CALL FLUSH(numout )
	926
	927	ENDIF
	928
	929	! 2. Start writing data
	930	! ---------------------
	931
[4292]	932	! ndex(1) est utilise ssi l'avant dernier argument est different de
[3]	933	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
	934	! donne le nombre d'elements, et ndex la liste des indices a sortir
	935
[11536]	936	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
[3]	937	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
	938	WRITE(numout,*) '~~~~~~ '
	939	ENDIF
	940
[6140]	941	IF( .NOT.ln_linssh ) THEN
[13237]	942	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! heat content
	943	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) * ze3t(:,:,:) , ndim_T , ndex_T ) ! salt content
	944	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
	945	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) * ze3t(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
[4292]	946	ELSE
[12377]	947	CALL histwrite( nid_T, "votemper", it, ts(:,:,:,jp_tem,Kmm) , ndim_T , ndex_T ) ! temperature
	948	CALL histwrite( nid_T, "vosaline", it, ts(:,:,:,jp_sal,Kmm) , ndim_T , ndex_T ) ! salinity
	949	CALL histwrite( nid_T, "sosstsst", it, ts(:,:,1,jp_tem,Kmm) , ndim_hT, ndex_hT ) ! sea surface temperature
	950	CALL histwrite( nid_T, "sosaline", it, ts(:,:,1,jp_sal,Kmm) , ndim_hT, ndex_hT ) ! sea surface salinity
[4292]	951	ENDIF
[6140]	952	IF( .NOT.ln_linssh ) THEN
[13237]	953	zw3d(:,:,:) = ( ( ze3t(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
	954	CALL histwrite( nid_T, "vovvle3t", it, ze3t (:,:,:) , ndim_T , ndex_T ) ! level thickness
	955	CALL histwrite( nid_T, "vovvldep", it, zgdept , ndim_T , ndex_T ) ! t-point depth
[4292]	956	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
	957	ENDIF
[12377]	958	CALL histwrite( nid_T, "sossheig", it, ssh(:,:,Kmm) , ndim_hT, ndex_hT ) ! sea surface height
[2528]	959	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
[4570]	960	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
[3625]	961	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
	962	! (includes virtual salt flux beneath ice
	963	! in linear free surface case)
[6140]	964	IF( ln_linssh ) THEN
[12377]	965	zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_tem,Kmm)
[4292]	966	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
[12377]	967	zw2d(:,:) = emp (:,:) * ts(:,:,1,jp_sal,Kmm)
[4292]	968	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
	969	ENDIF
[888]	970	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
[3]	971	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
[1585]	972	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
[3]	973	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
[1037]	974	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
[1649]	975	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
[12377]	976	!
	977	IF( ln_abl ) THEN
	978	ALLOCATE( zw3d_abl(jpi,jpj,jpka) )
	979	IF( ln_mskland ) THEN
	980	DO jk=1,jpka
	981	zw3d_abl(:,:,jk) = tmask(:,:,1)
	982	END DO
	983	ELSE
	984	zw3d_abl(:,:,:) = 1._wp
	985	ENDIF
	986	CALL histwrite( nid_A, "pblh" , it, pblh(:,:) *zw3d_abl(:,:,1 ), ndim_hA, ndex_hA ) ! pblh
	987	CALL histwrite( nid_A, "u_abl" , it, u_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! u_abl
	988	CALL histwrite( nid_A, "v_abl" , it, v_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! v_abl
	989	CALL histwrite( nid_A, "t_abl" , it, tq_abl (:,:,2:jpka,nt_n,1)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! t_abl
	990	CALL histwrite( nid_A, "q_abl" , it, tq_abl (:,:,2:jpka,nt_n,2)*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! q_abl
	991	CALL histwrite( nid_A, "tke_abl", it, tke_abl (:,:,2:jpka,nt_n )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! tke_abl
	992	CALL histwrite( nid_A, "avm_abl", it, avm_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avm_abl
	993	CALL histwrite( nid_A, "avt_abl", it, avt_abl (:,:,2:jpka )*zw3d_abl(:,:,2:jpka), ndim_A , ndex_A ) ! avt_abl
	994	#if defined key_si3
	995	CALL histwrite( nid_A, "oce_frac" , it, ato_i(:,:) , ndim_hA, ndex_hA ) ! ato_i
	996	#endif
	997	DEALLOCATE(zw3d_abl)
	998	ENDIF
	999	!
[3609]	1000	IF( ln_icebergs ) THEN
	1001	!
	1002	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
	1003	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
	1004	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
	1005	!
	1006	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
	1007	!
	1008	IF( ln_bergdia ) THEN
	1009	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
	1010	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
	1011	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
	1012	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
	1013	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
	1014	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
	1015	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
	1016	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
	1017	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
	1018	!
	1019	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
	1020	ENDIF
	1021	ENDIF
	1022
[11536]	1023	IF( ln_ssr ) THEN
[4990]	1024	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
	1025	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
[12377]	1026	zw2d(:,:) = erp(:,:) * ts(:,:,1,jp_sal,Kmm) * tmask(:,:,1)
[4990]	1027	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
	1028	ENDIF
	1029	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
	1030	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
[3]	1031
	1032	#if defined key_diahth
	1033	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
	1034	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
	1035	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
	1036	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
	1037	#endif
[888]	1038
[12377]	1039	CALL histwrite( nid_U, "vozocrtx", it, uu(:,:,:,Kmm) , ndim_U , ndex_U ) ! i-current
[888]	1040	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
[3]	1041
[12377]	1042	CALL histwrite( nid_V, "vomecrty", it, vv(:,:,:,Kmm) , ndim_V , ndex_V ) ! j-current
[888]	1043	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
[3]	1044
[11418]	1045	IF( ln_zad_Aimp ) THEN
[12377]	1046	CALL histwrite( nid_W, "vovecrtz", it, ww + wi , ndim_T, ndex_T ) ! vert. current
[11418]	1047	ELSE
[12377]	1048	CALL histwrite( nid_W, "vovecrtz", it, ww , ndim_T, ndex_T ) ! vert. current
[11418]	1049	ENDIF
[3]	1050	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
[9019]	1051	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
	1052	IF( ln_zdfddm ) THEN
	1053	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
[3]	1054	ENDIF
	1055
[7646]	1056	IF( ln_wave .AND. ln_sdw ) THEN
[9019]	1057	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
	1058	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
	1059	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
[7646]	1060	ENDIF
	1061
[1318]	1062	! 3. Close all files
	1063	! ---------------------------------------
[1561]	1064	IF( kt == nitend ) THEN
[3]	1065	CALL histclo( nid_T )
	1066	CALL histclo( nid_U )
	1067	CALL histclo( nid_V )
	1068	CALL histclo( nid_W )
[12377]	1069	IF(ln_abl) CALL histclo( nid_A )
[3]	1070	ENDIF
[2528]	1071	!
[9124]	1072	IF( ln_timing ) CALL timing_stop('dia_wri')
[3294]	1073	!
[3]	1074	END SUBROUTINE dia_wri
[1567]	1075	#endif
	1076
[12377]	1077	SUBROUTINE dia_wri_state( Kmm, cdfile_name )
[3]	1078	!!---------------------------------------------------------------------
	1079	!! * ROUTINE dia_wri_state *
	1080	!!
	1081	!! ** Purpose : create a NetCDF file named cdfile_name which contains
	1082	!! the instantaneous ocean state and forcing fields.
	1083	!! Used to find errors in the initial state or save the last
	1084	!! ocean state in case of abnormal end of a simulation
	1085	!!
	1086	!! ** Method : NetCDF files using ioipsl
	1087	!! File 'output.init.nc' is created if ninist = 1 (namelist)
	1088	!! File 'output.abort.nc' is created in case of abnormal job end
	1089	!!----------------------------------------------------------------------
[12377]	1090	INTEGER , INTENT( in ) :: Kmm ! time level index
[1334]	1091	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
[10425]	1092	!!
[12377]	1093	INTEGER :: inum, jk
[15540]	1094	REAL(dp), DIMENSION(jpi,jpj,jpk) :: ze3t, zgdept ! 3D workspace for qco substitution
[3]	1095	!!----------------------------------------------------------------------
[3294]	1096	!
[13237]	1097	IF(lwp) THEN
	1098	WRITE(numout,*)
	1099	WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
	1100	WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
	1101	WRITE(numout,*) ' and named :', cdfile_name, '...nc'
	1102	ENDIF
[12649]	1103	!
[13237]	1104	DO jk = 1, jpk
	1105	ze3t(:,:,jk) = e3t(:,:,jk,Kmm)
	1106	zgdept(:,:,jk) = gdept(:,:,jk,Kmm)
	1107	END DO
	1108	!
[12649]	1109	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
	1110	!
[12377]	1111	CALL iom_rstput( 0, 0, inum, 'votemper', ts(:,:,:,jp_tem,Kmm) ) ! now temperature
	1112	CALL iom_rstput( 0, 0, inum, 'vosaline', ts(:,:,:,jp_sal,Kmm) ) ! now salinity
	1113	CALL iom_rstput( 0, 0, inum, 'sossheig', ssh(:,:,Kmm) ) ! sea surface height
	1114	CALL iom_rstput( 0, 0, inum, 'vozocrtx', uu(:,:,:,Kmm) ) ! now i-velocity
	1115	CALL iom_rstput( 0, 0, inum, 'vomecrty', vv(:,:,:,Kmm) ) ! now j-velocity
[11418]	1116	IF( ln_zad_Aimp ) THEN
[12377]	1117	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww + wi ) ! now k-velocity
[11418]	1118	ELSE
[12377]	1119	CALL iom_rstput( 0, 0, inum, 'vovecrtz', ww ) ! now k-velocity
[11418]	1120	ENDIF
[12377]	1121	CALL iom_rstput( 0, 0, inum, 'risfdep', risfdep ) ! now k-velocity
[13237]	1122	CALL iom_rstput( 0, 0, inum, 'ht' , ht(:,:) ) ! now water column height
[12649]	1123	!
[12377]	1124	IF ( ln_isf ) THEN
	1125	IF (ln_isfcav_mlt) THEN
	1126	CALL iom_rstput( 0, 0, inum, 'fwfisf_cav', fwfisf_cav ) ! now k-velocity
	1127	CALL iom_rstput( 0, 0, inum, 'rhisf_cav_tbl', rhisf_tbl_cav ) ! now k-velocity
	1128	CALL iom_rstput( 0, 0, inum, 'rfrac_cav_tbl', rfrac_tbl_cav ) ! now k-velocity
[14219]	1129	CALL iom_rstput( 0, 0, inum, 'misfkb_cav', REAL(misfkb_cav,dp) ) ! now k-velocity
	1130	CALL iom_rstput( 0, 0, inum, 'misfkt_cav', REAL(misfkt_cav,dp) ) ! now k-velocity
	1131	CALL iom_rstput( 0, 0, inum, 'mskisf_cav', REAL(mskisf_cav,dp), ktype = jp_i1 )
[12377]	1132	END IF
	1133	IF (ln_isfpar_mlt) THEN
[14219]	1134	CALL iom_rstput( 0, 0, inum, 'isfmsk_par', REAL(mskisf_par,dp) ) ! now k-velocity
[12377]	1135	CALL iom_rstput( 0, 0, inum, 'fwfisf_par', fwfisf_par ) ! now k-velocity
	1136	CALL iom_rstput( 0, 0, inum, 'rhisf_par_tbl', rhisf_tbl_par ) ! now k-velocity
	1137	CALL iom_rstput( 0, 0, inum, 'rfrac_par_tbl', rfrac_tbl_par ) ! now k-velocity
[14219]	1138	CALL iom_rstput( 0, 0, inum, 'misfkb_par', REAL(misfkb_par,dp) ) ! now k-velocity
	1139	CALL iom_rstput( 0, 0, inum, 'misfkt_par', REAL(misfkt_par,dp) ) ! now k-velocity
	1140	CALL iom_rstput( 0, 0, inum, 'mskisf_par', REAL(mskisf_par,dp), ktype = jp_i1 )
[12377]	1141	END IF
	1142	END IF
[12649]	1143	!
[7646]	1144	IF( ALLOCATED(ahtu) ) THEN
[10425]	1145	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
	1146	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
[7646]	1147	ENDIF
	1148	IF( ALLOCATED(ahmt) ) THEN
[10425]	1149	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
	1150	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
[7646]	1151	ENDIF
[10425]	1152	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
	1153	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
	1154	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
	1155	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
	1156	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
	1157	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
[6140]	1158	IF( .NOT.ln_linssh ) THEN
[13237]	1159	CALL iom_rstput( 0, 0, inum, 'vovvldep', zgdept ) ! T-cell depth
	1160	CALL iom_rstput( 0, 0, inum, 'vovvle3t', ze3t ) ! T-cell thickness
[10425]	1161	END IF
[7646]	1162	IF( ln_wave .AND. ln_sdw ) THEN
[10425]	1163	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
	1164	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
	1165	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
[7646]	1166	ENDIF
[12377]	1167	IF ( ln_abl ) THEN
	1168	CALL iom_rstput ( 0, 0, inum, "uz1_abl", u_abl(:,:,2,nt_a ) ) ! now first level i-wind
	1169	CALL iom_rstput ( 0, 0, inum, "vz1_abl", v_abl(:,:,2,nt_a ) ) ! now first level j-wind
	1170	CALL iom_rstput ( 0, 0, inum, "tz1_abl", tq_abl(:,:,2,nt_a,1) ) ! now first level temperature
	1171	CALL iom_rstput ( 0, 0, inum, "qz1_abl", tq_abl(:,:,2,nt_a,2) ) ! now first level humidity
	1172	ENDIF
[14045]	1173	IF( ln_zdfosm ) THEN
	1174	CALL iom_rstput( 0, 0, inum, 'hbl', hbl*tmask(:,:,1) ) ! now boundary-layer depth
	1175	CALL iom_rstput( 0, 0, inum, 'hml', hml*tmask(:,:,1) ) ! now mixed-layer depth
	1176	CALL iom_rstput( 0, 0, inum, 'avt_k', avt_k*wmask ) ! w-level diffusion
	1177	CALL iom_rstput( 0, 0, inum, 'avm_k', avm_k*wmask ) ! now w-level viscosity
	1178	CALL iom_rstput( 0, 0, inum, 'ghamt', ghamt*wmask ) ! non-local t forcing
	1179	CALL iom_rstput( 0, 0, inum, 'ghams', ghams*wmask ) ! non-local s forcing
	1180	CALL iom_rstput( 0, 0, inum, 'ghamu', ghamu*umask ) ! non-local u forcing
	1181	CALL iom_rstput( 0, 0, inum, 'ghamv', ghamv*vmask ) ! non-local v forcing
	1182	IF( ln_osm_mle ) THEN
	1183	CALL iom_rstput( 0, 0, inum, 'hmle', hmle*tmask(:,:,1) ) ! now transition-layer depth
	1184	END IF
	1185	ENDIF
[12649]	1186	!
	1187	CALL iom_close( inum )
	1188	!
[10425]	1189	#if defined key_si3
	1190	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
[12649]	1191	CALL iom_open( TRIM(cdfile_name)//'_ice', inum, ldwrt = .TRUE., kdlev = jpl, cdcomp = 'ICE' )
[10425]	1192	CALL ice_wri_state( inum )
[12649]	1193	CALL iom_close( inum )
[1561]	1194	ENDIF
[12933]	1195	!
[1561]	1196	#endif
[3]	1197	END SUBROUTINE dia_wri_state
[6140]	1198
[3]	1199	!!======================================================================
	1200	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/branches/2021/dev_r14116_HPC-10_mcastril_Mixed_Precision_implementation/src/OCE/DIA/diawri.F90

Download in other formats: