da_psichi_to_uv_adj.inc

References to this file elsewhere.
1 subroutine da_psichi_to_uv_adj(u, v, coefx, coefy, psi, chi)
2 
3    !---------------------------------------------------------------------------
4    !  Purpose: Adjoint code of da_psichi_to_uv  
5    !
6    !  Method:  u = coefx * (-dpsi/dy + dchi/dx)
7    !           v = coefy * ( dpsi/dx + dchi/dy)
8    !
9    !  Assumptions: Unstaggered grid.
10    !               Lateral boundary conditions - dpsi/dn, dchi/dn = 0 (FCT)
11    !               grid size may or may not be equal
12    !
13    !---------------------------------------------------------------------------
14 
15    implicit none
16    
17    real, intent(inout):: u(ims:ime,jms:jme,kms:kme)   ! u wind comp.
18    real, intent(inout):: v(ims:ime,jms:jme,kms:kme)   ! v wind comp.
19    real, intent(inout):: psi(ims:ime,jms:jme,kms:kme) ! Stream function
20    real, intent(inout):: chi(ims:ime,jms:jme,kms:kme) ! Velocity potential
21    real, intent(in)   :: coefx(ims:ime,jms:jme)       ! Multiplicative coeff.
22    real, intent(in)   :: coefy(ims:ime,jms:jme)       ! Multiplicative coeff.
23    
24    integer            :: i, j, k                      ! Loop counters.
25    integer            :: is, ie                       ! 1st dim. end points.
26    integer            :: js, je                       ! 2nd dim. end points.
27 
28    !---------------------------------------------------------------------------
29    ! [5.0] For Global application, set Wast-Eest Periodic boundary
30    !---------------------------------------------------------------------------
31 
32    if (global) then
33       call da_set_boundary_3d(u) 
34       call da_set_boundary_3d(v) 
35    end if
36 
37    !---------------------------------------------------------------------------
38    ! Computation to check for edge of domain:
39    !---------------------------------------------------------------------------
40 
41    is = its-1
42    ie = ite+1
43    js = jts-1
44    je = jte+1
45    if (jts == jds) js = jds+1
46    if (jte == jde) je = jde-1
47 
48    if (.not. global) then
49       if (its == ids) is = ids+1
50       if (ite == ide) ie = ide-1
51 
52       do k = kts, kte
53 
54          !---------------------------------------------------------------------
55          ! [4.0] Corner points (assume average of surrounding points - poor?):
56          !---------------------------------------------------------------------
57 
58          ! [4.1] Bottom-left point:
59 
60          if (its == ids .AND. jts == jds) then
61             u(its+1,jts,k) = u(its+1,jts,k) + 0.5 * u(its,jts,k)
62             u(its,jts+1,k) = u(its,jts+1,k) + 0.5 * u(its,jts,k)
63             v(its+1,jts,k) = v(its+1,jts,k) + 0.5 * v(its,jts,k)
64             v(its,jts+1,k) = v(its,jts+1,k) + 0.5 * v(its,jts,k)
65          end if
66 
67         ! [4.2] Top-left point:
68 
69          if (ite == ide .AND. jts == jds) then
70             u(ite-1,jts,k) = u(ite-1,jts,k) + 0.5 * u(ite,jts,k)
71             u(ite,jts+1,k) = u(ite,jts+1,k) + 0.5 * u(ite,jts,k)
72             v(ite-1,jts,k) = v(ite-1,jts,k) + 0.5 * v(ite,jts,k)
73             v(ite,jts+1,k) = v(ite,jts+1,k) + 0.5 * v(ite,jts,k)
74          end if
75 
76          ! [4.3] Bottom-right point:
77 
78          if (its == ids .AND. jte == jde) then
79             u(its+1,jte,k) = u(its+1,jte,k) + 0.5 * u(its,jte,k)
80             u(its,jte-1,k) = u(its,jte-1,k) + 0.5 * u(its,jte,k)
81             v(its+1,jte,k) = v(its+1,jte,k) + 0.5 * v(its,jte,k)
82             v(its,jte-1,k) = v(its,jte-1,k) + 0.5 * v(its,jte,k)
83          end if
84 
85          ! [4.4] Top-right point:
86 
87          if (ite == ide .AND. jte == jde) then
88             u(ite-1,jte,k) = u(ite-1,jte,k) + 0.5 * u(ite,jte,k)
89             u(ite,jte-1,k) = u(ite,jte-1,k) + 0.5 * u(ite,jte,k)
90             v(ite-1,jte,k) = v(ite-1,jte,k) + 0.5 * v(ite,jte,k)
91             v(ite,jte-1,k) = v(ite,jte-1,k) + 0.5 * v(ite,jte,k)
92          end if
93       end do
94    end if
95 
96    ! [3.0] Compute u, v at domain boundaries:
97 
98    do k = kts, kte
99       if (.not. global) then
100          ! [3.4] Northern boundaries:
101 
102          if (jte == jde) then
103             j = jte
104 
105             do i = ie, is, -1
106                psi(i+1,j,k) = psi(i+1,j,k) + coefx(i,j) * v(i,j,k)
107                psi(i-1,j,k) = psi(i-1,j,k) - coefx(i,j) * v(i,j,k)
108                chi(i,j  ,k) = chi(i,j  ,k) + coefy(i,j) * v(i,j,k)
109                chi(i,j-2,k) = chi(i,j-2,k) - coefy(i,j) * v(i,j,k)
110 
111                psi(i,j  ,k) = psi(i,j  ,k) - coefy(i,j) * u(i,j,k)
112                psi(i,j-2,k) = psi(i,j-2,k) + coefy(i,j) * u(i,j,k)
113                chi(i+1,j,k) = chi(i+1,j,k) + coefx(i,j) * u(i,j,k)
114                chi(i-1,j,k) = chi(i-1,j,k) - coefx(i,j) * u(i,j,k)
115             end do
116          end if
117 
118          ! [3.3] Southern boundaries:
119 
120          if (jts == jds) then
121             j = jts
122 
123             do i = ie, is, -1
124 
125 
126                psi(i+1,j,k) = psi(i+1,j,k) + coefx(i,j) * v(i,j,k)
127                psi(i-1,j,k) = psi(i-1,j,k) - coefx(i,j) * v(i,j,k)
128                chi(i,j+2,k) = chi(i,j+2,k) + coefy(i,j) * v(i,j,k)
129                chi(i,j  ,k) = chi(i,j  ,k) - coefy(i,j) * v(i,j,k)
130 
131                psi(i,j+2,k) = psi(i,j+2,k) - coefy(i,j) * u(i,j,k)
132                psi(i,j  ,k) = psi(i,j  ,k) + coefy(i,j) * u(i,j,k)
133                chi(i+1,j,k) = chi(i+1,j,k) + coefx(i,j) * u(i,j,k)
134                chi(i-1,j,k) = chi(i-1,j,k) - coefx(i,j) * u(i,j,k)
135 
136             end do
137          end if
138 
139          ! [3.2] Eastern boundaries:
140 
141          if (ite == ide) then
142             i = ite
143             do j = je, js, -1
144                psi(i  ,j,k) = psi(i  ,j,k) + coefx(i,j) * v(i,j,k)
145                psi(i-2,j,k) = psi(i-2,j,k) - coefx(i,j) * v(i,j,k)
146                chi(i,j+1,k) = chi(i,j+1,k) + coefy(i,j) * v(i,j,k)
147                chi(i,j-1,k) = chi(i,j-1,k) - coefy(i,j) * v(i,j,k)
148 
149                psi(i,j+1,k) = psi(i,j+1,k) - coefy(i,j) * u(i,j,k)
150                psi(i,j-1,k) = psi(i,j-1,k) + coefy(i,j) * u(i,j,k)
151                chi(i  ,j,k) = chi(i  ,j,k) + coefx(i,j) * u(i,j,k)
152                chi(i-2,j,k) = chi(i-2,j,k) - coefx(i,j) * u(i,j,k)
153             end do
154          end if
155 
156          ! [3.1] Western boundaries:
157          if (its == ids) then
158             i = its
159 
160             do j = je, js, -1
161                psi(i+2,j,k) = psi(i+2,j,k) + coefx(i,j) * v(i,j,k)
162                psi(i  ,j,k) = psi(i  ,j,k) - coefx(i,j) * v(i,j,k)
163                chi(i,j+1,k) = chi(i,j+1,k) + coefy(i,j) * v(i,j,k)
164                chi(i,j-1,k) = chi(i,j-1,k) - coefy(i,j) * v(i,j,k)
165 
166                psi(i,j+1,k) = psi(i,j+1,k) - coefy(i,j) * u(i,j,k)
167                psi(i,j-1,k) = psi(i,j-1,k) + coefy(i,j) * u(i,j,k)
168                chi(i+2,j,k) = chi(i+2,j,k) + coefx(i,j) * u(i,j,k)
169                chi(i,  j,k) = chi(i,  j,k) - coefx(i,j) * u(i,j,k)
170             end do
171          end if
172       end if
173 
174       !------------------------------------------------------------------------
175       ! [2.0] Compute u, v at interior points (2nd order central finite diffs):
176       !------------------------------------------------------------------------
177 
178       do j = je, js, -1
179          do i = ie, is, -1
180             psi(i+1,j,k) = psi(i+1,j,k) + coefx(i,j) * v(i,j,k)
181             psi(i-1,j,k) = psi(i-1,j,k) - coefx(i,j) * v(i,j,k)
182             chi(i,j+1,k) = chi(i,j+1,k) + coefy(i,j) * v(i,j,k)
183             chi(i,j-1,k) = chi(i,j-1,k) - coefy(i,j) * v(i,j,k)
184 
185             psi(i,j+1,k) = psi(i,j+1,k) - coefy(i,j) * u(i,j,k)
186             psi(i,j-1,k) = psi(i,j-1,k) + coefy(i,j) * u(i,j,k)
187             chi(i+1,j,k) = chi(i+1,j,k) + coefx(i,j) * u(i,j,k)
188             chi(i-1,j,k) = chi(i-1,j,k) - coefx(i,j) * u(i,j,k)
189          end do
190       end do
191    end do    !  loop over levels
192 
193    !---------------------------------------------------------------------------
194    ! [5.0] For Global application, set Wast-Eest Periodic boundary
195    !---------------------------------------------------------------------------
196 
197    if (global) then
198        call da_set_boundary_3d(psi)
199        call da_set_boundary_3d(chi)
200    end if
201 
202 end subroutine da_psichi_to_uv_adj
203 
204