ROMS
Loading...
Searching...
No Matches
rbl4dvar_mod Module Reference

Functions/Subroutines

subroutine, public background_initialize (my_outer)
 
subroutine, public background (my_outer, runinterval)
 
subroutine, public increment (my_outer, runinterval)
 
subroutine, public analysis_initialize (my_outer)
 
subroutine, public analysis (my_outer, runinterval)
 
subroutine, public prior_error (ng)
 
subroutine, public posterior_error (runinterval)
 

Variables

logical ldone
 
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
 

Function/Subroutine Documentation

◆ analysis()

subroutine, public rbl4dvar_mod::analysis ( integer, intent(in) my_outer,
real(dp), intent(in) runinterval )

Definition at line 2372 of file rbl4dvar.F.

2373!
2374!=======================================================================
2375! !
2376! This routine computes 4D-Var data assimilation analysis, Xa. The !
2377! nonlinear model initial conditions are computed by adding the !
2378! 4D-Var increments to the current background: Xa = Xb + dXa. !
2379! !
2380! On Input: !
2381! !
2382! my_outer Outer-loop counter (integer) !
2383! RunInterval NLM kernel time stepping window (seconds) !
2384! !
2385!=======================================================================
2386!
2387! Imported variable declarations
2388!
2389 integer, intent(in) :: my_outer
2390!
2391 real(dp), intent(in) :: RunInterval
2392!
2393! Local variable declarations.
2394!
2395 logical :: DoneStepping
2396!
2397 integer :: i, lstr, ng
2398 integer :: Fcount
2399# ifdef PROFILE
2400 integer :: thread
2401# endif
2402# if defined MODEL_COUPLING && !defined MCT_LIB
2403 integer :: NstrStep, NendStep, extra
2404!
2405 real(dp) :: ENDtime, NEXTtime
2406# endif
2407!
2408 character (len=20) :: string
2409
2410 character (len=*), parameter :: MyFile = &
2411 & __FILE__//", analysis"
2412!
2413 sourcefile=myfile
2414!
2415# if !(defined MODEL_COUPLING && defined ESMF_LIB)
2416!
2417!-----------------------------------------------------------------------
2418! Set nonlinear model initial conditions.
2419!-----------------------------------------------------------------------
2420!
2421 CALL analysis_initialize (my_outer)
2422 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2423# endif
2424!
2425!-----------------------------------------------------------------------
2426! Run nonlinear model and compute a "new estimate" of the state
2427! trajectory, Xb(t)|n.
2428# if defined MODEL_COUPLING && !defined MCT_LIB
2429! Since the ROMS kernel has a delayed output and line diagnostics by
2430! one timestep, subtact an extra value to the report of starting and
2431! ending timestep for clarity. Usually, the model coupling interval
2432! is of the same size as ROMS timestep.
2433# endif
2434!-----------------------------------------------------------------------
2435
2436# ifdef PROFILE
2437!
2438! Start profile clock.
2439!
2440 DO ng=1,ngrids
2441 DO thread=thread_range
2442 CALL wclock_on (ng, inlm, 88, __line__, myfile)
2443 END DO
2444 END DO
2445# endif
2446!
2447! Initialize various parameters and switches.
2448!
2449 myruninterval=runinterval
2450!
2451 DO ng=1,ngrids
2452# ifdef AVERAGES
2453 idefavg(ng)=-1
2454 ldefavg(ng)=.true.
2455 lwrtavg(ng)=.true.
2456 WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
2457 lstr=len_trim(avg(ng)%name)
2458 avg(ng)%base=avg(ng)%name(1:lstr-3)
2459# endif
2460# ifdef DIAGNOSTICS
2461 idefdia(ng)=-1
2462 ldefdia(ng)=.true.
2463 lwrtdia(ng)=.true.
2464 WRITE (dia(ng)%name,10) trim(dia(ng)%head), my_outer
2465 lstr=len_trim(dia(ng)%name)
2466 dia(ng)%base=dia(ng)%name(1:lstr-3)
2467# endif
2468# if defined MODEL_COUPLING && !defined MCT_LIB
2469!
2470 nexttime=time(ng)+runinterval
2471 endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
2472 IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
2473 extra=0 ! last time interval
2474 ELSE
2475 extra=1
2476 END IF
2477 step_counter(ng)=0
2478 nstrstep=iic(ng)
2479 nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
2480 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, ninner, &
2481 & nstrstep, nendstep
2482# else
2483 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, ninner, &
2484 & ntstart(ng), ntend(ng)
2485# endif
2486 END DO
2487!
2488! Run nonlinear model kernel.
2489!
2490# ifdef SOLVE3D
2491 CALL main3d (runinterval)
2492# else
2493 CALL main2d (runinterval)
2494# endif
2495 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2496!
2497!-----------------------------------------------------------------------
2498! If completed stepping, reset switches and write out cost function.
2499# if defined MODEL_COUPLING && !defined MCT_LIB
2500! In coupled applications, RunInterval is much less than ntimes*dt,
2501! so we need to wait until the last coupling interval is finished.
2502! Otherwise, the control switches will be turned off prematurely.
2503# endif
2504!-----------------------------------------------------------------------
2505!
2506# if defined MODEL_COUPLING && !defined MCT_LIB
2507 IF (nendstep.eq.ntend(1)) THEN
2508 donestepping=.true.
2509 ELSE
2510 donestepping=.false.
2511 END IF
2512# else
2513 donestepping=.true.
2514# endif
2515!
2516 IF (donestepping) THEN
2517 DO ng=1,ngrids
2518# ifdef AVERAGES
2519 ldefavg(ng)=.false.
2520 lwrtavg(ng)=.false.
2521# endif
2522# ifdef DIAGNOSTICS
2523 ldefdia(ng)=.false.
2524 lwrtdia(ng)=.false.
2525# endif
2526 wrtnlmod(ng)=.false.
2527 wrttlmod(ng)=.false.
2528 END DO
2529
2530# ifdef RPCG
2531!
2532! Deactivate the NL forcing term for iau.
2533!
2534 DO ng=1,ngrids
2535 iauswitch(ng)=.false.
2536 END DO
2537# endif
2538!
2539! Set structure for the nonlinear forward trajectory to be processed
2540! by the tangent linear and adjoint models. Also, set switches to
2541! process the FWD structure in routine "check_multifile". Notice that
2542! it is possible to split solution into multiple NetCDF files to reduce
2543! their size.
2544!
2545 CALL edit_multifile ('HIS2FWD')
2546 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2547 DO ng=1,ngrids
2548 lreadfwd(ng)=.true.
2549 END DO
2550
2551# ifdef FORWARD_FLUXES
2552!
2553! Set the BLK structure to contain the nonlinear model surface fluxes
2554! needed by the tangent linear and adjoint models. Also, set switches
2555! to process that structure in routine "check_multifile". Notice that
2556! it is possible to split the solution into multiple NetCDF files to
2557! reduce their size.
2558!
2559 CALL edit_multifile ('QCK2BLK')
2560 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2561 DO ng=1,ngrids
2562 lreadblk(ng)=.true.
2563 END DO
2564# endif
2565!
2566! Report data penalty function.
2567!
2568 DO ng=1,ngrids
2569 IF (master) THEN
2570 DO i=0,nobsvar(ng)
2571 IF (i.eq.0) THEN
2572 string='Total'
2573 ELSE
2574 string=obsname(i)
2575 END IF
2576 IF (fourdvar(ng)%NLPenalty(i).ne.0.0_r8) THEN
2577 WRITE (stdout,30) my_outer, ninner, 'NLM', &
2578 & fourdvar(ng)%NLPenalty(i), &
2579 & trim(string)
2580 END IF
2581 END DO
2582 END IF
2583!
2584! Write out final data penalty function to NetCDF file.
2585!
2586 sourcefile=myfile
2587 SELECT CASE (dav(ng)%IOtype)
2588 CASE (io_nf90)
2589 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
2590 & 'NL_fDataPenalty', &
2591 & fourdvar(ng)%NLPenalty(0:), &
2592 & (/1,my_outer/), &
2593 & (/nobsvar(ng)+1,1/), &
2594 & ncid = dav(ng)%ncid)
2595
2596# if defined PIO_LIB && defined DISTRIBUTE
2597 CASE (io_pio)
2598 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
2599 & 'NL_fDataPenalty', &
2600 & fourdvar(ng)%NLPenalty(0:), &
2601 & (/1,my_outer/), &
2602 & (/nobsvar(ng)+1,1/), &
2603 & piofile = dav(ng)%pioFile)
2604# endif
2605 END SELECT
2606 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2607!
2608! Clean penalty array before next run of NL model.
2609!
2610 fourdvar(ng)%NLPenalty=0.0_r8
2611 END DO
2612 END IF
2613
2614# ifdef PROFILE
2615!
2616! Stop profile clock
2617!
2618 DO ng=1,ngrids
2619 DO thread=thread_range
2620 CALL wclock_off (ng, inlm, 88, __line__, myfile)
2621 END DO
2622 END DO
2623# endif
2624!
2625 10 FORMAT (a,'_outer',i0,'.nc')
2626 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
2627 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
2628 ', TimeSteps: ',i0,' - ',i0,')',/)
2629 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', &
2630 & 1p,e17.10,0p,t68,a)
2631!
2632 RETURN
subroutine edit_multifile(task)
subroutine main2d
Definition main2d.F:746
subroutine main3d(runinterval)
Definition main3d.F:4
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3

References analysis_initialize(), mod_iounits::avg, mod_iounits::dav, mod_iounits::dia, mod_scalars::dt, edit_multifile(), mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_scalars::iauswitch, mod_ncparam::idefavg, mod_ncparam::idefdia, mod_scalars::iic, mod_scalars::initime, mod_param::inlm, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_scalars::ldefavg, mod_scalars::ldefdia, mod_scalars::lreadblk, mod_scalars::lreadfwd, mod_scalars::lwrtavg, mod_scalars::lwrtdia, main2d(), main3d(), mod_parallel::master, mod_scalars::myruninterval, mod_param::ngrids, mod_scalars::ninner, mod_fourdvar::nobsvar, mod_scalars::noerror, mod_scalars::ntend, mod_scalars::ntimes, mod_scalars::ntstart, mod_fourdvar::obsname, mod_iounits::sourcefile, mod_iounits::stdout, mod_scalars::step_counter, mod_scalars::time, wclock_off(), wclock_on(), mod_fourdvar::wrtnlmod, and mod_fourdvar::wrttlmod.

Here is the call graph for this function:

◆ analysis_initialize()

subroutine, public rbl4dvar_mod::analysis_initialize ( integer, intent(in) my_outer)

Definition at line 2066 of file rbl4dvar.F.

2067!
2068!=======================================================================
2069! !
2070! This routine initializes the nonlinear kernel with the 4D-Var new !
2071! state estimate Xa = Xb + dXa. It is separated from the 'analysis' !
2072! 4D-Var phase in ESM coupling applications that use generic methods !
2073! for 'initialize', 'run', and 'finalize'. !
2074! !
2075! On Input: !
2076! !
2077! my_outer Outer-loop counter (integer) !
2078! !
2079!=======================================================================
2080!
2081! Imported variable declarations
2082!
2083 integer, intent(in) :: my_outer
2084!
2085! Local variable declarations.
2086!
2087 integer :: ifile, lstr, ng, tile
2088 integer :: Fcount
2089# ifdef PROFILE
2090 integer :: thread
2091# endif
2092!
2093 integer, dimension(Ngrids) :: indxSave
2094!
2095 character (len=10) :: suffix
2096
2097 character (len=*), parameter :: MyFile = &
2098 & __FILE__//", analysis_initialize"
2099!
2100 sourcefile=myfile
2101!
2102!-----------------------------------------------------------------------
2103! Initialize nonlinear model kernel with new 4D-Var state estimate.
2104!-----------------------------------------------------------------------
2105
2106# ifdef PROFILE
2107!
2108! Start profile clock.
2109!
2110 DO ng=1,ngrids
2111 DO thread=thread_range
2112 CALL wclock_on (ng, inlm, 88, __line__, myfile)
2113 END DO
2114 END DO
2115# endif
2116!
2117! Clear tangent arrays and the nonlinear model mixing arrays
2118! before running nonlinear model (important).
2119!
2120 DO ng=1,ngrids
2121 DO tile=first_tile(ng),last_tile(ng),+1
2122 CALL initialize_ocean (ng, tile, itlm)
2123 CALL initialize_forces (ng, tile, itlm)
2124# ifdef ADJUST_BOUNDARY
2125 CALL initialize_boundary (ng, tile, itlm)
2126# endif
2127 CALL initialize_mixing (ng, tile, inlm)
2128 END DO
2129 END DO
2130
2131# ifdef SPLIT_4DVAR
2132!
2133!-----------------------------------------------------------------------
2134! If split 4D-Var algorithm, set several variables that computed or
2135! assigned in other 4D-Var phase executable.
2136!-----------------------------------------------------------------------
2137!
2138! Set Nrun>1, to read in surface forcing and open boundary conditions
2139! increments in "initial", if appropriate.
2140!
2141 nrun=ninner+1
2142!
2143! Set ERstr=Nrun, to set the open boundary condition (OBC_time) and
2144! surface forcing (SF_time) adjustment times, if needed.
2145!
2146 erstr=nrun
2147# if defined MODEL_COUPLING && defined ESMF_LIB
2148 erend=nouter
2149# endif
2150!
2151! Open Nonlinear model initial conditions NetCDF file (INI struc) and
2152! inquire about its variable IDs. In particular, set the current value
2153! of INI(ng)%Rindex, which is needed in "initial".
2154!
2155 DO ng=1,ngrids
2156 ldefini(ng)=.false.
2157 CALL def_ini (ng)
2158 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2159 END DO
2160!
2161! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
2162!
2163 DO ng=1,ngrids
2164 ldefmod(ng)=.false.
2165 CALL def_mod (ng)
2166 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2167 END DO
2168!
2169! In the split 4D-Var algorithm, we need to read the values of ObsScale
2170! computed in the background phase from the DAV NetCDF file. Such data
2171! is in memory in the unsplit algorithm.
2172!
2173! HGA: What to do in 4D-Var with nested grids?
2174!
2175 ng=1
2176 SELECT CASE (dav(ng)%IOtype)
2177 CASE (io_nf90)
2178 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
2179 & vname(1,idobss), obsscale, &
2180 & ncid = dav(ng)%ncid)
2181
2182# if defined PIO_LIB && defined DISTRIBUTE
2183 CASE (io_pio)
2184 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
2185 & vname(1,idobss), obsscale, &
2186 & piofile = dav(ng)%pioFile)
2187# endif
2188 END SELECT
2189 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2190!
2191! If outer=Nouter, read initial NLmodVal ("NLmodel_initial") and load
2192! it in the scratch "NLincrement" array.
2193!
2194 IF (my_outer.eq.nouter) THEN
2195 ng=1
2196 SELECT CASE (dav(ng)%IOtype)
2197 CASE (io_nf90)
2198 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
2199 & vname(1,idnlmi), nlincrement, &
2200 & ncid = dav(ng)%ncid)
2201
2202# if defined PIO_LIB && defined DISTRIBUTE
2203 CASE (io_pio)
2204 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
2205 & vname(1,idnlmi), nlincrement, &
2206 & piofile = dav(ng)%pioFile)
2207# endif
2208 END SELECT
2209 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2210 END IF
2211
2212# ifdef FORWARD_FLUXES
2213!
2214! If not first NLM run, set BLK structure to the previous outer loop
2215! quicksave trajectory. There is a logic in "get_data" that reads
2216! from BLK.
2217! (HGA: This probably legacy code that it is no longer needed)
2218!
2219 DO ng=1,ngrids
2220 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
2221 lstr=len_trim(qck(ng)%name)
2222 qck(ng)%base=qck(ng)%name(1:lstr-3)
2223 END DO
2224!
2225 CALL edit_multifile ('QCK2BLK')
2226 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2227# endif
2228# endif
2229!
2230!-----------------------------------------------------------------------
2231! Initialize nonlinear model.
2232!-----------------------------------------------------------------------
2233!
2234! Set new basic state trajectory for next outer loop. Notice that the
2235! LreadFWD switch is turned off to suppress processing of the structure
2236! when "check_multifile" during nonlinear model execution.
2237!
2238 DO ng=1,ngrids
2239 idefhis(ng)=-1
2240 ldefhis(ng)=.true.
2241 lwrthis(ng)=.true.
2242 wrtnlmod(ng)=.true.
2243 wrttlmod(ng)=.false.
2244# ifdef FRC_FILES
2245 lreadfrc(ng)=.true.
2246# endif
2247 lreadfwd(ng)=.false.
2248 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer
2249 lstr=len_trim(his(ng)%name)
2250 his(ng)%base=his(ng)%name(1:lstr-3)
2251 IF (his(ng)%Nfiles.gt.1) THEN
2252 DO ifile=1,his(ng)%Nfiles
2253 WRITE (suffix,"('_',i4.4,'.nc')") ifile
2254 his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
2255 END DO
2256 his(ng)%name=trim(his(ng)%files(1))
2257 ELSE
2258 his(ng)%files(1)=trim(his(ng)%name)
2259 END IF
2260 END DO
2261!
2262! Set nonlinear output quicksave history file for the next outer loop.
2263! Notice that the LreadBLK switch is turned off to suppress processing
2264! of the structure when "check_multifile" during nonlinear model
2265! execution.
2266!
2267 DO ng=1,ngrids
2268 idefqck(ng)=-1
2269 ldefqck(ng)=.true.
2270 lwrtqck(ng)=.true.
2271# ifdef FORWARD_FLUXES
2272 lreadblk(ng)=.false.
2273# endif
2274 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer
2275 lstr=len_trim(qck(ng)%name)
2276 qck(ng)%base=qck(ng)%name(1:lstr-3)
2277 IF (qck(ng)%Nfiles.gt.1) THEN
2278 DO ifile=1,qck(ng)%Nfiles
2279 WRITE (suffix,"('_',i4.4,'.nc')") ifile
2280 qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
2281 END DO
2282 qck(ng)%name=trim(qck(ng)%files(1))
2283 ELSE
2284 qck(ng)%files(1)=trim(qck(ng)%name)
2285 END IF
2286 END DO
2287!
2288! If weak constraint, the impulses are time-interpolated at each
2289! time-steps.
2290!
2291 DO ng=1,ngrids
2292# ifdef SPLIT_4DVAR
2293 IF (nadj(ng).lt.ntimes(ng)) THEN
2294 SELECT CASE (adm(ng)%IOtype)
2295 CASE (io_nf90)
2296 CALL netcdf_get_dim (ng, inlm, adm(ng)%name, &
2297 & dimname = 'ocean_time', &
2298 & dimsize = frcrec(ng))
2299# if defined PIO_LIB && defined DISTRIBUTE
2300 CASE (io_pio)
2301 CALL pio_netcdf_get_dim (ng, inlm, adm(ng)%name, &
2302 & dimname = 'ocean_time', &
2303 & dimsize = frcrec(ng))
2304# endif
2305 END SELECT
2306 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2307 END IF
2308# endif
2309 IF (frcrec(ng).gt.3) THEN
2310 frequentimpulse(ng)=.true.
2311 END IF
2312 END DO
2313!
2314! Initialize nonlinear model INI(ng)%name file, record outer+2.
2315! Notice that NetCDF record index counter is saved because this
2316! counter is used to write initial conditions.
2317!
2318 DO ng=1,ngrids
2319 indxsave(ng)=ini(ng)%Rindex
2320 ini(ng)%Rindex=my_outer+2
2321 rst(ng)%Rindex=0
2322 fcount=rst(ng)%load
2323 rst(ng)%Nrec(fcount)=0
2324 END DO
2325
2326# ifdef RPCG
2327!
2328! Compute the incremental analysis update forcing for the NLM.
2329!
2330! Read the accumulated increment from Rec3 of the ITL file.
2331!
2332 DO ng=1,ngrids
2333 IF (timeiau(ng).gt.0.0_dp) THEN
2334 iauswitch(ng)=.true.
2335 CALL get_state (ng, iadm, 4, itl(ng), rec3, ltlm1)
2336 END IF
2337 END DO
2338!
2339! Set up the NL model forcing contribution for the incremental
2340! analysis update.
2341!
2342 DO ng=1,ngrids
2343 IF (timeiau(ng).gt.0.0_dp) THEN
2344 DO tile=first_tile(ng),last_tile(ng),+1
2345 CALL frc_iau_ini (ng, tile, ltlm1)
2346 END DO
2347 END IF
2348 END DO
2349# endif
2350!
2351 CALL initial
2352 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2353!
2354 DO ng=1,ngrids
2355 ini(ng)%Rindex=indxsave(ng)
2356 END DO
2357!
2358! Activate switch to write out final misfit between model and
2359! observations.
2360!
2361 IF (my_outer.eq.nouter) THEN
2362 DO ng=1,ngrids
2363 wrtmisfit(ng)=.true.
2364 END DO
2365 END IF
2366!
2367 10 FORMAT (a,'_outer',i0,'.nc')
2368!
2369 RETURN
subroutine initial
Definition initial.F:3

References mod_iounits::adm, mod_iounits::dav, def_ini_mod::def_ini(), def_mod_mod::def_mod(), edit_multifile(), mod_scalars::erend, mod_scalars::erstr, mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), frc_iau_mod::frc_iau_ini(), mod_scalars::frcrec, mod_scalars::frequentimpulse, mod_iounits::fwd, get_state_mod::get_state(), mod_iounits::his, mod_param::iadm, mod_scalars::iauswitch, mod_ncparam::idefhis, mod_ncparam::idefqck, mod_ncparam::idnlmi, mod_ncparam::idobss, mod_iounits::ini, initial(), mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_mixing::initialize_mixing(), mod_ocean::initialize_ocean(), mod_param::inlm, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_iounits::itl, mod_param::itlm, mod_parallel::last_tile, mod_scalars::ldefhis, mod_scalars::ldefini, mod_scalars::ldefmod, mod_scalars::ldefqck, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, ltlm1, mod_scalars::lwrthis, mod_scalars::lwrtqck, mod_scalars::nadj, mod_netcdf::netcdf_get_dim(), mod_param::ngrids, mod_scalars::ninner, mod_fourdvar::nlincrement, mod_scalars::noerror, mod_scalars::nouter, mod_scalars::nrun, mod_scalars::ntimes, mod_fourdvar::obsscale, mod_pio_netcdf::pio_netcdf_get_dim(), mod_iounits::qck, rec3, mod_iounits::rst, mod_iounits::sourcefile, mod_scalars::timeiau, mod_ncparam::vname, wclock_on(), mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, and mod_fourdvar::wrttlmod.

Referenced by analysis().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ background()

subroutine, public rbl4dvar_mod::background ( integer, intent(in) my_outer,
real(dp), intent(in) runinterval )

Definition at line 399 of file rbl4dvar.F.

400!
401!=======================================================================
402! !
403! This routine computes the backgound state trajectory, Xb_n-1(t), !
404! used to linearize the tangent linear and adjoint models in the !
405! inner loops. It interpolates the background at the observations !
406! locations, and computes the accept/reject quality control flag, !
407! ObsScale. !
408! !
409! On Input: !
410! !
411! my_outer Outer-loop counter (integer) !
412! RunInterval NLM kernel time stepping window (seconds) !
413! !
414!=======================================================================
415!
416! Imported variable declarations
417!
418 integer, intent(in) :: my_outer
419!
420 real(dp), intent(in) :: RunInterval
421!
422! Local variable declarations.
423!
424 logical :: DoneStepping
425!
426 integer :: i, lstr, ng, tile
427 integer :: Fcount, Tindex
428# ifdef PROFILE
429 integer :: thread
430# endif
431# if defined MODEL_COUPLING && !defined MCT_LIB
432 integer :: NstrStep, NendStep, extra
433!
434 real(dp) :: ENDtime, NEXTtime
435# endif
436!
437 character (len=20) :: string
438
439 character (len=*), parameter :: MyFile = &
440 & __FILE__//", background"
441!
442 sourcefile=myfile
443
444# if !(defined MODEL_COUPLING && defined ESMF_LIB)
445!
446!-----------------------------------------------------------------------
447! Set nonlinear model initial conditions.
448!-----------------------------------------------------------------------
449!
450 CALL background_initialize (my_outer)
451 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
452# endif
453!
454!-----------------------------------------------------------------------
455! Run nonlinear model and compute background state trajectory,
456! Xb_n-1(t), and the background values at the observation points
457! and times. It processes and writes the observations accept/reject
458! flag (ObsScale) once to allow background quality control, if any.
459# if defined MODEL_COUPLING && !defined MCT_LIB
460! Since the ROMS kernel has a delayed output and line diagnostics by
461! one timestep, subtact an extra value to the report of starting and
462! ending timestep for clarity. Usually, the model coupling interval
463! is of the same size as ROMS timestep.
464# endif
465!-----------------------------------------------------------------------
466
467# ifdef PROFILE
468!
469! Start profile clock.
470!
471 DO ng=1,ngrids
472 DO thread=thread_range
473 CALL wclock_on (ng, inlm, 86, __line__, myfile)
474 END DO
475 END DO
476# endif
477!
478! Initialize various parameters and switches.
479!
480 myruninterval=runinterval
481!
482 DO ng=1,ngrids
483# ifdef AVERAGES
484 ldefavg(ng)=.true.
485 lwrtavg(ng)=.true.
486 WRITE (avg(ng)%name,10) trim(avg(ng)%head), my_outer
487 lstr=len_trim(avg(ng)%name)
488 avg(ng)%base=avg(ng)%name(1:lstr-3)
489# endif
490# ifdef DIAGNOSTICS
491 ldefdia(ng)=.true.
492 lwrtdia(ng)=.true.
493 WRITE (dia(ng)%name,10) trim(dia(ng)%head), my_outer
494 lstr=len_trim(dia(ng)%name)
495 dia(ng)%base=dia(ng)%name(1:lstr-3)
496# endif
497 wrtmisfit(ng)=.true.
498 wrtobsscale(ng)=.true.
499 sporadicimpulse(ng)=.false.
500 frequentimpulse(ng)=.false.
501# if defined MODEL_COUPLING && !defined MCT_LIB
502!
503 nexttime=time(ng)+runinterval
504 endtime=initime(ng)+(ntimes(ng)-1)*dt(ng)
505 IF ((nexttime.eq.endtime).and.(ng.eq.1)) THEN
506 extra=0 ! last time interval
507 ELSE
508 extra=1
509 END IF
510 step_counter(ng)=0
511 nstrstep=iic(ng)
512 nendstep=nstrstep+int((myruninterval)/dt(ng))-extra
513 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0, &
514 & nstrstep, nendstep
515# else
516
517 IF (master) WRITE (stdout,20) 'NL', ng, my_outer, 0, &
518 & ntstart(ng), ntend(ng)
519# endif
520 END DO
521!
522! Run nonlinear model kernel.
523!
524# ifdef SOLVE3D
525 CALL main3d (runinterval)
526# else
527 CALL main2d (runinterval)
528# endif
529 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
530!
531!-----------------------------------------------------------------------
532! If completed stepping, reset switches and write out cost function.
533# if defined MODEL_COUPLING && !defined MCT_LIB
534! In coupled applications, RunInterval is much less than ntimes*dt,
535! so we need to wait until the last coupling interval is finished.
536! Otherwise, the control switches will be turned off prematurely.
537# endif
538!-----------------------------------------------------------------------
539!
540# if defined MODEL_COUPLING && !defined MCT_LIB
541 IF (nendstep.eq.ntend(1)) THEN
542 donestepping=.true.
543 ELSE
544 donestepping=.false.
545 END IF
546# else
547 donestepping=.true.
548# endif
549!
550 IF (donestepping) THEN
551 DO ng=1,ngrids
552 ldefqck(ng)=.false.
553 lwrtqck(ng)=.false.
554# ifdef AVERAGES
555 ldefavg(ng)=.false.
556 lwrtavg(ng)=.false.
557# endif
558# ifdef DIAGNOSTICS
559 ldefdia(ng)=.false.
560 lwrtdia(ng)=.false.
561# endif
562 wrtnlmod(ng)=.false.
563 wrtobsscale(ng)=.false.
564 END DO
565!
566! Report data penalty function.
567!
568 DO ng=1,ngrids
569 IF (master) THEN
570 DO i=0,nobsvar(ng)
571 IF (i.eq.0) THEN
572 string='Total'
573 ELSE
574 string=obsname(i)
575 END IF
576 IF (fourdvar(ng)%NLPenalty(i).ne.0.0_r8) THEN
577 WRITE (stdout,30) my_outer, 0, 'NLM', &
578 & fourdvar(ng)%NLPenalty(i), &
579 & trim(string)
580 END IF
581 END DO
582 END IF
583!
584! Write out initial data penalty function to NetCDF file.
585!
586 sourcefile=myfile
587 SELECT CASE (dav(ng)%IOtype)
588 CASE (io_nf90)
589 CALL netcdf_put_fvar (ng, inlm, dav(ng)%name, &
590 & 'NL_iDataPenalty', &
591 & fourdvar(ng)%NLPenalty(0:), &
592 & (/1/), (/nobsvar(ng)+1/), &
593 & ncid = dav(ng)%ncid)
594
595# if defined PIO_LIB && defined DISTRIBUTE
596 CASE (io_pio)
597 CALL pio_netcdf_put_fvar (ng, inlm, dav(ng)%name, &
598 & 'NL_iDataPenalty', &
599 & fourdvar(ng)%NLPenalty(0:), &
600 & (/1/), (/nobsvar(ng)+1/), &
601 & piofile = dav(ng)%pioFile)
602# endif
603 END SELECT
604 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
605!
606! Clean penalty array before next run of NL model.
607!
608 fourdvar(ng)%NLPenalty=0.0_r8
609 END DO
610 END IF
611
612# ifdef PROFILE
613!
614! Stop profile clock
615!
616 DO ng=1,ngrids
617 DO thread=thread_range
618 CALL wclock_off (ng, inlm, 86, __line__, myfile)
619 END DO
620 END DO
621# endif
622!
623 10 FORMAT (a,'_outer',i0,'.nc')
624 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
625 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
626 ', TimeSteps: ',i0,' - ',i0,')',/)
627 30 FORMAT (' (',i3.3,',',i3.3,'): ',a,' data penalty, Jdata = ', &
628 & 1p,e17.10,0p,t68,a)
629
630 RETURN

References mod_iounits::avg, background_initialize(), mod_iounits::dav, mod_iounits::dia, mod_scalars::dt, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_scalars::frequentimpulse, mod_scalars::iic, mod_scalars::initime, mod_param::inlm, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_scalars::ldefavg, mod_scalars::ldefdia, mod_scalars::ldefqck, mod_scalars::lwrtavg, mod_scalars::lwrtdia, mod_scalars::lwrtqck, main2d(), main3d(), mod_parallel::master, mod_scalars::myruninterval, mod_param::ngrids, mod_fourdvar::nobsvar, mod_scalars::noerror, mod_scalars::ntend, mod_scalars::ntimes, mod_scalars::ntstart, mod_fourdvar::obsname, mod_iounits::sourcefile, mod_scalars::sporadicimpulse, mod_iounits::stdout, mod_scalars::step_counter, mod_scalars::time, wclock_off(), wclock_on(), mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, and mod_fourdvar::wrtobsscale.

Here is the call graph for this function:

◆ background_initialize()

subroutine, public rbl4dvar_mod::background_initialize ( integer, intent(in) my_outer)

Definition at line 211 of file rbl4dvar.F.

212!
213!=======================================================================
214! !
215! This routine initializes ROMS nonlinear model trajectory Xb_n-1(0) !
216! for each (n-1) outer loops. It is separated from the 'background' !
217! 4D-Var phase in ESM coupling applications that use generic methods !
218! for 'initialize', 'run', and 'finalize'. !
219! !
220! On Input: !
221! !
222! my_outer Outer-loop counter (integer) !
223! !
224!=======================================================================
225!
226! Imported variable declarations
227!
228 integer, intent(in) :: my_outer
229!
230! Local variable declarations.
231!
232 integer :: lstr, ng, tile
233# ifdef STD_MODEL
234 integer :: Lstd
235# endif
236
237 integer :: Fcount, Tindex
238# ifdef PROFILE
239 integer :: thread
240# endif
241!
242 character (len=*), parameter :: MyFile = &
243 & __FILE__//", background_initialize"
244!
245 sourcefile=myfile
246!
247!-----------------------------------------------------------------------
248! Initialize nonlinear model kernel.
249!-----------------------------------------------------------------------
250
251# ifdef PROFILE
252!
253! Start profile clock.
254!
255 DO ng=1,ngrids
256 DO thread=thread_range
257 CALL wclock_on (ng, inlm, 86, __line__, myfile)
258 END DO
259 END DO
260# endif
261!
262! Initialize the switch to gather weak constraint forcing.
263!
264 DO ng=1,ngrids
265 wrtforce(ng)=.false.
266 END DO
267!
268! Clear TLM mixing arrays.
269!
270 DO ng=1,ngrids
271 DO tile=first_tile(ng),last_tile(ng),+1
272 CALL initialize_mixing (ng, tile, itlm)
273 END DO
274 END DO
275!
276! Get nonlinear model initial conditions.
277!
278 DO ng=1,ngrids
279 wrtnlmod(ng)=.true.
280 wrttlmod(ng)=.false.
281# ifdef FORWARD_FLUXES
282 lreadblk(ng)=.false.
283# endif
284# ifdef FRC_FILES
285 lreadfrc(ng)=.true.
286# endif
287 lreadfwd(ng)=.false.
288 rst(ng)%Rindex=0
289 fcount=rst(ng)%load
290 rst(ng)%Nrec(fcount)=0
291 END DO
292!
293 CALL initial
294 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
295!
296! Save nonlinear initial conditions (currently in time index 1,
297! background) into record "Lbck" of INI(ng)%name NetCDF file. The
298! record "Lbck" becomes the background state and "Lini" becomes
299! the current nonlinear initial conditions.
300!
301 tindex=1
302 DO ng=1,ngrids
303 ini(ng)%Rindex=1
304 fcount=ini(ng)%load
305 ini(ng)%Nrec(fcount)=1
306# ifdef DISTRIBUTE
307 CALL wrt_ini (ng, myrank, tindex, lbck)
308# else
309 CALL wrt_ini (ng, -1, tindex, lbck)
310# endif
311 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
312 END DO
313
314# ifdef STD_MODEL
315!
316! Estimate initial conditions standard deviation fields for the
317! background error covariance using the method of Mogensen et al.
318! (2012), which assumes that the background errors are proportional
319! to the vertical derivatives of the state. The error have simmilar
320! shape as the prior profile but iwth a displacement over the water
321! column.
322!
323 IF (my_outer.eq.0) THEN
324 lstd=1 ! time index in error arrays
325 DO ng=1,ngrids
326 DO tile=first_tile(ng),last_tile(ng),+1
327 CALL background_std (ng, tile, tindex, lstd)
328 END DO
329 END DO
330!
331! Write out estimated background standard deviation fields.
332!
333 DO ng=1,ngrids
334 IF (lwrtstd(ng)) THEN
335 std(5,ng)%Rindex=0
336 fcount=std(5,ng)%load
337 std(5,ng)%Nrec(fcount)=1
338# ifdef DISTRIBUTE
339 CALL wrt_std (ng, myrank, lstd)
340# else
341 CALL wrt_std (ng, -1, lstd)
342# endif
343 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
344 END IF
345 END DO
346 END IF
347# endif
348!
349! Set nonlinear output history file as the initial basic state
350! trajectory.
351!
352 DO ng=1,ngrids
353 ldefhis(ng)=.true.
354 lwrthis(ng)=.true.
355 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer
356 lstr=len_trim(his(ng)%name)
357 his(ng)%base=his(ng)%name(1:lstr-3)
358 END DO
359!
360! Set the nonlinear model to output the quicksave history file as a
361! function of the outer loop. It may be used as the basic state
362! trajectory for the surface fluxes (wind stress, shortwave, heat
363! flux, and E-P) because they can be saved frequently to resolve
364! the daily cycle while avoiding large files.
365!
366 DO ng=1,ngrids
367 ldefqck(ng)=.true.
368 lwrtqck(ng)=.true.
369 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer
370 lstr=len_trim(qck(ng)%name)
371 qck(ng)%base=qck(ng)%name(1:lstr-3)
372 END DO
373!
374! Define output 4D-Var NetCDF file (DAV struc) containing all
375! processed data at observation locations.
376!
377 DO ng=1,ngrids
378 ldefmod(ng)=.true.
379 CALL def_mod (ng)
380 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
381 END DO
382
383# ifdef PROFILE
384!
385! Stop profile clock
386!
387 DO ng=1,ngrids
388 DO thread=thread_range
389 CALL wclock_off (ng, inlm, 86, __line__, myfile)
390 END DO
391 END DO
392# endif
393!
394 10 FORMAT (a,'_outer',i0,'.nc')
395!
396 RETURN

References def_mod_mod::def_mod(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_iounits::fwd, mod_iounits::his, mod_iounits::ini, initial(), mod_mixing::initialize_mixing(), mod_param::inlm, mod_param::itlm, mod_parallel::last_tile, lbck, mod_scalars::ldefhis, mod_scalars::ldefmod, mod_scalars::ldefqck, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, mod_scalars::lwrthis, mod_scalars::lwrtqck, mod_parallel::myrank, mod_param::ngrids, mod_scalars::noerror, mod_iounits::qck, mod_iounits::rst, mod_iounits::sourcefile, mod_iounits::std, wclock_off(), wclock_on(), wrt_ini_mod::wrt_ini(), mod_fourdvar::wrtforce, mod_fourdvar::wrtnlmod, and mod_fourdvar::wrttlmod.

Referenced by background().

Here is the call graph for this function:
Here is the caller graph for this function:

◆ increment()

subroutine, public rbl4dvar_mod::increment ( integer, intent(in) my_outer,
real(dp), intent(in) runinterval )

Definition at line 633 of file rbl4dvar.F.

634!
635!=======================================================================
636! !
637! This routine computes the 4D-Var data assimilation state increment, !
638! dXa, by iterating the inner loops and minimizing the cost function. !
639! !
640! On Input: !
641! !
642! my_outer Outer-loop counter (integer) !
643! RunInterval TLM/ADM kernels time stepping window (seconds) !
644! !
645!=======================================================================
646
647# ifdef RPCG
648!
649 USE comp_jb0_mod, ONLY : aug_oper, comp_jb0
650 USE ini_adjust_mod, ONLY : ini_adjust
651 USE sum_grad_mod, ONLY : sum_grad
652 USE sum_imp_mod, ONLY : sum_imp
653# endif
654!
655! Imported variable declarations
656!
657 integer, intent(in) :: my_outer
658!
659 real(dp), intent(in) :: RunInterval
660!
661! Local variable declarations.
662!
663 logical :: Lcgini, Linner, Lposterior
664!
665 integer :: i, ifile, lstr, my_inner, ng, tile
666 integer :: Fcount, InpRec
667# ifdef RPCG
668 integer :: ADrec, nLAST
669 integer :: irec, jrec, jrec1, jrec2
670 integer :: LiNL
671# endif
672# ifdef PROFILE
673 integer :: thread
674# endif
675!
676 integer, dimension(Ngrids) :: Nrec
677# ifdef RPCG
678 integer, dimension(Ngrids) :: nADrec
679# endif
680# ifdef SPLIT_4DVAR
681!
682 real(dp) :: stime
683# endif
684!
685 character (len=8 ) :: driver = 'rbl4dvar'
686 character (len=10) :: suffix
687
688 character (len=*), parameter :: MyFile = &
689 & __FILE__//", increment"
690!
691 sourcefile=myfile
692!
693!=======================================================================
694! Compute 4D-Var increment.
695!=======================================================================
696
697# ifdef PROFILE
698!
699! Start profile clock.
700!
701 DO ng=1,ngrids
702 DO thread=thread_range
703 CALL wclock_on (ng, itlm, 87, __line__, myfile)
704 END DO
705 END DO
706# endif
707
708# ifdef SPLIT_4DVAR
709!
710!-----------------------------------------------------------------------
711! If split 4D-Var algorithm, set several variables that are computed
712! or assigned in other 4D-Var phase executable.
713!-----------------------------------------------------------------------
714!
715! Reset Nrun counter to a value greater than one.
716!
717 IF (my_outer.gt.1) THEN
718 nrun=1+(my_outer-1)*ninner
719 END IF
720!
721! Open Nonlinear model initial conditions NetCDF file (INI struc) and
722! inquire about its variables IDs.
723!
724 DO ng=1,ngrids
725 ldefini(ng)=.false.
726 CALL def_ini (ng)
727 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
728 END DO
729!
730! If outer>1, open tangent linear initial conditions NetCDF file
731! (ITL struc) and inquire about its variables IDs.
732!
733 IF (my_outer.gt.1) THEN
734 DO ng=1,ngrids
735 ldefitl(ng)=.false.
736 CALL tl_def_ini (ng)
737 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
738 END DO
739 END IF
740!
741! If outer>1, open impulse forcing NetCDF file (TLF struc) and inquire
742! about its variables IDs.
743!
744 IF (my_outer.gt.1) THEN
745 DO ng=1,ngrids
746 ldeftlf(ng)=.false.
747 CALL def_impulse (ng)
748 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
749 END DO
750 END IF
751!
752! If outer>1, open adjoint history NetCDF file (ADM struc) and inquire
753! about its variables IDs and set "FrcRec".
754!
755 IF (my_outer.gt.1) THEN
756 DO ng=1,ngrids
757 ldefadj(ng)=.false.
758 CALL ad_def_his (ng, ldefadj(ng))
759 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
760 fcount=adm(ng)%Fcount
761 frcrec(ng)=adm(ng)%Nrec(fcount)
762 END DO
763 END IF
764!
765! Open 4D-Var NetCDF file (DAV struc) and inquire about its variables.
766! Deactivate "haveNLmod" since its values are read below. Otherwise,
767! the TLM will read zero values when calling "obs_read" for outer>1.
768! The DAV file does not have values for NLmodel_value(:,outer) yet.
769!
770 DO ng=1,ngrids
771 ldefmod(ng)=.false.
772 CALL def_mod (ng)
773 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
774 havenlmod(ng)=.false. ! because is activated in "def_mod"
775 END DO
776!
777! If split 4D-Var and outer>1, read several variables from 4D-VAR
778! NetCDF (DAV struc) file needed in the conjugate gradient algorithm.
779! In the unsplit case, such values are available in memory.
780!
781 IF (my_outer.gt.1) THEN
782 DO ng=1,ngrids
783# ifdef RPCG
784 CALL cg_read_rpcg (ng, itlm, my_outer)
785# else
786 CALL cg_read_congrad (ng, itlm, my_outer)
787# endif
788 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
789 END DO
790 END IF
791!
792! In the split 4D-Var algorithm, we need to read in initial NLmodVal
793! ("NLmodel_initial") and ObsScale ("obs_scale"), computed in the
794! background phase, from the DAV NetCDF file. Such data is in memory
795! in the unsplit algorithm.
796!
797! HGA: What to do in 4D-Var with nested grids?
798!
799 IF (my_outer.eq.1) THEN ! initial values from "background"
800 ng=1
801 SELECT CASE (dav(ng)%IOtype)
802 CASE (io_nf90)
803 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
804 & vname(1,idnlmi), nlmodval, &
805 & ncid = dav(ng)%ncid)
806
807# if defined PIO_LIB && defined DISTRIBUTE
808 CASE (io_pio)
809 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
810 & vname(1,idnlmi), nlmodval, &
811 & piofile = dav(ng)%pioFile)
812
813# endif
814 END SELECT
815 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
816 END IF
817!
818 ng=1
819 SELECT CASE (dav(ng)%IOtype)
820 CASE (io_nf90)
821 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
822 & vname(1,idobss), obsscale, &
823 & ncid = dav(ng)%ncid)
824
825# if defined PIO_LIB && defined DISTRIBUTE
826 CASE (io_pio)
827 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
828 & vname(1,idobss), obsscale, &
829 & piofile = dav(ng)%pioFile)
830# endif
831 END SELECT
832 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
833!
834! In the split 4D-Var alorithm, we also need to read in ObsVal
835! ("obs_value") and ObsErr ("obs_error") from the observation file.
836! They are used in the initialization of the conjugate gradient
837! algorithm.
838!
839 ng=1
840 SELECT CASE (obs(ng)%IOtype)
841 CASE (io_nf90)
842 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
843 & vname(1,idoval), obsval)
844 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
845
846 CALL netcdf_get_fvar (ng, itlm, obs(ng)%name, &
847 & vname(1,idoerr), obserr)
848 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
849
850# if defined PIO_LIB && defined DISTRIBUTE
851 CASE (io_pio)
852 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
853 & vname(1,idoval), obsval)
854 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
855
856 CALL pio_netcdf_get_fvar (ng, itlm, obs(ng)%name, &
857 & vname(1,idoerr), obserr)
858 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
859# endif
860 END SELECT
861!
862! Set nonlinear output history file to be used as the basic state
863! trajectory. The 4D-Var increment phase is computed by a different
864! executable and needs to know some of the HIS structure information.
865!
866 DO ng=1,ngrids
867 WRITE (his(ng)%name,10) trim(fwd(ng)%head), my_outer-1
868 lstr=len_trim(his(ng)%name)
869 his(ng)%base=his(ng)%name(1:lstr-3)
870 IF (his(ng)%Nfiles.gt.1) THEN
871 DO ifile=1,his(ng)%Nfiles
872 WRITE (suffix,"('_',i4.4,'.nc')") ifile
873 his(ng)%files(ifile)=trim(his(ng)%base)//trim(suffix)
874 END DO
875 his(ng)%name=trim(his(ng)%files(1))
876 ELSE
877 his(ng)%files(1)=trim(his(ng)%name)
878 END IF
879 END DO
880!
881! Set the nonlinear model quicksave-history file as the basic state for
882! the surface fluxes computed in "bulk_flux", which may be available at
883! more frequent intervals while avoiding large files. Since the 4D-Var
884! increment phase is calculated by a different executable and needs to
885! know some of the QCK structure information.
886!
887 DO ng=1,ngrids
888 WRITE (qck(ng)%name,10) trim(qck(ng)%head), my_outer-1
889 lstr=len_trim(qck(ng)%name)
890 qck(ng)%base=qck(ng)%name(1:lstr-3)
891 IF (qck(ng)%Nfiles.gt.1) THEN
892 DO ifile=1,qck(ng)%Nfiles
893 WRITE (suffix,"('_',i4.4,'.nc')") ifile
894 qck(ng)%files(ifile)=trim(qck(ng)%base)//trim(suffix)
895 END DO
896 qck(ng)%name=trim(qck(ng)%files(1))
897 ELSE
898 qck(ng)%files(1)=trim(qck(ng)%name)
899 END IF
900 END DO
901!
902! Read in 4D-Var starting time (sec) from nonlinear trajectory.
903! Initialize "tday" which are needed to write the correct time in
904! the ITL NetCDF file. It is alse need for boundary and surface
905! forcing adjustments, if any.
906!
907 inprec=1
908 DO ng=1,ngrids
909 SELECT CASE (his(ng)%IOtype)
910 CASE (io_nf90)
911 CALL netcdf_get_fvar (ng, itlm, his(ng)%name, &
912 & vname(1,idtime), stime, &
913 & start = (/inprec/), &
914 & total = (/1/))
915
916# if defined PIO_LIB && defined DISTRIBUTE
917 CASE (io_pio)
918 CALL pio_netcdf_get_fvar (ng, itlm, his(ng)%name, &
919 & vname(1,idtime), stime, &
920 & start = (/inprec/), &
921 & total = (/1/))
922# endif
923 END SELECT
924 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
925 initime(ng)=stime
926
927# ifdef ADJUST_BOUNDARY
928!
929! Set time (sec) for the open boundary adjustment.
930!
931 obc_time(1,ng)=stime
932 DO i=2,nbrec(ng)
933 obc_time(i,ng)=obc_time(i-1,ng)+nobc(ng)*dt(ng)
934 END DO
935# endif
936
937# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
938!
939! Set time (sec) for the surface forcing adjustment.
940!
941 sf_time(1,ng)=stime
942 DO i=2,nfrec(ng)
943 sf_time(i,ng)=sf_time(i-1,ng)+nsff(ng)*dt(ng)
944 END DO
945# endif
946 END DO
947# endif
948!
949!-----------------------------------------------------------------------
950! Set few variables needed in the "increment" phase.
951!-----------------------------------------------------------------------
952!
953! Set structure for the nonlinear forward trajectory to be processed
954! by the tangent linear and adjoint models. Also, set switches to
955! process the FWD structure in routine "check_multifile". Notice that
956! it is possible to split solution into multiple NetCDF files to reduce
957! their size.
958!
959 CALL edit_multifile ('HIS2FWD')
960 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
961 DO ng=1,ngrids
962 lreadfwd(ng)=.true.
963 END DO
964
965# ifdef FORWARD_FLUXES
966!
967! Set the BLK structure to contain the nonlinear model surface fluxes
968! needed by the tangent linear and adjoint models. Also, set switches
969! to process that structure in routine "check_multifile". Notice that
970! it is possible to split the solution into multiple NetCDF files to
971! reduce their size.
972!
973! The switch LreadFRC is deactivated because all the atmospheric
974! forcing, including shortwave radiation, is read from the NLM
975! surface fluxes or is assigned during ESM coupling. Such fluxes
976! are available from the QCK structure. There is no need for reading
977! and processing from the FRC structure input forcing-files.
978!
979 CALL edit_multifile ('QCK2BLK')
980 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
981 DO ng=1,ngrids
982 lreadblk(ng)=.true.
983 lreadfrc(ng)=.false.
984 lreadqck(ng)=.false.
985 END DO
986# endif
987
988# ifdef RPCG
989!
990! If weak constraint, set the number of TLF records used in the
991! augmented solver.
992!
993 DO ng=1,ngrids
994 IF (nadj(ng).lt.ntimes(ng)) THEN
995 nadrec(ng)=2*(1+ntimes(ng)/nadj(ng))
996 ELSE
997 nadrec(ng)=0
998 END IF
999 END DO
1000# endif
1001!
1002!-----------------------------------------------------------------------
1003! On first pass (outer=1), create NetCDF files and initialize all the
1004! records in the ITL file to zero.
1005!-----------------------------------------------------------------------
1006!
1007 check_outer1 : IF (my_outer.eq.1) THEN
1008!
1009! If outer=1, define tangent linear initial conditions file.
1010!
1011 DO ng=1,ngrids
1012 ldefitl(ng)=.true.
1013 CALL tl_def_ini (ng)
1014 ldefitl(ng)=.false.
1015 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1016 END DO
1017
1018# ifdef RPCG
1019!
1020! If outer=1, Initialize all records of the ITL file to zero. The TLM
1021! state variables are zero when outer=1.
1022!
1023 DO ng=1,ngrids
1024 tdays(ng)=initime(ng)*sec2day
1025# ifdef DISTRIBUTE
1026 CALL tl_wrt_ini (ng, myrank, rec1, rec1)
1027# else
1028 CALL tl_wrt_ini (ng, -1, rec1, rec1)
1029# endif
1030 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1031# ifdef DISTRIBUTE
1032 CALL tl_wrt_ini (ng, myrank, rec1, rec2)
1033# else
1034 CALL tl_wrt_ini (ng, -1, rec1, rec2)
1035# endif
1036 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1037# ifdef DISTRIBUTE
1038 CALL tl_wrt_ini (ng, myrank, rec1, rec3)
1039# else
1040 CALL tl_wrt_ini (ng, -1, rec1, rec3)
1041# endif
1042 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1043# ifdef DISTRIBUTE
1044 CALL tl_wrt_ini (ng, myrank, rec1, rec4)
1045# else
1046 CALL tl_wrt_ini (ng, -1, rec1, rec4)
1047# endif
1048 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1049# ifdef DISTRIBUTE
1050 CALL tl_wrt_ini (ng, myrank, rec1, rec5)
1051# else
1052 CALL tl_wrt_ini (ng, -1, rec1, rec5)
1053# endif
1054 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1055!
1056 IF (nadj(ng).lt.ntimes(ng)) THEN
1057 nlast=rec5
1058 DO irec=1,nadrec(ng)
1059# ifdef DISTRIBUTE
1060 CALL tl_wrt_ini (ng, myrank, rec1, nlast+irec)
1061# else
1062 CALL tl_wrt_ini (ng, -1, rec1, nlast+irec)
1063# endif
1064 IF (founderror(exit_flag, noerror, &
1065 & __line__, myfile)) RETURN
1066 END DO
1067 END IF
1068 END DO
1069# endif
1070!
1071! If outer=1, define impulse forcing NetCDF file.
1072!
1073 DO ng=1,ngrids
1074 ldeftlf(ng)=.true.
1075 CALL def_impulse (ng)
1076 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1077 END DO
1078
1079# if defined POSTERIOR_EOFS || defined POSTERIOR_ERROR_I || \
1080 defined posterior_error_f
1081!
1082! If outer=1, define output Hessian NetCDF file that will eventually
1083! contain the intermediate posterior analysis error covariance matrix
1084! fields or its EOFs.
1085!
1086 DO ng=1,ngrids
1087 ldefhss(ng)=.true.
1088 CALL def_hessian (ng)
1089 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1090 END DO
1091# endif
1092
1093# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
1094!
1095! If outer=1, define output initial or final full posterior error
1096! covariance (diagonal) matrix NetCDF.
1097!
1098 DO ng=1,ngrids
1099 ldeferr(ng)=.true.
1100 CALL def_error (ng)
1101 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1102 END DO
1103# endif
1104
1105 END IF check_outer1
1106
1107# ifdef RPCG
1108!
1109!-----------------------------------------------------------------------
1110! If outer>1, integrate the tangent linear model to evaluate
1111! G[x(k)-x(k-1)], which is written to the DAV NetCDF in the
1112! variable "TLmodel_value". It is needed in routine "rpcg_lanczos".
1113! In the old unsplit algorithm, the TLM was run in the "analysis"
1114! phase. It is done here to allow lower resolution in the inner
1115! loops in the future.
1116!-----------------------------------------------------------------------
1117!
1118 check_outer2 : IF (my_outer.gt.1) THEN
1119!
1120! Clear tangent linear forcing arrays. This is very important since
1121! these arrays are non-zero and must be zero when running the tangent
1122! linear model. We also need to clear the nonlinear state arrays to
1123! make sure that the RHS terms rzeta, rubar, rvbar, ru, and rv are
1124! zero. Otherwise, those array will have the last computed values
1125! when running the nonlinear model if not processing the forward
1126! trajectory RHS terms. This needs to be done to get identical
1127! solutions with the split schemes.
1128!
1129 DO ng=1,ngrids
1130 DO tile=first_tile(ng),last_tile(ng),+1
1131 CALL initialize_ocean (ng, tile, inlm)
1132 CALL initialize_forces (ng, tile, itlm)
1133# ifdef ADJUST_BOUNDARY
1134 CALL initialize_boundary (ng, tile, itlm)
1135# endif
1136 END DO
1137 END DO
1138!
1139! Initialize tangent linear model from initial impulse which is now
1140! stored in file ITL(ng)%name.
1141!
1142 DO ng=1,ngrids
1143 wrtnlmod(ng)=.false.
1144 wrttlmod(ng)=.true.
1145 END DO
1146!
1147! If weak constraint, the impulses are time-interpolated at each
1148! time-steps.
1149!
1150 DO ng=1,ngrids
1151 IF (frcrec(ng).gt.3) THEN
1152 frequentimpulse(ng)=.true.
1153 END IF
1154 END DO
1155!
1156! Initialize tangent linear model from ITLname, record Rec3.
1157! The sum of the initial condition increments is in record 3.
1158!
1159 DO ng=1,ngrids
1160 itl(ng)%Rindex=rec3
1161 CALL tl_initial (ng)
1162 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1163 END DO
1164!
1165! Run tangent linear model.
1166!
1167 DO ng=1,ngrids
1168 IF (master) THEN
1169 WRITE (stdout,20) 'TL', ng, my_outer, 0, &
1170 & ntstart(ng), ntend(ng)
1171 END IF
1172 END DO
1173!
1174# ifdef SOLVE3D
1175 CALL tl_main3d (runinterval)
1176# else
1177 CALL tl_main2d (runinterval)
1178# endif
1179 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1180!
1181 DO ng=1,ngrids
1182 wrtnlmod(ng)=.false.
1183 wrttlmod(ng)=.false.
1184 END DO
1185 END IF check_outer2
1186# endif
1187!
1188!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1189! 4D-Var inner loops.
1190!:::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::::
1191!
1192! Clear tangent linear forcing arrays before entering inner-loop.
1193! This is very important since these arrays are non-zero and must
1194! be zero when running the tangent linear model.
1195!
1196 DO ng=1,ngrids
1197 DO tile=first_tile(ng),last_tile(ng),+1
1198 CALL initialize_ocean (ng, tile, inlm)
1199 CALL initialize_forces (ng, tile, itlm)
1200# ifdef ADJUST_BOUNDARY
1201 CALL initialize_boundary (ng, tile, itlm)
1202# endif
1203 END DO
1204 END DO
1205
1206# if defined BALANCE_OPERATOR && defined ZETA_ELLIPTIC
1207!
1208! Compute the reference zeta and biconjugate gradient arrays
1209! required for the balance of free surface.
1210!
1211 IF (balance(isfsur)) THEN
1212 DO ng=1,ngrids
1213 CALL get_state (ng, inlm, 2, ini(ng), lini, lini)
1214 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1215!
1216 DO tile=first_tile(ng),last_tile(ng),+1
1217 CALL balance_ref (ng, tile, lini)
1218 CALL biconj (ng, tile, inlm, lini)
1219 END DO
1220 wrtzetaref(ng)=.true.
1221 END DO
1222 END IF
1223# endif
1224!
1225! Inner loops iteration.
1226!
1227 inner_loop : DO my_inner=0,ninner
1228 inner=my_inner
1229
1230# ifdef RPCG
1231 IF (inner.ne.ninner) THEN
1232 linner=.true.
1233 ELSE
1234 linner=.false.
1235 END IF
1236!
1237! Retrieve TLmodVal and NLmodVal when inner=0 and outer>1 for use as
1238! Hbk and BCKmodVal, respectively.
1239!
1240 IF ((inner.eq.0).and.(my_outer.gt.1)) THEN
1241 DO ng=1,ngrids
1242 SELECT CASE (dav(ng)%IOtype)
1243 CASE (io_nf90)
1244 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1245 & 'TLmodel_value', tlmodval)
1246 IF (founderror(exit_flag, noerror, __line__, &
1247 & myfile)) RETURN
1248
1249 CALL netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1250 & 'NLmodel_value', nlmodval, &
1251 & ncid = dav(ng)%ncid, &
1252 & start = (/1,my_outer-1/), &
1253 & total = (/ndatum(ng),1/))
1254 IF (founderror(exit_flag, noerror, __line__, &
1255 & myfile)) RETURN
1256
1257# if defined PIO_LIB && defined DISTRIBUTE
1258 CASE (io_pio)
1259 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1260 & 'TLmodel_value', tlmodval)
1261 IF (founderror(exit_flag, noerror, __line__, &
1262 & myfile)) RETURN
1263
1264 CALL pio_netcdf_get_fvar (ng, itlm, dav(ng)%name, &
1265 & 'NLmodel_value', nlmodval, &
1266 & piofile = dav(ng)%pioFile, &
1267 & start = (/1,my_outer-1/), &
1268 & total = (/ndatum(ng),1/))
1269 IF (founderror(exit_flag, noerror, __line__, &
1270 & myfile)) RETURN
1271# endif
1272 END SELECT
1273 END DO
1274 END IF
1275!
1276 IF (inner.eq.0) THEN
1277 lcgini=.true.
1278 END IF
1279 DO ng=1,ngrids
1280 CALL rpcg_lanczos (ng, itlm, my_outer, inner, ninner, lcgini)
1281 END DO
1282# else
1283!
1284! Initialize conjugate gradient algorithm depending on hot start or
1285! outer loop index.
1286!
1287 IF (inner.eq.0) THEN
1288 lcgini=.true.
1289 DO ng=1,ngrids
1290 CALL congrad (ng, irpm, my_outer, inner, ninner, lcgini)
1291 END DO
1292 END IF
1293!
1294! If initialization step, skip the inner-loop computations.
1295!
1296 linner=.false.
1297 IF ((inner.ne.0).or.(nrun.ne.1)) THEN
1298 IF (((inner.eq.0).and.lhotstart).or.(inner.ne.0)) THEN
1299 linner=.true.
1300 END IF
1301 END IF
1302# endif
1303!
1304! Start inner loop computations.
1305!
1306 inner_compute : IF (linner) THEN
1307!
1308!-----------------------------------------------------------------------
1309! Integrate adjoint model forced with any vector PSI at the observation
1310! locations and generate adjoint trajectory, Lambda_n(t).
1311!-----------------------------------------------------------------------
1312!
1313! Initialize the adjoint model from rest.
1314!
1315 DO ng=1,ngrids
1316 CALL ad_initial (ng)
1317 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1318 wrtmisfit(ng)=.false.
1319 END DO
1320!
1321! Set adjoint history NetCDF parameters. Define adjoint history
1322! file only once to avoid opening too many files.
1323!
1324 DO ng=1,ngrids
1325# ifdef WEAK_NOINTERP
1326 wrtforce(ng)=.false.
1327# else
1328 wrtforce(ng)=.true.
1329# endif
1330 IF (nrun.gt.1) ldefadj(ng)=.false.
1331 fcount=adm(ng)%load
1332 adm(ng)%Nrec(fcount)=0
1333 adm(ng)%Rindex=0
1334 END DO
1335!
1336! Time-step adjoint model backwards forced with current PSI vector.
1337!
1338 DO ng=1,ngrids
1339 IF (master) THEN
1340 WRITE (stdout,20) 'AD', ng, my_outer, my_inner, &
1341 & ntstart(ng), ntend(ng)
1342 END IF
1343 END DO
1344!
1345# ifdef SOLVE3D
1346 CALL ad_main3d (runinterval)
1347# else
1348 CALL ad_main2d (runinterval)
1349# endif
1350 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1351!
1352! Activate "LwrtState2d" switch to write the correct ad_ubar, ad_vbar,
1353! ad_zeta here instead of ad_ubar_sol, ad_vbar_sol, and ad_zeta_sol in
1354! the calls to "ad_wrt_his" below (HGA).
1355! Here, ad_zeta(:,:,kout)=ad_zeta_sol. However, ad_ubar and ad_vbar
1356! not equal to ad_ubar_sol and ad_vbar_sol, respectively. It is
1357! irrelevant because ubar and vbar are not part of the state in
1358! 3D application. Notice that "LwrtState2d" will be turned off at
1359! the bottom of "error_covariance".
1360!
1361 DO ng=1,ngrids
1362 lwrtstate2d(ng)=.true.
1363 END DO
1364
1365# ifndef WEAK_NOINTERP
1366!
1367! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
1368! record into the adjoint history file. Note that the weak-constraint
1369! forcing is delayed by nADJ time-steps.
1370!
1371 DO ng=1,ngrids
1372# ifdef DISTRIBUTE
1373 CALL ad_wrt_his (ng, myrank)
1374# else
1375 CALL ad_wrt_his (ng, -1)
1376# endif
1377 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1378 END DO
1379# endif
1380!
1381! Write out adjoint initial condition record into the adjoint
1382! history file.
1383!
1384 DO ng=1,ngrids
1385 wrtforce(ng)=.false.
1386# ifdef DISTRIBUTE
1387 CALL ad_wrt_his (ng, myrank)
1388# else
1389 CALL ad_wrt_his (ng, -1)
1390# endif
1391 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1392 END DO
1393!
1394! Convolve adjoint trajectory with error covariances.
1395!
1396# ifdef POSTERIOR_ERROR_I
1397 lposterior=.true.
1398# else
1399 lposterior=.false.
1400# endif
1401# ifdef RPCG
1402!
1403! Set the flag that controls the augmentation of the model error
1404! forcing terms. This is ONLY done in the outer-loop so
1405! LaugWeak=.FALSE. here.
1406!
1407 laugweak=.false.
1408# endif
1409 CALL error_covariance (itlm, driver, my_outer, inner, &
1410 & lbck, lini, lold, lnew, &
1411 & rec1, rec2, lposterior)
1412 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1413!
1414! Convert the current adjoint solution in ADM(ng)%name to impulse
1415! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
1416! To facilitate the forcing to the TLM and RPM, the forcing is
1417! processed and written in increasing time coordinates (recall that
1418! the adjoint solution in ADM(ng)%name is backwards in time).
1419!
1420 IF (master) THEN
1421 WRITE (stdout,30) my_outer, inner
1422 END IF
1423 DO ng=1,ngrids
1424 tlf(ng)%Rindex=0
1425# ifdef DISTRIBUTE
1426 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1427# else
1428 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1429# endif
1430 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1431 END DO
1432!
1433!-----------------------------------------------------------------------
1434! Integrate tangent linear model forced by the convolved adjoint
1435! trajectory (impulse forcing) to compute R_n * PSI at observation
1436! points.
1437!-----------------------------------------------------------------------
1438!
1439! Initialize tangent linear model from initial impulse which is now
1440! stored in file ITL(ng)%name.
1441!
1442 DO ng=1,ngrids
1443 wrtnlmod(ng)=.false.
1444 wrttlmod(ng)=.true.
1445 END DO
1446!
1447! If weak constraint, the impulses are time-interpolated at each
1448! time-steps.
1449!
1450 DO ng=1,ngrids
1451 IF (frcrec(ng).gt.3) THEN
1452 frequentimpulse(ng)=.true.
1453 END IF
1454 END DO
1455!
1456! Initialize tangent linear model from ITL(ng)%name, record Rec1.
1457!
1458 DO ng=1,ngrids
1459 itl(ng)%Rindex=rec1
1460 CALL tl_initial (ng)
1461 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1462 END DO
1463!
1464! Activate switch to write out initial misfit between model and
1465! observations.
1466!
1467 IF ((my_outer.eq.1).and.(inner.eq.1)) THEN
1468 DO ng=1,ngrids
1469 wrtmisfit(ng)=.true.
1470 END DO
1471 END IF
1472!
1473! Run tangent linear model forward and force with convolved adjoint
1474! trajectory impulses. Compute (H M B M' H')_n * PSI at observation
1475! points which are used in the conjugate gradient algorithm.
1476!
1477 DO ng=1,ngrids
1478 IF (master) THEN
1479 WRITE (stdout,20) 'TL', ng, my_outer, my_inner, &
1480 & ntstart(ng), ntend(ng)
1481 END IF
1482 END DO
1483!
1484# ifdef SOLVE3D
1485 CALL tl_main3d (runinterval)
1486# else
1487 CALL tl_main2d (runinterval)
1488# endif
1489 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1490!
1491 DO ng=1,ngrids
1492 wrtnlmod(ng)=.false.
1493 wrttlmod(ng)=.false.
1494 END DO
1495
1496# ifdef POSTERIOR_ERROR_F
1497!
1498! Copy the final time tl_var(Lold) into ad_var(Lold) so that it can be
1499! written to the Hessian NetCDF file.
1500!
1501 add=.false.
1502 DO ng=1,ngrids
1503 DO tile=first_tile(ng),last_tile(ng),+1
1504 CALL load_tltoad (ng, tile, lold(ng), lold(ng), add)
1505 END DO
1506 END DO
1507!
1508! Write evolved tangent solution into hessian netcdf file for use
1509! later.
1510!
1511 IF (inner.ne.0) THEN
1512 DO ng=1,ngrids
1513# ifdef DISTRIBUTE
1514 CALL wrt_hessian (ng, tile, lold(ng), lold(ng))
1515# else
1516 CALL wrt_hessian (ng, -1, lold(ng), lold(ng))
1517# endif
1518 IF (founderror(exit_flag, noerror, &
1519 & __line__, myfile)) RETURN
1520 END DO
1521 END IF
1522# endif
1523!
1524 nrun=nrun+1
1525
1526# ifndef RPCG
1527!
1528! Use conjugate gradient algorithm to find a better approximation
1529! PSI to coefficients Beta_n.
1530!
1531 DO ng=1,ngrids
1532 lcgini=.false.
1533 CALL congrad (ng, itlm, my_outer, inner, ninner, lcgini)
1534 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1535 END DO
1536# else
1537 lcgini=.false.
1538# endif
1539
1540 END IF inner_compute
1541
1542 END DO inner_loop
1543!
1544!-----------------------------------------------------------------------
1545! Once the w_n, have been approximated with sufficient accuracy,
1546! compute estimates of Lambda_n and Xhat_n by carrying out one
1547! backward intergration of the adjoint model and one forward
1548! itegration of the nonlinear model.
1549!-----------------------------------------------------------------------
1550!
1551! Initialize the adjoint model always from rest.
1552!
1553 DO ng=1,ngrids
1554 CALL ad_initial (ng)
1555 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1556 END DO
1557!
1558! Set adjoint history NetCDF parameters. Define adjoint history
1559! file one to avoid opening to many files.
1560!
1561 DO ng=1,ngrids
1562# ifdef WEAK_NOINTERP
1563 wrtforce(ng)=.false.
1564# else
1565 wrtforce(ng)=.true.
1566# endif
1567 IF (nrun.gt.1) ldefadj(ng)=.false.
1568 fcount=adm(ng)%load
1569 adm(ng)%Nrec(fcount)=0
1570 adm(ng)%Rindex=0
1571 END DO
1572!
1573! Time-step adjoint model backwards forced with estimated coefficients,
1574! Beta_n.
1575!
1576 DO ng=1,ngrids
1577 IF (master) THEN
1578 WRITE (stdout,20) 'AD', ng, my_outer, ninner, &
1579 & ntstart(ng), ntend(ng)
1580 END IF
1581 END DO
1582!
1583# ifdef SOLVE3D
1584 CALL ad_main3d (runinterval)
1585# else
1586 CALL ad_main2d (runinterval)
1587# endif
1588 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1589
1590# ifndef WEAK_NOINTERP
1591!
1592! Write out last weak-constraint forcing (WRTforce is still .TRUE.)
1593! record into the adjoint history file. Note that the weak-constraint
1594! forcing is delayed by nADJ time-steps.
1595!
1596 DO ng=1,ngrids
1597# ifdef DISTRIBUTE
1598 CALL ad_wrt_his (ng, myrank)
1599# else
1600 CALL ad_wrt_his (ng, -1)
1601# endif
1602 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1603 END DO
1604# endif
1605!
1606! Write out adjoint initial condition record into the adjoint
1607! history file.
1608!
1609 DO ng=1,ngrids
1610 wrtforce(ng)=.false.
1611# ifdef DISTRIBUTE
1612 CALL ad_wrt_his (ng, myrank)
1613# else
1614 CALL ad_wrt_his (ng, -1)
1615# endif
1616 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1617 END DO
1618
1619# ifdef RPCG
1620!
1621! Get number of records in adjoint NetCDF.
1622! We need to do this here since ADM(ng)%Nrec is reset to zero in
1623! error_covariance.
1624!
1625 DO ng=1,ngrids
1626 fcount=adm(ng)%load
1627 nrec(ng)=adm(ng)%Nrec(fcount)
1628 END DO
1629# endif
1630!
1631! Convolve adjoint trajectory with error covariances.
1632!
1633 lposterior=.false.
1634
1635# ifdef RPCG
1636!
1637! Set the flag that controls the augmentation of the model error
1638! forcing terms. This is ONLY done in the outer-loop so
1639! LaugWeak=.TRUE. here.
1640!
1641 laugweak=.true.
1642 linl=my_outer+1
1643 CALL error_covariance (inlm, driver, my_outer, inner, &
1644 & linl, lini, lold, lnew, &
1645 & rec1, rec2, lposterior)
1646 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1647# else
1648 CALL error_covariance (inlm, driver, my_outer, inner, &
1649 & lbck, lini, lold, lnew, &
1650 & rec1, rec2, lposterior)
1651 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1652# endif
1653
1654# ifdef RPCG
1655!
1656! Augmented solver:
1657!
1658! NOTES: The ITL file contains 5 records -
1659! Rec2 = the new TL initial condition
1660! Rec3 = the sum of the TL initial conditions
1661! Rec4 = B^-1(xb-xk)=sum_j^k-1 G_j lambda_j
1662! Rec5 = the augmented correction to Rec2.
1663! Rec5+1 to Rec5+nADrec/2 = sum of the TLF forcing increments
1664! Rec5+nADrec/2+1 to Rec5+nADrec = B^-1 of the sum of the TLF
1665! forcing increments.
1666! Reset the flag LaugWeak flag.
1667!
1668 laugweak=.false.
1669
1670# else
1671!
1672! Write out nonlinear model initial conditions into INIname, record
1673! INI(ng)%Rindex.
1674!
1675 DO ng=1,ngrids
1676# ifdef DISTRIBUTE
1677 CALL wrt_ini (ng, myrank, lnew(ng))
1678# else
1679 CALL wrt_ini (ng, -1, lnew(ng))
1680# endif
1681 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1682
1683# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1684 defined adjust_boundary
1685# ifdef DISTRIBUTE
1686 CALL wrt_frc_ad (ng, myrank, lold(ng), ini(ng)%Rindex)
1687# else
1688 CALL wrt_frc_ad (ng, -1, lold(ng), ini(ng)%Rindex)
1689# endif
1690 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1691# endif
1692 END DO
1693# endif
1694
1695# ifdef RPCG
1696!
1697! Compute the augmented correction to the adjoint propagator.
1698! We need to use sum (x(k)-x(k-1)) before it is updated. This is
1699! in record 3 of the ITL file.
1700!
1701 DO ng=1,ngrids
1702 CALL get_state (ng, itlm, 4, itl(ng), rec3, ltlm1)
1703 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1704 END DO
1705 DO ng=1,ngrids
1706 DO tile=first_tile(ng),last_tile(ng),+1
1707 CALL aug_oper (ng, tile, ltlm1, ltlm2)
1708 END DO
1709 END DO
1710!
1711! Save this augmented piece in record 5 of the ITL file.
1712!
1713 DO ng=1,ngrids
1714# ifdef DISTRIBUTE
1715 CALL tl_wrt_ini (ng, myrank, ltlm2, rec5)
1716# else
1717 CALL tl_wrt_ini (ng, -1, ltlm2, rec5)
1718# endif
1719 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1720 END DO
1721!
1722! Complete the computation of the TL initial condition by adding the
1723! contribution from the augmented adjoint propagator.
1724!
1725 DO ng=1,ngrids
1726 CALL get_state (ng, itlm, 4, itl(ng), rec2, ltlm1)
1727 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1728 CALL get_state (ng, itlm, 4, itl(ng), rec5, ltlm2)
1729 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1730 END DO
1731!
1732 DO ng=1,ngrids
1733 DO tile=first_tile(ng),last_tile(ng),+1
1734 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1735 END DO
1736 END DO
1737!
1738! Write the final TL increment to Rec2 of the ITL file.
1739!
1740 DO ng=1,ngrids
1741# ifdef DISTRIBUTE
1742 CALL tl_wrt_ini (ng, myrank, ltlm2, rec2)
1743# else
1744 CALL tl_wrt_ini (ng, -1, ltlm2, rec2)
1745# endif
1746 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1747 END DO
1748!
1749! Now update the non-linear model initial conditions.
1750!
1751 linl=my_outer+1
1752 DO ng=1,ngrids
1753 CALL get_state (ng, inlm, 9, ini(ng), linl, lnew(ng))
1754 IF (timeiau(ng).eq.0.0_dp) THEN
1755 CALL get_state (ng, iadm, 4, itl(ng), rec2, ltlm1)
1756 END IF
1757 END DO
1758 DO ng=1,ngrids
1759 IF (timeiau(ng).eq.0.0_dp) THEN
1760 DO tile=first_tile(ng),last_tile(ng),+1
1761 CALL ini_adjust (ng, tile, ltlm1, lnew(ng))
1762 END DO
1763 END IF
1764 END DO
1765!
1766! Write out nonlinear model initial conditions into INIname, record
1767! INI(ng)%Rindex.
1768!
1769 DO ng=1,ngrids
1770# ifdef DISTRIBUTE
1771 CALL wrt_ini (ng, myrank, lnew(ng))
1772# else
1773 CALL wrt_ini (ng, -1, lnew(ng))
1774# endif
1775 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1776
1777# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1778 defined adjust_boundary
1779# ifdef DISTRIBUTE
1780 CALL wrt_frc_ad (ng, myrank, ltlm1, ini(ng)%Rindex)
1781# else
1782 CALL wrt_frc_ad (ng, -1, ltlm1, ini(ng)%Rindex)
1783# endif
1784 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1785# endif
1786 END DO
1787!
1788! Compute the new B^-1(x(k)-x(k-1)) term.
1789! Gather the final adjoint solutions increments, sum and
1790! save in record 4 of the ITL file. Use the tl arrays as temporary
1791! storage.
1792!
1793! First add the augmented piece which is computed from the previous
1794! sum.
1795!
1796 DO ng=1,ngrids
1797 CALL get_state (ng, itlm, 4, itl(ng), rec4, ltlm1)
1798 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1799 END DO
1800!
1801 DO ng=1,ngrids
1802 DO tile=first_tile(ng),last_tile(ng),+1
1803 CALL aug_oper (ng, tile, ltlm1, ltlm2)
1804 END DO
1805 END DO
1806!
1807 DO ng=1,ngrids
1808 DO tile=first_tile(ng),last_tile(ng),+1
1809 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1810 END DO
1811 END DO
1812!
1813 DO ng=1,ngrids
1814 adrec=nrec(ng)
1815 CALL get_state (ng, itlm, 4, adm(ng), adrec, ltlm1)
1816 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1817 END DO
1818!
1819 DO ng=1,ngrids
1820 DO tile=first_tile(ng),last_tile(ng),+1
1821 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1822 END DO
1823 END DO
1824!
1825! Write the current sum of adjoint solutions into record 4 of the ITL
1826! file.
1827!
1828 DO ng=1,ngrids
1829# ifdef DISTRIBUTE
1830 CALL tl_wrt_ini (ng, myrank, ltlm2, rec4)
1831# else
1832 CALL tl_wrt_ini (ng, -1, ltlm2, rec4)
1833# endif
1834 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1835 END DO
1836# endif
1837!
1838! Convert the current adjoint solution in ADM(ng)%name to impulse
1839! forcing. Write out impulse forcing into TLF(ng)%name NetCDF file.
1840! To facilitate the forcing to the TLM and RPM, the forcing is
1841! processed and written in increasing time coordinates (recall that
1842! the adjoint solution in ADM(ng)%name is backwards in time).
1843!
1844 IF (master) THEN
1845 WRITE (stdout,30) my_outer, inner
1846 END IF
1847 DO ng=1,ngrids
1848 tlf(ng)%Rindex=0
1849# ifdef DISTRIBUTE
1850 CALL wrt_impulse (ng, myrank, iadm, adm(ng)%name)
1851# else
1852 CALL wrt_impulse (ng, -1, iadm, adm(ng)%name)
1853# endif
1854 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1855 END DO
1856
1857# ifdef RPCG
1858!
1859! Now Compute the augmented corrections to the weak constraint
1860! forcing terms. The sums so far are in records 6 to
1861! 6+nADrec/2.
1862!
1863 DO ng=1,ngrids
1864 DO i=1,nadrec(ng)/2
1865 irec=i
1866 jrec=rec5+i
1867 CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
1868 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1869!
1870 DO tile=first_tile(ng),last_tile(ng),+1
1871 CALL aug_oper (ng, tile, ltlm1, ltlm2)
1872 END DO
1873!
1874! Complete the computation of the TLF forcing term by adding the
1875! contribution from the augmented adjoint propagator. Specify
1876! the value of 7 for the model variable since this the special
1877! case in get_state for reading the impulse forcing.
1878!
1879!! CALL get_state (ng, iTLM, 4, TLF(ng), irec, LTLM1) ! TEST
1880 CALL get_state (ng, 7, 4, tlf(ng), irec, ltlm1)
1881 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1882!
1883 DO tile=first_tile(ng),last_tile(ng),+1
1884!! CALL sum_grad (ng, tile, LTLM1, LTLM2) ! TEST
1885 CALL sum_imp (ng, tile, ltlm2)
1886 END DO
1887!
1888! Write the final forcing increment to he TLF file.
1889! Note the original TLF file is overwritten at this point.
1890!
1891 tlf(ng)%Rindex=0
1892# ifdef DISTRIBUTE
1893 CALL wrt_aug_imp (ng, myrank, itlm, ltlm2, i)
1894# else
1895 CALL wrt_aug_imp (ng, -1, itlm, ltlm2, i)
1896# endif
1897
1898 END DO
1899 END DO
1900!
1901! Gather the increments from the final inner-loop and
1902! save in record 3 of the ITL file. The current increment
1903! is in record 2 and the sum so far is in record 3.
1904!
1905 DO ng=1,ngrids
1906 CALL get_state (ng, itlm, 4, itl(ng), rec2, ltlm1)
1907 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1908 CALL get_state (ng, itlm, 4, itl(ng), rec3, ltlm2)
1909 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1910 END DO
1911!
1912 DO ng=1,ngrids
1913 DO tile=first_tile(ng),last_tile(ng),+1
1914 CALL sum_grad (ng, tile, ltlm1, ltlm2)
1915 END DO
1916 END DO
1917!
1918! Write the current sum into record 3 of the ITL file.
1919!
1920 DO ng=1,ngrids
1921# ifdef DISTRIBUTE
1922 CALL tl_wrt_ini (ng, myrank, ltlm2, rec3)
1923# else
1924 CALL tl_wrt_ini (ng, -1, ltlm2, rec3)
1925# endif
1926 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1927 END DO
1928
1929# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS || \
1930 defined adjust_boundary
1931!
1932! Write the current sum of the forcing and boundary increments
1933! into the INI file into index INI(ng)%Rindex-1 (i.e. we are
1934! overwriting the previous fields. Read the current sum into the
1935! adjoint arrays since this is what wrt_frc_AD uses.
1936!
1937 DO ng=1,ngrids
1938 CALL get_state (ng, iadm, 4, itl(ng), rec3, ltlm1)
1939 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1940 END DO
1941 DO ng=1,ngrids
1942# ifdef DISTRIBUTE
1943 CALL wrt_frc_ad (ng, myrank, ltlm1, ini(ng)%Rindex)
1944# else
1945 CALL wrt_frc_ad (ng, -1, ltlm1, ini(ng)%Rindex)
1946# endif
1947 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1948 END DO
1949# endif
1950!
1951! Gather the model error forcing increments and update the
1952! sum in records 6 to 6+nADrec/2.
1953!
1954 DO ng=1,ngrids
1955 DO i=1,nadrec(ng)/2
1956 irec=i
1957 jrec=rec5+i
1958 CALL get_state (ng, itlm, 4, itl(ng), jrec, ltlm1)
1959 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1960!! CALL get_state (ng, iTLM, 4, TLF(ng)%, irec, LTLM2) ! TEST
1961 CALL get_state (ng, 7, 4, tlf(ng), irec, ltlm2)
1962 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1963!
1964 DO tile=first_tile(ng),last_tile(ng),+1
1965!! CALL sum_grad (ng, tile, LTLM1, LTLM2) ! TEST
1966 CALL sum_imp (ng, tile, ltlm1)
1967 END DO
1968!
1969! Write the current sum into record jrec of the ITL file.
1970!
1971# ifdef DISTRIBUTE
1972!! CALL tl_wrt_ini (ng, MyRank, LTLM2, jrec) ! TEST
1973 CALL tl_wrt_ini (ng, myrank, ltlm1, jrec)
1974# else
1975!! CALL tl_wrt_ini (ng, -1, LTLM2, jrec) ! TEST
1976 CALL tl_wrt_ini (ng, -1, ltlm1, jrec)
1977# endif
1978 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1979 END DO
1980 END DO
1981!
1982! Now compute the background cost function Jb0.
1983! First compute the contribution from the increments in the
1984! initial conditions, surface forcing, and boundary conditions.
1985!
1986 jb0(my_outer)=0.0_r8
1987 DO ng=1,ngrids
1988 CALL get_state (ng, iadm, 8, itl(ng), rec4, ltlm1)
1989 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1990 CALL get_state (ng, itlm, 8, itl(ng), rec3, ltlm1)
1991 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1992 END DO
1993!
1994 DO ng=1,ngrids
1995 DO tile=first_tile(ng),last_tile(ng),+1
1996 CALL comp_jb0 (ng, tile, itlm, my_outer, ltlm1, ltlm1)
1997 END DO
1998 END DO
1999!
2000! Now compute the contribution of model error terms to the
2001! background cost function Jb0.
2002!
2003! NOTE: I THINK WE NEED A BETTER WAY OF COMPUTING THE OBS ERROR
2004! CONTRIBUTION TO Jb. THE FOLLOWING CALCULATION IS DANGEROUS BECAUSE
2005! IT RELIES ON THE FACT THAT THE SURFACING FORCING AND OBC TL ARRAYS
2006! READ FROM jrec2 ARE ZERO. THIS MEANS THAT ONLY THE MODEL STATE
2007! CONTRIBUTES TO Jb FOR THE MODEL ERROR TERM, WHICH OF COURSE IS
2008! WHAT SHOULD BE THE CASE. THE SURFACE FORCING AND OBC CONTRIBUTIONS
2009! ARE COMPUTED IN THE PREVIOUS CALL TO COMP_JB0.
2010!
2011 DO ng=1,ngrids
2012 DO irec=1,nadrec(ng)/2
2013 jrec1=rec5+irec
2014 jrec2=rec5+nadrec(ng)/2+irec
2015 CALL get_state (ng, iadm, 8, itl(ng), jrec1, ltlm1)
2016 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2017 CALL get_state (ng, itlm, 8, itl(ng), jrec2, ltlm1)
2018 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2019!
2020 DO tile=first_tile(ng),last_tile(ng),+1
2021 CALL comp_jb0 (ng, tile, itlm, my_outer, ltlm1, ltlm1)
2022 END DO
2023 END DO
2024 END DO
2025!
2026! Overwrite the TFL netcdf file with the sum of the model error
2027! forcing increments - required for the next run of the NLM and
2028! TLM.
2029!
2030 DO ng=1,ngrids
2031 DO irec=1,nadrec(ng)/2
2032 jrec=rec5+irec
2033 CALL get_state (ng, itlm, 8, itl(ng), jrec, ltlm1)
2034!
2035 tlf(ng)%Rindex=0
2036# ifdef DISTRIBUTE
2037 CALL wrt_aug_imp (ng, myrank, itlm, ltlm1, irec)
2038# else
2039 CALL wrt_aug_imp (ng, -1, itlm, ltlm1, irec)
2040# endif
2041 END DO
2042 END DO
2043# endif
2044
2045# ifdef PROFILE
2046!
2047! Stop profile clock
2048!
2049 DO ng=1,ngrids
2050 DO thread=thread_range
2051 CALL wclock_off (ng, itlm, 87, __line__, myfile)
2052 END DO
2053 END DO
2054# endif
2055!
2056 10 FORMAT (a,'_outer',i0,'.nc')
2057 20 FORMAT (/,1x,a,1x,'ROMS: started time-stepping:', &
2058 & ' (Grid: ',i0,', Outer: ',i2.2,', Inner: ',i3.3, &
2059 ', TimeSteps: ',i0,' - ',i0,')',/)
2060 30 FORMAT (/,' Converting Convolved Adjoint Trajectory to', &
2061 & ' Impulses: Outer = ',i3.3,' Inner = ',i3.3,/)
2062!
2063 RETURN
subroutine ad_initial(ng)
Definition ad_initial.F:4
subroutine ad_main2d
Definition ad_main2d.F:586
subroutine ad_main3d(runinterval)
Definition ad_main3d.F:4
subroutine, public comp_jb0(ng, tile, model, outloop, ltlm, ladj)
Definition comp_Jb0.F:49
subroutine, public aug_oper(ng, tile, linp, lout)
Definition comp_Jb0.F:359
subroutine, public load_tltoad(ng, tile, linp, lout, add)
subroutine, public ini_adjust(ng, tile, linp, lout)
Definition ini_adjust.F:69
subroutine, public sum_grad(ng, tile, linp, lout)
Definition sum_grad.F:26
subroutine, public sum_imp(ng, tile, lout)
Definition sum_imp.F:21
subroutine tl_initial(ng)
Definition tl_initial.F:4
subroutine tl_main2d
Definition tl_main2d.F:429
subroutine tl_main3d(runinterval)
Definition tl_main3d.F:4

References ad_def_his_mod::ad_def_his(), ad_initial(), ad_main2d(), ad_main3d(), ad_wrt_his_mod::ad_wrt_his(), mod_iounits::adm, comp_jb0_mod::aug_oper(), mod_scalars::balance, zeta_balance_mod::balance_ref(), zeta_balance_mod::biconj(), congrad_mod::cg_read_congrad(), rpcg_lanczos_mod::cg_read_rpcg(), comp_jb0_mod::comp_jb0(), congrad_mod::congrad(), mod_iounits::dav, def_error_mod::def_error(), def_hessian_mod::def_hessian(), def_impulse_mod::def_impulse(), def_ini_mod::def_ini(), def_mod_mod::def_mod(), mod_scalars::dt, edit_multifile(), convolve_mod::error_covariance(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), mod_scalars::frcrec, mod_scalars::frequentimpulse, mod_iounits::fwd, get_state_mod::get_state(), mod_fourdvar::havenlmod, mod_iounits::his, mod_param::iadm, mod_ncparam::idnlmi, mod_ncparam::idobss, mod_ncparam::idoerr, mod_ncparam::idoval, mod_ncparam::idtime, mod_iounits::ini, ini_adjust_mod::ini_adjust(), mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_ocean::initialize_ocean(), mod_scalars::initime, mod_param::inlm, mod_scalars::inner, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_param::irpm, mod_ncparam::isfsur, mod_iounits::itl, mod_param::itlm, mod_fourdvar::jb0, mod_parallel::last_tile, mod_fourdvar::laugweak, lbck, mod_scalars::ldefadj, mod_scalars::ldeferr, mod_scalars::ldefhss, mod_scalars::ldefini, mod_scalars::ldefitl, mod_scalars::ldefmod, mod_scalars::ldeftlf, mod_fourdvar::lhotstart, lini, mod_stepping::lnew, ini_adjust_mod::load_tltoad(), mod_stepping::lold, mod_scalars::lreadblk, mod_scalars::lreadfrc, mod_scalars::lreadfwd, mod_scalars::lreadqck, ltlm1, ltlm2, mod_scalars::lwrtstate2d, mod_parallel::master, mod_parallel::myrank, mod_scalars::nadj, mod_scalars::nbrec, mod_fourdvar::ndatum, mod_scalars::nfrec, mod_param::ngrids, mod_scalars::ninner, mod_fourdvar::nlmodval, mod_scalars::nobc, mod_scalars::noerror, mod_scalars::nrun, mod_scalars::nsff, mod_scalars::ntend, mod_scalars::ntimes, mod_scalars::ntstart, mod_scalars::obc_time, mod_iounits::obs, mod_fourdvar::obserr, mod_fourdvar::obsscale, mod_fourdvar::obsval, mod_iounits::qck, rec1, rec2, rec3, rec4, rec5, rpcg_lanczos_mod::rpcg_lanczos(), mod_scalars::sec2day, mod_scalars::sf_time, mod_iounits::sourcefile, mod_iounits::stdout, sum_grad_mod::sum_grad(), sum_imp_mod::sum_imp(), mod_scalars::tdays, mod_scalars::timeiau, tl_def_ini_mod::tl_def_ini(), tl_initial(), tl_main2d(), tl_main3d(), tl_wrt_ini_mod::tl_wrt_ini(), mod_iounits::tlf, mod_ncparam::vname, wclock_off(), wclock_on(), wrt_aug_imp_mod::wrt_aug_imp(), wrt_ini_mod::wrt_frc_ad(), wrt_hessian_mod::wrt_hessian(), wrt_impulse_mod::wrt_impulse(), wrt_ini_mod::wrt_ini(), mod_fourdvar::wrtforce, mod_fourdvar::wrtmisfit, mod_fourdvar::wrtnlmod, mod_fourdvar::wrttlmod, and mod_fourdvar::wrtzetaref.

Here is the call graph for this function:

◆ posterior_error()

subroutine, public rbl4dvar_mod::posterior_error ( real(dp), intent(in) runinterval)

Definition at line 2806 of file rbl4dvar.F.

2807!
2808!=======================================================================
2809! !
2810! This routine computes posterior analysis error covariance matrix !
2811! using the Lanczos vectors. !
2812! !
2813! On Input: !
2814! !
2815! RunInterval Kernel time stepping window (seconds) !
2816! !
2817! Warning: !
2818! !
2819! Currently, this code only works for a single outer-loop. !
2820! !
2821!=======================================================================
2822!
2823 USE random_ic_mod, ONLY : random_ic
2824!
2825! Imported variable declarations
2826!
2827 real(dp), intent(in) :: RunInterval
2828!
2829! Local variable declarations.
2830!
2831 logical :: Ltrace
2832!
2833 integer :: my_inner, my_outer, ng, tile
2834 integer :: Fcount, Rec
2835!
2836 character (len=8) :: driver = 'rbl4dvar'
2837
2838 character (len=*), parameter :: MyFile = &
2839 & __FILE__//", posterior_error"
2840!
2841 sourcefile=myfile
2842
2843# if defined POSTERIOR_ERROR_I || defined POSTERIOR_ERROR_F
2844!
2845!-----------------------------------------------------------------------
2846! Compute full (diagonal) posterior analysis error covariance matrix.
2847!-----------------------------------------------------------------------
2848!
2849! Clear tangent and adjoint arrays because they are used as
2850! work arrays below.
2851!
2852 DO ng=1,ngrids
2853 DO tile=first_tile(ng),last_tile(ng),+1
2854 CALL initialize_ocean (ng, tile, iadm)
2855 CALL initialize_ocean (ng, tile, itlm)
2856 CALL initialize_forces (ng, tile, iadm)
2857 CALL initialize_forces (ng, tile, itlm)
2858# ifdef ADJUST_BOUNDARY
2859 CALL initialize_boundary (ng, tile, iadm)
2860 CALL initialize_boundary (ng, tile, itlm)
2861# endif
2862 END DO
2863 END DO
2864!
2865! Compute the diagonal of the posterior/analysis error covariance
2866! matrix. The result is written to record 2 of the ITL netcdf file.
2867!
2868 var_oloop : DO my_outer=1,nouter
2869 outer=my_outer
2870 DO ng=1,ngrids
2871 DO tile=first_tile(ng),last_tile(ng),+1
2872 CALL posterior_var (ng, tile, itlm, outer)
2873 END DO
2874 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2875 END DO
2876 END DO var_oloop
2877!
2878! Write out the diagonal of the posterior/analysis covariance matrix
2879! which is in tl_var(Rec1) to 4DVar error NetCDF file.
2880!
2881 DO ng=1,ngrids
2882# ifdef DISTRIBUTE
2883 CALL wrt_error (ng, myrank, rec1, rec1)
2884# else
2885 CALL wrt_error (ng, -1, rec1, rec1)
2886# endif
2887 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2888 END DO
2889!
2890! Clear tangent and adjoint arrays because they are used as
2891! work arrays below.
2892!
2893 DO ng=1,ngrids
2894 DO tile=first_tile(ng),last_tile(ng),+1
2895 CALL initialize_ocean (ng, tile, iadm)
2896 CALL initialize_ocean (ng, tile, itlm)
2897 CALL initialize_forces (ng, tile, iadm)
2898 CALL initialize_forces (ng, tile, itlm)
2899# ifdef ADJUST_BOUNDARY
2900 CALL initialize_boundary (ng, tile, iadm)
2901 CALL initialize_boundary (ng, tile, itlm)
2902# endif
2903 END DO
2904 END DO
2905# endif
2906
2907# ifdef POSTERIOR_EOFS
2908!
2909!-----------------------------------------------------------------------
2910! Compute the posterior analysis error covariance matrix EOFs using a
2911! Lanczos algorithm.
2912!
2913! NOTE: Currently, this code only works for a single outer-loop.
2914!-----------------------------------------------------------------------
2915!
2916 IF (master) WRITE (stdout,10)
2917!
2918! Estimate first the trace of the posterior analysis error
2919! covariance matrix since the evolved and convolved Lanczos
2920! vectors stored in the Hessian NetCDF file will be destroyed
2921! later.
2922!
2923 ltrace=.true.
2924
2925 trace_oloop : DO my_outer=1,nouter
2926 outer=my_outer
2927 inner=0
2928
2929 trace_iloop : DO my_inner=1,nposti
2930 inner=my_inner
2931!
2932! Initialize the tangent linear variables with a random vector
2933! comprised of +1 and -1 elements randomly chosen.
2934!
2935 DO ng=1,ngrids
2936 DO tile=first_tile(ng),last_tile(ng),+1
2937 CALL random_ic (ng, tile, itlm, inner, outer, &
2938 & lold(ng), ltrace)
2939 END DO
2940 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2941 END DO
2942!
2943! Apply horizontal convolution.
2944!
2945 CALL convolve (driver, lini, lold, lnew)
2946 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2947!
2948! Compute Lanczos vector and eigenvectors of the posterior analysis
2949! error covariance matrix.
2950!
2951 DO ng=1,ngrids
2952 DO tile=first_tile(ng),last_tile(ng),+1
2953 CALL posterior (ng, tile, itlm, inner, outer, ltrace)
2954 END DO
2955 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2956 END DO
2957
2958 END DO trace_iloop
2959
2960 END DO trace_oloop
2961!
2962! Estimate posterior analysis error covariance matrix.
2963!
2964 ltrace=.false.
2965
2966 post_oloop : DO my_outer=1,nouter
2967 outer=my_outer
2968 inner=0
2969!
2970! The Lanczos algorithm requires to save all the Lanczos vectors.
2971! They are used to compute the posterior EOFs.
2972!
2973 DO ng=1,ngrids
2974 adm(ng)%Rindex=0
2975 fcount=adm(ng)%load
2976 adm(ng)%Nrec(fcount)=0
2977 END DO
2978
2979 post_iloop : DO my_inner=0,nposti
2980 inner=my_inner
2981!
2982! Read first record of ITL file and apply convolutions.
2983!
2984! NOTE: If inner=0, we would like to use a random starting vector.
2985! For now we can use what ever is in record 1.
2986!
2987 IF (inner.ne.0) THEN
2988 DO ng=1,ngrids
2989 rec=1
2990 CALL get_state (ng, itlm, 1, itl(ng), rec, lold(ng))
2991 IF (founderror(exit_flag, noerror, &
2992 & __line__, myfile)) RETURN
2993 END DO
2994 ELSE
2995
2996 DO ng=1,ngrids
2997 DO tile=first_tile(ng),last_tile(ng),+1
2998 CALL random_ic (ng, tile, itlm, inner, outer, &
2999 & lold(ng), ltrace)
3000 END DO
3001 IF (founderror(exit_flag, noerror, &
3002 & __line__, myfile)) RETURN
3003 END DO
3004 END IF
3005!
3006! Apply horizontal convolution.
3007!
3008 CALL convolve (driver, lini, lold, lnew)
3009 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3010!
3011! Compute Lanczos vector and eigenvectors of the posterior analysis
3012! error covariance matrix.
3013!
3014 DO ng=1,ngrids
3015 DO tile=first_tile(ng),last_tile(ng),+1
3016 CALL posterior (ng, tile, itlm, inner, outer, ltrace)
3017 END DO
3018 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3019 END DO
3020!
3021! Write the Lanczos vectors of the posterior error covariance
3022! to the adjoint NetCDF file.
3023!
3024 DO ng=1,ngrids
3025# if defined ADJUST_STFLUX || defined ADJUST_WSTRESS
3026 lfout(ng)=lnew(ng)
3027# endif
3028# ifdef ADJUST_BOUNDARY
3029 lbout(ng)=lnew(ng)
3030# endif
3031 kstp(ng)=lnew(ng)
3032# ifdef SOLVE3D
3033 nstp(ng)=lnew(ng)
3034# endif
3035 lwrtstate2d(ng)=.true.
3036# ifdef DISTRIBUTE
3037 CALL ad_wrt_his (ng, myrank)
3038# else
3039 CALL ad_wrt_his (ng, -1)
3040# endif
3041 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3042 lwrtstate2d(ng)=.false.
3043 END DO
3044!
3045! Write out tangent linear model initial conditions and tangent
3046! linear surface forcing adjustments for next inner
3047! loop into ITL(ng)%name (record Rec1).
3048!
3049 DO ng=1,ngrids
3050# ifdef DISTRIBUTE
3051 CALL tl_wrt_ini (ng, myrank, lnew(ng), rec1)
3052# else
3053 CALL tl_wrt_ini (ng, -1, lnew(ng), rec1)
3054# endif
3055 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
3056 END DO
3057
3058 END DO post_iloop
3059
3060 END DO post_oloop
3061
3062# endif
3063!
3064! Done. Set history file ID to closed state since we manipulated
3065! its indices with the forward file ID which was closed above.
3066!
3067 DO ng=1,ngrids
3068 IF (his(ng)%IOtype.eq.io_nf90) THEN
3069 his(ng)%ncid=-1
3070# if defined PIO_LIB && defined DISTRIBUTE
3071 ELSE IF (his(ng)%IOtype.eq.io_pio)
3072 his(ng)%pioFile%fh=-1
3073# endif
3074 END IF
3075 END DO
3076
3077# ifdef POSTERIOR_EOFS
3078!
3079 10 FORMAT (/,' <<<< Posterior Analysis Error Covariance Matrix', &
3080 & ' Estimation >>>>',/)
3081# endif
3082!
3083 RETURN
subroutine, public random_ic(ng, tile, model, innloop, outloop, lout, ltrace)
Definition random_ic.F:32

References ad_wrt_his_mod::ad_wrt_his(), mod_iounits::adm, convolve_mod::convolve(), mod_scalars::exit_flag, mod_parallel::first_tile, get_state_mod::get_state(), mod_iounits::his, mod_param::iadm, mod_boundary::initialize_boundary(), mod_forces::initialize_forces(), mod_ocean::initialize_ocean(), mod_scalars::inner, mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_iounits::itl, mod_param::itlm, mod_stepping::kstp, mod_parallel::last_tile, mod_stepping::lbout, mod_stepping::lfout, lini, mod_stepping::lnew, mod_stepping::lold, mod_scalars::lwrtstate2d, mod_parallel::master, mod_parallel::myrank, mod_param::ngrids, mod_scalars::noerror, mod_scalars::nouter, mod_fourdvar::nposti, mod_stepping::nstp, mod_scalars::outer, posterior_mod::posterior(), posterior_var_mod::posterior_var(), random_ic_mod::random_ic(), rec1, mod_iounits::sourcefile, mod_iounits::stdout, tl_wrt_ini_mod::tl_wrt_ini(), and wrt_error_mod::wrt_error().

Here is the call graph for this function:

◆ prior_error()

subroutine, public rbl4dvar_mod::prior_error ( integer, intent(in) ng)

Definition at line 2635 of file rbl4dvar.F.

2636!
2637!=======================================================================
2638! !
2639! This routine processes background prior error covariance standard !
2640! deviations and normalization coefficients. !
2641! !
2642! On Input: !
2643! !
2644! ng Nested grid number !
2645! !
2646!=======================================================================
2647!
2648 USE def_norm_mod, ONLY : def_norm
2650 USE strings_mod, ONLY : founderror
2651!
2652! Imported variable declarations
2653!
2654 integer, intent(in) :: ng
2655!
2656! Local variable declarations.
2657!
2658 integer :: tile
2659 integer :: NRMrec, STDrec, Tindex
2660!
2661 character (len=*), parameter :: MyFile = &
2662 & __FILE__//", prior_error"
2663!
2664 sourcefile=myfile
2665!
2666!-----------------------------------------------------------------------
2667! Set application grid, metrics, and associated variables and
2668! parameters.
2669!-----------------------------------------------------------------------
2670!
2671! The ROMS application grid configuration is done once. It is usually
2672! done in the "initial" kernel routine. However, since we are calling
2673! the "normalization" routine here, we need several grid variables and
2674! parameter. Also, if reading only water points, we need to know the
2675! land/sea mask arrays to unpack.
2676!
2677 IF (setgridconfig(ng)) THEN
2678 CALL set_grid (ng, inlm)
2679 END IF
2680!
2681!-----------------------------------------------------------------------
2682! Read in standard deviation factors for error covariance.
2683!-----------------------------------------------------------------------
2684!
2685! Initial conditions standard deviation. They are loaded in Tindex=1
2686! of the e_var(...,Tindex) state variables.
2687!
2688 IF (lreadstd(ng)) THEN
2689 stdrec=1
2690 tindex=1
2691# ifdef STD_MODEL
2692 CALL get_state (ng, 10, 10, std(5,ng), stdrec, tindex)
2693# else
2694 CALL get_state (ng, 10, 10, std(1,ng), stdrec, tindex)
2695# endif
2696 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2697 END IF
2698!
2699! Model error standard deviation. They are loaded in Tindex=2
2700! of the e_var(...,Tindex) state variables.
2701!
2702 stdrec=1
2703 tindex=2
2704 IF (nsa.eq.2) THEN
2705 CALL get_state (ng, 11, 11, std(2,ng), stdrec, tindex)
2706 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2707 END IF
2708
2709# ifdef ADJUST_BOUNDARY
2710!
2711! Open boundary conditions standard deviation.
2712!
2713 stdrec=1
2714 tindex=1
2715 CALL get_state (ng, 12, 12, std(3,ng), stdrec, tindex)
2716 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2717# endif
2718
2719# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2720!
2721! Surface forcing standard deviation.
2722!
2723 stdrec=1
2724 tindex=1
2725 CALL get_state (ng, 13, 13, std(4,ng), stdrec, tindex)
2726 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2727# endif
2728
2729# ifdef STD_MODEL
2730!
2731! If computing the standard deviation (STD) directly from the
2732! background (prior) field as an alternative to climatological
2733! read from the input NetCDF file, define the output STD NetCDF
2734! file.
2735!
2736 IF (ldefstd(ng)) THEN
2737 CALL def_std (ng)
2738 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2739 END IF
2740# endif
2741!
2742!-----------------------------------------------------------------------
2743! Error covariance normalization coefficients.
2744!-----------------------------------------------------------------------
2745!
2746! Compute or read in the error covariance normalization factors.
2747! If computing, write out factors to NetCDF. This is an expensive
2748! computation that needs to be computed only once for a particular
2749! application grid and decorrelation scales.
2750!
2751 IF (any(lwrtnrm(:,ng))) THEN
2752 CALL def_norm (ng, inlm, 1)
2753 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2754
2755 IF (nsa.eq.2) THEN
2756 CALL def_norm (ng, inlm, 2)
2757 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2758 END IF
2759
2760# ifdef ADJUST_BOUNDARY
2761 CALL def_norm (ng, inlm, 3)
2762 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2763# endif
2764
2765# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2766 CALL def_norm (ng, inlm, 4)
2767 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2768# endif
2769!
2770 DO tile=first_tile(ng),last_tile(ng),+1
2771 CALL normalization (ng, tile, 2)
2772 END DO
2773 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2774 ldefnrm(1:4,ng)=.false.
2775 lwrtnrm(1:4,ng)=.false.
2776
2777 ELSE
2778
2779 nrmrec=1
2780 CALL get_state (ng, 14, 14, nrm(1,ng), nrmrec, 1)
2781 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2782
2783 IF (nsa.eq.2) THEN
2784 CALL get_state (ng, 15, 15, nrm(2,ng), nrmrec, 2)
2785 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2786 END IF
2787
2788# ifdef ADJUST_BOUNDARY
2789 CALL get_state (ng, 16, 16, nrm(3,ng), nrmrec, 1)
2790 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2791# endif
2792
2793# if defined ADJUST_WSTRESS || defined ADJUST_STFLUX
2794 CALL get_state (ng, 17, 17, nrm(4,ng), nrmrec, 1)
2795 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2796# endif
2797
2798 END IF
2799!
2800 RETURN
subroutine, public def_norm(ng, model, ifile)
Definition def_norm.F:53
subroutine, public normalization(ng, tile, ifac)
logical function, public founderror(flag, noerr, line, routine)
Definition strings.F:52
subroutine set_grid(ng, model)
Definition set_grid.F:3

References def_norm_mod::def_norm(), mod_scalars::exit_flag, mod_parallel::first_tile, strings_mod::founderror(), get_state_mod::get_state(), mod_param::inlm, mod_parallel::last_tile, mod_scalars::ldefnrm, mod_scalars::lreadstd, mod_scalars::lwrtnrm, mod_scalars::noerror, normalization_mod::normalization(), mod_iounits::nrm, mod_param::nsa, posterior_mod::posterior_eofs(), set_grid(), mod_scalars::setgridconfig, mod_iounits::sourcefile, and mod_iounits::std.

Here is the call graph for this function:

Variable Documentation

◆ lbck

integer rbl4dvar_mod::lbck = 2

Definition at line 199 of file rbl4dvar.F.

199 integer :: Lbck = 2

Referenced by background_initialize(), and increment().

◆ ldone

logical rbl4dvar_mod::ldone

Definition at line 196 of file rbl4dvar.F.

196 logical :: Ldone

◆ lini

integer rbl4dvar_mod::lini = 1

Definition at line 200 of file rbl4dvar.F.

200 integer :: Lini = 1

Referenced by increment(), and posterior_error().

◆ ltlm1

integer rbl4dvar_mod::ltlm1 = 1

Definition at line 201 of file rbl4dvar.F.

201 integer :: LTLM1 = 1

Referenced by analysis_initialize(), and increment().

◆ ltlm2

integer rbl4dvar_mod::ltlm2 = 2

Definition at line 202 of file rbl4dvar.F.

202 integer :: LTLM2 = 2

Referenced by increment().

◆ rec1

integer rbl4dvar_mod::rec1 = 1

Definition at line 203 of file rbl4dvar.F.

203 integer :: Rec1 = 1

Referenced by increment(), and posterior_error().

◆ rec2

integer rbl4dvar_mod::rec2 = 2

Definition at line 204 of file rbl4dvar.F.

204 integer :: Rec2 = 2

Referenced by increment().

◆ rec3

integer rbl4dvar_mod::rec3 = 3

Definition at line 205 of file rbl4dvar.F.

205 integer :: Rec3 = 3

Referenced by analysis_initialize(), and increment().

◆ rec4

integer rbl4dvar_mod::rec4 = 4

Definition at line 206 of file rbl4dvar.F.

206 integer :: Rec4 = 4

Referenced by increment().

◆ rec5

integer rbl4dvar_mod::rec5 = 5

Definition at line 207 of file rbl4dvar.F.

207 integer :: Rec5 = 5

Referenced by increment().