text:
! This module stores the runtime parameters.  The probin_gmres_init() routine is
! used to initialize the runtime parameters

module probin_gmres_module

  use bl_types
  use bl_space

  implicit none

  ! For comments and instructions on how to set the input parameters see namelist section below
  !-------------------------------------------------------------
  integer   , save :: precon_type,visc_schur_approx
  integer   , save :: mg_verbose,cg_verbose,mg_max_vcycles,mg_minwidth
  integer   , save :: mg_bottom_solver,mg_nsmooths_down,mg_nsmooths_up,mg_nsmooths_bottom
  integer   , save :: mg_max_bottom_nlevels,stag_mg_verbosity,stag_mg_max_vcycles
  integer   , save :: stag_mg_minwidth,stag_mg_bottom_solver,stag_mg_nsmooths_down
  integer   , save :: stag_mg_nsmooths_up,stag_mg_nsmooths_bottom,stag_mg_smoother
  integer   , save :: stag_mg_max_bottom_nlevels,gmres_verbose,gmres_max_outer,gmres_max_inner
  integer   , save :: gmres_max_iter,gmres_min_iter
  real(dp_t), save :: p_norm_weight,scale_factor,mg_rel_tol,stag_mg_omega,stag_mg_rel_tol
  real(dp_t), save :: gmres_rel_tol,gmres_abs_tol

  !-------------------------------------------------------------
  ! Input parameters controlled via namelist input, with comments
  !-------------------------------------------------------------

  ! preconditioner type
  ! 1 = projection preconditioner
  !-1 = projection preconditioner with expensive pressure update
  ! 2 = lower triangular preconditioner
  !-2 = lower triangular preconditioner with negative sign
  ! 3 = upper triangular preconditioner
  !-3 = upper triangular preconditioner with negative sign
  ! 4 = Block diagonal preconditioner
  !-4 = Block diagonal preconditioner with negative sign
  namelist /probin_gmres/ precon_type

  ! use the viscosity-based BFBt Schur complement (from Georg Stadler)
  namelist /probin_gmres/ visc_schur_approx

  ! weighting of pressure when computing norms and inner products
  namelist /probin_gmres/ p_norm_weight

  ! scale theta_alpha, beta, gamma, and b_u by this, and then scale x_p by the inverse
  namelist /probin_gmres/ scale_factor

  ! MAC projection solver parameters:
  namelist /probin_gmres/ mg_verbose            ! multigrid verbosity
  namelist /probin_gmres/ cg_verbose            ! BiCGStab (mg_bottom_solver=1) verbosity
  namelist /probin_gmres/ mg_max_vcycles        ! maximum number of V-cycles
  namelist /probin_gmres/ mg_minwidth           ! length of box at coarsest multigrid level
  namelist /probin_gmres/ mg_bottom_solver      ! bottom solver type
                                                ! 0 = smooths only, controlled by mg_nsmooths_bottom
                                                ! 1 = BiCGStab
                                                ! 4 = Fancy bottom solve that coarsens down additionally
                                                !     and then applies mg_nsmooths_bottom smooths
  namelist /probin_gmres/ mg_nsmooths_down      ! number of smooths at each level on the way down
  namelist /probin_gmres/ mg_nsmooths_up        ! number of smooths at each level on the way up
  namelist /probin_gmres/ mg_nsmooths_bottom    ! number of smooths at the bottom (only if mg_bottom_solver=0)
  namelist /probin_gmres/ mg_max_bottom_nlevels ! for mg_bottom_solver=4, number of additional levels of multigrid
  namelist /probin_gmres/ mg_rel_tol            ! relative tolerance stopping criteria

  ! Staggered multigrid solver parameters
  namelist /probin_gmres/ stag_mg_verbosity       ! verbosity
  namelist /probin_gmres/ stag_mg_max_vcycles     ! max number of v-cycles
  namelist /probin_gmres/ stag_mg_minwidth        ! length of box at coarsest multigrid level
  namelist /probin_gmres/ stag_mg_bottom_solver   ! bottom solver type
                                                  ! 0 = smooths only, controlled by mg_nsmooths_bottom
                                                  ! 4 = Fancy bottom solve that coarsens additionally
                                                  !     and then applies stag_mg_nsmooths_bottom smooths
  namelist /probin_gmres/ stag_mg_nsmooths_down   ! number of smooths at each level on the way down
  namelist /probin_gmres/ stag_mg_nsmooths_up     ! number of smooths at each level on the way up
  namelist /probin_gmres/ stag_mg_nsmooths_bottom ! number of smooths at the bottom
  namelist /probin_gmres/ stag_mg_max_bottom_nlevels ! for stag_mg_bottom_solver=4, number of additional levels of multigrid
  namelist /probin_gmres/ stag_mg_omega           ! weighted-jacobi omega coefficient
  namelist /probin_gmres/ stag_mg_smoother        ! 0 = jacobi; 1 = 2*dm-color Gauss-Seidel
  namelist /probin_gmres/ stag_mg_rel_tol         ! relative tolerance stopping criteria

  ! GMRES solver parameters
  namelist /probin_gmres/ gmres_rel_tol         ! relative tolerance stopping criteria
  namelist /probin_gmres/ gmres_abs_tol         ! absolute tolerance stopping criteria
  namelist /probin_gmres/ gmres_verbose         ! gmres verbosity; if greater than 1, more residuals will be printed out
  namelist /probin_gmres/ gmres_max_outer       ! max number of outer iterations
  namelist /probin_gmres/ gmres_max_inner       ! max number of inner iterations, or restart number
  namelist /probin_gmres/ gmres_max_iter        ! max number of gmres iterations
  namelist /probin_gmres/ gmres_min_iter        ! min number of gmres iterations
  !-------------------------------------------------------------

contains

  subroutine probin_gmres_init()

    use f2kcli
    use parallel
    use bc_module
    use bl_IO_module
    use bl_prof_module
    use bl_error_module
    use bl_constants_module
    use cluster_module

    integer    :: narg, farg

    character(len=128) :: fname
    character(len=128) :: probin_env

    logical :: lexist
    logical :: need_inputs

    integer :: un, ierr

    narg = command_argument_count()

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! initialize the runtime parameters
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!

    ! Defaults
    precon_type = 1

    visc_schur_approx = 0

    p_norm_weight = 1.d0

    scale_factor = 1.d0

    mg_verbose = 0
    cg_verbose = 0
    mg_max_vcycles = 1
    mg_minwidth = 2
    mg_bottom_solver = 4
    mg_nsmooths_down = 2
    mg_nsmooths_up = 2
    mg_nsmooths_bottom = 8
    mg_max_bottom_nlevels = 10
    mg_rel_tol = 1.d-9

    stag_mg_verbosity = 0
    stag_mg_max_vcycles = 1
    stag_mg_minwidth = 2
    stag_mg_bottom_solver = 4
    stag_mg_nsmooths_down = 2
    stag_mg_nsmooths_up = 2
    stag_mg_nsmooths_bottom = 8
    stag_mg_max_bottom_nlevels = 10
    stag_mg_omega = 1.d0
    stag_mg_smoother = 1
    stag_mg_rel_tol = 1.d-9

    gmres_rel_tol = 1.d-9
    gmres_abs_tol = 0.d0
    gmres_verbose = 1
    gmres_max_outer = 20
    gmres_max_inner = 5
    gmres_max_iter = 100
    gmres_min_iter = 1

    need_inputs = .true.

    farg = 1
    if ( need_inputs .AND. narg >= 1 ) then
       call get_command_argument(farg, value = fname)
       inquire(file = fname, exist = lexist )
       if ( lexist ) then
          farg = farg + 1
          un = unit_new()
          open(unit=un, file = fname, status = 'old', action = 'read')
          read(unit=un, nml = probin_gmres)
          close(unit=un)
          need_inputs = .false.
       end if
    end if

    ! stuff that can be read in from the command line by appending, e.g., "--precon_type 1"
    do while ( farg <= narg )
       call get_command_argument(farg, value = fname)
       select case (fname)

       case ('--precon_type')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) precon_type

       case ('--visc_schur_approx')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) visc_schur_approx

       case ('--p_norm_weight')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) p_norm_weight

       case ('--scale_factor')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) scale_factor

       case ('--mg_verbose')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_verbose

       case ('--cg_verbose')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) cg_verbose

       case ('--mg_max_vcycles')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_max_vcycles

       case ('--mg_minwidth')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_minwidth

       case ('--mg_bottom_solver')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_bottom_solver

       case ('--mg_nsmooths_down')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_nsmooths_down

       case ('--mg_nsmooths_up')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_nsmooths_up

       case ('--mg_nsmooths_bottom')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_nsmooths_bottom

       case ('--mg_max_bottom_nlevels')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_max_bottom_nlevels

       case ('--mg_rel_tol')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) mg_rel_tol

       case ('--stag_mg_verbosity')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_verbosity

       case ('--stag_mg_max_vcycles')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_max_vcycles

       case ('--stag_mg_minwidth')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_minwidth

       case ('--stag_mg_bottom_solver')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_bottom_solver

       case ('--stag_mg_nsmooths_down')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_nsmooths_down

       case ('--stag_mg_nsmooths_up')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_nsmooths_up

       case ('--stag_mg_nsmooths_bottom')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_nsmooths_bottom

       case ('--stag_mg_max_bottom_nlevels')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_max_bottom_nlevels

       case ('--stag_mg_omega')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_omega

       case ('--stag_mg_smoother')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) stag_mg_smoother

       case ('--gmres_rel_tol')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_rel_tol

       case ('--gmres_abs_tol')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_rel_tol

       case ('--gmres_verbose')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_verbose

       case ('--gmres_max_outer')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_max_outer

       case ('--gmres_max_inner')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_max_inner

       case ('--gmres_max_iter')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_max_iter

      case ('--gmres_min_iter')
          farg = farg + 1
          call get_command_argument(farg, value = fname)
          read(fname, *) gmres_min_iter

       case ('--')
          farg = farg + 1
          exit

       case default

       end select

       farg = farg + 1
    end do

  end subroutine probin_gmres_init

end module probin_gmres_module

max_stars_repo_name: BoxLib-Codes/LowMachMixtures

max_issues_repo_name: BoxLib-Codes/LowMachMixtures

max_forks_repo_name: BoxLib-Codes/LowMachMixtures


text:
!=======================================================================
! Subroutines in this file have been adapted or rewritten
! for the PlanetSimulator of University Hamburg
! by Hui Wan (Max Planck Institute for Meteorology, 2009).
!=======================================================================
!
!****6***0*********0*********0*********0*********0*********0**********72
      subroutine tpcore(q,ps1,ps2,u,v,                &
                        dtoa,dtdx,                    &
                        iord,jord,kord,               &
                        nc,im,jm,nl,dap,dbk,          &
                        iml,j1,j2,js0,jn0,            &
                        cose,cosp,acosp,dlat,rcap,    &
                        cnst,deform,zcross,           &
                        fill,mfct,debug)
!****6***0*********0*********0*********0*********0*********0**********72
!
! Purpose: perform the transport of  3-D mixing ratio fields using
!          externally specified winds on the hybrid Eta-coordinate.
!          One call to tpcore updates the 3-D mixing ratio
!          fields for one time step. [vertical mass flux is computed
!          internally using a center differenced hydrostatic mass
!          continuity equation].
!
! Schemes: Multi-dimensional Flux Form Semi-Lagrangian (FFSL) schemes
!          (Lin and Rood 1996, MWR) with a modified MFCT option
!          (Zalesak 1979).
!
! Multitasking version: 6m
! History: Mriginal vertion by S.-J. Lin, Oct. 1, 1997
!          Modified for the PlanetSimulator by Hui Wan, May 2009
!
! Author of the algorithm: S.-J. Lin
! Address:
!                 Code 910.3, NASA/GSFC, Greenbelt, MD 20771
!                 Phone: 301-286-9540
!                 E-mail: lin@dao.gsfc.nasa.gov
!
! Affiliation:
!                 Joint Center for Earth Systems Technology
!                 The University of Maryland Baltimore County
!                 and  NASA - Goddard Space Flight Center
!
! The algorithm is based on the following papers:
!
! 1. Lin, S.-J., and R. B. Rood, 1996: Multidimensional flux form semi-
!    Lagrangian transport schemes. Mon. Wea. Rev., 124, 2046-2070.
!
! 2. Lin, S.-J., W. C. Chao, Y. C. Sud, and G. K. Walker, 1994: A class of
!    the van Leer-type transport schemes and its applications to the moist-
!    ure transport in a General Circulation Model. Mon. Wea. Rev., 122,
!    1575-1593.
!
! ======
! INPUT:
! ======
! Q(IM,JM,NL,NC): mixing ratios at current time (t)
! NC: total # of constituents
! IM: first (E-W) dimension; # of Grid intervals in E-W is IM
! JM: 2nd (N-S) dimension;   # of Gaussian latitudes is JM
! NL: 3rd dimension (# of layers); vertical index increases from 1 at
!       the model top to NL near the surface (see fig. below).
!       It is assumed that NL > 5.
!
! PS1(IM,JM): surface pressure at current time (t)
! PS2(IM,JM): surface pressure at next time level (t+dt)
!             predicted by the host model
! PS2 is replaced by the PS predicted in this routine (at t+dt) on output.
! Note: surface pressure can have any unit or can be multiplied by any
!       const.
!
! The hybrid ETA-coordinate:
!
! pressure at layer edges are defined as follows:
!
!        p(i,j,k) = AP(k)*PT  +  BP(k)*PS(i,j)          (1)
!
! Where PT is a constant having the same unit as PS.
! AP and BP are unitless constants given at layer edges.
! In all cases  BP(1) = 0., BP(NL+1) = 1.
! The pressure at the model top is PTOP = AP(1)*PT
!
! Assume the upper and lower boundaries of vertical layer k are
! layer edges k and k+1, respectively (see the sketch below).
! The layer thickness at grid point (i,j,k) reads
!
!   delp(i,j,k) = dap(k) + dbk(k)*PS(i,j)
!
! where    dap(k) = ( AP(k+1) - AP(k) )*PT
!          dbk(k) =   BP(k+1) - BP(k)
!
! *********************
! For pure sigma system
! *********************
! AP(k) = 1 for all k, PT = PTOP,
! BP(k) = sige(k) (sigma at edges), PS = Psfc - PTOP, where Psfc
! is the true surface pressure.
! In this implementation, we calculate the values of dap and dbk in
! the initialization phase of a model run, and pass them into this
! subroutine. Since the PlanetSimulator uses a pure sigma coordinate
! with PT = 0, we have
!
!         dap(:) = 0
!         dbk(:) = dsigma(:)
!
!
!
!                  /////////////////////////////////
!              / \ ------ Model top P=PTOP ---------  AP(1), BP(1)
!               |
!    delp(1)    |  ........... Q(i,j,1) ............
!               |
!     W(k=1)   \ / ---------------------------------  AP(2), BP(2)
!
!
!
!     W(k-1)   / \ ---------------------------------  AP(k), BP(k)
!               |
!    delp(K)    |  ........... Q(i,j,k) ............
!               |
!      W(k)    \ / ---------------------------------  AP(k+1), BP(k+1)
!
!
!
!              / \ ---------------------------------  AP(NL), BP(NL)
!               |
!    delp(NL)   |  ........... Q(i,j,NL) .........
!               |
!     W(NL)=0  \ / -----Earth's surface P=Psfc ------ AP(NL+1), BP(NL+1)
!                 //////////////////////////////////
!
! U(IM,JM,NL) & V(IM,JM,NL):winds (m/s) at mid-time-level (t+dt/2)
! Note that on return U and V are destroyed.
!
! dtoa(real): time step (in seconds) divided by the planet radius (in meters).
!      Suggested value for time step: 30 min. for 4x5, 15 min. for 2x2.5
!      (Lat-Lon) resolution. Smaller values maybe needed if the model
!      has a well-resolved stratosphere and Max(V) > 225 m/s
!
! dtdx: real array of shape (JM). time step (in seconds) divided by the
!       zonal grid size (in meters). The values are initialized in subroutine
!       tracer_ini.
!
! J1, J2 are the starting and ending indices of the Gaussian latitudes
!        outside the polar caps. Note that in the tracer transport routines
!        we count from south to north.
!        The South polar cap edge is located between the first and second
!        Gaussian latitudes from the south; The North polar cap edge is
!        located between the last two Gaussian latitudes.
!
! js0, jn0: [1, js0] and [jn0, JM] are the latitude ranges in which
!           the semi-Lagrangian method may be needed for transport in
!           the x-direction. Their values are initialized in subroutine
!           tracer_ini.
!
! cose,cosp,acosp,dlat,rcap: these parameters are related to the
!        locations of the Gaussiang latitudes. They are initialized in
!        subroutine tracer_ini.
!
! IORD, JORD, and KORD are integers controlling various options in E-W, N-S,
! and vertical transport, respectively.
!
!
!  _ORD=
!     1: 1st order upstream scheme (too diffusive, not a real option; it
!        can be used for debugging purposes; this is THE only known "linear"
!        monotonic advection scheme.).
!     2: 2nd order van Leer (full monotonicity constraint;
!        see Lin et al 1994, MWR)
!     3: monotonic PPM* (Collela & Woodward 1984)
!     4: semi-monotonic PPM (same as 3, but overshoots are allowed)
!     5: positive-definite PPM (constraint on the subgrid distribution is
!        only strong enough to prevent generation of negative values;
!        both overshoots & undershootes are possible).
!     6: un-constrained PPM (nearly diffusion free; faster but
!        positivity of the subgrid distribution is not quaranteed. Use
!        this option only when the fields and winds are very smooth or
!        when MFCT=.true.)
!
! *PPM: Piece-wise Parabolic Method
!
! Recommended values:
! IORD=JORD=3 for high horizontal resolution.
! KORD=6       if MFCT=.true.
! KORD=3 or 5  if MFCT=.false.
!
! The implicit numerical diffusion decreases as _ORD increases.
! DO not use option 4 or 5 for non-positive definite scalars
! (such as Ertel Potential Vorticity).
!
! If numerical diffusion is a problem (particularly at low horizontal
! resolution) then the following setup is recommended:
! IORD=JORD=KORD=6 and MFCT=.true.
!
! FILL (logical):   flag to do filling for negatives (see note below).
! MFCT (logical):   flag to do a Zalesak-type Multidimensional Flux
!                   correction. It shouldn't be necessary to call the
!                   filling routine when MFCT is true.
!
! ======
! Output
! ======
!
! Q: the updated mixing ratios at t+NDT (original values are over-written)
! W(;;NL): large-scale vertical mass flux as diagnosed from the hydrostatic
!          relationship. W will have the same unit as PS1 and PS2 (eg, mb).
!          W must be divided by dt to get the correct mass-flux unit.
!          The vertical Courant number C = W/delp_UPWIND, where delp_UPWIND
!          is the pressure thickness in the "upwind" direction. For example,
!          C(k) = W(k)/delp(k)   if W(k) > 0;
!          C(k) = W(k)/delp(k+1) if W(k) < 0.
!              ( W > 0 is downward, ie, toward surface)
! PS2: predicted PS at t+dt (original values are over-written)
!
! =====
! NOTES:
! =====
!
! This forward-in-time upstream-biased transport scheme degenerates to
! the 2nd order center-in-time center-in-space mass continuity eqn.
! if Q = 1 (constant fields will remain constant). This degeneracy ensures
! that the computed vertical velocity to be identical to GEOS-1 GCM
! for on-line transport.
!
! A larger polar cap is used if j1=3 (recommended for C-Grid winds or when
! winds are noisy near poles).
!
! PPM is 4th order accurate when grid spacing is uniform (x & y); 3rd
! order accurate for non-uniform grid (vertical sigma coord.).
!
! Time step is limitted only by transport in the meridional direction.
! (the FFSL scheme is not implemented in the meridional direction).
!
! Since only 1-D limiters are applied, negative values could
! potentially be generated when large time step is used and when the
! initial fields contain discontinuities.
! This does not necessarily imply the integration is unstable.
! These negatives are typically very small. A filling algorithm is
! activated if the user set "fill" to be true.
! Alternatively, one can use the MFCT option to enforce monotonicity.
!
      implicit none

! Input-Output arrays

      real,intent(inout) ::   q(im,jm,nl,nc)
      real,intent(inout) :: ps1(im,jm)
      real,intent(inout) :: ps2(im,jm)
      real,intent(inout) ::   u(im,jm,nl)
      real,intent(inout) ::   v(im,jm,nl)
      real,intent(in)    ::  dtdx(jm), dtoa

      integer,intent(in) :: iord,jord,kord,cnst
      integer,intent(in) :: nc,im,jm,nl,iml,j1,j2,js0,jn0
      real,intent(in)    :: cose(jm+1),cosp(jm), acosp(jm), dlat(jm), rcap
      real,intent(in)    :: dap(nl), dbk(nl)
      logical,intent(in) :: zcross, fill, mfct, deform, debug

! Local dynamic arrays

      integer :: js(nl),jn(nl)

      real :: dtdx5(jm)
      real ::   w(im,jm,nl)
      real :: crx(im,jm,nl),cry(im,jm,nl)
      real :: delp(im,jm,nl),delp1(im,jm,nl),delp2(im,jm,nl)
      real :: delp2dyn(im,jm,nl)
      real :: xmass(im,jm,nl),ymass(im,jm,nl)
      real ::   dg1(im),dg2(im,jm),dpi(im,jm,nl)
      real ::  qlow(im,jm,nl), dq(im,jm,nl)
      real ::    qz(im,jm,nl),qmax(im,jm,nl),qmin(im,jm,nl)
      real ::    wk(im,jm,nl),pu(im,jm,nl)
      real ::    fx(im+1,jm,nl),fy(im,jm,nl),fz(im,jm,nl+1)

! scalars

      integer :: imh, i,j,k,jt,ic
      real    :: d5, dtoa5, sum1, sum2

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

      imh      = im/2
      dtdx5(:) = 0.5*dtdx(:)
      dtoa5    = 0.5*dtoa

      if (debug) then
         write(*,*) '* Entered routine tpcore'
         do ic=1,nc
            write(*,*) '* tracer',ic,'* max. & min. q =', &
                       maxval(q(:,:,:,ic)),minval(q(:,:,:,ic))
         end do
      endif

! save the vertical layer thickness (at t+dt) predicted by
! the dynamical core

      do k=1,nl
         delp2dyn(:,:,k) = dap(k) + dbk(k)*ps2(:,:)
      end do

! estimate the value at t+0.5dt

      ps2(:,:) = 0.5*( ps1(:,:) + ps2(:,:) )

!****6***0*********0*********0*********0*********0*********0**********72
! Compute Courant number
!****6***0*********0*********0*********0*********0*********0**********72
! Convert winds on A-Grid to Courant # on C-Grid.

      do j=2,jm-1
      do i=2,im
         crx(i,j,:) = dtdx5(j)*(u(i,j,:)+u(i-1,j,:))
      end do
         crx(1,j,:) = dtdx5(j)*(u(1,j,:)+u(im,j,:))
      end do

      do j=2,jm
         cry(:,j,:) = dtoa5*(v(:,j,:)+v(:,j-1,:))
      enddo

!****6***0*********0*********0*********0*********0*********0**********72
! Find JN and JS
!****6***0*********0*********0*********0*********0*********0**********72
! [2,js(k)] and [jn(k),jm-1] are the latitudes at which semi-Lagrangian
! treatment is needed.

!MIC$ do all autoscope shared(JS,JN,CRX,CRY,PS2,U,V,DPI,ymass,delp2,PU)
!MIC$* shared(xmass)
!MIC$* private(i,j,k,sum1,sum2,D5)

      do k=1,nl

      js(k) = j1
      jn(k) = j2

      do 111 j=js0,j1+1,-1
      do 111 i=1,im
      if(abs(crx(i,j,k)) .gt. 1.) then
            js(k) = j
            go to 112
      endif
111   continue
112   continue

      do 122 j=jn0,j2-1
      do 122 i=1,im
      if(abs(crx(i,j,k)) .gt. 1.) then
            jn(k) = j
            go to 133
      endif
122   continue
133   continue

      enddo  !vertical layer loop

!****6***0*********0*********0*********0*********0*********0**********72
! ***** Compute horizontal mass fluxes *****
!****6***0*********0*********0*********0*********0*********0**********72
! for the polar caps: replace the values at the northest (southest)
! latitude with the zonal mean.

      sum1 = sum(ps1(:,1 ))/im
      sum2 = sum(ps1(:,jm))/im
      ps1(:,1)  = sum1
      ps1(:,jm) = sum2

      sum1 = sum(ps2(:,1 ))/im
      sum2 = sum(ps2(:,jm))/im
      ps2(:,1)  = sum1
      ps2(:,jm) = sum2

      do k=1,nl

! delp = pressure thickness: the psudo-density in a hydrostatic system.

      delp2(:,:,k) = dap(k) + dbk(k)*ps2(:,:)

! N-S componenet

      do j=j1,j2+1
         d5 = 0.5 * cose(j)
         ymass(:,j,k) = cry(:,j,k)*d5*(delp2(:,j,k)+delp2(:,j-1,k))
      enddo

      do j=j1,j2
         dpi(:,j,k) = (ymass(:,j,k)-ymass(:,j+1,k)) *acosp(j)/dlat(j)
      end do

! polar caps

      sum1 = -sum(ymass(:,j1  ,k))*rcap
      sum2 =  sum(ymass(:,j2+1,k))*rcap

      dpi(:, 1,k) = sum1
      dpi(:,jm,k) = sum2

! E-W component

      do j=j1,j2
      do i=2,im
         pu(i,j,k) = 0.5 * (delp2(i,j,k) + delp2(i-1,j,k))
      enddo
      enddo

      do j=j1,j2
         pu(1,j,k) = 0.5 * (delp2(1,j,k) + delp2(im,j,k))
      enddo

      do j=j1,j2
      do i=1,im
         xmass(i,j,k) = pu(i,j,k)*crx(i,j,k)
      enddo
      enddo

      do j=j1,j2
      do i=1,im-1
         dpi(i,j,k) = dpi(i,j,k) + xmass(i,j,k) - xmass(i+1,j,k)
      enddo
      enddo

      do j=j1,j2
         dpi(im,j,k) = dpi(im,j,k) + xmass(im,j,k) - xmass(1,j,k)
      enddo

      enddo ! vertical layer loop

!****6***0*********0*********0*********0*********0*********0**********72
! Compute Courant number at cell center (upwinding)
!****6***0*********0*********0*********0*********0*********0**********72

! E-W direction

      do k=1,nl

      do j=j1,j2
      do i=1,im-1
         if(crx(i,j,k)*crx(i+1,j,k) .gt. 0.) then
            if(crx(i,j,k) .gt. 0.) then
            u(i,j,k) = crx(i,j,k)
            else
            u(i,j,k) = crx(i+1,j,k)
            endif
         else
            u(i,j,k) = 0.
         endif
      enddo
      enddo

      i=im
      do j=j1,j2
      if(crx(i,j,k)*crx(1,j,k) .gt. 0.) then
         if(crx(i,j,k) .gt. 0.) then
         u(i,j,k) = crx(i,j,k)
         else
         u(i,j,k) = crx(1,j,k)
         endif
      else
         u(i,j,k) = 0.
      endif
      enddo

      enddo ! vertical layer loop

! N-S direction

      do k=1,nl

      do j=j1,j2
      do i=1,im
         if(cry(i,j,k)*cry(i,j+1,k) .gt. 0.) then
            if(cry(i,j,k) .gt. 0.) then
            v(i,j,k) = cry(i,j,k)/dlat(j)
            else
            v(i,j,k) = cry(i,j+1,k)/dlat(j)
            endif
         else
            v(i,j,k) = 0.
         endif
      enddo
      enddo

!++ to be checked ++
!     do 139 i=1,imh
!     v(i,     1,k) = 0.5*(cry(i,2,k)-cry(i+imh,2,k))
!     v(i+imh, 1,k) = -v(i,1,k)
!     v(i,    jm,k) = 0.5*(cry(i,jm,k)-cry(i+imh,jm1,k))
!139   v(i+imh,jm,k) = -v(i,jm,k)
!== to be checked ==

      enddo ! vertical layer loop

!****6***0*********0*********0*********0*********0*********0**********72
! Compute vertical mass flux (same dimensional unit as PS)
!****6***0*********0*********0*********0*********0*********0**********72

! compute total column mass CONVERGENCE.

!MIC$ do all autoscope shared(im,jm,DPI,PS1,PS2,W,DBK)
!MIC$* shared(DPI,PS1,PS2,W,DBK)
!MIC$* private(i,j,k,DG1)

      do j=1,jm

         dg1(1:im) = 0.
         do k=1,nl
            dg1(1:im) = dg1(1:im) + dpi(1:im,j,k)
         enddo

! Compute PS2 (PS at n+1) using the hydrostatic assumption.
! Changes (increases) to surface pressure = total column mass convergence

         ps2(1:im,j)  = ps1(1:im,j) + dg1(1:im)

! compute vertical mass flux from mass conservation principle:
! lower boundary of the first (uppermost) layer

         w(1:im,j,1) = dpi(1:im,j,1) - dbk(1)*dg1(1:im)

         do k=2,nl-1
         w(1:im,j,k) = w(1:im,j,k-1) + dpi(1:im,j,k) - dbk(k)*dg1(1:im)
         enddo

! Earth's surface

         w(1:im,j,nl) = 0.

      enddo

!MIC$ do all
!MIC$* shared(deform,NL,im,jm,delp,delp1,delp2,DPI,DAP,DBK,PS1,PS2)
!MIC$* private(i,j,k)

      do k=1,nl

         delp1(:,:,k) = dap(k) + dbk(k)*ps1(:,:)
         delp2(:,:,k) = dap(k) + dbk(k)*ps2(:,:)
         delp (:,:,k) = delp1(:,:,k) + dpi(:,:,k)

! Check deformation of the flow fields

        if(deform) then
          do j=1,jm
          do i=1,im
          if(delp(i,j,k) .le. 0.) then
             write(*,'(a)') '* FFSL transport'
             write(*,'(a,i3,a)') ' Vertical layer',k, &
                         'Noisy wind fields -> delp* is negative!'
             write(*,*) '* Smooth the wind fields or reduce time step'
             stop
          endif
          enddo
          enddo
        endif

      enddo !vertical layer loop

!****6***0*********0*********0*********0*********0*********0**********72
! Do transport one tracer at a time.
!****6***0*********0*********0*********0*********0*********0**********72

      DO 5000 ic=1,nc

!MIC$ do all autoscope
!MIC$* shared(q,DQ,delp1,U,V,j1,j2,JS,JN,im,jm,IML,IC,IORD,JORD)
!MIC$* shared(CRX,CRY,PU,xmass,ymass,fx,fy,acosp,rcap,qz)
!MIC$* private(i,j,k,jt,wk,DG2)

      do 2500 k=1,nl

! for the polar caps: replace the mixing ratio at the northest (southest)
! latitude with the zonal mean.

         sum1 = sum(q(:,1 ,k,ic))/im
         sum2 = sum(q(:,jm,k,ic))/im
         q(:,1 ,k,ic) = sum1
         q(:,jm,k,ic) = sum2

! Initialize DQ

         dq(:,:,k) = q(:,:,k,ic)*delp1(:,:,k)

! E-W advective cross term
      call xadv(im,jm,j1,j2,q(1,1,k,ic),u(1,1,k),js(k),jn(k),iml, &
                wk(1,1,1))

      wk(:,:,1) = q(:,:,k,ic) + 0.5*wk(:,:,1)

! N-S advective cross term
      do 66 j=j1,j2
      do 66 i=1,im
      jt = float(j) - v(i,j,k)
66    wk(i,j,2) = v(i,j,k) * (q(i,jt,k,ic) - q(i,jt+1,k,ic))

      do 77 j=j1,j2
      do 77 i=1,im
77    wk(i,j,2) = q(i,j,k,ic) + 0.5*wk(i,j,2)

!****6***0*********0*********0*********0*********0*********0**********72
! compute flux in  E-W direction
      call xtp(im,jm,iml,j1,j2,jn(k),js(k),pu(1,1,k),dq(1,1,k), &
               wk(1,1,2),crx(1,1,k),fx(1,1,k),xmass(1,1,k),iord)

! compute flux in  N-S direction
      call ytp(im,jm,j1,j2,acosp,dlat,rcap,dq(1,1,k),wk(1,1,1), &
               cry(1,1,k),dg2,ymass(1,1,k),wk(1,1,3),fy(1,1,k),jord)
!****6***0*********0*********0*********0*********0*********0**********72

      if(ZCROSS) then

! qz is the horizontal advection modified value for input to the
! vertical transport operator FZPPM
! Note: DQ contains only first order upwind contribution.

      do 88 j=1,JM
      do 88 i=1,IM
88    qz(i,j,k) = DQ(i,j,k) / delp(i,j,k)

      else

      do 99 j=1,JM
      do 99 i=1,IM
99    qz(i,j,k) = q(i,j,k,IC)

      endif

2500  continue     ! k-loop

!****6***0*********0*********0*********0*********0*********0**********72
! Compute fluxes in the vertical direction
      call FZPPM(qz,fz,IM,JM,NL,DQ,W,delp,KORD)
!****6***0*********0*********0*********0*********0*********0**********72

      if( MFCT ) then

      if (debug) write(*,*) '* mfct on'

! qlow is the low order "monotonic" solution

!MIC$ do all
!MIC$* shared(NL,im,jm,j1,jm1,qlow,DQ,delp2)
!MIC$* private(i,j,k)

      DO k=1,NL

      DO 560 j=1,JM
      DO 560 i=1,IM
560   qlow(i,j,k) = DQ(i,j,k) / delp2(i,j,k)

      if(j1.ne.2) then
      DO 561 i=1,IM
      qlow(i,   2,k) = qlow(i, 1,k)
      qlow(i,jm-1,k) = qlow(i,jm,k)
561   CONTINUE
      endif

      enddo

!****6***0*********0*********0*********0*********0*********0**********72
      call FCT3D(Q(1,1,1,IC),qlow,fx,fy,fz,IM,JM,NL,j1,j2,delp2, &
                 DPI,qz,wk,Qmax,Qmin,DG2,U,V,acosp,dlat,RCAP)
! Note: Q is destroyed!!!
!****6***0*********0*********0*********0*********0*********0**********72
      ENDIF

! Final update

!MIC$ do all autoscope
!MIC$* private(i,j,k,sum1,sum2)

      do 101 k=1,NL

      do 425 j=j1,j2
      do 425 i=1,IM
      DQ(i,j,k) = DQ(i,j,k) +  fx(i,j,k) - fx(i+1,j,k)                   &
                            + (fy(i,j,k) - fy(i,j+1,k))*acosp(j)/dlat(j) &
                            +  fz(i,j,k) - fz(i,j,k+1)
425   continue

! poles:
      sum1 = fy(IM,j1  ,k)
      sum2 = fy(IM,J2+1,k)

      do i=1,IM-1
      sum1 = sum1 + fy(i,j1  ,k)
      sum2 = sum2 + fy(i,J2+1,k)
      enddo

      DQ(1, 1,k) = DQ(1, 1,k) - sum1*RCAP + fz(1, 1,k) - fz(1, 1,k+1)
      DQ(1,JM,k) = DQ(1,JM,k) + sum2*RCAP + fz(1,JM,k) - fz(1,JM,k+1)

      do i=2,IM
      DQ(i, 1,k) = DQ(1, 1,k)
      DQ(i,JM,k) = DQ(1,JM,k)
      enddo

101   continue

!****6***0*********0*********0*********0*********0*********0**********72
      if(FILL) call qckxyz(DQ,DG2,IM,JM,NL,j1,j2,cosp,acosp,IC)
!****6***0*********0*********0*********0*********0*********0**********72

!MIC$ do all
!MIC$* shared(q,IC,NL,j1,im,jm,jm1,DQ,delp2)
!MIC$* private(i,j,k)

!******************************************************************
! finally, convert dq to q
!******************************************************************
! We have two options here:
! - to conserve the total mass of each tracer (cnst=2), use the
!   surface pressure predicted by the dynamical core
! - to preserve constant mixing ration, use the surface pressure
!   calculated in this subroutine (i.e., the updated ps2)

      select case (cnst)
      case(1)
        q(:,:,:,IC) = DQ(:,:,:) / delp2(:,:,:)
      case(2)
        q(:,:,:,IC) = DQ(:,:,:) / delp2dyn(:,:,:)
      end select

5000  continue !tracer loop

      if (debug) then
         write(*,*) '* Leaving routine tpcore'
         do ic=1,nc
            write(*,*) '* tracer',ic,'* max. & min. q =', &
                       maxval(q(:,:,:,ic)),minval(q(:,:,:,ic))
         end do
      endif

      RETURN
      END

!****6***0*********0*********0*********0*********0*********0**********72
      subroutine FCT3D(P,plow,fx,fy,fz,im,jm,km,j1,j2,delp,adx,ady, &
                       wk1,Qmax,Qmin,wkx,CRX,CRY,acosp,dlat,RCAP)
!****6***0*********0*********0*********0*********0*********0**********72

! MFCT Limiter
! plow: low order solution matrix
! P: current solution matrix
      integer im,jm,k,jm1,km,j1,j2,i,j,jt,it
      real ps1,ps2,pn1,pn2,ain,aou,bou,cin,cou,rcap,btop,bdon,bin
      real esl
      PARAMETER (esl = 1.E-30)
      REAL P(IM,JM,km),CRX(IM,JM,km),CRY(IM,JM,km),plow(IM,JM,km), &
           Qmax(IM,JM,km),Qmin(IM,JM,km),acosp(*),delp(im,jm,km), &
           adx(IM,JM,km),ady(IM,JM,km),fx(IM+1,JM,km), &
           fy(IM,JM,km),fz(im,jm,km+1),wk1(IM,JM,km), &
           wkx(im,jm),wkn(im,jm),dlat(jm)

      JM1 = JM-1

! Find local min/max of the low-order monotone solution
      call hilo3D(P,im,jm,km,j1,j2,adx,ady,Qmax,Qmin,wkx,wkn)
      call hilo3D(plow,im,jm,km,j1,j2,Qmax,Qmin,wk1,P,wkx,wkn)
! P is destroyed!

!     GOTO 123

!MIC$ do all autoscope
!MIC$* shared(im,j1,j2,km,CRX,CRY,adx,ady,Qmax,Qmin)
!MIC$* private(i,j,k,IT,JT,PS1,PS2,PN1,PN2)

      DO 1000 k=1,km
      do j=j1,j2
      DO i=1,IM

      IT = NINT( float(i) - CRX(i,j,k) )
! Wrap around in E-W
      if(IT .lt. 1) then
            IT = IM + IT
      elseif(IT .GT. IM) then
            IT = IT - IM
      endif

      JT = NINT( float(j) - CRY(i,j,k) )
      Qmax(i,j,k) = max(Qmax(i,j,k), adx(IT,JT,k))
      Qmin(i,j,k) = min(Qmin(i,j,k), ady(IT,JT,k))
      enddo
      enddo

! Poles:
      PS1 = max(Qmax(1, 1,k), adx(1, 1,k))
      PS2 = min(Qmin(1, 1,k), ady(1, 1,k))

      PN1 = max(Qmax(1,JM,k), adx(1,JM,k))
      PN2 = min(Qmin(1,JM,k), ady(1,JM,k))
      DO i=1,IM
      Qmax(i, 1,k) = PS1
      Qmin(i, 1,k) = PS2

      Qmax(i,JM,k) = PN1
      Qmin(i,JM,k) = PN2
      enddo
1000  continue

123   continue
! Flux Limiter

!MIC$ do all autoscope
!MIC$* shared(adx,ady,fx,fy,fz,plow,Qmax,Qmin,delp)
!MIC$* private(wkx,wkn)
!MIC$* private(i,j,k,ain,aou,bin,bou,cin,cou,btop,bdon)

      DO 2000 k=1,km

      DO j=j1,j2
      DO i=1,IM
         if(fx(i,j,k) .gt. 0.) then
         Ain = fx(i,j,k)
         Aou = 0.
         else
         Ain = 0.
         Aou = -fx(i,j,k)
         endif

         if(fx(i+1,j,k) .gt. 0.) then
         Aou = Aou + fx(i+1,j,k)
         else
         Ain = Ain - fx(i+1,j,k)
         endif

         if(fy(i,j,k) .gt. 0.) then
         Bin = fy(i,j,k)
         Bou = 0.
         else
         Bin = 0.
         Bou = -fy(i,j,k)
         endif

         if(fy(i,j+1,k) .gt. 0.) then
         Bou = Bou + fy(i,j+1,k)
         else
         Bin = Bin - fy(i,j+1,k)
         endif

         if(fz(i,j,k) .gt. 0.) then
         Cin = fz(i,j,k)
         Cou = 0.
         else
         Cin = 0.
         Cou = -fz(i,j,k)
         endif

         if(fz(i,j,k+1) .gt. 0.) then
         Cou = Cou + fz(i,j,k+1)
         else
         Cin = Cin - fz(i,j,k+1)
         endif

!HW 2009 -
!****6***0*********0*********0*********0*********0*********0**********72
         wkx(i,j) = Ain + Bin*acosp(j)/dlat(j) + Cin
         wkn(i,j) = Aou + Bou*acosp(j)/dlat(j) + Cou
!****6***0*********0*********0*********0*********0*********0**********72
!HW 2009 -
      enddo
      enddo

      DO j=j1,j2
      DO i=1,IM
      adx(i,j,k) = delp(i,j,k)*(Qmax(i,j,k)-plow(i,j,k))/(wkx(i,j)+esl)
      ady(i,j,k) = delp(i,j,k)*(plow(i,j,k)-Qmin(i,j,k))/(wkn(i,j)+esl)
      enddo
      enddo

! S Pole
      Ain = 0.
      Aou = 0.
      DO i=1,IM
      if(fy(i,j1,k).gt. 0.) then
           Aou = Aou + fy(i,j1,k)
      else
           Ain = Ain + fy(i,j1,k)
      endif
      enddo
      Ain = -Ain * RCAP
      Aou =  Aou * RCAP

! add vertical contribution...

      i=1
      j=1
      if(fz(i,j,k) .gt. 0.) then
      Cin = fz(i,j,k)
      Cou = 0.
      else
      Cin = 0.
      Cou = -fz(i,j,k)
      endif

      if(fz(i,j,k+1) .gt. 0.) then
      Cou = Cou + fz(i,j,k+1)
      else
      Cin = Cin - fz(i,j,k+1)
      endif

!****6***0*********0*********0*********0*********0*********0**********72
      btop = delp(1,1,k)*(Qmax(1,1,k)-plow(1,1,k))/(Ain+Cin+esl)
      bdon = delp(1,1,k)*(plow(1,1,k)-Qmin(1,1,k))/(Aou+Cou+esl)
!****6***0*********0*********0*********0*********0*********0**********72

      DO i=1,IM
      adx(i,j,k) = btop
      ady(i,j,k) = bdon
      enddo
! N Pole
      J=JM
      Ain = 0.
      Aou = 0.
      DO i=1,IM
      if(fy(i,j2+1,k).gt. 0.) then
           Ain = Ain + fy(i,j2+1,k)
      else
           Aou = Aou + fy(i,j2+1,k)
      endif
      enddo
      Ain =  Ain * RCAP
      Aou = -Aou * RCAP

! add vertical contribution...

      i=1
      if(fz(i,j,k) .gt. 0.) then
      Cin = fz(i,j,k)
      Cou = 0.
      else
      Cin = 0.
      Cou = -fz(i,j,k)
      endif

      if(fz(i,j,k+1) .gt. 0.) then
      Cou = Cou + fz(i,j,k+1)
      else
      Cin = Cin - fz(i,j,k+1)
      endif

!****6***0*********0*********0*********0*********0*********0**********72
      btop = delp(1,j,k)*(Qmax(1,j,k)-plow(1,j,k))/(Ain+Cin+esl)
      bdon = delp(1,j,k)*(plow(1,j,k)-Qmin(1,j,k))/(Aou+Cou+esl)
!****6***0*********0*********0*********0*********0*********0**********72

      DO i=1,IM
      adx(i,j,k) = btop
      ady(i,j,k) = bdon
      enddo

      if(j1 .ne. 2) then
      DO i=1,IM
! SP
      adx(i,2,k) = adx(i,1,k)
      ady(i,2,k) = ady(i,1,k)
! NP
      adx(i,JM1,k) = adx(i,JM,k)
      ady(i,JM1,k) = ady(i,JM,k)
      enddo
      endif
2000  continue

!MIC$ do all autoscope
!MIC$* shared(fz,adx,ady,im,jm,km)
!MIC$* private(i,j,k)

      DO 3000 k=1,km
      DO j=j1,j2
      do i=2,IM
      if(fx(i,j,k) .gt. 0.) then
      fx(i,j,k) = min(1.,ady(i-1,j,k),adx(i,j,k))*fx(i,j,k)
      else
      fx(i,j,k) = min(1.,adx(i-1,j,k),ady(i,j,k))*fx(i,j,k)
      endif
      enddo
      enddo

! For i=1
      DO j=j1,j2
      if(fx(1,j,k) .gt. 0.) then
      fx(1,j,k) = min(1.,ady(IM,j,k),adx(1,j,k))*fx(1,j,k)
      else
      fx(1,j,k) = min(1.,adx(IM,j,k),ady(1,j,k))*fx(1,j,k)
      endif
      fx(IM+1,j,k) = fx(1,j,k)
      enddo

      do j=j1,j2+1
      do i=1,IM
      if(fy(i,j,k) .gt. 0.) then
        fy(i,j,k) = min(1.,ady(i,j-1,k),adx(i,j,k))*fy(i,j,k)
      else
        fy(i,j,k) = min(1.,adx(i,j-1,k),ady(i,j,k))*fy(i,j,k)
      endif
      enddo
      enddo

      if(k .ne. 1) then
      do j=1,jm
      do i=1,im
      if(fz(i,j,k) .gt. 0.) then
        fz(i,j,k) = min(1.,ady(i,j,k-1),adx(i,j,k))*fz(i,j,k)
      else
        fz(i,j,k) = min(1.,adx(i,j,k-1),ady(i,j,k))*fz(i,j,k)
      endif
      enddo
      enddo
      endif

3000  continue

      return
      end

!****6***0*********0*********0*********0*********0*********0**********72
      subroutine ytp(IMR,JNP,j1,j2,acosp,dlat,RCAP,DQ,P,C,DC2 &
                    ,ymass,fy1,fy2,JORD)
!****6***0*********0*********0*********0*********0*********0**********72
!     do-loops re-written and "dlat" introduced by Hui Wan
!     for PlanetSimulator (2009-05).

      implicit none

      integer imr,jnp,j1,j2,jord
      REAL P(IMR,JNP),C(IMR,JNP),ymass(IMR,JNP),fy2(IMR,JNP), &
           DC2(IMR,JNP),DQ(IMR,JNP),acosp(JNP),dlat(JNP),rcap

! Work array
      REAL fy1(IMR,JNP)
      real sum1,sum2
      integer i,j,jt,jmr

      JMR = JNP - 1

      if(JORD.eq.1) then
!     upwind scheme

         do 1000 j=2,jnp
         do 1000 i=1,imr
            if ( c(i,j).le.0. ) then
               jt = j
            else
               jt = j-1
            endif
1000     fy1(i,j) = p(i,jt)

         do 1050 j=j1,j2+1
         do 1050 i=1,imr
1050     fy2(i,j) = 0.

      else

         call ymist(IMR,JNP,P,DC2)

         if(JORD.LE.0 .or. JORD.GE.3) then
!        PPM
!!           call fyppm(C,P,DC2,fy1,fy2,IMR,JNP,j1,j2,A6,AR,AL,JORD)
             write(*,*) 'PPM algorithm not yet adapted'
             stop

         else
!        van Leer
           do 1200 j=2,jnp
           do 1200 i=1,imr
              if ( c(i,j).le.0. ) then
                 jt = j
              else
                 jt = j-1
              endif
           fy1(i,j) = p(i,jt)
1200       fy2(i,j) = (sign(1.,C(i,j))-C(i,j)/dlat(jt))*DC2(i,jt)

         endif ! van Leer or PPM
      endif    ! upwind or higher order

      do 1300 j=2,jnp
      do 1300 i=1,imr
      fy1(i,j) = fy1(i,j)*ymass(i,j)
1300  fy2(i,j) = fy2(i,j)*ymass(i,j)

      DO 1400 j=j1,j2
      DO 1400 i=1,IMR
1400  DQ(i,j) = DQ(i,j) + (fy1(i,j) - fy1(i,j+1)) * acosp(j)/dlat(j)

! Poles
      sum1 = fy1(IMR,j1  )
      sum2 = fy1(IMR,J2+1)
      do i=1,IMR-1
      sum1 = sum1 + fy1(i,j1  )
      sum2 = sum2 + fy1(i,J2+1)
      enddo

      sum1 = DQ(1,  1) - sum1 * RCAP
      sum2 = DQ(1,JNP) + sum2 * RCAP
      do i=1,IMR
      DQ(i,  1) = sum1
      DQ(i,JNP) = sum2
      enddo

!     if(j1.ne.2) then
!     do i=1,IMR
!     DQ(i,  2) = sum1
!     DQ(i,JMR) = sum2
!     enddo
!     endif

      return
      end

!****6***0*********0*********0*********0*********0*********0**********72
      subroutine  ymist(IMR,JNP,P,DC)
!****6***0*********0*********0*********0*********0*********0**********72
!     re-written by Hui Wan for the PlanetSimulator (2009-05)

      implicit none

!     input and output

      integer imr,jnp
      real p(imr,jnp),dc(imr,jnp)

!     local variables

      integer imh,jmr,i,j
      real r24
      parameter ( r24 = 1./24. )
      real tmp,pmin,pmax
!
! 2nd order version for scalars
!
      imh = imr / 2
      jmr = jnp - 1

      do 10 j=2,jmr
      do 10 i=1,imr
      tmp = 0.25*(p(i,j+1) - p(i,j-1))
      pmax = max(p(i,j-1),p(i,j),p(i,j+1)) - p(i,j)
      pmin = p(i,j) - min(p(i,j-1),p(i,j),p(i,j+1))
      dc(i,j) = sign(min(abs(tmp),pmin,pmax),tmp)
10    continue
!
! Poles:
! Determine slopes in polar caps for scalars!

      do 20 i=1,IMH
! South
      tmp = 0.25*(p(i,2) - p(i+imh,2))
      Pmax = max(p(i,2),p(i,1), p(i+imh,2)) - p(i,1)
      Pmin = p(i,1) - min(p(i,2),p(i,1), p(i+imh,2))
      DC(i,1)=sign(min(abs(tmp),Pmax,Pmin),tmp)
! North.
      tmp = 0.25*(p(i+imh,JMR) - p(i,JMR))
      Pmax = max(p(i+imh,JMR),p(i,jnp), p(i,JMR)) - p(i,JNP)
      Pmin = p(i,JNP) - min(p(i+imh,JMR),p(i,jnp), p(i,JMR))
      DC(i,JNP) = sign(min(abs(tmp),Pmax,pmin),tmp)
20    continue

! Scalars:
      do 25 i=imh+1,IMR
      DC(i,  1) =  - DC(i-imh,  1)
      DC(i,JNP) =  - DC(i-imh,JNP)
25    continue

      return
      end

max_stars_repo_name: hadarohana/myCosmos

max_issues_repo_name: hadarohana/myCosmos

max_forks_repo_name: hadarohana/myCosmos


text:
module defns
  type point
     real :: x, y
  end type point
  type triangle
     type(point) :: a, b, c
  end type triangle
contains
  type(triangle) function tfunc(p)
    type(point) :: p
    tfunc%a = p
    tfunc%b = p
    tfunc%c = p
  end function tfunc
  real function rfunc(p)
    real :: p
    rfunc = p
  end function rfunc
end module defns

program geom
  use defns, pt=>point, tr=>triangle
  type point
     character(len=10) :: s
  end type point
  type triangle
     real :: x
  end type triangle
  type(pt) :: p1
  type(tr) :: t1
  type(triangle) :: bermuda
  type(point) :: sharp
  t1 = tfunc(p1)
  t1%a = p1
  sharp%s = "foo"
  bermuda%x = 1.23
  bermuda%x = rfunc(p1)  ! arg mismatch: actual pt/point, dummy real
  t1 = tfunc(sharp)      ! arg mismatch; actual local point, dummy module point
  t1%a = sharp           ! assign mismatch: local point to module point
end program geom

max_stars_repo_name: kurtsansom/ftnchek

max_issues_repo_name: kurtsansom/ftnchek

max_forks_repo_name: kurtsansom/ftnchek


text:
*----------------------------------------------------------------------*
      subroutine set_len_str2(len_str,nc,na,
     &                       len_str_array,nsym,mxidxms,
     &                       graph_c,idxms_c,gam_c,hpvx_c,
     &                       graph_a,idxms_a,gam_a,hpvx_a,
     &                       hpvxseq,resort)
*----------------------------------------------------------------------*
*     set string lengths
*     if (resort): resort to HPVX major sequence (default: C A major)
*----------------------------------------------------------------------*
      implicit none

      include 'opdim.h'

      logical, intent(in) ::
     &     resort
      integer, intent(in) ::
     &     nsym,mxidxms,
     &     len_str_array(nsym,mxidxms,*),
     &     nc, na,
     &     graph_c(nc),idxms_c(nc),gam_c(nc),hpvx_c(nc),
     &     graph_a(na),idxms_a(na),gam_a(na),hpvx_a(na),
     &     hpvxseq(ngastp)
      integer, intent(out) ::
     &     len_str(nc+na)

      integer ::
     &     idx_ca, hpvx, idx_hpvx, idx_c, idx_a, idx_high

      if (na+nc.eq.0) return
c dbg
c      print *,'gam_a: ',gam_a
c dbg

      if (resort) then

        idx_ca = 0
        idx_high = ngastp
        do idx_hpvx = 1, idx_high
          hpvx = hpvxseq(idx_hpvx)
          do idx_c = 1, nc
            if (hpvx_c(idx_c).eq.hpvx) then
              idx_ca = idx_ca + 1
              len_str(idx_ca) =
     &         len_str_array(gam_c(idx_c),idxms_c(idx_c),graph_c(idx_c))
            end if
          end do
          do idx_a = 1, na
            if (hpvx_a(idx_a).eq.hpvx) then
              idx_ca = idx_ca + 1
              len_str(idx_ca) =
     &         len_str_array(gam_a(idx_a),idxms_a(idx_a),graph_a(idx_a))
            end if
          end do
        end do

      else
        idx_ca = 0
        do idx_c = 1, nc
          idx_ca = idx_ca + 1
          len_str(idx_ca) =
     &         len_str_array(gam_c(idx_c),idxms_c(idx_c),graph_c(idx_c))
        end do
        do idx_a = 1, na
          idx_ca = idx_ca + 1
c dbg
c        print *,'ddd ',idx_ca,idx_a,nsym,mxidxms
c        print *,'ddd ',gam_a(idx_a),idxms_a(idx_a),graph_a(idx_a)
c        print *,'ddd ',
c     &   len_str_array(gam_a(idx_a),idxms_a(idx_a),graph_a(idx_a))
c dbg
          len_str(idx_ca) =
     &         len_str_array(gam_a(idx_a),idxms_a(idx_a),graph_a(idx_a))
        end do

      end if

      end

ext: f

max_stars_repo_name: ak-ustutt/GeCCo-public

max_issues_repo_name: ak-ustutt/GeCCo-public

max_forks_repo_name: ak-ustutt/GeCCo-public


text:
c
c     version 1.3  Aug 1992
c
c     Head wave and diffraction routines for TRAMP
c
c     ----------------------------------------------------------------
c
      subroutine hdwave(is,npt,ifam,nr,xshot,a1,nskip,itt,l,iszero,
     +    idot,iray,iamp,caust,spamp,nrskip,irayp,xsmax,isw,ampsmt,
     +    istep,anglew)
c
c     generate head waves at a spacing of hws along bottom of layer
c     and calculate corresponding amplitudes if iamp>0
c
      include 'tramp.par'
      real real,ls1,vave(ppray),amps(pnrayf),
     +     ls2,ls21,len,lenmin,lh,vpave(ppray),vsave(ppray),ab(ppray/2)
      complex q1,q2,z,zc,z1
      integer w,itt(1)
      include 'tramp.com'
c
      nhdw=1
      w=iwave
      on2=-omega/2.
      tb=0.
      nbnd=0
      iturn=0
      idh=id
      fidh=fid
      l=layer
      ib=iblk
      ilr=(id+3)/2
      ihamp=0
      irhbnd=ircbnd
      ichbnd=iccbnd
      nchbnd=nccbnd+1
      vpk=vp(npt,1)
      vsk=vs(npt,1)
      vk=vr(npt,1)
c
c     make preliminary head wave amplitude calculations
c
      if(iamp.gt.0) then
        if(idump.eq.1) write(14,5) ifam,nr
5       format(/'head wave:   ifam=',i3,'  ir=',i5//
     +'irt iw1 iw2  vp1     vs1     vp2     vs2    angle     zc',
     +5x,'phase')
        ihamp=1
        vp1=vp(1,2)
        vs1=vs(1,2)
        if(w.eq.1) then
          v1=vp1
        else
          v1=vs1
        end if
        ls1=sqrt(vk/v1)
        vpave(1)=(vp(1,2)+vp(2,1))/2.
        vsave(1)=(vs(1,2)+vs(2,1))/2.
        if(w.eq.1) then
          vave(1)=vpave(1)
        else
          vave(1)=vsave(1)
        end if
        if(npt.gt.2) then
          do 60 i=2,npt-1
             vpave(i)=(vp(i,2)+vp(i+1,1))/2.
             vsave(i)=(vs(i,2)+vs(i+1,1))/2.
             if(w.eq.1) then
               vave(i)=vpave(i)
             else
               vave(i)=vsave(i)
             end if
             ab(i)=atan(sh(i-1))
60        continue
          do 10 i=2,npt-1
             az1=ar(i,1)+ab(i)
             az2=ar(i,2)+ab(i)
             ls1=ls1*cos(az1)/cos(az2)
10        continue
        end if
c
        ls21=0.
        do 20 i=1,npt-1
           ls21=ls21+tr(i)*vave(i)
20      continue
c
        q1=cmplx(0.,1.)
        if(npt.gt.2) then
          do 30 i=2,npt-1
             az=abs(ar(i,1)+ab(i))
             if(az.gt.pi2) az=pi-az
             if(az.gt.1.5704) az=1.5704
             call zprtz(2,w,w,vpave(i-1),vpave(i),
     +                  vsave(i-1),vsave(i),az,zc)
             if(idump.eq.1) write(14,25) 2,w,w,vpave(i-1),
     +         vsave(i-1),vpave(i),vsave(i),az*pi18,cabs(zc),
     +         atan2(aimag(zc),real(zc))*pi18
25           format(i2,2i4,6f8.3,f8.2)
             q1=q1*cabs(zc)
30        continue
        end if
c
        call block(xr(npt),l+1,ib1)
        vpbnd=(vm(l+1,ib1,2)-vm(l+1,ib1,1))*(xr(npt)-xbnd(l+1,ib1,1))/
     +        (xbnd(l+1,ib1,2)-xbnd(l+1,ib1,1))+vm(l+1,ib1,1)
        vsbnd=vpbnd*vsvp(l+1,ib1)
        if(w.eq.1) then
          vbnd=vpbnd
        else
          vbnd=vsbnd
        end if
        if((vk/vbnd).ge.1.) then
          tc=pi2
        else
          tc=asin(vk/vbnd)
        end if
        call zprtz(2,w,w,vpk,vpbnd,vsk,vsbnd,tc,zc)
        if(idump.eq.1) write(14,25) 2,w,w,vpk,vsk,vpbnd,vsbnd,tc*pi18,
     +                 cabs(zc),atan2(aimag(zc),real(zc))*pi18
        zl=.5*cabs(zc)*tan(tc)/omega
        z1=zl*q1/ls1
        len=0.
        vhw=vbnd
        if(idump.eq.1) write(14,1)
1       format(' ')
      end if
c
      call ttime(isw,xshot,npt,nr,a1,ifam,itt,iszero,iamp,2,idump,
     +           omega)
      qh=qa
      if(iray.ne.0) call pltray(npt,nskip,idot,irayps,0,0.)
c
      t1=time
      l1=layer+1
      xr(1)=xr(npt)
      zr(1)=zr(npt)
      call block(xr(1),l1,ib1)
      vpbnd=(vm(l1,ib1,2)-vm(l1,ib1,1))*(xr(1)-xbnd(l1,ib1,1))/
     +      (xbnd(l1,ib1,2)-xbnd(l1,ib1,1))+vm(l1,ib1,1)
      vsbnd=vpbnd*vsvp(l1,ib1)
      if(w.eq.1) then
        vbnd=vpbnd
      else
        vbnd=vsbnd
      end if
      vr(1,1)=0.0
1020  vp(1,2)=(vm(layer,iblk,4)-vm(layer,iblk,3))*(xr(1)-
     +        xbnd(layer,iblk,1))/(xbnd(layer,iblk,2)-
     +        xbnd(layer,iblk,1))+vm(layer,iblk,3)
      if(vp(1,2).gt.0.) go to 1010
      layer=layer-1
      call block(xr(1),layer,iblk)
      go to 1020
1010  vs(1,2)=vp(1,2)*vsvp(layer,iblk)
      if(w.eq.1) then
        vr(1,2)=vp(1,2)
      else
        vr(1,2)=vs(1,2)
      end if
      qr(1)=q(layer,iblk,(3-iwave)/2)
      iflag2=-2
      abnd=atan(s(layer,iblk,2))
1000  if(ibsmth.eq.0) then
        alpha=atan(s(layer,iblk,2))
      else
        npl=ifix((xr(1)-xmin-.001)/xsinc)+1
        npr=npl+1
        xpl=xmin+float(npl-1)*xsinc
        alpha=(cosmth(layer+1,npr)-cosmth(layer+1,npl))*(xr(1)-xpl)
     +        /xsinc+cosmth(layer+1,npl)
      end if
      alpha2=atan(s(layer,iblk,2))
      time=t1+tb
      nbnd=0
      ircbnd=irhbnd
      iccbnd=ichbnd
      nccbnd=nchbnd
      iwave=w
      if(icbnd(iccbnd).eq.nccbnd) then
        iwave=-iwave
        iccbnd=iccbnd+1
        if(iwave.eq.1) then
          vr(1,2)=vp(1,2)
        else
          vr(1,2)=vs(1,2)
        end if
      end if
      if((vr(1,2)/vbnd).ge.1.) then
        a=pi2
      else
        a=pi-asin(vr(1,2)/vbnd)
      end if
      if(abs(abnd-alpha).le.angbnd) then
        idray(1)=l
        idray(2)=3
        ar(1,2)=fid*a-alpha
        npt=1
        call trace(npt,ifam,nr,iturn,iamp,ihamp,iflag)
c
        qa=qh
        call ttime(isw,xshot,npt,nr,a1,ifam,itt,iszero,iamp,
     +             iflag2,idump,omega)
        if((iray.eq.1.or.(iray.eq.2.and.vr(npt,2).gt.0.)).
     +    and.mod(nhdw-1,nrskip).eq.0)
     +    call pltray(npt,0,idot,irayps,istep,anglew)
c
c       complete head wave amplitude calculations
c
        if(iamp.gt.0.and.vr(npt,2).ne.0.) then
          if(len.lt.001) then
            xhmin=abs(xshot-xr(npt))
            lenmin=hedcut*xhmin
          end if
          if(len.lt.lenmin) then
            lh=lenmin
          else
            lh=len
          end if
c
          do 70 i=1,npt-1
             vpave(i)=(vp(i,2)+vp(i+1,1))/2.
             vsave(i)=(vs(i,2)+vs(i+1,1))/2.
             if(w.eq.1) then
               vave(i)=vpave(i)
             else
               vave(i)=vsave(i)
             end if
70        continue
          ls2=0.
          do 40 i=1,npt-1
             ls2=ls2+tr(i)*vave(i)
40        continue
          ls2=sqrt((lh*vhw+ls21)/v1)
c
          q2=cmplx(1.,0.)
          if(npt.gt.2) then
            do 50 i=2,npt-1
               az=abs(pi-fid1*(ar(i,1)+atan(sh(i-1))))
               if(az.gt.pi2) az=pi-az
               if(az.gt.1.5704) az=1.5704
               call zprtz(2,w,w,vpave(i-1),vpave(i),vsave(i-1),
     +                    vsave(i),az,zc)
               if(idump.eq.1) write(14,25) 2,w,w,vpave(i-1),
     +           vsave(i-1),vpave(i),vsave(i),az*pi18,cabs(zc),
     +           atan2(aimag(zc),real(zc))*pi18
               q2=q2*cabs(zc)
50          continue
          end if
c
          call zprtz(-2,w,w,vpbnd,vp(1,2),vsbnd,vs(1,2),1.570796,zc)
          if(idump.eq.1) write(14,25) 2,w,w,vpbnd,vsbnd,vp(1,2),
     +      vs(1,2),90.,cabs(zc),atan2(aimag(zc),real(zc))*pi18
          q2=q2*zc
c
          if(surcon.eq.1) then
            if(aamp(nbnd,1).gt.1.5704) aamp(nbnd,1)=1.5704
            call zprtz(3,w,icomp,vp(npt,1),.001,vs(npt,1),.001,
     +                 aamp(nbnd,1),zc)
            if(idump.eq.1) write(14,25) 3,w,w,vp(npt,1),vs(npt,1),
     +        0.,0.,aamp(nbnd,1)*pi18,cabs(zc),atan2(aimag(zc),
     +        real(zc))*pi18
            q2=q2*zc
          end if
c
          z=z1*(vhw/lh)**1.5*q2/(sqrt(vk)*ls2)
          if(surcon.eq.0) then
            if(icomp.eq.1) then
              if(w.eq.1) then
                z=z*abs(cos(ar(npt,1)))
              else
                z=z*abs(sin(ar(npt,1)))
              end if
            else
              if(w.eq.1) then
                z=z*abs(sin(ar(npt,1)))
              else
                z=z*abs(cos(ar(npt,1)))
              end if
            end if
          end if
          z=z*qa
          if(icbnd(1).eq.0) then
            z=z*spamp
          end if
          amp(ntt-1)=cabs(z)
          zimag=aimag(z)
          zreal=real(z)
          phase(ntt-1)=(caust+1)*pi2
          if(idump.eq.1) write(14,35) lh,vhw,qa
35        format(/'hdw bnd length=',f10.3,'  vhw=',f10.3,
     +            '  q fctr=',f10.4/)
        end if
      end if
      id=idh
      fid=fidh
      hwsl=hws
500   xrp=xr(1)+fid*hwsl*cos(alpha2)
      if(xrp.gt.xmin.and.xrp.lt.xmax) then
        vh1=vbnd
        la=l1
        if((fid*xrp).gt.(fid*xbnd(la,ib1,ilr))) then
          hl=abs(xbnd(la,ib1,ilr)-xr(1))/cos(alpha2)
          hwsl=hwsl-hl
          vpbnd=vm(la,ib1,ilr)
          vsbnd=vpbnd*vsvp(la,ib1)
          if(w.eq.1) then
            vbnd=vpbnd
          else
            vbnd=vsbnd
          end if
          vh2=vbnd
          if(abs(vh2-vh1).lt..001) then
            tb=tb+hl/vh1
          else
            tb=tb+hl*alog(vh2/vh1)/(vh2-vh1)
          end if
          if(iamp.gt.0) then
            vhw=(len*vhw+hl*(vh2+vh1)/2.)/(len+hl)
            len=len+hl
          end if
          xr(1)=xbnd(la,ib1,ilr)
          ib1=ib1+id
530       if(vm(la,ib1,1).eq.0.) then
            la=la+1
            if(la.gt.nlayer) go to 990
            call block(xr(1)+.001*fid,la,ib1)
            go to 530
          end if
          alpha=atan(s(la,ib1,1))
          vbnd=vm(la,ib1,-(id-3)/2)
          if(iq.eq.1) qh=qh*exp(on2*hl/
     +                   (((vh2+vh1)/2.)*q(la,ib1,(3-w)/2)))
          go to 500
        else
          xr(1)=xrp
          zr(1)=s(la,ib1,1)*xr(1)+b(la,ib1,1)
          if(xsmax.gt.0..and.abs(xshot-xr(1)).gt.xsmax) go to 990
          vr(1,1)=0.
          ls=l
520       call block(xr(1),ls,ib)
          if(vm(ls,ib,2).le.0.) then
            ls=ls-1
            if(ls.le.0) go to 990
            go to 520
          end if
510       vp(1,2)=(vm(ls,ib,4)-vm(ls,ib,3))*(xr(1)-xbnd(ls,ib,1))/
     +            (xbnd(ls,ib,2)-xbnd(ls,ib,1))+vm(ls,ib,3)
          vs(1,2)=vp(1,2)*vsvp(ls,ib)
          if(w.eq.1) then
            vr(1,2)=vp(1,2)
          else
            vr(1,2)=vs(1,2)
          end if
          la=l1
540       call block(xr(1),la,ib1)
          if(vm(la,ib1,1).eq.0.) then
            la=la+1
            if(la.gt.nlayer) go to 990
            go to 540
          end if
          vpbnd=(vm(la,ib1,2)-vm(la,ib1,1))*(xr(1)-xbnd(la,ib1,1))/
     +          (xbnd(la,ib1,2)-xbnd(la,ib1,1))+vm(la,ib1,1)
          vsbnd=vpbnd*vsvp(la,ib1)
          if(w.eq.1) then
            vbnd=vpbnd
          else
            vbnd=vsbnd
          end if
          if(vbnd.le..001) then
            write(11,45)
45          format('***  ray stopped - s-wave cannot propagate  ***')
            go to 990
          end if
          vh2=vbnd
          if(abs(vh2-vh1).lt..001) then
            tb=tb+hwsl/vh1
          else
            tb=tb+hwsl*alog(vh2/vh1)/(vh2-vh1)
          end if
          if(iamp.gt.0) then
            vhw=(len*vhw+hwsl*(vh2+vh1)/2.)/(len+hwsl)
            len=len+hwsl
          end if
          vratio=vr(1,2)/vbnd
          if(vratio.ge.1.) then
            a=pi2
          else
            a=pi-asin(vratio)
          end if
          if(iq.eq.1) qh=qh*exp(on2*hwsl/
     +                   (((vh2+vh1)/2.)*q(la,ib1,(3-w)/2)))
          layer=ls
          iblk=ib
          iflag2=-1
          nhdw=nhdw+1
          if(nhdw.gt.pnrayf) go to 990
          go to 1000
        end if
      end if
c
990   if(iamp.gt.0.and.itt(ifam).gt.0) then
        if(ampsmt.gt.0.) then
          nsmth=nint(ampsmt*(itt(ifam)-2))
          if(nsmth.gt.0) then
            do 9110 i=1,itt(ifam)
               if(amp(ntt-itt(ifam)-1+i).gt.0.) then
                 amps(i)=amp(ntt-itt(ifam)-1+i)
               else
                 amps(i)=0.
               end if
9110        continue
            do 9120 i=1,nsmth
               call smooth(amps,itt(ifam))
9120        continue
            do 9130 i=1,itt(ifam)
               amp(ntt-itt(ifam)-1+i)=amps(i)
9130        continue
          end if
        end if
c
        do 9140 i=1,itt(ifam)
           ipos=ntt-itt(ifam)-1+i
           write(13,15) isw,nr,range(ipos),amp(ipos),phase(ipos)*pi18,
     +                  rayid(ipos)
15         format(2i5,f10.3,e12.3,f10.2,f6.1)
9140    continue
      end if
c
      return
      end
c
c     ----------------------------------------------------------------
c
      subroutine difrct(is,npt,ifam,nr,xshot,a1,nskip,itt,iszero,
     +           idot,iray,nrskip,irayps,xsmax,isw,istep,anglew)
c
c     generate damped head waves at a spacing of hws along bottom of
c     layer to simulate diffracted arrivals associated with the top
c     of a low-velocity zone - amplitudes cannot be calculated
c
      include 'tramp.par'
      integer w,itt(1)
      include 'tramp.com'
c
      nhdw=1
      w=iwave
      difbnd=0
      tb=0.
      nbnd=0
      iturn=0
      idh=id
      fidh=fid
      l=layer
      ib=iblk
      ilr=(id+3)/2
      ihamp=0
      irhbnd=ircbnd
      ichbnd=iccbnd
c
      call ttime(isw,xshot,npt,nr,a1,ifam,itt,iszero,0,2,idump,
     +           omega)
      if(iray.ne.0) call pltray(npt,nskip,idot,irayps,0,0.)
c
      t1=time
      l1=layer+1
      xr(1)=xr(npt)
      zr(1)=zr(npt)
      call block(xr(1),l1,ib1)
      vpbnd=(vm(l1,ib1,2)-vm(l1,ib1,1))*(xr(1)-xbnd(l1,ib1,1))/
     +      (xbnd(l1,ib1,2)-xbnd(l1,ib1,1))+vm(l1,ib1,1)
      vsbnd=vpbnd*vsvp(l1,ib1)
      if(w.eq.1) then
        vbnd=vpbnd
      else
        vbnd=vsbnd
      end if
      abnd=atan(s(layer,iblk,2))
      a=pi-abs(ar(npt,1)+abnd)
      vr(1,1)=0.0
1020  vp(1,2)=(vm(layer,iblk,4)-vm(layer,iblk,3))*(xr(1)-
     +        xbnd(layer,iblk,1))/(xbnd(layer,iblk,2)-
     +        xbnd(layer,iblk,1))+vm(layer,iblk,3)
      if(vp(1,2).gt.0.) go to 1010
      layer=layer-1
      call block(xr(1),layer,iblk)
      go to 1020
1010  vs(1,2)=vp(1,2)*vsvp(layer,iblk)
      if(w.eq.1) then
        vr(1,2)=vp(1,2)
      else
        vr(1,2)=vs(1,2)
      end if
      iflag2=-2
1000  alpha=atan(s(layer,iblk,2))
      time=t1+tb
      nbnd=0
      ircbnd=irhbnd
      iccbnd=ichbnd
      if(abs(abnd-alpha).le.angbnd) then
        idray(1)=l
        idray(2)=3
        ar(1,2)=fid*a-alpha
        npt=1
        call trace(npt,ifam,nr,iturn,iamp,ihamp,iflag)
c
        call ttime(isw,xshot,npt,nr,a1,ifam,itt,iszero,0,iflag2,
     +             idump,omega)
        if((iray.eq.1.or.(iray.eq.2.and.vr(npt,2).gt.0.)).
     +    and.mod(nhdw-1,nrskip).eq.0)
     +    call pltray(npt,0,idot,irayps,istep,anglew)
c
        amp(ntt-1)=-1.
      end if
      id=idh
      fid=fidh
      hwsl=hws
500   xrp=xr(1)+fid*hwsl*cos(alpha2)
      if(xrp.gt.xmin.and.xrp.lt.xmax) then
        vh1=vbnd
        la=l1
        if((fid*xrp).gt.(fid*xbnd(la,ib1,ilr))) then
          hl=abs(xbnd(la,ib1,ilr)-xr(1))/cos(alpha2)
          hwsl=hwsl-hl
          vpbnd=vm(la,ib1,ilr)
          vsbnd=vpbnd*vsvp(la,ib1)
          if(w.eq.1) then
            vbnd=vpbnd
          else
            vbnd=vsbnd
          end if
          vh2=vbnd
          if(abs(vh2-vh1).lt..001) then
            tb=tb+hl/vh1
          else
            tb=tb+hl*alog(vh2/vh1)/(vh2-vh1)
          end if
          xr(1)=xbnd(la,ib1,ilr)
          ib1=ib1+id
530       if(vm(la,ib1,1).eq.0.) then
            la=la+1
            if(la.gt.nlayer) return
            call block(xr(1),la,ib1)
            go to 530
          end if
          alpha=atan(s(la,ib1,1))
          vbnd=vm(la,ib1,-(id-3)/2)
          go to 500
        else
          xr(1)=xrp
          zr(1)=s(la,ib1,1)*xr(1)+b(la,ib1,1)
          if(xsmax.gt.0..and.abs(xshot-xr(1)).gt.xsmax) return
          vr(1,1)=0.
          ls=l
520       call block(xr(1),ls,ib)
          if(vm(ls,ib,2).le.0.) then
            ls=ls-1
            if(ls.le.0) return
            go to 520
          end if
510       vp(1,2)=(vm(ls,ib,4)-vm(ls,ib,3))*(xr(1)-xbnd(ls,ib,1))/
     +            (xbnd(ls,ib,2)-xbnd(ls,ib,1))+vm(ls,ib,3)
          vs(1,2)=vp(1,2)*vsvp(ls,ib)
          if(w.eq.1) then
            vr(1,2)=vp(1,2)
          else
            vr(1,2)=vs(1,2)
          end if
          la=l1
540       call block(xr(1),la,ib1)
          if(vm(la,ib1,1).le.0.) then
            la=la+1
            if(la.gt.nlayer) return
            go to 540
          end if
          vpbnd=(vm(la,ib1,2)-vm(la,ib1,1))*(xr(1)-xbnd(la,ib1,1))/
     +          (xbnd(la,ib1,2)-xbnd(la,ib1,1))+vm(la,ib1,1)
          vsbnd=vpbnd*vsvp(la,ib1)
          if(w.eq.1) then
            vbnd=vpbnd
          else
            vbnd=vsbnd
          end if
          if(vbnd.le..001) then
            write(11,45)
45          format('***  ray stopped - s-wave cannot propagate  ***')
            return
          end if
          vh2=vbnd
          if(abs(vh2-vh1).lt..001) then
            tb=tb+hwsl/vh1
          else
            tb=tb+hwsl*alog(vh2/vh1)/(vh2-vh1)
          end if
          layer=ls
          iblk=ib
          iflag2=-1
          nhdw=nhdw+1
          if(nhdw.gt.pnrayf) return
          go to 1000
        end if
      end if
      return
      end


text:
!
! Copyright (c) 2014, NVIDIA CORPORATION.  All rights reserved.
!
! Licensed under the Apache License, Version 2.0 (the "License");
! you may not use this file except in compliance with the License.
! You may obtain a copy of the License at
!
!     http://www.apache.org/licenses/LICENSE-2.0
!
! Unless required by applicable law or agreed to in writing, software
! distributed under the License is distributed on an "AS IS" BASIS,
! WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
! See the License for the specific language governing permissions and
! limitations under the License.
!
! This tests line, file, date, and time directives,
! as well as the '#\n' (null directive)
!
#line 42 "universe.f90"
#
program p
    print *, __FILE__
    print *, __LINE__
    print *, __DATE__
    print *, __TIME__
    call check(.true., .true., 1)
end program

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
  CALL ConstructEmptyMatrix(matC, matA%actual_matrix_dimension, &
       & matA%process_grid, matA%is_complex)

  !$omp parallel
  !$omp do collapse(2)
  DO JJ = 1, matA%process_grid%number_of_blocks_columns
     DO II = 1, matA%process_grid%number_of_blocks_rows
        CALL PairwiseMultiplyMatrix(matA%LMAT(II,JJ), matB%LMAT(II,JJ), &
             & matC%LMAT(II,JJ))
     END DO
  END DO
  !$omp end do
  !$omp end parallel

max_stars_repo_name: Kokookster/NTPoly

max_issues_repo_name: omni-compiler/NTPoly

max_forks_repo_name: omni-compiler/NTPoly


text:
program testKernelLocalMem
!! Focal test program
!!

use Focal
use Focal_Test_Utils
use iso_fortran_env, only: sp=>real32, dp=>real64
implicit none

integer, parameter :: blockSize = 10

character(:), allocatable :: kernelSrc              ! Kernel source string
type(fclProgram) :: prog                            ! Focal program object
type(fclKernel) :: sumSqInt_k, sumSqFloat_k, sumSqDouble_k

real(sp), dimension(FCL_TEST_SIZE) :: hostReal32
real(dp), dimension(FCL_TEST_SIZE) :: hostReal64
integer, dimension(FCL_TEST_SIZE) :: hostInt32

type(fclDeviceFloat) :: deviceReal32, real32Sums
type(fclDeviceDouble) :: deviceReal64, real64Sums
type(fclDeviceInt32) :: deviceInt32, Int32Sums

integer :: i, nBlock

! --- Initialise ---
call fclTestInit()

! --- Initialise typed device buffers ---
call fclInitBuffer(deviceInt32,FCL_TEST_SIZE)
call fclInitBuffer(deviceReal32,FCL_TEST_SIZE)
call fclInitBuffer(deviceReal64,FCL_TEST_SIZE)

nBlock = FCL_TEST_SIZE/blockSize

call fclInitBuffer(real32Sums,nBLock)
call fclInitBuffer(real64Sums,nBLock)
call fclInitBuffer(int32Sums,nBLock)

! --- Initialise kernels ---
call fclGetKernelResource(kernelSrc)
prog = fclCompileProgram(kernelSrc)
sumSqInt_k = fclGetProgramKernel(prog,'sumSqInt32Test',[FCL_TEST_SIZE],[blockSize])
sumSqFloat_k = fclGetProgramKernel(prog,'sumSqFloatTest',[FCL_TEST_SIZE],[blockSize])
sumSqDouble_k = fclGetProgramKernel(prog,'sumSqDoubleTest',[FCL_TEST_SIZE],[blockSize])

! --- Initialise host arrays ---
hostInt32 = [(i,i=1,FCL_TEST_SIZE)]
hostReal32 = [(1.0*i,i=1,FCL_TEST_SIZE)]
hostReal64 = [(1.0d0*i,i=1,FCL_TEST_SIZE)]

! --- Transfer host to device ---
deviceInt32 = hostInt32
deviceReal32 = hostReal32
deviceReal64 = hostReal64

! --- Call kernels ---
call sumSqInt_k%launch(FCL_TEST_SIZE,blockSize,deviceInt32, &
                        fclLocalInt32(blockSize),int32Sums)
call sumSqFloat_k%launch(FCL_TEST_SIZE,blockSize,deviceReal32, &
                        fclLocalFloat(blockSize), real32Sums)
call sumSqDouble_k%launch(FCL_TEST_SIZE,blockSize,deviceReal64, &
                        fclLocalDouble(blockSize), real64Sums)

! --- Transfer device buffers to host ---
hostInt32(1:nBlock) = int32Sums
hostReal32(1:nBlock) = real32Sums
hostReal64(1:nBlock) = real64Sums

call fclWait()

! --- Check arrays ---
call fclTestAssert(sum(hostInt32(1:nBlock))==sum([(i*i,i=1,FCL_TEST_SIZE)]),'sumSqInt32Test')
call fclTestAssert(abs(sum(hostReal32(1:nBlock))-sum([(1.0*i*i,i=1,FCL_TEST_SIZE)])) < 1e-5,'sumSqReal32Test')
call fclTestAssert(abs(sum(hostReal64(1:nBlock))-sum([(1.0d0*i*i,i=1,FCL_TEST_SIZE)])) < 1e-10,'sumSqReal64Test')

call fclFreeBuffer(deviceInt32)
call fclFreeBuffer(deviceReal32)
call fclFreeBuffer(deviceReal64)

call fclTestFinish()

end program testKernelLocalMem
! -----------------------------------------------------------------------------

max_stars_repo_name: jborsen/test

max_issues_repo_name: jborsen/test

max_forks_repo_name: jborsen/test


text:
!**************************************************************************
!
! This code is copyright (c) 2009, Giovanni Bussi. When contacted by Miguel
! Caro via email (on 29/03/2021), regarding licensing, the author stated
! the following:
!
! "I don't have any specific requirement, it can be used at your will (even
! in a commercial code, and even if modified)."
!
! Therefore the code is redistributed according to these terms. The code
! has been slightly modified by Miguel Caro.
!
!**************************************************************************

module bussi

  contains

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! Stochastic velocity rescale, as described in
! Bussi, Donadio and Parrinello, J. Chem. Phys. 126, 014101 (2007)
!
! This subroutine implements Eq.(A7) and returns the new value for the kinetic energy,
! which can be used to rescale the velocities.
! The procedure can be applied to all atoms or to smaller groups.
! If it is applied to intersecting groups in sequence, the kinetic energy
! that is given as an input (kk) has to be up-to-date with respect to the previous rescalings.
!
! When applied to the entire system, and when performing standard molecular dynamics (fixed c.o.m. (center of mass))
! the degrees of freedom of the c.o.m. have to be discarded in the calculation of ndeg,
! and the c.o.m. momentum HAS TO BE SET TO ZERO.
! When applied to subgroups, one can chose to:
! (a) calculate the subgroup kinetic energy in the usual reference frame, and count the c.o.m. in ndeg
! (b) calculate the subgroup kinetic energy with respect to its c.o.m. motion, discard the c.o.m. in ndeg
!     and apply the rescale factor with respect to the subgroup c.o.m. velocity.
! They should be almost equivalent.
! If the subgroups are expected to move one respect to the other, the choice (b) should be better.
!
! If a null relaxation time is required (taut=0.0), the procedure reduces to an istantaneous
! randomization of the kinetic energy, as described in paragraph IIA.
!
! HOW TO CALCULATE THE EFFECTIVE-ENERGY DRIFT
! The effective-energy (htilde) drift can be used to check the integrator against discretization errors.
! The easiest recipe is:
! htilde = h + conint
! where h is the total energy (kinetic + potential)
! and conint is a quantity accumulated along the trajectory as minus the sum of all the increments of kinetic
! energy due to the thermostat.
!
function resamplekin(kk,T,ndeg,tau,dt)
  implicit none
  real*8               :: resamplekin
  real*8,  intent(in)  :: kk    ! present value of the kinetic energy of the atoms to be thermalized (in arbitrary units)
