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        1536 :   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        1536 :     call phys_getopts( history_scwaccm_forcing_out = history_scwaccm_forcing )
      87             : 
      88       24576 :     do i = 1,rxt_tag_cnt
      89       23040 :        pos = 0
      90       23040 :        pos = index(rxt_tag_lst(i),'tag_')
      91       23040 :        if (pos <= 0) pos = index(rxt_tag_lst(i),'usr_')
      92       23040 :        if (pos <= 0) pos = index(rxt_tag_lst(i),'cph_')
      93       19968 :        if (pos <= 0) pos = index(rxt_tag_lst(i),'ion_')
      94       19968 :        if (pos>0) then
      95        3072 :           name = 'r_'//trim(rxt_tag_lst(i)(5:))
      96             :        else
      97       19968 :           name = 'r_'//trim(rxt_tag_lst(i)(1:))
      98             :        endif
      99       23040 :        len = min(fieldname_len,len_trim(name))
     100       23040 :        rate_names(i) = trim(name(1:len))
     101       46080 :        call addfld(rate_names(i), (/ 'lev' /),'A', 'molecules/cm3/sec','reaction rate')
     102       24576 :        if (history_scwaccm_forcing .and. rate_names(i) == 'r_O1D_H2O') then
     103           0 :           call add_default( rate_names(i), 1, ' ')
     104             :        endif
     105             :     enddo
     106             : 
     107             :     ! parse the terms of the summations
     108        1536 :     grps_list => shr_exp_parse( rxn_rate_sums )
     109        1536 :     deallocate( rxn_rate_sums )
     110             : 
     111        1536 :     if (masterproc) write(iulog,*) 'rate_diags_init :'
     112             : 
     113        1536 :     grp => grps_list
     114        1536 :     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        1536 :     o3_ndx = get_spc_ndx('O3')
     131             : 
     132        1536 :   end subroutine rate_diags_init
     133             : 
     134             : !--------------------------------------------------------------------------------
     135             : !--------------------------------------------------------------------------------
     136     1489176 :   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     2978352 :     real(r8) :: group_rate(ncol,pver)
     147             :     type(shr_exp_item_t), pointer :: grp
     148             : 
     149     1489176 :     call set_rates( rxt_rates, vmr, ncol )
     150             : 
     151             :     ! output individual tagged rates
     152    23826816 :     do i = 1, rxt_tag_cnt
     153             :        ! convert from vmr/sec to molecules/cm3/sec
     154 34710095160 :        rxt_rates(:ncol,:,rxt_tag_map(i)) = rxt_rates(:ncol,:,rxt_tag_map(i)) * m(:ncol,:)
     155    23826816 :        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     1489176 :     grp => grps_list
     161     1489176 :     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     1489176 :   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        1536 :   subroutine rate_diags_final
     221             : 
     222        1536 :     if (associated(grps_list)) call shr_exp_list_destroy(grps_list)
     223             : 
     224        1536 :   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           0 :     findloop: do i = 1,rxt_tag_cnt
     243           0 :        if (trim(name) .eq. trim(rate_names(i)(3:))) then
     244           0 :           ndx = i
     245           0 :           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           0 :   end function lookup_tag_ndx
     254             : 
     255             : end module rate_diags