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