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

Context Navigation

source: branches/UKMO/AMM15_v3_6_STABLE_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 8059

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

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

Download in other formats: