text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!

!
! Test F2008 iso_fortran_env constants/kinds
!
MODULE m
  IMPLICIT NONE
  TYPE dt
    INTEGER, ALLOCATABLE :: i
    CONTAINS
      PROCEDURE :: write_dt
      GENERIC :: WRITE(formatted) => write_dt
  END TYPE

  CONTAINS
    SUBROUTINE write_dt(dtv, unit, iotype, v_list, iostat, iomsg)
      CLASS(dt), INTENT(IN) :: dtv
      INTEGER, INTENT(IN) :: unit
      CHARACTER(*), INTENT(IN) :: iotype
      INTEGER, INTENT(IN) :: v_list(:)
      INTEGER, INTENT(OUT) :: iostat
      CHARACTER(*), INTENT(INOUT) :: iomsg
      INQUIRE(unit, iostat = iostat)
    END SUBROUTINE
END MODULE m

program p
 use, intrinsic :: iso_fortran_env
 use m

  TYPE(dt) d
  CHARACTER(10) :: internal_file
  INTEGER :: iostat

  integer, parameter :: N = 21
  integer :: rslts(N)
  integer :: expect(N) = (/ 1, 2, 4, 8, 1, 2, 4, 8, 4, 8, 1, 2, 4, 8, 1, 2, 4, 8, 4, 8, 99/)

  integer(INT8) :: i1
  integer(INT16) :: i2
  integer(INT32) :: i4
  integer(INT64) :: i8

  logical(LOGICAL8) :: l1
  logical(LOGICAL16) :: l2
  logical(LOGICAL32) :: l4
  logical(LOGICAL64) :: l8

  real(REAL32) :: r32
  real(REAL64) :: r64

  rslts(1) = sizeof(i1)
  rslts(2) = sizeof(i2)
  rslts(3) = sizeof(i4)
  rslts(4) = sizeof(i8)
  rslts(5) = sizeof(l1)
  rslts(6) = sizeof(l2)
  rslts(7) = sizeof(l4)
  rslts(8) = sizeof(l8)
  rslts(9) = sizeof(r32)
  rslts(10) = sizeof(r64)

  rslts(11:14) = INTEGER_KINDS
  rslts(15:18) = LOGICAL_KINDS
  rslts(19:20) = REAL_KINDS(1:2)

  WRITE(internal_file, *, iostat = iostat) d
  rslts(21) = iostat

  call check(rslts, expect, N)
  !print*,rslts
end program

max_stars_repo_name: abrahamtovarmob/flang

max_issues_repo_name: abrahamtovarmob/flang

max_forks_repo_name: abrahamtovarmob/flang


text:
      double complex function a62(st,j1,j2,j3,j4,j5,j6,za,zb)
      implicit none
************************************************************************
*     Author: R.K. Ellis                                               *
*     July, 1998.                                                      *
*     Updated 8/16/01: UV subtraction included                         *
************************************************************************
*     implementation of Eqs. (2.7) and (2.8) of BDKW hep-ph/9610370
*     with ns=0
*     character string st can take the value pp or pm
      include 'constants.f'
      include 'zprods_decl.f'
      integer j1,j2,j3,j4,j5,j6
      character*2 st
      double complex a6,aa6sf,aa6tp,aa6uv

c----Includes ultraviolet subtraction aa6uv (with +ve sign since
c--- it appears in interference terms,e.g. 1234 x 3214)
      call a6routine(st,j1,j2,j3,j4,j5,j6,za,zb,aa6sf,aa6tp,aa6uv)
      a62=a6(st,j1,j2,j3,j4,j5,j6,za,zb)/xnsq
     & +(aa6sf*dble(nf)-aa6tp)/xn+aa6uv

      if (st .eq. 'pp') then
      a62=a62+a6('pm',j1,j3,j2,j4,j5,j6,za,zb)*(one+one/xnsq)
     &       -a6('sl',j2,j3,j1,j4,j5,j6,za,zb)/xnsq
      elseif (st .eq. 'pm') then
      a62=a62+a6('pp',j1,j3,j2,j4,j5,j6,za,zb)*(one+one/xnsq)
     &       +a6('sl',j3,j2,j1,j4,j5,j6,za,zb)/xnsq
      else
      write(6,*) 'Unimplemented st in a62',st
      stop
      endif
      return
      end

ext: f

max_stars_repo_name: tmartini/JHUGen

max_issues_repo_name: tmartini/JHUGen

max_forks_repo_name: tmartini/JHUGen


text:
      SUBROUTINE ADD_PATH( DIR_OUT, FILE_NAME )

!***********************************************************************
!* Adds a Directory Specification to a Path-Free File Name
!*
!* Language: Fortran
!*
!* Author: Stuart G. Mentzer
!*
!* Date: 1999/10/18
!***********************************************************************

      ! Headers
      INCLUDE 'platform.fi'
      INCLUDE 'uds_fxn.fi'

      ! Arguments ______________________________________________________

      CHARACTER*(*) DIR_OUT ! Output directory

      CHARACTER*(*) FILE_NAME ! File name

      ! Variables ______________________________________________________

      INTEGER LD

      CHARACTER FILE_TEMP*512

      ! Add DIR_OUT if FILE_NAME is path-free
      IF ( ( .NOT. BLANK( DIR_OUT ) ) .AND.
     & ( .NOT. ANY_CHARS( FILE_NAME, PATH_CHARS ) ) ) THEN
        FILE_TEMP = DIR_OUT
        LD = LEN_TRIM( DIR_OUT )
        IF ( .NOT. ANY_CHARS( DIR_OUT(LD:LD), PATH_CHARS ) ) THEN
          FILE_TEMP(LEN_TRIM(FILE_TEMP)+1:) = DIR_FILE_SEP
        END IF
        FILE_TEMP(LEN_TRIM(FILE_TEMP)+1:) = FILE_NAME
        FILE_NAME = FILE_TEMP
      END IF

      RETURN
      END

ext: for

max_stars_repo_name: DeadParrot/NHTSA-Tools

max_issues_repo_name: DeadParrot/NHTSA-Tools

max_forks_repo_name: DeadParrot/NHTSA-Tools


text:
      SUBROUTINE ztsrange (IFLTAB, CPATH, KSEARCH, CFIRSTP, CLASTP,
     *                     NFOUND)
