MODULE diamlr
!!======================================================================
!! *** MODULE diamlr ***
!! Management of the IOM context for multiple-linear-regression analysis
!!======================================================================
!! History : 4.0 ! 2019 (S. Mueller) Original code
!!----------------------------------------------------------------------
USE par_kind
USE par_oce , ONLY : wp, jpi, jpj
USE phycst , ONLY : rpi
USE dom_oce , ONLY : adatrj
USE tide_mod
!
USE in_out_manager , ONLY : lwp, numout, ln_timing
USE iom , ONLY : iom_put, iom_use, iom_update_file_name
USE timing , ONLY : timing_start, timing_stop
#if defined key_xios
USE xios
#endif
IMPLICIT NONE
PRIVATE
LOGICAL, PUBLIC :: lk_diamlr = .FALSE. !: ===>>> NOT a DOCTOR norm name : use l_diamlr
! lk_ is used only for logical controlled by a CPP key
PUBLIC :: dia_mlr_init, dia_mlr_iom_init, dia_mlr
!! * Substitutions
# include "do_loop_substitute.h90"
!!----------------------------------------------------------------------
!! NEMO/OCE 4.0 , NEMO Consortium (2019)
!! $Id$
!! Software governed by the CeCILL license (see ./LICENSE)
!!----------------------------------------------------------------------
CONTAINS
SUBROUTINE dia_mlr_init
!!----------------------------------------------------------------------
!! *** ROUTINE dia_mlr_init ***
!!
!! ** Purpose : initialisation of IOM context management for
!! multiple-linear-regression analysis
!!
!!----------------------------------------------------------------------
!
lk_diamlr = .TRUE.
!
IF(lwp) THEN
WRITE(numout, *)
WRITE(numout, *) 'dia_mlr_init : initialisation of IOM context management for'
WRITE(numout, *) '~~~~~~~~~~~~ multiple-linear-regression analysis'
END IF
!
END SUBROUTINE dia_mlr_init
SUBROUTINE dia_mlr_iom_init
!!----------------------------------------------------------------------
!! *** ROUTINE dia_mlr_iom_init ***
!!
!! ** Purpose : IOM context setup for multiple-linear-regression
!! analysis
!!
!!----------------------------------------------------------------------
#if defined key_xios
TYPE(xios_fieldgroup) :: slxhdl_fldgrp
TYPE(xios_filegroup) :: slxhdl_filgrp
TYPE(xios_field), ALLOCATABLE, DIMENSION(:) :: slxhdl_regs, slxhdl_flds
TYPE(xios_field) :: slxhdl_fld
TYPE(xios_file) :: slxhdl_fil
LOGICAL :: llxatt_enabled, llxatt_comment
CHARACTER(LEN=256) :: clxatt_expr, clxatt_comment
CHARACTER(LEN=32) :: clxatt_name1, clxatt_name2
CHARACTER(LEN=32) :: clxatt_gridref, clxatt_fieldref
INTEGER, PARAMETER :: jpscanmax = 999
INTEGER :: ireg, ifld
CHARACTER(LEN=3) :: cl3i
CHARACTER(LEN=6) :: cl6a
CHARACTER(LEN=7) :: cl7a
CHARACTER(LEN=1) :: clgt
CHARACTER(LEN=2) :: clgd
CHARACTER(LEN=25) :: clfloat
CHARACTER(LEN=32) :: clrepl
INTEGER :: jl, jm, jn
INTEGER :: itide ! Number of available tidal components
REAL(wp) :: ztide_phase ! Tidal-constituent phase at adatrj=0
CHARACTER (LEN=4), DIMENSION(jpmax_harmo) :: ctide_selected = 'n/a '
TYPE(tide_harmonic), DIMENSION(:), POINTER :: stideconst
IF(lwp) THEN
WRITE(numout, *)
WRITE(numout, *) 'dia_mlr_iom_init : IOM context setup for multiple-linear-regression'
WRITE(numout, *) '~~~~~~~~~~~~~~~~'
END IF
! Get handles to multiple-linear-regression analysis configuration (field
! group 'diamrl_fields' and file group 'diamlr_files'); if no suitable
! configuration is found, disable diamlr
IF ( lk_diamlr .AND. xios_is_valid_fieldgroup( "diamlr_fields" ) .AND. xios_is_valid_field( "diamlr_time" ) .AND. &
& xios_is_valid_filegroup( "diamlr_files" ) ) THEN
CALL xios_get_handle("diamlr_fields", slxhdl_fldgrp)
CALL xios_get_handle("diamlr_files", slxhdl_filgrp)
ELSE
IF (lwp) THEN
WRITE(numout, *) "diamlr: configuration not found or incomplete (field group 'diamlr_fields'"
WRITE(numout, *) " and/or file group 'diamlr_files' and/or field 'diamlr_time' missing);"
WRITE(numout, *) " disabling output for multiple-linear-regression analysis."
END IF
lk_diamlr = .FALSE.
END IF
! Set up IOM context for multiple-linear-regression analysis
IF ( lk_diamlr ) THEN
! Set up output files for grid types scalar, grid_T, grid_U, grid_V,
! and grid_W
DO jm = 1, 5
SELECT CASE( jm )
CASE( 1 )
cl6a = 'scalar'
CASE( 2 )
cl6a = 'grid_T'
CASE( 3 )
cl6a = 'grid_U'
CASE( 4 )
cl6a = 'grid_V'
CASE( 5 )
cl6a = 'grid_W'
END SELECT
CALL xios_add_child ( slxhdl_filgrp, slxhdl_fil, "diamlr_file_"//cl6a )
CALL xios_set_attr ( slxhdl_fil, name_suffix="_diamlr_"//cl6a, &
& description="Intermediate output for multiple-linear-regression analysis - "//cl6a )
CALL iom_update_file_name( "diamlr_file_"//cl6a )
END DO
! Compile lists of active regressors and of fields selected for
! analysis (fields "diamlr_r<nnn>" and "diamlr_f<nnn>", where <nnn> is
! a 3-digit integer); also carry out placeholder substitution of tidal
! parameters in regressor expressions
!
ALLOCATE( slxhdl_regs( jpscanmax ), slxhdl_flds( jpscanmax ) )
ireg = 0
ifld = 0
!
IF ( ln_tide ) THEN
! Retrieve information (frequency, phase, nodal correction) about all
! available tidal constituents for placeholder substitution below
! Warning: we must use the same character length in an array constructor (at least for gcc compiler)
ctide_selected(1:34) = (/ 'Mf ', 'Mm ', 'Ssa ', 'Mtm ', 'Msf ', &
& 'Msqm', 'Sa ', 'K1 ', 'O1 ', 'P1 ', &
& 'Q1 ', 'J1 ', 'S1 ', 'M2 ', 'S2 ', 'N2 ', &
& 'K2 ', 'nu2 ', 'mu2 ', '2N2 ', 'L2 ', &
& 'T2 ', 'eps2', 'lam2', 'R2 ', 'M3 ', &
& 'MKS2', 'MN4 ', 'MS4 ', 'M4 ', 'N4 ', &
& 'S4 ', 'M6 ', 'M8 ' /)
CALL tide_init_harmonics(ctide_selected, stideconst)
itide = size(stideconst)
ELSE
itide = 0
ENDIF
DO jm = 1, jpscanmax
WRITE (cl3i, '(i3.3)') jm
! Look for regressor
IF ( xios_is_valid_field( "diamlr_r"//cl3i ) ) THEN
CALL xios_get_handle( "diamlr_r"//cl3i, slxhdl_regs(ireg+1) )
! Retrieve pre-configured value of "enabled" attribute and
! regressor expression
CALL xios_get_attr ( slxhdl_regs(ireg+1), enabled=llxatt_enabled, expr=clxatt_expr )
! If enabled, keep handle in list of active regressors; also
! substitute placeholders for tidal frequencies, phases, and
! nodal corrections in regressor expressions
IF ( llxatt_enabled ) THEN
! Substitution of placeholders for tidal-constituent
! parameters (amplitudes, angular veloccities, nodal phase
! correction) with values that have been obtained from the
! tidal-forcing implementation (if enabled)
DO jn = 1, itide
! Compute phase of tidal constituent (incl. current nodal
! correction) at the start of the model run (i.e. for
! adatrj=0)
ztide_phase = MOD( stideconst(jn)%u + stideconst(jn)%v0 - adatrj * 86400.0_wp * stideconst(jn)%omega, &
& 2.0_wp * rpi )
clrepl = "__TDE_"//TRIM( stideconst(jn)%cname_tide )//"_omega__"
DO WHILE ( INDEX( clxatt_expr, TRIM( clrepl ) ) > 0 )
WRITE (clfloat, '(e25.18)') stideconst(jn)%omega
jl = INDEX( clxatt_expr, TRIM( clrepl ) )
clxatt_expr = clxatt_expr(1:jl - 1)//clfloat// &
& clxatt_expr(jl + LEN( TRIM( clrepl ) ):LEN( TRIM( clxatt_expr ) ))
END DO
clrepl = "__TDE_"//TRIM( stideconst(jn)%cname_tide )//"_phase__"
DO WHILE ( INDEX( clxatt_expr, TRIM( clrepl ) ) > 0 )
WRITE (clfloat, '(e25.18)') ztide_phase
jl = INDEX( clxatt_expr, TRIM( clrepl ) )
clxatt_expr = clxatt_expr(1:jl - 1)//clfloat// &
& clxatt_expr(jl + LEN( TRIM( clrepl ) ):LEN( TRIM( clxatt_expr ) ))
END DO
clrepl = "__TDE_"//TRIM( stideconst(jn)%cname_tide )//"_amplitude__"
DO WHILE (INDEX( clxatt_expr, TRIM( clrepl ) ) > 0 )
WRITE (clfloat, '(e25.18)') stideconst(jn)%f
jl = INDEX( clxatt_expr, TRIM( clrepl ) )
clxatt_expr = clxatt_expr(1:jl - 1)//clfloat// &
& clxatt_expr(jl + LEN( TRIM( clrepl ) ):LEN( TRIM( clxatt_expr ) ))
END DO
END DO
! Set standard value for comment attribute, including possible
! existing comment added in parantheses
CALL xios_is_defined_attr( slxhdl_regs(ireg+1), comment=llxatt_comment )
IF ( llxatt_comment ) THEN
CALL xios_get_attr( slxhdl_regs(ireg+1), comment=clxatt_comment )
clxatt_comment = "Regressor "//cl3i//" ("//TRIM( clxatt_comment )//") "
ELSE
clxatt_comment = "Regressor "//cl3i
END IF
! Set name attribute (and overwrite possible pre-configured
! name) with field id to enable id string retrieval from
! stored handle below, re-set expression with possible
! substitutions, and set or re-set comment attribute
CALL xios_set_attr ( slxhdl_regs(ireg+1), name="diamlr_r"//cl3i, expr=TRIM( clxatt_expr ), &
& comment=TRIM( clxatt_comment ) )
ireg = ireg + 1 ! Accept regressor in list of active regressors
END IF
END IF
! Look for field
IF ( xios_is_valid_field( "diamlr_f"//cl3i ) ) THEN
CALL xios_get_handle( "diamlr_f"//cl3i, slxhdl_flds(ifld+1) )
! Retrieve pre-configured value of "enabled" attribute
CALL xios_get_attr ( slxhdl_flds(ifld+1), enabled=llxatt_enabled )
! If enabled, keep handle in list of fields selected for analysis
IF ( llxatt_enabled ) THEN
! Set name attribute (and overwrite possible pre-configured name)
! with field id to enable id string retrieval from stored handle
! below
CALL xios_set_attr ( slxhdl_flds(ifld+1), name="diamlr_f"//cl3i )
ifld = ifld + 1 ! Accept field in list of fields selected for analysis
END IF
END IF
END DO
! Output number of active regressors and fields selected for analysis
IF ( lwp ) WRITE(numout,'(a,i3,a)' ) 'diamlr: ', ireg, ' active regressors found'
IF ( lwp ) WRITE(numout,'(a,i3,a)' ) 'diamlr: ', ifld, ' fields selected for analysis'
! Set up output of minimum, maximum, and average values of the time
! variable available for the computation of regressors (diamlr_time)
CALL xios_get_handle( "diamlr_file_scalar", slxhdl_fil )
CALL xios_add_child ( slxhdl_fil, slxhdl_fld, "diamlr_time_average" )
!$AGRIF_DO_NOT_TREAT
CALL xios_set_attr ( slxhdl_fld, standard_name="diamlr_time", &
& long_name="Elapsed model time at start of regression interval", &
& unit="s", operation="average", field_ref="diamlr_time", &
& grid_ref="diamlr_grid_2D_to_scalar" )
!$AGRIF_END_DO_NOT_TREAT
CALL xios_add_child ( slxhdl_fil, slxhdl_fld, "diamlr_time_minimum" )
!$AGRIF_DO_NOT_TREAT
CALL xios_set_attr ( slxhdl_fld, standard_name="diamlr_time", &
& long_name="Elapsed model time at start of regression interval", &
& unit="s", operation="minimum", field_ref="diamlr_time", &
& grid_ref="diamlr_grid_2D_to_scalar" )
!$AGRIF_END_DO_NOT_TREAT
CALL xios_add_child ( slxhdl_fil, slxhdl_fld, "diamlr_time_maximum" )
!$AGRIF_DO_NOT_TREAT
CALL xios_set_attr ( slxhdl_fld, standard_name="diamlr_time", &
& long_name="Elapsed model time at start of regression interval", &
& unit="s", operation="maximum", field_ref="diamlr_time", &
& grid_ref="diamlr_grid_2D_to_scalar" )
!$AGRIF_END_DO_NOT_TREAT
! For each active regressor:
DO jm = 1, ireg
! i) set up 2-dimensional and 3-dimensional versions of the
! regressors; explicitely set "enabled" attribute; note, while
! the scalar versions of regressors are part of the
! configuration, the respective 2-dimensional versions take
! over the defining expression, while the scalar and
! 3-dimensional versions are simply obtained via grid
! transformations from the 2-dimensional version.
CALL xios_get_attr ( slxhdl_regs( jm ), name=clxatt_name1, expr=clxatt_expr, &
& enabled=llxatt_enabled, comment=clxatt_comment )
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_T_2D" )
CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_T_2D", &
& field_ref="diamlr_time", enabled=llxatt_enabled )
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_U_2D" )
CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_U_2D", &
& field_ref="diamlr_time", enabled=llxatt_enabled )
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_V_2D" )
CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_V_2D", &
& field_ref="diamlr_time", enabled=llxatt_enabled )
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_W_2D" )
CALL xios_set_attr ( slxhdl_fld, expr=TRIM( clxatt_expr ), grid_ref="diamlr_grid_W_2D", &
& field_ref="diamlr_time", enabled=llxatt_enabled )
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_T_3D")
CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_T_3D", &
& field_ref=TRIM( clxatt_name1 )//"_grid_T_2D", enabled=llxatt_enabled)
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_U_3D")
CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_U_3D", &
& field_ref=TRIM( clxatt_name1 )//"_grid_U_2D", enabled=llxatt_enabled)
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_V_3D")
CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_V_3D", &
& field_ref=TRIM( clxatt_name1 )//"_grid_V_2D", enabled=llxatt_enabled)
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"_grid_W_3D")
CALL xios_set_attr ( slxhdl_fld, expr="this", grid_ref="diamlr_grid_2D_to_grid_W_3D", &
& field_ref=TRIM( clxatt_name1 )//"_grid_W_2D", enabled=llxatt_enabled)
CALL xios_set_attr ( slxhdl_regs(jm), expr="this", grid_ref="diamlr_grid_2D_to_scalar", &
& field_ref=TRIM( clxatt_name1 )//"_grid_T_2D", enabled=llxatt_enabled)
! ii) set up output of active regressors, including metadata
CALL xios_get_handle( "diamlr_file_scalar", slxhdl_fil )
! Add regressor to output file
CALL xios_add_child ( slxhdl_fil, slxhdl_fld, TRIM( clxatt_name1 ) )
CALL xios_set_attr ( slxhdl_fld, standard_name=TRIM( clxatt_comment ), long_name=TRIM( clxatt_expr ), &
& operation="average" )
! iii) set up the output of scalar products with itself and with
! other active regressors
CALL xios_get_attr ( slxhdl_regs(jm), name=clxatt_name1 )
DO jn = 1, jm
! Field for product between regressors
CALL xios_get_attr ( slxhdl_regs(jn), name=clxatt_name2 )
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ) )
! Set appropriate name attribute to avoid the possibility of
! using an inappropriate inherited name attribute as the variable
! name in the output file
CALL xios_set_attr ( slxhdl_fld, &
& name=TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ), &
& grid_ref="diamlr_grid_scalar", &
& expr="this * "//TRIM( clxatt_name2 ), &
& field_ref=TRIM( clxatt_name1 ), &
& enabled=llxatt_enabled, &
& long_name="Scalar product of regressor "//TRIM( clxatt_name1 )// &
& " and regressor "//TRIM( clxatt_name2 ), &
& standard_name=TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ), &
& operation="accumulate")
! Add regressor-product field to output file
CALL xios_add_child ( slxhdl_fil, slxhdl_fld, TRIM( clxatt_name1 )//"."//TRIM( clxatt_name2 ) )
END DO
! iv) set up definitions for the output of scalar products with
! fields selected for analysis
DO jn = 1, ifld
CALL xios_get_attr ( slxhdl_flds(jn), name=clxatt_name2, field_ref=clxatt_fieldref )
CALL xios_get_handle( TRIM( clxatt_fieldref ), slxhdl_fld )
CALL xios_get_attr ( slxhdl_fld, grid_ref=clxatt_gridref )
clgt="T"
IF ( INDEX( clxatt_gridref, "_U_" ) > 0 ) clgt="U"
IF ( INDEX( clxatt_gridref, "_V_" ) > 0 ) clgt="V"
IF ( INDEX( clxatt_gridref, "_W_" ) > 0 ) clgt="W"
clgd="2D"
cl7a=""
IF ( INDEX( clxatt_gridref, "_3D" ) > 0 ) THEN
clgd="3D"
ELSE
cl7a="diamlr_"
END IF
CALL xios_add_child ( slxhdl_fldgrp, slxhdl_fld, TRIM( clxatt_name2 )//"."//TRIM( clxatt_name1 ) )
! Set appropriate name attribute to avoid the possibility of
! using an inappropriate inherited name attribute as the variable
! name in the output file; use metadata (standard_name and
! long_name) to refer to the id of the analysed field
CALL xios_set_attr ( slxhdl_fld, &
& name=TRIM( clxatt_name2 )//"."//TRIM( clxatt_name1 ), &
& expr="this * "//TRIM( clxatt_fieldref ), &
& grid_ref=cl7a//"grid_"//clgt//"_"//clgd, &
& field_ref=TRIM( clxatt_name1 )//"_grid_"//clgt//"_"//clgd, &
& enabled=llxatt_enabled, &
& long_name="Scalar product of "//TRIM( clxatt_fieldref )// &
& " and regressor "//TRIM( clxatt_name1 ), &
& standard_name=TRIM( clxatt_fieldref )//"."//TRIM( clxatt_name1 ), &
& operation="accumulate" )
CALL xios_get_handle( "diamlr_file_grid_"//clgt, slxhdl_fil )
CALL xios_add_child ( slxhdl_fil, slxhdl_fld, TRIM( clxatt_name2 )//"."//TRIM( clxatt_name1 ) )
END DO
END DO
! Release list of active regressors and fields selected for analysis
DEALLOCATE( slxhdl_regs, slxhdl_flds )
END IF
#else
IF( .FALSE. ) write(numout,*) 'dia_mlr_iom_init: should not see this' ! useless statement to avoid compiler warnings
#endif
END SUBROUTINE dia_mlr_iom_init
SUBROUTINE dia_mlr
!!----------------------------------------------------------------------
!! *** ROUTINE dia_mlr ***
!!
!! ** Purpose : update time used in multiple-linear-regression analysis
!!
!!----------------------------------------------------------------------
REAL(wp), DIMENSION(jpi,jpj) :: zadatrj2d
!!----------------------------------------------------------------------
IF( ln_timing ) CALL timing_start('dia_mlr')
! Update time to the continuous time since the start of the model run
! (value of adatrj converted to time in units of seconds)
!
! A 2-dimensional field of constant value is sent, and subsequently used directly
! or transformed to a scalar or a constant 3-dimensional field as required.
zadatrj2d(:,:) = adatrj*86400.0_wp
IF ( iom_use('diamlr_time') ) CALL iom_put('diamlr_time', zadatrj2d)
!
IF( ln_timing ) CALL timing_stop('dia_mlr')
!
END SUBROUTINE dia_mlr
!!======================================================================
END MODULE diamlr