text:
      subroutine wet_pond

!!    ~ ~ ~ PURPOSE ~ ~ ~
!!    run wet pond processes

!!    ~ ~ ~ INCOMING VARIABLES ~ ~ ~ ~ ~ ~ ~
!!    name             |units        |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    wtp_evrsv        |none         |detention pond evaporation coefficient
!!    hhvaroute(2,:,:) |m^3 H2O      |water
!!    hhvaroute(3,:,:) |metric tons  |sediment or suspended solid load
!!    i_mo             |none         |current month of simulation
!!    sub_subp_dt(:,:)  |mm H2O      |precipitation for time step in subbasin
!!    wtp_pvol(:       |m^3 H2O      |volume of permanent pool including forebay
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ OUTGOING VARIABLES ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    name           |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    hhvaroute(2,:,:) |m^3 H2O    |water
!!    hhvaroute(3,:,:) |metric tons|sediment or suspended solid load
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ LOCAL DEFINITIONS ~ ~ ~
!!    name        |units         |definition
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~
!!    ii          |none          |time step counter
!!    k           |none          |weir stage counter
!!    titldum     |NA            |dummy string
!!    ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~ ~

!!    ~ ~ ~ SUBROUTINES/FUNCTIONS CALLED ~ ~ ~
!!    SWAT:  surf_seep,surf_evap,water_depth
!!    Intrinsic: Log,exp

!!    ~ ~ ~ ~ ~ ~ END SPECIFICATIONS ~ ~ ~ ~ ~ ~

      use parm
      implicit none

      character (len=80) :: titldum
      integer :: ii, k, sb
      real :: qin,qout,qpnd,sedin,sedout,sedpnd,spndconc
      real :: rf,imc,pndwdth,seep,evap,rain,seepa,surfa
      real :: a1,b1,qdepth,hdep_ext,usettle,tmpw,mu,qhyd,dp
      real :: alpha,beta,vol,hp,tvol,ht,decayexp,qweir,sedweir,coeff_R

      sb = inum1
      qout=0.; sedout=0.; qdepth=0.
      dp = wtp_dp(sb) * 0.1 !cm

      if (iyr<wtp_iyr(sb).or.(iyr==wtp_iyr(sb).and.i_mo<wtp_imo(sb)))
     &then
         return
      endif

!!    Estimate Pond Volume using Austin Design Manual if not entered by user
!!    COA Environmental Criteria Manual Chapter 1.6.6 and Table 1-9
      if (wtp_dim(sb)==0) then
         imc = subdr_ickm(ihout) / subdr_km(ihout) !fraction impervious cover
         rf = 0.5463 * imc ** 2 + 0.328 * imc + 0.0296
         wtp_pvol(sb) = 0.162 * rf * (subdr_km(ihout) * 247.11) !ac-feet
         wtp_pvol(sb) = wtp_pvol(sb) * 1233.4 !m3
         if (wtp_pvol(sb)<3000) wtp_pvol(sb) = 3000. !minimum area 21780ft2 * depth 5ft gives ~3000m3
         wtp_pdepth(sb) = 2. !m
         wtp_sdslope(sb) = 4. !h:v
         wtp_lenwdth(sb) = 2. !l:w
         wtp_pdia(sb) = 0.1524 !m (=6inch )
         wtp_plen(sb) = 5. !m
         wtp_pmann(sb) = 0.012 ! concrete surface
         wtp_ploss(sb) = 0.1 !minor loss

         call ext_dpth(wtp_extdepth(sb))
      end if

      alpha = wtp_lenwdth(sb)
      beta = wtp_sdslope(sb)
      hp = wtp_pdepth(sb)
      vol = wtp_pvol(sb)

      !pond width at the bottom
      a1 = 3. * beta * hp * (alpha + 1.)
      b1 = 12. * alpha * (3.* vol / hp - 4.* beta**2 * hp**2)
      pndwdth = (- a1 + sqrt(a1**2 + b1)) / (6.* alpha) !m

      !surface area of the permanent pool
      surfa = wtp_pvol(sb) / wtp_pdepth(sb) !m2,

      !pond volume (permanent pool + extended detention
      ht = hp + wtp_extdepth(sb)
      a1 = 2.* beta * ht
      tvol = (alpha * pndwdth **2 + (pndwdth + a1) * (alpha *
     & pndwdth + a1) + (pndwdth * (alpha * pndwdth + a1) + alpha *
     & pndwdth * (pndwdth + a1)) / 2.) * ht / 3.

!!    Get initial values from previous day
      qpnd = wtp_qi(sb) !m^3
      sedpnd = wtp_sedi(sb) * qpnd / 1.e6  !tons

!!	   iterate for subdaily flow/sediment routing
      do ii=1,nstep

         qout = 0.; spndconc = 0.
         qin = hhvaroute(2,ihout,ii)   !m^3
         qpnd =  qin + qpnd !m^3
         sedin = hhvaroute(3,ihout,ii) !tons

         !overflow to emergency outlet
         if (qpnd>tvol) then
            qweir = qpnd - tvol
            qpnd = tvol
            sedweir = sedin * qweir / qin
            sedpnd = sedpnd + sedin - sedweir
         else
            qweir = 0.
            sedweir = 0.
            sedpnd = sedin + sedpnd  !tons
         endif
         if (qpnd>1e-6) then
            spndconc = sedpnd / qpnd * 1.e6 !initial sed conc, mg/l
            if(spndconc<wtp_sede(sb)) spndconc = wtp_sede(sb)
         else
            cycle
         end if

!!       Estimate water depth, m
         call wpnd_depth(qpnd,pndwdth,beta,alpha,qdepth)
         hdep_ext = qdepth - hp
         if(hdep_ext<0) hdep_ext = 0.

!!       Calculate outflow
         if (hdep_ext>0.) then
!!          Use stage-discharge relationship if available
	      if (wtp_stagdis(sb)==1) then
	         select case(wtp_sdtype(sb))
	         case(1) !! 1 is exponential function
	            qout = wtp_sdc1(sb) * exp(wtp_sdexp(sb) * qdepth) +
     &   	            wtp_sdintc(sb)
	         case(2) !! 2 is Linear function
	            qout = wtp_sdc1(sb) * qdepth + wtp_sdintc(sb)
	         case(3) !! 3 is logarthmic function
	            qout = wtp_sdc1(sb) * log(qdepth) + wtp_sdintc(sb)
	         case(4) !! 4 is power function
	            qout = wtp_sdc1(sb) * (qdepth**3) + wtp_sdc2(sb) *
     &   	          (qdepth**2) + wtp_sdc3(sb) * qdepth + wtp_sdintc(sb)
               case(5)
                  qout = wtp_sdc1(sb) * (qdepth**wtp_sdexp(sb)) +
     &                wtp_sdintc(sb) !m3/s
	         end select
	      else
!!                Discharge out of extended detention storage through inverted PVC pipe
               call pipe_discharge(wtp_pdia(sb),wtp_plen(sb),
     &              hdep_ext,wtp_pmann(sb),wtp_ploss(sb),qout) !m3/s

	      end if
            qout = qout * idt * 60. !m3/s ->m^3
         !  no outflow from the permanent pool
		   if (qout>qpnd-wtp_pvol(sb)) qout = qpnd - wtp_pvol(sb)
             if (qout<0) qout = 0.
         else
!!          no discharge from the permanent pool
            qout = 0.
         end if

!!       Seepage, evaporation and rainfall
         a1 = wtp_lenwdth(sb)
         b1 = wtp_sdslope(sb) ** 2 + 1.
         seepa = a1 * pndwdth ** 2 + 2. * (a1 + 1) * pndwdth * qdepth *
     &      b1 ** 0.5 + 5.64 * qdepth ** 2 * b1 * wtp_sdslope(sb) !m2
         b1 = 2.* beta * qdepth                !jeong 05/22/14
         surfa = (a1 * pndwdth + b1) * (pndwdth + b1) !m2

         seep = wtp_k(sb) / 1000. / 60. * idt *  seepa !m3
         evap = wtp_evrsv(sb) * sub_pet(sb) / 1000. * surfa !m3
         rain = sub_subp_dt(sb,ii) / 1000. * surfa !m3

!!       Mass balance for flow
         qpnd = qpnd + rain - seep - evap - qout
         if (qpnd<0) qpnd = 0.

!----------------------------------------------------------------------------------
!!       Sediment removal
         tmpw = sub_hhwtmp(sb,ii)
         ! water viscosity (g/cm-s) using 3rd order polynomial interpolation
         mu = -3.e-6 * tmpw ** 3 + 0.0006 * tmpw ** 2 - 0.0469 *
     &        tmpw + 1.7517
         mu = mu * 1.e-2
         !settling velocity, cm/s
         usettle = 0.7 * 981. / 18. * 1.6 / mu * dp ** 2 !ro,s=1.6g/cm3, g=981cm/s
         b1 = 2.* beta * qdepth
         surfa = (a1 * pndwdth + b1) * (pndwdth + b1) !m2

         if(qout>0.) then
            !WERF equation
            qhyd = qout / idt / 60. / surfa * 100.  ! m3/s / m2 * 100 = cm/s
            coeff_R = (1.+ usettle /
     &        (qhyd * wtp_hydeff(sb))) ** (-wtp_hydeff(sb))
            spndconc = spndconc * coeff_R  ! mg/l
            sedpnd = spndconc * qpnd * 1e-6 !tons, amount sediment in the pond at the end of the time step
            !sediment conc no less than the minimum value
            if(spndconc<wtp_sede(sb)) spndconc = wtp_sede(sb)
         else
            !SWAT sediment equation for ponds
            decayexp = exp(-7.667e-3 * idt / 60. * dp)
            spndconc = (spndconc - wtp_sede(sb)) * decayexp +
     &         wtp_sede(sb)
            sedpnd = spndconc * qpnd * 1e-6 !tons
         endif

         !Sediment coming out of the pond
         sedout = spndconc * qout * 1e-6 !tons

         qout = qout + qweir
         sedout = sedout + sedweir

   	   !! Store flow/sediment out of the pond at the subbasin outlet
   	   hhvaroute(2,ihout,ii) = max(0.,qout)
   	   hhvaroute(3,ihout,ii) = max(0.,sedout)

      end do

      ! Store end-of-day values for next day
      wtp_qi(sb) = qpnd !m^3
      if (qpnd>0.1.and.sedpnd>0.000001) then
         wtp_sedi(sb) = sedpnd / qpnd * 1.e6 !tons
      else
         wtp_sedi(sb) = 0
      endif
	   if (wtp_sedi(sb)<wtp_sede(sb))  wtp_sedi(sb) = wtp_sede(sb)

      end subroutine

   !-------------------------------------------------------------------
      subroutine ext_dpth(hmax)
      use parm
      implicit none

      real,dimension(40) :: cumrain=(/0.,0.006,0.012,0.019,0.026,0.034,
     & 0.043,0.053,0.064,0.077,0.092,0.11,0.134,0.166,0.212,0.287,0.384,
     & 0.542,0.802,1.262,1.462,1.587,1.688,1.746,1.784,1.811,1.832,
     & 1.849,1.863,1.875,1.885,1.894,1.902,1.91,1.917,1.924,1.93,1.93,
     & 1.93,1.93/)
      real :: ia, ss, plen,inflow,outflow,pndvol,wdth,pndarea,vtmp,hdep
      real :: aa,bb,pdia,hvol,alpha
      real, dimension(40):: fa
      real, intent(out) :: hmax
      integer :: ii,sb

      sb = inum1
      fa=0; inflow=0.

      ss = 1000. / sub_cn2(sb) - 10.
      plen = wtp_plen(sb) * 3.2808 !ft
      pdia = wtp_pdia(sb) * 3.2808 !ft
      !volume of permanent pool
      pndvol = wtp_pvol(sb) * 3.2808 ** 3 !ft^3
      !pond width at the bottom of the pond
      wdth = ((140.3**2 - 4*(5828.2-0.072 * pndvol))**0.5 - 140.3) / 2. ! ft
      !pond area at the top of permanent pool
      pndarea = (wdth + 187.) * (2*wdth + 187.) ! ft^2
      alpha = wtp_lenwdth(sb)

      hmax=0.
      do ii=1,40
         if(hvol<0) hvol = 0.
         ! calculate surface runoff
         ia = 0.2 * ss
         if (ia>cumrain(ii)) then
            ia = cumrain(ii)
            fa(ii) = 0
         else
            fa(ii) = ss * (cumrain(ii) - ia) / (cumrain(ii) - ia + ss)
         endif
         if (ii>=2) inflow = fa(ii) - fa(ii-1) !inches

         ! calculate max water depth in extended detention
         inflow = inflow * subdr_km(sb) * 8.97e5 !inch-km2 ->ft^3
         hvol = hvol + inflow  !ft^3

         call wpnd_depth(hvol,wdth,wtp_sdslope(sb),alpha,hdep)
         if (hdep>hmax) hmax = hdep

         !compute discharge from outlet pipe
         call pipe_discharge(pdia,plen,hdep,wtp_pmann(sb),wtp_ploss(sb),
     &      outflow)

         outflow = outflow * 35.31 * 300.  !m3/s ->ft^3 per 5 minutes
         if (outflow>hvol) outflow = hvol
         hvol = hvol - outflow

      end do

      hmax = hmax / 3.2808 ! meter

      end subroutine
   !-------------------------------------------------------------------

      subroutine wpnd_depth(hvol,width,slp,lenwdth,hdep)
      !calculate ponding depth using Newton's method
      implicit none
      real, intent(in):: hvol,width,slp,lenwdth
      real, intent(out):: hdep
      real:: dfn,fn,alp,ll

      alp = lenwdth
      ll = width
      hdep = 0.1; fn = 100
      do while (abs(fn)>0.1)
         fn = alp * ll ** 2 * hdep + slp * ll * (1.+ alp) * hdep ** 2 +
     &     4. / 3.* slp ** 2 * hdep ** 3 - hvol
         dfn = alp * ll ** 2 + 2. * slp * ll * (1.+ alp) * hdep + 4.*
     &     slp ** 2 * hdep ** 2
         hdep = hdep - fn / dfn
      end do
      end subroutine

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

      subroutine pipe_discharge(pdia,plen,hdep,mann,mloss,outflow)
      ! calculate discharge from extended detention through pvc pipe,m3/s
      implicit none
      real, intent(in):: pdia,plen,hdep,mann,mloss
      real, intent(out):: outflow
      real:: parea,kf,rh

      parea = 3.14159 * (3.2808 * pdia) ** 2 / 4. !ft^2
      rh = 3.2808 * pdia / 4. !ft
      kf = 29. * plen * mann ** 2 / rh ** 1.33
      outflow = parea * ((64.4 * 3.2808 * hdep) / (1.+ kf + mloss))
     & ** 0.5 !cfs
      outflow = outflow / 35.31 !m3/s

      end subroutine


text:
!**********************************************************
!    RAMG tools, functions needed in ramg_driver and
!    ramg_interface.
!    Module defined in ramg_data.f
!
!    ramg_calcIAI
!    ramg_calcIvFC/CF
!    ramg_direct
!    ramg_allocate
!    ramg_deallocate
!    ramg_dump
!
!***********************************************************

!!***********************************************************
!      ramg_calcIvCF
!       calculate V coarse to fine,
!       v = I * VC
!***********************************************************
      subroutine ramg_calcIvCF(level1,level2,VC,vf,
     &           ilwork,BC,iBC,iper)
      use ramg_data
      include "common.h"
      include "mpif.h"

      !implicit none
      integer,intent(in)                :: level1,level2
      real(kind=8),intent(inout),dimension(amg_nshg(level1)) :: vf
      real(kind=8),intent(in),dimension(amg_nshg(level2)) :: VC

      integer,intent(in),dimension(nlwork)          :: ilwork
      integer,intent(in),dimension(nshg)            :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC

      integer                           :: i,j,k,p

      real(kind=8) :: cpusec(2)

      call cpu_time(cpusec(1))

      vf = 0
      do i=1,amg_nshg(level1)
         do k=I_cf_colm(level1)%p(i),I_cf_colm(level1)%p(i+1)-1
            j = I_cf_rowp(level1)%p(k)
            vf(i) = vf(i) + VC(j)*I_cf(level1)%p(k)
         enddo
      enddo

      if ( level1 .eq. -1 ) then
          call commIn(vf,ilwork,1,iper,iBC,BC)
          call commOut(vf,ilwork,1,iper,iBC,BC)
      end if

      call cpu_time(cpusec(2))
      ramg_time(13)=ramg_time(13)+cpusec(2)-cpusec(1)

      end subroutine ! ramg_calcIvCF

!***********************************************************
!      ramg_calcIvFC
!       calculate v fine to coarse,
!       VC = IT * v
!***********************************************************
      subroutine ramg_calcIvFC(level1,level2,vf,VC)
      use ramg_data
      implicit none
      integer,intent(in)                :: level1,level2
      real(kind=8),intent(in),dimension(amg_nshg(level1)) :: vf
      real(kind=8),intent(inout),dimension(amg_nshg(level2)) :: VC

      integer                           :: i,j,k,p
      real(kind=8) :: cpusec(2)

      call cpu_time(cpusec(1))

      VC = 0
      do i=1,amg_nshg(level2)
         do k=I_fc_colm(level1)%p(i),I_fc_colm(level1)%p(i+1)-1
            j = I_fc_rowp(level1)%p(k)
            VC(i) = VC(i) + vf(j)*I_fc(level1)%p(k)
         enddo
      enddo
      call cpu_time(cpusec(2))
      ramg_time(13)=ramg_time(13)+cpusec(2)-cpusec(1)

      end subroutine ! ramg_calcIvFC

!!***********************************************************
!      ramg_calcSvCF
!       calculate V coarse to fine,
!       v = S * VC
!***********************************************************

!***********************************************************
!      ramg_calcIvFC
!       calculate v fine to coarse,
!       VC = ST * v
!***********************************************************

!************************************************************
!      ramg_calcAv
!      calculate u = A*v
!************************************************************
      subroutine ramg_calcAv(gcolm,growp,glhs,gnshg,gnnz_tot,
     &                         u,v,gcomm)
      use ramg_data
      integer,intent(in)                 :: gnshg,gnnz_tot
      integer,intent(in),dimension(gnshg+1) :: gcolm
      integer,intent(in),dimension(gnnz_tot) :: growp
      real(kind=8),intent(in),dimension(gnnz_tot) :: glhs
      real(kind=8),intent(in),dimension(gnshg) :: v
      real(kind=8),intent(inout),dimension(gnshg) :: u
      integer,intent(in) :: gcomm

      integer                             :: i,j,k,p

      u = 0
      do i=1,gnshg
         do k=gcolm(i),gcolm(i+1)-1
            j = growp(k)
            u(i) = u(i)+glhs(k)*v(j)
         enddo
      enddo

!      if (gcomm.eq.1) then
!      end if

      end subroutine ! ramg_calcAv

!************************************************************
!      ramg_calcAv_g
!      calculate u = A*v
!      Globally: commout, do product, commin, (zeroout)
!************************************************************
      subroutine ramg_calcAv_g(level,u,v,colm,rowp,lhsK,lhsP,
     &           ilwork,BC,iBC,iper,gcomm)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer,intent(in),dimension(nlwork) :: ilwork
      integer,intent(in),dimension(nshg)              :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC)  :: BC
      integer,intent(in),dimension(nshg+1)       :: colm
      integer,intent(in),dimension(nnz_tot)      :: rowp

      real(kind=8),intent(in),dimension(9,nnz_tot)    :: lhsK
      real(kind=8),intent(in),dimension(4,nnz_tot)    :: lhsP
      integer,intent(in) :: level
      real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v
      real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u
      integer,intent(in) :: gcomm

      integer                             :: i,j,k,p,mi,mj
      real(kind=8)  :: cpusec(2)

      call cpu_time(cpusec(1))

      IF (level.eq.1) THEN
      !IF (.FALSE.) THEN
      call ramg_PPEAp(u,v,colm,rowp,lhsK,lhsP,ilwork,BC,iBC,iper)
      ELSE
      if (gcomm.eq.1) then
      call ramg_commOut(v,level,ilwork,1,iper,iBC,BC)
      endif
      u = 0
      do i=1,amg_nshg(level)
         mi = amg_paramap(level)%p(i)
         do k=amg_A_colm(level)%p(i),amg_A_colm(level)%p(i+1)-1
            j = amg_A_rowp(level)%p(k)
            mj = amg_paramap(level)%p(j)
            if (.not.( (mi.eq.mj).and.(mi.lt.0) )) then
                u(i) = u(i)+amg_A_lhs(level)%p(k,1)*v(j)
            endif
         enddo
      enddo
      if (gcomm.eq.1) then
      call ramg_commIn(u,level,ilwork,1,iper,iBC,BC)
      endif
      ENDIF
      call cpu_time(cpusec(2))
      if (level.eq.1) then
      ramg_time(11) = ramg_time(11) + cpusec(2)-cpusec(1)
      else
      ramg_time(12) = ramg_time(12) + cpusec(2)-cpusec(1)
      endif

      end subroutine ! ramg_calcAv_g

!************************************************************
!      ramg_calcAv_b: same as calcAv_g, but only apply
!                     on boundary nodes.
!      calculate u = A*v
!      Globally: commout, do product, commin, (zeroout)
!************************************************************
      subroutine ramg_calcAv_b(level,u,v,
     &           ilwork,BC,iBC,iper,gcomm,diag)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer,intent(in),dimension(nlwork) :: ilwork
      integer,intent(in),dimension(nshg)              :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC)  :: BC

      integer,intent(in) :: level
      real(kind=8),intent(inout),dimension(amg_nshg(level)) :: v
      real(kind=8),intent(inout),dimension(amg_nshg(level)) :: u
      integer,intent(in) :: gcomm
      integer,intent(in) :: diag !=1: NOT include diagonal A_(i,i)
                                 !=0: include diagonal

      integer                             :: i,j,k,p,mi,mj,mk
      real(kind=8)  :: cpusec(2)

      call cpu_time(cpusec(1))

      if (gcomm.eq.1) then
      call ramg_commOut(v,level,ilwork,1,iper,iBC,BC)
      endif
      u = 0
      do i=1,amg_nshg(level)
         mi = amg_paramap(level)%p(i)
         mk = amg_paraext(level)%p(i)
         IF (mk.ne.(myrank+1)) THEN ! only treat boundary nodes
         do k=amg_A_colm(level)%p(i)+diag,amg_A_colm(level)%p(i+1)-1
            j = amg_A_rowp(level)%p(k)
            mj = amg_paramap(level)%p(j)
            if (.not.( (mi.eq.mj).and.(mi.lt.0) )) then
                u(i) = u(i)+amg_A_lhs(level)%p(k,1)*v(j)
            endif
         enddo
         ELSE
             u(i) = v(i)
         ENDIF
      enddo
      if (gcomm.eq.1) then
      call ramg_commIn(u,level,ilwork,1,iper,iBC,BC)
      endif

      end subroutine ! ramg_calcAv_b

!*************************************************
!      check_paracomm
!*************************************************
      subroutine ramg_check_paracomm(ilwork,BC,iBC,iper)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer,intent(in),dimension(nlwork)        :: ilwork
      integer,intent(in),dimension(nshg)            :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC

      integer i,j
      type(r1d),dimension(ramg_levelx) :: tarray
      character :: fname*80

      do i=1,ramg_levelx
         allocate(tarray(i)%p(amg_nshg(i)))
         do j=1,amg_nshg(i)
            call random_number(tarray(i)%p)
         enddo
      enddo

      do i=1,ramg_levelx
         write(fname,'((A8)(I1))')'tarray_a',i
         call ramg_dump(tarray(i)%p,amg_nshg(i),fname)
      enddo

      do i=1,ramg_levelx
         call ramg_commIn(tarray(i)%p,i,ilwork,1,iper,iBC,BC)
         write(fname,'((A8)(I1))')'tarray_i',i
         call ramg_dump(tarray(i)%p,amg_nshg(i),fname)
      enddo

      do i=1,ramg_levelx
         call ramg_commOut(tarray(i)%p,i,ilwork,1,iper,iBC,BC)
         write(fname,'((A8)(I1))')'tarray_b',i
         call ramg_dump(tarray(i)%p,amg_nshg(i),fname)
      enddo

!      call commOut(tarray(1)%p,ilwork,1,iper,iBC,BC)

!      do i=1,ramg_levelx
!         write(fname,'((A8)(I1))')'tarray_i',i
!         call ramg_dump_rn_map(tarray(i)%p,amg_nshg(i),1,fname)
!      enddo

      do i=1,ramg_levelx
         deallocate(tarray(i)%p)
      enddo

      end subroutine ! ramg_check_paracomm

!**************************************************
!      ramg_zeroOut: zero out slave values
!**************************************************
      subroutine ramg_zeroOut(u,ilwork,BC,iBC,iper)
      include "common.h"
      include "auxmpi.h"
      include "mpif.h"

      real(kind=8),dimension(nshg),intent(inout)  :: u
      integer,intent(in),dimension(nlwork)        :: ilwork
      integer,intent(in),dimension(nshg)            :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC

      integer i,j,k,p
      integer            :: itag,iacc,iother,numseg,isgbeg,itkbeg
      integer            :: numtask,lenseg

      if (numpe.lt.2) return
      !call commOut(u,ilwork,1,iper,iBC,BC)
      numtask = ilwork(1)
      itkbeg = 1
      do i=1,numtask
         iacc = ilwork(itkbeg+2)
         numseg = ilwork(itkbeg+4)
         if (iacc.eq.0) then
         do j=1,numseg
            isgbeg = ilwork(itkbeg+3+j*2)
            lenseg = ilwork(itkbeg+4+j*2)
            isgend = isgbeg + lenseg -1
            u(isgbeg:isgend) = 0
         enddo
         endif
         itkbeg = itkbeg+4+2*numseg
      enddo

      end subroutine ! ramg_zeroOut

!************************************************
!      ramg_freeBC: set "fine" on boundary nodes
!************************************************
      subroutine ramg_freeBC(amg_F,ilwork,BC,iBC,iper)
      use ramg_data
      include "common.h"
      include "auxmpi.h"
      include "mpif.h"
      integer,intent(inout),dimension(nshg)  :: amg_F
      integer,intent(in),dimension(nlwork)     :: ilwork
      integer,intent(in),dimension(nshg)       :: iBC,iper
      real(kind=8),intent(in),dimension(nshg,ndofBC) :: BC
      integer      :: itag,iacc,iother,numseg,isgbeg,itkbeg,numtask
      integer      :: i,p
      integer,dimension(nshg,2)  :: tmpout

      character     :: filename*80
      integer       :: nentry,nshg,iglobal,ilocal,icf,iproc

      return

      if (numpe.ge.2) then

          numtask = ilwork(1)
          itkbeg = 1
          do i=1,numtask
             itag = ilwork(itkbeg+1)
             iacc = ilwork(itkbeg+2)
             iother = ilwork(itkbeg+3)
             numseg = ilwork(itkbeg+4)
             isgbeg = ilwork(itkbeg+5)
             do j=1,numseg
                p = ilwork(itkbeg+3+j*2)
                amg_F(p) = 2
             enddo
          enddo

      tmpout(:,1) = ncorp_map
      tmpout(:,2) = amg_F

      call ramg_dump_i(tmpout,nshg,2,'CFsplit   ')
      else  ! DEBUGGING PART, READ 2-PROC CASE INTO 1-proc
          filename = "cfs.dat"
          filename = trim(filename)
          open(unit=999,file=filename,status='unknown')
          read(999,*)nentry
          do i=1,nentry
             read(999,*)iglobal,ilocal,icf,iproc
             !amg_F(iglobal) = icf
          enddo
          close(999)
      end if

      end subroutine ! ramg_freeBC

!**********************************************
!      ramg_read_map:
!       read in ncorp array from local to global
!       mapping
!**********************************************
      subroutine ramg_prep_reduce_map
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      character*5 cname
      character*255 fname1
      integer igeom,map_nshg,i

      allocate(ncorp_map(nshg))
      if (numpe.lt.2) then! construct dummy-array
          do i=1,nshg
             ncorp_map(i) = i
          enddo
      else
      fname1='geombc.dat'
      fname1=trim(fname1)//cname(myrank+1)

      call openfile(fname1,"read",igeom)
      !write(*,*)'mcheck:',myrank,fname1,igeom

      fname1="mode number map from partition to global"
      ione=1
      call readheader(igeom,fname1,map_nshg,ione,"integer",iotype)

      call readdatablock(igeom,fname1,ncorp_map,map_nshg,
     &                   "integer",iotype)

      call closefile(igeom,"read")

      endif

      end subroutine ! ramg_read_map

!***********************************************
!      ramg_check_diag
!***********************************************
      subroutine ramg_check_diag(colm,rowp,lhs,nshg,nnz_tot)
      implicit none
      integer,intent(in)                   :: nshg,nnz_tot
      integer,intent(in),dimension(nshg+1) :: colm
      integer,intent(in),dimension(nnz_tot) :: rowp
      real(kind=8),intent(in),dimension(nnz_tot) :: lhs

      integer                              :: i,j,p
      real(kind=8)                         :: diagrow
      logical                              :: diagokay

      diagokay = .true.

      do i=1,nshg
         p = colm(i)
         if (rowp(p).ne.i) then
             write(*,*)'matrix first row entry is not diagonal',i
             diagokay = .false.
         end if
         diagrow = lhs(p)
         if (diagrow.lt.0) then
             write(*,*)'matrix diagonal < 0 at',i,diagrow
             diagokay = .false.
         end if
         do j = colm(i)+1,colm(i+1)-1
            p = rowp(j)
            if (lhs(p).gt.diagrow) then
                write(*,*)'matrix not diagonal dominant at row',i
                diagokay = .false.
                write(*,*)'diag:',diagrow,p,lhs(p)
                exit
            end if
         end do
      enddo

      if (diagokay) then
          write(*,*)'matrix check diagonal dominance okay'
      end if

      end subroutine

      subroutine ramg_dot_p(level,v,u,redun,norm)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer :: redun,level
      real(kind=8),intent(in),dimension(amg_nshg(level),redun) ::v,u
      real(kind=8),intent(inout)    :: norm
      integer                       :: i,k
      real(kind=8) :: normt
      normt = 0
      do i=1,amg_nshg(level)
         do k=1,redun
         normt = normt + v(i,k)*u(i,k)
         enddo
      enddo
      if (numpe.gt.1) then
      call MPI_AllReduce(normt,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM,
     & MPI_COMM_WORLD,ierr)
      else
          norm=normt
      endif
      end subroutine !ramg_dot_p

!************************************************************
!      ramg_L2_norm
!      calculate norm = | r |
!************************************************************
      subroutine ramg_L2_norm(nshg,v,norm)
      implicit none
      integer,intent(in)           :: nshg
      real(kind=8),intent(in),dimension(nshg) :: v
      real(kind=8),intent(inout)    :: norm
      integer                       :: i
      norm = 0
      do i=1,nshg
         norm = norm + v(i)*v(i)
      enddo
      norm = sqrt(norm)
      end subroutine !ramg_L2_norm

      subroutine ramg_L2_norm_p(level,v,vflow,norm)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer,intent(in) :: vflow,level
      real(kind=8),intent(in),dimension(amg_nshg(level),vflow) :: v
      real(kind=8),intent(inout)    :: norm
      integer                       :: i,k
      real(kind=8) :: normt
      normt = 0
      do i=1,amg_nshg(level)
         do k=1,vflow
         normt = normt + v(i,k)*v(i,k)
         enddo
      enddo
      if (numpe.gt.1) then
      call MPI_AllReduce(normt,norm,1,MPI_DOUBLE_PRECISION,MPI_SUM,
     & MPI_COMM_WORLD,ierr)
      else
          norm = normt
      endif
      norm = sqrt(norm)
      end subroutine !ramg_L2_norm_p

!************************************************************
!   ramg_allocate & deallocate
!    (de)allocate memory for level l, lhs matrix, rhs vector
!************************************************************
      subroutine ramg_allocate(
     &                level,lnshg,lnnz_tot,n_sol)
      use ramg_data
      implicit none

      integer,intent(in)             :: level,lnshg,lnnz_tot
      integer,intent(in)             :: n_sol
      integer                        :: mem_err,mem_err_s
      mem_err_s = 0
      if (lnnz_tot.ne.0) then ! zero means manual alloc later
      amg_nnz(level) = lnnz_tot
      allocate(amg_A_lhs(level)%p(amg_nnz(level),n_sol),
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      allocate(amg_A_rowp(level)%p(amg_nnz(level)),
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (lnshg.ne.0) then
      amg_nshg(level) = lnshg
      allocate(amg_A_colm(level)%p(amg_nshg(level)+1),
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      allocate(amg_A_rhs(level)%p(amg_nshg(level)),
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      allocate(amg_ppeDiag(level)%p(amg_nshg(level)),
     &        stat=mem_err)
      endif
      mem_err_s = mem_err_s + mem_err
      if (mem_err_s .ne. 0 ) then
          write(6,7001)level
      end if
7001  format(/' **** ramg: Allocation error at level',i5)
      ! zero out
      if (lnnz_tot.ne.0) then
      amg_A_lhs(level)%p(:,:)  = 0
      amg_A_rowp(level)%p(:) = 0
      endif
      if (lnshg.ne.0) then
      amg_A_colm(level)%p(:) = 0
      amg_A_rhs(level)%p(:)  = 0
      endif
      return

      end subroutine ! ramg_allocate

      subroutine ramg_deallocate(level)

      use ramg_data
      include "common.h"

      integer,intent(inout)            :: level
      integer                        :: mem_err,mem_err_s,i

      if (level.eq.1) then

      if (maxnev.gt.0) then
        deallocate(ramg_ggb_ev)
        deallocate(ramg_ggb_eA)
        deallocate(cmtxA)
        deallocate(cindx)
      endif

        if (iamg_reduce.gt.1) then
            deallocate(reducemap)
            deallocate(rmap1d)
        endif

      end if

      level = abs(level)

      mem_err_s = 0
      if (associated(amg_A_lhs(level)%p)) then
      deallocate(amg_A_lhs(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (associated(amg_A_rowp(level)%p)) then
      deallocate(amg_A_rowp(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (associated(amg_A_colm(level)%p)) then
      deallocate(amg_A_colm(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (associated(amg_ppeDiag(level)%p)) then
      deallocate(amg_ppeDiag(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (associated(amg_A_rhs(level)%p)) then
      deallocate(amg_A_rhs(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (associated(amg_ilwork(level)%p)) then
      deallocate(amg_ilwork(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      write(*,*)'mcheck deallocate ilwork,',level,myrank
      endif

      if (associated(amg_paramap(level)%p)) then
      deallocate(amg_paramap(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif
      if (associated(amg_paraext(level)%p)) then
      deallocate(amg_paraext(level)%p,
     &         stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      endif

      if (mem_err_s .ne. 0 ) then
          write(6,7002)level
      end if

      if (associated(I_fc_colm(level)%p)) then
      deallocate(CF_map(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(CF_revmap(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(I_cf_colm(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(I_cf_rowp(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(I_cf(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(I_fc_colm(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(I_fc_rowp(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      deallocate(I_fc(level)%p,stat=mem_err)
      mem_err_s = mem_err_s + mem_err
      end if
      amg_nnz(level) = 0
      amg_nshg(level) = 0
7002  format(/' **** ramg: Deallocation error at level',i5)

      return

      end subroutine ! ramg_deallocate

      subroutine ramg_readin_i(iarray,rn1,rfname)
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1
      integer,dimension(rn1) ::  iarray
      character*10     rfname
      character*5      cname

      character*20     mfname

      integer i,t
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG READ: ',mfname
      open(264,file=trim(mfname),status='unknown')
      do i=1,rn1
         read(264,*)t,iarray(i)
      enddo
      close(264)
      end subroutine ! ramg_readin_i

      subroutine ramg_dump(rarray,rn1,rfname)
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1
      real(kind=8),dimension(rn1) ::  rarray
      character*10     rfname
      character*5      cname

      character*20     mfname

      integer i
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG DUMP: ',mfname
      open(264,file=trim(mfname),status='unknown')
      do i=1,rn1
         !write(264,'((I10)(A)(D10.3))')i,'  ',rarray(i)
         write(264,*)i,rarray(i)
      enddo
      close(264)

      end subroutine ! ramg_dump

      subroutine ramg_dump_ic(rarray,rn1,rn2,rfname)
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1,rn2
      integer,dimension(rn1,rn2) ::  rarray
      character*10     rfname
      character*5      cname

      character*20     mfname

      integer i
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG DUMP: ',mfname
      open(264,file=trim(mfname),status='unknown')
      write(264,*)rn2
      do j=1,rn2
         write(264,*)j,(rarray(i,j),i=1,rn1)
      enddo
      close(264)
      end subroutine

      subroutine ramg_dump_i(rarray,rn1,rn2,rfname)
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1,rn2
      integer,dimension(rn1,rn2) ::  rarray
      character*10     rfname
      character*5      cname

      character*20     mfname

      integer i
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG DUMP: ',mfname
      open(264,file=trim(mfname),status='unknown')
      write(264,*)rn1
      do i=1,rn1
         write(264,*)i,(rarray(i,j),j=1,rn2)
      enddo
      close(264)

      end subroutine ! ramg_dump_i

      subroutine ramg_dump_ir(iarray,rarray,rn1,rn2,rfname)
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1,rn2
      integer,dimension(rn1) ::  iarray
      real(kind=8),dimension(rn1,rn2) ::  rarray
      character*10     rfname
      character*5      cname

      character*20     mfname

      character*20 mformat

      integer i
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(mformat,'((A6)(I1)(A7))')'((2I)(',rn2,'E14.3))'
      mformat = trim(mformat)
      write(*,*)'RAMG DUMP: ',mfname
      open(264,file=trim(mfname),status='unknown')
      do i=1,rn1
         write(264,mformat)
     &   i,iarray(i),(rarray(i,j),j=1,rn2)
      enddo
      close(264)

      end subroutine ! ramg_dump_ir

      subroutine ramg_dump_rn_map(rarray,rn1,rn2,rfname)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1,rn2
      real(kind=8),dimension(rn1,rn2) ::  rarray
      character*10     rfname
      character*5      cname

      character*20     mfname
      character*20 mformat

      integer i,j,ii
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG DUMP: ',mfname
      open(264,file=trim(mfname),status='unknown')
      write(mformat,'((A6)(I1)(A7))')'((1I)(',rn2,'E14.4))'
      mformat=trim(mformat)
      do i=1,rn1
         if (numpe.gt.1) then
             ii = ncorp_map(i)
         else
             ii = i
         endif
         write(264,mformat)ii,(rarray(i,j),j=1,rn2)
      enddo
      close(264)
      end subroutine ! ramg_dump_rn_map

      subroutine ramg_dump_rn(rarray,rn1,rn2,rfname)
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"
      integer rn1,rn2
      real(kind=8),dimension(rn1,rn2) ::  rarray
      character*10     rfname
      character*5      cname

      character*20     mfname
      character*20 mformat

      integer i,j
      mfname = trim(rfname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG DUMP: ',mfname
      open(264,file=trim(mfname),status='unknown')
      write(mformat,'((A6)(I1)(A7))')'((1I)(',rn2,'E14.4))'
      mformat=trim(mformat)
      do i=1,rn1
         write(264,mformat)i,(rarray(i,j),j=1,rn2)
      enddo
      close(264)

      end subroutine ! ramg_dump_rn

      subroutine ramg_dump_matlab_A(acolm,arowp,alhs,an,annz,redun,
     &           fname)

      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer :: an,annz,redun
      integer,dimension(an+1) :: acolm
      integer,dimension(annz) :: arowp
      real(kind=8),dimension(redun,annz) :: alhs
      character fname*10,mtxname*5
      character cname*5

      character mfname*15,mAname*5
      character mformat*20

      integer i,j,p,k

      mfname = trim(fname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG Dump to matlab  ',mfname,myrank
      open(264,file=trim(mfname),status='unknown')
      !write(264,*)an
      !write(264,*)annz
      !write(264,*)'1'
      write(mformat,'((A6)(I1)(A7))')'((2I)(',redun,'E14.4))'
      do i=1,an
         do p=acolm(i),acolm(i+1)-1
            j = arowp(p)
            !write(264,"((I6)(I7)(A)(E8.2))")i,j,' ',alhs(p)
!            if ( (amg_paramap(1)%p(i).eq.amg_paramap(2)%p(j)).and.
!           if   (iabs(amg_paramap(1)%p(i)).ne.(myrank+1)) then
!     &          (amg_paramap(1)%p(i).ne.1)) then
            write(264,mformat)i,j,(alhs(k,p),k=1,redun)
!            endif
         enddo
      enddo
      close(264)
      end subroutine

      subroutine ramg_dump_matlab_map(acolm,arowp,alhs,an,annz,redun,
     &           fname)
      use ramg_data
      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      integer :: an,annz,redun
      integer,dimension(an+1) :: acolm
      integer,dimension(annz) :: arowp
      real(kind=8),dimension(redun,annz) :: alhs
      !integer,dimension(an) :: ipmap
      character fname*10

      character mfname*20
      character cname*5
      character mformat*20

      integer i,j,p,k
      integer ii,jj

      if (numpe.eq.1) then
          call ramg_dump_matlab_A(acolm,arowp,alhs,an,annz,redun,
     &         fname)
         return;
      endif

      mfname = trim(fname)//'.dat'//cname(myrank+1)

      write(*,*)'RAMG Dump to matlab  ',mfname,nshg,myrank
      open(264,file=trim(mfname),status='unknown')
      !write(264,*)an
      !write(264,*)annz
      !write(264,*)'1'
      write(mformat,'((A6)(I1)(A7))')'((2I)(',redun,'E14.4))'
      !call ramg_dump_i(ncorp_map,nshg,1,'pam_corp  ')
      do i=1,an
         ii = ncorp_map(i)!ipmap(i))
         do p=acolm(i),acolm(i+1)-1
            j = arowp(p)
            jj = ncorp_map(j)!ipmap(j))
            !write(264,"((I6)(I7)(A)(E8.2))")i,j,' ',alhs(p)
            write(264,mformat)ii,jj,(alhs(k,p),k=1,redun)
         enddo
      enddo
      close(264)
      end subroutine

      subroutine ramg_dump_mupad_A(acolm,arowp,alhs,an,annz,
     &           fname,mtxname)
      integer :: an,annz
      integer,dimension(an+1) :: acolm
      integer,dimension(annz) :: arowp
      real(kind=8),dimension(annz) :: alhs
      character fname*10,mtxname*5

      character mfname*15,mAname*5

      integer i,j,p,k

      mfname = trim(fname)//'.mws'
      mAname = trim(mtxname)

      write(*,*)'RAMG Dump to Mupad  ',mfname
      open(264,file=trim(mfname),status='unknown')
      write(264,*)mAname,':= matrix(',an,',',an,'):'
      do i=1,an
         do p=acolm(i),acolm(i+1)-1
            j = arowp(p)
            write(264,*)mAname,'[',i,',',j,']:=',alhs(p),':'
         enddo
      enddo
      close(264)

      end subroutine

      subroutine ramg_print_time(str,v)

      include "common.h"
      include "mpif.h"
      include "auxmpi.h"

      character*(*) str
      real(kind=8) :: v,vave,vmax

      if (numpe.gt.1) then
          call MPI_AllReduce(v,vmax,1,
     &    MPI_DOUBLE_PRECISION,MPI_MAX,MPI_COMM_WORLD,ierr)
          call MPI_AllReduce(v,vave,1,
     &    MPI_DOUBLE_PRECISION,MPI_SUM,MPI_COMM_WORLD,ierr)
          vave = vave/numpe
      else
          vave = v
          vmax = v
      endif

      if ((iamg_verb.gt.1).and.(myrank.eq.master)) then
      write(*,7101)trim(str),vave,vmax
7101  format(T1,A,T40,f9.2,' s (ave), ',f9.2,' s (max)')
      endif

      end subroutine

      subroutine ramg_output_coarsening
      use readarrays
      use ramg_data
      include "common.h"
      include "mpif.h"

      character*20 mfname
      character*20 mformat

      integer i

      write(*,*)'ramg dump coarsening'
      mfname = 'amgcoarsen.dat'
      open(265,file=trim(mfname),status='unknown')
      write(265,*)nshg
      do i=1,nshg
      write(265,'((2I)(3E14.5))')i,amg_cfmap(i),
     & point2x(i,1),point2x(i,2),point2x(i,3)
      enddo
      close(265)

      end subroutine

!***********************************************************
!      <EOF> ramg TOOLS
!**********************************************************

ext: f


text:
C$Procedure      BODVAR ( Return values from the kernel pool )

      SUBROUTINE BODVAR ( BODY, ITEM, DIM, VALUES )

C$ Abstract
C
C     Deprecated: This routine has been superseded by BODVCD and
C     BODVRD.  This routine is supported for purposes of backward
C     compatibility only.
C
C     Return the values of some item for any body in the
C     kernel pool.
C
C$ Disclaimer
C
C     THIS SOFTWARE AND ANY RELATED MATERIALS WERE CREATED BY THE
C     CALIFORNIA INSTITUTE OF TECHNOLOGY (CALTECH) UNDER A U.S.
C     GOVERNMENT CONTRACT WITH THE NATIONAL AERONAUTICS AND SPACE
C     ADMINISTRATION (NASA). THE SOFTWARE IS TECHNOLOGY AND SOFTWARE
C     PUBLICLY AVAILABLE UNDER U.S. EXPORT LAWS AND IS PROVIDED "AS-IS"
C     TO THE RECIPIENT WITHOUT WARRANTY OF ANY KIND, INCLUDING ANY
C     WARRANTIES OF PERFORMANCE OR MERCHANTABILITY OR FITNESS FOR A
C     PARTICULAR USE OR PURPOSE (AS SET FORTH IN UNITED STATES UCC
C     SECTIONS 2312-2313) OR FOR ANY PURPOSE WHATSOEVER, FOR THE
C     SOFTWARE AND RELATED MATERIALS, HOWEVER USED.
C
C     IN NO EVENT SHALL CALTECH, ITS JET PROPULSION LABORATORY, OR NASA
C     BE LIABLE FOR ANY DAMAGES AND/OR COSTS, INCLUDING, BUT NOT
C     LIMITED TO, INCIDENTAL OR CONSEQUENTIAL DAMAGES OF ANY KIND,
C     INCLUDING ECONOMIC DAMAGE OR INJURY TO PROPERTY AND LOST PROFITS,
C     REGARDLESS OF WHETHER CALTECH, JPL, OR NASA BE ADVISED, HAVE
C     REASON TO KNOW, OR, IN FACT, SHALL KNOW OF THE POSSIBILITY.
C
C     RECIPIENT BEARS ALL RISK RELATING TO QUALITY AND PERFORMANCE OF
C     THE SOFTWARE AND ANY RELATED MATERIALS, AND AGREES TO INDEMNIFY
C     CALTECH AND NASA FOR ALL THIRD-PARTY CLAIMS RESULTING FROM THE
C     ACTIONS OF RECIPIENT IN THE USE OF THE SOFTWARE.
C
C$ Required_Reading
C
C     KERNEL
C
C$ Keywords
C
C     CONSTANTS
C
C$ Declarations

      INTEGER               BODY
      CHARACTER*(*)         ITEM
      INTEGER               DIM
      DOUBLE PRECISION      VALUES  ( * )

C$ Brief_I/O
C
C     VARIABLE  I/O  DESCRIPTION
C     --------  ---  --------------------------------------------------
C     BODY       I   ID code of body.
C     ITEM       I   Item for which values are desired. ('RADII',
C                    'NUT_PREC_ANGLES', etc. )
C     DIM        O   Number of values returned.
C     VALUES     O   Values.
C
C$ Detailed_Input
C
C     BODY       is the ID code of the body for which ITEM is
C                requested. Bodies are numbered according to the
C                standard NAIF numbering scheme.
C
C     ITEM       is the item to be returned. Together, the body and
C                item name combine to form a variable name, e.g.,
C
C                      'BODY599_RADII'
C                      'BODY401_POLE_RA'
C
C$ Detailed_Output
C
C     DIM        is the number of values associated with the variable.
C
C     VALUES     are the values associated with the variable.
C
C$ Parameters
C
C     None.
C
C$ Exceptions
C
C    None.
C
C$ Files
C
C     None.
C
C$ Particulars
C
C     None.
C
C$ Examples
C
C     The call
C
C         CALL BODVAR ( 399, 'RADII', DIM, VALUE )
C
C     returns the dimension and values associated with the variable
C     'BODY399_RADII', for example,
C
C          DIM      = 3
C          VALUE(1) = 6378.140
C          VALUE(2) = 6378.140
C          VALUE(3) = 6356.755
C
C$ Restrictions
C
C     1) If the requested item is not found, the error
C        SPICE(KERNELVARNOTFOUND) is signalled.
C
C$ Literature_References
C
C      1) Refer to the SPK required reading file for a complete list of
C         the NAIF integer ID codes for bodies.
C
C$ Author_and_Institution
C
C     W.L. Taber      (JPL)
C     I.M. Underwood  (JPL)
C
C$ Version
C
C-     SPICELIB Version 1.0.5, 18-MAY-2010 (BVS)
C
C        Index lines now state that this routine is deprecated.
C
C-     SPICELIB Version 1.0.4, 27-OCT-2005 (NJB)
C
C         Routine is now deprecated.
C
C-     SPICELIB Version 1.0.3, 08-JAN-2004 (EDW)
C
C         Trivial typo corrected.
C
C-     SPICELIB Version 1.0.2, 10-MAR-1992 (WLT)
C
C         Comment section for permuted index source lines was added
C         following the header.
C
C-     SPICELIB Version 1.0.1, 8-AUG-1990 (HAN)
C
C           Detailed Input section of the header was updated. The
C           description for the variable BODY was incorrect.
C
C-     SPICELIB Version 1.0.0, 31-JAN-1990  (WLT) (IMU)
C
C-&

C$ Index_Entries
C
C     DEPRECATED fetch constants for a body from the kernel pool
C     DEPRECATED physical constants for a body
C
C-&

C
C     SPICELIB functions
C
      LOGICAL               RETURN

C
C     Local variables
C
      CHARACTER*16          CODE
      CHARACTER*32          VARNAM
      LOGICAL               FOUND

C
C     Standard SPICE error handling.
C
      IF ( RETURN () ) THEN
        RETURN
      ELSE
        CALL CHKIN ( 'BODVAR' )
      END IF

C
C     Construct the variable name from BODY and ITEM.
C
      VARNAM = 'BODY'

      CALL INTSTR ( BODY, CODE )
      CALL SUFFIX ( CODE, 0, VARNAM )
      CALL SUFFIX ( '_',  0, VARNAM )
      CALL SUFFIX ( ITEM, 0, VARNAM )

C
C     Grab the items. Complain if they aren't there.
C
      CALL RTPOOL ( VARNAM, DIM, VALUES, FOUND )

      IF ( .NOT. FOUND ) THEN
         CALL SETMSG ( 'The variable # could not be found in the ' //
     .                 'kernel pool.'                             )
         CALL ERRCH  ( '#', VARNAM                 )
         CALL SIGERR ( 'SPICE(KERNELVARNOTFOUND)' )
      END IF

      CALL CHKOUT ( 'BODVAR' )
      RETURN
      END

max_stars_repo_name: agforero/FTFramework

max_issues_repo_name: agforero/fortran-testing-framework

max_forks_repo_name: agforero/fortran-testing-framework


! Copyright (c) 2015-2021, the ELSI team.
! All rights reserved.
!
! This file is part of ELSI and is distributed under the BSD 3-clause license,
! which may be found in the LICENSE file in the ELSI root directory.

!>
!! Provide implementations of the heapsort algorithm.
!!
module ELSI_SORT

   use ELSI_PRECISION, only: r8,i4,i8

   implicit none

   private

   public :: elsi_heapsort
   public :: elsi_permute
   public :: elsi_unpermute

   interface swap
      module procedure swap_i8
      module procedure swap_i4
      module procedure swap_r8
   end interface

   interface downheap
      module procedure downheap_i4
      module procedure downheap_i8
      module procedure downheap_r8
   end interface

   interface elsi_heapsort
      module procedure heapsort_i4
      module procedure heapsort_i8
      module procedure heapsort_r8
   end interface

   interface elsi_permute
      module procedure permute_i4
      module procedure permute_r8
      module procedure permute_c16
   end interface

   interface elsi_unpermute
      module procedure unpermute_r8
      module procedure unpermute_c16
   end interface

contains

!>
!! Sort a real(kind=r8) array by heapsort, returns an array of permutation.
!!
subroutine heapsort_r8(length,array,perm)

   implicit none

   integer(kind=i4), intent(in) :: length
   real(kind=r8), intent(inout) :: array(length)
   integer(kind=i4), intent(out) :: perm(length)

   integer(kind=i4) :: top
   integer(kind=i4) :: i

   do i = 1,length
      perm(i) = i
   end do

   top = length/2

   do i = top,1,-1
      call downheap(length,array,perm,i,length)
   end do

   i = length

   do while(i > 1)
      call swap(length,array,1,i)
      call swap(length,perm,1,i)

      i = i-1

      call downheap(length,array,perm,1,i)
   end do

end subroutine

!>
!! Restore the max-heap structure.
!!
subroutine downheap_r8(length,a_r8,b_i4,top,bottom)

   implicit none

   integer(kind=i4), intent(in) :: length
   real(kind=r8), intent(inout) :: a_r8(length)
   integer(kind=i4), intent(inout) :: b_i4(length)
   integer(kind=i4), intent(in) :: top
   integer(kind=i4), intent(in) :: bottom

   integer(kind=i4) :: v
   integer(kind=i4) :: w

   v = top
   w = 2*v

   do while(w <= bottom)
      if(w+1 <= bottom) then
         if(a_r8(w+1) > a_r8(w)) then
            w = w+1
         end if
      end if

      if(a_r8(v) >= a_r8(w)) then
         return
      else
         call swap(length,a_r8,v,w)
         call swap(length,b_i4,v,w)

         v = w
         w = 2*v
      end if
   end do

end subroutine

!>
!! Sort an integer(kind=i8) array by heapsort, returns an array of permutation.
!!
subroutine heapsort_i8(length,array,perm)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i8), intent(inout) :: array(length)
   integer(kind=i4), intent(out) :: perm(length)

   integer(kind=i4) :: top
   integer(kind=i4) :: i

   do i = 1,length
      perm(i) = i
   end do

   top = length/2

   do i = top,1,-1
      call downheap(length,array,perm,i,length)
   end do

   i = length

   do while(i > 1)
      call swap(length,array,1,i)
      call swap(length,perm,1,i)

      i = i-1

      call downheap(length,array,perm,1,i)
   end do

end subroutine

!>
!! Restore the max-heap structure.
!!
subroutine downheap_i8(length,a_i8,b_i4,top,bottom)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i8), intent(inout) :: a_i8(length)
   integer(kind=i4), intent(inout) :: b_i4(length)
   integer(kind=i4), intent(in) :: top
   integer(kind=i4), intent(in) :: bottom

   integer(kind=i4) :: v
   integer(kind=i4) :: w

   v = top
   w = 2*v

   do while(w <= bottom)
      if(w+1 <= bottom) then
         if(a_i8(w+1) > a_i8(w)) then
            w = w+1
         end if
      end if

      if(a_i8(v) >= a_i8(w)) then
         return
      else
         call swap(length,a_i8,v,w)
         call swap(length,b_i4,v,w)

         v = w
         w = 2*v
      end if
   end do

end subroutine

!>
!! Sort an integer(kind=i4) array by heapsort, returns an array of permutation.
!!
subroutine heapsort_i4(length,array,perm)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(inout) :: array(length)
   integer(kind=i4), intent(out) :: perm(length)

   integer(kind=i4) :: top
   integer(kind=i4) :: i

   do i = 1,length
      perm(i) = i
   end do

   top = length/2

   do i = top,1,-1
      call downheap(length,array,perm,i,length)
   end do

   i = length

   do while(i > 1)
      call swap(length,array,1,i)
      call swap(length,perm,1,i)

      i = i-1

      call downheap(length,array,perm,1,i)
   end do

end subroutine

!>
!! Restore the max-heap structure.
!!
subroutine downheap_i4(length,a_i4,b_i4,top,bottom)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(inout) :: a_i4(length)
   integer(kind=i4), intent(inout) :: b_i4(length)
   integer(kind=i4), intent(in) :: top
   integer(kind=i4), intent(in) :: bottom

   integer(kind=i4) :: v
   integer(kind=i4) :: w

   v = top
   w = 2*v

   do while(w <= bottom)
      if(w+1 <= bottom) then
         if(a_i4(w+1) > a_i4(w)) then
            w = w+1
         end if
      end if

      if(a_i4(v) >= a_i4(w)) then
         return
      else
         call swap(length,a_i4,v,w)
         call swap(length,b_i4,v,w)

         v = w
         w = 2*v
      end if
   end do

end subroutine

!>
!! Permute an integer(kind=i4) array.
!!
subroutine permute_i4(length,perm,array)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(in) :: perm(length)
   integer(kind=i4), intent(inout) :: array(length)

   integer(kind=i4) :: i
   integer(kind=i4) :: work(length)

   do i = 1,length
      work(i) = array(perm(i))
   end do

   array = work

end subroutine

!>
!! Permute a real(kind=r8) array.
!!
subroutine permute_r8(length,perm,array)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(in) :: perm(length)
   real(kind=r8), intent(inout) :: array(length)

   real(kind=r8) :: work(length)
   integer(kind=i4) :: i

   do i = 1,length
      work(i) = array(perm(i))
   end do

   array = work

end subroutine

!>
!! Permute a complex(kind=r8) array.
!!
subroutine permute_c16(length,perm,array)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(in) :: perm(length)
   complex(kind=r8), intent(inout) :: array(length)

   complex(kind=r8) :: work(length)
   integer(kind=i4) :: i

   do i = 1,length
      work(i) = array(perm(i))
   end do

   array = work

end subroutine

!>
!! Unpermute a real(kind=r8) array.
!!
subroutine unpermute_r8(length,perm,array)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(in) :: perm(length)
   real(kind=r8), intent(inout) :: array(length)

   real(kind=r8) :: work(length)
   integer(kind=i4) :: i

   do i = 1,length
      work(perm(i)) = array(i)
   end do

   array = work

end subroutine

!>
!! Unpermute a complex(kind=r8) array.
!!
subroutine unpermute_c16(length,perm,array)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(in) :: perm(length)
   complex(kind=r8), intent(inout) :: array(length)

   complex(kind=r8) :: work(length)
   integer(kind=i4) :: i

   do i = 1,length
      work(perm(i)) = array(i)
   end do

   array = work

end subroutine

!>
!! Swap two numbers in an integer(kind=i8) array.
!!
subroutine swap_i8(length,array,i,j)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i8), intent(inout) :: array(length)
   integer(kind=i4), intent(in) :: i
   integer(kind=i4), intent(in) :: j

   integer(kind=i8) :: tmp

   tmp = array(i)
   array(i) = array(j)
   array(j) = tmp

end subroutine

!>
!! Swap two numbers in an integer(kind=i4) array.
!!
subroutine swap_i4(length,array,i,j)

   implicit none

   integer(kind=i4), intent(in) :: length
   integer(kind=i4), intent(inout) :: array(length)
   integer(kind=i4), intent(in) :: i
   integer(kind=i4), intent(in) :: j

   integer(kind=i4) :: tmp

   tmp = array(i)
   array(i) = array(j)
   array(j) = tmp

end subroutine

!>
!! Swap two numbers in a real(kind=r8) array.
!!
subroutine swap_r8(length,array,i,j)

   implicit none

   integer(kind=i4), intent(in) :: length
   real(kind=r8), intent(inout) :: array(length)
   integer(kind=i4), intent(in) :: i
   integer(kind=i4), intent(in) :: j

   real(kind=r8) :: tmp

   tmp = array(i)
   array(i) = array(j)
   array(j) = tmp

end subroutine

end module ELSI_SORT


text:
!----------------------------------------------------------------------------
! This file is part of UCLALES.
!
! UCLALES 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; either version 3 of the License, or
! (at your option) any later version.
!
! UCLALES 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, see <http://www.gnu.org/licenses/>.
!
! Copyright 1999-2007, Bjorn B. Stevens, Dep't Atmos and Ocean Sci, UCLA
!----------------------------------------------------------------------------
!
!
module defs

  integer :: nv=1, nv1=2, mb=1
  real    :: totalpower
  real, parameter :: R      = 287.04
  real, parameter :: Rm     = 461.5
  real, parameter :: ep     = R/Rm
  real, parameter :: ep2    = Rm/R - 1.
  real, parameter :: cp     = 1005.
  real, parameter :: cv     = cp-R
  real, parameter :: rcp    = R/cp
  real, parameter :: cpr    = cp/R
  real, parameter :: g      = 9.8
  real, parameter :: p00    = 1.e+05
  real, parameter :: p00i   = 1./p00
  real, parameter :: omega  = 7.292e-05
  real, parameter :: alvl   = 2.5e+06  ! latent heat of vaporization
  real, parameter :: alvi   = 2.834e+06 ! latent heat of sublimation
  real, parameter :: rowt   = 1.e+3
  real, parameter :: roice  = 0.9e+3
  real, parameter :: vonk   = 0.40
  real, parameter :: stefan = 5.6696e-8
  real, parameter :: SolarConstant  = 1.365d+3
  real, parameter :: mair   = 28.967 ! molar mass of air
  real, parameter :: pi     = 3.14159265358979323846264338327
  real, parameter :: kb     = 1.3806e-23 !Boltzman constant

end module defs

max_stars_repo_name: huangynj/UCLALES-SALSA

max_issues_repo_name: huangynj/UCLALES-SALSA

max_forks_repo_name: huangynj/UCLALES-SALSA


text:
program raytr
  use, intrinsic :: iso_fortran_env, dp=>real64
  use geq_sphax_base
  include "gensrc/raytr-metricname.f90"
  use indiv_ngeod_solver
  implicit none
  real(dp), parameter :: pi = 3.1415192d0
  real(dp), dimension(0:3,0:3) :: g
  real(dp), dimension(1:3) :: pos, k, k0
  integer :: outputfile,logfile,metafile
  real(dp) :: dlmin, tol, sepgeod
  character(len=128) :: dataid, metricname
  integer :: view,i,imax,imin,ifactor

  selected_pot_func => get_potential
  selected_dpot_func => get_deriv_potential
  selected_get_metric => get_metric
  selected_geqrhs => geqrhs
  selected_surf => surface_sph

  indiv_ngeod_solver__logfile = logfile

  read(*,'(a)') dataid
  read(*,'(a)') metricname

  read(*,*) m
  read(*,*) a
  read(*,*) q
  read(*,*) asp_ratio
  read(*,*) view
  read(*,*) tol
  read(*,*) dlmin

  include "gensrc/raytr-paramtransf.f90"

  open(newunit=logfile, file="log.txt", status="replace", action="write")
  open(newunit=outputfile, file="data/raytr/"//trim(dataid)//"/"//trim(metricname)//".dat",&
    & form="unformatted", access="stream", status="replace")
  open(newunit=metafile, file="data/raytr/"//trim(dataid)//"/meta/"//trim(metricname)//".txt",&
    & status="replace", action="write")

  if (asp_ratio < 0 .and. asp_ratio >= -1) then
    !! Calculate aspect ratio, assuming solid body =====
    !!   >>> J2 = (Iz - Ix)/M*R^2 = Q/(M_g*R^2);  R=Req
    !!  If solid ellipsoid, Iz = 2/5*M*R^2, Ix = 1/5*M*(R^2 + Rp^2)
    !!      => (M_g/R) - (M_g/R)*Rp^2/R^2 = (Q/R^3)*5/2
    !!  /!\ The abs is added assuming that's always oblate!
    asp_ratio = sqrt(1-5*abs(q)/(2*m))
  end if

  surf_r_eq = 1d0
  coord_bl_a = a
  call set_coordinates(coord_bl)

  ! print*, get_potential([1d0,10d0,pi/2,0d0])

  include "print_art.f90"

  select case(view)
  case(1) !! xy plane
  !! Horizontal

    ifactor = 4
    imax = 640/ifactor
    imin = -640/ifactor

    do i = imin,imax
      write(*,'(a,a,f5.1,a,$)') char(13), "computing...", (imax+i)/(1.0*abs(imin)+abs(imax))*100.0, "%"

      sepgeod = 0.003125*ifactor
      pos = [15d0,sepgeod*real(i,dp),0d0]
      k = [-1d0,0d0,0d0]
      k0 = k

      !call ngeod_cart_rkf(pos,k,1d0,40d0,1d-5,2d0,1d-7,outputfile)
      ! call ngeod_cart_rkf(pos,k,1d0,50d0,1d-9,2d0,1d-11,outputfile)
      ! call ngeod_cart_rkf(pos,k,1d0,50d0,1d-8,5d0,1d-10,outputfile)
      call ngeod_cart_rkf(pos,k,1d0,50d0,dlmin,5d0,tol,outputfile)
      !write(logfile,*) i,0.25*i,k
      ! !  ngeod_cart_rkf(r  , k ,om,lmax,dlmin,dlmax,tol,outputfile)
      ! print*, acos(dot_product(k0,k)/sqrt(dot_product(k,k)*dot_product(k0,k0)))*180/pi
    end do

  case(2)
    !! Vertical

    ifactor = 64
    imax = 640/ifactor
    imin = -640/ifactor

    do i = imin,imax
      write(*,'(a,a,f5.1,a,$)') char(13), "computing...", (imax+i)/(1.0*abs(imin)+abs(imax))*100.0, "%"

      sepgeod = 0.003125*ifactor
      pos = [6d0,0.02d0,sepgeod*real(i,dp)]
      k = [-1d0,0d0,0d0]
      k0 = k
      ! call ngeod_cart_eul(pos,k,1d0,50d0,0.1d0,outputfile)

      ! call ngeod_cart_rkf(pos,k,1d0,40d0,1d-5,2d0,1d-7,outputfile)
      ! call ngeod_cart_rkf(pos,k,1d0,50d0,1d-8,2d0,1d-12,outputfile)

      call ngeod_cart_rkf(pos,k,1d0,20d0,dlmin,2d0,tol,outputfile)
      ! !  ngeod_cart_rkf(r  , k ,om,lmax,dlmin,dlmax,tol,outputfile)
      ! print*, acos(dot_product(k0,k)/sqrt(dot_product(k,k)*dot_product(k0,k0)))*180/pi
    end do
  end select

print*, ""
! print '(a)', "[Report]"
write(metafile,*) " dataid = ", dataid
write(metafile,*) " metricname = ", metricname
write(metafile,*) " m = ", m
write(metafile,*) " a = ", a
write(metafile,*) " q = ", q
write(metafile,*) " asp_ratio =", asp_ratio
write(metafile,*) " view = ", view
write(metafile,*) " dlmin = ", dlmin
write(metafile,*) " tol = ", tol
write(metafile,*) " imax = ", imax
write(metafile,*) " sepgeod = ", sepgeod
write(metafile,*) "8< - - - - - - - - - - - - - - -"

end program

max_stars_repo_name: gandreoliva/ujti

max_issues_repo_name: gandreoliva/ujti

max_forks_repo_name: gandreoliva/ujti


text:
      SUBROUTINE DLUM_3(LUM)
C     ****************************************************
C
C     Generated by MadGraph5_aMC@NLO v. %(version)s, %(date)s
C     By the MadGraph5_aMC@NLO Development Team
C     Visit launchpad.net/madgraph5 and amcatnlo.web.cern.ch
C     RETURNS PARTON LUMINOSITIES FOR MADFKS
C
C
C     Process: u g > t t~ u [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: c g > t t~ c [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: d g > t t~ d [ real = QCD QED ] QCD^2<=6 QED^2<=0
C     Process: s g > t t~ s [ real = QCD QED ] QCD^2<=6 QED^2<=0
C
C     ****************************************************
C
      IMPLICIT NONE
C
C     CONSTANTS
C
C
      INCLUDE 'genps.inc'
      INCLUDE 'nexternal.inc'
      DOUBLE PRECISION       CONV
      PARAMETER (CONV=389379660D0)  !CONV TO PICOBARNS
C
C     ARGUMENTS
C
C
      DOUBLE PRECISION LUM
C
C     LOCAL VARIABLES
C
C
      INTEGER I, ICROSS,LP
      DOUBLE PRECISION D1,U1,S1,C1
      DOUBLE PRECISION G2
C
C     EXTERNAL FUNCTIONS
C
C
      DOUBLE PRECISION PDG2PDF
C
C     GLOBAL VARIABLES
C
C
      INTEGER              IPROC
      DOUBLE PRECISION PD(0:MAXPROC)
      COMMON /SUBPROC/ PD, IPROC
      INCLUDE 'coupl.inc'
      INCLUDE 'run.inc'
      INTEGER IMIRROR
      COMMON/CMIRROR/IMIRROR
C
C     STUFF FOR DRESSED EE COLLISIONS
C
      INCLUDE 'eepdf.inc'
      DOUBLE PRECISION EE_COMP_PROD
      DOUBLE PRECISION DUMMY_COMPONENTS(N_EE)
      DOUBLE PRECISION D1_COMPONENTS(N_EE),U1_COMPONENTS(N_EE)
     $ ,S1_COMPONENTS(N_EE),C1_COMPONENTS(N_EE)
      DOUBLE PRECISION G2_COMPONENTS(N_EE)

      INTEGER I_EE
C
C
C
C     Common blocks
      CHARACTER*7         PDLABEL,EPA_LABEL
      INTEGER       LHAID
      COMMON/TO_PDF/LHAID,PDLABEL,EPA_LABEL
C
C     DATA
C
C
      DATA D1,U1,S1,C1/4*1D0/
      DATA G2/1*1D0/
      DATA ICROSS/1/
C     ----------
C
C     BEGIN CODE
C
C     ----------
C
      LUM = 0D0
      IF (ABS(LPP(1)) .GE. 1) THEN
        D1=PDG2PDF(LPP(1),1,1,XBK(1),DSQRT(Q2FACT(1)))
        IF ((ABS(LPP(1)).EQ.4.OR.ABS(LPP(1)).EQ.3)
     $   .AND.PDLABEL.NE.'none') D1_COMPONENTS(1:N_EE) =
     $    EE_COMPONENTS(1:N_EE)
        U1=PDG2PDF(LPP(1),2,1,XBK(1),DSQRT(Q2FACT(1)))
        IF ((ABS(LPP(1)).EQ.4.OR.ABS(LPP(1)).EQ.3)
     $   .AND.PDLABEL.NE.'none') U1_COMPONENTS(1:N_EE) =
     $    EE_COMPONENTS(1:N_EE)
        S1=PDG2PDF(LPP(1),3,1,XBK(1),DSQRT(Q2FACT(1)))
        IF ((ABS(LPP(1)).EQ.4.OR.ABS(LPP(1)).EQ.3)
     $   .AND.PDLABEL.NE.'none') S1_COMPONENTS(1:N_EE) =
     $    EE_COMPONENTS(1:N_EE)
        C1=PDG2PDF(LPP(1),4,1,XBK(1),DSQRT(Q2FACT(1)))
        IF ((ABS(LPP(1)).EQ.4.OR.ABS(LPP(1)).EQ.3)
     $   .AND.PDLABEL.NE.'none') C1_COMPONENTS(1:N_EE) =
     $    EE_COMPONENTS(1:N_EE)
      ENDIF
      IF (ABS(LPP(2)) .GE. 1) THEN
        G2=PDG2PDF(LPP(2),0,2,XBK(2),DSQRT(Q2FACT(2)))
        IF ((ABS(LPP(2)).EQ.4.OR.ABS(LPP(2)).EQ.3)
     $   .AND.PDLABEL.NE.'none') G2_COMPONENTS(1:N_EE) =
     $    EE_COMPONENTS(1:N_EE)
      ENDIF
      PD(0) = 0D0
      IPROC = 0
      IPROC=IPROC+1  ! u g > t t~ u
      PD(IPROC) = U1*G2
      IF (ABS(LPP(1)).EQ.ABS(LPP(2)).AND. (ABS(LPP(1))
     $ .EQ.3.OR.ABS(LPP(1)).EQ.4).AND.PDLABEL.NE.'none')PD(IPROC)
     $ =EE_COMP_PROD(U1_COMPONENTS,G2_COMPONENTS)
      IPROC=IPROC+1  ! c g > t t~ c
      PD(IPROC) = C1*G2
      IF (ABS(LPP(1)).EQ.ABS(LPP(2)).AND. (ABS(LPP(1))
     $ .EQ.3.OR.ABS(LPP(1)).EQ.4).AND.PDLABEL.NE.'none')PD(IPROC)
     $ =EE_COMP_PROD(C1_COMPONENTS,G2_COMPONENTS)
      IPROC=IPROC+1  ! d g > t t~ d
      PD(IPROC) = D1*G2
      IF (ABS(LPP(1)).EQ.ABS(LPP(2)).AND. (ABS(LPP(1))
     $ .EQ.3.OR.ABS(LPP(1)).EQ.4).AND.PDLABEL.NE.'none')PD(IPROC)
     $ =EE_COMP_PROD(D1_COMPONENTS,G2_COMPONENTS)
      IPROC=IPROC+1  ! s g > t t~ s
      PD(IPROC) = S1*G2
      IF (ABS(LPP(1)).EQ.ABS(LPP(2)).AND. (ABS(LPP(1))
     $ .EQ.3.OR.ABS(LPP(1)).EQ.4).AND.PDLABEL.NE.'none')PD(IPROC)
     $ =EE_COMP_PROD(S1_COMPONENTS,G2_COMPONENTS)
      DO I=1,IPROC
        IF (NINCOMING.EQ.2) THEN
          LUM = LUM + PD(I) * CONV
        ELSE
          LUM = LUM + PD(I)
        ENDIF
      ENDDO
      RETURN
      END


text:
      MODULE diagi_I
      INTERFACE
!...Generated by Pacific-Sierra Research 77to90  4.4G  10:47:08  03/09/06
      REAL(KIND(0.0D0)) FUNCTION diagi (IALPHA, IBETA, EIGA, XY, NMOS)
      USE vast_kind_param,ONLY: DOUBLE
      use meci_C, only : nmeci
      INTEGER, INTENT(IN) :: NMOS
      INTEGER, DIMENSION(NMOS), INTENT(IN) :: IALPHA
      INTEGER, DIMENSION(NMOS), INTENT(IN) :: IBETA
      REAL(DOUBLE), DIMENSION(NMECI), INTENT(IN) :: EIGA
      REAL(DOUBLE), DIMENSION(NMOS,NMOS,NMOS,NMOS), INTENT(IN) :: XY
      END FUNCTION
      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:
module fortapt_m
  !!
  !! A Test Anything Protocol (TAP) producer module inspired by Perl's
  !! Test::More module for planning and setting up tests.
  !!
  use test_base_m,      only: diag
  use test_base_m,      only: diag_unit
  use test_base_m,      only: fail
  use test_base_m,      only: note
  use test_base_m,      only: ok
  use test_base_m,      only: pass
  use test_base_m,      only: test_unit
  use test_base_m,      only: todo

  use test_planning_m,  only: plan
  use test_planning_m,  only: done_testing
  use test_planning_m,  only: skip_all
  use test_planning_m,  only: bail_out

  use test_more_m,      only: is
  use test_more_m,      only: isabs
  use test_more_m,      only: isrel
  use test_more_m,      only: isnear
  use test_more_m,      only: skip
  implicit none
  public
end module fortapt_m

max_stars_repo_name: jbdv-no/fortapt

max_issues_repo_name: jbdv-no/fortapt

max_forks_repo_name: jbdv-no/fortapt


text:
! Tabulated SFH defined in a file called sfh.dat that must reside in the data
! directory (or a file specified via the parameter sfh filename; see
! below). The file must contain three columns. The first column is time since
! the Big Bang in Gyr, the second is the SFR in units of solar masses per year,
! the third is the absolute metallicity. An example is provided in the data
! directory. The time grid in this file can be arbitrary (so long as the units
! are correct), but it is up to the user to ensure that the tabulated sfh is
! well-sampled so that the outputs are stable. Obviously, highly oscillatory
! data require dense sampling.
!
! SFRs are clipped to a minimum value of 1e-30 to avoid divide by zero errors.

subroutine setup_tabular_sfh(pset, nzin)

  use sps_vars, only: sfh_tab, ntabsfh, ntabmax, nz, &
                      tiny_number, tiny30, PARAMS, SPS_HOME
  implicit none
  type(PARAMS), intent(in) :: pset
  integer, intent(in) :: nzin
  integer :: stat, n

  IF (pset%sfh.EQ.2) THEN

     if (pset%sf_start.gt.tiny_number) then
        WRITE(*,*) 'COMPSP ERROR: Tabular sfh, but sf_start > 0'
        STOP
     endif

     ! Read the sfh file
     IF (TRIM(pset%sfh_filename).EQ.'') THEN
        OPEN(3,FILE=TRIM(SPS_HOME)//'/data/sfh.dat',ACTION='READ',STATUS='OLD')
     ELSE
        OPEN(3,FILE=TRIM(SPS_HOME)//'/data/'//TRIM(pset%sfh_filename),&
             ACTION='READ',STATUS='OLD')
     ENDIF
     DO n=1,ntabmax
        IF (nzin.EQ.nz) THEN
           READ(3,*,IOSTAT=stat) sfh_tab(1,n),sfh_tab(2,n),sfh_tab(3,n)
        ELSE
           READ(3,*,IOSTAT=stat) sfh_tab(1,n),sfh_tab(2,n)
           sfh_tab(3,n)=0.0
        ENDIF
        IF (stat.NE.0) GOTO 29
     ENDDO
     WRITE(*,*) 'COMPSP ERROR: didnt finish reading in the sfh file,'
     WRITE(*,*) '     increase ntabmax variable in sps_vars.f90 file'
     STOP
29   CONTINUE
     CLOSE(3)

     ntabsfh = n-1
     sfh_tab(1,1:ntabsfh) = sfh_tab(1,1:ntabsfh)*1E9 !convert to yrs

     !special switch to compute only the last time output
     !in the tabulated file
     ! IF (pset%tage.EQ.-99.) imin=imax

  ELSE IF (pset%sfh.EQ.3) THEN

     !sfh_tab array is supposed to already be filled in, check that it is
     IF (ntabsfh.EQ.0) THEN
        WRITE(*,*) 'COMPSP ERROR: sfh=3 but sfh_tab array not initialized!'
        STOP
     ENDIF

  ENDIF

  ! clip SFR to a minimum of 1e-30
  do n=1, ntabsfh
     sfh_tab(2, n) = max(sfh_tab(2, n), tiny30)
  enddo

end subroutine setup_tabular_sfh

max_stars_repo_name: christopherlovell/fsps

max_issues_repo_name: christopherlovell/fsps

max_forks_repo_name: christopherlovell/fsps


! { dg-do run }
! Tests the fix for PR42104 in which the call to the procedure pointer
! component caused an ICE because the "always_implicit flag was not used
! to force the passing of a descriptor for the array argument.
!
! Contributed by Martien Hulsen <m.a.hulsen@tue.nl>
!
module poisson_functions_m

  implicit none

contains

  function func ( nr, x )
    integer, intent(in) :: nr
    real, intent(in), dimension(:) :: x
    real :: func

    real :: pi

    pi = 4 * atan(1.)

    select case(nr)
      case(1)
        func = 0
      case(2)
        func = 1
      case(3)
        func = 1 + cos(pi*x(1))*cos(pi*x(2))
      case default
        write(*,'(/a,i0/)') 'Error func: wrong function number: ', nr
        stop
    end select

  end function func

end module poisson_functions_m

module element_defs_m

  implicit none

  abstract interface
    function dummyfunc ( nr, x )
      integer, intent(in) :: nr
      real, intent(in), dimension(:) :: x
      real :: dummyfunc
    end function dummyfunc
  end interface

  type function_p
    procedure(dummyfunc), nopass, pointer :: p => null()
  end type function_p

end module element_defs_m

program t

use poisson_functions_m
use element_defs_m

procedure(dummyfunc), pointer :: p => null()
type(function_p) :: funcp

p => func
funcp%p => func

print *, func(nr=3,x=(/0.1,0.1/))
print *, p(nr=3,x=(/0.1,0.1/))
print *, funcp%p(nr=3,x=(/0.1,0.1/))

end program t


text:
      SUBROUTINE CROTG(CA,CB,C,S)
*     .. Scalar Arguments ..
      COMPLEX CA,CB,S
      REAL C
*     ..
*
*  Purpose
*  =======
*
*  CROTG determines a complex Givens rotation.
*
*  =====================================================================
*
*     .. Local Scalars ..
      COMPLEX ALPHA
      REAL NORM,SCALE
*     ..
*     .. Intrinsic Functions ..
      INTRINSIC CABS,CONJG,SQRT
*     ..
      IF (CABS(CA).EQ.0.) THEN
         C = 0.
         S = (1.,0.)
         CA = CB
      ELSE
         SCALE = CABS(CA) + CABS(CB)
         NORM = SCALE*SQRT((CABS(CA/SCALE))**2+ (CABS(CB/SCALE))**2)
         ALPHA = CA/CABS(CA)
         C = CABS(CA)/NORM
         S = ALPHA*CONJG(CB)/NORM
         CA = ALPHA*NORM
      END IF
      RETURN
      END

max_stars_repo_name: NetlibRepository/lapack-www

max_issues_repo_name: NetlibRepository/lapack-www

max_forks_repo_name: NetlibRepository/lapack-www


text:
	SUBROUTINE GDXINP  ( gdfile, gdatim, gvcord, cxstns, gfunc,
     +			     cint, scale, line,
     +			     ptype, yaxis, border, gvect,
     +			     wind, refvec, skip, title,  clear,
     +                       device, text, panel, contur, fint, fline,
     +			     ctype, clrbar, ijskip, iret )
C************************************************************************
C* GDXINP								*
C*									*
C* This subroutine gets the input parameters for GDCROSS.		*
C*									*
C* GDXINP  ( GDFILE, GDATIM, GVCORD, CXSTNS, GFUNC, CINT, SCALE, LINE,	*
C*   	     PTYPE,  YAXIS, BORDER, GVECT, WIND, REFVEC, SKIP, TITLE,	*
C*	     CLEAR, DEVICE, TEXT, PANEL, CONTUR, FINT, FLINE,		*
C*	     CTYPE, CLRBAR, IJSKIP, IRET )				*
C*        								*
C**									*
C* Log:									*
C* K. F. BRill/GSC	 5/89   Created from GDPINP			*
C* S. Schotz/GSC	 7/90	Update for IN_LINE			*
C* K. Brill/NMC		01/92	Added CONTUR, FINT, FLINE, CTYPE	*
C* L. Sager/NMC		 7/93	Added REFVEC				*
C* S. Jacobs/EAI	 9/93	Added CLRBAR				*
C* T. Lee/SAIC		 3/02	Returned SKIP with ' '			*
C* K. Brill/HPC		12/02	Added IJSKIP				*
C************************************************************************
	CHARACTER*(*)	gdfile, gdatim, gvcord, cxstns, gfunc, gvect,
     +			     cint, scale, line,
     +			     ptype,  yaxis, border, wind, refvec,
     +			     skip, title, device, text, panel, contur,
     +			     fint, fline, ctype, clrbar, ijskip
	LOGICAL		clear
	INTEGER 	ier(25)
C------------------------------------------------------------------------
	skip = ' '
	CALL IP_STR  ( 'GDFILE',  gdfile,  ier(1) )
	CALL IP_STR  ( 'GDATTIM', gdatim,  ier(2) )
	CALL IP_STR  ( 'GVCORD',  gvcord,  ier(3) )
	CALL IP_STR  ( 'CXSTNS',  cxstns,  ier(4) )
	CALL IP_STR  ( 'GFUNC',   gfunc,   ier(5) )
	CALL IP_STR  ( 'CINT',    cint,    ier(6) )
	CALL IP_STR  ( 'SCALE',   scale,   ier(7) )
	CALL IP_STR  ( 'LINE',    line,    ier(8) )
	CALL IP_STR  ( 'PTYPE',   ptype,   ier(9) )
	CALL IP_STR  ( 'YAXIS',   yaxis,   ier(10) )
	CALL IP_STR  ( 'BORDER',  border,  ier(11) )
	CALL IP_STR  ( 'GVECT' ,  gvect,   ier(12) )
	CALL IP_STR  ( 'WIND',    wind,    ier(13) )
	CALL IP_STR  ( 'REFVEC',  refvec,  ier(14) )
 	CALL IP_STR  ( 'TITLE',   title,   ier(15) )
	CALL IP_LOG  ( 'CLEAR',   clear,   ier(16) )
	CALL IP_STR  ( 'DEVICE',  device,  ier(17) )
	CALL IP_STR  ( 'TEXT',    text,    ier(18) )
	CALL IP_STR  ( 'PANEL',   panel,   ier(19) )
	CALL IP_STR  ( 'CONTUR',  contur,  ier(20) )
	CALL IP_STR  ( 'FINT',    fint,    ier(21) )
	CALL IP_STR  ( 'FLINE',   fline,   ier(22) )
	CALL IP_STR  ( 'CTYPE',   ctype,   ier(23) )
	CALL IP_STR  ( 'CLRBAR',  clrbar,  ier(24) )
	CALL IP_STR  ( 'IJSKIP',  ijskip,  ier(25) )
C*
	iret = 0
	DO i = 1, 25
	  iret = iret + ier (i)
	END DO
C*
	IF  ( iret .ne. 0 )  iret = -2
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      subroutine read_xsect(ntit, titles, s02x, erelax, wpx, edge,
     1                      emu, gamma, ne, ne1, ik0,
     2                      er, ei, xsn, col4, col5)

      use json_module
      implicit double precision (a-h, o-z)

      include '../HEADERS/const.h'
      include '../HEADERS/dim.h'

      logical :: found
      type(json_file) :: json   !the JSON structure read from the file:
      double precision,dimension(:),allocatable :: dbpc1, dbpc2, dbpc3
      double precision,dimension(:),allocatable :: dbpc4, dbpc5
      character*80,dimension(:),allocatable :: strings

      character*80 titles(nheadx)
      double precision er(nex), ei(nex), xsn(nex)
      double precision col4(nex), col5(nex)

      call json%load_file('xsect.json')
      if (json_failed()) then   !if there was an error reading the file
         print *, "failed to read xsect.json"
         stop
      else
         call json%get('ntitle', ntit, found)
                   if (.not. found) call bailout('s02',   'xsect.json')
         call json%get('title',   strings, found)
                   if (.not. found) call bailout('title', 'xsect.json')
         do 10 i=1,ntit
            titles(i) = strings(i)
 10      continue

         call json%get('s02',    s02x, found)
                   if (.not. found) call bailout('s02',    'xsect.json')
         call json%get('erelax', erelax, found)
                   if (.not. found) call bailout('erelax', 'xsect.json')
         call json%get('wp',     wpx, found)
                   if (.not. found) call bailout('wp',     'xsect.json')
         call json%get('edge',   edge, found)
                   if (.not. found) call bailout('edge',   'xsect.json')
         call json%get('emu',    emu, found)
                   if (.not. found) call bailout('emu',    'xsect.json')
         call json%get('gamach', gamma, found)
                   if (.not. found) call bailout('gamach', 'xsect.json')
         call json%get('ne',     ne, found)
                   if (.not. found) call bailout('ne',     'xsect.json')
         call json%get('ne1',    ne1, found)
                   if (.not. found) call bailout('ne1',    'xsect.json')
         call json%get('ik0',    ik0, found)
                   if (.not. found) call bailout('ik0',    'xsect.json')

         call json%get('ereal',   dbpc1, found)
                   if (.not. found) call bailout('er',     'xsect.json')
         call json%get('eimag',   dbpc2, found)
                   if (.not. found) call bailout('ei',     'xsect.json')
         call json%get('xsnorm',  dbpc3, found)
                   if (.not. found) call bailout('xsn',    'xsect.json')
         call json%get('dum1',   dbpc4, found)
                   if (.not. found) call bailout('col4',   'xsect.json')
         call json%get('dum2',   dbpc5, found)
                   if (.not. found) call bailout('col5',   'xsect.json')

         do 20 i=1,ne
            er(i)   = dbpc1(i)
            ei(i) = dbpc2(i)
            xsn(i) = dbpc3(i)
            col4(i) = dbpc4(i)
            col5(i) = dbpc5(i)
 20      continue
         call json%destroy()
      end if

      return
      end

      subroutine read_titles(ntit, titles)
      use json_module
      implicit double precision (a-h, o-z)

      include '../HEADERS/const.h'
      include '../HEADERS/dim.h'

      logical :: found
      type(json_file) :: json   !the JSON structure read from the file:
      character*80,dimension(:),allocatable :: strings

      character*80 titles(nheadx)

      call json%load_file('xsect.json')
      if (json_failed()) then   !if there was an error reading the file
         print *, "failed to read xsect.json"
         stop
      else
         call json%get('ntitle', ntit, found)
                   if (.not. found) call bailout('s02',   'xsect.json')
         call json%get('title',   strings, found)
                   if (.not. found) call bailout('title', 'xsect.json')
         do 10 i=1,ntit
            titles(i) = strings(i)
 10      continue
         call json%destroy()
      end if

      return
      end


! { dg-do run }
!
! Ensure that move_alloc for CLASS resets the FROM variable's dynamic type
! to the declared one
!
implicit none
type t
end type t
type, extends(t) :: t2
end type t2

class(t), allocatable :: a, b, c
class(t), allocatable :: a2(:), b2(:), c2(:)
allocate (t2 :: a)
allocate (t2 :: a2(5))
call move_alloc (from=a, to=b)
call move_alloc (from=a2, to=b2)
!print *, same_type_as (a,c), same_type_as (a,b)
!print *, same_type_as (a2,c2), same_type_as (a2,b2)
if (.not. same_type_as (a,c) .or. same_type_as (a,b)) STOP 1
if (.not. same_type_as (a2,c2) .or. same_type_as (a2,b2)) STOP 2
end


text:
	SUBROUTINE MA_IFSP  ( rimnem, cimnem, iret )
C************************************************************************
C* MA_IFSP								*
C*									*
C* This subroutine initializes the interface mnemonic arrays and sets   *
C* the pointers within COMMON / RINTFP / and COMMON / CINTFP /.		*
C*									*
C* MA_IFSP  ( RIMNEM, CIMNEM, IRET )					*
C*									*
C* Output parameters:							*
C*	RIMNEM (*)	CHAR*		Interface mnemonics for reals   *
C*	CIMNEM (*)	CHAR*		Interface mnemonics for chars   *
C*	IRET		INTEGER		Return code 			*
C*					  0 = normal return 		*
C*					 -1 = One or more pointers	*
C*					      could not be set 		*
C*									*
C**									*
C* Log:									*
C* D. Kidwell/NCEP	10/97	Adapted from AF_IFSP			*
C* R. Hollern/NCEP	12/97	Replaced MA_IFFP, MA_IFMP calls with DC_*
C* D. Kidwell/NCEP	10/98	Added init of interface mnemonics;      *
C*				added intf mnemonics to calling sequence*
C* R. Hollern/NCEP	 7/99	Added 24 more mnemonics                 *
C* R. Hollern/NCEP	 9/99	Added drifting buoy section 4 mnemonics *
C* C. Caruso Magee/NCEP	 4/01	Modified for Coast Guard reports	*
C* D. Kidwell/NCEP	 3/02	Corrected PRMN to PSMN                  *
C* D. Kidwell/NCEP	 4/05	Added PKWD, PKWS, PKWT, PMN1, PMNT      *
C************************************************************************
	INCLUDE		'macmn.cmn'
C*
	CHARACTER*(*)	rimnem (*), cimnem (*)
C*
	LOGICAL		allok
	INTEGER		iploc (NRSLMN)
	CHARACTER	rifmn ( NRIMN )*8, cifmn ( NCIMN )*8
C*
C*	Establish array equivalence to COMMON / RINTFP /
C*
	EQUIVALENCE	( iploc (1), iryear )
C*
	DATA	( rifmn (i), i = 1, 58 ) /
     +		'YEAR'  , 'MNTH'  , 'DAYS'  , 'HOUR'  , 'MINU'  ,
     +		'SLAT'  , 'SLON'  , 'SELV'  ,
     +		'DRCT'  , 'SPED'  , 'GUM0'  , 'GUM2'  , 'GUM1'  ,
     +		'XS10'  , 'XS20'  , 'PRES'  , 'PMSL'  , 'P03D'  ,
     +		'CHPT'  , '3HPC'  , 'RELH'  , 'MSST'  , 'STMWC' ,
     +		'TMPC'  , 'DWPC'  , 'SSTC'  , 'TMWC'  , 'DTFV1' ,
     +		'MXTM'  , 'DTFV2' , 'MITM'  , 'QOPM'  , 'QCBH'  ,
     +		'QWTM'  , 'QATM'  , 'QBST'  , 'QCIL'  , 'MSDM'  ,
     +		'PXXM'  , 'P06M'  , 'P12M'  , 'P18M'  , 'P24M'  ,
     +		'P01M'  , 'P02M'  , 'P03M'  , 'P09M'  , 'P15M'  ,
     +		'INPC'  , 'ITSO'  , 'TOST'  , 'CBAS'  , 'VSBK'  ,
     +		'CFRT'  , 'WWMO'  , 'WWMA'  , 'PWWM'  , 'PWWA'  /
C*
	DATA	( rifmn (i), i = 59, NRSLMN ) /
     +		'PSW2'  , 'CFRL'  , 'CTYL'  , 'CLHL'  , 'CTYM'  ,
     +		'CLHM'  , 'CTYH'  , 'CLHH'  , 'TDMP'  , 'ASMP'  ,
     +		'WPER'  , 'POWW'  , 'HOWW'  , 'COIA'  , 'IDTH'  ,
     + 		'ROIA'  , '24PC'  , 'PWYR'  , 'PWMO'  , 'PWDY'  ,
     +		'PWHR'  , 'PWMN'  , 'PWDR'  , 'PWSP'  , 'DTFVM' ,
     +		'GUMS'  , 'WHGT'  , 'CORN'  , 'BPID'  , 'SUWS'  ,
     +		'QPOS'  , 'QTIM'  ,
     +		'QCLS'  , 'QDS1'  , 'QXS1'  , 'QDS2'  , 'QXS2'  ,
     +		'Q3D1'  , 'Q3D2'  , 'Q4CL'  , 'PWA2'  , 'VRTM'  ,
     +		'PSYR'  , 'PSMN'  , 'PSDY'  , 'PSHR'  , 'PSMI'  ,
     +		'DBVV'  , 'DBDD'  , 'BENG'  , 'DROT'  , 'CALN'  ,
     +		'DLATH' , 'DLONH' , 'SKNT'  , 'ALTI'  , 'TMPF'  ,
     +		'SSTF'  , 'VSBY'  , 'MXWH'  , 'TERC'  , 'GUST'  ,
     +		'PKWD'  , 'PKWS'  , 'PKWT'  , 'PMN1'  , 'PMNT'  /
C*
	DATA	( rifmn (i), i = NRSLMN + 1, NRIMN ) /
     +		'NDTS'  , 'DBSS'  , 'STMP'  , 'SALN'  ,
     +                    MDM1T3 * ' ',
     +          'NCLO'  , 'VSSO'  , 'CLAM'  , 'CLTP'  , 'HOCB'  ,
     +		          MCM1T4 * ' ' ,
     +		'NCL3'  , 'VSSO3' , 'CLAM3' , 'CLTP3' , 'HOCB3' ,
     +                    M3M1T4 * ' ',
     +          'NCWD'  , 'TPMI'  , 'WDRC'  , 'WDSC'  ,
     +		          MWM1T3 * ' ',
     +          'NSWV'  , 'DOSW'  , 'POSW'  , 'HOSW'  ,
     +		          MSM1T3 * ' '                          /
C*					Real interface mnemonics
C*
	DATA	( cifmn (i), i = 1, NCIMN ) 	/ 'STID' /
C*					Character interface mnemonics
C-----------------------------------------------------------------------
	iret = 0
C
C*	Initialize the interface mnemonics.
C
	DO i = 1, NRIMN
	    rimnem ( i ) = rifmn ( i )
	END DO
	DO i = 1, NCIMN
	    cimnem ( i ) = cifmn ( i )
	END DO
C
C*	The logical variable "allok" is initially set to .true. but will
C*	be reset to .false. if any of the pointers cannot be set.
C
	allok = .true.
C
C*	Set the pointers for the single-level interface mnemonics.
C
	DO i = 1, NRSLMN
	    CALL DC_IFFP ( rimnem ( i ), rimnem, NRIMN, allok,
     +                     iploc ( i ), ier )
	END DO
C
	CALL DC_IFFP  ( 'STID', cimnem, NCIMN, allok, icstid, ier )
C
C*	Set the pointers for depth/temperature/salinity data.
C
	CALL DC_IFFP  ( 'NDTS', rimnem, NRIMN, allok, irndts, ier )
	CALL DC_IFFP  ( 'DBSS', rimnem, NRIMN, allok, irdbss (1), ier )
	CALL DC_IFMP  ( 3, MXDLYR, irdbss, ier )
	CALL DC_IFFP  ( 'STMP', rimnem, NRIMN, allok, irstmp (1), ier )
	CALL DC_IFMP  ( 3, MXDLYR, irstmp, ier )
	CALL DC_IFFP  ( 'SALN', rimnem, NRIMN, allok, irsaln (1), ier )
	CALL DC_IFMP  ( 3, MXDLYR, irsaln, ier )
C
C*	Set the pointers for cloud data.
C
	CALL DC_IFFP  ( 'NCLO', rimnem, NRIMN, allok, irnclo, ier )
	CALL DC_IFFP  ( 'VSSO', rimnem, NRIMN, allok, irvsso (1), ier )
	CALL DC_IFMP  ( 4, MXCLYR, irvsso, ier )
	CALL DC_IFFP  ( 'CLAM', rimnem, NRIMN, allok, irclam (1), ier )
	CALL DC_IFMP  ( 4, MXCLYR, irclam, ier )
	CALL DC_IFFP  ( 'CLTP', rimnem, NRIMN, allok, ircltp (1), ier )
	CALL DC_IFMP  ( 4, MXCLYR, ircltp, ier )
	CALL DC_IFFP  ( 'HOCB', rimnem, NRIMN, allok, irhocb (1), ier )
	CALL DC_IFMP  ( 4, MXCLYR, irhocb, ier )
C
C*	Set the pointers for cloud3 data.
C
	CALL DC_IFFP  ( 'NCL3', rimnem, NRIMN, allok, irncl3, ier )
	CALL DC_IFFP  ( 'VSSO3', rimnem, NRIMN, allok, irvss3(1), ier )
	CALL DC_IFMP  ( 4, MX3LYR, irvss3, ier )
	CALL DC_IFFP  ( 'CLAM3', rimnem, NRIMN, allok, ircla3(1), ier )
	CALL DC_IFMP  ( 4, MX3LYR, ircla3, ier )
	CALL DC_IFFP  ( 'CLTP3', rimnem, NRIMN, allok, irclt3(1), ier )
	CALL DC_IFMP  ( 4, MX3LYR, irclt3, ier )
	CALL DC_IFFP  ( 'HOCB3', rimnem, NRIMN, allok, irhcb3(1), ier )
	CALL DC_IFMP  ( 4, MX3LYR, irhcb3, ier )
C
C*	Set the pointers for continuous wind data.
C
	CALL DC_IFFP  ( 'NCWD', rimnem, NRIMN, allok, irncwd, ier )
	CALL DC_IFFP  ( 'TPMI', rimnem, NRIMN, allok, irtpmi (1), ier )
	CALL DC_IFMP  ( 3, MXWLYR, irtpmi, ier )
	CALL DC_IFFP  ( 'WDRC', rimnem, NRIMN, allok, irwdrc (1), ier )
	CALL DC_IFMP  ( 3, MXWLYR, irwdrc, ier )
	CALL DC_IFFP  ( 'WDSC', rimnem, NRIMN, allok, irwdsc (1), ier )
	CALL DC_IFMP  ( 3, MXWLYR, irwdsc, ier )
C
C*	Set the pointers for wave data.
C
	CALL DC_IFFP  ( 'NSWV', rimnem, NRIMN, allok, irnswv, ier )
	CALL DC_IFFP  ( 'DOSW', rimnem, NRIMN, allok, irdosw (1), ier )
	CALL DC_IFMP  ( 3, MXSLYR, irdosw, ier )
	CALL DC_IFFP  ( 'POSW', rimnem, NRIMN, allok, irposw (1), ier )
	CALL DC_IFMP  ( 3, MXSLYR, irposw, ier )
	CALL DC_IFFP  ( 'HOSW', rimnem, NRIMN, allok, irhosw (1), ier )
	CALL DC_IFMP  ( 3, MXSLYR, irhosw, ier )
C
C*	Were all of the pointers properly set?
C
	IF  ( .not. allok ) THEN
	    iret = -1
	END IF
C*
	RETURN
	END