da_transfer_xatowrftl_adj.inc
References to this file elsewhere.
1 subroutine da_transfer_xatowrftl_adj(grid, config_flags, filnam)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Convert WRFTL variables to analysis increments
5 ! (inverse of the incremental part of xatowrf)
6 !---------------------------------------------------------------------------
7
8 implicit none
9
10 type(domain), intent(inout) :: grid
11 type(grid_config_rec_type), intent(inout) :: config_flags
12
13 character*4, intent(in) :: filnam
14
15 ! Local variables
16
17 integer :: i, j, k
18 real :: sdmd, s1md
19 real :: rho_cgrid
20
21 real, dimension(ims:ime,jms:jme,kms:kme) :: a_press
22 real, dimension(ims:ime,jms:jme,kms:kme) :: utmp
23 real, dimension(ims:ime,jms:jme,kms:kme) :: vtmp
24
25
26 integer ndynopt
27
28 if (trace_use) call da_trace_entry("da_transfer_xatowrftl_adj")
29
30 !---------------------------------------------------------------------------
31 ! [7.0] Adjoint of outPUT (inPUT)
32 !---------------------------------------------------------------------------
33
34 ndynopt = grid%dyn_opt
35 grid%dyn_opt = DYN_EM_AD
36 call nl_set_dyn_opt (1 , DYN_EM_AD)
37
38 call da_med_initialdata_input(grid , config_flags, filnam)
39
40 grid%dyn_opt = ndynopt
41 call nl_set_dyn_opt (1 , DYN_EM)
42
43 !---------------------------------------------------------------------------
44 ! [6.0] Adjoint of save OTHERinCREMENT
45 !---------------------------------------------------------------------------
46
47 do k=kts,kte+1
48 do j=jts,jte
49 do i=its,ite
50 grid%xa%w(i,j,k)=grid%xa%w(i,j,k)+grid%a_w_2(i,j,k)
51 end do
52 end do
53 end do
54
55 grid%a_w_2 = 0.0
56
57 #ifdef VAR4D_MICROPHYSICS
58 ! New code needed once we introduce the microphysics to 4dvar in 2008
59 if (size(grid%moist,dim=4) >= 4) then
60 do k=kts,kte
61 do j=jts,jte
62 do i=its,ite
63 grid%xa%qcw(i,j,k)=grid%a_moist(i,j,k,p_qcw)
64 grid%xa%qrn(i,j,k)=grid%a_moist(i,j,k,p_qrn)
65 end do
66 end do
67 end do
68 end if
69
70 if (size(grid%moist,dim=4) >= 6) then
71 do k=kts,kte
72 do j=jts,jte
73 do i=its,ite
74 grid%xa%qci(i,j,k)=grid%a_moist(i,j,k,p_qci)
75 grid%xa%qsn(i,j,k)=grid%a_moist(i,j,k,p_qsn)
76 end do
77 end do
78 end do
79 end if
80
81 if (size(grid%moist,dim=4) >= 7) then
82 do k=kts,kte
83 do j=jts,jte
84 do i=its,ite
85 grid%xa%qgr(i,j,k)=grid%a_moist(i,j,k,p_qgr)
86 end do
87 end do
88 end do
89 end if
90
91 #endif
92 !---------------------------------------------------------------------------
93 ! [5.0] Adjoint of CONVERT FROM A-GRID TO C-GRID
94 !---------------------------------------------------------------------------
95
96 ! Fill the halo region for a_u and a_v.
97 utmp=grid%xa%u
98 vtmp=grid%xa%v
99 grid%xa%u=grid%a_u_2
100 grid%xa%v=grid%a_v_2
101 #ifdef DM_PARALLEL
102 #include "HALO_PSICHI_UV_ADJ.inc"
103 #endif
104 grid%a_u_2=grid%xa%u
105 grid%a_v_2=grid%xa%v
106 grid%xa%u=utmp
107 grid%xa%v=vtmp
108 utmp=0.0
109 vtmp=0.0
110
111 do k=kts,kte
112 do j=jts,jte
113 do i=its,ite
114 utmp(i,j,k)=utmp(i,j,k)+0.5*(grid%a_u_2(i+1,j ,k)+grid%a_u_2(i,j,k))
115 vtmp(i,j,k)=vtmp(i,j,k)+0.5*(grid%a_v_2(i ,j+1,k)+grid%a_v_2(i,j,k))
116 end do
117 end do
118 end do
119
120 utmp(its-1,jts:jte,kts:kte)=utmp(its-1,jts:jte,kts:kte)+0.5*grid%a_u_2(its,jts:jte,kts:kte)
121 utmp(ite+1,jts:jte,kts:kte)=utmp(ite+1,jts:jte,kts:kte)+0.5*grid%a_u_2(ite+1,jts:jte,kts:kte)
122 vtmp(its:ite,jts-1,kts:kte)=vtmp(its:ite,jts-1,kts:kte)+0.5*grid%a_v_2(its:ite,jts,kts:kte)
123 vtmp(its:ite,jte+1,kts:kte)=vtmp(its:ite,jte+1,kts:kte)+0.5*grid%a_v_2(its:ite,jte+1,kts:kte)
124
125 ! The western boundary
126 if (its == grid%xp%ids) then
127 grid%xa%u(its ,jts:jte,kts:kte)=grid%xa%u(its ,jts:jte,kts:kte)+2.0*utmp(its-1,jts:jte,kts:kte)
128 grid%xa%u(its+1,jts:jte,kts:kte)=grid%xa%u(its+1,jts:jte,kts:kte)-utmp(its-1,jts:jte,kts:kte)
129 end if
130
131 ! The eastern boundary
132 if (ite == grid%xp%ide) then
133 grid%xa%u(ite ,jts:jte,kts:kte)=grid%xa%u(ite ,jts:jte,kts:kte)+2.0*utmp(ite+1,jts:jte,kts:kte)
134 grid%xa%u(ite-1,jts:jte,kts:kte)=grid%xa%u(ite-1,jts:jte,kts:kte)-utmp(ite+1,jts:jte,kts:kte)
135 end if
136
137 grid%xa%u=grid%xa%u+utmp
138
139 ! The southern boundary
140 if (jts == grid%xp%jds) then
141 grid%xa%v(its:ite,jts ,kts:kte)=grid%xa%v(its:ite,jts ,kts:kte)+2.0*vtmp(its:ite,jts-1,kts:kte)
142 grid%xa%v(its:ite,jts+1,kts:kte)=grid%xa%v(its:ite,jts+1,kts:kte)-vtmp(its:ite,jts-1,kts:kte)
143 end if
144
145 ! The northern boundary
146 if (jte == grid%xp%jde) then
147 grid%xa%v(its:ite,jte ,kts:kte)=grid%xa%v(its:ite,jte ,kts:kte)+2.0*vtmp(its:ite,jte+1,kts:kte)
148 grid%xa%v(its:ite,jte-1,kts:kte)=grid%xa%v(its:ite,jte-1,kts:kte)-vtmp(its:ite,jte+1,kts:kte)
149 end if
150
151 grid%xa%v=grid%xa%v+vtmp
152
153 grid%a_u_2 = 0.0
154 grid%a_v_2 = 0.0
155
156 !---------------------------------------------------------------------------
157 ! [4.0] Adjoint of CONVERT TEMPERATURE inCREMENTS inTO THETA inCREMENTS
158 ! EVALUATE ALSO THE inCREMENTS OF (1/rho) AND GEOPOTENTIAL
159 !---------------------------------------------------------------------------
160
161 a_press(its:ite,jts:jte,kts:kte+1)=0.0
162 do k=kte,kts,-1
163 do j=jts,jte
164 do i=its,ite
165 rho_cgrid=-(grid%em_ph_2(i,j,k+1)-grid%em_ph_2(i,j,k))*grid%a_ph_2(i,j,k+1)/grid%xb%rho(i,j,k)
166 a_press(i,j,k )=a_press(i,j,k )+grid%a_ph_2(i,j,k+1)/grid%xb%rho(i,j,k)
167 a_press(i,j,k+1)=a_press(i,j,k+1)-grid%a_ph_2(i,j,k+1)/grid%xb%rho(i,j,k)
168 grid%a_ph_2(i,j,k ) =grid%a_ph_2(i,j,k) +grid%a_ph_2(i,j,k+1)
169 grid%xa%q(i,j,k)=grid%xa%q(i,j,k)-grid%xb%rho(i,j,k)*0.61*rho_cgrid/(1.0+0.61*grid%xb%q(i,j,k))
170 grid%xa%t(i,j,k)=grid%xa%t(i,j,k)-grid%xb%rho(i,j,k)*rho_cgrid/grid%xb%t(i,j,k)
171 grid%xa%p(i,j,k)= grid%xa%p(i,j,k)+grid%xb%rho(i,j,k)*rho_cgrid/grid%xb%p(i,j,k)
172 end do
173 end do
174 end do
175
176 do k=kts,kte
177 do j=jts,jte
178 do i=its,ite
179 grid%xa%p(i,j,k)=grid%xa%p(i,j,k)-(t0+grid%em_t_2(i,j,k))*kappa*grid%a_t_2(i,j,k)/grid%xb%p(i,j,k)
180 grid%xa%t(i,j,k)=grid%xa%t(i,j,k)+(t0+grid%em_t_2(i,j,k))*grid%a_t_2(i,j,k)/grid%xb%t(i,j,k)
181 end do
182 end do
183 end do
184
185 grid%a_t_2 = 0.0
186 grid%a_ph_2 = 0.0
187
188 !---------------------------------------------------------------------------
189 ! [3.0] Adjoint of COMPUTE pressure increments (for computing theta increments)
190 !---------------------------------------------------------------------------
191
192 do k=kts,kte
193 do j=jts,jte
194 do i=its,ite
195 a_press(i,j,k+1)=a_press(i,j,k+1)+0.5*grid%xa%p(i,j,k)
196 a_press(i,j,k )=a_press(i,j,k )+0.5*grid%xa%p(i,j,k)
197 grid%xa%p(i,j,k)=0.0
198 grid%a_moist(i,j,k,P_A_QV)=grid%a_moist(i,j,k,P_A_QV)-(grid%em_mu_2(i,j)+grid%em_mub(i,j))*a_press(i,j,k)*grid%em_dn(k)
199 grid%a_mu_2(i,j)=grid%a_mu_2(i,j)-a_press(i,j,k)*(1.0+grid%moist(i,j,k,P_QV))*grid%em_dn(k)
200 a_press(i,j,k+1)=a_press(i,j,k+1)+a_press(i,j,k)
201 end do
202 end do
203 end do
204
205 !---------------------------------------------------------------------------
206 ! [2.0] Adjoint of COMPUTE increments of dry-column air mass per unit area
207 !---------------------------------------------------------------------------
208
209 do j=jts,jte
210 do i=its,ite
211 sdmd=0.0
212 s1md=0.0
213 do k=kts,kte
214 s1md=s1md+(1.0+grid%moist(i,j,k,P_QV))*grid%em_dnw(k)
215 end do
216 sdmd=sdmd-grid%xb%psac(i,j)*grid%a_mu_2(i,j)/s1md
217 grid%xa%psfc(i,j)=grid%xa%psfc(i,j)-grid%a_mu_2(i,j)/s1md
218 do k=kts,kte
219 grid%a_moist(i,j,k,P_A_QV)=grid%a_moist(i,j,k,P_A_QV)+sdmd*grid%em_dnw(k)
220 end do
221 end do
222 end do
223
224 grid%a_mu_2 = 0.0
225 !---------------------------------------------------------------------------
226 ! [1.0] Adjoint of Get the mixing ratio of moisture
227 !---------------------------------------------------------------------------
228 do k=kts,kte
229 do j=jts,jte
230 do i=its,ite
231 grid%xa%q(i,j,k)=grid%xa%q(i,j,k)+grid%a_moist(i,j,k,P_A_QV)/(1.0-grid%xb%q(i,j,k))**2
232 end do
233 end do
234 end do
235
236 grid%a_moist = 0.0
237
238 if (trace_use) call da_trace_exit("da_transfer_xatowrftl_adj")
239
240 end subroutine da_transfer_xatowrftl_adj
241
242