da_get_innov_vector_gpspw.inc

References to this file elsewhere.
1 subroutine da_get_innov_vector_gpspw( it, xb, xp, ob, iv)
2 
3    !----------------------------------------------------------------
4    ! Purpose: TBD
5    !----------------------------------------------------------------
6 
7    implicit none
8 
9    integer, intent(in)           :: it      ! External iteration.
10    type(xb_type), intent(in)    :: xb      ! first guess state.
11    type(xpose_type), intent(in) :: xp      ! Domain decomposition vars.
12    type(y_type),  intent(inout) :: ob      ! Observation structure.
13    type(ob_type), intent(inout) :: iv      ! O-B structure.
14 
15    integer                      :: n        ! Loop counter.
16    integer                      :: i, j     ! Index dimension.
17    real                         :: dx, dxm  ! Interpolation weights.
18    real                         :: dy, dym  ! Interpolation weights.
19    real                         :: model_tpw! Model value u at oblocation.
20    integer           :: ittpw,ittpwf
21 
22    !--------------------------------------------------------------------------
23    integer                :: k            ! Index dimension
24    real                   :: dpw, ddpw    ! adjustment pw [mm]
25    real                   :: obs_terr     ! real terrain height at GPS site [m]
26    real                   :: model_terr   ! model terrain height at GPS site[m]
27    real,dimension(xp%kts:xp%kte):: model_q ! model q at GPS site [kg/kg]
28    real,dimension(xp%kts:xp%kte):: model_z ! model z at GPS site [m]
29    real,dimension(xp%kts:xp%kte):: model_rho ! model rho at GPS site [kg/m^3]
30 
31    ! unit_gps = myproc + 140
32    !---------------------------------------------------------------------------
33 
34    !
35    ! GPS TPW Pseudo OBS test:
36    !
37    if (pseudo_var(1:3) == 'tpw' .and. num_pseudo > 0) then
38       ! Deallocate:
39       if (iv%num_gpspw > 0) then
40          write(unit=stdout, fmt='(a,i4)') &
41             'iv%num_gpspw =', iv%num_gpspw
42          deallocate(iv % gpspw)
43          deallocate(ob % gpspw)
44       end if
45 
46       Use_GpspwObs = .true.
47 
48       ! Allocate:
49       iv % num_gpspw  = num_pseudo
50       iv % ob_numb(1)%gpspw = num_pseudo
51       iv % num_pseudo = 0
52 
53       allocate(iv % gpspw(1:num_pseudo))
54       allocate(ob % gpspw(1:num_pseudo))
55 
56       write(unit=stdout,fmt='(a,i2)') &
57         '==> GPS TPW pseudo OBS test: num_pseudo=',num_pseudo
58 
59       iv % gpspw(1) % loc % x   = pseudo_x
60       iv % gpspw(1) % loc % y   = pseudo_y
61       iv % gpspw(1) % loc % i   = int(pseudo_x)
62       iv % gpspw(1) % loc % j   = int(pseudo_y)
63       iv % gpspw(1) % loc % dx  = pseudo_x-real(iv % gpspw(1) % loc % i)
64       iv % gpspw(1) % loc % dy  = pseudo_y-real(iv % gpspw(1) % loc % j)
65       iv % gpspw(1) % loc % dxm = 1.0 - iv % gpspw(1) % loc % dx
66       iv % gpspw(1) % loc % dym = 1.0 - iv % gpspw(1) % loc % dy
67 
68       iv % gpspw(1) % tpw % inv   = pseudo_val
69       iv % gpspw(1) % tpw % qc    = 0
70       iv % gpspw(1) % tpw % error = pseudo_err
71 
72       ! To consider the site elevation effect, set the model terrain height
73       ! to elevation for pseudo OBS:
74 
75       i = iv%gpspw(1)%loc%i
76       j = iv%gpspw(1)%loc%j
77       dx = iv%gpspw(1)%loc%dx
78       dy = iv%gpspw(1)%loc%dy
79       dxm = iv%gpspw(1)%loc%dxm
80       dym = iv%gpspw(1)%loc%dym
81 
82       iv%gpspw(1)%info%elv = dym*(dxm*xb%terr(i,j)   + dx*xb%terr(i+1,j)) + & 
83                              dy *(dxm*xb%terr(i,j+1) + dx*xb%terr(i+1,j+1))
84 
85       ! Set halo:
86       if ((iv%gpspw(1)%loc%i < xp%its-1) .or.(iv%gpspw(1)%loc%i > xp%ite) .or. & 
87           (iv%gpspw(1)%loc%j < xp%jts-1) .or.(iv%gpspw(1)%loc%j > xp%jte)) then 
88          call da_error(__FILE__,__LINE__,(/"Should never have obs outside halo by now"/))
89          ! iv%gpspw(1)%loc%proc_domain_with_halo = .false. 
90          iv%gpspw(1)%loc%proc_domain = .false. 
91       else 
92          ! iv%gpspw(1)%loc%proc_domain_with_halo = .true.  
93          iv%gpspw(1)%loc%proc_domain = .false.  
94 
95          if (iv%gpspw(1)%loc%i >= xp%its .and. iv%gpspw(1)%loc%i <= xp%ite .and. &  
96              iv%gpspw(1)%loc%j >= xp%jts .and. iv%gpspw(1)%loc%j <= xp%jte) then  
97             iv%gpspw(1)%loc%proc_domain = .true.  
98          end if  
99       end if 
100 
101       write(unit=stdout,fmt='(a4,2f15.5)') pseudo_var, pseudo_val, pseudo_err
102       write(unit=stdout,fmt='(3f15.5)')    pseudo_x,  pseudo_y,  pseudo_z
103       write(unit=stdout,fmt='(a,f8.2)')    'iv%gpspw: elv=',iv%gpspw(1)%info%elv
104    end if 
105 
106    if (iv % num_gpspw > 0) then
107 
108       ittpw   = 0
109       ittpwf  = 0
110 
111       ! write(unit=unit_gps,fmt='(3x,a3,10a10)') ' n ','     lat  ',  &
112       !   '     lon  ', '  obs ght ', '  mdl ght ',  &
113       !   ' obsh-mdlh', '   obs pw ', '  model pw',  &
114       !   '   O-B pw ', '    Dpw   ', '  O-B+Dpw '
115 
116       do n=iv%ob_numb(iv%current_ob_time-1)%gpspw + 1, iv%ob_numb(iv%current_ob_time)%gpspw
117 
118          ! [1.1] Get horizontal interpolation weights:
119 
120          i = iv%gpspw(n)%loc%i
121          j = iv%gpspw(n)%loc%j
122          dx = iv%gpspw(n)%loc%dx
123          dy = iv%gpspw(n)%loc%dy
124          dxm = iv%gpspw(n)%loc%dxm
125          dym = iv%gpspw(n)%loc%dym
126 
127          model_tpw = dym*(dxm*xb%tpw(i,j)   + dx*xb%tpw(i+1,j)) + &
128                      dy *(dxm*xb%tpw(i,j+1) + dx*xb%tpw(i+1,j+1))
129 
130          if (pseudo_var(1:3) == 'tpw' .and. num_pseudo > 0) then
131             ! To compute the 'ob':
132             ob % gpspw(n) % tpw = iv % gpspw(n) % tpw % inv + model_tpw
133          else
134             ! To compute the 'inv':
135             iv % gpspw(n) % tpw % inv = 0.0
136             if (ob % gpspw(n) % tpw > missing_r .AND. &
137                 iv % gpspw(n) % tpw % qc >= obs_qc_pointer) then
138                dpw = 0.0
139                obs_terr   = iv%gpspw(n)%info%elv
140                model_terr = dym*(dxm*xb%terr(i,j)   + dx*xb%terr(i+1,j)) + &
141                           dy *(dxm*xb%terr(i,j+1) + dx*xb%terr(i+1,j+1))
142 
143                if (obs_terr <= model_terr) then
144                   model_q(1) = dym*(dxm*xb%q(i,j,1)   + dx*xb%q(i+1,j,1)) + &
145                              dy *(dxm*xb%q(i,j+1,1) + dx*xb%q(i+1,j+1,1))
146                   model_rho(1) = dym*(dxm*xb%rho(i,j,1) + dx*xb%rho(i+1,j,1)) + &
147                              dy *(dxm*xb%rho(i,j+1,1) + dx*xb%rho(i+1,j+1,1))
148                   dpw = model_rho(1) * model_q(1) *( obs_terr - model_terr)
149                else
150                   model_z(1) = dym*(dxm*xb%hf(i,j,1)   + dx*xb%hf(i+1,j,1)) + &
151                             dy *(dxm*xb%hf(i,j+1,1) + dx*xb%hf(i+1,j+1,1))
152                   do k = xp%kts, xp%kte
153                      if (model_z(k) >= obs_terr) exit
154 
155                      model_z(k+1) = dym*(dxm*xb%hf(i,j,k+1)   + dx*xb%hf(i+1,j,k+1)) + &
156                                 dy *(dxm*xb%hf(i,j+1,k+1) + dx*xb%hf(i+1,j+1,k+1))
157                      model_q(k) = dym*(dxm*xb%q(i,j,k)   + dx*xb%q(i+1,j,k)) + &
158                               dy *(dxm*xb%q(i,j+1,k) + dx*xb%q(i+1,j+1,k))
159                      model_rho(k) = dym*(dxm*xb%rho(i,j,k)   + dx*xb%rho(i+1,j,k)) + &
160                                 dy *(dxm*xb%rho(i,j+1,k) + dx*xb%rho(i+1,j+1,k))
161 
162                      if (model_z(k+1) <= obs_terr) then
163                         ddpw = model_rho(k) * model_q(k) *( model_z(k+1) - model_z(k))
164                      else
165                         ddpw = model_rho(k) * model_q(k) *( obs_terr - model_z(k))
166                      end if
167 
168                      dpw = dpw + ddpw
169                   end do
170                end if
171 
172                iv % gpspw(n) % tpw % inv = ob % gpspw(n) % tpw - model_tpw &
173                                          + 0.1*dpw
174             end if
175          end if
176 
177          !------------------------------------------------------------------------
178          ! [5.0] Perform optional maximum error check:
179          !------------------------------------------------------------------------
180 
181          if ( check_max_iv) then
182             call da_check_max_iv_gpspw(it, iv % gpspw(n), ittpw,ittpwf)
183          end if
184 
185          ! write(unit=unit_gps, fmt='(i4,10f10.3)') n, &
186          !    iv%gpspw(n)%info%lat, iv%gpspw(n)%info%lon, obs_terr, &
187          !    model_terr, obs_terr - model_terr, ob%gpspw(n)%tpw,   &
188          !    model_tpw , ob%gpspw(n)%tpw-model_tpw, 0.1*dpw,       &
189          !    ob%gpspw(n)%tpw-model_tpw+0.1*dpw
190       end do
191 
192       if (rootproc .and. check_max_iv_print) then
193          write(unit = check_max_iv_unit, fmt ='(A,i5,A)')'For outer iteration ',it, &
194             ', Total Rejections for Gpspw follows:'
195          write(unit = check_max_iv_unit, fmt = '(/,10(2(A,I6),/))') &
196             'Number of failed TPW-observations:     ',ittpwf, ' on ',ittpw, &
197             'Finally all Gpspw rejections ',ittpwf,' on ',ittpw
198       end if
199    end if
200 
201 end subroutine da_get_innov_vector_gpspw
202 
203