Line data Source code
1 : module d2a3dikj_mod
2 :
3 : implicit none
4 : save
5 : private
6 : public :: d2a3dikj
7 :
8 : contains
9 :
10 : !-----------------------------------------------------------------------
11 : !BOP
12 : ! !ROUTINE: d2a3ikj -- Generalized 2nd order D-to-A grid transform (3D)
13 : ! Output array is i,k,j
14 : !
15 : ! !INTERFACE:
16 :
17 16128 : subroutine d2a3dikj(grid, am_geom_crrct, u, v, ua, va)
18 :
19 : ! !USES:
20 :
21 : use shr_kind_mod, only: r8 => shr_kind_r8
22 : use dynamics_vars, only : t_fvdycore_grid
23 :
24 : #if defined( SPMD )
25 : use parutilitiesmodule, only : parcollective, sumop
26 : use mod_comm, only: mp_send3d, mp_recv3d
27 : #endif
28 :
29 : use shr_reprosum_mod, only : shr_reprosum_calc, &
30 : shr_reprosum_tolExceeded, &
31 : shr_reprosum_reldiffmax, &
32 : shr_reprosum_recompute
33 : use cam_logfile, only : iulog
34 : use perf_mod
35 :
36 : implicit none
37 : ! !INPUT PARAMETERS:
38 : type (t_fvdycore_grid), intent(in) :: grid
39 : logical, intent(in) :: am_geom_crrct
40 : real(r8), intent(in) :: u(grid%ifirstxy:grid%ilastxy, &
41 : grid%jfirstxy:grid%jlastxy,grid%km) ! U-Wind
42 : real(r8), intent(in) :: v(grid%ifirstxy:grid%ilastxy, &
43 : grid%jfirstxy:grid%jlastxy,grid%km) ! V-Wind
44 :
45 : ! !INPUT/OUTPUT PARAMETERS:
46 : real(r8), intent(inout) :: ua(grid%ifirstxy:grid%ilastxy,grid%km, &
47 : grid%jfirstxy:grid%jlastxy) ! U-Wind
48 : real(r8), intent(inout) :: va(grid%ifirstxy:grid%ilastxy,grid%km, &
49 : grid%jfirstxy:grid%jlastxy) ! V-Wind
50 :
51 : ! !DESCRIPTION:
52 : !
53 : ! This routine performs a second order
54 : ! interpolation of three-dimensional wind
55 : ! fields on a D grid to an A grid. !
56 : !
57 : ! !REVISION HISTORY:
58 : ! WS 00.12.22 : Creation from d2a3d
59 : ! AAM 01.06.13 : Generalized to 2D decomposition
60 : ! WS 02.04.25 : New mod_comm interfaces
61 : ! WS 03.08.13 : Use unorth for ghosting U (aligned with d2a3dijk)
62 : ! WS 05.07.06 : Simplified interface with grid
63 : ! WS 06.09.08 : isolated magic numbers as F90 parameters
64 : ! PW 08.07.03 : introduced reprosum logic
65 : ! SS 12.10.29 : reprosum is now in csm_share
66 : !
67 : !EOP
68 : !-----------------------------------------------------------------------
69 : !BOC
70 : real(r8), parameter :: D0_0 = 0.0_r8
71 : real(r8), parameter :: D0_5 = 0.5_r8
72 :
73 : integer :: imh, i, j, k, m, itot, jtot, ltot, ik
74 32256 : real(r8) :: veast(grid%jfirstxy:grid%jlastxy,grid%km)
75 32256 : real(r8) :: unorth(grid%ifirstxy:grid%ilastxy,grid%km)
76 :
77 32256 : real(r8) :: uva(grid%ifirstxy:grid%ilastxy,grid%km,2)
78 32256 : real(r8) :: uvn(grid%km,2), uvs(grid%km,2)
79 32256 : real(r8) :: rel_diff(2,grid%km,2)
80 16128 : real(r8),allocatable :: uva_tmp(:)
81 :
82 : integer :: ifirstxy, ilastxy, jfirstxy, jlastxy, im, jm, km
83 : integer :: myidxy_y, myidxy_x, nprxy_x, iam
84 :
85 : logical :: write_warning
86 :
87 16128 : real(r8), pointer :: coslon(:),sinlon(:) ! Sine and cosine in longitude
88 :
89 : #if defined( SPMD )
90 : integer dest, src, incount, outcount
91 : #endif
92 :
93 : ! for AM correction
94 16128 : real(r8), pointer :: cosp(:), cose(:)
95 :
96 16128 : myidxy_y = grid%myidxy_y
97 16128 : myidxy_x = grid%myidxy_x
98 16128 : nprxy_x = grid%nprxy_x
99 16128 : iam = grid%iam
100 :
101 16128 : im = grid%im
102 16128 : jm = grid%jm
103 16128 : km = grid%km
104 :
105 16128 : ifirstxy = grid%ifirstxy
106 16128 : ilastxy = grid%ilastxy
107 16128 : jfirstxy = grid%jfirstxy
108 16128 : jlastxy = grid%jlastxy
109 :
110 16128 : itot = ilastxy-ifirstxy+1
111 16128 : jtot = jlastxy-jfirstxy+1
112 16128 : imh = im/2
113 :
114 16128 : coslon => grid%coslon
115 16128 : sinlon => grid%sinlon
116 16128 : cosp => grid%cosp
117 16128 : cose => grid%cose
118 :
119 : #if defined( SPMD )
120 : ! Set ua on A-grid
121 : call mp_send3d( grid%commxy, iam-nprxy_x, iam+nprxy_x, im, jm, km, &
122 : ifirstxy, ilastxy, jfirstxy, jlastxy, 1, km, &
123 16128 : ifirstxy, ilastxy, jfirstxy, jfirstxy, 1, km, u )
124 : call mp_recv3d( grid%commxy, iam+nprxy_x, im, jm, km, &
125 : ifirstxy, ilastxy, jlastxy+1, jlastxy+1, 1, km, &
126 16128 : ifirstxy, ilastxy, jlastxy+1, jlastxy+1, 1, km, unorth )
127 :
128 16128 : if ( jlastxy .lt. jm ) then
129 :
130 15876 : if (am_geom_crrct) then
131 : !$omp parallel do private(i, k)
132 0 : do k = 1, km
133 0 : do i = ifirstxy, ilastxy
134 0 : ua(i,k,jlastxy) = 0.5_r8*(u(i,jlastxy,k)*cose(jlastxy) + &
135 0 : unorth(i,k)*cose(jlastxy+1))/cosp(jlastxy)
136 : end do
137 : end do
138 : else
139 : !$omp parallel do private(i, k)
140 1127196 : do k = 1, km
141 27798876 : do i = ifirstxy, ilastxy
142 27783000 : ua(i,k,jlastxy) = 0.5_r8*(u(i,jlastxy,k) + unorth(i,k))
143 : end do
144 : end do
145 : end if
146 : end if
147 : #endif
148 :
149 16128 : if (am_geom_crrct) then
150 : !$omp parallel do private(i,j,k)
151 0 : do k = 1, km
152 0 : do j = jfirstxy, jlastxy-1
153 0 : do i = ifirstxy, ilastxy
154 0 : if (cosp(j) .ne. 0.0_r8) then
155 0 : ua(i,k,j) = 0.5_r8*(u(i,j,k)*cose(j) + u(i,j+1,k)*cose(j+1))/cosp(j) ! preserve curl-free flow
156 : else
157 0 : ua(i,k,j) = (u(i,j,k)*cose(j)+u(i,j+1,k)*cose(j+1))/(cose(j)+cose(j+1))
158 : end if
159 : end do
160 : end do
161 : end do
162 : else
163 : !$omp parallel do private(i,j,k)
164 1145088 : do k = 1, km
165 3403008 : do j = jfirstxy, jlastxy-1
166 57576960 : do i = ifirstxy, ilastxy
167 56448000 : ua(i,k,j) = 0.5_r8*(u(i,j,k) + u(i,j+1,k)) ! preserve solid-body flow
168 : end do
169 : end do
170 : end do
171 : end if
172 :
173 : ! Set va on A-grid
174 :
175 : !$omp parallel do private(j,k)
176 :
177 1145088 : do k = 1,km
178 4531968 : do j=jfirstxy,jlastxy
179 4515840 : veast(j,k) = v(ifirstxy,j,k)
180 : enddo
181 : enddo
182 :
183 : #if defined( SPMD )
184 16128 : if (itot .ne. im) then
185 16128 : dest = myidxy_y*nprxy_x + MOD(iam+nprxy_x-1,nprxy_x)
186 16128 : src = myidxy_y*nprxy_x + MOD(iam+1,nprxy_x)
187 : call mp_send3d( grid%commxy, dest, src, im, jm, km, &
188 : ifirstxy, ilastxy, jfirstxy, jlastxy, 1, km, &
189 16128 : ifirstxy, ifirstxy, jfirstxy, jlastxy, 1, km, v )
190 : call mp_recv3d( grid%commxy, src, im, jm, km, &
191 : ilastxy+1, ilastxy+1, jfirstxy, jlastxy, 1, km, &
192 16128 : ilastxy+1, ilastxy+1, jfirstxy, jlastxy, 1, km, veast )
193 : endif
194 : #endif
195 :
196 : !$omp parallel do private(i,j,k)
197 :
198 1145088 : do k=1,km
199 4531968 : do j=jfirstxy, jlastxy
200 81285120 : do i=ifirstxy,ilastxy-1
201 81285120 : va(i,k,j) = D0_5*(v(i,j,k) + v(i+1,j,k))
202 : enddo
203 4515840 : va(ilastxy,k,j) = D0_5*(v(ilastxy,j,k) + veast(j,k))
204 : enddo
205 : enddo
206 :
207 16128 : if (jfirstxy .eq. 1) then
208 : ! Was (something like) ...
209 : ! do i=1,imh
210 : ! us(k) = us(k) + (uvaglob(i+imh,k,1)-uvaglob(i,k,1))*sinlon(i) &
211 : ! + (uvaglob(i,k,2)-uvaglob(i+imh,k,2))*coslon(i)
212 : ! vs(k) = vs(k) + (uvaglob(i+imh,k,1)-uvaglob(i,k,1))*coslon(i) &
213 : ! + (uvaglob(i+imh,k,2)-uvaglob(i,k,2))*sinlon(i)
214 : ! enddo
215 :
216 : !$omp parallel do private(i,k)
217 17892 : do k = 1,km
218 229320 : do i=ifirstxy,min(imh,ilastxy)
219 211680 : uva(i,k,1) = -ua(i,k,2)*sinlon(i) + va(i,k,2)*coslon(i)
220 229320 : uva(i,k,2) = -ua(i,k,2)*coslon(i) - va(i,k,2)*sinlon(i)
221 : enddo
222 229572 : do i=max(imh+1,ifirstxy),ilastxy
223 211680 : uva(i,k,1) = ua(i,k,2)*sinlon(i-imh) - va(i,k,2)*coslon(i-imh)
224 229320 : uva(i,k,2) = ua(i,k,2)*coslon(i-imh) + va(i,k,2)*sinlon(i-imh)
225 : enddo
226 : enddo
227 :
228 252 : call t_startf("d2a3dikj_reprosum")
229 : call shr_reprosum_calc(uva, uvs, itot, itot, 2*km, gbl_count=im, &
230 252 : commid=grid%commxy_x, rel_diff=rel_diff)
231 252 : call t_stopf("d2a3dikj_reprosum")
232 :
233 : ! check that "fast" reproducible sum is accurate enough. If not, calculate
234 : ! using old method
235 252 : write_warning = .false.
236 252 : if (myidxy_x == 0) write_warning = .true.
237 252 : if ( shr_reprosum_tolExceeded('d2a3dikj/South Pole', 2*km, write_warning, &
238 : iulog, rel_diff) ) then
239 0 : if ( shr_reprosum_recompute ) then
240 0 : call t_startf("d2a3dikj_sumfix")
241 0 : allocate( uva_tmp(im) )
242 0 : do m = 1,2
243 0 : do k = 1,km
244 0 : if (rel_diff(1,k,m) > shr_reprosum_reldiffmax) then
245 0 : uva_tmp(:) = D0_0
246 0 : do i = ifirstxy,ilastxy
247 0 : uva_tmp(i) = uva(i,k,m)
248 : enddo
249 : #if defined(SPMD)
250 0 : call parcollective(grid%commxy_x,sumop,im,uva_tmp)
251 : #endif
252 0 : uvs(k,m) = D0_0
253 0 : do i = 1,im
254 0 : uvs(k,m) = uvs(k,m) + uva_tmp(i)
255 : enddo
256 : endif
257 : enddo
258 : enddo
259 0 : deallocate( uva_tmp )
260 0 : call t_stopf("d2a3dikj_sumfix")
261 : endif
262 : endif
263 :
264 : !$omp parallel do private(i,k)
265 17892 : do k = 1,km
266 17640 : uvs(k,1) = uvs(k,1)/im
267 17640 : uvs(k,2) = uvs(k,2)/im
268 229320 : do i=ifirstxy,min(imh,ilastxy)
269 211680 : ua(i,k,1) = -uvs(k,1)*sinlon(i) - uvs(k,2)*coslon(i)
270 229320 : va(i,k,1) = uvs(k,1)*coslon(i) - uvs(k,2)*sinlon(i)
271 : enddo
272 229572 : do i=max(imh+1,ifirstxy),ilastxy
273 211680 : ua(i,k,1) = uvs(k,1)*sinlon(i-imh) + uvs(k,2)*coslon(i-imh)
274 229320 : va(i,k,1) = -uvs(k,1)*coslon(i-imh) + uvs(k,2)*sinlon(i-imh)
275 : enddo
276 : enddo
277 :
278 : endif
279 :
280 16128 : if (jlastxy .eq. jm) then
281 : ! Was (something like) ...
282 : ! do i=1,imh
283 : ! un(k) = un(k) + (uaglob(i+imh,k,jm-1)-uaglob(i,k,jm-1))*sinlon(i) &
284 : ! + (vaglob(i+imh,k,jm-1)-vaglob(i,k,jm-1))*coslon(i)
285 : ! vn(k) = vn(k) + (uaglob(i,k,jm-1)-uaglob(i+imh,k,jm-1))*coslon(i) &
286 : ! + (vaglob(i+imh,k,jm-1)-vaglob(i,k,jm-1))*sinlon(i)
287 : ! enddo
288 :
289 : !$omp parallel do private(i,k)
290 17892 : do k = 1,km
291 229320 : do i=ifirstxy,min(imh,ilastxy)
292 211680 : uva(i,k,1) = -ua(i,k,jm-1)*sinlon(i) - va(i,k,jm-1)*coslon(i)
293 229320 : uva(i,k,2) = ua(i,k,jm-1)*coslon(i) - va(i,k,jm-1)*sinlon(i)
294 : enddo
295 229572 : do i=max(imh+1,ifirstxy),ilastxy
296 211680 : uva(i,k,1) = ua(i,k,jm-1)*sinlon(i-imh) + va(i,k,jm-1)*coslon(i-imh)
297 229320 : uva(i,k,2) = -ua(i,k,jm-1)*coslon(i-imh) + va(i,k,jm-1)*sinlon(i-imh)
298 : enddo
299 : enddo
300 :
301 252 : call t_startf("d2a3dikj_reprosum")
302 : call shr_reprosum_calc(uva, uvn, itot, itot, 2*km, gbl_count=im, &
303 252 : commid=grid%commxy_x, rel_diff=rel_diff)
304 252 : call t_stopf("d2a3dikj_reprosum")
305 :
306 : ! check that "fast" reproducible sum is accurate enough. If not, calculate
307 : ! using old method
308 252 : write_warning = .false.
309 252 : if (myidxy_x == 0) write_warning = .true.
310 252 : if ( shr_reprosum_tolExceeded('d2a3dikj/Nouth Pole', 2*km, write_warning, &
311 : iulog, rel_diff) ) then
312 0 : if ( shr_reprosum_recompute ) then
313 0 : call t_startf("d2a3dikj_sumfix")
314 0 : allocate( uva_tmp(im) )
315 0 : do m = 1,2
316 0 : do k = 1,km
317 0 : if (rel_diff(1,k,m) > shr_reprosum_reldiffmax) then
318 0 : uva_tmp(:) = D0_0
319 0 : do i = ifirstxy,ilastxy
320 0 : uva_tmp(i) = uva(i,k,m)
321 : enddo
322 : #if defined(SPMD)
323 0 : call parcollective(grid%commxy_x,sumop,im,uva_tmp)
324 : #endif
325 0 : uvn(k,m) = D0_0
326 0 : do i = 1,im
327 0 : uvn(k,m) = uvn(k,m) + uva_tmp(i)
328 : enddo
329 : endif
330 : enddo
331 : enddo
332 0 : deallocate( uva_tmp )
333 0 : call t_stopf("d2a3dikj_sumfix")
334 : endif
335 : endif
336 :
337 : !$omp parallel do private(i,k)
338 17892 : do k = 1,km
339 17640 : uvn(k,1) = uvn(k,1)/im
340 17640 : uvn(k,2) = uvn(k,2)/im
341 229320 : do i=ifirstxy,min(imh,ilastxy)
342 211680 : ua(i,k,jm) = -uvn(k,1)*sinlon(i) + uvn(k,2)*coslon(i)
343 229320 : va(i,k,jm) = -uvn(k,1)*coslon(i) - uvn(k,2)*sinlon(i)
344 : enddo
345 229572 : do i=max(imh+1,ifirstxy),ilastxy
346 211680 : ua(i,k,jm) = uvn(k,1)*sinlon(i-imh) - uvn(k,2)*coslon(i-imh)
347 229320 : va(i,k,jm) = uvn(k,1)*coslon(i-imh) + uvn(k,2)*sinlon(i-imh)
348 : enddo
349 : enddo
350 :
351 : endif
352 :
353 16128 : return
354 : !EOC
355 32256 : end subroutine d2a3dikj
356 : !-----------------------------------------------------------------------
357 : end module d2a3dikj_mod
|
