text:
program mpiio
  use mpi
  use, intrinsic :: iso_fortran_env, only : error_unit, output_unit
  implicit none

  integer, parameter :: datasize = 64
  integer :: rc, my_id, ntasks, localsize, i
  integer, dimension(:), allocatable :: localvector
  integer :: fh, dsize
  integer(kind=MPI_OFFSET_KIND) :: offset

  call mpi_init(rc)
  call mpi_comm_size(mpi_comm_world, ntasks, rc)
  call mpi_comm_rank(mpi_comm_world, my_id, rc)

  localsize = datasize / ntasks
  allocate(localvector(localsize))

  do i = 1, localsize
      localvector(i) = i + my_id * localsize
  end do

  call mpi_type_size(MPI_INTEGER, dsize, rc)
  offset = my_id * localsize * dsize

  call mpi_file_open(MPI_COMM_WORLD, 'output.dat', &
                     MPI_MODE_CREATE+MPI_MODE_WRONLY, MPI_INFO_NULL, fh, rc)
  call mpi_file_write_at_all(fh, offset, localvector, localsize, &
                             MPI_INTEGER, MPI_STATUS_IGNORE, rc)
  call mpi_file_close(fh, rc)

  deallocate(localvector)
  call mpi_finalize(rc)

end program mpiio

max_stars_repo_name: mrgransky/advanced-parallel-prog

max_issues_repo_name: mrgransky/advanced-parallel-prog

max_forks_repo_name: mrgransky/advanced-parallel-prog


MODULE trabbl
  USE oce
  USE dom_oce
  USE phycst
  USE eosbn2
  USE trd_oce
  USE trdtra
  USE iom
  USE in_out_manager
  USE lbclnk
  USE prtctl
  USE timing
  USE lib_fortran
  IMPLICIT NONE
  PRIVATE
  PUBLIC :: tra_bbl
  PUBLIC :: tra_bbl_init
  PUBLIC :: tra_bbl_dif
  PUBLIC :: tra_bbl_adv
  PUBLIC :: bbl
  LOGICAL, PUBLIC :: ln_trabbl
  INTEGER, PUBLIC :: nn_bbl_ldf
  INTEGER, PUBLIC :: nn_bbl_adv
  REAL(KIND = wp), PUBLIC :: rn_ahtbbl
  REAL(KIND = wp), PUBLIC :: rn_gambbl
  LOGICAL, PUBLIC :: l_bbl
  REAL(KIND = wp), ALLOCATABLE, SAVE, DIMENSION(:, :), PUBLIC :: utr_bbl, vtr_bbl
  REAL(KIND = wp), ALLOCATABLE, SAVE, DIMENSION(:, :), PUBLIC :: ahu_bbl, ahv_bbl
  INTEGER, ALLOCATABLE, SAVE, DIMENSION(:, :), PUBLIC :: mbku_d, mbkv_d
  INTEGER, ALLOCATABLE, SAVE, DIMENSION(:, :), PUBLIC :: mgrhu, mgrhv
  REAL(KIND = wp), ALLOCATABLE, SAVE, DIMENSION(:, :) :: ahu_bbl_0, ahv_bbl_0
  REAL(KIND = wp), ALLOCATABLE, SAVE, DIMENSION(:, :), PUBLIC :: e3u_bbl_0, e3v_bbl_0
  CONTAINS
  INTEGER FUNCTION tra_bbl_alloc()
    ALLOCATE(utr_bbl(jpi, jpj), ahu_bbl(jpi, jpj), mbku_d(jpi, jpj), mgrhu(jpi, jpj), vtr_bbl(jpi, jpj), ahv_bbl(jpi, jpj), mbkv_d(jpi, jpj), mgrhv(jpi, jpj), ahu_bbl_0(jpi, jpj), ahv_bbl_0(jpi, jpj), e3u_bbl_0(jpi, jpj), e3v_bbl_0(jpi, jpj), STAT = tra_bbl_alloc)
    IF (lk_mpp) CALL mpp_sum(tra_bbl_alloc)
    IF (tra_bbl_alloc > 0) CALL ctl_warn('tra_bbl_alloc: allocation of arrays failed.')
  END FUNCTION tra_bbl_alloc
  SUBROUTINE tra_bbl(kt)
    INTEGER, INTENT( IN ) :: kt
    REAL(KIND = wp), ALLOCATABLE, DIMENSION(:, :, :) :: ztrdt, ztrds
    IF (ln_timing) CALL timing_start('tra_bbl')
    IF (l_trdtra) THEN
      ALLOCATE(ztrdt(jpi, jpj, jpk), ztrds(jpi, jpj, jpk))
      !$ACC KERNELS
      ztrdt(:, :, :) = tsa(:, :, :, jp_tem)
      ztrds(:, :, :) = tsa(:, :, :, jp_sal)
      !$ACC END KERNELS
    END IF
    IF (l_bbl) CALL bbl(kt, nit000, 'TRA')
    IF (nn_bbl_ldf == 1) THEN
      CALL tra_bbl_dif(tsb, tsa, jpts)
      IF (ln_ctl) CALL prt_ctl(tab3d_1 = tsa(:, :, :, jp_tem), clinfo1 = ' bbl_ldf  - Ta: ', mask1 = tmask, tab3d_2 = tsa(:, :, :, jp_sal), clinfo2 = ' Sa: ', mask2 = tmask, clinfo3 = 'tra')
      CALL lbc_lnk_multi(ahu_bbl, 'U', 1., ahv_bbl, 'V', 1.)
      CALL iom_put("ahu_bbl", ahu_bbl)
      CALL iom_put("ahv_bbl", ahv_bbl)
    END IF
    IF (nn_bbl_adv /= 0) THEN
      CALL tra_bbl_adv(tsb, tsa, jpts)
      IF (ln_ctl) CALL prt_ctl(tab3d_1 = tsa(:, :, :, jp_tem), clinfo1 = ' bbl_adv  - Ta: ', mask1 = tmask, tab3d_2 = tsa(:, :, :, jp_sal), clinfo2 = ' Sa: ', mask2 = tmask, clinfo3 = 'tra')
      CALL lbc_lnk_multi(utr_bbl, 'U', 1., vtr_bbl, 'V', 1.)
      CALL iom_put("uoce_bbl", utr_bbl)
      CALL iom_put("voce_bbl", vtr_bbl)
    END IF
    IF (l_trdtra) THEN
      !$ACC KERNELS
      ztrdt(:, :, :) = tsa(:, :, :, jp_tem) - ztrdt(:, :, :)
      ztrds(:, :, :) = tsa(:, :, :, jp_sal) - ztrds(:, :, :)
      !$ACC END KERNELS
      CALL trd_tra(kt, 'TRA', jp_tem, jptra_bbl, ztrdt)
      CALL trd_tra(kt, 'TRA', jp_sal, jptra_bbl, ztrds)
      DEALLOCATE(ztrdt, ztrds)
    END IF
    IF (ln_timing) CALL timing_stop('tra_bbl')
  END SUBROUTINE tra_bbl
  SUBROUTINE tra_bbl_dif(ptb, pta, kjpt)
    INTEGER, INTENT(IN   ) :: kjpt
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpk, kjpt), INTENT(IN   ) :: ptb
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpk, kjpt), INTENT(INOUT) :: pta
    INTEGER :: ji, jj, jn
    INTEGER :: ik
    REAL(KIND = wp) :: zbtr
    REAL(KIND = wp), DIMENSION(jpi, jpj) :: zptb
    !$ACC KERNELS
    DO jn = 1, kjpt
      DO jj = 1, jpj
        DO ji = 1, jpi
          ik = mbkt(ji, jj)
          zptb(ji, jj) = ptb(ji, jj, ik, jn)
        END DO
      END DO
      DO jj = 2, jpjm1
        DO ji = 2, jpim1
          ik = mbkt(ji, jj)
          pta(ji, jj, ik, jn) = pta(ji, jj, ik, jn) + (ahu_bbl(ji, jj) * (zptb(ji + 1, jj) - zptb(ji, jj)) - ahu_bbl(ji - 1, jj) * (zptb(ji, jj) - zptb(ji - 1, jj)) + ahv_bbl(ji, jj) * (zptb(ji, jj + 1) - zptb(ji, jj)) - ahv_bbl(ji, jj - 1) * (zptb(ji, jj) - zptb(ji, jj - 1))) * r1_e1e2t(ji, jj) / e3t_n(ji, jj, ik)
        END DO
      END DO
    END DO
    !$ACC END KERNELS
  END SUBROUTINE tra_bbl_dif
  SUBROUTINE tra_bbl_adv(ptb, pta, kjpt)
    INTEGER, INTENT(IN   ) :: kjpt
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpk, kjpt), INTENT(IN   ) :: ptb
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpk, kjpt), INTENT(INOUT) :: pta
    INTEGER :: ji, jj, jk, jn
    INTEGER :: iis, iid, ijs, ijd
    INTEGER :: ikus, ikud, ikvs, ikvd
    REAL(KIND = wp) :: zbtr, ztra
    REAL(KIND = wp) :: zu_bbl, zv_bbl
    DO jn = 1, kjpt
      DO jj = 1, jpjm1
        DO ji = 1, jpim1
          IF (utr_bbl(ji, jj) /= 0.E0) THEN
            !$ACC KERNELS
            iid = ji + MAX(0, mgrhu(ji, jj))
            iis = ji + 1 - MAX(0, mgrhu(ji, jj))
            ikud = mbku_d(ji, jj)
            ikus = mbku(ji, jj)
            zu_bbl = ABS(utr_bbl(ji, jj))
            zbtr = r1_e1e2t(iis, jj) / e3t_n(iis, jj, ikus)
            ztra = zu_bbl * (ptb(iid, jj, ikus, jn) - ptb(iis, jj, ikus, jn)) * zbtr
            pta(iis, jj, ikus, jn) = pta(iis, jj, ikus, jn) + ztra
            DO jk = ikus, ikud - 1
              zbtr = r1_e1e2t(iid, jj) / e3t_n(iid, jj, jk)
              ztra = zu_bbl * (ptb(iid, jj, jk + 1, jn) - ptb(iid, jj, jk, jn)) * zbtr
              pta(iid, jj, jk, jn) = pta(iid, jj, jk, jn) + ztra
            END DO
            zbtr = r1_e1e2t(iid, jj) / e3t_n(iid, jj, ikud)
            ztra = zu_bbl * (ptb(iis, jj, ikus, jn) - ptb(iid, jj, ikud, jn)) * zbtr
            pta(iid, jj, ikud, jn) = pta(iid, jj, ikud, jn) + ztra
            !$ACC END KERNELS
          END IF
          IF (vtr_bbl(ji, jj) /= 0.E0) THEN
            !$ACC KERNELS
            ijd = jj + MAX(0, mgrhv(ji, jj))
            ijs = jj + 1 - MAX(0, mgrhv(ji, jj))
            ikvd = mbkv_d(ji, jj)
            ikvs = mbkv(ji, jj)
            zv_bbl = ABS(vtr_bbl(ji, jj))
            zbtr = r1_e1e2t(ji, ijs) / e3t_n(ji, ijs, ikvs)
            ztra = zv_bbl * (ptb(ji, ijd, ikvs, jn) - ptb(ji, ijs, ikvs, jn)) * zbtr
            pta(ji, ijs, ikvs, jn) = pta(ji, ijs, ikvs, jn) + ztra
            DO jk = ikvs, ikvd - 1
              zbtr = r1_e1e2t(ji, ijd) / e3t_n(ji, ijd, jk)
              ztra = zv_bbl * (ptb(ji, ijd, jk + 1, jn) - ptb(ji, ijd, jk, jn)) * zbtr
              pta(ji, ijd, jk, jn) = pta(ji, ijd, jk, jn) + ztra
            END DO
            zbtr = r1_e1e2t(ji, ijd) / e3t_n(ji, ijd, ikvd)
            ztra = zv_bbl * (ptb(ji, ijs, ikvs, jn) - ptb(ji, ijd, ikvd, jn)) * zbtr
            pta(ji, ijd, ikvd, jn) = pta(ji, ijd, ikvd, jn) + ztra
            !$ACC END KERNELS
          END IF
        END DO
      END DO
    END DO
  END SUBROUTINE tra_bbl_adv
  SUBROUTINE bbl(kt, kit000, cdtype)
    INTEGER, INTENT(IN   ) :: kt
    INTEGER, INTENT(IN   ) :: kit000
    CHARACTER(LEN = 3), INTENT(IN   ) :: cdtype
    INTEGER :: ji, jj
    INTEGER :: ik
    INTEGER :: iis, iid, ikus, ikud
    INTEGER :: ijs, ijd, ikvs, ikvd
    REAL(KIND = wp) :: za, zb, zgdrho
    REAL(KIND = wp) :: zsign, zsigna, zgbbl
    REAL(KIND = wp), DIMENSION(jpi, jpj, jpts) :: zts, zab
    REAL(KIND = wp), DIMENSION(jpi, jpj) :: zub, zvb, zdep
    IF (kt == kit000) THEN
      IF (lwp) WRITE(numout, FMT = *)
      IF (lwp) WRITE(numout, FMT = *) 'trabbl:bbl : Compute bbl velocities and diffusive coefficients in ', cdtype
      IF (lwp) WRITE(numout, FMT = *) '~~~~~~~~~~'
    END IF
    !$ACC KERNELS
    DO jj = 1, jpj
      DO ji = 1, jpi
        ik = mbkt(ji, jj)
        zts(ji, jj, jp_tem) = tsb(ji, jj, ik, jp_tem)
        zts(ji, jj, jp_sal) = tsb(ji, jj, ik, jp_sal)
        zdep(ji, jj) = gdept_n(ji, jj, ik)
        zub(ji, jj) = un(ji, jj, mbku(ji, jj))
        zvb(ji, jj) = vn(ji, jj, mbkv(ji, jj))
      END DO
    END DO
    !$ACC END KERNELS
    CALL eos_rab(zts, zdep, zab)
    IF (nn_bbl_ldf == 1) THEN
      !$ACC KERNELS
      DO jj = 1, jpjm1
        DO ji = 1, jpim1
          za = zab(ji + 1, jj, jp_tem) + zab(ji, jj, jp_tem)
          zb = zab(ji + 1, jj, jp_sal) + zab(ji, jj, jp_sal)
          zgdrho = (za * (zts(ji + 1, jj, jp_tem) - zts(ji, jj, jp_tem)) - zb * (zts(ji + 1, jj, jp_sal) - zts(ji, jj, jp_sal))) * umask(ji, jj, 1)
          zsign = SIGN(0.5, - zgdrho * REAL(mgrhu(ji, jj)))
          ahu_bbl(ji, jj) = (0.5 - zsign) * ahu_bbl_0(ji, jj)
          za = zab(ji, jj + 1, jp_tem) + zab(ji, jj, jp_tem)
          zb = zab(ji, jj + 1, jp_sal) + zab(ji, jj, jp_sal)
          zgdrho = (za * (zts(ji, jj + 1, jp_tem) - zts(ji, jj, jp_tem)) - zb * (zts(ji, jj + 1, jp_sal) - zts(ji, jj, jp_sal))) * vmask(ji, jj, 1)
          zsign = SIGN(0.5, - zgdrho * REAL(mgrhv(ji, jj)))
          ahv_bbl(ji, jj) = (0.5 - zsign) * ahv_bbl_0(ji, jj)
        END DO
      END DO
      !$ACC END KERNELS
    END IF
    IF (nn_bbl_adv /= 0) THEN
      SELECT CASE (nn_bbl_adv)
      CASE (1)
        !$ACC KERNELS
        DO jj = 1, jpjm1
          DO ji = 1, jpim1
            za = zab(ji + 1, jj, jp_tem) + zab(ji, jj, jp_tem)
            zb = zab(ji + 1, jj, jp_sal) + zab(ji, jj, jp_sal)
            zgdrho = (za * (zts(ji + 1, jj, jp_tem) - zts(ji, jj, jp_tem)) - zb * (zts(ji + 1, jj, jp_sal) - zts(ji, jj, jp_sal))) * umask(ji, jj, 1)
            zsign = SIGN(0.5, - zgdrho * REAL(mgrhu(ji, jj)))
            zsigna = SIGN(0.5, zub(ji, jj) * REAL(mgrhu(ji, jj)))
            utr_bbl(ji, jj) = (0.5 + zsigna) * (0.5 - zsign) * e2u(ji, jj) * e3u_bbl_0(ji, jj) * zub(ji, jj)
            za = zab(ji, jj + 1, jp_tem) + zab(ji, jj, jp_tem)
            zb = zab(ji, jj + 1, jp_sal) + zab(ji, jj, jp_sal)
            zgdrho = (za * (zts(ji, jj + 1, jp_tem) - zts(ji, jj, jp_tem)) - zb * (zts(ji, jj + 1, jp_sal) - zts(ji, jj, jp_sal))) * vmask(ji, jj, 1)
            zsign = SIGN(0.5, - zgdrho * REAL(mgrhv(ji, jj)))
            zsigna = SIGN(0.5, zvb(ji, jj) * REAL(mgrhv(ji, jj)))
            vtr_bbl(ji, jj) = (0.5 + zsigna) * (0.5 - zsign) * e1v(ji, jj) * e3v_bbl_0(ji, jj) * zvb(ji, jj)
          END DO
        END DO
        !$ACC END KERNELS
      CASE (2)
        !$ACC KERNELS
        zgbbl = grav * rn_gambbl
        DO jj = 1, jpjm1
          DO ji = 1, jpim1
            iid = ji + MAX(0, mgrhu(ji, jj))
            iis = ji + 1 - MAX(0, mgrhu(ji, jj))
            ikud = mbku_d(ji, jj)
            ikus = mbku(ji, jj)
            za = zab(ji + 1, jj, jp_tem) + zab(ji, jj, jp_tem)
            zb = zab(ji + 1, jj, jp_sal) + zab(ji, jj, jp_sal)
            zgdrho = 0.5 * (za * (zts(iid, jj, jp_tem) - zts(iis, jj, jp_tem)) - zb * (zts(iid, jj, jp_sal) - zts(iis, jj, jp_sal))) * umask(ji, jj, 1)
            zgdrho = MAX(0.E0, zgdrho)
            utr_bbl(ji, jj) = e2u(ji, jj) * e3u_bbl_0(ji, jj) * zgbbl * zgdrho * REAL(mgrhu(ji, jj))
            ijd = jj + MAX(0, mgrhv(ji, jj))
            ijs = jj + 1 - MAX(0, mgrhv(ji, jj))
            ikvd = mbkv_d(ji, jj)
            ikvs = mbkv(ji, jj)
            za = zab(ji, jj + 1, jp_tem) + zab(ji, jj, jp_tem)
            zb = zab(ji, jj + 1, jp_sal) + zab(ji, jj, jp_sal)
            zgdrho = 0.5 * (za * (zts(ji, ijd, jp_tem) - zts(ji, ijs, jp_tem)) - zb * (zts(ji, ijd, jp_sal) - zts(ji, ijs, jp_sal))) * vmask(ji, jj, 1)
            zgdrho = MAX(0.E0, zgdrho)
            vtr_bbl(ji, jj) = e1v(ji, jj) * e3v_bbl_0(ji, jj) * zgbbl * zgdrho * REAL(mgrhv(ji, jj))
          END DO
        END DO
        !$ACC END KERNELS
      END SELECT
    END IF
  END SUBROUTINE bbl
  SUBROUTINE tra_bbl_init
    INTEGER :: ji, jj
    INTEGER :: ii0, ii1, ij0, ij1, ios
    REAL(KIND = wp), DIMENSION(jpi, jpj) :: zmbku, zmbkv
    NAMELIST /nambbl/ ln_trabbl, nn_bbl_ldf, nn_bbl_adv, rn_ahtbbl, rn_gambbl
    REWIND(UNIT = numnam_ref)
    READ(numnam_ref, nambbl, IOSTAT = ios, ERR = 901)
901 IF (ios /= 0) CALL ctl_nam(ios, 'nambbl in reference namelist', lwp)
    REWIND(UNIT = numnam_cfg)
    READ(numnam_cfg, nambbl, IOSTAT = ios, ERR = 902)
902 IF (ios > 0) CALL ctl_nam(ios, 'nambbl in configuration namelist', lwp)
    IF (lwm) WRITE(numond, nambbl)
    l_bbl = .TRUE.
    IF (lwp) THEN
      WRITE(numout, FMT = *)
      WRITE(numout, FMT = *) 'tra_bbl_init : bottom boundary layer initialisation'
      WRITE(numout, FMT = *) '~~~~~~~~~~~~'
      WRITE(numout, FMT = *) '       Namelist nambbl : set bbl parameters'
      WRITE(numout, FMT = *) '          bottom boundary layer flag          ln_trabbl  = ', ln_trabbl
    END IF
    IF (.NOT. ln_trabbl) RETURN
    IF (lwp) THEN
      WRITE(numout, FMT = *) '          diffusive bbl (=1)   or not (=0)    nn_bbl_ldf = ', nn_bbl_ldf
      WRITE(numout, FMT = *) '          advective bbl (=1/2) or not (=0)    nn_bbl_adv = ', nn_bbl_adv
      WRITE(numout, FMT = *) '          diffusive bbl coefficient           rn_ahtbbl  = ', rn_ahtbbl, ' m2/s'
      WRITE(numout, FMT = *) '          advective bbl coefficient           rn_gambbl  = ', rn_gambbl, ' s'
    END IF
    IF (tra_bbl_alloc() /= 0) CALL ctl_stop('STOP', 'tra_bbl_init : unable to allocate arrays')
    IF (lwp) THEN
      IF (nn_bbl_adv == 1) WRITE(numout, FMT = *) '       * Advective BBL using upper velocity'
      IF (nn_bbl_adv == 2) WRITE(numout, FMT = *) '       * Advective BBL using velocity = F( delta rho)'
    END IF
    !$ACC KERNELS
    DO jj = 1, jpjm1
      DO ji = 1, jpim1
        mbku_d(ji, jj) = MAX(mbkt(ji + 1, jj), mbkt(ji, jj))
        mbkv_d(ji, jj) = MAX(mbkt(ji, jj + 1), mbkt(ji, jj))
      END DO
    END DO
    zmbku(:, :) = REAL(mbku_d(:, :), wp)
    zmbkv(:, :) = REAL(mbkv_d(:, :), wp)
    !$ACC END KERNELS
    CALL lbc_lnk_multi(zmbku, 'U', 1., zmbkv, 'V', 1.)
    !$ACC KERNELS
    mbku_d(:, :) = MAX(INT(zmbku(:, :)), 1)
    mbkv_d(:, :) = MAX(NINT(zmbkv(:, :)), 1)
    mgrhu(:, :) = 0
    mgrhv(:, :) = 0
    DO jj = 1, jpjm1
      DO ji = 1, jpim1
        IF (gdept_0(ji + 1, jj, mbkt(ji + 1, jj)) - gdept_0(ji, jj, mbkt(ji, jj)) /= 0._wp) THEN
          mgrhu(ji, jj) = INT(SIGN(1.E0, gdept_0(ji + 1, jj, mbkt(ji + 1, jj)) - gdept_0(ji, jj, mbkt(ji, jj))))
        END IF
        IF (gdept_0(ji, jj + 1, mbkt(ji, jj + 1)) - gdept_0(ji, jj, mbkt(ji, jj)) /= 0._wp) THEN
          mgrhv(ji, jj) = INT(SIGN(1.E0, gdept_0(ji, jj + 1, mbkt(ji, jj + 1)) - gdept_0(ji, jj, mbkt(ji, jj))))
        END IF
      END DO
    END DO
    DO jj = 1, jpjm1
      DO ji = 1, jpim1
        e3u_bbl_0(ji, jj) = MIN(e3u_0(ji, jj, mbkt(ji + 1, jj)), e3u_0(ji, jj, mbkt(ji, jj)))
        e3v_bbl_0(ji, jj) = MIN(e3v_0(ji, jj, mbkt(ji, jj + 1)), e3v_0(ji, jj, mbkt(ji, jj)))
      END DO
    END DO
    !$ACC END KERNELS
    CALL lbc_lnk_multi(e3u_bbl_0, 'U', 1., e3v_bbl_0, 'V', 1.)
    !$ACC KERNELS
    ahu_bbl_0(:, :) = rn_ahtbbl * e2_e1u(:, :) * e3u_bbl_0(:, :) * umask(:, :, 1)
    ahv_bbl_0(:, :) = rn_ahtbbl * e1_e2v(:, :) * e3v_bbl_0(:, :) * vmask(:, :, 1)
    !$ACC END KERNELS
  END SUBROUTINE tra_bbl_init
END MODULE trabbl


text:
  subroutine fermiBreakUpInterface (gsmObj, residINC, gsmRxn)

! ======================================================================
!
!    Fermi break-up calculation of nuclei with A<13 in Preco and Evap
!
!    Called from PRECOF
!
!    Written by K.K. Gudima, 06/23/06
!    Modified by SGM, 07/09/06
!    Edited by AJS, LANL T-2, December, 2011.
!    Edited by LMK, XCP-3, July 2013 (included error protection)
!
! ======================================================================

    use, intrinsic :: iso_fortran_env, only: int32, real64, int64
    use gsm_params, only: zro, one, thousand, twpi
    use fermiBreakupClass, only : &
         & fermiBreakUpResults, newFermiBreakUpResults, &
         & fermiBreakUpProgeny, fermiBreakUpNucleus

    implicit none
    class(GSM),   intent(inout) :: gsmObj
    class(GSMResidual), intent(in   ) :: residINC
    class(GSMReaction), intent(inout) :: gsmRxn

    integer(int32) :: iafr, ifr, izfr, jfr
    real(real64)   :: cosPhi, cosTheta, phi, totLinMome, &
         & sinPhi, sinTheta, theta, temp

    ! Results object:
    type(fermiBreakUpProgeny), dimension( nint(residINC%numBaryons) ) :: progenyBnk
    type(fermiBreakUpResults) :: fbuResults
    ! Nucleus type:
    type(fermiBreakUpNucleus) :: residual

    ! For photon emission:
    type(GSMResidual) :: photoResidual
    type(GSMResults), pointer :: results => NULL()

! ======================================================================

    ! Set up variables and nucleus:
    residual%numBaryons = residINC%numBaryons
    residual%numProtons = residINC%numProtons
    residual%kinEnergy  = residINC%kinEnergy/thousand
    residual%linearXMom = residINC%linearMom(1)
    residual%linearYMom = residINC%linearMom(2)
    residual%linearZMom = residINC%linearMom(3)
    gsmRxn%outData%ifermi = gsmRxn%outData%ifermi + 1

    ! Create a object for photon emission, if used:
    if ( gsmObj%usePhotonEmission ) then
       photoResidual%numBaryons    = residINC%numBaryons
       photoResidual%numProtons    = residINC%numProtons
       photoResidual%kinEnergy     = residINC%kinEnergy
       photoResidual%linearMom(1)  = residINC%linearMom(1)
       photoResidual%linearMom(2)  = residINC%linearMom(2)
       photoResidual%linearMom(3)  = residINC%linearMom(3)
    end if

    ! Construct the results object:
    fbuResults = newFermiBreakUpResults( progenyBnk )

    ! Perform Fermi Break-Up simulation
    call gsmObj%genModels%fbu%execute (residual, fbuResults)

    ! Now interface results:
    results => gsmRxn%results

    ! Set evaporation/fission flags
    results%info%fusion = zro
    results%info%wf = zro

    ! If fragments created, store data into secondary arrays of program
    progenyLoop: do ifr = 1,fbuResults%numProgeny

       ! Check for available particle storage:
       if ( results%numProgeny > results%maxProgenyM1 ) then
          write(gsmObj%io%message, 1000) fbuResults%numProgeny - ifr + 1
          call gsmObj%io%print(3, 3, gsmObj%io%message)
          exit progenyLoop
       end if

       ! Obtain interim results:
       totLinMome = sqrt( progenyBnk(ifr)%linearXMom**2 + progenyBnk(ifr)%linearYMom**2 + &
            & progenyBnk(ifr)%linearZMom**2 )
       cosTheta = progenyBnk(ifr)%linearZMom/totLinMome
       if ( abs(cosTheta) > one ) cosTheta = sign(one, cosTheta) ! Ensure cos(theta) in limits of [-1, 1]
       sinTheta = sqrt(one - cosTheta**2)
       if (sinTheta > zro) then
          ! Particle has forward movement (get phi component of direction)
          temp = totLinMome*sinTheta
          cosPhi = progenyBnk(ifr)%linearXMom/(temp) ! cos(phi) = px/[sin(theta) ptot]
          sinPhi = progenyBnk(ifr)%linearYMom/(temp) ! sin(phi) = py/[sin(theta) ptot]
       else
          ! Scatter is at 90 degrees (in theta) to incident beam
          cosPhi = one
          sinPhi = zro
       endif
       ! Obtain (phi), ensure within range [0, 2pi]
       phi = atan2 (sinPhi, cosPhi)
       if (phi < zro) phi = twpi + phi
       theta = atan2 (sinTheta, cosTheta)

       ! Obtain GSM fragment label:
       iafr = nint(progenyBnk(ifr)%numBaryons)
       izfr = nint(progenyBnk(ifr)%numProtons)
       jfr = zro
       if (progenyBnk(ifr)%numBaryons <= 4.1d0 .and. progenyBnk(ifr)%numProtons <= 2.1d0) then
          if (iafr == 1 .and. izfr == 0)  jfr = 1 ! n
          if (iafr == 1 .and. izfr == 1)  jfr = 2 ! p
          if (iafr == 2 .and. izfr == 1)  jfr = 3 ! d
          if (iafr == 3 .and. izfr == 1)  jfr = 4 ! t
          if (iafr == 3 .and. izfr == 2)  jfr = 5 ! He-3
          if (iafr == 4 .and. izfr == 2)  jfr = 6 ! He-4
       else
          jfr = thousand*izfr + (iafr - izfr) ! Type for all others
       endif
       if (jfr == 0) jfr = thousand*izfr + (iafr - izfr)

       ! Store progeny in particle bank:
       results%numProgeny = results%numProgeny + 1
       results%progenyBnk(results%numProgeny)%numBaryons = progenyBnk(ifr)%numBaryons
       results%progenyBnk(results%numProgeny)%numProtons = progenyBnk(ifr)%numProtons
       results%progenyBnk(results%numProgeny)%kinEnergy  = progenyBnk(ifr)%kinEnergy / thousand
       results%progenyBnk(results%numProgeny)%restMass   = progenyBnk(ifr)%restMass/thousand
       results%progenyBnk(results%numProgeny)%phi        = phi
       results%progenyBnk(results%numProgeny)%theta      = theta
       results%progenyBnk(results%numProgeny)%sinTheta   = sinTheta
       results%progenyBnk(results%numProgeny)%cosTheta   = cosTheta
       results%progenyBnk(results%numProgeny)%typeID     = jfr
       results%progenyBnk(results%numProgeny)%prodMech   = 1500

       ! Simulate photon emission:
       if ( gsmObj%usePhotonEmission ) then
          call gammaCascade(results%progenyBnk(results%numProgeny), photoResidual)
       end if

    end do progenyLoop

    return
! ======================================================================
1000 format("The GSM progeny array was exceeded. Cannot tally last ", i3, &
          & " fragments.")
! ======================================================================
  end subroutine fermiBreakUpInterface

max_stars_repo_name: lanl/generalized-spallation-model

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_name: lanl/generalized-spallation-model

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_name: lanl/generalized-spallation-model

max_forks_repo_licenses: Intel Unlicense


! { dg-do compile }
! { dg-options "-finit-local-zero -finit-derived -fdump-tree-original" }
!
! Make sure -finit-derived initializes components of local derived type
! variables to zero with -finit-local-zero.
!

subroutine dummy(i1,r1,c1,l1,i2,r2,c2,l2)
  implicit none
  integer, intent(in) :: i1
  real, intent(in) :: r1
  character, intent(in) :: c1
  logical, intent(in) :: l1
  integer, intent(out) :: i2
  real, intent(out) :: r2
  character, intent(out) :: c2
  logical, intent(out) :: l2
end subroutine

type t2
  integer i2
  real r2
  character c2
  logical l2
end type

type t1
  logical l1
  real r1
  character c1
  integer i1
  type (t2) y
end type

type (t1) :: x

call dummy (x%i1, x%r1, x%c1, x%l1, x%y%i2, x%y%r2, x%y%c2, x%y%l2)

end

