Rainout and transport rates are not accessible from the Python wrapper

[DanyIvan]

I would like to access rainout, transport rates, and distributed fluxes directly from the photochem Python wrapper. These rates are useful to construct pathways in photochem.

These variables could be included in the wrapper with these code additions: https://github.com/Nicholaswogan/photochem/compare/main...DanyIvan:photochem:chempath

I can create a pull request with these changes if you think this is ok.

[Nicholaswogan]

Daniel,

Thanks for these changes! They look great. Would it be possible to adjust the code for getting the transport_rates and distributed_fluxes so that it looks something like this just to avoid the same calculation being done multiple times?

Also, is the current transport_rates missing the surface and top-of-atmosphere transport rates?

    ! diffusion (interior grid points)
    do j = 2,var%nz-1
      do i = 1,dat%nq
        k = i + (j-1)*dat%nq
        advection_rates = ...
        diffusion_rates = ...
        transport rates = advection_rates + diffusion_rates

        wrk%transport_rates(i, j) = transport_rates
        wrk%advection_rates(i, j) =  advection_rates
        wrk%diffusion_rates(i, j) = transport rates

        rhs(k) = rhs(k) + transport_rates

      enddo
    enddo

[Nicholaswogan]

If it would work to open a PR, then make those minor adjustments, then it would be great to merge it so that there is better compatibility with chempath

[DanyIvan]

I will modify the code to avoid repeating calculations. Yes, transport_rates is missing the surface and top-of-atmosphere rates. I was unsure how to include these rates. Would something like this work?

! Lower boundary
do i = 1,dat%nq
      if (var%lowerboundcond(i) == VelocityBC .or. &
          var%lowerboundcond(i) == VelocityDistributedFluxBC) then
        transport_rates(i, 1) = wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) &
                        + wrk%DD(i,1)*wrk%usol(i,1) + wrk%ADD(i,1)*wrk%usol(i,1) &
                        - wrk%lower_vdep_copy(i)*wrk%usol(i,1)/var%dz(1)
        wrk%transport_rates(i, 1) = transport_rates(i, 1)
        rhs(i) = rhs(i) + transport_rates(i, 1) 
      elseif (var%lowerboundcond(i) == DensityBC .or. &
              var%lowerboundcond(i) == PressureBC) then
        rhs(i) = 0.0_dp
        wrk%transport_rates(i, 1) = 0.0_dp
      elseif (var%lowerboundcond(i) == FluxBC) then
        transport_rates(i, 1) = wrk%DU(i,1)*wrk%usol(i,2) + &
              wrk%ADU(i,1)*wrk%usol(i,2) + wrk%DD(i,1)*wrk%usol(i,1) + &
              wrk%ADD(i,1)*wrk%usol(i,1) + var%lower_flux(i)/var%dz(1)
        wrk%transport_rates(i, 1) = transport_rates(i, 1)
        rhs(i) = rhs(i) + transport_rates(i, 1) 
      ! Moses (2001) boundary condition for gas giants
      ! A deposition velocity controled by how quickly gases
      ! turbulantly mix vertically
      elseif (var%lowerboundcond(i) == MosesBC) then
        transport_rates(i, 1) = wrk%DU(i,1)*wrk%usol(i,2) + wrk%ADU(i,1)*wrk%usol(i,2) &
                        + wrk%DD(i,1)*wrk%usol(i,1) + wrk%ADD(i,1)*wrk%usol(i,1) &
                        - (var%edd(1)/wrk%scale_height(1))*wrk%usol(i,1)/var%dz(1)
        wrk%transport_rates(i, 1) = transport_rates(i, 1)
        rhs(i) = rhs(i) + transport_rates(i, 1) 
      endif
enddo

[Nicholaswogan]

This looks right, I think. Except the fixed boundary case. In this case the total sum of chemistry + tranpsort + other processes should cause a zero rate of change. So, the transport rates for this case might be something like the following:

photochem/src/evoatmosphere/photochem_evoatmosphere_utils.f90

Line 63 in b4d9110

! surface flux is molecules required to sustain the lower boundary

[DanyIvan]

Thanks,

I modified the code to avoid double calculations: main...DanyIvan:photochem:chempath

However, I am having problems testing these changes because I have not been able to install Photochem from source with python -m pip install --no-deps --no-build-isolation . I get an error related to CMake when I try this, probably related to my conda environment.

Do you have suggestions on how to test these changes?

[Nicholaswogan]

Ya, check out these instructions: https://github.com/Nicholaswogan/photochem/blob/docs/docs/docs/build_from_source.md

[DanyIvan]

Thanks. I made a couple of tests, and the changes seem to be working. I will go ahead with the pull request.

By the way, I had to use export CMAKE_POLICY_VERSION_MINIMUM=3.5 for the installation from source to work.