da_transform_xtoy_crtmk_f_adj.inc

References to this file elsewhere.
1 subroutine da_transform_xtoy_crtmk_f_adj ( iv, xp, jo_grad_y, jo_grad_x )
2 
3    !---------------------------------------------------------------------------
4    ! PURPOSE: transform gradient from obs space to model grid space.
5    !
6    ! METHOD:  jo_grad_x = H^T jo_grad_y =  - H^T R^-1 ( d - H delta_x )
7    !           1. input gradient in obs space and reference state of RTTOV
8    !           2. call adjoint of RTM
9    !           3. adjoint of interpolation from model grid to obs loc
10    !
11    !  HISTORY: 11/16/2006 - Creation            Zhiquan Liu
12    !
13    !---------------------------------------------------------------------------
14 
15    IMPLICIT NONE
16 
17    TYPE (x_type), INTENT(INOUT)   :: jo_grad_x ! 
18    TYPE (xpose_type), INTENT(IN)  :: xp        ! Domain decomposition vars.
19    TYPE (y_type),  INTENT(IN)     :: jo_grad_y ! H' delta_x
20    TYPE (ob_type), INTENT(IN)     :: iv        ! O-B structure.
21 
22    INTEGER                        :: l, i, j, k  ! Index dimension.
23    INTEGER                        :: num_rad  ! Number of radiance obs
24    REAL                           :: dx, dxm  ! Interpolation weights.
25    REAL                           :: dy, dym  ! Interpolation weights.
26    INTEGER                        :: inst, n
27    REAL, allocatable              :: q_ad(:),t_ad(:)
28    REAL                           :: ps_ad
29 
30    ! CRTM local varaibles and types
31    INTEGER :: Allocate_Status
32 !---------------------------------------------------------
33 
34 #if !defined(CRTM)
35     call da_error(__FILE__,__LINE__, &
36        (/"Must compile with $CRTM option for radiances"/))
37 #else
38    IF ( iv%num_inst < 1 ) return
39 
40    if (trace_use) call da_trace_entry("da_transform_xtoy_crtmk_f_adj")
41 
42 !-------------------------------------------------------------------------------
43 
44    do inst = 1, iv%num_inst                 ! loop for sensor
45       if ( iv%instid(inst)%num_rad < 1 ) cycle
46       num_rad = iv%instid(inst)%num_rad
47 
48   ! Allocate forward model solution RTSolution array to number of channels
49       ALLOCATE (t_ad(xp%kte-xp%kts+1))
50       ALLOCATE (q_ad(xp%kte-xp%kts+1))
51 
52 !----------------------------------------------------------------------------
53       do n=1,num_rad
54 
55       ! [1.5] Scale transformation and fill zero for no-control variable
56           ps_ad = 0.
57          do l=1, iv%instid(inst)%nchan
58           ps_ad = ps_ad + iv%instid(inst)%ps_jacobian(l,n)*jo_grad_y%instid(inst)%tb(l,n)
59          end do
60           ps_ad = ps_ad*0.01
61 
62          do k=xp%kts,xp%kte
63           t_ad(k) = 0.
64           q_ad(k) = 0.
65          do l=1, iv%instid(inst)%nchan
66           t_ad(k)  = t_ad(k) + iv%instid(inst)%t_jacobian(l,k,n)*jo_grad_y%instid(inst)%tb(l,n)
67           q_ad(k)  = q_ad(k) + iv%instid(inst)%q_jacobian(l,k,n)*jo_grad_y%instid(inst)%tb(l,n)
68          end do
69           q_ad(k) = 1000.*q_ad(k)
70          end do
71           
72             !-----------------------------------------------------
73             ! [1.6] Get horizontal interpolation weights:
74             !-----------------------------------------------------
75 
76             i = iv%instid(inst)%loc_i(n)
77             j = iv%instid(inst)%loc_j(n)
78             dx = iv%instid(inst)%loc_dx(n)
79             dy = iv%instid(inst)%loc_dy(n)
80             dxm = iv%instid(inst)%loc_dxm(n)
81             dym = iv%instid(inst)%loc_dym(n)
82 
83             ! [1.7] Adjoint of Interpolate horizontally from ob to grid:
84 
85             do k=xp%kts,xp%kte ! from bottem to top
86                if (iv%instid(inst)%pm(xp%kte-k+1,n) < 75.) q_ad(xp%kte-k+1) = 0.
87                call da_interp_lin_2d_adj( jo_grad_x%t(:,:,k), xp%ims, xp%ime, xp%jms, &
88                                           xp%jme, i, j, dx, dy, dxm, dym, &
89                                           t_ad(xp%kte-k+1))
90 
91                call da_interp_lin_2d_adj( jo_grad_x%q(:,:,k), xp%ims, xp%ime, xp%jms, & 
92                                            xp%jme, i, j, dx, dy, dxm, dym, &
93                                            q_ad(xp%kte-k+1) )
94             enddo
95 
96                call da_interp_lin_2d_adj(jo_grad_x% psfc, xp%ims, xp%ime, xp%jms, &
97                                          xp%jme, i, j, dx, dy, dxm, dym,  &
98                                          ps_ad )
99 
100                !if (n <=10 ) then
101                !   write(6,'(e15.5)') jo_grad_x% psfc(i,j)
102                !  do k=xp%kts,xp%kte
103                !   write(6,'(2e15.5)') jo_grad_x%t(i,j,k), jo_grad_x%q(i,j,k)
104                !  enddo
105                !endif
106 
107          end do       !  end loop for pixels
108 
109       !-------------------------------------------------------------------
110       ! [2.0] Deallocating CRTM structures
111       !-------------------------------------------------------------------
112 
113          deallocate( t_ad, q_ad, &
114                      STAT = Allocate_Status )
115          if ( Allocate_Status /= 0 ) THEN
116               call da_error(__FILE__,__LINE__, &
117                   (/"Error in deallocatting t_ad q_ad"/))
118          endif
119 
120    end do        ! end loop for sensor
121 
122    if (trace_use) call da_trace_exit("da_transform_xtoy_crtmk_f_adj")
123 #endif
124  
125 end subroutine da_transform_xtoy_crtmk_f_adj
126