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

Context Navigation

source: branches/2012/dev_r3337_NOCS10_ICB/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 3340

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

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

Download in other formats: