#include "cppdefs.h" #if (defined FOUR_DVAR || defined VERIFICATION) && defined OBSERVATIONS SUBROUTINE obs_read (ng, model, backward) ! !svn $Id: obs_read.F 906 2018-06-21 03:29:09Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2019 The ROMS/TOMS Group ! ! Licensed under a MIT/X style license ! ! See License_ROMS.txt ! !======================================================================= ! ! ! This subroutine reads in observations data when appropriate from ! ! observations input NetCDF file. The observations data is stored ! ! for use elsewhere. ! ! ! !======================================================================= ! USE mod_param USE mod_parallel USE mod_fourdvar USE mod_iounits USE mod_ncparam USE mod_netcdf USE mod_scalars ! USE dateclock_mod, ONLY : time_string USE strings_mod, ONLY : FoundError ! implicit none ! ! Imported variable declarations. ! logical, intent(in) :: backward integer, intent(in) :: ng, model ! ! Local variable declarations. ! logical :: readNLmod, readTLmod integer :: Mstr, Mend integer :: i, iobs, itrc, status character (len=22) :: t_code ! SourceFile=__FILE__ ! !--------------------------------------------------------------------- ! Read observation variables needed for interpolating the model ! state at the observation locations. !--------------------------------------------------------------------- ! IF (ProcessObs(ng)) THEN # if defined TLM_OBS # if defined IS4DVAR_SENSITIVITY readNLmod=.FALSE. # else readNLmod=.TRUE. # endif readTLmod=.TRUE. # else readNLmod=.FALSE. readTLmod=.FALSE. # endif ! ! Initialize observations processing counters. ! DO i=1,NobsVar(ng) FOURDVAR(ng)%ObsCount(i)=0 FOURDVAR(ng)%ObsReject(i)=0 END DO IF (backward) THEN ObsSurvey(ng)=ObsSurvey(ng)-1 ELSE ObsSurvey(ng)=ObsSurvey(ng)+1 END IF ! ! Set number of observations to process. ! Nobs(ng)=FOURDVAR(ng)%NobsSurvey(ObsSurvey(ng)) ! ! Set number of datum to process at current time-step. ! IF (backward) THEN NendObs(ng)=NstrObs(ng)-1 NstrObs(ng)=NstrObs(ng)-Nobs(ng) ELSE NstrObs(ng)=NendObs(ng)+1 NendObs(ng)=NstrObs(ng)+Nobs(ng)-1 END IF ! ! Set starting index of obervation vectors for reading. In weak ! constraint, the entire observation data is loaded. Otherwise, ! only the observartion for the current time window are loaded ! and started from vector index one. ! # ifdef WEAK_CONSTRAINT Mstr=NstrObs(ng) Mend=NendObs(ng) # else Mstr=1 Mend=Nobs(ng) # endif ! ! Read in observation type identifier. ! CALL netcdf_get_ivar (ng, model, OBS(ng)%name, Vname(1,idOtyp), & & ObsType(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN ! ! Read in observation provenance identifier. ! CALL netcdf_get_ivar (ng, model, OBS(ng)%name, Vname(1,idOpro), & & ObsProv(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN ! ! Read in observation time (days). ! CALL netcdf_get_time (ng, model, OBS(ng)%name, Vname(1,idObsT), & & Rclock%DateNumber, Tobs(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN ! ! Read in observation X-location (grid units). ! CALL netcdf_get_fvar (ng, model, OBS(ng)%name, Vname(1,idObsX), & & Xobs(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN ! ! Read in observation Y-location (grid units). ! CALL netcdf_get_fvar (ng, model, OBS(ng)%name, Vname(1,idObsY), & & Yobs(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN # ifdef SOLVE3D ! ! Read in observation depth, Zobs. If negative, depth is meter. If ! greater than zero, Zobs is in model fractional vertical levels ! (1 <= Zobs <= N). If Zobs < 0, its fractional level value is ! computed in routine "extract_obs3d" and over-written so it can ! be written into the observation NetCDF file for latter use. ! IF (wrote_Zobs(ng)) THEN CALL netcdf_get_fvar (ng, model, OBS(ng)%name, & & Vname(1,idObsZ), & & Zobs(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN ELSE CALL netcdf_get_fvar (ng, model, OBS(ng)%name, & & Vname(1,idObsD), & & Zobs(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN END IF Load_Zobs(ng)=.FALSE. IF ((MINVAL(Zobs).lt.0.0_r8).or. & & (MAXVAL(Zobs).lt.0.0_r8)) THEN Load_Zobs(ng)=.TRUE. END IF # ifdef DISTRIBUTE ! ! If distributed-memory and Zobs in meters (Zobs < 0), zero-out ! Zobs values in all nodes but itself to facilitate exchages between ! tiles latter before writting into observation NetCDF file. ! IF (.not.wrote_Zobs(ng)) THEN CALL obs_depth (ng, MyRank, model) END IF # endif # endif ! ! Read in observation values. ! CALL netcdf_get_fvar (ng, model, OBS(ng)%name, Vname(1,idOval), & & ObsVal(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN ! ! If appropriate, read in observation meta values. ! IF (haveObsMeta(ng)) THEN CALL netcdf_get_fvar (ng, model, OBS(ng)%name, & & Vname(1,idOmet), & & ObsMeta(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN END IF ! # ifdef WEAK_CONSTRAINT ! Read in observation error covariance. # else ! Read in observation error covariance. To avoid successive divisions, ! convert to inverse observation error covariance. # endif ! CALL netcdf_get_fvar (ng, model, OBS(ng)%name, Vname(1,idOerr), & & ObsErr(Mstr:), & & ncid = OBS(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN # ifndef WEAK_CONSTRAINT DO iobs=1,Nobs(ng) ObsErr(iobs)=1.0_r8/ObsErr(iobs) END DO # endif ! ! Read in nonlinear model values at observation locations. ! IF (readNLmod.and.haveNLmod(ng)) THEN # ifdef VERIFICATION CALL netcdf_get_fvar (ng, model, DAV(ng)%name, & & Vname(1,idNLmo), & & NLmodVal(Mstr:), & & ncid = DAV(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN # else CALL netcdf_get_fvar (ng, model, DAV(ng)%name, & & Vname(1,idNLmo), & & NLmodVal(Mstr:), & & ncid = DAV(ng)%ncid, & & start = (/NstrObs(ng),outer/), & & total = (/Nobs(ng),1/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN # endif # ifdef BGQC_NOT_NEEDED CALL netcdf_get_fvar (ng, model, DAV(ng)%name, & & Vname(1,idBgEr), & & BgErr(Mstr:), & & ncid = DAV(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN CALL netcdf_get_fvar (ng, model, DAV(ng)%name, & & Vname(1,idObsS), & & ObsScale(Mstr:), & & ncid = DAV(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN # endif END IF # if defined TLM_OBS && !defined WEAK_CONSTRAINT ! ! If adjoint pass and incremental 4DVar, read in tangent linear model ! values at observation locations. ! IF (readTLmod.and.haveTLmod(ng)) THEN CALL netcdf_get_fvar (ng, model, DAV(ng)%name, & & Vname(1,idTLmo), & & TLmodVal(Mstr:), & & ncid = DAV(ng)%ncid, & & start = (/NstrObs(ng)/), & & total = (/Nobs(ng)/)) IF (FoundError(exit_flag, NoError, __LINE__, & & __FILE__)) RETURN # if defined IS4DVAR ! ! Reset TLM values to zero at the first pass of the inner loop. ! IF (inner.eq.0) THEN DO iobs=1,Mobs TLmodVal(iobs)=0.0_r8 END DO END IF # endif END IF # endif ! !----------------------------------------------------------------------- ! Set counters for number of observations to processed for each state ! variable. !----------------------------------------------------------------------- ! DO iobs=Mstr,Mend IF (ObsType(iobs).eq.ObsState2Type(isFsur)) THEN FOURDVAR(ng)%ObsCount(isFsur)= & & FOURDVAR(ng)%ObsCount(isFsur)+1 ELSE IF (ObsType(iobs).eq.ObsState2Type(isUbar)) THEN FOURDVAR(ng)%ObsCount(isUbar)= & & FOURDVAR(ng)%ObsCount(isUbar)+1 ELSE IF (ObsType(iobs).eq.ObsState2Type(isVbar)) THEN FOURDVAR(ng)%ObsCount(isVbar)= & & FOURDVAR(ng)%ObsCount(isVbar)+1 # ifdef SOLVE3D ELSE IF (ObsType(iobs).eq.ObsState2Type(isUvel)) THEN FOURDVAR(ng)%ObsCount(isUvel)= & & FOURDVAR(ng)%ObsCount(isUvel)+1 ELSE IF (ObsType(iobs).eq.ObsState2Type(isVvel)) THEN FOURDVAR(ng)%ObsCount(isVvel)= & & FOURDVAR(ng)%ObsCount(isVvel)+1 ELSE IF (ObsType(iobs).eq.ObsState2Type(isRadial)) THEN haveObsMeta(ng)=.TRUE. FOURDVAR(ng)%ObsCount(isRadial)= & & FOURDVAR(ng)%ObsCount(isRadial)+1 ELSE DO itrc=1,NT(ng) IF (ObsType(iobs).eq.ObsState2Type(isTvar(itrc))) THEN i=isTvar(itrc) FOURDVAR(ng)%ObsCount(i)=FOURDVAR(ng)%ObsCount(i)+1 END IF END DO # endif END IF END DO ! !----------------------------------------------------------------------- ! If applicable, set next observation survey time to process. !----------------------------------------------------------------------- ! IF (Master) THEN CALL time_string (ObsTime(ng), t_code) WRITE (stdout,10) ObsTime(ng)*sec2day, t_code END IF IF (backward) THEN IF ((ObsSurvey(ng)-1).ge.1) THEN ObsTime(ng)=FOURDVAR(ng)%SurveyTime(ObsSurvey(ng)-1)*day2sec END IF ELSE IF ((ObsSurvey(ng)+1).le.Nsurvey(ng)) THEN ObsTime(ng)=FOURDVAR(ng)%SurveyTime(ObsSurvey(ng)+1)*day2sec END IF END IF END IF ! 10 FORMAT (/,' Number of State Observations Processed:',2x, & & 'ObsTime = ',f12.4,',',t68,a,/,/, & & 10x,'Variable',10x,'IstrObs',4x,'IendObs',6x,'Count', & & 3x,'Rejected',/) RETURN END SUBROUTINE obs_read #else SUBROUTINE obs_read RETURN END SUBROUTINE obs_read #endif