!  real*8,  intent(in)  :: sigma ! target average value of the kinetic energy (ndeg k_b T/2)  (in the same units as kk)
  integer, intent(in)  :: ndeg  ! number of degrees of freedom of the atoms to be thermalized
!
! This is the original definition. MAC introduced the time step explicitly, understanding that the thermostat
! is called at each individual step
!  real*8,  intent(in)  :: taut  ! relaxation time of the thermostat, in units of 'how often this routine is called'
!
! Use this instead
  real*8,  intent(in)  :: tau  ! relaxation time of the thermostat
  real*8,  intent(in)  :: dt    ! time step
  real*8,  intent(in)  :: T     ! temperature
  real*8 :: sigma, kB = 8.6173303d-5, taut
!
  real*8 :: factor,rr
!  real*8 :: gasdev

! MAC modified
  sigma = float(ndeg)*kB*T/2.0
  taut = tau/dt
!

  if(taut>0.1) then
    factor=exp(-1.0/taut)
  else
    factor=0.0
  end if
  rr = gasdev()
  resamplekin = kk + (1.0-factor)* (sigma*(sumnoises(ndeg-1)+rr**2)/ndeg-kk) &
               + 2.0*rr*sqrt(kk*sigma/ndeg*(1.0-factor)*factor)

end function

function sumnoises(nn)
  implicit none
  real*8 :: sumnoises
  integer, intent(in) :: nn
! returns the sum of n independent gaussian noises squared
! (i.e. equivalent to summing the square of the return values of nn calls to gasdev)
!  real*8 :: gamdev,gasdev
  if(nn==0) then
    sumnoises=0.0
  else if(nn==1) then
    sumnoises=gasdev()**2
  else if(modulo(nn,2)==0) then
    sumnoises=2.0*gamdev(nn/2)
  else
    sumnoises=2.0*gamdev((nn-1)/2) + gasdev()**2
  end if
end function

!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
! THE FOLLOWING ROUTINES ARE TRANSCRIBED FROM NUMERICAL RECIPES

function gamdev(ia)
! gamma-distributed random number, implemented as described in numerical recipes

implicit none
real*8 :: gamdev
integer, intent(in) :: ia
integer j
real*8 am,e,s,v1,v2,x,y
!real*8 :: ran1
!if(ia.lt.1)pause 'bad argument in gamdev'
if(ia.lt.1)then
  write(*,*) 'bad argument in gamdev'
  stop
end if
if(ia.lt.6)then
  x=1.
  do 11 j=1,ia
    x=x*ran1()
11  continue
  x=-log(x)
else
1 v1=2.*ran1()-1.
    v2=2.*ran1()-1.
  if(v1**2+v2**2.gt.1.)goto 1
    y=v2/v1
    am=ia-1
    s=sqrt(2.*am+1.)
    x=s*y+am
  if(x.le.0.)goto 1
    e=(1.+y**2)*exp(am*log(x/am)-s*y)
  if(ran1().gt.e)goto 1
endif
gamdev=x
end function

function gasdev()
! gaussian-distributed random number, implemented as described in numerical recipes

implicit none
real*8 :: gasdev
integer, save :: iset = 0
real*8, save :: gset
!real*8 :: ran1
real*8 fac,rsq,v1,v2
if(iset==0) then
1       v1=2.*ran1()-1.0d0
  v2=2.*ran1()-1.0d0
  rsq=v1**2+v2**2
  if(rsq.ge.1..or.rsq.eq.0.)goto 1
  fac=sqrt(-2.*log(rsq)/rsq)
  gset=v1*fac
  gasdev=v2*fac
  iset=1
else
  gasdev=gset
  iset=0
end if
end function

FUNCTION ran1()
! random number generator
INTEGER IA,IM,IQ,IR,NTAB,NDIV
REAL*8 ran1,AM,EPS,RNMX
PARAMETER (IA=16807,IM=2147483647,AM=1./IM,IQ=127773,IR=2836, &
  NTAB=32,NDIV=1+(IM-1)/NTAB,EPS=1.2e-7,RNMX=1.-EPS)
INTEGER j,k,iv(NTAB),iy
SAVE iv,iy
DATA iv /NTAB*0/, iy /0/
INTEGER, SAVE :: idum=0 !! ATTENTION: THE SEED IS HARDCODED
if (idum.le.0.or.iy.eq.0) then
  idum=max(-idum,1)
  do 11 j=NTAB+8,1,-1
    k=idum/IQ
    idum=IA*(idum-k*IQ)-IR*k
    if (idum.lt.0) idum=idum+IM
    if (j.le.NTAB) iv(j)=idum
11  continue
  iy=iv(1)
endif
k=idum/IQ
idum=IA*(idum-k*IQ)-IR*k
if (idum.lt.0) idum=idum+IM
j=1+iy/NDIV
iy=iv(j)
iv(j)=idum
ran1=min(AM*iy,RNMX)
return
      END FUNCTION

end module

max_stars_repo_name: richardjana/turbogap

max_issues_repo_name: richardjana/turbogap

max_forks_repo_name: richardjana/turbogap


text:
      SUBROUTINE DT_MODB(B,Nidx)

C***********************************************************************
C Sampling of impact parameter of collision.                           *
C    B          impact parameter    (output)                           *
C    NIDX       index of projectile/target material             (input)*
C Based on the original version by Shmakov et al.                      *
C This version dated 21.04.95 is revised by S. Roesler                 *
C Last change  5.5.2012 by S. Roesler.                                 *
C***********************************************************************

      IMPLICIT NONE
      DOUBLE PRECISION B , bb , DT_RNDM , ONE , ra , rb , TINY15 , TWO ,
     &                 x0 , x1 , x2 , y , y0 , y1 , y2 , ZERO
      INTEGER i0 , i1 , i2 , Nidx , ntarg
      SAVE

      INCLUDE 'inc/dtflka'

      PARAMETER (ZERO=0.0D0,TINY15=1.0D-15,ONE=1.0D0,TWO=2.0D0)

      LOGICAL left , lfirst

C central particle production, impact parameter biasing
      INCLUDE 'inc/dtimpa'
C emulsion treatment
      INCLUDE 'inc/dtcomp'
C Glauber formalism: parameters
      INCLUDE 'inc/dtglam'
C Glauber formalism: cross sections
      INCLUDE 'inc/dtglxs'

      DATA lfirst/.TRUE./

Cc                                                            mar'04 (26-02:10)
Cc -----------------------------------------------------------------------------
Cc  here we call into the pre-processed Glauber parametrization as derived for
Cc   dpmjet-II.5x - the called routine is a carbon-copy of this subroutine with
Cc  minor changes
Cc  also, by selecting KKMAT == -2 all the other routines in connection with
Cc  the Glauber formalism usually called are skipped
Cc
      IF ( Nidx.EQ.-2 ) THEN
         CALL FL_MODB(B,Nidx,y)
         RETURN
      END IF
Cc
Cc -----------------------------------------------------------------------------
Cc

      ntarg = ABS(Nidx)
      IF ( Nidx.LE.-1 ) THEN
         ra = RASh(1)
         rb = RBSh(ntarg)
      ELSE
         ra = RASh(ntarg)
         rb = RBSh(1)
      END IF

      IF ( ICEntr.NE.2 ) THEN
 50      y = DT_RNDM(bb)
         i0 = 1
         i2 = NSIteb
 100     i1 = (i0+i2)/2
         left = ((BSIte(0,1,ntarg,i0)-y)*(BSIte(0,1,ntarg,i1)-y))
     &          .LT.ZERO
         IF ( left ) THEN
            i2 = i1
         ELSE
            i0 = i1
         END IF
         IF ( i2-i0.LT.2 ) THEN
            i1 = i2 + 1
            IF ( i1.GT.NSIteb ) i1 = i0 - 1
         ELSE IF ( i2-i0.EQ.2 ) THEN
            i1 = i0 + 1
         ELSE
            GOTO 100
         END IF
         x0 = DBLE(i0-1)*BSTep(ntarg)
         x1 = DBLE(i1-1)*BSTep(ntarg)
         x2 = DBLE(i2-1)*BSTep(ntarg)
         y0 = BSIte(0,1,ntarg,i0)
         y1 = BSIte(0,1,ntarg,i1)
         y2 = BSIte(0,1,ntarg,i2)
C  80    CONTINUE
         B = x0*(y-y1)*(y-y2)/((y0-y1)*(y0-y2)+TINY15) + x1*(y-y0)
     &       *(y-y2)/((y1-y0)*(y1-y2)+TINY15) + x2*(y-y0)*(y-y1)
     &       /((y2-y0)*(y2-y1)+TINY15)
C*sr 5.4.98: shift B by half the bin width to be in agreement with BPROD
         B = B + 0.5D0*BSTep(ntarg)
         IF ( B.LT.ZERO ) B = x1
         IF ( B.GT.BMAx(ntarg) ) B = BMAx(ntarg)
         IF ( ICEntr.LT.0 ) THEN
            IF ( lfirst ) THEN
               lfirst = .FALSE.
               IF ( ICEntr.LE.-100 ) THEN
                  BIMin = 0.0D0
               ELSE
                  XSFrac = 0.0D0
               END IF
               CALL DT_GETBXS(XSFrac,BIMin,BIMax,ntarg)

               IF ( LPRi.GT.4 ) WRITE (LOUt,99010) RASh(1) , RBSh(ntarg)
     &              , BMAx(ntarg) , BIMin , BIMax , XSFrac*100.0D0 ,
     &              XSFrac*XSPro(1,1,ntarg)
99010          FORMAT (/,1X,'DT_MODB:      Biasing in impact parameter',
     &                 /,15X,'---------------------------'/,/,4X,
     &                 'average radii of proj / targ :',F10.3,' fm /',
     &                 F7.3,' fm',/,4X,'corresp. b_max (4*(r_p+r_t)) :',
     &                 F10.3,' fm',/,/,21X,'b_lo / b_hi :',F10.3,
     &                 ' fm /',F7.3,' fm',/,5X,'percentage of',
     &                 ' cross section :',F10.3,' %',/,5X,
     &                 'corresponding cross section :',F10.3,' mb',/)
            END IF
            IF ( ABS(BIMax-BIMin).LT.1.0D-3 ) THEN
               B = BIMin
            ELSE IF ( (B.LT.BIMin) .OR. (B.GT.BIMax) ) THEN
               GOTO 50
            END IF
         END IF
      ELSE IF ( ra.EQ.rb ) THEN
         bb = DT_RNDM(B)*(0.3D0*ra)**2
         B = SQRT(bb)
      ELSE IF ( ra.LT.rb ) THEN
         bb = DT_RNDM(B)*1.4D0*(rb-ra)**2
         B = SQRT(bb)
      ELSE IF ( ra.GT.rb ) THEN
         bb = DT_RNDM(B)*1.4D0*(ra-rb)**2
         B = SQRT(bb)
      END IF

      END SUBROUTINE

max_stars_repo_name: pzhristov/DPMJET

max_issues_repo_name: pzhristov/DPMJET

max_forks_repo_name: pzhristov/DPMJET


! written by jxzou at 20180331
! generate UHF natural orbtials (UNO) from UHF canonical orbitals

! updated by jxzou at 20180420: to support nalpha > nbeta
! updated by jxzou at 20180520: modify the intent of array alpha_coeff from (inout) to (in,copy)
! updated by jxzou at 20180825: fix the bug when only 1 pair
! updated by jxzou at 20191215: delete SVD generation of virtual/inactive MOs (PAO outside recommanded)
! updated by jxzou at 20200426: add NOON output
! updated by jxzou at 20210518: add an intent(in) parameter ON_thres

! This subroutine is designed to be imported as a module in Python.
!
! For GNU compiler, use
! --------------------------------------------
!  f2py -m uno -c pyuno.f90 --link-lapack_opt
! --------------------------------------------
!
! For Intel compiler, use
! -----------------------------------------------------------------------------------
!  f2py -m uno -c pyuno.f90 --link-lapack_opt --fcompiler=intelem --compiler=intelem
! -----------------------------------------------------------------------------------
!
!  to compile this file (a uno.so file will be generated). Then in Python
!  you can import the uno module.

subroutine uno(nbf, nif, nalpha, nbeta, alpha_coeff, beta_coeff, ao_ovlp, ON_thres, idx, noon, uno_coeff)
 implicit none
 integer, parameter :: iout = 6
 integer :: i, outid
 integer :: ndb, nact, nact0, nopen, nocc
 ! ndb: the number of doubly occupied MOs
 ! nact: the number of active occupied orbitals
 ! nact = nact0 + nopen
 ! nocc = nalpha + nact0

 integer :: nbf, nif, nalpha, nbeta
!f2py intent(in) :: nbf, nif, nalpha, nbeta
 ! nbf: the number of basis functions
 ! nif: the number of independent functions, i.e., the number of MOs
 ! nalpha: the number of alpha electrons
 ! nbeta: the number of beta electrons

 ! array idx is in Fortran convention, i.e., begins from 1, not 0
 integer :: idx(3) ! 1: start of pair, 2: start of virtual, 3: nopen
!f2py intent(out) :: idx
 integer, allocatable :: idx1(:), idx2(:)

 real(kind=8) :: maxv, abs_mean
 real(kind=8), parameter :: ON_criteria = 0.99999d0
 real(kind=8) :: ON_thres ! occupation number threshold
!f2py intent(in) :: ON_thres

 real(kind=8) :: ao_ovlp(nbf,nbf)
!f2py intent(in) :: ao_ovlp
!f2py depend(nbf) :: ao_ovlp

 real(kind=8) :: alpha_coeff(nbf,nif), beta_coeff(nbf,nif)
!f2py intent(in,copy) :: alpha_coeff
!f2py intent(in) :: beta_coeff
!f2py depend(nbf,nif) :: alpha_coeff, beta_coeff

 real(kind=8) :: noon(nif), uno_coeff(nbf,nif)
!f2py intent(out) :: noon, uno_coeff
!f2py depend(nif) :: noon
!f2py depend(nbf,nif) :: uno_coeff
 real(kind=8), allocatable :: mo_basis_ovlp(:,:), alpha_occ(:,:), beta_occ(:,:)
 real(kind=8), allocatable :: sv_occ0(:), sv_occ(:)
 character(len=7), parameter :: outname = 'uno.out'
 character(len=62), parameter :: on_warn1 = 'Warning in subroutine uno: ON_thres deviates from ON_criteria.'
 character(len=35), parameter :: on_warn2 = 'You better know what you are doing.'

 noon = 0d0
 nopen = nalpha - nbeta
 if(nopen < 0) then
  write(iout,'(A)') 'ERROR in subroutine uno: nalpha < nbeta.'
  write(iout,'(2(A,I0))') 'nalpha=', nalpha, ', nbeta=', nbeta
  stop
 end if

 open(newunit=outid,file=outname,status='replace')
 write(outid,'(A4,I5)') 'nbf=', nbf
 write(outid,'(A4,I5)') 'nif=', nif
 write(outid,'(A,F11.7)') 'ON_criteria=', ON_criteria
 write(outid,'(A,F11.7)') 'ON_thres=', ON_thres
 if(DABS(ON_thres - ON_criteria) > 1d-5) then
  write(outid,'(/,A)') on_warn1
  write(outid,'(A)') on_warn2
  write(iout,'(/,A)') on_warn1
  write(iout,'(A)') on_warn2
 end if

 ! check the orthonormality of initial Alpha and Beta MO, respectively
 allocate(mo_basis_ovlp(nif,nif))
 call get_mo_basis_ovlp(nif, nif, nbf, alpha_coeff, alpha_coeff, ao_ovlp, mo_basis_ovlp)
 call check_unity(nif, mo_basis_ovlp, maxv, abs_mean)
 write(outid,'(/,A)') 'The orthonormality of initial Alpha MO:'
 write(outid,'(A,F16.10)') 'maxv=', maxv
 write(outid,'(A,F16.10)') 'abs_mean=', abs_mean

 call get_mo_basis_ovlp(nif, nif, nbf, beta_coeff, beta_coeff, ao_ovlp, mo_basis_ovlp)
 call check_unity(nif, mo_basis_ovlp, maxv, abs_mean)
 deallocate(mo_basis_ovlp)
 write(outid,'(/,A)') 'The orthonormality of initial Beta MO:'
 write(outid,'(A,F16.10)') 'maxv=', maxv
 write(outid,'(A,F16.10)') 'abs_mean=', abs_mean
 ! check orthonormality done

 if(nbeta == 0) then ! no beta electrons, return
  forall(i = 1:nopen) noon(i) = 1d0
  uno_coeff = alpha_coeff
  idx = [1, nopen+1, nopen]
  write(outid,'(/,A6,I5)') 'ndb  =', 0
  write(outid,'(A6,I5)')   'nact =', nalpha
  write(outid,'(A6,I5)')   'nact0=', 0
  write(outid,'(A6,3I5)')  'idx  =', idx
  close(outid)
  return
 end if
 ! allocate arrays for alpha and beta occupied orbitals, prepare for SVD
 allocate(alpha_occ(nbf,nalpha), source=alpha_coeff(:,1:nalpha))
 allocate(beta_occ(nbf,nbeta), source=beta_coeff(:,1:nbeta))

 ! calculate the overlap between alpha and beta occupied spatial orbitals
 allocate(mo_basis_ovlp(nalpha,nbeta))
 call get_mo_basis_ovlp(nalpha, nbeta, nbf, alpha_occ, beta_occ, ao_ovlp, mo_basis_ovlp)
 ! calculate done

 ! do SVD on the alpha_beta_ovlp of occupied spatial orbitals
 allocate(sv_occ(nalpha))
 call svd_and_rotate(nalpha, nbeta, nbf, alpha_occ, beta_occ, mo_basis_ovlp, sv_occ, .false.)
 deallocate(mo_basis_ovlp)
 write(iout,'(/,A)') 'Singular values from SVD of Alpha/Beta MOs:'
 write(iout,'(5(1X,ES15.8))') (sv_occ(i), i=1,nalpha)
 ! SVD done in occ space

 allocate(sv_occ0(nalpha), source=sv_occ)
 nact = COUNT(sv_occ < ON_criteria)
 ndb = nalpha - nact
 nact0 = nact - nopen
 nocc = nalpha + nact0
 idx = [ndb+1, nocc+1, nopen]
 write(outid,'(/,A6,I5)') 'ndb  =', ndb
 write(outid,'(A6,I5)')   'nact =', nact
 write(outid,'(A6,I5)')   'nact0=', nact0
 write(outid,'(A6,3I5)')  'idx  =', idx

 ! generate NOON (Natural Orbital Occupation Number)
 forall(i = 1:nalpha) noon(i) = 1d0 + sv_occ(i)
 forall(i = 1:nact0) noon(nalpha+i) = 2d0 - noon(nbeta-i+1)
 deallocate(sv_occ)

 ! copy the doubly occupied MOs
 alpha_coeff(:,1:ndb) = alpha_occ(:,1:ndb)

 ! copy the singly occupied MO
 if(nopen > 0) alpha_coeff(:,nbeta+1:nalpha) = alpha_occ(:,nbeta+1:nalpha)

 ! transform the corresponding orbitals to UNOs in occ space
 allocate(idx1(nact0), idx2(nact0))
 forall(i = 1:nact0)
  idx1(i) = ndb + i
  idx2(i) = nocc + 1 - i
 end forall
 forall(i = 1:nact0)
  alpha_coeff(:,idx1(i)) = (alpha_occ(:,idx1(i)) + beta_occ(:,idx1(i)))/DSQRT(2d0*noon(idx1(i)))
  alpha_coeff(:,idx2(i)) = (alpha_occ(:,idx1(i)) - beta_occ(:,idx1(i)))/DSQRT(2d0*noon(idx2(i)))
 end forall
 deallocate(idx1, idx2, alpha_occ, beta_occ)
 ! done transform in occ space

 ! Set virtual MOs to be zero. They may be calculated in following PAO constructions
 if(nocc < nif) alpha_coeff(:,nocc+1:nif) = 0d0

 ! check the orthonormality of final Alpha MO
 allocate(mo_basis_ovlp(nocc,nocc))
 call get_mo_basis_ovlp(nocc, nocc, nbf, alpha_coeff(:,1:nocc), alpha_coeff(:,1:nocc), ao_ovlp, mo_basis_ovlp)
 call check_unity(nocc, mo_basis_ovlp, maxv, abs_mean)
 deallocate(mo_basis_ovlp)
 write(outid,'(/,A)') 'The orthonormality of final Alpha MO:'
 write(outid,'(A,F16.10)') 'maxv=', maxv
 write(outid,'(A,F16.10)') 'abs_mean=', abs_mean

 ! copy the UNO coefficients
 uno_coeff = alpha_coeff

 ! now let's update ndb, nact, nact0, idx according to input ON_thres
 nact = COUNT(sv_occ0 < ON_thres)
 ndb = nalpha - nact
 nact0 = nact - nopen
 nocc = nalpha + nact0
 idx = [ndb+1, nocc+1, nopen]
 write(outid,'(/,A6,I5)') 'ndb  =', ndb
 write(outid,'(A6,I5)')   'nact =', nact
 write(outid,'(A6,I5)')   'nact0=', nact0
 write(outid,'(A6,3I5)')  'idx  =', idx
 close(outid)
 deallocate(sv_occ0)
 return
end subroutine uno

subroutine get_mo_basis_ovlp(na, nb, nbf, c_alpha, c_beta, ao_ovlp, mo_basis_ovlp)
 implicit none
 integer, intent(in) :: na, nb, nbf
 real(kind=8), intent(in) :: c_alpha(nbf,na), c_beta(nbf,nb)
 real(kind=8), intent(in) :: ao_ovlp(nbf,nbf)
 real(kind=8), intent(out) :: mo_basis_ovlp(na,nb)
 real(kind=8) :: s_c_beta(nbf, nb)

 s_c_beta = 0d0
 mo_basis_ovlp = 0d0

 ! ?symm: Computes a matrix-matrix product where one input matrix is symmetric
 ! ?gemm: Computes a matrix-matrix product with general matrices
 ! Syntax FORTRAN 77:
 ! call dsymm(side, uplo, m, n, alpha, a, lda, b, ldb, beta, c, ldc)
 ! call dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
 call dsymm('L', 'U', nbf, nb, 1d0, ao_ovlp, nbf, c_beta, nbf, 0d0, s_c_beta, nbf)
 call dgemm('T', 'N', na, nb, nbf, 1d0, c_alpha, nbf, s_c_beta, nbf, 0d0, mo_basis_ovlp, na)
 return
end subroutine get_mo_basis_ovlp

! perform SVD on a given overlap matrix
subroutine svd_on_ovlp(m, n, a, u, vt, s)
 implicit none
 integer :: i, lwork
 integer, intent(in) :: m,n
 integer, parameter :: iout = 6
 real(kind=8), intent(in) :: a(m,n)
 real(kind=8), intent(out) :: u(m,m), vt(n,n), s(m)
 real(kind=8), allocatable :: work(:)
 real(kind=8) :: a_copy(m,n)

 a_copy = a
 u = 0d0; vt = 0d0; s = 0d0

 ! ?gesvd: Computes the singular value decomposition of a general rectangular matrix
 ! Syntax FORTRAN 77:
 ! call dgesvd(jobu, jobvt, m, n, a, lda, s, u, ldu, vt, ldvt, work, lwork, info)
 lwork = -1
 allocate(work(1))
 work = 0
 call dgesvd('A', 'A', m, n, a_copy, m, s, u, m, vt, n, work, lwork, i)
 lwork = CEILING(work(1))
 deallocate(work)
 allocate(work(lwork))
 call dgesvd('A', 'A', m, n, a_copy, m, s, u, m, vt, n, work, lwork, i)

 deallocate(work)
 if(i /= 0) then
  write(iout,'(A)') 'ERROR: info /= 0 in subroutine svd_on_ovlp. Please check why.'
  write(iout,'(A5,I0)') 'info=',i
  stop
 end if
 return
end subroutine svd_on_ovlp

subroutine svd_and_rotate(nalpha,nbeta,nbf,c_alpha,c_beta,alpha_beta_ovlp,sv,reverse)
 implicit none
 integer :: i
 integer, intent(in) :: nalpha, nbeta, nbf
 real(kind=8), intent(inout) :: c_alpha(nbf,nalpha), c_beta(nbf,nbeta)
 real(kind=8), intent(in) :: alpha_beta_ovlp(nalpha,nbeta)
 real(kind=8), intent(out) :: sv(nalpha)
 real(kind=8) :: u(nalpha,nalpha), vt(nbeta,nbeta)
 real(kind=8), allocatable :: sv2(:), u2(:,:), vt2(:,:)
 logical, intent(in) :: reverse

 sv = 0d0
 u = 0d0
 vt = 0d0

 ! perform SVD on alpha_beta_ovlp of occupied or virtual orbitals
 call svd_on_ovlp(nalpha, nbeta, alpha_beta_ovlp, u, vt, sv)

 ! rotate occupied or virtual orbitals
 if(reverse) then
  allocate(sv2(nalpha))
  forall(i = 1:nalpha) sv2(i) = sv(nalpha-i+1)
  sv = sv2
  deallocate(sv2)
  allocate(u2(nalpha,nalpha), vt2(nbeta,nbeta))
  forall(i = 1:nalpha) u2(:,i) = u(:,nalpha-i+1)
  u = u2
  forall(i = 1:nbeta) vt2(i,:) = vt(nbeta-i+1,:)
  vt = vt2
  deallocate(u2, vt2)
 end if

 ! ?gemm: Computes a matrix-matrix product with general matrices
 ! Syntax FORTRAN 77:
 ! call dgemm(transa, transb, m, n, k, alpha, a, lda, b, ldb, beta, c, ldc)
 allocate(u2(nbf,nalpha), source=0d0)
 allocate(vt2(nbf,nbeta), source=0d0)
 call dgemm('N', 'N', nbf, nalpha, nalpha, 1d0, c_alpha, nbf, u, nalpha, 0d0, u2, nbf)
 call dgemm('N', 'T', nbf, nbeta, nbeta, 1d0, c_beta, nbf, vt, nbeta, 0d0, vt2, nbf)

 c_alpha = u2
 c_beta = vt2
 deallocate(u2, vt2)
 return
end subroutine svd_and_rotate

subroutine check_unity(n, a, maxv, abs_mean)
 implicit none
 integer :: i, j
 integer, intent(in) :: n
 real(kind=8), intent(in) :: a(n,n)
 real(kind=8), intent(out) :: maxv, abs_mean
 real(kind=8) :: b(n,n)

 b = a
 forall(i = 1:n) b(i,i) = b(i,i) - 1d0
 forall(i=1:n, j=1:n) b(i,j) = DABS(b(i,j))

 maxv = MAXVAL(b)
 abs_mean = SUM(b)/DBLE(n*n)
 return
end subroutine check_unity

! calculated the SVD singular values of overlap of two sets of MOs
! Note: the input coeff1 and coeff2 must have the same dimension (nbf,nif)
subroutine svd_of_two_mo(nbf, nif, coeff1, coeff2, S)
 implicit none
 integer :: i
 integer :: nbf, nif
!f2py intent(in) :: nbf, nif
 integer, parameter :: iout = 6

 real(kind=8), intent(in) :: coeff1(nbf,nif), coeff2(nbf,nif), S(nbf,nbf)
!f2py intent(in) :: coeff1, coeff2, S
!f2py depend(nbf,nif) :: coeff1, coeff2
!f2py depend(nbf) :: S

 real(kind=8), allocatable :: mo_ovlp(:,:), u(:,:), vt(:,:), ev(:)

 allocate(mo_ovlp(nif,nif), source=0d0)
 call get_mo_basis_ovlp(nif, nif, nbf, coeff1, coeff2, S, mo_ovlp)

 allocate(u(nif,nif), vt(nif,nif), ev(nif))
 call svd_on_ovlp(nif, nif, mo_ovlp, u, vt, ev)
 deallocate(mo_ovlp, u, vt)

 write(iout,'(/,A)') 'SVD analysis of two sets of MOs:'
 write(iout,'(A,ES15.8)') 'The smallest singular value:', MINVAL(ev)
 i = COUNT(ev < 1d-1)
 write(iout,'(A,I0)') 'Number of singular values< 0.1: ', i
 i = COUNT(ev < 1d-2)
 write(iout,'(A,I0)') 'Number of singular values<0.01: ', i

 write(iout,'(A)') 'All singular values:'
 write(iout,'(5(1X,ES15.8))') (ev(i),i=1,nif)
 deallocate(ev)
 return
end subroutine svd_of_two_mo


text:
	SUBROUTINE realft_sp(data,isign,zdata)
	USE nrtype; USE nrutil, ONLY : assert,assert_eq,zroots_unity
	USE nr, ONLY : four1
	IMPLICIT NONE
	REAL(SP), DIMENSION(:), INTENT(INOUT) :: data
	INTEGER(I4B), INTENT(IN) :: isign
	COMPLEX(SPC), DIMENSION(:), OPTIONAL, TARGET :: zdata
	INTEGER(I4B) :: n,ndum,nh,nq
	COMPLEX(SPC), DIMENSION(size(data)/4) :: w
	COMPLEX(SPC), DIMENSION(size(data)/4-1) :: h1,h2
	COMPLEX(SPC), DIMENSION(:), POINTER :: cdata
	COMPLEX(SPC) :: z
	REAL(SP) :: c1=0.5_sp,c2
	n=size(data)
	call assert(iand(n,n-1)==0, 'n must be a power of 2 in realft_sp')
	nh=n/2
	nq=n/4
	if (present(zdata)) then
		ndum=assert_eq(n/2,size(zdata),'realft_sp')
		cdata=>zdata
		if (isign == 1) cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
	else
		allocate(cdata(n/2))
		cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
	end if
	if (isign == 1) then
		c2=-0.5_sp
		call four1(cdata,+1)
	else
		c2=0.5_sp
	end if
	w=zroots_unity(sign(n,isign),n/4)
	w=cmplx(-aimag(w),real(w),kind=spc)
	h1=c1*(cdata(2:nq)+conjg(cdata(nh:nq+2:-1)))
	h2=c2*(cdata(2:nq)-conjg(cdata(nh:nq+2:-1)))
	cdata(2:nq)=h1+w(2:nq)*h2
	cdata(nh:nq+2:-1)=conjg(h1-w(2:nq)*h2)
	z=cdata(1)
	if (isign == 1) then
		cdata(1)=cmplx(real(z)+aimag(z),real(z)-aimag(z),kind=spc)
	else
		cdata(1)=cmplx(c1*(real(z)+aimag(z)),c1*(real(z)-aimag(z)),kind=spc)
		call four1(cdata,-1)
	end if
	if (present(zdata)) then
		if (isign /= 1) then
			data(1:n-1:2)=real(cdata)
			data(2:n:2)=aimag(cdata)
		end if
	else
		data(1:n-1:2)=real(cdata)
		data(2:n:2)=aimag(cdata)
		deallocate(cdata)
	end if
	END SUBROUTINE realft_sp

	SUBROUTINE realft_dp(data,isign,zdata)
	USE nrtype; USE nrutil, ONLY : assert,assert_eq,zroots_unity
	USE nr, ONLY : four1
	IMPLICIT NONE
	REAL(DP), DIMENSION(:), INTENT(INOUT) :: data
	INTEGER(I4B), INTENT(IN) :: isign
	COMPLEX(DPC), DIMENSION(:), OPTIONAL, TARGET :: zdata
	INTEGER(I4B) :: n,ndum,nh,nq
	COMPLEX(DPC), DIMENSION(size(data)/4) :: w
	COMPLEX(DPC), DIMENSION(size(data)/4-1) :: h1,h2
	COMPLEX(DPC), DIMENSION(:), POINTER :: cdata
	COMPLEX(DPC) :: z
	REAL(DP) :: c1=0.5_dp,c2
	n=size(data)
	call assert(iand(n,n-1)==0, 'n must be a power of 2 in realft_dp')
	nh=n/2
	nq=n/4
	if (present(zdata)) then
		ndum=assert_eq(n/2,size(zdata),'realft_dp')
		cdata=>zdata
		if (isign == 1) cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
	else
		allocate(cdata(n/2))
		cdata=cmplx(data(1:n-1:2),data(2:n:2),kind=spc)
	end if
	if (isign == 1) then
		c2=-0.5_dp
		call four1(cdata,+1)
	else
		c2=0.5_dp
	end if
	w=zroots_unity(sign(n,isign),n/4)
	w=cmplx(-aimag(w),real(w),kind=dpc)
	h1=c1*(cdata(2:nq)+conjg(cdata(nh:nq+2:-1)))
	h2=c2*(cdata(2:nq)-conjg(cdata(nh:nq+2:-1)))
	cdata(2:nq)=h1+w(2:nq)*h2
	cdata(nh:nq+2:-1)=conjg(h1-w(2:nq)*h2)
	z=cdata(1)
	if (isign == 1) then
		cdata(1)=cmplx(real(z)+aimag(z),real(z)-aimag(z),kind=dpc)
	else
		cdata(1)=cmplx(c1*(real(z)+aimag(z)),c1*(real(z)-aimag(z)),kind=dpc)
		call four1(cdata,-1)
	end if
	if (present(zdata)) then
		if (isign /= 1) then
			data(1:n-1:2)=real(cdata)
			data(2:n:2)=aimag(cdata)
		end if
	else
		data(1:n-1:2)=real(cdata)
		data(2:n:2)=aimag(cdata)
		deallocate(cdata)
	end if
	END SUBROUTINE realft_dp

max_stars_repo_name: AntonBiryukovUofC/RayTracerFortran

max_issues_repo_name: AntonBiryukovUofC/RayTracerFortran

max_forks_repo_name: AntonBiryukovUofC/RayTracerFortran


text:
      INTEGER FUNCTION IDAMAX( N, X, INCX )

      INTEGER               INCX, N
      DOUBLE PRECISION      X(*)

      INTEGER               TEMP

      CALL IDAMAX_SUB( N, X, INCX, TEMP )
      IDAMAX = TEMP

      RETURN
      END


text:
 dog.love.art  http://www.clairebolding.com
Pets Pet Art Portraits by Users/ClaireBolding Claire Bolding

Email dog.love.art@gmail.com to request a commission

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search