Line data Source code
1 : module ssatcontrail
2 : ! contrail parameterization
3 : ! see Chen et al., 2012: Global contrail coverage simulated
4 : ! by CAM5 with the inventory of 2006 global aircraft emissions, JAMES
5 : ! https://doi.org/10.1029/2011MS000105
6 : use shr_kind_mod, only: r8 => shr_kind_r8
7 : use ppgrid, only: pcols, pver
8 : use physics_types, only: physics_state, physics_ptend, physics_ptend_init
9 : use physconst, only: cpair,mwdry,mwh2o, gravit, zvir, rair, pi, rearth, tmelt
10 : use physics_buffer, only: pbuf_get_index, pbuf_get_field, physics_buffer_desc, pbuf_old_tim_idx
11 : use constituents, only: cnst_get_ind, pcnst
12 : use phys_grid, only: get_wght_all_p
13 : use wv_saturation, only: qsat_water, qsat_ice
14 : use aircraft_emit, only: get_aircraft
15 :
16 : implicit none
17 : private
18 : save
19 :
20 : public ssatcontrail_d0
21 :
22 : ! Private data
23 : real(r8), parameter :: rhoi = 500.0_r8 ! density of ice (500 kg/m3)
24 : real(r8), parameter :: radius = 3.75e-6_r8 ! diameter of ice particle = 7.5 microns
25 :
26 :
27 : contains
28 :
29 0 : subroutine ssatcontrail_d0(state1,pbuf,dtime,ptend_loc)
30 : implicit none
31 :
32 : type(physics_state), intent(in) :: state1
33 : type(physics_ptend), intent(inout) :: ptend_loc
34 : type(physics_buffer_desc), pointer :: pbuf(:)
35 : real(r8), intent(in) :: dtime ! time step
36 : !------------------------Local storage------------------------------------------------------
37 : real(r8) :: Ma, Mh2o, epsi, Q, eta, p, G, T_contr, eslTc, eslT, RH_contr, qslTc, qslT
38 : real(r8) :: w, esiT, qsiT, ws, RH, ei
39 : integer :: i,k
40 : integer :: lchnk,ncol
41 : real(r8) :: contrail(pcols,pver), pcontrail(pcols,pver)
42 0 : real(r8), pointer, dimension(:,:) :: cld ! cloud fraction
43 0 : real(r8), pointer, dimension(:,:) :: ac_H2O
44 0 : real(r8), pointer, dimension(:,:) :: ac_SLANT_DIST
45 : integer :: itim, ifld
46 : integer :: ixcldice, ixcldliq ! indices for CLDICE and CLDLIQ
47 : integer :: ixnumice, ixnumliq
48 : real(r8):: zi, zm, rog
49 : logical :: has_aircraft_H2O
50 : logical :: has_aircraft_distance
51 : real(r8) :: hkl, hkk, tv
52 : real(r8) :: particle_mass
53 : real(r8) :: ICIWC0(pcols,pver), ICIWC, rho
54 : real(r8) :: qs
55 : real(r8) :: wght(pcols)
56 : real(r8) :: dz, ratio(pcols,pver)
57 : real(r8) :: dcld(pcols,pver)
58 : real(r8) :: ac_Q, ac_Q1, ac_Q2
59 : real(r8) :: RHcts(pcols,pver)
60 : logical :: lq(pcnst)
61 :
62 : integer :: yr, mon, day, ncsec
63 : real(r8) :: curr_factor
64 : integer :: aircraft_cnt
65 : character(len=16) :: spc_name_list(30)
66 :
67 0 : has_aircraft_H2O = .false.
68 0 : has_aircraft_distance = .false.
69 :
70 : ! Update constituents, all schemes use time split q: no tendency kept
71 0 : call cnst_get_ind('CLDICE', ixcldice, abort=.false.)
72 0 : call cnst_get_ind('CLDLIQ', ixcldliq, abort=.false.)
73 : ! Check for number concentration of cloud liquid and cloud ice (if not present)
74 : ! the indices will be set to -1)
75 0 : call cnst_get_ind('NUMICE', ixnumice, abort=.false.)
76 0 : call cnst_get_ind('NUMLIQ', ixnumliq, abort=.false.)
77 :
78 0 : call get_aircraft(aircraft_cnt, spc_name_list)
79 : !-----------------------------------------------------------------------------------------
80 : ! check if ac_H2O in namelist
81 : ! if not, then bypass this subroutine
82 : !-----------------------------------------------------------------------------------------
83 0 : if(aircraft_cnt>0) then
84 0 : do i=1,aircraft_cnt
85 0 : if(trim(spc_name_list(i)) == 'ac_H2O') then
86 0 : has_aircraft_H2O = .true.
87 : endif
88 0 : if(trim(spc_name_list(i)) == 'ac_SLANT_DIST' .or. trim(spc_name_list(i)) == 'ac_TRACK_DIST') then
89 0 : has_aircraft_distance = .true.
90 : endif
91 : enddo
92 : endif
93 : !------------------------------------------------------------------------------------------
94 0 : lq(:) = .FALSE.
95 0 : lq(1) = .TRUE.
96 0 : lq(ixcldice) = .TRUE.
97 0 : lq(ixnumice) = .TRUE.
98 :
99 0 : call physics_ptend_init(ptend_loc, state1%psetcols,'ssatcontrail',ls=.true.,lq = lq)
100 : !-----------------------------------------------------------------------------------------
101 0 : if(.not. has_aircraft_H2O) return
102 0 : if(.not. has_aircraft_distance) return
103 :
104 0 : particle_mass = 4._r8/3._r8*pi*rhoi*radius**3 ! mass of ice particle
105 :
106 0 : rog = rair/gravit ! Rd/Cp
107 :
108 :
109 0 : contrail(:,:) = 0.0_r8
110 0 : pcontrail(:,:) = 0.0_r8
111 0 : ICIWC0(:,:) = 0.0_r8
112 0 : RHcts(:,:) = 0.0_r8
113 :
114 0 : lchnk = state1%lchnk
115 0 : ncol = state1%ncol
116 :
117 0 : call get_wght_all_p(lchnk, ncol, wght)
118 :
119 0 : itim = pbuf_old_tim_idx()
120 0 : ifld = pbuf_get_index('CLD')
121 0 : call pbuf_get_field(pbuf, ifld, cld, (/1,1,itim/),(/pcols,pver,1/))
122 :
123 0 : ifld = pbuf_get_index('ac_H2O')
124 0 : call pbuf_get_field(pbuf,ifld,ac_H2O)
125 0 : ifld = pbuf_get_index('ac_SLANT_DIST')
126 0 : call pbuf_get_field(pbuf,ifld,ac_SLANT_DIST)
127 :
128 : ! adjust h2o to volume mixing ratio (mass adjustment and conversion from g/kg to kg/kg)
129 0 : Ma = mwdry
130 0 : Mh2o = mwh2o
131 :
132 : ! contrail paramter (G = CF*p/epi)
133 : ! and Schumann 1996 DOI: 10.1127/metz/5/1996/4, reprinted by Ponater 2002, JGR (eq 6-8) DOI: 10.1029/2011MS000105
134 :
135 0 : epsi = Mh2o/Ma
136 0 : ei = 1.21_r8 ! water vapor emision index (g) h2o per kg fuel (Schumann 96)
137 0 : Q = 43.e6_r8 ! specific combustion heat Schummann 1996, Q = 43 MJ/kg
138 0 : eta = 0.3_r8 ! propulsion effieciency (Ponater 2002)
139 :
140 0 : ratio = 0._r8
141 0 : dcld = 0._r8
142 :
143 0 : do i=1,ncol
144 0 : do k=1,pver
145 0 : p = (state1%pint(i,k)+state1%pint(i,k+1))/2.0_r8
146 :
147 0 : G = (ei*cpair*p)/(epsi*Q*(1.0_r8-eta)) ! eq 7, Ponater JGR 2002
148 :
149 0 : if( G > 0.053_r8 ) then
150 0 : T_contr = -46.46_r8+9.43_r8*log(G-0.053_r8)+0.72_r8*log(G-0.053_r8)*log(G-0.053_r8) ! eq 6, Ponater JGR 2002
151 0 : T_contr = T_contr + tmelt ! convert to Kelvin
152 :
153 : ! compute saturation pressure
154 0 : call qsat_water(T_contr, p, eslTc, qslTc)
155 0 : call qsat_water(state1%t(i,k), p, eslT, qslT)
156 :
157 0 : RH_contr = (G*(state1%t(i,k)-T_contr)+eslTc)/eslT
158 : ! RH_contr ranges between 0 and 1
159 0 : if(RH_contr>1.0_r8) RH_contr = 1.0_r8
160 0 : if(RH_contr<0.0_r8) RH_contr = 0.0_r8
161 :
162 0 : w = state1%q(i,k,1)/(1.0_r8-state1%q(i,k,1)) ! mixing ratio from specific humidity
163 0 : call qsat_ice(state1%t(i,k), p, esiT, qsiT)
164 0 : ws = epsi*esiT/(p-esiT) ! saturation mixing ration with respect to ice
165 0 : qs = ws/(1.0_r8+ws)
166 :
167 0 : RH = w/ws ! relative humidity with respect to ice
168 0 : if( RH>=1.0_r8 ) RHcts(i,k) = 1.0_r8
169 :
170 : ! Schumann, U. “Contrail Cirrus.” In Cirrus, edited by D. K. Lynch and others, 231–55. Oxford University Press, 2002
171 : ! IWC(g/m3) = exp(6.97+0.103*T(C))*1e-3
172 : ! IWC(kg/m3) = exp(6.97+0.103*T(C))*1e-6
173 :
174 0 : ICIWC0(i,k) = exp(6.97_r8+0.103_r8*(state1%t(i,k)-tmelt)) ! in mg/m3
175 0 : rho = p/(rair*state1%t(i,k))
176 0 : ICIWC = ICIWC0(i,k)/rho*1.0e-6_r8
177 :
178 :
179 : ! persistent contrail condition
180 0 : if( (state1%t(i,k)<T_contr).and.(RH>RH_contr).and.(RH>1.0_r8).and.(ac_H2O(i,k)>0.0_r8) ) then
181 :
182 : ! if persistent contrail, H2O emitted from aircraft turns into cloud ice
183 0 : dz = state1%zi(i,k)-state1%zi(i,k+1)
184 0 : ratio(i,k) = (ac_SLANT_DIST(i,k)*dtime*1.e4_r8)/(dz*rearth*rearth*wght(i))
185 :
186 0 : ac_Q = min(ac_H2O(i,k)*dtime + (state1%q(i,k,1)-qs)*ratio(i,k),ratio(i,k)*ICIWC)
187 0 : ptend_loc%q(i,k,ixcldice) = ac_Q/dtime
188 :
189 : ! take out water vapor from q
190 0 : ptend_loc%q(i,k,1) = -(ac_Q-ac_H2O(i,k)*dtime)/dtime
191 :
192 : ! modify cloud fraction
193 : ! by a prescribed ICIWC, we may deduce the new cloud fraction
194 :
195 0 : cld(i,k) = min(1._r8, cld(i,k)+ac_Q/ICIWC)
196 :
197 : ! modify cloud ice number concentration,
198 : ! by assuming the particle size, the number of ice particles may be obtained
199 0 : ptend_loc%q(i,k,ixnumice) = ac_Q/particle_mass/dtime
200 :
201 : else
202 : ! if not persistent contrail, just add ac_H2O to state1%q(1) (vapor phase)
203 :
204 0 : ptend_loc%q(i,k,1) = ac_H2O(i,k)
205 :
206 : endif
207 :
208 : else
209 0 : ptend_loc%q(i,k,1) = ac_H2O(i,k)
210 : end if
211 :
212 : enddo
213 :
214 : ! modify dry static energy if water field is added to any grid cell
215 : ! this bypasses geopotential_t which assumes dry static energy conservation
216 : ! water vapor added to the system is assumed to increase dry static energy
217 : ! conservation of dry static energy by geopotential_t will lower temperature to compensate
218 :
219 : zi = 0.0_r8
220 0 : do k=pver,1,-1
221 0 : hkl = state1%lnpint(i,k+1)-state1%lnpint(i,k)
222 0 : hkk = 1._r8 - state1%pint(i,k) * hkl * state1%rpdel(i,k)
223 :
224 0 : tv = state1%t(i,k) * (1._r8 + zvir*(state1%q(i,k,1)+ptend_loc%q(i,k,1)*dtime))
225 :
226 0 : zm = zi + rog * tv * hkk
227 0 : zi = zi + rog * tv * hkl
228 :
229 0 : ptend_loc%s(i,k) = (cpair*state1%t(i,k)+gravit*zm + state1%phis(i) - state1%s(i,k) )/dtime
230 : enddo
231 :
232 : enddo
233 :
234 0 : end subroutine ssatcontrail_d0
235 :
236 : end module ssatcontrail
|
