text:
c
c     ****************************************************************
c     *                                                              *
c     *                      subroutine incontact                    *
c     *                                                              *
c     *                       written by : ag                        *
c     *                                                              *
c     *                   last modified : 07/3/2017 rhd              *
c     *                                                              *
c     *        This routine processes the user input which describes *
c     *        the contact surfaces.                                 *
c     *                                                              *
c     ****************************************************************
c
      subroutine incontact(sbflg1,sbflg2)
      use global_data ! old common.main
      use contact
      implicit integer (a-z)
c
c
      logical :: numd, numi, sbflg1, sbflg2, point_set, norm_set,
     &          matchs, outside, inside
      character :: dums*1
      double precision ::
     &     zero, dot, stiff, fric, point(3,3),
     &     tol_val, rate_val(3), depth, dumd, radius, length,
     &     ten_billion, norm(3)
      real dumr
      data zero, tol_val, ten_billion /0.0, 0.0000001, 1.0e10 /
c
c           read in new line
c
 10   continue
      call readsc
 20   continue
c
      if ( matchs('shape',4)) goto 100
      if ( matchs('surface',4)) goto 100
      if ( matchs('clear',5)) goto 1000
      if ( matchs('dump',4)) goto 2000
c
      sbflg1         = .true.
      sbflg2         = .true.
c
      goto 9999
c
c          input for contact shape definition
c
 100  continue
c
c             read shape number
c
      if ( .not. numi(shape)) then
         call errmsg (294, dumi, dums, dumr, dumd)
         goto 10
      endif
c
      if ( shape .lt. 1 .or. shape .gt. maxcontact) then
         call errmsg (295, maxcontact, dums, dumr, dumd)
         goto 10
      endif
c
c            initialize temporary shape variables
c
      point_num = 0
      stiff = zero
      fric = zero
      depth = zero
      radius = zero
      length = zero
      point_set = .false.
      norm_set = .false.
      outside = .true.
      inside = .false.
      do i=1, 3
         rate_val(i) = zero
         do j=1, 3
            point(i,j) = zero
         enddo
      enddo
c
c            read shape type
c
      if ( matchs ('plane',4)) goto 200
      if ( matchs ('cylinder',3)) goto 300
      if ( matchs ('sphere',3)) goto 400
c
      call errmsg (296, dumi, dums, dumr, dumd)
      goto 10
c
c
c          Case 1: input for contact plane definition
c          -------
c
 200  continue
c
      call readsc
c
      if ( matchs('point',5) ) goto 210
      if ( matchs('contact',4) ) goto 200
      if ( matchs('stiffness',5) ) goto 220
      if ( matchs('friction',4) ) goto 230
      if ( matchs('rate',4) ) goto 240
      if ( matchs('velocity',3) ) goto 240
      if ( matchs('depth',5) ) goto 250
c
c            if we don't match any keywords, we assume that the
c            plane definition is complete.
c
      goto 290
c
c        Input point for plane definition
c
 210  continue
      point_num = point_num + 1
      if ( point_num .gt. 3) then
         call errmsg (297, dumi, dums, dumr, dumd)
         goto 200
      endif
c
      do i = 1, 3
         if ( .not. numd (point(point_num,i))) then
            call errmsg (298, dumi, dums, dumr, dumd)
            point_num = point_num - 1
            goto 200
         endif
      enddo
      goto 200
c
c        Input plane stiffness
c
 220  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(stiff))
     &     call errmsg (299, dumi, dums, dumr, dumd)
      goto 200
c
c        Input plane friction coeffficient
c
 230  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(fric))
     &     call errmsg (300, dumi, dums, dumr, dumd)
      goto 200
c
c        Input rate of plane motion
c
 240  continue
      if (matchs('=',1)) call splunj
      do i = 1, 3
         if ( .not. numd (rate_val(i))) then
            call errmsg (301, dumi, dums, dumr, dumd)
            rate_val(1:3) = zero
            goto 200
         endif
      enddo
      goto 200
c
c        Input depth of plane
c
 250  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(depth))
     &     call errmsg (302, dumi, dums, dumr, dumd)
      goto 200
c
c          save this data as part of contact plane
c
 290  continue
c
c             if we have too few points, then error and skip this plane.
c             also, if stiffness is negative or zero, then skip this
c             plane.
c
      if ( point_num .lt. 3 ) then
         call errmsg (303, dumi, dums, dumr, dumd)
         goto 20
      else if ( stiff .le. zero ) then
         call errmsg (304, dumi, dums, dumr, dumd)
         goto 20
      endif
c
c             compute v1 and v2, edge vectors of rectangle.  Check to make
c             sure the dot product is zero -- if not, produce error.
c
      dot = zero
      do i = 1, 3
         cshape_pnt(i,shape) = point(1,i)
         cplane_vec(i,1,shape) = point(2,i) - point(1,i)
         cplane_vec(i,2,shape) = point(3,i) - point(1,i)
         dot = dot + cplane_vec(i,1,shape) * cplane_vec(i,2,shape)
      enddo
      if ( abs(dot) .gt. tol_val) then
         call errmsg (305, dumi, dums, dumr, dumd)
         goto 20
      endif
c
c             compute v3 -- normal vector
c
      call cross_prod ( cplane_vec(1,1,shape),cplane_vec(1,2,shape),
     &     cshape_norm(1,shape) )
      call normalize (cshape_norm(1,shape), dumd)
c
c             store plane constats -- stiffness, friction, rate
c
      contact_shape(shape) = 1
      contact_stiff(shape) = stiff
      contact_fric(shape) = fric
      contact_outside(shape) = .true.
      if ( depth .ne. zero ) contact_depth(shape) = depth
      do i=1, 3
         cshape_rate(i,shape) = rate_val(i)
      enddo
      use_contact = .true.
c
      goto 20
c
c          Case 2: input for contact cylinder definition
c          -------
c
 300  continue
c
      call readsc
c
      if ( matchs('point',5) ) goto 310
      if ( matchs('contact',4) ) goto 300
      if ( matchs('center',4) ) goto 300
      if ( matchs('stiffness',5) ) goto 320
      if ( matchs('friction',4) ) goto 330
      if ( matchs('rate',4) ) goto 340
      if ( matchs('velocity',3) ) goto 340
      if ( matchs('radius',3) ) goto 350
      if ( matchs('length',3) ) goto 360
      if ( matchs('direction',3) ) goto 370
      if ( matchs('outside',3) ) go to  375
      if ( matchs('inside',3) ) go to 380
c
c            if we don't match any keywords, we assume that the
c            cylinder definition is complete.
c
      goto 390
c
c        Input point on line along center of cylinder
c
 310  continue
c
      point_set = .true.
      do i = 1, 3
         if ( .not. numd (point(i,1))) then
            call errmsg (298, dumi, dums, dumr, dumd)
            point_set = .false.
            goto 300
         endif
      enddo
      goto 300
c
c        Input cylinder stiffness
c
 320  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(stiff))
     &     call errmsg (299, dumi, dums, dumr, dumd)
      goto 300
c
c        Input cylinder friction coeffficient
c
 330  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(fric))
     &     call errmsg (300, dumi, dums, dumr, dumd)
      goto 300
c
c        Input rate of motion of cylinder -- only translation is allowed
c
 340  continue
      if (matchs('=',1)) call splunj
      do i = 1, 3
         if ( .not. numd (rate_val(i))) then
            call errmsg (301, dumi, dums, dumr, dumd)
            rate_val(1:3) = zero
            goto 300
         endif
      enddo
      goto 300
c
c        Input radius of cylinder
c
 350  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(radius))
     &     call errmsg (306, dumi, dums, dumr, dumd)
      goto 300
c
c        Input length of cylinder
c
 360  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(length))
     &     call errmsg (307, dumi, dums, dumr, dumd)
      goto 300
c
c        Input direction of center line of cylinder
c
 370  continue
      norm_set = .true.
      do i = 1, 3
         if ( .not. numd (norm(i))) then
            call errmsg (308, dumi, dums, dumr, dumd)
            norm_set = .false.
            goto 300
         endif
      enddo
      goto 300
c
c        enforce contact on outside (default) of cylinder or
c        inside
c
 375  continue
      outside = .true.
      inside = .false.
      go to 300
 380  continue
      inside = .true.
      outside = .false.
      go to 300
c
c          save this data as part of cylinder
c
 390  continue
c
c             if we don't have a valid point or normal, then skip this
c             shape.
c             also, if stiffness is negative or zero, then skip this
c             plane.
c
      if (.not. point_set .or. .not. norm_set ) then
         call errmsg (309, dumi, dums, dumr, dumd)
         goto 20
      else if ( stiff .le. zero ) then
         call errmsg (304, dumi, dums, dumr, dumd)
         goto 20
      else if ( radius .le. zero ) then
         call errmsg (306, dumi, dums, dumr, dumd)
         goto 20
      else if ( length .le. zero ) then
         call errmsg (307, dumi, dums, dumr, dumd)
         goto 20
      endif
c
c             store cylinder constants
c
      cshape_pnt (1:3,shape) = point(1:3,1)
      cshape_norm (1:3,shape) = norm(1:3)
      call normalize (cshape_norm (1:3,shape), dumd)
      cshape_param (1,shape) = radius
      cshape_param (2,shape) = length
      contact_shape(shape) = 2
      contact_stiff(shape) = stiff
      contact_fric(shape) = fric
      contact_outside(shape) = outside
      if( depth .ne. zero ) contact_depth(shape) = depth
      cshape_rate(1:3,shape) = rate_val(1:3)
      use_contact = .true.
c
      goto 20
c
c          Case 3: input for contact sphere definition
c          -------
c
 400  continue
c
      call readsc
c
      if ( matchs('point',5) ) goto 410
      if ( matchs('contact',4) ) goto 400
      if ( matchs('center',4) ) goto 410
      if ( matchs('stiffness',5) ) goto 420
      if ( matchs('friction',4) ) goto 430
      if ( matchs('rate',4) ) goto 440
      if ( matchs('velocity',3) ) goto 440
      if ( matchs('radius',3) ) goto 450
      if ( matchs('outside',3) ) go to  460
      if ( matchs('inside',3) ) go to 465
c
c            if we don't match any keywords, we assume that the
c            sphere definition is complete.
c
      goto 490
c
c        Input point on line along center of sphere
c
 410  continue
c
      point_set = .true.
      do i = 1, 3
         if ( .not. numd (point(i,1))) then
            call errmsg (297, dumi, dums, dumr, dumd)
            point_set = .false.
            goto 400
         endif
      enddo
      goto 400
c
c        Input sphere stiffness
c
 420  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(stiff))
     &     call errmsg (299, dumi, dums, dumr, dumd)
      goto 400
c
c        Input sphere friction coeffficient
c
 430  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(fric))
     &     call errmsg (300, dumi, dums, dumr, dumd)
      goto 400
c
c        Input rate of motion of sphere -- only translation is allowed
c
 440  continue
      if (matchs('=',1)) call splunj
      do i = 1, 3
         if ( .not. numd (rate_val(i))) then
            call errmsg (301, dumi, dums, dumr, dumd)
            rate_val(1:3) = zero
            goto 400
         endif
      enddo
      goto 400
c
c        Input radius of sphere
c
 450  continue
      if (matchs('=',1)) call splunj
      if ( .not. numd(radius))
     &     call errmsg (306, dumi, dums, dumr, dumd)
      goto 400
c
c        enforce contact on outside (default) of sphere or
c        inside
c
 460  continue
      outside = .true.
      inside = .false.
      go to 400
 465  continue
      inside = .true.
      outside = .false.
      go to 400
c
c          save this data as part of contact sphere
c
 490  continue
c
c             if we don't have a valid point or normal, then skip this
c             shape.
c             also, if stiffness is negative or zero, then skip this
c             sphere.
c
      if (.not. point_set )then
         call errmsg (309, dumi, dums, dumr, dumd)
         goto 20
      else if ( stiff .le. zero ) then
         call errmsg (304, dumi, dums, dumr, dumd)
         goto 20
      else if ( radius .le. zero ) then
         call errmsg (306, dumi, dums, dumr, dumd)
         goto 20
      endif
c
c             store sphere constants
c
      cshape_pnt (1:3,shape) = point(1:3,1)
      cshape_param (1,shape) = radius
      contact_shape(shape) = 3
      contact_stiff(shape) = stiff
      contact_fric(shape) = fric
      cshape_rate(1:3,shape) = rate_val(1:3)
      contact_outside(shape) = outside
      use_contact = .true.
c
      goto 20
c
c          keyword "clear" was found.  zero out the contact plane
c          definitions.
c
 1000 continue
      write (out,*) '>>> Clearing contact surfaces'
c
      call contact_remove(.true.)
c
      do i=1, maxcontact
         cplane_vec(1:3,1:2,i) = zero
         cshape_norm(1:3,i) = zero
         cshape_pnt(1:3,i) = zero
         cshape_rate(1:3,i) = zero
         cshape_param(1:maxcntprm,i) = zero
         contact_stiff(i) = zero
         contact_fric(i) = zero
         contact_force(i) = zero
         contact_shape(i) = 0
         contact_depth(i) = ten_billion
         contact_outside(i) = .true.
      end do
c
      use_contact = .false.
c
      goto 10
c
c
c          keyword "dump" was found.  zero out the contact plane
c          definitions.
c
 2000 continue
      write (out,*) '>>> Dumping contact surfaces...'
      do shape = 1, maxcontact
c
         if ( contact_shape(shape) .eq. 1 ) then
c
            write (out,*) ' --> Shape ',shape,' is a sub-plane'
            write (out,*) '    corner point:'
            write (out,'(7x,3e14.6)')(cshape_pnt(i,shape),i=1,3)
            write (out,*) '    edge vectors:'
            write (out,'(7x,3e14.6)')((cplane_vec(i,j,shape),i=1,3),
     &              j=1,2)
            write (out,*) '    normal vector:'
            write (out,'(7x,3e14.6)')(cshape_norm(i,shape),i=1,3)
            write (out,*) '    velocity:'
            write (out,'(7x,3e14.6)')(cshape_rate(i,shape),i=1,3)
            write (out,*) '    stiffness:', contact_stiff(shape)
            write (out,*) '    friction:', contact_fric(shape)
            write (out,*) '    depth:', contact_depth(shape)
c
         else if ( contact_shape(shape) .eq. 2 ) then
c
            write (out,*) ' --> Shape ',shape,' is a cylinder'
            write (out,*) '    center point:'
            write (out,'(7x,3e14.6)')(cshape_pnt(i,shape),i=1,3)
            write (out,*) '    center line direction:'
            write (out,'(7x,3e14.6)')(cshape_norm(i,shape),i=1,3)
            write (out,*) '    velocity:'
            write (out,'(7x,3e14.6)')(cshape_rate(i,shape),i=1,3)
            write (out,*) '    radius:', cshape_param(1,shape)
            write (out,*) '    length:', cshape_param(2,shape)
            write (out,*) '    stiffness:', contact_stiff(shape)
            write (out,*) '    friction:', contact_fric(shape)
            write (out,*) '    outside: ', contact_outside(shape)
c
         else if ( contact_shape(shape) .eq. 3 ) then
c
            write (out,*) ' --> Shape ',shape,' is a sphere'
            write (out,*) '    center point:'
            write (out,'(7x,3e14.6)')(cshape_pnt(i,shape),i=1,3)
            write (out,*) '    velocity:'
            write (out,'(7x,3e14.6)')(cshape_rate(i,shape),i=1,3)
            write (out,*) '    radius:', cshape_param(1,shape)
            write (out,*) '    stiffness:', contact_stiff(shape)
            write (out,*) '    friction:', contact_fric(shape)
            write (out,*) '    outside: ', contact_outside(shape)
c
         endif
c
      enddo
c
      goto 10
c
 9999 continue
      return
      end
c


text:
      program prog
      do i=1,2
        if (i .eq. 2) then
          stop
        endif
        write(6, '(a,i2)') 'iteration', i
      enddo
      end

ext: f


text:
      subroutine wrisub ( lu )

      include 'data.inc'
      integer     i, lu, iitem

c     ITEMS with defined ITEMGR are substances

          do i = 1,nitem
	    if (itemgr(i).ne.'') then
	        if (itemwk(i).eq.'x') then
      	        write (lu,1000)
     j        trim(itemid(i)),trim(itemnm(i)),trim(itemun(i))
	        else
      	        write (lu,1010)
     j        trim(itemid(i)),trim(itemnm(i)),trim(itemun(i))
	        endif
	    endif
	    enddo

c     ITEMS with undefined ITEMGR and defined ITEMDE are input

          do i = 1,nitem
	    if (itemgr(i).eq.''.and.itemde(i).ne.-999.) then
    	        write (lu,1020)
     j  trim(itemid(i)),trim(itemnm(i)),trim(itemun(i)),itemde(i)
	    endif
	    enddo

	    do i = 1,noutp
	        call zoek (outpit(i),nitem,itemid,10,iitem)
              if (iitem.le.0) stop 'Bug 17-01-2011'
	        write (lu,1030) trim(outpit(i)), trim(itemnm(iitem))
	    enddo

	    write (lu,1040)
	    write (lu,1041) (trim(procid(i)),trim(procnm(i)),i=1,nproc)
	    write (lu,1043)

 1000 format ('substance ''',a,''' active'/
     j        '   description        ''',a,''''/
     j        '   concentration-unit ''',a,''''/
     j        '   waste-load-unit    ''-'''/
     j        'end-substance')
 1010 format ('substance ''',a,''' inactive'/
     j        '   description        ''',a,''''/
     j        '   concentration-unit ''',a,''''/
     j        '   waste-load-unit    ''-'''/
     j        'end-substance')
 1020 format ('parameter ''',a,''''/
     j        '   description   ''',a,''''/
     j        '   unit          ''',a,''''/
     j        '   value          ',g15.7/
     j        'end-parameter')
 1030 format ('output ''',a,''''/
     j        '   description   ''',a,''''/
     j        'end-output')
 1040 format ('active-processes')
 1041 format ('   name  ''',a,''' ''',a,'''')
 1043 format ('end-active-processes')

      return
      end

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
C
C	$Id: fngbkd.f,v 1.5 2008-07-27 00:59:02 haley Exp $
C
C                Copyright (C)  2000
C        University Corporation for Atmospheric Research
C                All Rights Reserved
C
C The use of this Software is governed by a License Agreement.
C
      SUBROUTINE FNGBKD
C
C Calling this do-nothing subroutine forces "ld" to load the following
C block data routine (but only if they are in the same ".f" file).
C
        RETURN
C
      END
CNOSPLIT
      BLOCKDATA FNGBKDX
C
      PARAMETER ( NUM=71 , NUMA=468 )
      COMMON /SUBNMS/SUBNAM,ALLNAM
      CHARACTER*10 SUBNAM(NUM),ALLNAM(NUMA)
C
      DATA (SUBNAM(I),I=1,NUM)/
     -     'AXES'  ,'CONOP1','CONOP2','CONOP3','CONOP4',
     -     'CURVE' ,'CURVE3','CURVED','DASHLN','DASHD' ,
     -     'ENCODE','FENCE3','FL2INT','FLASH1','FLASH2','FLASH3',
     -     'FLASH4','FLUSH' ,'FLUSHB','FRAME' ,'FRST3' ,'FRSTD' ,
     -     'FRSTPT','GETCHR','GETCND','GETOPT','GETSET','GETSI' ,
     -     'INTT'  ,'JUSTFY','LINE'  ,'LINE3' ,'LINED' ,'LOC'   ,
     -     'MXMY'  ,'OPTION','OPTN'  ,'PACKUM','PERIM3','PERROR',
     -     'PLOTIT','POINT' ,'POINT3','POINTS','PORGN' ,'PREOUT',
     -     'PSCALE','PSYM'  ,'PSYM3' ,'PUT42' ,'PUTINS','PWRIT' ,
     -     'PWRITX','PWRITY','PWRM'  ,'PWRT'  ,'PWRX'  ,'PWRY'  ,
     -     'SET'   ,'SET3'  ,'SETCHR','SETCND','SETI'  ,'TICK3' ,
     -     'TICK43','TRANS' ,'ULIBER','VECT3' ,'VECTD' ,'VECTOR',
     -     'WRITEB'/
C
      DATA (ALLNAM(LL),LL=  1,108)  /
     -     'ADDCHR','AGAXIS','AGBACK','AGBNCH','AGCHAX','AGCHCU',
     -     'AGCHIL','AGCHNL','AGCTCS','AGCTKO','AGCURV','AGDASH',
     -     'AGDFLT','AGDLCH','AGDSHN','AGEXAX','AGEXUS','AGEZSU',
     -     'AGFPBN','AGFTOL','AGGETC','AGGETF','AGGETI','AGGETP',
     -     'AGGTCH','AGINIT','AGKURV','AGLBLS','AGMAXI','AGMINI',
     -     'AGNUMB','AGPPID','AGPWRT','AGQURV','AGRPCH','AGRSTR',
     -     'AGSAVE','AGSCAN','AGSETC','AGSETF','AGSETI','AGSETP',
     -     'AGSRCH','AGSTCH','AGSTUP','AGUTOL','ANOTAT','ARINT' ,
     -     'AXES'  ,'BD1'   ,'BD2'   ,'BLANK' ,'BOUND' ,'BUFF'  ,
     -     'C2W'   ,'CA2B'  ,'CALCNT','CARDIN','CCHECK','CCON'  ,
     -     'CERMSG','CFVLD' ,'CH8'   ,'CHKCYC','CLGEN' ,'CLSET' ,
     -     'CMP2'  ,'COMBIN','CONBDN','CONBND','CONCAL','CONCLD',
     -     'CONCLS','CONCOM','CONDET','CONDRW','CONDSD','CONECD',
     -     'CONGEN','CONINT','CONLCM','CONLIN','CONLOC','CONLOD',
     -     'CONOP1','CONOP2','CONOP3','CONOP4','CONOT2','CONOUT',
     -     'CONPDV','CONPMM','CONPMS','CONRAN','CONRAQ','CONRAS',
     -     'CONREC','CONREO','CONSET','CONSLD','CONSSD','CONSTP',
     -     'CONSWV','CONTLK','CONTNG','CONTOR','CONXCH','COPY'   /
      DATA (ALLNAM(LL),LL=109,216)  /
     -     'CREB15','CREBIN','CTCELL','CTRINT','CURVE' ,'CURVE3',
     -     'CURVED','CURVEW','CUTUP' ,'DANDR' ,'DASHBD','DASHD' ,
     -     'DASHLN','DASHMC','DAT1'  ,'DAT2'  ,'DAT3'  ,'DAT4'  ,
     -     'DAT5'  ,'DAT6'  ,'DAT7'  ,'DAT8'  ,'DCHECK','DCON'  ,
     -     'DGMOD' ,'DISPLA','DOTEDG','DPORT' ,'DRAWI' ,'DRAWMC',
     -     'DRAWPV','DRAWS' ,'DRAWT' ,'DRCNTR','CRDRLN','DRWBRB',
     -     'DRWSTR','DRWVEC','DTOA'  ,'E9RIN' ,'ENCD'  ,'ENCODE',
     -     'ENTSR' ,'EPRIN' ,'ERROF' ,'ETOA'  ,'EXPAND','EZBARB',
     -     'EZCNTR','EZHFTN','EZISOS','EZMXY' ,'EZMY'  ,'EZSRFC',
     -     'EZSTRM','EZVEC' ,'EZXY'  ,'EZY'   ,'FDUM'  ,'FDVDLD',
     -     'FENCE3','FILL'  ,'FILLBD','FILLCK','FILLIN','FILLNC',
     -     'FILLOP','FILLPA','FL2INT','FLASH1','FLASH2','FLASH3',
     -     'FLASH4','FLMHDR','FLSALL','FLSH'  ,'FLUSH' ,'FLUSHB',
     -     'FRADV' ,'FRAME' ,'FRST3' ,'FRSTC' ,'FRSTD' ,'FRSTPT',
     -     'FRSTS' ,'FRSTW' ,'FTOI'  ,'G16UPK','GAP'   ,'GASCCH',
     -     'GENCHR','GET16' ,'GETCHR','GETCND','GETDAT','GETFLT',
     -     'GETHOL','GETINT','GETOPT','GETRNG','GETSET','GETSI' ,
     -     'GETTYP','GNEWPT','GRAY'  ,'GRDINT','GRID'  ,'GRIDAL' /
      DATA (ALLNAM(LL),LL=217,324)  /
     -     'GRIDL' ,'GTCH'  ,'GTNUM' ,'GTNUMB','HAFTON','HALFAX',
     -     'HEADER','HFINIT','HILO'  ,'I8SAV' ,'IDICTL','IDIOT' ,
     -     'IMPCHK','INIT3D','INITZI','INITZS','INITZT','INSTP' ,
     -     'INTERT','INTT'  ,'ISOSRB','ISOSRF','JLM2'  ,'JUSTFY',
     -     'KURV1S','KURV2S','LABEL' ,'LABMOD','LASTD' ,'LGFILE',
     -     'LINE'  ,'LINE3' ,'LINE3W','LINEAR','LINED' ,'LINEW' ,
     -     'LOC'   ,'MAPBD' ,'MAPCEM','MAPCHI','MAPDRW','MAPEOS',
     -     'MAPFST','MAPGRD','MAPGTC','MAPGTI','MAPGTL','MAPGTR',
     -     'MAPINT','MAPIO' ,'MAPIQ' ,'MAPIT' ,'MAPLBL','MAPLMB',
     -     'MAPLOT','MAPPOS','MAPROJ','MAPRS' ,'MAPRST','MAPSAV',
     -     'MAPSET','MAPSTC','MAPSTI','MAPSTL','MAPSTR','MAPTRE',
     -     'MAPTRN','MAPTRP','MAPUSR','MAPVEC','MAPVP' ,'MARKL' ,
     -     'MAXMIN','MCSCAL','MCTR'  ,'MCTRAN','MCTRHX','MFCLOS',
     -     'MFCMDL','MFCPFR','MFED'  ,'MFEWT' ,'MFFSTF','MFGATH',
     -     'MFGETD','MFINIT','MFNREC','MFOPEN','MFPAR1','MFPAR2',
     -     'MFPAR4','MFPAR5','MFPAR6','MFPAR7','MFPARS','MFPERR',
     -     'MFREAD','MFRMBL','MFSORT','MFTYPE','MFWRIT','MINMAX',
     -     'MKMSK' ,'MMASK' ,'MOVCHR','MTITLE','MXMY'  ,'MXMYW'  /
      DATA (ALLNAM(LL),LL=325,432)  /
     -     'NAMEFR','NERRO' ,'NEWPEN','NUFRAM','OPTION','OPTN'  ,
     -     'OPTN2' ,'OUTZRO','PACKUM','PASS2' ,'PCRBIN','PERIM' ,
     -     'PERIM3','PERIML','PERROR','PLOT'  ,'PLOTIT','PLOTIW',
     -     'PLOTMC','PLOTS' ,'PLTGRD','PNUM'  ,'POINT' ,'POINT3',
     -     'POINTS','POINTW','PORGN' ,'POSITN','PRCARG','PREFIX',
     -     'PREOUT','PREP'  ,'PRGO'  ,'PRINT' ,'PRSET' ,'PRSIMX',
     -     'PRTOD' ,'PSCALE','PSTART','PSYM'  ,'PSYM3' ,'PUT42' ,
     -     'PUTINS','PWRIQ' ,'PWRIT' ,'PWRITV','PWRITW','PWRITX',
     -     'PWRITY','PWRM'  ,'PWRT'  ,'PWRX'  ,'PWRXBD','PWRY'  ,
     -     'PWRYBD','PWRYGT','PWRYMC','PWRYSO','PWRYW' ,'PWRZ'  ,
     -     'PWRZGI','PWRZGS','PWRZGT','PWRZI' ,'PWRZOI','PWRZOS',
     -     'PWRZOT','PWRZS' ,'PWRZT' ,'PWWBKD','PWWGET','PWWSOR',
     -     'QUICK' ,'RDAMRC','RDLINE','RDSYS' ,'RDWT'  ,'REMOVE',
     -     'REORD' ,'RESET' ,'RETSR' ,'RITE'  ,'SCALEB','SCAN'  ,
     -     'SCPLTW','SCPWTW','SCRLBD','SCROLL','SET'   ,'SET3'  ,
     -     'SET3D' ,'SETCHR','SETCND','SETER' ,'SETI'  ,'SETR'  ,
     -     'SETUM' ,'SETW'  ,'SKFILE','SKIP'  ,'SPACE' ,'SPLINT',
     -     'SRFABD','SRFACE','STCNTR','STLINE','STRMLN','SUPCON' /
      DATA (ALLNAM(LL),LL=433,468)  /
     -     'SUPMAP','T20'   ,'T268'  ,'T400'  ,'TERM'  ,'TEXT'  ,
     -     'THREBD','TICK3' ,'TICK4' ,'TICK43','TICKS' ,'TPGPHC',
     -     'TR32'  ,'TRANS' ,'TRIDIP','TRN32I','TRN32S','TRN32T',
     -     'UERRBD','ULIBER','UNHEX' ,'VECT3' ,'VECTD' ,'VECTOR',
     -     'VECTW' ,'VELDAT','VELVCT','VELVEC','VMAP'  ,'W2C'   ,
     -     'WINSYS','WNDBRB','WNDOUT','WRITEB','XTCH'  ,'XTCHDT' /
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
program standard_streams
  ! Example program of standard input and output streams.
  ! Listing 6.6 in the book.
  use iso_fortran_env, only: stdin => input_unit, &
                             stdout => output_unit, &
                             stderr => error_unit
  implicit none
  character(len=1000) :: text
  read(stdin, '(a)') text
  write(stdout, '(a)') trim(text)
  write(stderr, '(a)') 'This is an error message'
end program standard_streams

max_stars_repo_name: wilsonify/c-consumer

max_issues_repo_name: wilsonify/c-consumer

max_forks_repo_name: wilsonify/c-consumer


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!This file is "main.f90"
!This file contains mainroutine program
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130

Subroutine err_finalize(err_message)
  use Global_Variables
  use salmon_parallel, only: nproc_id_global, end_parallel
  use salmon_communication, only: comm_is_root
  implicit none
  character(*),intent(in) :: err_message
  if (comm_is_root(nproc_id_global)) then
    write(*,*) err_message
  endif
  call end_parallel

  stop
End Subroutine Err_finalize

subroutine arted
  use salmon_global,   only: use_ms_maxwell,use_geometry_opt,restart_option
  use control_sc,      only: tddft_sc
  use control_ms,      only: tddft_maxwell_ms
  use control_ms_raman,only: raman_maxwell_ms
  use optimization,    only: calc_opt_ground_state
  use md_ground_state, only: calc_md_ground_state

  use salmon_parallel
  use initialization
  use ground_state
  use io_gs_wfn_k
  use io_rt_wfn_k

  implicit none

  nproc_group_tdks = nproc_group_global
  nproc_id_tdks    = nproc_id_global
  nproc_size_tdks  = nproc_size_global

  call initialize

  if(restart_option == 'new')then
    select case(iflag_calc_mode)
    case(iflag_calc_mode_gs_rt)
      call calc_ground_state
      if(write_gs_wfn_k=='y') call read_write_gs_wfn_k(iflag_write)
    case(iflag_calc_mode_gs)
      call calc_ground_state
      call read_write_gs_wfn_k(iflag_write)
      if(use_geometry_opt=='y') call calc_opt_ground_state
      if(use_adiabatic_md=='y') call calc_md_ground_state
      return
    case(iflag_calc_mode_rt)
      if(use_ms_maxwell=='y' .and. read_gs_wfn_k_ms=='y' ) then
         call read_gs_wfn_k_ms_each_macro_grid
      else if(use_ms_maxwell=='y' .and. read_rt_wfn_k_ms=='y') then
         call read_write_rt_wfn_k_ms_each_macro_grid(iflag_read_rt)
      else
         call read_write_gs_wfn_k(iflag_read)
         if(read_rt_wfn_k=='y') call read_write_rt_wfn_k(iflag_read_rt)
      endif
    end select
  else if(restart_option == 'restart')then
  else
    call Err_finalize("Invalid restart_option!")
  end if

  select case(use_ms_maxwell)
  case ('y')
    if( use_potential_model == 'n'    ) call tddft_maxwell_ms
   !if( use_potential_model == 'Raman') call raman_maxwell_ms  !not open yet
    if(write_rt_wfn_k_ms=='y') call read_write_rt_wfn_k_ms_each_macro_grid(iflag_write_rt)
  case ('n')
    call tddft_sc
    if(write_rt_wfn_k=='y') call read_write_rt_wfn_k(iflag_write_rt)
  case default
    call Err_finalize("Invalid use_ms_maxwell parameter!")
  end select

  return
end subroutine arted

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
! This comment line below causes ROSE to crash

      SUBROUTINE DGEMM ( TRANSA )
      RETURN
      END

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
C ========================================================================
C                                                                        =
C  FDISPLAY -- Example fortran program showing the use of the Client     =
C  Display Library (CDL) Fortran interface for displaying images.  In    =
C  this simple program all input is prompted for on the command line.    =
C                                                                        =
C ========================================================================

	program fdisplay
	character*132  	imname
        character*132	imtdev

C 	--------------------------
C 	Initialize the CDL package
C 	--------------------------
        call getenv('IMTDEV',imtdev)
        call cfopen(imtdev, ier)
	if (ier .gt. 0) then
	    write (*,*) 'open: Error return from CDL'
	    goto 999
	endif

        write (*, "('Image Name: ', $)")
        read (5, *) imname
        write (*, "('Frame Number: ', $)")
        read (5, *) iframe
        write (*, "('Frame buffer configuration number: ', $)")
        read (5, *) ifb

C 	----------------------------------------------------------
C	If we've got a FITS format image, go ahead and display it.
C 	----------------------------------------------------------
	call cfisfits (imname, isfits)
	if (isfits .gt. 0) then
	    call cfdisplayfits (imname, iframe, ifb, 1, ier)
	    if (ier .gt. 0) then
	        write (*,*) 'displayFITS: Error return from CDL'
	        goto 999
	    endif
	else
C 	    --------------------------------------------------------
C	    We've got an IRAF format image, go ahead and display it.
C 	    --------------------------------------------------------
	    call cfisiraf (imname, isiraf)
	    if (isiraf .gt. 0) then
	        call cfdisplayiraf (imname, 1, iframe, ifb, 1, ier)
	        if (ier .gt. 0) then
	            write (*,*) 'displayIRAF: Error return from CDL'
	            goto 999
	        endif
	    else
C 	    	----------------------------------
C	    	Unrecognized image, punt and exit.
C 	    	----------------------------------
	        write (*,*) 'Unrecognized image format'
	    endif
	endif

C 	------------------
C 	Clean up and exit.
C 	------------------
999	continue
	call cfclose (ier)
	end

max_stars_repo_name: ysBach/irafdocgen

max_issues_repo_name: ysBach/irafdocgen

max_forks_repo_name: ysBach/irafdocgen


text:
      SUBROUTINE DT_PHOINI

C***********************************************************************
C Initialization PHOJET-event generator for nucleon-nucleon interact.  *
C This version dated 16.11.95 is written by S. Roesler                 *
C Last change: s.r. 21.01.01                                           *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION amp2 , ecmini , epf , etf , ONE , pp , ppf ,
     &                 ppftmp , pt , ptf , ptftmp , q2 , scpf , sigmax ,
     &                 TINY10 , ZERO
      INTEGER i , IDT_IPDGHA , ijp , ijt , irej1 , isav , k
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (TINY10=1.0D-10,ZERO=0.0D0,ONE=1.0D0)

C nucleon-nucleon event-generator
      INCLUDE 'inc/dtmodl'
C particle properties (BAMJET index convention)
      INCLUDE 'inc/dtpart'
C Lorentz-parameters of the current interaction
      INCLUDE 'inc/dtltra'
C properties of interacting particles
      INCLUDE 'inc/dtprta'
C properties of photon/lepton projectiles
      INCLUDE 'inc/dtgpro'
C emulsion treatment
      INCLUDE 'inc/dtcomp'
C VDM parameter for photon-nucleus interactions
      INCLUDE 'inc/dtvdmp'
C nuclear potential
      INCLUDE 'inc/dtnpot'
C Glauber formalism: flags and parameters for statistics
      INCLUDE 'inc/dtglgp'

C*anfe  Obsolete cascade modifications...
C* parameters for cascade calculations:
C* maximum mumber of PDF's which can be defined in phojet (limited
C* by the dimension of ipdfs in pho_setpdf)
C      PARAMETER (MAXPDF = 20)
C* PDF parametrization and number of set for the first 30 hadrons in
C* the bamjet-code list
C*   negative numbers mean that the PDF is set in phojet,
C*   zero stands for "not a hadron"
C      DIMENSION IPARPD(30),ISETPD(30)
C* PDF parametrization
C      DATA IPARPD /
C     &  -5,-5, 0, 0, 0, 0,-5,-5,-5, 0, 0, 5,-5,-5, 5, 5, 5, 5, 5, 5,
C     &   5, 5,-5, 5, 5, 0, 0, 0, 0, 0/
C* number of set
C      DATA ISETPD /
C     &  -6,-6, 0, 0, 0, 0,-3,-6,-6, 0, 0, 2,-2,-2, 2, 2, 6, 6, 2, 6,
C     &   6, 6,-2, 2, 2, 0, 0, 0, 0, 0/

C*PHOJET105a
C     COMMON /GLOCMS/ XECM,XPCM,PMASS(2),PVIRT(2),IFPAP(2),IFPAB(2)
C     PARAMETER ( MAXPRO = 16 )
C     PARAMETER ( MAXTAB = 20 )
C     COMMON /HAXSEC/ XSECTA(4,-1:MAXPRO,4,MAXTAB),XSECT(6,-1:MAXPRO),
C    &                MXSECT(0:4,-1:MAXPRO,4),ECMSH(4,MAXTAB),ISTTAB
C     CHARACTER*8 MDLNA
C     COMMON /MODELS/ MDLNA(50),ISWMDL(50),PARMDL(200),IPAMDL(100)
C     COMMON /PROCES/ IPROCE,IDNODF,IDIFR1,IDIFR2,IDDPOM,IPRON(15)
C*PHOJET110
C  current beam selection
      INCLUDE 'inc/pobeam'
C  global event kinematics and particle IDs
      INCLUDE 'inc/pogcms'
C  hard cross sections and MC selection weights
      INCLUDE 'inc/pohrcs'
C  model switches and parameters
      INCLUDE 'inc/pomdls'
C  general process information
      INCLUDE 'inc/poprcs'
C*
      DIMENSION pp(4) , pt(4)

      LOGICAL lstart
      DATA lstart/.TRUE./

      ijp = IJProj
      ijt = IJTarg
      q2 = VIRt
C lepton-projectiles: initialize real photon instead
      IF ( (ijp.EQ.3) .OR. (ijp.EQ.4) .OR. (ijp.EQ.10) .OR. (ijp.EQ.11)
     &     ) THEN
         ijp = 7
         q2 = ZERO
      END IF

C switch Reggeon off
C     IPAMDL(3)= 0
      irej1 = LPRi
      IF ( LPHoin ) CALL PHO_INIT(-1,LOUt,irej1)

      IF ( IP.EQ.1 ) THEN
         CALL PHO_SETPAR(1,IDT_IPDGHA(ijp),0,ZERO)
      ELSE
         CALL PHO_SETPAR(1,2212,0,ZERO)
      END IF
      IF ( IT.EQ.1 ) THEN
         CALL PHO_SETPAR(2,IDT_IPDGHA(ijt),0,ZERO)
      ELSE
         CALL PHO_SETPAR(2,2212,0,ZERO)
      END IF

      DO k = 1 , 4
         pp(k) = ZERO
         pt(k) = ZERO
      END DO
C get max. possible momenta of incoming particles to be used for PHOJET ini.
      ppf = ZERO
      ptf = ZERO
      scpf = 1.5D0
      IF ( UMO.GE.1.E5 ) scpf = 5.0D0
      IF ( NCOmpo.GT.0 ) THEN
         DO i = 1 , NCOmpo
            IF ( IT.GT.1 ) THEN
               CALL DT_NCLPOT(IEMuch(i),IEMuma(i),ITZ,IT,ZERO,ZERO,0)
            ELSE
               CALL DT_NCLPOT(IPZ,IP,IEMuch(i),IEMuma(i),ZERO,ZERO,0)
            END IF
            ppftmp = MAX(PFErmp(1),PFErmn(1))
            ptftmp = MAX(PFErmp(2),PFErmn(2))
            IF ( ppftmp.GT.ppf ) ppf = ppftmp
            IF ( ptftmp.GT.ptf ) ptf = ptftmp
         END DO
      ELSE
         CALL DT_NCLPOT(IPZ,IP,ITZ,IT,ZERO,ZERO,0)
         ppf = MAX(PFErmp(1),PFErmn(1))
         ptf = MAX(PFErmp(2),PFErmn(2))
      END IF
      ptf = -ptf
      ppf = scpf*ppf
      ptf = scpf*ptf
      IF ( ijp.EQ.7 ) THEN
         amp2 = SIGN(PMAss(1)**2,PMAss(1))
         pp(3) = PPCm
         pp(4) = SQRT(amp2+pp(3)**2)
      ELSE
         epf = SQRT(ppf**2+PMAss(1)**2)
         CALL DT_LTNUC(ppf,epf,pp(3),pp(4),2)
      END IF
      etf = SQRT(ptf**2+PMAss(2)**2)
      CALL DT_LTNUC(ptf,etf,pt(3),pt(4),3)
      ecmini = SQRT((pp(4)+pt(4))**2-(pp(1)+pt(1))**2-(pp(2)+pt(2))
     &         **2-(pp(3)+pt(3))**2)
      IF ( lstart ) THEN

         IF ( LPRi.GT.4 ) WRITE (LOUt,99010) IP , IPZ , scpf , ppf , pp
99010    FORMAT (
     &         ' DT_PHOINI:    PHOJET initialized for projectile A,Z = '
     &         ,I3,',',I2,/,F4.1,'xp_F(max) = ',E10.3,'  p(max) = ',
     &         4E10.3)
         IF ( NCOmpo.GT.0 ) THEN

            IF ( LPRi.GT.4 ) WRITE (LOUt,99020) scpf , ptf , pt
99020       FORMAT (
     &         ' DT_PHOINI:    PHOJET initialized for target emulsion  '
     &         ,/,F4.1,'xp_F(max) = ',E10.3,'  p(max) = ',4E10.3)
         ELSE

            IF ( LPRi.GT.4 ) WRITE (LOUt,99030) IT , ITZ , scpf , ptf ,
     &           pt
99030       FORMAT (
     &         ' DT_PHOINI:    PHOJET initialized for target     A,Z = '
     &         ,I3,',',I2,/,F4.1,'xp_F(max) = ',E10.3,'  p(max) = ',
     &         4E10.3)
         END IF

         IF ( LPRi.GT.4 ) WRITE (LOUt,99040) ecmini
99040    FORMAT (' E_cm = ',E10.3)

         IF ( LPRi.GT.4 .AND. ijp.EQ.8 ) WRITE (LOUt,99050)
99050    FORMAT (
     &         ' DT_PHOINI: warning! proton parameters used for neutron'
     &         ,' projectile')
         lstart = .FALSE.
      END IF
C switch off new diffractive cross sections at low energies for nuclei
C (temporary solution)
      IF ( (ISWmdl(30).NE.0) .AND. ((IP.GT.1) .OR. (IT.GT.1)) ) THEN

         IF ( LPRi.GT.4 ) WRITE (LOUt,'(1X,A)')
     &         ' DT_PHOINI: model-switch 30 for nuclei re-set !'
         CALL PHO_SETMDL(30,0,1)
      END IF
C
C     IF (IJP.EQ.7) THEN
C        AMP2  = SIGN(PMASS(1)**2,PMASS(1))
C        PP(3) = PPCM
C        PP(4) = SQRT(AMP2+PP(3)**2)
C     ELSE
C        PFERMX = ZERO
C        IF (IP.GT.1) PFERMX = 0.5D0
C        EFERMX = SQRT(PFERMX**2+PMASS(1)**2)
C        CALL DT_LTNUC(PFERMX,EFERMX,PP(3),PP(4),2)
C     ENDIF
C     PFERMX = ZERO
C     IF ((IT.GT.1).OR.(NCOMPO.GT.0)) PFERMX = -0.5D0
C     EFERMX = SQRT(PFERMX**2+PMASS(2)**2)
C     CALL DT_LTNUC(PFERMX,EFERMX,PT(3),PT(4),3)
C*sr 26.10.96
      isav = IPAmdl(13)
      IF ( (ISHad(2).EQ.1) .AND.
     &     ((IJProj.EQ.7) .OR. (IJProj.EQ.3) .OR. (IJProj.EQ.4) .OR.
     &     (IJProj.EQ.10) .OR. (IJProj.EQ.11)) ) IPAmdl(13) = 1
C*

C  Initialize phojet with protons and other hadrons on demand
      CALL PHO_SETPAR(1,2212,0,ZERO)
      CALL PHO_SETPAR(2,2212,0,ZERO)

      CALL PHO_EVENT(-1,pp,pt,sigmax,irej1)

C*sr 26.10.96
      IPAmdl(13) = isav
C*
C
C patch for cascade calculations:
C define parton distribution functions for other hadrons, i.e. other
C then defined already in phojet
C*anfe remove cascade mod, since the PDFs are defined in PHOJET
C       IF ((IOGLB.EQ.100).AND..FALSE.) THEN
C
C          IF (LPRI.GT.4)
C      .   WRITE(LOUT,1006)
C  1006    FORMAT(/,1X,'PHOINI: additional parton distribution functions',
C      &          ' assiged (ID,IPAR,ISET)',/)
C          NPDF = 0
C          DO 3 I=1,30
C             IF (IPARPD(I).NE.0) THEN
C                WRITE(LOUT,*) 'Cascade path overwrites PDF settings.'
C                NPDF = NPDF+1
C                IF (NPDF.GT.MAXPDF) STOP ' PHOINI: npdf > maxpdf !'
C                IF ((IPARPD(I).GT.0).AND.(ISETPD(I).GT.0)) THEN
C                   IDPDG = IDT_IPDGHA(I)
C                   IPAR  = IPARPD(I)
C                   ISET  = ISETPD(I)
C
C                   IF (LPRI.GT.4)
C      .            WRITE(LOUT,'(13X,A8,3I6)') ANAME(I),IDPDG,IPAR,ISET
C                   CALL PHO_SETPDF(IDPDG,IDUM,IPAR,ISET,0,0,-1)
C                ENDIF
C             ENDIF
C     3    CONTINUE
C       ENDIF

C     CALL PHO_PHIST(-1,SIGMAX)

      IF ( irej1.NE.0 ) THEN

         IF ( LPRi.GT.4 ) WRITE (LOUt,99060)
99060    FORMAT (1X,'PHOINI:   PHOJET event-initialization failed!')
         STOP
      END IF

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


text:
c=======================================================================
      subroutine dpmultmeta(nrec,nvar,y,sigma2e,
     &                      a0b0,m1rand,s2inv,s2invm2,
     &                      nu,s1ra nd,psiinv,
     &                      mcmc,nsave,
     &                      cpo,randsave,thetasave,
     &                      alpha,m1,s1,
     &                      ncluster,muclus,ss,ccluster,cstrt,
     &                      iflag,prob,sigma2ei,seed,
     &                      s1inv,s1invm1,
     &                      theta,workm1,workm2,workm3,
     &                      workmh1,workmh2,
     &                      workv1,workv2,ywork)
c=======================================================================
c
c     Copyright: Alejandro Jara, 2008-2010.
c
c     Version 1.0:
c
c     Last modification: 20-04-2008.
c
c      Alejandro Jara
c      Department of Statistics
c      Facultad de Matematicas
c      Pontificia Universidad Catolica de Chile
c      Casilla 306, Correo 22
c      Santiago
c      Chile
c      Voice: +56-2-3544506  URL  : http://www.mat.puc.cl/~ajara
c      Fax  : +56-2-3547729  Email: atjara@uc.cl
c
c=======================================================================

      implicit none

c+++++Data
      integer nrec,nvar
      real*8 y(nrec,nvar)
      real*8 sigma2e(nrec,nvar*(nvar+1)/2)

c+++++Prior
      integer nu,m1rand,s1rand
      real*8 aa0,ab0,a0b0(2)
      real*8 psiinv(nvar,nvar)
      real*8 s2inv(nvar,nvar),s2invm2(nvar)

c+++++MCMC parameters
      integer mcmc(3),nburn,nskip,nsave,ndisplay

c+++++Output
      real*8 cpo(nrec,2)
      real*8 randsave(nsave,(nrec+1)*nvar)
      real*8 thetasave(nsave,nvar+nvar*(nvar+1)/2+2)

c+++++Current values of the parameters
      integer ncluster,ss(nrec)
      real*8 alpha,m1(nvar),s1(nvar,nvar)
      real*8 muclus(nrec+1,nvar)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++External working space
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
      integer iflag(nvar)
      integer seed(2)
      real*8 s1inv(nvar,nvar)
      real*8 s1invm1(nvar)
      real*8 sigma2ei(nrec,nvar*(nvar+1)/2)
      real*8 theta(nvar)
      real*8 workm1(nvar,nvar),workm2(nvar,nvar),workm3(nvar,nvar)
      real*8 workmh1(nvar*(nvar+1)/2)
      real*8 workmh2(nvar*(nvar+1)/2)
      real*8 workv1(nvar),workv2(nvar)
      real*8 ywork(nvar)

c+++++DPM
      integer cstrt(nrec,nrec)
      integer ccluster(nrec)
      real*8 prob(nrec+1)

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c+++++Internal working space
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

c+++++General
      integer count,dispcount,evali
      integer i,ii,ihmssf,isave,iscan
      integer j,k,l,j1,j2
      integer ns,nscan,nuniqs,ok,sprint
      integer seed1,seed2,since,skipcount
      real*8 tmp1,tmp2

c+++++CPU time
      real*8 sec00,sec0,sec1,sec

c++++ Define parameters

      aa0=a0b0(1)
      ab0=a0b0(2)

      nburn=mcmc(1)
      nskip=mcmc(2)
      ndisplay=mcmc(3)

      nuniqs=nvar*(nvar+1)/2

c++++ set random number generator

      seed1=seed(1)
      seed2=seed(2)
      call setall(seed1,seed2)

c++++ cluster structure

      do i=1,nrec
         ccluster(ss(i))=ccluster(ss(i))+1
         cstrt(ss(i),ccluster(ss(i)))=i
         do j1=1,nvar
            do j2=1,nvar
               workm1(j1,j2)=sigma2e(i,ihmssf(j1,j2,nvar))
            end do
         end do
         call inverse(workm1,nvar,iflag)
         do j1=1,nvar
            do j2=j1,nvar
               sigma2ei(i,ihmssf(j1,j2,nvar))=workm1(j1,j2)
            end do
         end do
      end do

      do j1=1,nvar
         do j2=1,nvar
            s1inv(j1,j2)=s1(j1,j2)
         end do
      end do
      call inverse(s1inv,nvar,iflag)

      do j1=1,nvar
         tmp1=0.d0
         do j2=1,nvar
            tmp1=tmp1+s1inv(j1,j2)*m1(j2)
         end do
         s1invm1(j1)=tmp1
      end do

c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
c++++ start the MCMC algorithm
c+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++

      isave=0
      skipcount=0
      dispcount=0
      nscan=nburn+(nskip+1)*(nsave)

      call cpu_time(sec0)
      sec00=0.d0

      do iscan=1,nscan

c+++++++ check if the user has requested an interrupt
         call rchkusr()

c++++++++++++++++++++++++++++++++++
c+++++++ DP part
c++++++++++++++++++++++++++++++++++

c+++++++ a) Polya Urn based on a collapsed state

         do i=1,nrec

            ns=ccluster(ss(i))

            do j=1,nvar
               ywork(j)=y(i,j)
            end do

c++++++++++ subject in cluster with more than 1 observations

            if(ns.gt.1)then

               j=1
               ok=0
               do while(ok.eq.0.and.j.le.ns)
                  if(cstrt(ss(i),j).eq.i)ok=j
                  j=j+1
               end do

               do j=ok,ns-1
                  cstrt(ss(i),j)=cstrt(ss(i),j+1)
               end do

               ccluster(ss(i))=ccluster(ss(i))-1

               do j=1,ncluster
                  do k=1,nvar
                     do l=1,nvar
                        workm1(k,l)=s1inv(k,l)
                     end do
                     workv1(k)=s1invm1(k)
                  end do

                  do k=1,ccluster(j)
                     do j1=1,nvar
                        tmp1=0.d0
                        do j2=1,nvar
                           tmp1=tmp1+sigma2ei(cstrt(j,k),
     &                          ihmssf(j1,j2,nvar))*y(cstrt(j,k),j2)
                           workm1(j1,j2)=workm1(j1,j2)+
     &                          sigma2ei(cstrt(j,k),ihmssf(j1,j2,nvar))
                        end do
                        workv1(j1)=workv1(j1)+tmp1
                     end do
                  end do

                  call inverse(workm1,nvar,iflag)

                  do j1=1,nvar
                     tmp1=0.d0
                     do j2=1,nvar
                        tmp1=tmp1+workm1(j1,j2)*workv1(j2)
                        workm2(j1,j2)=workm1(j1,j2)+
     &                                sigma2e(i,ihmssf(j1,j2,nvar))
                     end do
                     workv2(j1)=tmp1
                  end do

                  call dmvnd(nvar,ywork,workv2,workm2,tmp2,iflag)
                  prob(j)=exp(log(dble(ccluster(j)))+tmp2)
               end do

               do j1=1,nvar
                  workv2(j1)=m1(j1)
                  do j2=1,nvar
                     workm2(j1,j2)=s1(j1,j2)+
     &                             sigma2e(i,ihmssf(j1,j2,nvar))
                  end do
               end do
               call dmvnd(nvar,ywork,workv2,workm2,tmp2,iflag)
               prob(ncluster+1)=exp(log(alpha)+tmp2)

               call simdisc(prob,nrec+1,ncluster+1,evali)

               ss(i)=evali

               ccluster(evali)=ccluster(evali)+1

               cstrt(evali,ccluster(evali))=i

               if(evali.gt.ncluster)then
                  ncluster=ncluster+1
               end if
            end if

c++++++++++ subject in cluster with only 1 observation

            if(ns.eq.1)then

               since=ss(i)

               if(since.lt.ncluster)then
                   call relabeldpm(i,since,nrec,1,ncluster,
     &                             ccluster,ss,cstrt)
               end if

               ccluster(ncluster)=ccluster(ncluster)-1
               ncluster=ncluster-1

               do j=1,ncluster
                  do k=1,nvar
                     do l=1,nvar
                        workm1(k,l)=s1inv(k,l)
                     end do
                     workv1(k)=s1invm1(k)
                  end do

                  do k=1,ccluster(j)
                     do j1=1,nvar
                        tmp1=0.d0
                        do j2=1,nvar
                           tmp1=tmp1+sigma2ei(cstrt(j,k),
     &                          ihmssf(j1,j2,nvar))*y(cstrt(j,k),j2)
                           workm1(j1,j2)=workm1(j1,j2)+
     &                          sigma2ei(cstrt(j,k),ihmssf(j1,j2,nvar))
                        end do
                        workv1(j1)=workv1(j1)+tmp1
                     end do
                  end do

                  call inverse(workm1,nvar,iflag)

                  do j1=1,nvar
                     tmp1=0.d0
                     do j2=1,nvar
                        tmp1=tmp1+workm1(j1,j2)*workv1(j2)
                        workm2(j1,j2)=workm1(j1,j2)+
     &                                sigma2e(i,ihmssf(j1,j2,nvar))
                     end do
                     workv2(j1)=tmp1
                  end do

                  call dmvnd(nvar,ywork,workv2,workm2,tmp2,iflag)
                  prob(j)=exp(log(dble(ccluster(j)))+tmp2)
               end do

               do j1=1,nvar
                  workv2(j1)=m1(j1)
                  do j2=1,nvar
                     workm2(j1,j2)=s1(j1,j2)+
     &                  sigma2e(i,ihmssf(j1,j2,nvar))
                  end do
               end do
               call dmvnd(nvar,ywork,workv2,workm2,tmp2,iflag)
               prob(ncluster+1)=exp(log(alpha)+tmp2)

               call simdisc(prob,nrec+1,ncluster+1,evali)

               ss(i)=evali

               ccluster(evali)=ccluster(evali)+1

               cstrt(evali,ccluster(evali))=i

               if(evali.gt.ncluster)then
                  ncluster=ncluster+1
               end if
            end if

         end do

c+++++++ b) Resampling step

         do ii=1,ncluster

c++++++++++ check if the user has requested an interrupt
            call rchkusr()

            ns=ccluster(ii)

            do k=1,nvar
               do l=1,nvar
                  workm1(k,l)=s1inv(k,l)
               end do
               workv1(k)=s1invm1(k)
            end do

            do k=1,ns
               do j1=1,nvar
                  tmp1=0.d0
                  do j2=1,nvar
                     tmp1=tmp1+sigma2ei(cstrt(ii,k),ihmssf(j1,j2,nvar))*
     &                    y(cstrt(ii,k),j2)
                     workm1(j1,j2)=workm1(j1,j2)+sigma2ei(cstrt(ii,k),
     &                             ihmssf(j1,j2,nvar))
                  end do
                  workv1(j1)=workv1(j1)+tmp1
               end do
            end do

            call inverse(workm1,nvar,iflag)

            do j1=1,nvar
               tmp1=0.d0
               do j2=1,nvar
                  tmp1=tmp1+workm1(j1,j2)*workv1(j2)
               end do
               workv2(j1)=tmp1
            end do

            call rmvnorm(nvar,workv2,workm1,workmh1,workv1,theta)

            do i=1,nvar
               muclus(ii,i)=theta(i)
            end do
         end do

c++++++++++++++++++++++++++++++++++
c+++++++ Base line distribution
c++++++++++++++++++++++++++++++++++

c+++++++ check if the user has requested an interrupt
         call rchkusr()

c+++++++ mean of the normal component

         if(m1rand.eq.1)then

            do i=1,nvar
               workv1(i)=s2invm2(i)
               do j=1,nvar
                  workm1(i,j)=s2inv(i,j)+dble(ncluster)*s1inv(i,j)
               end do
            end do
            call inverse(workm1,nvar,iflag)

            do i=1,ncluster
               do j=1,nvar
                  tmp1=0.d0
                  do k=1,nvar
                     tmp1=tmp1+s1inv(j,k)*muclus(i,k)
                  end do
                  workv1(j)=workv1(j)+tmp1
               end do
            end do

            do j1=1,nvar
               tmp1=0.d0
               do j2=1,nvar
                  tmp1=tmp1+workm1(j1,j2)*workv1(j2)
               end do
               workv2(j1)=tmp1
            end do

            call rmvnorm(nvar,workv2,workm1,workmh1,workv1,theta)

            do i=1,nvar
               m1(i)=theta(i)
            end do

         end if

c+++++++ scale matrix of the inverted-Wishart component

         if(s1rand.eq.1)then

            do i=1,nvar
               do j=1,nvar
                  workm1(i,j)=psiinv(i,j)
               end do
            end do

            do i=1,ncluster
               do j=1,nvar
                  do k=1,nvar
                     workm1(j,k)=workm1(j,k)+(muclus(i,j)-m1(j))*
     &                                       (muclus(i,k)-m1(k))
                  end do
               end do
            end do

            call riwishart(nvar,nu+ncluster,workm1,workm2,workm3,
     &                     workv1,workmh1,workmh2,iflag)

            do i=1,nvar
               do j=1,nvar
                  s1(i,j)=workm1(i,j)
                  s1inv(i,j)=workm2(i,j)
               end do
            end do
         end if

         do j1=1,nvar
            tmp1=0.d0
            do j2=1,nvar
               tmp1=tmp1+s1inv(j1,j2)*m1(j2)
            end do
            s1invm1(j1)=tmp1
         end do

c++++++++++++++++++++++++++++++++++
c+++++++ Precision parameter
c++++++++++++++++++++++++++++++++++

         if(aa0.gt.0.d0)then
            call samalph(alpha,aa0,ab0,ncluster,nrec)
         end if

c++++++++++++++++++++++++++++++++++
c+++++++ save samples
c++++++++++++++++++++++++++++++++++

         if(iscan.gt.nburn)then
            skipcount=skipcount+1
            if(skipcount.gt.nskip)then
               isave=isave+1
               dispcount=dispcount+1

               count=0

c+++++++++++++ normal baseline mean

               do i=1,nvar
                  count=count+1
                  thetasave(isave,count)=m1(i)
               end do

c+++++++++++++ IW baseline scale

               do i=1,nvar
                  do j=i,nvar
                     count=count+1
                     thetasave(isave,count)=s1(i,j)
                  end do
               end do

c+++++++++++++ cluster information

               count=count+1
               thetasave(isave,count)=ncluster
               count=count+1
               thetasave(isave,count)=alpha

c+++++++++++++ random effects
               count=0
               do i=1,nrec
                  do j=1,nvar
                     count=count+1
                     randsave(isave,count)=muclus(ss(i),j)
                  end do
               end do

c+++++++++++++ predictive information

               do i=1,ncluster
                  prob(i)=dble(ccluster(i))/(alpha+dble(nrec))
               end do
               prob(ncluster+1)=alpha/(alpha+dble(nrec))

               call simdisc(prob,nrec+1,ncluster+1,evali)

               if(evali.le.ncluster)then
                  do i=1,nvar
                     theta(i)=muclus(evali,i)
                  end do
                else
                  call rmvnorm(nvar,m1,s1,workmh1,workv1,theta)
               end if

               do i=1,nvar
                  count=count+1
                  randsave(isave,count)=theta(i)
               end do

c+++++++++++++ cpo and save samples

               do i=1,nrec
                  do j=1,nvar
                     ywork(j)=y(i,j)
                     workv1(j)=muclus(ss(i),j)
                     do k=1,nvar
                        workm1(j,k)=sigma2e(i,ihmssf(j,k,nvar))
                     end do
                  end do

                  call dmvnd(nvar,ywork,workv1,workm1,tmp1,iflag)

                  tmp1=exp(tmp1)
                  cpo(i,1)=cpo(i,1)+1.0d0/tmp1
                  cpo(i,2)=cpo(i,2)+tmp1
               end do

c+++++++++++++ print
               skipcount = 0
               if(dispcount.ge.ndisplay)then
                  call cpu_time(sec1)
                  sec00=sec00+(sec1-sec0)
                  sec=sec00
                  sec0=sec1
                  tmp1=sprint(isave,nsave,sec)
                  dispcount=0
               end if
            end if
         end if
      end do

      do i=1,nrec
         cpo(i,1)=dble(nsave)/cpo(i,1)
         cpo(i,2)=cpo(i,2)/dble(nsave)
      end do

      return
      end

max_stars_repo_name: kdevick/densitymediation

max_issues_repo_name: kdevick/densitymediation

max_forks_repo_name: kdevick/densitymediation


text:
      SUBROUTINE Q2BCD (EST,PLANAR,RMAT,ET,IERROR)
C
C     BASIC CALCULATIONS ARE PERFORMED FOR THE QDMEM2 ELEMENT IN THIS
C     ROUTINE (DOUBLE-PRECISION VERSION)
C
      LOGICAL          PLANAR
      REAL             EST(1)
      DOUBLE PRECISION MAG    ,D12(3) ,G1(3) ,IAREA  ,D13(3),GRID(3,5),
     1                 G2(3)  ,ITWOH  ,D24(3),VEC(3) ,G3(3) ,ET(3,3)  ,
     2                 G5(3)  ,G4(3)  ,DADOTB,RMAT(3,5)
      EQUIVALENCE      (GRID(1,1),G1(1)),(GRID(1,2),G2(1)),
     1                 (GRID(1,3),G3(1)),(GRID(1,4),G4(1)),
     2                 (GRID(1,5),G5(1))
C
C     MOVE GRID COORDINATES AND MAKE DOUBLE-PRECISION IF THIS IS THE
C     DOUBLE-PRECISION VERSION.
C
      DO 10 I = 1,3
      G1(I) = EST(I+10)
      G2(I) = EST(I+14)
      G3(I) = EST(I+18)
      G4(I) = EST(I+22)
   10 CONTINUE
C
C     FORM  D   , D   AND  D   VECTORS
C            13    24       12
C
      DO 20 I = 1,3
      D12(I) = G2(I) - G1(I)
      D13(I) = G3(I) - G1(I)
      D24(I) = G4(I) - G2(I)
   20 CONTINUE
C
C     NVEC = D13 CROSS D24 = K-VECTOR (UN-NORMALIZED)
C
      CALL DAXB (D13,D24,VEC)
      MAG   = DSQRT(DADOTB(VEC,VEC))
      IAREA = 0.5D0*MAG
C
C     NORMALIZE K-VECTOR
C
      IF (MAG) 100,100,30
   30 ET(1,3) = VEC(1)/MAG
      ET(2,3) = VEC(2)/MAG
      ET(3,3) = VEC(3)/MAG
C
C     H = .5 * (D   DOT K-VEC)
C                12
C
      ITWOH = DADOTB(D12,ET(1,3))
C
C     I-VECTOR (UN-NORMALIZED) = (D  ) - 2 H (K-VECTOR)
C                                  12
C
      DO 40 I = 1,3
      VEC(I) = D12(I) - ITWOH*ET(I,3)
   40 CONTINUE
      MAG = DSQRT(DADOTB(VEC,VEC))
C
C     NORMALIZE I-VECTOR
C
      IF (MAG) 100,100,50
   50 ET(1,1) = VEC(1)/MAG
      ET(2,1) = VEC(2)/MAG
      ET(3,1) = VEC(3)/MAG
C
C     JVEC = KVEC CROSS IVEC
C
      CALL DAXB (ET(1,3),ET(1,1),ET(1,2))
C
C     FILL THE SUB-TRIANGLE ELEMENT COORDINATE MATRIX
C
      DO 60 I = 1,3
      G5(I) = 0.25D0*(G1(I) + G2(I) + G3(I) + G4(I))
   60 CONTINUE
      RMAT(1,1) = 0.0D0
      RMAT(2,1) = 0.0D0
      RMAT(3,1) =-ITWOH/2.0D0
      DO 70 I = 2,5
      VEC(1) = GRID(1,I) - G1(1)
      VEC(2) = GRID(2,I) - G1(2)
      VEC(3) = GRID(3,I) - G1(3)
      CALL GMMATD (ET,3,3,0, VEC,3,1,0, RMAT(1,I))
      RMAT(1,I) = RMAT(1,I) + RMAT(1,1)
      RMAT(2,I) = RMAT(2,I) + RMAT(2,1)
      RMAT(3,I) = RMAT(3,I) + RMAT(3,1)
   70 CONTINUE
C
C     SET PLANAR FLAG .TRUE. OR .FALSE.
C
      IF ((ITWOH/2.0D0)**2/IAREA .LE. 0.01D0) GO TO 80
      PLANAR = .FALSE.
      GO TO 90
   80 PLANAR = .TRUE.
C
C     ALL BASIC CALCULATIONS NOW COMPLETE
C
   90 IERROR = 0
      RETURN
C
C     ERROR CONDITION, BAD ELEMENT GEOMETRY.
C
  100 IERROR = 1
      RETURN
      END


text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
MODULE scf_data
use salmon_global
use salmon_xc, only: xc_functional
use salmon_pp, only: pp_info,pp_grid
implicit none
!-------------------- Parameters
integer, parameter :: maxntmg=10

! Physical constants
real(8),parameter :: E2=14.4d0, H2M=3.81d0, a_B=0.529177d0
real(8),parameter :: Ry=13.6058d0, Pi=3.141592653589793d0
real(8),parameter :: fs2eVinv = 1.51925d0

! Coefficients of finite-difference
integer,parameter :: Nd=4       ! (2*Nd+1)-points finite-difference used in Taylor expansion
integer,parameter :: Ndh=4       ! (2*Nd+1)-points finite-difference used in Hartree routine
real(8),parameter :: cN0=-205.d0/72.d0                 ! for Laplacian
real(8),parameter :: cN1=8.d0/5.d0  , cN2=-1.d0/5.d0
real(8),parameter :: cN3=8.d0/315.d0, cN4=-1.d0/560.d0
real(8),parameter :: bN1=4.d0/5.d0  , bN2=-1.d0/5.d0   ! for Nabla
real(8),parameter :: bN3=4.d0/105.d0, bN4=-1.d0/280.d0

real(8),parameter :: cN0_Nd1=-2.d0
real(8),parameter :: cN1_Nd1=1.d0
real(8),parameter :: bN1_Nd1=1.d0/2.d0

real(8),parameter :: cN0_Nd2=-5.d0/2.d0
real(8),parameter :: cN1_Nd2=4.d0/3.d0, cN2_Nd2=-1.d0/12.d0
real(8),parameter :: bN1_Nd2=2.d0/3.d0, bN2_Nd2=-1.d0/12.d0

real(8),parameter :: cN0_Nd3=-49.d0/18.d0
real(8),parameter :: cN1_Nd3=3.d0/2.d0, cN2_Nd3=-3.d0/20.d0
real(8),parameter :: cN3_Nd3=1.d0/90.d0
real(8),parameter :: bN1_Nd3=3.d0/4.d0, bN2_Nd3=-3.d0/20.d0
real(8),parameter :: bN3_Nd3=1.d0/60.d0

real(8),parameter :: cN0_Nd5=-5269.d0/1800.d0
real(8),parameter :: cN1_Nd5=5.d0/3.d0, cN2_Nd5=-5.d0/21.d0
real(8),parameter :: cN3_Nd5=5.d0/126.d0, cN4_Nd5=-5.d0/1008.d0
real(8),parameter :: cN5_Nd5=1.d0/3150.d0

real(8),parameter :: cN0_Nd6=-5369.d0/1800.d0
real(8),parameter :: cN1_Nd6=12.d0/7.d0, cN2_Nd6=-15.d0/56.d0
real(8),parameter :: cN3_Nd6=10.d0/189.d0, cN4_Nd6=-1.d0/112.d0
real(8),parameter :: cN5_Nd6=2.d0/1925.d0, cN6_Nd6=-1.d0/16632.d0

real(8),parameter :: cN0_Nd7=-266681.d0/88200.d0
real(8),parameter :: cN1_Nd7=7.d0/4.d0, cN2_Nd7=-7.d0/24.d0
real(8),parameter :: cN3_Nd7=7.d0/108.d0, cN4_Nd7=-7.d0/528.d0
real(8),parameter :: cN5_Nd7=7.d0/3300.d0, cN6_Nd7=-7.d0/30888.d0
real(8),parameter :: cN7_Nd7=1.d0/84084.d0

real(8),parameter :: cN0_Nd8=-1077749.d0/352800.d0
real(8),parameter :: cN1_Nd8=16.d0/9.d0, cN2_Nd8=-14.d0/45.d0
real(8),parameter :: cN3_Nd8=112.d0/1485.d0, cN4_Nd8=-7.d0/396.d0
real(8),parameter :: cN5_Nd8=112.d0/32175.d0, cN6_Nd8=-2.d0/3861.d0
real(8),parameter :: cN7_Nd8=16.d0/315315.d0, cN8_Nd8=-1.d0/411840.d0

real(8),parameter :: cN0_Nd9=-9778141.d0/3175200.d0
real(8),parameter :: cN1_Nd9=9.d0/5.d0, cN2_Nd9=-18.d0/55.d0
real(8),parameter :: cN3_Nd9=14.d0/165.d0, cN4_Nd9=-63.d0/2860.d0
real(8),parameter :: cN5_Nd9=18.d0/3575.d0, cN6_Nd9=-2.d0/2145.d0
real(8),parameter :: cN7_Nd9=9.d0/70070.d0, cN8_Nd9=-9.d0/777920.d0
real(8),parameter :: cN9_Nd9=1.d0/1969110.d0

real(8),parameter :: cN0_Nd10=-1968329.d0/635040.d0
real(8),parameter :: cN1_Nd10=20.d0/11.d0, cN2_Nd10=-15.d0/44.d0
real(8),parameter :: cN3_Nd10=40.d0/429.d0, cN4_Nd10=-15.d0/572.d0
real(8),parameter :: cN5_Nd10=24.d0/3575.d0, cN6_Nd10=-5.d0/3432.d0
real(8),parameter :: cN7_Nd10=30.d0/119119.d0, cN8_Nd10=-5.d0/155584.d0
real(8),parameter :: cN9_Nd10=10.d0/3741309.d0, cN10_Nd10=-1.d0/9237800.d0

real(8),parameter :: cN0_Nd11=-239437889.d0/76839840.d0
real(8),parameter :: cN1_Nd11=11.d0/6.d0, cN2_Nd11=-55.d0/156.d0
real(8),parameter :: cN3_Nd11=55.d0/546.d0, cN4_Nd11=-11.d0/364.d0
real(8),parameter :: cN5_Nd11=11.d0/1300.d0, cN6_Nd11=-11.d0/5304.d0
real(8),parameter :: cN7_Nd11=55.d0/129948.d0, cN8_Nd11=-55.d0/806208.d0
real(8),parameter :: cN9_Nd11=11.d0/1360476.d0, cN10_Nd11=-11.d0/17635800.d0
real(8),parameter :: cN11_Nd11=1.d0/42678636.d0

real(8),parameter :: cN0_Nd12=-240505109.d0/76839840.d0
real(8),parameter :: cN1_Nd12=24.d0/13.d0, cN2_Nd12=-33.d0/91.d0
real(8),parameter :: cN3_Nd12=88.d0/819.d0, cN4_Nd12=-99.d0/2912.d0
real(8),parameter :: cN5_Nd12=396.d0/38675.d0, cN6_Nd12=-11.d0/3978.d0
real(8),parameter :: cN7_Nd12=132.d0/205751.d0, cN8_Nd12=-33.d0/268736.d0
real(8),parameter :: cN9_Nd12=44.d0/2380833.d0, cN10_Nd12=-3.d0/1469650.d0
real(8),parameter :: cN11_Nd12=12.d0/81800719.d0, cN12_Nd12=-1.d0/194699232.d0

real(8) :: cNmat(0:12,0:12),bNmat(0:12,0:12)

!-------------------- Global variables

! Exchange Correlation
type(xc_functional) :: xc_func

integer :: iflag_ps

integer :: inumcpu_check

integer,allocatable :: ob_sta_all_kgrid(:),ob_end_all_kgrid(:),iobnum_all_kgrid(:)

integer :: nproc_Mxin(3),nproc_Mxin_mul
integer :: nproc_ob_spin(2)
integer :: nproc_Mxin_s(3), nproc_Mxin_mul_s
integer :: nproc_Mxin_s_dm(3), nproc_Mxin_mul_s_dm

integer :: ista_Mx_ori(3),iend_Mx_ori(3),inum_Mx_ori(3)
integer,allocatable :: ista_Mxin(:,:),iend_Mxin(:,:),inum_Mxin(:,:)
integer,allocatable :: ista_Mxin_old(:,:),iend_Mxin_old(:,:),inum_Mxin_old(:,:)
integer,allocatable :: ista_Mxin_s(:,:),iend_Mxin_s(:,:),inum_Mxin_s(:,:)
integer :: max_lg_num

integer :: Miter       ! Total number of Iteration for SCF calculation
integer :: Miter_rt    ! Total number of Iteration for RT calculation

integer :: iflag_diisjump

integer :: iflag_writepsi
real(8),allocatable :: rto(:)
integer ,allocatable:: idip2int(:)
real(8),allocatable :: rto_ix(:,:)
real(8),allocatable :: rbox_array_dip2(:,:)
real(8),allocatable :: rbox_array2_dip2(:,:)
real(8),allocatable :: rbox_array_dip2q(:,:,:)
real(8),allocatable :: rbox_array2_dip2q(:,:,:)
real(8),allocatable :: rbox_array_dip2e(:)
real(8),allocatable :: rbox_array2_dip2e(:)

integer :: ilsda

integer :: iflag_stopt
integer :: iter_stopt
integer :: istopt
integer :: istopt_tranc

integer :: MST(2),ifMST(2),itotMST
integer :: itotfMST
integer :: MST0(2),itotMST0
integer :: Mx(3),Mxin(3),Mxin_old(3)

real(8) :: rnetot

character(8),allocatable :: AtomName(:)
integer,allocatable :: iAtomicNumber(:)
integer,allocatable :: istopt_a(:)

real(8) :: Hgs(3)        ! Grid spacing
real(8) :: Hold(3)     ! Grid spacing
real(8) :: Hvol
real(8) :: Harray(3,maxntmg)  ! Grid spacing
real(8) :: rLsize(3,maxntmg)    ! size of the box

integer :: maxMps

! Pseudopotential
type(pp_info) :: pp
type(pp_grid) :: ppg,ppg_all
integer,parameter :: Nrmax=3000,Lmax=4
integer,allocatable :: NRloc(:)
real(8),allocatable :: Rloc(:)
real(8),allocatable :: Zps(:)              ! Pseudo charge
real(8),allocatable :: Rps(:)              ! Core radius
real(8),allocatable :: Mass(:)             ! Atomic weight
integer,allocatable :: Jxyz_all(:,:,:),Mps_all(:),Jxxyyzz_all(:,:,:)
integer,allocatable :: Mps(:)
integer,allocatable :: Jxyz_tmp1(:,:,:)
integer,allocatable :: Jxyz_tmp2(:,:,:)
integer,allocatable :: Jxxyyzz_tmp1(:,:,:)
integer,allocatable :: Jxxyyzz_tmp2(:,:,:)
integer,allocatable :: Jxxyyzz(:,:,:)
integer :: Mlmps
integer :: Mlps(maxMKI),Lref(maxMKI)
real(8),allocatable :: Vpsl(:,:,:)                 ! Local pseudopotential
real(8),allocatable :: Vpsl_atom(:,:,:,:)
real(8),allocatable :: uV_all(:,:,:),uVu(:,:)          ! Non-local
!Nonlinear core correction
logical :: flag_nlcc = .false.

real(8),allocatable :: rocc(:,:)                    ! Occupation number
real(8),allocatable :: psi(:,:,:,:,:)              ! Single particle orbitals
real(8),allocatable :: psi_mesh(:,:,:,:,:)         ! Single particle orbitals
real(8),allocatable :: zpsi_mesh(:,:,:,:,:)        ! Single particle orbitals
real(8),allocatable :: psi_old(:,:,:,:,:)
complex(8),allocatable :: zpsi_old(:,:,:,:,:)

real(8),allocatable :: rho_in(:,:,:,:)
real(8),allocatable :: rho_out(:,:,:,:)
real(8),allocatable :: rho_s_in(:,:,:,:,:)
real(8),allocatable :: rho_s_out(:,:,:,:,:)

real(8),allocatable :: esp(:,:)         ! Single particle energy
real(8),allocatable :: esp2(:,:)        ! Single particle energy
real(8),allocatable :: rho(:,:,:)       ! Single particle density
real(8),allocatable :: rho0(:,:,:)      ! Single particle density
real(8),allocatable :: rho_diff(:,:,:)  ! Single particle density
real(8),allocatable :: rho_s(:,:,:,:)   ! Single particle density for each spin
real(8),allocatable :: Vh(:,:,:)        ! Hartree potential
real(8),allocatable :: Vxc(:,:,:)       ! Exchange-Correlation potential
real(8),allocatable :: Vxc_s(:,:,:,:)   ! Exchange-Correlation potential for each spin
integer :: ihpsieff
real(8),allocatable :: elf(:,:,:)
complex(8),allocatable :: zpsi(:,:,:,:,:)
complex(8),allocatable :: zpsi_in(:,:,:,:,:),zpsi_out(:,:,:,:,:)
complex(8),allocatable :: zpsi_t0(:,:,:,:,:)
complex(8),allocatable :: ttpsi(:,:,:)
integer :: iSCFRT

real(8),allocatable :: Vbox(:,:,:)

real(8),allocatable :: Ex_fast(:,:,:),Ec_fast(:,:,:)

real(8) :: elp3(3000)
real(8) :: elp5(3000)

integer, allocatable :: idiis_sd(:)

integer :: iwksta(3),iwkend(3),iwknum(3)
integer :: iwk2sta(3),iwk2end(3),iwk2num(3)
integer :: iwk3sta(3),iwk3end(3),iwk3num(3)
integer :: iwk_size

integer :: maxlm

integer :: imesh_oddeven(3)

integer :: version_num(2)

complex(8), allocatable :: zc(:)
real(8), allocatable :: Dp(:,:), Dp2(:,:,:)
real(8), allocatable :: Qp(:,:,:), Qp2(:,:,:,:)
real(8), allocatable :: rIe(:), rIe2(:,:)
real(8), allocatable :: tene(:)
real(8) :: vecDs(3)
real(8),allocatable :: vecDs2(:,:)
real(8) :: vecQs(3,3)
real(8),allocatable :: vecQs2(:,:,:)

integer :: num_pole

integer :: itotNtime

integer :: num_datafiles_OUT2

real(8),allocatable :: Gs(:,:,:),Gl(:,:,:)
complex(8),allocatable :: tx_exp(:,:),ty_exp(:,:),tz_exp(:,:)

integer :: lg_sta(3),lg_end(3),lg_num(3)
integer :: mg_sta(3),mg_end(3),mg_num(3)
integer :: ng_sta(3),ng_end(3),ng_num(3)

integer :: lg_old_sta(3),lg_old_end(3),lg_old_num(3)
integer :: mg_old_sta(3),mg_old_end(3),mg_old_num(3)
integer :: ng_old_sta(3),ng_old_end(3),ng_old_num(3)

real(8),allocatable :: gridcoo(:,:)

integer :: iobnum

integer :: k_sta,k_end,k_num

integer :: kx_hock_sta(3),kx_hock_end(3),kx_hock_num(3)

integer :: num_kpoints_3d(3)
integer :: num_kpoints_rd

real(8),allocatable :: wtk(:)

real(8),allocatable :: norm_diff_psi_stock(:,:)

real(8) :: Etot
real(8) :: Exc

integer :: imr(3),imrs(3),igroup

real(8),allocatable :: rho_stock(:,:,:,:)
real(8),allocatable :: Vlocal_stock(:,:,:,:)
integer,parameter :: num_rho_stock=21

integer :: iflag_subspace_diag

real(8),allocatable :: Vh_stock1(:,:,:)
real(8),allocatable :: Vh_stock2(:,:,:)

real(8),allocatable :: Vlocal(:,:,:,:)
real(8),allocatable :: Vlocal2(:,:,:,:)

! use for hartree routine
integer :: iterVh

real(8),allocatable :: rhobox(:,:,:)
real(8),allocatable :: rhobox_s(:,:,:,:)

integer :: lg_num_fmax(3)

real(8), allocatable :: rhobox1_all(:,:,:), rhobox2_all(:,:,:)

integer :: iDiter(maxntmg)

character(LEN=100) :: file_OUT_ini
integer :: num_mol
real(8) :: rlatcon
real(8), allocatable :: coo_mol_ini(:,:)

integer :: lg_sta_ini(3),lg_end_ini(3),lg_num_ini(3)
integer :: mg_sta_ini(3),mg_end_ini(3),mg_num_ini(3)
real(8) :: H_ini(3)
real(8) :: rLsize_ini(3)

integer :: img

integer :: itt
integer :: ikind_eext   !0:No external field, 1: dipoleApprox

character(3)  :: dir
character(2)  :: dir2

real(8) :: Fst,Fst2(2)
real(8) :: romega, romega2(2)
real(8) :: pulse_T, pulse_T2(2)
real(8) :: rlaser_I, rlaser_I2(2)
real(8) :: tau, tau2(2), delay, rcycle

character(2)  :: denplane  ! plane for writing density (xy, yz, xz)
integer       :: idensum   ! whether density is summed up along direction
                           ! perpendicular to the plane
                           ! (0: not summed, 1: summed)
real(8)       :: posplane  ! position of the plane
                           ! (only for idensum = 0)

character(1) :: circular

character(100):: rtOutFile
character(100):: rtDiffOutFile
character(100):: rtELFOutFile
character(100):: file_Projection
character(20):: fileNumber

integer,allocatable :: Jxyz(:,:,:)
real(8),allocatable :: uV(:,:,:)

integer,allocatable :: numatom_ps(:,:,:)
integer,allocatable :: iatomnum_ps(:,:,:,:)
integer :: maxMps_all

real(8), allocatable :: rho_n(:,:,:)
real(8), allocatable :: Vh_n(:,:,:)
real(8), allocatable :: Vh0(:,:,:)
complex(8), allocatable :: zpsi_n(:,:,:,:,:)

complex(8), allocatable :: Ex_static(:,:,:),Ey_static(:,:,:),Ez_static(:,:,:)

real(8),allocatable :: curr(:,:)
real(8),allocatable :: sumcurr(:,:)
real(8),allocatable :: rE_ind(:,:)

integer :: ilasbound_sta(3),ilasbound_end(3)
real(8) :: rlaser_center(3)

complex(8) :: cumnum

real(8) :: Eion

integer :: iblacsinit
integer :: CONTEXT, IAM, MYCOL, MYROW, NPCOL, NPROCS2, NPROW
integer :: DESCA( 50 ), DESCZ( 50 )

real(8) :: rho_region1(100)
real(8) :: rho_region2(100)
integer :: num_rho_region(100)
integer,allocatable :: rho_region_nx(:,:)
integer,allocatable :: rho_region_ny(:,:)
integer,allocatable :: rho_region_nz(:,:)

integer :: numspin

integer :: icalcforce
real(8),allocatable :: rforce(:,:)
integer :: iflag_md
real(8),allocatable :: dRion(:,:,:)
real(8),allocatable :: Rion_eq(:,:)
real(8),parameter :: umass=1822.9d0

integer :: wmaxMI

real(8) :: fcN(0:12)
real(8) :: fbN(0:12)

real(8),allocatable :: k_rd(:,:),ksquare(:)
real(8),allocatable :: k_rd0(:,:),ksquare0(:)

real(8),allocatable :: A_ext(:,:)
real(8),allocatable :: A_ind(:,:)
real(8),allocatable :: A_tot(:,:)
real(8),allocatable :: E_ext(:,:)
real(8),allocatable :: E_ind(:,:)
real(8),allocatable :: E_tot(:,:)

integer,allocatable :: oblist(:)

integer,allocatable::icoo1d(:,:)

integer :: MI_read

integer :: ik_oddeven

! variables for FFTE routine
integer,dimension(3) :: LNPU
integer :: NPUZ,NPUY,NPUW

real(8) :: absorption(0:100000)
real(8) :: absorption_d(0:100000)
real(8) :: absorption_id(0:100000)

integer :: iflag_dos
integer :: iflag_pdos

real(8) , allocatable :: rxk_ob(:,:,:,:),rhxk_ob(:,:,:,:),rgk_ob(:,:,:,:),rpk_ob(:,:,:,:)

complex(8) , allocatable :: zxk_ob(:,:,:,:),zhxk_ob(:,:,:,:),zgk_ob(:,:,:,:),zpk_ob(:,:,:,:)
complex(8) , allocatable :: zpko_ob(:,:,:,:),zhtpsi_ob(:,:,:,:)

integer :: iflag_ELF

integer :: iflag_indA

integer :: iflag_hartree

real(8),allocatable :: vonf_sd(:,:,:),eonf_sd(:,:,:,:)

!filename
character(100) :: file_OUT
character(100) :: file_IN
character(100) :: LDA_Info
character(100) :: file_eigen
character(100) :: file_RT
character(100) :: file_alpha_lr
character(100) :: file_alpha_pulse
character(100) :: file_RT_q
character(100) :: file_alpha_q
character(100) :: file_RT_e
character(100) :: file_RT_dip2
character(100) :: file_alpha_dip2
character(100) :: file_RT_dip2_q
character(100) :: file_alpha_dip2_q
character(100) :: file_RT_dip2_e
character(100) :: file_external
character(100) :: file_OUT_rt
character(100) :: file_IN_rt

!arrays for isend and irecv
real(8),allocatable :: srmatbox1_x_3d(:,:,:),srmatbox1_y_3d(:,:,:),srmatbox1_z_3d(:,:,:)
real(8),allocatable :: srmatbox2_x_3d(:,:,:),srmatbox2_y_3d(:,:,:),srmatbox2_z_3d(:,:,:)
real(8),allocatable :: srmatbox3_x_3d(:,:,:),srmatbox3_y_3d(:,:,:),srmatbox3_z_3d(:,:,:)
real(8),allocatable :: srmatbox4_x_3d(:,:,:),srmatbox4_y_3d(:,:,:),srmatbox4_z_3d(:,:,:)

complex(8),allocatable :: scmatbox1_x_3d(:,:,:),scmatbox1_y_3d(:,:,:),scmatbox1_z_3d(:,:,:)
complex(8),allocatable :: scmatbox2_x_3d(:,:,:),scmatbox2_y_3d(:,:,:),scmatbox2_z_3d(:,:,:)
complex(8),allocatable :: scmatbox3_x_3d(:,:,:),scmatbox3_y_3d(:,:,:),scmatbox3_z_3d(:,:,:)
complex(8),allocatable :: scmatbox4_x_3d(:,:,:),scmatbox4_y_3d(:,:,:),scmatbox4_z_3d(:,:,:)

real(8),allocatable :: srmatbox1_x_5d(:,:,:,:,:),srmatbox1_y_5d(:,:,:,:,:),srmatbox1_z_5d(:,:,:,:,:)
real(8),allocatable :: srmatbox2_x_5d(:,:,:,:,:),srmatbox2_y_5d(:,:,:,:,:),srmatbox2_z_5d(:,:,:,:,:)
real(8),allocatable :: srmatbox3_x_5d(:,:,:,:,:),srmatbox3_y_5d(:,:,:,:,:),srmatbox3_z_5d(:,:,:,:,:)
real(8),allocatable :: srmatbox4_x_5d(:,:,:,:,:),srmatbox4_y_5d(:,:,:,:,:),srmatbox4_z_5d(:,:,:,:,:)

complex(8),allocatable :: scmatbox1_x_5d(:,:,:,:,:),scmatbox1_y_5d(:,:,:,:,:),scmatbox1_z_5d(:,:,:,:,:)
complex(8),allocatable :: scmatbox2_x_5d(:,:,:,:,:),scmatbox2_y_5d(:,:,:,:,:),scmatbox2_z_5d(:,:,:,:,:)
complex(8),allocatable :: scmatbox3_x_5d(:,:,:,:,:),scmatbox3_y_5d(:,:,:,:,:),scmatbox3_z_5d(:,:,:,:,:)
complex(8),allocatable :: scmatbox4_x_5d(:,:,:,:,:),scmatbox4_y_5d(:,:,:,:,:),scmatbox4_z_5d(:,:,:,:,:)

CONTAINS

!=========================================================================
!-------------------------------------------------------------------------
!-------------------------------------------------------------------------
!=========================================================================
subroutine snum_procs
use salmon_parallel, only: nproc_size_global
implicit none
integer :: i,ibox,ipow

nproc_Mxin_mul=nproc_size_global/nproc_ob

ibox=1
ipow=0
do i=1,29
  if(ibox<nproc_Mxin_mul)then
    ipow=ipow+1
    ibox=ibox*2
  end if
end do

nproc_Mxin(1:3)=1
do i=1,ipow
  if(mod(i,3)==1)then
    nproc_Mxin(3)=nproc_Mxin(3)*2
  else if(mod(i,3)==2)then
    nproc_Mxin(2)=nproc_Mxin(2)*2
  else
    nproc_Mxin(1)=nproc_Mxin(1)*2
  end if
end do

ibox=1
ipow=0
do i=1,29
  if(ibox<nproc_size_global)then
    ipow=ipow+1
    ibox=ibox*2
  end if
end do

nproc_Mxin_s(1:3)=1
do i=1,ipow
  if(mod(i,3)==1)then
    nproc_Mxin_s(3)=nproc_Mxin_s(3)*2
  else if(mod(i,3)==2)then
    nproc_Mxin_s(2)=nproc_Mxin_s(2)*2
  else
    nproc_Mxin_s(1)=nproc_Mxin_s(1)*2
  end if
end do

end subroutine snum_procs

!======================================================================
subroutine init_mesh_s
use salmon_parallel, only: nproc_size_global, nproc_id_global
implicit none

nproc_Mxin_mul_s_dm=nproc_Mxin_s_dm(1)*nproc_Mxin_s_dm(2)*nproc_Mxin_s_dm(3)

allocate(ista_Mxin_s(3,0:nproc_size_global-1),iend_Mxin_s(3,0:nproc_size_global-1))
allocate(inum_Mxin_s(3,0:nproc_size_global-1))

call setng(ng_sta,ng_end,ng_num,ista_Mxin_s,iend_Mxin_s,inum_Mxin_s, &
           nproc_size_global,nproc_id_global,nproc_Mxin,nproc_Mxin_s_dm,ista_Mxin,iend_Mxin,isequential)

end subroutine init_mesh_s

!=====================================================================
subroutine make_iwksta_iwkend
use salmon_parallel, only: nproc_id_global
implicit none

if(iwk_size==1)then
  iwksta(1:3)=ista_Mxin(1:3,nproc_id_global)
  iwkend(1:3)=iend_Mxin(1:3,nproc_id_global)
  iwk2sta(1:3)=ista_Mxin(1:3,nproc_id_global)-Nd
  iwk2end(1:3)=iend_Mxin(1:3,nproc_id_global)+Nd
  iwk3sta(1:3)=ista_Mxin(1:3,nproc_id_global)
  iwk3end(1:3)=iend_Mxin(1:3,nproc_id_global)
else if(iwk_size==2)then
  iwksta(1:3)=ista_Mxin(1:3,nproc_id_global)-Nd
  iwkend(1:3)=iend_Mxin(1:3,nproc_id_global)+Nd
  iwk2sta(1:3)=ista_Mxin(1:3,nproc_id_global)-Nd
  iwk2end(1:3)=iend_Mxin(1:3,nproc_id_global)+Nd
  iwk3sta(1:3)=ista_Mxin(1:3,nproc_id_global)
  iwk3end(1:3)=iend_Mxin(1:3,nproc_id_global)
else if(iwk_size==11.or.iwk_size==31)then
  iwksta(1:3)=ista_Mxin_s(1:3,nproc_id_global)
  iwkend(1:3)=iend_Mxin_s(1:3,nproc_id_global)
  iwk2sta(1:3)=ista_Mxin_s(1:3,nproc_id_global)-Ndh
  iwk2end(1:3)=iend_Mxin_s(1:3,nproc_id_global)+Ndh
  iwk3sta(1:3)=ista_Mxin_s(1:3,nproc_id_global)
  iwk3end(1:3)=iend_Mxin_s(1:3,nproc_id_global)
else if(iwk_size==12.or.iwk_size==32)then
  iwksta(1:3)=ista_Mxin_s(1:3,nproc_id_global)-Ndh
  iwkend(1:3)=iend_Mxin_s(1:3,nproc_id_global)+Ndh
  iwk2sta(1:3)=ista_Mxin_s(1:3,nproc_id_global)-Ndh
  iwk2end(1:3)=iend_Mxin_s(1:3,nproc_id_global)+Ndh
  iwk3sta(1:3)=ista_Mxin_s(1:3,nproc_id_global)
  iwk3end(1:3)=iend_Mxin_s(1:3,nproc_id_global)
end if
iwknum(1:3)=iwkend(1:3)-iwksta(1:3)+1
iwk2num(1:3)=iwk2end(1:3)-iwk2sta(1:3)+1
iwk3num(1:3)=iwk3end(1:3)-iwk3sta(1:3)+1

end subroutine make_iwksta_iwkend

END MODULE scf_data


text:
C     LAST UPDATE 21/02/95
C+++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++
C
      REAL FUNCTION FALL(PALL)
      IMPLICIT NONE
C
C Purpose: This function allows the generation of profiles from
C          multiple lattices of different types.
C
C Calls:  3: FBRAG  , FCONT  , CHICAL
C
C-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-+-
C Parameters:
C
      INTEGER MAXALL,MAXLAT,MAXDAT
      PARAMETER(MAXALL=50,MAXLAT=3,MAXDAT=262144)
      INTEGER NLLPX
      PARAMETER(NLLPX=10000)
      INTEGER JWIDTHAV
      PARAMETER(JWIDTHAV=50)
      INTEGER IPOINTX
      PARAMETER(IPOINTX=JWIDTHAV*NLLPX)
      INTEGER IWIDTHAV,MPOINTX
      PARAMETER(IWIDTHAV=30,MPOINTX=IWIDTHAV*IPOINTX)
C
C Arguments:
C
      REAL PALL(MAXALL)
C
C Local variables:
C
      INTEGER I,J,NP,NPALL,IFLAG,NPAR
      REAL PR(14),SIG(NLLPX+8)
      REAL CHISQ,RFAC
C
C Common variables:
C
      REAL ALLINE(NLLPX+8),AD(MAXDAT),PSIG(MAXDAT)
      REAL ALLCEL(6,MAXLAT),LATLIM(6,MAXLAT),RADMM(MAXLAT)
      REAL SNORM(MAXLAT),SSCALE(NLLPX)
      REAL DELR,RMIN,DELZ,ZMIN
      REAL DMIN,DMAX,RADM,SN
      INTEGER*2 JMIN(NLLPX),JMAX(NLLPX)
      INTEGER*2 IMIN(IPOINTX),IMAX(IPOINTX)
      CHARACTER*1 IDCIN(MPOINTX)
      INTEGER L1,L2,NR1,NR2,NLDIV,NBESS
      INTEGER ISIZE(5),IHKL(4,NLLPX),INFO(3),MINFO(5,MAXLAT),
     &        MSIZE(5,MAXLAT),MBESS(MAXLAT),MLDIV(MAXLAT)
      INTEGER N1,N2,IBACK,NLLP,NPRD
      INTEGER MPTR,IPTR,LPTR,NPTR,NLAT
      LOGICAL SIGMA,POLAR
C
C Common blocks:
C
      COMMON /INSPEC/ DELR,RMIN,DELZ,ZMIN,POLAR
      COMMON /POINTS/ NLLP,NPRD,INFO,IHKL
      COMMON /COEFFS/ JMIN,JMAX,IMIN,IMAX,IDCIN,ISIZE
      COMMON /DLIMIT/ DMIN,DMAX,L1,L2,NR1,NR2,SN,NLDIV,NBESS,RADM
      COMMON /IMDATA/ N1,N2,IBACK,AD,PSIG,SIGMA
      COMMON /FITPAR/ ALLINE
      COMMON /POINTR/ MPTR,IPTR,LPTR,NPTR
      COMMON /LATICE/ NLAT,ALLCEL,MINFO,LATLIM,RADMM,MSIZE,MBESS,MLDIV
      COMMON /SCALES/ SNORM,SSCALE
C
C-----------------------------------------------------------------------
C
      IFLAG = 0
      NPALL = 0
      MPTR = 0
      IPTR = 0
      LPTR = 0
      NPTR = 0
      DO 40 I=1,NLAT
         WRITE(6,1000)I
         WRITE(4,1000)I
         NPRD = MINFO(4,I)
         DO 10 J=1,3
            INFO(J) = MINFO(J,I)
 10      CONTINUE
         IF(INFO(1).EQ.1)THEN
            NP = 14
         ELSE
            NP = 4
         ENDIF
         DO 20 J=1,NP
            PR(J) = PALL(NPALL+J)
 20      CONTINUE
         NPALL = NPALL + NP
         DO 30 J=1,5
            ISIZE(J) = MSIZE(J,I)
 30      CONTINUE
         NBESS = MBESS(I)
         RADM = RADMM(I)
         NLDIV = MLDIV(I)
         L1 = NINT(LATLIM(1,I))
         L2 = NINT(LATLIM(2,I))
         NR1 = NINT(LATLIM(3,I))
         NR2 = NINT(LATLIM(4,I))
         DMIN = LATLIM(5,I)
         DMAX = LATLIM(6,I)
         IF(INFO(1).EQ.0)THEN
            CALL FCONT(PR)
            MINFO(4,I) = NPRD
         ELSE
            CALL FBRAG(PR)
         ENDIF
         MINFO(5,I) = LPTR
         SNORM(I) = SN
 40   CONTINUE
      IF(MOD(IBACK/16,2).EQ.1)THEN
         NPAR = NLLP + 8
      ELSE
         NPAR = NLLP
      ENDIF
      CALL CHICAL(ALLINE,JMIN,JMAX,NLLP,NPAR,IMIN,IMAX,IPTR,IDCIN,
     &            MPTR,AD,PSIG,N1,N2,IBACK,CHISQ,RFAC,SIG,
     &            SIGMA,IFLAG)

      FALL = RFAC
      RETURN
 1000 FORMAT(/1X,'Calculating profiles for cell ',I1)
      END


text:
      subroutine muatco(xnval)
c               * angular coefficients *
c        sous programmes utilises  cwig3j
c
      implicit double precision (a-h,o-z)
      dimension xnval(30)
      common/itescf/testy,rap(2),teste,nz,norb,norbsc
      common/mulabk/afgk
      dimension afgk(30,30,0:3)
      common/ratom1/xnel(30),en(30),scc(30),scw(30),sce(30),
     1nq(30),kap(30),nmax(30)
c#mn
       external cwig3j

      do i=1,30
         do j=1,30
            do k=0,3
               afgk(i,j,k)=0.0d00
            enddo
         enddo
      enddo

      do i=1,norb
         li= abs(kap(i))*2-1
         do j=1,i
            lj= abs(kap(j))*2-1
            kmax=(li+lj)/2
            kmin= abs(li-lj)/2
            if ((kap(i)*kap(j)).lt.0) kmin=kmin+1
c calculate a_k(i,j)
            m=0
            if (j.eq.i .and. xnval(i).le.0.0d0) m=1
c           use to test SIC
c           if (j.eq.i) m=1

            afgk(j,i,0)=afgk(j,i,0)+xnel(i)*(xnel(j)-m)
            if (xnval(i).le.0.0d0 .or. xnval(j).le.0.0d0) then
c calculate b_k(i,j)
               b=afgk(j,i,0)
               if (j.eq.i .and. xnval(i).le.0.0d0) then
                  a=li
                  b=-b*(a+1.0d 00)/a
                  kmin = kmin+2
               endif
               do k = kmin, kmax,2
                  afgk(i,j,k/2)=b*(cwig3j(li,k*2,lj,1,0,2)**2)
               enddo
            endif
         enddo
      enddo
      return
      end


text:
*..............................................................

!-------------------------------------------------------------------------
!         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
!-------------------------------------------------------------------------
!BOP
!
! !ROUTINE:    ODS_GetList --- Reads a list from an ODS file
!
! !DESCRIPTION:
!    \label{ODS:GetList}
!     Reads the user-defined list or array of strings
!     (identified by the character string, ListName) from a
!     NetCDF data file (identified by the ODS file handle id).
!
!     Note: For a list of error codes, see Table~\ref{tab:errors}.
!
! !INTERFACE:
!
      subroutine ODS_GetList ( id, ListName, ListSz, List, ierr )
!
! !INPUT PARAMETERS:
      implicit         NONE
      integer          id          ! ODS file handle
      character  * (*) ListName    ! The name of the list to be
                                   !   obtained.  The case of each
                                   !   letter is significant.
!
! !INPUT/OUTPUT PARAMETERS:
      integer         ListSz       ! on input: maximum number of
                                   !   values (usually determined
                                   !   by the space allocated for
                                   !   storage)
                                   ! on output: the number of
                                   !   values obtained from the
                                   !   file
                                   ! note: If the input value is
                                   !   smaller than the number
                                   !   of values to be read, then
                                   !   the routine exits with a
                                   !   system error message and
                                   !   code of ODS_DimErr ( = -3 ).
                                   !   Also, nval is set to the
                                   !   minimum value required in
                                   !   order for the routine to
                                   !   execute successfully.
!
! !OUTPUT PARAMETERS:
      character  * (*) List  ( * ) ! The entries of the list.
      integer    ierr              ! Error code. If non-zero, an
                                   !  error has occurred. For a list
                                   !  of possible values, see the
                                   !  description section of this
                                   !  prologue.
!
!     NOTE:  Each file in the present version of ODS ( version 2.00 )
!            contains the lists with the following names
!
!            kt_names    Name of GEOS/DAS data types
!            kt_units    Units for each GEOS/DAS data type
!            kx_names    Name of GEOS/DAS data sources
!            kx_meta     kx specific metadata information.  This
!                        information is used to specify the meaning
!                        of the ODS variable "xm" or to specify the
!                        name of the external OMS file name.  This
!                        list is available only in versions 2.00
!                        or later.
!            qcx_names   information about the meaning of the ODS
!                        variable, "qcexcl".  This list is available
!                        only in versions 2.00 or later.
!
! !SEE ALSO:
!     ODS_PutList ( Put the user-defined list string or array
!                   of strings to the NetCDF data file. )
!
! !FILES USED:
!     netcdf.inc, a header file, for defining NetCDF library
!            parameters
!     ods_hdf.h, a header file, for defining hardwired constants
!            and defining global variables and setting up data
!            structures
!     ods_stdio.h, a header file, for defining standard input/output
!            unit numbers
!
! !LIBRARIES ACCESSED:
!     NetCDF
!
! !REVISION HISTORY:
!     08Mar1996   C. Redder   Original version
!     13Apr1998   C. Redder   Made routine to be an interface routine
!                             in ODS version 2.00.  Added checks made
!                             to input arguments.
!     20Apr1998   C. Redder   Fixed bug in bounds check
!     02Nov1999   C. Redder   Fixed bug in storing the list entries
!                             in the last block of data.
!     19Nov1999   C. Redder   Added a latex label in and moved the
!                             subroutine statement into the prologue.
!                             Modified the comments for the return
!                             status code.
!     06Dec1999   C. Redder   Corrections to the documentation in the
!                             prologue.
!     16Feb2000   R. Todling  Rename stdio.h to ods_stdio.h
!
!EOP
!-------------------------------------------------------------------------

      include 'netcdf.inc'
      include 'ods_hdf.h'
      include 'ods_stdio.h'

*     corners and edge lengths
*     ------------------------
      integer   corner ( MAXNCDIM )
      integer   edges  ( MAXNCDIM )

*     Other variables
*     ---------------
      character VarNam * ( MaxNCNam )
                            ! name of variable (i.e. list )
                            !   returned from NetCDF routine
                            !   ncvinq
      integer   nc_id       ! temporary storage for NetCDF
                            !   file id
      integer   varid       ! temporary storage for NetCDF
                            !   variable id
      integer   NC_VarType  ! type of NetCDF variable
      integer   NVDims      ! number of NetCDF dimensions
      integer   VDims    ( MaxVDims )
                            ! NetCDF variable dimension ids
      integer   DimSz       ! Dimension size
      integer   DimID       ! Dimension id
      integer   NVAtts      ! number of variable attribute
      integer   NC_NStr     ! the allocated number of list
                            !   entries allocated in the file
      integer   NC_StrSz    ! the allocated size of each string
                            !   or list entry allocated
      integer   iList       ! list entry index number
      integer   iChar0      ! location of substring in temporary
                            !   storage space
      integer   NBlocks     ! number of blocks or chunks necessary
                            !   to write the data to the file
                            !   by using in scratch space in
                            !   the variable, strtemp, ( defined
                            !   in the header file, ods_hdf.h )
      integer   iBlock      ! block index number for do loop
      integer   NChar       ! total number of characters in a
                            !   block of list entries
      integer   NSectors    ! number of sectors in the scatch
                            !   space.  Each sector is NC_StrSz
                            !   in length
      integer   iSector     ! sector index number for do loop

*     Default value for error code
*     ----------------------------
      ierr         = NCNoErr

*     Check to determine if the file handle id is valid
*     -------------------------------------------------
      if ( id            .lt. 1      .or.
     .     id            .gt. id_max .or.
     .     IOMode ( id ) .eq. CLOSED ) then
         write ( stderr, 901 )
         ierr = NCEBadID
         return
      end if

*     Set NetCDF file id
*     ------------------
      nc_id        = ncid  ( id )

*     Set NetCDF variable id
*     ----------------------
      varid        = ncvid ( nc_id, ListName, ierr )
      if ( ierr .ne. NCNoErr ) return

*     Get information about the NetCDF variable
*     -----------------------------------------
      call ncvinq ( nc_id,  varid,
     .              VarNam, NC_VarType,
     .              NVDims, VDims,
     .              NVAtts, ierr )
      if ( ierr .ne. NCNoErr ) return

*     Determine the hyperslab dimension in the NetCDF file
*     ----------------------------------------------------
      DimID        = VDims ( 1 )
      call NCDINQ ( nc_id, DimID, VarNam, DimSz, ierr )
      if ( ierr .ne. NCNoErr ) return
      NC_StrSz     = DimSz

      DimID        = VDims ( 2 )
      call NCDINQ ( nc_id, DimID, VarNam, DimSz, ierr )
      if ( ierr .ne. NCNoErr ) return
      NC_NStr      = DimSz

*     Return with an error message if the number of
*     observations to be read exceeds the user-defined maximum
*     --------------------------------------------------------
      if ( NC_NStr .gt. ListSz ) then
         write ( stderr, * )
         write ( stderr, 902 ) ListSz, NC_NStr
         ierr = ODS_DimErr
         return

      end if

*     Define a block of data and allocate
*     the work space, strtemp
*     -----------------------------------
      NChar        = NC_StrSz * NC_NStr
      NSectors     = min ( NC_StrSz * NC_NStr, MChar ) / NC_StrSz
      NBlocks      = NC_NStr / NSectors
      if ( NBlocks * NSectors .ne. NC_NStr )
     .   NBlocks   = NBlocks + 1

*     Initialize do loop index and pointers
*     -------------------------------------
      iList        = 0
      corner ( 1 ) = 1
      corner ( 2 ) = 1
      edges  ( 1 ) = NC_StrSz
      edges  ( 2 ) = 1

*     for each block ...
*     ------------------
      do 20, iBlock = 1, NBlocks

*        Define the hyperslab of the block of list entries
*        -------------------------------------------------
         edges ( 2 ) = NSectors
         if ( iBlock .eq. NBlocks ) then   ! for the last block
                                           ! ------------------
            edges ( 2 ) = mod ( NC_NStr - 1, NSectors ) + 1
            NSectors    = edges ( 2 )

         end if
         NChar       = edges ( 1 ) * edges ( 2 )

*        Write the block of entries to file
*        ----------------------------------
         call ncvgtc ( nc_id,   varid, corner, edges,
     .                 strtemp, NChar, ierr )
         if ( ierr .ne. NCNoErr ) return

*        Define the corner of the hyperslab
*        for the next block of entries
*        ----------------------------------
         corner ( 2 ) = corner ( 2 ) + edges ( 2 )

*        Copy the block of data in scatch space to the output array
*        ----------------------------------------------------------
         iChar0 = CLEAR
         do 10, iSector = 1, NSectors
            iList  = iList + 1
            List ( iList )
     .          = strtemp ( iChar0 + 1 : iChar0 + NC_StrSz )

            iChar0 = iChar0 + NC_StrSz

 10      continue

 20   continue

*     return the list size
*     --------------------
      ListSz = NC_NStr

      return
*     ------

 901  format ( /, ' ODS_GetList: File handle id number does not ',
     .         /, '              correspond to an opened ODS file' )
 902  format ( /, ' ODS_GetList: The number of list entries ',
     .         /, '              to be read exceeds the user ',
     .         /, '              defined maximum of ', i10, '.',
     .         /, '              Number of list entries to be ',
     .         /, '              read = ', i10, '.' )

      end

max_stars_repo_name: GEOS-ESM/AeroApps

max_forks_repo_name: GEOS-ESM/AeroApps


!! CONTROLA LAS CONDICIONES DE FRONTERA !!

subroutine boundaries_calc

    use arrays
    use global

    implicit none

    !modos R, L y extrapolacion de 2 orden (en la frontera)
    Rmode     = 0.5d0*(Pi - Psi)
    Rmode(0)  = 3.0d0*Rmode(1) - 3.0d0*Rmode(2) + Rmode(3)
    Rmode(Nx) = 3.0d0*Rmode(Nx-1) - 3.0d0*Rmode(Nx-2) + Rmode(Nx-3)
    Lmode     = 0.5d0*(Pi + Psi)
    Lmode(0)  = 3.0d0*Lmode(1) - 3.0d0*Lmode(2) + Lmode(3)
    Lmode(Nx) = 3.0d0*Lmode(Nx-1) - 3.0d0*Lmode(Nx-2) + Lmode(Nx-3)

    !Condiciones de frontera
    if (Boundaries.eq.0) then       !cd. onda saliente
        Psi(0)  = Lmode(0)
        Psi(Nx) = -Rmode(Nx)
        Pi(0)   = Lmode(0)
        Pi(Nx)  = Rmode(Nx)
    else if (Boundaries.eq.1) then  !cd. onda reflexiva
        Psi(0)  = -Psi(1)
        Psi(Nx) = -Psi(Nx-1)
        Pi(0)   = -Pi(1)
        Pi(Nx)  = -Pi(Nx-1)
    else if (Boundaries.eq.2) then  !cd. onda periodica
        Psi(0)  = Psi(Nx-1)
        Psi(Nx) = Psi(1)
        Pi(0)   = Pi(Nx-1)
        Pi(Nx)  = Pi(1)
    else if (Boundaries.eq.3) then  !cd. agujero negro (desaparece)
        Psi(0)  = 3.0d0*Psi(1) - 3.0d0*Psi(2) + Psi(3)
        Psi(Nx) = -Rmode(Nx)
        Pi(0)   = 3.0d0*Pi(1)  - 3.0d0*Pi(2)  + Pi(3)
        Pi(Nx)  = Rmode(Nx)
    else
        print*, 'NOT IMPLEMENTED YET'
    end if

end subroutine boundaries_calc


text:
      subroutine minpwr(pc,nd,npc,nt,x,pwrmin,nomit)

      integer nd,npc,i,j,k,nt,nomit
      real pwrmin,pwr
      complex pc(nd,2),x(nd,nt)
      pwravg = 0
      do i = 1,nd
         pwr = 0
         do j = 1,npc
            pwr = pwr + pc(i,j)*conjg(pc(i,j))
         enddo
         pwravg = pwravg+pwr
         if(pwr .lt. pwrmin) then
            do k = 1,nt
               x(j,k) = x(j,k)*1.e-10
            enddo
            nomit = nomit+1
         endif
      enddo

      write(*,*) 'Average Power= ',pwravg/nd,'    # omitted = ',nomit

      return
      end

max_stars_repo_name: simpeg-research/iris-mt-scratch

max_issues_repo_name: simpeg-research/iris-mt-scratch

max_forks_repo_name: simpeg-research/iris-mt-scratch