Line data Source code
1 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
2 : !! Begin GPU remap module !!
3 : !! by Rick Archibald, 2010 !!
4 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
5 : module vertremap_mod
6 :
7 : !**************************************************************************************
8 : !
9 : ! Purpose:
10 : ! Construct sub-grid-scale polynomials using piecewise spline method with
11 : ! monotone filters.
12 : !
13 : ! References: PCM - Zerroukat et al., Q.J.R. Meteorol. Soc., 2005. (ZWS2005QJR)
14 : ! PSM - Zerroukat et al., Int. J. Numer. Meth. Fluids, 2005. (ZWS2005IJMF)
15 : !
16 : !**************************************************************************************
17 :
18 : use shr_kind_mod, only: r8=>shr_kind_r8
19 : use dimensions_mod, only: np,nlev,qsize,nlevp,npsq,nc
20 : use element_mod, only: element_t
21 : use fvm_control_volume_mod, only: fvm_struct
22 : use perf_mod, only: t_startf, t_stopf ! _EXTERNAL
23 : use parallel_mod, only: parallel_t
24 : use cam_abortutils, only: endrun
25 :
26 : implicit none
27 :
28 : public remap1 ! remap any field, splines, monotone
29 : public remap1_nofilter ! remap any field, splines, no filter
30 : ! todo: tweak interface to match remap1 above, rename remap1_ppm:
31 : public remap_q_ppm ! remap state%Q, PPM, monotone
32 :
33 : contains
34 :
35 : !=======================================================================================================!
36 :
37 25974000 : subroutine remap1(Qdp,nx,qstart,qstop,qsize,dp1,dp2,ptop,identifier,Qdp_mass,kord)
38 : use fv_mapz, only: map1_ppm
39 : ! remap 1 field
40 : ! input: Qdp field to be remapped (NOTE: MASS, not MIXING RATIO)
41 : ! dp1 layer thickness (source)
42 : ! dp2 layer thickness (target)
43 : !
44 : ! output: remaped Qdp, conserving mass, monotone on Q=Qdp/dp
45 : !
46 : integer, intent(in) :: nx,qstart,qstop,qsize
47 : real (kind=r8), intent(inout) :: Qdp(nx,nx,nlev,qsize)
48 : real (kind=r8), intent(in) :: dp1(nx,nx,nlev),dp2(nx,nx,nlev)
49 : integer, intent(in) :: identifier !0: tracers, 1: T, -1: u,v
50 : real (kind=r8), intent(in) :: ptop
51 : logical, intent(in) :: Qdp_mass
52 : integer, intent(in) :: kord(qsize)
53 : ! ========================
54 : ! Local Variables
55 : ! ========================
56 51948000 : real (kind=r8) :: pe1(nx,nlev+1),pe2(nx,nlev+1),inv_dp(nx,nx,nlev),dp2_local(nx,nlev)
57 51948000 : real (kind=r8) :: tmp(nx,nlev), gz(nx)
58 : integer :: i,j,k,itrac
59 : logical :: logp
60 25974000 : integer :: kord_local(qsize)
61 :
62 129870000 : kord_local = kord
63 :
64 88311600 : if (any(kord(:) >= 0)) then
65 5194800 : if (.not.qdp_mass) then
66 62337600 : do itrac=1,qsize
67 62337600 : if (kord(itrac) >= 0) then
68 31428540000 : Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)*dp1(:,:,:)
69 : end if
70 : end do
71 : end if
72 5194800 : call remap_Q_ppm(qdp,nx,qstart,qstop,qsize,dp1,dp2,kord)
73 5194800 : if (.not.qdp_mass) then
74 62337600 : do itrac=1,qsize
75 62337600 : if (kord(itrac) >= 0) then
76 31428540000 : Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)/dp2(:,:,:)
77 : end if
78 : end do
79 : end if
80 : endif
81 25974000 : if (any(kord(:)<0)) then
82 : !
83 : ! check if remapping over p or log(p)
84 : !
85 : ! can not mix and match here (all kord's must >-20 or <=-20)
86 : !
87 25974000 : if (any(kord(:)>-20)) then
88 129870000 : kord_local = abs(kord)
89 : logp = .false.
90 : else
91 0 : kord_local = abs(kord/10)
92 0 : if (identifier==1) then
93 : logp = .true.
94 : else
95 0 : logp = .false.
96 : end if
97 : end if
98 : !
99 : ! modified FV3 vertical remapping
100 : !
101 25974000 : if (qdp_mass) then
102 20301278400 : inv_dp = 1.0_r8/dp1
103 46753200 : do itrac=1,qsize
104 46753200 : if (kord(itrac)<0) then
105 71054474400 : Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)*inv_dp(:,:,:)
106 : end if
107 : end do
108 : end if
109 25974000 : if (logp) then
110 0 : do j=1,nx
111 0 : pe1(:,1) = ptop
112 0 : pe2(:,1) = ptop
113 0 : do k=1,nlev
114 0 : do i=1,nx
115 0 : pe1(i,k+1) = pe1(i,k)+dp1(i,j,k)
116 0 : pe2(i,k+1) = pe2(i,k)+dp2(i,j,k)
117 : end do
118 : end do
119 0 : pe1(:,nlev+1) = pe2(:,nlev+1)
120 0 : do k=1,nlev+1
121 0 : do i=1,nx
122 0 : pe1(i,k) = log(pe1(i,k))
123 0 : pe2(i,k) = log(pe2(i,k))
124 : end do
125 : end do
126 :
127 0 : do itrac=1,qsize
128 0 : if (kord(itrac)<0) then
129 : call map1_ppm( nlev, pe1(:,:), Qdp(:,:,:,itrac), gz, &
130 : nlev, pe2(:,:), Qdp(:,:,:,itrac), &
131 0 : 1, nx, j, 1, nx, 1, nx, identifier, kord_local(itrac))
132 : end if
133 : end do
134 : ! call mapn_tracer(qsize, nlev, pe1, pe2, Qdp, dp2_local, kord, j, &
135 : ! 1, nx, 1, nx, 1, nx, 0.0_r8, fill)
136 : end do
137 : else
138 124675200 : do j=1,nx
139 477921600 : pe1(:,1) = ptop
140 477921600 : pe2(:,1) = ptop
141 9277912800 : do k=1,nlev
142 44545410000 : do i=1,nx
143 35267497200 : pe1(i,k+1) = pe1(i,k)+dp1(i,j,k)
144 44446708800 : pe2(i,k+1) = pe2(i,k)+dp2(i,j,k)
145 : end do
146 : end do
147 477921600 : pe1(:,nlev+1) = pe2(:,nlev+1)
148 483116400 : do itrac=1,qsize
149 457142400 : if (kord(itrac)<0) then
150 : call map1_ppm( nlev, pe1(:,:), Qdp(:,:,:,itrac), gz, &!phl
151 : nlev, pe2(:,:), Qdp(:,:,:,itrac), &
152 280519200 : 1, nx, j, 1, nx, 1, nx, identifier, kord_local(itrac))
153 : end if
154 : end do
155 : ! call mapn_tracer(qsize, nlev, pe1, pe2, Qdp, dp2_local, kord, j, &
156 : ! 1, nx, 1, nx, 1, nx, 0.0_r8, fill)
157 : end do
158 : end if
159 25974000 : if (qdp_mass) then
160 46753200 : do itrac=1,qsize
161 46753200 : if (kord(itrac)<0) then
162 71054474400 : Qdp(:,:,:,itrac) = Qdp(:,:,:,itrac)*dp2(:,:,:)
163 : end if
164 : end do
165 : end if
166 : end if
167 25974000 : end subroutine remap1
168 :
169 0 : subroutine remap1_nofilter(Qdp,nx,qsize,dp1,dp2)
170 : ! remap 1 field
171 : ! input: Qdp field to be remapped (NOTE: MASS, not MIXING RATIO)
172 : ! dp1 layer thickness (source)
173 : ! dp2 layer thickness (target)
174 : !
175 : ! output: remaped Qdp, conserving mass
176 : !
177 : implicit none
178 : integer, intent(in) :: nx,qsize
179 : real (kind=r8), intent(inout) :: Qdp(nx,nx,nlev,qsize)
180 : real (kind=r8), intent(in) :: dp1(nx,nx,nlev),dp2(nx,nx,nlev)
181 : ! ========================
182 : ! Local Variables
183 : ! ========================
184 :
185 : real (kind=r8), dimension(nlev+1) :: rhs,lower_diag,diag,upper_diag,q_diag,zgam,z1c,z2c,zv
186 : real (kind=r8), dimension(nlev) :: h,Qcol,za0,za1,za2,zarg,zhdp
187 : real (kind=r8) :: tmp_cal,zv1,zv2
188 : integer :: zkr(nlev+1),i,ilev,j,jk,k,q
189 : logical :: abort=.false.
190 : ! call t_startf('remap1_nofilter')
191 :
192 : #if (defined COLUMN_OPENMP)
193 : !$omp parallel do num_threads(tracer_num_threads) &
194 : !$omp private(q,i,j,z1c,z2c,zv,k,Qcol,zkr,ilev) &
195 : !$omp private(jk,zgam,zhdp,h,zarg,rhs,lower_diag,diag,upper_diag,q_diag,tmp_cal) &
196 : !$omp private(za0,za1,za2) &
197 : !$omp private(ip2,zv1,zv2)
198 : #endif
199 0 : do q=1,qsize
200 0 : do i=1,nx
201 0 : do j=1,nx
202 :
203 0 : z1c(1)=0 ! source grid
204 0 : z2c(1)=0 ! target grid
205 0 : do k=1,nlev
206 0 : z1c(k+1)=z1c(k)+dp1(i,j,k)
207 0 : z2c(k+1)=z2c(k)+dp2(i,j,k)
208 : enddo
209 :
210 0 : zv(1)=0
211 0 : do k=1,nlev
212 0 : Qcol(k)=Qdp(i,j,k,q)! *(z1c(k+1)-z1c(k)) input is mass
213 0 : zv(k+1) = zv(k)+Qcol(k)
214 : enddo
215 :
216 0 : if (ABS(z2c(nlev+1)-z1c(nlev+1)) >= 0.000001_r8) then
217 0 : write(6,*) 'SURFACE PRESSURE IMPLIED BY ADVECTION SCHEME'
218 0 : write(6,*) 'NOT CORRESPONDING TO SURFACE PRESSURE IN '
219 0 : write(6,*) 'DATA FOR MODEL LEVELS'
220 0 : write(6,*) 'PLEVMODEL=',z2c(nlev+1)
221 0 : write(6,*) 'PLEV =',z1c(nlev+1)
222 0 : write(6,*) 'DIFF =',z2c(nlev+1)-z1c(nlev+1)
223 0 : abort=.true.
224 : endif
225 :
226 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
227 : !! quadratic splies with UK met office monotonicity constraints !!
228 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
229 :
230 0 : zkr = 99
231 0 : ilev = 2
232 0 : zkr(1) = 1
233 0 : zkr(nlev+1) = nlev
234 0 : kloop: do k = 2,nlev
235 0 : do jk = ilev,nlev+1
236 0 : if (z1c(jk) >= z2c(k)) then
237 0 : ilev = jk
238 0 : zkr(k) = jk-1
239 0 : cycle kloop
240 : endif
241 : enddo
242 : enddo kloop
243 :
244 0 : zgam = (z2c(1:nlev+1)-z1c(zkr)) / (z1c(zkr+1)-z1c(zkr))
245 0 : zgam(1) = 0.0_r8
246 0 : zgam(nlev+1) = 1.0_r8
247 0 : zhdp = z1c(2:nlev+1)-z1c(1:nlev)
248 :
249 :
250 0 : h = 1/zhdp
251 0 : zarg = Qcol * h
252 0 : rhs = 0
253 0 : lower_diag = 0
254 0 : diag = 0
255 0 : upper_diag = 0
256 :
257 0 : rhs(1)=3*zarg(1)
258 0 : rhs(2:nlev) = 3*(zarg(2:nlev)*h(2:nlev) + zarg(1:nlev-1)*h(1:nlev-1))
259 0 : rhs(nlev+1)=3*zarg(nlev)
260 :
261 0 : lower_diag(1)=1
262 0 : lower_diag(2:nlev) = h(1:nlev-1)
263 0 : lower_diag(nlev+1)=1
264 :
265 0 : diag(1)=2
266 0 : diag(2:nlev) = 2*(h(2:nlev) + h(1:nlev-1))
267 0 : diag(nlev+1)=2
268 :
269 0 : upper_diag(1)=1
270 0 : upper_diag(2:nlev) = h(2:nlev)
271 0 : upper_diag(nlev+1)=0
272 :
273 0 : q_diag(1)=-upper_diag(1)/diag(1)
274 0 : rhs(1)= rhs(1)/diag(1)
275 :
276 0 : do k=2,nlev+1
277 0 : tmp_cal = 1/(diag(k)+lower_diag(k)*q_diag(k-1))
278 0 : q_diag(k) = -upper_diag(k)*tmp_cal
279 0 : rhs(k) = (rhs(k)-lower_diag(k)*rhs(k-1))*tmp_cal
280 : enddo
281 0 : do k=nlev,1,-1
282 0 : rhs(k)=rhs(k)+q_diag(k)*rhs(k+1)
283 : enddo
284 :
285 0 : za0 = rhs(1:nlev)
286 0 : za1 = -4*rhs(1:nlev) - 2*rhs(2:nlev+1) + 6*zarg
287 0 : za2 = 3*rhs(1:nlev) + 3*rhs(2:nlev+1) - 6*zarg
288 :
289 :
290 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
291 : !! start iteration from top to bottom of atmosphere !!
292 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
293 :
294 : zv1 = 0
295 0 : do k=1,nlev
296 0 : if (zgam(k+1)>1_r8) then
297 0 : WRITE(*,*) 'r not in [0:1]', zgam(k+1)
298 0 : abort=.true.
299 : endif
300 0 : zv2 = zv(zkr(k+1))+(za0(zkr(k+1))*zgam(k+1)+(za1(zkr(k+1))/2)*(zgam(k+1)**2)+ &
301 0 : (za2(zkr(k+1))/3)*(zgam(k+1)**3))*zhdp(zkr(k+1))
302 0 : Qdp(i,j,k,q) = (zv2 - zv1) ! / (z2c(k+1)-z2c(k) ) dont convert back to mixing ratio
303 0 : zv1 = zv2
304 : enddo
305 : enddo
306 : enddo
307 : enddo ! q loop
308 0 : if (abort) then
309 0 : call endrun('Bad levels in remap1_nofilter. usually CFL violatioin')
310 : end if
311 0 : end subroutine remap1_nofilter
312 :
313 : !=============================================================================!
314 :
315 : !This uses the exact same model and reference grids and data as remap_Q, but it interpolates
316 : !using PPM instead of splines.
317 5194800 : subroutine remap_Q_ppm(Qdp,nx,qstart,qstop,qsize,dp1,dp2,kord)
318 : ! remap 1 field
319 : ! input: Qdp field to be remapped (NOTE: MASS, not MIXING RATIO)
320 : ! dp1 layer thickness (source)
321 : ! dp2 layer thickness (target)
322 : !
323 : ! output: remaped Qdp, conserving mass
324 : !
325 : implicit none
326 : integer,intent(in) :: nx,qstart,qstop,qsize
327 : real (kind=r8), intent(inout) :: Qdp(nx,nx,nlev,qsize)
328 : real (kind=r8), intent(in) :: dp1(nx,nx,nlev),dp2(nx,nx,nlev)
329 : integer , intent(in) :: kord(qsize)
330 : ! Local Variables
331 : integer, parameter :: gs = 2 !Number of cells to place in the ghost region
332 : real(kind=r8), dimension( nlev+2 ) :: pio !Pressure at interfaces for old grid
333 : real(kind=r8), dimension( nlev+1 ) :: pin !Pressure at interfaces for new grid
334 : real(kind=r8), dimension( nlev+1 ) :: masso !Accumulate mass up to each interface
335 : real(kind=r8), dimension( 1-gs:nlev+gs) :: ao !Tracer value on old grid
336 : real(kind=r8), dimension( 1-gs:nlev+gs) :: dpo !change in pressure over a cell for old grid
337 : real(kind=r8), dimension( 1-gs:nlev+gs) :: dpn !change in pressure over a cell for old grid
338 : real(kind=r8), dimension(3, nlev ) :: coefs !PPM coefficients within each cell
339 : real(kind=r8), dimension( nlev ) :: z1, z2
340 : real(kind=r8) :: ppmdx(10,0:nlev+1) !grid spacings
341 : real(kind=r8) :: massn1, massn2, ext(2)
342 : integer :: i, j, k, q, kk, kid(nlev)
343 :
344 20779200 : do j = 1 , nx
345 67532400 : do i = 1 , nx
346 :
347 46753200 : pin(1)=0
348 46753200 : pio(1)=0
349 4394800800 : do k=1,nlev
350 4348047600 : dpn(k)=dp2(i,j,k)
351 4348047600 : dpo(k)=dp1(i,j,k)
352 4348047600 : pin(k+1)=pin(k)+dpn(k)
353 4394800800 : pio(k+1)=pio(k)+dpo(k)
354 : enddo
355 :
356 :
357 :
358 46753200 : pio(nlev+2) = pio(nlev+1) + 1._r8 !This is here to allow an entire block of k threads to run in the remapping phase.
359 : !It makes sure there's an old interface value below the domain that is larger.
360 46753200 : pin(nlev+1) = pio(nlev+1) !The total mass in a column does not change.
361 : !Therefore, the pressure of that mass cannot either.
362 : !Fill in the ghost regions with mirrored values. if vert_remap_q_alg is defined, this is of no consequence.
363 140259600 : do k = 1 , gs
364 93506400 : dpo(1 -k) = dpo( k)
365 140259600 : dpo(nlev+k) = dpo(nlev+1-k)
366 : enddo
367 :
368 : !Compute remapping intervals once for all tracers. Find the old grid cell index in which the
369 : !k-th new cell interface resides. Then integrate from the bottom of that old cell to the new
370 : !interface location. In practice, the grid never deforms past one cell, so the search can be
371 : !simplified by this. Also, the interval of integration is usually of magnitude close to zero
372 : !or close to dpo because of minimial deformation.
373 : !Numerous tests confirmed that the bottom and top of the grids match to machine precision, so
374 : !I set them equal to each other.
375 4394800800 : do k = 1 , nlev
376 4348047600 : kk = k !Keep from an order n^2 search operation by assuming the old cell index is close.
377 : !Find the index of the old grid cell in which this new cell's bottom interface resides.
378 : ! do while ( pio(kk) <= pin(k+1) )
379 : ! kk = kk + 1
380 : ! if(kk==nlev+2) exit
381 : ! enddo
382 : ! kk = kk - 1 !kk is now the cell index we're integrating over.
383 :
384 4348047600 : if (pio(kk) <= pin(k+1)) then
385 10744273667 : do while ( pio(kk) <= pin(k+1) )
386 6396229663 : kk = kk + 1
387 : enddo
388 4348044004 : kk = kk - 1 !kk is now the cell index we're integrating over.
389 : else
390 : call binary_search(pio, pin(k+1), kk)
391 : end if
392 4348047600 : if (kk == nlev+1) kk = nlev !This is to keep the indices in bounds.
393 : !Top bounds match anyway, so doesn't matter what coefficients are used
394 4348047600 : kid(k) = kk !Save for reuse
395 4348047600 : z1(k) = -0.5_R8 !This remapping assumes we're starting from the left interface of an old grid cell
396 : !In fact, we're usually integrating very little or almost all of the cell in question
397 4394800800 : z2(k) = ( pin(k+1) - ( pio(kk) + pio(kk+1) ) * 0.5_r8 ) / dpo(kk) !PPM interpolants are normalized to an independent
398 : !coordinate domain [-0.5,0.5].
399 : enddo
400 :
401 : !This turned out a big optimization, remembering that only parts of the PPM algorithm depends on the data, namely the
402 : !limiting. So anything that depends only on the grid is pre-computed outside the tracer loop.
403 46753200 : ppmdx(:,:) = compute_ppm_grids( dpo)
404 :
405 : !From here, we loop over tracers for only those portions which depend on tracer data, which includes PPM limiting and
406 : !mass accumulation
407 576622800 : do q = qstart, qstop
408 561038400 : if (kord(q) >= 0) then
409 : !Accumulate the old mass up to old grid cell interface locations to simplify integration
410 : !during remapping. Also, divide out the grid spacing so we're working with actual tracer
411 : !values and can conserve mass. The option for ifndef ZEROHORZ I believe is there to ensure
412 : !tracer consistency for an initially uniform field. I copied it from the old remap routine.
413 233766000 : masso(1) = 0._r8
414 :
415 21974004000 : do k = 1 , nlev
416 21740238000 : ao(k) = Qdp(i,j,k,q)
417 21740238000 : masso(k+1) = masso(k) + ao(k) !Accumulate the old mass. This will simplify the remapping
418 21974004000 : ao(k) = ao(k) / dpo(k) !Divide out the old grid spacing because we want the tracer mixing ratio, not mass.
419 : enddo
420 : !Fill in ghost values. Ignored if kord == 2
421 233766000 : if (kord(q) == 10) then
422 22207770000 : ext(1) = minval(ao(1:nlev))
423 22207770000 : ext(2) = maxval(ao(1:nlev))
424 233766000 : call linextrap(dpo(2), dpo(1), dpo(0), dpo(-1), ao(2), ao(1), ao(0), ao(-1), ext(1), ext(2))
425 : call linextrap(dpo(nlev-1), dpo(nlev), dpo(nlev+1), dpo(nlev+2), &
426 233766000 : ao(nlev-1), ao(nlev), ao(nlev+1), ao(nlev+2), ext(1), ext(2))
427 : else
428 0 : do k = 1 , gs
429 0 : ao(1 -k) = ao( k)
430 0 : ao(nlev+k) = ao(nlev+1-k)
431 : enddo
432 : end if
433 : !Compute monotonic and conservative PPM reconstruction over every cell
434 233766000 : coefs(:,:) = compute_ppm( ao , ppmdx, kord(q) )
435 : !Compute tracer values on the new grid by integrating from the old cell bottom to the new
436 : !cell interface to form a new grid mass accumulation. Taking the difference between
437 : !accumulation at successive interfaces gives the mass inside each cell. Since Qdp is
438 : !supposed to hold the full mass this needs no normalization.
439 233766000 : massn1 = 0._r8
440 21974004000 : do k = 1 , nlev
441 21740238000 : kk = kid(k)
442 21740238000 : massn2 = masso(kk) + integrate_parabola( coefs(:,kk) , z1(k) , z2(k) ) * dpo(kk)
443 21740238000 : Qdp(i,j,k,q) = massn2 - massn1
444 21974004000 : massn1 = massn2
445 : enddo
446 : end if
447 : enddo
448 : enddo
449 : enddo
450 : ! call t_stopf('remap_Q_ppm')
451 5194800 : end subroutine remap_Q_ppm
452 :
453 : ! Find k such that pio(k) <= pivot < pio(k+1). Provide a reasonable input
454 : ! value for k.
455 3596 : subroutine binary_search(pio, pivot, k)
456 : real(kind=r8), intent(in) :: pio(nlev+2), pivot
457 : integer, intent(inout) :: k
458 : integer :: lo, hi, mid
459 :
460 3596 : if (pio(k) > pivot) then
461 : lo = 1
462 : hi = k
463 : else
464 0 : lo = k
465 0 : hi = nlev+2
466 : end if
467 25394 : do while (hi > lo + 1)
468 21798 : k = (lo + hi)/2
469 25394 : if (pio(k) > pivot) then
470 : hi = k
471 : else
472 21798 : lo = k
473 : end if
474 : end do
475 3596 : k = lo
476 3596 : end subroutine binary_search
477 : !=======================================================================================================!
478 :
479 : !This compute grid-based coefficients from Collela & Woodward 1984.
480 46753200 : function compute_ppm_grids( dx ) result(rslt)
481 : implicit none
482 : real(kind=r8), intent(in) :: dx(-1:nlev+2) !grid spacings
483 : real(kind=r8) :: rslt(10,0:nlev+1) !grid spacings
484 : integer :: j
485 :
486 : !Calculate grid-based coefficients for stage 1 of compute_ppm
487 4488307200 : do j = 0 , nlev+1
488 4441554000 : rslt( 1,j) = dx(j) / ( dx(j-1) + dx(j) + dx(j+1) )
489 4441554000 : rslt( 2,j) = ( 2._r8*dx(j-1) + dx(j) ) / ( dx(j+1) + dx(j) )
490 4488307200 : rslt( 3,j) = ( dx(j) + 2._r8*dx(j+1) ) / ( dx(j-1) + dx(j) )
491 : enddo
492 :
493 : !Caculate grid-based coefficients for stage 2 of compute_ppm
494 4441554000 : do j = 0 , nlev
495 4394800800 : rslt( 4,j) = dx(j) / ( dx(j) + dx(j+1) )
496 26368804800 : rslt( 5,j) = 1._r8 / sum( dx(j-1:j+2) )
497 4394800800 : rslt( 6,j) = ( 2._r8 * dx(j+1) * dx(j) ) / ( dx(j) + dx(j+1 ) )
498 4394800800 : rslt( 7,j) = ( dx(j-1) + dx(j ) ) / ( 2._r8 * dx(j ) + dx(j+1) )
499 4394800800 : rslt( 8,j) = ( dx(j+2) + dx(j+1) ) / ( 2._r8 * dx(j+1) + dx(j ) )
500 4394800800 : rslt( 9,j) = dx(j ) * ( dx(j-1) + dx(j ) ) / ( 2._r8*dx(j ) + dx(j+1) )
501 4441554000 : rslt(10,j) = dx(j+1) * ( dx(j+1) + dx(j+2) ) / ( dx(j ) + 2._r8*dx(j+1) )
502 : enddo
503 46753200 : end function compute_ppm_grids
504 :
505 :
506 : !=======================================================================================================!
507 :
508 :
509 :
510 : !This computes a limited parabolic interpolant using a net 5-cell stencil, but the stages of computation are broken up into 3 stages
511 233766000 : function compute_ppm( a , dx , kord) result(coefs)
512 : implicit none
513 : real(kind=r8), intent(in) :: a ( -1:nlev+2) !Cell-mean values
514 : real(kind=r8), intent(in) :: dx (10, 0:nlev+1) !grid spacings
515 : integer, intent(in) :: kord
516 : real(kind=r8) :: coefs(0:2, nlev ) !PPM coefficients (for parabola)
517 : real(kind=r8) :: ai (0:nlev ) !fourth-order accurate, then limited interface values
518 : real(kind=r8) :: dma(0:nlev+1) !An expression from Collela's '84 publication
519 : real(kind=r8) :: da !Ditto
520 : ! Hold expressions based on the grid (which are cumbersome).
521 : real(kind=r8) :: al, ar !Left and right interface values for cell-local limiting
522 : integer :: j
523 : integer :: indB, indE
524 :
525 : ! Stage 1: Compute dma for each cell, allowing a 1-cell ghost stencil below and above the domain
526 22441536000 : do j = 0 , nlev+1
527 22207770000 : da = dx(1,j) * ( dx(2,j) * ( a(j+1) - a(j) ) + dx(3,j) * ( a(j) - a(j-1) ) )
528 >11103*10^7 : dma(j) = minval( (/ abs(da) , 2._r8 * abs( a(j) - a(j-1) ) , 2._r8 * abs( a(j+1) - a(j) ) /) ) * sign(1._r8,da)
529 22441536000 : if ( ( a(j+1) - a(j) ) * ( a(j) - a(j-1) ) <= 0._r8 ) dma(j) = 0._r8
530 : enddo
531 :
532 : ! Stage 2: Compute ai for each cell interface in the physical domain (dimension nlev+1)
533 22207770000 : do j = 0 , nlev
534 43948008000 : ai(j) = a(j) + dx(4,j) * ( a(j+1) - a(j) ) + dx(5,j) * ( dx(6,j) * ( dx(7,j) - dx(8,j) ) &
535 66155778000 : * ( a(j+1) - a(j) ) - dx(9,j) * dma(j+1) + dx(10,j) * dma(j) )
536 : enddo
537 :
538 : ! Stage 3: Compute limited PPM interpolant over each cell in the physical domain
539 : ! (dimension nlev) using ai on either side and ao within the cell.
540 21974004000 : do j = 1 , nlev
541 21740238000 : al = ai(j-1)
542 21740238000 : ar = ai(j )
543 21740238000 : if ( (ar - a(j)) * (a(j) - al) <= 0._r8 ) then
544 6319824757 : al = a(j)
545 6319824757 : ar = a(j)
546 : endif
547 21740238000 : if ( (ar - al) * (a(j) - (al + ar)/2._r8) > (ar - al)**2/6._r8 ) al = 3._r8*a(j) - 2._r8 * ar
548 21740238000 : if ( (ar - al) * (a(j) - (al + ar)/2._r8) < -(ar - al)**2/6._r8 ) ar = 3._r8*a(j) - 2._r8 * al
549 : !Computed these coefficients from the edge values and cell mean in Maple. Assumes normalized coordinates: xi=(x-x0)/dx
550 21740238000 : coefs(0,j) = 1.5_r8 * a(j) - ( al + ar ) / 4._r8
551 21740238000 : coefs(1,j) = ar - al
552 21974004000 : coefs(2,j) = 3._r8 * (-2._r8 * a(j) + ( al + ar ))
553 : enddo
554 :
555 : !If kord == 2, use piecewise constant in the boundaries, and don't use ghost cells.
556 233766000 : if (kord == 2) then
557 0 : coefs(0,1:2) = a(1:2)
558 0 : coefs(1:2,1:2) = 0._r8
559 0 : coefs(0,nlev-1:nlev) = a(nlev-1:nlev)
560 0 : coefs(1:2,nlev-1:nlev) = 0._r8
561 : endif
562 233766000 : end function compute_ppm
563 :
564 :
565 : !=======================================================================================================!
566 :
567 :
568 : !Simple function computes the definite integral of a parabola in normalized coordinates, xi=(x-x0)/dx,
569 : !given two bounds. Make sure this gets inlined during compilation.
570 21740238000 : function integrate_parabola( a , x1 , x2 ) result(mass)
571 : implicit none
572 : real(kind=r8), intent(in) :: a(0:2) !Coefficients of the parabola
573 : real(kind=r8), intent(in) :: x1 !lower domain bound for integration
574 : real(kind=r8), intent(in) :: x2 !upper domain bound for integration
575 : real(kind=r8) :: mass
576 21740238000 : mass = a(0) * (x2 - x1) + a(1) * (x2 ** 2 - x1 ** 2) / 0.2D1 + a(2) * (x2 ** 3 - x1 ** 3) / 0.3D1
577 21740238000 : end function integrate_parabola
578 :
579 :
580 : !=============================================================================================!
581 467532000 : subroutine linextrap(dx1,dx2,dx3,dx4,y1,y2,y3,y4,lo,hi)
582 : real(kind=r8), intent(in) :: dx1,dx2,dx3,dx4,y1,y2,lo,hi
583 : real(kind=r8), intent(out) :: y3,y4
584 :
585 : real(kind=r8), parameter :: half = 0.5_r8
586 :
587 : real(kind=r8) :: x1,x2,x3,x4,a
588 :
589 467532000 : x1 = half*dx1
590 467532000 : x2 = x1 + half*(dx1 + dx2)
591 467532000 : x3 = x2 + half*(dx2 + dx3)
592 467532000 : x4 = x3 + half*(dx3 + dx4)
593 467532000 : a = (x3-x1)/(x2-x1)
594 467532000 : y3 = (1.0_r8-a)*y1 + a*y2
595 467532000 : a = (x4-x1)/(x2-x1)
596 467532000 : y4 = (1.0_r8-a)*y1 + a*y2
597 467532000 : y3 = max(lo, min(hi, y3))
598 467532000 : y4 = max(lo, min(hi, y4))
599 467532000 : end subroutine linextrap
600 : end module vertremap_mod
601 :
602 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
603 : !! End GPU remap module !!
604 : !! by Rick Archibald, 2010 !!
605 : !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
|
