Line data Source code
1 : !------------------------------------------------------------------------------
2 : ! Harmonized Emissions Component (HEMCO) !
3 : !------------------------------------------------------------------------------
4 : !BOP
5 : !
6 : ! !MODULE: hcox_gc_RnPbBe_mod.F90
7 : !
8 : ! !DESCRIPTION: Defines the HEMCO extension for the GEOS-Chem Rn-Pb-Be
9 : ! specialty simulation.
10 : !\\
11 : !\\
12 : ! This extension parameterizes emissions of Rn and/or Pb based upon the
13 : ! literature given below. The emission fields become automatically added
14 : ! to the HEMCO emission array of the given species. It is possible to
15 : ! select only one of the two species (Rn or Pb) in the HEMCO configuration
16 : ! file. This may be useful if a gridded data inventory shall be applied to
17 : ! one of the species (through the standard HEMCO interface).
18 : !\\
19 : !\\
20 : ! !INTERFACE:
21 : !
22 : MODULE HCOX_GC_RnPbBe_Mod
23 : !
24 : ! !USES:
25 : !
26 : USE HCO_Error_Mod
27 : USE HCO_Diagn_Mod
28 : USE HCO_State_Mod, ONLY : HCO_State ! Derived type for HEMCO state
29 : USE HCOX_State_Mod, ONLY : Ext_State ! Derived type for External state
30 :
31 : IMPLICIT NONE
32 : PRIVATE
33 : !
34 : ! !PUBLIC MEMBER FUNCTIONS:
35 : !
36 : PUBLIC :: HcoX_GC_RnPbBe_Run
37 : PUBLIC :: HcoX_GC_RnPbBe_Init
38 : PUBLIC :: HcoX_Gc_RnPbBe_Final
39 : !
40 : ! !PRIVATE MEMBER FUNCTIONS:
41 : !
42 : PRIVATE :: Init_7Be_Emissions
43 : !
44 : ! !REMARKS:
45 : ! References:
46 : ! ============================================================================
47 : ! (1 ) Liu,H., D.Jacob, I.Bey, and R.M.Yantosca, Constraints from 210Pb
48 : ! and 7Be on wet deposition and transport in a global three-dimensional
49 : ! chemical tracer model driven by assimilated meteorological fields,
50 : ! JGR, 106, D11, 12,109-12,128, 2001.
51 : ! (2 ) Jacob et al.,Evaluation and intercomparison of global atmospheric
52 : ! transport models using Rn-222 and other short-lived tracers,
53 : ! JGR, 1997 (102):5953-5970
54 : ! (3 ) Dorothy Koch, JGR 101, D13, 18651, 1996.
55 : ! (4 ) Lal, D., and B. Peters, Cosmic ray produced radioactivity on the
56 : ! Earth. Handbuch der Physik, 46/2, 551-612, edited by K. Sitte,
57 : ! Springer-Verlag, New York, 1967.
58 : ! (5 ) Koch and Rind, Beryllium 10/beryllium 7 as a tracer of stratospheric
59 : ! transport, JGR, 103, D4, 3907-3917, 1998.
60 : !
61 : ! !REVISION HISTORY:
62 : ! 07 Jul 2014 - R. Yantosca - Initial version
63 : ! See https://github.com/geoschem/hemco for complete history
64 : !EOP
65 : !------------------------------------------------------------------------------
66 : !BOC
67 : !
68 : ! !PRIVATE TYPES:
69 : !
70 : TYPE :: MyInst
71 :
72 : ! Emissions indices etc.
73 : INTEGER :: Instance
74 : INTEGER :: ExtNr ! Main Extension number
75 : INTEGER :: ExtNrZhang ! ZHANG_Rn222 extension number
76 : INTEGER :: IDTRn222 ! Index # for Rn222
77 : INTEGER :: IDTBe7 ! Index # for Be7
78 : INTEGER :: IDTBe7s ! Index # for Be7s
79 : INTEGER :: IDTBe10 ! Index # for Be10
80 : INTEGER :: IDTBe10s ! Index # for Be10s
81 :
82 : ! For tracking Rn222, Be7, and Be10 emissions
83 : REAL(hp), POINTER :: EmissRn222(:,: )
84 : REAL(hp), POINTER :: EmissBe7 (:,:,:)
85 : REAL(hp), POINTER :: EmissBe7s (:,:,:)
86 : REAL(hp), POINTER :: EmissBe10 (:,:,:)
87 : REAL(hp), POINTER :: EmissBe10s(:,:,:)
88 :
89 : ! For Lal & Peters 7Be emissions input data
90 : REAL(hp), POINTER :: LATSOU(: ) ! Array for latitudes
91 : REAL(hp), POINTER :: PRESOU(: ) ! Array for pressures
92 : REAL(hp), POINTER :: BESOU (:,: ) ! Array for 7Be emissions
93 :
94 : TYPE(MyInst), POINTER :: NextInst => NULL()
95 : END TYPE MyInst
96 :
97 : ! Pointer to instances
98 : TYPE(MyInst), POINTER :: AllInst => NULL()
99 : !
100 : ! !DEFINED PARAMETERS:
101 : !
102 : ! To convert kg to atoms
103 : REAL*8, PARAMETER :: XNUMOL_Rn = ( 6.022140857d23 / 222.0d-3 )
104 : REAL*8, PARAMETER :: XNUMOL_Be7 = ( 6.022140857d23 / 7.0d-3 )
105 : REAL*8, PARAMETER :: XNUMOL_Be10 = ( 6.022140857d23 / 10.0d-3 )
106 :
107 : CONTAINS
108 : !EOC
109 : !------------------------------------------------------------------------------
110 : ! Harmonized Emissions Component (HEMCO) !
111 : !------------------------------------------------------------------------------
112 : !BOP
113 : !
114 : ! !IROUTINE: HCOX_Gc_RnPbBe_run
115 : !
116 : ! !DESCRIPTION: Subroutine HcoX\_Gc\_RnPbBe\_Run computes emissions of 222Rn,
117 : ! 7Be, and 10Be for the GEOS-Chem Rn-Pb-Be specialty simulation.
118 : !\\
119 : !\\
120 : ! !INTERFACE:
121 : !
122 0 : SUBROUTINE HCOX_Gc_RnPbBe_Run( ExtState, HcoState, RC )
123 : !
124 : ! !USES:
125 : !
126 : USE HCO_Calc_Mod, ONLY : HCO_EvalFld
127 : USE HCO_FluxArr_Mod, ONLY : HCO_EmisAdd
128 : !
129 : ! !INPUT PARAMETERS:
130 : !
131 : TYPE(Ext_State), POINTER :: ExtState ! Options for Rn-Pb-Be sim
132 : TYPE(HCO_State), POINTER :: HcoState ! HEMCO state
133 : !
134 : ! !INPUT/OUTPUT PARAMETERS:
135 : !
136 : INTEGER, INTENT(INOUT) :: RC ! Success or failure?
137 : !
138 : ! !REMARKS:
139 : ! This code is based on routine EMISSRnPbBe in prior versions of GEOS-Chem.
140 : !
141 : ! !REVISION HISTORY:
142 : ! 07 Jul 2014 - R. Yantosca - Initial version
143 : ! See https://github.com/geoschem/hemco for complete history
144 : !EOP
145 : !------------------------------------------------------------------------------
146 : !BOC
147 : !
148 : ! !LOCAL VARIABLES:
149 : !
150 :
151 : ! Scalars
152 : INTEGER :: I, J, L, N
153 : INTEGER :: HcoID
154 : REAL*8 :: A_CM2, ADD_Rn, Add_Be7, Add_Be10
155 : REAL*8 :: Rn_LAND, Rn_WATER, DTSRCE
156 : REAL*8 :: Rn_TMP, LAT, F_LAND
157 : REAL*8 :: F_WATER, F_BELOW_70, F_BELOW_60, F_ABOVE_60
158 : REAL*8 :: DENOM
159 : REAL(hp) :: LAT_TMP, P_TMP, Be_TMP
160 : CHARACTER(LEN=255):: MSG, LOC
161 :
162 : ! Pointers
163 : TYPE(MyInst), POINTER :: Inst
164 0 : REAL(hp), POINTER :: Arr2D(:,: )
165 0 : REAL(hp), POINTER :: Arr3D(:,:,:)
166 :
167 : !=======================================================================
168 : ! HCOX_GC_RnPbBe_RUN begins here!
169 : !=======================================================================
170 0 : LOC = 'HCOX_GC_RnPbBe_RUN (HCOX_GC_RNPBBE_MOD.F90)'
171 :
172 : ! Return if extension not turned on
173 0 : IF ( ExtState%GC_RnPbBe <= 0 ) RETURN
174 :
175 : ! Enter
176 0 : CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
177 0 : IF ( RC /= HCO_SUCCESS ) THEN
178 0 : CALL HCO_ERROR( 'ERROR 0', RC, THISLOC=LOC )
179 0 : RETURN
180 : ENDIF
181 :
182 : ! Set error flag
183 : !ERR = .FALSE.
184 :
185 : ! Get instance
186 0 : Inst => NULL()
187 0 : CALL InstGet ( ExtState%GC_RnPbBe, Inst, RC )
188 0 : IF ( RC /= HCO_SUCCESS ) THEN
189 0 : WRITE(MSG,*) 'Cannot find GC_RnPbBe instance Nr. ', ExtState%GC_RnPbBe
190 0 : CALL HCO_ERROR(MSG,RC)
191 0 : RETURN
192 : ENDIF
193 :
194 : ! Emission timestep [s]
195 0 : DTSRCE = HcoState%TS_EMIS
196 :
197 : ! Nullify
198 0 : Arr2D => NULL()
199 0 : Arr3D => NULL()
200 :
201 : !=======================================================================
202 : ! Compute 222Rn emissions [kg/m2/s], according to the following:
203 : !
204 : ! (1) 222Rn emission poleward of 70 degrees = 0.0 [atoms/cm2/s]
205 : !
206 : ! (2) For latitudes 70S-60S and 60N-70N (both land & ocean),
207 : ! 222Rn emission is 0.005 [atoms/cm2/s]
208 : !
209 : ! (3) For latitudes between 60S and 60N,
210 : ! 222Rn emission is 1 [atoms/cm2/s] over land or
211 : ! 0.005 [atoms/cm2/s] over oceans
212 : !
213 : ! (4) For grid boxes where the surface temperature is below
214 : ! 0 deg Celsius, reduce 222Rn emissions by a factor of 3.
215 : !
216 : ! Reference: Jacob et al.,Evaluation and intercomparison of
217 : ! global atmospheric transport models using Rn-222 and other
218 : ! short-lived tracers, JGR, 1997 (102):5953-5970
219 : !=======================================================================
220 0 : IF ( Inst%IDTRn222 > 0 ) THEN
221 :
222 0 : IF ( Inst%ExtNrZhang > 0 ) THEN
223 :
224 : !------------------------------------------------------------------
225 : ! Use Zhang et al Rn222 emissions
226 : ! cf https://doi.org/10.5194/acp-21-1861-2021
227 : !------------------------------------------------------------------
228 : CALL HCO_EvalFld( HcoState, 'ZHANG_Rn222_EMIS', &
229 0 : Inst%EmissRn222, RC )
230 0 : IF ( RC /= HCO_SUCCESS ) THEN
231 0 : CALL HCO_Error( 'Could not read ZHANG_Rn222_EMIS!', RC )
232 0 : RETURN
233 : ENDIF
234 :
235 : ELSE
236 :
237 : !------------------------------------------------------------------
238 : ! Use default Rn222 emissions, based on Jacob et al 1997
239 : !------------------------------------------------------------------
240 : !$OMP PARALLEL DO &
241 : !$OMP DEFAULT( SHARED ) &
242 : !$OMP PRIVATE( I, J, LAT, DENOM ) &
243 : !$OMP PRIVATE( F_BELOW_70, F_BELOW_60, F_ABOVE_60, Rn_LAND ) &
244 : !$OMP PRIVATE( Rn_WATER, F_LAND, F_WATER, ADD_Rn ) &
245 : !$OMP SCHEDULE( DYNAMIC )
246 0 : DO J = 1, HcoState%Ny
247 0 : DO I = 1, HcoState%Nx
248 :
249 : ! Get ABS( latitude ) of the grid box
250 0 : LAT = ABS( HcoState%Grid%YMID%Val( I, J ) )
251 :
252 : ! Zero for safety's sake
253 0 : F_BELOW_70 = 0d0
254 0 : F_BELOW_60 = 0d0
255 0 : F_ABOVE_60 = 0d0
256 :
257 : ! Baseline 222Rn emissions
258 : ! Rn_LAND [kg/m2/s] = [1 atom 222Rn/cm2/s] / [atoms/kg]
259 : ! * [1d4 cm2/m2]
260 0 : Rn_LAND = ( 1d0 / XNUMOL_Rn ) * 1d4
261 :
262 : ! Baseline 222Rn emissions over water or ice [kg]
263 0 : Rn_WATER = Rn_LAND * 0.005d0
264 :
265 : ! Fraction of grid box that is land
266 0 : F_LAND = ExtState%FRCLND%Arr%Val(I,J)
267 :
268 : ! Fraction of grid box that is water
269 0 : F_WATER = 1d0 - F_LAND
270 :
271 : !--------------------
272 : ! 90S-70S or 70N-90N
273 : !--------------------
274 0 : IF ( LAT >= 70d0 ) THEN
275 :
276 : ! 222Rn emissions are shut off poleward of 70 degrees
277 : ADD_Rn = 0.0d0
278 :
279 : !--------------------
280 : ! 70S-60S or 60N-70N
281 : !--------------------
282 0 : ELSE IF ( LAT >= 60d0 ) THEN
283 :
284 0 : IF ( LAT <= 70d0 ) THEN
285 :
286 : ! If the entire grid box lies equatorward of 70 deg,
287 : ! then 222Rn emissions here are 0.005 [atoms/cm2/s]
288 : ADD_Rn = Rn_WATER
289 :
290 : ELSE
291 :
292 : ! N-S extent of grid box [degrees]
293 0 : DENOM = HcoState%Grid%YMID%Val( I, J+1 ) &
294 0 : - HcoState%Grid%YMID%Val( I, J )
295 :
296 : ! Compute the fraction of the grid box below 70 degrees
297 0 : F_BELOW_70 = ( 70.0d0 - LAT ) / DENOM
298 :
299 : ! If the grid box straddles the 70S or 70N latitude
300 : ! line, then only count 222Rn emissions equatorward of
301 : ! 70 degrees. 222Rn emissions here are 0.005
302 : ! [atoms/cm2/s].
303 0 : ADD_Rn = F_BELOW_70 * Rn_WATER
304 :
305 : ENDIF
306 :
307 : ELSE
308 :
309 : !--------------------
310 : ! 70S-60S or 60N-70N
311 : !--------------------
312 0 : IF ( LAT > 60d0 ) THEN
313 :
314 : ! N-S extent of grid box [degrees]
315 0 : DENOM = HcoState%Grid%YMID%Val( I, J+1 ) &
316 0 : - HcoState%Grid%YMID%Val( I, J )
317 :
318 : ! Fraction of grid box with ABS( lat ) below 60 degrees
319 0 : F_BELOW_60 = ( 60.0d0 - LAT ) / DENOM
320 :
321 : ! Fraction of grid box with ABS( lat ) above 60 degrees
322 0 : F_ABOVE_60 = F_BELOW_60
323 :
324 : ADD_Rn = &
325 : ! Consider 222Rn emissions equatorward of
326 : ! 60 degrees for both land (1.0 [atoms/cm2/s])
327 : ! and water (0.005 [atoms/cm2/s])
328 : F_BELOW_60 * &
329 : ( Rn_LAND * F_LAND ) + &
330 : ( Rn_WATER * F_WATER ) + &
331 :
332 : ! If the grid box straddles the 60 degree boundary
333 : ! then also consider the emissions poleward of 60
334 : ! degrees. 222Rn emissions here are 0.005
335 : ! [atoms/cm2/s].
336 0 : F_ABOVE_60 * Rn_WATER
337 :
338 : !--------------------
339 : ! 60S-60N
340 : !--------------------
341 : ELSE
342 :
343 : ! Consider 222Rn emissions equatorward of 60 deg for
344 : ! land (1.0 [atoms/cm2/s]) and water (0.005 [atoms/cm2/s])
345 0 : ADD_Rn = ( Rn_LAND * F_LAND ) + ( Rn_WATER * F_WATER )
346 :
347 : ENDIF
348 : ENDIF
349 :
350 : ! For boxes below freezing, reduce 222Rn emissions by 3x
351 0 : IF ( ExtState%T2M%Arr%Val(I,J) < 273.15 ) THEN
352 0 : ADD_Rn = ADD_Rn / 3d0
353 : ENDIF
354 :
355 : ! Save 222Rn emissions into an array [kg/m2/s]
356 0 : Inst%EmissRn222(I,J) = ADD_Rn
357 : ENDDO
358 : ENDDO
359 : !$OMP END PARALLEL DO
360 :
361 : ENDIF
362 :
363 : !------------------------------------------------------------------------
364 : ! Add 222Rn emissions to HEMCO data structure & diagnostics
365 : !------------------------------------------------------------------------
366 :
367 : ! Add emissions
368 0 : Arr2D => Inst%EmissRn222(:,:)
369 : CALL HCO_EmisAdd( HcoState, Arr2D, Inst%IDTRn222, &
370 0 : RC, ExtNr=Inst%ExtNr )
371 0 : Arr2D => NULL()
372 0 : IF ( RC /= HCO_SUCCESS ) THEN
373 : CALL HCO_ERROR( &
374 0 : 'HCO_EmisAdd error: EmissRn222', RC )
375 0 : RETURN
376 : ENDIF
377 :
378 : ENDIF ! IDTRn222 > 0
379 :
380 : !=======================================================================
381 : ! Compute 7Be and 10Be emissions [kg/m2/s]
382 : !
383 : ! Original units of 7Be and 10Be emissions are [stars/g air/sec],
384 : ! where "stars" = # of nuclear disintegrations of cosmic rays
385 : !
386 : ! Now interpolate from 33 std levels onto GEOS-CHEM levels
387 : !
388 : ! 7Be and 10Be have identical source distributions (Koch and Rind, 1998)
389 : !=======================================================================
390 0 : IF ( Inst%IDTBe7 > 0 .or. Inst%IDTBe10 > 0 ) THEN
391 : !$OMP PARALLEL DO &
392 : !$OMP DEFAULT( SHARED ) &
393 : !$OMP PRIVATE( I, J, L, LAT_TMP, P_TMP, Be_TMP, ADD_Be7, ADD_Be10 ) &
394 : !$OMP SCHEDULE( DYNAMIC )
395 0 : DO L = 1, HcoState%Nz
396 0 : DO J = 1, HcoState%Ny
397 0 : DO I = 1, HcoState%Nx
398 :
399 : ! Get absolute value of latitude, since we will assume that
400 : ! the 7Be distribution is symmetric about the equator
401 0 : LAT_TMP = ABS( HcoState%Grid%YMID%Val( I, J ) )
402 :
403 : ! Pressure at (I,J,L) [hPa]
404 : ! Now calculate from edge points (ckeller, 10/06/1014)
405 0 : P_TMP = ( HcoState%Grid%PEDGE%Val(I,J,L) + &
406 0 : HcoState%Grid%PEDGE%Val(I,J,L+1) ) / 200.0_hp
407 :
408 : ! Interpolate 7Be [stars/g air/sec] to GEOS-Chem levels
409 : CALL SLQ( Inst%LATSOU, Inst%PRESOU, Inst%BESOU, 10, 33, &
410 0 : LAT_TMP, P_TMP, Be_TMP )
411 :
412 : ! Be_TMP = [stars/g air/s] * [0.045 atom/star] *
413 : ! [kg air] * [1e3 g/kg] = 7Be/10Be emissions [atoms/s]
414 0 : Be_TMP = Be_TMP * 0.045e+0_hp * ExtState%AIR%Arr%Val(I,J,L) * 1.e+3_hp
415 :
416 : ! ADD_Be = [atoms/s] / [atom/kg] / [m2] = 7Be/10Be emissions [kg/m2/s]
417 0 : ADD_Be7 = ( Be_TMP / XNUMOL_Be7 ) / HcoState%Grid%AREA_M2%Val(I,J)
418 0 : ADD_Be10 = ( Be_TMP / XNUMOL_Be10 ) / HcoState%Grid%AREA_M2%Val(I,J)
419 :
420 : ! Save emissions into an array for use below
421 0 : Inst%EmissBe7 (I,J,L) = ADD_Be7
422 0 : Inst%EmissBe10(I,J,L) = ADD_Be10
423 0 : IF ( L > ExtState%TropLev%Arr%Val(I,J) ) THEN
424 0 : IF ( Inst%IDTBe7s > 0 ) THEN
425 0 : Inst%EmissBe7s (I,J,L) = Add_Be7
426 : ENDIF
427 0 : IF ( Inst%IDTBe10s > 0 ) THEN
428 0 : Inst%EmissBe10s(I,J,L) = Add_Be10
429 : ENDIF
430 : ELSE
431 0 : IF ( Inst%IDTBe7s > 0 ) THEN
432 0 : Inst%EmissBe7s (I,J,L) = 0d0
433 : ENDIF
434 0 : IF ( Inst%IDTBe10s > 0 ) THEN
435 0 : Inst%EmissBe10s(I,J,L) = 0d0
436 : ENDIF
437 : ENDIF
438 :
439 : ENDDO
440 : ENDDO
441 : ENDDO
442 : !$OMP END PARALLEL DO
443 :
444 : !------------------------------------------------------------------------
445 : ! Add Be7 and Be10 emissions to HEMCO data structure & diagnostics
446 : !------------------------------------------------------------------------
447 :
448 : ! Add emissions
449 0 : IF ( Inst%IDTBe7 > 0 ) THEN
450 0 : Arr3D => Inst%EmissBe7(:,:,:)
451 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTBe7, &
452 0 : RC, ExtNr=Inst%ExtNr )
453 0 : Arr3D => NULL()
454 0 : IF ( RC /= HCO_SUCCESS ) THEN
455 : CALL HCO_ERROR( &
456 0 : 'HCO_EmisAdd error: EmissBe7', RC )
457 0 : RETURN
458 : ENDIF
459 : ENDIF
460 :
461 : ! Add emissions
462 0 : IF ( Inst%IDTBe7s > 0 ) THEN
463 0 : Arr3D => Inst%EmissBe7s(:,:,:)
464 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTBe7s, &
465 0 : RC, ExtNr=Inst%ExtNr )
466 0 : Arr3D => NULL()
467 0 : IF ( RC /= HCO_SUCCESS ) THEN
468 : CALL HCO_ERROR( &
469 0 : 'HCO_EmisAdd error: EmissBe7s', RC )
470 0 : RETURN
471 : ENDIF
472 : ENDIF
473 :
474 : ! Add emissions
475 0 : IF ( Inst%IDTBe10 > 0 ) THEN
476 0 : Arr3D => Inst%EmissBe10(:,:,:)
477 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTBe10, &
478 0 : RC, ExtNr=Inst%ExtNr )
479 0 : Arr3D => NULL()
480 0 : IF ( RC /= HCO_SUCCESS ) THEN
481 : CALL HCO_ERROR( &
482 0 : 'HCO_EmisAdd error: EmissBe10', RC )
483 0 : RETURN
484 : ENDIF
485 : ENDIF
486 :
487 : ! Add emissions
488 0 : IF ( Inst%IDTBe10s > 0 ) THEN
489 0 : Arr3D => Inst%EmissBe10s(:,:,:)
490 : CALL HCO_EmisAdd( HcoState, Arr3D, Inst%IDTBe10s, &
491 0 : RC, ExtNr=Inst%ExtNr )
492 0 : Arr3D => NULL()
493 0 : IF ( RC /= HCO_SUCCESS ) THEN
494 : CALL HCO_ERROR( &
495 0 : 'HCO_EmisAdd error: EmissBe10s', RC )
496 0 : RETURN
497 : ENDIF
498 : ENDIF
499 :
500 : ENDIF !IDTBe7 > 0 or IDTBe10 > 0
501 :
502 : !=======================================================================
503 : ! Cleanup & quit
504 : !=======================================================================
505 :
506 : ! Nullify pointers
507 0 : Inst => NULL()
508 :
509 : ! Return w/ success
510 0 : CALL HCO_LEAVE( HcoState%Config%Err,RC )
511 :
512 0 : END SUBROUTINE HCOX_Gc_RnPbBe_Run
513 : !EOC
514 : !------------------------------------------------------------------------------
515 : ! Harmonized Emissions Component (HEMCO) !
516 : !------------------------------------------------------------------------------
517 : !BOP
518 : !
519 : ! !IROUTINE: HCOX_Gc_RnPbBe_Init
520 : !
521 : ! !DESCRIPTION: Subroutine HcoX\_Gc\_RnPbBe\_Init initializes the HEMCO
522 : ! GC\_Rn-Pb-Be extension.
523 : !\\
524 : !\\
525 : ! !INTERFACE:
526 : !
527 0 : SUBROUTINE HCOX_Gc_RnPbBe_Init( HcoState, ExtName, ExtState, RC )
528 : !
529 : ! !USES:
530 : !
531 : USE HCO_ExtList_Mod, ONLY : GetExtNr
532 : USE HCO_ExtList_Mod, ONLY : GetExtOpt
533 : USE HCO_State_Mod, ONLY : HCO_GetExtHcoID
534 : !
535 : ! !INPUT PARAMETERS:
536 : !
537 : CHARACTER(LEN=*), INTENT(IN ) :: ExtName ! Extension name
538 : TYPE(Ext_State), POINTER :: ExtState ! Module options
539 : !
540 : ! !INPUT/OUTPUT PARAMETERS:
541 : !
542 : TYPE(HCO_State), POINTER :: HcoState ! Hemco state
543 : INTEGER, INTENT(INOUT) :: RC
544 :
545 : ! !REVISION HISTORY:
546 : ! 07 Jul 2014 - R. Yantosca - Initial version
547 : ! See https://github.com/geoschem/hemco for complete history
548 : !EOP
549 : !------------------------------------------------------------------------------
550 : !BOC
551 : !
552 : ! !LOCAL VARIABLES:
553 : !
554 : ! Scalars
555 : INTEGER :: N, nSpc, ExtNr, ExtNrZhang
556 : CHARACTER(LEN=255) :: MSG, LOC
557 :
558 : ! Arrays
559 0 : INTEGER, ALLOCATABLE :: HcoIDs(:)
560 0 : CHARACTER(LEN=31), ALLOCATABLE :: SpcNames(:)
561 :
562 : ! Pointers
563 : TYPE(MyInst), POINTER :: Inst
564 :
565 : !=======================================================================
566 : ! HCOX_GC_RnPbBe_INIT begins here!
567 : !=======================================================================
568 0 : LOC = 'HCOX_GC_RNPBBE_INIT (HCOX_GC_RNPBBE_MOD.F90)'
569 :
570 : ! Get the main extension number
571 0 : ExtNr = GetExtNr( HcoState%Config%ExtList, TRIM(ExtName) )
572 0 : IF ( ExtNr <= 0 ) RETURN
573 :
574 : ! Get the extension number for Zhang et al [2021] emissions
575 0 : ExtNrZhang = GetExtNr( HcoState%Config%ExtList, 'ZHANG_Rn222' )
576 :
577 : ! Enter
578 0 : CALL HCO_ENTER( HcoState%Config%Err, LOC, RC )
579 0 : IF ( RC /= HCO_SUCCESS ) THEN
580 0 : CALL HCO_ERROR( 'ERROR 1', RC, THISLOC=LOC )
581 0 : RETURN
582 : ENDIF
583 :
584 : ! Create Instance
585 0 : Inst => NULL()
586 0 : CALL InstCreate ( ExtNr, ExtState%GC_RnPbBe, Inst, RC )
587 0 : IF ( RC /= HCO_SUCCESS ) THEN
588 0 : CALL HCO_ERROR ( 'Cannot create GC_RnPbBe instance', RC )
589 0 : RETURN
590 : ENDIF
591 : ! Also fill the extension numbers in the Instance object
592 0 : Inst%ExtNr = ExtNr
593 0 : Inst%ExtNrZhang = ExtNrZhang
594 :
595 : ! Set HEMCO species IDs
596 0 : CALL HCO_GetExtHcoID( HcoState, Inst%ExtNr, HcoIDs, SpcNames, nSpc, RC )
597 0 : IF ( RC /= HCO_SUCCESS ) THEN
598 0 : CALL HCO_ERROR( 'Could not set HEMCO species IDs', RC )
599 0 : RETURN
600 : ENDIF
601 :
602 : ! Verbose mode
603 0 : IF ( HcoState%amIRoot ) THEN
604 :
605 : ! Write the name of the extension regardless of the verbose setting
606 0 : msg = 'Using HEMCO extension: GC_RnPbBe (radionuclide emissions)'
607 0 : IF ( HCO_IsVerb( HcoState%Config%Err ) ) THEN
608 0 : CALL HCO_Msg( HcoState%Config%Err, sep1='-' ) ! with separator
609 : ELSE
610 0 : CALL HCO_Msg( msg, verb=.TRUE. ) ! w/o separator
611 : ENDIF
612 :
613 : ! Write all other messages as debug printout only
614 0 : MSG = 'Use the following species (Name: HcoID):'
615 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
616 0 : DO N = 1, nSpc
617 0 : WRITE(MSG,*) TRIM(SpcNames(N)), ':', HcoIDs(N)
618 0 : CALL HCO_MSG(HcoState%Config%Err,MSG)
619 : ENDDO
620 : ENDIF
621 :
622 : ! Set up tracer and HEMCO indices
623 0 : DO N = 1, nSpc
624 0 : SELECT CASE( TRIM( SpcNames(N) ) )
625 : CASE( 'Rn', 'Rn222', '222Rn' )
626 0 : Inst%IDTRn222 = HcoIDs(N)
627 : CASE( 'Be', 'Be7', '7Be' )
628 0 : Inst%IDTBe7 = HcoIDs(N)
629 : CASE( 'Be7s', '7Bes' )
630 0 : Inst%IDTBe7s = HcoIDs(N)
631 : CASE( 'Be10', '10Be' )
632 0 : Inst%IDTBe10 = HcoIDs(N)
633 : CASE( 'Be10s', '10Bes' )
634 0 : Inst%IDTBe10s = HcoIDs(N)
635 : CASE DEFAULT
636 : ! Do nothing
637 : END SELECT
638 : ENDDO
639 :
640 : ! WARNING: Rn tracer is not found!
641 0 : IF ( Inst%IDTRn222 <= 0 .AND. HcoState%amIRoot ) THEN
642 : CALL HCO_WARNING( HcoState%Config%Err, &
643 0 : 'Cannot find Rn222 tracer in list of species!', RC )
644 : ENDIF
645 :
646 : ! WARNING: Be7 tracer is not found
647 0 : IF ( Inst%IDTBe7 <= 0 .AND. HcoState%amIRoot ) THEN
648 : CALL HCO_WARNING( HcoState%Config%Err, &
649 0 : 'Cannot find Be7 tracer in list of species!', RC )
650 : ENDIF
651 :
652 : ! WARNING: Be10 tracer is not found
653 0 : IF ( Inst%IDTBe10 <= 0 .AND. HcoState%amIRoot ) THEN
654 : CALL HCO_WARNING( HcoState%Config%Err, &
655 0 : 'Cannot find Be10 tracer in list of species!', RC )
656 : ENDIF
657 :
658 : ! ERROR: No tracer defined
659 0 : IF ( Inst%IDTRn222 <= 0 .AND. Inst%IDTBe7 <= 0 .AND. Inst%IDTBe10 <= 0) THEN
660 : CALL HCO_ERROR( &
661 0 : 'Cannot use RnPbBe extension: no valid species!', RC )
662 : ENDIF
663 :
664 : ! Activate met fields required by this extension
665 0 : ExtState%FRCLND%DoUse = .TRUE.
666 0 : ExtState%T2M%DoUse = .TRUE.
667 0 : ExtState%AIR%DoUse = .TRUE.
668 0 : ExtState%TropLev%DoUse = .TRUE.
669 :
670 : !=======================================================================
671 : ! Initialize data arrays
672 : !=======================================================================
673 :
674 0 : IF ( Inst%IDTRn222 > 0 ) THEN
675 0 : ALLOCATE( Inst%EmissRn222( HcoState%Nx, HcoState%NY ), STAT=RC )
676 0 : IF ( RC /= 0 ) THEN
677 : CALL HCO_ERROR ( &
678 0 : 'Cannot allocate EmissRn222', RC )
679 0 : RETURN
680 : ENDIF
681 : ENDIF
682 :
683 0 : IF ( Inst%IDTBe7 > 0 ) THEN
684 : ALLOCATE( Inst%EmissBe7( HcoState%Nx, HcoState%NY, HcoState%NZ ), &
685 0 : STAT=RC )
686 0 : IF ( RC /= 0 ) THEN
687 : CALL HCO_ERROR ( &
688 0 : 'Cannot allocate EmissBe7', RC )
689 0 : RETURN
690 : ENDIF
691 : IF ( RC /= 0 ) RETURN
692 :
693 : ! Array for latitudes (Lal & Peters data)
694 0 : ALLOCATE( Inst%LATSOU( 10 ), STAT=RC )
695 0 : IF ( RC /= 0 ) THEN
696 : CALL HCO_ERROR ( &
697 0 : 'Cannot allocate LATSOU', RC )
698 0 : RETURN
699 : ENDIF
700 :
701 : ! Array for pressures (Lal & Peters data)
702 0 : ALLOCATE( Inst%PRESOU( 33 ), STAT=RC )
703 0 : IF ( RC /= 0 ) THEN
704 : CALL HCO_ERROR ( &
705 0 : 'Cannot allocate PRESOU', RC )
706 0 : RETURN
707 : ENDIF
708 :
709 : ! Array for 7Be emissions ( Lal & Peters data)
710 0 : ALLOCATE( Inst%BESOU( 10, 33 ), STAT=RC )
711 0 : IF ( RC /= 0 ) THEN
712 : CALL HCO_ERROR ( &
713 0 : 'Cannot allocate BESOU', RC )
714 0 : RETURN
715 : ENDIF
716 :
717 : ! Initialize the 7Be emisisons data arrays
718 0 : CALL Init_7Be_Emissions( Inst )
719 : ENDIF
720 :
721 0 : IF ( Inst%IDTBe7s > 0 ) THEN
722 : ALLOCATE( Inst%EmissBe7s( HcoState%Nx, HcoState%NY, HcoState%NZ ), &
723 0 : STAT=RC )
724 0 : IF ( RC /= 0 ) THEN
725 : CALL HCO_ERROR ( &
726 0 : 'Cannot allocate EmissBe7s', RC )
727 0 : RETURN
728 : ENDIF
729 0 : Inst%EmissBe7s = 0.0_hp
730 : ENDIF
731 :
732 0 : IF ( Inst%IDTBe10 > 0 ) THEN
733 : ALLOCATE( Inst%EmissBe10( HcoState%Nx, HcoState%NY, HcoState%NZ ), &
734 0 : STAT=RC )
735 0 : IF ( RC /= 0 ) THEN
736 : CALL HCO_ERROR ( &
737 0 : 'Cannot allocate EmissBe10', RC )
738 0 : RETURN
739 : ENDIF
740 : ENDIF
741 :
742 0 : IF ( Inst%IDTBe10s > 0 ) THEN
743 : ALLOCATE( Inst%EmissBe10s( HcoState%Nx, HcoState%NY, HcoState%NZ ), &
744 0 : STAT=RC )
745 0 : IF ( RC /= 0 ) THEN
746 : CALL HCO_ERROR ( &
747 0 : 'Cannot allocate EmissBe10s', RC )
748 0 : RETURN
749 : ENDIF
750 0 : Inst%EmissBe10s = 0.0_hp
751 : ENDIF
752 :
753 : !=======================================================================
754 : ! Leave w/ success
755 : !=======================================================================
756 0 : IF ( ALLOCATED( HcoIDs ) ) DEALLOCATE( HcoIDs )
757 0 : IF ( ALLOCATED( SpcNames ) ) DEALLOCATE( SpcNames )
758 :
759 : ! Nullify pointers
760 0 : Inst => NULL()
761 :
762 0 : CALL HCO_LEAVE( HcoState%Config%Err,RC )
763 :
764 0 : END SUBROUTINE HCOX_Gc_RnPbBe_Init
765 : !EOC
766 : !------------------------------------------------------------------------------
767 : ! Harmonized Emissions Component (HEMCO) !
768 : !------------------------------------------------------------------------------
769 : !BOP
770 : !
771 : ! !IROUTINE: HCOX_Gc_RnPbBe_Final
772 : !
773 : ! !DESCRIPTION: Subroutine HcoX\_Gc\_RnPbBe\_Final finalizes the HEMCO
774 : ! extension for the GEOS-Chem Rn-Pb-Be specialty simulation. All module
775 : ! arrays will be deallocated.
776 : !\\
777 : !\\
778 : ! !INTERFACE:
779 : !
780 0 : SUBROUTINE HCOX_Gc_RnPbBe_Final( ExtState )
781 : !
782 : ! !INPUT PARAMETERS:
783 : !
784 : TYPE(Ext_State), POINTER :: ExtState ! Module options
785 : !
786 : ! !REVISION HISTORY:
787 : ! 13 Dec 2013 - C. Keller - Now a HEMCO extension
788 : ! See https://github.com/geoschem/hemco for complete history
789 : !EOP
790 : !------------------------------------------------------------------------------
791 : !BOC
792 :
793 : !=======================================================================
794 : ! HCOX_GC_RNPBBE_FINAL begins here!
795 : !=======================================================================
796 :
797 0 : CALL InstRemove ( ExtState%GC_RnPbBe )
798 :
799 0 : END SUBROUTINE HCOX_Gc_RnPbBe_Final
800 : !EOC
801 : !------------------------------------------------------------------------------
802 : ! Harmonized Emissions Component (HEMCO) !
803 : !------------------------------------------------------------------------------
804 : !BOP
805 : !
806 : ! !IROUTINE: Init_7Be_Emissions
807 : !
808 : ! !DESCRIPTION: Subroutine Init\_7Be\_Emissions initializes the 7Be emissions
809 : ! from Lal \& Peters on 33 pressure levels. This data used to be read from
810 : ! a file, but we have now hardwired it to facilitate I/O in the ESMF
811 : ! environment.
812 : !\\
813 : !\\
814 : ! !INTERFACE:
815 : !
816 0 : SUBROUTINE Init_7Be_Emissions( Inst )
817 : !
818 : ! !INPUT PARAMETERS:
819 : !
820 : TYPE(MyInst), POINTER :: Inst ! Instance
821 : !
822 : ! !REMARKS:
823 : ! (1) Reference: Lal, D., and B. Peters, Cosmic ray produced radioactivity
824 : ! on the Earth. Handbuch der Physik, 46/2, 551-612, edited by K. Sitte,
825 : ! Springer-Verlag, New York, 1967.
826 : ! .
827 : ! (2) In prior versions of GEOS-Chem, this routine was named READ_7BE, and
828 : ! it read the ASCII file "7Be.Lal". Because this data set is not placed
829 : ! on a lat/lon grid, ESMF cannot regrid it. To work around this, we now
830 : ! hardwire this data in module arrays rather than read it from disk.
831 : ! .
832 : ! (3) Units of 7Be emissions are [stars/g air/s].
833 : ! Here, "stars" = # of nuclear disintegrations of cosmic rays
834 : ! .
835 : ! (4) Original data from Lal & Peters (1967), w/ these modifications:
836 : ! (a) Replace data at (0hPa, 70S) following Koch 1996:
837 : ! (i ) old value = 3000
838 : ! (ii) new value = 1900
839 : ! (b) Copy data from 70S to 80S and 90S at all levels
840 : ! .
841 : ! !REVISION HISTORY:
842 : ! 07 Aug 2002 - H. Liu - Initial version
843 : ! See https://github.com/geoschem/hemco for complete history
844 : !EOP
845 : !------------------------------------------------------------------------------
846 : !BOC
847 :
848 : ! Define latitudes [degrees North]
849 : Inst%LATSOU = (/ 0.0_hp, 10.0_hp, 20.0_hp, 30.0_hp, &
850 : 40.0_hp, 50.0_hp, 60.0_hp, 70.0_hp, &
851 0 : 80.0_hp, 90.0_hp /)
852 :
853 : ! Define pressures [hPa]
854 : Inst%PRESOU = (/ 0.0_hp, 50.0_hp, 70.0_hp, 90.0_hp, &
855 : 110.0_hp, 130.0_hp, 150.0_hp, 170.0_hp, &
856 : 190.0_hp, 210.0_hp, 230.0_hp, 250.0_hp, &
857 : 270.0_hp, 290.0_hp, 313.0_hp, 338.0_hp, &
858 : 364.0_hp, 392.0_hp, 420.0_hp, 451.0_hp, &
859 : 485.0_hp, 518.0_hp, 555.0_hp, 592.0_hp, &
860 : 633.0_hp, 680.0_hp, 725.0_hp, 772.0_hp, &
861 : 822.0_hp, 875.0_hp, 930.0_hp, 985.0_hp, &
862 0 : 1030.0_hp /)
863 :
864 : ! Define 7Be emissions [stars/g air/s]
865 : ! 1 "star" = 1 nuclear disintegration via cosmic rays
866 : !
867 : ! NOTE: These statements were defined from printout of the file
868 : ! and need to be multiplied by 1d-5 below.
869 0 : Inst%BESOU(:,1) = (/ 150.0_hp, 156.0_hp, 188.0_hp, 285.0_hp, &
870 : 500.0_hp, 910.0_hp, 1700.0_hp, 1900.0_hp, &
871 0 : 1900.0_hp, 1900.0_hp /)
872 :
873 0 : Inst%BESOU(:,2) = (/ 280.0_hp, 310.0_hp, 390.0_hp, 590.0_hp, &
874 : 880.0_hp, 1390.0_hp, 1800.0_hp, 1800.0_hp, &
875 0 : 1800.0_hp, 1800.0_hp /)
876 :
877 0 : Inst%BESOU(:,3) = (/ 310.0_hp, 330.0_hp, 400.0_hp, 620.0_hp, &
878 : 880.0_hp, 1280.0_hp, 1450.0_hp, 1450.0_hp, &
879 0 : 1450.0_hp, 1450.0_hp /)
880 :
881 0 : Inst%BESOU(:,4) = (/ 285.0_hp, 310.0_hp, 375.0_hp, 570.0_hp, &
882 : 780.0_hp, 1100.0_hp, 1180.0_hp, 1180.0_hp, &
883 0 : 1180.0_hp, 1180.0_hp /)
884 :
885 0 : Inst%BESOU(:,5) = (/ 255.0_hp, 275.0_hp, 330.0_hp, 510.0_hp, &
886 : 680.0_hp, 950.0_hp, 1000.0_hp, 1000.0_hp, &
887 0 : 1000.0_hp, 1000.0_hp /)
888 :
889 0 : Inst%BESOU(:,6) = (/ 230.0_hp, 245.0_hp, 292.0_hp, 450.0_hp, &
890 : 600.0_hp, 820.0_hp, 875.0_hp, 875.0_hp, &
891 0 : 875.0_hp, 875.0_hp /)
892 :
893 0 : Inst%BESOU(:,7) = (/ 205.0_hp, 215.0_hp, 260.0_hp, 400.0_hp, &
894 : 530.0_hp, 730.0_hp, 750.0_hp, 750.0_hp, &
895 0 : 750.0_hp, 750.0_hp /)
896 :
897 0 : Inst%BESOU(:,8) = (/ 182.0_hp, 195.0_hp, 235.0_hp, 355.0_hp, &
898 : 480.0_hp, 630.0_hp, 650.0_hp, 650.0_hp, &
899 0 : 650.0_hp, 650.0_hp /)
900 :
901 0 : Inst%BESOU(:,9) = (/ 160.0_hp, 173.0_hp, 208.0_hp, 315.0_hp, &
902 : 410.0_hp, 543.0_hp, 550.0_hp, 550.0_hp, &
903 0 : 550.0_hp, 550.0_hp /)
904 :
905 0 : Inst%BESOU(:,10) = (/ 148.0_hp, 152.0_hp, 185.0_hp, 280.0_hp, &
906 : 370.0_hp, 480.0_hp, 500.0_hp, 500.0_hp, &
907 0 : 500.0_hp, 500.0_hp /)
908 :
909 0 : Inst%BESOU(:,11) = (/ 130.0_hp, 139.0_hp, 167.0_hp, 250.0_hp, &
910 : 320.0_hp, 425.0_hp, 430.0_hp, 430.0_hp, &
911 0 : 430.0_hp, 430.0_hp /)
912 :
913 0 : Inst%BESOU(:,12) = (/ 116.0_hp, 123.0_hp, 148.0_hp, 215.0_hp, &
914 : 285.0_hp, 365.0_hp, 375.0_hp, 375.0_hp, &
915 0 : 375.0_hp, 375.0_hp /)
916 :
917 0 : Inst%BESOU(:,13) = (/ 104.0_hp, 110.0_hp, 130.0_hp, 198.0_hp, &
918 : 250.0_hp, 320.0_hp, 330.0_hp, 330.0_hp, &
919 0 : 330.0_hp, 330.0_hp /)
920 :
921 0 : Inst%BESOU(:,14) = (/ 93.0_hp, 99.0_hp, 118.0_hp, 170.0_hp, &
922 : 222.0_hp, 280.0_hp, 288.0_hp, 288.0_hp, &
923 0 : 288.0_hp, 288.0_hp /)
924 :
925 0 : Inst%BESOU(:,15) = (/ 80.0_hp, 84.0_hp, 100.0_hp, 145.0_hp, &
926 : 190.0_hp, 235.0_hp, 250.0_hp, 250.0_hp, &
927 0 : 250.0_hp, 250.0_hp /)
928 :
929 0 : Inst%BESOU(:,16) = (/ 72.0_hp, 74.0_hp, 88.0_hp, 129.0_hp, &
930 : 168.0_hp, 210.0_hp, 218.0_hp, 218.0_hp, &
931 0 : 218.0_hp, 218.0_hp /)
932 :
933 0 : Inst%BESOU(:,17) = (/ 59.5_hp, 62.5_hp, 73.5_hp, 108.0_hp, &
934 : 138.0_hp, 171.0_hp, 178.0_hp, 178.0_hp, &
935 0 : 178.0_hp, 178.0_hp /)
936 :
937 0 : Inst%BESOU(:,18) = (/ 50.0_hp, 53.0_hp, 64.0_hp, 90.0_hp, &
938 : 115.0_hp, 148.0_hp, 150.0_hp, 150.0_hp, &
939 0 : 150.0_hp, 150.0_hp /)
940 :
941 0 : Inst%BESOU(:,19) = (/ 45.0_hp, 46.5_hp, 52.5_hp, 76.0_hp, &
942 : 98.0_hp, 122.0_hp, 128.0_hp, 128.0_hp, &
943 0 : 128.0_hp, 128.0_hp /)
944 :
945 0 : Inst%BESOU(:,20) = (/ 36.5_hp, 37.5_hp, 45.0_hp, 61.0_hp, &
946 : 77.0_hp, 98.0_hp, 102.0_hp, 102.0_hp, &
947 0 : 102.0_hp, 102.0_hp /)
948 :
949 0 : Inst%BESOU(:,21) = (/ 30.8_hp, 32.0_hp, 37.5_hp, 51.5_hp, &
950 : 65.0_hp, 81.0_hp, 85.0_hp, 85.0_hp, &
951 0 : 85.0_hp, 85.0_hp /)
952 :
953 0 : Inst%BESOU(:,22) = (/ 25.5_hp, 26.5_hp, 32.0_hp, 40.5_hp, &
954 : 54.0_hp, 67.5_hp, 69.5_hp, 69.5_hp, &
955 0 : 69.5_hp, 69.5_hp /)
956 :
957 0 : Inst%BESOU(:,23) = (/ 20.5_hp, 21.6_hp, 25.5_hp, 33.0_hp, &
958 : 42.0_hp, 53.5_hp, 55.0_hp, 55.0_hp, &
959 0 : 55.0_hp, 55.0_hp /)
960 :
961 0 : Inst%BESOU(:,24) = (/ 16.8_hp, 17.3_hp, 20.0_hp, 26.0_hp, &
962 : 33.5_hp, 41.0_hp, 43.0_hp, 43.0_hp, &
963 0 : 43.0_hp, 43.0_hp /)
964 :
965 0 : Inst%BESOU(:,25) = (/ 13.0_hp, 13.8_hp, 15.3_hp, 20.5_hp, &
966 : 26.8_hp, 32.5_hp, 33.5_hp, 33.5_hp, &
967 0 : 33.5_hp, 33.5_hp /)
968 :
969 0 : Inst%BESOU(:,26) = (/ 10.1_hp, 10.6_hp, 12.6_hp, 15.8_hp, &
970 : 20.0_hp, 24.5_hp, 25.8_hp, 25.8_hp, &
971 0 : 25.8_hp, 25.8_hp /)
972 :
973 0 : Inst%BESOU(:,27) = (/ 7.7_hp, 8.15_hp, 9.4_hp, 11.6_hp, &
974 : 14.8_hp, 17.8_hp, 18.5_hp, 18.5_hp, &
975 0 : 18.5_hp, 18.5_hp /)
976 :
977 0 : Inst%BESOU(:,28) = (/ 5.7_hp, 5.85_hp, 6.85_hp, 8.22_hp, &
978 : 11.0_hp, 13.1_hp, 13.2_hp, 13.2_hp, &
979 0 : 13.2_hp, 13.2_hp /)
980 :
981 0 : Inst%BESOU(:,29) = (/ 3.9_hp, 4.2_hp, 4.85_hp, 6.0_hp, &
982 : 7.6_hp, 9.0_hp, 9.2_hp, 9.2_hp, &
983 0 : 9.2_hp, 9.2_hp /)
984 :
985 0 : Inst%BESOU(:,30) = (/ 3.0_hp, 3.05_hp, 3.35_hp, 4.2_hp, &
986 : 5.3_hp, 5.9_hp, 6.25_hp, 6.25_hp, &
987 0 : 6.25_hp, 6.25_hp /)
988 :
989 0 : Inst%BESOU(:,31) = (/ 2.05_hp, 2.1_hp, 2.32_hp, 2.9_hp, &
990 : 3.4_hp, 3.9_hp, 4.1_hp, 4.1_hp, &
991 0 : 4.1_hp, 4.1_hp /)
992 :
993 0 : Inst%BESOU(:,32) = (/ 1.45_hp, 1.43_hp, 1.65_hp, 2.03_hp, &
994 : 2.4_hp, 2.75_hp, 2.65_hp, 2.65_hp, &
995 0 : 2.65_hp, 2.65_hp /)
996 :
997 0 : Inst%BESOU(:,33) = (/ 1.04_hp, 1.08_hp, 1.21_hp, 1.5_hp, &
998 : 1.68_hp, 1.8_hp, 1.8_hp, 1.8_hp, &
999 0 : 1.8_hp, 1.8_hp /)
1000 :
1001 : ! All the numbers of BESOU need to be multiplied by 1e-5 in order to put
1002 : ! them into the correct data range. NOTE: This multiplication statement
1003 : ! needs to be preserved here in order to ensure identical output to the
1004 : ! prior code! (bmy, 7/7/14)
1005 0 : Inst%BESOU = Inst%BESOU * 1.e-5_hp
1006 :
1007 0 : END SUBROUTINE Init_7Be_Emissions
1008 : !EOC
1009 : !------------------------------------------------------------------------------
1010 : ! Harmonized Emissions Component (HEMCO) !
1011 : !------------------------------------------------------------------------------
1012 : !BOP
1013 : !
1014 : ! !IROUTINE: SLQ
1015 : !
1016 : ! !DESCRIPTION: Subroutine SLQ is an interpolation subroutine from a
1017 : ! Chinese reference book (says Hongyu Liu).
1018 : !\\
1019 : !\\
1020 : ! !INTERFACE:
1021 : !
1022 0 : SUBROUTINE SLQ( X, Y, Z, N, M, U, V, W )
1023 : !
1024 : ! !INPUT PARAMETERS:
1025 : !
1026 : INTEGER :: N ! First dimension of Z
1027 : INTEGER :: M ! Second dimension of Z
1028 : REAL(hp) :: X(N) ! X-axis coordinate on original grid
1029 : REAL(hp) :: Y(M) ! Y-axis coordinate on original grid
1030 : REAL(hp) :: Z(N,M) ! Array of data on original grid
1031 : REAL(hp) :: U ! X-axis coordinate for desired interpolated value
1032 : REAL(hp) :: V ! Y-axis coordinate for desired interpolated value
1033 : !
1034 : ! !OUTPUT PARAMETERS:
1035 : !
1036 : REAL(hp) :: W ! Interpolated value of Z array, at coords (U,V)
1037 : !
1038 : ! !REMARKS:
1039 : ! This routine was taken from the old RnPbBe_mod.F.
1040 : !
1041 : ! !REVISION HISTORY:
1042 : ! 17 Mar 1998 - H. Liu - Initial version
1043 : ! See https://github.com/geoschem/hemco for complete history
1044 : !EOP
1045 : !------------------------------------------------------------------------------
1046 : !BOC
1047 : !
1048 : ! !LOCAL VARIABLES:
1049 : !
1050 : REAL(hp) :: B(3), HH
1051 : INTEGER :: NN, IP, I, J, L, IQ, K, MM
1052 :
1053 : !=======================================================================
1054 : ! SLQ begins here!
1055 : !=======================================================================
1056 0 : NN=3
1057 0 : IF(N.LE.3) THEN
1058 : IP=1
1059 : NN=N
1060 0 : ELSE IF (U.LE.X(2)) THEN
1061 : IP=1
1062 0 : ELSE IF (U.GE.X(N-1)) THEN
1063 0 : IP=N-2
1064 : ELSE
1065 : I=1
1066 : J=N
1067 0 : 10 IF (IABS(I-J).NE.1) THEN
1068 0 : L=(I+J)/2
1069 0 : IF (U.LT.X(L)) THEN
1070 : J=L
1071 : ELSE
1072 0 : I=L
1073 : END IF
1074 : GOTO 10
1075 : END IF
1076 0 : IF (ABS(U-X(I)).LT.ABS(U-X(J))) THEN
1077 0 : IP=I-1
1078 : ELSE
1079 : IP=I
1080 : END IF
1081 : END IF
1082 0 : MM=3
1083 0 : IF (M.LE.3) THEN
1084 : IQ=1
1085 : MM=N
1086 0 : ELSE IF (V.LE.Y(2)) THEN
1087 : IQ=1
1088 0 : ELSE IF (V.GE.Y(M-1)) THEN
1089 0 : IQ=M-2
1090 : ELSE
1091 : I=1
1092 : J=M
1093 0 : 20 IF (IABS(J-I).NE.1) THEN
1094 0 : L=(I+J)/2
1095 0 : IF (V.LT.Y(L)) THEN
1096 : J=L
1097 : ELSE
1098 0 : I=L
1099 : END IF
1100 : GOTO 20
1101 : END IF
1102 0 : IF (ABS(V-Y(I)).LT.ABS(V-Y(J))) THEN
1103 0 : IQ=I-1
1104 : ELSE
1105 : IQ=I
1106 : END IF
1107 : END IF
1108 0 : DO 50 I=1,NN
1109 0 : B(I)=0.0
1110 0 : DO 40 J=1,MM
1111 0 : HH=Z(IP+I-1,IQ+J-1)
1112 0 : DO 30 K=1,MM
1113 0 : IF (K.NE.J) THEN
1114 0 : HH=HH*(V-Y(IQ+K-1))/(Y(IQ+J-1)-Y(IQ+K-1))
1115 : END IF
1116 0 : 30 CONTINUE
1117 0 : B(I)=B(I)+HH
1118 0 : 40 CONTINUE
1119 0 : 50 CONTINUE
1120 0 : W=0.0
1121 0 : DO 70 I=1,NN
1122 0 : HH=B(I)
1123 0 : DO 60 J=1,NN
1124 0 : IF (J.NE.I) THEN
1125 0 : HH=HH*(U-X(IP+J-1))/(X(IP+I-1)-X(IP+J-1))
1126 : END IF
1127 0 : 60 CONTINUE
1128 0 : W=W+HH
1129 0 : 70 CONTINUE
1130 :
1131 0 : END SUBROUTINE SLQ
1132 : !EOC
1133 : !------------------------------------------------------------------------------
1134 : ! Harmonized Emissions Component (HEMCO) !
1135 : !------------------------------------------------------------------------------
1136 : !BOP
1137 : !
1138 : ! !IROUTINE: InstGet
1139 : !
1140 : ! !DESCRIPTION: Subroutine InstGet returns a poiner to the desired instance.
1141 : !\\
1142 : !\\
1143 : ! !INTERFACE:
1144 : !
1145 0 : SUBROUTINE InstGet ( Instance, Inst, RC, PrevInst )
1146 : !
1147 : ! !INPUT PARAMETERS:
1148 : !
1149 : INTEGER :: Instance
1150 : TYPE(MyInst), POINTER :: Inst
1151 : INTEGER :: RC
1152 : TYPE(MyInst), POINTER, OPTIONAL :: PrevInst
1153 : !
1154 : ! !REVISION HISTORY:
1155 : ! 18 Feb 2016 - C. Keller - Initial version
1156 : ! See https://github.com/geoschem/hemco for complete history
1157 : !EOP
1158 : !------------------------------------------------------------------------------
1159 : !BOC
1160 : TYPE(MyInst), POINTER :: PrvInst
1161 :
1162 : !=================================================================
1163 : ! InstGet begins here!
1164 : !=================================================================
1165 :
1166 : ! Get instance. Also archive previous instance.
1167 0 : PrvInst => NULL()
1168 0 : Inst => AllInst
1169 0 : DO WHILE ( ASSOCIATED(Inst) )
1170 0 : IF ( Inst%Instance == Instance ) EXIT
1171 0 : PrvInst => Inst
1172 0 : Inst => Inst%NextInst
1173 : END DO
1174 0 : IF ( .NOT. ASSOCIATED( Inst ) ) THEN
1175 0 : RC = HCO_FAIL
1176 0 : RETURN
1177 : ENDIF
1178 :
1179 : ! Pass output arguments
1180 0 : IF ( PRESENT(PrevInst) ) PrevInst => PrvInst
1181 :
1182 : ! Cleanup & Return
1183 0 : PrvInst => NULL()
1184 0 : RC = HCO_SUCCESS
1185 :
1186 : END SUBROUTINE InstGet
1187 : !EOC
1188 : !------------------------------------------------------------------------------
1189 : ! Harmonized Emissions Component (HEMCO) !
1190 : !------------------------------------------------------------------------------
1191 : !BOP
1192 : !
1193 : ! !IROUTINE: InstCreate
1194 : !
1195 : ! !DESCRIPTION: Subroutine InstCreate creates a new instance.
1196 : !\\
1197 : !\\
1198 : ! !INTERFACE:
1199 : !
1200 0 : SUBROUTINE InstCreate ( ExtNr, Instance, Inst, RC )
1201 : !
1202 : ! !INPUT PARAMETERS:
1203 : !
1204 : INTEGER, INTENT(IN) :: ExtNr
1205 : !
1206 : ! !OUTPUT PARAMETERS:
1207 : !
1208 : INTEGER, INTENT( OUT) :: Instance
1209 : TYPE(MyInst), POINTER :: Inst
1210 : !
1211 : ! !INPUT/OUTPUT PARAMETERS:
1212 : !
1213 : INTEGER, INTENT(INOUT) :: RC
1214 : !
1215 : ! !REVISION HISTORY:
1216 : ! 18 Feb 2016 - C. Keller - Initial version
1217 : ! See https://github.com/geoschem/hemco for complete history
1218 : !EOP
1219 : !------------------------------------------------------------------------------
1220 : !BOC
1221 : TYPE(MyInst), POINTER :: TmpInst
1222 : INTEGER :: nnInst
1223 :
1224 : !=================================================================
1225 : ! InstCreate begins here!
1226 : !=================================================================
1227 :
1228 : ! ----------------------------------------------------------------
1229 : ! Generic instance initialization
1230 : ! ----------------------------------------------------------------
1231 :
1232 : ! Initialize
1233 0 : Inst => NULL()
1234 :
1235 : ! Get number of already existing instances
1236 0 : TmpInst => AllInst
1237 0 : nnInst = 0
1238 0 : DO WHILE ( ASSOCIATED(TmpInst) )
1239 0 : nnInst = nnInst + 1
1240 0 : TmpInst => TmpInst%NextInst
1241 : END DO
1242 :
1243 : ! Create new instance
1244 0 : ALLOCATE(Inst)
1245 0 : Inst%Instance = nnInst + 1
1246 0 : Inst%ExtNr = ExtNr
1247 :
1248 : ! Attach to instance list
1249 0 : Inst%NextInst => AllInst
1250 0 : AllInst => Inst
1251 :
1252 : ! Update output instance
1253 0 : Instance = Inst%Instance
1254 :
1255 : ! ----------------------------------------------------------------
1256 : ! Type specific initialization statements follow below
1257 : ! ----------------------------------------------------------------
1258 :
1259 : ! Return w/ success
1260 0 : RC = HCO_SUCCESS
1261 :
1262 0 : END SUBROUTINE InstCreate
1263 : !EOC
1264 : !------------------------------------------------------------------------------
1265 : ! Harmonized Emissions Component (HEMCO) !
1266 : !------------------------------------------------------------------------------
1267 : !BOP
1268 : !BOP
1269 : !
1270 : ! !IROUTINE: InstRemove
1271 : !
1272 : ! !DESCRIPTION: Subroutine InstRemove creates a new instance.
1273 : !\\
1274 : !\\
1275 : ! !INTERFACE:
1276 : !
1277 0 : SUBROUTINE InstRemove ( Instance )
1278 : !
1279 : ! !INPUT PARAMETERS:
1280 : !
1281 : INTEGER :: Instance
1282 : !
1283 : ! !REVISION HISTORY:
1284 : ! 18 Feb 2016 - C. Keller - Initial version
1285 : ! See https://github.com/geoschem/hemco for complete history
1286 : !EOP
1287 : !------------------------------------------------------------------------------
1288 : !BOC
1289 : INTEGER :: RC
1290 : TYPE(MyInst), POINTER :: PrevInst
1291 : TYPE(MyInst), POINTER :: Inst
1292 :
1293 : !=================================================================
1294 : ! InstRemove begins here!
1295 : !=================================================================
1296 :
1297 : ! Init
1298 0 : PrevInst => NULL()
1299 0 : Inst => NULL()
1300 :
1301 : ! Get instance. Also archive previous instance.
1302 0 : CALL InstGet ( Instance, Inst, RC, PrevInst=PrevInst )
1303 :
1304 : ! Instance-specific deallocation
1305 0 : IF ( ASSOCIATED(Inst) ) THEN
1306 :
1307 : !---------------------------------------------------------------------
1308 : ! Deallocate fields of Inst before popping Inst off the list
1309 : ! in order to avoid memory leaks (Bob Yantosca, 17 Aug 2020)
1310 : !---------------------------------------------------------------------
1311 0 : IF ( ASSOCIATED( Inst%EmissRn222 ) ) THEN
1312 0 : DEALLOCATE( Inst%EmissRn222 )
1313 : ENDIF
1314 0 : Inst%EmissRn222 => NULL()
1315 :
1316 0 : IF ( ASSOCIATED( Inst%EmissBe7 ) ) THEN
1317 0 : DEALLOCATE( Inst%EmissBe7 )
1318 : ENDIF
1319 0 : Inst%EmissBe7 => NULL()
1320 :
1321 0 : IF ( ASSOCIATED( Inst%EmissBe7s ) ) THEN
1322 0 : DEALLOCATE( Inst%EmissBe7s )
1323 : ENDIF
1324 0 : Inst%EmissBe7s => NULL()
1325 :
1326 0 : IF ( ASSOCIATED( Inst%EmissBe10 ) ) THEN
1327 0 : DEALLOCATE(Inst%EmissBe10 )
1328 : ENDIF
1329 0 : Inst%EmissBe10 => NULL()
1330 :
1331 0 : IF ( ASSOCIATED( Inst%EmissBe10s ) ) THEN
1332 0 : DEALLOCATE( Inst%EmissBe10s )
1333 : ENDIF
1334 0 : Inst%EmissBe10s => NULL()
1335 :
1336 0 : IF ( ASSOCIATED( Inst%LATSOU ) ) THEN
1337 0 : DEALLOCATE( Inst%LATSOU )
1338 : ENDIF
1339 0 : Inst%LATSOU => NULL()
1340 :
1341 0 : IF ( ASSOCIATED( Inst%PRESOU ) ) THEN
1342 0 : DEALLOCATE(Inst%PRESOU )
1343 : ENDIF
1344 0 : Inst%PRESOU => NULL()
1345 :
1346 0 : IF ( ASSOCIATED( Inst%BESOU ) ) THEN
1347 0 : DEALLOCATE( Inst%BESOU )
1348 : ENDIF
1349 0 : Inst%BESOU => NULL()
1350 :
1351 : !---------------------------------------------------------------------
1352 : ! Pop off instance from list
1353 : !---------------------------------------------------------------------
1354 0 : IF ( ASSOCIATED(PrevInst) ) THEN
1355 0 : PrevInst%NextInst => Inst%NextInst
1356 : ELSE
1357 0 : AllInst => Inst%NextInst
1358 : ENDIF
1359 0 : DEALLOCATE(Inst)
1360 :
1361 : ENDIF
1362 :
1363 : ! Free pointers before exiting
1364 0 : PrevInst => NULL()
1365 0 : Inst => NULL()
1366 :
1367 0 : END SUBROUTINE InstRemove
1368 : !EOC
1369 0 : END MODULE HCOX_GC_RnPbBe_Mod
|
