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

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

diawri.F90 in branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/SAS_SRC – NEMO

Context Navigation

source: branches/2014/dev_r4621_NOC4_BDY_VERT_INTERP/NEMOGCM/NEMO/SAS_SRC/diawri.F90 @ 5901

Last change on this file since 5901 was 5901, checked in by jamesharle, 8 years ago
merging branch with head of the trunk
Property svn:keywords set to `Id`
File size: 22.0 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 sol_oce ! solver variables
29	USE sbc_oce ! Surface boundary condition: ocean fields
30	USE sbc_ice ! Surface boundary condition: ice fields
31	USE sbcssr ! restoring term toward SST/SSS climatology
32	USE phycst ! physical constants
33	USE zdfmxl ! mixed layer
34	USE dianam ! build name of file (routine)
35	USE zdfddm ! vertical physics: double diffusion
36	USE diahth ! thermocline diagnostics
37	USE lbclnk ! ocean lateral boundary conditions (or mpp link)
38	USE in_out_manager ! I/O manager
39	USE diadimg ! dimg direct access file format output
40	USE diaar5, ONLY : lk_diaar5
41	USE iom
42	USE ioipsl
43	#if defined key_lim2
44	USE limwri_2
45	#endif
46	USE lib_mpp ! MPP library
47	USE timing ! preformance summary
48	USE wrk_nemo ! working array
49
50	IMPLICIT NONE
51	PRIVATE
52
53	PUBLIC dia_wri ! routines called by step.F90
54	PUBLIC dia_wri_state
55	PUBLIC dia_wri_alloc ! Called by nemogcm module
56
57	INTEGER :: nid_T, nz_T, nh_T, ndim_T, ndim_hT ! grid_T file
58	INTEGER :: nid_U, nz_U, nh_U, ndim_U, ndim_hU ! grid_U file
59	INTEGER :: nid_V, nz_V, nh_V, ndim_V, ndim_hV ! grid_V file
60	INTEGER :: ndex(1) ! ???
61	INTEGER, SAVE, ALLOCATABLE, DIMENSION(:) :: ndex_hT, ndex_hU, ndex_hV
62
63	!! * Substitutions
64	# include "zdfddm_substitute.h90"
65	# include "domzgr_substitute.h90"
66	# include "vectopt_loop_substitute.h90"
67	!!----------------------------------------------------------------------
68	!! NEMO/OPA 3.3 , NEMO Consortium (2010)
69	!! $Id$
70	!! Software governed by the CeCILL licence (NEMOGCM/NEMO_CeCILL.txt)
71	!!----------------------------------------------------------------------
72	CONTAINS
73
74	INTEGER FUNCTION dia_wri_alloc()
75	!!----------------------------------------------------------------------
76	INTEGER :: ierr
77	!!----------------------------------------------------------------------
78	!
79	ALLOCATE( ndex_hT(jpijpj), ndex_hU(jpijpj), ndex_hV(jpi*jpj), STAT=dia_wri_alloc )
80	IF( lk_mpp ) CALL mpp_sum( dia_wri_alloc )
81	!
82	END FUNCTION dia_wri_alloc
83
84	#if defined key_dimgout
85	!!----------------------------------------------------------------------
86	!! 'key_dimgout' DIMG output file
87	!!----------------------------------------------------------------------
88	# include "diawri_dimg.h90"
89
90	#else
91	!!----------------------------------------------------------------------
92	!! Default option NetCDF output file
93	!!----------------------------------------------------------------------
94	# if defined key_iomput
95	!!----------------------------------------------------------------------
96	!! 'key_iomput' use IOM library
97	!!----------------------------------------------------------------------
98
99	SUBROUTINE dia_wri( kt )
100	!!---------------------------------------------------------------------
101	!! * ROUTINE dia_wri *
102	!!
103	!! ** Purpose : Standard output of opa: dynamics and tracer fields
104	!! NETCDF format is used by default
105	!! Standalone surface scheme
106	!!
107	!! ** Method : use iom_put
108	!!----------------------------------------------------------------------
109	!!
110	INTEGER, INTENT( in ) :: kt ! ocean time-step index
111	!!----------------------------------------------------------------------
112	!
113	!! no relevant 2D arrays to write in iomput case
114	!
115	END SUBROUTINE dia_wri
116
117	#else
118	!!----------------------------------------------------------------------
119	!! Default option use IOIPSL library
120	!!----------------------------------------------------------------------
121
122	SUBROUTINE dia_wri( kt )
123	!!---------------------------------------------------------------------
124	!! * ROUTINE dia_wri *
125	!!
126	!! ** Purpose : Standard output of opa: dynamics and tracer fields
127	!! NETCDF format is used by default
128	!!
129	!! ** Method : At the beginning of the first time step (nit000),
130	!! define all the NETCDF files and fields
131	!! At each time step call histdef to compute the mean if ncessary
132	!! Each nwrite time step, output the instantaneous or mean fields
133	!!----------------------------------------------------------------------
134	!!
135	INTEGER, INTENT( in ) :: kt ! ocean time-step index
136	!!
137	LOGICAL :: ll_print = .FALSE. ! =T print and flush numout
138	CHARACTER (len=40) :: clhstnam, clop, clmx ! local names
139	INTEGER :: inum = 11 ! temporary logical unit
140	INTEGER :: ji, jj, jk ! dummy loop indices
141	INTEGER :: ierr ! error code return from allocation
142	INTEGER :: iimi, iima, ipk, it, itmod, ijmi, ijma ! local integers
143	REAL(wp) :: zsto, zout, zmax, zjulian, zdt ! local scalars
144	!!----------------------------------------------------------------------
145	!
146	IF( nn_timing == 1 ) CALL timing_start('dia_wri')
147	!
148	! Output the initial state and forcings
149	IF( ninist == 1 ) THEN
150	CALL dia_wri_state( 'output.init', kt )
151	ninist = 0
152	ENDIF
153	!
154	! 0. Initialisation
155	! -----------------
156
157	! local variable for debugging
158	ll_print = .FALSE.
159	ll_print = ll_print .AND. lwp
160
161	! Define frequency of output and means
162	zdt = rdt
163	IF( nacc == 1 ) zdt = rdtmin
164	IF( ln_mskland ) THEN ; clop = "only(x)" ! put 1.e+20 on land (very expensive!!)
165	ELSE ; clop = "x" ! no use of the mask value (require less cpu time)
166	ENDIF
167	#if defined key_diainstant
168	zsto = nwrite * zdt
169	clop = "inst("//TRIM(clop)//")"
170	#else
171	zsto=zdt
172	clop = "ave("//TRIM(clop)//")"
173	#endif
174	zout = nwrite * zdt
175	zmax = ( nitend - nit000 + 1 ) * zdt
176
177	! Define indices of the horizontal output zoom and vertical limit storage
178	iimi = 1 ; iima = jpi
179	ijmi = 1 ; ijma = jpj
180	ipk = jpk
181
182	! define time axis
183	it = kt
184	itmod = kt - nit000 + 1
185
186
187	! 1. Define NETCDF files and fields at beginning of first time step
188	! -----------------------------------------------------------------
189
190	IF( kt == nit000 ) THEN
191
192	! Define the NETCDF files (one per grid)
193
194	! Compute julian date from starting date of the run
195	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian )
196	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
197	IF(lwp)WRITE(numout,*)
198	IF(lwp)WRITE(numout,*) 'Date 0 used :', nit000, ' YEAR ', nyear, &
199	& ' MONTH ', nmonth, ' DAY ', nday, 'Julian day : ', zjulian
200	IF(lwp)WRITE(numout,*) ' indexes of zoom = ', iimi, iima, ijmi, ijma, &
201	' limit storage in depth = ', ipk
202
203	! WRITE root name in date.file for use by postpro
204	IF(lwp) THEN
205	CALL dia_nam( clhstnam, nwrite,' ' )
206	CALL ctl_opn( inum, 'date.file', 'REPLACE', 'FORMATTED', 'SEQUENTIAL', -1, numout, lwp, narea )
207	WRITE(inum,*) clhstnam
208	CLOSE(inum)
209	ENDIF
210
211	! Define the T grid FILE ( nid_T )
212
213	CALL dia_nam( clhstnam, nwrite, 'grid_T' )
214	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
215	CALL histbeg( clhstnam, jpi, glamt, jpj, gphit, & ! Horizontal grid: glamt and gphit
216	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
217	& nit000-1, zjulian, zdt, nh_T, nid_T, domain_id=nidom, snc4chunks=snc4set )
218	CALL histvert( nid_T, "deptht", "Vertical T levels", & ! Vertical grid: gdept
219	& "m", ipk, gdept_1d, nz_T, "down" )
220	! ! Index of ocean points
221	CALL wheneq( jpi*jpj , tmask, 1, 1., ndex_hT, ndim_hT ) ! surface
222
223	! Define the U grid FILE ( nid_U )
224
225	CALL dia_nam( clhstnam, nwrite, 'grid_U' )
226	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam ! filename
227	CALL histbeg( clhstnam, jpi, glamu, jpj, gphiu, & ! Horizontal grid: glamu and gphiu
228	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
229	& nit000-1, zjulian, zdt, nh_U, nid_U, domain_id=nidom, snc4chunks=snc4set )
230	CALL histvert( nid_U, "depthu", "Vertical U levels", & ! Vertical grid: gdept
231	& "m", ipk, gdept_1d, nz_U, "down" )
232	! ! Index of ocean points
233	CALL wheneq( jpi*jpj , umask, 1, 1., ndex_hU, ndim_hU ) ! surface
234
235	! Define the V grid FILE ( nid_V )
236
237	CALL dia_nam( clhstnam, nwrite, 'grid_V' ) ! filename
238	IF(lwp) WRITE(numout,*) " Name of NETCDF file ", clhstnam
239	CALL histbeg( clhstnam, jpi, glamv, jpj, gphiv, & ! Horizontal grid: glamv and gphiv
240	& iimi, iima-iimi+1, ijmi, ijma-ijmi+1, &
241	& nit000-1, zjulian, zdt, nh_V, nid_V, domain_id=nidom, snc4chunks=snc4set )
242	CALL histvert( nid_V, "depthv", "Vertical V levels", & ! Vertical grid : gdept
243	& "m", ipk, gdept_1d, nz_V, "down" )
244	! ! Index of ocean points
245	CALL wheneq( jpi*jpj , vmask, 1, 1., ndex_hV, ndim_hV ) ! surface
246
247	! No W grid FILE
248
249	! Declare all the output fields as NETCDF variables
250
251	! !!! nid_T : 3D
252	CALL histdef( nid_T, "sst_m", "Sea Surface temperature" , "C" , & ! sst
253	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
254	CALL histdef( nid_T, "sss_m", "Sea Surface Salinity" , "PSU" , & ! sss
255	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
256	CALL histdef( nid_T, "sowaflup", "Net Upward Water Flux" , "Kg/m2/s", & ! (emp-rnf)
257	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
258	CALL histdef( nid_T, "sosfldow", "downward salt flux" , "PSU/m2/s", & ! (sfx)
259	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
260	CALL histdef( nid_T, "sohefldo", "Net Downward Heat Flux" , "W/m2" , & ! qns + qsr
261	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
262	CALL histdef( nid_T, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! qsr
263	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
264	CALL histdef( nid_T, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
265	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
266	CALL histdef( nid_T, "sowindsp", "wind speed at 10m" , "m/s" , & ! wndm
267	& jpi, jpj, nh_T, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
268
269	CALL histend( nid_T, snc4chunks=snc4set )
270
271	! !!! nid_U : 3D
272	CALL histdef( nid_U, "ssu_m", "Velocity component in x-direction", "m/s" , & ! ssu
273	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
274	CALL histdef( nid_U, "sozotaux", "Wind Stress along i-axis" , "N/m2" , & ! utau
275	& jpi, jpj, nh_U, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
276
277	CALL histend( nid_U, snc4chunks=snc4set )
278
279	! !!! nid_V : 3D
280	CALL histdef( nid_V, "ssv_m", "Velocity component in y-direction", "m/s", & ! ssv_m
281	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
282	CALL histdef( nid_V, "sometauy", "Wind Stress along j-axis" , "N/m2" , & ! vtau
283	& jpi, jpj, nh_V, 1 , 1, 1 , - 99, 32, clop, zsto, zout )
284
285	CALL histend( nid_V, snc4chunks=snc4set )
286
287	IF(lwp) WRITE(numout,*)
288	IF(lwp) WRITE(numout,*) 'End of NetCDF Initialization'
289	IF(ll_print) CALL FLUSH(numout )
290
291	ENDIF
292
293	! 2. Start writing data
294	! ---------------------
295
296	! ndex(1) est utilise ssi l'avant dernier argument est diffferent de
297	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
298	! donne le nombre d'elements, et ndex la liste des indices a sortir
299
300	IF( lwp .AND. MOD( itmod, nwrite ) == 0 ) THEN
301	WRITE(numout,*) 'dia_wri : write model outputs in NetCDF files at ', kt, 'time-step'
302	WRITE(numout,*) '~~~~~~ '
303	ENDIF
304
305	! Write fields on T grid
306	CALL histwrite( nid_T, "sst_m", it, sst_m, ndim_hT, ndex_hT ) ! sea surface temperature
307	CALL histwrite( nid_T, "sss_m", it, sss_m, ndim_hT, ndex_hT ) ! sea surface salinity
308	CALL histwrite( nid_T, "sowaflup", it, (emp - rnf ) , ndim_hT, ndex_hT ) ! upward water flux
309	CALL histwrite( nid_T, "sosfldow", it, sfx , ndim_hT, ndex_hT ) ! downward salt flux
310	! (includes virtual salt flux beneath ice
311	! in linear free surface case)
312
313	CALL histwrite( nid_T, "sohefldo", it, qns + qsr , ndim_hT, ndex_hT ) ! total heat flux
314	CALL histwrite( nid_T, "soshfldo", it, qsr , ndim_hT, ndex_hT ) ! solar heat flux
315	CALL histwrite( nid_T, "soicecov", it, fr_i , ndim_hT, ndex_hT ) ! ice fraction
316	CALL histwrite( nid_T, "sowindsp", it, wndm , ndim_hT, ndex_hT ) ! wind speed
317
318	! Write fields on U grid
319	CALL histwrite( nid_U, "ssu_m" , it, ssu_m , ndim_hU, ndex_hU ) ! i-current speed
320	CALL histwrite( nid_U, "sozotaux", it, utau , ndim_hU, ndex_hU ) ! i-wind stress
321
322	! Write fields on V grid
323	CALL histwrite( nid_V, "ssv_m" , it, ssv_m , ndim_hV, ndex_hV ) ! j-current speed
324	CALL histwrite( nid_V, "sometauy", it, vtau , ndim_hV, ndex_hV ) ! j-wind stress
325
326	! 3. Close all files
327	! ---------------------------------------
328	IF( kt == nitend ) THEN
329	CALL histclo( nid_T )
330	CALL histclo( nid_U )
331	CALL histclo( nid_V )
332	ENDIF
333	!
334	IF( nn_timing == 1 ) CALL timing_stop('dia_wri')
335	!
336	END SUBROUTINE dia_wri
337	# endif
338
339	#endif
340
341	SUBROUTINE dia_wri_state( cdfile_name, kt )
342	!!---------------------------------------------------------------------
343	!! * ROUTINE dia_wri_state *
344	!!
345	!! ** Purpose : create a NetCDF file named cdfile_name which contains
346	!! the instantaneous ocean state and forcing fields.
347	!! Used to find errors in the initial state or save the last
348	!! ocean state in case of abnormal end of a simulation
349	!!
350	!! ** Method : NetCDF files using ioipsl
351	!! File 'output.init.nc' is created if ninist = 1 (namelist)
352	!! File 'output.abort.nc' is created in case of abnormal job end
353	!!----------------------------------------------------------------------
354	CHARACTER (len=* ), INTENT( in ) :: cdfile_name ! name of the file created
355	INTEGER , INTENT( in ) :: kt ! ocean time-step index
356	!!
357	CHARACTER (len=32) :: clname
358	CHARACTER (len=40) :: clop
359	INTEGER :: id_i , nz_i, nh_i
360	INTEGER, DIMENSION(1) :: idex ! local workspace
361	REAL(wp) :: zsto, zout, zmax, zjulian, zdt
362	!!----------------------------------------------------------------------
363	!
364	IF( nn_timing == 1 ) CALL timing_start('dia_wri_state')
365
366	! 0. Initialisation
367	! -----------------
368
369	! Define name, frequency of output and means
370	clname = cdfile_name
371	IF( .NOT. Agrif_Root() ) clname = TRIM(Agrif_CFixed())//'_'//TRIM(clname)
372	zdt = rdt
373	zsto = rdt
374	clop = "inst(x)" ! no use of the mask value (require less cpu time)
375	zout = rdt
376	zmax = ( nitend - nit000 + 1 ) * zdt
377
378	IF(lwp) WRITE(numout,*)
379	IF(lwp) WRITE(numout,*) 'dia_wri_state : single instantaneous ocean state'
380	IF(lwp) WRITE(numout,*) '~~~~~~~~~~~~~ and forcing fields file created '
381	IF(lwp) WRITE(numout,*) ' and named :', clname, '.nc'
382
383
384	! 1. Define NETCDF files and fields at beginning of first time step
385	! -----------------------------------------------------------------
386
387	! Compute julian date from starting date of the run
388	CALL ymds2ju( nyear, nmonth, nday, rdt, zjulian ) ! time axis
389	zjulian = zjulian - adatrj ! set calendar origin to the beginning of the experiment
390	CALL histbeg( clname, jpi, glamt, jpj, gphit, &
391	1, jpi, 1, jpj, nit000-1, zjulian, zdt, nh_i, id_i, domain_id=nidom, snc4chunks=snc4set ) ! Horizontal grid : glamt and gphit
392	CALL histvert( id_i, "deptht", "Vertical T levels", & ! Vertical grid : gdept
393	"m", jpk, gdept_1d, nz_i, "down")
394
395	! Declare all the output fields as NetCDF variables
396
397	CALL histdef( id_i, "sowaflup", "Net Upward Water Flux" , "Kg/m2/S", & ! net freshwater
398	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
399	CALL histdef( id_i, "sohefldo", "Net Downward Heat Flux", "W/m2" , & ! net heat flux
400	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
401	CALL histdef( id_i, "soshfldo", "Shortwave Radiation" , "W/m2" , & ! solar flux
402	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
403	CALL histdef( id_i, "soicecov", "Ice fraction" , "[0,1]" , & ! fr_i
404	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
405	CALL histdef( id_i, "sozotaux", "Zonal Wind Stress" , "N/m2" , & ! i-wind stress
406	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
407	CALL histdef( id_i, "sometauy", "Meridional Wind Stress", "N/m2" , & ! j-wind stress
408	& jpi, jpj, nh_i, 1 , 1, 1 , -99 , 32, clop, zsto, zout )
409
410	#if defined key_lim2
411	CALL lim_wri_state_2( kt, id_i, nh_i )
412	#else
413	CALL histend( id_i, snc4chunks=snc4set )
414	#endif
415
416	! 2. Start writing data
417	! ---------------------
418	! idex(1) est utilise ssi l'avant dernier argument est diffferent de
419	! la taille du tableau en sortie. Dans ce cas , l'avant dernier argument
420	! donne le nombre d'elements, et idex la liste des indices a sortir
421	idex(1) = 1 ! init to avoid compil warning
422
423	! Write all fields on T grid
424	CALL histwrite( id_i, "sowaflup", kt, emp , jpi*jpj , idex ) ! freshwater budget
425	CALL histwrite( id_i, "sohefldo", kt, qsr + qns , jpi*jpj , idex ) ! total heat flux
426	CALL histwrite( id_i, "soshfldo", kt, qsr , jpi*jpj , idex ) ! solar heat flux
427	CALL histwrite( id_i, "soicecov", kt, fr_i , jpi*jpj , idex ) ! ice fraction
428	CALL histwrite( id_i, "sozotaux", kt, utau , jpi*jpj , idex ) ! i-wind stress
429	CALL histwrite( id_i, "sometauy", kt, vtau , jpi*jpj , idex ) ! j-wind stress
430
431	! 3. Close the file
432	! -----------------
433	CALL histclo( id_i )
434	#if ! defined key_iomput && ! defined key_dimgout
435	IF( ninist /= 1 ) THEN
436	CALL histclo( nid_T )
437	CALL histclo( nid_U )
438	CALL histclo( nid_V )
439	ENDIF
440	#endif
441
442	IF( nn_timing == 1 ) CALL timing_stop('dia_wri_state')
443	!
444
445	END SUBROUTINE dia_wri_state
446	!!======================================================================
447	END MODULE diawri

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

Download in other formats: