Context Navigation

diawri.F90 @ 15343

Last change on this file since 15343 was 15333, checked in by hadjt, 3 years ago

Adding Regional Means, but without XIOS or MLD.

Search on
!JT MLD
!JT IOM

in IOM/iom.F90 and DIA/diaregmean.F90

to see the code commented out.

File size: 51.1 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 dia25h ! 25h Mean output
49	!JT
50	USE diaregmean ! regionalmean
51	USE diapea ! pea
52	!JT
53	USE iom !
54	USE ioipsl !
55	USE eosbn2
56	#if defined key_si3
57	USE ice
58	USE icewri
59	#endif
60	USE lib_mpp ! MPP library
61	USE timing ! preformance summary
62	USE diurnal_bulk ! diurnal warm layer
63	USE cool_skin ! Cool skin
64
65	IMPLICIT NONE
66	PRIVATE
67
68	PUBLIC dia_wri ! routines called by step.F90
69	PUBLIC dia_wri_state
70	PUBLIC dia_wri_alloc ! Called by nemogcm module
71
72	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
73	INTEGER :: nb_T , ndim_bT ! grid_T file
74	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
75	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
76	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
77	INTEGER :: ndex(1) ! ???
78	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
79	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
80	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
81
82	!! * Substitutions
83	# include "vectopt_loop_substitute.h90"
84	!!----------------------------------------------------------------------
85	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
86	!! $Id$
87	!! Software governed by the CeCILL license (see ./LICENSE)
88	!!----------------------------------------------------------------------
89	CONTAINS
90
91	#if defined key_iomput
92	!!----------------------------------------------------------------------
93	!! 'key_iomput' use IOM library
94	!!----------------------------------------------------------------------
95	INTEGER FUNCTION dia_wri_alloc()
96	!
97	dia_wri_alloc = 0
98	!
99	END FUNCTION dia_wri_alloc
100
101
102	SUBROUTINE dia_wri( kt )
103	!!---------------------------------------------------------------------
104	!! * ROUTINE dia_wri *
105	!!
106	!! ** Purpose : Standard output of opa: dynamics and tracer fields
107	!! NETCDF format is used by default
108	!!
109	!! ** Method : use iom_put
110	!!----------------------------------------------------------------------
111	INTEGER, INTENT( in ) :: kt ! ocean time-step index
112	!!
113	INTEGER :: ji, jj, jk ! dummy loop indices
114	INTEGER :: ikbot ! local integer
115	REAL(wp):: zztmp , zztmpx ! local scalar
116	REAL(wp):: zztmp2, zztmpy ! - -
117	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
118	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
119	CHARACTER(len=4),SAVE :: ttype , stype ! temperature and salinity type
120	!!----------------------------------------------------------------------
121	!
122	IF( kt == nit000 ) THEN
123	IF( ln_TEOS10 ) THEN
124	IF ( iom_use("toce_pot") .OR. iom_use("soce_pra") .OR. iom_use("sst_pot") .OR. iom_use("sss_pra") &
125	& .OR. iom_use("sbt_pot") .OR. iom_use("sbs_pra") .OR. iom_use("sstgrad_pot") .OR. iom_use("sstgrad2_pot") &
126	& .OR. iom_use("tosmint_pot") .OR. iom_use("somint_pra")) THEN
127	CALL ctl_stop( 'diawri: potential temperature and practical salinity not available with ln_TEOS10' )
128	ELSE
129	ttype='con' ; stype='abs' ! teos-10 using conservative temperature and absolute salinity
130	ENDIF
131	ELSE IF( ln_EOS80 ) THEN
132	IF ( iom_use("toce_con") .OR. iom_use("soce_abs") .OR. iom_use("sst_con") .OR. iom_use("sss_abs") &
133	& .OR. iom_use("sbt_con") .OR. iom_use("sbs_abs") .OR. iom_use("sstgrad_con") .OR. iom_use("sstgrad2_con") &
134	& .OR. iom_use("tosmint_con") .OR. iom_use("somint_abs")) THEN
135	CALL ctl_stop( 'diawri: conservative temperature and absolute salinity not available with ln_EOS80' )
136	ELSE
137	ttype='pot' ; stype='pra' ! eos-80 using potential temperature and practical salinity
138	ENDIF
139	ELSE IF ( ln_SEOS) THEN
140	ttype='seos' ; stype='seos' ! seos using Simplified Equation of state
141	ENDIF
142	ENDIF
143
144	IF( ln_timing ) CALL timing_start('dia_wri')
145	!
146	! Output the initial state and forcings
147	IF( ninist == 1 ) THEN
148	CALL dia_wri_state( 'output.init' )
149	ninist = 0
150	ENDIF
151
152	! Output of initial vertical scale factor
153	CALL iom_put("e3t_0", e3t_0(:,:,:) )
154	CALL iom_put("e3u_0", e3u_0(:,:,:) )
155	CALL iom_put("e3v_0", e3v_0(:,:,:) )
156	!
157	CALL iom_put( "e3t" , e3t_n(:,:,:) )
158	CALL iom_put( "e3u" , e3u_n(:,:,:) )
159	CALL iom_put( "e3v" , e3v_n(:,:,:) )
160	CALL iom_put( "e3w" , e3w_n(:,:,:) )
161	IF( iom_use("e3tdef") ) &
162	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
163
164	IF( ll_wd ) THEN
165	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
166	ELSE
167	CALL iom_put( "ssh" , sshn ) ! sea surface height
168	ENDIF
169
170	IF( iom_use("wetdep") ) & ! wet depth
171	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
172
173	CALL iom_put( "toce_"//ttype, tsn(:,:,:,jp_tem) ) ! 3D temperature
174	CALL iom_put( "sst_"//ttype, tsn(:,:,1,jp_tem) ) ! surface temperature
175	IF ( iom_use("sbt_"//ttype) ) THEN
176	DO jj = 1, jpj
177	DO ji = 1, jpi
178	ikbot = mbkt(ji,jj)
179	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
180	END DO
181	END DO
182	CALL iom_put( "sbt_"//ttype, z2d ) ! bottom temperature
183	ENDIF
184
185	CALL iom_put( "soce_"//stype, tsn(:,:,:,jp_sal) ) ! 3D salinity
186	CALL iom_put( "sss_"//stype, tsn(:,:,1,jp_sal) ) ! surface salinity
187	IF ( iom_use("sbs_"//stype) ) THEN
188	DO jj = 1, jpj
189	DO ji = 1, jpi
190	ikbot = mbkt(ji,jj)
191	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
192	END DO
193	END DO
194	CALL iom_put( "sbs_"//stype, z2d ) ! bottom salinity
195	ENDIF
196
197	CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
198
199	IF ( iom_use("taubot") ) THEN ! bottom stress
200	zztmp = rau0 * 0.25
201	z2d(:,:) = 0._wp
202	DO jj = 2, jpjm1
203	DO ji = fs_2, fs_jpim1 ! vector opt.
204	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 &
205	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 &
206	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 &
207	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2
208	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
209	!
210	END DO
211	END DO
212	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
213	CALL iom_put( "taubot", z2d )
214	ENDIF
215
216	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
217	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
218	IF ( iom_use("sbu") ) THEN
219	DO jj = 1, jpj
220	DO ji = 1, jpi
221	ikbot = mbku(ji,jj)
222	z2d(ji,jj) = un(ji,jj,ikbot)
223	END DO
224	END DO
225	CALL iom_put( "sbu", z2d ) ! bottom i-current
226	ENDIF
227
228	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
229	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
230	IF ( iom_use("sbv") ) THEN
231	DO jj = 1, jpj
232	DO ji = 1, jpi
233	ikbot = mbkv(ji,jj)
234	z2d(ji,jj) = vn(ji,jj,ikbot)
235	END DO
236	END DO
237	CALL iom_put( "sbv", z2d ) ! bottom j-current
238	ENDIF
239
240	IF( ln_zad_Aimp ) wn = wn + wi ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
241	!
242	CALL iom_put( "woce", wn ) ! vertical velocity
243	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
244	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
245	z2d(:,:) = rau0 * e1e2t(:,:)
246	DO jk = 1, jpk
247	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
248	END DO
249	CALL iom_put( "w_masstr" , z3d )
250	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
251	ENDIF
252	!
253	IF( ln_zad_Aimp ) wn = wn - wi ! Remove implicit part of vertical velocity that was added for diagnostic output
254
255	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
256	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
257	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
258
259	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
260	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
261
262	IF ( iom_use("sstgrad_"//ttype) .OR. iom_use("sstgrad2_"//ttype) ) 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( 'diawri', z2d, 'T', 1. )
273	CALL iom_put( "sstgrad2_"//ttype, z2d ) ! square of module of sst gradient
274	z2d(:,:) = SQRT( z2d(:,:) )
275	CALL iom_put( "sstgrad_"//ttype , 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("eken") ) THEN
304	z3d(:,:,jpk) = 0._wp
305	DO jk = 1, jpkm1
306	DO jj = 2, jpjm1
307	DO ji = fs_2, fs_jpim1 ! vector opt.
308	zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
309	z3d(ji,jj,jk) = zztmp * ( un(ji-1,jj,jk)*2 e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) &
310	& + un(ji ,jj,jk)*2 e2u(ji ,jj) * e3u_n(ji ,jj,jk) &
311	& + vn(ji,jj-1,jk)*2 e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) &
312	& + vn(ji,jj ,jk)*2 e1v(ji,jj ) * e3v_n(ji,jj ,jk) )
313	END DO
314	END DO
315	END DO
316	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
317	CALL iom_put( "eken", z3d ) ! kinetic energy
318	ENDIF
319	!
320	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
321	!
322	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
323	z3d(:,:,jpk) = 0.e0
324	z2d(:,:) = 0.e0
325	DO jk = 1, jpkm1
326	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
327	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
328	END DO
329	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
330	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
331	ENDIF
332
333	IF( iom_use("u_heattr") ) THEN
334	z2d(:,:) = 0._wp
335	DO jk = 1, jpkm1
336	DO jj = 2, jpjm1
337	DO ji = fs_2, fs_jpim1 ! vector opt.
338	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
339	END DO
340	END DO
341	END DO
342	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
343	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
344	ENDIF
345
346	IF( iom_use("u_salttr") ) THEN
347	z2d(:,:) = 0.e0
348	DO jk = 1, jpkm1
349	DO jj = 2, jpjm1
350	DO ji = fs_2, fs_jpim1 ! vector opt.
351	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
352	END DO
353	END DO
354	END DO
355	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
356	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
357	ENDIF
358
359
360	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
361	z3d(:,:,jpk) = 0.e0
362	DO jk = 1, jpkm1
363	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
364	END DO
365	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
366	ENDIF
367
368	IF( iom_use("v_heattr") ) THEN
369	z2d(:,:) = 0.e0
370	DO jk = 1, jpkm1
371	DO jj = 2, jpjm1
372	DO ji = fs_2, fs_jpim1 ! vector opt.
373	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
374	END DO
375	END DO
376	END DO
377	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
378	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
379	ENDIF
380
381	IF( iom_use("v_salttr") ) THEN
382	z2d(:,:) = 0._wp
383	DO jk = 1, jpkm1
384	DO jj = 2, jpjm1
385	DO ji = fs_2, fs_jpim1 ! vector opt.
386	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
387	END DO
388	END DO
389	END DO
390	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
391	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
392	ENDIF
393
394	IF( iom_use("tosmint_"//ttype) ) THEN
395	z2d(:,:) = 0._wp
396	DO jk = 1, jpkm1
397	DO jj = 2, jpjm1
398	DO ji = fs_2, fs_jpim1 ! vector opt.
399	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
400	END DO
401	END DO
402	END DO
403	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
404	CALL iom_put( "tosmint_"//ttype, rau0 * z2d ) ! Vertical integral of temperature
405	ENDIF
406	IF( iom_use("somint_"//stype) ) THEN
407	z2d(:,:)=0._wp
408	DO jk = 1, jpkm1
409	DO jj = 2, jpjm1
410	DO ji = fs_2, fs_jpim1 ! vector opt.
411	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
412	END DO
413	END DO
414	END DO
415	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
416	CALL iom_put( "somint_"//stype, rau0 * z2d ) ! Vertical integral of salinity
417	ENDIF
418
419	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
420	!
421
422	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
423
424	!JT
425	if ( ln_pea ) THEN
426	CALL dia_pea( kt )
427	ENDIF
428
429	IF (ln_diaregmean) THEN
430	!write(,) kt,narea,'diawri before dia_regmean'
431	CALL dia_regmean( kt )
432	!write(,) kt,narea,'diawri after dia_regmean'
433	ENDIF
434
435	!JT
436
437	IF( ln_timing ) CALL timing_stop('dia_wri')
438	!
439	END SUBROUTINE dia_wri
440
441	#else
442	!!----------------------------------------------------------------------
443	!! Default option use IOIPSL library
444	!!----------------------------------------------------------------------
445
446	INTEGER FUNCTION dia_wri_alloc()
447	!!----------------------------------------------------------------------
448	INTEGER, DIMENSION(2) :: ierr
449	!!----------------------------------------------------------------------
450	IF( nn_write == -1 ) THEN
451	dia_wri_alloc = 0
452	ELSE
453	ierr = 0
454	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
455	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
456	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
457	!
458	dia_wri_alloc = MAXVAL(ierr)
459	CALL mpp_sum( 'diawri', dia_wri_alloc )
460	!
461	ENDIF
462	!
463	END FUNCTION dia_wri_alloc
464
465
466	SUBROUTINE dia_wri( kt )
467	!!---------------------------------------------------------------------
468	!! * ROUTINE dia_wri *
469	!!
470	!! ** Purpose : Standard output of opa: dynamics and tracer fields
471	!! NETCDF format is used by default
472	!!
473	!! ** Method : At the beginning of the first time step (nit000),
474	!! define all the NETCDF files and fields
475	!! At each time step call histdef to compute the mean if ncessary
476	!! Each nn_write time step, output the instantaneous or mean fields
477	!!----------------------------------------------------------------------
478	INTEGER, INTENT( in ) :: kt ! ocean time-step index
479	!
480	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
481	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
482	INTEGER :: inum = 11 ! temporary logical unit
483	INTEGER :: ji, jj, jk ! dummy loop indices
484	INTEGER :: ierr ! error code return from allocation
485	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
486	INTEGER :: jn, ierror ! local integers
487	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
488	!
489	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
490	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
491	!!----------------------------------------------------------------------
492	!
493	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
494	CALL dia_wri_state( 'output.init' )
495	ninist = 0
496	ENDIF
497	!
498	IF( nn_write == -1 ) RETURN ! we will never do any output
499	!
500	IF( ln_timing ) CALL timing_start('dia_wri')
501	!
502	! 0. Initialisation
503	! -----------------
504
505	ll_print = .FALSE. ! local variable for debugging
506	ll_print = ll_print .AND. lwp
507
508	! Define frequency of output and means
509	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
510	#if defined key_diainstant
511	zsto = nn_write * rdt
512	clop = "inst("//TRIM(clop)//")"
513	#else
514	zsto=rdt
515	clop = "ave("//TRIM(clop)//")"
516	#endif
517	zout = nn_write * rdt
518	zmax = ( nitend - nit000 + 1 ) * rdt
519
520	! Define indices of the horizontal output zoom and vertical limit storage
521	iimi = 1 ; iima = jpi
522	ijmi = 1 ; ijma = jpj
523	ipk = jpk
524
525	! define time axis
526	it = kt
527	itmod = kt - nit000 + 1
528
529
530	! 1. Define NETCDF files and fields at beginning of first time step
531	! -----------------------------------------------------------------
532
533	IF( kt == nit000 ) THEN
534
535	! Define the NETCDF files (one per grid)
536
537	! Compute julian date from starting date of the run
538	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
539	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
540	IF(lwp)WRITE(numout,*)
541	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
542	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
543	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
544	' limit storage in depth = ', ipk
545
546	! WRITE root name in date.file for use by postpro
547	IF(lwp) THEN
548	CALL dia_nam( clhstnam, nn_write,' ' )
549	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
550	WRITE(inum,*) clhstnam
551	CLOSE(inum)
552	ENDIF
553
554	! Define the T grid FILE ( nid_T )
555
556	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
557	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
558	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
559	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
560	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
561	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
562	& "m", ipk, gdept_1d, nz_T, "down" )
563	! ! Index of ocean points
564	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
565	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
566	!
567	IF( ln_icebergs ) THEN
568	!
569	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
570	!! that routine is called from nemogcm, so do it here immediately before its needed
571	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
572	CALL mpp_sum( 'diawri', ierror )
573	IF( ierror /= 0 ) THEN
574	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
575	RETURN
576	ENDIF
577	!
578	!! iceberg vertical coordinate is class number
579	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
580	& "number", nclasses, class_num, nb_T )
581	!
582	!! each class just needs the surface index pattern
583	ndim_bT = 3
584	DO jn = 1,nclasses
585	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
586	ENDDO
587	!
588	ENDIF
589
590	! Define the U grid FILE ( nid_U )
591
592	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
593	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
594	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
595	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
596	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
597	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
598	& "m", ipk, gdept_1d, nz_U, "down" )
599	! ! Index of ocean points
600	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
601	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
602
603	! Define the V grid FILE ( nid_V )
604
605	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
606	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
607	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
608	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
609	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
610	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
611	& "m", ipk, gdept_1d, nz_V, "down" )
612	! ! Index of ocean points
613	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
614	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
615
616	! Define the W grid FILE ( nid_W )
617
618	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
619	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
620	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
621	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
622	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
623	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
624	& "m", ipk, gdepw_1d, nz_W, "down" )
625
626
627	! Declare all the output fields as NETCDF variables
628
629	! !!! nid_T : 3D
630	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
631	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
632	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
633	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
634	IF( .NOT.ln_linssh ) THEN
635	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
636	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
637	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
638	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
639	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
640	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
641	ENDIF
642	! !!! nid_T : 2D
643	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
644	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
645	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
646	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
647	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
648	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
649	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
650	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
651	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
652	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
653	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
654	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
655	IF( ln_linssh ) THEN
656	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
657	& , "KgC/m2/s", & ! sosst_cd
658	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
659	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
660	& , "KgPSU/m2/s",& ! sosss_cd
661	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
662	ENDIF
663	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
664	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
665	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
666	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
667	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
668	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
669	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
670	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
671	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
672	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
673	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
674	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
675	!
676	IF( ln_icebergs ) THEN
677	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
678	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
679	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
680	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
681	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
682	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
683	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
684	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
685	IF( ln_bergdia ) THEN
686	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
687	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
688	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
689	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
690	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
691	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
692	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
693	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
694	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
695	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
696	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
697	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
698	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
699	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
700	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
701	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
702	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
703	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
704	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
705	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
706	ENDIF
707	ENDIF
708
709	IF( ln_ssr ) THEN
710	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
711	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
712	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
713	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
714	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
715	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
716	ENDIF
717
718	clmx ="l_max(only(x))" ! max index on a period
719	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
720	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
721	#if defined key_diahth
722	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
723	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
724	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
725	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
726	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
727	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
728	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
729	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
730	#endif
731
732	CALL histend( nid_T, snc4chunks=snc4set )
733
734	! !!! nid_U : 3D
735	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
736	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
737	IF( ln_wave .AND. ln_sdw) THEN
738	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
739	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
740	ENDIF
741	! !!! nid_U : 2D
742	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
743	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
744
745	CALL histend( nid_U, snc4chunks=snc4set )
746
747	! !!! nid_V : 3D
748	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
749	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
750	IF( ln_wave .AND. ln_sdw) THEN
751	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
752	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
753	ENDIF
754	! !!! nid_V : 2D
755	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
756	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
757
758	CALL histend( nid_V, snc4chunks=snc4set )
759
760	! !!! nid_W : 3D
761	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
762	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
763	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
764	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
765	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
766	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
767
768	IF( ln_zdfddm ) THEN
769	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
770	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
771	ENDIF
772
773	IF( ln_wave .AND. ln_sdw) THEN
774	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
775	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
776	ENDIF
777	! !!! nid_W : 2D
778	CALL histend( nid_W, snc4chunks=snc4set )
779
780	IF(lwp) WRITE(numout,*)
781	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
782	IF(ll_print) CALL FLUSH(numout )
783
784	ENDIF
785
786	! 2. Start writing data
787	! ---------------------
788
789	! ndex(1) est utilise ssi l'avant dernier argument est different de
790	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
791	! donne le nombre d'elements, et ndex la liste des indices a sortir
792
793	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
794	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
795	WRITE(numout,*) '~~~~~~ '
796	ENDIF
797
798	IF( .NOT.ln_linssh ) THEN
799	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
800	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
801	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
802	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
803	ELSE
804	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
805	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
806	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
807	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
808	ENDIF
809	IF( .NOT.ln_linssh ) THEN
810	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
811	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
812	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
813	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
814	ENDIF
815	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
816	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
817	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
818	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
819	! (includes virtual salt flux beneath ice
820	! in linear free surface case)
821	IF( ln_linssh ) THEN
822	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
823	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
824	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
825	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
826	ENDIF
827	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
828	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
829	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
830	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
831	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
832	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
833	!
834	IF( ln_icebergs ) THEN
835	!
836	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
837	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
838	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
839	!
840	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
841	!
842	IF( ln_bergdia ) THEN
843	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
844	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
845	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
846	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
847	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
848	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
849	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
850	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
851	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
852	!
853	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
854	ENDIF
855	ENDIF
856
857	IF( ln_ssr ) THEN
858	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
859	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
860	zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
861	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
862	ENDIF
863	! zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
864	! CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
865
866	#if defined key_diahth
867	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
868	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
869	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
870	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
871	#endif
872
873	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
874	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
875
876	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
877	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
878
879	IF( ln_zad_Aimp ) THEN
880	CALL histwrite( nid_W, "vovecrtz", it, wn + wi , ndim_T, ndex_T ) ! vert. current
881	ELSE
882	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
883	ENDIF
884	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
885	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
886	IF( ln_zdfddm ) THEN
887	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
888	ENDIF
889
890	IF( ln_wave .AND. ln_sdw ) THEN
891	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
892	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
893	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
894	ENDIF
895
896	! 3. Close all files
897	! ---------------------------------------
898	IF( kt == nitend ) THEN
899	CALL histclo( nid_T )
900	CALL histclo( nid_U )
901	CALL histclo( nid_V )
902	CALL histclo( nid_W )
903	ENDIF
904	!
905	IF( ln_timing ) CALL timing_stop('dia_wri')
906	!
907	END SUBROUTINE dia_wri
908	#endif
909
910	SUBROUTINE dia_wri_state( cdfile_name )
911	!!---------------------------------------------------------------------
912	!! * ROUTINE dia_wri_state *
913	!!
914	!! ** Purpose : create a NetCDF file named cdfile_name which contains
915	!! the instantaneous ocean state and forcing fields.
916	!! Used to find errors in the initial state or save the last
917	!! ocean state in case of abnormal end of a simulation
918	!!
919	!! ** Method : NetCDF files using ioipsl
920	!! File 'output.init.nc' is created if ninist = 1 (namelist)
921	!! File 'output.abort.nc' is created in case of abnormal job end
922	!!----------------------------------------------------------------------
923	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
924	!!
925	INTEGER :: inum
926	!!----------------------------------------------------------------------
927	!
928	IF(lwp) WRITE(numout,*)
929	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
930	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
931	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
932
933	#if defined key_si3
934	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
935	#else
936	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
937	#endif
938
939	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
940	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
941	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
942	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
943	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
944	IF( ln_zad_Aimp ) THEN
945	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi ) ! now k-velocity
946	ELSE
947	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
948	ENDIF
949	IF( ALLOCATED(ahtu) ) THEN
950	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
951	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
952	ENDIF
953	IF( ALLOCATED(ahmt) ) THEN
954	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
955	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
956	ENDIF
957	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
958	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
959	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
960	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
961	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
962	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
963	IF( .NOT.ln_linssh ) THEN
964	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
965	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
966	END IF
967	IF( ln_wave .AND. ln_sdw ) THEN
968	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
969	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
970	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
971	ENDIF
972	#if defined key_si3
973	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
974	CALL ice_wri_state( inum )
975	ENDIF
976	#endif
977	!
978	CALL iom_close( inum )
979	!
980	END SUBROUTINE dia_wri_state
981
982	!!======================================================================
983	END MODULE diawri

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

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

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_CO9_shelf_climate/src/OCE/DIA/diawri.F90 @ 15343

Download in other formats: