MODULE mod_fourdvar
!
!git $Id$
!svn $Id: mod_fourdvar.F 1202 2023-10-24 15:36:07Z arango $
!================================================== Hernan G. Arango ===
! Copyright (c) 2002-2023 The ROMS/TOMS Group !
! Licensed under a MIT/X style license !
! See License_ROMS.md !
!=======================================================================
! !
! Variational data assimilation variables: !
! !
! ADmodVal Adjoint model values at observation locations. !
! BackCost Current background cost function misfit (mean squared !
! difference) between model and background state, Jb. !
! CostFun Current iteration total cost function (background !
! plus observations), J. !
! CostFunOld Previous iteration total cost function (background !
! plus observations), J. !
! CostNorm Total cost function normalization scales !
! (minimization starting value, Jb+Jo). !
! Cost0 The total cost function at the beggining of the inner !
! loop (inner=0) for each outer loop. !
! DTsizeH Horizontal diffusion time-step size for spatial !
! convolutions. !
! DTsizeV Vertical diffusion time-step size for spatial !
! convolutions. !
! GradErr Upper bound on relatice error of the gradient. !
! HevecErr Maximum error bound on Hessian eigenvectors. !
! KhMax Maximum horizontal diffusion coefficient. !
! KhMin Minimum horizontal diffusion coefficient. !
! KvMax Maximum vertical diffusion coefficient. !
! KvMin Minimum vertical diffusion coefficient. !
! Load_Zobs Logical switch indicating that input Zobs is negative !
! and fractional vertical positions are computed in !
! "extract_obs3d" and written to observation NetCDF !
! file for latter use. !
! LhessianEV Logical switch to compute Hessian eigenvectors. !
! LhotStart Logical switch to activate hot start of subsquent !
! outer loops in the weak-constraint algorithms. !
! Lprecond Logical switch to activate conjugate gradient !
! preconditioning. !
! Lritz Logical switch to activw Ritz limited-memory !
! preconditioning. !
! LaugWeak Logical switch for computing augmented model error !
! terms in the routine "error_covariance". !
! nConvRitz Number of converged Ritz eigenvalues. !
! Nimpact Outer loop to consider for observations impact and !
! observations sensitivity. !
! NLmodVal Nonlinear model values at observation locations. !
! NHsteps Full number of horizontal diffusion steps for spatial !
! convolutions. !
! NLobsCost Current NLM observation cost function misfit (mean !
! squared difference) between NLM and observations. !
! NpostI Number of Lanczos iterations used for the posterior !
! analysis error covariance matrix estimation. !
! NritzEV If preconditioning, number of eigenpairs to use. !
! NVsteps Full number of vertical diffusion steps for spatial !
! convolutions. !
! ObsAngler Observation vector rotation angle at RHO-points. !
! ObsCost Current observation cost function misfit (mean !
! squared difference) between model and observations. !
! ObsCount Current observation counts per state variable. !
! ObsErr Observation error. !
! ObsMeta Observation meta values. !
! ObsName String describing the observation types. !
! ObsProv Observation provenance flags. !
! ObsReject Current rejected observation counts per state !
! variable. !
! ObsScale Observation screening and quality control scale, !
! 0: reject 1: accept. !
! ObsState2Type Mapping indices from state variable to observation !
! type. !
! ObsType2State Mapping indices from observation type to state !
! variable. !
! ObsType Observation type identifier. !
! ObsVal Observation values. !
! ObsVetting Processing flag used to reject (zero) or accept !
! (unity) observations. !
! Optimality normalized, optimal cost function minimum. !
! Ritz Ritz eigenvalues to compute approximated Hessian. !
! RitzMaxErr Ritz values maximum error limit. !
! TLmodVal Tangent linear model values at observation locations. !
! Vdecay Covariance vertical decorrelation scale (m). !
! Tobs Observations time (days). !
! Xobs Observations X-locations (grid coordinates). !
! Yobs Observations Y-locations (grid coordinates). !
! Zobs Observations Z-locations (grid coordinates). !
! cg_Gnorm Initial gradient vector normalization factor. !
! cg_Greduc Reduction in the gradient norm; excess cost function. !
! cg_QG Lanczos vector normalization factor. !
! cg_Ritz Eigenvalues of the Lanczos recurrence term Q(k)T(k). !
! cg_RitzErr Eigenvalues relative error. !
! cg_Tmatrix Lanczos recurrence symmetric, tridiagonal matrix. !
! cg_alpha Conjugate gradient coefficient. !
! cg_beta Conjugate gradient coefficient. !
! cg_delta Lanczos algorithm coefficient. !
! cg_gamma Lanczos algorithm coefficient. !
! cg_tau Conjugate gradient coefficient. !
! cg_zu Lanczos tridiagonal matrix, upper diagonal elements. !
! cg_zv Eigenvectors of Lanczos recurrence term Q(k)T(k). !
! haveADmod Logical switch indicating that there is representer !
! coefficients available to process in DAV(ng)%name !
! file. !
! haveNLmod Logical switch indicating that there is nonlinear !
! model data available to process in DAV(ng)%name !
! file. !
! haveTLmod Logical switch indicating that there is tangent !
! model data available to process in DAV(ng)%name !
! file. !
! haveObsMeta Logical switch indicating loading and processing of !
! observations meta field. !
! !
! Variables in observation space needed Desroziers et al. (2005) !
! statistics at observation locations: !
! !
! BgErr 4D-Var Background error standard deviation. !
! NLincrement 4D-Var NLM increment, analysis minus background. !
! innovation 4D-Var innovation, observations minus background. !
! residual 4D-Var residual, observation minus analysis. !
! !
! !
!=======================================================================
!
USE mod_param
!
implicit none
!
PUBLIC :: allocate_fourdvar
PUBLIC :: deallocate_fourdvar
PUBLIC :: initialize_fourdvar
!
!-----------------------------------------------------------------------
! Define T_FOURDVAR structure.
!-----------------------------------------------------------------------
!
TYPE T_FOURDVAR
integer , allocatable :: NobsSurvey(:)
integer , allocatable :: ObsCount(:)
integer , allocatable :: ObsReject(:)
real(r8), allocatable :: BackCost(:)
real(r8), allocatable :: Cost0(:)
real(r8), allocatable :: CostFun(:)
real(r8), allocatable :: CostFunOld(:)
real(r8), allocatable :: CostNorm(:)
real(r8), allocatable :: ObsCost(:)
real(r8), allocatable :: DataPenalty(:)
real(r8), allocatable :: NLPenalty(:)
real(r8), allocatable :: NLobsCost(:)
real(r8), allocatable :: SurveyTime(:)
real(r8), allocatable :: cg_pxout(:)
real(r8), allocatable :: cg_pxsave(:)
real(r8), allocatable :: cg_pxsum(:)
real(r8), allocatable :: cg_Dpxsum(:)
real(r8), allocatable :: tl_cg_pxsave(:)
real(r8), allocatable :: ad_cg_pxsave(:)
END TYPE T_FOURDVAR
!
TYPE (T_FOURDVAR), allocatable :: FOURDVAR(:)
!
!-----------------------------------------------------------------------
! Define other module variables.
!-----------------------------------------------------------------------
!
integer, allocatable :: ObsType(:)
integer, allocatable :: ObsState2Type(:)
integer, allocatable :: ObsType2State(:)
integer, allocatable :: ObsProv(:)
real(r8), allocatable :: ObsAngler(:)
real(r8), allocatable :: ObsErr(:)
real(r8), allocatable :: ObsMeta(:)
real(r8), allocatable :: ObsScale(:)
real(r8), allocatable :: ObsVetting(:)
real(r8), allocatable :: ObsVal(:)
real(dp), allocatable :: Tobs(:)
real(r8), allocatable :: Xobs(:)
real(r8), allocatable :: Yobs(:)
real(r8), allocatable :: Zobs(:)
real(r8), allocatable :: ADmodVal(:)
real(r8), allocatable :: NLmodVal(:)
real(r8), allocatable :: misfit(:)
real(r8), allocatable :: unvetted(:)
real(r8), allocatable :: uradial(:)
real(r8), allocatable :: vradial(:)
real(r8), allocatable :: uBgErr(:)
real(r8), allocatable :: vBgErr(:)
real(r8), allocatable :: TLmodVal(:)
real(r8), allocatable :: TLmodVal_S(:,:,:)
real(r8), allocatable :: BCKmodVal(:)
real(r8), allocatable :: Hbk(:,:)
real(r8), allocatable :: zcglwk(:,:,:)
real(r8), allocatable :: vcglwk(:,:,:)
real(r8), allocatable :: Jb0(:)
real(r8), allocatable :: vcglev(:,:,:)
real(r8), allocatable :: zgrad0(:,:)
real(r8), allocatable :: vgrad0(:)
real(r8), allocatable :: vgrad0s(:)
real(r8), allocatable :: gdgw(:)
real(r8), allocatable :: vgnorm(:)
real(r8), allocatable :: cg_innov(:)
real(r8), allocatable :: cg_dla(:,:)
!
!-----------------------------------------------------------------------
! Define observations parameters.
!-----------------------------------------------------------------------
!
! Switches indicating that at input observations vertical position were
! given in meters (Zobs < 0) so the fractional vertical level position
! is computed during extraction and written to Observation NetCDF file.
!
logical, allocatable :: Load_Zobs(:)
logical, allocatable :: wrote_Zobs(:)
!
! Maximum number of observations to process.
!
integer :: Mobs
!
! Number of model state variables to process.
!
integer, allocatable :: NstateVar(:)
!
! Number of observation types and names to process. If the observation
! operators for a specific application have a 1-to-1 association with
! the model state variables, NobsVar = NstateVar. However, if more
! that one state variable is needed to evaluate a particular type of
! observation (HF radials, travel time, pressure, etc.), a new class
! name is added to the state variable list to assess and differentiate
! its impact. Then, NobsVar = NstateVar + NextraObs.
!
integer, allocatable :: NobsVar(:)
character(len=40), allocatable :: ObsName(:)
!
! Number of extra-observation classes, extra-observation type indices,
! and extra-observation type names. It is used in observation operators
! that require more than one state variable or not directly associated
! with state variables.
!
integer :: NextraObs = 0
integer, allocatable :: ExtraIndex(:)
character(len=40), allocatable :: ExtraName(:)
!
! Number of interpolation weights and (I,J,K) indices offsets.
!
integer, parameter :: Nweights = 8
integer, parameter, dimension(Nweights) :: Ioffset = &
& (/ 0, 1, 0, 1, 0, 1, 0, 1 /)
integer, parameter, dimension(Nweights) :: Joffset = &
& (/ 0, 0, 1, 1, 0, 0, 1, 1 /)
integer, parameter, dimension(Nweights) :: Koffset = &
& (/ 0, 0, 0, 0, 1, 1, 1, 1 /)
!
! Size of observation NetCDF file unlimited dimension.
!
integer, allocatable :: Ndatum(:)
!
! Number of observations surveys available.
!
integer, allocatable :: Nsurvey(:)
!
! Observation surveys counter.
!
integer, allocatable :: ObsSurvey(:)
!
! Current number of observations processed.
!
integer, allocatable :: Nobs(:)
!
! Current starting and ending observation file index processed.
!
integer, allocatable :: NstrObs(:)
integer, allocatable :: NendObs(:)
!
! Background error covariance normalization method:
!
! [0] Exact, very expensive
! [1] Approximated, randomization
!
integer, allocatable :: Nmethod(:)
!
! Random number generation scheme for randomization:
!
! [0] Intrinsic F90 routine "randon_number"
! [1] Gaussian distributed deviates, numerical recipes
!
integer, allocatable :: Rscheme(:)
!
! Number of iterations in the randomization ensemble used to
! compute the background error covariance, B, normalization
! factors. This factors ensure that the diagonal elements of
! B are equal to unity (Fisher and Coutier, 1995).
!
integer :: Nrandom = 1000
!
! Number of Lanczos iterations used in the posterior analysis
! error covariance matrix estimation.
!
integer :: NpostI
!
! Outer loop to consider in the computatio of the observations
! impact or observations sensitivity, Nimpact =< Nouter. This
! facilitates the computations with multiple outer loop 4D-Var
! applications. The observation analysis needs to be computed
! separately for each outer loop. The full analysis for all
! outer loops are combined offline.
!
integer :: Nimpact
!
! Parameter to either process the eigenvector of the stabilized
! representer matrix to whe computing array modes (Nvct=1 less
! important eigenvector, Nvct=Ninner most important eigenvector)
! or cut-off eigenvector for the clipped analysis when (only
! eigenvector Nvct:Ninner to will be processed, others will be
! disgarded).
!
integer :: Nvct
!
! Optimal, normalized cost funtion minimum. If the statistical
! hypotheses between the background and observations errors
! is consistemt, the cost function value at the minimum, Jmin,
! is idealy equal to half the number of observations assimilated
! (Optimality=1=2*Jmin/Nobs), for a linear system.
!
real(r8), allocatable :: Optimality(:)
!
! Switch to activate the processing of model state at the observation
! locations.
!
logical, allocatable :: ProcessObs(:)
!
! Switch to activate writting of nonlinear model values at
! observations locations into observations NetCDF file.
!
logical, allocatable :: wrtNLmod(:)
!
! Sitch to process and write out observation accept/reject flag used
! for screening and quality control.
!
logical, allocatable :: wrtObsScale(:)
!
! Switch to activate writting of representer model values at
! observations locations into observations NetCDF file.
!
logical, allocatable :: wrtRPmod(:)
!
! Switch to activate writting of tangent linear model values at
! observations locations into observations NetCDF file.
!
logical, allocatable :: wrtTLmod(:)
!
! Switch to activate writting of initial and final model-observation
! misfit (innovation) vector into 4DVAR output NetCDF file.
!
logical, allocatable :: wrtMisfit(:)
!
! Swiches indicating that there is representer coeffiecients,
! nonlinear, and tangent linear model data available to process
! in 4DVAR NetCDF output file (DAV(ng)%name). At the beeginning,
! there is not data since this file has been just defined.
!
logical, allocatable :: haveADmod(:)
logical, allocatable :: haveNLmod(:)
logical, allocatable :: haveTLmod(:)
!
! Switch to indicate whether observations are present that
! require the specific meta field be loaded and processed.
!
logical, allocatable :: haveObsMeta(:)
!
!-----------------------------------------------------------------------
! Variables in observation space needed Desroziers et al. (2005)
! statistics at observation locations:
!-----------------------------------------------------------------------
!
! 4D-Var background error standard deviation [Mobs].
!
real(r8), allocatable :: BgErr(:)
!
! 4D-Var NLM increment vector, analysis minus background {Mobs].
!
real(r8), allocatable :: NLincrement(:)
!
! 4D-Var innovation vector, observations minus background [Mobs].
!
real(r8), allocatable :: innovation(:)
!
! 4D-Var residual vector, observation minus analysis [Mobs].
!
real(r8), allocatable :: residual(:)
!
!-----------------------------------------------------------------------
! Spatial convolutions parameters
!-----------------------------------------------------------------------
!
! Initial conditions, surface forcing and model error covariance: Full
! number of horizontal and vertical diffusion operator step for spatial
! convolutions.
!
integer, allocatable :: NHsteps(:,:)
integer, allocatable :: NVsteps(:,:)
!
! Boundary conditions error covariance: Full number of horizontal and
! vertical diffusion operator step for spatial convolutions.
!
integer, allocatable :: NHstepsB(:,:)
integer, allocatable :: NVstepsB(:,:)
!
! Initial conditions, surface forcing and model error covariance:
! Horizontal and vertical diffusion operator time-step size for
! spatial convolutions.
!
real(r8), allocatable :: DTsizeH(:,:)
real(r8), allocatable :: DTsizeV(:,:)
!
! Boundary conditions error covariance: Horizontal and vertical
! diffusion operator time-step size for spatial convolutions.
!
real(r8), allocatable :: DTsizeHB(:,:)
real(r8), allocatable :: DTsizeVB(:,:)
!
! Minimum and maximum Horizontal and vertical diffusion coefficients
! used in spatial convolutions.
!
real(r8), allocatable :: KhMin(:)
real(r8), allocatable :: KhMax(:)
real(r8), allocatable :: KvMin(:)
real(r8), allocatable :: KvMax(:)
!
!-----------------------------------------------------------------------
! Conjugate gradient parameters.
!-----------------------------------------------------------------------
!
! Conjugate gradient coefficients.
!
real(r8), allocatable :: cg_alpha(:,:)
real(r8), allocatable :: cg_beta(:,:)
real(r8), allocatable :: cg_tau(:,:)
!
! Ritz values maximum error limit.
!
real(r8) :: RitzMaxErr
!
! Converged Ritz eigenvalues used to approximate Hessian matrix during
! preconditioning.
!
real(r8), allocatable :: Ritz(:)
!
! Orthogonal eigenvectors of Lanczos recurrence term Q(k)T(k).
!
real(r8), allocatable :: cg_zv(:,:,:)
!
! Lanczos algorithm coefficients.
!
real(r8), allocatable :: cg_delta(:,:)
real(r8), allocatable :: cg_gamma(:,:)
!
! Initial gradient vector normalization factor.
!
real(r8), allocatable :: cg_Gnorm(:)
real(r8), allocatable :: cg_Gnorm_v(:)
real(r8), allocatable :: cg_Gnorm_y(:)
!
! Reduction in the gradient norm; excess cost function.
!
real(r8), allocatable :: cg_Greduc(:,:)
!
! Lanczos vector normalization factor.
!
real(r8), allocatable :: cg_QG(:,:)
!
! Lanczos recurrence symmetric, tridiagonal matrix, T(k).
!
real(r8), allocatable :: cg_Tmatrix(:,:)
real(r8), allocatable :: cg_zu(:,:)
!
! Eigenvalues of the Lanczos recurrence term Q(k)T(k)
! algorithm and their relative error (accuaracy).
!
real(r8), allocatable :: cg_Ritz(:,:)
real(r8), allocatable :: cg_RitzErr(:,:)
!
!-----------------------------------------------------------------------
! Descent algorithm parameters.
!-----------------------------------------------------------------------
!
! Switch to compute Hessian approximated eigenvalues and eigenvectors.
!
logical :: LhessianEV
!
! Switch switch to activate hot start of subsquent outerloops in the
! weak-constraint algorithms.
!
logical :: LhotStart
!
! Switch to activate conjugate gradient preconditioning.
!
logical :: Lprecond
!
! Switch to activate Ritz limited-memory preconditioning using the
! eigenpairs approximation of the Hessian matrix (Tshimanga et al.,
! 2008).
!
logical :: Lritz
!
! Switch that controls computation of the augumented contributions
! to the model error forcing terms in error_covariance.
!
logical :: LaugWeak=.FALSE.
!
!
! Number of converged Ritz eigenvalues. If input parameter NritzEV > 0,
! then nConvRitz=NritzEV. Therefore, only nConvRitz eigenpairs will
! used for preconditioning and the RitzMaxErr threshold is ignored.
!
integer :: NritzEV
integer :: nConvRitz = 0
!
! Weak contraint conjugate gradient norms.
!
real(r8) :: cg_gammam1 = 0.0_r8
real(r8) :: cg_sigmam1 = 0.0_r8
real(r8) :: cg_rnorm = 0.0_r8
!
! Upper bound on the relative error of the gradient when computing
! Hessian eigenvectors.
!
real(r8) :: GradErr
!
! Maximum error bound on Hessian eigenvectors. Note that even quite
! innacurate eigenvectors are usefull for pre-condtioning purposes.
!
real(r8) :: HevecErr
!
!------------------------------------------------------------------------
! Dot product parameters.
!------------------------------------------------------------------------
!
! Dot product between tangent linear and adjoint vectors.
!
real(r8) :: DotProduct
real(r8) :: adDotProduct
!
! Tangent linear model linearization check dot products.
!
integer :: ig1count ! counter for g1 dot product
integer :: ig2count ! counter for g2 dot product
real(r8), dimension(1000) :: g1
real(r8), dimension(1000) :: g2
!
!------------------------------------------------------------------------
! Weak constraint parameters.
!------------------------------------------------------------------------
!
! Switch to activate writing of weak constraint forings
! into the adjoint NetCDF file.
!
logical, allocatable :: WRTforce(:)
!
! Weak-constraint forcing time. A new variable is required since this
! forcing is delayed by nADJ time-steps.
!
real(r8), allocatable :: ForceTime(:)
!
CONTAINS
!
SUBROUTINE allocate_fourdvar
!
!=======================================================================
! !
! This routine allocates several variables in module "mod_fourdvar" !
! for all nested grids. !
! !
!=======================================================================
!
!-----------------------------------------------------------------------
! Allocate module variables due to nested grids.
!-----------------------------------------------------------------------
!
IF (.not.allocated(Load_Zobs)) THEN
allocate ( Load_Zobs(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(wrote_Zobs)) THEN
allocate ( wrote_Zobs(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(NstateVar)) THEN
allocate ( NstateVar(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(NobsVar)) THEN
allocate ( NobsVar(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(Ndatum)) THEN
allocate ( Ndatum(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(Nsurvey)) THEN
allocate ( Nsurvey(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(ObsSurvey)) THEN
allocate ( ObsSurvey(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(Nobs)) THEN
allocate ( Nobs(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(NstrObs)) THEN
allocate ( NstrObs(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(NendObs)) THEN
allocate ( NendObs(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(Nmethod)) THEN
allocate ( Nmethod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(Rscheme)) THEN
allocate ( Rscheme(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(Optimality)) THEN
allocate ( Optimality(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(ProcessObs)) THEN
allocate ( ProcessObs(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(wrtNLmod)) THEN
allocate ( wrtNLmod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(wrtObsScale)) THEN
allocate ( wrtObsScale(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(wrtRPmod)) THEN
allocate ( wrtRPmod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(wrtTLmod)) THEN
allocate ( wrtTLmod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(wrtMisfit)) THEN
allocate ( wrtMisfit(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(haveADmod)) THEN
allocate ( haveADmod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(haveNLmod)) THEN
allocate ( haveNLmod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(haveTLmod)) THEN
allocate ( haveTLmod(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(haveObsMeta)) THEN
allocate ( haveObsMeta(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(KhMin)) THEN
allocate ( KhMin(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(KhMax)) THEN
allocate ( KhMax(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(KvMin)) THEN
allocate ( KvMin(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(KvMax)) THEN
allocate ( KvMax(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(WRTforce)) THEN
allocate ( WRTforce(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
IF (.not.allocated(ForceTime)) THEN
allocate ( ForceTime(Ngrids) )
Dmem(1)=Dmem(1)+REAL(Ngrids,r8)
END IF
!
RETURN
END SUBROUTINE allocate_fourdvar
!
SUBROUTINE deallocate_fourdvar
!
!=======================================================================
! !
! This routine deallocates only the variables in module associated !
! with observation dimension parameters Ndatum, Nsurvey, and Mobs. !
! !
!=======================================================================
!
! Local variable declarations.
!
integer :: ng
!
!-----------------------------------------------------------------------
! Deallocate structure variables for each nested grids.
!-----------------------------------------------------------------------
!
DO ng=1,Ngrids
IF (allocated(FOURDVAR(ng) % NobsSurvey)) THEN
deallocate (FOURDVAR(ng) % NobsSurvey)
END IF
IF (allocated(FOURDVAR(ng) % SurveyTime)) THEN
deallocate (FOURDVAR(ng) % SurveyTime)
END IF
IF (allocated(FOURDVAR(ng) % cg_pxsave)) THEN
deallocate (FOURDVAR(ng) % cg_pxsave)
END IF
IF (allocated(FOURDVAR(ng) % cg_pxout)) THEN
deallocate (FOURDVAR(ng) % cg_pxout)
END IF
IF (allocated(FOURDVAR(ng) % cg_pxsum)) THEN
deallocate (FOURDVAR(ng) % cg_pxsum)
END IF
IF (allocated(FOURDVAR(ng) % cg_Dpxsum)) THEN
deallocate (FOURDVAR(ng) % cg_Dpxsum)
END IF
IF (allocated(FOURDVAR(ng) % tl_cg_pxsave)) THEN
deallocate (FOURDVAR(ng) % tl_cg_pxsave)
END IF
IF (allocated(FOURDVAR(ng) % ad_cg_pxsave)) THEN
deallocate (FOURDVAR(ng) % ad_cg_pxsave)
END IF
END DO
!
!-----------------------------------------------------------------------
! deallocate and initialize model and observation variables.
!-----------------------------------------------------------------------
!
IF (allocated(ObsAngler)) THEN
deallocate (ObsAngler)
END IF
IF (allocated(ObsType)) THEN
deallocate (ObsType)
END IF
IF (allocated(ObsProv)) THEN
deallocate (ObsProv)
END IF
IF (allocated(ObsErr)) THEN
deallocate (ObsErr)
END IF
IF (allocated(ObsMeta)) THEN
deallocate (ObsMeta)
END IF
IF (allocated(ObsScale)) THEN
deallocate (ObsScale)
END IF
IF (allocated(ObsVetting)) THEN
deallocate (ObsVetting)
END IF
IF (allocated(ObsVal)) THEN
deallocate (ObsVal)
END IF
IF (allocated(Tobs)) THEN
deallocate (Tobs)
END IF
IF (allocated(Xobs)) THEN
deallocate (Xobs)
END IF
IF (allocated(Yobs)) THEN
deallocate (Yobs)
END IF
IF (allocated(Zobs)) THEN
deallocate (Zobs)
END IF
IF (allocated(ADmodVal)) THEN
deallocate (ADmodVal)
END IF
IF (allocated(NLmodVal)) THEN
deallocate (NLmodVal)
END IF
IF (allocated(misfit)) THEN
deallocate (misfit)
END IF
IF (allocated(unvetted)) THEN
deallocate (unvetted)
END IF
IF (allocated(uradial)) THEN
deallocate (uradial)
END IF
IF (allocated(vradial)) THEN
deallocate (vradial)
END IF
IF (allocated(uBgErr)) THEN
deallocate (uBgErr)
END IF
IF (allocated(vBgErr)) THEN
deallocate (vBgErr)
END IF
IF (allocated(BgErr)) THEN
deallocate (BgErr)
END IF
IF (allocated(NLincrement)) THEN
deallocate (NLincrement)
END IF
IF (allocated(innovation)) THEN
deallocate (innovation)
END IF
IF (allocated(residual)) THEN
deallocate (residual)
END IF
IF (allocated(TLmodVal)) THEN
deallocate (TLmodVal)
END IF
IF (allocated(TLmodVal_S)) THEN
deallocate (TLmodVal_S)
END IF
IF (allocated(BCKmodVal)) THEN
deallocate (BCKmodVal)
END IF
IF (allocated(Hbk)) THEN
deallocate (Hbk)
END IF
!
! Deallocate and initialize weak constraint conjugate gradient vectors.
!
IF (allocated(zcglwk)) THEN
deallocate (zcglwk)
END IF
IF (allocated(vcglwk)) THEN
deallocate (vcglwk)
END IF
IF (allocated(vcglev)) THEN
deallocate (vcglev)
END IF
IF (allocated(zgrad0)) THEN
deallocate (zgrad0)
END IF
IF (allocated(vgrad0)) THEN
deallocate (vgrad0)
END IF
IF (allocated(vgrad0s)) THEN
deallocate (vgrad0s)
END IF
IF (allocated(cg_innov)) THEN
deallocate (cg_innov)
END IF
!
RETURN
END SUBROUTINE deallocate_fourdvar
!
SUBROUTINE initialize_fourdvar
!
!=======================================================================
! !
! This routine initializes several variables in module "mod_fourdvar" !
! for all nested grids. !
! !
!=======================================================================
!
USE mod_parallel
USE mod_scalars
USE mod_iounits
USE mod_ncparam
USE mod_netcdf
!
USE strings_mod, ONLY : FoundError
USE strings_mod, ONLY : uppercase
!
! Local variable declarations.
!
integer :: my_NobsVar, i, icount, ng
real(r8), parameter :: IniVal = 0.0_r8
!
character (len=*), parameter :: MyFile = &
& "ROMS/Modules/mod_fourdvar.F"//", initialize_fourdvar"
SourceFile=MyFile
!
!-----------------------------------------------------------------------
! Inquire observations NetCDF and determine the maximum dimension of
! several observations arrays.
!-----------------------------------------------------------------------
!
! Inquire about the size of the "datum" unlimitted dimension and
! "survey" dimension.
!
DO ng=1,Ngrids
SELECT CASE (OBS(ng)%IOtype)
CASE (io_nf90)
CALL netcdf_get_dim (ng, iNLM, TRIM(OBS(ng)%name))
CASE DEFAULT
IF (Master) WRITE (stdout,10) OBS(ng)%IOtype
exit_flag=3
END SELECT
IF (FoundError(exit_flag, NoError, 1642, MyFile)) RETURN
!
Ndatum(ng)=0
Nsurvey(ng)=0
DO i=1,n_dim
IF (TRIM(dim_name(i)).eq.'datum') then
Ndatum(ng)=dim_size(i)
ELSE IF (TRIM(dim_name(i)).eq.'survey') THEN
Nsurvey(ng)=dim_size(i)
ELSE IF (TRIM(dim_name(i)).eq.'state_variable') THEN
my_NobsVar=dim_size(i)
END IF
END DO
IF (Ndatum(ng).eq.0) THEN
WRITE (stdout,20) 'datum', TRIM(OBS(ng)%name)
exit_flag=2
RETURN
END IF
IF (Nsurvey(ng).eq.0) THEN
WRITE (stdout,20) 'survey', TRIM(OBS(ng)%name)
exit_flag=2
RETURN
END IF
END DO
!
!-----------------------------------------------------------------------
! Allocate module variables.
!-----------------------------------------------------------------------
!
!
! Allocate structure (1:Ngrids).
!
IF (.not.allocated( FOURDVAR)) THEN
allocate ( FOURDVAR(Ngrids) )
END IF
!
! Allocate variables in structure.
!
DO ng=1,Ngrids
!
! Number of state variables associated with data assimilation.
!
! zeta, ubar, vbar, Uvel, Vvel, Tvar(1:MT)
! Ustr, Vstr, Tsur(1:MT)
!
! The rest are ignored.
!
NstateVar(ng)=5+NT(ng)
NobsVar(ng)=NstateVar(ng)+NextraObs
IF (.not.allocated(FOURDVAR(ng) % NobsSurvey)) THEN
allocate ( FOURDVAR(ng) % NobsSurvey(Nsurvey(ng)) )
Dmem(ng)=Dmem(ng)+REAL(Nsurvey(ng),r8)
FOURDVAR(ng) % NobsSurvey = 0
END IF
IF (.not.allocated(FOURDVAR(ng) % SurveyTime)) THEN
allocate ( FOURDVAR(ng) % SurveyTime(Nsurvey(ng)) )
Dmem(ng)=Dmem(ng)+REAL(Nsurvey(ng),r8)
FOURDVAR(ng) % SurveyTime = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % BackCost)) THEN
allocate ( FOURDVAR(ng) % BackCost(0:NstateVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NstateVar(ng)+1,r8)
FOURDVAR(ng) % BackCost = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % Cost0)) THEN
allocate ( FOURDVAR(ng) % Cost0(Nouter) )
Dmem(ng)=Dmem(ng)+REAL(Nouter,r8)
FOURDVAR(ng) % Cost0 = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % CostFun)) THEN
allocate ( FOURDVAR(ng) % CostFun(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % CostFun = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % CostFunOld)) THEN
allocate ( FOURDVAR(ng) % CostFunOld(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % CostFunOld = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % CostNorm)) THEN
allocate ( FOURDVAR(ng) % CostNorm(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % CostNorm = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % ObsCost)) THEN
allocate ( FOURDVAR(ng) % ObsCost(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % ObsCost = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % DataPenalty)) THEN
allocate ( FOURDVAR(ng) % DataPenalty(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % DataPenalty = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % NLPenalty)) THEN
allocate ( FOURDVAR(ng) % NLPenalty(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % NLPenalty = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % NLobsCost)) THEN
allocate ( FOURDVAR(ng) % NLobsCost(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % NLobsCost = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % ObsCount)) THEN
allocate ( FOURDVAR(ng) % ObsCount(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % ObsCount = 0
END IF
IF (.not.allocated(FOURDVAR(ng) % ObsReject)) THEN
allocate ( FOURDVAR(ng) % ObsReject(0:NobsVar(ng)) )
Dmem(ng)=Dmem(ng)+REAL(NobsVar(ng)+1,r8)
FOURDVAR(ng) % ObsReject = 0
END IF
IF (.not.allocated(FOURDVAR(ng) % cg_pxout)) THEN
allocate ( FOURDVAR(ng) % cg_pxout(Nouter) )
Dmem(ng)=Dmem(ng)+REAL(Nouter)
FOURDVAR(ng) % cg_pxout = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % cg_pxsave)) THEN
allocate ( FOURDVAR(ng) % cg_pxsave(Ndatum(ng)+1) )
Dmem(ng)=Dmem(ng)+REAL(Ndatum(ng)+1,r8)
FOURDVAR(ng) % cg_pxsave = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % cg_pxsum)) THEN
allocate ( FOURDVAR(ng) % cg_pxsum(Ndatum(ng)+1) )
Dmem(ng)=Dmem(ng)+REAL(Ndatum(ng)+1,r8)
FOURDVAR(ng) % cg_pxsum = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % cg_Dpxsum)) THEN
allocate ( FOURDVAR(ng) % cg_Dpxsum(Ndatum(ng)+1) )
Dmem(ng)=Dmem(ng)+REAL(Ndatum(ng)+1,r8)
FOURDVAR(ng) % cg_Dpxsum = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % tl_cg_pxsave)) THEN
allocate ( FOURDVAR(ng) % tl_cg_pxsave(Ndatum(ng)+1) )
Dmem(ng)=Dmem(ng)+REAL(Ndatum(ng)+1,r8)
FOURDVAR(ng) % tl_cg_pxsave = IniVal
END IF
IF (.not.allocated(FOURDVAR(ng) % ad_cg_pxsave)) THEN
allocate ( FOURDVAR(ng) % ad_cg_pxsave(Ndatum(ng)+1) )
Dmem(ng)=Dmem(ng)+REAL(Ndatum(ng)+1,r8)
FOURDVAR(ng) % ad_cg_pxsave = IniVal
END IF
END DO
!
!-----------------------------------------------------------------------
! Read in number of observations available per survey at their times.
!-----------------------------------------------------------------------
!
DO ng=1,Ngrids
!
! Read in number of observations available per survey and the time of
! each observation survey.
!
SELECT CASE (OBS(ng)%IOtype)
CASE (io_nf90)
CALL netcdf_get_ivar (ng, iNLM, TRIM(OBS(ng)%name), &
& TRIM(Vname(1,idNobs)), &
& FOURDVAR(ng)%NobsSurvey, &
& start = (/1/), &
& total = (/Nsurvey(ng)/))
IF (FoundError(exit_flag, NoError, 1976, MyFile)) RETURN
!
CALL netcdf_get_fvar (ng, iNLM, TRIM(OBS(ng)%name), &
& TRIM(Vname(1,idOday)), &
& FOURDVAR(ng)%SurveyTime, &
& start = (/1/), &
& total = (/Nsurvey(ng)/))
IF (FoundError(exit_flag, NoError, 1983, MyFile)) RETURN
END SELECT
!
! Determine maximum size of observation arrays.
!
Mobs=MAXVAL(Ndatum)
END DO
!
!-----------------------------------------------------------------------
! Allocate and initialize model and observation variables.
!-----------------------------------------------------------------------
!
IF (.not.allocated(ObsAngler)) THEN
allocate ( ObsAngler(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsAngler = IniVal
IF (.not.allocated(ObsType)) THEN
allocate ( ObsType(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsType = 0
IF (.not.allocated(ObsState2Type)) THEN
allocate ( ObsState2Type(0:MAXVAL(NobsVar)) )
Dmem(1)=Dmem(1)+REAL(MAXVAL(NobsVar)+1,r8)
END IF
ObsState2Type = 0
IF (.not.allocated(ObsProv)) THEN
allocate ( ObsProv(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsProv = 0
IF (.not.allocated(ObsErr)) THEN
allocate ( ObsErr(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsErr = IniVal
IF (.not.allocated(ObsMeta)) THEN
allocate ( ObsMeta(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsMeta = IniVal
IF (.not.allocated(ObsScale)) THEN
allocate ( ObsScale(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsScale = IniVal
IF (.not.allocated(ObsVetting)) THEN
allocate ( ObsVetting(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsVetting = IniVal
IF (.not.allocated(ObsVal)) THEN
allocate ( ObsVal(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ObsVal = IniVal
IF (.not.allocated(Tobs)) THEN
allocate ( Tobs(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
Tobs = 0.0_dp
IF (.not.allocated(Xobs)) THEN
allocate ( Xobs(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
Xobs = IniVal
IF (.not.allocated(Yobs)) THEN
allocate ( Yobs(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
Yobs = IniVal
IF (.not.allocated(Zobs)) THEN
allocate ( Zobs(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
Zobs = IniVal
IF (.not.allocated(ADmodVal)) THEN
allocate ( ADmodVal(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
ADmodVal = IniVal
IF (.not.allocated(NLmodVal)) THEN
allocate ( NLmodVal(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
NLmodVal = IniVal
IF (.not.allocated(misfit)) THEN
allocate ( misfit(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
misfit = IniVal
IF (.not.allocated(unvetted)) THEN
allocate ( unvetted(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
unvetted = IniVal
IF (.not.allocated(uradial)) THEN
allocate ( uradial(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
uradial = IniVal
IF (.not.allocated(vradial)) THEN
allocate ( vradial(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
vradial = IniVal
IF (.not.allocated(uBgErr)) THEN
allocate ( uBgErr(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
uBgErr = IniVal
IF (.not.allocated(vBgErr)) THEN
allocate ( vBgErr(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
vBgErr = IniVal
IF (.not.allocated(BgErr)) THEN
allocate ( BgErr(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
BgErr = IniVal
IF (.not.allocated(NLincrement)) THEN
allocate ( NLincrement(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
NLincrement = IniVal
IF (.not.allocated(innovation)) THEN
allocate ( innovation(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
innovation = IniVal
IF (.not.allocated(residual)) THEN
allocate ( residual(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
residual = IniVal
IF (.not.allocated(TLmodVal)) THEN
allocate ( TLmodVal(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
TLmodVal = IniVal
!
! The following arrays are only needed and used by the master node.
! In order to avoid hogging memory penalties in the other nodes, they
! are only allocated by the master node.
!
IF (.not.allocated(TLmodVal_S)) THEN
allocate ( TLmodVal_S(Mobs,Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Mobs*Ninner*Nouter,r8)
END IF
TLmodVal_S = IniVal
IF (.not.allocated(BCKmodVal)) THEN
allocate ( BCKmodVal(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
BCKmodVal = IniVal
IF (.not.allocated(Hbk)) THEN
allocate ( Hbk(Mobs,Nouter) )
Dmem(1)=Dmem(1)+REAL(Mobs*Nouter,r8)
END IF
Hbk = IniVal
!
! Allocate and initialize observation types names and indices.
! Notice that a mapping from state-to-type (ObsState2Type) and its
! inverse type-to-state (ObsType2State) indices are needed because
! the User is allowed to add extra-observation operators with
! nonsequential type enumeration.
!
! Both mapping arrays ObsState2Type and ObsType2State have a zero
! array element to allow applications with no extra-observations
! to work with their zero associated state index (as initialized
! in mod_ncparam.F). For example, if the index "isRadial" is not
! redefined below, the following assignment
!
! ObsState2Type(isRadial)=ObsState2Type(0)=0
! ObsType2State(isRadial)=ObsType2State(0)=0
!
! is still legal with isRadial=0. It avoids a Fortran segmentation
! violation (i.e., subscript #1 of the array ObsState2Type has
! value 0 which is less than the lower bound of 1). Sorry for
! the awkward logic but we need a generic way to specify extra-
! observation operators.
!
IF (.not.allocated(ObsName)) THEN
allocate ( ObsName(MAXVAL(NobsVar)) )
Dmem(1)=Dmem(1)+REAL(MAXVAL(NobsVar),r8)
END IF
icount=MAXVAL(NstateVar)
ObsState2Type(0)=0
DO i=1,icount ! 5+NT
ObsState2Type(i)=i
ObsName(i)=TRIM(Vname(1,idSvar(i)))
END DO
IF (NextraObs.gt.0) THEN
DO i=1,NextraObs
icount=icount+1
ObsName(icount)=TRIM(ExtraName(i))
SELECT CASE (TRIM(uppercase(ExtraName(i))))
CASE ('RADIAL')
isRadial=icount
ObsState2Type(icount)=ExtraIndex(i)
END SELECT
END DO
END IF
!
! Set inverse mapping type-to-state.
!
IF (.not.allocated(ObsType2State)) THEN
allocate ( ObsType2State(0:MAXVAL(ObsState2Type)) )
Dmem(1)=Dmem(1)+REAL(MAXVAL(ObsState2Type)+1,r8)
END IF
ObsType2State(0)=0
ObsType2State(1:MAXVAL(ObsState2Type))=-1
DO i=1,MAXVAL(NstateVar)
ObsType2State(i)=i
END DO
IF (NextraObs.gt.0) THEN
DO i=1,NextraObs
icount=ExtraIndex(i)
SELECT CASE (TRIM(uppercase(ExtraName(i))))
CASE ('RADIAL')
ObsType2State(icount)=isRadial
END SELECT
END DO
END IF
!
! Allocate and initialize weak constraint conjugate gradient vectors.
!
IF (.not.allocated(zcglwk)) THEN
allocate ( zcglwk(Mobs+1,Ninner+1,Nouter) )
Dmem(1)=Dmem(1)+REAL((Mobs+1)*(Ninner+1)*Nouter,r8)
END IF
zcglwk = IniVal
IF (.not.allocated(vcglwk)) THEN
allocate ( vcglwk(Mobs+1,Ninner+1,Nouter) )
Dmem(1)=Dmem(1)+REAL((Mobs+1)*(Ninner+1)*Nouter,r8)
END IF
vcglwk = IniVal
IF (.not.allocated(Jb0)) THEN
allocate ( Jb0(0:Nouter) )
Dmem(1)=Dmem(1)+REAL(Nouter+1,r8)
END IF
Jb0 = IniVal
IF (.not.allocated(vcglev)) THEN
allocate ( vcglev(Mobs+1,Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL((Mobs+1)*Ninner*Nouter,r8)
END IF
vcglev = IniVal
IF (.not.allocated(zgrad0)) THEN
allocate ( zgrad0(Mobs+1,Nouter) )
Dmem(1)=Dmem(1)+REAL((Mobs+1)*Nouter,r8)
END IF
zgrad0 = IniVal
IF (.not.allocated(vgrad0)) THEN
allocate ( vgrad0(Mobs+1) )
Dmem(1)=Dmem(1)+REAL(Mobs+1,r8)
END IF
vgrad0 = IniVal
IF (.not.allocated(vgrad0s)) THEN
allocate ( vgrad0s(Mobs+1) )
Dmem(1)=Dmem(1)+REAL(Mobs+1,r8)
END IF
vgrad0s = IniVal
IF (.not.allocated(gdgw)) THEN
allocate ( gdgw(Mobs) )
Dmem(1)=Dmem(1)+REAL(Mobs,r8)
END IF
gdgw = IniVal
IF (.not.allocated(vgnorm)) THEN
allocate ( vgnorm(Nouter) )
Dmem(1)=Dmem(1)+REAL(Nouter,r8)
END IF
vgnorm = IniVal
IF (.not.allocated(cg_innov)) THEN
allocate ( cg_innov(Mobs+1) )
Dmem(1)=Dmem(1)+REAL(Mobs+1,r8)
END IF
cg_innov = IniVal
IF (.not.allocated(cg_dla)) THEN
allocate ( cg_dla(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_dla = IniVal
!
! Allocate convolution parameters.
!
IF (.not.allocated(NHsteps)) THEN
allocate ( NHsteps(2,MstateVar) )
Dmem(1)=Dmem(1)+2.0_r8*REAL(MstateVar,r8)
END IF
NHsteps = 0
IF (.not.allocated(NVsteps)) THEN
allocate ( NVsteps(2,MstateVar) )
Dmem(1)=Dmem(1)+2.0_r8*REAL(MstateVar,r8)
END IF
NVsteps = 0
IF (.not.allocated(DTsizeH)) THEN
allocate ( DTsizeH(2,MstateVar) )
Dmem(1)=Dmem(1)+2.0_r8*REAL(MstateVar,r8)
END IF
DTsizeH = IniVal
IF (.not.allocated(DTsizeV)) THEN
allocate ( DTsizeV(2,MstateVar) )
Dmem(1)=Dmem(1)+2.0_r8*REAL(MstateVar,r8)
END IF
DTsizeV = IniVal
IF (.not.allocated(NVstepsB)) THEN
allocate ( NVstepsB(4,MstateVar) )
Dmem(1)=Dmem(1)+4.0_r8*REAL(MstateVar,r8)
END IF
NVstepsB = 0
IF (.not.allocated(NHstepsB)) THEN
allocate ( NHstepsB(4,MstateVar) )
Dmem(1)=Dmem(1)+4.0_r8*REAL(MstateVar,r8)
END IF
NHstepsB = 0
IF (.not.allocated(DTsizeHB)) THEN
allocate ( DTsizeHB(4,MstateVar) )
Dmem(1)=Dmem(1)+4.0_r8*REAL(MstateVar,r8)
END IF
DTsizeHB = IniVal
IF (.not.allocated(DTsizeVB)) THEN
allocate ( DTsizeVB(4,MstateVar) )
Dmem(1)=Dmem(1)+4.0_r8*REAL(MstateVar,r8)
END IF
DTsizeVB = IniVal
!
! Allocate conjugate gradient variables.
!
IF (.not.allocated(cg_alpha)) THEN
allocate ( cg_alpha(0:Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL((Ninner+1)*Nouter,r8)
END IF
cg_alpha = IniVal
IF (.not.allocated(cg_beta)) THEN
allocate ( cg_beta(Ninner+1,Nouter) )
Dmem(1)=Dmem(1)+REAL((Ninner+1)*Nouter,r8)
END IF
cg_beta = IniVal
IF (.not.allocated(cg_tau)) THEN
allocate ( cg_tau(0:Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL((Ninner+1)*Nouter,r8)
END IF
cg_tau = IniVal
IF (.not.allocated(Ritz)) THEN
allocate ( Ritz(Ninner+1) )
Dmem(1)=Dmem(1)+REAL(Ninner+1,r8)
END IF
Ritz = IniVal
IF (.not.allocated(cg_zv)) THEN
allocate ( cg_zv(Ninner,Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Ninner*Nouter,r8)
END IF
cg_zv = IniVal
IF (.not.allocated(cg_delta)) THEN
allocate ( cg_delta(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_delta = IniVal
IF (.not.allocated(cg_gamma)) THEN
allocate ( cg_gamma(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_gamma = IniVal
IF (.not.allocated(cg_Gnorm)) THEN
allocate ( cg_Gnorm(Nouter) )
Dmem(1)=Dmem(1)+REAL(Nouter,r8)
END IF
cg_Gnorm = IniVal
IF (.not.allocated(cg_Gnorm_v)) THEN
allocate ( cg_Gnorm_v(Nouter) )
Dmem(1)=Dmem(1)+REAL(Nouter,r8)
END IF
cg_Gnorm_v = IniVal
IF (.not.allocated(cg_Gnorm_y)) THEN
allocate ( cg_Gnorm_y(Nouter) )
Dmem(1)=Dmem(1)+REAL(Nouter,r8)
END IF
cg_Gnorm_y = IniVal
IF (.not.allocated(cg_Greduc)) THEN
allocate ( cg_Greduc(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_Greduc = IniVal
IF (.not.allocated(cg_QG)) THEN
allocate ( cg_QG(Ninner+1,Nouter) )
Dmem(1)=Dmem(1)+REAL((Ninner+1)*Nouter,r8)
END IF
cg_QG = IniVal
IF (.not.allocated(cg_Tmatrix)) THEN
allocate ( cg_Tmatrix(Ninner,3) )
Dmem(1)=Dmem(1)+3.0_r8*REAL(Ninner,r8)
END IF
cg_Tmatrix = IniVal
IF (.not.allocated(cg_Ritz)) THEN
allocate ( cg_Ritz(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_Ritz = IniVal
IF (.not.allocated(cg_RitzErr)) THEN
allocate ( cg_RitzErr(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_RitzErr = IniVal
IF (.not.allocated(cg_zu)) THEN
allocate ( cg_zu(Ninner,Nouter) )
Dmem(1)=Dmem(1)+REAL(Ninner*Nouter,r8)
END IF
cg_zu = IniVal
!
!-----------------------------------------------------------------------
! Initialize various variables.
!-----------------------------------------------------------------------
!
DO ng=1,Ngrids
Load_Zobs(ng)=.FALSE.
ProcessObs(ng)=.FALSE.
haveObsMeta(ng)=.FALSE.
wrote_Zobs(ng)=.FALSE.
wrtMisfit(ng)=.FALSE.
wrtNLmod(ng)=.FALSE.
wrtObsScale(ng)=.FALSE.
wrtRPmod(ng)=.FALSE.
wrtTLmod(ng)=.FALSE.
KhMin(ng)=1.0_r8
KhMax(ng)=1.0_r8
KvMin(ng)=1.0_r8
KvMax(ng)=1.0_r8
Optimality(ng)=0.0_r8
ForceTime(ng)=0.0_r8
END DO
!
10 FORMAT (/,' INITIALIZE_FOURDVAR - Illegal output type,', &
& ' io_type = ',i0)
20 FORMAT (/,' INITIALIZE_FOURDVAR - error inquiring dimension:', &
& 1x,a,2x,'in input NetCDF file: ',a)
RETURN
END SUBROUTINE initialize_fourdvar
END MODULE mod_fourdvar