1!/*****************************************************************************/ 2! * 3! * Elmer, A Finite Element Software for Multiphysical Problems 4! * 5! * Copyright 1st April 1995 - , CSC - IT Center for Science Ltd., Finland 6! * 7! * This library is free software; you can redistribute it and/or 8! * modify it under the terms of the GNU Lesser General Public 9! * License as published by the Free Software Foundation; either 10! * version 2.1 of the License, or (at your option) any later version. 11! * 12! * This library is distributed in the hope that it will be useful, 13! * but WITHOUT ANY WARRANTY; without even the implied warranty of 14! * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU 15! * Lesser General Public License for more details. 16! * 17! * You should have received a copy of the GNU Lesser General Public 18! * License along with this library (in file ../LGPL-2.1); if not, write 19! * to the Free Software Foundation, Inc., 51 Franklin Street, 20! * Fifth Floor, Boston, MA 02110-1301 USA 21! * 22! *****************************************************************************/ 23! 24!/****************************************************************************** 25! * 26! * Authors: Juha Ruokolainen 27! * Email: Juha.Ruokolainen@csc.fi 28! * Web: http://www.csc.fi/elmer 29! * Address: CSC - IT Center for Science Ltd. 30! * Keilaranta 14 31! * 02101 Espoo, Finland 32! * 33! * Original Date: 01 Oct 1996 34! * 35! *****************************************************************************/ 36 37#include "../config.h" 38 39!> \ingroup ElmerLib 40!> \} 41 42!----------------------------------------------------------------------------- 43!> Miscellaneous utilities. 44!----------------------------------------------------------------------------- 45MODULE GeneralUtils 46 47USE Types 48USE LoadMod 49 50#ifdef HAVE_LUA 51USE, INTRINSIC :: ISO_C_BINDING 52#endif 53 54IMPLICIT NONE 55 56INTERFACE AllocateVector 57 MODULE PROCEDURE AllocateRealVector, AllocateIntegerVector, & 58 AllocateComplexVector, AllocateLogicalVector, & 59 AllocateElementVector 60END INTERFACE 61 62INTERFACE AllocateArray 63 MODULE PROCEDURE AllocateRealArray, AllocateIntegerArray, & 64 AllocateComplexArray, AllocateLogicalArray 65END INTERFACE 66 67INTERFACE ComponentName 68 MODULE PROCEDURE ComponentNameStr, ComponentNameVar 69END INTERFACE 70 71 REAL(KIND=dp), PRIVATE :: AdvanceTime1, AdvanceTime2 72 73CONTAINS 74 75!------------------------------------------------------------------------------ 76 SUBROUTINE StartAdvanceOutput( SolverName, OutputType ) 77!------------------------------------------------------------------------------ 78 CHARACTER(LEN=*) :: SolverName, OutputType 79!------------------------------------------------------------------------------ 80 AdvanceTime1 = RealTime() 81 AdvanceTime2 = RealTime() 82 CALL Info( SolverName, OutputType, Level=5 ) 83!------------------------------------------------------------------------------ 84 END SUBROUTINE StartAdvanceOutput 85!------------------------------------------------------------------------------ 86 87 88!------------------------------------------------------------------------------ 89 SUBROUTINE AdvanceOutput(t,n,dot_t,percent_t) 90!------------------------------------------------------------------------------ 91 IMPLICIT NONE 92 INTEGER :: t,n 93 REAL(KIND=dp), OPTIONAL :: dot_t,percent_t 94!------------------------------------------------------------------------------ 95 INTEGER :: i 96 REAL(KIND=dp) :: d_t, p_t 97!------------------------------------------------------------------------------ 98 d_t = 1._dp 99 p_t = 20._dp 100 IF ( PRESENT(dot_t) ) d_t = dot_t 101 IF ( PRESENT(percent_t) ) p_t = percent_t 102 103 IF ( RealTime() - AdvanceTime1 > d_t ) THEN 104 CALL Info( '', '.', Level=5, noAdvance=.TRUE. ) 105 106 IF ( RealTime() - AdvanceTime2 > p_t ) THEN 107 i = NINT(t*100.0/n) 108 WRITE(Message, '(i3,a)' ) i, '%' 109 CALL Info( '', Message, Level=5 ) 110 AdvanceTime2 = RealTime() 111 END IF 112 AdvanceTime1 = RealTime() 113 END IF 114!------------------------------------------------------------------------------ 115 END SUBROUTINE AdvanceOutput 116!------------------------------------------------------------------------------ 117 118 119!------------------------------------------------------------------------------ 120 PURE FUNCTION lentrim(str) RESULT(n) 121!------------------------------------------------------------------------------ 122 CHARACTER(LEN=*), INTENT(IN) :: str 123 INTEGER :: n 124 DO n=LEN(str),1,-1 125 IF ( str(n:n) /= ' ' ) EXIT 126 END DO 127!------------------------------------------------------------------------------ 128 END FUNCTION lentrim 129!------------------------------------------------------------------------------ 130 131 132!------------------------------------------------------------------------------ 133!> Compare equality of start of s1 to (in most uses string literal) s2. 134!------------------------------------------------------------------------------ 135!------------------------------------------------------------------------------ 136 PURE FUNCTION SEQL(s1,s2) RESULT(L) 137!------------------------------------------------------------------------------ 138 LOGICAL :: L 139 CHARACTER(LEN=*), INTENT(IN) :: s1,s2 140!------------------------------------------------------------------------------ 141 INTEGER :: n 142!------------------------------------------------------------------------------ 143 L = .FALSE. 144 n = LEN(s2) 145 IF(LEN(s1) < n) RETURN 146 IF (s1(1:n)==s2) L=.TRUE. 147!------------------------------------------------------------------------------ 148 END FUNCTION SEQL 149!------------------------------------------------------------------------------ 150 151 152!------------------------------------------------------------------------------ 153!> Converts integer to string. Handy when writing output with integer data. 154!------------------------------------------------------------------------------ 155 PURE FUNCTION i2s(ival) RESULT(str) 156!------------------------------------------------------------------------------ 157 INTEGER, INTENT(in) :: ival 158 CHARACTER(LEN=12) :: str 159!------------------------------------------------------------------------------ 160 INTEGER :: i,j,n,m,t,v 161 CHARACTER, PARAMETER :: DIGITS(0:9)=['0','1','2','3','4','5','6','7','8','9'] 162!------------------------------------------------------------------------------ 163 str = ' ' 164 165 IF ( ival >= 0 ) THEN 166 j=0 167 v=ival 168 ELSE 169 str(1:1)='-' 170 j=1 171 v=-ival 172 END IF 173 174 IF (v<10) THEN 175 str(j+1:j+1)=DIGITS(v) 176 ELSE 177 n=2 178 m=10 179 DO WHILE(10*m<=v) 180 n=n+1 181 m=m*10 182 END DO 183 184 DO i=j+1,j+n 185 t = v / m 186 str(i:i) = DIGITS(t) 187 v = v - t*m 188 m = m / 10 189 END DO 190 END IF 191!------------------------------------------------------------------------------ 192 END FUNCTION i2s 193!------------------------------------------------------------------------------ 194 195 196!------------------------------------------------------------------------------ 197!> Converts string of length n to an integer number. 198!------------------------------------------------------------------------------ 199 PURE FUNCTION s2i(str,n) RESULT(ival) 200!------------------------------------------------------------------------------ 201 INTEGER, INTENT(IN) :: n 202 INTEGER :: ival 203 CHARACTER(LEN=n), INTENT(IN) :: str 204!------------------------------------------------------------------------------ 205 LOGICAL :: neg 206 INTEGER :: j,k 207 INTEGER, PARAMETER :: ic0 = ICHAR('0') 208 209 neg = str(1:1)=='-' 210 k=1 211 IF ( neg ) k=2 212 213 ival = 0 214 DO j=k,n 215 ival = 10*ival + ICHAR(str(j:j)) - ic0 216 END DO 217 IF(neg) ival=-ival 218 END FUNCTION s2i 219!------------------------------------------------------------------------------ 220 221 222!------------------------------------------------------------------------------ 223!> Converts a string into a number of integer numbers 224!> It is assumed that the integers may also be separated by 225!> the given separator. 226!------------------------------------------------------------------------------ 227 FUNCTION str2ints(str,ints,sep) RESULT(n) 228!------------------------------------------------------------------------------ 229 INTEGER, INTENT(out) :: ints(:) 230 CHARACTER(LEN=*), INTENT(in) :: str 231 CHARACTER, OPTIONAL, INTENT(in) :: sep 232 233 INTEGER :: i,k,l,m,n,ic, icsep 234 INTEGER, PARAMETER :: ic0 = ICHAR('0'), ic9 = ICHAR('9'), icm = ICHAR('-'), & 235 ics = ICHAR(' ') 236 237 IF( PRESENT( sep ) ) THEN 238 icsep = ICHAR(sep) 239 ELSE 240 icsep = ics 241 END IF 242 243 244 k = LEN_TRIM(str) 245 l = 1 246 n = 0 247 DO WHILE(l<=k.AND.n<SIZE(ints)) 248 DO WHILE(l<=k) 249 ic = ICHAR(str(l:l)) 250 IF( ic == ics .OR. ic == icsep ) THEN 251 CONTINUE 252 ELSE 253 EXIT 254 END IF 255 l=l+1 256 END DO 257 IF(l>k) EXIT 258 IF(.NOT.(ic==icm .OR. ic>=ic0 .AND. ic<=ic9)) EXIT 259 260 m = l+1 261 DO WHILE(m<=k) 262 ic = ICHAR(str(m:m)) 263 IF(ic<ic0 .OR. ic>ic9) EXIT 264 m=m+1 265 END DO 266 267 n = n + 1 268 ints(n) = s2i(str(l:m-1),m-l) 269 l = m 270 END DO 271 END FUNCTION str2ints 272!------------------------------------------------------------------------------ 273 274 275 276 277!------------------------------------------------------------------------------ 278 SUBROUTINE SystemCommand( cmd ) 279!------------------------------------------------------------------------------ 280 CHARACTER(LEN=*) :: cmd 281 CALL SystemC( TRIM(cmd) // CHAR(0) ) 282!------------------------------------------------------------------------------ 283 END SUBROUTINE SystemCommand 284!------------------------------------------------------------------------------ 285 286 287!------------------------------------------------------------------------------ 288 PURE FUNCTION FileNameQualified(file) RESULT(L) 289!------------------------------------------------------------------------------ 290 LOGICAL :: L 291 CHARACTER(*), INTENT(IN) :: file 292!------------------------------------------------------------------------------ 293 L = INDEX(file,':')>0 .OR. file(1:1)=='/' .OR. file(1:1)==Backslash 294!------------------------------------------------------------------------------ 295 END FUNCTION FileNameQualified 296!------------------------------------------------------------------------------ 297 298!------------------------------------------------------------------------------ 299 PURE FUNCTION LittleEndian() RESULT(L) 300!------------------------------------------------------------------------------ 301 LOGICAL :: L 302!------------------------------------------------------------------------------ 303 INTEGER(1) :: s(2) 304 INTEGER(2), PARAMETER :: t = 256*7+8 305!------------------------------------------------------------------------------ 306 s = TRANSFER(t,s) 307 L=s(1)==8 308!------------------------------------------------------------------------------ 309 END FUNCTION LittleEndian 310!------------------------------------------------------------------------------ 311 312!------------------------------------------------------------------------------ 313 FUNCTION FormatDate() RESULT( date ) 314!------------------------------------------------------------------------------ 315 CHARACTER( LEN=20 ) :: date 316 INTEGER :: dates(8) 317 318 CALL DATE_AND_TIME( VALUES=dates ) 319 WRITE( date, & 320 '(I4,"/",I2.2,"/",I2.2," ",I2.2,":",I2.2,":",I2.2)' ) & 321 dates(1),dates(2),dates(3),dates(5),dates(6),dates(7) 322!------------------------------------------------------------------------------ 323 END FUNCTION FormatDate 324!------------------------------------------------------------------------------ 325 326 327!------------------------------------------------------------------------------ 328!> Sort an array of integer values. 329!------------------------------------------------------------------------------ 330 PURE SUBROUTINE Sort( n,a ) 331!------------------------------------------------------------------------------ 332 INTEGER, INTENT(in) :: n 333 INTEGER, INTENT(inout) :: a(:) 334!------------------------------------------------------------------------------ 335 336 INTEGER :: i,j,l,ir,ra 337!------------------------------------------------------------------------------ 338 339 IF ( n <= 1 ) RETURN 340 341 l = n / 2 + 1 342 ir = n 343 DO WHILE( .TRUE. ) 344 IF ( l > 1 ) THEN 345 l = l - 1 346 ra = a(l) 347 ELSE 348 ra = a(ir) 349 a(ir) = a(1) 350 ir = ir - 1 351 IF ( ir == 1 ) THEN 352 a(1) = ra 353 RETURN 354 END IF 355 END IF 356 i = l 357 j = l + l 358 DO WHILE( j <= ir ) 359 IF ( j<ir ) THEN 360 IF ( a(j)<a(j+1) ) j = j+1 361 END IF 362 363 IF ( ra<a(j) ) THEN 364 a(i) = a(j) 365 i = j 366 j = j + i 367 ELSE 368 j = ir + 1 369 END IF 370 a(i) = ra 371 END DO 372 END DO 373 374!------------------------------------------------------------------------------ 375 END SUBROUTINE Sort 376!------------------------------------------------------------------------------ 377 378!------------------------------------------------------------------------------ 379!> Sort an integer array a, together with an another integer array. 380!------------------------------------------------------------------------------ 381 PURE SUBROUTINE SortI( n,a,b ) 382!------------------------------------------------------------------------------ 383 INTEGER, INTENT(in) :: n 384 INTEGER, INTENT(inout) :: a(:),b(:) 385!------------------------------------------------------------------------------ 386 387 INTEGER :: i,j,l,ir,ra,rb 388!------------------------------------------------------------------------------ 389 390 IF ( n <= 1 ) RETURN 391 392 l = n / 2 + 1 393 ir = n 394 DO WHILE( .TRUE. ) 395 IF ( l > 1 ) THEN 396 l = l - 1 397 ra = a(l) 398 rb = b(l) 399 ELSE 400 ra = a(ir) 401 rb = b(ir) 402 a(ir) = a(1) 403 b(ir) = b(1) 404 ir = ir - 1 405 IF ( ir == 1 ) THEN 406 a(1) = ra 407 b(1) = rb 408 RETURN 409 END IF 410 END IF 411 i = l 412 j = l + l 413 DO WHILE( j <= ir ) 414 IF ( j<ir ) THEN 415 IF ( a(j)<a(j+1) ) j = j+1 416 END IF 417 IF ( ra<a(j) ) THEN 418 a(i) = a(j) 419 b(i) = b(j) 420 i = j 421 j = j + i 422 ELSE 423 j = ir + 1 424 END IF 425 a(i) = ra 426 b(i) = rb 427 END DO 428 END DO 429 430!------------------------------------------------------------------------------ 431 END SUBROUTINE SortI 432!------------------------------------------------------------------------------ 433 434 435!------------------------------------------------------------------------------ 436!> Sort an index array, and change the order of an real array accordingly. 437!------------------------------------------------------------------------------ 438 PURE SUBROUTINE SortF( n,a,b ) 439!------------------------------------------------------------------------------ 440 INTEGER, INTENT(in) :: n 441 INTEGER, INTENT(inout) :: a(:) 442 REAL(KIND=dp), INTENT(inout) :: b(:) 443!------------------------------------------------------------------------------ 444 445 INTEGER :: i,j,l,ir,ra 446 REAL(KIND=dp) :: rb 447!------------------------------------------------------------------------------ 448 449 IF ( n <= 1 ) RETURN 450 451 l = n / 2 + 1 452 ir = n 453 DO WHILE( .TRUE. ) 454 455 IF ( l > 1 ) THEN 456 l = l - 1 457 ra = a(l) 458 rb = b(l) 459 ELSE 460 ra = a(ir) 461 rb = b(ir) 462 a(ir) = a(1) 463 b(ir) = b(1) 464 ir = ir - 1 465 IF ( ir == 1 ) THEN 466 a(1) = ra 467 b(1) = rb 468 RETURN 469 END IF 470 END IF 471 i = l 472 j = l + l 473 DO WHILE( j <= ir ) 474 IF ( j<ir ) THEN 475 IF ( a(j)<a(j+1) ) j = j+1 476 END IF 477 IF ( ra<a(j) ) THEN 478 a(i) = a(j) 479 b(i) = b(j) 480 i = j 481 j = j + i 482 ELSE 483 j = ir + 1 484 END IF 485 a(i) = ra 486 b(i) = rb 487 END DO 488 END DO 489 490!------------------------------------------------------------------------------ 491 END SUBROUTINE SortF 492!------------------------------------------------------------------------------ 493 494 495!------------------------------------------------------------------------------ 496!> Sort an real array, and change the order of an index array accordingly. 497!------------------------------------------------------------------------------ 498 PURE SUBROUTINE SortD( n,a,b ) 499!------------------------------------------------------------------------------ 500 INTEGER, INTENT(in) :: n 501 INTEGER, INTENT(inout) :: b(:) 502 REAL(KIND=dp), INTENT(inout) :: a(:) 503!------------------------------------------------------------------------------ 504 505 INTEGER :: i,j,l,ir,rb 506 REAL(KIND=dp) :: ra 507!------------------------------------------------------------------------------ 508 509 IF ( n <= 1 ) RETURN 510 511 l = n / 2 + 1 512 ir = n 513 DO WHILE( .TRUE. ) 514 515 IF ( l > 1 ) THEN 516 l = l - 1 517 ra = a(l) 518 rb = b(l) 519 ELSE 520 ra = a(ir) 521 rb = b(ir) 522 a(ir) = a(1) 523 b(ir) = b(1) 524 ir = ir - 1 525 IF ( ir == 1 ) THEN 526 a(1) = ra 527 b(1) = rb 528 RETURN 529 END IF 530 END IF 531 i = l 532 j = l + l 533 DO WHILE( j <= ir ) 534 IF ( j<ir ) THEN 535 IF ( a(j)<a(j+1) ) j = j+1 536 END IF 537 IF ( ra<a(j) ) THEN 538 a(i) = a(j) 539 b(i) = b(j) 540 i = j 541 j = j + i 542 ELSE 543 j = ir + 1 544 END IF 545 a(i) = ra 546 b(i) = rb 547 END DO 548 END DO 549 550!------------------------------------------------------------------------------ 551 END SUBROUTINE SortD 552!------------------------------------------------------------------------------ 553 554 555!------------------------------------------------------------------------------ 556!> Sort an complex array, and organize an index table accordingly. 557!------------------------------------------------------------------------------ 558 PURE SUBROUTINE SortC( n,a,b ) 559!------------------------------------------------------------------------------ 560 INTEGER, INTENT(in) :: n 561 INTEGER, INTENT(inout) :: b(:) 562 COMPLEX(KIND=dp), INTENT(inout) :: a(:) 563!------------------------------------------------------------------------------ 564 565 INTEGER :: i,j,l,ir,rb 566 COMPLEX(KIND=dp) :: ra 567!------------------------------------------------------------------------------ 568 569 IF ( n <= 1 ) RETURN 570 571 l = n / 2 + 1 572 ir = n 573 DO WHILE( .TRUE. ) 574 IF ( l > 1 ) THEN 575 l = l - 1 576 ra = a(l) 577 rb = b(l) 578 ELSE 579 ra = a(ir) 580 rb = b(ir) 581 a(ir) = a(1) 582 b(ir) = b(1) 583 ir = ir - 1 584 IF ( ir == 1 ) THEN 585 a(1) = ra 586 b(1) = rb 587 RETURN 588 END IF 589 END IF 590 i = l 591 j = l + l 592 DO WHILE( j <= ir ) 593 IF ( j<ir ) THEN 594 IF ( ABS(a(j))<ABS(a(j+1)) ) j = j+1 595 END IF 596 IF ( ABS(ra)<ABS(a(j)) ) THEN 597 a(i) = a(j) 598 b(i) = b(j) 599 i = j 600 j = j + i 601 ELSE 602 j = ir + 1 603 END IF 604 a(i) = ra 605 b(i) = rb 606 END DO 607 END DO 608 609!------------------------------------------------------------------------------ 610 END SUBROUTINE SortC 611!------------------------------------------------------------------------------ 612 613 614!------------------------------------------------------------------------------ 615!> Order real components in b in a decreasing order and return the new order 616!> of indexes in a. 617!------------------------------------------------------------------------------ 618 PURE SUBROUTINE SortR( n,a,b ) 619!------------------------------------------------------------------------------ 620 INTEGER, INTENT(in) :: n 621 INTEGER, INTENT(inout) :: a(:) 622 REAL(KIND=dp), INTENT(inout) :: b(:) 623!------------------------------------------------------------------------------ 624 625 INTEGER :: i,j,l,ir,ra 626 REAL(KIND=dp) :: rb 627!------------------------------------------------------------------------------ 628 629 IF ( n <= 1 ) RETURN 630 631 l = n / 2 + 1 632 ir = n 633 DO WHILE( .TRUE. ) 634 635 IF ( l > 1 ) THEN 636 l = l - 1 637 ra = a(l) 638 rb = b(l) 639 ELSE 640 ra = a(ir) 641 rb = b(ir) 642 a(ir) = a(1) 643 b(ir) = b(1) 644 ir = ir - 1 645 IF ( ir == 1 ) THEN 646 a(1) = ra 647 b(1) = rb 648 RETURN 649 END IF 650 END IF 651 i = l 652 j = l + l 653 DO WHILE( j <= ir ) 654 IF ( j<ir ) THEN 655 IF ( b(j) > b(j+1) ) j = j+1 656 END IF 657 IF ( rb > b(j) ) THEN 658 a(i) = a(j) 659 b(i) = b(j) 660 i = j 661 j = j + i 662 ELSE 663 j = ir + 1 664 END IF 665 a(i) = ra 666 b(i) = rb 667 END DO 668 END DO 669 670!------------------------------------------------------------------------------ 671 END SUBROUTINE SortR 672!------------------------------------------------------------------------------ 673 674!------------------------------------------------------------------------------ 675!> Search an integer value in an ordered array. 676!------------------------------------------------------------------------------ 677 PURE FUNCTION SearchI( N,Array,Val ) RESULT ( Idx ) 678!------------------------------------------------------------------------------ 679 INTEGER, INTENT(in) :: N,Val,Array(:) 680!------------------------------------------------------------------------------ 681 INTEGER :: Lower, Upper,Lou,Idx 682!------------------------------------------------------------------------------ 683 684 Idx = 0 685 Upper = N 686 Lower = 1 687 688 ! Handle the special case 689 690 IF ( Upper == 0 ) RETURN 691 692 DO WHILE( .TRUE. ) 693 IF ( Array(Lower) == Val) THEN 694 Idx = Lower 695 EXIT 696 ELSE IF ( Array(Upper) == Val ) THEN 697 Idx = Upper 698 EXIT 699 END IF 700 701 IF ( (Upper-Lower)>1 ) THEN 702 Lou = ISHFT((Upper + Lower), -1) 703 IF ( Array(Lou) < Val ) THEN 704 Lower = Lou 705 ELSE 706 Upper = Lou 707 END IF 708 ELSE 709 EXIT 710 END IF 711 END DO 712 713 RETURN 714 715!------------------------------------------------------------------------------ 716 END FUNCTION SearchI 717!------------------------------------------------------------------------------ 718 719 720 721!------------------------------------------------------------------------------ 722!> Search a real value in an ordered array. 723!------------------------------------------------------------------------------ 724 PURE FUNCTION SearchR( N,Array,Val ) RESULT ( Idx ) 725!------------------------------------------------------------------------------ 726 727 INTEGER, INTENT(in) :: N 728 INTEGER :: Idx 729 REAL(KIND=dP), INTENT(in) :: Val,Array(:) 730!------------------------------------------------------------------------------ 731 INTEGER :: Lower, Upper,Lou 732!------------------------------------------------------------------------------ 733 734 Idx = 0 735 Upper = N 736 Lower = 1 737 738 ! Handle the special case 739 IF ( Upper == 0 ) RETURN 740 741 DO WHILE( .TRUE. ) 742 IF ( ABS( Array(Lower) - Val) < TINY(Val) ) THEN 743 Idx = Lower 744 EXIT 745 ELSE IF ( ABS( Array(Upper) - Val ) < TINY(Val) ) THEN 746 Idx = Upper 747 EXIT 748 END IF 749 750 IF ( (Upper-Lower) > 1 ) THEN 751 Lou = ISHFT((Upper + Lower), -1) 752 IF ( Array(Lou) < Val ) THEN 753 Lower = Lou 754 ELSE 755 Upper = Lou 756 END IF 757 ELSE 758 EXIT 759 END IF 760 END DO 761 762 RETURN 763 764!------------------------------------------------------------------------------ 765 END FUNCTION SearchR 766!------------------------------------------------------------------------------ 767 768 769!------------------------------------------------------------------------------ 770 SUBROUTINE OpenIncludeFile( Unit, FileName, IncludePath ) 771!------------------------------------------------------------------------------ 772 INTEGER :: Unit 773 CHARACTER(LEN=*) :: FileName, IncludePath 774!------------------------------------------------------------------------------ 775 INTEGER :: i,j,k,k0,k1,l,iostat 776 CHARACTER(LEN=1024) :: name, TmpName 777!------------------------------------------------------------------------------ 778 779 i = 1 780 name = FileName 781 DO WHILE( name(i:i) == ' ' .OR. name(i:i)=='"') 782 i = i + 1 783 END DO 784 j = LEN_TRIM(name) 785 IF ( name(j:j) == '"' ) j=j-1 786 name = TRIM(name(i:j)) 787 788 IF ( INDEX(name,':') == 0 .AND. name(1:1) /= '/' .AND. & 789 name(1:1) /= Backslash ) THEN 790 k0 = 1 791 DO WHILE( IncludePath(k0:k0) == '"' ) 792 k0 = k0+1 793 END DO 794 k1 = INDEX( IncludePath, ';' ) 795 796 DO WHILE( k1 >= k0 ) 797 DO k = k1-1,k0,-1 798 IF ( IncludePath(k:k) /= ' ' .AND. IncludePath(k:k)/='"' ) EXIT 799 END DO 800 IF ( IncludePath(k:k) == '"' ) k=k-1 801 IF ( k >= k0 ) THEN 802 WRITE( tmpName,'(a,a,a)' ) IncludePath(k0:k), '/', TRIM(name) 803 OPEN( Unit, FILE=TRIM(tmpName), STATUS='OLD',ERR=10 ) 804 RETURN 805 END IF 80610 CONTINUE 807 k0 = k1+1 808 k1 = INDEX( IncludePath(k0:), ';' ) + k0 - 1 809 END DO 810 811 IF ( LEN_TRIM(IncludePath(k0:))>0 ) THEN 812 k1 = INDEX( IncludePath(k0:), '"' ) + k0 - 2 813 IF ( k1 < k0 ) k1=LEN_TRIM(IncludePath) 814 tmpName = TRIM(IncludePath(k0:k1)) // '/' // TRIM(name) 815 OPEN( Unit, FILE=TRIM(TmpName), STATUS='OLD',ERR=20 ) 816 RETURN 817 END IF 818 81920 CONTINUE 820 OPEN( Unit, FILE=TRIM(name), STATUS='OLD',IOSTAT=iostat ) 821 ELSE 822 OPEN( Unit, FILE=TRIM(name), STATUS='OLD',IOSTAT=iostat ) 823 END IF 824 825 IF( iostat /= 0 ) THEN 826 CALL Fatal('OpenIncludeFile','Cannot open include file: '//TRIM(Name)) 827 END IF 828 829 830!------------------------------------------------------------------------------ 831 END SUBROUTINE OpenIncludeFile 832!------------------------------------------------------------------------------ 833 834 835!------------------------------------------------------------------------------ 836!> Read a (logical) line from FORTRAN device Unit and remove leading, trailing, 837!> and multiple blanks between words. Also convert uppercase characters to 838!> lowercase.The logical line can continue the several physical lines by adding 839!> the backslash (\) mark at the end of a physical line. 840!------------------------------------------------------------------------------ 841 RECURSIVE FUNCTION ReadAndTrim( Unit,str,echo,literal,noeval ) RESULT(l) 842!------------------------------------------------------------------------------ 843 INTEGER :: Unit !< Fortran unit number to read from 844 CHARACTER(LEN=:), ALLOCATABLE :: str !< The string read from the file 845 LOGICAL, OPTIONAL :: Echo 846 LOGICAL, OPTIONAL :: literal 847 LOGICAL, OPTIONAL :: noeval 848 LOGICAL :: l !< Success of the read operation 849!------------------------------------------------------------------------------ 850 INTEGER, PARAMETER :: MAXLEN = 16384 851 852 CHARACTER(LEN=:), ALLOCATABLE :: temp 853 CHARACTER(LEN=12) :: tmpstr 854 CHARACTER(LEN=MAXLEN) :: readstr = ' ', copystr = ' ', matcstr=' ' , IncludePath=' ' 855 856 LOGICAL :: InsideQuotes, OpenSection=.FALSE., DoEval 857 INTEGER :: i,j,k,m,ValueStarts=0,inlen,ninlen,outlen,IncludeUnit=28,IncludeUnitBase=28 858 859 CHARACTER(LEN=MAX_NAME_LEN) :: Prefix = ' ' 860 861 INTEGER, PARAMETER :: A=ICHAR('A'),Z=ICHAR('Z'),U2L=ICHAR('a')-ICHAR('A'),Tab=9 862 CHARACTER(LEN=MAXLEN) :: tmatcstr, tcmdstr 863 INTEGER :: tninlen 864 865 SAVE ReadStr, ValueStarts, Prefix, OpenSection 866 867 IF ( PRESENT(literal) ) literal=.FALSE. 868 l = .TRUE. 869 870 ! Optionally do not expand the MATC and LUA expressions. 871 DoEval = .TRUE. 872 IF( PRESENT( NoEval ) ) THEN 873 DoEval = .NOT. NoEval 874 END IF 875 876 IF(.NOT.ALLOCATED(str)) ALLOCATE(CHARACTER(512)::str) 877 outlen = LEN(str) 878 879 IF ( ValueStarts==0 .AND. OpenSection ) THEN 880 str = 'end' 881 ValueStarts = 0 882 OpenSection = .FALSE. 883 RETURN 884 END IF 885 886 IF ( ValueStarts == 0 ) THEN 887 tmpstr = ' ' 888 DO WHILE( .TRUE. ) 889 IF ( IncludeUnit < IncludeUnitBase ) THEN 890 READ( IncludeUnit,'(A)',END=1,ERR=1 ) readstr 891 GO TO 2 8921 CLOSE(IncludeUnit) 893 IncludeUnit = IncludeUnit+1 894 READ( Unit,'(A)',END=10,ERR=10 ) readstr 8952 CONTINUE 896 ELSE 897 READ( Unit,'(A)',END=10,ERR=10 ) readstr 898 END IF 899 900 readstr = ADJUSTL(readstr) 901 902 DO k=1,12 903 j = ICHAR(readstr(k:k)) 904 IF ( j >= A .AND. j<= Z ) THEN 905 Tmpstr(k:k) = CHAR(j+U2L) 906 ELSE 907 tmpstr(k:k) = readstr(k:k) 908 END IF 909 END DO 910 911 IF ( SEQL(Tmpstr, 'include path') ) THEN 912 k = LEN_TRIM(readstr) 913 IncludePath(1:k-13) = readstr(14:k) 914 Tmpstr = '' 915 ELSE 916 EXIT 917 END IF 918 END DO 919 920 IF ( SEQL(tmpstr, 'include ') ) THEN 921 IncludeUnit = IncludeUnit-1 922 923 CALL Info('OpenIncludeFile','Trying to include file: '//TRIM(readstr(9:)),Level=10) 924 925 inlen = LEN_TRIM(readstr) 926 927 i = INDEX( readstr(1:inlen), '$' ) 928 IF ( i>0 .AND. i<inlen ) THEN 929 CALL TrimMatcExpression() 930 CALL Info('ReadAndTrim','Include file after MATC trimming: '& 931 //TRIM(readstr(9:)),Level=6) 932 END IF 933 934 i = INDEX( readstr(1:inlen),'#') 935 IF( i>0 .AND. i<inlen ) THEN 936#ifdef HAVE_LUA 937 CALL TrimLuaExpression() 938 CALL Info('ReadAndTrim','Include file after LUA trimming: '& 939 //TRIM(readstr(9:)),Level=6) 940#else 941 CALL Fatal('ReadAndTrim','LUA not included, cannot continue') 942#endif 943 END IF 944 945 CALL OpenIncludeFile( IncludeUnit, TRIM(readstr(9:)), IncludePath ) 946 947 READ( IncludeUnit,'(A)',END=3,ERR=3 ) readstr 948 GO TO 4 9493 CLOSE(IncludeUnit) 950 IncludeUnit = IncludeUnit+1 951 READ( Unit,'(A)',END=10,ERR=10 ) readstr 9524 CONTINUE 953 END IF 954 ninlen = LEN_TRIM(readstr) 955 ELSE 956 inlen = LEN_TRIM(readstr) 957 ninlen = inlen-ValueStarts+1 958 IF ( Prefix == ' ' ) THEN 959 readstr = readstr(ValueStarts:inlen) 960 ELSE IF ( Prefix == '::' ) THEN 961 readstr = readstr(ValueStarts:inlen) 962 OpenSection = .TRUE. 963 Prefix = ' ' 964 ELSE 965 DO i=ValueStarts,inlen 966 IF ( readstr(i:i) == ')' ) THEN 967 readstr(i:i) = ' ' 968 EXIT 969 ELSE IF ( readstr(i:i) == ',' ) THEN 970 readstr(i:i) = ' ' 971 END IF 972 END DO 973 ninlen = ninlen + LEN_TRIM(Prefix) + 1 974 readstr = TRIM(Prefix) // ' ' // readstr(ValueStarts:inlen) 975 END IF 976 END IF 977 978 ValueStarts = 0 979 InsideQuotes = .FALSE. 980 981 i = INDEX( readstr(1:ninlen), '!' ) 982 IF ( i>0 ) ninlen=i-1 983 984 i = 1 985 inlen = ninlen 986 DO WHILE( i <= inlen ) 987 IF ( readstr(i:i) == '"' ) InsideQuotes = .NOT.InsideQuotes 988 IF ( .NOT. InsideQuotes .AND. readstr(i:i) == CHAR(92) .AND. i==inlen ) THEN 989 readstr(i:i) = ' ' 990 IF ( IncludeUnit < IncludeUnitBase ) THEN 991 READ( IncludeUnit,'(A)',END=10,ERR=10 ) readstr(i+1:MAXLEN) 992 ELSE 993 READ( Unit,'(A)',END=10,ERR=10 ) readstr(i+1:MAXLEN) 994 END IF 995 DO j=LEN(readstr),i+1,-1 996 IF ( readstr(j:j) /= ' ' ) EXIT 997 END DO 998 inlen = inlen + j-i 999 END IF 1000 i = i + 1 1001 END DO 1002 1003 i = INDEX( readstr(1:inlen), '#' ) 1004 IF ( i>0 .AND. i<inlen .AND. DoEval ) THEN 1005#ifdef HAVE_LUA 1006 CALL TrimLuaExpression() 1007#else 1008 CALL Fatal('ReadAndTrim','LUA not included, cannot continue') 1009#endif 1010 END IF 1011 1012 i = INDEX( readstr(1:inlen), '$' ) 1013 IF ( i>0 .AND. i<inlen .AND. DoEval ) THEN 1014 CALL TrimMatcExpression() 1015 END IF 1016 1017 IF ( PRESENT( Echo ) ) THEN 1018 IF ( Echo ) WRITE( 6, '(a)' ) readstr(1:inlen) 1019 END IF 1020 1021 i = 1 1022 DO WHILE(i <= inlen ) 1023 IF (readstr(i:i) /= ' ' .AND. ICHAR(readstr(i:i))/=Tab ) EXIT 1024 i = i + 1 1025 END DO 1026 1027 InsideQuotes = .FALSE. 1028 1029 IF ( PRESENT(literal) ) THEN 1030 IF ( readstr(i:i) == '"' ) literal=.TRUE. 1031 END IF 1032 1033 k = 1 1034 DO WHILE( i<=inlen ) 1035 IF ( readstr(i:i) == '"' ) THEN 1036 InsideQuotes = .NOT.InsideQuotes 1037 i=i+1 1038 IF ( i>inlen ) EXIT 1039 END IF 1040 1041 IF ( .NOT.InsideQuotes ) THEN 1042 IF ( readstr(i:i) == '!' .OR. readstr(i:i) == '#' .OR. & 1043 readstr(i:i) == '=' .OR. readstr(i:i) == '(' .OR. & 1044 readstr(i:i) == ';' .OR. readstr(i:i+1) == '::' ) EXIT 1045 IF (ICHAR( readstr(i:i))<32.AND.ICHAR(readstr(i:i))/=Tab) EXIT 1046 END IF 1047 1048 DO WHILE( i <= inlen ) 1049 IF ( readstr(i:i) == '"' ) THEN 1050 InsideQuotes = .NOT.InsideQuotes 1051 i=i+1 1052 IF ( i>inlen ) EXIT 1053 END IF 1054 1055 IF ( .NOT.InsideQuotes ) THEN 1056 IF ( readstr(i:i) == ' ' .OR. readstr(i:i) == '=' .OR. & 1057 readstr(i:i) == ';' .OR. readstr(i:i) == '(' .OR. & 1058 readstr(i:i+1) == '::' ) EXIT 1059 IF ( ICHAR( readstr(i:i))<32 ) EXIT 1060 END IF 1061 1062 IF ( k>outlen ) THEN 1063 temp = str 1064 DEALLOCATE(str) 1065 outlen=LEN(temp)+512 1066 ALLOCATE(CHARACTER(outlen)::str) 1067 str(1:LEN(temp))=temp; str(LEN(temp)+1:)='' 1068 DEALLOCATE(temp) 1069 END IF 1070 1071 j = ICHAR( readstr(i:i) ) 1072 IF ( .NOT.InsideQuotes .AND. j>=A .AND. j<=Z ) THEN 1073 str(k:k) = CHAR(j+U2L) 1074 ELSE IF ( .NOT.InsideQuotes .AND. j==Tab ) THEN 1075 str(k:k) = ' ' 1076 ELSE 1077 str(k:k) = readstr(i:i) 1078 ENDIF 1079 1080 i = i + 1 1081 k = k + 1 1082 END DO 1083 1084 IF ( k <= outlen ) str(k:k) = ' ' 1085 k = k + 1 1086 1087 DO WHILE( i<=inlen ) 1088 IF ( readstr(i:i) /= ' ' .AND. ICHAR(readstr(i:i))/=Tab ) EXIT 1089 i = i + 1 1090 END DO 1091 END DO 1092 str(k:)=' ' 1093 1094 IF ( i <= inlen ) THEN 1095 Prefix = ' ' 1096 IF ( ReadStr(i:i) == '=' ) THEN 1097 ValueStarts = i + 1 1098 ELSE IF ( ReadStr(i:i) == ';' ) THEN 1099 ValueStarts = i + 1 1100 ELSE IF ( ReadStr(i:i) == '(' ) THEN 1101 ValueStarts = i + 1 1102 Prefix = 'Size' 1103 ELSE IF ( ReadStr(i:i+1) == '::' ) THEN 1104 ValueStarts = i + 2 1105 Prefix = '::' 1106 ELSE IF ( ICHAR(readstr(i:i)) < 32 ) THEN 1107 DO WHILE( i <= inlen ) 1108 IF ( ICHAR(readstr(i:i)) >= 32 ) EXIT 1109 i = i + 1 1110 END DO 1111 IF ( i <= inlen ) ValueStarts = i 1112 END IF 1113 END IF 1114 RETURN 111510 CONTINUE 1116 l = .FALSE. 1117!------------------------------------------------------------------------------ 1118 1119 CONTAINS 1120 1121 SUBROUTINE TrimMatcExpression() 1122 1123 i = INDEX( readstr(1:inlen), '$' ) 1124 IF ( i>0 .AND. i<inlen ) THEN 1125 m = i 1126 copystr(i:inlen) = readstr(i:inlen) 1127 DO WHILE(i<=inlen) 1128 IF ( copystr(i:i) == '$' ) THEN 1129 DO j=i+1,inlen-1 1130 IF ( copystr(j:j) == '$' ) EXIT 1131 END DO 1132 ninlen = j - i 1133 1134 ! Initialize variables for each copy of MATC separately 1135 1136 !$OMP PARALLEL DEFAULT(NONE) & 1137 !$OMP SHARED(copystr, i, matcstr, ninlen, inlen) & 1138 !$OMP PRIVATE(tcmdstr, tmatcstr, tninlen) 1139 1140 tninlen = ninlen 1141 tcmdstr = copystr(i+1:inlen) 1142 CALL MATC( tcmdstr, tmatcstr, tninlen ) 1143 !$OMP BARRIER 1144 1145 !$OMP SINGLE 1146 matcstr(1:tninlen) = tmatcstr(1:tninlen) 1147 ninlen = tninlen 1148 !$OMP END SINGLE 1149 1150 !$OMP END PARALLEL 1151 1152 DO k=1,ninlen 1153 readstr(m:m) = matcstr(k:k) 1154 m = m + 1 1155 END DO 1156 i = j+1 1157 ELSE 1158 readstr(m:m) = copystr(i:i) 1159 i = i + 1 1160 m = m + 1 1161 END IF 1162 END DO 1163 IF ( m <= inlen ) readstr(m:inlen) = ' ' 1164 inlen = m-1 1165 END IF 1166 1167 END SUBROUTINE TrimMatcExpression 1168 1169#ifdef HAVE_LUA 1170 SUBROUTINE TrimLuaExpression() 1171 1172 INTEGER :: lstat 1173 character(kind=c_char, len=:), pointer :: lua_result 1174 integer :: result_len 1175 logical :: closed_region, first_bang 1176 closed_region = .false. 1177 first_bang = .true. 1178 m = i 1179 copystr(i:inlen) = readstr(i:inlen) 1180 DO WHILE(i<=inlen) 1181 IF ( copystr(i:i) == '#' ) THEN 1182 DO j=i+1,inlen-1 1183 IF ( copystr(j:j) == '#' ) EXIT 1184 END DO 1185 ninlen = j - i 1186 1187 ! Initialize variables for each copy of Lua interpreter separately 1188 1189 !$OMP PARALLEL DEFAULT(NONE) & 1190 !$OMP SHARED(copystr, i, matcstr, ninlen, inlen, closed_region, first_bang, j) & 1191 !$OMP PRIVATE(tcmdstr, tninlen, lstat, result_len, lua_result) 1192 1193 tninlen = ninlen 1194 tcmdstr = copystr(i+1:inlen) 1195 1196 IF(tcmdstr(tninlen:tninlen) == '#') then 1197 closed_region = .TRUE. 1198 ELSE 1199 closed_region = .FALSE. 1200 END IF 1201 1202 IF(closed_region) THEN 1203 lstat = lua_dostring( LuaState, & 1204 'return tostring('// tcmdstr(1:tninlen-1) // ')'//c_null_char, 1) 1205 ELSE 1206 IF (i == 1 .and. first_bang .and. j == inlen) THEN ! ' # <luacode>' case, do not do 'return tostring(..)'. 1207 ! Instead, just execute the line in the lua interpreter 1208 lstat = lua_dostring( LuaState, tcmdstr(1:tninlen) // c_null_char, 1) 1209 ELSE ! 'abc = # <luacode>' case, oneliners only 1210 lstat = lua_dostring( LuaState, & 1211 'return tostring('// tcmdstr(1:tninlen) // ')'//c_null_char, 1) 1212 END IF 1213 END IF 1214 lua_result => lua_popstring(LuaState, result_len) 1215 1216 !$OMP SINGLE 1217 matcstr(1:result_len) = lua_result(1:result_len) 1218 ninlen = result_len 1219 !$OMP END SINGLE 1220 1221 !$OMP END PARALLEL 1222 1223 DO k=1,ninlen 1224 readstr(m:m) = matcstr(k:k) 1225 m = m + 1 1226 END DO 1227 i = j+1 1228 ELSE 1229 readstr(m:m) = copystr(i:i) 1230 i = i + 1 1231 m = m + 1 1232 END IF 1233 first_bang = .false. 1234 END DO 1235 IF ( m <= inlen ) readstr(m:inlen) = ' ' 1236 inlen = m-1 1237 END SUBROUTINE TrimLuaExpression 1238#endif 1239 1240 END FUNCTION ReadAndTrim 1241!------------------------------------------------------------------------------ 1242 1243 1244!------------------------------------------------------------------------------ 1245 PURE FUNCTION GetVarName(Var) RESULT(str) 1246!------------------------------------------------------------------------------ 1247 TYPE(Variable_t), INTENT(in) :: Var 1248 CHARACTER(LEN=Var % NameLen) :: str 1249!------------------------------------------------------------------------------ 1250 str = Var % Name(1:Var % NameLen) 1251!------------------------------------------------------------------------------ 1252 END FUNCTION GetVarName 1253!------------------------------------------------------------------------------ 1254 1255 1256!------------------------------------------------------------------------------ 1257 FUNCTION ComponentNameVar( Var, Component ) RESULT(str) 1258!------------------------------------------------------------------------------ 1259 TYPE(Variable_t),INTENT(in) :: Var 1260 INTEGER, OPTIONAL,INTENT(in) :: Component 1261!------------------------------------------------------------------------------ 1262 CHARACTER(LEN=MAX_NAME_LEN) :: str 1263!------------------------------------------------------------------------------ 1264 IF ( Var % Name(1:Var % NameLen) == 'flow solution' ) THEN 1265 str='flow solution' 1266 IF ( .NOT. PRESENT(Component) ) RETURN 1267 IF ( Component == Var % DOFs ) THEN 1268 str = 'pressure' 1269 RETURN 1270 ELSE 1271 str = 'velocity ' // TRIM(i2s(Component)) 1272 END IF 1273 ELSE 1274 str = ComponentName(Var % Name, Component) 1275 END IF 1276!------------------------------------------------------------------------------ 1277END FUNCTION ComponentNameVar 1278!------------------------------------------------------------------------------ 1279 1280 1281!------------------------------------------------------------------------------ 1282 FUNCTION ComponentNameStr( BaseName, Component_arg ) RESULT(str) 1283!------------------------------------------------------------------------------ 1284 INTEGER, OPTIONAL, INTENT(in) :: Component_arg 1285 CHARACTER(LEN=*), INTENT(in) :: BaseName 1286!------------------------------------------------------------------------------ 1287 INTEGER :: ind, ind1, DOFsTot, DOFs, Component 1288 CHARACTER(LEN=MAX_NAME_LEN) :: str 1289!------------------------------------------------------------------------------ 1290 ind = INDEX( BaseName,'[' ) 1291 1292 Component = 0 1293 IF ( PRESENT(Component_arg) ) Component=Component_arg 1294 1295 IF ( ind<=0 ) THEN 1296 str = BaseName 1297 IF ( Component > 0 ) THEN 1298 str = TRIM(str) // ' ' // TRIM(i2s(Component) ) 1299 END IF 1300 ELSE IF( Component == 0 ) THEN 1301 str = BaseName(1:ind-1) 1302 ELSE 1303 DOFsTot = 0 1304 DO WHILE( .TRUE. ) 1305 ind1 = INDEX( BaseName(ind+1:),':' )+ind 1306 IF ( ind1 <= ind ) THEN 1307 CALL Fatal( 'ComponentName', 'Syntax error in variable definition.' ) 1308 END IF 1309 READ(BaseName(ind1+1:),'(i1)') DOFs 1310 DOFsTot = DOFsTot+DOFs 1311 IF ( DOFsTot>=Component ) EXIT 1312 ind = ind1+2 1313 END DO 1314 str = BaseName(ind+1:ind1-1) 1315 IF ( DOFs>1 ) THEN 1316 DOFs = Component - DOFsTot + DOFs 1317 str = TRIM(str) // ' ' // TRIM(i2s(DOFs) ) 1318 END IF 1319 END IF 1320!------------------------------------------------------------------------------ 1321 END FUNCTION ComponentNameStr 1322!------------------------------------------------------------------------------ 1323 1324 1325!------------------------------------------------------------------------------ 1326!> Solves a tridiagonal linear system. 1327!------------------------------------------------------------------------------ 1328 PURE SUBROUTINE SolveTriDiag( n, y, h, r ) 1329!------------------------------------------------------------------------------ 1330 INTEGER, INTENT(in) :: n 1331 REAL(KIND=dp), INTENT(out) :: r(:) 1332 REAL(KIND=dp), INTENT(in) :: y(:), h(:) 1333 1334 REAL(KIND=dp) :: s,b(n) 1335 INTEGER :: i 1336 1337 DO i=2,n-1 1338 b(i) = 2 * ( h(i-1) + h(i) ) 1339 r(i) = 3 * ( h(i) * ( y(i)-y(i-1) ) / h(i-1) + & 1340 h(i-1) * ( y(i+1)-y(i) ) / h(i) ) 1341 END DO 1342 1343 r(2) = r(2) - h(2) * r(1) 1344 DO i=2,n-2 1345 s = -h(i+1) / b(i) 1346 r(i+1) = r(i+1) + s * r(i) 1347 b(i+1) = b(i+1) + s * h(i-1) 1348 END DO 1349 1350 DO i=n-1,2,-1 1351 r(i) = (r(i) - h(i-1) * r(i+1)) / b(i) 1352 END DO 1353!------------------------------------------------------------------------------ 1354 END SUBROUTINE SolveTriDiag 1355!------------------------------------------------------------------------------ 1356 1357 1358 FUNCTION CheckMonotone(n,x) RESULT ( Monotone ) 1359 REAL(KIND=dp), INTENT(in) :: x(:) 1360 INTEGER, INTENT(in) :: n 1361 LOGICAL :: Monotone 1362 1363 INTEGER :: i 1364 1365 Monotone = .TRUE. 1366 DO i=1,n-1 1367 IF( x(i+1) <= x(i) ) THEN 1368 Monotone = .FALSE. 1369 WRITE (Message,'(E14.7,A,E14.7)') x(i),'>=',x(i+1) 1370 CALL WARN('CheckMonotone', Message) 1371 EXIT 1372 END IF 1373 END DO 1374 1375 END FUNCTION CheckMonotone 1376 1377!------------------------------------------------------------------------------ 1378!> Solver for the coefficients of a cubic spline. 1379!------------------------------------------------------------------------------ 1380 PURE SUBROUTINE CubicSpline( n,x,y,r, monotone ) 1381!------------------------------------------------------------------------------ 1382 REAL(KIND=dp), INTENT(in) :: x(:),y(:) 1383 REAL(KIND=dp), INTENT(out) :: r(:) 1384 INTEGER, INTENT(in) :: n 1385 LOGICAL, OPTIONAL, INTENT(in) :: monotone 1386 1387 REAL(KIND=dp) :: t,h(n),tau, alpha, beta 1388 INTEGER :: i 1389 LOGICAL :: mono 1390 1391 DO i=1,n-1 1392 h(i) = x(i+1) - x(i) 1393 END DO 1394 1395 r(1) = (y(2) - y(1) ) / h(1) 1396 r(n) = (y(n) - y(n-1) ) / h(n-1) 1397 1398 mono = .FALSE. 1399 IF(PRESENT(monotone)) mono = Monotone 1400 1401 IF (mono) THEN 1402 DO i=1,n-1 1403 h(i) = (y(i+1) - y(i) ) / h(i) 1404 END DO 1405 1406 DO i=2,n-1 1407 r(i) = (h(i-1) + h(i))/2 1408 END DO 1409 1410 DO i=1,n-1 1411 IF(ABS(h(i))<10*AEPS) THEN 1412 r(i) = 0._dp; r(i+1) = 0._dp 1413 CYCLE 1414 END IF 1415 1416 alpha = r(i) / h(i); beta = r(i+1) / h(i) 1417 IF ( alpha < 0._dp .OR. beta < 0._dp ) THEN 1418 r(i) = 0._dp; 1419 CYCLE 1420 END IF 1421 1422 tau = SQRT(alpha**2 + beta**2) 1423 IF(tau > 3) THEN 1424 tau = 3._dp / tau 1425 r(i) = alpha*tau*h(i); r(i+1)=beta*tau*h(i) 1426 END IF 1427 END DO 1428 ELSE 1429 CALL SolveTriDiag( n,y,h,r ) 1430 END IF 1431!------------------------------------------------------------------------------ 1432 END SUBROUTINE CubicSpline 1433!------------------------------------------------------------------------------ 1434 1435!------------------------------------------------------------------------------ 1436!> Evalulate a cubic spline. 1437!------------------------------------------------------------------------------ 1438 PURE FUNCTION CubicSplineVal(x,y,r,t) RESULT(s) 1439!------------------------------------------------------------------------------ 1440 REAL(KIND=dp) :: s 1441 REAL(KIND=dp), INTENT(in) :: x(:),y(:),r(:),t 1442!------------------------------------------------------------------------------ 1443 REAL(KIND=dp) :: a,b,c,d,h,lt 1444 1445 h = x(2)-x(1) 1446 a = -2 * ( y(2) - y(1) ) + ( r(1) + r(2) ) * h 1447 b = 3 * ( y(2) - y(1) ) - ( 2*r(1) + r(2) ) * h 1448 c = r(1) * h 1449 d = y(1) 1450 1451 lt = (t - x(1)) / h 1452 s = ((a*lt + b) * lt + c) * lt + d 1453!------------------------------------------------------------------------------ 1454 END FUNCTION CubicSplineVal 1455!------------------------------------------------------------------------------ 1456 1457 1458!------------------------------------------------------------------------------ 1459!> Evalulate derivative of cubic spline. 1460!------------------------------------------------------------------------------ 1461 PURE FUNCTION CubicSplinedVal(x,y,r,t) RESULT(s) 1462!------------------------------------------------------------------------------ 1463 REAL(KIND=dp) :: s 1464 REAL(KIND=dp), INTENT(in) :: x(:),y(:),r(:),t 1465!------------------------------------------------------------------------------ 1466 REAL(KIND=dp) :: a,b,c,h,lt 1467 1468 h = x(2)-x(1) 1469 a = -2 * ( y(2) - y(1) ) + ( r(1) + r(2) ) * h 1470 b = 3 * ( y(2) - y(1) ) - ( 2*r(1) + r(2) ) * h 1471 c = r(1) * h 1472 1473 lt = (t - x(1)) / h 1474 s = ((3*a*lt + 2*b) * lt + c)/h 1475!------------------------------------------------------------------------------ 1476 END FUNCTION CubicSplinedVal 1477!------------------------------------------------------------------------------ 1478 1479 1480!------------------------------------------------------------------------------ 1481!> Search array index such that tval(i) <= t < tval(i+1) 1482!------------------------------------------------------------------------------ 1483 PURE FUNCTION SearchInterval( tval, t ) RESULT(i) 1484!------------------------------------------------------------------------------ 1485 INTEGER :: i 1486 REAL(KIND=dp), INTENT(in) :: tval(:), t 1487!------------------------------------------------------------------------------ 1488 INTEGER :: n,n0,n1 1489!------------------------------------------------------------------------------ 1490 1491 n = SIZE(tval) 1492 1493 IF (t < tval(2)) THEN 1494 i = 1 1495 ELSE IF (t>=tval(n-1)) THEN 1496 i = n-1 1497 ELSE 1498 n0 = 1 1499 n1 = n 1500 i = (n0+n1)/2 1501 DO WHILE(.TRUE.) 1502 IF ( tval(i) <= t .AND. tval(i+1)>t ) EXIT 1503 1504 IF ( tval(i) > t ) THEN 1505 n1 = i-1 1506 ELSE 1507 n0 = i+1 1508 END IF 1509 i = (n0+n1)/2 1510 END DO 1511 END IF 1512 IF(i>n-1) i=n-1 1513 1514!------------------------------------------------------------------------------ 1515 END FUNCTION SearchInterval 1516!------------------------------------------------------------------------------ 1517 1518!------------------------------------------------------------------------------ 1519!> As SearchInterval, but doesn't assume we'll find the value in the interval 1520 !------------------------------------------------------------------------------ 1521 PURE FUNCTION SearchIntPosition( tval, t ) RESULT(i) 1522 !------------------------------------------------------------------------------ 1523 INTEGER :: i 1524 INTEGER, INTENT(in) :: tval(:), t 1525 !------------------------------------------------------------------------------ 1526 INTEGER :: n,n0,n1 1527 !------------------------------------------------------------------------------ 1528 1529 n = SIZE(tval) 1530 1531 IF (t < tval(1)) THEN 1532 i = 0 1533 ELSE IF (t>=tval(n)) THEN 1534 i = n 1535 ELSE 1536 n0 = 1 1537 n1 = n 1538 i = (n0+n1)/2 1539 DO WHILE(.TRUE.) 1540 IF ( tval(i) <= t .AND. tval(i+1)>t ) EXIT 1541 1542 IF ( tval(i) > t ) THEN 1543 n1 = i-1 1544 ELSE 1545 n0 = i+1 1546 END IF 1547 i = (n0+n1)/2 1548 END DO 1549 END IF 1550 IF(i>n) i=n 1551 1552 !------------------------------------------------------------------------------ 1553 END FUNCTION SearchIntPosition 1554 !------------------------------------------------------------------------------ 1555 1556 1557!------------------------------------------------------------------------------ 1558!> Interpolate values in a curve given by linear table or splines. 1559!------------------------------------------------------------------------------ 1560 PURE FUNCTION InterpolateCurve( TValues,FValues,T, CubicCoeff) RESULT( F ) 1561!------------------------------------------------------------------------------ 1562 REAL(KIND=dp) :: F 1563 REAL(KIND=dp), INTENT(iN) :: TValues(:),FValues(:),T 1564 REAL(KIND=dp), OPTIONAL, POINTER, INTENT(in) :: CubicCoeff(:) 1565!------------------------------------------------------------------------------ 1566 INTEGER :: i,n 1567 LOGICAL :: Cubic 1568!------------------------------------------------------------------------------ 1569 1570 n = SIZE(TValues) 1571 1572 ! This is a misuse of the interpolation in case of standard dependency 1573 ! of type y=a*x. 1574 IF( n == 1 ) THEN 1575 F = FValues(1) * T 1576 RETURN 1577 END IF 1578 1579 i = SearchInterval( Tvalues, t ) 1580 1581 Cubic = PRESENT(CubicCoeff) 1582 Cubic = Cubic .AND. T>=Tvalues(1) .AND. T<=Tvalues(n) 1583 IF ( Cubic ) Cubic = Cubic.AND.ASSOCIATED(CubicCoeff) 1584 1585 IF ( Cubic ) THEN 1586 F = CubicSplineVal(Tvalues(i:i+1),FValues(i:i+1),CubicCoeff(i:i+1),T) 1587 ELSE 1588 F = (T-TValues(i)) / (TValues(i+1)-TValues(i)) 1589 F = (1-F)*FValues(i) + F*FValues(i+1) 1590 END IF 1591 END FUNCTION InterpolateCurve 1592!------------------------------------------------------------------------------ 1593 1594 1595!------------------------------------------------------------------------------ 1596!> Derivate a curve given by linear table or splines. 1597!------------------------------------------------------------------------------ 1598 PURE FUNCTION DerivateCurve( TValues,FValues,T,CubicCoeff ) RESULT( F ) 1599!------------------------------------------------------------------------------ 1600 REAL(KIND=dp) :: F 1601 REAL(KIND=dp), INTENT(in) :: TValues(:),FValues(:),T 1602 REAL(KIND=dp), OPTIONAL, POINTER, INTENT(in) :: CubicCoeff(:) 1603!------------------------------------------------------------------------------ 1604 INTEGER :: i,n 1605 LOGICAL :: Cubic 1606!------------------------------------------------------------------------------ 1607 n = SIZE(TValues) 1608 1609 i = SearchInterval( Tvalues, t ) 1610 1611 Cubic = PRESENT(CubicCoeff) 1612 Cubic = Cubic .AND. T>=Tvalues(1) .AND. T<=Tvalues(n) 1613 IF ( Cubic ) Cubic = Cubic.AND.ASSOCIATED(CubicCoeff) 1614 1615 IF (Cubic) THEN 1616 F = CubicSplinedVal(Tvalues(i:i+1),FValues(i:i+1),CubicCoeff(i:i+1),T) 1617 ELSE 1618 F = (FValues(i+1)-FValues(i)) / (TValues(i+1)-TValues(i)) 1619 END IF 1620!------------------------------------------------------------------------------ 1621 END FUNCTION DerivateCurve 1622!------------------------------------------------------------------------------ 1623 1624 1625!------------------------------------------------------------------------------ 1626!> Integrate a curve given by linear table or splines. 1627!------------------------------------------------------------------------------ 1628 PURE SUBROUTINE CumulativeIntegral(TValues,FValues,CubicCoeff,Cumulative) 1629!------------------------------------------------------------------------------ 1630 REAL(KIND=dp), INTENT(in) :: TValues(:),FValues(:) 1631 REAL(KIND=dp), INTENT(out) :: Cumulative(:) 1632 REAL(KIND=dp), OPTIONAL, POINTER, INTENT(in) :: CubicCoeff(:) 1633!------------------------------------------------------------------------------ 1634 INTEGER :: i,n 1635 LOGICAL :: Cubic 1636 REAL(KIND=dp) :: t(2), y(2), r(2), h, a, b, c, d 1637!------------------------------------------------------------------------------ 1638 n = SIZE(TValues) 1639 1640 Cubic = PRESENT(CubicCoeff) 1641 IF ( Cubic ) Cubic = Cubic.AND.ASSOCIATED(CubicCoeff) 1642 1643 ! here only complete intervals: 1644 ! ----------------------------- 1645 Cumulative(1) = 0._dp 1646 IF ( Cubic ) THEN 1647 DO i=1,n-1 1648 t(1) = Tvalues(i) 1649 t(2) = Tvalues(i+1) 1650 1651 y(1) = FValues(i) 1652 y(2) = FValues(i+1) 1653 1654 r(1) = CubicCoeff(i) 1655 r(2) = CubicCoeff(i+1) 1656 1657 h = t(2) - t(1) 1658 1659 a = (-2 * (y(2) - y(1)) + ( r(1) + r(2)) * h)/4 1660 b = ( 3 * (y(2) - y(1)) - (2*r(1) + r(2)) * h)/3 1661 c = (r(1) * h)/2 1662 d = y(1) 1663 Cumulative(i+1) = Cumulative(i) + h*(a+b+c+d) 1664 END DO 1665 ELSE 1666 DO i=1,n-1 1667 t(1) = Tvalues(i) 1668 t(2) = Tvalues(i+1) 1669 1670 y(1) = FValues(i) 1671 y(2) = FValues(i+1) 1672 1673 h = t(2) - t(1) 1674 c = (y(2)-y(1))/2 1675 d = y(1) 1676 Cumulative(i+1) = Cumulative(i) + h*(c+d) 1677 END DO 1678 END IF 1679!------------------------------------------------------------------------------ 1680 END SUBROUTINE CumulativeIntegral 1681!------------------------------------------------------------------------------ 1682 1683!------------------------------------------------------------------------------ 1684!> Integrate a curve given by linear table or splines. 1685!------------------------------------------------------------------------------ 1686 PURE FUNCTION IntegrateCurve(TValues,FValues,CubicCoeff,T0,T1,Cumulative) RESULT(sumf) 1687!------------------------------------------------------------------------------ 1688 REAL(KIND=dp) :: sumf 1689 1690 REAL(KIND=dp), INTENT(in) :: TValues(:),FValues(:) 1691 REAL(KIND=dp), OPTIONAL, INTENT(in) :: T0, T1 1692 REAL(KIND=dp), OPTIONAL, INTENT(in) :: Cumulative(:) 1693 REAL(KIND=dp), OPTIONAL, POINTER, INTENT(in) :: CubicCoeff(:) 1694!------------------------------------------------------------------------------ 1695 INTEGER :: i,n,i0,i1 1696 LOGICAL :: Cubic 1697 REAL(KIND=dp) :: t(2), y(2), r(2), h, a, b, c, d, s0, s1, tt0, tt1 1698!------------------------------------------------------------------------------ 1699 n = SIZE(TValues) 1700 1701 tt0 = TValues(1) 1702 IF(PRESENT(t0)) tt0=t0 1703 1704 tt1 = TValues(n) 1705 IF(PRESENT(t1)) tt1=t1 1706 1707 sumf = 0._dp 1708 IF(tt0>=tt1) RETURN 1709 1710 ! t0 < first, t1 <= first 1711 IF(tt1<=Tvalues(1)) THEN 1712 t(1) = Tvalues(1) 1713 t(2) = Tvalues(2) 1714 1715 y(1) = FValues(1) 1716 y(2) = FValues(2) 1717 1718 h = t(2) - t(1) 1719 s0 = (tt0 - t(1)) / h 1720 s1 = (tt1 - t(1)) / h 1721 c = (y(2) - y(1)) / 2 1722 d = y(1) 1723 sumf = sumf + h * ((c*s1 + d)*s1 - (c*s0 + d)*s0) 1724 RETURN 1725 END IF 1726 1727 ! t0 >= last, t1 > last 1728 IF(tt0>=Tvalues(n)) THEN 1729 t(1) = Tvalues(n-1) 1730 t(2) = Tvalues(n) 1731 1732 y(1) = FValues(n-1) 1733 y(2) = FValues(n) 1734 1735 h = t(2) - t(1) 1736 s0 = (tt0 - t(1)) / h 1737 s1 = (tt1 - t(1)) / h 1738 c = (y(2) - y(1)) / 2 1739 d = y(1) 1740 sumf = sumf + h * ((c*s1 + d)*s1 - (c*s0 + d)*s0) 1741 RETURN 1742 END IF 1743 1744 ! first interval outside 1745 IF(tt0<Tvalues(1)) THEN 1746 t(1) = Tvalues(1) 1747 t(2) = Tvalues(2) 1748 1749 y(1) = FValues(1) 1750 y(2) = FValues(2) 1751 1752 h = t(2) - t(1) 1753 s0 = (tt0 - t(1)) / h 1754 c = (y(2) - y(1)) / 2 1755 d = y(1) 1756 sumf = sumf - h * (c*s0 + d)*s0 1757 tt0 = Tvalues(1) 1758 END IF 1759 1760 ! last interval outside 1761 IF(tt1>Tvalues(n)) THEN 1762 t(1) = Tvalues(n-1) 1763 t(2) = Tvalues(n) 1764 1765 y(1) = FValues(n-1) 1766 y(2) = FValues(n) 1767 1768 h = t(2) - t(1) 1769 s1 = (tt1 - t(1)) / h 1770 c = (y(2) - y(1)) / 2 1771 d = y(1) 1772 sumf = sumf + h * ( (c*s1 + d)*s1 - (c+d) ) 1773 tt1 = Tvalues(n) 1774 END IF 1775 1776 IF(tt0 >= tt1) RETURN 1777 1778 Cubic = PRESENT(CubicCoeff) 1779 IF ( Cubic ) Cubic = Cubic.AND.ASSOCIATED(CubicCoeff) 1780 1781 i0 = SearchInterval( Tvalues, tt0 ) 1782 1783 ! first (possibly partial) interval: 1784 ! ------------------------------------- 1785 t(1) = Tvalues(i0) 1786 t(2) = Tvalues(i0+1) 1787 1788 h = t(2) - t(1) 1789 s0 = (tt0-t(1))/h 1790 s1 = MIN((tt1-t(1))/h,1._dp) 1791 1792 IF(s0>0 .OR. s1<1) THEN 1793 y(1) = FValues(i0) 1794 y(2) = FValues(i0+1) 1795 1796 IF(Cubic) THEN 1797 r(1) = CubicCoeff(i0) 1798 r(2) = CubicCoeff(i0+1) 1799 1800 a = (-2 * (y(2) - y(1)) + ( r(1) + r(2)) * h)/4 1801 b = ( 3 * (y(2) - y(1)) - (2*r(1) + r(2)) * h)/3 1802 c = (r(1) * h)/2 1803 d = y(1) 1804 sumf = sumf + h * ( (((a*s1 + b)*s1 + c)*s1 + d)*s1 - & 1805 (((a*s0 + b)*s0 + c)*s0 + d)*s0 ) 1806 1807 ELSE 1808 c = (y(2)-y(1))/2 1809 d = y(1) 1810 sumf = sumf + h * ( (c*s1 + d)*s1 - (c*s0 + d)*s0 ) 1811 END IF 1812 i0 = i0 + 1 1813 tt0 = Tvalues(i0) 1814 IF(tt0 >= tt1) RETURN 1815 END IF 1816 1817 i1 = SearchInterval( Tvalues, tt1 ) 1818 1819 ! last (possibly partial) interval: 1820 ! ------------------------------------ 1821 t(1) = Tvalues(i1) 1822 t(2) = Tvalues(i1+1) 1823 1824 h = t(2) - t(1) 1825 1826 s0 = MAX((tt0-t(1))/h, 0.0_dp) 1827 s1 = (tt1-t(1))/h 1828 1829 IF(s0>0 .OR. s1<1) THEN 1830 y(1) = FValues(i1) 1831 y(2) = FValues(i1+1) 1832 1833 IF(Cubic) THEN 1834 r(1) = CubicCoeff(i1) 1835 r(2) = CubicCoeff(i1+1) 1836 1837 a = (-2 * (y(2) - y(1)) + ( r(1) + r(2)) * h)/4 1838 b = ( 3 * (y(2) - y(1)) - (2*r(1) + r(2)) * h)/3 1839 c = (r(1) * h)/2 1840 d = y(1) 1841 sumf = sumf + h * ( (((a*s1 + b)*s1 + c)*s1 + d)*s1 - & 1842 (((a*s0 + b)*s0 + c)*s0 + d)*s0 ) 1843 ELSE 1844 c = (y(2)-y(1))/2 1845 d = y(1) 1846 sumf = sumf + h * ( (c*s1 + d)*s1 - (c*s0 + d)*s0 ) 1847 END IF 1848 i1 = i1 - 1 1849 tt1 = Tvalues(i1+1) 1850 IF(tt0 >= tt1) RETURN 1851 END IF 1852 1853 ! here only complete intervals: 1854 ! ----------------------------- 1855 1856 IF(PRESENT(Cumulative)) THEN 1857 sumf = sumf + Cumulative(i1+1) - Cumulative(i0) 1858 RETURN 1859 END IF 1860 1861 IF ( Cubic ) THEN 1862 DO i=i0,i1 1863 t(1) = Tvalues(i) 1864 t(2) = Tvalues(i+1) 1865 1866 y(1) = FValues(i) 1867 y(2) = FValues(i+1) 1868 1869 r(1) = CubicCoeff(i) 1870 r(2) = CubicCoeff(i+1) 1871 1872 h = t(2) - t(1) 1873 1874 a = (-2 * (y(2) - y(1)) + ( r(1) + r(2)) * h)/4 1875 b = ( 3 * (y(2) - y(1)) - (2*r(1) + r(2)) * h)/3 1876 c = (r(1) * h)/2 1877 d = y(1) 1878 sumf = sumf + h * (a+b+c+d) 1879 END DO 1880 ELSE 1881 DO i=i0,i1 1882 t(1) = Tvalues(i) 1883 t(2) = Tvalues(i+1) 1884 1885 y(1) = FValues(i) 1886 y(2) = FValues(i+1) 1887 1888 h = t(2) - t(1) 1889 c = (y(2)-y(1))/2 1890 d = y(1) 1891 sumf = sumf + h * (c+d) 1892 END DO 1893 END IF 1894!------------------------------------------------------------------------------ 1895 END FUNCTION IntegrateCurve 1896!------------------------------------------------------------------------------ 1897 1898 1899!------------------------------------------------------------------------------ 1900!> Solves a 2 x 2 linear system. 1901!------------------------------------------------------------------------------ 1902 SUBROUTINE SolveLinSys2x2( A, x, b ) 1903!------------------------------------------------------------------------------ 1904 REAL(KIND=dp), INTENT(out) :: x(:) 1905 REAL(KIND=dp), INTENT(in) :: A(:,:),b(:) 1906!------------------------------------------------------------------------------ 1907 REAL(KIND=dp) :: detA 1908!------------------------------------------------------------------------------ 1909 detA = A(1,1) * A(2,2) - A(1,2) * A(2,1) 1910 1911 IF ( detA == 0.0d0 ) THEN 1912 WRITE( Message, * ) 'Singular matrix, sorry!' 1913 CALL Error( 'SolveLinSys2x2', Message ) 1914 RETURN 1915 END IF 1916 1917 detA = 1.0d0 / detA 1918 x(1) = detA * (A(2,2) * b(1) - A(1,2) * b(2)) 1919 x(2) = detA * (A(1,1) * b(2) - A(2,1) * b(1)) 1920!------------------------------------------------------------------------------ 1921 END SUBROUTINE SolveLinSys2x2 1922!------------------------------------------------------------------------------ 1923 1924 1925!------------------------------------------------------------------------------ 1926!> Solves a 3 x 3 linear system. 1927!------------------------------------------------------------------------------ 1928 SUBROUTINE SolveLinSys3x3( A, x, b ) 1929!------------------------------------------------------------------------------ 1930 REAL(KIND=dp), INTENT(out) :: x(:) 1931 REAL(KIND=dp), INTENT(in) :: A(:,:),b(:) 1932!------------------------------------------------------------------------------ 1933 REAL(KIND=dp) :: C(2,2),y(2),g(2),s,t,q 1934!------------------------------------------------------------------------------ 1935 1936 IF ( ABS(A(1,1))>ABS(A(1,2)) .AND. ABS(A(1,1))>ABS(A(1,3)) ) THEN 1937 q = 1.0d0 / A(1,1) 1938 s = q * A(2,1) 1939 t = q * A(3,1) 1940 C(1,1) = A(2,2) - s * A(1,2) 1941 C(1,2) = A(2,3) - s * A(1,3) 1942 C(2,1) = A(3,2) - t * A(1,2) 1943 C(2,2) = A(3,3) - t * A(1,3) 1944 1945 g(1) = b(2) - s * b(1) 1946 g(2) = b(3) - t * b(1) 1947 CALL SolveLinSys2x2( C,y,g ) 1948 1949 x(2) = y(1) 1950 x(3) = y(2) 1951 x(1) = q * ( b(1) - A(1,2) * x(2) - A(1,3) * x(3) ) 1952 ELSE IF ( ABS(A(1,2)) > ABS(A(1,3)) ) THEN 1953 q = 1.0d0 / A(1,2) 1954 s = q * A(2,2) 1955 t = q * A(3,2) 1956 C(1,1) = A(2,1) - s * A(1,1) 1957 C(1,2) = A(2,3) - s * A(1,3) 1958 C(2,1) = A(3,1) - t * A(1,1) 1959 C(2,2) = A(3,3) - t * A(1,3) 1960 1961 g(1) = b(2) - s * b(1) 1962 g(2) = b(3) - t * b(1) 1963 CALL SolveLinSys2x2( C,y,g ) 1964 1965 x(1) = y(1) 1966 x(3) = y(2) 1967 x(2) = q * ( b(1) - A(1,1) * x(1) - A(1,3) * x(3) ) 1968 ELSE 1969 q = 1.0d0 / A(1,3) 1970 s = q * A(2,3) 1971 t = q * A(3,3) 1972 C(1,1) = A(2,1) - s * A(1,1) 1973 C(1,2) = A(2,2) - s * A(1,2) 1974 C(2,1) = A(3,1) - t * A(1,1) 1975 C(2,2) = A(3,2) - t * A(1,2) 1976 1977 g(1) = b(2) - s * b(1) 1978 g(2) = b(3) - t * b(1) 1979 CALL SolveLinSys2x2( C,y,g ) 1980 1981 x(1) = y(1) 1982 x(2) = y(2) 1983 x(3) = q * ( b(1) - A(1,1) * x(1) - A(1,2) * x(2) ) 1984 END IF 1985!------------------------------------------------------------------------------ 1986 END SUBROUTINE SolveLinSys3x3 1987!------------------------------------------------------------------------------ 1988 1989 1990 1991!------------------------------------------------------------------------------ 1992 SUBROUTINE ClearMatrix( Matrix ) 1993 TYPE(Matrix_t), POINTER, INTENT(in) :: Matrix 1994INCLUDE "mpif.h" 1995 1996 Matrix % FORMAT = MATRIX_CRS 1997 1998 NULLIFY( Matrix % Child ) 1999 NULLIFY( Matrix % Parent ) 2000 NULLIFY( Matrix % EMatrix ) 2001 NULLIFY( Matrix % ConstraintMatrix ) 2002 2003 NULLIFY( Matrix % Perm ) 2004 NULLIFY( Matrix % InvPerm ) 2005 2006 NULLIFY( Matrix % Cols ) 2007 NULLIFY( Matrix % Rows ) 2008 NULLIFY( Matrix % Diag ) 2009 2010 NULLIFY( Matrix % RHS ) 2011 NULLIFY( Matrix % Force ) 2012 NULLIFY( Matrix % RHS_im ) 2013 2014 NULLIFY( Matrix % Values ) 2015 NULLIFY( Matrix % ILUValues ) 2016 NULLIFY( Matrix % MassValues ) 2017 NULLIFY( Matrix % DampValues ) 2018 2019 NULLIFY( Matrix % BulkRHS ) 2020 NULLIFY( Matrix % BulkValues ) 2021 2022 NULLIFY( Matrix % BulkResidual ) 2023 2024 NULLIFY( Matrix % ILUCols ) 2025 NULLIFY( Matrix % ILURows ) 2026 NULLIFY( Matrix % ILUDiag ) 2027 2028 NULLIFY( Matrix % CRHS ) 2029 NULLIFY( Matrix % CForce ) 2030 2031 NULLIFY( Matrix % ParMatrix ) 2032 2033 NULLIFY( Matrix % CValues ) 2034 NULLIFY( Matrix % CILUValues ) 2035 NULLIFY( Matrix % CMassValues ) 2036 NULLIFY( Matrix % CDampValues ) 2037 2038! NULLIFY( Matrix % GRows ) 2039! NULLIFY( Matrix % RowOwner ) 2040 NULLIFY( Matrix % GOrder ) 2041 NULLIFY( Matrix % EPerm ) 2042 2043 NULLIFY( Matrix % ParMatrix ) 2044 2045 NULLIFY( Matrix % ParallelInfo ) 2046#ifdef HAVE_UMFPACK 2047 Matrix % UMFPack_Numeric = 0 2048#endif 2049 2050 Matrix % Cholesky = .FALSE. 2051 Matrix % Lumped = .FALSE. 2052 Matrix % Ordered = .FALSE. 2053 Matrix % COMPLEX = .FALSE. 2054 Matrix % Symmetric = .FALSE. 2055 Matrix % SolveCount = 0 2056 Matrix % NumberOfRows = 0 2057 2058 Matrix % ProjectorBC = 0 2059 Matrix % ProjectorType = PROJECTOR_TYPE_DEFAULT 2060 2061 Matrix % Solver => NULL() 2062 2063 Matrix % DGMatrix = .FALSE. 2064 Matrix % Comm = ELMER_COMM_WORLD 2065 2066 END SUBROUTINE ClearMatrix 2067 2068 2069!------------------------------------------------------------------------------ 2070 FUNCTION AllocateMatrix() RESULT(Matrix) 2071!------------------------------------------------------------------------------ 2072 TYPE(Matrix_t), POINTER :: Matrix 2073!------------------------------------------------------------------------------ 2074 ALLOCATE( Matrix ) 2075 CALL ClearMatrix( Matrix ) 2076!------------------------------------------------------------------------------ 2077 END FUNCTION AllocateMatrix 2078!------------------------------------------------------------------------------ 2079 2080 2081!------------------------------------------------------------------------------ 2082 RECURSIVE SUBROUTINE FreeQuadrantTree( Root ) 2083!------------------------------------------------------------------------------ 2084 TYPE(Quadrant_t), POINTER :: Root 2085 2086 INTEGER :: i 2087 2088 IF ( .NOT. ASSOCIATED( Root ) ) RETURN 2089 2090 IF ( ASSOCIATED(Root % Elements) ) DEALLOCATE( Root % Elements ) 2091 2092 IF ( ASSOCIATED( Root % ChildQuadrants ) ) THEN 2093 DO i=1,SIZE(Root % ChildQuadrants) 2094 CALL FreeQuadrantTree( Root % ChildQuadrants(i) % Quadrant ) 2095 END DO 2096 DEALLOCATE( Root % ChildQuadrants ) 2097 NULLIFY( Root % ChildQuadrants ) 2098 END IF 2099 2100 DEALLOCATE( Root ) 2101!------------------------------------------------------------------------------ 2102 END SUBROUTINE FreeQuadrantTree 2103!------------------------------------------------------------------------------ 2104 2105 2106!------------------------------------------------------------------------------ 2107 SUBROUTINE AllocateRealVector( F, n, From, FailureMessage ) 2108!------------------------------------------------------------------------------ 2109 REAL(KIND=dp), POINTER :: F(:) 2110 INTEGER :: n 2111 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2112!------------------------------------------------------------------------------ 2113 INTEGER :: istat 2114!------------------------------------------------------------------------------ 2115 2116 istat = -1 2117 ALLOCATE( F(n), STAT=istat ) 2118 IF ( istat /= 0 ) THEN 2119 IF ( PRESENT( FailureMessage ) ) THEN 2120 WRITE( Message, * )'Unable to allocate ', n, ' element real array.' 2121 CALL Error( 'AllocateRealVector', Message ) 2122 IF ( PRESENT( From ) ) THEN 2123 WRITE( Message, * )'Requested From: ', TRIM(From) 2124 CALL Error( 'AllocateRealVector', Message ) 2125 END IF 2126 IF ( PRESENT( FailureMessage ) ) THEN 2127 CALL Fatal( 'AllocateRealVector', FailureMessage ) 2128 END IF 2129 END IF 2130 END IF 2131!------------------------------------------------------------------------------ 2132 END SUBROUTINE AllocateRealVector 2133!------------------------------------------------------------------------------ 2134 2135 2136!------------------------------------------------------------------------------ 2137 SUBROUTINE AllocateComplexVector( f, n, From, FailureMessage ) 2138!------------------------------------------------------------------------------ 2139 COMPLEX(KIND=dp), POINTER :: f(:) 2140 INTEGER :: n 2141 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2142!------------------------------------------------------------------------------ 2143 INTEGER :: istat 2144!------------------------------------------------------------------------------ 2145 2146 istat = -1 2147 ALLOCATE( f(n), STAT=istat ) 2148 IF ( istat /= 0 ) THEN 2149 IF ( PRESENT( FailureMessage ) ) THEN 2150 WRITE( Message, * )'Unable to allocate ', n, ' element real array.' 2151 CALL Error( 'AllocateComplexVector', Message ) 2152 IF ( PRESENT( From ) ) THEN 2153 WRITE( Message, * )'Requested From: ', TRIM(From) 2154 CALL Error( 'AllocateComplexVector', Message ) 2155 END IF 2156 IF ( PRESENT( FailureMessage ) ) THEN 2157 CALL Fatal( 'AllocateComplexVector', FailureMessage ) 2158 END IF 2159 END IF 2160 END IF 2161!------------------------------------------------------------------------------ 2162 END SUBROUTINE AllocateComplexVector 2163!------------------------------------------------------------------------------ 2164 2165 2166!------------------------------------------------------------------------------ 2167 SUBROUTINE AllocateIntegerVector( f, n, From, FailureMessage ) 2168!------------------------------------------------------------------------------ 2169 INTEGER, POINTER :: f(:) 2170 INTEGER :: n 2171 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2172!------------------------------------------------------------------------------ 2173 INTEGER :: istat 2174!------------------------------------------------------------------------------ 2175 2176 istat = -1 2177 ALLOCATE( f(n), STAT=istat ) 2178 IF ( istat /= 0 ) THEN 2179 IF ( PRESENT( FailureMessage ) ) THEN 2180 WRITE( Message, * )'Unable to allocate ', n, ' element integer array.' 2181 CALL Error( 'AllocateIntegerVector', Message ) 2182 IF ( PRESENT( From ) ) THEN 2183 WRITE( Message, * )'Requested From: ', TRIM(From) 2184 CALL Error( 'AllocateIntegerVector', Message ) 2185 END IF 2186 IF ( PRESENT( FailureMessage ) ) THEN 2187 CALL Fatal( 'AllocateIntegerVector', FailureMessage ) 2188 END IF 2189 END IF 2190 END IF 2191!------------------------------------------------------------------------------ 2192 END SUBROUTINE AllocateIntegerVector 2193!------------------------------------------------------------------------------ 2194 2195 2196!------------------------------------------------------------------------------ 2197 SUBROUTINE AllocateLogicalVector( f, n, From, FailureMessage ) 2198!------------------------------------------------------------------------------ 2199 LOGICAL, POINTER :: f(:) 2200 INTEGER :: n 2201 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2202!------------------------------------------------------------------------------ 2203 INTEGER :: istat 2204!------------------------------------------------------------------------------ 2205 2206 istat = -1 2207 ALLOCATE( f(n), STAT=istat ) 2208 IF ( istat /= 0 ) THEN 2209 IF ( PRESENT( FailureMessage ) ) THEN 2210 WRITE( Message, * )'Unable to allocate ', n, ' element integer array.' 2211 CALL Error( 'AllocateLogicalVector', Message ) 2212 IF ( PRESENT( From ) ) THEN 2213 WRITE( Message, * )'Requested From: ', TRIM(From) 2214 CALL Error( 'AllocateLogicalVector', Message ) 2215 END IF 2216 IF ( PRESENT( FailureMessage ) ) THEN 2217 CALL Fatal( 'AllocateLogicalVector', FailureMessage ) 2218 END IF 2219 END IF 2220 END IF 2221!------------------------------------------------------------------------------ 2222 END SUBROUTINE AllocateLogicalVector 2223!------------------------------------------------------------------------------ 2224 2225 2226!------------------------------------------------------------------------------ 2227 SUBROUTINE AllocateElementVector( f, n, From, FailureMessage ) 2228!------------------------------------------------------------------------------ 2229 TYPE(Element_t), POINTER :: f(:) 2230 INTEGER :: n 2231 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2232!------------------------------------------------------------------------------ 2233 INTEGER :: istat 2234!------------------------------------------------------------------------------ 2235 2236 istat = -1 2237 ALLOCATE( f(n), STAT=istat ) 2238 IF ( istat /= 0 ) THEN 2239 WRITE( Message, * )'Unable to allocate ', n, ' element integer array.' 2240 CALL Error( 'AllocateElementVector', Message ) 2241 IF ( PRESENT( From ) ) THEN 2242 WRITE( Message, * )'Requested From: ', TRIM(From) 2243 CALL Error( 'AllocateElementVector', Message ) 2244 END IF 2245 IF ( PRESENT( FailureMessage ) ) THEN 2246 CALL Fatal( 'AllocateElementVector', FailureMessage ) 2247 END IF 2248 END IF 2249!------------------------------------------------------------------------------ 2250 END SUBROUTINE AllocateElementVector 2251!------------------------------------------------------------------------------ 2252 2253 2254!------------------------------------------------------------------------------ 2255 SUBROUTINE AllocateRealArray( f, n1, n2, From, FailureMessage ) 2256!------------------------------------------------------------------------------ 2257 REAL(KIND=dp), POINTER :: f(:,:) 2258 INTEGER :: n1,n2 2259 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2260!------------------------------------------------------------------------------ 2261 INTEGER :: istat 2262!------------------------------------------------------------------------------ 2263 2264 istat = -1 2265 ALLOCATE( f(n1,n2), STAT=istat ) 2266 if ( istat /= 0 ) THEN 2267 WRITE( Message, * )'Unable to allocate ', n1, ' by ', n2, ' element real matrix.' 2268 CALL Error( 'AllocateRealArray', Message ) 2269 IF ( PRESENT( From ) ) THEN 2270 WRITE( Message, * )'Requested From: ', TRIM(From) 2271 CALL Error( 'AllocateRealArray', Message ) 2272 END IF 2273 IF ( PRESENT( FailureMessage ) ) THEN 2274 CALL Fatal( 'AllocateRealArray', FailureMessage ) 2275 END IF 2276 END IF 2277!------------------------------------------------------------------------------ 2278 END SUBROUTINE AllocateRealArray 2279!------------------------------------------------------------------------------ 2280 2281!------------------------------------------------------------------------------ 2282 SUBROUTINE AllocateComplexArray( f, n1, n2, From, FailureMessage ) 2283!------------------------------------------------------------------------------ 2284 COMPLEX(KIND=dp), POINTER :: f(:,:) 2285 INTEGER :: n1,n2 2286 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2287!------------------------------------------------------------------------------ 2288 INTEGER :: istat 2289!------------------------------------------------------------------------------ 2290 2291 istat = -1 2292 ALLOCATE( f(n1,n2), STAT=istat ) 2293 if ( istat /= 0 ) THEN 2294 WRITE( Message, * )'Unable to allocate ', n1, ' by ', n2, ' element real matrix.' 2295 CALL Error( 'AllocateComplexArray', Message ) 2296 IF ( PRESENT( From ) ) THEN 2297 WRITE( Message, * )'Requested From: ', TRIM(From) 2298 CALL Error( 'AllocateComplexArray', Message ) 2299 END IF 2300 IF ( PRESENT( FailureMessage ) ) THEN 2301 CALL Fatal( 'AllocateComplexArray', FailureMessage ) 2302 END IF 2303 END IF 2304!------------------------------------------------------------------------------ 2305 END SUBROUTINE AllocateComplexArray 2306!------------------------------------------------------------------------------ 2307 2308 2309!------------------------------------------------------------------------------ 2310 SUBROUTINE AllocateIntegerArray( f, n1, n2, From, FailureMessage ) 2311!------------------------------------------------------------------------------ 2312 INTEGER, POINTER :: f(:,:) 2313 INTEGER :: n1,n2 2314 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2315!------------------------------------------------------------------------------ 2316 INTEGER :: istat 2317!------------------------------------------------------------------------------ 2318 2319 istat = -1 2320 IF ( n1 > 0 .AND. n2 > 0 ) THEN 2321 ALLOCATE( f(n1,n2), STAT=istat ) 2322 END IF 2323 IF ( istat /= 0 ) THEN 2324 WRITE( Message, * )'Unable to allocate ', n1, ' by ', n2, ' element integer matrix.' 2325 CALL Error( 'AllocateIntegerArray', Message ) 2326 IF ( PRESENT( From ) ) THEN 2327 WRITE( Message, * )'Requested From: ', TRIM(From) 2328 CALL Error( 'AllocateIntegerArray', Message ) 2329 END IF 2330 IF ( PRESENT( FailureMessage ) ) THEN 2331 CALL Fatal( 'AllocateIntegerArray', FailureMessage ) 2332 END IF 2333 END IF 2334!------------------------------------------------------------------------------ 2335 END SUBROUTINE AllocateIntegerArray 2336!------------------------------------------------------------------------------ 2337 2338 2339 2340!------------------------------------------------------------------------------ 2341 SUBROUTINE AllocateLogicalArray( f, n1, n2, From, FailureMessage ) 2342!------------------------------------------------------------------------------ 2343 LOGICAL, POINTER :: f(:,:) 2344 INTEGER :: n1,n2 2345 CHARACTER(LEN=*), OPTIONAL :: From, FailureMessage 2346!------------------------------------------------------------------------------ 2347 INTEGER :: istat 2348!------------------------------------------------------------------------------ 2349 2350 istat = -1 2351 IF ( n1 > 0 .AND. n2 > 0 ) THEN 2352 ALLOCATE( f(n1,n2), STAT=istat ) 2353 END IF 2354 IF ( istat /= 0 ) THEN 2355 WRITE( Message, * )'Unable to allocate ', n1, ' by ', n2, ' element integer matrix.' 2356 CALL Error( 'AllocateLogicalArray', Message ) 2357 IF ( PRESENT( From ) ) THEN 2358 WRITE( Message, * )'Requested From: ', TRIM(From) 2359 CALL Error( 'AllocateLogicalArray', Message ) 2360 END IF 2361 IF ( PRESENT( FailureMessage ) ) THEN 2362 CALL Fatal( 'AllocateLogicalArray', FailureMessage ) 2363 END IF 2364 END IF 2365!------------------------------------------------------------------------------ 2366 END SUBROUTINE AllocateLogicalArray 2367!------------------------------------------------------------------------------ 2368 2369 ! Pad given integer value to be the next largest multiple of nbyte 2370 FUNCTION IntegerNBytePad(val, nbyte) RESULT(padval) 2371 IMPLICIT NONE 2372 2373 INTEGER, INTENT(IN) :: val, nbyte 2374 INTEGER :: padval 2375 ! Parameters and variables 2376 INTEGER, PARAMETER :: bytesinint = KIND(val) 2377 INTEGER :: nbytesinint 2378 2379 ! Compute number of nbytes in int 2380 nbytesinint = nbyte/bytesinint 2381 ! Compute value padded to multiples of n-byte 2382 padval=((val-1)/nbytesinint)*nbytesinint+nbytesinint 2383 END FUNCTION IntegerNBytePad 2384 2385 ! Pad given value to be the next largest multiple of nbyte 2386 FUNCTION NBytePad(val, bytesinelem, nbyte) RESULT(padval) 2387 IMPLICIT NONE 2388 2389 INTEGER, INTENT(IN) :: val, bytesinelem, nbyte 2390 INTEGER :: padval 2391 ! Variables 2392 INTEGER :: nbytesinelem 2393 2394 ! Compute number of nbytes in a single element 2395 nbytesinelem = nbyte/bytesinelem 2396 ! Compute value padded to multiples of n-byte 2397 padval=((val-1)/nbytesinelem)*nbytesinelem+nbytesinelem 2398 END FUNCTION NBytePad 2399 2400!------------------------------------------------------------------------------ 2401!> Given the filename0 (and suffix0) find the 1st free filename 2402!> that does not exist in the current working directory 2403!------------------------------------------------------------------------------ 2404 2405 FUNCTION NextFreeFilename(Filename0,Suffix0,LastExisting) RESULT (Filename) 2406 2407 CHARACTER(LEN=MAX_NAME_LEN) :: Filename0 2408 CHARACTER(LEN=MAX_NAME_LEN), OPTIONAL :: Suffix0 2409 LOGICAL, OPTIONAL :: LastExisting 2410 CHARACTER(LEN=MAX_NAME_LEN) :: Filename 2411 CHARACTER(LEN=MAX_NAME_LEN) :: Prefix, Suffix, PrevFilename 2412 LOGICAL :: FileIs 2413 INTEGER :: No, ind, len 2414 2415 ind = INDEX( FileName0,'.',.TRUE. ) 2416 len = LEN_TRIM(Filename0) 2417 IF(ind > 0) THEN 2418 Prefix = Filename0(1:ind-1) 2419 Suffix = Filename0(ind:len) 2420 ELSE 2421 Prefix = Filename0(1:len) 2422 IF(PRESENT(Suffix0)) THEN 2423 Suffix = '.'//TRIM(Suffix0) 2424 ELSE 2425 Suffix = '.dat' 2426 END IF 2427 END IF 2428 2429 DO No = 1,9999 2430 IF( No > 0 ) PrevFilename = Filename 2431 IF( No < 10) THEN 2432 WRITE( FileName,'(A,I1,A)') TRIM(Prefix),No,TRIM(Suffix) 2433 ELSE IF( No < 100) THEN 2434 WRITE( FileName,'(A,I2,A)') TRIM(Prefix),No,TRIM(Suffix) 2435 ELSE IF( No < 1000) THEN 2436 WRITE( FileName,'(A,I3,A)') TRIM(Prefix),No,TRIM(Suffix) 2437 ELSE IF( No < 10000) THEN 2438 WRITE( FileName,'(A,I4,A)') TRIM(Prefix),No,TRIM(Suffix) 2439 END IF 2440 INQUIRE( FILE=Filename, EXIST=FileIs ) 2441 IF(.NOT. FileIs) EXIT 2442 END DO 2443 2444 IF( PRESENT(LastExisting)) THEN 2445 IF( LastExisting ) Filename = PrevFilename 2446 END IF 2447 2448!------------------------------------------------------------------------------ 2449 END FUNCTION NextFreeFilename 2450!------------------------------------------------------------------------------ 2451 2452!------------------------------------------------------------------------------ 2453!> Given the filename0 add a string related to the partitioning. 2454!------------------------------------------------------------------------------ 2455 2456 FUNCTION AddFilenameParSuffix(Filename0,Suffix0,Parallel,MyPe) RESULT (Filename) 2457 2458 CHARACTER(LEN=MAX_NAME_LEN) :: Filename0 2459 CHARACTER(LEN=*), OPTIONAL :: Suffix0 2460 LOGICAL :: Parallel 2461 INTEGER :: MyPe 2462 CHARACTER(LEN=MAX_NAME_LEN) :: Filename 2463 CHARACTER(LEN=MAX_NAME_LEN) :: Prefix, Suffix 2464 INTEGER :: No, ind, len 2465 2466 2467 ind = INDEX( FileName0,'.',.TRUE. ) 2468 len = LEN_TRIM(Filename0) 2469 2470 ! If the only dot is the first one it only related to the current working directory. 2471 IF(ind > 1) THEN 2472 Prefix = Filename0(1:ind-1) 2473 Suffix = Filename0(ind:len) 2474 ELSE 2475 Prefix = Filename0(1:len) 2476 IF(PRESENT(Suffix0)) THEN 2477 Suffix = '.'//TRIM(Suffix0) 2478 ELSE 2479 Suffix = '.dat' 2480 END IF 2481 END IF 2482 2483 IF( Parallel ) THEN 2484 No = MyPe + 1 2485 IF( No < 10000 ) THEN 2486 WRITE( FileName,'(A,I4.4,A)') TRIM(Prefix),No,TRIM(Suffix) 2487 ELSE 2488 WRITE( FileName,'(A,I0,A)') TRIM(Prefix),No,TRIM(Suffix) 2489 END IF 2490 ELSE 2491 FileName = TRIM(Prefix)//TRIM(Suffix) 2492 END IF 2493 2494!------------------------------------------------------------------------------ 2495 END FUNCTION AddFilenameParSuffix 2496!------------------------------------------------------------------------------ 2497 2498 2499 !---------------------------------------------------------< 2500 !> Returns values from a normal distribution to be used in 2501 !> thermal velocity distribution, for example. 2502 !--------------------------------------------------------- 2503 FUNCTION NormalRandom() RESULT ( normalrand ) 2504 2505 REAL(KIND=dp) :: normalrand,mean 2506 INTEGER :: flag = 0 2507 REAL(KIND=dp) :: fac,gsave,rsq,r1,r2 2508 2509 SAVE flag,gsave 2510 2511 IF (flag == 0) THEN 2512 rsq=2.0_dp 2513 2514 DO WHILE(rsq >= 1.0_dp .OR. rsq == 0.0_dp ) 2515 CALL RANDOM_NUMBER(r1) 2516 CALL RANDOM_NUMBER(r2) 2517 r1 = 2.0_dp * r1 - 1.0_dp 2518 r2 = 2.0_dp * r2 - 1.0_dp 2519 rsq = r1*r1 + r2*r2 2520 ENDDO 2521 2522 fac = SQRT(-2.0_dp * LOG(rsq) / rsq) 2523 gsave = r1 * fac 2524 normalrand = r2 * fac 2525 flag = 1 2526 ELSE 2527 normalrand = gsave 2528 flag = 0 2529 ENDIF 2530 2531 END FUNCTION NormalRandom 2532 2533 2534 !--------------------------------------------------------- 2535 !> Returns values from a even distribution [0,1] 2536 !--------------------------------------------------------- 2537 FUNCTION EvenRandom() RESULT ( rand ) 2538 REAL(KIND=dp) :: rand 2539 CALL RANDOM_NUMBER(rand) 2540 END FUNCTION EvenRandom 2541 2542 2543 SUBROUTINE ForceLoad 2544 CALL MPI_SEND() 2545 END SUBROUTINE ForceLoad 2546 2547 2548END MODULE GeneralUtils 2549 2550 2551!--------------------------------------------------------- 2552!> Module mainly for writing xml based vtk files. 2553!> The idea is that same routines save both the ascii 2554!> and binary format. 2555!--------------------------------------------------------- 2556MODULE AscBinOutputUtils 2557 2558 2559 USE Types 2560 IMPLICIT NONE 2561 2562 LOGICAL, PRIVATE :: AsciiOutput, SinglePrec, CalcSum = .FALSE. 2563 INTEGER, PRIVATE :: VtuUnit = 0, BufferSize = 0 2564 REAL, POINTER, PRIVATE :: FVals(:) 2565 REAL(KIND=dp), POINTER, PRIVATE :: DVals(:) 2566 INTEGER, POINTER, PRIVATE :: IVals(:) 2567 INTEGER, PRIVATE :: INoVals, NoVals 2568 REAL(KIND=dp) :: RSum = 0.0_dp, Isum = 0.0_dp 2569 INTEGER :: Rcount = 0, Icount = 0, Scount = 0, Ssum = 0 2570 2571 SAVE :: AsciiOutput, SinglePrec, VtuUnit, BufferSize, & 2572 FVals, DVals, IVals, CalcSum, Rsum, Isum, Ssum, RCount, Icount, Scount 2573 2574 2575 2576CONTAINS 2577 2578 2579 ! Initialize the buffer for writing, choose mode etc. 2580 !----------------------------------------------------------------- 2581 SUBROUTINE AscBinWriteInit( IsAscii, IsSingle, UnitNo, BufSize ) 2582 2583 LOGICAL :: IsAscii, IsSingle 2584 INTEGER :: UnitNo, BufSize 2585 2586 AsciiOutput = IsAscii 2587 SinglePrec = IsSingle 2588 VtuUnit = UnitNo 2589 BufferSize = BufSize 2590 2591 CALL Info('AscBinWriteInit','Initializing buffered ascii/binary writing',Level=8) 2592 IF( AsciiOutput ) THEN 2593 CALL Info('AscBinWriteInit','Writing in ascii',Level=10) 2594 ELSE 2595 CALL Info('AscBinWriteInit','Writing in binary',Level=10) 2596 END IF 2597 2598 IF( SinglePrec ) THEN 2599 CALL Info('AscBinWriteInit','Writing in single precision',Level=10) 2600 ELSE 2601 CALL Info('AscBinWriteInit','Writing in double precision',Level=10) 2602 END IF 2603 2604 WRITE(Message,'(A,I0)') 'Writing to unit number: ',VtuUnit 2605 CALL Info('AscBinWriteInit',Message,Level=10) 2606 2607 IF(.NOT. AsciiOutput ) THEN 2608 WRITE(Message,'(A,I0)') 'Size of buffer is: ',BufferSize 2609 CALL Info('AscBinWriteInit',Message,Level=10) 2610 2611 ALLOCATE( Ivals( BufferSize ) ) 2612 IF( SinglePrec ) THEN 2613 ALLOCATE( FVals( BufferSize ) ) 2614 ELSE 2615 ALLOCATE( Dvals( BufferSize ) ) 2616 END IF 2617 2618 INoVals = 0 2619 NoVals = 0 2620 END IF 2621 2622 END SUBROUTINE AscBinWriteInit 2623 2624 2625 ! Free the buffer, next buffer can be different in size and type 2626 !--------------------------------------------------------------- 2627 SUBROUTINE AscBinWriteFree() 2628 2629 CALL Info('AscBinWriteFree','Terminating buffered ascii/binary writing',Level=10) 2630 2631 IF( AsciiOutput ) RETURN 2632 2633 IF( SinglePrec ) THEN 2634 DEALLOCATE( FVals ) 2635 ELSE 2636 DEALLOCATE( DVals ) 2637 END IF 2638 DEALLOCATE( IVals ) 2639 2640 BufferSize = 0 2641 VtuUnit = 0 2642 2643 END SUBROUTINE AscBinWriteFree 2644 2645 2646 2647 ! The writing of xml strings is done here to allow easier modification 2648 ! of output strategies. 2649 !------------------------------------------------------------------------- 2650 SUBROUTINE AscBinStrWrite( Str ) 2651 2652 CHARACTER(LEN=1024) :: Str 2653 INTEGER, PARAMETER :: VtuUnit = 58 2654 2655 WRITE( VtuUnit ) TRIM(Str) 2656 IF( CalcSum ) THEN 2657 Scount = Scount + 1 2658 Ssum = Ssum + len_trim( Str ) 2659 END IF 2660 2661 2662 END SUBROUTINE AscBinStrWrite 2663 2664 2665 ! Write a binary value, either in single or double precision 2666 !------------------------------------------------------------------------ 2667 SUBROUTINE AscBinRealWrite( val, EmptyBuffer ) 2668 2669 INTEGER, PARAMETER :: VtuUnit = 58 2670 REAL(KIND=dp) :: val 2671 LOGICAL, OPTIONAL :: EmptyBuffer 2672 LOGICAL :: Empty 2673 CHARACTER(LEN=1024) :: Str 2674 2675 IF( VtuUnit == 0 ) THEN 2676 CALL Fatal('AscBinRealWrite','Buffer not initialized for writing') 2677 END IF 2678 2679 IF( PRESENT( EmptyBuffer ) ) THEN 2680 Empty = EmptyBuffer 2681 ELSE 2682 Empty = .FALSE. 2683 END IF 2684 2685 ! Ascii output is not buffered 2686 IF( AsciiOutput ) THEN 2687 IF( Empty ) RETURN 2688 IF( ABS( val ) <= TINY ( val ) ) THEN 2689 WRITE(Str,'(A)') " 0.0" 2690 ELSE IF( SinglePrec ) THEN 2691 WRITE( Str,'(ES12.3E3)') val 2692 ELSE 2693 WRITE( Str,'(ES16.7E3)') val 2694 END IF 2695 WRITE( VtuUnit ) TRIM(Str) 2696 IF( CalcSum ) THEN 2697 Rcount = Rcount + 1 2698 Rsum = Rsum + val 2699 END IF 2700 RETURN 2701 END IF 2702 2703 ! Buffered binary output 2704 IF( Empty .OR. NoVals == BufferSize ) THEN 2705 IF( NoVals == 0 ) THEN 2706 RETURN 2707 ELSE IF( SinglePrec ) THEN 2708 WRITE( VtuUnit ) Fvals(1:NoVals) 2709 ELSE 2710 WRITE( VtuUnit ) DVals(1:NoVals) 2711 END IF 2712 2713 IF( CalcSum ) THEN 2714 Rcount = Rcount + NoVals 2715 IF( SinglePrec ) THEN 2716 Rsum = 1.0_dp * SUM( Fvals(1:NoVals) ) 2717 ELSE 2718 Rsum = SUM( DVals(1:NoVals) ) 2719 END IF 2720 END IF 2721 2722 NoVals = 0 2723 IF( Empty ) RETURN 2724 END IF 2725 2726 ! Save values in the buffer (either single or double prec.) 2727 NoVals = NoVals + 1 2728 IF( SinglePrec ) THEN 2729 Fvals(NoVals) = val 2730 ELSE 2731 DVals(NoVals) = val 2732 END IF 2733 2734 2735 END SUBROUTINE AscBinRealWrite 2736 2737 2738 ! Write an integer value 2739 !------------------------------------------------- 2740 SUBROUTINE AscBinIntegerWrite( ival, EmptyBuffer ) 2741 2742 INTEGER, PARAMETER :: VtuUnit = 58 2743 INTEGER :: ival 2744 LOGICAL, OPTIONAL :: EmptyBuffer 2745 LOGICAL :: Empty 2746 CHARACTER(LEN=1024) :: Str 2747 2748 IF( VtuUnit == 0 ) THEN 2749 CALL Fatal('AscBinIntegerWrite','Buffer not initialized for writing') 2750 END IF 2751 2752 IF( PRESENT( EmptyBuffer ) ) THEN 2753 Empty = EmptyBuffer 2754 ELSE 2755 Empty = .FALSE. 2756 END IF 2757 2758 IF( AsciiOutput ) THEN 2759 IF( Empty ) RETURN 2760 WRITE( Str, '(" ",I0)') ival 2761 WRITE( VtuUnit ) TRIM(Str) 2762 IF( CalcSum ) Isum = Isum + ival 2763 RETURN 2764 END IF 2765 2766 IF( Empty .OR. INoVals == BufferSize ) THEN 2767 IF( INoVals == 0 ) THEN 2768 RETURN 2769 ELSE 2770 WRITE( VtuUnit ) Ivals(1:INoVals) 2771 IF( CalcSum ) THEN 2772 Icount = Icount + 1 2773 Isum = Isum + SUM( Ivals(1:INoVals) ) 2774 END IF 2775 END IF 2776 INoVals = 0 2777 IF( Empty ) RETURN 2778 END IF 2779 2780 INoVals = INoVals + 1 2781 Ivals(INoVals) = ival 2782 2783 END SUBROUTINE AscBinIntegerWrite 2784 2785 2786 SUBROUTINE AscBinInitNorm(CalcNorm) 2787 LOGICAL :: CalcNorm 2788 2789 CalcSum = CalcNorm 2790 RSum = 0.0_dp 2791 ISum = 0.0_dp 2792 Ssum = 0 2793 Rcount = 0 2794 Icount = 0 2795 Scount = 0 2796 2797 END SUBROUTINE AscBinInitNorm 2798 2799 2800 FUNCTION AscBinCompareNorm(RefResults) RESULT ( RelativeNorm ) 2801 REAL(KIND=dp), DIMENSION(*) :: RefResults 2802 REAL(KIND=dp) :: RelativeNorm 2803 REAL(KIND=dp) :: ThisResults(6) 2804 REAL(KIND=dp) :: c 2805 INTEGER :: i 2806 2807 ThisResults(1) = scount 2808 ThisResults(2) = icount 2809 ThisResults(3) = rcount 2810 ThisResults(4) = ssum 2811 ThisResults(5) = isum ! we use real for isum since it could be huge too! 2812 ThisResults(6) = rsum 2813 2814 PRINT *,'Checksums for file output:' 2815 PRINT *,'RefResults:',NINT(RefResults(1:4)),RefResults(5:6) 2816 PRINT *,'ThisResults:',NINT(ThisResults(1:4)),ThisResults(5:6) 2817 2818 ! We have a special relative pseunonorm that should be one! 2819 RelativeNorm = 0.0 2820 DO i=1,6 2821 c = ThisResults(i)/RefResults(i) 2822 RelativeNorm = RelativeNorm + MAX( c, 1.0_dp /c ) / 6 2823 END DO 2824 2825 END FUNCTION AscBinCompareNorm 2826 2827 2828END MODULE AscBinOutputUtils 2829 2830 2831!> \} 2832 2833