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

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

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

Context Navigation

source: branches/UKMO/dev_r5518_GO6_package/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 6793

Last change on this file since 6793 was 6498, checked in by timgraham, 8 years ago
Merge head of nemo_v3_6_STABLE into package branch
File size: 57.4 KB

Line
1	MODULE diawri
2	!!======================================================================
3	!! * MODULE diawri *
4	!! Ocean diagnostics : write ocean output files
5	!!=====================================================================
6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
7	!! 4.0 ! 1991-11 (G. Madec)
8	!! ! 1992-06 (M. Imbard) correction restart file
9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
10	!! ! 1993-03 (M. Imbard) suppress writibm
11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
15	!! - ! 2002-09 (G. Madec) F90: Free form and module
16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
19	!!----------------------------------------------------------------------
20
21	!!----------------------------------------------------------------------
22	!! dia_wri : create the standart output files
23	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24	!!----------------------------------------------------------------------
25	USE oce ! ocean dynamics and tracers
26	USE dom_oce ! ocean space and time domain
27	USE dynadv, ONLY: ln_dynadv_vec
28	USE zdf_oce ! ocean vertical physics
29	USE ldftra_oce ! ocean active tracers: lateral physics
30	USE ldfdyn_oce ! ocean dynamics: lateral physics
31	USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
32	USE sol_oce ! solver variables
33	USE sbc_oce ! Surface boundary condition: ocean fields
34	USE sbc_ice ! Surface boundary condition: ice fields
35	USE icb_oce ! Icebergs
36	USE icbdia ! Iceberg budgets
37	USE sbcssr ! restoring term toward SST/SSS climatology
38	USE phycst ! physical constants
39	USE zdfmxl ! mixed layer
40	USE dianam ! build name of file (routine)
41	USE zdftke ! vertical physics: one-equation scheme
42	USE zdfddm ! vertical physics: double diffusion
43	USE diahth ! thermocline diagnostics
44	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
45	USE in_out_manager ! I/O manager
46	USE diadimg ! dimg direct access file format output
47	USE iom
48	USE ioipsl
49	USE dynspg_oce, ONLY: un_adv, vn_adv ! barotropic velocities
50
51	#if defined key_lim2
52	USE limwri_2
53	#elif defined key_lim3
54	USE limwri
55	#endif
56	USE lib_mpp ! MPP library
57	USE timing ! preformance summary
58	USE wrk_nemo ! working array
59
60	IMPLICIT NONE
61	PRIVATE
62
63	PUBLIC dia_wri ! routines called by step.F90
64	PUBLIC dia_wri_state
65	PUBLIC dia_wri_alloc ! Called by nemogcm module
66
67	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
68	INTEGER :: nb_T , ndim_bT ! grid_T file
69	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
70	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
71	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
72	INTEGER :: ndex(1) ! ???
73	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
74	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
75	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
76
77	!! * Substitutions
78	# include "zdfddm_substitute.h90"
79	# include "domzgr_substitute.h90"
80	# include "vectopt_loop_substitute.h90"
81	!!----------------------------------------------------------------------
82	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
83	!! $Id$
84	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
85	!!----------------------------------------------------------------------
86	CONTAINS
87
88	INTEGER FUNCTION dia_wri_alloc()
89	!!----------------------------------------------------------------------
90	INTEGER, DIMENSION(2) :: ierr
91	!!----------------------------------------------------------------------
92	ierr = 0
93	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
94	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
95	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
96	!
97	dia_wri_alloc = MAXVAL(ierr)
98	IF( lk_mpp ) CALL mpp_sum( dia_wri_alloc )
99	!
100	END FUNCTION dia_wri_alloc
101
102	#if defined key_dimgout
103	!!----------------------------------------------------------------------
104	!! 'key_dimgout' DIMG output file
105	!!----------------------------------------------------------------------
106	# include "diawri_dimg.h90"
107
108	#else
109	!!----------------------------------------------------------------------
110	!! Default option NetCDF output file
111	!!----------------------------------------------------------------------
112	# if defined key_iomput
113	!!----------------------------------------------------------------------
114	!! 'key_iomput' use IOM library
115	!!----------------------------------------------------------------------
116
117	SUBROUTINE dia_wri( kt )
118	!!---------------------------------------------------------------------
119	!! * ROUTINE dia_wri *
120	!!
121	!! ** Purpose : Standard output of opa: dynamics and tracer fields
122	!! NETCDF format is used by default
123	!!
124	!! ** Method : use iom_put
125	!!----------------------------------------------------------------------
126	!!
127	INTEGER, INTENT( in ) :: kt ! ocean time-step index
128	!!
129	INTEGER :: ji, jj, jk ! dummy loop indices
130	INTEGER :: jkbot !
131	REAL(wp) :: zztmp, zztmpx, zztmpy !
132	!!
133	REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace
134	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace
135	!!----------------------------------------------------------------------
136	!
137	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
138	!
139	CALL wrk_alloc( jpi , jpj , z2d )
140	CALL wrk_alloc( jpi , jpj, jpk , z3d )
141	!
142	! Output the initial state and forcings
143	IF( ninist == 1 ) THEN
144	CALL dia_wri_state( 'output.init', kt )
145	ninist = 0
146	ENDIF
147
148	! Output of initial vertical scale factor
149	CALL iom_put("e3t_0", e3t_0(:,:,:) )
150	CALL iom_put("e3u_0", e3t_0(:,:,:) )
151	CALL iom_put("e3v_0", e3t_0(:,:,:) )
152	!
153	CALL iom_put( "e3t" , fse3t_n(:,:,:) )
154	CALL iom_put( "e3u" , fse3u_n(:,:,:) )
155	CALL iom_put( "e3v" , fse3v_n(:,:,:) )
156	CALL iom_put( "e3w" , fse3w_n(:,:,:) )
157	IF( iom_use("e3tdef") ) &
158	CALL iom_put( "e3tdef" , ( ( fse3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
159
160
161	CALL iom_put( "ssh" , sshn ) ! sea surface height
162
163	CALL iom_put( "toce", tsn(:,:,:,jp_tem) ) ! 3D temperature
164	CALL iom_put( "sst", tsn(:,:,1,jp_tem) ) ! surface temperature
165	IF ( iom_use("sbt") ) THEN
166	DO jj = 1, jpj
167	DO ji = 1, jpi
168	jkbot = mbkt(ji,jj)
169	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_tem)
170	END DO
171	END DO
172	CALL iom_put( "sbt", z2d ) ! bottom temperature
173	ENDIF
174
175	CALL iom_put( "soce", tsn(:,:,:,jp_sal) ) ! 3D salinity
176	CALL iom_put( "sss", tsn(:,:,1,jp_sal) ) ! surface salinity
177	IF ( iom_use("sbs") ) THEN
178	DO jj = 1, jpj
179	DO ji = 1, jpi
180	jkbot = mbkt(ji,jj)
181	z2d(ji,jj) = tsn(ji,jj,jkbot,jp_sal)
182	END DO
183	END DO
184	CALL iom_put( "sbs", z2d ) ! bottom salinity
185	ENDIF
186
187	IF ( iom_use("taubot") ) THEN ! bottom stress
188	z2d(:,:) = 0._wp
189	DO jj = 2, jpjm1
190	DO ji = fs_2, fs_jpim1 ! vector opt.
191	zztmpx = ( bfrua(ji ,jj) * un(ji ,jj,mbku(ji ,jj)) &
192	& + bfrua(ji-1,jj) * un(ji-1,jj,mbku(ji-1,jj)) )
193	zztmpy = ( bfrva(ji, jj) * vn(ji,jj ,mbkv(ji,jj )) &
194	& + bfrva(ji,jj-1) * vn(ji,jj-1,mbkv(ji,jj-1)) )
195	z2d(ji,jj) = rau0 * SQRT( zztmpx * zztmpx + zztmpy * zztmpy ) * tmask(ji,jj,1)
196	!
197	ENDDO
198	ENDDO
199	CALL lbc_lnk( z2d, 'T', 1. )
200	CALL iom_put( "taubot", z2d )
201	ENDIF
202
203	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
204	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
205	IF ( iom_use("sbu") ) THEN
206	DO jj = 1, jpj
207	DO ji = 1, jpi
208	jkbot = mbku(ji,jj)
209	z2d(ji,jj) = un(ji,jj,jkbot)
210	END DO
211	END DO
212	CALL iom_put( "sbu", z2d ) ! bottom i-current
213	ENDIF
214	#if defined key_dynspg_ts
215	CALL iom_put( "ubar", un_adv(:,:) ) ! barotropic i-current
216	#else
217	CALL iom_put( "ubar", un_b(:,:) ) ! barotropic i-current
218	#endif
219
220	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
221	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
222	IF ( iom_use("sbv") ) THEN
223	DO jj = 1, jpj
224	DO ji = 1, jpi
225	jkbot = mbkv(ji,jj)
226	z2d(ji,jj) = vn(ji,jj,jkbot)
227	END DO
228	END DO
229	CALL iom_put( "sbv", z2d ) ! bottom j-current
230	ENDIF
231	#if defined key_dynspg_ts
232	CALL iom_put( "vbar", vn_adv(:,:) ) ! barotropic j-current
233	#else
234	CALL iom_put( "vbar", vn_b(:,:) ) ! barotropic j-current
235	#endif
236
237	CALL iom_put( "woce", wn ) ! vertical velocity
238	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
239	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
240	z2d(:,:) = rau0 * e12t(:,:)
241	DO jk = 1, jpk
242	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
243	END DO
244	CALL iom_put( "w_masstr" , z3d )
245	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
246	ENDIF
247
248	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
249	CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef.
250	IF( lk_zdftke ) THEN
251	CALL iom_put( "tke" , en ) ! TKE budget: Turbulent Kinetic Energy
252	CALL iom_put( "tke_niw" , e_niw ) ! TKE budget: Near-inertial waves
253	ENDIF
254	CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef. (useful only with key_zdfddm)
255	! Log of eddy diff coef
256	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt (:,:,:) ) ) )
257	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, fsavs(:,:,:) ) ) )
258
259	IF ( iom_use("sstgrad") .OR. iom_use("sstgrad2") ) THEN
260	DO jj = 2, jpjm1 ! sst gradient
261	DO ji = fs_2, fs_jpim1 ! vector opt.
262	zztmp = tsn(ji,jj,1,jp_tem)
263	zztmpx = ( tsn(ji+1,jj ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) / e1u(ji-1,jj )
264	zztmpy = ( tsn(ji ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) / e2v(ji ,jj-1)
265	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
266	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
267	END DO
268	END DO
269	CALL lbc_lnk( z2d, 'T', 1. )
270	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
271	z2d(:,:) = SQRT( z2d(:,:) )
272	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
273	ENDIF
274
275	! clem: heat and salt content
276	IF( iom_use("heatc") ) THEN
277	z2d(:,:) = 0._wp
278	DO jk = 1, jpkm1
279	DO jj = 1, jpj
280	DO ji = 1, jpi
281	z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
282	END DO
283	END DO
284	END DO
285	CALL iom_put( "heatc", (rau0 * rcp) * z2d ) ! vertically integrated heat content (J/m2)
286	ENDIF
287
288	IF( iom_use("saltc") ) THEN
289	z2d(:,:) = 0._wp
290	DO jk = 1, jpkm1
291	DO jj = 1, jpj
292	DO ji = 1, jpi
293	z2d(ji,jj) = z2d(ji,jj) + fse3t(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
294	END DO
295	END DO
296	END DO
297	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
298	ENDIF
299	!
300	IF ( iom_use("eken") ) THEN
301	rke(:,:,jk) = 0._wp ! kinetic energy
302	DO jk = 1, jpkm1
303	DO jj = 2, jpjm1
304	DO ji = fs_2, fs_jpim1 ! vector opt.
305	zztmp = 1._wp / ( e1e2t(ji,jj) * fse3t(ji,jj,jk) )
306	zztmpx = 0.5 * ( un(ji-1,jj,jk) * un(ji-1,jj,jk) * e2u(ji-1,jj) * fse3u(ji-1,jj,jk) &
307	& + un(ji ,jj,jk) * un(ji ,jj,jk) * e2u(ji ,jj) * fse3u(ji ,jj,jk) ) &
308	& * zztmp
309	!
310	zztmpy = 0.5 * ( vn(ji,jj-1,jk) * vn(ji,jj-1,jk) * e1v(ji,jj-1) * fse3v(ji,jj-1,jk) &
311	& + vn(ji,jj ,jk) * vn(ji,jj ,jk) * e1v(ji,jj ) * fse3v(ji,jj ,jk) ) &
312	& * zztmp
313	!
314	rke(ji,jj,jk) = 0.5_wp * ( zztmpx + zztmpy )
315	!
316	ENDDO
317	ENDDO
318	ENDDO
319	CALL lbc_lnk( rke, 'T', 1. )
320	CALL iom_put( "eken", rke )
321	ENDIF
322	!
323	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
324	!
325	IF( iom_use("u_masstr") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
326	z3d(:,:,jpk) = 0.e0
327	DO jk = 1, jpkm1
328	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk) * umask(:,:,jk)
329	END DO
330	CALL iom_put( "u_masstr", z3d ) ! mass transport in i-direction
331	ENDIF
332
333	IF( iom_use("u_heattr") ) THEN
334	z2d(:,:) = 0.e0
335	DO jk = 1, jpkm1
336	DO jj = 2, jpjm1
337	DO ji = fs_2, fs_jpim1 ! vector opt.
338	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
339	END DO
340	END DO
341	END DO
342	CALL lbc_lnk( z2d, 'U', -1. )
343	CALL iom_put( "u_heattr", (0.5 * rcp) * z2d ) ! heat transport in i-direction
344	ENDIF
345
346	IF( iom_use("u_salttr") ) THEN
347	z2d(:,:) = 0.e0
348	DO jk = 1, jpkm1
349	DO jj = 2, jpjm1
350	DO ji = fs_2, fs_jpim1 ! vector opt.
351	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
352	END DO
353	END DO
354	END DO
355	CALL lbc_lnk( z2d, 'U', -1. )
356	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
357	ENDIF
358
359
360	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
361	z3d(:,:,jpk) = 0.e0
362	DO jk = 1, jpkm1
363	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk) * vmask(:,:,jk)
364	END DO
365	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
366	ENDIF
367
368	IF( iom_use("v_heattr") ) THEN
369	z2d(:,:) = 0.e0
370	DO jk = 1, jpkm1
371	DO jj = 2, jpjm1
372	DO ji = fs_2, fs_jpim1 ! vector opt.
373	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
374	END DO
375	END DO
376	END DO
377	CALL lbc_lnk( z2d, 'V', -1. )
378	CALL iom_put( "v_heattr", (0.5 * rcp) * z2d ) ! heat transport in j-direction
379	ENDIF
380
381	IF( iom_use("v_salttr") ) THEN
382	z2d(:,:) = 0.e0
383	DO jk = 1, jpkm1
384	DO jj = 2, jpjm1
385	DO ji = fs_2, fs_jpim1 ! vector opt.
386	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
387	END DO
388	END DO
389	END DO
390	CALL lbc_lnk( z2d, 'V', -1. )
391	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
392	ENDIF
393	!
394	CALL wrk_dealloc( jpi , jpj , z2d )
395	CALL wrk_dealloc( jpi , jpj, jpk , z3d )
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.

Download in other formats: