text:
      SUBROUTINE MP_BORN_AMPS_AND_WFS(P)
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C
C     Computes all the AMP and WFS in quadruple precision for the
C     phase space point P(0:3,NEXTERNAL)
C
C     Process: d d~ > t t~ QED=0 QCD=2 [ virt = QCD ] WEIGHTED=6
C
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=1)
      INTEGER    NLOOPAMPS, NCTAMPS
      PARAMETER (NLOOPAMPS=43, NCTAMPS=29)
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=4)
      INTEGER    NWAVEFUNCS
      PARAMETER (NWAVEFUNCS=6)
      INTEGER    NCOMB
      PARAMETER (NCOMB=16)
      REAL*16     ZERO
      PARAMETER (ZERO=0E0_16)
      COMPLEX*32 IMAG1
      PARAMETER (IMAG1=(0E0_16,1E0_16))

C
C     ARGUMENTS
C
      REAL*16 P(0:3,NEXTERNAL)
C
C     LOCAL VARIABLES
C
      INTEGER I,J,H
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      DATA IC/NEXTERNAL*1/

C
C     GLOBAL VARIABLES
C
      INCLUDE 'mp_coupl_same_name.inc'

      INTEGER NTRY
      LOGICAL CHECKPHASE,HELDOUBLECHECKED
      REAL*8 REF
      COMMON/INIT/NTRY,CHECKPHASE,HELDOUBLECHECKED,REF

      LOGICAL GOODHEL(NCOMB)
      LOGICAL GOODAMP(NLOOPAMPS,NCOMB)
      COMMON/FILTERS/GOODAMP,GOODHEL

      INTEGER HELPICKED
      COMMON/HELCHOICE/HELPICKED

      COMPLEX*32 AMP(NBORNAMPS,NCOMB)
      COMMON/MP_AMPS/AMP
      COMPLEX*16 DPAMP(NBORNAMPS,NCOMB)
      COMMON/AMPS/DPAMP
      COMPLEX*32 W(20,NWAVEFUNCS,NCOMB)
      COMMON/MP_WFS/W

      COMPLEX*32 AMPL(3,NCTAMPS)
      COMMON/MP_AMPL/AMPL

      COMPLEX*16 DPW(20,NWAVEFUNCS,NCOMB)
      COMMON/WFCTS/DPW

      COMPLEX*16 DPAMPL(3,NLOOPAMPS)
      LOGICAL S(NLOOPAMPS)
      COMMON/AMPL/DPAMPL,S

      INTEGER HELC(NEXTERNAL,NCOMB)
      COMMON/HELCONFIGS/HELC

      LOGICAL MP_DONE_ONCE
      COMMON/MP_DONE_ONCE/MP_DONE_ONCE

C     ----------
C     BEGIN CODE
C     ---------

      MP_DONE_ONCE=.TRUE.

C     To be on the safe side, we always update the MP params here.
C     It can be redundant as this routine can be called a couple of
C      times for the same PS point during the stability checks.
C     But it is really not time consuming and I would rather be safe.
      CALL MP_UPDATE_AS_PARAM()

      DO H=1,NCOMB
        IF ((HELPICKED.EQ.H).OR.((HELPICKED.EQ.-1).AND.((CHECKPHASE.OR.
     $   .NOT.HELDOUBLECHECKED).OR.GOODHEL(H)))) THEN
          DO I=1,NEXTERNAL
            NHEL(I)=HELC(I,H)
          ENDDO
          CALL MP_IXXXXX(P(0,1),ZERO,NHEL(1),+1*IC(1),W(1,1,H))
          CALL MP_OXXXXX(P(0,2),ZERO,NHEL(2),-1*IC(2),W(1,2,H))
          CALL MP_OXXXXX(P(0,3),MDL_MT,NHEL(3),+1*IC(3),W(1,3,H))
          CALL MP_IXXXXX(P(0,4),MDL_MT,NHEL(4),-1*IC(4),W(1,4,H))
          CALL MP_FFV1P0_3(W(1,1,H),W(1,2,H),GC_5,ZERO,ZERO,W(1,5,H))
C         Amplitude(s) for born diagram with ID 1
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),GC_5,AMP(1,H))
          CALL MP_FFV1P0_3(W(1,4,H),W(1,3,H),GC_5,ZERO,ZERO,W(1,6,H))
C         Counter-term amplitude(s) for loop diagram number 2
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,6,H),R2_GGQ,AMPL(1,1))
C         Counter-term amplitude(s) for loop diagram number 5
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQQ_1EPS
     $     ,AMPL(2,2))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQQ_1EPS
     $     ,AMPL(2,3))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQQ_1EPS
     $     ,AMPL(2,4))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQQ_1EPS
     $     ,AMPL(2,5))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQB,AMPL(1,6))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQQ_1EPS
     $     ,AMPL(2,7))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQT,AMPL(1,8))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQQ_1EPS
     $     ,AMPL(2,9))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),UV_GQQG_1EPS
     $     ,AMPL(2,10))
          CALL MP_FFV1_0(W(1,1,H),W(1,2,H),W(1,6,H),R2_GQQ,AMPL(1,11))
C         Counter-term amplitude(s) for loop diagram number 7
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,6,H),R2_GGQ,AMPL(1,12))
C         Counter-term amplitude(s) for loop diagram number 8
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,6,H),R2_GGQ,AMPL(1,13))
C         Counter-term amplitude(s) for loop diagram number 9
          CALL MP_R2_GG_1_0(W(1,5,H),W(1,6,H),R2_GGQ,AMPL(1,14))
C         Counter-term amplitude(s) for loop diagram number 10
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,5,H),W(1,6,H),R2_GGQ,R2_GGB
     $     ,AMPL(1,15))
C         Counter-term amplitude(s) for loop diagram number 11
          CALL MP_R2_GG_1_R2_GG_3_0(W(1,5,H),W(1,6,H),R2_GGQ,R2_GGT
     $     ,AMPL(1,16))
C         Counter-term amplitude(s) for loop diagram number 12
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQQ_1EPS
     $     ,AMPL(2,17))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQQ_1EPS
     $     ,AMPL(2,18))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQQ_1EPS
     $     ,AMPL(2,19))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQQ_1EPS
     $     ,AMPL(2,20))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQB,AMPL(1,21))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQQ_1EPS
     $     ,AMPL(2,22))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQT,AMPL(1,23))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQQ_1EPS
     $     ,AMPL(2,24))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),UV_GQQG_1EPS
     $     ,AMPL(2,25))
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),R2_GQQ,AMPL(1,26))
C         Counter-term amplitude(s) for loop diagram number 14
          CALL MP_R2_GG_1_R2_GG_2_0(W(1,5,H),W(1,6,H),R2_GGG_1
     $     ,R2_GGG_2,AMPL(1,27))
C         Amplitude(s) for UVCT diagram with ID 16
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),GC_5,AMPL(1,28))
          AMPL(1,28)=AMPL(1,28)*(2.0D0*UVWFCT_T_0)
C         Amplitude(s) for UVCT diagram with ID 17
          CALL MP_FFV1_0(W(1,4,H),W(1,3,H),W(1,5,H),GC_5,AMPL(2,29))
          AMPL(2,29)=AMPL(2,29)*(2.0D0*UVWFCT_B_0_1EPS)
C         Copy the qp wfs to the dp ones as they are used to setup the
C          CT calls.
          DO I=1,NWAVEFUNCS
            DO J=1,20
              DPW(J,I,H)=W(J,I,H)
            ENDDO
          ENDDO
C         Same for the counterterms amplitudes
          DO I=1,NCTAMPS
            DO J=1,3
              DPAMPL(J,I)=AMPL(J,I)
              S(I)=.TRUE.
            ENDDO
          ENDDO
          DO I=1,NBORNAMPS
            DPAMP(I,H)=AMP(I,H)
          ENDDO
        ENDIF
      ENDDO

      END

ext: f


text:
      SUBROUTINE FRSTNZ(IAPADR,NSAMPS,NELEM,ISCR)
C     FRSTNZ RETURNS THE ELEMENT COUNT OF THE FIRST NONZERO ELEMENT OF AN ARRAY
C  STARTING AT AP LOCATION IAPADR.
C
C  ARGUMENTS:
C  IAPADR - THE AP ADDRESS OF THE BEGINNING OF THE ARRAY TO BE SCANNED.
C  NSAMPS - THE NUMBER OF SAMPLES OR ELEMENTS TO BE SCANNED.
C  NELEM  - THE ELEMENT COUNT COMPUTED BY FRSTNZ.  THE FIRST ELEMENT
C           (AT AP ADDRESS IAPADR) IS 1.
C  ISCR   - THE AP ADDRESS OF AN AP SCRATCH ARRAY AT LEAST NSAMPS+2 LONG.
C
C  WRITTEN AND COPYRIGHTED BY
C  PAUL HENKART, SCRIPPS INSTITUTION OF OCEANOGRAPHY, 2 JANUARY 1980
C  ALL RIGHTS ARE RESERVED BY THE AUTHOR.  PERMISSION TO COPY OR REPRODUCE THIS
C  SUBROUTINE, BY COMPUTER OR OTHER MEANS, MAY BE OBTAINED ONLY FROM THE AUTHOR.
C
      CALL LVNOT(IAPADR,1,ISCR,1,NSAMPS)                                ! /* ZEROES=1 AND NONZEROES=0
      ISCR1=ISCR+NSAMPS                                                 ! /* GET ADDITIONAL SCRATCH SPACE (2 WORDS)
      CALL MINV(ISCR,1,ISCR1,NSAMPS)                                    ! /* FIND THE MINIMUM VALUE (THE FIRST ZERO)
      ISCR2=ISCR1+1                                                     ! /* THIS POINTS TO THE AP ADDRESS OF THE MINIMUM
      CALL APWR                                                         ! /* WAIT FOR RESULTS TO BE COMPLETED
      CALL APGET(A,ISCR2,1,2)                                           ! /* THE ADDRESS IS FLOATING POINT!
      CALL APWD                                                         ! /* WAIT FOR DATA TRANSMISSION COMPLETION
      NELEM=A                                                           ! /* COMPUTE THE ELEMENT COUNT FROM THE BEGINNING OF THE ARRAY
      RETURN
      END

max_stars_repo_name: schwehr/sioseis

max_issues_repo_name: schwehr/sioseis

max_forks_repo_name: schwehr/sioseis


text:
module mod_latlon

 !===================================================================================85
 !---------------------------- DESCRIPTION ------------------------------------------85
 !
 !
 !
 !          Module  name:
 !                         'mod_latlon'
 !
 !          Purpose:
 !                    Latitude-longtitude geographical coordinates.
 !                    This class is simple wrapper which conveniently
 !                    handles creation of geographiv coordinates(lat-lon).
 !
 !
 !          History:
 !                        Date: 15-09-2017
 !                        Time: 09:36 GMT+2
 !
 !          Version:
 !
 !                      Major: 1
 !                      Minor: 0
 !                      Micro: 0
 !
 !          Author:  Bernard Gingold
 !
 !
 !
 !
 !
 !          E-mail:
 !
 !                      beniekg@gmail.com
 !==================================================================================85

 ! Tab:5 col - Type and etc.. definitions
 ! Tab:10,11 col - Type , function and subroutine code blocks.

    implicit none
    private ! everything is private beside variables and procedures deemed public
    use ISO_FORTRAN_ENV, only : stderr=>ERROR_UNIT , &
                                stdout=>OUTPUT_UNIT
    use IFPORT, only : TRACEBACKQQ
    use module_class_error_check, only : array2D_not_alloc

    use module_kinds
    use mod_constants
    use module_logger
    use mod_code_timing

    !=====================================================59
    !  File and module information:
    !  version,creation and build date, author,description
    !=====================================================59

    ! Version major
    integer(I32P), parameter, public :: MOD_LATLON_MAJOR = 1

    ! Version minor
    integer(I32P), parameter, public :: MOD_LATLON_MINOR = 0

    ! Version micro
    integer(I32P), parameter, public :: MOD_LATLON_MICRO = 0

    ! Module/file full version
    integer(I32P), parameter, public :: MOD_LATLON_FULLVER = 1000*MOD_LATLON_MAJOR+100*MOD_LATLON_MINOR+ &
                                                             10*MOD_LATLON_MICRO

    ! Module creation date
    character(*),  parameter, public :: MOD_LATLON_CREATE_DATE = "15-09-2017 14:09 +00200 (FRI 15 SEP 2017 GMT+2)"

    ! Module build date (should be set after successful build date/time)
    character(*),  parameter, public :: MOD_LATLON_BUILD_DATE = " "

    ! Module author info.
    character(*),  parameter, public :: MOD_LATLON_AUTHOR = "Programmer: Bernard Gingold e-mail: beniekg@gmail.com"

    ! Module short description.
    character(*),  parameter, public :: MOD_LATLON_DESCRIPT = "Latitude-longtitude geographical coordinates."

    !======================================================60
    !  Type: LatLon
    !======================================================60

    type, public :: LatLon

         ! public by default.
         ! indices
         integer(I64P) :: m_ims,m_ime,m_jms,m_jme, &
                          m_ids,m_ide,m_jds,m_jde, &
                          m_its,m_ite,m_jts,m_jte

        ! Lat-lon coordinates
        ! m_lat --  latitude, south is negative (degree)
        ! m_long -- longitude, west is negative (degree)
        real(R64P), allocatable, dimension(:,:) :: m_lat,m_long

!DIR$   ATTRIBUTES ALIGN : 32 :: m_lat
!DIR$   ATTRIBUTES ALIGN : 32 :: m_long

        ! built indicator (logical)
        logical(I32P) :: m_isbuilt

        contains

        !===========================================55
        ! Constructors and Destructor procedures.
        !===========================================55

        procedure, pass(this), public :: init

        procedure, pass(this), public :: copy

        procedure, pass(this), public :: destroy

        !============================================55
        ! Read/Write procedures
        !============================================55

        procedure, nopass, public :: read_latlon

        procedure, nopass, public :: write_latlon

    end type LatLon

       !=================================================60
       !  Module operators
       !=================================================60

       interface assignment (=)
            module procedure assign_latlon
       end interface

       interface operator (/=)
            module procedure latlon_neq_latlon
       end interface

       interface operator (==)
            module procedure latlon_eq_latlon
       end interface

       interface operator (>)
            module procedure latlon_gt_latlon
       end interface

       interface operator (<)
             module procedure latlon_lt_latlon
       end interface

    contains

    !=================================================60
    !              Implementation
    !=================================================60

    !======================================================60
    ! subroutine: init
    ! Initialziation by user(caller) passed array which
    ! usually contain specific random scalar distribution.
    !   @Warning:
    !            Upon non-fatal early exit parameter
    !            'err' will be set to -1
    !======================================================60
    subroutine init(this,indices,lat,long,logging,    &
                    profiling,filename,append,dbg,err,qpctimer)
          implicit none
          class(LatLon),                intent(inout) :: this
          integer(I64P), dimension(12), intent(in)    :: indices
          real(R64P), dimension(:,:),   intent(in)    :: lat,long
!DIR$     ASSUME_ALIGNED lat:32,long:32
          logical(I32P),                intent(in)    :: logging,profiling
          character(len=*),             intent(in)    :: filename
          logical(I32P),                intent(in)    :: append,dbg
          integer(I32P),                intent(inout) :: err
          type(QPCTimer_t),             intent(inout) :: qpctimer
          ! Locals
          character(len=40)  :: dstr,tstr
          character(len=256) :: emsg
          integer(I64P)      :: j,i
          integer(BOOL)      :: ifail
          logical(I32P)      :: bfail
!DIR$  IF (USE_LOOP_BLOCKING .EQ. 1)
          integer(I64P)      :: jj,ii
!DIR$  ENDIF
          integer(I32P)      :: aerr
          ! Start of executable statements
          if(err.LT.0) err = 0
          if(array2D_not_alloc(lat) .OR. &
             array2D_not_alloc(long)) then
              err = -1
              return
          end if
          if(this%m_isbuilt .EQ. .true.) then
              if(logging) then
                  call log_startup(filename,append)
                  call log_UsrMsg("logger:182, In->mod_latlon/init: LatLon already initialized!!")
                  call log_shutdown()
              else
                  call DATE_AND_TIME(date=dstr,time=tstr)
                  write(stderr,*) "===========================NON-FATAL=========================="
                  write(sdterr,*) " ( mod_latlon/init:182, LatLon already initialized!)"
                  write(stderr,*) " ( Non-Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
                  write(stderr,*) "===========================NON-FATAL=========================="
              end if
              err = -1
              return
          end if
          ! Begin construction
          this%m_ims=indices(1)
          this%m_ime=indices(2)
          this%m_jms=indices(3)
          this%m_jme=indices(4)
          this%m_ids=indices(5)
          this%m_ide=indices(6)
          this%m_jds=indices(7)
          this%m_jde=indices(8)
          this%m_its=indices(9)
          this%m_ite=indices(10)
          this%m_jts=indices(11)
          this%m_jte=indices(12)
          ! Array alocation
          associate(d1s=>this%m_ims, &
                    d1e=>this%m_ime, &
                    d2s=>this%m_jms, &
                    d2e=>this%m_jme )
              allocate(this%m_lat(d1s:d1e,d2s:d2e),   &
                       this%m_long(d1s:d1e,d2s:d2e),  &
                       STAT=aerr,                     &
                       ERRMSG=emsg )
          end associate
          if(aerr.NE.0) then
              if(logging) then
                  call log_startup(filename,append)
                  call log_UsrMsg("logger:221, In->mod_latlon/init: Failed to allocate: [m_coord] ")
                  call log_shutdown()
              else
                  call DATE_AND_TIME(date=dstr,time=tstr)
                  write(ERROR_UNIT,*) "===========================FATAL========================="
                  write(ERROR_UNIT,*) "   (mod_latlon/init:221, Failed to allocate: [m_coord)"
                  write(ERROR_UNIT,*) "   (System message:)", emsg
                  write(ERROR_UNIT,*) "   ( Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
                  write(ERROR_UNIT,*) "===========================FATAL=========================="
              end if
!DIR$         IF (SHOW_CALLSTACK .EQ. 1)
              call TRACEBACKQQ(STRING="FATAL-ERROR", USER_EXIT_CODE= -1)
!DIR$         ENDIF
              ERROR STOP "mod_latlon/init:221 -> [FATAL-ERROR]: Terminating execution!!"
          end if
          !Array initilaization
          if(profiling) then
              call qpctimer_start(qpctimer,ifail)
              if(ifail.EQ.0) then
                  write(stderr,*) "mod_latlon/init: qpctimer_start failed to query performance frequency counter!"
              end if
          end if
!DIR$  IF (USE_LOOP_BLOCKING .EQ. 1)
          do j = this%m_jms, this%m_jme, DEFAULT_BLOCK_SIZE
              do i = this%m_ims, this%m_kme, DEFAULT_BLOCK_SIZE
                  do jj = j, DEFAULT_BLOCK_SIZE
!DIR$                 SIMD VECTORLENGTHFOR(REAL(KIND=8))
                      do ii = i, DEFAULT_BLOCK_SIZE
!DIR$           IF (USE_SOFT_PREFETCH .EQ. 1)
                          call MM_PREFETCH(lat(ii+1,jj),1)
                          call MM_PREFETCH(long(ii+1,jj),1)
!DIR$           ENDIF
                          this%m_lat(ii,jj) = lat(ii,jj)
                          this%m_long(ii,jj) = long(ii,jj)
                      end do
                  end do
              end do
          end do
!DIR$ ELSE
          do j = this%m_jms, this%m_jme
!DIR$         SIMD VECTORLENGTHFOR(REAL(KIND=8))
              do i = this%m_ims, this%m_ime
!DIR$         IF (USE_SOFT_PREFETCH .EQ. 1)
                  call MM_PREFETCH(lat(i+1,j),1)
                  call MM_PREFETCH(long(i+1,j),1)
!DIR$         ENDIF
                  this%m_lat(i,j) = lat(i,j)
                  this%m_long(i,j) = long(i,j)
              end do
          end do
!DIR$     ENDIF
          if(profiling) then
                if(ifail.EQ.0) then
                    call qpctimer_stop(qpctimer,ifail)
                    call qpctimer_delta(qpctimer,bfail)
                    if(bfail .EQ. .false.) then
                         call qpctimer_print(qpctimer)
                    else
                         write(stderr,*) "mod_latlon/init: qpctimer_delta: failed to compute delta measurement!!"
                    end if
                else
                    write(stderr,*) "mod_latlon/init: Unable to read performance counter -- fatal!!"
                end if
          end if
          this%m_isbuilt = .true.
          if(dbg) then
              print*, " Latitude:  ", this%m_lat
              print*, " Longtitude:", this%m_long
          end if
    end subroutine

    !======================================================60
    ! subroutine: copy
    ! Remark:
    !           Upon detecting non-fatal error inout
    !           integer indicator 'err' will be set to -1.
    !======================================================60
    subroutine copy(this,other,logging,filename,append,err)
          implicit none
          class(LatLon),    intent(inout) :: this
          class(LatLon),    intent(in)    :: other
          logical(I32P),    intent(in)    :: logging
          character(len=*), intent(in)    :: filename
          logical(I32P),    intent(in)    :: append
          integer(I32P),    intent(inout) :: err
          ! Locals
          character(len=40) :: dstr,tstr
          ! Start of executable statements
          if(err.LT.0) err = 0
          if(this%m_isbuilt .EQ. .true.) then
              if(logging) then
                  call log_startup(filename,append)
                  call log_UsrMsg("logger:330, In->mod_latlon/copy: LatLon already initilaized!!")
                  call log_shutdown()
              else
                  call DATE_AND_TIME(date=dstr,time=tstr)
                  write(stderr,*) "===========================NON-FATAL=========================="
                  write(sdterr,*) " ( mod_latlon/copy:330, LatLon already initialized!)"
                  write(stderr,*) " ( Non-Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
                  write(stderr,*) "===========================NON-FATAL=========================="
              end if
              err = -1
              return
          end if
          ! Begin construction
          this%m_ims=other%m_ims
          this%m_ime=other%m_ime
          this%m_jms=other%m_jms
          this%m_jme=other%m_jme
          this%m_ids=other%m_ids
          this%m_ide=other%m_ide
          this%m_jds=other%m_jds
          this%m_jde=other%m_jde
          this%m_its=other%m_its
          this%m_ite=other%m_ite
          this%m_jts=other%m_jts
          this%m_jte=other%m_jte
          this%m_lat=other%m_lat
          this%m_long=other%m_long
          this%m_isbuilt = .true.
    end subroutine

    !======================================================60
    !  subroutine: destroy
    !  Remark:
    !           Upon detecting non-fatal error inout
    !           integer indicator 'err' will be set to -1.
    !======================================================60
    subroutine destroy(this,logging,filename,append,err)
          implicit none
          class(LatLon),    intent(inout) :: this
          logical(I32P),    intent(in)    :: logging
          character(len=*), intent(in)    :: filename
          logical(I32P),    intent(in)    :: append
          integer(I32P),    intent(inout) :: err
          ! Locals
          character(len=40)  :: dstr,tstr
          character(len=256) :: emsg
          integer(I32P)      :: derr
          ! Start of executable statemetns
          if(err.LT.0) err = 0
          if(this%m_isbuilt .EQ. .false.) then
              if(logging) then
                  call log_startup(filename,append)
                  call log_UsrMsg("logger:384, In->mod_latlon/destroy: LatLon already destroyed!!")
                  call log_shutdown()
              else
                  call DATE_AND_TIME(date=dstr,time=tstr)
                  write(stderr,*) "===========================NON-FATAL=========================="
                  write(sdterr,*) " ( mod_latlon/destroy:330, LatLon already destroyed!)"
                  write(stderr,*) " ( Non-Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
                  write(stderr,*) "===========================NON-FATAL=========================="
              end if
              err = -1
              return
          end if
          this%m_ims=0
          this%m_ime=0
          this%m_jms=0
          this%m_jme=0
          this%m_ids=0
          this%m_ide=0
          this%m_jds=0
          this%m_jde=0
          this%m_its=0
          this%m_ite=0
          this%m_jts=0
          this%m_jte=0
          ! Arrays deallocation
          if(ALLOCATED(this%m_lat)) then
              deallocate(this%m_lat, &
                         STAT=derr,  &
                         ERRMSG=emsg)
               if(derr.NE.0) then
                   if(logging) then
                       call log_startup(filename,append)
                       call log_UsrMsg("logger:418, In->mod_latlon/destroy: Failed to deallocate: [m_lat] !!")
                       call log_shutdown()
                   else
                       call DATE_AND_TIME(date=dstr,time=tstr)
                       write(ERROR_UNIT,*) "===========================FATAL========================="
                       write(ERROR_UNIT,*) "   (mod_latlon/destroy:418, Failed deallocate: [m_lat]!!)"
                       write(ERROR_UNIT,*) "   (System message:)", emsg
                       write(ERROR_UNIT,*) "   ( Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
                       write(ERROR_UNIT,*) "===========================FATAL=========================="
                end if
!DIR$         IF (SHOW_CALLSTACK .EQ. 1)
              call TRACEBACKQQ(STRING="FATAL-ERROR", USER_EXIT_CODE= -1)
!DIR$         ENDIF
              ERROR STOP "mod_latlon/destroy:418 -> [FATAL-ERROR]: Terminating execution!!"
              end if
          end if
          if(ALLOCATED(this%m_long)) then
              deallocate(this%m_long, &
                         STAT=derr,   &
                         ERRMSG=emsg)
              if(derr.NE.0) then
                  if(logging) then
                      call log_startup(filename,append)
                      call log_UsrMsg("logger:443, In->mod_latlon/destroy: Failed to deallocate: [m_long] !!")
                      call log_shutdown()
                  else
                       call DATE_AND_TIME(date=dstr,time=tstr)
                       write(ERROR_UNIT,*) "===========================FATAL========================="
                       write(ERROR_UNIT,*) "   (mod_latlon/destroy:443, Failed deallocate: [m_long]!!)"
                       write(ERROR_UNIT,*) "   (System message:)", emsg
                       write(ERROR_UNIT,*) "   ( Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
                       write(ERROR_UNIT,*) "===========================FATAL=========================="
                end if
!DIR$         IF (SHOW_CALLSTACK .EQ. 1)
              call TRACEBACKQQ(STRING="FATAL-ERROR", USER_EXIT_CODE= -1)
!DIR$         ENDIF
              ERROR STOP "mod_latlon/destroy:443 -> [FATAL-ERROR]: Terminating execution!!"
              end if
          end if
          this%m_isbuilt = .false.
    end subroutine

    !======================================================60
    !  subroutine: read_latlon
    !  Remark:
    !           Upon detecting non-fatal error inout
    !           integer indicator 'err' will be set to -1.
    !======================================================60
    subroutine read_latlon(this,unit,ioerr,err)
          implicit none
          class(LatLon), intent(in)    :: this
          integer,       intent(in)    :: unit
          integer(I32P), intent(inout) :: ioerr,err
          ! Strat of executable statements
          if(err.LT.0) err = 0
          if(this%m_isbuilt .EQ. .false.) then
              err = -1
              return
          end if
          READ(unit,iostat=ioerr) this
    end subroutine

    !======================================================60
    !  subroutine: write_latlon
    !  Remark:
    !           Upon detecting non-fatal error inout
    !           integer indicator 'err' will be set to -1.
    !======================================================60
    subroutine write_latlon(this,unit,ioerr,err)
          implicit none
          class(LatLon), intent(in)    :: this
          integer,       intent(in)    :: unit
          integer(I32P), intent(inout) :: ioerr,err
          ! Start of executable statements
          if(err.LT.0) err = 0
          if(this%m_isbuilt .EQ. .false.) then
              err = -1
              return
          end if
          WRITE(unit,iostat=ioerr) this
    end subroutine

    !======================================================60
    !                   Module operators
    !======================================================60

    !=======================================================60
    !  subroutine: overloaded assignment (=)
    !=======================================================60
    subroutine assign_latlon(this,other)
          implicit none
          type(LatLon), intent(inout) :: this
          type(LatLon), intent(in)    :: other
          ! Start of executable statements
          if(LOC(this).EQ.LOC(other)) then
               call DATE_AND_TIME(date=dstr,time=tstr)
               write(stderr,*) "===========================NON-FATAL=========================="
               write(sdterr,*) " ( mod_latlon/assignment(=):519, Attempted self_assignment!)"
               write(stderr,*) " ( Non-Fatal Error at:) ",  &
                                     sdate(1:4),"-",sdate(5:6),"-",sdate(7:8), " ", &
                                     stime(1:2),":",stime(3:4),":",stime(5:6)
               write(stderr,*) "===========================NON-FATAL=========================="
               return
          end if
          this%m_ims=other%m_ims
          this%m_ime=other%m_ime
          this%m_jms=other%m_jms
          this%m_jme=other%m_jme
          this%m_ids=other%m_ids
          this%m_ide=other%m_ide
          this%m_jds=other%m_jds
          this%m_jde=other%m_jde
          this%m_its=other%m_its
          this%m_ite=other%m_ite
          this%m_jts=other%m_jts
          this%m_jte=other%m_jte
          this%m_lat=other%m_lat
          this%m_long=other%m_long
          this%m_isbuilt = .true.
    end subroutine

    !======================================================60
    ! function: latlon_neq_latlon i.e. overloaded
    !           operator (/=)
    ! Warning:  No error checking is made!!
    !           Both arguments should be the same i.e.
    !           Same arrays size.
    !======================================================60
    function latlon_neq_latlon(c1,c2) result(neq)
         use module_helper_fields_equality
         implicit none
         type(LatLon), intent(in) :: c1,c2
         type(BoolCompField2D)    :: neq
         ! Start of executable statemetns
         neq%e1   = c1%m_lat  /= c2%m_lat
         neq%e2 = c1%m_long /= c2%m_long
    end function

    !======================================================60
    ! function: latlon_eq_latlon i.e. overloaded
    !           operator (==)
    ! Warning:  No error checking is made!!
    !           Both arguments should be the same i.e.
    !           Same arrays size.
    !======================================================60
    function latlon_eq_latlon(c1,c2) result(eq)
         use module_helper_fields_equality
         implicit none
         type(LatLon), intent(in) :: c1,c2
         type(BoolCompField2D)    :: eq
         ! Start of executable statements
         eq%e1 = c1%m_lat  == c2%m_lat
         eq%e2 = c1%m_long == c2%m_long
    end function

    !======================================================60
    ! function: latlon_gt_latlon i.e. overloaded
    !           operator (>)
    ! Warning:  No error checking is made!!
    !           Both arguments should be the same i.e.
    !           Same arrays size.
    !======================================================60
    function latlon_gt_latlon(c1,c2) result(gt)
          use module_helper_fields_equality
          implicit none
          type(LatLon), intent(in) :: c1,c2
          type(BoolCompField2D)    :: eq
          ! Start of executable stateemtns
          gt%e1 = c1%m_lat  > c2%m_lat
          gt%e2 = c1%m_long > c2%m_long
    end function

    !======================================================60
    ! function: latlon_lt_latlon i.e. overloaded
    !           operator (<)
    ! Warning:  No error checking is made!!
    !           Both arguments should be the same i.e.
    !           Same arrays size.
    !======================================================60
    function latlon_gt_latlon(c1,c2) result(lt)
          use module_helper_fields_equality
          implicit none
          type(LatLon), intent(in) :: c1,c2
          type(BoolCompField2D)    :: eq
          ! Start of executable statements
          lt%eq1 = c1%m_lat  < c2%m_lat
          lt%eq2 = c1%m_long < c2%m_long
    end function
end module mod_latlon

max_stars_repo_name: bgin/LibAtmosModel

max_issues_repo_name: bgin/LibAtmosModel

max_forks_repo_name: bgin/LibAtmosModel


text:
module type_face_info
    use mod_kinds,  only: ik

    !> Container for locating a given face in the mesh data, via indices.
    !!
    !!
    !!  @author Nathan A. Wukie
    !!  @date   2/1/2016
    !!
    !!
    !----------------------------------------------------------------------------------
    type, public :: face_info_t

        integer(ik)     :: idomain_g
        integer(ik)     :: idomain_l
        integer(ik)     :: ielement_g
        integer(ik)     :: ielement_l
        integer(ik)     :: iface        ! Element-local face index
        integer(ik)     :: dof_start

    end type face_info_t
    !**********************************************************************************

    interface face_info
        module procedure face_info_constructor
    end interface

contains

    !>
    !!
    !!  @author Nathan A. Wukie (AFRL)
    !!  @date   7/21/2017
    !!
    !----------------------------------------------------------------------------------
    function face_info_constructor(idomain_g, idomain_l, ielement_g, ielement_l, iface, dof_start) result(face_info)
        integer(ik),    intent(in)  :: idomain_g
        integer(ik),    intent(in)  :: idomain_l
        integer(ik),    intent(in)  :: ielement_g
        integer(ik),    intent(in)  :: ielement_l
        integer(ik),    intent(in)  :: iface
        integer(ik),    intent(in)  :: dof_start

        type(face_info_t)   :: face_info

        face_info%idomain_g  = idomain_g
        face_info%idomain_l  = idomain_l
        face_info%ielement_g = ielement_g
        face_info%ielement_l = ielement_l
        face_info%iface      = iface
        face_info%dof_start  = dof_start

    end function face_info_constructor
    !**********************************************************************************

end module type_face_info

max_stars_repo_name: nwukie/ChiDG

max_issues_repo_name: nwukie/ChiDG

max_forks_repo_name: nwukie/ChiDG


text:
module screening_module

  implicit none

  private
  public :: screen5, screenz, add_screening_factor, screening_init, &
            plasma_state, fill_plasma_state

  ! Constants
  double precision, parameter :: ZERO    =  0.0d0
  double precision, parameter :: ONE     =  1.0d0
  double precision, parameter :: TWO     =  2.0d0
  double precision, parameter :: THREE   =  3.0d0
  double precision, parameter :: FOUR    =  4.0d0
  double precision, parameter :: FIVE    =  5.0d0
  double precision, parameter :: SIX     =  6.0d0
  double precision, parameter :: SEVEN   =  7.0d0
  double precision, parameter :: EIGHT   =  8.0d0
  double precision, parameter :: NINE    =  9.0d0
  double precision, parameter :: TEN     = 10.0d0

  double precision, parameter :: ELEVEN  = 11.0d0
  double precision, parameter :: TWELVE  = 12.0d0
  double precision, parameter :: FIFTEEN = 15.0d0
  double precision, parameter :: SIXTEEN = 16.0d0

  double precision, parameter :: HALF    = 0.5d0
  double precision, parameter :: THIRD   = ONE/THREE
  double precision, parameter :: FOURTH  = 0.25d0
  double precision, parameter :: FIFTH   = ONE/FIVE
  double precision, parameter :: SIXTH   = ONE/SIX
  double precision, parameter :: SEVENTH = ONE/SEVEN
  double precision, parameter :: EIGHTH  = 0.125d0
  double precision, parameter :: NINETH  = ONE/NINE
  double precision, parameter :: TENTH   = 0.10d0
  double precision, parameter :: TWELFTH = ONE/TWELVE

  double precision, parameter :: TWO3RD    = TWO/THREE
  double precision, parameter :: FOUR3RD   = FOUR/THREE
  double precision, parameter :: FIVE3RD   = FIVE/THREE
  double precision, parameter :: FIVE6TH   = FIVE/SIX

  double precision, parameter :: THREE4TH  = 0.75d0

  double precision, parameter :: FIVE12TH = FIVE/TWELVE
  double precision, parameter :: SEVEN12TH = SEVEN/TWELVE

  double precision, parameter :: FIVE32ND = FIVE/32.0d0

  !! Pi
  double precision, parameter :: M_PI    = &
       3.141592653589793238462643383279502884197d0
  double precision, parameter :: M_SQRT_PI  = &
       1.772453850905516027298167483341145182798d0

  !! Roots
  double precision, parameter :: M_SQRT_2  = &
       1.414213562373095048801688724209698078570d0

  integer, parameter :: nscreen_max = 500
  integer            :: nscreen = 0

  double precision, parameter :: fact       = 1.25992104989487d0
  double precision, parameter :: co2        = THIRD * 4.248719d3
  double precision, parameter :: gamefx     = 0.3d0
  double precision, parameter :: gamefs     = 0.8d0
  double precision, parameter :: blend_frac = 0.05d0

  double precision :: z1scr(nscreen_max)
  double precision :: z2scr(nscreen_max)
  double precision :: a1scr(nscreen_max)
  double precision :: a2scr(nscreen_max)

  ! zs13    = (z1+z2)**(1./3.)
  ! zhat    = combination of z1 and z2 raised to the 5/3 power
  ! zhat2   = combination of z1 and z2 raised to the 5/12 power
  ! lzav    = log of effective charge
  ! aznut   = combination of a1,z1,a2,z2 raised to 1/3 power

  double precision :: zs13(nscreen_max)
  double precision :: zs13inv(nscreen_max)
  double precision :: zhat(nscreen_max)
  double precision :: zhat2(nscreen_max)
  double precision :: lzav(nscreen_max)
  double precision :: aznut(nscreen_max)

  type :: plasma_state

     double precision :: qlam0z
     double precision :: qlam0zdt
     double precision :: qlam0zdd

     double precision :: taufac
     double precision :: taufacdt

     double precision :: aa
     double precision :: daadt
     double precision :: daadd

  end type plasma_state

contains

  subroutine screening_init()
    ! This routine assumes that we have already filled z1scr, z2scr,
    ! a1scr, and a2scr.

    integer :: i

    do i = 1, nscreen

       zs13(i)    = (z1scr(i) + z2scr(i))**THIRD
       zs13inv(i) = ONE/zs13(i)
       zhat(i)    = (z1scr(i) + z2scr(i))**FIVE3RD  - z1scr(i)**FIVE3RD - z2scr(i)**FIVE3RD
       zhat2(i)   = (z1scr(i) + z2scr(i))**FIVE12TH - z1scr(i)**FIVE12TH -z2scr(i)**FIVE12TH
       lzav(i)    = FIVE3RD * log(z1scr(i)*z2scr(i)/(z1scr(i) + z2scr(i)))
       aznut(i)   = (z1scr(i)**2 * z2scr(i)**2 * a1scr(i)*a2scr(i) / (a1scr(i) + a2scr(i)))**THIRD

    enddo

  end subroutine screening_init

  subroutine add_screening_factor(z1, a1, z2, a2)

    ! this is only called at initialization

    double precision :: z1, a1, z2, a2

    nscreen = nscreen + 1

    z1scr(nscreen) = z1
    a1scr(nscreen) = a1
    z2scr(nscreen) = z2
    a2scr(nscreen) = a2

  end subroutine add_screening_factor

  subroutine fill_plasma_state(state, temp, dens, y)

    use network, only: nspec, zion

    ! Input variables

    type (plasma_state) :: state
    double precision :: temp, dens, y(nspec)

    ! Local variables

    double precision :: abar, zbar, z2bar
    double precision :: ytot, rr, tempi, dtempi, deni
    double precision :: pp, qq, dppdt, xni
!    double precision :: dppdd

    abar   = ONE / sum(y)
    zbar   = sum(zion * y) * abar
    z2bar  = sum(zion**2 * y) * abar

    ytot             = ONE / abar
    rr               = dens * ytot
    tempi            = ONE / temp
    dtempi           = -tempi * tempi
    deni             = ONE / dens

    pp               = sqrt(rr*tempi*(z2bar + zbar))
    qq               = HALF/pp *(z2bar + zbar)
    dppdt            = qq*rr*dtempi
    !dppdd            = qq * ytot * tempi

    state % qlam0z   = 1.88d8 * tempi * pp
    state % qlam0zdt = 1.88d8 * (dtempi*pp + tempi*dppdt)
    !state % qlam0zdd = 1.88d8 * tempi * dppdd

    state % taufac   = co2 * tempi**THIRD
    state % taufacdt = -THIRD * state % taufac * tempi

    qq               = rr * zbar
    xni              = qq**THIRD
    !dxnidd           = THIRD * xni * deni

    state % aa       = 2.27493d5 * tempi * xni
    state % daadt    = 2.27493d5 * dtempi * xni
    !state % daadd    = 2.27493d5 * tempi * dxnidd

  end subroutine fill_plasma_state

  subroutine screen5(state,jscreen,scor,scordt,scordd)

    implicit none

    ! this subroutine calculates screening factors and their derivatives
    ! for nuclear reaction rates in the weak, intermediate and strong regimes.
    ! based on graboske, dewit, grossman and cooper apj 181 457 1973 for
    ! weak screening. based on alastuey and jancovici apj 226 1034 1978,
    ! with plasma parameters from itoh et al apj 234 1079 1979, for strong
    ! screening.

    ! input:
    ! state   = plasma state (T, rho, abar, zbar, etc.)
    ! jscreen = counter of which reaction is being calculated

    ! output:
    ! scor    = screening correction
    ! scordt  = derivative of screening correction with temperature
    ! scordd  = derivative of screening correction with density

    ! declare the pass
    integer             :: jscreen
    type (plasma_state) :: state
    double precision    :: scor, scordt, scordd

    ! local variables
    double precision :: z1, a1, z2, a2

    double precision :: bb,cc,dccdt, &
                        qq,dqqdt,rr,drrdt, &
                        ss,dssdt,tt,dttdt,uu,duudt, &
                        vv,dvvdt,a3,da3, &
                        h12w,dh12wdt,h12,dh12dt, &
                        h12x,dh12xdt,alfa,beta, &
                        gamp,gampdt, &
                        gamef,gamefdt, &
                        tau12,tau12dt,alph12,alph12dt, &
                        xlgfac,dxlgfacdt, &
                        gamp14,gamp14dt
!    double precision :: dccdd,dqqdd,dvvdd,drrdd,dssdd,dttdd,duudd
!    double precision :: dh12dd,dh12wdd,dh12xdd,alph12dd
!    double precision :: gampdd,gamefdd,dxlgcfacdd,gamp14dd

    ! Get the ion data based on the input index

    z1 = z1scr(jscreen)
    a1 = a1scr(jscreen)
    z2 = z2scr(jscreen)
    a2 = a2scr(jscreen)

    ! calculate individual screening factors
    bb       = z1 * z2
    gamp     = state % aa
    gampdt   = state % daadt
    !gampdd   = state % daadd

    qq       = fact * bb * zs13inv(jscreen)
    gamef    = qq * gamp
    gamefdt  = qq * gampdt
    !gamefdd  = qq * gampdd

    tau12    = state % taufac * aznut(jscreen)
    tau12dt  = state % taufacdt * aznut(jscreen)

    qq       = ONE/tau12
    alph12   = gamef * qq
    alph12dt = (gamefdt - alph12*tau12dt) * qq
    !alph12dd = gamefdd * qq

    ! limit alph12 to 1.6 to prevent unphysical behavior.
    ! this should really be replaced by a pycnonuclear reaction rate formula
    if (alph12 .gt. 1.6) then
       alph12   = 1.6d0
       alph12dt = ZERO
       !alph12dd = ZERO

       gamef    = 1.6d0 * tau12
       gamefdt  = 1.6d0 * tau12dt
       !gamefdd  = ZERO

       qq       = zs13(jscreen)/(fact * bb)
       gamp     = gamef * qq
       gampdt   = gamefdt * qq
       !gampdd   = ZERO
    end if

    ! weak screening regime
    h12w    = bb * state % qlam0z
    dh12wdt = bb * state % qlam0zdt
    !dh12wdd = bb * qlam0zdd

    h12     = h12w
    dh12dt  = dh12wdt
    !dh12dd  = dh12wdd

    ! intermediate and strong sceening regime
    if (gamef .gt. gamefx) then

       gamp14   = gamp**FOURTH
       rr       = ONE/gamp
       qq       = 0.25d0*gamp14*rr
       gamp14dt = qq * gampdt
       !gamp14dd = qq * gampdd

       cc       =   0.896434d0 * gamp * zhat(jscreen) &
            - 3.44740d0  * gamp14 * zhat2(jscreen) &
            - 0.5551d0   * (log(gamp) + lzav(jscreen)) &
            - 2.996d0

       dccdt    =   0.896434d0 * gampdt * zhat(jscreen) &
            - 3.44740d0  * gamp14dt * zhat2(jscreen) &
            - 0.5551d0*rr*gampdt

       !dccdd    =   0.896434d0 * gampdd * zhat(jscreen) &
       !     - 3.44740d0  * gamp14dd * zhat2(jscreen) &
       !     - 0.5551d0*rr*gampdd

       a3     = alph12 * alph12 * alph12
       da3    = 3.0d0 * alph12 * alph12

       qq     = 0.014d0 + 0.0128d0*alph12
       dqqdt  = 0.0128d0*alph12dt
       !dqqdd  = 0.0128d0*alph12dd

       rr     = FIVE32ND - alph12*qq
       drrdt  = -(alph12dt*qq + alph12*dqqdt)
       !drrdd  = -(alph12dd*qq + alph12*dqqdd)

       ss     = tau12*rr
       dssdt  = tau12dt*rr + tau12*drrdt
       !dssdd  = tau12*drrdd

       tt     =  -0.0098d0 + 0.0048d0*alph12
       dttdt  = 0.0048d0*alph12dt
       !dttdd  = 0.0048d0*alph12dd

       uu     =  0.0055d0 + alph12*tt
       duudt  = alph12dt*tt + alph12*dttdt
       !duudd  = alph12dd*tt + alph12*dttdd

       vv   = gamef * alph12 * uu
       dvvdt= gamefdt*alph12*uu + gamef*alph12dt*uu + gamef*alph12*duudt
       !dvvdd= gamefdd*alph12*uu + gamef*alph12dd*uu + gamef*alph12*duudd

       h12     = cc - a3 * (ss + vv)
       rr      = da3 * (ss + vv)
       dh12dt  = dccdt - rr*alph12dt - a3*(dssdt + dvvdt)
       !dh12dd  = dccdd - rr*alph12dd - a3*(dssdd + dvvdd)

       rr     =  ONE - 0.0562d0*a3
       ss     =  -0.0562d0*da3
       drrdt  = ss*alph12dt
       !drrdd  = ss*alph12dd

       if (rr .ge. 0.77d0) then
          xlgfac    = rr
          dxlgfacdt = drrdt
          !dxlgfacdd = drrdd
       else
          xlgfac    = 0.77d0
          dxlgfacdt = ZERO
          !dxlgfacdd = ZERO
       end if

       h12    = log(xlgfac) + h12
       rr     = ONE/xlgfac
       dh12dt = rr*dxlgfacdt + dh12dt
       !dh12dd = rr*dxlgfacdd + dh12dd

       if (gamef .le. gamefs) then
          rr     =  2.0d0*(gamefs - gamef)
          drrdt  = -2.0d0*gamefdt
          !drrdd  = -2.0d0*gamefdd

          ss     = 2.0d0*(gamef - gamefx)
          dssdt  = 2.0d0*gamefdt
          !dssdd  = 2.0d0*gamefdd

          ! store current values for possible blending
          h12x    = h12
          dh12xdt = dh12dt
          !dh12xdd = dh12dd

          vv     = h12
          h12    = h12w*rr + vv*ss
          dh12dt = dh12wdt*rr + h12w*drrdt + dh12dt*ss + vv*dssdt
          !dh12dd = dh12wdd*rr + h12w*drrdd + dh12dd*ss + vv*dssdd

          ! blend the transition region - from bill paxton
          if (gamefs - gamef .lt. blend_frac*(gamefs - gamefx)) then
             alfa   = (gamefs - gamef) / (blend_frac*(gamefs - gamefx))
             alfa   = HALF * (ONE - cos(M_PI*alfa))
             beta   = ONE - alfa
             h12    = alfa * h12 + beta * h12x
             dh12dt = alfa * dh12dt + beta * dh12xdt
             !dh12dd = alfa * dh12dd + beta * dh12xdd
          end if
       end if

       ! end of intermediate and strong screening if
    end if

    ! machine limit the output
    ! further limit to avoid the pycnonuclear regime
    h12    = max(min(h12, 30.0d0), ZERO)
    scor   = exp(h12)
    if (h12 .eq. 30.0d0) then
       scordt = ZERO
       !scordd = ZERO
    else
       scordt = scor * dh12dt
       !scordd = scor * dh12dd
    end if

  end subroutine screen5

  subroutine screenz (t,d,z1,z2,a1,a2,ymass,scfac,dscfacdt)

    use network, only: aion, zion, nspec

    implicit none

    double precision :: t, d, z1, z2, a1, a2
    double precision :: ymass(nspec)
    double precision :: scfac
    double precision :: dscfacdt

    ! this subroutine calculates screening factors for nuclear reaction
    ! rates in the weak, intermediate , and strong regimes given the
    ! temperature (t--degk), the density (d--g/cc), the atomic numbers
    ! and weights of the elements in the reaction channel with the
    ! largest coulomb barrier (z1,z2,a1,a2), and the mean plasma
    ! parameters calculated in main and passed over in common aver:
    ! (mean atomic number--zbar, mean square of the atomic
    ! number--z2bar, mean atomic weight--abar, and total number of moles
    ! of nuclei per gram--ytot1).  the unscreened rate is to be
    ! multiplied by the dimensionless the treatment is based on
    ! graboske, dewit, grossman, and cooper ap j. 181,457 (1973) for
    ! weak screening and on alastuey and jancovici, ap.j. 226, 1034,
    ! 1978, with plasma parameters from itoh, totsuji, setsuo, and
    ! dewitt, ap.j. 234, 1079,1979, for strong screening (rkw
    ! modification).

    !.... last revision 15 nov 1982

    double precision abar, zbar, ytot1, z2bar, theta

    integer iy

    double precision qlam0, ztilda, qlam0z, gamp, taufac
    double precision dqlam0dt, dqlam0zdt, dgampdt, dtaufacdt
    double precision zhat, zhat2, gamef, tau12, alph12
    double precision dgamefdt, dtau12dt, dalph12dt
    double precision h12w, h12, c, h12fac
    double precision dh12wdt, dh12dt, dcdt

    ! calculate averages for screening routine
    ! nb  y = x/a with x the mass fraction
    ! zi and ai are the nuclear chage and atomic mass number
    ! respectively

    ! this part came in through a common block in Kepler -- do it
    ! directly here
    abar=0.d0
    zbar=0.d0
    ytot1=0.d0
    z2bar=0.d0

    do iy = 1, nspec
       ytot1 = ytot1 + ymass(iy)
       z2bar = z2bar + zion(iy)**2*ymass(iy)
       abar = abar + aion(iy)*ymass(iy)
       zbar = zbar + zion(iy)*ymass(iy)
    enddo

    z2bar=z2bar/ytot1
    abar=abar/ytot1
    zbar=zbar/ytot1

    ! resume original Kepler screen...
    theta=1.d0
    ytot1=1.d0/abar

    !.... calculate average plasma parameters
    !....
    if ((z1*z2) > 0.d0) then

       qlam0=1.88d+8*sqrt(d/(abar*t**3))
       dqlam0dt=-1.5d+0*qlam0 / t

       ztilda=sqrt(z2bar+zbar*theta)

       qlam0z=qlam0*ztilda
       dqlam0zdt=ztilda*dqlam0dt

       gamp=2.27493d+5*(d*zbar*ytot1)**THIRD/t
       dgampdt=-gamp/t

       taufac=4.248719d+3/t**THIRD
       dtaufacdt=-THIRD*taufac/t

       !.... calculate screening factor
       !.... approx. for strong screening only good for alpha .lt. 1.6

       zhat=(z1+z2)**FIVE3RD-z1**FIVE3RD-z2**FIVE3RD
       zhat2=(z1+z2)**FIVE12TH-z1**FIVE12TH-z2**FIVE12TH

       gamef=2.d0**THIRD*gamp*z1*z2/(z1+z2)**THIRD
       dgamefdt=gamef*dgampdt/gamp

       tau12=taufac*(z1**2*z2**2*a1*a2/(a1+a2))**THIRD
       dtau12dt=tau12*dtaufacdt/taufac

       alph12=3.d0*gamef/tau12
       dalph12dt=alph12*(dgamefdt/gamef - dtau12dt/tau12)

       !....
       !.... limit alph12 to 1.6 to prevent unphysical behavior
       !.... (h dec. as rho inc.) at high rho.  this should really
       !.... be replaced by a pycnonuclear reaction rate formula.
       !....
       if (alph12 > 1.6d0) then

          alph12=1.6d0
          dalph12dt=0.0d0

          gamef=1.6d0*tau12/3.d0
          dgamefdt=gamef*dtau12dt/tau12

          gamp=gamef*(z1+z2)**THIRD/(2.d0**THIRD*z1*z2)
          dgampdt=gamp*dgamefdt/gamef

       endif

       h12w=z1*z2*qlam0z
       dh12wdt=h12w*dqlam0zdt/qlam0z

       h12=h12w
       dh12dt=dh12wdt

       if (gamef > 0.3d0) then

          c=0.896434d0*gamp*zhat-3.44740d0*gamp**FOURTH*zhat2- &
               0.5551d0*(log(gamp)+FIVE3RD*log(z1*z2/(z1+z2)))-2.996d0

          dcdt=0.896434d0*dgampdt*zhat- &
               3.44740d0*FOURTH*gamp**(FOURTH-1.0d0)*zhat2*dgampdt- &
               0.5551d0*dgampdt/gamp

          h12=c-(tau12/3.d0)*(5.d0*alph12**3/32.d0-0.014d0*alph12**4 &
               -0.0128d0*alph12**5)-gamef*(0.0055d0*alph12**4 &
               -0.0098d0*alph12**5+0.0048d0*alph12**6)

          dh12dt=dcdt - ((dtau12dt*alph12**3 + 3.0d0*tau12*alph12**2* &
               dalph12dt)*(5.d0/32.d0 - 0.014d0*alph12 - &
               0.0128d0*alph12**2) + tau12*alph12**3*dalph12dt*(-0.014d0 &
               - 2.d0*0.0128d0*alph12))/3.d0 -(dgamefdt*alph12**4 + 4.d0 &
               *gamef*alph12**3*dalph12dt)*(0.0055d0 - 0.0098d0*alph12 - &
               0.0048d0*alph12**2) - gamef*alph12**4*dalph12dt*(-0.0098d0 &
               + 2.d0*0.0048d0*alph12)

          h12fac=0.77d0

          h12=log(max(1.d+0-0.0562d+0*alph12**3,h12fac))+h12
          if (1.d+0-0.0562d+0*alph12**3 .gt. h12fac) then
             dh12dt=(-3.d0*0.0562d0*alph12**2*dalph12dt)/ &
                  (1.d0-0.0562d0*alph12**3) + dh12dt
          endif

          if(gamef <= 0.8d0) then

             h12=h12w*((0.8d0-gamef)/0.5d0)+h12*((gamef-0.3d0)/0.5d0)
             dh12dt=((dh12wdt*(0.8d0-gamef) - h12w*dgamefdt + dh12dt* &
                  (gamef-0.3d0) + h12*dgamefdt)/0.5d0)
          endif
       endif

       if (h12.gt.300.d0) then
          h12=300.d0
          dh12dt=0.d0
       endif

       if (h12.lt.0.d0) then
          h12=0.d0
          dh12dt=0.d0
       endif

       scfac=exp(h12)
       dscfacdt=scfac*dh12dt

    else
       scfac=1.d0
       dscfacdt=0.d0
    endif

    return

  end subroutine screenz

end module screening_module

max_stars_repo_name: zingale/pynucastro

max_issues_repo_name: zingale/pynucastro

max_forks_repo_name: zingale/pynucastro


! { dg-do compile }
!
! Tests correction to implementation of pointer function assignments.
!
! Contributed by Mikael Morin  <mikael.morin@sfr.fr>
!
program p
  integer, target :: a(3) = 2
  integer :: b(3, 3) = 1
  integer :: c

  c = 3
  func (b(2, 2)) = b ! { dg-error "Different ranks" }
  func (c) = b       ! { dg-error "Different ranks" }

contains
  function func(arg) result(r)
    integer, pointer :: r(:)
    integer :: arg

    if (arg == 1) then
      r => a
    else
      r => null()
    end if
  end function func
end program p


text:
      program test_hf_eop
      use hfeop_xyu      !module that contains data and calcuation routines.
! Program to read in HF-EOP files in IERS format, compute HF-EOP for a two week period and write out the results.
! The models are presumed to be in the directory ../models
!    John Gipson
!    2016Oct20
!    301-614-6876
!   john.m.gipson@nasa.gov
!

     implicit none
!
! functions
      integer*4 julday           !julian date
      character*21 jd_to_date_1  !convert fjday to ascii string. Assumes FJDAY is UTC.

      character*128 lhfeop_file     !input file that contains model information
      character*128 leop_out       !output file
      real*8 delta_T
      real*8 MJDAY_UTC         !UTC time tag
      real*8 FJDAY_UTC         !Julian version of above
      real*8 MJDAY_TT          !TT tag
      real*8 EOP(4,2)          !XP,YP,UT1, LOD  and their rates. Units are (uas,uas,uts, uas/Day)

      integer i                !counter
      integer j                !counter over models
      integer ind

      Delta_T=  68.931118         !Difference on 2017Nov28

! Read in several models and compute HF-EOP for each one.
! k1_xyu and s2_xyu just have a single tide and are for testing purposes.

      do j=1,6
      if(j .eq. 1) then
        lhfeop_file= "../models/k1_xyu.txt"
        leop_out="k1_xyu.out"
      else if(j .eq. 2) then
        lhfeop_file= "../models/s2_xyu.txt"
        leop_out="s2_xyu.out"
      else if(j .eq. 3) then
        lhfeop_file="../models/abn_comb_2012_xyu.txt"
        leop_out="abn_comb_2012_xyu.out"
      else if(j .eq. 4) then
        lhfeop_file="../models/2017a_astro_xyu.txt"
        leop_out="2017a_astro_xyu.out"
      else if(j .eq. 5) then
        lhfeop_file="../models/fes2012_xyu.txt"
        leop_out="fes2012_xyu.out"
      else if(j .eq. 6) then
        lhfeop_file= "../models/madzak_xyu.txt"
        leop_out="madzak_xyu.out"
      endif

      write(*,*) "Doing "//trim(lhfeop_file)

      call import_tides_xyu(lhfeop_file)
!
      open(2,file=leop_out)
! write out the results for the CONT17 time period. 2017NOV28 through 2017DEC12
      write(2,*) "HFEOP predictions from model "//trim(lhfeop_file)
      write(2,*) "X, Y, UT1, LOD units= (uas,uas,us, us/day) "
      write(2,*) "Time derivatives multiplied by 1000"

      write(2,'("    MJD          Time            |     X      Y       UT1   LOD     |    Xdot    Ydot   UT1dot  LODdot |")')
      MJDAY_UTC=58085           ! 2017NOV28
! write out the EOP for two weeks
      do  i=1, 15*24
         MJDAY_TT=MJDAY_UTC+Delta_T/86400.d0
         call calc_hf_eop_xyu(MJDAY_TT, Delta_T, eop)
         fjday_UTC=mjday_UTC+2400000.5                !convert to modified Julian day  for use by jd_to_date_1
         write(2,'(f10.3,1x,a," | ",2(4f8.2," | ") )') mjday_UTC, jd_to_date_1(fjday_UTC) , &
      &       eop(1:4,1), eop(1:4,2)*1000.d0

         mjday_UTC=mjday_UTC+1.d0/24.d0
      end do
      close(2)

      end do      !end looping over models.
      end program

max_stars_repo_name: gahjelle/where

max_issues_repo_name: gahjelle/where

max_forks_repo_name: gahjelle/where


text:
        program main
        include 'mpif.h'
        integer count, errcnt, size, rank, ierr, i
        integer comm
        logical fnderr
        integer max_size
        integer world_rank
        parameter (max_size=100)
        integer intin(max_size), intout(max_size), intsol(max_size)
        real    realin(max_size), realout(max_size), realsol(max_size)
        double precision dblein(max_size), dbleout(max_size),
     *                   dblesol(max_size)
        complex cplxin(max_size), cplxout(max_size), cplxsol(max_size)
        logical login(max_size), logout(max_size), logsol(max_size)
C
C
C
C       Declare work areas
C
        call MPI_INIT( ierr )

        errcnt = 0
        comm = MPI_COMM_WORLD
        call MPI_COMM_RANK( comm, rank, ierr )
        world_rank = rank
        call MPI_COMM_SIZE( comm, size, ierr )
        count = 10

C Test sum
        if (world_rank .eq. 0) print *, ' MPI_SUM'

       fnderr = .false.
       do 23000 i=1,count
        intin(i) = i
        intsol(i) = i*size
        intout(i) = 0
23000   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_SUM, comm, ierr )
              do 23001 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23001   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_SUM'
        endif

       fnderr = .false.
       do 23002 i=1,count
        realin(i) = i
        realsol(i) = i*size
        realout(i) = 0
23002   continue
       call MPI_Allreduce( realin, realout, count,
     *      MPI_REAL, MPI_SUM, comm, ierr )
              do 23003 i=1,count
        if (realout(i).ne.realsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23003   continue
        if (fnderr) then
          print *, 'Error for type MPI_REAL and op MPI_SUM'
        endif

       fnderr = .false.
       do 23004 i=1,count
        dblein(i) = i
        dblesol(i) = i*size
        dbleout(i) = 0
23004   continue
       call MPI_Allreduce( dblein, dbleout, count,
     *      MPI_DOUBLE_PRECISION, MPI_SUM, comm, ierr )
              do 23005 i=1,count
        if (dbleout(i).ne.dblesol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23005   continue
        if (fnderr) then
          print *, 'Error for type MPI_DOUBLE_PRECISION and op MPI_SUM'
        endif

       fnderr = .false.
       do 23006 i=1,count
        cplxin(i) = i
        cplxsol(i) = i*size
        cplxout(i) = 0
23006   continue
       call MPI_Allreduce( cplxin, cplxout, count,
     *      MPI_COMPLEX, MPI_SUM, comm, ierr )
              do 23007 i=1,count
        if (cplxout(i).ne.cplxsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23007   continue
        if (fnderr) then
          print *, 'Error for type MPI_COMPLEX and op MPI_SUM'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank, ' for MPI_SUM'
        endif
        errcnt = 0

C Test product
        if (world_rank .eq. 0) print *, ' MPI_PROD'

       fnderr = .false.
       do 23008 i=1,count
        intin(i) = i
        intsol(i) = (i)**(size)
        intout(i) = 0
23008   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_PROD, comm, ierr )
              do 23009 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23009   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_PROD'
        endif

       fnderr = .false.
       do 23010 i=1,count
        realin(i) = i
        realsol(i) = (i)**(size)
        realout(i) = 0
23010   continue
       call MPI_Allreduce( realin, realout, count,
     *      MPI_REAL, MPI_PROD, comm, ierr )
              do 23011 i=1,count
        if (realout(i).ne.realsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23011   continue
        if (fnderr) then
          print *, 'Error for type MPI_REAL and op MPI_PROD'
        endif

       fnderr = .false.
       do 23012 i=1,count
        dblein(i) = i
        dblesol(i) = (i)**(size)
        dbleout(i) = 0
23012   continue
       call MPI_Allreduce( dblein, dbleout, count,
     *      MPI_DOUBLE_PRECISION, MPI_PROD, comm, ierr )
              do 23013 i=1,count
        if (dbleout(i).ne.dblesol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23013   continue
        if (fnderr) then
          print *, 'Error for type MPI_DOUBLE_PRECISION and op MPI_PROD'
        endif

       fnderr = .false.
       do 23014 i=1,count
        cplxin(i) = i
        cplxsol(i) = (i)**(size)
        cplxout(i) = 0
23014   continue
       call MPI_Allreduce( cplxin, cplxout, count,
     *      MPI_COMPLEX, MPI_PROD, comm, ierr )
              do 23015 i=1,count
        if (cplxout(i).ne.cplxsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23015   continue
        if (fnderr) then
          print *, 'Error for type MPI_COMPLEX and op MPI_PROD'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank, ' for MPI_PROD'
        endif
        errcnt = 0

C  Test max
        if (world_rank .eq. 0) print *, ' MPI_MAX'

       fnderr = .false.
       do 23016 i=1,count
        intin(i) = (rank + i)
        intsol(i) = (size - 1 + i)
        intout(i) = 0
23016   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_MAX, comm, ierr )
              do 23017 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23017   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_MAX'
        endif

       fnderr = .false.
       do 23018 i=1,count
        realin(i) = (rank + i)
        realsol(i) = (size - 1 + i)
        realout(i) = 0
23018   continue
       call MPI_Allreduce( realin, realout, count,
     *      MPI_REAL, MPI_MAX, comm, ierr )
              do 23019 i=1,count
        if (realout(i).ne.realsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23019   continue
        if (fnderr) then
          print *, 'Error for type MPI_REAL and op MPI_MAX'
        endif

       fnderr = .false.
       do 23020 i=1,count
        dblein(i) = (rank + i)
        dblesol(i) = (size - 1 + i)
        dbleout(i) = 0
23020   continue
       call MPI_Allreduce( dblein, dbleout, count,
     *      MPI_DOUBLE_PRECISION, MPI_MAX, comm, ierr )
              do 23021 i=1,count
        if (dbleout(i).ne.dblesol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23021   continue
        if (fnderr) then
          print *, 'Error for type MPI_DOUBLE_PRECISION and op MPI_MAX'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank, ' for MPI_MAX'
        endif
        errcnt = 0

C Test min
        if (world_rank .eq. 0) print *, ' MPI_MIN'

       fnderr = .false.
       do 23022 i=1,count
        intin(i) = (rank + i)
        intsol(i) = i
        intout(i) = 0
23022   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_MIN, comm, ierr )
              do 23023 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23023   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_MIN'
        endif

       fnderr = .false.
       do 23024 i=1,count
        realin(i) = (rank + i)
        realsol(i) = i
        realout(i) = 0
23024   continue
       call MPI_Allreduce( realin, realout, count,
     *      MPI_REAL, MPI_MIN, comm, ierr )
              do 23025 i=1,count
        if (realout(i).ne.realsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23025   continue
        if (fnderr) then
          print *, 'Error for type MPI_REAL and op MPI_MIN'
        endif

       fnderr = .false.
       do 23026 i=1,count
        dblein(i) = (rank + i)
        dblesol(i) = i
        dbleout(i) = 0
23026   continue
       call MPI_Allreduce( dblein, dbleout, count,
     *      MPI_DOUBLE_PRECISION, MPI_MIN, comm, ierr )
              do 23027 i=1,count
        if (dbleout(i).ne.dblesol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23027   continue
        if (fnderr) then
          print *, 'Error for type MPI_DOUBLE_PRECISION and op MPI_MIN'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank, ' for MPI_MIN'
        endif
        errcnt = 0

C Test LOR
        if (world_rank .eq. 0) print *, ' MPI_LOR'

       fnderr = .false.
       do 23028 i=1,count
        login(i) = (mod(rank,2) .eq. 1)
        logsol(i) = (size .gt. 1)
        logout(i) = .FALSE.
23028   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LOR, comm, ierr )
              do 23029 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23029   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LOR'
        endif

        if (errcnt .gt. 0) then
           print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_LOR(0)'
        endif
        errcnt = 0

       fnderr = .false.
       do 23030 i=1,count
        login(i) = .false.
        logsol(i) = .false.
        logout(i) = .FALSE.
23030   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LOR, comm, ierr )
              do 23031 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23031   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LOR'
        endif

        if (errcnt .gt. 0) then
           print *, 'Found ', errcnt, ' errors on ', rank,
     *              ' for MPI_LOR(1)'
        endif
        errcnt = 0

C Test LXOR
        if (world_rank .eq. 0) print *, ' MPI_LXOR'

       fnderr = .false.
       do 23032 i=1,count
        login(i) = (rank .eq. 1)
        logsol(i) = (size .gt. 1)
        logout(i) = .FALSE.
23032   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LXOR, comm, ierr )
              do 23033 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23033   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LXOR'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ',errcnt,' errors on ', rank, ' for MPI_LXOR'
        endif
        errcnt = 0

       fnderr = .false.
       do 23034 i=1,count
        login(i) = .false.
        logsol(i) = .false.
        logout(i) = .FALSE.
23034   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LXOR, comm, ierr )
              do 23035 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23035   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LXOR'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ',errcnt,' errors on ',rank,' for MPI_LXOR(0)'
        endif
        errcnt = 0

       fnderr = .false.
       do 23036 i=1,count
        login(i) = .true.
        logsol(i) = mod(size,2) .ne. 0
        logout(i) = .FALSE.
23036   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LXOR, comm, ierr )
              do 23037 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23037   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LXOR'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ',errcnt,' errors on ',rank,' for MPI_LXOR(1-0)'
        endif
        errcnt = 0

C Test LAND
        if (world_rank .eq. 0) print *, ' MPI_LAND'

       fnderr = .false.
       do 23038 i=1,count
        login(i) = (mod(rank,2) .eq. 1)
        logsol(i) = .false.
        logout(i) = .FALSE.
23038   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LAND, comm, ierr )
              do 23039 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23039   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LAND'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank, ' for MPI_LAND'
        endif
        errcnt = 0

       fnderr = .false.
       do 23040 i=1,count
        login(i) = .true.
        logsol(i) = .true.
        logout(i) = .FALSE.
23040   continue
       call MPI_Allreduce( login, logout, count,
     *      MPI_LOGICAL, MPI_LAND, comm, ierr )
              do 23041 i=1,count
        if (logout(i).neqv.logsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23041   continue
        if (fnderr) then
      print *, 'Error for type MPI_LOGICAL and op MPI_LAND'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ',errcnt,' errors on ',rank,
     *      ' for MPI_LAND(true)'
        endif
        errcnt = 0

C Test BOR
        if (world_rank .eq. 0) print *, ' MPI_BOR'
        if (size .lt. 3) then

       fnderr = .false.
       do 23042 i=1,count
        intin(i) = mod(rank,4)
        intsol(i) = size - 1
        intout(i) = 0
23042   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BOR, comm, ierr )
              do 23043 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23043   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BOR'
        endif

        else

       fnderr = .false.
       do 23044 i=1,count
        intin(i) = mod(rank,4)
        intsol(i) = 3
        intout(i) = 0
23044   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BOR, comm, ierr )
              do 23045 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23045   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BOR'
        endif

        endif
        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *           ' for MPI_BOR(1)'
        endif
        errcnt = 0

C Test BAND
        if (world_rank .eq. 0) print *, ' MPI_BAND'
C See bottom for function definitions

       fnderr = .false.
       do 23046 i=1,count
        intin(i) = ibxandval(rank,size,i)
        intsol(i) = i
        intout(i) = 0
23046   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BAND, comm, ierr )
              do 23047 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23047   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BAND'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_BAND(1)'
        endif
        errcnt = 0

       fnderr = .false.
       do 23048 i=1,count
        intin(i) = ibxandval1(rank,size,i)
        intsol(i) = 0
        intout(i) = 0
23048   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BAND, comm, ierr )
              do 23049 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23049   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BAND'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_BAND(0)'
        endif
        errcnt = 0

C Test BXOR
        if (world_rank .eq. 0) print *, ' MPI_BXOR'
C See below for function definitions

       fnderr = .false.
       do 23050 i=1,count
        intin(i) = ibxorval1(rank)
        intsol(i) = ibxorsol1(size)
        intout(i) = 0
23050   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BXOR, comm, ierr )
              do 23051 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23051   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BXOR'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_BXOR(1)'
        endif
        errcnt = 0

       fnderr = .false.
       do 23052 i=1,count
        intin(i) = 0
        intsol(i) = 0
        intout(i) = 0
23052   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BXOR, comm, ierr )
              do 23053 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23053   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BXOR'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_BXOR(0)'
        endif
        errcnt = 0

C Assumes -1 == all bits set

       fnderr = .false.
       do 23054 i=1,count
        intin(i) = (-1)
        if (mod(size,2) .eq. 0) then
            intsol(i) = 0
        else
            intsol(i) = -1
        endif
        intout(i) = 0
23054   continue
       call MPI_Allreduce( intin, intout, count,
     *      MPI_INTEGER, MPI_BXOR, comm, ierr )
              do 23055 i=1,count
        if (intout(i).ne.intsol(i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23055   continue
        if (fnderr) then
          print *, 'Error for type MPI_INTEGER and op MPI_BXOR'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_BXOR(1-0)'
        endif
        errcnt = 0

C Test Maxloc
        if (world_rank .eq. 0) print *, ' MPI_MAXLOC'

        fnderr = .false.
        do 23056 i=1, count
           intin(2*i-1) = (rank + i)
           intin(2*i)   = rank
           intsol(2*i-1) = (size - 1 + i)
           intsol(2*i) = (size-1)
           intout(2*i-1) = 0
           intout(2*i)   = 0
23056   continue
                call MPI_Allreduce( intin, intout, count,
     *      MPI_2INTEGER, MPI_MAXLOC, comm, ierr )
        do 23057 i=1, count
        if (intout(2*i-1) .ne. intsol(2*i-1) .or.
     *      intout(2*i) .ne. intsol(2*i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23057   continue
        if (fnderr) then
        print *, 'Error for type MPI_2INTEGER and op MPI_MAXLOC'
        endif

        fnderr = .false.
        do 23058 i=1, count
           realin(2*i-1) = (rank + i)
           realin(2*i)   = rank
           realsol(2*i-1) = (size - 1 + i)
           realsol(2*i) = (size-1)
           realout(2*i-1) = 0
           realout(2*i)   = 0
23058   continue
                call MPI_Allreduce( realin, realout, count,
     *      MPI_2REAL, MPI_MAXLOC, comm, ierr )
        do 23059 i=1, count
        if (realout(2*i-1) .ne. realsol(2*i-1) .or.
     *      realout(2*i) .ne. realsol(2*i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23059   continue
        if (fnderr) then
        print *, 'Error for type MPI_2REAL and op MPI_MAXLOC'
        endif

        fnderr = .false.
        do 23060 i=1, count
           dblein(2*i-1) = (rank + i)
           dblein(2*i)   = rank
           dblesol(2*i-1) = (size - 1 + i)
           dblesol(2*i) = (size-1)
           dbleout(2*i-1) = 0
           dbleout(2*i)   = 0
23060   continue
                call MPI_Allreduce( dblein, dbleout, count,
     *      MPI_2DOUBLE_PRECISION, MPI_MAXLOC, comm, ierr )
        do 23061 i=1, count
        if (dbleout(2*i-1) .ne. dblesol(2*i-1) .or.
     *      dbleout(2*i) .ne. dblesol(2*i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23061   continue
        if (fnderr) then
           print *,
     *     'Error for type MPI_2DOUBLE_PRECISION and op MPI_MAXLOC'

        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_MAXLOC'
        endif
        errcnt = 0

C Test minloc
        if (world_rank .eq. 0) print *, ' MPI_MINLOC'

        fnderr = .false.
        do 23062 i=1, count
           intin(2*i-1) = (rank + i)
           intin(2*i)   = rank
           intsol(2*i-1) = i
           intsol(2*i) = 0
           intout(2*i-1) = 0
           intout(2*i)   = 0
23062   continue
                call MPI_Allreduce( intin, intout, count,
     *      MPI_2INTEGER, MPI_MINLOC, comm, ierr )
        do 23063 i=1, count
        if (intout(2*i-1) .ne. intsol(2*i-1) .or.
     *      intout(2*i) .ne. intsol(2*i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23063   continue
        if (fnderr) then
        print *, 'Error for type MPI_2INTEGER and op MPI_MINLOC'
        endif

        fnderr = .false.
        do 23064 i=1, count
           realin(2*i-1) = (rank + i)
           realin(2*i)   = rank
           realsol(2*i-1) = i
           realsol(2*i) = 0
           realout(2*i-1) = 0
           realout(2*i)   = 0
23064   continue
                call MPI_Allreduce( realin, realout, count,
     *      MPI_2REAL, MPI_MINLOC, comm, ierr )
        do 23065 i=1, count
        if (realout(2*i-1) .ne. realsol(2*i-1) .or.
     *      realout(2*i) .ne. realsol(2*i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23065   continue
        if (fnderr) then
        print *, 'Error for type MPI_2REAL and op MPI_MINLOC'
        endif

        fnderr = .false.
        do 23066 i=1, count
           dblein(2*i-1) = (rank + i)
           dblein(2*i)   = rank
           dblesol(2*i-1) = i
           dblesol(2*i) = 0
           dbleout(2*i-1) = 0
           dbleout(2*i)   = 0
23066   continue
                call MPI_Allreduce( dblein, dbleout, count,
     *      MPI_2DOUBLE_PRECISION, MPI_MINLOC, comm, ierr )
        do 23067 i=1, count
        if (dbleout(2*i-1) .ne. dblesol(2*i-1) .or.
     *      dbleout(2*i) .ne. dblesol(2*i)) then
            errcnt = errcnt + 1
            fnderr = .true.
        endif
23067   continue
        if (fnderr) then
           print *,
     *      'Error for type MPI_2DOUBLE_PRECISION and op MPI_MINLOC'
        endif

        if (errcnt .gt. 0) then
        print *, 'Found ', errcnt, ' errors on ', rank,
     *          ' for MPI_MINLOC'
        endif
        errcnt = 0

        call MPI_Finalize( ierr )
        end

        integer function ibxorval1( ir )
        ibxorval1 = 0
        if (ir .eq. 1) ibxorval1 = 16+32+64+128
        return
        end

        integer function ibxorsol1( is )
        ibxorsol1 = 0
        if (is .gt. 1) ibxorsol1 = 16+32+64+128
        return
        end

C
C       Assumes -1 == all bits set
        integer function ibxandval( ir, is, i )
        integer ir, is, i
        ibxandval = -1
        if (ir .eq. is - 1) ibxandval = i
        return
        end
C
        integer function ibxandval1( ir, is, i )
        integer ir, is, i
        ibxandval1 = 0
        if (ir .eq. is - 1) ibxandval1 = i
        return
        end

max_stars_repo_name: scottkwarren/config-db

max_issues_repo_name: scottkwarren/config-db

max_forks_repo_name: scottkwarren/config-db


text:
!-----------------------------------------------------------------------
      subroutine w3utcdat(idat)
!$$$   SUBPROGRAM  DOCUMENTATION  BLOCK
!
! SUBPROGRAM: W3UTCDAT       RETURN THE UTC DATE AND TIME
!   AUTHOR: MARK IREDELL     ORG: WP23       DATE: 98-01-05
!
! ABSTRACT: THIS SUBPROGRAM RETURNS THE UTC (GREENWICH) DATE AND TIME
!   IN THE NCEP ABSOLUTE DATE AND TIME DATA STRUCTURE.
!
! PROGRAM HISTORY LOG:
!   98-01-05  MARK IREDELL
! 1999-04-28  Gilbert         - added a patch to check for the proper
!                               UTC offset.  Needed until the IBM bug
!                               in date_and_time is fixed.  The patch
!                               can then be removed.  See comments in
!                               the section blocked with "&&&&&&&&&&&".
! 1999-08-12  Gilbert         - Changed so that czone variable is saved
!                               and the system call is only done for
!                               first invocation of this routine.
!
! USAGE:  CALL W3UTCDAT(IDAT)
!
!   OUTPUT VARIABLES:
!     IDAT       INTEGER (8) NCEP ABSOLUTE DATE AND TIME
!                (YEAR, MONTH, DAY, TIME ZONE,
!                 HOUR, MINUTE, SECOND, MILLISECOND)
!
! SUBPROGRAMS CALLED:
!     DATE_AND_TIME  FORTRAN 90 SYSTEM DATE INTRINSIC
!     IW3JDN         COMPUTE JULIAN DAY NUMBER
!     W3FS26         YEAR, MONTH, DAY FROM JULIAN DAY NUMBER
!
! ATTRIBUTES:
!   LANGUAGE: FORTRAN 90
!
!$$$
      integer idat(8)
      character cdate*8,ctime*10,czone*5
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  get local date and time but use the character time zone
      call date_and_time(cdate,ctime,czone,idat)
      read(czone,'(i5)') idat(4)
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
!  convert to hours and minutes to UTC time
!  and possibly adjust the date as well
      idat(6)=idat(6)-mod(idat(4),100)
      idat(5)=idat(5)-idat(4)/100
      idat(4)=0
      if(idat(6).lt.00) then
        idat(6)=idat(6)+60
        idat(5)=idat(5)-1
      elseif(idat(6).ge.60) then
        idat(6)=idat(6)-60
        idat(5)=idat(5)+1
      endif
      if(idat(5).lt.00) then
        idat(5)=idat(5)+24
        jldayn=iw3jdn(idat(1),idat(2),idat(3))-1
        call w3fs26(jldayn,idat(1),idat(2),idat(3),idaywk,idayyr)
      elseif(idat(5).ge.24) then
        idat(5)=idat(5)-24
        jldayn=iw3jdn(idat(1),idat(2),idat(3))+1
        call w3fs26(jldayn,idat(1),idat(2),idat(3),idaywk,idayyr)
      endif
! - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -
      end

max_stars_repo_name: GMAO-SI-Team/jedi-stack

max_issues_repo_name: GMAO-SI-Team/jedi-stack


text:
Name: Kaethe Sullivan

Phone Number: (530) 7586000

Office: Downtown, Lyon Real Estate, The Tiki Hut Davis Tiki Hut

Personality: Fun, hard working and local yocal...Realtor, Oceangoer, Inventor

Website: http://www.davistikihut.com  http://www.sullivanspecialistgroup.com

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search