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 NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/UKMO/NEMO_4.0.4_GO8_package/src/OCE/DIA/diawri.F90 @ 14342

Last change on this file since 14342 was 14342, checked in by davestorkey, 3 years ago
UKMO/NEMO_4.0.4_GO8_package: Add Gent-McWilliams? variants.
File size: 51.2 KB

Line
1	MODULE diawri
2	!!======================================================================
3	!! * MODULE diawri *
4	!! Ocean diagnostics : write ocean output files
5	!!=====================================================================
6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
7	!! 4.0 ! 1991-11 (G. Madec)
8	!! ! 1992-06 (M. Imbard) correction restart file
9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
10	!! ! 1993-03 (M. Imbard) suppress writibm
11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
15	!! - ! 2002-09 (G. Madec) F90: Free form and module
16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
19	!! 3.7 ! 2014-01 (G. Madec) remove eddy induced velocity from no-IOM output
20	!! ! change name of output variables in dia_wri_state
21	!!----------------------------------------------------------------------
22
23	!!----------------------------------------------------------------------
24	!! dia_wri : create the standart output files
25	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
26	!!----------------------------------------------------------------------
27	USE oce ! ocean dynamics and tracers
28	USE dom_oce ! ocean space and time domain
29	USE phycst ! physical constants
30	USE dianam ! build name of file (routine)
31	USE diahth ! thermocline diagnostics
32	USE dynadv , ONLY: ln_dynadv_vec
33	USE icb_oce ! Icebergs
34	USE icbdia ! Iceberg budgets
35	USE ldftra ! lateral physics: eddy diffusivity coef.
36	USE ldfdyn ! lateral physics: eddy viscosity coef.
37	USE sbc_oce ! Surface boundary condition: ocean fields
38	USE sbc_ice ! Surface boundary condition: ice fields
39	USE sbcssr ! restoring term toward SST/SSS climatology
40	USE sbcwave ! wave parameters
41	USE wet_dry ! wetting and drying
42	USE zdf_oce ! ocean vertical physics
43	USE zdfdrg ! ocean vertical physics: top/bottom friction
44	USE zdfmxl ! mixed layer
45	!
46	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
47	USE in_out_manager ! I/O manager
48	USE dia25h ! 25h Mean output
49	USE iom !
50	USE ioipsl !
51	USE eosbn2
52	#if defined key_si3
53	USE ice
54	USE icewri
55	#endif
56	USE lib_mpp ! MPP library
57	USE timing ! preformance summary
58	USE diurnal_bulk ! diurnal warm layer
59	USE cool_skin ! Cool skin
60
61	IMPLICIT NONE
62	PRIVATE
63
64	PUBLIC dia_wri ! routines called by step.F90
65	PUBLIC dia_wri_state
66	PUBLIC dia_wri_alloc ! Called by nemogcm module
67
68	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
69	INTEGER :: nb_T , ndim_bT ! grid_T file
70	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
71	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
72	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
73	INTEGER :: ndex(1) ! ???
74	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
75	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
76	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_bT
77
78	!! * Substitutions
79	# include "vectopt_loop_substitute.h90"
80	!!----------------------------------------------------------------------
81	!! NEMO/OCE 4.0 , NEMO Consortium (2018)
82	!! $Id$
83	!! Software governed by the CeCILL license (see ./LICENSE)
84	!!----------------------------------------------------------------------
85	CONTAINS
86
87	#if defined key_iomput
88	!!----------------------------------------------------------------------
89	!! 'key_iomput' use IOM library
90	!!----------------------------------------------------------------------
91	INTEGER FUNCTION dia_wri_alloc()
92	!
93	dia_wri_alloc = 0
94	!
95	END FUNCTION dia_wri_alloc
96
97
98	SUBROUTINE dia_wri( kt )
99	!!---------------------------------------------------------------------
100	!! * ROUTINE dia_wri *
101	!!
102	!! ** Purpose : Standard output of opa: dynamics and tracer fields
103	!! NETCDF format is used by default
104	!!
105	!! ** Method : use iom_put
106	!!----------------------------------------------------------------------
107	INTEGER, INTENT( in ) :: kt ! ocean time-step index
108	!!
109	INTEGER :: ji, jj, jk ! dummy loop indices
110	INTEGER :: ikbot ! local integer
111	REAL(wp):: zztmp , zztmpx ! local scalar
112	REAL(wp):: zztmp2, zztmpy ! - -
113	REAL(wp), DIMENSION(jpi,jpj) :: z2d ! 2D workspace
114	REAL(wp), DIMENSION(jpi,jpj,jpk) :: z3d ! 3D workspace
115	CHARACTER(len=4),SAVE :: ttype , stype ! temperature and salinity type
116	!!----------------------------------------------------------------------
117	!
118	IF( kt == nit000 ) THEN
119	IF( ln_TEOS10 ) THEN
120	IF ( iom_use("toce_pot") .OR. iom_use("soce_pra") .OR. iom_use("sst_pot") .OR. iom_use("sss_pra") &
121	& .OR. iom_use("sbt_pot") .OR. iom_use("sbs_pra") .OR. iom_use("sstgrad_pot") .OR. iom_use("sstgrad2_pot") &
122	& .OR. iom_use("tosmint_pot") .OR. iom_use("somint_pra")) THEN
123	CALL ctl_stop( 'diawri: potential temperature and practical salinity not available with ln_TEOS10' )
124	ELSE
125	ttype='con' ; stype='abs' ! teos-10 using conservative temperature and absolute salinity
126	ENDIF
127	ELSE IF( ln_EOS80 ) THEN
128	IF ( iom_use("toce_con") .OR. iom_use("soce_abs") .OR. iom_use("sst_con") .OR. iom_use("sss_abs") &
129	& .OR. iom_use("sbt_con") .OR. iom_use("sbs_abs") .OR. iom_use("sstgrad_con") .OR. iom_use("sstgrad2_con") &
130	& .OR. iom_use("tosmint_con") .OR. iom_use("somint_abs")) THEN
131	CALL ctl_stop( 'diawri: conservative temperature and absolute salinity not available with ln_EOS80' )
132	ELSE
133	ttype='pot' ; stype='pra' ! eos-80 using potential temperature and practical salinity
134	ENDIF
135	ELSE IF ( ln_SEOS) THEN
136	ttype='seos' ; stype='seos' ! seos using Simplified Equation of state
137	ENDIF
138	ENDIF
139
140	IF( ln_timing ) CALL timing_start('dia_wri')
141	!
142	! Output the initial state and forcings
143	IF( ninist == 1 ) THEN
144	CALL dia_wri_state( 'output.init' )
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", e3u_0(:,:,:) )
151	CALL iom_put("e3v_0", e3v_0(:,:,:) )
152	!
153	CALL iom_put( "e3t" , e3t_n(:,:,:) )
154	CALL iom_put( "e3u" , e3u_n(:,:,:) )
155	CALL iom_put( "e3v" , e3v_n(:,:,:) )
156	CALL iom_put( "e3w" , e3w_n(:,:,:) )
157	IF( iom_use("e3tdef") ) &
158	CALL iom_put( "e3tdef" , ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2 )
159
160	IF( ll_wd ) THEN
161	CALL iom_put( "ssh" , (sshn+ssh_ref)*tmask(:,:,1) ) ! sea surface height (brought back to the reference used for wetting and drying)
162	ELSE
163	CALL iom_put( "ssh" , sshn ) ! sea surface height
164	ENDIF
165
166	IF( iom_use("wetdep") ) & ! wet depth
167	CALL iom_put( "wetdep" , ht_0(:,:) + sshn(:,:) )
168
169	CALL iom_put( "toce_"//ttype, tsn(:,:,:,jp_tem) ) ! 3D temperature
170	CALL iom_put( "sst_"//ttype, tsn(:,:,1,jp_tem) ) ! surface temperature
171	IF ( iom_use("sbt_"//ttype) ) THEN
172	DO jj = 1, jpj
173	DO ji = 1, jpi
174	ikbot = mbkt(ji,jj)
175	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_tem)
176	END DO
177	END DO
178	CALL iom_put( "sbt_"//ttype, z2d ) ! bottom temperature
179	ENDIF
180
181	CALL iom_put( "soce_"//stype, tsn(:,:,:,jp_sal) ) ! 3D salinity
182	CALL iom_put( "sss_"//stype, tsn(:,:,1,jp_sal) ) ! surface salinity
183	IF ( iom_use("sbs_"//stype) ) THEN
184	DO jj = 1, jpj
185	DO ji = 1, jpi
186	ikbot = mbkt(ji,jj)
187	z2d(ji,jj) = tsn(ji,jj,ikbot,jp_sal)
188	END DO
189	END DO
190	CALL iom_put( "sbs_"//stype, z2d ) ! bottom salinity
191	ENDIF
192
193	CALL iom_put( "rhop", rhop(:,:,:) ) ! 3D potential density (sigma0)
194
195	IF ( iom_use("taubot") ) THEN ! bottom stress
196	zztmp = rau0 * 0.25
197	z2d(:,:) = 0._wp
198	DO jj = 2, jpjm1
199	DO ji = fs_2, fs_jpim1 ! vector opt.
200	zztmp2 = ( ( rCdU_bot(ji+1,jj)+rCdU_bot(ji ,jj) ) * un(ji ,jj,mbku(ji ,jj)) )**2 &
201	& + ( ( rCdU_bot(ji ,jj)+rCdU_bot(ji-1,jj) ) * un(ji-1,jj,mbku(ji-1,jj)) )**2 &
202	& + ( ( rCdU_bot(ji,jj+1)+rCdU_bot(ji,jj ) ) * vn(ji,jj ,mbkv(ji,jj )) )**2 &
203	& + ( ( rCdU_bot(ji,jj )+rCdU_bot(ji,jj-1) ) * vn(ji,jj-1,mbkv(ji,jj-1)) )**2
204	z2d(ji,jj) = zztmp * SQRT( zztmp2 ) * tmask(ji,jj,1)
205	!
206	END DO
207	END DO
208	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
209	CALL iom_put( "taubot", z2d )
210	ENDIF
211
212	CALL iom_put( "uoce", un(:,:,:) ) ! 3D i-current
213	CALL iom_put( "ssu", un(:,:,1) ) ! surface i-current
214	IF ( iom_use("sbu") ) THEN
215	DO jj = 1, jpj
216	DO ji = 1, jpi
217	ikbot = mbku(ji,jj)
218	z2d(ji,jj) = un(ji,jj,ikbot)
219	END DO
220	END DO
221	CALL iom_put( "sbu", z2d ) ! bottom i-current
222	ENDIF
223
224	CALL iom_put( "voce", vn(:,:,:) ) ! 3D j-current
225	CALL iom_put( "ssv", vn(:,:,1) ) ! surface j-current
226	IF ( iom_use("sbv") ) THEN
227	DO jj = 1, jpj
228	DO ji = 1, jpi
229	ikbot = mbkv(ji,jj)
230	z2d(ji,jj) = vn(ji,jj,ikbot)
231	END DO
232	END DO
233	CALL iom_put( "sbv", z2d ) ! bottom j-current
234	ENDIF
235
236	IF( ln_zad_Aimp ) wn = wn + wi ! Recombine explicit and implicit parts of vertical velocity for diagnostic output
237	!
238	CALL iom_put( "woce", wn ) ! vertical velocity
239	IF( iom_use('w_masstr') .OR. iom_use('w_masstr2') ) THEN ! vertical mass transport & its square value
240	! Caution: in the VVL case, it only correponds to the baroclinic mass transport.
241	z2d(:,:) = rau0 * e1e2t(:,:)
242	DO jk = 1, jpk
243	z3d(:,:,jk) = wn(:,:,jk) * z2d(:,:)
244	END DO
245	CALL iom_put( "w_masstr" , z3d )
246	IF( iom_use('w_masstr2') ) CALL iom_put( "w_masstr2", z3d(:,:,:) * z3d(:,:,:) )
247	ENDIF
248	!
249	IF( ln_zad_Aimp ) wn = wn - wi ! Remove implicit part of vertical velocity that was added for diagnostic output
250
251	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
252	CALL iom_put( "avs" , avs ) ! S vert. eddy diff. coef.
253	CALL iom_put( "avm" , avm ) ! T vert. eddy visc. coef.
254
255	IF( iom_use('logavt') ) CALL iom_put( "logavt", LOG( MAX( 1.e-20_wp, avt(:,:,:) ) ) )
256	IF( iom_use('logavs') ) CALL iom_put( "logavs", LOG( MAX( 1.e-20_wp, avs(:,:,:) ) ) )
257
258	IF ( iom_use("sstgrad_"//ttype) .OR. iom_use("sstgrad2_"//ttype) ) THEN
259	DO jj = 2, jpjm1 ! sst gradient
260	DO ji = fs_2, fs_jpim1 ! vector opt.
261	zztmp = tsn(ji,jj,1,jp_tem)
262	zztmpx = ( tsn(ji+1,jj,1,jp_tem) - zztmp ) * r1_e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) * r1_e1u(ji-1,jj)
263	zztmpy = ( tsn(ji,jj+1,1,jp_tem) - zztmp ) * r1_e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) * r1_e2v(ji,jj-1)
264	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
265	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
266	END DO
267	END DO
268	CALL lbc_lnk( 'diawri', z2d, 'T', 1. )
269	CALL iom_put( "sstgrad2_"//ttype, z2d ) ! square of module of sst gradient
270	z2d(:,:) = SQRT( z2d(:,:) )
271	CALL iom_put( "sstgrad_"//ttype , z2d ) ! module of sst gradient
272	ENDIF
273
274	! heat and salt contents
275	IF( iom_use("heatc") ) THEN
276	z2d(:,:) = 0._wp
277	DO jk = 1, jpkm1
278	DO jj = 1, jpj
279	DO ji = 1, jpi
280	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem) * tmask(ji,jj,jk)
281	END DO
282	END DO
283	END DO
284	CALL iom_put( "heatc", rau0_rcp * z2d ) ! vertically integrated heat content (J/m2)
285	ENDIF
286
287	IF( iom_use("saltc") ) THEN
288	z2d(:,:) = 0._wp
289	DO jk = 1, jpkm1
290	DO jj = 1, jpj
291	DO ji = 1, jpi
292	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal) * tmask(ji,jj,jk)
293	END DO
294	END DO
295	END DO
296	CALL iom_put( "saltc", rau0 * z2d ) ! vertically integrated salt content (PSU*kg/m2)
297	ENDIF
298	!
299	IF ( iom_use("eken") ) THEN
300	z3d(:,:,jpk) = 0._wp
301	DO jk = 1, jpkm1
302	DO jj = 2, jpjm1
303	DO ji = fs_2, fs_jpim1 ! vector opt.
304	zztmp = 0.25_wp * r1_e1e2t(ji,jj) / e3t_n(ji,jj,jk)
305	z3d(ji,jj,jk) = zztmp * ( un(ji-1,jj,jk)*2 e2u(ji-1,jj) * e3u_n(ji-1,jj,jk) &
306	& + un(ji ,jj,jk)*2 e2u(ji ,jj) * e3u_n(ji ,jj,jk) &
307	& + vn(ji,jj-1,jk)*2 e1v(ji,jj-1) * e3v_n(ji,jj-1,jk) &
308	& + vn(ji,jj ,jk)*2 e1v(ji,jj ) * e3v_n(ji,jj ,jk) )
309	END DO
310	END DO
311	END DO
312	CALL lbc_lnk( 'diawri', z3d, 'T', 1. )
313	CALL iom_put( "eken", z3d ) ! kinetic energy
314	ENDIF
315	!
316	CALL iom_put( "hdiv", hdivn ) ! Horizontal divergence
317	!
318	IF( iom_use("u_masstr") .OR. iom_use("u_masstr_vint") .OR. iom_use("u_heattr") .OR. iom_use("u_salttr") ) THEN
319	z3d(:,:,jpk) = 0.e0
320	z2d(:,:) = 0.e0
321	DO jk = 1, jpkm1
322	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * e3u_n(:,:,jk) * umask(:,:,jk)
323	z2d(:,:) = z2d(:,:) + z3d(:,:,jk)
324	END DO
325	CALL iom_put( "u_masstr" , z3d ) ! mass transport in i-direction
326	CALL iom_put( "u_masstr_vint", z2d ) ! mass transport in i-direction vertical sum
327	ENDIF
328
329	IF( iom_use("u_heattr") ) THEN
330	z2d(:,:) = 0._wp
331	DO jk = 1, jpkm1
332	DO jj = 2, jpjm1
333	DO ji = fs_2, fs_jpim1 ! vector opt.
334	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
335	END DO
336	END DO
337	END DO
338	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
339	CALL iom_put( "u_heattr", 0.5rcp z2d ) ! heat transport in i-direction
340	ENDIF
341
342	IF( iom_use("u_salttr") ) THEN
343	z2d(:,:) = 0.e0
344	DO jk = 1, jpkm1
345	DO jj = 2, jpjm1
346	DO ji = fs_2, fs_jpim1 ! vector opt.
347	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji+1,jj,jk,jp_sal) )
348	END DO
349	END DO
350	END DO
351	CALL lbc_lnk( 'diawri', z2d, 'U', -1. )
352	CALL iom_put( "u_salttr", 0.5 * z2d ) ! heat transport in i-direction
353	ENDIF
354
355
356	IF( iom_use("v_masstr") .OR. iom_use("v_heattr") .OR. iom_use("v_salttr") ) THEN
357	z3d(:,:,jpk) = 0.e0
358	DO jk = 1, jpkm1
359	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * e3v_n(:,:,jk) * vmask(:,:,jk)
360	END DO
361	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
362	ENDIF
363
364	IF( iom_use("v_heattr") ) THEN
365	z2d(:,:) = 0.e0
366	DO jk = 1, jpkm1
367	DO jj = 2, jpjm1
368	DO ji = fs_2, fs_jpim1 ! vector opt.
369	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
370	END DO
371	END DO
372	END DO
373	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
374	CALL iom_put( "v_heattr", 0.5rcp z2d ) ! heat transport in j-direction
375	ENDIF
376
377	IF( iom_use("v_salttr") ) THEN
378	z2d(:,:) = 0._wp
379	DO jk = 1, jpkm1
380	DO jj = 2, jpjm1
381	DO ji = fs_2, fs_jpim1 ! vector opt.
382	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * ( tsn(ji,jj,jk,jp_sal) + tsn(ji,jj+1,jk,jp_sal) )
383	END DO
384	END DO
385	END DO
386	CALL lbc_lnk( 'diawri', z2d, 'V', -1. )
387	CALL iom_put( "v_salttr", 0.5 * z2d ) ! heat transport in j-direction
388	ENDIF
389
390	IF( (.NOT.l_ldfeiv_time) .AND. ( iom_use('RossRad') .OR. iom_use('RossRadlim') &
391	& .OR. iom_use('TclinicR') .OR. iom_use('RR_GG') &
392	& .OR. iom_use('aeiu_2d') .OR. iom_use('aeiv_2d') ) ) THEN
393	CALL ldf_eiv(kt, 75.0, z2d, z3d(:,:,1))
394	CALL iom_put('aeiu_2d', z2d)
395	CALL iom_put('aeiv_2d', z3d(:,:,1))
396	ENDIF
397
398	IF( iom_use("tosmint_"//ttype) ) THEN
399	z2d(:,:) = 0._wp
400	DO jk = 1, jpkm1
401	DO jj = 2, jpjm1
402	DO ji = fs_2, fs_jpim1 ! vector opt.
403	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_tem)
404	END DO
405	END DO
406	END DO
407	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
408	CALL iom_put( "tosmint_"//ttype, rau0 * z2d ) ! Vertical integral of temperature
409	ENDIF
410	IF( iom_use("somint_"//stype) ) THEN
411	z2d(:,:)=0._wp
412	DO jk = 1, jpkm1
413	DO jj = 2, jpjm1
414	DO ji = fs_2, fs_jpim1 ! vector opt.
415	z2d(ji,jj) = z2d(ji,jj) + e3t_n(ji,jj,jk) * tsn(ji,jj,jk,jp_sal)
416	END DO
417	END DO
418	END DO
419	CALL lbc_lnk( 'diawri', z2d, 'T', -1. )
420	CALL iom_put( "somint_"//stype, rau0 * z2d ) ! Vertical integral of salinity
421	ENDIF
422
423	CALL iom_put( "bn2", rn2 ) ! Brunt-Vaisala buoyancy frequency (N^2)
424	!
425
426	IF (ln_dia25h) CALL dia_25h( kt ) ! 25h averaging
427
428	IF( ln_timing ) CALL timing_stop('dia_wri')
429	!
430	END SUBROUTINE dia_wri
431
432	#else
433	!!----------------------------------------------------------------------
434	!! Default option use IOIPSL library
435	!!----------------------------------------------------------------------
436
437	INTEGER FUNCTION dia_wri_alloc()
438	!!----------------------------------------------------------------------
439	INTEGER, DIMENSION(2) :: ierr
440	!!----------------------------------------------------------------------
441	IF( nn_write == -1 ) THEN
442	dia_wri_alloc = 0
443	ELSE
444	ierr = 0
445	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
446	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
447	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
448	!
449	dia_wri_alloc = MAXVAL(ierr)
450	CALL mpp_sum( 'diawri', dia_wri_alloc )
451	!
452	ENDIF
453	!
454	END FUNCTION dia_wri_alloc
455
456
457	SUBROUTINE dia_wri( kt )
458	!!---------------------------------------------------------------------
459	!! * ROUTINE dia_wri *
460	!!
461	!! ** Purpose : Standard output of opa: dynamics and tracer fields
462	!! NETCDF format is used by default
463	!!
464	!! ** Method : At the beginning of the first time step (nit000),
465	!! define all the NETCDF files and fields
466	!! At each time step call histdef to compute the mean if ncessary
467	!! Each nn_write time step, output the instantaneous or mean fields
468	!!----------------------------------------------------------------------
469	INTEGER, INTENT( in ) :: kt ! ocean time-step index
470	!
471	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
472	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
473	INTEGER :: inum = 11 ! temporary logical unit
474	INTEGER :: ji, jj, jk ! dummy loop indices
475	INTEGER :: ierr ! error code return from allocation
476	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
477	INTEGER :: jn, ierror ! local integers
478	REAL(wp) :: zsto, zout, zmax, zjulian ! local scalars
479	!
480	REAL(wp), DIMENSION(jpi,jpj) :: zw2d ! 2D workspace
481	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! 3D workspace
482	!!----------------------------------------------------------------------
483	!
484	IF( ninist == 1 ) THEN !== Output the initial state and forcings ==!
485	CALL dia_wri_state( 'output.init' )
486	ninist = 0
487	ENDIF
488	!
489	IF( nn_write == -1 ) RETURN ! we will never do any output
490	!
491	IF( ln_timing ) CALL timing_start('dia_wri')
492	!
493	! 0. Initialisation
494	! -----------------
495
496	ll_print = .FALSE. ! local variable for debugging
497	ll_print = ll_print .AND. lwp
498
499	! Define frequency of output and means
500	clop = "x" ! no use of the mask value (require less cpu time and otherwise the model crashes)
501	#if defined key_diainstant
502	zsto = nn_write * rdt
503	clop = "inst("//TRIM(clop)//")"
504	#else
505	zsto=rdt
506	clop = "ave("//TRIM(clop)//")"
507	#endif
508	zout = nn_write * rdt
509	zmax = ( nitend - nit000 + 1 ) * rdt
510
511	! Define indices of the horizontal output zoom and vertical limit storage
512	iimi = 1 ; iima = jpi
513	ijmi = 1 ; ijma = jpj
514	ipk = jpk
515
516	! define time axis
517	it = kt
518	itmod = kt - nit000 + 1
519
520
521	! 1. Define NETCDF files and fields at beginning of first time step
522	! -----------------------------------------------------------------
523
524	IF( kt == nit000 ) THEN
525
526	! Define the NETCDF files (one per grid)
527
528	! Compute julian date from starting date of the run
529	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
530	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
531	IF(lwp)WRITE(numout,*)
532	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
533	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
534	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
535	' limit storage in depth = ', ipk
536
537	! WRITE root name in date.file for use by postpro
538	IF(lwp) THEN
539	CALL dia_nam( clhstnam, nn_write,' ' )
540	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
541	WRITE(inum,*) clhstnam
542	CLOSE(inum)
543	ENDIF
544
545	! Define the T grid FILE ( nid_T )
546
547	CALL dia_nam( clhstnam, nn_write, 'grid_T' )
548	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
549	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
550	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
551	& nit000-1, zjulian, rdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
552	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
553	& "m", ipk, gdept_1d, nz_T, "down" )
554	! ! Index of ocean points
555	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
556	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
557	!
558	IF( ln_icebergs ) THEN
559	!
560	!! allocation cant go in dia_wri_alloc because ln_icebergs is only set after
561	!! that routine is called from nemogcm, so do it here immediately before its needed
562	ALLOCATE( ndex_bT(jpijpjnclasses), STAT=ierror )
563	CALL mpp_sum( 'diawri', ierror )
564	IF( ierror /= 0 ) THEN
565	CALL ctl_stop('dia_wri: failed to allocate iceberg diagnostic array')
566	RETURN
567	ENDIF
568	!
569	!! iceberg vertical coordinate is class number
570	CALL histvert( nid_T, "class", "Iceberg class", & ! Vertical grid: class
571	& "number", nclasses, class_num, nb_T )
572	!
573	!! each class just needs the surface index pattern
574	ndim_bT = 3
575	DO jn = 1,nclasses
576	ndex_bT((jn-1)jpijpj+1:jnjpijpj) = ndex_hT(1:jpi*jpj)
577	ENDDO
578	!
579	ENDIF
580
581	! Define the U grid FILE ( nid_U )
582
583	CALL dia_nam( clhstnam, nn_write, 'grid_U' )
584	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
585	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
586	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
587	& nit000-1, zjulian, rdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
588	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
589	& "m", ipk, gdept_1d, nz_U, "down" )
590	! ! Index of ocean points
591	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
592	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
593
594	! Define the V grid FILE ( nid_V )
595
596	CALL dia_nam( clhstnam, nn_write, 'grid_V' ) ! filename
597	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
598	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
599	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
600	& nit000-1, zjulian, rdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
601	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
602	& "m", ipk, gdept_1d, nz_V, "down" )
603	! ! Index of ocean points
604	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
605	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
606
607	! Define the W grid FILE ( nid_W )
608
609	CALL dia_nam( clhstnam, nn_write, 'grid_W' ) ! filename
610	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
611	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
612	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
613	& nit000-1, zjulian, rdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
614	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
615	& "m", ipk, gdepw_1d, nz_W, "down" )
616
617
618	! Declare all the output fields as NETCDF variables
619
620	! !!! nid_T : 3D
621	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
622	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
623	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
624	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
625	IF( .NOT.ln_linssh ) THEN
626	CALL histdef( nid_T, "vovvle3t", "Level thickness" , "m" ,& ! e3t_n
627	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
628	CALL histdef( nid_T, "vovvldep", "T point depth" , "m" ,& ! e3t_n
629	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
630	CALL histdef( nid_T, "vovvldef", "Squared level deformation" , "%^2" ,& ! e3t_n
631	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
632	ENDIF
633	! !!! nid_T : 2D
634	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
635	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
636	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
637	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
638	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
639	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
640	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
641	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
642	CALL histdef( nid_T, "sorunoff", "River runoffs" , "Kg/m2/s", & ! runoffs
643	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
644	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! sfx
645	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
646	IF( ln_linssh ) THEN
647	CALL histdef( nid_T, "sosst_cd", "Concentration/Dilution term on temperature" & ! emp * tsn(:,:,1,jp_tem)
648	& , "KgC/m2/s", & ! sosst_cd
649	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
650	CALL histdef( nid_T, "sosss_cd", "Concentration/Dilution term on salinity" & ! emp * tsn(:,:,1,jp_sal)
651	& , "KgPSU/m2/s",& ! sosss_cd
652	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
653	ENDIF
654	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
655	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
656	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
657	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
658	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
659	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
660	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
661	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
662	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
663	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
664	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
665	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
666	!
667	IF( ln_icebergs ) THEN
668	CALL histdef( nid_T, "calving" , "calving mass input" , "kg/s" , &
669	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
670	CALL histdef( nid_T, "calving_heat" , "calving heat flux" , "XXXX" , &
671	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
672	CALL histdef( nid_T, "berg_floating_melt" , "Melt rate of icebergs + bits" , "kg/m2/s", &
673	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
674	CALL histdef( nid_T, "berg_stored_ice" , "Accumulated ice mass by class" , "kg" , &
675	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
676	IF( ln_bergdia ) THEN
677	CALL histdef( nid_T, "berg_melt" , "Melt rate of icebergs" , "kg/m2/s", &
678	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
679	CALL histdef( nid_T, "berg_buoy_melt" , "Buoyancy component of iceberg melt rate" , "kg/m2/s", &
680	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
681	CALL histdef( nid_T, "berg_eros_melt" , "Erosion component of iceberg melt rate" , "kg/m2/s", &
682	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
683	CALL histdef( nid_T, "berg_conv_melt" , "Convective component of iceberg melt rate", "kg/m2/s", &
684	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
685	CALL histdef( nid_T, "berg_virtual_area" , "Virtual coverage by icebergs" , "m2" , &
686	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
687	CALL histdef( nid_T, "bits_src" , "Mass source of bergy bits" , "kg/m2/s", &
688	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
689	CALL histdef( nid_T, "bits_melt" , "Melt rate of bergy bits" , "kg/m2/s", &
690	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
691	CALL histdef( nid_T, "bits_mass" , "Bergy bit density field" , "kg/m2" , &
692	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
693	CALL histdef( nid_T, "berg_mass" , "Iceberg density field" , "kg/m2" , &
694	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
695	CALL histdef( nid_T, "berg_real_calving" , "Calving into iceberg class" , "kg/s" , &
696	& jpi, jpj, nh_T, nclasses , 1, nclasses , nb_T , 32, clop, zsto, zout )
697	ENDIF
698	ENDIF
699
700	IF( ln_ssr ) THEN
701	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
702	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
703	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
704	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
705	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
706	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
707	ENDIF
708
709	clmx ="l_max(only(x))" ! max index on a period
710	! CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
711	! & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
712	#if defined key_diahth
713	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
714	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
715	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
716	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
717	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
718	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
719	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "J/m2" , & ! htc3
720	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
721	#endif
722
723	CALL histend( nid_T, snc4chunks=snc4set )
724
725	! !!! nid_U : 3D
726	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
727	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
728	IF( ln_wave .AND. ln_sdw) THEN
729	CALL histdef( nid_U, "sdzocrtx", "Stokes Drift Zonal Current" , "m/s" , & ! usd
730	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
731	ENDIF
732	! !!! nid_U : 2D
733	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
734	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
735
736	CALL histend( nid_U, snc4chunks=snc4set )
737
738	! !!! nid_V : 3D
739	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
740	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
741	IF( ln_wave .AND. ln_sdw) THEN
742	CALL histdef( nid_V, "sdmecrty", "Stokes Drift Meridional Current" , "m/s" , & ! vsd
743	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
744	ENDIF
745	! !!! nid_V : 2D
746	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
747	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
748
749	CALL histend( nid_V, snc4chunks=snc4set )
750
751	! !!! nid_W : 3D
752	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
753	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
754	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
755	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
756	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avm
757	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
758
759	IF( ln_zdfddm ) THEN
760	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
761	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
762	ENDIF
763
764	IF( ln_wave .AND. ln_sdw) THEN
765	CALL histdef( nid_W, "sdvecrtz", "Stokes Drift Vertical Current" , "m/s" , & ! wsd
766	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
767	ENDIF
768	! !!! nid_W : 2D
769	CALL histend( nid_W, snc4chunks=snc4set )
770
771	IF(lwp) WRITE(numout,*)
772	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
773	IF(ll_print) CALL FLUSH(numout )
774
775	ENDIF
776
777	! 2. Start writing data
778	! ---------------------
779
780	! ndex(1) est utilise ssi l'avant dernier argument est different de
781	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
782	! donne le nombre d'elements, et ndex la liste des indices a sortir
783
784	IF( lwp .AND. MOD( itmod, nn_write ) == 0 ) THEN
785	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
786	WRITE(numout,*) '~~~~~~ '
787	ENDIF
788
789	IF( .NOT.ln_linssh ) THEN
790	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! heat content
791	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) * e3t_n(:,:,:) , ndim_T , ndex_T ) ! salt content
792	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface heat content
793	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) * e3t_n(:,:,1) , ndim_hT, ndex_hT ) ! sea surface salinity content
794	ELSE
795	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem) , ndim_T , ndex_T ) ! temperature
796	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal) , ndim_T , ndex_T ) ! salinity
797	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem) , ndim_hT, ndex_hT ) ! sea surface temperature
798	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal) , ndim_hT, ndex_hT ) ! sea surface salinity
799	ENDIF
800	IF( .NOT.ln_linssh ) THEN
801	zw3d(:,:,:) = ( ( e3t_n(:,:,:) - e3t_0(:,:,:) ) / e3t_0(:,:,:) * 100 * tmask(:,:,:) ) ** 2
802	CALL histwrite( nid_T, "vovvle3t", it, e3t_n (:,:,:) , ndim_T , ndex_T ) ! level thickness
803	CALL histwrite( nid_T, "vovvldep", it, gdept_n(:,:,:) , ndim_T , ndex_T ) ! t-point depth
804	CALL histwrite( nid_T, "vovvldef", it, zw3d , ndim_T , ndex_T ) ! level thickness deformation
805	ENDIF
806	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
807	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
808	CALL histwrite( nid_T, "sorunoff", it, rnf , ndim_hT, ndex_hT ) ! river runoffs
809	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
810	! (includes virtual salt flux beneath ice
811	! in linear free surface case)
812	IF( ln_linssh ) THEN
813	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_tem)
814	CALL histwrite( nid_T, "sosst_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sst
815	zw2d(:,:) = emp (:,:) * tsn(:,:,1,jp_sal)
816	CALL histwrite( nid_T, "sosss_cd", it, zw2d, ndim_hT, ndex_hT ) ! c/d term on sss
817	ENDIF
818	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
819	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
820	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
821	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
822	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
823	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
824	!
825	IF( ln_icebergs ) THEN
826	!
827	CALL histwrite( nid_T, "calving" , it, berg_grid%calving , ndim_hT, ndex_hT )
828	CALL histwrite( nid_T, "calving_heat" , it, berg_grid%calving_hflx , ndim_hT, ndex_hT )
829	CALL histwrite( nid_T, "berg_floating_melt" , it, berg_grid%floating_melt, ndim_hT, ndex_hT )
830	!
831	CALL histwrite( nid_T, "berg_stored_ice" , it, berg_grid%stored_ice , ndim_bT, ndex_bT )
832	!
833	IF( ln_bergdia ) THEN
834	CALL histwrite( nid_T, "berg_melt" , it, berg_melt , ndim_hT, ndex_hT )
835	CALL histwrite( nid_T, "berg_buoy_melt" , it, buoy_melt , ndim_hT, ndex_hT )
836	CALL histwrite( nid_T, "berg_eros_melt" , it, eros_melt , ndim_hT, ndex_hT )
837	CALL histwrite( nid_T, "berg_conv_melt" , it, conv_melt , ndim_hT, ndex_hT )
838	CALL histwrite( nid_T, "berg_virtual_area" , it, virtual_area , ndim_hT, ndex_hT )
839	CALL histwrite( nid_T, "bits_src" , it, bits_src , ndim_hT, ndex_hT )
840	CALL histwrite( nid_T, "bits_melt" , it, bits_melt , ndim_hT, ndex_hT )
841	CALL histwrite( nid_T, "bits_mass" , it, bits_mass , ndim_hT, ndex_hT )
842	CALL histwrite( nid_T, "berg_mass" , it, berg_mass , ndim_hT, ndex_hT )
843	!
844	CALL histwrite( nid_T, "berg_real_calving" , it, real_calving , ndim_bT, ndex_bT )
845	ENDIF
846	ENDIF
847
848	IF( ln_ssr ) 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	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	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
865	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
866
867	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
868	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
869
870	IF( ln_zad_Aimp ) THEN
871	CALL histwrite( nid_W, "vovecrtz", it, wn + wi , ndim_T, ndex_T ) ! vert. current
872	ELSE
873	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
874	ENDIF
875	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
876	CALL histwrite( nid_W, "votkeavm", it, avm , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
877	IF( ln_zdfddm ) THEN
878	CALL histwrite( nid_W, "voddmavs", it, avs , ndim_T, ndex_T ) ! S vert. eddy diff. coef.
879	ENDIF
880
881	IF( ln_wave .AND. ln_sdw ) THEN
882	CALL histwrite( nid_U, "sdzocrtx", it, usd , ndim_U , ndex_U ) ! i-StokesDrift-current
883	CALL histwrite( nid_V, "sdmecrty", it, vsd , ndim_V , ndex_V ) ! j-StokesDrift-current
884	CALL histwrite( nid_W, "sdvecrtz", it, wsd , ndim_T , ndex_T ) ! StokesDrift vert. current
885	ENDIF
886
887	! 3. Close all files
888	! ---------------------------------------
889	IF( kt == nitend ) THEN
890	CALL histclo( nid_T )
891	CALL histclo( nid_U )
892	CALL histclo( nid_V )
893	CALL histclo( nid_W )
894	ENDIF
895	!
896	IF( ln_timing ) CALL timing_stop('dia_wri')
897	!
898	END SUBROUTINE dia_wri
899	#endif
900
901	SUBROUTINE dia_wri_state( cdfile_name )
902	!!---------------------------------------------------------------------
903	!! * ROUTINE dia_wri_state *
904	!!
905	!! ** Purpose : create a NetCDF file named cdfile_name which contains
906	!! the instantaneous ocean state and forcing fields.
907	!! Used to find errors in the initial state or save the last
908	!! ocean state in case of abnormal end of a simulation
909	!!
910	!! ** Method : NetCDF files using ioipsl
911	!! File 'output.init.nc' is created if ninist = 1 (namelist)
912	!! File 'output.abort.nc' is created in case of abnormal job end
913	!!----------------------------------------------------------------------
914	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
915	!!
916	INTEGER :: inum
917	!!----------------------------------------------------------------------
918	!
919	IF(lwp) WRITE(numout,*)
920	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
921	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
922	IF(lwp) WRITE(numout,*) ' and named :', cdfile_name, '...nc'
923
924	#if defined key_si3
925	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE., kdlev = jpl )
926	#else
927	CALL iom_open( TRIM(cdfile_name), inum, ldwrt = .TRUE. )
928	#endif
929
930	CALL iom_rstput( 0, 0, inum, 'votemper', tsn(:,:,:,jp_tem) ) ! now temperature
931	CALL iom_rstput( 0, 0, inum, 'vosaline', tsn(:,:,:,jp_sal) ) ! now salinity
932	CALL iom_rstput( 0, 0, inum, 'sossheig', sshn ) ! sea surface height
933	CALL iom_rstput( 0, 0, inum, 'vozocrtx', un ) ! now i-velocity
934	CALL iom_rstput( 0, 0, inum, 'vomecrty', vn ) ! now j-velocity
935	IF( ln_zad_Aimp ) THEN
936	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn + wi ) ! now k-velocity
937	ELSE
938	CALL iom_rstput( 0, 0, inum, 'vovecrtz', wn ) ! now k-velocity
939	ENDIF
940	IF( ALLOCATED(ahtu) ) THEN
941	CALL iom_rstput( 0, 0, inum, 'ahtu', ahtu ) ! aht at u-point
942	CALL iom_rstput( 0, 0, inum, 'ahtv', ahtv ) ! aht at v-point
943	ENDIF
944	IF( ALLOCATED(ahmt) ) THEN
945	CALL iom_rstput( 0, 0, inum, 'ahmt', ahmt ) ! ahmt at u-point
946	CALL iom_rstput( 0, 0, inum, 'ahmf', ahmf ) ! ahmf at v-point
947	ENDIF
948	CALL iom_rstput( 0, 0, inum, 'sowaflup', emp - rnf ) ! freshwater budget
949	CALL iom_rstput( 0, 0, inum, 'sohefldo', qsr + qns ) ! total heat flux
950	CALL iom_rstput( 0, 0, inum, 'soshfldo', qsr ) ! solar heat flux
951	CALL iom_rstput( 0, 0, inum, 'soicecov', fr_i ) ! ice fraction
952	CALL iom_rstput( 0, 0, inum, 'sozotaux', utau ) ! i-wind stress
953	CALL iom_rstput( 0, 0, inum, 'sometauy', vtau ) ! j-wind stress
954	IF( .NOT.ln_linssh ) THEN
955	CALL iom_rstput( 0, 0, inum, 'vovvldep', gdept_n ) ! T-cell depth
956	CALL iom_rstput( 0, 0, inum, 'vovvle3t', e3t_n ) ! T-cell thickness
957	END IF
958	IF( ln_wave .AND. ln_sdw ) THEN
959	CALL iom_rstput( 0, 0, inum, 'sdzocrtx', usd ) ! now StokesDrift i-velocity
960	CALL iom_rstput( 0, 0, inum, 'sdmecrty', vsd ) ! now StokesDrift j-velocity
961	CALL iom_rstput( 0, 0, inum, 'sdvecrtz', wsd ) ! now StokesDrift k-velocity
962	ENDIF
963	#if defined key_si3
964	IF( nn_ice == 2 ) THEN ! condition needed in case agrif + ice-model but no-ice in child grid
965	CALL ice_wri_state( inum )
966	ENDIF
967	#endif
968	!
969	CALL iom_close( inum )
970	!
971	END SUBROUTINE dia_wri_state
972
973	!!======================================================================
974	END MODULE diawri

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

Download in other formats: