part_3_code/box_model.F90

Go to the documentation of this file.

program box_model
 
  use camp_camp_core
  use camp_camp_state
  use camp_chem_spec_data
  use camp_constants
 
  !! [NO2 photolysis modules]
  use camp_mechanism_data
  use camp_rxn_data
  use camp_rxn_photolysis
  use camp_rxn_factory
 
  !! [NO2 photolysis modules]
 
#ifdef CAMP_USE_MPI
  use mpi
  use camp_mpi
#endif
 
  implicit none
 
  type(camp_core_t), pointer :: camp_core
  type(camp_state_t), pointer :: camp_state
 
  integer(kind=i_kind) :: idx_o3, idx_no, idx_no2, idx_o2
  type(chem_spec_data_t), pointer :: chem_spec_data
 
  character(len=*), parameter :: fmt_hdr = "(A10,',',A10,',',A10,',',A10,',',A10)"
  character(len=*), parameter :: fmt_dat = "(ES10.4,',',ES10.4,',',ES10.4,',',ES10.4,',',ES10.4)"
 
  integer(kind=i_kind) :: i_time
 
  !! [NO2 photolysis variables]
  integer(kind=i_kind) :: i_rxn
  character(len=:), allocatable :: photo_label
  type(mechanism_data_t), pointer :: mechanism
  type(rxn_factory_t) :: rxn_factory
  class(rxn_data_t), pointer :: photo_rxn
  type(rxn_update_data_photolysis_t) :: no2_photolysis
 
  !! [NO2 photolysis variables]
 
  ! allows example to run with MPI supprt compiled in
#ifdef CAMP_USE_MPI
  call camp_mpi_init( )
#endif
 
  camp_core => camp_core_t( "my_config_file.json" )
  call camp_core%initialize( )
 
  if( .not.camp_core%get_chem_spec_data( chem_spec_data ) ) then
    write(*,*) "Something's gone wrong!"
    stop 3
  end if
 
  idx_o3  = chem_spec_data%gas_state_id( "O3"  )
  idx_no  = chem_spec_data%gas_state_id( "NO"  )
  idx_no2 = chem_spec_data%gas_state_id( "NO2" )
  idx_o2  = chem_spec_data%gas_state_id( "O2"  )
  if( idx_o3.eq.0 .or. idx_no2.eq.0 .or.idx_o2.eq.0 ) then
    write(*,*) "Missing species!"
    stop 3
  end if
 
  !! [Find NO2 photolysis]
  if( .not.camp_core%get_mechanism( "my simple mechanism", mechanism ) ) then
    write(*,*) "Missing mechanism!"
    stop 3
  end if
 
  do i_rxn = 1, mechanism%size( )
    photo_rxn => mechanism%get_rxn( i_rxn )
    select type( photo_rxn )
      class is( rxn_photolysis_t )
        if( photo_rxn%property_set%get_string( "my photo label", photo_label ) ) then
          if( photo_label .eq. "NO2 photolysis" ) then
            call camp_core%initialize_update_object( photo_rxn, no2_photolysis )
          end if
        end if
    end select
  end do
  !! [Find NO2 photolysis]
 
  call camp_core%solver_initialize( )
  camp_state => camp_core%new_state( )
 
  call camp_state%env_states(1)%set_temperature_K(  275.4_dp )
  call camp_state%env_states(1)%set_pressure_Pa( 101532.2_dp )
 
  camp_state%state_var( idx_o3   ) = 0.13  ! [O3] in ppm
  camp_state%state_var( idx_no   ) = 0.02  ! [NO] in ppm
  camp_state%state_var( idx_no2  ) = 0.053 ! [NO2] in ppm
  camp_state%state_var( idx_o2   ) = 2.1e5 ! [O2] in ppm
 
  !! [Set NO2 photolysis]
  call no2_photolysis%set_rate( 12.2d0 ) ! rate in s-1
  call camp_core%update_data( no2_photolysis )
  !! [Set NO2 photolysis]
 
  write(*,fmt_hdr) "time", "O3", "NO", "NO2", "O2"
  do i_time = 1, 100
    call camp_core%solve( camp_state, 1.0d-15 ) ! time step in s
    write(*,fmt_dat) i_time*1.0e-15, &
                     camp_state%state_var( idx_o3  ), &
                     camp_state%state_var( idx_no  ), &
                     camp_state%state_var( idx_no2 ), &
                     camp_state%state_var( idx_o2  )
  end do
 
  deallocate( camp_core )
  deallocate( camp_state )
 
#ifdef CAMP_USE_MPI
  call camp_mpi_finalize( )
#endif
 
end program box_model