Line data    Source code

       1             : module charge_neutrality
       2             : 
       3             :   use shr_kind_mod, only : r8 => shr_kind_r8
       4             :   use ppgrid,       only : pcols, pver
       5             :   use mo_chem_utls, only : get_spc_ndx
       6             : 
       7             :   implicit none
       8             : 
       9             :   private
      10             :   public :: charge_balance
      11             : 
      12             :   interface charge_balance
      13             :      module procedure charge_fix_vmr
      14             :      module procedure charge_fix_mmr   ! for fixing charge balance after vertical diffusion
      15             :   end interface
      16             : 
      17             :   integer, parameter :: pos_ion_n = 22
      18             :   character(len=16), parameter :: pos_ion_names(pos_ion_n) = (/ &
      19             :        'Np              ','N2p             ','Op              ','O2p             ','NOp             ', &
      20             :        'O4p             ','O2p_H2O         ','Hp_H2O          ','Hp_2H2O         ','Hp_3H2O         ', &
      21             :        'Hp_4H2O         ','Hp_5H2O         ','H3Op_OH         ','Hp_3N1          ','Hp_4N1          ', &
      22             :        'NOp_H2O         ','NOp_2H2O        ','NOp_3H2O        ','NOp_CO2         ','NOp_N2          ', &
      23             :        'Op2P            ','Op2D            ' /)
      24             : 
      25             :   integer, parameter :: neg_ion_n = 21
      26             :   character(len=16), parameter :: neg_ion_names(neg_ion_n) = (/ &
      27             :        'Om              ','O2m             ','O3m             ','O4m             ','OHm             ', &
      28             :        'CO3m            ','CO4m            ','NO2m            ','NO3m            ','HCO3m           ', &
      29             :        'CLm             ','CLOm            ','CLm_H2O         ','CLm_HCL         ','CO3m_H2O        ', &
      30             :        'NO3m_H2O        ','CO3m2H2O        ','NO2m_H2O        ','NO3m2H2O        ','NO3mHNO3        ', &
      31             :        'NO3m_HCL        ' /)
      32             : 
      33             : contains
      34             : 
      35             :   !-----------------------------------------------------------------------      
      36             :   !        ... force ion/electron balance
      37             :   !-----------------------------------------------------------------------      
      38     2978352 :   subroutine charge_fix_vmr( ncol, vmr )
      39             : 
      40             :     !-----------------------------------------------------------------------      
      41             :     !        ... dummy arguments
      42             :     !-----------------------------------------------------------------------      
      43             :     integer,  intent(in)    :: ncol
      44             :     real(r8), intent(inout) :: vmr(:,:,:)         ! concentration
      45             : 
      46             :     !-----------------------------------------------------------------------      
      47             :     !        ... local variables
      48             :     !-----------------------------------------------------------------------      
      49             :     integer  :: i, n
      50             :     integer  :: elec_ndx
      51     5956704 :     real(r8) :: wrk(ncol,pver)
      52             : 
      53     2978352 :     elec_ndx = get_spc_ndx('e')
      54             : 
      55             :     !--------------------------------------------------------------------
      56             :     ! If electrons are in the chemistry add up charges to get electrons
      57             :     !--------------------------------------------------------------------  
      58     2978352 :     if( elec_ndx > 0 ) then
      59           0 :        wrk(:,:) = 0._r8
      60             : 
      61           0 :        do i = 1,pos_ion_n
      62           0 :           n = get_spc_ndx(pos_ion_names(i))
      63           0 :           if (n>0) then
      64           0 :              wrk(:ncol,:) = wrk(:ncol,:) + vmr(:ncol,:,n)
      65             :           endif
      66             :        enddo
      67           0 :        do i = 1,neg_ion_n
      68           0 :           n = get_spc_ndx(neg_ion_names(i))
      69           0 :           if (n>0) then 
      70           0 :              wrk(:ncol,:) = wrk(:ncol,:) - vmr(:ncol,:,n)
      71             :           endif
      72             :        enddo
      73             : 
      74           0 :        where ( wrk(:,:)<0._r8 )
      75             :           wrk(:,:)=0._r8
      76             :        end where
      77             : 
      78           0 :        vmr(:ncol,:,elec_ndx) = wrk(:ncol,:)
      79             :       
      80             :     end if
      81             : 
      82     2978352 :   end subroutine charge_fix_vmr
      83             : 
      84             :   !-----------------------------------------------------------------------
      85             :   !        ... force ion/electron balance
      86             :   !-----------------------------------------------------------------------
      87     1489176 :   subroutine charge_fix_mmr(state, pbuf)
      88             : 
      89             :     use constituents,        only : cnst_get_ind
      90             :     use air_composition,     only : mbarv                       ! Constituent dependent mbar
      91             :     use short_lived_species, only : slvd_index,slvd_pbf_ndx => pbf_idx ! Routines to access short lived species in pbuf
      92             :     use chem_mods,           only : adv_mass
      93             :     use physics_buffer,      only : pbuf_get_field,physics_buffer_desc ! Needed to get variables from physics buffer
      94             :     use physics_types,       only : physics_state
      95             : 
      96             :     !-----------------------------------------------------------------------      
      97             :     !        ... dummy arguments
      98             :     !-----------------------------------------------------------------------      
      99             :     type(physics_state), intent(inout), target :: state
     100             :     type(physics_buffer_desc), pointer :: pbuf(:)    ! physics buffer
     101             : 
     102             :     !-----------------------------------------------------------------------      
     103             :     !        ... local variables
     104             :     !-----------------------------------------------------------------------      
     105             :     integer  :: i, n, ns, nc
     106             :     integer  :: elec_ndx
     107             :     integer  :: lchnk                 !Chunk number from state structure
     108             :     integer  :: ncol                  !Number of columns in this chunk from state structure
     109             : 
     110     1489176 :     real(r8), dimension(:,:,:), pointer :: q         ! model mass mixing ratios
     111     1489176 :     real(r8), dimension(:,:),   pointer :: qs        ! Pointer to access fields in pbuf
     112             : 
     113             :     character(len=16) :: name
     114     2978352 :     real(r8) :: vmr(state%ncol,pver)  
     115     2978352 :     real(r8) :: wrk(state%ncol,pver)
     116             : 
     117             :     !-----------------------------------------------------------------------      
     118     1489176 :     elec_ndx = get_spc_ndx('e')
     119             : 
     120             :     !--------------------------------------------------------------------
     121             :     ! If electrons are simulated enforce charge neutrality ...
     122             :     !--------------------------------------------------------------------  
     123     1489176 :     if( elec_ndx > 0 ) then
     124           0 :        lchnk = state%lchnk
     125           0 :        ncol  = state%ncol
     126           0 :        q => state%q
     127           0 :        wrk(:,:) = 0._r8
     128             : 
     129           0 :        do i = 1,pos_ion_n+neg_ion_n
     130           0 :           if (i .le. pos_ion_n) then
     131           0 :              name = pos_ion_names(i)
     132             :           else
     133           0 :              name = neg_ion_names(i-pos_ion_n)
     134             :           endif
     135           0 :           n = get_spc_ndx(name) 
     136             : 
     137           0 :           if (n>0) then
     138           0 :              call cnst_get_ind( name, nc, abort=.false. )
     139           0 :              if (nc>0) then
     140           0 :                 vmr(:ncol,:) = mbarv(:ncol,:,lchnk) * q(:ncol,:,nc) / adv_mass(n)
     141             :              else
     142             :                 ! not transported
     143           0 :                 ns = slvd_index( name )
     144           0 :                 if (ns>0) then
     145           0 :                    call pbuf_get_field(pbuf, slvd_pbf_ndx, qs, start=(/1,1,ns/), kount=(/pcols,pver,1/) )
     146           0 :                    vmr(:ncol,:) = mbarv(:ncol,:,lchnk) * qs(:ncol,:) / adv_mass(n)
     147             :                 endif
     148             :              endif
     149           0 :              if (i .le. pos_ion_n) then
     150           0 :                 wrk(:ncol,:) = wrk(:ncol,:) + vmr(:ncol,:)
     151             :              else
     152           0 :                 wrk(:ncol,:) = wrk(:ncol,:) - vmr(:ncol,:)
     153             :              endif
     154             :           end if
     155             :        end do
     156             : 
     157           0 :        where ( wrk(:,:)<0._r8 )
     158             :           wrk(:,:)=0._r8
     159             :        end where
     160             : 
     161           0 :        call cnst_get_ind( 'e', nc, abort=.false. )  
     162             : 
     163           0 :        if (nc>0) then 
     164           0 :           q(:ncol,:,nc) = adv_mass(elec_ndx) * wrk(:ncol,:) / mbarv(:ncol,:,lchnk)
     165             :        else
     166             :           ! not transported
     167           0 :           ns = slvd_index( 'e' )
     168           0 :           call pbuf_get_field(pbuf, slvd_pbf_ndx, qs, start=(/1,1,ns/), kount=(/pcols,pver,1/) )
     169           0 :           qs(:ncol,:) = adv_mass(elec_ndx) * wrk(:ncol,:) / mbarv(:ncol,:,lchnk)
     170             :        endif
     171             : 
     172             :     endif
     173             : 
     174     1489176 :   end subroutine charge_fix_mmr
     175             : 
     176             : end module charge_neutrality