text:
Minutes from Senate meetings are typically available in SGAO and at http://asucd.ucdavis.edu asucd.ucdavis.edu.  Any questions/concerns with minutes should be forwarded to SGAO or to the current ASUCD Recorder.  Officers reports and meeting notes will be put on file in SGAO on the third floor of the MU without the actual minutes until they have been completed. Minutes from the previous weeks meeting are approved every Senate meeting and are found at http://asucd.ucdavis.edu/resources/governmentdocuments/minutes/ this website shortly after.

Minutes are downloadable and are stored in the .doc format, best used with Microsoft Word or equivalent software.

A recent shift in the style of the minutes is due to Jasmine Chen, who had been recorder for close to three years, recently finishing her term alongside President Rebecca Sterling at the end of the 2013 Winter Quarter.

The California Aggie writes Senate Briefs in their weekly edition. They are notable in their interesting prioritization of which details make print and tone, as seen in http://www.theaggie.org/2013/05/29/asucdsenatebriefs5/ this May 2013 example.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


text:
      program dipole

c*********************************************************************
c
c     dl_poly program to calculate cluster dipoles
c     for selected atom clusters from dl_poly HISTORY file
c
c     copyright daresbury laboratory 1996
c     author  w.smith oct 1996
c
c*********************************************************************

      implicit real*8(a-h,o-z)

      parameter (mxatms=1080,mxcls=2)
      character*1 fmt
      character*40 fname
      character*80 title
      character*8 atmnam(mxatms)
      dimension xxx(mxatms),yyy(mxatms),zzz(mxatms)
      dimension weight(mxatms),chge(mxatms),cell(9)
      dimension cmass(3,mxcls),dipole(3,mxcls)
      dimension ibc(mxcls),iec(mxcls)

c     open the I/O files

      open(5,file='dip_input')
      open(6,file='dip_output')

      write(6,'(a)')
     x  '# Cluster Dipole Program'

c     read the control variables

c     number of clusters
      read(5,*)ncls

c     specify atom clusters

      do i=1,ncls
        read(5,*)ibc(i),iec(i)
      enddo

c     name of selected HISTORY file
      read(5,*)fname

c     HISTORY file formatted or unformatted
      read(5,*)fmt

c     total number of atoms in a HISTORY file configuration
      read(5,*)natms

c     max number of configurations to sample
      read(5,*)nconf

c     read time interval between configurations
      read(5,*)tstep

c     check on specified control variables

      write(6,'(a,a40)')'# name of target HISTORY file   : ',fname
      write(6,'(a,i8)')'# number of atomic clusters     : ',ncls
      do i=1,ncls
        write(6,'(a,2i8)')'# start and end atoms   : ',ibc(i),iec(i)
      enddo
      write(6,'(a,i8)')'# total no. of atoms in config  : ',natms
      write(6,'(a,i8)')'# number of configurations      : ',nconf
      write(6,'(a,f12.6)')'# time interval between configs : ',tstep

      if(natms.gt.mxatms)then

        write(6,'(a,2i5)')
     x    '# error - too many atoms in system'
        stop

      endif

c     set default cell properties

      do i=1,9
        cell(i)=0.d0
      enddo
      cell(1)=1.d0
      cell(5)=1.d0
      cell(9)=1.d0

      iflg=0
      keytrj=0

      do iconf=0,nconf-1

        if(fmt.eq."U".or.fmt.eq."u")then

          call uread
     x      (fname,title,atmnam,iflg,imcon,keytrj,natms,nstep,timstp,
     x       cell,chge,weight,xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz)

        else

          call hread
     x      (fname,title,atmnam,iflg,imcon,keytrj,natms,nstep,timstp,
     x       cell,chge,weight,xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz)

        endif

        if(imcon.ne.0)then

          write(6,*)'error - periodic boundary condition in operation'
          stop

        endif

        if(iflg.lt.0)go to 100

c     initialise dipole and centre of mass arrays

        do i=1,mxcls

          cmass(1,i)=0.d0
          cmass(2,i)=0.d0
          cmass(3,i)=0.d0
          dipole(1,i)=0.d0
          dipole(2,i)=0.d0
          dipole(3,i)=0.d0

        enddo

        do j=1,ncls

          huge=0.d0

          do i=ibc(j),iec(j)

            huge=huge+weight(i)

            cmass(1,j)=weight(i)*xxx(i)+cmass(1,j)
            cmass(2,j)=weight(i)*yyy(i)+cmass(2,j)
            cmass(3,j)=weight(i)*zzz(i)+cmass(3,j)

          enddo

          cmass(1,j)=cmass(1,j)/huge
          cmass(2,j)=cmass(2,j)/huge
          cmass(3,j)=cmass(3,j)/huge

          do i=ibc(j),iec(j)

            dipole(1,j)=chge(i)*(xxx(i)-cmass(1,j))+dipole(1,j)
            dipole(2,j)=chge(i)*(yyy(i)-cmass(2,j))+dipole(2,j)
            dipole(3,j)=chge(i)*(zzz(i)-cmass(3,j))+dipole(3,j)

          enddo

        enddo

        time=tstep*dble(iconf)
        write(6,'(1p,10e12.4)')time,dipole

      enddo

  100 continue

      stop

      end
      subroutine uread
     x  (history,cfgname,atmnam,iflg,imcon,keytrj,natms,nstep,tstep,
     x   cell,chge,weight,xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz)

c
c***********************************************************************
c
c     dl_poly subroutine for reading unformatted history files
c
c     double precision, single processor version
c
c     copyright - daresbury laboratory 1996
c     author    - w. smith jan 1996.
c
c***********************************************************************
c

      implicit real*8(a-h,o-z)

      logical new

      character*80 cfgname
      character*40 history
      character*8 atmnam(*)

      dimension cell(9)
      dimension chge(*),weight(*)
      dimension xxx(*),yyy(*),zzz(*)
      dimension vxx(*),vyy(*),vzz(*)
      dimension fxx(*),fyy(*),fzz(*)

      save new

      data new/.true./,nhist/77/

      iflg=0

c     open the history file if new job

      if(new)then

        open(nhist,file=history,form='unformatted',
     x       status='old',err=100)

        read(nhist,err=200) cfgname
        write(*,'(a,a)')'# History file header: ',cfgname
        read(nhist,end=200) datms
        if(natms.ne.nint(datms))then

          matms=nint(datms)
          write(*,'(a)')'# error - incorrect number of atoms in file'
          write(*,'(a,i6,a)')'# file contains',matms,' atoms'
          stop

        endif
        read(nhist,end=200) (atmnam(i),i=1,natms)
        read(nhist,end=200) (weight(i),i=1,natms)
        read(nhist,end=200) (chge(i),i=1,natms)

      endif

      read(nhist,end=200)dstep,datms,trjkey,dimcon
      nstep=nint(dstep)
      ktrj=nint(trjkey)
      imcon=nint(dimcon)
      if(keytrj.gt.ktrj)then

        if(ktrj.eq.0)write(*,'(a)')'# error - no velocities in file'
        if(keytrj.gt.1)write(*,'(a)')'# error - no forces in file'
        stop

      endif

      if(imcon.gt.0) read(nhist,end=200) cell

      read(nhist,end=200) (xxx(i),i = 1,natms)
      read(nhist,end=200) (yyy(i),i = 1,natms)
      read(nhist,end=200) (zzz(i),i = 1,natms)

      if(keytrj.ge.1)then
        read(nhist,end=200) (vxx(i),i = 1,natms)
        read(nhist,end=200) (vyy(i),i = 1,natms)
        read(nhist,end=200) (vzz(i),i = 1,natms)
      else if(ktrj.ge.1)then
        read(nhist,end=200)
        read(nhist,end=200)
        read(nhist,end=200)
      endif
      if(keytrj.ge.2)then
        read(nhist,end=200) (fxx(i),i = 1,natms)
        read(nhist,end=200) (fyy(i),i = 1,natms)
        read(nhist,end=200) (fzz(i),i = 1,natms)
      else if(ktrj.ge.2)then
        read(nhist,end=200)
        read(nhist,end=200)
        read(nhist,end=200)
      endif

      iflg=1

      return

  100 continue

      write(*,'(a)')'# error - History file not found'
      stop

  200 continue
      write(*,'(a)')'# warning - end of History file encountered'
      close (nhist)
      iflg=-1

      return
      end
      subroutine hread
     x  (history,cfgname,atmnam,iflg,imcon,keytrj,natms,nstep,tstep,
     x   cell,chge,weight,xxx,yyy,zzz,vxx,vyy,vzz,fxx,fyy,fzz)

c
c***********************************************************************
c
c     dl_poly subroutine for reading the formatted history file
c
c     copyright - daresbury laboratory 1996
c     author    - w. smith jan 1996.
c
c     single processor version
c
c***********************************************************************
c

      implicit real*8(a-h,o-z)

      logical new

      character*80 cfgname
      character*40 history
      character*8 atmnam(*),step

      dimension cell(9)
      dimension chge(*),weight(*)
      dimension xxx(*),yyy(*),zzz(*)
      dimension vxx(*),vyy(*),vzz(*)
      dimension fxx(*),fyy(*),fzz(*)

      save new

      data new/.true./,nhist/77/

      iflg=0

c     open history file if new job

      if(new)then

        open(nhist,file=history,status='old',err=100)

        read(nhist,'(a80)',err=200) cfgname
        write(*,'(a,a)')'# History file header: ',cfgname
        read(nhist,'(2i10)',end=200) ktrj,imcon
        if(keytrj.gt.ktrj)then

          if(ktrj.eq.0)write(*,'(a)')'# error - no velocities in file'
          if(keytrj.gt.1)write(*,'(a)')'# error - no forces in file'
          stop

        endif

        new=.false.

      endif

      read(nhist,*,end=200)
     x     step,nstep,matms,ktrj,imcon,tstep

      if(natms.ne.matms)then

        write(*,'(a)')'# error - incorrect number of atoms in file'
        write(*,'(a,i6,a)')'# file contains',matms,' atoms'
        stop

      endif

      if(imcon.gt.0) read(nhist,'(3g12.4)',end=200) cell

      do i = 1,natms
        read(nhist,'(a8,i10,2f12.6)',end=200)
     x    atmnam(i),j,weight(i),chge(i)
        read(nhist,'(1p,3e12.4)',end=200) xxx(i),yyy(i),zzz(i)
        if(keytrj.ge.1)then
          read(nhist,'(1p,3e12.4)',end=200) vxx(i),vyy(i),vzz(i)
        else if(ktrj.ge.1)then
          read(nhist,'(1p,3e12.4)',end=200) vx,vy,vz
        endif
        if(keytrj.ge.2)then
          read(nhist,'(1p,3e12.4)',end=200) fxx(i),fyy(i),fzz(i)
        else if(ktrj.ge.2)then
          read(nhist,'(1p,3e12.4)',end=200) fx,fy,fz
        endif
      enddo

      iflg=1

      return

  100 continue

      write(*,'(a)')'# error - History file not found'
      stop

  200 continue
      write(*,'(a)')'# warning - end of History file encountered'
      close (nhist)
      iflg=-1

      return
      end


text:
C     Check that a new module can be added to a workspace
C     when it is named in tpips

      program add_module00

      end


text:
$TB_SRC_PATH/os-vvm/SortListPkg_int.vhd \
$TB_SRC_PATH/os-vvm/RandomBasePkg.vhd \
$TB_SRC_PATH/os-vvm/RandomPkg.vhd \
$TB_SRC_PATH/os-vvm/CoveragePkg.vhd \

max_stars_repo_name: synvue/generic-fifo

max_issues_repo_name: synvue/generic-fifo

max_forks_repo_name: synvue/generic-fifo


text:
subroutine csp_gen(mass_ssp, lbol_ssp, spec_ssp, &
     pset, tage, nzin, mass_csp, lbol_csp, spec_csp, &
     mdust_csp,emlin_ssp,emlin_csp)
  !
  ! Return the spectrum (and mass and lbol) of a composite stellar population.
  !
  ! Inputs
  ! ---------
  !
  ! mass_ssp, lbol_ssp, spec_ssp:
  !   The (surviving) stellar masses, bolometric luminosities, and spectra of
  !   the SSPs.
  !
  ! pset:
  !   A `PARAMS` structure containing the SFH parameters
  !
  ! tage:
  !   The age (in Gyr of forward time) at which the spectrum is desired.  Note
  !   that this can be different than pset%tage if the latter is 0.
  !
  ! Outputs
  ! ---------
  !
  ! mass_csp, lbol_csp, spec_csp:
  !   The (surviving) stellar masses, bolometric luminosity, and spectrum of
  !   the composite stellar population at tage, normalized to 1 M_sun *formed*.

  use sps_vars, only: ntfull, nspec, time_full, tiny_number, tiny_logt, &
                      zlegend, nz, sfh_tab, ntabsfh, compute_light_ages, &
                      SFHPARAMS, PARAMS, SP, nemline, dust_type
  use sps_utils, only: locate, sfh_weight, sfhinfo, add_dust
  implicit none

  real(SP), intent(in), dimension(ntfull, nzin) :: mass_ssp, lbol_ssp
  real(SP), intent(in), dimension(nspec, ntfull, nzin) :: spec_ssp
  type(PARAMS), intent(in) :: pset
  real(SP), intent(in) :: tage
  integer, intent(in) :: nzin

  real(SP), intent(out) :: mass_csp, lbol_csp, mdust_csp
  real(SP), intent(out), dimension(nspec) :: spec_csp

  real(SP), DIMENSION(nemline, ntfull, nzin), intent(in) :: emlin_ssp
  real(SP), DIMENSION(nemline), intent(out) :: emlin_csp

  real(SP), dimension(nspec) :: csp1, csp2, lw_age
  real(SP), dimension(nemline) :: ncsp1, ncsp2, nlw_age
  real(SP), dimension(ntfull, nzin) :: total_weights
  real(SP), dimension(ntfull) :: w1=0., w2=0.
  integer :: i, j, k, imin, imax, i_tesc
  type(SFHPARAMS) :: sfhpars
  real(SP) :: m1, m2, frac_linear, mfrac, sfr, fburst
  real(SP) :: t1, t2, dt, zbin, dz  ! for tabular calculations
  real(SP) :: lbol_age, mass_age ! for mass and lbol weighted ages

  ! ------- Setup ----------

  ! Build a structure containing useful units, numbers, and switches for the
  ! weight calculations. The units of the parameters in the `sfhparams`
  ! structure are years of lookback time.
  call convert_sfhparams(pset, tage, sfhpars)
  ! Only calculate SFH weights for SSPs up to tage
  ! (plus the next couple, to bracket and be safe).
  imin = 0
  imax = min(max(locate(time_full, log10(sfhpars%tage)) + 2, 1), ntfull)

  ! ----- Get SFH weights -----

  total_weights = 0.

  ! SSP.
  if (pset%sfh.eq.0) then
     ! Make sure to use SSP weighting scheme
     sfhpars%type = -1
     ! Use tage as the burst lookback time, instead of tage-tburst.
     sfhpars%tb = sfhpars%tage
     ! Only need weights at two SSP points. Though in practice this doesn't
     ! matter, as the appropriate ages are located within sfh_weight.
     ! But it speeds up the matrix multiply later
     imin = max(imax - 2, 1)
     ! These come pre-normalized to 1 Msun
     total_weights(:, 1) = sfh_weight(sfhpars, imin, imax)
  endif

  ! Tau and delayed-tau.
  if ((pset%sfh.eq.1).or.(pset%sfh.eq.4)) then
     imin = 0
     sfhpars%type = pset%sfh
     total_weights(:, 1) = sfh_weight(sfhpars, imin, imax)
     ! Could save some loops by having proper normalization analytically from sfh_weight
     m1 = sum(total_weights(1:imax, 1))
     if (m1.lt.tiny_number) m1 = 1.0
     total_weights(:, 1) = total_weights(:, 1) / m1
  endif

  ! Add constant and burst weights for SFH=1,4
  if (((pset%sfh.eq.1).or.(pset%sfh.eq.4)).and.&
       ((pset%const.gt.0).or.(pset%fburst.gt.tiny_number))) then
     imin = 0
     ! Constant
     sfhpars%type = 0
     w1 = sfh_weight(sfhpars, imin, imax)
     m1 = sum(w1(1:imax))
     ! Burst.  These weights come pre-normalized to 1 Msun.
     ! If burst happens after age of system then we kill it entirely.
     sfhpars%type = -1
     if (sfhpars%tb.lt.0) then
        w2 = 0.
        fburst = 0.
     else
        w2 = sfh_weight(sfhpars, imin, imax)
        fburst = pset%fburst
        ! We'll need to add any early bursts when summing later
        imax = max(imax, min(max(locate(time_full, log10(sfhpars%tb)) + 2, 1), ntfull))
     endif
     ! Sum with proper relative normalization.  Beware divide by zero.
     if (m1.lt.tiny_number) m1 = 1.0
     total_weights(:,1) = (1 - pset%const - fburst) * total_weights(:,1) + &
                          pset%const * (w1 / m1) + &
                          fburst * w2
  endif

  ! Simha
  if (pset%sfh.eq.5) then
     imin = 0
     ! Delayed-tau model portion, could set imin here to be just before
     ! sfhpars%tq, for small speed increase.
     ! imin = min(max(locate(time_full, log10(sfhpars%tq)), 0), ntfull-1)
     sfhpars%type = 4
     w1 = sfh_weight(sfhpars, imin, imax)
     m1 = sum(w1(1:imax))
     ! Linear portion.  Need to set `use_simha_limits` flag to get correct
     ! integration limits. Could set imax here to be just after sfhpars%tq,
     ! but with imin=0, for small speed increase.
     ! imax = imin + 1
     ! imin = 0
     sfhpars%type = 5
     sfhpars%use_simha_limits = 1
     w2 = sfh_weight(sfhpars, imin, imax)
     sfhpars%use_simha_limits = 0
     m2 = sum(w2(1:imax))
     ! Normalize and sum.  need to be careful of divide by zero here, if all
     ! linear or all delay-tau s.t. weights sum to zero for a component.
     if (m1.lt.tiny_number) m1 = 1.0
     if (m2.lt.tiny_number) m2 = 1.0
     ! need to get the fraction of the mass formed in the linear portion
     call sfhinfo(pset, tage, mfrac, sfr, frac_linear)
     total_weights(:, 1) = (w1 / m1) * (1 - frac_linear) + frac_linear * (w2 / m2)
     ! imax = min(max(locate(time_full, log10(sfhpars%tage)) + 2, 1), ntfull)
  endif

  ! Tabular.  Time units in sfh_tab are assumed to be linear years of time
  ! since big bang (forward time).  We are going to treat this as a sum of
  ! linear SFHs, one for each bin in the table
  if (pset%sfh.eq.2.or.pset%sfh.eq.3) then
     total_weights = 0.

     ! Assume linear SFH within the bins
     sfhpars%type = 5
     ! Loop over each bin in the table.
     do j=1, ntabsfh-1
        ! Edges of the bin in lookback time. Note that the order of sfhtab gets
        ! flipped, since it is given in forward time and then we convert to
        ! lookback time.  So j=0 is the `oldest` in terms of lookback time
        t1 = tage*1e9 - sfh_tab(1, j+1)
        t2 = tage*1e9 - sfh_tab(1, j)
        if (t2.lt.0) then
           ! Entire bin is in the future, skip
           cycle
        endif
        ! Metallicity of the bin. Just a straight average.
        zbin = (sfh_tab(3, j) + sfh_tab(3, j+1)) / 2

        ! Linear slope.  Positive should be sfr *decreasing* in time since big bang.
        sfhpars%sf_slope = -(sfh_tab(2, j+1) - sfh_tab(2, j)) / (t2 - t1) / sfh_tab(2, j+1)
        ! Set integration limits using bin edges clipped to valid times.
        ! That is, don't include any portion of a bin that goes to negative
        ! time, or beyond the oldest isochrone.
        sfhpars%tq = min(max(t1, 10**tiny_logt), 10**time_full(ntfull))  ! lower limit (in lookback time)
        sfhpars%tage = min(max(t2, 10**tiny_logt), 10**time_full(ntfull)) ! upper limit
        sfhpars%sf_trunc = sfhpars%tage - sfhpars%tq
        ! Mass that formed within these valid times.
        dt = (sfhpars%tage - sfhpars%tq)
        m2 = sfh_tab(2, j+1) * (1 + sfhpars%sf_slope/2. * (sfhpars%tage + sfhpars%tq - 2*t1)) * dt
        ! min and max ssps to consider, being conservative.
        imin = min(max(locate(time_full, log10(t1)) - 1, 0), ntfull)
        imax = min(max(locate(time_full, log10(t2)) + 2, 0), ntfull)

        ! Get the weights for this bin in the tabulated sfh and add to the
        ! total weight, after normalizing.
        w1 = sfh_weight(sfhpars, imin, imax)
        m1 = sum(w1)
        if (m1.lt.tiny_number) m1 = 1.0
        ! This is where we'd assign to specific metallicities, if taking that
        ! into account.  This scheme assumes entire bin is at average of the
        ! two enclosing Z values.
        if (nzin.gt.1) then
           k = max(min(locate(zlegend, zbin), nz-1), 1)
           dz = (log10(zbin) - log10(zlegend(k))) / &
                (log10(zlegend(k+1)) - log10(zlegend(k)))
           dz = max(min(dz, 1.0), -1.0) !don't extrapolate
           total_weights(:, k) = total_weights(:, k) + (1-dz) * w1 * (m2 / m1)
           total_weights(:, k+1) = total_weights(:, k+1) + dz * w1 * (m2 / m1)
        else
           total_weights(:, 1) = total_weights(:, 1) + w1 * (m2 / m1)
        endif
     enddo
     ! Reset imin and imax for the spectral sum.
     imin = 0
     imax = ntfull
  endif

  ! ----- Combine SSPs with dust -------
  ! Now weight each SSP by `total_weight`, assign to young or old, and feed to
  ! add_dust, which does the final sum as well as attenuating.
  ! This matrix multiply could probably be optimized!!!!

  csp1  = 0. ! young
  csp2  = 0. ! old
  ncsp1 = 0. ! young nebular
  ncsp2 = 0. ! old nebular
  if (compute_light_ages.eq.1) then
     nlw_age  = 0.
     lw_age   = 0.
     lbol_age = 0.
     mass_age = 0.
  endif

  ! Dust treatment is not strictly correct, since the age bin older than
  ! dust_tesc will include some contribution from young star dust due to the
  ! interpolation, and changing dust_tesc by values smaller than half the ssp
  ! age grid resolution will have no effect on the output.
  i_tesc = locate(time_full, pset%dust_tesc)
  do i=max(imin, 1), imax
     do k=1,nzin
        if (total_weights(i, k).gt.tiny_number) then
           if (i.le.i_tesc) then
              csp1  = csp1  + total_weights(i, k) * spec_ssp(:, i, k)
              ncsp1 = ncsp1 + total_weights(i, k) * emlin_ssp(:, i, k)
           else
              csp2  = csp2  + total_weights(i, k) * spec_ssp(:, i, k)
              ncsp2 = ncsp2 + total_weights(i, k) * emlin_ssp(:, i, k)
           endif
           ! Now do numerator in case of light and mass weighted ages.
           if (compute_light_ages.eq.1) then
              nlw_age = nlw_age + total_weights(i, k) * emlin_ssp(:, i, k) &
                       * 10**(time_full(i)-9)
              lw_age = lw_age + total_weights(i, k) * spec_ssp(:, i, k) &
                       * 10**(time_full(i)-9)
              lbol_age = lbol_age + 10**lbol_ssp(i, k) * total_weights(i, k) &
                         * 10**(time_full(i)-9)
              mass_age = mass_age + mass_ssp(i, k) * total_weights(i, k) &
                         * 10**(time_full(i)-9)
           endif
        endif
     enddo
  enddo
  mass_csp = sum(mass_ssp * total_weights)
  lbol_csp = log10(sum(10**lbol_ssp * total_weights))

  ! Here we add young and old spectra with dust.
  if (((pset%dust1.gt.tiny_number).or.(pset%dust2.gt.tiny_number).or.(dust_type.eq.3))&
       .and.(compute_light_ages.eq.0)) then
     call add_dust(pset, csp1, csp2, spec_csp, mdust_csp, ncsp1, ncsp2, emlin_csp)

  else
     spec_csp  = csp1 + csp2
     emlin_csp = ncsp1 + ncsp2
     mdust_csp = 0.0
  endif

  if (compute_light_ages.eq.1) then
     spec_csp  = lw_age / spec_csp
     emlin_csp = nlw_age / emlin_csp
     lbol_csp  = lbol_age / 10**lbol_csp
     mass_csp  = mass_age / mass_csp
  endif

end subroutine csp_gen

subroutine convert_sfhparams(pset, tage, sfh)
  ! Convert the pset values to yrs, subtract sf_start from the relevant times,
  ! pre-calculate some useful lookback times, and store in a structure.  Note
  ! that this subroutine uses but does not alter the pset values.
  !
  ! Inputs
  ! ----------
  !
  ! pset:
  !    The parameter set.
  !
  ! tage:
  !    The age (in forward time, not lookback time) at which you are
  !    calculating the composite spectrum, in Gyr.
  !
  ! Outputs
  ! ----------
  ! sfh:
  !    An SFHPARAMS structure, containing SFH parameters in yrs, with sf_start
  !    subtracted, and some useful lookback times.
  !       - `tq` is the truncation time, in lookback time
  !       - `t0` is the zero crossing time for a linear SFH, in lookback time.
  !       - `tb` is the burst time, in lookback time.
  !       - `sf_slope` is the fractional change in the SFR in inverse years.  It is
  !          positive for SFR that increases with *lookback* time.
  !
  use sps_vars, only: tiny_number, SFHPARAMS, PARAMS, SP
  implicit none

  type(PARAMS), intent(in) :: pset
  real(SP), intent(in) :: tage

  type(SFHPARAMS), intent(inout) :: sfh

  real(SP) :: start

  ! Define a starting time iff SFH=1,4,5
  if ((pset%sfh.eq.1).or.(pset%sfh.eq.4).or.(pset%sfh.eq.5)) then
     start = pset%sf_start * 1e9
  else
     start = 0.
  endif
  ! Convert units from Gyr to yr, subtract sf_start
  sfh%tage = tage * 1e9 - start
  sfh%tburst = pset%tburst * 1e9 - start
  sfh%sf_trunc = pset%sf_trunc * 1e9 - start
  sfh%tau = pset%tau * 1e9
  ! Note the sign flip here!  pset%sf_slope is positive for sfr that increases
  ! with forward time, sfh%sf_slope is positive for sfr that increases with
  ! *lookback time*.  Yeah, awesome!
  sfh%sf_slope = -pset%sf_slope / 1e9

  ! convert tburst to lookback time
  sfh%tb = sfh%tage - sfh%tburst

  ! convert sf_trunc to to lookback time
  if ((sfh%sf_trunc.le.0).or.(sfh%sf_trunc.gt.sfh%tage)) then
     sfh%tq = 0.
  else
     sfh%tq = sfh%tage - sfh%sf_trunc
  endif
  ! For simha get zero crossing time (in lookback time), avoiding divison by zero
  ! Note that only positive slopes have a chance to hit zero SFR.
  if (sfh%sf_slope.gt.tiny_number) then
     sfh%t0 = sfh%tq - 1. / sfh%sf_slope
  else
     sfh%t0 = 0.
  endif
  ! If the zero crossing time is outside the linear regime, set it to zero.
  if ((sfh%t0.gt.sfh%tq).or.(sfh%t0.le.0)) then
     sfh%t0 = 0.
  endif

end subroutine convert_sfhparams

max_stars_repo_name: christopherlovell/fsps

max_issues_repo_name: christopherlovell/fsps

max_forks_repo_name: christopherlovell/fsps


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

module mod_Kosselmaster
  !! author: MDG
  !! version: 1.0
  !! date: 03/25/20
  !!
  !! class definition for the EMKosselmaster program

use mod_kinds
use mod_global

IMPLICIT NONE

private :: CalcKthick, CalcKint

! namelist for the EMKosselmaster program
type, public :: KosselmasterNameListType
  integer(kind=irg)       :: numthick
  integer(kind=irg)       :: npx
  integer(kind=irg)       :: nthreads
  real(kind=sgl)          :: voltage
  real(kind=sgl)          :: dmin
  real(kind=sgl)          :: startthick
  real(kind=sgl)          :: thickinc
  real(kind=sgl)          :: tfraction
  character(6)            :: Kosselmode
  character(fnlen)        :: xtalname
  character(fnlen)        :: outname
  character(fnlen)        :: BetheParametersFile
end type KosselmasterNameListType

! class definition
type, public :: Kosselmaster_T
private
  character(fnlen)                :: nmldeffile = 'EMKosselmaster.nml'
  type(KosselmasterNameListType)  :: nml

contains
private
  procedure, pass(self) :: readNameList_
  procedure, pass(self) :: writeHDFNameList_
  procedure, pass(self) :: getNameList_
  procedure, pass(self) :: Kosselmaster_
  procedure, pass(self) :: get_numthick_
  procedure, pass(self) :: get_npx_
  procedure, pass(self) :: get_nthreads_
  procedure, pass(self) :: get_voltage_
  procedure, pass(self) :: get_dmin_
  procedure, pass(self) :: get_startthick_
  procedure, pass(self) :: get_thickinc_
  procedure, pass(self) :: get_tfraction_
  procedure, pass(self) :: get_Kosselmode_
  procedure, pass(self) :: get_xtalname_
  procedure, pass(self) :: get_outname_
  procedure, pass(self) :: get_BetheParametersFile_
  procedure, pass(self) :: set_numthick_
  procedure, pass(self) :: set_npx_
  procedure, pass(self) :: set_nthreads_
  procedure, pass(self) :: set_voltage_
  procedure, pass(self) :: set_dmin_
  procedure, pass(self) :: set_startthick_
  procedure, pass(self) :: set_thickinc_
  procedure, pass(self) :: set_tfraction_
  procedure, pass(self) :: set_Kosselmode_
  procedure, pass(self) :: set_xtalname_
  procedure, pass(self) :: set_outname_
  procedure, pass(self) :: set_BetheParametersFile_

  generic, public :: getNameList => getNameList_
  generic, public :: writeHDFNameList => writeHDFNameList_
  generic, public :: readNameList => readNameList_
  generic, public :: Kosselmaster => Kosselmaster_
  generic, public :: get_numthick => get_numthick_
  generic, public :: get_npx => get_npx_
  generic, public :: get_nthreads => get_nthreads_
  generic, public :: get_voltage => get_voltage_
  generic, public :: get_dmin => get_dmin_
  generic, public :: get_startthick => get_startthick_
  generic, public :: get_thickinc => get_thickinc_
  generic, public :: get_tfraction => get_tfraction_
  generic, public :: get_Kosselmode => get_Kosselmode_
  generic, public :: get_xtalname => get_xtalname_
  generic, public :: get_outname => get_outname_
  generic, public :: get_BetheParametersFile => get_BetheParametersFile_
  generic, public :: set_numthick => set_numthick_
  generic, public :: set_npx => set_npx_
  generic, public :: set_nthreads => set_nthreads_
  generic, public :: set_voltage => set_voltage_
  generic, public :: set_dmin => set_dmin_
  generic, public :: set_startthick => set_startthick_
  generic, public :: set_thickinc => set_thickinc_
  generic, public :: set_tfraction => set_tfraction_
  generic, public :: set_Kosselmode => set_Kosselmode_
  generic, public :: set_xtalname => set_xtalname_
  generic, public :: set_outname => set_outname_
  generic, public :: set_BetheParametersFile => set_BetheParametersFile_
end type Kosselmaster_T

! the constructor routine for this class
interface Kosselmaster_T
  module procedure Kosselmaster_constructor
end interface Kosselmaster_T

contains

!--------------------------------------------------------------------------
type(Kosselmaster_T) function Kosselmaster_constructor( nmlfile ) result(Kosselmaster)
!DEC$ ATTRIBUTES DLLEXPORT :: Kosselmaster_constructor
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! constructor for the Kosselmaster_T Class; reads the name list

IMPLICIT NONE

character(fnlen), OPTIONAL   :: nmlfile

call Kosselmaster%readNameList(nmlfile)

end function Kosselmaster_constructor

!--------------------------------------------------------------------------
subroutine Kosselmaster_destructor(self)
!DEC$ ATTRIBUTES DLLEXPORT ::  Kosselmaster_destructor
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! destructor for the Kosselmaster_T Class

IMPLICIT NONE

type(Kosselmaster_T), INTENT(INOUT)  :: self

call reportDestructor('Kosselmaster_T')

end subroutine Kosselmaster_destructor

!--------------------------------------------------------------------------
subroutine readNameList_(self, nmlfile, initonly)
!DEC$ ATTRIBUTES DLLEXPORT :: readNameList_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! read the namelist from an nml file for the Kosselmaster_T Class

use mod_io
use mod_EMsoft

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT) :: self
character(fnlen),INTENT(IN)          :: nmlfile
 !! full path to namelist file
logical,OPTIONAL,INTENT(IN)          :: initonly
 !! fill in the default values only; do not read the file

type(EMsoft_T)                       :: EMsoft
type(IO_T)                           :: Message
logical                              :: skipread = .FALSE.

integer(kind=irg)       :: numthick
integer(kind=irg)       :: npx
integer(kind=irg)       :: nthreads
real(kind=sgl)          :: voltage
real(kind=sgl)          :: dmin
real(kind=sgl)          :: startthick
real(kind=sgl)          :: thickinc
real(kind=sgl)          :: tfraction
character(6)            :: Kosselmode
character(fnlen)        :: xtalname
character(fnlen)        :: outname
character(fnlen)        :: BetheParametersFile

namelist /Kosselmasterlist/ xtalname, voltage, dmin, nthreads, BetheParametersFile, &
                            startthick, thickinc, numthick, tfraction, outname, npx, Kosselmode

! set the input parameters to default values (except for xtalname, which must be present)
numthick = 10                   ! number of increments
npx = 256                       ! output arrays will have size npix x npix
nthreads = 4                    ! default number of threads for OpenMP
voltage = 20.0                  ! acceleration voltage [kV]
dmin = 0.05                     ! smallest d-spacing to include in dynamical matrix [nm]
startthick = 10.0               ! starting thickness [nm]
thickinc = 10.0                 ! thickness increment
xtalname = 'undefined'          ! initial value to check that the keyword is present in the nml file
outname = 'Kosselout.data'      ! output filename
BetheParametersFile = 'BetheParameters.nml' ! file name for Bethe parameters
Kosselmode = 'normal'           ! 'thicks' for thickness determination, 'normal' for normal plot
tfraction = 0.1                 ! thickness fraction for 'thicks' mode

if (present(initonly)) then
  if (initonly) skipread = .TRUE.
end if

if (.not.skipread) then
! read the namelist file
 open(UNIT=dataunit,FILE=trim(nmlfile),DELIM='apostrophe',STATUS='old')
 read(UNIT=dataunit,NML=Kosselmasterlist)
 close(UNIT=dataunit,STATUS='keep')

! check for required entries
 if (trim(xtalname).eq.'undefined') then
  call Message%printError('readNameList:',' structure file name is undefined in '//nmlfile)
 end if
end if

! if we get here, then all appears to be ok, and we need to fill in the knl fields
self%nml%numthick = numthick
self%nml%npx = npx
self%nml%nthreads = nthreads
self%nml%voltage = voltage
self%nml%dmin = dmin
self%nml%startthick = startthick
self%nml%thickinc = thickinc
self%nml%tfraction = tfraction
self%nml%Kosselmode = Kosselmode
self%nml%xtalname = xtalname
self%nml%outname = outname
self%nml%BetheParametersFile = BetheParametersFile

end subroutine readNameList_

!--------------------------------------------------------------------------
function getNameList_(self) result(nml)
!DEC$ ATTRIBUTES DLLEXPORT :: getNameList_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! pass the namelist for the Kosselmaster_T Class to the calling program

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)          :: self
type(KosselmasterNameListType)                :: nml

nml = self%nml

end function getNameList_

!--------------------------------------------------------------------------
recursive subroutine writeHDFNameList_(self, HDF, HDFnames)
!DEC$ ATTRIBUTES DLLEXPORT :: writeHDFNameList_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! write namelist to HDF file

use mod_HDFsupport
use mod_HDFnames
use stringconstants

use ISO_C_BINDING

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)    :: self
type(HDF_T), INTENT(INOUT)              :: HDF
type(HDFnames_T), INTENT(INOUT)         :: HDFnames

integer(kind=irg),parameter             :: n_int = 3, n_real = 5
integer(kind=irg)                       :: hdferr,  io_int(n_int)
real(kind=sgl)                          :: io_real(n_real)
character(20)                           :: intlist(n_int), reallist(n_real)
character(fnlen)                        :: dataset, sval(1),groupname
character(fnlen,kind=c_char)            :: line2(1)

associate( knl => self%nml )

! create the group for this namelist
groupname = trim(HDFnames%get_NMLlist())
hdferr = HDF%createGroup(groupname)

! write all the single integers
io_int = (/ knl%numthick, knl%npx, knl%nthreads /)
intlist(1) = 'numthick'
intlist(2) = 'npx'
intlist(3) = 'nthreads'
call HDF%writeNMLintegers(io_int, intlist, n_int)

! write all the single reals
io_real = (/ knl%voltage, knl%dmin, knl%startthick, knl%thickinc, knl%tfraction /)
reallist(1) = 'voltage'
reallist(2) = 'dmin'
reallist(3) = 'startthick'
reallist(4) = 'thickinc'
reallist(5) = 'tfraction'
call HDF%writeNMLreals(io_real, reallist, n_real)

! write all the strings
dataset = SC_Kosselmode
line2(1) = knl%Kosselmode
hdferr = HDF%writeDatasetStringArray(dataset, line2, 1)
if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create Kosselmode dataset',hdferr)

dataset = SC_xtalname
line2(1) = knl%xtalname
hdferr = HDF%writeDatasetStringArray(dataset, line2, 1)
if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create xtalname dataset',hdferr)

dataset = SC_outname
line2(1) = knl%outname
hdferr = HDF%writeDatasetStringArray(dataset, line2, 1)
if (hdferr.ne.0) call HDF%error_check('writeHDFNameList: unable to create outname dataset',hdferr)

! and pop this group off the stack
call HDF%pop()

end associate

end subroutine writeHDFNameList_

!--------------------------------------------------------------------------
function get_numthick_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_numthick_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get numthick from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
integer(kind=irg)                        :: out

out = self%nml%numthick

end function get_numthick_

!--------------------------------------------------------------------------
subroutine set_numthick_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_numthick_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set numthick in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
integer(kind=irg), INTENT(IN)            :: inp

self%nml%numthick = inp

end subroutine set_numthick_

!--------------------------------------------------------------------------
function get_npx_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_npx_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get npx from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
integer(kind=irg)                        :: out

out = self%nml%npx

end function get_npx_

!--------------------------------------------------------------------------
subroutine set_npx_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_npx_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set npx in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
integer(kind=irg), INTENT(IN)            :: inp

self%nml%npx = inp

end subroutine set_npx_

!--------------------------------------------------------------------------
function get_nthreads_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_nthreads_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get nthreads from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
integer(kind=irg)                        :: out

out = self%nml%nthreads

end function get_nthreads_

!--------------------------------------------------------------------------
subroutine set_nthreads_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_nthreads_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set nthreads in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
integer(kind=irg), INTENT(IN)            :: inp

self%nml%nthreads = inp

end subroutine set_nthreads_

!--------------------------------------------------------------------------
function get_voltage_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_voltage_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get voltage from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl)                           :: out

out = self%nml%voltage

end function get_voltage_

!--------------------------------------------------------------------------
subroutine set_voltage_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_voltage_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set voltage in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl), INTENT(IN)               :: inp

self%nml%voltage = inp

end subroutine set_voltage_

!--------------------------------------------------------------------------
function get_dmin_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_dmin_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get dmin from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl)                           :: out

out = self%nml%dmin

end function get_dmin_

!--------------------------------------------------------------------------
subroutine set_dmin_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_dmin_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set dmin in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl), INTENT(IN)               :: inp

self%nml%dmin = inp

end subroutine set_dmin_

!--------------------------------------------------------------------------
function get_startthick_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_startthick_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get startthick from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl)                           :: out

out = self%nml%startthick

end function get_startthick_

!--------------------------------------------------------------------------
subroutine set_startthick_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_startthick_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set startthick in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl), INTENT(IN)               :: inp

self%nml%startthick = inp

end subroutine set_startthick_

!--------------------------------------------------------------------------
function get_thickinc_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_thickinc_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get thickinc from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl)                           :: out

out = self%nml%thickinc

end function get_thickinc_

!--------------------------------------------------------------------------
subroutine set_thickinc_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_thickinc_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set thickinc in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl), INTENT(IN)               :: inp

self%nml%thickinc = inp

end subroutine set_thickinc_

!--------------------------------------------------------------------------
function get_tfraction_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_tfraction_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get tfraction from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl)                           :: out

out = self%nml%tfraction

end function get_tfraction_

!--------------------------------------------------------------------------
subroutine set_tfraction_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_tfraction_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set tfraction in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
real(kind=sgl), INTENT(IN)               :: inp

self%nml%tfraction = inp

end subroutine set_tfraction_

!--------------------------------------------------------------------------
function get_Kosselmode_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_Kosselmode_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get Kosselmode from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(6)                             :: out

out = self%nml%Kosselmode

end function get_Kosselmode_

!--------------------------------------------------------------------------
subroutine set_Kosselmode_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_Kosselmode_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set Kosselmode in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(6), INTENT(IN)                 :: inp

self%nml%Kosselmode = inp

end subroutine set_Kosselmode_

!--------------------------------------------------------------------------
function get_xtalname_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_xtalname_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get xtalname from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(fnlen)                         :: out

out = self%nml%xtalname

end function get_xtalname_

!--------------------------------------------------------------------------
subroutine set_xtalname_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_xtalname_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set xtalname in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(fnlen), INTENT(IN)             :: inp

self%nml%xtalname = inp

end subroutine set_xtalname_

!--------------------------------------------------------------------------
function get_outname_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_outname_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get outname from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(fnlen)                         :: out

out = self%nml%outname

end function get_outname_

!--------------------------------------------------------------------------
subroutine set_outname_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_outname_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set outname in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(fnlen), INTENT(IN)             :: inp

self%nml%outname = inp

end subroutine set_outname_

!--------------------------------------------------------------------------
function get_BetheParametersFile_(self) result(out)
!DEC$ ATTRIBUTES DLLEXPORT :: get_BetheParametersFile_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! get BetheParametersFile from the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(fnlen)                         :: out

out = self%nml%BetheParametersFile

end function get_BetheParametersFile_

!--------------------------------------------------------------------------
subroutine set_BetheParametersFile_(self,inp)
!DEC$ ATTRIBUTES DLLEXPORT :: set_BetheParametersFile_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! set BetheParametersFile in the Kosselmaster_T class

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)     :: self
character(fnlen), INTENT(IN)             :: inp

self%nml%BetheParametersFile = inp

end subroutine set_BetheParametersFile_

!--------------------------------------------------------------------------
subroutine Kosselmaster_(self, EMsoft, progname)
!DEC$ ATTRIBUTES DLLEXPORT :: Kosselmaster_
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! perform the computations

use mod_EMsoft
use mod_initializers
use mod_symmetry
use mod_crystallography
use mod_gvectors
use mod_kvectors
use mod_io
use mod_diffraction
use mod_timing
use mod_Lambert
use mod_MPfiles
use mod_math
use HDF5
use mod_HDFsupport
use mod_HDFnames
use ISO_C_BINDING
use omp_lib
use mod_OMPsupport
use stringconstants

IMPLICIT NONE

class(Kosselmaster_T), INTENT(INOUT)    :: self
type(EMsoft_T), INTENT(INOUT)           :: EMsoft
character(fnlen), INTENT(INOUT)         :: progname

type(Cell_T)                    :: cell
type(DynType)                   :: Dyn
type(Timing_T)                  :: timer
type(IO_T)                      :: Message
type(Lambert_T)                 :: L
type(HDF_T)                     :: HDF
type(SpaceGroup_T)              :: SG
type(Diffraction_T)             :: Diff
type(MPfile_T)                  :: MPFT
type(kvectors_T)                :: kvec
type(gvectors_T)                :: reflist
type(HDFnames_T)                :: HDFnames

real(kind=dbl)                  :: ctmp(192,3), arg
integer(HSIZE_T)                :: dims3(3), cnt3(3), offset3(3)
integer(HSIZE_T)                :: dims2(2), cnt2(2), offset2(2)
integer(kind=irg)               :: isym,i,j,ik,npy,ipx,ipy,ipz,debug,izz, izzmax, iequiv(3,48), nequiv, num_el, MCnthreads, & ! counters
                                   SamplingType,numk,numthick,  nthreads, & ! number of independent incident beam directions
                                   ir,kk(3), npyhex, skip, ijmax, one, NUMTHREADS, TID, hdferr, tickstart, &
                                   n,ix,iy, io_int(6), nns, nnw, nref, nix, niy, nixp, niyp, ierr, &
                                   istat,gzero,ic,ip,ikk, totstrong, totweak     ! counters
real(kind=dbl)                  :: tpi,Znsq, kkl, DBWF, kin, xy(2), dc(3), edge, scl, tmp, dx, dxm, dy, dym, xyz(3), Radius !!
real(kind=sgl)                  :: io_real(5), selE, kn, FN(3), kkk(3), bp(4), tstop, xyzs(3), tav
complex(kind=dbl)               :: czero
real(kind=sgl),allocatable      :: mLPNH(:,:,:), mLPSH(:,:,:), Iz(:), thick(:), trange(:,:), masterSPNH(:,:,:), &
                                   masterSPSH(:,:,:), trangeSP(:,:)
real(kind=sgl),allocatable      :: auxNH(:,:,:), auxSH(:,:,:), auxtrange(:,:)
logical                         :: usehex, switchmirror, verbose, insert=.TRUE., overwrite=.TRUE., silent=.TRUE., g_exists
character(11)                   :: dstr
character(15)                   :: tstrb
character(15)                   :: tstre
character(fnlen)                :: oldprogname, groupname, energyfile, outname, attributename, datagroupname, HDF_FileVersion
character(fnlen, KIND=c_char),allocatable,TARGET :: stringarray(:)
character(fnlen,kind=c_char)                     :: line2(1)

type(reflisttype), pointer      :: firstw
type(kvectorlist), pointer      :: ktmp
type(gnode),save                :: rlp
type(BetheParameterType)        :: BetheParameters
real(kind=sgl),allocatable      :: karray(:,:)
integer(kind=irg),allocatable   :: kij(:,:)
complex(kind=dbl),allocatable   :: DynMat(:,:)
character(fnlen)                :: dataset, instring

!$OMP THREADPRIVATE(rlp)

call openFortranHDFInterface()

associate( kmnl=>self%nml )

! set the HDF group names for this program
HDF = HDF_T()
HDFnames = HDFnames_T()

! initialize the timing routines
timer = Timing_T()
tstrb = timer%getTimeString()

tpi = 2.D0*cPi
czero = cmplx(0.D0,0.D0)
npy = kmnl%npx
gzero = 1  ! index of incident beam

!=============================================
!=============================================
! crystallography section
verbose = .TRUE.

call cell%setFileName(kmnl%xtalname)
call Diff%setrlpmethod('WK')

call Diff%setV(dble(kmnl%voltage))
call Initialize_Cell(cell, Diff, SG, Dyn, EMsoft, kmnl%dmin, verbose, useHDF=HDF)

! check the crystal system and setting; abort the program for trigonal with rhombohedral setting with
! an explanation for the user

if ((SG%getSpaceGroupXtalSystem().eq.5).and.(cell%getLatParm('b').eq.cell%getLatParm('c'))) then
    call Message%printMessage( (/ &
    '                                                                         ', &
    ' ========Program Aborted========                                         ', &
    ' The EBSD master pattern simulation for rhombohedral/trigonal structures ', &
    ' requires that the structure be described using the hexagonal reference  ', &
    ' frame.  Please re-enter the crystal structure in this setting and re-run', &
    ' the Monte Carlo calculation and this master pattern program.            '/) )
    stop
end if

! determine the point group number
 j=0
 do i=1,32
  if (SGPG(i).le.SG%getSpaceGroupNumber()) j=i
 end do
 isym = j

! Here, we encode isym into a new number that describes the sampling scheme; the new schemes are
! described in detail in the EBSD manual pdf file.

SamplingType = PGSamplingType(isym)

! next, intercept the special cases (hexagonal vs. rhombohedral cases that require special treatment)
if ((SamplingType.eq.-1).or.(isym.eq.14).or.(isym.eq.26)) then
  SamplingType = SG%getHexvsRho(isym)
end if

! if the point group is trigonal or hexagonal, we need to switch usehex to .TRUE. so that
! the program will use the hexagonal sampling method
usehex = .FALSE.
if ((SG%getSpaceGroupXtalSystem().eq.4).or.(SG%getSpaceGroupXtalSystem().eq.5)) usehex = .TRUE.

! ---------- end of symmetry and crystallography section
!=============================================
!=============================================

!=============================================
!=============================================
! this is where we determine the values for the thicknesses
if (kmnl%Kosselmode.eq.'thicks') then
  numthick = 1
else
  numthick = kmnl%numthick
end if

allocate(thick(numthick),Iz(numthick))

do i=1,numthick
  thick(i) = kmnl%startthick + float(i-1)*kmnl%thickinc
end do
!=============================================
!=============================================

!=============================================
!=============================================
! ---------- allocate memory for the master pattern
! we need to sample the stereographic projection Northern hemisphere or a portion
! thereoff, depending on the current symmetry.
if (kmnl%Kosselmode.eq.'normal') then
  allocate(mLPNH(-kmnl%npx:kmnl%npx,-npy:npy,1:kmnl%numthick),stat=istat)
  allocate(mLPSH(-kmnl%npx:kmnl%npx,-npy:npy,1:kmnl%numthick),stat=istat)
  allocate(masterSPNH(-kmnl%npx:kmnl%npx,-npy:npy,1:kmnl%numthick))
  allocate(masterSPSH(-kmnl%npx:kmnl%npx,-npy:npy,1:kmnl%numthick))

! set various arrays to zero
   mLPNH = 0.0
   mLPSH = 0.0
else
! Kosselmode must be 'thicks'
  allocate(trange(-kmnl%npx:kmnl%npx,-npy:npy),trangeSP(-kmnl%npx:kmnl%npx,-npy:npy),stat=istat)
  trange = 0.0
  trangeSP = 0.0
end if
! ---------- end allocate memory for the master pattern
!=============================================
!=============================================

! force dynamical matrix routine to read new Bethe parameters from file
  call Diff%SetBetheParameters(EMsoft, .FALSE., kmnl%BetheParametersFile)

!=============================================
! create the HDF5 output file
!=============================================
  HDF = HDF_T()
  call HDFnames%set_ProgramData(SC_Kosselmaster)
  call HDFnames%set_NMLlist(SC_KosselmasterNameList)
  call HDFnames%set_NMLfilename(SC_KosselmasterNML)

  outname = EMsoft%generateFilePath('EMdatapathname',trim(kmnl%outname))

! Create a new file using the default properties.
  hdferr =  HDF%createFile(outname)

! write the EMheader to the file
groupname = trim(HDFnames%get_ProgramData())
  call HDF%writeEMheader(EMsoft,dstr, tstrb, tstre, progname, groupname)

! add the CrystalData group at the top level of the file
call cell%addXtalDataGroup(SG, EMsoft, HDF)

! create a namelist group to write all the namelist files into
  hdferr = HDF%createGroup(HDFnames%get_NMLfiles())

! read the text file and write the array to the file
dataset = trim(HDFnames%get_NMLfilename())
  hdferr = HDF%writeDatasetTextFile(dataset, EMsoft%nmldeffile)

! leave this group
  call HDF%pop()

! create a namelist group to write all the namelist files into
  hdferr = HDF%createGroup(HDFnames%get_NMLparameters())
  call self%writeHDFNameList(HDF, HDFnames)
  call Diff%writeBetheparameterNameList(HDF)

! leave this group
  call HDF%pop()

! then the remainder of the data in a EMData group
  hdferr = HDF%createGroup(HDFnames%get_EMData())
  hdferr = HDF%createGroup(HDFnames%get_ProgramData())

! add a HDF_FileVersion attribute
  HDF_FileVersion = '4.0'
  HDF_FileVersion = cstringify(HDF_FileVersion)
  attributename = SC_HDFFileVersion
  hdferr = HDF%addStringAttributeToGroup(attributename, HDF_FileVersion)

dataset = SC_xtalname
  allocate(stringarray(1))
  stringarray(1)= trim(kmnl%xtalname)
  call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
  if (g_exists) then
    hdferr = HDF%writeDatasetStringArray(dataset, stringarray, 1, overwrite)
  else
    hdferr = HDF%writeDatasetStringArray(dataset, stringarray, 1)
  end if

dataset = SC_npx
  hdferr = HDF%writeDatasetInteger(dataset, kmnl%npx)

dataset = SC_BetheParameters
  bp(1) = Diff%getBetheParameter('c1')
  bp(2) = Diff%getBetheParameter('c2')
  bp(3) = Diff%getBetheParameter('c3')
  bp(4) = Diff%getBetheParameter('sg')
  call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
  if (g_exists) then
    hdferr = HDF%writeDatasetFloatArray(dataset, bp, 4, overwrite)
  else
    hdferr = HDF%writeDatasetFloatArray(dataset, bp, 4)
  end if

dataset = SC_numthick
  call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
  if (g_exists) then
    hdferr = HDF%writeDatasetInteger(dataset, kmnl%numthick, overwrite)
  else
    hdferr = HDF%writeDatasetInteger(dataset, kmnl%numthick)
  end if

dataset = SC_startthick
  call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
  if (g_exists) then
    hdferr = HDF%writeDatasetFloat(dataset, kmnl%startthick, overwrite)
  else
    hdferr = HDF%writeDatasetFloat(dataset, kmnl%startthick)
  end if

dataset = SC_thickinc
  call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
  if (g_exists) then
    hdferr = HDF%writeDatasetFloat(dataset, kmnl%thickinc, overwrite)
  else
    hdferr = HDF%writeDatasetFloat(dataset, kmnl%thickinc)
  end if

dataset = SC_Kosselmode
  stringarray(1)= trim(kmnl%Kosselmode)
  hdferr = HDF%writeDatasetStringArray(dataset, stringarray, 1)

! create the hyperslabs and write zeroes to them for now
  if (kmnl%Kosselmode.eq.'normal') then
dataset = SC_mLPNH
    dims3 = (/  2*kmnl%npx+1, 2*kmnl%npx+1, numthick /)
    cnt3 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1, numthick /)
    offset3 = (/ 0, 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, mLPNH, dims3, offset3, cnt3)

dataset = SC_mLPSH
    dims3 = (/  2*kmnl%npx+1, 2*kmnl%npx+1, numthick /)
    cnt3 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1, numthick /)
    offset3 = (/ 0, 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, mLPSH, dims3, offset3, cnt3)
  else
dataset = SC_trange
    dims2 = (/  2*kmnl%npx+1, 2*kmnl%npx+1 /)
    cnt2 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1 /)
    offset2 = (/ 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, trange, dims2, offset2, cnt2)

dataset = SC_trangeSP
    dims2 = (/  2*kmnl%npx+1, 2*kmnl%npx+1 /)
    cnt2 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1 /)
    offset2 = (/ 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, trangeSP, dims2, offset2, cnt2)
  end if

  call HDF%pop(.TRUE.)
!=============================================
! completes the HDF5 output file
!=============================================

!=============================================
!=============================================
! print a message to indicate we're starting the computation
call Message%printMessage('Starting computation', frm = "(/A)")

!=============================================
! ---------- create the incident beam directions list
! determine all independent incident beam directions (use a linked list starting at khead)
! numk is the total number of k-vectors to be included in this computation;
! note that this needs to be redone for each energy, since the wave vector changes with energy
kvec = kvectors_T()   ! initialize the wave vector list
call kvec%set_kinp( (/ 0.D0, 0.D0, 1.D0 /) )
call kvec%set_ktmax( 0.D0 )
call kvec%set_SamplingType( SamplingType )

call kvec%set_mapmode('RoscaLambert')
if (usehex) then
  call kvec%Calckvectors(cell, SG, Diff, (/ 0.D0, 0.D0, 0.D0 /),kmnl%npx,npy, ijmax,usehex)
else
  call kvec%Calckvectors(cell, SG, Diff, (/ 0.D0, 0.D0, 0.D0 /),kmnl%npx,npy, ijmax,usehex)
end if
numk = kvec%get_numk()
io_int(1)=numk
call Message%WriteValue('# independent beam directions to be considered = ', io_int, 1, "(I8)")

! convert part of the kvector linked list into arrays for OpenMP
allocate(karray(4,numk), kij(3,numk),stat=istat)
! point to the first beam direction
ktmp => kvec%get_ListHead()
! and loop through the list, keeping k, kn, and i,j
karray(1:3,1) = sngl(ktmp%k(1:3))
karray(4,1) = sngl(ktmp%kn)
kij(1:3,1) = (/ ktmp%i, ktmp%j, ktmp%hs /)
 do ik=2,numk
   ktmp => ktmp%next
   karray(1:3,ik) = sngl(ktmp%k(1:3))
   karray(4,ik) = sngl(ktmp%kn)
   kij(1:3,ik) = (/ ktmp%i, ktmp%j, ktmp%hs /)
 end do
! and remove the linked list
call kvec%Delete_kvectorlist()

verbose = .FALSE.
totstrong = 0
totweak = 0

! ---------- end of "create the incident beam directions list"
!=============================================

! here's where we introduce the OpenMP calls, to speed up the overall calculations...
call timer%Time_tick(1)

! set the number of OpenMP threads
call OMP_setNThreads(kmnl%nthreads)

! use OpenMP to run on multiple cores ...
!$OMP PARALLEL COPYIN(rlp) &
!$OMP& PRIVATE(DynMat,ik,FN,TID,kn,ipx,ipy,ix,iequiv,nequiv,reflist,firstw) &
!$OMP& PRIVATE(kkk,nns,nnw,nref,io_int,Iz)

NUMTHREADS = OMP_GET_NUM_THREADS()
TID = OMP_GET_THREAD_NUM()

!$OMP DO SCHEDULE(DYNAMIC,100)
! ---------- and here we start the beam direction loop
beamloop:do ik = 1,numk
!=============================================
! ---------- create the master reflection list for this beam direction
! Then we must determine the masterlist of reflections (also a linked list);
! This list basically samples a large reciprocal space volume; it does not
! distinguish between zero and higher order Laue zones, since that
! distinction becomes meaningless when we consider the complete
! reciprocal lattice.
  reflist = gvectors_T()
  kkk = karray(1:3,ik)
  FN = kkk

  call reflist%Initialize_ReflectionList(cell, SG, Diff, FN, kkk, self%nml%dmin, verbose)
  nref = reflist%get_nref()
! ---------- end of "create the master reflection list"
!=============================================

! determine strong and weak reflections
  nullify(firstw)
  nns = 0
  nnw = 0
  call reflist%Apply_BethePotentials(Diff, firstw, nns, nnw)

! generate the dynamical matrix
  if (allocated(DynMat)) deallocate(DynMat)
  allocate(DynMat(nns,nns))
  call reflist%GetDynMat(cell, Diff, firstw, DynMat, nns, nnw)
  totstrong = totstrong + nns
  totweak = totweak + nnw

! solve the dynamical eigenvalue equation for this beam direction
  kn = karray(4,ik)
  if (self%nml%Kosselmode.eq.'thicks') then
    call CalcKthick(DynMat,kn,nns,self%nml%tfraction,Iz)
  else
    call CalcKint(DynMat,kn,nns,numthick,thick,Iz)
  end if
  deallocate(DynMat)

! and store the resulting svals values, applying point group symmetry where needed.
  ipx = kij(1,ik)
  ipy = kij(2,ik)
  ipz = kij(3,ik)
!
  if (usehex) then
   call L%Apply3DPGSymmetry(cell,SG,ipx,ipy,ipz,self%nml%npx,iequiv,nequiv,usehex)
  else
   if ((SG%getSpaceGroupNumber().ge.195).and.(SG%getSpaceGroupNumber().le.230)) then
     call L%Apply3DPGSymmetry(cell,SG,ipx,ipy,ipz,self%nml%npx,iequiv,nequiv,cubictype=SamplingType)
   else
     call L%Apply3DPGSymmetry(cell,SG,ipx,ipy,ipz,self%nml%npx,iequiv,nequiv)
   end if
  end if
!$OMP CRITICAL
  if (self%nml%Kosselmode.eq.'normal') then
    do ix=1,nequiv
      if (iequiv(3,ix).eq.-1) mLPSH(iequiv(1,ix),iequiv(2,ix),1:numthick) = Iz(1:numthick)
      if (iequiv(3,ix).eq.1) mLPNH(iequiv(1,ix),iequiv(2,ix),1:numthick) = Iz(1:numthick)
    end do
  else
    do ix=1,nequiv
      trange(iequiv(1,ix),iequiv(2,ix)) = Iz(1)
    end do
  end if
!$OMP END CRITICAL

  if (mod(ik,5000).eq.0) then
   io_int(1) = ik
   io_int(2) = numk
   call Message%WriteValue('  completed beam direction ',io_int, 2, "(I8,' of ',I8)")
  end if

  call reflist%Delete_gvectorlist()

end do beamloop

! end of OpenMP portion
!$OMP END PARALLEL

deallocate(karray, kij)

if (usehex) then
! and finally, we convert the hexagonally sampled array to a square Lambert projection which will be used
! for all EBSD pattern interpolations;  we need to do this for both the Northern and Southern hemispheres

! we begin by allocating auxiliary arrays to hold copies of the hexagonal data; the original arrays will
! then be overwritten with the newly interpolated data.
  if (self%nml%Kosselmode.eq.'normal') then
    allocate(auxNH(-kmnl%npx:kmnl%npx,-npy:npy,1:kmnl%numthick),stat=istat)
    allocate(auxSH(-kmnl%npx:kmnl%npx,-npy:npy,1:kmnl%numthick),stat=istat)
    auxNH = mLPNH
    auxSH = mLPSH
  else
    allocate(auxtrange(-kmnl%npx:kmnl%npx,-npy:npy),stat=istat)
    auxtrange = trange
  end if

  edge = 1.D0 / dble(kmnl%npx)
  scl = float(kmnl%npx)
  do i=-kmnl%npx,kmnl%npx
    do j=-npy,npy
! determine the spherical direction for this point
      L = Lambert_T( xyd = (/ dble(i), dble(j) /) * edge )
      ierr = L%LambertSquareToSphere(dc)
! convert direction cosines to hexagonal Lambert projections
      L = Lambert_T( xyzd = dc )
      ierr = L%LambertSphereToHex(xy)
      xy = xy * scl
! interpolate intensity from the neighboring points
      if (ierr.eq.0) then
        nix = floor(xy(1))
        niy = floor(xy(2))
        nixp = nix+1
        niyp = niy+1
        if (nixp.gt.kmnl%npx) nixp = nix
        if (niyp.gt.kmnl%npx) niyp = niy
        dx = xy(1) - nix
        dy = xy(2) - niy
        dxm = 1.D0 - dx
        dym = 1.D0 - dy
        if (kmnl%Kosselmode.eq.'normal') then
          mLPNH(i,j,1:numthick) = auxNH(nix,niy,1:numthick)*dxm*dym + auxNH(nixp,niy,1:numthick)*dx*dym + &
                                  auxNH(nix,niyp,1:numthick)*dxm*dy + auxNH(nixp,niyp,1:numthick)*dx*dy
          mLPSH(i,j,1:numthick) = auxSH(nix,niy,1:numthick)*dxm*dym + auxSH(nixp,niy,1:numthick)*dx*dym + &
                                  auxSH(nix,niyp,1:numthick)*dxm*dy + auxSH(nixp,niyp,1:numthick)*dx*dy
        else
          trange(i,j) = auxtrange(nix,niy)*dxm*dym + auxtrange(nixp,niy)*dx*dym + &
                        auxtrange(nix,niyp)*dxm*dy + auxtrange(nixp,niyp)*dx*dy
        end if
      end if
    end do
  end do
  deallocate(auxNH, auxSH)
end if

! make sure that the outer pixel rim of the mLPSH patterns is identical to
! that of the mLPNH array.
if (kmnl%Kosselmode.eq.'normal') then
  mLPSH(-kmnl%npx,-kmnl%npx:kmnl%npx,1:numthick) = mLPNH(-kmnl%npx,-kmnl%npx:kmnl%npx,1:numthick)
  mLPSH( kmnl%npx,-kmnl%npx:kmnl%npx,1:numthick) = mLPNH( kmnl%npx,-kmnl%npx:kmnl%npx,1:numthick)
  mLPSH(-kmnl%npx:kmnl%npx,-kmnl%npx,1:numthick) = mLPNH(-kmnl%npx:kmnl%npx,-kmnl%npx,1:numthick)
  mLPSH(-kmnl%npx:kmnl%npx, kmnl%npx,1:numthick) = mLPNH(-kmnl%npx:kmnl%npx, kmnl%npx,1:numthick)

! get stereographic projections (summed over the atomic positions)
  Radius = 1.0
  do i=-kmnl%npx,kmnl%npx
    do j=-kmnl%npx,kmnl%npx
      L = Lambert_T( xyd = (/ dble(i), dble(j) /) / dble(kmnl%npx) )
      ierr = L%StereoGraphicInverse( xyz, Radius )
      xyz = xyz/vecnorm(xyz)
      if (ierr.ne.0) then
        masterSPNH(i,j,1) = 0.0
        masterSPSH(i,j,1) = 0.0
      else
        masterSPNH(i,j,1:numthick) = InterpolateLambert(xyz, mLPNH, kmnl%npx, numthick)
        masterSPSH(i,j,1:numthick) = InterpolateLambert(xyz, mLPSH, kmnl%npx, numthick)
      end if
    end do
  end do
else
  ! set the background to the average intensity level
  tav = sum(trange)/float(2*kmnl%npx+1)**2
  ! get stereographic projections (summed over the atomic positions)
  Radius = 1.0
  do i=-kmnl%npx,kmnl%npx
    do j=-kmnl%npx,kmnl%npx
      L = Lambert_T( xyd = (/ dble(i), dble(j) /) / dble(kmnl%npx) )
      ierr = L%StereoGraphicInverse( xyz, Radius )
      xyz = xyz/vecnorm(xyz)
      if (ierr.ne.0) then
        trangeSP(i,j) = tav
      else
        xyzs = sngl(xyz)
        trangeSP(i,j) = InterpolateLambert(xyzs, trange, kmnl%npx)
      end if
    end do
  end do
end if

io_int(1) = nint(float(totstrong)/float(numk))
call Message%WriteValue(' -> Average number of strong reflections = ',io_int, 1, "(I5)")
io_int(1) = nint(float(totweak)/float(numk))
call Message%WriteValue(' -> Average number of weak reflections   = ',io_int, 1, "(I5)")

call timer%makeTimeStamp()
dstr = timer%getDateString()
tstre = timer%getTimeString()

!===================================
!===================================
! and write the final output arrays to the HDF5 file
datagroupname = trim(HDFnames%get_ProgramData())

! open the existing file using the default properties.
hdferr =  HDF%openFile(outname)

! update the time string
hdferr = HDF%openGroup(HDFnames%get_EMheader())
hdferr = HDF%openGroup(datagroupname)

dataset = SC_StopTime
  call timer%Time_tock(1)
  tstop = timer%getInterval(1)
  call timer%Time_reset(1)
  line2(1) = dstr//', '//tstre
  hdferr = HDF%writeDatasetStringArray(dataset, line2, 1, overwrite)

  io_int(1) = tstop
  call Message%WriteValue(' Execution time [s]: ',io_int,1)

dataset = SC_Duration
  call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
  if (g_exists) then
    hdferr = HDF%writeDatasetFloat(dataset, tstop, overwrite)
  else
    hdferr = HDF%writeDatasetFloat(dataset, tstop)
  end if

  call HDF%pop()
  call HDF%pop()

  hdferr = HDF%openGroup(HDFnames%get_EMData())
  hdferr = HDF%openGroup(datagroupname)

! add data to the hyperslab
if (kmnl%Kosselmode.eq.'normal') then
  dataset = SC_mLPNH
    dims3 = (/  2*kmnl%npx+1, 2*kmnl%npx+1, numthick /)
    cnt3 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1, numthick /)
    offset3 = (/ 0, 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, mLPNH, dims3, offset3, cnt3, insert)

  dataset = SC_mLPSH
    hdferr = HDF%writeHyperslabFloatArray(dataset, mLPSH, dims3, offset3, cnt3, insert)

  dataset = SC_masterSPNH
    call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
    if (g_exists) then
      hdferr = HDF%writeHyperslabFloatArray(dataset, masterSPNH, dims3, offset3, cnt3, insert)
    else
      hdferr = HDF%writeHyperslabFloatArray(dataset, masterSPNH, dims3, offset3, cnt3)
    end if

  dataset = SC_masterSPSH
    call H5Lexists_f(HDF%getobjectID(),trim(dataset),g_exists, hdferr)
    if (g_exists) then
      hdferr = HDF%writeHyperslabFloatArray(dataset, masterSPSH, dims3, offset3, cnt3, insert)
    else
      hdferr = HDF%writeHyperslabFloatArray(dataset, masterSPSH, dims3, offset3, cnt3)
    end if

else
  dataset = SC_trange
    dims2 = (/  2*kmnl%npx+1, 2*kmnl%npx+1 /)
    cnt2 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1 /)
    offset2 = (/ 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, trange, dims2, offset2, cnt2, insert)

  dataset = SC_trangeSP
    dims2 = (/  2*kmnl%npx+1, 2*kmnl%npx+1 /)
    cnt2 = (/ 2*kmnl%npx+1, 2*kmnl%npx+1 /)
    offset2 = (/ 0, 0 /)
    hdferr = HDF%writeHyperslabFloatArray(dataset, trangeSP, dims2, offset2, cnt2, insert)
end if

call HDF%pop(.TRUE.)

! and close the fortran hdf interface
call closeFortranHDFInterface()

call Message%printMessage('Final data stored in file '//trim(kmnl%outname), frm = "(A/)")

end associate

end subroutine Kosselmaster_

!--------------------------------------------------------------------------
recursive subroutine CalcKthick(DynMat,kn,nn,thresh,Iz)
!DEC$ ATTRIBUTES DLLEXPORT :: CalcKthick
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! compute the thickness for which the Kossel intensity drops below a treshold

use mod_io
use mod_diffraction
use mod_kvectors
use mod_gvectors

IMPLICIT NONE

integer(kind=irg),INTENT(IN)    :: nn                   !< number of strong beams
complex(kind=dbl),INTENT(IN)    :: DynMat(nn,nn)
real(kind=sgl),INTENT(IN)       :: kn
real(kind=sgl),INTENT(IN)       :: thresh               !< thickness fraction parameter
real(kind=sgl),INTENT(INOUT)    :: Iz(1)                !< output (thickness)
!f2py intent(in,out) ::  Iz

type(diffraction_T)             :: Diff
integer(kind=irg)               :: j, IPIV(nn), k
complex(kind=dbl)               :: CGinv(nn,nn), Minp(nn,nn), Wloc(nn), lCG(nn,nn), lW(nn), lalpha(nn)
real(kind=dbl)                  :: s(nn), q(nn), t, ss

! compute the eigenvalues and eigenvectors
 Minp = DynMat
 IPIV = 0
 call Diff%BWsolve(Minp,Wloc,lCG,CGinv,nn,IPIV)

! the alpha coefficients are in the first column of the inverse matrix
 lW = cPi*Wloc/cmplx(kn,0.0)
 lalpha(1:nn) = CGinv(1:nn,1)

! make sure the alpha excitation coefficients are normalized
 ss = sum(abs(lalpha(1:nn))**2)
 if (ss.ne.1.D0) then
  ss = cmplx(1.D0/dsqrt(ss),0.D0)
  lalpha = lalpha*ss
 endif

! compute the thickness value in steps of 0.25 nm until less than thresh
 do j=1,nn
  q(j) = -4.D0*cPi*aimag(lW(j))
  s(j) = abs(lalpha(j))**2
 end do
 t = 0.D0
 do
   t = t+0.25D0
   ss = sum(s*dexp(t*q))
   if (ss.le.thresh) EXIT
 end do
 Iz(1) = t

end subroutine CalcKthick

!--------------------------------------------------------------------------
recursive subroutine CalcKint(DynMat,kn,nn,nt,thick,Iz)
!DEC$ ATTRIBUTES DLLEXPORT :: CalcKint
!! author: MDG
!! version: 1.0
!! date: 03/25/20
!!
!! compute the Kossel intensities for a range of thicknesses

use mod_io
use mod_diffraction
use mod_kvectors
use mod_gvectors

IMPLICIT NONE

integer(kind=irg),INTENT(IN)    :: nn                   !< number of strong beams
complex(kind=dbl),INTENT(IN)    :: DynMat(nn,nn)
real(kind=sgl),INTENT(IN)       :: kn
integer(kind=irg),INTENT(IN)    :: nt                   !< number of thickness values
real(kind=sgl),INTENT(IN)       :: thick(nt)            !< thickness array
real(kind=sgl),INTENT(INOUT)    :: Iz(nt)               !< output intensities
!f2py intent(in,out) ::  Iz

type(diffraction_T)             :: Diff
integer(kind=irg)               :: j, IPIV(nn), k
complex(kind=dbl)               :: CGinv(nn,nn), Minp(nn,nn), Wloc(nn), lCG(nn,nn), lW(nn), lalpha(nn)
real(kind=dbl)                  :: s, q, t

! compute the eigenvalues and eigenvectors
 Minp = DynMat
 IPIV = 0
 call Diff%BWsolve(Minp,Wloc,lCG,CGinv,nn,IPIV)

! the alpha coefficients are in the first column of the inverse matrix
 lW = cPi*Wloc/cmplx(kn,0.0)
 lalpha(1:nn) = CGinv(1:nn,1)

! make sure the alpha excitation coefficients are normalized
 s = sum(abs(lalpha(1:nn))**2)
 if (s.ne.1.D0) then
  s = cmplx(1.D0/dsqrt(s),0.D0)
  lalpha = lalpha*s
 endif

! compute the thickness array
 Iz = 0.D0
 do j=1,nn
    q = -4.D0*cPi*aimag(lW(j))
    s = abs(lalpha(j))**2
    do k=1,nt
      t = q*thick(k)
      if (abs(t).lt.30.D0) Iz(k) = Iz(k) +  s * exp(t)
    end do
 end do

end subroutine CalcKint

end module mod_Kosselmaster

max_stars_repo_name: EMsoft-org/EMsoftOO

max_issues_repo_name: EMsoft-org/EMsoftOO

max_forks_repo_name: EMsoft-org/EMsoftOO


text:
Katie was a teenage mother at four months old, and was taken in to a shelter.  A UC Davis veternary student found her in that shelter and, since Katie was sick and going to be euthanised, that student rescued her from the shelter and cared for her at home.  One day, while studying for a music history test, Margie went over to this students house and met Katie for the first time, and it was love at first sight.  Several weeks later, when Katie was healthy enough to leave, Margie adopted Katie and theyve been best friends ever since.

Katie has lived with Users/AllisonMears and Users/NatalieKitchiner, among others who are not on DavisWiki.  People have described Katie to be a bit of a pest, but very cute and sweet.

When Katie was a kitten, her favorite thing to do was eat homework and bank statements.  This did not bode well with her owner or roommates.

Currently, Katie is living in Sunnyvale with Margie and her family, including three other cats by the names of Dotty, Daisy, and Stanley. Katie will move out when Margie does and live with her.   Katie can not wait to move out, because she hates Stanley, Dotty and Daisy, as they are mean and sort of smelly.

Katies full name is Her Royal Highness and Majesty Our Reverend Lady Katherine Galinda Sunshine Pussywilwo Douglass Shelby Dingleberry Halloran III Esq. PhD, which was a group effort by Margie and her friends.  Dont ask Margie what a pussywilwo is; she made it up one day whilst playing Oregon Trail and it stuck.

  Question: What did Her Royal Highness and Majesty Our Reverend Lady Katherine Galinda Sunshine Pussywilwo Douglass Shelby Dingleberry Halloran III Esq. PhD get her doctorate in?   Users/JeffTolentino

Katie received her PhD in Advanced Rug Peeing in 2007.

Katie likes to play with dental floss, chase laser pointers, look out the window, and sleep in the crook of peoples knees.  Her new favorite thing since moving to Sunnyvale is to run and frolick outside, and Margie hopes she can get a place with a yard so that Katie does not have to remain indoors.


text:
      SUBROUTINE simp3(a,mp,np,i1,k1,ip,kp)
      INTEGER i1,ip,k1,kp,mp,np
      REAL a(mp,np)
      INTEGER ii,kk
      REAL piv
      piv=1./a(ip+1,kp+1)
      do 12 ii=1,i1+1
        if(ii-1.ne.ip)then
          a(ii,kp+1)=a(ii,kp+1)*piv
          do 11 kk=1,k1+1
            if(kk-1.ne.kp)then
              a(ii,kk)=a(ii,kk)-a(ip+1,kk)*a(ii,kp+1)
            endif
11        continue
        endif
12    continue
      do 13 kk=1,k1+1
        if(kk-1.ne.kp)a(ip+1,kk)=-a(ip+1,kk)*piv
13    continue
      a(ip+1,kp+1)=piv
      return
      END

ext: for

max_stars_repo_path: NR-Functions/Numerical Recipes- Example & Functions/Functions/simp3.for


text:
module constants
!$$$   module documentation block
!                .      .    .                                       .
! module:    constants
!   prgmmr: treadon          org: np23                date: 2003-09-25
!
! abstract:  This module contains the definition of various constants
!            used in the gsi code
!
! program history log:
!   2003-09-25  treadon - original code
!   2004-03-02  treadon - allow global and regional constants to differ
!   2004-06-16  treadon - update documentation
!   2004-10-28  treadon - replace parameter tiny=1.e-12 with tiny_r_kind
!                         and tiny_single
!   2004-11-16  treadon - add huge_r_kind parameter
!
! Subroutines Included:
!   sub init_constants  - compute derived constants, set regional/global constants
!
! Variable Definitions:
!   see below
!
! attributes:
!   language: f90
!   machine:  ibm RS/6000 SP
!
!$$$
  use kinds, only: r_single,r_kind
  implicit none

! Declare constants
  integer izero
  real(r_kind) rearth,grav,omega,rd,rv,cp,cv,cvap,cliq
  real(r_kind) csol,hvap,hfus,psat,t0c,ttp,jcal,cp_mass
  real(r_kind) fv,deg2rad,rad2deg,pi,tiny_r_kind,huge_r_kind
  real(r_kind) ozcon,rozcon,tpwcon,rd_over_g,rd_over_cp,g_over_rd
  real(r_kind) amsua_clw_d1,amsua_clw_d2,constoz,zero,one,two,four
  real(r_kind) one_tenth,quarter,three,five
  real(r_kind) rearth_equator,stndrd_atmos_ps
  real(r_kind) semi_major_axis,semi_minor_axis,n_a,n_b
  real(r_kind) eccentricity,grav_polar,grav_ratio
  real(r_kind) grav_equator,earth_omega,grav_constant
  real(r_kind) flattening,eccentricity_linear,somigliana
  real(r_kind) dldt,dldti,hsub,psatk,tmix,xa,xai,xb,xbi
  real(r_kind) eps,epsm1,omeps
  real(r_kind) elocp,cpr,el2orc,cclimit,climit,epsq
  real(r_kind) pcpeff0,pcpeff1,pcpeff2,pcpeff3,rcp,c0,delta
  real(r_kind) h1000,factor1,factor2,rhcbot,rhctop,dx_max,dx_min,dx_inv
  real(r_kind) h300,half,cmr,cws,ke2,row,rrow
  real(r_single) zero_single,tiny_single,huge_single

! Define constants common to global and regional applications
!           name     value                  description                     units
!           ----     -----                  -----------                     -----
  parameter(rearth_equator= 6.37813662e6_r_kind) ! equatorial earth radius (m)
  parameter(omega  = 7.2921e-5_r_kind)  !  angular velocity of earth       (1/s)
  parameter(cp     = 1.0046e+3_r_kind)  !  specific heat of air @pressure  (J/kg/K)
  parameter(cvap   = 1.8460e+3_r_kind)  !  specific heat of h2o vapor      (J/kg/K)
  parameter(csol   = 2.1060e+3_r_kind)  !  specific heat of solid h2o (ice)(J/kg/K)
  parameter(hvap   = 2.5000e+6_r_kind)  !  latent heat of h2o condensation (J/kg)
  parameter(hfus   = 3.3358e+5_r_kind)  !  latent heat of h2o fusion       (J/kg)
  parameter(psat   = 6.1078e+2_r_kind)  !  pressure at h2o triple point    (Pa)
  parameter(t0c    = 2.7315e+2_r_kind)  !  temperature at zero celsius     (K)
  parameter(ttp    = 2.7316e+2_r_kind)  !  temperature at h2o triple point (K)
  parameter(jcal   = 4.1855e+0_r_kind)  !  joules per calorie              ()
  parameter(stndrd_atmos_ps = 1013.25e2_r_kind) ! 1976 US standard atmosphere ps   (Pa)

! Numeric constants
  parameter(izero  = 0)
  parameter(zero_single = 0.0_r_single)
  parameter(zero   = 0.0_r_kind)
  parameter(one_tenth  = 0.10_r_kind)
  parameter(quarter= 0.25_r_kind)
  parameter(one    = 1.0_r_kind)
  parameter(two    = 2.0_r_kind)
  parameter(three  = 3.0_r_kind)
  parameter(four   = 4.0_r_kind)
  parameter(five   = 5.0_r_kind)

! Constants for gps refractivity
  parameter(n_a=77.6_r_kind) !K/mb
  parameter(n_b=3.73e+5_r_kind) !K^2/mb

! Parameters below from WGS-84 model software inside GPS receivers.
  parameter(semi_major_axis = 6378.1370e3_r_kind)    !                     (m)
  parameter(semi_minor_axis = 6356.7523142e3_r_kind) !                     (m)
  parameter(grav_polar = 9.8321849378_r_kind)        !                     (m/s2)
  parameter(grav_equator = 9.7803253359_r_kind)      !                     (m/s2)
  parameter(earth_omega = 7.292115e-5_r_kind)        !                     (rad/s)
  parameter(grav_constant = 3.986004418e14_r_kind)   !                     (m3/s2)

! Derived geophysical constants
  parameter(flattening = (semi_major_axis-semi_minor_axis)/semi_major_axis)!()
  parameter(somigliana = &
       (semi_minor_axis/semi_major_axis) * (grav_polar/grav_equator) - one)!()
  parameter(grav_ratio = (earth_omega*earth_omega * &
       semi_major_axis*semi_major_axis * semi_minor_axis) / grav_constant) !()

! Derived thermodynamic constants
  parameter ( dldti = cvap-csol )
  parameter ( hsub = hvap+hfus )
  parameter ( psatk = psat*0.001_r_kind )
  parameter ( tmix = ttp-20._r_kind )
  parameter ( elocp = hvap/cp )
  parameter ( rcp  = one/cp )

! Constants used in GFS moist physics
  parameter ( h300 = 300._r_kind )
  parameter ( half = 0.5_r_kind )
  parameter ( cclimit = 0.001_r_kind )
  parameter ( climit = 1.e-20_r_kind)
  parameter ( epsq = 2.e-12_r_kind )
  parameter ( h1000 = 1000.0_r_kind)
  parameter ( rhcbot=0.85_r_kind )
  parameter ( rhctop=0.85_r_kind )
  parameter ( dx_max=-8.8818363_r_kind )
  parameter ( dx_min=-5.2574954_r_kind )
  parameter ( dx_inv=one/(dx_max-dx_min) )
  parameter ( c0=0.002_r_kind )
  parameter ( delta=0.6077338_r_kind )
  parameter ( pcpeff0=1.591_r_kind )
  parameter ( pcpeff1=-0.639_r_kind )
  parameter ( pcpeff2=0.0953_r_kind )
  parameter ( pcpeff3=-0.00496_r_kind )
  parameter ( cmr = one/0.0003_r_kind )
  parameter ( cws = 0.025_r_kind )
  parameter ( ke2 = 0.00002_r_kind )
  parameter ( row = 1000._r_kind )
  parameter ( rrow = one/row )

! Constant used to process ozone
  parameter ( constoz = 604229.0_r_kind)

! Constants used in cloud liquid water correction for AMSU-A
! brightness temperatures
  parameter ( amsua_clw_d1 = 0.754_r_kind )
  parameter ( amsua_clw_d2 = -2.265_r_kind )

contains
  subroutine init_constants(regional)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    init_constants     set derived, regional, global constants
!     prgmmr:    treadon          org: np23           date: 2004-03-02
!
! abstract:  This routine sets derived constants as well as those that
!            are specific to regional or global applications of the
!            gsi
!
! program history log:
!   2004-03-02  treadon
!   2004-06-16  treadon, documentation
!   2004-10-28  treadon - use intrinsic TINY function to set value
!                         for smallest machine representable positive
!                         number
!
!   input argument list:
!     regional - if .true., set regional gsi constants;
!                otherwise (.false.), use global constants
!
!   output argument list:
!
! attributes:
!   language: f90
!   machine:  ibm rs/6000 sp
!
!$$$
    logical regional
    real(r_kind) reradius,g,r_d,r_v,cliq_wrf

!   Trigonometric constants
    pi      = acos(-one)
    deg2rad = pi/180.0_r_kind
    rad2deg = one/deg2rad
    tiny_r_kind = tiny(zero)
    huge_r_kind = huge(zero)
    tiny_single = tiny(zero_single)
    huge_single = huge(zero_single)

!   Geophysical parameters used in conversion of geopotential to
!   geometric height
    eccentricity_linear = sqrt(semi_major_axis**2 - semi_minor_axis**2)
    eccentricity = eccentricity_linear / semi_major_axis

!   Define regional constants here
    if (regional) then

!      Name given to WRF constants
       reradius = one/6370.e03_r_kind
       g        = 9.81_r_kind
       r_d      = 287.04_r_kind
       r_v      = 461.6_r_kind
       cliq_wrf = 4190.0_r_kind
       cp_mass  = 1004.67_r_kind

!      Transfer WRF constants into unified GSI constants
       rearth = one/reradius
       grav   = g
       rd     = r_d
       rv     = r_v
       cv     = cp-r_d
       cliq   = cliq_wrf

!   Define global constants here
    else
       rearth = 6.3712e+6_r_kind
       grav   = 9.80665e+0_r_kind
       rd     = 2.8705e+2_r_kind
       rv     = 4.6150e+2_r_kind
       cv     = 7.1760e+2_r_kind
       cliq   = 4.1855e+3_r_kind
       cp_mass= zero
    endif

!   Now define derived constants which depend on constants
!   which differ between global and regional applications.

!   Constants related to ozone assimilation
    ozcon = grav*21.4e-9_r_kind
    rozcon= one/ozcon

!   Constant used in vertical integral for precipitable water
    tpwcon = 100.0_r_kind/grav

!   Derived atmospheric constants
    fv         = rv/rd-one    ! used in virtual temperature equation
    dldt       = cvap-cliq
    xa         = -(dldt/rv)
    xai        = -(dldti/rv)
    xb         = xa+hvap/(rv*ttp)
    xbi        = xai+hsub/(rv*ttp)
    eps        = rd/rv
    epsm1      = rd/rv-one
    omeps      = one-eps
    factor1    = (cvap-cliq)/rv
    factor2    = hvap/rv-factor1*t0c
    cpr        = cp*rd
    el2orc     = hvap*hvap/(rv*cp)
    rd_over_g  = rd/grav
    rd_over_cp = rd/cp
    g_over_rd  = grav/rd

    return
  end subroutine init_constants

end module constants

max_stars_repo_name: hsbadr/crtm

max_issues_repo_name: hsbadr/crtm

max_forks_repo_name: hsbadr/crtm