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

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

diawri.F90 in NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2020/dev_r14116_HPC-04_mcastril_Mixed_Precision_implementation_final/src/OCE/DIA/diawri.F90 @ 14286

Last change on this file since 14286 was 14286, checked in by mcastril, 3 years ago
Reformatting and allowing to use key_qco
Property svn:keywords set to `Id`
File size: 64.5 KB

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

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

Download in other formats: