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_GO6_package_FOAMv14_ssh_dt/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_FOAMv14_ssh_dt/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 13074

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

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

Download in other formats: