!----------------------------------------------------------------------------------------------- ! chemistry driver; uses either the Rosenbrock or Mozart solver method !----------------------------------------------------------------------------------------------- module chemistry_driver use kinetics_module, only : kinetics_type use ccpp_kinds, only : kind_phys use ODE_solver, only : baseOdeSolver use Rosenbrock_Solver, only : RosenbrockSolver use Mozart_Solver, only : MozartSolver implicit none private public :: chemistry_driver_init public :: chemistry_driver_run integer :: icntrl(20) ! integer control array for ODE solver real(kind_phys) :: absTol=-huge(1._kind_phys), relTol=-huge(1._kind_phys) real(kind_phys) :: rcntrl(20) ! real control array for ODE solver character(len=80) :: Solver_method = 'NONE' contains !> \section arg_table_chemistry_driver_init Argument Table !! \htmlinclude chemistry_driver_init.html !! subroutine chemistry_driver_init(nSpecies, nkRxt, njRxt, TimeStart, TimeEnd, dt, options_filepath, print_log_message, errmsg, errflg) implicit none !----------------------------------------------------------- ! these dimension parameters will be set by the cafe/configurator !----------------------------------------------------------- real(kind_phys), intent(in) :: TimeStart, TimeEnd real(kind_phys), intent(in) :: dt integer, intent(in) :: nSpecies ! number prognostic constituents integer, intent(in) :: nkRxt ! number gas phase reactions integer, intent(in) :: njRxt ! number of photochemical reactions character(len=*), intent(in) :: options_filepath logical, intent(in) :: print_log_message character(len=512),intent(out) :: errmsg integer, intent(out) :: errflg ! error index from CPF ! real(kind_phys) :: Hstart real(kind_phys) :: Abs_tol(nSpecies), Rel_tol(nSpecies) namelist /micm_solv_opts/ Solver_method namelist /micm_solv_err_cntrl/ absTol, relTol ! namelist /timeCntrl/ Hstart !----------------------------------------------------------- ! get and set the Solver method !----------------------------------------------------------- open(unit=10,file=options_filepath) read(unit=10,nml=micm_solv_opts) read(unit=10,nml=micm_solv_err_cntrl) ! read(unit=10,nml=timeCntrl) close(unit=10) !--- initialize CCPP error handling variables errmsg = '' errflg = 0 !----------------------------------------------------------- ! initialize method control arrays !----------------------------------------------------------- icntrl(:) = 0 rcntrl(:) = 0._kind_phys select case(theSolver) case( 'rosenbrock', 'Rosenbrock', 'ROSENBROCK' ) !----------------------------------------------------------- ! set ode solver "control" variables for Rosenbrock solver !----------------------------------------------------------- icntrl(1) = 1 ! autonomous, F depends only on Y icntrl(3) = 2 ! ros3 solver rcntrl(2) = dt ! Hmax rcntrl(3) = .01_kind_phys*dt ! Hstart case( 'mozart', 'Mozart', 'MOZART' ) !----------------------------------------------------------- ! set ode solver "control" variables for MOZART solver !----------------------------------------------------------- icntrl(1) = 1 ! autonomous, F depends only on Y rcntrl(2) = dt ! Hmax, max time step rcntrl(3) = .01_kind_phys*dt ! Hstart, initial time step case default write(errmsg,*) 'chemistry_driver: solver method must be {Mozart,Rosenbrock}; not ',trim(Solver_method) errflg = -1 return end select if (print_log_message) then Abs_tol(:) = absTol Rel_tol(:) = relTol write(*,*) ' ' write(*,*) 'icntrl settings' write(*,'(10i6)') icntrl(1:10) write(*,*) 'rcntrl settings' write(*,'(1p,10(1x,g10.4))') rcntrl(1:10) write(*,*) 'Absolute error tolerances' write(*,*) Abs_tol(:) write(*,*) 'Relative error tolerances' write(*,*) Rel_tol(:) write(*,*) ' ' end if end subroutine chemistry_driver_init !> \section arg_table_chemistry_driver_run Argument Table !! \htmlinclude chemistry_driver_run.html !! subroutine chemistry_driver_run(vmr, TimeStart, TimeEnd, j_rateConst, k_rateConst, number_density_air, errmsg, errflg) use kinetics_utilities, only: kinetics_init, kinetics_final implicit none !----------------------------------------------------------- ! these dimension parameters will be set by the cafe/configurator !----------------------------------------------------------- real(kind=kind_phys), intent(inout) :: vmr(:) ! "working" concentration passed thru CPF real(kind_phys), intent(in) :: TimeStart, TimeEnd real(kind=kind_phys), intent(in) :: j_rateConst(:) ! host model provides photolysis rates for now real(kind=kind_phys), intent(in) :: k_rateConst(:) ! host model provides photolysis rates for now real(kind=kind_phys), intent(in) :: number_density_air ! number density character(len=512), intent(out) :: errmsg integer, intent(out) :: errflg ! error index from CPF real(kind=kind_phys) :: number_density(size(vmr)) ! "working" number density of each molecule integer :: nTotRxt ! total number of chemical reactions integer :: nSpecies ! number prognostic constituents real(kind_phys) :: Abs_tol(size(vmr)), Rel_tol(size(vmr)) class(baseOdeSolver), pointer :: theSolver => null() type(RosenbrockSolver), target :: aRosenbrockSolver type(MozartSolver), target :: aMozartSolver type(kinetics_type) :: theKinetics !--- initialize CCPP error handling variables errmsg = '' errflg = 0 select case( Solver_method ) case( 'mozart', 'Mozart', 'MOZART' ) theSolver => aMozartSolver case( 'rosenbrock', 'Rosenbrock', 'ROSENBROCK' ) theSolver => aRosenbrockSolver case default write(errmsg,*) 'chemistry_driver: solver method must be {Mozart,Rosenbrock}; not ',trim(Solver_method) errflg = -1 return end select !----------------------------------------------------------- ! initialize the solver !----------------------------------------------------------- Abs_tol(:) = absTol Rel_tol(:) = relTol theSolver%print_log_message = .false. call theSolver%Initialize( Tstart=TimeStart, Tend=TimeEnd, AbsTol=Abs_tol, RelTol=Rel_tol, & ICNTRL=icntrl, RCNTRL=rcntrl, Ierr=errflg ) !----------------------------------------------------------- ! initialize the kinetics type !----------------------------------------------------------- nTotRxt = size(j_rateConst) + size(k_rateConst) nSpecies = size(vmr) call thekinetics%init( nTotRxt, nSpecies ) call kinetics_init(vmr, number_density, number_density_air) !----------------------------------------------------------- ! update the kinetics !----------------------------------------------------------- call theKinetics%rateConst_update( k_rateConst, j_rateConst, number_density_air ) !----------------------------------------------------------- ! solve the current timestep's chemistry !----------------------------------------------------------- call theSolver%Run( Tstart=TimeStart, Tend=TimeEnd, y=number_density, theKinetics=theKinetics, Ierr=errflg ) call kinetics_final(vmr, number_density, number_density_air) if (errflg /= 0) then errmsg = 'chemistry_driver_run: ERROR in theSolver%Run' end if end subroutine chemistry_driver_run !> \section arg_table_chemistry_driver_finalize Argument Table !! \htmlinclude chemistry_driver_finalize.html !! subroutine chemistry_driver_finalize(errmsg,errflg) character(len=512), intent(out) :: errmsg integer, intent(out) :: errflg ! error index from CPF errmsg = '' errflg = 0 end subroutine chemistry_driver_finalize end module chemistry_driver