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