1 2! KGEN-generated Fortran source file 3! 4! Filename : prim_advection_mod.F90 5! Generated at: 2015-02-24 15:34:48 6! KGEN version: 0.4.4 7 8 9 10 MODULE vertremap_mod 11 !************************************************************************************** 12 ! 13 ! Purpose: 14 ! Construct sub-grid-scale polynomials using piecewise spline method with 15 ! monotone filters. 16 ! 17 ! References: PCM - Zerroukat et al., Q.J.R. Meteorol. Soc., 2005. (ZWS2005QJR) 18 ! PSM - Zerroukat et al., Int. J. Numer. Meth. Fluids, 2005. (ZWS2005IJMF) 19 ! 20 !************************************************************************************** 21 USE kinds, ONLY: real_kind 22 USE dimensions_mod, ONLY: nlev 23 USE perf_mod, ONLY: t_startf 24 USE perf_mod, ONLY: t_stopf ! _EXTERNAL 25 INTEGER, PARAMETER :: kgen_dp = selected_real_kind(15, 307) 26 PUBLIC remap1 27 type, public :: check_t 28 logical :: Passed 29 integer :: numFatal 30 integer :: numTotal 31 integer :: numIdentical 32 integer :: numWarning 33 integer :: VerboseLevel 34 real(kind=kgen_dp) :: tolerance 35 end type check_t 36 ! remap any field, splines, monotone 37 ! remap any field, splines, no filter 38 ! todo: tweak interface to match remap1 above, rename remap1_ppm: 39 PUBLIC remap_q_ppm ! remap state%Q, PPM, monotone 40 CONTAINS 41 subroutine kgen_init_check(check,tolerance) 42 type(check_t), intent(inout) :: check 43 real(kind=kgen_dp), intent(in), optional :: tolerance 44 check%Passed = .TRUE. 45 check%numFatal = 0 46 check%numWarning = 0 47 check%numTotal = 0 48 check%numIdentical = 0 49 check%VerboseLevel = 1 50 if(present(tolerance)) then 51 check%tolerance = tolerance 52 else 53 check%tolerance = 1.E-14 54 endif 55 end subroutine kgen_init_check 56 subroutine kgen_print_check(kname, check) 57 character(len=*) :: kname 58 type(check_t), intent(in) :: check 59 write (*,*) 60 write (*,*) TRIM(kname),' KGENPrtCheck: Tolerance for normalized RMS: ',check%tolerance 61 write (*,*) TRIM(kname),' KGENPrtCheck: Number of variables checked: ',check%numTotal 62 write (*,*) TRIM(kname),' KGENPrtCheck: Number of Identical results: ',check%numIdentical 63 write (*,*) TRIM(kname),' KGENPrtCheck: Number of warnings detected: ',check%numWarning 64 write (*,*) TRIM(kname),' KGENPrtCheck: Number of fatal errors detected: ', check%numFatal 65 if (check%numFatal> 0) then 66 write(*,*) TRIM(kname),' KGENPrtCheck: verification FAILED' 67 else 68 write(*,*) TRIM(kname),' KGENPrtCheck: verification PASSED' 69 endif 70 end subroutine kgen_print_check 71 !=======================================================================================================! 72 !remap_calc_grids computes the vertical pressures and pressure differences for one vertical column for the reference grid 73 !and for the deformed Lagrangian grid. This was pulled out of each routine since it was a repeated task. 74 75 !=======================================================================================================! 76 77 SUBROUTINE remap1(nx, qsize, qdp, dp1, dp2, kgen_unit) 78 ! remap 1 field 79 ! input: Qdp field to be remapped (NOTE: MASS, not MIXING RATIO) 80 ! dp1 layer thickness (source) 81 ! dp2 layer thickness (target) 82 ! 83 ! output: remaped Qdp, conserving mass, monotone on Q=Qdp/dp 84 ! 85 IMPLICIT NONE 86 integer, intent(in) :: kgen_unit 87 88 ! read interface 89 interface kgen_read_var 90 procedure read_var_real_real_kind_dim4 91 end interface kgen_read_var 92 93 94 95 ! verification interface 96 interface kgen_verify_var 97 procedure verify_var_logical 98 procedure verify_var_integer 99 procedure verify_var_real 100 procedure verify_var_character 101 procedure verify_var_real_real_kind_dim4 102 end interface kgen_verify_var 103 104 INTEGER*8 :: kgen_intvar, start_clock, stop_clock, rate_clock 105 TYPE(check_t):: check_status 106 REAL(KIND=kgen_dp) :: tolerance 107 INTEGER, intent(in) :: nx 108 INTEGER, intent(in) :: qsize 109 REAL(KIND=real_kind), intent(inout) :: qdp(nx,nx,nlev,qsize) 110 REAL(KIND=real_kind), allocatable :: ref_qdp(:,:,:,:) 111 REAL(KIND=real_kind), intent(in) :: dp1(nx,nx,nlev) 112 REAL(KIND=real_kind), intent(in) :: dp2(nx,nx,nlev) 113 ! ======================== 114 ! Local Variables 115 ! ======================== 116 tolerance = 1.E-14 117 CALL kgen_init_check(check_status, tolerance) 118 ! None 119 call kgen_read_var(ref_qdp, kgen_unit) 120 ! call to kernel 121 CALL remap_q_ppm(qdp, nx, qsize, dp1, dp2) 122 ! kernel verification for output variables 123 call kgen_verify_var("qdp", check_status, qdp, ref_qdp) 124 CALL kgen_print_check("remap_q_ppm", check_status) 125 CALL system_clock(start_clock, rate_clock) 126 DO kgen_intvar=1,10 127 CALL remap_q_ppm(qdp, nx, qsize, dp1, dp2) 128 END DO 129 CALL system_clock(stop_clock, rate_clock) 130 WRITE(*,*) 131 PRINT *, "Elapsed time (sec): ", (stop_clock - start_clock)/REAL(rate_clock*10) 132 ! q loop 133 CONTAINS 134 135 ! read subroutines 136 subroutine read_var_real_real_kind_dim4(var, kgen_unit) 137 integer, intent(in) :: kgen_unit 138 real(kind=real_kind), intent(out), dimension(:,:,:,:), allocatable :: var 139 integer, dimension(2,4) :: kgen_bound 140 logical is_save 141 142 READ(UNIT = kgen_unit) is_save 143 if ( is_save ) then 144 READ(UNIT = kgen_unit) kgen_bound(1, 1) 145 READ(UNIT = kgen_unit) kgen_bound(2, 1) 146 READ(UNIT = kgen_unit) kgen_bound(1, 2) 147 READ(UNIT = kgen_unit) kgen_bound(2, 2) 148 READ(UNIT = kgen_unit) kgen_bound(1, 3) 149 READ(UNIT = kgen_unit) kgen_bound(2, 3) 150 READ(UNIT = kgen_unit) kgen_bound(1, 4) 151 READ(UNIT = kgen_unit) kgen_bound(2, 4) 152 ALLOCATE(var(kgen_bound(2, 1) - kgen_bound(1, 1) + 1, kgen_bound(2, 2) - kgen_bound(1, 2) + 1, kgen_bound(2, 3) - kgen_bound(1, 3) + 1, kgen_bound(2, 4) - kgen_bound(1, 4) + 1)) 153 READ(UNIT = kgen_unit) var 154 end if 155 end subroutine 156 157 subroutine verify_var_logical(varname, check_status, var, ref_var) 158 character(*), intent(in) :: varname 159 type(check_t), intent(inout) :: check_status 160 logical, intent(in) :: var, ref_var 161 162 check_status%numTotal = check_status%numTotal + 1 163 IF ( var .eqv. ref_var ) THEN 164 check_status%numIdentical = check_status%numIdentical + 1 165 if(check_status%verboseLevel > 1) then 166 WRITE(*,*) 167 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 168 endif 169 ELSE 170 if(check_status%verboseLevel > 1) then 171 WRITE(*,*) 172 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 173 if(check_status%verboseLevel > 2) then 174 WRITE(*,*) "KERNEL: ", var 175 WRITE(*,*) "REF. : ", ref_var 176 endif 177 endif 178 check_status%numFatal = check_status%numFatal + 1 179 END IF 180 end subroutine 181 182 subroutine verify_var_integer(varname, check_status, var, ref_var) 183 character(*), intent(in) :: varname 184 type(check_t), intent(inout) :: check_status 185 integer, intent(in) :: var, ref_var 186 187 check_status%numTotal = check_status%numTotal + 1 188 IF ( var == ref_var ) THEN 189 check_status%numIdentical = check_status%numIdentical + 1 190 if(check_status%verboseLevel > 1) then 191 WRITE(*,*) 192 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 193 endif 194 ELSE 195 if(check_status%verboseLevel > 0) then 196 WRITE(*,*) 197 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 198 if(check_status%verboseLevel > 2) then 199 WRITE(*,*) "KERNEL: ", var 200 WRITE(*,*) "REF. : ", ref_var 201 endif 202 endif 203 check_status%numFatal = check_status%numFatal + 1 204 END IF 205 end subroutine 206 207 subroutine verify_var_real(varname, check_status, var, ref_var) 208 character(*), intent(in) :: varname 209 type(check_t), intent(inout) :: check_status 210 real, intent(in) :: var, ref_var 211 212 check_status%numTotal = check_status%numTotal + 1 213 IF ( var == ref_var ) THEN 214 check_status%numIdentical = check_status%numIdentical + 1 215 if(check_status%verboseLevel > 1) then 216 WRITE(*,*) 217 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 218 endif 219 ELSE 220 if(check_status%verboseLevel > 0) then 221 WRITE(*,*) 222 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 223 if(check_status%verboseLevel > 2) then 224 WRITE(*,*) "KERNEL: ", var 225 WRITE(*,*) "REF. : ", ref_var 226 endif 227 endif 228 check_status%numFatal = check_status%numFatal + 1 229 END IF 230 end subroutine 231 232 subroutine verify_var_character(varname, check_status, var, ref_var) 233 character(*), intent(in) :: varname 234 type(check_t), intent(inout) :: check_status 235 character(*), intent(in) :: var, ref_var 236 237 check_status%numTotal = check_status%numTotal + 1 238 IF ( var == ref_var ) THEN 239 check_status%numIdentical = check_status%numIdentical + 1 240 if(check_status%verboseLevel > 1) then 241 WRITE(*,*) 242 WRITE(*,*) trim(adjustl(varname)), " is IDENTICAL( ", var, " )." 243 endif 244 ELSE 245 if(check_status%verboseLevel > 0) then 246 WRITE(*,*) 247 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 248 if(check_status%verboseLevel > 2) then 249 WRITE(*,*) "KERNEL: ", var 250 WRITE(*,*) "REF. : ", ref_var 251 end if 252 end if 253 check_status%numFatal = check_status%numFatal + 1 254 END IF 255 end subroutine 256 257 subroutine verify_var_real_real_kind_dim4(varname, check_status, var, ref_var) 258 character(*), intent(in) :: varname 259 type(check_t), intent(inout) :: check_status 260 real(kind=real_kind), intent(in), dimension(:,:,:,:) :: var 261 real(kind=real_kind), intent(in), allocatable, dimension(:,:,:,:) :: ref_var 262 real(kind=real_kind) :: nrmsdiff, rmsdiff 263 real(kind=real_kind), allocatable :: temp(:,:,:,:), temp2(:,:,:,:) 264 integer :: n 265 266 267 IF ( ALLOCATED(ref_var) ) THEN 268 check_status%numTotal = check_status%numTotal + 1 269 allocate(temp(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3),SIZE(var,dim=4))) 270 allocate(temp2(SIZE(var,dim=1),SIZE(var,dim=2),SIZE(var,dim=3),SIZE(var,dim=4))) 271 IF ( ALL( var == ref_var ) ) THEN 272 check_status%numIdentical = check_status%numIdentical + 1 273 if(check_status%verboseLevel > 1) then 274 WRITE(*,*) 275 WRITE(*,*) "All elements of ", trim(adjustl(varname)), " are IDENTICAL." 276 !WRITE(*,*) "KERNEL: ", var 277 !WRITE(*,*) "REF. : ", ref_var 278 IF ( ALL( var == 0 ) ) THEN 279 if(check_status%verboseLevel > 2) then 280 WRITE(*,*) "All values are zero." 281 end if 282 END IF 283 end if 284 ELSE 285 n = count(var/=ref_var) 286 where(ref_var .NE. 0) 287 temp = ((var-ref_var)/ref_var)**2 288 temp2 = (var-ref_var)**2 289 elsewhere 290 temp = (var-ref_var)**2 291 temp2 = temp 292 endwhere 293 nrmsdiff = sqrt(sum(temp)/real(n)) 294 rmsdiff = sqrt(sum(temp2)/real(n)) 295 296 if(check_status%verboseLevel > 0) then 297 WRITE(*,*) 298 WRITE(*,*) trim(adjustl(varname)), " is NOT IDENTICAL." 299 WRITE(*,*) count( var /= ref_var), " of ", size( var ), " elements are different." 300 if(check_status%verboseLevel > 1) then 301 WRITE(*,*) "Average - kernel ", sum(var)/real(size(var)) 302 WRITE(*,*) "Average - reference ", sum(ref_var)/real(size(ref_var)) 303 endif 304 WRITE(*,*) "RMS of difference is ",rmsdiff 305 WRITE(*,*) "Normalized RMS of difference is ",nrmsdiff 306 end if 307 308 if (nrmsdiff > check_status%tolerance) then 309 check_status%numFatal = check_status%numFatal+1 310 else 311 check_status%numWarning = check_status%numWarning+1 312 endif 313 END IF 314 deallocate(temp,temp2) 315 END IF 316 end subroutine 317 318 END SUBROUTINE remap1 319 320 !=======================================================================================================! 321 !This uses the exact same model and reference grids and data as remap_Q, but it interpolates 322 !using PPM instead of splines. 323 324 SUBROUTINE remap_q_ppm(qdp, nx, qsize, dp1, dp2) 325 ! remap 1 field 326 ! input: Qdp field to be remapped (NOTE: MASS, not MIXING RATIO) 327 ! dp1 layer thickness (source) 328 ! dp2 layer thickness (target) 329 ! 330 ! output: remaped Qdp, conserving mass 331 ! 332 USE control_mod, ONLY: vert_remap_q_alg 333 IMPLICIT NONE 334 INTEGER, intent(in) :: nx, qsize 335 REAL(KIND=real_kind), intent(inout) :: qdp(nx,nx,nlev,qsize) 336 REAL(KIND=real_kind), intent(in) :: dp1(nx,nx,nlev), dp2(nx,nx,nlev) 337 ! Local Variables 338 INTEGER, parameter :: gs = 2 !Number of cells to place in the ghost region 339 REAL(KIND=real_kind), dimension(nlev+2) :: pio !Pressure at interfaces for old grid 340 REAL(KIND=real_kind), dimension(nlev+1) :: pin !Pressure at interfaces for new grid 341 REAL(KIND=real_kind), dimension(nlev+1) :: masso !Accumulate mass up to each interface 342 REAL(KIND=real_kind), dimension(1-gs:nlev+gs) :: ao !Tracer value on old grid 343 REAL(KIND=real_kind), dimension(1-gs:nlev+gs) :: dpo !change in pressure over a cell for old grid 344 REAL(KIND=real_kind), dimension(1-gs:nlev+gs) :: dpn !change in pressure over a cell for old grid 345 REAL(KIND=real_kind), dimension(3, nlev) :: coefs !PPM coefficients within each cell 346 REAL(KIND=real_kind), dimension( nlev ) :: z1, z2 347 REAL(KIND=real_kind) :: ppmdx(10,0:nlev+1) !grid spacings 348 REAL(KIND=real_kind) :: mymass, massn1, massn2 349 INTEGER :: i, j, k, q, kk, kid(nlev) 350 CALL t_startf('remap_Q_ppm') 351 DO j = 1 , nx 352 DO i = 1 , nx 353 pin(1) = 0 354 pio(1) = 0 355 DO k=1,nlev 356 dpn(k) = dp2(i,j,k) 357 dpo(k) = dp1(i,j,k) 358 pin(k+1) = pin(k)+dpn(k) 359 pio(k+1) = pio(k)+dpo(k) 360 END DO 361 pio(nlev+2) = pio(nlev+1) + 1. !This is here to allow an entire block of k threads to run in the remapping phase. 362 !It makes sure there's an old interface value below the domain that is larger. 363 pin(nlev+1) = pio(nlev+1) !The total mass in a column does not change. 364 !Therefore, the pressure of that mass cannot either. 365 !Fill in the ghost regions with mirrored values. if vert_remap_q_alg is defined, this is of no consequence. 366 DO k = 1 , gs 367 dpo(1 -k) = dpo( k) 368 dpo(nlev+k) = dpo(nlev+1-k) 369 END DO 370 !Compute remapping intervals once for all tracers. Find the old grid cell index in which the 371 !k-th new cell interface resides. Then integrate from the bottom of that old cell to the new 372 !interface location. In practice, the grid never deforms past one cell, so the search can be 373 !simplified by this. Also, the interval of integration is usually of magnitude close to zero 374 !or close to dpo because of minimial deformation. 375 !Numerous tests confirmed that the bottom and top of the grids match to machine precision, so 376 !I set them equal to each other. 377 DO k = 1 , nlev 378 kk = k !Keep from an order n^2 search operation by assuming the old cell index is close. 379 !Find the index of the old grid cell in which this new cell's bottom interface resides. 380 DO while (pio(kk) <= pin(k+1)) 381 kk = kk + 1 382 END DO 383 kk = kk - 1 !kk is now the cell index we're integrating over. 384 IF (kk == nlev+1) kk = nlev !This is to keep the indices in bounds. 385 !Top bounds match anyway, so doesn't matter what coefficients are used 386 kid(k) = kk !Save for reuse 387 z1(k) = -0.5d0 !This remapping assumes we're starting from the left interface of an old grid cell 388 !In fact, we're usually integrating very little or almost all of the cell in question 389 z2(k) = (pin(k+1) - ( pio(kk) + pio(kk+1) ) * 0.5) / dpo(kk) !PPM interpolants are normalized to an independent 390 !coordinate domain [-0.5,0.5]. 391 END DO 392 !This turned out a big optimization, remembering that only parts of the PPM algorithm depends on the data, 393 ! namely the 394 !limiting. So anything that depends only on the grid is pre-computed outside the tracer loop. 395 ppmdx(:,:) = compute_ppm_grids( dpo ) 396 !From here, we loop over tracers for only those portions which depend on tracer data, which includes PPM 397 ! limiting and 398 !mass accumulation 399 DO q = 1 , qsize 400 !Accumulate the old mass up to old grid cell interface locations to simplify integration 401 !during remapping. Also, divide out the grid spacing so we're working with actual tracer 402 !values and can conserve mass. The option for ifndef ZEROHORZ I believe is there to ensure 403 !tracer consistency for an initially uniform field. I copied it from the old remap routine. 404 masso(1) = 0. 405 DO k = 1 , nlev 406 ao(k) = qdp(i,j,k,q) 407 masso(k+1) = masso(k) + ao(k) !Accumulate the old mass. This will simplify the remapping 408 ao(k) = ao(k) / dpo(k) !Divide out the old grid spacing because we want the tracer mixing ratio, not mass. 409 END DO 410 !Fill in ghost values. Ignored if vert_remap_q_alg == 2 411 DO k = 1 , gs 412 ao(1 -k) = ao( k) 413 ao(nlev+k) = ao(nlev+1-k) 414 END DO 415 !Compute monotonic and conservative PPM reconstruction over every cell 416 coefs(:,:) = compute_ppm(ao , ppmdx) 417 !Compute tracer values on the new grid by integrating from the old cell bottom to the new 418 !cell interface to form a new grid mass accumulation. Taking the difference between 419 !accumulation at successive interfaces gives the mass inside each cell. Since Qdp is 420 !supposed to hold the full mass this needs no normalization. 421 massn1 = 0. 422 DO k = 1 , nlev 423 kk = kid(k) 424 massn2 = masso(kk) + integrate_parabola(coefs(:,kk) , z1(k) , z2(k)) * dpo(kk) 425 qdp(i,j,k,q) = massn2 - massn1 426 massn1 = massn2 427 END DO 428 END DO 429 END DO 430 END DO 431 CALL t_stopf('remap_Q_ppm') 432 END SUBROUTINE remap_q_ppm 433 !=======================================================================================================! 434 !THis compute grid-based coefficients from Collela & Woodward 1984. 435 436 FUNCTION compute_ppm_grids(dx) RESULT ( rslt ) 437 USE control_mod, ONLY: vert_remap_q_alg 438 IMPLICIT NONE 439 REAL(KIND=real_kind), intent(in) :: dx(-1:nlev+2) !grid spacings 440 REAL(KIND=real_kind) :: rslt(10,0:nlev+1) !grid spacings 441 INTEGER :: j 442 INTEGER :: indb, inde 443 !Calculate grid-based coefficients for stage 1 of compute_ppm 444 IF (vert_remap_q_alg == 2) THEN 445 indb = 2 446 inde = nlev-1 447 ELSE 448 indb = 0 449 inde = nlev+1 450 END IF 451 DO j = indb , inde 452 rslt(1,j) = dx(j) / (dx(j-1) + dx(j) + dx(j+1)) 453 rslt(2,j) = (2.*dx(j-1) + dx(j)) / (dx(j+1) + dx(j)) 454 rslt(3,j) = (dx(j) + 2.*dx(j+1)) / (dx(j-1) + dx(j)) 455 END DO 456 !Caculate grid-based coefficients for stage 2 of compute_ppm 457 IF (vert_remap_q_alg == 2) THEN 458 indb = 2 459 inde = nlev-2 460 ELSE 461 indb = 0 462 inde = nlev 463 END IF 464 DO j = indb , inde 465 rslt(4,j) = dx(j) / (dx(j) + dx(j+1)) 466 rslt(5,j) = 1. / sum(dx(j-1:j+2)) 467 rslt(6,j) = (2. * dx(j+1) * dx(j)) / (dx(j) + dx(j+1 )) 468 rslt(7,j) = (dx(j-1) + dx(j )) / (2. * dx(j ) + dx(j+1)) 469 rslt(8,j) = (dx(j+2) + dx(j+1)) / (2. * dx(j+1) + dx(j )) 470 rslt(9,j) = dx(j ) * (dx(j-1) + dx(j )) / (2.*dx(j ) + dx(j+1)) 471 rslt(10,j) = dx(j+1) * (dx(j+1) + dx(j+2)) / (dx(j ) + 2.*dx(j+1)) 472 END DO 473 END FUNCTION compute_ppm_grids 474 !=======================================================================================================! 475 !This computes a limited parabolic interpolant using a net 5-cell stencil, but the stages of computation are broken up 476 ! into 3 stages 477 478 FUNCTION compute_ppm(a, dx) RESULT ( coefs ) 479 USE control_mod, ONLY: vert_remap_q_alg 480 IMPLICIT NONE 481 REAL(KIND=real_kind), intent(in) :: a (-1:nlev+2) !Cell-mean values 482 REAL(KIND=real_kind), intent(in) :: dx (10, 0:nlev+1) !grid spacings 483 REAL(KIND=real_kind) :: coefs(0:2, nlev) !PPM coefficients (for parabola) 484 REAL(KIND=real_kind) :: ai (0:nlev) !fourth-order accurate, then limited interface values 485 REAL(KIND=real_kind) :: dma(0:nlev+1) !An expression from Collela's '84 publication 486 REAL(KIND=real_kind) :: da !Ditto 487 ! Hold expressions based on the grid (which are cumbersome). 488 REAL(KIND=real_kind) :: dx1, dx2, dx3, dx4, dx5, dx6, dx7, dx8, dx9, dx10 489 REAL(KIND=real_kind) :: al, ar !Left and right interface values for cell-local limiting 490 INTEGER :: j 491 INTEGER :: indb, inde 492 ! Stage 1: Compute dma for each cell, allowing a 1-cell ghost stencil below and above the domain 493 IF (vert_remap_q_alg == 2) THEN 494 indb = 2 495 inde = nlev-1 496 ELSE 497 indb = 0 498 inde = nlev+1 499 END IF 500 DO j = indb , inde 501 da = dx(1,j) * (dx(2,j) * ( a(j+1) - a(j) ) + dx(3,j) * ( a(j) - a(j-1) )) 502 dma(j) = minval((/ abs(da) , 2. * abs( a(j) - a(j-1) ) , 2. * abs( a(j+1) - a(j) ) /)) * sign(1.d0,da) 503 IF (( a(j+1) - a(j) ) * ( a(j) - a(j-1) ) <= 0.) dma(j) = 0. 504 END DO 505 ! Stage 2: Compute ai for each cell interface in the physical domain (dimension nlev+1) 506 IF (vert_remap_q_alg == 2) THEN 507 indb = 2 508 inde = nlev-2 509 ELSE 510 indb = 0 511 inde = nlev 512 END IF 513 DO j = indb , inde 514 ai(j) = a(j) + dx(4,j) * (a(j+1) - a(j)) + dx(5,j) * (dx(6,j) * ( dx(7,j) - dx(8,j) ) * ( a(j+1) - a(j) )& 515 - dx(9,j) * dma(j+1) + dx(10,j) * dma(j)) 516 END DO 517 ! Stage 3: Compute limited PPM interpolant over each cell in the physical domain 518 ! (dimension nlev) using ai on either side and ao within the cell. 519 IF (vert_remap_q_alg == 2) THEN 520 indb = 3 521 inde = nlev-2 522 ELSE 523 indb = 1 524 inde = nlev 525 END IF 526 DO j = indb , inde 527 al = ai(j-1) 528 ar = ai(j ) 529 IF ((ar - a(j)) * (a(j) - al) <= 0.) THEN 530 al = a(j) 531 ar = a(j) 532 END IF 533 IF ((ar - al) * (a(j) - (al + ar)/2.) > (ar - al)**2/6.) al = 3.*a(j) - 2. * ar 534 IF ((ar - al) * (a(j) - (al + ar)/2.) < -(ar - al)**2/6.) ar = 3.*a(j) - 2. * al 535 !Computed these coefficients from the edge values and cell mean in Maple. Assumes normalized coordinates: xi=( 536 ! x-x0)/dx 537 coefs(0,j) = 1.5 * a(j) - (al + ar) / 4. 538 coefs(1,j) = ar - al 539 coefs(2,j) = -6. * a(j) + 3. * (al + ar) 540 END DO 541 !If we're not using a mirrored boundary condition, then make the two cells bordering the top and bottom 542 !material boundaries piecewise constant. Zeroing out the first and second moments, and setting the zeroth 543 !moment to the cell mean is sufficient to maintain conservation. 544 IF (vert_remap_q_alg == 2) THEN 545 coefs(0,1:2) = a(1:2) 546 coefs(1:2,1:2) = 0. 547 coefs(0,nlev-1:nlev) = a(nlev-1:nlev) 548 coefs(1:2,nlev-1:nlev) = 0.d0 549 END IF 550 END FUNCTION compute_ppm 551 !=======================================================================================================! 552 !Simple function computes the definite integral of a parabola in normalized coordinates, xi=(x-x0)/dx, 553 !given two bounds. Make sure this gets inlined during compilation. 554 555 FUNCTION integrate_parabola(a, x1, x2) RESULT ( mass ) 556 IMPLICIT NONE 557 REAL(KIND=real_kind), intent(in) :: a(0:2) !Coefficients of the parabola 558 REAL(KIND=real_kind), intent(in) :: x1 !lower domain bound for integration 559 REAL(KIND=real_kind), intent(in) :: x2 !upper domain bound for integration 560 REAL(KIND=real_kind) :: mass 561 mass = a(0) * (x2 - x1) + a(1) * (x2 ** 2 - x1 ** 2) / 0.2d1 + a(2) * (x2 ** 3 - x1 ** 3) / 0.3d1 562 END FUNCTION integrate_parabola 563 !=============================================================================================! 564 END MODULE vertremap_mod 565 566 567 568 569 570 571 572 573 574 575 576 577 578 579 580 581 582 583 584 585 586 587 588 589 590 591 592 593 594