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/2014/dev_r4627_COMODO_UPW2D/NEMOGCM/CONFIG/UPW2D_eo_NOadv_NOvvl/MY_SRC – NEMO

Context Navigation

source: branches/2014/dev_r4627_COMODO_UPW2D/NEMOGCM/CONFIG/UPW2D_eo_NOadv_NOvvl/MY_SRC/diawri.F90 @ 4648

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

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

Download in other formats: