Context Navigation

diawri.F90 @ 12918

Last change on this file since 12918 was 12918, checked in by petesykes, 4 years ago
Land masking for operational diagnostics
File size: 60.0 KB

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

Download in other formats: