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 @ 13020

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

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

Download in other formats: