Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diawri.F90 @ 12721

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

Download in other formats: