MODULE rbl4dvar_mod ! !git $Id$ !svn $Id: rbl4dvar.F 1151 2023-02-09 03:08:53Z arango $ !================================================== Hernan G. Arango === ! Copyright (c) 2002-2023 The ROMS/TOMS Group Andrew M. Moore ! ! Licensed under a MIT/X style license ! ! See License_ROMS.md ! !======================================================================= ! ! ! This module splits the RBL4D-Var data assimilation algorithm into ! ! its logical components routines. ! ! ! ! background_initialize: ! ! ! ! Initializes the nonlinear model kernel. It is separated from ! ! the 'background' phase to allow ESM coupling. ! ! ! ! background: ! ! ! ! Timesteps the nonlinear model to compute the basic state Xb(t) ! ! used to linearize the tangent linear and adjoint models. It ! ! interpolates the nonlinear model trajectory to the observations ! ! locations. ! ! ! ! increment: ! ! ! ! Minimizes of the 4D-Var cost function over Ninner inner loops ! ! iterations to compute the data assimilation increment, dXa. It ! ! also computes the nonlinear model initial conditions by adding ! ! the increment to the background, Xa = Xb(t=0) + dXa. The new ! ! NLM initial conditions is written inot INI(ng)%name NetCDF file.! ! ! ! analysis_initialize: ! ! ! ! Initializes the nonlinear model with the 4D-Var analysis, Xa. ! ! It is separated from 'analysis' phase to allow ESM coupling. ! ! ! ! analysis: ! ! ! ! Timesteps the nonlinear model to compute new state trajectory ! ! that becomes the basic state for the next outer loop. It also ! ! interpolates the solution at the observation locations. ! ! ! ! prior_error: ! ! ! ! Processes the background and model prior error covariance and ! ! its normalization coefficients. ! ! ! ! posterior_error: ! ! ! ! Computes the posterior error analysis. ! ! ! ! ! ! References: ! ! ! ! Moore, A.M., H.G. Arango, G. Broquet, B.S. Powell, A.T. Weaver, ! ! and J. Zavala-Garay, 2011: The Regional Ocean Modeling System ! ! (ROMS) 4-dimensional variational data assimilations systems, ! ! Part I - System overview and formulation, Prog. Oceanogr., 91, ! ! 34-49, doi:10.1016/j.pocean.2011.05.004. ! ! ! ! Moore, A.M., H.G. Arango, G. Broquet, C. Edward, M. Veneziani, ! ! B. Powell, D. Foley, J.D. Doyle, D. Costa, and P. Robinson, ! ! 2011: The Regional Ocean Modeling System (ROMS) 4-dimensional ! ! variational data assimilations systems, Part II - Performance ! ! and application to the California Current System, Prog. ! ! Oceanogr., 91, 50-73, doi:10.1016/j.pocean.2011.05.003. ! ! ! ! Gurol, S., A.T. Weaver, A.M. Moore, A. Piacentini, H.G. Arango, ! ! S. Gratton, 2014: B-preconditioned minimization algorithms for ! ! data assimilation with the dual formulation, QJRMS, 140, ! ! 539-556. ! ! ! !======================================================================= ! USE mod_kinds USE mod_param USE mod_parallel USE mod_fourdvar USE mod_iounits USE mod_ncparam USE mod_netcdf USE mod_scalars USE mod_stepping ! USE mod_forces, ONLY : initialize_forces USE mod_mixing, ONLY : initialize_mixing USE mod_ocean, ONLY : initialize_ocean ! USE ad_wrt_his_mod, ONLY : ad_wrt_his USE convolve_mod, ONLY : error_covariance USE def_impulse_mod, ONLY : def_impulse USE def_mod_mod, ONLY : def_mod USE def_norm_mod, ONLY : def_norm USE get_state_mod, ONLY : get_state USE normalization_mod, ONLY : normalization USE rpcg_lanczos_mod, ONLY : rpcg_lanczos USE strings_mod, ONLY : FoundError USE tl_def_ini_mod, ONLY : tl_def_ini USE tl_wrt_ini_mod, ONLY : tl_wrt_ini USE wrt_aug_imp_mod, ONLY : wrt_aug_imp USE wrt_impulse_mod, ONLY : wrt_impulse USE wrt_ini_mod, ONLY : wrt_ini ! implicit none ! PUBLIC :: background_initialize PUBLIC :: background PUBLIC :: increment PUBLIC :: analysis_initialize PUBLIC :: analysis PUBLIC :: prior_error ! ! Set module internal parameters. ! integer :: Lbck = 2 integer :: Lini = 1 integer :: LTLM1 = 1 integer :: LTLM2 = 2 integer :: Rec1 = 1 integer :: Rec2 = 2 integer :: Rec3 = 3 integer :: Rec4 = 4 integer :: Rec5 = 5 ! CONTAINS ! SUBROUTINE background_initialize (my_outer) ! !======================================================================= ! ! ! This routine initializes ROMS nonlinear model trajectory Xb_n-1(0) ! ! for each (n-1) outer loops. It is separated from the 'background' ! ! 4D-Var phase in ESM coupling applications that use generic methods ! ! for 'initialize', 'run', and 'finalize'. ! ! ! ! On Input: ! ! ! ! my_outer Outer-loop counter (integer) ! ! ! !======================================================================= ! ! Imported variable declarations ! integer, intent(in) :: my_outer ! ! Local variable declarations. ! integer :: lstr, ng, tile integer :: Fcount, Tindex integer :: thread ! character (len=*), parameter :: MyFile = & & "ROMS/Drivers/rbl4dvar.F"//", background_initialize" ! SourceFile=MyFile ! !----------------------------------------------------------------------- ! Initialize nonlinear model kernel. !----------------------------------------------------------------------- ! ! Start profile clock. ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_on (ng, iNLM, 86, 242, MyFile) END DO END DO ! ! Initialize the switch to gather weak constraint forcing. ! DO ng=1,Ngrids WRTforce(ng)=.FALSE. END DO ! ! Clear TLM mixing arrays. ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL initialize_mixing (ng, tile, iTLM) END DO END DO ! ! Get nonlinear model initial conditions. ! DO ng=1,Ngrids wrtNLmod(ng)=.TRUE. wrtTLmod(ng)=.FALSE. LreadBLK(ng)=.FALSE. LreadFWD(ng)=.FALSE. RST(ng)%Rindex=0 Fcount=RST(ng)%load RST(ng)%Nrec(Fcount)=0 END DO ! CALL initial IF (FoundError(exit_flag, NoError, 279, MyFile)) RETURN ! ! Save nonlinear initial conditions (currently in time index 1, ! background) into record "Lbck" of INI(ng)%name NetCDF file. The ! record "Lbck" becomes the background state and "Lini" becomes ! the current nonlinear initial conditions. ! Tindex=1 DO ng=1,Ngrids INI(ng)%Rindex=1 Fcount=INI(ng)%load INI(ng)%Nrec(Fcount)=1 CALL wrt_ini (ng, MyRank, Tindex, Lbck) IF (FoundError(exit_flag, NoError, 296, MyFile)) RETURN END DO ! ! Set nonlinear output history file as the initial basic state ! trajectory. ! DO ng=1,Ngrids LdefHIS(ng)=.TRUE. LwrtHIS(ng)=.TRUE. WRITE (HIS(ng)%name,10) TRIM(FWD(ng)%head), my_outer lstr=LEN_TRIM(HIS(ng)%name) HIS(ng)%base=HIS(ng)%name(1:lstr-3) END DO ! ! Set the nonlinear model to output the quicksave history file as a ! function of the outer loop. It may be used as the basic state ! trajectory for the surface fluxes (wind stress, shortwave, heat ! flux, and E-P) because they can be saved frequently to resolve ! the daily cycle while avoiding large files. ! DO ng=1,Ngrids LdefQCK(ng)=.TRUE. LwrtQCK(ng)=.TRUE. WRITE (QCK(ng)%name,10) TRIM(QCK(ng)%head), my_outer lstr=LEN_TRIM(QCK(ng)%name) QCK(ng)%base=QCK(ng)%name(1:lstr-3) END DO ! ! Define output 4D-Var NetCDF file (DAV struc) containing all ! processed data at observation locations. ! DO ng=1,Ngrids LdefMOD(ng)=.TRUE. CALL def_mod (ng) IF (FoundError(exit_flag, NoError, 330, MyFile)) RETURN END DO ! ! Stop profile clock ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_off (ng, iNLM, 86, 339, MyFile) END DO END DO ! 10 FORMAT (a,'_outer',i0,'.nc') ! RETURN END SUBROUTINE background_initialize ! SUBROUTINE background (my_outer, RunInterval) ! !======================================================================= ! ! ! This routine computes the backgound state trajectory, Xb_n-1(t), ! ! used to linearize the tangent linear and adjoint models in the ! ! inner loops. It interpolates the background at the observations ! ! locations, and computes the accept/reject quality control flag, ! ! ObsScale. ! ! ! ! On Input: ! ! ! ! my_outer Outer-loop counter (integer) ! ! RunInterval NLM kernel time stepping window (seconds) ! ! ! !======================================================================= ! ! Imported variable declarations ! integer, intent(in) :: my_outer ! real(dp), intent(in) :: RunInterval ! ! Local variable declarations. ! logical :: DoneStepping ! integer :: i, lstr, ng, tile integer :: Fcount, Tindex integer :: thread ! character (len=20) :: string character (len=*), parameter :: MyFile = & & "ROMS/Drivers/rbl4dvar.F"//", background" ! SourceFile=MyFile ! !----------------------------------------------------------------------- ! Set nonlinear model initial conditions. !----------------------------------------------------------------------- ! CALL background_initialize (my_outer) IF (FoundError(exit_flag, NoError, 401, MyFile)) RETURN ! !----------------------------------------------------------------------- ! Run nonlinear model and compute background state trajectory, ! Xb_n-1(t), and the background values at the observation points ! and times. It processes and writes the observations accept/reject ! flag (ObsScale) once to allow background quality control, if any. !----------------------------------------------------------------------- ! ! Start profile clock. ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_on (ng, iNLM, 86, 423, MyFile) END DO END DO ! ! Initialize various parameters and switches. ! MyRunInterval=RunInterval ! DO ng=1,Ngrids wrtMisfit(ng)=.TRUE. wrtObsScale(ng)=.TRUE. SporadicImpulse(ng)=.FALSE. FrequentImpulse(ng)=.FALSE. IF (Master) WRITE (stdout,20) 'NL', ng, my_outer, 0, & & ntstart(ng), ntend(ng) END DO ! ! Run nonlinear model kernel. ! CALL main3d (RunInterval) IF (FoundError(exit_flag, NoError, 479, MyFile)) RETURN ! !----------------------------------------------------------------------- ! If completed stepping, reset switches and write out cost function. !----------------------------------------------------------------------- ! DoneStepping=.TRUE. ! IF (DoneStepping) THEN DO ng=1,Ngrids LdefQCK(ng)=.FALSE. LwrtQCK(ng)=.FALSE. wrtNLmod(ng)=.FALSE. wrtObsScale(ng)=.FALSE. END DO ! ! Report data penalty function. ! DO ng=1,Ngrids IF (Master) THEN DO i=0,NobsVar(ng) IF (i.eq.0) THEN string='Total' ELSE string=ObsName(i) END IF IF (FOURDVAR(ng)%NLPenalty(i).ne.0.0_r8) THEN WRITE (stdout,30) my_outer, 0, 'NLM', & & FOURDVAR(ng)%NLPenalty(i), & & TRIM(string) END IF END DO END IF ! ! Write out initial data penalty function to NetCDF file. ! SourceFile=MyFile SELECT CASE (DAV(ng)%IOtype) CASE (io_nf90) CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, & & 'NL_iDataPenalty', & & FOURDVAR(ng)%NLPenalty(0:), & & (/1/), (/NobsVar(ng)+1/), & & ncid = DAV(ng)%ncid) END SELECT IF (FoundError(exit_flag, NoError, 554, MyFile)) RETURN ! ! Clean penalty array before next run of NL model. ! FOURDVAR(ng)%NLPenalty=0.0_r8 END DO END IF ! ! Stop profile clock ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_off (ng, iNLM, 86, 568, MyFile) END DO END DO ! 10 FORMAT (a,'_outer',i0,'.nc') 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', & & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, & ', TimeSteps: ',i0,' - ',i0,')',/) 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', & & 1p,e17.10,0p,t68,a) RETURN END SUBROUTINE background ! SUBROUTINE increment (my_outer, RunInterval) ! !======================================================================= ! ! ! This routine computes the 4D-Var data assimilation state increment, ! ! dXa, by iterating the inner loops and minimizing the cost function. ! ! ! ! On Input: ! ! ! ! my_outer Outer-loop counter (integer) ! ! RunInterval TLM/ADM kernels time stepping window (seconds) ! ! ! !======================================================================= ! USE comp_Jb0_mod, ONLY : aug_oper, comp_Jb0 USE ini_adjust_mod, ONLY : ini_adjust USE sum_grad_mod, ONLY : sum_grad USE sum_imp_mod, ONLY : sum_imp ! ! Imported variable declarations ! integer, intent(in) :: my_outer ! real(dp), intent(in) :: RunInterval ! ! Local variable declarations. ! logical :: Lcgini, Linner, Lposterior ! integer :: i, ifile, lstr, my_inner, ng, tile integer :: Fcount, InpRec integer :: ADrec, nLAST integer :: irec, jrec, jrec1, jrec2 integer :: LiNL integer :: thread ! integer, dimension(Ngrids) :: Nrec integer, dimension(Ngrids) :: nADrec ! character (len=8 ) :: driver = 'rbl4dvar' character (len=10) :: suffix character (len=*), parameter :: MyFile = & & "ROMS/Drivers/rbl4dvar.F"//", increment" ! SourceFile=MyFile ! !======================================================================= ! Compute 4D-Var increment. !======================================================================= ! ! Start profile clock. ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_on (ng, iTLM, 87, 653, MyFile) END DO END DO ! !----------------------------------------------------------------------- ! Set few variables needed in the "increment" phase. !----------------------------------------------------------------------- ! ! Set structure for the nonlinear forward trajectory to be processed ! by the tangent linear and adjoint models. Also, set switches to ! process the FWD structure in routine "check_multifile". Notice that ! it is possible to split solution into multiple NetCDF files to reduce ! their size. ! CALL edit_multifile ('HIS2FWD') IF (FoundError(exit_flag, NoError, 910, MyFile)) RETURN DO ng=1,Ngrids LreadFWD(ng)=.TRUE. END DO ! ! Set the BLK structure to contain the nonlinear model surface fluxes ! needed by the tangent linear and adjoint models. Also, set switches ! to process that structure in routine "check_multifile". Notice that ! it is possible to split the solution into multiple NetCDF files to ! reduce their size. ! ! The switch LreadFRC is deactivated because all the atmospheric ! forcing, including shortwave radiation, is read from the NLM ! surface fluxes or is assigned during ESM coupling. Such fluxes ! are available from the QCK structure. There is no need for reading ! and processing from the FRC structure input forcing-files. ! CALL edit_multifile ('QCK2BLK') IF (FoundError(exit_flag, NoError, 930, MyFile)) RETURN DO ng=1,Ngrids LreadBLK(ng)=.TRUE. LreadFRC(ng)=.FALSE. LreadQCK(ng)=.FALSE. END DO ! ! If weak constraint, set the number of TLF records used in the ! augmented solver. ! DO ng=1,Ngrids IF (nADJ(ng).lt.ntimes(ng)) THEN nADrec(ng)=2*(1+ntimes(ng)/nADJ(ng)) ELSE nADrec(ng)=0 END IF END DO ! !----------------------------------------------------------------------- ! On first pass (outer=1), create NetCDF files and initialize all the ! records in the ITL file to zero. !----------------------------------------------------------------------- ! CHECK_OUTER1 : IF (my_outer.eq.1) THEN ! ! If outer=1, define tangent linear initial conditions file. ! DO ng=1,Ngrids LdefITL(ng)=.TRUE. CALL tl_def_ini (ng) LdefITL(ng)=.FALSE. IF (FoundError(exit_flag, NoError, 965, MyFile)) RETURN END DO ! ! If outer=1, Initialize all records of the ITL file to zero. The TLM ! state variables are zero when outer=1. ! DO ng=1,Ngrids tdays(ng)=INItime(ng)*sec2day CALL tl_wrt_ini (ng, MyRank, Rec1, Rec1) IF (FoundError(exit_flag, NoError, 980, MyFile)) RETURN CALL tl_wrt_ini (ng, MyRank, Rec1, Rec2) IF (FoundError(exit_flag, NoError, 986, MyFile)) RETURN CALL tl_wrt_ini (ng, MyRank, Rec1, Rec3) IF (FoundError(exit_flag, NoError, 992, MyFile)) RETURN CALL tl_wrt_ini (ng, MyRank, Rec1, Rec4) IF (FoundError(exit_flag, NoError, 998, MyFile)) RETURN CALL tl_wrt_ini (ng, MyRank, Rec1, Rec5) IF (FoundError(exit_flag, NoError, 1004, MyFile)) RETURN ! IF (nADJ(ng).lt.ntimes(ng)) THEN nLAST=Rec5 DO irec=1,nADrec(ng) CALL tl_wrt_ini (ng, MyRank, Rec1, nLAST+irec) IF (FoundError(exit_flag, NoError, & & 1015, MyFile)) RETURN END DO END IF END DO ! ! If outer=1, define impulse forcing NetCDF file. ! DO ng=1,Ngrids LdefTLF(ng)=.TRUE. CALL def_impulse (ng) IF (FoundError(exit_flag, NoError, 1026, MyFile)) RETURN END DO END IF CHECK_OUTER1 ! !----------------------------------------------------------------------- ! If outer>1, integrate the tangent linear model to evaluate ! G[x(k)-x(k-1)], which is written to the DAV NetCDF in the ! variable "TLmodel_value". It is needed in routine "rpcg_lanczos". ! In the old unsplit algorithm, the TLM was run in the "analysis" ! phase. It is done here to allow lower resolution in the inner ! loops in the future. !----------------------------------------------------------------------- ! CHECK_OUTER2 : IF (my_outer.gt.1) THEN ! ! Clear tangent linear forcing arrays. This is very important since ! these arrays are non-zero and must be zero when running the tangent ! linear model. We also need to clear the nonlinear state arrays to ! make sure that the RHS terms rzeta, rubar, rvbar, ru, and rv are ! zero. Otherwise, those array will have the last computed values ! when running the nonlinear model if not processing the forward ! trajectory RHS terms. This needs to be done to get identical ! solutions with the split schemes. ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL initialize_ocean (ng, tile, iNLM) CALL initialize_forces (ng, tile, iTLM) END DO END DO ! ! Initialize tangent linear model from initial impulse which is now ! stored in file ITL(ng)%name. ! DO ng=1,Ngrids wrtNLmod(ng)=.FALSE. wrtTLmod(ng)=.TRUE. END DO ! ! If weak constraint, the impulses are time-interpolated at each ! time-steps. ! DO ng=1,Ngrids IF (FrcRec(ng).gt.3) THEN FrequentImpulse(ng)=.TRUE. END IF END DO ! ! Initialize tangent linear model from ITLname, record Rec3. ! The sum of the initial condition increments is in record 3. ! DO ng=1,Ngrids ITL(ng)%Rindex=Rec3 CALL tl_initial (ng) IF (FoundError(exit_flag, NoError, 1112, MyFile)) RETURN END DO ! ! Run tangent linear model. ! DO ng=1,Ngrids IF (Master) THEN WRITE (stdout,20) 'TL', ng, my_outer, 0, & & ntstart(ng), ntend(ng) END IF END DO ! CALL tl_main3d (RunInterval) IF (FoundError(exit_flag, NoError, 1129, MyFile)) RETURN ! DO ng=1,Ngrids wrtNLmod(ng)=.FALSE. wrtTLmod(ng)=.FALSE. END DO END IF CHECK_OUTER2 ! !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! 4D-Var inner loops. !::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::: ! ! Clear tangent linear forcing arrays before entering inner-loop. ! This is very important since these arrays are non-zero and must ! be zero when running the tangent linear model. ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL initialize_ocean (ng, tile, iNLM) CALL initialize_forces (ng, tile, iTLM) END DO END DO ! ! Inner loops iteration. ! INNER_LOOP : DO my_inner=0,Ninner inner=my_inner IF (inner.ne.Ninner) THEN Linner=.TRUE. ELSE Linner=.FALSE. END IF ! ! Retrieve TLmodVal and NLmodVal when inner=0 and outer>1 for use as ! Hbk and BCKmodVal, respectively. ! IF ((inner.eq.0).and.(my_outer.gt.1)) THEN DO ng=1,Ngrids SELECT CASE (DAV(ng)%IOtype) CASE (io_nf90) CALL netcdf_get_fvar (ng, iTLM, DAV(ng)%name, & & 'TLmodel_value', TLmodVal) IF (FoundError(exit_flag, NoError, 1196, & & MyFile)) RETURN CALL netcdf_get_fvar (ng, iTLM, DAV(ng)%name, & & 'NLmodel_value', NLmodVal, & & ncid = DAV(ng)%ncid, & & start = (/1,my_outer-1/), & & total = (/Ndatum(ng),1/)) IF (FoundError(exit_flag, NoError, 1204, & & MyFile)) RETURN END SELECT END DO END IF ! IF (inner.eq.0) THEN Lcgini=.TRUE. END IF DO ng=1,Ngrids CALL rpcg_lanczos (ng, iTLM, my_outer, inner, Ninner, Lcgini) END DO ! ! Start inner loop computations. ! INNER_COMPUTE : IF (Linner) THEN ! !----------------------------------------------------------------------- ! Integrate adjoint model forced with any vector PSI at the observation ! locations and generate adjoint trajectory, Lambda_n(t). !----------------------------------------------------------------------- ! ! Initialize the adjoint model from rest. ! DO ng=1,Ngrids CALL ad_initial (ng) IF (FoundError(exit_flag, NoError, 1267, MyFile)) RETURN wrtMisfit(ng)=.FALSE. END DO ! ! Set adjoint history NetCDF parameters. Define adjoint history ! file only once to avoid opening too many files. ! DO ng=1,Ngrids WRTforce(ng)=.TRUE. IF (Nrun.gt.1) LdefADJ(ng)=.FALSE. Fcount=ADM(ng)%load ADM(ng)%Nrec(Fcount)=0 ADM(ng)%Rindex=0 END DO ! ! Time-step adjoint model backwards forced with current PSI vector. ! DO ng=1,Ngrids IF (Master) THEN WRITE (stdout,20) 'AD', ng, my_outer, my_inner, & & ntstart(ng), ntend(ng) END IF END DO ! CALL ad_main3d (RunInterval) IF (FoundError(exit_flag, NoError, 1300, MyFile)) RETURN ! ! Activate "LwrtState2d" switch to write the correct ad_ubar, ad_vbar, ! ad_zeta here instead of ad_ubar_sol, ad_vbar_sol, and ad_zeta_sol in ! the calls to "ad_wrt_his" below (HGA). ! Here, ad_zeta(:,:,kout)=ad_zeta_sol. However, ad_ubar and ad_vbar ! not equal to ad_ubar_sol and ad_vbar_sol, respectively. It is ! irrelevant because ubar and vbar are not part of the state in ! 3D application. Notice that "LwrtState2d" will be turned off at ! the bottom of "error_covariance". ! DO ng=1,Ngrids LwrtState2d(ng)=.TRUE. END DO ! ! Write out last weak-constraint forcing (WRTforce is still .TRUE.) ! record into the adjoint history file. Note that the weak-constraint ! forcing is delayed by nADJ time-steps. ! DO ng=1,Ngrids CALL ad_wrt_his (ng, MyRank) IF (FoundError(exit_flag, NoError, 1327, MyFile)) RETURN END DO ! ! Write out adjoint initial condition record into the adjoint ! history file. ! DO ng=1,Ngrids WRTforce(ng)=.FALSE. CALL ad_wrt_his (ng, MyRank) IF (FoundError(exit_flag, NoError, 1341, MyFile)) RETURN END DO ! ! Convolve adjoint trajectory with error covariances. ! Lposterior=.FALSE. ! ! Set the flag that controls the augmentation of the model error ! forcing terms. This is ONLY done in the outer-loop so ! LaugWeak=.FALSE. here. ! LaugWeak=.FALSE. CALL error_covariance (iTLM, driver, my_outer, inner, & & Lbck, Lini, Lold, Lnew, & & Rec1, Rec2, Lposterior) IF (FoundError(exit_flag, NoError, 1362, MyFile)) RETURN ! ! Convert the current adjoint solution in ADM(ng)%name to impulse ! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file. ! To facilitate the forcing to the TLM and RPM, the forcing is ! processed and written in increasing time coordinates (recall that ! the adjoint solution in ADM(ng)%name is backwards in time). ! IF (Master) THEN WRITE (stdout,30) my_outer, inner END IF DO ng=1,Ngrids TLF(ng)%Rindex=0 CALL wrt_impulse (ng, MyRank, iADM, ADM(ng)%name) IF (FoundError(exit_flag, NoError, 1380, MyFile)) RETURN END DO ! !----------------------------------------------------------------------- ! Integrate tangent linear model forced by the convolved adjoint ! trajectory (impulse forcing) to compute R_n * PSI at observation ! points. !----------------------------------------------------------------------- ! ! Initialize tangent linear model from initial impulse which is now ! stored in file ITL(ng)%name. ! DO ng=1,Ngrids wrtNLmod(ng)=.FALSE. wrtTLmod(ng)=.TRUE. END DO ! ! If weak constraint, the impulses are time-interpolated at each ! time-steps. ! DO ng=1,Ngrids IF (FrcRec(ng).gt.3) THEN FrequentImpulse(ng)=.TRUE. END IF END DO ! ! Initialize tangent linear model from ITL(ng)%name, record Rec1. ! DO ng=1,Ngrids ITL(ng)%Rindex=Rec1 CALL tl_initial (ng) IF (FoundError(exit_flag, NoError, 1411, MyFile)) RETURN END DO ! ! Activate switch to write out initial misfit between model and ! observations. ! IF ((my_outer.eq.1).and.(inner.eq.1)) THEN DO ng=1,Ngrids wrtMisfit(ng)=.TRUE. END DO END IF ! ! Run tangent linear model forward and force with convolved adjoint ! trajectory impulses. Compute (H M B M' H')_n * PSI at observation ! points which are used in the conjugate gradient algorithm. ! DO ng=1,Ngrids IF (Master) THEN WRITE (stdout,20) 'TL', ng, my_outer, my_inner, & & ntstart(ng), ntend(ng) END IF END DO ! CALL tl_main3d (RunInterval) IF (FoundError(exit_flag, NoError, 1439, MyFile)) RETURN ! DO ng=1,Ngrids wrtNLmod(ng)=.FALSE. wrtTLmod(ng)=.FALSE. END DO ! Nrun=Nrun+1 Lcgini=.FALSE. END IF INNER_COMPUTE END DO INNER_LOOP ! !----------------------------------------------------------------------- ! Once the w_n, have been approximated with sufficient accuracy, ! compute estimates of Lambda_n and Xhat_n by carrying out one ! backward intergration of the adjoint model and one forward ! itegration of the nonlinear model. !----------------------------------------------------------------------- ! ! Initialize the adjoint model always from rest. ! DO ng=1,Ngrids CALL ad_initial (ng) IF (FoundError(exit_flag, NoError, 1505, MyFile)) RETURN END DO ! ! Set adjoint history NetCDF parameters. Define adjoint history ! file one to avoid opening to many files. ! DO ng=1,Ngrids WRTforce(ng)=.TRUE. IF (Nrun.gt.1) LdefADJ(ng)=.FALSE. Fcount=ADM(ng)%load ADM(ng)%Nrec(Fcount)=0 ADM(ng)%Rindex=0 END DO ! ! Time-step adjoint model backwards forced with estimated coefficients, ! Beta_n. ! DO ng=1,Ngrids IF (Master) THEN WRITE (stdout,20) 'AD', ng, my_outer, Ninner, & & ntstart(ng), ntend(ng) END IF END DO ! CALL ad_main3d (RunInterval) IF (FoundError(exit_flag, NoError, 1538, MyFile)) RETURN ! ! Write out last weak-constraint forcing (WRTforce is still .TRUE.) ! record into the adjoint history file. Note that the weak-constraint ! forcing is delayed by nADJ time-steps. ! DO ng=1,Ngrids CALL ad_wrt_his (ng, MyRank) IF (FoundError(exit_flag, NoError, 1552, MyFile)) RETURN END DO ! ! Write out adjoint initial condition record into the adjoint ! history file. ! DO ng=1,Ngrids WRTforce(ng)=.FALSE. CALL ad_wrt_his (ng, MyRank) IF (FoundError(exit_flag, NoError, 1566, MyFile)) RETURN END DO ! ! Get number of records in adjoint NetCDF. ! We need to do this here since ADM(ng)%Nrec is reset to zero in ! error_covariance. ! DO ng=1,Ngrids Fcount=ADM(ng)%load Nrec(ng)=ADM(ng)%Nrec(Fcount) END DO ! ! Convolve adjoint trajectory with error covariances. ! Lposterior=.FALSE. ! ! Set the flag that controls the augmentation of the model error ! forcing terms. This is ONLY done in the outer-loop so ! LaugWeak=.TRUE. here. ! LaugWeak=.TRUE. LiNL=my_outer+1 CALL error_covariance (iNLM, driver, my_outer, inner, & & LiNL, Lini, Lold, Lnew, & & Rec1, Rec2, Lposterior) IF (FoundError(exit_flag, NoError, 1596, MyFile)) RETURN ! ! Augmented solver: ! ! NOTES: The ITL file contains 5 records - ! Rec2 = the new TL initial condition ! Rec3 = the sum of the TL initial conditions ! Rec4 = B^-1(xb-xk)=sum_j^k-1 G_j lambda_j ! Rec5 = the augmented correction to Rec2. ! Rec5+1 to Rec5+nADrec/2 = sum of the TLF forcing increments ! Rec5+nADrec/2+1 to Rec5+nADrec = B^-1 of the sum of the TLF ! forcing increments. ! Reset the flag LaugWeak flag. ! LaugWeak=.FALSE. ! ! Compute the augmented correction to the adjoint propagator. ! We need to use sum (x(k)-x(k-1)) before it is updated. This is ! in record 3 of the ITL file. ! DO ng=1,Ngrids CALL get_state (ng, iTLM, 4, ITL(ng), Rec3, LTLM1) IF (FoundError(exit_flag, NoError, 1653, MyFile)) RETURN END DO DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL aug_oper (ng, tile, LTLM1, LTLM2) END DO END DO ! ! Save this augmented piece in record 5 of the ITL file. ! DO ng=1,Ngrids CALL tl_wrt_ini (ng, MyRank, LTLM2, Rec5) IF (FoundError(exit_flag, NoError, 1669, MyFile)) RETURN END DO ! ! Complete the computation of the TL initial condition by adding the ! contribution from the augmented adjoint propagator. ! DO ng=1,Ngrids CALL get_state (ng, iTLM, 4, ITL(ng), Rec2, LTLM1) IF (FoundError(exit_flag, NoError, 1677, MyFile)) RETURN CALL get_state (ng, iTLM, 4, ITL(ng), Rec5, LTLM2) IF (FoundError(exit_flag, NoError, 1679, MyFile)) RETURN END DO ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL sum_grad (ng, tile, LTLM1, LTLM2) END DO END DO ! ! Write the final TL increment to Rec2 of the ITL file. ! DO ng=1,Ngrids CALL tl_wrt_ini (ng, MyRank, LTLM2, Rec2) IF (FoundError(exit_flag, NoError, 1696, MyFile)) RETURN END DO ! ! Now update the non-linear model initial conditions. ! LiNL=my_outer+1 DO ng=1,Ngrids CALL get_state (ng, iNLM, 9, INI(ng), LiNL, Lnew(ng)) CALL get_state (ng, iADM, 4, ITL(ng), Rec2, LTLM1) END DO DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL ini_adjust (ng, tile, LTLM1, Lnew(ng)) END DO END DO ! ! Write out nonlinear model initial conditions into INIname, record ! INI(ng)%Rindex. ! DO ng=1,Ngrids CALL wrt_ini (ng, MyRank, Lnew(ng)) IF (FoundError(exit_flag, NoError, 1721, MyFile)) RETURN END DO ! ! Compute the new B^-1(x(k)-x(k-1)) term. ! Gather the final adjoint solutions increments, sum and ! save in record 4 of the ITL file. Use the tl arrays as temporary ! storage. ! ! First add the augmented piece which is computed from the previous ! sum. ! DO ng=1,Ngrids CALL get_state (ng, iTLM, 4, ITL(ng), Rec4, LTLM1) IF (FoundError(exit_flag, NoError, 1744, MyFile)) RETURN END DO ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL aug_oper (ng, tile, LTLM1, LTLM2) END DO END DO ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL sum_grad (ng, tile, LTLM1, LTLM2) END DO END DO ! DO ng=1,Ngrids ADrec=Nrec(ng) CALL get_state (ng, iTLM, 4, ADM(ng), ADrec, LTLM1) IF (FoundError(exit_flag, NoError, 1762, MyFile)) RETURN END DO ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL sum_grad (ng, tile, LTLM1, LTLM2) END DO END DO ! ! Write the current sum of adjoint solutions into record 4 of the ITL ! file. ! DO ng=1,Ngrids CALL tl_wrt_ini (ng, MyRank, LTLM2, Rec4) IF (FoundError(exit_flag, NoError, 1780, MyFile)) RETURN END DO ! ! Convert the current adjoint solution in ADM(ng)%name to impulse ! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file. ! To facilitate the forcing to the TLM and RPM, the forcing is ! processed and written in increasing time coordinates (recall that ! the adjoint solution in ADM(ng)%name is backwards in time). ! IF (Master) THEN WRITE (stdout,30) my_outer, inner END IF DO ng=1,Ngrids TLF(ng)%Rindex=0 CALL wrt_impulse (ng, MyRank, iADM, ADM(ng)%name) IF (FoundError(exit_flag, NoError, 1800, MyFile)) RETURN END DO ! ! Now Compute the augmented corrections to the weak constraint ! forcing terms. The sums so far are in records 6 to ! 6+nADrec/2. ! DO ng=1,Ngrids DO i=1,nADrec(ng)/2 irec=i jrec=Rec5+i CALL get_state (ng, iTLM, 4, ITL(ng), jrec, LTLM1) IF (FoundError(exit_flag, NoError, 1814, MyFile)) RETURN ! DO tile=first_tile(ng),last_tile(ng),+1 CALL aug_oper (ng, tile, LTLM1, LTLM2) END DO ! ! Complete the computation of the TLF forcing term by adding the ! contribution from the augmented adjoint propagator. Specify ! the value of 7 for the model variable since this the special ! case in get_state for reading the impulse forcing. ! CALL get_state (ng, 7, 4, TLF(ng), irec, LTLM1) IF (FoundError(exit_flag, NoError, 1827, MyFile)) RETURN ! DO tile=first_tile(ng),last_tile(ng),+1 CALL sum_imp (ng, tile, LTLM2) END DO ! ! Write the final forcing increment to he TLF file. ! Note the original TLF file is overwritten at this point. ! TLF(ng)%Rindex=0 CALL wrt_aug_imp (ng, MyRank, iTLM, LTLM2, i) END DO END DO ! ! Gather the increments from the final inner-loop and ! save in record 3 of the ITL file. The current increment ! is in record 2 and the sum so far is in record 3. ! DO ng=1,Ngrids CALL get_state (ng, iTLM, 4, ITL(ng), Rec2, LTLM1) IF (FoundError(exit_flag, NoError, 1853, MyFile)) RETURN CALL get_state (ng, iTLM, 4, ITL(ng), Rec3, LTLM2) IF (FoundError(exit_flag, NoError, 1855, MyFile)) RETURN END DO ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL sum_grad (ng, tile, LTLM1, LTLM2) END DO END DO ! ! Write the current sum into record 3 of the ITL file. ! DO ng=1,Ngrids CALL tl_wrt_ini (ng, MyRank, LTLM2, Rec3) IF (FoundError(exit_flag, NoError, 1872, MyFile)) RETURN END DO ! ! Gather the model error forcing increments and update the ! sum in records 6 to 6+nADrec/2. ! DO ng=1,Ngrids DO i=1,nADrec(ng)/2 irec=i jrec=Rec5+i CALL get_state (ng, iTLM, 4, ITL(ng), jrec, LTLM1) IF (FoundError(exit_flag, NoError, 1905, MyFile)) RETURN CALL get_state (ng, 7, 4, TLF(ng), irec, LTLM2) IF (FoundError(exit_flag, NoError, 1908, MyFile)) RETURN ! DO tile=first_tile(ng),last_tile(ng),+1 CALL sum_imp (ng, tile, LTLM1) END DO ! ! Write the current sum into record jrec of the ITL file. ! CALL tl_wrt_ini (ng, MyRank, LTLM1, jrec) IF (FoundError(exit_flag, NoError, 1924, MyFile)) RETURN END DO END DO ! ! Now compute the background cost function Jb0. ! First compute the contribution from the increments in the ! initial conditions, surface forcing, and boundary conditions. ! Jb0(my_outer)=0.0_r8 DO ng=1,Ngrids CALL get_state (ng, iADM, 8, ITL(ng), Rec4, LTLM1) IF (FoundError(exit_flag, NoError, 1935, MyFile)) RETURN CALL get_state (ng, iTLM, 8, ITL(ng), Rec3, LTLM1) IF (FoundError(exit_flag, NoError, 1937, MyFile)) RETURN END DO ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL comp_Jb0 (ng, tile, iTLM, my_outer, LTLM1, LTLM1) END DO END DO ! ! Now compute the contribution of model error terms to the ! background cost function Jb0. ! ! NOTE: I THINK WE NEED A BETTER WAY OF COMPUTING THE OBS ERROR ! CONTRIBUTION TO Jb. THE FOLLOWING CALCULATION IS DANGEROUS BECAUSE ! IT RELIES ON THE FACT THAT THE SURFACING FORCING AND OBC TL ARRAYS ! READ FROM jrec2 ARE ZERO. THIS MEANS THAT ONLY THE MODEL STATE ! CONTRIBUTES TO Jb FOR THE MODEL ERROR TERM, WHICH OF COURSE IS ! WHAT SHOULD BE THE CASE. THE SURFACE FORCING AND OBC CONTRIBUTIONS ! ARE COMPUTED IN THE PREVIOUS CALL TO COMP_JB0. ! DO ng=1,Ngrids DO irec=1,nADrec(ng)/2 jrec1=Rec5+irec jrec2=Rec5+nADrec(ng)/2+irec CALL get_state (ng, iADM, 8, ITL(ng), jrec1, LTLM1) IF (FoundError(exit_flag, NoError, 1962, MyFile)) RETURN CALL get_state (ng, iTLM, 8, ITL(ng), jrec2, LTLM1) IF (FoundError(exit_flag, NoError, 1964, MyFile)) RETURN ! DO tile=first_tile(ng),last_tile(ng),+1 CALL comp_Jb0 (ng, tile, iTLM, my_outer, LTLM1, LTLM1) END DO END DO END DO ! ! Overwrite the TFL netcdf file with the sum of the model error ! forcing increments - required for the next run of the NLM and ! TLM. ! DO ng=1,Ngrids DO irec=1,nADrec(ng)/2 jrec=Rec5+irec CALL get_state (ng, iTLM, 8, ITL(ng), jrec, LTLM1) ! TLF(ng)%Rindex=0 CALL wrt_aug_imp (ng, MyRank, iTLM, LTLM1, irec) END DO END DO ! ! Stop profile clock ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_off (ng, iTLM, 87, 1997, MyFile) END DO END DO ! 10 FORMAT (a,'_outer',i0,'.nc') 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', & & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, & ', TimeSteps: ',i0,' - ',i0,')',/) 30 FORMAT (/,' Converting Convolved Adjoint Trajectory to', & & ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/) ! RETURN END SUBROUTINE increment ! SUBROUTINE analysis_initialize (my_outer) ! !======================================================================= ! ! ! This routine initializes the nonlinear kernel with the 4D-Var new ! ! state estimate Xa = Xb + dXa. It is separated from the 'analysis' ! ! 4D-Var phase in ESM coupling applications that use generic methods ! ! for 'initialize', 'run', and 'finalize'. ! ! ! ! On Input: ! ! ! ! my_outer Outer-loop counter (integer) ! ! ! !======================================================================= ! ! Imported variable declarations ! integer, intent(in) :: my_outer ! ! Local variable declarations. ! integer :: ifile, lstr, ng, tile integer :: Fcount integer :: thread ! integer, dimension(Ngrids) :: indxSave ! character (len=10) :: suffix character (len=*), parameter :: MyFile = & & "ROMS/Drivers/rbl4dvar.F"//", analysis_initialize" ! SourceFile=MyFile ! !----------------------------------------------------------------------- ! Initialize nonlinear model kernel with new 4D-Var state estimate. !----------------------------------------------------------------------- ! ! Start profile clock. ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_on (ng, iNLM, 88, 2058, MyFile) END DO END DO ! ! Clear tangent arrays and the nonlinear model mixing arrays ! before running nonlinear model (important). ! DO ng=1,Ngrids DO tile=first_tile(ng),last_tile(ng),+1 CALL initialize_ocean (ng, tile, iTLM) CALL initialize_forces (ng, tile, iTLM) CALL initialize_mixing (ng, tile, iNLM) END DO END DO ! !----------------------------------------------------------------------- ! Initialize nonlinear model. !----------------------------------------------------------------------- ! ! Set new basic state trajectory for next outer loop. Notice that the ! LreadFWD switch is turned off to suppress processing of the structure ! when "check_multifile" during nonlinear model execution. ! DO ng=1,Ngrids idefHIS(ng)=-1 LdefHIS(ng)=.TRUE. LwrtHIS(ng)=.TRUE. wrtNLmod(ng)=.TRUE. wrtTLmod(ng)=.FALSE. LreadFWD(ng)=.FALSE. WRITE (HIS(ng)%name,10) TRIM(FWD(ng)%head), my_outer lstr=LEN_TRIM(HIS(ng)%name) HIS(ng)%base=HIS(ng)%name(1:lstr-3) IF (HIS(ng)%Nfiles.gt.1) THEN DO ifile=1,HIS(ng)%Nfiles WRITE (suffix,"('_',i4.4,'.nc')") ifile HIS(ng)%files(ifile)=TRIM(HIS(ng)%base)//TRIM(suffix) END DO HIS(ng)%name=TRIM(HIS(ng)%files(1)) ELSE HIS(ng)%files(1)=TRIM(HIS(ng)%name) END IF END DO ! ! Set nonlinear output quicksave history file for the next outer loop. ! Notice that the LreadBLK switch is turned off to suppress processing ! of the structure when "check_multifile" during nonlinear model ! execution. ! DO ng=1,Ngrids idefQCK(ng)=-1 LdefQCK(ng)=.TRUE. LwrtQCK(ng)=.TRUE. LreadBLK(ng)=.FALSE. WRITE (QCK(ng)%name,10) TRIM(QCK(ng)%head), my_outer lstr=LEN_TRIM(QCK(ng)%name) QCK(ng)%base=QCK(ng)%name(1:lstr-3) IF (QCK(ng)%Nfiles.gt.1) THEN DO ifile=1,QCK(ng)%Nfiles WRITE (suffix,"('_',i4.4,'.nc')") ifile QCK(ng)%files(ifile)=TRIM(QCK(ng)%base)//TRIM(suffix) END DO QCK(ng)%name=TRIM(QCK(ng)%files(1)) ELSE QCK(ng)%files(1)=TRIM(QCK(ng)%name) END IF END DO ! ! If weak constraint, the impulses are time-interpolated at each ! time-steps. ! DO ng=1,Ngrids IF (FrcRec(ng).gt.3) THEN FrequentImpulse(ng)=.TRUE. END IF END DO ! ! Initialize nonlinear model INI(ng)%name file, record outer+2. ! Notice that NetCDF record index counter is saved because this ! counter is used to write initial conditions. ! DO ng=1,Ngrids indxSave(ng)=INI(ng)%Rindex INI(ng)%Rindex=my_outer+2 RST(ng)%Rindex=0 Fcount=RST(ng)%load RST(ng)%Nrec(Fcount)=0 END DO ! CALL initial IF (FoundError(exit_flag, NoError, 2273, MyFile)) RETURN ! DO ng=1,Ngrids INI(ng)%Rindex=indxSave(ng) END DO ! ! Activate switch to write out final misfit between model and ! observations. ! IF (my_outer.eq.Nouter) THEN DO ng=1,Ngrids wrtMisfit(ng)=.TRUE. END DO END IF ! 10 FORMAT (a,'_outer',i0,'.nc') ! RETURN END SUBROUTINE analysis_initialize ! SUBROUTINE analysis (my_outer, RunInterval) ! !======================================================================= ! ! ! This routine computes 4D-Var data assimilation analysis, Xa. The ! ! nonlinear model initial conditions are computed by adding the ! ! 4D-Var increments to the current background: Xa = Xb + dXa. ! ! ! ! On Input: ! ! ! ! my_outer Outer-loop counter (integer) ! ! RunInterval NLM kernel time stepping window (seconds) ! ! ! !======================================================================= ! ! Imported variable declarations ! integer, intent(in) :: my_outer ! real(dp), intent(in) :: RunInterval ! ! Local variable declarations. ! logical :: DoneStepping ! integer :: i, lstr, ng integer :: Fcount integer :: thread ! character (len=20) :: string character (len=*), parameter :: MyFile = & & "ROMS/Drivers/rbl4dvar.F"//", analysis" ! SourceFile=MyFile ! ! !----------------------------------------------------------------------- ! Set nonlinear model initial conditions. !----------------------------------------------------------------------- ! CALL analysis_initialize (my_outer) IF (FoundError(exit_flag, NoError, 2343, MyFile)) RETURN ! !----------------------------------------------------------------------- ! Run nonlinear model and compute a "new estimate" of the state ! trajectory, Xb(t)|n. !----------------------------------------------------------------------- ! ! Start profile clock. ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_on (ng, iNLM, 88, 2363, MyFile) END DO END DO ! ! Initialize various parameters and switches. ! MyRunInterval=RunInterval ! DO ng=1,Ngrids IF (Master) WRITE (stdout,20) 'NL', ng, my_outer, Ninner, & & ntstart(ng), ntend(ng) END DO ! ! Run nonlinear model kernel. ! CALL main3d (RunInterval) IF (FoundError(exit_flag, NoError, 2416, MyFile)) RETURN ! !----------------------------------------------------------------------- ! If completed stepping, reset switches and write out cost function. !----------------------------------------------------------------------- ! DoneStepping=.TRUE. ! IF (DoneStepping) THEN DO ng=1,Ngrids wrtNLmod(ng)=.FALSE. wrtTLmod(ng)=.FALSE. END DO ! ! Set structure for the nonlinear forward trajectory to be processed ! by the tangent linear and adjoint models. Also, set switches to ! process the FWD structure in routine "check_multifile". Notice that ! it is possible to split solution into multiple NetCDF files to reduce ! their size. ! CALL edit_multifile ('HIS2FWD') IF (FoundError(exit_flag, NoError, 2458, MyFile)) RETURN DO ng=1,Ngrids LreadFWD(ng)=.TRUE. END DO ! ! Set the BLK structure to contain the nonlinear model surface fluxes ! needed by the tangent linear and adjoint models. Also, set switches ! to process that structure in routine "check_multifile". Notice that ! it is possible to split the solution into multiple NetCDF files to ! reduce their size. ! CALL edit_multifile ('QCK2BLK') IF (FoundError(exit_flag, NoError, 2472, MyFile)) RETURN DO ng=1,Ngrids LreadBLK(ng)=.TRUE. END DO ! ! Report data penalty function. ! DO ng=1,Ngrids IF (Master) THEN DO i=0,NobsVar(ng) IF (i.eq.0) THEN string='Total' ELSE string=ObsName(i) END IF IF (FOURDVAR(ng)%NLPenalty(i).ne.0.0_r8) THEN WRITE (stdout,30) my_outer, Ninner, 'NLM', & & FOURDVAR(ng)%NLPenalty(i), & & TRIM(string) END IF END DO END IF ! ! Write out final data penalty function to NetCDF file. ! SourceFile=MyFile SELECT CASE (DAV(ng)%IOtype) CASE (io_nf90) CALL netcdf_put_fvar (ng, iNLM, DAV(ng)%name, & & 'NL_fDataPenalty', & & FOURDVAR(ng)%NLPenalty(0:), & & (/1,my_outer/), & & (/NobsVar(ng)+1,1/), & & ncid = DAV(ng)%ncid) END SELECT IF (FoundError(exit_flag, NoError, 2518, MyFile)) RETURN ! ! Clean penalty array before next run of NL model. ! FOURDVAR(ng)%NLPenalty=0.0_r8 END DO END IF ! ! Stop profile clock ! DO ng=1,Ngrids DO thread=MyRank,MyRank CALL wclock_off (ng, iNLM, 88, 2532, MyFile) END DO END DO ! 10 FORMAT (a,'_outer',i0,'.nc') 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', & & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, & ', TimeSteps: ',i0,' - ',i0,')',/) 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', & & 1p,e17.10,0p,t68,a) ! RETURN END SUBROUTINE analysis ! SUBROUTINE prior_error (ng) ! !======================================================================= ! ! ! This routine processes background prior error covariance standard ! ! deviations and normalization coefficients. ! ! ! ! On Input: ! ! ! ! ng Nested grid number ! ! ! !======================================================================= ! USE def_norm_mod, ONLY : def_norm USE normalization_mod, ONLY : normalization USE strings_mod, ONLY : FoundError ! ! Imported variable declarations ! integer, intent(in) :: ng ! ! Local variable declarations. ! integer :: tile integer :: NRMrec, STDrec, Tindex ! character (len=*), parameter :: MyFile = & & "ROMS/Drivers/rbl4dvar.F"//", prior_error" ! SourceFile=MyFile ! !----------------------------------------------------------------------- ! Set application grid, metrics, and associated variables and ! parameters. !----------------------------------------------------------------------- ! ! The ROMS application grid configuration is done once. It is usually ! done in the "initial" kernel routine. However, since we are calling ! the "normalization" routine here, we need several grid variables and ! parameter. Also, if reading only water points, we need to know the ! land/sea mask arrays to unpack. ! IF (SetGridConfig(ng)) THEN CALL set_grid (ng, iNLM) END IF ! !----------------------------------------------------------------------- ! Read in standard deviation factors for error covariance. !----------------------------------------------------------------------- ! ! Initial conditions standard deviation. They are loaded in Tindex=1 ! of the e_var(...,Tindex) state variables. ! STDrec=1 Tindex=1 CALL get_state (ng, 10, 10, STD(1,ng), STDrec, Tindex) IF (FoundError(exit_flag, NoError, 2603, MyFile)) RETURN ! ! Model error standard deviation. They are loaded in Tindex=2 ! of the e_var(...,Tindex) state variables. ! STDrec=1 Tindex=2 IF (NSA.eq.2) THEN CALL get_state (ng, 11, 11, STD(2,ng), STDrec, Tindex) IF (FoundError(exit_flag, NoError, 2612, MyFile)) RETURN END IF ! !----------------------------------------------------------------------- ! Error covariance normalization coefficients. !----------------------------------------------------------------------- ! ! Compute or read in the error covariance normalization factors. ! If computing, write out factors to NetCDF. This is an expensive ! computation that needs to be computed only once for a particular ! application grid and decorrelation scales. ! IF (ANY(LwrtNRM(:,ng))) THEN CALL def_norm (ng, iNLM, 1) IF (FoundError(exit_flag, NoError, 2646, MyFile)) RETURN IF (NSA.eq.2) THEN CALL def_norm (ng, iNLM, 2) IF (FoundError(exit_flag, NoError, 2650, MyFile)) RETURN END IF ! DO tile=first_tile(ng),last_tile(ng),+1 CALL normalization (ng, tile, 2) END DO IF (FoundError(exit_flag, NoError, 2666, MyFile)) RETURN LdefNRM(1:4,ng)=.FALSE. LwrtNRM(1:4,ng)=.FALSE. ELSE NRMrec=1 CALL get_state (ng, 14, 14, NRM(1,ng), NRMrec, 1) IF (FoundError(exit_flag, NoError, 2674, MyFile)) RETURN IF (NSA.eq.2) THEN CALL get_state (ng, 15, 15, NRM(2,ng), NRMrec, 2) IF (FoundError(exit_flag, NoError, 2678, MyFile)) RETURN END IF END IF ! RETURN END SUBROUTINE prior_error END MODULE rbl4dvar_mod