Context Navigation

diawri.F90 @ 6448

Last change on this file since 6448 was 6448, checked in by dancopsey, 8 years ago
Merged in changeset 5239 of dev_r5021_nn_etau_revision
File size: 57.7 KB

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

Download in other formats: