! 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_l.F90 1147 2023-03-12 16:14:34Z lnerger $ !> PDAF-OMI routines for local observations !! !! This module contains generic routines for several observation-related !! operations for local filters. The routines are !! !! * PDAFomi_set_debug_flag \n !! Set or unset the debugging flag for PDAFomi routines !! * PDAFomi_init_dim_obs_l \n !! Initialize dimension of local obs. vetor and arrays for !! local observations !! * PDAFomi_cnt_dim_obs_l \n !! Set dimension of local obs. vector with isotropic localization !! * PDAFomi_cnt_dim_obs_l_noniso \n !! Set dimension of local obs. vector with nonisotropic localization !! * PDAFomi_init_obsarrays_l \n !! Initialize arrays for the index of a local observation in !! the full observation vector and its corresponding distance. !! * PDAFomi_init_obsarrays_l_noniso \n !! Initialize arrays for the index of a local observation in !! the full observation vector and its corresponding distance !! with onoisotrppic localization. !! * PDAFomi_g2l_obs \n !! Initialize local observation vector from full observation vector !! * PDAFomi_init_obs_l \n !! Initialize the local vector of observations !! * PDAFomi_prodRinvA_l \n !! Multiply an intermediate matrix of the local filter analysis !! with the inverse of the observation error covariance matrix !! and apply observation localization !! * PDAFomi_prodRinvA_hyb_l \n !! Multiply an intermediate matrix of the local filter analysis !! with the inverse of the observation error covariance matrix !! and apply observation localization. In addition apply the !! hybrid weight !! * PDAFomi_init_obsvar_l \n !! Compute mean observation error variance !! * PDAFomi_likelihood_l \n !! Compute local likelihood for an ensemble member !! * PDAFomi_likelihood_hyb_l \n !! Compute local likelihood for an ensemble member taking into !! account a hybrid weight for tempering !! * PDAFomi_localize_covar_iso \n !! Apply covariance isotropic localization in LEnKF !! * PDAFomi_localize_covar_noniso \n !! Apply non-isotropic covariance localization in LEnKF !! * PDAFomi_g2l_obs_internal \n !! Internal routine to initialize local observation vector from full !! observation vector (used by PDAFomi_init_obs_l and PDAFomi_g2l_obs) !! * PDAFomi_comp_dist2 \n !! Compute squared distance !! * PDAFomi_check_dist2 \n !! Compute and check distance for isotropic localization !! * PDAFomi_check_dist2_noniso \n !! Compute and check distance for non-isotropic localization !! * PDAFomi_weights_l \n !! Compute a vector of localization weights !! * PDAFomi_deallocate_obs \n !! Deallocate arrays in observation type !! * PDAFomi_dealloc \n !! Deallocate arrays in all observation types !! * PDAFomi_omit_by_innovation_l \ !! Exclude observations if innovation is too large (thisobs%inno_exclude) !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code !! * Later revisions - see repository log !! MODULE PDAFomi_obs_l USE PDAFomi_obs_f, ONLY: obs_f, r_earth, pi, debug, n_obstypes, error USE PDAF_mod_filter, ONLY: screen, obs_member USE PDAF_mod_filtermpi, ONLY: mype, npes_filter IMPLICIT NONE SAVE ! *** Module internal variables ! Data type to define the local observations by internally shared variables of the module TYPE obs_l INTEGER :: dim_obs_l !< number of local observations INTEGER :: off_obs_l !< Offset of this observation in overall local obs. vector INTEGER, ALLOCATABLE :: id_obs_l(:) !< Indices of local observations in full obs. vector REAL, ALLOCATABLE :: distance_l(:) !< Distances of local observations REAL, ALLOCATABLE :: cradius_l(:) !< directional cut-off radii of local observations REAL, ALLOCATABLE :: sradius_l(:) !< directional support radii of local observations REAL, ALLOCATABLE :: ivar_obs_l(:) !< Inverse variance of local observations REAL, ALLOCATABLE :: dist_l_v(:) !< Vertical distances of local observations INTEGER :: locweight !< Specify localization function INTEGER :: locweight_v=0 !< Specify localization function for vertical direction INTEGER :: nradii !< Length of CRADIUS and SRADIUS REAL, ALLOCATABLE :: cradius(:) !< Localization cut-off radius (single value or vector) REAL, ALLOCATABLE :: sradius(:) !< support radius for localization function (single value or vector) END TYPE obs_l TYPE obs_arr_l ! Type for pointer array over all observation types TYPE(obs_l), POINTER :: ptr END TYPE obs_arr_l TYPE(obs_arr_l), ALLOCATABLE :: obs_l_all(:) ! Declare pointer array INTEGER :: firstobs = 0 ! Flag for very first call to init_dim_obs_l INTEGER :: offset_obs_l = 0 ! offset of current observation in overall local obs. vector !$OMP THREADPRIVATE(obs_l_all, firstobs, offset_obs_l) INTERFACE PDAFomi_init_dim_obs_l MODULE PROCEDURE PDAFomi_init_dim_obs_l_iso MODULE PROCEDURE PDAFomi_init_dim_obs_l_noniso MODULE PROCEDURE PDAFomi_init_dim_obs_l_noniso_locweights END INTERFACE INTERFACE PDAFomi_localize_covar MODULE PROCEDURE PDAFomi_localize_covar_iso MODULE PROCEDURE PDAFomi_localize_covar_noniso END INTERFACE !------------------------------------------------------------------------------- CONTAINS !!> Set debugging flag !! !! This routine sets the debug flag for PDAF-OMI. !! One can set the flag dependent on the local analysis !! domain, the MPI rank, or the OpenMP thread ID, or !! and combination of them. !! !! For debugval>0 additional information is written by !! the OMI routine to stdout. One should activate the !! debugging before calling some selected routine(s) and !! deactivate it with debugval=0 afterwards. This allows !! for a targeted checking of the functionality. !! !! __Revision history:__ !! * 2019-09 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_set_debug_flag(debugval) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: debugval !< Value for debugging flag debug = debugval ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug set_debug_flag: mype_filter', mype, 'activate', debug END IF END SUBROUTINE PDAFomi_set_debug_flag !------------------------------------------------------------------------------- !> Set dimension of local obs. vector and local obs. arrays !! !! This routine sets the number of local observations for the !! current observation type for the local analysis domain !! with coordinates COORD_l and localization cut-off radius CRADIUS. !! Further the routine initializes arrays for the index of a !! local observation in the full observation vector and its !! corresponding distance. !! The operation are performed by calling the routines !! cnt_dim_obs_l and init_obsarrays_l. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_dim_obs_l_iso(thisobs_l, thisobs, coords_l, locweight, cradius, & sradius, cnt_obs_l) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation TYPE(obs_l), TARGET, INTENT(inout) :: thisobs_l !< Data type with local observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current analysis domain INTEGER, INTENT(in) :: locweight !< Type of localization function REAL, INTENT(in) :: cradius !< Localization cut-off radius (single or vector) REAL, INTENT(in) :: sradius !< Support radius of localization function (single or vector) INTEGER, INTENT(inout) :: cnt_obs_l !< Local dimension of current observation vector ! *** Local variables *** REAL :: maxcoords_l, mincoords_l ! Min/Max domain coordinates to check geographic coords REAL :: maxocoords_l, minocoords_l ! Min/Max observation coordinates to check geographic coords doassim: IF (thisobs%doassim == 1) THEN ! *********************************************** ! *** Check offset in full observation vector *** ! *********************************************** IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_dim_obs_l -- START' IF (thisobs%ncoord/=3 .AND. thisobs%disttype>=10) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: factorized 2+1D localization can only be used for thisobs%ncoord=3' error = 14 END IF ! Store ID of first observation type that call the routine ! This is reset in PDAFomi_deallocate_obs IF (firstobs == 0) THEN firstobs = thisobs%obsid END IF ! Reset offset of currrent observation in overall local obs. vector IF (thisobs%obsid == firstobs) THEN offset_obs_l = 0 cnt_obs_l = 0 END IF ! ************************************** ! *** Store localization information *** ! ************************************** thisobs_l%locweight = locweight IF (ALLOCATED(thisobs_l%cradius)) DEALLOCATE(thisobs_l%cradius) ALLOCATE(thisobs_l%cradius(1)) IF (ALLOCATED(thisobs_l%sradius)) DEALLOCATE(thisobs_l%sradius) ALLOCATE(thisobs_l%sradius(1)) thisobs_l%nradii = 1 thisobs_l%cradius(1) = cradius thisobs_l%sradius(1) = sradius ! Store offset thisobs_l%off_obs_l = offset_obs_l ! ************************************************** ! *** Initialize local observation pointer array *** ! ************************************************** ! Initialize pointer array IF (thisobs%obsid == firstobs) THEN IF (ALLOCATED(obs_l_all)) DEALLOCATE(obs_l_all) ALLOCATE(obs_l_all(n_obstypes)) END IF ! Set pointer to current observation obs_l_all(thisobs%obsid)%ptr => thisobs_l ! ********************************************** ! *** Initialize local observation dimension *** ! ********************************************** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_init_dim_obs_l -- count local observations' IF (thisobs%obsid == firstobs) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, ' Re-init dim_obs_l=0' END IF WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, ' coords_l', coords_l ! For geographic coordinates check whether their range is reasonable IF (thisobs%disttype==2 .OR. thisobs%disttype==3 .OR. thisobs%disttype==12 .OR. thisobs%disttype==13) THEN maxcoords_l = MAXVAL(coords_l) mincoords_l = MINVAL(coords_l) maxocoords_l = MAXVAL(thisobs%ocoord_f(1:2, :)) minocoords_l = MINVAL(thisobs%ocoord_f(1:2, :)) IF (maxcoords_l>2.0*pi .OR. mincoords_l<-pi .OR. maxocoords_l>2.0*pi .OR. minocoords_l<-pi) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, & ' WARNING: The unit for geographic coordinates is radian, thus range (0,2*pi) or (-pi,pi)!' END IF END IF WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, & ' Note: Please ensure that coords_l and observation coordinates have the same unit' END IF CALL PDAFomi_cnt_dim_obs_l(thisobs_l, thisobs, coords_l) ! Store number of local module-type observations for output cnt_obs_l = cnt_obs_l + thisobs_l%dim_obs_l ! ************************************************************ ! *** Initialize internal local arrays for local distances *** ! *** and indices of local obs. in full obs. vector *** ! ************************************************************ IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_init_dim_obs_l -- initialize local observation arrays' ! Initialize ID_OBS_L and DISTANCE_L and increment offsets CALL PDAFomi_init_obsarrays_l(thisobs_l, thisobs, coords_l, offset_obs_l) ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, ' thisobs_l%dim_obs_l', thisobs_l%dim_obs_l IF (thisobs_l%dim_obs_l>0) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, ' thisobs_l%id_obs_l', thisobs_l%id_obs_l WRITE (*,*) '++ OMI-debug init_dim_obs_l:', debug, ' thisobs_l%distance_l', thisobs_l%distance_l END IF WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_dim_obs_l -- END' END IF ELSE doassim cnt_obs_l = cnt_obs_l + 0 END IF doassim END SUBROUTINE PDAFomi_init_dim_obs_l_iso !------------------------------------------------------------------------------- !> Set dimension of local obs. vector and local obs. arrays (non-isotropic) !! !! This routine sets the number of local observations for the !! current observation type for the local analysis domain !! with coordinates COORD_l and a vector of localization cut-off !! radii CRADIUS. !! Further the routine initializes arrays for the index of a !! local observation in the full observation vector and its !! corresponding distance. !! The operation are performed by calling the routines !! cnt_dim_obs_l and init_obsarrays_l. !! !! __Revision history:__ !! * 2024-02 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_dim_obs_l_noniso(thisobs_l, thisobs, coords_l, locweight, cradius, & sradius, cnt_obs_l) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation TYPE(obs_l), TARGET, INTENT(inout) :: thisobs_l !< Data type with local observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current analysis domain INTEGER, INTENT(in) :: locweight !< Type of localization function REAL, INTENT(in) :: cradius(:) !< Vector of localization cut-off radii REAL, INTENT(in) :: sradius(:) !< Vector of support radii of localization function INTEGER, INTENT(inout) :: cnt_obs_l !< Local dimension of current observation vector ! *** Local variables *** REAL :: maxcoords_l, mincoords_l ! Min/Max domain coordinates to check geographic coords REAL :: maxocoords_l, minocoords_l ! Min/Max observation coordinates to check geographic coords doassim: IF (thisobs%doassim == 1) THEN ! *********************************************** ! *** Check offset in full observation vector *** ! *********************************************** IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_dim_obs_l_noniso -- START' ! Check consistency of dimensions IF (SIZE(cradius) /= thisobs%ncoord) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: non-isotropic localization: Size of CRADIUS /= thisobs%ncoord' error = 12 END IF IF (SIZE(sradius) /= thisobs%ncoord) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: non-isotropic localization: Size of SRADIUS /= thisobs%ncoord' error = 13 END IF IF (thisobs%ncoord/=3 .AND. thisobs%disttype>=10) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: factorized 2+1D localization can only be used for thisobs%ncoord=3' error = 14 END IF ! Store ID of first observation type that call the routine ! This is reset in PDAFomi_deallocate_obs IF (firstobs == 0) THEN firstobs = thisobs%obsid END IF ! Reset offset of currrent observation in overall local obs. vector IF (thisobs%obsid == firstobs) THEN offset_obs_l = 0 cnt_obs_l = 0 END IF ! ************************************** ! *** Store localization information *** ! ************************************** thisobs_l%locweight = locweight ! Allocate vectors for localization radii and store their values IF (ALLOCATED(thisobs_l%cradius)) DEALLOCATE(thisobs_l%cradius) ALLOCATE(thisobs_l%cradius(thisobs%ncoord)) IF (ALLOCATED(thisobs_l%sradius)) DEALLOCATE(thisobs_l%sradius) ALLOCATE(thisobs_l%sradius(thisobs%ncoord)) thisobs_l%nradii = thisobs%ncoord thisobs_l%cradius(:) = cradius(:) thisobs_l%sradius(:) = sradius(:) ! Store offset thisobs_l%off_obs_l = offset_obs_l ! ************************************************** ! *** Initialize local observation pointer array *** ! ************************************************** ! Initialize pointer array IF (thisobs%obsid == firstobs) THEN IF (ALLOCATED(obs_l_all)) DEALLOCATE(obs_l_all) ALLOCATE(obs_l_all(n_obstypes)) END IF ! Set pointer to current observation obs_l_all(thisobs%obsid)%ptr => thisobs_l ! ********************************************** ! *** Initialize local observation dimension *** ! ********************************************** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_init_dim_obs_l_noniso -- count local observations' IF (thisobs%obsid == firstobs) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, ' Re-init dim_obs_l=0' END IF WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, ' coords_l', coords_l ! For geographic coordinates check whether their range is reasonable IF (thisobs%disttype==2 .OR. thisobs%disttype==3 .OR. thisobs%disttype==12 .OR. thisobs%disttype==13) THEN maxcoords_l = MAXVAL(coords_l) mincoords_l = MINVAL(coords_l) maxocoords_l = MAXVAL(thisobs%ocoord_f(1:2, :)) minocoords_l = MINVAL(thisobs%ocoord_f(1:2, :)) IF (maxcoords_l>2.0*pi .OR. mincoords_l<-pi .OR. maxocoords_l>2.0*pi .OR. minocoords_l<-pi) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, & ' WARNING: The unit for geographic coordinates is radian, thus range (0,2*pi) or (-pi,pi)!' END IF END IF WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, & ' Note: Please ensure that coords_l and observation coordinates have the same unit' END IF CALL PDAFomi_cnt_dim_obs_l_noniso(thisobs_l, thisobs, coords_l) ! Store number of local module-type observations for output cnt_obs_l = cnt_obs_l + thisobs_l%dim_obs_l ! ************************************************************ ! *** Initialize internal local arrays for local distances *** ! *** and indices of local obs. in full obs. vector *** ! ************************************************************ IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_init_dim_obs_l_noniso -- initialize local observation arrays' ! Initialize ID_OBS_L and DISTANCE_L and increment offsets CALL PDAFomi_init_obsarrays_l_noniso(thisobs_l, thisobs, coords_l, offset_obs_l) ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, ' thisobs_l%dim_obs_l', thisobs_l%dim_obs_l IF (thisobs_l%dim_obs_l>0) THEN WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, ' thisobs_l%id_obs_l', thisobs_l%id_obs_l WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso:', debug, ' thisobs_l%distance_l', thisobs_l%distance_l END IF WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_dim_obs_l_noniso -- END' END IF ELSE doassim cnt_obs_l = cnt_obs_l + 0 END IF doassim END SUBROUTINE PDAFomi_init_dim_obs_l_noniso !------------------------------------------------------------------------------- !> Set dimension of local obs. vector and local obs. arrays !! !! This routine is a variant of PDAFomi_init_dim_obs_l_noniso with !! support for a vector of localization weights. This is used !! to specify different localization functions for the vertical and !! horizontal directions. The routine only stores the value of !! locweights(2) for the vertical and calls PDAFomi_init_dim_obs_l_iso. !! !! __Revision history:__ !! * 2024-04 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_dim_obs_l_noniso_locweights(thisobs_l, thisobs, coords_l, locweights, cradius, & sradius, cnt_obs_l) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation TYPE(obs_l), TARGET, INTENT(inout) :: thisobs_l !< Data type with local observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current analysis domain INTEGER, INTENT(in) :: locweights(:) !< Types of localization function REAL, INTENT(in) :: cradius(:) !< Vector of localization cut-off radii REAL, INTENT(in) :: sradius(:) !< Vector of support radii of localization function INTEGER, INTENT(inout) :: cnt_obs_l !< Local dimension of current observation vector ! *** Store vertical locweight and call standard routine IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_dim_obs_l_noniso_locweights -- START' WRITE (*,*) '++ OMI-debug init_dim_obs_l_noniso_locweights:', debug, ' locweights', locweights END IF ! Check consistency of dimensions IF (SIZE(locweights) /= 2) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: Input for locweights in horizontal and vertical needs size 2' error = 15 END IF IF (thisobs%ncoord /= 3) THEN WRITE (*,*) '+++++ WARNING PDAF-OMI: separate locweight for vertical is only utilized if thisobs%ncoord=3' END IF IF (thisobs%ncoord == 3) THEN ! locweight for the vertical is treated separately thisobs_l%locweight_v = locweights(2) END IF ! Call to usual routine that handles a single locweight setting CALL PDAFomi_init_dim_obs_l_noniso(thisobs_l, thisobs, coords_l, locweights(1), cradius, & sradius, cnt_obs_l) IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_dim_obs_l_noniso_locweights -- END' END SUBROUTINE PDAFomi_init_dim_obs_l_noniso_locweights !------------------------------------------------------------------------------- !> Set dimension of local observation vector !! !! This routine sets the number of local observations for the !! current observation type for the local analysis domain !! with coordinates COORDS_L and localization cut-off radius CRADIUS. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_cnt_dim_obs_l(thisobs_l, thisobs, coords_l) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current analysis domain (thisobs%ncoord) ! *** Local variables *** INTEGER :: i, cnt ! Counters REAL :: cradius ! localization cut-off radius REAL :: distance2 ! squared distance REAL :: sradius ! support radius LOGICAL :: checkdist ! Flag whether distance nis not larger than cut-off radius ! ********************************************** ! *** Initialize local observation dimension *** ! ********************************************** ! Count local observations thisobs_l%dim_obs_l = 0 cnt = 0 IF (debug>0) THEN WRITE (*,*) '++ OMI-debug cnt_dim_obs_l: ', debug, ' thisobs%ncoord', thisobs%ncoord WRITE (*,*) '++ OMI-debug cnt_dim_obs_l: ', debug, ' thisobs_l%cradius', thisobs_l%cradius WRITE (*,*) '++ OMI-debug cnt_dim_obs_l: ', debug, ' Check for observations within radius' END IF scancount: DO i = 1, thisobs%dim_obs_f CALL PDAFomi_check_dist2(thisobs, thisobs_l, coords_l, thisobs%ocoord_f(1:thisobs%ncoord, i), distance2, & checkdist, i-1, cnt) ! If distance below limit, add observation to local domain IF (checkdist) THEN IF (debug>0) THEN WRITE (*,*) '++ OMI-debug cnt_dim_obs_l: ', debug, & ' valid observation with coordinates', thisobs%ocoord_f(1:thisobs%ncoord, i) END IF thisobs_l%dim_obs_l = thisobs_l%dim_obs_l + 1 END IF END DO scancount END SUBROUTINE PDAFomi_cnt_dim_obs_l !------------------------------------------------------------------------------- !> Set dimension of local observation vector for nonisotropic localization !! !! This routine sets the number of local observations for the !! current observation type for the local analysis domain !! with coordinates COORDS_L and localization cut-off radius CRADIUS. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_cnt_dim_obs_l_noniso(thisobs_l, thisobs, coords_l) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current analysis domain (thisobs%ncoord) ! *** Local variables *** INTEGER :: i, cnt ! Counters REAL :: cradius ! localization cut-off radius REAL :: distance2 ! squared distance REAL :: sradius ! support radius LOGICAL :: checkdist ! Flag whether distance nis not larger than cut-off radius REAL :: dists(thisobs%ncoord) ! Distance vector between analysis point and observation ! ********************************************** ! *** Initialize local observation dimension *** ! ********************************************** ! Count local observations thisobs_l%dim_obs_l = 0 cnt = 0 IF (debug>0) THEN WRITE (*,*) '++ OMI-debug cnt_dim_obs_l_noniso: ', debug, ' thisobs%ncoord', thisobs%ncoord WRITE (*,*) '++ OMI-debug cnt_dim_obs_l_noniso: ', debug, ' thisobs_l%cradius', thisobs_l%cradius WRITE (*,*) '++ OMI-debug cnt_dim_obs_l_noniso: ', debug, ' Check for observations within radius' END IF IF (thisobs_l%nradii==1) THEN ! 1D but with radius specified as array scancount: DO i = 1, thisobs%dim_obs_f CALL PDAFomi_check_dist2(thisobs, thisobs_l, coords_l, thisobs%ocoord_f(1:thisobs%ncoord, i), distance2, & checkdist, i-1, cnt) ! If distance below limit, add observation to local domain IF (checkdist) THEN IF (debug>0) THEN WRITE (*,*) '++ OMI-debug cnt_dim_obs_l_noniso: ', debug, & ' valid observation with coordinates', thisobs%ocoord_f(1:thisobs%ncoord, i) END IF thisobs_l%dim_obs_l = thisobs_l%dim_obs_l + 1 END IF END DO scancount ELSEIF (thisobs_l%nradii==2 .OR. thisobs_l%nradii==3) THEN ! Nonisotropic in 2 or 3 dimensions scancountB: DO i = 1, thisobs%dim_obs_f CALL PDAFomi_check_dist2_noniso(thisobs, thisobs_l, coords_l, thisobs%ocoord_f(1:thisobs%ncoord, i), distance2, & dists, cradius, sradius, checkdist, i-1, cnt) ! If distance below limit, add observation to local domain IF (checkdist .AND. debug>0) THEN WRITE (*,*) '++ OMI-debug cnt_dim_obs_l_noniso: ', debug, & ' valid observation with coordinates', thisobs%ocoord_f(1:thisobs%ncoord, i) WRITE (*,*) '++ OMI-debug cnt_dim_obs_l_noniso: ', debug, & ' valid observation distance, cradius, sradius', SQRT(distance2), cradius, sradius END IF END DO scancountB thisobs_l%dim_obs_l = cnt ELSE WRITE (*,*) '+++++ ERROR PDAF-OMI: nonisotropic localization is only possible in 1, 2 or 3 dimensions' error = 10 END IF END SUBROUTINE PDAFomi_cnt_dim_obs_l_noniso !------------------------------------------------------------------------------- !> Initialize local arrays for an observation !! !! This routine has to initialize for the current !! observation type the indices of the local observations !! in the full observation vector and the corresponding !! distances from the local analysis domain. The offset !! of the observation type in the local onbservation !! vector is given by OFF_OBS_L_ALL. !! !! The routine has also to return OFF_OBS_L_ALL incremented !! by the number of initialized local observations. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obsarrays_l(thisobs_l, thisobs, coords_l, off_obs_l_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current water column (thisobs%ncoord) INTEGER, INTENT(inout) :: off_obs_l_all !< input: offset of current obs. in local obs. vector !< output: input + thisobs_l%dim_obs_l ! *** Local variables *** INTEGER :: i, off_obs ! Counters REAL :: cradius ! localization cut-off radius REAL :: distance2 ! squared distance REAL :: sradius ! support radius LOGICAL :: checkdist ! Flag whether distance nis not larger than cut-off radius ! ********************************************** ! *** Initialize local observation dimension *** ! ********************************************** ! Allocate module-internal index array for indices in module-type observation vector IF (ALLOCATED(thisobs_l%id_obs_l)) DEALLOCATE(thisobs_l%id_obs_l) IF (ALLOCATED(thisobs_l%distance_l)) DEALLOCATE(thisobs_l%distance_l) IF (ALLOCATED(thisobs_l%cradius_l)) DEALLOCATE(thisobs_l%cradius_l) IF (ALLOCATED(thisobs_l%sradius_l)) DEALLOCATE(thisobs_l%sradius_l) IF (thisobs_l%dim_obs_l>0) THEN ALLOCATE(thisobs_l%id_obs_l(thisobs_l%dim_obs_l)) ALLOCATE(thisobs_l%distance_l(thisobs_l%dim_obs_l)) ALLOCATE(thisobs_l%cradius_l(thisobs_l%dim_obs_l)) ALLOCATE(thisobs_l%sradius_l(thisobs_l%dim_obs_l)) ELSE ALLOCATE(thisobs_l%id_obs_l(1)) ALLOCATE(thisobs_l%distance_l(1)) ALLOCATE(thisobs_l%cradius_l(1)) ALLOCATE(thisobs_l%sradius_l(1)) END IF off_obs = 0 ! Count local observations and initialize index and distance arrays IF (thisobs_l%dim_obs_l>0) THEN scancount: DO i = 1, thisobs%dim_obs_f CALL PDAFomi_check_dist2(thisobs, thisobs_l, coords_l, thisobs%ocoord_f(1:thisobs%ncoord, i), distance2, & checkdist, i-1, off_obs) ! If distance below limit, add observation to local domain IF (checkdist) THEN ! For internal storage (use in prodRinvA_l) thisobs_l%id_obs_l(off_obs) = i ! node index thisobs_l%distance_l(off_obs) = SQRT(distance2) ! distance thisobs_l%cradius_l(off_obs) = thisobs_l%cradius(1) ! isotropic cut-off radius thisobs_l%sradius_l(off_obs) = thisobs_l%sradius(1) ! isotropic support radius END IF END DO scancount ! Count overall local observations off_obs_l_all = off_obs_l_all + off_obs ! dimension END IF END SUBROUTINE PDAFomi_init_obsarrays_l !------------------------------------------------------------------------------- !> Initialize local arrays for an observation for nonisotropic localization !! !! This routine has to initialize for the current !! observation type the indices of the local observations !! in the full observation vector and the corresponding !! distances from the local analysis domain. The offset !! of the observation type in the local onbservation !! vector is given by OFF_OBS_L_ALL. !! !! The routine has also to return OFF_OBS_L_ALL incremented !! by the number of initialized local observations. !! !! __Revision history:__ !! * 2024-04 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obsarrays_l_noniso(thisobs_l, thisobs, coords_l, off_obs_l_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: coords_l(:) !< Coordinates of current water column (thisobs%ncoord) INTEGER, INTENT(inout) :: off_obs_l_all !< input: offset of current obs. in local obs. vector !< output: input + thisobs_l%dim_obs_l ! *** Local variables *** INTEGER :: i, off_obs ! Counters REAL :: cradius ! localization cut-off radius REAL :: distance2 ! squared distance REAL :: sradius ! support radius LOGICAL :: checkdist ! Flag whether distance nis not larger than cut-off radius REAL :: dists(thisobs%ncoord) ! Distance vector between analysis point and observation ! ********************************************** ! *** Initialize local observation dimension *** ! ********************************************** ! Allocate module-internal index array for indices in module-type observation vector IF (ALLOCATED(thisobs_l%id_obs_l)) DEALLOCATE(thisobs_l%id_obs_l) IF (ALLOCATED(thisobs_l%distance_l)) DEALLOCATE(thisobs_l%distance_l) IF (ALLOCATED(thisobs_l%cradius_l)) DEALLOCATE(thisobs_l%cradius_l) IF (ALLOCATED(thisobs_l%sradius_l)) DEALLOCATE(thisobs_l%sradius_l) IF (thisobs_l%dim_obs_l>0) THEN ALLOCATE(thisobs_l%id_obs_l(thisobs_l%dim_obs_l)) ALLOCATE(thisobs_l%distance_l(thisobs_l%dim_obs_l)) ALLOCATE(thisobs_l%cradius_l(thisobs_l%dim_obs_l)) ALLOCATE(thisobs_l%sradius_l(thisobs_l%dim_obs_l)) ELSE ALLOCATE(thisobs_l%id_obs_l(1)) ALLOCATE(thisobs_l%distance_l(1)) ALLOCATE(thisobs_l%cradius_l(1)) ALLOCATE(thisobs_l%sradius_l(1)) END IF ! Allocate array for vertical distance in factorized localization (disttype>=10) IF (thisobs_l%locweight_v>0) THEN IF (ALLOCATED(thisobs_l%dist_l_v)) DEALLOCATE(thisobs_l%dist_l_v) IF (thisobs_l%dim_obs_l>0) THEN ALLOCATE(thisobs_l%dist_l_v(thisobs_l%dim_obs_l)) ELSE ALLOCATE(thisobs_l%dist_l_v(1)) END IF END IF off_obs = 0 ! Count local observations and initialize index and distance arrays IF (thisobs_l%nradii==1) THEN ! 1D but with radius specified as array IF (thisobs_l%dim_obs_l>0) THEN scancount: DO i = 1, thisobs%dim_obs_f CALL PDAFomi_check_dist2(thisobs, thisobs_l, coords_l, thisobs%ocoord_f(1:thisobs%ncoord, i), distance2, & checkdist, i-1, off_obs) ! If distance below limit, add observation to local domain IF (checkdist) THEN ! For internal storage (use in prodRinvA_l) thisobs_l%id_obs_l(off_obs) = i ! node index thisobs_l%distance_l(off_obs) = SQRT(distance2) ! distance thisobs_l%cradius_l(off_obs) = thisobs_l%cradius(1) ! isotropic cut-off radius thisobs_l%sradius_l(off_obs) = thisobs_l%sradius(1) ! isotropic support radius END IF END DO scancount ! Count overall local observations off_obs_l_all = off_obs_l_all + off_obs ! dimension END IF ELSEIF (thisobs_l%nradii==2 .OR. thisobs_l%nradii==3) THEN ! Nonisotropic in 2 or 3 dimensions IF (thisobs_l%dim_obs_l>0) THEN off_obs = 0 scancountB: DO i = 1, thisobs%dim_obs_f CALL PDAFomi_check_dist2_noniso(thisobs, thisobs_l, coords_l, thisobs%ocoord_f(1:thisobs%ncoord, i), distance2, & dists, cradius, sradius, checkdist, i-1, off_obs) ! If distance below limit, add observation to local domain IF (checkdist) THEN ! For internal storage (use in prodRinvA_l) thisobs_l%id_obs_l(off_obs) = i ! node index thisobs_l%distance_l(off_obs) = SQRT(distance2) ! distance thisobs_l%cradius_l(off_obs) = cradius ! directional cut-off radius thisobs_l%sradius_l(off_obs) = sradius ! directional support radius IF (thisobs_l%locweight_v>0 .AND. thisobs_l%nradii==3) THEN thisobs_l%dist_l_v(off_obs) = dists(3) ! distance in vertical direction END if END IF END DO scancountB ! Count overall local observations off_obs_l_all = off_obs_l_all + off_obs ! dimension END IF ELSE WRITE (*,*) '+++++ ERROR PDAF-OMI: nonisotropic localization is only possible in 1, 2 or 3 dimensions' error = 11 END IF END SUBROUTINE PDAFomi_init_obsarrays_l_noniso !------------------------------------------------------------------------------- !> Initialize local observation vector !! !! This routine has to initialize the part of the !! overall local observation vector corresponding !! to the current observation type. The offset of !! the current observation type in the local obs. !! vector is given by OFFSET_OBS_l_ALL. !! This routine uses a shortened interface and just !! passed the operation to the actually routine. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_g2l_obs(thisobs_l, thisobs, obs_f_all, obs_l_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: obs_f_all(:) !< Full obs. vector of current obs. for all variables REAL, INTENT(inout) :: obs_l_all(:) !< Local observation vector for all variables ! ******************************************* ! *** Initialize local observation vector *** ! ******************************************* doassim: IF (thisobs%doassim == 1) THEN IF (debug>0) THEN IF (obs_member==0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_g2l_obs -- START Get local observed ensemble mean' ELSE WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_g2l_obs -- START Get local observed ensemble member', obs_member END IF END IF CALL PDAFomi_g2l_obs_internal(thisobs_l, & obs_f_all(thisobs%off_obs_f+1:thisobs%off_obs_f+thisobs%dim_obs_f), & thisobs_l%off_obs_l, obs_l_all) IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_g2l_obs -- END' END IF doassim END SUBROUTINE PDAFomi_g2l_obs !------------------------------------------------------------------------------- !> Initialize local observation vector and inverse error variance !! !! This routine has to initialize the part of the !! overall local observation vector corresponding !! to the current observation type. The offset of !! the current observation type in the local obs. !! vector is given by OFFSET_OBS_l_ALL. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obs_l(thisobs_l, thisobs, obs_l_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(inout) :: obs_l_all(:) !< Local observation vector for all variables ! ******************************************* ! *** Initialize local observation vector *** ! ******************************************* doassim: IF (thisobs%doassim == 1) THEN IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_obs_l -- START' WRITE (*,*) '++ OMI-debug init_obs_l: ', debug, ' thisobs_l%dim_obs_l', thisobs_l%dim_obs_l WRITE (*,*) '++ OMI-debug init_obs_l: ', debug, ' thisobs_l%off_obs_l', thisobs_l%off_obs_l WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_init_obs_l -- Get local vector of observations' END IF ! Initialize local observations CALL PDAFomi_g2l_obs_internal(thisobs_l, thisobs%obs_f, thisobs_l%off_obs_l, obs_l_all) ! Initialize local inverse variances for current observation ! they will be used in prodRinva_l IF (ALLOCATED(thisobs_l%ivar_obs_l)) DEALLOCATE(thisobs_l%ivar_obs_l) IF (thisobs_l%dim_obs_l>0) THEN ALLOCATE(thisobs_l%ivar_obs_l(thisobs_l%dim_obs_l)) ELSE ALLOCATE(thisobs_l%ivar_obs_l(1)) END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & ' PDAFomi_init_obs_l -- Get local vector of inverse obs. variances' END IF CALL PDAFomi_g2l_obs_internal(thisobs_l, thisobs%ivar_obs_f, 0, thisobs_l%ivar_obs_l) ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_init_obs_l -- END' END IF END IF doassim END SUBROUTINE PDAFomi_init_obs_l !------------------------------------------------------------------------------- !> 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 the loop over all !! local analysis domains during each analysis !! by the routine PDAF_set_forget_local that !! estimates a local adaptive forgetting factor. !! The routine has to initialize the mean observation !! error variance for the current local analysis !! domain. (See init_obsvar_f for a global variant) !! !! The routine assumed 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 first the variance sum is computed by !! multiplying with the observation counter. !! !! __Revision history:__ !! * 2019-09 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_init_obsvar_l(thisobs_l, thisobs, meanvar_l, cnt_obs_l) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(inout) :: meanvar_l !< Mean variance INTEGER, INTENT(inout) :: cnt_obs_l !< Observation counter ! Local variables INTEGER :: i ! Counter ! *********************************** ! *** Compute local mean variance *** ! *********************************** doassim: IF (thisobs%doassim == 1) THEN IF (cnt_obs_l==0) THEN ! Reset mean variance meanvar_l = 0.0 ELSE ! Compute sum of variances from mean variance meanvar_l = meanvar_l * REAL(cnt_obs_l) END IF ! Add observation error variances DO i = 1, thisobs_l%dim_obs_l meanvar_l = meanvar_l + 1.0 / thisobs_l%ivar_obs_l(i) END DO ! Increment observation count cnt_obs_l = cnt_obs_l + thisobs_l%dim_obs_l ! Compute updated mean variance meanvar_l = meanvar_l / REAL(cnt_obs_l) ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug init_obsvar_l: ', debug, 'thisobs_l%dim_obs_l', thisobs_l%dim_obs_l WRITE (*,*) '++ OMI-debug init_obsvar_l: ', debug, 'cnt_obs_l', cnt_obs_l WRITE (*,*) '++ OMI-debug init_obsvar_l: ', debug, 'thisobs_l%ivar_obs_l', thisobs_l%ivar_obs_l(:) WRITE (*,*) '++ OMI-debug init_obsvar_l: ', debug, 'meanvar_l', meanvar_l END IF END IF doassim END SUBROUTINE PDAFomi_init_obsvar_l !------------------------------------------------------------------------------- !> Compute product of inverse of R with some matrix !! !! The routine is called during the analysis step !! on each local analysis domain. It has to !! compute the product of the inverse of the local !! observation error covariance matrix with !! the matrix of locally observed ensemble !! perturbations. !! !! Next to computing the product, a localizing !! weighting ("observation localization") can be !! applied to matrix A. !! !! This implementation assumes a diagonal observation !! error covariance matrix, and supports varying !! observation error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_prodRinvA_l(thisobs_l, thisobs, nobs_all, ncols, & A_l, C_l, verbose) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nobs_all !< Dimension of local obs. vector (all obs. types) INTEGER, INTENT(in) :: ncols !< Rank of initial covariance matrix REAL, INTENT(inout) :: A_l(:, :) !< Input matrix (thisobs_l%dim_obs_l, ncols) REAL, INTENT(out) :: C_l(:, :) !< Output matrix (thisobs_l%dim_obs_l, ncols) INTEGER, INTENT(in) :: verbose !< Verbosity flag ! *** local variables *** INTEGER :: i, j ! Index of observation component REAL, ALLOCATABLE :: weight(:) ! Localization weights REAL, ALLOCATABLE :: weight_v(:) ! Localization weights for vertical (for locweight_v>0) INTEGER :: idummy ! Dummy to access nobs_all INTEGER :: off ! row offset in A_l and C_l ! ********************** ! *** INITIALIZATION *** ! ********************** doassim: IF (thisobs%doassim == 1) THEN ! Initialize dummy to prevent compiler warning idummy = nobs_all ! Initialize offset off = thisobs_l%off_obs_l ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_prodrinva_l -- START Multiply with inverse R and and apply localization' WRITE (*,*) '++ OMI-debug prodrinva_l: ', debug, ' thisobs_l%locweight', thisobs_l%locweight IF (thisobs_l%locweight_v>0) & WRITE (*,*) '++ OMI-debug prodrinva_l: ', debug, ' thisobs_l%locweight_v', thisobs_l%locweight_v WRITE (*,*) '++ OMI-debug prodRinvA_l: ', debug, ' thisobs%dim_obs_l', thisobs_l%dim_obs_l WRITE (*,*) '++ OMI-debug prodRinvA_l: ', debug, ' thisobs%ivar_obs_l', thisobs_l%ivar_obs_l WRITE (*,*) '++ OMI-debug prodRinvA_l: ', debug, ' Input matrix A_l', A_l END IF IF (verbose == 1) THEN WRITE (*, '(a, 5x, a, 1x, i3)') & 'PDAFomi', '--- Domain localization for obs. type ID',thisobs%obsid IF (thisobs_l%nradii==1) THEN IF (thisobs%disttype<10) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- isotropic localization' ELSE WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- isotropic localization factorized in 2+1 dimensions' END IF WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Localization cut-off radius', thisobs_l%cradius WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Support radius', thisobs_l%sradius ELSE IF (thisobs_l%nradii==2) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- non-isotropic localization' WRITE (*, '(a, 8x, a, 1x, 2es11.3)') & 'PDAFomi', '--- Localization cut-off radii', thisobs_l%cradius WRITE (*, '(a, 8x, a, 1x, 2es11.3)') & 'PDAFomi', '--- Support radii', thisobs_l%sradius ELSE IF (thisobs_l%nradii==3) THEN IF (thisobs%disttype<10) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- non-isotropic localization in 3 dimensions' ELSE WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- non-isotropic localization factorized in 2+1 dimensions' END IF IF (thisobs_l%locweight_v>0) & WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- use separate localization function in vertical direction' WRITE (*, '(a, 8x, a, 1x, 3es11.3)') & 'PDAFomi', '--- Localization cut-off radii', thisobs_l%cradius WRITE (*, '(a, 8x, a, 1x, 3es11.3)') & 'PDAFomi', '--- Support radii', thisobs_l%sradius END IF ENDIF ! *********************************************** ! *** Apply a weight matrix with correlations *** ! *** of compact support to matrix A or the *** ! *** observation error covariance matrix. *** ! *********************************************** ! *** Initialize weight array ALLOCATE(weight(thisobs_l%dim_obs_l)) CALL PDAFomi_weights_l(verbose, thisobs_l%dim_obs_l, ncols, thisobs_l%locweight, & thisobs_l%cradius_l, thisobs_l%sradius_l, & A_l, thisobs_l%ivar_obs_l, thisobs_l%distance_l, weight) ! *** For factorized 2+1D localization use product of horizontal and vertical weights IF (thisobs_l%locweight_v>0) then IF (verbose == 1) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- initialize also weight function for vertical direction' END IF ALLOCATE(weight_v(thisobs_l%dim_obs_l)) CALL PDAFomi_weights_l_sgnl(verbose, thisobs_l%dim_obs_l, ncols, thisobs_l%locweight_v, & thisobs_l%cradius(3), thisobs_l%sradius(3), & A_l, thisobs_l%ivar_obs_l, thisobs_l%dist_l_v, weight_v) DO i = 1, thisobs_l%dim_obs_l weight(i) = weight(i) * weight_v(i) END DO DEALLOCATE(weight_v) END IF ! *** Handling of special weighting types *** lw2: IF (thisobs_l%locweight == 26) THEN ! Use square-root of 5th-order polynomial on A DO i = 1, thisobs_l%dim_obs_l ! Check if weight >0 (Could be <0 due to numerical precision) IF (weight(i) > 0.0) THEN weight(i) = SQRT(weight(i)) ELSE weight(i) = 0.0 END IF END DO END IF lw2 ! *** Apply weight doweighting: IF (thisobs_l%locweight >= 11) THEN ! *** Apply weight to matrix A DO j = 1, ncols DO i = 1, thisobs_l%dim_obs_l A_l(i+off, j) = weight(i) * A_l(i+off, j) END DO END DO ! *** -1 ! *** C = R A DO j = 1, ncols DO i = 1, thisobs_l%dim_obs_l C_l(i+off, j) = thisobs_l%ivar_obs_l(i) * A_l(i+off, j) END DO END DO ELSE doweighting ! *** Apply weight to matrix R only DO j = 1, ncols DO i = 1, thisobs_l%dim_obs_l C_l(i+off, j) = thisobs_l%ivar_obs_l(i) * weight(i) * A_l(i+off, j) END DO END DO END IF doweighting ! *** Clean up *** DEALLOCATE(weight) IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_prodrinva_l -- END' ENDIF doassim END SUBROUTINE PDAFomi_prodRinvA_l !------------------------------------------------------------------------------- !> Compute product of inverse of R with some matrix and hybrid weight !! !! The routine is called during the analysis step !! on each local analysis domain. It has to !! compute the product of the inverse of the local !! observation error covariance matrix with !! the matrix of locally observed ensemble !! perturbations. !! !! Next to computing the product, a localizing !! weighting ("observation localization") can be !! applied to matrix A. In addition the hybrid !! weight alpha is applied. !! !! This implementation assumes a diagonal observation !! error covariance matrix, and supports varying !! observation error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. !! !! __Revision history:__ !! * 2022-03 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_prodRinvA_hyb_l(thisobs_l, thisobs, nobs_all, ncols, & gamma, A_l, C_l, verbose) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: nobs_all !< Dimension of local obs. vector (all obs. types) INTEGER, INTENT(in) :: ncols !< Rank of initial covariance matrix REAL, INTENT(in) :: gamma !< Hybrid weight REAL, INTENT(inout) :: A_l(:, :) !< Input matrix (thisobs_l%dim_obs_l, ncols) REAL, INTENT(out) :: C_l(:, :) !< Output matrix (thisobs_l%dim_obs_l, ncols) INTEGER, INTENT(in) :: verbose !< Verbosity flag ! *** local variables *** INTEGER :: i, j ! Index of observation component REAL, ALLOCATABLE :: weight(:) ! Localization weights REAL, ALLOCATABLE :: weight_v(:) ! Localization weights for vertical (for locweight_v>0) INTEGER :: idummy ! Dummy to access nobs_all INTEGER :: off ! row offset in A_l and C_l ! ********************** ! *** INITIALIZATION *** ! ********************** doassim: IF (thisobs%doassim == 1) THEN ! Initialize dummy to prevent compiler warning idummy = nobs_all ! Initialize offset off = thisobs_l%off_obs_l ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_prodrinva_hyb_l -- START Multiply with inverse R and and apply localization' WRITE (*,*) '++ OMI-debug prodrinva_hyb_l: ', debug, ' thisobs_l%locweight', thisobs_l%locweight WRITE (*,*) '++ OMI-debug prodrinva_hyb_l: ', debug, 'thisobs%dim_obs_f', thisobs_l%dim_obs_l WRITE (*,*) '++ OMI-debug prodrinva_hyb_l: ', debug, 'thisobs%ivar_obs_f', thisobs_l%ivar_obs_l WRITE (*,*) '++ OMI-debug prodrinva_hyb_l: ', debug, 'Input matrix A_l', A_l END IF IF (verbose == 1) THEN WRITE (*,'(a, 5x, a, f12.5)') 'PDAFomi', '--- hybrid gamma=', gamma WRITE (*, '(a, 5x, a, 1x)') & 'PDAFomi', '--- Domain localization' WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Localization cut-off radius', thisobs_l%cradius WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Support radius', thisobs_l%sradius ENDIF ! *********************************************** ! *** Apply a weight matrix with correlations *** ! *** of compact support to matrix A or the *** ! *** observation error covariance matrix. *** ! *********************************************** ! *** Initialize weight array ALLOCATE(weight(thisobs_l%dim_obs_l)) CALL PDAFomi_weights_l(verbose, thisobs_l%dim_obs_l, ncols, thisobs_l%locweight, & thisobs_l%cradius_l, thisobs_l%sradius_l, & A_l, thisobs_l%ivar_obs_l, thisobs_l%distance_l, weight) ! *** For factorized 2+1D localization use product of horizontal and vertical weights IF (thisobs_l%locweight_v>0) then IF (verbose == 1) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- initialize also weight function for vertical direction' END IF ALLOCATE(weight_v(thisobs_l%dim_obs_l)) CALL PDAFomi_weights_l_sgnl(verbose, thisobs_l%dim_obs_l, ncols, thisobs_l%locweight_v, & thisobs_l%cradius(3), thisobs_l%sradius(3), & A_l, thisobs_l%ivar_obs_l, thisobs_l%dist_l_v, weight_v) DO i = 1, thisobs_l%dim_obs_l weight(i) = weight(i) * weight_v(i) END DO DEALLOCATE(weight_v) END IF ! *** Handling of special weighting types *** lw2: IF (thisobs_l%locweight == 26) THEN ! Use square-root of 5th-order polynomial on A DO i = 1, thisobs_l%dim_obs_l ! Check if weight >0 (Could be <0 due to numerical precision) IF (weight(i) > 0.0) THEN weight(i) = SQRT(weight(i)) ELSE weight(i) = 0.0 END IF END DO END IF lw2 ! *** Apply weight doweighting: IF (thisobs_l%locweight >= 11) THEN ! *** Apply weight to matrix A DO j = 1, ncols DO i = 1, thisobs_l%dim_obs_l A_l(i+off, j) = weight(i) * A_l(i+off, j) END DO END DO ! *** -1 ! *** C = R A DO j = 1, ncols DO i = 1, thisobs_l%dim_obs_l C_l(i+off, j) = gamma * thisobs_l%ivar_obs_l(i) * A_l(i+off, j) END DO END DO ELSE doweighting ! *** Apply weight to matrix R only DO j = 1, ncols DO i = 1, thisobs_l%dim_obs_l C_l(i+off, j) = gamma * thisobs_l%ivar_obs_l(i) * weight(i) * A_l(i+off, j) END DO END DO END IF doweighting ! *** Clean up *** DEALLOCATE(weight) IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_prodrinva_hyb_l -- END' ENDIF doassim END SUBROUTINE PDAFomi_prodRinvA_hyb_l !------------------------------------------------------------------------------- !> Compute local likelihood for an ensemble member !! !! The routine is called during the analysis step !! of the localized NETF. !! 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 addition, a localizing weighting of the !! inverse of R by expotential decrease or a 5-th order !! polynomial of compact support can be applied. This is !! defined by the variables 'locweight', 'cradius', !! 'cradius2' and 'sradius' in the main program. !! !! In general this routine is similar to the routine !! prodRinvA_l 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. !! !! This implementation assumes a diagonal observation !! error covariance matrix, and supports varying !! observation error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_likelihood_l(thisobs_l, thisobs, resid_l, lhood_l, verbose) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(inout) :: resid_l(:) !< Input vector of residuum REAL, INTENT(inout) :: lhood_l !< Output vector - log likelihood INTEGER, INTENT(in) :: verbose !< Verbosity flag ! *** local variables *** INTEGER :: i ! Index of observation component REAL, ALLOCATABLE :: weight(:) ! Localization weights REAL, ALLOCATABLE :: weight_v(:) ! Localization weights for vertical (for locweight_v>0) REAL, ALLOCATABLE :: resid_obs(:,:) ! Array for a single row of resid_l REAL, ALLOCATABLE :: Rinvresid_l(:) ! R^-1 times residual REAL :: lhood_one ! Likelihood for this observation doassim: IF (thisobs%doassim == 1) THEN ! ********************** ! *** INITIALIZATION *** ! ********************** ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_likelihood_l -- START localization and likelihood, member', obs_member WRITE (*,*) '++ OMI-debug likelihood_l: ', debug, ' thisobs_l%locweight', thisobs_l%locweight END IF ! Screen output IF (verbose == 1) THEN IF (thisobs%obs_err_type==0) THEN WRITE (*, '(a, 5x, a)') & 'PDAFomi', '--- Assume Gaussian observation errors' ELSE WRITE (*, '(a, 5x, a)') & 'PDAFomi', '--- Assume double-exponential observation errors' END IF WRITE (*, '(a, 5x, a, 1x)') & 'PDAFomi', '--- Domain localization' WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Localization cut-off radius', thisobs_l%cradius WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Support radius', thisobs_l%sradius ENDIF ! *********************************************** ! *** Before computing the likelihood, apply *** ! *** the localization weight and scale by *** ! *** observation error variance *** ! *** -1 *** ! *** Rinvresid = Weight R resid *** ! *** *** ! *** Apply a weight matrix with correlations *** ! *** of compact support to residual or the *** ! *** observation error variance. *** ! *********************************************** ! *** Initialize weight array ALLOCATE(weight(thisobs_l%dim_obs_l)) ALLOCATE(resid_obs(thisobs_l%dim_obs_l,1)) resid_obs(:,1) = resid_l(:) CALL PDAFomi_weights_l(verbose, thisobs_l%dim_obs_l, 1, thisobs_l%locweight, & thisobs_l%cradius_l, thisobs_l%sradius_l, & resid_obs, thisobs_l%ivar_obs_l, thisobs_l%distance_l, weight) ! *** For factorized 2+1D localization use product of horizontal and vertical weights IF (thisobs_l%locweight_v>0) then IF (verbose == 1) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- initialize also weight function for vertical direction' END IF ALLOCATE(weight_v(thisobs_l%dim_obs_l)) CALL PDAFomi_weights_l_sgnl(verbose, thisobs_l%dim_obs_l, 1, thisobs_l%locweight_v, & thisobs_l%cradius(3), thisobs_l%sradius(3), & resid_obs, thisobs_l%ivar_obs_l, thisobs_l%dist_l_v, weight_v) DO i = 1, thisobs_l%dim_obs_l weight(i) = weight(i) * weight_v(i) END DO DEALLOCATE(weight_v) END IF DEALLOCATE(resid_obs) ! *** Handling of special weighting types *** lw2: IF (thisobs_l%locweight == 26) THEN ! Use square-root of 5th-order polynomial on A DO i = 1, thisobs_l%dim_obs_l ! Check if weight >0 (Could be <0 due to numerical precision) IF (weight(i) > 0.0) THEN weight(i) = SQRT(weight(i)) ELSE weight(i) = 0.0 END IF END DO END IF lw2 ! *** Apply weight ALLOCATE(Rinvresid_l(thisobs_l%dim_obs_l)) DO i = 1, thisobs_l%dim_obs_l Rinvresid_l(i) = thisobs_l%ivar_obs_l(i) * weight(i) * resid_l(thisobs_l%off_obs_l+i) END DO ! ******************************** ! *** Compute local 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_l>0.0) lhood_l = - LOG(lhood_l) lhood_one = 0.0 DO i = 1, thisobs_l%dim_obs_l lhood_one = lhood_one + 0.5*resid_l(thisobs_l%off_obs_l+i)*Rinvresid_l(i) END DO lhood_l = EXP(-(lhood_l + lhood_one)) ELSE ! Double-exponential errors ! Calculate exp(-SUM(ABS(resid))) ! Transform pack to log likelihood to increment its values IF (lhood_l>0.0) lhood_l = - LOG(lhood_l) lhood_one = 0.0 DO i = 1, thisobs_l%dim_obs_l lhood_one = lhood_one + ABS(Rinvresid_l(i)) END DO lhood_l = EXP(-(lhood_l + lhood_one)) END IF ! *** Clean up *** DEALLOCATE(weight, Rinvresid_l) ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug likelihood_l: ', debug, ' accumulated likelihood', lhood_l WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_likelihood_l -- END' END IF END IF doassim END SUBROUTINE PDAFomi_likelihood_l !------------------------------------------------------------------------------- !> Compute local likelihood for an ensemble member using hybrid weight !! !! The routine is called during the analysis step !! of the localized NETF. !! 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 addition, a localizing weighting of the !! inverse of R by expotential decrease or a 5-th order !! polynomial of compact support can be applied. This is !! defined by the variables 'locweight', 'cradius', !! 'cradius2' and 'sradius' in the main program. !! A tempering is appply by using the hybrid weight 'gamma'. !! !! In general this routine is similar to the routine !! prodRinvA_l 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. !! !! This implementation assumes a diagonal observation !! error covariance matrix, and supports varying !! observation error variances. !! !! The routine can be applied with either all observations !! of different types at once, or separately for each !! observation type. !! !! __Revision history:__ !! * 2022-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_likelihood_hyb_l(thisobs_l, thisobs, resid_l, gamma, lhood_l, verbose) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(inout) :: resid_l(:) !< Input vector of residuum REAL, INTENT(inout) :: lhood_l !< Output vector - log likelihood REAL, INTENT(in) :: gamma !< Hybrid weight INTEGER, INTENT(in) :: verbose !< Verbosity flag ! *** local variables *** INTEGER :: i ! Index of observation component REAL, ALLOCATABLE :: weight(:) ! Localization weights REAL, ALLOCATABLE :: weight_v(:) ! Localization weights for vertical (for locweight_v>0) REAL, ALLOCATABLE :: resid_obs(:,:) ! Array for a single row of resid_l REAL, ALLOCATABLE :: Rinvresid_l(:) ! R^-1 times residual REAL :: lhood_one ! Likelihood for this observation doassim: IF (thisobs%doassim == 1) THEN ! ********************** ! *** INITIALIZATION *** ! ********************** ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, & 'PDAFomi_likelihood_hyb_l -- START localization and likelihood, member', obs_member WRITE (*,*) '++ OMI-debug likelihood_hyb_l: ', debug, ' thisobs_l%locweight', thisobs_l%locweight END IF ! Screen output IF (verbose == 1) THEN IF (thisobs%obs_err_type==0) THEN WRITE (*, '(a, 5x, a)') & 'PDAFomi', '--- Assume Gaussian observation errors' ELSE WRITE (*, '(a, 5x, a)') & 'PDAFomi', '--- Assume double-exponential observation errors' END IF WRITE (*, '(a, 5x, a, f12.5)') & 'PDAFomi', '--- Apply tempering with 1.0-gamma= ', 1.0 - gamma WRITE (*, '(a, 5x, a, 1x)') & 'PDAFomi', '--- Domain localization' WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Localization cut-off radius', thisobs_l%cradius WRITE (*, '(a, 8x, a, 1x, es11.3)') & 'PDAFomi', '--- Support radius', thisobs_l%sradius ENDIF ! *********************************************** ! *** Before computing the likelihood, apply *** ! *** the localization weight and scale by *** ! *** observation error variance *** ! *** -1 *** ! *** Rinvresid = Weight R resid *** ! *** *** ! *** Apply a weight matrix with correlations *** ! *** of compact support to residual or the *** ! *** observation error variance. *** ! *********************************************** ! *** Initialize weight array ALLOCATE(weight(thisobs_l%dim_obs_l)) ALLOCATE(resid_obs(thisobs_l%dim_obs_l,1)) resid_obs(:,1) = resid_l(:) CALL PDAFomi_weights_l(verbose, thisobs_l%dim_obs_l, 1, thisobs_l%locweight, & thisobs_l%cradius_l, thisobs_l%sradius_l, & resid_obs, thisobs_l%ivar_obs_l, thisobs_l%distance_l, weight) ! *** For factorized 2+1D localization use product of horizontal and vertical weights IF (thisobs_l%locweight_v>0) then IF (verbose == 1) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- initialize also weight function for vertical direction' END IF ALLOCATE(weight_v(thisobs_l%dim_obs_l)) CALL PDAFomi_weights_l_sgnl(verbose, thisobs_l%dim_obs_l, 1, thisobs_l%locweight_v, & thisobs_l%cradius(3), thisobs_l%sradius(3), & resid_obs, thisobs_l%ivar_obs_l, thisobs_l%dist_l_v, weight_v) DO i = 1, thisobs_l%dim_obs_l weight(i) = weight(i) * weight_v(i) END DO DEALLOCATE(weight_v) END IF DEALLOCATE(resid_obs) ! *** Handling of special weighting types *** lw2: IF (thisobs_l%locweight == 26) THEN ! Use square-root of 5th-order polynomial on A DO i = 1, thisobs_l%dim_obs_l ! Check if weight >0 (Could be <0 due to numerical precision) IF (weight(i) > 0.0) THEN weight(i) = SQRT(weight(i)) ELSE weight(i) = 0.0 END IF END DO END IF lw2 ! *** Apply weight ALLOCATE(Rinvresid_l(thisobs_l%dim_obs_l)) DO i = 1, thisobs_l%dim_obs_l Rinvresid_l(i) = (1.0-gamma) * thisobs_l%ivar_obs_l(i) * weight(i) * resid_l(i) END DO ! ******************************** ! *** Compute local likelihood *** ! ******************************** IF (thisobs%obs_err_type == 0) THEN ! Gaussian errors ! Calculate exp(-0.5*resid^T*R^-1*resid) ! Transform pack to log likelihood to increment its values IF (lhood_l>0.0) lhood_l = - LOG(lhood_l) CALL dgemv('t', thisobs_l%dim_obs_l, 1, 0.5, resid_l, & thisobs_l%dim_obs_l, Rinvresid_l, 1, 0.0, lhood_one, 1) lhood_l = EXP(-(lhood_l + lhood_one)) ELSE ! Double-exponential errors ! Calculate exp(-SUM(ABS(resid))) ! Transform pack to log likelihood to increment its values IF (lhood_l>0.0) lhood_l = - LOG(lhood_l) lhood_one = 0.0 DO i = 1, thisobs_l%dim_obs_l lhood_one = lhood_one + ABS(Rinvresid_l(i)) END DO lhood_l = EXP(-(lhood_l + lhood_one)) END IF ! *** Clean up *** DEALLOCATE(weight, Rinvresid_l) ! Screen output IF (debug>0) THEN WRITE (*,*) '++ OMI-debug likelihood_hyb_l: ', debug, ' accumulated likelihood', lhood_l WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_likelihood_hyb)l -- END' END IF END IF doassim END SUBROUTINE PDAFomi_likelihood_hyb_l !------------------------------------------------------------------------------- !> Apply covariance localization !! !! This routine applies a localization matrix B !! to the matrices HP and HPH^T of the localized EnKF. !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_localize_covar_iso(thisobs, dim, locweight, cradius, sradius, & coords, HP, HPH) USE PDAFomi_obs_f, & ONLY: obsdims USE PDAF_mod_filtermpi, & ONLY: npes IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: dim !< State dimension INTEGER, INTENT(in) :: locweight !< Localization weight type REAL, INTENT(in) :: cradius !< localization radius REAL, INTENT(in) :: sradius !< support radius for weight functions REAL, INTENT(in) :: coords(:,:) !< Coordinates of state vector elements REAL, INTENT(inout) :: HP(:, :) !< Matrix HP, dimension (nobs, dim) REAL, INTENT(inout) :: HPH(:, :) !< Matrix HPH, dimension (nobs, nobs) ! *** local variables *** INTEGER :: i, j, pe, cnt ! counters INTEGER :: ncoord ! Number of coordinates REAL :: distance ! Distance between points in the domain REAL :: weight ! Localization weight REAL, ALLOCATABLE :: weights(:) ! Localization weights array REAL :: tmp(1,1)= 1.0 ! Temporary, but unused array INTEGER :: wtype ! Type of weight function INTEGER :: rtype ! Type of weight regulation REAL, ALLOCATABLE :: co(:), oc(:) ! Coordinates of single point 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 INTEGER, ALLOCATABLE :: obs_map(:) ! Mapping indiced for observations doassim: IF (thisobs%doassim == 1) THEN ! ********************** ! *** INITIALIZATION *** ! ********************** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_localize_covar -- START' WRITE (*,*) '++ OMI-debug localize_covar:', debug, 'thisobs%off_obs_g', thisobs%off_obs_g WRITE (*,*) '++ OMI-debug localize_covar:', debug, 'thisobs%dim_obs_g', thisobs%dim_obs_g END IF ! Screen output IF (screen > 0 .AND. mype==0) THEN WRITE (*,'(a, 8x, a, 1x, i3)') & 'PDAFomi', '--- Apply covariance localization, obs. type ID', thisobs%obsid WRITE (*, '(a, 12x, a, 1x, f12.2)') & 'PDAFomi', '--- Local influence radius', cradius IF (locweight == 0) THEN WRITE (*, '(a, 12x, a)') & 'PDAFomi', '--- Use uniform weight' ELSE IF (locweight == 1) THEN WRITE (*, '(a, 12x, a)') & 'PDAFomi', '--- Use exponential distance-dependent weight' ELSE IF (locweight == 2) THEN WRITE (*, '(a, 12x, a)') & 'PDAFomi', '--- Use distance-dependent weight by 5th-order polynomial' END IF ENDIF ! Set ncoord locally for compact code ncoord = thisobs%ncoord ! *** Initialize mapping indices ALLOCATE(id_start(npes), id_end(npes)) ! thisobs%ocoord_f is global for all observations of one type ! while HP and HPH are ordered obstype-first. Thus the observations ! of all types of one sub-domain are combined. This mapping ! ensures that the correct indices are used in HP and HPH. 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 ALLOCATE(obs_map(thisobs%dim_obs_g)) cnt = 1 DO pe = 1, npes DO i = id_start(pe), id_end(pe) obs_map(cnt) = i cnt = cnt + 1 END DO END DO ! ************************** ! *** Apply localization *** ! ************************** ! Set parameters for weight calculation IF (locweight == 0) THEN ! Uniform (unit) weighting wtype = 0 rtype = 0 ELSE IF (locweight == 1) THEN ! Exponential weighting wtype = 1 rtype = 0 ELSE IF (locweight == 2) THEN ! 5th-order polynomial (Gaspari&Cohn, 1999) wtype = 2 rtype = 0 END IF ALLOCATE(oc(ncoord)) ALLOCATE(co(ncoord)) ! *** Localize HP *** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug localize_covar:', debug, ' localize matrix HP' END IF ALLOCATE(weights(thisobs%dim_obs_g)) DO i = 1, dim ! Initialize coordinate co(1:ncoord) = coords(1:thisobs%ncoord, i) DO j = 1, thisobs%dim_obs_g ! Initialize coordinate oc(1:ncoord) = thisobs%ocoord_f(1:thisobs%ncoord, j) ! Compute distance CALL PDAFomi_comp_dist2(thisobs, co, oc, distance, (i*j)-1) distance = SQRT(distance) ! Compute weight CALL PDAF_local_weight(wtype, rtype, cradius, sradius, distance, & 1, 1, tmp, 1.0, weights(j), 0) END DO IF (debug==i) THEN WRITE (*,*) '++ OMI-debug localize_covar: ', debug, 'weights for row in HP', weights END IF DO j = 1, thisobs%dim_obs_g ! Apply localization HP(obs_map(j), i) = weights(j) * HP(obs_map(j), i) END DO END DO DEALLOCATE(weights) ! *** Localize HPH^T *** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug localize_covar:', debug, ' localize matrix HPH^T' END IF DO i = 1, thisobs%dim_obs_g ! Initialize coordinate co(1:ncoord) = thisobs%ocoord_f(1:thisobs%ncoord, i) DO j = 1, thisobs%dim_obs_g ! Initialize coordinate oc(1:ncoord) = thisobs%ocoord_f(1:thisobs%ncoord, j) ! Compute distance CALL PDAFomi_comp_dist2(thisobs, co, oc, distance, (i*j)-1) distance = SQRT(distance) ! Compute weight CALL PDAF_local_weight(wtype, rtype, cradius, sradius, distance, & 1, 1, tmp, 1.0, weight, 0) ! Apply localization HPH(obs_map(j), obs_map(i)) = weight * HPH(obs_map(j), obs_map(i)) END DO END DO ! clean up DEALLOCATE(co, oc) DEALLOCATE(id_start, id_end, obs_map) IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_localize_covar -- END' END IF doassim END SUBROUTINE PDAFomi_localize_covar_iso !------------------------------------------------------------------------------- !> Apply covariance localization: 2+1D factorized with vertical localization weight !! !! This routine is a variant of PDAFomi_localize_covar_noniso with !! support for a vector of localization weights. This is used !! to specify different localization functions for the vertical and !! horizontal directions. The routine only stores the value of !! locweights(2) for the vertical and calls PDAFomi_init_dim_obs_l_noniso. !! !! __Revision history:__ !! * 2024-04 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_localize_covar_noniso_locweights(thisobs, dim, locweights, cradius, sradius, & coords, HP, HPH) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: dim !< State dimension INTEGER, INTENT(in) :: locweights(:) !< Types of localization function REAL, INTENT(in) :: cradius(:) !< Vector of localization cut-off radii REAL, INTENT(in) :: sradius(:) !< Vector of support radii of localization function REAL, INTENT(in) :: coords(:,:) !< Coordinates of state vector elements REAL, INTENT(inout) :: HP(:, :) !< Matrix HP, dimension (nobs, dim) REAL, INTENT(inout) :: HPH(:, :) !< Matrix HPH, dimension (nobs, nobs) ! *** Store vertical locweight and call standard routine IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_localize_covar_noniso_locweights -- START' WRITE (*,*) '++ OMI-debug PDAFomi_localize_covar_noniso_locweights:', debug, ' locweights', locweights END IF ! Check consistency of dimensions IF (SIZE(locweights) /= 2) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: Input for locweights in horizontal and vertical needs size 2' error = 15 END IF IF (thisobs%ncoord /= 3) THEN WRITE (*,*) '+++++ WARNING PDAF-OMI: separate locweight for vertical is only utilized if thisobs%ncoord=3' END IF IF (thisobs%ncoord == 3) THEN ! locweight for the vertical is treated separately thisobs%locweight_v = locweights(2) END IF CALL PDAFomi_localize_covar_noniso(thisobs, dim, locweights(1), cradius, sradius, & coords, HP, HPH) END SUBROUTINE PDAFomi_localize_covar_noniso_locweights !------------------------------------------------------------------------------- !> Apply covariance localization !! !! This routine applies a localization matrix B !! to the matrices HP and HPH^T of the localized EnKF. !! This variant is for non-iceotropic localization !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_localize_covar_noniso(thisobs, dim, locweight, cradius, sradius, & coords, HP, HPH) USE PDAFomi_obs_f, & ONLY: obsdims USE PDAF_mod_filtermpi, & ONLY: npes IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation INTEGER, INTENT(in) :: dim !< State dimension INTEGER, INTENT(in) :: locweight !< Localization weight type REAL, INTENT(in) :: cradius(:) !< Vector of localization cut-off radii REAL, INTENT(in) :: sradius(:) !< Vector of support radii of localization function REAL, INTENT(in) :: coords(:,:) !< Coordinates of state vector elements REAL, INTENT(inout) :: HP(:, :) !< Matrix HP, dimension (nobs, dim) REAL, INTENT(inout) :: HPH(:, :) !< Matrix HPH, dimension (nobs, nobs) ! *** local variables *** INTEGER :: i, j, pe, cnt ! counters INTEGER :: ncoord ! Number of coordinates REAL :: distance ! Distance between points in the domain REAL :: weight ! Localization weight REAL, ALLOCATABLE :: weights(:) ! Localization weights array REAL :: weight_v ! Weights in vertical direction REAL :: tmp(1,1)= 1.0 ! Temporary, but unused array INTEGER :: wtype ! Type of weight function INTEGER :: rtype ! Type of weight regulation REAL :: srad, crad ! localization cut-off radius REAL, ALLOCATABLE :: co(:), oc(:) ! Coordinates of single point 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 INTEGER, ALLOCATABLE :: obs_map(:) ! Mapping indiced for observations LOGICAL :: checkdist ! Flag whether distance nis not larger than cut-off radius REAL :: dists(thisobs%ncoord) ! Distance vector between analysis point and observation TYPE(obs_l) :: thisobs_l ! local observation doassim: IF (thisobs%doassim == 1) THEN ! ********************** ! *** INITIALIZATION *** ! ********************** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_localize_covar_noniso -- START' WRITE (*,*) '++ OMI-debug localize_covar_noniso:', debug, 'thisobs%off_obs_g', thisobs%off_obs_g WRITE (*,*) '++ OMI-debug localize_covar_noniso:', debug, 'thisobs%dim_obs_g', thisobs%dim_obs_g END IF ! Check consistency of dimensions IF (SIZE(cradius) /= thisobs%ncoord) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: non-isotropic localization: Size of CRADIUS /= thisobs%ncoord' error = 12 END IF IF (SIZE(sradius) /= thisobs%ncoord) THEN WRITE (*,*) '+++++ ERROR PDAF-OMI: non-isotropic localization: Size of SRADIUS /= thisobs%ncoord' error = 13 END IF ! Screen output IF (screen > 0 .AND. mype==0) THEN WRITE (*,'(a, 8x, a)') & 'PDAFomi', '--- Apply non-isotropic covariance localization' IF (thisobs%ncoord==1) THEN WRITE (*, '(a, 12x, a, 1x, es11.3)') & 'PDAFomi', '--- Local influence radius', cradius ELSEIF (thisobs%ncoord==2) THEN WRITE (*, '(a, 12x, a, 1x, 2es11.3)') & 'PDAFomi', '--- Local influence radii', cradius ELSE WRITE (*, '(a, 12x, a, 1x, 3f11.3)') & 'PDAFomi', '--- Local influence radii', cradius END IF IF (locweight == 0) THEN WRITE (*, '(a, 12x, a)') & 'PDAFomi', '--- Use uniform weight' ELSE IF (locweight == 1) THEN WRITE (*, '(a, 12x, a)') & 'PDAFomi', '--- Use exponential distance-dependent weight' ELSE IF (locweight == 2) THEN WRITE (*, '(a, 12x, a)') & 'PDAFomi', '--- Use distance-dependent weight by 5th-order polynomial' END IF ENDIF ! Set ncoord locally for compact code ncoord = thisobs%ncoord ! *** Initialize mapping indices ALLOCATE(id_start(npes), id_end(npes)) ! thisobs%ocoord_f is global for all observations of one type ! while HP and HPH are ordered obstype-first. Thus the observations ! of all types of one sub-domain are combined. This mapping ! ensures that the correct indices are used in HP and HPH. 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 ALLOCATE(obs_map(thisobs%dim_obs_g)) cnt = 1 DO pe = 1, npes DO i = id_start(pe), id_end(pe) obs_map(cnt) = i cnt = cnt + 1 END DO END DO ! Allocate vectors for localization radii and store their values IF (ALLOCATED(thisobs_l%cradius)) DEALLOCATE(thisobs_l%cradius) ALLOCATE(thisobs_l%cradius(thisobs%ncoord)) IF (ALLOCATED(thisobs_l%sradius)) DEALLOCATE(thisobs_l%sradius) ALLOCATE(thisobs_l%sradius(thisobs%ncoord)) ! Set number of localization radii thisobs_l%nradii = thisobs%ncoord thisobs_l%cradius(:) = cradius(:) thisobs_l%sradius(:) = sradius(:) ! ************************** ! *** Apply localization *** ! ************************** ! Set parameters for weight calculation IF (locweight == 0) THEN ! Uniform (unit) weighting wtype = 0 rtype = 0 ELSE IF (locweight == 1) THEN ! Exponential weighting wtype = 1 rtype = 0 ELSE IF (locweight == 2) THEN ! 5th-order polynomial (Gaspari&Cohn, 1999) wtype = 2 rtype = 0 END IF ALLOCATE(oc(ncoord)) ALLOCATE(co(ncoord)) ! *** Localize HP *** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug localize_covar:', debug, ' localize matrix HP' END IF ALLOCATE(weights(thisobs%dim_obs_g)) cnt=0 DO i = 1, dim ! Initialize coordinate co(1:ncoord) = coords(1:thisobs%ncoord, i) DO j = 1, thisobs%dim_obs_g ! Initialize coordinate oc(1:ncoord) = thisobs%ocoord_f(1:thisobs%ncoord, j) ! Compute distance CALL PDAFomi_check_dist2_noniso(thisobs, thisobs_l, co, oc, distance, & dists, crad, srad, checkdist, (i*j)-1, cnt) distance = SQRT(distance) ! Compute weight CALL PDAF_local_weight(wtype, rtype, crad, srad, distance, & 1, 1, tmp, 1.0, weights(j), 0) ! Compute separate weight for vertical direction (for factorized 2+1D localization) IF (thisobs%locweight_v>0 .AND. thisobs%ncoord==3) THEN CALL PDAF_local_weight(wtype, rtype, cradius(3), sradius(3), dists(3), & 1, 1, tmp, 1.0, weight_v, 0) weights(j) = weights(j) * weight_v END IF END DO IF (debug==i) THEN WRITE (*,*) '++ OMI-debug localize_covar: ', debug, 'weights for row in HP', weights END IF DO j = 1, thisobs%dim_obs_g ! Apply localization HP(obs_map(j), i) = weights(j) * HP(obs_map(j), i) END DO END DO DEALLOCATE(weights) ! *** Localize HPH^T *** IF (debug>0) THEN WRITE (*,*) '++ OMI-debug localize_covar:', debug, ' localize matrix HPH^T' END IF cnt=0 DO i = 1, thisobs%dim_obs_g ! Initialize coordinate co(1:ncoord) = thisobs%ocoord_f(1:thisobs%ncoord, i) DO j = 1, thisobs%dim_obs_g ! Initialize coordinate oc(1:ncoord) = thisobs%ocoord_f(1:thisobs%ncoord, j) ! Compute distance CALL PDAFomi_check_dist2_noniso(thisobs, thisobs_l, co, oc, distance, & dists, crad, srad, checkdist, (i*j)-1, cnt) distance = SQRT(distance) ! Compute weight CALL PDAF_local_weight(wtype, rtype, crad, srad, distance, & 1, 1, tmp, 1.0, weight, 0) ! Compute separate weight for vertical direction (for factorized 2+1D localization) IF (thisobs%locweight_v>0 .AND. thisobs%ncoord==3) THEN CALL PDAF_local_weight(wtype, rtype, cradius(3), sradius(3), dists(3), & 1, 1, tmp, 1.0, weight_v, 0) weights(j) = weights(j) * weight_v END IF ! Apply localization HPH(obs_map(j), obs_map(i)) = weight * HPH(obs_map(j), obs_map(i)) END DO END DO ! clean up DEALLOCATE(co, oc) DEALLOCATE(id_start, id_end, obs_map) IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_localize_covar -- END' END IF doassim END SUBROUTINE PDAFomi_localize_covar_noniso !------------------------------------------------------------------------------- !> Initialize local observation vector !! !! This routine has to initialize the part of the !! overall local observation vector corresponding !! to the current observation type. The offset of !! the current observation type in the local obs. !! vector is given by OFFSET_OBS_l_ALL. !! This routine is both used directly to initialize !! the local part of a vector in observation space !! (PDAFomi_g2l_obs) and it's called from !! PDAFomi_init_obs_l to initialize the local part !! of the vector of observations and the corresponding !! vector of variances. !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_g2l_obs_internal(thisobs_l, obs_f_one, offset_obs_l_all, obs_l_all) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation REAL, INTENT(in) :: obs_f_one(:) !< Full obs. vector of current obs. type (nobs_f_one) REAL, INTENT(inout) :: obs_l_all(:) !< Local observation vector for all variables (nobs_l_all) INTEGER, INTENT(in) :: offset_obs_l_all !< Offset of current observation in obs_l_all and ivar_l_all ! *** Local variables *** INTEGER :: i ! Counter ! ******************************************* ! *** Initialize local observation vector *** ! ******************************************* ! Local observations DO i = 1, thisobs_l%dim_obs_l obs_l_all(i+offset_obs_l_all) = obs_f_one(thisobs_l%id_obs_l(i)) ENDDO ! Print debug information IF (debug>0) THEN IF (thisobs_l%dim_obs_l>0) THEN WRITE (*,*) '++ OMI-debug g2l_obs: ', debug, ' thisobs_l%id_obs_l', thisobs_l%id_obs_l WRITE (*,*) '++ OMI-debug g2l_obs: ', debug, ' obs_l', & obs_l_all(1+offset_obs_l_all:offset_obs_l_all+thisobs_l%dim_obs_l) ELSE WRITE (*,*) '++ OMI-debug g2l_obs: ', debug, ' no local observations present' END IF END IF END SUBROUTINE PDAFomi_g2l_obs_internal !------------------------------------------------------------------------------- !> Compute square distance between two locations !! !! This routine computes the distance between two locations. !! The computation can be for Cartesian grids with and without !! periodicity and for geographic coordinates. For Cartesian !! grids, the coordinats can be in any unit, while geographic !! coordinates must be provided in radians and the resulting !! distance will be in meters. !! !! Choices for distance computation - disttype: !! 0: Cartesian distance in ncoord dimensions !! 1: Cartesian distance in ncoord dimensions with periodicity !! (Needs specification of domsize(ncoord)) !! 2: Aproximate geographic distance with horizontal coordinates in radians (-pi/+pi) !! 3: Geographic distance computation using haversine formula !! !! __Revision history:__ !! * 2019-06 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_comp_dist2(thisobs, coordsA, coordsB, distance2, verbose) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation REAL, INTENT(in) :: coordsA(:) !< Coordinates of current analysis domain (ncoord) REAL, INTENT(in) :: coordsB(:) !< Coordinates of observation (ncoord) REAL, INTENT(out) :: distance2 !< Squared distance INTEGER, INTENT(in) :: verbose !< Control screen output ! *** Local variables *** INTEGER :: k ! Counters REAL :: dists(thisobs%ncoord) ! Distance vector between analysis point and observation REAL :: slon, slat ! sine of distance in longitude or latitude INTEGER :: domsize ! Flag whether domainsize is set ! ************************ ! *** Compute distance *** ! ************************ IF (.NOT.ALLOCATED(thisobs%domainsize)) THEN domsize = 0 ELSE domsize = 1 END IF norm: IF ((thisobs%disttype==0 .OR. thisobs%disttype==10) .OR. & ((thisobs%disttype==1 .OR. thisobs%disttype==11) .AND. domsize==0)) THEN ! *** Compute Cartesian distance *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug comp_dist2: ', debug, ' compute Cartesian distance' END IF ! Distance per direction DO k = 1, thisobs%ncoord dists(k) = ABS(coordsA(k) - coordsB(k)) END DO ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSEIF ((thisobs%disttype==1 .OR. thisobs%disttype==11) .AND. domsize==1) THEN norm ! *** Compute periodic Cartesian distance *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug comp_dist2: ', debug, ' compute periodic Cartesian distance' END IF ! Distance per direction DO k = 1, thisobs%ncoord IF (thisobs%domainsize(k)<=0.0) THEN dists(k) = ABS(coordsA(k) - coordsB(k)) ELSE dists(k) = MIN(ABS(coordsA(k) - coordsB(k)), & ABS(ABS(coordsA(k) - coordsB(k))-thisobs%domainsize(k))) END IF END DO ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF ELSEIF (thisobs%disttype==2 .OR. thisobs%disttype==12) THEN norm ! *** Compute distance from geographic coordinates *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug comp_dist2: ', debug, ' compute geographic distance' END IF ! approximate distances in longitude and latitude dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (thisobs%ncoord>2) dists(3) = ABS(coordsA(3) - coordsB(3)) ! full squared distance in meters distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSEIF (thisobs%disttype==3 .OR. thisobs%disttype==13) THEN norm ! *** Compute distance from geographic coordinates with haversine formula *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug comp_dist2: ', debug, & ' compute geographic distance using haversine function' END IF slon = SIN((coordsA(1) - coordsB(1))/2) slat = SIN((coordsA(2) - coordsB(2))/2) dists(2) = SQRT(slat*slat + COS(coordsA(2))*COS(coordsB(2))*slon*slon) IF (dists(2)<=1.0) THEN dists(2) = 2.0 * r_earth* ASIN(dists(2)) ELSE dists(2) = r_earth* pi END IF IF (thisobs%ncoord>2) dists(3) = ABS(coordsA(3) - coordsB(3)) ! full squared distance in meters distance2 = 0.0 DO k = 2, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF norm END SUBROUTINE PDAFomi_comp_dist2 !------------------------------------------------------------------------------- !> Check distance in case of isotropic localization !! !! This routine computes the distance between two locations. !! The computation can be for Cartesian grids with and without !! periodicity and for geographic coordinates. For Cartesian !! grids, the coordinates can be in any unit, while geographic !! coordinates must be provided in radians and the resulting !! distance will be in meters. Finally, the routine checks !! whether the distance is not larger than the cut-off radius. !! !! Choices for distance computation - disttype: !! 0: Cartesian distance in ncoord dimensions !! 1: Cartesian distance in ncoord dimensions with periodicity !! (Needs specification of domsize(ncoord)) !! 2: Aproximate geographic distance with horizontal coordinates in radians (-pi/+pi) !! 3: Geographic distance computation using haversine formula !! 10-13: Variants of distance types 0-3, but particularly for 3 dimensions in which !! a 2+1 dimensional localization is applied (distance weighting only in the horizontal) !! !! __Revision history:__ !! * 2024-04 - Lars Nerger - Initial code based on PDAFomi_comp_dist2 !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_check_dist2(thisobs, thisobs_l, coordsA, coordsB, distance2, & checkdist, verbose, cnt_obs) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation TYPE(obs_l), INTENT(in) :: thisobs_l !< Data type with local observation REAL, INTENT(in) :: coordsA(:) !< Coordinates of current analysis domain (ncoord) REAL, INTENT(in) :: coordsB(:) !< Coordinates of observation (ncoord) REAL, INTENT(out) :: distance2 !< Squared distance LOGICAL, INTENT(out) :: checkdist !< Flag whether distance is within cut-off radius INTEGER, INTENT(in) :: verbose !< Control screen output INTEGER, INTENT(inout) :: cnt_obs !< Count number of local observations ! *** Local variables *** INTEGER :: k ! Counters REAL :: dists(thisobs%ncoord) ! Distance vector between analysis point and observation REAL :: slon, slat ! sine of distance in longitude or latitude REAL :: cradius2 ! squared localization cut-off radius INTEGER :: domsize ! Flag whether domainsize is set LOGICAL :: distflag ! Flag whether distance in a coordinate direction is within cradius ! ********************** ! *** Initialization *** ! ********************** ! Initialize distance flag checkdist = .FALSE. distflag = .TRUE. ! ************************ ! *** Compute distance *** ! ************************ IF (.NOT.ALLOCATED(thisobs%domainsize)) THEN domsize = 0 ELSE domsize = 1 END IF norm: IF ((thisobs%disttype==0 .OR. thisobs%disttype==10) .OR. & ((thisobs%disttype==1 .OR. thisobs%disttype==11) .AND. domsize==0)) THEN ! *** Compute Cartesian distance *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2: ', debug, ' compute Cartesian distance' END IF IF (thisobs%ncoord>=3) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) IF (dists(3)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(2) = ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(1) = ABS(coordsA(1) - coordsB(1)) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSEIF (thisobs%ncoord==2) THEN dists(2) = ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(1) = ABS(coordsA(1) - coordsB(1)) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF (thisobs%ncoord==1) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF ((thisobs%disttype==1 .OR. thisobs%disttype==11) .AND. domsize==1) THEN norm ! *** Compute periodic Cartesian distance *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2: ', debug, ' compute periodic Cartesian distance' END IF IF (thisobs%ncoord>=3) THEN IF (thisobs%domainsize(3)<=0.0) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) ELSE dists(3) = MIN(ABS(coordsA(3) - coordsB(3)), & ABS(ABS(coordsA(3) - coordsB(3))-thisobs%domainsize(3))) END IF IF (dists(3)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE IF (thisobs%domainsize(2)<=0.0) THEN dists(2) = ABS(coordsA(2) - coordsB(2)) ELSE dists(2) = MIN(ABS(coordsA(2) - coordsB(2)), & ABS(ABS(coordsA(2) - coordsB(2))-thisobs%domainsize(2))) END IF IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE IF (thisobs%domainsize(1)<=0.0) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) ELSE dists(1) = MIN(ABS(coordsA(1) - coordsB(1)), & ABS(ABS(coordsA(1) - coordsB(1))-thisobs%domainsize(1))) END IF IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSEIF (thisobs%ncoord==2) THEN IF (thisobs%domainsize(2)<=0.0) THEN dists(2) = ABS(coordsA(2) - coordsB(2)) ELSE dists(2) = MIN(ABS(coordsA(2) - coordsB(2)), & ABS(ABS(coordsA(2) - coordsB(2))-thisobs%domainsize(2))) END IF IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE IF (thisobs%domainsize(1)<=0.0) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) ELSE dists(1) = MIN(ABS(coordsA(1) - coordsB(1)), & ABS(ABS(coordsA(1) - coordsB(1))-thisobs%domainsize(1))) END IF IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF (thisobs%ncoord==1) THEN IF (thisobs%domainsize(1)<=0.0) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) ELSE dists(1) = MIN(ABS(coordsA(1) - coordsB(1)), & ABS(ABS(coordsA(1) - coordsB(1))-thisobs%domainsize(1))) END IF IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF (thisobs%disttype==2 .OR. thisobs%disttype==12) THEN norm ! *** Compute distance from geographic coordinates *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2: ', debug, ' compute geographic distance' END IF IF (thisobs%ncoord==3) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) IF (dists(3)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF ELSEIF (thisobs%disttype==3 .OR. thisobs%disttype==13) THEN norm ! *** Compute distance from geographic coordinates with haversine formula *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2: ', debug, & ' compute geographic distance using haversine function' END IF IF (thisobs%ncoord==3) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) IF (dists(3)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! Haversine formula slon = SIN((coordsA(1) - coordsB(1))/2) slat = SIN((coordsA(2) - coordsB(2))/2) dists(2) = SQRT(slat*slat + COS(coordsA(2))*COS(coordsB(2))*slon*slon) IF (dists(2)<=1.0) THEN dists(2) = 2.0 * r_earth* ASIN(dists(2)) ELSE dists(2) = r_earth* pi END IF IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 2, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 2, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! Haversine formula slon = SIN((coordsA(1) - coordsB(1))/2) slat = SIN((coordsA(2) - coordsB(2))/2) dists(2) = SQRT(slat*slat + COS(coordsA(2))*COS(coordsB(2))*slon*slon) IF (dists(2)<=1.0) THEN dists(2) = 2.0 * r_earth* ASIN(dists(2)) ELSE dists(2) = r_earth* pi END IF IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF norm IF (distflag) THEN cradius2 = thisobs_l%cradius(1)*thisobs_l%cradius(1) IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. ! Increment counter cnt_obs = cnt_obs + 1 END IF END IF END SUBROUTINE PDAFomi_check_dist2 !------------------------------------------------------------------------------- !> Check distance in case of nonisotropic localization !! !! This routine computes the distance between the observation and a !! model grid point and the cut-off radius of an ellipse (in 2D) !! or ellipsoid (in 3D) in the direction of the distance. Finally, !! the routine checks whether the distance is not larger than the !! cut-off radius. !! !! Choices for distance computation - disttype: !! 0: Cartesian distance in ncoord dimensions !! 1: Cartesian distance in ncoord dimensions with periodicity !! (Needs specification of domsize(ncoord)) !! 2: Aproximate geographic distance with horizontal coordinates in radians (-pi/+pi) !! 3: Geographic distance computation using haversine formula !! 10-13: Variants of distance types 0-3, but particularly for 3 dimensions in which !! a 2+1 dimensional localization is applied (distance weighting only in the horizontal) !! !! __Revision history:__ !! * 2024-02 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_check_dist2_noniso(thisobs, thisobs_l, coordsA, coordsB, distance2, & dists, cradius, sradius, checkdist, verbose, cnt_obs) IMPLICIT NONE ! *** Arguments *** TYPE(obs_f), INTENT(in) :: thisobs !< Data type with full observation TYPE(obs_l), INTENT(in) :: thisobs_l !< Data type with local observation REAL, INTENT(in) :: coordsA(:) !< Coordinates of current analysis domain (ncoord) REAL, INTENT(in) :: coordsB(:) !< Coordinates of observation (ncoord) REAL, INTENT(out) :: distance2 !< Squared distance REAL, INTENT(inout) :: dists(:) !< Vector of distance in each coordinate direction REAL, INTENT(out) :: cradius !< Directional cut-off radius REAL, INTENT(inout) :: sradius !< Directional support radius LOGICAL, INTENT(out) :: checkdist !< Flag whether distance is within cut-off radius INTEGER, INTENT(in) :: verbose !< Control screen output INTEGER, INTENT(inout) :: cnt_obs !< Count number of local observations ! *** Local variables *** INTEGER :: k ! Counters REAL :: slon, slat ! sine of distance in longitude or latitude INTEGER :: domsize ! Flag whether domainsize is set REAL :: cradius2 ! cut-off radius on ellipse or ellipsoid REAL :: phi, theta ! Angles in ellipse or ellipsoid REAL :: dist_xy ! Distance in xy-plan in 3D case LOGICAL :: distflag ! Flag whether distance in a coordinate direction is within cradius ! ********************** ! *** Initialization *** ! ********************** ! Initialize distance flag checkdist = .FALSE. ! Whether an observation lies within the local box distflag = .TRUE. ! Whether an observation lies within the local radius (ellipse, ellipsoid) ! ************************ ! *** Compute distance *** ! ************************ IF (.NOT.ALLOCATED(thisobs%domainsize)) THEN domsize = 0 ELSE domsize = 1 END IF ! Debug output IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' use non-isotropic localization' END IF norm: IF ((thisobs%disttype==0 .OR. thisobs%disttype==10) .OR. & ((thisobs%disttype==1 .OR. thisobs%disttype==11) .AND. domsize==0)) THEN ! *** Compute Cartesian distance *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' compute Cartesian distance' END IF IF (thisobs%ncoord==3) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) IF (dists(3)>thisobs_l%cradius(3)) THEN distflag = .FALSE. ELSE dists(2) = ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE dists(1) = ABS(coordsA(1) - coordsB(1)) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSEIF (thisobs%ncoord==2) THEN dists(2) = ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE dists(1) = ABS(coordsA(1) - coordsB(1)) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF (thisobs%ncoord==1) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF ((thisobs%disttype==1 .OR. thisobs%disttype==11) .AND. domsize==1) THEN norm ! *** Compute periodic Cartesian distance *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' compute periodic Cartesian distance' END IF IF (thisobs%ncoord==3) THEN IF (thisobs%domainsize(3)<=0.0) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) ELSE dists(3) = MIN(ABS(coordsA(3) - coordsB(3)), & ABS(ABS(coordsA(3) - coordsB(3))-thisobs%domainsize(3))) END IF IF (dists(3)>thisobs_l%cradius(3)) THEN distflag = .FALSE. ELSE IF (thisobs%domainsize(2)<=0.0) THEN dists(2) = ABS(coordsA(2) - coordsB(2)) ELSE dists(2) = MIN(ABS(coordsA(2) - coordsB(2)), & ABS(ABS(coordsA(2) - coordsB(2))-thisobs%domainsize(2))) END IF IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE IF (thisobs%domainsize(1)<=0.0) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) ELSE dists(1) = MIN(ABS(coordsA(1) - coordsB(1)), & ABS(ABS(coordsA(1) - coordsB(1))-thisobs%domainsize(1))) END IF IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSEIF (thisobs%ncoord==2) THEN IF (thisobs%domainsize(2)<=0.0) THEN dists(2) = ABS(coordsA(2) - coordsB(2)) ELSE dists(2) = MIN(ABS(coordsA(2) - coordsB(2)), & ABS(ABS(coordsA(2) - coordsB(2))-thisobs%domainsize(2))) END IF IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE IF (thisobs%domainsize(1)<=0.0) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) ELSE dists(1) = MIN(ABS(coordsA(1) - coordsB(1)), & ABS(ABS(coordsA(1) - coordsB(1))-thisobs%domainsize(1))) END IF IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF (thisobs%ncoord==1) THEN IF (thisobs%domainsize(1)<=0.0) THEN dists(1) = ABS(coordsA(1) - coordsB(1)) ELSE dists(1) = MIN(ABS(coordsA(1) - coordsB(1)), & ABS(ABS(coordsA(1) - coordsB(1))-thisobs%domainsize(1))) END IF IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF ELSEIF (thisobs%disttype==2 .OR. thisobs%disttype==12) THEN norm ! *** Compute distance from geographic coordinates *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' compute geographic distance' END IF IF (thisobs%ncoord==3) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) IF (dists(3)>thisobs_l%cradius(3)) THEN distflag = .FALSE. ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 1, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) IF (dists(1)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF ELSEIF (thisobs%disttype==3 .OR. thisobs%disttype==13) THEN norm ! *** Compute distance from geographic coordinates with haversine formula *** IF (debug>0 .AND. verbose==0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, & ' compute geographic distance using haversine function' END IF IF (thisobs%ncoord==3) THEN dists(3) = ABS(coordsA(3) - coordsB(3)) IF (dists(3)>thisobs_l%cradius(3)) THEN distflag = .FALSE. ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) ! Haversine formula slon = SIN((coordsA(1) - coordsB(1))/2) slat = SIN((coordsA(2) - coordsB(2))/2) dists(2) = SQRT(slat*slat + COS(coordsA(2))*COS(coordsB(2))*slon*slon) IF (dists(2)<=1.0) THEN dists(2) = 2.0 * r_earth* ASIN(dists(2)) ELSE dists(2) = r_earth* pi END IF IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 IF (thisobs%disttype<10) THEN ! full 3D localization DO k = 2, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO ELSE ! factorized 2+1D localization DO k = 2, thisobs%ncoord-1 distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF ELSE dists(2) = r_earth * ABS(coordsA(2) - coordsB(2)) IF (dists(2)>thisobs_l%cradius(2)) THEN distflag = .FALSE. ELSE dists(1) = r_earth * MIN( ABS(coordsA(1) - coordsB(1))* COS(coordsA(2)), & ABS(ABS(coordsA(1) - coordsB(1)) - 2.0*pi) * COS(coordsA(2))) ! Haversine formula slon = SIN((coordsA(1) - coordsB(1))/2) slat = SIN((coordsA(2) - coordsB(2))/2) dists(2) = SQRT(slat*slat + COS(coordsA(2))*COS(coordsB(2))*slon*slon) IF (dists(2)<=1.0) THEN dists(2) = 2.0 * r_earth* ASIN(dists(2)) ELSE dists(2) = r_earth* pi END IF IF (dists(2)>thisobs_l%cradius(1)) THEN distflag = .FALSE. ELSE ! full squared distance distance2 = 0.0 DO k = 1, thisobs%ncoord distance2 = distance2 + dists(k)*dists(k) END DO END IF END IF END IF END IF norm ! *************************************************************************** ! *** Compute directional cut-off and support radii and set distance flag *** ! *************************************************************************** dflag: IF (distflag) THEN nrad: IF (thisobs_l%nradii == 2 .OR. (thisobs_l%nradii == 3 .AND. thisobs%disttype >= 10)) THEN IF ((thisobs_l%cradius(1) == thisobs_l%cradius(2)) .OR. & (thisobs_l%sradius(1) == thisobs_l%sradius(2))) THEN ! 2D isotropic case cradius2 = thisobs_l%cradius(1) * thisobs_l%cradius(1) IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. cnt_obs = cnt_obs + 1 cradius = thisobs_l%cradius(1) sradius = thisobs_l%sradius(1) IF (debug>0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, & ' 2D isotropic with separately specified, but equal, radii' WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' theta, cradius, sradius', & theta*180/pi, cradius, sradius END IF END IF ELSE ! *** 2D anisotropic case: Use polar radius of ellipse in 2 dimensions *** ! Compute angle IF (dists(1) /= 0.0) THEN theta = ATAN(dists(2) / dists(1)) ELSE theta = pi / 2.0 END IF ! Compute radius in direction of theta IF (thisobs_l%cradius(1)>0.0 .OR. thisobs_l%cradius(2)>0.0) THEN cradius = thisobs_l%cradius(1) * thisobs_l%cradius(2) / & SQRT( (thisobs_l%cradius(2)*COS(theta))**2 & + (thisobs_l%cradius(1)*SIN(theta))**2 ) ELSE cradius = 0.0 END IF cradius2 = cradius * cradius IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. cnt_obs = cnt_obs + 1 ! Compute support radius in direction of theta IF (thisobs_l%sradius(1)>0.0 .OR. thisobs_l%sradius(2)>0.0) THEN sradius = thisobs_l%sradius(1) * thisobs_l%sradius(2) / & SQRT( (thisobs_l%sradius(2)*COS(theta))**2 & + (thisobs_l%sradius(1)*SIN(theta))**2 ) ELSE sradius = 0.0 END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, & ' 2D nonisotropic localization' WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' theta, cradius, sradius', & theta*180/pi, cradius, sradius END IF END IF END IF ELSE IF (thisobs_l%nradii == 3 .AND. thisobs%disttype < 10) THEN nrad ! To save computing time, we here distinguish whether ! - the horizontal radii are equal and only direction 3 has a different radius ! - whether all radii are equal (isotropic but specified with separate radii) ! - the anisotropy is in all 3 dimensions (all radii different) aniso: IF ((thisobs_l%cradius(1) == thisobs_l%cradius(2)) .AND. & (thisobs_l%cradius(1) /= thisobs_l%cradius(3)) .AND. & (thisobs_l%sradius(1) == thisobs_l%sradius(2))) THEN ! *** Isotropic in horizontal direction, distinct radius in the third direction (vertical) *** dist_xy = SQRT(dists(1)*dists(1) + dists(2)*dists(2)) ! 2D anisotropy: Polar radius of ellipse in 2 dimensions ! Compute angle IF (dist_xy /= 0.0) THEN theta = ATAN(dists(3) / dist_xy) ELSE theta = pi / 2.0 END IF ! Compute radius in direction of theta IF (thisobs_l%cradius(1)>0.0 .OR. thisobs_l%cradius(3)>0.0) THEN cradius = thisobs_l%cradius(1) * thisobs_l%cradius(3) / & SQRT( (thisobs_l%cradius(3)*COS(theta))**2 & + (thisobs_l%cradius(1)*SIN(theta))**2 ) ELSE cradius = 0.0 END IF cradius2 = cradius * cradius IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. cnt_obs = cnt_obs + 1 ! Compute support radius in direction of theta IF (thisobs_l%sradius(1)>0.0 .OR. thisobs_l%sradius(3)>0.0) THEN sradius = thisobs_l%sradius(1) * thisobs_l%sradius(3) / & SQRT( (thisobs_l%sradius(3)*COS(theta))**2 & + (thisobs_l%sradius(1)*SIN(theta))**2 ) ELSE sradius = 0.0 END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, & ' 3D: isotropic in directions 1 and 2, nonisotropic in direction 3' WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' theta, cradius, sradius', & theta*180/pi, cradius, sradius END IF END IF ELSEIF ((thisobs_l%cradius(1) == thisobs_l%cradius(2)) .AND. & (thisobs_l%cradius(1) == thisobs_l%cradius(3)) .AND. & (thisobs_l%sradius(1) == thisobs_l%sradius(2)) .AND. & (thisobs_l%sradius(2) == thisobs_l%sradius(3))) THEN aniso ! *** 3D isotropic case (all radii equal) *** cradius = thisobs_l%cradius(1) cradius2 = thisobs_l%cradius(1) * thisobs_l%cradius(1) sradius = thisobs_l%sradius(1) IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. cnt_obs = cnt_obs + 1 END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, & ' 3D isotropic case specified with vector of radii' WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' theta, cradius, sradius', & theta*180/pi, cradius, sradius END IF ELSE aniso ! *** general 3D anisotropic case *** ! Polar radius of ellipsoid in 3 dimensions ! Compute angle in x-y direction IF (dists(1) /= 0.0) THEN theta = ATAN(dists(2) / dists(1)) ELSE theta = pi / 2.0 END IF ! Distance in xy-plane dist_xy = SQRT(dists(1)**2 + dists(2)**2) ! Compute angle of xy-plane to z direction IF (dist_xy /= 0.0) THEN phi = ATAN(dists(3) / dist_xy) ELSE phi = 0.0 END IF ! Compute radius in direction of theta IF (thisobs_l%cradius(1)>0.0 .OR. thisobs_l%cradius(2)>0.0 .OR. thisobs_l%cradius(3)>0.0) THEN cradius = thisobs_l%cradius(1) * thisobs_l%cradius(2) * thisobs_l%cradius(3) / & SQRT( (thisobs_l%cradius(2)*thisobs_l%cradius(3)*COS(phi)*COS(theta))**2 & + (thisobs_l%cradius(1)*thisobs_l%cradius(3)*COS(phi)*SIN(theta))**2 & + (thisobs_l%cradius(1)*thisobs_l%cradius(2)*SIN(phi))**2 ) ELSE cradius = 0.0 END IF cradius2 = cradius * cradius IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. cnt_obs = cnt_obs + 1 ! Compute support radius in direction of theta IF (thisobs_l%sradius(1)>0.0 .OR. thisobs_l%sradius(2)>0.0 .OR. thisobs_l%sradius(3)>0.0) THEN sradius = thisobs_l%sradius(1) * thisobs_l%sradius(2) * thisobs_l%sradius(3) / & SQRT( (thisobs_l%sradius(2)*thisobs_l%sradius(3)*COS(phi)*COS(theta))**2 & + (thisobs_l%sradius(1)*thisobs_l%sradius(3)*COS(phi)*SIN(theta))**2 & + (thisobs_l%sradius(1)*thisobs_l%sradius(2)*SIN(phi))**2 ) ELSE sradius = 0.0 END IF IF (debug>0) THEN WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, & ' 3D nonisotropic localization' WRITE (*,*) '++ OMI-debug check_dist2_noniso: ', debug, ' theta, phi, distance, cradius, sradius', & theta*180/pi, phi*180/pi, SQRT(distance2), cradius, sradius END IF END IF END IF aniso ELSEIF (thisobs_l%nradii == 1) THEN nrad cradius = thisobs_l%cradius(1) cradius2 = thisobs_l%cradius(1) * thisobs_l%cradius(1) sradius = thisobs_l%sradius(1) IF (distance2 <= cradius2) THEN ! Set flag for valid observation checkdist = .TRUE. cnt_obs = cnt_obs + 1 END IF END IF nrad END IF dflag END SUBROUTINE PDAFomi_check_dist2_noniso !------------------------------------------------------------------------------- !> Compute weight vector for localization !! !! The routine computes a weight vector according to the !! distances of observations from the local analysis !! domain. !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_weights_l(verbose, nobs_l, ncols, locweight, cradius, sradius, & matA, ivar_obs_l, dist_l, weight_l) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: verbose !< Verbosity flag INTEGER, INTENT(in) :: nobs_l !< Number of local observations INTEGER, INTENT(in) :: ncols !< INTEGER, INTENT(in) :: locweight !< Localization weight type REAL, INTENT(in) :: cradius(:) !< Localization cut-off radius REAL, INTENT(in) :: sradius(:) !< support radius for weight functions REAL, INTENT(in) :: matA(:,:) !< REAL, INTENT(in) :: ivar_obs_l(:) !< Local vector of inverse obs. variances (nobs_l) REAL, INTENT(in) :: dist_l(:) !< Local vector of obs. distances (nobs_l) REAL, INTENT(out) :: weight_l(:) !< Output: vector of weights ! *** local variables *** INTEGER :: i ! Index of observation component INTEGER :: verbose_w ! Verbosity flag for weight computation INTEGER :: wtype ! Type of weight function INTEGER :: rtype ! Type of weight regulation REAL :: var_obs_l ! Variance of observation error REAL, ALLOCATABLE :: A_obs(:,:) ! Array for a single row of matA ! ********************** ! *** INITIALIZATION *** ! ********************** ! Screen output IF (verbose == 1) THEN IF (locweight == 1 .OR. locweight == 2 .OR. locweight == 3 & .OR. locweight == 4 .OR. locweight == 5. .OR. locweight == 15 & .OR. locweight == 16) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use distance-dependent weight for observation errors' IF (locweight == 3 .OR. locweight == 15) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use regulated weight with mean error variance' ELSE IF (locweight == 4 .OR. locweight == 16) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use regulated weight with single-point error variance' END IF ELSE IF (locweight == 11 .OR. locweight == 26 .OR. locweight == 27) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use distance-dependent weight for observed ensemble' END IF ENDIF ! *** Initialize weight vector IF (locweight == 0) THEN ! Uniform (unit) weighting wtype = 0 rtype = 0 ELSE IF (locweight == 1 .OR. locweight == 11) THEN ! Exponential weighting wtype = 1 rtype = 0 ELSE IF (locweight == 2 .OR. locweight == 3 .OR. locweight == 4 & .OR. locweight == 16 .OR. locweight == 17) THEN ! 5th-order polynomial (Gaspari&Cohn, 1999) wtype = 2 IF (locweight == 3 .OR. locweight == 4) THEN ! Use regulated weight rtype = 1 ELSE ! No regulated weight rtype = 0 END IF ELSE IF (locweight == 5 .OR. locweight == 15 .OR. locweight ==16) THEN ! 5th-order polynomial (Gaspari&Cohn, 1999) wtype = 3 IF (locweight == 15 .OR. locweight == 16) THEN ! Use regulated weight rtype = 1 ELSE ! No regulated weight rtype = 0 END IF END IF IF (locweight == 4 .OR. locweight == 16) THEN ! Allocate array for single observation point ALLOCATE(A_obs(1, ncols)) END IF ! Control verbosity of PDAF_local_weight IF (verbose==1) verbose_w = 1 IF (locweight /= 4) THEN ! All localizations except regulated weight based on variance at ! single observation point DO i = 1, nobs_l ! set observation variance value var_obs_l = 1.0 / ivar_obs_l(i) CALL PDAF_local_weight(wtype, rtype, cradius(i), sradius(i), dist_l(i), & nobs_l, ncols, matA, var_obs_l, weight_l(i), verbose_w) verbose_w = 0 END DO ELSE ! Regulated weight using variance at single observation point DO i = 1, nobs_l ! set observation variance value var_obs_l = 1.0 / ivar_obs_l(i) A_obs(1,:) = matA(i,:) CALL PDAF_local_weight(wtype, rtype, cradius(i), sradius(i), dist_l(i), & 1, ncols, A_obs, var_obs_l, weight_l(i), verbose_w) verbose_w = 0 END DO END IF ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'thisobs_l%dim_obs_l', nobs_l WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'thisobs%distance_l', dist_l(1:nobs_l) WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'weight_l', weight_l(1:nobs_l) WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'thisobs_l%ivar_obs_l', ivar_obs_l(1:nobs_l) END IF IF (locweight == 4) DEALLOCATE(A_obs) END SUBROUTINE PDAFomi_weights_l !------------------------------------------------------------------------------- !> Compute weight vector for localization !! !! The routine computes a weight vector according to the !! distances of observations from the local analysis !! domain. !! !! __Revision history:__ !! * 2020-03 - Lars Nerger - Initial code from restructuring observation routines !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_weights_l_sgnl(verbose, nobs_l, ncols, locweight, cradius, sradius, & matA, ivar_obs_l, dist_l, weight_l) IMPLICIT NONE ! *** Arguments *** INTEGER, INTENT(in) :: verbose !< Verbosity flag INTEGER, INTENT(in) :: nobs_l !< Number of local observations INTEGER, INTENT(in) :: ncols !< INTEGER, INTENT(in) :: locweight !< Localization weight type REAL, INTENT(in) :: cradius !< Localization cut-off radius REAL, INTENT(in) :: sradius !< support radius for weight functions REAL, INTENT(in) :: matA(:,:) !< REAL, INTENT(in) :: ivar_obs_l(:) !< Local vector of inverse obs. variances (nobs_l) REAL, INTENT(in) :: dist_l(:) !< Local vector of obs. distances (nobs_l) REAL, INTENT(out) :: weight_l(:) !< Output: vector of weights ! *** local variables *** INTEGER :: i ! Index of observation component INTEGER :: verbose_w ! Verbosity flag for weight computation INTEGER :: wtype ! Type of weight function INTEGER :: rtype ! Type of weight regulation REAL :: var_obs_l ! Variance of observation error REAL, ALLOCATABLE :: A_obs(:,:) ! Array for a single row of matA ! ********************** ! *** INITIALIZATION *** ! ********************** ! Screen output IF (verbose == 1) THEN IF (locweight == 1 .OR. locweight == 2 .OR. locweight == 3 & .OR. locweight == 4 .OR. locweight == 5. .OR. locweight == 15 & .OR. locweight == 16) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use distance-dependent weight for observation errors' IF (locweight == 3 .OR. locweight == 15) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use regulated weight with mean error variance' ELSE IF (locweight == 4 .OR. locweight == 16) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use regulated weight with single-point error variance' END IF ELSE IF (locweight == 11 .OR. locweight == 26 .OR. locweight == 27) THEN WRITE (*, '(a, 8x, a)') & 'PDAFomi', '--- Use distance-dependent weight for observed ensemble' END IF ENDIF ! *** Initialize weight vector IF (locweight == 0) THEN ! Uniform (unit) weighting wtype = 0 rtype = 0 ELSE IF (locweight == 1 .OR. locweight == 11) THEN ! Exponential weighting wtype = 1 rtype = 0 ELSE IF (locweight == 2 .OR. locweight == 3 .OR. locweight == 4 & .OR. locweight == 16 .OR. locweight == 17) THEN ! 5th-order polynomial (Gaspari&Cohn, 1999) wtype = 2 IF (locweight == 3 .OR. locweight == 4) THEN ! Use regulated weight rtype = 1 ELSE ! No regulated weight rtype = 0 END IF ELSE IF (locweight == 5 .OR. locweight == 15 .OR. locweight ==16) THEN ! 5th-order polynomial (Gaspari&Cohn, 1999) wtype = 3 IF (locweight == 15 .OR. locweight == 16) THEN ! Use regulated weight rtype = 1 ELSE ! No regulated weight rtype = 0 END IF END IF IF (locweight == 4 .OR. locweight == 16) THEN ! Allocate array for single observation point ALLOCATE(A_obs(1, ncols)) END IF ! Control verbosity of PDAF_local_weight IF (verbose==1) verbose_w = 1 IF (locweight /= 4) THEN ! All localizations except regulated weight based on variance at ! single observation point DO i = 1, nobs_l ! set observation variance value var_obs_l = 1.0 / ivar_obs_l(i) CALL PDAF_local_weight(wtype, rtype, cradius, sradius, dist_l(i), & nobs_l, ncols, matA, var_obs_l, weight_l(i), verbose_w) verbose_w = 0 END DO ELSE ! Regulated weight using variance at single observation point DO i = 1, nobs_l ! set observation variance value var_obs_l = 1.0 / ivar_obs_l(i) A_obs(1,:) = matA(i,:) CALL PDAF_local_weight(wtype, rtype, cradius, sradius, dist_l(i), & 1, ncols, A_obs, var_obs_l, weight_l(i), verbose_w) verbose_w = 0 END DO END IF ! Print debug information IF (debug>0) THEN WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'thisobs_l%dim_obs_l', nobs_l WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'thisobs%distance_l', dist_l(1:nobs_l) WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'weight_l', weight_l(1:nobs_l) WRITE (*,*) '++ OMI-debug weights_l: ', debug, 'thisobs_l%ivar_obs_l', ivar_obs_l(1:nobs_l) END IF IF (locweight == 4) DEALLOCATE(A_obs) END SUBROUTINE PDAFomi_weights_l_sgnl !------------------------------------------------------------------------------- !> Deallocate arrays in observation type !! !! This routine deallocates arrays in the data type THISOBS. !! The routine mainly operates on the full observation type. !! It is included here to avoid cross-dependences between !! PDAFomi_obs_f and PDAFomi_obs_l. !! !! The routine is called by all filter processes. !! !! __Revision history:__ !! * 2019-10 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_deallocate_obs(thisobs) USE PDAFomi_obs_f, & ONLY: obs_f, n_obstypes, obscnt, offset_obs, obs_f_all, & offset_obs_g, obsdims, map_obs_id IMPLICIT NONE ! *** Arguments TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation ! *** Perform deallocation *** IF (ALLOCATED(thisobs%obs_f)) DEALLOCATE(thisobs%obs_f) IF (ALLOCATED(thisobs%ocoord_f)) DEALLOCATE(thisobs%ocoord_f) IF (ALLOCATED(thisobs%id_obs_p)) DEALLOCATE(thisobs%id_obs_p) IF (ALLOCATED(thisobs%ivar_obs_f)) DEALLOCATE(thisobs%ivar_obs_f) IF (ALLOCATED(thisobs%icoeff_p)) DEALLOCATE(thisobs%icoeff_p) IF (ALLOCATED(thisobs%domainsize)) DEALLOCATE(thisobs%domainsize) IF (ALLOCATED(thisobs%id_obs_f_lim)) DEALLOCATE(thisobs%id_obs_f_lim) ! Reset assim flag thisobs%doassim = 0 IF (ALLOCATED(obs_f_all)) DEALLOCATE(obs_f_all) IF (ALLOCATED(obsdims)) DEALLOCATE(obsdims) IF (ALLOCATED(map_obs_id)) DEALLOCATE(map_obs_id) ! Reset counters over all observation types n_obstypes = 0 obscnt = 0 offset_obs = 0 offset_obs_g = 0 ! Reset flag for first call to local observations firstobs = 0 END SUBROUTINE PDAFomi_deallocate_obs !------------------------------------------------------------------------------- !> Exclude observations for too high innovation !! !! The routine is called during the analysis step !! on each local analysis domain. 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:__ !! * 2022-12 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_omit_by_inno_l(thisobs_l, thisobs, inno_l, obs_l_all, obsid, cnt_all, verbose) IMPLICIT NONE ! *** Arguments *** TYPE(obs_l), INTENT(inout) :: thisobs_l !< Data type with local observation TYPE(obs_f), INTENT(inout) :: thisobs !< Data type with full observation REAL, INTENT(in) :: inno_l(:) !< Input vector of observation innovation REAL, INTENT(in) :: obs_l_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 INTEGER, INTENT(in) :: verbose !< Verbosity flag ! *** 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_l -- START obs-ID', obsid WRITE (*,*) '++ OMI-debug omit_by_inno_l:', debug, 'limit for innovation', & thisobs%inno_omit WRITE (*,*) '++ OMI-debug omit_by_inno_l:', debug, 'inno_omit_ivar', & thisobs%inno_omit_ivar WRITE (*,*) '++ OMI-debug omit_by_inno_l:', debug, 'innovation_l', inno_l ENDIF ! Squared limit factor limit2 = thisobs%inno_omit * thisobs%inno_omit ! Check for observations to be excluded cnt = 0 DO i = 1, thisobs_l%dim_obs_l ! Squared innovation inno2 = inno_l(i + thisobs_l%off_obs_l)* inno_l(i + thisobs_l%off_obs_l) IF (inno2 > limit2 * 1.0/thisobs_l%ivar_obs_l(i)) THEN IF (debug>0) THEN WRITE (*,*) '++ OMI-debug omit_by_inno_l:', debug, 'omit: innovation:', & inno_l(i + thisobs_l%off_obs_l), 'observation:', obs_l_all(i + thisobs_l%off_obs_l) END IF ! Exclude observation by increased its observation error thisobs_l%ivar_obs_l(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_l:', debug, 'count of excluded obs.: ', cnt WRITE (*,*) '++ OMI-debug omit_by_inno_l:', debug, 'updated thisobs_l%ivar_obs_l ', & thisobs_l%ivar_obs_l ENDIF IF (debug>0) & WRITE (*,*) '++ OMI-debug: ', debug, 'PDAFomi_omit_by_inno_l -- END obs-ID', obsid cnt_all = cnt_all + cnt END IF ENDIF doassim END SUBROUTINE PDAFomi_omit_by_inno_l !------------------------------------------------------------------------------- !> Get statistics on local observations !! !! The routine is called in the update routine of !! local filters and write statistics on locally !! used and excluded observations. !! !! __Revision history:__ !! * 2023-03 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_obsstats_l(screen) USE MPI USE PDAFomi_obs_f, ONLY: ostats_omit USE PDAF_mod_filtermpi, & ONLY: COMM_filter, MPIerr 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, 5, MPI_INTEGER, MPI_SUM, & 0, COMM_filter, MPIerr) CALL MPI_Reduce(ostats_omit(6:7), ostats_omit_g(6:7), 2, MPI_INTEGER, MPI_MAX, & 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, 5x, a)') 'PDAFomi', '--- Global statistics for locally omitted observations:' WRITE (*, '(a, 8x, a, i10)') & 'PDAFomi', 'Local domains with omitted observations: ', ostats_omit_g(1) WRITE (*, '(a, 8x, a, i10)') & 'PDAFomi', 'Local domains without omitted observations: ', ostats_omit_g(2) WRITE (*, '(a, 8x, a, i10)') & 'PDAFomi', 'Maximum number of locally omitted observations: ', ostats_omit_g(6) WRITE (*, '(a, 8x, a, i10)') & 'PDAFomi', 'Maximum number of locally used observations: ', ostats_omit_g(7) IF (ostats_omit_g(2) > 0) THEN WRITE (*, '(a, 8x, a, f9.1)') & 'PDAFomi', 'Total avg. locally omitted observations:', & REAL(ostats_omit_g(3)) / REAL(ostats_omit_g(1) + ostats_omit_g(2)) WRITE (*, '(a, 8x, a, f9.1)') & 'PDAFomi', 'Total avg. locally used observations: ', & REAL(ostats_omit_g(4)) / REAL(ostats_omit_g(1) + ostats_omit_g(2)) WRITE (*, '(a, 8x, a, f9.1)') & 'PDAFomi', 'Avg. omitted for domains with omitted observations: ', & REAL(ostats_omit_g(3)) / REAL(ostats_omit_g(1)) WRITE (*, '(a, 8x, a, f9.1)') & 'PDAFomi', 'Avg. used for domains with omitted observations: ', & REAL(ostats_omit_g(5)) / REAL(ostats_omit_g(1)) END IF ELSEIF (mype == 0 .AND. screen > 0 .AND. ostats_omit_g(2)>0) THEN WRITE (*, '(a, 5x, a)') 'PDAFomi', '--- Global statistics for locally omitted observations:' WRITE (*, '(a, 9x, a)') & 'PDAFomi', 'Zero observations omitted' END IF END SUBROUTINE PDAFomi_obsstats_l !------------------------------------------------------------------------------- !> Deallocate arrays in all observation types !! !! This routine deallocates arrays in all observation types. !! The routine is only called internally in PDAF. !! !! The routine is called by all filter processes. !! !! __Revision history:__ !! * 2021-04 - Lars Nerger - Initial code !! * Later revisions - see repository log !! SUBROUTINE PDAFomi_dealloc() USE PDAFomi_obs_f, & ONLY: obs_f, n_obstypes, obscnt, offset_obs, obs_f_all, & offset_obs_g, obsdims, map_obs_id ! *** Local variables INTEGER :: i ! *** Perform deallocation *** IF (n_obstypes>0) THEN DO i=1, n_obstypes IF (ALLOCATED(obs_f_all(i)%ptr%obs_f)) DEALLOCATE(obs_f_all(i)%ptr%obs_f) IF (ALLOCATED(obs_f_all(i)%ptr%ocoord_f)) DEALLOCATE(obs_f_all(i)%ptr%ocoord_f) IF (ALLOCATED(obs_f_all(i)%ptr%id_obs_p)) DEALLOCATE(obs_f_all(i)%ptr%id_obs_p) IF (ALLOCATED(obs_f_all(i)%ptr%ivar_obs_f)) DEALLOCATE(obs_f_all(i)%ptr%ivar_obs_f) IF (ALLOCATED(obs_f_all(i)%ptr%icoeff_p)) DEALLOCATE(obs_f_all(i)%ptr%icoeff_p) IF (ALLOCATED(obs_f_all(i)%ptr%ivar_obs_f)) DEALLOCATE(obs_f_all(i)%ptr%ivar_obs_f) IF (ALLOCATED(obs_f_all(i)%ptr%icoeff_p)) DEALLOCATE(obs_f_all(i)%ptr%icoeff_p) IF (ALLOCATED(obs_f_all(i)%ptr%domainsize)) DEALLOCATE(obs_f_all(i)%ptr%domainsize) IF (ALLOCATED(obs_f_all(i)%ptr%id_obs_f_lim)) DEALLOCATE(obs_f_all(i)%ptr%id_obs_f_lim) ! Reset assim flag obs_f_all(i)%ptr%doassim = 0 END DO IF (ALLOCATED(obs_f_all)) DEALLOCATE(obs_f_all) IF (ALLOCATED(obsdims)) DEALLOCATE(obsdims) IF (ALLOCATED(map_obs_id)) DEALLOCATE(map_obs_id) ! Reset counters over all observation types n_obstypes = 0 obscnt = 0 offset_obs = 0 offset_obs_g = 0 ! Reset flag for first call to local observations firstobs = 0 END IF END SUBROUTINE PDAFomi_dealloc END MODULE PDAFomi_obs_l