Line data Source code
1 :
2 : module pfixer
3 :
4 : !-----------------------------------------------------------------------
5 : !
6 : ! BOP
7 : !
8 : ! !MODULE: pfixer
9 : !
10 : ! !DESCRIPTION
11 : ! Corrects (or fixes) mass fluxes and edge pressures to be consistent
12 : ! with offline surface pressure at next model time step.
13 : !
14 : ! !USES
15 : use shr_kind_mod, only: r8 => shr_kind_r8
16 : use cam_abortutils, only: endrun
17 : use dynamics_vars, only: T_FVDYCORE_GRID
18 : use cam_logfile, only: iulog
19 : use metdata, only: met_rlx
20 :
21 : ! !PUBLIC MEMBER FUNCTIONS
22 : public :: adjust_press
23 :
24 : ! !REVISION HISTORY:
25 : ! 04.01.30 Stacy Walters Creation
26 : ! 04.02.15 F Vitt Fixed bug in edge pressure corrections
27 : ! 04.08.27 F Vitt Added ability to handle 2D decomposition
28 : ! 05.07.06 Sawyer Simplified interfaces with grid
29 : !
30 : ! EOP
31 : !-----------------------------------------------------------------------
32 :
33 : private
34 : real(r8), parameter :: D0_0 = 0.0_r8
35 : real(r8), parameter :: D0_5 = 0.5_r8
36 : real(r8), parameter :: D1_0 = 1.0_r8
37 : real(r8), parameter :: D100_0 = 100.0_r8
38 : contains
39 :
40 : !-----------------------------------------------------------------------
41 : ! ... adjust mass fluxes and pressures for lin-rood transport
42 : !-----------------------------------------------------------------------
43 :
44 64512 : subroutine adjust_press( grid, ps_mod, ps_obs, mfx, mfy, pexy )
45 :
46 : #if defined( SPMD )
47 : use mod_comm, only : mp_send3d, mp_recv3d, &
48 : mp_sendirr, mp_recvirr
49 : use parutilitiesmodule, only: parcollective, SUMOP
50 : #endif
51 : use time_manager, only : get_nstep
52 : !!$ use history, only: outfld
53 :
54 : implicit none
55 :
56 : !-----------------------------------------------------------------------
57 : ! ... dummy arguments
58 : !-----------------------------------------------------------------------
59 : type (T_FVDYCORE_GRID), intent(in) :: grid
60 : real(r8), intent(in) :: ps_obs(grid%im,grid%jfirst:grid%jlast) ! surface pressure at t(n) (Pa)
61 : real(r8), intent(in) :: ps_mod(grid%im,grid%jfirst:grid%jlast) ! surface pressure at t(n) (Pa)
62 : real(r8), intent(inout) :: mfx(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast) ! zonal mass flux
63 : real(r8), intent(inout) :: mfy(grid%im,grid%jfirst:grid%jlast+1,grid%kfirst:grid%klast) ! meridional mass flux
64 : real(r8), intent(inout) :: pexy(grid%ifirstxy:grid%ilastxy,grid%km+1,grid%jfirstxy:grid%jlastxy)
65 :
66 : !-----------------------------------------------------------------------
67 : ! ... local variables
68 : !-----------------------------------------------------------------------
69 : integer, parameter :: nstep0 = -1
70 :
71 : integer :: i, j, k, km1
72 : integer :: nstep
73 : #if defined( SPMD )
74 : integer :: dest, src
75 : #endif
76 : integer :: ndx(2)
77 64512 : real(r8) :: dpi(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
78 64512 : real(r8) :: dpixy(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,grid%km)
79 64512 : real(r8) :: dpi_in(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
80 64512 : real(r8) :: dpi_inxy(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy,1:grid%km)
81 64512 : real(r8) :: dpi_c(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
82 64512 : real(r8) :: dps (grid%im,grid%jfirst:grid%jlast)
83 64512 : real(r8) :: dps_in(grid%im,grid%jfirst:grid%jlast)
84 64512 : real(r8) :: dps_inxy(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
85 64512 : real(r8) :: ps_diffxy(grid%ifirstxy:grid%ilastxy,grid%jfirstxy:grid%jlastxy)
86 :
87 64512 : real(r8) :: dmfx(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast) ! original zonnal mass flux
88 64512 : real(r8) :: dmfy(grid%im,grid%jfirst:grid%jlast+1,grid%kfirst:grid%klast) ! original meridional mass flux
89 64512 : real(r8) :: emfx(grid%im,grid%jfirst:grid%jlast) ! zonal mass flux error
90 64512 : real(r8) :: emfy(grid%im,grid%jfirst:grid%jlast+1) ! meridional mass flux error
91 :
92 : real(r8) :: tmp2d(grid%im,grid%kfirst:grid%klast)
93 : logical :: debug = .false.
94 : logical :: method1 = .true.
95 :
96 : integer :: im, jm, km, ifirstxy, ilastxy, jfirstxy, jlastxy
97 : integer :: jfirst, jlast, kfirst, klast
98 :
99 32256 : im = grid%im
100 32256 : jm = grid%jm
101 32256 : km = grid%km
102 :
103 32256 : ifirstxy = grid%ifirstxy
104 32256 : ilastxy = grid%ilastxy
105 32256 : jfirstxy = grid%jfirstxy
106 32256 : jlastxy = grid%jlastxy
107 :
108 32256 : jfirst = grid%jfirst
109 32256 : jlast = grid%jlast
110 32256 : kfirst = grid%kfirst
111 32256 : klast = grid%klast
112 :
113 : #if defined( SPMD )
114 : ! Send one latitude of mfy to the south
115 32256 : if( mod(grid%iam,grid%npr_y) /= 0 ) then
116 31752 : dest = grid%iam-1
117 : else
118 504 : dest = -1
119 : end if
120 32256 : if( mod(grid%iam+1,grid%npr_y) /= 0 ) then
121 31752 : src = grid%iam+1
122 : else
123 504 : src = -1
124 : end if
125 : call mp_send3d( grid%commxy, dest, src, im, jm, km, &
126 : 1, im, jfirst, jlast+1, kfirst, klast, &
127 32256 : 1, im, jfirst, jfirst, kfirst, klast, mfy)
128 : #endif
129 :
130 129024 : do j = jfirst,jlast
131 27998208 : dps(:,j) = ps_obs(:,j) - ps_mod(:,j)
132 : end do
133 :
134 : ! ghost mfy
135 : #if defined( SPMD )
136 : call mp_recv3d( grid%commxy, src, im, jm, km, &
137 : 1, im, jfirst, jlast+1, kfirst, klast, &
138 32256 : 1, im, jlast+1, jlast+1, kfirst, klast, mfy)
139 : #endif
140 :
141 32256 : nstep = get_nstep()
142 : !-----------------------------------------------------------------------
143 : ! ... store incoming mass fluxes
144 : !-----------------------------------------------------------------------
145 32256 : if (debug) then
146 0 : do k = kfirst,klast
147 0 : do j = jfirst,jlast
148 0 : dmfx(:,j,k) = mfx(:,j,k)
149 0 : dmfy(:,j,k) = mfy(:,j,k)
150 : end do
151 0 : if( jlast /= jm ) then
152 0 : dmfy(:,jlast+1,k) = mfy(:,jlast+1,k)
153 : end if
154 : end do
155 : endif
156 :
157 : !-----------------------------------------------------------------------
158 : ! ... incoming mass flux divergence
159 : !-----------------------------------------------------------------------
160 32256 : call calc_divergence( grid, mfx, mfy, dpi_in )
161 :
162 : !-----------------------------------------------------------------------
163 : ! ... surface pressure from mass flux divergence
164 : !-----------------------------------------------------------------------
165 : ! Two different methods to compute change in ps give differnt
166 : ! results if 2D decomp is used (round off error). Method 1 gives
167 : ! identical 1D vs 2D decomposition results.
168 : !-----------------------------------------------------------------------
169 32256 : if (method1) then
170 :
171 : ! xfer dpi_in to dpi_inxy
172 32256 : if (grid%twod_decomp .eq. 1) then
173 : #if defined (SPMD)
174 : call mp_sendirr( grid%commxy, grid%ijk_yz_to_xy%SendDesc, &
175 : grid%ijk_yz_to_xy%RecvDesc, dpi_in, dpi_inxy, &
176 32256 : modc=grid%modc_dynrun )
177 : call mp_recvirr( grid%commxy, grid%ijk_yz_to_xy%SendDesc, &
178 : grid%ijk_yz_to_xy%RecvDesc, dpi_in, dpi_inxy, &
179 32256 : modc=grid%modc_dynrun )
180 : #endif
181 : else
182 0 : dpi_inxy(:,:,:) = dpi_in(:,:,:)
183 : endif
184 :
185 : ! vertical sum
186 129024 : do j = jfirstxy,jlastxy
187 2451456 : do i = ifirstxy,ilastxy
188 132475392 : dps_inxy(i,j) = sum( dpi_inxy(i,j,1:km) )
189 : end do
190 : end do
191 :
192 : ! xfer dps_inxy to dps_in
193 : ! Embed in 3D array since transpose machinery cannot handle 2D arrays
194 32256 : if (grid%twod_decomp .eq. 1) then
195 : #if defined (SPMD)
196 1838592 : do k = 1,km
197 7257600 : do j = jfirstxy,jlastxy
198 137281536 : do i = ifirstxy,ilastxy
199 135475200 : dpixy(i,j,k) = dps_inxy(i,j)
200 : enddo
201 : enddo
202 : enddo
203 :
204 : call mp_sendirr( grid%commxy, grid%ijk_xy_to_yz%SendDesc, &
205 : grid%ijk_xy_to_yz%RecvDesc, dpixy, dpi, &
206 32256 : modc=grid%modc_dynrun )
207 : call mp_recvirr( grid%commxy, grid%ijk_xy_to_yz%SendDesc, &
208 : grid%ijk_xy_to_yz%RecvDesc, dpixy, dpi, &
209 32256 : modc=grid%modc_dynrun )
210 :
211 129024 : do j = jfirst,jlast
212 27998208 : do i = 1,im
213 27965952 : dps_in(i,j) = dpi(i,j,kfirst)
214 : enddo
215 : enddo
216 : #endif
217 : else
218 0 : dps_in(:,:) = dps_inxy(:,:)
219 : endif
220 :
221 : else ! method1
222 :
223 : ! this method does not give identical results as the above method
224 : ! when two dimensional decomposition is used
225 :
226 0 : do j = jfirst,jlast
227 0 : do i = 1,im
228 0 : dps_in(i,j) = sum( dpi_in(i,j,kfirst:klast) )
229 : end do
230 : end do
231 :
232 : #if ( defined SPMD )
233 0 : if (grid%twod_decomp .eq. 1) then
234 0 : call parcollective( grid%comm_z, SUMOP, im, jlast-jfirst+1, dps_in )
235 : endif
236 : #endif
237 :
238 : endif ! method1
239 :
240 : !-----------------------------------------------------------------------
241 : ! ... modify (fix) mass fluxes
242 : !-----------------------------------------------------------------------
243 32256 : call do_press_fix_llnl( grid, dps, dps_in, mfx, mfy )
244 :
245 : !-----------------------------------------------------------------------
246 : ! ... modified mass flux divergence
247 : !-----------------------------------------------------------------------
248 32256 : call calc_divergence( grid, mfx, mfy, dpi_c )
249 :
250 : !-----------------------------------------------------------------------
251 : ! ... differential mass flux divergence
252 : !-----------------------------------------------------------------------
253 182784 : do k = kfirst,klast
254 634368 : do j = jfirst,jlast
255 130658304 : dpi(:,j,k) = dpi_c(:,j,k) - dpi_in(:,j,k)
256 : end do
257 : end do
258 :
259 32256 : if (grid%twod_decomp .eq. 1) then
260 : #if defined (SPMD)
261 : call mp_sendirr( grid%commxy, grid%ijk_yz_to_xy%SendDesc, &
262 : grid%ijk_yz_to_xy%RecvDesc, dpi, dpixy, &
263 32256 : modc=grid%modc_dynrun )
264 : call mp_recvirr( grid%commxy, grid%ijk_yz_to_xy%SendDesc, &
265 : grid%ijk_yz_to_xy%RecvDesc, dpi, dpixy, &
266 32256 : modc=grid%modc_dynrun )
267 : #endif
268 : else
269 0 : dpixy(:,:,:) = dpi(:,:,:)
270 : endif
271 :
272 :
273 : !-----------------------------------------------------------------------
274 : ! ... modify pe
275 : !-----------------------------------------------------------------------
276 :
277 32256 : if (debug) then
278 0 : write(iulog,*) ' '
279 0 : write(iulog,*) 'adjust_press: max pe diff % @ nstep,ifirstxy,ilastxy,jfirstxy,jlastxy = ',&
280 0 : nstep,ifirstxy,ilastxy,jfirstxy,jlastxy
281 : endif
282 :
283 1838592 : do k = 1+1,km+1
284 1806336 : km1 = k - 1
285 :
286 1806336 : if (debug) then
287 0 : do j = jfirstxy,jlastxy
288 0 : do i = ifirstxy,ilastxy
289 0 : ps_diffxy(i,j) = sum( dpixy(i,j,1:km1) )/ pexy(i,k,j )
290 : end do
291 : end do
292 : endif
293 :
294 1806336 : if( nstep > nstep0 ) then
295 7225344 : do j = jfirstxy,jlastxy
296 137281536 : do i = ifirstxy,ilastxy
297 3842076672 : pexy(i,k,j) = pexy(i,k,j) + sum( dpixy(i,j,1:km1) )
298 : end do
299 : end do
300 : end if
301 1838592 : if (debug) then
302 :
303 0 : ndx(:) = maxloc( abs( ps_diffxy(:,:) ) )
304 :
305 0 : ndx(1) = ndx(1) + ifirstxy - 1
306 0 : ndx(2) = ndx(2) + jfirstxy - 1
307 :
308 : write(iulog,'("pfixer press change error (% error,press adjmnt,new pe)",1x,3i5,1p,3g15.7)') &
309 0 : k,ndx(:),D100_0*abs( ps_diffxy(ndx(1),ndx(2)) ), &
310 0 : dpixy(ndx(1),ndx(2),km1),pexy(ndx(1),k,ndx(2))
311 :
312 : endif
313 : end do
314 :
315 32256 : if (debug) then
316 0 : write(iulog,*) ' '
317 0 : write(iulog,*) 'adjust_press: max mass flux error @ nstep,jfirst,jlast,kfirst,klast = ',&
318 0 : nstep,jfirst,jlast,kfirst,klast
319 :
320 0 : do k = kfirst,klast
321 :
322 0 : do j=jfirst,jlast
323 0 : do i=1,im
324 0 : emfx(i,j) = ( mfx(i,j,k)-dmfx(i,j,k) )
325 : enddo
326 : enddo
327 :
328 0 : ndx(:) = maxloc( abs( emfx(:,:) ) )
329 0 : ndx(2) = ndx(2) + jfirst - 1
330 :
331 : write(iulog,'("pfixer max x flux error (diff,fixed,orig) ",1x,3i5,1p,3g15.7)') &
332 0 : k,ndx(:), emfx( ndx(1),ndx(2) ) , &
333 0 : mfx(ndx(1),ndx(2),k), dmfx(ndx(1),ndx(2),k)
334 :
335 0 : do j=jfirst,jlast+1
336 0 : do i=1,im
337 0 : emfy(i,j) = ( mfy(i,j,k)-dmfy(i,j,k) )
338 : enddo
339 : enddo
340 :
341 0 : ndx(:) = maxloc( abs( emfy(:,:) ) )
342 0 : ndx(2) = ndx(2) + jfirst - 1
343 :
344 : write(iulog,'("pfixer max y flux error (diff,fixed,orig) ",1x,3i5,1p,3g15.7)') &
345 0 : k,ndx(:), emfy( ndx(1),ndx(2) ) , &
346 0 : mfy(ndx(1),ndx(2),k), dmfy(ndx(1),ndx(2),k)
347 :
348 : enddo
349 : endif
350 :
351 32256 : end subroutine adjust_press
352 :
353 : !-----------------------------------------------------------------------
354 : ! ... calculate horizontal mass flux divergence
355 : !-----------------------------------------------------------------------
356 64512 : subroutine calc_divergence( grid, mfx, mfy, dpi )
357 :
358 : implicit none
359 :
360 : !-----------------------------------------------------------------------
361 : ! ... dummy arguments
362 : !-----------------------------------------------------------------------
363 : type (T_FVDYCORE_GRID), intent(in) :: grid
364 : real(r8), intent(in) :: mfx(grid%im,grid%jfirst:grid%jlast, &
365 : grid%kfirst:grid%klast) ! zonal mass flux
366 : real(r8), intent(in) :: mfy(grid%im,grid%jfirst:grid%jlast+1, &
367 : grid%kfirst:grid%klast) ! meridional mass flux
368 : real(r8), intent(inout) :: dpi(grid%im,grid%jfirst:grid%jlast, &
369 : grid%kfirst:grid%klast) ! horizontal mass flux divergence
370 :
371 : !-----------------------------------------------------------------------
372 : ! ... local variables
373 : !-----------------------------------------------------------------------
374 : integer :: i, j, k, js2g0, jn2g0
375 : real(r8) :: sum1
376 : integer :: im, jm, km, jfirst, jlast, kfirst, klast
377 :
378 64512 : im = grid%im
379 64512 : jm = grid%jm
380 64512 : km = grid%km
381 64512 : jfirst= grid%jfirst
382 64512 : jlast = grid%jlast
383 64512 : kfirst= grid%kfirst
384 64512 : klast = grid%klast
385 :
386 64512 : js2g0 = max( 2,jfirst )
387 64512 : jn2g0 = min( jm-1,jlast )
388 :
389 : !$omp parallel do private( j, k, sum1 )
390 365568 : do k = kfirst,klast
391 : !-----------------------------------------------------------------------
392 : ! ... north-south component
393 : !-----------------------------------------------------------------------
394 1194816 : do j = js2g0,jn2g0
395 258597696 : dpi(:,j,k) = (mfy(:,j,k) - mfy(:,j+1,k)) * grid%acosp(j)
396 : end do
397 : !-----------------------------------------------------------------------
398 : ! ... east-west component
399 : !-----------------------------------------------------------------------
400 1194816 : do j = js2g0,jn2g0
401 257402880 : dpi(:im-1,j,k) = dpi(:im-1,j,k) + mfx(:im-1,j,k) - mfx(2:im,j,k)
402 1194816 : dpi(im,j,k) = dpi(im,j,k) + mfx(im,j,k) - mfx(1,j,k)
403 : end do
404 : !-----------------------------------------------------------------------
405 : ! ... poles
406 : !-----------------------------------------------------------------------
407 301056 : if( jfirst == 1 ) then
408 1359456 : sum1 = -sum( mfy(:,2,k) )*grid%rcap
409 1359456 : dpi(:,1,k) = sum1
410 : end if
411 365568 : if( jlast == jm ) then
412 1359456 : sum1 = sum( mfy(:,jm,k) ) * grid%rcap
413 1359456 : dpi(:,jm,k) = sum1
414 : end if
415 : end do
416 : !$omp end parallel do
417 :
418 32256 : end subroutine calc_divergence
419 :
420 : !-----------------------------------------------------------------------
421 : ! ... fix the mass fluxes to match the met field pressure tendency
422 : ! See: http://asd.llnl.gov/pfix/index.html
423 : !-----------------------------------------------------------------------
424 32256 : subroutine do_press_fix_llnl( grid, dps, dps_ctm, mfx, mfy )
425 :
426 : use commap, only : gw => w
427 : #ifdef SPMD
428 : use mpishorthand, only : mpicom, mpi_double_precision, mpi_success
429 : use spmd_dyn, only : compute_gsfactors
430 : #endif
431 : use spmd_utils, only : npes
432 : use hycoef, only : hybd
433 : implicit none
434 :
435 : !-----------------------------------------------------------------------
436 : ! ... dummy arguments
437 : !-----------------------------------------------------------------------
438 : type (T_FVDYCORE_GRID), intent(in) :: grid
439 : ! surface pressure change from met fields
440 : real(r8), intent(in) :: dps(grid%im,grid%jfirst:grid%jlast)
441 : ! vert. sum of dpi from original mass fluxes
442 : real(r8), intent(in) :: dps_ctm(grid%im,grid%jfirst:grid%jlast)
443 : ! zonal mass flux
444 : real(r8), intent(inout) :: mfx(grid%im,grid%jfirst:grid%jlast,grid%kfirst:grid%klast)
445 : ! meridional mass flux
446 : real(r8), intent(inout) :: mfy(grid%im,grid%jfirst:grid%jlast+1,grid%kfirst:grid%klast)
447 :
448 : !-----------------------------------------------------------------------
449 : ! ... local variables
450 : !-----------------------------------------------------------------------
451 : integer :: i, j, jglob, k, astat, ierr
452 : integer :: jn2g0, js2g0, jn2g1
453 : integer :: cnt
454 : #ifdef SPMD
455 64512 : integer :: numrecv(0:npes-1)
456 64512 : integer :: displs(0:npes-1)
457 : #endif
458 : real(r8) :: dpress_g ! global pressure error
459 : real(r8) :: fxmean, factor
460 64512 : real(r8) :: ddps(grid%im,grid%jfirst:grid%jlast) ! surface pressure change error
461 64512 : real(r8) :: dpresslat(grid%jm)
462 64512 : real(r8) :: mmfd(grid%jm)
463 64512 : real(r8) :: mmf(grid%jm+1)
464 64512 : real(r8) :: fxintegral(grid%im+1)
465 64512 : real(r8) :: xcolmass_fix(grid%im,grid%jfirst:grid%jlast)
466 64512 : real(r8) :: temp(grid%jfirst:grid%jlast)
467 :
468 : integer :: im, jm, km, jfirst, jlast, kfirst, klast
469 :
470 32256 : im = grid%im
471 32256 : jm = grid%jm
472 32256 : km = grid%km
473 32256 : jfirst= grid%jfirst
474 32256 : jlast = grid%jlast
475 32256 : kfirst= grid%kfirst
476 32256 : klast = grid%klast
477 :
478 32256 : js2g0 = max( 2,jfirst )
479 32256 : jn2g0 = min( jm-1,jlast )
480 32256 : jn2g1 = min( jm-1,jlast+1 )
481 :
482 129024 : do j = jfirst,jlast
483 27998208 : ddps(:,j) = dps(:,j) - dps_ctm(:,j)
484 : end do
485 32256 : factor = D0_5/im
486 129024 : do j = jfirst,jlast
487 27998208 : dpresslat(j) = sum( ddps(:,j) ) * gw(j) * factor
488 : end do
489 :
490 : #ifdef SPMD
491 32256 : call compute_gsfactors( 1, cnt, numrecv, displs )
492 161280 : temp(jfirst:jlast) = dpresslat(jfirst:jlast)
493 : call mpi_allgatherv( temp(jfirst:jlast), cnt, mpi_double_precision, &
494 32256 : dpresslat, numrecv, displs, mpi_double_precision, mpicom, ierr )
495 32256 : if( ierr /= mpi_success ) then
496 0 : write(iulog,*) 'do_press_fix_llnl: mpi_allgatherv failed; error code = ',ierr
497 0 : call endrun
498 : end if
499 : #endif
500 :
501 6225408 : dpress_g = sum( dpresslat(:) )
502 32256 : if( grid%iam == 0 ) then
503 42 : write(iulog,*) 'do_press_fix_llnl: dpress_g = ',dpress_g
504 : end if
505 :
506 : !-----------------------------------------------------------------------
507 : ! calculate mean meridional flux divergence (df/dy).
508 : ! note that mmfd is actually the zonal mean pressure change,
509 : ! which is related to df/dy by geometrical factors.
510 : !-----------------------------------------------------------------------
511 : !-----------------------------------------------------------------------
512 : ! ... handle non-pole regions.
513 : !-----------------------------------------------------------------------
514 32256 : factor = D1_0/im
515 129024 : do j = jfirst,jlast
516 27998208 : mmfd(j) = dpress_g - sum( ddps(:,j) ) * factor
517 : end do
518 :
519 : #ifdef SPMD
520 32256 : cnt = jlast - jfirst + 1
521 129024 : temp(jfirst:jlast) = mmfd(jfirst:jlast)
522 : call mpi_allgatherv( temp(jfirst:jlast), cnt, mpi_double_precision, &
523 32256 : mmfd, numrecv, displs, mpi_double_precision, mpicom, ierr )
524 32256 : if( ierr /= mpi_success ) then
525 0 : write(iulog,*) 'do_press_fix_llnl: mpi_allgatherv failed; error code = ',ierr
526 0 : call endrun
527 : end if
528 : #endif
529 :
530 : !-----------------------------------------------------------------------
531 : ! calculate mean meridional fluxes (cosp*fy).
532 : ! nb: this calculation is being done for global lats, i.e., (1,jm)
533 : !-----------------------------------------------------------------------
534 32256 : mmf(2) = mmfd(1) / (grid%rcap*im)
535 6160896 : do j = 2,jm-1
536 6160896 : mmf(j+1) = mmf(j) + mmfd(j) * grid%cosp(j)
537 : end do
538 :
539 : !-----------------------------------------------------------------------
540 : ! fix latitude bands.
541 : ! note that we do not need to worry about geometry here because
542 : ! all boxes in a latitude band are identical.
543 : ! note also that fxintegral(im+1) should equal fxintegral(1),
544 : ! i.e., zero.
545 : !-----------------------------------------------------------------------
546 : !$omp parallel do private( i, j, k, fxmean, fxintegral )
547 128016 : do j = js2g0,jn2g0
548 95760 : fxintegral(1) = D0_0
549 27674640 : do i = 1,im
550 27674640 : fxintegral(i+1) = fxintegral(i) - (ddps(i,j) - dpress_g) - mmfd(j)
551 : end do
552 95760 : fxintegral(1) = fxintegral(im+1)
553 27770400 : fxmean = sum( fxintegral(:im) ) * factor
554 27706896 : xcolmass_fix(:im,j) = fxintegral(:im) - fxmean
555 : end do
556 : !$omp end parallel do
557 :
558 : !-----------------------------------------------------------------------
559 : ! ... distribute colmass_fix in vertical
560 : !-----------------------------------------------------------------------
561 : !$omp parallel do private( j, k )
562 182784 : do k = kfirst,klast
563 597408 : do j = js2g0,jn2g0
564 129298848 : mfx(:,j,k) = mfx(:,j,k) + met_rlx(k) * xcolmass_fix(:,j) * hybd(k)
565 : end do
566 745584 : do j = js2g0,jn2g1
567 172121712 : mfy(:,j,k) = mfy(:,j,k) + met_rlx(k) * mmf(j) * hybd(k)
568 : end do
569 182784 : if( jlast == jm ) then
570 679728 : mfy(:,jm,k) = mfy(:,jm,k) + met_rlx(k) * mmf(jm) * hybd(k)
571 : end if
572 : end do
573 : !$omp end parallel do
574 :
575 32256 : end subroutine do_press_fix_llnl
576 :
577 : end module pfixer
|
