Context Navigation

← Previous Revision
Latest Revision
Next Revision →
Blame
Revision Log

diawri.F90 @ 9330

Last change on this file since 9330 was 7590, checked in by kingr, 7 years ago
Added diagnostics from operational version: Kara ML, diaopfoam, diurnal.
File size: 56.9 KB

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

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

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

Context Navigation

source: branches/UKMO/CO6_KD490_amm7_oper_fabm_hem08/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 9330

Download in other formats: