# Presentation: 045424

## GSI Code

- Russ Treadon
- NOAA/NWS/NCEP/EMC

- GSI User Orientation
- 4-5 January 2005

- Room 209, WWB
- Camp Springs, MD

## Caveats / Reminders

- Describing code as it is,
- not necessarily as it should or will be
- Development discussion is encouraged  later, not now
- All GSI Orientation presentations will be available on the web
- There’s a lot to cover in a short time
- Can’t cover everything
- Information content will be too much / too little for various users

## GSI Code - General Flow

- User input &
- initializations

- Read in &
- distribute
- observations

- Additional
- initializations

- Outer loop
- (Read/distribute) guess & derived fields
- (Read/distribute) background error
- Compute/qc O-G & prepare for inner loop
- Call inner loop
- Compute gradient information
- Apply background error
- Compute search direction
- Compute step size
- Update analysis increment

- Write analysis & related output

## GSI Main driver

- Code level:  gsimain.F90
- Start up
- Start MPI
- Allocate arrays, initialize variables/arrays
- Read namelist input
- Namelist dependent allocations/initializations
- Pass execution onto lower level routines  gsisub.f90
- Shut down
- Deallocate arrays
- Stop MPI

- User input &
- initializations

## GSI namelists
currently three
SETUP
OBS_INPUT
SINGLEOB_TEST
future code updates will better organize namelist elements into logical groups
e.g., 
quality control parameters
grid information
background error parameters
diagnostic output parameters
... 

- currently three
- SETUP
- OBS_INPUT
- SINGLEOB_TEST
- future code updates will better organize namelist elements into logical groups
- e.g.,
- quality control parameters
- grid information
- background error parameters
- diagnostic output parameters
- ...

- GSI namelists

## Contents
Grid specifications
Background error parameters
Iteration and diagnostic control
(nonlinear) Quality control parameters
Miscellaneous parameters

- Grid specifications
- Background error parameters
- Iteration and diagnostic control
- (nonlinear) Quality control parameters
- Miscellaneous parameters

- GSI namelist SETUP

## GSI namelist SETUP

- Grid information
- jcap:  spectral resolution (global only)
- nlat,nlon,nsig:  x-, y-, z- dimensions
- NOTE:  regional gets grid dimension from guess file
- hybrid:  .*true*.=hybrid z coordinate
- regional:  .*true*.=regional run, .*false*.=global

- &SETUP
-    ***JCAP***=$JCAP,***NLAT***=$NLAT,***NLON***=$LONA,***nsig***=$LEVS,
-    miter=2,niter(1)=100,niter(2)=100,niter_no_qc(1)=30,niter_no_qc(2)=0,
-    write_diag(1)=.true.,write_diag(2)=.false.,write_diag(3)=.true.,
-    as=0.55,0.55,0.85,0.95,1.20,1.20,vs=0.6,hzscl=0.75,1.50,3.00,bw=0.0,
-    diag_rad=.true.,diag_pcp=.false.,diag_ozone=.true.,diag_conv=.true.,
-    ermint=1.3,ermaxt=5.6,grosst=10.0,b_t=10.0,pg_t=0.007,gencode=78,
-    dfact=0.1,factqmin=0.8,factqmax=1.2,deltim=$DELTIM,
-    ndat=29,jpch=493,npred=5,iguess=-1,
-    ***hybrid***=.false.,***regional***=.false.,oneobtest=.false.,
- /

## GSI namelist SETUP

- Background error parameters
- as:  background error variance factors
- as(1,2,3,4,5,6,7,8) = , , ln(ps), Tv, q/qsg, rlyr_oz, SST, rcld
- hzscl: horizontal scale stretching coefficients
- hzwgt: horizontal scale weighting factors
- vs: vertical scale stretching coefficients

- &SETUP
-    JCAP=$JCAP,NLAT=$NLAT,NLON=$LONA,nsig=$LEVS,
-    miter=2,niter(1)=100,niter(2)=100,niter_no_qc(1)=30,niter_no_qc(2)=0,
-    write_diag(1)=.true.,write_diag(2)=.false.,write_diag(3)=.true.,
-    ***as***=0.55,0.55,0.85,0.95,1.20,1.20,***vs***=0.6,***hzscl***=0.75,1.50,3.00,bw=0.0,
-    diag_rad=.true.,diag_pcp=.false.,diag_ozone=.true.,diag_conv=.true.,
-    ermint=1.3,ermaxt=5.6,grosst=10.0,b_t=10.0,pg_t=0.007,gencode=78,
-    dfact=0.1,factqmin=0.8,factqmax=1.2,deltim=$DELTIM,
-    ndat=29,jpch=493,npred=5,iguess=-1,
-    hybrid=.false.,regional=.false.,oneobtest=.false.,
- /

## GSI namelist SETUP

- Iteration and diagnostic control
- miter:  number of outer iterations
- niter:  number of inner iterations per outer loop
- diag_conv, diag_ozone, diag_pcp, diag_rad:  logical to turn on (.*true*.  default) / off (.*false*.) generation of diagnostic output files

- &SETUP
-    JCAP=$JCAP,NLAT=$NLAT,NLON=$LONA,nsig=$LEVS,
-    ***miter***=2,***niter(1)***=100,***niter(2)***=100,niter_no_qc(1)=30,niter_no_qc(2)=0,
-    write_diag(1)=.true.,write_diag(2)=.false.,write_diag(3)=.true.,
-    as=0.55,0.55,0.85,0.95,1.20,1.20,vs=0.6,hzscl=0.75,1.50,3.00,bw=0.0,
-    ***diag_rad***=.true.,***diag_pcp***=.false.,***diag_ozone***=.true.,***diag_conv***=.true.,
-    ermint=1.3,ermaxt=5.6,grosst=10.0,b_t=10.0,pg_t=0.007,gencode=78,
-    dfact=0.1,factqmin=0.8,factqmax=1.2,deltim=$DELTIM,
-    ndat=29,jpch=493,npred=5,iguess=-1,
-    hybrid=.false.,regional=.false.,oneobtest=.false.,
- /

## GSI namelist SETUP

- Iteration and diagnostic control (continued)
- write_diag: .*true*. = write diagnostic file on indicated iteration
- setting write_diag(miter+1)=.*true*. generates O-A statistics

- &SETUP
-    JCAP=$JCAP,NLAT=$NLAT,NLON=$LONA,nsig=$LEVS,
-    miter=2,niter(1)=100,niter(2)=100,niter_no_qc(1)=30,niter_no_qc(2)=0,
-    ***write_diag(1)***=.true.,***write_diag(2)***=.false.,***write_diag(3)***=.true.,
-    as=0.55,0.55,0.85,0.95,1.20,1.20,vs=0.6,hzscl=0.75,1.50,3.00,bw=0.0,
-    diag_rad=.true.,diag_pcp=.false.,diag_ozone=.true.,diag_conv=.true.,
-    ermint=1.3,ermaxt=5.6,grosst=10.0,b_t=10.0,pg_t=0.007,gencode=78,
-    dfact=0.1,factqmin=0.8,factqmax=1.2,deltim=$DELTIM,
-    ndat=29,jpch=493,npred=5,iguess=-1,
-    hybrid=.false.,regional=.false.,oneobtest=.false.,
- /

## GSI namelist SETUP

- Quality control parameters (Tv example)
- ermint,ermaxt: min/max values for Tv observation error
- grosst:  (Obs-Ges)/Error threshold above which observation will be tossed by gross check

- &SETUP
-    JCAP=$JCAP,NLAT=$NLAT,NLON=$LONA,nsig=$LEVS,
-    miter=2,niter(1)=100,niter(2)=100,niter_no_qc(1)=30,niter_no_qc(2)=0,
-    write_diag(1)=.true.,write_diag(2)=.false.,write_diag(3)=.true.,
-    as=0.55,0.55,0.85,0.95,1.20,1.20,vs=0.6,hzscl=0.75,1.50,3.00,bw=0.0,
-    diag_rad=.true.,diag_pcp=.false.,diag_ozone=.true.,diag_conv=.true.,
-    ***ermint***=1.3,***ermaxt***=5.6,***grosst***=10.0,b_t=10.0,pg_t=0.007,gencode=78,
-    dfact=0.1,factqmin=0.8,factqmax=1.2,deltim=$DELTIM,
-    ndat=29,jpch=493,npred=5,iguess=-1,
-    hybrid=.false.,regional=.false.,oneobtest=.false.,
- /

## GSI namelist SETUP

- Nonlinear quality control parameters (Tv example)
- b_t:  range for gross errors normalized by observation error
- pg_t:  probability of gross errors for Tv
- niter_no_qc:  number of inner loop iterations without nonlinear quality control

- &SETUP
-    JCAP=$JCAP,NLAT=$NLAT,NLON=$LONA,nsig=$LEVS,
-    miter=2,niter(1)=100,niter(2)=100,***niter_no_qc(1)***=30,***niter_no_qc(2)***=0,
-    write_diag(1)=.true.,write_diag(2)=.false.,write_diag(3)=.true.,
-    as=0.55,0.55,0.85,0.95,1.20,1.20,vs=0.6,hzscl=0.75,1.50,3.00,bw=0.0,
-    diag_rad=.true.,diag_pcp=.false.,diag_ozone=.true.,diag_conv=.true.,
-    ermint=1.3,ermaxt=5.6,grosst=10.0,***b_t***=10.0,***pg_t***=0.007,gencode=78,
-    dfact=0.1,factqmin=0.8,factqmax=1.2,deltim=$DELTIM,
-    ndat=29,jpch=493,npred=5,iguess=-1,
-    hybrid=.false.,regional=.false.,oneobtest=.false.,
- /

## GSI namelist SETUP

- Miscellaneous parameters (just a few listed below)
- factqmin,factqmax:  scaling factors for negative / supersaturation q constraint (0.0 = no constraint)
- ndat:  number of observation files to read
- jpch:  number of channels (Tb assimilation) over all satellites/sensors  number of entries in satinfo file
- iguess:  flag for use/generation of guess solution file

- &SETUP
-    JCAP=$JCAP,NLAT=$NLAT,NLON=$LONA,nsig=$LEVS,
-    miter=2,niter(1)=100,niter(2)=100,niter_no_qc(1)=30,niter_no_qc(2)=0,
-    write_diag(1)=.true.,write_diag(2)=.false.,write_diag(3)=.true.,
-    as=0.55,0.55,0.85,0.95,1.20,1.20,vs=0.6,hzscl=0.75,1.50,3.00,bw=0.0,
-    diag_rad=.true.,diag_pcp=.false.,diag_ozone=.true.,diag_conv=.true.,
-    ermint=1.3,ermaxt=5.6,grosst=10.0,b_t=10.0,pg_t=0.007,gencode=78,
-    dfact=0.1,***factqmin***=0.8,***factqmax***=1.2,deltim=$DELTIM,
-    ***ndat***=29,***jpch***=493,npred=5,***iguess***=-1,
-    hybrid=.false.,regional=.false.,oneobtest=.false.,
- /

## Contents – applicable to all obs types
dfile:  local name of data file
dtype:  character string indicating obs type
ps, t, q, pw, uv, spd, rw, srw, dw, gps_ref, sst, sbuv2, pcp_ssm/I, pcp_tmi, amsua, amsub, hirs/2, hirs/3, msu, goes, goesimg, airs, eos_amsua, hsb, avhrr_navy, avhrr  see read_obs.f90
id:        satellite id (Tb and pcp obs)
dload:  computational load factor

- dfile:  local name of data file
- dtype:  character string indicating obs type
- ps, t, q, pw, uv, spd, rw, srw, dw, gps_ref, sst, sbuv2, pcp_ssm/I, pcp_tmi, amsua, amsub, hirs/2, hirs/3, msu, goes, goesimg, airs, eos_amsua, hsb, avhrr_navy, avhrr  see read_obs.f90
- id:        satellite id (Tb and pcp obs)
- dload:  computational load factor

- GSI namelist OBS_INPUT

-  &OBS_INPUT
-    dmesh(1)=145.0,dmesh(2)=240.,dmesh(3)=180.,
-    ***dfile(1)***='amsuabufr', dtype(1)='amsua',  id(1)=15, dload(1)=225.0, dval(1)=10.0, dthin(1)=1,
- ...
-    dfile(8)='tmi',       ***dtype(8)***='pcp_tmi',id(8)=211,dload(8)=300.0, dval(8)=1.0,  dthin(8)=-1,
- ...
-    dfile(11)='hirs3bufr',dtype(11)='hirs/3',***id(11)***=16,dload(11)=285.0,dval(11)=6.0, dthin(11)=3,
- ...
-    dfile(21)='prepqc',   dtype(21)='uv',    id(21)=0, ***dload(21)***=3.478,dval(21)=1.0, dthin(21)=0,
-    ...
-  /

## Contents – specific to Tb obs
dthin:  group indicator 
dval:    relative weighting of profile within dthin group
dmesh:  thinning grid mesh size (km) for each dthin group

- dthin:  group indicator
- dval:    relative weighting of profile within dthin group
- dmesh:  thinning grid mesh size (km) for each dthin group

- GSI namelist OBS_INPUT

-  &OBS_INPUT
-    ***dmesh(1)=145.0***,dmesh(2)=240.,**dmesh(3)=180.**,
-    dfile(1)='amsuabufr',dtype(1)='amsua',  id(1)=15, dload(1)=225.0, ***dval(1)=10.0,  dthin(1)=1,***
-    dfile(4)='msubufr',  dtype(4)='msu',    id(4)=14, dload(4)=60.0,  ***dval(4)=2.0,   dthin(4)=1,***
-    dfile(19)='airs', dtype(19)='eos_amsua',id(19)=49,dload(19)=225.0,***dval(19)=10.0, dthin(19)=1***, 				...
-    dfile(3)='hirs2bufr',dtype(3)='hirs/2', id(3)=14, dload(3)=285.0, **dval(3)=6.0,   dthin(3)=3,**
-    dfile(5)='prepqc',   dtype(5)='goes',   id(5)=10, dload(5)=270.0, **dval(5)=6.0,   dthin(5)=3,**
**   **dfile(18)='airs',    dtype(18)='airs',  id(18)=49,dload(18)=4215.,**dval(18)=20.0, dthin(18)=3,**
- /

## To turn on, include oneobtest=.true. in namelist SETUP
Contents
maginnov:  magnitude of O-G
magoberr:  assigned observation error
oneob_type:  observation type
currently must be ps, t, q, u, or v

- Contents
- maginnov:  magnitude of O-G
- magoberr:  assigned observation error
- oneob_type:  observation type
- currently must be ps, t, q, u, or v

- GSI namelist SINGLEOB_TEST

-  &SINGLEOB_TEST
-    ***maginnov***=0.1,***magoberr***=0.1,***oneob***_type='t',
-    oblat=45.,oblon=180.,obpres=1000.,obdattim=2004041512,
-    obhourset=0.,
-  /

## Contents (continued)
oblat, oblon, obpres: observation latitude (°), longitude (°), and pressure (hPa) 
obdattim, obhourset:
Analysis date/time (YYYYMMDDHH) 
Observation time in hours relative to analysis time

- oblat, oblon, obpres: observation latitude (°), longitude (°), and pressure (hPa)
- obdattim, obhourset:
- Analysis date/time (YYYYMMDDHH)
- Observation time in hours relative to analysis time

- GSI namelist SINGLEOB_TEST

-  &SINGLEOB_TEST
-    maginnov=0.1,magoberr=0.1,oneob_type='t',
-    ***oblat***=45.,***oblon***=180.,***obpres***=1000.,***obdattim***=2004041512,
-    ***obhourset***=0.,
-  /

## Observation read/distribute

- gsisub.f90
- Performs one-time only functions
- Code flow
- Various initializations
- Variables & arrays
- Get analysis time, grid information, and type/error/usage information for radiance, precipitation, and ozone assimilation

- User input &
- initializations

- Read in &
- distribute
- observations

## Code level:  gsisub.f90
read_obs.f90
Read ndat observation files listed in OBS_INPUT namelist
Loop over ndat data files
MPI tasks 0 through ndat-1 read data
If ndat > npe, MPI tasks read more than one observation file
Each task reads an observation type / data file
Limited quality control (some gross checks)
Horizontal thinning of brightness temperature data
Each task writes data retained to a unique local file (e.g., obsfil1000, obsfile1001, obsfil1002, ...)
Routines called:  read_airs, read_avhrr_navy, read_bufrtovs, read_goesimg, read_goesndr, read_gps_ref, read_ieeetovs, read_lidar, read_ozone, read_pcp, read_prepbufr, read_radar, read_superwinds

- read_obs.f90
- Read ndat observation files listed in OBS_INPUT namelist
- Loop over ndat data files
- MPI tasks 0 through ndat-1 read data
- If ndat > npe, MPI tasks read more than one observation file
- Each task reads an observation type / data file
- Limited quality control (some gross checks)
- Horizontal thinning of brightness temperature data
- Each task writes data retained to a unique local file (e.g., obsfil1000, obsfile1001, obsfil1002, ...)
- Routines called:  read_airs, read_avhrr_navy, read_bufrtovs, read_goesimg, read_goesndr, read_gps_ref, read_ieeetovs, read_lidar, read_ozone, read_pcp, read_prepbufr, read_radar, read_superwinds

- Read observations

## Code level:  gsisub.f90
deter_subdomain.f90
Divide analysis domain into npe (number of MPI tasks) geographic subdomains
Define subdomains sizes such that computational load is “equally” spread over all tasks

- deter_subdomain.f90
- Divide analysis domain into npe (number of MPI tasks) geographic subdomains
- Define subdomains sizes such that computational load is “equally” spread over all tasks

- Define subdomains

_**dtype**_	_**dload**_
- ps	1.000
- t	2.483
- uv	3.478
- spd	3.478
- q	3.991
- msu	60.0

_**dtype**_	_**dload**_
- ozone	75.0
- amsub	75.0
- amsua	225.0
- hirs	285.0
- pcp	300.0
- airs	4215.0

## Define subdomains: example

- 48 MPI tasks

- Computational load

## Define subdomains: example

- Remove SSM/I and TRMM TMI rain rates

- 48 MPI tasks

- Computational load

## Define subdomains: example

- Set dload = 1.0 for all obs types

- 48 MPI tasks

- Computational load

## Code level:  gsisub.f90
obs_para.f90
Distribute observations to appropriate geographic subdomains
Each MPI tasks reads observation files (obsfil1000, ...)
Observation belonging to geographic subdomain represented by given task is written to task specific observation file  obs_s1000, obs_s1001, ... 
Each observation assigned to a single MPI task   luse
Observations along subdomain edges placed on more than one MPI task  need to do this for int*f90 routines  

- obs_para.f90
- Distribute observations to appropriate geographic subdomains
- Each MPI tasks reads observation files (obsfil1000, ...)
- Observation belonging to geographic subdomain represented by given task is written to task specific observation file  obs_s1000, obs_s1001, ...
- Each observation assigned to a *single* MPI task   luse
- Observations along subdomain edges placed on more than one MPI task  need to do this for int*f90 routines

- Distribute obs to subdomains

## Loops driver/Analysis output

- User input &
- initializations

- Read in &
- distribute
- observations

- Additional
- initializations

- Outer loop
- (Read/distribute) guess & derived fields
- (Read/distribute) background error
- Compute/qc O-G & prepare for inner loop
- Call inner loop
- Compute gradient information
- Apply background error
- Compute search direction
- Compute step size
- Update analysis increment

- Write analysis & related output

## Loops driver/Analysis output

- Code level:  glbsoi.f90
- Determine how many guess files to read (read_*files.f90)
- Allocate/initialize arrays for
- guess fields - create_ges_bias, create_sfc_grids
- minimization - set_pointer, create_jfunc
- background error - create_berror_vars
- radiative transfer model - initialize_rtm
- emissivity model - initialize_irsse_model
*Optionally* read guess solution (iguess)
- read_guess_solution

## Loops driver/Analysis output

- Code level:  glbsoi.f90 (continued)
- “Drive” inner loop miter times
- Write analysis  and other output files:  write_all.f90
- Deallocate arrays:  destroy_sfc_grids, destroy_berror_vars, destroy_ges_bias_grids, destroy_jfunc, destroy_rtm, destroy_irsse_model
- Return execution to calling routines
- gsisub.f90  gsimain.F90

## Cycle through the following 7 steps miter times
Allocate local arrays
create_atm_grids, create_pcp_grids
“Read” guess fields:  read_*guess.f90
1st outer loop only  read guess fields 
All outer loops  compute derived fields
Horizontal derivatives - used in precipitation assimilation
Saturation specific humidity
Geopotential heights
Tropopause pressure

- Allocate local arrays
- create_atm_grids, create_pcp_grids
- “Read” guess fields:  read_*guess.f90
- 1st outer loop only  read guess fields
- All outer loops  compute derived fields
- Horizontal derivatives - used in precipitation assimilation
- Saturation specific humidity
- Geopotential heights
- Tropopause pressure

- Inner loop driver  Outer loop

## Outer loop code flow (continued)
Load arrays containing background error
1st outer loop only
prewgt.f90 or prewgt_reg.f90
Compute/qc O-G & prepare for inner loop (setuprhsall.f90)
Routines involved:  
setup*.f90:  process unique observation types
prep*.f90:    data processing, some quality control
spr*.f90:      gross check, load arrays for inner loop
statsconv.f90, statsoz.f90, statspcp.f90, statsrad.f90:  generate ASCII output of O-G statistics

- Load arrays containing background error
- 1st outer loop only
- prewgt.f90 or prewgt_reg.f90
- Compute/qc O-G & prepare for inner loop (setuprhsall.f90)
- Routines involved:
- setup*.f90:  process unique observation types
- prep*.f90:    data processing, some quality control
- spr*.f90:      gross check, load arrays for inner loop
- statsconv.f90, statsoz.f90, statspcp.f90, statsrad.f90:  generate ASCII output of O-G statistics

- Inner loop driver  Outer loop

## Outer loop code flow (continued)
Deallocate local arrays
destroy_atm_grids, destroy_pcp_grids
Set background error variances for Tb bias correction and ozone
1st outer loop only
set_predictors_vars, set_ozone_var
Call inner (minimization) loop  pcgsoi.f90

- Deallocate local arrays
- destroy_atm_grids, destroy_pcp_grids
- Set background error variances for Tb bias correction and ozone
- 1st outer loop only
- set_predictors_vars, set_ozone_var
- Call inner (minimization) loop  pcgsoi.f90

- Inner loop driver  Outer loop

## Each MPI task
Read observations from task specific file (obs_s1000, obs_s1001, etc.)
Generate simulated observation from guess fields
Apply gross and other quality control decisions
Load arrays containing observation information for use in minimization
Brightness temperature and precipitation information written to scratch files (obs_o1000,obs_o1001, etc.)
Optionally write diagnostic files with innovation information
Local filename OBS_JITER.0000, OBS_JITER.0001, ...
OBS  conv, sbuv2, pcp_tmi, ... etc.
JITER  01, 02, ... 

- Read observations from task specific file (obs_s1000, obs_s1001, etc.)
- Generate simulated observation from guess fields
- Apply gross and other quality control decisions
- Load arrays containing observation information for use in minimization
- Brightness temperature and precipitation information written to scratch files (obs_o1000,obs_o1001, etc.)
*Optionally* write diagnostic files with innovation information
- Local filename* OBS*_*JITER*.0000, *OBS*_*JITER*.0001, ...
*OBS*  conv, sbuv2, pcp_tmi, ... etc.
*JITER*  01, 02, ...

- Step 4 - setuprhsall.f90

## General pcgsoi.f90 flow:
SETUP
Initialize variables and arrays
Read brightness temperature and precipitation observation innovation information from scratch files (obs_o1000, obs_o1001, etc.)
INNER LOOP
Apply preconditioned conjugate gradient descent algorithm to minimize cost function
Loop up to niter times
Exit loop early if 
Magnitude[gradient] reduced more than 10 orders
Total penalty < 1.0e-10

- SETUP
- Initialize variables and arrays
- Read brightness temperature and precipitation observation innovation information from scratch files (obs_o1000, obs_o1001, etc.)
- INNER LOOP
- Apply preconditioned conjugate gradient descent algorithm to minimize cost function
- Loop up to niter times
- Exit loop early if
- Magnitude[gradient] reduced more than 10 orders
- Total penalty < 1.0e-10

- Step 6 - Inner loop (general)

## General pcgsoi.f90 flow:
UPDATE & Exit
Deallocate local arrays 
Update bias correction coefficients
Update guess with current analysis increment

- UPDATE & Exit
- Deallocate local arrays
- Update bias correction coefficients
- Update guess with current analysis increment

- Step 6 - Inner loop (general)

## Cycle through the following 6 steps niter times
Obtain gradient information:  intall.F90 or intall_qc.F90
Separate routine for each observation type:  intt, intps, intq, ...
Apply background error:  bkerror.f90
Compute search direction
Estimate step size:  stpcalc.f90 or stpcalc_qc.f90 
Separate routine for each observation type:  stpt, stpps, stpq, ...
Update solution
Perform diagnostic calculations
Write penalty, gradient, step size information
Check for possible early inner loop exit

- Obtain gradient information:  intall.F90 or intall_qc.F90
- Separate routine for each observation type:  intt, intps, intq, ...
- Apply background error:  bkerror.f90
- Compute search direction
- Estimate step size:  stpcalc.f90 or stpcalc_qc.f90
- Separate routine for each observation type:  stpt, stpps, stpq, ...
- Update solution
- Perform diagnostic calculations
- Write penalty, gradient, step size information
- Check for possible early inner loop exit

- Step 6 – pcgsoi.f90

## Loops driver/Analysis output

- User input &
- initializations

- Read in &
- distribute
- observations

- Additional
- initializations

- Outer loop
- (Read/distribute) guess & derived fields
- (Read/distribute) background error
- Compute/qc O-G & prepare for inner loop
- Call inner loop
- Compute gradient information
- Apply background error
- Compute search direction
- Compute step size
- Update analysis increment

- Write analysis & related output

## Last section of glbsoi.f90
Code flow:
Optionally generate O-A statistics / output files
Write analysis files:  write_all.f90
Atmospheric analysis
Surface analysis (only update Tskin)
Optionally write 
atmospheric bias correction fields
guess (iteration) solution
Deallocate arrays  
destroy_sfc_grids, destroy_berror_vars, destroy_ges_bias_grids, destroy_jfunc, destroy_rtm, destroy_irsse_model

- Code flow:
*Optionally* generate O-A statistics / output files
- Write analysis files:  write_all.f90
- Atmospheric analysis
- Surface analysis (only update Tskin)
*Optionally* write
- atmospheric bias correction fields
- guess (iteration) solution
- Deallocate arrays
- destroy_sfc_grids, destroy_berror_vars, destroy_ges_bias_grids, destroy_jfunc, destroy_rtm, destroy_irsse_model

- Analysis output