MODULE obs_write_mod
!
!git $Id$
!svn $Id: obs_write.F 1189 2023-08-15 21:26:58Z 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 interpolates nonlinear (background) and/or tangent !
! linear model (increments) state at observations location when !
! appropriate. !
! !
! It writes the model values at observation locations using the !
! standard NetCDF library or the Parallel-IO (PIO) library. !
! !
!=======================================================================
!
USE mod_param
USE mod_parallel
USE mod_fourdvar
USE mod_grid
USE mod_iounits
USE mod_ncparam
USE mod_ocean
USE mod_scalars
USE mod_stepping
!
USE distribute_mod, ONLY : mp_collect
USE extract_obs_mod, ONLY : extract_obs2d
USE extract_obs_mod, ONLY : extract_obs3d
USE strings_mod, ONLY : FoundError
!
implicit none
!
PUBLIC :: obs_write
PRIVATE :: obs_operator
PRIVATE :: obs_write_nf90
!
CONTAINS
!
!***********************************************************************
SUBROUTINE obs_write (ng, tile, model)
!***********************************************************************
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
!
! Local variable declarations.
!
integer :: Mstr, Mend
integer :: LBi, UBi, LBj, UBj
integer :: ObsSum, ObsVoid
integer :: i, ie, is, iobs, itrc
!
character (len=50) :: string
character (len=*), parameter :: MyFile = &
& "ROMS/Utility/obs_write.F"
!
!-----------------------------------------------------------------------
! Extract state vector at the observation locations in space and
! time.
!-----------------------------------------------------------------------
!
H_OPERATOR : IF (ProcessObs(ng)) THEN
LBi=BOUNDS(ng)%LBi(tile)
UBi=BOUNDS(ng)%UBi(tile)
LBj=BOUNDS(ng)%LBj(tile)
UBj=BOUNDS(ng)%UBj(tile)
!
CALL obs_operator (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& Mstr, Mend)
!
!-----------------------------------------------------------------------
! Write out state vector to output Data Assimilation Variables (DAV)
! NetCDF file.
!-----------------------------------------------------------------------
!
SELECT CASE (DAV(ng)%IOtype)
CASE (io_nf90)
CALL obs_write_nf90 (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& Mstr, Mend)
CASE DEFAULT
IF (Master) WRITE (stdout,10) DAV(ng)%IOtype
exit_flag=3
END SELECT
IF (FoundError(exit_flag, NoError, 110, MyFile)) RETURN
!
!-----------------------------------------------------------------------
! Set counters for the number of rejected observations for each state
! variable.
!-----------------------------------------------------------------------
!
DO iobs=Mstr,Mend
IF (ObsScale(iobs).lt.1.0) THEN
IF (ObsType(iobs).eq.ObsState2Type(isFsur)) THEN
FOURDVAR(ng)%ObsReject(isFsur)= &
& FOURDVAR(ng)%ObsReject(isFsur)+1
ELSE IF (ObsType(iobs).eq.ObsState2Type(isUbar)) THEN
FOURDVAR(ng)%ObsReject(isUbar)= &
& FOURDVAR(ng)%ObsReject(isUbar)+1
ELSE IF (ObsType(iobs).eq.ObsState2Type(isVbar)) THEN
FOURDVAR(ng)%ObsReject(isVbar)= &
& FOURDVAR(ng)%ObsReject(isVbar)+1
ELSE IF (ObsType(iobs).eq.ObsState2Type(isUvel)) THEN
FOURDVAR(ng)%ObsReject(isUvel)= &
& FOURDVAR(ng)%ObsReject(isUvel)+1
ELSE IF (ObsType(iobs).eq.ObsState2Type(isVvel)) THEN
FOURDVAR(ng)%ObsReject(isVvel)= &
& FOURDVAR(ng)%ObsReject(isVvel)+1
ELSE IF (ObsType(iobs).eq.ObsState2Type(isRadial)) THEN
FOURDVAR(ng)%ObsReject(isRadial)= &
& FOURDVAR(ng)%ObsReject(isRadial)+1
ELSE
DO itrc=1,NT(ng)
IF (ObsType(iobs).eq.ObsState2Type(isTvar(itrc))) THEN
i=isTvar(itrc)
FOURDVAR(ng)%ObsReject(i)=FOURDVAR(ng)%ObsReject(i)+1
END IF
END DO
END IF
END IF
END DO
!
! Load total available and rejected observations into structure
! array.
!
DO i=1,NobsVar(ng)
FOURDVAR(ng)%ObsCount(0)=FOURDVAR(ng)%ObsCount(0)+ &
& FOURDVAR(ng)%ObsCount(i)
FOURDVAR(ng)%ObsReject(0)=FOURDVAR(ng)%ObsReject(0)+ &
& FOURDVAR(ng)%ObsReject(i)
END DO
!
!-----------------------------------------------------------------------
! Report observation processing information.
!-----------------------------------------------------------------------
!
IF (Master) THEN
ObsSum=0
ObsVoid=0
is=NstrObs(ng)
DO i=1,NobsVar(ng)
IF (FOURDVAR(ng)%ObsCount(i).gt.0) THEN
ie=is+FOURDVAR(ng)%ObsCount(i)-1
WRITE (stdout,20) TRIM(ObsName(i)), is, ie, &
& ie-is+1, FOURDVAR(ng)%ObsReject(i)
is=ie+1
ObsSum=ObsSum+FOURDVAR(ng)%ObsCount(i)
ObsVoid=ObsVoid+FOURDVAR(ng)%ObsReject(i)
END IF
END DO
WRITE (stdout,30) ObsSum, ObsVoid, &
& FOURDVAR(ng)%ObsCount(0), &
& FOURDVAR(ng)%ObsReject(0)
END IF
!
IF (wrtNLmod(ng)) THEN
string='Wrote NLM state at observation locations, '
ELSE IF (wrtTLmod(ng)) THEN
string='Wrote TLM state at observation locations, '
ELSE IF (wrtRPmod(ng)) THEN
string='Wrote RPM state at observation locations, '
END IF
IF (Master) THEN
IF (wrtNLmod(ng).or.wrtTLmod(ng).or.wrtRPmod(ng)) THEN
WRITE (stdout,40) TRIM(string), NstrObs(ng), NendObs(ng)
END IF
END IF
END IF H_OPERATOR
!
10 FORMAT (' OBS_WRITE - Illegal output file type, io_type = ',i0, &
& /,13x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
20 FORMAT (10x,a,t25,4(1x,i10))
30 FORMAT (/,10x,'Total',t47,2(1x,i10), &
& /,10x,'Obs Tally',t47,2(1x,i10),/)
40 FORMAT (1x,a,' datum = ',i0,' - ',i0,/)
!
RETURN
END SUBROUTINE obs_write
!
!***********************************************************************
SUBROUTINE obs_operator (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& Mstr, Mend)
!***********************************************************************
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(out) :: Mstr, Mend
!
! Local variable declarations.
!
integer :: i, ic, iobs, itrc
integer :: j, k
integer :: Ncollect
!
real(r8), parameter :: IniVal = 0.0_r8
real(r8) :: angle
real(r8) :: misfit(Mobs), unvetted(Mobs)
real(r8) :: uradial(Mobs), vradial(Mobs)
real(r8) :: uBgErr(Mobs), vBgErr(Mobs)
!
character (len=*), parameter :: MyFile = &
& "ROMS/Utility/obs_write.F"//", obs_operator"
!
!-----------------------------------------------------------------------
! Set lower and upper tile bounds and staggered variables bounds for
! this horizontal domain partition. Notice that if tile=-1, it will
! set the values for the global grid.
!-----------------------------------------------------------------------
!
integer :: Istr, IstrB, IstrP, IstrR, IstrT, IstrM, IstrU
integer :: Iend, IendB, IendP, IendR, IendT
integer :: Jstr, JstrB, JstrP, JstrR, JstrT, JstrM, JstrV
integer :: Jend, JendB, JendP, JendR, JendT
integer :: Istrm3, Istrm2, Istrm1, IstrUm2, IstrUm1
integer :: Iendp1, Iendp2, Iendp2i, Iendp3
integer :: Jstrm3, Jstrm2, Jstrm1, JstrVm2, JstrVm1
integer :: Jendp1, Jendp2, Jendp2i, Jendp3
!
Istr =BOUNDS(ng) % Istr (tile)
IstrB =BOUNDS(ng) % IstrB (tile)
IstrM =BOUNDS(ng) % IstrM (tile)
IstrP =BOUNDS(ng) % IstrP (tile)
IstrR =BOUNDS(ng) % IstrR (tile)
IstrT =BOUNDS(ng) % IstrT (tile)
IstrU =BOUNDS(ng) % IstrU (tile)
Iend =BOUNDS(ng) % Iend (tile)
IendB =BOUNDS(ng) % IendB (tile)
IendP =BOUNDS(ng) % IendP (tile)
IendR =BOUNDS(ng) % IendR (tile)
IendT =BOUNDS(ng) % IendT (tile)
Jstr =BOUNDS(ng) % Jstr (tile)
JstrB =BOUNDS(ng) % JstrB (tile)
JstrM =BOUNDS(ng) % JstrM (tile)
JstrP =BOUNDS(ng) % JstrP (tile)
JstrR =BOUNDS(ng) % JstrR (tile)
JstrT =BOUNDS(ng) % JstrT (tile)
JstrV =BOUNDS(ng) % JstrV (tile)
Jend =BOUNDS(ng) % Jend (tile)
JendB =BOUNDS(ng) % JendB (tile)
JendP =BOUNDS(ng) % JendP (tile)
JendR =BOUNDS(ng) % JendR (tile)
JendT =BOUNDS(ng) % JendT (tile)
!
Istrm3 =BOUNDS(ng) % Istrm3 (tile) ! Istr-3
Istrm2 =BOUNDS(ng) % Istrm2 (tile) ! Istr-2
Istrm1 =BOUNDS(ng) % Istrm1 (tile) ! Istr-1
IstrUm2=BOUNDS(ng) % IstrUm2(tile) ! IstrU-2
IstrUm1=BOUNDS(ng) % IstrUm1(tile) ! IstrU-1
Iendp1 =BOUNDS(ng) % Iendp1 (tile) ! Iend+1
Iendp2 =BOUNDS(ng) % Iendp2 (tile) ! Iend+2
Iendp2i=BOUNDS(ng) % Iendp2i(tile) ! Iend+2 interior
Iendp3 =BOUNDS(ng) % Iendp3 (tile) ! Iend+3
Jstrm3 =BOUNDS(ng) % Jstrm3 (tile) ! Jstr-3
Jstrm2 =BOUNDS(ng) % Jstrm2 (tile) ! Jstr-2
Jstrm1 =BOUNDS(ng) % Jstrm1 (tile) ! Jstr-1
JstrVm2=BOUNDS(ng) % JstrVm2(tile) ! JstrV-2
JstrVm1=BOUNDS(ng) % JstrVm1(tile) ! JstrV-1
Jendp1 =BOUNDS(ng) % Jendp1 (tile) ! Jend+1
Jendp2 =BOUNDS(ng) % Jendp2 (tile) ! Jend+2
Jendp2i=BOUNDS(ng) % Jendp2i(tile) ! Jend+2 interior
Jendp3 =BOUNDS(ng) % Jendp3 (tile) ! Jend+3
!
SourceFile=MyFile
!
!=======================================================================
! Interpolate model state at observation locations.
!=======================================================================
!
!
! Set starting and ending indices of observations to process for the
! current time survey. In the dual formulation, the entire observation
! vector for the assimilation window is maintained and the data array
! indexes are global.
!
Mstr=NstrObs(ng)
Mend=NendObs(ng)
!
! Initialize observation reject/accept processing flag.
!
DO iobs=Mstr,Mend
unvetted(iobs)=IniVal
ObsVetting(iobs)=IniVal
END DO
!
! Some entries are not computed in the extraction routine. Set values
! to zero to avoid problems when writing non initialized values.
! Notice that only the appropriate indices are zero-out to facilate
! collecting all the extrated data as sum between all nodes.
!
IF (wrtNLmod(ng)) THEN
DO iobs=Mstr,Mend
NLmodVal(iobs)=IniVal
BgErr(iobs)=IniVal
uBgErr(iobs)=IniVal
vBgErr(iobs)=IniVal
END DO
END IF
IF (wrtTLmod(ng).or.wrtRPmod(ng)) THEN
DO iobs=Mstr,Mend
TLmodVal(iobs)=IniVal
END DO
END IF
!
!-----------------------------------------------------------------------
! Free-surface operator.
!-----------------------------------------------------------------------
!
! Nonlinear free-surface at observation locations.
!
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isFsur).gt.0)) THEN
CALL extract_obs2d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isFsur), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%zeta(:,:,kstp(ng)), &
& GRID(ng)%rmask, &
& NLmodVal)
!
! Free-surface backgound error standard deviation at observation
! locations.
!
CALL extract_obs2d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isFsur), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%e_zeta(:,:,1), &
& GRID(ng)%rmask, &
& BgErr)
END IF
!
! Tangent linear free-surface at observation locations.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isFsur).gt.0)) THEN
CALL extract_obs2d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isFsur), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%tl_zeta(:,:,kstp(ng)), &
& GRID(ng)%rmask, &
& TLmodVal)
END IF
!
!-----------------------------------------------------------------------
! Vertically integrated u-velocity operator.
!-----------------------------------------------------------------------
!
! Nonlinear 2D u-velocity at observation locations.
!
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isUbar).gt.0)) THEN
CALL extract_obs2d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isUbar), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%ubar(:,:,kstp(ng)), &
& GRID(ng)%umask, &
& NLmodVal)
!
! 2D u-velocity backgound error standard deviation at observation
! locations.
!
CALL extract_obs2d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isUbar), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%e_ubar(:,:,1), &
& GRID(ng)%umask, &
& BgErr)
END IF
!
! Tangent linear 2D u-velocity at observation locations.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isUbar).gt.0)) THEN
CALL extract_obs2d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isUbar), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%tl_ubar(:,:,kstp(ng)), &
& GRID(ng)%umask, &
& TLmodVal)
END IF
!
!-----------------------------------------------------------------------
! Vertically integrated v-velocity observations.
!-----------------------------------------------------------------------
!
! Nonlinear 2D v-velocity at observation locations.
!
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isVbar).gt.0)) THEN
CALL extract_obs2d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isVbar), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%vbar(:,:,kstp(ng)), &
& GRID(ng)%vmask, &
& NLmodVal)
!
! 2D v-velocity backgound error standard deviation at observation
! locations.
!
CALL extract_obs2d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isVbar), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%e_vbar(:,:,1), &
& GRID(ng)%vmask, &
& BgErr)
END IF
!
! Tangent linear 2D v-velocity at observation locations.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isVbar).gt.0)) THEN
CALL extract_obs2d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isVbar), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& OCEAN(ng)%tl_vbar(:,:,kstp(ng)), &
& GRID(ng)%vmask, &
& TLmodVal)
END IF
!
!-----------------------------------------------------------------------
! 3D u-velocity observations.
!-----------------------------------------------------------------------
!
! Nonlinear 3D u-velocity at observation locations.
!
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isUvel).gt.0)) THEN
DO k=1,N(ng)
DO j=Jstr-1,Jend+1
DO i=IstrU-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i-1,j,k)+ &
& GRID(ng)%z_r(i ,j,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isUvel), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%u(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%umask, &
& NLmodVal)
!
! 3D u-velocity backgound error standard deviation at observation
! locations.
!
CALL extract_obs3d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isUvel), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%e_u(:,:,:,1), &
& GRID(ng)%z_v, &
& GRID(ng)%umask, &
& BgErr)
END IF
!
! Tangent linear 3D u-velocity at observation locations.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isUvel).gt.0)) THEN
DO k=1,N(ng)
DO j=Jstr-1,Jend+1
DO i=IstrU-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i-1,j,k)+ &
& GRID(ng)%z_r(i ,j,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isUvel), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%tl_u(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%umask, &
& TLmodVal)
END IF
!
!-----------------------------------------------------------------------
! 3D v-velocity observations.
!-----------------------------------------------------------------------
!
! Nonlinear 3D v-velocity at observation locations.
!
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isVvel).gt.0)) THEN
DO k=1,N(ng)
DO j=JstrV-1,Jend+1
DO i=Istr-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i,j-1,k)+ &
& GRID(ng)%z_r(i,j ,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isVvel), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%v(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%vmask, &
& NLmodVal)
!
! 3D v-velocity backgound error standard deviation at observation
! locations.
!
CALL extract_obs3d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isVvel), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%e_v(:,:,:,1), &
& GRID(ng)%z_v, &
& GRID(ng)%vmask, &
& BgErr)
END IF
!
! Tangent linear 3D v-velocity at observation locations.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isVvel).gt.0)) THEN
DO k=1,N(ng)
DO j=JstrV-1,Jend+1
DO i=Istr-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i,j-1,k)+ &
& GRID(ng)%z_r(i,j ,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isVvel), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%tl_v(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%vmask, &
& TLmodVal)
END IF
!
!-----------------------------------------------------------------------
! Radial Velocity. The observations are in terms of radial speed and
! angle (stored in obs_meta). The observation angle converts the
! velocity components to geographical EAST and North components.
! By default, the radial velocity observations are processed as
! magnitude and heading angle (obs_meta; radians) in the navigation
! convention: an azimuth that is clockwise from TRUE North.
!
! In curvilinear coordinates, the radial forward problem is:
!
! radial = u * SIN(obs_meta + angler) + v * COS(obs_meta + angler)
!-----------------------------------------------------------------------
!
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isRadial).gt.0)) THEN
DO iobs=Mstr,Mend
uradial(iobs)=IniVal
vradial(iobs)=IniVal
END DO
!
! If curvilinear coordinates, interpolate grid rotation angle (radians)
! for each observation location.
!
CALL extract_obs2d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& GRID(ng)%angler, &
& GRID(ng)%rmask, &
& ObsAngler)
!
! Radial nonlinear u-component at observation locations.
!
DO k=1,N(ng)
DO j=Jstr-1,Jend+1
DO i=IstrU-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i-1,j,k)+ &
& GRID(ng)%z_r(i ,j,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%u(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%umask, &
& uradial)
!
! Radial u-component backgound error standard deviation at observation
! locations.
!
CALL extract_obs3d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%e_u(:,:,:,1), &
& GRID(ng)%z_v, &
& GRID(ng)%umask, &
& uBgErr)
!
! Radial nonlinear v-component at observation locations.
!
DO k=1,N(ng)
DO j=JstrV-1,Jend+1
DO i=Istr-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i,j-1,k)+ &
& GRID(ng)%z_r(i,j ,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%v(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%vmask, &
& vradial)
!
! Radial v-component backgound error standard deviation at observation
! locations.
!
CALL extract_obs3d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%e_v(:,:,:,1), &
& GRID(ng)%z_v, &
& GRID(ng)%vmask, &
& vBgErr)
!
! Compute nonlinear radial velocity.
!
DO iobs=Mstr,Mend
IF (ObsType(iobs).eq.ObsState2Type(isRadial)) THEN
angle=ObsMeta(iobs)+ObsAngler(iobs)
NLmodVal(iobs)=uradial(iobs)*SIN(angle)+ &
& vradial(iobs)*COS(angle)
BgErr(iobs)=MAX(uBgErr(iobs), vBgErr(iobs))
END IF
END DO
END IF
!
! Tangent linear radial velocities.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isRadial).gt.0)) THEN
DO iobs=Mstr,Mend
uradial(iobs)=IniVal
vradial(iobs)=IniVal
END DO
!
! If curvilinear coordinates, interpolate grid rotation angle (radians)
! for each observation location.
!
CALL extract_obs2d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, &
& GRID(ng)%angler, &
& GRID(ng)%rmask, &
& ObsAngler)
!
! Tangent linear radial u-component at the observation locations.
!
DO k=1,N(ng)
DO j=Jstr-1,Jend+1
DO i=IstrU-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i-1,j,k)+ &
& GRID(ng)%z_r(i ,j,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 1, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& uXmin(ng), uXmax(ng), &
& uYmin(ng), uYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%tl_u(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%umask, &
& uradial)
!
! Tangent linear radial v-component at the observation locations.
!
DO k=1,N(ng)
DO j=JstrV-1,Jend+1
DO i=Istr-1,Iend+1
GRID(ng)%z_v(i,j,k)=0.5_r8*(GRID(ng)%z_r(i,j-1,k)+ &
& GRID(ng)%z_r(i,j ,k))
END DO
END DO
END DO
CALL extract_obs3d (ng, 0, Lm(ng)+1, 1, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isRadial), &
& Mobs, Mstr, Mend, &
& vXmin(ng), vXmax(ng), &
& vYmin(ng), vYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%tl_v(:,:,:,nrhs(ng)), &
& GRID(ng)%z_v, &
& GRID(ng)%vmask, &
& vradial)
!
! Compute tangent linear radial velocity.
!
DO iobs=Mstr,Mend
IF (ObsType(iobs).eq.ObsState2Type(isRadial)) THEN
angle=ObsMeta(iobs)+ObsAngler(iobs)
TLmodVal(iobs)=uradial(iobs)*SIN(angle)+ &
& vradial(iobs)*COS(angle)
END IF
END DO
END IF
!
!-----------------------------------------------------------------------
! Tracers variable observations.
!-----------------------------------------------------------------------
!
DO itrc=1,NT(ng)
!
! Nonlinear tracer at the observation locations.
IF (wrtNLmod(ng).and. &
& (FOURDVAR(ng)%ObsCount(isTvar(itrc)).gt.0)) THEN
CALL extract_obs3d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isTvar(itrc)), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%t(:,:,:,nrhs(ng),itrc), &
& GRID(ng)%z_r, &
& GRID(ng)%rmask, &
& NLmodVal)
!
! Tracer backgound error standard deviation at observation locations.
!
CALL extract_obs3d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isTvar(itrc)), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%e_t(:,:,:,1,itrc), &
& GRID(ng)%z_r, &
& GRID(ng)%rmask, &
& BgErr)
END IF
!
! Tangent linear tracer at the observation locations.
!
IF ((wrtTLmod(ng).or.(wrtRPmod(ng))).and. &
& (FOURDVAR(ng)%ObsCount(isTvar(itrc)).gt.0)) THEN
CALL extract_obs3d (ng, 0, Lm(ng)+1, 0, Mm(ng)+1, &
& LBi, UBi, LBj, UBj, 1, N(ng), &
& ObsState2Type(isTvar(itrc)), &
& Mobs, Mstr, Mend, &
& rXmin(ng), rXmax(ng), &
& rYmin(ng), rYmax(ng), &
& time(ng), dt(ng), &
& ObsType, ObsVetting, &
& Tobs, Xobs, Yobs, Zobs, &
& OCEAN(ng)%tl_t(:,:,:,nrhs(ng),itrc), &
& GRID(ng)%z_r, &
& GRID(ng)%rmask, &
& TLmodVal)
END IF
END DO
!
!-----------------------------------------------------------------------
! If processing nonlinear trajectory, load observation reject/accept
! flag into screening variable. This needs to be done once and it is
! controlled by the main driver. Notice that in routine "extract_obs"
! the observations are only accepted if they are at the appropriate
! time window and grid bounded at water points.
!-----------------------------------------------------------------------
!
IF (wrtObsScale(ng).and.wrtNLmod(ng)) THEN
DO iobs=Mstr,Mend
ObsScale(iobs)=ObsVetting(iobs)
END DO
END IF
!
!-----------------------------------------------------------------------
! Collect screening variable for all extracted data.
!-----------------------------------------------------------------------
!
Ncollect=Mend-Mstr+1
IF (wrtObsScale(ng).and.wrtNLmod(ng)) THEN
CALL mp_collect (ng, model, Ncollect, IniVal, &
& ObsScale(Mstr:))
END IF
!
!-----------------------------------------------------------------------
! Scale extracted data at observations location. Save NLM unvetted
! values before they are rejected by quality control. The unvetted
! vector is used only for output purposes to analyze why a datum was
! rejected.
!-----------------------------------------------------------------------
!
! Nonlinear state vector.
!
IF (wrtNLmod(ng)) THEN
DO iobs=Mstr,Mend
unvetted(iobs)=NLmodVal(iobs)
NLmodVal(iobs)=ObsScale(iobs)*NLmodVal(iobs)
END DO
END IF
!
! Tangent linear state vector.
!
IF (wrtTLmod(ng).or.wrtRPmod(ng)) THEN
DO iobs=Mstr,Mend
TLmodVal(iobs)=ObsScale(iobs)*TLmodVal(iobs)
END DO
END IF
!
!-----------------------------------------------------------------------
! Collect extracted data.
!-----------------------------------------------------------------------
!
!
! Nonlinear state vector.
!
IF (wrtNLmod(ng)) THEN
CALL mp_collect (ng, model, Ncollect, IniVal, &
& NLmodVal(Mstr:))
!
CALL mp_collect (ng, model, Ncollect, IniVal, &
& unvetted(Mstr:))
END IF
!
! Background error standard deviation state vector.
!
IF (wrtObsScale(ng).and.wrtNLmod(ng)) THEN
CALL mp_collect (ng, model, Ncollect, IniVal, &
& BgErr(Mstr:))
END IF
!
! Tangent linear state vector.
!
IF (wrtTLmod(ng).or.wrtRPmod(ng)) THEN
CALL mp_collect (ng, model, Ncollect, IniVal, &
& TLmodVal(Mstr:))
END IF
!
! State vector depths.
!
IF (Load_Zobs(ng)) THEN
CALL mp_collect (ng, model, Ncollect, IniVal, &
& Zobs(Mstr:))
END IF
RETURN
END SUBROUTINE obs_operator
!
!***********************************************************************
SUBROUTINE obs_write_nf90 (ng, tile, model, &
& LBi, UBi, LBj, UBj, &
& Mstr, Mend)
!***********************************************************************
!
USE mod_netcdf
!
! Imported variable declarations.
!
integer, intent(in) :: ng, tile, model
integer, intent(in) :: LBi, UBi, LBj, UBj
integer, intent(in) :: Mstr, Mend
!
! Local variable declarations.
!
logical :: Lwrote
!
integer :: Tindex, iobs, status
!
real(dp) :: scale
!
character (len=*), parameter :: MyFile = &
& "ROMS/Utility/obs_write.F"//", obs_write_nf90"
!
SourceFile=MyFile
!
!-----------------------------------------------------------------------
! Compute and write initial and final model-observation misfit
! (innovation) vector for output purposes only. Write also initial
! nonlinear model at observation locations.
!-----------------------------------------------------------------------
!
IF (wrtMisfit(ng)) THEN
DO iobs=Mstr,Mend
misfit(iobs)=ObsScale(iobs)/SQRT(ObsErr(iobs))* &
& (NLmodVal(iobs)-ObsVal(iobs))
END DO
!
! Write out initial or final values into DAV NetCDF file.
!
IF (Nrun.eq.1) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idNLmp), &
& unvetted(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idNLmp))
IF (FoundError(exit_flag, NoError, 1369, MyFile)) RETURN
!
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idNLmi), &
& NLmodVal(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idNLmi))
IF (FoundError(exit_flag, NoError, 1377, MyFile)) RETURN
!
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idMOMi), &
& misfit(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idMOMi))
IF (FoundError(exit_flag, NoError, 1385, MyFile)) RETURN
ELSE
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idMOMf), &
& misfit(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idMOMf))
IF (FoundError(exit_flag, NoError, 1393, MyFile)) RETURN
END IF
END IF
!
!-----------------------------------------------------------------------
! Write out variables needed to compute Desroziers et al. (2005)
! statistics.
!-----------------------------------------------------------------------
!
IF (wrtNLmod(ng)) THEN
!
! 4D-Var background error standard deviation at observation locations.
!
IF (wrtObsScale(ng)) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idBgEr), &
& BgErr(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idBgEr))
IF (FoundError(exit_flag, NoError, 1431, MyFile)) RETURN
END IF
!
! 4D-Var innovation vector: observation minus background. Save initial
! background in the NLincrement vector.
!
IF (Nrun.eq.1) THEN
DO iobs=Mstr,Mend
innovation(iobs)=ObsScale(iobs)* &
& (ObsVal(iobs)-NLmodVal(iobs))
NLincrement(iobs)=NLmodVal(iobs)
END DO
!
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idInno), &
& innovation(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idInno))
IF (FoundError(exit_flag, NoError, 1450, MyFile)) RETURN
END IF
!
! 4D-Var increment (analysis minus background) and residual
! (observation minus analysis).
!
IF (outer.eq.Nouter) THEN
DO iobs=Mstr,Mend
NLincrement(iobs)=ObsScale(iobs)* &
& (NLmodVal(iobs)-NLincrement(iobs))
residual(iobs)=ObsScale(iobs)* &
& (ObsVal(iobs)-NLmodVal(iobs))
END DO
!
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idIncr), &
& NLincrement(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idIncr))
IF (FoundError(exit_flag, NoError, 1470, MyFile)) RETURN
!
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idResi), &
& residual(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idResi))
IF (FoundError(exit_flag, NoError, 1478, MyFile)) RETURN
END IF
END IF
!
!-----------------------------------------------------------------------
! Write out other variable into output DAC NetCDF file.
!-----------------------------------------------------------------------
!
! Current outer and inner loop.
!
IF (wrtNLmod(ng).or.wrtTLmod(ng).or.wrtRPmod(ng)) THEN
CALL netcdf_put_ivar (ng, model, DAV(ng)%name, 'outer', &
& outer, (/0/), (/0/), &
& ncid = DAV(ng)%ncid)
IF (FoundError(exit_flag, NoError, 1496, MyFile)) RETURN
!
CALL netcdf_put_ivar (ng, model, DAV(ng)%name, 'inner', &
& inner, (/0/), (/0/), &
& ncid = DAV(ng)%ncid)
IF (FoundError(exit_flag, NoError, 1501, MyFile)) RETURN
END IF
!
! Observation screening flag: accept or reject.
!
IF (wrtObsScale(ng).and.wrtNLmod(ng)) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idObsS), &
& ObsScale(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idObsS))
IF (FoundError(exit_flag, NoError, 1514, MyFile)) RETURN
END IF
!
! Nonlinear model or first guess (background) state at observation
! locations. If 4D-Var, the values are written for every outer loop.
! Also, the final value is written as individual variable to
! facilitate ERDDAP processing. Notice that it is always overwritten
! to avoid complicated logic with all the 4D-Var drivers. The last
! one written correspont to the final values.
!
IF (wrtNLmod(ng).and.outer.gt.0) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idNLmo), &
& NLmodVal(Mstr:), &
& (/NstrObs(ng),outer/), (/Nobs(ng),1/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idNLmo))
IF (FoundError(exit_flag, NoError, 1545, MyFile)) RETURN
END IF
!
IF (wrtNLmod(ng).and.outer.gt.0) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idNLmf), &
& NLmodVal(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idNLmf))
IF (FoundError(exit_flag, NoError, 1557, MyFile)) RETURN
haveNLmod(ng)=.TRUE.
END IF
!
! Write out unvetted nonlinear model at the observation locations per
! outer loop.
!
IF (wrtNLmod(ng).and.outer.gt.0) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idNLmu), &
& unvetted(Mstr:), &
& (/NstrObs(ng),outer/), (/Nobs(ng),1/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idNLmu))
IF (FoundError(exit_flag, NoError, 1575, MyFile)) RETURN
END IF
!
! Tangent linear model state increments at observation locations.
!
IF (wrtTLmod(ng).or.wrtRPmod(ng)) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idTLmo), &
& TLmodVal(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idTLmo))
IF (FoundError(exit_flag, NoError, 1590, MyFile)) RETURN
haveTLmod(ng)=.TRUE.
END IF
!
! Now, that Zobs has all the Z-locations in grid coordinates, write
! "obs_Zgrid" into data assimilation output file.
!
IF ((Nrun.eq.1).and. &
& (wrtNLmod(ng).or.wrtTLmod(ng).or.wrtRPmod(ng))) THEN
CALL netcdf_put_fvar (ng, model, DAV(ng)%name, &
& Vname(1,idObsZ), &
& Zobs(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = DAV(ng)%ncid, &
& varid = DAV(ng)%Vid(idObsZ))
IF (FoundError(exit_flag, NoError, 1682, MyFile)) RETURN
END IF
!
! Write Z-location of observation in grid coordinates, if applicable.
! This values are needed elsewhere when using the interpolation
! weight matrix. Recall that the depth of observations can be in
! meters or grid coordinates. Notice that since the model levels
! evolve in time, the fractional level coordinate is unknow during
! the processing of the observations.
!
IF ((Nrun.eq.1).and. &
& (wrtNLmod(ng).or.wrtTLmod(ng).or.wrtRPmod(ng))) THEN
CALL netcdf_put_fvar (ng, model, OBS(ng)%name, &
& Vname(1,idObsZ), &
& Zobs(Mstr:), &
& (/NstrObs(ng)/), (/Nobs(ng)/), &
& ncid = OBS(ng)%ncid, &
& varid = OBS(ng)%Vid(idObsZ))
IF (FoundError(exit_flag, NoError, 1700, MyFile)) RETURN
IF (model.eq.iADM) THEN
Lwrote=ObsSurvey(ng).eq.1
ELSE
Lwrote=ObsSurvey(ng).eq.Nsurvey(ng)
END IF
IF (Lwrote) wrote_Zobs(ng)=.TRUE.
END IF
!
!-----------------------------------------------------------------------
! Synchronize observations NetCDF file to disk.
!-----------------------------------------------------------------------
!
IF (wrtNLmod(ng).or.wrtTLmod(ng).or.wrtRPmod(ng)) THEN
CALL netcdf_sync (ng, model, DAV(ng)%name, DAV(ng)%ncid)
IF (FoundError(exit_flag, NoError, 1737, MyFile)) RETURN
IF (Nrun.eq.1) THEN
CALL netcdf_sync (ng, model, OBS(ng)%name, OBS(ng)%ncid)
IF (FoundError(exit_flag, NoError, 1742, MyFile)) RETURN
END IF
END IF
!
RETURN
END SUBROUTINE obs_write_nf90
END MODULE obs_write_mod