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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_text_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 10774

Last change on this file since 10774 was 10774, checked in by andmirek, 5 years ago
GMED 450 add flush after prints
File size: 59.0 KB

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

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

Download in other formats: