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 branches/UKMO/dev_r5518_fix_diag_bitcomp/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_fix_diag_bitcomp/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 9502

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

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

Download in other formats: