Context Navigation

diawri.F90 @ 7367

Last change on this file since 7367 was 7367, checked in by deazer, 8 years ago
Contains merged code for CO5 reference.
File size: 65.3 KB

Line
1	MODULE diawri
2	!!======================================================================
3	!! * MODULE diawri *
4	!! Ocean diagnostics : write ocean output files
5	!!=====================================================================
6	!! History : OPA ! 1991-03 (M.-A. Foujols) Original code
7	!! 4.0 ! 1991-11 (G. Madec)
8	!! ! 1992-06 (M. Imbard) correction restart file
9	!! ! 1992-07 (M. Imbard) split into diawri and rstwri
10	!! ! 1993-03 (M. Imbard) suppress writibm
11	!! ! 1998-01 (C. Levy) NETCDF format using ioipsl INTERFACE
12	!! ! 1999-02 (E. Guilyardi) name of netCDF files + variables
13	!! 8.2 ! 2000-06 (M. Imbard) Original code (diabort.F)
14	!! NEMO 1.0 ! 2002-06 (A.Bozec, E. Durand) Original code (diainit.F)
15	!! - ! 2002-09 (G. Madec) F90: Free form and module
16	!! - ! 2002-12 (G. Madec) merge of diabort and diainit, F90
17	!! ! 2005-11 (V. Garnier) Surface pressure gradient organization
18	!! 3.2 ! 2008-11 (B. Lemaire) creation from old diawri
19	!!----------------------------------------------------------------------
20
21	!!----------------------------------------------------------------------
22	!! dia_wri : create the standart output files
23	!! dia_wri_state : create an output NetCDF file for a single instantaeous ocean state and forcing fields
24	!!----------------------------------------------------------------------
25	USE oce ! ocean dynamics and tracers
26	USE dom_oce ! ocean space and time domain
27	USE zdf_oce ! ocean vertical physics
28	USE ldftra_oce ! ocean active tracers: lateral physics
29	USE ldfdyn_oce ! ocean dynamics: lateral physics
30	USE traldf_iso_grif, ONLY : psix_eiv, psiy_eiv
31	USE sol_oce ! solver variables
32	USE sbc_oce ! Surface boundary condition: ocean fields
33	USE sbc_ice ! Surface boundary condition: ice fields
34	USE sbcssr ! restoring term toward SST/SSS climatology
35	USE phycst ! physical constants
36	USE zdfmxl ! mixed layer
37	USE dianam ! build name of file (routine)
38	USE zdfddm ! vertical physics: double diffusion
39	USE diahth ! thermocline diagnostics
40	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
41	USE in_out_manager ! I/O manager
42	USE diadimg ! dimg direct access file format output
43	USE diaar5, ONLY : lk_diaar5
44	USE iom
45	USE ioipsl
46	USE diafoam, ONLY: dia_wri_foam
47	!CEOD USE insitu_tem, ONLY: insitu_t, theta2t
48	USE bartrop_uv, ONLY: un_dm, vn_dm, bartrop_vel
49	#if defined key_lim2
50	USE limwri_2
51	USE ice_2 ! LIM_2 ice model variables
52	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
53	#endif
54	#if defined key_lim3
55	USE ice_3 ! LIM_3 ice model variables
56	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
57	#endif
58	USE daymod ! calendar
59	USE insitu_tem, ONLY: insitu_t, theta2t
60	#if defined key_top
61	USE par_trc ! biogeochemical variables
62	USE trc
63	#endif
64	#if defined key_spm
65	USE spm_con, ONLY: Eps0XS
66	#endif
67	#if defined key_zdftke
68	USE zdftke, ONLY: en
69	#endif
70	USE zdf_oce, ONLY: avt, avm
71	#if defined key_zdfgls
72	USE zdfgls, ONLY: mxln, en
73	#endif
74	USE lib_mpp ! MPP library
75	USE timing ! preformance summary
76	USE wrk_nemo ! working array
77
78	IMPLICIT NONE
79	PRIVATE
80
81	PUBLIC dia_wri_tmb_init ! Called by nemogcm module
82	PUBLIC dia_wri ! routines called by step.F90
83	PUBLIC dia_wri_state
84	PUBLIC dia_wri_alloc ! Called by nemogcm module
85	PUBLIC dia_wri_tide_init ! Called by nemogcm module
86
87	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
88	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
89	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
90	INTEGER :: nid_W, nz_W, nh_W ! grid_W file
91	INTEGER :: ndex(1) ! ???
92	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
93	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_T, ndex_U, ndex_V
94
95	!! * variables for calculating 25-hourly means
96	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: tn_25h , sn_25h , insitu_t_25h
97	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:) :: sshn_25h, hmld_kara_25h
98	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: un_25h , vn_25h , wn_25h
99	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: avt_25h , avm_25h
100	#if defined key_zdfgls \|\| key_zdftke
101	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: en_25h
102	#endif
103	#if defined key_zdfgls
104	REAL(wp),SAVE, ALLOCATABLE, DIMENSION(:,:,:) :: mxln_25h
105	#endif
106	INTEGER, SAVE :: cnt_25h ! Counter for 25 hour means
107
108
109
110	!! * Substitutions
111	# include "zdfddm_substitute.h90"
112	# include "domzgr_substitute.h90"
113	# include "vectopt_loop_substitute.h90"
114	!!----------------------------------------------------------------------
115	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
116	!! $Id$
117	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
118	!!----------------------------------------------------------------------
119	CONTAINS
120
121	INTEGER FUNCTION dia_wri_alloc()
122	!!----------------------------------------------------------------------
123	INTEGER, DIMENSION(2) :: ierr
124	!!----------------------------------------------------------------------
125	!
126	ierr = 0
127	!
128	ALLOCATE( ndex_hT(jpijpj) , ndex_T(jpijpj*jpk) , &
129	& ndex_hU(jpijpj) , ndex_U(jpijpj*jpk) , &
130	& ndex_hV(jpijpj) , ndex_V(jpijpj*jpk) , STAT=ierr(1) )
131	!
132	dia_wri_alloc = MAXVAL(ierr)
133	IF( lk_mpp ) CALL mpp_sum( dia_wri_alloc )
134	!
135	END FUNCTION dia_wri_alloc
136
137	#if defined key_dimgout
138	!!----------------------------------------------------------------------
139	!! 'key_dimgout' DIMG output file
140	!!----------------------------------------------------------------------
141	# include "diawri_dimg.h90"
142
143	#else
144	!!----------------------------------------------------------------------
145	!! Default option NetCDF output file
146	!!----------------------------------------------------------------------
147	# if defined key_iomput
148	!!----------------------------------------------------------------------
149	!! 'key_iomput' use IOM library
150	!!----------------------------------------------------------------------
151
152	SUBROUTINE dia_wri( kt )
153	!!---------------------------------------------------------------------
154	!! * ROUTINE dia_wri *
155	!!
156	!! ** Purpose : Standard output of opa: dynamics and tracer fields
157	!! NETCDF format is used by default
158	!!
159	!! ** Method : use iom_put
160	!!----------------------------------------------------------------------
161	!!
162	INTEGER, INTENT( in ) :: kt ! ocean time-step index
163	!!
164	INTEGER :: ji, jj, jk ! dummy loop indices
165	REAL(wp) :: zztmp, zztmpx, zztmpy !
166	!!
167	REAL(wp), POINTER, DIMENSION(:,:) :: z2d ! 2D workspace
168	REAL(wp), POINTER, DIMENSION(:,:,:) :: z3d ! 3D workspace
169	REAL(wp), POINTER, DIMENSION(:,:,:) :: zwtmb ! temporary workspace for tmb
170	!!----------------------------------------------------------------------
171	!
172	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
173	!
174	CALL wrk_alloc( jpi , jpj , z2d )
175	CALL wrk_alloc( jpi , jpj, jpk , z3d )
176	CALL wrk_alloc( jpi , jpj, 3 , zwtmb )
177	!
178	! Output the initial state and forcings
179	IF( ninist == 1 ) THEN
180	CALL dia_wri_state( 'output.init', kt )
181	ninist = 0
182	ENDIF
183
184	IF (ln_diatide) THEN
185	CALL dia_wri_tide(kt)
186	ENDIF
187
188	CALL iom_put( "toce" , tsn(:,:,:,jp_tem) ) ! temperature
189	CALL theta2t ! in-situ temperature conversion
190	!CEOD CALL iom_put( "tinsitu", insitu_t(:,:,:) ) ! in-situ temperature
191	CALL iom_put( "soce" , tsn(:,:,:,jp_sal) ) ! salinity
192	CALL iom_put( "sst" , tsn(:,:,1,jp_tem) ) ! sea surface temperature
193	CALL iom_put( "sst2" , tsn(:,:,1,jp_tem) * tsn(:,:,1,jp_tem) ) ! square of sea surface temperature
194	CALL iom_put( "sss" , tsn(:,:,1,jp_sal) ) ! sea surface salinity
195	CALL iom_put( "sss2" , tsn(:,:,1,jp_sal) * tsn(:,:,1,jp_sal) ) ! square of sea surface salinity
196	CALL iom_put( "uoce" , un ) ! i-current
197	CALL iom_put( "voce" , vn ) ! j-current
198	CALL iom_put( "ssu" , un(:,:,1) ) ! sea surface U velocity
199	CALL iom_put( "ssv" , vn(:,:,1) ) ! sea surface V velocity
200	IF( cp_cfg == "natl" .OR. cp_cfg == "ind12" ) CALL bartrop_vel ! barotropic velocity conversion
201	!These don't exist independently in this branch so we remove them to get a CO5
202	!that works on the Cray
203	!CEOD CALL iom_put( "uocebt" , un_dm ) ! barotropic i-current
204	!CEOD CALL iom_put( "vocebt" , vn_dm ) ! barotropic j-current
205	CALL iom_put( "avt" , avt ) ! T vert. eddy diff. coef.
206	CALL iom_put( "avm" , avmu ) ! T vert. eddy visc. coef.
207	IF( lk_zdfddm ) THEN
208	CALL iom_put( "avs" , fsavs(:,:,:) ) ! S vert. eddy diff. coef.
209	ENDIF
210	!
211	! If we want tmb values
212
213	IF (ln_diatmb) THEN
214	CALL dia_wri_calctmb( tsn(:,:,:,jp_tem),zwtmb )
215	!ssh already output but here we output it masked
216	CALL iom_put( "sshnmasked" , sshn(:,:)tmask(:,:,1)+missing_val(1-tmask(:,:,1 ) ) ) ! tmb Temperature
217	CALL iom_put( "top_temp" , zwtmb(:,:,1) ) ! tmb Temperature
218	CALL iom_put( "mid_temp" , zwtmb(:,:,2) ) ! tmb Temperature
219	CALL iom_put( "bot_temp" , zwtmb(:,:,3) ) ! tmb Temperature
220	! CALL iom_put( "sotrefml" , hmld_tref(:,:) ) ! "T criterion Mixed Layer Depth
221
222	CALL dia_wri_calctmb( tsn(:,:,:,jp_sal),zwtmb )
223	CALL iom_put( "top_sal" , zwtmb(:,:,1) ) ! tmb Salinity
224	CALL iom_put( "mid_sal" , zwtmb(:,:,2) ) ! tmb Salinity
225	CALL iom_put( "bot_sal" , zwtmb(:,:,3) ) ! tmb Salinity
226
227	CALL dia_wri_calctmb( un(:,:,:),zwtmb )
228	CALL iom_put( "top_u" , zwtmb(:,:,1) ) ! tmb U Velocity
229	CALL iom_put( "mid_u" , zwtmb(:,:,2) ) ! tmb U Velocity
230	CALL iom_put( "bot_u" , zwtmb(:,:,3) ) ! tmb U Velocity
231	!Called in dynspg_ts.F90 CALL iom_put( "baro_u" , un_b ) ! Barotropic U Velocity
232
233	CALL dia_wri_calctmb( vn(:,:,:),zwtmb )
234	CALL iom_put( "top_v" , zwtmb(:,:,1) ) ! tmb V Velocity
235	CALL iom_put( "mid_v" , zwtmb(:,:,2) ) ! tmb V Velocity
236	CALL iom_put( "bot_v" , zwtmb(:,:,3) ) ! tmb V Velocity
237	!Called in dynspg_ts.F90 CALL iom_put( "baro_v" , vn_b ) ! Barotropic V Velocity
238	ENDIF
239
240	DO jj = 2, jpjm1 ! sst gradient
241	DO ji = fs_2, fs_jpim1 ! vector opt.
242	zztmp = tsn(ji,jj,1,jp_tem)
243	zztmpx = ( tsn(ji+1,jj ,1,jp_tem) - zztmp ) / e1u(ji,jj) + ( zztmp - tsn(ji-1,jj ,1,jp_tem) ) / e1u(ji-1,jj )
244	zztmpy = ( tsn(ji ,jj+1,1,jp_tem) - zztmp ) / e2v(ji,jj) + ( zztmp - tsn(ji ,jj-1,1,jp_tem) ) / e2v(ji ,jj-1)
245	z2d(ji,jj) = 0.25 * ( zztmpx * zztmpx + zztmpy * zztmpy ) &
246	& * umask(ji,jj,1) * umask(ji-1,jj,1) * vmask(ji,jj,1) * umask(ji,jj-1,1)
247	END DO
248	END DO
249	CALL lbc_lnk( z2d, 'T', 1. )
250	CALL iom_put( "sstgrad2", z2d ) ! square of module of sst gradient
251	!CDIR NOVERRCHK
252	z2d(:,:) = SQRT( z2d(:,:) )
253	CALL iom_put( "sstgrad" , z2d ) ! module of sst gradient
254
255	IF( lk_diaar5 ) THEN
256	z3d(:,:,jpk) = 0.e0
257	DO jk = 1, jpkm1
258	z3d(:,:,jk) = rau0 * un(:,:,jk) * e2u(:,:) * fse3u(:,:,jk)
259	END DO
260	CALL iom_put( "u_masstr", z3d ) ! mass transport in i-direction
261	zztmp = 0.5 * rcp
262	z2d(:,:) = 0.e0
263	DO jk = 1, jpkm1
264	DO jj = 2, jpjm1
265	DO ji = fs_2, fs_jpim1 ! vector opt.
266	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji+1,jj,jk,jp_tem) )
267	END DO
268	END DO
269	END DO
270	CALL lbc_lnk( z2d, 'U', -1. )
271	CALL iom_put( "u_heattr", z2d ) ! heat transport in i-direction
272	DO jk = 1, jpkm1
273	z3d(:,:,jk) = rau0 * vn(:,:,jk) * e1v(:,:) * fse3v(:,:,jk)
274	END DO
275	CALL iom_put( "v_masstr", z3d ) ! mass transport in j-direction
276	z2d(:,:) = 0.e0
277	DO jk = 1, jpkm1
278	DO jj = 2, jpjm1
279	DO ji = fs_2, fs_jpim1 ! vector opt.
280	z2d(ji,jj) = z2d(ji,jj) + z3d(ji,jj,jk) * zztmp * ( tsn(ji,jj,jk,jp_tem) + tsn(ji,jj+1,jk,jp_tem) )
281	END DO
282	END DO
283	END DO
284	CALL lbc_lnk( z2d, 'V', -1. )
285	CALL iom_put( "v_heattr", z2d ) ! heat transport in i-direction
286	ENDIF
287	!
288	CALL wrk_dealloc( jpi , jpj , z2d )
289	CALL wrk_dealloc( jpi , jpj, jpk , z3d )
290	!
291	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
292	!
293	END SUBROUTINE dia_wri
294
295	#else
296	!!----------------------------------------------------------------------
297	!! Default option use IOIPSL library
298	!!----------------------------------------------------------------------
299
300	SUBROUTINE dia_wri( kt )
301	!!---------------------------------------------------------------------
302	!! * ROUTINE dia_wri *
303	!!
304	!! ** Purpose : Standard output of opa: dynamics and tracer fields
305	!! NETCDF format is used by default
306	!!
307	!! ** Method : At the beginning of the first time step (nit000),
308	!! define all the NETCDF files and fields
309	!! At each time step call histdef to compute the mean if ncessary
310	!! Each nwrite time step, output the instantaneous or mean fields
311	!!----------------------------------------------------------------------
312	!!
313	INTEGER, INTENT( in ) :: kt ! ocean time-step index
314	!!
315	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
316	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
317	INTEGER :: inum = 11 ! temporary logical unit
318	INTEGER :: ji, jj, jk ! dummy loop indices
319	INTEGER :: ierr ! error code return from allocation
320	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
321	REAL(wp) :: zsto, zout, zmax, zjulian, zdt ! local scalars
322	!!
323	REAL(wp), POINTER, DIMENSION(:,:) :: zw2d ! 2D workspace
324	REAL(wp), POINTER, DIMENSION(:,:,:) :: zw3d ! 3D workspace
325	!!----------------------------------------------------------------------
326	!
327	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
328	!
329	CALL wrk_alloc( jpi , jpj , zw2d )
330	IF ( ln_traldf_gdia ) call wrk_alloc( jpi , jpj , jpk , zw3d )
331	!
332	! Output the initial state and forcings
333	IF( ninist == 1 ) THEN
334	CALL dia_wri_state( 'output.init', kt )
335	ninist = 0
336	ENDIF
337	!
338	! -1. Alternative routine
339	!------------------------
340	IF (ln_diafoam) THEN
341	CALL dia_wri_foam( kt )
342	RETURN
343	END IF
344	!
345	! 0. Initialisation
346	! -----------------
347
348	! local variable for debugging
349	ll_print = .FALSE.
350	ll_print = ll_print .AND. lwp
351
352	! Define frequency of output and means
353	zdt = rdt
354	IF( nacc == 1 ) zdt = rdtmin
355	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
356	ELSE ; clop = "x" ! no use of the mask value (require less cpu time)
357	ENDIF
358	#if defined key_diainstant
359	zsto = nwrite * zdt
360	clop = "inst("//TRIM(clop)//")"
361	#else
362	zsto=zdt
363	clop = "ave("//TRIM(clop)//")"
364	#endif
365	zout = nwrite * zdt
366	zmax = ( nitend - nit000 + 1 ) * zdt
367
368	! Define indices of the horizontal output zoom and vertical limit storage
369	iimi = 1 ; iima = jpi
370	ijmi = 1 ; ijma = jpj
371	ipk = jpk
372
373	! define time axis
374	it = kt
375	itmod = kt - nit000 + 1
376
377
378	! 1. Define NETCDF files and fields at beginning of first time step
379	! -----------------------------------------------------------------
380
381	IF( kt == nit000 ) THEN
382
383	! Define the NETCDF files (one per grid)
384
385	! Compute julian date from starting date of the run
386	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
387	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
388	IF(lwp)WRITE(numout,*)
389	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
390	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
391	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
392	' limit storage in depth = ', ipk
393
394	! WRITE root name in date.file for use by postpro
395	IF(lwp) THEN
396	CALL dia_nam( clhstnam, nwrite,' ' )
397	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
398	WRITE(inum,*) clhstnam
399	CLOSE(inum)
400	ENDIF
401
402	! Define the T grid FILE ( nid_T )
403
404	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
405	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
406	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
407	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
408	& nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
409	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
410	& "m", ipk, gdept_0, nz_T, "down" )
411	! ! Index of ocean points
412	CALL wheneq( jpijpjipk, tmask, 1, 1., ndex_T , ndim_T ) ! volume
413	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
414
415	! Define the U grid FILE ( nid_U )
416
417	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
418	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
419	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
420	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
421	& nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
422	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
423	& "m", ipk, gdept_0, nz_U, "down" )
424	! ! Index of ocean points
425	CALL wheneq( jpijpjipk, umask, 1, 1., ndex_U , ndim_U ) ! volume
426	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
427
428	! Define the V grid FILE ( nid_V )
429
430	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
431	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
432	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
433	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
434	& nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
435	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
436	& "m", ipk, gdept_0, nz_V, "down" )
437	! ! Index of ocean points
438	CALL wheneq( jpijpjipk, vmask, 1, 1., ndex_V , ndim_V ) ! volume
439	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
440
441	! Define the W grid FILE ( nid_W )
442
443	CALL dia_nam( clhstnam, nwrite, 'grid_W' ) ! filename
444	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
445	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
446	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
447	& nit000-1, zjulian, zdt, nh_W, nid_W, domain_id=nidom, snc4chunks=snc4set )
448	CALL histvert( nid_W, "depthw", "Vertical W levels", & ! Vertical grid: gdepw
449	& "m", ipk, gdepw_0, nz_W, "down" )
450
451
452	! Declare all the output fields as NETCDF variables
453
454	! !!! nid_T : 3D
455	CALL histdef( nid_T, "votemper", "Temperature" , "C" , & ! tn
456	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
457	CALL histdef( nid_T, "vosaline", "Salinity" , "PSU" , & ! sn
458	& jpi, jpj, nh_T, ipk, 1, ipk, nz_T, 32, clop, zsto, zout )
459	! !!! nid_T : 2D
460	CALL histdef( nid_T, "sosstsst", "Sea Surface temperature" , "C" , & ! sst
461	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
462	CALL histdef( nid_T, "sosaline", "Sea Surface Salinity" , "PSU" , & ! sss
463	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
464	CALL histdef( nid_T, "sossheig", "Sea Surface Height" , "m" , & ! ssh
465	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
466	!!$#if defined key_lim3 \|\| defined key_lim2
467	!!$ ! sowaflup = sowaflep + sorunoff + sowafldp + a term associated to
468	!!$ ! internal damping to Levitus that can be diagnosed from others
469	!!$ ! sowaflcd = sowaflep + sorunoff + sowafldp + iowaflup
470	!!$ CALL histdef( nid_T, "iowaflup", "Ice=>ocean net freshwater" , "kg/m2/s", & ! fsalt
471	!!$ & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
472	!!$ CALL histdef( nid_T, "sowaflep", "atmos=>ocean net freshwater" , "kg/m2/s", & ! fmass
473	!!$ & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
474	!!$#endif
475	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
476	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
477	!!$ CALL histdef( nid_T, "sorunoff", "Runoffs" , "Kg/m2/s", & ! runoffs
478	!!$ & jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
479	CALL histdef( nid_T, "sowaflcd", "concentration/dilution water flux" , "kg/m2/s", & ! (emps-rnf)
480	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
481	CALL histdef( nid_T, "sosalflx", "Surface Salt Flux" , "Kg/m2/s", & ! (emps-rnf) * sn
482	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
483	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
484	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
485	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
486	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
487	CALL histdef( nid_T, "somixhgt", "Turbocline Depth" , "m" , & ! hmld
488	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
489	CALL histdef( nid_T, "somxl010", "Mixed Layer Depth 0.01" , "m" , & ! hmlp
490	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
491	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
492	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
493	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
494	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
495	#if ! defined key_coupled
496	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
497	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
498	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
499	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
500	CALL histdef( nid_T, "sosafldp", "Surface salt flux: damping" , "Kg/m2/s", & ! erp * sn
501	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
502	#endif
503
504
505
506	#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
507	CALL histdef( nid_T, "sohefldp", "Surface Heat Flux: Damping" , "W/m2" , & ! qrp
508	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
509	CALL histdef( nid_T, "sowafldp", "Surface Water Flux: Damping" , "Kg/m2/s", & ! erp
510	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
511	CALL histdef( nid_T, "sosafldp", "Surface salt flux: Damping" , "Kg/m2/s", & ! erp * sn
512	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
513	#endif
514	clmx ="l_max(only(x))" ! max index on a period
515	CALL histdef( nid_T, "sobowlin", "Bowl Index" , "W-point", & ! bowl INDEX
516	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clmx, zsto, zout )
517	#if defined key_diahth
518	CALL histdef( nid_T, "sothedep", "Thermocline Depth" , "m" , & ! hth
519	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
520	CALL histdef( nid_T, "so20chgt", "Depth of 20C isotherm" , "m" , & ! hd20
521	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
522	CALL histdef( nid_T, "so28chgt", "Depth of 28C isotherm" , "m" , & ! hd28
523	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
524	CALL histdef( nid_T, "sohtc300", "Heat content 300 m" , "W" , & ! htc3
525	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
526	#endif
527
528	#if defined key_coupled
529	# if defined key_lim3
530	Must be adapted to LIM3
531	# else
532	CALL histdef( nid_T,"soicetem" , "Ice Surface Temperature" , "K" , & ! tn_ice
533	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
534	CALL histdef( nid_T,"soicealb" , "Ice Albedo" , "[0,1]" , & ! alb_ice
535	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
536	# endif
537	#endif
538
539	CALL histend( nid_T, snc4chunks=snc4set )
540
541	! !!! nid_U : 3D
542	CALL histdef( nid_U, "vozocrtx", "Zonal Current" , "m/s" , & ! un
543	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
544	IF( ln_traldf_gdia ) THEN
545	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
546	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
547	ELSE
548	#if defined key_diaeiv
549	CALL histdef( nid_U, "vozoeivu", "Zonal EIV Current" , "m/s" , & ! u_eiv
550	& jpi, jpj, nh_U, ipk, 1, ipk, nz_U, 32, clop, zsto, zout )
551	#endif
552	END IF
553	! !!! nid_U : 2D
554	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
555	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
556
557	CALL histend( nid_U, snc4chunks=snc4set )
558
559	! !!! nid_V : 3D
560	CALL histdef( nid_V, "vomecrty", "Meridional Current" , "m/s" , & ! vn
561	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
562	IF( ln_traldf_gdia ) THEN
563	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
564	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
565	ELSE
566	#if defined key_diaeiv
567	CALL histdef( nid_V, "vomeeivv", "Meridional EIV Current" , "m/s" , & ! v_eiv
568	& jpi, jpj, nh_V, ipk, 1, ipk, nz_V, 32, clop, zsto, zout )
569	#endif
570	END IF
571	! !!! nid_V : 2D
572	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
573	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
574
575	CALL histend( nid_V, snc4chunks=snc4set )
576
577	! !!! nid_W : 3D
578	CALL histdef( nid_W, "vovecrtz", "Vertical Velocity" , "m/s" , & ! wn
579	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
580	IF( ln_traldf_gdia ) THEN
581	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
582	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
583	ELSE
584	#if defined key_diaeiv
585	CALL histdef( nid_W, "voveeivw", "Vertical EIV Velocity" , "m/s" , & ! w_eiv
586	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
587	#endif
588	END IF
589	CALL histdef( nid_W, "votkeavt", "Vertical Eddy Diffusivity" , "m2/s" , & ! avt
590	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
591	CALL histdef( nid_W, "votkeavm", "Vertical Eddy Viscosity" , "m2/s" , & ! avmu
592	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
593
594	IF( lk_zdfddm ) THEN
595	CALL histdef( nid_W,"voddmavs","Salt Vertical Eddy Diffusivity" , "m2/s" , & ! avs
596	& jpi, jpj, nh_W, ipk, 1, ipk, nz_W, 32, clop, zsto, zout )
597	ENDIF
598	! !!! nid_W : 2D
599	#if defined key_traldf_c2d
600	CALL histdef( nid_W, "soleahtw", "lateral eddy diffusivity" , "m2/s" , & ! ahtw
601	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
602	# if defined key_traldf_eiv
603	CALL histdef( nid_W, "soleaeiw", "eddy induced vel. coeff. at w-point", "m2/s", & ! aeiw
604	& jpi, jpj, nh_W, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
605	# endif
606	#endif
607
608	CALL histend( nid_W, snc4chunks=snc4set )
609
610	IF(lwp) WRITE(numout,*)
611	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
612	IF(ll_print) CALL FLUSH(numout )
613
614	ENDIF
615
616	! 2. Start writing data
617	! ---------------------
618
619	! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
620	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
621	! donne le nombre d'elements, et ndex la liste des indices a sortir
622
623	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
624	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
625	WRITE(numout,*) '~~~~~~ '
626	ENDIF
627
628	! Write fields on T grid
629	CALL histwrite( nid_T, "votemper", it, tsn(:,:,:,jp_tem), ndim_T , ndex_T ) ! temperature
630	CALL histwrite( nid_T, "vosaline", it, tsn(:,:,:,jp_sal), ndim_T , ndex_T ) ! salinity
631	CALL histwrite( nid_T, "sosstsst", it, tsn(:,:,1,jp_tem), ndim_hT, ndex_hT ) ! sea surface temperature
632	CALL histwrite( nid_T, "sosaline", it, tsn(:,:,1,jp_sal), ndim_hT, ndex_hT ) ! sea surface salinity
633	CALL histwrite( nid_T, "sossheig", it, sshn , ndim_hT, ndex_hT ) ! sea surface height
634	!!$#if defined key_lim3 \|\| defined key_lim2
635	!!$ CALL histwrite( nid_T, "iowaflup", it, fsalt(:,:) , ndim_hT, ndex_hT ) ! ice=>ocean water flux
636	!!$ CALL histwrite( nid_T, "sowaflep", it, fmass(:,:) , ndim_hT, ndex_hT ) ! atmos=>ocean water flux
637	!!$#endif
638	CALL histwrite( nid_T, "sowaflup", it, ( emp-rnf ) , ndim_hT, ndex_hT ) ! upward water flux
639	!!$ CALL histwrite( nid_T, "sorunoff", it, runoff , ndim_hT, ndex_hT ) ! runoff
640	CALL histwrite( nid_T, "sowaflcd", it, ( emps-rnf ) , ndim_hT, ndex_hT ) ! c/d water flux
641	zw2d(:,:) = ( emps(:,:) - rnf(:,:) ) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
642	CALL histwrite( nid_T, "sosalflx", it, zw2d , ndim_hT, ndex_hT ) ! c/d salt flux
643	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
644	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
645	CALL histwrite( nid_T, "somixhgt", it, hmld , ndim_hT, ndex_hT ) ! turbocline depth
646	CALL histwrite( nid_T, "somxl010", it, hmlp , ndim_hT, ndex_hT ) ! mixed layer depth
647	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
648	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
649	#if ! defined key_coupled
650	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
651	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
652	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
653	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
654	#endif
655	#if ( defined key_coupled && ! defined key_lim3 && ! defined key_lim2 )
656	CALL histwrite( nid_T, "sohefldp", it, qrp , ndim_hT, ndex_hT ) ! heat flux damping
657	CALL histwrite( nid_T, "sowafldp", it, erp , ndim_hT, ndex_hT ) ! freshwater flux damping
658	IF( ln_ssr ) zw2d(:,:) = erp(:,:) * tsn(:,:,1,jp_sal) * tmask(:,:,1)
659	CALL histwrite( nid_T, "sosafldp", it, zw2d , ndim_hT, ndex_hT ) ! salt flux damping
660	#endif
661	zw2d(:,:) = FLOAT( nmln(:,:) ) * tmask(:,:,1)
662	CALL histwrite( nid_T, "sobowlin", it, zw2d , ndim_hT, ndex_hT ) ! ???
663
664	#if defined key_diahth
665	CALL histwrite( nid_T, "sothedep", it, hth , ndim_hT, ndex_hT ) ! depth of the thermocline
666	CALL histwrite( nid_T, "so20chgt", it, hd20 , ndim_hT, ndex_hT ) ! depth of the 20 isotherm
667	CALL histwrite( nid_T, "so28chgt", it, hd28 , ndim_hT, ndex_hT ) ! depth of the 28 isotherm
668	CALL histwrite( nid_T, "sohtc300", it, htc3 , ndim_hT, ndex_hT ) ! first 300m heaat content
669	#endif
670
671	#if defined key_coupled
672	# if defined key_lim3
673	Must be adapted for LIM3
674	CALL histwrite( nid_T, "soicetem", it, tn_ice , ndim_hT, ndex_hT ) ! surf. ice temperature
675	CALL histwrite( nid_T, "soicealb", it, alb_ice , ndim_hT, ndex_hT ) ! ice albedo
676	# else
677	CALL histwrite( nid_T, "soicetem", it, tn_ice(:,:,1) , ndim_hT, ndex_hT ) ! surf. ice temperature
678	CALL histwrite( nid_T, "soicealb", it, alb_ice(:,:,1), ndim_hT, ndex_hT ) ! ice albedo
679	# endif
680	#endif
681	! Write fields on U grid
682	CALL histwrite( nid_U, "vozocrtx", it, un , ndim_U , ndex_U ) ! i-current
683	IF( ln_traldf_gdia ) THEN
684	IF (.not. ALLOCATED(psix_eiv))THEN
685	ALLOCATE( psix_eiv(jpi,jpj,jpk) , psiy_eiv(jpi,jpj,jpk) , STAT=ierr )
686	IF( lk_mpp ) CALL mpp_sum ( ierr )
687	IF( ierr > 0 ) CALL ctl_stop('STOP', 'diawri: unable to allocate psi{x,y}_eiv')
688	psix_eiv(:,:,:) = 0.0_wp
689	psiy_eiv(:,:,:) = 0.0_wp
690	ENDIF
691	DO jk=1,jpkm1
692	zw3d(:,:,jk) = (psix_eiv(:,:,jk+1) - psix_eiv(:,:,jk))/fse3u(:,:,jk) ! u_eiv = -dpsix/dz
693	END DO
694	zw3d(:,:,jpk) = 0._wp
695	CALL histwrite( nid_U, "vozoeivu", it, zw3d, ndim_U , ndex_U ) ! i-eiv current
696	ELSE
697	#if defined key_diaeiv
698	CALL histwrite( nid_U, "vozoeivu", it, u_eiv, ndim_U , ndex_U ) ! i-eiv current
699	#endif
700	ENDIF
701	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
702
703	! Write fields on V grid
704	CALL histwrite( nid_V, "vomecrty", it, vn , ndim_V , ndex_V ) ! j-current
705	IF( ln_traldf_gdia ) THEN
706	DO jk=1,jpk-1
707	zw3d(:,:,jk) = (psiy_eiv(:,:,jk+1) - psiy_eiv(:,:,jk))/fse3v(:,:,jk) ! v_eiv = -dpsiy/dz
708	END DO
709	zw3d(:,:,jpk) = 0._wp
710	CALL histwrite( nid_V, "vomeeivv", it, zw3d, ndim_V , ndex_V ) ! j-eiv current
711	ELSE
712	#if defined key_diaeiv
713	CALL histwrite( nid_V, "vomeeivv", it, v_eiv, ndim_V , ndex_V ) ! j-eiv current
714	#endif
715	ENDIF
716	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
717
718	! Write fields on W grid
719	CALL histwrite( nid_W, "vovecrtz", it, wn , ndim_T, ndex_T ) ! vert. current
720	IF( ln_traldf_gdia ) THEN
721	DO jk=1,jpk-1
722	DO jj = 2, jpjm1
723	DO ji = fs_2, fs_jpim1 ! vector opt.
724	zw3d(ji,jj,jk) = (psiy_eiv(ji,jj,jk) - psiy_eiv(ji,jj-1,jk))/e2v(ji,jj) + &
725	& (psix_eiv(ji,jj,jk) - psix_eiv(ji-1,jj,jk))/e1u(ji,jj) ! w_eiv = dpsiy/dy + dpsiy/dx
726	END DO
727	END DO
728	END DO
729	zw3d(:,:,jpk) = 0._wp
730	CALL histwrite( nid_W, "voveeivw", it, zw3d , ndim_T, ndex_T ) ! vert. eiv current
731	ELSE
732	# if defined key_diaeiv
733	CALL histwrite( nid_W, "voveeivw", it, w_eiv , ndim_T, ndex_T ) ! vert. eiv current
734	# endif
735	ENDIF
736	CALL histwrite( nid_W, "votkeavt", it, avt , ndim_T, ndex_T ) ! T vert. eddy diff. coef.
737	CALL histwrite( nid_W, "votkeavm", it, avmu , ndim_T, ndex_T ) ! T vert. eddy visc. coef.
738	IF( lk_zdfddm ) THEN
739	CALL histwrite( nid_W, "voddmavs", it, fsavs(:,:,:), ndim_T, ndex_T ) ! S vert. eddy diff. coef.
740	ENDIF
741	#if defined key_traldf_c2d
742	CALL histwrite( nid_W, "soleahtw", it, ahtw , ndim_hT, ndex_hT ) ! lateral eddy diff. coef.
743	# if defined key_traldf_eiv
744	CALL histwrite( nid_W, "soleaeiw", it, aeiw , ndim_hT, ndex_hT ) ! EIV coefficient at w-point
745	# endif
746	#endif
747
748	! 3. Close all files
749	! ---------------------------------------
750	IF( kt == nitend ) THEN
751	CALL histclo( nid_T )
752	CALL histclo( nid_U )
753	CALL histclo( nid_V )
754	CALL histclo( nid_W )
755	ENDIF
756	!
757	CALL wrk_dealloc( jpi , jpj , zw2d )
758	IF ( ln_traldf_gdia ) call wrk_dealloc( jpi , jpj , jpk , zw3d )
759	!
760	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
761	!
762	END SUBROUTINE dia_wri
763	# endif
764
765	#endif
766
767	SUBROUTINE dia_wri_calctmb( infield,outtmb )
768	!!---------------------------------------------------------------------
769	!! * ROUTINE dia_tmb *
770	!!
771	!! ** Purpose : Write diagnostics for Top,Mid, and Bottom of water Column
772	!!
773	!! ** Method :
774	!! use mbathy to find surface, mid and bottom of model levels
775	!!
776	!! History :
777	!! 3.4 ! 04-13 (E. O'Dea) Routine taken from old dia_wri_foam
778	!!----------------------------------------------------------------------
779	!! * Modules used
780
781	! Routine to map 3d field to top, middle, bottom
782	IMPLICIT NONE
783
784	! Routine arguments
785	REAL(wp), DIMENSION(jpi, jpj, jpk), INTENT(IN ) :: infield ! Input 3d field and mask
786	REAL(wp), DIMENSION(jpi, jpj, 3 ), INTENT( OUT) :: outtmb ! Output top, middle, bottom
787
788	! Local variables
789	INTEGER :: ji,jj,jk ! Dummy loop indices
790
791	! Calculate top
792	outtmb(:,:,1) = infield(:,:,1)tmask(:,:,1) + missing_val(1-tmask(:,:,1))
793
794	! Calculate middle
795	DO ji = 1,jpi
796	DO jj = 1,jpj
797	jk = max(1,mbathy(ji,jj)/2)
798	outtmb(ji,jj,2) = infield(ji,jj,jk)tmask(ji,jj,jk) + missing_val(1-tmask(ji,jj,jk))
799	END DO
800	END DO
801
802	! Calculate bottom
803	DO ji = 1,jpi
804	DO jj = 1,jpj
805	jk = max(1,mbathy(ji,jj) - 1)
806	outtmb(ji,jj,3) = infield(ji,jj,jk)tmask(ji,jj,jk) + missing_val(1-tmask(ji,jj,jk))
807	END DO
808	END DO
809
810	END SUBROUTINE dia_wri_calctmb
811
812	SUBROUTINE dia_wri_tmb_init
813	!!---------------------------------------------------------------------------
814	!! * ROUTINE dia_wri_tmb_init *
815	!!
816	!! ** Purpose: Initialization of tmb namelist
817	!!
818	!! ** Method : Read namelist
819	!! History
820	!! 3.4 ! 04-13 (E. O'Dea) Routine to initialize dia_wri_tmb
821	!!---------------------------------------------------------------------------
822	!!
823	INTEGER :: ierror ! local integer
824	!!
825	NAMELIST/nam_diatmb/ ln_diatmb
826	!!
827	!!----------------------------------------------------------------------
828	!
829	REWIND ( numnam ) ! Read Namelist nam_diatmb
830	READ ( numnam, nam_diatmb )
831	!
832	IF(lwp) THEN ! Control print
833	WRITE(numout,*)
834	WRITE(numout,*) 'dia_wri_tmb_init : Output Top, Middle, Bottom Diagnostics'
835	WRITE(numout,*) '~~~~~~~~~~~~'
836	WRITE(numout,*) ' Namelist nam_diatmb : set tmb outputs '
837	WRITE(numout,*) ' Switch for TMB diagnostics (T) or not (F) ln_diatmb = ', ln_diatmb
838	ENDIF
839
840	END SUBROUTINE dia_wri_tmb_init
841
842
843	SUBROUTINE dia_wri_state( cdfile_name, kt )
844	!!---------------------------------------------------------------------
845	!! * ROUTINE dia_wri_state *
846	!!
847	!! ** Purpose : create a NetCDF file named cdfile_name which contains
848	!! the instantaneous ocean state and forcing fields.
849	!! Used to find errors in the initial state or save the last
850	!! ocean state in case of abnormal end of a simulation
851	!!
852	!! ** Method : NetCDF files using ioipsl
853	!! File 'output.init.nc' is created if ninist = 1 (namelist)
854	!! File 'output.abort.nc' is created in case of abnormal job end
855	!!----------------------------------------------------------------------
856	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
857	INTEGER , INTENT( in ) :: kt ! ocean time-step index
858	!!
859	CHARACTER (len=32) :: clname
860	CHARACTER (len=40) :: clop
861	INTEGER :: id_i , nz_i, nh_i
862	INTEGER, DIMENSION(1) :: idex ! local workspace
863	REAL(wp) :: zsto, zout, zmax, zjulian, zdt
864	!!----------------------------------------------------------------------
865	!
866	IF( nn_timing == 1 ) CALL timing_start('dia_wri_state')
867
868	! 0. Initialisation
869	! -----------------
870
871	! Define name, frequency of output and means
872	clname = cdfile_name
873	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
874	zdt = rdt
875	zsto = rdt
876	clop = "inst(x)" ! no use of the mask value (require less cpu time)
877	zout = rdt
878	zmax = ( nitend - nit000 + 1 ) * zdt
879
880	IF(lwp) WRITE(numout,*)
881	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
882	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
883	IF(lwp) WRITE(numout,*) ' and named :', clname, '.nc'
884
885
886	! 1. Define NETCDF files and fields at beginning of first time step
887	! -----------------------------------------------------------------
888
889	! Compute julian date from starting date of the run
890	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis
891	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
892	CALL histbeg( clname, jpi, glamt, jpj, gphit, &
893	1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
894	CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept
895	"m", jpk, gdept_0, nz_i, "down")
896
897	! Declare all the output fields as NetCDF variables
898
899	CALL histdef( id_i, "vosaline", "Salinity" , "PSU" , & ! salinity
900	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
901	CALL histdef( id_i, "votemper", "Temperature" , "C" , & ! temperature
902	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
903	CALL histdef( id_i, "sossheig", "Sea Surface Height" , "m" , & ! ssh
904	& jpi, jpj, nh_i, 1 , 1, 1 , nz_i, 32, clop, zsto, zout )
905	CALL histdef( id_i, "vozocrtx", "Zonal Current" , "m/s" , & ! zonal current
906	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
907	CALL histdef( id_i, "vomecrty", "Meridional Current" , "m/s" , & ! meridonal current
908	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
909	CALL histdef( id_i, "vovecrtz", "Vertical Velocity" , "m/s" , & ! vertical current
910	& jpi, jpj, nh_i, jpk, 1, jpk, nz_i, 32, clop, zsto, zout )
911	CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S", & ! net freshwater
912	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
913	CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2" , & ! net heat flux
914	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
915	CALL histdef( id_i, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! solar flux
916	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
917	CALL histdef( id_i, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
918	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
919	CALL histdef( id_i, "sozotaux", "Zonal Wind Stress" , "N/m2" , & ! i-wind stress
920	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
921	CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress
922	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
923
924	#if defined key_lim2
925	CALL lim_wri_state_2( kt, id_i, nh_i )
926	#else
927	CALL histend( id_i, snc4chunks=snc4set )
928	#endif
929
930	! 2. Start writing data
931	! ---------------------
932	! idex(1) est utilise ssi l'avant dernier argument est diffferent de
933	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
934	! donne le nombre d'elements, et idex la liste des indices a sortir
935	idex(1) = 1 ! init to avoid compil warning
936
937	! Write all fields on T grid
938	CALL histwrite( id_i, "votemper", kt, tsn(:,:,:,jp_tem), jpijpjjpk, idex ) ! now temperature
939	CALL histwrite( id_i, "vosaline", kt, tsn(:,:,:,jp_sal), jpijpjjpk, idex ) ! now salinity
940	CALL histwrite( id_i, "sossheig", kt, sshn , jpi*jpj , idex ) ! sea surface height
941	CALL histwrite( id_i, "vozocrtx", kt, un , jpijpjjpk, idex ) ! now i-velocity
942	CALL histwrite( id_i, "vomecrty", kt, vn , jpijpjjpk, idex ) ! now j-velocity
943	CALL histwrite( id_i, "vovecrtz", kt, wn , jpijpjjpk, idex ) ! now k-velocity
944	CALL histwrite( id_i, "sowaflup", kt, (emp-rnf ) , jpi*jpj , idex ) ! freshwater budget
945	CALL histwrite( id_i, "sohefldo", kt, qsr + qns , jpi*jpj , idex ) ! total heat flux
946	CALL histwrite( id_i, "soshfldo", kt, qsr , jpi*jpj , idex ) ! solar heat flux
947	CALL histwrite( id_i, "soicecov", kt, fr_i , jpi*jpj , idex ) ! ice fraction
948	CALL histwrite( id_i, "sozotaux", kt, utau , jpi*jpj , idex ) ! i-wind stress
949	CALL histwrite( id_i, "sometauy", kt, vtau , jpi*jpj , idex ) ! j-wind stress
950
951	! 3. Close the file
952	! -----------------
953	CALL histclo( id_i )
954	#if ! defined key_iomput && ! defined key_dimgout
955	IF( ninist /= 1 ) THEN
956	CALL histclo( nid_T )
957	CALL histclo( nid_U )
958	CALL histclo( nid_V )
959	CALL histclo( nid_W )
960	ENDIF
961	#endif
962
963	IF( nn_timing == 1 ) CALL timing_stop('dia_wri_state')
964	!
965
966	END SUBROUTINE dia_wri_state
967	!!======================================================================
968	!!======================================================================
969
970	SUBROUTINE dia_wri_tide( kt )
971	!!---------------------------------------------------------------------
972	!! * ROUTINE dia_tide *
973	!!
974	!! ** Purpose : Write diagnostics with M2/S2 tide removed
975	!!
976	!! ** Method :
977	!! 25hr mean outputs for shelf seas
978	!!
979	!! History :
980	!! ?.0 ! 07-04 (A. Hines) New routine, developed from dia_wri_foam
981	!! 3.4 ! 02-13 (J. Siddorn) Routine taken from old dia_wri_foam
982	!!----------------------------------------------------------------------
983	!! * Modules used
984
985	IMPLICIT NONE
986
987	!! * Arguments
988	INTEGER, INTENT( in ) :: kt ! ocean time-step index
989
990
991	!! * Local declarations
992	INTEGER :: ji, jj, jk
993
994	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
995	REAL(wp) :: zsto, zout, zmax, zjulian, zdt, zmdi ! temporary integers
996	INTEGER :: i_steps ! no of timesteps per hour
997	REAL(wp), DIMENSION(jpi,jpj ) :: zw2d, un_dm, vn_dm ! temporary workspace
998	REAL(wp), DIMENSION(jpi,jpj,jpk) :: zw3d ! temporary workspace
999	REAL(wp), DIMENSION(jpi,jpj,3) :: zwtmb ! temporary workspace
1000	INTEGER :: nyear0, nmonth0,nday0 ! start year,month,day
1001	!#if defined key_top
1002	! CHARACTER (len=20) :: cltra, cltrau
1003	! CHARACTER (len=80) :: cltral
1004	! INTEGER :: jn, jl
1005	!#endif
1006	!#if defined key_spm
1007	! ! variables needed to calculate visibility field from sediment fields
1008	! REAL(wp), DIMENSION(jpi,jpj,jpk) :: vis3d ! derived 3D visibility field
1009	! REAL(wp) :: epsessX = 0.07d-03 ! attenuation coefficient applied to the sediment (as used in ERSEM)
1010	! REAL(wp) :: tiny = 1.0d-15 ! to prevent division by zero in visibility calculation
1011	!#endif
1012
1013	!!----------------------------------------------------------------------
1014
1015	! 0. Initialisation
1016	! -----------------
1017	! Define frequency of summing to create 25 h mean
1018	zdt = rdt
1019	IF( nacc == 1 ) zdt = rdtmin
1020
1021	IF( MOD( 3600,INT(zdt) ) == 0 ) THEN
1022	i_steps = 3600/INT(zdt)
1023	ELSE
1024	CALL ctl_stop('STOP', 'dia_wri_tide: timestep must give MOD(3600,rdt) = 0 otherwise no hourly values are possible')
1025	ENDIF
1026
1027	#if defined key_lim3 \|\| defined key_lim2
1028	CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
1029	#endif
1030	#if defined key_spm \|\| defined key_MOersem
1031	CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ERSEM')
1032	#endif
1033
1034	! local variable for debugging
1035	ll_print = ll_print .AND. lwp
1036
1037	! Sum of 25 hourly instantaneous values to give a 25h mean from 24hours
1038	! every day
1039	IF( MOD( kt, i_steps ) == 0 .and. kt .ne. nn_it000 ) THEN
1040
1041	IF (lwp) THEN
1042	WRITE(numout,*) 'dia_wri_tide : Summing instantaneous hourly diagnostics at timestep ',kt
1043	WRITE(numout,*) '~~~~~~~~~~~~ '
1044	ENDIF
1045
1046	tn_25h(:,:,:) = tn_25h(:,:,:) + tsn(:,:,:,jp_tem)
1047	sn_25h(:,:,:) = sn_25h(:,:,:) + tsn(:,:,:,jp_sal)
1048	CALL theta2t
1049	insitu_t_25h(:,:,:) = insitu_t_25h(:,:,:) + insitu_t(:,:,:)
1050	sshn_25h(:,:) = sshn_25h(:,:) + sshn (:,:)
1051	! hmld_kara_25h(:,:) = hmld_kara_25h(:,:) + hmld_kara(:,:)
1052	#if defined key_lim3 \|\| defined key_lim2
1053	hsnif_25h(:,:) = hsnif_25h(:,:) + hsnif(:,:)
1054	hicif_25h(:,:) = hicif_25h(:,:) + hicif(:,:)
1055	frld_25h(:,:) = frld_25h(:,:) + frld(:,:)
1056	#endif
1057	#if defined key_spm \|\| defined key_MOersem
1058	trn_25h(:,:,:,:) = trn_25h(:,:,:,:) + trn (:,:,:,:)
1059	trc3d_25h(:,:,:,:) = trc3d_25h(:,:,:,:) + trc3d(:,:,:,:)
1060	trc2d_25h(:,:,:) = trc2d_25h(:,:,:) + trc2d(:,:,:)
1061	#endif
1062	un_25h(:,:,:) = un_25h(:,:,:) + un(:,:,:)
1063	vn_25h(:,:,:) = vn_25h(:,:,:) + vn(:,:,:)
1064	wn_25h(:,:,:) = wn_25h(:,:,:) + wn(:,:,:)
1065	avt_25h(:,:,:) = avt_25h(:,:,:) + avt(:,:,:)
1066	avm_25h(:,:,:) = avm_25h(:,:,:) + avm(:,:,:)
1067	# if defined key_zdfgls \|\| defined key_zdftke
1068	en_25h(:,:,:) = en_25h(:,:,:) + en(:,:,:)
1069	#endif
1070	# if defined key_zdfgls
1071	mxln_25h(:,:,:) = mxln_25h(:,:,:) + mxln(:,:,:)
1072	#endif
1073	cnt_25h = cnt_25h + 1
1074
1075	IF (lwp) THEN
1076	WRITE(numout,*) 'dia_tide : Summed the following number of hourly values so far',cnt_25h
1077	ENDIF
1078
1079	ENDIF ! MOD( kt, i_steps ) == 0
1080
1081	! Write data for 25 hour mean output streams
1082	IF( cnt_25h .EQ. 25 .AND. MOD( kt, i_steps*24) == 0 .AND. kt .NE. nn_it000 ) THEN
1083
1084	IF(lwp) THEN
1085	WRITE(numout,*) 'dia_wri_tide : Writing 25 hour mean tide diagnostics at timestep', kt
1086	WRITE(numout,*) '~~~~~~~~~~~~ '
1087	ENDIF
1088
1089	tn_25h(:,:,:) = tn_25h(:,:,:) / 25.0_wp
1090	sn_25h(:,:,:) = sn_25h(:,:,:) / 25.0_wp
1091	insitu_t_25h(:,:,:) = insitu_t_25h(:,:,:) / 25.0_wp
1092	sshn_25h(:,:) = sshn_25h(:,:) / 25.0_wp
1093	! hmld_kara_25h(:,:) = hmld_kara_25h(:,:) / 25.0_wp
1094	#if defined key_lim3 \|\| defined key_lim2
1095	hsnif_25h(:,:) = hsnif_25h(:,:) / 25.0_wp
1096	hicif_25h(:,:) = hicif_25h(:,:) / 25.0_wp
1097	frld_25h(:,:) = frld_25h(:,:) / 25.0_wp
1098	#endif
1099	#if defined key_spm \|\| defined key_MOersem
1100	trn_25h(:,:,:,:) = trn_25h(:,:,:,:) / 25.0_wp
1101	trc3d_25h(:,:,:,:) = trc3d_25h(:,:,:,:) / 25.0_wp
1102	trc2d_25h(:,:,:) = trc2d_25h(:,:,:) / 25.0_wp
1103	#endif
1104	un_25h(:,:,:) = un_25h(:,:,:) / 25.0_wp
1105	vn_25h(:,:,:) = vn_25h(:,:,:) / 25.0_wp
1106	wn_25h(:,:,:) = wn_25h(:,:,:) / 25.0_wp
1107	avt_25h(:,:,:) = avt_25h(:,:,:) / 25.0_wp
1108	avm_25h(:,:,:) = avm_25h(:,:,:) / 25.0_wp
1109	# if defined key_zdfgls \|\| defined key_zdftke
1110	en_25h(:,:,:) = en_25h(:,:,:) / 25.0_wp
1111	#endif
1112	# if defined key_zdfgls
1113	mxln_25h(:,:,:) = mxln_25h(:,:,:) / 25.0_wp
1114	#endif
1115
1116	IF (lwp) WRITE(numout,*) 'dia_wri_tide : Mean calculated by dividing 25 hour sums and writing output'
1117
1118	zmdi=missing_val ! for masking
1119	! write tracers (instantaneous)
1120	zw3d(:,:,:) = tn_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1121	CALL iom_put("temper25h", zw3d) ! potential temperature
1122	CALL theta2t ! calculate insitu temp
1123	zw3d(:,:,:) = insitu_t_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1124	CALL iom_put("tempis25h", zw3d) ! in-situ temperature
1125	zw3d(:,:,:) = sn_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1126	CALL iom_put( "salin25h", zw3d ) ! salinity
1127	zw2d(:,:) = sshn_25h(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
1128	CALL iom_put( "ssh25h", zw2d ) ! sea surface
1129	! zw2d(:,:) = hmld_kara_25h(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
1130	! CALL iom_put( "kara25h", zw2d ) ! mixed layer
1131
1132	#if defined key_lim3 \|\| defined key_lim2
1133	! Write ice model variables (instantaneous)
1134	zw2d(:,:) = hsnif_25h(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
1135	CALL iom_put("isnowthi", zw2d ) ! ice thickness
1136	zw2d(:,:) = hicif_25h(:,:)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
1137	CALL iom_put("iicethic", zw2d ) ! ice thickness
1138	zw2d(:,:) = (1.0-frld_25h(:,:))tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
1139	CALL iom_put("iiceconc", zw2d ) ! ice concetration
1140	#endif
1141	#if defined key_spm \|\| key_MOersem
1142	! output biogeochemical variables:
1143	! output main tracers
1144	DO jn = 1, jptra
1145	cltra = ctrcnm(jn) ! short title for tracer
1146	zw3d(:,:,:) = trn_25h(:,:,:,jn)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1147	IF( lutsav(jn) ) CALL iom_put( cltra, zw3d ) ! temperature
1148	END DO
1149	! more 3D horizontal arrays from diagnostics
1150	DO jl = 1, jpdia3d
1151	cltra = ctrc3d(jl) ! short title for 3D diagnostic
1152	zw3d(:,:,:) = trc3d_25h(:,:,:,jl)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1153	CALL iom_put( cltra, zw3d )
1154	END DO
1155	! more 2D horizontal arrays from diagnostics
1156	DO jl = 1, jpdia2d
1157	cltra = ctrc2d(jl) ! short title for 2D diagnostic
1158	zw2d(:,:) = trc2d_25h(:,:,jl)tmask(:,:,1) + zmdi(1.0-tmask(:,:,1))
1159	CALL iom_put(cltra, zw2d )
1160	END DO
1161	#endif
1162
1163	! Write velocities (instantaneous)
1164	zw3d(:,:,:) = un_25h(:,:,:)umask(:,:,:) + zmdi(1.0-umask(:,:,:))
1165	CALL iom_put("vozocrtx25h", zw3d) ! i-current
1166	zw3d(:,:,:) = vn_25h(:,:,:)vmask(:,:,:) + zmdi(1.0-vmask(:,:,:))
1167	CALL iom_put("vomecrty25h", zw3d ) ! j-current
1168
1169	zw3d(:,:,:) = wn_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1170	CALL iom_put("vomecrtz25h", zw3d ) ! k-current
1171	zw3d(:,:,:) = avt_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1172	CALL iom_put("avt25h", zw3d ) ! diffusivity
1173	zw3d(:,:,:) = avm_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1174	CALL iom_put("avm25h", zw3d) ! viscosity
1175	#if defined key_zdftke \|\| defined key_zdfgls
1176	zw3d(:,:,:) = en_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1177	CALL iom_put("tke25h", zw3d) ! tke
1178	#endif
1179	#if defined key_zdfgls
1180	zw3d(:,:,:) = mxln_25h(:,:,:)tmask(:,:,:) + zmdi(1.0-tmask(:,:,:))
1181	CALL iom_put( "mxln25h",zw3d)
1182	#endif
1183
1184	! After the write reset the values to cnt=1 and sum values equal current value
1185	tn_25h(:,:,:) = tsn(:,:,:,jp_tem)
1186	sn_25h(:,:,:) = tsn(:,:,:,jp_sal)
1187	CALL theta2t
1188	insitu_t_25h(:,:,:) = insitu_t(:,:,:)
1189	sshn_25h(:,:) = sshn (:,:)
1190	! hmld_kara_25h(:,:) = hmld_kara(:,:)
1191	#if defined key_lim3 \|\| defined key_lim2
1192	hsnif_25h(:,:) = hsnif(:,:)
1193	hicif_25h(:,:) = hicif(:,:)
1194	frld_25h(:,:) = frld(:,:)
1195	#endif
1196	#if defined key_spm \|\| defined key_MOersem
1197	trn_25h(:,:,:,:) = trn (:,:,:,:)
1198	trc3d_25h(:,:,:,:) = trc3d(:,:,:,:)
1199	trc2d_25h(:,:,:) = trc2d(:,:,:)
1200	#endif
1201	un_25h(:,:,:) = un(:,:,:)
1202	vn_25h(:,:,:) = vn(:,:,:)
1203	wn_25h(:,:,:) = wn(:,:,:)
1204	avt_25h(:,:,:) = avt(:,:,:)
1205	avm_25h(:,:,:) = avm(:,:,:)
1206	# if defined key_zdfgls \|\| defined key_zdftke
1207	en_25h(:,:,:) = en(:,:,:)
1208	#endif
1209	# if defined key_zdfgls
1210	mxln_25h(:,:,:) = mxln(:,:,:)
1211	#endif
1212	cnt_25h = 1
1213	IF (lwp) WRITE(numout,*) 'dia_wri_tide : After 25hr mean write, reset sum to current value and cnt_25h to one for overlapping average',cnt_25h
1214
1215	ENDIF ! cnt_25h .EQ. 25 .AND. MOD( kt, i_steps * 24) == 0 .AND. kt .NE. nn_it000
1216
1217	END SUBROUTINE dia_wri_tide
1218	!!======================================================================
1219
1220	SUBROUTINE dia_wri_tide_init
1221	!!---------------------------------------------------------------------------
1222	!! * ROUTINE dia_wri_tide_init *
1223	!!
1224	!! ** Purpose: Initialization of 25hour mean variables for detided output
1225	!!
1226	!! ** Method : Read namelist, allocate and assign initial values
1227	!! History
1228	!! 3.4 ! 03-13 (E. O'Dea) Routine to initialize dia_wri_tide
1229	!!---------------------------------------------------------------------------
1230	!!
1231	INTEGER :: ierror ! local integer
1232	!!
1233	NAMELIST/nam_diatide/ ln_diatide
1234	!!
1235	!!----------------------------------------------------------------------
1236	!
1237	REWIND ( numnam ) ! Read Namelist nam_tiatide
1238	READ ( numnam, nam_diatide )
1239	!
1240	IF(lwp) THEN ! Control print
1241	WRITE(numout,*)
1242	WRITE(numout,*) 'dia_wri_tide_init : Output 25 hour Mean Diagnostics'
1243	WRITE(numout,*) '~~~~~~~~~~~~'
1244	WRITE(numout,*) ' Namelist nam_diatide : set 25hour mean outputs '
1245	WRITE(numout,*) ' Switch for 25 hour mean diagnostics (T) or not (F) ln_diatide = ', ln_diatide
1246	ENDIF
1247	IF( .NOT. ln_diatide ) RETURN
1248
1249	! ------------------- !
1250	! 1 - Allocate memory !
1251	! ------------------- !
1252	ALLOCATE( tn_25h(jpi,jpj,jpk), STAT=ierror )
1253	IF( ierror > 0 ) THEN
1254	CALL ctl_stop( 'dia_tide: unable to allocate tn_25h' ) ; RETURN
1255	ENDIF
1256	ALLOCATE( sn_25h(jpi,jpj,jpk), STAT=ierror )
1257	IF( ierror > 0 ) THEN
1258	CALL ctl_stop( 'dia_tide: unable to allocate sn_25h' ) ; RETURN
1259	ENDIF
1260	ALLOCATE( insitu_t_25h(jpi,jpj,jpk), STAT=ierror )
1261	IF( ierror > 0 ) THEN
1262	CALL ctl_stop( 'dia_tide: unable to allocate insitu_t_25h' ) ; RETURN
1263	ENDIF
1264	ALLOCATE( un_25h(jpi,jpj,jpk), STAT=ierror )
1265	IF( ierror > 0 ) THEN
1266	CALL ctl_stop( 'dia_tide: unable to allocate un_25h' ) ; RETURN
1267	ENDIF
1268	ALLOCATE( vn_25h(jpi,jpj,jpk), STAT=ierror )
1269	IF( ierror > 0 ) THEN
1270	CALL ctl_stop( 'dia_tide: unable to allocate vn_25h' ) ; RETURN
1271	ENDIF
1272	ALLOCATE( wn_25h(jpi,jpj,jpk), STAT=ierror )
1273	IF( ierror > 0 ) THEN
1274	CALL ctl_stop( 'dia_tide: unable to allocate wn_25h' ) ; RETURN
1275	ENDIF
1276	ALLOCATE( avt_25h(jpi,jpj,jpk), STAT=ierror )
1277	IF( ierror > 0 ) THEN
1278	CALL ctl_stop( 'dia_tide: unable to allocate avt_25h' ) ; RETURN
1279	ENDIF
1280	ALLOCATE( avm_25h(jpi,jpj,jpk), STAT=ierror )
1281	IF( ierror > 0 ) THEN
1282	CALL ctl_stop( 'dia_tide: unable to allocate avm_25h' ) ; RETURN
1283	ENDIF
1284	# if defined key_zdfgls \|\| defined key_zdftke
1285	ALLOCATE( en_25h(jpi,jpj,jpk), STAT=ierror )
1286	IF( ierror > 0 ) THEN
1287	CALL ctl_stop( 'dia_tide: unable to allocate en_25h' ) ; RETURN
1288	ENDIF
1289	#endif
1290	# if defined key_zdfgls
1291	ALLOCATE( mxln_25h(jpi,jpj,jpk), STAT=ierror )
1292	IF( ierror > 0 ) THEN
1293	CALL ctl_stop( 'dia_tide: unable to allocate mxln_25h' ) ; RETURN
1294	ENDIF
1295	#endif
1296	ALLOCATE( sshn_25h(jpi,jpj), STAT=ierror )
1297	IF( ierror > 0 ) THEN
1298	CALL ctl_stop( 'dia_tide: unable to allocate sshn_25h' ) ; RETURN
1299	ENDIF
1300	ALLOCATE( hmld_kara_25h(jpi,jpj), STAT=ierror )
1301	IF( ierror > 0 ) THEN
1302	CALL ctl_stop( 'dia_tide: unable to allocate hmld_kara_25h' ) ; RETURN
1303	ENDIF
1304	! ------------------------- !
1305	! 2 - Assign Initial Values !
1306	! ------------------------- !
1307	cnt_25h = 1 ! sets the first value of sum at timestep 1 (note - should strictly be at timestep zero so before values used where possible)
1308	tn_25h(:,:,:) = tsb(:,:,:,jp_tem)
1309	sn_25h(:,:,:) = tsb(:,:,:,jp_sal)
1310	CALL theta2t
1311	insitu_t_25h(:,:,:) = insitu_t(:,:,:)
1312	sshn_25h(:,:) = sshb(:,:)
1313	! hmld_kara_25h(:,:) = hmld_kara(:,:)
1314	un_25h(:,:,:) = ub(:,:,:)
1315	vn_25h(:,:,:) = vb(:,:,:)
1316	wn_25h(:,:,:) = wn(:,:,:)
1317	avt_25h(:,:,:) = avt(:,:,:)
1318	avm_25h(:,:,:) = avm(:,:,:)
1319	# if defined key_zdfgls \|\| defined key_zdftke
1320	en_25h(:,:,:) = en(:,:,:)
1321	#endif
1322	# if defined key_zdfgls
1323	mxln_25h(:,:,:) = mxln(:,:,:)
1324	#endif
1325	#if defined key_lim3 \|\| defined key_lim2
1326	CALL ctl_stop('STOP', 'dia_wri_tide not setup yet to do tidemean ice')
1327	#endif
1328
1329	! -------------------------- !
1330	! 3 - Return to dia_wri_tide !
1331	! -------------------------- !
1332
1333
1334	END SUBROUTINE dia_wri_tide_init
1335
1336
1337	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/CO5_package_branch/NEMOGCM/NEMO/OPA_SRC/DIA/diawri.F90 @ 7367

Download in other formats: