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 43776 : 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 87552 : real(r8) :: wrk(ncol,pver)
52 :
53 43776 : elec_ndx = get_spc_ndx('e')
54 :
55 : !--------------------------------------------------------------------
56 : ! If electrons are in the chemistry add up charges to get electrons
57 : !--------------------------------------------------------------------
58 43776 : if( elec_ndx > 0 ) then
59 74025216 : wrk(:,:) = 0._r8
60 :
61 1006848 : do i = 1,pos_ion_n
62 963072 : n = get_spc_ndx(pos_ion_names(i))
63 1006848 : if (n>0) then
64 1628554752 : wrk(:ncol,:) = wrk(:ncol,:) + vmr(:ncol,:,n)
65 : endif
66 : enddo
67 963072 : do i = 1,neg_ion_n
68 919296 : n = get_spc_ndx(neg_ion_names(i))
69 963072 : if (n>0) then
70 1554529536 : wrk(:ncol,:) = wrk(:ncol,:) - vmr(:ncol,:,n)
71 : endif
72 : enddo
73 :
74 74025216 : where ( wrk(:,:)<0._r8 )
75 : wrk(:,:)=0._r8
76 : end where
77 :
78 74025216 : vmr(:ncol,:,elec_ndx) = wrk(:ncol,:)
79 :
80 : end if
81 :
82 43776 : end subroutine charge_fix_vmr
83 :
84 : !-----------------------------------------------------------------------
85 : ! ... force ion/electron balance
86 : !-----------------------------------------------------------------------
87 43776 : 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 43776 : real(r8), dimension(:,:,:), pointer :: q ! model mass mixing ratios
111 43776 : real(r8), dimension(:,:), pointer :: qs ! Pointer to access fields in pbuf
112 :
113 : character(len=16) :: name
114 87552 : real(r8) :: vmr(state%ncol,pver)
115 87552 : real(r8) :: wrk(state%ncol,pver)
116 :
117 : !-----------------------------------------------------------------------
118 43776 : elec_ndx = get_spc_ndx('e')
119 :
120 : !--------------------------------------------------------------------
121 : ! If electrons are simulated enforce charge neutrality ...
122 : !--------------------------------------------------------------------
123 43776 : if( elec_ndx > 0 ) then
124 43776 : lchnk = state%lchnk
125 43776 : ncol = state%ncol
126 43776 : q => state%q
127 74025216 : wrk(:,:) = 0._r8
128 :
129 1926144 : do i = 1,pos_ion_n+neg_ion_n
130 1882368 : if (i .le. pos_ion_n) then
131 963072 : name = pos_ion_names(i)
132 : else
133 919296 : name = neg_ion_names(i-pos_ion_n)
134 : endif
135 1882368 : n = get_spc_ndx(name)
136 :
137 1926144 : if (n>0) then
138 1882368 : call cnst_get_ind( name, nc, abort=.false. )
139 1882368 : if (nc>0) then
140 0 : vmr(:ncol,:) = mbarv(:ncol,:,lchnk) * q(:ncol,:,nc) / adv_mass(n)
141 : else
142 : ! not transported
143 1882368 : ns = slvd_index( name )
144 1882368 : if (ns>0) then
145 7529472 : call pbuf_get_field(pbuf, slvd_pbf_ndx, qs, start=(/1,1,ns/), kount=(/pcols,pver,1/) )
146 3183084288 : vmr(:ncol,:) = mbarv(:ncol,:,lchnk) * qs(:ncol,:) / adv_mass(n)
147 : endif
148 : endif
149 1882368 : if (i .le. pos_ion_n) then
150 1628554752 : wrk(:ncol,:) = wrk(:ncol,:) + vmr(:ncol,:)
151 : else
152 1554529536 : wrk(:ncol,:) = wrk(:ncol,:) - vmr(:ncol,:)
153 : endif
154 : end if
155 : end do
156 :
157 74025216 : where ( wrk(:,:)<0._r8 )
158 : wrk(:,:)=0._r8
159 : end where
160 :
161 43776 : call cnst_get_ind( 'e', nc, abort=.false. )
162 :
163 43776 : if (nc>0) then
164 0 : q(:ncol,:,nc) = adv_mass(elec_ndx) * wrk(:ncol,:) / mbarv(:ncol,:,lchnk)
165 : else
166 : ! not transported
167 43776 : ns = slvd_index( 'e' )
168 175104 : call pbuf_get_field(pbuf, slvd_pbf_ndx, qs, start=(/1,1,ns/), kount=(/pcols,pver,1/) )
169 74025216 : qs(:ncol,:) = adv_mass(elec_ndx) * wrk(:ncol,:) / mbarv(:ncol,:,lchnk)
170 : endif
171 :
172 : endif
173 :
174 87552 : end subroutine charge_fix_mmr
175 :
176 : end module charge_neutrality
|
