! Copyright (c) 2004-2024 Lars Nerger ! ! This file is part of PDAF. ! ! PDAF is free software: you can redistribute it and/or modify ! it under the terms of the GNU Lesser General Public License ! as published by the Free Software Foundation, either version ! 3 of the License, or (at your option) any later version. ! ! PDAF is distributed in the hope that it will be useful, ! but WITHOUT ANY WARRANTY; without even the implied warranty of ! MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the ! GNU Lesser General Public License for more details. ! ! You should have received a copy of the GNU Lesser General Public ! License along with PDAF. If not, see <http://www.gnu.org/licenses/>. ! !$Id$ !> PDAF-OMI observation operators !! !! This module contains generic routines for several observation !! operators to be used after preparation with init_dim_obs_f !! !! The observation operators are: !! !! * PDAFomi_obs_op_gridpoint\n !! Observation operator for data at grid points. The routine !! selects values of the state vector according to an index array !! * PDAFomi_obs_op_gridavg\n !! Observation operator for the case that the observations are the !! average of grid point values. The routine computes these !! averages according to an index array. !! * PDAFomi_obs_op_interp_lin\n !! Observation operator for the case that the observations are !! linear interpolated from the grid points. The interpolation !! coefficients are pre-computed. !! * PDAFomi_obs_op_gatheronly\n !! Observation operator for the case of strongly coupled assimilation !! to gather an observation which only exists in other compartments. !! !! Adjoint observation operators: !! !! * PDAFomi_obs_op_adj_gridpoint\n !! Adjoint observation operator for data at grid points. The routine !! selects values of the state vector according to an index array !! * PDAFomi_obs_op_adj_gridavg\n !! Adjoint observation operator for the case that the observations !! are the average of grid point values. The routine computes these !! averages according to an index array. !! * PDAFomi_obs_op_adj_interp_lin\n !! Adjoint observation operator for the case that the observations !! are linear interpolatied from the grid points. The interpolation !! coefficients are pre-computed. !! * PDAFomi_obs_op_adj_gatheronly\n !! Adjoint observation operator for the case of strongly coupled assimilation !! to gather an observation which only exists in other compartments. !! !! Helper routines for the operators: !! * PDAFomi_get_interp_coeff_tri \n !! Routine to compute interpolation coefficients for triangular !! interpolation from barycentric coordinates. !! * PDAFomi_get_interp_coeff_lin1D \n !! Routine to comput linear interpolation in 1D !! * PDAFomi_get_interp_coeff_lin \n !! Routine to compute interpolation coefficients for linear !! interpolations (linear, bi-linear, tri-linear) !! MODULE PDAFomi_obs_op USE PDAFomi_obs_f, ONLY: obs_f, PDAFomi_gather_obsstate, debug CONTAINS !------------------------------------------------------------------------------- !> observation operator for data at grid points !! !! Application of observation operator for the case that !! model variables are observerved at model grid points. !! !! For this case INIT_DIM_OBS_F will prepare the index !! array thisobs%id_obs_p containing the information which !! elements of the PE-local state vector contain the !! observed values. !! !! The routine is called by all filter processes. It first !! selects the observed elements for a PE-local domain. !! Afterwards, the values are gathered into the full vector !! using PDAFomi_gather_obsstate. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_gridpoint(thisobs, state_p, obs_f_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: state_p(:) !< PE-local model state (dim_p) REAL, INTENT(inout) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) ! *** Local variables *** INTEGER :: i ! Counter REAL, ALLOCATABLE :: ostate_p(:) ! local observed part of state vector ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* doassim: IF (thisobs%doassim == 1) THEN ! Consistency check IF (.NOT.ALLOCATED(thisobs%id_obs_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_gridpoint - thisobs%id_obs_p is not allocated' END IF ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gridpoint -- START' WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_obs_op_gridpoint -- Process-local selection' WRITE (*,*) '++ OMI-debug obs_op_gridpoint:', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug obs_op_gridpoint:', debug, 'thisobs%id_obs_p', thisobs%id_obs_p END IF ! *** PE-local: Initialize observed part state vector IF (thisobs%dim_obs_p>0) THEN ALLOCATE(ostate_p(thisobs%dim_obs_p)) ELSE ALLOCATE(ostate_p(1)) END IF DO i = 1, thisobs%dim_obs_p ostate_p(i) = state_p(thisobs%id_obs_p(1, i)) ENDDO ! *** Global: Gather full observed state vector CALL PDAFomi_gather_obsstate(thisobs, ostate_p, obs_f_all) ! *** Clean up DEALLOCATE(ostate_p) ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gridpoint -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_gridpoint !------------------------------------------------------------------------------- !> Observation operator for averaging grid point values !! !! Application of observation operator for the case that !! the observation value is given as the average of model !! grid point values. !! For this case INIT_DIM_OBS_F will prepare the index !! array thisobs%id_obs_p that contains several rows holding !! the indices of the state vector elements which are to !! be averaged to represent an observation. !! !! The routine is called by all filter processes, !! and the operation has to be performed by each !! these processes for its PE-local domain before the !! information from all PEs is gathered. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_gridavg(thisobs, nrows, state_p, obs_f_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nrows !< Number of values to be averaged REAL, INTENT(in) :: state_p(:) !< PE-local model state (dim_p) REAL, INTENT(inout) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) ! *** Local variables *** INTEGER :: i, row ! Counter REAL, ALLOCATABLE :: ostate_p(:) ! local observed part of state vector REAL :: rrows ! Real-value for nrows ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gridavg -- START' WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_obs_op_gridavg -- Process-local averaging' WRITE (*,*) '++ OMI-debug obs_op_gridavg:', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug obs_op_gridavg:', debug, 'number of points to average', nrows WRITE (*,*) '++ OMI-debug obs_op_gridavg:', debug, 'thisobs%id_obs_p', thisobs%id_obs_p END IF ! Consistency check IF (.NOT.ALLOCATED(thisobs%id_obs_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_gridavg - thisobs%id_obs_p is not allocated' END IF ! *** PE-local: Initialize observed part state vector by averaging IF (thisobs%dim_obs_p>0) THEN ALLOCATE(ostate_p(thisobs%dim_obs_p)) ELSE ALLOCATE(ostate_p(1)) END IF rrows = REAL(nrows) DO i = 1, thisobs%dim_obs_p ostate_p(i) = 0.0 DO row = 1, nrows ostate_p(i) = ostate_p(i) + state_p(thisobs%id_obs_p(row,i)) END DO ostate_p(i) = ostate_p(i) / rrows ENDDO ! *** Global: Gather full observed state vector CALL PDAFomi_gather_obsstate(thisobs, ostate_p, obs_f_all) ! *** Clean up DEALLOCATE(ostate_p) ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gridavg -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_gridavg !------------------------------------------------------------------------------- !> Observation operator for linear interpolation !! !! Application of observation operator for the case that the !! observation value is given as the interpolation using !! pre-computed coefficients. !! !! For this case INIT_DIM_OBS_F will prepare the index !! array thisobs%id_obs_p containing the information which !! elements of the PE-local state vector contain the !! observed values. Further the array thisobs%icoeff_p is !! prepared which contains the interpolation coefficients. !! This can be prepared using a help routine like !! get_interp_coeff_tri. !! !! The routine is called by all filter processes. It first !! selects the observed elements for a PE-local domain. !! Afterwards, the values are gathered into the full vector !! using PDAFomi_gather_obsstate. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2019-12 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_interp_lin(thisobs, nrows, state_p, obs_f_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nrows !< Number of values to be averaged REAL, INTENT(in) :: state_p(:) !< PE-local model state (dim_p) REAL, INTENT(inout) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) ! *** Local variables *** INTEGER :: i, row ! Counters REAL, ALLOCATABLE :: ostate_p(:) ! local observed part of state vector ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_interp_lin -- START' WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_obs_op_interp_lin -- Process-local interpolation' WRITE (*,*) '++ OMI-debug obs_op_interp_lin:', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug obs_op_interp_lin:', debug, 'number of points in interpolation', nrows WRITE (*,*) '++ OMI-debug obs_op_interp_lin:', debug, 'thisobs%id_obs_p', thisobs%id_obs_p WRITE (*,*) '++ OMI-debug obs_op_interp_lin:', debug, 'thisobs%icoeff_p', thisobs%icoeff_p END IF ! Check if required arrays are allocated (assuming that they are initialzed in this case) IF (.NOT.ALLOCATED(thisobs%id_obs_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_interp_lin - thisobs%id_obs_p is not allocated' END IF IF (.NOT.ALLOCATED(thisobs%icoeff_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_interp_lin - thisobs%icoeff_p is not allocated' END IF ! *** PE-local: Initialize observed part state vector by weighted averaging IF (thisobs%dim_obs_p>0) THEN ALLOCATE(ostate_p(thisobs%dim_obs_p)) ELSE ALLOCATE(ostate_p(1)) END IF DO i = 1, thisobs%dim_obs_p ostate_p(i) = 0.0 DO row = 1, nrows ostate_p(i) = ostate_p(i) + thisobs%icoeff_p(row,i)*state_p(thisobs%id_obs_p(row,i)) END DO ENDDO ! *** Global: Gather full observed state vector CALL PDAFomi_gather_obsstate(thisobs, ostate_p, obs_f_all) ! *** Clean up DEALLOCATE(ostate_p) ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_interp_lin -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_interp_lin !------------------------------------------------------------------------------- !> observation operator for the case that observations belong to other compartment !! !! Application of observation operator for the case that !! model variables are observerved in another compartment !! only. Thus DIM_OBS_P of the current compartment is 0. !! Accordingly, this observation operator only performs !! the gather operation to obtain the full observations. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2020-04 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_gatheronly(thisobs, state_p, obs_f_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: state_p(:) !< PE-local model state (dim_p) REAL, INTENT(inout) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) ! *** Local variables *** REAL, ALLOCATABLE :: ostate_p(:) ! local observed part of state vector REAL :: rdummy ! dummy variable to prevent compiler warning ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* ! Initialize dummy to prevent compiler warning rdummy = state_p(1) doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gatheronly -- START' ! *** PE-local: Nothing to be done! ALLOCATE(ostate_p(1)) ostate_p = 0.0 ! *** Global: Gather full observed state vector CALL PDAFomi_gather_obsstate(thisobs, ostate_p, obs_f_all) ! *** Clean up DEALLOCATE(ostate_p) ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gatheronly -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_gatheronly !------------------------------------------------------------------------------- !> Helper routine: Initialize interpolation coefficients in triangle !! !! The routine computes the coefficients for triangular interpolation !! as barycentric coordinates. !! The computation is done for one observation given the !! observation coordinates (OC) as well as the coordinates of the !! grid points (GPC). In GPC each row contains the coordinates !! for one grid point. Thus the first index determines the grid point, !! while the second the coordinates of this grid point !! !! __Revision history:__ !! * 2019-12 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_get_interp_coeff_tri(gpc, oc, icoeff) IMPLICIT NONE ! *** Arguments *** REAL, INTENT(in) :: gpc(:,:) !< Coordinates of grid points; dim(3,2) ! 3 rows; each containing lon and lat coordinates REAL, INTENT(in) :: oc(:) !< Coordinates of observation; dim(2) REAL, INTENT(inout) :: icoeff(:) !< Interpolation coefficients; dim(3) ! *** Local variables *** REAL :: denum ! denumerator ! ****************************************** ! *** Compute interpolation coefficients *** ! *** as barycentric coordinates *** ! ****************************************** ! common denumerator for coefficients 1 and 2 denum = (gpc(2,2) - gpc(3,2)) * (gpc(1,1) - gpc(3,1)) + (gpc(3,1) - gpc(2,1)) * (gpc(1,2) - gpc(3,2)) ! compute coefficients icoeff(1) = (gpc(2,2) - gpc(3,2)) * (oc(1) - gpc(3,1)) + (gpc(3,1) - gpc(2,1)) * (oc(2) - gpc(3,2)) icoeff(1) = icoeff(1) / denum icoeff(2) = (gpc(3,2) - gpc(1,2)) * (oc(1) - gpc(3,1)) + (gpc(1,1) - gpc(3,1)) * (oc(2) - gpc(3,2)) icoeff(2) = icoeff(2) / denum icoeff(3) = 1.0 - icoeff(1) - icoeff(2) END SUBROUTINE PDAFomi_get_interp_coeff_tri !------------------------------------------------------------------------------- !> Helper routine: Initialize linear interpolation coefficients in 1D !! !! The routine computes the coefficients for linear interpolation !! in 1 dimensions. !! The computation is done for one observation given the !! observation coordinates (OC) as well as the coordinates of the !! grid points (GPC). !! !! __Revision history:__ !! * 2019-12 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_get_interp_coeff_lin1D(gpc, oc, icoeff) IMPLICIT NONE ! *** Arguments *** REAL, INTENT(in) :: gpc(:) !< Coordinates of grid points (dim=2) REAL, INTENT(in) :: oc !< Coordinates of observation REAL, INTENT(inout) :: icoeff(:) !< Interpolation coefficients (dim=2) ! **************************************************************** ! *** Compute linear interpolation coefficients in 1 dimension *** ! **************************************************************** icoeff(1) = (gpc(2) - oc) / (gpc(2) - gpc(1)) icoeff(2) = (oc - gpc(1)) / (gpc(2) - gpc(1)) END SUBROUTINE PDAFomi_get_interp_coeff_lin1D !------------------------------------------------------------------------------- !> Helper routine: Initialize linear interpolation coefficients !! !! The routine computes the coefficients for linear interpolation !! in 1, 2, or 3 dimensions. !! The computation is done for one observation given the !! observation coordinates (OC) as well as the coordinates of the !! grid points (GPC). In GPC each row contains the coordinates !! for one grid point. !! !! Setup of GPC: !! The first index is specifies the grid point, while the second the !! coordinate !! * For n_dim=X only the first X coordinate values are used !! * The ordering of the coordinates and coefficient is the following: !! !! (7)------(8) !! /| /| with !! (5)+-----(6)| - column 1 !! | | | | / column 2 !! |(3)-----+(4) | column 3 !! |/ |/ !! (1) ---- (2) !! !! thus gpc(1,1)/=gpc(2,1), gpc(1,2)/=gpc(3,2), gpc(1,3)/=gpc(5,3) !! but gpc(1,1)=gpc(3,1)=gpc(5,1), gpc(1,2)=gpc(2,2)=gpc(5,2), !! gpc(1,3)=gpc(2,3)=gpc(3,3) !! * For bi-linear interpolation only the coordinates for grid !! points 1, 2, and 3 are used to compute the coefficients !! * For tri-linear interpolation only the coordinates for grid !! points 1, 2, 3, and 5 are used to compute the coefficients !! * (for bi-linear interpolation gpc only needs to have length 3 !! for tri-linear the length 5) !! * num_gp=2 for n_dim=1; num_gp=4 for n_dim=2; num_gp=8 for n_dim=3 !! is required !! !! __Revision history:__ !! * 2019-12 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_get_interp_coeff_lin(num_gp, n_dim, gpc, oc, icoeff) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: num_gp !< Length of icoeff INTEGER, INTENT(in) :: n_dim !< Number of dimensions in interpolation REAL, INTENT(in) :: gpc(:,:) !< Coordinates of grid points REAL, INTENT(in) :: oc(:) !< Coordinates of observation REAL, INTENT(inout) :: icoeff(:) !< Interpolation coefficients (num_gp) ! *** Local variables *** REAL :: denum ! denumerator IF (n_dim == 1) THEN ! **************************************************************** ! *** Compute linear interpolation coefficients in 1 dimension *** ! **************************************************************** ! Checks IF (num_gp /= 2) WRITE (*,'(a,3x,a)') & 'PDAFomi', 'ERROR: get_interp_coeff_lin - NUM_GP=2 required!' IF (gpc(2,1) == gpc(1,1)) WRITE (*,'(a,3x,a)') & 'PDAFomi', 'ERROR: get_interp_coeff_lin - wrong setting of coordinates!' ! Compute coefficients icoeff(1) = (gpc(2,1) - oc(1)) / (gpc(2,1) - gpc(1,1)) icoeff(2) = (oc(1) - gpc(1,1)) / (gpc(2,1) - gpc(1,1)) ELSE IF (n_dim == 2) THEN ! ******************************************************** ! *** Compute coefficients for bi-linear interpolation *** ! *** Order of coefficients: (3) ---- (4) *** ! *** | | *** ! *** (1) ---- (2) *** ! ******************************************************** ! Checks IF (num_gp /= 4) WRITE (*,'(a,5x,a)') & 'PDAFomi', 'ERROR: get_interp_coeff_lin - NUM_GP=4 required!' IF (gpc(2,1) == gpc(1,1) .OR. gpc(3,2) == gpc(1,2)) WRITE (*,'(a,3x,a)') & 'PDAFomi', 'ERROR: get_interp_coeff_lin - wrong setting of coordinates!' ! Compute coefficients denum = (gpc(2,1) - gpc(1,1)) * (gpc(3,2) - gpc(1,2)) icoeff(1) = (gpc(2,1) - oc(1)) * (gpc(3,2) - oc(2)) / denum icoeff(2) = (oc(1) - gpc(1,1)) * (gpc(3,2) - oc(2)) / denum icoeff(3) = (gpc(2,1) - oc(1)) * (oc(2) - gpc(1,2)) / denum icoeff(4) = (oc(1) - gpc(1,1)) * (oc(2) - gpc(1,2)) / denum ELSE IF (n_dim == 3) THEN ! ********************************************************* ! *** Compute coefficients for tri-linear interpolation *** ! *** Order of coefficients: (7)------(8) *** ! *** /| /| *** ! *** (5)+-----(6)| *** ! *** | | | | *** ! *** |(3)-----+(4) *** ! *** |/ |/ *** ! *** (1) ---- (2) *** ! ********************************************************* ! Checks IF (num_gp /= 8) WRITE (*,'(a,5x,a)') & 'PDAFomi', 'ERROR: get_interp_coeff_lin - NUM_GP=8 required!' IF (gpc(2,1) == gpc(1,1) .OR. gpc(3,2) == gpc(1,2) .OR. gpc(5,3) == gpc(1,3)) & WRITE (*,'(a,3x,a)') & 'PDAFomi', 'ERROR: get_interp_coeff_lin - wrong setting of coordinates!' ! Compute coefficients denum = (gpc(2,1) - gpc(1,1)) * (gpc(3,2) - gpc(1,2)) * (gpc(5,3) - gpc(1,3)) icoeff(1) = (gpc(2,1) - oc(1)) * (gpc(3,2) - oc(2)) * (gpc(5,3) - oc(3)) / denum icoeff(2) = (oc(1) - gpc(1,1)) * (gpc(3,2) - oc(2)) * (gpc(5,3) - oc(3)) / denum icoeff(3) = (gpc(2,1) - oc(1)) * (oc(2) - gpc(1,2)) * (gpc(5,3) - oc(3)) / denum icoeff(4) = (oc(1) - gpc(1,1)) * (oc(2) - gpc(1,2)) * (gpc(5,3) - oc(3)) / denum icoeff(5) = (gpc(2,1) - oc(1)) * (gpc(3,2) - oc(2)) * (oc(3) - gpc(1,3)) / denum icoeff(6) = (oc(1) - gpc(1,1)) * (gpc(3,2) - oc(2)) * (oc(3) - gpc(1,3)) / denum icoeff(7) = (gpc(2,1) - oc(1)) * (oc(2) - gpc(1,2)) * (oc(3) - gpc(1,3)) / denum icoeff(8) = (oc(1) - gpc(1,1)) * (oc(2) - gpc(1,2)) * (oc(3) - gpc(1,3)) / denum END IF END SUBROUTINE PDAFomi_get_interp_coeff_lin !------------------------------------------------------------------------------- !> Adjoint observation operator for data at grid points !! !! Application of adjoint observation operator for the case !! that model variables are observerved at model grid points. !! !! For this case INIT_DIM_OBS_F will prepare the index !! array thisobs%id_obs_p containing the information which !! elements of the PE-local state vector contain the !! observed values. !! !! The routine is called by all filter processes. It first !! selects the observed elements for a PE-local domain. !! Afterwards, the values are gathered into the full vector !! using PDAFomi_gather_obsstate. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2021-04 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_adj_gridpoint(thisobs, obs_f_all, state_p) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) REAL, INTENT(inout) :: state_p(:) !< PE-local model state (dim_p) ! *** Local variables *** INTEGER :: i ! Counter ! ************************************************** ! *** Perform application of adjoint observation *** ! *** operator H^T on vector or matrix column *** ! ************************************************** doassim: IF (thisobs%doassim == 1) THEN ! Consistency check IF (.NOT.ALLOCATED(thisobs%id_obs_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_adj_gridpoint - thisobs%id_obs_p is not allocated' END IF ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_adj_gridpoint -- START' WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_obs_op_adj_gridpoint -- Process-local selection' WRITE (*,*) '++ OMI-debug obs_op_adj_gridpoint:', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug obs_op_adj_gridpoint:', debug, 'thisobs%id_obs_p', thisobs%id_obs_p WRITE (*,*) '++ OMI-debug obs_op_adj_gridpoint:', debug, 'thisobs%off_obs_f', thisobs%off_obs_f END IF ! *** PE-local: Apply adjoint observation operator DO i = 1, thisobs%dim_obs_p state_p(thisobs%id_obs_p(1, i)) & = state_p(thisobs%id_obs_p(1, i)) + obs_f_all(thisobs%off_obs_f+i) ENDDO ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_adj_gridpoint -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_adj_gridpoint !------------------------------------------------------------------------------- !> Adjoint observation operator for averaging grid point values !! !! Application of adjoint observation operator for the case !! that the observation value is given as the average of !! model grid point values. !! !! For this case INIT_DIM_OBS_F will prepare the index !! array thisobs%id_obs_p that contains several rows holding !! the indices of the state vector elements which are to !! be averaged to represent an observation. !! !! The routine is called by all filter processes, !! and the operation has to be performed by each !! these processes for its PE-local domain before the !! information from all PEs is gathered. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2021-12 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_adj_gridavg(thisobs, nrows, obs_f_all, state_p) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nrows !< Number of values to be averaged REAL, INTENT(in) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) REAL, INTENT(inout) :: state_p(:) !< PE-local model state (dim_p) ! *** Local variables *** INTEGER :: i, row ! Counter REAL :: rrows ! Real-value for nrows ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_adj_gridavg -- START' WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_obs_op_adj_gridavg -- Process-local averaging' WRITE (*,*) '++ OMI-debug obs_op_adj_gridpoint:', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug obs_op_adj_gridpoint:', debug, 'number of points to average', nrows WRITE (*,*) '++ OMI-debug obs_op_adj_gridpoint:', debug, 'thisobs%id_obs_p', thisobs%id_obs_p END IF ! Consistency check IF (.NOT.ALLOCATED(thisobs%id_obs_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_adj_gridavg - thisobs%id_obs_p is not allocated' END IF ! *** PE-local: Initialize observed part state vector by averaging rrows = REAL(nrows) DO i = 1, thisobs%dim_obs_p DO row = 1, nrows state_p(thisobs%id_obs_p(row, i)) & = state_p(thisobs%id_obs_p(row, i)) + obs_f_all(thisobs%off_obs_f+i) / rrows end DO ENDDO ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_adj_gridavg -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_adj_gridavg !------------------------------------------------------------------------------- !> Observation operator for linear interpolation !! !! Application of adjoint observation operator for the case !! that the observation value is given as the interpolation !! using pre-computed coefficients. !! !! For this case INIT_DIM_OBS_F will prepare the index !! array thisobs%id_obs_p containing the information which !! elements of the PE-local state vector contain the !! observed values. Further the array thisobs%icoeff_p is !! prepared which contains the interpolation coefficients. !! This can be prepared using a help routine like !! get_interp_coeff_tri. !! !! The routine is called by all filter processes. It first !! selects the observed elements for a PE-local domain. !! Afterwards, the values are gathered into the full vector !! using PDAFomi_gather_obsstate. !! !! The routine has to fill the part of the full observation !! vector OBS_F_ALL that represents the current observation !! type. The routine first applied the observation operator !! for the current observation type and the calls !! PDAFomi_gather_obsstate to gather the observation over !! all processes and fills OBS_F_ALL. !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2021-12 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_adj_interp_lin(thisobs, nrows, obs_f_all, state_p) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nrows !< Number of values to be averaged REAL, INTENT(in) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) REAL, INTENT(inout) :: state_p(:) !< PE-local model state (dim_p) ! *** Local variables *** INTEGER :: i, row ! Counters ! ************************************************** ! *** Perform application of adjoint observation *** ! *** operator H^T on vector or matrix column *** ! ************************************************** doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_adj_interp_lin -- START' WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_obs_op_adj_interp_lin -- Process-local interpolation' WRITE (*,*) '++ OMI-debug obs_op_adj_interp_lin:', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug obs_op_adj_interp_lin:', debug, 'number of points in interpolation', nrows WRITE (*,*) '++ OMI-debug obs_op_adj_interp_lin:', debug, 'thisobs%id_obs_p', thisobs%id_obs_p WRITE (*,*) '++ OMI-debug obs_op_adj_interp_lin:', debug, 'thisobs%icoeff_p', thisobs%icoeff_p END IF ! Check if required arrays are allocated (assuming that they are initialzed in this case) IF (.NOT.ALLOCATED(thisobs%id_obs_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_adj_interp_lin - thisobs%id_obs_p is not allocated' END IF IF (.NOT.ALLOCATED(thisobs%icoeff_p)) THEN WRITE (*,*) 'ERROR: PDAFomi_obs_op_adj_interp_lin - thisobs%icoeff_p is not allocated' END IF ! *** PE-local: Apply adjoint observation operator DO i = 1, thisobs%dim_obs_p DO row = 1, nrows state_p(thisobs%id_obs_p(row, i)) & = state_p(thisobs%id_obs_p(row, i)) + thisobs%icoeff_p(row,i)*obs_f_all(thisobs%off_obs_f+i) end DO ENDDO ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_adj_interp_lin -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_adj_interp_lin !------------------------------------------------------------------------------- !> adjoint observation operator for the case that observations belong !! to other compartment !! !! Application of adjoint observation operator for the case that !! model variables are observerved in another compartment !! only. Thus DIM_OBS_P of the current compartment is 0. !! Accordingly, this observation operator only performs !! the gather operation to obtain the full observations. !! !! The routine does nothing, since there are no observations !! on the compartment !! !! The routine has to be called by all filter processes. !! !! __Revision history:__ !! * 2021-12 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obs_op_adj_gatheronly(thisobs, obs_f_all, state_p) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: state_p(:) !< PE-local model state (dim_p) REAL, INTENT(inout) :: obs_f_all(:) !< Full observed state for all observation types (nobs_f_all) ! *** Local variables *** ! ********************************************* ! *** Perform application of measurement *** ! *** operator H on vector or matrix column *** ! ********************************************* doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gatheronly -- START' ! *** Nothing to be done! ! Print debug information IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_obs_op_gatheronly -- END' END IF doassim END SUBROUTINE PDAFomi_obs_op_adj_gatheronly END MODULE PDAFomi_obs_op