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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diawri.F90 in branches/UKMO/CO6_KD490/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/CO6_KD490/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6332

Last change on this file since 6332 was 6332, checked in by deazer, 8 years ago
Tested Initial run one day physics only in rose suite.
File size: 56.8 KB

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

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

Download in other formats: