da_psichi_to_uv.inc
References to this file elsewhere.
1 subroutine da_psichi_to_uv(psi, chi, coefx,coefy, u, v)
2
3 !---------------------------------------------------------------------------
4 ! Purpose: Calculate wind components u and v from psi and chi.
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 implicit none
15
16 real, intent(inout) :: psi(ims:ime,jms:jme,kms:kme) ! Stream function
17 real, intent(inout) :: chi(ims:ime,jms:jme,kms:kme) ! Velocity potential
18 real, intent(in) :: coefx(ims:ime,jms:jme) ! Multiplicative coeff.
19 real, intent(in) :: coefy(ims:ime,jms:jme) ! Multiplicative coeff.
20 real, intent(out) :: u(ims:ime,jms:jme,kms:kme) ! u wind comp.
21 real, intent(out) :: v(ims:ime,jms:jme,kms:kme) ! v wind comp.
22
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 if (trace_use) call da_trace_entry("da_psichi_to_uv")
29
30 !---------------------------------------------------------------------------
31 ! [1.0] For Global application, set Wast-Eest Periodic boundary
32 !---------------------------------------------------------------------------
33
34 ! Computation to check for edge of domain:
35 is = its
36 ie = ite
37 js = jts
38 je = jte
39 if (jts == jds) js = jds+1
40 if (jte == jde) je = jde-1
41
42 if (global) then
43 call da_set_boundary_3d(psi)
44 call da_set_boundary_3d(chi)
45 else
46 if (its == ids) is = ids+1
47 if (ite == ide) ie = ide-1
48 end if
49
50 do k = kts, kte
51 !------------------------------------------------------------------------
52 ! [2.0] Compute u, v at interior points (2nd order central finite diffs):
53 !------------------------------------------------------------------------
54
55 do j = js, je
56 do i = is, ie
57 u(i,j,k) = -coefy(i,j)*(psi(i ,j+1,k) - psi(i ,j-1,k)) + &
58 coefx(i,j)*(chi(i+1,j ,k) - chi(i-1,j ,k))
59
60 v(i,j,k) = coefx(i,j)*(psi(i+1,j ,k) - psi(i-1,j ,k)) + &
61 coefy(i,j)*(chi(i ,j+1,k) - chi(i ,j-1,k))
62 end do
63 end do
64
65 if (global) cycle
66 !------------------------------------------------------------------------
67 ! [3.0] Compute u, v at domain boundaries:
68 !------------------------------------------------------------------------
69
70 ! [3.1] Western boundaries:
71
72 if (its == ids) then
73 i = its
74 do j = js, je
75 u(i,j,k) = -coefy(i,j)*(psi(i ,j+1,k) - psi(i,j-1,k)) + &
76 coefx(i,j)*(chi(i+2,j ,k) - chi(i,j ,k))
77
78 v(i,j,k) = coefx(i,j)*(psi(i+2,j ,k) - psi(i,j ,k)) + &
79 coefy(i,j)*(chi(i ,j+1,k) - chi(i,j-1,k))
80 end do
81 end if
82
83 ! [3.2] Eastern boundaries:
84
85 if (ite == ide) then
86 i = ite
87 do j = js, je
88 u(i,j,k) = -coefy(i,j)*(psi(i,j+1,k) - psi(i ,j-1,k)) + &
89 coefx(i,j)*(chi(i,j ,k) - chi(i-2,j ,k))
90
91 v(i,j,k) = coefx(i,j)*(psi(i,j ,k) - psi(i-2,j ,k))+ &
92 coefy(i,j)*(chi(i,j+1,k) - chi(i ,j-1,k))
93 end do
94 end if
95
96 ! [3.3] Southern boundaries:
97
98 if (jts == jds) then
99 j = jts
100 do i = is, ie
101 u(i,j,k) = -coefy(i,j)*(psi(i ,j+2,k) - psi(i ,j,k)) + &
102 coefx(i,j)*(chi(i+1,j ,k) - chi(i-1,j,k))
103
104 v(i,j,k) = coefx(i,j)*(psi(i+1,j ,k) - psi(i-1,j,k)) + &
105 coefy(i,j)*(chi(i ,j+2,k) - chi(i ,j,k))
106 end do
107 end if
108
109 ! [3.4] Northern boundaries:
110
111 if (jte == jde) then
112 j = jte
113 do i = is, ie
114 u(i,j,k) = -coefy(i,j)*(psi(i ,j,k) - psi(i ,j-2,k)) + &
115 coefx(i,j)*(chi(i+1,j,k) - chi(i-1,j ,k))
116
117 v(i,j,k) = coefx(i,j)*(psi(i+1,j,k) - psi(i-1,j ,k))+ &
118 coefy(i,j)*(chi(i ,j,k) - chi(i ,j-2,k))
119 end do
120 end if
121
122 !------------------------------------------------------------------------
123 ! [4.0] Corner points (assume average of surrounding points - poor?):
124 !------------------------------------------------------------------------
125
126 ! [4.1] Bottom-left point:
127
128 if (its == ids .AND. jts == jds) then
129 u(its,jts,k) = 0.5 * (u(its+1,jts,k) + u(its,jts+1,k))
130 v(its,jts,k) = 0.5 * (v(its+1,jts,k) + v(its,jts+1,k))
131 end if
132
133 ! [4.2] Top-left point:
134
135 if (ite == ide .AND. jts == jds) then
136 u(ite,jts,k) = 0.5 * (u(ite-1,jts,k) + u(ite,jts+1,k))
137 v(ite,jts,k) = 0.5 * (v(ite-1,jts,k) + v(ite,jts+1,k))
138 end if
139
140 ! [4.3] Bottom-right point:
141
142 if (its == ids .AND. jte == jde) then
143 u(its,jte,k) = 0.5 * (u(its+1,jte,k) + u(its,jte-1,k))
144 v(its,jte,k) = 0.5 * (v(its+1,jte,k) + v(its,jte-1,k))
145 end if
146
147 ! [4.4] Top-right point:
148
149 if (ite == ide .AND. jte == jde) then
150 u(ite,jte,k) = 0.5 * (u(ite-1,jte,k) + u(ite,jte-1,k))
151 v(ite,jte,k) = 0.5 * (v(ite-1,jte,k) + v(ite,jte-1,k))
152 end if
153 end do
154
155 !---------------------------------------------------------------------------
156 ! [5.0] For Global application, set Wast-Eest Periodic boundary
157 !---------------------------------------------------------------------------
158
159 if (global) then
160 call da_set_boundary_3d(u)
161 call da_set_boundary_3d(v)
162 end if
163
164 if (trace_use) call da_trace_exit("da_psichi_to_uv")
165
166 end subroutine da_psichi_to_uv
167
168