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/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/CANAL/MY_SRC – NEMO

Context Navigation

source: branches/2017/dev_merge_2017/NEMOGCM/CONFIG/TEST_CASES/CANAL/MY_SRC/diawri.F90 @ 9403

Last change on this file since 9403 was 9403, checked in by smasson, 6 years ago
dev_merge_2017: update test_cases CANAL
File size: 60.0 KB

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

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

Download in other formats: