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

Since March 2022 along with NEMO 4.2 release, the code development moved to a self-hosted GitLab.
This present forge is now archived and remained online for history.

diawri.F90 in branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/DIA – NEMO

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package_asm_surf_bgc/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8485

Last change on this file since 8485 was 8400, checked in by timgraham, 7 years ago

GMED ticket 335:

Merge dev_r5518_GO6_package_inc_asm into package branch to make everything easier for data assimilation
No effect on configs without data assimilation

File size: 58.9 KB

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

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

Download in other formats: