Context Navigation

diawri.F90 @ 12344

Last change on this file since 12344 was 12340, checked in by acc, 4 years ago

Branch 2019/dev_r11943_MERGE_2019. This commit introduces basic do loop macro
substitution to the 2019 option 1, merge branch. These changes have been SETTE
tested. The only addition is the do_loop_substitute.h90 file in the OCE directory but
the macros defined therein are used throughout the code to replace identifiable, 2D-
and 3D- nested loop opening and closing statements with single-line alternatives. Code
indents are also adjusted accordingly.

The following explanation is taken from comments in the new header file:

This header file contains preprocessor definitions and macros used in the do-loop
substitutions introduced between version 4.0 and 4.2. The primary aim of these macros
is to assist in future applications of tiling to improve performance. This is expected
to be achieved by alternative versions of these macros in selected locations. The
initial introduction of these macros simply replaces all identifiable nested 2D- and
3D-loops with single line statements (and adjusts indenting accordingly). Do loops
are identifiable if they comform to either:

DO jk = ....

DO jj = .... DO jj = ...

DO ji = .... DO ji = ...
. OR .
. .

END DO END DO

END DO END DO

END DO

and white-space variants thereof.

Additionally, only loops with recognised jj and ji loops limits are treated; these are:
Lower limits of 1, 2 or fs_2
Upper limits of jpi, jpim1 or fs_jpim1 (for ji) or jpj, jpjm1 or fs_jpjm1 (for jj)

The macro naming convention takes the form: DO_2D_BT_LR where:

B is the Bottom offset from the PE's inner domain;
T is the Top offset from the PE's inner domain;
L is the Left offset from the PE's inner domain;
R is the Right offset from the PE's inner domain

So, given an inner domain of 2,jpim1 and 2,jpjm1, a typical example would replace:

DO jj = 2, jpj

DO ji = 1, jpim1
.
.

END DO

END DO

with:

DO_2D_01_10
.
.
END_2D

similar conventions apply to the 3D loops macros. jk loop limits are retained
through macro arguments and are not restricted. This includes the possibility of
strides for which an extra set of DO_3DS macros are defined.

In the example definition below the inner PE domain is defined by start indices of
(kIs, kJs) and end indices of (kIe, KJe)

#define DO_2D_00_00 DO jj = kJs, kJe ; DO ji = kIs, kIe
#define END_2D END DO ; END DO

TO DO:

Only conventional nested loops have been identified and replaced by this step. There are constructs such as:

DO jk = 2, jpkm1

z2d(:,:) = z2d(:,:) + e3w(:,:,jk,Kmm) * z3d(:,:,jk) * wmask(:,:,jk)

END DO

which may need to be considered.

Property svn:keywords set to Id

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

Download in other formats: