text:
	SUBROUTINE GSCNAM  ( icolr, color, iret )
C************************************************************************
C* GSCNAM								*
C*									*
C* This subroutine defines the color corresponding to a color 		*
C* number by specifying a color name.  Maximum characters of		*
C* the color name = 40.  The colors are stored in coltbl.tbl.		*
C*									*
C* GSCNAM  ( ICOLR, COLOR, IRET )					*
C*									*
C* Input parameters:							*
C* 	ICOLR		INTEGER		Color number			*
C*	COLOR		CHAR*		Color name			*
C*									*
C* Output parameters:							*
C*	IRET		INTEGER	 	Return code			*
C**									*
C* Log:									*
C* I. Graffman/RDS	 8/85						*
C* M. desJardins/GSFC	 5/88	Documentation				*
C* D. Keiser/GSC	12/95	Removed reference to GEMTBL in comments	*
C* L. Williams/EAI	 3/96   Removed kolor and Call to ST_LCUC	*
C************************************************************************
	INCLUDE 	'ERROR.PRM'
C
	CHARACTER*(*) 	color
C-----------------------------------------------------------------------
C*	Set the color by name.
C
	CALL DSCNAM  ( icolr, color, iret )
C*
	RETURN
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
-makelib ies_lib/xil_defaultlib -sv \
  "G:/Xilinx/Vivado/2018.3/data/ip/xpm/xpm_memory/hdl/xpm_memory.sv" \
-endlib
-makelib ies_lib/xpm \
  "G:/Xilinx/Vivado/2018.3/data/ip/xpm/xpm_VCOMP.vhd" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../ipstatic/hdl/csi2_parser.vhd" \
  "../../../ipstatic/hdl/lane_align.vhd" \
  "../../../ipstatic/hdl/lane_merge.vhd" \
  "../../../ipstatic/hdl/csi_to_axis_v1_0.vhd" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/src/csi_to_axis_0/sim/csi_to_axis_0.vhd" \
  "../../../ipstatic/hdl/phy_clock_system.vhd" \
  "../../../ipstatic/hdl/line_if.vhd" \
  "../../../ipstatic/hdl/csi2_d_phy_rx.vhd" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/src/csi2_d_phy_rx_0/sim/csi2_d_phy_rx_0.vhd" \
-endlib
-makelib ies_lib/blk_mem_gen_v8_4_2 \
  "../../../ipstatic/simulation/blk_mem_gen_v8_4.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/src/blk_mem_gen_0/sim/blk_mem_gen_0.v" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/src/blk_mem_gen_1/sim/blk_mem_gen_1.v" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/sim/Driver_Bayer_To_RGB.v" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/sim/Driver_Csi_To_Dvp.v" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/sim/Driver_MIPI.v" \
  "../../../../Digital-Recognition.srcs/sources_1/ip/Driver_MIPI_0/sim/Driver_MIPI_0.v" \
-endlib
-makelib ies_lib/xil_defaultlib \
  glbl.v
-endlib


text:
!> Module holding basic KIND parameters
!> This module imitates the naming of dl_esm_inf
!> It uses C interoperable types for safety with Regent.
module kind_params_mod
  use iso_c_binding
  implicit none

  public

  !> Douple precision kind parameter
  integer, parameter :: GO_WP = c_double

  ! Kind type for double precision
  integer, parameter :: GO_DP = c_double

end module kind_params_mod

max_stars_repo_name: stfc/PSycloneBench

max_issues_repo_name: stfc/PSycloneBench

max_forks_repo_name: stfc/PSycloneBench


text:
C    Copyright(C) 1999-2020 National Technology & Engineering Solutions
C    of Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C    NTESS, the U.S. Government retains certain rights in this software.
C
C    See packages/seacas/LICENSE for details

      SUBROUTINE DONRM2 (COORD, LTNESS, MAP, DIRCOS, TEMP,
     *    NSEG, NUMNIQ, NUMNP)
      DIMENSION COORD(NUMNP, *), LTNESS(2,*), MAP(*), DIRCOS(4,*),
     *    TEMP(2,*)

      DO 10 ISEG = 1, NSEG
          XI = COORD( LTNESS(1,ISEG),1 )
          YI = COORD( LTNESS(1,ISEG),2 )

          XJ = COORD( LTNESS(2,ISEG),1 )
          YJ = COORD( LTNESS(2,ISEG),2 )

          DX = XI - XJ
          DY = YI - YJ
          RMAG = SQRT ( DX**2 + DY**2)

          TEMP(1,ISEG) = -DY / RMAG
          TEMP(2,ISEG) =  DX / RMAG

   10 CONTINUE

      DO 20 I=1,NUMNIQ
          DIRCOS(1,I) = 0.0
          DIRCOS(2,I) = 0.0
   20 CONTINUE

      DO 40 ISEG = 1, NSEG
          DO 30 J = 1, 2
              MISEG = MAP( 2 * (ISEG-1) + J )
              DIRCOS(1,MISEG) = DIRCOS(1,MISEG) + TEMP(1,ISEG)
              DIRCOS(2,MISEG) = DIRCOS(2,MISEG) + TEMP(2,ISEG)
   30     CONTINUE
   40 CONTINUE

C ... NORMALIZE ALL DIRECTION COSINES

      DO 50 I = 1, NUMNIQ
          A = DIRCOS(1,I)
          B = DIRCOS(2,I)
          R = SQRT(A**2 + B**2)
          DIRCOS(1,I) = A/R
          DIRCOS(2,I) = B/R
   50 CONTINUE
      RETURN
      END

max_stars_repo_name: jschueller/seacas

max_issues_repo_name: jschueller/seacas

max_forks_repo_name: jschueller/seacas


text:
!	Define the latent vector dimension dimlatin and the transformed latent vector dimension dimlatout.,
!	custom is for a special linear transformation.
subroutine getdim(dimlatin,dimlatout,custom)
implicit none
integer::io
integer,intent(out)::dimlatin,dimlatout
logical,intent(out)::custom
namelist/dimension/dimlatin,dimlatout,custom
dimlatin=1
dimlatout=1

custom=.false.
read(*,nml=dimension,iostat=io)
if(io/=0) stop "Dimension or integration data not read successfully."
dimlatin=max(1,dimlatin)
dimlatout=max(1,dimlatout)
return
end subroutine getdim

max_stars_repo_name: EducationalTestingService/MIRT

max_issues_repo_name: EducationalTestingService/MIRT

max_forks_repo_name: EducationalTestingService/MIRT


text:
      WRITE (6,1) (I,I = 1,10)
    1 FORMAT (4(1X,I0,","),1X,I0)
      END

max_issues_repo_name: seanwallawalla-forks/RosettaCodeData

max_forks_repo_name: seanwallawalla-forks/RosettaCodeData


text:
      subroutine tffscalc(kdp,df,iqcol,lfp,
     $     nqcola,nqcola1,ibegin,
     $     r,rp,rstab,nstab,residual1,
     $     zcal,wcal,parallel,lout,error)
      use tfstk
      use ffs, only:ndim,nlat,flv,maxcond,ffs_bound,nvevx,nelvx
      use ffs_flag
      use ffs_pointer
      use ffs_fit
      use ffs_wake
      use tffitcode
      use tfshare
      use tfcsi,only:lfno
      use eeval
      use macmath
      use mathfun
      implicit none
c      include 'DEBUG.inc'
      type (sad_descriptor) kx
      type (ffs_bound) fbound,ibound
      integer*4 ,intent(in):: nqcola,nqcola1,
     $     kdp(maxcond),iqcol(maxcond),lfp(2,maxcond)
      integer*4 ,intent(out):: ibegin
      integer*4 i1,i2,i3,i,ii,j,iter,kt,iq,l,maxf,
     $     ie,ie1,iv,nstab,lout,irtc
      real*8 ,intent(out):: df(maxcond),r,rp,residual1(-ndimmax:ndimmax)
      logical*4 ,intent(in):: parallel
      logical*4 ,intent(out):: error,zcal,wcal
      real*8 anux0,anuy0,anux0h,anuy0h,anuxi,anuyi,anuxih,anuyih,
     $     rw,drw,rstab,wi,anusumi,anusum0,anudiffi,anudiff0
      logical*4 fam,beg,zerores
      integer*4 irw,isw,ipr,ifb,ife,idir,
     $     jjfam(-nfam:nfam),ifpe,ntfun
      integer*4, external :: fork_worker,waitpid,
     $     itfdownlevel,itfuplevel,itgetfpe
      integer*8 iutm,jb
      integer*4 , parameter:: ivoid=9999
      integer*8 ,save :: iprolog=0,iepilog=0,imr=0,inr=0,isl=0
      associate (nvar=>flv%nvar)
c     begin initialize for preventing compiler warning
      anux0=0.d0
      anuy0=0.d0
      anux0h=0.d0
      anuy0h=0.d0
      anusum0=0.d0
      anudiff0=0.d0
      iutm=0
c     end   initialize for preventing compiler warning
      if(iprolog .eq. 0)then
        iprolog=ktfsymbolz('OpticsProlog',12)
        iepilog=ktfsymbolz('OpticsEpilog',12)
        imr=    ktfsymbolz('MatchingResidual',16)-4
        inr=    ktfsymbolz('NetResidual',11)-4
        isl=    ktfsymbolz('StabilityLevel',14)-4
      endif
      ifpe=itgetfpe()
      call tclrfpe
c      call tfmemcheckprint('tffscalc-before-prolog',.true.,irtc)
      l=itfuplevel()
      kx=tfeeval(dfromk(ktfsymbol+iprolog),.true.,irtc)
      l=itfdownlevel()
      if(irtc .ne. 0)then
        if(irtc .gt. 0)then
          if(ierrorprint .ne. 0)then
            call tfaddmessage(' ',0,lfno)
          endif
          call tclrfpe
          call termes(lfno,'Error in OpticsProlog.',' ')
        endif
        error=.true.
        return
      endif
      call tsetfpe(ifpe)
      if(zcal)then
        call tfgeo(latt,geomet .or. .not. fitflg)
      endif
      call tffsbound(fbound)
      optstat(nfam1:nfam)%over=.false.
      optstat(nfam1:nfam)%stabx=.true.
      optstat(nfam1:nfam)%staby=.true.
      optstat(nfam1:nfam)%stabz=.true.
      jjfam(nfam1:nfam)=ivoid
      ntfun=merge(ntwissfun,mfitdetr,orbitcal .or. calc6d)
      beg=ibegin .gt. fbound%lb
      if(beg)then
        twiss(ibegin,0,1:ntfun)
     $       =utwiss(1:ntfun,0,itwissp(ibegin))
        fbound%fb=0.d0
      else
        call twmov(1,twiss,nlat,ndim,.true.)
        twiss(1,0,mfitnx)=0.d0
        twiss(1,0,mfitny)=0.d0
        if(orbitcal)then
          twiss(1,0,mfitddp)=twiss(1,0,mfitddp)+dp(0)
        endif
        if(fbound%lb .gt. 1)then
          twiss(fbound%lb,0,1:ntfun)=twiss(1,0,1:ntfun)
        endif
        ibegin=fbound%lb
      endif
      ibound=fbound
      ibound%lb=ibegin
      if(wake)then
        call tffswake(ibound,beg)
      else
        call qcell1(ibound,0,optstat(0),.false.,chgini,lout)
        call tffssetutwiss(0,nlat,fbound,beg,.true.,.true.)
        if(cell)then
          anux0=aint(twiss(nlat,0,mfitnx)/pi2)
          anuy0=aint(twiss(nlat,0,mfitny)/pi2)
          anux0h=aint(twiss(nlat,0,mfitnx)/pi)
          anuy0h=aint(twiss(nlat,0,mfitny)/pi)
          anusum0=aint((twiss(nlat,0,mfitnx)+twiss(nlat,0,mfitny))/pi2)
          anudiff0=tfloor(
     $         (twiss(nlat,0,mfitnx)-twiss(nlat,0,mfitny))/pi2)
c          anudiff0=twiss(nlat,0,mfitnx)/pi2-
c     $         aint(twiss(nlat,0,mfitnx)/pi2)-
c     $         twiss(nlat,0,mfitny)/pi2+
c     $         aint(twiss(nlat,0,mfitny)/pi2)
          if(optstat(0)%stabx .and. optstat(0)%staby
     $         .and. optstat(0)%stabz .or. chgini)then
            call twmov(1,twiss,nlat,ndim,.false.)
          endif
        endif
        if(nfam .ne. 0)then
          if(parallel)then
            irtc=0
            iutm=ktfallocshared((2*nfam+1)*4)
c            iutm=mapalloc8(rlist(1),(2*nfam+1)*4,8,irtc)
            if(irtc .eq. 0)then
              ipr=fork_worker()
            else
              ipr=-1
            endif
          else
            ipr=-1
          endif
          if(ipr .gt. 0)then
            idir=-1
            ifb=-1
            ife=-nfr
          else
            idir=1
            ifb=1
            ife=nfr
          endif
 2        i1=0
          i2=0
          i3=0
          fam=.false.
 1        do ii=ifb,ife,idir
            if(fam)then
              if(kfam(ii) .eq. 0)then
                cycle
              elseif(ipr .gt. 0 .and. jfam(ii) .ge. 0 .or.
     $               ipr .eq. 0 .and. jfam(ii) .lt. 0)then
                cycle
              endif
              i3=i2
              i2=jfam(ii)
              i1=ii
            else
              i3=i2
              i2=i1
              i1=ii
            endif
            if(optstat(i2)%over)then
              utwiss(1,ii,1:nut)=utwiss(1,i2,1:nut)
              utwiss(mfitddp,ii,1:nut)=utwiss(mfitddp,ii,1:nut)+dp(ii)
              optstat(ii)%over=.true.
              optstat(ii)%stabx=.false.
              optstat(ii)%staby=.false.
              optstat(ii)%stabz=.false.
              optstat(ii)%tracex=optstat(i2)%tracex
              optstat(ii)%tracey=optstat(i2)%tracey
              optstat(ii)%tracez=optstat(i2)%tracez
            else
              if(beg)then
                twiss(ibegin,1,1:ntfun)
     $               =utwiss(1:ntfun,ii,itwissp(ibegin))
              else
                twiss(1,1,:)=utwiss(:,0,1)
                if(inicond)then
                  if(uini(mfitbx,ii) .gt. 0.d0)then
                    twiss(1,1,:)=uini(:,ii)
                  endif
                  twiss(1,1,mfitdx:mfitddp)=
     $                 utwiss(mfitdx:mfitddp,0,1)+
     $                 uini(mfitdx:mfitddp,ii)
                else
                  if(fam)then
                    twiss(1,1,mfitdx:mfitdpy )=
     $                   utwiss(mfitdx:mfitdpy ,i2,1)+dfam(1:4,ii)
                  else
c                    call tgetphysdispu(utwiss(1,i2,1),physd)
c                    physd1=(utwiss(mfitdx:mfitdpy,i2,1)
c     $                     -utwiss(mfitdx:mfitdpy,i3,1))/(dp(i2)-dp(i3))
c                    do i=1,4
c                      if(abs(physd(i)) .gt. abs(physd1(i)))then
c                        physd(i)=physd1(i)
c                      endif
c                    enddo
                    twiss(1,1,mfitdx:mfitdpy)=
     $                   utwiss(mfitdx:mfitdpy,i2,1)
c     $                   +(dp(ii)-dp(i2))*physd
c                    if(ii .eq. nfr)then
c                      write(*,'(a,1p10g12.4)')'tffscalc ',
c     $                     twiss(1,1,mfitdx:mfitdpy),
c     $                     utwiss(mfitex:mfitepy,i2,1),dp(ii),dp(i2)
c                    endif
                  endif
                  twiss(1,1,mfitdz )=utwiss(mfitdz ,i2,1)
                  twiss(1,1,mfitddp)=twiss(1,0,mfitddp)+dp(ii)
                endif
                twiss(1,1,mfitnx)=0.d0
                twiss(1,1,mfitny)=0.d0
                if(fbound%lb .gt. 1)then
                  twiss(fbound%lb,1,1:ntfun)=twiss(1,1,1:ntfun)
                endif
              endif
              call qcell1(ibound,1,optstat(ii),fam,chgini,lout)
              call tffssetutwiss(ii,nlat,fbound,beg,.true.,.true.)
            endif
            if(cell)then
              anuxih=aint(twiss(nlat,1,mfitnx)/pi)
              anuyih=aint(twiss(nlat,1,mfitny)/pi)
              anuxi=aint(twiss(nlat,1,mfitnx)/pi2)
              anuyi=aint(twiss(nlat,1,mfitny)/pi2)
              anusumi=aint((twiss(nlat,1,mfitnx)
     $             +twiss(nlat,1,mfitny))/pi2)
              anudiffi=tfloor((twiss(nlat,1,mfitnx)
     $             -twiss(nlat,1,mfitny))/pi2)
c              anudiffi=twiss(nlat,1,mfitnx)/pi2-
c     $             aint(twiss(nlat,1,mfitnx)/pi2)-
c     $             twiss(nlat,1,mfitny)/pi2+
c     $             aint(twiss(nlat,1,mfitny)/pi2)
c              write(*,*)'tffscalc ',anudiffi,anudiff0
              optstat(ii)%stabx=optstat(ii)%stabx .and. (fam .or.
     $             (.not. intres .or. anuxi .eq. anux0) .and.
     $             (.not. halfres .or. anuxih .eq. anux0h) .and.
     $             (.not. sumres .or. anusumi .eq. anusum0) .and.
     $             (.not. diffres .or. anudiffi .eq. anudiff0))
              optstat(ii)%staby=optstat(ii)%staby .and. (fam .or.
     $             (.not. intres .or. anuyi .eq. anuy0) .and.
     $             (.not. halfres .or. anuyih .eq. anuy0h) .and.
     $             (.not. sumres .or. anusumi .eq. anusum0) .and.
     $             (.not. diffres .or. anudiffi .eq. anudiff0))
            endif
          enddo
          if(ipr .eq. -1 .and. .not. fam .and. idir .eq. 1)then
            idir=-1
            ifb=-1
            ife=-nfr
            go to 2
          endif
          if(.not. fam .and. nfam .gt. nfr)then
            fam=.true.
            ifb=nfam1
            ife=nfam
            idir=1
            go to 1
          endif
          if(ipr .gt. 0)then
            irw=0
            do while(irw .ne. ipr)
              irw=waitpid(-1,isw)
            enddo
            if(isw .ne. 0)then
              call termes(lfno,
     1             '?Error in parallel process.',' ')
            endif
            do i=nfam1,nfam
              if(i .gt. 0 .and. i .le. nfr .or.
     $             kfam(i) .ne. 0 .and. jfam(i) .ge. 0)then
                jb=iutm+(i+nfam)*4
                optstat(i)%stabx=ilist(1,jb+1)
                optstat(i)%staby=ilist(1,jb+2)
                optstat(i)%tracex=rlist(jb+3)
                optstat(i)%tracey=rlist(jb+4)
              endif
            enddo
          elseif(ipr .eq. 0)then
            do i=nfam1,nfam
              if(i .gt. 0 .and. i .le. nfr .or.
     $             kfam(i) .ne. 0 .and. jfam(i) .ge. 0)then
                jb=iutm+(i+nfam)*4
                ilist(1,jb+1)=optstat(i)%stabx
                ilist(1,jb+2)=optstat(i)%staby
                rlist(jb+3)=optstat(i)%tracex
                rlist(jb+4)=optstat(i)%tracey
              endif
            enddo
            stop
          endif
          if(ipr .gt. 0)then
            call tfreeshared(iutm)
          endif
        endif
      endif
      call tdfun(iqcol,lfp,nqcola,nqcola1,kdp,df,error)
      if(error)then
        call termes(lfno,
     1         '?Too many fit conditions.',' ')
        return
      endif
      if(wcal)then
        iter=1
        zcal=.false.
        do i=1,nvar
          kt=idtypec(nelvx(nvevx(i)%ivarele)%klp)
          if(kt .ne. 6 .and. kt .ne. 8 .and. kt .ne. 10
     $         .and. kt .ne. 12)then
            zcal=.true.
            exit
          endif
        enddo
        if(fitflg)then
          if(nvar .le. 0)then
            call termes(lfno,'?No variable.',' ')
            error=.true.
          endif
          if(.not. cell .and. .not. geomet)then
            do i=ibegin,nlat-1
              ii=icomp(i)
              ie=iele1(ii)
              ie1=iele1(i)
              do j=1,nvar
                iv=nvevx(j)%ivvar
                if(iv .eq. nelvx(ie)%ival .and. nvevx(j)%ivarele .eq. ie
     $               .and. (nvevx(j)%ivcomp .eq. 0 .or.
     $               nvevx(j)%ivcomp .eq. ii))then
                  ibegin=i
                  go to 1023
                elseif(iv .ne. nelvx(ie)%ival .and.
     $                 nvevx(j)%ivarele .eq. ie1
     $                 .and. (nvevx(j)%ivcomp .eq. 0 .or.
     $                 nvevx(j)%ivcomp .eq. i))then
                  ibegin=i
                  go to 1023
                endif
                if(nvevx(j)%ivarele .gt. ie)then
                  exit
                endif
              enddo
            enddo
            ibegin=nlat
          endif
        endif
1023    etamax=1.d-2
        if(cell)then
          maxf=fbound%le
        else
          maxf=fbound%lb
          do i=1,nfcol
            maxf=max(maxf,
     $           flv%ifitp(flv%kfitp(i)),flv%ifitp1(flv%kfitp(i)))
          enddo
          if(pos(maxf) .le. pos(1))then
            maxf=fbound%le
          endif
        endif
        do l=fbound%lb,maxf
          etamax=max(abs(twiss(l,0,mfitex)),
     $         abs(twiss(l,0,mfitey)),etamax)
        enddo
        avebeta=(pos(maxf)-pos(1))/
     $       max(twiss(maxf,0,mfitnx),twiss(maxf,0,mfitny))
      endif
      call twfit(flv%kfit,
     1     flv%ifitp,flv%kfitp,kdp,nqcola,iqcol,maxf,wcal)
      wcal=.false.
      rw=0.d0
      residual1(nfam1:nfam)=0.d0
      wsum=0.d0
      do i=1,nqcola
        if(kdp(i) .ne. 0)then
          iq=iqcol(i)
          wi=(offmw/2.d0/
     $         sqrt(dble(max(1,abs(flv%mfitp(flv%kfitp(iq)))))))**2
          wsum=wsum+wi
          wiq(i)=wiq(i)*wi**(1.d0/wexponent)
        else
          wsum=wsum+1.d0
        endif
        df(i)=df(i)*wiq(i)
        if(df(i) .ne. 0.d0)then
          drw=min(1.d50,max(1.d-50,abs(df(i))))**wexponent
          rw=rw+drw
          residual1(kdp(i))=residual1(kdp(i))+drw
        endif
      enddo
      r=merge(wsum*(max(rw,1.d-50)/wsum)**(2.d0/wexponent),
     $     0.d0,rw .gt. 0.d0)
      rp=r
      nstab=0
      if(cell)then
        if(rstab .eq. 0.d0)then
          rstab=10.d0
          do while(rstab .lt. r)
            rstab=rstab*10.d0
          enddo
        endif
        cellstab=.true.
        zerores=.true.
        do i=nfam1,nfam
          if(residual1(i) .ne. 0.d0)then
            zerores=.false.
            if(.not. optstat(i)%stabx)then
              nstab=nstab+1
              residual1(i)=residual1(i)+rstab
              r=r+rstab
              cellstab=.false.
            endif
            if(.not. optstat(i)%staby)then
              nstab=nstab+1
              residual1(i)=residual1(i)+rstab
              r=r+rstab
              cellstab=.false.
            endif
          endif
        enddo
        if(zerores)then
          do i=nfam1,nfam
            cellstab=cellstab .and.
     $           optstat(i)%stabx .and. optstat(i)%staby
          enddo
        endif
      else
        do i=nfam1,nfam
          if(.not. optstat(i)%stabx)then
            nstab=nstab+2
            residual1(i)=residual1(i)+20.d0
          endif
        enddo
      endif
      rlist(imr)=r
      rlist(inr)=rp
      rlist(isl)=dble(nstab)
      ifpe=itgetfpe()
      call tclrfpe
      l=itfuplevel()
      kx=tfeeval(dfromk(ktfsymbol+iepilog),.true.,irtc)
      l=itfdownlevel()
      if(irtc .ne. 0)then
        if(irtc .gt. 0)then
          if(ierrorprint .ne. 0)then
            call tfaddmessage(' ',0,lfno)
          endif
          call tclrfpe
          call termes(lfno,'Error in OpticsEpilog.',' ')
        endif
        error=.true.
      endif
c      call tfevals('Print["PROF: ",LINE["PROFILE","Q1"]]',kxx,irtc)
      call tsetfpe(ifpe)
      return
      end associate
      end

      subroutine twfit(kfit,
     1     ifitp,kfitp,kdp,nqcola,iqcol,maxf,wcal)
      use tfstk
      use ffs, only:emx,emy,dpmax,coumin,emminv
      use ffs_pointer
      use ffs_fit
      use ffs_flag, only:cell
      use tffitcode
      use eeval
      implicit none
      integer*4 ,intent(in):: nqcola,maxf,
     $     kfit(*),ifitp(*),kfitp(*),kdp(*),iqcol(nqcola)
      integer*4 i,j,k,iq
      real*8 coum,emxx,emyy,dpm,coup,em
      integer*4 itfuplevel, level,irtc,idp
      character*16 name
      logical*4 ,intent(in):: wcal
      type (sad_descriptor) kx
      type (sad_descriptor) ,save::kfv
      data kfv%k /0/
      type (sad_dlist), pointer , save::klv
      type (sad_rlist), pointer , save::klid
      integer*8 , save:: ifvloc,ifvfun
      real*8 , parameter :: almin=1.d0
      if(kfv%k .eq. 0)then
        kfv=kxadaloc(0,4,klv)
        klv%head=dtfcopy(kxsymbolz('`FitWeight',10))
        ifvloc=ktsalocb(0,'                ',MAXPNAME+8)
        ifvfun=ktsalocb(0,'        ',MAXPNAME)
        klv%body(1)=ktfstring+ifvloc
        klv%body(2)=ktfstring+ifvfun
        klv%dbody(3)=kxraaloc(0,2,klid)
        klv%body(4)=0
      endif
      em=max(emminv,abs(emx)+abs(emy))
      coum=min(1.d0,
     $     max(coumin,0.01d0,
     $     1.d0/(abs(emx/max(emminv,emy))+abs(emy/max(emminv,emx)))))
      coum=coum/(1.d0+coum)
      emxx=max(emx,coum*em)
      emyy=max(emy,coum*em)
      dpm=max(0.001d0,dpmax)
      coup=sqrt(abs(emxx/emyy)+abs(emyy/emxx))
      if(wcal)then
        do i=1,nfcol
          k=kfit(kfitp(i))
          j=ifitp(kfitp(i))
          select case(k)
          case (mfitex,mfitpex)
            wfit(i)=dpm/sqrt(emxx*twiss(j,0,mfitbx))
          case (mfitepx,mfitpepx)
            wfit(i)=dpm*sqrt(twiss(j,0,mfitbx)
     1           /(1.d0+twiss(j,0,mfitax)**2)/emxx)
          case (mfitey,mfitpey)
            wfit(i)=dpm/sqrt(emyy*twiss(j,0,mfitby))
          case (mfitepy,mfitpepy)
            wfit(i)=dpm*sqrt(twiss(j,0,mfitby)
     1           /(1.d0+twiss(j,0,mfitay)**2)/emyy)
          case (mfitr1)
            wfit(i)=coup*sqrt(twiss(j,0,mfitbx)/twiss(j,0,mfitby))
          case (mfitr2)
            wfit(i)=coup/sqrt(twiss(j,0,mfitbx)*twiss(j,0,mfitby))
          case (mfitr3)
            wfit(i)=coup*sqrt(twiss(j,0,mfitbx)*twiss(j,0,mfitby))
          case (mfitr4)
            wfit(i)=coup*sqrt(twiss(j,0,mfitby)/twiss(j,0,mfitbx))
          case (mfitdx)
            wfit(i)=1.d0/sqrt(twiss(j,0,mfitbx)*em)
          case (mfitdpx)
            wfit(i)=sqrt(twiss(j,0,mfitbx)/
     $           (1.d0+twiss(j,0,mfitax)**2)/em)
          case (mfitdy)
            wfit(i)=1.d0/sqrt(twiss(j,0,mfitby)*em)
          case (mfitdpy)
            wfit(i)=sqrt(twiss(j,0,mfitby)/
     $           (1.d0+twiss(j,0,mfitay)**2)/em)
          case (mfitleng)
            wfit(i)=1.d0/(pos(maxf)-pos(1))*
     $           max(twiss(maxf,0,mfitnx),twiss(maxf,0,mfitny))
          case (mfitdz)
            wfit(i)=0.01d0*sqrt(
     $           max(twiss(maxf,0,mfitnx),twiss(maxf,0,mfitny))
     $           /em/max(almin,pos(maxf)-pos(1)))
          case (mfitgx,mfitgy,mfitgz)
            wfit(i)=0.01d0*sqrt(
     $           max(twiss(maxf,0,mfitnx),twiss(maxf,0,mfitny))
     $           /em/max(almin,pos(maxf)-pos(1)))
          case (mfitchi1,mfitchi2,mfitchi3)
            wfit(i)=1.d3
          case default
            wfit(i)=1.d0
          end select
          if(cell .and.
     $         kfitp(i) .gt. nfc0 .and. kfitp(i) .le. nfc0+4)then
            wfit(i)=wfit(i)*0.7d0
          endif
        enddo
      endif
      do iq=1,nqcola
        i=iqcol(iq)
        if(i .gt. 0)then
          k=kfit(kfitp(i))
          j=ifitp(kfitp(i))
          idp=kdp(iq)
          wiq(iq)=wfit(i)
          call elname(j,name)
          call tfpadstr(name,ifvloc+1,len_trim(name))
          ilist(1,ifvloc)=len_trim(name)
          call tfpadstr(nlist(k),ifvfun+1,len_trim(nlist(k)))
          ilist(1,ifvfun)=len_trim(nlist(k))
          klid%rbody(1)=merge(dble(iuid(idp)),dble(kfam(idp)),
     $         inicond)
          klid%rbody(2)=dp(idp)
          klv%rbody(4)=wfit(i)
          call tclrfpe
          level=itfuplevel()
          kx=tfleval(klv,.true.,irtc)
          call tfconnect(kx,irtc)
          if(irtc .ne. 0)then
            if(ierrorprint .ne. 0)then
              call tfaddmessage(' ',2,6)
            endif
            call termes(6,'Error in FitWeight '//
     $           nlist(k)//' at '//name,' ')
          elseif(ktfrealq(kx,wiq(iq)))then
          endif
        else
          wiq(iq)=1.d0
        endif
      enddo
      return
      end

      subroutine tffsbound(fbound)
      use tfstk
      use ffs, only:nlat,ffs_bound
      use ffs_pointer
      use mackw
      implicit none
      type (ffs_bound) ,intent(out):: fbound
      call tffsbound1(1,nlat,fbound)
      return
      end

      subroutine tffsbound1(lb,le,fbound)
      use tfstk
      use ffs, only:nlat,ffs_bound
      use ffs_pointer
      use mackw
      implicit none
      type (ffs_bound) ,intent(out):: fbound
      integer*4 ,intent(in):: lb,le
      integer*4 le1
      real*8 xnlat,offset,tffsmarkoffset
      xnlat=nlat
      offset=tffsmarkoffset(lb)
      if(offset .ne. 0.d0)then
        offset=max(1.d0,min(xnlat,lb+offset))
        fbound%lb=int(offset)
        fbound%fb=offset-fbound%lb
      else
        fbound%lb=lb
        fbound%fb=0.d0
      endif
      le1=le
      if(le .eq. nlat)then
        if(idtypec(nlat-1) .eq. icMARK)then
          le1=le1-1
        else
          fbound%le=le1
          fbound%fe=0.d0
          return
        endif
      endif
      offset=tffsmarkoffset(le1)
      if(offset .ne. 0.d0)then
        offset=max(1.d0,min(xnlat,le1+offset))
        fbound%le=int(offset)
        fbound%fe=offset-fbound%le
      else
        fbound%le=le1
        fbound%fe=0.d0
      endif
      return
      end

      real*8 function tffsmarkoffset(lp)
      use tfstk
      use ffs_pointer
      use mackw
      use tmacro, only:nlat
      implicit none
      type (sad_comp), pointer :: cmp
      integer*4 ,intent(in):: lp
      integer*4 k
      if(lp .ge. nlat)then
        tffsmarkoffset=0.d0
        return
      endif
      call compelc(lp,cmp)
      k=kytbl(kwOFFSET,idtype(cmp%id))
      if(k .gt. 0)then
        tffsmarkoffset=cmp%value(k)
        if(tffsmarkoffset .ne. 0.d0)then
          tffsmarkoffset=(tffsmarkoffset-.5d0)*direlc(lp)+.5d0
        endif
      else
        tffsmarkoffset=0
      endif
      return
      end

      real*8 function tffselmoffset(l)
      use tfstk
      use ffs, only:nlat
      use ffs_pointer
      use mackw
      use tfcsi, only:icslfno
      implicit none
      integer*4 ,intent(in):: l
      integer*4 lm,nm,lx,nmmax
      parameter (nmmax=256)
      real*8 offset,xp,xe,tffsmarkoffset
      nm=0
      xp=l
      if(l .ne. nlat .and. kytbl(kwOFFSET,idtypec(l)) .ne. 0)then
        xe=nlat
        lm=l
        do
          offset=tffsmarkoffset(lm)
          if(offset .ne. 0.d0)then
            xp=offset+lm
            if(xp .ge. 1.d0 .and. xp .le. xe)then
              lx=int(xp)
              if(idtypec(lx) .eq. icMARK)then
                nm=nm+1
                if(nm .lt. nmmax)then
                  lm=lx
                  cycle
                else
                  call termes(icslfno(),
     $                 '?Recursive OFFSET in',pnamec(l))
                endif
              endif
            endif
          endif
          exit
        enddo
      endif
      tffselmoffset=xp
      return
      end

      subroutine tffssetutwiss(idp,nlat,fbound,beg,begin,end)
      use tfstk
      use ffs, only:ffs_bound
      use ffs_pointer
      use tffitcode
      use mackw
      implicit none
      type (ffs_bound) ,intent(in):: fbound
      integer*4 ,intent(in):: nlat,idp
      integer*4 jdp,j,jp,le1
      logical*4 ,intent(in):: beg,begin,end
      jdp=min(1,abs(idp))
      do concurrent (j=fbound%lb:fbound%le)
        jp=itwissp(j)
        if(jp .gt. 0)then
          utwiss(:,idp,jp)=twiss(j,jdp,:)
        endif
      enddo
      jp=itwissp(1)
      if(begin)then
        if(jp .gt. 0 .and. .not. beg)then
          twiss(1,jdp,:)=twiss(fbound%lb,jdp,:)
          utwiss(:,idp,jp)=twiss(fbound%lb,jdp,:)
        endif
      endif
      if(end)then
        if(idtypec(nlat-1) .eq. icMARK)then
          le1=merge(fbound%le,fbound%le+1,fbound%fe .eq. 0.d0)
          twiss(nlat-1,jdp,:)=twiss(le1,jdp,:)
          twiss(nlat,jdp,:)=twiss(le1,jdp,:)
          jp=itwissp(nlat-1)
          if(jp .gt. 0)then
            utwiss(:,idp,jp)=twiss(le1,jdp,:)
          endif
          jp=itwissp(nlat)
          if(jp .gt. 0)then
            utwiss(:,idp,jp)=twiss(le1,jdp,:)
          endif
        endif
      endif
      return
      end

max_stars_repo_name: noboruatkek/SAD

max_issues_repo_name: noboruatkek/SAD

max_forks_repo_name: noboruatkek/SAD


text:
c this file contains the following subroutines:
c
c readharmfile: reads a file with the harmonic expansion of the potential
c               and density.
c all the relevant stuff ends up in a common block that only needs to be seen by
c the routines in this file (I think!). This routine must be called first.
c
c dens(r,z): density at point (r,z)
c
c densrpsi(r,psi): density as a function of r and potential psi.
c
c pot(r,z): potential at point (r,z)
c
c all the rest is (slightly nobbled in places) numerical recipes routines.
c

      subroutine readharmfile(filename,ibuf1)
      parameter(pi=3.1415926535)
      parameter(nrsize=20000)
      common /potconstants/ apot(20,0:nrsize), frad(20,0:nrsize),
     +     dr, nr, lmax, potcor
      common /legendre/ plcon(0:40)
      common /gparameters/  a, b, c, v0, q, psi00,
     +                      psiout, rho1, sigbulge2,
     +                      rmdisk, rdisk, zdisk, outdisk, drtrunc,
     +                      potcor1
      common /moreconstants/ v02, v03, rdisk2, diskconst, bulgea
      common /diskpars/ sigr0, disksr, nrdisk
      common /flags/ idiskflag, ibulgeflag, ihaloflag

      real adens(20,0:nrsize)
      integer*4 ibuf(18),ibufb(3), ibuf1(18)
c this allows the gparameters to be passes as a C structure
      character filename*60
      equivalence (a, ibuf),(idiskflag,ibufb)
      external alldensnohalo,alldensnobulge,zerofunc

c  constants for legendre polynomials
      do l=0,40
         plcon(l) = sqrt((2*l+1)/(4.0*pi))
      enddo
c
      open(14,file=filename,status='old')
      read(14,*)
      read(14,'(2x,7g15.5,i6,i4)') psi00,v0,q,a,b,c,dr,nr,lmax
      read(14,*)
      read(14,'(2x,3g15.5)') psiout,rho1,sigbulge
      read(14,*)
      read(14,'(2x,5g15.5)') rmdisk, rdisk, zdisk, outdisk, drtrunc
      read(14,'(2x,3i5)') idiskflag, ibulgeflag, ihaloflag
      read(14,*)
!      write(0,*) ' this model had psi00,v0,q,a,b,c,dr,nr,lmax='
!      write(0,'(x,7g15.5,i6,i4)') psi00,v0,q,a,b,c,dr,nr,lmax

c these are needed for the bulgeconstants common block
c
      if(ibulgeflag.NE.3) then
      sigbulge2=sigbulge*sigbulge
      bulgea = rho1*exp((psi00 - psiout)/sigbulge2)
      endif
c
c additional disk constants for efficiency - common diskconstants
c
      nrdisk= int((outdisk + 2.0*drtrunc)/dr) + 10
      rdisk2 = rdisk*rdisk
      diskconst = rmdisk/(4.0*pi*rdisk2*zdisk)
c
c and some more constants...
c
      if(ihaloflag.NE.3) then
      v02=v0*v0
      v03=v0*v02
      endif
! in case ihaloflag=3: v0=rscl,q=rv,a=rhs ...
! in case ibulgeflag=3: psiout=bulgerscl, rho1=bulgerv,sigbulge=bulgerhs....

      if(ihaloflag.eq.3) call sethalodf(v0,q,a)
 !     print*,'##',v0,q,a,ihaloflag
      if(ibulgeflag.eq.3)
     &      call setbulgedf(psiout,rho1,sigbulge,halormin(),halormax())
! ugly
!      print*,psiout,rho1,sigbulge,halormax()
!      print*,v0,q,a

      if(ihaloflag.eq.3) then
        call inithalodf(alldensnohalo)
      endif

      if(ibulgeflag.eq.3) then
        call initbulgedf(alldensnobulge)
      endif

      write(0,*) 'Reading harmonics file '//trim(filename)
      do ir=0,nr
         read(14,'(8g16.8)') rr,(adens(l/2+1,ir),l=0,lmax,2)
         enddo
      read(14,*)
      read(14,*)
      do ir=0,nr
         read(14,'(8g16.8)') rr,(apot(l/2+1,ir),l=0,lmax,2)
         enddo
      read(14,*)
      read(14,*)
      do ir=0,nr
         read(14,'(8g16.8)') rr,(frad(l/2+1,ir),l=0,lmax,2)
         enddo
      close(14)
c      write(0,*) 'have read density and potential'
c      write(*,*) 'lmax =',lmax
c      write(*,*) 'apot starts as:'
c      do ll=1,3
c         write(*,*) (apot(ll,i),i=1,10)
c      enddo
c
c
c
c calculate potential correction
c
!!!!
      redge = nr*dr
      potcor = 0
      do l=0,lmax,2
          potcor = potcor + apot(l/2+1,nr) + frad(l/2+1,nr)*redge/(l+1)
      enddo
      potcor = potcor*plcon(0)
      potcor1 = potcor
!!!
c
c transfer gparams buffer for output
c

      do i=1,15
          ibuf1(i) = ibuf(i)
      enddo
      do i=1,3
          ibuf1(i+15) = ibufb(i)
      enddo

      return
      end

max_stars_repo_name: rknop/amuse

max_issues_repo_name: rknop/amuse

max_forks_repo_name: rieder/amuse


text:
C Test case for PR fortran/3743
C Origin: David Billinghurst <David.Billinghurst@riotinto.com>
C
C { dg-do link }
      integer   i
      i = bit_size(i)
      end


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2015 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine e_c3d_em(co,konl,lakonl,s,sm,
     &  ff,nelem,nmethod,ielmat,
     &  ntmat_,t1,ithermal,vold,idist,
     &  matname,mi,mass,rhsi,
     &  ncmat_,elcon,nelcon,h0,iactive,
     &  alcon,nalcon,istartset,iendset,ialset)
!
!     computation of the element matrix and rhs for the element with
!     the topology in konl
!
!     ff: rhs without temperature and eigenstress contribution
!
!     nmethod=0: check for positive Jacobian
!     nmethod=1: stiffness matrix + right hand side
!     nmethod=2: stiffness matrix + mass matrix
!     nmethod=3: static stiffness + buckling stiffness
!     nmethod=4: stiffness matrix + mass matrix
!
      implicit none
!
      logical mass,rhsi
!
      character*8 lakonl
      character*80 matname(*),amat
!
      integer konl(26),ifaceq(8,6),nelem,nmethod,iactive(3),
     &  ithermal,idist,i,j,k,i1,m,one,ii,jj,id,ipointer,ig,kk,mi(*),
     &  ielmat(mi(3),*),ntmat_,nope,nopes,imat,mint2d,mint3d,
     &  ifacet(6,4),nopev,ifacew(8,5),ipointeri,ipointerj,iflag,
     &  nelcon(2,*),ncmat_,nalcon(2,*),iel,ii1,ilength,istart,iset,
     &  isurf,jj1,istartset(*),iendset(*),ialset(*),three,nfaces
!
      real*8 co(3,*),xl(3,26),shp(4,26),s(100,100),ff(100),xs2(3,7),
     &  t1(*),h0(3,*),xl2(3,9),xsj2(3),shp2(7,9),vold(0:mi(2),*),
     &  xi,et,ze,xsj,sm(100,100),t1l,weight,
     &  elcon(0:ncmat_,ntmat_,*),elconloc(21),sigma,
     &  h0l(3),h0l2(3,8),alcon(0:6,ntmat_,*),diag,um,uminv,alpha(6),
     &  d(3,3)
!
      include "gauss.f"
!
      data ifaceq /4,3,2,1,11,10,9,12,
     &            5,6,7,8,13,14,15,16,
     &            1,2,6,5,9,18,13,17,
     &            2,3,7,6,10,19,14,18,
     &            3,4,8,7,11,20,15,19,
     &            4,1,5,8,12,17,16,20/
      data ifacet /1,3,2,7,6,5,
     &             1,2,4,5,9,8,
     &             2,3,4,6,10,9,
     &             1,4,3,8,10,7/
      data ifacew /1,3,2,9,8,7,0,0,
     &             4,5,6,10,11,12,0,0,
     &             1,2,5,4,7,14,10,13,
     &             2,3,6,5,8,15,11,14,
     &             4,6,3,1,12,15,9,13/
      data d /1.d0,0.d0,0.d0,0.d0,1.d0,0.d0,0.d0,0.d0,1.d0/
!
      one=1
      three=3
      iflag=3
!
      imat=ielmat(1,nelem)
      amat=matname(imat)
!
      if(lakonl(4:5).eq.'20') then
         nope=20
         nopev=8
         nopes=8
      elseif(lakonl(4:4).eq.'8') then
         nope=8
         nopev=8
         nopes=4
      elseif(lakonl(4:5).eq.'10') then
         nope=10
         nopev=4
         nopes=6
      elseif(lakonl(4:4).eq.'4') then
         nope=4
         nopev=4
         nopes=3
      elseif(lakonl(4:5).eq.'15') then
         nope=15
         nopev=6
      elseif(lakonl(4:4).eq.'6') then
         nope=6
         nopev=6
      endif
!
      if(lakonl(4:5).eq.'8R') then
         mint2d=1
         mint3d=1
      elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R')) then
         mint2d=4
         mint3d=8
      elseif(lakonl(4:4).eq.'2') then
         mint2d=9
         mint3d=27
      elseif(lakonl(4:5).eq.'10') then
         mint2d=3
         mint3d=4
      elseif(lakonl(4:4).eq.'4') then
         mint2d=1
         mint3d=1
      elseif(lakonl(4:5).eq.'15') then
         mint3d=9
      elseif(lakonl(4:4).eq.'6') then
         mint3d=2
      else
         mint3d=0
      endif
!
!     computation of the coordinates of the local nodes
!
      do i=1,nope
        do j=1,3
          xl(j,i)=co(j,konl(i))
        enddo
      enddo
!
!       initialisation for distributed forces
!
      if(rhsi) then
        if(idist.ne.0) then
          do i=1,5*nope
            ff(i)=0.d0
          enddo
        endif
      endif
!
!     initialisation of sm
!
      if(mass) then
        do i=1,5*nope
          do j=i,5*nope
            sm(i,j)=0.d0
          enddo
        enddo
      endif
!
!     initialisation of s
!
      do i=1,5*nope
        do j=i,5*nope
          s(i,j)=0.d0
        enddo
      enddo
!
!     computation of the matrix: loop over the Gauss points
!
      do kk=1,mint3d
         if(lakonl(4:5).eq.'8R') then
            xi=gauss3d1(1,kk)
            et=gauss3d1(2,kk)
            ze=gauss3d1(3,kk)
            weight=weight3d1(kk)
         elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R'))
     &           then
            xi=gauss3d2(1,kk)
            et=gauss3d2(2,kk)
            ze=gauss3d2(3,kk)
            weight=weight3d2(kk)
         elseif(lakonl(4:4).eq.'2') then
            xi=gauss3d3(1,kk)
            et=gauss3d3(2,kk)
            ze=gauss3d3(3,kk)
            weight=weight3d3(kk)
         elseif(lakonl(4:5).eq.'10') then
            xi=gauss3d5(1,kk)
            et=gauss3d5(2,kk)
            ze=gauss3d5(3,kk)
            weight=weight3d5(kk)
         elseif(lakonl(4:4).eq.'4') then
            xi=gauss3d4(1,kk)
            et=gauss3d4(2,kk)
            ze=gauss3d4(3,kk)
            weight=weight3d4(kk)
         elseif(lakonl(4:5).eq.'15') then
            xi=gauss3d8(1,kk)
            et=gauss3d8(2,kk)
            ze=gauss3d8(3,kk)
            weight=weight3d8(kk)
         else
            xi=gauss3d7(1,kk)
            et=gauss3d7(2,kk)
            ze=gauss3d7(3,kk)
            weight=weight3d7(kk)
         endif
!
!     calculation of the shape functions and their derivatives
!     in the gauss point
!
         if(nope.eq.20) then
            call shape20h(xi,et,ze,xl,xsj,shp,iflag)
         elseif(nope.eq.8) then
            call shape8h(xi,et,ze,xl,xsj,shp,iflag)
         elseif(nope.eq.10) then
            call shape10tet(xi,et,ze,xl,xsj,shp,iflag)
         elseif(nope.eq.4) then
            call shape4tet(xi,et,ze,xl,xsj,shp,iflag)
         elseif(nope.eq.15) then
            call shape15w(xi,et,ze,xl,xsj,shp,iflag)
         else
            call shape6w(xi,et,ze,xl,xsj,shp,iflag)
         endif
!
!           check the jacobian determinant
!
         if(xsj.lt.1.d-20) then
            write(*,*) '*ERROR in e_c3d_th: nonpositive jacobian'
            write(*,*) '       determinant in element',nelem
            write(*,*)
            xsj=dabs(xsj)
            nmethod=0
         endif
!
!        calculating the temperature and the magnetic intensity
!        in the integration point (one order lower than the
!        interpolation of the primary unknowns)
!
         do j=1,3
            h0l(j)=0.d0
         enddo
         t1l=0.d0
!
!        calculating the magnetic intensity
!
         if(lakonl(4:5).eq.'8 ') then
            do i1=1,nope
               do j=1,3
                  h0l(j)=h0l(j)+h0(j,konl(i1))/8.d0
               enddo
            enddo
         elseif(lakonl(4:6).eq.'20 ')then
            call linvec(h0,konl,nope,kk,h0l,one,three)
         else
            do i1=1,nope
               do j=1,3
                  h0l(j)=h0l(j)+shp(4,i1)*h0(j,konl(i1))
               enddo
            enddo
         endif
!
!        calculating the temperature
!
         if(ithermal.eq.1) then
            if(lakonl(4:5).eq.'8 ') then
               do i1=1,nope
                  t1l=t1l+t1(konl(i1))/8.d0
               enddo
            elseif(lakonl(4:6).eq.'20 ')then
               call linscal(t1,konl,nope,kk,t1l,one)
            else
               do i1=1,nope
                  t1l=t1l+shp(4,i1)*t1(konl(i1))
               enddo
            endif
         elseif(ithermal.ge.2) then
            if(lakonl(4:5).eq.'8 ') then
               do i1=1,nope
                  t1l=t1l+vold(0,konl(i1))/8.d0
               enddo
            elseif(lakonl(4:6).eq.'20 ')then
               call linscal(vold,konl,nope,kk,t1l,mi(2))
            else
               do i1=1,nope
                  t1l=t1l+shp(4,i1)*vold(0,konl(i1))
               enddo
            endif
         endif
!
!        material data (electric conductivity and
!        magnetic permeability)
!
         call materialdata_em(elcon,nelcon,alcon,nalcon,
     &        imat,ntmat_,t1l,elconloc,ncmat_,alpha)
!
         weight=weight*xsj
!
         uminv=weight/elconloc(1)
         um=weight*elconloc(1)
         sigma=weight*alpha(1)
!
         jj1=0
         do jj=1,nope
            ii1=0
            do ii=1,jj
!
               if((int(elconloc(2)).eq.2).or.(int(elconloc(2)).eq.3))
     &             then
!
!                 K_AA matrix
!
                  diag=shp(1,ii)*shp(1,jj)+shp(2,ii)*shp(2,jj)+
     &                 shp(3,ii)*shp(3,jj)
                  do k=1,3
                     do m=1,3
                        s(ii1+k,jj1+m)=s(ii1+k,jj1+m)+
     &                    (diag*d(k,m)-shp(m,ii)*shp(k,jj)
     &                                +shp(k,ii)*shp(m,jj))*uminv
                     enddo
                  enddo
!
!                 M_AA matrix
!
                  if(mass) then
                     do k=1,3
                        sm(ii1+k,jj1+k)=sm(ii1+k,jj1+k)+
     &                       sigma*shp(4,ii)*shp(4,jj)
                     enddo
                  endif
                  if((int(elconloc(2)).eq.2).and.mass) then
!
!                    M_AV and M_VA matrix
!
                     do k=1,3
                        sm(ii1+k,jj1+4)=sm(ii1+k,jj1+4)+
     &                       sigma*shp(4,ii)*shp(k,jj)
                        sm(ii1+4,jj1+k)=sm(ii1+4,jj1+k)+
     &                       sigma*shp(k,ii)*shp(4,jj)
                     enddo
!
!                    M_VV matrix
!
                     sm(ii1+4,jj1+4)=sm(ii1+4,jj1+4)+
     &                    sigma*(shp(1,ii)*shp(1,jj)+
     &                    shp(2,ii)*shp(2,jj)+
     &                    shp(3,ii)*shp(3,jj))
                  endif
               else if(int(elconloc(2)).eq.1) then
!
!                 K_phiphi matrix
!
                  s(ii1+5,jj1+5)=s(ii1+5,jj1+5)-
     &                     um*(shp(1,ii)*shp(1,jj)+
     &                    shp(2,ii)*shp(2,jj)+
     &                    shp(3,ii)*shp(3,jj))
               endif
!
               ii1=ii1+5
            enddo
!
!           F_phi matrix
!
            if((rhsi).and.(int(elconloc(2)).eq.1)) then
               ff(jj1+5)=ff(jj1+5)-um*(shp(1,jj)*h0l(1)+
     &                                 shp(2,jj)*h0l(2)+
     &                                 shp(3,jj)*h0l(3))
            endif
!
            jj1=jj1+5
         enddo
      enddo
!
!     surface integrals
!
!     determining the number of faces per element
!
      if((lakonl(4:4).eq.'8').or.(lakonl(4:4).eq.'2')) then
         nfaces=6
      elseif((lakonl(4:4).eq.'6').or.(lakonl(4:5).eq.'15')) then
         nfaces=5
      elseif((lakonl(4:4).eq.'4').or.(lakonl(4:5).eq.'10')) then
         nfaces=4
      endif
!
      m=int(elconloc(2))
      if(iactive(m).eq.0) return
         iset=iactive(m)
         istart=istartset(iset)
         ilength=iendset(iset)-istart+1
!
         isurf=10*nelem+nfaces
         call nident(ialset(istart),isurf,ilength,id)
!
         do
            if(id.eq.0) exit
            isurf=ialset(istart+id-1)
            iel=int(isurf/10.d0)
            if(iel.ne.nelem) exit
            ig=isurf-10*iel
!
!     treatment of wedge faces
!
            if(lakonl(4:4).eq.'6') then
               mint2d=1
               if(ig.le.2) then
                  nopes=3
               else
                  nopes=4
               endif
            endif
            if(lakonl(4:5).eq.'15') then
               if(ig.le.2) then
                  mint2d=3
                  nopes=6
               else
                  mint2d=4
                  nopes=8
               endif
            endif
!
            if((nope.eq.20).or.(nope.eq.8)) then
               do i=1,nopes
                  do j=1,3
                     h0l2(j,i)=h0(j,konl(ifaceq(i,ig)))
                     xl2(j,i)=co(j,konl(ifaceq(i,ig)))
                  enddo
               enddo
            elseif((nope.eq.10).or.(nope.eq.4)) then
               do i=1,nopes
                  do j=1,3
                     h0l2(j,i)=h0(j,konl(ifacet(i,ig)))
                     xl2(j,i)=co(j,konl(ifacet(i,ig)))
                  enddo
               enddo
            else
               do i=1,nopes
                  do j=1,3
                     h0l2(j,i)=h0(j,konl(ifacew(i,ig)))
                     xl2(j,i)=co(j,konl(ifacew(i,ig)))
                  enddo
               enddo
            endif
!
            do i=1,mint2d
!
               if((lakonl(4:5).eq.'8R').or.
     &              ((lakonl(4:4).eq.'6').and.(nopes.eq.4))) then
                  xi=gauss2d1(1,i)
                  et=gauss2d1(2,i)
                  weight=weight2d1(i)
               elseif((lakonl(4:4).eq.'8').or.(lakonl(4:6).eq.'20R').or.
     &                 ((lakonl(4:5).eq.'15').and.(nopes.eq.8))) then
                  xi=gauss2d2(1,i)
                  et=gauss2d2(2,i)
                  weight=weight2d2(i)
               elseif(lakonl(4:4).eq.'2') then
                  xi=gauss2d3(1,i)
                  et=gauss2d3(2,i)
                  weight=weight2d3(i)
               elseif((lakonl(4:5).eq.'10').or.
     &                 ((lakonl(4:5).eq.'15').and.(nopes.eq.6))) then
                  xi=gauss2d5(1,i)
                  et=gauss2d5(2,i)
                  weight=weight2d5(i)
               elseif((lakonl(4:4).eq.'4').or.
     &                 ((lakonl(4:4).eq.'6').and.(nopes.eq.3))) then
                  xi=gauss2d4(1,i)
                  et=gauss2d4(2,i)
                  weight=weight2d4(i)
               endif
!
               if(nopes.eq.8) then
                  call shape8q(xi,et,xl2,xsj2,xs2,shp2,iflag)
               elseif(nopes.eq.4) then
                  call shape4q(xi,et,xl2,xsj2,xs2,shp2,iflag)
               elseif(nopes.eq.6) then
                  call shape6tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
               else
                  call shape3tri(xi,et,xl2,xsj2,xs2,shp2,iflag)
               endif
!
               do k=1,3
                  h0l(k)=0.d0
                  do j=1,nopes
                     h0l(k)=h0l(k)+h0l2(k,j)*shp2(4,j)
                  enddo
               enddo
!
               if((rhsi).and.(m.gt.1)) then
                  do k=1,nopes
                     if((nope.eq.20).or.(nope.eq.8)) then
                        ipointer=5*(ifaceq(k,ig)-1)
                     elseif((nope.eq.10).or.(nope.eq.4)) then
                        ipointer=5*(ifacet(k,ig)-1)
                     else
                        ipointer=5*(ifacew(k,ig)-1)
                     endif
!
!                    F_A vector
!
                     ff(ipointer+1)=ff(ipointer+1)+
     &                  shp2(4,k)*(h0l(2)*xsj2(3)-h0l(3)*xsj2(2))*weight
                     ff(ipointer+2)=ff(ipointer+2)+
     &                  shp2(4,k)*(h0l(3)*xsj2(1)-h0l(1)*xsj2(3))*weight
                     ff(ipointer+3)=ff(ipointer+3)+
     &                  shp2(4,k)*(h0l(1)*xsj2(2)-h0l(2)*xsj2(1))*weight
!
                  enddo
               endif
!
               do ii=1,nopes
                  if((nope.eq.20).or.(nope.eq.8)) then
                     ipointeri=5*(ifaceq(ii,ig)-1)
                  elseif((nope.eq.10).or.(nope.eq.4)) then
                     ipointeri=5*(ifacet(ii,ig)-1)
                  else
                     ipointeri=5*(ifacew(ii,ig)-1)
                  endif
                  do jj=1,nopes
                     if((nope.eq.20).or.(nope.eq.8)) then
                        ipointerj=5*(ifaceq(jj,ig)-1)
                     elseif((nope.eq.10).or.(nope.eq.4)) then
                        ipointerj=5*(ifacet(jj,ig)-1)
                     else
                        ipointerj=5*(ifacew(jj,ig)-1)
                     endif
!
                     if(m.gt.1) then
!
!                       K_Aphi matrix
!
                        if((ipointeri+1).gt.ipointerj+5) cycle
                        s(ipointeri+1,ipointerj+5)=
     &                       s(ipointeri+1,ipointerj+5)
     &                       +shp2(4,jj)*(shp2(3,ii)*xsj2(2)
     &                                   -shp2(2,ii)*xsj2(3))
     &                       *weight
!
                        if((ipointeri+2).gt.ipointerj+5) cycle
                        s(ipointeri+2,ipointerj+5)=
     &                       s(ipointeri+2,ipointerj+5)
     &                       +shp2(4,jj)*(shp2(1,ii)*xsj2(3)
     &                                   -shp2(3,ii)*xsj2(1))
     &                       *weight
!
                        if((ipointeri+3).gt.ipointerj+5) cycle
                        s(ipointeri+3,ipointerj+5)=
     &                       s(ipointeri+3,ipointerj+5)
     &                       +shp2(4,jj)*(shp2(2,ii)*xsj2(1)
     &                                   -shp2(1,ii)*xsj2(2))
     &                       *weight
!
!                       K_phiA matrix
!
                        if(ipointeri+5.gt.(ipointerj+1)) cycle
                        s(ipointeri+5,ipointerj+1)=
     &                       s(ipointeri+5,ipointerj+1)
     &                       +shp2(4,ii)*(shp2(3,jj)*xsj2(2)
     &                                   -shp2(2,jj)*xsj2(3))
     &                       *weight
!
                        if(ipointeri+5.gt.(ipointerj+2)) cycle
                        s(ipointeri+5,ipointerj+2)=
     &                       s(ipointeri+5,ipointerj+2)
     &                       +shp2(4,ii)*(shp2(1,jj)*xsj2(3)
     &                                   -shp2(3,jj)*xsj2(1))
     &                       *weight
!
                        if(ipointeri+5.gt.(ipointerj+3)) cycle
                        s(ipointeri+5,ipointerj+3)=
     &                       s(ipointeri+5,ipointerj+3)
     &                       +shp2(4,ii)*(shp2(2,jj)*xsj2(1)
     &                                   -shp2(1,jj)*xsj2(2))
     &                       *weight
                     endif
                  enddo
               enddo
!
            enddo
!
            id=id-1
         enddo
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
module index_space_test
  use legion_fortran
  use iso_c_binding
  implicit none

  integer(c_int), parameter :: TOP_LEVEL_TASK_ID = 0
  integer(c_int), parameter :: INDEX_SPACE_TASK_ID = 1

contains
  subroutine index_space_task(tdata, tdatalen, userdata, userlen, p)
    implicit none

    type(c_ptr), value, intent(in) :: tdata
    integer(c_size_t), value, intent(in) :: tdatalen
    type(c_ptr), value, intent(in) ::userdata
    integer(c_size_t), value, intent(in) :: userlen
    integer(c_long_long), value, intent(in) :: p

    type(FRuntime) :: runtime
    type(FContext) :: ctx
    type(FTask) :: task
    type(FPhysicalRegionList) :: pr_list

    integer(c_size_t) :: arglen
    integer(kind=4), pointer :: task_arg
    type(c_ptr) :: task_arg_ptr
    integer(kind=4), target :: input

    call legion_task_prolog(tdata, tdatalen, userdata, userlen, p, &
                            task, pr_list, &
                            ctx, runtime)
    arglen = task%get_local_arglen()
    task_arg_ptr = task%get_local_args()
    call c_f_pointer(task_arg_ptr, task_arg)
    input = task_arg
    print *, "Index Task", input
    input = input * 2
    call legion_task_epilog(runtime, ctx, c_loc(input), c_sizeof(input))
  end subroutine index_space_task

  subroutine top_level_task(tdata, tdatalen, userdata, userlen, p)
    implicit none

    type(c_ptr), value, intent(in) :: tdata
    integer(c_size_t), value, intent(in) :: tdatalen
    type(c_ptr), value, intent(in) ::userdata
    integer(c_size_t), value, intent(in) :: userlen
    integer(c_long_long), value, intent(in) :: p

    type(FRuntime) :: runtime
    type(FContext) :: ctx
    type(FTask) :: task
    type(FPhysicalRegionList) :: pr_list
    type(FArgumentMap) :: arg_map
    type(FIndexLauncher) :: index_launcher
    type(FRect1D) :: launch_bounds
    type(FIndexSpace) :: launch_is
    type(FFutureMap) :: fm
    type(FFuture) :: f

    integer :: num_points = 4

    integer :: i = 0
    integer(kind=4), target :: input = 0
    integer(kind=4) :: received, expected
    logical :: all_passed = .true.

    Print *, "TOP Level Task!"

    call legion_task_prolog(tdata, tdatalen, userdata, userlen, p, &
                            task, pr_list, &
                            ctx, runtime)

    launch_bounds = FRect1D(0, num_points-1)
    launch_is = runtime%create_index_space(ctx, launch_bounds)

    arg_map = FArgumentMap()
    do i = 0, num_points-1
      input = i + 10
      call arg_map%set_point(i, FTaskArgument(c_loc(input), c_sizeof(input)))
    end do

    index_launcher = FIndexLauncher(INDEX_SPACE_TASK_ID, launch_is, &
                                    FTaskArgument(), arg_map)
    fm = runtime%execute_index_space(ctx, index_launcher)
    call fm%wait_all_results()
    do i = 0, num_points-1
      f = fm%get_future(i)
      call f%get_result(received)
      print *, "result", received
      expected = 2*(i+10)
      if (expected .ne. received) then
        print *, "Check failed for point", i, expected, received
      end if
    end do

    if (all_passed) then
      print *, "All checks passed!"
    end if

    call index_launcher%destroy()
    call arg_map%destroy()

    call legion_task_epilog(runtime, ctx)
  end subroutine top_level_task
end module index_space_test

Program index_space_main
  use iso_c_binding
  use legion_fortran
  use index_space_test
  implicit none

  type(FTaskVariantRegistrar) :: registrar_top, registrar_index_space

  call set_top_level_task_id(TOP_LEVEL_TASK_ID)

  registrar_top = FTaskVariantRegistrar(TOP_LEVEL_TASK_ID)
  call registrar_top%add_constraint(FProcessorConstraint(LOC_PROC))
  call preregister_task_variant(top_level_task, registrar_top, "top_level_task")
  call registrar_top%destroy()

  registrar_index_space = FTaskVariantRegistrar(INDEX_SPACE_TASK_ID)
  call registrar_index_space%add_constraint(FProcessorConstraint(LOC_PROC))
  call registrar_index_space%set_leaf()
  call preregister_task_variant(index_space_task, registrar_index_space, "index_space_task")
  call registrar_index_space%destroy()

  call legion_runtime_start()
End Program index_space_main

max_stars_repo_name: karasevb/legion

max_issues_repo_name: karasevb/legion

max_forks_repo_name: karasevb/legion


text:
      SUBROUTINE balanc(a,n,np)
      INTEGER n,np
      REAL a(np,np),RADIX,SQRDX
      PARAMETER (RADIX=2.,SQRDX=RADIX**2)
      INTEGER i,j,last
      REAL c,f,g,r,s
1     continue
        last=1
        do 14 i=1,n
          c=0.
          r=0.
          do 11 j=1,n
            if(j.ne.i)then
              c=c+abs(a(j,i))
              r=r+abs(a(i,j))
            endif
11        continue
          if(c.ne.0..and.r.ne.0.)then
            g=r/RADIX
            f=1.
            s=c+r
2           if(c.lt.g)then
              f=f*RADIX
              c=c*SQRDX
            goto 2
            endif
            g=r*RADIX
3           if(c.gt.g)then
              f=f/RADIX
              c=c/SQRDX
            goto 3
            endif
            if((c+r)/f.lt.0.95*s)then
              last=0
              g=1./f
              do 12 j=1,n
                a(i,j)=a(i,j)*g
12            continue
              do 13 j=1,n
                a(j,i)=a(j,i)*f
13            continue
            endif
          endif
14      continue
      if(last.eq.0)goto 1
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/balanc.for


text:
	FUNCTION PR_WVSW  ( dosw, posw, hosw )
C************************************************************************
C* PR_WVSW								*
C*									*
C* This function computes WVSW, the combined swell wave direction in    *
C* tens of degrees, the swell wave period in seconds, and the swell     *
C* wave height in half meters, from DOSW, POSW, and HOSW.               *
C* The following equation is used:                                      *
C*									*
C*       WVSW = ( DOSW / 10 ) * 10000 + PR_WVPH ( POSW, 2 * HOSW, 0 )   *
C*									*
C* REAL PR_WVSW  ( DOSW, POSW, HOSW )                                   *
C*									*
C* Input parameters:							*
C*	DOSW		REAL	Swell wave direction in degrees         *
C*	POSW		REAL	Swell wave period in seconds            *
C*	HOSW 		REAL	Swell wave height in meters             *
C*									*
C* Output parameters:							*
C*	PR_WVSW		REAL	Combined direction, period and height   *
C**									*
C* Log:									*
C* D. Kidwell/NCEP	 2/99	                                        *
C************************************************************************
        INCLUDE   'GEMPRM.PRM'
        INCLUDE   'ERMISS.FNC'
C------------------------------------------------------------------------
	PR_WVSW = RMISSD
C
	IF ( ( .not. ERMISS ( dosw ) ) .and.
     +       ( .not. ERMISS ( posw ) ) .and.
     +       ( .not. ERMISS ( hosw ) ) ) THEN
	    pphh = PR_WVPH ( posw, 2. * hosw, 0. )
	    IF ( .not. ERMISS ( pphh ) ) THEN
	        PR_WVSW = NINT ( dosw / 10. ) * 10000 + NINT ( pphh )
	    END IF
	END IF
C*
	RETURN
	END


text:
C
C**************************  TPNODI2D **********************************
C
C  ripartisco ai nodi le quantita' che sono date per
C  elementi e che devono andare ai nodi
C
C***********************************************************************
C
      SUBROUTINE TPNODI2d(nnod,ntri,TRIANG,z,TP2d,eltria)

      IMPLICIT NONE
      include 'CATHY.H'
      INTEGER  K,II,INOD,ISTOP
      INTEGER  nnod,ntri
      INTEGER  TRIANG(4,*),TP2d(*)
      REAL*8   z(*),eltria(*)
!CMS
!CMS      REAL*8   ZTEMP(2856)
!CMS      INTEGER  jj,kk,kz
!CMS

      do k=1,nnod
        tp2d(k)=0
        z(k)=0.0d0
      end do

      DO K=1,NTRI
         DO II=1,3
            INOD=TRIANG(II,K)
            z(INOD)=z(INOD) + eltria(K)
            TP2d(INOD)=TP2d(INOD)+1
         END DO
      END DO
      ISTOP=0
      DO K=1,NNOD
         IF (TP2d(K) .EQ. 0) THEN
            WRITE(6,1000) K
            ISTOP=1
         ELSE
            z(K)=z(K)/TP2d(K)
         END IF
      END DO
!CMS
!CMS      OPEN(98,file='./grid2d')
!CMS      DO kk=1,NNOD
!CMS          READ(98,*) ZTEMP(kk)
!CMS          z(kk) = ZTEMP(kk)
!CMS      END DO
!CMS      CLOSE(98)
!CMS
      IF (ISTOP .NE. 0) THEN
         STOP
      END IF
C
      RETURN
 1000 FORMAT(/,' ERROR IN MESH CONSTRUCTION:    NODE ',I7,
     1         ' IS NOT CONNECTED TO ANY ELEMENTS')
      END


!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    get_wrf_mass_ensperts  read wrf arw model ensemble members
!   prgmmr: mtong           org: np22                date: 2010-06-28
!
! abstract: dummy for read ensemble members from the wrf arw model in both
! binary and netcdf
!             format, for use with hybrid ensemble option. ensemble spread is
!             also
!             written out as a byproduct for diagnostic purposes.
!
!
! program history log:
!
!   input argument list:
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$ end documentation block

module get_wrf_mass_ensperts_mod
use abstract_get_wrf_mass_ensperts_mod
  type, extends(abstract_get_wrf_mass_ensperts_class) :: get_wrf_mass_ensperts_class
  contains
    procedure, pass(this) :: get_wrf_mass_ensperts => get_wrf_mass_ensperts_dummy
    procedure, pass(this) :: ens_spread_dualres_regional => ens_spread_dualres_regional_dummy
!   procedure, pass(this) :: get_wrf_mass_ensperts_netcdf_dummy
  end type get_wrf_mass_ensperts_class

contains
  subroutine ens_spread_dualres_regional_dummy(this,mype,en_perts,nelen,en_bar)
    use kinds, only: i_kind
    use gsi_bundlemod, only: gsi_bundle
    implicit none
    class(get_wrf_mass_ensperts_class), intent(inout) :: this
    integer(i_kind),intent(in):: mype
    type(gsi_bundle),allocatable, intent(in   ) :: en_perts(:,:)
    integer(i_kind), intent(in   ):: nelen
    type(gsi_bundle),OPTIONAL,intent(in):: en_bar
    write(6,*)'get_wrf_mass_ensperts:  ***WARNING*** dummy call ... does nothing!'
  return
  end subroutine ens_spread_dualres_regional_dummy

  subroutine get_wrf_mass_ensperts_dummy(this,en_perts,nelen,ps_bar)
    use kinds, only: i_kind,r_single
    use gsi_bundlemod, only: gsi_bundle
    implicit none
    class(get_wrf_mass_ensperts_class), intent(inout) :: this
    type(gsi_bundle),allocatable, intent(inout) :: en_perts(:,:)
    integer(i_kind), intent(in   ):: nelen
    real(r_single),dimension(:,:,:),allocatable:: ps_bar
    write(6,*)'get_wrf_mass_ensperts:  ***WARNING*** dummy call ... does nothing!'
  return
  end subroutine get_wrf_mass_ensperts_dummy
end module get_wrf_mass_ensperts_mod