module control_vectors
!$$$ module documentation block
! . . . .
! module: control_vectors
! prgmmr: tremolet
!
! abstract: define control vectors and basic operators
!
! program history log:
! 2007-04-13 tremolet - initial code
! 2007-08-03 todling - using get_lun for file unit definition
! 2008-01-04 tremolet - improve allocate/deallocate
! 2008-12-29 todling - add omp to various loops
! 2009-01-27 todling - rename prt_norms to prevent IBM compiler confusion
! 2009-08-12 lueken - updated documentation
! 2009-09-20 parrish - add pointer variable a_en to definition of type control_state
! also, add module variables n_ens
! 2009-11-10 todling - remove redundant dot products
! 2010-02-17 treadon - fix in predictor part of random vector
! 2010-02-20 parrish - add functions dplevs_ens for use with dual-resolution hybrid ensemble option.
! 2010-03-11 zhu - add changes for generalizing control variables,i.e.,
! init_anacv,nvars,nrf*,allocate_cs,deallocate_cs,
! assign_cs2array,assign_array2cs
! 2010-05-01 todling - introduce gsi_bundle
! 2010-05-05 derber - omp commands removed
! 2010-05-22 todling - add a wired-in set of variables composing a motley (not fully part of CVector)
! 2010-05-28 todling - remove all nrf2/3_VAR-specific "pointers"
! 2011-07-04 todling - fixes to run either single or double precision
! 2013-05-20 zhu - add aircraft temperature bias correction coefficients as control variables
! 2016-02-15 Johnson, Y. Wang, X. Wang - add variables to control reading
! state variables for radar DA. POC: xuguang.wang@ou.edu
!
! subroutines included:
! sub init_anacv
! sub final_anacv
! sub setup_control_vectors
! sub allocate_cv
! sub deallocate_cv
! sub assign_scalar2cv
! sub assign_cv2cv
! sub assign_array2cv
! sub assign_cv2array
! sub qdot_prod_vars
! sub prt_norms
! sub prt_norms_vars
! sub axpy
! sub random_cv
! sub write_cv
! sub read_cv
! sub inquire_cv
!
! functions included:
! dplevs_ens
! qdot_prod_sub
! dot_prod_cv
! qdot_prod_cv
! qdot_prod_cv_eb
! maxval_cv
! qdot_product
!
! variable definitions:
! def n_ens - number of ensemble perturbations (=0 except when hybrid ensemble option turned on)
! def imp_physics - type of microphysics used in the GFS. 99: Zhao-Carr, 11: GFDL
! def lupp - if T, UPP is used and additional variables are output
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use kinds, only: r_kind,r_double,r_single,i_kind,r_quad
use mpimod, only: mpi_comm_world,mpi_max,mpi_rtype,mype,npe,ierror
use constants, only: zero, one, two, three, zero_quad, tiny_r_kind
use gsi_4dvar, only: iadatebgn
use file_utility, only : get_lun
use mpl_allreducemod, only: mpl_allreduce
use hybrid_ensemble_parameters, only: beta_s0,l_hyb_ens
use hybrid_ensemble_parameters, only: grd_ens
use constants, only : max_varname_length
use m_rerank, only : rerank
use GSI_BundleMod, only : GSI_BundleCreate
use GSI_BundleMod, only : GSI_BundleSet
use GSI_BundleMod, only : GSI_Bundle
use GSI_BundleMod, only : GSI_BundleGetPointer
use GSI_BundleMod, only : dplevs => GSI_BundleDpLevs
use GSI_BundleMod, only : GSI_BundleUnSet
use GSI_BundleMod, only : GSI_BundleDestroy
use GSI_BundleMod, only : GSI_Grid
use GSI_BundleMod, only : GSI_GridCreate
use mpeu_util, only: gettablesize
use mpeu_util, only: gettable
implicit none
save
private
!
! Public functions
!
public control_vector
public allocate_cv
public deallocate_cv
public assignment(=)
public dot_product
public prt_control_norms, axpy, random_cv, setup_control_vectors, &
write_cv, read_cv, inquire_cv, maxval, qdot_prod_sub, init_anacv, &
final_anacv
!
! Public variables
!
public cvars2d, cvars3d, cvarsmd, evars2d, evars3d, nrf_var
public nc2d ! number of 2d static control fields
public nc3d ! number of 3d static control fields
public mvars ! number of motley fields
public nrf ! total number of static control fields
public nvars ! total number of static plus motley fields
public nrf_3d ! when .t., indicates 3d-fields
public as3d ! normalized scale factor for background error 3d-variables
public as2d ! normalized scale factor for background error 2d-variables
public atsfc_sdv ! standard deviation of surface temperature error over (1) land (and (2) ice
public an_amp0 ! multiplying factors on reference background error variances
public imp_physics ! type of GFS microphysics
public lupp ! when .t., UPP is used and extra variables are output
public nrf2_loc,nrf3_loc,nmotl_loc ! what are these for??
public ntracer
public :: w_exist ! w will be used in the control variables ,only for
! wrf_mass_region =.true.
public :: dbz_exist ! dbz will be used in the control variables ,only for
! wrf_mass_region =.true.
type control_vector
integer(i_kind) :: lencv
real(r_kind), pointer :: values(:) => NULL()
type(GSI_Grid) :: grid_step
type(GSI_Bundle), pointer :: step(:)
type(GSI_Bundle), pointer :: motley(:)
type(GSI_Grid) :: grid_aens
type(GSI_Bundle), pointer :: aens(:,:)
real(r_kind), pointer :: predr(:) => NULL()
real(r_kind), pointer :: predp(:) => NULL()
real(r_kind), pointer :: predt(:) => NULL()
logical :: lallocated = .false.
end type control_vector
character(len=*),parameter:: myname='control_vectors'
integer(i_kind) :: nclen,nclen1,nsclen,npclen,ntclen,nrclen,nsubwin,nval_len
integer(i_kind) :: latlon11,latlon1n,lat2,lon2,nsig,n_ens
integer(i_kind) :: nval_lenz_en,imp_physics
logical :: lsqrtb,lupp
integer(i_kind) :: m_vec_alloc, max_vec_alloc, m_allocs, m_deallocs
logical,allocatable,dimension(:):: nrf_3d
integer(i_kind),allocatable,dimension(:):: nrf2_loc,nrf3_loc,nmotl_loc
integer(i_kind) nrf,nvars
integer(i_kind) ntracer
integer(i_kind) :: nc2d,nc3d,mvars
character(len=max_varname_length),allocatable,dimension(:) :: nrf_var ! names of all variables
character(len=max_varname_length),allocatable,dimension(:) :: cvars2d ! 2-d fields for static CV
character(len=max_varname_length),allocatable,dimension(:) :: cvars3d ! 3-d fields for static CV
character(len=max_varname_length),allocatable,dimension(:) :: evars2d ! 2-d fields for ensemble CV
character(len=max_varname_length),allocatable,dimension(:) :: evars3d ! 3-d fields for ensemble CV
character(len=max_varname_length),allocatable,dimension(:) :: cvarsmd ! motley variable names
real(r_kind) ,allocatable,dimension(:) :: as3d
real(r_kind) ,allocatable,dimension(:) :: as2d
real(r_kind) ,allocatable,dimension(:) :: atsfc_sdv
real(r_kind) ,allocatable,dimension(:) :: an_amp0
integer(i_kind) mx_clouds
parameter(mx_clouds=10)
character(len=10) cloud_var(mx_clouds)
logical w_exist
logical dbz_exist
logical :: llinit = .false.
! ----------------------------------------------------------------------
INTERFACE ASSIGNMENT (=)
MODULE PROCEDURE assign_scalar2cv, assign_array2cv, assign_cv2array, assign_cv2cv
END INTERFACE
INTERFACE DOT_PRODUCT
MODULE PROCEDURE dot_prod_cv,qdot_prod_cv,qdot_prod_cv_eb
END INTERFACE
INTERFACE MAXVAL
MODULE PROCEDURE maxval_cv
END INTERFACE
INTERFACE PRT_CONTROL_NORMS
MODULE PROCEDURE prt_norms
END INTERFACE
! ----------------------------------------------------------------------
contains
! ----------------------------------------------------------------------
subroutine setup_control_vectors(ksig,klat,klon,katlon11,katlon1n, &
ksclen,kpclen,ktclen,kclen,ksubwin,kval_len,ldsqrtb,k_ens,&
kval_lenz_en)
!$$$ subprogram documentation block
! . . . .
! subprogram: setup_control_vectors
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
! 2009-09-20 parrish - add optional input variable k_ens, which communicates size
! of ensemble used when hybrid ensemble option is turned on.
! 2010-05-19 todling - k_ens no longer optional
! 2010-05-23 todling - move lev pointer init from jfunc here
!
! input argument list:
! ksig
! klat,klon
! katlon11
! katlon1n
! ksclen
! kpclen
! ktclen
! kclen
! ksubwin
! kval_len
! ldsqrtb
! k_ens - if applicable, size of ensemble used in hybrid ensemble option
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind) , intent(in ) :: ksig,klat,klon,katlon11,katlon1n, &
ksclen,kpclen,ktclen,kclen,ksubwin,kval_len,k_ens,&
kval_lenz_en
logical , intent(in ) :: ldsqrtb
nsig=ksig
lat2=klat
lon2=klon
latlon11=katlon11
latlon1n=katlon1n
nsclen=ksclen
npclen=kpclen
ntclen=ktclen
nrclen=nsclen+npclen+ntclen
nclen =kclen
nclen1=nclen-nrclen
nsubwin=ksubwin
nval_len=kval_len
lsqrtb=ldsqrtb
n_ens=k_ens
nval_lenz_en=kval_lenz_en
llinit = .true.
m_vec_alloc=0
max_vec_alloc=0
m_allocs=0
m_deallocs=0
call inquire_cv
return
end subroutine setup_control_vectors
! ----------------------------------------------------------------------
subroutine init_anacv
!$$$ subprogram documentation block
! . . . .
! subprogram: init_anacv
! prgmmr: zhu org: np23 date: 2008-03-29
!
! abstract: read in control variables information
!
! program history log:
! 2010-03-11 zhu - initial code
! 2010-05-30 todling - revamp initial code
! 2014-02-11 todling - rank-2 must have lev=1, anything else is rank-3
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine: ibm rs/6000 sp
!
implicit none
!character(len=*),parameter:: rcname='anavinfo.txt'
character(len=*),parameter:: rcname='anavinfo' ! filename should have extension
character(len=*),parameter:: tbname='control_vector::'
character(len=256),allocatable,dimension(:):: utable
character(len=20) var,source,funcof
character(len=*),parameter::myname_=myname//'*init_anacv'
integer(i_kind) luin,ii,ntot
integer(i_kind) ilev, itracer
real(r_kind) aas,amp
integer(i_kind) icloud
! load file
luin=get_lun()
open(luin,file=rcname,form='formatted')
! Scan file for desired table first
! and get size of table
call gettablesize(tbname,luin,ntot,nvars)
! Get contents of table
allocate(utable(nvars))
call gettable(tbname,luin,ntot,nvars,utable)
! release file unit
close(luin)
! Retrieve each token of interest from table and define
! variables participating in state vector
! Count variables first
nc3d=0; nc2d=0;mvars=0
do ii=1,nvars
read(utable(ii),*) var, ilev, itracer, aas, amp, source, funcof
if(trim(adjustl(source))=='motley') then
mvars=mvars+1
else
if(ilev==1) then
nc2d=nc2d+1
else
nc3d=nc3d+1
endif
endif
enddo
allocate(nrf_var(nvars),cvars3d(nc3d),cvars2d(nc2d))
allocate(as3d(nc3d),as2d(nc2d))
allocate(cvarsmd(mvars))
allocate(atsfc_sdv(mvars))
allocate(an_amp0(nvars))
! want to rid code from the following ...
nrf=nc2d+nc3d
allocate(nrf_3d(nrf),nrf2_loc(nc2d),nrf3_loc(nc3d),nmotl_loc(max(1,mvars)))
! Now load information from table
nc3d=0;nc2d=0;mvars=0
nrf_3d=.false.
do ii=1,nvars
read(utable(ii),*) var, ilev, itracer, aas, amp, source, funcof
if(trim(adjustl(source))=='motley') then
mvars=mvars+1
cvarsmd(mvars)=trim(adjustl(var))
nmotl_loc(mvars)=ii
atsfc_sdv(mvars)=aas
else
if(ilev==1) then
nc2d=nc2d+1
cvars2d(nc2d)=trim(adjustl(var))
nrf2_loc(nc2d)=ii ! rid of soon
as2d(nc2d)=aas
else
nc3d=nc3d+1
cvars3d(nc3d)=trim(adjustl(var))
nrf3_loc(nc3d)=ii ! rid of soon
nrf_3d(ii)=.true.
as3d(nc3d)=aas
endif
endif
nrf_var(ii)=trim(adjustl(var))
if(amp>zero) then
an_amp0(ii)=amp
else
an_amp0(ii)=one/three
endif
enddo
deallocate(utable)
! right now, ens is made ideantical to static CV
allocate(evars2d(nc2d),evars3d(nc3d))
evars2d=cvars2d
evars3d=cvars3d
if (mype==0) then
write(6,*) myname_,': 2D-CONTROL VARIABLES ARE ', cvars2d
write(6,*) myname_,': 3D-CONTROL VARIABLES ARE ', cvars3d
write(6,*) myname_,': MOTLEY CONTROL VARIABLES ', cvarsmd
write(6,*) myname_,': ALL CONTROL VARIABLES ', nrf_var
end if
w_exist=.false.
dbz_exist=.false.
cloud_var="none"
icloud=0
do ii=1,nc3d
if(mype == 0 ) write(6,*)"anacv cvars3d is ",cvars3d(ii)
if(trim(cvars3d(ii)) == 'w'.or.trim(cvars3d(ii))=='W') w_exist=.true.
if(trim(cvars3d(ii))=='dbz'.or.trim(cvars3d(ii))=='DBZ') then
dbz_exist=.true.
endif
if(trim(cvars3d(ii))=='qr'.or.trim(cvars3d(ii))=='QR') then
icloud=icloud+1
cloud_var(icloud)="qr"
endif
if(trim(cvars3d(ii))=='ql'.or.trim(cvars3d(ii))=='QL') then
icloud=icloud+1
cloud_var(icloud)="ql"
endif
if(trim(cvars3d(ii))=='qs'.or.trim(cvars3d(ii))=='QS') then
icloud=icloud+1
cloud_var(icloud)="qs"
endif
if(trim(cvars3d(ii))=='qi'.or.trim(cvars3d(ii))=='QI') then
icloud=icloud+1
cloud_var(icloud)="qi"
endif
if(trim(cvars3d(ii))=='qg'.or.trim(cvars3d(ii))=='QG') then
icloud=icloud+1
cloud_var(icloud)="qg"
endif
enddo ! ii
imp_physics=99
lupp = .false.
end subroutine init_anacv
subroutine final_anacv
implicit none
deallocate(nrf_var)
deallocate(nrf_3d,nrf2_loc,nrf3_loc,nmotl_loc)
deallocate(as3d,as2d)
deallocate(an_amp0)
deallocate(atsfc_sdv)
deallocate(cvarsmd)
deallocate(cvars2d,cvars3d)
deallocate(evars2d,evars3d)
end subroutine final_anacv
! ----------------------------------------------------------------------
subroutine allocate_cv(ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: allocate_cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
! 2009-09-20 parrish - add optional allocation of hybrid ensemble control variable a_en
! 2010-02-20 parrish - add structure variable grd_ens as part of changes for dual-resolution
! hybrid ensemble system.
! 2010-02-25 zhu - use nrf_var and nrf_3d to specify the order control variables
! 2010-05-01 todling - update to use gsi_bundle
! 2010-05-17 todling - add back ens control; w/ a twist from original
! 2010-05-22 todling - add support for motley variables
!
! input argument list:
!
! output argument list:
! ycv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use hybrid_ensemble_parameters, only: grd_ens
implicit none
type(control_vector), intent( out) :: ycv
integer(i_kind) :: ii,jj,nn,ndim,ierror,n_step,n_aens
character(len=256)::bname
character(len=max_varname_length)::ltmp(1)
type(gsi_grid) :: grid_motley
if (ycv%lallocated) then
write(6,*)'allocate_cv: vector already allocated'
call stop2(108)
end if
ycv%lallocated=.true.
ycv%lencv = nclen
ALLOCATE(ycv%values(ycv%lencv))
! If so, define grid of regular control vector
n_step=0
! if (beta_s0>tiny_r_kind) then
ALLOCATE(ycv%step(nsubwin))
call GSI_GridCreate(ycv%grid_step,lat2,lon2,nsig)
if (lsqrtb) then
n_step=nval_len
else
n_step=nc3d*latlon1n+nc2d*latlon11
endif
if(mvars>0) then
ALLOCATE(ycv%motley(nsubwin))
call GSI_GridCreate(grid_motley,lat2,lon2,-1) ! this is sort of wired-in
endif
! endif
! If so, define grid of ensemble control vector
n_aens=0
if (l_hyb_ens) then
ALLOCATE(ycv%aens(nsubwin,n_ens))
call GSI_GridCreate(ycv%grid_aens,grd_ens%lat2,grd_ens%lon2,grd_ens%nsig)
if (lsqrtb) then
n_aens=nval_lenz_en
else
n_aens=grd_ens%latlon11*grd_ens%nsig
endif
endif
! Loop over sub-windows in time
ii=0; ndim=0
do jj=1,nsubwin
! Set static part of control vector (non-ensemble-based)
! if (beta_s0>tiny_r_kind) then
ycv%step(jj)%values => ycv%values(ii+1:ii+n_step)
write(bname,'(a,i3.3)') 'Static Control Bundle subwin-',jj
call GSI_BundleSet(ycv%step(jj),ycv%grid_step,bname,ierror,names2d=cvars2d,names3d=cvars3d,bundle_kind=r_kind)
if (ierror/=0) then
write(6,*)'allocate_cv: error alloc(static bundle)'
call stop2(109)
endif
ndim=ndim+ycv%step(jj)%ndim
ii=ii+n_step
! If motley variables needed (these don't contribute to CV size)
! Presently, there are thankfully only 2d-fields in this category
if(mvars>0) then
write(bname,'(a,i3.3,a,i4.4)') 'Motley Control Bundle subwin-',jj
call GSI_BundleCreate(ycv%motley(jj),grid_motley,bname,ierror,names2d=cvarsmd,bundle_kind=r_kind)
if (ierror/=0) then
write(6,*)'allocate_cv: error alloc(motley bundle)'
call stop2(109)
endif
endif
! endif ! beta_s0
! Set ensemble-based part of control vector
if (l_hyb_ens) then
ltmp(1)='a_en'
do nn=1,n_ens
ycv%aens(jj,nn)%values => ycv%values(ii+1:ii+n_aens)
write(bname,'(a,i3.3,a,i4.4)') 'Ensemble Control Bundle subwin-',jj,' and member-',nn
call GSI_BundleSet(ycv%aens(jj,nn),ycv%grid_aens,bname,ierror,names3d=ltmp,bundle_kind=r_kind)
if (ierror/=0) then
write(6,*)'allocate_cv: error alloc(ensemble bundle)'
call stop2(109)
endif
ndim=ndim+ycv%aens(jj,nn)%ndim
ii=ii+n_aens
enddo
endif
enddo
ycv%predr => ycv%values(ii+1:ii+nsclen)
ii=ii+nsclen
ycv%predp => ycv%values(ii+1:ii+npclen)
ii=ii+npclen
if (ntclen>0) then
ycv%predt => ycv%values(ii+1:ii+ntclen)
ii=ii+ntclen
end if
if (ii/=nclen) then
write(6,*)'allocate_cv: error length',ii,nclen
call stop2(109)
end if
! Construct a list of integer pointers to the static and motley part of the
! control vector (this is to support operations on current (mpi) distribution)
! it's enough to get these pointers to nsubwin=1, since they only serve the
! purpose of indexing arrays operating on single window (e.g. see grid2sub)
! allocate(ycv%ivalues(n3d+n2d+mvars))
! ii=0
! do i=1,nc3d
! call gsi_bundlegetpointer (ycv%step(1),cvars3d,iptr,istatus,ival=ival)
! if(istatus==0) then
! ii=ii+1
! ycv%ivalues(ii)=ival
! endif
! enddo
! do i=1,nc2d
! call gsi_bundlegetpointer (ycv%step(1),cvars2d,iptr,istatus,ival=ival)
! if(istatus==0) then
! ii=ii+1
! ycv%ivalues(ii)=ival
! endif
! enddo
! do i=1,mvars
! call gsi_bundlegetpointer (ycv%motley(1),cvarsmd,iptr,istatus,ival=ival)
! if(istatus==0) then
! ii=ii+1
! ycv%ivalues(ii)=ival
! endif
! enddo
m_allocs=m_allocs+1
m_vec_alloc=m_vec_alloc+1
max_vec_alloc=MAX(max_vec_alloc,m_vec_alloc)
return
end subroutine allocate_cv
! ----------------------------------------------------------------------
subroutine deallocate_cv(ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: deallocate_cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
! 2009-09-20 parrish - add optional removal of pointer to hybrid ensemble control variable a_en
! 2010-02-25 zhu - add flexibility to control variable
! 2010-05-01 todling - update to use gsi_bundle
! 2010-05-17 todling - add back ens control; w/ a twist from original
! 2010-05-22 todling - add support for motley variables
! 2011-09-05 todling - replace destroy calls with unset calls
!
! input argument list:
! ycv
!
! output argument list:
! ycv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(inout) :: ycv
integer(i_kind) :: ii,nn,ierror
if (ycv%lallocated) then
do ii=1,nsubwin
if (l_hyb_ens) then
do nn=n_ens,1,-1
call GSI_BundleUnset(ycv%aens(ii,nn),ierror)
enddo
endif
! if (beta_s0>tiny_r_kind) then
if(mvars>0) then
call GSI_BundleDestroy(ycv%motley(ii),ierror)
endif
call GSI_BundleUnset(ycv%step(ii),ierror)
! endif ! beta_s0
end do
NULLIFY(ycv%predr)
NULLIFY(ycv%predp)
NULLIFY(ycv%predt)
if(l_hyb_ens) DEALLOCATE(ycv%aens)
if(mvars>0) DEALLOCATE(ycv%motley)
DEALLOCATE(ycv%step)
DEALLOCATE(ycv%values)
ycv%lallocated=.false.
m_deallocs=m_deallocs+1
m_vec_alloc=m_vec_alloc-1
else
if (mype==0) write(6,*)'deallocate_cv warning: vector not allocated'
endif
return
end subroutine deallocate_cv
! ----------------------------------------------------------------------
subroutine assign_scalar2cv(ycv,pval)
!$$$ subprogram documentation block
! . . . .
! subprogram: assign_scalar2cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! ycv
! pval
!
! output argument list:
! ycv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(inout) :: ycv
real(r_kind) , intent(in ) :: pval
integer(i_kind) :: ii
DO ii=1,ycv%lencv
ycv%values(ii)=pval
ENDDO
return
end subroutine assign_scalar2cv
! ----------------------------------------------------------------------
subroutine assign_cv2cv(ycv,xcv)
!$$$ subprogram documentation block
! . . . .
! subprogram: assign_cv2cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! ycv
! xcv
!
! output argument list:
! ycv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(inout) :: ycv
type(control_vector), intent(in ) :: xcv
integer(i_kind) :: ii
if (xcv%lencv/=ycv%lencv) then
write(6,*)'assign_cv2cv: error length',xcv%lencv,ycv%lencv
call stop2(110)
end if
DO ii=1,ycv%lencv
ycv%values(ii)=xcv%values(ii)
ENDDO
return
end subroutine assign_cv2cv
! ----------------------------------------------------------------------
subroutine assign_array2cv(ycv,parray)
!$$$ subprogram documentation block
! . . . .
! subprogram: assign_array2cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! ycv
! parray
!
! output argument list:
! ycv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(inout) :: ycv
real(r_kind) , intent(in ) :: parray(:)
integer(i_kind) :: ii
if (size(parray)/=ycv%lencv) then
write(6,*)'assign_array2cv: array wrong length',size(parray),ycv%lencv
call stop2(111)
end if
DO ii=1,ycv%lencv
ycv%values(ii)=parray(ii)
ENDDO
return
end subroutine assign_array2cv
! ----------------------------------------------------------------------
subroutine assign_cv2array(parray,ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: assign_array2cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! ycv
!
! output argument list:
! parray
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
real(r_kind) , intent( out) :: parray(:)
type(control_vector), intent(in ) :: ycv
integer(i_kind) :: ii
if (size(parray)/=ycv%lencv) then
write(6,*)'assign_cv2array: array wrong length',size(parray),ycv%lencv
call stop2(112)
end if
DO ii=1,ycv%lencv
parray(ii)=ycv%values(ii)
ENDDO
return
end subroutine assign_cv2array
! ----------------------------------------------------------------------
real(r_quad) function qdot_prod_sub(xcv,ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: qdot_prod_sub
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
! 2009-09-20 parrish - add hybrid ensemble control variable a_en contribution to dot product
! 2010-05-17 todling - update to use bundle
! 2010-06-02 parrish - add contribution from ensemble control variable to dot product
! 2010-10-08 derber - optimize and clean up
! 2011-09-05 todling - add hybrid to sqrt part of code
!
! input argument list:
! xcv,ycv
!
! output argument list:
! prods
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(in ) :: xcv, ycv
integer(i_kind) :: ii,nn,m3d,m2d,i,j,itot
real(r_quad),allocatable,dimension(:) :: partsum
qdot_prod_sub=zero_quad
! Independent part of vector
if (lsqrtb) then
do ii=1,nsubwin
qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%step(ii)%values(:) ,ycv%step(ii)%values(:) )
end do
if(l_hyb_ens) then
do nn=1,n_ens
do ii=1,nsubwin
qdot_prod_sub=qdot_prod_sub+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) )
end do
end do
endif
else
do ii=1,nsubwin
m3d=xcv%step(ii)%n3d
m2d=xcv%step(ii)%n2d
itot=max(m3d,0)+max(m2d,0)
if(l_hyb_ens)itot=itot+n_ens
allocate(partsum(itot))
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,m3d
partsum(i) = dplevs(xcv%step(ii)%r3(i)%q,ycv%step(ii)%r3(i)%q,ihalo=1)
enddo
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,m2d
partsum(m3d+i) = dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1)
enddo
if(l_hyb_ens) then
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,n_ens
partsum(m3d+m2d+i) = dplevs(xcv%aens(ii,i)%r3(1)%q,ycv%aens(ii,i)%r3(1)%q,ihalo=1)
end do
end if
do i=1,itot
qdot_prod_sub = qdot_prod_sub + partsum(i)
end do
deallocate(partsum)
end do
end if
! Duplicated part of vector
if(mype == 0)then
do j=nclen1+1,nclen
qdot_prod_sub=qdot_prod_sub+xcv%values(j)*ycv%values(j)
end do
end if
return
end function qdot_prod_sub
! ----------------------------------------------------------------------
subroutine qdot_prod_vars_eb(xcv,ycv,prods,eb)
!$$$ subprogram documentation block
! . . . .
! subprogram: qdot_prod_vars_eb copy of qdot_prod_vars for J_ens
! prgmmr: parrish org: np22 date: 2009-09-20
!
! abstract: copy qdot_prod_vars and add extra code to compute J_b or
! J_ens when running in hybrid_ensemble mode. extra input
! character string eb is used to signal if J_b or J_ens is
! to be computed.
!
! program history log:
! 2009-09-20 parrish - initial documentation
! 2010-05-17 todling - update to use bundle
! 2010-10-08 derber - optimize and clean up
! 2011-09-05 todling - add hybrid to sqrt part of code
!
! input argument list:
! xcv,ycv
! eb - eb= 'cost_b' then return J_b in prods
! = 'cost_e' then return J_ens in prods
!
! output argument list:
! prods
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(in ) :: xcv, ycv
character(len=*) , intent(in ) :: eb
real(r_quad) , intent( out) :: prods(nsubwin+1)
real(r_quad) :: zz(nsubwin)
integer(i_kind) :: ii,i,nn,m3d,m2d
real(r_quad),allocatable,dimension(:) :: partsum
prods(:)=zero_quad
zz(:)=zero_quad
! Independent part of vector
if (lsqrtb) then
if(trim(eb) == 'cost_b') then
do ii=1,nsubwin
zz(ii)=zz(ii)+qdot_product( xcv%step(ii)%values(:) ,ycv%step(ii)%values(:) )
end do
endif
if(trim(eb) == 'cost_e') then
do nn=1,n_ens
do ii=1,nsubwin
zz(ii)=zz(ii)+qdot_product( xcv%aens(ii,nn)%values(:) ,ycv%aens(ii,nn)%values(:) )
end do
end do
endif
else
if(trim(eb) == 'cost_b') then
do ii=1,nsubwin
m3d=xcv%step(ii)%n3d
m2d=xcv%step(ii)%n2d
allocate(partsum(m2d+m3d))
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,m3d
partsum(i)= dplevs(xcv%step(ii)%r3(i)%q,ycv%step(ii)%r3(i)%q,ihalo=1)
enddo
!$omp parallel do schedule(dynamic,1) private(i)
do i = 1,m2d
partsum(m3d+i)= dplevs(xcv%step(ii)%r2(i)%q,ycv%step(ii)%r2(i)%q,ihalo=1)
enddo
do i = 1,m2d+m3d
zz(ii)=zz(ii) + partsum(i)
end do
deallocate(partsum)
end do
end if
if(trim(eb) == 'cost_e') then
do ii=1,nsubwin ! RTod: somebody could work in opt/zing this ...
allocate(partsum(n_ens))
!$omp parallel do schedule(dynamic,1) private(nn,m3d,m2d)
do nn=1,n_ens
partsum(nn) = zero_quad
m3d=xcv%aens(ii,nn)%n3d
do i = 1,m3d
partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r3(i)%q,ycv%aens(ii,nn)%r3(i)%q,ihalo=1)
enddo
m2d=xcv%aens(ii,nn)%n2d
do i = 1,m2d
partsum(nn)= partsum(nn) + dplevs(xcv%aens(ii,nn)%r2(i)%q,ycv%aens(ii,nn)%r2(i)%q,ihalo=1)
enddo
enddo
do nn=1,n_ens
zz(ii)=zz(ii)+partsum(nn)
end do
deallocate(partsum)
end do
end if
end if
call mpl_allreduce(nsubwin,qpvals=zz)
prods(1:nsubwin) = zz(1:nsubwin)
! Duplicated part of vector
if(trim(eb) == 'cost_b') then
if (nsclen>0) then
prods(nsubwin+1) = prods(nsubwin+1) + qdot_product(xcv%predr(:),ycv%predr(:))
endif
if (npclen>0) then
prods(nsubwin+1) = prods(nsubwin+1) + qdot_product(xcv%predp(:),ycv%predp(:))
endif
if (ntclen>0) then
prods(nsubwin+1) = prods(nsubwin+1) + qdot_product(xcv%predt(:),ycv%predt(:))
endif
end if
return
end subroutine qdot_prod_vars_eb
! ----------------------------------------------------------------------
real(r_kind) function dot_prod_cv(xcv,ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: dot_prod_cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! xcv,ycv
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(in ) :: xcv, ycv
! local variables
real(r_quad) :: dd(1)
if (xcv%lencv/=ycv%lencv) then
write(6,*)'dot_prod_cv: error length',xcv%lencv,ycv%lencv
call stop2(113)
end if
dd(1) = qdot_prod_sub(xcv,ycv)
call mpl_allreduce(1,qpvals=dd)
dot_prod_cv = dd(1)
return
end function dot_prod_cv
! ----------------------------------------------------------------------
real(r_quad) function qdot_prod_cv(xcv,ycv,mold)
!$$$ subprogram documentation block
! . . . .
! subprogram: qdot_prod_cv
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
! 2017-03-06 todling - rename interface variable to mold
!
! input argument list:
! xcv,ycv
! mold - interface device
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind) , intent(in ) :: mold
type(control_vector), intent(in ) :: xcv, ycv
! local variables
real(r_quad) :: dd(1)
if (xcv%lencv/=ycv%lencv) then
write(6,*)'qdot_prod_cv: error length',xcv%lencv,ycv%lencv
call stop2(114)
end if
dd(1) = qdot_prod_sub(xcv,ycv)
call mpl_allreduce(1,qpvals=dd)
qdot_prod_cv = dd(1)
return
end function qdot_prod_cv
! ----------------------------------------------------------------------
real(r_quad) function qdot_prod_cv_eb(xcv,ycv,mold,eb)
!$$$ subprogram documentation block
! . . . .
! subprogram: qdot_prod_cv_eb copy of qdot_prod_cv for J_ens
! prgmmr: parrish org: np22 date: 2009-09-20
!
! abstract: copy qdot_prod_cv and add extra code to compute J_b or
! J_ens when running in hybrid_ensemble mode. extra input
! character string eb is used to signal if J_b or J_ens is
! to be computed.
!
! program history log:
! 2009-09-20 parrish - initial documentation
! 2017-03-06 todling - rename interface variable to mold
!
! input argument list:
! xcv,ycv
! mold - interface device
! eb - eb= 'cost_b' then return J_b in prods
! = 'cost_e' then return J_ens in prods
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
integer(i_kind) , intent(in ) :: mold
character(len=*) , intent(in ) :: eb
type(control_vector), intent(in ) :: xcv, ycv
! local variables
real(r_quad) :: zz(nsubwin+1)
integer(i_kind) :: ii
if (xcv%lencv/=ycv%lencv) then
write(6,*)'qdot_prod_cv_eb: error length',xcv%lencv,ycv%lencv
call stop2(114)
end if
call qdot_prod_vars_eb(xcv,ycv,zz,eb)
qdot_prod_cv_eb= zero_quad
do ii=1,nsubwin+1
qdot_prod_cv_eb = qdot_prod_cv_eb + zz(ii)
enddo
return
end function qdot_prod_cv_eb
! ----------------------------------------------------------------------
subroutine prt_norms(xcv,sgrep)
!$$$ subprogram documentation block
! . . . .
! subprogram: prt_norms
! prgmmr: org: date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! xcv
! sgrep
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(in ) :: xcv
character(len=*) , intent(in ) :: sgrep
real(r_quad) :: zt
zt = qdot_prod_cv(xcv,xcv,r_quad)
zt=sqrt(zt)
if (mype==0) then
write(6,*)sgrep,' global norm =',real(zt,r_kind)
endif
!_RT call prt_norms_vars(xcv,sgrep) --->> this routine is hanging
return
end subroutine prt_norms
! ----------------------------------------------------------------------
subroutine prt_norms_vars(xcv,sgrep)
!$$$ subprogram documentation block
! . . . .
! subprogram: prt_norms_vars
! prgmmr: Jing Guo org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
! 2010-05-06 todling- update to use gsi_bundle-like vector
!
! input argument list:
! xcv
! sgrep
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
use m_stats,only : stats_sum,stats_allreduce
implicit none
type(control_vector), intent(in ) :: xcv
character(len=*) , intent(in ) :: sgrep
real (r_kind),allocatable,dimension(:) :: vdot,vsum,vmin,vmax
integer(i_kind),allocatable,dimension(:) :: vnum
integer :: jj,nsw,iv,nv
real(r_kind),pointer,dimension(:) :: piv
nsw=size(xcv%step)
! process 3d variables first
allocate(vdot(nc3d),vsum(nc3d),vmin(nc3d),vmax(nc3d),vnum(nc3d))
do iv=1,nc3d
do jj=1,nsw
piv => rerank(xcv%step(jj)%r3(iv)%q)
call stats_sum(piv, &
vdot(iv),vsum(iv),vmin(iv),vmax(iv),vnum(iv),add=jj>1)
enddo
call stats_allreduce(vdot(iv),vsum(iv),vmin(iv),vmax(iv), &
vnum(iv),MPI_comm_world)
nv=max(vnum(iv),1)
if(mype==0) then
write(6,'(2(1x,a),4(1x,ES20.12),1x,i10)') &
sgrep,cvars3d(iv),sqrt(vdot(iv)/nv),vsum(iv)/nv, &
vmin(iv),vmax(iv),vnum(iv)
endif
end do
deallocate(vdot,vsum,vmin,vmax,vnum)
! process 2d now
allocate(vdot(nc2d),vsum(nc2d),vmin(nc2d),vmax(nc2d),vnum(nc2d))
do iv=1,nc2d
do jj=1,nsw
piv => rerank(xcv%step(jj)%r2(iv)%q)
call stats_sum(piv, &
vdot(iv),vsum(iv),vmin(iv),vmax(iv),vnum(iv),add=jj>1)
enddo
call stats_allreduce(vdot(iv),vsum(iv),vmin(iv),vmax(iv), &
vnum(iv),MPI_comm_world)
nv=max(vnum(iv),1)
if(mype==0) then
write(6,'(2(1x,a),4(1x,ES20.12),1x,i10)') &
sgrep,cvars2d(iv),sqrt(vdot(iv)/nv),vsum(iv)/nv, &
vmin(iv),vmax(iv),vnum(iv)
endif
end do
deallocate(vdot,vsum,vmin,vmax,vnum)
! release pointer
piv => null()
end subroutine prt_norms_vars
! ----------------------------------------------------------------------
subroutine axpy(alpha,xcv,ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: axpy
! prgmmr: todling org: gmao date:
!
! abstract: similar to BLAS axpy
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! xcv
! alpha
!
! output argument list:
! ycv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
real(r_kind) , intent(in ) :: alpha
type(control_vector), intent(in ) :: xcv
type(control_vector), intent(inout) :: ycv
integer(i_kind) :: ii
if (xcv%lencv/=ycv%lencv) then
write(6,*)'axpy: error length',xcv%lencv,ycv%lencv
call stop2(115)
end if
DO ii=1,ycv%lencv
ycv%values(ii) = ycv%values(ii) + alpha * xcv%values(ii)
ENDDO
return
end subroutine axpy
! ----------------------------------------------------------------------
subroutine random_cv(ycv,kseed)
!$$$ subprogram documentation block
! . . . .
! subprogram: random_cv
! prgmmr: tremolet org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! ycv
! kseed
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector) , intent(inout) :: ycv
integer(i_kind), optional, intent(in ) :: kseed
integer(i_kind):: ii,jj,nn,iseed
integer, allocatable :: nseed(:) ! Intentionaly default integer
real(r_kind), allocatable :: zz(:)
iseed=iadatebgn
if (present(kseed)) iseed=iseed+kseed
call random_seed(size=jj)
allocate(nseed(jj))
nseed(1:jj)=iseed
! The following because we don't want all procs to get
! exactly the same sequence (which would be repeated in
! the then not so random vector) but it makes the test
! not reproducible if the number of procs is changed.
nseed(1)=iseed+mype
call random_seed(put=nseed)
deallocate(nseed)
allocate(zz(nval_len))
do jj=1,nsubwin
call random_number(zz)
do ii=1,nval_len
ycv%step(jj)%values(ii) = two*zz(ii)-one
enddo
enddo
deallocate(zz)
if (nval_lenz_en>0) then
allocate(zz(nval_lenz_en))
do nn=1,n_ens
do jj=1,nsubwin
call random_number(zz)
do ii=1,nval_lenz_en
ycv%aens(jj,nn)%values(ii) = two*zz(ii)-one
enddo
enddo
enddo
deallocate(zz)
endif
if (nsclen>0) then
allocate(zz(nsclen))
call random_number(zz)
do ii=1,nsclen
ycv%predr(ii) = two*zz(ii)-one
enddo
deallocate(zz)
endif
if (npclen>0) then
allocate(zz(npclen))
call random_number(zz)
do ii=1,npclen
ycv%predp(ii) = two*zz(ii)-one
enddo
deallocate(zz)
endif
if (ntclen>0) then
allocate(zz(ntclen))
call random_number(zz)
do ii=1,ntclen
ycv%predt(ii) = two*zz(ii)-one
enddo
deallocate(zz)
endif
return
end subroutine random_cv
! ----------------------------------------------------------------------
subroutine write_cv(xcv,cdfile)
!$$$ subprogram documentation block
! . . . .
! subprogram: write_cv
! prgmmr: tremolet org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! xcv
! cdfile
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(in ) :: xcv
character(len=*) , intent(in ) :: cdfile
character(len=100) :: clfile
character(len=5) :: clmype
integer(i_kind):: iunit
iunit=get_lun()
clmype='.YYYY'
write(clmype(2:5),'(I4.4)')mype
clfile=trim(cdfile)//clmype
if (mype==0) write(6,*)'Writing control vector to file ',clfile
open(iunit,file=trim(clfile),form='unformatted')
write(iunit)xcv%lencv
write(iunit)xcv%values(1:xcv%lencv)
close(iunit)
return
end subroutine write_cv
! ----------------------------------------------------------------------
subroutine read_cv(xcv,cdfile)
!$$$ subprogram documentation block
! . . . .
! subprogram: read_cv
! prgmmr: tremolet org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! xcv
! cdfile
!
! output argument list:
! xcv
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(inout) :: xcv
character(len=*) , intent(in ) :: cdfile
character(len=100) :: clfile
character(len=5) :: clmype
integer(i_kind):: iunit,ilen
iunit=get_lun()
clmype='.YYYY'
write(clmype(2:5),'(I4.4)')mype
clfile=trim(cdfile)//clmype
if (mype==0) write(6,*)'Reading control vector from file ',clfile
open(iunit,file=trim(clfile),form='unformatted')
read(iunit)ilen
if (ilen/=xcv%lencv) then
write(6,*)'read_cv: wrong length',ilen,xcv%lencv
call stop2(116)
end if
read(iunit)xcv%values(1:xcv%lencv)
close(iunit)
return
end subroutine read_cv
! ----------------------------------------------------------------------
subroutine inquire_cv
!$$$ subprogram documentation block
! . . . .
! subprogram: inquire_cv
! prgmmr: tremolet org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
real(r_kind) :: zz
if (mype==0) then
zz=real(max_vec_alloc*nclen,r_kind)*8.0_r_kind/1.048e6_r_kind
write(6,*)'control_vectors: length=',nclen
write(6,*)'control_vectors: currently allocated=',m_vec_alloc
write(6,*)'control_vectors: maximum allocated=',max_vec_alloc
write(6,*)'control_vectors: number of allocates=',m_allocs
write(6,*)'control_vectors: number of deallocates=',m_deallocs
write(6,'(A,F8.1,A)')'control_vectors: Estimated max memory used= ',zz,' Mb'
endif
end subroutine inquire_cv
! ----------------------------------------------------------------------
real(r_kind) function maxval_cv(ycv)
!$$$ subprogram documentation block
! . . . .
! subprogram: maxval_cv
! prgmmr: tremolet org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! ycv
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
type(control_vector), intent(in ) :: ycv
real(r_kind) :: zloc(1),zglo(1)
zloc(1)=maxval(ycv%values(:))
call mpi_allreduce(zloc,zglo,1,mpi_rtype,mpi_max,mpi_comm_world,ierror)
if (ierror/=0) then
write(6,*)'maxval_cv: MPI error',ierror
call stop2(117)
end if
maxval_cv=zglo(1)
return
end function maxval_cv
! ----------------------------------------------------------------------
real(r_quad) function qdot_product(x,y)
!$$$ subprogram documentation block
! . . . .
! subprogram: qdot_product
! prgmmr: todling org: gmao date:
!
! abstract:
!
! program history log:
! 2009-08-04 lueken - added subprogram doc block
!
! input argument list:
! x,y
!
! output argument list:
!
! attributes:
! language: f90
! machine:
!
!$$$ end documentation block
implicit none
real(r_kind),intent(in ) :: x(:),y(:)
real(r_quad):: zz
integer(i_kind) :: nx,ny,i
nx=size(x)
ny=size(y)
if(nx/=ny) then
write(6,*)'qdot_product: inconsistent dims',nx,ny
call stop2(118)
end if
zz=zero_quad
do i=1,nx
zz = zz + x(i)*y(i)
enddo
qdot_product=zz
end function qdot_product
! ----------------------------------------------------------------------
end module control_vectors