! 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: PDAFomi_obs_f.F90 1147 2023-03-12 16:14:34Z lnerger $ !> PDAF-OMI routines for full observations !! !! This module contains subroutines to handle full observations. !! Further, it contains routines to restrict the global full vector of observations !! to those observations that are relevant for a process-local model subdomain. !! The routines are !! !! * PDAFomi_gather_obs \n !! Gather full observation information !! * PDAFomi_gather_obsstate \n !! Gather a full observed state vector (used in observation operators) !! * PDAFomi_init_obs_f \n !! Initialize full vector of observations for adaptive forgetting factor !! * PDAFomi_init_obsvar_f \n !! Compute mean observation error variance for adaptive forgetting factor !! * PDAFomi_prodRinvA \n !! Multiply an intermediate matrix of the global filter analysis !! with the inverse of the observation error covariance matrix !! * PDAFomi_likelihood \n !! Compute likelihood for an ensemble member !! * PDAFomi_add_obs_err \n !! Add observation error to some matrix !! * PDAFomi_init_obscovar \n !! Initialize global observation error covariance matrix !! * PDAFomi_set_domain_limits \n !! Set min/max coordinate locations of decomposed grid !! * PDAFomi_get_domain_limits_unstr \n !! Find min/max coordinate locations in unstructured grid !! * PDAFomi_get_local_ids_obs_f \n !! Find observations inside or close to process domain !! * PDAFomi_limit_obs_f \n !! Reduce full observation vector to part relevant for local process domain !! !! The routine PDAFomi_get_domain_limits_unstr assumed geographic coordinates in radians !! and within the range -pi to +pi for longitude (- is westward) and -pi/2 to +pi/2 for !! latitude. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code !! * Later revisions - see repository log !! MODULE PDAFomi_obs_f USE mpi USE PDAF_mod_filtermpi, & ONLY: mype, npes, COMM_FILTER, MPIerr USE PDAF_mod_filter, & ONLY: screen, obs_member, filterstr, dim_p IMPLICIT NONE SAVE ! *** Module internal variables INTEGER :: debug=0 !< Debugging flag INTEGER :: error=0 !< Error flag REAL, ALLOCATABLE :: domain_limits(:) !< Limiting coordinates (NSWE) for process domain REAL, PARAMETER :: r_earth=6.3675e6 !< Earth radius in meters REAL, PARAMETER :: pi=3.141592653589793 !< Pi ! *** Data type to define the full observations by internally shared variables of the module TYPE obs_f ! ---- Mandatory variables to be set in INIT_DIM_OBS ---- INTEGER :: doassim=0 !< Whether to assimilate this observation type INTEGER :: disttype !< Type of distance computation to use for localization !< (0) Cartesian, (1) Cartesian periodic !< (2) simplified geographic, (3) geographic haversine function !< (10,11,12,13) factorized 2+1D localization with distance !< calculation from (0)-(3); obs. weighting is only done with !< horizontal distance, which vertical uses only cut-off radius INTEGER :: ncoord !< Number of coordinates use for distance computation INTEGER, ALLOCATABLE :: id_obs_p(:,:) !< Indices of process-local observed field in state vector ! ---- Optional variables - they can be set in INIT_DIM_OBS ---- REAL, ALLOCATABLE :: icoeff_p(:,:) !< Interpolation coefficients for obs. operator (optional) REAL, ALLOCATABLE :: domainsize(:) !< Size of domain for periodicity (<=0 for no periodicity) (optional) ! ---- Variables with predefined values - they can be changed in INIT_DIM_OBS ---- INTEGER :: obs_err_type=0 !< Type of observation error: (0) Gauss, (1) Laplace INTEGER :: use_global_obs=1 !< Whether to use (1) global full obs. !< or (0) obs. restricted to those relevant for a process domain REAL :: inno_omit=0.0 !< Omit obs. if squared innovation larger this factor times !< observation variance (only active for >0) REAL :: inno_omit_ivar=1.0e-12 !< Value of inverse variance to omit observation !< (should be much larger than actual observation error variance) ! ---- The following variables are set in the routine PDAFomi_gather_obs --- INTEGER :: dim_obs_p !< number of PE-local observations INTEGER :: dim_obs_f !< number of full observations INTEGER :: dim_obs_g !< global number of observations INTEGER :: off_obs_f !< Offset of this observation in overall full obs. vector INTEGER :: off_obs_g !< Offset of this observation in overall global obs. vector INTEGER :: obsid !< Index of observation over all assimilated observations REAL, ALLOCATABLE :: obs_f(:) !< Full observed field REAL, ALLOCATABLE :: ocoord_f(:,:) !< Coordinates of full observation vector REAL, ALLOCATABLE :: ivar_obs_f(:) !< Inverse variance of full observations INTEGER, ALLOCATABLE :: id_obs_f_lim(:) !< Indices of domain-relevant full obs. in global vector of obs. ! ---- Other internal variables --- INTEGER :: locweight_v !< Type of localization function in vertical direction (for disttype>=10) END TYPE obs_f INTEGER :: n_obstypes = 0 ! Number of observation types INTEGER :: obscnt = 0 ! current ID of observation type INTEGER :: offset_obs = 0 ! offset of current observation in overall observation vector INTEGER :: offset_obs_g = 0 ! offset of current observation in global observation vector LOGICAL :: omit_obs = .FALSE. ! Flag whether observations are omitted for large innovation LOGICAL :: omi_was_used = .FALSE. ! Flag whether OMI was used INTEGER, ALLOCATABLE :: obsdims(:,:) ! Observation dimensions over all types and process sub-domains INTEGER, ALLOCATABLE :: map_obs_id(:) ! Index array to map obstype-first index to domain-first index INTEGER :: ostats_omit(7) ! PE-local statistics ! ostats_omit(1): Number of local domains with excluded observations ! ostats_omit(2): Number of local domains without excluded observations ! ostats_omit(3): Sum of all excluded observations for all domains ! ostats_omit(4): Sum of all used observations for all domains ! ostats_omit(5): Sum of all used observations for all domains ! ostats_omit(6): Maximum count of excluded observations over all domains ! ostats_omit(7): Maximum count of used observations over all domains TYPE obs_arr_f TYPE(obs_f), POINTER :: ptr END TYPE obs_arr_f TYPE(obs_arr_f), ALLOCATABLE :: obs_f_all(:) !$OMP THREADPRIVATE(debug) !------------------------------------------------------------------------------- CONTAINS !> Gather full observational information !! !! This routine uses different gather routines from PDAF to obtain !! combined full observational information. !! !! The observation-type specific variables that are initialized here are !! * thisobs\%dim_obs_p - PE-local number of module-type observations !! * thisobs\%dim_obs_f - full number of module-type observations !! * thisobs\%off_obs_f - Offset of full module-type observation in overall full obs. vector !! * thisobs\%off_obs_g - Offset of global module-type observation in overall full obs. vector !! * thisobs\%obs_f - full vector of module-type observations !! * thisobs\%ocoord_f - coordinates of observations in OBS_MOD_F !! * thisobs\%ivar_obs_f - full vector of inverse obs. error variances of module-type !! !! If the full observations are restricted to those relevant for a process domain, !! there are further initialized !! * thisobs\%dim_obs_g - Number of global observations !! * thisobs\%id_obs_f_lim - Ids of full observations in global observations !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_gather_obs(thisobs, dim_obs_p, obs_p, ivar_obs_p, ocoord_p, & ncoord, lradius, dim_obs_f) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: dim_obs_p !< Number of process-local observation REAL, INTENT(in) :: obs_p(:) !< Vector of process-local observations REAL, INTENT(in) :: ivar_obs_p(:) !< Vector of process-local inverse observation error variance REAL, INTENT(in) :: ocoord_p(:,:) !< Array of process-local observation coordinates INTEGER, INTENT(in) :: ncoord !< Number of rows of coordinate array REAL, INTENT(in) :: lradius !< Localization radius (the maximum radius used in this process domain) INTEGER, INTENT(out) :: dim_obs_f !< Full number of observations ! *** Local variables *** INTEGER :: i ! Counter REAL, ALLOCATABLE :: obs_g(:) ! Global full observation vector (used in case of limited obs.) REAL, ALLOCATABLE :: ivar_obs_g(:) ! Global full inverse variances (used in case of limited obs.) REAL, ALLOCATABLE :: ocoord_g(:,:) ! Global full observation coordinates (used in case of limited obs.) INTEGER :: status ! Status flag for PDAF gather operation INTEGER :: localfilter ! Whether the filter is domain-localized INTEGER :: globalobs ! Whether the filter needs global observations INTEGER :: maxid ! maximum index in thisobs%id_obs_p ! ********************** ! *** Initialization *** ! ********************** ! Increment counter of observation types n_obstypes = n_obstypes + 1 ! Set observation ID thisobs%obsid = n_obstypes ! Set flag that OMI was used omi_was_used = .TRUE. IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 5x, a, 1x, i3)') 'PDAFomi', '--- Initialize observation type ID', thisobs%obsid ! ************************************** ! *** Gather full observation arrays *** ! ************************************** ! Consistency check maxid = MAXVAL(thisobs%id_obs_p) IF (maxid > dim_p .AND. dim_obs_p>0) THEN ! Maximum value of id_obs_p point to outside of state vector WRITE (*,'(a)') 'PDAFomi - ERROR: thisobs%id_obs_p too large - index points to outside of state vector !!!' error = 1 END IF ! Check whether the filter is domain-localized CALL PDAF_get_localfilter(localfilter) ! Check whether the filter needs global observations CALL PDAF_get_globalobs(globalobs) ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_gather_obs -- START Gather full observation vector' IF (localfilter==1) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'domain localized filter' ELSE WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'filter without domain-localization' END IF IF (globalobs==1) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'filter uses global observations' END IF END IF lfilter: IF (localfilter==1 .OR. globalobs==1) THEN ! For domain-localized filters: gather full observations fullobs: IF (thisobs%use_global_obs==1 .OR. globalobs==1) THEN ! *** Use global full observations *** IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 5x, a)') 'PDAFomi', '--- Use global full observations' ! *** Initialize global dimension of observation vector *** CALL PDAF_gather_dim_obs_f(dim_obs_p, dim_obs_f) ! Store full and PE-local observation dimensions in module variables thisobs%dim_obs_p = dim_obs_p thisobs%dim_obs_f = dim_obs_f IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 8x, a, i7)') 'PDAFomi', & '--- Number of full observations ', dim_obs_f IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_f', thisobs%dim_obs_f WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'obs_p', obs_p WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'ocoord_p', ocoord_p WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%disttype', thisobs%disttype END IF ! *** Gather full observation vector and corresponding coordinates *** ! Allocate full observation arrays ! The arrays are deallocated in PDAFomi_deallocate_obs in PDAFomi_obs_l IF (dim_obs_f > 0) THEN ALLOCATE(thisobs%obs_f(dim_obs_f)) ALLOCATE(thisobs%ivar_obs_f(dim_obs_f)) ALLOCATE(thisobs%ocoord_f(ncoord, dim_obs_f)) ELSE ALLOCATE(thisobs%obs_f(1)) ALLOCATE(thisobs%ivar_obs_f(1)) ALLOCATE(thisobs%ocoord_f(ncoord, 1)) END IF CALL PDAFomi_gather_obs_f_flex(dim_obs_p, obs_p, thisobs%obs_f, status) CALL PDAFomi_gather_obs_f_flex(dim_obs_p, ivar_obs_p, thisobs%ivar_obs_f, status) CALL PDAFomi_gather_obs_f2_flex(dim_obs_p, ocoord_p, thisobs%ocoord_f, ncoord, status) ELSE fullobs ! *** Use full observations limited to those relevant for a process domain *** ! *** This can be more efficient as in the local analysis loop less *** ! *** observations have a be checked for each analysis domain *** IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 5x, a)') 'PDAFomi', '--- Use limited full observations' ! *** Initialize global dimension of observation vector *** CALL PDAF_gather_dim_obs_f(dim_obs_p, thisobs%dim_obs_g) ! *** First gather global observation vector and corresponding coordinates *** ! Allocate global observation arrays IF (thisobs%dim_obs_g > 0) THEN ALLOCATE(obs_g(thisobs%dim_obs_g)) ALLOCATE(ivar_obs_g(thisobs%dim_obs_g)) ALLOCATE(ocoord_g(ncoord, thisobs%dim_obs_g)) ELSE ALLOCATE(obs_g(1)) ALLOCATE(ivar_obs_g(1)) ALLOCATE(ocoord_g(ncoord, 1)) END IF CALL PDAFomi_gather_obs_f_flex(dim_obs_p, obs_p, obs_g, status) CALL PDAFomi_gather_obs_f_flex(dim_obs_p, ivar_obs_p, ivar_obs_g, status) CALL PDAFomi_gather_obs_f2_flex(dim_obs_p, ocoord_p, ocoord_g, ncoord, status) ! *** Now restrict the global observation arrays to the process-relevant parts *** ! Get number of full observation relevant for the process domain ! and corresponding indices in global observation vector IF (thisobs%dim_obs_g > 0) THEN ALLOCATE(thisobs%id_obs_f_lim(thisobs%dim_obs_g)) ELSE ALLOCATE(thisobs%id_obs_f_lim(1)) END IF IF (ALLOCATED(thisobs%domainsize)) THEN CALL PDAFomi_get_local_ids_obs_f(thisobs%dim_obs_g, lradius, ocoord_g, dim_obs_f, & thisobs%id_obs_f_lim, thisobs%disttype, thisobs%domainsize) ELSE CALL PDAFomi_get_local_ids_obs_f(thisobs%dim_obs_g, lradius, ocoord_g, dim_obs_f, & thisobs%id_obs_f_lim, thisobs%disttype) END IF ! Store full and PE-local observation dimensions in module variables thisobs%dim_obs_p = dim_obs_p thisobs%dim_obs_f = dim_obs_f IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_f', thisobs%dim_obs_f WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'obs_p', obs_p WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'ocoord_p', ocoord_p END IF ! Allocate global observation arrays ! The arrays are deallocated in PDAFomi_deallocate_obs in PDAFomi_obs_l IF (dim_obs_f > 0) THEN ALLOCATE(thisobs%obs_f(dim_obs_f)) ALLOCATE(thisobs%ivar_obs_f(dim_obs_f)) ALLOCATE(thisobs%ocoord_f(ncoord, dim_obs_f)) ! Get process-relevant full observation arrays CALL PDAFomi_limit_obs_f(thisobs, 0, obs_g, thisobs%obs_f) CALL PDAFomi_limit_obs_f(thisobs, 0, ivar_obs_g, thisobs%ivar_obs_f) DO i = 1, ncoord CALL PDAFomi_limit_obs_f(thisobs, 0, ocoord_g(i,:), thisobs%ocoord_f(i,:)) END DO ELSE ALLOCATE(thisobs%obs_f(1)) ALLOCATE(thisobs%ivar_obs_f(1)) ALLOCATE(thisobs%ocoord_f(ncoord, 1)) END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_gather_obs -- Limited full observations' WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_g', thisobs%dim_obs_g WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'obs_g', obs_g END IF DEALLOCATE(obs_g, ivar_obs_g, ocoord_g) END IF fullobs ELSE lfilter ! *** For global filters use process-local observations without gathering *** IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 5x, a)') 'PDAFomi', '--- Use process-local observations for global filters' ! *** Initialize global dimension of observation vector *** dim_obs_f = dim_obs_p ! *** Initialize global dimension of observation vector *** CALL PDAF_gather_dim_obs_f(dim_obs_p, thisobs%dim_obs_g) IF (TRIM(filterstr)=='ENKF' .OR. TRIM(filterstr)=='LENKF') THEN IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 8x, a, i7)') 'PDAFomi', & '--- Number of global observations ', thisobs%dim_obs_g ELSE IF (mype == 0 .AND. screen > 0) & WRITE (*, '(a, 8x, a, i7)') 'PDAFomi', & '--- Number of full observations ', dim_obs_f END IF ! *** Gather full observation vector and corresponding coordinates *** ! Allocate full observation arrays ! The arrays are deallocated in PDAFomi_deallocate_obs in PDAFomi_obs_l IF (dim_obs_f > 0) THEN ALLOCATE(thisobs%obs_f(dim_obs_f)) IF (TRIM(filterstr)=='ENKF' .OR. TRIM(filterstr)=='LENKF') THEN ! The LEnKF needs the global array ivar_obs_f ALLOCATE(thisobs%ivar_obs_f(thisobs%dim_obs_g)) ALLOCATE(thisobs%ocoord_f(ncoord, thisobs%dim_obs_g)) ELSE ALLOCATE(thisobs%ivar_obs_f(dim_obs_f)) ALLOCATE(thisobs%ocoord_f(ncoord, dim_obs_f)) END IF ELSE ALLOCATE(thisobs%obs_f(1)) IF (thisobs%dim_obs_g>0 .AND. (TRIM(filterstr)=='ENKF' .OR. TRIM(filterstr)=='LENKF')) THEN ! The LEnKF needs the global array ivar_obs_f ! Here dim_obs_f=0, but dim_obs_g>0 is possible in case of domain-decomposition IF (thisobs%dim_obs_g>0) THEN ALLOCATE(thisobs%ivar_obs_f(thisobs%dim_obs_g)) ALLOCATE(thisobs%ocoord_f(ncoord, thisobs%dim_obs_g)) ELSE ALLOCATE(thisobs%ivar_obs_f(1)) ALLOCATE(thisobs%ocoord_f(ncoord, 1)) END IF ELSE ALLOCATE(thisobs%ivar_obs_f(1)) ALLOCATE(thisobs%ocoord_f(ncoord, 1)) END IF END IF thisobs%obs_f = obs_p IF (TRIM(filterstr)=='ENKF' .OR. TRIM(filterstr)=='LENKF') THEN ! The EnKF and LEnKF need the global array ivar_obs_f CALL PDAFomi_gather_obs_f_flex(dim_obs_p, ivar_obs_p, & thisobs%ivar_obs_f, status) CALL PDAFomi_gather_obs_f2_flex(dim_obs_p, ocoord_p, & thisobs%ocoord_f, ncoord, status) ELSE thisobs%ivar_obs_f = ivar_obs_p END IF ! Store full and PE-local observation dimensions in module variables thisobs%dim_obs_p = dim_obs_p thisobs%dim_obs_f = dim_obs_f IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_f', thisobs%dim_obs_f IF (TRIM(filterstr)=='ENKF' .OR. TRIM(filterstr)=='LENKF') & WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%dim_obs_g', thisobs%dim_obs_g WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'obs_p', obs_p END IF END IF lfilter ! set observation offset thisobs%off_obs_f = offset_obs offset_obs = offset_obs + thisobs%dim_obs_f ! set global observation offset for EnKF/LEnKF IF (TRIM(filterstr)=='ENKF' .OR. TRIM(filterstr)=='LENKF') THEN thisobs%off_obs_g = offset_obs_g offset_obs_g = offset_obs_g + thisobs%dim_obs_g IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%off_obs_g', thisobs%off_obs_g END IF END IF ! Set general flag if observations are omitted for large innovation IF (thisobs%inno_omit > 0.0) omit_obs = .TRUE. ! Initialize statistics for observations omitted for large innovation ostats_omit = 0 ! Screen output IF (mype == 0 .AND. screen > 0.AND. thisobs%inno_omit > 0.0) THEN WRITE (*, '(a, 5x, a, 1x, i3, 1x , a, f8.2,a)') & 'PDAFomi', '--- Exclude obs. type ID', n_obstypes, ' if innovation^2 > ', & thisobs%inno_omit,' times obs. error variance' END IF ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%obs_f', thisobs%obs_f WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%ivar_obs_f', thisobs%ivar_obs_f WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'thisobs%ocoord_f', thisobs%ocoord_f WRITE (*,*) '++ OMI-debug gather_obs: ', debug, 'initialized obs. ID', thisobs%obsid WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_gather_obs -- END' END IF END SUBROUTINE PDAFomi_gather_obs !------------------------------------------------------------------------------- !> Gather full observational information !! !! This routine uses PDAFomi_gather_obs_f_flex to obtain !! a full observed state vector. The routine is usually called !! the observation operators. !! !! __Revision history:__ !! * 2020-05 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_gather_obsstate(thisobs, obsstate_p, obsstate_f) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), TARGET, INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: obsstate_p(:) !< Vector of process-local observed state REAL, INTENT(inout) :: obsstate_f(:) !< Full observed vector for all types ! *** Local variables *** INTEGER :: status ! Status flag for PDAF gather operation INTEGER :: localfilter ! Whether the filter is domain-localized INTEGER :: globalobs ! Whether the filter needs global observations REAL, ALLOCATABLE :: obsstate_tmp(:) ! Temporary vector of globally full observations ! ************************************** ! *** Gather full observation arrays *** ! ************************************** ! Check whether the filter is domain-localized CALL PDAF_get_localfilter(localfilter) ! Check whether the filter needs global observations CALL PDAF_get_globalobs(globalobs) ! Print debug information IF (debug>0) THEN IF (obs_member==0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_gather_obsstate -- START Gather full observed ensemble mean' ELSE WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_gather_obsstate -- START Gather full observed ensemble state', obs_member END IF WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, 'observation ID', thisobs%obsid WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, 'thisobs%dim_obs_p', thisobs%dim_obs_p WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, 'thisobs%dim_obs_f', thisobs%dim_obs_f WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, 'thisobs%off_obs_f', thisobs%off_obs_f IF (thisobs%use_global_obs==0) & WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, 'thisobs%dim_obs_g', thisobs%dim_obs_g WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, 'obsstate_p', obsstate_p END IF lfilter: IF (localfilter==1 .OR. globalobs==1) THEN ! For domain-localized filters: gather full observations fullobs: IF (thisobs%use_global_obs==1 .OR. globalobs==1) THEN ! *** Gather global full observation vector *** CALL PDAFomi_gather_obs_f_flex(thisobs%dim_obs_p, obsstate_p, & obsstate_f(thisobs%off_obs_f+1 : thisobs%off_obs_f+thisobs%dim_obs_f), status) ELSE fullobs ! *** Use full observations limited to those relevant for a process domain *** ! *** This can be more efficient as in the local analysis loop less *** ! *** observations have a be checked for each analysis domain *** IF (thisobs%dim_obs_g>0) THEN ALLOCATE(obsstate_tmp(thisobs%dim_obs_g)) ELSE ALLOCATE(obsstate_tmp(1)) END IF ! *** Gather observation vector *** CALL PDAFomi_gather_obs_f_flex(thisobs%dim_obs_p, obsstate_p, & obsstate_tmp, status) ! Now restrict observation vector to process-relevant part CALL PDAFomi_limit_obs_f(thisobs, thisobs%off_obs_f, obsstate_tmp, obsstate_f) DEALLOCATE(obsstate_tmp) END IF fullobs ELSE lfilter ! *** For global filters use process-local observations without gathering *** ! In case of a global filter store process-local observed state obsstate_f(thisobs%off_obs_f+1 : thisobs%off_obs_f+thisobs%dim_obs_p) & = obsstate_p(1:thisobs%dim_obs_p) END IF lfilter IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, & 'obsstate_f', obsstate_f(thisobs%off_obs_f+1 : thisobs%off_obs_f+thisobs%dim_obs_p) END IF ! Initialize pointer array IF (obscnt == 0) THEN IF (.NOT.ALLOCATED(obs_f_all)) ALLOCATE(obs_f_all(n_obstypes)) obscnt = 1 IF (debug>0) THEN WRITE (*,*) '++ OMI-debug gather_obsstate: ', debug, & 'initialize pointer array for ', n_obstypes, 'observation types' END IF END IF ! Set pointer to current observation obs_f_all(thisobs%obsid)%ptr => thisobs IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, ' PDAFomi_gather_obsstate -- END' END IF END SUBROUTINE PDAFomi_gather_obsstate !------------------------------------------------------------------------------- !> Initialize full vector of observations !! !! This routine initializes the part of the full vector of !! observations for the current observation type. !! It has to fill the observations to obsstate_f from !! position OFFSET_OBS+1. For the return value OFFSET_OBS !! has to be incremented by the number of added observations. !! The routine will only be called if the adaptive forgetting !! factor is used. !! !! __Revision history:__ !! * 2019-09 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obs_f(thisobs, dim_obs_f, obsstate_f, offset) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: dim_obs_f !< Dimension of full observed state (all observed fields) REAL, INTENT(inout) :: obsstate_f(:) !< Full observation vector (dim_obs_f) INTEGER, INTENT(inout) :: offset !< input: offset of module-type observations in obsstate_f !< output: input + number of added observations ! ****************************************** ! *** Initialize full observation vector *** ! ****************************************** ! Consistency check IF (dim_obs_f < offset+thisobs%dim_obs_f) THEN WRITE (*,'(a)') 'PDAFomi - ERROR: PDAFomi_init_obs_f - dim_obs_f is too small !!!' error = 2 END IF doassim: IF (thisobs%doassim == 1) THEN ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_init_obs_f -- START Initialize observation vector' WRITE (*,*) '++ OMI-debug init_obs_f: ', debug, 'observation ID', thisobs%obsid WRITE (*,*) '++ OMI-debug init_obs_f: ', debug, 'thisobs%dim_obs_f', thisobs%dim_obs_f WRITE (*,*) '++ OMI-debug init_obs_f: ', debug, 'thisobs%obs_f', thisobs%obs_f END IF ! Fill part of full observation vector obsstate_f(offset+1 : offset+thisobs%dim_obs_f) = thisobs%obs_f(1 : thisobs%dim_obs_f) ! Increment offset offset = offset + thisobs%dim_obs_f IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_init_obs_f -- END' END IF END IF doassim END SUBROUTINE PDAFomi_init_obs_f !------------------------------------------------------------------------------- !> Compute mean observation error variance !! !! This routine will only be called, if the adaptive !! forgetting factor feature is used. Please note that !! this is an experimental feature. !! !! The routine is called in global filters (like ESTKF) !! during the analysis or in local filters (e.g. LESTKF) !! before the loop over local analysis domains !! by the routine PDAF_set_forget that estimates an !! adaptive forgetting factor. The routine has to !! initialize the mean observation error variance. !! For global filters this should be the global mean, !! while for local filters it should be the mean for the !! PE-local sub-domain. (init_obsvar_l_TYPE is the !! localized variant for local filters) !! !! The routine assumes a diagonal observation error !! covariance matrix. !! !! If the observation counter is zero the computation !! of the mean variance is initialized. The output is !! always the mean variance. If the observation counter !! is >0 the computed was initialized for another !! observation. To proceed the computation, meanvar is !! first multiplied with the observation counter to !! obtain the variance sum. Then the computation of the !! mean is continued. !! !! __Revision history:__ !! * 2019-09 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obsvar_f(thisobs, meanvar, cnt_obs) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(inout) :: meanvar !< Mean variance INTEGER, INTENT(inout) :: cnt_obs !< Observation counter ! Local variables *** INTEGER :: i ! Counter ! *********************************** ! *** Compute local mean variance *** ! *********************************** doassim: IF (thisobs%doassim == 1) THEN IF (cnt_obs==0) THEN ! Reset mean variance meanvar = 0.0 ELSE ! Compute sum of variances from mean variance meanvar = meanvar * REAL(cnt_obs) END IF ! Add observation error variances DO i = 1, thisobs%dim_obs_f meanvar = meanvar + 1.0 / thisobs%ivar_obs_f(i) END DO ! Increment observation count cnt_obs = cnt_obs + thisobs%dim_obs_f ! Compute updated mean variance meanvar = meanvar / REAL(cnt_obs) END IF doassim END SUBROUTINE PDAFomi_init_obsvar_f !------------------------------------------------------------------------------- !> Compute product of inverse of R with some matrix !! !! The routine is called during the analysis step !! of the global square-root filters. It has to !! compute the product of the inverse of the !! process-local observation error covariance matrix !! with the matrix of process-local observed ensemble !! perturbations. !! !! This routine assumes a diagonal observation error !! covariance matrix, but allows for varying observation !! error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. The operation is done with all !! process-local observations !! !! __Revision history:__ !! * 2019-12 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_prodRinvA(thisobs, ncols, A_p, C_p) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: ncols !< Number of columns in A_p and C_p REAL, INTENT(in) :: A_p(:, :) !< Input matrix (nobs_f, ncols) REAL, INTENT(out) :: C_p(:, :) !< Output matrix (nobs_f, ncols) ! *** local variables *** INTEGER :: i, j ! index of observation component INTEGER :: off ! row offset in A_l and C_l ! ************************************* ! *** -1 *** ! *** C = R A *** ! *** *** ! *** The inverse observation error *** ! *** covariance matrix is not *** ! *** computed explicitely. *** ! ************************************* doassim: IF (thisobs%doassim == 1) THEN ! Check process-local observation dimension IF (thisobs%dim_obs_p /= thisobs%dim_obs_f) THEN ! This error usually happens when localfilter=1 WRITE (*,'(a)') 'PDAFomi - ERROR: PDAFomi_prodRinvA - INCONSISTENT value for DIM_OBS_P !!!' error = 3 END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_prodRinvA -- START Multiply with inverse observation variance' WRITE (*,*) '++ OMI-debug prodRinvA: ', debug, 'observation ID', thisobs%obsid WRITE (*,*) '++ OMI-debug prodRinvA: ', debug, 'thisobs%dim_obs_f', thisobs%dim_obs_f WRITE (*,*) '++ OMI-debug prodRinvA: ', debug, 'thisobs%ivar_obs_f', thisobs%ivar_obs_f WRITE (*,*) '++ OMI-debug prodRinvA: ', debug, 'Input matrix A_p', A_p END IF ! Initialize offset off = thisobs%off_obs_f DO j = 1, ncols DO i = 1, thisobs%dim_obs_f C_p(i+off, j) = thisobs%ivar_obs_f(i) * A_p(i+off, j) END DO END DO IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_prodRinvA -- END' END IF END IF doassim END SUBROUTINE PDAFomi_prodRinvA !------------------------------------------------------------------------------- !> Compute likelihood for an ensemble member !! !! The routine is called during the analysis step !! of the NETF or a particle filter. !! It has to compute the likelihood of the !! ensemble according to the difference from the !! observation (residual) and the error distribution !! of the observations. !! !! In general this routine is similar to the routine !! prodRinvA used for ensemble square root Kalman !! filters. As an addition to this routine, we here have !! to evaluate the likelihood weight according the !! assumed observation error statistics. !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_likelihood(thisobs, resid, lhood) USE PDAF_mod_filter, & ONLY: obs_member IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: resid(:) ! Input vector of residuum REAL, INTENT(inout) :: lhood ! Output vector - log likelihood ! *** local variables *** INTEGER :: i ! index of observation component REAL, ALLOCATABLE :: Rinvresid(:) ! R^-1 times residual REAL :: lhood_one ! Likelihood for this observation doassim: IF (thisobs%doassim == 1) THEN ! **************************************** ! *** First scale by observation error *** ! *** -1 *** ! *** Rinvresid = R resid *** ! *** *** ! *** We assume a diagonal matrix R *** ! **************************************** ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_likelihood -- START Compute likelihood, member', obs_member IF (thisobs%obs_err_type==0) THEN WRITE (*,*) '++ OMI-debug likelihood: ', debug, & ' Assume Gaussian observation errors' ELSE WRITE (*,*) '++ OMI-debug likelihood: ', debug, & ' Assume double-exponential observation errors' END IF END IF ! Compute product of R^-1 with residuum ALLOCATE(Rinvresid(thisobs%dim_obs_f)) DO i = 1, thisobs%dim_obs_f Rinvresid(i) = thisobs%ivar_obs_f(i) * resid(thisobs%off_obs_f+i) END DO ! ****************************** ! *** Compute log likelihood *** ! ****************************** IF (thisobs%obs_err_type==0) THEN ! Gaussian errors ! Calculate exp(-0.5*resid^T*R^-1*resid) ! Transform back to log likelihood to increment its values IF (lhood>0.0) lhood = - LOG(lhood) lhood_one = 0.0 DO i = 1, thisobs%dim_obs_f lhood_one = lhood_one + 0.5*resid(thisobs%off_obs_f+i)*Rinvresid(i) END DO lhood = EXP(-(lhood + lhood_one)) ELSE ! Double-exponential errors ! Calculate exp(-SUM(ABS(resid))) ! Transform pack to log likelihood to increment its values IF (lhood>0.0) lhood = - LOG(lhood) lhood_one = 0.0 DO i = 1, thisobs%dim_obs_f lhood_one = lhood_one + ABS(Rinvresid(i)) END DO lhood = EXP(-(lhood + lhood_one)) END IF ! *** Clean up *** DEALLOCATE(Rinvresid) ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug likelihood: ', debug, ' accumulated likelihood', lhood WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_likelihood -- END' END IF END IF doassim END SUBROUTINE PDAFomi_likelihood !------------------------------------------------------------------------------- !> Add observation error to some matrix !! !! The routine is called during the analysis step !! of the stochastic EnKF. It it provided with a !! matrix in observation space and has to add the !! observation error covariance matrix. !! !! This routine assumes a diagonal observation error !! covariance matrix, but allows for varying observation !! error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. The operation is done with all !! process-local observations !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_add_obs_error(thisobs, nobs_all, matC) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nobs_all !< Number of observations REAL, INTENT(inout) :: matC(:, :) !< Input/Output matrix (nobs_f, rank) ! *** local variables *** INTEGER :: i, pe, cnt ! Counters INTEGER :: idummy ! Dummy to access nobs_all INTEGER, ALLOCATABLE :: id_start(:) ! Start index of obs. type in global averall obs. vector INTEGER, ALLOCATABLE :: id_end(:) ! End index of obs. type in global averall obs. vector doassim: IF (thisobs%doassim == 1) THEN ! Initialize dummy to prevent compiler warning idummy = nobs_all ! ************************************************* ! *** Check process-local observation dimension *** ! ************************************************* IF (thisobs%dim_obs_p /= thisobs%dim_obs_f) THEN ! This error usually happens when localfilter=1 WRITE (*,'(a)') 'PDAFomi - ERROR: PDAFomi_add_obs_error - INCONSISTENT VALUE for DIM_OBS_P !!!' error = 4 END IF ! ******************************************************** ! *** Initialize indices of observation type in global *** ! *** state vector over all observation types. *** ! ******************************************************** ALLOCATE(id_start(npes), id_end(npes)) pe = 1 id_start(1) = 1 IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 DO pe = 2, npes id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 END DO ! ************************************* ! *** Add observation error *** ! *** *** ! *** Measurements are uncorrelated *** ! *** here, thus R is diagonal *** ! ************************************* cnt = 1 DO pe = 1, npes DO i = id_start(pe), id_end(pe) matC(i, i) = matC(i, i) + 1.0/thisobs%ivar_obs_f(cnt) cnt = cnt + 1 ENDDO ENDDO DEALLOCATE(id_start, id_end) END IF doassim END SUBROUTINE PDAFomi_add_obs_error !------------------------------------------------------------------------------- !> Initialize global observation error covariance matrix !! !! The routine is called during the analysis !! step when an ensemble of observations is !! generated by PDAF_enkf_obs_ensemble. !! It has to initialize the global observation !! error covariance matrix. !! !! This routine assumes a diagonal observation error !! covariance matrix, but allows for varying observation !! error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. The operation is done with all !! process-local observations !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obscovar(thisobs, nobs_all, covar, isdiag) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nobs_all !< Number of observations REAL, INTENT(inout) :: covar(:, :) !< Input/Output matrix (nobs_all, nobs_all) !< (needs to be set =0 before calling the routine) LOGICAL, INTENT(out) :: isdiag !< Whether matrix R is diagonal ! *** local variables *** INTEGER :: i, pe, cnt ! Counters INTEGER :: idummy ! Dummy to access nobs_all INTEGER, ALLOCATABLE :: id_start(:) ! Start index of obs. type in global averall obs. vector INTEGER, ALLOCATABLE :: id_end(:) ! End index of obs. type in global averall obs. vector doassim: IF (thisobs%doassim == 1) THEN ! Initialize dummy to prevent compiler warning idummy = nobs_all ! ************************************************* ! *** Initialize indices of observation type in *** ! *** global state vector over all obsservation *** ! *** types and corresponding mapping vector. *** ! ************************************************* ALLOCATE(id_start(npes), id_end(npes)) ! Initialize indices pe = 1 id_start(1) = 1 IF (thisobs%obsid>1) id_start(1) = id_start(1) + sum(obsdims(1, 1:thisobs%obsid-1)) id_end(1) = id_start(1) + obsdims(1,thisobs%obsid) - 1 DO pe = 2, npes id_start(pe) = id_start(pe-1) + SUM(obsdims(pe-1,thisobs%obsid:)) IF (thisobs%obsid>1) id_start(pe) = id_start(pe) + sum(obsdims(pe,1:thisobs%obsid-1)) id_end(pe) = id_start(pe) + obsdims(pe,thisobs%obsid) - 1 END DO ! Initialize mapping vector (to be used in PDAF_enkf_obs_ensemble) cnt = 1 IF (thisobs%obsid-1 > 0) cnt = cnt+ SUM(obsdims(:,1:thisobs%obsid-1)) DO pe = 1, npes DO i = id_start(pe), id_end(pe) map_obs_id(i) = cnt cnt = cnt + 1 END DO END DO ! ************************************* ! *** Initialize covariances *** ! *** *** ! *** Measurements are uncorrelated *** ! *** here, thus R is diagonal *** ! ************************************* cnt = 1 DO pe = 1, npes DO i = id_start(pe), id_end(pe) covar(i, i) = covar(i, i) + 1.0/thisobs%ivar_obs_f(cnt) cnt = cnt + 1 ENDDO ENDDO ! The matrix is diagonal ! This setting avoids the computation of the SVD of COVAR ! in PDAF_enkf_obs_ensemble isdiag = .TRUE. DEALLOCATE(id_start, id_end) END IF doassim END SUBROUTINE PDAFomi_init_obscovar !------------------------------------------------------------------------------- !> Initialize vector of observation errors for generating synthetic obs. !! !! This routine initializes a vector of observation errors !! used to perturb an observe dmodel state when generating !! synthetic observations. !! !! __Revision history:__ !! * 2020-05 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obserr_f(thisobs, obserr_f) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation REAL, INTENT(inout) :: obserr_f(:) !< Full vector of observation errors ! *** Local variables *** INTEGER :: i ! Counter ! ***************************************************************************** ! *** Initialize vector of observation errors for generating synthetic obs. *** ! ***************************************************************************** doassim: IF (thisobs%doassim == 1) THEN DO i = 1, thisobs%dim_obs_f IF (thisobs%ivar_obs_f(i)>0.0) THEN obserr_f(i + thisobs%off_obs_f) = 1.0/SQRT(thisobs%ivar_obs_f(i)) ELSE obserr_f(i + thisobs%off_obs_f) = 1.0e12 END IF END DO END IF doassim END SUBROUTINE PDAFomi_init_obserr_f !------------------------------------------------------------------------------- !> Set min/max coordinate locations of a decomposed grid !! !! This routine sets the limiting coordinates of a !! process domain, i.e. the northern-, southern-, !! eastern-, and western-most coordinate. The !! information can be used to restrict the full !! observations for PDAF to those that might be !! used for the local analysis. !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_set_domain_limits(lim_coords) IMPLICIT NONE ! *** Arguments *** REAL, INTENT(in) :: lim_coords(2,2) !< geographic coordinate array (1: longitude, 2: latitude) !< ranges: longitude (-pi, pi), latitude (-pi/2, pi/2) ! Store domain limiting coordinates in module array IF (.NOT.ALLOCATED(domain_limits)) ALLOCATE(domain_limits(4)) domain_limits(1) = lim_coords(2,1) ! Northern edge domain_limits(2) = lim_coords(2,2) ! Southern edge domain_limits(3) = lim_coords(1,1) ! Western edge domain_limits(4) = lim_coords(1,2) ! Eastern edge IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_set_domain_limits -- START' WRITE (*,*) '++ OMI-debug set_domain_limits: ', debug, 'domain limits', domain_limits WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_set_domain_limits -- END' END IF END SUBROUTINE PDAFomi_set_domain_limits !------------------------------------------------------------------------------- !> Find min/max coordinate locations in unstructured grid !! !! This routine finds the limiting coordinates of a !! process domain, i.e. the northern-, southern-, !! eastern-, and western-most coordinate. The !! information can be used to restrict the full !! observations for PDAF to those that might be !! used for the local analysis. Usually this is used !! for unstructured grids, like in FESOM, because the !! grid point indices do not contain information on !! coordinates in this case. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_get_domain_limits_unstr(npoints_p, coords_p) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: npoints_p !< number of process-local grid points REAL, INTENT(in) :: coords_p(:,:) !< geographic coordinate array (row 1: longitude, 2: latitude) !< ranges: longitude (-pi, pi), latitude (-pi/2, pi/2) ! *** Local variables *** INTEGER :: i ! Counter REAL :: nlimit, slimit, elimit, wlimit ! Limiting coordinates REAL :: abslonmin ! absolute minimum longitude ! *** Determine limiting coordinates *** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_get_domain_limits_unstr -- START' END IF ! Initialize limiting values nlimit = -100.0 slimit = 100.0 wlimit = 100.0 elimit = -100.0 abslonmin = 100.0 DO i=1, npoints_p ! Get North/South Limits IF (coords_p(2,i) < slimit) slimit = coords_p(2,i) IF (coords_p(2,i) > nlimit) nlimit = coords_p(2,i) ! Get East/West Limits IF (coords_p(1,i) < wlimit) wlimit = coords_p(1,i) IF (coords_p(1,i) > elimit) elimit = coords_p(1,i) IF (ABS(coords_p(1,i)) < abslonmin) THEN abslonmin = ABS(coords_p(1,i)) END IF ENDDO IF (elimit*wlimit<0.0) THEN ! Domain crosses prime meridian or date line IF (wlimit<-3.1 .AND. elimit>3.1 .AND. abslonmin>0.5) THEN ! If the domain crosses the date line, we have to search the longitudinal limits differently elimit = -100.0 wlimit = 100.0 DO i=1, npoints_p IF (coords_p(1,i)<0.0 .AND. coords_p(1,i)>elimit) elimit = coords_p(1,i) IF (coords_p(1,i)>0.0 .AND. coords_p(1,i)<wlimit) wlimit = coords_p(1,i) END DO IF (debug>0) THEN WRITE (*,*) '++ OMI-debug get_domain_limits_unstr: ', debug, & 'limits (crossing date line) ', nlimit, slimit, wlimit, elimit END IF ELSE ! In this case the domain crosses the prime meridian IF (debug>0) THEN WRITE (*,*) '++ OMI-debug get_domain_limits_unstr: ', debug, & 'limits (cossing prime meridian) ', nlimit, slimit, wlimit, elimit END IF END IF ELSE ! Standard case IF (debug>0) THEN WRITE (*,*) '++ OMI-debug get_domain_limits_unstr: ', debug, & 'limits (standard case) ', nlimit, slimit, wlimit, elimit END IF END IF ! Store domain limiting coordinates in module array IF (.NOT.ALLOCATED(domain_limits)) ALLOCATE(domain_limits(4)) domain_limits(1) = nlimit domain_limits(2) = slimit domain_limits(3) = wlimit domain_limits(4) = elimit IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_get_domain_limits_unstr -- END' END IF END SUBROUTINE PDAFomi_get_domain_limits_unstr !------------------------------------------------------------------------------- !> Find observations inside or close to process domain !! !! This routine finds observations that lie inside the !! local process sub-domain or within the distance !! LRADIUS around it. The observations are counted and !! an index array is initialized storing the indices !! of the process-local relevant full observations in the !! global full observation vector. !! !! The routine has to be called by all filter processes !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_get_local_ids_obs_f(dim_obs_g, lradius, oc_f, cnt_lim, id_lim, disttype, domainsize) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: dim_obs_g !< Global full number of observations REAL, INTENT(in) :: lradius !< Localization radius (used is a constant one here) REAL, INTENT(in) :: oc_f(:,:) !< observation coordinates (radians), row 1: lon, 2: lat !< ranges: longitude (-pi, pi), latitude (-pi/2, pi/2) INTEGER, INTENT(out) :: cnt_lim !< Number of full observation for local process domain INTEGER, INTENT(out) :: id_lim(:) !< Indices of process-local full obs. in global full vector !< It has to be allocated sufficiently large INTEGER, INTENT(in) :: disttype !< type of distance computation REAL, INTENT(in), OPTIONAL :: domainsize(:) !< Global size of model domain ! *** Local variables *** INTEGER :: i ! Counter INTEGER :: flag ! Counting flag REAL :: limdist ! Limit distance normalized by r_earth INTEGER :: cnt_lim_max, cnt_lim_min ! min/max number over all domains REAL :: maxlat ! Highest latitude of a domain REAL :: period(2) ! Size of domain in case of periodicity ! ********************** ! *** Initialization *** ! ********************** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_get_local_ids_obs_f -- START Limit observations to process sub-domains' WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'domain limiting coordinates', domain_limits END IF IF (.NOT.ALLOCATED(domain_limits)) THEN WRITE (*,'(a)') 'PDAFomi - ERROR: PDAFomi_get_local_ids_obs_f - DOMAIN_LIMITS is not initialized !!!' error = 5 END IF IF ((disttype==1 .OR. disttype==11) .AND. .NOT.PRESENT(domainsize)) THEN WRITE (*,'(a)') 'PDAFomi - ERROR: PDAFomi_get_local_ids_obs_f - THISOBS%DOMAINSIZE is not initialized !!!' error = 6 END IF ! initialize index array id_lim = 0 ! *************************************** ! *** Find relevant full observations *** ! *************************************** cnt_lim = 0 dtype: IF (disttype==2 .OR. disttype==3 .OR. disttype==12 .OR. disttype==13) THEN ! Limit distance around the domain limdist = lradius / r_earth IF (debug==1) THEN WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'limit for geographic coordinates' WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'limiting distance (m)', limdist END IF fullobsloop: DO i = 1, dim_obs_g ! Init flag for latitudinal check flag = 0 ! First check in latitudinal direction checklat: IF (oc_f(2,i)<=domain_limits(1) .AND. oc_f(2,i)>=domain_limits(2)) THEN ! inside domain north-south extent flag=1 ELSEIF (oc_f(2,i)>domain_limits(1)) THEN ! north of the domain IF (ABS(oc_f(2,i)-domain_limits(1)) <= limdist) flag=1 ELSEIF (oc_f(2,i)<domain_limits(2)) THEN ! south of the domain IF (ABS(oc_f(2,i)-domain_limits(2)) <= limdist) flag=1 END IF checklat ! Store highest latitude maxlat = MAX(ABS(domain_limits(1)), ABS(domain_limits(2))) ! if observation fits in the latitudinal direction check longitudinal direction lat_ok: IF (flag==1) THEN lontypes: IF (domain_limits(4)>=0.0 .OR. (domain_limits(4)<0.0 .AND. domain_limits(3)<0.0)) THEN IF (oc_f(1,i)>=domain_limits(3) .AND. oc_f(1,i)<=domain_limits(4)) THEN ! fully inside domain extent cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i ELSEIF (oc_f(1,i)<domain_limits(3)) THEN ! west of the domain IF (ABS(COS(maxlat)*(oc_f(1,i)-domain_limits(3))) <= limdist .OR. & (ABS(COS(maxlat)*(oc_f(1,i)-domain_limits(3)))-2.0*pi) <= limdist ) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i END IF ELSEIF (oc_f(1,i)>domain_limits(4)) THEN ! east of the domain IF (ABS(COS(maxlat)*(oc_f(1,i)-domain_limits(4))) <= limdist .OR. & (ABS(COS(maxlat)*(oc_f(1,i)-domain_limits(4)))-2.0*pi) <= limdist ) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i END IF ENDIF ELSE lontypes IF ((oc_f(1,i)>=domain_limits(3) .AND. oc_f(1,i)<=pi) .OR. & (oc_f(1,i)<=domain_limits(4)) .AND. oc_f(1,i)>=-pi) THEN ! fully inside domain extent cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i ELSEIF (oc_f(1,i)<domain_limits(3) .AND. oc_f(1,i)>=0.0) THEN ! east of the domain IF (ABS(COS(maxlat)*(oc_f(1,i)-domain_limits(3))) <= limdist) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i ENDIF ELSEIF (oc_f(1,i)>domain_limits(4)) THEN ! west of the domain IF (ABS(COS(maxlat)*(oc_f(1,i)-domain_limits(4))) <= limdist) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i ENDIF ENDIF ENDIF lontypes ENDIF lat_ok END DO fullobsloop ELSE IF (disttype==0 .OR. disttype==10) THEN ! *** Check Cartesian coordinates without periodicity *** limdist = lradius IF (debug==1) THEN WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'limit for Cartesian coordinates' WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'limiting distance', limdist END IF fullobsloopB: DO i = 1, dim_obs_g ! Init flag for latitudinal check flag = 0 ! First check in latitudinal direction checklatC: IF (oc_f(2,i)<=domain_limits(1) .AND. oc_f(2,i)>=domain_limits(2)) THEN ! inside domain north-south extent flag=1 ELSEIF (oc_f(2,i)>domain_limits(1)) THEN ! north of the domain IF (ABS(oc_f(2,i)-domain_limits(1)) <= limdist) flag=1 ELSEIF (oc_f(2,i)<domain_limits(2)) THEN ! south of the domain IF (ABS(oc_f(2,i)-domain_limits(2)) <= limdist) flag=1 END IF checklatC ! if observation fits in the latitudinal direction check longitudinal direction lat_okB: IF (flag==1) THEN IF (oc_f(1,i)>=domain_limits(3) .AND. oc_f(1,i)<=domain_limits(4)) THEN ! fully inside domain extent cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i ELSEIF (oc_f(1,i)<domain_limits(3)) THEN ! West of the domain IF (ABS(oc_f(1,i)-domain_limits(3)) <= limdist) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i END IF ELSEIF (oc_f(1,i)>domain_limits(4)) THEN ! East of the domain IF (ABS(oc_f(1,i)-domain_limits(4)) <= limdist) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i END IF ENDIF END IF lat_okB END DO fullobsloopB ELSE IF (disttype==1 .OR. disttype==11) THEN ! *** Check Cartesian coordinates with periodicity *** limdist = lradius IF (debug==1) THEN WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'limit for periodic Cartesian coordinates' WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'limiting distance', limdist WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'thisobs%domainsize', domainsize END IF fullobsloopC: DO i = 1, dim_obs_g ! Init flag for latitudinal check flag = 0 IF (domainsize(1)>0.0) THEN period(1) = domainsize(1) ELSE period(1) = 0.0 END IF IF (domainsize(2)>0.0) THEN period(2) = domainsize(2) ELSE period(2) = 0.0 END IF ! First check in latitudinal direction checklatB: IF (oc_f(2,i)<=domain_limits(1) .AND. oc_f(2,i)>=domain_limits(2)) THEN ! inside domain north-south extent flag=1 ELSEIF (oc_f(2,i)>domain_limits(1)) THEN ! north of the domain IF ((ABS(oc_f(2,i)-domain_limits(1)) <= limdist) .OR. & (ABS(oc_f(2,i)-domain_limits(1) - period(2)) <= limdist)) flag=1 ELSEIF (oc_f(2,i)<domain_limits(2)) THEN ! south of the domain IF ((ABS(oc_f(2,i)-domain_limits(2)) <= limdist) .OR. & (ABS(oc_f(2,i)-domain_limits(2) + period(2)) <= limdist)) flag=1 END IF checklatB ! if observation fits in the latitudinal direction check longitudinal direction lat_okC: IF (flag==1) THEN IF (oc_f(1,i)>=domain_limits(3) .AND. oc_f(1,i)<=domain_limits(4)) THEN ! fully inside domain extent cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i ELSEIF (oc_f(1,i)<domain_limits(3)) THEN ! West of the domain IF (ABS(oc_f(1,i)-domain_limits(3)) <= limdist .OR. & (ABS(oc_f(1,i)-domain_limits(3)) - period(1)) <= limdist) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i END IF ELSEIF (oc_f(1,i)>domain_limits(4)) THEN ! East of the domain IF (ABS(oc_f(1,i)-domain_limits(4)) <= limdist .OR. & (ABS(oc_f(1,i)-domain_limits(4)) - period(1)) <= limdist) THEN cnt_lim = cnt_lim+1 id_lim(cnt_lim) = i END IF ENDIF END IF lat_okC END DO fullobsloopC END IF dtype IF (debug>0) THEN WRITE (*,*) '++ OMI-debug get_local_ids_obs_f: ', debug, 'obs. ids for process domains', id_lim(1:cnt_lim) WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_get_local_ids_obs_f -- END' END IF ! Get number of min/max process-local full observation dimensions CALL MPI_Allreduce (cnt_lim, cnt_lim_max, 1, MPI_INTEGER, MPI_MAX, & COMM_filter, MPIerr) CALL MPI_Allreduce (cnt_lim, cnt_lim_min, 1, MPI_INTEGER, MPI_MIN, & COMM_filter, MPIerr) IF (mype == 0 .AND. screen > 0) THEN WRITE (*,'(a,8x,a,i8)') 'PDAFomi','--- global observation dimension', dim_obs_g WRITE (*,'(a,8x,a,i7,1x,i7)') 'PDAFomi','--- process-local min/max full obs. dimensions', & cnt_lim_min, cnt_lim_max END IF END SUBROUTINE PDAFomi_get_local_ids_obs_f !------------------------------------------------------------------------------- !> Reduce full observation vector to part relevant for local process domain !! !! This routine initializes a full vector of observations that only !! contains those full observations that are relevant for a process !! subdomain. The indices of these observations were determined !! using get_local_ids_obs_f. !! !! __Revision history:__ !! * 2019-07 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_limit_obs_f(thisobs, offset, obs_f_one, obs_f_lim) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: obs_f_one(:) !< Global full observation vector (nobs_f) REAL, INTENT(out) :: obs_f_lim(:) !< full observation vector for process domains (nobs_lim) INTEGER, INTENT(in) :: offset !< offset of this observation in obs_f_lim ! *** Local variables *** INTEGER :: i ! Counter ! ******************************************** ! *** Initialize process-local full vector *** ! ******************************************** IF (.NOT.ALLOCATED(thisobs%id_obs_f_lim)) THEN WRITE (*,'(a)') 'PDAFomi - ERROR: PDAFomi_limit_obs_f - thisobs%id_obs_f_lim is not allocated !!!' error = 7 END IF DO i = 1, thisobs%dim_obs_f obs_f_lim(i+offset) = obs_f_one(thisobs%id_obs_f_lim(i)) END DO END SUBROUTINE PDAFomi_limit_obs_f SUBROUTINE PDAFomi_gather_obs_f_flex(dim_obs_p, obs_p, obs_f, status) ! !DESCRIPTION: ! If the local filter is used with a domain-decomposed model, ! the observational information from different sub-domains ! has to be combined into the full observation vector. ! In this routine the process-local parts of the observation ! vector are gathered into a full observation vector. ! The routine requires that PDAF_gather_dim_obs_f was executed ! before, because this routine initializes dimensions that are ! used here. ! The routine can also be used to gather full arrays of coordinates. ! It is however, only usable if the coordinates are stored row- ! wise, i.e. each row represents the set of coordinates for one ! observation point. It has to be called separately for each column. ! A better alternative is the row-wise storage of coordinates. In this ! case the routine PDAF_gather_dim_obs_f allows the gather the full ! coordinate array in one step. ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2019-03 - Lars Nerger - Initial code ! Later revisions - see svn log ! ! !USES: ! Include definitions for real type of different precision ! (Defines BLAS/LAPACK routines and MPI_REALTYPE) #include "typedefs.h" USE PDAF_mod_filtermpi, & ONLY: COMM_filter, MPIerr, mype_filter, npes_filter IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: dim_obs_p ! PE-local observation dimension REAL, INTENT(in) :: obs_p(:) ! PE-local vector REAL, INTENT(out) :: obs_f(:) ! Full gathered vector INTEGER, INTENT(out) :: status ! Status flag: (0) no error ! !CALLING SEQUENCE: ! Called by: user code ! Calls: MPI_Allreduce ! Calls: MPI_Allgather ! Calls: MPI_AllgatherV !EOP ! Local variables INTEGER :: i ! Counter INTEGER :: dimobs_f ! full dimension of observation vector obtained from allreduce INTEGER, ALLOCATABLE :: all_dim_obs_p(:) ! PE-Local observation dimensions INTEGER, ALLOCATABLE :: all_dis_obs_p(:) ! PE-Local observation displacements ! ********************************************************** ! *** Compute global sum of local observation dimensions *** ! ********************************************************** IF (npes_filter>1) THEN CALL MPI_Allreduce(dim_obs_p, dimobs_f, 1, MPI_INTEGER, MPI_SUM, & COMM_filter, MPIerr) ELSE dimobs_f = dim_obs_p END IF ! **************************************************************************** ! *** Gather and store array of process-local dimensions and displacements *** ! **************************************************************************** ALLOCATE(all_dim_obs_p(npes_filter)) ALLOCATE(all_dis_obs_p(npes_filter)) IF (npes_filter>1) THEN CALL MPI_Allgather(dim_obs_p, 1, MPI_INTEGER, all_dim_obs_p, 1, & MPI_INTEGER, COMM_filter, MPIerr) ! Init array of displacements for observation vector all_dis_obs_p(1) = 0 DO i = 2, npes_filter all_dis_obs_p(i) = all_dis_obs_p(i-1) + all_dim_obs_p(i-1) END DO ELSE all_dim_obs_p = dim_obs_p all_dis_obs_p = 0 END IF ! ********************************************************** ! *** Gather full observation vector *** ! ********************************************************** IF (npes_filter>1) THEN CALL MPI_AllGatherV(obs_p, all_dim_obs_p(mype_filter+1), MPI_REALTYPE, & obs_f, all_dim_obs_p, all_dis_obs_p, MPI_REALTYPE, & COMM_filter, MPIerr) status = MPIerr ELSE obs_f = obs_p status = 0 END IF ! **************** ! *** Clean up *** ! **************** DEALLOCATE(all_dim_obs_p, all_dis_obs_p) END SUBROUTINE PDAFomi_gather_obs_f_flex SUBROUTINE PDAFomi_gather_obs_f2_flex(dim_obs_p, coords_p, coords_f, & nrows, status) ! !DESCRIPTION: ! If the local filter is used with a domain-decomposed model, ! the observational information from different sub-domains ! has to be combined into the full observation vector. ! In this routine the process-local parts of a coordinate array ! accompanying the observation vector are gathered into a full ! array of coordinates. ! The routine is for the case that the observation coordinates ! are stored column-wise, i.e. each column is the set of coordinates ! for one observation. This should be the usual case, as in this ! case the set of coordinates of one observations are stored ! next to each other in memory. If the coordinates are stored row- ! wise, the routine PDAF_gather_obs_f can be used, but has to be ! called separately for each column. ! ! ! This is a core routine of PDAF and ! should not be changed by the user ! ! ! !REVISION HISTORY: ! 2019-03 - Lars Nerger - Initial code ! Later revisions - see svn log ! ! !USES: ! Include definitions for real type of different precision ! (Defines BLAS/LAPACK routines and MPI_REALTYPE) #include "typedefs.h" USE PDAF_mod_filtermpi, & ONLY: COMM_filter, MPIerr, mype_filter, npes_filter IMPLICIT NONE ! !ARGUMENTS: INTEGER, INTENT(in) :: dim_obs_p ! PE-local observation dimension INTEGER, INTENT(in) :: nrows ! Number of rows in array REAL, INTENT(in) :: coords_p(:,:) ! PE-local array REAL, INTENT(out) :: coords_f(:,:) ! Full gathered array INTEGER, INTENT(out) :: status ! Status flag: (0) no error ! !CALLING SEQUENCE: ! Called by: user code ! Calls: MPI_Allreduce ! Calls: MPI_Allgather ! Calls: MPI_AllgatherV !EOP ! local variables INTEGER :: i ! Counter INTEGER :: dimobs_f ! full dimension of observation vector obtained from allreduce INTEGER, ALLOCATABLE :: all_dim_obs_p(:) ! PE-Local observation dimensions INTEGER, ALLOCATABLE :: all_dim_obs_p2(:) ! local-dims for multi-row array INTEGER, ALLOCATABLE :: all_dis_obs_p2(:) ! displacements to gather multi-row array ! ********************************************************** ! *** Compute global sum of local observation dimensions *** ! ********************************************************** IF (npes_filter>1) THEN CALL MPI_Allreduce(dim_obs_p, dimobs_f, 1, MPI_INTEGER, MPI_SUM, & COMM_filter, MPIerr) ELSE dimobs_f = dim_obs_p END IF ! **************************************************************************** ! *** Gather and store array of process-local dimensions and displacements *** ! **************************************************************************** ALLOCATE(all_dim_obs_p(npes_filter)) IF (npes_filter>1) THEN CALL MPI_Allgather(dim_obs_p, 1, MPI_INTEGER, all_dim_obs_p, 1, & MPI_INTEGER, COMM_filter, MPIerr) ELSE all_dim_obs_p = dim_obs_p END IF ! ********************************************************** ! *** Gather full observation coordinates array *** ! ********************************************************** IF (npes_filter>1) THEN ALLOCATE(all_dis_obs_p2(npes_filter)) ALLOCATE(all_dim_obs_p2(npes_filter)) ! Init array of local dimensions do i = 1, npes_filter all_dim_obs_p2(i) = nrows * all_dim_obs_p(i) end do ! Init array of displacements for observation vector all_dis_obs_p2(1) = 0 DO i = 2, npes_filter all_dis_obs_p2(i) = all_dis_obs_p2(i-1) + all_dim_obs_p2(i-1) END DO CALL MPI_AllGatherV(coords_p, all_dim_obs_p2(mype_filter+1), MPI_REALTYPE, & coords_f, all_dim_obs_p2, all_dis_obs_p2, MPI_REALTYPE, & COMM_filter, MPIerr) DEALLOCATE(all_dim_obs_p2, all_dis_obs_p2) status = MPIerr ELSE coords_f = coords_p status = 0 END IF END SUBROUTINE PDAFomi_gather_obs_f2_flex !------------------------------------------------------------------------------- !> Exclude observations for too high innovation !! !! The routine is called during the analysis step !! of a global filter. It checks the size of the !! innovation and sets the observation error to a !! high value if the squared innovation exceeds a !! limit relative to the observation error variance. !! !! __Revision history:__ !! * 2024-02 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_omit_by_inno(thisobs, inno_f, obs_f_all, obsid, cnt_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: inno_f(:) !< Input vector of observation innovation REAL, INTENT(in) :: obs_f_all(:) !< Input vector of local observations INTEGER, INTENT(in) :: obsid !< ID of observation type INTEGER, INTENT(inout) :: cnt_all !< Count of omitted observation over all types ! *** local variables *** INTEGER :: i ! Index of observation component INTEGER :: cnt ! Counter REAL :: inno2 ! Squared innovation REAL :: limit2 ! Squared limit ! ********************** ! *** INITIALIZATION *** ! ********************** doassim: IF (thisobs%doassim == 1) THEN IF (thisobs%inno_omit > 0.0) THEN IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_omit_by_inno -- START obs-ID', obsid WRITE (*,*) '++ OMI-debug omit_by_inno:', debug, 'limit for innovation', & thisobs%inno_omit WRITE (*,*) '++ OMI-debug omit_by_inno:', debug, 'inno_omit_ivar', & thisobs%inno_omit_ivar WRITE (*,*) '++ OMI-debug omit_by_inno:', debug, 'innovation_f(1:10)', inno_f(1:10) ENDIF ! Squared limit factor limit2 = thisobs%inno_omit * thisobs%inno_omit ! Check for observations to be excluded cnt = 0 DO i = 1, thisobs%dim_obs_f ! Squared innovation inno2 = inno_f(i + thisobs%off_obs_f)* inno_f(i + thisobs%off_obs_f) IF (inno2 > limit2 * 1.0/thisobs%ivar_obs_f(i)) THEN IF (debug>0) THEN WRITE (*,*) '++ OMI-debug omit_by_inno:', debug, 'omit: innovation:', & inno_f(i + thisobs%off_obs_f), 'observation:', obs_f_all(i + thisobs%off_obs_f) END IF ! Exclude observation by increased its observation error thisobs%ivar_obs_f(i) = thisobs%inno_omit_ivar ! Count excluded obs cnt = cnt + 1 END IF ENDDO IF (debug>0 .and. cnt>0) THEN WRITE (*,*) '++ OMI-debug omit_by_inno:', debug, 'count of excluded obs.: ', cnt WRITE (*,*) '++ OMI-debug omit_by_inno:', debug, 'updated thisobs_f%ivar_obs_f ', & thisobs%ivar_obs_f ENDIF IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_omit_by_inno_f -- END obs-ID', obsid cnt_all = cnt_all + cnt END IF ENDIF doassim END SUBROUTINE PDAFomi_omit_by_inno !------------------------------------------------------------------------------- !> Get statistics on local observations !! !! The routine is called in the update routine of !! global filters and writes statistics on !! used and excluded observations. !! !! __Revision history:__ !! * 2023-03 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obsstats(screen) USE MPI USE PDAF_mod_filtermpi, & ONLY: COMM_filter, MPIerr, npes_filter IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: screen !< Verbosity flag ! *** Local variables *** INTEGER :: ostats_omit_g(7) IF (npes_filter>1) THEN CALL MPI_Reduce(ostats_omit, ostats_omit_g, 7, MPI_INTEGER, MPI_SUM, & 0, COMM_filter, MPIerr) ELSE ! This is a work around for working with nullmpi.F90 ostats_omit_g = ostats_omit END IF IF (mype == 0 .AND. screen > 0 .AND. ostats_omit_g(1)>0) THEN WRITE (*, '(a, 9x, a)') 'PDAFomi', 'Global statistics for omitted observations:' WRITE (*, '(a, 11x, a, i10)') & 'PDAFomi', 'Global number of omitted observations: ', ostats_omit_g(6) WRITE (*, '(a, 11x, a, i10)') & 'PDAFomi', 'Global number of used observations: ', ostats_omit_g(7) ELSEIF (mype == 0 .AND. screen > 0) THEN WRITE (*, '(a, 9x, a)') 'PDAFomi', 'Global statistics for omitted observations:' WRITE (*, '(a, 11x, a)') & 'PDAFomi', 'Zero observations omitted' END IF END SUBROUTINE PDAFomi_obsstats !------------------------------------------------------------------------------- !> Gather global observation dimension information !! !! This routine gathers the information about !! the full dimension of each observation type !! in each process-local subdomain. !! SUBROUTINE PDAFomi_gather_obsdims() IMPLICIT NONE ! *** local variables *** INTEGER :: i ! Loop counter INTEGER :: dim_obs_all ! Full number of global observations ! ***************************************** ! *** Gather all observation dimensions *** ! ***************************************** ALLOCATE(obsdims(npes,n_obstypes)) DO i = 1, n_obstypes CALL MPI_Allgather(obs_f_all(i)%ptr%dim_obs_f, 1, MPI_INTEGER, obsdims(:,i), 1, & MPI_INTEGER, COMM_filter, MPIerr) END DO ! Determine overall number of observations dim_obs_all = SUM(obsdims) ! Allocate mapping vector ALLOCATE(map_obs_id(dim_obs_all)) END SUBROUTINE PDAFomi_gather_obsdims !------------------------------------------------------------------------------- !> Check error flag !! !! This routine returns the value of the PDAF-OMI internal error flag. !! !! __Revision history:__ !! * 2024-06 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_check_error(flag) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(inout) :: flag !< Error flag ! Set error flag flag = error END SUBROUTINE PDAFomi_check_error END MODULE PDAFomi_obs_f