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.

diamlr.F90 in NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DIA – NEMO

Context Navigation

source: NEMO/branches/2019/dev_r11879_ENHANCE-05_SimonM-Harmonic_Analysis/src/OCE/DIA/diamlr.F90 @ 12012

Last change on this file since 12012 was 12012, checked in by smueller, 4 years ago
Further modification that completes the enabling of the grid-type propagation started with [12010] (ticket #2175)
File size: 22.8 KB

Line
1	MODULE diamlr
2	!!======================================================================
3	!! * MODULE diamlr *
4	!! Management of the IOM context for multiple-linear-regression analysis
5	!!======================================================================
6	!! History : ! 2019 (S. Mueller)
7	!!----------------------------------------------------------------------
8
9	USE par_oce , ONLY : wp, jpi, jpj
10	USE phycst , ONLY : rpi
11	USE in_out_manager , ONLY : lwp, numout, ln_timing
12	USE iom , ONLY : iom_put, iom_use, iom_update_file_name
13	USE dom_oce , ONLY : adatrj
14	USE timing , ONLY : timing_start, timing_stop
15	USE xios
16	USE tide_mod , ONLY : tide_harmo, jpmax_harmo, Wave
17
18	IMPLICIT NONE
19	PRIVATE
20
21	LOGICAL, PUBLIC :: lk_diamlr = .FALSE.
22
23	PUBLIC :: dia_mlr_init, dia_mlr_iom_init, dia_mlr
24
25	!!----------------------------------------------------------------------
26	!! NEMO/OCE 4.0 , NEMO Consortium (2019)
27	!! $Id$
28	!! Software governed by the CeCILL license (see ./LICENSE)
29	!!----------------------------------------------------------------------
30	CONTAINS
31
32	SUBROUTINE dia_mlr_init
33	!!----------------------------------------------------------------------
34	!! * ROUTINE dia_mlr_init *
35	!!
36	!! ** Purpose : initialisation of IOM context management for
37	!! multiple-linear-regression analysis
38	!!
39	!!----------------------------------------------------------------------
40
41	lk_diamlr = .TRUE.
42
43	IF(lwp) THEN
44	WRITE(numout, *)
45	WRITE(numout, *) 'dia_mlr_init : initialisation of IOM context management for'
46	WRITE(numout, *) '~~~~~~~~~~~~ multiple-linear-regression analysis'
47	END IF
48
49	END SUBROUTINE dia_mlr_init
50
51	SUBROUTINE dia_mlr_iom_init
52	!!----------------------------------------------------------------------
53	!! * ROUTINE dia_mlr_iom_init *
54	!!
55	!! ** Purpose : IOM context setup for multiple-linear-regression
56	!! analysis
57	!!
58	!!----------------------------------------------------------------------
59
60	TYPE(xios_fieldgroup) :: slxhdl_fldgrp
61	TYPE(xios_filegroup) :: slxhdl_filgrp
62	TYPE(xios_field), ALLOCATABLE, DIMENSION(:) :: slxhdl_regs, slxhdl_flds
63	TYPE(xios_field) :: slxhdl_fld
64	TYPE(xios_file) :: slxhdl_fil
65	LOGICAL :: llxatt_enabled, llxatt_comment
66	CHARACTER(LEN=256) :: clxatt_expr, clxatt_comment
67	CHARACTER(LEN=32) :: clxatt_name1, clxatt_name2
68	CHARACTER(LEN=32) :: clxatt_gridref, clxatt_fieldref
69	INTEGER, PARAMETER :: jpscanmax = 999
70	INTEGER :: ireg, ifld
71	CHARACTER(LEN=3) :: cl3i
72	CHARACTER(LEN=6) :: cl6a
73	CHARACTER(LEN=7) :: cl7a
74	CHARACTER(LEN=1) :: clgt
75	CHARACTER(LEN=2) :: clgd
76	CHARACTER(LEN=25) :: clfloat
77	CHARACTER(LEN=32) :: clrepl
78	INTEGER :: jl, jm, jn
79	INTEGER :: itide ! Number of available tidal components
80	INTEGER, ALLOCATABLE, DIMENSION(:) :: itide_const ! Index list of selected tidal constituents
81	REAL(wp), ALLOCATABLE, DIMENSION(:) :: ztide_omega, ztide_u, & ! Tidal frequency, phase, nodal correction
82	& ztide_v, ztide_f
83	REAL(wp) :: ztide_phase ! Tidal-constituent phase at adatrj=0
84
85	IF(lwp) THEN
86	WRITE(numout, *)
87	WRITE(numout, *) 'dia_mlr_iom_init : IOM context setup for multiple-linear-regression'
88	WRITE(numout, *) '~~~~~~~~~~~~~~~~'
89	END IF
90
91	! Get handles to multiple-linear-regression analysis configuration (field
92	! group 'diamrl_fields' and file group 'diamlr_files'); if no suitable
93	! configuration is found, disable diamlr
94	IF ( lk_diamlr .AND. xios_is_valid_fieldgroup( "diamlr_fields" ) .AND. xios_is_valid_field( "diamlr_time" ) .AND. &
95	& xios_is_valid_filegroup( "diamlr_files" ) ) THEN
96	CALL xios_get_handle("diamlr_fields", slxhdl_fldgrp)
97	CALL xios_get_handle("diamlr_files", slxhdl_filgrp)
98	ELSE
99	IF (lwp) THEN
100	WRITE(numout, *) "diamlr: configuration not found or icomplete (field group 'diamlr_fields'"
101	WRITE(numout, *) " and/or file group 'diamlr_files' and/or field 'diamlr_time' missing);"
102	WRITE(numout, *) " disabling output for multiple-linear-regression analysis."
103	END IF
104	lk_diamlr = .FALSE.
105	END IF
106
107	! Set up IOM context for multiple-linear-regression analysis
108	IF ( lk_diamlr ) THEN
109
110	! Set up output files for grid types scalar, grid_T, grid_U, grid_V,
111	! and grid_W
112	DO jm = 1, 5
113	SELECT CASE( jm )
114	CASE( 1 )
115	cl6a = 'scalar'
116	CASE( 2 )
117	cl6a = 'grid_T'
118	CASE( 3 )
119	cl6a = 'grid_U'
120	CASE( 4 )
121	cl6a = 'grid_V'
122	CASE( 5 )
123	cl6a = 'grid_W'
124	END SELECT
125	CALL xios_add_child ( slxhdl_filgrp, slxhdl_fil, "diamlr_file_"//cl6a )
126	CALL xios_set_attr ( slxhdl_fil, name_suffix="_diamlr_"//cl6a, &
127	& description="Intermediate output for multiple-linear-regression analysis - "//cl6a )
128	CALL iom_update_file_name( "diamlr_file_"//cl6a )
129	END DO
130
131	! Compile lists of active regressors and of fields selected for
132	! analysis (fields "diamlr_r<nnn>" and "diamlr_f<nnn>", where <nnn> is
133	! a 3-digit integer); also carry out placeholder substitution of tidal
134	! parameters in regressor expressions
135	!
136	ALLOCATE( slxhdl_regs( jpscanmax ), slxhdl_flds( jpscanmax ) )
137	ireg = 0
138	ifld = 0
139	!
140	! Retrieve information (frequency, phase, nodal correction) about all
141	! available tidal constituents for placeholder substitution below
142	itide = jpmax_harmo
143	ALLOCATE(itide_const(itide), ztide_omega(itide), ztide_u(itide), ztide_v(itide), ztide_f(itide))
144	DO jn = 1, itide
145	itide_const(jn) = jn ! Select all available tidal constituents
146	END DO
147	CALL tide_harmo( ztide_omega, ztide_v, ztide_u, ztide_f, itide_const, itide )
148
149	DO jm = 1, jpscanmax
150	WRITE (cl3i, '(i3.3)') jm
151
152	! Look for regressor
153	IF ( xios_is_valid_field( "diamlr_r"//cl3i ) ) THEN
154
155	CALL xios_get_handle( "diamlr_r"//cl3i, slxhdl_regs(ireg+1) )
156	! Retrieve pre-configured value of "enabled" attribute and
157	! regressor expression
158	CALL xios_get_attr ( slxhdl_regs(ireg+1), enabled=llxatt_enabled, expr=clxatt_expr )
159	! If enabled, keep handle in list of active regressors; also
160	! substitute placeholders for tidal frequencies, phases, and
161	! nodal corrections in regressor expressions
162	IF ( llxatt_enabled ) THEN
163
164	! Substitution of placeholders for tidal-constituent
165	! parameters (amplitudes, angular veloccities, nodal phase
166	! correction) with values that have been obtained from the
167	! tidal-forcing implementation
168	DO jn = 1, itide
169	! Compute phase of tidal constituent (incl. current nodal
170	! correction) at the start of the model run (i.e. for
171	! adatrj=0)
172	ztide_phase = MOD( ztide_u(jn) + ztide_v(jn) - adatrj * 86400.0_wp * ztide_omega(jn), 2.0_wp * rpi )
173	clrepl = "__TDE_"//TRIM( Wave(jn)%cname_tide )//"_omega__"
174	DO WHILE ( INDEX( clxatt_expr, TRIM( clrepl ) ) > 0 )
175	WRITE (clfloat, '(e25.18)') ztide_omega(jn)
176	jl = INDEX( clxatt_expr, TRIM( clrepl ) )
177	clxatt_expr = clxatt_expr(1:jl - 1)//clfloat// &
178	& clxatt_expr(jl + LEN( TRIM( clrepl ) ):LEN( TRIM( clxatt_expr ) ))
179	END DO
180	clrepl = "__TDE_"//TRIM( Wave(jn)%cname_tide )//"_phase__"
181	DO WHILE ( INDEX( clxatt_expr, TRIM( clrepl ) ) > 0 )
182	WRITE (clfloat, '(e25.18)') ztide_phase
183	jl = INDEX( clxatt_expr, TRIM( clrepl ) )
184	clxatt_expr = clxatt_expr(1:jl - 1)//clfloat// &
185	& clxatt_expr(jl + LEN( TRIM( clrepl ) ):LEN( TRIM( clxatt_expr ) ))
186	END DO
187	clrepl = "__TDE_"//TRIM( Wave(jn)%cname_tide )//"_amplitude__"
188	DO WHILE (INDEX( clxatt_expr, TRIM( clrepl ) ) > 0 )
189	WRITE (clfloat, '(e25.18)') ztide_f(jn)
190	jl = INDEX( clxatt_expr, TRIM( clrepl ) )
191	clxatt_expr = clxatt_expr(1:jl - 1)//clfloat// &
192	& clxatt_expr(jl + LEN( TRIM( clrepl ) ):LEN( TRIM( clxatt_expr ) ))
193	END DO
194	END DO
195
196	! Set standard value for comment attribute, including possible
197	! existing comment added in parantheses
198	CALL xios_is_defined_attr( slxhdl_regs(ireg+1), comment=llxatt_comment )
199	IF ( llxatt_comment ) THEN
200	CALL xios_get_attr( slxhdl_regs(ireg+1), comment=clxatt_comment )
201	clxatt_comment = "Regressor "//cl3i//" ("//TRIM( clxatt_comment )//") "
202	ELSE
203	clxatt_comment = "Regressor "//cl3i
204	END IF
205
206	! Set name attribute (and overwrite possible pre-configured
207	! name) with field id to enable id string retrieval from
208	! stored handle below, re-set expression with possible
209	! substitutions, and set or re-set comment attribute
210	CALL xios_set_attr ( slxhdl_regs(ireg+1), name="diamlr_r"//cl3i, expr=TRIM( clxatt_expr ), &
211	& comment=TRIM( clxatt_comment ) )
212
213	ireg = ireg + 1 ! Accept regressor in list of active regressors
214
215	END IF
216	END IF
217
218	! Look for field
219	IF ( xios_is_valid_field( "diamlr_f"//cl3i ) ) THEN
220
221	CALL xios_get_handle( "diamlr_f"//cl3i, slxhdl_flds(ifld+1) )
222	! Retrieve pre-configured value of "enabled" attribute
223	CALL xios_get_attr ( slxhdl_flds(ifld+1), enabled=llxatt_enabled )
224	! If enabled, keep handle in list of fields selected for analysis
225	IF ( llxatt_enabled ) THEN
226
227	! Set name attribute (and overwrite possible pre-configured name)
228	! with field id to enable id string retrieval from stored handle
229	! below
230	CALL xios_set_attr ( slxhdl_flds(ifld+1), name="diamlr_f"//cl3i )
231
232	ifld = ifld + 1 ! Accept field in list of fields selected for analysis
233
234	END IF
235	END IF
236
237	END DO
238
239	! Release tidal data
240	DEALLOCATE( itide_const, ztide_omega, ztide_u, ztide_v, ztide_f )
241
242	! Output number of active regressors and fields selected for analysis
243	IF ( lwp ) WRITE(numout,'(a,i3,a)' ) 'diamlr: ', ireg, ' active regressors found'
244	IF ( lwp ) WRITE(numout,'(a,i3,a)' ) 'diamlr: ', ifld, ' fields selected for analysis'
245
246	! Set up output of minimum, maximum, and average values of the time
247	! variable available for the computation of regressors (diamlr_time)
248	CALL xios_get_handle( "diamlr_file_scalar", slxhdl_fil )
249	CALL xios_add_child ( slxhdl_fil, slxhdl_fld, "diamlr_time_average" )
250	CALL xios_set_attr ( slxhdl_fld, standard_name="diamlr_time", &
251	& long_name="Elapsed model time at start of regression interval", &
252	& unit="s", operation="average", field_ref="diamlr_time", &
253	& grid_ref="diamlr_grid_2D_to_scalar" )
254	CALL xios_add_child ( slxhdl_fil, slxhdl_fld, "diamlr_time_minimum" )
255	CALL xios_set_attr ( slxhdl_fld, standard_name="diamlr_time", &
256	& long_name="Elapsed model time at start of regression interval", &
257	& unit="s", operation="minimum", field_ref="diamlr_time", &
258	& grid_ref="diamlr_grid_2D_to_scalar" )
259	CALL xios_add_child ( slxhdl_fil, slxhdl_fld, "diamlr_time_maximum" )
260	CALL xios_set_attr ( slxhdl_fld, standard_name="diamlr_time", &
261	& long_name="Elapsed model time at start of regression interval", &
262	& unit="s", operation="maximum", field_ref="diamlr_time", &
263	& grid_ref="diamlr_grid_2D_to_scalar" )
264
265	! For each active regressor:
266	DO jm = 1, ireg
267
268	! i) set up 2-dimensional and 3-dimensional versions of the
269	! regressors; explicitely set "enabled" attribute; note, while
270	! the scalar versions of regressors are part of the
271	! configuration, the respective 2-dimensional versions take
272	! over the defining expression, while the scalar and
273	! 3-dimensional versions are simply obtained via grid
274	! transformations from the 2-dimensional version.
275	CALL xios_get_attr ( slxhdl_regs( jm ), name=clxatt_name1, expr=clxatt_expr, &
276	& enabled=llxatt_enabled, comment=clxatt_comment )
277	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_T_2D" )
278	CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_T_2D", &
279	& field_ref="diamlr_time", enabled=llxatt_enabled )
280	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_U_2D" )
281	CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_U_2D", &
282	& field_ref="diamlr_time", enabled=llxatt_enabled )
283	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_V_2D" )
284	CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_V_2D", &
285	& field_ref="diamlr_time", enabled=llxatt_enabled )
286	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_W_2D" )
287	CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_W_2D", &
288	& field_ref="diamlr_time", enabled=llxatt_enabled )
289	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_T_3D")
290	CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_T_3D", &
291	& field_ref=TRIM( clxatt_name1 )//"_grid_T_2D", enabled=llxatt_enabled)
292	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_U_3D")
293	CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_U_3D", &
294	& field_ref=TRIM( clxatt_name1 )//"_grid_U_2D", enabled=llxatt_enabled)
295	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_V_3D")
296	CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_V_3D", &
297	& field_ref=TRIM( clxatt_name1 )//"_grid_V_2D", enabled=llxatt_enabled)
298	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_W_3D")
299	CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_W_3D", &
300	& field_ref=TRIM( clxatt_name1 )//"_grid_W_2D", enabled=llxatt_enabled)
301	CALL xios_set_attr ( slxhdl_regs(jm), expr="this", grid_ref="diamlr_grid_2D_to_scalar", &
302	& field_ref=TRIM( clxatt_name1 )//"_grid_T_2D", enabled=llxatt_enabled)
303
304	! ii) set up output of active regressors, including metadata
305	CALL xios_get_handle( "diamlr_file_scalar", slxhdl_fil )
306	! Add regressor to output file
307	CALL xios_add_child ( slxhdl_fil, slxhdl_fld, TRIM( clxatt_name1 ) )
308	CALL xios_set_attr ( slxhdl_fld, standard_name=TRIM( clxatt_comment ), long_name=TRIM( clxatt_expr ), &
309	& operation="average" )
310
311	! iii) set up the output of scalar products with itself and with
312	! other active regressors
313	CALL xios_get_attr ( slxhdl_regs(jm), name=clxatt_name1 )
314	DO jn = 1, jm
315	! Field for product between regressors
316	CALL xios_get_attr ( slxhdl_regs(jn), name=clxatt_name2 )
317	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ) )
318	! Set appropriate name attribute to avoid the possibility of
319	! using an inappropriate inherited name attribute as the variable
320	! name in the output file
321	CALL xios_set_attr ( slxhdl_fld, &
322	& name=TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ), &
323	& grid_ref="diamlr_grid_scalar", &
324	& expr="this * "//TRIM( clxatt_name2 ), &
325	& field_ref=TRIM( clxatt_name1 ), &
326	& enabled=llxatt_enabled, &
327	& long_name="Scalar product of regressor "//TRIM( clxatt_name1 )// &
328	& " and regressor "//TRIM( clxatt_name2 ), &
329	& standard_name=TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ), &
330	& operation="accumulate")
331	! Add regressor-product field to output file
332	CALL xios_add_child ( slxhdl_fil, slxhdl_fld, TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ) )
333	END DO
334
335	! iv) set up definitions for the output of scalar products with
336	! fields selected for analysis
337	DO jn = 1, ifld
338	CALL xios_get_attr ( slxhdl_flds(jn), name=clxatt_name2, field_ref=clxatt_fieldref )
339	CALL xios_get_handle( TRIM( clxatt_fieldref ), slxhdl_fld )
340	CALL xios_get_attr ( slxhdl_fld, grid_ref=clxatt_gridref )
341	clgt="T"
342	IF ( INDEX( clxatt_gridref, "_U_" ) > 0 ) clgt="U"
343	IF ( INDEX( clxatt_gridref, "_V_" ) > 0 ) clgt="V"
344	IF ( INDEX( clxatt_gridref, "_W_" ) > 0 ) clgt="W"
345	clgd="2D"
346	cl7a=""
347	IF ( INDEX( clxatt_gridref, "_3D" ) > 0 ) THEN
348	clgd="3D"
349	ELSE
350	cl7a="diamlr_"
351	END IF
352	CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name2 )//"."//TRIM( clxatt_name1 ) )
353	! Set appropriate name attribute to avoid the possibility of
354	! using an inappropriate inherited name attribute as the variable
355	! name in the output file; use metadata (standard_name and
356	! long_name) to refer to the id of the analysed field
357	CALL xios_set_attr ( slxhdl_fld, &
358	& name=TRIM( clxatt_name2 )//"."//TRIM( clxatt_name1 ), &
359	& expr="this * "//TRIM( clxatt_fieldref ), &
360	& grid_ref=cl7a//"grid_"//clgt//"_"//clgd, &
361	& field_ref=TRIM( clxatt_name1 )//"_grid_"//clgt//"_"//clgd, &
362	& enabled=llxatt_enabled, &
363	& long_name="Scalar product of "//TRIM( clxatt_fieldref )// &
364	& " and regressor "//TRIM( clxatt_name1 ), &
365	& standard_name=TRIM( clxatt_fieldref )//"."//TRIM( clxatt_name1 ), &
366	& operation="accumulate" )
367	CALL xios_get_handle( "diamlr_file_grid_"//clgt, slxhdl_fil )
368	CALL xios_add_child ( slxhdl_fil, slxhdl_fld, TRIM( clxatt_name2 )//"."//TRIM( clxatt_name1 ) )
369	END DO
370
371	END DO
372
373	! Release list of active regressors and fields selected for analysis
374	DEALLOCATE( slxhdl_regs, slxhdl_flds )
375
376	END IF
377
378	END SUBROUTINE dia_mlr_iom_init
379
380	SUBROUTINE dia_mlr
381	!!----------------------------------------------------------------------
382	!! * ROUTINE dia_mlr *
383	!!
384	!! ** Purpose : update time used in multiple-linear-regression analysis
385	!!
386	!!----------------------------------------------------------------------
387
388	REAL(wp), DIMENSION(jpi,jpj) :: zadatrj2d
389
390	IF( ln_timing ) CALL timing_start('dia_mlr')
391
392	! Update time to the continuous time since the start of the model run
393	! (value of adatrj converted to time in units of seconds)
394	!
395	! A 2-dimensional field of constant value is sent, and subsequently used
396	! directly or transformed to a scalar or a constant 3-dimensional field as
397	! required.
398	zadatrj2d(:,:) = adatrj*86400.0_wp
399	IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d)
400
401	IF( ln_timing ) CALL timing_stop('dia_mlr')
402
403	END SUBROUTINE dia_mlr
404
405	END MODULE diamlr

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

Download in other formats: