Docstoc

4DVAR

Document Sample
4DVAR Powered By Docstoc
					                          M o d e li n                                                           -F o llo w i
            al   Ocean                   g S                                             r rai n              ng
                                                                                    Te
     i   on                                  ys
                                                t
 g




                                                em
Re




                                                                                                                   m
                                                                                   ce




                                                                               O
                                                                                        an                         te
            R es
                   e a rc h C o mm u ni
                                           ty                                                M o d e lin g S y s




                                                4D Variational Data Assimilation
                                                    Observation Operators



                                                       Hernan G. Arango
                 ROMS 4DVAR ALGORITHMS

•   Strong Constraint
      Conventional (S4DVAR): outer loop, NLM, ADM
      Incremental (IS4DVAR): inner and outer loops, NLM, TLM, ADM (Courtier et al.,
        1994)
           IS4DVAR_OLD: inefficient old conjugate gradient algorithm
            (is4dvar_ocean_old.h, descent.F)
           IS4DVAR: new conjugate gradient algorithm, CONGRAD, Fisher 1997
            (is4dvar_ocean.h, cgradient.h)
           IS4DVAR, LANCZOS: conjugate gradient and Lanczos algorithm, Fisher
            1997 (is4dvar_ocean_lanczos.h, cgradient_lanczos.h)
•   Weak Constraint
      Indirect Representer Method (W4DVAR): inner and outer loops, NLM, TLM,
        RPM, ADM (Egbert et al., 1994; Bennett et al, 1997)
      Physical Space Statistical Analysis (W4DPSAS): inner and outer loops, NLM,
        TLM, ADM (Courtier, 1997)
              Strong Constraint, Incremental 4DVAR

Let’s introduce a new minimization variable v, such that:
   xk = B-1/2(xk + xk-1 – xb)
   xk = B1/2vk + xk-1 – xb
yielding

   J(vk) = ½(vk)Tvk + ½(Hxk – dk-1)TO-1(Hxk – dk-1)

The gradient of J in minimization-space, denoted v J, is given by:

   v J = vk + BT/2HTO-1(Hxk – dk-1) = vk + BT/2x Jo => vk + W -1/2LT/2GS

The background-error covariance matrix can be factored as:

   B = SCS => S(GL1/2W -1/2)(W -1/2LT/2G)S

where S is the background-error standard deviations, C is the
background-error correlations which can be factorized as C = C1/2CT/2,
G is the normalization matrix which ensures that the diagonal
elements of C are equal to unity, L is a 3D self-adjoint filtering
operator, and W is the grid cell area or volume.
                                  Basic IS4DVAR Procedure

             (1)           Choose an x(0) = xb(0)

             (2)           Integrate NLROMS t  [0, ] and save x(t)                     (NLM at OBS)

                           (a) Choose a x(0)

                           (b) Integrate TLROMS t  [0, ] and compute J                 (TLM at OBS)

                           (c) Integrate ADROMS t  [0, ] to yield                      (ADM forcing at OBS)
Outer Loop

              Inner Loop




                           (d) Compute

                           (e) Use a descent algorithm to determine a “down gradient”
                              correction to x(0) that will yield a smaller value of J

                           (f) Back to (b) until converged

             (3)           Compute new x(0) = x(0) + x(0) and back to (2) until converged
        Incremental, Strong Constraint 4DVar
                     (IS4DVAR)




• Given a first guess (a forward trajectory)…
• And given the available data…
       Incremental, Strong Constraint 4DVar
                    (IS4DVAR)




• Given a first guess (a forward trajectory)…
• And given the available data…
• IS4DVAR computes the changes (or increments) to
 the initial conditions so that the forward model fits
 the observations.
                 4DVAR Observations
                    NetCDF File




•   Utility/obs_initial.F
•   Utility/obs_read.F
•   Utility/obs_write.F
•   Utility/obs_scale.F
•   Utility/obs_depth.F
•   Utility/extract_obs.F
•   Adjoint/ad_extract_obs.F
•   Adjoint/ad_misfit.F
                      Metadata

Dimensions:

survey                Number of different time
weight                Number of interpolation weight
datum                 Observations counter, unlimited dimension

Variables:

Nobs(survey)          Number of observations per time survey
survey_time(survey)   Survey time (days)
obs_type(datum)       State variable ID associated with observation
obs_time(datum)       Time of observation (days)
obs_lon(datum)        Longitude of observation (degrees_east)
obs_lat(datum)        Latitude of observation (degrees_north)
obs_depth(datum)      Depth of observation (meters or level)
obs_Xgrid(datum)      X-grid observation location (nondimensional)
obs_Ygrid(datum)      Y-grid observation location (nondimensional)
obs_Zgrid(datum)      Z-grid observation location (nondimensional)
obs_error(datum)      Observation error, assigned weight
obs_value(datum)      Observation value
dimensions:
                                                                       Observations NetCDF
        survey = 1 ;
        tate_variable = 7 ;
        datum = UNLIMITED ; // (79416 currently)
                                                                                                        Variable   ID
                                                                                                            ζ
variables:
                                                                                                                    1
        char spherical ;
                spherical:long_name = "grid type logical switch" ;                                         u        2
        int Nobs(survey) ;                                                                                 v        3
                Nobs:long_name = "number of observations with the same survey time" ;
        double survey_time(survey) ;                                                                       u        4
                survey_time:long_name = "survey time" ;
                survey_time:units = "days since 2000-01-01 00:00:00" ;                                     v        5
                survey_time:calendar = "365.25 days in every year" ;
        double obs_variance(state_variable) ;                                                            temp       6
                obs_variance:long_name = "global (time and space) observation variance" ;                 salt      7
                obs_variance:units = "squared state variable units" ;
        int obs_type(datum) ;
                obs_type:long_name = "model state variable associated with observation" ;
                obs_type:units = "nondimensional" ;
        double obs_time(datum) ;
                obs_time:long_name = "time of observation" ;
                obs_time:units = "days since 2000-01-01 00:00:00" ;
                obs_time:calendar = "365.25 days in every year" ;
        double obs_depth(datum) ;
                obs_depth:long_name = "depth of observation" ;
                obs_depth:units = "meter" ;
        double obs_Xgrid(datum) ;
                obs_Xgrid:long_name = "x-grid observation location" ;
                obs_Xgrid:units = "nondimensional" ;                                                        8           (i2,j2,k2)
        double obs_Ygrid(datum) ;
                                                                                                                   7
                                                                                                    5
                obs_Ygrid:long_name = "y-grid observation location" ;
                obs_Ygrid:units = "nondimensional" ;
        double obs_Zgrid(datum) ;
                obs_Zgrid:long_name = "z-grid observation location" ;
                                                                                                        6
                                                                                                            4
                obs_Zgrid:units = "nondimensional" ;
        double obs_error(datum) ;
                obs_error:long_name = "observation error, assigned weight, inverse variance" ;
                obs_error:units = "inverse squared state variable units" ;
                                                                                                                   3
        double obs_value(datum) ;
                obs_value:long_name = "observation value" ;                            (i1,j1,k1)
                                                                                                    1
                obs_value:units = "state variable units" ;
                                                                                                        2
                   Processing
• Use hindices, try_range and inside routines to
  transform (lon,lat) to (,)
• Find how many survey times occur within the data
  set (survey dimension)
• Count observations available per survey (Nobs) and
  assign their times (survey_time)
• Sort the observation in ascending time order and
  observation time for efficiency
• Save a copy of the observation file
• Several matlab scripts to process observations
              ROMS GRID



• Horizontal curvilinear orthogonal coordinates
  on an Arakawa C-grid
• Terrain-following coordinates on a staggered
  vertical grid
Curvilinear Transformation
Staggered C-Grid, RHO-points
Staggered C-Grid, U-points
Staggered C-Grid, V-points
      Vertical Terrain-following Coordinates
Vieste                   Dubrovnik
(Italy)                  (Croatia)




Depth
 (m)




            Longitude
Parabolic Splines Reconstruction
      PARALLEL TILE PARTITIONS




8x8




                    Ny
                         } Nx
            East-West MPI Communications
                                          AsendE    ArecvE                 ArecvW    AsendW
                      TILE L                                                                                 TILE R

                                                               receive
               Jend                                                                                                   Jend

Nonlinear                                                                                                                     With Respect
                                                                                                                              To Tile R
                                         i-2 i-1    i   i+1                i-2 i-1   i     i+1
               Jstr                                                                                                   Jstr

                                                                 send

                            Istr             Iend                                   Istr                 Iend
                                 L         R                                               R           L
                               Vi =      Vi                                              V i-2 =     V i-2
                                 L         R                                               R           L
                               V i+1 =   V i+1                                           V i-1 =     V i-1

                      TILE L              AsendE    ArecvE                 ArecvW    AsendW
                                                                                                             TILE R

                                                               ad_send

               Jend                                                                                                   Jend


Adjoint                                                                                                                       With Respect
                                                                                                                              To Tile R
                                         i-2 i-1    i   i+1                i-2 i-1   i     i+1
               Jstr                                                                                                   Jstr

                                                              ad_receive

                            Istr           Iend                                     Istr                Iend
                 L            L          R         R                     R              R                 L            L
            ad_V i-2 =   ad_V i-2 + ad_V i-2 ; ad_Vi-2 = 0          ad_V i     =   ad_V i        +        i     ; ad_V i     =0
                 L          L          R         R                       R          R          L          L
            ad_V i-1 = ad_V i-1 + ad_V i-1 ; ad_Vi-1 = 0            ad_V i+1 = ad_V i+1 + ad_V i+1 ; ad_V i+1 = 0
                     North-South MPI Communications
                                 B       T                                      T          T          B          B
                               V j+2 = V j+2                               ad_V j+2 = ad_V j+2 + ad_V j+2 ; ad_V j+2 = 0
                                 B       T                                      T          T          B          B
                               V j+1 = V j+1                               ad_V j+1 = ad_V j+1 + ad_V j+1 ; ad_V j+1 = 0

                                 TILE T                                                            TILE T
                                                                With Respect
                                                                  to Tile B
            Jend                                                                Jend




                          j+2                                                                j+2

                          j+1                         AsendS                                 j+1                       AsendS
            Jstr                                                                Jstr
                           j                                                                  j
Nonlinear




                          j-1                          ArecvS                                j-1                       ArecvS




                                                                                                                                Adjoint
                          Istr            Iend                                              Istr        Iend
                   send                          receive                       ad_receive                      ad_send
                                 TILE B                                                            TILE B
                          j+2                                                                j+2

                          j+1                          ArecvN                                j+1                       ArecvN

                           j                                                                  j
            Jend                                                                Jend
                          j-1                          AsendN                                j-1                       AsendN




            Jstr                                                                Jstr




                          Istr            Iend                                              Istr        Iend
                                 T      B                                       B             B         T          T
                                Vj   = Vj                                  ad_V j      = ad_V j     +   j   ; ad_V j     =0
                                  T       B                                     B          B          T          T
                                V j-1 = V j-1                              ad_V j-1 = ad_V j-1 + ad_V j-1 ; ad_V j-1 = 0
Observations
                                     THE STATE UNIVERSITY OF NEW JERSEY


                                     RUTGERS
                                                                                        LEO-15                                 LEO              NJSOS
LEO-15 Research Stations                                                                                         Kilometers
                                                                                                     Beach
  M                                                                                                  Haven
      ul                                                                                                         0    1   2
           lic
                 a
                     Ri
                          ve
                               r
                                                                                            Long Beach
                                   Great
                                                         Field
                                                        Station
                                                       Field Station
                                    Bay                                                                              Longterm Ecosystem Observatory
                                                                       Little Egg
                                                                          Inlet
                                                                          C
                                                                              ab
                                                                                   le
 Depth
                                                                                                             LEO-15
 Land and
 Wetlands
                                                                                        Node A
                                                                                                 Node B      Research
 1 meter                                                   Brigantine                                         Area
 3 meters                                                     Inlet
 6 meters
 10 meters
 14 meters
 16 meters
 18 meters
 22 meters
                                        Brigantine                                                                                          New Jersey Shelf
                                                                                                                                            Observing
                                                                                                                                            System
                                                   3km x 3km
                                                  1996-Present


                                                                                                                              30km x 30km
                                                                                                                               1998-2001

                                                                                                                                                                              Satellites,
                                                                                                                                                                       Aircraft, Surface
                                                                                                                                                                    RADAR, Glider AUVs
                                                                                                                                                         300km x 300km
                                                                                                                                                         Beginning 2001
                   Assumptions
• All scalar observations are assumed to be at
 RHO-points.
• All vector observations are assumed to be rotated
 to curvilinear grid, if applicable. Vector
 observations are always measured at the same
 location.
• All observation horizontal locations are assumed
 to be in fractional curvilinear grid coordinates.
• Vertical locations can be in fractional levels (1:N)
 or actual depths (negative values).
• Removal of tidal signal?
• Filtering of non-resolved processes?
              Observation Operators

• Currently, all observations must be in terms of
 model state variables (same units):
  – 2D configuration: zeta, ubar, vbar
  – 3D configuration: zeta, u, v, T, S, …
• There is not a time interpolation of model solution
 at observation times:
       time - 0.5*dt < ObsTime < time + 0.5*dt
• There is not observations quality control
 (screening) inside ROMS, except ObsScale.
• No observation constraints are implemented
 (Satellite SST measurements)
           Observation Interpolation

• Only spatial linear interpolation is coded.
• If land/sea masking, the interpolation
 coefficients are weighted by the mask.
• If shallower than z_r(:,:,N), observations are
 assigned to the surface level.
• If deeper than z_r(:,:,1), observations are assigned
 to bottom level.
               Recommedations

• Create a NetCDF file for each observation type.
• Use a processing program to meld NetCDF
 observation files (nc_4dvar_meld.m).
• Keep a master copy of each observation file, in
 case that you are running your application at
 different resolutions.
• Decimation of observations. Finite volume
 representation?
   BACKGOUND
ERROR COVARIANCE
    Model/Background Error Covariance, B
• Use a generalized diffusion squared-root operator
  (symmetric) as in Weaver et al. (2003):
                                 1/2       -1/2        -1/2   T/2
              B = S C S = S (G L       W          ) (W        L     G)
• The normalization matrix, G, ensure that the diagonal
  elements of the correlation matrix, C, are equal to unity.
  They are computed using the exact (expensive) or
  randomization (cheaper) methods.
• The spatial convolution of the self-adjoint filtering operator,
   1/2
  L , is split in horizontal and vertical components and
  discretized both explicitly and implicitly.
• The model/background standard deviation matrix, S, is
  computed from long (monthly, seasonal) simulations.
                                                    -1/2
• The grid cell area or volume matrix, W                   , is assumed to be
  time invariant.
      Model/Background Error Correlation (C)

  Horizontal                 Vertical (implicit)




                              Vdecay = 100 m



Hdecay = 100 km
      Model/Background Error Correlation Normalization
                     Coefficients (G)

            SSH                        Temperature

EAC                             EAC




                                       Bottom Level
                             job_is4dvar.sh
#! / bi n/ csh - f

          ent
# I ncr em al S4DVAR j ob scr i pt :
#
# Thi s scr i pt NEEDS t o be r un bef or e any r un:
#
#    ( 1) I t copi es a new cl ean nonl i near m    odel i ni t i al condi t i ons
#         f i l e. The nonl i near model i s i ni t i al i zed f r om t he
#         backgr ound or r ef er ence st at e.
#    ( 2) Speci f y t he backgr ound- er r or st andar d devi at i on,
#                al
          nor m i zat i on, obser vat i on, Hessi an vect or s, and out put
#         4DVAR f i l es.
#                                                     pl
     ( 3) Cr eat e 4DVAR i nput scr i pt f r om t em at e " s4dvar . i n" and
#         r epl ace t he val ues of I O f i l es wi t h t hose speci f i ed i n ( 2) .
#

                e/                   s/
 set Di r =/ hom ar ango/ ocean/ t om r eposi t or y/ Pr oj ect s/ doubl e_gyr e

# Copy nonl i near model i ni t i al condi t i ons f i l e.

 cp - p ${ Di r } / Dat a/ gyr e3d_bck. nc gyr e3d_i ni . nc               # 10 year mean

# Set backgr ound- er r or covar i ance st andar d devi at i ons f i l e.
           e=.
 set STDnam . / Dat a/ gyr e3d_st d. nc

                                             al
# Set backgr ound- er r or covar i ance nor m i zat i on f act or f i l e

        nam
 set NRM e=. . / Dat a/ gyr e3d_nr m_100k100i . nc

# Set obser vat i ons f i l e.
#set OBSname=gyr e3d_zSST_obs. nc
 set OBSname=gyr e3d_zuvTS_obs. nc

    et
# G a cl ean copy of t he obser vat i on f i l e.         Thi s i s r eal l y
     por                               odi
# i m t ant si nce t hi s f i l e i s m f i ed.

                     BS/ BSnam .
 cp - p ${ Di r } / O $O      e
   odi             pl
# M f y S4DVAR t em at e i nput scr i pt and speci f y above f i l es.
                                                            build.sh
#! / bi n/ csh - f
#
# svn $I d: bui l d. sh 677 2008- 08- 05 20: 17: 30Z ar ango $
                                                                                                                                   i
#: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : John W l ki n : : :
# Copyr i ght ( c) 2002- 2008 The RO S/ TO S G oup             M         M r                                                                   :::
#                                          I
       Li censed under a M T/ X st yl e l i cense                                                                                              :::
#      See Li cense_RO S. t xt     M                                                                                                           :::
                                                                                                                             .
#: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : Her nan G Ar ango : : :
#                                                                                                                                              :::
       M         M           pi
# RO S/ TO S Com l i ng Scr i pt                                                                                                               :::
#                                                                                                                                              :::
                             pi
# Scr i pt t o com l e an user appl i cat i on wher e t he appl i cat i on- speci f i c : : :
# f i l es ar e kept separ at e f r om t he RO S sour ce code.           M                                                                     :::
#                                                                                                                                              :::
     :
# Q How/ why does t hi s scr i pt wor k?                                                                                                       :::
#                                                                                                                                              :::
# A: The RO S m      M akef i l e conf i gur es user - def i ned opt i ons wi t h a set of                                                     :::
#                                        M
         f l ags such as RO S_APPLI CATI O Br owse t he m          N.                            akef i l e t o see t hese. : : :
#        I f an opt i on i n t he m                akef i l e uses t he synt ax ?= i n set t i ng t he                                         :::
#        def aul t , t hi s m          eans t hat m          ake wi l l check whet her an envi r onm                             ent           :::
#                                                e
         var i abl e by t hat nam i s set i n t he shel l t hat cal l s m                                            ake. I f so               :::
#        t he envi r onm         ent var i abl e val ue over r i des t he def aul t ( and t he                                                 :::
#                                      ai
         user need not m nt ai n separ at e m                              akef i l es, or f r equent l y edi t                                :::
#        t he m    akef i l e, t o r un separ at e appl i cat i ons) .                                                                         :::
#                                                                                                                                              :::
# Usage:                                                                                                                                       :::
#                                                                                                                                              :::
#        . / bui l d. sh [ opt i ons]                                                                                                          :::
#                                                                                                                                              :::
     pt
# O i ons:                                                                                                                                     :::
#                                                                                                                                              :::
#        - j [ N]                      pi
                                 Com l e i n par al l el usi ng N CPUs                                                                         :::
#                                        i
                                     om t ar gum         ent f or al l avai l abl e CPUs                                                       :::
#        - nocl ean              Do not cl ean al r eady com l ed obj ect s      pi                                                            :::
#                                                                                                                                              :::
                                 et es
# Not i ce t hat som i m t he par al l el com l at i on f ai l t o f i nd M    pi                                            PI                :::
# i ncl ude f i l e " m f . h" . pi                                                                                                            :::
#                                                                                                                                              :::
#: : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : : :
set par al l el = 0
                                           s4dvar.in
                   i
! 4DVar assi m l at i on i nput par am er s.et
!
! svn $I d: s4dvar . i n 503 2008- 01- 10 00: 11: 51Z ar ango $
! ========================================================= Her nan G Ar ango ===      .
                                            M
! Copyr i ght ( c) 2002- 2008 The RO S/ TO S G oup M r                                            !
!                              I
        Li censed under a M T/ X st yl e l i cense                                                !
!                        M
        See Li cense_RO S. t xt                                                                   !
! ==============================================================================
!                                                                                                 !
                  et
! I nput par am er s can be ent er ed i n ANY or der , pr ovi ded t hat t he par am er     et     !
         O
! KEYW RD ( usual l y, upper case) i s t yped cor r ect l y f ol l owed by " =" or " =="          !
        bol             m
! sym s. Any com ent l i nes ar e al l owed and m                                         at
                                                             ust begi n wi t h an excl am i on    !
    ar                   n
! m k ( ! ) i n col um one. Com ent s m m         ay appear t o t he r i ght of a par am eret     !
                            pr            ent                m
! speci f i cat i on t o i m ove docum at i on. Com ent s wi l l be i gnor ed dur i ng            !
! r eadi ng. Bl ank l i nes ar e al so al l owed and i gnor ed. Cont i nuat i on l i nes i n      !
              et
! a par am er speci f i cat i on ar e al l owed and m       ust be pr eceded by a backsl ash      !
                   e                 or
! ( \ ) . I n som i nst ances, m e t han one val ue i s r equi r ed f or a par am er .    et      !
! I f f ewer val ues ar e pr ovi ded, t he l ast val ue i s assi gned f or t he ent i r e         !
           et                    ul
! par am er ar r ay. The m t i pl i cat i on sym         bol ( * ) , wi t hout bl ank spaces i n  !
                                            et                                    pl
! bet ween, i s al l owed f or a par am er speci f i cat i on. For exam e, i n a t wo             !
! gr i ds nest ed appl i cat i on:                                                                !
!                                                                                                 !
!       AKT_BAK == 2* 1. 0d- 6 2* 5. 0d- 6                           2/
                                                                 ! m s                            !
!                                                                                                 !
! i ndi cat es t hat t he f i r st t wo ent r i es of ar r ay AKT_BAK, i n f or t r an col um  n- !
    aj                                    e
! m or or der , wi l l have t he sam val ue of " 1. 0d- 6" f or gr i d 1, wher eas t he !
                                                e
! next t wo ent r i es wi l l have t he sam val ue of " 5. 0d- 6" f or gr i d 2.                  !
!                                                                                                 !
         ul                                           ul                         ai
! I n m t i pl e l evel s of nest i ng and/ or m t i pl e connect ed dom ns st ep- ups, !
                                                    e                     et
! " Ngr i ds" ent r i es ar e expect ed f or som of t hese par am er s. I n such case, !
                                                   et                 el     por
! t he or der of t he ent r i es f or a par am er i s ext r em y i m t ant . I t m            ust !
                      e
! f ol l ow t he sam or der ( 1: Ngr i ds) as i n t he st at e var i abl e decl ar at i on. The !
            ay
! USER m f ol l ow t he above gui del i nes f or speci f yi ng hi s/ her val ues. These !
           et           ar
! par am er s ar e m ked by " ==" pl ur al sym                                 O
                                                         bol af t er t he KEYW RD.                !
!                                                                                                 !
! ==============================================================================
!

! Swi t ch t o r ecycl e r ecor ds i n i ni t i al f i l e, [ 1: Ngr i ds] .
      Lcycl eI NI == T