Context Navigation

← Previous Revision
Next Revision →
Blame
Revision Log

diawri.F90

Last change on this file was 12555, checked in by charris, 4 years ago

Changes from GO6 package branch (GMED ticket 450):

svn merge -r 11035:11101 svn+ssh://charris@forge.ipsl.jussieu.fr/ipsl/forge/projets/nemo/svn/branches/UKMO/dev_r5518_GO6_package

File size: 58.9 KB

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

Download in other formats: