Line data Source code
1 : module diag_module
2 :
3 : ! diagnostic calcs
4 : !
5 : ! REVISION HISTORY:
6 : ! 05.09.10 Rasch Creation of compute_vdot_gradp
7 : ! 05.10.18 Sawyer Revisions for 2D decomp, placed in module
8 : ! 07.01.29 Chen Removed pft2d calculation for OMGA (is in cd_core)
9 :
10 : use shr_kind_mod, only: r8 => shr_kind_r8
11 : use spmd_utils, only: masterproc
12 : use pmgrid, only: plon, plev, plevp
13 : use cam_logfile, only: iulog
14 : use cam_history, only: addfld, outfld, add_default, horiz_only
15 : use cam_abortutils, only: endrun
16 :
17 : implicit none
18 : private
19 : save
20 :
21 : public :: &
22 : fv_diag_init, &
23 : fv_diag_am_calc, &
24 : compute_vdot_gradp
25 :
26 : real(r8) :: rplon
27 : real(r8) :: iref_p(plevp) ! interface reference pressure for vertical interpolation
28 : integer :: ip_b ! level index where hybrid levels become purely pressure
29 : integer :: zm_limit
30 :
31 : !========================================================================================
32 : CONTAINS
33 : !========================================================================================
34 :
35 0 : subroutine fv_diag_init()
36 :
37 : use hycoef, only : hyai, hybi, ps0
38 :
39 : ! local variables
40 : integer :: k
41 : logical :: history_waccm
42 : !---------------------------------------------------------------------------
43 :
44 0 : rplon = 1._r8/plon
45 0 : zm_limit = plon/3
46 :
47 : !-------------------------------------------------------------
48 : ! Calculate reference pressure
49 : !-------------------------------------------------------------
50 0 : do k = 1, plevp
51 0 : iref_p(k) = (hyai(k) + hybi(k)) * ps0
52 : end do
53 0 : if( masterproc ) then
54 0 : write(iulog,*) 'fv_diag_inti: iref_p'
55 0 : write(iulog,'(1p5g15.7)') iref_p(:)
56 : end if
57 :
58 : !-------------------------------------------------------------
59 : ! Find level where hybrid levels become purely pressure
60 : !-------------------------------------------------------------
61 0 : ip_b = -1
62 0 : do k = 1,plev
63 0 : if( hybi(k) == 0._r8 ) ip_b = k
64 : end do
65 :
66 : ! Fields for diagnosing angular momentum conservation. They are supplemental
67 : ! to the fields computed by do_circulation_diags
68 :
69 0 : call addfld ('dUzm' ,(/ 'ilev' /),'A','M/S','Zonal-Mean U dyn increm - defined on ilev', gridname='fv_centers_zonal' )
70 0 : call addfld ('dVzm' ,(/ 'ilev' /),'A','M/S','Zonal-Mean V dyn increm - defined on ilev', gridname='fv_centers_zonal' )
71 0 : call addfld ('dUazm',(/ 'ilev' /),'A','M/S','Zonal-Mean U adv increm - defined on ilev', gridname='fv_centers_zonal' )
72 0 : call addfld ('dVazm',(/ 'ilev' /),'A','M/S','Zonal-Mean V adv increm - defined on ilev', gridname='fv_centers_zonal' )
73 0 : call addfld ('dUfzm',(/ 'ilev' /),'A','M/S','Zonal-Mean U fixer incr - defined on ilev', gridname='fv_centers_zonal' )
74 0 : call addfld ('dU', (/ 'lev' /), 'A','K', 'U dyn increm', gridname='fv_centers' )
75 0 : call addfld ('dV', (/ 'lev' /), 'A','K', 'V dyn increm', gridname='fv_centers' )
76 0 : call addfld ('dUa', (/ 'lev' /), 'A','K', 'U adv increm', gridname='fv_centers' )
77 0 : call addfld ('dVa', (/ 'lev' /), 'A','K', 'V adv increm', gridname='fv_centers' )
78 0 : call addfld ('dUf', (/ 'lev' /), 'A','K', 'U fixer incr', gridname='fv_centers' )
79 :
80 0 : call add_default ('dUzm' ,1, ' ')
81 0 : call add_default ('dVzm' ,1, ' ')
82 0 : call add_default ('dUazm' ,1, ' ')
83 0 : call add_default ('dVazm' ,1, ' ')
84 0 : call add_default ('dUfzm' ,1, ' ')
85 0 : call add_default ('dU' , 1, ' ')
86 0 : call add_default ('dV' , 1, ' ')
87 0 : call add_default ('dUa', 1, ' ')
88 0 : call add_default ('dVa', 1, ' ')
89 0 : call add_default ('dUf', 1, ' ')
90 :
91 0 : end subroutine fv_diag_init
92 :
93 : !========================================================================================
94 :
95 0 : subroutine fv_diag_am_calc(grid, ps, pe, du3, dv3, dua3, dva3, duf3)
96 :
97 : ! Compute fields for diagnosing angular momentum conservation. They are supplemental
98 : ! to the fields computed by do_circulation_diags
99 :
100 : use dynamics_vars, only : T_FVDYCORE_GRID
101 : use interpolate_data, only : vertinterp
102 : #ifdef SPMD
103 : use mpishorthand, only : mpilog, mpiint
104 : use parutilitiesmodule, only : pargatherint
105 : #endif
106 :
107 : !-------------------------------------------------------------
108 : ! ... dummy arguments
109 : !-------------------------------------------------------------
110 : type(T_FVDYCORE_GRID), intent(in) :: grid ! FV Dynamics grid
111 : real(r8), intent(in) :: ps(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy) ! surface pressure (pa)
112 : real(r8), intent(in) :: pe(grid%ifirstxy:grid%ilastxy,plevp,grid%jfirstxy:grid%jlastxy) ! interface pressure (pa)
113 : real(r8), intent(in) :: du3 (grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! U increment, total (m/s/timestep)
114 : real(r8), intent(in) :: dv3 (grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! V increment, total (m/s/timestep)
115 : real(r8), intent(in) :: dua3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! U increment, advec (m/s/timestep)
116 : real(r8), intent(in) :: dva3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! V increment, advec (m/s/timestep)
117 : real(r8), intent(in) :: duf3(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! U increment, fixer (m/s/timestep)
118 :
119 : !-------------------------------------------------------------
120 : ! ... local variables
121 : !-------------------------------------------------------------
122 :
123 : real(r8) :: pinterp
124 0 : real(r8) :: pm(grid%ifirstxy:grid%ilastxy,plev,grid%jfirstxy:grid%jlastxy) ! mid-point pressure
125 : real(r8) :: psurf
126 :
127 0 : real(r8) :: dui (grid%ifirstxy:grid%ilastxy,plevp) ! interp. zonal mean U increment, total FV
128 0 : real(r8) :: dvi (grid%ifirstxy:grid%ilastxy,plevp) ! interp. zonal mean V increment, total FV
129 0 : real(r8) :: duai(grid%ifirstxy:grid%ilastxy,plevp) ! interp. zonal mean U increment, advection
130 0 : real(r8) :: dvai(grid%ifirstxy:grid%ilastxy,plevp) ! interp. zonal mean V increment, advection
131 0 : real(r8) :: dufi(grid%ifirstxy:grid%ilastxy,plevp) ! interp. zonal mean U increment, fixer
132 :
133 : real(r8) :: dum (plevp) ! zonal mean U increment, total FV
134 : real(r8) :: dvm (plevp) ! zonal mean V increment, total FV
135 : real(r8) :: duam(plevp) ! zonal mean U increment, advection
136 : real(r8) :: dvam(plevp) ! zonal mean V increment, advection
137 : real(r8) :: dufm(plevp) ! zonal mean U increment, fixer
138 :
139 : real(r8) :: rdiv(plevp)
140 :
141 0 : integer :: ip_gm1g(plon,grid%jfirstxy:grid%jlastxy) ! contains level index-1 where blocked points begin
142 : integer :: zm_cnt(plevp) ! counter
143 : integer :: i,j,k
144 : integer :: nlons
145 :
146 0 : logical :: has_zm(plevp,grid%jfirstxy:grid%jlastxy) ! .true. the (z,y) point is a valid zonal mean
147 0 : integer :: ip_gm1(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy) ! contains level index-1 where blocked points begin
148 :
149 0 : real(r8) :: du2d (plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged U dycore increment
150 0 : real(r8) :: dv2d (plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged V dycore increment
151 0 : real(r8) :: dua2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged U advect increment
152 0 : real(r8) :: dva2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged V advect increment
153 0 : real(r8) :: duf2d(plevp,grid%jfirstxy:grid%jlastxy) ! zonally averaged U fixer increment
154 :
155 : real(r8) :: tmp2(grid%ifirstxy:grid%ilastxy)
156 0 : real(r8) :: tmp3(grid%ifirstxy:grid%ilastxy,plevp)
157 0 : real(r8) :: tmph(grid%ifirstxy:grid%ilastxy,plev)
158 :
159 : integer :: beglat, endlat ! begin,end latitude indicies
160 : integer :: beglon, endlon ! begin,end longitude indicies
161 :
162 0 : beglon = grid%ifirstxy
163 0 : endlon = grid%ilastxy
164 0 : beglat = grid%jfirstxy
165 0 : endlat = grid%jlastxy
166 :
167 : !omp parallel do private (i,j,k,psurf)
168 : lat_loop1 : &
169 0 : do j = beglat, endlat
170 0 : do k = 1, plev
171 0 : do i = beglon, endlon
172 0 : pm(i,k,j) = 0.5_r8 * ( pe(i,k,j) + pe(i,k+1,j) )
173 : end do
174 : end do
175 : !-------------------------------------------------------------
176 : ! Keep track of where the bottom is in each column
177 : ! (i.e., largest index for which P(k) <= PS)
178 : !-------------------------------------------------------------
179 0 : ip_gm1(:,j) = plevp
180 0 : do i = beglon, endlon
181 0 : psurf = ps(i,j)
182 0 : do k = ip_b+1, plevp
183 0 : if( iref_p(k) <= psurf ) then
184 0 : ip_gm1(i,j) = k
185 : end if
186 : end do
187 : end do
188 : end do lat_loop1
189 :
190 0 : nlons = endlon - beglon + 1
191 :
192 : #ifdef SPMD
193 0 : if( grid%twod_decomp == 1 ) then
194 0 : if (grid%iam .lt. grid%npes_xy) then
195 0 : call pargatherint( grid%commxy_x, 0, ip_gm1, grid%strip2dx, ip_gm1g )
196 : endif
197 : else
198 0 : ip_gm1g(:,:) = ip_gm1(:,:)
199 : end if
200 : #else
201 : ip_gm1g(:,:) = ip_gm1(:,:)
202 : #endif
203 : #ifdef SPMD
204 0 : if( grid%myidxy_x == 0 ) then
205 : #endif
206 : lat_loop2 : &
207 0 : do j = beglat, endlat
208 0 : zm_cnt(:ip_b) = plon
209 0 : do k = ip_b+1, plevp
210 0 : zm_cnt(k) = count( ip_gm1g(:,j) >= k )
211 : end do
212 0 : has_zm(:ip_b,j) = .true.
213 0 : do k = ip_b+1, plevp
214 0 : has_zm(k,j) = zm_cnt(k) >= zm_limit
215 : end do
216 : end do lat_loop2
217 : #ifdef SPMD
218 : end if
219 0 : if( grid%twod_decomp == 1 ) then
220 0 : call mpibcast( has_zm, plevp*(endlat-beglat+1), mpilog, 0, grid%commxy_x )
221 0 : call mpibcast( ip_gm1g, plon*(endlat-beglat+1), mpiint, 0, grid%commxy_x )
222 : end if
223 : #endif
224 :
225 : lat_loop3 : &
226 0 : do j = beglat, endlat
227 : !-------------------------------------------------------------
228 : ! Vertical interpolation
229 : !-------------------------------------------------------------
230 0 : do k = 1, plevp
231 0 : pinterp = iref_p(k)
232 : !-------------------------------------------------------------
233 : ! Zonal & meridional velocity increments
234 : !-------------------------------------------------------------
235 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
236 0 : du3(beglon,1,j) , dui (beglon,k) )
237 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
238 0 : dv3(beglon,1,j) , dvi (beglon,k) )
239 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
240 0 : dua3(beglon,1,j), duai(beglon,k) )
241 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
242 0 : dva3(beglon,1,j), dvai(beglon,k) )
243 0 : call vertinterp( nlons, nlons, plev, pm(beglon,1,j), pinterp, &
244 0 : duf3(beglon,1,j), dufi(beglon,k) )
245 : end do
246 :
247 : !-------------------------------------------------------------
248 : ! Calculate zonal averages
249 : !-------------------------------------------------------------
250 0 : do k = ip_b+1, plevp
251 0 : where( ip_gm1(beglon:endlon,j) < k )
252 0 : dui (beglon:endlon,k)= 0._r8
253 : dvi (beglon:endlon,k)= 0._r8
254 : duai(beglon:endlon,k)= 0._r8
255 : dvai(beglon:endlon,k)= 0._r8
256 : dufi(beglon:endlon,k)= 0._r8
257 : endwhere
258 : end do
259 :
260 0 : call par_xsum(grid, dui, plevp, dum)
261 0 : call par_xsum(grid, dvi, plevp, dvm)
262 0 : call par_xsum(grid, duai, plevp, duam)
263 0 : call par_xsum(grid, dvai, plevp, dvam)
264 0 : call par_xsum(grid, dufi, plevp, dufm)
265 :
266 0 : do k = 1, ip_b
267 0 : du2d(k,j) = dum(k) * rplon
268 0 : dv2d(k,j) = dvm(k) * rplon
269 0 : dua2d(k,j)= duam(k)* rplon
270 0 : dva2d(k,j)= dvam(k)* rplon
271 0 : duf2d(k,j)= dufm(k)* rplon
272 : end do
273 :
274 0 : do k = ip_b+1, plevp
275 0 : if( has_zm(k,j) ) then
276 0 : rdiv(k) = 1._r8/count( ip_gm1g(:,j) >= k )
277 0 : du2d(k,j) = dum(k) * rdiv(k)
278 0 : dv2d(k,j) = dvm(k) * rdiv(k)
279 0 : dua2d(k,j)= duam(k)* rdiv(k)
280 0 : dva2d(k,j)= dvam(k)* rdiv(k)
281 0 : duf2d(k,j)= dufm(k)* rdiv(k)
282 : else
283 0 : du2d(k,j) = 0._r8
284 0 : dv2d(k,j) = 0._r8
285 0 : dua2d(k,j)= 0._r8
286 0 : dva2d(k,j)= 0._r8
287 0 : duf2d(k,j)= 0._r8
288 : end if
289 : end do
290 :
291 : end do lat_loop3
292 :
293 : !-------------------------------------------------------------
294 : ! Do the output
295 : !-------------------------------------------------------------
296 0 : latloop: do j = beglat,endlat
297 :
298 : !-------------------------------------------------------------
299 : ! zonal-mean output
300 : !-------------------------------------------------------------
301 :
302 0 : do k = 1,plevp
303 0 : tmp3(grid%ifirstxy,k) = du2d(k,j)
304 : enddo
305 0 : call outfld( 'dUzm', tmp3(grid%ifirstxy,:), 1, j )
306 :
307 0 : do k = 1,plevp
308 0 : tmp3(grid%ifirstxy,k) = dv2d(k,j)
309 : enddo
310 0 : call outfld( 'dVzm', tmp3(grid%ifirstxy,:), 1, j )
311 :
312 0 : do k = 1,plevp
313 0 : tmp3(grid%ifirstxy,k) = dua2d(k,j)
314 : enddo
315 0 : call outfld( 'dUazm', tmp3(grid%ifirstxy,:), 1, j )
316 :
317 0 : do k = 1,plevp
318 0 : tmp3(grid%ifirstxy,k) = dva2d(k,j)
319 : enddo
320 0 : call outfld( 'dVazm', tmp3(grid%ifirstxy,:), 1, j )
321 :
322 0 : do k = 1,plevp
323 0 : tmp3(grid%ifirstxy,k) = duf2d(k,j)
324 : enddo
325 0 : call outfld( 'dUfzm', tmp3(grid%ifirstxy,:), 1, j )
326 :
327 : !-------------------------------------------------------------
328 : ! 3D output
329 : !-------------------------------------------------------------
330 :
331 0 : do k = 1,plev
332 0 : do i = beglon,endlon
333 0 : tmph(i,k) = du3(i,k,j)
334 : enddo
335 : enddo
336 0 : call outfld( 'dU', tmph, nlons, j )
337 :
338 0 : do k = 1,plev
339 0 : do i = beglon,endlon
340 0 : tmph(i,k) = dv3(i,k,j)
341 : enddo
342 : enddo
343 0 : call outfld( 'dV', tmph, nlons, j )
344 :
345 0 : do k = 1,plev
346 0 : do i = beglon,endlon
347 0 : tmph(i,k) = dua3(i,k,j)
348 : enddo
349 : enddo
350 0 : call outfld( 'dUa', tmph, nlons, j )
351 :
352 0 : do k = 1,plev
353 0 : do i = beglon,endlon
354 0 : tmph(i,k) = dva3(i,k,j)
355 : enddo
356 : enddo
357 0 : call outfld( 'dVa', tmph, nlons, j )
358 :
359 0 : do k = 1,plev
360 0 : do i = beglon,endlon
361 0 : tmph(i,k) = duf3(i,k,j)
362 : enddo
363 : enddo
364 0 : call outfld( 'dUf', tmph, nlons, j )
365 :
366 : enddo latloop
367 :
368 0 : end subroutine fv_diag_am_calc
369 :
370 : !========================================================================================
371 :
372 16128 : subroutine compute_vdot_gradp(grid, dt, frac, cx, cy, pexy, omgaxy )
373 :
374 : use shr_kind_mod, only : r8 => shr_kind_r8
375 : use dynamics_vars, only : T_FVDYCORE_GRID
376 : #if defined( SPMD )
377 : use mod_comm, only: mp_send3d, mp_recv3d, &
378 : mp_sendirr, mp_recvirr
379 : #endif
380 :
381 : implicit none
382 :
383 : ! !INPUT PARAMETERS:
384 : type (T_FVDYCORE_GRID), intent(in) :: grid
385 : real(r8), intent(in):: dt
386 : real(r8), intent(in):: frac
387 :
388 : real(r8), intent(in):: cx(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
389 : real(r8), intent(in):: cy(grid%im,grid%jfirst:grid%jlast+1,grid%kfirst:grid%klast)
390 : real(r8), target, intent(in):: &
391 : pexy(grid%ifirstxy:grid%ilastxy,grid%km+1,grid%jfirstxy:grid%jlastxy) ! P (pascal) at layer edges
392 : real(r8), target, intent(inout):: &
393 : omgaxy(grid%ifirstxy:grid%ilastxy,grid%km,grid%jfirstxy:grid%jlastxy) ! vert. press. velocity (pa/sec)
394 :
395 : ! Local
396 : integer :: im ! dimension in east-west
397 : integer :: jm ! dimension in North-South
398 : integer :: km ! number of Lagrangian layers
399 : integer :: jfirst ! starting latitude index for MPI
400 : integer :: jlast ! ending latitude index for MPI
401 : integer :: kfirst ! starting level index for MPI
402 : integer :: klast ! ending level index for MPI
403 : integer :: js2g0 ! j==1 not included
404 : integer :: jn2g0 ! j==jm not included
405 :
406 32256 : real(r8) :: pm(grid%im, grid%jfirst-1:grid%jlast+1)
407 32256 : real(r8) :: grad(grid%im, grid%jfirst:grid%jlast+1)
408 : real(r8) :: fac, sum1
409 :
410 16128 : real(r8), pointer :: pe(:,:,:) ! YZ version of edge pressures
411 16128 : real(r8), pointer :: omga(:,:,:) ! YZ version of vert. vel.
412 :
413 : real(r8), parameter :: half = 0.5_r8
414 : real(r8), parameter :: zero = 0.0_r8
415 :
416 : integer :: i,j,k
417 :
418 : #if defined( SPMD )
419 : integer :: iam, dest, src, npr_y, npes_yz
420 32256 : real(r8) :: penorth(grid%im, grid%kfirst:grid%klast+1)
421 32256 : real(r8) :: pesouth(grid%im, grid%kfirst:grid%klast+1)
422 : #endif
423 :
424 16128 : im = grid%im
425 16128 : jm = grid%jm
426 16128 : km = grid%km
427 16128 : jfirst = grid%jfirst
428 16128 : jlast = grid%jlast
429 16128 : kfirst = grid%kfirst
430 16128 : klast = grid%klast
431 16128 : js2g0 = grid%js2g0
432 16128 : jn2g0 = grid%jn2g0
433 :
434 16128 : fac = half / (dt * frac)
435 :
436 : #if defined( SPMD )
437 16128 : if ( grid%twod_decomp == 1 ) then
438 80640 : allocate(pe(im,kfirst:klast+1,jfirst:jlast))
439 80640 : allocate(omga(im,kfirst:klast,jfirst:jlast))
440 : call mp_sendirr( grid%commxy, grid%ikj_xy_to_yz%SendDesc, &
441 : grid%ikj_xy_to_yz%RecvDesc, omgaxy, omga, &
442 16128 : modc=grid%modc_dynrun )
443 : call mp_recvirr( grid%commxy, grid%ikj_xy_to_yz%SendDesc, &
444 : grid%ikj_xy_to_yz%RecvDesc, omgaxy, omga, &
445 16128 : modc=grid%modc_dynrun )
446 : call mp_sendirr( grid%commxy, grid%pexy_to_pe%SendDesc, &
447 : grid%pexy_to_pe%RecvDesc, pexy, pe, &
448 16128 : modc=grid%modc_dynrun )
449 : call mp_recvirr( grid%commxy, grid%pexy_to_pe%SendDesc, &
450 : grid%pexy_to_pe%RecvDesc, pexy, pe, &
451 16128 : modc=grid%modc_dynrun )
452 : else
453 0 : pe => pexy
454 0 : omga => omgaxy
455 : endif
456 16128 : iam = grid%iam
457 16128 : npes_yz = grid%npes_yz
458 16128 : if (iam .lt. npes_yz) then
459 16128 : npr_y = grid%npr_y
460 16128 : dest = iam+1
461 16128 : src = iam-1
462 16128 : if ( mod(iam+1,npr_y) == 0 ) dest = -1
463 16128 : if ( mod(iam,npr_y) == 0 ) src = -1
464 :
465 : !
466 : ! Have to give more thought to the source and destination for 2D
467 : !
468 : call mp_send3d(grid%commyz, dest, src, im, km+1, jm, &
469 : 1, im, kfirst, klast+1, jfirst, jlast, &
470 16128 : 1, im, kfirst, klast+1, jlast, jlast, pe)
471 : call mp_recv3d(grid%commyz, src, im, km+1, jm, &
472 : 1, im, kfirst, klast+1, jfirst-1, jfirst-1, &
473 16128 : 1, im, kfirst, klast+1, jfirst-1, jfirst-1, pesouth)
474 : call mp_send3d(grid%commyz, src, dest, im, km+1, jm, &
475 : 1, im, kfirst, klast+1, jfirst, jlast, &
476 16128 : 1, im, kfirst, klast+1, jfirst, jfirst, pe)
477 : call mp_recv3d(grid%commyz, dest, im, km+1, jm, &
478 : 1, im, kfirst, klast+1, jlast+1, jlast+1, &
479 16128 : 1, im, kfirst, klast+1, jlast+1, jlast+1, penorth)
480 : end if ! (iam .lt. npes_yz)
481 : #else
482 : pe => pexy
483 : omga => omgaxy
484 : #endif
485 :
486 : !$omp parallel do private(i,j,k,pm,grad, sum1)
487 59136 : do k=kfirst,klast
488 :
489 : ! Compute layer mean p
490 172032 : do j=jfirst,jlast
491 37330944 : do i=1,im
492 37287936 : pm(i,j) = half * ( pe(i,k,j) + pe(i,k+1,j) )
493 : enddo
494 : enddo
495 :
496 : #if defined( SPMD )
497 43008 : if ( jfirst/=1 ) then
498 12235104 : do i=1,im
499 12235104 : pm(i,jfirst-1) = half * ( pesouth(i,k) + pesouth(i,k+1))
500 : enddo
501 : endif
502 :
503 43008 : if ( jlast/=jm ) then
504 12235104 : do i=1,im
505 12235104 : pm(i,jlast+1) = half * ( penorth(i,k) + penorth(i,k+1))
506 : enddo
507 : endif
508 : #endif
509 :
510 170688 : do j=js2g0,jn2g0
511 127680 : i=1
512 127680 : grad(i,j) = fac * cx(i,j,k) * (pm(i,j)-pm(im,j))
513 36814848 : do i=2,im
514 36771840 : grad(i,j) = fac * cx(i,j,k) * (pm(i,j)-pm(i-1,j))
515 : enddo
516 : enddo
517 :
518 170688 : do j=js2g0,jn2g0
519 36771840 : do i=1,im-1
520 36771840 : omga(i,k,j) = omga(i,k,j) + grad(i,j) + grad(i+1,j)
521 : enddo
522 127680 : i=im
523 170688 : omga(i,k,j) = omga(i,k,j) + grad(i,j) + grad(1,j)
524 : enddo
525 :
526 213696 : do j=js2g0,min(jm,jlast+1)
527 49371840 : do i=1,im
528 49328832 : grad(i,j) = fac * cy(i,j,k) * (pm(i,j)-pm(i,j-1))
529 : enddo
530 : enddo
531 :
532 170688 : do j=js2g0,jn2g0
533 36942528 : do i=1,im
534 36899520 : omga(i,k,j) = omga(i,k,j) + grad(i,j) + grad(i,j+1)
535 : enddo
536 : enddo
537 :
538 : ! Note: Since V*grad(P) at poles are harder to compute accurately we use the average of sourding points
539 : ! to be used as input to physics.
540 :
541 43008 : if ( jfirst==1 ) then
542 672 : sum1 = zero
543 194208 : do i=1,im
544 194208 : sum1 = sum1 + omga(i,k,2)
545 : enddo
546 672 : sum1 = sum1 / real(im,r8)
547 194208 : do i=1,im
548 194208 : omga(i,k,1) = sum1
549 : enddo
550 : endif
551 :
552 59136 : if ( jlast==jm ) then
553 672 : sum1 = zero
554 194208 : do i=1,im
555 194208 : sum1 = sum1 + omga(i,k,jm-1)
556 : enddo
557 672 : sum1 = sum1 / real(im,r8)
558 194208 : do i=1,im
559 194208 : omga(i,k,jm) = sum1
560 : enddo
561 : endif
562 : enddo
563 :
564 : #if defined( SPMD)
565 16128 : if ( grid%twod_decomp == 1 ) then
566 : !
567 : ! Transpose back to XY (if 1D, the changes to omgaxy were made in place)
568 : !
569 : call mp_sendirr( grid%commxy, grid%ikj_yz_to_xy%SendDesc, &
570 : grid%ikj_yz_to_xy%RecvDesc, omga, omgaxy, &
571 16128 : modc=grid%modc_dynrun )
572 : call mp_recvirr( grid%commxy, grid%ikj_yz_to_xy%SendDesc, &
573 : grid%ikj_yz_to_xy%RecvDesc, omga, omgaxy, &
574 16128 : modc=grid%modc_dynrun )
575 16128 : deallocate( pe )
576 16128 : deallocate( omga )
577 : endif
578 : #endif
579 :
580 16128 : end subroutine compute_vdot_gradp
581 :
582 : end module diag_module
|
