Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diawri.F90 @ 10207

Last change on this file since 10207 was 10207, checked in by cmao, 6 years ago
remove svn keyword
File size: 58.3 KB

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

Download in other formats: