ROMS
Loading...
Searching...
No Matches
Functions/Subroutines
rpcg_lanczos_mod Module Reference

Functions/Subroutines

subroutine, public rpcg_lanczos (ng, model, outloop, innloop, ninnloop, lcgini)
 
subroutine, private rpevecs (ng, outloop, ninnloop)
 
subroutine, private rprecond (ng, lscale, laug, outloop, ninnloop, py)
 
subroutine, private cg_write_rpcg (ng, model, innloop, outloop, jf, jdata, jmod, jopt, jb, jobs, jact, preducv, preducy)
 
subroutine, private cg_write_rpcg_nf90 (ng, model, innloop, outloop, jf, jdata, jmod, jopt, jb, jobs, jact, preducv, preducy)
 
subroutine, private cg_write_rpcg_pio (ng, model, innloop, outloop, jf, jdata, jmod, jopt, jb, jobs, jact, preducv, preducy)
 
subroutine, public cg_read_rpcg (ng, model, outloop)
 
subroutine, private cg_read_rpcg_nf90 (ng, model, outloop)
 
subroutine, private cg_read_rpcg_pio (ng, model, outloop)
 

Function/Subroutine Documentation

◆ cg_read_rpcg()

subroutine, public rpcg_lanczos_mod::cg_read_rpcg ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) outloop )

Definition at line 2034 of file rpcg_lanczos.F.

2035!***********************************************************************
2036! !
2037! This routine reads conjugate gradient variables for previous outer !
2038! loops using either the standard NetCDF library or the Parallel-IO !
2039! (PIO) library. !
2040! !
2041!=======================================================================
2042!
2043! Imported variable declarations
2044!
2045 integer, intent(in) :: ng, model, outLoop
2046!
2047! Local variable declarations.
2048!
2049 character (len=*), parameter :: MyFile = &
2050 & __FILE__//", cg_read_rpcg"
2051!
2052 sourcefile=myfile
2053!
2054!-----------------------------------------------------------------------
2055! Read in conjugate gradient variables.
2056!-----------------------------------------------------------------------
2057!
2058 SELECT CASE (dav(ng)%IOtype)
2059 CASE (io_nf90)
2060 CALL cg_read_rpcg_nf90 (ng, model, outloop)
2061
2062# if defined PIO_LIB && defined DISTRIBUTE
2063 CASE (io_pio)
2064 CALL cg_read_rpcg_pio (ng, model, outloop)
2065# endif
2066 CASE DEFAULT
2067 IF (master) WRITE (stdout,10) dav(ng)%IOtype
2068 exit_flag=3
2069 END SELECT
2070 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2071!
2072 10 FORMAT (' CG_READ_RPCG - Illegal output type, io_type = ',i0)
2073!
2074 RETURN

References cg_read_rpcg_nf90(), cg_read_rpcg_pio(), mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_parallel::master, mod_scalars::noerror, mod_iounits::sourcefile, and mod_iounits::stdout.

Referenced by rbl4dvar_mod::increment().

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

◆ cg_read_rpcg_nf90()

subroutine, private rpcg_lanczos_mod::cg_read_rpcg_nf90 ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) outloop )
private

Definition at line 2078 of file rpcg_lanczos.F.

2079!***********************************************************************
2080!
2081 USE mod_netcdf
2082!
2083! Imported variable declarations
2084!
2085 integer, intent(in) :: ng, model, outLoop
2086!
2087! Local variable declarations.
2088!
2089 integer :: status
2090!
2091 character (len=*), parameter :: MyFile = &
2092 & __FILE__//", cg_read_rpcg_nf90"
2093!
2094 sourcefile=myfile
2095!
2096!-----------------------------------------------------------------------
2097! If split 4D-Var and outer>1, Read in conjugate gradient variables
2098! four outerloop restart.
2099!-----------------------------------------------------------------------
2100!
2101! Open DAV NetCDF file for reading.
2102!
2103 IF (dav(ng)%ncid.eq.-1) THEN
2104 CALL netcdf_open (ng, model, dav(ng)%name, 1, dav(ng)%ncid)
2105 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2106 END IF
2107!
2108! Read in number of converget Ritz eigenvalues
2109!
2110 CALL netcdf_get_ivar (ng, model, dav(ng)%name, &
2111 & 'nConvRitz', nconvritz, &
2112 & ncid = dav(ng)%ncid, &
2113 & start = (/outloop-1/), &
2114 & total = (/1/))
2115 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2116!
2117! Read in Ritz eigenvalues.
2118!
2119 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2120 & 'Ritz', ritz, &
2121 & ncid = dav(ng)%ncid, &
2122 & start = (/1,outloop-1/), &
2123 & total = (/ninner,1/))
2124 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2125!
2126! Read in conjugate gradient "beta" coefficients.
2127!
2128 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2129 & 'cg_beta', cg_beta, &
2130 & ncid = dav(ng)%ncid, &
2131 & start = (/1,1/), &
2132 & total = (/ninner+1,outloop-1/))
2133 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2134!
2135! Read in conjugate gradient "delta" coefficients.
2136!
2137 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2138 & 'cg_delta', cg_delta, &
2139 & ncid = dav(ng)%ncid, &
2140 & start = (/1,1/), &
2141 & total = (/ninner,outloop-1/))
2142 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2143!
2144! Read in normalization coefficients for Lanczos vectors.
2145!
2146 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2147 & 'cg_dla', cg_dla, &
2148 & ncid = dav(ng)%ncid, &
2149 & start = (/1,1/), &
2150 & total = (/ninner,outloop-1/))
2151 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2152!
2153! Read in initial gradient normalization factor.
2154!
2155 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2156 & 'cg_Gnorm_y', cg_gnorm_y, &
2157 & ncid = dav(ng)%ncid, &
2158 & start = (/1/), &
2159 & total = (/outloop-1/))
2160 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2161!
2162 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2163 & 'cg_Gnorm_v', cg_gnorm_v, &
2164 & ncid = dav(ng)%ncid, &
2165 & start = (/1/), &
2166 & total = (/outloop-1/))
2167 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2168!
2169! Read in Lanczos vector normalization factor.
2170!
2171 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2172 & 'cg_QG', cg_qg, &
2173 & ncid = dav(ng)%ncid, &
2174 & start = (/1,1/), &
2175 & total = (/ninner+1,outloop-1/))
2176 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2177!
2178! Read in eigenvalues of Lanczos recurrence relationship.
2179!
2180 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2181 & 'cg_Ritz', cg_ritz, &
2182 & ncid = dav(ng)%ncid, &
2183 & start = (/1,1/), &
2184 & total = (/ninner,outloop-1/))
2185 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2186!
2187! Read eigenvalues relative error.
2188!
2189 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2190 & 'cg_RitzErr', cg_ritzerr, &
2191 & ncid = dav(ng)%ncid, &
2192 & start = (/1,1/), &
2193 & total = (/ninner,outloop-1/))
2194 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2195!
2196! Read in eigenvectors of Lanczos recurrence relationship.
2197!
2198 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2199 & 'cg_zv', cg_zv, &
2200 & ncid = dav(ng)%ncid, &
2201 & start = (/1,1,1/), &
2202 & total = (/ninner,ninner,outloop-1/))
2203 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2204!
2205! Read in initial gradient for minimization.
2206!
2207 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2208 & 'zgrad0', zgrad0, &
2209 & ncid = dav(ng)%ncid, &
2210 & start = (/1,1/), &
2211 & total = (/ndatum(ng)+1,outloop-1/))
2212 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2213!
2214 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2215 & 'vgrad0', vgrad0, &
2216 & ncid = dav(ng)%ncid, &
2217 & start = (/1/), &
2218 & total = (/ndatum(ng)+1/))
2219 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2220!
2221 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2222 & 'Hbk', hbk, &
2223 & ncid = dav(ng)%ncid, &
2224 & start = (/1,1/), &
2225 & total = (/ndatum(ng),outloop-1/))
2226 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2227!
2228 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2229 & 'Jb0', jb0(0:), &
2230 & ncid = dav(ng)%ncid, &
2231 & start = (/1/), &
2232 & total = (/nouter+1/))
2233 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2234!
2235! Read in Lanczos vectors.
2236!
2237 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2238 & 'zcglwk', zcglwk, &
2239 & ncid = dav(ng)%ncid, &
2240 & start = (/1,1,1/), &
2241 & total = (/ndatum(ng)+1,ninner+1,outloop-1/))
2242 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2243!
2244 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2245 & 'vcglwk', vcglwk, &
2246 & ncid = dav(ng)%ncid, &
2247 & start = (/1,1,1/), &
2248 & total = (/ndatum(ng)+1,ninner+1,outloop-1/))
2249 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2250!
2251! Read in converged Lanczos eigenvectors.
2252!
2253 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2254 & 'vcglev', vcglev, &
2255 & ncid = dav(ng)%ncid, &
2256 & start = (/1,1,1/), &
2257 & total = (/ndatum(ng)+1,ninner,outloop-1/))
2258 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2259!
2260! Read in TLmodVal_S.
2261!
2262 CALL netcdf_get_fvar (ng, model, dav(ng)%name, &
2263 & 'TLmodVal_S', tlmodval_s, &
2264 & ncid = dav(ng)%ncid, &
2265 & start = (/1,1,1/), &
2266 & total = (/ndatum(ng),ninner,outloop-1/))
2267 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2268!
2269 RETURN
mod_netcdf::netcdf_get_fvar
Definition mod_netcdf.F:45
mod_netcdf::netcdf_get_ivar
Definition mod_netcdf.F:59
mod_netcdf
Definition mod_netcdf.F:2
mod_netcdf::netcdf_open
subroutine, public netcdf_open(ng, model, ncname, omode, ncid)
Definition mod_netcdf.F:9293

References mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_dla, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_gnorm_y, mod_fourdvar::cg_qg, mod_fourdvar::cg_ritz, mod_fourdvar::cg_ritzerr, mod_fourdvar::cg_zv, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::jb0, mod_fourdvar::nconvritz, mod_fourdvar::ndatum, mod_netcdf::netcdf_open(), mod_scalars::ninner, mod_scalars::noerror, mod_scalars::nouter, mod_fourdvar::ritz, mod_iounits::sourcefile, mod_fourdvar::vcglev, mod_fourdvar::vcglwk, mod_fourdvar::vgrad0, mod_fourdvar::zcglwk, and mod_fourdvar::zgrad0.

Referenced by cg_read_rpcg().

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

◆ cg_read_rpcg_pio()

subroutine, private rpcg_lanczos_mod::cg_read_rpcg_pio ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) outloop )
private

Definition at line 2275 of file rpcg_lanczos.F.

2276!***********************************************************************
2277!
2278 USE mod_pio_netcdf
2279!
2280! Imported variable declarations
2281!
2282 integer, intent(in) :: ng, model, outLoop
2283!
2284! Local variable declarations.
2285!
2286 integer :: status
2287!
2288 character (len=*), parameter :: MyFile = &
2289 & __FILE__//", cg_read_rpcg_pio"
2290!
2291 sourcefile=myfile
2292!
2293!-----------------------------------------------------------------------
2294! If split 4D-Var and outer>1, Read in conjugate gradient variables
2295! four outerloop restart.
2296!-----------------------------------------------------------------------
2297!
2298! Open DAV NetCDF file for reading.
2299!
2300 IF (dav(ng)%pioFile%fh.eq.-1) THEN
2301 CALL pio_netcdf_open (ng, model, dav(ng)%name, 1, &
2302 & dav(ng)%pioFile)
2303 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2304 END IF
2305!
2306! Read in number of converget Ritz eigenvalues
2307!
2308 CALL pio_netcdf_get_ivar (ng, model, dav(ng)%name, &
2309 & 'nConvRitz', nconvritz, &
2310 & piofile = dav(ng)%pioFile, &
2311 & start = (/outloop-1/), &
2312 & total = (/1/))
2313 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2314!
2315! Read in Ritz eigenvalues.
2316!
2317 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2318 & 'Ritz', ritz, &
2319 & piofile = dav(ng)%pioFile, &
2320 & start = (/1,outloop-1/), &
2321 & total = (/ninner,1/))
2322 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2323!
2324! Read in conjugate gradient "beta" coefficients.
2325!
2326 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2327 & 'cg_beta', cg_beta, &
2328 & piofile = dav(ng)%pioFile, &
2329 & start = (/1,1/), &
2330 & total = (/ninner+1,outloop-1/))
2331 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2332!
2333! Read in conjugate gradient "delta" coefficients.
2334!
2335 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2336 & 'cg_delta', cg_delta, &
2337 & piofile = dav(ng)%pioFile, &
2338 & start = (/1,1/), &
2339 & total = (/ninner,outloop-1/))
2340 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2341!
2342! Read in normalization coefficients for Lanczos vectors.
2343!
2344 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2345 & 'cg_dla', cg_dla, &
2346 & piofile = dav(ng)%pioFile, &
2347 & start = (/1,1/), &
2348 & total = (/ninner,outloop-1/))
2349 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2350!
2351! Read in initial gradient normalization factor.
2352!
2353 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2354 & 'cg_Gnorm_y', cg_gnorm_y, &
2355 & piofile = dav(ng)%pioFile, &
2356 & start = (/1/), &
2357 & total = (/outloop-1/))
2358 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2359!
2360 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2361 & 'cg_Gnorm_v', cg_gnorm_v, &
2362 & piofile = dav(ng)%pioFile, &
2363 & start = (/1/), &
2364 & total = (/outloop-1/))
2365 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2366!
2367! Read in Lanczos vector normalization factor.
2368!
2369 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2370 & 'cg_QG', cg_qg, &
2371 & piofile = dav(ng)%pioFile, &
2372 & start = (/1,1/), &
2373 & total = (/ninner+1,outloop-1/))
2374 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2375!
2376! Read in eigenvalues of Lanczos recurrence relationship.
2377!
2378 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2379 & 'cg_Ritz', cg_ritz, &
2380 & piofile = dav(ng)%pioFile, &
2381 & start = (/1,1/), &
2382 & total = (/ninner,outloop-1/))
2383 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2384!
2385! Read eigenvalues relative error.
2386!
2387 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2388 & 'cg_RitzErr', cg_ritzerr, &
2389 & piofile = dav(ng)%pioFile, &
2390 & start = (/1,1/), &
2391 & total = (/ninner,outloop-1/))
2392 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2393!
2394! Read in eigenvectors of Lanczos recurrence relationship.
2395!
2396 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2397 & 'cg_zv', cg_zv, &
2398 & piofile = dav(ng)%pioFile, &
2399 & start = (/1,1,1/), &
2400 & total = (/ninner,ninner,outloop-1/))
2401 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2402!
2403! Read in initial gradient for minimization.
2404!
2405 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2406 & 'zgrad0', zgrad0, &
2407 & piofile = dav(ng)%pioFile, &
2408 & start = (/1,1/), &
2409 & total = (/ndatum(ng)+1,outloop-1/))
2410 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2411!
2412 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2413 & 'vgrad0', vgrad0, &
2414 & piofile = dav(ng)%pioFile, &
2415 & start = (/1/), &
2416 & total = (/ndatum(ng)+1/))
2417 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2418!
2419 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2420 & 'Hbk', hbk, &
2421 & piofile = dav(ng)%pioFile, &
2422 & start = (/1,1/), &
2423 & total = (/ndatum(ng),outloop-1/))
2424 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2425!
2426 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2427 & 'Jb0', jb0(0:), &
2428 & piofile = dav(ng)%pioFile, &
2429 & start = (/1/), &
2430 & total = (/nouter+1/))
2431 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2432!
2433! Read in Lanczos vectors.
2434!
2435 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2436 & 'zcglwk', zcglwk, &
2437 & piofile = dav(ng)%pioFile, &
2438 & start = (/1,1,1/), &
2439 & total = (/ndatum(ng)+1,ninner+1, &
2440 & outloop-1/))
2441 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2442!
2443 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2444 & 'vcglwk', vcglwk, &
2445 & piofile = dav(ng)%pioFile, &
2446 & start = (/1,1,1/), &
2447 & total = (/ndatum(ng)+1,ninner+1, &
2448 & outloop-1/))
2449 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2450!
2451! Read in converged Lanczos eigenvectors.
2452!
2453 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2454 & 'vcglev', vcglev, &
2455 & piofile = dav(ng)%pioFile, &
2456 & start = (/1,1,1/), &
2457 & total = (/ndatum(ng)+1,ninner, &
2458 & outloop-1/))
2459 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2460!
2461! Read in TLmodVal_S.
2462!
2463 CALL pio_netcdf_get_fvar (ng, model, dav(ng)%name, &
2464 & 'TLmodVal_S', tlmodval_s, &
2465 & piofile = dav(ng)%pioFile, &
2466 & start = (/1,1,1/), &
2467 & total = (/ndatum(ng),ninner,outloop-1/))
2468 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2469!
2470 RETURN
mod_pio_netcdf::pio_netcdf_get_fvar
Definition mod_pio_netcdf.F:52
mod_pio_netcdf::pio_netcdf_get_ivar
Definition mod_pio_netcdf.F:66
mod_pio_netcdf
Definition mod_pio_netcdf.F:2
mod_pio_netcdf::pio_netcdf_open
subroutine, public pio_netcdf_open(ng, model, ncname, omode, piofile)
Definition mod_pio_netcdf.F:9148

References mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_dla, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_gnorm_y, mod_fourdvar::cg_qg, mod_fourdvar::cg_ritz, mod_fourdvar::cg_ritzerr, mod_fourdvar::cg_zv, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::jb0, mod_fourdvar::nconvritz, mod_fourdvar::ndatum, mod_scalars::ninner, mod_scalars::noerror, mod_scalars::nouter, mod_pio_netcdf::pio_netcdf_open(), mod_fourdvar::ritz, mod_iounits::sourcefile, mod_fourdvar::vcglev, mod_fourdvar::vcglwk, mod_fourdvar::vgrad0, mod_fourdvar::zcglwk, and mod_fourdvar::zgrad0.

Referenced by cg_read_rpcg().

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

◆ cg_write_rpcg()

subroutine, private rpcg_lanczos_mod::cg_write_rpcg ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) innloop,
integer, intent(in) outloop,
real(dp), intent(in) jf,
real(dp), intent(in) jdata,
real(dp), intent(in) jmod,
real(dp), intent(in) jopt,
real(dp), intent(in) jb,
real(dp), intent(in) jobs,
real(dp), intent(in) jact,
real(dp), intent(in) preducv,
real(dp), intent(in) preducy )
private

Definition at line 1413 of file rpcg_lanczos.F.

1416!***********************************************************************
1417! !
1418! This routine writes conjugate gradient vectors into DAV file !
1419! using either the standard NetCDF library or the Parallel-IO (PIO) !
1420! library. It saves all the variables needed for split or restart !
1421! schemes. !
1422! !
1423!***********************************************************************
1424!
1425! Imported variable declarations
1426!
1427 integer, intent(in) :: ng, model, innLoop, outLoop
1428
1429 real(dp), intent(in) :: Jf, Jdata, Jmod, Jopt, preducv, preducy
1430 real(dp), intent(in) :: Jb, Jobs, Jact
1431!
1432! Local variable declarations.
1433!
1434 character (len=*), parameter :: MyFile = &
1435 & __FILE__//", cg_write_rpcg"
1436!
1437 sourcefile=myfile
1438!
1439!-----------------------------------------------------------------------
1440! Write out conjugate gradient variables.
1441!-----------------------------------------------------------------------
1442!
1443 SELECT CASE (dav(ng)%IOtype)
1444 CASE (io_nf90)
1445 CALL cg_write_rpcg_nf90 (ng, model, innloop, outloop, &
1446 & jf, jdata, jmod, jopt, jb, jobs, &
1447 & jact, preducv, preducy)
1448
1449# if defined PIO_LIB && defined DISTRIBUTE
1450 CASE (io_pio)
1451 CALL cg_write_rpcg_pio (ng, model, innloop, outloop, &
1452 & jf, jdata, jmod, jopt, jb, jobs, &
1453 & jact, preducv, preducy)
1454# endif
1455 CASE DEFAULT
1456 IF (master) WRITE (stdout,10) dav(ng)%IOtype
1457 exit_flag=3
1458 END SELECT
1459 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1460!
1461 10 FORMAT (' CG_WRITE_RPCG - Illegal output file type, io_type = ', &
1462 & i0,/,17x,'Check KeyWord ''OUT_LIB'' in ''roms.in''.')
1463!
1464 RETURN

References cg_write_rpcg_nf90(), cg_write_rpcg_pio(), mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_ncparam::io_nf90, mod_ncparam::io_pio, mod_fourdvar::jobs, mod_parallel::master, mod_scalars::noerror, mod_iounits::sourcefile, and mod_iounits::stdout.

Referenced by rpcg_lanczos().

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

◆ cg_write_rpcg_nf90()

subroutine, private rpcg_lanczos_mod::cg_write_rpcg_nf90 ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) innloop,
integer, intent(in) outloop,
real(dp), intent(in) jf,
real(dp), intent(in) jdata,
real(dp), intent(in) jmod,
real(dp), intent(in) jopt,
real(dp), intent(in) jb,
real(dp), intent(in) jobs,
real(dp), intent(in) jact,
real(dp), intent(in) preducv,
real(dp), intent(in) preducy )
private

Definition at line 1468 of file rpcg_lanczos.F.

1471!***********************************************************************
1472!
1473 USE mod_netcdf
1474!
1475! Imported variable declarations
1476!
1477 integer, intent(in) :: ng, model, innLoop, outLoop
1478
1479 real(dp), intent(in) :: Jf, Jdata, Jmod, Jopt, preducv, preducy
1480 real(dp), intent(in) :: Jb, Jobs, Jact
1481!
1482! Local variable declarations.
1483!
1484 integer :: Nconv, status
1485!
1486 character (len=*), parameter :: MyFile = &
1487 & __FILE__//", cg_write_rpcg_nf90"
1488!
1489 sourcefile=myfile
1490!
1491!-----------------------------------------------------------------------
1492! Write out conjugate gradient variables.
1493!-----------------------------------------------------------------------
1494!
1495! Write out outer and inner iteration.
1496!
1497 CALL netcdf_put_ivar (ng, model, dav(ng)%name, &
1498 & 'outer', outer, &
1499 & (/0/), (/0/), &
1500 & ncid = dav(ng)%ncid)
1501 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1502
1503 CALL netcdf_put_ivar (ng, model, dav(ng)%name, &
1504 & 'inner', inner, &
1505 & (/0/), (/0/), &
1506 & ncid = dav(ng)%ncid)
1507 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1508!
1509! Write out number of converged Ritz eigenvalues.
1510!
1511 IF (innloop.eq.ninner) THEN
1512 CALL netcdf_put_ivar (ng, model, dav(ng)%name, &
1513 & 'nConvRitz', nconvritz, &
1514 & (/outloop/), (/1/), &
1515 & ncid = dav(ng)%ncid)
1516 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1517 END IF
1518!
1519! Write out converged Ritz eigenvalues.
1520!
1521 IF (innloop.eq.ninner) THEN
1522 nconv=max(1,nconvritz)
1523 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1524 & 'Ritz', ritz(1:nconv), &
1525 & (/1,outloop/), (/nconv,1/), &
1526 & ncid = dav(ng)%ncid)
1527 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1528 END IF
1529!
1530! Write out conjugate gradient norm.
1531!
1532 IF (innloop.gt.0) THEN
1533 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1534 & 'cg_beta', cg_beta, &
1535 & (/1,1/), (/ninner+1,nouter/), &
1536 & ncid = dav(ng)%ncid)
1537 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1538 END IF
1539!
1540! Write out Lanczos algorithm coefficients.
1541!
1542 IF (innloop.gt.0) THEN
1543 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1544 & 'cg_delta', cg_delta, &
1545 & (/1,1/), (/ninner,nouter/), &
1546 & ncid = dav(ng)%ncid)
1547 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1548 END IF
1549!
1550! Write normalization coefficients for Lanczos vectors.
1551!
1552 IF (innloop.gt.0) THEN
1553 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1554 & 'cg_dla', cg_dla, &
1555 & (/1,1/), (/ninner,nouter/), &
1556 & ncid = dav(ng)%ncid)
1557 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1558 END IF
1559!
1560! Initial gradient normalization factor.
1561!
1562 IF (innloop.eq.1) THEN
1563 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1564 & 'cg_Gnorm_y', cg_gnorm_y, &
1565 & (/1/), (/nouter/), &
1566 & ncid = dav(ng)%ncid)
1567 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1568
1569 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1570 & 'cg_Gnorm_v', cg_gnorm_v, &
1571 & (/1/), (/nouter/), &
1572 & ncid = dav(ng)%ncid)
1573 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1574 END IF
1575!
1576! Lanczos vector normalization factor.
1577!
1578 IF (innloop.gt.0) THEN
1579 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1580 & 'cg_QG', cg_qg, &
1581 & (/1,1/), (/ninner+1,nouter/), &
1582 & ncid = dav(ng)%ncid)
1583 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1584 END IF
1585!
1586! Reduction in the gradient norm.
1587!
1588 IF (innloop.gt.0) THEN
1589 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1590 & 'cg_Greduc_y', preducy, &
1591 & (/innloop,outloop/), (/1,1/), &
1592 & ncid = dav(ng)%ncid)
1593 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1594
1595 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1596 & 'cg_Greduc_v', preducv, &
1597 & (/innloop,outloop/), (/1,1/), &
1598 & ncid = dav(ng)%ncid)
1599 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1600 END IF
1601!
1602! Eigenvalues of Lanczos recurrence relationship.
1603!
1604 IF (innloop.gt.0) THEN
1605 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1606 & 'cg_Ritz', cg_ritz, &
1607 & (/1,1/), (/ninner,nouter/), &
1608 & ncid = dav(ng)%ncid)
1609 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1610 END IF
1611!
1612! Eigenvalues relative error.
1613!
1614 IF (innloop.gt.0) THEN
1615 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1616 & 'cg_RitzErr', cg_ritzerr, &
1617 & (/1,1/), (/ninner,nouter/), &
1618 & ncid = dav(ng)%ncid)
1619 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1620 END IF
1621!
1622! Eigenvectors of Lanczos recurrence relationship.
1623!
1624 IF (innloop.gt.0) THEN
1625 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1626 & 'cg_zv', cg_zv(:,:,outloop), &
1627 & (/1,1,outloop/), (/ninner,ninner,1/), &
1628 & ncid = dav(ng)%ncid)
1629 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1630 END IF
1631!
1632! First guess initial data misfit.
1633!
1634 IF (innloop.gt.0) THEN
1635 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1636 & 'Jf', jf, &
1637 & (/innloop+1,outloop/), (/1,1/), &
1638 & ncid = dav(ng)%ncid)
1639 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1640!
1641! State estimate data misfit.
1642!
1643 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1644 & 'Jdata', jdata, &
1645 & (/innloop+1,outloop/), (/1,1/), &
1646 & ncid = dav(ng)%ncid)
1647 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1648!
1649! Model penalty function.
1650!
1651 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1652 & 'Jmod', jmod, &
1653 & (/innloop+1,outloop/), (/1,1/), &
1654 & ncid = dav(ng)%ncid)
1655 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1656!
1657! Optimal penalty function.
1658!
1659 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1660 & 'Jopt', jopt, &
1661 & (/innloop+1,outloop/), (/1,1/), &
1662 & ncid = dav(ng)%ncid)
1663 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1664!
1665! Actual model penalty function.
1666!
1667 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1668 & 'Jb', jb, &
1669 & (/innloop+1,outloop/), (/1,1/), &
1670 & ncid = dav(ng)%ncid)
1671 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1672!
1673! Actual data penalty function.
1674!
1675 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1676 & 'Jobs', jobs, &
1677 & (/innloop+1,outloop/), (/1,1/), &
1678 & ncid = dav(ng)%ncid)
1679 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1680!
1681! Actual total penalty function.
1682!
1683 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1684 & 'Jact', jact, &
1685 & (/innloop+1,outloop/), (/1,1/), &
1686 & ncid = dav(ng)%ncid)
1687 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1688 END IF
1689!
1690! Initial gradient for minimization.
1691!
1692 IF (innloop.eq.1) THEN
1693 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1694 & 'zgrad0', zgrad0(:,outloop), &
1695 & (/1,outloop/), (/ndatum(ng)+1,1/), &
1696 & ncid = dav(ng)%ncid)
1697 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1698
1699 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1700 & 'vgrad0', vgrad0, &
1701 & (/1/), (/ndatum(ng)+1/), &
1702 & ncid = dav(ng)%ncid)
1703 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1704
1705 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1706 & 'Hbk', hbk(:,outloop), &
1707 & (/1,outloop/), (/ndatum(ng),1/), &
1708 & ncid = dav(ng)%ncid)
1709 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1710 END IF
1711!
1712! Lanczos vectors.
1713!
1714 IF (innloop.gt.0) THEN
1715 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1716 & 'zcglwk', zcglwk(:,innloop,outloop), &
1717 & (/1,innloop,outloop/), &
1718 & (/ndatum(ng)+1,1,1/), &
1719 & ncid = dav(ng)%ncid)
1720 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1721
1722 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1723 & 'vcglwk', vcglwk(:,innloop,outloop), &
1724 & (/1,innloop,outloop/), &
1725 & (/ndatum(ng)+1,1,1/), &
1726 & ncid = dav(ng)%ncid)
1727 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1728 END IF
1729!
1730! Saved TLmodVal_S.
1731!
1732 CALL netcdf_put_fvar (ng, model, dav(ng)%name, &
1733 & 'TLmodVal_S', &
1734 & tlmodval_s(:,innloop,outloop), &
1735 & (/1,innloop,outloop/), (/ndatum(ng),1,1/), &
1736 & ncid = dav(ng)%ncid)
1737 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1738!
1739!-----------------------------------------------------------------------
1740! Synchronize model/observation NetCDF file to disk.
1741!-----------------------------------------------------------------------
1742!
1743 CALL netcdf_sync (ng, model, dav(ng)%name, dav(ng)%ncid)
1744!
1745 RETURN
mod_netcdf::netcdf_put_fvar
Definition mod_netcdf.F:82
mod_netcdf::netcdf_put_ivar
Definition mod_netcdf.F:96
mod_netcdf::netcdf_sync
subroutine, public netcdf_sync(ng, model, ncname, ncid)
Definition mod_netcdf.F:9505

References mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_dla, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_gnorm_y, mod_fourdvar::cg_qg, mod_fourdvar::cg_ritz, mod_fourdvar::cg_ritzerr, mod_fourdvar::cg_zv, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_scalars::inner, mod_fourdvar::jobs, mod_scalars::nconv, mod_fourdvar::nconvritz, mod_fourdvar::ndatum, mod_netcdf::netcdf_sync(), mod_scalars::ninner, mod_scalars::noerror, mod_scalars::nouter, mod_scalars::outer, mod_fourdvar::ritz, mod_iounits::sourcefile, mod_fourdvar::vcglwk, mod_fourdvar::vgrad0, mod_fourdvar::zcglwk, and mod_fourdvar::zgrad0.

Referenced by cg_write_rpcg().

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

◆ cg_write_rpcg_pio()

subroutine, private rpcg_lanczos_mod::cg_write_rpcg_pio ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) innloop,
integer, intent(in) outloop,
real(dp), intent(in) jf,
real(dp), intent(in) jdata,
real(dp), intent(in) jmod,
real(dp), intent(in) jopt,
real(dp), intent(in) jb,
real(dp), intent(in) jobs,
real(dp), intent(in) jact,
real(dp), intent(in) preducv,
real(dp), intent(in) preducy )
private

Definition at line 1751 of file rpcg_lanczos.F.

1754!***********************************************************************
1755!
1756 USE mod_pio_netcdf
1757!
1758! Imported variable declarations
1759!
1760 integer, intent(in) :: ng, model, innLoop, outLoop
1761
1762 real(dp), intent(in) :: Jf, Jdata, Jmod, Jopt, preducv, preducy
1763 real(dp), intent(in) :: Jb, Jobs, Jact
1764!
1765! Local variable declarations.
1766!
1767 integer :: Nconv, status
1768!
1769 character (len=*), parameter :: MyFile = &
1770 & __FILE__//", cg_write_rpcg_pio"
1771!
1772 sourcefile=myfile
1773!
1774!-----------------------------------------------------------------------
1775! Write out conjugate gradient variables.
1776!-----------------------------------------------------------------------
1777!
1778! Write out outer and inner iteration.
1779!
1780 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, &
1781 & 'outer', outer, &
1782 & (/0/), (/0/), &
1783 & piofile = dav(ng)%pioFile)
1784 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1785
1786 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, &
1787 & 'inner', inner, &
1788 & (/0/), (/0/), &
1789 & piofile = dav(ng)%pioFile)
1790 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1791!
1792! Write out number of converged Ritz eigenvalues.
1793!
1794 IF (innloop.eq.ninner) THEN
1795 CALL pio_netcdf_put_ivar (ng, model, dav(ng)%name, &
1796 & 'nConvRitz', nconvritz, &
1797 & (/outloop/), (/1/), &
1798 & piofile = dav(ng)%pioFile)
1799 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1800 END IF
1801!
1802! Write out converged Ritz eigenvalues.
1803!
1804 IF (innloop.eq.ninner) THEN
1805 nconv=max(1,nconvritz)
1806 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1807 & 'Ritz', ritz(1:nconv), &
1808 & (/1,outloop/), (/nconv,1/), &
1809 & piofile = dav(ng)%pioFile)
1810 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1811 END IF
1812!
1813! Write out conjugate gradient norm.
1814!
1815 IF (innloop.gt.0) THEN
1816 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1817 & 'cg_beta', cg_beta, &
1818 & (/1,1/), (/ninner+1,nouter/), &
1819 & piofile = dav(ng)%pioFile)
1820 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1821 END IF
1822!
1823! Write out Lanczos algorithm coefficients.
1824!
1825 IF (innloop.gt.0) THEN
1826 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1827 & 'cg_delta', cg_delta, &
1828 & (/1,1/), (/ninner,nouter/), &
1829 & piofile = dav(ng)%pioFile)
1830 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1831 END IF
1832!
1833! Write normalization coefficients for Lanczos vectors.
1834!
1835 IF (innloop.gt.0) THEN
1836 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1837 & 'cg_dla', cg_dla, &
1838 & (/1,1/), (/ninner,nouter/), &
1839 & piofile = dav(ng)%pioFile)
1840 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1841 END IF
1842!
1843! Initial gradient normalization factor.
1844!
1845 IF (innloop.eq.1) THEN
1846 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1847 & 'cg_Gnorm_y', cg_gnorm_y, &
1848 & (/1/), (/nouter/), &
1849 & piofile = dav(ng)%pioFile)
1850 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1851
1852 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1853 & 'cg_Gnorm_v', cg_gnorm_v, &
1854 & (/1/), (/nouter/), &
1855 & piofile = dav(ng)%pioFile)
1856 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1857 END IF
1858!
1859! Lanczos vector normalization factor.
1860!
1861 IF (innloop.gt.0) THEN
1862 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1863 & 'cg_QG', cg_qg, &
1864 & (/1,1/), (/ninner+1,nouter/), &
1865 & piofile = dav(ng)%pioFile)
1866 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1867 END IF
1868!
1869! Reduction in the gradient norm.
1870!
1871 IF (innloop.gt.0) THEN
1872 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1873 & 'cg_Greduc_y', preducy, &
1874 & (/innloop,outloop/), (/1,1/), &
1875 & piofile = dav(ng)%pioFile)
1876 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1877
1878 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1879 & 'cg_Greduc_v', preducv, &
1880 & (/innloop,outloop/), (/1,1/), &
1881 & piofile = dav(ng)%pioFile)
1882 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1883 END IF
1884!
1885! Eigenvalues of Lanczos recurrence relationship.
1886!
1887 IF (innloop.gt.0) THEN
1888 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1889 & 'cg_Ritz', cg_ritz, &
1890 & (/1,1/), (/ninner,nouter/), &
1891 & piofile = dav(ng)%pioFile)
1892 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1893 END IF
1894!
1895! Eigenvalues relative error.
1896!
1897 IF (innloop.gt.0) THEN
1898 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1899 & 'cg_RitzErr', cg_ritzerr, &
1900 & (/1,1/), (/ninner,nouter/), &
1901 & piofile = dav(ng)%pioFile)
1902 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1903 END IF
1904!
1905! Eigenvectors of Lanczos recurrence relationship.
1906!
1907 IF (innloop.gt.0) THEN
1908 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1909 & 'cg_zv', cg_zv(:,:,outloop), &
1910 & (/1,1,outloop/), (/ninner,ninner,1/), &
1911 & piofile = dav(ng)%pioFile)
1912 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1913 END IF
1914!
1915! First guess initial data misfit.
1916!
1917 IF (innloop.gt.0) THEN
1918 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1919 & 'Jf', jf, &
1920 & (/innloop+1,outloop/), (/1,1/), &
1921 & piofile = dav(ng)%pioFile)
1922 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1923!
1924! State estimate data misfit.
1925!
1926 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1927 & 'Jdata', jdata, &
1928 & (/innloop+1,outloop/), (/1,1/), &
1929 & piofile = dav(ng)%pioFile)
1930 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1931!
1932! Model penalty function.
1933!
1934 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1935 & 'Jmod', jmod, &
1936 & (/innloop+1,outloop/), (/1,1/), &
1937 & piofile = dav(ng)%pioFile)
1938 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1939!
1940! Optimal penalty function.
1941!
1942 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1943 & 'Jopt', jopt, &
1944 & (/innloop+1,outloop/), (/1,1/), &
1945 & piofile = dav(ng)%pioFile)
1946 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1947!
1948! Actual model penalty function.
1949!
1950 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1951 & 'Jb', jb, &
1952 & (/innloop+1,outloop/), (/1,1/), &
1953 & piofile = dav(ng)%pioFile)
1954 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1955!
1956! Actual data penalty function.
1957!
1958 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1959 & 'Jobs', jobs, &
1960 & (/innloop+1,outloop/), (/1,1/), &
1961 & piofile = dav(ng)%pioFile)
1962 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1963!
1964! Actual total penalty function.
1965!
1966 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1967 & 'Jact', jact, &
1968 & (/innloop+1,outloop/), (/1,1/), &
1969 & piofile = dav(ng)%pioFile)
1970 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1971 END IF
1972!
1973! Initial gradient for minimization.
1974!
1975 IF (innloop.eq.1) THEN
1976 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1977 & 'zgrad0', zgrad0(:,outloop), &
1978 & (/1,outloop/), (/ndatum(ng)+1,1/), &
1979 & piofile = dav(ng)%pioFile)
1980 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1981
1982 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1983 & 'vgrad0', vgrad0, &
1984 & (/1/), (/ndatum(ng)+1/), &
1985 & piofile = dav(ng)%pioFile)
1986 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1987
1988 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1989 & 'Hbk', hbk(:,outloop), &
1990 & (/1,outloop/), (/ndatum(ng),1/), &
1991 & piofile = dav(ng)%pioFile)
1992 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1993 END IF
1994!
1995! Lanczos vectors.
1996!
1997 IF (innloop.gt.0) THEN
1998 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
1999 & 'zcglwk', zcglwk(:,innloop,outloop), &
2000 & (/1,innloop,outloop/), &
2001 & (/ndatum(ng)+1,1,1/), &
2002 & piofile = dav(ng)%pioFile)
2003 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2004
2005 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2006 & 'vcglwk', vcglwk(:,innloop,outloop), &
2007 & (/1,innloop,outloop/), &
2008 & (/ndatum(ng)+1,1,1/), &
2009 & piofile = dav(ng)%pioFile)
2010 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2011 END IF
2012!
2013! Saved TLmodVal_S.
2014!
2015 CALL pio_netcdf_put_fvar (ng, model, dav(ng)%name, &
2016 & 'TLmodVal_S', &
2017 & tlmodval_s(:,innloop,outloop), &
2018 & (/1,innloop,outloop/), &
2019 & (/ndatum(ng),1,1/), &
2020 & piofile = dav(ng)%pioFile)
2021 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
2022!
2023!-----------------------------------------------------------------------
2024! Synchronize model/observation NetCDF file to disk.
2025!-----------------------------------------------------------------------
2026!
2027 CALL pio_netcdf_sync (ng, model, dav(ng)%name, dav(ng)%pioFile)
2028!
2029 RETURN
mod_pio_netcdf::pio_netcdf_put_fvar
Definition mod_pio_netcdf.F:94
mod_pio_netcdf::pio_netcdf_put_ivar
Definition mod_pio_netcdf.F:108
mod_pio_netcdf::pio_netcdf_sync
subroutine, public pio_netcdf_sync(ng, model, ncname, piofile)
Definition mod_pio_netcdf.F:9275

References mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_dla, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_gnorm_y, mod_fourdvar::cg_qg, mod_fourdvar::cg_ritz, mod_fourdvar::cg_ritzerr, mod_fourdvar::cg_zv, mod_iounits::dav, mod_scalars::exit_flag, strings_mod::founderror(), mod_scalars::inner, mod_fourdvar::jobs, mod_scalars::nconv, mod_fourdvar::nconvritz, mod_fourdvar::ndatum, mod_scalars::ninner, mod_scalars::noerror, mod_scalars::nouter, mod_scalars::outer, mod_pio_netcdf::pio_netcdf_sync(), mod_fourdvar::ritz, mod_iounits::sourcefile, mod_fourdvar::vcglwk, mod_fourdvar::vgrad0, mod_fourdvar::zcglwk, and mod_fourdvar::zgrad0.

Referenced by cg_write_rpcg().

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

◆ rpcg_lanczos()

subroutine, public rpcg_lanczos_mod::rpcg_lanczos ( integer, intent(in) ng,
integer, intent(in) model,
integer, intent(in) outloop,
integer, intent(in) innloop,
integer, intent(in) ninnloop,
logical, intent(in) lcgini )

Definition at line 71 of file rpcg_lanczos.F.

73!***********************************************************************
74!
75! Imported variable declarations
76!
77 logical, intent(in) :: Lcgini
78 integer, intent(in) :: ng, model, outLoop, innLoop, NinnLoop
79!
80! Local variable declarations.
81!
82 logical :: Ltrans, Laug
83!
84 integer :: i, ic, j, iobs, ivec, Lscale, info
85!
86 real(dp) :: zbet, eps, preducv, preducy
87 real(dp) :: Jopt, Jf, Jmod, Jdata, Jb, Jobs, Jact, cff
88!
89 real(dp), dimension(NinnLoop) :: zu, zgam, zfact, ztemp3
90 real(dp), dimension(NinnLoop) :: ztemp1, ztemp2, zu1, zu2
91 real(dp), dimension(NinnLoop) :: ztemp4
92 real(dp), dimension(Ndatum(ng)+1) :: px, pgrad, zw, zt
93 real(dp), dimension(Ndatum(ng)+1) :: zhv, zht, zd
94 real(dp), dimension(Ninner,3) :: zwork
95 real(dp), dimension(2*(NinnLoop-1)) :: work
96 real(dp), dimension(Ninner,Ninner) :: zgv
97!
98 character (len=13) :: string
99
100 character (len=*), parameter :: MyFile = &
101 & __FILE__
102!
103!=======================================================================
104! Weak constraint 4D-Var conjugate gradient, Lanczos-based algorithm.
105!=======================================================================
106
107# ifdef PROFILE
108!
109! Turn on time clock.
110!
111 CALL wclock_on (ng, model, 85, __line__, myfile)
112# endif
113!
114! This conjugate gradient algorithm is not run in parallel since the
115! weak constraint is done in observation space. Mostly all the
116! variables are 1D arrays. Therefore, in parallel applications (only
117! distributed-memory is possible) the master node does the computations
118! and then broadcasts results to remaining nodes in the communicator.
119!
120! This algorithm solves A(x+x0)=b for x, by minimizing
121! J=0.5*x'Ax-x'(b+Ax0), where x0 is a first-guess estimate of the
122! representer coefficients from the previous outer-loop.
123! For the first outer-loop, x0=0. In the code x0=cg_pxsave and
124! x=px.
125!
126 ltrans=.false.
127!
128! Initialize internal parameters. This is needed here for output
129! reasons.
130!
131 jopt=0.0_dp
132 jf=0.0_dp
133 jdata=0.0_dp
134 jmod=0.0_dp
135 jb=jb0(outloop-1)
136 jobs=0.0_dp
137 jact=0.0_dp
138 ritzmaxerr=hevecerr
139 eps=0.0_dp
140 preducv=0.0_dp
141 preducy=0.0_dp
142 zwork=0.0_dp
143 zgv=0.0_dp
144
145!
146 master_thread : IF (master) THEN
147!
148! Multiply TLmodVal by ObsScale to effectively remove any rejected
149! observations from the analysis.
150!
151 IF (innloop.gt.0) THEN
152 DO iobs=1,ndatum(ng)
153 tlmodval(iobs)=obsscale(iobs)*tlmodval(iobs)
154 END DO
155 END IF
156!
157!-----------------------------------------------------------------------
158! Initialization step, innloop = 0.
159!-----------------------------------------------------------------------
160!
161 minimize : IF (innloop.eq.0) THEN
162
163# if defined RBL4DVAR || \
164 defined rbl4dvar_ana_sensitivity || \
165 defined rbl4dvar_fct_sensitivity || \
166 defined tl_rbl4dvar
167!
168 DO iobs=1,ndatum(ng)
169 bckmodval(iobs)=nlmodval(iobs)
170 END DO
171# endif
172!
173! For the first outer-loop, x0=0.
174!
175 IF (outloop.eq.1) THEN
176 DO iobs=1,ndatum(ng)+1
177 px(iobs)=0.0_dp
178 END DO
179 DO iobs=1,ndatum(ng)
180 hbk(iobs,outloop)=0.0_dp
181 END DO
182 ELSE
183!! NOTE: CODE WILL NOT WORK FOR R4D-Var AT THE MOMENT SINCE THIS IS
184!! THE WRONG TLmodVal HERE!!!
185 DO iobs=1,ndatum(ng)
186 hbk(iobs,outloop)=-tlmodval(iobs)
187 END DO
188 IF (lprecond) THEN
189 DO iobs=1,ndatum(ng)+1
190 fourdvar(ng)%cg_Dpxsum(iobs)=fourdvar(ng)%cg_pxsum(iobs)
191 END DO
192 lscale=1 ! Spectral LMP
193 laug=.false.
194 CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, &
195 & fourdvar(ng)%cg_Dpxsum)
196 DO iobs=1,ndatum(ng)+1
197 fourdvar(ng)%cg_Dpxsum(iobs)= &
198 & -fourdvar(ng)%cg_Dpxsum(iobs)+fourdvar(ng)%cg_pxsum(iobs)
199 END DO
200 END IF
201 END IF
202!
203! If this is the first pass of the inner loop, set up the vectors and
204! store the first guess. The initial starting guess is assumed to be
205! zero in which case the residual is just: (obs-model).
206!
207 DO iobs=1,ndatum(ng)
208!
209! Selime code: r=b.
210!
211# if defined RBL4DVAR || \
212 defined rbl4dvar_ana_sensitivity || \
213 defined rbl4dvar_fct_sensitivity || \
214 defined tl_rbl4dvar
215 pgrad(iobs)=obsscale(iobs)* &
216 & (obsval(iobs)-bckmodval(iobs))
217# else
218 pgrad(iobs)=obsscale(iobs)* &
219 & (obsval(iobs)-tlmodval(iobs))
220# endif
221 IF (obserr(iobs).ne.0.0_r8) THEN
222 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
223 END IF
224 vgrad0(iobs)=pgrad(iobs)
225 vgrad0s(iobs)=vgrad0(iobs)
226 END DO
227!
228! AUGMENTED
229! Augment the residual vector.
230!
231 pgrad(ndatum(ng)+1)=1.0_dp
232 vgrad0(ndatum(ng)+1)=pgrad(ndatum(ng)+1)
233 vgrad0s(ndatum(ng)+1)=vgrad0(ndatum(ng)+1)
234!
235! If preconditioning, convert pgrad from v-space to y-space.
236!
237! Selime code: z1=G*r.
238!
239 IF (lprecond.and.(outloop.gt.1)) THEN
240 lscale=1 ! Spectral LMP
241 laug=.true.
242 CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, &
243 & pgrad)
244 END IF
245!
246! Selime code: zgrad0=z1.
247!
248! Selime code: vgrad0=r.
249!
250 DO iobs=1,ndatum(ng)+1
251 zgrad0(iobs,outloop)=pgrad(iobs)
252 END DO
253!
254! Selime code: ADmodVal=z1.
255!
256 DO iobs=1,ndatum(ng)
257 admodval(iobs)=zgrad0(iobs,outloop)
258 END DO
259!
260 preducv=1.0_dp
261 preducy=1.0_dp
262!
263!-----------------------------------------------------------------------
264! Minimization step, innLoop > 0.
265!-----------------------------------------------------------------------
266!
267 ELSE
268!
269! COMPLETE THE CALCULATION OF PREVIOUS LANCZOS VECTOR.
270!
271! Selime's code: t1=GDG'*z1=TLmodVal.
272!
273! Complete the computation of the Lanczos vector from previous inner-loop
274!
275 IF (innloop.EQ.1) THEN
276 cg_gnorm_v(outloop)=0.0_dp
277 DO iobs=1,ndatum(ng)
278 cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+ &
279 & tlmodval(iobs)*vgrad0(iobs)
280 END DO
281!
282! AUGMENTED. Here Hbk=H(x0-xb). Hbk=0 and Jb0=0 when outer=1.
283! Jb0 is computed using sum_i (G'*w)_i, i.e. the sum of the outputs'
284! of the adjoint model at the end of AD_LOOP2.
285!
286 DO iobs=1,ndatum(ng)
287 cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+ &
288 & hbk(iobs,outloop)* &
289 & zgrad0(ndatum(ng)+1,outloop)* &
290 & vgrad0(iobs)+ &
291 & hbk(iobs,outloop)* &
292 & vgrad0(ndatum(ng)+1)* &
293 & zgrad0(iobs,outloop)
294 END DO
295 cg_gnorm_v(outloop)=cg_gnorm_v(outloop)+ &
296 & jb0(outloop-1)*vgrad0(ndatum(ng)+1)* &
297 & zgrad0(ndatum(ng)+1,outloop)
298!
299 cg_gnorm_v(outloop)=sqrt(cg_gnorm_v(outloop))
300!
301! Selime's code: beta_first=sqrt(r't1)=cg_Gnorm_v.
302! cg_QG(1)=beta_first
303!
304 cg_qg(1,outloop)=cg_gnorm_v(outloop)
305!
306! Normalize TLmodVal here and save in TLmodVal_S.
307!
308! AUGMENTED.
309!
310 DO iobs=1,ndatum(ng)+1
311!
312! Selime code: v1=r/beta_first=zcglwk(1).
313!
314 zcglwk(iobs,1,outloop)=vgrad0(iobs)/ &
315 & cg_gnorm_v(outloop)
316!
317! Selime code: z1=z1/beta_first=vcglwk(1).
318!
319 vcglwk(iobs,1,outloop)=zgrad0(iobs,outloop)/ &
320 & cg_gnorm_v(outloop)
321 END DO
322!
323! AUGMENTED
324!
325 DO iobs=1,ndatum(ng)
326!
327! Selime code: t1=t1/beta_first=TLmodVal
328!
329 tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
330 tlmodval(iobs)=tlmodval(iobs)/cg_gnorm_v(outloop)
331 END DO
332!
333 preducv=1.0_dp
334 preducy=1.0_dp
335!
336! Compute the actual penalty function terms.
337!
338 jobs=0.0_dp
339!
340! NO AUGMENTED TERM HERE.
341!
342 jact=jb
343 DO iobs=1,ndatum(ng)
344 IF (obserr(iobs).ne.0.0_r8) THEN
345 jact=jact+obserr(iobs)*vgrad0s(iobs)*vgrad0s(iobs)
346 END IF
347 END DO
348 jobs=jact-jb
349!
350 ELSE
351!
352! First grab the non-normalized Lanczos vector.
353!
354! AUGMENTED.
355!
356 DO iobs=1,ndatum(ng)+1
357!
358! Selime code: w=pgrad, t1=TLmodVal.
359!
360 pgrad(iobs)=zcglwk(iobs,innloop,outloop)
361 END DO
362!
363! Selime's code: beta=pgrad*HBH'*pgrad=cg_beta
364!
365 cg_beta(innloop,outloop)=0.0_dp
366 DO iobs=1,ndatum(ng)
367 cg_beta(innloop,outloop)=cg_beta(innloop,outloop)+ &
368 & pgrad(iobs)*tlmodval(iobs)
369 END DO
370!
371! AUGMENTED.
372!
373 DO iobs=1,ndatum(ng)
374 cg_beta(innloop,outloop)=cg_beta(innloop,outloop)+ &
375 & hbk(iobs,outloop)* &
376 & vcglwk(ndatum(ng)+1,innloop, &
377 & outloop)* &
378 & pgrad(iobs)+ &
379 & hbk(iobs,outloop)* &
380 & pgrad(ndatum(ng)+1)* &
381 & vcglwk(iobs,innloop,outloop)
382 END DO
383 cg_beta(innloop,outloop)=cg_beta(innloop,outloop)+ &
384 & jb0(outloop-1)* &
385 & pgrad(ndatum(ng)+1)* &
386 & vcglwk(ndatum(ng)+1,innloop, &
387 & outloop)
388!
389 cg_beta(innloop,outloop)=sqrt(cg_beta(innloop,outloop))
390!
391! Normalize TLmodVal here and save in TLmodVal_S.
392!
393! AUGMENTED
394!
395 DO iobs=1,ndatum(ng)+1
396!
397! Selime code: v1=w/beta=zcglwk
398!
399 zcglwk(iobs,innloop,outloop)=pgrad(iobs)/ &
400 & cg_beta(innloop,outloop)
401!
402! Selime code: z1=w/beta=vcglwk
403!
404 vcglwk(iobs,innloop,outloop)=vcglwk(iobs,innloop,outloop)/&
405 & cg_beta(innloop,outloop)
406 END DO
407!
408! Selime code: t1=t1/beta=TLmodVal
409!
410! AUGMENTED
411!
412 DO iobs=1,ndatum(ng)
413 tlmodval_s(iobs,innloop,outloop)=tlmodval(iobs)
414 tlmodval(iobs)=tlmodval(iobs)/cg_beta(innloop,outloop)
415 END DO
416!
417! Selime's code: cg_QG=zgrad0*HBH'*zcglwk
418!
419 cg_qg(innloop,outloop)=0.0_dp
420 DO iobs=1,ndatum(ng)
421 cg_qg(innloop,outloop)=cg_qg(innloop,outloop)+ &
422 & tlmodval(iobs)* &
423 & vgrad0(iobs)
424 END DO
425!
426! AUGMENTED.
427!
428 DO iobs=1,ndatum(ng)
429 cg_qg(innloop,outloop)=cg_qg(innloop,outloop)+ &
430 & hbk(iobs,outloop)* &
431 & vcglwk(ndatum(ng)+1,innloop, &
432 & outloop)* &
433 & vgrad0(iobs)+ &
434 & hbk(iobs,outloop)* &
435 & vcglwk(iobs,innloop,outloop)* &
436 & vgrad0(ndatum(ng)+1)
437 END DO
438 cg_qg(innloop,outloop)=cg_qg(innloop,outloop)+ &
439 & jb0(outloop-1)* &
440 & vcglwk(ndatum(ng)+1,innloop, &
441 & outloop)* &
442 & vgrad0(ndatum(ng)+1)
443!TEST
444!TEST IF (innLoop.gt.1) cg_QG(innLoop,outLoop)=0.0_dp
445!TEST
446!
447! Calculate the new solution based upon the new, orthonormalized
448! Lanczos vector. First, the tridiagonal system is solved by
449! decomposition and forward/back substitution.
450!
451 zbet=cg_delta(1,outloop)
452 zu(1)=cg_qg(1,outloop)/zbet
453 DO ivec=2,innloop-1
454 zgam(ivec)=cg_beta(ivec,outloop)/zbet
455 zbet=cg_delta(ivec,outloop)-cg_beta(ivec,outloop)* &
456 & zgam(ivec)
457 zu(ivec)=(cg_qg(ivec,outloop)-cg_beta(ivec,outloop)* &
458 & zu(ivec-1))/zbet
459 END DO
460 zwork(innloop-1,3)=zu(innloop-1)
461
462 DO ivec=innloop-2,1,-1
463 zu(ivec)=zu(ivec)-zgam(ivec+1)*zu(ivec+1)
464 zwork(ivec,3)=zu(ivec)
465 END DO
466!
467! AUGMENTED.
468!
469 DO iobs=1,ndatum(ng)+1
470 zw(iobs)=zgrad0(iobs,outloop)+ &
471 & cg_beta(innloop,outloop)* &
472 & vcglwk(iobs,innloop,outloop)* &
473 & zwork(innloop-1,3)
474 END DO
475
476!
477! AUGMENTED.
478!
479 DO iobs=1,ndatum(ng)+1
480 px(iobs)=0.0_dp
481 DO ivec=1,innloop-1
482 px(iobs)=px(iobs)+ &
483 & zcglwk(iobs,ivec,outloop)*zwork(ivec,3)
484 zw(iobs)=zw(iobs)+ &
485 & zcglwk(iobs,ivec,outloop)*cg_qg(ivec,outloop)
486 END DO
487 END DO
488!
489! If preconditioning, convert px from y-space to v-space.
490! We will always keep px in v-space.
491!
492 IF (lprecond.and.(outloop.gt.1)) THEN
493 lscale=1 ! Spectral LMP
494 laug=.true.
495 CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, px)
496 END IF
497!
498! Compute the reduction in the gradient Ax-b (in y-space).
499!
500!AMM preducy=0.0_dp
501!AMM DO iobs=1,Ndatum(ng)
502!AMM preducy=preducy+zw(iobs)*zw(iobs)
503!AMM END DO
504!AMM preducy=SQRT(preducy)/cg_Gnorm_y(outLoop)
505!
506! Compute the reduction in the gradient Ax-b (in v-space).
507!
508 IF (lprecond.and.(outloop.gt.1)) THEN
509 lscale=1 ! Spectral LMP
510 laug=.true.
511 CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, zw)
512 END IF
513!
514 preducv=0.0_dp
515!
516! AUGMENTED.
517!
518 DO iobs=1,ndatum(ng)+1
519 preducv=preducv+zw(iobs)*zw(iobs)
520 END DO
521 preducv=sqrt(preducv)/cg_gnorm_v(outloop)
522!
523! Estimate the residual (in y-space) by:
524!
525! ||Ax-b|| = cg_beta(k+1,outLoop)*zwork(innLoop,3)
526!
527 eps=abs(cg_beta(innloop,outloop)*zwork(innloop-1,3))
528!
529! Estimate the various penalty function components based on
530! Bennett (2002), page 22 and pages 42-44:
531!
532! Jopt = value of the penalty function when true px
533! has been identified.
534! Jf = initial data misfit for the first-guess
535! Jdata = data misfit for the state estimate
536! Jmod = the model penalty function (the sum of the
537! dynamical, initial and boundary penalties)
538! Jb = the ACTUAL model penalty function.
539! Jobs = the ACTUAL data penalty function.
540! Jact = the ACTUAL total penalty function.
541!
542! NOTE: Unlike I4D-Var where we compute the actual cost function
543! for each iteration that decreases with increasing
544! number of iterations, the various penalty terms
545! represent the optimal values (i.e. the values when
546! optimal state estimate has been identified). Therefore
547! the various penalty terms will only be correct when the
548! optimal state has been found. Therefore, as the
549! R4D-Var proceeds, Jopt, Jdata and Jmod will asymptote
550! to their final values and WILL NOT necessary decrease
551! with increasing iteration number.
552!
553 IF (outloop.eq.1) THEN
554 DO iobs=1,ndatum(ng)
555 IF (obserr(iobs).ne.0.0_r8) THEN
556! Jopt=Jopt+px(iobs)*vgrad0(iobs)/SQRT(ObsErr(iobs))
557 jopt=jopt+px(iobs)*tlmodval_s(iobs,1,outloop)
558! Jf=Jf+vgrad0(iobs)*vgrad0(iobs)/ObsErr(iobs)
559 jf=jf+vgrad0(iobs)*tlmodval_s(iobs,1,outloop)
560 END IF
561!AMM: For Jdata we need GDG'*px but we don't have this except at the end
562!AMM Jdata=Jdata+px(iobs)*px(iobs)
563 END DO
564!AMM Jmod=Jopt-Jdata
565 ELSE
566 DO iobs=1,ndatum(ng)
567!AMM Jopt=Jopt- &
568!AMM & (px(iobs)+cg_pxsave(iobs))*cg_innov(iobs)
569 IF (obserr(iobs).ne.0.0_r8) THEN
570 jf=jf+cg_innov(iobs)*cg_innov(iobs)
571 END IF
572!AMM Jdata=Jdata+ &
573!AMM & (px(iobs)+cg_pxsave(iobs))* &
574!AMM & (px(iobs)+cg_pxsave(iobs))
575 END DO
576 jmod=jopt-jdata
577 END IF
578!
579! Compute the actual penalty function terms.
580!
581 DO ivec=1,innloop-1
582 IF (ivec.eq.1) THEN
583 zfact(ivec)=1.0_dp/cg_gnorm_v(outloop)
584 ELSE
585 zfact(ivec)=1.0_dp/cg_beta(ivec,outloop)
586 ENDIF
587 END DO
588!
589! AUGMENTED.
590!
591 DO iobs=1,ndatum(ng)+1
592 zd(iobs)=vgrad0s(iobs)
593 END DO
594!
595! AUGMENTED.
596!
597 DO iobs=1,ndatum(ng)+1
598 zhv(iobs)=zd(iobs)
599 END DO
600 DO ivec=1,innloop-1
601 ztemp1(ivec)=0.0_dp
602 ztemp3(ivec)=0.0_dp
603 ztemp4(ivec)=0.0_dp
604 DO iobs=1,ndatum(ng)
605 ztemp1(ivec)=ztemp1(ivec)+ &
606 & tlmodval_s(iobs,ivec,outloop)* &
607 & zhv(iobs)* &
608 & zfact(ivec)
609 ztemp3(ivec)=ztemp3(ivec)+ &
610 & vcglwk(iobs,ivec,outloop)* &
611 & hbk(iobs,outloop)
612 ztemp4(ivec)=ztemp4(ivec)+ &
613 & tlmodval_s(iobs,ivec,outloop)* &
614 & zgrad0(iobs,outloop)* &
615 & zfact(ivec)
616 END DO
617!
618! AUGMENTED
619! Note the factor zfact is not needed in the next loop since it has
620! already been applied to vcglwk above.
621!
622 DO iobs=1,ndatum(ng)
623 ztemp1(ivec)=ztemp1(ivec)+ &
624 & zhv(iobs)* &
625 & hbk(iobs,outloop)* &
626 & vcglwk(ndatum(ng)+1,ivec,outloop)+ &
627 & hbk(iobs,outloop)* &
628 & vcglwk(iobs,ivec,outloop)* &
629 & zhv(ndatum(ng)+1)
630 ztemp4(ivec)=ztemp4(ivec)+ &
631 & zgrad0(iobs,outloop)* &
632 & hbk(iobs,outloop)* &
633 & vcglwk(ndatum(ng)+1,ivec,outloop)+ &
634 & hbk(iobs,outloop)* &
635 & vcglwk(iobs,ivec,outloop)* &
636 & zgrad0(ndatum(ng)+1,outloop)
637 END DO
638 ztemp1(ivec)=ztemp1(ivec)+ &
639 & zhv(ndatum(ng)+1)* &
640 & jb0(outloop-1)* &
641 & vcglwk(ndatum(ng)+1,ivec,outloop)
642 ztemp3(ivec)=ztemp3(ivec)+ &
643 & jb0(outloop-1)* &
644 & vcglwk(ndatum(ng)+1,ivec,outloop)
645 ztemp4(ivec)=ztemp4(ivec)+ &
646 & zgrad0(ndatum(ng)+1,outloop)* &
647 & jb0(outloop-1)* &
648 & vcglwk(ndatum(ng)+1,ivec,outloop)
649!TEST
650 ztemp4(ivec)=0.0_dp
651!TEST
652 END DO
653 zbet=cg_delta(1,outloop)
654!TEST zu1(1)=ztemp1(1)/zbet
655!TEST zu1(1)=ztemp4(1)/zbet
656 zu1(1)=cg_qg(1,outloop)/zbet
657!TEST
658 DO ivec=2,innloop-1
659 zgam(ivec)=cg_beta(ivec,outloop)/zbet
660 zbet=cg_delta(ivec,outloop)- &
661 & cg_beta(ivec,outloop)*zgam(ivec)
662!TEST zu1(ivec)=(ztemp1(ivec)- &
663 zu1(ivec)=(ztemp4(ivec)- &
664 & cg_beta(ivec,outloop)*zu1(ivec-1))/zbet
665 END DO
666 ztemp2(innloop-1)=zu1(innloop-1)
667
668 DO ivec=innloop-2,1,-1
669 zu1(ivec)=zu1(ivec)-zgam(ivec+1)*zu1(ivec+1)
670 ztemp2(ivec)=zu1(ivec)
671 END DO
672!
673 jb=jb0(outloop-1)
674 jact=jb
675 DO iobs=1,ndatum(ng)
676 IF (obserr(iobs).ne.0.0_r8) THEN
677 jact=jact+obserr(iobs)*vgrad0s(iobs)*vgrad0s(iobs)
678 END IF
679 END DO
680 DO ivec=1,innloop-1
681 jb=jb+ztemp2(ivec)*ztemp2(ivec)
682!TEST Jact=Jact-ztemp1(ivec)*ztemp2(ivec)+ &
683! & 2.0_dp*ztemp2(ivec)*ztemp3(ivec)
684 jact=jact-ztemp1(ivec)*ztemp2(ivec)
685 END DO
686 DO iobs=1,ndatum(ng)
687 jb=jb-2.0_dp*px(iobs)*hbk(iobs,outloop)
688!TEST
689 jact=jact+2.0_dp*px(iobs)*hbk(iobs,outloop)
690!TEST
691 END DO
692 jb=jb-2.0_dp*px(ndatum(ng)+1)*jb0(outloop-1)
693!TEST
694 jact=jact+2.0_dp*px(ndatum(ng)+1)*jb0(outloop-1)
695!TEST
696 jobs=jact-jb
697
698 END IF
699!
700! NEXT CALCULATION
701!
702! After the initialization, for every other inner loop, calculate a
703! new Lanczos vector, store it in the matrix, and update search.
704!
705! AUGMENTED.
706!
707 DO iobs=1,ndatum(ng)
708!
709! Selime code: w=t1.
710!
711 pgrad(iobs)=tlmodval(iobs)+ &
712 & hbk(iobs,outloop)* &
713 & vcglwk(ndatum(ng)+1,innloop,outloop)
714 END DO
715 DO iobs=1,ndatum(ng)
716!
717! Convert gradient contribution from x-space to v-space.
718!
719 IF (obserr(iobs).ne.0.0_r8) THEN
720!
721! Selime code: w=t1/R. No augmented term here since Raug^-1=[R^-1 0].
722!
723 pgrad(iobs)=pgrad(iobs)/obserr(iobs)
724 END IF
725 END DO
726 pgrad(ndatum(ng)+1)=0.0_dp
727!
728! Selime code: w=t1/R+z1
729!
730 DO iobs=1,ndatum(ng)
731 pgrad(iobs)=pgrad(iobs)+vcglwk(iobs,innloop,outloop)
732 END DO
733 pgrad(ndatum(ng)+1)=pgrad(ndatum(ng)+1)+ &
734 & vcglwk(ndatum(ng)+1,innloop,outloop)
735!
736 IF (innloop.gt.1) THEN
737!
738! Selime code: w=w-v0*beta
739!
740! AUGMENTED.
741!
742 DO iobs=1,ndatum(ng)+1
743 pgrad(iobs)=pgrad(iobs)- &
744 & cg_beta(innloop,outloop)* &
745 & zcglwk(iobs,innloop-1,outloop)
746 END DO
747 END IF
748!
749 cg_delta(innloop,outloop)=0.0_dp
750!
751! Selime's code: alpha=w'*t1=cg_delta
752!
753! AUGMENTED.
754!
755 DO iobs=1,ndatum(ng)
756 cg_delta(innloop,outloop)=cg_delta(innloop,outloop)+ &
757 & tlmodval(iobs)* &
758 & pgrad(iobs)+ &
759 & hbk(iobs,outloop)* &
760 & vcglwk(ndatum(ng)+1,innloop, &
761 & outloop)* &
762 & pgrad(iobs)+ &
763 & hbk(iobs,outloop)* &
764 & vcglwk(iobs,innloop,outloop)* &
765 & pgrad(ndatum(ng)+1)
766 END DO
767 cg_delta(innloop,outloop)=cg_delta(innloop,outloop)+ &
768 & jb0(outloop-1)* &
769 & vcglwk(ndatum(ng)+1,innloop, &
770 & outloop)* &
771 & pgrad(ndatum(ng)+1)
772!
773! Exit, if not a positive definite algorithm.
774!
775 IF (cg_delta(innloop,outloop).le.0.0_dp) THEN
776 WRITE (stdout,20) cg_delta(innloop,outloop), outloop, &
777 & innloop
778 exit_flag=8
779 GO TO 10
780 END IF
781!
782! Selime code: w=w-alpha*v1
783!
784! AUGMENTED.
785!
786 DO iobs=1,ndatum(ng)+1
787 pgrad(iobs)=pgrad(iobs)-cg_delta(innloop,outloop)* &
788 & zcglwk(iobs,innloop,outloop)
789 END DO
790!
791! Orthonormalize against previous directions.
792!
793! Identify the appropriate Lanczos vector normalization coefficient.
794!
795 DO ivec=1,innloop
796 IF (ivec.eq.1) THEN
797 zfact(ivec)=1.0_dp/cg_gnorm_v(outloop)
798 ELSE
799 zfact(ivec)=1.0_dp/cg_beta(ivec,outloop)
800 ENDIF
801 END DO
802!
803! Selime's code: cg_dla=pgrad*HBH'*v1=pgrad*TLmodVal_S/zfact
804!
805! AUGMENTED
806!
807 DO ivec=innloop,1,-1
808 cg_dla(ivec,outloop)=0.0_dp
809 DO iobs=1,ndatum(ng)
810 cg_dla(ivec,outloop)=cg_dla(ivec,outloop)+ &
811 & pgrad(iobs)* &
812 & tlmodval_s(iobs,ivec,outloop)* &
813 & zfact(ivec)+pgrad(iobs)* &
814 & hbk(iobs,outloop)* &
815 & vcglwk(ndatum(ng)+1,ivec,outloop)+ &
816 & hbk(iobs,outloop)* &
817 & vcglwk(iobs,ivec,outloop)* &
818 & pgrad(ndatum(ng)+1)
819 END DO
820 cg_dla(ivec,outloop)=cg_dla(ivec,outloop)+ &
821 & jb0(outloop-1)* &
822 & vcglwk(ndatum(ng)+1,ivec,outloop)* &
823 & pgrad(ndatum(ng)+1)
824 DO iobs=1,ndatum(ng)+1
825 pgrad(iobs)=pgrad(iobs)- &
826 & cg_dla(ivec,outloop)* &
827 & vcglwk(iobs,ivec,outloop)
828 END DO
829 END DO
830!
831! Save the non-normalized Lanczos vector. zcglwk is used as temporary
832! storage. The calculation will be completed at the start of the next
833! inner-loop.
834!
835! Selime code: w=zcglwk.
836!
837! AUGMENTED
838!
839 DO iobs=1,ndatum(ng)+1
840 zcglwk(iobs,innloop+1,outloop)=pgrad(iobs)
841 END DO
842!
843! If preconditioning, convert pgrad from v-space to y-space.
844!
845! Selime code: z1=Gw.
846!
847 IF (lprecond.and.(outloop.gt.1)) THEN
848 lscale=1 ! Spectral LMP
849 laug=.true.
850 CALL rprecond (ng, lscale, laug, outloop, ninnloop-1, pgrad)
851 END IF
852!
853! Selime code: z1=vcglwk.
854!
855! AUGMENTED.
856!
857 DO iobs=1,ndatum(ng)+1
858 vcglwk(iobs,innloop+1,outloop)=pgrad(iobs)
859 END DO
860!
861! Put the new trial solution into the adjoint vector for the next
862! loop or put the final solution into the adjoint vector on the
863! final inner-loop.
864!
865 IF (innloop.eq.ninnloop) THEN
866!
867! Note: px is already in v-space so there is no need to convert
868! if preconditioning; cg_pxsave is also in v-space.
869!
870 DO iobs=1,ndatum(ng)
871 admodval(iobs)=px(iobs)
872 END DO
873 DO iobs=1,ndatum(ng)+1
874 fourdvar(ng)%cg_pxsum(iobs)=fourdvar(ng)%cg_pxsum(iobs)+ &
875 & fourdvar(ng)%cg_pxsave(iobs)
876 fourdvar(ng)%cg_pxsave(iobs)=px(iobs)
877!AMM FOURDVAR(ng)%cg_pxsum(iobs)=FOURDVAR(ng)%cg_pxsum(iobs)+ &
878!AMM & px(iobs)
879 END DO
880 ELSE
881!
882! Selime code: z1=ADmodVal.
883!
884 DO iobs=1,ndatum(ng)
885 admodval(iobs)=vcglwk(iobs,innloop+1,outloop)
886 END DO
887!
888 END IF
889!
890 END IF minimize
891!
892!-----------------------------------------------------------------------
893! Report minimization progress.
894!-----------------------------------------------------------------------
895!
896 IF (innloop.gt.0) THEN
897 WRITE (stdout,30) ndatum(ng), &
898 & outloop, innloop-1, eps, &
899 & outloop, innloop-1, preducy, &
900 & outloop, innloop-1, preducv, &
901 & outloop, innloop-1, jf, &
902 & outloop, innloop-1, jdata, &
903 & outloop, innloop-1, jmod, &
904 & outloop, innloop-1, jopt, &
905 & outloop, innloop-1, jb, &
906 & outloop, innloop-1, jobs, &
907 & outloop, innloop-1, jact, &
908 & outloop, innloop-1
909 END IF
910 DO ivec=1,innloop
911 WRITE (stdout,40) ivec, cg_delta(ivec,outloop), &
912 & cg_beta(ivec,outloop), zwork(ivec,3)
913 END DO
914 WRITE (stdout,50) innloop+1, cg_beta(innloop+1,outloop)
915!
916! If preconditioning, report the Ritz eigenvalues used.
917!
918 IF ((lhessianev.or.lprecond).and.(outloop.gt.1)) THEN
919 IF (nritzev.gt.0) THEN
920 WRITE (stdout,60) outloop, innloop, nconvritz
921 ic=0
922 DO ivec=1,ninnloop-1
923 IF (ivec.gt.(ninnloop-1-nconvritz)) THEN
924 string=' '
925 ic=ic+1
926 WRITE (stdout,70) ivec, cg_ritz(ivec,outloop-1), &
927 & cg_ritzerr(ivec,outloop-1), &
928 & trim(adjustl(string)), &
929 & 'Used=', ic
930 ELSE
931 string=' '
932 WRITE (stdout,80) ivec, cg_ritz(ivec,outloop-1), &
933 & cg_ritzerr(ivec,outloop-1), &
934 & trim(adjustl(string))
935 END IF
936 END DO
937 ELSE
938 WRITE (stdout,90) outloop, innloop, ritzmaxerr
939 ic=0
940 DO ivec=1,ninnloop
941 IF (cg_ritzerr(ivec,outloop-1).le.ritzmaxerr) THEN
942 string='converged'
943 ic=ic+1
944 WRITE (stdout,70) ivec, cg_ritz(ivec,outloop-1), &
945 & cg_ritzerr(ivec,outloop-1), &
946 & trim(adjustl(string)), &
947 & 'Good=', ic
948 ELSE
949 string='not converged'
950 WRITE (stdout,80) ivec, cg_ritz(ivec,outloop-1), &
951 & cg_ritzerr(ivec,outloop-1), &
952 & trim(adjustl(string))
953 END IF
954 END DO
955 END IF
956 END IF
957!
958!-----------------------------------------------------------------------
959! If last inner loop, innLoop = NinnLoop.
960!-----------------------------------------------------------------------
961!
962! Compute the eigenvalues and eigenvectors of the tridiagonal matrix.
963!
964 IF (innloop.eq.ninnloop) THEN
965 IF (lhessianev.or.lprecond) THEN
966 DO ivec=1,innloop-1
967 cg_ritz(ivec,outloop)=cg_delta(ivec,outloop)
968 END DO
969 DO ivec=1,innloop-2
970 zwork(ivec,1)=cg_beta(ivec+1,outloop)
971 END DO
972!
973! Use the LAPACK routine DSTEQR to compute the eigenvectors and
974! eigenvalues of the tridiagonal matrix. If applicable, the
975! eigenpairs are computed by master thread only.
976!
977 CALL dsteqr ('I', innloop-1, cg_ritz(1,outloop), zwork(1,1),&
978 & zgv, ninner, work, info)
979 IF (info.ne.0) THEN
980 WRITE (stdout,*) ' RPCG_LANCZOS - Error in DSTEQR:', &
981 & ' info = ',info
982 exit_flag=8
983 GO TO 10
984 END IF
985!
986 DO j=1,ninner
987 DO i=1,ninner
988 cg_zv(i,j,outloop)=zgv(i,j)
989 END DO
990 END DO
991!
992! Estimate the Ritz value error bounds.
993!
994 DO ivec=1,innloop-1
995 cg_ritzerr(ivec,outloop)=abs(cg_beta(innloop,outloop)* &
996 & cg_zv(innloop-1,ivec,outloop))
997 END DO
998!
999! Check for exploding or negative Ritz values.
1000!
1001 DO ivec=1,innloop-1
1002 IF (cg_ritzerr(ivec,outloop).lt.0.0_dp) THEN
1003 WRITE (stdout,*) ' RPCG_LANCZOS - negative Ritz value', &
1004 & ' found.'
1005 exit_flag=8
1006 GO TO 10
1007 END IF
1008 END DO
1009!
1010! Change the error bounds for the acceptable eigenvectors.
1011!
1012 DO ivec=1,ninnloop-1
1013 cg_ritzerr(ivec,outloop)=cg_ritzerr(ivec,outloop)/ &
1014 & cg_ritz(ninnloop-1,outloop)
1015 END DO
1016!
1017! Report Eigenvalues and their relative accuracy.
1018!
1019 WRITE (stdout,100) outloop, innloop, ritzmaxerr
1020 ic=0
1021 DO ivec=1,ninnloop-1
1022 IF (cg_ritzerr(ivec,outloop).le.ritzmaxerr) THEN
1023 string='converged'
1024 ic=ic+1
1025 WRITE (stdout,70) ivec, cg_ritz(ivec,outloop), &
1026 & cg_ritzerr(ivec,outloop), &
1027 & trim(adjustl(string)), &
1028 & 'Good=', ic
1029 ELSE
1030 string='not converged'
1031 WRITE (stdout,80) ivec, cg_ritz(ivec,outloop), &
1032 & cg_ritzerr(ivec,outloop), &
1033 & trim(adjustl(string))
1034 END IF
1035 END DO
1036!
1037! Calculate the converged eigenvectors (vcglev) of the [R_n + Cobs].
1038!
1039 CALL rpevecs (ng, outloop, ninnloop-1)
1040 END IF
1041 END IF
1042
1043 END IF master_thread
1044!
1045!-----------------------------------------------------------------------
1046! Come here if not a possite definite algorithm.
1047!-----------------------------------------------------------------------
1048!
1049 10 CONTINUE
1050# ifdef DISTRIBUTE
1051 CALL mp_bcasti (ng, model, exit_flag)
1052# endif
1053 IF (founderror(exit_flag, noerror, __line__, myfile)) RETURN
1054
1055# ifdef DISTRIBUTE
1056!
1057!-----------------------------------------------------------------------
1058! Broadcast new solution to other nodes.
1059!-----------------------------------------------------------------------
1060!
1061 CALL mp_bcastf (ng, model, admodval)
1062# if defined RBL4DVAR || defined R4DVAR || \
1063 defined sensitivity_4dvar || defined tl_rbl4dvar || \
1064 defined tl_r4dvar
1065!! CALL mp_bcastf (ng, model, TLmodVal_S(:,:,outLoop)) ! not needed
1066 CALL mp_bcastf (ng, model, hbk(:,outloop))
1067 CALL mp_bcastf (ng, model, jb0(outloop-1))
1068# endif
1069 CALL mp_bcasti (ng, model, info)
1070 CALL mp_bcastf (ng, model, cg_beta(:,outloop))
1071 CALL mp_bcastf (ng, model, cg_qg(:,outloop))
1072 CALL mp_bcastf (ng, model, cg_delta(:,outloop))
1073 CALL mp_bcastf (ng, model, cg_dla(:,outloop))
1074 CALL mp_bcastf (ng, model, cg_gnorm_y(:))
1075 CALL mp_bcastf (ng, model, cg_gnorm_v(:))
1076 CALL mp_bcastf (ng, model, fourdvar(ng)%cg_pxsave(:))
1077!AMM CALL mp_bcastf (ng, model, cg_pxsave(:))
1078 CALL mp_bcastf (ng, model, cg_innov(:))
1079 CALL mp_bcastf (ng, model, zgrad0(:,outloop))
1080 CALL mp_bcastf (ng, model, zcglwk(:,:,outloop))
1081 CALL mp_bcastf (ng, model, vcglwk(:,:,outloop))
1082 IF ((lhessianev.or.lprecond).and.(innloop.eq.ninnloop)) THEN
1083 CALL mp_bcastf (ng, model, cg_ritz(:,outloop))
1084 CALL mp_bcastf (ng, model, cg_ritzerr(:,outloop))
1085 CALL mp_bcastf (ng, model, cg_zv(:,:,outloop))
1086 CALL mp_bcastf (ng, model, vcglev(:,:,outloop))
1087 END IF
1088 CALL mp_bcastf (ng, model, jf)
1089 CALL mp_bcastf (ng, model, jdata)
1090 CALL mp_bcastf (ng, model, jmod)
1091 CALL mp_bcastf (ng, model, jopt)
1092 CALL mp_bcastf (ng, model, jobs)
1093 CALL mp_bcastf (ng, model, jact)
1094 CALL mp_bcastf (ng, model, preducv)
1095 CALL mp_bcastf (ng, model, preducy)
1096# endif
1097!
1098!-----------------------------------------------------------------------
1099! Write out conjugate gradient vectors into 4D-Var NetCDF file.
1100!-----------------------------------------------------------------------
1101!
1102# ifdef BGQC
1103!
1104! Save the cost function values for background quality control decisions
1105!
1106 jact_s(innloop)=jact
1107!
1108# endif
1109 IF (innloop.gt.0) THEN
1110 CALL cg_write_rpcg (ng, model, innloop, outloop, &
1111 & jf, jdata, jmod, jopt, jb, jobs, jact, &
1112 & preducv, preducy)
1113 END IF
1114
1115# ifdef PROFILE
1116!
1117! Turn off time clock.
1118!
1119 CALL wclock_off (ng, model, 85, __line__, myfile)
1120# endif
1121!
1122 20 FORMAT (/,' RPCG_LANCZOS - Fatal error, not possitive definite', &
1123 & ' algorithm:',/, &
1124 & /,11x,'cg_delta = ',1p,e15.8,0p,3x,'(',i3.3,', ',i3.3,')')
1125 30 FORMAT (/,' RPCG_LANCZOS - Conjugate Gradient Information:',/, &
1126 & /,11x,'Ndatum = ',i0,/, /, &
1127 & 1x,'(',i3.3,',',i3.3,'): ', &
1128 & 'Residual estimate, eps = ', &
1129 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1130 & 'Reduction in gradient norm, Greduc y-space = ', &
1131 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1132 & 'Reduction in gradient norm, Greduc v-space = ', &
1133 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1134 & 'First guess initial data misfit, Jf = ', &
1135 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1136 & 'State estimate data misfit, Jdata = ', &
1137 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1138 & 'Model penalty function, Jmod = ', &
1139 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1140 & 'Optimal penalty function, Jopt = ', &
1141 & 1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ', &
1142 & 'Actual Model penalty function, Jb = ', &
1143 & 1p,e14.7,/,1x,'(',i3.3,',',i3.3,'): ', &
1144 & 'Actual data penalty function, Jobs = ', &
1145 & 1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ', &
1146 & 'Actual total penalty function, Jact = ', &
1147 & 1p,e14.7,/,/,1x,'(',i3.3,',',i3.3,'): ', &
1148 & 'Lanczos vectors - cg_delta, cg_beta, zwork:',/)
1149 40 FORMAT (6x,i3.3,4x,3(1p,e15.8,5x))
1150 50 FORMAT (6x,i3.3,24x,1p,e15.8)
1151 60 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1152 & 'Ritz eigenvalues used in preconditioning, ', &
1153 & 'nConvRitz = ',i3.3,/)
1154 70 FORMAT (6x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a,5x,'('a,i3.3,')')
1155 80 FORMAT (6x,i3.3,2x,1p,e14.7,2x,1p,e14.7,2x,a)
1156 90 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1157 & 'Ritz eigenvalues used in preconditioning, ', &
1158 & 'RitzMaxErr = ',1p,e12.5,/)
1159100 FORMAT (/,1x,'(',i3.3,',',i3.3,'): ', &
1160 & 'New Ritz eigenvalues and their accuracy, ', &
1161 & 'RitzMaxErr = ',1p,e12.5,/)
1162!
1163 RETURN
wclock_off
recursive subroutine wclock_off(ng, model, region, line, routine)
Definition timers.F:148
wclock_on
recursive subroutine wclock_on(ng, model, region, line, routine)
Definition timers.F:3

References mod_fourdvar::admodval, mod_fourdvar::cg_beta, mod_fourdvar::cg_delta, mod_fourdvar::cg_dla, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_gnorm_y, mod_fourdvar::cg_innov, mod_fourdvar::cg_qg, mod_fourdvar::cg_ritz, mod_fourdvar::cg_ritzerr, cg_write_rpcg(), mod_fourdvar::cg_zv, lapack_mod::dsteqr(), mod_scalars::exit_flag, strings_mod::founderror(), mod_fourdvar::fourdvar, mod_fourdvar::hevecerr, mod_fourdvar::jact_s, mod_fourdvar::jb0, mod_fourdvar::jobs, mod_fourdvar::lhessianev, mod_fourdvar::lprecond, mod_parallel::master, mod_fourdvar::nconvritz, mod_fourdvar::ndatum, mod_scalars::ninner, mod_fourdvar::nlmodval, mod_scalars::noerror, mod_fourdvar::nritzev, mod_fourdvar::obserr, mod_fourdvar::obsscale, mod_fourdvar::obsval, mod_fourdvar::ritzmaxerr, rpevecs(), rprecond(), mod_iounits::stdout, mod_fourdvar::vcglev, mod_fourdvar::vcglwk, mod_fourdvar::vgrad0, mod_fourdvar::vgrad0s, wclock_off(), wclock_on(), mod_fourdvar::zcglwk, and mod_fourdvar::zgrad0.

Referenced by rbl4dvar_mod::increment().

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

◆ rpevecs()

subroutine, private rpcg_lanczos_mod::rpevecs ( integer, intent(in) ng,
integer, intent(in) outloop,
integer, intent(in) ninnloop )
private

Definition at line 1167 of file rpcg_lanczos.F.

1168!***********************************************************************
1169! !
1170! This routine computes the converged eigenvectors of (R_n+Cobs). !
1171! !
1172!***********************************************************************
1173!
1174! Imported variable declarations
1175!
1176 integer, intent(in) :: ng, outLoop, NinnLoop
1177!
1178! Local variable declarations.
1179!
1180 integer :: iobs, ivec, jn, jk, ij, jkk, ingood
1181!
1182 real(dp) :: dla, zfact
1183 real(dp), dimension(Ndatum(ng)+1) :: zt
1184!
1185!-----------------------------------------------------------------------
1186! Compute and store converged eigenvectors of (R_n+Cobs).
1187!-----------------------------------------------------------------------
1188!
1189! Count and collect the converged eigenvalues.
1190!
1191 ingood=0
1192 DO ivec=ninnloop,1,-1
1193 ingood=ingood+1
1194 ritz(ingood)=cg_ritz(ivec,outloop)
1195 END DO
1196!
1197! Calculate the converged eigenvectors of (R_n+Cobs).
1198!
1199 ij=0
1200 DO jn=ninnloop,1,-1
1201 ij=ij+1
1202! AUGMENTED
1203 DO iobs=1,ndatum(ng)
1204 vcglev(iobs,ij,outloop)=0.0_dp
1205 END DO
1206 DO jk=1,ninnloop
1207 DO iobs=1,ndatum(ng)
1208 vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)+ &
1209 & (vcglwk(iobs,jk,outloop)- &
1210 & fourdvar(ng)%cg_pxsum(iobs)* &
1211 & vcglwk(ndatum(ng)+1,jk,outloop))* &
1212 & cg_zv(jk,jn,outloop)
1213 END DO
1214 END DO
1215!
1216! Orthonormalize.
1217!
1218 DO jk=1,ij-1
1219 DO iobs=1,ndatum(ng)
1220 zt(iobs)=0.0_dp
1221 END DO
1222 DO iobs=1,ndatum(ng)
1223 DO jkk=1,ninnloop
1224 IF (jkk.eq.1) THEN
1225 zfact=1.0_dp/cg_gnorm_v(outloop)
1226 ELSE
1227 zfact=1.0_dp/cg_beta(jkk,outloop)
1228 ENDIF
1229 zt(iobs)=zt(iobs)+ &
1230 & (tlmodval_s(iobs,jkk,outloop)* &
1231 & zfact+ &
1232 & hbk(iobs,outloop)* &
1233 & vcglwk(ndatum(ng)+1,jkk,outloop))* &
1234 & cg_zv(jkk,jk,outloop)
1235 END DO
1236 END DO
1237 dla=0.0_dp
1238 DO iobs=1,ndatum(ng)
1239 dla=dla+zt(iobs)*vcglev(iobs,ij,outloop)
1240 END DO
1241 DO iobs=1,ndatum(ng)
1242 vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)- &
1243 & dla*zt(iobs)
1244 END DO
1245 END DO
1246 dla=0.0_dp
1247! AUGMENTED
1248 DO iobs=1,ndatum(ng)
1249 dla=dla+vcglev(iobs,ij,outloop)*vcglev(iobs,ij,outloop)
1250 END DO
1251! AUGMENTED+1
1252 DO iobs=1,ndatum(ng)
1253 vcglev(iobs,ij,outloop)=vcglev(iobs,ij,outloop)/sqrt(dla)
1254 END DO
1255 END DO
1256!
1257 RETURN

References mod_fourdvar::cg_beta, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_ritz, mod_fourdvar::cg_zv, mod_fourdvar::fourdvar, mod_fourdvar::ndatum, mod_fourdvar::ritz, mod_fourdvar::vcglev, and mod_fourdvar::vcglwk.

Referenced by rpcg_lanczos().

Here is the caller graph for this function:

◆ rprecond()

subroutine, private rpcg_lanczos_mod::rprecond ( integer, intent(in) ng,
integer, intent(in) lscale,
logical, intent(in) laug,
integer, intent(in) outloop,
integer, intent(in) ninnloop,
real(dp), dimension(ndatum(ng)+1), intent(inout) py )
private

Definition at line 1261 of file rpcg_lanczos.F.

1262!***********************************************************************
1263! !
1264! This routine is the preconditioner in observation space. !
1265! !
1266!***********************************************************************
1267!
1268! Imported variable declarations
1269!
1270 logical, intent(in) :: Laug
1271 integer, intent(in) :: ng, Lscale, outLoop, NinnLoop
1272!
1273 real(dp), intent(inout) :: py(Ndatum(ng)+1) ! AUGMENTED
1274!
1275! Local variable declarations.
1276!
1277 integer :: is, ie, inc, iss, i, jk
1278 integer :: nol, nols, nole, ninc
1279 integer :: ingood
1280 integer :: iobs, nvec, ivec
1281!
1282 real(dp) :: dla, dlar, fac, facritz, zfact
1283 real(dp), dimension(NinnLoop) :: zvect
1284 real(dp), dimension(Ndatum(ng)+1) :: zt
1285!
1286! Set the do-loop indices for the sequential preconditioner
1287! loop.
1288!
1289 nols=1
1290 nole=outloop-1
1291 ninc=1
1292!
1293! Apply the preconditioners derived from all previous outer-loops
1294! sequentially.
1295!
1296 DO nol=nols,nole,ninc
1297!
1298! Determine the number of Ritz vectors to use.
1299! For Lritz=.TRUE., choose HvecErr to be larger.
1300!
1301 ingood=0
1302 DO i=1,ninnloop
1303 IF (cg_ritzerr(i,nol).le.ritzmaxerr) THEN
1304 ingood=ingood+1
1305 END IF
1306 END DO
1307 IF (nritzev.gt.0) THEN
1308 nconvritz=nritzev
1309 ingood=nconvritz
1310 ELSE
1311 nconvritz=ingood
1312 END IF
1313!
1314 IF (lscale.gt.0) THEN
1315 is=1
1316 ie=ingood
1317 inc=1
1318 ELSE
1319 is=ingood
1320 ie=1
1321 inc=-1
1322 END IF
1323!
1324! Lscale determines the form of the preconditioner:
1325!
1326! 1= LMP
1327! -1= Inverse LMP
1328!
1329 DO nvec=is,ie,inc
1330 DO iobs=1,ndatum(ng)
1331 zt(iobs)=0.0_dp
1332 END DO
1333 DO iobs=1,ndatum(ng)
1334 DO jk=1,ninnloop
1335 IF (jk.eq.1) THEN
1336 zfact=1.0_dp/cg_gnorm_v(nol)
1337 ELSE
1338 zfact=1.0_dp/cg_beta(jk,nol)
1339 ENDIF
1340 zt(iobs)=zt(iobs)+ &
1341 & (tlmodval_s(iobs,jk,nol)* &
1342 & zfact+ &
1343 & hbk(iobs,nol)* &
1344 & vcglwk(ndatum(ng)+1,jk,nol))* &
1345 & cg_zv(jk,nvec,nol)
1346 END DO
1347 END DO
1348!
1349 dla=0.0_dp
1350 DO iobs=1,ndatum(ng)
1351 dla=dla+zt(iobs)*py(iobs)
1352 END DO
1353 IF (lritz) THEN
1354 dlar=0.0_dp
1355 DO iobs=1,ndatum(ng)
1356 dlar=dlar+ &
1357 & (tlmodval_s(iobs,ninnloop+1,nol)/ &
1358 & cg_beta(ninnloop+1,nol)+ &
1359 & hbk(iobs,nol)* &
1360 & vcglwk(ndatum(ng)+1,ninnloop+1,nol))* &
1361 & py(iobs)
1362 END DO
1363 END IF
1364!
1365! NOTE: Lscale=1 and Lscale=-1 cases are not complete yet.
1366!
1367 IF (lscale.eq.-1) THEN
1368 fac=(cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
1369 ELSE
1370 fac=(1.0_dp/cg_ritz(ninnloop+1-nvec,nol)-1.0_dp)*dla
1371 END IF
1372!
1373 DO iobs=1,ndatum(ng)
1374 py(iobs)=py(iobs)+fac*vcglev(iobs,nvec,nol)
1375 END DO
1376!
1377 IF (lritz) THEN
1378!
1379! NOTE: We still need the code for Lscale=-1.
1380!
1381 IF (lscale.eq.1) THEN
1382 facritz=cg_beta(ninnloop+1,nol)* &
1383 & cg_zv(ninnloop,ninnloop+1-nvec,nol)/ &
1384 & cg_ritz(ninnloop+1-nvec,nol)
1385 DO iobs=1,ndatum(ng)
1386 py(iobs)=py(iobs)+ &
1387 & facritz*facritz*vcglev(iobs,nvec,nol)- &
1388 & facritz*dlar* &
1389 & (vcglwk(iobs,ninnloop+1,nol)- &
1390 & fourdvar(ng)%cg_pxsum(iobs)* &
1391 & vcglwk(ndatum(ng)+1,ninnloop+1,nol))
1392 END DO
1393 DO iobs=1,ndatum(ng)
1394 py(iobs)=py(iobs)-facritz*dla*vcglev(iobs,nvec,nol)
1395 END DO
1396 END IF
1397 END IF
1398!
1399 END DO
1400 END DO
1401!
1402 IF (laug) THEN
1403 DO iobs=1,ndatum(ng)
1404 py(iobs)=py(iobs)+ &
1405 & fourdvar(ng)%cg_Dpxsum(iobs)*py(ndatum(ng)+1)
1406 END DO
1407 END IF
1408!
1409 RETURN

References mod_fourdvar::cg_beta, mod_fourdvar::cg_gnorm_v, mod_fourdvar::cg_ritz, mod_fourdvar::cg_ritzerr, mod_fourdvar::cg_zv, mod_fourdvar::fourdvar, mod_fourdvar::lritz, mod_fourdvar::nconvritz, mod_fourdvar::ndatum, mod_fourdvar::nritzev, mod_fourdvar::ritzmaxerr, mod_fourdvar::vcglev, and mod_fourdvar::vcglwk.

Referenced by rpcg_lanczos().

Here is the caller graph for this function: