obs_initial.F

Go to the documentation of this file.

#include "cppdefs.h"
      MODULE obs_initial_mod
#if (defined FOUR_DVAR || defined VERIFICATION) && defined OBSERVATIONS
!
!git $Id$
!================================================== Hernan G. Arango ===
!  Copyright (c) 2002-2025 The ROMS Group                              !
!    Licensed under a MIT/X style license                              !
!    See License_ROMS.md                                               !
!=======================================================================
!                                                                      !
!  This module opens and inquires observations file using the standard !
!  NetCDF library or the Parallel-IO (PIO) library.                    !
!                                                                      !
!  It sets various variables needed for the processing of the state    !
!  solution at observations locations in variational data assimilation.!
!                                                                      !
!=======================================================================
!
      USE mod_param
      USE mod_parallel
      USE mod_fourdvar
      USE mod_iounits
      USE mod_ncparam
      USE mod_scalars
!
      USE strings_mod, ONLY : founderror
!
      implicit none
!
      PUBLIC  :: obs_initial
      PRIVATE :: obs_initial_nf90
# if defined PIO_LIB && defined DISTRIBUTE
      PRIVATE :: obs_initial_pio
# endif
!
      CONTAINS
!
!***********************************************************************
      SUBROUTINE obs_initial (ng, model, backward)
!***********************************************************************
!
!  Imported variable declarations.
!
      logical, intent(in) :: backward
!
      integer, intent(in) :: ng, model
!
!  Local variable declarations.
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__
!
!-----------------------------------------------------------------------
!  Write out time-averaged fields according to IO type.
!-----------------------------------------------------------------------
!
      SELECT CASE (obs(ng)%IOtype)
        CASE (io_nf90)
          CALL obs_initial_nf90 (ng, model, backward)
 
# if defined PIO_LIB && defined DISTRIBUTE
        CASE (io_pio)
          CALL obs_initial_pio (ng, model, backward)
# endif
        CASE DEFAULT
          IF (master) WRITE (stdout,10) obs(ng)%IOtype
          exit_flag=2
      END SELECT
      IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
  10  FORMAT (' OBS_INITIAL - Illegal input file type, io_type = ',i0,  &
     &        /,15x,'Check KeyWord ''INP_LIB'' in ''roms.in''.')
!
      RETURN
      END SUBROUTINE obs_initial
!
!***********************************************************************
      SUBROUTINE obs_initial_nf90 (ng, model, backward)
!***********************************************************************
!
      USE mod_netcdf
!
!  Imported variable declarations.
!
      logical, intent(in) :: backward
!
      integer, intent(in) :: ng, model
!
!  Local variable declarations.
!
      logical, dimension(NV) :: got_var(nv)
!
      integer :: ifirst, i, nvd, recdim, status
      integer :: vsize(4)
!
# ifdef WEAK_CONSTRAINT
      real(r8), parameter :: inival = 0.0_r8
# endif
      real(r8) :: tend
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//"obs_initial_nf90"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Inquire about the contents of observation NetCDF file:  Inquire about
!  the dimensions and variables.
!-----------------------------------------------------------------------
!
      query : IF (obs(ng)%ncid.eq.-1) THEN
!
!  Open observations NetCDF file.
!
        CALL netcdf_open (ng, model, obs(ng)%name, 1, obs(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          WRITE (stdout,10) trim(obs(ng)%name)
          RETURN
        END IF
!
!  Inquire about the variables.
!
        CALL netcdf_inq_var (ng, model, obs(ng)%name,                   &
     &                       obs(ng)%ncid)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize logical switches.
!
        DO i=1,nv
          got_var(i)=.false.
        END DO
!
!  Scan variable list from observation NetCDF and activate switches for
!  required variables.
!
        DO i=1,n_var
          IF (trim(var_name(i)).eq.trim(vname(1,idoday))) THEN
            got_var(idoday)=.true.
            obs(ng)%Vid(idoday)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idnobs))) THEN
            got_var(idnobs)=.true.
            obs(ng)%Vid(idnobs)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idotyp))) THEN
            got_var(idotyp)=.true.
            obs(ng)%Vid(idotyp)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobst))) THEN
            got_var(idobst)=.true.
            obs(ng)%Vid(idobst)=var_id(i)
# ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsd))) THEN
            got_var(idobsd)=.true.
            obs(ng)%Vid(idobsd)=var_id(i)
# endif
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsx))) THEN
            got_var(idobsx)=.true.
            obs(ng)%Vid(idobsx)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsy))) THEN
            got_var(idobsy)=.true.
            obs(ng)%Vid(idobsy)=var_id(i)
# ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsz))) THEN
            got_var(idobsz)=.true.
            obs(ng)%Vid(idobsz)=var_id(i)
# endif
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idoerr))) THEN
            got_var(idoerr)=.true.
            obs(ng)%Vid(idoerr)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idoval))) THEN
            got_var(idoval)=.true.
            obs(ng)%Vid(idoval)=var_id(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idomet))) THEN
            got_var(idomet)=.true.
            haveobsmeta(ng)=.true.
            obs(ng)%Vid(idomet)=var_id(i)
          END IF
        END DO
!
!  Check if needed obsrvation variables are available.
!
        IF (.not.got_var(idoday)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idoday)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idnobs)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idnobs)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idotyp)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idotyp)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idobst)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobst)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
# ifdef SOLVE3D
        IF (.not.got_var(idobsd)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsd)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
# endif
        IF (.not.got_var(idobsx)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsx)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idobsy)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsy)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
# ifdef SOLVE3D
        IF (.not.got_var(idobsz)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsz)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
# endif
        IF (.not.got_var(idoerr)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idoerr)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idoval)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idoval)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
      END IF query
!
!-----------------------------------------------------------------------
!  Set observation processing variables.
!-----------------------------------------------------------------------
!
!  Determine if there is any data available in the model time
!  window.  Set first survey record to process.
!
      IF (backward) THEN
        ifirst=0
# ifdef GENERIC_DSTART
        tend=(time(ng)-ntimes(ng)*dt(ng))*sec2day
# else
        tend=(time(ng)-(ntstart(ng)-1)*dt(ng))*sec2day
# endif
        DO i=1,nsurvey(ng)
          IF ((tend.le.fourdvar(ng)%SurveyTime(i)).and.                 &
     &        (fourdvar(ng)%SurveyTime(i).le.tdays(ng))) THEN
            ifirst=max(ifirst,i)
          END IF
        END DO
# ifndef SP4DVAR
        IF (ifirst.eq.0) THEN
          IF (master) WRITE (stdout,30) tend, tdays(ng)
          exit_flag=2
          RETURN
        END IF
# endif
      ELSE
        ifirst=nsurvey(ng)
        tend=(time(ng)+ntimes(ng)*dt(ng))*sec2day
        DO i=1,nsurvey(ng)
          IF ((tdays(ng).le.fourdvar(ng)%SurveyTime(i)).and.            &
     &        (fourdvar(ng)%SurveyTime(i).le.tend)) THEN
            ifirst=min(ifirst,i)
          END IF
        END DO
# ifndef SP4DVAR
        IF (ifirst.eq.0) THEN
          IF (master) WRITE (stdout,30) tdays(ng), tend
          exit_flag=2
          RETURN
        END IF
# endif
      END IF
      obstime(ng)=fourdvar(ng)%SurveyTime(ifirst)*day2sec
!
!  Initialize observation survey counter.  This is the counter of data
!  assimilation cycles or observations survey times on which the model
!  state is extracted (interpolated) at the observation locations.
!
      IF (backward) THEN
        obssurvey(ng)=ifirst+1
      ELSE
        obssurvey(ng)=ifirst-1
      END IF
!
!  Initialize time switch to process model state at observation
!  locations.
!
      processobs(ng)=.false.
 
# ifdef I4DVAR
!
!  Initialize cost function misfit between model and observations.
!  The IF statement is to avoid rewritting its value before it is
!  written into the initial NetCDF file.
!
      IF (.not.backward) THEN
        DO i=0,nobsvar(ng)
          fourdvar(ng)%ObsCost(i)=0.0_r8
        END DO
      END IF
# endif
!
!  Set staring and ending observation indices.
!
      IF (backward) THEN
        nstrobs(ng)=0
        nendobs(ng)=0
        DO i=1,ifirst
          nstrobs(ng)=nstrobs(ng)+fourdvar(ng)%NobsSurvey(i)
        END DO
        nstrobs(ng)=nstrobs(ng)+1
      ELSE
        IF (ifirst.eq.1) THEN
          nstrobs(ng)=0
          nendobs(ng)=0
        ELSE
          nstrobs(ng)=0
          nendobs(ng)=0
          DO i=1,ifirst-1
            nendobs(ng)=nendobs(ng)+fourdvar(ng)%NobsSurvey(i)
          END DO
        END IF
      END IF
!
!  Initialize total observation counter in structure array.  Notice
!  that the zero index carries the total summation.
!
      fourdvar(ng)%ObsCount(0)=0
      fourdvar(ng)%ObsReject(0)=0
 
# ifdef WEAK_CONSTRAINT
!
!  Initialize model values at observation locations.
!
      DO i=1,mobs
        nlmodval(i)=inival
        obsvetting(i)=inival
#  ifndef SP4DVAR
        tlmodval(i)=inival
#  endif
#  if defined SP4DVAR           && defined DISJOINTED && \
      defined concurrent_kernel
        tlmodval(i)=inival
#  endif
      END DO
#  if defined SP4DVAR           && defined DISJOINTED && \
     !defined CONCURRENT_KERNEL
      IF (model.ne.iadm) THEN
        DO i=1,mobs
          tlmodval(i)=inival
        END DO
      END IF
#  endif
# endif
!
  10  FORMAT (/,' OBS_INITIAL_NF90 - unable to open input NetCDF',      &
     &        ' file: ',a)
  20  FORMAT (/,' OBS_INITIAL_NF90 - unable to find model variable:',   &
     &        1x,a,/,20x,'in input NetCDF file: ',a)
  30  FORMAT (/,' OBS_INITIAL_NF90 - No are observations available',    &
     &        ' for time window (days): ',/,12x,f12.4,' - ',f12.4,/)
!
      RETURN
      END SUBROUTINE obs_initial_nf90
 
# if defined PIO_LIB && defined DISTRIBUTE
!
!***********************************************************************
      SUBROUTINE obs_initial_pio (ng, model, backward)
!***********************************************************************
!
      USE mod_pio_netcdf
!
!  Imported variable declarations.
!
      logical, intent(in) :: backward
!
      integer, intent(in) :: ng, model
!
!  Local variable declarations.
!
      logical, dimension(NV) :: got_var(nv)
!
      integer :: ifirst, i, nvd, recdim, status
      integer :: vsize(4)
!
#  ifdef WEAK_CONSTRAINT
      real(r8), parameter :: inival = 0.0_r8
#  endif
      real(r8) :: tend
!
      character (len=*), parameter :: myfile =                          &
     &  __FILE__//"obs_initial_nf90"
!
      sourcefile=myfile
!
!-----------------------------------------------------------------------
!  Inquire about the contents of observation NetCDF file:  Inquire about
!  the dimensions and variables.
!-----------------------------------------------------------------------
!
      query : IF (obs(ng)%ncid.eq.-1) THEN
