program hmipt_matsuara_dmft_g0
  USE DMFT_IPT
  USE IOTOOLS
  USE ERROR
  implicit none
  logical                :: converged,check
  real(8)                :: n,z
  integer                :: i,iloop
  complex(8)             :: zeta
  complex(8),allocatable :: fg(:),fg0(:),sigma(:),GFold(:)
  real(8),allocatable    :: wm(:),sigt(:)
  real(8) :: C0,C1,n0

  call read_input("inputIPT.in")

  !allocate functions:
  allocate(fg(L))
  allocate(sigma(L))
  allocate(fg0(L))
  allocate(sigt(0:L))

  !build freq. array
  allocate(wm(L))
  wm(:)  = pi/beta*real(2*arange(1,L)-1,8)

  !set the energy scale:
  D=2.d0*ts

  !get or read first sigma
  call  get_initial_weiss_field(fg0,"g0.restart")
  call splot("guessG0_iw.ipt",wm,fg0)

  !dmft loop:
  iloop=0 ; converged=.false.
  do while(.not.converged.AND.iloop<nloop)
     iloop=iloop+1
     write(*,"(A,i5)",advance="no")"DMFT-loop",iloop
     !IMPURITY SOLVER
     sigma= solve_ipt_matsubara(fg0)
     !
     sigma=xi*dimag(sigma)!ACTHUNG
     !
     !SELF-CONSISTENCY:
     do i=1,L
        zeta  = xi*wm(i) - sigma(i)
        fg(i) = gfbethe(wm(i),zeta,D)
     enddo
     n   = get_local_density(fg,beta)
     !update weiss-field:
     fg0 = one/(one/fg + sigma)
     !if(iloop>1)fg0 = weight*fg0 + (1.d0-weight)*GFold
     !
     converged=check_convergence(fg0,dmft_error,nsuccess,nloop)
     !GET OBSERVABLES
     z=1.d0 - dimag(sigma(1))/wm(1);z=1.d0/z
     call splot("observables_all.ipt",iloop,beta,xmu,u,z,n,append=.true.)
  enddo
  call splot("G_iw.ipt",wm,fg)
  call splot("G0_iw.ipt",wm,fg0)
  call splot("Sigma_iw.ipt",wm,sigma)
  call splot("observables.ipt",iloop,u,beta,xmu,z,n)

  n0=n/2.d0
  C0=U*(n0-0.5d0)
  C1=U**2*n0*(1.d0-n0)
  print*,n0,C0,C1
  sigma=sigma- C0 - C1/(xi*wm)
  call fftgf_iw2tau(sigma,sigt(0:),beta,notail=.true.)
  sigt=sigt-C1*0.5d0
  open(100,file="fft_sigma_iw.ipt")
  do i=0,L
     write(100,*)i*beta/dble(L),sigt(i)
  enddo
  close(100)

contains

  subroutine get_initial_weiss_field(self,file)
    complex(8),dimension(:)       :: self
    real(8),dimension(size(self)) :: wm_
    character(len=*)              :: file
    logical                       :: check
    inquire(file=file,exist=check)
    if(check)then
       print*,'Reading G0'
       call sread(file,wm_,self)
    else
       print*,"Using non-interacting GF as guess:"
       call bethe_guess_g0(self,D,beta,0.d0)
    endif
  end subroutine get_initial_weiss_field

end program hmipt_matsuara_dmft_g0


PROGRAM test_cycle
INTEGER :: i
DO i = 1, 5
   IF ( i == 3 ) CYCLE
   WRITE (*,*) i
END DO
WRITE (*,*) 'End of loop!'
END PROGRAM test_cycle


text:
! ###########################################################
! ###########################################################
module CN_Thermal
    use CN_BaseClass
    implicit none
!!!        ---------------------------------------
    type,extends( CN_Base):: ThermalClass
        real*8 , ALLOCATABLE :: rThermal_DijSij(:)

        contains
!!!!!!!!!        procedure ::ThermalClass
        procedure ::CN_Thermal_CleanMem
!!!!!!!!!            Interface CN_Thermal_CleanMem;                Module Procedure CN_Thermal_CleanMem ;          End Interface
    end type ThermalClass
!            Interface init_ThermalClass;                       Module Procedure init_ThermalClass ;            End Interface
!!!        ---------------------------------------
!!!        ---------------------------------------
        contains
! ###########################################################
! ###########################################################
        subroutine CN_Thermal_CleanMem (tCN_object)
            class ( ThermalClass ), intent(inout) :: tCN_object

            call tCN_object%CleanMem()
            IF (ALLOCATED (tCN_object%rThermal_DijSij))           DEALLOCATE (tCN_object%rThermal_DijSij)

        end subroutine CN_Thermal_CleanMem
end module CN_Thermal
! ###########################################################
! ###########################################################
type ( ThermalClass )function init_ThermalClass ( mNodCountB,mElemCountB )
    use CN_Thermal
    integer, intent(in) :: mNodCountB,mElemCountB

        init_ThermalClass%CN_Base=init_CN_Base(mNodCountB,mElemCountB)
        init_ThermalClass%iSolutionActive=1
        Call init_ThermalClass%AllocateMem()

        allocate(init_ThermalClass%rThermal_DijSij(mElemCountB))

end function init_ThermalClass

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


text:
      MODULE ENV_VARS

C*************************************************************************
C
C  FUNCTION:  To get environment variables
C
C  PRECONDITIONS: None
C
C  KEY SUBROUTINES/FUNCTIONS CALLED: None
C
C  REVISION HISTORY: Prototype created by Jerry Gipson, July, 1999
C  Modified by Chris Nolte December 2008 to add ALLVARS flag
C
C*************************************************************************

      INTEGER       LOGUNIT         ! Output log unit number

      INTEGER       LOCOL           ! Lowest column of block
      INTEGER       LOROW           ! Lowest row of block
      INTEGER       LOLEV           ! Lowest level of block

      INTEGER       HICOL           ! Highest column of block
      INTEGER       HIROW           ! Highest row of block
      INTEGER       HILEV           ! Highest level of block

      INTEGER       NVAR            ! Number of variables to output

      INTEGER       SDATE           ! Start date
      INTEGER       STIME           ! Start time

      INTEGER       NSTEPS          ! No. of steps to output

c      CHARACTER*16  SUMORAVG        ! Output sums or avvgs ( default=avg)

      CHARACTER*16  OUTFORMAT       ! Output time format

      CHARACTER*16  TIME_CONV       ! Time zone to convert to

      CHARACTER*256 OUT_FNAME       ! Name of output data file

      CHARACTER*16, ALLOCATABLE :: VNAME( : ) ! Variable names

      LOGICAL       ALLVARS         ! flag if output all vars--specify 'ALL' in specfile

      CONTAINS

         SUBROUTINE GET_ENVS
C*************************************************************************
C
C  FUNCTION:  To get environment variables
C
C  PRECONDITIONS: None
C
C  KEY SUBROUTINES/FUNCTIONS CALLED: None
C
C  REVISION HISTORY: Prototype created by Jerry Gipson, July, 1999
C
C*************************************************************************
         USE M3UTILIO
         IMPLICIT NONE

C..ARGUMENTS: None

C..PARAMETERS:
         INTEGER    N_VALTZ
         PARAMETER( N_VALTZ = 2 )

C..SAVED LOCAL VARIABLES: None

C..SCRATCH LOCAL VARIABLES:
         CHARACTER*16    ENV_DFLT         ! Environment variable default value
         CHARACTER*16    OUTFILE          ! Env var for output file name
         CHARACTER*16    PNAME            ! Program Name
         CHARACTER*80    ENV_DESC         ! Environment variable description
         CHARACTER*80    MSG              ! Error message
         CHARACTER*256   RET_VAL          ! Returned value of env var
         CHARACTER*9     FORMAT

         CHARACTER*3     VALTZ( N_VALTZ )      ! Valid time zone codes

         INTEGER   N                      ! Loop index
         INTEGER   STATUS                 ! Status code

         LOGICAL   LERROR                 ! Error flag

C**********************************************************************
         DATA  PNAME           / 'GET_ENVS'        /
         DATA  OUTFILE         / 'OUTFILE'         /
         DATA  FORMAT          / 'OUTFORMAT' /
         DATA VALTZ / 'GMT', 'EST' /

         LOGUNIT = INIT3()

         LERROR = .FALSE.

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the variables to retrieve
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         ALLOCATE ( VNAME( MXVARS3 ) )

         OPEN( UNIT = 9, FILE='specfile' )

         DO N = 1, MXVARS3
            READ( 9, '( A )', END = 100 ) VNAME( N )
            NVAR = NVAR + 1
         ENDDO

  100    CONTINUE

         CLOSE( 9 )

         ALLVARS = ( VNAME( 1 ) .EQ. 'ALL' )

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the start date and time
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         SDATE = 0
         ENV_DESC = 'Start date for output'
         SDATE = ENVINT( 'SDATE', ENV_DESC, SDATE, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: SDATE not specified; using file(s) start date'
            CALL M3MESG( MSG )
         ENDIF

         STIME = 0
         ENV_DESC = 'Start time for output'
         STIME = ENVINT( 'STIME', ENV_DESC, STIME, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: STIME not specified; using 000000'
            CALL M3MESG( MSG )
         ENDIF

         NSTEPS = 0
         ENV_DESC = 'Start time for output'
         NSTEPS = ENVINT( 'NSTEPS', ENV_DESC, NSTEPS, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: NSTEPS not specified; using file end time'
            CALL M3MESG( MSG )
         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the block definition
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         LOCOL = 0
         ENV_DESC = 'Lowest column number'
         LOCOL = ENVINT( 'LOCOL', ENV_DESC, LOCOL, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: LOCOL not specified; will be set to 1'
            CALL M3MESG( MSG )
         ENDIF

         HICOL = 0
         ENV_DESC = 'Highest column number'
         HICOL = ENVINT( 'HICOL', ENV_DESC, HICOL, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: HICOL not specified; will be set to NCOLS'
            CALL M3MESG( MSG )
         ENDIF

         LOROW = 0
         ENV_DESC = 'Lowest row number'
         LOROW = ENVINT( 'LOROW', ENV_DESC, LOROW, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: LOROW not specified; will be set to 1'
            CALL M3MESG( MSG )
         ENDIF

         HIROW = 0
         ENV_DESC = 'Highest row number'
         HIROW = ENVINT( 'HIROW', ENV_DESC, HIROW, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: HIROW not specified; will be set to NROWS'
            CALL M3MESG( MSG )
         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the levels to use
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         LOLEV = 0
         ENV_DESC = 'Lowest level number'
         LOLEV = ENVINT( 'LOLEV', ENV_DESC, LOLEV, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: LOLEV not specified; wuill be set to 1'
            CALL M3MESG( MSG )
         ENDIF

         HILEV = 0
         ENV_DESC = 'Highest level'
         HILEV = ENVINT( 'HILEV', ENV_DESC, HILEV, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: HILEV not specified; will be set to NLAYS'
            CALL M3MESG( MSG )
         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the name of the output file
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         ENV_DFLT = 'OUTFILE'
         ENV_DESC = 'Output file name'
         CALL ENVSTR( OUTFILE, ENV_DESC, ENV_DFLT, OUT_FNAME, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'ERROR: no file name assigned to OUTFILE'
            CALL M3MESG( MSG )
            LERROR = .TRUE.
         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the time convention to use on output
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         TIME_CONV = 'GMT'
         ENV_DFLT  = 'GMT'
         ENV_DESC  = 'Time convention to use for output data'
         CALL ENVSTR( 'TIME_ZONE', ENV_DESC, ENV_DFLT, TIME_CONV, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: No time zone assigned; using default GMT'
            CALL M3MESG( MSG )
         ENDIF

         IF( INDEX1( TIME_CONV, N_VALTZ, VALTZ ) .EQ. 0 ) THEN
            MSG = 'ERROR: Invalid time zone for output file'
            CALL M3MESG( MSG )
            LERROR = .TRUE.
         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get the format to use for output
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
         ENV_DFLT  = 'IOAPI'
         ENV_DESC  = 'Time format to use for output data'

         CALL ENVSTR( FORMAT, ENV_DESC, ENV_DFLT, RET_VAL, STATUS)

         IF( STATUS .NE. 0 ) THEN
            MSG = 'WARNING: No output time format assigned; using ' //
     &            ' IOAPI convention'
            CALL M3MESG( MSG )
            OUTFORMAT = 'IOAPI'
         ELSE
            OUTFORMAT = RET_VAL
         ENDIF

         IF( OUTFORMAT .NE. 'SAS' .AND. OUTFORMAT .NE. 'IOAPI') THEN
            MSG = 'WARNING: Invalid output time format requested: ' //
     &            'using IOAPI convention'
            CALL M3MESG( MSG )
            OUTFORMAT = 'IOAPI'
         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Get sum or avg flag
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c         ENV_DFLT  = 'AVG'
c         ENV_DESC  = 'Flag to output sum or average'
c         CALL ENVSTR( 'SUMORAVG', ENV_DESC, ENV_DFLT, RET_VAL, STATUS)

c         IF( STATUS .NE. 0 ) THEN
c            MSG = 'WARNING: Sums or averages not specified; averages ' //
c     &            ' output'
c            CALL M3MESG( MSG )
c            SUMORAVG = 'AVG'
c         ELSE
c            SUMORAVG = RET_VAL
c         ENDIF

c         IF( SUMORAVG .NE. 'AVG' .AND. SUMORAVG .NE. 'SUM') THEN
c            MSG = 'WARNING: Invalid SUMORAVG value: using AVG'
c            CALL M3MESG( MSG )
c            SUMORAVG = 'AVG'
c         ENDIF

cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc
c  Stop if errors detected
cccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccccc

         IF( LERROR ) THEN
            MSG = 'Error in assigning environment variables -- stopping'
            CALL M3ERR( PNAME, 0, 0, MSG, .TRUE. )
         ENDIF

         RETURN

         END SUBROUTINE GET_ENVS

      END MODULE ENV_VARS

ext: f

max_stars_repo_name: Simeng-unique/CMAQ-changed

max_issues_repo_name: Simeng-unique/CMAQ-changed

max_forks_repo_name: Simeng-unique/CMAQ-changed


text:
Rising Star Stables Inc. is an equestrian facility offering boarding and breeding of Arabians, straight Egyptian and Egyptianbred horses. Facilities include both an arena and pastureland; boarders are turned out daily and fed twice a day.

To learn more about equestrian resources in or near Davis, please visit our Horses page.

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
	SUBROUTINE TG_RANG  ( gdattm, ntimf, filtim, ntime, times,
     +			      iret )
C************************************************************************
C* TG_RANG								*
C*									*
C* This subroutine converts the user input for a grid time range	*
C* into a list of times.						*
C*									*
C* TG_RANG  ( GDATTM, NTIMF, FILTIM, NTIME, TIMES, IRET )		*
C*									*
C* Input parameters:							*
C*	GDATTM		CHAR*		Input grid time			*
C*	NTIMF		INTEGER		Number of times in file		*
C*	FILTIM (NTIMF)	CHAR*		Times in file			*
C*									*
C* Output parameters:							*
C*	NTIME		INTEGER		Number of times selected	*
C*	TIMES  (NTIME)	CHAR*		Selected times			*
C*	IRET		INTEGER		Return code			*
C*					  0 = normal return		*
C*					 -4 = no times in file		*
C*					 -5 = invalid time range	*
C*					 -6 = single time invalid	*
C*					 -7 = invalid forecast type	*
C*					 -8 = invalid first time	*
C*					 -9 = invalid last time		*
C*					-10 = first time after last	*
C*					-11 = first and last are same	*
C*					-12 = forecast types not same	*
C*					-13 = forecast times not same	*
C*					-14 = invalid time increment	*
C*					-15 = too many times in list	*
C*					-16 = no times in range		*
C**									*
C* Log:									*
C* M. desJardins/GSFC	 4/89						*
C* M. desJardins/GSFC	 7/89	Added iret = -16 for no times in range	*
C* G. Huffman/USRA	10/89	Make TG_RNOI, TG_RINC return IRET	*
C* K. Brill/NMC		11/91	Remove string index range from sub call *
C* D. Kidwell/NCEP	 3/99	Fixed for Y2K                 		*
C* S. Jacobs/NCEP	12/00	Removed checks for begin time = end time*
C************************************************************************
	CHARACTER*(*)	gdattm, filtim (*), times (*)
C*
	CHARACTER	firstm*20, lasttm*20, tstart*20, tstop*20,
     +			tinc*20, tbegin*20, tend*20, ctype*1, ginput*48,
     +			ccc*48, cbegin*2, cend*2
C------------------------------------------------------------------------
	iret   = 0
	ntime  = 0
	CALL ST_LCUC  ( gdattm, ginput, ier )
C
C*	Make sure there are times in the file and that an output time
C*	has been requested.
C
	IF  ( ntimf .lt. 1 )  THEN
	    iret = -4
	    RETURN
	  ELSE IF  ( ginput .eq. ' ' )  THEN
	    iret = -5
	    RETURN
	END IF
C
C*	Get first and last time in file.
C
	firstm = filtim (1)
	lasttm = filtim (ntimf)
C
C*	Break range into parts.
C
	CALL ST_RANG  ( ginput, tstart, tstop, tinc, itype, ier )
C
C*	If this is not a range, check for "ALL".
C
	IF  ( itype .eq. 0 )  THEN
	    indall = INDEX ( ginput, 'ALL' )
C
C*	    If ALL is not included, this is not a valid time range.
C
	    IF  ( indall .eq. 0 )  THEN
		iret = -6
		RETURN
C
C*		If the entire string is "ALL", set first and last times.
C
	      ELSE IF  ( ginput .eq. 'ALL' )  THEN
		tstart = firstm
		tstop  = lasttm
		itype  = 1
C
C*		If ALL is first, but not entire string, add first and
C*		last times.
C
	      ELSE IF  ( indall .eq. 1 )  THEN
		tstart = firstm ( :11 ) // ginput ( 4: )
		tstop  = lasttm ( :11 ) // ginput ( 4: )
		itype  = 1
C
C*		Check for FALL, AALL, GALL.
C
	      ELSE
		ctype = ginput ( indall-1 : indall-1 )
		IF  ( ( ctype .eq. 'F' ) .or. ( ctype .eq. 'A' ) .or.
     +		      ( ctype .eq. 'G' ) .or. ( ctype .eq. 'I' ) )  THEN
		    ccc = ginput ( :indall-1 )
		    CALL TG_FULL  ( ccc, firstm, lasttm, tbegin, ier )
		    IF  ( ier .ne. 0 )  THEN
			iret = -8
			RETURN
		    END IF
		    tbegin ( 12:12 ) = ctype
		    DO  i = 1, ntimf
			IF  ( ( filtim (i)(:12) .eq. tbegin (:12) ) .or.
     +			      ( ( filtim (i)(:11) .eq. tbegin (:11) )
     +						.and.
     +				( filtim (i)(12:12) .eq. ' ' ) .and.
     +				( tbegin (12:12) .eq. 'A' ) ) )  THEN
			    ntime = ntime + 1
			    times ( ntime ) = filtim (i)
			END IF
		    END DO
		    IF  ( ntime .eq. 0 )  iret = -16
		    RETURN
		  ELSE
		    iret = -7
		    RETURN
		END IF
	    END IF
	END IF
C
C*	Retrieve the first time.
C
	CALL TG_FULL  ( tstart, firstm, lasttm, tbegin, ier1 )
C
C*	Get the last time.
C
	IF  (  INDEX ( tstop, 'LAST' ) .eq. 0 )  THEN
	    CALL TG_FULL  ( tstop, firstm, tbegin, tend, ier2 )
	  ELSE
	    CALL TG_FULL  ( tstop, firstm, lasttm, tend, ier2 )
	END IF
C
C*	Check for errors.
C
	IF  ( ier1 .ne. 0 )  THEN
	    iret = -8
	    RETURN
	  ELSE IF  ( ier2 .ne. 0 )  THEN
	    iret = -9
	    RETURN
	END IF
C
C*	If first time is A00 time in an analysis range, add A00.
C
	IF  ( ( tbegin (12:12) .eq. ' ' ) .and.
     +	      ( tend   (12:12) .ne. ' ' ) )  tbegin (12: ) = 'A00'
C
C*	Check that first time is before last time.
C
	IF  ( tbegin .gt. tend )  THEN
C
C*	    Check for new century.
C
	    CALL TI_CCNT ( tbegin, cbegin, ier )
	    CALL TI_CCNT ( tend, cend, ier )
	    IF ( cbegin .ge. cend ) THEN
	        iret = -10
	        RETURN
	    END IF
	END IF
C
C*	Check that either time or forecast field is not changed.
C
	IF  ( tbegin ( : 11 ) .eq. tend ( : 11 ) )  THEN
	    IF  ( tbegin (12:12) .eq. tend (12:12) )  THEN
		itftyp = 1
	      ELSE
		iret = -12
		RETURN
	    END IF
	  ELSE IF  ( tbegin (12: ) .eq. tend (12: ) )  THEN
	    itftyp = 2
	  ELSE
	    iret = -13
	    RETURN
	END IF
C
C*	Call proper subroutine according to type.
C
	IF  ( itype .eq. 1 )  THEN
	    CALL TG_RNOI  ( tbegin, tend, itftyp, ntimf, filtim,
     +			    ntime, times, iret )
	  ELSE
	    CALL TG_RINC  ( tbegin, tend, tinc, itftyp, ntime, times,
     +			    iret )
	END IF
	IF  ( ntime .eq. 0 )  iret = -16
C*
	RETURN
	END

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
C    @(#)ldptigen.f	20.17 5/3/00
C****************************************************************
C
C     File: ldptigen.f
C
C     Purpose: Routine to load PTI generator data from raw data file
C
c     Return code:  n = 0 : Success
c                   n = 1 " Error
c
C     Author: Walt Powell  Date: 21 May 1996
C     Called by: load_pti.f
C
C****************************************************************
      integer function ldptigen (xbuf, options, numver, error)
      integer error, options(*), numver
      character *(*) xbuf

      include 'ipfinc/parametr.inc'

      include 'ipfinc/prt.inc'
      include 'ipfinc/pti_data.inc'
      include 'ipfinc/blank.inc'
      include 'ipfinc/alt_case.inc'
      include 'ipfinc/alpha.inc'
      include 'ipfinc/bus.inc'
      include 'ipfinc/cbus.inc'

      character word(30)*10, busname*8, type*1, cbtyp*1, cbown*2,
     &          cbkyr*2
      logical first
      integer find_bus, ftn_atoi, status, status1, gtptinam

      ldptigen = 0
      error = 0
      last = index (xbuf, '/') - 1
      if (last .le. 0) last = lastch(xbuf)
      call ptiscan (xbuf(1:last), word, nwrd, ',',  ' ')
c
c     Add defaults
c
      if (word(2) .eq. ' ') word(2) = '1'
      if (word(5) .eq. ' ') word(5) = '9999'
      if (word(6) .eq. ' ') word(6) = '-9999'
      if (word(7) .eq. ' ') word(7) = '1.000'
      if (word(15) .eq. ' ') word(15) = '1'

      status1 = gtptinam (word(1), numpti, num1, num2, busname, basekv)
      status = ftn_atoi(word(13))
      if (numpti .eq. 0 .and. busname .eq. ' ') then
        ldptigen = 1
      else if (status .eq. 0) then
      else
        numarea = 0
        numzone = 0
        if (numpti .ne. 0 .and. istatus .eq. 1) then
          write (errbuf(1), 100)
  100     format (' Generator record is not preceded with a Bus record
     &in raw data file.')
          errbuf(2)= ' (' // xbuf(1:60) // ')'
          call prterx ('W',2)
          error = 1
          go to 900
        endif
        num = find_bus (busname, basekv)
        if (num .lt. 0) then
          write (errbuf(1), 110) busname, basekv
  110     format (' Bus (', a8, f7.1, ') in Generator record is not in

     & system.')
          call prterx ('W', 1)
          error = 1
          go to 900
        endif
        if (kbsdta(1,num) .ne. 2 .and. kbsdta(1,num) .ne. 3 .and.
     &             kbsdta(1,num) .ne. 8) then
          call typno (type, kbsdta(1,num))
          write (errbuf(1), 120) busname, basekv, type
  120     format (' Generator bus (', a8, f7.1, ') is type B', a)
          call prterx ('W', 1)
          error = 1
          go to 900
        endif
        if (kbsdta(1,num) .ne. 3) then
          kbsdta(1,num) = 8                  ! Change to type BG
        endif
        busdta(7,num) = busdta(7,num) + ftn_atof(word(17))    ! Pmax
        if (word(2) .eq. '1') then
          busdta(8,num) = busdta(8,num) + ftn_atof(word(3))   ! PG
          busdta(9,num) = busdta(9,num) + ftn_atof(word(5))   ! Qmax
          busdta(10,num) = busdta(10,num) + ftn_atof(word(6)) ! Qmin
        else
          if (ntot2 .ge. MAXCBS) then
            write (errbuf(1), 130) MAXCBS
  130       format(' More than ', i5, ' + bus records entities.')
            errbuf(2) = ' (' // xbuf(1:60) // ')'
            call prterx ('W',2)
            error = 1
            ntot2 = 1
            go to 900
          endif
          ntot2 = ntot2 + 1
          kbctbl(1,ntot2) = num
          call linkcbus (ntot2, error)
          bctbl(2,ntot2) = 0.0
          bctbl(3,ntot2) = 0.0
          bctbl(4,ntot2) = 0.0
          bctbl(5,ntot2) = 0.0
          bctbl(6,ntot2) = ftn_atof(word(3))             ! PG
          kbctbl(7,ntot2) = 0
          cbown = owner(num)
          cbtyp = 'A'
          cbkyr = word(2)
          call putchr (1, cbtyp, kbctbl(8,ntot2))
          call putchr (2, cbkyr, kbctbl(9,ntot2))
          call putchr (3, cbown, kbctbl(10,ntot2))
          bctbl(11,ntot2) = ftn_atof(word(5))            ! Qmax
          bctbl(12,ntot2) = ftn_atof(word(6))            ! Qmin
        endif
        if (word(5) .eq. '9999') then
          if (kbsdta(1,num) .ne. 3) then
            kbsdta(1,num) = 2
          endif
        endif
        if (word(7) .ne. ' ') then
          busdta(11,num) = ftn_atof(word(7))             ! Vsched
          busdta(12,num) = busdta(11,num)                ! Vsched
        endif
        if (word(8) .ne. ' ') then
          istatus = gtptinam (word(8), numrpti, num1, num2, busname,
     &                        basekv)
          if (istatus .eq. 1 .or. numrpti .eq. 0) then
            nb = 0
          else
            nb = find_bus (busname, basekv)
            if (nb .lt. 0) nb = 0
          endif
          kbsdta(13,num) = nb                  ! Controlled bus
        else if (kbsdta(1,num) .eq. 8) then
          kbsdta(1,num) = 7                    ! Change to type BQ
          busdta(12,num) = 0.0
        endif
        pct = ftn_atof(word(16))                         ! %-var
        if (pct .gt. 0.0 .and. kbsdta(1,num) .ne. 8) then
          call typno (type, kbsdta(1,num))
          write (errbuf(1), 140) busname, basekv, type, pct
  140     format (' Generator (', a8, f7.1, ') is type B', a,
     &      ' and not BG. Percentage vars (', f7.1, ') ignored ')
          call prterx ('W', 1)
          error = 1
        else
          busdta(14,num) = pct
        endif
        call vltlim(num, vlimn(num), vlimx(num), vstart(num))
      endif
  900 continue
      return
      end

max_stars_repo_name: mbheinen/bpa-ipf-tsp

max_issues_repo_name: cuihantao/bpa-ipf-tsp

max_forks_repo_name: mbheinen/bpa-ipf-tsp


text:
cdis    Open Source License/Disclaimer, Forecast Systems Laboratory
cdis    NOAA/OAR/FSL, 325 Broadway Boulder, CO 80305
cdis
cdis    This software is distributed under the Open Source Definition,
cdis    which may be found at http://www.opensource.org/osd.html.
cdis
cdis    In particular, redistribution and use in source and binary forms,
cdis    with or without modification, are permitted provided that the
cdis    following conditions are met:
cdis
cdis    - Redistributions of source code must retain this notice, this
cdis    list of conditions and the following disclaimer.
cdis
cdis    - Redistributions in binary form must provide access to this
cdis    notice, this list of conditions and the following disclaimer, and
cdis    the underlying source code.
cdis
cdis    - All modifications to this software must be clearly documented,
cdis    and are solely the responsibility of the agent making the
cdis    modifications.
cdis
cdis    - If significant modifications or enhancements are made to this
cdis    software, the FSL Software Policy Manager
cdis    (softwaremgr@fsl.noaa.gov) should be notified.
cdis
cdis    THIS SOFTWARE AND ITS DOCUMENTATION ARE IN THE PUBLIC DOMAIN
cdis    AND ARE FURNISHED "AS IS."  THE AUTHORS, THE UNITED STATES
cdis    GOVERNMENT, ITS INSTRUMENTALITIES, OFFICERS, EMPLOYEES, AND
cdis    AGENTS MAKE NO WARRANTY, EXPRESS OR IMPLIED, AS TO THE USEFULNESS
cdis    OF THE SOFTWARE AND DOCUMENTATION FOR ANY PURPOSE.  THEY ASSUME
cdis    NO RESPONSIBILITY (1) FOR THE USE OF THE SOFTWARE AND
cdis    DOCUMENTATION; OR (2) TO PROVIDE TECHNICAL SUPPORT TO USERS.
cdis
cdis
      subroutine ofm (kk,laps_p,laps_t,laps_q,laps_sfc_t,
     1     psfc, jday, lat, ZA,
     1                                Tbest,
     1     radest,
     1     sec_za,              !optran 90 variable
     1     sfc_emis,            !optran 90 variable
     1     sfc_refl,            !optran 90 variable
     1     sec_solar            !optran 90 variable
     1     )

c   This routine interfaces GOES 8/10 satellite broadcast network data (and
c   local GVAR data) to the LAPS moisture analysis.  In 1999, this routine
c   was modified from an earlier version that used the University of
c   Wisconsin -- Madison's forward model to a newer model developed at
c   NESDIS.  OPTRAN (optical transmittance) forward model was developed by
c   Thomas Kleespies (NESDIS) and questions about this model should be
c   directed to him.  Forecast Systems Laboratory does not in any way
c   guarantee the validity of OPTRAN and distributes this software on an
c   as-is basis.  MOREOVER, FSL HAS PERMISSION TO DISTRIBUTE OPTRAN AS PART
c   OF LAPS TO FEDERAL AGENCIES.  NON-FEDERAL ENTITIES NEED TO INQUIRE WITH
c   NESDIS TO ESTABLISH THEIR RIGHTS AND OBLIGATIONS WITH REGARD TO OPTRAN.
c
c   The version of OPTRAN with which this software is used, has been
c   modified by FSL to include both sounder and imager channels for a
c   particular satellite in one call to the routine.  Thus a user only need
c   to setup OPTRAN for a particular satellite.  After doing such, either
c   the imager or sounding instrument can be used with the software without
c   further recompilation.

c      USE module_sfc_structure

      Implicit None

      save

      Include '../include/constants_optran.inc'

      include '../include/trans.inc'

      integer Mchan
      parameter (Mchan=18) ! adjusted for optran 90.

      integer kk                ! kk the laps profile dim,
      integer Nk                ! the size of the composite vectors
      integer start_level       !lowest level of climo used
      real laps_p(kk),laps_t(kk),laps_q(kk), lat, laps_sfc_t, psfc
      real ZA                   !zenith angle (degrees)
      real sec_za               !secant of ZA (for optran 90)
      real sfc_emis             !surface emissivity for optran 90
      real sfc_refl             !surface reflectivity for optran 90
      real sec_solar            !local secant of solar angle for optran 90
      real, dimension (:,:), allocatable :: tau90,flux_tau,
     1     solar_tau            !tau for optran 90 return variable
      real, dimension (:), allocatable ::  upwelling_radiance,
     1     brightness_temperature !upwelling radiance
      integer compute_rtm      ! optran 90 wrapper function is type int
      integer jday, istatus
      logical first_call
      data first_call /.true./

c     optran 90 variables

      real, dimension (0:Nlevel,1)::  level_p,level_t, level_w, level_o
      real, dimension (Nlevel,1)  ::  layer_p,layer_t, layer_w,layer_o
      integer nk_90 ! optran 90 reduced nk counter

c     climo variables
      real c_p (40)
      real c_t (40)
      real c_m (40)

      real conversion
      real mdf

c     satellite number

      integer goes_number
      common /sat_id/ goes_number

      Real O(Nlevel)

      Real Tbest(Nchan)         ! Brightness temps from estimated tau
      real radest(nchan)

      integer ichan

      Integer Channels(Mchan)
      data channels/1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18/

      integer i

c     optran 90 coefficient data information

      character*256 sndr_coeff, sndr_trans
      integer sndr_coeff_len, sndr_trans_len
      common /optn_coef/ sndr_coeff, sndr_trans, sndr_coeff_len,
     1     sndr_trans_len

c     ----------------------do ozone first call ----

      if (first_call) then
         first_call = .false.

c     determine climo start level

         call climate_sm(lat, jday, c_p, c_t, c_m, istatus)

         start_level = 40       ! (high c_p)
         do while (  c_p(start_level) .ge. laps_p(kk) )
            start_level = start_level - 1
         enddo

c     fill pressure level

         nk = 0

         do i = 1,start_level
            nk = nk +1
            p(nk)   = c_p(i)
         enddo

         do i = kk,1,-1
            nk      = nk+1
            p(nk)   = laps_p(i)
         enddo

         write(6,*) 'Ozone levels used = ',nk

c     grab ozone

         call oh_zone (p,o,nk,1,istatus)

c     convert (g/g) to (g/kg)

         conversion    = 1.0e+3

c     ozone called only one time

      endif

c------------------------------------------------

c     00000000000000000000000000000
c     assemble vectors module 1
c000000000000000000000000000000

c     build vectors, count total levels (not including surface)

c     determine climo start level

      start_level = 40          ! (high c_p)
      do while (  c_p(start_level) .ge. laps_p(kk) )
         start_level = start_level - 1
      enddo

      nk = 0

      do i = 1,start_level
         nk = nk +1
         p(nk)   = c_p(i)
         t(nk,1) = c_t(i)
         q(nk,1) = c_m(i)       ! note this is g/kg mixing ratio
      enddo

      do i = kk,1, -1
         if(laps_p(i) .lt. psfc) then ! accept
            nk      = nk+1
            p(nk)   = laps_p(i)
            t(nk,1) = laps_t(i)
            q(nk,1) = laps_q(i) * conversion ! note this is g/kg mixing ratio
         endif
      enddo

c     now at sfc or cloud top

      nk = nk + 1
      p(nk) = psfc
      t(nk,1) = laps_sfc_t
      q(nk,1) = q(nk-1,1)       ! approximation for now

      if(nk.gt.Nlevel) then
         write(6,*) 'Array dimension error'
         write(6,*) 'Module ofm.f'
         write(6,*) 'Parameter Nlevel too small'
      endif

c0000000000000000000000
c     end assembling vectors
c00000000000000000000000

cc insert here new code for OPTRAN 90 INTERFACE  dB  2002
c upper level layer computations.

      do i = 1,nk
         level_p (i-1,1) = p(i)
         level_t (i-1,1) = t(i,1)
         level_w (i-1,1) = q(i,1)
         level_o (i-1,1) = o(i)
      enddo
      nk_90 = nk - 1            ! (start counter at zero for layers)

c     generate layers

      call make_optran_layers ( level_t, level_w, level_o,
     1     layer_t, layer_w, layer_o, level_p,layer_p,
     1     nk_90, istatus)

c surface parameters

c     call to optran 90

c      if (nk .eq. 10000) then

         allocate (tau90 (1:nk_90, 1:mchan))
         allocate (flux_tau (1:nk_90,1:mchan))
         allocate (solar_tau (1:nk_90,1:mchan))
         allocate (upwelling_radiance (1:mchan))
         allocate (brightness_temperature (1:mchan))

         tau90 = 0.0
         flux_tau = 0.0
         solar_tau = 0.0
         upwelling_radiance = 0.0
         brightness_temperature = 0.0

         call optran90_fm (     ! call to optran 90 forward model
     1        nk_90,level_p, layer_p, layer_t, layer_w, layer_o,
     1        laps_sfc_t,
     1        sfc_emis,
     1        sfc_refl,
     1        sec_za,
     1        sec_solar,
     1        Mchan,
     1        tau90,
     1        flux_tau,
     1        solar_tau,
     1        upwelling_radiance,
     1        brightness_temperature,
     1        sndr_coeff,
     1        sndr_trans,
     1        sndr_coeff_len,
     1        sndr_trans_len
     1        )

         Do Ichan = 1 , Mchan

            TbEst (Ichan) = brightness_temperature (ichan)
            radest (ichan) = upwelling_radiance(ichan)

         EndDo

         deallocate (tau90,flux_tau)
         deallocate (upwelling_radiance, brightness_temperature)
         deallocate (solar_tau)

      return

      End

max_stars_repo_name: maxinye/laps-mirror

max_stars_repo_licenses: Intel Unlicense OLDAP-2.2.1 NetCDF

max_issues_repo_name: longwosion/laps-mirror

max_issues_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense

max_forks_repo_name: longwosion/laps-mirror

max_forks_repo_licenses: Intel NetCDF OLDAP-2.2.1 Unlicense


text:
C     Check that J is not privatized in the first loop, but privatized
C     in the second one. This is the same vase as privatization01.f, but
C     with a proper initialization of J

      subroutine privatization08
      integer a(10)

      j = 0
      do i = 1, 10
         j = j + 1
         a(i) = j
      enddo

      do i = 1, 10
         j = 1
         a(i) = j
      enddo
      end


text:
C Copyright (c) 2008 Sandia Corporation.  Under the terms of Contract
C DE-AC04-94AL85000 with Sandia Corporation, the U.S. Government retains
C certain rights in this software
C
C Redistribution and use in source and binary forms, with or without
C modification, are permitted provided that the following conditions are
C met:
C
C     * Redistributions of source code must retain the above copyright
C       notice, this list of conditions and the following disclaimer.
C
C     * Redistributions in binary form must reproduce the above
C       copyright notice, this list of conditions and the following
C       disclaimer in the documentation and/or other materials provided
C       with the distribution.
C
C     * Neither the name of Sandia Corporation nor the names of its
C       contributors may be used to endorse or promote products derived
C       from this software without specific prior written permission.
C
C THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
C "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
C LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR
C A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT
C OWNER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL,
C SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT
C LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE,
C DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY
C THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT
C (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE
C OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE.
C

C=======================================================================
      SUBROUTINE MUNELB (NELBLK, ISTAT, NUMEL,
     &   IDELB, NUMELB, NUMLNK, NUMATR,
     &   LINK, ATRIB, LINKX, ATRIBX, IXEL, IXELB, NELBX,
     $   ISCR, NAMELB, SCRSTR, LLINK, LATRIB)
C=======================================================================
C $Id: munelb.f,v 1.3 1999/03/22 15:53:00 gdsjaar Exp $
C $Log: munelb.f,v $
C Revision 1.3  1999/03/22 15:53:00  gdsjaar
C Fix problem with too many element blocks. GJoin has fixed size of 256
C element blocks, but didn't check whether maximum was exceeded.
C Changed so element block names are dynamically dimensioned.
C
C Revision 1.2  1999/01/27 15:12:36  gdsjaar
C Fixed handling of element order map. For gjoin, we force a 1..numel
C map on the combined mesh.
C
C Added writing of coordinate names
C
C Fixed error handling to not overwrite error code before call to exerr
C
C Revision 1.1  1999/01/18 19:21:23  gdsjaar
C ExodusII version of gjoin, needs testing and syncing with exodus 1 version, but is being committed to permit easier testing and modifications.  This was created by Dave Fry at Goodyear
C
c Revision 1.1.1.1  1998/11/05  16:23:26  a294617
c Initial import == gjoin 1.36
c
C Revision 1.2  1997/06/02 19:39:13  gdsjaar
C Fixed problem with link and attribute array sizes. The munelb routine
C now returns the new size of the arrays following element deletion and
C the arrays are resized using mdlong. Later calls to find the size of
C the link and attrib arrays using mdfind() now get the correct size.
C
C Revision 1.1.1.1  1990/11/12 14:35:10  gdsjaar
C GJOIN - X1.00.40 - 7/17/90
C
c Revision 1.1  90/11/12  14:35:08  gdsjaar
c Initial revision
c

C   --*** MUNELB *** (GJOIN) Compress and rearrange element blocks
C   --   Written by Amy Gilkey - revised 09/29/87
C   --   Modified by Greg Sjaardema, 07/11/90
C   --      Added element block names
C   --
C   --MUNELB processes the element blocks according to the block status.
C   --Blocks may be combined or deleted.
C   --
C   --Parameters:
C   --   NELBLK - IN/OUT - the number of element blocks
C   --   ISTAT - IN - the status of each block:
C   --      0 = same
C   --      - = delete
C   --      n = combine with block n
C   --   NUMEL - IN/OUT - the number of elements
C   --   IDELB - IN/OUT - the element block IDs for each block
C   --   NUMELB - IN/OUT - the number of elements in each block
C   --   NUMLNK - IN/OUT - the number of nodes per element in each block
C   --   NUMATR - IN/OUT - the number of attributes in each block
C   --   LINK - IN/OUT - the connectivity for each block
C   --   ATRIB - IN/OUT - the attributes for each block
C   --   LINKX - SCRATCH - sized to hold the new connectivity
C   --   ATRIBX - SCRATCH - sized to hold the new attributes
C   --   IXEL - OUT - the new element number for each element
C   --   IXELB - SCRATCH - size = NELBLK
C   --   NELBX - SCRATCH - size = NELBLK
C   --   ISCR - SCRATCH - size = NELBLK
C   --   NAMELB - IN/OUT - the names of the element blocks
C   --   SCRSTR - SCRATCH - size = size of NAMELB

      include 'exodusII.inc'

      INTEGER ISTAT(*)
      INTEGER IDELB(*)
      INTEGER NUMELB(*)
      INTEGER NUMLNK(*)
      INTEGER NUMATR(*)
      INTEGER LINK(*), LINKX(*)
      REAL ATRIB(*), ATRIBX(*)
      INTEGER IXEL(*)
      INTEGER IXELB(*)
      INTEGER NELBX(*)
      INTEGER ISCR(*)
      character*(MXSTLN) namelb(*), scrstr(*)

      DO 100 I = 1, NUMEL
         IXEL(I) = 0
  100 CONTINUE

      JEL0 = 0
      JLNK = 1
      JATR = 1

      JBLK = 0
      DO 140 IELB = 1, NELBLK
         IF (ISTAT(IELB) .EQ. 0) THEN
            NINSET = 1
            ISCR(NINSET) = IELB
         ELSE IF (ISTAT(IELB) .EQ. IELB) THEN
            CALL GETALL (IELB, NELBLK, ISTAT, NINSET, ISCR)
         ELSE
            NINSET = 0
         END IF

         IF (NINSET .GT. 0) THEN
            JBLK = JBLK + 1
            IXELB(JBLK) = IELB
            NELBX(JBLK) = 0
         END IF
         DO 130 ISET = 1, NINSET
            IBLK = ISCR(ISET)
            IEL0 = 0
            ILNK = 1
            IATR = 1
            DO 110 I = 1, IBLK-1
               IEL0 = IEL0 + NUMELB(I)
               ILNK = ILNK + NUMLNK(I) * NUMELB(I)
               IATR = IATR + NUMATR(I) * NUMELB(I)
  110       CONTINUE
            DO 120 I = 1, NUMELB(IBLK)
               IXEL(IEL0+I) = JEL0 + I
  120       CONTINUE
            CALL MOVINT (NUMELB(IBLK) * NUMLNK(IBLK),
     &         LINK(ILNK), LINKX(JLNK))
            CALL MOVREA (NUMELB(IBLK) * NUMATR(IBLK),
     &         ATRIB(IATR), ATRIBX(JATR))
            NELBX(JBLK) = NELBX(JBLK) + NUMELB(IBLK)
            JEL0 = JEL0 + NUMELB(IBLK)
            JLNK = JLNK + NUMLNK(IBLK) * NUMELB(IBLK)
            JATR = JATR + NUMATR(IBLK) * NUMELB(IBLK)
  130    CONTINUE
  140 CONTINUE

      CALL ORDIX  (JBLK, IXELB, NELBLK, IDELB, ISCR, IDELB)
      CALL MOVINT (JBLK, NELBX, NUMELB)
      CALL ORDIX  (JBLK, IXELB, NELBLK, NUMLNK, ISCR, NUMLNK)
      CALL ORDIX  (JBLK, IXELB, NELBLK, NUMATR, ISCR, NUMATR)
      CALL ORDSTR (JBLK, IXELB, NELBLK, NAMELB, SCRSTR, NAMELB)
      NELBLK = JBLK
      NUMEL  = JEL0
      LLINK  = JLNK-1
      LATRIB = JATR-1

      CALL MOVINT (LLINK,  LINKX,  LINK)
      CALL MOVREA (LATRIB, ATRIBX, ATRIB)

      RETURN
      END

max_stars_repo_name: milljm/seacas

max_issues_repo_name: milljm/seacas

max_forks_repo_name: milljm/seacas


module type_RealHarmonicCoefficients

    use spherepack_precision, only: &
        wp, & ! working precision
        ip ! integer precision

    ! Explicit typing only
    implicit none

    ! Everything is private unless stated otherwise
    private
    public :: RealHarmonicCoefficients

    type, public :: RealHarmonicCoefficients
        ! Type components
        logical,               public :: initialized = .false.
        integer(ip),           public :: NUMBER_OF_LONGITUDES = 0
        integer(ip),           public :: NUMBER_OF_LATITUDES = 0
        integer(ip),           public :: NUMBER_OF_SYNTHESES = 0
        real(wp), allocatable, public :: real_component(:, :)!, :)
        real(wp), allocatable, public :: imaginary_component(:, :)!, :)
    contains
        ! Type-bound procedures
        procedure, public  :: create => create_real_harmonic_coefficients
        procedure, public  :: destroy => destroy_real_harmonic_coefficients
        procedure, private :: copy_real_harmonic_coefficients
        ! Generic type-bound procedures
        generic, public :: assignment(=) => copy_real_harmonic_coefficients
    end type RealHarmonicCoefficients

    ! Declare user-defined constructor
    interface RealHarmonicCoefficients
        module procedure real_harmonic_coefficients_constructor
    end interface

contains

    function real_harmonic_coefficients_constructor(nlat, nlon, nt) &
        result (return_value)

        ! Dummy arguments
        integer(ip),           intent(in) :: nlat ! Number of latitudinal points 0 <= theta <= pi
        integer(ip),           intent(in) :: nlon ! Number of longitudinal points 0 <= phi <= 2*pi
        integer(ip), optional, intent(in) :: nt ! Number of syntheses
        type(RealHarmonicCoefficients)    :: return_value

        ! Local variables
        integer(ip) :: nt_op

        ! Address optional argument
        if (present(nt)) then
            nt_op = nt
        else
            nt_op = 1
        end if

        call return_value%create(nlat, nlon, nt_op)

    end function real_harmonic_coefficients_constructor

    subroutine create_real_harmonic_coefficients(self, nlat, nlon, nt)

        ! Dummy arguments
        class(RealHarmonicCoefficients), intent(inout) :: self
        integer(ip),                     intent(in)    :: nlat
        integer(ip),                     intent(in)    :: nlon
        integer(ip),                     intent(in)    :: nt

        ! Local variables
        integer(ip) :: mdab, ndab

        ! Ensure that object is usable
        call self%destroy()

        !  Set constants
        self%NUMBER_OF_LATITUDES = nlat
        self%NUMBER_OF_LONGITUDES = nlon
        self%NUMBER_OF_SYNTHESES = nt

                !  Set upper limit for vector m subscript
        !        select case(mod(nlon, 2))
        !            case(0)
        !                mdab = min(nlat, (nlon + 2)/2)
        !            case default
        !                mdab = min(nlat, (nlon + 1)/2)
        !        end select
        mdab = nlat
        ndab = nlat

        !  Allocate memory
        allocate (self%real_component(mdab, ndab))!, nt))
        allocate (self%imaginary_component(mdab, ndab))!, nt))

        ! Set flag
        self%initialized = .true.

    end subroutine create_real_harmonic_coefficients

    subroutine copy_real_harmonic_coefficients(self, other)

        ! Dummy arguments
        class(RealHarmonicCoefficients), intent(out) :: self
        class(RealHarmonicCoefficients), intent(in)  :: other

        ! Check if object is usable
        if (.not.other%initialized) then
            error stop 'Uninitialized object of class(RealHarmonicCoefficients): '&
                //'in assignment(=) '
        end if

        !  Make copies
        self%initialized = other%initialized
        self%NUMBER_OF_LATITUDES = other%NUMBER_OF_LATITUDES
        self%NUMBER_OF_LONGITUDES = other%NUMBER_OF_LONGITUDES
        self%NUMBER_OF_SYNTHESES = other%NUMBER_OF_SYNTHESES
        self%real_component = other%real_component
        self%imaginary_component = other%imaginary_component

    end subroutine copy_real_harmonic_coefficients

    subroutine destroy_real_harmonic_coefficients(self)

        ! Dummy arguments
        class(RealHarmonicCoefficients), intent(inout) :: self

        ! Check flag
        if (.not.self%initialized) return

        if (allocated(self%real_component)) then
            deallocate (self%real_component)
        end if

        if (allocated(self%imaginary_component)) then
            deallocate (self%imaginary_component)
        end if

        !  Reset constants
        self%NUMBER_OF_LONGITUDES = 0
        self%NUMBER_OF_LATITUDES = 0
        self%NUMBER_OF_SYNTHESES = 0

        ! Reset initialization flag
        self%initialized = .false.

    end subroutine destroy_real_harmonic_coefficients

end module type_RealHarmonicCoefficients


text:
!
!     CalculiX - A 3-dimensional finite element program
!              Copyright (C) 1998-2018 Guido Dhondt
!
!     This program is free software; you can redistribute it and/or
!     modify it under the terms of the GNU General Public License as
!     published by the Free Software Foundation(version 2);
!
!
!     This program is distributed in the hope that it will be useful,
!     but WITHOUT ANY WARRANTY; without even the implied warranty of
!     MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!     GNU General Public License for more details.
!
!     You should have received a copy of the GNU General Public License
!     along with this program; if not, write to the Free Software
!     Foundation, Inc., 675 Mass Ave, Cambridge, MA 02139, USA.
!
      subroutine creep(decra,deswa,statev,serd,ec,esw,p,qtild,
     &  temp,dtemp,predef,dpred,time,dtime,cmname,leximp,lend,
     &  coords,nstatv,noel,npt,layer,kspt,kstep,kinc)
!
!     user creep routine
!
!     INPUT (general):
!
!     statev(1..nstatv)  internal variables
!     serd               not used
!     ec(1)              equivalent creep at the start of the increment
!     ec(2)              not used
!     esw(1..2)          not used
!     p                  not used
!     temp               temperature at the end of the increment
!     dtemp              not used
!     predef             not used
!     dpred              not used
!     time(1)            value of the step time at the end of the increment
!     time(2)            value of the total time at the end of the increment
!     dtime              time increment
!     cmname             material name
!     leximp             not used
!     lend               if = 2: isotropic creep
!                        if = 3: anisotropic creep
!     coords(1..3)       coordinates of the current integration point
!     nstatv             number of internal variables
!     noel               element number
!     npt                integration point number
!     layer              not used
!     kspt               not used
!     kstep              not used
!     kinc               not used
!
!    INPUT only for elastic isotropic materials:
!     qtild              von Mises stress
!
!    INPUT only for elastic anisotropic materials:
!     decra(1)           equivalent deviatoric creep strain increment
!
!
!     OUTPUT (general):
!
!     decra(5)           derivative of the equivalent deviatoric
!                        creep strain increment w.r.t. the von Mises
!                        stress
!
!     OUTPUT only for elastic isotropic materials:
!     decra(1)           equivalent deviatoric creep strain increment
!
!     OUTPUT only for elastic anisotropic materials:
!     qtild              von Mises stress
!
      implicit none
!
      character*80 cmname
!
      integer leximp,lend,nstatv,noel,npt,layer,kspt,kstep,kinc
!
      real*8 decra(5),deswa(5),statev(*),serd,ec(2),esw(2),p,qtild,
     &  temp,dtemp,predef(*),dpred(*),time(2),dtime,coords(*)
!
      qtild=(1.d10*decra(1)/dtime)**0.2d0
      decra(5)=5.d-10*dtime*qtild**4
!
      return
      end

max_stars_repo_name: alleindrach/calculix-desktop

max_issues_repo_name: alleindrach/calculix-desktop

max_forks_repo_name: alleindrach/calculix-desktop


text:
program AlphaPhase
    use ProgramParametersModule
    use alphaFullChromModule
    use CompatibilityModule
    implicit none

    type(ProgramParameters) :: params

    ! Main is a seperate subroutine and called here due to compiler issues (Roberto / 16.0.3)
    call readParametersFile(params)

      if (params%plinkinputfile /= "") then
        call runPlink(params%plinkinputfile, params, runAlphaPhase)

    else
        call runAlphaPhase(params)
    endif
    contains

        subroutine readParametersFile(params)
            use ProgramParametersModule
            use InputOutput
            use AlphaPhaseResultsModule
            use PedigreeModule
            use HaplotypeLibraryModule
            use ProgramParametersModule
            use InputOutput
            use AlphaPhaseFunctions
            use HaplotypeModule

            type(ProgramParameters), intent(out) :: params
            type(PedigreeHolder) :: p
            integer, dimension(:,:), pointer :: CoreIndex
            character(len=4096) :: specfile, cmd
            integer :: specloc, nAnisG,nAnisP
            type(AlphaPhaseResults) :: results

            if (Command_Argument_Count() > 0) then
                call get_command_argument(1,cmd)
                if (cmd(1:2) .eq. "-v") then
                    call PrintVersion
                    call exit(0)
                end if
                if ((cmd(1:2) .eq. "-s") .or. (cmd(1:2) .eq. "-c")) then
                    call Titles
                    if (Command_Argument_Count() > 1) then
                        call Get_Command_Argument(2,specfile)
                    else
                        specfile="AlphaPhaseSpec.txt"
                    end if
                    params = ReadInParameterFile(specfile)
                    if (cmd(1:2) .eq. "-s") then
                        call CountInData(nAnisP, nAnisG, params)
                    CoreIndex => calculateCores(params%nSnp,params%params%Jump,params%params%offset)

                        print *
                        call printCoreInfo(CoreIndex)

                        ! Doesn't make sense to use Plink + specifing cores to run (at least for now) so set "plink" output directory
                        ! tot he current one
                        params%outputParams%resultFolderPath = "."
                        call makeDirectories(params%outputParams)

                        call writeCoreIndex(params%outputParams, size(CoreIndex,1), nAnisG, params%nSnp, CoreIndex(:,1), CoreIndex(:,2))
                        print *
                        print *, "Setup complete"
                        call exit(0)
                    end if
                    if (cmd(1:2) .eq. "-c") then
                        params%outputParams%outputPerCore = .false.
                        params%outputParams%outputCombined = .true.

                        ! The following are already done on a per core basis so will already have been created and aren't created right anyway.
                        ! Bit hacky and should probably change in future but works for now.
                        params%outputParams%outputHaplotypeLibraryText = .false.
                        params%outputParams%outputHaplotypeLibraryBinary = .false.
                        params%outputParams%outputSurrogatesSummary = .false.
                        params%outputParams%outputSurrogates = .false.
                        params%outputParams%outputHapCommonality = .false.
                        params%outputParams%outputIndivMistakes = .false.
                        params%outputParams%outputIndivMistakesPercent = .false.
                        params%outputParams%outputCoreMistakesPercent = .false.
                        ! As above
                        params%outputParams%resultFolderPath = "."

                        call ParsePedigreeAndGenotypeData(params, p)

                        ! Hack copied from below...
                        if (p%nHd == 0) then
                            p%nHd = p%nGenotyped
                            p%hdMap = p%genotypeMap

                            if (.not. allocated(p%hdDictionary)) then
                                allocate(p%hdDictionary)
                            endif
                            p%hdDictionary = p%genotypeDictionary
                        endif

                        call readInPerCoreResults(results, params%outputParams, p)
                        call writeAlphaPhaseResults(results, p, params%outputParams)
                        print *, "Results combined"
                        call exit(0)
                    end if
                end if
            end if

            call Titles

            if (Command_Argument_Count() > 0) then
                if (cmd(1:2) .eq. "-r") then
                    specloc = 3
                else
                    specloc = 1
                end if
                if (specloc <= Command_Argument_Count()) then
                    call Get_Command_Argument(specloc,specfile)
                else
                    specfile="AlphaPhaseSpec.txt"
                end if
            else
                specfile="AlphaPhaseSpec.txt"
            end if
            params = ReadInParameterFile(specfile)

        end subroutine readParametersFile

        subroutine runAlphaPhase(in, pedIn)
            use HaplotypeLibraryModule
            use ProgramParametersModule
            use InputOutput
            use baseSpecFileModule
            use AlphaPhaseResultsModule
            use AlphaPhaseFunctions
            use InputOutput
            use PedigreeModule
            use HaplotypeModule

            implicit none
            class(baseSpecFile),target :: in
            type(ProgramParameters),pointer :: params
            type(PedigreeHolder), target, optional :: pedIn
            type(PedigreeHolder) :: p
            type(Haplotype), pointer, dimension(:,:) :: Phase, TruePhase
            type(HaplotypeLibrary), dimension(:), pointer :: existingLibraries
            integer :: nAnisG

            logical :: notPrephased
            integer, dimension(:,:), pointer :: CoreIndex
            type(AlphaPhaseResults) :: results

            !Linux max path length is 4096 which is more than windows or mac (all according to google)

            character(len=4096) :: cmd
            character(len=4096) :: limit

            	select type(in)

                type is (ProgramParameters)
                params => in
                class default
                write(error_unit, *) "ERROR: AlphaPhase given incorrect object type as input"
                call abort()
            end select

            if (.not. present(pedIn)) then
                call ParsePedigreeAndGenotypeData(params, p)
            else
                p = pedIn
            endif

            if (p%nHd == 0) then
				p%nHd = p%nGenotyped
				p%hdMap = p%genotypeMap

				if (.not. allocated(p%hdDictionary)) then
					allocate(p%hdDictionary)
				endif
				p%hdDictionary = p%genotypeDictionary
			endif

            if (params%prePhaseFile .ne. "None") then
                call p%addPhaseInformationFromFile(params%prePhaseFile, p%nsnpsPopulation, p%nGenotyped)
            end if

            if (Command_Argument_Count() > 0) then
                call get_command_argument(1,cmd)
                if (cmd(1:2) .eq. "-r") then
                    call get_command_argument(2,limit)
                    if (index(limit,"-") > 0) then
                        params%params%startCoreChar = limit(1:index(limit,"-")-1)
                        params%params%endCoreChar = trim(limit(index(limit,"-")+1:4096))
                    else
                        params%params%startCoreChar = limit
                        params%params%endCoreChar = limit
                    end if
                    params%outputParams%outputCombined = .false.
                    params%outputParams%outputPerCore = .true.
                    params%outputParams%outputTimer = .false.
                    params%outputParams%outputGlobalCoreMistakesPercent = .false.
                end if
            end if

            nAnisG = p%nHd

            notPrephased = (params%GenotypeFileFormat /= 2)

            if (params%CoreFile .ne. "None") then
                CoreIndex => readInCores(params%CoreFile)
            end if

            if (notPrephased) then
                if (params%Simulation) then
                    TruePhase => ParsePhaseData(params%TruePhaseFile,nAnisG,params%nSnp)
                    if (params%Library .ne. "None") then
                        existingLibraries => getHaplotypeLibraries(params%library)
                        results = phaseAndCreateLibraries(p, params%params, existingLibraries, TruePhase, quiet = .false.)
                    else
                        if (params%CoreFile .ne. "None") then
                            results = phaseAndCreateLibraries(p, params%params, truePhase = TruePhase, userCoreIndex = coreIndex, quiet = .false.)
                        else
                            results = phaseAndCreateLibraries(p, params%params, truePhase = TruePhase, quiet = .false.)
                        end if
                    end if
                else
                    if (params%Library .ne. "None") then
                        existingLibraries => getHaplotypeLibraries(params%library)
                        results = phaseAndCreateLibraries(p, params%params, existingLibraries, quiet = .false.)
                    else
                        if (params%CoreFile .ne. "None") then
                            results = phaseAndCreateLibraries(p, params%params, userCoreIndex = coreIndex, quiet = .false.)
                        else
                            results = phaseAndCreateLibraries(p, params%params, quiet = .false.)
                        end if
                    end if
                end if
            else
                Phase => ParsePhaseData(params%GenotypeFile,nAnisG,params%nSnp)
                if (params%Library .ne. "None") then
                    existingLibraries => getHaplotypeLibraries(params%library)
                    results = createLibraries(Phase, params%params, existingLibraries)
                else
                    if (params%CoreFile .ne. "None") then
                        results = createLibraries(Phase, params%params, userCoreIndex = coreIndex)
                    else
                        results = createLibraries(Phase, params%params)
                    end if
                end if
            end if

            ! Copy the path from the Plink chromosome stuff to output params
            params%outputParams%resultFolderPath = params%resultFolderPath

            call writeAlphaPhaseResults(results, p, params%outputParams)

            call PrintTimerTitles(params%outputParams)
        end subroutine runAlphaPhase

    end program AlphaPhase

max_stars_repo_name: AlphaGenes/alphaphase

max_issues_repo_name: AlphaGenes/alphaphase

max_forks_repo_name: AlphaGenes/alphaphase


text:
program GenerateEvents
    !Run from the command line
    ! ./GenerateEvents M_X SIGMA_SI RHO
    ! where
    !   - M_X is the WIMP mass in GeV
    !   - SIGMA_SI is the WIMP-nucleon xsec in cm^2
    !   - RHO is the local density in GeV/cm^3

use ranlib_poisson
use LabFuncs
use DDrate
use modulation
use utils
use expt !E_min, E_max, t_start(), t_end(), A_det, m_det, resolution

implicit none

double precision :: N_exp
integer :: i, N_obs

double precision :: m_x, sigma_SI, rho

double precision, allocatable :: events_E(:), events_t(:)
double precision, allocatable :: reslist(:)

!RNG stuff
integer :: values(1:8), k
integer, dimension(:), allocatable :: seed
double precision, allocatable :: x(:)
double precision :: r

!Stuff for doing the inverse transform sampling
integer :: NE = 200
integer :: Nt
integer :: i_t, i_E
integer :: i_cut

!Lists for generating the inverse CDFs in (t, E)
double precision, allocatable :: tlist(:), Elist(:)
double precision, allocatable :: rate(:), cum_rate(:)

!Load in cross section and density from the command line
character(len=100) :: input_str
call getarg(1, input_str)
read(input_str, *) m_x

call getarg(2, input_str)
read(input_str,*) sigma_SI

call getarg(3, input_str)
read(input_str,*) rho

write(*,*) "Generating events using:"
write(*,*) "    m_x [GeV]:", sigma_SI
write(*,*) "    sigma_SI [cm^2]:", sigma_SI
write(*,*) "    rho [GeV/cm^3]:", rho

! Initialise seed for rng
call initialize()
call date_and_time(values=values)
call random_seed(size=k)
allocate(seed(1:k))
seed(:) = values(8)
call random_seed(put=seed)
call random_number(r)
call set_initial_seed ( floor(1d6*r), floor(1d7*r) )

!Read in the data tables for the daily modulations
call initialise_modulation

! Calculate total number of events
N_exp = (rho/rho0)*m_det*Nevents(E_min, E_max, t_start(), t_end(), A_det, m_x, sigma_SI)
write(*,*) "    N_exp:", N_exp
N_obs = IGNPOI(real(N_exp))
write(*,*) "    N_obs:", N_obs

allocate(events_E(N_obs))
allocate(events_t(N_obs))
allocate(x(N_obs))

!----------------
write(*,*) "    Sampling event times..."

!Calculate the rate using roughly 100 time steps per day
Nt = 99*nint(t_end() - t_start()) + 1

allocate(tlist(Nt))
allocate(rate(Nt))
allocate(cum_rate(Nt))

! Sampling event rate as a function of time
! integrated over energy
tlist = linspace(t_start(), t_end(), Nt)
do i_t = 1, Nt
    rate(i_t) = (rho/rho0)*m_det*Nevents_fixedt(E_min, E_max, tlist(i_t), A_det, m_x, sigma_SI)
end do

! Inverse transform sample to get event times
cum_rate =  cumtrapz(tlist, rate)
cum_rate = cum_rate/cum_rate(Nt)
call random_number(x)

events_t = linterp(cum_rate, tlist, x)

deallocate(tlist)
deallocate(rate)
deallocate(cum_rate)

!----------------
write(*,*) "    Sampling event energies..."

allocate(Elist(NE))
allocate(reslist(NE))
allocate(rate(NE))
allocate(cum_rate(NE))

Elist = 10**(linspace(log10(E_th), log10(E_max), NE))

do i = 1, N_obs
    do i_E = 1, NE
        rate(i_E) = (rho/rho0)*m_det*dRdE_res(Elist(i_E), events_t(i), A_det, m_x, sigma_SI)
    end do

    !Find the point where the rate drops to zero,
    !and cut there, otherwise the interpolation
    !can get screwed up...
    i_cut = 1
    do while (rate(i_cut) > 1d-30)
        i_cut = i_cut + 1
        if (i_cut > NE) then
            exit
        end if
    end do
    !In fact, cut one step earlier, just to be sure...
    i_cut = i_cut - 1

    cum_rate(1:i_cut) =  cumtrapz(Elist(1:i_cut), rate(1:i_cut))
    cum_rate(1:i_cut) = cum_rate(1:i_cut)/cum_rate(i_cut)

    call random_number(r)
    events_E(i) = linterp_scalar(cum_rate(1:i_cut), Elist(1:i_cut), r)

end do

!Save to file
write(*,*) "    Saving to file <events.txt>..."

open (unit = 2, file = "events.txt")
do i = 1, N_obs
    write (2, *) events_t(i), events_E(i)
end do
close(2)

end program GenerateEvents

max_stars_repo_name: bradkav/EarthScatterLikelihood

max_issues_repo_name: bradkav/EarthScatterLikelihood

max_forks_repo_name: bradkav/EarthScatterLikelihood


text:
C sketch
C        metropolis hasting monte carlo
C        hamiltonian monte carlo
C not really gonna use it, consider there are
C MCMCpack and STAN
C rather for fun and jot down thoughts
C
C   arXiv:1808.08490v2
C   MikeXL/bayes/R/mcmc.r
C
C metropolis hasting, the simplest one
C target the logged posterior function
C nmc number iterations
C nbi burnin
C chain the Markov chain
C23456789
      subroutine metrop(target, nmc, nbi, chain)
        call random_number(chain(1, :))
        do 71 i = 2, nmc
          x = chain(i-1, :)
          call random_number(chain(i, :)
          call random_number(metropolis)
          h = exp(target(chain(i, :)) - target(x))
          if h > metropolis then
C accept
          else
C reject
            chain(i, :) = x
          end
 71     continue
      end subroutine metrop
C
C Hamiltonian, put out your sketchbook
C derivatives
C
      subroutine ham()
        do 99 i = 1, nmc
C where the magic happens
          call molecular_dynamics(..., Hi, Hf)
          h = exp(Hf - Hi)
          call random_number(metropolis)
          if h > metropolis then
C
          else
C
          end if
 99     continue
      end subroutine ham

max_stars_repo_name: MikeXL/bayes

max_issues_repo_name: MikeXL/bayes

max_forks_repo_name: MikeXL/bayes


text:
      SUBROUTINE SB03SX( TRANA, UPLO, LYAPUN, N, XANORM, T, LDT, U, LDU,
     $                   R, LDR, FERR, IWORK, DWORK, LDWORK, INFO )
C
C     SLICOT RELEASE 5.7.
C
C     Copyright (c) 2002-2020 NICONET e.V.
C
C     PURPOSE
C
C     To estimate a forward error bound for the solution X of a real
C     discrete-time Lyapunov matrix equation,
C
C            op(A)'*X*op(A) - X = C,
C
C     where op(A) = A or A' (A**T) and C is symmetric (C = C**T). The
C     matrix A, the right hand side C, and the solution X are N-by-N.
C     An absolute residual matrix, which takes into account the rounding
C     errors in forming it, is given in the array R.
C
C     ARGUMENTS
C
C     Mode Parameters
C
C     TRANA   CHARACTER*1
C             Specifies the form of op(A) to be used, as follows:
C             = 'N':  op(A) = A    (No transpose);
C             = 'T':  op(A) = A**T (Transpose);
C             = 'C':  op(A) = A**T (Conjugate transpose = Transpose).
C
C     UPLO    CHARACTER*1
C             Specifies which part of the symmetric matrix R is to be
C             used, as follows:
C             = 'U':  Upper triangular part;
C             = 'L':  Lower triangular part.
C
C     LYAPUN  CHARACTER*1
C             Specifies whether or not the original Lyapunov equations
C             should be solved, as follows:
C             = 'O':  Solve the original Lyapunov equations, updating
C                     the right-hand sides and solutions with the
C                     matrix U, e.g., X <-- U'*X*U;
C             = 'R':  Solve reduced Lyapunov equations only, without
C                     updating the right-hand sides and solutions.
C
C     Input/Output Parameters
C
C     N       (input) INTEGER
C             The order of the matrices A and R.  N >= 0.
C
C     XANORM  (input) DOUBLE PRECISION
C             The absolute (maximal) norm of the symmetric solution
C             matrix X of the Lyapunov equation.  XANORM >= 0.
C
C     T       (input) DOUBLE PRECISION array, dimension (LDT,N)
C             The leading N-by-N upper Hessenberg part of this array
C             must contain the upper quasi-triangular matrix T in Schur
C             canonical form from a Schur factorization of A.
C
C     LDT     INTEGER
C             The leading dimension of array T.  LDT >= MAX(1,N).
C
C     U       (input) DOUBLE PRECISION array, dimension (LDU,N)
C             The leading N-by-N part of this array must contain the
C             orthogonal matrix U from a real Schur factorization of A.
C             If LYAPUN = 'R', the array U is not referenced.
C
C     LDU     INTEGER
C             The leading dimension of array U.
C             LDU >= 1,        if LYAPUN = 'R';
C             LDU >= MAX(1,N), if LYAPUN = 'O'.
C
C     R       (input/output) DOUBLE PRECISION array, dimension (LDR,N)
C             On entry, if UPLO = 'U', the leading N-by-N upper
C             triangular part of this array must contain the upper
C             triangular part of the absolute residual matrix R, with
C             bounds on rounding errors added.
C             On entry, if UPLO = 'L', the leading N-by-N lower
C             triangular part of this array must contain the lower
C             triangular part of the absolute residual matrix R, with
C             bounds on rounding errors added.
C             On exit, the leading N-by-N part of this array contains
C             the symmetric absolute residual matrix R (with bounds on
C             rounding errors added), fully stored.
C
C     LDR     INTEGER
C             The leading dimension of array R.  LDR >= MAX(1,N).
C
C     FERR    (output) DOUBLE PRECISION
C             An estimated forward error bound for the solution X.
C             If XTRUE is the true solution, FERR bounds the magnitude
C             of the largest entry in (X - XTRUE) divided by the
C             magnitude of the largest entry in X.
C             If N = 0 or XANORM = 0, FERR is set to 0, without any
C             calculations.
C
C     Workspace
C
C     IWORK   INTEGER array, dimension (N*N)
C
C     DWORK   DOUBLE PRECISION array, dimension (LDWORK)
C
C     LDWORK  INTEGER
C             The length of the array DWORK.
C             LDWORK >= 0,            if N = 0;
C             LDWORK >= MAX(3,2*N*N), if N > 0.
C
C     Error Indicator
C
C     INFO    INTEGER
C             = 0:  successful exit;
C             < 0:  if INFO = -i, the i-th argument had an illegal
C                   value;
C             = N+1:  if T has almost reciprocal eigenvalues; perturbed
C                   values were used to solve Lyapunov equations (but
C                   the matrix T is unchanged).
C
C     METHOD
C
C     The forward error bound is estimated using a practical error bound
C     similar to the one proposed in [1], based on the 1-norm estimator
C     in [2].
C
C     REFERENCES
C
C     [1] Higham, N.J.
C         Perturbation theory and backward error for AX-XB=C.
C         BIT, vol. 33, pp. 124-136, 1993.
C
C     [2] Higham, N.J.
C         FORTRAN codes for estimating the one-norm of a real or
C         complex matrix, with applications to condition estimation.
C         ACM Trans. Math. Softw., 14, pp. 381-396, 1988.
C
C     NUMERICAL ASPECTS
C                               3
C     The algorithm requires 0(N ) operations.
C
C     FURTHER COMMENTS
C
C     The option LYAPUN = 'R' may occasionally produce slightly worse
C     or better estimates, and it is much faster than the option 'O'.
C     The routine can be also used as a final step in estimating a
C     forward error bound for the solution of a discrete-time algebraic
C     matrix Riccati equation.
C
C     CONTRIBUTOR
C
C     V. Sima, Research Institute for Informatics, Bucharest, Romania,
C     Oct. 1998. Partly based on DDLSVX (and then SB03SD) by P. Petkov,
C     Tech. University of Sofia, March 1998 (and December 1998).
C
C     REVISIONS
C
C     February 6, 1999, V. Sima, Katholieke Univ. Leuven, Belgium.
C     V. Sima, Research Institute for Informatics, Bucharest, Oct. 2004,
C     May 2020.
C
C     KEYWORDS
C
C     Lyapunov equation, orthogonal transformation, real Schur form.
C
C     ******************************************************************
C
C     .. Parameters ..
      DOUBLE PRECISION   ZERO, ONE, HALF
      PARAMETER          ( ZERO = 0.0D+0, ONE = 1.0D+0, HALF = 0.5D+0 )
C     ..
C     .. Scalar Arguments ..
      CHARACTER          LYAPUN, TRANA, UPLO
      INTEGER            INFO, LDR, LDT, LDU, LDWORK, N
      DOUBLE PRECISION   FERR, XANORM
C     ..
C     .. Array Arguments ..
      INTEGER            IWORK( * )
      DOUBLE PRECISION   DWORK( * ), R( LDR, * ), T( LDT, * ),
     $                   U( LDU, * )
C     ..
C     .. Local Scalars ..
      LOGICAL            LOWER, NOTRNA, UPDATE
      CHARACTER          TRANAT, UPLOW
      INTEGER            I, IJ, INFO2, ITMP, J, KASE, NN
      DOUBLE PRECISION   EST, SCALE, TEMP
C     ..
C     .. Local Arrays ..
      INTEGER            ISAVE( 3 )
C     ..
C     .. External Functions ..
      LOGICAL            LSAME
      DOUBLE PRECISION   DLANSY
      EXTERNAL           DLANSY, LSAME
C     ..
C     .. External Subroutines ..
      EXTERNAL           DLACN2, DSCAL, MA02ED, MB01RU, SB03MX, XERBLA
C     ..
C     .. Intrinsic Functions ..
      INTRINSIC          MAX
C     ..
C     .. Executable Statements ..
C
C     Decode and Test input parameters.
C
      NOTRNA = LSAME( TRANA,  'N' )
      UPDATE = LSAME( LYAPUN, 'O' )
C
      NN   = N*N
      INFO = 0
      IF( .NOT.( NOTRNA .OR. LSAME( TRANA, 'T' ) .OR.
     $                       LSAME( TRANA, 'C' ) ) ) THEN
         INFO = -1
      ELSE IF( .NOT.( LSAME( UPLO, 'L' ) .OR. LSAME( UPLO, 'U' ) ) )
     $   THEN
         INFO = -2
      ELSE IF( .NOT.( UPDATE .OR. LSAME( LYAPUN, 'R' ) ) ) THEN
         INFO = -3
      ELSE IF( N.LT.0 ) THEN
         INFO = -4
      ELSE IF( XANORM.LT.ZERO ) THEN
         INFO = -5
      ELSE IF( LDT.LT.MAX( 1, N ) ) THEN
         INFO = -7
      ELSE IF( LDU.LT.1 .OR. ( UPDATE .AND. LDU.LT.N ) ) THEN
         INFO = -9
      ELSE IF( LDR.LT.MAX( 1, N ) ) THEN
         INFO = -11
      ELSE IF( LDWORK.LT.0 .OR.
     $       ( LDWORK.LT.MAX( 3, 2*NN ) .AND. N.GT.0 ) ) THEN
         INFO = -15
      END IF
C
      IF( INFO.NE.0 ) THEN
         CALL XERBLA( 'SB03SX', -INFO )
         RETURN
      END IF
C
C     Quick return if possible.
C
      FERR = ZERO
      IF( N.EQ.0 .OR. XANORM.EQ.ZERO )
     $   RETURN
C
      ITMP = NN + 1
C
      IF( NOTRNA ) THEN
         TRANAT = 'T'
      ELSE
         TRANAT = 'N'
      END IF
C
C     Fill in the remaining triangle of the symmetric residual matrix.
C
      CALL MA02ED( UPLO, N, R, LDR )
C
      KASE = 0
C
C     REPEAT
   10 CONTINUE
      CALL DLACN2( NN, DWORK( ITMP ), DWORK, IWORK, EST, KASE, ISAVE )
      IF( KASE.NE.0 ) THEN
C
C        Select the triangular part of symmetric matrix to be used.
C
         IF( DLANSY( '1-norm', 'Upper', N, DWORK, N, DWORK( ITMP ) )
     $       .GE.
     $       DLANSY( '1-norm', 'Lower', N, DWORK, N, DWORK( ITMP ) )
     $     ) THEN
            UPLOW = 'U'
            LOWER = .FALSE.
         ELSE
            UPLOW = 'L'
            LOWER = .TRUE.
         END IF
C
         IF( KASE.EQ.2 ) THEN
            IJ = 0
            IF( LOWER ) THEN
C
C              Scale the lower triangular part of symmetric matrix
C              by the residual matrix.
C
               DO 30 J = 1, N
                  DO 20 I = J, N
                     IJ = IJ + 1
                     DWORK( IJ ) = DWORK( IJ )*R( I, J )
   20             CONTINUE
                  IJ = IJ + J
   30          CONTINUE
            ELSE
C
C              Scale the upper triangular part of symmetric matrix
C              by the residual matrix.
C
               DO 50 J = 1, N
                  DO 40 I = 1, J
                     IJ = IJ + 1
                     DWORK( IJ ) = DWORK( IJ )*R( I, J )
   40             CONTINUE
                  IJ = IJ + N - J
   50          CONTINUE
            END IF
         END IF
C
         IF( UPDATE ) THEN
C
C           Transform the right-hand side: RHS := U'*RHS*U.
C
            CALL MB01RU( UPLOW, 'Transpose', N, N, ZERO, ONE, DWORK, N,
     $                   U, LDU, DWORK, N, DWORK( ITMP ), NN, INFO2 )
            CALL DSCAL( N, HALF, DWORK, N+1 )
         END IF
         CALL MA02ED( UPLOW, N, DWORK, N )
C
         IF( KASE.EQ.2 ) THEN
C
C           Solve op(T)'*Y*op(T) - Y = scale*RHS.
C
            CALL SB03MX( TRANA, N, T, LDT, DWORK, N, SCALE,
     $                   DWORK( ITMP ), INFO2 )
         ELSE
C
C           Solve op(T)*W*op(T)' - W = scale*RHS.
C
            CALL SB03MX( TRANAT, N, T, LDT, DWORK, N, SCALE,
     $                   DWORK( ITMP ), INFO2 )
         END IF
C
         IF( INFO2.GT.0 )
     $      INFO = N + 1
C
         IF( UPDATE ) THEN
C
C           Transform back to obtain the solution: Z := U*Z*U', with
C           Z = Y or Z = W.
C
            CALL MB01RU( UPLOW, 'No transpose', N, N, ZERO, ONE, DWORK,
     $                   N, U, LDU, DWORK, N, DWORK( ITMP ), NN, INFO2 )
            CALL DSCAL( N, HALF, DWORK, N+1 )
         END IF
C
         IF( KASE.EQ.1 ) THEN
            IJ = 0
            IF( LOWER ) THEN
C
C              Scale the lower triangular part of symmetric matrix
C              by the residual matrix.
C
               DO 70 J = 1, N
                  DO 60 I = J, N
                     IJ = IJ + 1
                     DWORK( IJ ) = DWORK( IJ )*R( I, J )
   60             CONTINUE
                  IJ = IJ + J
   70          CONTINUE
            ELSE
C
C              Scale the upper triangular part of symmetric matrix
C              by the residual matrix.
C
               DO 90 J = 1, N
                  DO 80 I = 1, J
                     IJ = IJ + 1
                     DWORK( IJ ) = DWORK( IJ )*R( I, J )
   80             CONTINUE
                  IJ = IJ + N - J
   90          CONTINUE
            END IF
         END IF
C
C        Fill in the remaining triangle of the symmetric matrix.
C
         CALL MA02ED( UPLOW, N, DWORK, N )
         GO TO 10
      END IF
C
C     UNTIL KASE = 0
C
C     Compute the estimate of the relative error.
C
      TEMP = XANORM*SCALE
      IF( TEMP.GT.EST ) THEN
         FERR = EST / TEMP
      ELSE
         FERR = ONE
      END IF
C
      RETURN
C
C *** Last line of SB03SX ***
      END

max_stars_repo_name: bnavigator/SLICOT-Reference

max_issues_repo_name: RJHKnight/slicotr

max_forks_repo_name: RJHKnight/slicotr


!----------------------------------------------------------------------
!----------------------------------------------------------------------
!   ezp :    Complex bifurcation in a boundary value problem
!----------------------------------------------------------------------
!----------------------------------------------------------------------

      SUBROUTINE FUNC(NDIM,U,ICP,PAR,IJAC,F,DFDU,DFDP)
!     ---------- ----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM, ICP(*), IJAC
      COMPLEX(KIND(1.0D0)), INTENT(IN) :: U(NDIM)
      DOUBLE PRECISION, INTENT(IN) :: PAR(*)
      COMPLEX(KIND(1.0D0)), INTENT(OUT) :: F(NDIM)
      COMPLEX(KIND(1.0D0)), INTENT(INOUT) :: DFDU(NDIM,NDIM), DFDP(NDIM,*)

      COMPLEX(KIND(1.0D0)) :: U1, U2, RL, E

       U1=U(1)
       U2=U(2)

       RL=PAR(1)

       E=CDEXP(U1)
       F(1)=U2
       F(2)=-RL*E

      END SUBROUTINE FUNC

      SUBROUTINE STPNT(NDIM,U,PAR,T)
!     ---------- -----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM
      COMPLEX(KIND(1.0D0)), INTENT(INOUT) :: U(NDIM)
      DOUBLE PRECISION, INTENT(INOUT) :: PAR(*)
      DOUBLE PRECISION, INTENT(IN) :: T

       U(1)=0
       U(2)=0

      END SUBROUTINE STPNT

      SUBROUTINE BCND(NDIM,PAR,ICP,NBC,U0,U1,F,IJAC,DBC)
!     ---------- ----

      IMPLICIT NONE
      INTEGER, INTENT(IN) :: NDIM, ICP(*), NBC, IJAC
      DOUBLE PRECISION, INTENT(IN) :: PAR(*)
      COMPLEX(KIND(1.0D0)), INTENT(IN) :: U0(NDIM), U1(NDIM)
      COMPLEX(KIND(1.0D0)), INTENT(OUT) :: F(NBC)
      COMPLEX(KIND(1.0D0)), INTENT(INOUT) :: DBC(NBC,*)

       F(1)=U0(1)
       F(2)=U1(1)

      END SUBROUTINE BCND

      SUBROUTINE ICND
      END SUBROUTINE ICND

      SUBROUTINE FOPT
      END SUBROUTINE FOPT

      SUBROUTINE PVLS
      END SUBROUTINE PVLS


text:
************************************************************************
* Build B-matrices
*
* Build the B-matrix blocks in the system matrix by exact evaluation.
* This is possible due to the fact that the pressure basis functions
* are piecewise constant and the velocity basis functions
* piecewise linear. So there is an exact quadrature rule!
*
* Works only if Q1~ is used!
*
* In:
*   NEL,
*   NVT,
*   NMT,
*   KVERT,
*   KMID,
*   KADJ,
*   DCORVG  - usual geometry information
*   B1      - array [1..2*NMT] of double
*             Memory for first B-matrix block
*   B2      - array [1..2*NMT] of double
*             Memory for 2nd B-matrix block
*   KCOLB   - array [1..2*NMT] of double
*             Memory for column structure
*   KLDB    - array [1..NMT+1] of double
*             Memory for row structure structure
*
* Out:
*   NB      - Number of entries in the B-matrices
*   B1      - The first NB entries are filled with data
*   B2      - The first NB entries are filled with data
*   KCOL    - The first NB entries are filled with data
*   B1      - Is filled with data
*
* After calling this routine, if the caller allocated memory for the
* B-matrices dynamically, the memory for B1,B2,KCOL can be released
* except for the actual NB entries.
************************************************************************

      SUBROUTINE BBUILD(KVERT,KMID,KADJ,DCORVG,B1,B2,KCOLB,KLDB,
     *                  NB,NEL,NVT,NMT)

      IMPLICIT NONE

C Common blocks

      INCLUDE 'cbasictria.inc'
      INCLUDE 'cout.inc'

C parameters

      INTEGER KVERT(NNVE,*),KMID(NNVE,*),KADJ(NNVE,*)
      DOUBLE PRECISION DCORVG(2,*)
      INTEGER KCOLB(*),KLDB(*)

      INTEGER NB,NEL,NVT,NMT

C local variables

      INTEGER ILD, IADJ, IVEN, IMID, IVT, IVTN, IEL, IVE

      DOUBLE PRECISION B1(*),B2(*)
      DOUBLE PRECISION PX, PY, PXN, PYN, DN1, DN2

      ILD=0

C Loop through all elements

      DO IEL=1,NEL

C On each element loop through the vertices

        DO IVE=1,NNVE

C What's the neighbour element next to the edge
C starting with our current vertex?

          IADJ=KADJ(IVE,IEL)
          IF ((IADJ.LE.0).OR.(IADJ.GE.IEL)) THEN

C What's the next vertex = endpoint of our current edge?

            IVEN=IVE+1
            IF (IVEN.EQ.5) IVEN=1

C Calculate the number of our current vertex, the neighbour
C vertex and the midpoint inbetween...

            IMID=KMID (IVE ,IEL)-NVT
            IVT =KVERT(IVE ,IEL)
            IVTN=KVERT(IVEN,IEL)

C Then the coordinates of the two vertices

            PX =DCORVG(1,IVT)
            PY =DCORVG(2,IVT)
            PXN=DCORVG(1,IVTN)
            PYN=DCORVG(2,IVTN)

C and the length of the edge in X- and Y-direction, resp.

            DN1=-PYN+PY
            DN2= PXN-PX

C calculate the value in the B-matrix in the current row
C by exact quadrature

            ILD=ILD+1
            KLDB (IMID)=ILD
            KCOLB(ILD )=IEL
            B1   (ILD )=DN1
            B2   (ILD )=DN2

            IF (IADJ.GT.0) THEN
              ILD=ILD+1
              KCOLB(ILD)= IADJ
              B1   (ILD)=-DN1
              B2   (ILD)=-DN2
            END IF

          END IF ! ((IADJ.LE.0).OR.(IADJ.GE.IEL))

        END DO ! IVE

      END DO ! IEL

      NB         =ILD
      KLDB(NMT+1)=ILD+1

      END

************************************************************************
* Matrix multiplication routine with B-matrix, matrix-free
*
* Build the matrix-vector-product with the B-matrices elementwise,
* i.e. without existing B-matrices. This realizes IPRECB=2.
*
* Can only be used with Q1~/Q0 !
*
* In:
*   NEL,NVT,NMT,
*   KVERT,
*   KMID,
*   KADJ,
*   DCORVG - usual geometry information
*   DP     - array [1..NEL] of double
*            Pressure vector
*   DF1    - array [1..NMT] of double
*            Velocity/RHS vector, 1st component
*   DF2    - array [1..NMT] of double
*            Velocity/RHS vector, 2nd component
*   A1     - multiplication factor for DF1
*   A2     - multiplication factor for DF2
*
* Out:
*   DFx = A1*Bx*DP + A2*DFx
************************************************************************

      SUBROUTINE BMUL1 (KVERT,KMID,KADJ,DCORVG,DP,DF1,DF2,
     *                  NEL,NVT,NMT,A1,A2)

      IMPLICIT NONE

C Common blocks

      INCLUDE 'cbasictria.inc'

C parameters

      INTEGER KVERT(NNVE,*),KMID(NNVE,*),KADJ(NNVE,*)
      DOUBLE PRECISION DCORVG(2,*)
      INTEGER NEL,NVT,NMT
      DOUBLE PRECISION DP(NEL),DF1(NMT),DF2(NMT),A1,A2

C local variables

      INTEGER IEL, IVE, IVEN, IMID, IVT, IVTN, IADJ
      DOUBLE PRECISION PX, PY, PXN, PYN, DN1, DN2

C     Loop over the elements

      DO IEL=1,NEL

C       Loop over the vertices

        DO IVE=1,NNVE

          IADJ=KADJ(IVE,IEL)

          IF ((IADJ.LE.0).OR.(IADJ.GE.IEL)) THEN

C           Build the entry of the B-matrix by exact quadrature
C           (which luckily exists for Q1~/Q0), like in BBUILD above:

            IVEN=IVE+1
            IF (IVEN.EQ.5) IVEN=1

            IMID=KMID (IVE ,IEL)-NVT
            IVT =KVERT(IVE ,IEL)
            IVTN=KVERT(IVEN,IEL)

            PX =DCORVG(1,IVT)
            PY =DCORVG(2,IVT)
            PXN=DCORVG(1,IVTN)
            PYN=DCORVG(2,IVTN)

            DN1=-PYN+PY
            DN2= PXN-PX

C           Calculate the matrix-vector product

            IF (IADJ.GT.0) THEN
              DF1(IMID)=A2*DF1(IMID)+A1*DN1*(DP(IEL)-DP(IADJ))
              DF2(IMID)=A2*DF2(IMID)+A1*DN2*(DP(IEL)-DP(IADJ))
            ELSE
              DF1(IMID)=A2*DF1(IMID)+A1*DN1*DP(IEL)
              DF2(IMID)=A2*DF2(IMID)+A1*DN2*DP(IEL)
            ENDIF

          END IF ! ((IADJ.LE.0).OR.(IADJ.GE.IEL))

        END DO ! IVE

      END DO ! IEL

      END

************************************************************************
* Matrix multiplication routine. Multiply with B^T, matrix-free.
*
* Build the matrix-vector-product with the B^T-matrices elementwise,
* i.e. without existing B-matrices. This realizes IPRECB=2.
*
* Can only be used with Q1~/Q0 !
*
* In:
*   NEL,NVT,NMT,
*   KVERT,
*   KMID,
*   KADJ,
*   DCORVG - usual geometry information
*   DU1    - array [1..NMT] of double
*            Velocity/RHS vector, 1st component
*   DU2    - array [1..NMT] of double
*            Velocity/RHS vector, 2nd component
*   A1     - multiplication factor for DF1
*
* Out:
*   DFP    - array [1..NEL] of double
*            DFP  =  A1*B^T*DU  =  A1 * (B1^T*DU1 + B2^T*DU2)
************************************************************************

      SUBROUTINE BTMUL1(KVERT,KMID,KADJ,DCORVG,DU1,DU2,DFP,
     *                  NEL,NVT,NMT,A1)

      IMPLICIT NONE

C Common blocks

      INCLUDE 'cbasictria.inc'

C parameters

      INTEGER KVERT(NNVE,*),KMID(NNVE,*),KADJ(NNVE,*)
      DOUBLE PRECISION DCORVG(2,*)
      INTEGER NEL,NVT,NMT
      DOUBLE PRECISION DU1(NMT),DU2(NMT),DFP(NEL),A1

C local variables

      INTEGER IEL, IVE, IVEN, IMID, IVT, IVTN
      DOUBLE PRECISION PX, PY, PXN, PYN, DN1, DN2

C     Loop over the elements

      DO IEL=1,NEL

        DFP(IEL)=0D0

C       Loop over the vertices

        DO IVE=1,NNVE

C         Build the matrix entry as in BBUILD:

          IVEN=IVE+1
          IF (IVEN.EQ.5) IVEN=1

          IMID=KMID (IVE ,IEL)-NVT
          IVT =KVERT(IVE ,IEL)
          IVTN=KVERT(IVEN,IEL)

          PX =DCORVG(1,IVT)
          PY =DCORVG(2,IVT)
          PXN=DCORVG(1,IVTN)
          PYN=DCORVG(2,IVTN)

          DN1=-PYN+PY
          DN2= PXN-PX

          DFP(IEL)=DFP(IEL)+A1*(DN1*DU1(IMID)+DN2*DU2(IMID))

        END DO ! IVE

      END DO ! IEL

      END

max_stars_repo_licenses: Intel Unlicense

max_issues_repo_licenses: Intel Unlicense

max_forks_repo_licenses: Intel Unlicense


text:
      SUBROUTINE RSURFPLT
C
C CHANGE RECORD
C **  SUBROUTINE SURFPLT WRITES FILES TO CONTOUR FREE SURFACE
C **  ELEVATION
C
      USE GLOBAL
      CHARACTER*80 TITLE
      IF(JSRPPH.NE.1) GOTO 300
      OPEN(10,FILE='RSURFCN.OUT',STATUS='UNKNOWN')
      CLOSE(10,STATUS='DELETE')
      OPEN(10,FILE='RSURFCN.OUT',STATUS='UNKNOWN')
      TITLE='INSTANTANEOUS SURFACE ELEVATION CONTOURS'
      LINES=LA-1
      LEVELS=1
      DBS=0.
      WRITE (10,99) TITLE
      WRITE (10,100)LINES,LEVELS
      WRITE (10,250)DBS
      CLOSE(10)
      JSRPPH=0
  300 CONTINUE
      IF(ISDYNSTP.EQ.0)THEN
        TIME=DT*FLOAT(N)+TCON*TBEGIN
        TIME=TIME/TCON
      ELSE
        TIME=TIMESEC/TCON
      ENDIF
      OPEN(10,FILE='RSURFCN.OUT',POSITION='APPEND',STATUS='UNKNOWN')
      WRITE (10,100)N,TIME
      DO L=2,LA
        SURFEL=HLPF(L)+BELV(L)
        WRITE(10,200)IL(L),JL(L),DLON(L),DLAT(L),SURFEL
      ENDDO
      CLOSE(10)
   99 FORMAT(A80)
  100 FORMAT(I10,F12.4)
  101 FORMAT(2I10)
  200 FORMAT(2I5,1X,6E14.6)
  250 FORMAT(12E12.4)
      RETURN
      END

max_stars_repo_name: fearghalodonncha/DeepCurrent

max_issues_repo_name: fearghalodonncha/DeepCurrent

max_forks_repo_name: fearghalodonncha/DeepCurrent


text:
!!  Copyright (C)  Stichting Deltares, 2012-2016.
!!
!!  This program is free software: you can redistribute it and/or modify
!!  it under the terms of the GNU General Public License version 3,
!!  as published by the Free Software Foundation.
!!
!!  This program is distributed in the hope that it will be useful,
!!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
!!  GNU General Public License for more details.
!!
!!  You should have received a copy of the GNU General Public License
!!  along with this program. If not, see <http://www.gnu.org/licenses/>.
!!
!!  contact: delft3d.support@deltares.nl
!!  Stichting Deltares
!!  P.O. Box 177
!!  2600 MH Delft, The Netherlands
!!
!!  All indications and logos of, and references to registered trademarks
!!  of Stichting Deltares remain the property of Stichting Deltares. All
!!  rights reserved.

      subroutine varoxy ( pmsa   , fl     , ipoint , increm , noseg  ,
     &                    noflux , iexpnt , iknmrk , noq1   , noq2   ,
     &                    noq3   , noq4   )
!>\file
!>       Variation of oxygen due to variation in primary production within day

!
!     Description of the module :
!
! Name    T   L I/O   Description                                   Units
! ----    --- -  -    -------------------                            ----

!     Logical Units : -

!     Modules called : -

!     Name     Type   Library
!     ------   -----  ------------

      REAL     PMSA  ( * ) , FL    (*)
      INTEGER  IPOINT( * ) , INCREM(*) , NOSEG , NOFLUX,
     +         IEXPNT(4,*) , IKNMRK(*) , NOQ1, NOQ2, NOQ3, NOQ4

      INTEGER  IP1 , IP2 , IP3 , IP4 , IP5 , IP6 , IP7 , IP8 , IP9 ,
     J         IP10, IP11, IP12, I
      REAL     TIMSIM, DELTAT, TIMNUL, T1MXPP, T2MXPP, DAYLEN, FPPTOT,
     J         FRESPI, DEPTHW, T1    , T2    , PPMAX , TRISE ,
     J         TSET  , TOTAL , V1    , V2
      REAL     INTEGR(0:12*24), PPLAST, RELAST
      LOGICAL  FIRST
      SAVE     FIRST, PPLAST, RELAST, INTEGR
      DATA     FIRST /.TRUE./
      DATA     PPLAST, RELAST /-999.,-999./

      IP1  = IPOINT( 1)
      IP2  = IPOINT( 2)
      IP3  = IPOINT( 3)
      IP4  = IPOINT( 4)
      IP5  = IPOINT( 5)
      IP6  = IPOINT( 6)
      IP7  = IPOINT( 7)
      IP8  = IPOINT( 8)
      IP9  = IPOINT( 9)
      IP10 = IPOINT(10)
      IP11 = IPOINT(11)
      IP12 = IPOINT(12)

!     Check whether certain input parameters are independent of X

      IF (FIRST) THEN
          FIRST = .FALSE.
          IF ( (INCREM(1) .GT. 0) .OR.
     J         (INCREM(2) .GT. 0) .OR.
     J         (INCREM(3) .GT. 0) .OR.
     J         (INCREM(4) .GT. 0) .OR.
     J         (INCREM(5) .GT. 0) .OR.
     J         (INCREM(6) .GT. 0) .OR.
     J         (INCREM(7) .GT. 0) ) THEN

              WRITE (*,*)
     J        ' VAROXY: Time parameters function(x) not ALLOWED'
              CALL SRSTOP(1)
          ENDIF
      ENDIF

      IFLUX = 1
      DO 9000 ISEG = 1 , NOSEG

          TIMSIM = PMSA(IP1)/PMSA(IP3)
          DELTAT = PMSA(IP4)
          TIMNUL = PMSA(IP2)
          T1MXPP = PMSA(IP5)
          T2MXPP = PMSA(IP6)
          DAYLEN = PMSA(IP7)*24.
          FPPTOT = PMSA(IP8)
          FRESPI = PMSA(IP9)
          DEPTHW = PMSA(IP10)
          TRISE  = 12.0-0.5*DAYLEN
          TSET   = 12.0+0.5*DAYLEN

!         Initialize light variation curve for present cycle
!         ONLY if fluxes have changed

          IF ( ISEG .EQ. 1 ) THEN

              IF ( (ABS(FPPTOT-PPLAST) .GT. 1E-3) .OR.
     J             (ABS(FRESPI-RELAST) .GT. 1E-3) ) THEN

!                 WRITE (*,*) ' INITIALIZE ', TIMSIM
                  PPLAST = FPPTOT
                  RELAST = FRESPI

!                 Check on conditions for daylength

                  IF ( T1MXPP .LT. TRISE .OR.
     J                 T2MXPP .GT. TSET ) THEN
                      WRITE (*,*)
     J                ' VAROXY: Illegal definition of T1MXPP/T2MXPP'
                      CALL SRSTOP(1)
                  ENDIF

                  PPMAX = 48.0/(T2MXPP-T1MXPP+DAYLEN)
!                 PPMAX = 48.0*(FPPTOT+FRESPI)/(T2MXPP-T1MXPP+DAYLEN)

!                 Compute normalized integral Flux.dt in of (gC/m2/d)*h
!                 from t=0 to t=T every 5 minutes

                  TOTAL = 0.0
                  INTEGR(0) = 0.0
                  T1 = 0.0
                  V1 = 0.0
                  DO 100 I = 1,12*24
                      T2 = REAL(I)/12.
                      IF ( T2 .LE. TRISE .OR. T2 .GE. TSET ) THEN
                          V2 = 0.0
                      ELSEIF ( T2.GT.TRISE .AND. T2.LT.T1MXPP ) THEN
                          V2 = PPMAX*(T2-TRISE)/(T1MXPP-TRISE)
                      ELSEIF ( T2.GE.T1MXPP .AND. T2.LE.T2MXPP ) THEN
                          V2 = PPMAX
                      ELSEIF ( T2.GT.T2MXPP .AND. T2.LT.TSET ) THEN
                          V2 = PPMAX*(1.0-(T2-T2MXPP)/(TSET- T2MXPP) )
                      ENDIF
                      TOTAL = TOTAL + ((V1+V2)/2.0) * (T2-T1)
                      INTEGR(I) = TOTAL
                      V1 = V2
                      T1 = T2
  100             CONTINUE
              ENDIF
          ENDIF

!!        CALL DHKMRK(1,IKNMRK(ISEG),IKMRK1)
!!        IF (IKMRK1.EQ.1) THEN
          IF (BTEST(IKNMRK(ISEG),0)) THEN
!
!            Compute FLUX only if SWITCH is 1.0
!
             IF ( PMSA(IP11) .GT. 0.5 ) THEN

!               Compute relative time within day of time step to come

                T1 = (TIMSIM-INT(TIMSIM))*24.0 + TIMNUL
                IF ( T1 .GE. 24.0 ) T1 = T1 - 24.0
                T2 = T1 + DELTAT*24.0
                IF ( T2 .GT. 24.001 ) THEN
!                   WRITE (*,*) ' VAROXY: Illegal time frame'
!                   CALL SRSTOP(1)
                    T2 = 24.
                    T1 = T2 - DELTAT*24.0
                ENDIF
                PMSA(IP12) = T1

!               Compute flux for interval [T1:T2] by subtracting integrals
!               for both times and dividing by time interval

!               FL(IFLUX)  = (( INTEGR(NINT(T2*12.0))
!    J                         -INTEGR(NINT(T1*12.0)) )
!    J                       / (T2-T1)- FRESPI ) /DEPTHW
                FL(IFLUX)  = (( INTEGR(NINT(T2*12.0))
     J                         -INTEGR(NINT(T1*12.0)) )
     J                       / (T2-T1)* (FPPTOT+FRESPI)
     J                       - FRESPI ) / DEPTHW
!                IF (ISEG.LE.2)
!    J           WRITE (*,*) ISEG,FL(IFLUX),DEPTHW

             ELSE
                FL(IFLUX)  = 0.0
             ENDIF
          ENDIF
!
          IFLUX = IFLUX + NOFLUX
!
          IP1  = IP1  + INCREM( 1)
          IP2  = IP2  + INCREM( 2)
          IP3  = IP3  + INCREM( 3)
          IP4  = IP4  + INCREM( 4)
          IP5  = IP5  + INCREM( 5)
          IP6  = IP6  + INCREM( 6)
          IP7  = IP7  + INCREM( 7)
          IP8  = IP8  + INCREM( 8)
          IP9  = IP9  + INCREM( 9)
          IP10 = IP10 + INCREM(10)
          IP11 = IP11 + INCREM(11)
          IP12 = IP12 + INCREM(12)
!
 9000 CONTINUE
!
      RETURN
!
      END

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd