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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7651

Last change on this file since 7651 was 7651, checked in by timgraham, 7 years ago
Added iom_put call for wpt_dep (used for CMIP6 zhalf) as this came into the data request quite late
File size: 58.7 KB

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

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

Download in other formats: