text:
      SUBROUTINE IFP1B
C
C     THIS ROUTINE DETERMINES THE LOOP CONDITIONS AND CASE CONTROL
C     REQUEST CHNGES
C     LOOP$ -- THE CURRENT PROBLEM WILL LOOP
C
C     LOOP1$-- THE OLD PROBLEM WAS A LOOP AND CASE CONTROL IS CHANGED
C     IN LENGTH
C
C     COMMENTS FROM G.C.  10/92
C     IWORD AND IBIT 200 WORDS EACH CORRESPOND TO 200 WORDS IN CASECC
C     ZERO IN IWORD MEANS NO FURTHER CHECKING
C     INTEGER VALUE IN IBIT POINTS TO RESTART BIT POSITION, AND WILL BE
C     SAVED IN BITS(17) AND BITS(18), BITS FOR LCC. (LBD = 16)
C
C     LAST REVISED  7/91, BY G.CHAN/UNISYS, TO ALLOW HUGE THRU-RANGE ON
C     SET IN CASE CONTROL SECTION FOR PRINTOUT OR PLOTTING
C
      EXTERNAL        ANDF,ORF
      LOGICAL         NEW,DEBUG
      INTEGER         NAME(2),OPTP,CASECC,BITS,TWO1,CORE(2),CASE,CC,SS,
     1                ORF,ANDF,COREY(401)
      DIMENSION       ICASE(200,2),IWORD(200),IBIT(200)
      CHARACTER       UFM*23,UWM*25,UIM*29,SFM*25
      COMMON /XMSSG / UFM,UWM,UIM,SFM
      COMMON /TWO   / TWO1(32)
      COMMON /SYSTEM/ IBUF,NOUT
      COMMON /IFPX0 / LBD,LCC,BITS(1)
      COMMON /XIFP1 / IBLANK
      COMMON /ZZZZZZ/ COREX(1)
      COMMON /IFP1A / SCR1,CASECC,IS,NWPC,NCPW4,NMODES,ICC,NSET,NSYM,
     1                ZZZZBB,ISTR,ISUB,LENCC,IBEN,EQUAL,IEOR
      EQUIVALENCE    (COREX(1),COREY(1),ICASE(1,1)),(CORE(1),COREY(401))
      DATA    NAME  / 4HIFP1,  4HB          /
      DATA    CASE  , CC     / 4HCASE,4HCC  /
      DATA    SS    / 4HSS   /
      DATA    OPTP  / 4HOPTP /
      DATA    IWORD /
     1      -1,01,01,01,01,01,01,01,01,00,01,01,01,01,01,-1,00,01,01,00,
     2      01,01,00,01,01,00,01,01,00,01,01,00,01,01,00,01,01,01,-1,-1,
     3      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
     4      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
     5      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
     6      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,
     7      -1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,00,-1,-1,01,01,01,
     8      01,01,01,01,00,00,-1,01,01,-1,00,01,01,00,01,01,00,01,01,01,
     9      -1,-1,01,01,01,-1,00,01,01,00,01,01,00,01,01,00,01,01,01,00,
     X      01,01,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1,-1/
      DATA    IBIT  /
     1      00,02,03,04,05,06,07,08,09,10,10,10,13,14,15,00,18,18,18,18,
     2      18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,18,17,00,00,
     3      00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,
     4      00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,
     5      00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,
     6      00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,
     7      00,00,00,00,00,00,00,00,00,00,00,00,00,00,16,00,00,20,21,21,
     8      22,22,23,23,18,18,00,24,25,00,10,10,10,10,10,10,10,10,10,27,
     9      00,00,30,26,29,00,18,18,18,18,18,18,10,10,10,10,10,10,33,18,
     X      18,18,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00,00/
      DATA    NEW,DEBUG   / 2*.FALSE. /
C
      K  = LBD + 1
      IFIROD = 0
      IOLOOP = 0
      ILOOP  = 0
      IFIRST = 0
      IEOPTP = 0
C
C     ALLOCATE GINO BUFFERS
C
      NZ     = KORSZ(CORE)
      IBUF1  = NZ - IBUF + 1
      IBUF2  = IBUF1 - IBUF
      NZ     = NZ - 2*IBUF
      ICRQ   =-NZ
      IF (NZ .LE. 0) GO TO 700
      IECASE = 0
C
C     TRY TO FIND CASECC ON OPTP - TRY TO ASSUME PROPER POSITION
C
      CALL OPEN (*560,OPTP,CORE(IBUF1),2)
      IOPN = 0
C
C     FIND CASECC
C
   10 CALL READ (*600,*610,OPTP,CORE(1),2,1,IFLAG)
      IF (CORE(1).EQ.CASE .AND. CORE(2).EQ.CC) GO TO 30
      IF (CORE(1).EQ.CASE .AND. CORE(2).EQ.SS) GO TO 20
      IF (IOPN .EQ. 0) CALL REWIND (OPTP)
      IOPN = 1
      CALL SKPFIL (OPTP,1)
      GO TO 10
C
C     CASESS FOUND ON OPTP - SKIP TO CASECC
C
   20 CALL READ (*600,*610,OPTP,CORE(1),2,1,IFLAG)
      IF (CORE(1).NE.CASE .OR. CORE(2).NE.CC) GO TO 20
C
C     CASECC FOUND ON OLD PROB TAPE
C
C     OPEN CASECC AND SKIP CASESS IF PRESENT
C
   30 CALL OPEN (*650,CASECC,CORE(IBUF2),0)
   40 CALL READ (*670,*680,CASECC,CORE(1),2,1,IFLAG)
      IF (CORE(1).NE.CASE .OR. CORE(2).NE.CC) GO TO 40
      ASSIGN 50 TO IHOP
   50 CALL READ (*550,*610,OPTP,ICASE(1,2),LENCC,0,IFLAG)
      IF (ICASE(16,2) .EQ. 0) GO TO 60
      CALL FWDREC (*600,OPTP)
      GO TO 50
   60 IF (ICASE(LENCC,2) .EQ. 0) GO TO 70
      LSYM = ICASE(LENCC,2)
      CALL READ (*600,*610,OPTP,CORE(1),-LSYM,0,IFLAG )
   70 CALL READ (*600,*80 ,OPTP,CORE(1),   NZ,1,IFOPTP)
      ICRQ = NZ
      GO TO 700
   80 CALL READ (*510,*680,CASECC,ICASE(1,1),LENCC,0,IFLAG)
      IF (ICASE(16,1) .EQ. 0) GO TO 90
      CALL FWDREC (*670,CASECC)
      GO TO 80
   90 IF (ICASE(LENCC,1) .EQ. 0) GO TO 100
      LSYM = ICASE(LENCC,1)
      CALL READ (*670,*680,CASECC,CORE(IFOPTP+1),-LSYM,0,IFLAG)
  100 CALL READ (*670,*110,CASECC,CORE(IFOPTP+1),NZ-IFOPTP,1,IFCASE)
      ICRQ = NZ - IFOPTP
      GO TO 700
C
C     CHECK FOR LOOPING PROBLEM
C
  110 IF (IFIRST .NE. 0) GO TO 120
      IFIRST= 1
      ISPC  = ICASE(  3,1)
      IMPC  = ICASE(  2,1)
      IMTD  = ICASE(  5,1)
      IFREQ = ICASE( 14,1)
      ITFL  = ICASE( 15,1)
      IK1   = ICASE(139,1)
      IK2   = ICASE(140,1)
      IM1   = ICASE(141,1)
      IM2   = ICASE(142,1)
      IB1   = ICASE(143,1)
      IB2   = ICASE(144,1)
      IF (ICASE(165,1).GT.0 .OR. ICASE(164,1).GT.0) GO TO 130
      GO TO 140
  120 IF (ICASE(  3,1) .NE. ISPC) GO TO 130
      IF (ICASE(  2,1) .NE. IMPC) GO TO 130
      IF (ICASE(  5,1) .NE. IMTD) GO TO 130
      IF (ICASE(139,1).NE.IK1 .OR. ICASE(140,1).NE.IK2) GO TO 130
      IF (ICASE(141,1).NE.IM1 .OR. ICASE(142,1).NE.IM2) GO TO 130
      IF (ICASE(143,1).NE.IB1 .OR. ICASE(144,1).NE.IB2) GO TO 130
      IF (ICASE( 15,1) .NE. ITFL ) GO TO 130
      IF (ICASE( 14,1) .NE. IFREQ) GO TO 130
      IF (ICASE(138,1) .GT. 0) GO TO 130
      IF (ICASE( 38,1) .NE. 0) GO TO 130
      GO TO 140
C
C     SET LOOP$
C
  130 BITS(K) = ORF(BITS(K),TWO1(11))
      ILOOP = 1
  140 CONTINUE
C
C     DETERMINE IF OLD PROBLEM WOULD HAVE LOOPED
C
      IF (IFIROD .NE. 0) GO TO 150
      IFIROD= 1
      ISPC1 = ICASE(  3,2)
      IMPC1 = ICASE(  2,2)
      IMTD1 = ICASE(  5,2)
      IK11  = ICASE(139,2)
      IK21  = ICASE(140,2)
      IM11  = ICASE(141,2)
      IM21  = ICASE(142,2)
      IB11  = ICASE(143,2)
      IB21  = ICASE(144,2)
      ITFL1 = ICASE( 15,2)
      IFREQ1= ICASE( 14,2)
      IF (ICASE(164,2).GT.0 .OR. ICASE(165,2).GT.0) GO TO 160
      GO TO 170
C
C     SECOND RECORD APPLY LOOP RULES
C
  150 IF (ICASE(  3,2) .NE. ISPC1) GO TO 160
      IF (ICASE(  2,2) .NE. IMPC1) GO TO 160
      IF (ICASE(  5,2) .NE. IMTD1) GO TO 160
      IF (ICASE(139,2).NE.IK11 .OR. ICASE(140,2).NE.IK21) GO TO 160
      IF (ICASE(141,2).NE.IM11 .OR. ICASE(142,2).NE.IM21) GO TO 160
      IF (ICASE(143,2).NE.IB11 .OR. ICASE(144,2).NE.IB21) GO TO 160
      IF (ICASE(138,2) .GT. 0) GO TO 160
      IF (ICASE( 38,2) .NE. 0) GO TO 160
      IF (ICASE( 15,2) .NE. ITFL1 ) GO TO 160
      IF (ICASE( 14,2) .NE. IFREQ1) GO TO 160
      GO TO 170
  160 IOLOOP = 1
  170 CONTINUE
      IF (IECASE .NE. 1) GO TO 180
      IF (IOLOOP .EQ. 1) GO TO 530
      GO TO 520
C
C     CHECK FOR CHANGES -
C
  180 IF (IEOPTP .EQ. 1) IEOPTP = 2
      DO 500 I = 1,LENCC
      IF (IBIT(I) .EQ. 0) GO TO 500
      L = IBIT(I)
      IF (L.LE.32 .AND. ANDF(BITS(K),TWO1(L)).NE.0) GO TO 500
      IF (IWORD(I) .EQ. 0) GO TO 210
      IF (ICASE(I,1) .EQ. ICASE(I,2)) GO TO 500
  190 IF (L .GT. 32) GO TO 200
      BITS(K) = ORF(BITS(K),TWO1(L))
      GO TO 500
C
C     SECOND CASECC WORD
C
  200 L = L - 31
      BITS(K+1) = ORF(BITS(K+1),TWO1(L))
      GO TO 500
C
C     CHECK FOR PRESENCE OF PRINT AND PLOT REQUESTS
C
  210 IF (I.NE.135 .AND. ICASE(I,1).EQ.0) GO TO 500
      IF (I .NE. 135) GO TO 220
      IF (ICASE(I,1).NE.0 .OR. ICASE(I,2).NE.0) GO TO 190
      GO TO 500
  220 IF (IBIT(I) .EQ. 18) BITS(K+1) = ORF(BITS(K+1),TWO1(3))
      IF (IBIT(I) .EQ. 10) BITS(K+1) = ORF(BITS(K+1),TWO1(4))
      IF (IEOPTP  .EQ.  2) GO TO 190
      IF (ICASE(I,1).LT.0 .AND. ICASE(I,2).LT.0) GO TO 500
      IF (ICASE(I,1).LT.0 .AND. ICASE(I,2).GE.0) GO TO 190
      IF (ICASE(I,1).GT.0 .AND. ICASE(I,2).LE.0) GO TO 190
      IPCASE = IFOPTP + 1
  230 IF (IPCASE .GT. IFOPTP+IFCASE) GO TO 690
      IF (CORE(IPCASE) .EQ. ICASE(I,1)) GO TO 240
      IPCASE = IPCASE + CORE(IPCASE+1) + 2
      GO TO 230
  240 IPOPTP = 1
  250 IF (IPOPTP .GT. IFOPTP) GO TO 620
      IF (CORE(IPOPTP) .EQ. ICASE(I,2)) GO TO 260
      IPOPTP = IPOPTP + CORE(IPOPTP+1) + 2
      GO TO 250
  260 IQCASE = IFOPTP + IFCASE + 1
      IX = IPCASE
      IY = IQCASE
      ASSIGN 280 TO JUMP
      IF (DEBUG) WRITE (NOUT,270)
  270 FORMAT (/,' ------ NPTP PASS ------')
      GO TO 360
  280 IQOPTP = IY
      IX = IPOPTP
      ASSIGN 300 TO JUMP
      IF (DEBUG) WRITE (NOUT,290)
  290 FORMAT (/,' ------ OPTP PASS ------')
      GO TO 360
  300 LENG1 = IQOPTP - IQCASE
      LENG2 = IY - IQOPTP
      IF (DEBUG) WRITE (NOUT,310) CORE(IPCASE),LENG1,LENG2,IY,IQOPTP,
     1                            IQCASE
  310 FORMAT (//,' IFP1B/@310  CHECKING SETS',I9,' FROM NPTP AND OPTP',
     1        /5X,'LENG1,LENG2, IY,IQOPTP,IQCASE =', 2I5,3I7)
      IF (LENG1 .NE. LENG2) GO TO 340
      DO 320 MM = 1,LENG1
      IF (CORE(IQCASE+MM-1) .NE. CORE(IQOPTP+MM-1)) GO TO 340
  320 CONTINUE
      IF (DEBUG) WRITE (NOUT,330) CORE(IPCASE)
  330 FORMAT (' ... NO DIFFERENCES IN SET',I8)
      GO TO 500
  340 WRITE  (NOUT,350) UIM,CORE(IPCASE)
  350 FORMAT (A29,', SET',I9,' DEFINITION HAS BEEN CHANGED IN RESTART')
      GO TO 190
C
C     A NEW NON-EXPANDING METHOD IS IMPLEMENTED HERE BY  G.CAHN/UNISYS
C     8/91, IN CASE THE ORIGINAL LOGIC RUNS OUT OF CORE SPACE
C
C     THE NEW METHOD WILL CONCATINATE VARIATIONS OF SET DEFINITION TO
C     THE SIMPLEST FORM.  E.G. THE NEXT 3 LINES SPECIFY THE SAME SET
C     10 THRU 400000  (THIS IS THE SIMPLEST FORM)
C     10, 11, 12 THRU 400008, 400009, 400000
C     10 THRU 20, 21, 22, 23 THRU 200, 201 THRU 500, 501 502 THRU 400000
C
  360 IF (NEW) GO TO 420
      IN = CORE(IX+1)
      IX = IX + 2
      M  = 0
  370 M  = M + 1
      IF (M-IN) 380,400,490
  380 IF (CORE(IX+M) .GT. 0) GO TO 400
      M1 = CORE(IX+M-1)
      M2 =-CORE(IX+M  )
      ICRQ = IY + M2 - M1 - NZ
      IF (ICRQ .GT. 0) GO TO 410
      DO 390 MM = M1,M2
      CORE(IY) = MM
      IY = IY + 1
  390 CONTINUE
      M = M + 1
      GO TO 370
  400 ICRQ = IY - NZ
      IF (IY .GT. NZ) GO TO 700
      CORE(IY) = CORE(IX+M-1)
      IY = IY + 1
      GO TO 370
C
C     INSUFFICIENT CORE SPACE, SWITCH TO NEW METHOD
C
  410 NEW = .TRUE.
      GO TO 260
C
C     NEW LOGIC WITHOUT THRU RANGE EXPANSION
C
  420 IN = CORE(IX+1)
      IX = IX + 2
      M0 = IY
      CORE(IY) = CORE(IX)
      IY = IY + 1
      IF (IN .EQ. 1) GO TO 490
      CORE(IY) = CORE(IX+1)
      IF (CORE(IY) .EQ. CORE(IY-1)+1) CORE(IY) = -CORE(IY)
      IY = IY + 1
      M  = 1
  430 M  = M + 1
      IF (M .GE. IN) GO TO 470
      M1 = CORE(IX+M)
      M2 = IABS(M1)
      IF (DEBUG) WRITE (NOUT,440) M,IN,IX,IY,M1,CORE(IY-1)
  440 FORMAT (' @440   M,IN,IX,IY,M1,CORE(IY-1) =',6I8)
      IF (M1 .LT. 0) GO TO 450
      IF (M1 .NE. 1-CORE(IY-1)) GO TO 460
      CORE(IY-1) = -M2
      GO TO 430
  450 IF (CORE(IY-1) .GT. 0) GO TO 460
      CORE(IY-1) = -M2
      GO TO 430
  460 CORE(IY) = M1
      IF (M1 .EQ. CORE(IY-1)+1) CORE(IY) = -M2
      IY = IY + 1
      GO TO 430
  470 ICRQ = IY - NZ
      IF (IY .GT. NZ) GO TO 700
      M1 = IY - 1
      IF (DEBUG) WRITE (NOUT,480) CORE(IX-2),(CORE(J),J=M0,M1)
  480 FORMAT (/,' IFP1B/@480    SET',I8, /,(2X,15I8))
C
  490 GO TO JUMP, (280,300)
C
  500 CONTINUE
      GO TO IHOP, (50,80)
C
C     EOF ON CASECC
C
  510 CALL CLOSE (CASECC,1)
      IF (IEOPTP .NE. 0) GO TO 530
      IECASE = 1
      GO TO 150
  520 CALL READ (*530,*610,OPTP,ICASE(1,2),LENCC,1,IFLAG)
      IF (ICASE(16,2) .NE. 0) GO TO 520
      GO TO 150
  530 CALL CLOSE (OPTP,2)
      IF (IEOPTP.EQ.1 .OR. IOLOOP.EQ.0) GO TO 540
C
C     SET LOOP1  THIS SHOULD REEXECUTE THE ENTIRE LOOP
C
      BITS(K) = ORF(BITS(K),TWO1(12))
C
C     CHECK FOR LOOP$ IF NOT ON SET NOLOOP$
C
  540 IF (ILOOP .EQ. 0) BITS(K) = ORF(BITS(K),TWO1(32))
      RETURN
C
C     EOF ON  OPTP
C
  550 ASSIGN 80 TO IHOP
      IEOPTP = 1
      GO TO 80
C
C     ERROR MESSAGES
C
  560 IP1 = -1
  570 IP2 = OPTP
  580 CALL MESAGE (IP1,IP2,NAME)
      RETURN
C
  600 IP1 = -2
      GO TO 570
  610 IP1 = -3
      GO TO 570
  620 CORE(1) = OPTP
      CORE(2) = IBLANK
  630 WRITE  (NOUT,640) SFM,CORE(1),CORE(2)
  640 FORMAT (A25,' 651, LOGIC ERROR IN SUBROUTINE IFP1B WHILE ',
     1       'PROCESSING SET DATA ON ',2A4,' FILE.')
      IP1 = -37
      GO TO 580
  650 IP1 = -1
  660 IP2 = CASECC
      GO TO 580
  670 IP1 = -2
      GO TO 660
  680 IP1 = -3
      GO TO 660
  690 CORE(1) = CASE
      CORE(2) = CC
      GO TO 630
  700 IP1 = -8
      IP2 = ICRQ
      GO TO 580
      END


text:
Part pilot.  All geek.  Aviator, foodie, wino, computer guy.  Not necessarily in that order.

20110423 12:29:06 nbsp Welcome to the Wiki! Users/TomGarberson

20120324 20:31:44 nbsp Why did you remove the emails from The Ramble Apartments?  Do you have reason to know that those people dont want their email addresses on the page? Users/CovertProfessor

20120324 21:42:57 nbsp Much as you seem to have a preference for some anonymity, so do they.  Please respect it. Users/pilotgeek

20120324 21:51:23 nbsp Youll note that I have been respecting it, since I did not revert your change.  I was simply confirming that you were not trying to protect people who had not asked to be protected.  Try checking your attitude next time. Users/CovertProfessor

20120324 21:57:14 nbsp The question was shared by others... with no reason given, we dont know if are you removing those emails they added because they wish it, or because you are a competitor, or because you dont like the at symbol.  With no answer, I still dont know what to think about the edit. You may well have a great reason.  Could you please share it? Users/JabberWokky Evan JabberWokky Edwards

20120324 22:35:50 nbsp @cp That is always an issue with the electronic medium.  I intended clarify, you read attitude (none intended). Users/pilotgeek

20120324 22:46:16 nbsp @evan What question?  Who others?  Society of Wiki Editors? :) Seriously, I was asked to remove the addresses by the manager, as they prefer not to have their personal emails on a web page. Users/pilotgeek

20120324 23:27:04 nbsp Well, at least two people, cp and myself.  Thats a we.  I was just pointing out that s/he wasnt alone in wondering.  This seems somehow awfully formal and antagonistic for a simple query about an edit.  You removed stuff... people asked why.  Its not exactly an interrogation or asking you to defend your existence.  Just wondering why you deleted some stuff.  Thats all.  People ask me the same thing.  Heck, Im pretty sure Ive asked cp in the past why something was deleted. Users/JabberWokky

20120325 10:14:51 nbsp I didnt read in any attitude.  When you reminded me of my own anonymity, you were reminding me to be consistent, as though I might not be.  When you said, Please respect it, you implied that I might not respect their anonymity, so you had to request that I do so  even though I had not undone your edits and had specifically asked you if they had requested privacy.  In other words, you assumed bad behavior on my part even though all evidence suggested good behavior.  Thats insulting.  As JW notes, people remove information from the wiki all the time (something that is normally frowned upon), and not all of their reasons are good ones.  For example, within the last month someone changed the phone number of a Chinese restaurant to the number of a different Chinese restaurant in town.  So, we (meaning, we who care about the wiki remaining a useful resource) have to ask when information is changed or removed.

Look, I get it.  The Ramble Apartments opened up with a lot of fanfare, but based on the longwinded yet articulate comments on the page, it hasnt lived up to the hype.  I can certainly understand why, under those circumstances, people would not want their email addresses on the page.  Now, Im sure youll tell me that that isnt the reason, but that is how it appears from the outside.  They might want to consider that. Users/CovertProfessor

20120330 13:41:07 nbsp @cp OK.. I read this and gave myself some time to not be annoyed.  Seriously.. there was no attitude in my original comment.  I wrote that below and meant it.  I do care about the wiki as a useful resource, as I use it all the time and contribute occasionally.  There is nothing sinister going on in my removing the addresses, just contact preferences.

Since you mentioned the information on the Ramble page, I will say that the editorial at the top of all pages should be mostly factual, as opposed to opinion.  And at the ramble that is simply not the case, but I am refraining from editing that, as it will ultimately take care of itself (as the facility gets out of its infancy).  I have seen numerous times where tilted opinions are removed from the informative section, both by you.. and frankly by me as well.  Of course, comments are comments  opinions  and they belong in that location on the page.  I would never suggest ones freedom of speech should revoked. Users/pilotgeek
 I will say that the editorial at the top of all pages should be mostly factual, as opposed to opinion. You have not been paying attention. The informative/comment section distinction you perceive is an unfortunate artifact of people being more willing to use the comment box than the edit button. There are pages loaded with opinion that dont have a comment barthat was the default state of pages on this wiki for a long time in fact. Users/WilliamLewis

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
c     Subroutine convert date time string to jtime
c     AJ_Kettle, 11Feb2021

      SUBROUTINE get_jtime_from_datetime(s_obs_date_time,
     +     d_obs_jtime)

      IMPLICIT NONE
c************************************************************************
c     Variables passed into subroutine

      CHARACTER(*)        :: s_obs_date_time

      DOUBLE PRECISION    :: d_obs_jtime

c*****
c     Variables used within subroutine

      INTEGER             :: i,j,k,ii,jj,kk

      CHARACTER(LEN=4)    :: s_year
      CHARACTER(LEN=2)    :: s_month
      CHARACTER(LEN=2)    :: s_day
      CHARACTER(LEN=2)    :: s_hour
      CHARACTER(LEN=2)    :: s_minute
      CHARACTER(LEN=2)    :: s_second

      CHARACTER(LEN=10)   :: s_date_new
      CHARACTER(LEN=8)    :: s_time_new
c************************************************************************

      s_year    =s_obs_date_time(1:4)
      s_month   =s_obs_date_time(6:7)
      s_day     =s_obs_date_time(9:10)
      s_hour    =s_obs_date_time(12:13)
      s_minute  =s_obs_date_time(15:16)
      s_second  =s_obs_date_time(18:19)

      s_date_new=s_day//'/'//s_month//'/'//s_year
      s_time_new=s_hour//'/'//s_minute//'/'//s_second

      CALL str_to_dt_pvwave3(d_obs_jtime,
     +  s_date_new,s_time_new,2,-1)

c      print*,'s_date_new=',s_date_new
c      print*,'s_time_new=',s_time_new
c      print*,'d_obs_jtime=',d_obs_jtime

      RETURN
      END

max_stars_repo_name: glamod/glamod-nuim

max_issues_repo_name: glamod/glamod-nuim

max_forks_repo_name: glamod/glamod-nuim


        program testjdqz
        use myjdqz_mod, only :  myjdqz_setup
        implicit none

        integer, parameter :: r8 = kind(1.0d0)

        integer :: lwork
        double complex, allocatable, dimension(:,:) :: zwork
        double complex, allocatable, dimension(:) :: alpha,beta
        double complex, allocatable, dimension(:) :: v,Av,Bv,residu
        double complex, allocatable, dimension(:,:) :: eivec
        double complex :: ztarget

        real*8, allocatable, dimension(:,:,:) :: egn_vectors
        real*8, allocatable, dimension(:) :: rho, omega
	integer, allocatable, dimension(:) :: im_col, in_col
	integer :: ns, mn_col, i, k, jj, m
	real*8 :: enrg1, enrg2, center_freq

        integer :: kmax
        integer :: neq, jmin,jmax,maxstep, testspace
        integer :: order, method, gmres_m, cgstab_l, maxnmv
        real(kind=r8) :: tol, lock
        logical :: wanted

        real(kind=8), external :: dznrm2

        integer :: j,ierr
        character*80 fileA, fileB, arg1
        real(kind=r8) :: time_1, time_2, elapse
        CALL getarg(1, arg1)
        read(arg1,'(e12.4)') center_freq
	write(*,'("freq. = ",e12.4)') center_freq
        open(unit=15,file="jdqz.dat",status="unknown")
	read(15,*) ns, mn_col
c	write(*,*) ns,mn_col
	allocate(im_col(mn_col), in_col(mn_col), stat=ierr)
	allocate(rho(ns), stat=ierr)
	do i=1,mn_col
	 read(15,*) im_col(i), in_col(i)
c	 write(*,*) im_col(i), in_col(i)
	end do
	do i=1,ns
	 read(15,*) rho(i)
c	 write(*,*) rho(i)
	end do
	close(unit=15)
!       ztarget = (0.0d0,0.0d0)
!       ztarget = dcmplx(400.0d0,0.0d0)
        ztarget = dcmplx(center_freq**2,0.0d0)
!        ztarget = dcmplx(0.05*center_freq,center_freq)               !for TAEFL test problem
	write(*,'("eigenvalue target = ",e15.7)') real(ztarget)
!	write(*,'("eigenvalue target = ",2(e15.7,2x))') ztarget       !for TAEFL test problem
        fileA = 'a_matrix.dat'
        fileB = 'b_matrix.dat'
        call myjdqz_setup(ztarget, fileA,fileB,neq)

      tol = 1.d-9

!...  kmax is number of eigenvalues we want
      kmax = 40
      allocate(omega(kmax), stat=ierr)

!...  jmin,jmax size of search space
      jmax = min(150,max(30,3*kmax))
      jmin = min(jmax,2*kmax)

      maxstep = 200
!...   take it small to avoid missing eigensolutions
      lock = 1.d-9
!...     order =  0: nearest to target
!...     order = -1: smallest real part
!...     order =  1: largest real part
!...     order = -2: smallest complex part
!...     order =  2: largest complex part
      order = 0
!...     method = 1: gmres(m)
!...     method = 2: cgstab(l)
      method = 1
!...     for gmres(m):
      gmres_m = 99
!...     for cgstab(l):
      cgstab_l= 4
!...     maximum number of matvecs in cgstab or gmres
      maxnmv = 200

!...  testspace = 3 if a reasonable value for target is known, else
!...  set testspace = 2
!...  testspace = 1 seems to better match results in matlab
      testspace = 3
!...     Testspace 1: w = "Standard Petrov" * v (Section 3.1.1)
!...     Testspace 2: w = "Standard 'variable' Petrov" * v (Section 3.1.2)
!...     Testspace 3: w = "Harmonic Petrov" * v (Section 3.5.1)
!
      if ( method .eq. 1 ) then
!...   ---------------------------
!...   storage if GMRES(M) is used
!...   ---------------------------
         lwork =  4 +  gmres_m  + 5*jmax + 3*kmax
      else
!...    --------------------------------
!...    storage if BiCGSTAB(ell) is used
!...    --------------------------------
         lwork = 10 + 6*cgstab_l + 5*jmax + 3*kmax
      end if
!... converged eigenvectors if wanted = .true.
!...  eivec are converged eigenvectors if wanted = .true.,
!...  else converged Schur vectors
      wanted = .true.
!

      allocate( eivec(neq,kmax), alpha(jmax),beta(jmax),
     &          v(neq),Av(neq),Bv(neq),residu(neq),stat=ierr )
      if (ierr.ne.0) then
              write(*,*) 'allocate(eivec), neq,kmax,ierr ',
     &                    neq,kmax,ierr
              stop '** error ** '
      endif

      allocate( zwork(neq,lwork),stat=ierr)
      if (ierr.ne.0) then
              write(*,*) 'allocate(zwork), neq,lwork,ierr',
     &                     neq,lwork,ierr
              stop '** error ** '
      endif

      call cpu_time(time_1)
      call jdqz(alpha, beta, eivec, wanted, neq, ztarget, tol,
     &     kmax, jmax, jmin,
     &     method, gmres_m, cgstab_l, maxnmv, maxstep,
     &     lock, order, testspace, zwork, lwork )
!
      call cpu_time(time_2)
      elapse = time_2 - time_1

      deallocate( zwork, stat=ierr)
      if (ierr.ne.0) then
              write(*,*) 'deallocate(zwork), ierr ',ierr
              stop '** error ** '
      endif

!
!...     Compute the norms of the residuals:
      do j = 1, kmax
         v(1:neq) = eivec(1:neq,j)
         call amul(neq, v, Av )
         call bmul(neq, v, Bv )
         residu(1:neq) = beta(j)*Av(1:neq) - alpha(j)*Bv(1:neq)
         print '("lambda(",i2,"): (",1p,e11.4,",",e11.4,
     &        " )")', j,alpha(j)/beta(j)
         print '(a30,d13.6)', '||beta Ax - alpha Bx||:',
     &          dznrm2( neq, residu, 1 )
      end do
      write(*,10)  elapse, elapse
!
   10 format(1x,'END JDQZ AFTER ',f6.2,' SEC. CPU-TIME AND ', f6.2,
     &       ' SEC. ELAPSED TIME' )

      allocate(egn_vectors(kmax,mn_col,ns), stat=ierr)
      do k=1,kmax
       if(beta(k) .ne. 0.) omega(k)=sqrt(real(alpha(k)/beta(k)))
       jj = 0
       do i=1,ns
        do j = 1, mn_col
         jj = jj + 1
	 egn_vectors(k,j,i) = real(eivec(jj,k))
	end do
       end do
      end do

      write(*,'("Freq. Range: ",e12.4," to " e12.4)')
     >   minval(omega),maxval(omega)
      open(unit=33,file="egn_mode_asci.dat",status="unknown")
      open(unit=34,file="egn_values.dat",status="unknown")
      enrg1 = 1.; enrg2 = 1.   !dummy values for now
       write(33,*) kmax
       write(33,*) mn_col
       write(33,*) ns
       do i=1,mn_col
        write(33,*) im_col(i)
        write(33,*) in_col(i)
       end do
       do i=1,kmax
        write(33,*) real(alpha(i)/beta(i))
        write(34,'(e22.14,2(2x,e15.7))')
     >       dble(alpha(i)/beta(i)), enrg1, enrg2
       end do
       do i=1,ns
        write(33,*) rho(i)
       end do
       do j=1,kmax
        do m = 1,ns
         do i=1,mn_col
          write(33,*) egn_vectors(j,i,m)
          end do
         end do
       end do
       close(unit=33)

      stop
      end program

        subroutine precon(neq,q)
        use myjdqz_mod, only : m,n, L,D,ipiv,U
        use symtrd_mod, only : symtrd_solve
!       ----------------------------------
!       subroutine to compute q = inv(K)*q
!       ----------------------------------
        integer, parameter :: r8 = kind(1.0d0)

        integer neq
        double complex q(neq)

        integer, parameter :: nrhs = 2
        real(kind=r8) :: Brhs(m,n,nrhs)
        integer :: ip,i,k, info

        logical :: isok

        isok = (neq .eq. m*n)
        if (.not.isok) then
                write(*,*) 'precon:invalid neq ',neq
        endif

!       -------------------------------------
!       convert q(1:(m*n)) to Brhs(1:m,1:n,1)
!       -------------------------------------
        ip = 1
        do k=1,n
        do i=1,m
           ip =  i + (k-1)*m
           Brhs(i,k,1) = dble(q(ip))
           Brhs(i,k,2) = dimag(q(ip))
        enddo
        enddo
        call symtrd_solve( m,n,nrhs, L,D,ipiv,U, Brhs,info)
        if (info.ne.0) then
                write(*,*) 'precon: symtrd_solve return info',info
        endif

!       ------------------------------------------
!       convert Brhs(1:m,1:n,1) back to q(1:(m*n))
!       ------------------------------------------

        ip = 1
        do k=1,n
        do i=1,m
          ip = i + (k-1)*m
          q(ip) = dcmplx(Brhs(i,k,1),Brhs(i,k,2))
        enddo
        enddo
        return
        end subroutine precon
        subroutine amul( neq, q, r )
        use myjdqz_mod, only : m,n, Adiag,Asub,Asup
        use symtrd_mod, only : symtrd_matmul
!       -----------------------------
!       subroutine to compute r = A*q
!       -----------------------------
        integer, parameter :: r8 = kind(1.0d0)

        integer neq
        double complex q(neq),r(neq)

        integer, parameter :: nrhs = 2
        real(kind=r8), dimension(m,n,nrhs) :: X,Y
        integer :: i,k,ip

!       ---------------------------------------
!       convert q(1:(m*n)) to X(1:m,1:n,1:nrhs)
!       ---------------------------------------
        ip = 1
        do k=1,n
        do i=1,m
           ip = i + (k-1)*m
           X(i,k,1) = dble(q(ip))
           X(i,k,2) = dimag(q(ip))
        enddo
        enddo

        call symtrd_matmul( m,n,nrhs, Adiag,Asub,Asup, X,Y )

!       ---------------------------------------
!       convert Y(1:m,1:n,1:nrhs) to r(1:(m*n))
!       ---------------------------------------

         ip = 1
         do k=1,n
         do i=1,m
           ip = i + (k-1)*m
           r(ip) = dcmplx(Y(i,k,1),Y(i,k,2))
         enddo
         enddo

         return
         end subroutine amul
        subroutine bmul( neq, q, r )
        use myjdqz_mod, only : m,n, Bdiag,Bsub,Bsup
        use symtrd_mod, only : symtrd_matmul
!       -----------------------------
!       subroutine to compute r = B*q
!       -----------------------------
        integer, parameter :: r8 = kind(1.0d0)
        integer neq
        double complex q(neq),r(neq)

        integer, parameter :: nrhs = 2
        real(kind=8), dimension(m,n,nrhs) :: X,Y
        integer :: i,k,ip

!       ---------------------------------------
!       convert q(1:(m*n)) to X(1:m,1:n,1:nrhs)
!       ---------------------------------------
        ip = 1
        do k=1,n
        do i=1,m
           ip = i + (k-1)*m
           X(i,k,1) = dble(q(ip))
           X(i,k,2) = dimag(q(ip))
        enddo
        enddo

        call symtrd_matmul( m,n,nrhs, Bdiag,Bsub,Bsup, X,Y )

!       ---------------------------------------
!       convert Y(1:m,1:n,1:nrhs) to r(1:(m*n))
!       ---------------------------------------

         ip = 1
         do k=1,n
         do i=1,m
           ip = i + (k-1)*m
           r(ip) = dcmplx(Y(i,k,1),Y(i,k,2))
         enddo
         enddo

         return
         end subroutine bmul


text:
../../src/cacode.v
../../src/gps.v
../../src/pcode.v
../../../aes/src/aes_128.v
../../../aes/src/aes_192.v
../../../aes/src/aes_top_axi4lite.v
../../../aes/src/aes_top_wb.v
../../../aes/src/round.v
../../../aes/src/table.v

max_stars_repo_name: animeshbchowdhury/robust-pnr-time

max_issues_repo_name: animeshbchowdhury/robust-pnr-time

max_forks_repo_name: ieee-ceda-datc/datc-rdf


text:
!-----------------------------------------------------------------------
! Basic parallel PIC library for MPI communications with OpenMP
! pplib2.f90 contains basic communications procedures for 1d partitions:
! PPINIT2 initializes parallel processing for Fortran90, returns
!         number of processors and processor id.
! PPEXIT terminates parallel processing.
! PPABORT aborts parallel processing.
! PWTIMERA performs parallel local wall clock timing.
! PPSUM performs parallel sum of a real vector.
! PPDSUM performs parallel sum of a double precision vector.
! PPIMAX performs parallel maximum of an integer vector.
! PPDMAX performs parallel maximum of a double precision vector.
! PPNCGUARD2L copies data to guard cells in y for scalar data, linear
!             interpolation, and distributed data with non-uniform
!             partition.
! PPNAGUARD2L adds guard cells in y for scalar array, linear
!             interpolation, and distributed data with non-uniform
!             partition.
! PPNACGUARD2L adds guard cells in y for vector array, linear
!              interpolation, and distributed data with non-uniform
!              partition.
! PPTPOSE performs a transpose of a complex scalar array, distributed
!         in y, to a complex scalar array, distributed in x.
! PPNTPOSE performs a transpose of an n component complex vector array,
!          distributed in y, to an n component complex vector array,
!          distributed in x.
! PPPMOVE2 moves particles into appropriate spatial regions for tiled
!          distributed data.
! written by viktor k. decyk, ucla
! copyright 1995, regents of the university of california
! update: may 9, 2015
      module mpplib2
      use mpi
      implicit none
!
! common data for parallel processing
! lstat = length of status array
      integer, parameter :: lstat = MPI_STATUS_SIZE
! nproc = number of real or virtual processors obtained
! lgrp = current communicator
! mreal = default datatype for reals
! mint = default datatype for integers
! mcplx = default datatype for complex type
! mdouble = default double precision type
! lworld = MPI_COMM_WORLD communicator
      integer :: nproc, lgrp, mreal, mint, mcplx, mdouble, lworld
! msum = MPI_SUM
! mmax = MPI_MAX
      integer :: msum, mmax
      save
!
      private
      public :: PPINIT2, PPEXIT, PPABORT, PWTIMERA
      public :: PPSUM, PPDSUM, PPIMAX, PPDMAX
      public :: PPNCGUARD2L, PPNAGUARD2L, PPNACGUARD2L
      public :: PPTPOSE, PPNTPOSE, PPPMOVE2
!
      contains
!
!-----------------------------------------------------------------------
      subroutine PPINIT2(idproc,nvp)
! this subroutine initializes parallel processing
! lgrp communicator = MPI_COMM_WORLD
! output: idproc, nvp
! idproc = processor id in lgrp communicator
! nvp = number of real or virtual processors obtained
      implicit none
      integer, intent(inout) :: idproc, nvp
! nproc = number of real or virtual processors obtained
! lgrp = current communicator
! mreal = default datatype for reals
! mint = default datatype for integers
! mcplx = default datatype for complex type
! mdouble = default double precision type
! lworld = MPI_COMM_WORLD communicator
! msum = MPI_SUM
! mmax = MPI_MAX
! local data
      integer :: ierror, ndprec, idprec
      integer :: iprec
      logical :: flag
      real :: prec
! ndprec = (0,1) = (no,yes) use (normal,autodouble) precision
      if (digits(prec) > 24) then
         ndprec = 1
      else
         ndprec = 0
      endif
! idprec = (0,1) = (no,yes) use (normal,autodouble) integer precision
      if (digits(iprec) > 31) then
         idprec = 1
      else
         idprec = 0
      endif
! this segment is used for mpi computers
! indicate whether MPI_INIT has been called
      call MPI_INITIALIZED(flag,ierror)
      if (.not.flag) then
! initialize the MPI execution environment
         call MPI_INIT(ierror)
         if (ierror /= 0) stop
      endif
      lworld = MPI_COMM_WORLD
      lgrp = lworld
! determine the rank of the calling process in the communicator
      call MPI_COMM_RANK(lgrp,idproc,ierror)
! determine the size of the group associated with a communicator
      call MPI_COMM_SIZE(lgrp,nproc,ierror)
! set default datatypes
      mint = MPI_INTEGER
      mdouble = MPI_DOUBLE_PRECISION
! single precision real
      if (ndprec==0) then
         mreal = MPI_REAL
         mcplx = MPI_COMPLEX
! double precision real
      else
         mreal = MPI_DOUBLE_PRECISION
         mcplx = MPI_DOUBLE_COMPLEX
      endif
! single precision integer
!     if (idprec==0) then
!        mint = MPI_INTEGER
! double precision integer
!     else
!        mint = MPI_INTEGER8
!     endif
! operators
      msum = MPI_SUM
      mmax = MPI_MAX
      nvp = nproc
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPEXIT()
! this subroutine terminates parallel processing
      implicit none
! lworld = MPI_COMM_WORLD communicator
! local data
      integer :: ierror
      logical :: flag
! indicate whether MPI_INIT has been called
      call MPI_INITIALIZED(flag,ierror)
      if (flag) then
! synchronize processes
         call MPI_BARRIER(lworld,ierror)
! terminate MPI execution environment
         call MPI_FINALIZE(ierror)
      endif
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPABORT()
! this subroutine aborts parallel processing
      implicit none
! lworld = MPI_COMM_WORLD communicator
! local data
      integer :: errorcode, ierror
      logical :: flag
! indicate whether MPI_INIT has been called
      call MPI_INITIALIZED(flag,ierror)
      if (flag) then
         errorcode = 1
! terminate MPI execution environment
         call MPI_ABORT(lworld,errorcode,ierror)
      endif
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PWTIMERA(icntrl,time,dtime)
! this subroutine performs local wall clock timing
! input: icntrl, dtime
! icntrl = (-1,0,1) = (initialize,ignore,read) clock
! clock should be initialized before it is read!
! time = elapsed time in seconds
! dtime = current time
! written for mpi
      implicit none
      integer, intent(in) :: icntrl
      real, intent(inout) :: time
      double precision, intent(inout) :: dtime
! local data
      double precision :: jclock
! initialize clock
      if (icntrl==(-1)) then
         dtime = MPI_WTIME()
! read clock and write time difference from last clock initialization
      else if (icntrl==1) then
         jclock = dtime
         dtime = MPI_WTIME()
         time = real(dtime - jclock)
      endif
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPSUM(f,g,nxp)
! this subroutine performs a parallel sum of a vector, that is:
! f(j,k) = sum over k of f(j,k)
! at the end, all processors contain the same summation.
! f = input and output real data
! g = scratch real array
! nxp = number of data values in vector
      implicit none
      integer, intent(in) :: nxp
      real, dimension(nxp), intent(inout) :: f, g
! lgrp = current communicator
! mreal = default datatype for reals
! msum = MPI_SUM
! local data
      integer :: j, ierr
! perform sum
      call MPI_ALLREDUCE(f,g,nxp,mreal,msum,lgrp,ierr)
! copy output from scratch array
      do j = 1, nxp
         f(j) = g(j)
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPDSUM(f,g,nxp)
! this subroutine performs a parallel sum of a vector, that is:
! f(j,k) = sum over k of f(j,k)
! at the end, all processors contain the same summation.
! f = input and output double precision data
! g = scratch double precision array
! nxp = number of data values in vector
      implicit none
      integer, intent(in) :: nxp
      double precision, dimension(nxp), intent(inout) :: f, g
! lgrp = current communicator
! mdouble = default double precision type
! msum = MPI_SUM
! local data
      integer :: j, ierr
! perform sum
      call MPI_ALLREDUCE(f,g,nxp,mdouble,msum,lgrp,ierr)
! copy output from scratch array
      do j = 1, nxp
         f(j) = g(j)
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPIMAX(if,ig,nxp)
! this subroutine finds parallel maximum for each element of a vector
! that is, if(j,k) = maximum as a function of k of if(j,k)
! at the end, all processors contain the same maximum.
! if = input and output integer data
! ig = scratch integer array
! nxp = number of data values in vector
      implicit none
      integer, intent(in) :: nxp
      integer, dimension(nxp), intent(inout) :: if, ig
! lgrp = current communicator
! mint = default datatype for integers
! mmax = MPI_MAX
! local data
      integer :: j, ierr
! find maximum
      call MPI_ALLREDUCE(if,ig,nxp,mint,mmax,lgrp,ierr)
! copy output from scratch array
      do j = 1, nxp
         if(j) = ig(j)
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPDMAX(f,g,nxp)
! this subroutine finds parallel maximum for each element of a vector
! that is, f(j,k) = maximum as a function of k of f(j,k)
! at the end, all processors contain the same maximum.
! f = input and output double precision data
! g = scratch double precision array
! nxp = number of data values in vector
      implicit none
      integer, intent(in) :: nxp
      double precision, dimension(nxp), intent(inout) :: f, g
! lgrp = current communicator
! mdouble = default double precision type
! mmax = MPI_MAX
! local data
      integer j, ierr
! find maximum
      call MPI_ALLREDUCE(f,g,nxp,mdouble,mmax,lgrp,ierr)
! copy output from scratch array
      do j = 1, nxp
         f(j) = g(j)
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNCGUARD2L(f,nyp,kstrt,nvp,nxv,nypmx)
! this subroutine copies data to guard cells in non-uniform partitions
! f(j,k) = real data for grid j,k in particle partition.
! the grid is non-uniform and includes one extra guard cell.
! output: f
! nyp = number of primary gridpoints in field partition
! it is assumed the nyp > 0.
! kstrt = starting data block number
! nvp = number of real or virtual processors
! nxv = first dimension of f, must be >= nx
! nypmx = maximum size of field partition, including guard cell.
! linear interpolation, for distributed data
      implicit none
      integer, intent(in) :: nyp, kstrt, nvp, nxv, nypmx
      real, dimension(nxv,nypmx), intent(inout) :: f
! lgrp = current communicator
! mreal = default datatype for reals
! local data
      integer :: j, ks, moff, kl, kr
      integer :: msid, ierr
      integer, dimension(lstat) :: istatus
! special case for one processor
      if (nvp==1) then
         do j = 1, nxv
            f(j,nyp+1) = f(j,1)
         enddo
         return
      endif
      ks = kstrt - 1
      moff = nypmx*nvp + 2
! copy to guard cells
      kr = ks + 1
      if (kr >= nvp) kr = kr - nvp
      kl = ks - 1
      if (kl < 0)  kl = kl + nvp
      ks = nyp + 1
! this segment is used for mpi computers
      call MPI_IRECV(f(1,ks),nxv,mreal,kr,moff,lgrp,msid,ierr)
      call MPI_SEND(f,nxv,mreal,kl,moff,lgrp,ierr)
      call MPI_WAIT(msid,istatus,ierr)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNAGUARD2L(f,scr,nyp,nx,kstrt,nvp,nxv,nypmx)
! this subroutine adds data from guard cells in non-uniform partitions
! f(j,k) = real data for grid j,k in particle partition.
! the grid is non-uniform and includes one extra guard cell.
! output: f, scr
! scr(j) = scratch array for particle partition
! nyp = number of primary gridpoints in particle partition
! it is assumed the nyp > 0.
! kstrt = starting data block number
! nvp = number of real or virtual processors
! nx = system length in x direction
! nxv = first dimension of f, must be >= nx
! nypmx = maximum size of field partition, including guard cells.
! linear interpolation, for distributed data
      implicit none
      integer, intent(in) :: nyp, kstrt, nvp, nx, nxv, nypmx
      real, dimension(nxv,nypmx), intent(inout) :: f
      real, dimension(nxv), intent(inout) :: scr
! lgrp = current communicator
! mreal = default datatype for reals
! local data
      integer :: j, nx1, ks, moff, kl, kr
      integer :: msid, ierr
      integer, dimension(lstat) :: istatus
      nx1 = nx + 1
! special case for one processor
      if (nvp==1) then
         do j = 1, nx1
            f(j,1) = f(j,1) + f(j,nyp+1)
            f(j,nyp+1) = 0.
         enddo
         return
      endif
      ks = kstrt - 1
      moff = nypmx*nvp + 1
! add guard cells
      kr = ks + 1
      if (kr >= nvp) kr = kr - nvp
      kl = ks - 1
      if (kl < 0) kl = kl + nvp
      ks = nyp + 1
! this segment is used for mpi computers
      call MPI_IRECV(scr,nxv,mreal,kl,moff,lgrp,msid,ierr)
      call MPI_SEND(f(1,ks),nxv,mreal,kr,moff,lgrp,ierr)
      call MPI_WAIT(msid,istatus,ierr)
! add up the guard cells
      do j = 1, nx1
         f(j,1) = f(j,1) + scr(j)
         f(j,nyp+1) = 0.0
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNACGUARD2L(f,scr,nyp,nx,ndim,kstrt,nvp,nxv,nypmx)
! this subroutine adds data from guard cells in non-uniform partitions
! f(ndim,j,k) = real data for grid j,k in particle partition.
! the grid is non-uniform and includes one extra guard cell.
! output: f, scr
! scr(ndim,j) = scratch array for particle partition
! nyp = number of primary gridpoints in particle partition
! it is assumed the nyp > 0.
! kstrt = starting data block number
! nvp = number of real or virtual processors
! nx = system length in x direction
! ndim = leading dimension of array f
! nxv = first dimension of f, must be >= nx
! nypmx = maximum size of field partition, including guard cells.
! linear interpolation, for distributed data
      implicit none
      integer, intent(in) ::  nyp, kstrt, nvp, nx, ndim, nxv, nypmx
      real, dimension(ndim,nxv,nypmx), intent(inout) :: f
      real, dimension(ndim,nxv), intent(inout) :: scr
! lgrp = current communicator
! mreal = default datatype for reals
! local data
      integer :: j, n, nx1, ks, moff, kl, kr
      integer :: nnxv
      integer :: msid, ierr
      integer, dimension(lstat) :: istatus
      nx1 = nx + 1
! special case for one processor
      if (nvp==1) then
         do j = 1, nx1
            do n = 1, ndim
               f(n,j,1) = f(n,j,1) + f(n,j,nyp+1)
               f(n,j,nyp+1) = 0.0
            enddo
         enddo
         return
      endif
      ks = kstrt - 1
      moff = nypmx*nvp + 1
      nnxv = ndim*nxv
! add guard cells
      kr = ks + 1
      if (kr >= nvp) kr = kr - nvp
      kl = ks - 1
      if (kl < 0) kl = kl + nvp
      ks = nyp + 1
! this segment is used for mpi computers
      call MPI_IRECV(scr,nnxv,mreal,kl,moff,lgrp,msid,ierr)
      call MPI_SEND(f(1,1,ks),nnxv,mreal,kr,moff,lgrp,ierr)
      call MPI_WAIT(msid,istatus,ierr)
! add up the guard cells
      do j = 1, nx1
         do n = 1, ndim
            f(n,j,1) = f(n,j,1) + scr(n,j)
            f(n,j,nyp+1) = 0.0
         enddo
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPTPOSE(f,g,s,t,nx,ny,kxp,kyp,kstrt,nvp,nxv,nyv,kxpd,  &
     &kypd)
! this subroutine performs a transpose of a matrix f, distributed in y,
! to a matrix g, distributed in x, that is,
! g(k+kyp*(m-1),j,l) = f(j+kxp*(l-1),k,m), where
! 1 <= j <= kxp, 1 <= k <= kyp, 1 <= l <= nx/kxp, 1 <= m <= ny/kyp
! and where indices l and m can be distributed across processors.
! this subroutine sends and receives one message at a time, either
! synchronously or asynchronously. it uses a minimum of system resources
! f = complex input array
! g = complex output array
! s, t = complex scratch arrays
! nx/ny = number of points in x/y
! kxp/kyp = number of data values per block in x/y
! kstrt = starting data block number
! nvp = number of real or virtual processors
! nxv/nyv = first dimension of f/g
! kypd/kxpd = second dimension of f/g
      implicit none
      integer, intent(in) :: nx, ny, kxp, kyp, kstrt, nvp, nxv, nyv
      integer, intent(in) :: kxpd, kypd
      complex, dimension(nxv,kypd), intent(in) :: f
      complex, dimension(nyv,kxpd), intent(inout) :: g
      complex, dimension(kxp*kyp), intent(inout) :: s, t
! lgrp = current communicator
! mcplx = default datatype for complex
! local data
      integer :: n, j, k, ks, kxps, kyps, kxyp, id, joff, koff, ld
      integer :: ierr, msid
      integer, dimension(lstat) :: istatus
      ks = kstrt - 1
      kxps = min(kxp,max(0,nx-kxp*ks))
      kyps = min(kyp,max(0,ny-kyp*ks))
      kxyp = kxp*kyp
! special case for one processor
      if (nvp==1) then
!$OMP PARALLEL DO PRIVATE(j,k)
         do k = 1, kyp
            do j = 1, kxp
               g(k,j) = f(j,k)
            enddo
         enddo
!$OMP END PARALLEL DO
         return
      endif
! this segment is used for shared memory computers
!     do m = 1, min(ny,nvp)
!        koff = kyp*(m - 1)
!        do k = 1, min(kyp,max(0,ny-koff))
!           do l = 1, min(nx,nvp)
!              joff = kxp*(l - 1)
!              do j = 1, min(kxp,max(0,nx-joff))
!                 g(k+koff,j+joff) = f(j+joff,k+koff)
!              enddo
!           enddo
!        enddo
!     enddo
! this segment is used for mpi computers
      do n = 1, nvp
         id = n - ks - 1
         if (id < 0) id = id + nvp
! extract data to send
         joff = kxp*id
         ld = min(kxp,max(0,nx-joff))
!$OMP PARALLEL DO PRIVATE(j,k)
         do k = 1, kyps
            do j = 1, ld
               s(j+ld*(k-1)) = f(j+joff,k)
            enddo
         enddo
!$OMP END PARALLEL DO
         ld = ld*kyps
! post receive
         call MPI_IRECV(t,kxyp,mcplx,id,n,lgrp,msid,ierr)
! send data
         call MPI_SEND(s,ld,mcplx,id,n,lgrp,ierr)
! receive data
         call MPI_WAIT(msid,istatus,ierr)
! insert data received
         koff = kyp*id
         ld = min(kyp,max(0,ny-koff))
!$OMP PARALLEL DO PRIVATE(j,k)
         do k = 1, ld
            do j = 1, kxps
               g(k+koff,j) = t(j+kxps*(k-1))
            enddo
         enddo
!$OMP END PARALLEL DO
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNTPOSE(f,g,s,t,nx,ny,kxp,kyp,kstrt,nvp,ndim,nxv,nyv, &
     &kxpd,kypd)
! this subroutine performs a transpose of a matrix f, distributed in y,
! to a matrix g, distributed in x, that is,
! g(1:ndim,k+kyp*(m-1),j,l) = f(1:ndim,j+kxp*(l-1),k,m), where
! 1 <= j <= kxp, 1 <= k <= kyp, 1 <= l <= nx/kxp, 1 <= m <= ny/kyp
! and where indices l and m can be distributed across processors.
! this subroutine sends and receives one message at a time, either
! synchronously or asynchronously. it uses a minimum of system resources
! f = complex input array
! g = complex output array
! s, t = complex scratch arrays
! nx/ny = number of points in x/y
! kxp/kyp = number of data values per block in x/y
! kstrt = starting data block number
! nvp = number of real or virtual processors
! ndim = leading dimension of arrays f and g
! nxv/nyv = first dimension of f/g
! kypd/kxpd = second dimension of f/g
      implicit none
      integer, intent(in) :: nx, ny, kxp, kyp, kstrt, nvp, ndim
      integer, intent(in) :: nxv, nyv, kxpd, kypd
      complex, dimension(ndim,nxv,kypd), intent(in) :: f
      complex, dimension(ndim,nyv,kxpd), intent(inout) :: g
      complex, dimension(ndim,kxp*kyp), intent(inout) :: s, t
! lgrp = current communicator
! mcplx = default datatype for complex
! local data
      integer :: i, n, j, k, ks, kxps, kyps, kxyp, id, joff, koff, ld
      integer :: ierr, msid
      integer, dimension(lstat) :: istatus
      ks = kstrt - 1
      kxps = min(kxp,max(0,nx-kxp*ks))
      kyps = min(kyp,max(0,ny-kyp*ks))
      kxyp = ndim*kxp*kyp
! special case for one processor
      if (nvp==1) then
!$OMP PARALLEL DO PRIVATE(i,j,k)
         do k = 1, kyp
            do j = 1, kxp
               do i = 1, ndim
                  g(i,k,j) = f(i,j,k)
               enddo
            enddo
         enddo
!$OMP END PARALLEL DO
         return
      endif
! this segment is used for shared memory computers
!     do m = 1, min(ny,nvp)
!        koff = kyp*(m - 1)
!        do k = 1, min(kyp,max(0,ny-koff))
!           do l = 1, min(nx,nvp)
!              joff = kxp*(l - 1)
!              do j = 1, min(kxp,max(0,nx-joff))
!                 do i = 1, ndim
!                    g(i,k+koff,j+joff) = f(i,j+joff,k+koff)
!                 enddo
!              enddo
!           enddo
!        enddo
!     enddo
! this segment is used for mpi computers
      do n = 1, nvp
         id = n - ks - 1
         if (id.lt.0) id = id + nvp
! extract data to send
         joff = kxp*id
         ld = min(kxp,max(0,nx-joff))
