| 1 | |
|---|
| 2 | ! |
|---|
| 3 | ! This program may be freely redistributed under the |
|---|
| 4 | ! condition that the copyright notices (including this |
|---|
| 5 | ! entire header) are not removed, and no compensation |
|---|
| 6 | ! is received through use of the software. Private, |
|---|
| 7 | ! research, and institutional use is free. You may |
|---|
| 8 | ! distribute modified versions of this code UNDER THE |
|---|
| 9 | ! CONDITION THAT THIS CODE AND ANY MODIFICATIONS MADE |
|---|
| 10 | ! TO IT IN THE SAME FILE REMAIN UNDER COPYRIGHT OF THE |
|---|
| 11 | ! ORIGINAL AUTHOR, BOTH SOURCE AND OBJECT CODE ARE |
|---|
| 12 | ! MADE FREELY AVAILABLE WITHOUT CHARGE, AND CLEAR |
|---|
| 13 | ! NOTICE IS GIVEN OF THE MODIFICATIONS. Distribution |
|---|
| 14 | ! of this code as part of a commercial system is |
|---|
| 15 | ! permissible ONLY BY DIRECT ARRANGEMENT WITH THE |
|---|
| 16 | ! AUTHOR. (If you are not directly supplying this |
|---|
| 17 | ! code to a customer, and you are instead telling them |
|---|
| 18 | ! how they can obtain it for free, then you are not |
|---|
| 19 | ! required to make any arrangement with me.) |
|---|
| 20 | ! |
|---|
| 21 | ! Disclaimer: Neither I nor: Columbia University, the |
|---|
| 22 | ! National Aeronautics and Space Administration, nor |
|---|
| 23 | ! the Massachusetts Institute of Technology warrant |
|---|
| 24 | ! or certify this code in any way whatsoever. This |
|---|
| 25 | ! code is provided "as-is" to be used at your own risk. |
|---|
| 26 | ! |
|---|
| 27 | ! |
|---|
| 28 | |
|---|
| 29 | ! |
|---|
| 30 | ! PQM.f90: a 1d slope-limited, piecewise quartic recon. |
|---|
| 31 | ! |
|---|
| 32 | ! Darren Engwirda |
|---|
| 33 | ! 08-Sep-2016 |
|---|
| 34 | ! de2363 [at] columbia [dot] edu |
|---|
| 35 | ! |
|---|
| 36 | ! |
|---|
| 37 | |
|---|
| 38 | ! White, L. and Adcroft, A., A high-order finite volume |
|---|
| 39 | ! remapping scheme for nonuniform grids: The piecewise |
|---|
| 40 | ! quartic method (PQM), J. Comp. Phys., 227 (15), 2008, |
|---|
| 41 | ! 7394-7422, https://doi.org/10.1016/j.jcp.2008.04.026. |
|---|
| 42 | ! |
|---|
| 43 | |
|---|
| 44 | pure subroutine pqm(npos,nvar,ndof,delx, & |
|---|
| 45 | & fdat,fhat,edge,dfdx, & |
|---|
| 46 | & oscl,dmin,ilim,wlim, & |
|---|
| 47 | & halo) |
|---|
| 48 | |
|---|
| 49 | ! |
|---|
| 50 | ! NPOS no. edges over grid. |
|---|
| 51 | ! NVAR no. state variables. |
|---|
| 52 | ! NDOF no. degrees-of-freedom per grid-cell. |
|---|
| 53 | ! DELX grid-cell spacing array. LENGTH(DELX) == +1 if |
|---|
| 54 | ! spacing is uniform . |
|---|
| 55 | ! FDAT grid-cell moments array. FDAT is an array with |
|---|
| 56 | ! SIZE = NDOF-by-NVAR-by-NPOS-1 . |
|---|
| 57 | ! FHAT grid-cell re-con. array. FHAT is an array with |
|---|
| 58 | ! SIZE = MDOF-by-NVAR-by-NPOS-1 . |
|---|
| 59 | ! EDGE edge-centred interp. for function-value. EDGE |
|---|
| 60 | ! is an array with SIZE = NVAR-by-NPOS . |
|---|
| 61 | ! DFDX edge-centred interp. for 1st-derivative. DFDX |
|---|
| 62 | ! is an array with SIZE = NVAR-by-NPOS . |
|---|
| 63 | ! OSCL grid-cell oscil. dof.'s. OSCL is an array with |
|---|
| 64 | ! SIZE = +2 -by-NVAR-by-NPOS-1 . |
|---|
| 65 | ! DMIN min. grid-cell spacing thresh . |
|---|
| 66 | ! ILIM cell slope-limiting selection . |
|---|
| 67 | ! WLIM wall slope-limiting selection . |
|---|
| 68 | ! HALO width of re-con. stencil, symmetric about mid. . |
|---|
| 69 | ! |
|---|
| 70 | |
|---|
| 71 | implicit none |
|---|
| 72 | |
|---|
| 73 | !------------------------------------------- arguments ! |
|---|
| 74 | integer, intent(in) :: npos,nvar,ndof |
|---|
| 75 | integer, intent(in) :: ilim,wlim,halo |
|---|
| 76 | real*8 , intent(in) :: dmin |
|---|
| 77 | real*8 , intent(out) :: fhat(:,:,:) |
|---|
| 78 | real*8 , intent(in) :: oscl(:,:,:) |
|---|
| 79 | real*8 , intent(in) :: delx(:) |
|---|
| 80 | real*8 , intent(in) :: fdat(:,:,:) |
|---|
| 81 | real*8 , intent(in) :: edge(:,:) |
|---|
| 82 | real*8 , intent(in) :: dfdx(:,:) |
|---|
| 83 | |
|---|
| 84 | !------------------------------------------- variables ! |
|---|
| 85 | integer :: ipos,ivar,iill,iirr,head,tail |
|---|
| 86 | real*8 :: ff00,ffll,ffrr,hh00,hhll,hhrr |
|---|
| 87 | real*8 :: xhat |
|---|
| 88 | integer :: mono |
|---|
| 89 | real*8 :: fell,ferr |
|---|
| 90 | real*8 :: dell,derr |
|---|
| 91 | real*8 :: dfds(-1:+1) |
|---|
| 92 | real*8 :: uhat(+1:+5) |
|---|
| 93 | real*8 :: lhat(+1:+5) |
|---|
| 94 | real*8 :: wval(+1:+2) |
|---|
| 95 | |
|---|
| 96 | head = +1; tail = npos - 1 |
|---|
| 97 | |
|---|
| 98 | if (npos.le.2) then |
|---|
| 99 | !----- default to reduced order if insufficient points ! |
|---|
| 100 | do ivar = +1, nvar |
|---|
| 101 | fhat(1,ivar,+1) = fdat(1,ivar,+1) |
|---|
| 102 | fhat(2,ivar,+1) = 0.d0 |
|---|
| 103 | fhat(3,ivar,+1) = 0.d0 |
|---|
| 104 | fhat(4,ivar,+1) = 0.d0 |
|---|
| 105 | fhat(5,ivar,+1) = 0.d0 |
|---|
| 106 | end do |
|---|
| 107 | end if |
|---|
| 108 | |
|---|
| 109 | if (npos.le.2) return |
|---|
| 110 | |
|---|
| 111 | !------------------- reconstruct function on each cell ! |
|---|
| 112 | |
|---|
| 113 | do ipos = +1 , npos-1 |
|---|
| 114 | |
|---|
| 115 | iill = max(head,ipos-1) |
|---|
| 116 | iirr = min(tail,ipos+1) |
|---|
| 117 | |
|---|
| 118 | do ivar = +1 , nvar-0 |
|---|
| 119 | |
|---|
| 120 | !----------------------------- cell mean + edge values ! |
|---|
| 121 | |
|---|
| 122 | ff00 = fdat(1,ivar,ipos) |
|---|
| 123 | ffll = fdat(1,ivar,iill) |
|---|
| 124 | ffrr = fdat(1,ivar,iirr) |
|---|
| 125 | |
|---|
| 126 | fell = edge(ivar,ipos+0) |
|---|
| 127 | ferr = edge(ivar,ipos+1) |
|---|
| 128 | |
|---|
| 129 | !----------------------------- calc. LL/00/RR gradient ! |
|---|
| 130 | |
|---|
| 131 | if (size(delx).gt.+1) then |
|---|
| 132 | |
|---|
| 133 | hh00 = delx(ipos) |
|---|
| 134 | hhll = delx(iill) |
|---|
| 135 | hhrr = delx(iirr) |
|---|
| 136 | |
|---|
| 137 | xhat = delx(ipos+0)*.5d+0 |
|---|
| 138 | |
|---|
| 139 | call plsv (ffll,hhll,ff00, & |
|---|
| 140 | & hh00,ffrr,hhrr, & |
|---|
| 141 | & dfds) |
|---|
| 142 | else |
|---|
| 143 | |
|---|
| 144 | xhat = delx( +1)*.5d+0 |
|---|
| 145 | |
|---|
| 146 | call plsc (ffll,ff00,ffrr, & |
|---|
| 147 | & dfds) |
|---|
| 148 | |
|---|
| 149 | end if |
|---|
| 150 | |
|---|
| 151 | dell = dfdx (ivar,ipos+0) |
|---|
| 152 | dell = dell * xhat |
|---|
| 153 | |
|---|
| 154 | derr = dfdx (ivar,ipos+1) |
|---|
| 155 | derr = derr * xhat |
|---|
| 156 | |
|---|
| 157 | !----------------------------- calc. cell-wise profile ! |
|---|
| 158 | |
|---|
| 159 | select case(ilim) |
|---|
| 160 | case (null_limit) |
|---|
| 161 | |
|---|
| 162 | !----------------------------- calc. unlimited profile ! |
|---|
| 163 | |
|---|
| 164 | call pqmfn(ff00,ffll,ffrr, & |
|---|
| 165 | & fell,ferr,dell, & |
|---|
| 166 | & derr,dfds,uhat, & |
|---|
| 167 | & lhat,mono) |
|---|
| 168 | |
|---|
| 169 | !----------------------------- pref. unlimited profile ! |
|---|
| 170 | |
|---|
| 171 | wval(1) = +1.d+0 |
|---|
| 172 | wval(2) = +0.d+0 |
|---|
| 173 | |
|---|
| 174 | case (mono_limit) |
|---|
| 175 | |
|---|
| 176 | !----------------------------- calc. monotonic profile ! |
|---|
| 177 | |
|---|
| 178 | call pqmfn(ff00,ffll,ffrr, & |
|---|
| 179 | & fell,ferr,dell, & |
|---|
| 180 | & derr,dfds,uhat, & |
|---|
| 181 | & lhat,mono) |
|---|
| 182 | |
|---|
| 183 | !----------------------------- pref. monotonic profile ! |
|---|
| 184 | |
|---|
| 185 | wval(1) = +0.d+0 |
|---|
| 186 | wval(2) = +1.d+0 |
|---|
| 187 | |
|---|
| 188 | case (weno_limit) |
|---|
| 189 | |
|---|
| 190 | !----------------------------- calc. monotonic profile ! |
|---|
| 191 | |
|---|
| 192 | call pqmfn(ff00,ffll,ffrr, & |
|---|
| 193 | & fell,ferr,dell, & |
|---|
| 194 | & derr,dfds,uhat, & |
|---|
| 195 | & lhat,mono) |
|---|
| 196 | |
|---|
| 197 | if (mono.gt.+0) then |
|---|
| 198 | |
|---|
| 199 | !----------------------------- calc. WENO-type weights ! |
|---|
| 200 | |
|---|
| 201 | call wenoi(npos,delx,oscl, & |
|---|
| 202 | & ipos,ivar,halo, & |
|---|
| 203 | & wlim,wval) |
|---|
| 204 | |
|---|
| 205 | else |
|---|
| 206 | |
|---|
| 207 | !----------------------------- pref. unlimited profile ! |
|---|
| 208 | |
|---|
| 209 | wval(1) = +1.d+0 |
|---|
| 210 | wval(2) = +0.d+0 |
|---|
| 211 | |
|---|
| 212 | end if |
|---|
| 213 | |
|---|
| 214 | end select |
|---|
| 215 | |
|---|
| 216 | !----------------------------- blend "null" and "mono" ! |
|---|
| 217 | |
|---|
| 218 | fhat(1,ivar,ipos) = & |
|---|
| 219 | & wval(1) * uhat(1) + & |
|---|
| 220 | & wval(2) * lhat(1) |
|---|
| 221 | fhat(2,ivar,ipos) = & |
|---|
| 222 | & wval(1) * uhat(2) + & |
|---|
| 223 | & wval(2) * lhat(2) |
|---|
| 224 | fhat(3,ivar,ipos) = & |
|---|
| 225 | & wval(1) * uhat(3) + & |
|---|
| 226 | & wval(2) * lhat(3) |
|---|
| 227 | fhat(4,ivar,ipos) = & |
|---|
| 228 | & wval(1) * uhat(4) + & |
|---|
| 229 | & wval(2) * lhat(4) |
|---|
| 230 | fhat(5,ivar,ipos) = & |
|---|
| 231 | & wval(1) * uhat(5) + & |
|---|
| 232 | & wval(2) * lhat(5) |
|---|
| 233 | |
|---|
| 234 | end do |
|---|
| 235 | |
|---|
| 236 | end do |
|---|
| 237 | |
|---|
| 238 | return |
|---|
| 239 | |
|---|
| 240 | end subroutine |
|---|
| 241 | |
|---|
| 242 | !----------- assemble piecewise quartic reconstruction ! |
|---|
| 243 | |
|---|
| 244 | pure subroutine pqmfn(ff00,ffll,ffrr,fell, & |
|---|
| 245 | & ferr,dell,derr,dfds, & |
|---|
| 246 | & uhat,lhat,mono) |
|---|
| 247 | |
|---|
| 248 | ! |
|---|
| 249 | ! FF00 centred grid-cell mean. |
|---|
| 250 | ! FFLL left -biased grid-cell mean. |
|---|
| 251 | ! FFRR right-biased grid-cell mean. |
|---|
| 252 | ! FELL left -biased edge interp. |
|---|
| 253 | ! FERR right-biased edge interp. |
|---|
| 254 | ! DELL left -biased edge df//dx. |
|---|
| 255 | ! DERR right-biased edge df//dx. |
|---|
| 256 | ! DFDS piecewise linear gradients in local co-ord.'s. |
|---|
| 257 | ! DFDS(+0) is a centred, slope-limited estimate, |
|---|
| 258 | ! DFDS(-1), DFDS(+1) are left- and right-biased |
|---|
| 259 | ! estimates (unlimited). |
|---|
| 260 | ! UHAT unlimited PPM reconstruction coefficients . |
|---|
| 261 | ! LHAT monotonic PPM reconstruction coefficients . |
|---|
| 262 | ! MONO slope-limiting indicator, MONO > +0 if some |
|---|
| 263 | ! limiting has occured . |
|---|
| 264 | ! |
|---|
| 265 | |
|---|
| 266 | implicit none |
|---|
| 267 | |
|---|
| 268 | !------------------------------------------- arguments ! |
|---|
| 269 | real*8 , intent(in) :: ff00 |
|---|
| 270 | real*8 , intent(in) :: ffll,ffrr |
|---|
| 271 | real*8 , intent(inout) :: fell,ferr |
|---|
| 272 | real*8 , intent(inout) :: dell,derr |
|---|
| 273 | real*8 , intent(in) :: dfds(-1:+1) |
|---|
| 274 | real*8 , intent(out) :: uhat(+1:+5) |
|---|
| 275 | real*8 , intent(out) :: lhat(+1:+5) |
|---|
| 276 | integer, intent(out) :: mono |
|---|
| 277 | |
|---|
| 278 | !------------------------------------------- variables ! |
|---|
| 279 | integer :: turn |
|---|
| 280 | real*8 :: grad, iflx(+1:+2) |
|---|
| 281 | logical :: haveroot |
|---|
| 282 | |
|---|
| 283 | !-------------------------------- "null" slope-limiter ! |
|---|
| 284 | |
|---|
| 285 | mono = 0 |
|---|
| 286 | |
|---|
| 287 | uhat(1) = & |
|---|
| 288 | & + (30.d+0 / 16.d+0) * ff00 & |
|---|
| 289 | & - ( 7.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 290 | & + ( 1.d+0 / 16.d+0) *(derr-dell) |
|---|
| 291 | uhat(2) = & |
|---|
| 292 | & + ( 3.d+0 / 4.d+0) *(ferr-fell) & |
|---|
| 293 | & - ( 1.d+0 / 4.d+0) *(derr+dell) |
|---|
| 294 | uhat(3) = & |
|---|
| 295 | & - (30.d+0 / 8.d+0) * ff00 & |
|---|
| 296 | & + (15.d+0 / 8.d+0) *(ferr+fell) & |
|---|
| 297 | & - ( 3.d+0 / 8.d+0) *(derr-dell) |
|---|
| 298 | uhat(4) = & |
|---|
| 299 | & - ( 1.d+0 / 4.d+0) *(ferr-fell & |
|---|
| 300 | & -derr-dell) |
|---|
| 301 | uhat(5) = & |
|---|
| 302 | & + (30.d+0 / 16.d+0) * ff00 & |
|---|
| 303 | & - (15.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 304 | & + ( 5.d+0 / 16.d+0) *(derr-dell) |
|---|
| 305 | |
|---|
| 306 | !-------------------------------- "mono" slope-limiter ! |
|---|
| 307 | |
|---|
| 308 | if((ffrr - ff00) * & |
|---|
| 309 | & (ff00 - ffll) .le. 0.d+0) then |
|---|
| 310 | |
|---|
| 311 | !----------------------------------- "flatten" extrema ! |
|---|
| 312 | |
|---|
| 313 | mono = +1 |
|---|
| 314 | |
|---|
| 315 | lhat(1) = ff00 |
|---|
| 316 | lhat(2) = 0.d0 |
|---|
| 317 | lhat(3) = 0.d0 |
|---|
| 318 | lhat(4) = 0.d0 |
|---|
| 319 | lhat(5) = 0.d0 |
|---|
| 320 | |
|---|
| 321 | return |
|---|
| 322 | |
|---|
| 323 | end if |
|---|
| 324 | |
|---|
| 325 | !----------------------------------- limit edge values ! |
|---|
| 326 | |
|---|
| 327 | if((ffll - fell) * & |
|---|
| 328 | & (fell - ff00) .le. 0.d+0) then |
|---|
| 329 | |
|---|
| 330 | mono = +1 |
|---|
| 331 | |
|---|
| 332 | fell = ff00 - dfds(0) |
|---|
| 333 | |
|---|
| 334 | end if |
|---|
| 335 | |
|---|
| 336 | if (dell * dfds(0) .lt. 0.d+0) then |
|---|
| 337 | |
|---|
| 338 | mono = +1 |
|---|
| 339 | |
|---|
| 340 | dell = dfds(0) |
|---|
| 341 | |
|---|
| 342 | end if |
|---|
| 343 | |
|---|
| 344 | if((ffrr - ferr) * & |
|---|
| 345 | & (ferr - ff00) .le. 0.d+0) then |
|---|
| 346 | |
|---|
| 347 | mono = +1 |
|---|
| 348 | |
|---|
| 349 | ferr = ff00 + dfds(0) |
|---|
| 350 | |
|---|
| 351 | end if |
|---|
| 352 | |
|---|
| 353 | if (derr * dfds(0) .lt. 0.d+0) then |
|---|
| 354 | |
|---|
| 355 | mono = +1 |
|---|
| 356 | |
|---|
| 357 | derr = dfds(0) |
|---|
| 358 | |
|---|
| 359 | end if |
|---|
| 360 | |
|---|
| 361 | !----------------------------------- limit cell values ! |
|---|
| 362 | |
|---|
| 363 | lhat(1) = & |
|---|
| 364 | & + (30.d+0 / 16.d+0) * ff00 & |
|---|
| 365 | & - ( 7.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 366 | & + ( 1.d+0 / 16.d+0) *(derr-dell) |
|---|
| 367 | lhat(2) = & |
|---|
| 368 | & + ( 3.d+0 / 4.d+0) *(ferr-fell) & |
|---|
| 369 | & - ( 1.d+0 / 4.d+0) *(derr+dell) |
|---|
| 370 | lhat(3) = & |
|---|
| 371 | & - (30.d+0 / 8.d+0) * ff00 & |
|---|
| 372 | & + (15.d+0 / 8.d+0) *(ferr+fell) & |
|---|
| 373 | & - ( 3.d+0 / 8.d+0) *(derr-dell) |
|---|
| 374 | lhat(4) = & |
|---|
| 375 | & - ( 1.d+0 / 4.d+0) *(ferr-fell & |
|---|
| 376 | & -derr-dell) |
|---|
| 377 | lhat(5) = & |
|---|
| 378 | & + (30.d+0 / 16.d+0) * ff00 & |
|---|
| 379 | & - (15.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 380 | & + ( 5.d+0 / 16.d+0) *(derr-dell) |
|---|
| 381 | |
|---|
| 382 | !------------------ calc. inflexion via 2nd-derivative ! |
|---|
| 383 | |
|---|
| 384 | call roots_2(12.d+0 * lhat(5), & |
|---|
| 385 | & 6.d+0 * lhat(4), & |
|---|
| 386 | & 2.d+0 * lhat(3), & |
|---|
| 387 | & iflx , haveroot ) |
|---|
| 388 | |
|---|
| 389 | if (haveroot) then |
|---|
| 390 | |
|---|
| 391 | turn = +0 |
|---|
| 392 | |
|---|
| 393 | if ( ( iflx(1) .gt. -1.d+0 ) & |
|---|
| 394 | & .and. ( iflx(1) .lt. +1.d+0 ) ) then |
|---|
| 395 | |
|---|
| 396 | !------------------ check for non-monotonic inflection ! |
|---|
| 397 | |
|---|
| 398 | grad = lhat(2) & |
|---|
| 399 | &+ (iflx(1)**1) * 2.d+0* lhat(3) & |
|---|
| 400 | &+ (iflx(1)**2) * 3.d+0* lhat(4) & |
|---|
| 401 | &+ (iflx(1)**3) * 4.d+0* lhat(5) |
|---|
| 402 | |
|---|
| 403 | if (grad * dfds(0) .lt. 0.d+0) then |
|---|
| 404 | |
|---|
| 405 | if (abs(dfds(-1)) & |
|---|
| 406 | & .lt. abs(dfds(+1)) ) then |
|---|
| 407 | |
|---|
| 408 | turn = -1 ! modify L |
|---|
| 409 | |
|---|
| 410 | else |
|---|
| 411 | |
|---|
| 412 | turn = +1 ! modify R |
|---|
| 413 | |
|---|
| 414 | end if |
|---|
| 415 | |
|---|
| 416 | end if |
|---|
| 417 | |
|---|
| 418 | end if |
|---|
| 419 | |
|---|
| 420 | if ( ( iflx(2) .gt. -1.d+0 ) & |
|---|
| 421 | & .and. ( iflx(2) .lt. +1.d+0 ) ) then |
|---|
| 422 | |
|---|
| 423 | !------------------ check for non-monotonic inflection ! |
|---|
| 424 | |
|---|
| 425 | grad = lhat(2) & |
|---|
| 426 | &+ (iflx(2)**1) * 2.d+0* lhat(3) & |
|---|
| 427 | &+ (iflx(2)**2) * 3.d+0* lhat(4) & |
|---|
| 428 | &+ (iflx(2)**3) * 4.d+0* lhat(5) |
|---|
| 429 | |
|---|
| 430 | if (grad * dfds(0) .lt. 0.d+0) then |
|---|
| 431 | |
|---|
| 432 | if (abs(dfds(-1)) & |
|---|
| 433 | & .lt. abs(dfds(+1)) ) then |
|---|
| 434 | |
|---|
| 435 | turn = -1 ! modify L |
|---|
| 436 | |
|---|
| 437 | else |
|---|
| 438 | |
|---|
| 439 | turn = +1 ! modify R |
|---|
| 440 | |
|---|
| 441 | end if |
|---|
| 442 | |
|---|
| 443 | end if |
|---|
| 444 | |
|---|
| 445 | end if |
|---|
| 446 | |
|---|
| 447 | !------------------ pop non-monotone inflexion to edge ! |
|---|
| 448 | |
|---|
| 449 | if (turn .eq. -1) then |
|---|
| 450 | |
|---|
| 451 | !------------------ pop inflection points onto -1 edge ! |
|---|
| 452 | |
|---|
| 453 | mono = +2 |
|---|
| 454 | |
|---|
| 455 | derr = & |
|---|
| 456 | &- ( 5.d+0 / 1.d+0) * ff00 & |
|---|
| 457 | &+ ( 3.d+0 / 1.d+0) * ferr & |
|---|
| 458 | &+ ( 2.d+0 / 1.d+0) * fell |
|---|
| 459 | dell = & |
|---|
| 460 | &+ ( 5.d+0 / 3.d+0) * ff00 & |
|---|
| 461 | &- ( 1.d+0 / 3.d+0) * ferr & |
|---|
| 462 | &- ( 4.d+0 / 3.d+0) * fell |
|---|
| 463 | |
|---|
| 464 | if (dell*dfds(+0) .lt. 0.d+0) then |
|---|
| 465 | |
|---|
| 466 | dell = 0.d+0 |
|---|
| 467 | |
|---|
| 468 | ferr = & |
|---|
| 469 | &+ ( 5.d+0 / 1.d+0) * ff00 & |
|---|
| 470 | &- ( 4.d+0 / 1.d+0) * fell |
|---|
| 471 | derr = & |
|---|
| 472 | &+ (10.d+0 / 1.d+0) * ff00 & |
|---|
| 473 | &- (10.d+0 / 1.d+0) * fell |
|---|
| 474 | |
|---|
| 475 | else & |
|---|
| 476 | & if (derr*dfds(+0) .lt. 0.d+0) then |
|---|
| 477 | |
|---|
| 478 | derr = 0.d+0 |
|---|
| 479 | |
|---|
| 480 | fell = & |
|---|
| 481 | &+ ( 5.d+0 / 2.d+0) * ff00 & |
|---|
| 482 | &- ( 3.d+0 / 2.d+0) * ferr |
|---|
| 483 | dell = & |
|---|
| 484 | &- ( 5.d+0 / 3.d+0) * ff00 & |
|---|
| 485 | &+ ( 5.d+0 / 3.d+0) * ferr |
|---|
| 486 | |
|---|
| 487 | end if |
|---|
| 488 | |
|---|
| 489 | lhat(1) = & |
|---|
| 490 | &+ (30.d+0 / 16.d+0) * ff00 & |
|---|
| 491 | &- ( 7.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 492 | &+ ( 1.d+0 / 16.d+0) *(derr-dell) |
|---|
| 493 | lhat(2) = & |
|---|
| 494 | &+ ( 3.d+0 / 4.d+0) *(ferr-fell) & |
|---|
| 495 | &- ( 1.d+0 / 4.d+0) *(derr+dell) |
|---|
| 496 | lhat(3) = & |
|---|
| 497 | &- (30.d+0 / 8.d+0) * ff00 & |
|---|
| 498 | &+ (15.d+0 / 8.d+0) *(ferr+fell) & |
|---|
| 499 | &- ( 3.d+0 / 8.d+0) *(derr-dell) |
|---|
| 500 | lhat(4) = & |
|---|
| 501 | &- ( 1.d+0 / 4.d+0) *(ferr-fell & |
|---|
| 502 | & -derr-dell) |
|---|
| 503 | lhat(5) = & |
|---|
| 504 | &+ (30.d+0 / 16.d+0) * ff00 & |
|---|
| 505 | &- (15.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 506 | &+ ( 5.d+0 / 16.d+0) *(derr-dell) |
|---|
| 507 | |
|---|
| 508 | end if |
|---|
| 509 | |
|---|
| 510 | if (turn .eq. +1) then |
|---|
| 511 | |
|---|
| 512 | !------------------ pop inflection points onto -1 edge ! |
|---|
| 513 | |
|---|
| 514 | mono = +2 |
|---|
| 515 | |
|---|
| 516 | derr = & |
|---|
| 517 | &- ( 5.d+0 / 3.d+0) * ff00 & |
|---|
| 518 | &+ ( 4.d+0 / 3.d+0) * ferr & |
|---|
| 519 | &+ ( 1.d+0 / 3.d+0) * fell |
|---|
| 520 | dell = & |
|---|
| 521 | &+ ( 5.d+0 / 1.d+0) * ff00 & |
|---|
| 522 | &- ( 2.d+0 / 1.d+0) * ferr & |
|---|
| 523 | &- ( 3.d+0 / 1.d+0) * fell |
|---|
| 524 | |
|---|
| 525 | if (dell*dfds(+0) .lt. 0.d+0) then |
|---|
| 526 | |
|---|
| 527 | dell = 0.d+0 |
|---|
| 528 | |
|---|
| 529 | ferr = & |
|---|
| 530 | &+ ( 5.d+0 / 2.d+0) * ff00 & |
|---|
| 531 | &- ( 3.d+0 / 2.d+0) * fell |
|---|
| 532 | derr = & |
|---|
| 533 | &+ ( 5.d+0 / 3.d+0) * ff00 & |
|---|
| 534 | &- ( 5.d+0 / 3.d+0) * fell |
|---|
| 535 | |
|---|
| 536 | else & |
|---|
| 537 | & if (derr*dfds(+0) .lt. 0.d+0) then |
|---|
| 538 | |
|---|
| 539 | derr = 0.d+0 |
|---|
| 540 | |
|---|
| 541 | fell = & |
|---|
| 542 | &+ ( 5.d+0 / 1.d+0) * ff00 & |
|---|
| 543 | &- ( 4.d+0 / 1.d+0) * ferr |
|---|
| 544 | dell = & |
|---|
| 545 | &- (10.d+0 / 1.d+0) * ff00 & |
|---|
| 546 | &+ (10.d+0 / 1.d+0) * ferr |
|---|
| 547 | |
|---|
| 548 | end if |
|---|
| 549 | |
|---|
| 550 | lhat(1) = & |
|---|
| 551 | &+ (30.d+0 / 16.d+0) * ff00 & |
|---|
| 552 | &- ( 7.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 553 | &+ ( 1.d+0 / 16.d+0) *(derr-dell) |
|---|
| 554 | lhat(2) = & |
|---|
| 555 | &+ ( 3.d+0 / 4.d+0) *(ferr-fell) & |
|---|
| 556 | &- ( 1.d+0 / 4.d+0) *(derr+dell) |
|---|
| 557 | lhat(3) = & |
|---|
| 558 | &- (30.d+0 / 8.d+0) * ff00 & |
|---|
| 559 | &+ (15.d+0 / 8.d+0) *(ferr+fell) & |
|---|
| 560 | &- ( 3.d+0 / 8.d+0) *(derr-dell) |
|---|
| 561 | lhat(4) = & |
|---|
| 562 | &- ( 1.d+0 / 4.d+0) *(ferr-fell & |
|---|
| 563 | & -derr-dell) |
|---|
| 564 | lhat(5) = & |
|---|
| 565 | &+ (30.d+0 / 16.d+0) * ff00 & |
|---|
| 566 | &- (15.d+0 / 16.d+0) *(ferr+fell) & |
|---|
| 567 | &+ ( 5.d+0 / 16.d+0) *(derr-dell) |
|---|
| 568 | |
|---|
| 569 | end if |
|---|
| 570 | |
|---|
| 571 | end if ! haveroot |
|---|
| 572 | |
|---|
| 573 | return |
|---|
| 574 | |
|---|
| 575 | end subroutine |
|---|
| 576 | |
|---|
| 577 | |
|---|
| 578 | |
|---|
