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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_for_static_diagnostics/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7583

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

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

Download in other formats: