text:
!***************************************************************************************************************************
! This is the code from the section (DYNAMIC == RUNINIT) lines 1479 - 1822 of the original CSCAS code. The names of the
! dummy arguments are the same as in the original CSCAS code and the call statement and are declared here. The variables
! that are not arguments are declared in module YCA_First_Trans_m. Unless identified as by MF, all comments are those of
! the original CSCAS.FOR code.
!
! Subroutine YCA_RUNINIT sets parameters, gets command line arguments, creates file names and their logical unit
! numbers, writes the headers for the main ooutput files.
!***************************************************************************************************************************

    SUBROUTINE YCA_RunInit (&
        CN          , DOY         , FILEIOIN    , FROP        , IDETL       , ISWNIT      , ON          , RN          , &
        RNMODE      , RUN         , SN          , TN          , YEAR        &
        )

        USE ModuleDefs
        USE YCA_First_Trans_m

        IMPLICIT NONE

        INTEGER CN          , DOY         , FROP        , ON          , RN          , RUN         , SN          , TN
        INTEGER YEAR
        !INTEGER VERSIONCSCAS
        INTEGER TVILENT                                                                       ! Integer function call

        CHARACTER(LEN=1)   IDETL       , ISWNIT      , ISWWAT      , RNMODE
        CHARACTER(LEN=120) FILEIOIN
        CHARACTER(LEN=10)  TL10FROMI                                                          ! Character function call

        !-----------------------------------------------------------------------
        !         The body of YCA_RunInit
        !-----------------------------------------------------------------------

        IF (RUNCRP <= 0) THEN
            MODNAME(1:8) = 'CSYCA' // ModelVerTxt
            VERSIONCSCAS = 010114                                                             ! MF 15SE14 Changed from VERSION, conflict with ModuleDefs
            GENFLCHK(1:15) = 'CSYCA047.041621'

            !-----------------------------------------------------------------------
            !         Set parameters (Most should be placed in input files!)
            !-----------------------------------------------------------------------

            ! Physical constants
            MJPERE = 220.0*1.0E-3  ! MJ per Einstein at 540 nm
            PARMJFAC = 0.5         ! PAR in SRAD (fr)

            ! Standard parameters
            STDAY = 20.0      ! TT in standard day
            STEPNUM = 1       ! Step number per day set to 1

            ! Methods that not in control file
            ! (Impementation needs a re-compile)
            MERNU = 'CSM'     ! Root N uptake
            MEDEV = 'DEVU'    ! Control of development (Lnum or devu) LPM 04MAR15 change from LNUM to DEVU

            ! Control flags that not in control file
            ! (Impementation needs a  re-compile)
            ! (Only useful when examining how the model is working)
            ISWWATCROP = 'Y'    ! N=no water stress effect on crop growth
            CFLAFLF = 'Y'       ! Assimilate effect on leaf size
            CFLTFG = 'Y'        ! Temperature effect on leaf size

            !-----------------------------------------------------------------------
            !         Read command line arguments for model name and path (Cropsim)
            !-----------------------------------------------------------------------

            arg = ' '
            tvi2 = 0
            tvi3 = 0
            tvi4 = 0
            ! Following not good for all platforms.Changed for portability
            ! CALL GETARG (0,arg,arglen)
            CALL GETARG(0,arg)
            arglen = len_trim(arg)
            DO tvi1 = 1,arglen
                IF (arg(tvi1:tvi1) == Slash) tvi2=tvi1
                IF (arg(tvi1:tvi1) == '.') tvi3=tvi1
                IF (arg(tvi1:tvi1) == ' ' .AND. tvi4 == 0) tvi4=tvi1
            ENDDO
            IF (TVI3 == 0 .AND. TVI4 > 0) THEN
                tvi3 = tvi4
            ELSEIF (TVI3 == 0 .AND. TVI4 == 0) THEN
                tvi3 = arglen+1
            ENDIF
            MODEL = ARG(TVI2+1:TVI3-1)
            CALL UCASE(MODEL)

            !-----------------------------------------------------------------------
            !         Set configuration file name (Cropsim)
            !-----------------------------------------------------------------------

            IF (FILEIOT(1:2) /= 'DS') THEN
                CFGDFILE = ' '
                IF (TVI2 > 1) THEN
                    CFGDFILE = ARG(1:TVI2)//'CROPSIM.CFG'
                ELSE
                    CFGDFILE(1:12) = 'CROPSIM.CFG '
                ENDIF
            ENDIF

            !-----------------------------------------------------------------------
            !         Determine input file type (Dssat or X-file) and check if there
            !-----------------------------------------------------------------------

            TVI1 = TVILENT(FILEIOIN)
            IF (FILEIOIN(TVI1-2:TVI1) == 'INP') THEN
                FILEIOIN(TVI1:TVI1) = 'H'
                FILEIOT = 'DS4'
            ELSE
                FILEIOT = 'XFL'
            ENDIF
            FILEIO = ' '
            FILEIO(1:TVI1) = FILEIOIN(1:TVI1)
            INQUIRE (FILE = FILEIO,EXIST = FFLAG)
            IF (.NOT.(FFLAG)) THEN
                CALL GETLUN ('ERROR.OUT',FNUMERR)
                OPEN (UNIT = FNUMERR,FILE = 'ERROR.OUT')
                WRITE(fnumerr,*) ' '
                WRITE(fnumerr,*) 'Input file not found!     '
                WRITE(fnumerr,*) 'File sought was:          '
                WRITE(fnumerr,*) Fileio(1:78)
                WRITE(fnumerr,*) 'Please check'
                WRITE(*,*) ' Input file not found!     '
                WRITE(*,*) 'File sought was:          '
                WRITE(*,*) Fileio(1:78)
                WRITE(*,*) ' Program will have to stop'
                CLOSE (FNUMERR)
                STOP ' '
            ENDIF

            !-----------------------------------------------------------------------
            !         Create output file extensions (For different components)
            !-----------------------------------------------------------------------

            CNCHAR = ' '
            CNCHAR2 = '  '
            IF (CN == 1.OR.CN == 0) THEN
                OUT = 'OUT'
                CNCHAR2= '1 '
            ELSE
                CNCHAR = TL10FROMI(CN)
                OUT = 'OU'//CNCHAR(1:1)
                CNCHAR2(1:1) = CNCHAR(1:1)
            ENDIF

            !-----------------------------------------------------------------------
            !         Set output flags to agree with run modes and control switches
            !-----------------------------------------------------------------------

            IF (FILEIOT == 'XFL') THEN
                IF (RNMODE == 'I'.OR.RNMODE == 'E'.OR.RNMODE == 'A') THEN
                    IDETD = 'M'
                ELSEIF (RNMODE == 'B'.OR.RNMODE == 'N'.OR.RNMODE == 'Q' .OR. RNMODE .EQ. 'Y')THEN
                    IDETD = 'S'
                ENDIF
            ELSE
                IDETD = 'N'
            ENDIF
            FROPADJ = FROP
            IF (RNMODE == 'T') FROPADJ = 1
            IF (IDETL == 'D'.OR.IDETL == 'A') FROPADJ = 1

            !-----------------------------------------------------------------------
            !         Set file names and determine file unit numbers
            !-----------------------------------------------------------------------

            ! DATA FILES
            CALL GETLUN ('FILET',FNUMT)

            ! WORK,ERROR,AND TEMPORARY FILES
            CALL GETLUN ('ERROR.OUT',FNUMERR)
            CALL GETLUN ('FNAMETMP',FNUMTMP)

            ! IDETG FILES
            ! Check if need to change file names
            CALL XREADC (FILEIO,TN,RN,SN,ON,CN,'FNAME',fname)
            IF (FNAME == 'Y') THEN   ! File name change required.
                CALL XREADT (FILEIO,TN,RN,SN,ON,CN,'EXPER',excode)
                !NB. Renaming of Plantgro and Plantn handled by CSM
                OUTPG = 'PlantGro.'//OUT
                OUTPN = 'PlantN.'//OUT
                OUTPG2 = EXCODE(1:8)//'.OP2'
                OUTPG3 = 'FreshWt.'//OUT
                OUTPGF = EXCODE(1:8)//'.OPF'
            ELSE
                OUTPG = 'PlantGro.'//OUT
                OUTPN = 'PlantN.'//OUT
                OUTPG2 = 'PlantGr2.'//OUT
                OUTPG3 = 'FreshWt.'//OUT
                OUTPGF = 'PlantGrf.'//OUT
            ENDIF
            CALL GETLUN (OUTPG,NOUTPG)
            CALL GETLUN (OUTPG2,NOUTPG2)
            CALL GETLUN (OUTPG3,NOUTPF)
            CALL GETLUN (OUTPGF,NOUTPGF)
            CALL GETLUN (OUTPN,NOUTPN)

            ! IDETO FILES
            ! NB. Renaming of Overview and Evaluate handled by CSM
            FNAMEOV = 'Overview.'//out
            FNAMEEVAL = 'Evaluate.'//out
            FNAMEMEAS = 'Measured.'//out
            CALL GETLUN (FNAMEEVAL,fnumeval)
            CALL GETLUN (FNAMEOV,fnumov)
            CALL GETLUN (FNAMEMEAS,fnummeas)

            ! IDETS FILES
            FNAMEPSUM(1:12)   = 'Plantsum.'//OUT
            CALL GETLUN (FNAMEPSUM,  fnumpsum)

            ! RESPONSE FILES
            FNAMEPRES(1:12)   = 'Plantres.'//out
            FNAMEPREM(1:12) = 'Plantrem.'//out
            CALL GETLUN (FNAMEPRES,  fnumpres)
            CALL GETLUN (FNAMEPREM,fnumprem)

            ! LEAVES FILES
            FNAMELEAVES(1:10) = 'Leaves.'//OUT
            CALL GETLUN (FNAMELEAVES,fnumlvs)

            ! PHENOL FILES
            FNAMEPHASES(1:9) = 'Tiers.'//out
            FNAMEPHENOLS(1:11) = 'Phenols.'//out
            FNAMEPHENOLM(1:11) = 'Phenolm.'//out
            CALL GETLUN (FNAMEPHASES,fnumpha)
            CALL GETLUN (FNAMEPHENOLS,fnumphes)
            CALL GETLUN (FNAMEPHENOLM,fnumphem)

            ! ERROR FILES
            FNAMEERA(1:12) = 'Plantera.'//out
            FNAMEERT(1:12) = 'Plantert.'//out
            CALL GETLUN (FNAMEERT,fnumert)
            CALL GETLUN (FNAMEERA,fnumera)

            !-----------------------------------------------------------------------
            !         Open and write main headers to output files
            !-----------------------------------------------------------------------

            ! WARNING AND WORK FILES
! FO/LPM/GH/CHP - 12-04-2020 - READS.out file removed from CSM output.
!            IF (IDETL .EQ. 'Y') THEN
!              ! Close and re-open Work file
!              CALL Getlun('READS.OUT',fnumrea)
!              ! Close and re-open Reads file
!              CLOSE (FNUMREA, STATUS = 'DELETE')
!              OPEN (UNIT = FNUMREA,FILE = 'READS.OUT', STATUS = 'NEW', ACTION = 'READWRITE')
!              WRITE(fnumrea,*)' '
!              WRITE(fnumrea,*)' File closed and re-opened to avoid generating huge file'
!            ENDIF

            IF (RUN == 1) THEN
                ! IDETG FILES
                OPEN (UNIT = NOUTPG, FILE = OUTPG)
                WRITE (NOUTPG,'(A27)')'$GROWTH ASPECTS OUTPUT FILE'
                CLOSE (NOUTPG)
                OPEN (UNIT = NOUTPG2, FILE = OUTPG2)
                WRITE (NOUTPG2,'(A38)')'$GROWTH ASPECTS SECONDARY OUTPUTS FILE'
                CLOSE (NOUTPG2)
                OPEN (UNIT = NOUTPF, FILE = OUTPG3)
                WRITE(NOUTPF,'("*Fresh Weight Output File")')
                CLOSE (NOUTPF)
                OPEN (UNIT = NOUTPGF, FILE = OUTPGF)
                WRITE (NOUTPGF,'(A27)')'$GROWTH FACTOR OUTPUTS FILE'
                CLOSE (NOUTPGF)
                IF (ISWNIT /= 'N') THEN
                    OPEN (UNIT = NOUTPN, FILE = OUTPN)
                    WRITE (NOUTPN,'(A35)')'$PLANT NITROGEN ASPECTS OUTPUT FILE'
                    CLOSE (NOUTPN)
                ELSE
                    INQUIRE (FILE = OUTPN,EXIST = FEXIST)
                    IF (FEXIST) THEN
                        OPEN (UNIT = NOUTPN, FILE = OUTPN, STATUS='UNKNOWN',POSITION = 'APPEND')
                        CLOSE (UNIT=NOUTPN, STATUS = 'DELETE')
                    ENDIF
                ENDIF

                ! IDETO FILES
                OPEN (UNIT = FNUMOV, FILE = FNAMEOV)
                WRITE(FNUMOV,'(A20)') '$SIMULATION_OVERVIEW'
                CLOSE(FNUMOV)
                OPEN (UNIT = FNUMEVAL, FILE = FNAMEEVAL)
                WRITE(FNUMEVAL,'(A17)') '$PLANT_EVALUATION'
                CLOSE(FNUMEVAL)
                OPEN (UNIT = FNUMMEAS,FILE = FNAMEMEAS)
                WRITE (FNUMMEAS,'(A22)') '$TIME_COURSE(MEASURED)'
                CLOSE(FNUMMEAS)

                ! IDETS FILES
                OPEN (UNIT = FNUMPSUM,FILE = FNAMEPSUM)
                WRITE (FNUMPSUM,'(A27)') '$PLANT_SUMMARY             '
                CLOSE(FNUMPSUM)

                OPEN(UNIT=FNUMLVS,FILE=FNAMELEAVES)
                WRITE (FNUMLVS,'(A11)') '$LEAF_SIZES'
                CLOSE(FNUMLVS)

                OPEN(UNIT=FNUMPHA,FILE=FNAMEPHASES)
                WRITE (FNUMPHA,'(A16)') '$TIER_CONDITIONS'
                CLOSE(FNUMPHA)
                OPEN(UNIT=FNUMPHES,FILE=FNAMEPHENOLS)
                WRITE (FNUMPHES,'(A27)') '$PHENOLOGY_DATES(SIMULATED)'
                CLOSE(FNUMPHES)
                OPEN(UNIT=FNUMPHEM,FILE=FNAMEPHENOLM)
                WRITE (FNUMPHEM,'(A27)') '$PHENOLOGY_DATES(MEASURED) '
                CLOSE(FNUMPHEM)

                OPEN (UNIT = FNUMPRES,FILE = FNAMEPRES,STATUS = 'UNKNOWN')
                WRITE (FNUMPRES,'(A27)') '$PLANT_RESPONSES(SIMULATED)'
                CLOSE(FNUMPRES)
                OPEN (UNIT = FNUMPREM,FILE = FNAMEPREM,STATUS = 'UNKNOWN')
                WRITE (FNUMPREM,'(A26)') '$PLANT_RESPONSES(MEASURED)'
                CLOSE(FNUMPREM)

                ! ERROR FILES
                INQUIRE (FILE = FNAMEERA,EXIST = FFLAG)
                OPEN (UNIT = FNUMERA,FILE = FNAMEERA,STATUS = 'UNKNOWN')
                WRITE (FNUMERA,'(A27)') '$ERRORS   As % of measured '
                CLOSE(FNUMERA)
                OPEN (UNIT = FNUMERT,FILE = FNAMEERT,STATUS = 'UNKNOWN')
                WRITE (FNUMERT,'(A20)') '$ERRORS(TIME_COURSE)'
                WRITE (FNUMERT,*)' '
                WRITE (FNUMERT,'(A25)')'! Errors as % of measured'
                CLOSE(FNUMERT)

                ! Initialize 'previous' variables
                CROPPREV = '  '
                VARNOPREV = ' '
                CUDIRFLPREV = ' '
                ECONOPREV = ' '
                ECDIRFLPREV = ' '
                SPDIRFLPREV = ' '
            ENDIF

        ENDIF ! End of first time through stuff

        !-----------------------------------------------------------------------
        !       Record/set starting information
        !-----------------------------------------------------------------------

        YEARSIM = YEAR*1000 + DOY
        SEASENDOUT = 'N'  ! Season end outputs flag

        !-----------------------------------------------------------------------
        !       Increment counter for re-initialized runs within module
        !-----------------------------------------------------------------------

        RUNCRP = RUNCRP + 1

    END SUBROUTINE YCA_RunInit

max_stars_repo_name: Hafmic/dssat-csm-os

max_issues_repo_name: Hafmic/dssat-csm-os

max_forks_repo_name: Hafmic/dssat-csm-os


text:
! difference operator, cell to node
module diff_cn_op
implicit none
contains

subroutine diff_cn(df, f, i, a, i1, i2, x, dx)
real, intent(out) :: df(:,:,:)
real, intent(in) :: f(:,:,:,:), x(:,:,:,:), dx(3)
integer, intent(in) :: i, a, i1(3), i2(3)
real :: h, b
integer :: j, k, l, b, c

if (any(i1 > i2)) return

h = sign(1.0 / 12.0, product(dx))
b = modulo(a, 3) + 1
c = modulo(a + 1, 3) + 1
!$omp parallel do schedule(static) private(j, k, l)
do l = i1(3), i2(3)
do k = i1(2), i2(2)
do j = i1(1), i2(1)
df(j,k,l) = h * f(j-1,k-1,l-1,i) * &
( (x(j-1,k,l,b) - x(j,k-1,l-1,b)) * (x(j,k,l-1,c) - x(j,k-1,l,c)) + x(j-1,k,l,b) * (x(j-1,k,l-1,c) - x(j-1,k-1,l,c)) &
+ (x(j,k-1,l,b) - x(j-1,k,l-1,b)) * (x(j-1,k,l,c) - x(j,k,l-1,c)) + x(j,k-1,l,b) * (x(j-1,k-1,l,c) - x(j,k-1,l-1,c)) &
+ (x(j,k,l-1,b) - x(j-1,k-1,l,b)) * (x(j,k-1,l,c) - x(j-1,k,l,c)) + x(j,k,l-1,b) * (x(j,k-1,l-1,c) - x(j-1,k,l-1,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h  * f(j,k-1,l-1,i) * &
( (x(j+1,k,l,b) - x(j,k-1,l-1,b)) * (x(j,k-1,l,c) - x(j,k,l-1,c)) + x(j+1,k,l,b) * (x(j+1,k-1,l,c) - x(j+1,k,l-1,c)) &
+ (x(j,k-1,l,b) - x(j+1,k,l-1,b)) * (x(j,k,l-1,c) - x(j+1,k,l,c)) + x(j,k-1,l,b) * (x(j,k-1,l-1,c) - x(j+1,k-1,l,c)) &
+ (x(j,k,l-1,b) - x(j+1,k-1,l,b)) * (x(j+1,k,l,c) - x(j,k-1,l,c)) + x(j,k,l-1,b) * (x(j+1,k,l-1,c) - x(j,k-1,l-1,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h * f(j-1,k-1,l,i) * &
( (x(j,k,l+1,b) - x(j-1,k-1,l,b)) * (x(j-1,k,l,c) - x(j,k-1,l,c)) + x(j,k,l+1,b) * (x(j-1,k,l+1,c) - x(j,k-1,l+1,c)) &
+ (x(j-1,k,l,b) - x(j,k-1,l+1,b)) * (x(j,k-1,l,c) - x(j,k,l+1,c)) + x(j-1,k,l,b) * (x(j-1,k-1,l,c) - x(j-1,k,l+1,c)) &
+ (x(j,k-1,l,b) - x(j-1,k,l+1,b)) * (x(j,k,l+1,c) - x(j-1,k,l,c)) + x(j,k-1,l,b) * (x(j,k-1,l+1,c) - x(j-1,k-1,l,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h * f(j,k-1,l,i) * &
( (x(j,k-1,l,b) - x(j+1,k,l+1,b)) * (x(j+1,k,l,c) - x(j,k,l+1,c)) + x(j,k-1,l,b) * (x(j+1,k-1,l,c) - x(j,k-1,l+1,c)) &
+ (x(j,k,l+1,b) - x(j+1,k-1,l,b)) * (x(j,k-1,l,c) - x(j+1,k,l,c)) + x(j,k,l+1,b) * (x(j,k-1,l+1,c) - x(j+1,k,l+1,c)) &
+ (x(j+1,k,l,b) - x(j,k-1,l+1,b)) * (x(j,k,l+1,c) - x(j,k-1,l,c)) + x(j+1,k,l,b) * (x(j+1,k,l+1,c) - x(j+1,k-1,l,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h * f(j-1,k,l-1,i) * &
( (x(j,k+1,l,b) - x(j-1,k,l-1,b)) * (x(j,k,l-1,c) - x(j-1,k,l,c)) + x(j,k+1,l,b) * (x(j,k+1,l-1,c) - x(j-1,k+1,l,c)) &
+ (x(j,k,l-1,b) - x(j-1,k+1,l,b)) * (x(j-1,k,l,c) - x(j,k+1,l,c)) + x(j,k,l-1,b) * (x(j-1,k,l-1,c) - x(j,k+1,l-1,c)) &
+ (x(j-1,k,l,b) - x(j,k+1,l-1,b)) * (x(j,k+1,l,c) - x(j,k,l-1,c)) + x(j-1,k,l,b) * (x(j-1,k+1,l,c) - x(j-1,k,l-1,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h * f(j,k,l-1,i) * &
( (x(j,k,l-1,b) - x(j+1,k+1,l,b)) * (x(j,k+1,l,c) - x(j+1,k,l,c)) + x(j,k,l-1,b) * (x(j,k+1,l-1,c) - x(j+1,k,l-1,c)) &
+ (x(j+1,k,l,b) - x(j,k+1,l-1,b)) * (x(j,k,l-1,c) - x(j,k+1,l,c)) + x(j+1,k,l,b) * (x(j+1,k,l-1,c) - x(j+1,k+1,l,c)) &
+ (x(j,k+1,l,b) - x(j+1,k,l-1,b)) * (x(j+1,k,l,c) - x(j,k,l-1,c)) + x(j,k+1,l,b) * (x(j+1,k+1,l,c) - x(j,k+1,l-1,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h * f(j-1,k,l,i) * &
( (x(j-1,k,l,b) - x(j,k+1,l+1,b)) * (x(j,k,l+1,c) - x(j,k+1,l,c)) + x(j-1,k,l,b) * (x(j-1,k,l+1,c) - x(j-1,k+1,l,c)) &
+ (x(j,k+1,l,b) - x(j-1,k,l+1,b)) * (x(j-1,k,l,c) - x(j,k,l+1,c)) + x(j,k+1,l,b) * (x(j-1,k+1,l,c) - x(j,k+1,l+1,c)) &
+ (x(j,k,l+1,b) - x(j-1,k+1,l,b)) * (x(j,k+1,l,c) - x(j-1,k,l,c)) + x(j,k,l+1,b) * (x(j,k+1,l+1,c) - x(j-1,k,l+1,c)) )
end do
do j = i1(1), i2(1)
df(j,k,l) = df(j,k,l) + h * f(j,k,l,i) * &
( (x(j+1,k,l,b) - x(j,k+1,l+1,b)) * (x(j,k+1,l,c) - x(j,k,l+1,c)) + x(j+1,k,l,b) * (x(j+1,k+1,l,c) - x(j+1,k,l+1,c)) &
+ (x(j,k+1,l,b) - x(j+1,k,l+1,b)) * (x(j,k,l+1,c) - x(j+1,k,l,c)) + x(j,k+1,l,b) * (x(j,k+1,l+1,c) - x(j+1,k+1,l,c)) &
+ (x(j,k,l+1,b) - x(j+1,k+1,l,b)) * (x(j+1,k,l,c) - x(j,k+1,l,c)) + x(j,k,l+1,b) * (x(j+1,k,l+1,c) - x(j,k+1,l+1,c)) )
end do
end do
end do
!$omp end parallel do

end subroutine

end module

max_stars_repo_name: gely/coseis

max_issues_repo_name: gely/coseis

max_forks_repo_name: gely/coseis


text:
	PROGRAM SFDELT
C************************************************************************
C* SFDELT								*
C*									*
C* This program deletes data from a surface data file.			*
C*									*
C* Log:									*
C* I. Graffman/RDS	12/85						*
C* M. desJardins/GSFC	10/86	Added GEMPAK parameter names		*
C* M. desJardins/GSFC	 6/88	Rewrote from SFDELTM			*
C* M. desJardins/GSFC	 1/89	Added declaration for variable TTT	*
C* L. Williams/EAI	 3/94	Clean up declarations of user input	*
C*				variables				*
C* S. Jacobs/NMC         4/94   Removed unused variables        	*
C* L. Williams/EAI	 7/94	Removed call to SFDUPD			*
C* K. Tyle/GSC	 	 8/96	Added FL_MFIL to search for file type	*
C* S. Maxwell/GSC        7/97   Increased input character length        *
C* A. Hardy/GSC		 3/99   Added priority parameter to SF_SNXT     *
C* A. Hardy/GSC		 3/99   Removed ispri = 0			*
C* T. Piper/SAIC	 4/02	Fixed UMR; initialized sffcur		*
C************************************************************************
	INCLUDE		'GEMPRM.PRM'
C*
	CHARACTER 	sffile*(LLMXLN), dattim*(LLMXLN), area*(LLMXLN)
C*
	LOGICAL		respnd, done, proces, allflg
	CHARACTER	timdst (LLMXTM)*20, times  (LLMXTM)*20
	CHARACTER	ttt*20, sffcur*72, filnam*72
	CHARACTER	stn*8
C-----------------------------------------------------------------------
	sffcur = ' '
C
C*	Initilaize user interface.
C
	CALL IP_INIT  ( respnd, iperr )
	IF  ( iperr .ne. 0 )  THEN
	    CALL ER_WMSG  ( 'SFDELT', iperr, ' ', ier )
	    CALL SS_EXIT
	END IF
	CALL IP_IDNT  ( 'SFDELT', ier )
C
C*	Main loop.
C
	done   = .false.
	isffln = 0
	DO WHILE (.not. done)
C
C*	    Get user input and exit if there is an error.
C
	    CALL SFDINP  ( sffile, dattim, area, iperr )
	    IF  ( iperr .ne. 0 )  THEN
		CALL ER_WMSG  ( 'SFDELT', iperr, ' ', ier )
		CALL SS_EXIT
	    END IF
	    proces = .true.
C
C*	    Open the surface data file.
C
	    CALL FL_MFIL ( sffile, ' ', filnam, iret )
	    IF ( iret .ne. 0 ) CALL ER_WMSG ( 'FL', iret, ' ', ier )
	    CALL SFDFIL  ( filnam, sffcur, isffln, ntdset, timdst, ier )
	    IF  ( ier .ne. 0 )  proces = .false.
C
C*	    Find the times to delete.
C
	    IF  ( proces )  THEN
		CALL TI_FIND  ( dattim, ntdset, timdst, ttt,
     +				ntime,  times,  ier )
		IF  ( ier .ne. 0 )  proces = .false.
	    END IF
C
C*	    Check area for DSET or ALL.  Otherwise, set area.
C
	    IF  ( proces )  THEN
		CALL ST_LCUC  ( area, area, ier )
		IF (( area .eq. 'DSET' ) .or. ( area .eq. 'ALL' )) THEN
		    allflg = .true.
		  ELSE
		    allflg = .false.
		    CALL LC_SARE  ( area, isffln, stn, ier )
		    IF  ( ier .ne. 0 )  proces = .false.
		END IF
	    END IF
C
C*	    Loop through all times.
C
	    itime = 1
	    DO WHILE  ( proces .and. ( itime .le. ntime ) )
C
C*		Give the user a chance to exit.
C
		CALL SFDOPT  ( times (itime), filnam, area,
     +			       ier )
C
C*		If user wants to exit, do not process any more times.
C
		IF  ( ier .ne. 0 )  THEN
		    proces = .false.
C
C*		    Otherwise, process time.
C
		  ELSE IF  ( allflg )  THEN
C
C*		    Delete all data for this time.
C
		    CALL SF_DTIM  ( isffln, times (itime), ier )
		    IF  ( ier .ne. 0 )  THEN
			CALL ER_WMSG ('SF', ier, times (itime), ierr)
			CALL ER_WMSG ('SFDELT', -3, times (itime), ier)
		    END IF
		  ELSE
C
C*		    Set the time.  Then loop through deleting data from
C*		    stations.
C
		    CALL SF_STIM  ( isffln, times (itime), ier )
		    ier = 0
		    DO WHILE  ( ier .eq. 0 )
			CALL SF_SNXT  ( isffln, stn, istnm, slat,
     +					slon, selv, ispri, ier )
			IF  ( ier .eq. 0 )  CALL SF_DDAT ( isffln, i )
		    END DO
		END IF
C
C*		Go to next time.
C
		itime = itime + 1
	    END DO
C
C*	    Call the dynamic tutor.
C
	    CALL IP_DYNM  ( done, ier )
	END DO
C
C*	Close any file that is open.
C
	CALL SF_CLOS  ( isffln, ier )
C
C*	Exit.
C
  	CALL IP_EXIT  ( iret )
	END

ext: f

max_stars_repo_name: oxelson/gempak

max_issues_repo_name: oxelson/gempak

max_forks_repo_name: oxelson/gempak


text:
      SUBROUTINE MP_HELAS_CALLS_AMPB_1(P,NHEL,H,IC)
C
      USE POLYNOMIAL_CONSTANTS
      IMPLICIT NONE
C
C     CONSTANTS
C
      INTEGER    NEXTERNAL
      PARAMETER (NEXTERNAL=3)
      INTEGER    NCOMB
      PARAMETER (NCOMB=12)

      INTEGER NBORNAMPS
      PARAMETER (NBORNAMPS=1)
      INTEGER    NLOOPS, NLOOPGROUPS, NCTAMPS
      PARAMETER (NLOOPS=1, NLOOPGROUPS=1, NCTAMPS=1)
      INTEGER    NLOOPAMPS
      PARAMETER (NLOOPAMPS=2)
      INTEGER    NWAVEFUNCS,NLOOPWAVEFUNCS
      PARAMETER (NWAVEFUNCS=3,NLOOPWAVEFUNCS=3)
      REAL*16     ZERO
      PARAMETER (ZERO=0.0E0_16)
      COMPLEX*32     IZERO
      PARAMETER (IZERO=CMPLX(0.0E0_16,0.0E0_16,KIND=16))
C     These are constants related to the split orders
      INTEGER    NSO, NSQUAREDSO, NAMPSO
      PARAMETER (NSO=1, NSQUAREDSO=1, NAMPSO=2)
C
C     ARGUMENTS
C
      REAL*16 P(0:3,NEXTERNAL)
      INTEGER NHEL(NEXTERNAL), IC(NEXTERNAL)
      INTEGER H
C
C     LOCAL VARIABLES
C
      INTEGER I,J,K
      COMPLEX*32 COEFS(MAXLWFSIZE,0:VERTEXMAXCOEFS-1,MAXLWFSIZE)
C
C     GLOBAL VARIABLES
C
      INCLUDE 'mp_coupl_same_name.inc'

      INTEGER GOODHEL(NCOMB)
      LOGICAL GOODAMP(NSQUAREDSO,NLOOPGROUPS)
      COMMON/FILTERS/GOODAMP,GOODHEL

      INTEGER SQSO_TARGET
      COMMON/SOCHOICE/SQSO_TARGET

      LOGICAL UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE,CT_REQ_SO_DONE
     $ ,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE,MP_LOOP_REQ_SO_DONE
     $ ,CTCALL_REQ_SO_DONE,FILTER_SO
      COMMON/SO_REQS/UVCT_REQ_SO_DONE,MP_UVCT_REQ_SO_DONE
     $ ,CT_REQ_SO_DONE,MP_CT_REQ_SO_DONE,LOOP_REQ_SO_DONE
     $ ,MP_LOOP_REQ_SO_DONE,CTCALL_REQ_SO_DONE,FILTER_SO

      COMPLEX*32 AMP(NBORNAMPS)
      COMMON/MP_AMPS/AMP
      COMPLEX*32 W(20,NWAVEFUNCS)
      COMMON/MP_W/W

      COMPLEX*32 WL(MAXLWFSIZE,0:LOOPMAXCOEFS-1,MAXLWFSIZE
     $ ,0:NLOOPWAVEFUNCS)
      COMPLEX*32 PL(0:3,0:NLOOPWAVEFUNCS)
      COMMON/MP_WL/WL,PL

      COMPLEX*32 AMPL(3,NCTAMPS)
      COMMON/MP_AMPL/AMPL

C
C     ----------
C     BEGIN CODE
C     ----------

C     The target squared split order contribution is already reached
C      if true.
      IF (FILTER_SO.AND.MP_CT_REQ_SO_DONE) THEN
        GOTO 1001
      ENDIF

      CALL MP_OXXXXX(P(0,1),ZERO,NHEL(1),-1*IC(1),W(1,1))
      CALL MP_IXXXXX(P(0,2),ZERO,NHEL(2),+1*IC(2),W(1,2))
      CALL MP_VXXXXX(P(0,3),MDL_MW,NHEL(3),+1*IC(3),W(1,3))
C     Amplitude(s) for born diagram with ID 1
      CALL MP_FFV2_0(W(1,2),W(1,1),W(1,3),GC_47,AMP(1))
C     Counter-term amplitude(s) for loop diagram number 2
      CALL MP_FFV2_0(W(1,2),W(1,1),W(1,3),R2_SXCW,AMPL(1,1))
C     At this point, all CT amps needed for (QCD=2), i.e. of split
C      order ID=1, are computed.
      IF(FILTER_SO.AND.SQSO_TARGET.EQ.1) GOTO 2000

      GOTO 1001
 2000 CONTINUE
      MP_CT_REQ_SO_DONE=.TRUE.
 1001 CONTINUE
      END

max_stars_repo_name: madnklo/madnklo

max_issues_repo_name: madnklo/madnklo

max_forks_repo_name: madnklo/madnklo


text:
      subroutine cubpol (x1,x2,y1,y2,dy1,dy2,b)
C...subroutine     cubpol    01/09/73
      dimension b(4)
      h0=1/(x2-x1)
      h1=h0*h0
      dy=y2-y1
      s1=x2+2*x1
      s2=x1+2*x2
      b(1)=h1*(2*h0*x1*x2*(x1*y2-x2*y1)+y1*x2**2+y2*x1**2-x1*x2*(dy1*x2+
     1dy2*x1))
      b(2)=h1*(-6*h0*x1*x2*dy+dy1*x2*s1+dy2*x1*s2)
      b(3)=3*h0*h1*(dy*(x2+x1))-h1*(dy1*s2+dy2*s1)
      b(4)=h1*(dy2+dy1-2*h0*dy)
      return
      end

max_stars_repo_name: sergev/grafor

max_issues_repo_name: sergev/grafor

max_forks_repo_name: sergev/grafor


!@(#)  show basic line drawing, text and (if applicable) color.
!
!  As none of the projection routines have been called we
!  move and draw in the initial coordinate system -1.0 to 1.0.
!
program fsimple
!(LICENSE:PD)
use M_draw
integer,parameter :: BLACK = 0, GREEN = 2
character(len=50) :: device
character(len=80) :: fname
character(len=11) :: p
data p/'Hello world'/
   print*,'Enter output device:'
   read(*,'(a)') device
   print*,'Enter a font name:'
   read(*,'(a)') fname
!  set up device
   call vinit(device)
!  change font to the argument
   call font(fname)
!  set current color
   call color(BLACK)
!  clear screen to current color
   call clear
   call color(GREEN)
!  2 d move to start where we want drawstr to start
   call move2(-0.9, 0.9)
!  draw string in current color
   call drawstr('A Simple Example')
!  the next four lines draw the x
   call move2(0.0, 0.0)
   call draw2(0.76, 0.76)
   call move2(0.0, 0.76)
   call draw2(0.76, 0.0)

   call move2(0.0, 0.5)
   call drawstr('x done')
   call drawstr('next sentence')

   call move2(0.0, 0.1)
   do i = 1, 11
      call drawchar(p(i:i))
   enddo
!  the next five lines draw the square
   call move2(0.0, 0.0)
   call draw2(0.76, 0.0)
   call draw2(0.76, 0.76)
   call draw2(0.0, 0.76)
   call draw2(0.0, 0.0)
!  wait for some input
   idum=getkey()
!  set the screen back to its original state
   call vexit
end program fsimple


text:
      SUBROUTINE SGMTML(X,I,J,K)

!Purpose:      Form product of operators to generate the full group

      INCLUDE       '../INCLDS/COPYRIGT.FOR'

C       This program was developed for
C                    The Division of Chemistry
C                               of
C               The National Research Council of Canada
C                               by
C       Allen C. Larson, 14 Cerrado Loop, Santa Fe, NM 87505, USA

      REAL*4        X(5,4,25)
      INTEGER*4     I                   !Input matrix number
      INTEGER*4     J                   !Input matrix number
      INTEGER*4     K                   !Output matrix number

      DO L=1,4
        DO M=1,4
          X(L,M,K) = 0.0
          DO N=1,4
            X(L,M,K) = X(L,M,K)+X(L,N,J)*X(N,M,I)
          END DO
        END DO
      END DO

      X(1,4,K) = MOD(NINT((7.0+X(1,4,K))*12)/12.0,1.0)            !Force the translations to be in the cell
      X(2,4,K) = MOD(NINT((7.0+X(2,4,K))*12)/12.0,1.0)            !Also reset them to the value nearest to n/12
      X(3,4,K) = MOD(NINT((7.0+X(3,4,K))*12)/12.0,1.0)

      X(5,1,K) = 81*(2*X(1,1,K)+3*X(1,2,K)+4*X(1,3,K))            !Calculate a matrix flag
     1  +9*(2*X(2,1,K)+3*X(2,2,K)+4*X(2,3,K))
     1  +2*X(3,1,K)+3*X(3,2,K)+4*X(3,3,K)
      X(5,2,K) = 1728*X(1,4,K)+144*X(2,4,K)+12*X(3,4,K)            !Calculate the translation flag
      X(5,2,K) = NINT(X(5,2,K))                              !These should be whole numbers
      X(5,3,K) = IEOR(NINT(X(5,3,J)),NINT(X(5,3,I)))                  !Note the generator matrix number
      X(5,4,K) = 0.0

      RETURN
      END

ext: for

max_stars_repo_name: pedrobcst/Xerus

max_issues_repo_name: pedrobcst/Xerus

max_forks_repo_name: pedrobcst/Xerus


text:
      subroutine total_Qe(y, z, ix)
      use EC_Objects_manager

!!      integer, parameter :: nume       = 40000
	  common/bk08/kprint,nstep,ite,ilimit,newstf
	  common/wblock8/  abc(573,nume,4), his(573,nume,4)
!!!	  common /stressflag/ ElemFractCode(nume),ElemDecayCount(nume)
	  common/meshnum/ numnpo, numelto
	  common /t_Qe/ totqe, totqe_last, ele_area

	  dimension y(*), z(*), ix(4,*)

	  integer numelto, ElemFractCode, nstep, ele
	  real totqe, totqe_last, ele_area
	  real y3, z3, y1, z1, ds0

	  if(nstep==0) then            ! only for square element
	      y3=y(ix(3,1))
	      z3=z(ix(3,1))
	      y1=y(ix(1,1))
	      z1=z(ix(1,1))
	      ds0=sqrt((y3-y1)**2+(z3-z1)**2)
		  ele_area=(ds0/sqrt(2.0))**2
	  end if

	  totqe_last=totqe
	  totqe=0.0
	  do ele=1,numelto
		  if(EC_GetElemSplit(ele)==0) then !-- already decayed
!!!	      if(ElemFractCode(ele)<2) then
		      totqe=totqe+abc(3,ele,1)*ele_area
		  end if
      end do

	  write(959,10) nstep, totqe, totqe-totqe_last

   10 format(i8,2(e20.10))

      end

max_stars_repo_name: imohame/LabCode

max_issues_repo_name: imohame/LabCode

max_forks_repo_name: imohame/LabCode


!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!                                                                            !!!!!!!!
!!!!!!!!                  MAIN FILE WITH SIMULATION CODE                            !!!!!!!!
!!!!!!!!                                                                            !!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!============================ files with subroutines ======================================!
	INCLUDE 'subroutines/modules.f90'		!modules
	INCLUDE 'subroutines/allocation.f90'		!allocation of vector sizes

	INCLUDE 'config.f90'				!configuration - parameters
	INCLUDE 'initial_conditions.f90'		!values for initial conditions
	INCLUDE 'initial_conditions2.f90'		!values for initial conditions

	INCLUDE 'subroutines/conditions.f90'		!sequence of three cases
	INCLUDE 'subroutines/simulation.f90'		!main simulation
	INCLUDE 'subroutines/LIF.f90'			!leaky integrate-and-fire implementation
	INCLUDE 'subroutines/input_step.f90'		!inputs from step
!==========================================================================================!
	PROGRAM COMPLEMENTARY_ISTDP
	USE PARAMETERS
	USE VARIABLES
	USE TEMP
	IMPLICIT NONE
	INTEGER			::	k
	CHARACTER(len=1)	:: arg
!======================== getting value from user for sims ================================!
	CALL get_command_argument(1, arg)
	READ(arg,'(I1)')sims_type
!======================== initialising subroutines ========================================!
	CALL config()					!call config file with parameters
	CALL allocation()				!allocate vectors and matrices based on parameters
	CALL initial_conditions()			!implement initial conditions
!==========================================================================================!
!========================== three conditions for sims =====================================!
	DO k = 1,3
		CALL conditions(k)	!1=control case; 2=pop 1 OFF; 3=pop 2 OFF
	END DO
!==========================================================================================!
	END PROGRAM
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!


text:
subroutine rdheat(lunmd     ,lundia    ,error     ,nrrec       ,mdfrec    , &
                & runid     ,filtem    ,fmttem    ,ktemp       ,dt        , &
                & itstrt    ,itfinish  ,ivapop    ,solrad_read ,gdp       )
!----- GPL ---------------------------------------------------------------------
!
!  Copyright (C)  Stichting Deltares, 2011-2016.
!
!  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 3.
!
!  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, "Delft3D" and "Deltares"
!  are registered trademarks of Stichting Deltares, and remain the property of
!  Stichting Deltares. All rights reserved.
!
!-------------------------------------------------------------------------------
!  $Id: rdheat.f90 5717 2016-01-12 11:35:24Z mourits $
!  $HeadURL: https://svn.oss.deltares.nl/repos/delft3d/tags/6686/src/engines_gpl/flow2d3d/packages/io/src/preprocessor/rdheat.f90 $
!!--description-----------------------------------------------------------------
!
!    Function: - Reads the time dependent boundary data for the
!                heat model directly from the MD-file or indi-
!                rectly from the attribute file Filtmp.
!              - Tests the file or data consistency.
!              - Checks whether the file exists or the required
!                data is not empty.
!              - An essential assumption is that the data has to
!                be specified sequentially in time. This imply
!                that NT times Background temperature records
!                should exist in the file (NT unrestricted).
!              - Writes the data to an unformatted file
! Method used:
!
!!--pseudo code and references--------------------------------------------------
! NONE
!!--declarations----------------------------------------------------------------
    use precision
    use properties
    !
    use globaldata
    use string_module
    !
    implicit none
    !
    type(globdat),target :: gdp
    !
    ! The following list of pointer parameters is used to point inside the gdp structure
    !
    integer , pointer :: itis
!
! Global variables
!
    integer                                      :: itfinish    !  Description and declaration in inttim.igs
    integer                                      :: itstrt      !  Description and declaration in inttim.igs
    integer                        , intent(in)  :: ivapop      !  Description and declaration in heat.igs
    integer                        , intent(in)  :: ktemp       !  Description and declaration in tricom.igs
    integer                                      :: lundia      !  Description and declaration in inout.igs
    integer                                      :: lunmd       !  Description and declaration in inout.igs
    integer                                      :: nrrec       !!  Pointer to the record number in the MD-file
    logical                                      :: error       !!  Flag=TRUE if an error is encountered
    logical                        , intent(in)  :: solrad_read !  Description and declaration in heat.igs
    real(fp)                                     :: dt          !  Description and declaration in esm_alloc_real.f90
    character(*)                                 :: filtem      !!  File name for the time varying heat model file
    character(*)                                 :: mdfrec      !!  Standard rec. length in MD-file (300)
    character(*)                                 :: runid       !!  Run identification code for the current simulation (used to determine
                                                                !!  the names of the in- /output files used by the system)
    character(2)                                 :: fmttem      !!  File format for the time varying heat model file
!
! Local variables
!
    integer                      :: iocond   ! IO status for reading
    integer                      :: itold    ! Help var. to store last read time to test accending order
    integer                      :: ittdep   ! Help var. for the time read (now defined as multiples of DT, but in future it may take any value)
    integer                      :: l        ! Help var.
    integer                      :: lenc     ! Help var. (length of var. cvar to be looked for in the MD-file)
    integer                      :: lf
    integer                      :: lkw      ! Length (in characters) of keyword
    integer                      :: lrid     ! Length of character string runid
    integer                      :: lunout   ! Unit number for unformatted FLOW help file between TDATOM and TRISIM
    integer                      :: n        ! Help var.
    integer                      :: newlun
    integer                      :: nlook    ! Help var.: nr. of data to look for in the MD-file
    integer                      :: nrval    ! Number of values to read from file
    integer                      :: ntrec    ! Help. var to keep track of NRREC
    logical                      :: dtn
    logical                      :: ex       ! Flag to test if file exists
    logical                      :: rec1st   ! Flag set to TRUE if the record read is the first record
    real(fp)                     :: rdummy
    real(fp)      , dimension(5) :: rval     ! Help array (real) where the data, recently read from the MD-file, are stored temporarily
    character(11)                :: fmtdef   ! Default file format (usually = blank)
    character(11)                :: fmttmp   ! Help variable for file format
    character(12)                :: fildef   ! Default file name (usually = blank)
    character(256)               :: filout   ! Help variable for file name
    character(6)                 :: keyw     ! Name of record to look for in the MD-file (usually KEYWRD or RECNAM)
    character(300)               :: message
!
!! executable statements -------------------------------------------------------
!
    itis  => gdp%gdrdpara%itis
    !
    ! initialize global parameters
    !
    itold  = -1
    rec1st = .true.
    nlook  = 1
    fildef = ' '
    filout = ' '
    fmtdef = 'FRformatted'
    fmttmp = ' '
    !
    lunout = 8
    !
    ! Locate 'Filtem' record for time varying heat model data in separate input file containing time-varying data
    ! Keyword does not have to exist, temperature might be specified space-varying via the meteo-module
    !
    filtem = fildef
    call prop_get_string(gdp%mdfile_ptr, '*', 'Filtmp', filtem)
    !
    ! Time varying heat module data in file?
    !
    if (filtem /= fildef) then
       !
       ! It is assumed the file is free-formatted
       !
       fmttem = 'FR'
       fmttmp = 'formatted'
       error  = .false.
       !
       ! Check filename "filtem" <> TMP file or
       ! "filtem" = TMP file and access is unformatted
       ! Define length of RUNID
       ! open output file (ONLY VERSION 2.48 or lower) +
       ! Set name for the constituents, Regenerated locally
       !
       call remove_leading_spaces(runid     ,lrid      )
       filout = 'TMP_' // runid(:lrid) // '.tem'
       !
       ! define length of file name
       !
       call remove_leading_spaces(filtem, lf)
       !
       ! open unformatted tem-file
       !
       lunout = newlun(gdp)
       inquire (file = filout(:8 + lrid), exist = ex)
       if (ex) then
          open (lunout, file = filout(:8 + lrid))
          close (lunout, status = 'delete')
       endif
       open (lunout, file = filout(:8 + lrid), form = 'unformatted',      &
           & status = 'unknown')
       !
       write (message, '(2a)') 'Reading Heat module file ', filtem(:lf)
       call prterr(lundia, 'G051', trim(message))
       nrval = 1
       if (ktemp /= 3) then
          nrval = 3
          if (ivapop == 1 .or. solrad_read) then
             nrval = 4
          endif

       endif
       call rdtdf(lundia    ,lunout    ,error     ,filtem    ,fmttmp    , &
                & nrval     ,rval      ,dt        ,itstrt    ,itfinish  , &
                & gdp       )
    else
       !
       ! time varying heat module data in a .tem file
       !
       ! Remove file TMP_<runid>.tem if it exists
       !
       call remove_leading_spaces(runid, lrid)
       filout = 'TMP_' // runid(:lrid) // '.tem'
       !
       ! Open file
       !
       lunout = newlun(gdp)
       inquire (file = filout(:8 + lrid), exist = ex)
       if (ex) then
          open (lunout, file = filout(:8 + lrid))
          close (lunout, status = 'delete')
       endif
       !
       ! Check whether heat parameters are handled by the (new) meteo module
       !
       filtem = fildef
       ex     = .false.
                     call prop_get_string(gdp%mdfile_ptr, '*', 'Fwndgr', filtem)
       if (filtem /= fildef) ex = .true.
       if (.not. ex) call prop_get_string(gdp%mdfile_ptr, '*', 'Fwndgt', filtem)
       if (filtem /= fildef) ex = .true.
       if (.not. ex) call prop_get_string(gdp%mdfile_ptr, '*', 'Fwndgc', filtem)
       if (filtem /= fildef) ex = .true.
       if (.not. ex) call prop_get_string(gdp%mdfile_ptr, '*', 'Filwr' , filtem)
       if (filtem /= fildef) ex = .true.
       if (.not. ex) call prop_get_string(gdp%mdfile_ptr, '*', 'Filwt' , filtem)
       if (filtem /= fildef) ex = .true.
       if (.not. ex) call prop_get_string(gdp%mdfile_ptr, '*', 'Filwc' , filtem)
       if (ex) then
          !
          ! Reading/checking heat data is done by the meteo module,
          ! when starting the calculation
          ! Skip this subroutine without error
          !
           error  = .false.
           lunout = 8
           goto 9999
       endif
       !
       ! -->
       !
  110  continue
       !
       ! If an NaN is read -> error
       !
       do n = 1, nrval
           if (isnan(rval(n))) then
              write(message,'(2a)') 'NaN in ', trim(filtem)
              call prterr(lundia    ,'P004'    ,message      )
              !
              error = .true.
              goto 9999
           endif
       enddo
       !
       ittdep = nint(rval(1)/dt)
       if (dtn(ittdep, rval(1), dt)) then
          error = .true.
          call prterr(lundia    ,'U044'    ,'Tstmp'   )
       !
       endif
       !
       ! test times and define minimum time
       !
       if (rec1st) then
          if (ittdep > itstrt) then
             call prterr(lundia    ,'U041'    ,'First time Tstmp  >' )
             error = .true.
          endif
          rec1st = .false.
       endif
       !
       if (ittdep <= itold) then
          call prterr(lundia    ,'U060'    ,'Tstmp'   )
          error = .true.
       endif
       !
       ! writing to LUNOUT
       !
       write (lunout) (rval(l), l = 1, nlook)
       !
       ! next time to read
       !
       itold = ittdep
       !
       ! <--
       !
       goto 110
       !
       ! stop reading
       !
       !
       ! define actual number of times for time varying heat data
       ! and define maximum time
       !
  500  continue
       if (itold/= - 1) then
          if (itold < itfinish) then
             write(message,'(3a)') 'Last time in file ', trim(filtem), ' <'
             call prterr(lundia    ,'U042'    ,message)
             error = .true.
             goto 9999
          endif
       endif
    endif
    !
    ! close files
    !
 9999 continue
    if (filtem == fildef) then
       error = .true.
       write (message, '(a)') 'No heat module file(s) found'
       call prterr(lundia, 'P004', message)
    endif
    if (lunout /= 8) then
       if (error) then
          close (lunout, status = 'delete')
       else
          close (lunout)
       endif
    endif
end subroutine rdheat

max_stars_repo_name: liujiamingustc/phd

max_issues_repo_name: liujiamingustc/phd

max_forks_repo_name: liujiamingustc/phd


text:
C--------------------------------------------------------------------
C Called by MONGO to decode the symbols as in NEWPLOT :
C Input : ISYMB, SIZE
C Output : NSIDES, ISTYLE, ANGLE, EXPAND
C ISTYLE : 0=OPEN, 1=SKELETAL (VERTICES CONNECTED TO THE CENTRE)
C	   2=STARRED, 3=FILLED IN
C--------------------------------------------------------------------
	SUBROUTINE MONGO_SYMBOL(ISYMB,SIZE,NSIDES,ISTYLE,ANGLE,EXPAND)

	EXPAND=SIZE/0.15
	ANGLE=0.

C Small dot :
	IF(ISYMB.EQ.1)THEN
	NSIDES=0
	ISTYLE=0
	ENDIF

C White triangle :
	IF(ISYMB.EQ.2)THEN
	NSIDES=3
	ISTYLE=0
	ENDIF

C Black triangle :
	IF(ISYMB.EQ.3)THEN
	NSIDES=3
	ISTYLE=3
	ENDIF

C Cross + :
	IF(ISYMB.EQ.4)THEN
	NSIDES=4
	ISTYLE=1
	ANGLE=45.
	ENDIF

C Cross X :
	IF(ISYMB.EQ.5)THEN
	NSIDES=4
	ISTYLE=1
	ENDIF

C White square :
	IF(ISYMB.EQ.6)THEN
	NSIDES=4
	ISTYLE=0
	ENDIF

C Black square :
	IF(ISYMB.EQ.7)THEN
	NSIDES=4
	ISTYLE=3
	ENDIF

C White circle :
	IF(ISYMB.EQ.8)THEN
	NSIDES=20
	ISTYLE=0
	ENDIF

C Black circle :
	IF(ISYMB.EQ.9)THEN
	NSIDES=20
	ISTYLE=3
	ENDIF

	RETURN
	END