! We expect to see each component initialized exactly once in MAIN.
! NB. the "once" qualifier also tests that the dummy variables aren't
! given an extraneous initializer.
! { dg-final { scan-tree-dump-times "i1= *0" 1 "original" } }
! { dg-final { scan-tree-dump-times "r1= *0" 1 "original" } }
! { dg-final { scan-tree-dump-times "c1= *\"\"" 1 "original" } }
! { dg-final { scan-tree-dump-times "l1= *0" 1 "original" } }
! { dg-final { scan-tree-dump-times "i2= *0" 1 "original" } }
! { dg-final { scan-tree-dump-times "r2= *0" 1 "original" } }
! { dg-final { scan-tree-dump-times "c2= *\"\"" 1 "original" } }
! { dg-final { scan-tree-dump-times "l2= *0" 1 "original" } }


text:
      module mg_param
c
c     ******************************************************************
c     *                                                                *
c     *   parameters controlling the multigrid process                 *
c     *                                                                *
c     ******************************************************************
c
      integer   :: kode,mode,lcyc
      integer   :: k1,kx,kw

      real      :: fcoll,fadd,fbc
      real      :: cflf,cflc,cfl0,hmf,hmc
      real      :: epsf,epsc

      end module mg_param

ext: f

max_stars_repo_name: AlexT-L/RANS

max_issues_repo_name: AlexT-L/RANS

max_forks_repo_name: AlexT-L/RANS


text:
../tb/utils.sv
../tb/log_shift_tester.sv
../rtl/log_shifter.sv


text:
!------------------------------------------------------------------------------
!P+
! NAME:
!       SRF_netCDF_IO_Test
!
! PURPOSE:
!       Program to test the SRF netCDF format I/O functions.
!
! CATEGORY:
!       Instrument Information : SRF
!
! LANGUAGE:
!       Fortran-95
!
! MODULES:
!       Type_Kinds:                 Module containing definitions for kinds
!                                   of variable types.
!
!       Message_Handler:            Module to define simple error codes and
!                                   handle error conditions
!                                   USEs: FILE_UTILITY module
!
!       SRF_Define:                 Module defining the SRF data structure
!                                   and containing routines to manipulate it.
!                                   USEs: TYPE_KINDS module
!                                         Message_Handler module
!
!       SRF_netCDF_IO:              Module containing routines to read and
!                                   write SRF netCDF format files.
!                                   USEs: TYPE_KINDS module
!                                         Message_Handler module
!                                         SRF_DEFINE module
!                                         NETCDF module
!                                         NETCDF_UTILITY module
!
! CONTAINS:
!       None.
!
! INCLUDE FILES:
!       None.
!
! EXTERNALS:
!       None.
!
! COMMON BLOCKS:
!       None.
!
! FILES ACCESSED:
!       netCDF SRF data file
!
! SIDE EFFECTS:
!       None.
!
! RESTRICTIONS:
!       None.
!
! CREATION HISTORY:
!       Written by:     Paul van Delst, CIMSS/SSEC 16-Apr-2002
!                       paul.vandelst@ssec.wisc.edu
!
!  Copyright (C) 2002 Paul van Delst
!
!P-
!------------------------------------------------------------------------------

PROGRAM SRF_netCDF_IO_Test

  ! ------------
  ! Module usage
  ! ------------

  USE Type_Kinds
  USE Message_Handler

  USE SRF_Define
  USE SRF_netCDF_IO

  ! -----------------------
  ! Disable implicit typing
  ! -----------------------

  IMPLICIT NONE

  ! ----------
  ! Parameters
  ! ----------

  CHARACTER( * ), PARAMETER :: PROGRAM_NAME = 'SRF_netCDF_IO_Test'
  CHARACTER( * ), PARAMETER :: PROGRAM_RCS_ID = &
  CHARACTER( * ), PARAMETER :: PROGRAM_HEADER = &
  '**********************************************************'

  CHARACTER( * ), PARAMETER :: OUTPUT_SRF_FILENAME = 'Test.srf.nc'

  ! ---------
  ! Variables
  ! ---------

  CHARACTER( 256 ) :: Message

  INTEGER         :: pn_pos
  CHARACTER( 80 ) :: pn_fmt

  CHARACTER( 256 ) :: SRF_Filename

  INTEGER :: Error_Status
  INTEGER :: Allocate_Status

  ! -- SRF dimensions
  INTEGER :: n_Channels, l

  ! -- SRF global attributes
  CHARACTER( 256 ) :: Title
  CHARACTER( 256 ) :: History
  CHARACTER( 256 ) :: Sensor_Name
  CHARACTER( 256 ) :: Platform_Name
  CHARACTER( 256 ) :: Comment

  ! -- SRF sensor IDs
  INTEGER :: NCEP_Sensor_ID
  INTEGER :: WMO_Satellite_ID
  INTEGER :: WMO_Sensor_ID

  ! -- SRF channel list data array
  INTEGER, DIMENSION( : ), ALLOCATABLE :: Channel_List

  ! -- SRF data structure
  TYPE( SRF_type ) :: SRF

  !#----------------------------------------------------------------------------#
  !#                       -- OUTPUT DESCRIPTIVE HEADER --                      #
  !#----------------------------------------------------------------------------#

  pn_pos = ( LEN( PROGRAM_HEADER ) / 2 ) - &
           ( LEN( PROGRAM_NAME ) / 2 )
  pn_pos = MAX( pn_pos, 0 ) + 5
  WRITE( pn_fmt, '( "( ",i2,"x, a )" )' ) pn_pos

  WRITE( *, '(/5x, a)' ) PROGRAM_HEADER
  WRITE( *, FMT = TRIM( pn_fmt ) ) PROGRAM_NAME
  WRITE( *, * )
  WRITE( *, '( 5x, " Program to test read netCDF SRF data files.")' )
  WRITE( *, * )
  WRITE( *, '( 5x, " $Revision: 1.7 $")' )
  WRITE( *, '( 5x, a)' ) PROGRAM_HEADER

  !#----------------------------------------------------------------------------#
  !#                     -- ENTER A NETCDF SRF FILENAME --                      #
  !#----------------------------------------------------------------------------#

  WRITE( *, FMT     = '( /5x, "Enter a netCDF SRF data file : " )', &
            ADVANCE = 'NO' )
  READ( *, FMT = '( a )' ) SRF_Filename
  SRF_Filename = ADJUSTL( SRF_Filename )

  !#----------------------------------------------------------------------------#
  !#                        -- INQUIRE THE SRF FILE --                          #
  !#----------------------------------------------------------------------------#

  ! ---------------------------------------------
  ! Get the file dimensions and global attributes
  ! ---------------------------------------------

  Error_Status = Inquire_SRF_netCDF( TRIM( SRF_Filename ), &
                                     n_Channels       = n_Channels, &
                                     NCEP_Sensor_ID   = NCEP_Sensor_ID, &
                                     WMO_Satellite_ID = WMO_Satellite_ID, &
                                     WMO_Sensor_ID    = WMO_Sensor_ID, &
                                     Title            = Title, &
                                     History          = History, &
                                     Sensor_Name      = Sensor_Name, &
                                     Platform_Name    = Platform_Name, &
                                     Comment          = Comment )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL display_message( PROGRAM_NAME, &
                          'Error inquiring the netCDF SRF file '//&
                          TRIM( SRF_Filename ), &
                          Error_Status )
    STOP
  END IF

  ! ----------------
  ! Output some info
  ! ----------------

  WRITE( *, '( /5x, "SRF data set attributes:" )' )
  WRITE( *, '(  5x, "TITLE         : ", a )' ) TRIM( Title         )
  WRITE( *, '(  5x, "HISTORY       : ", a )' ) TRIM( History       )
  WRITE( *, '(  5x, "SENSOR_NAME   : ", a )' ) TRIM( Sensor_Name   )
  WRITE( *, '(  5x, "PLATFORM_NAME : ", a )' ) TRIM( Platform_Name )
  WRITE( *, '(  5x, "COMMENT       : ", a )' ) TRIM( Comment       )

  ! -------------------------------
  ! Allocate the channel list array
  ! -------------------------------

  ALLOCATE( Channel_List( n_Channels ), &
            STAT = Allocate_Status )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL display_message( PROGRAM_NAME, &
                          'Error allocating CHANNEL_LIST or SRF arrays.', &
                          Error_Status )
    STOP
  END IF

  ! ---------------------
  ! Read the channel list
  ! ---------------------

  Error_Status = Inquire_SRF_netCDF( TRIM( SRF_Filename ), &
                                     Channel_List = Channel_List )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL display_message( PROGRAM_NAME, &
                          'Error reading CHANNEL_LIST from netCDF SRF file '//&
                          TRIM( SRF_Filename ), &
                          Error_Status )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#                      -- CREATE A FILE FOR OUTPUT --                        #
  !#----------------------------------------------------------------------------#

  Error_Status = Create_SRF_netCDF( OUTPUT_SRF_FILENAME, &
                                    Channel_List, &
                                    NCEP_Sensor_ID   = NCEP_Sensor_ID, &
                                    WMO_Satellite_ID = WMO_Satellite_ID, &
                                    WMO_Sensor_ID    = WMO_Sensor_ID, &
                                    Title            = TRIM( Title ), &
                                    History          = PROGRAM_RCS_ID//';'//TRIM( History ), &
                                    Sensor_Name      = TRIM( Sensor_Name ), &
                                    Platform_Name    = TRIM( Platform_Name ), &
                                    Comment          = TRIM( Comment ) )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL Display_Message( PROGRAM_NAME, &
                          'Error creating '//TRIM( OUTPUT_SRF_FILENAME ), &
                          Error_Status )
    STOP
  END IF

  !#----------------------------------------------------------------------------#
  !#                         -- LOOP OVER CHANNELS --                           #
  !#----------------------------------------------------------------------------#

  WRITE( *, '( /5x, "Channel",4x,"n_Points",6x,"Frequency limits" )' )

  Channel_Read_Loop: DO l = 1, n_Channels

    ! ---------------------------------
    ! Read the current SRF channel data
    ! ---------------------------------

    Error_Status = Read_SRF_netCDF( TRIM( SRF_Filename ), &  ! Input
                                    Channel_List( l ),    &  ! Input
                                    SRF                   )  ! Output

    IF ( Error_Status /= SUCCESS ) THEN
      WRITE( Message, '( "Error reading channel #", i5, " SRF from ", a )' ) &
                      Channel_List( l ), TRIM( SRF_Filename )
      CALL display_message( PROGRAM_NAME, &
                            TRIM( Message ), &
                            Error_Status )
      STOP
    END IF

    ! ----------------------------------------------
    ! Ouput some information about the SRF data read
    ! ----------------------------------------------

    WRITE( *, '( 5x, i5, 5x, i8, 5x,  "[",f9.4,",",f9.4,"]" )' ) &
              SRF%Channel, SRF%n_Points, &
              SRF%Begin_Frequency, SRF%End_Frequency

    IF ( MOD( l, 30 ) == 0 ) THEN
      WRITE( *, '( 5x, "Press <ENTER> to continue..." )' )
      READ( *, * )
    END IF

    ! ---------------------------------------
    ! Write the SRF to the output netCDF file
    ! ---------------------------------------

    Error_Status = Write_SRF_netCDF( OUTPUT_SRF_FILENAME, &
                                     SRF )

    IF ( Error_Status /= SUCCESS ) THEN
      WRITE( Message, '( "Error writing channel #", i5, " SRF to ", a )' ) &
                      Channel_List( l ), OUTPUT_SRF_FILENAME
      CALL display_message( PROGRAM_NAME, &
                            TRIM( Message ), &
                            Error_Status )
      STOP
    END IF

    ! ------------------------------------------------
    ! Destroy the SRF data structure for the next read
    ! ------------------------------------------------

    Error_Status = Destroy_SRF( SRF )

    IF ( Error_Status /= SUCCESS ) THEN
      WRITE( Message, '( "Error occurred destroying SRF structure ", &
                        &"after channel #", i4, " read." )' ) &
                      Channel_List( l )
      CALL display_message( PROGRAM_NAME,    &
                            TRIM( Message ), &
                            Error_Status     )
      STOP
    END IF

  END DO Channel_Read_Loop

  !#----------------------------------------------------------------------------#
  !#                    -- DEALLOCATE THE CHANNEL LIST --                       #
  !#----------------------------------------------------------------------------#

  DEALLOCATE( Channel_List, &
              STAT = Allocate_Status )

  IF ( Error_Status /= SUCCESS ) THEN
    CALL display_message( PROGRAM_NAME, &
                          'Error deallocating CHANNEL_LIST array after netCDF test read.', &
                          WARNING )
  END IF

END PROGRAM SRF_netCDF_IO_Test

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm


      INTEGER N_HK, N_RT, N_TREH
      REAL, DIMENSION(1,6) :: HK_RANGE
      REAL, DIMENSION(1,5) :: RT_RANGE
      REAL, DIMENSION(1,3) :: TREH_RANGE
      REAL, DIMENSION(6,5,3) :: CEC_TABLE

      N_HK   = 6
      N_RT   = 5
      N_TREH = 3

      HK_RANGE(1,1) = 7.0E0/5.0E0
      HK_RANGE(1,2) = 8.0E0/5.0E0
      HK_RANGE(1,3) = 9.0E0/5.0E0
      HK_RANGE(1,4) = 2.0E0
      HK_RANGE(1,5) = 1.1E1/5.0E0
      HK_RANGE(1,6) = 1.2E1/5.0E0

      RT_RANGE(1,1) = 8.0E2
      RT_RANGE(1,2) = 1.6E3
      RT_RANGE(1,3) = 3.2E3
      RT_RANGE(1,4) = 6.4E3
      RT_RANGE(1,5) = 1.28E4

      TREH_RANGE(1,2) = 5.0E1
      TREH_RANGE(1,3) = 1.0E2

      CEC_TABLE(1,1,2) = 2.872680733645279E-1
      CEC_TABLE(1,1,3) = 7.749716029388073E-1
      CEC_TABLE(1,2,2) = 2.039644535775086E-1
      CEC_TABLE(1,2,3) = 5.924318712770024E-1
      CEC_TABLE(1,3,2) = 1.443700895511554E-1
      CEC_TABLE(1,3,3) = 4.633605499494066E-1
      CEC_TABLE(1,4,2) = 1.004911816583118E-1
      CEC_TABLE(1,4,3) = 3.699178360235265E-1
      CEC_TABLE(1,5,2) = 6.733821510084832E-2
      CEC_TABLE(1,5,3) = 3.001150710959782E-1
      CEC_TABLE(2,1,2) = 2.15679590997346E-1
      CEC_TABLE(2,1,3) = 7.032901821333088E-1
      CEC_TABLE(2,2,2) = 1.38804898233508E-1
      CEC_TABLE(2,2,3) = 5.190092839824512E-1
      CEC_TABLE(2,3,2) = 8.540581665774229E-2
      CEC_TABLE(2,3,3) = 3.955784955811628E-1
      CEC_TABLE(2,4,2) = 4.687087515554092E-2
      CEC_TABLE(2,4,3) = 3.083346385006999E-1
      CEC_TABLE(2,5,2) = 1.823657536308939E-2
      CEC_TABLE(2,5,3) = 2.440957040061711E-1
      CEC_TABLE(3,1,2) = 1.351465804725244E-1
      CEC_TABLE(3,1,3) = 5.944539142560507E-1
      CEC_TABLE(3,2,2) = 6.948210676888233E-2
      CEC_TABLE(3,2,3) = 4.287398277454012E-1
      CEC_TABLE(3,3,2) = 2.446511411233687E-2
      CEC_TABLE(3,3,3) = 3.186820261578073E-1
      CEC_TABLE(3,4,3) = 2.413596674968999E-1
      CEC_TABLE(3,5,3) = 1.847296324018285E-1
      CEC_TABLE(4,1,2) = 5.497866179885413E-2
      CEC_TABLE(4,1,3) = 4.803653864888552E-1
      CEC_TABLE(4,2,2) = 1.021202577184958E-3
      CEC_TABLE(4,2,3) = 3.363622877434947E-1
      CEC_TABLE(4,3,3) = 2.409694933687837E-1
      CEC_TABLE(4,4,3) = 1.741794817845234E-1
      CEC_TABLE(4,5,3) = 1.254740804951762E-1
      CEC_TABLE(5,1,3) = 3.674090874850944E-1
      CEC_TABLE(5,2,3) = 2.450395028559258E-1
      CEC_TABLE(5,3,3) = 1.641449222712615E-1
      CEC_TABLE(5,4,3) = 1.077346391909527E-1
      CEC_TABLE(5,5,3) = 6.682986734578829E-2
      CEC_TABLE(6,1,3) = 2.566658710465144E-1
      CEC_TABLE(6,2,3) = 1.54999205732387E-1
      CEC_TABLE(6,3,3) = 8.805729195880023E-2
      CEC_TABLE(6,4,3) = 4.168825812144024E-2
      CEC_TABLE(6,5,3) = 8.364384007724771E-3


text:
!*************************************************************************
subroutine fan_gmres(nd)
  use global_com,only: dp,nitermax,precis
  use global_dim
  implicit none
  integer,parameter :: m = 30
  integer::itertotal,iterout,n,nd,k
  complex(kind=dp)::cx(nd)
  complex(kind=dp)::cb(m+1)
  complex(kind=dp)::cy(nd,m+1),ch(m+1,m)
  complex(kind=dp)::cs(m)
  real(kind=dp)::rc(m)
  complex(kind=dp)::ctemp
  real(kind=dp)::ay0, be, bea
  real(kind=dp)::xtime,xtim
  real(kind=dp),external:: scnrm22

  print *, '.......in fan_gmres algorithm ......'
  open(unit=161,file='con_fangmres.out')
  PRINT*, "Starting iteration of solution..."
  call cpu_time(xtime)

  cx    = (0.0_dp,0.0_dp)
  ay0   = scnrm22(nd,crhs)
  itertotal = 0

!print*,"nd",nd

!************************************************************************
  do iterout = 1, nitermax
!************************************************************************
     itertotal = itertotal +1
!====================================================================
     call            convol(cx(1:nd),cy(1:nd,m+1),nd)
     cy(1:nd,m+1)  = crhs(1:nd) -cy(1:nd,m+1)
     cb(1)         = scnrm22(nd,cy(1:nd,m+1))
     be            = cb(1)/ay0
!====================================================================
     !write(*,*) 'restart', iterout,itertotal,'  truebe',be
     !write(161,*) itertotal,be
     if(be.le.precis.or.itertotal>nitermax) exit
     cy(1:nd,1) =  cy(1:nd,m+1)/cb(1)
!************************************************************************
     do n = 1, m
!************************************************************************
        itertotal = itertotal +1
!========================step1=======================================
        call convol(cy(1:nd,n),cy(1:nd,n+1),nd)
!========================step2=======================================
  	do k=1,n
           ch(k,n) = dot_product(cy(1:nd,k),cy(1:nd,n+1))
           cy(1:nd,n+1) = cy(1:nd,n+1) - cy(1:nd,k)*ch(k,n)
        enddo
        ch(n+1,n)  =  scnrm22(nd,cy(1:nd,n+1))
        if(n.lt.m)  cy(1:nd,n+1) = cy(1:nd,n+1)/ch(n+1,n)
        if(n.ne.1) then
           do k=1,n-1
              ctemp       = rc(k)*ch(k+1,n) -        cs(k) *ch(k,n)
              ch(k,n)     = rc(k)*ch(k,n)   + conjg (cs(k))*ch(k+1,n)
              ch(k+1,n)   = ctemp
           enddo
        endif
!========================step3=======================================
        call givens_loc(ch(n,n),ch(n+1,n),rc(n),cs(n))
!========================step4=======================================
        cb(n+1) =  -cs(n)*cb(n)
        cb(n)   =   rc(n)*cb(n)
        ch(n,n)   = rc(n)*ch(n,n)+ conjg (cs(n))*ch(n+1,n)
        ch(n+1,n) = (0.0_dp,0.0_dp)
!========================step5=======================================
        bea = abs(cb(n+1))/ay0
        !write(*,*) itertotal,bea
        !write(161,*) itertotal,bea
        if(n.eq.m.or.bea.le.precis) then
           call cutrisub(n, ch(1:n,1:n), cb(1:n), cb(1:n))
           cy(1:nd,m+1) = matmul(cy(1:nd,1:n),cb(1:n))
           cx = cx + cy(1:nd,m+1)
           exit
        endif
!===================================================================
     enddo   !n }
  enddo !iterout}
!===================================================================
  if(be>precis) then
     write( *,*) 'convergence is not achieved'
     write(161,*) 'convergence is not achieved'
  else
     write( *,*) 'convergence achieved'
     write(161,*) 'convergence achieved'
  endif
  call cpu_time(xtim)
  write(*, *) "gmres itr time...", xtim-xtime
  write(161,*) "gmres itr time...", xtim-xtime

!------------------------solution------------------------
  rj = cx

  close(161)

  return
end subroutine fan_gmres
!*************************************************************************

!*************************************************************************
subroutine givens_loc(z1,z2,c,s)
  use global_com,only: dp
  implicit none
  complex(kind=dp)::  z1,z2,s
  real(kind=dp)::      c,vnormz

  vnormz = sqrt ( abs(z1)**2 + abs(z2)**2)
  if (abs(z1).ne.0.) then
     c =     abs(z1)/vnormz
     s =     z2/z1* c
  else if (abs(z2).ne.0.) then
     c = 0.
     s = z2/ abs(z2)
  else
     c = 1.
     s = 0.
  endif

end subroutine givens_loc
!*************************************************************************

!*************************************************************************
function  scnrm22(nd,cy9)
  use global_com,only: dp
  implicit none
  integer nd,i
  complex(kind=dp)::  cy9(nd)
  real(kind=dp):: scnrm22

  scnrm22=0.0_dp
  do i=1,nd
     scnrm22= scnrm22+( conjg( cy9(i) )*cy9(i))
  end do
  scnrm22=sqrt(scnrm22)
  !	write(*,*) scnrm22
  !	pause

  return
end function scnrm22
!*************************************************************************

!*************************************************************************
subroutine cutrisub(n, cma, cvx, cvy)
  use global_com,only: dp
  implicit none
  integer n
  complex(kind=dp):: cma(n, n), cvx(n), cvy(n)
  integer i, j

  do i = n, 1, - 1
     cvy(i) = cvx(i)
     do j = i + 1, n
        cvy(i) = cvy(i) - cma(i, j) * cvy(j)
     enddo
     cvy(i) = cvy(i) / cma(i, i)
  enddo

end subroutine cutrisub
!*************************************************************************

!*************************************************************************
subroutine convol(cin,cout,nd)  !!!convolnear
  use global_com,only: dp,ntriangle,nipp
  use global_dim,only: zzsparse
  implicit none
  integer::nd,i,kk1,kk2
  complex(kind=dp):: cin(nd),cout(nd)

  cout  = (0.0_dp,0.0_dp)

  do i = 1, ntriangle*nipp
     kk1 = (i-1)*ntriangle*nipp+1
     kk2 = i*ntriangle*nipp
     cout(i) = dot_product(conjg(zzsparse(kk1:kk2)),cin(1:ntriangle*nipp))
  enddo

end subroutine convol
!*************************************************************************

max_stars_repo_name: wawando/stars-h

max_issues_repo_name: wawando/stars-h

max_forks_repo_name: wawando/stars-h


text:
      MODULE interp_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:21  03/09/06
      subroutine interp(np, nq, mode, e, fp, cp)
      USE vast_kind_param,ONLY: DOUBLE
      use molkst_C, only : norbs
      integer  :: np
      integer  :: nq
      integer , intent(inout) :: mode
      real(double) , intent(in) :: e
      real(double) , intent(in) :: fp((norbs*(norbs + 1))/2)
      real(double) , intent(inout) :: cp(norbs,norbs)
      END SUBROUTINE
      END INTERFACE
      END MODULE

max_stars_repo_name: openmopac/MOPAC-archive

max_issues_repo_name: openmopac/MOPAC-archive

max_forks_repo_name: openmopac/MOPAC-archive


text:
C
C***********************************************************************
C
      SUBROUTINE DLARAN(N, X)
C
C  THIS SUBROUTINE SETS THE VECTOR X TO RANDOM NUMBERS
C
C  FORMAL PARAMETERS
C
      INTEGER N
      DOUBLE PRECISION X(N)
C
C  LOCAL VARIABLES
C
      INTEGER I, IURAND
C
C  FUNCTIONS CALLED
C
      REAL URAND
      DOUBLE PRECISION DBLE
C
C  SUBROUTINES CALLED
C
C     NONE
C
C  INITIALIZE SEED
C
      DATA IURAND /0/
C
      DO 10 I = 1, N
         X(I) = DBLE(URAND(IURAND)) - 0.5D0
   10 CONTINUE
      RETURN
      END

max_stars_repo_name: eile/ITK

max_issues_repo_name: GEHC-Surgery/ITK

max_forks_repo_name: GEHC-Surgery/ITK


text:
C
C $Id: gasetr.f,v 1.8 2008-07-27 00:17:14 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 GASETR (PNAM,RVAL)
C
        CHARACTER*(*) PNAM
C
C The subroutine GASETR may be used to set GRIDAL parameters which have
C values of type REAL.
C
C Declare the common block containing real and integer parameters.
C
        COMMON /GAREIN/ ICAX,ICLB,ICMJ,ICMN,ILTY,IORX,NCFX,NCFY,RCWX,
     +                  RCWY,RDCX,RDCY,RMJX,RMJY,RMNX,RMNY,RWAX,RWLB,
     +                  RWMJ,RWMN
        SAVE   /GAREIN/
C
C Do a call forcing a BLOCKDATA to be loaded from a binary library.
C
        CALL GABLDT
C
C Check for an uncleared prior error.
C
        IF (ICFELL('GASETR - UNCLEARED PRIOR ERROR',1).NE.0) RETURN
C
C Set the selected parameter.
C
        IF      (PNAM(1:3).EQ.'CAX'.OR.PNAM(1:3).EQ.'cax') THEN
          ICAX=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'CLB'.OR.PNAM(1:3).EQ.'clb') THEN
          ICLB=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'CMJ'.OR.PNAM(1:3).EQ.'cmj') THEN
          ICMJ=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'CMN'.OR.PNAM(1:3).EQ.'cmn') THEN
          ICMN=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'LTY'.OR.PNAM(1:3).EQ.'lty') THEN
          ILTY=MAX(0,MIN(2,INT(RVAL)))
        ELSE IF (PNAM(1:3).EQ.'WAX'.OR.PNAM(1:3).EQ.'wax') THEN
          RWAX=RVAL
        ELSE IF (PNAM(1:3).EQ.'WLB'.OR.PNAM(1:3).EQ.'wlb') THEN
          RWLB=RVAL
        ELSE IF (PNAM(1:3).EQ.'WMJ'.OR.PNAM(1:3).EQ.'wmj') THEN
          RWMJ=RVAL
        ELSE IF (PNAM(1:3).EQ.'WMN'.OR.PNAM(1:3).EQ.'wmn') THEN
          RWMN=RVAL
        ELSE IF (PNAM(1:3).EQ.'XLL'.OR.PNAM(1:3).EQ.'xll') THEN
          NCFX=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'XLO'.OR.PNAM(1:3).EQ.'xlo') THEN
          RDCY=RVAL
        ELSE IF (PNAM(1:3).EQ.'XLS'.OR.PNAM(1:3).EQ.'xls') THEN
          RCWX=RVAL
        ELSE IF (PNAM(1:3).EQ.'XMJ'.OR.PNAM(1:3).EQ.'xmj') THEN
          RMJX=RVAL
        ELSE IF (PNAM(1:3).EQ.'XMN'.OR.PNAM(1:3).EQ.'xmn') THEN
          RMNX=RVAL
        ELSE IF (PNAM(1:3).EQ.'XOR'.OR.PNAM(1:3).EQ.'xor') THEN
          IORX=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'YLL'.OR.PNAM(1:3).EQ.'yll') THEN
          NCFY=INT(RVAL)
        ELSE IF (PNAM(1:3).EQ.'YLO'.OR.PNAM(1:3).EQ.'ylo') THEN
          RDCX=RVAL
        ELSE IF (PNAM(1:3).EQ.'YLS'.OR.PNAM(1:3).EQ.'yls') THEN
          RCWY=RVAL
        ELSE IF (PNAM(1:3).EQ.'YMJ'.OR.PNAM(1:3).EQ.'ymj') THEN
          RMJY=RVAL
        ELSE IF (PNAM(1:3).EQ.'YMN'.OR.PNAM(1:3).EQ.'ymn') THEN
          RMNY=RVAL
        ELSE
          CALL SETER ('GASETR - UNRECOGNIZED PARAMETER NAME',2,1)
        END IF
C
C Done.
C
        RETURN
C
      END

max_stars_repo_name: tenomoto/ncl

max_issues_repo_name: tenomoto/ncl

max_forks_repo_name: tenomoto/ncl


text:
\ path for include

: find-last-/  ( a n -- a n' )
    begin dup while
        1- 2dup + c@ '/' = if  1+ exit  then
    repeat ;

\ return the path up to and including the last /
\ returns ./ if no / in the path
: dirname ( a n -- a n' )
    find-last-/  dup 0= if  2drop s" ./"  then ;

\ returns the filename after the last / (which may be empty)
\ return original input if no /
: basename ( a n -- a' n' )
    2dup find-last-/  nip /string ;

\ true if filename starts with "./"
: ./?  ( a n -- a n f )
    over 2 s" ./" compare 0=  over 1 > and ;

\ append ./basename to dirname of the current source file
: append-to-dir ( a n -- a' n' )
    source-name @ count dirname
    >r pad r@ move ( a n )
    2 /string  pad r@ + swap  dup>r move
    pad r> r> + ;

\ If the file starts with "./", try to include it from
\ the same directory as the current source file.
: open-on-path ( a n -- a' n' fid )
    ./? source-file @ file? and if
        2dup append-to-dir
        2dup r/o open-file 0= if  >r rot drop rot drop r> exit  then  drop
    then
    \ todo: search path goes here
    2dup r/o open-file abort" file not found" ;

: included  ( str len -- )  open-on-path include-file ;

: include   parse-name included ;

\S
: included? ( name len -- f )  2>r includes
    begin @ dup while
        dup cell+ count  2r@ compare 0=
    until then  r>drop r>drop ;

\ : required ( a n -- )  2dup included? if 2drop else included then ;
\ : require  parse-name required ;

: .included ( -- ) \ first to last
    includes begin dup @ ?dup 0= until
    begin dup includes = not while cell+ count type cr repeat drop ;
\    includes begin @ ?dup while dup cell+ count type cr repeat ;

max_stars_repo_name: amckewan/mcforth

max_issues_repo_name: amckewan/mcforth

max_forks_repo_name: amckewan/mcforth


text:
program eval_entropy
!
! Purpose: To evaluate the entropy of a given probability
! distribution function. The binning should be linearly
! graduated.
!
! Record of revisions:
!  Date         Who                Description
! ==========   ===============    =====================================
! 9-12-2011     ThMosqueiro        Original code.
!
!

  implicit none

  character(len=100)     :: bugger, inputfile
  real(8)                :: Entropy, dx, logbase

! Volatile variables
  real(8)                :: x, prob

! Purely auxiliary variables
  real(8)                :: x0, x1, prob0, prob1
  integer                :: Reason

! Reading input file name into inputfile variable.
  call getarg(1, bugger)
  read(bugger,*) inputfile

! This makes the correct base for entropy
  logbase = 1.d0/dlog(2.d0)

! Openning as read-only the inputfile
  open(unit=27, file=inputfile, action="read", status="old")

! Initial condition for variables
  Entropy = 0.d0

  read(27,*) x0, prob0
  read(27,*) x1, prob1

  dx = dabs(x0 - x1)

  if ( prob0 .ne. 0.d0 ) then
     Entropy = - prob0*dx*dlog(prob0*dx)*logbase + Entropy
  end if
  if ( prob1 .ne. 0.d0 ) then
     Entropy = - prob1*dx*dlog(prob1*dx)*logbase + Entropy
  end if

! Now read every line of the file and evaluate its contribution to
! the entropy.
  do
     read(27, *, IOSTAT=Reason) x, prob

     if ( Reason .lt. 0 ) then
        close(27)
        exit
     else

        if ( prob .ne. 0.d0 ) then
           Entropy = - prob*dx*dlog(prob*dx)*logbase + Entropy
        end if

     end if

  end do

  close(27)

  write(*,*) Entropy

  stop
end program eval_entropy

max_stars_repo_name: thmosqueiro/Probab-Greenberg-Hastings

max_issues_repo_name: thmosqueiro/Probab-Greenberg-Hastings

max_forks_repo_name: thmosqueiro/Probab-Greenberg-Hastings


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!
! Test the fix for segment fault about procedure pointer when all the following
! three conditions are satisfied:
! 1) module procedue used as proc-interface
! 2) keyword RESULT appears
! 3) none dummy argument

module m
  procedure(g), pointer :: f
contains
  function g() result(z)
    integer :: z(3)
    z = (/1, 2, 3/)
  end function
end module

program test
  use m
  implicit none
  integer :: a(3)
  f => g
  a = f()
  if (.not. all(a .eq. (/1, 2, 3/))) STOP 1

  print *, 'PASS'
end program

max_stars_repo_name: vigbalu/flang

max_issues_repo_name: vigbalu/flang

max_forks_repo_name: vigbalu/flang


text:
      SUBROUTINE ZHST01( N, ILO, IHI, A, LDA, H, LDH, Q, LDQ, WORK,
     $                   LWORK, RWORK, RESULT )
*
*  -- LAPACK test routine (version 3.1) --
*     Univ. of Tennessee, Univ. of California Berkeley and NAG Ltd..
*     November 2006
*
*     .. Scalar Arguments ..
      INTEGER            IHI, ILO, LDA, LDH, LDQ, LWORK, N
*     ..
*     .. Array Arguments ..
      DOUBLE PRECISION   RESULT( 2 ), RWORK( * )
      COMPLEX*16         A( LDA, * ), H( LDH, * ), Q( LDQ, * ),
     $                   WORK( LWORK )
*     ..
*
*  Purpose
*  =======
*
*  ZHST01 tests the reduction of a general matrix A to upper Hessenberg
*  form:  A = Q*H*Q'.  Two test ratios are computed;
*
*  RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
*  RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
*
*  The matrix Q is assumed to be given explicitly as it would be
*  following ZGEHRD + ZUNGHR.
*
*  In this version, ILO and IHI are not used, but they could be used
*  to save some work if this is desired.
*
*  Arguments
*  =========
*
*  N       (input) INTEGER
*          The order of the matrix A.  N >= 0.
*
*  ILO     (input) INTEGER
*  IHI     (input) INTEGER
*          A is assumed to be upper triangular in rows and columns
*          1:ILO-1 and IHI+1:N, so Q differs from the identity only in
*          rows and columns ILO+1:IHI.
*
*  A       (input) COMPLEX*16 array, dimension (LDA,N)
*          The original n by n matrix A.
*
*  LDA     (input) INTEGER
*          The leading dimension of the array A.  LDA >= max(1,N).
*
*  H       (input) COMPLEX*16 array, dimension (LDH,N)
*          The upper Hessenberg matrix H from the reduction A = Q*H*Q'
*          as computed by ZGEHRD.  H is assumed to be zero below the
*          first subdiagonal.
*
*  LDH     (input) INTEGER
*          The leading dimension of the array H.  LDH >= max(1,N).
*
*  Q       (input) COMPLEX*16 array, dimension (LDQ,N)
*          The orthogonal matrix Q from the reduction A = Q*H*Q' as
*          computed by ZGEHRD + ZUNGHR.
*
*  LDQ     (input) INTEGER
*          The leading dimension of the array Q.  LDQ >= max(1,N).
*
*  WORK    (workspace) COMPLEX*16 array, dimension (LWORK)
*
*  LWORK   (input) INTEGER
*          The length of the array WORK.  LWORK >= 2*N*N.
*
*  RWORK   (workspace) DOUBLE PRECISION array, dimension (N)
*
*  RESULT  (output) DOUBLE PRECISION array, dimension (2)
*          RESULT(1) = norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
*          RESULT(2) = norm( I - Q'*Q ) / ( N * EPS )
*
*  =====================================================================
*
*     .. Parameters ..
      DOUBLE PRECISION   ONE, ZERO
      PARAMETER          ( ONE = 1.0D+0, ZERO = 0.0D+0 )
*     ..
*     .. Local Scalars ..
      INTEGER            LDWORK
      DOUBLE PRECISION   ANORM, EPS, OVFL, SMLNUM, UNFL, WNORM
*     ..
*     .. External Functions ..
      DOUBLE PRECISION   DLAMCH, ZLANGE
      EXTERNAL           DLAMCH, ZLANGE
*     ..
*     .. External Subroutines ..
      EXTERNAL           DLABAD, ZGEMM, ZLACPY, ZUNT01
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC          DCMPLX, MAX, MIN
*     ..
*     .. Executable Statements ..
*
*     Quick return if possible
*
      IF( N.LE.0 ) THEN
         RESULT( 1 ) = ZERO
         RESULT( 2 ) = ZERO
         RETURN
      END IF
*
      UNFL = DLAMCH( 'Safe minimum' )
      EPS = DLAMCH( 'Precision' )
      OVFL = ONE / UNFL
      CALL DLABAD( UNFL, OVFL )
      SMLNUM = UNFL*N / EPS
*
*     Test 1:  Compute norm( A - Q*H*Q' ) / ( norm(A) * N * EPS )
*
*     Copy A to WORK
*
      LDWORK = MAX( 1, N )
      CALL ZLACPY( ' ', N, N, A, LDA, WORK, LDWORK )
*
*     Compute Q*H
*
      CALL ZGEMM( 'No transpose', 'No transpose', N, N, N,
     $            DCMPLX( ONE ), Q, LDQ, H, LDH, DCMPLX( ZERO ),
     $            WORK( LDWORK*N+1 ), LDWORK )
*
*     Compute A - Q*H*Q'
*
      CALL ZGEMM( 'No transpose', 'Conjugate transpose', N, N, N,
     $            DCMPLX( -ONE ), WORK( LDWORK*N+1 ), LDWORK, Q, LDQ,
     $            DCMPLX( ONE ), WORK, LDWORK )
*
      ANORM = MAX( ZLANGE( '1', N, N, A, LDA, RWORK ), UNFL )
      WNORM = ZLANGE( '1', N, N, WORK, LDWORK, RWORK )
*
*     Note that RESULT(1) cannot overflow and is bounded by 1/(N*EPS)
*
      RESULT( 1 ) = MIN( WNORM, ANORM ) / MAX( SMLNUM, ANORM*EPS ) / N
*
*     Test 2:  Compute norm( I - Q'*Q ) / ( N * EPS )
*
      CALL ZUNT01( 'Columns', N, N, Q, LDQ, WORK, LWORK, RWORK,
     $             RESULT( 2 ) )
*
      RETURN
*
*     End of ZHST01
*
      END

max_stars_repo_name: nashmit/lapack-release

max_issues_repo_name: nashmit/lapack-release

max_forks_repo_name: nashmit/lapack-release


text:
      subroutine Bamp_pmm(q,mc,ms,Bpmm)
c--- u + g  ->  c + s + d  (t-channel single-charm)
************************************************************************
*                                                                      *
* AUTHORS: R. FREDERIX AND F. TRAMONTANO                               *
* DATE  : 12/17/2008                                                    *
*                                                                      *
************************************************************************
      implicit none
      include 'constants.f'
      include 'zprods_com.f'
      include 'epinv.f'
      include 'stopf1inc.f'
      double precision q(mxpart,4),dot,cDs,gDs,cDg,mc,ms,
     . mc2,ms2,qsq,s,t,u,xsn,xsd,xs
      double complex trc,trg,trs,trsgc,zpmm,Bpmm

      mc2=mc**2
      ms2=ms**2

      cDs=dot(q,3,4)+mc2*dot(q,4,2)/2d0/dot(q,3,2)
     . +ms2*dot(q,3,2)/2d0/dot(q,4,2)
      cDg=dot(q,3,2)
      gDs=dot(q,4,2)
      qsq=mc2+ms2+2d0*cDs+2d0*cDg+2d0*gDs
      s=ms2+2d0*gDs
      t=mc2+2d0*cDg
      u=mc2+ms2+2d0*cDs

      xsn=(1d0-dsqrt(1d0-4d0*ms*mc/(u-(ms-mc)**2)))
      xsd=(1d0+dsqrt(1d0-4d0*ms*mc/(u-(ms-mc)**2)))
      xs=-xsn/xsd

      trg=2d0*za(5,2)*zb(2,1)
      trs=2d0*za(5,4)*zb(4,1)+ms**2*za(5,2)*zb(2,1)/gDs
      trc=2d0*za(5,3)*zb(3,1)+mc**2*za(5,2)*zb(2,1)/cDg
      trsgc=2d0*zb(1,4)*za(4,2)*zb(2,3)*za(3,5)
      zpmm=za(2,4)*za(2,3)*zb(4,3)**2

      Bpmm = ms*(gDs+cDg)*(mc2*gDs**3*(cDs*(4*trc*gDs**2-2*(trsgc+mc2*t
     &   rg-2*ms2*trc)*gDs-ms2*(trsgc+3*mc2*trg))+mc2*ms2*(2*(trg-trs)*g
     &   Ds-trsgc+3*ms2*trg)-4*trg*cDs**2*(gDs+ms2))-cDg**2*gDs*(2*cDs**
     &   2*(4*trc*gDs**2-2*trsgc*gDs+2*ms2*(trs+trg+3*trc)*gDs-2*ms2*trs
     &   gc+ms2**2*trg)+ms2*cDs*(-4*trc*gDs**2+2*(trsgc+mc2*trg-2*ms2*tr
     &   c)*gDs+ms2*(trsgc+3*mc2*trg))+mc2*ms2**2*(-2*(trs+trg)*gDs+trsg
     &   c-ms2*trg))+2*cDg*gDs**2*(cDs*(mc2*((2*gDs+3*ms2)*(trc*gDs-ms2*
     &   trg)+2*ms2*trs*gDs)+mc2*trg*cDs*(2*gDs+3*ms2)-4*trc*cDs*gDs*(gD
     &   s+ms2)+4*trg*cDs**2*(gDs+ms2))+trsgc*(2*cDs**2*(gDs+ms2)-mc2*cD
     &   s*gDs-mc2*ms2**2))+2*ms2*cDg**3*cDs*(trc*gDs*(2*gDs+3*ms2)-trsg
     &   c*(gDs+ms2)))*tr5Xs/(cDg*(cDs-mc*ms)*(cDs+mc*ms)*gDs**3)/2.0d+0

      Bpmm = Bpmm-mc2*ms*(gDs+cDg)*(2*mc2**2*trg*cDs*gDs**3+cDg**3*(cD
     &   s*(-4*mc2*trs*gDs-2*ms2*trc*gDs+ms2*trsgc+3*mc2*ms2*trg)-4*trg*
     &   cDs**3-2*(trsgc+ms2*trg)*cDs**2+mc2*ms2*(trsgc+ms2*trg))+2*cDg*
     &   *2*gDs*(mc2**2*(trs*gDs-ms2*trg)+cDs**2*(2*trc*gDs-trsgc+2*mc2*
     &   trg)+2*trg*cDs**3+mc2*trsgc*cDs)-mc2*cDg*gDs**2*(cDs*(2*trc*gDs
     &   -trsgc+mc2*trg)+6*trg*cDs**2+mc2*(trsgc-ms2*trg))+cDg**4*(4*trs
     &   *cDs**2-2*mc2*ms2*trs))*tr5Xc/(cDg**3*(cDs**2-mc2*ms2)*gDs)/2.0
     &   d+0

      Bpmm = ms*(cDg**2*gDs*(2*cDs**2*(4*trc*gDs**2-2*trsgc*gDs+2*ms2*(
     &   trs+trg+3*trc)*gDs-2*ms2*trsgc+ms2**2*trg)+ms2*cDs*(-4*trc*gDs*
     &   *2+2*(trsgc+mc2*trg-2*ms2*trc)*gDs+ms2*(trsgc+3*mc2*trg))+mc2*m
     &   s2**2*(-2*(trs+trg)*gDs+trsgc-ms2*trg))+mc2*gDs**3*(cDs*(-4*trc
     &   *gDs**2+2*(trsgc+mc2*trg-2*ms2*trc)*gDs+ms2*(trsgc+3*mc2*trg))+
     &   mc2*ms2*(2*trs*gDs-2*trg*gDs+trsgc-3*ms2*trg)+4*trg*cDs**2*(gDs
     &   +ms2))+2*cDg*gDs**2*(trsgc*(-2*cDs**2*(gDs+ms2)+mc2*cDs*gDs+mc2
     &   *ms2**2)-cDs*(mc2*((2*gDs+3*ms2)*(trc*gDs-ms2*trg)+2*ms2*trs*gD
     &   s)+mc2*trg*cDs*(2*gDs+3*ms2)-4*trc*cDs*gDs*(gDs+ms2)+4*trg*cDs*
     &   *2*(gDs+ms2)))+2*ms2*cDg**3*cDs*(trsgc*(gDs+ms2)-trc*gDs*(2*gDs
     &   +3*ms2)))*tr4Xs/(cDg*(cDs-mc*ms)*(cDs+mc*ms)*gDs**2)/2.0d+0+Bp
     &   mm

      Bpmm = mc2*ms*(2*mc2**2*trg*cDs*gDs**3+cDg**3*(cDs*(-4*mc2*trs*gD
     &   s-2*ms2*trc*gDs+ms2*trsgc+3*mc2*ms2*trg)-4*trg*cDs**3-2*(trsgc+
     &   ms2*trg)*cDs**2+mc2*ms2*(trsgc+ms2*trg))+2*cDg**2*gDs*(mc2**2*(
     &   trs*gDs-ms2*trg)+cDs**2*(2*trc*gDs-trsgc+2*mc2*trg)+2*trg*cDs**
     &   3+mc2*trsgc*cDs)-mc2*cDg*gDs**2*(cDs*(2*trc*gDs-trsgc+mc2*trg)+
     &   6*trg*cDs**2+mc2*(trsgc-ms2*trg))+cDg**4*(4*trs*cDs**2-2*mc2*ms
     &   2*trs))*tr4Xc/(cDg**2*(cDs**2-mc2*ms2)*gDs)/2.0d+0+Bpmm

      Bpmm = Bpmm-ms*(-2*mc2**3*ms2*trg*gDs**3-2*cDg**2*gDs*(mc2*ms2*(
     &   -13*(trg+trc)*gDs**2+2*(trsgc+mc2*trg)*gDs+mc2*(trsgc+2*ms2*trg
     &   ))+mc2*cDs*(mc2*trs*gDs+2*ms2*trc*gDs-ms2*trsgc+mc2*ms2*trg)+cD
     &   s**2*gDs*(13*(trg+trc)*gDs-2*(trsgc+mc2*trg)))+cDg**3*(6*(2*trg
     &   -3*trs)*(cDs**2-mc**2*ms**2)*gDs**2+2*((trsgc+ms2*trg)*cDs**2+m
     &   c2*ms2*(-trsgc-ms2*trg+ms2*trc)+2*mc2**2*ms2*trs)*gDs+mc2*ms2*(
     &   trsgc*cDs+ms2*trg*cDs-ms2*trsgc+mc2*ms2*trg))-2*cDg**4*(2*mc2*m
     &   s2*(trs*gDs-ms2*trg)+2*cDs**2*(ms2*trg-trs*gDs)+mc2*ms2*trs*cDs
     &   )+mc2*cDg*gDs**2*(2*(mc2*ms2*(trg+trc)-trg*cDs**2)*gDs+mc2*(trs
     &   gc*cDs+5*ms2*trg*cDs-ms2*trsgc+mc2*ms2*trg)))*tr3Xs/(cDg**2*(cD
     &   s**2-mc**2*ms**2)*gDs)/2.0d+0

      Bpmm = ms*(2*cDg**2*gDs*(mc2*cDs*(8*trc*gDs**2-4*trsgc*gDs+12*ms2
     &   *trc*gDs-3*ms2*trsgc+ms2**2*trg)+mc2*ms2*(3*trc*gDs**2-(trsgc+1
     &   3*ms2*trc)*gDs+ms2*(trsgc+3*mc2*trg))+cDs**2*gDs*(-7*trc*gDs+3*
     &   trsgc+2*mc2*trg+9*ms2*trc)+2*(trs+trg)*cDs**3*gDs)-2*mc2*cDg*gD
     &   s**2*(mc2*(4*trc*gDs**2+6*ms2*trc*gDs+trsgc*(-2*gDs-ms2)+ms2**2
     &   *trg)+cDs*(-8*trc*gDs**2+4*(trsgc+mc2*trg-2*ms2*trc)*gDs+2*ms2*
     &   (trsgc+3*mc2*trg))+cDs**2*(4*(trs+trg)*gDs+trsgc+ms2*trg))+mc2*
     &   gDs**3*(mc2*(-8*trc*gDs**2+(4*trsgc+4*mc2*trg-11*ms2*trc)*gDs+m
     &   s2*(2*trsgc+6*mc2*trg+3*ms2*trc))+2*mc2*cDs*(2*(trs+trg)*gDs+tr
     &   sgc+ms2*trg)+3*trc*cDs**2*(gDs-ms2))+4*cDg**3*(mc2*ms2*(gDs*(2*
     &   (trs+trg)*gDs+ms2*(trg-3*trc))+trsgc*(2*gDs+ms2))-cDs**2*gDs*(2
     &   *(trs+trg+trc)*gDs+trsgc+ms2*trg))-8*trg*cDg**4*(cDs-mc*ms)*(cD
     &   s+mc*ms)*gDs)*tr3Xc/(cDg*(cDs-mc*ms)*(cDs+mc*ms)*gDs**2)/4.0d+0
     &   +Bpmm

      Bpmm = ms*(-mc2*cDg**3*gDs*(mc2*ms2*((36*trs-24*trg+5*trc)*gDs**2
     &   +(ms2*(2*trs-2*trg+trc)+8*mc2*trs)*gDs+2*ms2**2*(trs+trg))-cDs*
     &   *2*(3*(12*trs-8*trg+7*trc)*gDs**2-8*trsgc*gDs+ms2*(2*trs-2*trg+
     &   9*trc)*gDs+2*ms2**2*(trs+trg))-2*mc2*ms2*cDs*(4*(-trs-trg+trc)*
     &   gDs-3*trsgc+ms2*(trs+3*trg))+2*(ms2*(trs+3*trg)-3*trsgc)*cDs**3
     &   )+2*cDg**4*(mc2**2*ms2*(8*trs*gDs**2-(trsgc+ms2*(trs+7*trg+2*tr
     &   c))*gDs+ms2*(trsgc+ms2*trs))+mc2*cDs**2*(ms2*(trs*gDs+7*trg*gDs
     &   -trsgc)+gDs*(trsgc-8*trs*gDs)-ms2**2*(trs+3*trc))+cDs**4*(2*(tr
     &   s+trg)*gDs+3*ms2*trc)+mc2*ms2*cDs*gDs*(3*trc*gDs-trsgc+2*mc2*tr
     &   s-2*mc2*trg-11*ms2*trc)+cDs**3*gDs*(-7*trc*gDs+3*trsgc+2*mc2*tr
     &   g+9*ms2*trc))+mc2*cDg**2*gDs**2*(cDs**2*(52*(trg+trc)*gDs**2+(-
     &   8*trsgc+ms2*(2*trs-trg+trc)+4*mc2*trs-8*mc2*trg)*gDs-2*(3*mc2*t
     &   rsgc+ms2**2*(trs+trg)-mc2*ms2*trs))+mc2*ms2*(-52*(trg+trc)*gDs*
     &   *2+(8*trsgc+ms2*(-2*trs+trg-trc)-4*mc2*(trc-3*trg))*gDs+2*(3*mc
     &   2*trsgc+ms2**2*(trs+trg)-mc2*ms2*trs))-mc2*cDs*(8*trc*gDs**2+(5
     &   *ms2*trc-4*(trsgc+mc2*trs))*gDs+ms2**2*(2*(trs+trg)+3*trc))+cDs
     &   **3*(trc*gDs+ms2*(2*(trs+trg)+3*trc)))+2*mc2*(cDs-mc*ms)*(cDs+m
     &   c*ms)*gDs**3*(trs*(cDs+ms2)*(4*gDs**2+ms2*(4*(gDs+cDs)-mc2)+8*c
     &   Ds*gDs-2*mc2*gDs+4*cDs**2-2*mc2*cDs+ms2**2+mc2**2)+ms2*(trg+trc
     &   )*gDs*(2*gDs+2*cDs+ms2-mc2))-2*cDg**5*(cDs-mc*ms)*(cDs+mc*ms)*(
     &   cDs*(2*(trs+trg+2*trc)*gDs+2*trsgc+2*ms2*trg-3*ms2*trc)+gDs*(7*
     &   trc*gDs-3*trsgc-2*mc2*trg-9*ms2*trc))+2*mc2*cDg*(cDs-mc*ms)*(cD
     &   s+mc*ms)*gDs**3*((ms2*(4*trs+trg+trc)+2*mc2*trg)*gDs+2*(2*ms2+m
     &   c2)*trs*cDs+3*ms2*(trg+trc)*cDs+(2*ms2**2+mc2**2)*trs)-4*cDg**6
     &   *(cDs-mc*ms)*(cDs+mc*ms)*(2*(trs+trg+trc)*gDs+2*trg*cDs+trsgc+m
     &   s2*trg)-8*trg*cDg**7*(cDs-mc*ms)*(cDs+mc*ms))*lVs/(mc2*cDg**3*(
     &   mc2*ms2-cDs**2)*gDs**2)/8.0d+0+Bpmm

      Bpmm = (-2*cDg**4*(cDs**2*(4*trs*gDs**3+ms2*(6*trc-5*trg)*gDs**2-
     &   ms2*(3*trsgc+2*mc2*(trs+trg)+ms2*(2*trg+13*trc))*gDs+ms2**2*(2*
     &   trsgc-mc2*(trs+trg+3*trc)))+mc2*ms2*(-4*trs*gDs**3+ms2*(5*trg-2
     &   *trc)*gDs**2+ms2*(trsgc+2*ms2*trg+2*mc2*trg+15*ms2*trc)*gDs+ms2
     &   **2*(mc2*(trs+trg)-2*trsgc))+cDs**3*gDs*(4*trs*gDs+ms2*(3*trc-2
     &   *(trs+5*trg)))-mc2*ms2*cDs*gDs*(4*trs*gDs+ms2*(trc-8*trg))+3*ms
     &   2*trc*cDs**4)-2*cDg**3*gDs*(cDs**2*(6*(2*trg-3*trs)*gDs**3+2*(t
     &   rsgc-4*mc2*trs+ms2*trg)*gDs**2+mc2*ms2*(5*trs+6*trg-4*trc)*gDs+
     &   mc2*ms2**2*(trs+trg))-mc2*ms2*(6*(2*trg-3*trs)*gDs**3+2*(trsgc-
     &   4*mc2*trs+ms2*trg)*gDs**2+mc2*ms2*(trs+2*trg)*gDs+mc2*ms2**2*(t
     &   rs+trg))+cDs**3*(6*(2*trg-3*trs)*gDs**2+(-2*trsgc+6*ms2*trg+8*m
     &   s2*trc)*gDs-ms2*(4*trsgc+mc2*trs+2*mc2*trg))+mc2*ms2*cDs*(-2*(-
     &   9*trs+6*trg+4*trc)*gDs**2+(6*trsgc+4*mc2*trs-6*ms2*(trg+2*trc))
     &   *gDs+ms2*(4*trsgc+mc2*(trs+2*trg))))+cDg**2*gDs**2*(cDs**2*(52*
     &   (trg+trc)*gDs**3-8*(trsgc+mc2*trg)*gDs**2+mc2*(-8*trsgc+8*mc2*t
     &   rs+ms2*(-2*trs+trg+10*trc))*gDs-mc2*ms2*(4*trsgc+ms2*(3*trc-2*(
     &   trs+trg))+2*mc2*(trs+trg)))+mc2*ms2*(-52*(trg+trc)*gDs**3+8*(tr
     &   sgc+mc2*(trg-trc))*gDs**2-mc2*(4*(mc2*trs-3*trsgc)+ms2*(-2*trs+
     &   trg+14*trc))*gDs+mc2*ms2*(4*trsgc+ms2*(3*trc-2*(trs+trg))+2*mc2
     &   *(trs+trg)))+cDs**3*(52*(trg+trc)*gDs**2-8*(trsgc+mc2*trg)*gDs-
     &   2*mc2*ms2*(trs+trg+3*trc))+2*mc2*ms2*cDs*(-26*(trg+trc)*gDs**2+
     &   (4*trsgc+2*mc2*(trs+3*trg-trc))*gDs+mc2*ms2*(trs+trg+3*trc)))-2
     &   *mc2*(cDs-mc*ms)*(cDs+mc*ms)*gDs**3*(trs*(cDs+ms2)*(4*gDs**2+ms
     &   2*(4*(gDs+cDs)-mc2)+8*cDs*gDs-2*mc2*gDs+4*cDs**2-2*mc2*cDs+ms2*
     &   *2+mc2**2)+ms2*(trg+trc)*gDs*(2*gDs+2*cDs+ms2-mc2))-2*mc2*cDg*(
     &   cDs-mc*ms)*(cDs+mc*ms)*gDs**3*(-2*trg*gDs**2+cDs*(-2*trg*gDs+ms
     &   2*(4*trs+3*(trg+trc))+2*mc2*trs)+ms2*(4*trs+trg+trc)*gDs+mc2*tr
     &   sgc+2*ms2**2*trs+mc2*ms2*trg)+2*ms2*cDg**5*(mc2*ms2-cDs**2)*(2*
     &   (3*trs*gDs+2*(trg+trc)*gDs+trsgc+ms2*trg)+3*trc*cDs)+8*ms2*trg*
     &   cDg**6*(mc2*ms2-cDs**2))*lVc/(cDg**3*(mc2*ms*ms2-ms*cDs**2)*gDs
     &   **2)/8.0d+0+Bpmm

      Bpmm = xs*cDs*(2*trc*cDg*gDs-mc2*trg*gDs-trsgc*cDg)*(mc2*gDs**2-2
     &   *cDg*cDs*gDs+ms2*cDg**2)*lRcs/(mc*(xs**2-1)*cDg**2*gDs**2)+Bpm
     &   m

      Bpmm = (mc2*cDg**4*(-288*trs*gDs**6-2*cDs*(320*trs*gDs**5-4*ms2*(
     &   -68*trs+88*trg+199*trc)*gDs**4-2*ms2*(-180*trsgc+mc2*(2*trs-132
     &   *trg+51*trc)-34*ms2*trs+24*ms2*trg+298*ms2*trc)*gDs**3+ms2**2*(
     &   252*trsgc+6*ms2*trs+mc2*(-66*trs+104*trg-123*trc)+64*ms2*trg-90
     &   *ms2*trc)*gDs**2-2*ms2**2*(mc2*(9*trsgc+2*ms2*(8*(trs+trg)+9*tr
     &   c))-18*ms2*trsgc)*gDs-9*mc2*ms2**3*(trsgc+ms2*trg))+16*ms2*(-17
     &   5*trs+132*trg+6*trc)*gDs**5+2*ms2*(144*trsgc-652*ms2*trs-4*mc2*
     &   trs+756*ms2*trg+2*mc2*trg+186*ms2*trc+153*mc2*trc)*gDs**4-ms2*(
     &   mc2*(216*trsgc+104*ms2*trs+854*ms2*trg+785*ms2*trc)-12*ms2**2*(
     &   3*trs+21*trg+14*trc)+144*mc2**2*trg)*gDs**3-6*cDs**2*gDs*(2*gDs
     &   +ms2)*(24*trs*gDs**2+ms2*(24*(trs+trg)+13*trc)*gDs+ms2*(-3*trsg
     &   c+6*ms2*trg+3*mc2*(trc-2*trg)-22*ms2*trc))-ms2**2*(mc2*(144*trs
     &   gc+ms2*(50*trs+500*trg+937*trc))-12*ms2*(ms2*(trs+trg)-6*trsgc)
     &   +4*mc2**2*(5*trs+23*trg))*gDs**2+36*ms2*(trs+trg)*cDs**3*gDs*(2
     &   *gDs+ms2)-2*mc2*ms2**3*(-9*trsgc+3*ms2*(trs+7*trg+39*trc)+5*mc2
     &   *trs-13*mc2*trg)*gDs+18*mc2*ms2**4*(trsgc+mc2*trg))-2*mc2*cDg**
     &   3*gDs*(216*(2*trg-3*trs)*gDs**6+cDs*(432*(2*trg-3*trs)*gDs**5+8
     &   *(ms2*(-81*trs+104*trg+50*trc)-28*mc2*trs)*gDs**4-2*ms2*(72*trs
     &   gc+mc2*(10*trs+274*trg+229*trc)+6*ms2*trs-94*ms2*trg-85*ms2*trc
     &   )*gDs**3-ms2*(6*ms2*(12*trsgc+ms2*(trs+trg+3*trc))+mc2*(ms2*(-5
     &   2*trs+210*trg+319*trc)-180*trsgc)+mc2**2*(26*trs-92*trg))*gDs**
     &   2+mc2*ms2**2*(126*trsgc+ms2*(3*trs+32*trg-45*trc)-13*mc2*trs+82
     &   *mc2*trg)*gDs+18*mc2*ms2**3*(trsgc+mc2*trg))+4*(18*(trsgc-3*mc2
     &   *trs)+ms2*(-81*trs+536*trg+464*trc))*gDs**5+2*ms2*(-126*trsgc+m
     &   c2*(656*trs-611*trg+217*trc)-6*ms2*trs+485*ms2*trg+461*ms2*trc)
     &   *gDs**4+cDs**2*(2*gDs+ms2)*(108*(2*trg-3*trs)*gDs**3-3*(12*(trs
     &   gc-mc2*trs)+ms2*(12*trg+trc))*gDs**2+ms2*(19*ms2*trc-12*mc2*(tr
     &   c-3*(trs+trg)))*gDs+9*mc2*ms2*(trsgc+ms2*trg))+ms2*(-3*ms2*(48*
     &   trsgc+ms2*(-2*trs-5*trg+trc))+mc2*(ms2*(680*trs-1001*trg+154*tr
     &   c)-396*trsgc)+6*mc2**2*(29*trs-16*trg+12*trc))*gDs**3+ms2*(-3*m
     &   c2*ms2*(54*trsgc+7*ms2*trs+67*ms2*trg+29*ms2*trc)+mc2**2*(18*tr
     &   sgc+109*ms2*trs+100*ms2*trg+72*ms2*trc)+6*ms2**3*(trs+trg))*gDs
     &   **2+8*ms2*trc*cDs**3*gDs*(2*gDs+ms2)+mc2*ms2**2*(mc2*(27*trsgc+
     &   ms2*(11*trs+92*trg+18*trc))-3*ms2*(ms2*(trs+trg+9*trc)-6*trsgc)
     &   )*gDs+9*mc2**2*ms2**3*(trsgc+ms2*trg))+mc2*cDg**2*gDs**2*((2*gD
     &   s+ms2)*(936*(trg+trc)*gDs**5-36*(4*trsgc-9*mc2*trs+10*mc2*trg)*
     &   gDs**4+4*(-mc2*(45*trsgc+12*ms2*trs+506*ms2*trg+425*ms2*trc)+3*
     &   ms2**2*(4*trs+trg+trc)+72*mc2**2*trs)*gDs**3+ms2*(3*mc2*(72*trs
     &   gc+ms2*(-6*trs-trg+trc))+24*ms2**2*trs+mc2**2*(-462*trs+419*trg
     &   -229*trc))*gDs**2+mc2*ms2*(3*mc2*(78*trsgc+4*ms2*trs+70*ms2*trg
     &   -3*ms2*trc)-6*ms2**2*(trs+trg)+mc2**2*(46*trg-62*trs))*gDs+27*m
     &   c2**2*ms2**3*trc)+2*cDs*gDs*(1872*(trg+trc)*gDs**4+72*(-4*trsgc
     &   +9*mc2*trs+13*ms2*(trg+trc)-10*mc2*trg)*gDs**3+4*(-mc2*(36*trsg
     &   c+ms2*(-75*trs+245*trg+119*trc))-36*ms2*trsgc+3*ms2**2*(4*trs+3
     &   *(trg+trc))+108*mc2**2*trs)*gDs**2+ms2*(mc2*(72*trsgc-ms2*(6*tr
     &   s+304*trg+223*trc))+6*mc2**2*(33*trs+43*trg+5*trc)+6*ms2**2*(4*
     &   trs+3*(trg+trc)))*gDs+3*mc2*ms2**2*(24*trsgc+ms2*(trs+trg+3*trc
     &   )+mc2*(-3*trs+43*trg+5*trc)))+cDs**2*(2*gDs+ms2)*(936*(trg+trc)
     &   *gDs**3-36*(4*trsgc-9*mc2*trs+10*mc2*trg)*gDs**2+3*mc2*(12*trsg
     &   c+ms2*(12*trg+trc))*gDs-19*mc2*ms2**2*trc)-8*mc2*ms2*trc*cDs**3
     &   *(2*gDs+ms2))+2*ms2*cDg**5*(mc2*(6*(32*trs+101*trc)*gDs**4+(-14
     &   4*trsgc+8*mc2*(9*trs+19*trc)+84*ms2*trs-264*ms2*trg-491*ms2*trc
     &   )*gDs**3+(2*mc2*(45*trsgc+62*ms2*trs+107*ms2*trg+46*ms2*trc)-ms
     &   2*(36*trsgc+6*ms2*trs+132*ms2*trg+865*ms2*trc))*gDs**2+ms2*(mc2
     &   *(99*trsgc+ms2*(28*trs+181*trg-12*trc))-6*ms2*(ms2*(trs+trg+39*
     &   trc)-9*trsgc))*gDs+ms2**2*(mc2*(27*trsgc-8*ms2*trs+37*ms2*trg-1
     &   0*ms2*trc)+18*ms2*trsgc))+cDs*(-4*(4*mc2*(6*trs+10*trg+3*trc)-9
     &   *ms2*trc)*gDs**3+2*mc2*(54*trsgc+40*ms2*trs+60*ms2*trg-9*mc2*(t
     &   rc-4*trg)+207*ms2*trc)*gDs**2+mc2*ms2*(90*trsgc+64*ms2*trs+136*
     &   ms2*trg-9*mc2*(trc-4*trg)+219*ms2*trc)*gDs+18*mc2*ms2**2*(trsgc
     &   +ms2*trg))-6*cDs**2*(2*gDs+ms2)*(21*trc*gDs**2+(-9*trsgc-6*mc2*
     &   trg-22*ms2*trc+9*mc2*trc)*gDs+3*mc2*(trsgc+ms2*trg))+36*(trs+tr
     &   g)*cDs**3*gDs*(2*gDs+ms2))+2*mc2*cDg*gDs**3*(2*gDs+ms2)*(18*mc2
     &   *(15*trg+13*trc)*gDs**4+cDs*(12*(2*ms2*trs+mc2*(-2*trs+45*trg+3
     &   9*trc))*gDs**3-12*(mc2*(6*trsgc+ms2*(7*trs+trg+trc))-ms2**2*(6*
     &   trs+trg+trc)+mc2**2*(6*trg-trs))*gDs**2+(-3*mc2**2*(18*trsgc-6*
     &   ms2*trs+41*ms2*trg+23*ms2*trc)+30*ms2**3*trs+70*mc2**3*trs+9*mc
     &   2*ms2**2*(-4*trs+trg+trc))*gDs+18*mc2**2*ms2*(trsgc+mc2*trg))-1
     &   2*(mc2*(3*trsgc+ms2*(2*trs+trg+trc))-ms2**2*(2*trs+trg+trc)+3*m
     &   c2**2*trg)*gDs**3+3*(-3*mc2**2*(6*trsgc+29*ms2*trg+15*ms2*trc)+
     &   2*ms2**3*(4*trs+trg+trc)+12*mc2**3*trs+mc2*ms2**2*(-8*trs-3*(tr
     &   g+trc)))*gDs**2+6*cDs**2*gDs*((8*ms2*trs+mc2*(-8*trs+45*trg+39*
     &   trc))*gDs-2*mc2*(3*trsgc+5*ms2*trs)+mc2**2*(2*trs-6*trg)+8*ms2*
     &   *2*trs)-24*(mc-ms)*(ms+mc)*trs*cDs**3*gDs+ms2*(6*ms2**3-6*mc2*m
     &   s2**2+6*mc2**2*ms2-43*mc2**3)*trs*gDs+18*mc2**3*ms2*(trsgc+ms2*
     &   trg))+2*mc2**2*gDs**4*(2*gDs+ms2)*(18*mc2*trg*gDs**3+12*(ms2*tr
     &   s-mc2*trs+3*mc2*trg)*cDs*gDs**2+6*ms2*(ms2-mc2)*(2*trs+trg+trc)
     &   *gDs**2+6*cDs**2*(3*mc2*trg*gDs+(mc-ms)*(ms+mc)*trs*(-4*gDs-4*m
     &   s2+mc2))-9*mc2**2*trsgc*(gDs+cDs+ms2)+6*(mc-ms)*(ms+mc)*(mc2*tr
     &   s-ms2*(6*trs+trg+trc))*cDs*gDs+mc2**2*ms2*(6*trs+2*trg+29*trc)*
     &   gDs+3*ms2**3*(4*trs+trg+trc)*gDs-6*mc2*ms2**2*(3*trs+trg+trc)*g
     &   Ds-12*(mc-ms)*(ms+mc)*trs*cDs**3+15*ms2**3*trs*cDs-24*mc2*ms2**
     &   2*trs*cDs+5*mc2**3*trs*cDs-3*mc2**2*ms2*(3*trg-4*trs)*cDs-mc2**
     &   3*ms2*(11*trs+9*trg)+3*ms2**4*trs-6*mc2*ms2**3*trs+6*mc2**2*ms2
     &   **2*trs)-2*ms2*cDg**6*(2*gDs+ms2)*(-mc2*((72*(trs+2*trg)+89*trc
     &   )*gDs**2+(117*trsgc+ms2*(108*trs+234*trg+97*trc)+90*mc2*trg)*gD
     &   s+2*ms2*(27*trsgc-8*ms2*trs+19*ms2*trg+18*mc2*trg-10*ms2*trc))+
     &   6*cDs*(42*trc*gDs**2+(-18*trsgc+3*mc2*(3*trs-trg+5*trc)-49*ms2*
     &   trc)*gDs+6*mc2*(trsgc+ms2*trg))+36*cDs**2*(2*trc*gDs+trsgc+(ms2
     &   +mc2)*trg))-36*ms2*cDg**7*(2*gDs+ms2)*(7*trc*gDs**2+cDs*((6*(tr
     &   s+trg)+8*trc)*gDs+4*(trsgc+(ms2+mc2)*trg))+(2*mc2*(trs-4*trg+tr
     &   c)-3*(trsgc+3*ms2*trc))*gDs+4*trg*cDs**2+mc2*(trsgc-3*ms2*trg))
     &   -72*ms2*cDg**8*(2*gDs+ms2)*(2*(trs+trg+trc)*gDs+trg*(4*cDs+ms2+
     &   mc2)+trsgc)-144*ms2*trg*cDg**9*(2*gDs+ms2))/(mc2*ms*cDg**3*(2*c
     &   Dg+mc2)*gDs**3*(2*gDs+ms2))/2.4d+1+Bpmm

      Bpmm = B0cgsf*(cDg**4*(-2*mc2*(4*trs*gDs**4+ms2*(3*trs-5*trg-11*t
     &   rc)*gDs**3+ms2*(2*mc2*(trs+trg)+ms2*(trs+5*trc))*gDs**2+ms2**2*
     &   (2*trsgc+ms2*(trs+13*trc))*gDs-ms2**3*(trsgc+mc2*trg))+cDs*(-2*
     &   (4*mc2*trs+7*ms2*trc)*gDs**3+ms2*(6*(trsgc+3*ms2*trc)+mc2*(-16*
     &   trs+16*trg+7*trc))*gDs**2+4*mc2*ms2**2*(2*trs+3*trg+2*trc)*gDs+
     &   2*mc2*ms2**2*(trsgc+ms2*trg))+2*ms2*cDs**2*gDs*((2*(trs+trg)-7*
     &   trc)*gDs+3*trsgc+2*mc2*trg+9*ms2*trc-3*mc2*trc)+4*ms2*(trs+trg)
     &   *cDs**3*gDs)+mc2*cDg**2*gDs**2*(4*(9*trs+7*trg+13*trc)*gDs**4-(
     &   12*trsgc-16*mc2*trs+55*ms2*trg+8*mc2*trg+51*ms2*trc)*gDs**3+cDs
     &   *gDs*(4*(18*trs+trg+13*trc)*gDs**2-(8*(trsgc+mc2*(trg-3*trs))+m
     &   s2*(23*trg+22*trc))*gDs+ms2*(8*trsgc+ms2*(2*(trs+trg)+3*trc)-4*
     &   mc2*trs+12*mc2*trg))+(-mc2*(8*trsgc+ms2*(40*trs-24*trg+7*trc))+
     &   ms2*(8*trsgc+ms2*(2*trs-trg+trc))+8*mc2**2*trs)*gDs**2+cDs**2*(
     &   4*gDs*(9*trs*gDs-6*trg*gDs+trsgc)+4*ms2*trg*gDs-2*ms2**2*trc)-2
     &   *ms2*(2*mc2*(ms2*(-trs-3*trg+trc)-2*trsgc)+ms2**2*(trs+trg)+4*m
     &   c2**2*trs)*gDs-ms2*trc*cDs**3+3*mc2*ms2**3*trc)+mc2*cDg**3*gDs*
     &   (-4*(6*trg-7*trs)*gDs**4-2*cDs*(2*(6*trg-5*trs)*gDs**3-(2*trsgc
     &   -5*ms2*trs+11*ms2*trg+6*ms2*trc)*gDs**2-ms2*(trsgc-ms2*trs+2*mc
     &   2*trs-2*ms2*trg+5*ms2*trc)*gDs+2*ms2**2*(trsgc+mc2*trg))+(-4*tr
     &   sgc+16*mc2*trs+ms2*(-38*trs+31*trg+8*trc))*gDs**3-cDs**2*(8*trs
     &   *gDs**2+ms2*(8*trs-8*trg-3*trc)*gDs+2*ms2*(trsgc+ms2*trg))+ms2*
     &   (8*trsgc+2*ms2*trs-4*mc2*(trg+2*trc)+10*ms2*trg+9*ms2*trc)*gDs*
     &   *2+2*ms2*(mc2*(trsgc-ms2*(2*trs+9*trg+2*trc))+ms2**2*(trs+trg))
     &   *gDs-2*mc2*ms2**2*(trsgc+ms2*trg))+2*mc2*cDg*gDs**3*(26*(trg+tr
     &   c)*gDs**4+cDs*(52*(trg+trc)*gDs**3-2*(4*trsgc+3*mc2*trg)*gDs**2
     &   +(-2*mc2*(2*trsgc+ms2*(trs+7*trg+5*trc))+ms2**2*(4*trs+3*(trg+t
     &   rc))+6*mc2**2*trs)*gDs+2*mc2*ms2*(trsgc+mc2*trg))-2*(2*trsgc+mc
     &   2*trg)*gDs**3+(-2*mc2*(2*trsgc+ms2*(2*trs+14*trg+11*trc))+ms2**
     &   2*(4*trs+trg+trc)+4*mc2**2*trs)*gDs**2+2*cDs**2*gDs*(13*(trg+tr
     &   c)*gDs-2*(trsgc+mc2*trg))+(2*mc2*ms2*(trsgc-ms2*trs)+mc2**2*(ms
     &   2*(trg-2*trs)-trsgc)+2*ms2**3*trs)*gDs+2*mc2**2*ms2*(trsgc+ms2*
     &   trg))+2*mc2*gDs**4*(2*mc2*trg*gDs**3+4*(ms2*trs+mc2*(trg-trs))*
     &   cDs*gDs**2+2*ms2*(ms2-mc2)*(2*trs+trg+trc)*gDs**2+2*cDs**2*(mc2
     &   *trg*gDs+(mc-ms)*(ms+mc)*trs*(-4*gDs-4*ms2+mc2))-mc2**2*trsgc*(
     &   gDs+cDs+ms2)+2*(mc-ms)*(ms+mc)*(mc2*trs-ms2*(6*trs+trg+trc))*cD
     &   s*gDs+ms2**3*(4*trs+trg+trc)*gDs-2*mc2*ms2**2*(3*trs+trg+trc)*g
     &   Ds+mc2**2*ms2*(2*trs+trg+4*trc)*gDs+4*(ms2-mc2)*trs*cDs**3+5*ms
     &   2**3*trs*cDs-8*mc2*ms2**2*trs*cDs-mc2**2*ms2*(trg-4*trs)*cDs-mc
     &   2**3*ms2*(2*trs+trg)+ms2**4*trs-2*mc2*ms2**3*trs+2*mc2**2*ms2**
     &   2*trs)+ms2*cDg**5*(-14*trc*gDs**3+(6*(trsgc+3*ms2*trc)+mc2*(-4*
     &   trs+12*trg+9*trc))*gDs**2-2*cDs*gDs*(2*(trs+trg+9*trc)*gDs-4*tr
     &   sgc+2*ms2*trg-4*mc2*trg-18*ms2*trc+3*mc2*trc)-4*cDs**2*(2*trc*g
     &   Ds+trsgc+ms2*trg)+2*mc2*(3*trsgc+ms2*(6*trs+11*trg+trc))*gDs-2*
     &   mc2*ms2*(ms2*(trs-2*trg+trc)-3*trsgc))-2*ms2*cDg**6*((4*(trs+tr
     &   g)+11*trc)*gDs**2+2*cDs*((3*trs+5*trg+4*trc)*gDs+2*(trsgc+ms2*t
     &   rg))-(trsgc-2*ms2*trg+6*mc2*trg+9*ms2*trc)*gDs+4*trg*cDs**2-4*m
     &   c2*ms2*trg)-4*ms2*cDg**7*(2*(trs+2*trg+trc)*gDs+4*trg*cDs+trsgc
     &   +ms2*trg)-8*ms2*trg*cDg**8)/(mc2*ms*cDg**3*gDs**3)/8.0d+0+Bpmm

      Bpmm = tr3s00ft*(cDg**2*(12*(2*trg-3*trs)*gDs**4+cDs*(24*(2*trg-3
     &   *trs)*gDs**3+37*ms2*(trg+trc)*gDs**2-2*ms2*(4*trsgc-9*mc2*trs+8
     &   *mc2*trg)*gDs+2*mc2*ms2*(trsgc+ms2*trg))+(4*trsgc-16*mc2*trs+10
     &   5*ms2*trg+101*ms2*trc)*gDs**3+2*ms2*(mc2*(27*trs-19*trg+trc)-8*
     &   trsgc)*gDs**2+4*cDs**2*gDs*(-9*trs*gDs+6*trg*gDs-trsgc-ms2*trg)
     &   +2*mc2*ms2*(-7*trsgc+4*mc2*trs-7*ms2*trg-2*ms2*trc)*gDs+2*mc2*m
     &   s2**2*(trsgc+mc2*trg))-2*cDg*gDs*(26*(trg+trc)*gDs**4+cDs*(52*(
     &   trg+trc)*gDs**3-8*(trsgc+mc2*trg)*gDs**2+mc2*(-4*trsgc+2*mc2*tr
     &   s-17*ms2*trg-13*ms2*trc)*gDs+2*mc2*ms2*(trsgc+mc2*trg))-4*(trsg
     &   c+mc2*trg)*gDs**3+mc2*(-4*trsgc+4*mc2*trs-27*ms2*trg-19*ms2*trc
     &   )*gDs**2+2*cDs**2*gDs*(13*(trg+trc)*gDs-2*(trsgc+mc2*trg))+mc2*
     &   ms2*(2*trsgc-9*mc2*trs+mc2*trg)*gDs+2*mc2**2*ms2*(trsgc+ms2*trg
     &   ))+cDg**3*(8*trs*gDs**3+2*cDs*(16*trs*gDs**2+ms2*(15*trs-4*trg+
     &   4*trc)*gDs+ms2*(4*trsgc-2*mc2*trs+5*ms2*trg+3*ms2*trc))+ms2*(78
     &   *trs-59*trg-3*trc)*gDs**2+8*trs*cDs**2*gDs-2*ms2*(4*trsgc-7*mc2
     &   *trs+6*ms2*trg+4*ms2*trc)*gDs+2*ms2**2*(2*trsgc+5*mc2*trg))+2*m
     &   c2*gDs**2*(-2*trg*gDs**3-4*trg*cDs*gDs**2+(mc2*(trsgc-2*ms2*trg
     &   -5*ms2*trc)-2*trg*cDs**2)*gDs+mc2*((trsgc+ms2*trg)*cDs+ms2*(trs
     &   gc+3*mc2*trs+mc2*trg)))+2*ms2*cDg**4*(-5*trs*gDs-4*trs*cDs+7*ms
     &   2*trg+ms2*trc))/(ms*cDg**3*gDs)/2.0d+0+Bpmm

      Bpmm = lc*(cDg**2*(52*(trg+trc)*gDs**4+2*cDs*(26*(trg+trc)*gDs**3
     &   +(-4*trsgc+9*mc2*trs-10*mc2*trg)*gDs**2+mc2*(trsgc-7*ms2*trg)*g
     &   Ds+2*mc2**2*ms2*trg)-2*(4*trsgc-9*mc2*trs+10*mc2*trg)*gDs**3+2*
     &   mc2*(-5*trsgc+8*mc2*trs-ms2*trg)*gDs**2-6*mc2**2*ms2*trg*gDs+mc
     &   2**2*ms2*(trsgc+ms2*trg))-2*cDg**3*(6*(2*trg-3*trs)*gDs**3+2*(t
     &   rg*(6*cDs+ms2)-3*trs*(3*cDs+mc2)+trsgc)*gDs**2-(mc2*(ms2*(trg+t
     &   rc)-2*trs*cDs)+2*(trsgc-3*ms2*trg)*cDs)*gDs+mc2*ms2*(-2*trg*(4*
     &   cDs+ms2)-trsgc+mc2*trs))+2*mc2*cDg*gDs*((15*trg+13*trc)*gDs**3+
     &   cDs*((15*trg+13*trc)*gDs**2-2*(trsgc+mc2*trg)*gDs-2*mc2*ms2*trg
     &   )-2*(trsgc+mc2*trg)*gDs**2+mc2*(-3*trsgc+2*mc2*trs+ms2*trg)*gDs
     &   -mc2**2*ms2*trg)+mc2**2*gDs**2*(2*trg*gDs**2+2*trg*cDs*gDs-mc2*
     &   trsgc+mc2*ms2*trg)+2*cDg**4*(ms2*((5*trg+trc)*gDs+2*trg*(4*cDs+
     &   ms2))-trs*(4*gDs*(gDs+cDs)+3*mc2*ms2))-4*ms2*trs*cDg**5)/(ms*cD
     &   g*(2*cDg+mc2)**2*gDs)/2.0d+0+Bpmm

      Bpmm = LsB2*ms*(2*mc2**3*trg*cDs*gDs**4+mc2*cDg**2*gDs**2*(mc2*(2
     &   *(mc2*trs+ms2*trc)*gDs-ms2*(trsgc+3*mc2*trg))+cDs**2*(4*trc*gDs
     &   -2*trsgc+6*mc2*trg)+12*trg*cDs**3+3*mc2*(trsgc-ms2*trg)*cDs)+cD
     &   g**3*gDs*(3*mc2*cDs*(ms2*(trsgc+3*mc2*trg)-2*(mc2*trs+ms2*trc)*
     &   gDs)-8*trg*cDs**4-12*mc2*trg*cDs**3+2*mc2*(ms2*trg-3*trsgc)*cDs
     &   **2+mc2**2*ms2*(3*trsgc+ms2*trg))-mc2**2*cDg*gDs**3*(cDs*(2*trc
     &   *gDs-trsgc+mc2*trg)+8*trg*cDs**2+mc2*(trsgc-ms2*trg))+cDg**4*(2
     &   *mc2*(6*trs*cDs**2-3*mc2*ms2*trs+ms2**2*trc)*gDs+8*trg*cDs**4+4
     &   *(trsgc+ms2*trg)*cDs**3-8*mc2*ms2*trg*cDs**2-3*mc2*ms2*(trsgc+m
     &   s2*trg)*cDs+mc2*ms2**2*(mc2*trg-trsgc))+cDg**5*(6*mc2*ms2*trs*c
     &   Ds-8*trs*cDs**3))/(cDg**3*(cDs**2-mc2*ms2)*gDs)+Bpmm

      Bpmm = Bpmm-ls*ms*(-2*cDg**2*(mc2*(-8*trs*gDs**3+2*(trsgc-2*ms2*
     &   trs+ms2*trg)*gDs**2+ms2*(3*trsgc-2*ms2*trs-3*ms2*trg)*gDs+ms2**
     &   2*(trsgc-ms2*trg))+cDs*gDs*(-14*trc*gDs**2+(6*trsgc+4*mc2*trg+5
     &   *ms2*trc)*gDs+ms2*(3*trsgc+2*mc2*trg+9*ms2*trc))+2*(trs+trg)*cD
     &   s**2*gDs*(2*gDs+ms2))-mc2*cDg*gDs*(10*trc*gDs**3+(8*(-2*trs*cDs
     &   +4*trg*cDs+trsgc)+ms2*trc)*gDs**2+(4*ms2*(-2*trs*cDs+11*trg*cDs
     &   +3*trsgc)-4*trsgc*cDs-3*ms2**2*trc)*gDs+2*ms2*(-trsgc*cDs+7*ms2
     &   *trg*cDs+2*ms2*trsgc))+mc2*gDs**2*(trc*cDs*(-2*gDs**2+3*ms2*gDs
     &   +3*ms2**2)+2*mc2*(2*gDs+ms2)*(-2*trs*gDs+4*trg*gDs-trsgc+3*ms2*
     &   trg))+2*cDg**3*(2*gDs+ms2)*(2*cDs*((trs+trg+2*trc)*gDs+trsgc+ms
     &   2*trg)+gDs*(7*trc*gDs-3*trsgc-2*mc2*trg-9*ms2*trc))+4*cDg**4*(2
     &   *gDs+ms2)*(2*(trs+trg+trc)*gDs+2*trg*cDs+trsgc+ms2*trg)+8*trg*c
     &   Dg**5*(2*gDs+ms2))/(mc2*cDg*gDs*(2*gDs+ms2)**2)/4.0d+0

      Bpmm = ms*tr3c00fs*(-2*cDg**2*(mc2*(17*trc*gDs**3-2*(2*trsgc+2*mc
     &   2*trg+9*ms2*trc)*gDs**2-13*ms2**2*trc*gDs+ms2**2*(trsgc+mc2*trg
     &   ))+mc2*cDs*(4*(trg+2*trc)*gDs**2+2*ms2*(6*trs+7*trg+2*trc)*gDs+
     &   ms2*(trsgc+ms2*trg))+cDs**2*gDs*(-7*trc*gDs+3*trsgc+2*mc2*trg)+
     &   2*(trs+trg)*cDs**3*gDs)+2*cDg**3*(-mc2*(4*(trs+2*trg+2*trc)*gDs
     &   **2+(5*trsgc+6*ms2*trs+13*ms2*trg)*gDs+ms2*(3*trsgc-3*ms2*trs-2
     &   *ms2*trc))+cDs*gDs*(14*trc*gDs-6*trsgc-4*mc2*trg-9*ms2*trc)+2*c
     &   Ds**2*(2*trc*gDs+trsgc+ms2*trg))+2*mc2*cDg*gDs*(mc2*(4*(trg+trc
     &   )*gDs**2+(trsgc+3*ms2*trs+6*ms2*trg+2*ms2*trc)*gDs+ms2*(trsgc+m
     &   s2*trg))+cDs*(-17*trc*gDs**2+4*(trsgc+mc2*trg-2*ms2*trc)*gDs+2*
     &   ms2*(trsgc+mc2*trg))+cDs**2*(4*(trs+trg)*gDs+trsgc+ms2*trg))+mc
     &   2**2*gDs**2*(20*trc*gDs**2+(5*ms2*trc-2*(2*(trs+trg)*cDs+trsgc+
     &   2*mc2*trg))*gDs-2*ms2*(trg*(cDs+mc2)+trsgc)-2*trsgc*cDs-3*ms2**
     &   2*trc)+2*cDg**4*(7*trc*gDs**2+cDs*((6*(trs+trg)+8*trc)*gDs+4*(t
     &   rsgc+ms2*trg))-(3*trsgc+10*mc2*trg+9*ms2*trc)*gDs+4*trg*cDs**2-
     &   4*mc2*ms2*trg)+4*cDg**5*(2*(trs+trg+trc)*gDs+4*trg*cDs+trsgc+ms
     &   2*trg)+8*trg*cDg**6)/(mc2*cDg*gDs**3)/2.0d+0+Bpmm

      Bpmm = BfunX*(-ms2*cDg**4*gDs*(mc2*((2*(trs+trg)+3*trc)*gDs**2+ms
     &   2*(2*(trs+trg)+trc)*gDs-2*ms2**2*(trs+trg))+cDs*(12*trc*gDs**2+
     &   mc2*(4*trs+9*trc)*gDs+4*mc2*ms2*(trs+trg))+6*trc*cDs**2*(4*gDs+
     &   mc2)+12*trc*cDs**3)-mc2*ms2*cDg**2*gDs**2*(gDs*((20*trs+3*trg+7
     &   *trc)*gDs**2+ms2*(6*trs-5*trg+trc)*gDs-2*ms2**2*(trs+trg))+cDs*
     &   *2*((4*(trs+trg)+11*trc)*gDs+ms2*trc)+cDs*gDs*(24*trs*gDs+7*trg
     &   *gDs+17*trc*gDs-2*ms2*trs-2*ms2*trg+6*ms2*trc)+trc*cDs**3)-ms2*
     &   cDg**5*(-mc2*(trc*gDs**2+2*ms2*(2*(trs+trg)+trc)*gDs+2*ms2**2*(
     &   trs+trg+trc))+6*trc*cDs*gDs*(4*gDs+mc2)+24*trc*cDs**2*gDs)-mc2*
     &   ms2*cDg**3*gDs**2*((6*trs+trg+7*trc)*gDs**2+2*cDs*((trs-trg+5*t
     &   rc)*gDs+ms2*(trs-trg))+2*ms2*(trs+trg+2*trc)*gDs-(4*trs+8*trg+3
     &   *trc)*cDs**2+2*ms2**2*(trs+trg))-2*mc2*cDg*gDs**4*(2*(gDs+cDs)+
     &   ms2)*(2*trs*(cDs+2*ms2)*(2*(gDs+cDs)+ms2)+3*ms2*(trg+trc)*(gDs+
     &   cDs))-2*mc2*gDs**4*(2*(gDs+cDs)+ms2)**2*(trs*(cDs+ms2)*(2*(gDs+
     &   cDs)+ms2)+ms2*(trg+trc)*gDs)-12*ms2*trc*cDg**6*cDs*gDs)/(mc2*ms
     &   *cDg**3*gDs**3)/8.0d+0+Bpmm

      Bpmm = Bpmm-LsB1*ms*(cDg**2*gDs**2*(mc2*cDs*(12*trc*gDs**2+2*(-3
     &   *trsgc+3*ms2*trs-ms2*trg+5*ms2*trc)*gDs+ms2*(trsgc-5*ms2*trg))+
     &   4*cDs**2*(-2*trc*gDs**2+(trsgc+mc2*trg-2*ms2*trc)*gDs+ms2*(trsg
     &   c+2*mc2*trg))+mc2*ms2*(-4*trc*gDs**2+2*(trsgc+mc2*trg-2*ms2*trc
     &   )*gDs+ms2*(mc2*trg-trsgc))+cDs**3*(8*(trg+trc)*gDs-4*trsgc+8*ms
     &   2*trg))-cDg**3*gDs*(2*cDs**2*(4*trc*gDs**2-2*trsgc*gDs+2*ms2*(t
     &   rs+trg+4*trc)*gDs-3*ms2*trsgc+ms2**2*trg)+mc2*ms2*(4*trc*gDs**2
     &   +2*(ms2*(-trs-trg+trc)-trsgc)*gDs+ms2*(trsgc-ms2*trg))+ms2*cDs*
     &   (-4*trc*gDs**2+2*(trsgc+mc2*trg-2*ms2*trc)*gDs+ms2*(trsgc+3*mc2
     &   *trg)))-mc2*cDg*gDs**3*(mc2*(4*trc*gDs**2+2*(ms2*(3*trs-trg+trc
     &   )-trsgc)*gDs+ms2*(trsgc-5*ms2*trg))+cDs*(-12*trc*gDs**2+6*(trsg
     &   c+mc2*trg-2*ms2*trc)*gDs+ms2*(3*trsgc+5*mc2*trg))+4*cDs**2*(2*t
     &   rg*(gDs+ms2)+trc*gDs)+4*trg*cDs**3)+mc2*gDs**4*(-4*mc2*trc*gDs*
     &   *2+2*mc2*((trs+trg)*cDs+trsgc+mc2*trg-2*ms2*trc)*gDs+2*(mc2*trg
     &   -trsgc)*cDs**2+mc2*(trsgc+ms2*trg)*cDs+mc2*ms2*(3*trsgc+mc2*trg
     &   ))+2*ms2*cDg**4*cDs*(trc*gDs*(2*gDs+3*ms2)-trsgc*(gDs+ms2)))/(c
     &   Dg*(cDs-mc*ms)*(cDs+mc*ms)*gDs**3)

      Bpmm = Bpmm-B0csf*ms*(mc2*gDs*(mc2*(2*trc*gDs**2+(-trsgc-mc2*trs
     &   +ms2*(trg-trs))*gDs+ms2*(ms2+mc2)*trg)+cDs*(2*trc*gDs**2+(-trsg
     &   c+mc2*(-2*trs-trg+trc)+ms2*trc)*gDs+2*mc2*ms2*trg)+cDs**2*(2*tr
     &   c*gDs-trg*(2*gDs+ms2+mc2))-2*trg*cDs**3)+cDg*(2*mc2*cDs*(-2*trc
     &   *gDs**2+(trsgc+ms2*(trs-trg)+mc2*trs)*gDs-ms2*(trsgc+(ms2+mc2)*
     &   trg))+cDs**2*(-4*trc*gDs**2+2*(trsgc+mc2*(2*trs+trg-trc)-ms2*tr
     &   c)*gDs+mc2*(trsgc-4*ms2*trg)+ms2*trsgc)+2*cDs**3*(trg*(2*gDs+ms
     &   2+mc2)-2*trc*gDs+trsgc)+4*trg*cDs**4-mc2*ms2*(ms2+mc2)*trsgc)+c
     &   Dg**2*(ms2*cDs*(2*trc*gDs-trsgc+mc2*(2*trs-trg+trc)+ms2*trc)+mc
     &   2*ms2*(2*trc*gDs-trsgc+ms2*(trs+trg)+mc2*trs)-4*trs*cDs**3+2*(m
     &   s2*(-trs-trg+trc)-mc2*trs)*cDs**2))/(cDg*(cDs-mc*ms)*(cDs+mc*ms
     &   )*gDs)/2.0d+0

      Bpmm = tr3s002ft*(6*trs*(2*cDg+mc2)*cDs*((2*cDg+mc2)**2*gDs**2+ms
     &   2*cDg**3)+3*ms2*(trg+trc)*cDg**2*gDs*(mc2*(gDs+6*cDs)+10*cDg*cD
     &   s))/(ms*cDg**3*gDs)/2.0d+0+Bpmm

      Bpmm = 3.0d+0*ms*tr3c002fs*(trc*cDg**3*(gDs**2+2*ms2*gDs+2*ms2**2
     &   )+cDg*cDs*(3*trc*cDs-4*mc2*(trs+trg))*gDs**2-cDg**2*gDs*(3*trc*
     &   cDs*gDs+2*mc2*ms2*(trs+trg))+gDs**2*(2*mc2**2*(trs+trg)*gDs-trc
     &   *cDs**3))/(2.0d+0*cDg*gDs**3)+Bpmm

      Bpmm = 3.0d+0*ms*tr3c001fs*(-mc2*ms2*trc*cDs**2*gDs**2+2*mc2*cDg*
     &   gDs**2*(mc2*(trs+trg)*(2*gDs+ms2)-ms2*trc*cDs)-cDg**2*gDs*(8*mc
     &   2*(trs+trg)*cDs*(gDs+ms2)+mc2*ms2*trc*gDs-6*ms2*trc*cDs**2)+2*m
     &   s2*cDg**3*(3*trc*cDs*gDs+mc2*ms2*(trs+trg)))/(2.0d+0*mc2*cDg*gD
     &   s**3)+Bpmm

      Bpmm = epinv*(4*lp*xs*cDs+mc*(ms-ms*xs**2))*(2*trc*cDg*gDs-mc2*tr
     &   g*gDs-trsgc*cDg)*(mc2*gDs**2-2*cDg*cDs*gDs+ms2*cDg**2)/(mc*(xs*
     &   *2-1)*cDg**2*gDs**2)/4.0d+0+Bpmm

      Bpmm = ms*tr2fu*cDs*(2*trc*cDg*gDs-mc2*trg*gDs-trsgc*cDg)*(mc2*gD
     &   s**2-2*cDg*cDs*gDs+ms2*cDg**2)/(cDg**2*gDs**2)+Bpmm

      Bpmm = tr3s001ft*(3*ms*(trg+trc)*(2*ms2*cDg**3*cDs-gDs**2*(2*mc2*
     &   cDg*(gDs-3*cDs)+cDg**2*(gDs-7*cDs)+2*mc2**2*gDs))+6*ms*trs*gDs*
     &   (cDg**3*(gDs+5*cDs)+3*mc2*cDg**2*(gDs+cDs)+3*mc2**2*cDg*gDs+mc2
     &   **3*gDs))/(cDg**3*gDs)/2.0d+0+Bpmm

      Bpmm=Bpmm/zpmm

      return
      end

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen