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/dev_r5518_GO6_package_XIOS_read/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_XIOS_read/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7923

Last change on this file since 7923 was 7923, checked in by andmirek, 7 years ago
merge changes up to 7573
File size: 58.6 KB

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

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

Download in other formats: