text:
      PROGRAM FM045

C     COMMENT SECTION                                                   00010045
C                                                                       00020045
C     FM045                                                             00030045
C                                                                       00040045
C         THIS ROUTINE TESTS ARITHMETIC ASSIGNMENTS USING INTEGER       00050045
C     VARIABLES CONNECTED BY A SERIES OF ARITHMETIC OPERATORS.          00060045
C     DIFFERENT COMBINATIONS OF PARENTHETICAL NOTATION ARE EXERCIZED.   00070045
C                                                                       00080045
C                                                                       00090045
C      REFERENCES                                                       00100045
C        AMERICAN NATIONAL STANDARD PROGRAMMING LANGUAGE FORTRAN,       00110045
C              X3.9-1978                                                00120045
C                                                                       00130045
C        SECTION 4.3, INTEGER TYPE                                      00140045
C        SECTION 4.3.1, INTEGER CONSTANT                                00150045
C        SECTION 6.1, ARITHMETIC EXPRESSIONS                            00160045
C        SECTION 6.6, EVALUATION OF EXPRESSIONS                         00170045
C        SECTION 10.1, ARITHMETIC ASSIGNMENT STATEMENT                  00180045
C                                                                       00190045
C                                                                       00200045
C                                                                       00210045
C      **********************************************************       00220045
C                                                                       00230045
C         A COMPILER VALIDATION SYSTEM FOR THE FORTRAN LANGUAGE         00240045
C     BASED ON SPECIFICATIONS AS DEFINED IN AMERICAN NATIONAL STANDARD  00250045
C     PROGRAMMING LANGUAGE FORTRAN X3.9-1978, HAS BEEN DEVELOPED BY THE 00260045
C     FEDERAL COBOL COMPILER TESTING SERVICE.  THE FORTRAN COMPILER     00270045
C     VALIDATION SYSTEM (FCVS) CONSISTS OF AUDIT ROUTINES, THEIR RELATED00280045
C     DATA, AND AN EXECUTIVE SYSTEM.  EACH AUDIT ROUTINE IS A FORTRAN   00290045
C     PROGRAM, SUBPROGRAM OR FUNCTION WHICH INCLUDES TESTS OF SPECIFIC  00300045
C     LANGUAGE ELEMENTS AND SUPPORTING PROCEDURES INDICATING THE RESULT 00310045
C     OF EXECUTING THESE TESTS.                                         00320045
C                                                                       00330045
C         THIS PARTICULAR PROGRAM/SUBPROGRAM/FUNCTION CONTAINS FEATURES 00340045
C     FOUND ONLY IN THE SUBSET AS DEFINED IN X3.9-1978.                 00350045
C                                                                       00360045
C         SUGGESTIONS AND COMMENTS SHOULD BE FORWARDED TO -             00370045
C                                                                       00380045
C              NATIONAL INSTITUTE OF STANDARDS AND TECHNOLOGY           00390045
C                   SOFTWARE STANDARDS VALIDATION GROUP                 00400045
C                          BUILDING 225  RM A266                        00410045
C                         GAITHERSBURG, MD  20899                       00420045
C      **********************************************************       00430045
C                                                                       00440045
C                                                                       00450045
C                                                                       00460045
C     INITIALIZATION SECTION                                            00470045
C                                                                       00480045
C     INITIALIZE CONSTANTS                                              00490045
C      **************                                                   00500045
C     I01 CONTAINS THE LOGICAL UNIT NUMBER FOR THE CARD READER.         00510045
      I01 = 5                                                           00520045
C     I02 CONTAINS THE LOGICAL UNIT NUMBER FOR THE PRINTER.             00530045
      I02 = 6                                                           00540045
C     SYSTEM ENVIRONMENT SECTION                                        00550045
C                                                                       00560045
CX010    THIS CARD IS REPLACED BY CONTENTS OF FEXEC X-010 CONTROL CARD. 00570045
C     THE CX010 CARD IS FOR OVERRIDING THE PROGRAM DEFAULT I01 = 5      00580045
C     (UNIT NUMBER FOR CARD READER).                                    00590045
CX011    THIS CARD IS REPLACED BY CONTENTS OF FEXEC X-011 CONTROL CARD. 00600045
C     THE CX011 CARD IS FOR SYSTEMS WHICH REQUIRE ADDITIONAL            00610045
C     FORTRAN STATEMENTS FOR FILES ASSOCIATED WITH CX010 ABOVE.         00620045
C                                                                       00630045
CX020    THIS CARD IS REPLACED BY CONTENTS OF FEXEC X-020 CONTROL CARD. 00640045
C     THE CX020 CARD IS FOR OVERRIDING THE PROGRAM DEFAULT I02 = 6      00650045
C     (UNIT NUMBER FOR PRINTER).                                        00660045
CX021    THIS CARD IS REPLACED BY CONTENTS OF FEXEC X-021 CONTROL CARD. 00670045
C     THE CX021 CARD IS FOR SYSTEMS WHICH REQUIRE ADDITIONAL            00680045
C     FORTRAN STATEMENTS FOR FILES ASSOCIATED WITH CX020 ABOVE.         00690045
C                                                                       00700045
      IVPASS=0                                                          00710045
      IVFAIL=0                                                          00720045
      IVDELE=0                                                          00730045
      ICZERO=0                                                          00740045
C                                                                       00750045
C     WRITE PAGE HEADERS                                                00760045
      WRITE (I02,90000)                                                 00770045
      WRITE (I02,90001)                                                 00780045
      WRITE (I02,90002)                                                 00790045
      WRITE (I02, 90002)                                                00800045
      WRITE (I02,90003)                                                 00810045
      WRITE (I02,90002)                                                 00820045
      WRITE (I02,90004)                                                 00830045
      WRITE (I02,90002)                                                 00840045
      WRITE (I02,90011)                                                 00850045
      WRITE (I02,90002)                                                 00860045
      WRITE (I02,90002)                                                 00870045
      WRITE (I02,90005)                                                 00880045
      WRITE (I02,90006)                                                 00890045
      WRITE (I02,90002)                                                 00900045
C                                                                       00910045
C                                                                       00920045
C     TEST SECTION                                                      00930045
C                                                                       00940045
C         ARITHMETIC ASSIGNMENT STATEMENT                               00950045
C                                                                       00960045
C                                                                       00970045
C     TESTS 747 THROUGH 755 USE THE SAME STRING OF VARIABLES AND        00980045
C     OPERATORS, BUT USE DIFFERENT COMBINATIONS OF PARENTHETICAL        00990045
C     NOTATION  TO ALTER PRIORITIES IN ORDER OF EVALUATION.             01000045
C                                                                       01010045
C     TESTS 756 THROUGH 759 CHECK THE CAPABILITY TO ENCLOSE THE ENTIRE  01020045
C     RIGHT HAND SIDE OF AN ASSIGNMENT STATEMENT IN PARENTHESES OR SETS 01030045
C     OF NESTED PARENTHESES.                                            01040045
C                                                                       01050045
C                                                                       01060045
C                                                                       01070045
C                                                                       01080045
C                                                                       01090045
C                                                                       01100045
C                                                                       01110045
      IVTNUM = 747                                                      01120045
C                                                                       01130045
C      ****  TEST 747  ****                                             01140045
C                                                                       01150045
      IF (ICZERO) 37470, 7470, 37470                                    01160045
 7470 CONTINUE                                                          01170045
      IVON01 = 15                                                       01180045
      IVON02 =  9                                                       01190045
      IVON03 =  4                                                       01200045
      IVON04 = 18                                                       01210045
      IVON05 =  6                                                       01220045
      IVON06 =  2                                                       01230045
      IVCOMP = IVON01 + IVON02 - IVON03 * IVON04 / IVON05 ** IVON06     01240045
      GO TO 47470                                                       01250045
37470 IVDELE = IVDELE + 1                                               01260045
      WRITE (I02,80003) IVTNUM                                          01270045
      IF (ICZERO) 47470, 7481, 47470                                    01280045
47470 IF (IVCOMP - 22) 27470,17470,27470                                01290045
17470 IVPASS = IVPASS + 1                                               01300045
      WRITE (I02,80001) IVTNUM                                          01310045
      GO TO 7481                                                        01320045
27470 IVFAIL = IVFAIL + 1                                               01330045
      IVCORR = 22                                                       01340045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          01350045
 7481 CONTINUE                                                          01360045
      IVTNUM = 748                                                      01370045
C                                                                       01380045
C      ****  TEST 748  ****                                             01390045
C                                                                       01400045
      IF (ICZERO) 37480, 7480, 37480                                    01410045
 7480 CONTINUE                                                          01420045
      IVON01 = 15                                                       01430045
      IVON02 =  9                                                       01440045
      IVON03 =  4                                                       01450045
      IVON04 = 18                                                       01460045
      IVON05 =  6                                                       01470045
      IVON06 =  2                                                       01480045
      IVCOMP = ((((IVON01 + IVON02) - IVON03) * IVON04) / IVON05)       01490045
     *         ** IVON06                                                01500045
      GO TO 47480                                                       01510045
37480 IVDELE = IVDELE + 1                                               01520045
      WRITE (I02,80003) IVTNUM                                          01530045
      IF (ICZERO) 47480, 7491, 47480                                    01540045
47480 IF (IVCOMP - 3600) 27480,17480,27480                              01550045
17480 IVPASS = IVPASS + 1                                               01560045
      WRITE (I02,80001) IVTNUM                                          01570045
      GO TO 7491                                                        01580045
27480 IVFAIL = IVFAIL + 1                                               01590045
      IVCORR = 3600                                                     01600045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          01610045
 7491 CONTINUE                                                          01620045
      IVTNUM = 749                                                      01630045
C                                                                       01640045
C      ****  TEST 749  ****                                             01650045
C                                                                       01660045
      IF (ICZERO) 37490, 7490, 37490                                    01670045
 7490 CONTINUE                                                          01680045
      IVON01 = 15                                                       01690045
      IVON02 =  9                                                       01700045
      IVON03 =  4                                                       01710045
      IVON04 = 36                                                       01720045
      IVON05 =  6                                                       01730045
      IVON06 =  2                                                       01740045
      IVCOMP = (IVON01 + IVON02 - IVON03) * (IVON04 / IVON05 ** IVON06) 01750045
      GO TO 47490                                                       01760045
37490 IVDELE = IVDELE + 1                                               01770045
      WRITE (I02,80003) IVTNUM                                          01780045
      IF (ICZERO) 47490, 7501, 47490                                    01790045
47490 IF (IVCOMP - 20) 27490,17490,27490                                01800045
17490 IVPASS = IVPASS + 1                                               01810045
      WRITE (I02,80001) IVTNUM                                          01820045
      GO TO 7501                                                        01830045
27490 IVFAIL = IVFAIL + 1                                               01840045
      IVCORR = 20                                                       01850045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          01860045
 7501 CONTINUE                                                          01870045
      IVTNUM = 750                                                      01880045
C                                                                       01890045
C      ****  TEST 750  ****                                             01900045
C                                                                       01910045
      IF (ICZERO) 37500, 7500, 37500                                    01920045
 7500 CONTINUE                                                          01930045
      IVON01 = 15                                                       01940045
      IVON02 =  9                                                       01950045
      IVON03 =  4                                                       01960045
      IVON04 = 36                                                       01970045
      IVON05 =  6                                                       01980045
      IVON06 =  2                                                       01990045
      IVCOMP = (IVON01 + IVON02) - (IVON03 * IVON04) / (IVON05 **       02000045
     *         IVON06)                                                  02010045
      GO TO 47500                                                       02020045
37500 IVDELE = IVDELE + 1                                               02030045
      WRITE (I02,80003) IVTNUM                                          02040045
      IF (ICZERO) 47500, 7511, 47500                                    02050045
47500 IF (IVCOMP - 20) 27500,17500,27500                                02060045
17500 IVPASS = IVPASS + 1                                               02070045
      WRITE (I02,80001) IVTNUM                                          02080045
      GO TO 7511                                                        02090045
27500 IVFAIL = IVFAIL + 1                                               02100045
      IVCORR = 20                                                       02110045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          02120045
 7511 CONTINUE                                                          02130045
      IVTNUM = 751                                                      02140045
C                                                                       02150045
C      ****  TEST 751  ****                                             02160045
C                                                                       02170045
      IF (ICZERO) 37510, 7510, 37510                                    02180045
 7510 CONTINUE                                                          02190045
      IVON01 = 15                                                       02200045
      IVON02 =  9                                                       02210045
      IVON03 =  4                                                       02220045
      IVON04 = 36                                                       02230045
      IVON05 =  6                                                       02240045
      IVON06 =  2                                                       02250045
      IVCOMP = ((IVON01 + IVON02) - (IVON03 * IVON04)) / (IVON05 **     02260045
     *         IVON06)                                                  02270045
      GO TO 47510                                                       02280045
37510 IVDELE = IVDELE + 1                                               02290045
      WRITE (I02,80003) IVTNUM                                          02300045
      IF (ICZERO)  47510, 7521, 47510                                   02310045
47510 IF (IVCOMP + 3)  27510,17510,27510                                02320045
17510 IVPASS = IVPASS + 1                                               02330045
      WRITE (I02,80001) IVTNUM                                          02340045
      GO TO 7521                                                        02350045
27510 IVFAIL = IVFAIL + 1                                               02360045
      IVCORR = -3                                                       02370045
C     ACTUAL ANSWER IS  -3.333333...     TRUNCATION IS NECESSARY        02380045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          02390045
 7521 CONTINUE                                                          02400045
      IVTNUM = 752                                                      02410045
C                                                                       02420045
C      ****  TEST 752  ****                                             02430045
C                                                                       02440045
      IF (ICZERO) 37520, 7520, 37520                                    02450045
 7520 CONTINUE                                                          02460045
      IVON01 = 15                                                       02470045
      IVON02 =  9                                                       02480045
      IVON03 =  4                                                       02490045
      IVON04 = 36                                                       02500045
      IVON05 =  6                                                       02510045
      IVON06 =  2                                                       02520045
      IVCOMP = (IVON01 + IVON02) - (IVON03 * IVON04 / IVON05) ** IVON06 02530045
      GO TO 47520                                                       02540045
37520 IVDELE = IVDELE + 1                                               02550045
      WRITE (I02,80003) IVTNUM                                          02560045
      IF (ICZERO) 47520, 7531, 47520                                    02570045
47520 IF (IVCOMP + 552) 27520,17520,27520                               02580045
17520 IVPASS = IVPASS + 1                                               02590045
      WRITE (I02,80001) IVTNUM                                          02600045
      GO TO 7531                                                        02610045
27520 IVFAIL = IVFAIL + 1                                               02620045
      IVCORR = -552                                                     02630045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          02640045
 7531 CONTINUE                                                          02650045
      IVTNUM = 753                                                      02660045
C                                                                       02670045
C      ****  TEST 753  ****                                             02680045
C                                                                       02690045
      IF (ICZERO) 37530, 7530, 37530                                    02700045
 7530 CONTINUE                                                          02710045
      IVON01 = 15                                                       02720045
      IVON02 =  9                                                       02730045
      IVON03 =  4                                                       02740045
      IVON04 = 36                                                       02750045
      IVON05 =  6                                                       02760045
      IVON06 =  2                                                       02770045
      IVCOMP = IVON01 + (IVON02 - IVON03 * IVON04) / IVON05 ** IVON06   02780045
      GO TO 47530                                                       02790045
37530 IVDELE = IVDELE + 1                                               02800045
      WRITE (I02,80003) IVTNUM                                          02810045
      IF (ICZERO) 47530, 7541, 47530                                    02820045
47530 IF (IVCOMP - 12) 27530,17530,27530                                02830045
17530 IVPASS = IVPASS + 1                                               02840045
      WRITE (I02,80001) IVTNUM                                          02850045
      GO TO 7541                                                        02860045
27530 IVFAIL = IVFAIL + 1                                               02870045
      IVCORR = 12                                                       02880045
C     ACTUAL ANSWER IS  11.25            TRUNCATION IS NECESSARY        02890045
C                                        DURING AN INTERMEDIATE STEP    02900045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          02910045
 7541 CONTINUE                                                          02920045
      IVTNUM = 754                                                      02930045
C                                                                       02940045
C      ****  TEST 754  ****                                             02950045
C                                                                       02960045
      IF (ICZERO) 37540, 7540, 37540                                    02970045
 7540 CONTINUE                                                          02980045
      IVON01 = 15                                                       02990045
      IVON02 =  9                                                       03000045
      IVON03 =  4                                                       03010045
      IVON04 = 36                                                       03020045
      IVON05 =  6                                                       03030045
      IVON06 =  2                                                       03040045
      IVCOMP = IVON01 + (IVON02 - IVON03) * (IVON04 / IVON05) ** IVON06 03050045
      GO TO 47540                                                       03060045
37540 IVDELE = IVDELE + 1                                               03070045
      WRITE (I02,80003) IVTNUM                                          03080045
      IF (ICZERO) 47540, 7551, 47540                                    03090045
47540 IF (IVCOMP - 195) 27540,17540,27540                               03100045
17540 IVPASS = IVPASS + 1                                               03110045
      WRITE (I02,80001) IVTNUM                                          03120045
      GO TO 7551                                                        03130045
27540 IVFAIL = IVFAIL + 1                                               03140045
      IVCORR = 195                                                      03150045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          03160045
 7551 CONTINUE                                                          03170045
      IVTNUM = 755                                                      03180045
C                                                                       03190045
C      ****  TEST 755  ****                                             03200045
C                                                                       03210045
      IF (ICZERO) 37550, 7550, 37550                                    03220045
 7550 CONTINUE                                                          03230045
      IVON01 = 15                                                       03240045
      IVON02 =  9                                                       03250045
      IVON03 =  4                                                       03260045
      IVON04 = 36                                                       03270045
      IVON05 =  6                                                       03280045
      IVON06 =  2                                                       03290045
      IVCOMP = ((IVON01 + (IVON02 - IVON03) * IVON04) / IVON05) **      03300045
     *         IVON06                                                   03310045
      GO TO 47550                                                       03320045
37550 IVDELE = IVDELE + 1                                               03330045
      WRITE (I02,80003) IVTNUM                                          03340045
      IF (ICZERO) 47550, 7561, 47550                                    03350045
47550 IF (IVCOMP - 1024)  27550,17550,27550                             03360045
17550 IVPASS = IVPASS + 1                                               03370045
      WRITE (I02,80001) IVTNUM                                          03380045
      GO TO 7561                                                        03390045
27550 IVFAIL = IVFAIL + 1                                               03400045
      IVCORR = 1024                                                     03410045
C     ACTUAL ANSWER IS  1056.25         TRUNCATION IS NECESSARY         03420045
C                                       DURING AN INTERMEDIATE STEP     03430045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          03440045
 7561 CONTINUE                                                          03450045
      IVTNUM = 756                                                      03460045
C                                                                       03470045
C      ****  TEST 756  ****                                             03480045
C          SINGLE PARENTHESES                                           03490045
C                                                                       03500045
      IF (ICZERO) 37560, 7560, 37560                                    03510045
 7560 CONTINUE                                                          03520045
      IVON01 = 13                                                       03530045
      IVON02 = 37                                                       03540045
      IVCOMP = (IVON01 + IVON02)                                        03550045
      GO TO 47560                                                       03560045
37560 IVDELE = IVDELE + 1                                               03570045
      WRITE (I02,80003) IVTNUM                                          03580045
      IF (ICZERO) 47560, 7571, 47560                                    03590045
47560 IF (IVCOMP - 50) 27560,17560,27560                                03600045
17560 IVPASS = IVPASS + 1                                               03610045
      WRITE (I02,80001) IVTNUM                                          03620045
      GO TO 7571                                                        03630045
27560 IVFAIL = IVFAIL + 1                                               03640045
      IVCORR = 50                                                       03650045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          03660045
 7571 CONTINUE                                                          03670045
      IVTNUM = 757                                                      03680045
C                                                                       03690045
C      ****  TEST 757  ****                                             03700045
C          NESTED PARENTHESES (TWO SETS)                                03710045
C                                                                       03720045
      IF (ICZERO) 37570, 7570, 37570                                    03730045
 7570 CONTINUE                                                          03740045
      IVON01 = 13                                                       03750045
      IVON02 = 37                                                       03760045
      IVCOMP = ((IVON01 - IVON02))                                      03770045
      GO TO 47570                                                       03780045
37570 IVDELE = IVDELE + 1                                               03790045
      WRITE (I02,80003) IVTNUM                                          03800045
      IF (ICZERO) 47570, 7581, 47570                                    03810045
47570 IF (IVCOMP + 24) 27570,17570,27570                                03820045
17570 IVPASS = IVPASS + 1                                               03830045
      WRITE (I02,80001) IVTNUM                                          03840045
      GO TO 7581                                                        03850045
27570 IVFAIL = IVFAIL + 1                                               03860045
      IVCORR = -24                                                      03870045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          03880045
 7581 CONTINUE                                                          03890045
      IVTNUM = 758                                                      03900045
C                                                                       03910045
C      ****  TEST 758  ****                                             03920045
C          NESTED PARENTHESES (21 SETS - SAME LINE)                     03930045
C                                                                       03940045
      IF (ICZERO) 37580, 7580, 37580                                    03950045
 7580 CONTINUE                                                          03960045
      IVON01 = 13                                                       03970045
      IVON02 = 37                                                       03980045
      IVCOMP = (((((((((((((((((((((IVON01 * IVON02)))))))))))))))))))))03990045
      GO TO 47580                                                       04000045
37580 IVDELE = IVDELE + 1                                               04010045
      WRITE (I02,80003) IVTNUM                                          04020045
      IF (ICZERO) 47580, 7591, 47580                                    04030045
47580 IF (IVCOMP - 481) 27580,17580,27580                               04040045
17580 IVPASS = IVPASS + 1                                               04050045
      WRITE (I02,80001) IVTNUM                                          04060045
      GO TO 7591                                                        04070045
27580 IVFAIL = IVFAIL + 1                                               04080045
      IVCORR = 481                                                      04090045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          04100045
 7591 CONTINUE                                                          04110045
      IVTNUM = 759                                                      04120045
C                                                                       04130045
C      ****  TEST 759  ****                                             04140045
C          NESTED PARENTHESES (57 SETS - MULTIPLE LINES)                04150045
C                                                                       04160045
      IF (ICZERO) 37590, 7590, 37590                                    04170045
 7590 CONTINUE                                                          04180045
      IVON01 = 13                                                       04190045
      IVON02 = 37                                                       04200045
      IVCOMP = (((((((((((((((((((((((((((((((((((((((((((((((((((((((((04210045
     *         IVON01 / IVON02                                          04220045
     *         )))))))))))))))))))))))))))))))))))))))))))))))))))))))))04230045
      GO TO 47590                                                       04240045
37590 IVDELE = IVDELE + 1                                               04250045
      WRITE (I02,80003) IVTNUM                                          04260045
      IF (ICZERO) 47590, 7601, 47590                                    04270045
47590 IF (IVCOMP) 27590,17590,27590                                     04280045
17590 IVPASS = IVPASS + 1                                               04290045
      WRITE (I02,80001) IVTNUM                                          04300045
      GO TO 7601                                                        04310045
27590 IVFAIL = IVFAIL + 1                                               04320045
      IVCORR = 0                                                        04330045
      WRITE (I02,80004) IVTNUM, IVCOMP, IVCORR                          04340045
 7601 CONTINUE                                                          04350045
C                                                                       04360045
C     WRITE PAGE FOOTINGS AND RUN SUMMARIES                             04370045
99999 CONTINUE                                                          04380045
      WRITE (I02,90002)                                                 04390045
      WRITE (I02,90006)                                                 04400045
      WRITE (I02,90002)                                                 04410045
      WRITE (I02,90002)                                                 04420045
      WRITE (I02,90007)                                                 04430045
      WRITE (I02,90002)                                                 04440045
      WRITE (I02,90008)  IVFAIL                                         04450045
      WRITE (I02,90009) IVPASS                                          04460045
      WRITE (I02,90010) IVDELE                                          04470045
C                                                                       04480045
C                                                                       04490045
C     TERMINATE ROUTINE EXECUTION                                       04500045
      STOP                                                              04510045
C                                                                       04520045
C     FORMAT STATEMENTS FOR PAGE HEADERS                                04530045
90000 FORMAT ("1")                                                      04540045
90002 FORMAT (" ")                                                      04550045
90001 FORMAT (" ",10X,"FORTRAN COMPILER VALIDATION SYSTEM" )            04560045
90003 FORMAT (" ",21X,"VERSION 2.1" )                                   04570045
90004 FORMAT (" ",10X,"FOR OFFICIAL USE ONLY - COPYRIGHT 1978" )        04580045
90005 FORMAT (" ",5X,"TEST",5X,"PASS/FAIL", 5X,"COMPUTED",8X,"CORRECT") 04590045
90006 FORMAT (" ",5X,"----------------------------------------------" ) 04600045
90011 FORMAT (" ",18X,"SUBSET LEVEL TEST" )                             04610045
C                                                                       04620045
C     FORMAT STATEMENTS FOR RUN SUMMARIES                               04630045
90008 FORMAT (" ",15X,I5," ERRORS ENCOUNTERED" )                        04640045
90009 FORMAT (" ",15X,I5," TESTS PASSED" )                              04650045
90010 FORMAT (" ",15X,I5," TESTS DELETED" )                             04660045
C                                                                       04670045
C     FORMAT STATEMENTS FOR TEST RESULTS                                04680045
80001 FORMAT (" ",4X,I5,7X,"PASS")                                      04690045
80002 FORMAT (" ",4X,I5,7X,"FAIL")                                      04700045
80003 FORMAT (" ",4X,I5,7X,"DELETED")                                   04710045
80004 FORMAT (" ",4X,I5,7X,"FAIL",10X,I6,9X,I6)                         04720045
80005 FORMAT (" ",4X,I5,7X,"FAIL",4X,E12.5,3X,E12.5)                    04730045
C                                                                       04740045
90007 FORMAT (" ",20X,"END OF PROGRAM FM045" )                          04750045
      END                                                               04760045

max_stars_repo_name: MrSidims/llvm-test-suite

max_issues_repo_name: MrSidims/llvm-test-suite

max_forks_repo_name: MrSidims/llvm-test-suite


text:
      subroutine zsphbesjh(lmax,x,bj,hm,eiz,job)
c_______________________________________________________________________
c calculate spherical bessel and hankel functions
c input: lmax integer scalar, max l
c        x    complex*16 scalar, argument of the bessel functions
c        job  integer scalar. job=+1: retarded hankel function is calculated
c             job=-1: advanced hankel function is calculated
c returns: bj complex*16 array of (0:lmax), j_l
c          hm complex*16 array of (0:lmax), h(+)_l (retarded hankel function)
c             this is scaled out of a factor of exp(i*x)/x
c          eiz complex*16 scalar, factor exp(i*x)
c calls zsphbesj
c  computes spherical bessel/hankel functions j and h+
c    h+ = jl + i nl = sinx/x - i cosx/x
c    h- = jl - i nl = sinx/x + i cosx/x
c
c   xgz ornl 1994
      implicit real*8 (a-h,o-z)
      complex*16 ci
      parameter(ci=(0.0d0,1.0d0))
      complex*16 bj(0:lmax),hm(0:lmax)
      complex*16 x,x1,eiz,fll,ac,as
c
      if(x.eq.0.0d0) goto 200

c use hm for scratch in zsphbesj
c zsphbesj returns x1=1/x, as=sin(x),ac=cos(x)
      call zsphbesj(lmax,x,as,ac,x1,bj,hm)
c pull out factor of e^ix/x from h
c eiz = job*exp(job*ci*x)
      eiz=job*exp(job*ci*x)

c now use hm for hankel hm = jl + i nl
      hm(0)=-ci
c hm(1) = -1.d0 - ci*x1
      hm(1)=-dcmplx(1.0d0-dimag(x1),dreal(x1))

      do l=2,lmax
        hm(l)=(2*l-1)*hm(l-1)*x1-hm(l-2)
      enddo
      return

200   write(6,*)'error bessel: zero argument.'
      stop
      end

ext: f

max_stars_repo_name: rdietric/lsms

max_issues_repo_name: rdietric/lsms

max_forks_repo_name: hkershaw-brown/MuST


text:
c=======================================================================

      subroutine write_illuminated_earth(day_night)

      implicit none
      include 'parameters'
      include 'constants'
      double precision alpha_grid, delta_grid
      double precision day_night(1:n_alpha, 1:n_delta)
      integer i,j

      open(unit = 21, file = '../OUTPUT/illum_earth.out')

      do i = 1, n_alpha
         alpha_grid = 2.d0*pi*dble(i-1)/dble(n_alpha-1)

         do j = 1, n_delta
            delta_grid = -pi/2.d0 + dble(j-1)*pi/dble(n_delta-1)

            write(21,*) alpha_grid*RAD, delta_grid*RAD,
     &           day_night(i,j)*albedo
            write(21,*)
         enddo
      enddo

      close(21)

      return
      end

c=======================================================================

      subroutine write_seen(day_night)

      implicit none
      include 'parameters'
      include 'constants'
      double precision alpha_grid, delta_grid
      double precision day_night(1:n_alpha, 1:n_delta)
      integer i,j

      open(unit = 21, file = '../OUTPUT/seen.out')

      do i = 1, n_alpha
         alpha_grid = 2.d0*pi*dble(i-1)/dble(n_alpha-1)

         do j = 1, n_delta
            delta_grid = -pi/2.d0 + dble(j-1)*pi/dble(n_delta-1)

            write(21,*) alpha_grid*RAD, delta_grid*RAD,
     &           day_night(i,j)*albedo
            write(21,*)
         enddo
      enddo

      close(21)

      return
      end

c=======================================================================

      subroutine write_angles(day_night)

      implicit none
      include 'parameters'
      include 'constants'
      double precision alpha_grid, delta_grid
      double precision day_night(1:n_alpha, 1:n_delta)
      integer i,j

      open(unit = 21, file = '../OUTPUT/angles.out')

      do i = 1, n_alpha
         alpha_grid = 2.d0*pi*dble(i-1)/dble(n_alpha-1)

         do j = 1, n_delta
            delta_grid = -pi/2.d0 + dble(j-1)*pi/dble(n_delta-1)

            write(21,*) alpha_grid*RAD, delta_grid*RAD,
     &           day_night(i,j)*albedo
            write(21,*)
         enddo
      enddo

      close(21)

      return
      end

c=======================================================================

      subroutine write_contamination(day_night)

      implicit none
      include 'parameters'
      include 'constants'
      double precision alpha_grid, delta_grid
      double precision day_night(1:n_alpha, 1:n_delta)
      integer i,j

      open(unit = 21, file = '../OUTPUT/contamination.out')

      do i = 1, n_alpha
         alpha_grid = 2.d0*pi*dble(i-1)/dble(n_alpha-1)

         do j = 1, n_delta
            delta_grid = -pi/2.d0 + dble(j-1)*pi/dble(n_delta-1)

            write(21,*) alpha_grid*RAD, delta_grid*RAD,
     &           day_night(i,j)*albedo
            write(21,*)
         enddo
      enddo

      close(21)

      return
      end

c=======================================================================

      subroutine write_star(ra,dec)

      implicit none
      include 'parameters'
      include 'constants'
      double precision ra,dec

      open(unit = 21, file = '../OUTPUT/star.out')

      write(21,*) ra, dec

      close(21)

      return
      end

c=======================================================================

      subroutine write_error_targets(ra,dec)

      implicit none
      include 'parameters'
      include 'constants'
      double precision ra,dec, JD

      call READ_DATE_IN_JD(JD)

      open(unit = 21, file = '../OUTPUT/errors.out', ACCESS='APPEND')

      write(21,*) JD, ra, dec

      close(21)

      return
      end

max_stars_repo_name: kuntzer/SALSA-public

max_issues_repo_name: kuntzer/SALSA-public

max_forks_repo_name: kuntzer/SALSA-public


text:
      FUNCTION LJUST(STR)

C$$$  SUBPROGRAM DOCUMENTATION BLOCK
C
C SUBPROGRAM:    LJUST
C   PRGMMR: WOOLLEN          ORG: NP20       DATE: 1998-07-08
C
C ABSTRACT: THIS FUNCTION REMOVES ALL LEADING BLANKS FROM A CHARACTER
C   STRING (LEFT JUSTIFIES).  IT IS CONSIDERED OBSOLETE AND MAY BE
C   REMOVED FROM THE BUFR ARCHIVE LIBRARY IN A FUTURE VERSION.  USERS
C   SHOULD INSTEAD MIGRATE TO THE USE OF SUBROUTINE JSTCHR.
C
C PROGRAM HISTORY LOG:
C 1998-07-08  J. WOOLLEN -- ORIGINAL AUTHOR
C 2003-11-04  S. BENDER  -- ADDED REMARKS/BUFRLIB ROUTINE
C                           INTERDEPENDENCIES
C 2003-11-04  D. KEYSER  -- UNIFIED/PORTABLE FOR WRF; ADDED
C                           DOCUMENTATION (INCLUDING HISTORY)
C 2007-01-19  J. ATOR    -- MARKED AS OBSOLETE AND ADDED PRINT
C                           NOTIFICATION
C
C USAGE:    LJUST (STR)
C   INPUT ARGUMENT LIST:
C     STR      - CHARACTER*(*): STRING TO BE LEFT-JUSTIFED
C
C   OUTPUT ARGUMENT LIST:
C     STR      - CHARACTER*(*): LEFT-JUSTIFED STRING (UNLESS INPUT
C                VALUE CONTAINED ALL BLANKS - THEN INPUT VALUE IS
C                RETURNED WITH NO CHANGES)
C     LJUST    - INTEGER: ALWAYS RETURNED AS 0 (DUMMY)
C
C REMARKS:
C    THIS ROUTINE CALLS:        None
C    THIS ROUTINE IS CALLED BY: None at this time
C                               Normally not called by any application
C                               Programs but it could be.
C
C ATTRIBUTES:
C   LANGUAGE: FORTRAN 77
C   MACHINE:  PORTABLE TO ALL PLATFORMS
C
C$$$

      CHARACTER*(*) STR

      COMMON /QUIET / IPRT

      DATA IFIRST/0/

      SAVE IFIRST

C-----------------------------------------------------------------------
C-----------------------------------------------------------------------

      IF(IFIRST.EQ.0) THEN
         IF(IPRT.GE.0) THEN
      PRINT*
      PRINT*,'+++++++++++++++++BUFR ARCHIVE LIBRARY++++++++++++++++++++'
      PRINT 101
101   FORMAT(' BUFRLIB: LJUST - THIS FUNCTION IS NOW OBSOLETE; ',
     . 'USE SUBROUTINE JSTCHR INSTEAD')
      PRINT*,'+++++++++++++++++BUFR ARCHIVE LIBRARY++++++++++++++++++++'
      PRINT*
         ENDIF
         IFIRST = 1
      ENDIF

      LJUST = 0

      IF(STR.EQ.' ') GOTO 100

      LSTR = LEN(STR)
      DO I=1,LSTR
      DO WHILE(STR(I:I).EQ.' ' .AND. STR(I+1:LSTR).NE.' ')
         STR(I:LSTR) = STR(I+1:LSTR)
      ENDDO
      ENDDO

C  EXIT
C  ----

100   RETURN
      END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


    Program zgeevx_example

!     ZGEEVX Example Program Text

!     Copyright (c) 2018, Numerical Algorithms Group (NAG Ltd.)
!     For licence see
!       https://github.com/numericalalgorithmsgroup/LAPACK_Examples/blob/master/LICENCE.md

!     .. Use Statements ..
      Use lapack_interfaces, Only: zgeevx
      Use lapack_precision, Only: dp
!     .. Implicit None Statement ..
      Implicit None
!     .. Parameters ..
      Integer, Parameter :: nb = 64, nin = 5, nout = 6
!     .. Local Scalars ..
      Real (Kind=dp) :: abnrm, eps, tol
      Integer :: i, ihi, ilo, info, j, lda, ldvl, ldvr, lwork, n
!     .. Local Arrays ..
      Complex (Kind=dp), Allocatable :: a(:, :), vl(:, :), vr(:, :), w(:), &
        work(:)
      Complex (Kind=dp) :: dummy(1)
      Real (Kind=dp), Allocatable :: rconde(:), rcondv(:), rwork(:), scale(:)
!     .. Intrinsic Procedures ..
      Intrinsic :: epsilon, max, nint, real
!     .. Executable Statements ..
      Write (nout, *) 'ZGEEVX Example Program Results'
!     Skip heading in data file
      Read (nin, *)
      Read (nin, *) n
      lda = n
      ldvl = n
      ldvr = n
      Allocate (a(lda,n), vl(ldvl,n), vr(ldvr,n), w(n), rconde(n), rcondv(n), &
        rwork(2*n), scale(n))

!     Use routine workspace query to get optimal workspace.
      lwork = -1
      Call zgeevx('Balance', 'Vectors (left)', 'Vectors (right)', &
        'Both reciprocal condition numbers', n, a, lda, w, vl, ldvl, vr, ldvr, &
        ilo, ihi, scale, abnrm, rconde, rcondv, dummy, lwork, rwork, info)

!     Make sure that there is enough workspace for block size nb.
      lwork = max((nb+1)*n, nint(real(dummy(1))))
      Allocate (work(lwork))

!     Read the matrix A from data file

      Read (nin, *)(a(i,1:n), i=1, n)

!     Solve the eigenvalue problem

      Call zgeevx('Balance', 'Vectors (left)', 'Vectors (right)', &
        'Both reciprocal condition numbers', n, a, lda, w, vl, ldvl, vr, ldvr, &
        ilo, ihi, scale, abnrm, rconde, rcondv, work, lwork, rwork, info)

      If (info==0) Then

!       Compute the machine precision

        eps = epsilon(1.0E0_dp)
        tol = eps*abnrm

!       Print the eigenvalues and vectors, and associated condition
!       number and bounds

        Write (nout, *)
        Write (nout, *) 'Eigenvalues'
        Write (nout, *)
        Write (nout, *) '         Eigenvalue           rcond    error'

        Do j = 1, n

!         Print information on j-th eigenvalue

          If (rconde(j)>0.0_dp) Then
            If (tol/rconde(j)<10.0_dp*eps) Then
              Write (nout, 100) j, w(j), rconde(j), '-'
            Else
              Write (nout, 110) j, w(j), rconde(j), tol/rconde(j)
            End If
          Else
            Write (nout, 100) j, w(j), rconde(j), 'Inf'
          End If

        End Do

        Write (nout, *)
        Write (nout, *) 'Eigenvectors'
        Write (nout, *)
        Write (nout, *) '         Eigenvector          rcond    error'

        Do j = 1, n

!         Print information on j-th eigenvector

          Write (nout, *)

!         Make first real part component be positive
          If (real(vr(1,j))<0.0_dp) Then
            vr(1:n, j) = -vr(1:n, j)
          End If
          If (rcondv(j)>0.0_dp) Then
            If (tol/rcondv(j)<10.0_dp*eps) Then
              Write (nout, 100) j, vr(1, j), rcondv(j), '-'
            Else
              Write (nout, 110) j, vr(1, j), rcondv(j), tol/rcondv(j)
            End If
          Else
            Write (nout, 100) j, vr(1, j), rcondv(j), 'Inf'
          End If

          Write (nout, 120) vr(2:n, j)

        End Do
        Write (nout, *)
        Write (nout, *) 'Errors below 10*machine precision are not displayed'
      Else
        Write (nout, *)
        Write (nout, 130) 'Failure in ZGEEVX. INFO =', info
      End If

100   Format (1X, I2, 1X, '(', 1P, E11.4, ',', E11.4, ')', 1X, 0P, F7.4, 4X, &
        A)
110   Format (1X, I2, 1X, '(', 1P, E11.4, ',', E11.4, ')', 1X, 0P, F7.4, 1X, &
        1P, E8.1)
120   Format (1X, 3X, '(', 1P, E11.4, ',', E11.4, ')')
130   Format (1X, A, I4)

    End Program


text:
      MODULE dex2_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:13  03/09/06
      REAL(KIND(0.0D0)) FUNCTION dex2 (M)
      integer, INTENT(IN) :: M
      END FUNCTION
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis

	function esilo(t)
c
c   this function returns the saturation vapor pressure over ice
c   esilo (millibars) given the temperature t (celsius).
c
c	baker,schlatter	17-may-1982	original version
c
c   the formula is due to lowe, paul r., 1977: an approximating polynomial for
c   the computation of saturation vapor pressure, journal of applied
c   meteorology, vol. 16, no. 1 (january), pp. 100-103.
c   the polynomial coefficients are a0 through a6.
c
	data a0,a1,a2,a3,a4,a5,a6
     1	/6.109177956,     5.034698970e-01, 1.886013408e-02,
     2	 4.176223716e-04, 5.824720280e-06, 4.838803174e-08,
     3	 1.838826904e-10/
	if (t.le.0.) go to 5
	esilo = 9999.
	write(6,3)esilo
    3	format(' saturation vapor pressure over ice is undefined for',
     1	/' temperature > 0c. esilo =',f6.0)
	return
    5	continue
	esilo = a0+t*(a1+t*(a2+t*(a3+t*(a4+t*(a5+a6*t)))))
	return
	end

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
program test_goto
21  GO TO 90
90  CONTINUE
end program

max_stars_repo_name: OpenFortranProject/ofp-sdf

max_issues_repo_name: OpenFortranProject/ofp-sdf

max_forks_repo_name: OpenFortranProject/ofp-sdf


text:
subroutine vdist3d(vx,vy,vz,lenv,nx,ny,nz,dfn)
   implicit none
   integer, intent(in) :: lenv, nx, ny, nz
   double precision, intent(in), dimension(lenv) :: vx, vy, vz
   double precision, intent(out), dimension(nx,ny,nz) :: dfn
   integer :: i,j,k,l
   double precision :: dvx, dvy, dvz
   double precision, dimension(6) :: mm

   mm=(/ maxval(vx) , minval(vx) &
        ,maxval(vy) , minval(vy) &
        ,maxval(vz) , minval(vz) /)

   dvx = (mm(1) - mm(2))/nx
   dvy = (mm(3) - mm(4))/ny
   dvz = (mm(5) - mm(6))/nz

   write(*,*) dvx,dvy,dvz

   do i=1,lenv
      j=floor((vx(i)-mm(2))/dvx + 0.5)
      k=floor((vy(i)-mm(4))/dvy + 0.5)
      l=floor((vz(i)-mm(6))/dvz + 0.5)
      !write(*,*) 'j,k,l',j,k,l
      dfn(j,k,l) = dfn(j,k,l) + 1
   enddo
end subroutine vdist3d

max_stars_repo_name: ahmadryan/TurbAn

max_issues_repo_name: ahmadryan/TurbAn

max_forks_repo_name: ahmadryan/TurbAn


text:
module q3_minima
  implicit none

  interface
    function pR2R(x)
      !! Real -> Real
      real, intent(in) :: x
      real :: pR2R
    end function pR2R
  end interface

contains

  subroutine q3()
    !! Find minima via golden section search
    print "(a)", "3. Minima of (x^2 - 1)^2 + x:"
    print "(f12.6)", golden(f3, x_min=-2.0, x_max=0.0, tol=1e-6)
    print "(f12.6)", golden(f3, x_min= 0.0, x_max=1.0, tol=1e-6)
    print *
  end subroutine

  function golden(f, x_min, x_max, tol)
    !! Find a minimum of f(x) within interval via golden section search
    !! This implementation doubles the necessary calls to f; but it is simpler
    procedure(pR2R) :: f
    real, intent(in) :: x_min, x_max, tol
    real :: golden, a, b, c, d
    real, parameter :: phi = 1.6180339887498948482045868343656381177203091798057

    a = x_min
    b = x_max
    c = b - (b - a) / phi
    d = a + (b - a) / phi

    do while (tol < abs(c - d))
      if (f(c) < f(d)) then
        b = d
      else
        a = c
      end if

      c = b - (b - a) / phi
      d = a + (b - a) / phi
    end do

    golden = 0.5 * (b + a)
  end function

  pure function f3(x)
    !! Function for question 3
    real, intent(in) :: x
    real :: f3

    f3 = (x**2 - 1.0)**2 + x
  end function

end module q3_minima

max_stars_repo_name: YodaEmbedding/experiments

max_issues_repo_name: YodaEmbedding/experiments

max_forks_repo_name: YodaEmbedding/experiments


text:
Describe Users/AdamD here.

20100708 09:35:46 nbsp Welcome to the Wiki. Dont use tiny url links on the wiki. They trigger spam alerts. Users/JasonAller

20100708 09:54:34 nbsp Spammers suck. Have you thought about what a terrible person you are, abusing community resources in an attempt to send traffic to your site and make money? Users/WilliamLewis

http://img818.imageshack.us/img818/6392/whywecanthavenicethings.jpg Users/hankim

