text:
>Image(toplogo.jpg, noborder)<

Founded in 1987, Pilipino Culture Night (PCN) is a production hosted by Mga Kapatid that showcases singing, acting, and traditional Pilipino dances. It is a yearlong process that builds community between the seven Pilipino American organizations, as well as those who are excited to learn about and share our culture. Dedication, an enthusiastic spirit, and openness to new experiences are the characteristics needed to be a part of this landmark event.

PCN welcomes and encourages people of all backgrounds to come and learn about the Pilipino culture.

MISSION STATEMENT
We aspire to promote Pilipino culture and awareness by educating the participants and the attending audience of Pilipino Culture Night, as well as build a stronger, more unified family among the Pilipino and Pilipino American community by producing a great show.

The 20132014 school year will hold the 27th Annual PCN. This years PCN Coordinators are Joseph Fernandez, Ruth Genota, and Kevin Gonzales.

Previous PCN Themes
 20122013  Walang Kapalit
 20112012  Suntok sa Buwan
 20102011  Ang Tanging Pagmamahal / The Greatest Love of All
 20092010  Mga Kababayan Ko
 20082009  Ang Mamatay ng Dahil Sayo
 20072008  A Snapshot of Paradise, Wish You Were Here
 20062007  Kayumanging Anino: Brown on Both Sides
 20052006  Wansaponatym: Simply Epic
 20042005  Finding History: Revealing Echoes
 20032004  Where I Belong, Here All Along
 20022003  Tsismis / Gossip
 20012002  One of a Kind
 20002001  Wake the Sun
 19992000  Bridge to Glory
 19981999  Bound by Pride
 19971998  Ang Puso Ng Pilipinas / Heart of the Philippines
 19961997  Harana / Serenade
 19951996  Kung Hindi Kita Nakilala / If I Never Knew You
 19941995  Sunding Ang Puso Mo / Fold Your Heart
 19931994  Tayong Lahat, Mga Kapatid / Brothers and Sisters
 19921993  Mga Bayani Sa Aming Piling / Heroes Among Us
 19911992  Ang Daan Sa Bagong Mundo / A Road to the New World
 19901991  Ganito Tayo Ngayon / This Is Who We Are Now
 19891990  Ay Pagibig / My Love
 19881989  Ako Ay Pilipino / I Am Filipino
 19871988  Sino Ba Ako? / Who Am I?

How to Get Involved
Check the PCN website. http://www.mgakapatid.org/pcn
Find us on Facebook: https://www.facebook.com/groups/180400765353986/

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
include "sparseMod.f03"
      program sparse01
      USE sparseMod
!
!     This program will serve as a sandbox for learning about sparse vector and
!     sparse matrix algebra.
!
!
      implicit none
      integer::i,j
      integer::nDim,nDimSparse
      real,dimension(:),allocatable::vec1
      integer,dimension(:),allocatable::indexVec1
      real,dimension(:),allocatable::denseVec1
!
!     Let's get some values from the user...
!
      write(*,*)' What is the dimension of the vector space?'
      read(*,*) nDim
      write(*,*)' How many non-zero elements are there in the sparse vector?'
      read(*,*) nDimSparse
      ALLOCATE(vec1(nDimSparse),indexVec1(nDimSparse))
      write(*,*)'Please enter the non-zero vector elements and their values in pairs.'
      do i = 1,nDimSparse
        read(*,*) indexVec1(i),vec1(i)
      endDo
!
!     Echo the user input...
!
      write(*,*)' nDim = ',nDim
      write(*,*)' nDimSparse = ',nDimSparse
      do i = 1,nDimSparse
        write(*,*) 'i=',indexVec1(i),'  vec1(i) = ',vec1(i)
      endDo
!
!     Allocate space for a full-memory version of the sparse vector and then use
!     a scatter function.
!
      ALLOCATE(denseVec1(nDim))
      denseVec1 = scatter_realVector(vec1,indexVec1,nDim)
      do i = 1,nDim
        write(*,*) i,denseVec1(i)
      endDo
!
!     Scan over the values of the dense array and build a new index before we
!     run gather and rebuild the sparse array.
!
      j = 0
      do i = 1,nDim
        if (abs(denseVec1(i)).gt.0.001) then
          j = j+1
          if(j.gt.nDimSparse) STOP
          indexVec1(j) = i
        endIf
      endDo
      if(j.ne.nDimSparse) then
        write(*,*)' j = ',j,'  BUT nDimSparse = ',nDimSparse
        write(*,*)' YIKES!!!'
        STOP
      endIf
!
!     Test out the gather function.
!
      vec1 = 999.123123
      vec1 = gather_realVector(indexVec1,denseVec1)
      write(*,*)
      write(*,*)' Here''s the sparse vector after running the gather function...'
      do i = 1,nDimSparse
        write(*,*) 'i=',indexVec1(i),'  vec1(i) = ',vec1(i)
      endDo
!
      end program sparse01

max_stars_repo_name: MQCPack/summerCodingWorkshop

max_issues_repo_name: MQCPack/summerCodingWorkshop

max_forks_repo_name: MQCPack/summerCodingWorkshop


text:
Jump to Timeline #Navigation Navigation

  DresbachHuntBoyer House William Dresbach essentially founds Davis by leasing the Yolo House Jerome C. Davis house and maintaining surrounding buildings in the area of the anticipated railroad junction.

ext: f


text:
cdt this version had code added to calculate rain infiltration as
c constant head infiltration for a period. I have to finish it to allow final
c flux to just finish up the water.
CDT  9/10/2014 Added CO2 as a weather variable and took it out of the
CDT   initials file.

C**  MSW1-  switch to indicate if daily wet bulb temperatures are     **
C**         available (=1 if YES).                                    **
C**                                                                   **
C**  MSW2-  switch to indicate if hourly wind is available (=1 if YES).**
C**  MSW3-  Daily rain intensities                                    **
C** cdt don't think these two variables (above) are needed  for hourly data
C**  MSW4-  Chemical concentrations in rain water                     **
C**  MSW5-  1 for flood irrigation                                    **
C**  MSW6-  1 if relative humidity is available                       **
C**  MSW7-  1 if Daily CO2 Measurements are available                 **
C**                                                                   **
C**  BSOLAR-factor for changing solar radiation units. Equals         **
C**         radiation in J m-2 divided by radiation in units used in  **
C**         weather data file.                                        **
C**                                                                   **
C**  BTEMP- factor for changing temperature units. Equals change in   **
C**         temperature, in units used in weather data file, that is  **
C**         equivalent to a one oC change.                            **
C**                                                                   **
C**  ATEMP- factor for changing temperature units. Equals temperature **
C**         used in weather data file that is equivalent to zero oC.  **
C**                                                                   **
C**  ERAIN- factor for changing rain units. Equals rainfall in mm     **
C**         day-1 divided by rainfall in units used in weather data   **
C**         file.                                                     **
C**                                                                   **
C**  BWIND- factor for changing wind units. Equals windspeed in       **
C**         km hr-1 divided by windspeed in units used in weather     **
C**         data file.                                                **
C**                                                                   **
C**                                                                   **
C**                                                                   **
C**  WIND-  if windspeed is unavailable read in average value for     **
C**         site. windspeed at 2 meters (km hr-1).                    **
C**                                                                   **
C** CO2   -  If daily values are not available read average           **
C**                                                                   **
C**                                                                   **
C**  (File 1, Line 3 to EOF, free format)                             **
C**      (WIND, TWET and TDRY are optional)                           **
C**                                                                   **
C**  JDAY-  day of year- Julian Date.  Only used for management       **
C**         of the weather file. Program converts calendar date       **
C**         to a day of year relative to 3/1/1900                     **
C**                                                                   **
C**                                                                   **
C**  RI-    daily solar radiation integral (J m-2).                   **
C**                                                                   **
C**  TMAX-  maximum air temperature during the day (oC).              **
C**                                                                   **
C**  TMIN-  minimum air temperature during the day (oC).              **
C**                                                                   **
C**  RAIN-  rainfall (or irrigation) (mm day-1).                      **
C**                                                                   **
C**  WIND-  windspeed at 2 meters (km hr-1).                          **
C**                                                                   **
C**  TWET-  wet bulb temperature (oC).                                **
C**                                                                   **
C**  TDRY- dry bulb temperature corresponding to TWET (oC).           **
C**                                                                   **
C**  Hourly values for radiation will be converted to watts m-2 as
C**   soon as they are read in.
C**                                                                   **

c inputs hourly data
      Subroutine SetSurfaceH()
      Include 'Public.ins'
      Include 'Puplant.ins'
      Include 'Puweath.ins'
      Parameter (PERIOD =1./24.)
      integer jday,m,ThisYear,DayOfYear,CurYear
      double precision St,t
      real Interval
      character*10 date
      integer HOUR
      Common /weather/ il,im,HRAIN(24),HSR(24),HTEMP(24),HTEMPY(24),
     &     HWIND(24), Rel_Humid(24),AVP(24),isol,Date1,ModNum,
     &     Interval

      Dimension CLIMAT(20),SDERP(9),SINALT(24),SINAZI(24),HRANG(24),
     &           DIFWAT(24),SARANG(24),DIFINT(24),ROWINC(24),DIRINT(24),
     &           SHADOW(24),SOLALT(24),SOLAZI(24),
     &           HTM(24)
      Dimension xS(NumBPD),iS(NumBPD),kS(NumBPD)
      Data SDERP/0.3964E-0,0.3631E1,0.3838E-1,0.7659E-1,0.0000E0,
     & -0.2297E2,-0.3885E0,-0.1587E-0,-0.1021E-1/,PI /3.1415926/,
     &  TPI /6.2831852/,DEGRAD/0.017453293/,IPERD/24/

CYAP
       If(NSurf.eq.0) then
           RowSP=1.
          height=0.
        endif
CYAPEND
      t=Time
      St=0.1D0*Step
      If (lInput.eq.0) goto 11
C
C  First and last days
C
      JDFRST=idint(t+St)
      JDLAST=idint(tFin+St)
c
        NumMod=NumMod+1
        ModNum=NumMod
        tNext(ModNum)=time
        Interval=1/24.

c
      im=160
      il=0
C      Open (5,file='Weatherhrly.dat',status='old',ERR=10)
      Open (5,file=ClimateFile,status='old',ERR=10)
cdt
c       Do i=1,NumNP
c        PcodeW(i)=CodeW(i)
c         Enddo
C
C  Read descriptors and conversion constants for weather data
C
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10) LATUDE
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10) MSW1,MSW2,MSW3,MSW4,MSW5,MSW6,MSW7
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
       ISOL=NumSol*Movers(2)
       Read (5,*,ERR=10) BSOLAR,BTEMP,ATEMP,ERAIN,BWIND,BIR
cdt 10/25/2007 changed from 1-MSW4
      ISOL=NumSol*(MSW4)*Movers(2)
cdt  9/10/2014 added CO2 as a weather variable. Made this line consistent
Cdt   with the same line from the daily weather file
CDT    TODO this needs to be tested fully
      NCD=4-MSW2-MSW3+ISOL+(NumG+1)*Movers(4)-MSW7
C      NCD=1-MSW2+ISOL+(NumG+1)*Movers(4)
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read(5,*,ERR=10) (CLIMAT(i),i=1,NCD)
      IF(MSW2.eq.0) WINDA=CLIMAT(1)
C     IF(MSW4.eq.0.AND.NumSol.gt.0) then
Cdt 02/03/2009 changed from 0-MSW4
      IF(MSW4.eq.1.AND.NumSol.gt.0) then
        Do i=1,NumSol
          CPREC(i)=CLIMAT(2-MSW2+i)
        Enddo
      Endif
      If(NumG.ne.0) then
        PG=CLIMAT(3-MSW2+ISOL)
        Do i=1,NumG
          GAIR(i)=CLIMAT(3-MSW2-MSW3+ISOL+i)
        Enddo
      Endif
      IF(MSW7.eq.0) then
         CO2=CLIMAT(NCD)
       Endif

      close(5)
C
C Total number of weather data
C
      NCD=3+2*MSW1+MSW2+ISOL*MSW4+MSW6+MSW7
C
C     Nodal numbers of furrow nodes
C
      im1=im
      il1=il
      im=180
      il=0
      If(MSW5.gt.0) then
        Open(6,file='Furnod.dat',status='old',ERR=10)
        im=im+1
        il=il+1
        Read(6,*,ERR=10)
        im=im+1
        il=il+1
        Read(6,*,ERR=10)
        im=im+1
        il=il+1
        Read(6,*,ERR=10) NumFP
        im=im+1
        il=il+1
        Read(6,*,ERR=10)
        im=im+1
        il=il+1
        Read(6,*,ERR=10) (NumF(i),i=1,NumFP)
        im=im+1
        il=il+1
        Read(6,*,ERR=10)
        im=im+1
        il=il+1
        Read(6,*,ERR=10) (hFur(i),i=1,NumFP)
        IFUR=0
      Endif
      im=im1
      il=il1
C
C FIND CORRECT DAY IN WEATHER FILE
C
      Open (5,file=WeatherFile,status='old',ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
      im=im+1
      il=il+1
      Read (5,*,ERR=10)
1111  il=il+1
      Read (5,*,ERR=10) iDum, date, iDum
      MDAY=julday(date)
C       MDAY = 38065
cdt changed from LE to LT since jday was always gt. time
      If (MDAY.LT.JDFRST) GO TO 1111
      backspace(5)

C
C   CALCULATE A CLOUD COVER FACTOR FOR THIS LATITUDE
C
      If (LATUDE.LE.25.0) THEN
          CLDFAC = 0.45 - (LATUDE*0.004)
      Else
          CLDFAC = 0.30 + (LATUDE*0.002)
      Endif
c
c..................... Routine calculations

c..................... Input daily data
   11 If((linput.eq.1).OR.idint(t+St).eq.(JDAY+1)) then

      tAtm=t          !insure we call the weather code the next step

c       Do i=1,NumNP
c        CodeW(i)=PcodeW(i)
c         Enddo

c
c ZEROING DAILY ARRAYS
c
        Do i=1,24
          RINT(i)=0.
          RADINT(i)=0.
          WATTSM(i)=0.
          il=il-1
        Enddo
c
c Read meteodata
c
        Do m=1,24
          il=il+1
C          Read(5,*,ERR=10) JDAY,HOUR, DATE1, (CLIMAT(i),i=1,NCD)
          Read(5,*,ERR=10) JDAY, date, HOUR,  (CLIMAT(i),i=1,NCD)
           JDAY=julday(date)
c    SInce some routines need the day of the calendar year we need a variable to hold this
C    since the julian day is referenced to a time longer in the past
           ThisYear=CurYear(date)
           write (date,'("01/01/",i4.4)') ThisYear
           DayOfYear=JDAY-julday(date)

           HSR(M)=climat(1)*BSOLAR/3600   ! convert to watts m-2 hard
                                          !coded for hourly now

           if (lInput.ne.1) HTEMPY(M)=HTEMP(M)  ! save to yesterday's
                                                !temperature
           HTEMP(M)=(climat(2)-ATEMP)/BTEMP     ! convert to celcius
           HRAIN(M) = CLIMAT(3)*erain         ! convert to cm rain in
                                              ! this hour
           If (MSW2.gt.0) HWIND(M)=CLIMAT(4)
           If(MSW6.gt.0) then
             Rel_Humid(m)=Min(Climat(3+2*MSW1+MSW2+ISOL+MSW6)
     #              /100.,0.98)
            endif
         enddo

       If(MSW7.gt.0) then
          CO2=Climat(4+2*MSW1+MSW2+MSW3+ISOL*MSW4+MSW6+MSW7)
       EndIf

C
C      ADJUST UNITS FOR HOURLY DATA
c      HSR(M) is converted to an hours worth of Watts m-2 (Joules m-2 h-1)
C
202          RI = 0.0
          SUM24 = 0.0
          TMAX = 0.0
          TMIN = 1000.0
          DO 1015 M = 1,24
c RI should be total J m-2 -- work in Watts * time = total energy
             RI = RI + amax1(0.0,HSR(M))*3600
             SUM24 = SUM24 + HRAIN(M)
             IF(TMAX .LT. HTEMP(M)) THEN
               TMAX = HTEMP(M)
             ENDIF
             IF(TMIN .GT. HTEMP(M)) THEN
               TMIN = HTEMP(M)
             ENDIF
             IF(MSW2 .GT. 0) THEN
               WIND = HWIND(M) * BWIND
             Else
               WIND=WINDA*BWIND
             ENDIF
1015      CONTINUE
          RAIN = SUM24
        If (MSW1.GT.0) then
          TWET = (CLIMAT(5+MSW2) - ATEMP)/BTEMP
          TDRY = (CLIMAT(6+MSW2) - ATEMP)/BTEMP
        Endif
        If (MSW4.GT.0.AND.NumSol.ne.0) then
          Do i=1,NumSol
            CPREC(i)=CLIMAT(4+MSW2+2*MSW1+MSW3+i)
          Enddo
        Endif

c......................... Radiation submodel
C
C  CALCULATE SOLAR DECLINATION
C
        XLAT = LATUDE*DEGRAD
        DEC = SDERP(1)
        Do I = 2,5
          N = I - 1
          J = I + 4
          D11 = N*0.01721*JDAY
          DEC = DEC + SDERP(I)*SIN(D11) + SDERP(J)*COS(D11)
        Enddo
        DEC = DEC*DEGRAD
C
C  CALCULATE SOME PARTS OF THE SOLAR RADIATION EQUATIONS THAT
C         OCCUR REPEATEDLY LATER
C
        D12 = SIN(XLAT)*SIN(DEC)
        D13 = COS(XLAT)*COS(DEC)
        D14 = D12 + D13
C
C  CALCULATE DAYLENGTH
C
        DAYLNG = ACOS((-0.014544 - D12)/D13)*7.6394
C       7.6394 = 180/3.1416/360*24*2
C
C  CALCULATE SOLAR RADIATION INCIDENT ON TOP OF THE EARTH'S
C         ATMOSPHERE AT SOLAR NOON
C
        RADVEC = 1 + (0.01674*SIN((DayOfYear - 93.5)*0.9863*DEGRAD))
        WATATM = 1325.4*D14/(RADVEC*RADVEC)
C
C  CALCULATE AN ATMOSPHERIC TRANSMISSION COEFFICIENT FOR THIS
C         LATITUDE AND TIME OF YEAR. ***
C
        If (DayOfYear.LT.145) then
          ATRANS = 0.68 + ((5.25E-5*LATUDE) - 0.1E-3)*(145 - DayOfYear)
        Else If (DayOfYear.LE.237) then
          ATRANS = 0.68
        Else
          If (LATUDE.LE.30.0) then
            D15 = (LATUDE*5.5E-5) - 0.1E-3
          Else
            D15 = 0.65E-3 + (LATUDE*3.0E-5)
          Endif
          ATRANS = 0.68 + D15*(DayOfYear - 237)
        Endif
C
C  CALCULATE POTENTIAL DIRECT + DIfFUSE RADIATION INCIDENT ON
C         CROP AT SOLAR NOON
C
        WATPOT = WATATM*0.5*(0.93 - (0.02/D14) + ATRANS**(1/D14))

C
C  CALCULATE ACTUAL RADIATION INCIDENT ON CROP AT SOLAR NOON
C  GIVEN DAILY INTEGRAL AND ASSUMING RADIATION FLUX DENSITY
C  VARIES AS A HALF SINE WAVE OVER THE PHOTOPERIOD
C
        WATACT = RI*4.363E-4/DAYLNG     ! Watts per m2
C       4.363E-4 = 3.1416/3600/2
C
C  CALCULATE CLOUD COVER. ***
C
        WATRAT = WATACT/WATPOT
        If (WATRAT.GE.1.0) then
          CLOUD = 0.0
        Else
          CLOUD =(CLDFAC - SQRT((CLDFAC*CLDFAC) + (1.52*(1 - WATRAT))))
     &   /(-0.76)
C        1.52 = 4*0.38
C        0.76 = 2*0.38
          If (CLOUD.GE.1.0) CLOUD = 1.0
        Endif
C
C       DUSK AND DAWN
C
        DAWN = 12.0 - (DAYLNG/2.0)
        DUSK = 12.0 + (DAYLNG/2.0)
C
C  CALCULATE TOTAL RADIATION INCIDENT ON THE CROP AT EACH TIME
C
C  DIVIDE PHOTOPERIOD INTO IPERD EQUAL INCREMENTS AND CALCULATE
C  HOUR ANGLES BETWEEN THE MIDPOINTS AND SOLAR NOON
C
        IHPERD = IPERD/2
        Do 20, J = 1,IHPERD
          HRANG(J) = PI/12.0*(12 - J + 0.5)
 20     HRANG(IPERD - J + 1) = HRANG(J)
        DDIf = DAYLNG - (2.0*IfIX(DAYLNG/2))
        IUP = 13 - IfIX(DAYLNG/2)
        IDN = 12 + (13 - IUP)
        If (.NOT.((DDIf.GT.0.0).OR.(DDIf.LT.0.0)))then
          IDAWN = IUP
          IDUSK = IDN
        Else
          IDAWN = IUP - 1
          IDUSK = IDN + 1
        ENDIf
C
C  CALCULATE TOTAL RADIATION INCIDENT ON THE CROP AT EACH TIME
C
        Do I = 1,IPERD
          TIMH = I - 0.5
          WATTSM(I) = 0.0
          If(I.GE.IUP.AND.I.LE.IDN)
cdbg check the units here
     &      WATTSM(I) = HSR(I)
        Enddo
        If ((DDIf.GT.0.0).OR.(DDIf.LT.0.0)) then
           WATTSM(IDAWN) = HSR(IDAWN)
           WATTSM(IDUSK) = HSR(IDUSK)
        ENDIf

c..................... Temperature and vapour pressure submodel

C
C  ASSIGN AIR TEMPERATURE AT EACH TIME
C
C FIND THE TIME AFTER DAWN IN HOURS WHEN MAXIMUM
C  TEMPERATURE IS REACHED
C
        tMax=-1e6
        Do 70, m = 1,IPERD
           TAIR(M) = (HTEMP(M)-ATEMP)/BTEMP
           if(m.eq.IDUSK) TDUSK=TAIR(m)
           if (TAIR(M).GT.TMAX) then
             TMAX=TAIR(m)
             TMAXHR=m
            End If
 70   CONTINUE

C
C  If DAILY WET BULB TEMPERATURES ARE NOT AVAILABLE, ASSUME
C  MINIMUM TEMPERATURE = DEW POINT TEMPERATURE AND CALCULATE
C  WATER VAPOR PRESSURE.  ASSUME A VALUE FOR THE PSYCHROMETRIC
C  "CONSTANT"
C
cdt added msw6 to if statement here

        If (MSW1.NE.1.and.MSW6.eq.0) then
          Do i=1,iperd
            AVP(i) = 0.61*EXP((17.27*TMIN)/(TMIN + 237.3))
           enddo
          GAMMA = 0.0645
        Else
cdt
           if(MSW1.eq.1) then
C
C  SINCE DAILY WET BULB TEMPERATURES ARE AVAILABLE CALCULATE
C  SATURATION VAPOR PRESSURE AT THE WET BULB TEMPERATURES
C    http://www.pmel.org/HandBook/HBpage21.htm see this site for more info
C    on using Dewpoint to calculate RH
C
          SVPW = 0.61*EXP((17.27*TWET)/(TWET + 237.3))
C
C  CALCULATE THE HUMIDITY RATIO, D31, LATENT HEAT OF EVAPORATION,
C  D32 AND THE PSYCHROMETRIC "CONSTANT"
C
          D31 = 0.622*(SVPW/(101.3 - SVPW))
          D32 = 2500.8 - (2.37*TWET)
          GAMMA = 0.62*(1.006 + (1.846*D31))
     &   /((0.622 + D31)*(0.622 + D31)*D32)*101.3
C
C  CALCULATE ACTUAL WATER VAPOR PRESSURE
C
          Do i=1,iperd
             AVP(i) = SVPW - (GAMMA*(TDRY - TWET))
          enddo
        Endif

CDT
            if(MSW6.eq.1) then
             Do i=1,iperd
                TMean=(TMax+TMin)/2.0
                SVPW= 0.61*EXP((17.27*HTEMP(i))/(HTEMP(i) + 237.3))
                D31 = 0.622*(SVPW/(101.3 - SVPW))
                D32 = 2500.8 - (2.37*HTEMP(i))
                GAMMA = 0.62*(1.006 + (1.846*D31))
     &          /((0.622 + D31)*(0.622 + D31)*D32)*101.3
                AVP(i)= Rel_Humid(i)*SVPW
              enddo
          endif

        Endif
CDT
C
C  CALCULATE WATER STAURATION VAPOR PRESSURE AT AIR TEMPERATURE
C  AND THE SLOPE OF THE RELATIONSHIP AT THAT POINT
C
        Do 80, I = 1,IPERD
          SVPA = 0.61*EXP((17.27*TAIR(I))/(TAIR(I) + 237.3))
          DEL(I) = (0.61*EXP((17.27*(TAIR(I) + 1.))
     &   /(TAIR(I) + 1.0 + 237.3))) - SVPA
C
C  CALCULATE WATER VAPOR PRESSURE DEFICIT
C
          VPD(I) = SVPA - AVP(I)
          If (VPD(I).LE.0.0) VPD(I) = 0.0
 80     Continue

c....................Light interception submodel

        If(NShoot.ne.0) then
C
C  CALCULATE SOLAR ALTITUDE AND SOLAR AZIMUTH FOR EACH TIME
C
          Do 30 I = 1,IPERD
            SINALT(I) = D12 + (D13*COS(HRANG(I)))
            SOLALT(I) = ASIN(SINALT(I))
            SINAZI(I) = -COS(DEC)*SIN(HRANG(I))/COS(SOLALT(I))
 30       Continue
          D18 = 0.0
          Do 60, I = 1,IHPERD
            SOLAZI(I) = PI + ASIN(SINAZI(I))
            If (.NOT.((D18.GT.0.0).OR.(D18.LT.0.0))) then
              D19 = HRANG(I) - 0.01
              TESAZI = PI + ((ASIN(-COS(DEC)*SIN(D19))
     &       /(COS(ASIN(D12 + (D13*COS(D19)))))))
              If (TESAZI.LT.SOLAZI(I)) then
                SOLAZI(I) = PI - SOLAZI(I)
              Else
                D18 = 1.0
              Endif
            Endif
            SOLAZI(IPERD - I + 1) = TPI - SOLAZI(I)
 60       Continue
C
C  CALCULATE PROPORTION OF TOTAL RADIATION THAT IS DIFFUSE AT
C  EACH TIME
C
          Do 40 I = 1,IPERD
            D16 = SINALT(I)
            If (D16.LT.0.01) then
              DIFWAT(I) = 1.0
            Else
              D17 = (0.93 - (0.02/D16) - ATRANS**(1/D16))
     &       /(0.93 - (0.02/D16) + ATRANS**(1/D16))
              If (CLOUD.GT.0.0) then
                DIFWAT(I) = 1.0 - ((1.0 - D17)*(1.0 - CLOUD)
     &         /(1.0 - ((CLDFAC+(0.38*CLOUD))*CLOUD)))
              Else
                DIFWAT(I) = D17
              Endif
              If (DIFWAT(I).LT.0.0) DIFWAT(I) = 0.0
              If (DIFWAT(I).GT.1.0) DIFWAT(I) = 1.0
            Endif
C
C     *** CALCULATE PHOTOSYNTHETICALLY ACTIVE RADIATION (PAR)
C         INCIDENT ON THE CROP ASSUMING THAT
C         PAR = 0.43 DIRECT + 0.57 DIFFUSE RADIATION. ***
C
           PAR(I) = WATTSM(I)
     &          *((DIFWAT(I)*0.57) + ((1.0 - DIFWAT(I))*0.43))
40         CONTINUE
C
C
C  If CLOUD COVER IS COMPLETE, OMIT CALCULATIONS PERTAINING
C  TO DIRECT RADIATION INTERCEPTION
C

          If (CLOUD.GE.1.0) then
            Do I = 1,IPERD
              DIRINT(I) = 0.0
            Enddo
          Endif
C
C  CALCULATE ANGLE BETWEEN ROW AND SOLAR AZIMUTH FOR EACH TIME
C
          Do I = 1,IPERD
            SARANG(I) = ABS(SOLAZI(I) - (ROWANG*DEGRAD))
            If (SARANG(I).GT.3.*PI/2.) SARANG(I)
     &       = ABS(SARANG(I) - TPI)
            If (SARANG(I).GT.PI/2.) SARANG(I)
     &       = ABS(SARANG(I) - PI)
C
C  CALCULATE PROPORTION OF DIRECT RADIATION INTERCEPTED BY ROWS
C  OF PLANTS ASSUMING THEY ARE OPAQUE CYLINDERS
C
            SHADOW(I) = HEIGHT/SIN(ATAN(TAN(SOLALT(I))
     &     /SIN(SARANG(I))))
            DIRINT(I) = SHADOW(I)/ROWSP
            If (DIRINT(I).GT.1.0) DIRINT(I) = 1.0
            If (DIRINT(I).LT.0.0) DIRINT(I) = 0.0
          Enddo

C
C  CALCULATE PROPORTION OF DIfFUSE RADIATION INTERCEPTED BY ROWS
C  OF PLANTS ASSUMING THEY ARE OPAQUE CYLINDERS
C  FIRST DIVIDE ROWSPACING INTO 20 EQUAL PARTS AND FIND THE
C  MIDPOINTS FROM ROW TO MID ROW
C
          DIfFIN = 0.0
          Do I = 1,10
            ROWINC(I) = ROWSP/20.0*(I - 0.5)
C
C  CALCULATE PROPORTION OF SKY OBSCURED BY "OPAQUE" ROWS AT EACH
C  POSITION FROM ROW TO MID-ROW
C
            DIfINT(I) = (ATAN(HEIGHT/2.0/(ROWSP - ROWINC(I)))
     &     +ATAN(HEIGHT/2.0/ROWINC(I)))/1.5708
            If (DIfINT(I).GT.1.0) DIfINT(I) = 1.0
C
C  INTEGRATE DIfFUSE RADIATION INTERCEPTION ACROSS THE ROW. ***
C
            DIfFIN = DIfFIN + DIfINT(I)
          Enddo
c
          DIfFIN = DIfFIN/10.0
C
C  CALCULATE LEAF AREA PER UNIT GROUND AREA COVERED BY CROP
C  CANOPY
C
C

C  CALCULATE AN EFFECTIVE LCAI ALLOWING FOR THE FACT THAT LIGHT
C  AT A LOW ANGLE TRAVERSES MORE LEAF LAYERS TO REACH THE SOIL
C

           LCAI=LAI

          Do 50, I = 1,IPERD
            If (SOLALT(I).GT.0.0) then
              ELCAI = LCAI/(SIN(ACOS(COS(SARANG(I))
     &       *COS(SOLALT(I)))))
              If (SHADOW(I).GT.ROWSP)
     &           ELCAI = ELCAI*SHADOW(I)/ROWSP
C
C  CALCULATE BEERS LAW CORRECTION FOR RADIATION INTERCEPTION
C
              If((ELCAI*CEC).LE.88.0) then
                BEERS = (1 - EXP(-ELCAI*CEC))
              Else
                BEERS=1.0
              Endif

C  CALCULATE PROPORTION OF PHOTOSYNTHETICALLY ACTIVE RADIATION
C   (PAR) INTERCEPTED BY THE CROP CANOPY ASSUMING THAT
C    PAR = 0.43 DIRECT + 0.57 DIFFUSE RADIATION. ***

      PARINT(I) = WATTSM(I)*((DIFFIN*0.57*DIFWAT(I))
     &     + (DIRINT(I)*0.43*(1.0 - DIFWAT(I))))
     &     /PAR(I)*BEERS
C
C  CALCULATE PROPORTION OF TOTAL RADIATION INTERCEPTED BY THE
C  CROP CANOPY
C
              RADINT(I) = ((DIfFIN*DIfWAT(I)) + (DIRINT(I)
     &         *(1.0 - DIfWAT(I))))*BEERS
            Else If (WATTSM(I).GT.0.0) then
              BEERS = (1 - EXP(-LCAI*CEC))
              PARINT(I) = COVER*BEERS
              RADINT(I) = COVER*BEERS
            Else
               PARINT(I) = 0.0
               RADINT(I) = 0.0
            Endif
 50       Continue

C
c..................End of radiation interception submodel

        Endif  !nshoot >0
C
c..................Daily calculations for the precipitation
C
          Do i=1,24
            RINT(i)=HRain(i)*24.0  ! scale hourly rainfall to a daily rate
          Enddo

        If(Rain.lt.0.) then
          QF=QF-Rain
          Do i=1,NumFP
            k=NumF(i)
            CodeW(k)=1
            hNew(k)=hFur(i)
          Enddo
          IfUR=1
        Endif
c........................End of daily calculations
      Endif
C
C     Further we have hourly calculations

C
55    If(Abs(Time-tNext(ModNum)).lt.0.01*Step.or.lInput.eq.1) then
       ITIME=Idint(dMOD(t+St,1.D0)/PERIOD+1)
       if(itime.eq.10) then
        iii=1;
        endif

       Do i=1, NumBP
         do j=1, 3
           VarBW(i,j)=0.0
         Enddo
        Enddo

c       Do i=1,NumNP
c        CodeW(i)=PCodeW(i)
c         Enddo

C
C  FOR DAYLIGHT PERIODS:
C  CALCULATE NET UPWARD LONG-WAVE RADIATION ***
C
1000    Continue
        If (WATTSM(ITIME).GT.0.0) then
          RNLU = (1.6E-3*WATRAT*(100.0 - TAIR(ITIME)))*697.6
        Else
          RNLU=0.
        ENDIf
C
C  NO SHADE IF NO PLANT
C
        If(NShoot.eq.0) then
          SHADE = 0.
          COVER = 0.
        ENDIf
C
C  CALCULATE POTENTIAL WATER EVAPORATION RATE FROM THE SOIL
C  SURFACE
C
C  THE TOTAL RADIATION FALLING ON THE SOIL IS ASSUMED TO BE
C  CONCENTRATED ON THE EXPOSED SOIL CELLS AND AN EQUIVALENT
C  TOTAL RADIATION IS CALCULATED
C
        If (Cover.lt.0.97) then
C
          D11 = WATTSM(ITIME)

        Else
          D11=0.
        Endif
c
c  MARCH ALONG THE SURFACE TO FIND SOIL-ATMOSPHERE BOUNDARY NODES
c  AND THEIR LATERAL COORDINATES
c
        ic=0
        Do k=1,NumBP
          i=KXB(k)
          If(abs(CodeW(i)).eq.4) then
            ic=ic+1
            xS(ic)=x(i)
            iS(ic)=i
            kS(ic)=k
          Endif
        Enddo
        Call SORT03(ic,xS,iS,kS)
C
        Do 90 i=1,ic
          PSh=FSh(i,xBSTEM,SHADE,xS,ic)
c
c  IF THE NODE IS EXPOSED THEN
c
          If(PSh.gt.0.) then
C
C  CALCULATE THE ALBEDO OF THE EXPOSED SOIL CELLS
C
cdt changed 0.5 to 0.1
              LAMDAS = 0.25 - 0.05*ThNew(iS(i))
C
C  CALCULATE NET RADIATION ON THE EXPOSED SOIL CELLS
C
              RNS = ((1.0 - LAMDAS)*D11) - RNLU
              If(RNS.LE.0.0) RNS = 0.0
C
C  CALCULATE POTENTIAL EVAPORATION RATE FROM EXPOSED SOIL CELLS
C
              D12 = max(1.0,2.0 - (2.0*COVER))
              If (D12.GE.1.0) D12 = 1.0
              ESO = ((DEL(ITIME)/GAMMA*RNS*3600.0/(2500.8
     &     - (2.3668*TAIR(ITIME))))
     &     + (VPD(ITIME)*109.375*(1.0 + (0.149*WIND*D12))))
     &      /((DEL(ITIME)/GAMMA) + 1.0)
C
C   PROPORTIONAL TO THE AREA EXPOSED
C

              VarBW(kS(i),2)=PSh*24.*ESO/10000.
c
c  ELSE IF THE NODE IS COVERED
c
          Else
            VarBW(kS(i),2)=0.
          Endif
c                Endif for the surface nodes
 90     Continue
C
C........................ End evaporation
C
C........................ Transpiration
C
        If (NShoot.ne.0) then
C
C  CALCULATE A CROP SURFACE ROUGHNESS PARAMETER
C
        If (ITIME.EQ.1) then
           ROUGH = 1.0/(ABS((HEIGHT/ROWSP) - 0.5) + 0.5)
           If (ROUGH.LT.1.0) ROUGH = 1.0
        Endif
C
C   ON STILL HOT DAYS, CONVECTION CURRENTS AT CROP LEVEL ASSURE
C   SOME AIR MOVEMENT
C
        WINDL = WIND
        If (TAIR(ITIME).GT.25.0.AND.WINDL.LE.0.36) WINDL = 0.36
C
C   CALCULATE POTENTIAL TRANSPIRATION RATE FOR THE CROP
C
C   ASSUME A VALUE FOR THE ALBEDO OF THE CROP
C
        LAMDAC = 0.23
C
C   THE TOTAL RADIATION INTERCEPTED BY THE CROP IS ASSUMED TO BE
C   SPREAD UNIfORMLY OVER THE AREA COVERED BY THE CROP AND AN
C   EQUIVALENT TOTAL RADIATION IS CALCULATED
C
C   CODE TO ACCOUNT FOR THE CASE WHERE THE PLANT
C   HAS NOT YET EMERGED AND COVER =0
        If (COVER.le.0.) then
          WATPL=0.0
        Else
          WATPL = WATTSM(ITIME) !*RADINT(ITIME)/COVER /dt don't need when used in 2dsoil
        Endif
C
C   CALCULATE NET RADIATION ON THE CROP
C
        RNC = ((1.0 - LAMDAC)*WATPL) - RNLU
        If (RNC.LE.0.0) RNC = 0.0
C
C   CALCULATE POTENTIAL TRANSPIRATION RATE FOR CROP ALLOWING FOR
C   INCOMPLETE GROUND COVER
C
        EPO = ((DEL(ITIME)/GAMMA*RNC*3600.0/(2500.8
     &   - (2.3668*TAIR(ITIME))))
     &   + (VPD(ITIME)*109.375*(ROUGH +(0.149*WINDL))))
     &   /((DEL(ITIME)/GAMMA) + 1.0)
C
C   CODE ADDED TO PREVENT DIVISION BY o
C
        If (COVER.le.0.0.OR.EPO.LE.0.0) EPO=1.0E-4
C
C     *** CALCULATE POTENTIAL TRANSPIRATION RATE FOR AN AVERAGE PLANT **
C
cdt I took out etcorr, the ET correction factor - units are per plant
        EO = 24.*EPO*cover/(poprow/rowsp*100)
        If(COVER.le.0.0) then
          TPot=0.
        Else
          TPot=EO/COVER
        Endif
      Endif
c................... End Transpiration .................................

c................... Precipitation and irrigation

      Do i=1,NumBP
      n=KXB(i)
      k=CodeW(n)
      If(K.eq.4.or.K.eq.-4) then
         VarBW(i,1)=RINT(ITIME)
       If(NumSol.ne.0) then
          do j=1,NumSol
cdt was varbs(n,j)
          VarBS(i,j)=CPREC(j)
          Enddo
       Endif
*
         VarBW(i,3)=VarBW(i,2)-VarBW(i,1)
         Q(n)=-Width(i)*VarBW(i,3)
         If (Q(n).gt.0.0) then
          CodeW(n)=-4
          endif
       Endif
       Enddo
c................... End of the precipitation and irrigation
C ............... Heat Balance is calculated here
C
C
CYAP
C FCSH is sensible heat flux when the ground is hotter than the air
       if(rowsp.le.0) then
         FCSH=4.0E-3+1.39E-3*Wind
        else
         FCSH=4.0E-3+1.39E-3*Wind*max(1.0,(2.-(2.*height/rowsp)))
       endif

c      FCSH=0.02
cdt
* FELWR is reradiation when soil is warmer than the air.
      FCSH=FCSH*4.1856*60.0*24.0
      FELWR=6.8E-03+9.0E-05*Tair(Itime)
      FELWR=FELWR*4.1856*60.0*24.0
* units are Joules cm-2 d-1 C-1
CYAPEND
C units are changed from cal cm-2 min-1 C-1 -->J cm-2 d-1 C-1
      cWat=4.1856
      TairN=Tair(Itime)
      Do i=1,NumBP
        n=kxb(i)
        k=codeT(n)
        Varbt(i,4)=0.0
        VarBT(i,3)=0.0
        VarBT(i,2)=0.0
        If (k.eq.-4) then
          If(abs(x(n)-xBSTEM).gt.SHADE.or.NShoot.eq.0) then
* If evaporation
* calc heat flux due to latent heat of evaporation
* given evaporation at it's potential rate
* units are J day-1
* RNS*3600*24/10000 is total energy going to soil in J day-1 cm-2
*  40.6 for mississippi
c           FELWR=max(0.0,(tmpr(i)-TairN)*FELWR)
c           If (VarBW(i,2).ge.0.0.and.RNS.gt.0.0) then
                VarBT(i,4)=RNS*3600.0*24.0/10000.0

*  calculate convected heat  if soil is warmer than air

              VarBT(i,3)=(FCSH+FELWR)*TAirN
              VarBT(i,2)=FCSH+FELWR
         else

*  the amount of heat added to the soil
*     is determined by air temperature.
*  units are millcal cm-2 min-1 C-1
               VarBT(i,2)=((0.058+1.7e-4*TAirN)+
     !                  (0.052*EXP(.058*TAirN)))
* the first term is thermal conductivity of the air.
* the second term is the conductivity of water vapor
* Change units to Joules cm-2 d-1
cdt change this for a test
               VarBT(i,2)=VarBT(i,2)*4.1856*60.0*24.0/1000
cdt /1000.0
               VarBT(i,3)=VarBT(i,2)*TairN
*  units are now J cm-2 d-1
*  heat into soil with rain
* end if for crop canopy
        endif

* add heat input with rain
c       If(VarBW(i,3).lt.0.0) then
c         VarBT(i,2)=VarBT(i,2)+varBW(i,3)*cWat
c         VarBT(i,3)=VarBT(i,3)+VarBW(i,3)*cWat*TAirN

* should this be + varbw*width?
c               Endif

       VarBT(i,1)=TAirN
        Endif
       Enddo
C ...............End of heat balance
c................... Gas movement
      Do i=1,NumBP
      n=KXB(i)
      k=CodeW(n)
      If(K.eq.4.or.K.eq.-4) then
        do jjj=1,NumG
          VarBG(i,jjj,2)=PG*Width(i)
          VarBG(i,jjj,3)=PG*GAIR(jjj)*Width(i)
        Enddo
      Endif
      Enddo
c................... End gas movement
c................... Furrow irrigation
c
      If(IfUR.eq.1) then
        Do i=1,NumFP
        k=NumF(i)
        QF=QF-Q(k)*Step
        Enddo
        If(QF.lt.0.0) then
          QF=0.
          IfUR=0
          Do i=1,NumFP
           k=NumF(i)
           CodeW(k)=-4
          Enddo
        Endif
      Endif
c................... End of the furrow irrigation

c................... This is the end of hourly calculations

      tNext(ModNum)=Time+period
      Endif      ! hourly loop
c

      RETURN
10    call errmes(im,il)
      END
C

ext: for

max_stars_repo_path: Soil Source/hourwea.for

max_stars_repo_name: cropbox-sk/MAIZSIM

max_issues_repo_path: Soil Source/hourwea.for

max_issues_repo_name: cropbox-sk/MAIZSIM

max_forks_repo_path: Soil Source/hourwea.for

max_forks_repo_name: cropbox-sk/MAIZSIM


      subroutine path_cha_res_read

      use constituent_mass_module
      use input_file_module
      use maximum_data_module
      use channel_data_module
      use hydrograph_module
      use sd_channel_module
      use organic_mineral_mass_module

      character (len=80) :: titldum
      character (len=80) :: header
      integer :: ipathi
      integer :: eof, imax
      logical :: i_exist              !none       |check to determine if file exists

      eof = 0

      !read all export coefficient data
      inquire (file=in_init%path_water, exist=i_exist)
      if (i_exist .or. in_init%path_water /= "null") then
        do
          open (107,file=in_init%path_water)
          read (107,*,iostat=eof) titldum
          if (eof < 0) exit
          read (107,*,iostat=eof) header
            if (eof < 0) exit
          imax = 0
          do while (eof == 0)
            read (107,*,iostat=eof) titldum     !name
            if (eof < 0) exit
            do ipath = 1, cs_db%num_paths
              read (107,*,iostat=eof) titldum   !water
              if (eof < 0) exit
              read (107,*,iostat=eof) titldum   !benthic
              if (eof < 0) exit
            end do
            imax = imax + 1
          end do

          db_mx%pathw_ini = imax

          allocate (path_water_ini(imax))
          allocate (path_init_name(imax))

          do ipathi = 1, imax
            allocate (path_water_ini(ipathi)%water(cs_db%num_paths))
            allocate (path_water_ini(ipathi)%benthic(cs_db%num_paths))
          end do
          rewind (107)
          read (107,*,iostat=eof) titldum
          if (eof < 0) exit

          do ipathi = 1, imax
            read (107,*,iostat=eof) header
            if (eof < 0) exit
            read (107,*,iostat=eof) path_init_name(ipathi)
            if (eof < 0) exit
            read (107,*,iostat=eof) titldum, path_water_ini(ipathi)%water
            if (eof < 0) exit
            read (107,*,iostat=eof) titldum, path_water_ini(ipathi)%benthic
            if (eof < 0) exit
          end do
          close (107)
          exit
        end do
      end if

      return
      end subroutine path_cha_res_read


text:
!*robodoc*u* tests/testACE
!  NAME
!    testACE
!  SYNOPSIS

!$Id: testACE.f90 526 2018-03-25 23:44:51Z mexas $

program testACE

!  PURPOSE
!    Test cgca_fwci/cgca_m2out - a debugging output routine.
!  DESCRIPTION
!    cgca_fwci/cgca_m2out collects coarray data from all images and
!    writes it out from image 1 in a formatted text file. A single
!    line contains the state of a single cells with its full local
!    global coordinates, i.e. its position in the coarray on this
!    image, and the image position in the coarray grid.
!    The data is dumped out as soon as the solidificaiton is finished.
!
!    Both grain and fracture layers are used and dumped.
!  NOTE
!    The program uses cgca_ins RND seed, to obtain a reproducible
!    model data. If this test is re-run on the same platform with
!    the same number of images, the results must be the same.
!  AUTHOR
!    Anton Shterenlikht
!  COPYRIGHT
!    See LICENSE
!  USES
!    cgca testaux
!  USED BY
!    Part of CGPACK test suite
!  SOURCE

use testaux

implicit none

logical( kind=ldef ), parameter :: yesdebug=.true., nodebug=.false.,   &
 periodicbc=.true.,  noperiodicbc=.false.

real, parameter :: gigabyte=real(2**30)

integer( kind=idef ) :: ir(3), img, nimgs,                             &
  ng,  & ! number of grains in the whole model
  c(3)   ! coarray dimensions

integer( kind=iarr ), allocatable :: space(:,:,:,:) [:,:,:]

integer( kind=ilrg ) :: icells, mcells

real( kind=rdef ) ::    &
 qual,                  & ! quality
 bsz0(3),               & ! the given "box" size
 bsz(3),                & ! updated "box" size
 dm,                    & ! mean grain size, linear dim, phys units
 lres,                  & ! linear resolution, cells per unit of length
 res                      ! resolutions, cells per grain
real( kind=rdef ), allocatable :: grt(:,:,:)[:,:,:]

logical( kind=ldef ) :: solid

!*********************************************************************72
! first executable statement

! physical dimensions of the box, assume mm
bsz0 = (/ 1.0, 1.0, 1.0 /)

! mean grain size, linear dimension, e.g. mean grain diameter, also mm
dm = 5.0e-1

! resolution
res = 1.0e5

! In this test set the number of images via the env var
! the code must be able to cope with any value >= 1.
   img = this_image()
 nimgs = num_images()

! do a check on image 1
if ( img .eq. 1 ) then

 ! print a banner
 call banner("ACE")

 ! print the parameter values
 call cgca_pdmp
 write (*,'(a,i0,a)') "running on ", nimgs, " images in a 3D grid"

end if

! want to sync here to make sure the banner is
! printed before the rest.
sync all

! each image calculates the coarray grid dimensions
call cgca_gdim( nimgs, ir, qual )

! calculate the resolution and the actual phys dimensions
! of the box
! subroutine cgca_cadim( bsz, res, dm, ir, c, lres, ng )
! c - coarray sizes
! ir - coarray grid sizes
bsz = bsz0
call cgca_cadim( bsz, res, dm, ir, c, lres, ng )

! total number of cells in a coarray
icells = int( c(1), kind=ilrg ) *    &
         int( c(2), kind=ilrg ) *    &
         int( c(3), kind=ilrg )

! total number of cells in the model
mcells = icells * int( nimgs, kind=ilrg )

if ( img .eq. 1 ) then
  write ( *, "(9(a,i0), 2(tr1, es9.2), 3(a, es9.2), a)" )              &
    "img: ", img  , " nimgs: ", nimgs, " (", c(1) ,                    &
    ","    , c(2) , ","       , c(3) , ")[", ir(1),                    &
    ","    , ir(2), ","       , ir(3), "] ", ng   ,                    &
    qual, lres,                                                        &
    " (", bsz(1), ",", bsz(2), ",", bsz(3), ")"
  write (*,*) "dataset sizes for ParaView", c*ir
  write (*,'(a,i0)') "Cells on each image: ", icells
  write (*,"(a, es10.2, a, i0)") "Total cells in the model (real): ",  &
    product( real(c) * real(ir) ), " (int): ", mcells
end if

! allocate space coarray with two layers, implicit SYNC ALL inside
call cgca_as(1, c(1), 1, c(2), 1, c(3), 1, ir(1), 1, ir(2), 1, 2, space)

! initialise the reproducible random number seed
call cgca_ins( yesdebug )

! Set initial values to all layers of space
space( :, :, :, cgca_state_type_grain ) = cgca_liquid_state
space( :, :, :, cgca_state_type_frac  ) = cgca_intact_state

! Allocate rotation tensors, implicit SYNC ALL inside
call cgca_art( 1, ng, 1, ir(1), 1, ir(2), 1, grt )

! Image 1 sets crack nuclei
if ( img .eq. 1 ) then

  ! set a single crack nucleus in the centre of the x3=max(x3) face
  space( c(1)/2, c(2)/2, c(3), cgca_state_type_frac )                  &
        [ ir(1)/2, ir(2)/2, ir(3) ] = cgca_clvg_state_100_edge
end if

! Set grain nuclei, SYNC ALL inside
call cgca_nr( space, ng, yesdebug )

! assign rotation tensors, SYNC ALL inside
call cgca_rt( grt )

! solidify, implicit SYNC ALL inside
!subroutine cgca_sld( coarray, periodicbc, iter, heartbeat, solid )
call cgca_sld( space, noperiodicbc, 0, 10, solid )

! Initiate grain boundaries. cgca_igb has no remote comms. Halo
! exchange is needed to update other images.
call cgca_igb( space )
sync all
call cgca_hxi( space )
sync all

! Smoothen the GB, several iterations.
! cgca_gbs has no remote comms.
call cgca_gbs( space )
sync all
call cgca_hxi( space )
sync all
call cgca_gbs( space )
sync all
call cgca_hxi( space )
sync all

if ( img .eq. 1 ) write (*,*) "dumping model to files"

  ! dump space arrays to files, only image 1 does it, all others
  ! wait at the barrier, hence sync needed
  call cgca_fwci( space, cgca_state_type_grain, "zg0.raw" )
  call cgca_fwci( space, cgca_state_type_grain, "zf0.raw" )

if ( img .eq. 1) write (*,*) "finished dumping model to files"

sync all

! deallocate all coaarrays, implicit sync all.
call cgca_ds( space )
call cgca_drt( grt )

if ( img .eq. 1 ) write (*,*) "Test ACE completed sucessfully"

end program testACE

!*roboend*

max_stars_repo_name: lcebaman/casup

max_issues_repo_name: lcebaman/casup

max_forks_repo_name: lcebaman/casup


text:
C Copyright(C) 2009-2017 National Technology & Engineering Solutions of
C Sandia, LLC (NTESS).  Under the terms of Contract DE-NA0003525 with
C NTESS, the U.S. Government retains 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     * Neither the name of NTESS 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============================================================================
      SUBROUTINE STCLST(MODULE)
C============================================================================

C   --*** STCLST ***  (BLOT) set color list mapping
C
C --Parameters:
C
C --  MODULE - character name of the subroutine module setting the colors
C
C  PURPOSE:  SETS THE DESIRED MAPPING OF THE DEFAULT COLORS.  THE ROUTINE
C            WAS WRITTEN TO ALLOW THE USER TO CHANGE THE MAPPING FOR
C            DIFFERENT PROGRAMS, SPECIFICALLY MAPPING ONE WAY FOR CONTOURS
C            AND PAINTING AND ANOTHER FOR PLOTTING.
C
      CHARACTER*(*) MODULE
      include 'params.blk'
      include 'cmap-lst.blk'
      CHARACTER*1 LSTMOD
      SAVE LSTMOD
      LOGICAL PLTICL

      CHARACTER*(MXSTLN) LIST1(NCOLOR), LIST2(NCOLOR)

      DATA LSTMOD / ' '/
      DATA LIST1 / 'BLACK   ', 'WHITE   ', 'RED     ',
     &   'GREEN   ', 'YELLOW  ', 'BLUE    ',
     &   'CYAN    ', 'MAGENTA ', '        ' /
      DATA LIST2 / 'BLACK   ', 'WHITE   ', 'BLUE    ',
     &   'MAGENTA ', 'CYAN    ', 'GREEN   ',
     &   'YELLOW  ', 'RED     ', '        ' /

C
C SET COLLST TO EITHER LIST DEPENDING ON THE VALUE OF MODULE
C
      IF(LSTMOD .NE. MODULE(1:1)) THEN
         LSTMOD = MODULE(1:1)

         IF(MODULE(1:1) .EQ. 'D' .OR. MODULE(1:1) .EQ. 'd') THEN
            DO 10 I = 1, NCOLOR
               COLLST(I) = LIST2(I)
10          CONTINUE
         ELSE
            DO 20 I = 1, NCOLOR
               COLLST(I) = LIST1(I)
20          CONTINUE
         END IF

         DO 30 I = 1, NCOLOR-2
            IF (PLTICL (COLLST(I+2), RCOLOR)) THEN
               COLMAP(I) = RCOLOR
            ELSE
               COLMAP(I) = -1
            ENDIF
30       CONTINUE
      END IF

      RETURN
      END

max_stars_repo_name: brian-kelley/seacas

max_issues_repo_name: brian-kelley/seacas

max_forks_repo_name: brian-kelley/seacas


text:
Name: Ellas Bates

Phone Number: (530) 8675309

Office: Davis, CA and San Jose, Alajuela, CR

Personality: Computer and political geek by trade. Well Wiki World Travelers traveled and read. Aviation buff. Loves motorcycles. Programs in a dead language.

I play RPG World Online when I get a chance...www.RPGWO.com... 2D graphics but great gameplay.

20070331 01:20:00 nbsp You arent by any chance a wiki:WikiPedia:Bo_Diddley Bo Diddley fan? Users/JasonAller

20070404 18:09:40 nbsp He is Otha Ellas Bates...no relationship Users/EllasBates


text:
      SUBROUTINE ASCM04 (NAME,IPHASE,ISOL,NOGO)
C
C     REDUCE COMMAND DMAP DATA
C
      INTEGER        COMND(6,1),XTRA(4),SUBNAM(2),ISAVE(54),
     1               RDMAP(18,23),RDMAP1(18,9),RDMAP2(18,9),
     2               RDMAP3(18,5),OCT(3,16),OCT1(3,16),PTBS(7,67),
     3               PTBS1(7,18),PTBS2(7,18),PTBS3(7,18),PTBS4(7,13)
      COMMON /ASDBD/ IRDM,NRDM,IXTRA,NXTRA,IOCT,NOCT,IPTBS,NPTBS,
     1               IPH,NPH,IDAT(935)
      EQUIVALENCE    (RDMAP1(1,1),RDMAP(1,1)),(RDMAP2(1,1),RDMAP(1,10)),
     1               (RDMAP3(1,1),RDMAP(1,19)),(OCT1(1,1),OCT(1,1)),
     2               (PTBS 1(1,1),PTBS(1, 1)),(PTBS2(1,1),PTBS(1,19)),
     3               (PTBS 3(1,1),PTBS(1,37)),(PTBS4(1,1),PTBS(1,55))
      DATA COMND   /
     1               4HREDU    , 23    ,  4    , 16    , 67    ,  0  /
      DATA SLASH   / 1H/       /
      DATA ISAVE   /
     1    1,14,1,  3,12,1,  3,16,3,  4, 5,1,  8, 8,2,  8,11,1,  9, 7,1,
     2   10,10,3, 12, 7,1, 14, 7,1, 16, 7,2, 16,10,1, 19, 7,2, 19,10,1,
     3   20, 7,2, 20,10,1, 21, 7,1, 23, 6,2    /
      DATA RDMAP 1 /
     1 4HREDU,4HCE  ,4H  CA,4HSECC,4H,GEO,4HM4/P,4HVNOA,4H,USS,4HTP,I,
     * 4HNSTP,4H/STP,4H/S,N,4H,DRY,4H!*PV,4HEC* ,4H$   ,4H    ,4H    ,
     2 4HCOND,4H    ,4H  LB,4HRSTP,4H,DRY,4H $  ,12*4H    ,
     3 4HSOFI,4H    ,4H  /K,4HNOA,,4HMNOA,4H,PNO,4HA,BN,4HOA,K,4H4NOA,
     * 4H/S,N,4H,DRY,4H!*NA,4HME00,4H0A*/,4H*KMT,4HX*!*,4HMMTX,4H*/  ,
     4 4H    ,4H    ,4H  *P,4HVEC*,4H!*BM,4HTX*/,4H*K4M,4HX* $,4H    ,
     * 4H    , 8*4H    ,
     5 4HCOND,4H    ,4H  LB,4HRSTP,4H,DRY,4H $  ,12*4H    ,
     6 4HSMP1,4H    ,4H  US,4HSTP,,4HKNOA,4H,,,/,4HGONO,4HA,KN,4HOB,K,
     * 4HONOA,4H,LON,4HOA,,,4H,,, ,4H$   , 4*4H    ,
     7 4HMERG,4HE   ,4H  GO,4HNOA,,4HINST,4HP,,,,4H,PVN,4HOA/G,4HNOA/,
     * 4H1/TY,4HP/2 ,4H$   , 6*4H    ,
     8 4HSOFO,4H    ,4H  ,G,4HNOA,,4HLONO,4HA,,,,4H//DR,4HY!*N,4HAME0,
     * 4H00A*,4H!*HO,4HRG*/,4H*LMT,4HX* $, 4*4H    ,
     9 4HSOFO,4H    ,4H  ,K,4HNOB,,4H,,,/,4H/DRY,4H!*NA,4HME00,4H0B*/,
     * 4H*KMT,4HX* $, 7*4H       /
      DATA RDMAP 2 /
     O 4HSOFI,4H    ,4H  /G,4HNOA,,4H,,,/,4HS,N,,4HDRY/,4H*NAM,4HE000,
     * 4HA*!*,4HHORG,4H* $ , 6*4H    ,
     1 4HMPY3,4H    ,4H  GN,4HOA,M,4HNOA,,4H/MNO,4HB/0/,4H0 $ ,4H    ,
     * 4H    , 8*4H    ,
     2 4HSOFO,4H    ,4H  ,M,4HNOB,,4H,,,/,4H/DRY,4H!*NA,4HME00,4H0B*/,
     * 4H*MMT,4HX* $, 7*4H    ,
     3 4HMPY3,4H    ,4H  GN,4HOA,B,4HNOA,,4H/BNO,4HB/0/,4H0 $ ,4H    ,
     * 4H    , 8*4H    ,
     4 4HSOFO,4H    ,4H  ,B,4HNOB,,4H,,,/,4H/DRY,4H!*NA,4HME00,4H0B*/,
     * 4H*BMT,4HX* $, 7*4H    ,
     5 4HMPY3,4H    ,4H  GN,4HOA,K,4H4NOA,4H,/K4,4HNOB/,4H0/0 ,4H$   ,
     * 4H    , 8*4H    ,
     6 4HSOFO,4H    ,4H  ,K,4H4NOB,4H,,,,,4H//DR,4HY!*N,4HAME0,4H00B*,
     * 4H!*K4,4HMX* ,4H$   , 6*4H    ,
     7 4HPART,4HN   ,4H  PN,4HOA,,,4HPVNO,4HA/PO,4HNOA,,4H,,/1,4H/1/2,
     * 4H $  , 8*4H    ,
     8 4HMPYA,4HD   ,4H  GN,4HOA,P,4HNOA,,4H/PNO,4HB/1/,4H1/0/,4H1 $ ,
     * 4H    , 8*4H        /
      DATA RDMAP 3 /
     9 4HSOFO,4H    ,4H  ,P,4HONOA,4H,,,,,4H//DR,4HY!*N,4HAME0,4H00A*,
     * 4H!*PO,4HVE* ,4H$   , 6*4H    ,
     O 4HSOFO,4H    ,4H  ,P,4HVNOA,4H,,,,,4H//DR,4HY!*N,4HAME0,4H00A*,
     * 4H!*UP,4HRT* ,4H$   , 6*4H    ,
     1 4HSOFO,4H    ,4H  ,P,4HNOB,,4H,,,/,4H/DRY,4H!*NA,4HME00,4H0B*/,
     * 4H*PVE,4HC* $, 7*4H    ,
     2 4HLABE,4HL   ,4H  LB,4HRSTP,4H $  ,13*4H    ,
     3 4HLODA,4HPP  ,4H  PN,4HOB,P,4HONOA,4H/!*N,4HAME0,4H00B*,4H/S,N,
     * 4H,DRY,4H $  , 7*4H      /
      DATA XTRA    /
     1               4HOUTP,4HNAME,4HBOUN,4HRSAV        /
      DATA OCT 1   /
     1                    6    ,         0    ,         1  ,
     2                    7    ,         0    ,         1  ,
     3                    8    ,         0    ,         1  ,
     4                    9    ,         0    ,         1  ,
     5                   10    ,         1    ,        62  ,
     6                   11    ,         0    ,         2  ,
     7                   12    ,         0    ,         2  ,
     8                   13    ,         0    ,        16  ,
     9                   14    ,         0    ,        16  ,
     O                   15    ,         0    ,        32  ,
     1                   16    ,         0    ,        32  ,
     2                   17    ,         0    ,        12  ,
     3                   18    ,         0    ,        12  ,
     4                   19    ,         0    ,        12  ,
     5                   21    ,         0    ,        12  ,
     6                   23    ,         0    ,         8  /
      DATA PTBS 1  /
     1          1  , 24  , 26  ,  3  ,4HNONA  ,         0  ,  0  ,
     2          1  , 32  , 32  ,  3  ,4HSTEP  ,         0  ,  0  ,
     3          1  , 38  , 38  ,  3  ,4HSTEP  ,         0  ,  0  ,
     4          1  , 41  , 42  ,  3  ,4HSTEP  ,         0  ,  0  ,
     5          1  , 53  , 55  ,  4  ,4HPITM  ,         0  ,  0  ,
     6          2  , 14  , 14  ,  3  ,4HSTEP  ,         0  ,  0  ,
     7          3  , 12  , 13  ,  3  ,4HNONA  ,         1  , -1  ,
     8          3  , 17  , 18  ,  3  ,4HNONA  ,         2  , -1  ,
     9          3  , 22  , 23  ,  3  ,4HNONA  ,        12  , -1  ,
     O          3  , 27  , 28  ,  3  ,4HNONA  ,        16  , -1  ,
     1          3  , 32  , 34  ,  3  ,4HNONA  ,        32  , -1  ,
     2          3  , 45  , 47  ,  8  ,4HNAMA  ,         0  ,  0  ,
     3          4  , 11  , 12  ,  4  ,4HPITM  ,         0  ,  0  ,
     4          5  , 14  , 14  ,  3  ,4HSTEP  ,         0  ,  0  ,
     5          6  , 11  , 13  ,  3  ,4HSTEP  ,         0  ,  0  ,
     6          6  , 17  , 18  ,  3  ,4HNONA  ,         0  ,  0  ,
     7          6  , 25  , 27  ,  3  ,4HNONA  ,         0  , -1  ,
     8          6  , 31  , 32  ,  3  ,4HNONB  ,         0  , -1  /
      DATA PTBS 2  /
     1          6  , 36  , 38  ,  3  ,4HNONA  ,         0  ,  0  ,
     2          6  , 42  , 44  ,  3  ,4HNONA  ,         0  ,  0  ,
     3          7  , 11  , 13  ,  3  ,4HNONA  ,         0  ,  0  ,
     4          7  , 17  , 19  ,  3  ,4HSTEP  ,         0  ,  0  ,
     5          7  , 26  , 28  ,  3  ,4HNONA  ,         0  ,  0  ,
     6          7  , 32  , 33  ,  3  ,4HNONA  ,         0  , -1  ,
     7          7  , 38  , 39  ,  3  ,4HPREC  ,         0  ,  0  ,
     8          8  , 12  , 13  ,  3  ,4HNONA  ,         0  ,  0  ,
     9          8  , 17  , 19  ,  3  ,4HNONA  ,         0  ,  0  ,
     O          8  , 17  , 21  ,  0  ,4HRSAV  ,         0  ,  0  ,
     1          8  , 30  , 32  ,  8  ,4HNAMA  ,         0  ,  0  ,
     2          8  , 48  , 54  ,  0  ,4HRSAV  ,         0  ,  0  ,
     3          9  , 12  , 13  ,  3  ,4HNONB  ,         0  ,  0  ,
     4          9  , 25  , 27  ,  8  ,4HNAMB  ,         0  ,  0  ,
     5         10  , 12  , 13  ,  3  ,4HNONA  ,         0  , -1  ,
     6         10  , 28  , 30  ,  8  ,4HNAMA  ,         0  ,  0  ,
     7         11  , 11  , 12  ,  3  ,4HNONA  ,         0  ,  0  ,
     8         11  , 16  , 17  ,  3  ,4HNONA  ,         0  ,  0  /
      DATA PTBS 3  /
     1         11  , 22  , 23  ,  3  ,4HNONB  ,         0  , -1  ,
     2         12  , 12  , 13  ,  3  ,4HNONB  ,         0  ,  0  ,
     3         12  , 25  , 27  ,  8  ,4HNAMB  ,         0  ,  0  ,
     4         13  , 11  , 12  ,  3  ,4HNONA  ,         0  ,  0  ,
     5         13  , 16  , 17  ,  3  ,4HNONA  ,         0  ,  0  ,
     6         13  , 22  , 23  ,  3  ,4HNONB  ,         0  , -1  ,
     7         14  , 12  , 13  ,  3  ,4HNONB  ,         0  ,  0  ,
     8         14  , 25  , 27  ,  8  ,4HNAMB  ,         0  ,  0  ,
     9         15  , 11  , 12  ,  3  ,4HNONA  ,         0  ,  0  ,
     O         15  , 16  , 18  ,  3  ,4HNONA  ,         0  ,  0  ,
     1         15  , 23  , 25  ,  3  ,4HNONB  ,         0  , -1  ,
     2         16  , 12  , 14  ,  3  ,4HNONB  ,         0  ,  0  ,
     3         16  , 26  , 28  ,  8  ,4HNAMB  ,         0  ,  0  ,
     4         17  , 11  , 12  ,  3  ,4HNONA  ,         0  ,  0  ,
     5         17  , 17  , 19  ,  3  ,4HNONA  ,         0  ,  0  ,
     6         17  , 23  , 25  ,  3  ,4HNONA  ,         0  ,  0  ,
     7         18  , 11  , 12  ,  3  ,4HNONA  ,         0  ,  0  ,
     8         18  , 16  , 17  ,  3  ,4HNONA  ,         0  ,  0  /
      DATA PTBS 4  /
     1         18  , 22  , 23  ,  3  ,4HNONB  ,         0  , -1  ,
     2         19  , 12  , 14  ,  3  ,4HNONA  ,         0  ,  0  ,
     3         19  , 26  , 28  ,  8  ,4HNAMA  ,         0  ,  0  ,
     4         19  , 37  , 39  ,  4  ,4HPOIT  ,         0  ,  0  ,
     5         20  , 12  , 14  ,  3  ,4HNONA  ,         0  ,  0  ,
     6         20  , 26  , 28  ,  8  ,4HNAMA  ,         0  ,  0  ,
     7         21  , 12  , 13  ,  3  ,4HNONB  ,         0  ,  1  ,
     8         21  , 25  , 27  ,  8  ,4HNAMB  ,         0  ,  0  ,
     9         21  , 36  , 38  ,  4  ,4HPITM  ,         0  ,  0  ,
     O         22  , 14  , 14  ,  3  ,4HSTEP  ,         0  ,  0  ,
     1         23  , 11  , 12  ,  3  ,4HNONB  ,         0  ,  1  ,
     2         23  , 16  , 18  ,  3  ,4HNONA  ,         0  ,  1  ,
     3         23  , 22  , 24  ,  8  ,4HNAMB  ,         0  ,  0  /
C
      DATA SUBNAM  / 4HASCM,2H04  /
C
C     RESTORE TO ORIGINAL DATA BY REPLACEING !* BY /* IN RDMAP ARRAY
C     (SEE ASCM01 FOR EXPLANATION)
C
      DO 20 L = 1,51,3
      I = ISAVE(L+1)
      J = ISAVE(L  )
      K = ISAVE(L+2)
      RDMAP(I,J) = KHRFN1(RDMAP(I,J),K,SLASH,1)
   20 CONTINUE
C
C     VALIDATE COMMAND AND SET POINTERS
C
      IF (NAME .NE. COMND(1,1)) GO TO 70
      ICOMND = 1
      IRDM   = 1
      NRDM   = COMND(2,ICOMND)
      IXTRA  = IRDM  + 18*NRDM
      NXTRA  = COMND(3,ICOMND)
      IOCT   = IXTRA + NXTRA
      NOCT   = COMND(4,ICOMND)
      IPTBS  = IOCT  + 3*NOCT
      NPTBS  = COMND(5,ICOMND)
      IPH    = IPTBS + 7*NPTBS
      NPH    = COMND(6,ICOMND)
C
C     MOVE RDMAP DATA
C
      K = 0
      IF (NRDM .EQ. 0) GO TO 35
      DO 30 J = 1,NRDM
      DO 30 I = 1,18
      K = K + 1
   30 IDAT(K) = RDMAP(I,J)
   35 CONTINUE
C
C     MOVE XTRA DATA
C
      IF (NXTRA .EQ. 0) GO TO 45
      DO 40 I = 1,NXTRA
      K = K + 1
   40 IDAT(K) = XTRA(I)
   45 CONTINUE
C
C     MOVE OCT DATA
C
      IF (NOCT .EQ. 0) GO TO 55
      DO 50 J = 1,NOCT
      DO 50 I = 1,3
      K = K + 1
   50 IDAT(K) = OCT(I,J)
   55 CONTINUE
C
C     MOVE PTBS DATA
C
      IF (NPTBS .EQ. 0) GO TO 65
      DO 60 J = 1,NPTBS
      DO 60 I = 1,7
      K = K + 1
   60 IDAT(K) = PTBS(I,J)
   65 CONTINUE
C
      RETURN
C
C     INPUT ERROR
C
   70 CALL MESAGE (7,0,SUBNAM)
      NOGO = 1
      RETURN
C
      END


text:
PROGRAM F033

    ! Copyright 2021 Melwyn Francis Carlo

    IMPLICIT NONE

    INTEGER, DIMENSION(4, 2) :: PRODUCTS_LIST

    CHARACTER (LEN=10) :: TEMP_STR

    REAL    :: VALUE_1, VALUE_2

    INTEGER :: COUNT_VAL = 1

    INTEGER :: I, J, TEMP_NUM_11, TEMP_NUM_12, TEMP_NUM_21,                    &
                     TEMP_NUM_22, TEMP_NUM_3,  TEMP_NUM_4

    LOGICAL :: REDUCED

    DO I = 10, 50

        DO J = 10, 99

            IF ((I >= J) .OR. (MOD(I, 10) == 0) .OR. (MOD(J, 10) == 0)) CYCLE

            WRITE(TEMP_STR, '(I10)') I
            TEMP_STR = ADJUSTL(TEMP_STR)

            TEMP_NUM_11 = IACHAR(TEMP_STR(1:1)) - IACHAR('0')
            TEMP_NUM_12 = IACHAR(TEMP_STR(2:2)) - IACHAR('0')

            WRITE(TEMP_STR, '(I10)') J
            TEMP_STR = ADJUSTL(TEMP_STR)

            TEMP_NUM_21 = IACHAR(TEMP_STR(1:1)) - IACHAR('0')
            TEMP_NUM_22 = IACHAR(TEMP_STR(2:2)) - IACHAR('0')

            TEMP_NUM_3 = 99999
            TEMP_NUM_4 = 99999

             IF  ((TEMP_NUM_11 == TEMP_NUM_21)                                 &
            .AND. (TEMP_NUM_12 /= TEMP_NUM_22)) THEN

                TEMP_NUM_3 = TEMP_NUM_12
                TEMP_NUM_4 = TEMP_NUM_22

            ELSE IF ((TEMP_NUM_11 == TEMP_NUM_22)                              &
            .AND.    (TEMP_NUM_12 /= TEMP_NUM_21)) THEN

                TEMP_NUM_3 = TEMP_NUM_12
                TEMP_NUM_4 = TEMP_NUM_21

            ELSE IF ((TEMP_NUM_12 == TEMP_NUM_21)                              &
            .AND.    (TEMP_NUM_11 /= TEMP_NUM_22)) THEN

                TEMP_NUM_3 = TEMP_NUM_11
                TEMP_NUM_4 = TEMP_NUM_22

            ELSE IF ((TEMP_NUM_12 == TEMP_NUM_22)                              &
            .AND.    (TEMP_NUM_11 /= TEMP_NUM_21)) THEN

                TEMP_NUM_3 = TEMP_NUM_11
                TEMP_NUM_4 = TEMP_NUM_21

            END IF

            IF (TEMP_NUM_3 /= 99999) THEN

                VALUE_1 = (FLOOR(REAL(I) * 1.0E5 / REAL(J))) / 1.0E5;
                VALUE_2 = (FLOOR(REAL(TEMP_NUM_3) * 1.0E5 /                    &
                           REAL(TEMP_NUM_4))) / 1.0E5;

                IF  ((VALUE_1 - VALUE_2) <  1.0E-5                             &
               .AND. (VALUE_1 - VALUE_2) > -1.0E-5) THEN

                    PRODUCTS_LIST(COUNT_VAL, 1) = I
                    PRODUCTS_LIST(COUNT_VAL, 2) = J
                    COUNT_VAL = COUNT_VAL + 1
                END IF

            END IF

        END DO

    END DO

    PRODUCTS_LIST(1, 1) = PRODUCTS_LIST(1, 1) * PRODUCTS_LIST(2, 1)            &
                        * PRODUCTS_LIST(3, 1) * PRODUCTS_LIST(4, 1)

    PRODUCTS_LIST(1, 2) = PRODUCTS_LIST(1, 2) * PRODUCTS_LIST(2, 2)            &
                        * PRODUCTS_LIST(3, 2) * PRODUCTS_LIST(4, 2)

    DO WHILE ((PRODUCTS_LIST(1, 1) /= 1) .AND. (PRODUCTS_LIST(1, 2) /= 1))

        I = 2
        REDUCED = .FALSE.

        DO WHILE (.NOT. REDUCED)

            IF ((MOD(PRODUCTS_LIST(1, 1), I) == 0)                             &
          .AND. (MOD(PRODUCTS_LIST(1, 2), I) == 0)) THEN

                PRODUCTS_LIST(1, 1) = PRODUCTS_LIST(1, 1) / I
                PRODUCTS_LIST(1, 2) = PRODUCTS_LIST(1, 2) / I
                REDUCED = .TRUE.

            ELSE
                I = I + 1
            END IF

        END DO

    END DO

    PRINT ('(I0)'), PRODUCTS_LIST(1, 2)

END PROGRAM F033

max_stars_repo_name: melwyncarlo/ProjectEuler

max_issues_repo_name: melwyncarlo/ProjectEuler

max_forks_repo_name: melwyncarlo/ProjectEuler


text:
   module netcdf_utils

   use anthro_types, only : glb_att

   implicit none

   integer :: ngatts
   type(glb_att), allocatable :: attrs(:)

   private
   public :: get_glb_atts, set_glb_atts, handle_ncerr
   public :: ngatts

!---------------------------------------------------------------------
!  include files
!---------------------------------------------------------------------
   include 'netcdf.inc'

   CONTAINS

   subroutine get_glb_atts( ncid )
!---------------------------------------------------------------------
!   read the global attributes
!---------------------------------------------------------------------

   implicit none

   integer, intent(in) :: ncid

!---------------------------------------------------------------------
!   local variables
!---------------------------------------------------------------------
   integer              :: m
   integer              :: astat
   integer              :: attr_len
   character(len=132)   :: message
   character(len=132)   :: attr_name

!---------------------------------------------------------------------
!   get global attr count
!---------------------------------------------------------------------
   message = 'glb_attr: Failed to get glb attr count'
   call handle_ncerr( nf_inq_natts( ncid, ngatts ), message )
!---------------------------------------------------------------------
!   allocate variables
!---------------------------------------------------------------------
   if( ngatts > 0 ) then
     if( allocated( attrs ) ) then
       call dealloc_glb_atts
     endif
     allocate( attrs(ngatts),stat=astat )
     if( astat /= 0 ) then
       write(*,*) 'glb_attr: failed to allocate type glb_att'
       stop 'Alloc err'
     endif
     attrs(:)%name = ' '
   endif
!---------------------------------------------------------------------
!   loop over glb attributes
!---------------------------------------------------------------------
glb_attr_loop : &
   do m = 1,ngatts
     write(message,*) 'glb_attr: Failed to get glb attr # ',m,' name'
     call handle_ncerr( nf_inq_attname( ncid, nf_global, m, attr_name ), message )
     attrs(m)%name = attr_name
     write(message,*) 'glb_attr: Failed to get glb attr # ',m,' type,len'
     call handle_ncerr( nf_inq_att( ncid, nf_global, trim(attr_name), attrs(m)%type, attr_len ), message )
     attrs(m)%len = attr_len
     message = 'glb_attr: Failed to get ' // trim(attr_name)
     select case( attrs(m)%type )
       case( nf_byte )
         allocate( attrs(m)%attr_byte(attr_len),stat=astat )
         if( astat /= 0 ) then
           write(*,*) 'glb_attr: failed to allocate attr_byte'
           stop 'Alloc err'
         endif
         call handle_ncerr( nf_get_att_int1( ncid, nf_global, trim(attr_name), attrs(m)%attr_byte ), message )
       case( nf_char )
         attrs(m)%attr_char = ' '
         call handle_ncerr( nf_get_att_text( ncid, nf_global, trim(attr_name), attrs(m)%attr_char ), message )
       case( nf_short )
         allocate( attrs(m)%attr_short(attr_len),stat=astat )
         if( astat /= 0 ) then
           write(*,*) 'glb_attr: failed to allocate attr_short'
           stop 'Alloc err'
         endif
         call handle_ncerr( nf_get_att_int2( ncid, nf_global, trim(attr_name), attrs(m)%attr_short ), message )
       case( nf_int )
         allocate( attrs(m)%attr_int(attr_len),stat=astat )
         if( astat /= 0 ) then
           write(*,*) 'glb_attr: failed to allocate attr_int'
           stop 'Alloc err'
         endif
         call handle_ncerr( nf_get_att_int( ncid, nf_global, trim(attr_name), attrs(m)%attr_int ), message )
       case( nf_float )
         allocate( attrs(m)%attr_real(attr_len),stat=astat )
         if( astat /= 0 ) then
           write(*,*) 'glb_attr: failed to allocate attr_real'
           stop 'Alloc err'
         endif
         call handle_ncerr( nf_get_att_real( ncid, nf_global, trim(attr_name), attrs(m)%attr_real ), message )
       case( nf_double )
         allocate( attrs(m)%attr_dbl(attr_len),stat=astat )
         if( astat /= 0 ) then
           write(*,*) 'glb_attr: failed to allocate attr_dbl'
           stop 'Alloc err'
         endif
         call handle_ncerr( nf_get_att_double( ncid, nf_global, trim(attr_name), attrs(m)%attr_dbl ), message )
     end select
   end do glb_attr_loop

   end subroutine get_glb_atts

   subroutine set_glb_atts( ncid )
!---------------------------------------------------------------------
!   set the global attributes
!---------------------------------------------------------------------

   implicit none

   integer, intent(in) :: ncid

!---------------------------------------------------------------------
!   local variables
!---------------------------------------------------------------------
   integer              :: m
   integer              :: ncd_err
   integer              :: attr_len
   integer              :: attr_xtype
   integer              :: slen
   character(len=132)   :: message
   character(len=132)   :: attr_name

!---------------------------------------------------------------------
!   loop over glb attributes
!---------------------------------------------------------------------
glb_attr_loop : &
   do m = 1,ngatts
     attr_name = trim(attrs(m)%name)
     if( trim(attr_name) == 'TITLE' .or. trim(attr_name) == 'START_DATE' .or. &
         trim(attr_name) == 'SIMULATION_START_DATE' ) then
       cycle
     endif
     slen      = len_trim(attr_name)
     write(message,*) 'set_glb_att: Failed to define glb att ',trim(attr_name)
     attr_len   = attrs(m)%len
     attr_xtype = attrs(m)%type
     select case( attrs(m)%type )
       case( nf_byte )
         ncd_err = nf_put_att_int1( ncid, nf_global, attr_name(:slen), attr_xtype, attr_len, attrs(m)%attr_byte )
       case( nf_char )
         ncd_err = nf_put_att_text( ncid, nf_global, attr_name(:slen), attr_len, attrs(m)%attr_char )
       case( nf_short )
         ncd_err = nf_put_att_int2( ncid, nf_global, attr_name(:slen), attr_xtype, attr_len, attrs(m)%attr_short )
       case( nf_int )
         ncd_err = nf_put_att_int( ncid, nf_global, attr_name(:slen), attr_xtype, attr_len, attrs(m)%attr_int )
       case( nf_float )
         ncd_err = nf_put_att_real( ncid, nf_global, attr_name(:slen), attr_xtype, attr_len, attrs(m)%attr_real )
       case( nf_double )
         ncd_err = nf_put_att_double( ncid, nf_global, attr_name(:slen), attr_xtype, attr_len, attrs(m)%attr_dbl )
     end select
     call handle_ncerr( ncd_err, message )
   end do glb_attr_loop

   end subroutine set_glb_atts

   subroutine dealloc_glb_atts
!---------------------------------------------------------------------
!   deallocate variables
!---------------------------------------------------------------------

   integer :: m

   if( allocated( attrs ) ) then
     do m = 1,ngatts
       select case( attrs(m)%type )
         case( nf_byte )
           deallocate( attrs(m)%attr_byte )
         case( nf_short )
           deallocate( attrs(m)%attr_short )
         case( nf_int )
           deallocate( attrs(m)%attr_int )
         case( nf_float )
           deallocate( attrs(m)%attr_real )
         case( nf_double )
           deallocate( attrs(m)%attr_dbl )
       end select
     end do
     deallocate( attrs )
   endif

   end subroutine dealloc_glb_atts

   subroutine handle_ncerr( ret, mes )
!---------------------------------------------------------------------
!	... netcdf error handling routine
!---------------------------------------------------------------------

!---------------------------------------------------------------------
!	... dummy arguments
!---------------------------------------------------------------------
   integer, intent(in) :: ret
   character(len=*), intent(in) :: mes

   if( ret /= nf_noerr ) then
      write(*,*) nf_strerror( ret )
      stop 'netcdf error'
   endif

   end subroutine handle_ncerr

   end module netcdf_utils

max_stars_repo_name: mattldawson/IPT

max_issues_repo_name: mattldawson/IPT

max_forks_repo_name: mattldawson/IPT


text:
!! Regularized Optimization for Hypers-spectral Analysis (ROHSA)
program ROHSA

  use mod_constants
  use mod_start
  use mod_inout
  use mod_rohsa

  implicit none

  logical :: noise           !! if false --> STD map computed by ROHSA with lstd and ustd (if true given by the user)
  logical :: regul           !! if true --> activate regulation
  logical :: descent         !! if true --> activate hierarchical descent to initiate the optimization
  logical :: save_grid       !! save grid of fitted parameters at each step of the multiresolution process
  logical :: lym             !! if true --> activate 2-Gaussian decomposition for Lyman alpha nebula emission
  logical :: init_spec       !! if true --> use params mean spectrum with input
  logical :: init_grid       !! if true --> use fileinit to give the initialization of the last grid

  integer :: n_gauss         !! number of gaussian to fit
  integer :: n_gauss_add     !! number of gaussian to add at each step
  integer :: m               !! number of corrections used in the limited memory matrix by LBFGS-B
  integer :: lstd            !! lower bound to compute the standard deviation map of the cube (if noise .eq. false)
  integer :: ustd            !! upper bound to compute the standrad deviation map of the cube (if noise .eq. false)
  integer :: iprint          !! print option
  integer :: iprint_init     !! print option init
  integer :: maxiter         !! max iteration for L-BFGS-B alogorithm
  integer :: maxiter_init    !! max iteration for L-BFGS-B alogorithm (init mean spectrum)

  real(xp) :: lambda_amp     !! lambda for amplitude parameter
  real(xp) :: lambda_mu      !! lamnda for mean position parameter
  real(xp) :: lambda_sig     !! lambda for dispersion parameter

  real(xp) :: lambda_var_amp !! lambda for variance amplitude parameter
  real(xp) :: lambda_var_mu  !! lambda for variance mean position parameter
  real(xp) :: lambda_var_sig !! lambda for variance dispersion parameter

  real(xp) :: lambda_lym_sig !! lambda for variance dispersion parameter

  real(xp) :: amp_fact_init  !! times max amplitude of additional Gaussian
  real(xp) :: sig_init       !! dispersion of additional Gaussian
  real(xp) :: lb_sig_init    !! lower bound sigma init
  real(xp) :: ub_sig_init    !! upper bound sigma init
  real(xp) :: lb_sig         !! lower bound sigma
  real(xp) :: ub_sig         !! upper bound sigma

  character(len=512) :: filename_parameters !! name of the parameters file (default parameters.txt)
  character(len=512) :: filename            !! name of the data file
  character(len=512) :: fileout             !! name of the output result
  character(len=512) :: timeout             !! name of the output result
  character(len=512) :: filename_noise      !! name of the file with STD map (if noise .eq. true)
  character(len=512) :: filename_init_spec
  character(len=512) :: fileinit           !! name of the file with init last grid
  character(len=8)   :: init_option !!Init ROHSA with the mean or the std spectrum

  real(xp) :: start, finish

  real(xp), dimension(:,:,:), allocatable :: data        !! initial fits data
  real(xp), dimension(:,:,:), allocatable :: data_init   !! initial fits data init full grid
  real(xp), dimension(:,:), allocatable   :: std_cube    !! standard deviation map fo the cube is given by the user
  real(xp), dimension(:), allocatable     :: params_init     !! standard deviation map fo the cube is given by the user

  call cpu_time(start)

  !Print header and get filename in argument
  call get_command_argument(1, filename_parameters)

  !Default user parameters
  n_gauss = 6

  lambda_amp = 1._xp
  lambda_mu = 1._xp
  lambda_sig = 1._xp

  lambda_var_amp = 0._xp
  lambda_var_mu = 0._xp
  lambda_var_sig = 1._xp

  lambda_lym_sig = 0._xp

  amp_fact_init = 2._xp/3._xp
  sig_init = 5._xp
  lb_sig_init = 0.001_xp
  ub_sig_init = 100._xp
  lb_sig = 0.001_xp
  ub_sig = 100._xp
  maxiter_init = 15000
  maxiter = 800
  m = 10
  noise = .false.
  regul = .true.
  descent = .false.
  lstd = 1; ustd = 20
  init_option = "mean"
  iprint = -1
  iprint_init = -1
  save_grid = .true.
  lym = .false.
  init_grid = .false.
  init_spec = .false.

  !Read parameters
  call read_parameters(filename_parameters, filename, fileout, timeout, filename_noise, filename_init_spec, &
       n_gauss, lambda_amp, lambda_mu, lambda_sig, lambda_var_amp, lambda_var_mu, lambda_var_sig, lambda_lym_sig, &
       amp_fact_init, sig_init, lb_sig_init, ub_sig_init, lb_sig, ub_sig, init_option, maxiter_init, maxiter, &
       m, noise, regul, descent, lstd, ustd, iprint, iprint_init, save_grid, lym, init_grid, fileinit, init_spec)

  !Call header
  call header()

  print*, "filename = '",trim(filename),"'"

  !Load data
  call read_cube(filename, data)

  if (init_spec .eqv. .true.) then
     call read_array(filename_init_spec, params_init)
  end if

  !Load grid init
  if (init_grid .eqv. .true.) then
     call read_output(fileinit, data_init, n_gauss, shape(data))
  end if

  if (noise .eqv. .true.) then
     if (filename_noise == " ") then
        print*, "--> noise = .true. (no input rms map)"
     end if
     call read_map(filename_noise, std_cube)
  end if

  !Call ROHSA subroutine
  call main_rohsa(data, std_cube, fileout, timeout, n_gauss, lambda_amp, lambda_mu, lambda_sig, &
       lambda_var_amp, lambda_var_mu, lambda_var_sig, lambda_lym_sig, amp_fact_init, sig_init, lb_sig_init, &
       ub_sig_init, lb_sig, ub_sig, maxiter_init, maxiter, m, noise, regul, descent, lstd, ustd, init_option, &
       iprint, iprint_init, save_grid, lym, init_grid, fileinit, data_init, params_init, init_spec)

  call ender()

  call cpu_time(finish)
  print '("Time = ",f6.3," seconds.")',finish-start

end program ROHSA

max_stars_repo_name: antoinemarchal/ROHSA

max_issues_repo_name: antoinemarchal/ROHSA

max_forks_repo_name: antoinemarchal/ROHSA


!***********************************************************************************************************************************
!**                                        S U B R O U T I N E   H E A T  E X C H A N G E                                         **
!***********************************************************************************************************************************

SUBROUTINE HEAT_EXCHANGE
  USE GLOBAL; USE GDAYC; USE SURFHE; USE TVDC; USE SHADEC; USE PREC
  IMPLICIT NONE

! Type declaration
  REAL     :: JDAY, LOCAL
  REAL(R8) :: TSUR
  REAL(R8) :: MPS_TO_MPH,W_M2_TO_BTU_FT2_DAY,FLUX_BR_TO_FLUX_SI,BOWEN_CONSTANT
  REAL(R8) :: DEG_F,DEG_C,BTU_FT2_DAY_TO_W_M2
  REAL(R8) :: STANDARD,HOUR,TAUD,SINAL,A0,TDEW_F,TAIR_F,SRO_BR,WIND_MPH,WIND2M,ACONV,BCONV
  REAL(R8) :: TSTAR,BETA,FW,RA,ETP,EA,ES,TAIRV,DTV,DTVL
  INTEGER :: IDAY,J
  !LOGICAL :: BIRDM

  REAL(R8)  :: r, solarzen, RToD, DToR, solaralt, Uo   ! SW 9/14/2018
  Real      :: directHz, diffuseHz, globalHz
  !REAL      :: albedo, dust, At, Ba, K1, AOD380nm,  AOD500nm
  REAL(R8)  :: B0, B1, BB, Beta, jme
  REAL(R8)  :: A, B, jd, t, LO, M, Mrad, sinM, sin2M, sin3M, c, TLO, v
  REAL(R8)  :: omega, lamda, second, e0, e, tananum, tanadenom, alpha, sint, theta, ep
  REAL(R8)  :: y, sin2LO, cos2LO, sin4LO, Etime, EQT
    real ::         day_decimal,Jdaym   ! ,HOUR
  REAL(R8)  :: Zone, daySavings, hh, mm, ss, timenow
  REAL(R8)  :: solarTimeFix, trueSolarTime, hourangle, harad, csz, zenith, azDenom, azRad
  REAL(R8)  :: azimuth, exoatmElevation, step1, step2, step3, refractionCorrection
  SAVE

! Data declaration

  DATA MPS_TO_MPH          /2.23714D0/, W_M2_TO_BTU_FT2_DAY /7.60796D0/, FLUX_BR_TO_FLUX_SI /0.23659D0/
  DATA BTU_FT2_DAY_TO_W_M2 /0.1314D0/
  DATA BOWEN_CONSTANT      /0.47D0/

  DATA PI/3.14159265358979323846/
  DATA SROCON/1367.00/      !W/m2

  !BIRDM= .FALSE.
  !INQUIRE(FILE='BirdMeeseSR.npt',EXIST=BIRDM)     ! SW 9/19/2018
  !IF(BIRDM)THEN
  !
  !    OPEN(9990,FILE='BirdMeeseSR.npt',STATUS='OLD')
  !    READ(9990,*)
  !    READ(9990,*)ALBEDO
  !    READ(9990,*)
  !    READ(9990,*)DUST
  !    READ(9990,*)
  !    READ(9990,*)ALBEDO
  !    READ(9990,*)
  !    READ(9990,*)AT
  !    READ(9990,*)
  !    READ(9990,*)BA
  !    READ(9990,*)
  !    READ(9990,*)K1
  !    READ(9990,*)
  !    READ(9990,*)AOD380nm
  !    READ(9990,*)
  !    READ(9990,*)AOD500nm
  !    CLOSE(9990)
  !
  !ENDIF

RETURN

!***********************************************************************************************************************************
!**                                            S H O R T  W A V E  R A D I A T I O N                                              **
!***********************************************************************************************************************************

ENTRY SHORT_WAVE_RADIATION (JDAY)
  IF(BIRDM)THEN    ! BIRD/MEESE SHORT WAVE SOLAR FORMULAE    SW 9/14/2018

    RToD = (180.0/PI)
    DToR = (PI/180.0)
    TZ(i)= TZ(i)*-1.0

! Calculates Solar position and radiation from Meeus code, currently set up for clear skies only
    CLOUD=0.0
!    year=yeard(i,k)
    year = 2004
    LEAP_YEAR             = MOD(year,4) == 0
    Jdaym=Jdayd(i,k)+TZ(i)/24.0
    call GREGORIAN_DATE(i,k)

    HOUR  =  (Jdaym-INT(Jdaym))*24
    day_decimal = day + HOUR/24.0

    if(month <= 2)then
     year = year - 1
     month = month + 12
    end if

    A = FLOOR(year/100.0)
    B = 2 - A + FLOOR(A/4)
    jd = FLOOR(365.25*(year+4716.0)) + FLOOR(30.6001*(month+1))+day_decimal + B - 1524.5         ! Eq 7.1, pg 61

! gp (Greg Pellitier) put the year and month back where they belong
    if(month == 13)then
     month = 1
     year = year + 1
    end if
    if(month == 14)then
     month = 2
     year = year + 1
    end if

! time in Julian centuries from epoch J2000.0
      t = (jd-2451545.0)/36525.0                                                               ! Eq 25.1, pg 163
! julian ephemeris millennium
      jme  = t/10.0

! calculate the Geometric Mean Longitude of the Sun, degrees
      LO = 280.46646+t*(36000.76983+0.0003032*t)                                               ! Eq 25.2, pg 163
!      LO = 280.4664567+36000.76982779*t+0.0003032028*(t**2)+(t**3)/49931-(t**4)/15300-(t**5)/2000000 ! Eq 28.2, pg 183, difference with line above is negligible
      DO WHILE(LO > 360 .or. LO < 0)
       If(LO > 360.0) LO = LO - 360
       If(LO < 0.0)   LO = LO + 360
      END DO
! calculate the Geometric Mean Anomaly of the Sun, degrees
      M = 357.52911 + t*(35999.05029-0.0001537*t)                                              ! Eq 25.3, pg 163
! calculate the eccentricity of earth's orbit
      e = 0.016708634 - t*(0.000042037+0.0000001267*t)                                         ! Eq 25.4, pg 163
! calculate the equation of center for the sun, degrees
      Mrad = DToR*M
      sinM = Sin(Mrad)
      sin2M = Sin(2*Mrad)
      sin3M = Sin(3*Mrad)
      c = sinM*(1.914602-t*(0.004817+0.000014*t))+sin2M*(0.019993-0.000101*t)+sin3M*0.000289   ! Between Eq 25.4 and Eq 25.5, pg 164
! calculate the true longitude of the sun, degrees
      TLO = LO + c                                                                             ! Between Eq 25.4 and Eq 25.5, pg 164
! calculate the true anomaly of the sun, degrees
      v = M + c                                                                                ! Between Eq 25.4 and Eq 25.5, pg 164
! calculate the distance to the sun in AU
      r = (1.000001018*(1 - e*e))/(1 + e*Cos(DToR*v))                                          ! Eq 25.5, pg 164
! calculate the apparent longitude of the sun, degrees
      omega = 125.04 - 1934.136*t                                                              ! Between Eq 25.5 and Eq 25.6, pg 164
!      omega  = 125.04452 - 1934.136261*t+0.0020708*(t**2)+(t**3)/45000                        ! Equation on page 144, difference with line above is negligible
      lambda = TLO - 0.00569 - 0.00478*Sin(DToR*omega)                                         ! Between Eq 25.5 and Eq 25.6, pg 164
! calculate the mean obliquity of the ecliptic, degrees
      seconds = 21.448 - t*(46.8150 + t*(0.00059 - 0.001813*t))                                ! Eq 22.2, pg 147
      e0 = 23.0 + (26.0 + (seconds/60.0))/60.0                                                 ! Eq 22.2, pg 147
! calculate the corrected obliquity of the ecliptic, degrees
      ep = e0 + 0.00256 * Cos(DToR*omega)
!      ep = e0 + (9.20/3600.0)*Cos(DToR*omega)+(0.57/3600.0)*cos(DToR*2*LO)                    ! Eq on page 144, still truncated though here in code, difference with line above is negligible

! earth heliocentric latitude, BB, degrees, Equations 9 to 12, NREL/TP-560-34302
      B0   = 280.0*Cos(3.199+84334.662*jme)+102.0*Cos(5.422+5507.553*jme)+  &
              80.0*Cos(3.880+ 5223.690*jme)+ 44.0*Cos(3.700+2352.870*jme)+  &
              32.0*Cos(4.000+ 1577.340*jme)
      B1   =   9.0*Cos(3.900+ 5507.550*jme)+  6.0*Cos(1.730+5223.690*jme)
      BB   = ((B0+B1*jme)/100000000.0)*RToD
! geocentric latitude, Beta, degrees
      Beta = -1.0*BB
! calculate the right ascension of the sun, degrees
      tananum = Cos(DToR*ep) * Sin(DToR*lambda)
!      tananum   = Cos(DToR*ep)*Sin(DToR*lambda)-tan(DToR*Beta)*Sin(DToR*ep)                   ! Eq 13.3 on page 93, difference with line above is negligible
      tanadenom = Cos(DToR*lambda)                                                             ! Eq 13.3 on page 93
! Javascript and FORTRAN use same format, Atan2(Y,X)
      alpha = RtoD*atan2(tananum, tanadenom)                                                   ! Eq 13.3 on page 93

! calculate the declination of the sun, degrees
      sint = Sin(DToR*ep) * Sin(DToR*lambda)
!      sint = Sin(DToR*ep)*Sin(DToR*lambda)*Cos(DToR*Beta)+Cos(DToR*ep)*Sin(DToR*Beta)         ! Eq 13.4 on page 93, difference with line above is negligible
      theta = RToD*Asin(sint)                                                                  ! Eq 13.4 on page 93

! calculate the equation of time, the difference between true solar time
! and mean solar time, in minutes of time
      y = (Tan(DToR*ep/2.0))**2
      sin2LO = Sin(2.0*DToR*LO)
      cos2LO = Cos(2.0*DToR*LO)
      sin4LO = Sin(4.0*DToR*LO)
      Etime = y*sin2LO - 2.0*e*sinM + 4.0*e*y*sinM*cos2LO - 0.5*y*y*sin4LO - 1.25*e*e*sin2M    ! Eq 28.3 on page 185
      EQT = RToD*Etime*4.0

! change sign convention for longitude from negative to positive in western hemisphere
      if(Long(i) < 0.0)then
       Long2(i) = Long(i)*-1.0
      else
       Long2(i) = Long(i)
      end if
      if(Lat(i) > 89.8)  Lat(i) = 89.8
      if(Lat(i) < -89.8) Lat(i) = -89.8

      solarTimeFix = EQT - 4.0*Long2(i) !+ 60.0*TZ(i)  ! in minutes
      trueSolarTime = HOUR*60.0 + solarTimeFix   !in minutes

      if(trueSolarTime > 1440)then
       trueSolarTime = trueSolarTime - 1440
      end if

      hourangle = trueSolarTime/4.0 - 180.0

!Thanks to Louis Schwarzmayr for the next line:
      if(hourangle < -180.0)hourangle = hourangle + 360.0

!Solar zenith (Solar altitude)
      harad = DToR*hourangle
      csz = Sin(DToR*Lat(i))*Sin(DToR*theta) + Cos(DToR*Lat(i))*Cos(DToR*theta)*Cos(harad)
      if(csz > 1.0) Then
       csz = 1.0
      else if(csz < -1.0) Then
       csz = -1.0
      end if
      zenith = RToD*Acos(csz)

!Solar azimuth
      azDenom = Cos(DToR*Lat(i))*Sin(DToR*zenith)
      if(Abs(azDenom) > 0.001)then
       azRad = ((Sin(DToR*Lat(i))*Cos(DToR*zenith)) - Sin(DToR*theta))/azDenom
       if(Abs(azRad) > 1.0)then
        if(azRad < 0.0)then
         azRad = -1.0
        else
         azRad = 1.0
        end if
       end if
       azimuth = 180.0 - RToD*Acos(azRad)
       if(hourangle > 0.0)then
        azimuth = azimuth*-1.0
       end if
      else
       if (Lat(i) > 0.0)then
        azimuth = 180.0
       else
        azimuth = 0.0
       end if
      end if

      if(azimuth < 0.0)then
       azimuth = azimuth + 360.0
      end if

! Correction for the effect of atmospheric refraction, NOAA website (SElev is Solar Altitude)
      SElev = 90.0 - zenith
      if(SElev > 85.0)then
       RC = 0.0
      else
       Te = Tan(DToR*SElev)
       if(SElev > 5.0 .and. SElev <= 85.0)   RC = (58.1/Te - 0.07/(Te**3) + 0.000086/(Te**5))/3600.0
       if(SElev > -0.575 .and. SElev <= 5.0) RC = (1735-518.2*SElev+103.4*SElev**2-12.79*SElev**3+0.711*SElev**4)/3600.0
       if(SElev < -0.575)                    RC = (-20.774/Te)/3600.0
      end if

       solarzen     = zenith - RC                    ! solar zenith, corrected, degrees
       solarazimuth = azimuth                        ! solar azimuth, degrees
       solaralt     = 90.0 - solarzen                ! solar altitude, degrees
       solarzen     = solarzen*DToR                  ! solar zenith, radians

       !Call Ozone(i,k,Jdaym)
          if(Lat(i) >= 0.0)then
      A = 150.0; C = 40.0; F = -30.000; H = 3.0; B = 1.28
      if(Long(i) >= 0.0)then
       P = 20.0
      else
       P =  0.0
      end if
    else
      A = 100.0; C = 30.0; F = 152.625; H = 2.0; B = 1.50; P = -75.0
    end if
! D converts the number of days to a fraction part of 360 degrees, in degrees
!    D  = 360.0/365.25, 0.98562628
! Uo is calculated in atm-cm
    Uo  = (235+(A+C*Sin(DToR*0.9856*(int(Jdaym)+F))+20.0*Sin(DToR*H*(Long(i)+P)))*(sin(DToR*B*Lat(i))**2))/1000.0
       !Call BirdModel(i,k,Jdaym)
       ! Atmospheric Pressure
     P         = 1013
! Ratio of forward-scattered irradiance to the total scattered, see SERI pages 8 & 9
!     Ba        = 0.84
! Constant used in Bird Model associated with aerosol absorptance, see SERI pages 7 & 9
!     K1        = 0.10
! Used same albdeo from other code for Methods 1 to 3
!     albedo    = 0.10
! NREL solar constant of 1367 W/m^2
!     SROCON    = 1367.00
! Bras/TVA (Io=Wo(sin(solaralt))/r2, Equation 2.9 on page 26) adjustment of ETR normal to beam for earth sun distance
     birdETR   = SROCON/r**2
! broadband aerosol optical depth from surface in a vertical path (broadband turbidity)
     tauA      = 0.2758*AOD380nm + 0.35*AOD500nm
!     tauA      = 0.05
    if(RToD*solarzen < 89.0)then
! Relative Optical Air Mass, SERI, pg 8
     AM      = (((288-0.0065*Elev(i))/288)**5.256)/(Sin(DToR*solaralt)+0.1500*(solaralt+3.885)**-1.253)                    ! Equation 23, optical air mass
! Pressure corrected air mass
     AMP       = AM*P/1013
! Total amount of ozone in a slanted path (cm)
     Xo        = Uo*AM
! Amount of Precipitable water in a vertical column  from surface (cm), used 1.50 cm for Mt Vernon, range from MLS and USS models 1.42 to 2.93
     Uw        = exp(-0.0592+0.06912*DP(i,k))                                                                                ! Equation 27, mean hourly precipitable water content, cm
! Total amount of precipitable water in a slanted path (cm)
     Xw        = Uw*AM
! Transmittance of Rayleigh Scattering, SERI, pg 8
     TR        = Exp(-0.0903*(AMP**0.84)*(1+AMP-AMP**1.01))
! Transmittance of ozone absorptance, SERI, pg 8
     Toh       = 1-0.1611*Xo*(1+139.48*Xo)**-0.3035-0.002715*Xo/(1+0.044*Xo+0.0003*Xo**2)
! Transmittance of absorptance of uniformly mixed gases (carbon dioxide and oxygen)
     TUM       = Exp(-0.0127*AMP**0.26)
! Transmittance of water vapor absoptance (1-aw)
     Tw        = 1-2.4959*Xw/((1 + 79.034*Xw)**0.6828+6.385*Xw)
! Transmittance of aerosol aborptance and scattering
     TA        = Exp(-(tauA**0.873)*(1+tauA-tauA**0.7088)*AM**0.9108)
! Transmittance of aerosol aborptance
     TAA       = 1-K1*(1-AM+AM**1.06)*(1-TA)
! Sky, or atmospheric, albedo
     rs        = 0.0685+(1-Ba)*(1-(TA/TAA))
! Direct solar irradiance on a horizontal ground surface (W/m^2)
     if(RToD*solarzen < 90.0)then
      directHz = birdETR*Cos(solarzen)*0.9662*TA*Tw*TUM*Toh*TR
     else
      directHz = 0.0
     end if
! Solar irradiance on a horizontal surface from atmospheric scattering (W/m2)
     numI      = 0.5*(1-TR)+ Ba * (1-(TA/TAA))
     denI      = 1-AM+AM**1.02
     Ias       = birdETR*Cos(solarzen)*0.79*Toh*TUM*Tw*TAA*(numI/DenI)
   else
     AM = 0.0; TR = 0.0; Toh = 0.0; TUM = 0.0; Tw = 0.0; TA = 0.0; TAA = 0.0; rs = 0.0; Id = 0.0; Ias = 0.0; AMP = 0.0
     Xw = 0.0; Xo = 0.0; mp = 0.0
   end if

!global radiation on a horizontal ground surface (W/m^2)
   if(AM > 0.0)then
    globalHz   = (directHz+Ias)/(1-albedo*rs)
   else
    globalHz   = 0.0
   end if

!diffuse radiation on a horizontal ground surface (W/m^2)
   diffuseHz = globalHz - directHz
       SRON4(i,k)   = globalHz
       TZ(i)= TZ(i)*-1.0

  ELSE
  LOCAL    =  LONGIT(JW)
  STANDARD =  15.0*INT(LONGIT(JW)/15.0)
  HOUR     = (JDAY-INT(JDAY))*24.0
  IDAY     =  JDAY-((INT(JDAY/365))*365)
  IDAY     =  IDAY+INT(INT(JDAY/365)/4)
  TAUD     = (2.*PI*(IDAY-1))/365.
  EQTNEW   =  0.170*SIN(4.*PI*(IDAY-80)/373.)-0.129*SIN(2.*PI*(IDAY-8)/355.)
  HH(JW)   =  0.261799*(HOUR-(LOCAL-STANDARD)*0.0666667+EQTNEW-12.0)
  DECL(JW) =  0.006918-0.399912*COS(TAUD)+0.070257*SIN(TAUD)-0.006758*COS(2.*TAUD)+0.000907*SIN(2.*TAUD)-0.002697*COS(3.*TAUD)        &
              +0.001480*SIN(3.*TAUD)
  SINAL    =  SIN(LAT(JW)*.0174533)*SIN(DECL(JW))+COS(LAT(JW)*.0174533)*COS(DECL(JW))*COS(HH(JW))
  A00(JW)  =  57.2957795*ASIN(SINAL)
  A0       =  A00(JW)
  IF (A0 > 0.0) THEN
    SRON(JW) = (1.0-0.0065*CLOUD(JW)**2)*24.0*(2.044*A0+0.1296*A0**2-1.941E-3*A0**3+7.591E-6*A0**4)*BTU_FT2_DAY_TO_W_M2
  ELSE
    SRON(JW) = 0.0
  END IF
  ENDIF
RETURN

!***********************************************************************************************************************************
!**                                          E Q U I L I B R I U M  T E M P E R A T U R E                                         **
!***********************************************************************************************************************************

ENTRY EQUILIBRIUM_TEMPERATURE

! British units

  TDEW_F   = DEG_F(TDEW(JW))
  TAIR_F   = DEG_F(TAIR(JW))
  SRO_BR   = SRON(JW)*W_M2_TO_BTU_FT2_DAY*SHADE(I)
  WIND_MPH = WIND(JW)*WSC(I)*MPS_TO_MPH
  WIND2M   = WIND_MPH*DLOG(2.0D0/Z0(JW))/DLOG(WINDH(JW)/Z0(JW))+NONZERO     ! SW 11/28/07  old version z0=0.003
  ACONV    = W_M2_TO_BTU_FT2_DAY
  IF (CFW(JW) == 1.0) BCONV = 3.401062
  IF (CFW(JW) == 2.0) BCONV = 1.520411

! Equilibrium temperature and heat exchange coefficient

  ET(I)   =  TDEW_F
  TSTAR   = (ET(I)+TDEW_F)*0.5
  BETA    =  0.255-(8.5E-3*TSTAR)+(2.04E-4*TSTAR*TSTAR)
  FW      =  ACONV*AFW(JW)+BCONV*BFW(JW)*WIND2M**CFW(JW)
  CSHE(I) =  15.7+(0.26+BETA)*FW
  RA      =  3.1872E-08*(TAIR_F+459.67)**4
  ETP     = (SRO_BR+RA-1801.0)/CSHE(I)+(CSHE(I)-15.7)*(0.26*TAIR_F+BETA*TDEW_F)/(CSHE(I)*(0.26+BETA))
  J       =  0
  DO WHILE (ABS(ETP-ET(I)) > 0.05 .AND. J < 10)
    ET(I)   =  ETP
    TSTAR   = (ET(I)+TDEW_F)*0.5
    BETA    =  0.255-(8.5E-3*TSTAR)+(2.04E-4*TSTAR*TSTAR)
    CSHE(I) =  15.7+(0.26+BETA)*FW
    ETP     = (SRO_BR+RA-1801.0)/CSHE(I)+(CSHE(I)-15.7)*(0.26*TAIR_F+BETA*TDEW_F)/(CSHE(I)*(0.26+BETA))
    J       =  J+1
  END DO

! SI units

  ET(I)   = DEG_C(ET(I))
  CSHE(I) = CSHE(I)*FLUX_BR_TO_FLUX_SI/RHOWCP
RETURN

!***********************************************************************************************************************************
!**                                                   S U R F A C E   T E R M S                                                   **
!***********************************************************************************************************************************

ENTRY SURFACE_TERMS (TSUR)

! Partial water vapor pressure of air (mm hg)

!  EA = EXP(2.3026*(9.5*TDEW(JW)/(TDEW(JW)+265.5)+0.6609))         ! SW 6/10/2011
!  IF (TDEW(JW) > 0.0) EA = EXP(2.3026*(7.5*TDEW(JW)/(TDEW(JW)+237.3)+0.6609))
  EA = DEXP(2.3026D0*(7.5D0*TDEW(JW)/(TDEW(JW)+237.3D0)+0.6609D0))

! Partial water vapor pressure at the water surface

  IF(TSUR<0.0)THEN
  ES = DEXP(2.3026D0*(9.5D0*TSUR/(TSUR+265.5D0)+0.6609D0))
  ELSE
  ES = DEXP(2.3026D0*(7.5D0*TSUR/(TSUR+237.3D0)+0.6609D0))
  ENDIF

! Wind function

  IF (RH_EVAP(JW)) THEN
    TAIRV = (TAIR(JW)+273.0D0)/(1.0D0-0.378D0*EA/760.0D0)
    DTV   = (TSUR+273.0D0)/(1.0D0-0.378D0*ES/760.0D0)-TAIRV
    DTVL  =  0.0084D0*WIND2(I)**3
    IF (DTV < DTVL) DTV = DTVL
    FW = (3.59D0*DTV**0.3333D0+4.26D0*WIND2(I))
  ELSE
    FW = AFW(JW)+BFW(JW)*WIND2(I)**CFW(JW)
  END IF

! Evaporative flux

  RE(I) = FW*(ES-EA)
  !IF(RE(I) < 0.0)RE(I)=0.0     ! SW 6/22/2016  SHOULD WE USE THIS AS AN ANALOG FOR CONDENSATION WHEN LESS THAN 0? TVA(1972) SUGGESTS YOU CAN - SEE SECTION 4 AND 5 but Ryan, Harleman suggest setting RE=0 if less than zero since the condensation coefficients are unknown...

! Conductive flux

  RC(I) = FW*BOWEN_CONSTANT*(TSUR-TAIR(JW))

! Back radiation flux

  RB(I) = 5.51D-8*(TSUR+273.15D0)**4
  RETURN
END SUBROUTINE HEAT_EXCHANGE

! Function declaration

   REAL(R8) FUNCTION DEG_F(X)
   USE PREC
   REAL(R8) :: X
   DEG_F =  X*1.8+32.0
   END FUNCTION DEG_F
   REAL(R8) FUNCTION DEG_C(X)
   USE PREC
   REAL(R8) :: X
   DEG_C = (X-32.0)*5.0/9.0
   END FUNCTION DEG_C