!$OMP PARALLEL DO PRIVATE(i,j,k)
         do k = 1, kyps
            do j = 1, ld
               do i = 1, ndim
                  s(i,j+ld*(k-1)) = f(i,j+joff,k)
               enddo
            enddo
         enddo
!$OMP END PARALLEL DO
         ld = ndim*ld*kyps
! post receive
         call MPI_IRECV(t,kxyp,mcplx,id,n,lgrp,msid,ierr)
! send data
         call MPI_SEND(s,ld,mcplx,id,n,lgrp,ierr)
! receive data
         call MPI_WAIT(msid,istatus,ierr)
! insert data received
         koff = kyp*id
         ld = min(kyp,max(0,ny-koff))
!$OMP PARALLEL DO PRIVATE(i,j,k)
         do k = 1, ld
            do j = 1, kxps
               do i = 1, ndim
                  g(i,k+koff,j) = t(i,j+kxps*(k-1))
               enddo
            enddo
         enddo
!$OMP END PARALLEL DO
      enddo
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPPMOVE2(sbufr,sbufl,rbufr,rbufl,ncll,nclr,mcll,mclr,  &
     &kstrt,nvp,idimp,nbmax,mx1)
! this subroutine moves particles into appropriate spatial regions
! for distributed data, with 1d domain decomposition in y.
! tiles are assumed to be arranged in 2D linear memory
! output: rbufr, rbufl, mcll, mclr
! sbufl = buffer for particles being sent to lower processor
! sbufr = buffer for particles being sent to upper processor
! rbufl = buffer for particles being received from lower processor
! rbufr = buffer for particles being received from upper processor
! ncll = particle number being sent to lower processor
! nclr = particle number being sent to upper processor
! mcll = particle number being received from lower processor
! mclr = particle number being received from upper processor
! kstrt = starting data block number
! nvp = number of real or virtual processors
! idimp = size of phase space = 4 or 5
! nbmax =  size of buffers for passing particles between processors
! mx1 = (system length in x direction - 1)/mx + 1
      implicit none
      integer, intent(in) :: kstrt, nvp, idimp, nbmax, mx1
      real, dimension(idimp,nbmax), intent(in) :: sbufl, sbufr
      real, dimension(idimp,nbmax), intent(inout) :: rbufl, rbufr
      integer, dimension(3,mx1), intent(in) :: ncll, nclr
      integer, dimension(3,mx1), intent(inout) :: mcll, mclr
! lgrp = current communicator
! mint = default datatype for integers
! mreal = default datatype for reals
! local data
      integer :: ierr, ks, kl, kr, i, j, jsl, jsr
      integer :: nbsize, ncsize
      integer, dimension(8) :: msid
      integer, dimension(4) :: itg
      integer, dimension(lstat) :: istatus
      data itg /3,4,5,6/
      ks = kstrt - 1
      nbsize = idimp*nbmax
      ncsize = 3*mx1
! copy particle buffers: update rbufl, rbufr, mcll, mclr
! special case for one processor
      if (nvp==1) then
         do j = 1, mx1
            do i = 1, 3
               mcll(i,j) = nclr(i,j)
            enddo
         continue
         enddo
         do j = 1, mx1
            do i = 1, 3
               mclr(i,j) = ncll(i,j)
            enddo
         enddo
         do j = 1, nclr(3,mx1)
            do i = 1, idimp
               rbufl(i,j) = sbufr(i,j)
            enddo
         enddo
         do j = 1, ncll(3,mx1)
            do i = 1, idimp
               rbufr(i,j) = sbufl(i,j)
            enddo
         enddo
! this segment is used for mpi computers
      else
! get particles from below and above
         kr = ks + 1
         if (kr >= nvp) kr = kr - nvp
         kl = ks - 1
         if (kl < 0) kl = kl + nvp
! post receives
         call MPI_IRECV(mcll,ncsize,mint,kl,itg(1),lgrp,msid(1),ierr)
         call MPI_IRECV(mclr,ncsize,mint,kr,itg(2),lgrp,msid(2),ierr)
         call MPI_IRECV(rbufl,nbsize,mreal,kl,itg(3),lgrp,msid(3),ierr)
         call MPI_IRECV(rbufr,nbsize,mreal,kr,itg(4),lgrp,msid(4),ierr)
! send particle number offsets
         call MPI_ISEND(nclr,ncsize,mint,kr,itg(1),lgrp,msid(5),ierr)
         call MPI_ISEND(ncll,ncsize,mint,kl,itg(2),lgrp,msid(6),ierr)
         call MPI_WAIT(msid(1),istatus,ierr)
         call MPI_WAIT(msid(2),istatus,ierr)
! send particles
         jsr = idimp*nclr(3,mx1)
         call MPI_ISEND(sbufr,jsr,mreal,kr,itg(3),lgrp,msid(7),ierr)
         jsl = idimp*ncll(3,mx1)
         call MPI_ISEND(sbufl,jsl,mreal,kl,itg(4),lgrp,msid(8),ierr)
         call MPI_WAIT(msid(3),istatus,ierr)
         call MPI_WAIT(msid(4),istatus,ierr)
      endif
! make sure sbufr, sbufl, ncll, and nclr have been sent
      if (nvp /= 1) then
         do i = 1, 4
            call MPI_WAIT(msid(i+4),istatus,ierr)
         enddo
      endif
      end subroutine
!
      end module
!
! Make functions callable by Fortran77
!
!-----------------------------------------------------------------------
      subroutine PPINIT2(idproc,nvp)
      use mpplib2, only: SUB => PPINIT2
      implicit none
      integer, intent(inout) :: idproc, nvp
      call SUB(idproc,nvp)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPEXIT
      use mpplib2, only: SUB => PPEXIT
      implicit none
      call SUB()
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPABORT
      use mpplib2, only: SUB => PPABORT
      implicit none
      call SUB()
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PWTIMERA(icntrl,time,dtime)
      use mpplib2, only: SUB => PWTIMERA
      implicit none
      integer, intent(in) :: icntrl
      real, intent(inout) :: time
      double precision, intent(inout) :: dtime
      call SUB(icntrl,time,dtime)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPSUM(f,g,nxp)
      use mpplib2, only: SUB => PPSUM
      implicit none
      integer, intent(in) :: nxp
      real, dimension(nxp), intent(inout) :: f, g
      call SUB(f,g,nxp)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPDSUM(f,g,nxp)
      use mpplib2, only: SUB => PPDSUM
      implicit none
      integer, intent(in) :: nxp
      double precision, dimension(nxp), intent(inout) :: f, g
      call SUB(f,g,nxp)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPIMAX(if,ig,nxp)
      use mpplib2, only: SUB => PPIMAX
      implicit none
      integer, intent(in) :: nxp
      integer, dimension(nxp), intent(inout) :: if, ig
      call SUB(if,ig,nxp)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPDMAX(f,g,nxp)
      use mpplib2, only: SUB => PPDMAX
      implicit none
      integer, intent(in) :: nxp
      double precision, dimension(nxp), intent(inout) :: f, g
      call SUB(f,g,nxp)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNCGUARD2L(f,nyp,kstrt,nvp,nxv,nypmx)
      use mpplib2, only: SUB => PPNCGUARD2L
      implicit none
      integer, intent(in) :: nyp, kstrt, nvp, nxv, nypmx
      real, dimension(nxv,nypmx), intent(inout) :: f
      call SUB(f,nyp,kstrt,nvp,nxv,nypmx)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNAGUARD2L(f,scr,nyp,nx,kstrt,nvp,nxv,nypmx)
      use mpplib2, only: SUB => PPNAGUARD2L
      implicit none
      integer, intent(in) :: nyp, kstrt, nvp, nx, nxv, nypmx
      real, dimension(nxv,nypmx), intent(inout) :: f
      real, dimension(nxv), intent(inout) :: scr
      call SUB(f,scr,nyp,nx,kstrt,nvp,nxv,nypmx)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNACGUARD2L(f,scr,nyp,nx,ndim,kstrt,nvp,nxv,nypmx)
      use mpplib2, only: SUB => PPNACGUARD2L
      implicit none
      integer, intent(in) :: nyp, kstrt, nvp, nx, ndim, nxv, nypmx
      real, dimension(ndim,nxv,nypmx), intent(inout) :: f
      real, dimension(ndim,nxv), intent(inout) :: scr
      call SUB(f,scr,nyp,nx,ndim,kstrt,nvp,nxv,nypmx)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPTPOSE(f,g,s,t,nx,ny,kxp,kyp,kstrt,nvp,nxv,nyv,kxpd,  &
     &kypd)
      use mpplib2, only: SUB => PPTPOSE
      implicit none
      integer, intent(in) :: nx, ny, kxp, kyp, kstrt, nvp, nxv, nyv
      integer, intent(in) :: kxpd, kypd
      complex, dimension(nxv,kypd), intent(in) :: f
      complex, dimension(nyv,kxpd), intent(inout) :: g
      complex, dimension(kxp*kyp), intent(inout) :: s, t
      call SUB(f,g,s,t,nx,ny,kxp,kyp,kstrt,nvp,nxv,nyv,kxpd,kypd)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPNTPOSE(f,g,s,t,nx,ny,kxp,kyp,kstrt,nvp,ndim,nxv,nyv, &
     &kxpd,kypd)
      use mpplib2, only: SUB => PPNTPOSE
      implicit none
      integer, intent(in) :: nx, ny, kxp, kyp, kstrt, nvp, ndim
      integer, intent(in) :: nxv, nyv, kxpd, kypd
      complex, dimension(ndim,nxv,kypd), intent(in) :: f
      complex, dimension(ndim,nyv,kxpd), intent(inout) :: g
      complex, dimension(ndim,kxp*kyp), intent(inout) :: s, t
      call SUB(f,g,s,t,nx,ny,kxp,kyp,kstrt,nvp,ndim,nxv,nyv,kxpd,kypd)
      end subroutine
!
!-----------------------------------------------------------------------
      subroutine PPPMOVE2(sbufr,sbufl,rbufr,rbufl,ncll,nclr,mcll,mclr,  &
     &kstrt,nvp,idimp,nbmax,mx1)
      use mpplib2, only: SUB => PPPMOVE2
      implicit none
      integer, intent(in) :: kstrt, nvp, idimp, nbmax, mx1
      real, dimension(idimp,nbmax), intent(in) :: sbufl, sbufr
      real, dimension(idimp,nbmax), intent(inout) :: rbufl, rbufr
      integer, dimension(3,mx1), intent(in) :: ncll, nclr
      integer, dimension(3,mx1), intent(inout) :: mcll, mclr
      call SUB(sbufr,sbufl,rbufr,rbufl,ncll,nclr,mcll,mclr,kstrt,nvp,   &
     &idimp,nbmax,mx1)
      end subroutine

max_stars_repo_name: gcasabona/cuda

max_issues_repo_name: gcasabona/cuda

max_forks_repo_name: gcasabona/cuda


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 insert_radar(i4time,cldcv,cld_hts
     1         ,temp_3d,temp_sfc_k,td_sfc_k                          ! I
     1         ,grid_spacing_m,ni,nj,nk,kcloud                       ! I
     1         ,cloud_base,ref_base                                  ! I
!    1         ,tb8_k,sh_3d                                          ! I
     1         ,topo,solar_alt,r_missing_data                        ! I
     1         ,grid_ra_ref,dbz_max_2d                               ! I/O
     1         ,vis_radar_thresh_dbz                                 ! I
     1         ,l_unresolved                                         ! O
     1         ,heights_3d                                           ! I
     1         ,istatus)                                             ! O

!       Insert radar data into cloud grid and reconcile cloud/radar fields.

!       For wideband radar, "resolved" echoes can be added in as cloud if the
!       echo top is >2000m AGL and the echo is located above an existing
!       "nominal" cloud base. The nominal cloud base is taken (if present) from
!       the pre-existing cloud field, either at the horizontal grid location
!       in question or from a neighboring grid point that has existing cloud.

!       High echo tops with consistent cloud bases are treated as "resolved".
!       Low echo tops are treated as "unresolved" as are echoes that lack
!       pre-existing cloud.

!       For all "unresolved" echoes no cloud is added and the echo is blanked
!       out.

!       THIS PARAGRAPH SHOULD BE REVIEWED TO ENSURE IT IS CURRENT
!       Narrowband radar is treated similarly, except that the echo top
!       threshold does not apply. Instead, this subroutine will be actively
!       adding clouds (between the pre-existing cloud layers) only for
!       "trusted" narrowband radars such as NOWRAD. For other narrowband
!       radars, the radar is blanked out if there is no pre-existing cloud
!       in the calling routine 'laps_cloud_sub'.

!       Further reconciling of radar and satellite/cloud information is done
!       later on in subroutine 'insert_vis'.

!       1992 - 2003        Steve Albers

        use mem_namelist, ONLY: echotop_thr_a

        real temp_3d(ni,nj,nk)
        real heights_3d(ni,nj,nk)
        real temp_sfc_k(ni,nj)
        real td_sfc_k(ni,nj)

        real cldcv(ni,nj,kcloud)
        real cld_hts(kcloud)

        real grid_ra_ref(ni,nj,nk)
        real dbz_max_2d(ni,nj)
        real dbz_max_thr(ni,nj)
        real cloud_base(ni,nj)
        real cloud_base_buf(ni,nj) ! Lowest SAO/IR base within search radius
        real topo(ni,nj)
        real solar_alt(ni,nj)

        real echo_top(ni,nj)           ! L
        real echo_top_agl(ni,nj)       ! L
        real echo_top_temp(ni,nj)      ! L
        real echo_agl_thr(ni,nj)       ! L
        real lcl_agl(ni,nj)            ! L

        integer idebug_a(ni,nj)        ! L

!       Cloud not filled in unless radar echo is higher than base calculated
!       with THIS threshold.
        real     thresh_cvr
        parameter (thresh_cvr = 0.10)

        logical l_below_base, l_inserted
        logical l_resolved(ni,nj)
        logical l_unresolved(ni,nj)

!       Calculate Cloud Bases
        unlimited_ceiling = 200000.

!       Echo top threshold
!       echo_agl_thr = 2000.

!       Search radius for cloud bases
        search_radius_m = 60000.

        call get_ref_base_useable(ref_base_useable,istatus)
        if(istatus .ne. 1)return

        itest = min(13,ni)

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

!           Approximation based on a 4 deg F dewpoint depression for each
!           1000 ft of cloud base above the ground
            lcl_agl(i,j) = (temp_sfc_k(i,j) - td_sfc_k(i,j)) * 137.16 ! meters

            if(i .eq. (i/5)*5 .and. j .eq. (j/5)*5)then
                idebug_a(i,j) = 1
            else
                idebug_a(i,j) = 0
            endif

            if(solar_alt(i,j) .lt. 15. .and.
     1          temp_sfc_k(i,j) .gt. echotop_thr_a(3))then ! warmer evenings
                echo_agl_thr(i,j) = echotop_thr_a(2)
                dbz_max_thr(i,j) = 20.
            else                                           ! colder evenings
                echo_agl_thr(i,j) = echotop_thr_a(1)
                dbz_max_thr(i,j) = ref_base_useable
            endif

            cloud_base(i,j) = unlimited_ceiling
            echo_top(i,j) = r_missing_data
            echo_top_agl(i,j) = r_missing_data

            do k = kcloud-1,1,-1
                if(cldcv(i,j,k  ) .lt. thresh_cvr .and.
     1             cldcv(i,j,k+1) .ge. thresh_cvr       )then
                    cloud_base(i,j) = 0.5 * (cld_hts(k) + cld_hts(k+1))
                endif
            enddo ! k

            cloud_base_buf(i,j) = cloud_base(i,j)

            l_resolved(i,j) = .true.

        enddo ! i
        enddo ! j

        icount_below = 0
        isearch_base = 0
        insert_count_tot = 0
        iwrite_inserted = 0

        isearch_radius = nint(search_radius_m / grid_spacing_m)
        intvl = max((isearch_radius / 4),1)

        write(6,*)' isearch_radius/intvl = ',isearch_radius,intvl

        kp1 = nk

        do k = nk,1,-1 ! Essential that this go downward to detect radar tops
                       ! in time to search for a new cloud base

            icount_radar_lvl = 0
            insert_count_lvl = 0

            klow = 10000
            khigh = -1

           !Find range of cloud grid levels for this LAPS grid level
           !Find lowest level
            do kk = 1,kcloud
c               write(6,*)k,kk,height_to_zcoord(cld_hts(kk),istatus)
                if(nint(height_to_zcoord(cld_hts(kk),istatus)) .eq. k
     1                                                        )then
c                   write(6,*)' klow = ',kk
                    klow = kk
                    goto400
                endif
            enddo
400         do kk = kcloud,1,-1
                if(nint(height_to_zcoord(cld_hts(kk),istatus)) .eq. k
     1                                                        )then
c                   write(6,*)' khigh = ',kk
                    khigh = kk
                    goto500
                endif
            enddo

            if(istatus .ne. 1)then
                write(6,*)' Error in insert_radar'
                return
            endif

            if(klow .eq. 10000 .or. khigh .eq. -1)then
                write(6,*)' Skip looping, no cloud levels'
                goto600
            endif

500         write(6,*)' Looping through level ',k,klow,khigh

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

                l_inserted = .false.

                if(temp_3d(i,j,k) .lt. 278.)then
                    ref_thresh = ref_base_useable
                else
                    ref_thresh = ref_base_useable ! Lowered from 20. and 13. 9/94
!                   ref_thresh = 20.              ! Raised back to 13. in 2012
                endif

!               if(idebug_a(i,j) .eq. 1 .and. k .eq. 29)then
!               if(idebug_a(i,j) .eq. 1)then
!                   write(6,*)'etop 1:',i,j,k
!    1                    ,grid_ra_ref(i,j,k),grid_ra_ref(i,j,kp1)
!               endif

                if(grid_ra_ref(i,j,k) .ge. ref_thresh)then
                    icount_radar_lvl = icount_radar_lvl + 1
                    l_below_base = .false.

!                   if(idebug_a(i,j) .eq. 1)then
!                       write(6,*)'etop 2:',i,j,k
!    1                        ,grid_ra_ref(i,j,k),grid_ra_ref(i,j,kp1)
!                   endif

!                   Test if we are at echo top
                    if(k .eq. nk   .or.
     1                 grid_ra_ref(i,j,kp1) .lt. ref_thresh)then

                        echo_top(i,j) = heights_3d(i,j,k)
                        echo_top_agl(i,j) = echo_top(i,j) - topo(i,j)
                        echo_top_temp(i,j) = temp_3d(i,j,k)

!                       if(idebug_a(i,j) .eq. 1)then
!                           write(6,*)'Echo top found:',i,j,k
!    1                                                  ,echo_top(i,j)
!                       endif

!                       Test if we are below the cloud base
                        if(echo_top(i,j) .lt. cloud_base(i,j))then

!                           Radar Echo Top is below analyzed cloud base
!                           Various tests attempting to resolve...
!                           Search for Modified Cloud Base, i.e. other neighboring
!                           bases lower than the current buffer value

                            ilow =  max(i-isearch_radius,1)
                            ihigh = min(i+isearch_radius,ni)
                            jlow =  max(j-isearch_radius,1)
                            jhigh = min(j+isearch_radius,nj)

                            do jj = jlow,jhigh,intvl
                            do ii = ilow,ihigh,intvl
                                base_buf_lcl = max(cloud_base(ii,jj)
     1                                            ,lcl_agl(i,j))
                                cloud_base_buf(i,j)
     1                    = min(base_buf_lcl,cloud_base_buf(i,j))
                            enddo ! ii
                            enddo ! jj

                            if(cloud_base_buf(i,j) .lt. echo_top(i,j)
     1                                       .AND.
     1                         echo_top_agl(i,j) .gt. echo_agl_thr(i,j)
     1                                       .AND.
     1                         echo_top_temp(i,j) .le. 278.15   ! Cold Precip
     1                                       .AND.
     1                         dbz_max_2d(i,j) .ge. dbz_max_thr(i,j)
     1                                                             )then

                              isearch_base = isearch_base + 1
                              if(idebug_a(i,j) .eq. 1)then
                                write(6,71)i,j,k
     1                                    ,nint(echo_top(i,j))
     1                                    ,nint(cloud_base(i,j))
     1                                    ,nint(cloud_base_buf(i,j))
71                              format(' Rdr Top > Bse ',2i4,i3,3i7
     1                                ,' Resolved')
                              endif

!                             l_resolved is still TRUE

                            else ! Unresolved base
                              l_resolved(i,j) = .false.

                            endif ! Resolved Base tests (nearby base, high echo
                                  !                      top, and cold)

                        endif ! Below Cloud Base (try to resolve)

                    endif ! At Echo Top

!                   Loop through range of cloud grid levels for this LAPS level
                    do kk = klow,khigh
!                       Insert radar if we are above cloud base
                        if(cld_hts(kk) .gt. cloud_base_buf(i,j)
     1                                .and. (l_resolved(i,j).eqv..true.)
     1                                                           )then
                            cldcv(i,j,kk) = 1.0
                            insert_count_lvl = insert_count_lvl + 1
                            insert_count_tot = insert_count_tot + 1
                            l_inserted = .true.
                        else ! Radar Echo below cloud base
                            l_below_base = .true.
                        endif
                    enddo ! kk

                    if(l_below_base)then
                        icount_below = icount_below + 1
                        if(idebug_a(i,j) .eq. 1)then
                            write(6,81)i,j,k,nint(cld_hts(klow))
     1                                  ,nint(cloud_base_buf(i,j))
81                          format(' Rdr     < Bse ',2i4,i3,2i7)
                        endif
                    endif ! below base

                endif ! Reflectivity > thresh

                if( (l_inserted .eqv. .true.) .AND.
     1              (idebug_a(i,j) .eq. 1)          )then
                    write(6,591)i,j,k,l_resolved(i,j)
     1                         ,dbz_max_2d(i,j),echo_top(i,j)
     1                         ,echo_top_temp(i,j)
591                 format(' Inserted radar',2i4,i3,l2
     1                                      ,1x,f8.1,f8.0,f8.2)
                    if(echo_top(i,j) .eq. r_missing_data)then
                        write(6,*)grid_ra_ref(i,j,:)
                    endif
                    iwrite_inserted = iwrite_inserted + 1
                endif

            enddo ! j
            enddo ! i

            write(6,592)k,klow,khigh
     1          ,icount_radar_lvl,insert_count_lvl,insert_count_tot
592         format(' Inserted radar',3i3,3i8)

            kp1 = k

600         continue

        enddo ! k (LAPS grid level)

        write(6,*)' Total cloud grid points modified by radar = '
     1                                          ,insert_count_tot

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

            if(echo_top(i,j) .ne. r_missing_data)then
                if(echo_top_agl(i,j) .lt. echo_agl_thr(i,j) .and.
     1             (l_resolved(i,j).eqv..true.)       )then

!                   Should we set these to unresolved here or better yet above?
                    if(echo_top_temp(i,j) .ge. 278.15)then
                      iwarm_resolved = 1
                      l_resolved(i,j) = .false.
                    else
                      iwarm_resolved = 0
                    endif

                    if(idebug_a(i,j) .eq. 1)then
                      write(6,610)i,j
     1                        ,nint(echo_top_agl(i,j))
     1                        ,nint(echo_agl_thr(i,j))
     1                        ,nint(echo_top(i,j))
     1                        ,nint(cloud_base(i,j))
     1                        ,nint(cloud_base_buf(i,j))
     1                        ,echo_top_temp(i,j),iwarm_resolved
 610                  format('CLD_RDR - Low echo top yet resolved '
     1                    ,3i5,' ET',2i5,' Base',2i7,f7.1,i3)
                    endif
                endif
            endif

            if(l_resolved(i,j) .eqv. .false.)then

!               Reconcile radar and satellite
!               (Block out "unresolved" radar echoes)

                if(dbz_max_2d(i,j) .lt. vis_radar_thresh_dbz)then
!                   Blank out radar
                    if(idebug_a(i,j) .eq. 1)then
                      write(6,601)i,j,dbz_max_2d(i,j)
     1                         ,nint(echo_top_agl(i,j))
     1                         ,nint(echo_agl_thr(i,j))
     1                         ,echo_top_temp(i,j)
601                   format('CLD_RDR - insert_radar: '
     1                         ,'Blank out radar < '
     1                         ,2x,2i4,f6.1,' dbz',i6,1x,i6,' agl',f8.1)
                    endif

                    if(.true.)then ! Block out the radar
                        do k = 1,nk
                            grid_ra_ref(i,j,k) = ref_base
                        enddo ! k
                        dbz_max_2d(i,j) = ref_base
                    endif

                else

!                   This situation is undesirable because we want to
!                   avoid cases where the radar reflectivity is high
!                   but is "unresolved" (e.g. analyzed (SAO/IR) clouds at
!                   this point in) the analysis. We should check as to why
!                   it is unresolved (e.g. no clouds analyzed). For now we
!                   blank out the radar, but we could also create a cloud at
!                   the radar echo location.

                    if(idebug_a(i,j) .eq. 1)then
                      write(6,602)i,j,dbz_max_2d(i,j)
     1                         ,nint(echo_top_agl(i,j))
     1                         ,nint(echo_agl_thr(i,j))
     1                         ,echo_top_temp(i,j)
602                   format('CLD_RDR - insert_radar: '
     1                         ,'Blank out radar > *'
     1                         ,1x,2i4,f6.1,' dbz',i6,1x,i6,' agl',f8.1)
                    endif

                    if(.true.)then ! Block out the radar
                        do k = 1,nk
                            grid_ra_ref(i,j,k) = ref_base
                        enddo ! k
                        dbz_max_2d(i,j) = ref_base
                    endif

                endif

            else
                if(idebug_a(i,j) .eq. 1)then
                    if(echo_top(i,j) .ne. r_missing_data)then
                        write(6,612)i,j,nint(cloud_base(i,j))
     1                             ,nint(echo_top(i,j))
612                     format('Resolved ',2i4,1x,2i6)
                    endif
                endif

            endif ! l_resolved

            if(echo_top(i,j) .ne. r_missing_data)then
                if(echo_top_temp(i,j) .ge. 278.15  .AND.
     1             dbz_max_2d(i,j)    .ge. ref_base_useable)then
                    if(idebug_a(i,j) .eq. 1)then
                        write(6,*)' Warm echo remaining at',i,j
     1                           ,dbz_max_2d(i,j),echo_top_agl(i,j)
     1                           ,echo_top_temp(i,j)
                    endif
                endif
            endif

        enddo ! j
        enddo ! i

!       Blank out echoes less than "useable"
        where (grid_ra_ref .lt. ref_base_useable)grid_ra_ref = ref_base
        where (dbz_max_2d  .lt. ref_base_useable)dbz_max_2d  = ref_base

999     l_unresolved = .not. (l_resolved)

        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 bench_inline
use gpfunc, only: build_K, kern_sqexp
implicit none

integer, parameter :: nx = 4096
integer, parameter :: nd = 4

real(8), allocatable :: x(:, :), y(:), l(:)
real(8), allocatable :: K(:, :)
integer :: kx, kd
integer(8) :: count, count_rate, count_max
real(8) :: tic, toc

integer :: info

external :: dpotrf
external :: dtrmv

allocate(K(nx, nx), x(nd, nx), l(nd), y(nx))

l = 1d0

do kx = 1, nx
    call random_number(y(kx))
    do kd = 1, nd
        call random_number(x(kd, kx))
    end do
end do

!-------------------------------------------------------------------------------
! Row-wise kernel
K = 0d0
call system_clock(count, count_rate, count_max)
tic = count*1d3/count_rate

call build_K(nd, nx, nx, x, x, l, K, kern_sqexp)

call system_clock(count, count_rate, count_max)
toc = count*1d3/count_rate
print *, 'Time build_K_vec:', toc - tic, 'ms'
print *, K(5,5), K(nx,nx-5)
print *, ''

! Nugget regularization
do kx = 1, nx
    K(kx, kx) = K(kx, kx) + 1d-10
end do

call system_clock(count, count_rate, count_max)
tic = count*1d3/count_rate

! This doesn't touch anything above diagonal
print *, K(1,1:3)
print *, K(2,1:3)
print *, K(3,1:3)
call dpotrf('L', nx, K, nx, info)
print *, '------------------'
print *, K(1,1:3)
print *, K(2,1:3)
print *, K(3,1:3)

call system_clock(count, count_rate, count_max)

toc = count*1d3/count_rate
print *, 'Result Cholesky:', info
print *, 'Time Cholesky:', toc - tic, 'ms'
print *, ''

call system_clock(count, count_rate, count_max)
tic = count*1d3/count_rate

call dtrmv('L', 'N', 'N', nx, K, nx, y, 1)

call system_clock(count, count_rate, count_max)
toc = count*1d3/count_rate
print *, 'Time matmul:', toc - tic, 'ms'

deallocate(K, x, y, l)

end program bench_inline

max_stars_repo_name: krystophny/profit

max_issues_repo_name: krystophny/profit

max_forks_repo_name: krystophny/profit


text:
!==========================================================================
pure subroutine gsw_geo_strf_dyn_height_pc (sa, ct, delta_p, &
                                            geo_strf_dyn_height_pc, p_mid)
!==========================================================================
!
!  Calculates dynamic height anomaly as the integral of specific volume
!  anomaly from the the sea surface pressure (0 Pa) to the pressure p.
!  This function, gsw_geo_strf_dyn_height_pc, is to used when the
!  Absolute Salinity and Conservative Temperature are piecewise constant in
!  the vertical over sucessive pressure intervals of delta_p (such as in
!  a forward "z-coordinate" ocean model).  "geo_strf_dyn_height_pc" is
!  the dynamic height anomaly with respect to the sea surface.  That is,
!  "geo_strf_dyn_height_pc" is the geostrophic streamfunction for the
!  difference between the horizontal velocity at the pressure concerned, p,
!  and the horizontal velocity at the sea surface.  Dynamic height anomaly
!  is the geostrophic streamfunction in an isobaric surface.  The reference
!  values used for the specific volume anomaly are SA = SSO = 35.16504 g/kg
!  and CT = 0 deg C.  The output values of geo_strf_dyn_height_pc are
!  given at the mid-point pressures, p_mid, of each layer in which SA and
!  CT are vertically piecewice constant (pc).  This function calculates
!  enthalpy using the computationally-efficient 75-term expression for
!  specific volume of Roquet et al., (2015).
!
!  SA       =  Absolute Salinity                                   [ g/kg ]
!  CT       =  Conservative Temperature (ITS-90)                  [ deg C ]
!  delta_p  =  difference in sea pressure between the deep and     [ dbar ]
!              shallow extents of each layer in which SA and CT
!              are vertically constant. delta_p must be positive.
!
!  Note. sea pressure is absolute pressure minus 10.1325 dbar.
!
!  geo_strf_dyn_height_pc =  dynamic height anomaly             [ m^2/s^2 ]
!  p_mid                  =  mid-point pressure in each layer      [ dbar ]
!--------------------------------------------------------------------------

use gsw_mod_toolbox, only : gsw_enthalpy_diff, gsw_enthalpy_sso_0

use gsw_mod_error_functions, only : gsw_error_code

use gsw_mod_kinds

implicit none

real (r8), intent(in) :: sa(:), ct(:), delta_p(:)
real (r8), intent(out) :: geo_strf_dyn_height_pc(:), p_mid(:)

integer :: i, np
real (r8), allocatable :: delta_h(:), delta_h_half(:), dyn_height_deep(:)
real (r8), allocatable :: p_deep(:), p_shallow(:)

character (*), parameter :: func_name = "gsw_geo_strf_dyn_height_pc"

if (any(delta_p .lt. 0.0_r8)) then
    geo_strf_dyn_height_pc = gsw_error_code(1,func_name)
    p_mid = geo_strf_dyn_height_pc
    return
end if

np = size(delta_p)
allocate (delta_h(np),delta_h_half(np),dyn_height_deep(np))
allocate (p_deep(np),p_shallow(np))

p_deep(1) = delta_p(1)
do i = 2, np
   p_deep(i) = p_deep(i-1) + delta_p(i)
end do
p_shallow = p_deep - delta_p

delta_h = gsw_enthalpy_diff(sa,ct,p_shallow,p_deep)

dyn_height_deep(1) = -delta_h(1)
do i = 2, np
   dyn_height_deep(i) = dyn_height_deep(i-1) - delta_h(i)
end do
!            This is Phi minus Phi_0 of Eqn. (3.32.2) of IOC et al. (2010).

p_mid = 0.5_r8*(p_shallow  + p_deep)
delta_h_half = gsw_enthalpy_diff(sa,ct,p_mid,p_deep)

geo_strf_dyn_height_pc = gsw_enthalpy_sso_0(p_mid) + &
                           dyn_height_deep + delta_h_half

return
end subroutine

!--------------------------------------------------------------------------

max_stars_repo_name: minsukji/ci-debug

max_issues_repo_name: minsukji/ci-debug

max_forks_repo_name: minsukji/ci-debug


text:
      subroutine pgmast(isp1,kx,irtc)
      use tfstk
      use ffs
      use tffitcode
      use ffs_pointer
      implicit none
      integer*8 kx,kax
      integer*4 isp1,irtc,itfmessage
      integer*4 narg

      narg=isp-isp1
      if(narg .gt. 1) go to 9010
      irtc=0
      kax=ktavaloc(-1,nlat)
      call pgmast1(rlist(kax+1:kax+nlat),icomp(1:nlat),
     $     rlist(ifpos:ifpos+nlat-1))
      kx=ktflist+kax
      return
 9010 irtc=itfmessage(9,'General::narg','"1"')
      return
      end
c
      subroutine pgmast1(rmaster,icomp,pos)
      use tfstk
      use ffs
      use tffitcode
      use ffs_pointer, only:idtypec
      implicit none
      integer*4 icomp(nlat)
      real*8 pos(nlat),rmaster(nlat)
      integer*4 i,ie,id,j
      do 10 i=1,nlat
        rmaster(i)=0.d0
 10   continue
      ie=0
      do 1010 i=1,nlat-1
        if(i .le. ie)then
          go to 1010
        endif
        id=idtypec(i)
c        if(id .eq. 20)then
c          if(rmaster(i) .eq. 0.d0)then
c            if(rlist(latt(2,i)+8) .ne. 0.d0)then
c              do 1030 j=i+1,nlat-1
c                if(idtype(latt(1,j)) .eq. 20)then
c                  if(rlist(latt(2,j)+8) .ne. 0.d0)then
c                    rmaster(i)=j
c                    rmaster(j)=-i
c                    ie=i
c                    write(*,*)'b ',i,j,rmaster(i),rmaster(j)
c                    go to 1010
c                  endif
c                endif
c 1030         continue
c              call termes(6,'?Missing end of solenoid',
c     1             pname(latt(1,i)))
c              write(*,*)'a ',i,rmaster(i)
c              ie=i
c            endif
c          endif
c          go to 1010
c        endif
c     if(id .gt. 8 .and. id .ne. 31 .and. id .eq. 22
c     1     .or. id .eq. 1)then
        if(id .gt. 8 .and. id .ne. 31 .and. id .ne. 22
     $       .and. id .ne. icSOL
     1       .or. id .eq. 1)then
          go to 1010
        endif
        ie=i
        do 1020 j=i+1,nlat-1
          if(icomp(j) .eq. icomp(i))then
            ie=j
            rmaster(j)=-i
          elseif(pos(j) .eq. pos(j+1))then
          else
            go to 1021
          endif
 1020   continue
 1021   rmaster(i)=ie
 1010 continue
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD