Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hcox_gc_POPs_mod.F90
7 : !
8 : ! !DESCRIPTION: Defines the HEMCO extension for the GEOS-Chem persistent
9 : ! organic pollutants (POPs) specialty simulation.
10 : !\\
11 : !\\
12 : ! !INTERFACE:
13 : !
14 : MODULE HCOX_GC_POPs_Mod
15 : !
16 : ! !USES:
17 : !
18 : USE HCO_Error_Mod
19 : USE HCO_Diagn_Mod
20 : USE HCO_State_Mod, ONLY : HCO_State ! Derived type for HEMCO state
21 : USE HCOX_State_Mod, ONLY : Ext_State ! Derived type for External state
22 :
23 : IMPLICIT NONE
24 : PRIVATE
25 : !
26 : ! !PUBLIC MEMBER FUNCTIONS:
27 : !
28 : PUBLIC :: HcoX_GC_POPs_Run
29 : PUBLIC :: HcoX_GC_POPs_Init
30 : PUBLIC :: HcoX_Gc_POPs_Final
31 : !
32 : ! !PRIVATE MEMBER FUNCTIONS:
33 : !
34 : PRIVATE :: VEGEMISPOP
35 : PRIVATE :: LAKEEMISPOP
36 : PRIVATE :: SOILEMISPOP
37 : PRIVATE :: IS_LAND
38 : PRIVATE :: IS_ICE
39 : !
40 : ! !REMARKS:
41 : !
42 : ! References:
43 : ! ============================================================================
44 : ! (1 ) Zhang, Y., and S. Tao, Global atmospheric emission inventory of
45 : ! polycyclic aromatic hydrocarbons (PAHs) for 2004. Atm Env, 43, 812-819,
46 : ! 2009.
47 : ! (2 ) Friedman, C.L, and N.E. Selin, Long-Range Atmospheric Transport of
48 : ! Polycyclic Aromatic Hydrocarbons: A Global 3-D Model Analysis
49 : ! Including Evaluation of Arctic Sources, Environ. Sci. Technol., 46(17),
50 : ! 9501-9510, 2012.
51 : ! (3 ) Friedman, C.L., Y. Zhang, and N.E. Selin, Climate change and
52 : ! emissions impacts on atmospheric PAH transport to the Arctic, Environ.
53 : ! Sci. Technol., 48, 429-437, 2014.
54 : ! (4 ) Friedman, C.L., J.R. Pierce, and N.E. Selin, Assessing the influence of
55 : ! secondary organic versus primary carbonaceous aerosols on long-range
56 : ! atmospheric polycyclic aromatic hydrocarbon transport, Environ. Sci.
57 : ! Technol., 48(6), 3293-3302, 2014.
58 : !
59 : ! !REVISION HISTORY:
60 : ! 20 Sep 2010 - N.E. Selin - Initial Version
61 : ! See https://github.com/geoschem/hemco for complete history
62 : !EOP
63 : !------------------------------------------------------------------------------
64 : !BOC
65 : !
66 : ! !DEFINED PARAMETERS:
67 : !
68 : REAL(hp), PARAMETER :: SMALLNUM = 1e-20_hp
69 : !
70 : ! !PRIVATE TYPES:
71 : !
72 : TYPE :: MyInst
73 : ! Fields required by module
74 : INTEGER :: Instance
75 : INTEGER :: ExtNr ! HEMCO Extension number
76 : INTEGER :: IDTPOPPOCPO ! Index # for POPPOC tracer
77 : INTEGER :: IDTPOPPBCPO ! Index # for POPPBC tracer
78 : INTEGER :: IDTPOPG ! Index # for POPG tracer
79 :
80 : ! Pointers to emission arrays read from disk
81 : REAL(sp), POINTER :: POP_TOT_EM(:,:) => NULL() ! [kg/m2/s]
82 : REAL(sp), POINTER :: POP_SURF(:,:) => NULL() ! [kg]
83 : REAL(sp), POINTER :: C_OC(:,:,:) => NULL() ! [kg]
84 : REAL(sp), POINTER :: C_BC(:,:,:) => NULL() ! [kg]
85 : REAL(sp), POINTER :: F_OC_SOIL(:,:) => NULL() ! [kg/m2]
86 :
87 : ! Calculated emissions of OC-phase, BC-phase, and gas-phase POPs [kg/m2/s]
88 : REAL(sp), POINTER :: EmisPOPPOCPO(:,:,:)
89 : REAL(sp), POINTER :: EmisPOPPBCPO(:,:,:)
90 : REAL(sp), POINTER :: EmisPOPG (:,:,:)
91 :
92 : ! For diagnostics
93 : REAL(sp), POINTER :: EmisPOPGfromSoil (:,:)
94 : REAL(sp), POINTER :: EmisPOPGfromLake (:,:)
95 : REAL(sp), POINTER :: EmisPOPGfromLeaf (:,:)
96 : REAL(sp), POINTER :: EmisPOPGfromSnow (:,:)
97 : REAL(sp), POINTER :: EmisPOPGfromOcean (:,:)
98 : REAL(sp), POINTER :: FluxPOPGfromSoilToAir(:,:)
99 : REAL(sp), POINTER :: FluxPOPGfromAirToSoil(:,:)
100 : REAL(sp), POINTER :: FluxPOPGfromLakeToAir(:,:)
101 : REAL(sp), POINTER :: FluxPOPGfromAirToLake(:,:)
102 : REAL(sp), POINTER :: FluxPOPGfromLeafToAir(:,:)
103 : REAL(sp), POINTER :: FluxPOPGfromAirtoLeaf(:,:)
104 : REAL(sp), POINTER :: FugacitySoilToAir (:,:)
105 : REAL(sp), POINTER :: FugacityLakeToAir (:,:)
106 : REAL(sp), POINTER :: FugacityLeafToAir (:,:)
107 :
108 : TYPE(MyInst), POINTER :: NextInst => NULL()
109 : END TYPE MyInst
110 :
111 : ! Pointer to instances
112 : TYPE(MyInst), POINTER :: AllInst => NULL()
113 :
114 : CONTAINS
115 : !EOC
116 : !------------------------------------------------------------------------------
117 : ! Harmonized Emissions Component (HEMCO) !
118 : !------------------------------------------------------------------------------
119 : !BOP
120 : !
121 : ! !IROUTINE: HCOX_GC_POPs_run
122 : !
123 : ! !DESCRIPTION: Subroutine HcoX\_Gc\_POPs\_Run computes emissions of OC-phase,
124 : ! BC-phase, and gas-phase POPs for the GEOS-Chem POPs specialty simulation.
125 : !\\
126 : !\\
127 : ! !INTERFACE:
128 : !
129 0 : SUBROUTINE HCOX_GC_POPs_Run( ExtState, HcoState, RC )
130 : !
131 : ! !USES:
132 : !
133 : USE HCO_EmisList_Mod, ONLY : HCO_GetPtr
134 : USE HCO_FluxArr_Mod, ONLY : HCO_EmisAdd
135 : USE HCO_Clock_Mod, ONLY : HcoClock_First
136 : !
137 : ! !INPUT PARAMETERS:
138 : !
139 : TYPE(Ext_State), POINTER :: ExtState ! Options for POPs sim
140 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
141 : !
142 : ! !INPUT/OUTPUT PARAMETERS:
143 : !
144 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
145 : !
146 : ! !REMARKS:
147 : ! This code is based on routine EMISSPOPS in prior versions of GEOS-Chem.
148 : !
149 : ! !REVISION HISTORY:
150 : ! 20 Sep 2010 - N.E. Selin - Initial Version based on EMISSMERCURY
151 : ! See https://github.com/geoschem/hemco for complete history
152 : !EOP
153 : !------------------------------------------------------------------------------
154 : !BOC
155 : !
156 : ! !DEFINED PARAMETERS:
157 : !
158 : ! Universal gas constant for adjusting KOA for temp: 8.3144598 [J/mol/K]
159 : REAL(hp), PARAMETER :: R = 8.3144598e+0_hp
160 :
161 : ! Density of octanol, needed for partitioning into OC: 820 [kg/m^3]
162 : REAL(hp), PARAMETER :: DENS_OCT = 82e+1_hp
163 :
164 : ! Density of BC, needed for partitioning onto BC: 1 [kg/L] or 1000 [kg/m^3]
165 : ! From Lohmann and Lammel, Environ. Sci. Technol., 2004, 38:3793-3803.
166 : REAL(hp), PARAMETER :: DENS_BC = 1e+3_hp
167 : !
168 : ! !LOCAL VARIABLES:
169 : !
170 : INTEGER :: I, J, L
171 : INTEGER :: PBL_MAX
172 : INTEGER :: MONTH, YEAR
173 : REAL(hp) :: F_OF_PBL, TK
174 : REAL(hp) :: T_POP
175 : REAL(hp) :: C_OC1, C_BC1
176 : REAL(hp) :: C_OC2, C_BC2
177 : REAL(hp) :: F_POP_OC, F_POP_BC
178 : REAL(hp) :: F_POP_G, AIR_VOL
179 : REAL(hp) :: KOA_T, KBC_T
180 : REAL(hp) :: KOC_BC_T, KBC_OC_T
181 : REAL(hp) :: VR_OC_AIR, VR_OC_BC
182 : REAL(hp) :: VR_BC_AIR, VR_BC_OC
183 : REAL(hp) :: SUM_F
184 : REAL(hp) :: OC_AIR_RATIO, OC_BC_RATIO
185 : REAL(hp) :: BC_AIR_RATIO, BC_OC_RATIO
186 : REAL(hp) :: FRAC_SNOW_OR_ICE, FRAC_SNOWFREE_LAND
187 : REAL(hp) :: FRAC_LEAF, FRAC_LAKE, FRAC_SOIL
188 : ! LOGICAL, SAVE :: FIRST = .TRUE.
189 : LOGICAL :: aIR
190 : LOGICAL :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE
191 : CHARACTER(LEN=255) :: MSG, LOC
192 :
193 : ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer
194 : ! from gas phase to OC. For now we use Delta H for phase transfer
195 : ! from the gas phase to the pure liquid state.
196 : ! For PHENANTHRENE:
197 : ! this is taken as the negative of the Delta H for phase transfer
198 : ! from the pure liquid state to the gas phase (Schwarzenbach,
199 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol].
200 : ! For PYRENE:
201 : ! this is taken as the negative of the Delta H for phase transfer
202 : ! from the pure liquid state to the gas phase (Schwarzenbach,
203 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol].
204 : ! For BENZO[a]PYRENE:
205 : ! this is also taken as the negative of the Delta H for phase transfer
206 : ! from the pure liquid state to the gas phase (Schwarzenbach,
207 : ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol]
208 : REAL(hp) :: DEL_H
209 :
210 : ! KOA_298 for partitioning of gas phase POP to atmospheric OC
211 : ! KOA_298 = Cpop in octanol/Cpop in atmosphere at 298 K
212 : ! For PHENANTHRENE:
213 : ! log KOA_298 = 7.64, or 4.37*10^7 [unitless]
214 : ! For PYRENE:
215 : ! log KOA_298 = 8.86, or 7.24*10^8 [unitless]
216 : ! For BENZO[a]PYRENE:
217 : ! log KOA_298 = 11.48, or 3.02*10^11 [unitless]
218 : ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825).
219 : REAL(hp) :: KOA_298
220 :
221 : ! KBC_298 for partitioning of gas phase POP to atmospheric BC
222 : ! KBC_298 = Cpop in black carbon/Cpop in atmosphere at 298 K
223 : ! For PHENANTHRENE:
224 : ! log KBC_298 = 10.0, or 1.0*10^10 [unitless]
225 : ! For PYRENE:
226 : ! log KBC_298 = 11.0, or 1.0*10^11 [unitless]
227 : ! For BENZO[a]PYRENE:
228 : ! log KBC_298 = 13.9, or 7.94*10^13 [unitless]
229 : ! (Lohmann and Lammel, EST, 2004, 38:3793-3802)
230 : REAL(hp) :: KBC_298
231 :
232 : ! Pointers
233 0 : REAL(sp), POINTER :: Arr3D(:,:,:)
234 : TYPE(MyInst), POINTER :: Inst
235 :
236 : !=======================================================================
237 : ! HCOX_GC_POPs_RUN begins here!
238 : !=======================================================================
239 0 : LOC = 'HCOX_GC_POPs_RUN (HCOX_GC_POPS_MOD.F90)'
240 :
241 : ! Return if extension not turned on
242 0 : IF ( ExtState%GC_POPs <= 0 ) RETURN
243 :
244 : ! Enter
245 0 : CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
246 0 : IF ( RC /= HCO_SUCCESS ) THEN
247 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
248 0 : RETURN
249 : ENDIF
250 :
251 : ! Get instance
252 0 : Inst => NULL()
253 0 : CALL InstGet ( ExtState%GC_POPs, Inst, RC )
254 0 : IF ( RC /= HCO_SUCCESS ) THEN
255 0 : WRITE(MSG,*) 'Cannot find GC_POPs instance Nr. ', ExtState%GC_POPs
256 0 : CALL HCO_ERROR(MSG,RC)
257 0 : RETURN
258 : ENDIF
259 :
260 0 : DEL_H = ExtState%POP_DEL_H
261 0 : KOA_298 = ExtState%POP_KOA
262 0 : KBC_298 = ExtState%POP_KBC
263 0 : Arr3D => NULL()
264 :
265 : !=======================================================================
266 : ! Get pointers to gridded data imported through config. file
267 : !=======================================================================
268 0 : IF ( HcoClock_First(HcoState%Clock,.TRUE.) ) THEN
269 :
270 0 : CALL HCO_GetPtr( HcoState, 'TOT_POP', Inst%POP_TOT_EM, RC )
271 0 : IF ( RC /= HCO_SUCCESS ) THEN
272 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
273 0 : RETURN
274 : ENDIF
275 :
276 0 : CALL HCO_GetPtr( HcoState, 'GLOBAL_OC', Inst%C_OC, RC )
277 0 : IF ( RC /= HCO_SUCCESS ) THEN
278 0 : CALL HCO_ERROR( 'ERROR 2', RC, THISLOC=LOC )
279 0 : RETURN
280 : ENDIF
281 :
282 0 : CALL HCO_GetPtr( HcoState, 'GLOBAL_BC', Inst%C_BC, RC )
283 0 : IF ( RC /= HCO_SUCCESS ) THEN
284 0 : CALL HCO_ERROR( 'ERROR 3', RC, THISLOC=LOC )
285 0 : RETURN
286 : ENDIF
287 :
288 0 : CALL HCO_GetPtr( HcoState, 'SURF_POP', Inst%POP_SURF, RC )
289 0 : IF ( RC /= HCO_SUCCESS ) THEN
290 0 : CALL HCO_ERROR( 'ERROR 4', RC, THISLOC=LOC )
291 0 : RETURN
292 : ENDIF
293 :
294 0 : CALL HCO_GetPtr( HcoState, 'SOIL_CARBON', Inst%F_OC_SOIL, RC )
295 0 : IF ( RC /= HCO_SUCCESS ) THEN
296 0 : CALL HCO_ERROR( 'ERROR 5', RC, THISLOC=LOC )
297 0 : RETURN
298 : ENDIF
299 :
300 : ! Convert F_OC_SOIL from kg/m2 to fraction
301 0 : DO J=1, HcoState%NY
302 0 : DO I=1, HcoState%NX
303 :
304 : ! Assume most of carbon mass extends to 5 cm and calculate
305 : ! concentration in kg/kg
306 : ! For now, assume a mean soil bulk density of 1300 kg/m3 similar to
307 : ! McLachlan 2002 to calculate a dry weight fraction
308 0 : Inst%F_OC_SOIL(I,J) = Inst%F_OC_SOIL(I,J) / 30e-2_hp / 13e+2_hp
309 :
310 : ENDDO
311 : ENDDO
312 :
313 : ! FIRST = .FALSE.
314 :
315 : ENDIF
316 :
317 : ! Maximum extent of the PBL [model level]
318 0 : IF ( .NOT. ASSOCIATED(ExtState%PBL_MAX) ) THEN
319 0 : CALL HCO_ERROR ( 'PBL_MAX not defined in ExtState!', RC )
320 0 : RETURN
321 : ELSE
322 0 : PBL_MAX = DBLE( ExtState%PBL_MAX )
323 : ENDIF
324 :
325 : !=================================================================
326 : ! Call emissions routines for revolatilization fluxes from surfaces
327 : ! Assume all re-emited POPs are in the gas phase until partitioning
328 : ! to ambient OC and BC in the boundary layer.
329 : ! Re-emission flux/mass depends on type of POP
330 : ! First draft, CLF, 28 Aug 2012
331 : !=================================================================
332 0 : CALL SOILEMISPOP( Inst%POP_SURF, Inst%F_OC_SOIL, HcoState, ExtState, Inst )
333 0 : CALL LAKEEMISPOP( Inst%POP_SURF, HcoState, ExtState, Inst )
334 0 : CALL VEGEMISPOP ( Inst%POP_SURF, HcoState, ExtState, Inst )
335 :
336 0 : Inst%EmisPOPGfromSnow = 0e+0_hp
337 0 : Inst%EmisPOPGfromOcean = 0e+0_hp
338 :
339 : ! Loop over grid boxes
340 0 : DO J = 1, HcoState%Ny
341 0 : DO I = 1, HcoState%Nx
342 :
343 0 : F_OF_PBL = 0e+0_hp
344 0 : T_POP = 0e+0_hp
345 :
346 : ! Here, save the total from the emissions array
347 : ! into the T_POP variable [kg/m2/s]
348 0 : T_POP = Inst%POP_TOT_EM(I,J)
349 :
350 : ! Now add revolatilization (secondary) emissions to primary [kg/m2/s]
351 0 : T_POP = T_POP + Inst%EmisPOPGfromSoil(I,J) + &
352 0 : Inst%EmisPOPGfromLake(I,J) + &
353 0 : Inst%EmisPOPGfromLeaf(I,J) !+ &
354 : !Inst%EmisPOPGfromSnow(I,J) + &
355 : !Inst%EmisPOPGfromOcean(I,J)
356 :
357 : !====================================================================
358 : ! Apportion total POPs emitted to gas phase, OC-bound, and BC-bound
359 : ! emissions (clf, 2/1/2011)
360 : ! Then partition POP throughout PBL; store into STT [kg]
361 : ! Now make sure STT does not underflow (cdh, bmy, 4/6/06; eck 9/20/10)
362 : !====================================================================
363 :
364 : ! Loop up to max PBL level
365 0 : DO L = 1, PBL_MAX
366 :
367 : ! Get temp [K]
368 0 : TK = ExtState%TK%Arr%Val(I,J,L)
369 :
370 : ! Define temperature-dependent partition coefficients:
371 : ! KOA_T, the octanol-air partition coeff at temp T [unitless]
372 0 : KOA_T = KOA_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK) - (1e+0_hp/298e+0_hp)))
373 :
374 : ! Define KBC_T, the BC-air partition coeff at temp T [unitless]
375 : ! TURN OFF TEMPERATURE DEPENDENCY FOR SENSITIVITY ANALYSIS
376 0 : KBC_T = KBC_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK) - (1e+0_hp/298e+0_hp)))
377 :
378 : ! Define KOC_BC_T, theoretical OC-BC part coeff at temp T [unitless]
379 0 : KOC_BC_T = KOA_T / KBC_T
380 :
381 : ! Define KBC_OC_T, theoretical BC_OC part coeff at temp T [unitless]
382 0 : KBC_OC_T = 1d0 / KOC_BC_T
383 :
384 : ! Get monthly mean OC and BC concentrations [kg/box]
385 0 : C_OC1 = Inst%C_OC(I,J,L)
386 0 : C_BC1 = Inst%C_BC(I,J,L)
387 :
388 : ! Make sure OC is not negative
389 0 : C_OC1 = MAX( C_OC1, 0e+0_hp )
390 :
391 : ! Convert C_OC and C_BC units to volume per box
392 : ! [m^3 OC or BC/box]
393 0 : C_OC2 = C_OC1 / DENS_OCT
394 0 : C_BC2 = C_BC1 / DENS_BC
395 :
396 : ! Get air volume (m^3)
397 0 : AIR_VOL = ExtState%AIRVOL%Arr%Val(I,J,L)
398 :
399 : ! Define volume ratios:
400 : ! VR_OC_AIR = volume ratio of OC to air [unitless]
401 0 : VR_OC_AIR = C_OC2 / AIR_VOL
402 :
403 : ! VR_OC_BC = volume ratio of OC to BC [unitless]
404 0 : VR_OC_BC = C_OC2 / C_BC2
405 :
406 : ! VR_BC_AIR = volume ratio of BC to air [unitless]
407 0 : VR_BC_AIR = VR_OC_AIR / VR_OC_BC
408 :
409 : ! VR_BC_OC = volume ratio of BC to OC [unitless]
410 : !VR_BC_OC(I,J,L) = 1d0 / VR_OC_BC(I,J,L)
411 0 : VR_BC_OC = 1d0 / VR_OC_BC
412 :
413 : ! Redefine fractions of total POPs in box (I,J,L) that are OC-phase,
414 : ! BC-phase, and gas phase with new time step (should only change if
415 : ! temp changes or OC/BC concentrations change)
416 0 : OC_AIR_RATIO = 1e+0_hp / (KOA_T * VR_OC_AIR)
417 0 : OC_BC_RATIO = 1e+0_hp / (KOC_BC_T * VR_OC_BC)
418 :
419 0 : BC_AIR_RATIO = 1e+0_hp / (KBC_T * VR_BC_AIR)
420 0 : BC_OC_RATIO = 1e+0_hp / (KBC_OC_T * VR_BC_OC)
421 :
422 : ! If there are zeros in OC or BC concentrations, make sure they
423 : ! don't cause problems with phase fractions
424 0 : IF ( C_OC1 > SMALLNUM .and. C_BC1 > SMALLNUM ) THEN
425 0 : F_POP_OC = 1e+0_hp / (1e+0_hp + OC_AIR_RATIO + OC_BC_RATIO)
426 0 : F_POP_BC = 1e+0_hp / (1e+0_hp + BC_AIR_RATIO + BC_OC_RATIO)
427 :
428 0 : ELSE IF ( C_OC1 > SMALLNUM .and. C_BC1 .le. SMALLNUM ) THEN
429 0 : F_POP_OC = 1e+0_hp / (1e+0_hp + OC_AIR_RATIO)
430 0 : F_POP_BC = SMALLNUM
431 :
432 0 : ELSE IF ( C_OC1 .le. SMALLNUM .and. C_BC1 > SMALLNUM ) THEN
433 0 : F_POP_OC = SMALLNUM
434 0 : F_POP_BC = 1e+0_hp / (1e+0_hp + BC_AIR_RATIO)
435 :
436 0 : ELSE IF ( C_OC1 .le. SMALLNUM .and. C_BC1 .le. SMALLNUM) THEN
437 0 : F_POP_OC = SMALLNUM
438 0 : F_POP_BC = SMALLNUM
439 : ENDIF
440 :
441 : ! Gas-phase:
442 0 : F_POP_G = 1e+0_hp - F_POP_OC - F_POP_BC
443 :
444 : ! Check that sum of fractions equals 1
445 0 : SUM_F = F_POP_OC + F_POP_BC + F_POP_G
446 :
447 : ! Fraction of PBL that box (I,J,L) makes up [unitless]
448 0 : F_OF_PBL = ExtState%FRAC_OF_PBL%Arr%Val(I,J,L)
449 :
450 : ! Calculate rates of POP emissions in each phase [kg/m2/s]
451 : ! OC-phase:
452 0 : Inst%EmisPOPPOCPO(I,J,L) = F_POP_OC * F_OF_PBL * T_POP
453 :
454 : ! BC-phase
455 0 : Inst%EmisPOPPBCPO(I,J,L) = F_POP_BC * F_OF_PBL * T_POP
456 :
457 : ! Gas-phase
458 0 : Inst%EmisPOPG(I,J,L) = F_POP_G * F_OF_PBL * T_POP
459 :
460 : ENDDO
461 :
462 : ENDDO
463 : ENDDO
464 :
465 : !=======================================================================
466 : ! Add POPs emissions to HEMCO data structure & diagnostics
467 : !=======================================================================
468 :
469 : !----------------------
470 : ! OC-PHASE EMISSIONS
471 : !----------------------
472 0 : IF ( Inst%IDTPOPPOCPO > 0 ) THEN
473 :
474 : ! Add flux to emissions array
475 0 : Arr3D => Inst%EmisPOPPOCPO(:,:,:)
476 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPPOCPO, &
477 0 : RC, ExtNr=Inst%ExtNr )
478 0 : Arr3D => NULL()
479 0 : IF ( RC /= HCO_SUCCESS ) THEN
480 0 : CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPPOCPO', RC )
481 0 : RETURN
482 : ENDIF
483 : ENDIF
484 :
485 : !----------------------
486 : ! BC-PHASE EMISSIONS
487 : !----------------------
488 0 : IF ( Inst%IDTPOPPBCPO > 0 ) THEN
489 :
490 : ! Add flux to emissions array
491 0 : Arr3D => Inst%EmisPOPPBCPO(:,:,:)
492 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPPBCPO, &
493 0 : RC, ExtNr=Inst%ExtNr )
494 0 : Arr3D => NULL()
495 0 : IF ( RC /= HCO_SUCCESS ) THEN
496 0 : CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPPBCPO', RC )
497 0 : RETURN
498 : ENDIF
499 : ENDIF
500 :
501 : !----------------------
502 : ! GASEOUS EMISSIONS
503 : !----------------------
504 0 : IF ( Inst%IDTPOPG > 0 ) THEN
505 :
506 : ! Add flux to emissions array
507 0 : Arr3D => Inst%EmisPOPG(:,:,:)
508 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTPOPG, &
509 0 : RC, ExtNr=Inst%ExtNr )
510 0 : Arr3D => NULL()
511 0 : IF ( RC /= HCO_SUCCESS ) THEN
512 0 : CALL HCO_ERROR( 'HCO_EmisAdd error: EmisPOPG', RC )
513 0 : RETURN
514 : ENDIF
515 :
516 : ENDIF
517 :
518 : !=======================================================================
519 : ! Cleanup & quit
520 : !=======================================================================
521 :
522 : ! Nullify pointers
523 0 : Inst => NULL()
524 :
525 : ! Return w/ success
526 0 : CALL HCO_LEAVE( HcoState%Config%Err,RC )
527 :
528 0 : END SUBROUTINE HCOX_GC_POPs_Run
529 : !EOC
530 : !------------------------------------------------------------------------------
531 : ! GEOS-Chem Global Chemical Transport Model !
532 : !------------------------------------------------------------------------------
533 : !BOP
534 : !
535 : ! !IROUTINE: soilemispop
536 : !
537 : ! !DESCRIPTION: Subroutine SOILEMISPOP is the subroutine for secondary
538 : ! POP emissions from soils.
539 : !\\
540 : !\\
541 : ! !INTERFACE:
542 : !
543 0 : SUBROUTINE SOILEMISPOP( POP_SURF, F_OC_SOIL, HcoState, ExtState, Inst )
544 : !
545 : ! !INPUT PARAMETERS:
546 : !
547 : REAL(sp), DIMENSION(:,:), INTENT(IN) :: POP_SURF ! POP sfc conc [kg]
548 : REAL(sp), DIMENSION(:,:), INTENT(IN) :: F_OC_SOIL ! Frac C in soil [g/g]
549 : TYPE(HCO_STATE), POINTER :: HcoState ! Hemco state
550 : TYPE(Ext_State), POINTER :: ExtState ! Module options
551 : TYPE(MyInst), POINTER :: Inst ! Instance
552 : !
553 : ! !REVISION HISTORY:
554 : ! 21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD
555 : ! See https://github.com/geoschem/hemco for complete history
556 : !EOP
557 : !------------------------------------------------------------------------------
558 : !BOC
559 : !
560 : ! !LOCAL VARIABLES:
561 : !
562 : LOGICAL :: IS_SNOW_OR_ICE
563 : LOGICAL :: IS_LAND_OR_ICE
564 : INTEGER :: I, J, L
565 : REAL(hp) :: POPG
566 : REAL(hp) :: TK_SURF
567 : REAL(hp) :: KSA_T, FLUX, F_OC
568 : REAL(hp) :: SOIL_CONC, DTSRCE, KOA_T
569 : REAL(hp) :: DIFF, FSOIL, FAIR, DS
570 : REAL(hp) :: DSA, DAD, DWD, PL
571 : REAL(hp) :: ZSOIL, ZAIR, TK
572 : REAL(hp) :: DTCHEM, NEWSOIL, E_KDEG
573 : REAL(hp) :: FRAC_LAKE, FRAC_SOIL
574 : REAL(hp) :: FUG_R
575 : REAL(hp) :: AREA_M2
576 : !
577 : ! !DEFINED PARAMETERS:
578 : !
579 : ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer
580 : ! from gas phase to OC. For now we use Delta H for phase transfer
581 : ! from the gas phase to the pure liquid state.
582 : ! For PHENANTHRENE:
583 : ! this is taken as the negative of the Delta H for phase transfer
584 : ! from the pure liquid state to the gas phase (Schwarzenbach,
585 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol].
586 : ! For PYRENE:
587 : ! this is taken as the negative of the Delta H for phase transfer
588 : ! from the pure liquid state to the gas phase (Schwarzenbach,
589 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol].
590 : ! For BENZO[a]PYRENE:
591 : ! this is also taken as the negative of the Delta H for phase transfer
592 : ! from the pure liquid state to the gas phase (Schwarzenbach,
593 : ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol]
594 : REAL(hp) :: DEL_H
595 :
596 : ! R = universal gas constant for adjusting KOA for temp:
597 : ! 8.3144598 [J/mol/K OR m3*Pa/K/mol]
598 : REAL(hp), PARAMETER :: R = 8.3144598d0
599 :
600 : ! Molecular weight
601 : ! For phe, 0.17823 kg/mol
602 : ! For pyr, 0.20225 kg/mol
603 : ! For BaP, 0,25231 kg/mol
604 : REAL(hp) :: MW
605 :
606 : ! Molecular weight of air
607 : REAL(hp), PARAMETER :: MWAIR = 28.97d0 ! g/mol
608 :
609 : ! For PHENANTHRENE:
610 : ! log KOA_298 = 7.64, or 4.37*10^7 [unitless]
611 : ! For PYRENE:
612 : ! log KOA_298 = 8.86, or 7.24*10^8 [unitless]
613 : ! For BENZO[a]PYRENE:
614 : ! log KOA_298 = 11.48, or 3.02*10^11 [unitless]
615 : ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825).
616 : REAL(hp) :: KOA_298
617 :
618 : ! Set transfer velocity and diffusion coefficient values
619 : REAL(hp), PARAMETER :: KSA = 1d0 ![m/h]
620 : REAL(hp), PARAMETER :: BA = 0.04d0 ![m2/h]
621 : REAL(hp), PARAMETER :: BW = 4d-6 ![m2/h]
622 :
623 : ! Set soil degradation rate
624 : REAL(hp), PARAMETER :: DEGR = 3.5d-5 ![/h]
625 :
626 : !=================================================================
627 : ! SOILEMISPOP begins here!
628 : !=================================================================
629 :
630 : ! Copy values from ExtState
631 0 : DEL_H = ExtState%POP_DEL_H
632 0 : KOA_298 = ExtState%POP_KOA
633 0 : MW = ExtState%POP_XMW
634 :
635 : ! Emission timestep [s]
636 0 : DTSRCE = HcoState%TS_EMIS
637 :
638 : ! Chemistry timestep [h]
639 0 : DTCHEM = HcoState%TS_CHEM / 60e+0_hp
640 :
641 0 : DO J=1, HcoState%NY
642 0 : DO I=1, HcoState%NX
643 :
644 : ! Set logicals
645 : ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE)
646 : ! IS_LAND will return non-ocean boxes but may still contain lakes
647 : ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE)
648 : IS_LAND_OR_ICE = ( ( IS_LAND(I,J,ExtState) ) .OR. &
649 0 : ( IS_ICE (I,J,ExtState) ) )
650 : IS_SNOW_OR_ICE = ( ( IS_ICE (I,J,ExtState) ) .OR. &
651 0 : ( IS_LAND(I,J,ExtState) .AND. &
652 0 : ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) )
653 :
654 : ! Do soils routine only if we are on land that is not covered with
655 : ! snow or ice
656 0 : IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN
657 :
658 : ! Get fraction of grid box covered by lake surface area
659 0 : FRAC_LAKE = ExtState%FRLAKE%Arr%Val(I,J)
660 :
661 : ! Get fraction of land remaining
662 : ! Assume the remaining land is soil and get OC content.
663 : ! If remaining land is not soil (e.g., desert), there
664 : ! should be a characteristically low OC content
665 : ! that will have little capacity to store POPs
666 : ! ONLY SUBTRACT FRAC LAKE NOW
667 0 : FRAC_SOIL = MAX(1e+0_hp - FRAC_LAKE, 0e+0_hp)
668 :
669 : ! Get surface skin temp [K]
670 0 : TK_SURF = ExtState%TSKIN%Arr%Val(I,J)
671 :
672 : ! Get air temp [K]
673 0 : TK = ExtState%TK%Arr%Val(I,J,1)
674 :
675 : ! Get gas phase air POP concentration at surface in mol/m3
676 : ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15)
677 0 : POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM )
678 :
679 : ! (kg trc/kg dry air) / (0.178 kg trc/mol) * (kg dry air/m3)
680 0 : POPG = POPG / MW * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3
681 :
682 : ! Grid box surface area [m2]
683 0 : AREA_M2 = HcoState%Grid%AREA_M2%Val(I,J)
684 :
685 : ! Get soil concentration in top 5 cm of soil
686 : ! (following Howsam et al 2000)
687 : ! From Howsam et al, soil burdens are equal to
688 : ! 2.6 years deposition for PHE,
689 : ! 10 years for PYR, and 9.4 years for BaP
690 : ! Convert to mol/m3
691 : ! 2.6 yrs * kg deposited to soil in 1 yr * / 0.178 kg/mol
692 : ! / area grid box (m2) / 0.05 m
693 0 : SOIL_CONC = 10e+0_hp * POP_SURF(I,J) / MW / &
694 0 : AREA_M2 / 5e-2_hp ! mol/m3
695 :
696 : ! Get rid of mass due to degradation
697 : ! Use rate constant for BaP from Mackay and Paterson 1991:
698 : ! 3.5*10^-5 /h
699 : ! Calculate exponential factor
700 0 : E_KDEG = EXP (-DEGR * DTCHEM)
701 :
702 : ! Adjust conc
703 0 : NEWSOIL = SOIL_CONC * E_KDEG
704 :
705 : ! Get foc from GTMM saved files
706 0 : F_OC = F_OC_SOIL(I,J)
707 :
708 : ! Define temperature-dependent KOA:
709 : KOA_T = KOA_298 * EXP((-DEL_H/R) * ( ( 1e+0_hp / TK_SURF ) - &
710 0 : ( 1e+0_hp / 298e+0_hp ) ))
711 :
712 : ! Dimensionless coefficient (mol/m3 soil / mol/m3 air)
713 : ! KSA = 1.5 (fTOC)*Koa
714 0 : KSA_T = 1.5 * F_OC * KOA_T
715 0 : KSA_T = MAX( KSA_T, SMALLNUM )
716 :
717 : ! Calculate fugacities from concentrations by dividing by "Z"
718 : ! values, or the fugacity capacity in mol/m3*Pa following
719 : ! Mackay and Paterson, 1991
720 :
721 : ! fsoil = Csoil * R * T / KSA [Pa]
722 : ! where Csoil is in mol/m3, R is in Pa * m3 / mol K
723 : ! T is in K and KSA is dimensionless
724 0 : FSOIL = NEWSOIL * R * TK_SURF / KSA_T
725 :
726 : ! fair = Cair * R * T [Pa]
727 : ! where Cair is in mol/m3, R is in Pa * m3 / mol K and T is in K
728 0 : FAIR = POPG * R * TK
729 :
730 : ! Calculate the fugacity gradient [Pa]
731 : ! If the gradient is negative, fair is larger and the POP will
732 : ! diffuse from air to soil
733 : ! If the gradient is positive, fsoil is larger and POP will
734 : ! diffuse from soil to air
735 0 : DIFF = FSOIL - FAIR
736 0 : FUG_R = FSOIL/FAIR
737 :
738 : ! Calculate "Z" values from fugacities.
739 : ! Z is the fugacity capacity in mol/m3*Pa. C = Z*f, so Z = C/f
740 0 : ZAIR = POPG / FAIR ! (mol/m3) / (Pa)
741 0 : ZSOIL = NEWSOIL / FSOIL ! (mol/m3) / (Pa)
742 :
743 : ! Calculate the "D" value, or the transfer coefficient that
744 : ! describes the movement of POP between phases (Mackay and
745 : ! Paterson, 1991). [mol/h*Pa]
746 : ! The D value for soil-air diffusion is given by
747 : ! Ds = 1 / (1/Dsa + 1/(Dad + Dwd))
748 : ! Dsa is the air-side boundary layer diffusion parameter [mol/h*Pa]
749 : ! Dad is the diffusion parameter between soil particles and
750 : ! "soil air" [mol/h*Pa]
751 : ! Dwd is the diffusion parameter between soil particles and
752 : ! porewater [mol/h*Pa]
753 : ! Dsa is in series with soil-air and soil-water diffusion,
754 : ! which are in parallel
755 :
756 : ! Need to define each D value
757 : ! DSA = kSA * Zair
758 : ! where kSA is a mass transfer coefficient [m/h],
759 : ! Zair is the air fugacity capacity [mol/m3*Pa]
760 : ! ***********
761 : ! DAD = BA * Zair
762 : ! where BA is the molecular diffusivity in air [m2/h]
763 : ! ***********
764 : ! DWD = BW * Zsoil
765 : ! where BW is the molecular diffusivity in water [m2/h]
766 : ! **** PL = the soil diffusion pathlength, set to half the
767 : ! soil depth (0.025 m)
768 0 : DSA = KSA * ZAIR ! (m/h) * (mol/m3*Pa) = mol/m2*h*Pa
769 0 : DAD = BA* ZAIR ! (m2/h) * (mol/m3*Pa) = mol/m*h*Pa
770 0 : DWD = BW * ZSOIL ! (m2/h) * (mol/m3*Pa) = mol/m*h*Pa
771 0 : PL = 0.025e+0_hp
772 0 : DS = 1e+0_hp / ( 1e+0_hp/DSA + PL/(DAD+DWD) ) ! mol/(m2*h*Pa) [* m3*Pa/K/mol = m/h/K
773 :
774 : ! Calculate Flux in mol/m2/h
775 0 : FLUX = DS * DIFF
776 :
777 : ! Change to units of ng/m2/d for storage
778 0 : FLUX = FLUX * 24e+0_hp * MW * 1e+12_hp
779 :
780 : ! Kludge soil emissions from poles for now
781 : ! Bug somewhere that allows GCAP versions to think some high polar
782 : ! boxes during some months are land rather than ice - results in
783 : ! extremely high fluxes
784 0 : IF ( HcoState%Grid%YMID%Val(I,J) > 60 .OR. &
785 : HcoState%Grid%YMID%Val(I,J) < -60 ) THEN
786 0 : FLUX = 0e+0_hp
787 0 : DIFF = 0e+0_hp
788 : ENDIF
789 :
790 : ! Convert to an emission rate in kg/m2/s for returning to
791 : ! HcoX_GC_POPs_Run
792 0 : Inst%EmisPOPGfromSoil(I,J) = MAX( FLUX / 24e+0_hp / 3600e+0_hp / &
793 0 : 1e+12_hp, 0e+0_hp )
794 :
795 : ! Multiply the emissions by the fraction of land that is soil
796 0 : Inst%EmisPOPGfromSoil(I,J) = Inst%EmisPOPGfromSoil(I,J) * FRAC_SOIL
797 :
798 : ! If the flux is positive, then the direction will be from the
799 : ! soil to the air.
800 : ! If the flux is zero or negative, store it in a separate array.
801 0 : IF ( FLUX > 0e+0_hp ) THEN
802 :
803 : ! Store positive flux
804 0 : Inst%FluxPOPGfromSoilToAir(I,J) = FLUX
805 :
806 : ! Make sure negative flux diagnostic has nothing added to it
807 0 : Inst%FluxPOPGfromAirToSoil(I,J) = 0e+0_hp
808 :
809 : ! Store the soil/air fugacity ratio
810 0 : Inst%FugacitySoilToAir(I,J) = FUG_R
811 :
812 0 : ELSE IF ( FLUX <= 0e+0_hp ) THEN
813 :
814 : ! Store the negative flux
815 0 : Inst%FluxPOPGfromAirToSoil(I,J) = FLUX
816 :
817 : ! Add nothing to positive flux
818 0 : Inst%FluxPOPGfromSoilToAir(I,J) = 0e+0_hp
819 :
820 : ! Continue to store the fugacity ratio
821 0 : Inst%FugacitySoilToAir(I,J) = FUG_R
822 :
823 : ENDIF
824 :
825 : ELSE
826 :
827 : ! We are not on land or the land is covered with ice or snow
828 0 : FLUX = 0e+0_hp
829 0 : Inst%EmisPOPGfromSoil(I,J) = 0e+0_hp
830 0 : Inst%FluxPOPGfromSoilToAir(I,J) = 0e+0_hp
831 0 : Inst%FluxPOPGfromAirToSoil(I,J) = 0e+0_hp
832 0 : Inst%FugacitySoilToAir(I,J) = 0e+0_hp
833 :
834 : ENDIF
835 :
836 : ENDDO
837 : ENDDO
838 :
839 0 : END SUBROUTINE SOILEMISPOP
840 : !EOC
841 :
842 : !------------------------------------------------------------------------------
843 : ! GEOS-Chem Global Chemical Transport Model !
844 : !------------------------------------------------------------------------------
845 : !BOP
846 : !
847 : ! !IROUTINE: lakeemispop
848 : !
849 : ! !DESCRIPTION: Subroutine LAKEEMISPOP is the subroutine for secondary
850 : ! POP emissions from lakes.
851 : !\\
852 : !\\
853 : ! !INTERFACE:
854 : !
855 0 : SUBROUTINE LAKEEMISPOP( POP_SURF, HcoState, ExtState, Inst )
856 : !
857 : ! !INPUT PARAMETERS:
858 : !
859 : REAL(sp), DIMENSION(:,:), INTENT(IN) :: POP_SURF ! POP sfc conc [kg]
860 : TYPE(HCO_STATE), POINTER :: HcoState ! Hemco state
861 : TYPE(Ext_State), POINTER :: ExtState ! Module options
862 : TYPE(MyInst), POINTER :: Inst ! Instance
863 : !
864 : ! !REVISION HISTORY:
865 : ! 21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD
866 : ! See https://github.com/geoschem/hemco for complete history
867 : !EOP
868 : !------------------------------------------------------------------------------
869 : !BOC
870 : !
871 : ! !LOCAL VARIABLES:
872 : !
873 : INTEGER :: I, J, L
874 : REAL(hp) :: TK_SURF
875 : REAL(hp) :: KAW_T, FLUX, KOL_T
876 : REAL(hp) :: DTSRCE
877 : REAL(hp) :: KA_H2O, KA_POP, KW_CO2, KW_POP
878 : REAL(hp) :: DA_H2O, DA_POP, DW_CO2, DW_POP
879 : REAL(hp) :: SCH_CO2, SCH_POP
880 : REAL(hp) :: TK, PRESS
881 : REAL(hp) :: C_DISS
882 : REAL(hp) :: POPG, U10M, ALPHA
883 : REAL(hp) :: FRAC_LAKE, VISC_H2O
884 : REAL(hp) :: SFCWINDSQR
885 : REAL(hp) :: AREA_M2
886 : LOGICAL :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE
887 : !
888 : ! !DEFINED PARAMETERS:
889 : !
890 : ! Delta H for POP:
891 : ! For PHENANTHRENE:
892 : ! this is the Delta H for phase transfer
893 : ! from air to water (Schwarzenbach,
894 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or 47 [kJ/mol].
895 : ! For PYRENE: 43000 [kJ/mol]
896 : REAL(hp) :: DEL_HW
897 :
898 : ! R = universal gas constant for adjusting KOA for temp:
899 : ! 8.3144598d-3 [kJ/mol/K]
900 : REAL(hp), PARAMETER :: R = 8.3144598d-3
901 :
902 : ! Molecular weight
903 : ! For phe, 0.17823 kg/mol
904 : REAL(hp) :: MWPOP
905 :
906 : ! Molecular weight of air
907 : REAL(hp), PARAMETER :: MWAIR = 28.97d0 ! g/mol
908 :
909 : ! Molecular weight of water
910 : REAL(hp), PARAMETER :: MWH2O = 18.1d0 ! g/mol
911 :
912 : ! Molar volumes calculated following Abraham and McGowan 1987 as
913 : ! summarized by Schwarzenbach et al. 2003.
914 : ! Each element is assigned a characteristic atomic volume, and an atomic
915 : ! volume of 6.56 cm3/mol is subtracted for each bond, no matter whether
916 : ! single, double, or triple
917 : ! C = 16.35, H = 8.71, O = 12.43, N = 14.39, P = 24.87, F = 10.48
918 : ! Br = 26.21, I = 34.53, S = 22.91, Si = 26.83
919 :
920 : ! Molar volume of water
921 : ! 2*(8.71) + 12.43 - 2*(6.56) = 16.73
922 : REAL(hp), PARAMETER :: V_H2O = 16.73d0 ! cm3/mol
923 :
924 : ! Molar volume of CO2
925 : ! 16.35 + 2*(12.43) - 2*(6.56) = 28.1
926 : REAL(hp), PARAMETER :: V_CO2 = 28.1d0 ! cm3/mol
927 :
928 : ! Molar volume of POP
929 : ! For PHE (C16H10):
930 : ! 16*(16.35) + 10*(8.71) - 29*(6.56) =
931 : REAL(hp), PARAMETER :: V_POP = 538.94d0 ! cm3/mol
932 :
933 : ! Molar volume of air - average of gases in air
934 : REAL(hp), PARAMETER :: V_AIR = 20.1d0 ! cm3/mol
935 :
936 : ! For PHENANTHRENE:
937 : ! log KAW_298 = -2.76, or 1.74*10-3 [unitless]
938 : ! For PYRENE:
939 : ! log KAW_298 = -3.27, or 5.37*10-4 [unitless]
940 : REAL(hp) :: KAW_298
941 :
942 : ! Set the kinematic viscosity of freshwater at 20C
943 : !REAL(hp), PARAMETER :: VISC_H2O ! = 1d0 ![cm2/s]
944 :
945 : ! Set aqueous degradation rate
946 : !REAL(hp), PARAMETER :: DEGR != 3.5d-5 ![/h]
947 :
948 : !=================================================================
949 : ! LAKEEMISPOP begins here!
950 : !=================================================================
951 :
952 : ! Copy values from ExtState
953 0 : DEL_HW = ExtState%POP_DEL_HW
954 0 : MWPOP = ExtState%POP_XMW
955 0 : KAW_298 = ExtState%POP_HSTAR
956 :
957 : ! Emission timestep [s]
958 0 : DTSRCE = HcoState%TS_EMIS
959 :
960 0 : DO J=1, HcoState%NY
961 0 : DO I=1, HcoState%NX
962 :
963 : ! Set logicals
964 : ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE)
965 : ! IS_LAND will return non-ocean boxes but may still contain lakes
966 : ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE)
967 : IS_LAND_OR_ICE = ( ( IS_LAND(I,J,ExtState) ) .OR. &
968 0 : ( IS_ICE (I,J,ExtState) ) )
969 : IS_SNOW_OR_ICE = ( ( IS_ICE (I,J,ExtState) ) .OR. &
970 0 : ( IS_LAND(I,J,ExtState) .AND. &
971 0 : ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) )
972 :
973 : ! Do soils routine only if we are on land that is not covered with
974 : ! snow or ice
975 0 : IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN
976 :
977 : ! Get fraction of grid box covered by lake surface area
978 0 : FRAC_LAKE = ExtState%FRLAKE%Arr%Val(I,J)
979 :
980 0 : IF ( FRAC_LAKE > 0e+0_hp ) THEN
981 :
982 : ! Get surface skin temp [K]
983 0 : TK_SURF = ExtState%TSKIN%Arr%Val(I,J)
984 :
985 : ! Get air temp [K]
986 0 : TK = ExtState%TK%Arr%Val(I,J,1)
987 :
988 : ! Get surface pressure at end of dynamic time step [hPa]
989 0 : PRESS = ExtState%PSC2_WET%Arr%Val(I,J)
990 :
991 : ! Convert to units of atm
992 0 : PRESS = PRESS / 1013.25e+0_hp
993 :
994 : ! Get gas phase air POP concentration at surface in mol/m3
995 : ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15)
996 0 : POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM )
997 :
998 : ! (kg trc/kg dry air) / (kg trc/mol) * (kg dry air/m3)
999 0 : POPG = POPG / MWPOP * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3
1000 :
1001 : ! Grid box surface area [m2]
1002 0 : AREA_M2 = HcoState%Grid%AREA_M2%Val(I,J)
1003 :
1004 : ! Get the dissolved POP concentration at lake surface in mol/m3
1005 : ! Distribute the total deposited mass to a volume that best
1006 : ! matches observed dissolved concentrations
1007 : ! Future versions should consider aqueous particle concentrations
1008 : ! and sinking rates and photolytic/microbial degradation
1009 :
1010 : ! Start with 1 m - scale by 100
1011 0 : C_DISS = POP_SURF(I,J) / MWPOP / AREA_M2 / 100e+0_hp
1012 :
1013 : ! Wind speed at 10m altitude [m/s]
1014 0 : SFCWINDSQR = ExtState%U10M%Arr%Val(I,J)**2 &
1015 0 : + ExtState%V10M%Arr%Val(I,J)**2
1016 0 : U10M = SQRT( SFCWINDSQR )
1017 :
1018 : ! Need to calculate water-side and air-side mass transfer
1019 : ! coefficients
1020 : ! Start with air-side
1021 : ! Relate air-side MTC of POP to that of H2O
1022 : ! First, calculate air-side MTC for water following
1023 : ! Schwarzenbach, Gschwend, Imboden 2003
1024 0 : KA_H2O = 0.2e+0_hp * U10M + 0.3e+0_hp ! cm/s
1025 :
1026 : ! Relate air-side MTC for water to that of POP via diffusivities
1027 : ! following Schwarzenbach et al 2003
1028 : ! Calculate temperature-dependent diffusivities in air first
1029 : ! folling Fuller et al. 1966 (summarized by Schwarzenbach
1030 : ! et al. 2003)
1031 : DA_POP = 1e-3_hp * TK**1.75e+0_hp * ( (1e+0_hp/MWAIR) + &
1032 : (1e+0_hp / (MWPOP*1e+3_hp) ) )**0.5e+0_hp / &
1033 : ( PRESS * ( V_AIR**(1e+0_hp/3e+0_hp) + &
1034 0 : V_POP**(1e+0_hp/3e+0_hp))**2e+0_hp ) ! cm2/s
1035 :
1036 : DA_H2O = 1e-3_hp * TK**1.75e+0_hp * ( (1e+0_hp/MWAIR) + &
1037 : (1e+0_hp / MWH2O ) )**0.5e+0_hp / &
1038 : ( PRESS * ( V_AIR**(1e+0_hp/3e+0_hp) + &
1039 0 : V_H2O**(1e+0_hp/3e+0_hp))**2e+0_hp ) ! cm2/s
1040 :
1041 : ! Relate POP and H2O air-side MTCs
1042 0 : KA_POP = KA_H2O * ( DA_POP / DA_H2O )**( 0.67e+0_hp ) ! cm/s
1043 :
1044 : ! Now calculate water-side MTCs
1045 : ! Start with calculating the water side MTC of CO2
1046 : ! This depends on wind speed (Schwarzenbach et al 2003)
1047 : ! Three different scenarios for the water surface under
1048 : ! different wind speeds are considered:
1049 : ! Smooth Surface Regime (SSR): u10 <= 4.2 m/s
1050 : ! Rough Surface Regime (RSR): 4.2 m/s < u10 <= 13 m/s
1051 : ! Breaking Wave Regime (BWR): u10 > 13 m/s
1052 : ! Alpha, the exponent in the relationship for the water side MTC,
1053 : ! also depends on the wind speed. Set this as well
1054 0 : IF ( U10M <= 4.2e+0_hp ) THEN
1055 : KW_CO2 = 0.65e-3_hp ! cm/s
1056 : ALPHA = 0.67e+0_hp
1057 0 : ELSE IF ( U10M > 4.2e+0_hp .AND. U10M <= 13e+0_hp ) THEN
1058 0 : KW_CO2 = ( 0.79e+0_hp * U10M - 2.68e+0_hp ) * 1e-3_hp
1059 0 : ALPHA = 0.5e+0_hp
1060 0 : ELSE IF (U10M > 13e+0_hp) THEN
1061 0 : KW_CO2 = ( 1.64e+0_hp * U10M - 13.69e+0_hp ) * 1e-3_hp
1062 0 : ALPHA = 0.50e+0_hp
1063 : ENDIF
1064 :
1065 : ! Get the temperature-dependent kinematic viscosity of water
1066 0 : IF (TK_SURF <= 273.15 ) THEN
1067 : VISC_H2O = 1.787e-2_hp ! [cm2/s]
1068 0 : ELSE IF (TK_SURF > 273.15 .AND. TK_SURF <= 278.15 ) THEN
1069 : VISC_H2O = 1.518e-2_hp
1070 0 : ELSE IF (TK_SURF > 258.15 .AND. TK_SURF <= 283.15 ) THEN
1071 : VISC_H2O = 1.307e-2_hp
1072 0 : ELSE IF (TK_SURF > 283.15 .AND. TK_SURF <= 287.15 ) THEN
1073 : VISC_H2O = 1.139e-2_hp
1074 0 : ELSE IF (TK_SURF > 287.15 .AND. TK_SURF <= 293.15 ) THEN
1075 : VISC_H2O = 1.002e-2_hp
1076 0 : ELSE IF (TK_SURF > 293.15 .AND. TK_SURF <= 298.15 ) THEN
1077 : VISC_H2O = 0.89e-2_hp
1078 0 : ELSE IF (TK_SURF > 298.15 ) THEN
1079 0 : VISC_H2O = 0.797e-2_hp
1080 : ENDIF
1081 :
1082 : ! Calculate the diffusivites of CO2 and POP in water
1083 : DW_CO2 = ( 13.26 * 1e-5_hp ) / &
1084 : (( VISC_H2O*1e+2_hp )**1.14e+0_hp * &
1085 0 : (V_CO2)**0.589e+0_hp ) ! [cm2/s]
1086 :
1087 : DW_POP = ( 13.26 * 1e-5_hp ) / &
1088 : (( VISC_H2O*1e+2_hp )**1.14e+0_hp * &
1089 0 : (V_POP)**0.589e+0_hp ) ! [cm2/s]
1090 :
1091 : ! Calculate the Schmidt numbers for CO2 and POP
1092 0 : SCH_CO2 = VISC_H2O / DW_CO2 ! [unitless]
1093 0 : SCH_POP = VISC_H2O / DW_POP ! [unitless]
1094 :
1095 : ! Calculate the water-side MTC for POP
1096 0 : KW_POP = KW_CO2 * ( SCH_POP / SCH_CO2 ) ** (-ALPHA) ! [cm/s]
1097 :
1098 : ! Calculate the temperature-dependent dimensionless Henry's Law
1099 : ! constant
1100 : KAW_T = KAW_298 * EXP((-DEL_HW/R) * ((1e+0_hp/TK_SURF) - &
1101 0 : (1e+0_hp/298e+0_hp))) ! [unitless]
1102 :
1103 : ! Now calculate the overall air-water MTC
1104 : KOL_T = 1e+0_hp / ( 1e+0_hp/KW_POP + &
1105 0 : 1e+0_hp / (KA_POP*KAW_T) ) ! [cm/s]
1106 :
1107 : ! Calculate Flux in ng/m2/day !
1108 : FLUX = KOL_T * 3600e+0_hp * 24e+0_hp * ( C_DISS - &
1109 0 : POPG/KAW_T ) * MWPOP * 1e+12_hp / 100e+0_hp
1110 :
1111 : ! Convert to an emission rate in kg/m2/s for returning to
1112 : ! HcoX_GC_POPs_Run
1113 : ! Only return it if it's positive
1114 0 : Inst%EmisPOPGfromLake(I,J) = MAX(FLUX / 24e+0_hp / 3600e+0_hp / &
1115 0 : 1e+12_hp, 0e+0_hp )
1116 :
1117 : ! Multiply the emissions by the fraction of land that is water
1118 0 : Inst%EmisPOPGfromLake(I,J) = Inst%EmisPOPGfromLake(I,J) * &
1119 0 : FRAC_LAKE
1120 :
1121 : ! If the flux is positive, then the direction will be from the
1122 : ! soil to the air.
1123 : ! If the flux is zero or negative, store it in a separate array.
1124 0 : IF ( FLUX > 0e+0_hp ) THEN
1125 :
1126 : ! Store positive flux
1127 0 : Inst%FluxPOPGfromLakeToAir(I,J) = + FLUX
1128 :
1129 : ! Make sure negative flux diagnostic has nothing added to it
1130 0 : Inst%FluxPOPGfromAirToLake(I,J) = 0e+0_hp
1131 :
1132 : ! Store the soil/air fugacity ratio
1133 0 : Inst%FugacityLakeToAir(I,J) = C_DISS / (POPG/KAW_T)
1134 :
1135 0 : ELSE IF ( FLUX <= 0e+0_hp ) THEN
1136 :
1137 : ! Store the negative flux
1138 0 : Inst%FluxPOPGfromAirToLake(I,J) = FLUX
1139 :
1140 : ! Add nothing to positive flux
1141 0 : Inst%FluxPOPGfromLakeToAir(I,J) = 0e+0_hp
1142 :
1143 : ! Continue to store the fugacity ratio
1144 0 : Inst%FugacityLakeToAir(I,J) = C_DISS / (POPG/KAW_T)
1145 :
1146 : ENDIF
1147 :
1148 : ENDIF
1149 :
1150 : ELSE
1151 :
1152 : ! We are not on land or the land is covered with ice or snow
1153 : ! or we are land but there is no water
1154 0 : FLUX = 0e+0_hp
1155 0 : Inst%EmisPOPGfromLake = 0e+0_hp
1156 0 : Inst%FluxPOPGfromLakeToAir(I,J) = 0e+0_hp
1157 0 : Inst%FluxPOPGfromAirToLake(I,J) = 0e+0_hp
1158 0 : Inst%FugacityLakeToAir(I,J) = 0e+0_hp
1159 :
1160 : ENDIF
1161 :
1162 : ENDDO
1163 : ENDDO
1164 :
1165 0 : END SUBROUTINE LAKEEMISPOP
1166 : !EOC
1167 : !------------------------------------------------------------------------------
1168 : ! GEOS-Chem Global Chemical Transport Model !
1169 : !------------------------------------------------------------------------------
1170 : !BOP
1171 : !
1172 : ! !IROUTINE: vegemispop
1173 : !
1174 : ! !DESCRIPTION: Subroutine VEGEMISPOP is the subroutine for secondary
1175 : ! POP emissions from leaves.
1176 : !\\
1177 : !\\
1178 : ! !INTERFACE:
1179 : !
1180 0 : SUBROUTINE VEGEMISPOP( POP_SURF, HcoState, ExtState, Inst )
1181 : !
1182 : ! !INPUT PARAMETERS:
1183 : !
1184 : REAL(sp), DIMENSION(:,:), INTENT(IN) :: POP_SURF ! POP sfc conc [kg]
1185 : TYPE(HCO_State), POINTER :: HcoState ! Hemco state
1186 : TYPE(Ext_State), POINTER :: ExtState ! Module options
1187 : TYPE(MyInst), POINTER :: Inst ! Instance
1188 : !
1189 : ! !REVISION HISTORY:
1190 : ! 21 Aug 2012 - C.L. Friedman - Initial version based on LAND_MERCURY_MOD
1191 : ! See https://github.com/geoschem/hemco for complete history
1192 : !EOP
1193 : !------------------------------------------------------------------------------
1194 : !BOC
1195 : !
1196 : ! !LOCAL VARIABLES:
1197 : !
1198 : INTEGER :: I, J, L
1199 : REAL(hp) :: POPG, POPG_GL
1200 : REAL(hp) :: FRAC_SNOWFREE_LAND, TK_SURF
1201 : REAL(hp) :: KLA_T, FLUX, K_MT
1202 : REAL(hp) :: LEAF_CONC, DTSRCE, KOA_T
1203 : REAL(hp) :: DIFF, FLEAF, FAIR, DS
1204 : REAL(hp) :: DAB_F, DC, DAL, PC, UC
1205 : REAL(hp) :: ZLEAF, ZAIR, TK
1206 : REAL(hp) :: NEWLEAF
1207 : REAL(hp) :: LAI, KAW_T, KOW_T
1208 : REAL(hp) :: FUG_R, NEW_LEAF, DLA
1209 : REAL(hp) :: AREA_M2
1210 : LOGICAL :: IS_SNOW_OR_ICE, IS_LAND_OR_ICE
1211 : LOGICAL :: IS_SNOWFREE_LAND
1212 :
1213 : !
1214 : ! !DEFINED PARAMETERS:
1215 : !
1216 : ! Delta H for POP [kJ/mol]. Delta H is enthalpy of phase transfer
1217 : ! from gas phase to OC. For now we use Delta H for phase transfer
1218 : ! from the gas phase to the pure liquid state.
1219 : ! For PHENANTHRENE:
1220 : ! this is taken as the negative of the Delta H for phase transfer
1221 : ! from the pure liquid state to the gas phase (Schwarzenbach,
1222 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -74000 [J/mol].
1223 : ! For PYRENE:
1224 : ! this is taken as the negative of the Delta H for phase transfer
1225 : ! from the pure liquid state to the gas phase (Schwarzenbach,
1226 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or -87000 [J/mol].
1227 : ! For BENZO[a]PYRENE:
1228 : ! this is also taken as the negative of the Delta H for phase transfer
1229 : ! from the pure liquid state to the gas phase (Schwarzenbach,
1230 : ! Gschwend, Imboden, 2003, pg 452, Prob 11.1), or -110,000 [J/mol]
1231 : REAL(hp) :: DEL_H
1232 :
1233 : ! Delta H for POP [kJ/mol].
1234 : ! For PHENANTHRENE:
1235 : ! this is the Delta H for phase transfer
1236 : ! from air to water (Schwarzenbach,
1237 : ! Gschwend, Imboden, 2003, pg 200, Table 6.3), or 47000 [J/mol].
1238 : ! For PYRENE: 43000 [J/mol]
1239 : REAL(hp) :: DEL_HW
1240 :
1241 : ! R = universal gas constant for adjusting KOA for temp:
1242 : ! 8.314 [J/mol/K OR m3*Pa/K/mol]
1243 : REAL(hp), PARAMETER :: R = 8.3144598e+0_hp
1244 :
1245 : ! Molecular weight
1246 : ! For phe, 0.17823 kg/mol
1247 : REAL(hp) :: MW
1248 :
1249 : ! Molecular weight of air
1250 : REAL(hp), PARAMETER :: MWAIR = 28.97d0 ! g/mol
1251 :
1252 : ! For PHENANTHRENE:
1253 : ! log KOA_298 = 7.64, or 4.37*10^7 [unitless]
1254 : ! For PYRENE:
1255 : ! log KOA_298 = 8.86, or 7.24*10^8 [unitless]
1256 : ! For BENZO[a]PYRENE:
1257 : ! log KOA_298 = 11.48, or 3.02*10^11 [unitless]
1258 : ! (Ma et al., J. Chem. Eng. Data, 2010, 55:819-825).
1259 : REAL(hp) :: KOA_298
1260 :
1261 : ! For PHENANTHRENE:
1262 : ! log KAW_298 = -2.76, or 1.74*10-3 [unitless]
1263 : ! For PYRENE:
1264 : ! log KAW_298 = -3.27, or 5.37*10-4 [unitless]
1265 : REAL(hp) :: KAW_298
1266 :
1267 : ! Set volume fractions of octanol and water in surface and reservoir
1268 : ! leaf compartments [unitless]
1269 : REAL(hp), PARAMETER :: OCT_SURF = 0.8e+0_hp
1270 : REAL(hp), PARAMETER :: OCT_RES = 0.02e+0_hp
1271 : REAL(hp), PARAMETER :: H2O_RES = 0.7e+0_hp
1272 :
1273 : ! Set thickness of different leaf compartments. Volumes calculated by
1274 : ! multiplying thicknesses by leaf area index
1275 : REAL(hp), PARAMETER :: SURF_THICK = 2e-6_hp ! m
1276 : REAL(hp), PARAMETER :: RES_THICK = 250e-6_hp ! m
1277 :
1278 : ! Set transfer velocity and diffusion coefficient values
1279 : REAL(hp), PARAMETER :: UAB_F = 9e+0_hp ![m/h]
1280 :
1281 : ! Set soil degradation rate
1282 : REAL(hp), PARAMETER :: DEGR = 3.5e-5_hp ![/h]
1283 :
1284 : !=================================================================
1285 : ! VEGEMISPOP begins here!
1286 : !=================================================================
1287 :
1288 : ! Copy values from ExtState
1289 0 : DEL_H = ExtState%POP_DEL_H
1290 0 : KOA_298 = ExtState%POP_KOA
1291 0 : DEL_HW = ExtState%POP_DEL_Hw
1292 0 : KAW_298 = ExtState%POP_HSTAR
1293 0 : MW = ExtState%POP_XMW
1294 :
1295 : ! Emission timestep [s]
1296 0 : DTSRCE = HcoState%TS_EMIS
1297 :
1298 0 : DO J=1, HcoState%NY
1299 0 : DO I=1, HcoState%NX
1300 :
1301 : ! Set logicals
1302 : ! Is grid box covered by land/ice or by water? (IS_LAND_OR_ICE)
1303 : ! IS_LAND will return non-ocean boxes but may still contain lakes
1304 : ! If land, is it covered by snow/ice? (IS_SNOW_OR_ICE)
1305 : IS_LAND_OR_ICE = ( (IS_LAND(I,J,ExtState)) .OR. &
1306 0 : (IS_ICE (I,J,ExtState)) )
1307 : IS_SNOW_OR_ICE = ( (IS_ICE (I,J,ExtState)) .OR. &
1308 0 : (IS_LAND(I,J,ExtState) .AND. &
1309 0 : ExtState%SNOWHGT%Arr%Val(I,J) > 10e+0_hp ) )
1310 :
1311 : ! Do soils routine only if we are on land that is not covered with
1312 : ! snow or ice
1313 0 : IF ((IS_LAND_OR_ICE) .AND. .NOT. ( IS_SNOW_OR_ICE )) THEN
1314 :
1315 : ! Get fraction of grid box covered by leaf surface area
1316 : ! Do not consider different vegetation types for now
1317 0 : LAI = ExtState%LAI%Arr%Val(I,J)
1318 :
1319 0 : IF ( LAI > 0e+0_hp ) THEN
1320 :
1321 : ! Get surface skin temp [K]
1322 0 : TK_SURF = ExtState%TSKIN%Arr%Val(I,J)
1323 :
1324 : ! Get air temp [K]
1325 0 : TK = ExtState%TK%Arr%Val(I,J,1)
1326 :
1327 : ! Get gas phase air POP concentration at surface in mol/m3
1328 : ! ExtState%POPG is now in units of kg/kg dry air (ewl, 10/2/15)
1329 0 : POPG = MAX( ExtState%POPG%Arr%Val(I,J,1), SMALLNUM )
1330 :
1331 : ! (kg trc/kg dry air) / (kg trc/mol) * (kg dry air/m3)
1332 0 : POPG = POPG / MW * ExtState%AIRDEN%Arr%Val(I,J,1) ! mol/m3
1333 :
1334 : ! Grid box surface area [m2]
1335 0 : AREA_M2 = HcoState%Grid%AREA_M2%Val(I,J)
1336 :
1337 : ! Only consider partitioning into leaf surface (not reservoir)
1338 : ! for now following Mackay et al 2006 Environ Sci & Pollut Res
1339 : ! Include reservoir when land-atm models become dynamic
1340 :
1341 : ! Assume that all leaf surfaces contain an average lipid content
1342 : ! of 80% (Mackay et al 2006)
1343 :
1344 : ! Get leaf concentration
1345 : ! Convert to mol/m3
1346 : ! kg deposited to leaf in 1 yr * / 0.178 kg/mol
1347 : ! / area grid box (m2) / surface thickness m
1348 0 : LEAF_CONC = POP_SURF(I,J) / MW / AREA_M2 / SURF_THICK ! mol/m3
1349 :
1350 : ! Check concentration in leaves by assuming a density similar
1351 : ! to water (1 g/cm3)
1352 :
1353 : ! No degradation/metabolism for now. Just scale leaf
1354 : ! concentrations to match flux observations
1355 0 : NEWLEAF = LEAF_CONC/1e+4_hp !SCALING FACTOR
1356 :
1357 : ! Define temperature-dependent KOA:
1358 : KOA_T = KOA_298 * EXP((-DEL_H/R) * ((1e+0_hp/TK_SURF) - &
1359 0 : (1e+0_hp/298e+0_hp)))
1360 :
1361 : ! Calculate the temperature-dependent dimensionless Henry's Law
1362 : ! constant
1363 : KAW_T = KAW_298 * EXP((-DEL_HW/R) * ((1e+0_hp/TK_SURF) - &
1364 0 : (1e+0_hp/298e+0_hp))) ! [unitless]
1365 :
1366 : ! Estimate the temperature-dependent dimensionless octanol-water
1367 : ! constant
1368 0 : KOW_T = KOA_T * KAW_T
1369 :
1370 : ! Define dimensionless leaf surface-air partition coefficient
1371 : ! (mol/m3 leaf / mol/m3 air)
1372 : ! KLA = foct_surf * Koa
1373 0 : KLA_T = OCT_SURF * KOA_T
1374 : ! KLA_T = MAX( KLA_T, SMALLNUM )
1375 :
1376 : ! Calculate fugacities from concentrations by dividing by "Z"
1377 : ! values, or the fugacity capacity in mol/m3*Pa following Mackay
1378 : ! and Paterson, 1991
1379 :
1380 : ! fleaf = Cleaf * R * T / KSA [Pa]
1381 : ! where Cleaf is in mol/m3, R is in Pa * m3 / mol K
1382 : ! T is in K and KLA is dimensionless
1383 0 : FLEAF = NEWLEAF * R * TK_SURF / KLA_T
1384 :
1385 : ! fair = Cair * R * T [Pa]
1386 : ! where Cair is in mol/m3, R is in Pa * m3 / mol K and T is in K
1387 0 : FAIR = POPG * R * TK
1388 :
1389 : ! Calculate the fugacity gradient [Pa]
1390 : ! If the gradient is negative, fair is larger and the POP will
1391 : ! diffuse from air to soil
1392 : ! If the gradient is positive, fsoil is larger and POP will
1393 : ! diffuse from soil to air
1394 0 : DIFF = FLEAF - FAIR
1395 0 : FUG_R = FLEAF/FAIR
1396 :
1397 : ! Calculate "Z" values from fugacities.
1398 : ! Z is the fugacity capacity in mol/m3*Pa. C = Z*f, so Z = C/f
1399 0 : ZAIR = POPG / FAIR ! (mol/m3) / (Pa)
1400 0 : ZLEAF = NEWLEAF / FLEAF ! (mol/m3) / (Pa)
1401 :
1402 : ! Calculate the "D" value, or the transfer coefficient that
1403 : ! describes the movement of POP between phases (Mackay and
1404 : ! Paterson, 1991, Cousins and Mackay 2000, internal report).
1405 : ! [mol/h*Pa]
1406 : ! The D value for leaf surface-air gas diffusion is given by
1407 : ! Dla = 1 / (1/Dc + 1/(Dab-f)) [mol/(Pa*h)]
1408 : ! where Dab-f is the boundary layer diffusion [mol/h*Pa]
1409 : ! given by Dab-f = As * L * Uab-f * Za
1410 : ! where As is the area of the land surface [m2], L is the leaf
1411 : ! area index [m2/m2],
1412 : ! Uab-f is a mass transfer coefficient for surface-air boundary
1413 : ! layer diffusion [m/h],
1414 : ! and Za is the fugacity capacity of the air [mol/(m3*Pa)]
1415 : ! Dc is the cuticle diffusion, given by
1416 : ! Dc = As * L * Uc * Zf
1417 : ! where As and L are as above, Uc is the cuticle mass transfer
1418 : ! coefficient [m/h],
1419 : ! and Zf is the fugacity capacity of the leaf surface
1420 : ! (mol/(m3*Pa))
1421 :
1422 : ! Uc is given by
1423 : ! Uc = 3600 * Pc * 1/Kaw
1424 : ! where Pc is the cuticle permeance (m/s) and Kaw is the
1425 : ! dimensionless air-water partition coefficient.
1426 : ! Pc is given by
1427 : ! Log Pc = ((0.704 * log Kow - 11.2) +
1428 : ! (-3.47 - 2.79 * logMW + 0.970 log Kow)) / 2
1429 : ! (an average of two equations)
1430 :
1431 : ! Need to define each D value
1432 : ! DAB_F:
1433 : ! m/h * mol/(m3*Pa) = (mol/h*Pa*m2)
1434 0 : DAB_F = UAB_F * ZAIR ! mol/(h*Pa*m2)
1435 :
1436 : ! Calculate PC and then Uc in order to calculate Dc
1437 : ! PC, UC are calculated according to Cousins and Mackay,
1438 : ! Chemosphere, 2001, Table 2
1439 : PC = 10** (( 0.704e+0_hp * LOG (KOW_T) - 11.2e+0_hp ) + &
1440 : ( -3.47e+0_hp -2.79e+0_hp* LOG(MW*1000d0) + 0.97e+0_hp &
1441 0 : * LOG(KOW_T))/2e+0_hp) ![m/s]
1442 :
1443 0 : UC = 3600e+0_hp * PC * 1e+0_hp/KAW_T ! [m/h]
1444 :
1445 0 : DC = UC * ZLEAF ! mol/(h*Pa*m2)
1446 :
1447 : ! Now calculate overall transfer ! mol/(h*Pa*m2)
1448 0 : DLA = 1e+0_hp / (1e+0_hp/DC + 1e+0_hp/DAB_F)
1449 :
1450 : ! Calculate Flux in mol/h/m2
1451 0 : FLUX = DLA * DIFF
1452 :
1453 : ! Change to units of ng/m2/d for storage
1454 0 : FLUX = FLUX * 24e+0_hp * MW * 1e+12_hp
1455 :
1456 : ! Convert to an emission rate in kg/m2/s for returning to
1457 : ! HcoX_GC_POPs_Run
1458 : ! Only want to add rates that are positive
1459 0 : Inst%EmisPOPGfromLeaf(I,J) = MAX(FLUX * LAI / 24e+0_hp / &
1460 0 : 3600e+0_hp / 1e+12_hp, 0e+0_hp)
1461 :
1462 : ! If the flux is positive, then the direction will be from the
1463 : ! soil to the air.
1464 : ! If the flux is zero or negative, store it in a separate array.
1465 0 : IF ( FLUX > 0e+0_hp ) THEN
1466 :
1467 : ! Store positive flux
1468 0 : Inst%FluxPOPGfromLeafToAir(I,J) = FLUX
1469 :
1470 : ! Make sure negative flux diagnostic has nothing added to it
1471 0 : Inst%FluxPOPGfromAirtoLeaf = 0e+0_hp
1472 :
1473 : ! Store the soil/air fugacity ratio
1474 0 : Inst%FugacityLeafToAir(I,J) = FUG_R
1475 :
1476 0 : ELSE IF ( FLUX <= 0e+0_hp ) THEN
1477 :
1478 : ! Store the negative flux
1479 0 : Inst%FluxPOPGfromAirtoLeaf(I,J) = FLUX
1480 :
1481 : ! Add nothing to positive flux
1482 0 : Inst%FluxPOPGfromLeafToAir(I,J) = 0e+0_hp
1483 :
1484 : ! Continue to store the fugacity ratio
1485 0 : Inst%FugacityLeafToAir(I,J) = FUG_R
1486 :
1487 : ENDIF
1488 :
1489 : ELSE
1490 :
1491 : ! We are not on land or the land is covered with ice or snow
1492 0 : FLUX = 0e+0_hp
1493 0 : Inst%EmisPOPGfromLeaf(I,J) = 0e+0_hp
1494 0 : Inst%FluxPOPGfromLeafToAir(I,J) = 0e+0_hp
1495 0 : Inst%FluxPOPGfromAirtoLeaf(I,J) = 0e+0_hp
1496 0 : Inst%FugacityLeafToAir(I,J) = 0e+0_hp
1497 :
1498 : ENDIF
1499 :
1500 : ENDIF
1501 :
1502 : ENDDO
1503 : ENDDO
1504 :
1505 0 : END SUBROUTINE VEGEMISPOP
1506 : !EOC
1507 : !------------------------------------------------------------------------------
1508 : ! GEOS-Chem Global Chemical Transport Model !
1509 : !------------------------------------------------------------------------------
1510 : !BOP
1511 : !
1512 : ! !IROUTINE: is_land
1513 : !
1514 : ! !DESCRIPTION: Function IS\_LAND returns TRUE if surface grid box (I,J) is
1515 : ! a land box.
1516 : !\\
1517 : !\\
1518 : ! !INTERFACE:
1519 : !
1520 0 : FUNCTION IS_LAND( I, J, ExtState ) RESULT ( LAND )
1521 : !
1522 : ! !USES:
1523 : !
1524 : USE HCO_GeoTools_Mod, ONLY : HCO_LANDTYPE
1525 : !
1526 : ! !INPUT PARAMETERS:
1527 : !
1528 : INTEGER, INTENT(IN) :: I ! Longitude index of grid box
1529 : INTEGER, INTENT(IN) :: J ! Latitude index of grid box
1530 : TYPE(Ext_State), POINTER :: ExtState ! Module options
1531 : !
1532 : ! !RETURN VALUE:
1533 : !
1534 : LOGICAL :: LAND ! =T if it is a land box
1535 : !
1536 : ! !REVISION HISTORY:
1537 : ! 26 Jun 2000 - R. Yantosca - Initial version
1538 : ! See https://github.com/geoschem/hemco for complete history
1539 : !EOP
1540 : !------------------------------------------------------------------------------
1541 : !BOC
1542 0 : LAND = HCO_LANDTYPE( ExtState%FRLAND%Arr%Val(I,J), &
1543 0 : ExtState%FRLANDIC%Arr%Val(I,J), &
1544 0 : ExtState%FROCEAN%Arr%Val(I,J), &
1545 0 : ExtState%FRSEAICE%Arr%Val(I,J), &
1546 0 : ExtState%FRLAKE%Arr%Val(I,J) ) == 1
1547 :
1548 0 : END FUNCTION IS_LAND
1549 : !EOC
1550 : !------------------------------------------------------------------------------
1551 : ! GEOS-Chem Global Chemical Transport Model !
1552 : !------------------------------------------------------------------------------
1553 : !BOP
1554 : !
1555 : ! !IROUTINE: is_ice
1556 : !
1557 : ! !DESCRIPTION: Function IS\_ICE returns TRUE if surface grid box (I,J)
1558 : ! contains either land-ice or sea-ice.
1559 : !\\
1560 : !\\
1561 : ! !INTERFACE:
1562 : !
1563 0 : FUNCTION IS_ICE( I, J, ExtState ) RESULT ( ICE )
1564 : !
1565 : ! !USES:
1566 : !
1567 : USE HCO_GeoTools_Mod, ONLY : HCO_LANDTYPE
1568 : !
1569 : ! !INPUT PARAMETERS:
1570 : !
1571 : INTEGER, INTENT(IN) :: I ! Longitude index of grid box
1572 : INTEGER, INTENT(IN) :: J ! Latitude index of grid box
1573 : TYPE(Ext_State), POINTER :: ExtState ! Module options
1574 : !
1575 : ! !RETURN VALUE:
1576 : !
1577 : LOGICAL :: ICE ! =T if this is an ice box
1578 : !
1579 : ! !REVISION HISTORY:
1580 : ! 09 Aug 2005 - R. Yantosca - Initial version
1581 : ! See https://github.com/geoschem/hemco for complete history
1582 : !EOP
1583 : !------------------------------------------------------------------------------
1584 : !BOC
1585 :
1586 0 : ICE = HCO_LANDTYPE( ExtState%FRLAND%Arr%Val(I,J), &
1587 0 : ExtState%FRLANDIC%Arr%Val(I,J), &
1588 0 : ExtState%FROCEAN%Arr%Val(I,J), &
1589 0 : ExtState%FRSEAICE%Arr%Val(I,J), &
1590 0 : ExtState%FRLAKE%Arr%Val(I,J) ) == 2
1591 :
1592 0 : END FUNCTION IS_ICE
1593 : !EOC
1594 : !------------------------------------------------------------------------------
1595 : ! Harmonized Emissions Component (HEMCO) !
1596 : !------------------------------------------------------------------------------
1597 : !BOP
1598 : !
1599 : ! !IROUTINE: HCOX_GC_POPs_Init
1600 : !
1601 : ! !DESCRIPTION: Subroutine HcoX\_GC\_POPs\_Init initializes the HEMCO
1602 : ! GC\_POPs extension.
1603 : !\\
1604 : !\\
1605 : ! !INTERFACE:
1606 : !
1607 0 : SUBROUTINE HCOX_GC_POPs_Init( HcoState, ExtName, ExtState, RC )
1608 : !
1609 : ! !USES:
1610 : !
1611 : USE HCO_ExtList_Mod, ONLY : GetExtNr
1612 : USE HCO_STATE_MOD, ONLY : HCO_GetExtHcoID
1613 : !
1614 : ! !INPUT PARAMETERS:
1615 : !
1616 : CHARACTER(LEN=*), INTENT(IN ) :: ExtName ! Extension name
1617 : TYPE(Ext_State), POINTER :: ExtState ! Module options
1618 : TYPE(HCO_State), POINTER :: HcoState ! Hemco state
1619 : !
1620 : ! !INPUT/OUTPUT PARAMETERS:
1621 : !
1622 : INTEGER, INTENT(INOUT) :: RC
1623 :
1624 : ! !REVISION HISTORY:
1625 : ! 19 Aug 2014 - M. Sulprizio- Initial version
1626 : ! See https://github.com/geoschem/hemco for complete history
1627 : !EOP
1628 : !------------------------------------------------------------------------------
1629 : !BOC
1630 : !
1631 : ! !LOCAL VARIABLES:
1632 : !
1633 : ! Scalars
1634 : INTEGER :: I, N, nSpc, ExtNr
1635 : CHARACTER(LEN=255) :: MSG, LOC
1636 : CHARACTER(LEN=255) :: DiagnName
1637 : CHARACTER(LEN=255) :: OutUnit
1638 :
1639 : ! Arrays
1640 0 : INTEGER, ALLOCATABLE :: HcoIDs(:)
1641 0 : CHARACTER(LEN=31), ALLOCATABLE :: SpcNames(:)
1642 :
1643 : ! Pointers
1644 : TYPE(MyInst), POINTER :: Inst
1645 0 : REAL(sp), POINTER :: Target2D(:,:)
1646 :
1647 : !=======================================================================
1648 : ! HCOX_GC_POPs_INIT begins here!
1649 : !=======================================================================
1650 0 : LOC = 'HCOX_GC_POPs_INIT (HCOX_GC_POPS_MOD.F90)'
1651 :
1652 : ! Get the extension number
1653 0 : ExtNr = GetExtNr( HcoState%Config%ExtList, TRIM(ExtName) )
1654 0 : IF ( ExtNr <= 0 ) RETURN
1655 :
1656 : ! Enter HEMCO
1657 0 : CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
1658 0 : IF ( RC /= HCO_SUCCESS ) THEN
1659 0 : CALL HCO_ERROR( 'ERROR 6', RC, THISLOC=LOC )
1660 0 : RETURN
1661 : ENDIF
1662 :
1663 : ! Create Instance
1664 0 : Inst => NULL()
1665 0 : CALL InstCreate ( ExtNr, ExtState%GC_POPs, Inst, RC )
1666 0 : IF ( RC /= HCO_SUCCESS ) THEN
1667 0 : CALL HCO_ERROR ( 'Cannot create GC_POPs instance', RC )
1668 0 : RETURN
1669 : ENDIF
1670 : ! Also fill ExtNrSS - this is the same as the parent ExtNr
1671 0 : Inst%ExtNr = ExtNr
1672 :
1673 : ! Set species IDs
1674 0 : CALL HCO_GetExtHcoID( HcoState, Inst%ExtNr, HcoIDs, SpcNames, nSpc, RC )
1675 0 : IF ( RC /= HCO_SUCCESS ) THEN
1676 0 : CALL HCO_ERROR( 'ERROR 7', RC, THISLOC=LOC )
1677 0 : RETURN
1678 : ENDIF
1679 :
1680 : ! Verbose mode
1681 0 : IF ( HcoState%amIRoot ) THEN
1682 :
1683 : ! Write the name of the extension regardless of the verbose setting
1684 0 : msg = 'Using HEMCO extension: GC_POPs (POPs emissions)'
1685 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
1686 0 : CALL HCO_Msg( HcoState%Config%Err, msg, sep1='-' ) ! with separator
1687 : ELSE
1688 0 : CALL HCO_Msg( msg, verb=.TRUE. ) ! w/o separator
1689 : ENDIF
1690 :
1691 : ! Write all other messages as debug printout only
1692 0 : MSG = 'Use the following species (Name: HcoID):'
1693 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1694 0 : DO N = 1, nSpc
1695 0 : WRITE(MSG,*) TRIM(SpcNames(N)), ':', HcoIDs(N)
1696 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
1697 : ENDDO
1698 : ENDIF
1699 :
1700 : ! Set up tracer indices
1701 0 : DO N = 1, nSpc
1702 0 : SELECT CASE( TRIM( SpcNames(N) ) )
1703 : CASE( 'POPG', 'POPG_BaP', 'POPG_PHE', 'POPG_PYR' )
1704 0 : Inst%IDTPOPG = HcoIDs(N)
1705 : CASE( 'POPPOCPO', 'POPPOCPO_BaP', 'POPPOCPO_PHE', 'POPPOCPO_PYR' )
1706 0 : Inst%IDTPOPPOCPO = HcoIDs(N)
1707 : CASE( 'POPPBCPO', 'POPPBCPO_BaP', 'POPPBCPO_PHE', 'POPPBCPO_PYR' )
1708 0 : Inst%IDTPOPPBCPO = HcoIDs(N)
1709 : CASE DEFAULT
1710 : ! Do nothing
1711 : END SELECT
1712 : ENDDO
1713 :
1714 : ! ERROR: POPG tracer is not found!
1715 0 : IF ( Inst%IDTPOPG <= 0 ) THEN
1716 0 : RC = HCO_FAIL
1717 0 : CALL HCO_ERROR( 'Cannot find POPG tracer in list of species!', RC )
1718 0 : RETURN
1719 : ENDIF
1720 :
1721 : ! ERROR! POPPOCPO tracer is not found
1722 0 : IF ( Inst%IDTPOPPOCPO <= 0 ) THEN
1723 0 : RC = HCO_FAIL
1724 0 : CALL HCO_ERROR( 'Cannot find POPPOCPO tracer in list of species!', RC )
1725 0 : RETURN
1726 : ENDIF
1727 :
1728 : ! ERROR! POPPBCPO tracer is not found
1729 0 : IF ( Inst%IDTPOPPBCPO <= 0 ) THEN
1730 0 : RC = HCO_FAIL
1731 0 : CALL HCO_ERROR( 'Cannot find POPPBCPO tracer in list of species!', RC )
1732 0 : RETURN
1733 : ENDIF
1734 :
1735 : !=======================================================================
1736 : ! Activate this module and the fields of ExtState that it uses
1737 : !=======================================================================
1738 :
1739 : ! Activate met fields required by this extension
1740 0 : ExtState%POPG%DoUse = .TRUE.
1741 0 : ExtState%AIRVOL%DoUse = .TRUE.
1742 0 : ExtState%AIRDEN%DoUse = .TRUE.
1743 0 : ExtState%FRAC_OF_PBL%DoUse = .TRUE.
1744 0 : ExtState%LAI%DoUse = .TRUE.
1745 0 : ExtState%PSC2_WET%DoUse = .TRUE.
1746 0 : ExtState%SNOWHGT%DoUse = .TRUE.
1747 0 : ExtState%TK%DoUse = .TRUE.
1748 0 : ExtState%TSKIN%DoUse = .TRUE.
1749 0 : ExtState%U10M%DoUse = .TRUE.
1750 0 : ExtState%V10M%DoUse = .TRUE.
1751 0 : ExtState%FRLAND%DoUse = .TRUE.
1752 0 : ExtState%FRLANDIC%DoUse = .TRUE.
1753 0 : ExtState%FROCEAN%DoUse = .TRUE.
1754 0 : ExtState%FRSEAICE%DoUse = .TRUE.
1755 0 : ExtState%FRLAKE%DoUse = .TRUE.
1756 :
1757 : !=======================================================================
1758 : ! Initialize data arrays
1759 : !=======================================================================
1760 :
1761 0 : ALLOCATE( Inst%EmisPOPG( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC )
1762 0 : IF ( RC /= 0 ) THEN
1763 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPG', RC )
1764 0 : RETURN
1765 : ENDIF
1766 0 : Inst%EmisPOPG = 0.0e0_hp
1767 :
1768 0 : ALLOCATE( Inst%EmisPOPPOCPO( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC )
1769 0 : IF ( RC /= 0 ) THEN
1770 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPPOCPO', RC )
1771 0 : RETURN
1772 : ENDIF
1773 0 : Inst%EmisPOPPOCPO = 0.0e0_hp
1774 :
1775 0 : ALLOCATE( Inst%EmisPOPPBCPO( HcoState%NX, HcoState%NY, HcoState%NZ ), STAT=RC )
1776 0 : IF ( RC /= 0 ) THEN
1777 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPPBCPO', RC )
1778 0 : RETURN
1779 : ENDIF
1780 0 : Inst%EmisPOPPBCPO = 0.0e0_hp
1781 :
1782 0 : ALLOCATE( Inst%EmisPOPGfromSoil( HcoState%NX, HcoState%NY ), STAT=RC )
1783 0 : IF ( RC /= 0 ) THEN
1784 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromSoil', RC )
1785 0 : RETURN
1786 : ENDIF
1787 0 : Inst%EmisPOPGfromSoil = 0.0e0_hp
1788 :
1789 0 : ALLOCATE( Inst%EmisPOPGfromLake( HcoState%NX, HcoState%NY ), STAT=RC )
1790 0 : IF ( RC /= 0 ) THEN
1791 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromLake', RC )
1792 0 : RETURN
1793 : ENDIF
1794 0 : Inst%EmisPOPGfromLake = 0.0e0_hp
1795 :
1796 0 : ALLOCATE( Inst%EmisPOPGfromLeaf( HcoState%NX, HcoState%NY ), STAT=RC )
1797 0 : IF ( RC /= 0 ) THEN
1798 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromLeaf', RC )
1799 0 : RETURN
1800 : ENDIF
1801 0 : Inst%EmisPOPGfromLeaf = 0.0e0_hp
1802 :
1803 0 : ALLOCATE( Inst%EmisPOPGfromSnow( HcoState%NX, HcoState%NY ), STAT=RC )
1804 0 : IF ( RC /= 0 ) THEN
1805 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromSnow', RC )
1806 0 : RETURN
1807 : ENDIF
1808 0 : Inst%EmisPOPGfromSnow = 0.0e0_hp
1809 :
1810 0 : ALLOCATE( Inst%EmisPOPGfromOcean( HcoState%NX, HcoState%NY ), STAT=RC )
1811 0 : IF ( RC /= 0 ) THEN
1812 0 : CALL HCO_ERROR ( 'Cannot allocate EmisPOPGfromOcean', RC )
1813 0 : RETURN
1814 : ENDIF
1815 0 : Inst%EmisPOPGfromOcean = 0.0e0_hp
1816 :
1817 0 : ALLOCATE( Inst%FluxPOPGfromSoilToAir( HcoState%NX, HcoState%NY ), STAT=RC )
1818 0 : IF ( RC /= 0 ) THEN
1819 0 : CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromSoilToAir', RC )
1820 0 : RETURN
1821 : ENDIF
1822 0 : Inst%FluxPOPGfromSoilToAir = 0.0e0_hp
1823 :
1824 0 : ALLOCATE( Inst%FluxPOPGfromAirToSoil( HcoState%NX, HcoState%NY ), STAT=RC )
1825 0 : IF ( RC /= 0 ) THEN
1826 0 : CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToSoil', RC )
1827 0 : RETURN
1828 : ENDIF
1829 0 : Inst%FluxPOPGfromAirToSoil = 0.0e0_hp
1830 :
1831 0 : ALLOCATE( Inst%FluxPOPGfromLakeToAir( HcoState%NX, HcoState%NY ), STAT=RC )
1832 0 : IF ( RC /= 0 ) THEN
1833 0 : CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromLakeToAir', RC )
1834 0 : RETURN
1835 : ENDIF
1836 0 : Inst%FluxPOPGfromLakeToAir = 0.0e0_hp
1837 :
1838 0 : ALLOCATE( Inst%FluxPOPGfromAirToLake( HcoState%NX, HcoState%NY ), STAT=RC )
1839 0 : IF ( RC /= 0 ) THEN
1840 0 : CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToLake', RC )
1841 0 : RETURN
1842 : ENDIF
1843 0 : Inst%FluxPOPGfromAirToLake = 0.0e0_hp
1844 :
1845 0 : ALLOCATE( Inst%FluxPOPGfromLeafToAir( HcoState%NX, HcoState%NY ), STAT=RC )
1846 0 : IF ( RC /= 0 ) THEN
1847 0 : CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromLeafToAir', RC )
1848 0 : RETURN
1849 : ENDIF
1850 0 : Inst%FluxPOPGfromLeafToAir = 0.0e0_hp
1851 :
1852 0 : ALLOCATE( Inst%FluxPOPGfromAirToLeaf( HcoState%NX, HcoState%NY ), STAT=RC )
1853 0 : IF ( RC /= 0 ) THEN
1854 0 : CALL HCO_ERROR ( 'Cannot allocate FluxPOPGfromAirToLeaf', RC )
1855 0 : RETURN
1856 : ENDIF
1857 0 : Inst%FluxPOPGfromAirToLeaf = 0.0e0_hp
1858 :
1859 0 : ALLOCATE( Inst%FugacitySoilToAir( HcoState%NX, HcoState%NY ), STAT=RC )
1860 0 : IF ( RC /= 0 ) THEN
1861 0 : CALL HCO_ERROR ( 'Cannot allocate FugacitySoilToAir', RC )
1862 0 : RETURN
1863 : ENDIF
1864 0 : Inst%FugacitySoilToAir = 0.0e0_hp
1865 :
1866 0 : ALLOCATE( Inst%FugacityLakeToAir( HcoState%NX, HcoState%NY ), STAT=RC )
1867 0 : IF ( RC /= 0 ) THEN
1868 0 : CALL HCO_ERROR ( 'Cannot allocate FugacityLakeToAir', RC )
1869 0 : RETURN
1870 : ENDIF
1871 0 : Inst%FugacityLakeToAir = 0.0e0_hp
1872 :
1873 0 : ALLOCATE( Inst%FugacityLeafToAir( HcoState%NX, HcoState%NY ), STAT=RC )
1874 0 : IF ( RC /= 0 ) THEN
1875 0 : CALL HCO_ERROR ( 'Cannot allocate FugacityLeafToAir', RC )
1876 0 : RETURN
1877 : ENDIF
1878 0 : Inst%FugacityLeafToAir = 0.0e0_hp
1879 :
1880 : !=======================================================================
1881 : ! Create manual diagnostics
1882 : !=======================================================================
1883 0 : DO I = 1,12
1884 :
1885 : ! Define diagnostic names. These have to match the names
1886 : ! in module HEMCO/Extensions/hcox_gc_POPs_mod.F90.
1887 : IF ( I == 1 ) THEN
1888 0 : DiagnName = 'EmisPOPGfromSoil'
1889 0 : OutUnit = 'kg/m2/s'
1890 0 : Target2D => Inst%EmisPOPGfromSoil
1891 : ELSEIF ( I == 2 ) THEN
1892 0 : DiagnName = 'EmisPOPGfromLake'
1893 0 : OutUnit = 'kg/m2/s'
1894 0 : Target2D => Inst%EmisPOPGfromLake
1895 : ELSEIF ( I == 3 ) THEN
1896 0 : DiagnName = 'EmisPOPGfromLeaf'
1897 0 : OutUnit = 'kg/m2/s'
1898 0 : Target2D => Inst%EmisPOPGfromLeaf
1899 : ELSEIF ( I == 4 ) THEN
1900 0 : DiagnName = 'FluxPOPGfromSoilToAir'
1901 0 : OutUnit = 'ng/m2/day'
1902 0 : Target2D => Inst%FluxPOPGfromSoilToAir
1903 : ELSEIF ( I == 5 ) THEN
1904 0 : DiagnName = 'FluxPOPGfromAirToSoil'
1905 0 : OutUnit = 'ng/m2/day'
1906 0 : Target2D => Inst%FluxPOPGfromAirToSoil
1907 : ELSEIF ( I == 6 ) THEN
1908 0 : DiagnName = 'FluxPOPGfromLakeToAir'
1909 0 : OutUnit = 'ng/m2/day'
1910 0 : Target2D => Inst%FluxPOPGfromLakeToAir
1911 : ELSEIF ( I == 7 ) THEN
1912 0 : DiagnName = 'FluxPOPGfromAirToLake'
1913 0 : OutUnit = 'ng/m2/day'
1914 0 : Target2D => Inst%FluxPOPGfromAirToLake
1915 : ELSEIF ( I == 8 ) THEN
1916 0 : DiagnName = 'FluxPOPGfromLeafToAir'
1917 0 : OutUnit = 'ng/m2/day'
1918 0 : Target2D => Inst%FluxPOPGfromLeafToAir
1919 : ELSEIF ( I == 9 ) THEN
1920 0 : DiagnName = 'FluxPOPGfromAirToLeaf'
1921 0 : OutUnit = 'ng/m2/day'
1922 0 : Target2D => Inst%FluxPOPGfromAirToLeaf
1923 : ELSEIF ( I == 10 ) THEN
1924 0 : DiagnName = 'FugacitySoilToAir'
1925 0 : OutUnit = '1'
1926 0 : Target2D => Inst%FugacitySoilToAir
1927 : ELSEIF ( I == 11 ) THEN
1928 0 : DiagnName = 'FugacityLakeToAir'
1929 0 : OutUnit = '1'
1930 0 : Target2D => Inst%FugacityLakeToAir
1931 : ELSEIF ( I == 12 ) THEN
1932 0 : DiagnName = 'FugacityLeafToAir'
1933 0 : OutUnit = '1'
1934 0 : Target2D => Inst%FugacityLeafToAir
1935 : ENDIF
1936 :
1937 : ! Create manual diagnostics
1938 : CALL Diagn_Create( HcoState = HcoState, &
1939 : cName = TRIM( DiagnName ), &
1940 : ExtNr = ExtNr, &
1941 : Cat = -1, &
1942 : Hier = -1, &
1943 : HcoID = -1, &
1944 : SpaceDim = 2, &
1945 : OutUnit = OutUnit, &
1946 : AutoFill = 0, &
1947 : Trgt2D = Target2D, &
1948 0 : RC = RC )
1949 0 : IF ( RC /= HCO_SUCCESS ) THEN
1950 0 : CALL HCO_ERROR( 'ERROR 8', RC, THISLOC=LOC )
1951 0 : RETURN
1952 : ENDIF
1953 0 : Target2D => NULL()
1954 :
1955 : ENDDO
1956 :
1957 : !=======================================================================
1958 : ! Leave w/ success
1959 : !=======================================================================
1960 0 : IF ( ALLOCATED( HcoIDs ) ) DEALLOCATE( HcoIDs )
1961 0 : IF ( ALLOCATED( SpcNames ) ) DEALLOCATE( SpcNames )
1962 :
1963 : ! Nullify pointers
1964 0 : Inst => NULL()
1965 :
1966 0 : CALL HCO_LEAVE( HcoState%Config%Err,RC )
1967 :
1968 0 : END SUBROUTINE HCOX_GC_POPs_Init
1969 : !EOC
1970 : !------------------------------------------------------------------------------
1971 : ! Harmonized Emissions Component (HEMCO) !
1972 : !------------------------------------------------------------------------------
1973 : !BOP
1974 : !
1975 : ! !IROUTINE: HCOX_GC_POPs_Final
1976 : !
1977 : ! !DESCRIPTION: Subroutine HcoX\_GC\_POPs\_Final finalizes the HEMCO
1978 : ! extension for the GEOS-Chem POPs specialty simulation. All module
1979 : ! arrays will be deallocated.
1980 : !\\
1981 : !\\
1982 : ! !INTERFACE:
1983 : !
1984 0 : SUBROUTINE HCOX_GC_POPs_Final( ExtState )
1985 : !
1986 : ! !INPUT PARAMETERS:
1987 : !
1988 : TYPE(Ext_State), POINTER :: ExtState ! Module options
1989 : !
1990 : ! !REVISION HISTORY:
1991 : ! 19 Aug 2014 - M. Sulprizio- Initial version
1992 : ! See https://github.com/geoschem/hemco for complete history
1993 : !EOP
1994 : !------------------------------------------------------------------------------
1995 : !BOC
1996 :
1997 : !=======================================================================
1998 : ! HCOX_GC_POPs_FINAL begins here!
1999 : !=======================================================================
2000 :
2001 0 : CALL InstRemove( ExtState%GC_POPs )
2002 :
2003 0 : END SUBROUTINE HCOX_GC_POPs_Final
2004 : !EOC
2005 : !------------------------------------------------------------------------------
2006 : ! Harmonized Emissions Component (HEMCO) !
2007 : !------------------------------------------------------------------------------
2008 : !BOP
2009 : !
2010 : ! !IROUTINE: InstGet
2011 : !
2012 : ! !DESCRIPTION: Subroutine InstGet returns a poiner to the desired instance.
2013 : !\\
2014 : !\\
2015 : ! !INTERFACE:
2016 : !
2017 0 : SUBROUTINE InstGet ( Instance, Inst, RC, PrevInst )
2018 : !
2019 : ! !INPUT PARAMETERS:
2020 : !
2021 : INTEGER :: Instance
2022 : TYPE(MyInst), POINTER :: Inst
2023 : INTEGER :: RC
2024 : TYPE(MyInst), POINTER, OPTIONAL :: PrevInst
2025 : !
2026 : ! !REVISION HISTORY:
2027 : ! 18 Feb 2016 - C. Keller - Initial version
2028 : ! See https://github.com/geoschem/hemco for complete history
2029 : !EOP
2030 : !------------------------------------------------------------------------------
2031 : !BOC
2032 : TYPE(MyInst), POINTER :: PrvInst
2033 :
2034 : !=================================================================
2035 : ! InstGet begins here!
2036 : !=================================================================
2037 :
2038 : ! Get instance. Also archive previous instance.
2039 0 : PrvInst => NULL()
2040 0 : Inst => AllInst
2041 0 : DO WHILE ( ASSOCIATED(Inst) )
2042 0 : IF ( Inst%Instance == Instance ) EXIT
2043 0 : PrvInst => Inst
2044 0 : Inst => Inst%NextInst
2045 : END DO
2046 0 : IF ( .NOT. ASSOCIATED( Inst ) ) THEN
2047 0 : RC = HCO_FAIL
2048 0 : RETURN
2049 : ENDIF
2050 :
2051 : ! Pass output arguments
2052 0 : IF ( PRESENT(PrevInst) ) PrevInst => PrvInst
2053 :
2054 : ! Cleanup & Return
2055 0 : PrvInst => NULL()
2056 0 : RC = HCO_SUCCESS
2057 :
2058 : END SUBROUTINE InstGet
2059 : !EOC
2060 : !------------------------------------------------------------------------------
2061 : ! Harmonized Emissions Component (HEMCO) !
2062 : !------------------------------------------------------------------------------
2063 : !BOP
2064 : !
2065 : ! !IROUTINE: InstCreate
2066 : !
2067 : ! !DESCRIPTION: Subroutine InstCreate creates a new instance.
2068 : !\\
2069 : !\\
2070 : ! !INTERFACE:
2071 : !
2072 0 : SUBROUTINE InstCreate ( ExtNr, Instance, Inst, RC )
2073 : !
2074 : ! !INPUT PARAMETERS:
2075 : !
2076 : INTEGER, INTENT(IN) :: ExtNr
2077 : !
2078 : ! !OUTPUT PARAMETERS:
2079 : !
2080 : INTEGER, INTENT( OUT) :: Instance
2081 : TYPE(MyInst), POINTER :: Inst
2082 : !
2083 : ! !INPUT/OUTPUT PARAMETERS:
2084 : !
2085 : INTEGER, INTENT(INOUT) :: RC
2086 : !
2087 : ! !REVISION HISTORY:
2088 : ! 18 Feb 2016 - C. Keller - Initial version
2089 : ! See https://github.com/geoschem/hemco for complete history
2090 : !EOP
2091 : !------------------------------------------------------------------------------
2092 : !BOC
2093 : TYPE(MyInst), POINTER :: TmpInst
2094 : INTEGER :: nnInst
2095 :
2096 : !=================================================================
2097 : ! InstCreate begins here!
2098 : !=================================================================
2099 :
2100 : ! ----------------------------------------------------------------
2101 : ! Generic instance initialization
2102 : ! ----------------------------------------------------------------
2103 :
2104 : ! Initialize
2105 0 : Inst => NULL()
2106 :
2107 : ! Get number of already existing instances
2108 0 : TmpInst => AllInst
2109 0 : nnInst = 0
2110 0 : DO WHILE ( ASSOCIATED(TmpInst) )
2111 0 : nnInst = nnInst + 1
2112 0 : TmpInst => TmpInst%NextInst
2113 : END DO
2114 :
2115 : ! Create new instance
2116 0 : ALLOCATE(Inst)
2117 0 : Inst%Instance = nnInst + 1
2118 0 : Inst%ExtNr = ExtNr
2119 :
2120 : ! Attach to instance list
2121 0 : Inst%NextInst => AllInst
2122 0 : AllInst => Inst
2123 :
2124 : ! Update output instance
2125 0 : Instance = Inst%Instance
2126 :
2127 : ! ----------------------------------------------------------------
2128 : ! Type specific initialization statements follow below
2129 : ! ----------------------------------------------------------------
2130 :
2131 : ! Return w/ success
2132 0 : RC = HCO_SUCCESS
2133 :
2134 0 : END SUBROUTINE InstCreate
2135 : !EOC
2136 : !------------------------------------------------------------------------------
2137 : ! Harmonized Emissions Component (HEMCO) !
2138 : !------------------------------------------------------------------------------
2139 : !BOP
2140 : !
2141 : ! !IROUTINE: InstRemove
2142 : !
2143 : ! !DESCRIPTION: Subroutine InstRemove creates a new instance.
2144 : !\\
2145 : !\\
2146 : ! !INTERFACE:
2147 : !
2148 0 : SUBROUTINE InstRemove ( Instance )
2149 : !
2150 : ! !INPUT PARAMETERS:
2151 : !
2152 : INTEGER :: Instance
2153 : !
2154 : ! !REVISION HISTORY:
2155 : ! 18 Feb 2016 - C. Keller - Initial version
2156 : ! See https://github.com/geoschem/hemco for complete history
2157 : !EOP
2158 : !------------------------------------------------------------------------------
2159 : !BOC
2160 : INTEGER :: RC
2161 : TYPE(MyInst), POINTER :: PrevInst
2162 : TYPE(MyInst), POINTER :: Inst
2163 :
2164 : !=================================================================
2165 : ! InstRemove begins here!
2166 : !=================================================================
2167 :
2168 : ! Init
2169 0 : PrevInst => NULL()
2170 0 : Inst => NULL()
2171 :
2172 : ! Get instance. Also archive previous instance.
2173 0 : CALL InstGet ( Instance, Inst, RC, PrevInst=PrevInst )
2174 :
2175 : ! Instance-specific deallocation
2176 0 : IF ( ASSOCIATED(Inst) ) THEN
2177 :
2178 : !---------------------------------------------------------------------
2179 : ! Deallocate fields of Inst before popping off from the list
2180 : ! in order to avoid memory leaks (Bob Yantosca (17 Aug 2022)
2181 : !---------------------------------------------------------------------
2182 0 : IF ( ASSOCIATED( Inst%EmisPOPPOCPO ) ) THEN
2183 0 : DEALLOCATE( Inst%EmisPOPPOCPO )
2184 : ENDIF
2185 0 : Inst%EmisPOPPOCPO => NULL()
2186 :
2187 0 : IF ( ASSOCIATED( Inst%EmisPOPPBCPO ) ) THEN
2188 0 : DEALLOCATE( Inst%EmisPOPPBCPO )
2189 : ENDIF
2190 0 : Inst%EmisPOPPBCPO => NULL()
2191 :
2192 0 : IF ( ASSOCIATED( Inst%EmisPOPG ) ) THEN
2193 0 : DEALLOCATE( Inst%EmisPOPG )
2194 : ENDIF
2195 0 : Inst%EmisPOPG => NULL()
2196 :
2197 0 : IF ( ASSOCIATED( Inst%EmisPOPGfromSoil ) ) THEN
2198 0 : DEALLOCATE( Inst%EmisPOPGfromSoil )
2199 : ENDIF
2200 0 : Inst%EmisPOPGfromSoil => NULL()
2201 :
2202 0 : IF ( ASSOCIATED( Inst%EmisPOPGfromLake ) ) THEN
2203 0 : DEALLOCATE( Inst%EmisPOPGfromLake )
2204 : ENDIF
2205 0 : Inst%EmisPOPGfromLake => NULL()
2206 :
2207 0 : IF ( ASSOCIATED( Inst%EmisPOPGfromLeaf ) ) THEN
2208 0 : DEALLOCATE( Inst%EmisPOPGfromLeaf )
2209 : ENDIF
2210 0 : Inst%EmisPOPGfromLeaf => NULL()
2211 :
2212 0 : IF ( ASSOCIATED( Inst%EmisPOPGfromSnow ) ) THEN
2213 0 : DEALLOCATE( Inst%EmisPOPGfromSnow )
2214 : ENDIF
2215 0 : Inst%EmisPOPGfromSnow => NULL()
2216 :
2217 0 : IF ( ASSOCIATED( Inst%EmisPOPGfromOcean ) ) THEN
2218 0 : DEALLOCATE( Inst%EmisPOPGfromOcean )
2219 : ENDIF
2220 0 : Inst%EmisPOPGfromOcean => NULL()
2221 :
2222 0 : IF ( ASSOCIATED( Inst%FluxPOPGfromSoilToAir ) ) THEN
2223 0 : DEALLOCATE( Inst%FluxPOPGfromSoilToAir )
2224 : ENDIF
2225 0 : Inst%FluxPOPGfromSoilToAir => NULL()
2226 :
2227 0 : IF ( ASSOCIATED( Inst%FluxPOPGfromAirToSoil ) ) THEN
2228 0 : DEALLOCATE( Inst%FluxPOPGfromAirToSoil )
2229 : ENDIF
2230 0 : Inst%FluxPOPGfromAirToSoil => NULL()
2231 :
2232 0 : IF ( ASSOCIATED( Inst%FluxPOPGfromLakeToAir ) ) THEN
2233 0 : DEALLOCATE( Inst%FluxPOPGfromLakeToAir )
2234 : ENDIF
2235 0 : Inst%FluxPOPGfromLakeToAir => NULL()
2236 :
2237 0 : IF ( ASSOCIATED( Inst%FluxPOPGfromAirToLake ) ) THEN
2238 0 : DEALLOCATE( Inst%FluxPOPGfromAirToLake )
2239 : ENDIF
2240 0 : Inst%FluxPOPGfromAirToLake => NULL()
2241 :
2242 0 : IF ( ASSOCIATED( Inst%FluxPOPGfromLeafToAir ) ) THEN
2243 0 : DEALLOCATE( Inst%FluxPOPGfromLeafToAir )
2244 : ENDIF
2245 0 : Inst%FluxPOPGfromLeafToAir => NULL()
2246 :
2247 0 : IF ( ASSOCIATED( Inst%FluxPOPGfromAirtoLeaf ) ) THEN
2248 0 : DEALLOCATE( Inst%FluxPOPGfromAirtoLeaf )
2249 : ENDIF
2250 0 : Inst%FluxPOPGfromAirtoLeaf => NULL()
2251 :
2252 0 : IF ( ASSOCIATED( Inst%FugacitySoilToAir ) ) THEN
2253 0 : DEALLOCATE( Inst%FugacitySoilToAir )
2254 : ENDIF
2255 0 : Inst%FugacitySoilToAir => NULL()
2256 :
2257 0 : IF ( ASSOCIATED( Inst%FugacityLakeToAir ) ) THEN
2258 0 : DEALLOCATE( Inst%FugacityLakeToAir )
2259 : ENDIF
2260 0 : Inst%FugacityLakeToAir => NULL()
2261 :
2262 0 : IF ( ASSOCIATED( Inst%FugacityLeafToAir ) ) THEN
2263 0 : DEALLOCATE( Inst%FugacityLeafToAir )
2264 : ENDIF
2265 0 : Inst%FugacityLeafToAir => NULL()
2266 :
2267 : !---------------------------------------------------------------------
2268 : ! Pop off instance from list
2269 : !---------------------------------------------------------------------
2270 0 : IF ( ASSOCIATED(PrevInst) ) THEN
2271 0 : PrevInst%NextInst => Inst%NextInst
2272 : ELSE
2273 0 : AllInst => Inst%NextInst
2274 : ENDIF
2275 0 : DEALLOCATE(Inst)
2276 : ENDIF
2277 :
2278 : ! Free pointers before exiting
2279 0 : PrevInst => NULL()
2280 0 : Inst => NULL()
2281 :
2282 0 : END SUBROUTINE InstRemove
2283 : !EOC
2284 0 : END MODULE HCOX_GC_POPs_Mod
|