!
!  Open observations NetCDF file.
!
        CALL pio_netcdf_open (ng, model, obs(ng)%name, 1,               &
     &                        obs(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) THEN
          WRITE (stdout,10) trim(obs(ng)%name)
          RETURN
        END IF
!
!  Inquire about the variables.
!
        CALL pio_netcdf_inq_var (ng, model, obs(ng)%name,               &
     &                           obs(ng)%pioFile)
        IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
!
!  Initialize logical switches.
!
        DO i=1,nv
          got_var(i)=.false.
        END DO
!
!  Scan variable list from observation NetCDF and activate switches for
!  required variables.
!
        DO i=1,n_var
          IF (trim(var_name(i)).eq.trim(vname(1,idoday))) THEN
            got_var(idoday)=.true.
            obs(ng)%pioVar(idoday)%vd=var_desc(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idnobs))) THEN
            got_var(idnobs)=.true.
            obs(ng)%pioVar(idnobs)%vd=var_desc(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idotyp))) THEN
            got_var(idotyp)=.true.
            obs(ng)%pioVar(idotyp)%vd=var_desc(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobst))) THEN
            got_var(idobst)=.true.
            obs(ng)%pioVar(idobst)%vd=var_desc(i)
#  ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsd))) THEN
            got_var(idobsd)=.true.
            obs(ng)%pioVar(idobsd)%vd=var_desc(i)
#  endif
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsx))) THEN
            got_var(idobsx)=.true.
            obs(ng)%pioVar(idobsx)%vd=var_desc(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsy))) THEN
            got_var(idobsy)=.true.
            obs(ng)%pioVar(idobsy)%vd=var_desc(i)
#  ifdef SOLVE3D
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idobsz))) THEN
            got_var(idobsz)=.true.
            obs(ng)%pioVar(idobsz)%vd=var_desc(i)
#  endif
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idoerr))) THEN
            got_var(idoerr)=.true.
            obs(ng)%pioVar(idoerr)%vd=var_desc(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idoval))) THEN
            got_var(idoval)=.true.
            obs(ng)%pioVar(idoval)%vd=var_desc(i)
          ELSE IF (trim(var_name(i)).eq.trim(vname(1,idomet))) THEN
            got_var(idomet)=.true.
            haveobsmeta(ng)=.true.
            obs(ng)%pioVar(idomet)%vd=var_desc(i)
          END IF
        END DO
!
!  Check if needed obsrvation variables are available.
!
        IF (.not.got_var(idoday)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idoday)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idnobs)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idnobs)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idotyp)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idotyp)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idobst)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobst)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
#  ifdef SOLVE3D
        IF (.not.got_var(idobsd)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsd)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
#  endif
        IF (.not.got_var(idobsx)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsx)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idobsy)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsy)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
#  ifdef SOLVE3D
        IF (.not.got_var(idobsz)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idobsz)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
#  endif
        IF (.not.got_var(idoerr)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idoerr)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
        IF (.not.got_var(idoval)) THEN
          IF (master) WRITE (stdout,20) trim(vname(1,idoval)),          &
     &                                  trim(obs(ng)%name)
          exit_flag=2
          RETURN
        END IF
      END IF query
!
!-----------------------------------------------------------------------
!  Set observation processing variables.
!-----------------------------------------------------------------------
!
!  Determine if there is any data available in the model time
!  window.  Set first survey record to process.
!
      IF (backward) THEN
        ifirst=0
#  ifdef GENERIC_DSTART
        tend=(time(ng)-ntimes(ng)*dt(ng))*sec2day
#  else
        tend=(time(ng)-(ntstart(ng)-1)*dt(ng))*sec2day
#  endif
        DO i=1,nsurvey(ng)
          IF ((tend.le.fourdvar(ng)%SurveyTime(i)).and.                 &
     &        (fourdvar(ng)%SurveyTime(i).le.tdays(ng))) THEN
            ifirst=max(ifirst,i)
          END IF
        END DO
#  ifndef SP4DVAR
        IF (ifirst.eq.0) THEN
          IF (master) WRITE (stdout,30) tend, tdays(ng)
          exit_flag=2
          RETURN
        END IF
#  endif
      ELSE
        ifirst=nsurvey(ng)
        tend=(time(ng)+ntimes(ng)*dt(ng))*sec2day
        DO i=1,nsurvey(ng)
          IF ((tdays(ng).le.fourdvar(ng)%SurveyTime(i)).and.            &
     &        (fourdvar(ng)%SurveyTime(i).le.tend)) THEN
            ifirst=min(ifirst,i)
          END IF
        END DO
#  ifndef SP4DVAR
        IF (ifirst.eq.0) THEN
          IF (master) WRITE (stdout,30) tdays(ng), tend
          exit_flag=2
          RETURN
        END IF
#  endif
      END IF
      obstime(ng)=fourdvar(ng)%SurveyTime(ifirst)*day2sec
!
!  Initialize observation survey counter.  This is the counter of data
!  assimilation cycles or observations survey times on which the model
!  state is extracted (interpolated) at the observation locations.
!
      IF (backward) THEN
        obssurvey(ng)=ifirst+1
      ELSE
        obssurvey(ng)=ifirst-1
      END IF
!
!  Initialize time switch to process model state at observation
!  locations.
!
      processobs(ng)=.false.
 
#  ifdef I4DVAR
!
!  Initialize cost function misfit between model and observations.
!  The IF statement is to avoid rewritting its value before it is
!  written into the initial NetCDF file.
!
      IF (.not.backward) THEN
        DO i=0,nobsvar(ng)
          fourdvar(ng)%ObsCost(i)=0.0_r8
        END DO
      END IF
#  endif
!
!  Set staring and ending observation indices.
!
      IF (backward) THEN
        nstrobs(ng)=0
        nendobs(ng)=0
        DO i=1,ifirst
          nstrobs(ng)=nstrobs(ng)+fourdvar(ng)%NobsSurvey(i)
        END DO
        nstrobs(ng)=nstrobs(ng)+1
      ELSE
        IF (ifirst.eq.1) THEN
          nstrobs(ng)=0
          nendobs(ng)=0
        ELSE
          nstrobs(ng)=0
          nendobs(ng)=0
          DO i=1,ifirst-1
            nendobs(ng)=nendobs(ng)+fourdvar(ng)%NobsSurvey(i)
          END DO
        END IF
      END IF
!
!  Initialize total observation counter in structure array.  Notice
!  that the zero index carries the total summation.
!
      fourdvar(ng)%ObsCount(0)=0
      fourdvar(ng)%ObsReject(0)=0
 
#  ifdef WEAK_CONSTRAINT
!
!  Initialize model values at observation locations.
!
      DO i=1,mobs
        nlmodval(i)=inival
        obsvetting(i)=inival
#   ifndef SP4DVAR
        tlmodval(i)=inival
#   endif
#   if defined SP4DVAR           && defined DISJOINTED && \
      defined concurrent_kernel
        tlmodval(i)=inival
#   endif
      END DO
#   if defined SP4DVAR           && defined DISJOINTED && \
     !defined CONCURRENT_KERNEL
      IF (model.ne.iadm) THEN
        DO i=1,mobs
          tlmodval(i)=inival
        END DO
      END IF
#   endif
#  endif
!
  10  FORMAT (/,' OBS_INITIAL_PIO - unable to open input NetCDF',       &
     &        ' file: ',a)
  20  FORMAT (/,' OBS_INITIAL_PIO - unable to find model variable:',    &
     &        1x,a,/,20x,'in input NetCDF file: ',a)
  30  FORMAT (/,' OBS_INITIAL_PIO - No are observations available',     &
     &        ' for time window (days): ',/,12x,f12.4,' - ',f12.4,/)
!
      RETURN
      END SUBROUTINE obs_initial_pio
# endif
#endif
      END MODULE obs_initial_mod