20100709 08:22:14 nbsp We dont take kindly to your kind round here.  Users/MasonMurray

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
!##############################################################################
!# ****************************************************************************
!# <name> heatcond_method5 </name>
!# ****************************************************************************
!#
!# <purpose>
!# This module is a demonstation program how to solve a simple
!# heat conduction problem with constant coefficients
!# on a simple domain.
!#
!# The routine splits up the tasks of reading the domain, creating
!# triangulations, discretisation, solving, postprocessing and creanup into
!# different subroutines. The communication between these subroutines
!# is done using an application-specific structure saving problem data
!# as well as a collection structure for the communication with callback
!# routines.
!#
!# The different tasks of setting up the triangulation, discretisation, solver,
!# etc. are separated into different sub-modules as follows:
!#
!# - heatcond_basic
!#   -> Defines a basic problem structure that contains all parameters for
!#      the simulation
!#
!# - heatcond_callback
!#   -> Contains user-defined routines for setting up the RHS or for
!#      configuring values on the boundary.
!#
!# - heatcond_partridiscr
!#   -> Contains routines to read in the parametrisation, triangulation and
!#      to set up the discretisation (element, cubature rules,...=
!#
!# - heatcond_matvec
!#   -> Contains routines to initialise matrices and vectors.
!#
!# - heatcond_bounadrycondition
!#   -> Contains routines to set up analytical and discrete boundary conditions.
!#
!# - heatcond_solver
!#   -> Configures the linear solver for the problem.
!#
!# - heatcond_timeloop
!#   -> Contains the actual time stepping algorithm and nonstationary solver.
!# </purpose>
!##############################################################################

module heatcond_method5

  use fsystem
  use genoutput
  use storage
  use boundary
  use cubature
  use derivatives
  use matrixfilters
  use vectorfilters
  use discretebc
  use bcassembly
  use sortstrategy
  use triangulation
  use element
  use spatialdiscretisation
  use coarsegridcorrection
  use filtersupport
  use linearsystemscalar
  use linearsystemblock
  use scalarpde
  use bilinearformevaluation
  use linearformevaluation
  use multilevelprojection
  use linearsolver
  use discretebc
  use ucd

  use collection
  use paramlist

  use heatcond_callback

  use heatcond_basic
  use heatcond_matvec
  use heatcond_boundarycondition
  use heatcond_partridiscr
  use heatcond_solver
  use heatcond_timeloop

  implicit none

contains

  ! ***************************************************************************

!<subroutine>

  subroutine hc5_initparameters (rparams,rproblem)

!<description>
  ! Reads the DAT file from disc into the parameter list rparams and
  ! initialises basic variables (number of levels, time stepping technique)
  ! in rproblem according to these settings.
!</description>

!<inputoutput>
  ! A parameter list structure accepting the parameters from the DAT file.
  type(t_parlist), intent(inout) :: rparams

  ! A problem structure saving problem-dependent information.
  type(t_problem), intent(inout), target :: rproblem
!</inputoutput>

!</subroutine>

  ! local variables
  integer :: cscheme,niterations
  real(DP) :: dtheta,dtstep,dtimemin,dtimemax
  character(len=SYS_STRLEN) :: smaster,sstring

    ! Read the parameters from disc and put a reference to it
    ! to the collection
    call sys_getcommandLineArg(1,smaster,sdefault="./data/heatcond.dat")
    call parlst_readfromfile (rparams, smaster)

    ! We want to solve our Laplace problem on level...
    call parlst_getvalue_int (rparams, "GENERAL", "NLMIN", rproblem%ilvmin, 7)
    call parlst_getvalue_int (rparams, "GENERAL", "NLMAX", rproblem%ilvmax, 7)

    ! Get the parameters for the time stepping scheme from the parameter list
    call parlst_getvalue_int (rparams, "TIMESTEPPING", "CSCHEME", cscheme, 0)
    call parlst_getvalue_int (rparams, "TIMESTEPPING", "NITERATIONS", niterations, 1000)
    call parlst_getvalue_double (rparams, "TIMESTEPPING", "DTHETA", dtheta, 1.0_DP)
    call parlst_getvalue_double (rparams, "TIMESTEPPING", "DTSTEP", dtstep, 0.1_DP)
    call parlst_getvalue_double (rparams, "TIMESTEPPING", "DTIMEMIN", dtimemin, 0.0_DP)
    call parlst_getvalue_double (rparams, "TIMESTEPPING", "DTIMEMAX", dtimemax, 1.0_DP)

    ! Get the path where to write gmv`s to.
    call parlst_getvalue_string (rparams, "GENERAL", &
                                 "sucddir", sstring,"""./pre/QUAD.prm""")
    read(sstring,*) rproblem%sucddir

    ! Get the path of the prm/tri files.
    call parlst_getvalue_string (rparams, "GENERAL", &
                                 "sprmfile", sstring,"""./pre/QUAD.prm""")
    read(sstring,*) rproblem%sprmfile

    call parlst_getvalue_string (rparams, "GENERAL", &
                                 "strifile", sstring,"""./pre/QUAD.tri""")
    read(sstring,*) rproblem%strifile

    ! Initialise the time stepping in the problem structure
    call timstp_init (rproblem%rtimedependence%rtimestepping, &
                      cscheme, dtimemin, dtstep, dtheta)

    rproblem%rtimedependence%niterations = niterations

    rproblem%rtimedependence%dtimemin = dtimemin
    rproblem%rtimedependence%dtime = dtimemin
    rproblem%rtimedependence%dtimemax = dtimemax

  end subroutine

  ! ***************************************************************************

!<subroutine>

  subroutine heatcond5

!<description>
  ! This is a "separated" heatcond solver for solving a nonstationary heat
  ! conduction problem. The different tasks of the problem are separated into
  ! subroutines. The problem uses a problem-specific structure for the
  ! communication: All subroutines add their generated information to the
  ! structure, so that the other subroutines can work with them.
  ! (This is somehow a cleaner implementation than using a collection!).
  ! For the communication to callback routines of black-box subroutines
  ! (matrix-assembly), a collection is used.
  !
  ! The following tasks are performed by the subroutines:
  !
  ! 1.) Read in parametrisation
  ! 2.) Read in triangulation
  ! 3.) Set up RHS
  ! 4.) Set up matrix
  ! 5.) Create solver structure
  ! 6.) Solve the problem
  ! 7.) Write solution to VTK file
  ! 8.) Release all variables, finish
!</description>

!</subroutine>

    ! A paramlist structure with parameters from the dat file
    type(t_parlist) :: rparams

    ! A problem structure for our problem
    type(t_problem), target :: rproblem

    ! An initial RHS vector and a solution vector
    type(t_vectorBlock) :: rrhs,rvector

    integer :: i

    ! Initialise the parameter list
    call parlst_init(rparams)

    ! Initialise the collection.
    call collct_init (rproblem%rcollection)
    call collct_setvalue_parlst (rproblem%rcollection, "PARAMS", rparams, .true.)

    ! Read in the parameters from the DAT file and initialise the basic
    ! structures with these.
    call hc5_initparameters (rparams,rproblem)

    ! Add space for level information in the collection
    do i=1,rproblem%ilvmax
      call collct_addlevel (rproblem%rcollection)
    end do

    ! Allocate memory for the level information
    allocate (rproblem%RlevelInfo(rproblem%ilvmax))

    ! So now the different steps - one after the other.
    !
    ! Initialisation
    call hc5_initParamTriang (rproblem%ilvmin,rproblem%ilvmax,rproblem)
    call hc5_initDiscretisation (rproblem)
    call hc5_initMatVec (rproblem,rparams)

    ! Use the auxiliary RHS vector on the finest level to create an
    ! initial RHS and solution vector, which we pass later to the timeloop.
    call lsysbl_createVecBlockIndirect (rproblem%rrhs,rrhs,.false.)
    call lsysbl_createVecBlockIndirect (rproblem%rrhs,rvector,.true.)

    ! Calculate the initial RHS, do not incorporate any BC`s.
    call hc5_calcRHS (rproblem,rrhs)

    ! Discretise the boundary conditions
    call hc5_initDiscreteBC (rproblem)

    ! Implement them into the initial solution vector, as we have a zero
    ! vector as initial solution.
    rvector%p_rdiscreteBC => rproblem%rrhs%p_rdiscreteBC
    call vecfil_discreteBCsol (rvector)

    ! Initialise the solver
    call hc5_initSolver (rproblem)

    ! Call the timeloop to solve the problem
    call hc5_timeloop (rproblem,rvector,rrhs)

    ! Release the solver, we don't need it anymore
    call hc5_doneSolver (rproblem)

    ! Cleanup
    call hc5_doneMatVec (rproblem)
    call hc5_doneBC (rproblem)
    call hc5_doneDiscretisation (rproblem)
    call hc5_doneParamTriang (rproblem)

    ! Release memory for level information
    deallocate (rproblem%RlevelInfo)

    ! Release parameter list
    call collct_deletevalue (rproblem%rcollection,"PARAMS")
    call parlst_done (rparams)

    ! Release RHS and solution vector
    call lsysbl_releaseVector (rvector)
    call lsysbl_releaseVector (rrhs)

    ! Print some statistical data about the collection - anything forgotten?
    call output_lbrk ()
    call output_line ("Remaining collection statistics:")
    call output_line ("--------------------------------")
    call output_lbrk ()
    call collct_printStatistics (rproblem%rcollection)

    ! Finally release the collection.
    call collct_done (rproblem%rcollection)

  end subroutine

end module

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c        ACKNOWLEDGE THE USE OF THIS PACKAGE
c         WHEN YOU PUBLISH YOUR RESULTS !!!
c>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>
c
c************* Sample main program #10 *****************
c  1d Heisenberg antiferromagnet with 15 spins
c Degeneracy check by various initial vectors
c***********************************************
c
      parameter (n=15,idim=6435,ibond=n,ic=ibond+1)
      implicit real*8 (a-h,o-z)
      dimension E(4)
      dimension elemnt(idim,ic),loc(idim,ic)
      dimension list1(idim),list2(2,0:2**15)
      dimension bondwt(ibond),ipair(2*ibond),zrtio(ibond)
      dimension wk(idim,4)
      dimension x(idim)
      dimension v(idim,3),norm(3)
c
      data bondwt/ibond*-1.0d0/
      data zrtio/ibond*1.0d0/
      data ipair/1,2, 2,3,  3,4,  4,5,  5,6,  6,7,  7,8, 8,9,
     &     9,10, 10,11, 11,12, 12,13, 13,14, 14,15, 15,1/
      nvec=1
      iv=idim/5
c
      call sz(n,idim,0.5d0,list1,list2)
c- You may alternatively use szdy or sztn for faster processing -
c      call szdy(n,idim,0.5d0,list1,list2)
c  or
c      call sztn(n,idim,0.5d0,list1,list2)
c------------------------------------------------------------
      call elm2(n,idim,ipair,bondwt,zrtio,ibond,
     &          elemnt,loc,idim,ic,list1,list2)
      k=0
c *** Two different initial conditions
      do 10 iv=91,idim,idim/3
        k=k+1
        call lnc2(elemnt,loc,idim,idim,ic,nvec,iv,E,itr,wk)
        print 100,k,e
  100    format(/' #',i2,' [Eigenvalues]  '/2x,4f14.8)
        call inv2(elemnt,loc,idim,idim,ic,iv,E(1),x,wk)
        print *,'[Eigenvector components (selected)]'
        print 120,(x(j),j=13,idim,idim/20)
  120    format(4d18.9)
        call check2(elemnt,loc,idim,idim,ic,x,wk,Hexpec)
        do 20 j=1,idim
  20     v(j,k)=x(j)
  10   continue
c *** Degeneracy check
      call orthg(idim,idim,v,norm,idgn,3)
      print 110,idgn
  110  format(/' [Degeneracy]'/i6)
      end

max_stars_repo_name: wistaria/rokko

max_issues_repo_name: wistaria/rokko

max_forks_repo_name: wistaria/rokko


text:
!--------------------------------------------------------------------------------
! Copyright (c) 2016 Peter Grünberg Institut, Forschungszentrum Jülich, Germany
! This file is part of FLEUR and available as free software under the conditions
! of the MIT license as expressed in the LICENSE file in more detail.
!--------------------------------------------------------------------------------

      module m_wann_rad_twf
      use m_juDFT
      contains
      subroutine wann_rad_twf(
     >               nwfs,jmtd,natd,ind,rwf,zona,regio,
     >               us,dus,uds,duds,ff,gg,lmaxd,ikpt,
     >               ntypd,ntp,jri,rmsh,dx,
     >               nlod,flo,llo,nlo,
     <               rads)
c*********************************************************************
c..this subroutine calculates on the radial grid the radial function
c..corresponding to the twf, its mt sphere parameters, diffusivity
c..and region
c
c..Y. Mokrousov
c
c..tidied-up && extension of rwf-parameter
c..Frank Freimuth
c*********************************************************************

      USE m_constants
      use m_intgr, only : intgr3

      implicit none

      integer, intent (in)  :: nwfs,natd,jmtd,ntypd,lmaxd,ikpt,nlod
      integer, intent (in)  :: rwf(nwfs),ind(nwfs),ntp(natd),jri(ntypd)
      real,    intent (in)  :: zona(nwfs),regio(nwfs),rmsh(jmtd,ntypd)
      real,    intent (in)  :: us(0:lmaxd,ntypd),dus(0:lmaxd,ntypd)
      real,    intent (in)  :: uds(0:lmaxd,ntypd),duds(0:lmaxd,ntypd)
      real,    intent (in)  :: ff(ntypd,jmtd,2,0:lmaxd)
      real,    intent (in)  :: gg(ntypd,jmtd,2,0:lmaxd),dx(ntypd)
      real,    intent (in)  :: flo(ntypd,jmtd,2,nlod)
      integer, intent (in)  :: llo(nlod,ntypd)
      integer, intent (in)  :: nlo(ntypd)
      real,    intent (out) :: rads(nwfs,0:3,jmtd,2)

      integer :: nwf,nat,ntyp,j,n,l,lo
      real    :: rr,alpha,const,fact,wronk,radi
      real    :: acft,bcft,aa,bb,a1,b1,rho,radf(jmtd)

      do nwf = 1,nwfs

       if (abs(regio(nwf)-1.00000).ge.0.00001) then
          CALL juDFT_error("no projections outside muffin tins",calledby
     +         ="wann_rad_twf")
       endif

       nat = ind(nwf)
       ntyp = ntp(nat)
       alpha = zona(nwf)

       rads (nwf,:,:,:) = 0.

       if (rwf(nwf).eq.-5) then
c*********************************************************************
c..For rwf.eq.-5 the following way of implementing is realized.
c..Suppose, the twf has its center on some atom. The
c..self-consistent potential V_l(r) and the energy parameters
c..E_l are known. A linear combination of the u_l and \dot{u}_l
c..is constructed, with the coefficients A and B such that
c..on the mt boundary A*u_l + B*\dot{u}_l is smoothly continuous
c..with the 1s function 2*((zona)**(1.5))*exp(-r*zona)
c*********************************************************************
         aa = 2.*(zona(nwf)**(1.5))
         bb = zona(nwf)

         a1 = aa*exp(-bb*rmsh(jri(ntyp),ntyp))
         b1 = -aa*bb*exp(-bb*rmsh(jri(ntyp),ntyp))

         do l = 0,3

            wronk = us(l,ntyp)*duds(l,ntyp)-uds(l,ntyp)*dus(l,ntyp)

            acft = (a1*duds(l,ntyp) - b1*uds(l,ntyp))/wronk
            bcft = (b1*us(l,ntyp) - a1*dus(l,ntyp))/wronk

            do j = 1,jri(ntyp)
               rads(nwf,l,j,:) = acft*ff(ntyp,j,:,l) +
     +                           bcft*gg(ntyp,j,:,l)
            enddo

         enddo
       elseif(rwf(nwf).eq.0) then
c***********************************************************
c..For rwf.eq.0 the zona parameter mixes the radial function
c..and its energy derivative linearly.
c***********************************************************
         do l = 0,3
            do j = 1,jri(ntyp)
               rads(nwf,l,j,:) = ff(ntyp,j,:,l)*(1.-abs(zona(nwf))) +
     +                         zona(nwf)*gg(ntyp,j,:,l)
            enddo
         enddo
c**************************************
c..project onto local orbitals
c**************************************
       elseif(rwf(nwf).eq.-6)then
         IF(nlod<1) CALL juDFT_error("nlod<1",calledby="wann_rad_twf"
     +         )
          do l=0,3
           do j=1,jri(ntyp)
            rads(nwf,l,j,:)=flo(ntyp,j,:,1)
           enddo!j
          enddo!l
       elseif(rwf(nwf).eq.-7)then
         IF(nlod<2) CALL juDFT_error("nlod<2",calledby="wann_rad_twf"
     +         )
          do l=0,3
           do j=1,jri(ntyp)
            rads(nwf,l,j,:)=flo(ntyp,j,:,2)
           enddo!j
          enddo!l
       elseif(rwf(nwf).eq.-8)then
          IF(nlod<3) CALL juDFT_error("nlod<3",calledby ="wann_rad_twf"
     +         )
          do l=0,3
           do j=1,jri(ntyp)
            rads(nwf,l,j,:)=flo(ntyp,j,:,3)
           enddo!j
          enddo!l
       elseif(rwf(nwf).eq.-9)then
          rads(nwf,:,:,:)=0
          do lo=1,nlo(ntyp)
           l=llo(lo,ntyp)
           do j=1,jri(ntyp)
            rads(nwf,l,j,:)=flo(ntyp,j,:,lo)
           enddo!j
          enddo !lo

       elseif((-5.lt.rwf(nwf)).and.(rwf(nwf).lt.0)) then
c************************************************************
c..use the radial with l=abs(rwf)-1 for all components of the
c..trial wave function
c************************************************************
           do l = 0,3
            do j = 1,jri(ntyp)
               rads(nwf,l,j,:) =
     =             ff(ntyp,j,:,abs(rwf(nwf))-1)*(1.-abs(zona(nwf))) +
     +             zona(nwf)*gg(ntyp,j,:,abs(rwf(nwf))-1)
            enddo
           enddo

c**************************************************************
c..Hydrogen-like test orbitals.
c..here we have tabulated the radial functions from orbitron.
c..for 1 <= rwf <= 6 different components of the test orbital
c..are constructed from different radials.
c  For Hybrides (e.g. sp3) it might be better to use the same
c  radials for all components. See below (7 <= rwf <= 12)
c**************************************************************
       elseif (rwf(nwf).eq.1) then
c..n=1
          do j = 1,jri(ntyp)
             rho = 2.*alpha*rmsh(j,ntyp)
             rads(nwf,0,j,1) = 2.*((alpha)**(1.5))*
     *            exp(-rho/2.)
          enddo

       elseif (rwf(nwf).eq.2) then
c..n=2
          do j = 1,jri(ntyp)
             rho = alpha*(rmsh(j,ntyp))
             rads(nwf,0,j,1) = (1./(2.*sqrt(2.)))*
     *              (2.-rho)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,1,j,1) = (1./(2.*sqrt(6.)))*
     *              rho*
     *              ((alpha)**(1.5))*exp(-rho/2.)
          enddo

       elseif (rwf(nwf).eq.3) then
c..n=3
          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/3.
             rads(nwf,0,j,1) = (1./(9.*sqrt(3.)))*
     *              (6. - 6.*rho + rho**2)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,1,j,1) = (1./(9.*sqrt(6.)))*
     *              (4. - rho)*rho*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,2,j,1) = (1./(9.*sqrt(30.)))*
     *               rho*rho*
     *              ((alpha)**(1.5))*exp(-rho/2.)
          enddo

       elseif (rwf(nwf).eq.4) then
c..n=4
          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/4.
             rads(nwf,0,j,1) = (1./96.)*
     *              (24. - 36.*rho + 12.*(rho**2) - (rho**3))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,1,j,1) = (1./(32.*sqrt(15.)))*
     *              rho*(20. - 10.*rho + rho**2)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,2,j,1) = (1./(96.*sqrt(5.)))*
     *              rho*rho*(6. - rho)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,3,j,1) = (1./(96.*sqrt(35.)))*
     *              rho*rho*rho*
     *              ((alpha)**(1.5))*exp(-rho/2.)
          enddo

       elseif (rwf(nwf).eq.5) then
          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/5.
             rads(nwf,0,j,1) = (1./(300.*sqrt(5.)))*
     *              (120. - 240.*rho + 120.*(rho**2) - 20.*(rho**3)
     +                            +(rho**4))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,1,j,1) = (1./(150.*sqrt(30.)))*
     *              rho*(120. - 90.*rho + 18.*(rho**2) - (rho**3))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,2,j,1) = (1./(150.*sqrt(70.)))*
     *              rho*rho*(42. - 14.*rho + (rho*rho))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,3,j,1) = (1./(300.*sqrt(70.)))*
     *              rho*rho*rho*(8. - rho)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
          enddo
        elseif (rwf(nwf).eq.6) then

          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/6.
             rads(nwf,0,j,1) = (1./2160.*(sqrt(6.)))*
     *              (720. - 1800.*rho + 1200.*rho*rho-
     -               300.*rho*rho*rho + 30.*(rho**4) - (rho**5))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             rads(nwf,1,j,1) = (1./432.*(sqrt(210.)))*
     *              rho*(840. - 840.*rho + 252.*rho*rho-
     -               28.*rho*rho*rho + (rho**4))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
          enddo

        elseif (rwf(nwf).eq.7) then
c..n=1
          do j = 1,jri(ntyp)
             rho = 2.*alpha*rmsh(j,ntyp)
             do l=0,3
             rads(nwf,l,j,1) = 2.*((alpha)**(1.5))*
     *            exp(-rho/2.)
             enddo
          enddo

        elseif (rwf(nwf).eq.8) then
c..n=2
          do j = 1,jri(ntyp)
             rho = alpha*(rmsh(j,ntyp))
             do l=0,3
             rads(nwf,l,j,1) = (1./(2.*sqrt(2.)))*
     *              (2.-rho)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             enddo
          enddo

        elseif (rwf(nwf).eq.9) then
c..n=3
          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/3.
             do l=0,3
             rads(nwf,l,j,1) = (1./(9.*sqrt(3.)))*
     *              (6. - 6.*rho + rho**2)*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             enddo
          enddo

        elseif (rwf(nwf).eq.10) then
c..n=4
          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/4.
             do l=0,3
             rads(nwf,l,j,1) = (1./96.)*
     *              (24. - 36.*rho + 12.*(rho**2) - (rho**3))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             enddo
          enddo
        elseif (rwf(nwf).eq.11) then

          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/5.
             do l=0,3
             rads(nwf,l,j,1) = (1./(300.*sqrt(5.)))*
     *              (120. - 240.*rho + 120.*(rho**2) - 20.*(rho**3)
     +                            +(rho**4))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             enddo
          enddo

        elseif (rwf(nwf).eq.12) then

          do j = 1,jri(ntyp)
             rho = 2*alpha*rmsh(j,ntyp)/6.
             do l=0,3
             rads(nwf,0,j,1) = 0*(1./2160.*(sqrt(6.)))*
     *              (720. - 1800.*rho + 1200.*rho*rho-
     -               300.*rho*rho*rho + 30.*(rho**4) - (rho**5))*
     *              ((alpha)**(1.5))*exp(-rho/2.)
             enddo
          enddo

        else

           CALL juDFT_error("radial function is not tabulated" ,calledby
     +          ="wann_rad_twf")

       endif

       if (ikpt.eq.1) then
        do l = 0,3
         do j = 1,jri(ntyp)
          radf(j) = rmsh(j,ntyp)*rmsh(j,ntyp)*rads(nwf,l,j,1)*
     *                rads(nwf,l,j,1)
c         radf(j) = rads(nwf,l,j,1)*rads(nwf,l,j,1)
         enddo
         call intgr3(radf,rmsh(1,ntyp),dx(ntyp),jri(ntyp),radi)
         write (oUnit,*)
         write (oUnit,*) 'Wannier Function N:',nwf
         write (oUnit,*) 'angular momentum',l
         write (oUnit,*) 'radial function at the MT boundary:',
     &                rads(nwf,l,jri(ntyp),1)
         write (oUnit,*) 'norma =',radi
        enddo
       endif

      enddo ! by twfs

      return

      end subroutine wann_rad_twf
      end module m_wann_rad_twf

max_stars_repo_name: MRedies/FLEUR

max_issues_repo_name: MRedies/FLEUR

max_forks_repo_name: MRedies/FLEUR


text:
	SUBROUTINE RS_RTIM  ( isdtar, day, hour, data, defalt, obfreq,
     +			      irdtar, rtime, ihhmm, iret )
C************************************************************************
C* RS_RTIM								*
C*									*
C* This subroutine combines an integer system time, the bulletin	*
C* time containing the day and hour and the report day, hour, and	*
C* minute into an observation time.  It is assumed that the time in	*
C* isdtar accurately reflects the year and month of the observation	*
C* and is later than that time.						*
C*									*
C* RS_RTIM  ( ISDTAR, DAY, HOUR, DATA, DEFALT, OBFREQ, IRDTAR, RTIME,	*
C*	      IHHMM, IRET )						*
C*									*
C* Input parameters:							*
C*	ISDTAR (5)	INTEGER		System time			*
C*	DAY		INTEGER		Nominal day of report		*
C*	HOUR		INTEGER		Nominal hour of report		*
C*	DATA (*)	INTEGER		Decoded data from RS_DECO	*
C*	DEFALT		INTEGER		Missing data flag		*
C*	OBFREQ		INTEGER		Interval (min) between ob times	*
C*									*
C* Output parameters:							*
C*	IRDTAR (5)	INTEGER		Integer report nominal time	*
C*	RTIME		CHAR*		Report nominal date/time	*
C*	IHHMM		INTEGER		Report actual time		*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -1 = invalid time		*
C*					 -2 = time conflicts with obfreq*
C**									*
C* Log:									*
C* J. Nielsen		 2/92	Loosely dapted from RA_RTIM		*
C* D. Kidwell/NCEP	 9/98	Fixed bugs in month, year changeovers   *
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	INTEGER		isdtar (*), day, hour, data (*), obfreq, defalt
	INTEGER		irdtar (*), ihhmm
	CHARACTER*(*)	rtime
C*
	LOGICAL		subday, addday
C*
	INTEGER		iodtar (5)
	EQUIVALENCE	( iodtar (1), ioyear ), ( iodtar (2), iomnth ),
     +			( iodtar (3), ioday  ), ( iodtar (4), iohour ),
     +			( iodtar (5), iomin  )
C------------------------------------------------------------------------
	iret   = 0
	addday = .false.
	subday = .false.
C
C*	Use the system time for the observation time for the year and
C*	month.
C
	ioyear = isdtar (1)
	iomnth = isdtar (2)
C
C*	Get times from report, from header, or from system time, in
C*	order of preference.  Fix date if we have an ob time of 00Z
C*	and a nominal time of prior to 00Z.  Use minutes from report
C*	if given, and assume zero otherwise.
C
	IF  (( data (91) .ne. defalt ) .and.
     +	     ( data (91) .le. 2359 ))  THEN
C
C*	    Hour and minute of ob reported
C
	    iomin = MOD ( data(91), 100 )
	    iohour = ( data(91) - iomin) / 100
	    IF  ( iomin .ge. 60 )  iomin = 0
	    IF  ( iomin .lt. 0 )  iomin = 0
	    IF  (( data (3) .ne. defalt ) .and.
     +		 ( data (3) .le. 31 ))  THEN
	        ioday = data (3)
		IF  ( iohour - data (4) .eq. 23 )  subday = .true.
	    ELSE IF  ( day .le. 31 )  THEN
	        ioday = day
		IF  ( iohour - hour .eq. 23 )  subday = .true.
	    ELSE
	        ioday = isdtar (3)
	    ENDIF
	ELSE IF  (( data(4) .ne. defalt ) .and.
     +		  ( data(4) .le. 23 ))  THEN
C
C*	    Day and hour of ob reported
C
	    iohour = data(4)
	    iomin = 0
	    IF  (( data (3) .ne. defalt ) .and.
     +		 ( data (3) .le. 31 ))  THEN
	        ioday = data (3)
	    ELSE IF  ( day .le. 31 )  THEN
	        ioday = day
		IF  ( iohour - hour .eq. 23 )  subday = .true.
	    ELSE
	        ioday = isdtar (3)
	    ENDIF
	ELSE IF  ( hour .le. 23 ) THEN
C
C*	    Get day and hour from header info
C
	    iohour = hour
	    iomin = 0
	    IF  ( day .le. 31 )  THEN
	        ioday = day
	    ELSE
	        ioday = isdtar (3)
	    ENDIF
	ELSE
C
C*	    Use system day and hour
C
	    iohour = isdtar (4)
	    iomin = 0
	    ioday = isdtar (3)
	ENDIF
C
C*	Make sure the report and bulletin time are within one hour.
C
	ihdiff = ABS  ( hour - iohour )
	IF  ( ( ihdiff .gt. 1 ) .and. ( ihdiff .ne. 23 ) ) THEN
C
C*	    Bulletin and report time do not match.
C
	    iret  = -1
	    rtime = ' '
	    RETURN
	END IF
C
C*	Before adjusting time to correspond to obfreq, save ihhmm
C
	ihhmm = iohour * 100 + iomin
C
C*	Establish the proper minutes using obfreq
C
	itmin = iohour * 60 + iomin
	IF  ( obfreq .gt. 60 )  THEN
	    ifreq = 60
	ELSE
	    ifreq = obfreq
	END IF
	IF  ( obfreq .gt. 0 )  THEN
	    ileft = MOD ( itmin, obfreq )
	    IF  ( ileft .ge. ifreq/2 )  THEN
		IF  ( ( obfreq - ileft ) .gt. ifreq/2 )  THEN
		    iret = -2
		    RETURN
		END IF
		itmin = itmin - ileft + obfreq
	    ELSE
		itmin = itmin - ileft
	    ENDIF
	ENDIF
	IF  ( itmin .ge. 24*60 )  THEN
	    addday = .true.
	    itmin = itmin - 24*60
	ENDIF
	iomin = MOD ( itmin, 60 )
	iohour = ( itmin - iomin ) / 60

C
C*	Check that the month has not changed.
C
	IF  ( ( isdtar (3) .eq. 1 ) .and. ( iodtar (3) .gt. 27 ) )
     +								THEN
	    iomnth = iomnth - 1
	    IF  ( iomnth .eq. 0 )  THEN
		iomnth = 12
		ioyear = ioyear - 1
	    END IF
	ELSE IF  ( ( isdtar (3) .gt. 27 ) .and.
     +		   ( iodtar (3) .eq. 1 ) )  THEN
	    iomnth = iomnth + 1
	    IF  ( iomnth .eq. 13 )  THEN
		iomnth = 1
		ioyear = ioyear + 1
	    END IF
	END IF
C
C*	Add or subtract days if necessary.
C
	IF  ( subday )  CALL TI_SUBD  ( iodtar, iodtar, ier )
	IF  ( addday )  CALL TI_ADDD  ( iodtar, iodtar, ier )
C
C*	Convert to a GEMPAK time.
C
	CALL TI_ITOC  ( iodtar, rtime, iret )
	DO  i = 1, 5
	    irdtar (i) = iodtar (i)
	END DO
C*
	RETURN
	END


text:
This editor can edit this entry and tell us a bit about themselves by clicking the Edit icon.

20110512 18:07:30 nbsp Welcome to the Wiki. How about introducing yourself and expanding your edits to share your knowledge of Davis beyond the pages for Cascade Apartments, Viking Apartments, Sherwood & Forest Arms Apartments, and StarCrossed Properties? Users/JasonAller


text:
cdis
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
cdis
cdis
cdis

        subroutine vert_wind(uanl,vanl                              ! I ! /O
     1          ,u_sfc,v_sfc,ni,nj,nk                               ! I
     1          ,wanl                                               ! O
     1          ,topo,lat,lon,grid_spacing_m                        ! I
     1          ,rk_terrain,r_missing_data,l_grid_north             ! I
     1          ,istatus)                                           ! O

!      ~1990        Steve Albers  Orig version
!       1997 Jun    Ken Dritz     Added ni and nj as dummy
!                                 arguments, thus making wsum and one and
!                                 several other arrays automatic.
!       1997 Jun    Ken Dritz     Initialize wsum and one dynamically, instead
!                                 of by DATA statements.
!       1997 Jun    Ken Dritz     Added r_missing_data as dummy argument.
!       1997 Jun    Ken Dritz     Removed include of 'lapsparms.for'.
!       1997 Oct    Steve Albers  Pass lon to fflxc. Misc Cleanup.
!       1997 Dec    Steve Albers  Changed NX_L_MAX/NY_L_MAX to ni/nj
!       1998 Nov    Steve Albers  Change to M, hopefully cancels out change to
!                                 fflxc/sigma

        include 'trigd.inc'
        real m ! Grid points per meter

        real wsum(ni,nj)
        real  one(ni,nj)

        logical l_grid_north ! Sfc & 3D winds

        real rk_terrain(ni,nj)
        integer k_terrain(ni,nj)

        DATA scale/1./

        real uanl(ni,nj,nk),vanl(ni,nj,nk)
        real wanl(ni,nj,nk) ! omega (pascals/second)
        real terrain_w(ni,nj),conv(ni,nj)
        real u_sfc(ni,nj),v_sfc(ni,nj)
        real topo_pa(ni,nj)

        real lat(ni,nj)
        real lon(ni,nj)
        real topo(ni,nj)

        real flu(ni,nj)
        real flv(ni,nj)
        real sigma(ni,nj)

        real u_sfc_grid(ni,nj),
     1         v_sfc_grid(ni,nj)
        real beta_factor(ni,nj)

        real radius_earth
        parameter (radius_earth = 6371e3)

        character*6 c6_maproj

        do i=1,ni
           do j=1,nj
              wsum(i,j) = 0.0
              one(i,j) = 1.0
           enddo
        enddo

        ierrcnt = 0
        istatus = 1

!       grid_spacing_m = sqrt(
!       1                      (  lat(1,2) - lat(1,1)                  )**2
!       1                    + ( (lon(1,2) - lon(1,1))*cosd(lat(1,1))  )**2
!       1                                       )    * 111317. ! Grid spacing m

        call get_standard_latitudes(slat1,slat2,istatus)
        if(istatus .ne. 1)then
            return
        endif

        call get_c6_maproj(c6_maproj,istatus)
        if(istatus .ne. 1)then
            return
        endif

        if(c6_maproj .eq. 'plrstr')then
            call get_ps_parms(slat1,slat2,grid_spacing_m,phi0
     1                                 ,grid_spacing_proj_m)
        else
            grid_spacing_proj_m = grid_spacing_m
        endif

        write(6,*)' Grid spacings (m) = ',grid_spacing_m
     1                                   ,grid_spacing_proj_m

        m = 1.0 / grid_spacing_proj_m

        imaxm1 = ni - 1
        jmaxm1 = nj - 1

!       Calculate surface pressure array
        do j=1,nj
        do i=1,ni
            topo_pa(i,j) = ztopsa(topo(i,j)) * 100.
        enddo
        enddo

!       Rotate Sfc Winds
        if(.not. l_grid_north)then
          write(6,*)' Rotating SFC winds to calculate terrain forcing'
          write(6,*)' Calculating Beta Factor'
          do j=1,nj
          do i=1,ni
            if(u_sfc(i,j) .ne. r_missing_data
     1                   .and. v_sfc(i,j) .ne. r_missing_data)then

                call uvtrue_to_uvgrid(u_sfc(i,j),
     1                              v_sfc(i,j),
     1                              u_sfc_grid(i,j),
     1                              v_sfc_grid(i,j),
     1                              lon(i,j))

            endif

            beta_factor(i,j) = sind(lat(i,j)) / radius_earth

          enddo
          enddo

        else ! Winds are already rotated
          do j=1,nj
          do i=1,ni
              u_sfc_grid(i,j) = u_sfc(i,j)
              v_sfc_grid(i,j) = v_sfc(i,j)
              beta_factor(i,j) = 0.
          enddo
          enddo

        endif

        if(c6_maproj .eq. 'latlon')then
            write(6,*)' latlon grid, skipping omega...'
            wanl = r_missing_data
            return
        endif

!       Calculate terrain induced omega
        do j=2,jmaxm1
        do i=2,imaxm1

!           Units of dterdx are pascals/meter
            dterdx=(topo_pa(i+1,j) - topo_pa(i-1,j)) * .5/grid_spacing_m
            dterdy=(topo_pa(i,j+1) - topo_pa(i,j-1)) * .5/grid_spacing_m

!           Units of ubar,vbar are m/s
            ubar=(u_sfc_grid(i-1,j-1) + u_sfc_grid(i,j-1)) *.5
            vbar=(v_sfc_grid(i-1,j-1) + v_sfc_grid(i-1,j)) *.5

!           Units of terrain_w are pascals/second
            terrain_w(i,j) = UBAR*DTERDX+VBAR*DTERDY

        enddo ! i
        enddo ! j

!       Fill in edges
        do i = 1,ni
            terrain_w(i,   1) = terrain_w(i,     2)
            terrain_w(i,nj)   = terrain_w(i,jmaxm1)
        enddo ! i

        do j = 1,nj
            terrain_w(1,   j) = terrain_w(2,     j)
            terrain_w(ni,j)   = terrain_w(imaxm1,j)
        enddo ! j

        do j = 1,nj
        do i = 1,ni
            k_terrain(i,j) = max( nint(rk_terrain(i,j)) ,1 )
        enddo ! i
        enddo ! j

        write(6,*)'            terrain       conv       omega      beta
     1corr'

        do k = 1,nk
            call FFLXC(ni,nj,M,SCALE
     1                ,uanl(1,1,k),vanl(1,1,k),one,conv,lat,lon
     1                ,flu,flv,sigma,r_missing_data)

            if(k.gt.1)then
                z_interval=abs(zcoord_of_level(k)-zcoord_of_level(k-1))
            else
                z_interval=abs(zcoord_of_level(2)-zcoord_of_level(2-1))
            endif

            do j = 1,nj
            do i = 1,ni

                k_terr = k_terrain(i,j)

                if(k .lt. k_terr)then !
                    wanl(i,j,k) = r_missing_data
!                   uanl(i,j,k) = 0.
!                   vanl(i,j,k) = 0.

                elseif(k .eq. k_terr)then !
                    wsum(i,j)   = terrain_w(i,j)
!    1       - (conv(i,j) + vanl(i,j,k) * beta_factor(i,j)) * z_interval
                    wanl(i,j,k) = wsum(i,j)

                else ! k .gt. k_terr
                    wsum(i,j)   = wsum(i,j)
     1       - (conv(i,j) + vanl(i,j,k) * beta_factor(i,j)) * z_interval
                    wanl(i,j,k) = wsum(i,j)

                endif

                if(abs(conv(i,j)) .gt. 1.0)then
                    ierrcnt = ierrcnt + 1
                    if(ierrcnt .lt. 20)then
                        write(6,*)' Error: Large Convergence'
     1                          ,i,j,k,k_terrain(i,j),conv(i,j)
                    endif

                    istatus = 0
                endif

                if(j .eq. 29 .and. i .eq. 29)then
                    write(6,111)i,j,k,terrain_w(i,j)
     1              ,conv(i,j),wanl(i,j,k),beta_factor(i,j)*vanl(i,j,k)
111                 format(3i3,4e12.3)
                endif

            enddo ! i
            enddo ! j

        enddo ! k

        return
        end


text:
!-----------------------------------------
subroutine vortex_force_calc(U_k, U_j, vort_k, vort_j, res, order, l, m, n)

  implicit none

  integer, intent(in):: l, m, n, order
  real*8, dimension(l, m, n), intent(in) :: U_k, U_j , vort_k, vort_j
  real*8, dimension(l, m, n), intent(inout) :: res

  !f2py intent(inplace):: vort_k, vort_j, res, U_k, U_j

  integer:: i, j, k
  integer, parameter :: out_unit=30
  real*8:: c1,c2,c3, b1, b2, b3, b4, b5
  real*8:: UU_k, up_k, um_k, qm_k, qp_k, U_k_interp, UU_j, up_j, um_j, qm_j, qp_j, U_j_interp, U_omega

  c1 = -1./6.
  c2 =  5./6.
  c3 =  2./6.

  b1 = 2./60.
  b2 = -13./60.
  b3 = 47./60.
  b4 = 27./60.
  b5 = -3./60.

  do k = 1,l
    do j = 1,m
      do i = 1,n

        !w_j U_k
        if ((k.gt.1).and.(j.lt.m)) then
          U_k_interp = 0.25 * ( U_k(k, j, i) + U_k(k-1, j, i) + U_k(k-1, j+1, i) + U_k(k, j+1, i) )
        elseif (j.lt.m) then
          U_k_interp = 0.25 * ( U_k(k, j+1, i) + U_k(k, j, i)) !no normal velocity
        else
        endif

        UU_k = abs(U_k_interp)
        up_k = 0.5*(U_k_interp+UU_k) ! right-going flux
        um_k = 0.5*(U_k_interp-UU_k) ! left-going flux

        !write(out_unit,*) "up_k ", up_k, "um_k", um_k

        if ((i.gt.2).and.(order.eq.5)) then
           qp_j = b1*vort_j(k,j,i-2) + b2*vort_j(k,j,i-1) + b3*vort_j(k,j,i) + b4*vort_j(k,j,i+1) + b5*vort_j(k,j,i+2)
        elseif ((i.gt.1).and.(order.ge.3)) then
           qp_j = c1*vort_j(k,j,i-1) + c2*vort_j(k,j,i) + c3*vort_j(k,j,i+1)
        else
           qp_j = vort_j(k,j,i)
        endif

        if ((i.le.(n-3)).and.(order.eq.5)) then
           qm_j = b5*vort_j(k,j,i-1) + b4*vort_j(k,j,i) + b3*vort_j(k,j,i+1) + b2*vort_j(k,j,i+2) + b1*vort_j(k,j,i+3)
        elseif ((i.le.(n-2)).and.(order.ge.3)) then
           qm_j = c3*vort_j(k,j,i) + c2*vort_j(k,j,i+1) + c1*vort_j(k,j,i+2)
        elseif (i.le.(n-1)) then
           qm_j = vort_j(k,j,i+1)
        else
           ! no flux through the boundary
           ! TODO: change this in case of inflow/outflow
           ! cf Winters 2012
           qm_j = 0.
        endif

        !w_k U_j
        if ((k.gt.1).and.(i.lt.m)) then
          U_j_interp = 0.25 * ( U_j(k, j, i) + U_j(k-1, j, i) + U_j(k-1, j, i+1) + U_j(k, j, i+1) )
        elseif (k.lt.l) then
          U_j_interp = 0.25 * ( U_j(k+1, j, i) + U_j(k, j, i)) !no normal velocity
        else
        endif

        UU_j = abs(U_j_interp)
        up_j = 0.5*(U_j_interp+UU_j) ! right-going flux
        um_j = 0.5*(U_j_interp-UU_j) ! left-going flux

        !write(out_unit,*) "up_j ", up_j, "um_j", um_j

        if ((j.gt.2).and.(order.eq.5)) then
           qp_k = b1*vort_k(k,j-2,i) + b2*vort_k(k,j-1,i) + b3*vort_k(k,j,i) + b4*vort_k(k,j+1,i) + b5*vort_k(k,j+2,i)
        elseif ((j.gt.1).and.(order.ge.3)) then
           qp_k = c1*vort_k(k,j-1,i) + c2*vort_k(k,j,i) + c3*vort_k(k,j+1,i)
        else
           qp_k = vort_k(k,j,i)
        endif

        if ((j.le.(m-3)).and.(order.eq.5)) then
           qm_k = b5*vort_k(k,j-1,i) + b4*vort_k(k,j,i) + b3*vort_k(k,j+1,i) + b2*vort_k(k,j+2,i) + b1*vort_k(k,j+3,i)
        elseif ((j.le.(m-2)).and.(order.ge.3)) then
           qm_k = c3*vort_k(k,j,i) + c2*vort_k(k,j+1,i) + c1*vort_k(k,j+2,i)
        elseif (j.lt.m) then
           qm_k = vort_k(k,j+1,i)
        else
           ! no flux through the boundary
           ! TODO: change this in case of inflow/outflow
           ! cf Winters 2012
           qm_k = 0.
        endif

        !write(out_unit,*), "qm_j", qm_j, "qp_j", qp_j, "qp_k", qp_k, "qm_k", qm_k

        ! this term is U_j * omega_k

        U_omega =  (up_k*qp_j + um_k*qm_j) - ( up_j * qp_k + um_j * qm_k)

        !write(out_unit,*) "U_omega ", U_omega

        res(k,j,i) = res(k,j,i) - U_omega

      end do
    end do
  end do

end subroutine vortex_force_calc

max_stars_repo_name: pvthinker/Nyles

max_issues_repo_name: pvthinker/Nyles

max_forks_repo_name: pvthinker/Nyles