#848 |
|
arango
|
Fixed
|
VERY IMPORTANT: Corrected bugs in tracer advection
|
Description |
Corrected critical bug in some of the tracer advection options:
- Corrected a bug in the HSIMT advection. Previously, I fixed the issue that three ghost points are needed for this scheme. At that time, I missed making a similar change in the computation of the inverse thickness (oHz). Therefore, in step3d_t.F we need to have
logical :: LapplySrc, Lhsimt, Lmpdata
...
Lhsimt =ANY(Hadvection(:,ng)%HSIMT).and. &
& ANY(Vadvection(:,ng)%HSIMT)
Lmpdata=ANY(Hadvection(:,ng)%MPDATA).and. &
& ANY(Vadvection(:,ng)%MPDATA)
...
!
! Compute reciprocal thickness, 1/Hz.
!
IF (Lmpdata.or.Lhsimt) THEN
DO k=1,N(ng)
DO j=Jstrm2,Jendp2
DO i=Istrm2,Iendp2
oHz(i,j,k)=1.0_r8/Hz(i,j,k)
END DO
END DO
END DO
ELSE
DO k=1,N(ng)
DO j=Jstr,Jend
DO i=Istr,Iend
oHz(i,j,k)=1.0_r8/Hz(i,j,k)
END DO
END DO
END DO
END IF
I have been hunting on and off for a bug in the HSIMT scheme for awhile. Many thanks to Pierre St-Laurent for reporting this issue. It is always helpful when another set of eyes look into a code.
- Corrected a weird bug in the tangent linear advection switches. In inp_par.F we need to have instead:
#if defined TANGENT || defined TL_IOMS
!
! Set tracer advection scheme switches for the tangent linear models
! (TLM and RPM) to the same values as the adjoint model.
!
DO ng=1,Ngrids
DO i=1,NT(ng)
tl_Hadvection(i,ng)%AKIMA4 = ad_Hadvection(i,ng)%AKIMA4
tl_Hadvection(i,ng)%CENTERED2 = ad_Hadvection(i,ng)%CENTERED2
tl_Hadvection(i,ng)%CENTERED4 = ad_Hadvection(i,ng)%CENTERED4
tl_Hadvection(i,ng)%HSIMT = ad_Hadvection(i,ng)%HSIMT
tl_Hadvection(i,ng)%MPDATA = ad_Hadvection(i,ng)%MPDATA
tl_Hadvection(i,ng)%SPLINES = ad_Hadvection(i,ng)%SPLINES
tl_Hadvection(i,ng)%SPLIT_U3 = ad_Hadvection(i,ng)%SPLIT_U3
tl_Hadvection(i,ng)%UPSTREAM3 = ad_Hadvection(i,ng)%UPSTREAM3
!
tl_Vadvection(i,ng)%AKIMA4 = ad_Vadvection(i,ng)%AKIMA4
tl_Vadvection(i,ng)%CENTERED2 = ad_Vadvection(i,ng)%CENTERED2
tl_Vadvection(i,ng)%CENTERED4 = ad_Vadvection(i,ng)%CENTERED4
tl_Vadvection(i,ng)%HSIMT = ad_Vadvection(i,ng)%HSIMT
tl_Vadvection(i,ng)%MPDATA = ad_Vadvection(i,ng)%MPDATA
tl_Vadvection(i,ng)%SPLINES = ad_Vadvection(i,ng)%SPLINES
tl_Vadvection(i,ng)%SPLIT_U3 = ad_Vadvection(i,ng)%SPLIT_U3
tl_Vadvection(i,ng)%UPSTREAM3 = ad_Vadvection(i,ng)%UPSTREAM3
END DO
END DO
#endif
Before, the assigment was to use the nonlinear values instead. As a consequence, the 4D-Var cost function was not decreasing correctly and the minimization failed because the broken symmetry between tangent linear (TLM) and adjoint (ADM) model operators. We need to have the same advection scheme in both TLM and ADM. Recall that in roms.in, the user specify just the values for the adjoint model.
Many thanks to Julia Levin for reporting this problem and Andy Moore for his help in debugging.
|
#849 |
|
arango
|
Done
|
Miscellaneous Update
|
Description |
This update includes a couple of changes:
- Fixed the tracer advection logic in inp_par.F and tadv.F for the shallow-water (SOLVE3D undefined) applications not needing tracer advection.
!
! Set switch for three ghost-points in the halo region.
!
#ifdef SOLVE3D
ThreeGhostPoints=ANY(Hadvection(:,:)%MPDATA).or. &
& ANY(Hadvection(:,:)%HSIMT)
#endif
#ifdef UV_VIS4
ThreeGhostPoints=.TRUE.
#endif
Many thanks to Guangyu Xu for reporting this issue.
- Updated inp_decode.F function load_s1d to have the following interface:
INTERFACE load_s1d
MODULE PROCEDURE load_s1d1 ! 1D structrure, S(:)
MODULE PROCEDURE load_s1d2 ! 2D structrure, S(Ie,:)
END INTERFACE load_s1d
so we have specific funtions for 1D (load_s1d1) and 2D (load_s1d2) TYPE_IO structures. For example, in read_asspar.F we have:
CASE ('STDnameI')
label='STD - initial conditions standard deviation'
Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid, &
& Ngrids, Nfiles, 4, 1, STD)
CASE ('STDnameM')
label='STD - model error standard deviation'
Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid, &
& Ngrids, Nfiles, 4, 2, STD)
CASE ('STDnameB')
label='STD - boundary conditions standard deviation'
Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid, &
& Ngrids, Nfiles, 4, 3, STD)
CASE ('STDnameF')
label='STD - surface forcing standard deviation'
Npts=load_s1d(Nval, Cval, Cdim, line, label, igrid, &
& Ngrids, Nfiles, 4, 4, STD)
to process each element of the STD(1:4,Ngrids) standard deviation structure used in 4D-Var.
|
#851 |
|
arango
|
Done
|
VERY IMPORTANT: Consolidating Research Repositories with Released Version
|
Description |
I have several research repositories for several ongoing new algorithms that will be available in the future. I have so many exceptional flags to update the codes that sometimes one fails, and undesired changes appear in public released version of ROMS and erase features in the research versions. To avoid doing that in the future, I am consolidating the research 4DVAR repository with the public version of the code.
For the last two months, we have been developing the saddle-point 4D-Var. By far, one of the most technically tricky algorithms that Andy and I have put our hands on. As formulated, the current available 4D-Var algorithms are sequential, and its inner loops cannot be parallelized in time to accelerate the analysis. The conjugate gradient solver used to minimize the cost function J, depends sequentially on the previous iteration.
Fisher and Gürol (2017) proposed a new weak constraint 4D-Var algorithm based on the saddle-point formulation that can be parallelized in time. The time trajectory can be split into the solution of Nsaddle intervals that can run concurrently to reduce the 4D-Var computational cost.
To achieve these technically challenging computations, we need to make changes to several ROMS routines to:
- To split the MPI communicator to carry out disjointed and concurrent runs of the nonlinear, tangent linear, and adjoint ROMS kernels. Therefore, several of the routines of distribute.F have an optional argument to pass the desired MPI communicators that are defined in mod_parallel.F.
- Add optional arguments to output NetCDF routines nf_fwrite2d.F, nf_fwrite2d_bry.F, nf_fwrite3d.F, nf_fwrite3d_bry.F, and nf_fwrite4d.F to get the minimun and maximun values ot the written fields.
- Include several new generic routines to manipulate ROMS state vector: ADfromTL.F, def_state.F, state_join.F, state_read.F, and wrt_state.F.
- No change is made to ROMS nonlinear, tangent linear, and adjoint time-stepping kernels.
- Added several routines from the LAPACK library that are needed for the Generalized Minimal Residual (GMRES) minimization solver.
- The saddle-point algorithm with be released in the future when it is fully tested and we have the proper documentation and relevant publication.
|