Context Navigation

diawri.F90 @ 7033

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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_MEDUSA_extra_CMIP6_diags/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7033

Download in other formats: