text:
!
!  Copyright 2017-2019 SALMON developers
!
!  Licensed under the Apache License, Version 2.0 (the "License");
!  you may not use this file except in compliance with the License.
!  You may obtain a copy of the License at
!
!      http://www.apache.org/licenses/LICENSE-2.0
!
!  Unless required by applicable law or agreed to in writing, software
!  distributed under the License is distributed on an "AS IS" BASIS,
!  WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
!  See the License for the specific language governing permissions and
!  limitations under the License.
!
!This file is "md_gs.f90"
!--------10--------20--------30--------40--------50--------60--------70--------80--------90--------100-------110-------120-------130

!-----------------------------------------------------------------
! Optimization in the ground state by conjugate gradient method
!-----------------------------------------------------------------
module md_ground_state
  use Global_Variables
  use salmon_parallel, only: nproc_group_global, nproc_id_global
  use salmon_communication, only: comm_bcast, comm_sync_all, comm_is_root
  use salmon_parallel, only: nproc_id_global
  implicit none
contains

subroutine calc_md_ground_state
  use timer
  use opt_variables
  use performance_analyzer
  use misc_routines, only: get_wtime
  use salmon_global, only: format3d, out_dns, out_dns_rt, out_dns_rt_step
  use ground_state
  use io_gs_wfn_k
  implicit none
  integer :: entrance_it,it,ia
  real(8) :: Vion_gs, Ework, Ework_integ_fdR_gs, dt_h, aforce(3,NI), Eall00
  real(8) :: Temperature_ion, mass_au, Rion_prev(3,NI)
  real(8) :: Htot, Enh, Enh_gkTlns, gkT, Qnh
  character(100) :: comment_line

  iflag_gs_init_wf = 1  !flag to skip giving randam number initial guess
  PrLv_scf = 0

  position_option='rewind'
  entrance_it=-1
  Rion_update_rt = rion_update_on

  ! Export electronic density
  if(out_dns == 'y') call write_density(it,'gs')

  ! Export to file_trj (initial step)
  if (out_rvf_rt=='y')then
       call Ion_Force_omp(Rion_update_rt,calc_mode_gs)
       write(comment_line,110) -1, 0.0d0
       if(ensemble=="NVT" .and. thermostat=="nose-hoover") &
       &  write(comment_line,112) trim(comment_line), xi_nh
       call write_xyz(comment_line,"new","rvf")
       call write_xyz(comment_line,"add","rvf")
  endif

  it=entrance_it+1
  dRion(:,:,it-1) = dRion(:,:,it) - velocity(:,:)*dt  !only for euler
  dt_h       = dt*0.5d0
  Eall00     = Eall0
  Enh_gkTlns = 0d0
  Enh        = 0d0

  if(ensemble=="NVT" .and. thermostat=="nose-hoover") then
     gkT = 3d0*NI * kB/hartree2J*temperature0_ion
     Qnh = gkT * thermostat_tau**2d0
  endif

  if(comm_is_root(nproc_id_global)) then
     write(*,'(a)') "# time, Tion, Vion_gs, Ework, Enh, Eall, Eall-Eall0, Htot, Temperature_ion [all in a.u.]"
  endif

  call comm_sync_all

  ! === Main loop of time step ===
  do it = entrance_it+1, Nt

     ! Velocity Verlet integrator

     !NHC act on velocity with dt/2
     if(ensemble=="NVT" .and. thermostat=="nose-hoover")then
        call apply_nose_hoover_velocity(dt_h)
     endif

     do ia=1,NI
        mass_au = umass*Mass(Kion(ia))
        velocity(:,ia) = velocity(:,ia) + force(:,ia)/mass_au * dt_h
     enddo

     if(step_velocity_scaling>=1 .and. mod(it,step_velocity_scaling)==0) then
        call cal_Tion_Temperature_ion(Tion,Temperature_ion,velocity)
        call apply_velocity_scaling_ion(Temperature_ion,velocity)
     endif

     Rion_prev(:,:) = Rion(:,:)
     do ia=1,NI
        Rion(:,ia) = Rion(:,ia) + velocity(:,ia)*dt
     enddo

     !put SHAKE here in future

     !NHC act on thermostat with dt
     if(ensemble=="NVT" .and. thermostat=="nose-hoover")then
        call cal_Tion_Temperature_ion(Tion,Temperature_ion,velocity)
        call apply_nose_hoover_thermostat(Temperature_ion,dt)
        Enh_gkTlns = Enh_gkTlns + gkT * xi_nh*dt
        Enh        = Enh_gkTlns + 0.5d0 * Qnh * xi_nh*xi_nh
     endif

     !update force (electric state) with updated coordinate
     aforce(:,:) = force(:,:)
     call prep_ps_periodic('update_all       ')
     call read_write_gs_wfn_k(iflag_read)
     call calc_ground_state
     call read_write_gs_wfn_k(iflag_write)

     call Ion_Force_omp(Rion_update_rt,calc_mode_gs)
     call Total_Energy_omp(Rion_update_rt,calc_mode_gs)

     aforce(:,:) = 0.5d0*( aforce(:,:) + force(:,:) )

     Ework = 0d0
     do ia=1,NI
        mass_au = umass*Mass(Kion(ia))
        velocity(:,ia) = velocity(:,ia) + force(:,ia)/mass_au * dt_h
        Ework = Ework - sum(aforce(:,ia)*(Rion(:,ia)-Rion_prev(:,ia)))
     enddo

     !NHC act on velocity with dt/2
     if(ensemble=="NVT" .and. thermostat=="nose-hoover")then
        call apply_nose_hoover_velocity(dt_h)
     endif

     !remove system momentum
     if(stop_system_momt=='y') call remove_system_momentum(0)

     call cal_Tion_Temperature_ion(Tion,Temperature_ion,velocity)
     Vion_gs  = Eall
     Eall     = Eall + Tion
     Htot     = Eall + Enh
     Ework_integ_fdR_gs = Ework_integ_fdR_gs + Ework

     !---Write section---
     ! Export to file_rt_data
     call write_md_gs_data(it,Vion_gs,Tion,Temperature_ion,Eall,Eall00,Ework_integ_fdR_gs,Enh,Htot)

     ! Export to standard log file
     if(comm_is_root(nproc_id_global)) then
        write(*,120) it*dt,Tion,Vion_gs,Ework_integ_fdR_gs,Enh,Eall,Eall-Eall00,Htot,Temperature_ion
120     format(1x,f10.4, 7e20.10E3,f12.3)
     endif

     ! Export to file_trj
     if (out_rvf_rt=='y' .and. mod(it,out_rvf_rt_step)==0)then
        write(comment_line,110) it, it*dt
110     format("#md-gs  step=",i8,"   time",e16.6)
        if(ensemble=="NVT" .and. thermostat=="nose-hoover") &
        &  write(comment_line,112) trim(comment_line), xi_nh
112     format(a,"  xi_nh=",e18.10)
        call write_xyz(comment_line,"add","rvf")
     endif

     ! Export electronic density (cube or vtk)
     if(out_dns_rt=='y' .and. mod(it,out_dns_rt_step)==0) then
        call write_density(it,'rt')
     end if

  enddo  !main loop of it

  call print_restart_data_md_gs
  call comm_sync_all

  if(comm_is_root(nproc_id_global)) write(*,*) 'This calculation is shutdown successfully!'

end subroutine calc_md_ground_state

  subroutine apply_nose_hoover_velocity(dt_h) !assuming half dt
    implicit none
    integer :: ia
    real(8) :: tmp_exp,dt_h

    tmp_exp = exp(-xi_nh*dt_h)
    do ia=1,NI
       velocity(:,ia) = velocity(:,ia) * tmp_exp
    enddo

    return
  end subroutine apply_nose_hoover_velocity

  subroutine apply_nose_hoover_thermostat(Temp_ion,dt_f) !assuming full dt
    implicit none
    real(8) :: Temp_ion,dt_f

    xi_nh = xi_nh + (Temp_ion/temperature0_ion-1d0)/(thermostat_tau**2d0)*dt_f

    return
  end subroutine apply_nose_hoover_thermostat

  subroutine write_md_gs_data(it,Vion_gs_t,Tion_gs_t,Temp_ion_t,Eall_gs_t,Eall_gs0,Ework_gs_t,Enh_t,Htot_gs_t)
    use inputoutput, only: t_unit_time,t_unit_energy
    implicit none
    integer :: fh_rt=202, it
    real(8) :: Vion_gs_t,Tion_gs_t,Temp_ion_t,Eall_gs_t,Ework_gs_t,Eall_gs0,Enh_t,Htot_gs_t

    if(comm_is_root(nproc_id_global)) then
      if(it==0) then
         open(fh_rt,file=trim(file_rt_data),status="unknown")

         write(fh_rt, '("#",1X,A)') "Real time calculation (Ground state MD)"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Tion", "Kinetic energy of ions"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Vion", "Potential energy of ions"
         write(fh_rt, '("#",1X,A,":",1X,A)') "E_work", "Work energy of ions(sum f*dr)"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Enh", "Energy of NH thermostat"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Eall", "Total energy"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Eall0", "Initial energy"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Hnvt", "Hamiltonian with NH thermostat"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Hnvt'", "Hnvt using E_work"
         write(fh_rt, '("#",1X,A,":",1X,A)') "Temperature_ion", "Temperature of ions"
         write(fh_rt, '("#",99(1X,I0,":",A,"[",A,"]"))') &
              &  1, "Time",   trim(t_unit_time%name), &
              &  2, "Tion",   trim(t_unit_energy%name), &
              &  3, "Vion",   trim(t_unit_energy%name), &
              &  4, "E_work", trim(t_unit_energy%name), &
              &  5, "Enh",    trim(t_unit_energy%name), &
              &  6, "Eall",   trim(t_unit_energy%name), &
              &  7, "Eall-Eall0", trim(t_unit_energy%name), &
              &  8, "Hnvt",   trim(t_unit_energy%name), &
              &  9, "Hnvt'",  trim(t_unit_energy%name), &
              & 10, "Temperature_ion", "K"
      endif

      write(fh_rt, "(F16.8,99(1X,E23.15E3))") &
          & it * dt    * t_unit_time%conv, &
          & Tion_gs_t  * t_unit_energy%conv, &
          & Vion_gs_t  * t_unit_energy%conv, &
          & Ework_gs_t * t_unit_energy%conv, &
          & Enh_t      * t_unit_energy%conv, &
          & Eall_gs_t  * t_unit_energy%conv, &
          & (Eall_gs_t-Eall_gs0) * t_unit_energy%conv, &
          & Htot_gs_t  * t_unit_energy%conv, &
          & (Tion_gs_t+Ework_gs_t+Enh_t) * t_unit_energy%conv, &
          & Temp_ion_t
      flush(fh_rt)

      if(it==Nt) close(fh_rt)
    end if

    call comm_sync_all
    return
  end subroutine write_md_gs_data

  subroutine cal_Tion_Temperature_ion(Ene_ion,Temp_ion, vel)
     implicit none
     integer :: ia
     real(8) :: Ene_ion,Temp_ion, vel(3,NI)

     Ene_ion=0.d0
     do ia=1,NI
        Ene_ion = Ene_ion + 0.5d0 * umass*Mass(Kion(ia)) * sum(vel(:,ia)**2d0)
     enddo
     Temp_ion = Ene_ion * 2d0 / (3d0*NI) / (kB/hartree2J)

     return
  end subroutine cal_Tion_Temperature_ion

  subroutine apply_velocity_scaling_ion(Temp_ion,vel)
     implicit none
     real(8) :: Temp_ion, vel(3,NI), fac_vscaling

     fac_vscaling = sqrt(temperature0_ion/Temp_ion)
     vel(:,:) = vel(:,:) * fac_vscaling

     return
  end subroutine apply_velocity_scaling_ion

  subroutine print_restart_data_md_gs
    use salmon_communication, only: comm_is_root
    use salmon_parallel, only: nproc_id_global
    use inputoutput, only: au_length_aa
    implicit none
    integer :: ia,unit_atomic_coor_tmp=201,ik,j
    real(8) :: tmpr(3), uconv
    character(100)  :: char_atom

    if(comm_is_root(nproc_id_global)) then

       select case(iflag_atom_coor)
       case(ntype_atom_coor_cartesian)
          open(unit_atomic_coor_tmp,file='.atomic_coor.tmp',status='old')
       case(ntype_atom_coor_reduced)
          open(unit_atomic_coor_tmp,file='.atomic_red_coor.tmp',status='old')
       end select

       if(unit_length=='AA')then
          uconv = au_length_aa
       else  !au
          uconv = 1d0
       endif

       write(*,*)
       write(*,9000) "##------ Restarting Data for MD-GS -------"
       write(*,9000) "# Coordinate in specified unit in input"
       write(*,9000) "# Copy into new input file"
       write(*,9000) "&atomic_coor"
       do ia = 1,NI
          read(unit_atomic_coor_tmp,*) char_atom, (tmpr(j),j=1,3),ik
          write(*,7000) trim(char_atom), Rion(1:3,ia)*uconv, ik
7000      format("     '",a,"'  ",3f18.10,i4)
       enddo
       write(*,9000) "/"

       write(*,*)
       write(*,9000) "# Velocity (atoms, thermostat if nose-hoover ooption) in [au]"
       write(*,9000) "# Copy to separated file used in file_ini_vel option with set_ini_velocity='r'"
       do ia = 1,NI
          write(*,8000) velocity(1:3,ia)
       enddo
8000   format(3e18.10)
       if(ensemble=="NVT" .and. thermostat=="nose-hoover") then
          write(*,8000) xi_nh
       endif
       write(*,*)

       close(unit_atomic_coor_tmp)
9000   format(a)
    endif

    return
  end subroutine print_restart_data_md_gs

  subroutine remove_system_momentum(flag_print_check)
    ! remove center of mass and momentum of whole system
    use salmon_communication, only: comm_is_root
    use salmon_parallel, only: nproc_id_global
    implicit none
    integer :: ia, flag_print_check
    real(8) :: v_com(3), sum_mass, mass_au

    !velocity of center of mass is removed
    v_com(:)=0d0
    sum_mass=0d0
    do ia=1,NI
       mass_au = umass*Mass(Kion(ia))
       v_com(:) = v_com(:) + mass_au * velocity(:,ia)
       sum_mass = sum_mass + mass_au
    enddo
    v_com(:) = v_com(:)/sum_mass
    do ia=1,NI
       velocity(:,ia) = velocity(:,ia) - v_com(:)
    enddo

    !rotation of system is removed
    !(this is only for isolated system)--- do nothing

    !(check velocity of center of mass)
    if(flag_print_check==1) then

       v_com(:)=0d0
       do ia=1,NI
          v_com(:) = v_com(:) + umass*Mass(Kion(ia)) * velocity(:,ia)
       enddo
       v_com(:) = v_com(:)/sum_mass
       if (comm_is_root(nproc_id_global)) &
       write(*,*) "    v_com =",real(v_com(:))

    endif

  end subroutine remove_system_momentum

  subroutine write_xyz(comment,action,rvf)
  ! Write xyz in xyz format but also velocity and force are printed if necessary
  ! (these can be used for restart of opt and md)
    use Global_Variables
    use inputoutput, only: au_length_aa
    use salmon_global, only: SYSname,iflag_atom_coor,ntype_atom_coor_cartesian,ntype_atom_coor_reduced
    use salmon_parallel, only: nproc_id_global
    use salmon_communication, only: comm_is_root
    implicit none
    integer :: ia,unit_xyz=200
    character(3) :: action,rvf
    character(1024) :: file_trj
    character(*) :: comment

    if(.not. comm_is_root(nproc_id_global)) return

    if(action=='new') then

       file_trj=trim(SYSname)//'_trj.xyz'
       open(unit_xyz,file=trim(file_trj),status="unknown")

    else if(action=='add') then

       write(unit_xyz,*) NI
       write(unit_xyz,*) trim(comment)
       do ia=1,NI
          if(      rvf=="r  " ) then
             write(unit_xyz,100) trim(atom_name(ia)),Rion(1:3,ia)*au_length_aa
          else if( rvf=="rv " ) then
             write(unit_xyz,110) trim(atom_name(ia)),Rion(1:3,ia)*au_length_aa,velocity(1:3,ia)
          else if( rvf=="rvf" ) then
             write(unit_xyz,120) trim(atom_name(ia)),Rion(1:3,ia)*au_length_aa,velocity(1:3,ia),force(1:3,ia)
          endif
       enddo

    else if(action=='end') then
       close(unit_xyz)
    endif

100 format(a2,3f18.10)
110 format(a2,3f18.10, "  #v=",3f18.10)
120 format(a2,3f18.10, "  #v=",3f18.10, "  #f=",3f18.10)

  end subroutine write_xyz

end module md_ground_state

max_stars_repo_name: shunsuke-sato/SALMON

max_issues_repo_name: shunsuke-sato/SALMON

max_forks_repo_name: shunsuke-sato/SALMON


text:
David Berman went to Davis Senior High School DHS, and currently attends http://en.wikipedia.org/wiki/UCLA UCLA, one of the most amazing universities on this Earth.
Overlooking all of these faults, he tends to be somewhat amusingly sarcastic.

This particularly unattractive photo of this alleged former Madrigal Choir Madrigal was taken in Wales.

ext: f

max_stars_repo_name: voflo/Search

max_issues_repo_name: voflo/Search

max_forks_repo_name: voflo/Search


module mod_types
  ! Andrew Valentine
  ! Universiteit Utrect
  ! a.p.valentine@uu.nl
  ! 2008-2014

  ! Module providing various compound types and functions to operate upon them
  ! Incorporates code from earlier modules:
  ! - datetime
  ! - ah

  use mod_log
  implicit none
  private

  public :: k_i0,k_rs,k_rd,date,julian_date,location,moment_tensor,assignment(=),operator(+),operator(-)
  public :: pi, twopi, sqrt2, euler, pi_d, twopi_d, fourpi_d, bigg

  integer,parameter :: k_i0=kind(1)
  integer,parameter :: k_rs=kind(1.0)  ! single precision
  integer,parameter :: k_rd=kind(1.d0) ! double precision

  ! Numerical constants (based on Numerical Recipes via DA-A)
  real(k_rs), parameter :: pi=3.141592653589793238462643383279502884197_k_rs
  real(k_rs), parameter :: twopi=6.283185307179586476925286766559005768394_k_rs
  real(k_rs), parameter :: sqrt2=1.41421356237309504880168872420969807856967_k_rs
  real(k_rs), parameter :: euler=0.5772156649015328606065120900824024310422_k_rs
  real(k_rd), parameter :: pi_d=3.141592653589793238462643383279502884197_k_rd
  real(k_rd), parameter :: twopi_d=6.283185307179586476925286766559005768394_k_rd
  real(k_rd), parameter :: fourpi_d=12.56637061_k_rd
  real(k_rd), parameter :: bigg=6.6723e-11_k_rd ! this was changed from 6.673e-11_k_rd
  !real(k_rd), parameter :: bigg=6.6723e-21_k_rd ! this was changed from 6.673e-11_k_rd

  type date
    sequence
    integer(k_i0) :: year
    integer(k_i0) :: month
    integer(k_i0) :: day
    integer(k_i0) :: hour
    integer(k_i0) :: minute
    real(k_rs)    :: second
  end type date

  type julian_date
    sequence
    integer(k_i0) :: year
    integer(k_i0) :: day
    integer(k_i0) :: hour
    integer(k_i0) :: minute
    real(k_rs)    :: second
  end type julian_date

  type location
    sequence
    real(k_rs) :: lat
    real(k_rs) :: lon
    real(k_rs) :: z
  end type location

  type moment_tensor
    sequence
    real(k_rs) :: mrr
    real(k_rs) :: mtt
    real(k_rs) :: mpp
    real(k_rs) :: mrt
    real(k_rs) :: mrp
    real(k_rs) :: mtp
    integer(k_i0) :: exponent
  end type moment_tensor

  interface assignment(=)
    module procedure d2j !date to julian
    module procedure j2d !julian to date
  end interface

  interface operator(+) !Add fixed number of seconds to a given date
    module procedure jpdbl !julian plus double
    module procedure jprea !julian plus single
    module procedure jpint !julian plus integer
    module procedure dpdbl !date plus double
    module procedure dprea !etc
    module procedure dpint
  end interface

  interface operator(-) ! Subtract seconds from date; compute seconds between two dates
    module procedure jmdbl !julian minus double
    module procedure jmrea !etc
    module procedure jmint
    module procedure jmj !julian minus julian
    module procedure jmd !etc
    module procedure dmdbl
    module procedure dmrea
    module procedure dmint
    module procedure dmd
    module procedure dmj
  end interface

  integer,dimension(12,0:1),parameter :: daylist=reshape((/0,31,59,90,120,151,181,212,243,273,304,334 &
                                                          ,0,31,60,91,121,152,182,213,244,274,305,335/),(/12,2/))

contains
  function lpyr(year) result(l)
    !Test whether 'year' is  leap year.
    !A year is a leap year if: divisible by 400
    !                          divisible by 4 but not 100 (except as above)
    integer,intent(in) :: year
    integer ::l !We want integer not logical to allow simple lookup in 'daylist'

    if(mod(year,400).eq.0) then
      l=1
    else if(mod(year,4).eq.0.and.mod(year,100).ne.0) then
      l=1
    else
      l=0
    end if
  end function lpyr

  subroutine d2j(myjul,mydate)
    type(date),intent(in) :: mydate
    type(julian_date),intent(out) :: myjul

    myjul%hour=mydate%hour
    myjul%minute=mydate%minute
    myjul%second=mydate%second
    myjul%year=mydate%year
    myjul%day=mydate%day+daylist(mydate%month,lpyr(mydate%year))
  end subroutine d2j

  subroutine j2d(mydate,myjul)
    type(julian_date),intent(in):: myjul
    type(date),intent(out):: mydate

    integer :: isleap

    mydate%hour=myjul%hour
    mydate%minute=myjul%minute
    mydate%second=myjul%second
    mydate%year=myjul%year

    isleap=lpyr(myjul%year)
    mydate%month=min(12,1+myjul%day/30)
    mydate%day=myjul%day-daylist(mydate%month,isleap)
    if(mydate%day.le.0) then
      mydate%month=mydate%month-1
      mydate%day=myjul%day-daylist(mydate%month,isleap)
    end if
  end subroutine j2d

  function jpdbl(myjul,second) result(j)
    type(julian_date),intent(in) :: myjul
    type(julian_date)::j
    real(k_rd),intent(in):: second

    real(k_rd) :: jps
    integer::jpints,minuteadd,houradd,dayadd,iflp

    j%year =myjul%year
    j%day  =myjul%day
    j%hour=myjul%hour
    j%minute =myjul%minute
    jps=real(myjul%second,k_rd)+second

    if(jps.ge.60.d0) then
      jpints=int(jps)
      minuteadd=int(jpints/60)
      j%minute=j%minute+minuteadd
      j%second=real((jps-minuteadd*60.d0),k_rs)
      if(j%second.lt.0.) then
        j%second=j%second+60.
        j%minute=j%minute-1
      end if
      if(j%second.ge.60.) then
        j%second=j%second-60.
        j%minute=j%minute+1
      end if
      if(j%second.lt.0..or.j%second.ge.60.) call fatal("mod_types","jpdbl","Logic error: unexpected value in 'j%second'.")
      if(j%minute.ge.60) then
        houradd=int(j%minute/60)
        j%hour=j%hour+houradd
        j%minute=j%minute-(houradd*60)
        if(j%hour.ge.24) then
          dayadd=int(j%hour/24)
          j%day=j%day+dayadd
          j%hour=j%hour-(24*dayadd)
          fixyr: do
            iflp=lpyr(j%year)
            if(j%day.gt.365+iflp) then
              j%day=j%day-365-iflp
              j%year=j%year+1
              cycle fixyr
            end if
            exit fixyr
          end do fixyr
        end if
      end if
    else if(jps.lt.0.d0) then
      minuteadd=int((-jps/60.d0)+1.d0)
      j%minute=j%minute-minuteadd
      j%second=real((jps+(real(minuteadd,k_rd)*60.d0)),k_rs)
      if(j%second.lt.0.) then
        j%second=j%second+60.
        j%minute=j%minute-1
      end if
      if(j%second.ge.60.) then
        j%second=j%second-60.
        j%minute=j%minute+1
      end if
      if(j%second.lt.0..or.j%second.ge.60.) call fatal("mod_types","jpdbl","Logic error: Unexpected value in 'j%second'.")
      if(j%minute.lt.0) then
        houradd=int((-j%minute-1)/60)+1
        j%hour=j%hour-houradd
        j%minute=j%minute+(60*houradd)
        if(j%minute.lt.0.or.j%minute.ge.60) call fatal("mod_types","jpdbl","Logic error: Unexpected value in 'j%minute'.")
        if(j%hour.lt.0) then
          dayadd=(-j%hour-1)/24+1
          j%day=j%day-dayadd
          j%hour=j%hour+(24*dayadd)
          if(j%hour.lt.0.or.j%hour.ge.24) call fatal("mod_types","jpdbl","Logic error: Unexpected value in 'j%hour'.")
          fixyr2: do
            if(j%day.le.0) then
              j%year=j%year-1
              iflp=lpyr(j%year)
              j%day=j%day+365+iflp
              cycle fixyr2
            end if
            exit fixyr2
          end do fixyr2
        end if
      end if
    else
      j%second=real(jps,k_rs)
    end if
  end function jpdbl

  function jmdbl(myjul,second) result(j)
    type(julian_date),intent(in) :: myjul
    type(julian_date) :: j
    real(k_rd),intent(in) :: second
    j=jpdbl(myjul,-second)
  end function jmdbl

  function jprea(myjul,second) result(j)
    type(julian_date),intent(in) :: myjul
    type(julian_date) :: j
    real(k_rs),intent(in) :: second
    j=jpdbl(myjul,real(second,k_rd))
  end function jprea

  function jmrea(myjul,second) result(j)
    type(julian_date),intent(in) :: myjul
    type(julian_date):: j
    real(k_rs),intent(in) :: second
    j=jpdbl(myjul,-real(second,k_rd))
  end function jmrea

  function jpint(myjul,second) result(j)
    type(julian_date),intent(in) :: myjul
    type(julian_date) ::j
    integer,intent(in) :: second
    j=jpdbl(myjul,real(second,k_rd))
  end function jpint

  function jmint(myjul,second) result(j)
    type(julian_date),intent(in) :: myjul
    type(julian_date)::j
    integer,intent(in)          :: second
    j=jpdbl(myjul,-real(second,k_rd))
  end function jmint

  function jmj(jul1,jul2) result(s)
    !Returns jul1-jul2 in seconds
    type(julian_date),intent(in) :: jul1,jul2
    real(k_rd)::s
    integer :: min,daydiff,i,islp

    daydiff=0
    if(jul1%year.ne.jul2%year) then
      do i=min(jul1%year,jul2%year),max(jul1%year,jul2%year)
        islp=lpyr(i)
        daydiff=daydiff+365+islp
      end do
      if(jul1%year.lt.jul2%year) daydiff=-daydiff
    end if
    s=jul1%second-jul2%second+60.d0*(real(jul1%minute-jul2%minute,k_rd)+60.d0*(real(jul1%hour-jul2%hour,k_rd)+ &
    24.d0*real(daydiff+jul1%day-jul2%day,k_rd)))
  end function jmj

  function jmd(jul1,day2) result(s)
    type(julian_date),intent(in) :: jul1
    real(k_rd)::s
    type(date),intent(in) :: day2
    type(julian_date) :: jul2
    jul2=day2 ! Does type conversion
    s=jul1-jul2
  end function jmd

  function dpdbl(mydat,second) result(d)
    type(date),intent(in) :: mydat
    type(date)::d
    real(k_rd),intent(in) :: second
    type(julian_date) :: myjul
    myjul=mydat
    myjul=myjul+second
    d=myjul
  end function dpdbl

  function dmdbl(mydat,second) result(d)
    type(date),intent(in)   :: mydat
    type(date):: d
    real(k_rd),intent(in) :: second
    d=dpdbl(mydat,-second)
  end function dmdbl

  function dprea(mydat,second) result(d)
    type(date),intent(in) :: mydat
    type(date)::d
    real(k_rs),intent(in) :: second
    type(julian_date) :: myjul
    myjul=mydat
    myjul=myjul+second
    d=myjul
  end function dprea

  function dmrea(mydat,second) result(d)
    type(date),intent(in)   :: mydat
    type(date)::d
    real(k_rs),intent(in) :: second
    d=dpdbl(mydat,-real(second,k_rd))
  end function dmrea

  function dpint(mydat,second) result(d)
    type(date),intent(in) :: mydat
    type(date)::d
    integer,intent(in) :: second

    type(julian_date) :: myjul
    myjul=mydat
    myjul=myjul+second
    d=myjul
  end function dpint

  function dmint(mydat,second) result(d)
    type(date),intent(in)   :: mydat
    integer,intent(in)      :: second
    type(date)::d
    d=mydat+(-real(second,k_rd))
  end function dmint

  function dmd(day1,day2) result(s)
    type(date),intent(in) :: day1
    type(date),intent(in) :: day2
    real(k_rd)::s
    type(julian_date) :: jul1,jul2

    jul1=day1
    jul2=day2

    s=jul1-jul2
  end function dmd

  function dmj(day1,jul2) result(s)
    type(date),intent(in)        :: day1
    type(julian_date),intent(in) :: jul2
    real(k_rd)::s
    type(julian_date) :: jul1
    jul1=day1
    s=jul1-jul2
  end function dmj
end module mod_types


! MIT License
!
! Copyright (c) 2020 SHEMAT-Suite
!
! Permission is hereby granted, free of charge, to any person obtaining a copy
! of this software and associated documentation files (the "Software"), to deal
! in the Software without restriction, including without limitation the rights
! to use, copy, modify, merge, publish, distribute, sublicense, and/or sell
! copies of the Software, and to permit persons to whom the Software is
! furnished to do so, subject to the following conditions:
!
! The above copyright notice and this permission notice shall be included in all
! copies or substantial portions of the Software.
!
! THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS OR
! IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
! FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL THE
! AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
! LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING FROM,
! OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
! SOFTWARE.

!>    @brief returns the begin and the end of the argument of the key in the line string
!>    @param[in] key search string
!>    @param[in] line current character line (find here the key)
!>    @param[out] first begin index of the argument
!>    @param[out] last end index of the argument
      SUBROUTINE get_arg(key,line,first,last)
        IMPLICIT NONE
        character (len=*) :: key, line
        INTEGER first, last, locstr, cfirst
        EXTERNAL locstr, cfirst
        INTRINSIC achar

        last = 0
        first = locstr(line,key)
        IF (first==0) RETURN
!
        first = first + len(key)
        first = first - 1 + cfirst(line(first:len(line)))
!
        DO last = first, len(line) - 1
          IF (line(last+1:last+1)==' ' .OR. line(last+1:last+1)==',' &
            .OR. line(last+1:last+1)==';' .OR. &
            line(last+1:last+1)==':' .OR. line(last+1:last+1)=='=' &
            .OR. line(last+1:last+1)==achar(9)) RETURN
        END DO
        last = len(line)
!
        RETURN
      END


text:
module pnmio

implicit none

integer, parameter ::    &
    PNM_BW_ASCII    = 1, &
    PNM_GRAY_ASCII  = 2, &
    PNM_RGB_ASCII   = 3, &
    PNM_BW_BINARY   = 4, &
    PNM_GRAY_BINARY = 5, &
    PNM_RGB_BINARY  = 6

contains

!=======================================================================

integer function writepnm(frm, b, fname, header)

! frm
!     One of the 6 parameters above indicating black and white, grayscale, or
!     RGB, as well as ASCII or binary
!
! b(:,:)
!     Rank-two byte array of pixel values.  The RGB versions have the RGB
!     triplets stretched out along the first array dimension.  The binary B&W
!     version has 8 bits of B&W pixels packed into a single byte of the array
!     along the first dimension.
!
! fname
!     File name without extension

character :: fname*(*), dat*256
character(len = :), allocatable :: str, frmstr

integer :: frm, iunit
character, allocatable :: b(:,:)

logical, optional, intent(in) :: header
logical :: headerl

writepnm = -1

if (present(header)) then
  headerl = header
else
  headerl = .true.
end if

if (frm == PNM_BW_ASCII) then

  ! B & W (0 and 1)
  open(newunit = iunit, file = trim(fname)//'.pbm')
  write(iunit, '(a)') 'P1'

else if (frm == PNM_GRAY_ASCII) then

  ! Grayscale (0 - 255)
  open(newunit = iunit, file = trim(fname)//'.pgm')
  write(iunit, '(a)') 'P2'

else if (frm == PNM_RGB_ASCII) then

  ! RGB (0 - 255 triplets)
  open(newunit = iunit, file = trim(fname)//'.ppm')
  write(iunit, '(a)') 'P3'

else if (frm == PNM_BW_BINARY) then

  ! B & W (0 and 1)
  open(newunit = iunit, file = trim(fname)//'.pbm', form = 'unformatted', access = 'stream')
  write(iunit) 'P4', char(10)

else if (frm == PNM_GRAY_BINARY) then

  ! Grayscale (0 - 255)
  open(newunit = iunit, file = trim(fname)//'.pgm', form = 'unformatted', access = 'stream')
  write(iunit) 'P5', char(10)

else if (frm == PNM_RGB_BINARY) then

  ! RGB (0 - 255 triplets)
  open(newunit = iunit, file = trim(fname)//'.ppm', form = 'unformatted', access = 'stream')
  write(iunit) 'P6', char(10)

else

  write(*,*)
  write(*,*) 'ERROR in writepnm'
  write(*,*) 'Invalid format'
  write(*,*) 'Valid values frm = 1 - 6'
  write(*,*)
  writepnm = -2
  return

end if

if (headerl) then
  ! This header string could be corrected to reflect the current function name,
  ! but it will require updated the "expected output" in parent repo tests.
  str = '# Pixel bit map in netpbm format, generated by writepbm'
  if (frm >= PNM_BW_ASCII .and. frm <= PNM_RGB_ASCII) then
    write(iunit, '(a)') str
  else
    write(iunit) str, char(10)
  end if
end if

if (frm >= PNM_BW_ASCII .and. frm <= PNM_RGB_ASCII) then
  if (frm == PNM_BW_ASCII .or. frm == PNM_GRAY_ASCII) then
    write(iunit, '(i0, a, i0)')  ubound(b, 1) - lbound(b, 1) + 1,      ' ', ubound(b, 2) - lbound(b, 2) + 1
  else
    write(iunit, '(i0, a, i0)') (ubound(b, 1) - lbound(b, 1) + 1) / 3, ' ', ubound(b, 2) - lbound(b, 2) + 1
  end if
else
  if (frm == PNM_BW_BINARY) then
    ! 8 bits packed into a byte
    write(dat, '(i0, a, i0)') (ubound(b, 1) - lbound(b, 1) + 1)*8,   ' ', ubound(b, 2) - lbound(b, 2) + 1
  else if (frm == PNM_GRAY_BINARY) then
    write(dat, '(i0, a, i0)')  ubound(b, 1) - lbound(b, 1) + 1,      ' ', ubound(b, 2) - lbound(b, 2) + 1
  else
    ! RGB triples
    write(dat, '(i0, a, i0)') (ubound(b, 1) - lbound(b, 1) + 1) / 3, ' ', ubound(b, 2) - lbound(b, 2) + 1
  end if
  write(iunit) trim(dat), char(10)
end if

if (frm >= PNM_BW_ASCII .and. frm <= PNM_RGB_ASCII) then
  if (frm == PNM_GRAY_ASCII .or. frm == PNM_RGB_ASCII) then
    write(iunit, '(i0)') 255
  end if
else
  if (frm == PNM_GRAY_BINARY .or. frm == PNM_RGB_BINARY) then
    write(dat, '(i0)') 255
    write(iunit) trim(dat), char(10)
  end if
end if

if (frm >= PNM_BW_ASCII .and. frm <= PNM_RGB_ASCII) then

  write(dat, '(a,i0,a)') '(', ubound(b,1) - lbound(b,1) + 1, '(i0,x))'
  frmstr = trim(dat)
  write(iunit, frmstr) ichar(b(:, ubound(b,2): lbound(b,2): -1))

else
  write(iunit) b(:, ubound(b,2): lbound(b,2): -1)
end if

close(iunit)

writepnm = 0

end function writepnm

!=======================================================================

end module pnmio

max_stars_repo_name: JeffIrwin/pnmio

max_issues_repo_name: JeffIrwin/pnmio

max_forks_repo_name: JeffIrwin/pnmio


text:
!=======================================================================
!
!  Analyse the formation/destruction routes of the specified species
!  that are dominant in the given particle by examining the reaction
!  network and all rates involved. Save a breakdown of these results
!  to the specified output file.
!
!-----------------------------------------------------------------------
SUBROUTINE ANALYSE_CHEMISTRY(PREFIX,SUFFIX,NPART,NSPEC,NREAC,TIME, &
                           & SPECIES,REACTANT,PRODUCT,PARTICLE, &
                           & AV_LIST,SPECIES_LIST)

   USE HEALPIX_TYPES
   USE PARTICLE_MODULE
   USE NUM2STR_FUNCTION
   USE GLOBAL_MODULE, ONLY : nH,nH2,nHe,nelect

   IMPLICIT NONE

   CHARACTER(LEN=*),    INTENT(IN) :: PREFIX,SUFFIX
   INTEGER(KIND=I4B),   INTENT(IN) :: NPART,NSPEC,NREAC
   REAL(KIND=DP),       INTENT(IN) :: TIME,AV_LIST(:)
   CHARACTER(LEN=*),    INTENT(IN) :: SPECIES(:),SPECIES_LIST(:)
   CHARACTER(LEN=*),    INTENT(IN) :: REACTANT(:,:),PRODUCT(:,:)
   TYPE(PARTICLE_TYPE), INTENT(IN) :: PARTICLE(:)

   INTEGER(KIND=I4B)  :: I,J,K,L,M,N,P,INDEX,IP,ID,II
   INTEGER(KIND=I4B)  :: NPR(1:NREAC),NDR(1:NREAC)
   INTEGER(KIND=I4B)  :: LENGTH(4),NMAX(1),MULTIPLE,TOTAL
   REAL(KIND=DP)      :: ABUNDANCE1,ABUNDANCE2,PERCENTAGE
   REAL(KIND=DP)      :: PRODUCTION_RATE(1:NREAC),TOTAL_PRODUCTION
   REAL(KIND=DP)      :: DESTRUCTION_RATE(1:NREAC),TOTAL_DESTRUCTION
   CHARACTER(LEN=20)  :: LHS,RHS
   CHARACTER(LEN=256) :: FILENAME

   INTEGER(KIND=I4B), SAVE :: UNIT=1

!  Specify the output file name
   FILENAME=TRIM(ADJUSTL(PREFIX))//TRIM(ADJUSTL('.chem'))//TRIM(ADJUSTL(SUFFIX))
   WRITE(6,*) 'Writing file [',TRIM(FILENAME),']'
   OPEN(UNIT=UNIT,FILE=FILENAME,STATUS='REPLACE')

   DO INDEX=LBOUND(AV_LIST,1),UBOUND(AV_LIST,1) ! Loop over list of visual extinctions to analyse

!  Find index P of the particle with AV closest to
!  the current entry in the visual extinction list
   P=1
   DO I=1,NPART
      IF(ABS(PARTICLE(I)%AV(0)-AV_LIST(INDEX)).LT.ABS(PARTICLE(P)%AV(0)-AV_LIST(INDEX))) P=I
   END DO

   DO L=LBOUND(SPECIES_LIST,1),UBOUND(SPECIES_LIST,1) ! Loop over list of species to analyse

!     Find index I of the species that corresponds
!     to the current entry L in the species list
      DO I=1,NSPEC
         IF(SPECIES(I).EQ.SPECIES_LIST(L)) EXIT
      END DO

!     Check that species I corresponds to the desired species
!     Continue to the next entry in the species list if not
      IF(SPECIES(I).NE.SPECIES_LIST(L)) CYCLE

      WRITE(UNIT,1) TRIM(NUM2STR(P)),TIME,PARTICLE(P)%AV(0),PARTICLE(P)%GAS_DENSITY, &
                  & TRIM(NUM2STR(PARTICLE(P)%GAS_TEMPERATURE,1)),TRIM(SPECIES(I))

!     Reset the formation/destruction rate counters
      IP=0
      ID=0

!     Reset the total formation/destruction rates
      TOTAL_PRODUCTION=0.0D0
      TOTAL_DESTRUCTION=0.0D0

!     Check all reactions to find relevant ones
      DO J=1,NREAC

!        Reset the reactant abundances
         ABUNDANCE1=0.0D0
         ABUNDANCE2=0.0D0

!-----------------------------------------------------------------------
!        Formation routes for species I
!-----------------------------------------------------------------------

         IF(ANY(PRODUCT(J,:).EQ.SPECIES(I))) THEN

            IP=IP+1

!           Store the reactant abundances
            DO N=1,NSPEC
               IF(SPECIES(N).EQ.REACTANT(J,1)) ABUNDANCE1=PARTICLE(P)%ABUNDANCE(N)
               IF(SPECIES(N).EQ.REACTANT(J,2)) ABUNDANCE2=PARTICLE(P)%ABUNDANCE(N)
            END DO

            MULTIPLE=0
            IF(PRODUCT(J,1).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1
            IF(PRODUCT(J,2).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1
            IF(PRODUCT(J,3).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1
            IF(PRODUCT(J,4).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1

!           Calculate the reaction rate
            IF(REACTANT(J,2).EQ."PHOTON" .OR. &
             & REACTANT(J,2).EQ."CRP   " .OR. &
             & REACTANT(J,2).EQ."CRPHOT" .OR. &
             & REACTANT(J,2).EQ."XRAY  " .OR. &
             & REACTANT(J,2).EQ."FREEZE" .OR. &
             & REACTANT(J,2).EQ."CRH   " .OR. &
             & REACTANT(J,2).EQ."PHOTD " .OR. &
             & REACTANT(J,2).EQ."THERM ") THEN
               PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1*MULTIPLE
            ELSE IF(REACTANT(J,2).EQ."XRSEC ") THEN
               IF(REACTANT(J,1).EQ."H2 ") THEN
                  PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(41.9D0*(1.0D0+6.72D0*(1.83D0*PARTICLE(P)%ABUNDANCE(nelect)/(1.0D0+0.83D0*PARTICLE(P)%ABUNDANCE(nelect)))**0.824D0)*(PARTICLE(P)%ABUNDANCE(nH2)+0.53D0*PARTICLE(P)%ABUNDANCE(nH)))
               ELSE IF(REACTANT(J,1).EQ."He ") THEN
                  PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(487.D0*(1.0D0+12.5D0*PARTICLE(P)%ABUNDANCE(nelect)**0.994D0)*(PARTICLE(P)%ABUNDANCE(nHe)))
               ELSE
                  PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(39.8D0*(1.0D0+12.2D0*PARTICLE(P)%ABUNDANCE(nelect)**0.866D0)*(PARTICLE(P)%ABUNDANCE(nH)+1.89D0*PARTICLE(P)%ABUNDANCE(nH2)))
               END IF
            ELSE IF(REACTANT(J,2).EQ."XRLYA ") THEN
               PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(39.8D0*(1.0D0+12.2D0*PARTICLE(P)%ABUNDANCE(nelect)**0.866D0)*(PARTICLE(P)%ABUNDANCE(nH)+1.89D0*PARTICLE(P)%ABUNDANCE(nH2)))
            ELSE IF(REACTANT(J,2).EQ."XRPHOT") THEN
               PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(41.9D0*(1.0D0+6.72D0*(1.83D0*PARTICLE(P)%ABUNDANCE(nelect)/(1.0D0+0.83D0*PARTICLE(P)%ABUNDANCE(nelect)))**0.824D0)*(PARTICLE(P)%ABUNDANCE(nH2)+0.53D0*PARTICLE(P)%ABUNDANCE(nH)))
            ELSE
               PRODUCTION_RATE(IP)=PARTICLE(P)%RATE(J)*ABUNDANCE1*ABUNDANCE2*PARTICLE(P)%GAS_DENSITY*MULTIPLE
            END IF

            TOTAL_PRODUCTION=TOTAL_PRODUCTION+PRODUCTION_RATE(IP)
            NPR(IP)=J
         END IF

!-----------------------------------------------------------------------
!        Destruction routes for species I
!-----------------------------------------------------------------------

         IF(ANY(REACTANT(J,:).EQ.SPECIES(I))) THEN

            ID=ID+1

!           Store the reactant abundances
            DO N=1,NSPEC
               IF(SPECIES(N).EQ.REACTANT(J,1)) ABUNDANCE1=PARTICLE(P)%ABUNDANCE(N)
               IF(SPECIES(N).EQ.REACTANT(J,2)) ABUNDANCE2=PARTICLE(P)%ABUNDANCE(N)
            END DO

            MULTIPLE=0
            IF(REACTANT(J,1).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1
            IF(REACTANT(J,2).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1
            IF(REACTANT(J,3).EQ.SPECIES(I)) MULTIPLE=MULTIPLE+1

!           Calculate the reaction rate
            IF(REACTANT(J,2).EQ."PHOTON" .OR. &
             & REACTANT(J,2).EQ."CRP   " .OR. &
             & REACTANT(J,2).EQ."CRPHOT" .OR. &
             & REACTANT(J,2).EQ."XRAY  " .OR. &
             & REACTANT(J,2).EQ."FREEZE" .OR. &
             & REACTANT(J,2).EQ."CRH   " .OR. &
             & REACTANT(J,2).EQ."PHOTD " .OR. &
             & REACTANT(J,2).EQ."THERM ") THEN
               DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1*MULTIPLE
            ELSE IF(REACTANT(J,2).EQ."XRSEC ") THEN
               IF(REACTANT(J,1).EQ."H2 ") THEN
                  DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(41.9D0*(1.0D0+6.72D0*(1.83D0*PARTICLE(P)%ABUNDANCE(nelect)/(1.0D0+0.83D0*PARTICLE(P)%ABUNDANCE(nelect)))**0.824D0)*(PARTICLE(P)%ABUNDANCE(nH2)+0.53D0*PARTICLE(P)%ABUNDANCE(nH)))
               ELSE IF(REACTANT(J,1).EQ."He ") THEN
                  DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(487.D0*(1.0D0+12.5D0*PARTICLE(P)%ABUNDANCE(nelect)**0.994D0)*(PARTICLE(P)%ABUNDANCE(nHe)))
               ELSE
                  DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(39.8D0*(1.0D0+12.2D0*PARTICLE(P)%ABUNDANCE(nelect)**0.866D0)*(PARTICLE(P)%ABUNDANCE(nH)+1.89D0*PARTICLE(P)%ABUNDANCE(nH2)))
               END IF
            ELSE IF(REACTANT(J,2).EQ."XRLYA ") THEN
               DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(39.8D0*(1.0D0+12.2D0*PARTICLE(P)%ABUNDANCE(nelect)**0.866D0)*(PARTICLE(P)%ABUNDANCE(nH)+1.89D0*PARTICLE(P)%ABUNDANCE(nH2)))
            ELSE IF(REACTANT(J,2).EQ."XRPHOT") THEN
               DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1/(41.9D0*(1.0D0+6.72D0*(1.83D0*PARTICLE(P)%ABUNDANCE(nelect)/(1.0D0+0.83D0*PARTICLE(P)%ABUNDANCE(nelect)))**0.824D0)*(PARTICLE(P)%ABUNDANCE(nH2)+0.53D0*PARTICLE(P)%ABUNDANCE(nH)))
            ELSE
               DESTRUCTION_RATE(ID)=PARTICLE(P)%RATE(J)*ABUNDANCE1*ABUNDANCE2*PARTICLE(P)%GAS_DENSITY*MULTIPLE
            END IF

            TOTAL_DESTRUCTION=TOTAL_DESTRUCTION+DESTRUCTION_RATE(ID)
            NDR(ID)=J
         END IF

      END DO ! End of loop over reactions

!     Prevent divide-by-zero errors by setting minimum finite
!     values for the total formation and destruction rates
      IF(TOTAL_PRODUCTION.LT.1.0D-99) TOTAL_PRODUCTION=1.0D-99
      IF(TOTAL_DESTRUCTION.LT.1.0D-99) TOTAL_DESTRUCTION=1.0D-99

!-----------------------------------------------------------------------
!     Find the lengths of the longest reactant and product names
!-----------------------------------------------------------------------

      LENGTH=0
      DO N=1,IP
!        Only consider reactions that represent > 0.5% of the total rate
         IF(1.0D2*(PRODUCTION_RATE(N)/TOTAL_PRODUCTION).GT.0.5D0) THEN
            LENGTH(1)=MAX(LENGTH(1),LEN_TRIM(REACTANT(NPR(N),1)),LEN_TRIM(SPECIES(I)))
            LENGTH(2)=MAX(LENGTH(2),LEN_TRIM(TRIM(ADJUSTL(REACTANT(NPR(N),2)))//'   '// &
                                           & TRIM(ADJUSTL(REACTANT(NPR(N),3)))))
         END IF
      END DO
      DO N=1,ID
!        Only consider reactions that represent > 0.5% of the total rate
         IF(1.0D2*(DESTRUCTION_RATE(N)/TOTAL_DESTRUCTION).GT.0.5D0) THEN
            LENGTH(2)=MAX(LENGTH(2),LEN_TRIM(TRIM(ADJUSTL(REACTANT(NDR(N),2)))//'   '// &
                                           & TRIM(ADJUSTL(REACTANT(NDR(N),3)))))
            LENGTH(3)=MAX(LENGTH(3),LEN_TRIM(PRODUCT(NDR(N),1)),LEN_TRIM(SPECIES(I)))
         END IF
      END DO

!-----------------------------------------------------------------------
!     Output the formation reactions and their rates
!-----------------------------------------------------------------------

      NMAX=0
      TOTAL=0

!     List the formation reactions in order of decreasing importance
      DO M=1,IP

!        Find the location of the maximum value
         NMAX=MAXLOC(PRODUCTION_RATE(1:IP))
         N=NMAX(1)

!        Exit the loop once the reaction rates reach zero
         IF(PRODUCTION_RATE(N).EQ.0.0D0) EXIT

!        Calculate the percentage of the total formation
!        rate that is contributed by the current reaction
         PERCENTAGE=1.0D2*(PRODUCTION_RATE(N)/TOTAL_PRODUCTION)
         IF(PERCENTAGE.LT.0.5D0) PERCENTAGE=1.0D0
         TOTAL=TOTAL+NINT(PERCENTAGE)

!        Find the first instance of the species of interest in the list of products
!         NMAX=MAXLOC(PRODUCT(NPR(N),:),(PRODUCT(NPR(N),:).EQ.SPECIES(I)))
         DO II=1,4
            IF (PRODUCT(NPR(N),II) .EQ. SPECIES(I)) THEN
               NMAX = II
               EXIT
            END IF
         END DO
         K=NMAX(1)

!        Create the left- and right-hand sides of the reaction string
         LHS=REACTANT(NPR(N),1)
         DO J=2,3 ! Loop over reactants
            IF(REACTANT(NPR(N),J).NE." ") THEN
               IF(LEN_TRIM(LHS).EQ.LEN_TRIM(REACTANT(NPR(N),1))) THEN
                  LHS=LHS(1:LENGTH(1))//' + '//REACTANT(NPR(N),J)
               ELSE
                  LHS=TRIM(ADJUSTL(LHS))//' + '//REACTANT(NPR(N),J)
               END IF
            END IF
         END DO ! End of loop over reactants

         RHS=PRODUCT(NPR(N),K)
         DO J=1,4 ! Loop over products
            IF(J.EQ.K) CYCLE ! Skip this one since it has already been added to the string
            IF(PRODUCT(NPR(N),J).NE." ") THEN
               IF(LEN_TRIM(RHS).EQ.LEN_TRIM(PRODUCT(NPR(N),K))) THEN
                  RHS=RHS(1:LENGTH(3))//' + '//PRODUCT(NPR(N),J)
               ELSE
                  RHS=TRIM(ADJUSTL(RHS))//' + '//PRODUCT(NPR(N),J)
               END IF
            END IF
         END DO ! End of loop over products

!        Right-align the left-hand side of the reaction string
         LHS=REPEAT(' ',LEN(LHS)-LENGTH(1)-LENGTH(2)-3)//ADJUSTL(LHS)

!        Left-align the right-hand side of the reaction string
         RHS=ADJUSTL(RHS)

!        Write the reaction and its percentage of the total formation rate
         WRITE(UNIT,2) LHS,RHS,NINT(PERCENTAGE)

!        Exit the loop once the sum of the reaction rates reaches 100%
         IF(TOTAL.GE.100) EXIT

!        Set the formation rate of this reaction to zero
!        to prevent it from being included more than once
         PRODUCTION_RATE(N)=0.0D0

      END DO ! End of loop over formation reactions
      WRITE(UNIT,*)

!-----------------------------------------------------------------------
!     Output the destruction reactions and their rates
!-----------------------------------------------------------------------

      NMAX=0
      TOTAL=0

!     List the destruction reactions in order of decreasing importance
      DO M=1,ID

!        Find the location of the maximum value
         NMAX=MAXLOC(DESTRUCTION_RATE(1:ID))
         N=NMAX(1)

!        Exit the loop once the reaction rates reach zero
         IF(DESTRUCTION_RATE(N).EQ.0.0D0) EXIT

!        Calculate the percentage of the total destruction
!        rate that is contributed by the current reaction
         PERCENTAGE=1.0D2*(DESTRUCTION_RATE(N)/TOTAL_DESTRUCTION)
         IF(PERCENTAGE.LT.0.5D0) PERCENTAGE=1.0D0
         TOTAL=TOTAL+NINT(PERCENTAGE)

!        Find the first instance of the species of interest in the list of reactants
!         NMAX=MAXLOC(REACTANT(NDR(N),:),(REACTANT(NDR(N),:).EQ.SPECIES(I)))
         DO II=1,4
            IF (REACTANT(NDR(N),II) .EQ. SPECIES(I)) THEN
               NMAX = II
               EXIT
            END IF
         END DO
         K=NMAX(1)

!        Create the left- and right-hand sides of the reaction string
         LHS=REACTANT(NDR(N),K)
         DO J=1,3 ! Loop over reactants
            IF(J.EQ.K) CYCLE ! Skip this one since it has already been added to the string
            IF(REACTANT(NDR(N),J).NE." ") THEN
               IF(LEN_TRIM(LHS).EQ.LEN_TRIM(REACTANT(NDR(N),K))) THEN
                  LHS=LHS(1:LENGTH(1))//' + '//REACTANT(NDR(N),J)
               ELSE
                  LHS=TRIM(ADJUSTL(LHS))//' + '//REACTANT(NDR(N),J)
               END IF
            END IF
         END DO ! End of loop over reactants

         RHS=PRODUCT(NDR(N),1)
         DO J=2,4 ! Loop over products
            IF(PRODUCT(NDR(N),J).NE." ") THEN
               IF(LEN_TRIM(RHS).EQ.LEN_TRIM(PRODUCT(NDR(N),1))) THEN
                  RHS=RHS(1:LENGTH(3))//' + '//PRODUCT(NDR(N),J)
               ELSE
                  RHS=TRIM(ADJUSTL(RHS))//' + '//PRODUCT(NDR(N),J)
               END IF
            END IF
         END DO ! End of loop over products

!        Right-align the left-hand side of the reaction string
         LHS=REPEAT(' ',LEN(LHS)-LENGTH(1)-LENGTH(2)-3)//ADJUSTL(LHS)

!        Left-align the right-hand side of the reaction string
         RHS=ADJUSTL(RHS)

!        Write the reaction and its percentage of the total destruction rate
         WRITE(UNIT,2) LHS,RHS,-NINT(PERCENTAGE)

!        Exit the loop once the sum of the reaction rates reaches 100%
         IF(TOTAL.GE.100) EXIT

!        Set the destruction rate of this reaction to zero
!        to prevent it from being included more than once
         DESTRUCTION_RATE(N)=0.0D0

      END DO ! End of loop over destruction reactions
      WRITE(UNIT,*)

!-----------------------------------------------------------------------
!     Output the abundance and total formation/destruction rates
!-----------------------------------------------------------------------

      WRITE(UNIT,3) PARTICLE(P)%ABUNDANCE(I),TOTAL_PRODUCTION,TOTAL_DESTRUCTION
      IF(L.LT.UBOUND(SPECIES_LIST,1)) WRITE(UNIT,4)

   END DO ! End of loop over species list
   WRITE(UNIT,5)

   END DO ! End of loop over visual extinction list
   CLOSE(UNIT)

 1 FORMAT('Particle = ',A,', t =',ES7.1E1,' yr, A_V =',ES7.1E1,' mag',/, &
        & 'n_H =',ES7.1E1,' cm^-3, T_gas = ',A,' K',/, &
        & 'Species = ',A,/)
 2 FORMAT(A20,' --> ',A20,'Rate:',SP,I4,'%')
 3 FORMAT('Total Abundance  =',ES10.3,/, &
        & 'Formation Rate   =',ES10.3,' s^-1',/, &
        & 'Destruction Rate =',ES10.3,' s^-1',/)
 4 FORMAT(80('-'),/)
 5 FORMAT(80('='),/)

   RETURN
END SUBROUTINE ANALYSE_CHEMISTRY
!=======================================================================

max_stars_repo_name: uclchem/uclpdr

max_issues_repo_name: uclchem/uclpdr

max_forks_repo_name: uclchem/uclpdr


text:
      subroutine blkind(dvar,cvar,cfmt)
      implicit none
c
      double precision dvar
      character        cvar*6,cfmt*(*)
c
      integer       lp
      common/linepr/lp
c
c     read in one d.p. value from stdin
c
      character*6 cvarin
c
      read(*,*) dvar,cvarin
      write(lp,cfmt) cvarin,dvar
      call flush(lp)
c
      if     (cvar.ne.cvarin) then
        write(lp,*)
        write(lp,*) 'error in blkind - input ',cvarin,
     +                      ' but should be ',cvar
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
      end
      subroutine blkinr(rvar,cvar,cfmt)
      implicit none
c
      real      rvar
      character cvar*6,cfmt*(*)
c
      integer       lp
      common/linepr/lp
c
c     read in one real value from stdin
c
      character*6 cvarin
c
      read(*,*) rvar,cvarin
      write(lp,cfmt) cvarin,rvar
      call flush(lp)
c
      if     (cvar.ne.cvarin) then
        write(lp,*)
        write(lp,*) 'error in blkinr - input ',cvarin,
     +                      ' but should be ',cvar
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
      end
      subroutine blkinr2(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2)
      implicit none
c
      real      rvar
      integer   nvar
      character cvar1*6,cvar2*6,cfmt1*(*),cfmt2*(*)
c
      integer       lp
      common/linepr/lp
c
c     read in one real value from stdin,
c     identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2)
c
      character*6 cvarin
c
      read(*,*) rvar,cvarin
c
      if     (cvar1.eq.cvarin) then
        nvar = 1
        write(lp,cfmt1) cvarin,rvar
        call flush(lp)
      elseif (cvar2.eq.cvarin) then
        nvar = 2
        write(lp,cfmt2) cvarin,rvar
        call flush(lp)
      else
        write(lp,*)
        write(lp,*) 'error in blkinr2 - input ',cvarin,
     +                      ' but should be ',cvar1,' or ',cvar2
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
      end
      subroutine blkinr3(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3)
      implicit none
c
      real      rvar
      integer   nvar
      character cvar1*6,cvar2*6,cvar3*6,cfmt1*(*),cfmt2*(*),cfmt3*(*)
c
      integer       lp
      common/linepr/lp
c
c     read in one of three possible real values from stdin,
c     identified as either cvar1 (return nvar=1) or
c                          cvar2 (return nvar=2) or
c                          cvar3 (return nvar=3)
c
      call blkinr9(rvar,nvar,cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3,
     &             'XXXXXX','("blkinr: ",a6," =",f11.4," ")',
     &             'XXXXXX','("blkinr: ",a6," =",f11.4," ")',
     &             'XXXXXX','("blkinr: ",a6," =",f11.4," ")',
     &             'XXXXXX','("blkinr: ",a6," =",f11.4," ")',
     &             'XXXXXX','("blkinr: ",a6," =",f11.4," ")',
     &             'XXXXXX','("blkinr: ",a6," =",f11.4," ")')
      return
      end
      subroutine blkinr9(rvar,nvar,
     &                   cvar1,cfmt1,cvar2,cfmt2,cvar3,cfmt3,
     &                   cvar4,cfmt4,cvar5,cfmt5,cvar6,cfmt6,
     &                   cvar7,cfmt7,cvar8,cfmt8,cvar9,cfmt9)
      implicit none
c
      real      rvar
      integer   nvar
      character cvar1*6,cvar2*6,cvar3*6,cfmt1*(*),cfmt2*(*),cfmt3*(*)
      character cvar4*6,cvar5*6,cvar6*6,cfmt4*(*),cfmt5*(*),cfmt6*(*)
      character cvar7*6,cvar8*6,cvar9*6,cfmt7*(*),cfmt8*(*),cfmt9*(*)
c
      integer       lp
      common/linepr/lp
c
c     read in one of nine possible real values from stdin,
c     identified as either cvar1 (return nvar=1) or
c                          cvar2 (return nvar=2) or
c                          ...
c                          cvar9 (return nvar=9)
c
      character*6 cvarin
c
      read(*,*) rvar,cvarin
c
      if     (cvar1.eq.cvarin) then
        nvar = 1
        write(lp,cfmt1) cvarin,rvar
        call flush(lp)
      elseif (cvar2.eq.cvarin) then
        nvar = 2
        write(lp,cfmt2) cvarin,rvar
        call flush(lp)
      elseif (cvar3.eq.cvarin) then
        nvar = 3
        write(lp,cfmt3) cvarin,rvar
        call flush(lp)
      elseif (cvar4.eq.cvarin) then
        nvar = 4
        write(lp,cfmt4) cvarin,rvar
        call flush(lp)
      elseif (cvar5.eq.cvarin) then
        nvar = 5
        write(lp,cfmt5) cvarin,rvar
        call flush(lp)
      elseif (cvar6.eq.cvarin) then
        nvar = 6
        write(lp,cfmt6) cvarin,rvar
        call flush(lp)
      elseif (cvar7.eq.cvarin) then
        nvar = 7
        write(lp,cfmt7) cvarin,rvar
        call flush(lp)
      elseif (cvar8.eq.cvarin) then
        nvar = 8
        write(lp,cfmt8) cvarin,rvar
        call flush(lp)
      elseif (cvar9.eq.cvarin) then
        nvar = 9
        write(lp,cfmt9) cvarin,rvar
        call flush(lp)
      else
        write(lp,*)
        write(lp,*) 'error in blkinr9 - input ',cvarin,
     +              ' but should be:'
        write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or'
        write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or'
        write(lp,*) cvar7,' or ',cvar8,' or ',cvar9
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
      end
      subroutine blkini(ivar,cvar)
      implicit none
c
      integer     ivar
      character*6 cvar
c
      integer       lp
      common/linepr/lp
c
c     read in one integer value from stdin
c
      character*6 cvarin
c
      read(*,*) ivar,cvarin
      write(lp,6000) cvarin,ivar
      call flush(lp)
c
      if     (cvar.ne.cvarin) then
        write(lp,*)
        write(lp,*) 'error in blkini - input ',cvarin,
     +                      ' but should be ',cvar
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
 6000 format('blkini: ',a6,' =',i6)
      end
      subroutine blkini2(ivar,nvar,cvar1,cvar2)
      implicit none
c
      integer     ivar,nvar
      character*6 cvar1,cvar2
c
      integer       lp
      common/linepr/lp
c
c     read in one integer value from stdin
c     identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2)
c
      character*6 cvarin
c
      read(*,*) ivar,cvarin
      write(lp,6000) cvarin,ivar
      call flush(lp)
c
      if     (cvarin.eq.cvar1) then
        nvar = 1
      elseif (cvarin.eq.cvar2) then
        nvar = 2
      else
        write(lp,*)
        write(lp,*) 'error in blkini2 - input ',cvarin,
     +                      ' but should be ',cvar1,' or ',cvar2
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
 6000 format('blkini: ',a6,' =',i6)
      end
      subroutine blkini3(ivar,nvar,cvar1,cvar2,cvar3)
      implicit none
c
      integer     ivar,nvar
      character*6 cvar1,cvar2,cvar3
c
      integer       lp
      common/linepr/lp
c
c     read in one of three possible integer values from stdin
c     identified as either cvar1 (return nvar=1) or
c                          cvar2 (return nvar=2) or
c                          cvar3 (return nvar=3)
c
      call blkini9(ivar,nvar,cvar1,cvar2,cvar3,
     &                       'XXXXXX','XXXXXX','XXXXXX',
     &                       'XXXXXX','XXXXXX','XXXXXX')
      return
      end
      subroutine blkini4(ivar,nvar,cvar1,cvar2,cvar3,cvar4)
      implicit none
c
      integer     ivar,nvar
      character*6 cvar1,cvar2,cvar3,cvar4
c
      integer       lp
      common/linepr/lp
c
c     read in one of four possible integer values from stdin
c     identified as either cvar1 (return nvar=1) or
c                          cvar2 (return nvar=2) or
c                          cvar3 (return nvar=3) or
c                          cvar4 (return nvar=4)
c
      call blkini9(ivar,nvar,cvar1,   cvar2,   cvar3,
     &                       cvar4,   'XXXXXX','XXXXXX',
     &                       'XXXXXX','XXXXXX','XXXXXX')
      return
      end
      subroutine blkini5(ivar,nvar,cvar1,cvar2,cvar3,cvar4,cvar5)
      implicit none
c
      integer     ivar,nvar
      character*6 cvar1,cvar2,cvar3,cvar4,cvar5
c
      integer       lp
      common/linepr/lp
c
c     read in one of five possible integer values from stdin
c     identified as either cvar1 (return nvar=1) or
c                          cvar2 (return nvar=2) or
c                          cvar3 (return nvar=3) or
c                          cvar4 (return nvar=4) or
c                          cvar5 (return nvar=5)
c
      call blkini9(ivar,nvar,cvar1,   cvar2,   cvar3,
     &                       cvar4,   cvar5,   'XXXXXX',
     &                       'XXXXXX','XXXXXX','XXXXXX')
      return
      end
      subroutine blkini9(ivar,nvar,cvar1,cvar2,cvar3,
     &                             cvar4,cvar5,cvar6,
     &                             cvar7,cvar8,cvar9)
      implicit none
c
      integer     ivar,nvar
      character*6 cvar1,cvar2,cvar3
      character*6 cvar4,cvar5,cvar6
      character*6 cvar7,cvar8,cvar9
c
      integer       lp
      common/linepr/lp
c
c     read in one of nine possible integer values from stdin
c     identified as either cvar1 (return nvar=1) or
c                          cvar2 (return nvar=2) or
c                          ...
c                          cvar9 (return nvar=9)
c                          cvar3 (return nvar=3) or
c
      character*6 cvarin
c
      read(*,*) ivar,cvarin
      write(lp,6000) cvarin,ivar
      call flush(lp)
c
      if     (cvarin.eq.cvar1) then
        nvar = 1
      elseif (cvarin.eq.cvar2) then
        nvar = 2
      elseif (cvarin.eq.cvar3) then
        nvar = 3
      elseif (cvar4.eq.cvarin) then
        nvar = 4
      elseif (cvar5.eq.cvarin) then
        nvar = 5
      elseif (cvar6.eq.cvarin) then
        nvar = 6
      elseif (cvar7.eq.cvarin) then
        nvar = 7
      elseif (cvar8.eq.cvarin) then
        nvar = 8
      elseif (cvar9.eq.cvarin) then
        nvar = 9
      else
        write(lp,*)
        write(lp,*) 'error in blkini9 - input ',cvarin,
     +              ' but should be:'
        write(lp,*) cvar1,' or ',cvar2,' or ',cvar3,' or'
        write(lp,*) cvar4,' or ',cvar5,' or ',cvar6,' or'
        write(lp,*) cvar7,' or ',cvar8,' or ',cvar9
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
 6000 format('blkini: ',a6,' =',i6)
      end
      subroutine blkinl(lvar,cvar)
      implicit none
c
      logical     lvar
      character*6 cvar
c
      integer       lp
      common/linepr/lp
c
c     read in one logical value from stdin
c     due to a SGI bug for logical I/O: read in an integer 0=F,1=T
c
      character*6 cvarin
      integer     ivar
c
      read(*,*) ivar,cvarin
      lvar = ivar .ne. 0
      write(lp,6000) cvarin,lvar
      call flush(lp)
c
      if     (cvar.ne.cvarin) then
        write(lp,*)
        write(lp,*) 'error in blkinl - input ',cvarin,
     +                      ' but should be ',cvar
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
 6000 format('blkinl: ',a6,' =',l6)
      end
      subroutine blkinl_99(lvar,cvar)
      implicit none
c
      logical     lvar
      character*6 cvar
c
      integer       lp
      common/linepr/lp
c
c     read in one logical value from unit 99
c     due to a SGI bug for logical I/O: read in an integer 0=F,1=T
c
      character*6 cvarin
      integer     ivar
c
      read(99,*) ivar,cvarin
      lvar = ivar .ne. 0
      write(lp,6000) cvarin,lvar
      call flush(lp)
c
      if     (cvar.ne.cvarin) then
        write(lp,*)
        write(lp,*) 'error in blkinl_99 - input ',cvarin,
     +                      ' but should be ',cvar
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
 6000 format('blkinl: ',a6,' =',l6)
      end
      subroutine blkinl2_99(lvar,nvar,cvar1,cvar2)
      implicit none
c
      logical     lvar
      integer     nvar
      character*6 cvar1,cvar2
c
      integer       lp
      common/linepr/lp
c
c     read in one logical value from unit 99
c     identified as either cvar1 (return nvar=1) or cvar2 (return nvar=2)
c     due to a SGI bug for logical I/O: read in an integer 0=F,1=T
c
      character*6 cvarin
      integer     ivar
c
      read(99,*) ivar,cvarin
      lvar = ivar .ne. 0
      write(lp,6000) cvarin,lvar
      call flush(lp)
c
      if     (cvarin.eq.cvar1) then
        nvar = 1
      elseif (cvarin.eq.cvar2) then
        nvar = 2
      else
        write(lp,*)
        write(lp,*) 'error in blkinl2_99 - input ',cvarin,
     +                      ' but should be ',cvar1,' or ',cvar2
        write(lp,*)
        call flush(lp)
        stop
      endif
      return
 6000 format('blkinl: ',a6,' =',l6)
      end

max_stars_repo_name: TillRasmussen/HYCOM-tools

max_issues_repo_name: TillRasmussen/HYCOM-tools

max_forks_repo_name: TillRasmussen/HYCOM-tools


text:
*DECK ZASYI
      SUBROUTINE ZASYI (ZR, ZI, FNU, KODE, N, YR, YI, NZ, RL, TOL, ELIM,
     +   ALIM)
C***BEGIN PROLOGUE  ZASYI
C***SUBSIDIARY
C***PURPOSE  Subsidiary to ZBESI and ZBESK
C***LIBRARY   SLATEC
C***TYPE      ALL (CASYI-A, ZASYI-A)
C***AUTHOR  Amos, D. E., (SNL)
C***DESCRIPTION
C
C     ZASYI COMPUTES THE I BESSEL FUNCTION FOR REAL(Z).GE.0.0 BY
C     MEANS OF THE ASYMPTOTIC EXPANSION FOR LARGE ABS(Z) IN THE
C     REGION ABS(Z).GT.MAX(RL,FNU*FNU/2). NZ=0 IS A NORMAL RETURN.
C     NZ.LT.0 INDICATES AN OVERFLOW ON KODE=1.
C
C***SEE ALSO  ZBESI, ZBESK
C***ROUTINES CALLED  D1MACH, ZABS, ZDIV, ZEXP, ZMLT, ZSQRT
C***REVISION HISTORY  (YYMMDD)
C   830501  DATE WRITTEN
C   910415  Prologue converted to Version 4.0 format.  (BAB)
C   930122  Added ZEXP and ZSQRT to EXTERNAL statement.  (RWC)
C***END PROLOGUE  ZASYI
C     COMPLEX AK1,CK,CONE,CS1,CS2,CZ,CZERO,DK,EZ,P1,RZ,S2,Y,Z
      DOUBLE PRECISION AA, AEZ, AK, AK1I, AK1R, ALIM, ARG, ARM, ATOL,
     * AZ, BB, BK, CKI, CKR, CONEI, CONER, CS1I, CS1R, CS2I, CS2R, CZI,
     * CZR, DFNU, DKI, DKR, DNU2, ELIM, EZI, EZR, FDN, FNU, PI, P1I,
     * P1R, RAZ, RL, RTPI, RTR1, RZI, RZR, S, SGN, SQK, STI, STR, S2I,
     * S2R, TOL, TZI, TZR, YI, YR, ZEROI, ZEROR, ZI, ZR, D1MACH, ZABS
      INTEGER I, IB, IL, INU, J, JL, K, KODE, KODED, M, N, NN, NZ
      DIMENSION YR(N), YI(N)
      EXTERNAL ZABS, ZEXP, ZSQRT
      DATA PI, RTPI  /3.14159265358979324D0 , 0.159154943091895336D0 /
      DATA ZEROR,ZEROI,CONER,CONEI / 0.0D0, 0.0D0, 1.0D0, 0.0D0 /
C***FIRST EXECUTABLE STATEMENT  ZASYI
      NZ = 0
      AZ = ZABS(ZR,ZI)
      ARM = 1.0D+3*D1MACH(1)
      RTR1 = SQRT(ARM)
      IL = MIN(2,N)
      DFNU = FNU + (N-IL)
C-----------------------------------------------------------------------
C     OVERFLOW TEST
C-----------------------------------------------------------------------
      RAZ = 1.0D0/AZ
      STR = ZR*RAZ
      STI = -ZI*RAZ
      AK1R = RTPI*STR*RAZ
      AK1I = RTPI*STI*RAZ
      CALL ZSQRT(AK1R, AK1I, AK1R, AK1I)
      CZR = ZR
      CZI = ZI
      IF (KODE.NE.2) GO TO 10
      CZR = ZEROR
      CZI = ZI
   10 CONTINUE
      IF (ABS(CZR).GT.ELIM) GO TO 100
      DNU2 = DFNU + DFNU
      KODED = 1
      IF ((ABS(CZR).GT.ALIM) .AND. (N.GT.2)) GO TO 20
      KODED = 0
      CALL ZEXP(CZR, CZI, STR, STI)
      CALL ZMLT(AK1R, AK1I, STR, STI, AK1R, AK1I)
   20 CONTINUE
      FDN = 0.0D0
      IF (DNU2.GT.RTR1) FDN = DNU2*DNU2
      EZR = ZR*8.0D0
      EZI = ZI*8.0D0
C-----------------------------------------------------------------------
C     WHEN Z IS IMAGINARY, THE ERROR TEST MUST BE MADE RELATIVE TO THE
C     FIRST RECIPROCAL POWER SINCE THIS IS THE LEADING TERM OF THE
C     EXPANSION FOR THE IMAGINARY PART.
C-----------------------------------------------------------------------
      AEZ = 8.0D0*AZ
      S = TOL/AEZ
      JL = RL+RL + 2
      P1R = ZEROR
      P1I = ZEROI
      IF (ZI.EQ.0.0D0) GO TO 30
C-----------------------------------------------------------------------
C     CALCULATE EXP(PI*(0.5+FNU+N-IL)*I) TO MINIMIZE LOSSES OF
C     SIGNIFICANCE WHEN FNU OR N IS LARGE
C-----------------------------------------------------------------------
      INU = FNU
      ARG = (FNU-INU)*PI
      INU = INU + N - IL
      AK = -SIN(ARG)
      BK = COS(ARG)
      IF (ZI.LT.0.0D0) BK = -BK
      P1R = AK
      P1I = BK
      IF (MOD(INU,2).EQ.0) GO TO 30
      P1R = -P1R
      P1I = -P1I
   30 CONTINUE
      DO 70 K=1,IL
        SQK = FDN - 1.0D0
        ATOL = S*ABS(SQK)
        SGN = 1.0D0
        CS1R = CONER
        CS1I = CONEI
        CS2R = CONER
        CS2I = CONEI
        CKR = CONER
        CKI = CONEI
        AK = 0.0D0
        AA = 1.0D0
        BB = AEZ
        DKR = EZR
        DKI = EZI
        DO 40 J=1,JL
          CALL ZDIV(CKR, CKI, DKR, DKI, STR, STI)
          CKR = STR*SQK
          CKI = STI*SQK
          CS2R = CS2R + CKR
          CS2I = CS2I + CKI
          SGN = -SGN
          CS1R = CS1R + CKR*SGN
          CS1I = CS1I + CKI*SGN
          DKR = DKR + EZR
          DKI = DKI + EZI
          AA = AA*ABS(SQK)/BB
          BB = BB + AEZ
          AK = AK + 8.0D0
          SQK = SQK - AK
          IF (AA.LE.ATOL) GO TO 50
   40   CONTINUE
        GO TO 110
   50   CONTINUE
        S2R = CS1R
        S2I = CS1I
        IF (ZR+ZR.GE.ELIM) GO TO 60
        TZR = ZR + ZR
        TZI = ZI + ZI
        CALL ZEXP(-TZR, -TZI, STR, STI)
        CALL ZMLT(STR, STI, P1R, P1I, STR, STI)
        CALL ZMLT(STR, STI, CS2R, CS2I, STR, STI)
        S2R = S2R + STR
        S2I = S2I + STI
   60   CONTINUE
        FDN = FDN + 8.0D0*DFNU + 4.0D0
        P1R = -P1R
        P1I = -P1I
        M = N - IL + K
        YR(M) = S2R*AK1R - S2I*AK1I
        YI(M) = S2R*AK1I + S2I*AK1R
   70 CONTINUE
      IF (N.LE.2) RETURN
      NN = N
      K = NN - 2
      AK = K
      STR = ZR*RAZ
      STI = -ZI*RAZ
      RZR = (STR+STR)*RAZ
      RZI = (STI+STI)*RAZ
      IB = 3
      DO 80 I=IB,NN
        YR(K) = (AK+FNU)*(RZR*YR(K+1)-RZI*YI(K+1)) + YR(K+2)
        YI(K) = (AK+FNU)*(RZR*YI(K+1)+RZI*YR(K+1)) + YI(K+2)
        AK = AK - 1.0D0
        K = K - 1
   80 CONTINUE
      IF (KODED.EQ.0) RETURN
      CALL ZEXP(CZR, CZI, CKR, CKI)
      DO 90 I=1,NN
        STR = YR(I)*CKR - YI(I)*CKI
        YI(I) = YR(I)*CKI + YI(I)*CKR
        YR(I) = STR
   90 CONTINUE
      RETURN
  100 CONTINUE
      NZ = -1
      RETURN
  110 CONTINUE
      NZ=-2
      RETURN
      END


text:
c     =====================================================
      subroutine qinit(meqn,mbc,mx,my,xlower,ylower,
     &     dx,dy,q,maux,aux)
c     =====================================================
c
c     # Set initial conditions for q.
c     # Sample scalar equation with data that is piecewise constant with
c     # q = 1.0  if  0.1 < x < 0.6   and   0.1 < y < 0.6
c     #     0.1  otherwise
c
      implicit double precision (a-h,o-z)
      dimension q(meqn, 1-mbc:mx+mbc, 1-mbc:my+mbc)
      dimension aux(maux, 1-mbc:mx+mbc, 1-mbc:my+mbc)

c     # set concentration profile
c     ---------------------------

      qin = 1.d0    ! value inside blob
      qout = 0.d0   ! value outside

       do 50 i=1,mx
          xlow = xlower + (i-1)*dx
          do 20 j=1,my
c            # set (xlow,ylow) to lower left corner of grid cell:
             ylow = ylower + (j-1)*dy
             call cellave(xlow,ylow,dx,dy,win)
c            # win is now the fraction of the cell that lies inside the circle
             do 10 m=1,meqn
                q(m,i,j) = win*qin + (1.d0-win)*qout
  10            continue
  20         continue
  50      continue

      return
      end


text:
!===============================================================================
!  S H E B A - Shear-wave Birefringence Analysis
!===============================================================================
!  This software is distributed in the hope that it will be useful,
!  but WITHOUT ANY WARRANTY; without even the implied warranty of
!  MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.
!===============================================================================
!
!  James Wookey, School of Earth Sciences, University of Bristol
!
!===============================================================================
   program sheba_main
!===============================================================================
      use f90sac  ! use the f90sac module
      use sheba_config ! use the config module
!-------------------------------------------------------------------------------
      implicit none
         real, parameter :: sheba_version = 1.01

! ** print copyright information
     call sheba_banner()

!  ** initialise configuration
      call sheba_config_init()

!  ** get input file information
      call read_config()

!  ** run the analysis
      call sheba()

!  ** finish up

      stop
   end program sheba_main
!===============================================================================

max_stars_repo_name: jwookey/sheba

max_issues_repo_name: jwookey/sheba

max_forks_repo_name: jwookey/sheba


! Tests for correct parsing of complicated initialization
      MODULE funcon
      IMPLICIT NONE
      REAL(8),DIMENSION(3,3),PARAMETER :: &
         imat= reshape((/1.d0,0.d0,0.d0, &
                         0.d0,1.d0,0.d0, &
                         0.d0,0.d0,1.d0/),(/3,3/))
      REAL(8),DIMENSION(6,6),PARAMETER :: &
         imat6= reshape((/1.d0,0.d0,0.d0,0.d0,0.d0,0.d0, &
                          0.d0,1.d0,0.d0,0.d0,0.d0,0.d0, &
		 	  0.d0,0.d0,1.d0,0.d0,0.d0,0.d0, &
                          0.d0,0.d0,0.d0,1.d0,0.d0,0.d0, &
			  0.d0,0.d0,0.d0,0.d0,1.d0,0.d0, &
                          0.d0,0.d0,0.d0,0.d0,0.d0,1.d0/),(/6,6/))

      END MODULE funcon


text:
module initdata_module

end module initdata_module

max_stars_repo_name: pgrete/MAESTROeX

max_issues_repo_name: pgrete/MAESTROeX

max_forks_repo_name: pgrete/MAESTROeX


subroutine read_aerosol(nread,ndata,nodata,jsatid,infile,gstime,lunout, &
           obstype,twind,sis,ithin,rmesh, &
           mype_root,mype_sub,npe_sub,mpi_comm_sub,nobs)
!$$$  subprogram documentation block
!                .      .    .                                       .
! subprogram:    read_aerosol                    read aerosol data
!   prgmmr: hchuang     org: np23                date: 2009-01-26
!
! abstract:  This routine reads MODIS aerosol total column AOD observations.
!            ONLY total column values are read in.  The routine has
!            the ability to read both IEEE and BUFR format MODIS
!            as well as BUFR format VIIRS
!            aerosol data files.
!
!            When running the gsi in regional mode, the code only
!            retains those observations that fall within the regional
!            domain
!
! program history log:
!   2009-04-08  Huang   - modified from read_ozone to read in MODIS AEROSOL data
!   2010-10-20  hclin   - modified for total aod in channels
!   2011-01-05  hclin   - added three more BUFR records (STYP DBCF QAOD)
!   2011-08-01  lueken  - changed F90 to f90 (no machine logic)
!   2013-01-26  parrish - change from grdcrd to grdcrd1 (to allow successful debug compile on WCOSS)
!   2015-02-23  Rancic/Thomas - add thin4d to time window logical
!   2015-10-01  guo      - calc ob location once in deg
!   2018-05-21  j.jin    - added time-thinning. Moved the checking of thin4d into satthin.F90.
!   2019-03-22  martin  - add VIIRS BUFR capability based on code from S-W Wei and  Q. Zhao
!
!   input argument list:
!     obstype  - observation type to process
!     jsatid   - satellite id to read
!     infile   - unit from which to read aerosol data
!     gstime   - analysis time in minutes from reference date
!     lunout   - unit to which to write data for further processing
!     obstype  - observation type to process
!     twind    - input group time window (hours)
!     sis      - satellite/instrument/sensor indicator
!     ithin    - flag to thin data
!     rmesh    - thinning mesh size (km)
!     mype     - mpi task id
!     mype_root - "root" task for sub-communicator
!     mype_sub - mpi task id within sub-communicator
!     npe_sub  - number of data read tasks
!     mpi_comm_sub - sub-communicator for data read
!
!   output argument list:
!     nread    - number of modis aerosol observations read
!     ndata    - number of modis aerosol profiles retained for further processing
!     nodata   - number of modis aerosol observations retained for further processing
!     nobs     - array of observations on each subdomain for each processor
!
! remarks:
!
! attributes:
!   language: f90
!   machine:  IBM AIX Cirrus
!
!$$$
  use kinds,     only: r_kind, r_double, i_kind
  use gridmod,   only: nlat, nlon, regional, tll2xy, rlats, rlons
  use chemmod,   only: aod_qa_limit, luse_deepblue
  use constants, only: deg2rad, zero, one, two, three, four, five, r0_01, r60inv
  use obsmod,    only: rmiss_single
  use gsi_4dvar, only: l4dvar,l4densvar,iwinbgn,winlen
  use satthin,   only: itxmax,makegrids,destroygrids,checkob, &
      finalcheck,map2tgrid,score_crit
  use satthin,   only: radthin_time_info,tdiff2crit
  use obsmod,    only: time_window_max
  use mpimod, only: npe
  implicit none
!
! Declare local parameters
  real(r_kind), parameter :: r360 = 360.0_r_kind
!
! Declare passed variables
!
  character(len=*),intent(in)    :: obstype, infile, jsatid
  character(len=20),intent(in)   :: sis
  integer(i_kind), intent(in)    :: lunout, ithin
  integer(i_kind), intent(inout) :: nread
  integer(i_kind),dimension(npe), intent(inout) :: nobs
  integer(i_kind), intent(inout) :: ndata, nodata
  integer(i_kind) ,intent(in)    :: mype_root
  integer(i_kind) ,intent(in)    :: mype_sub
  integer(i_kind) ,intent(in)    :: npe_sub
  integer(i_kind) ,intent(in)    :: mpi_comm_sub
  real(r_kind),    intent(in)    :: gstime, twind, rmesh
!
! Declare local variables
!
  logical :: outside, iuse

  character (len= 8) :: subset
  character (len=10) :: date

  integer(i_kind) :: naerodat
  integer(i_kind) :: idate, jdate, ksatid, iy, iret, im, ihh, idd
  integer(i_kind) :: lunin = 10
  integer(i_kind) :: nmind, i, n
  integer(i_kind) :: k, ilat, ilon, nreal, nchanl
  integer(i_kind) :: kidsat
  integer(i_kind), dimension(5) :: idate5
!
!| NC008041 | SAID    AEROSOL  CLONH   CLATH YYMMDD  HHMMSS  SOZA  SOLAZI       |
!| NC008041 | SCATTA  OPTD  AEROTP                                              |
!
!| YYMMDD   | YEAR    MNTH    DAYS                                              |
!|          |                                                                   |
!| HHMMSS   | HOUR    MINU    SECO                                              |
!
!    SAID    Satellite identifier code table (eg, 783 == 'TERRA')
!    AEROSOL Aerosol Optical Depth (AOD) source code table (eg, 5 == 'AATSR' )
!    YEAR    Year
!    MNTH    Month
!    DAYS    Day
!    HOUR    Hour
!    MINU    Minute
!    SECO    Second
!    CLATH   Latitude (high accuracy)     degree (5 decimal precision)
!    CLONH   Longitude (high accuracy)    degree (5 decimal precision)
!    SOLAZI  Solar azimuth                degree (2 decimal precision)
!    SOZA    Solar zenith angle           degree (2 decimal precision)
!    OPTD    Optical depth                numeric
!    SCATTA  Scattering angle             degree (2 decimal precsion)
!    AEROTP  Aerosol type land            code table (eg, 1 == 'DUST', 2 == 'SULFATE')
!
!    0-15-195 - AEROTP (Aerosol land type)
!
!    CODE  DESCRIPTION
!    ====  ===========
!    0     Mixed
!    1     Dust
!    2     Sulfate
!    3     Smoke
!    4     Heavy absorbing smoke
!    5-14  Reserved
!    15    Missing value
!
  character (len= 4) :: aerostr  = 'OPTD'
  character (len=53) :: aerogstr = &
      'SAID CLATH CLONH YEAR MNTH DAYS HOUR MINU SOZA SOLAZI'

! VIIRS AOD code
  character (len= 9) :: vaodchstr  = 'CHWL AOTH'
  character (len=69) :: vaodgstr = &
      'SAID CLATH CLONH YEAR MNTH DAYS HOUR MINU SOZA SOLAZI RSST VAOTQ QPLR'
  integer(i_kind), parameter :: mxib  = 20
  integer(i_kind) :: nib
  integer(i_kind) :: ibit(mxib)

  integer(i_kind) :: itx, itt, irec

  real(r_kind) :: tdiff, sstime, dlon, dlat, t4dv, crit1, dist1
  real(r_kind) :: slons0, slats0, rsat, solzen, azimuth, dlat_earth, dlon_earth
  real(r_kind) :: dlat_earth_deg, dlon_earth_deg
  real(r_kind) :: styp, dbcf, qaod, smask, qcall
  real(r_kind) :: qcall_limit  ! qcall >= qcall_limit will be retained

  real(r_kind),dimension(0:6):: rlndsea

  real(r_kind), allocatable, dimension(:,:) :: aeroout
  real(r_kind), allocatable, dimension(:)   :: dataaod
  integer(i_kind),allocatable,dimension(:)  :: nrec
  real(r_double), dimension( 10) :: hdraerog
  real(r_double)                 :: aod_550
  real(r_double), dimension(3)   :: aod_flags
  real(r_kind)    :: ptime,timeinflat,crit0
  integer(i_kind) :: ithin_time,n_tbin,it_mesh
! for VIIRS
  real(r_double),dimension(13)   :: hdrvaodg
  real(r_double),dimension(2,12) :: vaodch
  real(r_double)                 :: aod_lb,aod_ub

!**************************************************************************
! Set constants.  Initialize variables
  rsat=999._r_kind
  ! output position of LON and LAT
  ilon=3
  ilat=4
  nread = 0
  ndata = 0
  nodata = 0

  ! Set rlndsea for types we would prefer selecting
  rlndsea(0) = zero        ! styp 0: water
  rlndsea(1) = 15._r_kind  ! styp 1: coast
  rlndsea(2) = 20._r_kind  ! styp 2: desert
  rlndsea(3) = 10._r_kind  ! styp 3: land
  rlndsea(4) = 25._r_kind  ! styp 4: deep blue
  rlndsea(5) = 30._r_kind  ! styp 5: nnr ocean
  rlndsea(6) = 35._r_kind  ! styp 6: nnr land

  call radthin_time_info(obstype, jsatid, sis, ptime, ithin_time)
  if( ptime > 0.0_r_kind) then
     n_tbin=nint(2*time_window_max/ptime)
  else
     n_tbin=1
  endif
! Make thinning grids
  call makegrids(rmesh,ithin,n_tbin=n_tbin)

  if ( obstype == 'modis_aod' ) then
!
     open(lunin,file=trim(infile),form='unformatted')
     call openbf(lunin,'IN',lunin)
     call datelen(10)
     call readmg(lunin,subset,idate,iret)

     if ( iret == 0 ) then
!
        if (subset == 'NC008041') then
           write(6,*)'READ_AEROSOL: MODIS data type, subset = ',subset
           !          Set dependent variables and allocate arrays
           nreal=11   !9
           nchanl=20 ! 19 + 1 additional vis channel in CRTM coeff file
           naerodat=nreal+nchanl
           allocate (aeroout(naerodat,itxmax),nrec(itxmax))
           allocate (dataaod(nchanl))

           iy = 0
           im = 0
           idd= 0
           ihh= 0
           write(date,'( i10)') idate
           read (date,'(i4,3i2)') iy,im,idd,ihh
           write(6,'(''READ_AEROSOL: aerosol bufr file '',a,''  date is '',i4,4i2.2,a)')trim(infile),iy,im,idd,ihh

           nrec=999999
           irec=0
           read_modis: do
              irec=irec+1
              call readsb(lunin,iret)
              if (iret/=0) then
                 call readmg(lunin,subset,jdate,iret)
                 if (iret/=0) exit read_modis
                 cycle read_modis
              endif

              !    extract header information
              call ufbint(lunin,hdraerog,10,1,iret,aerogstr)
              rsat = hdraerog(1); ksatid=rsat

              if ( jsatid == 'terra' ) kidsat = 783
              if ( jsatid == 'aqua'  ) kidsat = 784

              if ( ksatid /= kidsat  ) cycle read_modis

              !    Convert observation location to radians
              slats0= hdraerog(2)
              slons0= hdraerog(3)
              if(slons0< zero) slons0=slons0+r360
              if(slons0>=r360) slons0=slons0-r360
              dlat_earth_deg = slats0
              dlon_earth_deg = slons0
              dlat_earth = slats0 * deg2rad
              dlon_earth = slons0 * deg2rad

              if(regional)then
                 call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
                 if(outside) cycle read_modis
              else
                 dlat = dlat_earth
                 dlon = dlon_earth
                 call grdcrd1(dlat,rlats,nlat,1)
                 call grdcrd1(dlon,rlons,nlon,1)
              endif

              solzen  = hdraerog(9)
              azimuth = hdraerog(10)

              !    Convert observation time to relative time
              idate5(1) = hdraerog(4)  !year
              idate5(2) = hdraerog(5)  !month
              idate5(3) = hdraerog(6)  !day
              idate5(4) = hdraerog(7)  !hour
              idate5(5) = hdraerog(8)  !minute

              !    extract total column aod 1 value 'OPTD' as defined in aerostr
              call ufbint(lunin,aod_550,1,1,iret,aerostr)

              call w3fs21(idate5,nmind)
              t4dv=real((nmind-iwinbgn),r_kind)*r60inv
              sstime=real(nmind,r_kind)
              tdiff=(sstime-gstime)*r60inv
              if (l4dvar.or.l4densvar) then
                 if(t4dv<zero .OR. t4dv>winlen) cycle read_modis
              else
                 if ( abs(tdiff) > twind ) cycle read_modis
              end if

              nread = nread + 1   !nread = nread + nchanl

              if ( aod_550 > 1.0e+10_r_double ) cycle read_modis

              ! extract STYP, DBCF, and QAOD
              ! these are missing from the 008041 bufr files
              styp = rmiss_single
              dbcf = rmiss_single
              qaod = zero

              if ( .not. luse_deepblue .and. nint(styp)==4 ) cycle read_modis
              if ( qaod > aod_qa_limit ) cycle read_modis

              ! Map obs to thinning grid
              crit0 = 0.01_r_kind
              timeinflat=two
              call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh)
              call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis,it_mesh=it_mesh)
              if ( .not. iuse ) cycle read_modis

              if ( (styp > rmiss_single) .and. (styp >= zero .and. styp <= four) ) then
                 crit1 = crit1 + rlndsea(nint(styp))
              end if
              if ( (qaod > rmiss_single) .and. (qaod >= aod_qa_limit .and. qaod <= three) ) then
                 crit1 = crit1 + 10.0_r_kind*(four-qaod)
              end if
              call checkob(dist1,crit1,itx,iuse)
              if ( .not. iuse ) cycle read_modis

              ! Compute "score" for observation.  All scores>=0.0.  Lowest score is "best"
              call finalcheck(dist1,crit1,itx,iuse)
              if ( .not. iuse ) cycle read_modis

              dataaod = rmiss_single
              dataaod(4) = aod_550

              aeroout( 1,itx) = rsat
              aeroout( 2,itx) = tdiff
              aeroout( 3,itx) = dlon               ! grid relative longitude
              aeroout( 4,itx) = dlat               ! grid relative latitude
              aeroout( 5,itx) = dlon_earth_deg     ! earth relative longitude (degrees)
              aeroout( 6,itx) = dlat_earth_deg     ! earth relative latitude (degrees)
              aeroout( 7,itx) = qaod               ! total column AOD error flag
              aeroout( 8,itx) = solzen             ! solar zenith angle
              aeroout( 9,itx) = azimuth            ! solar azimuth angle
              aeroout(10,itx) = styp               ! surface type
              aeroout(11,itx) = dbcf               ! deep blue confidence flag
              do i = 1, nchanl
                 aeroout(i+nreal,itx) = dataaod(i)
              end do
              nrec(itx)=irec

           end do read_modis

           call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
              naerodat,itxmax,nread,ndata,aeroout,score_crit,nrec)

           if ( mype_sub == mype_root ) then
              do n = 1, ndata
                 do i = 1, nchanl
                    if ( aeroout(i+nreal,n) > rmiss_single ) nodata = nodata + 1
                 end do
              end do
              ! Write final set of "best" observations to output file
              call count_obs(ndata,naerodat,ilat,ilon,aeroout,nobs)
              write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
              write(lunout) ((aeroout(k,n),k=1,naerodat),n=1,ndata)
           end if

           ! Deallocate local arrays
           deallocate(aeroout)
           deallocate(dataaod)

           ! End of MODIS bufr block
        else       ! subset /= NC008041
           write(6,*)'READ_AEROSOL:  *** WARNING: unknown aerosol data type, subset=',subset
           write(6,*)' infile=',infile, ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid
           write(6,*)' SKIP PROCESSING OF THIS MODIS FILE'
        endif

     else          ! read subset iret /= 0
        write(6,*)'READ_AEROSOL:  *** WARNING: read subset error, obstype=',obstype,', iret=',iret
     end if
     call closbf(lunin)
     close(lunin)
  else if ( obstype == 'viirs_aod' ) then
     open(lunin,file=trim(infile),form='unformatted')
     call openbf(lunin,'IN',lunin)
     call datelen(10)
     call readmg(lunin,subset,idate,iret)

     if ( iret == 0 ) then

        if (subset == 'NC008043') then
           write(6,*)'READ_AEROSOL: VIIRS AOD data type, subset = ',subset
           !          Set dependent variables and allocate arrays
           nreal=10
           nchanl=11
           naerodat=nreal+nchanl
           allocate (aeroout(naerodat,itxmax),nrec(itxmax))
           allocate (dataaod(nchanl))

           iy = 0
           im = 0
           idd= 0
           ihh= 0
           write(date,'( i10)') idate
           read (date,'(i4,3i2)') iy,im,idd,ihh
           write(6,'(''READ_AEROSOL: aerosol bufr file '',a,''  date is '',i4,3i2.2)') trim(infile),iy,im,idd,ihh

!          set qcall_limit
           if (idate >= 2018021300) then
              qcall_limit = aod_qa_limit + r0_01 ! for the viirs data after 2018/02/13
           else
              qcall_limit = aod_qa_limit - r0_01
           end if

!          set valid range of AOD to ingest
           aod_lb = zero
           aod_ub = five
           nrec=999999
           irec=0
           read_viirs: do
              irec=irec+1
              call readsb(lunin,iret)
              if (iret/=0) then
                 call readmg(lunin,subset,jdate,iret)
                 if (iret/=0) exit read_viirs
                 cycle read_viirs
              endif

              !    extract header information
              call ufbint(lunin,hdrvaodg,13,1,iret,vaodgstr)
              rsat = hdrvaodg(1); ksatid=rsat

              if ( jsatid == 'NPP' .or. jsatid == 'npp' ) kidsat = 225

              if ( ksatid /= kidsat  ) cycle read_viirs

              !    Convert observation location to radians
              slats0= hdrvaodg(2)
              slons0= hdrvaodg(3)
              if(slons0< zero) slons0=slons0+r360
              if(slons0>=r360) slons0=slons0-r360
              dlat_earth_deg = slats0
              dlon_earth_deg = slons0
              dlat_earth = slats0 * deg2rad
              dlon_earth = slons0 * deg2rad

              if(regional)then
                 call tll2xy(dlon_earth,dlat_earth,dlon,dlat,outside)
                 if(outside) cycle read_viirs
              else
                 dlat = dlat_earth
                 dlon = dlon_earth
                 call grdcrd1(dlat,rlats,nlat,1)
                 call grdcrd1(dlon,rlons,nlon,1)
              endif

              solzen  = hdrvaodg(9)
              azimuth = hdrvaodg(10)

              smask = zero                          ! over water
              if (nint(hdrvaodg(11)) > 0) then      ! over land
                 smask = one                        ! dark surface
                 call upftbv(lunin,"VAOTQ",hdrvaodg(12),mxib,ibit,nib)
                 if (nib > 0) then
                    if(any(ibit(1:nib) == 6)) then
                       smask = two                  !  bright surface
                    endif
                 endif
              endif

              qcall   = hdrvaodg(13)

              !    Convert observation time to relative time
              idate5(1) = hdrvaodg(4)  !year
              idate5(2) = hdrvaodg(5)  !month
              idate5(3) = hdrvaodg(6)  !day
              idate5(4) = hdrvaodg(7)  !hour
              idate5(5) = hdrvaodg(8)  !minute

              call w3fs21(idate5,nmind)
              t4dv=real((nmind-iwinbgn),r_kind)*r60inv
              sstime=real(nmind,r_kind)
              tdiff=(sstime-gstime)*r60inv

              if (l4dvar.or.l4densvar) then
                 if(t4dv<zero .OR. t4dv>winlen) cycle read_viirs
              else
                 if ( abs(tdiff) > twind ) cycle read_viirs
              end if

              nread = nread + 1   !nread = nread + nchanl

              if (idate >= 2018021300) then
                 if ( qcall > qcall_limit ) cycle read_viirs
              else
                 if ( qcall < qcall_limit ) cycle read_viirs
              end if

              !    extract VAODCH pairs 'CHWL AOPT' as defined in vaodchstr
              call ufbrep(lunin,vaodch,2,12,iret,vaodchstr)
              aod_550 = vaodch(2,12)

              if ( aod_550 < aod_lb .OR. aod_550 >= aod_ub ) cycle read_viirs

              ! Map obs to thinning grid
              crit0 = 0.01_r_kind
              timeinflat=two
              call tdiff2crit(tdiff,ptime,ithin_time,timeinflat,crit0,crit1,it_mesh)
              call map2tgrid(dlat_earth,dlon_earth,dist1,crit1,itx,ithin,itt,iuse,sis)

              if ( .not. iuse ) cycle read_viirs

              crit1 = crit1 + 10.0_r_kind*smask
              ! is below needed now because of the change in QC flags? CRM
              if (idate >= 2018021300) then
                crit1 = crit1 + 10.0_r_kind*(four+qcall)
              else
                crit1 = crit1 + 10.0_r_kind*(four-qcall)
              end if
              call checkob(dist1,crit1,itx,iuse)

              if ( .not. iuse ) cycle read_viirs

              ! Compute "score" for observation.  All scores>=0.0.  Lowest score
              ! is "best"
              call finalcheck(dist1,crit1,itx,iuse)

              if ( .not. iuse ) cycle read_viirs

              dataaod = rmiss_single
              dataaod(4) = aod_550

              aeroout( 1,itx) = rsat
              aeroout( 2,itx) = tdiff
              aeroout( 3,itx) = dlon               ! grid relative longitude
              aeroout( 4,itx) = dlat               ! grid relative latitude
              aeroout( 5,itx) = dlon_earth_deg     ! earth relative longitude (degrees)
              aeroout( 6,itx) = dlat_earth_deg     ! earth relative latitude (degrees)
              aeroout( 7,itx) = qcall              ! total column AOD error flag
              aeroout( 8,itx) = solzen             ! solar zenith angle
              aeroout( 9,itx) = azimuth            ! solar azimuth angle
              aeroout(10,itx) = smask              ! surface type mask
              do i = 1, nchanl
                 aeroout(i+nreal,itx) = dataaod(i)
              enddo
              nrec(itx)=irec

           end do read_viirs

           call combine_radobs(mype_sub,mype_root,npe_sub,mpi_comm_sub,&
              naerodat,itxmax,nread,ndata,aeroout,score_crit,nrec)

           if ( mype_sub == mype_root ) then
              do n = 1, ndata
                 do i = 1, nchanl
                    if ( aeroout(i+nreal,n) > rmiss_single ) nodata = nodata + 1
                 end do
              end do
              ! Write final set of "best" observations to output file
              call count_obs(ndata,naerodat,ilat,ilon,aeroout,nobs)
              write(lunout) obstype,sis,nreal,nchanl,ilat,ilon
              write(lunout) ((aeroout(k,n),k=1,naerodat),n=1,ndata)
           end if

           ! Deallocate local arrays
           deallocate(aeroout)

           ! End of VIIRS AOD bufr block

        else       ! subset /= NC008043
           write(6,*)'READ_AEROSOL:  *** WARNING: unknown aerosol data type, subset=',subset
           write(6,*)' infile=',infile, ', lunin=',lunin, ', obstype=',obstype,', jsatid=',jsatid
           write(6,*)' SKIP PROCESSING OF THIS VIIRS_AOD FILE'
        endif

     else          ! read subset iret /= 0
        write(6,*)'READ_AEROSOL:  *** WARNING: read subset error, obstype=',obstype,', iret=',iret
     end if
     call closbf(lunin)
     close(lunin)

  else             ! obstype /= 'modis' or 'viirs'
     write(6,*)'READ_AEROSOL:  *** WARNING: unknown aerosol input type, obstype=',obstype
  endif

  ! Deallocate satthin arrays
  call destroygrids

end subroutine read_aerosol


text:
! computes using Simpson's method
program CompSimp
implicit none
double precision x0,xn,x,h,Simp,exact,diff,f
integer j,n,m
write(*,*) 'enter n (must be even)'
read(*,*) n
m=n/2
x0=0.d0
xn=1.d0
h=(xn-x0)/dfloat(n)
Simp=f(x0)+f(xn)
do j=1,m-1
  x=x0+dfloat(2*j)*h
  Simp=Simp+2.d0*f(x)
enddo
do j=1,m
  x=x0+dfloat(2*j-1)*h
  Simp=Simp+4.d0*f(x)
enddo
Simp=h/3.d0 *Simp
exact=0.528408
diff=dabs(Simp-exact)
write(*,*)'exact=',exact
write(*,*)'Simp=',Simp
write(*,*)'|Simp-exact|=',diff
stop
end program CompSimp

!!!! 0.63629436111989061883446424291635313615100026

function f(x)
double precision x,f
f=x**(-1/4) *dsin(x)
return
end function f

max_stars_repo_name: whisker-rebellion/Fortran

max_issues_repo_name: whisker-rebellion/Fortran

max_forks_repo_name: whisker-rebellion/Fortran


text:
program ex5

    implicit none

    integer, parameter :: rk = selected_real_kind(15,300)

    real(rk) :: A(20,20)
    real(rk) :: B(20,20)
    real(rk) :: C(20,20)
    integer  :: i, j

    A(10,:) = 1.0_rk
    A(:,10) = 2.0_rk

    B = 0.0_rk

    do i=1,20
        do j=1,20
            C(i,j) = 42.0_rk
        end do
    end do

end program ex5

max_stars_repo_name: jonaslindemann/compute-course-public

max_issues_repo_name: jonaslindemann/compute-course-public

max_forks_repo_name: jonaslindemann/compute-course-public


! Daniel R. Reynolds
! SMU Mathematics
! Math 4370/6370
! 7 February 2015
!-----------------------------------------------------------------
! Description:
!    Evolves the first-order 2D wave equations in time.
!=================================================================

program advection
  !===============================================================
  ! inclusions
  implicit none

  ! interfaces
  interface
     subroutine output(u,t,nx,ny,noutput)
       integer, intent(in) :: nx, ny, noutput
       real*8,  intent(in) :: t, u(nx,ny)
     end subroutine output
     subroutine initialize(u,v1,v2,v3,c,dx,dy,nx,ny)
       integer, intent(in)  :: nx, ny
       real*8,  intent(in)  :: c, dx, dy
       real*8, dimension(nx,ny), intent(out) :: u, v1, v2, v3
     end subroutine initialize
  end interface

  ! declarations
  integer :: nx, ny, nt, i, j, noutput, it
  real*8 :: tstop, c, dtoutput, dx, dy, t, toutput, dt
  real*8 :: stime, ftime, v2_E, v2_W, v3_N, v3_S, v1_W, v1_E, v1_S, v1_N
  real*8 :: runtime, iotime, inittime
  real*8, dimension(:,:), allocatable :: u, v1, v2, v3
  namelist /inputs/ nx, ny, nt, tstop, c, dtoutput

  ! internals

  call get_time(stime)
  ! read problem parameters from input file (should be in this order):
  !    nx - number of nodes in x-direction
  !    ny - number of nodes in y-direction
  !    nt - number of time steps
  !    tstop - final time (will stop at nt or stop, whichever is 1st)
  !    c - wave speed
  !    dtoutput - time frequency for outputting solutions
  open(100,file='input.txt')
  read(100,inputs)
  close(100)
  print *, '  '
  print *, 'Running wave problem:'
  print *, '   nx =',nx,',  ny =',ny
  print *, '   nt =',nt,',  tstop =',tstop
  print *, '   c =',c
  print *, '   dtoutput =',dtoutput
  print *, '  '
  call get_time(ftime)
  iotime = ftime-stime

  ! allocate arrays
  call get_time(stime)
  allocate(u(nx,ny),v1(nx,ny),v2(nx,ny),v3(nx,ny))

  ! set grid spacing
  dx = 1.d0/nx
  dy = 1.d0/ny

  ! set initial conditions
  call initialize(u,v1,v2,v3,c,dx,dy,nx,ny)
  call get_time(ftime)
  inittime = ftime-stime

  ! set initial time, output initial solution
  t = 0.d0
  toutput = 0.d0
  noutput = 0
  print '(A,i4,A,i8,A,es9.2)', 'writing output file ',noutput, &
       ', step = ',0,',  t =',t
  call get_time(stime)
  call output(u,t,nx,ny,noutput)
  call get_time(ftime)
  iotime = iotime+ftime-stime

  ! set time step
  dt = min(dx/c/5.d1, dy/c/5.d1)

  ! start time stepping
  runtime = 0.d0
  do it=1,nt

     ! start timer
     call get_time(stime)

     ! first update v1 to get to half time step
     !    get relevant values for this location
     do j=1,ny
        do i=1,nx

           ! access relevant components of v2 and v3
           if (i==nx) then
              v2_E = v2(1,j)
           else
              v2_E = v2(i+1,j)
           end if
           v2_W = v2(i,j)
           if (j==ny) then
              v3_N = v3(i,1)
           else
              v3_N = v3(i,j+1)
           end if
           v3_S = v3(i,j)

           ! update v1
           v1(i,j) = v1(i,j) + c*dt/dx*(v2_E - v2_W) + c*dt/dy*(v3_N - v3_S)

        end do
     end do

     ! next update v2 & v3 to get to full time step
     !    get relevant values for this location
     do j=1,ny
        do i=1,nx

           ! access relevant components of v1
           if (i==1) then
              v1_W = v1(nx,j)
           else
              v1_W = v1(i-1,j)
           end if
           v1_E = v1(i,j)
           if (j==1) then
              v1_S = v1(i,ny)
           else
              v1_S = v1(i,j-1)
           end if
           v1_N = v1(i,j)

           ! update v2 and v3
           v2(i,j) = v2(i,j) + c*dt/dx*(v1_E - v1_W)
           v3(i,j) = v3(i,j) + c*dt/dy*(v1_N - v1_S)

        end do
     end do

     ! update and plot solution
     u = u + dt*v1

     ! update time
     t = t + dt
     call get_time(ftime)
     runtime = runtime+ftime-stime

     ! stop simulation if we've reached tstop
     if (t >= tstop)  exit

     ! output solution periodically
     if (abs(t-toutput-dtoutput) <= 1.e-14) then
        call get_time(stime)
        toutput = t
        noutput = noutput+1
        print '(A,i4,A,i8,A,es9.2)', 'writing output file ',noutput, &
             ', step = ',it,',  t =',t
        call output(u,t,nx,ny,noutput)
        call get_time(ftime)
        iotime = iotime+ftime-stime
     end if

  end do

  ! output final solution
  call get_time(stime)
  toutput = t
  noutput = noutput+1
  print '(A,i4,A,i8,A,es9.2)', 'writing output file ',noutput, &
       ', step = ',it,',  t =',t
  call output(u,t,nx,ny,noutput)
  call get_time(ftime)
  iotime = iotime+ftime-stime

  ! output runtime
  print *, ' total initialization time = ', inittime
  print *, ' total input/output time   = ', iotime
  print *, ' total simulation time     = ', runtime

  ! free memory
  deallocate(u,v1,v2,v3)

  ! end program
end program advection
!=================================================================