MODULE extract_obs_mod ! !git $Id$ !svn $Id: extract_obs.F 1151 2023-02-09 03:08:53Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! This routine extracts model data at the requested observations ! ! positions (Xobs,Yobs,Zobs). The extraction is done via linear ! ! interpolation. The (Xobs,Yobs) positions must be in fractional ! ! grid coordinates. Zobs can be in fractional grid coordinates ! ! (Zobs >= 0) or actual depths (Zobs < 0), if applicable. ! ! ! ! On Input: ! ! ! ! ng Nested grid number. ! ! Imin Global I-coordinate lower bound threshold for ! ! requested state field type. ! ! Imax Global I-coordinate upper bound threshold for ! ! requested state field type. ! ! Jmin Global J-coordinate lower bound threshold for ! ! requested state field type. ! ! Jmax Global J-coordinate upper bound threshold for ! ! requested state field type. ! ! LBi I-dimension Lower bound. ! ! UBi I-dimension Upper bound. ! ! LBj J-dimension Lower bound. ! ! UBj J-dimension Upper bound. ! ! LBk K-dimension Lower bound. ! ! UBk K-dimension Upper bound. ! ! ifield State field identification to process. ! ! Mobs Observation dimension in the calling program. ! ! NobsSTR Starting observation to process. ! ! NobsEND Last observations to process. ! ! Xmin Global minimum fractional I-coordinate to consider. ! ! Xmax Global maximum fractional I-coordinate to consider. ! ! Ymin Global minimum fractional J-coordinate to consider. ! ! Ymax Global maximum fractional J-coordinate to consider. ! ! time Current model time (seconds). ! ! dt Model baroclinic time-step (seconds). ! ! ObsType Observations type. ! ! Tobs Observations time (days). ! ! Xobs Observations X-locations (grid coordinates). ! ! Yobs Observations Y-locations (grid coordinates). ! ! Zobs Observations Z-locations (grid coordinates or meters).! ! A Model array (2D or 3D) to process. ! ! Adepth Depths (meter; negative). ! ! Amask Land-sea masking. ! ! ! ! On Output: ! ! ! ! ObsVetting Observation screenning flag, 0: reject or 1: accept. ! ! Aobs Extracted model values at observation positions. ! ! Zobs Observations Z-locations (grid coordinates). ! ! ! ! The interpolation weights matrix, Hmat(1:8), is as follows: ! ! ! ! 8____________7 ! ! /. /| (i2,j2,k2) ! ! / . / | ! ! 5/___________/6 | ! ! | . | | ! ! | . | | Grid Cell ! ! | 4.........|..|3 ! ! | . | / ! ! |. | / ! ! (i1,j1,k1) |___________|/ ! ! 1 2 ! ! ! ! Notice that the indices i2 and j2 are reset when observations are ! ! located exactly at the eastern and/or northern boundaries. This is ! ! needed to avoid out-of-range array computations. ! ! ! ! All the observations are assumed to in fractional coordinates with ! ! respect to RHO-points: ! ! ! ! ! ! M r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r ! ! : : ! ! Mm+.5 v p++v++p++v++p++v++p++v++p++v++p++v++p++v++p++v++p v ! ! : + | | | | | | | + : ! ! Mm r u r u r u r u r u r u r u r u r u r ! ! : + | | | | | | | + : ! ! Mm-.5 v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v ! ! : + | | | | | | | + : ! ! r u r u r u r u r u r u r u r u r u r ! ! : + | | | | | | | + : ! ! v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v ! ! : + | | | | | | | + : ! ! r u r u r u r u r u r u r u r u r u r ! ! : + | | | | | | | + : ! ! 2.5 v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v ! ! : + | | | | | | | + : ! ! 2.0 r u r u r u r u r u r u r u r u r u r ! ! : + | | | | | | | + : ! ! 1.5 v p--v--p--v--p--v--p--v--p--v--p--v--p--v--p--v--p v ! ! : + | | | | | | | + : ! ! 1.0 r u r u r u r u r u r u r u r u r u r ! ! : + | | | | | | | + : ! ! 0.5 v p++v++p++v++p++v++p++v++p++v++p++v++p++v++p++v++p v ! ! : : ! ! 0.0 r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r..u..r ! ! ! ! 0.5 1.5 2.5 Lm-.5 Lm+.5 ! ! ! ! 0.0 1.0 2.0 Lm L ! ! ! !======================================================================= ! USE mod_kinds implicit none PUBLIC extract_obs2d PUBLIC extract_obs3d CONTAINS ! !*********************************************************************** SUBROUTINE extract_obs2d (ng, Imin, Imax, Jmin, Jmax, & & LBi, UBi, LBj, UBj, & & ifield, Mobs, NobsSTR, NobsEND, & & Xmin, Xmax, Ymin, Ymax, & & time, dt, & & ObsType, ObsVetting, & & Tobs, Xobs, Yobs, & & A, & & Amask, & & Aobs) !*********************************************************************** ! USE mod_ncparam, ONLY : isUbar, isVbar USE mod_fourdvar, ONLY : ObsState2Type ! ! Imported variable declarations. ! integer, intent(in) :: ng, Imin, Imax, Jmin, Jmax integer, intent(in) :: LBi, UBi, LBj, UBj integer, intent(in) :: ifield, Mobs, NobsSTR, NobsEND ! real(r8), intent(in) :: Xmin, Xmax, Ymin, Ymax real(dp), intent(in) :: time, dt ! integer, intent(in) :: ObsType(:) real(dp), intent(in) :: Tobs(:) real(r8), intent(in) :: Xobs(:) real(r8), intent(in) :: Yobs(:) real(r8), intent(in) :: A(LBi:,LBj:) real(r8), intent(in) :: Amask(LBi:,LBj:) real(r8), intent(inout) :: ObsVetting(:) real(r8), intent(inout) :: Aobs(:) ! ! Local variable declarations. ! integer :: ic, iobs, i1, i2, j1, j2 real(dp) :: TimeLB, TimeUB real(r8) :: p1, p2, q1, q2, wsum real(r8), dimension(8) :: Hmat ! !----------------------------------------------------------------------- ! Interpolate from requested 2D state field when appropriate. !----------------------------------------------------------------------- ! TimeLB=(time-0.5_dp*dt)/86400.0_dp TimeUB=(time+0.5_dp*dt)/86400.0_dp ! DO iobs=NobsSTR,NobsEND IF ((ObsType(iobs).eq.ifield).and. & & ((TimeLB.le.Tobs(iobs)).and.(Tobs(iobs).lt.TimeUB)).and. & & ((Xmin.le.Xobs(iobs)).and.(Xobs(iobs).lt.Xmax)).and. & & ((Ymin.le.Yobs(iobs)).and.(Yobs(iobs).lt.Ymax))) THEN IF (ObsType(iobs).eq.ObsState2Type(isUbar)) THEN i1=INT(Xobs(iobs)+0.5_r8) ! 2D U-grid type variable j1=INT(Yobs(iobs)) ELSE IF (ObsType(iobs).eq.ObsState2Type(isVbar)) THEN i1=INT(Xobs(iobs)) ! 2D V-grid type variable j1=INT(Yobs(iobs)+0.5_r8) ELSE i1=INT(Xobs(iobs)) ! 2D RHO-grid type variable j1=INT(Yobs(iobs)) END IF i2=i1+1 j2=j1+1 IF (i2.gt.Imax) THEN i2=i1 ! Observation at the eastern boundary END IF IF (j2.gt.Jmax) THEN j2=j1 ! Observation at the northern boundary END IF p2=REAL(i2-i1,r8)*(Xobs(iobs)-REAL(i1,r8)) q2=REAL(j2-j1,r8)*(Yobs(iobs)-REAL(j1,r8)) p1=1.0_r8-p2 q1=1.0_r8-q2 Hmat(1)=p1*q1 Hmat(2)=p2*q1 Hmat(3)=p2*q2 Hmat(4)=p1*q2 Hmat(1)=Hmat(1)*Amask(i1,j1) Hmat(2)=Hmat(2)*Amask(i2,j1) Hmat(3)=Hmat(3)*Amask(i2,j2) Hmat(4)=Hmat(4)*Amask(i1,j2) wsum=0.0_r8 DO ic=1,4 wsum=wsum+Hmat(ic) END DO IF (wsum.gt.0.0_r8) THEN wsum=1.0_r8/wsum DO ic=1,4 Hmat(ic)=Hmat(ic)*wsum END DO END IF Aobs(iobs)=Hmat(1)*A(i1,j1)+ & & Hmat(2)*A(i2,j1)+ & & Hmat(3)*A(i2,j2)+ & & Hmat(4)*A(i1,j2) IF (wsum.gt.0.0_r8) ObsVetting(iobs)=1.0_r8 END IF END DO RETURN END SUBROUTINE extract_obs2d ! !*********************************************************************** SUBROUTINE extract_obs3d (ng, Imin, Imax, Jmin, Jmax, & & LBi, UBi, LBj, UBj, LBk, UBk, & & ifield, Mobs, NobsSTR, NobsEND, & & Xmin, Xmax, Ymin, Ymax, & & time, dt, & & ObsType, ObsVetting, & & Tobs, Xobs, Yobs, Zobs, & & A, Adepth, & & Amask, & & Aobs) !*********************************************************************** ! USE mod_param ! USE mod_ncparam, ONLY : isUvel, isVvel USE mod_fourdvar, ONLY : ObsState2Type ! ! Imported variable declarations. ! integer, intent(in) :: ng, Imin, Imax, Jmin, Jmax integer, intent(in) :: LBi, UBi, LBj, UBj, LBk, UBk integer, intent(in) :: ifield, Mobs, NobsSTR, NobsEND ! real(r8), intent(in) :: Xmin, Xmax, Ymin, Ymax real(dp), intent(in) :: time, dt ! integer, intent(in) :: ObsType(:) real(dp), intent(in) :: Tobs(:) real(r8), intent(in) :: Xobs(:) real(r8), intent(in) :: Yobs(:) real(r8), intent(in) :: A(LBi:,LBj:,LBk:) real(r8), intent(in) :: Adepth(LBi:,LBj:,LBk:) real(r8), intent(in) :: Amask(LBi:,LBj:) real(r8), intent(inout) :: ObsVetting(:) real(r8), intent(inout) :: Zobs(:) real(r8), intent(inout) :: Aobs(:) ! ! Local variable declarations. ! integer :: i, ic, iobs, i1, i2, j1, j2, k, k1, k2 real(dp) :: TimeLB, TimeUB real(r8) :: Zbot, Ztop, dz, p1, p2, q1, q2, r1, r2 real(r8) :: w11, w12, w21, w22, wsum real(r8), dimension(8) :: Hmat ! !----------------------------------------------------------------------- ! Interpolate from requested 3D state field. !----------------------------------------------------------------------- ! TimeLB=(time-0.5_dp*dt)/86400.0_dp TimeUB=(time+0.5_dp*dt)/86400.0_dp ! DO iobs=NobsSTR,NobsEND IF ((ObsType(iobs).eq.ifield).and. & & ((TimeLB.le.Tobs(iobs)).and.(Tobs(iobs).lt.TimeUB)).and. & & ((Xmin.le.Xobs(iobs)).and.(Xobs(iobs).lt.Xmax)).and. & & ((Ymin.le.Yobs(iobs)).and.(Yobs(iobs).lt.Ymax))) THEN IF (ObsType(iobs).eq.ObsState2Type(isUvel)) THEN i1=INT(Xobs(iobs)+0.5_r8) ! 3D U-grid type variable j1=INT(Yobs(iobs)) ELSE IF (ObsType(iobs).eq.ObsState2Type(isVvel)) THEN i1=INT(Xobs(iobs)) ! 3D V-grid type variable j1=INT(Yobs(iobs)+0.5_r8) ELSE i1=INT(Xobs(iobs)) ! 3D RHO-grid type variable j1=INT(Yobs(iobs)) END IF i2=i1+1 j2=j1+1 IF (i2.gt.Imax) THEN i2=i1 ! Observation at the eastern boundary END IF IF (j2.gt.Jmax) THEN j2=j1 ! Observation at the northern boundary END IF p2=REAL(i2-i1,r8)*(Xobs(iobs)-REAL(i1,r8)) q2=REAL(j2-j1,r8)*(Yobs(iobs)-REAL(j1,r8)) p1=1.0_r8-p2 q1=1.0_r8-q2 w11=p1*q1 w21=p2*q1 w22=p2*q2 w12=p1*q2 IF (Zobs(iobs).gt.0.0_r8) THEN k1=MAX(1,INT(Zobs(iobs))) ! Positions in fractional k2=MIN(INT(Zobs(iobs))+1,N(ng)) ! levels r2=REAL(k2-k1,r8)*(Zobs(iobs)-REAL(k1,r8)) r1=1.0_r8-r2 ELSE Ztop=Adepth(i1,j1,N(ng)) Zbot=Adepth(i1,j1,1 ) IF (Zobs(iobs).ge.Ztop) THEN k1=N(ng) ! If shallower, assign to k2=N(ng) ! top grid cell. The r1=1.0_r8 ! observation is located r2=0.0_r8 ! on the upper cell half Zobs(iobs)=REAL(N(ng),r8) ! above its middle depth. ELSE IF (Zbot.ge.Zobs(iobs)) THEN r1=0.0_r8 ! If deeper, ignore. r2=0.0_r8 ObsVetting(iobs)=0.0_r8 ELSE DO k=N(ng),2,-1 ! Otherwise, interpolate Ztop=Adepth(i1,j1,k ) ! to fractional level Zbot=Adepth(i1,j1,k-1) IF ((Ztop.gt.Zobs(iobs)).and.(Zobs(iobs).ge.Zbot)) THEN k1=k-1 k2=k END IF END DO dz=Adepth(i1,j1,k2)-Adepth(i1,j1,k1) r2=(Zobs(iobs)-Adepth(i1,j1,k1))/dz r1=1.0_r8-r2 Zobs(iobs)=REAL(k1,r8)+r2 ! overwrite END IF END IF IF ((r1+r2).gt.0.0_r8) THEN Hmat(1)=w11*r1 Hmat(2)=w21*r1 Hmat(3)=w22*r1 Hmat(4)=w12*r1 Hmat(5)=w11*r2 Hmat(6)=w21*r2 Hmat(7)=w22*r2 Hmat(8)=w12*r2 Hmat(1)=Hmat(1)*Amask(i1,j1) Hmat(2)=Hmat(2)*Amask(i2,j1) Hmat(3)=Hmat(3)*Amask(i2,j2) Hmat(4)=Hmat(4)*Amask(i1,j2) Hmat(5)=Hmat(5)*Amask(i1,j1) Hmat(6)=Hmat(6)*Amask(i2,j1) Hmat(7)=Hmat(7)*Amask(i2,j2) Hmat(8)=Hmat(8)*Amask(i1,j2) wsum=0.0_r8 DO ic=1,8 wsum=wsum+Hmat(ic) END DO IF (wsum.gt.0.0_r8) THEN wsum=1.0_r8/wsum DO ic=1,8 Hmat(ic)=Hmat(ic)*wsum END DO END IF Aobs(iobs)=Hmat(1)*A(i1,j1,k1)+ & & Hmat(2)*A(i2,j1,k1)+ & & Hmat(3)*A(i2,j2,k1)+ & & Hmat(4)*A(i1,j2,k1)+ & & Hmat(5)*A(i1,j1,k2)+ & & Hmat(6)*A(i2,j1,k2)+ & & Hmat(7)*A(i2,j2,k2)+ & & Hmat(8)*A(i1,j2,k2) IF (wsum.gt.0.0_r8) ObsVetting(iobs)=1.0_r8 ! ! Reject observations that lie in the lower bottom grid cell (k=1) to ! avoid clustering due shallowing of bathymetry during smoothing and ! coarse level half-thickness (-h < Zobs < Adepth(:,:,1)) in deep ! water. ! IF ((Zobs(iobs).gt.0.0_r8).and.(Zobs(iobs).le.1.0_r8)) THEN ObsVetting(iobs)=0.0_r8 END IF END IF END IF END DO ! RETURN END SUBROUTINE extract_obs3d END MODULE extract_obs_mod