LCOV - code coverage report
Current view: | top level - chemistry/mozart - rate_diags.F90 (source / functions) |  |  | Hit | Total | Coverage
Test: | coverage.info |  | Lines: | 12 | 66 | 18.2 %
Date: | 2024-12-17 17:57:11 |  | Functions: | 1 | 6 | 16.7 %

              Line data    Source code
           1             : !--------------------------------------------------------------------------------
           2             : ! Manages writing reaction rates to history
           3             : !--------------------------------------------------------------------------------
           4             : module rate_diags
           5             :
           6             :   use shr_kind_mod,     only : r8 => shr_kind_r8
           7             :   use shr_kind_mod,     only : CL => SHR_KIND_CL
           8             :   use cam_history,      only : fieldname_len
           9             :   use cam_history,      only : addfld, add_default
          10             :   use cam_history,      only : outfld
          11             :   use chem_mods,        only : rxt_tag_cnt, rxt_tag_lst, rxt_tag_map
          12             :   use ppgrid,           only : pver
          13             :   use spmd_utils,       only : masterproc
          14             :   use cam_logfile,      only : iulog
          15             :   use cam_abortutils,   only : endrun
          16             :   use shr_expr_parser_mod , only : shr_exp_parse, shr_exp_item_t, shr_exp_list_destroy
          17             :
          18             :   implicit none
          19             :   private
          20             :   public :: rate_diags_init
          21             :   public :: rate_diags_calc
          22             :   public :: rate_diags_readnl
          23             :   public :: rate_diags_o3s_loss
          24             :   public :: rate_diags_final
          25             :
          26             :   character(len=fieldname_len) :: rate_names(rxt_tag_cnt)
          27             :
          28             :   type(shr_exp_item_t), pointer :: grps_list => null()
          29             :
          30             :   integer, parameter :: maxlines = 200
          31             :   character(len=CL), allocatable :: rxn_rate_sums(:)
          32             :
          33             :   integer :: o3_ndx = -1
          34             :
          35             : contains
          36             :
          37             : !-------------------------------------------------------------------
          38             : !-------------------------------------------------------------------
          39        1536 :   subroutine rate_diags_readnl(nlfile)
          40             :
          41             :     use namelist_utils,  only: find_group_name
          42             :     use units,           only: getunit, freeunit
          43             :     use spmd_utils,      only: mpicom, mpi_character, masterprocid
          44             :
          45             :     ! args
          46             :     character(len=*), intent(in) :: nlfile  ! filepath for file containing namelist input
          47             :
          48             :     ! Local variables
          49             :     integer :: unitn, ierr
          50             :
          51             :     namelist /rxn_rate_diags_nl/ rxn_rate_sums
          52             :
          53        1536 :     allocate( rxn_rate_sums( maxlines ) )
          54      308736 :     rxn_rate_sums(:) = ' '
          55             :
          56             :     ! Read namelist
          57        1536 :     if (masterproc) then
          58           2 :        unitn = getunit()
          59           2 :        open( unitn, file=trim(nlfile), status='old' )
          60           2 :        call find_group_name(unitn, 'rxn_rate_diags_nl', status=ierr)
          61           2 :        if (ierr == 0) then
          62           0 :           read(unitn, rxn_rate_diags_nl, iostat=ierr)
          63           0 :           if (ierr /= 0) then
          64           0 :              call endrun('rate_diags_readnl:: ERROR reading namelist')
          65             :           end if
          66             :        end if
          67           2 :        close(unitn)
          68           2 :        call freeunit(unitn)
          69             :     end if
          70             :
          71             :     ! Broadcast namelist variables
          72        1536 :     call mpi_bcast(rxn_rate_sums,len(rxn_rate_sums(1))*maxlines, mpi_character, masterprocid, mpicom, ierr)
          73             :
          74        1536 :   end subroutine rate_diags_readnl
          75             : !--------------------------------------------------------------------------------
          76             : !--------------------------------------------------------------------------------
          77           0 :   subroutine rate_diags_init
          78             :     use phys_control, only : phys_getopts
          79             :     use mo_chem_utls, only : get_spc_ndx
          80             :
          81             :     integer :: i,j, len, pos
          82             :     character(len=64) :: name
          83             :     logical :: history_scwaccm_forcing
          84             :     type(shr_exp_item_t), pointer :: grp
          85             :
          86           0 :     call phys_getopts( history_scwaccm_forcing_out = history_scwaccm_forcing )
          87             :
          88             :     do i = 1,rxt_tag_cnt
          89             :        pos = 0
          90             :        pos = index(rxt_tag_lst(i),'tag_')
          91             :        if (pos <= 0) pos = index(rxt_tag_lst(i),'usr_')
          92             :        if (pos <= 0) pos = index(rxt_tag_lst(i),'cph_')
          93             :        if (pos <= 0) pos = index(rxt_tag_lst(i),'ion_')
          94             :        if (pos>0) then
          95             :           name = 'r_'//trim(rxt_tag_lst(i)(5:))
          96             :        else
          97             :           name = 'r_'//trim(rxt_tag_lst(i)(1:))
          98             :        endif
          99             :        len = min(fieldname_len,len_trim(name))
         100             :        rate_names(i) = trim(name(1:len))
         101             :        call addfld(rate_names(i), (/ 'lev' /),'A', 'molecules/cm3/sec','reaction rate')
         102             :        if (history_scwaccm_forcing .and. rate_names(i) == 'r_O1D_H2O') then
         103             :           call add_default( rate_names(i), 1, ' ')
         104             :        endif
         105             :     enddo
         106             :
         107             :     ! parse the terms of the summations
         108           0 :     grps_list => shr_exp_parse( rxn_rate_sums )
         109           0 :     deallocate( rxn_rate_sums )
         110             :
         111           0 :     if (masterproc) write(iulog,*) 'rate_diags_init :'
         112             :
         113           0 :     grp => grps_list
         114           0 :     do while(associated(grp))
         115             :
         116           0 :        if (masterproc) then
         117           0 :           write(iulog,*) ' grp name : ',trim(grp%name)
         118             :
         119           0 :           do j = 1, grp%n_terms
         120           0 :              write(iulog,'(f12.4,a,a)') grp%coeffs(j),' * ',trim(grp%vars(j))
         121             :           end do
         122             :        end if
         123             :
         124           0 :        call addfld( grp%name, (/ 'lev' /),'A', 'molecules/cm3/sec','reaction rate group')
         125             :
         126           0 :        grp => grp%next_item
         127             :
         128             :     enddo
         129             :
         130           0 :     o3_ndx = get_spc_ndx('O3')
         131             :
         132           0 :   end subroutine rate_diags_init
         133             :
         134             : !--------------------------------------------------------------------------------
         135             : !--------------------------------------------------------------------------------
         136           0 :   subroutine rate_diags_calc( rxt_rates, vmr, m, ncol, lchnk )
         137             :
         138             :     use mo_rxt_rates_conv, only: set_rates
         139             :
         140             :     real(r8), intent(inout) :: rxt_rates(:,:,:) ! 'molec/cm3/sec'
         141             :     real(r8), intent(in)    :: vmr(:,:,:)
         142             :     real(r8), intent(in)    :: m(:,:)           ! air density (molecules/cm3)
         143             :     integer,  intent(in)    :: ncol, lchnk
         144             :
         145             :     integer :: i, j, ndx
         146           0 :     real(r8) :: group_rate(ncol,pver)
         147             :     type(shr_exp_item_t), pointer :: grp
         148             :
         149           0 :     call set_rates( rxt_rates, vmr, ncol )
         150             :
         151             :     ! output individual tagged rates
         152             :     do i = 1, rxt_tag_cnt
         153             :        ! convert from vmr/sec to molecules/cm3/sec
         154             :        rxt_rates(:ncol,:,rxt_tag_map(i)) = rxt_rates(:ncol,:,rxt_tag_map(i)) * m(:ncol,:)
         155             :        call outfld( rate_names(i), rxt_rates(:ncol,:,rxt_tag_map(i)), ncol, lchnk )
         156             :     enddo
         157             :
         158             :     ! output rate groups ( or families )
         159             :
         160           0 :     grp => grps_list
         161           0 :     do while(associated(grp))
         162             :
         163           0 :        group_rate(:,:) = 0._r8
         164           0 :        do j = 1, grp%n_terms
         165           0 :          ndx = lookup_tag_ndx(grp%vars(j))
         166           0 :          group_rate(:ncol,:) = group_rate(:ncol,:) + grp%coeffs(j)*rxt_rates(:ncol,:,ndx)
         167             :        enddo
         168           0 :        call outfld( grp%name, group_rate(:ncol,:), ncol, lchnk )
         169             :
         170           0 :        grp => grp%next_item
         171             :
         172             :     end do
         173             :
         174           0 :   end subroutine rate_diags_calc
         175             :
         176             : !--------------------------------------------------------------------------------
         177             : !--------------------------------------------------------------------------------
         178           0 :   function rate_diags_o3s_loss( rxt_rates, vmr, ncol ) result(o3s_loss)
         179             :     use mo_rxt_rates_conv, only: set_rates
         180             :     use chem_mods,         only: rxntot
         181             :
         182             :     real(r8), intent(in) :: rxt_rates(:,:,:)
         183             :     real(r8), intent(in) :: vmr(:,:,:)
         184             :     integer,  intent(in) :: ncol
         185             :
         186             :     real(r8) :: o3s_loss(ncol,pver) ! /sec
         187             :
         188             :     integer :: j, ndx
         189           0 :     real(r8) :: group_rate(ncol,pver)
         190           0 :     real(r8) :: lcl_rxt_rates(ncol,pver,rxntot)
         191             :     type(shr_exp_item_t), pointer :: grp
         192             :
         193           0 :     o3s_loss(:,:) = 0._r8
         194             :
         195           0 :     if (o3_ndx>0) then
         196           0 :        lcl_rxt_rates(:ncol,:,:) = rxt_rates(:ncol,:,:)
         197           0 :        call set_rates( lcl_rxt_rates, vmr, ncol )
         198             :
         199           0 :        grp => grps_list
         200           0 :        loop: do while(associated(grp))
         201             :
         202           0 :           if (trim(grp%name)=='O3S_Loss') then
         203           0 :              group_rate(:,:) = 0._r8
         204           0 :              do j = 1, grp%n_terms
         205           0 :                 ndx = lookup_tag_ndx(grp%vars(j))
         206           0 :                 group_rate(:ncol,:) = group_rate(:ncol,:) + grp%coeffs(j)*lcl_rxt_rates(:ncol,:,ndx)
         207             :              enddo
         208           0 :              o3s_loss(:ncol,:) = group_rate(:ncol,:)/vmr(:ncol,:,o3_ndx)
         209           0 :              exit loop
         210             :           endif
         211           0 :           grp => grp%next_item
         212             :
         213             :        end do loop
         214             :     endif
         215             :
         216           0 :   end function rate_diags_o3s_loss
         217             :
         218             : !--------------------------------------------------------------------------------
         219             : !--------------------------------------------------------------------------------
         220           0 :   subroutine rate_diags_final
         221             :
         222           0 :     if (associated(grps_list)) call shr_exp_list_destroy(grps_list)
         223             :
         224           0 :   end subroutine rate_diags_final
         225             :
         226             : !-------------------------------------------------------------------
         227             : ! Private routines :
         228             : !-------------------------------------------------------------------
         229             : !-------------------------------------------------------------------
         230             :
         231             : !-------------------------------------------------------------------
         232             : ! finds the index corresponging to a given reacton name
         233             : !-------------------------------------------------------------------
         234           0 :   function lookup_tag_ndx( name ) result( ndx )
         235             :     character(len=*) :: name
         236             :     integer :: ndx
         237             :
         238             :     integer :: i
         239             :
         240           0 :     ndx = -1
         241             :
         242             :     findloop: do i = 1,rxt_tag_cnt
         243             :        if (trim(name) .eq. trim(rate_names(i)(3:))) then
         244             :           ndx = i
         245             :           return
         246             :        endif
         247             :     end do findloop
         248             :
         249             :     if (ndx<0) then
         250           0 :        call endrun('rate_diags: not able to find rxn tag name: '//trim(name))
         251             :     endif
         252             :
         253             :   end function lookup_tag_ndx
         254             :
         255             : end module rate_diags

Generated by: LCOV version 1.14