C
C
C     Given a time series pathname, this routine will return the
C     pathname of the first and last data set in that series
C     (varying D (date) part).
C
C     The input pathname may or may not have a valid D (date part)
C
C     This routine can do a quick search or an exhaustive search.
C
C     The quick search will find a valid record and search forwards
C     and backwards until up to ISEARCH records are not found
C     (for example, for a monthly block record, this will check
C     records forwards until a check fails.  At that point the
C     date will be incremented for up to ISEARCH values (e.g., 10)
C     checking for the records existence.  If no more are found, the
C     last one that was found is returned.
C     For most normal database files, this is fast and satisfactory.
C
C     The exhaustive search searches every record in the database,
C     then returns the first and last in the series.  This method
C     will always be accurate, but may take considerable time for
C     larger database files.
C
C     Typically data blocks of one day (10 minute interval or less)
C     use the exhaustive search, while larger block use the quick
C     search.
C
C
C     IFLTAB - (input) Array from the zopen6 call.
C     CPATH  - (input) Pathname to start from.  It does not
C              have to have a valid (or known) D part.
C     ISEARCH - (input)  The number of records to check beyond
C              the last one found (typically 10), OR set to
C              zero (0) to do an exhaustive search.
C     CFIRSTP - (output)  The first pathname in the series.
C     CLASTP  - (output)  The last pathname in the series.
C     NFOUND  - (output)  The number of valid (checked) records found.
C              (this is set to -1 on an error)
C
      INTEGER IFLTAB(*)
      CHARACTER CPATH*(*), CFIRSTP*(*), CLASTP*(*)
      INTEGER ISEARCH, NFOUND
C
      LOGICAL LFOUND, LOPEN
      CHARACTER CPATH1*392, CPATH2*392, CPATH3*392
      CHARACTER CLINE*392
      INTEGER IBPART(6), IEPART(6), ILPART(6)
      INTEGER IDX, KSEARCH
C
C     Standard quick search
      ISEARCH = KSEARCH
      CFIRSTP = ' '
      CLASTP = ' '
      IF (ISEARCH.GT.0) THEN
C        Does the record exist?
         CALL CHRLNB (CPATH, NPATH)
         CALL zcheck (IFLTAB, CPATH, NPATH, NHEAD, NDATA, LFOUND)
C
C        If not, search for a valid record
         IF (.NOT.LFOUND) THEN
            CALL ztssrch (IFLTAB, CPATH, 0, CPATH1, NFOUND)
            IF (NFOUND.LE.0) THEN
               ISEARCH = 0
               GO TO 200
            ENDIF
         ELSE
            CPATH1 = CPATH
         ENDIF
C
C        Now loop backwards in time
         CPATH2 = CPATH1
         call strcpy(CFIRSTP, CPATH1)
         CALL CHRLNB (CPATH1, NPATH)
         NFOUND = 1
C
         JSEARCH = 0
 40      CONTINUE
            CALL znextts (IFLTAB, CPATH2, CPATH3, .FALSE., ISTAT)
            IF (ISTAT.NE.0) GO TO 900
            CALL zcheck (IFLTAB, CPATH3, NPATH, NHEAD, NDATA, LFOUND)
            IF (LFOUND) THEN
               NFOUND = NFOUND + 1
               JSEARCH = 0
               call strcpy(CFIRSTP, CPATH3)
            ELSE
               JSEARCH = JSEARCH + 1
               IF (JSEARCH.GE.ISEARCH) GO TO 100
            ENDIF
            CPATH2 = CPATH3
            GO TO 40
C
 100     CONTINUE
C        Search forwards in time
         JSEARCH = 0
         CPATH2 = CPATH1
         call strcpy(CLASTP, CPATH1)
C
 140     CONTINUE
            CALL znextts (IFLTAB, CPATH2, CPATH3, .TRUE., ISTAT)
            IF (ISTAT.NE.0) GO TO 900
            CALL zcheck (IFLTAB, CPATH3, NPATH, NHEAD, NDATA, LFOUND)
            IF (LFOUND) THEN
               NFOUND = NFOUND + 1
               JSEARCH = 0
               call strcpy(CLASTP, CPATH3)
            ELSE
               JSEARCH = JSEARCH + 1
               IF (JSEARCH.GE.ISEARCH) THEN
                  IF (NFOUND.EQ.0) THEN
                     ISEARCH = 0
                     GO TO 200
                  ENDIF
                  GO TO 800
               ENDIF
            ENDIF
            CPATH2 = CPATH3
            GO TO 140
C
C           can not ever get here!
C
      ENDIF

 200  CONTINUE
      IF (ISEARCH.EQ.0) THEN
C        Exhaustive search
C        Try to find a scrach unit to open
         JUNIT = 80
         DO 160 I=1,30
            INQUIRE (UNIT=JUNIT, OPENED=LOPEN)
            IF (.NOT.LOPEN) GO TO 180
            JUNIT = JUNIT + 1
 160     CONTINUE
C        Shouldn't get here
         GO TO 920
C
 180     CONTINUE
         CALL TEMPNAME (CPATH3, JUNIT)
         CALL CHRLNB (CPATH3, N)
         OPEN (UNIT=JUNIT, FILE=CPATH3(1:N), ERR=920)
         CALL ztssrch (IFLTAB, CPATH, JUNIT, CPATH1, NFOUND)
         IF (NFOUND.LE.0) THEN
            CLOSE (UNIT=JUNIT, STATUS='DELETE')
            GO TO 900
         ENDIF
C
         REWIND (UNIT=JUNIT)
         DO 240 I=1,NFOUND
            READ (JUNIT,220,END=300) CLINE
 220        FORMAT (A)
            IDX = INDEX(CLINE, '/')
            IF (IDX.GT.0) THEN
              CPATH1 = CLINE(IDX:)
            ELSE
              CPATH1 = CLINE
            ENDIF
            IF (I.EQ.1) THEN
               CPATH2 = CPATH1
C              Get the location of the "D" part
               CALL zupath (CPATH1, IBPART, IEPART, ILPART, ISTAT)
               IF (ISTAT.NE.0) GO TO 900
               I4TH = IBPART(4)
               I5TH = IEPART(4)
            ENDIF
            CALL DATJUL (CPATH1(I4TH:I5TH), JULIAN, ISTAT)
            IF (ISTAT.NE.0) GO TO 900
            IF (I.EQ.1) THEN
               JULS = JULIAN
               JULE = JULIAN
            ELSE
               IF (JULIAN.LT.JULS) JULS = JULIAN
               IF (JULIAN.GT.JULE) JULE = JULIAN
            ENDIF
 240     CONTINUE
C
 300     CONTINUE
         CLOSE (UNIT=JUNIT, STATUS='DELETE')
         call strcpy(CFIRSTP, CPATH2)
         call strcpy(CLASTP, CPATH2)
         CALL JULDAT (JULS, 104, CFIRSTP(I4TH:I5TH), N)
         CALL JULDAT (JULE, 104, CLASTP(I4TH:I5TH), N)
      ENDIF
C
C
 800  CONTINUE
      RETURN
C
 900  CONTINUE
      NFOUND = -1
      GO TO 800
C
 920  CONTINUE
      NFOUND = -2
      GO TO 800
C
      END

max_stars_repo_name: HydrologicEngineeringCenter/heclib

max_issues_repo_name: HydrologicEngineeringCenter/heclib

max_forks_repo_name: HydrologicEngineeringCenter/heclib


!
!
!                             AMG4PSBLAS version 1.0
!    Algebraic Multigrid Package
!               based on PSBLAS (Parallel Sparse BLAS version 3.7)
!
!    (C) Copyright 2021
!
!        Salvatore Filippone
!        Pasqua D'Ambra
!        Fabio Durastante
!
!    Redistribution and use in source and binary forms, with or without
!    modification, are permitted provided that the following conditions
!    are met:
!      1. Redistributions of source code must retain the above copyright
!         notice, this list of conditions and the following disclaimer.
!      2. Redistributions in binary form must reproduce the above copyright
!         notice, this list of conditions, and the following disclaimer in the
!         documentation and/or other materials provided with the distribution.
!      3. The name of the AMG4PSBLAS group or the names of its contributors may
!         not be used to endorse or promote products derived from this
!         software without specific written permission.
!
!    THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
!    ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT LIMITED
!    TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A PARTICULAR
!    PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE AMG4PSBLAS GROUP OR ITS CONTRIBUTORS
!    BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, SPECIAL, EXEMPLARY, OR
!    CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT LIMITED TO, PROCUREMENT OF
!    SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, DATA, OR PROFITS; OR BUSINESS
!    INTERRUPTION) HOWEVER CAUSED AND ON ANY THEORY OF LIABILITY, WHETHER IN
!    CONTRACT, STRICT LIABILITY, OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE)
!    ARISING IN ANY WAY OUT OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE
!    POSSIBILITY OF SUCH DAMAGE.
!
!
! File: amg_zaggrmat_smth_bld.F90
!
! Subroutine: amg_zaggrmat_smth_bld
! Version:    complex
!
!  This routine builds a coarse-level matrix A_C from a fine-level matrix A
!  by using the Galerkin approach, i.e.
!
!                               A_C = P_C^T A P_C,
!
!  where P_C is a prolongator from the coarse level to the fine one.
!
!  The prolongator P_C is built according to a smoothed aggregation algorithm,
!  i.e. it is obtained by applying a damped Jacobi smoother to the piecewise
!  constant interpolation operator P corresponding to the fine-to-coarse level
!  mapping built by the amg_aggrmap_bld subroutine:
!
!                            P_C = (I - omega*D^(-1)A) * P,
!
!  where D is the diagonal matrix with main diagonal equal to the main diagonal
!  of A, and omega is a suitable smoothing parameter. An estimate of the spectral
!  radius of D^(-1)A, to be used in the computation of omega, is provided,
!  according to the value of p%parms%aggr_omega_alg, specified by the user
!  through amg_zprecinit and amg_zprecset.
!
!  The coarse-level matrix A_C is distributed among the parallel processes or
!  replicated on each of them, according to the value of p%parms%coarse_mat,
!  specified by the user through amg_zprecinit and amg_zprecset.
!  On output from this routine the entries of AC, op_prol, op_restr
!  are still in "global numbering" mode; this is fixed in the calling routine
!  aggregator%mat_bld.
!
!
! Arguments:
!    a          -  type(psb_zspmat_type), input.
!                  The sparse matrix structure containing the local part of
!                  the fine-level matrix.
!    desc_a     -  type(psb_desc_type), input.
!                  The communication descriptor of the fine-level matrix.
!    p          -  type(amg_z_onelev_type), input/output.
!                  The 'one-level' data structure that will contain the local
!                  part of the matrix to be built as well as the information
!                  concerning the prolongator and its transpose.
!    parms      -   type(amg_dml_parms), input
!                  Parameters controlling the choice of algorithm
!    ac         -  type(psb_zspmat_type), output
!                  The coarse matrix on output
!
!    ilaggr     -  integer, dimension(:), input
!                  The mapping between the row indices of the coarse-level
!                  matrix and the row indices of the fine-level matrix.
!                  ilaggr(i)=j means that node i in the adjacency graph
!                  of the fine-level matrix is mapped onto node j in the
!                  adjacency graph of the coarse-level matrix. Note that the indices
!                  are assumed to be shifted so as to make sure the ranges on
!                  the various processes do not   overlap.
!    nlaggr     -  integer, dimension(:) input
!                  nlaggr(i) contains the aggregates held by process i.
!    op_prol    -  type(psb_zspmat_type), input/output
!                  The tentative prolongator on input, the computed prolongator on output
!
!    op_restr    -  type(psb_zspmat_type), output
!                  The restrictor operator; normally, it is the transpose of the prolongator.
!
!    info       -  integer, output.
!                  Error code.
!
subroutine amg_zaggrmat_smth_bld(a,desc_a,ilaggr,nlaggr,parms,&
     & ac,desc_ac,op_prol,op_restr,t_prol,info)
  use psb_base_mod
  use amg_base_prec_type
  use amg_z_inner_mod, amg_protect_name => amg_zaggrmat_smth_bld
  use amg_z_base_aggregator_mod

  implicit none

  ! Arguments
  type(psb_zspmat_type), intent(in)      :: a
  type(psb_desc_type), intent(inout)       :: desc_a
  integer(psb_lpk_), intent(inout)         :: ilaggr(:), nlaggr(:)
  type(amg_dml_parms), intent(inout)      :: parms
  type(psb_zspmat_type), intent(inout)      :: op_prol,ac,op_restr
  type(psb_lzspmat_type), intent(inout)    :: t_prol
  type(psb_desc_type), intent(inout)       :: desc_ac
  integer(psb_ipk_), intent(out)           :: info

  ! Local variables
  integer(psb_lpk_) :: nrow, nglob, ncol, ntaggr, ip, &
       & naggr, nzl,naggrm1,naggrp1, i, j, k, jd, icolF, nrw
  integer(psb_ipk_) :: inaggr, nzlp
  type(psb_ctxt_type) :: ctxt
  integer(psb_ipk_) :: np, me
  character(len=20) :: name
  type(psb_lz_coo_sparse_mat) :: tmpcoo
  type(psb_z_coo_sparse_mat) :: coo_prol, coo_restr
  type(psb_z_csr_sparse_mat) :: acsr1,  acsrf, csr_prol, acsr
  complex(psb_dpk_), allocatable :: adiag(:)
  real(psb_dpk_), allocatable :: arwsum(:)
  integer(psb_ipk_)  :: ierr(5)
  logical            :: filter_mat
  integer(psb_ipk_)            :: debug_level, debug_unit, err_act
  integer(psb_ipk_), parameter :: ncmax=16
  real(psb_dpk_)     :: anorm, omega, tmp, dg, theta
  logical, parameter :: debug_new=.false.
  character(len=80) :: filename

  name='amg_aggrmat_smth_bld'
  info=psb_success_
  call psb_erractionsave(err_act)
  if (psb_errstatus_fatal()) then
    info = psb_err_internal_error_; goto 9999
  end if
  debug_unit  = psb_get_debug_unit()
  debug_level = psb_get_debug_level()

  ctxt = desc_a%get_context()

  call psb_info(ctxt, me, np)

  nglob = desc_a%get_global_rows()
  nrow  = desc_a%get_local_rows()
  ncol  = desc_a%get_local_cols()

  theta = parms%aggr_thresh

  naggr  = nlaggr(me+1)
  ntaggr = sum(nlaggr)

  naggrm1 = sum(nlaggr(1:me))
  naggrp1 = sum(nlaggr(1:me+1))
  filter_mat = (parms%aggr_filter == amg_filter_mat_)

  !
  ! naggr: number of local aggregates
  ! nrow: local rows.
  !

  ! Get the diagonal D
  adiag = a%get_diag(info)
  if (info == psb_success_) &
       & call psb_realloc(ncol,adiag,info)
  if (info == psb_success_) &
       & call psb_halo(adiag,desc_a,info)
  if (info == psb_success_) call a%cp_to(acsr)

  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='sp_getdiag')
    goto 9999
  end if

  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & ' Initial copies done.'

  call acsr%cp_to_fmt(acsrf,info)

  if (filter_mat) then
    !
    ! Build the filtered matrix Af from A
    !

    do i=1, nrow
      tmp = zzero
      jd  = -1
      do j=acsrf%irp(i),acsrf%irp(i+1)-1
        if (acsrf%ja(j) == i) jd = j
        if (abs(acsrf%val(j)) < theta*sqrt(abs(adiag(i)*adiag(acsrf%ja(j))))) then
          tmp=tmp+acsrf%val(j)
          acsrf%val(j)=zzero
        endif

      enddo
      if (jd == -1) then
        write(0,*) 'Wrong input: we need the diagonal!!!!', i
      else
        acsrf%val(jd)=acsrf%val(jd)-tmp
      end if
    enddo
    ! Take out zeroed terms
    call acsrf%clean_zeros(info)
  end if

  do i=1,size(adiag)
    if (adiag(i) /= zzero) then
      adiag(i) = zone / adiag(i)
    else
      adiag(i) = zone
    end if
  end do

  if (parms%aggr_omega_alg == amg_eig_est_) then

    if (parms%aggr_eig == amg_max_norm_) then
      allocate(arwsum(nrow))
      call acsr%arwsum(arwsum)
      anorm = maxval(abs(adiag(1:nrow)*arwsum(1:nrow)))
      call psb_amx(ctxt,anorm)
      omega = 4.d0/(3.d0*anorm)
      parms%aggr_omega_val = omega

    else
      info = psb_err_internal_error_
      call psb_errpush(info,name,a_err='invalid amg_aggr_eig_')
      goto 9999
    end if

  else if (parms%aggr_omega_alg == amg_user_choice_) then

    omega = parms%aggr_omega_val

  else if (parms%aggr_omega_alg /= amg_user_choice_) then
    info = psb_err_internal_error_
    call psb_errpush(info,name,a_err='invalid amg_aggr_omega_alg_')
    goto 9999
  end if

  call acsrf%scal(adiag,info)
  if (info /= psb_success_) goto 9999

  call t_prol%mv_to(tmpcoo)
  inaggr = naggr
  call psb_cdall(ctxt,desc_ac,info,nl=inaggr)
  nzlp = tmpcoo%get_nzeros()
  call desc_ac%indxmap%g2lip_ins(tmpcoo%ja(1:nzlp),info)
  call tmpcoo%set_ncols(desc_ac%get_local_cols())
  call tmpcoo%mv_to_ifmt(csr_prol,info)

  call psb_cdasb(desc_ac,info)
  call psb_cd_reinit(desc_ac,info)
  !
  ! Build the smoothed prolongator using either A or Af
  !    acsr1 = (I-w*D*A) Prol      acsr1 = (I-w*D*Af) Prol
  ! This is always done through the variable acsrf which
  ! is a bit less readable, but saves space and one matrix copy
  !
  call omega_smooth(omega,acsrf)
  call psb_par_spspmm(acsrf,desc_a,csr_prol,acsr1,desc_ac,info)
  if(info /= psb_success_) then
    call psb_errpush(psb_err_from_subroutine_,name,a_err='spspmm 1')
    goto 9999
  end if

  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done SPSPMM 1'
  nzl = acsr1%get_nzeros()
  call acsr1%mv_to_coo(coo_prol,info)

  call amg_ptap_bld(acsr,desc_a,nlaggr,parms,ac,&
       & coo_prol,desc_ac,coo_restr,info)

  call op_prol%mv_from(coo_prol)
  call op_restr%mv_from(coo_restr)

  if (debug_level >= psb_debug_outer_) &
       & write(debug_unit,*) me,' ',trim(name),&
       & 'Done smooth_aggregate '
  call psb_erractionrestore(err_act)
  return

9999 continue
  call psb_errpush(info,name)
  call psb_error_handler(err_act)
  return

contains

  subroutine omega_smooth(omega,acsr)
    implicit none
    real(psb_dpk_),intent(in) :: omega
    type(psb_z_csr_sparse_mat), intent(inout) :: acsr
    !
    integer(psb_lpk_) :: i,j
    do i=1,acsr%get_nrows()
      do j=acsr%irp(i),acsr%irp(i+1)-1
        if (acsr%ja(j) == i) then
          acsr%val(j) = zone - omega*acsr%val(j)
        else
          acsr%val(j) = - omega*acsr%val(j)
        end if
      end do
    end do
  end subroutine omega_smooth

end subroutine amg_zaggrmat_smth_bld


text:
module my_collective

    interface myy_alltoallv
    module procedure imyall2allv
    module procedure rmyall2allv
    module procedure dmyall2allv
    end interface

    integer ncalls,nsend,nget
    double precision t1,t2,btime,ttime
    data ncalls,nsend,nget / 0,      0,     0 /
    data t1,t2,btime,ttime / 0.0d0 , 0.0d0, 0.0d0, 0.0d0 /
    private ncalls,nsend,nget,t1,t2,btime,ttime

contains

    subroutine dump_stat(base)
    implicit none
    character(len=*) , optional :: base
    character (len=80) temp
    integer myid, ierr
    include "mpif.h"
    call MPI_COMM_RANK( MPI_COMM_WORLD, myid, ierr )
    if (present(base) )then
    	write(temp,"(a,i5.5)")trim(base),myid
    else
    	write(temp,"(a,i5.5)")trim("myall2allv"),myid
    endif
    open(58,file=temp,access="append")
    write(58,'(i5,i12,i12,4f10.5)')ncalls,nsend,nget,btime,ttime
    close(58)
    end subroutine

! myall2allv is a replacement for the MPI_ALLTOALLV
! i could not get the mpi version to work correctly
! myall2allv is based on the hypercube algorithm presented
! in class but it does not require a power of 2 processors
subroutine imyall2allv(sendbuf,sendcounts,sdispls,it1,&
                       recvbuf,recvcounts,rdispls,it2,&
                       comm,ierr)
  include "mpif.h"

  integer sendbuf(*),recvbuf(*)
  integer sendcounts(0:*),sdispls(0:*),recvcounts(0:*),rdispls(0:*)
  integer comm,ierr

  integer myid,numnodes,status(MPI_STATUS_SIZE)
  integer i,n2,xchng,my_tag

  call MPI_COMM_RANK( comm, myid, ierr )
  call MPI_COMM_SIZE( comm, numnodes, ierr )
  my_tag=49256
  if ( it1 .ne. MPI_INTEGER .or. it2 .ne. MPI_INTEGER)then
  	write(*,*)"type mismatch in myall2allv",it1,it2," not",MPI_INTEGER
  endif

  ncalls=ncalls+1
  nsend=nsend+sum(sendcounts(0:numnodes-1))
  nget=nget+sum(recvcounts(0:numnodes-1))
  t1=MPI_Wtime()

! find the power of two >= numnodes
  n2=nint(log((1.0d0*numnodes))/log(2.0d0))
  if(2**n2 < numnodes)n2=n2+1
  n2=2**n2
  n2=n2-1
  do i=1,n2
! do xor to find the processor xchng
      xchng=ieor(i,myid)
      if(xchng .le. (numnodes-1))then
          if(myid < xchng)then
          	if ( sendcounts(xchng) .gt. 0)then
              call MPI_SEND(sendbuf(sdispls(xchng)+1),sendcounts(xchng),&
                  MPI_INTEGER,xchng,my_tag,comm,ierr)
            endif
            if ( recvcounts(xchng) .gt. 0) then
              call MPI_RECV(recvbuf(rdispls(xchng)+1),recvcounts(xchng),&
                  MPI_INTEGER,xchng,my_tag,comm,status,ierr)
            endif
          else
            if ( recvcounts(xchng) .gt. 0) then
              call MPI_RECV(recvbuf(rdispls(xchng)+1),recvcounts(xchng),&
                  MPI_INTEGER,xchng,my_tag,comm,status,ierr)
            endif
            if ( sendcounts(xchng) .gt. 0) then
              call MPI_SEND(sendbuf(sdispls(xchng)+1),sendcounts(xchng),&
                  MPI_INTEGER,xchng,my_tag,comm,ierr)
            endif
          endif
      endif
  enddo
! exchange data with myself
  if( sendcounts(myid) .gt. 0)then
  do i=1,sendcounts(myid)
      recvbuf(rdispls(myid)+i)=sendbuf(sdispls(myid)+i)
  enddo
  endif

! gather stats part two
    t2=MPI_Wtime()
    ttime=ttime+(t2-t1)
    t1=t2
    call MPI_Barrier(comm,ierr)
    t2=MPI_Wtime()
    btime=btime+(t2-t1)
! gather stats part two

  return
end subroutine imyall2allv
subroutine rmyall2allv(sendbuf,sendcounts,sdispls,it1,&
                       recvbuf,recvcounts,rdispls,it2,&
                       comm,ierr)
  include "mpif.h"

  real sendbuf(*),recvbuf(*)
  integer sendcounts(0:*),sdispls(0:*),recvcounts(0:*),rdispls(0:*)
  integer comm,ierr

  integer myid,numnodes,status(MPI_STATUS_SIZE)
  integer i,n2,xchng,my_tag

  call MPI_COMM_RANK( comm, myid, ierr )
  call MPI_COMM_SIZE( comm, numnodes, ierr )
  my_tag=49258
  if ( it1 .ne. MPI_REAL .or. it2 .ne. MPI_REAL)then
  	write(*,*)"type mismatch in myall2allv",it1,it2," not",MPI_REAL
  endif

  ncalls=ncalls+1
  nsend=nsend+sum(sendcounts(0:numnodes-1))
  nget=nget+sum(recvcounts(0:numnodes-1))
  t1=MPI_Wtime()

! find the power of two >= numnodes
  n2=nint(log((1.0d0*numnodes))/log(2.0d0))
  if(2**n2 < numnodes)n2=n2+1
  n2=2**n2
  n2=n2-1
  do i=1,n2
! do xor to find the processor xchng
      xchng=ieor(i,myid)
      if(xchng .le. (numnodes-1))then
          if(myid < xchng)then
          	if ( sendcounts(xchng) .gt. 0)then
              call MPI_SEND(sendbuf(sdispls(xchng)+1),sendcounts(xchng),&
                  MPI_REAL,xchng,my_tag,comm,ierr)
            endif
            if ( recvcounts(xchng) .gt. 0) then
              call MPI_RECV(recvbuf(rdispls(xchng)+1),recvcounts(xchng),&
                  MPI_REAL,xchng,my_tag,comm,status,ierr)
            endif
          else
            if ( recvcounts(xchng) .gt. 0) then
              call MPI_RECV(recvbuf(rdispls(xchng)+1),recvcounts(xchng),&
                  MPI_REAL,xchng,my_tag,comm,status,ierr)
            endif
            if ( sendcounts(xchng) .gt. 0) then
              call MPI_SEND(sendbuf(sdispls(xchng)+1),sendcounts(xchng),&
                  MPI_REAL,xchng,my_tag,comm,ierr)
            endif
          endif
      endif
  enddo
! exchange data with myself
  if( sendcounts(myid) .gt. 0)then
  do i=1,sendcounts(myid)
      recvbuf(rdispls(myid)+i)=sendbuf(sdispls(myid)+i)
  enddo
  endif

! gather stats part two
    t2=MPI_Wtime()
    ttime=ttime+(t2-t1)
    t1=t2
    call MPI_Barrier(comm,ierr)
    t2=MPI_Wtime()
    btime=btime+(t2-t1)
! gather stats part two

  return
end subroutine rmyall2allv
subroutine dmyall2allv(sendbuf,sendcounts,sdispls,it1,&
                       recvbuf,recvcounts,rdispls,it2,&
                       comm,ierr)
  include "mpif.h"

  double precision sendbuf(*),recvbuf(*)
  integer sendcounts(0:*),sdispls(0:*),recvcounts(0:*),rdispls(0:*)
  integer comm,ierr

  integer myid,numnodes,status(MPI_STATUS_SIZE)
  integer i,n2,xchng,my_tag

  call MPI_COMM_RANK( comm, myid, ierr )
  call MPI_COMM_SIZE( comm, numnodes, ierr )
  my_tag=49257
  if ( it1 .ne. MPI_REAL8 .or. it2 .ne. MPI_REAL8)then
  	write(*,*)"type mismatch in myall2allv",it1,it2," not",MPI_REAL8
  endif

  ncalls=ncalls+1
  nsend=nsend+sum(sendcounts(0:numnodes-1))
  nget=nget+sum(recvcounts(0:numnodes-1))
  t1=MPI_Wtime()

! find the power of two >= numnodes
  n2=nint(log((1.0d0*numnodes))/log(2.0d0))
  if(2**n2 < numnodes)n2=n2+1
  n2=2**n2
  n2=n2-1
  do i=1,n2
! do xor to find the processor xchng
      xchng=ieor(i,myid)
      if(xchng .le. (numnodes-1))then
          if(myid < xchng)then
          	if ( sendcounts(xchng) .gt. 0)then
              call MPI_SEND(sendbuf(sdispls(xchng)+1),sendcounts(xchng),&
                  MPI_REAL8,xchng,my_tag,comm,ierr)
            endif
            if ( recvcounts(xchng) .gt. 0) then
              call MPI_RECV(recvbuf(rdispls(xchng)+1),recvcounts(xchng),&
                  MPI_REAL8,xchng,my_tag,comm,status,ierr)
            endif
          else
            if ( recvcounts(xchng) .gt. 0) then
              call MPI_RECV(recvbuf(rdispls(xchng)+1),recvcounts(xchng),&
                  MPI_REAL8,xchng,my_tag,comm,status,ierr)
            endif
            if ( sendcounts(xchng) .gt. 0) then
              call MPI_SEND(sendbuf(sdispls(xchng)+1),sendcounts(xchng),&
                  MPI_REAL8,xchng,my_tag,comm,ierr)
            endif
          endif
      endif
  enddo
! exchange data with myself
  if( sendcounts(myid) .gt. 0)then
  do i=1,sendcounts(myid)
      recvbuf(rdispls(myid)+i)=sendbuf(sdispls(myid)+i)
  enddo
  endif

! gather stats part two
    t2=MPI_Wtime()
    ttime=ttime+(t2-t1)
    t1=t2
    call MPI_Barrier(comm,ierr)
    t2=MPI_Wtime()
    btime=btime+(t2-t1)
! gather stats part two

  return
end subroutine dmyall2allv

end module

      module fmpi
      include "mpif.h"
      end module
! This program shows how to use fmpi_Alltoallv.  Each processor
! send/rec a different and random amount of data to/from other
! processors.
! We use fmpi_Alltoall to tell how much data is going to be sent.

module global
    integer numnodes,myid,mpi_err
    integer, parameter :: my_root=0
end module
subroutine init
    use fmpi
    use global
    implicit none
! do the mpi init stuff
    call MPI_INIT( mpi_err )
    call MPI_COMM_SIZE( MPI_COMM_WORLD, numnodes, mpi_err )
    call MPI_Comm_rank(MPI_COMM_WORLD, myid, mpi_err)
end subroutine init

program test1
! poe a.out -procs 3 -rmpool 1
    use fmpi
    use global
    use my_collective
    implicit none
    real, allocatable :: sray(:),rray(:)
    integer, allocatable :: sdisp(:),scounts(:),rdisp(:),rcounts(:)
    integer ssize,rsize,i,k,j
    real z
    integer mytype
    call init
! counts and displacement arrays
    allocate(scounts(0:numnodes-1))
    allocate(rcounts(0:numnodes-1))
    allocate(sdisp(0:numnodes-1))
    allocate(rdisp(0:numnodes-1))
! set counts to send to other processors
! seed the random number generator with a
! different number on each processor
    call seed_random
    do i=0,numnodes-1
        call random_number(z)
        scounts(i)=nint(10.0*z)+1
    enddo
! tell the other processors how much data is coming
    write(*,*)"myid= ",myid," scounts= ",scounts
    call MPI_alltoall(  scounts,1,MPI_INTEGER, &
                        rcounts,1,MPI_INTEGER, &
                        MPI_COMM_WORLD,mpi_err)
    write(*,*)"myid= ",myid," rcounts= ",rcounts
! calculate displacements and the size of the arrays
    sdisp(0)=0
    do i=1,numnodes-1
        sdisp(i)=scounts(i-1)+sdisp(i-1)
    enddo
    rdisp(0)=0
    do i=1,numnodes-1
        rdisp(i)=rcounts(i-1)+rdisp(i-1)
    enddo
    ssize=sum(scounts)
    rsize=sum(rcounts)
! allocate send and rec arrays
    allocate(sray(0:ssize-1))
    allocate(rray(0:rsize-1))
    sray=myid
! send/rec different amounts of data to/from each processor
    mytype=MPI_REAL
    call MPI_alltoallv( sray,scounts,sdisp,mytype, &
                        rray,rcounts,rdisp,mytype, &
                        MPI_COMM_WORLD,mpi_err)
    write(*,*)"myid= ",myid,"    rray= ",rray
    call MYY_Alltoallv( sray,scounts,sdisp,mytype,  &
                        rray,rcounts,rdisp,mytype,  &
                        MPI_COMM_WORLD,mpi_err)

    write(*,*)"myid= ",myid,"    rray= ",rray
    call dump_stat()
    call mpi_finalize(mpi_err)
end program

!typical output
! myid=  0  scounts=  1 5 2
! myid=  0  rcounts=  1 1 1
! myid=  0     rray=  0 1 2

! myid=  1  scounts=  1 7 3
! myid=  1  rcounts=  5 7 5
! myid=  1     rray=  0 0 0 0 0 1 1 1 1 1 1 1 2 2 2 2 2

! myid=  2  scounts=  1 5 10
! myid=  2  rcounts=  2 3 10
! myid=  2     rray=  0 0 1 1 1 2 2 2 2 2 2 2 2 2 2

subroutine seed_random
    use global
    implicit none
    integer the_size,j
    integer, allocatable :: seed(:)
    real z
    call random_seed(size=the_size) ! how big is the intrisic seed?
    allocate(seed(the_size))        ! allocate space for seed
    do j=1,the_size                 ! create the seed
        seed(j)=abs(myid*10)+(j*myid*myid)+100  ! abs is generic
    enddo
    call random_seed(put=seed)      ! assign the seed
    deallocate(seed)
end subroutine

max_stars_repo_name: timkphd/examples

max_issues_repo_name: timkphd/examples

max_forks_repo_name: timkphd/examples


text:
c=======================================================================
c
c     subroutine RNLSSX
c
c     Robust Non-Linear Least Square optimization, smooth solution
c
c------------------------------------------------------------------------
      SUBROUTINE rnlssx ( simext, mxpar, xini, dxmin, df1, epsabs,
     &                    mode, maxitr, maxsim, binf, bsup,
     &                    nobs, moyobs, covobs, skobs, spadim, grid,
     &                    musmo, matsmo,
     &                    liudat, iudata, ldudat, dudata, iwork, dwork,
     &                    xsol, funct, grad,
     &                    totitr, totsim, info )
c------------------------------------------------------------------------
c     INPUT
c            simext : entry point of an external subroutine
c                     provided by the user
c            mxpar  : number of parameters to find               integer
c            xini   : initial solution +Lagrang. vector(mxpar+1)  double
c            dxmin  : precision on  parameters +Lagrang.,
c                     vector(mxpar+1)                             double
c            df1    : estimation of the first decrease of the criter
c                     used if mode=1 computed if =<0
c                     must be nearly abs( f(xvar)-f(xopt) ) )     double
c            epsabs : precision stop test                         double
c            mode   : type of initialization                     integer
c                     = 1 : nothing to initialize
c                     = else : see n2qn1 documentation
c            maxitr : maximum of iterations authorized           integer
c            maxsim : maximum of calls of simul authorized       integer
c                     3*maxitr if largely good
c            binf   : inferior bounds on parameters +Lagrang.,
c                     vector(mxpar+1)                             double
c            bsup   : superior bounds on parameters +Lagrang.,
c                     vector(mxpar+1)                             double
c            nobs   : number of functions (observations)         integer
c            moyobs : mean observations  vector(nobs)             double
c            covobs : covariance matrix of observations
c                     vectorized matrix(nobs*nobs)                double
c            skobs  : stress coefficient of observations          double
c            spadim : space dimension                            integer
c            grid   : grid     vector(nobs*spadim)                double
c            musmo  : smoothing coefficient                       double
c            matsmo : smoothing matrix  vector(mxpar*mxpar)       double
c            liudat : size of integer user data for SIMEXT       integer
c            iudata : integer user data for SIMEXT               integer
c                     vector(liudat)
c            ldudat : size of double prec. user data for SIMEXT  integer
c            dudata : double precision user data for SIMEXT       double
c                     vector(ldudat)
c
c     WORKSPACE
c            iwork  : workspace                                  integer
c                     vector( 2*mxopt + 4 + 12*nobs + liudat )
c            dwork  : double precision  workspace                 double
c                     vector(  mxopt*(mxopt+9)/2
c                            + (mxpar*(mxpar+2))
c                            + nobs*(7*nobs+mxpar+spadim+33)
c                            + 2 + ldudat )
c                     with mxopt = mxpar + 1
c
c     OUTPUT
c            xsol   : solution +Lagrang.    vector(mxpar+1)       double
c            funct  : value of the function at the minimum        double
c            grad   : gradient solution +Lagrang. at the minimum,
c                     vector(mxpar+1)                             double
c            totitr : total of iterations                        integer
c            totsim : total of calls of simul                    integer
c            info   : = 0 successful exit                        integer
c
c     CALL
c        INIRNLS : initial computing of RNLS
c        YV      : copy a vector in a vector
c        IVX     : initialization at a scalar of a vector
c        PMX     : computing M*X = matrix
c                  ( M matrix(n*m), X scalar, gives M*X matrix(n*m) )
c        YVI     : copy a vector of integers in a vector of integers
c        BFGSBOX : quasi-Newton optimizer with BFGS method
c                  Expert version
c        SIMRNLSS: subroutine used by BFGS to compute the value of the
c                  function and its gradient, smooth solution
c
c-----------------------------------------------------------------------
c
c     SIMEXT SIGNATURE
c
c     INPUT
c            nobs   : number of functions                        integer
c            mxopt  : size of the vector solution                integer
c            xopt   : solution  vector(mxopt)                     double
c            spadim : space dimension                            integer
c            grid   : grid     vector(nobs*spadim)                double
c            liudat : size of integers data provided by the user integer
c            iusdat : integers user data provided by the user
c                     vector(liudat)                             integer
c            ldudat : size of doubles data provided by the user  integer
c            dusdat : doubles user data provided by the user
c                     vector(ldudat)                              double
c
c     OUTPUT
c            funct  : functions   vector(nobs)                    double
c            jacob  : jacobian    matrix(nobs*mx)                 double
c            info   : = 0 successful exit                        integer
c
c-----------------------------------------------------------------------
c
      IMPLICIT NONE
c
c     external functions
      EXTERNAL simext, simrnlss, YV, YVI, PMX, inirnls, bfgsbox
c
c     arguments i/o
      INTEGER mxpar, mode, maxitr, maxsim, nobs, spadim, totitr, totsim,
     &        liudat, ldudat, info
      INTEGER iudata(*)
      DOUBLE PRECISION df1, epsabs, skobs, musmo, funct
      DOUBLE PRECISION xini(*), dxmin(*), binf(*), bsup(*), grad(*),
     &                 moyobs(*), grid(*), covobs(nobs,*), xsol(*),
     &                 dudata(*), matsmo(*)
c
c     workspaces
      INTEGER iwork(*)
      DOUBLE PRECISION dwork(*)
c
c     local variables
      INTEGER mxopt
      INTEGER piwbf, pnobs, pspad, pliu, pldu, pinfo, piwu, piwini,
     &        pdwini,
     &        pdwbf, pwsnls, pwsext, pgrid, psk, palpha, pbeta, pvp,
     &        prspec, prcov, psmo, pdwu
      DOUBLE PRECISION eigmax
c
c-----------------------------------------------------------------------
c
c     initializations
      info = 0
      mxopt = mxpar + 1
c
c     pointers for integer work space  : iwork
c
      piwini = 1
c     piwini : pointer for INIRNLS who needs (12*nobs), so 12*nobs more
      piwbf = piwini + 12*nobs
c     piwbf  : pointer for BFGSBOX who needs (2*mxopt + 1),
c             so (2*mxopt + 1) more
      pnobs = piwbf + (2*mxopt + 1)
c     pnobs  : pointer for the number of functions, communication
c              with SIMUL ( simulator of user ), so 1 more
      pspad = pnobs + 1
c     pspad  : pointer for the space dimension, so 1 more
      pliu   = pspad + 1
c     pliu   : pointer for size of integers data provided by the user,
c              ( communication with SIMEXT, user simulator )
      pldu   = pliu + 1
c     pldu   : pointer for size of doubles data provided by the user,
c              ( communication with SIMEXT, user simulator )
      pinfo  = pldu + 1
c     pinfo  : pointer for info (errors of simnls), so 1 more
      piwu = pinfo + 1
c     piwu  : pointer for user SIMEXT workspace, so (liudat) more
c
c     so size of iwork array = ( 2*mxopt + 6 + 12*nobs + liudat )
c
c
c     pointers for double precision work space  : dwork
c
      prcov = 1
c     prcov : pointer for the reverse covariance matrix,
c             so (nobs*nobs) more
      pdwini = prcov + nobs*nobs
c     pdwini: pointer for INIRNLS who needs nobs*(2*nobs+27)
      pdwbf =  pdwini + nobs*(2*nobs+27)
c     pdwbf : pointer for BFGSBOX who needs ( mxopt*(mxopt+9)/2 )
      pwsnls = pdwbf + mxopt*(mxopt+9)/2
c     pwsnls: pointer for internal workspace of SIMRNLSS,
c             so (nobs*(3*nobs+mxpar+3)+2*mxpar) more
      pwsext = pwsnls + (nobs*(3*nobs+mxpar+4)+2*mxpar)
c     pwsext: pointer for workspace ( communication with SIMNLS,
c             and SIMEXT, BFGSBOX simulator ),
c             so ( nobs*(spadim+nobs+2) + 2 + ldudat ) more
      pgrid = pwsext
c     pgrid : pointer for grid matrix ( communication with SIMEXT,
c             user simulator ), so (nobs*spadim) more
      psk = pgrid + (nobs*spadim)
c     psk   : pointer for the stress coefficient k ( communication
c             with SIMRNLSS, user simulator ), so 1 more
      palpha = psk + 1
c     palpha: pointer for the coefficient alpha ( communication
c             with SIMRNLSS, user simulator ), so 1 more
      pbeta = palpha + 1
c     pbeta : pointer for the vector beta ( communication with SIMRNLSS,
c             user simulator ), so (nobs) more
      pvp = pbeta + nobs
c     pvp   : pointer for the eigenvestors matrix ( communication
c             with SIMRNLSS user simulator ), so (nobs*nobs) more
      prspec = pvp + nobs*nobs
c     prspec: pointer for eigenvalues reverse vector ( communication
c             with SIMRNLSS, user simulator ), so nobs more
      psmo = prspec + nobs
c     psmo  : pointer for the smoothing matrix ( communication
c             with SIMNLSS, BFGSBOX simulator ), so (mxpar*mxpar) more
      pdwu = psmo + (mxpar*mxpar)
c     pdwu   : pointer for user SIMEXT workspace, so (ldudat) more
c
c        so size of dwork array = nobs*nobs + nobs*(2*nobs+27)
c                               + mxopt*(mxopt+9)/2
c                               + nobs*(3*nobs+mxpar+4)
c                               + (mxpar*mxpar) + 2*mxpar
c                               + nobs*(spadim+nobs+2) + 2 + ldudat
c        with mxopt = mxpar + 1
c                      = mxopt*(mxopt+11)/2 + (mxpar*(mxpar+2))
c                      + nobs*(7*nobs+mxpar+spadim+33) + 2 + ldudat
c
c-----------------------------------------------------------------------
c
c     initial computing for SIMRNLSS
      CALL inirnls ( nobs, moyobs, covobs, skobs,
     &               iwork(piwini), dwork(pdwini),
     &               dwork(palpha), dwork(pbeta), dwork(pvp),
     &               dwork(prspec), dwork(prcov), eigmax, info )
      IF ( info .NE. 0 ) RETURN
c
c     initialization of solution
      CALL YV ( mxopt, xini, xsol )
c
c     compute smoothing matrix musmo*matsmo
      CALL PMX ( mxpar, mxpar, matsmo, musmo, dwork(psmo) )
c
c     communication with SIMEXT (user simulator)
      iwork(pnobs) = nobs
      iwork(pspad) = spadim
      iwork(pliu) = liudat
      iwork(pldu) = ldudat
      iwork(pinfo)= 0
      CALL YV ( nobs*spadim, grid, dwork(pgrid) )
      CALL YV ( ldudat, dudata, dwork(pdwu) )
      CALL YVI ( liudat, iudata, iwork(piwu) )
      dwork(psk)= skobs
c
c     optimization BFGS
      CALL bfgsbox ( simrnlss, simext, mxopt, xsol,
     &               dxmin, df1, epsabs, mode, maxitr, maxsim,
     &               binf, bsup, iwork(piwbf), dwork(pdwbf),
     &               xsol, funct, grad, totitr, totsim, info)
      IF (iwork(pinfo) .NE. 0) THEN
         info = iwork(pinfo)
         RETURN
      ENDIF
c
      RETURN
      END

max_stars_repo_name: alpgodev/numx

max_issues_repo_name: alpgodev/numx

max_forks_repo_name: alpgodev/numx


text:
!
! Part of the LLVM Project, under the Apache License v2.0 with LLVM Exceptions.
! See https://llvm.org/LICENSE.txt for license information.
! SPDX-License-Identifier: Apache-2.0 WITH LLVM-exception
!

!          THIS CODE AND INFORMATION ARE PROVIDED "AS IS" WITHOUT
!   WARRANTY OF ANY KIND, EITHER EXPRESSED OR IMPLIED, INCLUDING BUT
!   NOT LIMITED TO THE IMPLIED WARRANTIES OF MERCHANTABILITY AND/OR
!   FITNESS FOR A PARTICULAR PURPOSE.
!

! iso_fortran_env.f90
! 32/64 bit  linux and windows.  Add further targets as required.

        module ISO_FORTRAN_ENV

	public

	integer CHARACTER_STORAGE_SIZE
	parameter (CHARACTER_STORAGE_SIZE = 8)
	integer ERROR_UNIT
	parameter (ERROR_UNIT = 0)
	integer FILE_STORAGE_SIZE
	parameter (FILE_STORAGE_SIZE = 8)
	integer INPUT_UNIT
	parameter (INPUT_UNIT = 5)
	integer IOSTAT_END
	parameter (IOSTAT_END = -1)
	integer IOSTAT_EOR
	parameter (IOSTAT_EOR = -2)
	integer NUMERIC_STORAGE_SIZE
	parameter (NUMERIC_STORAGE_SIZE = 32)
	integer OUTPUT_UNIT
	parameter (OUTPUT_UNIT = 6)

        integer IOSTAT_INQUIRE_INTERNAL_UNIT
        parameter (IOSTAT_INQUIRE_INTERNAL_UNIT=99)

	integer INT8
	parameter (INT8 = 1)
	integer INT16
	parameter (INT16 = 2)
	integer INT32
	parameter (INT32 = 4)
	integer INT64
	parameter (INT64 = 8)
	integer LOGICAL8
	parameter (LOGICAL8 = 1)
	integer LOGICAL16
	parameter (LOGICAL16 = 2)
	integer LOGICAL32
	parameter (LOGICAL32 = 4)
	integer LOGICAL64
	parameter (LOGICAL64 = 8)
	integer REAL32
	parameter (REAL32 = 4)
	integer REAL64
	parameter (REAL64 = 8)
	integer REAL128
	parameter (REAL128 = 16)

        integer INTEGER_KINDS(4)
        parameter (INTEGER_KINDS = (/INT8, INT16, INT32, INT64/))
        integer LOGICAL_KINDS(4)
        parameter (LOGICAL_KINDS = (/LOGICAL8, LOGICAL16, LOGICAL32, LOGICAL64/))
        integer REAL_KINDS(3)
        parameter (REAL_KINDS = (/REAL32, REAL64, REAL128/))

        end module  ISO_FORTRAN_ENV

max_stars_repo_name: vigbalu/flang

max_issues_repo_name: vigbalu/flang

max_forks_repo_name: vigbalu/flang


text:
program sin_03
use iso_fortran_env, only: dp => real64
implicit none
real(dp) :: x
x = sin(1.5_dp)
print *, x
end program

max_stars_repo_name: Thirumalai-Shaktivel/lfortran

max_issues_repo_name: Thirumalai-Shaktivel/lfortran

max_forks_repo_name: Thirumalai-Shaktivel/lfortran


program main
    use Spoof
    implicit none

    type(CSRMatrix)         :: csr_matrix
    type(DiagonalMatrix)    :: dia_matrix
    double precision        :: dense(10,10)
    double precision        :: vec(10)

    dense = 42.0D+00
    vec = 10.0D+00

    ! Create the matrices
    csr_matrix = CSRMatrix(dense)
    dia_matrix = DiagonalMatrix(dense)

    ! Multiply them by vec
    vec = csr_matrix%multiply(vec)              ! vec is filled with 4200
    vec = dia_matrix%multiply(vec)              ! vec is filled with 176400

    print *, sum(vec)

end program


text:
! Copyright (c) 2011, 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.
!

module a
  implicit none
  type a_t
  logical :: result = .false.
   contains
     procedure :: say_hello
  end type a_t
  private
  public :: a_t
contains
  subroutine say_hello (this)
    class(a_t), intent(inout) :: this
    print *,'Hello from a'
    this%result = .true.
  end subroutine say_hello
end module a

module b
  use a
  implicit none
  type b_t
     type(a_t) :: a
   contains
     procedure :: say_hello
  end type b_t
contains
  subroutine say_hello (this)
    class(b_t), intent(inout) :: this
    call this%a%say_hello()
  end subroutine
end module b

program p
USE CHECK_MOD
   use b
   logical results(1)
   logical expect(1)
   type(b_t) :: bt
   results = .false.
   expect = .true.
   call bt%say_hello()
   results(1) = bt%a%result
   call check(results,expect,1)
end

max_stars_repo_name: kammerdienerb/flang

max_issues_repo_name: kammerdienerb/flang

max_forks_repo_name: kammerdienerb/flang


text:
! illegal-forward-ref-type.f90
! An illegal forward reference to a derived type causes the front end to crash.
! Specifically, the Java VM detects a SIGSEGV in libstdc++ in '__dynamic_cast'.

program p

  type(t2) :: y  ! illegal forward reference

  type :: t2
    real ::  x
  end type

end program

max_stars_repo_name: maurizioabba/rose

max_issues_repo_name: sujankh/rose-matlab

max_forks_repo_name: sujankh/rose-matlab


text:
Amazing Tans is a chain retailer that provides Beauty Salons and Barber Shops tanning services. They have a variety of tanning booth styles including facial tanning beds and stand up options. Amazing Tans caters to any skin type and provides Sun Protection sunless tanning as well as UV. An array of tanning lotions are available for purchase. Services can be purchased as part of package deals, or on a membership basis.

  To learn about other businesses in Davis that provide professional aesthetic services, visit our Beauty Salons and Barber Shops page.

20060310 18:29:34 nbsp a little on the pricy side, but its definitely worth it. i have been tanning for years in my home town and usually start out at 10 minutes my first time (im really fair skinned). well, today i went for the first time to amazing tans and they suggested i start out at 6 minutes. at first i was a little skeptical since i have always started at 10 and never really saw any difference. wow was i wrong! their beds are amazing. i have never gotten so much color in a first session. the girl at the front counter was also really friendly and helpful. i would recommend them! Users/CorinneEpley

20060409 17:16:16 nbsp This place is really nice.  They have monthly plans and tan packages.  In the summer, they sometimes have a deal where you just pay a flat amount and can tan (at level 1) unlimitedly until the end of they year.  The people working there are usually very friendly and informative.  I would definitely recommend this place. Users/CatLa

20070316 13:47:57 nbsp Amazing tans now has a sign in the window asking patrons to notify the manager if the store is closed during posted business hours.  Anyone know if this is in response to recent problems or complaints? Users/DavidGrundler
  Not in response to anything in particular.  Just a standard business sign.

20070405 13:22:01 nbsp WOW!!! The girls there know their stuff.  I never knew there was so much to learn about tanning beds, lotions etc.  Great customer service! A bit pricey, but I guess you get what you pay for!  Users/Jessica

20080308 09:55:06 nbsp Great customer service. On the days I just wanted to get in the tanning beds they did just that. On the days I had questions, they answered them better than I expected. They were very helpful, but it was the most I ever paid to tanning. If you are pale, it sucks because its cheaper if you pay for a certain amount of tans. So someone who tans for 15 minutes pays the same price as someone tanning for 5. I got my base before, so it wasnt too bad Users/lv

20080428 10:26:29 nbsp The place is pretty cool, and they treat dudes fine (I am not the only one I see there). You have to really work their deals though, they are extremly tricky and will NOT give you any kind of price list (cell phone camera anyone? Its the only way. Ive talked to their corporate office and its the same runaround). As with all tanning salons, its really only worth it to go with the most expensive beds, otherwise you will just get burned. Also, DONT get pressured into buying their lotions! What costs 100 bucks there is less than 40 on amazon, 25 on eBay, and so on. All in all though I cant complain, it is very convenient and my GF is not complaining either.  Users/YesItsMe
  20090721 18:35:10 nbsp I just went in to start tanning.  The girl showed me all the beds and how long it would take to get a base tan, etc.  They had a price list on the counter. I didnt get tricked into anything.  She was very straightforward and VERY nice.  No sales pressure at all. Glad I chose to go here.  Users/EricaMacGregor

20080516 13:57:18 nbsp Girls here are far nicer than at the other tanning place in Davis and more helpful. Wait is never too long for a bed. Yes, it is expensive, but worth it.  Users/tgdavis

20080928 13:10:42 nbsp not really sure why people saw this place is expensive!! it is far more reasonable than any other place in town! 5 bucks for 1 tan!!!! cant get better then that and memberships at 20$ a month unlimited...what is expensive?!?! love the employees, love the services! Users/sara87

20081031 14:50:11 nbsp I love the girls at amazing tans... great service, great beds, they arent too pushy, and overall a great tanning salon... Users/meaghan

20090309 00:02:09 nbsp Love this place! In response to the note on lotions, take a look at the ingredients on the back, you get what you pay for! When youre drying out your skin and aging yourself in the tanning process, dont skimp out on lotion. It will help save your tan and keep your skin soft in the process! Many of the lotions also have a ton of bronzers (DHA)! Save your skin people! Get a great lotion! Im loving mine!:) Users/allsmiles

20091011 10:54:23 nbsp Amazing Tans is great! I love this place. I believe they have the best prices in town, and it may be one of their guarantees to beat any other tanning salon prices (dont quote me on that). The girls that work there are all super knowledgeable, friendly and nice (with the exception with one girl I had to deal with.) If you plan on tanning at least once a week, I recommend getting a VIP membership. I pay $24.88 a month for unlimited tans in a level 1 bed and upgrades to higher level beds are cheap! The lotions are super expensive but they do have sales on them occasionally and the bottles are huge and will last for months. Friendly staff, relaxing atmosphere, best place to tan!  Users/MillaMi

20091028 15:12:57 nbsp Was about to go to prom last year...went here to get a spray on tan. Have had great first time experiences at other places, so I wasnt worried. The tan turned out awful, terribly streaky, and there were drips down my legs. Same experience for my two friends. We went back to try and get our money back, they wouldnt return it and told us that it happens every first time, and we had to come back multiple times to have it look good. Thanks for telling us! Wouldnt recommend the spray tans here.  Users/KB@CSUMB

20110129 12:35:10 nbsp AVOID AT ALL COSTS!!! Everything is nice at first, besides being wayyy overpriced, but attempting to end a membership is IMPOSSIBLE. First, there is a rule that you have to come in and quit in writing by signing a form. Then the policy is to charge you for a full month after you originally want to quit. However, I was suckered into the Stay and Save Plan, the manager told me if I signed a paper and paid $5 a month, then if I ever wanted to come back I wouldnt have to pay a new sign up fee (another one of their millions of fees). She told me there was no downside and I could end the Stay and Save program at any time. So I signed it, paid my last month fee (a month more than I wanted to stay, because of the policy) and paid the $5 a month for about 5 months. After this time, money was tight and I wanted out all together. To my despair, I found out that in the fine print it says that to get off the Stay and Save Plan, you have to pay another full price month (about $70). If I knew that I would not have consented to the Stay and Save plan. There was a change in management and the new manager said that even though the old manager had LIED to me and said that there would be no additional fees and that my last month was paid for back when I first signed the Stay and Save plan, she told me there is nothing she can do. I swallowed the bullet  resigned to pay. (and I promised myself I would never sign anything again without reading every word and never trust anyone who gives me their word on what it says) So I went in to try to pay and quit but the employee there told me the manager was not in and that she had taken my forms home with her. Long story short, I am going in today for the 3rd time to try to pay and quit, praying I will be successful. By the way the manager left me a voicemail that I have to make it in by 3 pm today because that is what time she leaves. So I have to adjust my schedule and make it in or I will be locked into paying Amazing Tans huge amounts of money forever. (So considerate, accommodating, and great customer service policy)
DO NOT GO HERE UNLESS YOU LIKE TO BE LIED TO AND SCAMMED. Users/EllieR

20110401 13:44:07 nbsp I LOVE Amazing Tans Fluffy Donuts!!

Definitely the best donuts Ive ever had! Users/OrofinJackson

20131009 18:51:44 nbsp The location on 5th Street is closed as of 9/27/13.   Users/CalamityJanie

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
*----------------------------------------------------------------------*
*                          PHENOLOGY1                                  *
*                          **********                                  *
*                                                                      *
* Used for grasses and crops. Lai being controlled by the available    *
* npp storage.                                                         *
*                                                                      *
*----------------------------------------------------------------------*
      SUBROUTINE PHENOLOGY1(ftagh,ftdth,
     &                      bbm,bb0,bbmax,bblim,ssm,sss,sslim,
     &nppstore,leafnpp,stemnpp,rootnpp,tran,rlai,lai,lairat,rem,leafls,
     &stemls,rootls,leafmol,respref,soil2g,half,tmem,leafv,stemv,rootv,
     &daysoff,laimax,leaflit,stemlit,rootlit,mnth,day,s_ln,s_sr,s_sn,
     &s_rr,s_rn,bb,ss,bbgs,dsbb,nppstorx,nppstor2,daynpp,
     &maxlai,wtfc,yield,resp,sm_trig,suma,tsumam,stemfr,lmor_sc,
     &chill,dschill,s_lr)
*----------------------------------------------------------------------*
      IMPLICIT NONE
      REAL*8 ftagh,nppstore,rootnpp,tran,rlai,lairat,leafmol,respref
      REAL*8 soil2g,half,rem,leaflit,leafv(3600),stemv(1000),rootv(1000)
      REAL*8 leafnpp,stemnpp,laiinc,tmem(200),resp,nppstor2,s_ln,s_mr
      REAL*8 daysoff,laimax,stemnppr,rootnppr,nppstorx,wtfc,s_lr
      REAL*8 s_sn,s_rr,s_rn,stemlit,rootlit,daynpp,maxlai,bb0,bbmax,s_sr
      REAL*8 bblim,sslim,bbsum,yield,sm_trig(30),smtrig,xdaynpp
      REAL*8 suma(360),tsuma,tsumam,stemfr,maint,yy,lmor_sc(3600)
      INTEGER lai,leafls,stemls,rootls,mnth,day,ij,i,bb,gs,bbgs,sssum
      INTEGER bb2bbmin,bb2bbmax,bbm,ssm,sss,ss,dsbb
      INTEGER ftdth,harvest,chill,dschill
      INCLUDE 'param.inc'
      real sumrr,sumsr,sumlr,summr
      save sumrr,sumsr,sumlr,summr

      leaflit = 0.0d0
      stemlit = 0.0d0
      rootlit = 0.0d0
      yield = 0.0d0

      bb2bbmin = 285
      bb2bbmax = 390
      gs = 60

      rem = rlai - int(rlai)
      lai = int(rlai) + 1
      resp = 0.0d0

      harvest = 0

      ss = 0

      xdaynpp = daynpp

*----------------------------------------------------------------------*
* Compute active root and stem npp
*----------------------------------------------------------------------*
      stemnppr = 0.0d0
      rootnppr = 0.0d0
      DO ij=1,stemls
        stemnppr = stemnppr + stemv(ij)
      ENDDO
      DO ij=1,rootls
        rootnppr = rootnppr + rootv(ij)
      ENDDO

*----------------------------------------------------------------------*
* Pay for days roots if veg exists
*----------------------------------------------------------------------*
*      IF (tran/13.4/12.0d0*p_rootsc.LT.0.001d0*nppstore) THEN
*        yy = tran/13.4/12.0d0*p_rootsc
*      ELSE
*        yy = 0.001d0*nppstore
*      ENDIF
*      IF (respref.LT.0.00001d0) THEN
*        yy = 0.0d0
*      ENDIF

      IF ((nppstore.GT.0.0d0).AND.(daynpp.GT.0.0d0)) THEN
        yy = nppstore*p_rootfr
      ELSE
        yy = 0.0d0
      ENDIF

      nppstore = nppstore - yy
      xdaynpp = xdaynpp - yy

      rootnpp = rootnpp + yy
      s_rn = s_rn + yy

*----------------------------------------------------------------------*
* Age roots and set daily vector
*----------------------------------------------------------------------*
      DO ij=1,rootls-1
        rootv(rootls+1-ij) = rootv(rootls-ij)
      ENDDO
      rootv(1) = yy
      IF (rootv(1).LT.0.0d0) THEN
        print*,'rootv negativ: rootv=',rootv(1),'tran=',tran
      ENDIF

*----------------------------------------------------------------------*
* Root resperation
*----------------------------------------------------------------------*
      DO ij=1,rootls-1
        resp = resp + respref*rootv(rootls-ij)
        s_rr = s_rr + respref*rootv(rootls-ij)
        rootnpp = rootnpp - respref*rootv(rootls-ij)
        rootv(rootls-ij) = rootv(rootls-ij)*(1.0 - respref)
      ENDDO

*----------------------------------------------------------------------*
* Check for chilling.                                                  *
*----------------------------------------------------------------------*
      IF (chill.EQ.0) THEN
        bbsum = 0.0d0
        DO i=1,20
          IF (tmem(i).LT.-5.0d0)  bbsum = bbsum + max(-10.0d0,tmem(i)+
     &5.0d0)
        ENDDO
        IF (bbsum.LT.-100) THEN
          chill = 1
          dschill = 1
        ENDIF
      ENDIF
      IF (chill.EQ.1) THEN
        dschill = dschill + 1
      ENDIF
      IF (dschill.GT.260) THEN
        chill = 0
        dschill = 0
      ENDIF

*      print*,mnth,day,bbsum,chill,dschill,' chill'

*----------------------------------------------------------------------*
* Bubburst, if no budburst set, and sufficient soil moisture, then     *
* check for  budburst.                                                 *
*----------------------------------------------------------------------*
      IF (((bb.eq.0).AND.(soil2g.GT.half+0.25d0*(wtfc-half))).OR.
     &((dsbb.GT.bb2bbmax).AND.(soil2g.GT.half+0.1d0*(wtfc-half)))) THEN

        smtrig = 0.0d0
        DO i=1,30
          smtrig = smtrig + sm_trig(i)
        ENDDO

*        print*,mnth,day,smtrig
        IF ((smtrig.GT.30.0d0).OR.(dsbb.GT.bb2bbmax)) THEN

*----------------------------------------------------------------------*
* Check for budburst using degree days.                                *
*----------------------------------------------------------------------*
          bbsum = 0.0d0
          DO i=1,bbm
            IF (tmem(i).GT.bb0)  bbsum = bbsum + min(bbmax,tmem(i)-bb0)
          ENDDO
*          print*,mnth,day,bbsum,bblim,dschill

          IF ((real(bbsum).GE.real(bblim)*exp(-0.01d0*real(dschill)))
     &.OR.(dsbb.GT.bb2bbmax)) THEN
*----------------------------------------------------------------------*
* Adjust proportion of gpp going into stem production based on suma.   *
* This is essentially the LAI control.                                 *
*----------------------------------------------------------------------*
            tsuma = 0.0d0
            DO i=1,360
              tsuma = tsuma + suma(i)
            ENDDO
*            print*,'tsuma ',tsuma,stemfr,nppstore
            maint = max(1.0d0,(real(leafls)/360.0d0))
            tsuma = tsuma - leafmol*1.25d0/maint*p_opt
            IF (tsuma.GT.leafmol*1.25d0/maint*p_opt)
     &tsuma = leafmol*1.25d0/maint*p_opt
            IF (tsuma.LT.-leafmol*1.25d0/maint*p_opt)
     &tsuma = -leafmol*1.25d0/maint*p_opt
            tsumam = 0.0d0*tsumam + 1.0d0*tsuma
            stemfr = stemfr + tsumam*p_laimem
            IF (stemfr.LT.10.0d0) stemfr = 10.0d0

*----------------------------------------------------------------------*
* Bud burst occurance
*----------------------------------------------------------------------*
            bb = (mnth-1)*30 + day
            bbgs = 0
            dsbb = 0
            IF (stemfr.LT.0.8d0*nppstore) THEN
              nppstor2 = nppstore - stemfr
            ELSE
              IF (stemfr.LT.0.75d0*nppstore) THEN
                nppstor2 = nppstore - stemfr
              ELSE
                nppstor2 = nppstore*0.25d0
              ENDIF
              stemfr = stemfr*0.8d0
            ENDIF
            laiinc = (nppstore - 0.0*nppstor2)/leafmol/1.25d0
            nppstorx = nppstore
*            print*,'laiinc ',laiinc,nppstore*12.0
          ENDIF
        ENDIF
      ENDIF

*----------------------------------------------------------------------*
* Compute length of growing season, and set to zero when equal to      *
* growing season.                                                      *
*----------------------------------------------------------------------*
      IF (bb.GT.0)  bbgs = bbgs + 1
      IF (bbgs-gs.gt.bb2bbmin) THEN
        bb = 0
        bbgs = 0
      ENDIF

      IF (dsbb.LT.500) dsbb = dsbb + 1

*----------------------------------------------------------------------*
* Set LAI increase.                                                    *
*----------------------------------------------------------------------*
      IF ((bb.GT.0).AND.(bbgs.LT.gs).AND.(nppstore.GT.0.0d0)) THEN
        laiinc = lairat*(nppstorx - 0.0d0*nppstor2)/leafmol/1.0d0
*        IF (rlai+laiinc.LT.0.5d0)  laiinc = 0.5d0 - rlai
        IF (rlai+laiinc.GT.maxlai)  laiinc = maxlai - rlai
        IF (rlai+laiinc.GT.11.5d0)  laiinc = 11.5d0 - rlai
        IF ((rlai.GT.0).and.(nppstore.LT.0.0d0)) laiinc = 0.0d0
      ELSE
        laiinc = 0.0d0
      ENDIF

*----------------------------------------------------------------------*
* Senescence, if rlai is greater than zero, compute senescence.        *
*----------------------------------------------------------------------*
      IF (rlai.GT.1.0e-6) THEN
        IF (abs(ftagh).LT.1.0E-6) THEN
*----------------------------------------------------------------------*
* Grass senescence
*----------------------------------------------------------------------*
          IF (soil2g.LT.half*0.5d0) THEN
            laiinc = -rlai
            ss = day + (mnth - 1)*30
          ELSEIF (bbgs.GT.100) THEN
*----------------------------------------------------------------------*
* Check for senescence, senescence occurs there are 'sss' days colder  *
* than 'sslim' out of the last 'ssm' days.                             *
*----------------------------------------------------------------------*
            sssum = 0
            DO i=1,ssm
              IF (tmem(i).LT.sslim)  sssum = sssum + 1
            ENDDO
            IF (sssum.GE.sss) THEN
*----------------------------------------------------------------------*
* Senescence event due to temperature.                                 *
*----------------------------------------------------------------------*
              laiinc =-rlai
              ss = day + (mnth - 1)*30
            ENDIF
          ENDIF

        ELSE
*----------------------------------------------------------------------*
* Crop senescence
*----------------------------------------------------------------------*
          IF ((bbgs.GT.leafls).OR.(bbgs.GT.ftdth)) THEN
            harvest = 1
            laiinc =-rlai
            ss = day + (mnth - 1)*30
          ENDIF

        ENDIF

      ENDIF

*----------------------------------------------------------------------*
* Pay for new leaves.
*----------------------------------------------------------------------*
      IF (laiinc.GT.0.0d0) THEN
        nppstore = nppstore - laiinc*leafmol*1.25d0
        nppstorx = nppstorx - laiinc*leafmol*1.25d0
        leafnpp = leafnpp + laiinc*leafmol
        resp = resp + 0.25d0*laiinc*leafmol
        xdaynpp = xdaynpp - laiinc*leafmol*1.25d0

        s_ln = s_ln + laiinc*leafmol*1.25d0
*        print*,rlai,laiinc,s_ln*12.0
      ENDIF
      s_mr = 0.25d0*laiinc*leafmol

*----------------------------------------------------------------------*
* Age leaves by one day, kill any which have died of old age,
* adjust by laiinc (+ or -), and then sum to get 'rlai'.
* 'leafls' is an integer variable of leaf lifespan in days.
*----------------------------------------------------------------------*
      CALL LAIALT(leafv,rlai,laiinc,leafls,leaflit)

      IF (harvest.EQ.1) THEN
        yield = leaflit*ftagh
        leaflit = (1.0d0 - ftagh)*leaflit
      ELSE
        yield = 0.0d0
      ENDIF

      rem = rlai - int(rlai)
      lai = int(rlai) + 1

*----------------------------------------------------------------------*
* leaf death not through age mortality
*----------------------------------------------------------------------*
*      DO i=1,leafls
*        rlai = rlai - leafv(i)*(1.0d0 - lmor**(1.0d0/dble(leafls)))
*        leaflit = leaflit + leafv(i)*(1.0d0 -
*     &lmor**(1.0d0/dble(leafls)))
*        leafv(i)=leafv(i)*lmor**(1.0d0/real(leafls))
*      ENDDO
*----------------------------------------------------------------------*
* leaf death not through age mortality
*----------------------------------------------------------------------*
*      DO i=1,leafls
*        xx = (real(leafls - i)/real(leafls))**lmor
*        rlai = rlai - leafv(i)*(1.0d0 - xx)
*        leaflit = leaflit + leafv(i)*(1.0d0 - xx)
*        leafv(i)=leafv(i)*xx
*      ENDDO
*----------------------------------------------------------------------*
* leaf death not through age mortality
*----------------------------------------------------------------------*
      rlai = rlai + leaflit
      DO i=1,leafls
        IF (leafv(i).GT.0.0d0) THEN
          IF (leafv(i).GT.1.0e-6) THEN
            leaflit = leaflit + leafv(i)
            leafv(i)=leafv(i)*lmor_sc(i)
            leaflit = leaflit - leafv(i)
          ELSE
            leaflit = leaflit + leafv(i)
            leafv(i) = 0.0d0
          ENDIF
        ENDIF
      ENDDO
      rlai = rlai - leaflit

*----------------------------------------------------------------------*
* Stem resperation, and NPP                                            *
*----------------------------------------------------------------------*
*      nppstore = nppstore - stemnppr*respref
*      resp = resp + stemnppr*respref
*      s_sr = s_sr + stemnppr*respref
*      xdaynpp = xdaynpp - stemnppr*respref

      s_sr = resp
      DO ij=1,stemls-1
        resp = resp + respref*stemv(stemls-ij)
        stemnpp = stemnpp - respref*stemv(stemls-ij)
        stemv(stemls-ij) = stemv(stemls-ij)*(1.0 - respref)
      ENDDO
      s_sr = resp - s_sr

      DO ij=1,stemls-1
        stemv(stemls+1-ij) = stemv(stemls-ij)
      ENDDO

      IF ((nppstore.GT.0.0d0).and.(daynpp.GT.0.0d0)) THEN
        yy = nppstore*p_stemfr
        stemv(1) = yy
        s_sn = s_sn + yy
        stemnpp = stemnpp + yy
        nppstore = nppstore - yy
      ELSE
        stemv(1) = 0.0d0
      ENDIF

      IF (stemnppr.GT.0d0) THEN
        IF ((stemv(1).LT.p_stmin)) then
          stemnpp = stemnpp + p_stmin - stemv(1)
          nppstore = nppstore - p_stmin + stemv(1)
          s_sn = s_sn + p_stmin - stemv(1)
          stemv(1) = p_stmin
        ENDIF
      ENDIF
*----------------------------------------------------------------------*

      IF (rlai.GT.laimax)  laimax = rlai
      IF (.NOT.(rlai.GT.0.0d0))  daysoff = daysoff + 1.0d0

      IF (nppstore.LT.0.0d0) THEN
        DO i=1,leafls
          leaflit = leaflit + leafv(i)
          leafv(i) = 0.0d0
        ENDDO
c added by Ghislain 07/10/03
        stemlit = 0.0d0
        DO i=1,stemls
          stemlit = stemlit + stemv(i)
          stemv(i) = 0.0d0
        ENDDO
c added by Ghislain 07/10/03
        rootlit = 0.0d0
        DO i=1,rootls
          rootlit = rootlit + rootv(i)
          rootv(i) = 0.0d0
        ENDDO

        nppstore = 0.0d0
        rlai = 0.0d0

      ELSE
        stemlit = 0.0d0
        rootlit = 0.0d0
      ENDIF

      if ((day == 1).and.(mnth == 1)) then
        sumsr=0.0
        sumrr=0.0
        sumlr=0.0
        summr=0.0
      endif
      sumsr = sumsr + 12.0*s_sr
      sumrr = sumrr + 12.0*s_rr
      sumlr = sumlr + 12.0*s_lr
      summr = summr + 12.0*s_mr
      if ((day == -30).and.(mnth == 12))
     &print'(4f10.3)',sumsr,sumrr,sumlr,summr

      RETURN
      END

*----------------------------------------------------------------------*
*                          PHENOLOGY2                                  *
*                          **********                                  *
*                                                                      *
* Used for trees. Lai being controlled by a weighted memory of npp     *
* storage over the previous 5 years.                                   *
*                                                                      *
*----------------------------------------------------------------------*
      SUBROUTINE PHENOLOGY2(bbm,bb0,bbmax,bblim,ssm,sss,sslim,
     &nppstore,leafnpp,stemnpp,rootnpp,tran,rlai,lai,lairat,rem,leafls,
     &stemls,rootls,leafmol,respref,soil2g,half,tmem,leafv,stemv,rootv,
     &daysoff,laimax,leaflit,stemlit,rootlit,mnth,day,s_ln,s_sr,s_sn,
     &s_rr,s_rn,bb,ss,bbgs,dsbb,nppstorx,nppstor2,daynpp,
     &maxlai,wtfc,yield,resp,sm_trig,suma,tsumam,stemfr,lmor_sc,
     &chill,dschill,s_lr)
*----------------------------------------------------------------------*
      IMPLICIT NONE
      REAL*8 nppstore,rootnpp,tran,rlai,lairat,leafmol,respref,soil2g
      REAL*8 half,rem,leaflit,leafv(3600),stemv(1000),rootv(1000)
      REAL*8 leafnpp,stemnpp,laiinc,tmem(200),resp,nppstor2,s_ln
      REAL*8 daysoff,laimax,stemnppr,rootnppr,nppstorx,wtfc,s_lr,s_mr
      REAL*8 s_sn,s_rr,s_rn,stemlit,rootlit,daynpp,maxlai,bb0,bbmax,s_sr
      REAL*8 bblim,sslim,bbsum,yield,sm_trig(30),smtrig,xdaynpp,ans
      REAL*8 suma(360),tsuma,tsumam,stemfr,maint,yy,lmor_sc(3600)
      real sumrr,sumsr,sumlr,summr,resp_r,resp_s,resp_m,resp_l
      SAVE sumrr,sumsr,sumlr,summr
      INTEGER lai,leafls,stemls,rootls,mnth,day,ij,i,bb,gs,bbgs,sssum
      INTEGER bb2bbmin,bb2bbmax,bbm,ssm,sss,ss,dsbb
      INTEGER chill,dschill
      INCLUDE 'param.inc'

      bb2bbmin = 315
      bb2bbmax = 375
      gs = 30

      rem = rlai - int(rlai)
      lai = int(rlai) + 1
      resp = 0.0d0

      yield = 0.0d0

      ss = 0

*      print*
*      print*,'day npp ',mnth,day,daynpp*12
*      print*,nppstore*12,nppstorx*12,nppstor2*12
*      print*,'stemfr lai ',stemfr*12.0,rlai

      xdaynpp = daynpp

*----------------------------------------------------------------------*
* Compute active root and stem npp
*----------------------------------------------------------------------*
      stemnppr = 0.0d0
      rootnppr = 0.0d0
      DO ij=1,stemls
        stemnppr = stemnppr + stemv(ij)
      ENDDO
      DO ij=1,rootls
        rootnppr = rootnppr + rootv(ij)
      ENDDO

*----------------------------------------------------------------------*
* Pay for days roots if veg exists
*----------------------------------------------------------------------*
*      IF (tran/13.4/12.0d0*p_rootsc.LT.0.01d0*nppstore) THEN
*        yy = tran/13.4/12.0d0*p_rootsc
*      ELSE
*        yy = 0.01d0*nppstore
*      ENDIF
*      IF (respref.LT.0.00001d0) THEN
*        yy = 0.0d0
*      ENDIF

      IF ((nppstore.GT.0.0d0).AND.(daynpp.GT.0.0d0)) THEN
        yy = nppstore*p_rootfr
      ELSE
        yy = 0.0d0
      ENDIF

      nppstore = nppstore - yy
      xdaynpp = xdaynpp - yy

      rootnpp = rootnpp + yy
      s_rn = s_rn + yy

*----------------------------------------------------------------------*
* Age roots and set daily vector
*----------------------------------------------------------------------*
      DO ij=1,rootls-1
        rootv(rootls+1-ij) = rootv(rootls-ij)
      ENDDO
*      rootv(1) = (tran/13.4d0)/12.0d0*p_rootsc
      rootv(1) = yy
      IF (rootv(1).LT.0.0d0) THEN
        print*,'rootv negativ: rootv=',rootv(1),'tran=',tran
      ENDIF

*----------------------------------------------------------------------*
* Root resperation
*----------------------------------------------------------------------*
      resp_r = resp
      DO ij=1,rootls-1
        resp = resp + respref*rootv(rootls-ij)
        s_rr = s_rr + respref*rootv(rootls-ij)
        rootnpp = rootnpp - respref*rootv(rootls-ij)
        rootv(rootls-ij) = rootv(rootls-ij)*(1.0 - respref)
      ENDDO
      resp_r = resp - resp_r

*      print*,'root resp ',resp*12.0

*      nppstore = nppstore - rootnppr*respref
*      resp = resp + rootnppr*respref
*      xdaynpp = xdaynpp - rootnppr*respref
*      s_rr = s_rr + rootnppr*respref

*----------------------------------------------------------------------*
* Check for chilling.                                                  *
*----------------------------------------------------------------------*
      IF (chill.EQ.0) THEN
        bbsum = 0.0d0
        DO i=1,20
          IF (tmem(i).LT.-5.0d0)  bbsum = bbsum + max(-10.0d0,tmem(i)+
     &5.0d0)
        ENDDO
        IF (bbsum.LT.-100) THEN
          chill = 1
          dschill = 1
        ENDIF
      ENDIF
      IF (chill.EQ.1) THEN
        dschill = dschill + 1
      ENDIF
      IF (dschill.GT.260) THEN
        chill = 0
        dschill = 0
      ENDIF

*----------------------------------------------------------------------*
* Bubburst, if no budburst set, and sufficient soil moisture, then     *
* check for  budburst.                                                 *
*----------------------------------------------------------------------*
*      IF ((bb.eq.0).AND.(soil2g.GT.half+0.0*(wtfc-half)).AND.
*     &(day+mnth*30.eq.60)) THEN
*      print'(5i4,3f8.2)',mnth,day,bb,dsbb,bb2bbmax,smtrig,soil2g,half
      IF (((bb.eq.0).AND.(soil2g.GT.half+0.5d0*(wtfc-half))).OR.
     &((dsbb.GT.bb2bbmax).AND.(soil2g.GT.half+0.1d0*(wtfc-half)))) THEN
*      IF ((bb.eq.0).AND.(soil2g.GT.half+xxx*(wtfc-half)).AND.
*     &(dsw.GT.dswmin).AND.(dsw.LT.dswmax)) THEN

        smtrig = 0.0d0
        DO i=1,30
          smtrig = smtrig + sm_trig(i)
        ENDDO

        IF ((smtrig.GT.30.0d0).OR.(dsbb.GT.bb2bbmax)) THEN

*----------------------------------------------------------------------*
* Check for budburst using degree days.                                *
*----------------------------------------------------------------------*
          bbsum = 0.0d0
          DO i=1,bbm
            IF (tmem(i).GT.bb0)  bbsum = bbsum + min(bbmax,tmem(i)-bb0)
          ENDDO

          IF ((real(bbsum).GE.real(bblim)*exp(-0.01d0*real(dschill)))
     &.OR.(dsbb.GT.bb2bbmax)) THEN
*----------------------------------------------------------------------*
* Adjust proportion of gpp going into stem production based on suma.   *
* This is essentially the LAI control.                                 *
*----------------------------------------------------------------------*
            tsuma = 0.0d0
            DO i=1,360
              tsuma = tsuma + suma(i)
            ENDDO
!            print*,'tsuma ',tsuma
            maint = max(1.0d0,(real(leafls)/360.0d0)*1.0d0)
            tsuma = tsuma - leafmol*1.25d0/maint*p_opt
            IF (tsuma.GT.leafmol*1.25d0/maint*p_opt)
     &tsuma = leafmol*1.25d0/maint*p_opt
            IF (tsuma.LT.-leafmol*1.25d0/maint*p_opt)
     &tsuma = -leafmol*1.25d0/maint*p_opt
            tsumam = 0.0d0*tsumam + 1.0d0*tsuma
            stemfr = stemfr + tsumam*p_laimem

            IF (stemfr.LT.10.0d0) stemfr = 10.0d0

*----------------------------------------------------------------------*
* Bud burst occurance.                                                 *
*----------------------------------------------------------------------*
            bb = (mnth-1)*30 + day
            bbgs = 0
            dsbb = 0
            IF (stemfr.LT.0.75d0*nppstore) THEN
              nppstor2 = nppstore - stemfr
            ELSE
              IF (stemfr.LT.0.75d0*nppstore) THEN
                nppstor2 = nppstore - stemfr
              ELSE
                nppstor2 = nppstore*0.25d0
              ENDIF
              stemfr = stemfr*0.95
            ENDIF
            laiinc = (nppstore - nppstor2)/leafmol/1.25d0
            nppstorx = nppstore
!            print*,'laiinc ',laiinc,tsuma,stemfr*12.0
          ENDIF
        ENDIF
      ENDIF

*----------------------------------------------------------------------*
* Compute length of growing season, and set to zero when equal to      *
* growing season.                                                      *
*----------------------------------------------------------------------*
      IF (bb.GT.0)  bbgs = bbgs + 1
      IF (bbgs-gs.gt.bb2bbmin) THEN
        bb = 0
        bbgs = 0
      ENDIF

      IF (dsbb.LT.500) dsbb = dsbb + 1

*----------------------------------------------------------------------*
* Set LAI increase.                                                    *
*----------------------------------------------------------------------*
      IF ((bb.GT.0).AND.(bbgs.LT.gs).AND.(nppstore.GT.1.0d0)) THEN
        laiinc = lairat*(nppstorx - nppstor2)/leafmol/1.25d0
*        IF (rlai+laiinc.LT.0.5d0)  laiinc = 0.5d0 - rlai
        IF (rlai+laiinc.GT.maxlai)  laiinc = maxlai - rlai
        IF (rlai+laiinc.GT.11.5d0)  laiinc = 11.5d0 - rlai
        IF ((rlai.GT.0).and.(nppstore.LT.5.0)) laiinc = 0.0d0
      ELSE
        laiinc = 0.0d0
      ENDIF

*----------------------------------------------------------------------*
* Senescence, if rlai is greater than zero, compute senescence.        *
*----------------------------------------------------------------------*
      IF (rlai.GT.1.0e-6) THEN
        IF (soil2g.LT.half*0.0d0) THEN
*----------------------------------------------------------------------*
* Senescence event due to soil moisture.                               *
*----------------------------------------------------------------------*
          laiinc = -rlai
          ss = day + (mnth - 1)*30
*          print*,'senescence soil'
        ELSEIF (bbgs.GT.100) THEN
*----------------------------------------------------------------------*
* Check for senescence, senescence occurs there are 'sss' days colder  *
* than 'sslim' out of the last 'ssm' days.                             *
*----------------------------------------------------------------------*
          sssum = 0
          DO i=1,ssm
            IF (tmem(i).LT.sslim)  sssum = sssum + 1
          ENDDO
          IF (sssum.GE.sss) THEN
*----------------------------------------------------------------------*
* Senescence event due to temperature.                                 *
*----------------------------------------------------------------------*
            laiinc =-rlai
            ss = day + (mnth - 1)*30
*          print*,'senescence temp'
*          print'(100f6.1)',(tmem(i),i=1,ssm)
          ENDIF
        ENDIF
      ENDIF

*----------------------------------------------------------------------*
* Pay for new leaves.
*----------------------------------------------------------------------*
      resp_m = resp - resp_m
      IF (laiinc.GT.0.0d0) THEN
        nppstore = nppstore - laiinc*leafmol*1.25d0
        nppstorx = nppstorx - laiinc*leafmol*1.25d0
        leafnpp = leafnpp + laiinc*leafmol
        resp = resp + 0.25d0*laiinc*leafmol
        xdaynpp = xdaynpp - laiinc*leafmol*1.25d0

        s_ln = s_ln + laiinc*leafmol*1.25d0
      ENDIF
*      print*,'maint & lai inc ',12.0*0.25d0*laiinc*leafmol,laiinc
      s_mr = 0.25d0*laiinc*leafmol
      resp_m = resp - resp_m

*      print*,'plus leaf resp ',resp*12.0

*----------------------------------------------------------------------*
* Age leaves by one day, kill any which have died of old age,
* adjust by laiinc (+ or -), and then sum to get 'rlai'.
* 'leafls' is an integer variable of leaf lifespan in days.
*----------------------------------------------------------------------*
      CALL LAIALT(leafv,rlai,laiinc,leafls,leaflit)

      rem = rlai - int(rlai)
      lai = int(rlai) + 1

*----------------------------------------------------------------------*
* leaf death not through age mortality
*----------------------------------------------------------------------*
      rlai = rlai + leaflit
      DO i=1,leafls
        IF (leafv(i).GT.0.0d0) THEN
          IF (leafv(i).GT.1.0e-6) THEN
            leaflit = leaflit + leafv(i)
            leafv(i)=leafv(i)*lmor_sc(i)
            leaflit = leaflit - leafv(i)
          ELSE
            leaflit = leaflit + leafv(i)
            leafv(i) = 0.0d0
          ENDIF
        ENDIF
      ENDDO
      rlai = rlai - leaflit

*----------------------------------------------------------------------*
* Stem resperation, and NPP                                            *
*----------------------------------------------------------------------*
*      nppstore = nppstore - stemnppr*respref
*      resp = resp + stemnppr*respref
*      s_sr = s_sr + stemnppr*respref
*      xdaynpp = xdaynpp - stemnppr*respref

      resp_s = resp - resp_s
      s_sr = resp
      DO ij=1,stemls-1
        resp = resp + respref*stemv(stemls-ij)
        stemnpp = stemnpp - respref*stemv(stemls-ij)
        stemv(stemls-ij) = stemv(stemls-ij)*(1.0 - respref)
      ENDDO
      s_sr = resp - s_sr
      resp_s = resp - resp_s

      DO ij=1,stemls-1
        stemv(stemls+1-ij) = stemv(stemls-ij)
      ENDDO

      IF ((nppstore.GT.0.0d0).and.(daynpp.GT.0.0d0)) THEN
        yy = nppstore*p_stemfr
        stemv(1) = yy
        s_sn = s_sn + yy
        stemnpp = stemnpp + yy
        nppstore = nppstore - yy
      ELSE
        stemv(1) = 0.0d0
      ENDIF

      IF (stemnppr.GT.0d0) THEN
        IF ((stemv(1).LT.p_stmin)) then
          stemnpp = stemnpp + p_stmin - stemv(1)
          nppstore = nppstore - p_stmin + stemv(1)
          s_sn = s_sn + p_stmin - stemv(1)
          stemv(1) = p_stmin
        ENDIF
      ENDIF

*      print*,mnth,day,stemv(1),stemnppr*respref

*----------------------------------------------------------------------*

      IF (rlai.GT.laimax)  laimax = rlai
      IF (.NOT.(rlai.GT.0.0d0))  daysoff = daysoff + 1.0d0

      IF (nppstore.LT.0.0d0) THEN
        DO i=1,leafls
          leaflit = leaflit + leafv(i)
          leafv(i) = 0.0d0
        ENDDO
c added by Ghislain 07/10/03
        stemlit = 0.0d0
        DO i=1,stemls
          stemlit = stemlit + stemv(i)
          stemv(i) = 0.0d0
        ENDDO
c added by Ghislain 07/10/03
        rootlit = 0.0d0
        DO i=1,rootls
          rootlit = rootlit + rootv(i)
          rootv(i) = 0.0d0
        ENDDO

        nppstore = 0.0d0
        rlai = 0.0d0

      ELSE
        stemlit = 0.0d0
        rootlit = 0.0d0
      ENDIF

      ans = 0.0d0
      DO i=1,stemls
        ans = ans + stemv(i)
      ENDDO

!      print*,'plus stem resp ',resp*12.0
!      print*,nppstore*12,nppstorx*12,nppstor2*12
!      print*,'stemfr lai ',12.0*stemfr,rlai

!      print*,mnth,day,s_sr*12.0,s_rr*12.0

      if ((day == 1).and.(mnth == 1)) then
        sumsr=0.0
        sumrr=0.0
        sumlr=0.0
        summr=0.0
      endif
      sumsr = sumsr + 12.0*s_sr
      sumrr = sumrr + 12.0*s_rr
      sumlr = sumlr + 12.0*s_lr
      summr = summr + 12.0*s_mr
      if ((day == -30).and.(mnth == 12))
     &print'(4f10.3)',sumsr,sumrr,sumlr,summr

      RETURN
      END

*----------------------------------------------------------------------*
*                          LAIALT                                      *
*----------------------------------------------------------------------*
      SUBROUTINE LAIALT(leafv,rlai,laiinc,leafls,leaflit)
*----------------------------------------------------------------------*
      IMPLICIT NONE
      REAL*8 leafv(3600),rlai,laiinc,sum,leaflit
      INTEGER leafls,i

      leaflit = leafv(leafls)
      DO i=1,leafls-1
        leafv(leafls+1-i) = leafv(leafls-i)
      ENDDO
      leafv(1) = 0.0d0

      IF (laiinc.GT.0.0d0) THEN
        leafv(1) = laiinc
      ELSEIF (laiinc+leaflit.LT.0.0d0) THEN
        sum = laiinc+leaflit
        i = leafls
20      CONTINUE
          sum = sum + leafv(i)
          IF (sum.GT.1.0e-6) then
            leaflit = leaflit + leafv(i) - sum
            leafv(i) = sum
          ELSE
            leaflit = leaflit + leafv(i)
            leafv(i) = 0.0d0
            i = i - 1
          ENDIF
        IF ((.NOT.(sum.GT.0.0d0)).AND.(i.GT.0))  GOTO 20
      ENDIF

      rlai = 0.0d0
      DO i=1,leafls
        rlai = rlai +  leafv(i)
      ENDDO

      RETURN
      END


text:
	SUBROUTINE GROAM  ( ityp, sys, x, y, iret )
C************************************************************************
C* GROAM								*
C* 									*
C* This subroutine moves the current roam window to the specified       *
C* position in any coordinate system except 'S'.  The base point of the *
C* roam can be the upper left of the screen or the center of the screen.*
C*									*
C* GROAM  ( ITYP, SYS, X, Y, IRET )					*
C*									*
C* Input parameters:							*
C* 	ITYP		INTEGER		The base point of roam 		*
C*                                        0 = upper left screen corner  *
C*                                        1 = center of the screen 	*
C*                                        2 = delta_x, delta_y		*
C*      SYS             CHAR*           Coordinate system               *
C*                                        'S' = screen coordinates      *
C*                                        'D' = device coordinates      *
C*                                        'N' = normalized coordinates  *
C*                                        'V' = view coordinates        *
C*                                        'P' = plot coordinates        *
C*                                        'M' = map coordinates         *
C*                                        'G' = grid coordinates        *
C* 	X		REAL		Upper left x coordinate 	*
C* 	Y		REAL		Upper left y coordinate 	*
C*									*
C* Output parameters:							*
C* 	IRET		INTEGER		Return code			*
C**									*
C* Log:									*
C* C. Lin/EAI	 	6/97						*
C* A. Hardy/GSC         6/98		Cleaned up prolog               *
C************************************************************************
	INCLUDE		'FUNCCODE.PRM'
	INCLUDE		'ERROR.PRM'
	INCLUDE		'GEMPRM.PRM'
C*
        CHARACTER*(*)   sys
C*
	INTEGER		isend (4)
	REAL		rsend (2)
C------------------------------------------------------------------------
	iret = NORMAL
C
C*      Check validity of the coordinate system.
C
        isys = INDEX ( syslo, sys ) + INDEX ( sysup, sys )
        IF  ( isys .eq. 0 ) THEN
            iret = NOCORD
            RETURN
        END IF
C
C*	Load input parameters into buffer and write them to the mailbox.
C
	isend (1) = 6
	isend (2) = FROAM
	isend (3) = ityp
	isend (4) = isys
C
	CALL GPUT  ( isend, 4, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
	rsend (1) = x
	rsend (2) = y
C
	CALL GPUTR  ( rsend, 2, iret )
	IF  ( iret .ne. NORMAL )  RETURN
C
        CALL GGET  ( iret, 1, ierr )
        IF  ( ierr .ne. NORMAL )  iret = ierr
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


module procmod

integer :: counter
integer :: counter_mutex

contains

    subroutine example(arg)
    implicit none
    integer :: arg
    integer :: info

    integer :: i
    integer :: lcnt

    do i=1,100
        call forthread_mutex_lock(counter_mutex,info)
        lcnt = counter
        lcnt = lcnt + 1
        counter = lcnt
        call forthread_mutex_unlock(counter_mutex,info)
    end do
    end subroutine example
end module procmod

program test04

use procmod

!include '../src/iface/ciface.h'

integer, parameter       :: n = 100
integer, dimension(n)    :: tid
integer                  :: arg
integer                  :: ret_val
integer                  :: thread_id
integer                  :: aid
integer                  :: i
integer                  :: info

aid = -1
arg = 0
print *, '====================== test04 ======================'
print *,'initializing threads'
call forthread_init(info)
if (info.ne.0) then
    print *,'error initializing '
endif

call forthread_mutex_init(counter_mutex,aid,info)

do i=1,n
  call forthread_create(thread_id,aid,example,arg,info)
  tid(i) = thread_id
  if (info.ne.0) then
    print *,'error creating thread ',i
  endif
enddo

do i=1,n
  call forthread_join(tid(i),ret_val,info)
  if (info.ne.0) then
    print *,'error joining thread ',i
  endif
enddo

call forthread_mutex_destroy(counter_mutex,info)
call forthread_destroy(info)
if (info.ne.0) then
    print *,'error destroying.'
endif

print *,'counter (10000):',counter

print *, '==================== end test04 ===================='

end program test04


text:
subroutine random(idum,prandom)
   implicit none
   INTEGER,intent(inout) :: idum
   double precision,intent(out) :: prandom
   integer,parameter :: IA=(16807),IM=(2147483647),IQ=(127773),  &
      IR=(2836),MASK=(123459876)
   double precision :: AM
   INTEGER :: k
   AM=1./IM
   idum=ieor(idum,MASK)
   k=idum/IQ
   idum=IA*(idum-k*IQ)-IR*k
   if (idum.lt.0) idum=idum+IM
   prandom=AM*idum
   idum=ieor(idum,MASK)
end subroutine random

subroutine init_rand(r)
   implicit none
   integer,intent(out) :: r
   double precision :: r1

   call random(r,r1)

end subroutine init_rand

subroutine gaussian(r,grandom)
   implicit none
   integer,intent(inout) :: r
   double precision,intent(out) :: grandom

   double precision :: u1, v1, u2, v2, s

   do
      call random(r,u1)
      call random(r,u2)
      v1 = 2 * u1 - 1.d0
      v2 = 2 * u2 - 1.d0
      s = v1**2 + v2**2
      if ((s.gt.0.d0).and.(s .lt. 1.d0))exit
   enddo

   grandom = dsqrt(-2.d0*dlog(s)/s)*v1

end subroutine gaussian

subroutine sde_euler(n,x,f_func,g_func,r,dt)
   implicit none
   integer, intent(in) :: n
   double precision, intent(in) :: dt
   double precision,intent(inout) :: x(n)
   integer,intent(inout) :: r
   external :: f_func, g_func

   double precision :: f(n), g(n), v(n), dt2
   integer :: i

   dt2 = dsqrt(dt)
   call f_func(n,x,f)
   call g_func(n,x,g)
   do i = 1, n
      call gaussian(r,v(i))
   enddo

   x = x + f * dt + g*v * dt2

end subroutine sde_euler

subroutine sde_milstein(n,x,f_func,g_func,dg_func,r,dt)
   implicit none
   integer, intent(in) :: n
   double precision, intent(in) :: dt
   double precision,intent(inout) :: x(n)
   integer,intent(inout) :: r
   external :: f_func, g_func, dg_func

   double precision :: f(n), g(n), dg(n), v(n), dt2
   integer :: i

   dt2 = dsqrt(dt)
   call f_func(n,x,f)
   call g_func(n,x,g)
   call dg_func(n,x,dg)
   do i = 1, n
      call gaussian(r,v(i))
   enddo

   x = x + f * dt + g*v * dt2 + &
      0.5d0* g * dg * (v*v-1.d0) * dt

end subroutine sde_milstein

subroutine sde_rk(n,x,f_func,g_func,r,dt)
   implicit none
   integer, intent(in) :: n
   double precision, intent(in) :: dt
   double precision,intent(inout) :: x(n)
   integer,intent(inout) :: r
   external :: f_func, g_func

   double precision :: f(n), g(n), v(n), dt2, x1(n), g1(n), x2(n)
   integer :: i

   dt2 = dsqrt(dt)
   call f_func(n,x,f)
   call g_func(n,x,g)
   do i = 1, n
      call gaussian(r,v(i))
   enddo

   x2 = x + f * dt !+ g*v * dt2
   x1 = x2 + g * dt2

   call g_func(n,x1,g1)

   x = x2 + g*v*dt2 + 0.5d0 *(g1 - g) * (v*v-1.d0) * dt2

end subroutine sde_rk

function hill_inh(x, x0, n)
   implicit none
   integer, intent(in) :: n
   double precision,intent(in) :: x, x0
   double precision :: hill_inh

   double precision :: a

   a = (x/x0)**n
   hill_inh = 1/(1+a)

end function hill_inh

subroutine f_ts(n, x, f)
   implicit none
   integer, intent(in) :: n
   double precision,intent(in) :: x(n)
   double precision,intent(out) :: f(n)
   double precision :: g0 = (10.d0), g1 = (40.d0), x_hill = (100.d0), k = (0.1d0)
   integer :: n_hill = (4)
   double precision :: hill_inh

   f(1) =  g0 + g1*hill_inh(x(2), x_hill, n_hill) - k * x(1)
   f(2) =  g0 + g1*hill_inh(x(1), x_hill, n_hill) - k * x(2)

end subroutine f_ts

subroutine g_ts(n, x, g)
   implicit none
   integer, intent(in) :: n
   double precision,intent(in) :: x(n)
   double precision,intent(out) :: g(n)

   g(1:2) = 20.d0

end subroutine g_ts

subroutine sde_simulation(iseed, t_total, kappa, nrun)
   implicit none
   integer,intent(in) :: iseed
   double precision,intent(in) :: t_total
   double precision,intent(out) :: kappa(2)
   integer, intent(out) :: nrun(2)

   double precision :: t, r1, r2, dt, t0, trun, x(2)
   integer :: n, state, r, t_relax
   double precision :: mfpt(2)
   logical :: start
   external :: f_ts, g_ts

   n = 2
   dt = 0.01d0
   t_relax = 100.d0
   r = iseed
   call init_rand(r)  ! initiate random number generator

   call random(r, r1)
   call random(r, r2)
   x(:) = (/r1*300.d0, r2*300.d0/)   ! the initial condition

   t = 0
   t0 = 0
   start = .false.
   nrun(:) = 0
   mfpt(:) = 0.d0
   do
      call sde_euler(n,x,f_ts,g_ts,r,dt)
      t = t + dt
      if(t > t_total)exit
      if(.not.start) then
        if(t - t0 > t_relax)then
          if(x(1) < x(2))then
            state = 1
          else
            state = 2
          endif
          trun = 0.d0
          start = .true.
        endif
      else
        if((state.eq.1).and.(X(1) > X(2)))then
          nrun(1) = nrun(1) + 1
          mfpt(1) = mfpt(1) + trun
          t0 = t
          start = .false.
          trun = 0.d0
          state = 2
        elseif((state.eq.2).and.(X(1) <= X(2)))then
          nrun(2) = nrun(2) + 1
          mfpt(2) = mfpt(2) + trun
          t0 = t
          start = .false.
          trun = 0.d0
          state = 1
        endif
      endif

      trun = trun + dt
   enddo

   mfpt(:) = mfpt(:)/nrun(:)
   kappa(:) = 1.d0/2.d0/mfpt(:)

end subroutine sde_simulation

max_stars_repo_name: lusystemsbio/numericalR

max_issues_repo_name: lusystemsbio/numericalR

max_forks_repo_name: lusystemsbio/numericalR