da_uv_to_divergence_adj.inc
References to this file elsewhere.
1 subroutine da_uv_to_divergence_adj(grid, u, v, div)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Adjoint of the subroutine da_uv_to_divergence
5 ! d U d V
6 ! Div = m^2 *[---(---) + ---(---) ]
7 ! dx m dy M
8 !---------------------------------------------------------------------------
9
10 implicit none
11
12 type (domain), intent(inout) :: grid
13
14 real, intent(out) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp.
15 real, intent(out) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp.
16 real, intent(inout):: div(ims:ime,jms:jme,kms:kme) ! Divergence.
17
18 integer :: i, j, k ! Loop counters.
19 integer :: is, ie ! 1st dim. end points.
20 integer :: js, je ! 2nd dim. end points.
21 real :: one_third ! 1/3.
22
23 real :: coeff, inv_2ds
24 real :: um(ims:ime,jms:jme) ! Temp. storage of u/m.
25 real :: vm(ims:ime,jms:jme) ! Temp. storage of v/m
26
27 if (trace_use) call da_trace_entry("da_uv_to_divergence_adj")
28
29 !---------------------------------------------------------------------------
30 ! [1.0] Initialise:
31 !---------------------------------------------------------------------------
32
33 one_third = 1.0 / 3.0
34
35 if (test_wrfvar) then
36 is = its; ie = ite; js = jts; je = jte
37 else
38 is = its - 1; ie = ite + 1; js = jts - 1; je = jte + 1
39 end if
40
41 if (.not. global .and. its == ids) is = ids+1
42 if (.not. global .and. ite == ide) ie = ide-1
43 if (jts == jds) js = jds+1; if (jte == jde) je = jde-1
44
45 !---------------------------------------------------------------------------
46 ! [2.0] Calculate divergence:
47 !---------------------------------------------------------------------------
48
49 if (.not. global) then
50 inv_2ds = 0.5 / grid%xb%ds
51 do k =kds, kde
52 um(ims:ime,jms:jme) = 0.0
53 vm(ims:ime,jms:jme) = 0.0
54 ! [2.2] Impose zero divergence gradient condition at boundaries:
55
56 ! [2.2.4] Right boundaries:
57
58 if (jte == jde) then
59 j = jte
60 do i = its, ite ! This is different to original
61 div(i,j-1,k)=div(i,j-1,k)+div(i,j,k)*one_third*4.0
62 div(i,j-2,k)=div(i,j-2,k)-div(i,j,k)*one_third
63 div(i,j,k)=0.0
64 end do
65 end if
66
67 ! [2.2.3] Left boundaries:
68
69 if (jts == jds) then
70 j = jts
71 do i = its, ite ! This is different to original
72 div(i,j+1,k)=div(i,j+1,k)+div(i,j,k)*one_third*4.0
73 div(i,j+2,k)=div(i,j+2,k)-div(i,j,k)*one_third
74 div(i,j,k)=0.0
75 end do
76 end if
77
78 ! [2.2.2] Top boundaries:
79
80 if (ite == ide) then
81 i = ite
82 do j = jts, jte
83 div(i-1,j,k)=div(i-1,j,k)+div(i,j,k)*one_third*4.0
84 div(i-2,j,k)=div(i-2,j,k)-div(i,j,k)*one_third
85 div(i,j,k)=0.0
86 end do
87 end if
88
89 ! [2.2.1] Bottom boundaries:
90
91 if (its == ids) then
92 i = its
93 do j = jts, jte
94 div(i+1,j,k)=div(i+1,j,k)+div(i,j,k)*one_third*4.0
95 div(i+2,j,k)=div(i+2,j,k)-div(i,j,k)*one_third
96 div(i,j,k)=0.0
97 end do
98 end if
99
100 ! [2.1] Compute fd divergence at interior points:
101 ! Computation to check for edge of domain:
102 ! This is only for adjoint, as we have to cross the processor boundary
103 ! to get the contribution.
104
105 if (.not. test_wrfvar) then
106 grid%xp%vxy(its:ite, jts:jte) = div(its:ite, jts:jte, k)
107 #ifdef DM_PARALLEL
108 #include "HALO_2D_WORK.inc"
109 #endif
110
111 div(is, js:je, k) = grid%xp%vxy(is, js:je)
112 div(ie, js:je, k) = grid%xp%vxy(ie, js:je)
113 div(is:ie, js, k) = grid%xp%vxy(is:ie, js)
114 div(is:ie, je, k) = grid%xp%vxy(is:ie, je)
115 end if
116
117 do j = js, je
118 do i = is, ie
119 coeff = grid%xb%map_factor(i,j) * grid%xb%map_factor(i,j) * inv_2ds
120 um(i+1,j)=um(i+1,j)+div(i,j,k)*coeff
121 um(i-1,j)=um(i-1,j)-div(i,j,k)*coeff
122 vm(i,j+1)=vm(i,j+1)+div(i,j,k)*coeff
123 vm(i,j-1)=vm(i,j-1)-div(i,j,k)*coeff
124 end do
125 end do
126
127 u(is-1:ie+1,js-1:je+1,k) = u(is-1:ie+1,js-1:je+1,k) +um(is-1:ie+1,js-1:je+1) / grid%xb%map_factor(is-1:ie+1,js-1:je+1)
128 v(is-1:ie+1,js-1:je+1,k) = v(is-1:ie+1,js-1:je+1,k) +vm(is-1:ie+1,js-1:je+1) / grid%xb%map_factor(is-1:ie+1,js-1:je+1)
129 end do
130
131 else ! global
132 call da_set_boundary_3d(div)
133
134 do k =kds, kde
135 !-------------------------------------------------------------------------
136 ! [2.1] Compute fd divergence at interior points:
137 !-------------------------------------------------------------------------
138
139 do j = je, js, -1
140 do i = ie, is, -1
141 u(i+1,j,k) = u(i+1,j,k) + grid%xb%coefx(i,j) * div(i,j,k)
142 u(i-1,j,k) = u(i-1,j,k) - grid%xb%coefx(i,j) * div(i,j,k)
143 v(i,j+1,k) = v(i,j+1,k) + grid%xb%coefy(i,j) * div(i,j,k)
144 v(i,j-1,k) = v(i,j-1,k) - grid%xb%coefy(i,j) * div(i,j,k)
145 end do
146 end do
147 end do
148 end if
149
150 div = 0.0
151
152 if (trace_use) call da_trace_exit("da_uv_to_divergence_adj")
153
154 end subroutine da_uv_to_divergence_adj
155
156