File: C:\NOAA\NEMS_11731\src\chem\gocart\src\Components\GOCART_GridComp\BC_GridComp\BC_GridCompMod.F90

1     #ifdef GEOS5
2     #include "MAPL_Generic.h"
3     #endif
4     
5     !-------------------------------------------------------------------------
6     !         NASA/GSFC, Data Assimilation Office, Code 910.3, GEOS/DAS      !
7     !-------------------------------------------------------------------------
8     !BOP
9     !
10     ! !MODULE:  BC_GridCompMod --- BC Grid Component Class
11     !
12     ! !INTERFACE:
13     !
14     
15        module  BC_GridCompMod
16     
17     ! !USES:
18     
19     #ifdef GEOS5
20        USE ESMF_Mod
21        USE MAPL_Mod
22     #endif
23     
24        use Chem_Mod              ! Chemistry Base Class
25        use Chem_StateMod         ! Chemistry State
26        use Chem_DepositionMod    ! Aerosol Deposition
27        use Chem_ConstMod, only: grav, von_karman, cpd, &
28                                 undefval => undef         ! Constants !
29        use Chem_UtilMod          ! I/O
30        use Chem_MieMod           ! Aerosol LU Tables, calculator
31        use m_inpak90             ! Resource file management
32        use m_die, only: die
33     
34        implicit none
35     
36     ! !PUBLIC TYPES:
37     !
38        PRIVATE
39        PUBLIC  BC_GridComp       ! The BC object 
40     
41     !
42     ! !PUBLIC MEMBER FUNCTIONS:
43     !
44     
45        PUBLIC  BC_GridCompInitialize
46        PUBLIC  BC_GridCompRun
47        PUBLIC  BC_GridCompFinalize
48     
49     !
50     ! !DESCRIPTION:
51     !
52     !  This module implements the (pre-ESMF) BC Grid Component. 
53     !
54     ! !REVISION HISTORY:
55     !
56     !  16Sep2003 da Silva  First crack.
57     !
58     !EOP
59     !-------------------------------------------------------------------------
60     
61       type BC_GridComp
62             character(len=255) :: name
63             type(Chem_Mie), pointer :: mie_tables  ! aod LUTs
64             real, pointer :: biofuel_src(:,:)
65             real, pointer :: biomass_src_(:,:) ! before diurnal
66             real, pointer :: biomass_src(:,:)
67             real, pointer :: ebcant1_src(:,:)  ! level 1
68             real, pointer :: ebcant2_src(:,:)  ! level 2
69             real, pointer :: bc_ship_src(:,:)
70             real :: fHydrophobic         ! Fraction of emissions hydrophobic
71             real :: eBiofuel             ! Emission factor of Biofuel to BC aerosol
72             real :: eBiomassBurning      ! Emission factor of Biomass Burning to BC
73             integer :: nymd   ! date of last emissions/prodction
74             character(len=255) :: bb_srcfilen
75             character(len=255) :: bf_srcfilen
76             character(len=255) :: ebcant1_srcfilen
77             character(len=255) :: ebcant2_srcfilen
78     	character(len=255) :: bc_ship_srcfilen
79             integer :: nymd_bf
80             integer :: nymd_ebcant1
81             integer :: nymd_ebcant2
82     	integer :: nymd_bc_ship
83       end type BC_GridComp
84     
85       real, parameter :: OCEAN=0.0, LAND = 1.0, SEA_ICE = 2.0
86     
87     CONTAINS
88     
89     !-------------------------------------------------------------------------
90     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
91     !-------------------------------------------------------------------------
92     !BOP
93     !
94     ! !IROUTINE:  BC_GridCompInitialize --- Initialize BC_GridComp
95     !
96     ! !INTERFACE:
97     !
98     
99        subroutine BC_GridCompInitialize ( gcBC, w_c, impChem, expChem, &
100                                           nymd, nhms, cdt, rc )
101     
102     ! !USES:
103     
104       implicit NONE
105     
106     ! !INPUT PARAMETERS:
107     
108        type(Chem_Bundle), intent(inout) :: w_c        ! Chemical tracer fields      
109        integer, intent(in) :: nymd, nhms           ! time
110        real, intent(in) :: cdt                     ! chemistry timestep (secs)
111     
112     ! !OUTPUT PARAMETERS:
113     
114        type(BC_GridComp), intent(inout) :: gcBC   ! Grid Component
115        type(ESMF_State), intent(inout)  :: impChem  ! Import State
116        type(ESMF_State), intent(inout)  :: expChem  ! Export State
117        integer, intent(out) ::  rc                  ! Error return code:
118                                                     !  0 - all is well
119                                                     !  1 - 
120     
121     ! !DESCRIPTION: Initializes the BC Grid Component. It primarily sets
122     !               the import state for each active constituent package.
123     !
124     ! !REVISION HISTORY:
125     !
126     !  18Sep2003 da Silva  First crack.
127     !
128     !EOP
129     !-------------------------------------------------------------------------
130     
131        character(len=*), parameter :: myname = 'BC_GridCompInitialize'
132     
133     
134        character(len=255) :: rcfilen = 'BC_GridComp.rc'
135        integer :: ios, n
136        integer :: i1, i2, im, j1, j2, jm, nbins, nbeg, nend, nbins_rc
137        integer :: nTimes, begTime, incSecs
138        integer, allocatable :: ier(:)
139        real, allocatable :: buffer(:,:)
140        real :: qmax, qmin
141     
142     
143        gcBC%name = 'BC Constituent Package'
144     
145     !  Initialize local variables
146     !  --------------------------
147        rc = 0
148        i1 = w_c%grid%i1; i2 = w_c%grid%i2; im = w_c%grid%im
149        j1 = w_c%grid%j1; j2 = w_c%grid%j2; jm = w_c%grid%jm
150        nbins = w_c%reg%n_BC
151        nbeg  = w_c%reg%i_BC
152        nend  = w_c%reg%j_BC
153     
154        call init_()
155        if ( rc /= 0 ) return
156     
157     
158     !                       -------------------
159     !                       Parse resource file
160     !                       -------------------
161     
162     !  Load resource file
163     !  ------------------
164        call i90_loadf ( rcfilen, ier(1) )
165        if ( ier(1) .ne. 0 ) then
166           call final_(10)
167           return
168        end if
169     
170        call i90_label ( 'number_bc_classes:', ier(1) )
171        nbins_rc = i90_gint ( ier(2) )
172        if ( any(ier(1:2) /= 0) ) then
173           call final_(20)
174           return
175        end if
176        if ( nbins_rc /= nbins ) then
177           call final_(25)
178           return
179        end if
180     
181     
182     !  BC sources files
183     !  ---------------------
184        call i90_label ( 'bb_srcfilen:', ier(1) )
185        if ( ier(1) /= 0 ) then
186           call final_(30)
187           return
188        else
189           call i90_gtoken ( gcBC%bb_srcfilen, ier(1) )
190           if ( ier(1) /= 0 ) then
191              call final_(40)
192              return
193           end if
194        end if
195     
196        call i90_label ( 'bf_srcfilen:', ier(1) )
197        if ( ier(1) /= 0 ) then
198           call final_(30)
199           return
200        else
201           call i90_gtoken ( gcBC%bf_srcfilen, ier(1) )
202           if ( ier(1) /= 0 ) then
203              call final_(40)
204              return
205           end if
206        end if
207     
208        call i90_label ( 'ebcant1_srcfilen:', ier(1) )
209        if ( ier(1) /= 0 ) then
210           call final_(30)
211           return
212        else
213           call i90_gtoken ( gcBC%ebcant1_srcfilen, ier(1) )
214           if ( ier(1) /= 0 ) then
215              call final_(40)
216              return
217           end if
218        end if
219     
220        call i90_label ( 'ebcant2_srcfilen:', ier(1) )
221        if ( ier(1) /= 0 ) then
222           call final_(30)
223           return
224        else
225           call i90_gtoken ( gcBC%ebcant2_srcfilen, ier(1) )
226           if ( ier(1) /= 0 ) then
227              call final_(40)
228              return
229           end if
230        end if
231     
232        call i90_label ( 'bc_ship_srcfilen:', ier(1) )
233        if ( ier(1) /= 0 ) then
234           call final_(30)
235           return
236        else
237           call i90_gtoken ( gcBC%bc_ship_srcfilen, ier(1) )
238           if ( ier(1) /= 0 ) then
239              call final_(40)
240              return
241           end if
242        end if
243     
244     !  Hydrophobic fraction
245     !  ---------------
246        call i90_label ( 'hydrophobic_fraction:', ier(1) )
247        gcBC%fHydrophobic = i90_gfloat ( ier(2) )
248        if ( any(ier(1:2) /= 0) ) then
249           call final_(50)
250           return
251        end if
252     
253     
254     !  Biofuel Emission Factor
255     !  ---------------
256        call i90_label ( 'biofuel_emission_factor:', ier(1) )
257        gcBC%eBiofuel = i90_gfloat ( ier(2) )
258        if ( any(ier(1:2) /= 0) ) then
259           call final_(50)
260           return
261        end if
262     
263     
264     !  Biomass Burning Emission Factor
265     !  ---------------
266        call i90_label ( 'biomass_burning_emission_factor:', ier(1) )
267        gcBC%eBiomassBurning = i90_gfloat ( ier(2) )
268        if ( any(ier(1:2) /= 0) ) then
269           call final_(50)
270           return
271        end if
272     
273     
274     !  Scavenging Efficiency
275     !  To be used in convtran.F90, this parameter
276     !  is the scavenging efficiency of the tracer [km -1]
277     !  ---------------
278        call i90_label ( 'fscav:', ier(1) )
279        do n = 1, nbins
280           w_c%reg%fscav(nbeg+n-1) = i90_gfloat ( ier(n+1) )
281        end do
282        if ( any(ier(1:nbins+1) /= 0) ) then
283           call final_(50)
284           return
285        end if
286     !                          -------
287     
288     !  Check initial date of inventory emission/oxidant files
289     !  ------------------------------------------------------
290     !  The intent here is that these files are valid for a particular
291     !  YYYY or YYYYMMDD (if 1x year in file).  We need to request
292     !  the correct date
293        call Chem_UtilGetTimeInfo( gcBC%bf_srcfilen, gcBC%nymd_bf, &
294                                   begTime, nTimes, incSecs )
295        call Chem_UtilGetTimeInfo( gcBC%ebcant1_srcfilen, gcBC%nymd_ebcant1, &
296                                   begTime, nTimes, incSecs )
297        call Chem_UtilGetTimeInfo( gcBC%ebcant2_srcfilen, gcBC%nymd_ebcant2, &
298                                   begTime, nTimes, incSecs )
299        call Chem_UtilGetTimeInfo( gcBC%bc_ship_srcfilen, gcBC%nymd_bc_ship, &
300                                   begTime, nTimes, incSecs )
301        ier(1) = gcBC%nymd_bf
302        ier(2) = gcBC%nymd_ebcant1
303        ier(3) = gcBC%nymd_ebcant2
304        ier(4) = gcBC%nymd_bc_ship
305        if( any(ier(1:4) < 0 ) ) then
306          call final_(60)
307          return
308        endif
309     
310     
311     !  Initialize date for BCs
312     !  -----------------------
313        gcBC%nymd = -1   ! nothing read yet
314     
315     #ifndef GEOS5
316     
317     !\/--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---\/ 
318     
319     !  Set which fvGCM fields are needed
320     !  ---------------------------------
321        call Chem_StateSetNeeded ( impChem, iFRACLAKE, .true., ier(1) )
322        call Chem_StateSetNeeded ( impChem, iGWETTOP,  .true., ier(2) )
323        call Chem_StateSetNeeded ( impChem, iORO,      .true., ier(3) )
324        call Chem_StateSetNeeded ( impChem, iU10M,     .true., ier(4) )
325        call Chem_StateSetNeeded ( impChem, iV10M,     .true., ier(5) )
326        call Chem_StateSetNeeded ( impChem, iLAI,      .true., ier(6) )
327        call Chem_StateSetNeeded ( impChem, iUSTAR,    .true., ier(7) )
328        call Chem_StateSetNeeded ( impChem, iPRECC,    .true., ier(8) )
329        call Chem_StateSetNeeded ( impChem, iPRECL,    .true., ier(9) )
330        call Chem_StateSetNeeded ( impChem, iDQCOND,   .true., ier(10) )
331        call Chem_StateSetNeeded ( impChem, iT,        .true., ier(11) )
332        call Chem_StateSetNeeded ( impChem, iAIRDENS,  .true., ier(12) )
333        call Chem_StateSetNeeded ( impChem, iU,        .true., ier(13) )
334        call Chem_StateSetNeeded ( impChem, iV,        .true., ier(14) )
335        call Chem_StateSetNeeded ( impChem, iPBLH,     .true., ier(15) )
336        call Chem_StateSetNeeded ( impChem, iSHFX,     .true., ier(16) )
337        call Chem_StateSetNeeded ( impChem, iZ0H,      .true., ier(17) )
338        call Chem_StateSetNeeded ( impChem, iHSURF,    .true., ier(18) )
339        call Chem_StateSetNeeded ( impChem, iHGHTE,    .true., ier(19) )
340     
341        if ( any(ier(1:19) /= 0) ) then
342             call final_(60)
343             return
344        endif
345     
346     !  Select fields to be produced in the export state
347     !  ------------------------------------------------
348        n = nbins
349     
350     !  Emission Flux
351        if(n>0) call Chem_StateSetNeeded ( expChem, iBCEM001, .true., ier(1) )
352        if(n>1) call Chem_StateSetNeeded ( expChem, iBCEM002, .true., ier(2) )
353        if(n>2) call Chem_StateSetNeeded ( expChem, iBCEM003, .true., ier(3) )
354        if(n>3) call Chem_StateSetNeeded ( expChem, iBCEM004, .true., ier(4) )
355        if(n>4) call Chem_StateSetNeeded ( expChem, iBCEM005, .true., ier(5) )
356        if(n>5) call Chem_StateSetNeeded ( expChem, iBCEM006, .true., ier(6) )
357        if(n>6) call Chem_StateSetNeeded ( expChem, iBCEM007, .true., ier(7) )
358        if(n>7) call Chem_StateSetNeeded ( expChem, iBCEM008, .true., ier(8) )
359        if(n>8) ier(9) = 1 ! not enough bins - need to change mod_diag.F
360     
361        if ( any(ier(1:9) /= 0) ) then
362             call final_(70)
363             return
364        endif
365     
366     !  Dry Deposition Flux
367        if(n>0) call Chem_StateSetNeeded ( expChem, iBCDP001, .true., ier(1) )
368        if(n>1) call Chem_StateSetNeeded ( expChem, iBCDP002, .true., ier(2) )
369        if(n>2) call Chem_StateSetNeeded ( expChem, iBCDP003, .true., ier(3) )
370        if(n>3) call Chem_StateSetNeeded ( expChem, iBCDP004, .true., ier(4) )
371        if(n>4) call Chem_StateSetNeeded ( expChem, iBCDP005, .true., ier(5) )
372        if(n>5) call Chem_StateSetNeeded ( expChem, iBCDP006, .true., ier(6) )
373        if(n>6) call Chem_StateSetNeeded ( expChem, iBCDP007, .true., ier(7) )
374        if(n>7) call Chem_StateSetNeeded ( expChem, iBCDP008, .true., ier(8) )
375        if(n>8) ier(9) = 1 ! not enough bins - need to change mod_diag.F
376     
377        if ( any(ier(1:9) /= 0) ) then
378             call final_(70)
379             return
380        endif
381     
382     !  Wet Deposition Flux
383        if(n>0) call Chem_StateSetNeeded ( expChem, iBCWT001, .true., ier(1) )
384        if(n>1) call Chem_StateSetNeeded ( expChem, iBCWT002, .true., ier(2) )
385        if(n>2) call Chem_StateSetNeeded ( expChem, iBCWT003, .true., ier(3) )
386        if(n>3) call Chem_StateSetNeeded ( expChem, iBCWT004, .true., ier(4) )
387        if(n>4) call Chem_StateSetNeeded ( expChem, iBCWT005, .true., ier(5) )
388        if(n>5) call Chem_StateSetNeeded ( expChem, iBCWT006, .true., ier(6) )
389        if(n>6) call Chem_StateSetNeeded ( expChem, iBCWT007, .true., ier(7) )
390        if(n>7) call Chem_StateSetNeeded ( expChem, iBCWT008, .true., ier(8) )
391        if(n>8) ier(9) = 1 ! not enough bins - need to change mod_diag.F
392     
393        if ( any(ier(1:9) /= 0) ) then
394             call final_(70)
395             return
396        endif
397     
398     
399     !  Other diagnostics
400        call Chem_StateSetNeeded ( expChem, iBCSMASS, .true., ier(1) )
401        call Chem_StateSetNeeded ( expChem, iBCCMASS, .true., ier(2) )
402        call Chem_StateSetNeeded ( expChem, iBCMASS, .true., ier(3) )
403        call Chem_StateSetNeeded ( expChem, iBCEXTTAU, .true., ier(4) )
404        call Chem_StateSetNeeded ( expChem, iBCSCATAU, .true., ier(5) )
405        call Chem_StateSetNeeded ( expChem, iBCEMAN, .true., ier(6))
406        call Chem_StateSetNeeded ( expChem, iBCEMBB, .true., ier(7))
407        call Chem_StateSetNeeded ( expChem, iBCEMBF, .true., ier(8))
408        call Chem_StateSetNeeded ( expChem, iBCHYPHIL, .true., ier(9))
409     
410        if ( any(ier(1:9) /= 0) ) then
411             call final_(70)
412             return
413        endif
414     
415     !/\--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---/\
416     
417     #endif
418     
419     
420     !  All done
421     !  --------
422        call i90_release()
423        deallocate(ier)
424     
425        return
426     
427     
428     CONTAINS
429     
430        subroutine init_()
431        integer ios, nerr
432        nerr = max ( 32, nbins+1 )
433        allocate ( gcBC%biomass_src(i1:i2,j1:j2), gcBC%biofuel_src(i1:i2,j1:j2), &
434                   gcBC%biomass_src_(i1:i2,j1:j2), &
435                   gcBC%ebcant1_src(i1:i2,j1:j2), gcBC%ebcant2_src(i1:i2,j1:j2), &
436                   gcBC%bc_ship_src(i1:i2,j1:j2), ier(nerr), stat=ios )
437        if ( ios /= 0 ) rc = 100
438        end subroutine init_
439     
440        subroutine final_(ierr)
441        integer :: ierr
442        integer ios
443        deallocate ( gcBC%biomass_src, gcBC%biofuel_src, gcBC%bc_ship_src, &
444                     gcBC%biomass_src_, &
445                     gcBC%ebcant1_src, gcBC%ebcant2_src, ier, stat=ios )
446        call i90_release()
447        rc = ierr
448        end subroutine final_
449     
450        end subroutine BC_GridCompInitialize
451     
452     !-------------------------------------------------------------------------
453     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
454     !-------------------------------------------------------------------------
455     !BOP
456     !
457     ! !IROUTINE:  BC_GridCompRun --- The Chem Driver 
458     !
459     ! !INTERFACE:
460     !
461     
462        subroutine BC_GridCompRun ( gcBC, w_c, impChem, expChem, &
463                                    nymd, nhms, cdt, rc )
464     
465     ! !USES:
466     
467       implicit NONE
468     
469     ! !INPUT/OUTPUT PARAMETERS:
470     
471        type(BC_GridComp), intent(inout) :: gcBC   ! Grid Component
472        type(Chem_Bundle), intent(inout) :: w_c      ! Chemical tracer fields   
473     
474     ! !INPUT PARAMETERS:
475     
476        type(ESMF_State), intent(inout) :: impChem ! Import State
477        integer, intent(in) :: nymd, nhms          ! time
478        real, intent(in) :: cdt                    ! chemistry timestep (secs)
479     
480     ! !OUTPUT PARAMETERS:
481     
482        type(ESMF_State), intent(inout) :: expChem   ! Export State
483        integer, intent(out) ::  rc                  ! Error return code:
484                                                     !  0 - all is well
485                                                     !  1 -
486      
487     ! !DESCRIPTION: This routine implements the so-called BC Driver. That 
488     !               is, adds chemical tendencies to each of the constituents,
489     !  Note: water wapor, the first constituent is not considered a chemical
490     !  constituents.
491     !
492     ! !REVISION HISTORY:
493     !
494     !  18Sep2003 da Silva  First crack.
495     !
496     !EOP
497     !-------------------------------------------------------------------------
498     
499        character(len=*), parameter :: myname = 'BC_GridCompRun'
500        character(len=*), parameter :: Iam = myname
501     
502        integer :: ier(20), idiag
503        integer :: i1, i2, im, j1, j2, jm, nbins, nbeg, nend, km, n, ios
504        integer :: i, j, k, nymd1, nhms1, ijl, ijkl
505        real :: qmax, qmin
506        real :: qUpdate, delq
507        real, pointer :: BC_radius(:), BC_rhop(:)
508     
509     
510     !  Input fields from fvGCM
511     !  -----------------------
512        real, pointer, dimension(:,:) :: fraclake, gwettop, oro, u10m, v10m, &
513                                         xlai, ustar, precc, precl,          &
514                                         pblh, shflux, z0h, hsurf
515        real, pointer, dimension(:,:,:) ::  dqcond, tmpu, rhoa, u, v, hghte
516     
517     #ifdef GEOS5 
518     
519     #define EXPORT     expChem
520     
521     #define ptrBCWT       BC_wet
522     #define ptrBCEM       BC_emis
523     #define ptrBCDP       BC_dep
524     
525     #define ptrBCMASS     BC_mass
526     #define ptrBCEMAN     BC_emisAN
527     #define ptrBCEMBB     BC_emisBB
528     #define ptrBCEMBF     BC_emisBF
529     #define ptrBCHYPHIL   BC_toHydrophilic
530     #define ptrBCSMASS    BC_sfcmass
531     #define ptrBCCMASS    BC_colmass
532     #define ptrBCEXTTAU   BC_exttau
533     #define ptrBCSCATAU   BC_scatau
534     
535        integer :: STATUS
536     
537     #include "BC_GetPointer___.h"
538     
539     #else
540     
541     !\/--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---\/ 
542     
543     !  Quantities to be exported
544     !  -------------------------
545        type(Chem_Array), pointer :: BC_emis(:), BC_dep(:), BC_wet(:), &
546                                     BC_sfcmass, BC_colmass, BC_mass, BC_exttau, &
547                                     BC_scatau, BC_emisAN, BC_emisBB, BC_emisBF, &
548                                     BC_toHydrophilic
549     
550     !/\--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---/\
551     
552     #endif
553     
554     !  Initialize local variables
555     !  --------------------------
556        rc = 0
557        i1 = w_c%grid%i1; i2 = w_c%grid%i2; im = w_c%grid%im
558        j1 = w_c%grid%j1; j2 = w_c%grid%j2; jm = w_c%grid%jm
559     
560        km = w_c%grid%km
561        nbins = w_c%reg%n_BC
562        nbeg  = w_c%reg%i_BC
563        nend  = w_c%reg%j_BC
564     
565        ijl  = ( i2 - i1 + 1 ) * ( j2 - j1 + 1 )
566        ijkl = ijl * km
567     
568     
569     
570     ! Update emissions/production if necessary (daily)
571     !  -----------------------------------------------
572        if(gcBC%nymd .ne. nymd) then
573     
574     !   Biomass Burning -- select on known inventories
575     !   ----------------------------------------------
576     
577     !   Daily files (e.g., MODIS) or GFED v.2 (1997 - 2005 valid)
578         if (  index(gcBC%bb_srcfilen,'%') .gt. 0 .or. &
579               index(gcBC%bb_srcfilen,'gfed') .gt. 0 ) then  
580            nymd1 = nymd
581            nhms1 = 120000
582     
583     !   Assume GFED climatology or Martin (Duncan) climatology
584         else                                            
585            nymd1 = 1971*10000 + mod ( nymd, 10000 )  ! assumes 1971
586     !       nymd1 = nymd
587            nhms1 = 120000
588         end if
589     
590         call Chem_UtilMPread ( gcBC%bb_srcfilen, 'biomass', nymd1, nhms1, &
591                                i1, i2, 0, im, j1, j2, 0, jm, 0, &
592                                var2d=gcBC%biomass_src, cyclic=.true., &
593                                grid = w_c%grid_esmf  )
594     
595     
596     !   Biofuel and anthropogenic emissions (inventories)
597     !   -------------------------------------------------
598         nymd1 = (gcBC%nymd_bf/10000)*10000 + mod ( nymd, 10000 )
599         nhms1 = 120000
600         call Chem_UtilMPread ( gcBC%bf_srcfilen, 'biofuel', nymd1, nhms1, &
601                                i1, i2, 0, im, j1, j2, 0, jm, 0, &
602                                var2d=gcBC%biofuel_src, cyclic=.true., &
603                                grid = w_c%grid_esmf  )
604     
605         nymd1 = gcBC%nymd_ebcant1
606         nhms1 = 120000
607         call Chem_UtilMPread ( gcBC%ebcant1_srcfilen, 'antebc1', nymd1, nhms1, &
608                                i1, i2, 0, im, j1, j2, 0, jm, 0, &
609                                var2d=gcBC%ebcant1_src, cyclic=.true., &
610                                grid = w_c%grid_esmf  )
611     
612     !   Functionality for only a single layer of emissions
613         if( index(gcBC%ebcant2_srcfilen,'--') .gt. 0) then
614          gcBC%ebcant2_src(i1:i2,j1:j2) = 0.
615         else
616          nymd1 = gcBC%nymd_ebcant2
617          nhms1 = 120000
618          call Chem_UtilMPread ( gcBC%ebcant2_srcfilen, 'antebc2', nymd1, nhms1, &
619                                 i1, i2, 0, im, j1, j2, 0, jm, 0, &
620                                 var2d=gcBC%ebcant2_src, grid = w_c%grid_esmf  )
621     
622         endif
623     
624     !   Ship based BC emissions
625         nymd1 = gcBC%nymd_bc_ship
626         nhms1 = 120000
627         call Chem_UtilMPread ( gcBC%bc_ship_srcfilen, 'bc_ship', nymd1, nhms1, &
628                                i1, i2, 0, im, j1, j2, 0, jm, 0, &
629                                var2d=gcBC%bc_ship_src, cyclic=.true., &
630                                grid = w_c%grid_esmf  )
631     
632     
633     
634     !   As a safety check, where value is undefined set to 0
635         do j = j1, j2
636          do i = i1, i2
637           if(1.01*gcBC%biomass_src(i,j) .gt. undefval) gcBC%biomass_src(i,j) = 0.
638           if(1.01*gcBC%biofuel_src(i,j) .gt. undefval) gcBC%biofuel_src(i,j) = 0.
639           if(1.01*gcBC%ebcant1_src(i,j) .gt. undefval) gcBC%ebcant1_src(i,j) = 0.
640           if(1.01*gcBC%ebcant2_src(i,j) .gt. undefval) gcBC%ebcant2_src(i,j) = 0.
641           if(1.01*gcBC%bc_ship_src(i,j) .gt. undefval) gcBC%bc_ship_src(i,j) = 0.
642          enddo
643         enddo
644     
645     
646     #ifdef DEBUG
647         call pmaxmin ( 'BC: biomass', gcBC%biomass_src, qmin, qmax, ijl,1, 1. )
648         call pmaxmin ( 'BC: biofuel', gcBC%biofuel_src, qmin, qmax, ijl,1, 1. )
649         call pmaxmin ( 'BC: ebcant1', gcBC%ebcant1_src, qmin, qmax, ijl,1,1.)
650         call pmaxmin ( 'BC: ebcant2', gcBC%ebcant2_src, qmin, qmax, ijl,1,1.)
651         call pmaxmin ( 'BC: bc_ship', gcBC%bc_ship_src, qmin, qmax, ijl,1,1.)
652     #endif
653     
654     !   Save this in case we need to apply diurnal cycle
655     !   ------------------------------------------------
656        if ( w_c%diurnal_bb ) then
657             gcBC%biomass_src_(:,:) = gcBC%biomass_src(:,:)
658        end if
659     
660         gcBC%nymd = nymd
661     
662        endif
663     
664     
665     !  Apply diurnal cycle if so desired
666     !  ---------------------------------
667        if ( w_c%diurnal_bb ) then
668           call Chem_BiomassDiurnal ( gcBC%biomass_src, gcBC%biomass_src_,   &
669                                      w_c%grid%lon(:), w_c%grid%lat(:), nhms, cdt )      
670        end if
671     
672     #ifndef GEOS5
673     
674     !\/--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---\/ 
675     
676     !  Work space for holding bc output
677     !  ----------------------------------
678        allocate ( BC_emis(nbins),  BC_dep(nbins), BC_wet(nbins), &
679                   BC_sfcmass, BC_colmass, BC_mass, BC_exttau, BC_scatau, &
680                   BC_emisAN, BC_emisBB, BC_emisBF, BC_toHydrophilic, &
681                   stat = ios )
682        if ( ios /= 0 ) then
683           rc = 1
684           return
685        end if
686     
687     !/\--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---/\
688     
689     #endif
690     
691     #ifdef DEBUG
692        do n = nbeg, nend
693           call pmaxmin ( 'BC: q_beg', w_c%qa(n)%data3d(i1:i2,j1:j2,1:km), qmin, qmax, &
694                        ijl, km, 1. )
695        end do
696     #endif
697     
698     !  BC particle radius [m] and density [kg m-3]
699     !  ---------------------------------------------
700     !  For now these are dummy values and unused below
701        allocate( BC_radius(nbins), BC_rhop(nbins) )
702        BC_radius(:) = 0.5e-6
703        BC_rhop(:)   = 1000.
704     
705     
706     #ifdef GEOS5
707     
708     !  Get 2D Imports
709     !  --------------
710        call MAPL_GetPointer ( impChem, fraclake, 'FRLAKE',  rc=ier(1) )
711        call MAPL_GetPointer ( impChem, gwettop,  'WET1',    rc=ier(2) )
712        call MAPL_GetPointer ( impChem, oro,      'LWI',     rc=ier(3) )
713        call MAPL_GetPointer ( impChem, u10m,     'U10M',    rc=ier(4) )
714        call MAPL_GetPointer ( impChem, v10m,     'V10M',    rc=ier(5) )
715        call MAPL_GetPointer ( impChem, xlai,     'LAI',     rc=ier(6) )
716        call MAPL_GetPointer ( impChem, ustar,    'USTAR',   rc=ier(7) )
717        call MAPL_GetPointer ( impChem, precc,    'CN_PRCP', rc=ier(8) )
718        call MAPL_GetPointer ( impChem, precl,    'NCN_PRCP',   rc=ier(9) )
719        call MAPL_GetPointer ( impChem, pblh,     'ZPBL',    rc=ier(10) )
720        call MAPL_GetPointer ( impChem, shflux,   'SH',      rc=ier(11) )
721        call MAPL_GetPointer ( impChem, z0h,      'Z0H',     rc=ier(12) )
722        ier(13) = 0 ! see below for hsurf
723     
724     !  Get 3D Imports
725     !  --------------
726        call MAPL_GetPointer ( impChem, dqcond, 'DQDT',    rc=ier(14) )
727        call MAPL_GetPointer ( impChem, tmpu,   'T',       rc=ier(15) )
728        call MAPL_GetPointer ( impChem, rhoa,   'AIRDENS', rc=ier(16) )
729        call MAPL_GetPointer ( impChem, u,      'U',       rc=ier(17) )
730        call MAPL_GetPointer ( impChem, v,      'V',       rc=ier(18) )
731        call MAPL_GetPointer ( impChem, hghte,  'ZLE',     rc=ier(19) )
732     
733     !  Unlike GEOS-4 hghte is defined for km+1
734     !  ---------------------------------------
735        hsurf => hghte(i1:i2,j1:j2,km) ! Recall: GEOS-5 has edges with k in [0,km]
736         
737     
738        if ( .not. associated(gwettop) ) &
739             call write_parallel('gwettop NOT associated')
740     
741     #else
742     
743     !\/--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---\/ 
744     
745     !  Get input fvGCM 2D diagnostics
746     !  ------------------------------
747        call Chem_StateGetArray2D ( impChem, iFRACLAKE, fraclake, ier(1) )
748        call Chem_StateGetArray2D ( impChem, iGWETTOP,  gwettop,  ier(2) )
749        call Chem_StateGetArray2D ( impChem, iORO,      oro,      ier(3) )
750        call Chem_StateGetArray2D ( impChem, iU10M,     u10m,     ier(4) )
751        call Chem_StateGetArray2D ( impChem, iV10M,     v10m,     ier(5) )
752        call Chem_StateGetArray2D ( impChem, iLAI,      xlai,     ier(6) )
753        call Chem_StateGetArray2D ( impChem, iUSTAR,    ustar,    ier(7) )
754        call Chem_StateGetArray2D ( impChem, iPRECC,    precc,    ier(8) )
755        call Chem_StateGetArray2D ( impChem, iPRECL,    precl,    ier(9) )
756        call Chem_StateGetArray2D ( impChem, iPBLH,     pblh,     ier(10) )
757        call Chem_StateGetArray2D ( impChem, iSHFX,     shflux,   ier(11) )
758        call Chem_StateGetArray2D ( impChem, iZ0H,      z0h,      ier(12) )
759        call Chem_StateGetArray2D ( impChem, iHSURF,    hsurf,      ier(13) )
760     
761     !  Get input fvGCM 3D diagnostics
762     !  ------------------------------
763        call Chem_StateGetArray3D ( impChem, iDQCOND,   dqcond,   ier(14) )
764        call Chem_StateGetArray3D ( impChem, iT,        tmpu,     ier(15) )
765        call Chem_StateGetArray3D ( impChem, iAIRDENS,  rhoa,     ier(16) )
766        call Chem_StateGetArray3D ( impChem, iU,        u,        ier(17) )
767        call Chem_StateGetArray3D ( impChem, iV,        v,        ier(18) )
768        call Chem_StateGetArray3D ( impChem, iHGHTE,    hghte,    ier(19) )
769     
770     !/\--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---/\
771     
772     #endif
773     
774        if ( any(ier(1:19) /= 0) ) then
775             rc = 10 
776             return
777        end if
778     
779        if ( any(ier(1:19) /= 0) ) then
780             rc = 10 
781             return
782        end if
783     
784     !  Make sure LAI has values over ocean
785     !  -----------------------------------
786        where ( oro /= LAND  )  xlai = 0.0
787     
788     #ifdef DEBUG
789     
790        call pmaxmin('BC: fraclake   ', fraclake, qmin, qmax, ijl,1, 1. )
791        call pmaxmin('BC: gwtop      ', gwettop , qmin, qmax, ijl,1, 1. )
792        call pmaxmin('BC: oro        ', oro     , qmin, qmax, ijl,1, 1. )
793        call pmaxmin('BC: u10m       ', u10m    , qmin, qmax, ijl,1, 1. )
794        call pmaxmin('BC: v10m       ', v10m    , qmin, qmax, ijl,1, 1. )
795        call pmaxmin('BC: xlai       ', xlai    , qmin, qmax, ijl,1, 1. )
796        call pmaxmin('BC: ustar      ', ustar   , qmin, qmax, ijl,1, 1. )
797        call pmaxmin('BC: precc      ', precc   , qmin, qmax, ijl,1, 1. )
798        call pmaxmin('BC: precl      ', precl   , qmin, qmax, ijl,1, 1. )
799        call pmaxmin('BC: pblh       ', pblh    , qmin, qmax, ijl,1, 1. )
800        call pmaxmin('BC: shfflux    ', shflux  , qmin, qmax, ijl,1, 1. )
801        call pmaxmin('BC: z0h        ', z0h     , qmin, qmax, ijl,1, 1. )
802        call pmaxmin('BC: hsurf      ', hsurf   , qmin, qmax, ijl,1, 1. )
803     
804        call pmaxmin('BC: dqcond     ', dqcond  , qmin, qmax, ijkl,1, 1. )
805        call pmaxmin('BC: tmpu       ', tmpu    , qmin, qmax, ijkl,1, 1. )
806        call pmaxmin('BC: rhoa       ', rhoa    , qmin, qmax, ijkl,1, 1. )
807        call pmaxmin('BC: u          ', u       , qmin, qmax, ijkl,1, 1. )
808        call pmaxmin('BC: v          ', v       , qmin, qmax, ijkl,1, 1. )
809        call pmaxmin('BC: hghte      ', hghte   , qmin, qmax, ijkl,1, 1. )
810     
811     #endif
812     
813     #ifndef GEOS5
814     
815     !\/--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---\/ 
816     
817     
818     !  Get pointers to export state
819     !  ----------------------------
820        do n = 1, nbins
821           idiag = iBCEM001 + n - 1
822           call Chem_StateGetArray2D ( expChem, idiag, BC_emis(n)%data2d, ier(n) )
823        end do
824        if ( any(ier(1:nbins) /= 0) ) then
825             rc = 15 
826             return
827        end if
828     
829        do n = 1, nbins
830           idiag = iBCDP001 + n - 1
831           call Chem_StateGetArray2D ( expChem, idiag, BC_dep(n)%data2d, ier(n) )
832        end do
833        if ( any(ier(1:nbins) /= 0) ) then
834             rc = 15 
835             return
836        end if
837     
838        do n = 1, nbins
839           idiag = iBCWT001 + n - 1
840           call Chem_StateGetArray2D ( expChem, idiag, BC_wet(n)%data2d, ier(n) )
841        end do
842        if ( any(ier(1:nbins) /= 0) ) then
843             rc = 15 
844             return
845        end if
846     
847        idiag = iBCSMASS
848        call Chem_StateGetArray2D ( expChem, idiag, BC_sfcmass%data2d, ier(1) )
849        if ( ier(1) /= 0 ) then
850             rc = 15 
851             return
852        end if
853     
854        idiag = iBCCMASS
855        call Chem_StateGetArray2D ( expChem, idiag, BC_colmass%data2d, ier(1) )
856        if ( ier(1) /= 0 ) then
857             rc = 15 
858             return
859        end if
860     
861        idiag = iBCMASS
862        call Chem_StateGetArray3D ( expChem, idiag, BC_mass%data3d, ier(1) )
863        if ( ier(1) /= 0 ) then
864             rc = 15 
865             return
866        end if
867     
868        idiag = iBCEXTTAU
869        call Chem_StateGetArray2D ( expChem, idiag, BC_exttau%data2d, ier(1) )
870        if ( ier(1) /= 0 ) then
871             rc = 15 
872             return
873        end if
874     
875        idiag = iBCSCATAU
876        call Chem_StateGetArray2D ( expChem, idiag, BC_scatau%data2d, ier(1) )
877        if ( ier(1) /= 0 ) then
878             rc = 15 
879             return
880        end if
881     
882        idiag = iBCEMAN
883        call Chem_StateGetArray2D ( expChem, idiag, BC_emisAN%data2d, ier(1) )
884        if ( ier(1) /= 0 ) then
885             rc = 15 
886             return
887        end if
888     
889        idiag = iBCEMBF
890        call Chem_StateGetArray2D ( expChem, idiag, BC_emisBF%data2d, ier(1) )
891        if ( ier(1) /= 0 ) then
892             rc = 15 
893             return
894        end if
895     
896        idiag = iBCEMBB
897        call Chem_StateGetArray2D ( expChem, idiag, BC_emisBB%data2d, ier(1) )
898        if ( ier(1) /= 0 ) then
899             rc = 15 
900             return
901        end if
902     
903        idiag = iBCHYPHIL
904        call Chem_StateGetArray2D ( expChem, idiag, BC_toHydrophilic%data2d, ier(1) )
905        if ( ier(1) /= 0 ) then
906             rc = 15 
907             return
908        end if
909     
910     !/\--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---/\
911     
912     #endif
913     
914     
915     !  BC Source
916     !  -----------
917        call BC_Emission ( i1, i2, j1, j2, km, nbins, cdt, gcBC, w_c, &
918                           pblh, tmpu, rhoa, BC_emis, &
919                           BC_emisAN, BC_emisBB, BC_emisBF, rc )
920     
921     #ifdef DEBUG
922        do n = nbeg, nend
923           call pmaxmin('BC: q_emi', w_c%qa(n)%data3d(i1:i2,j1:j2,1:km), &
924                         qmin, qmax, ijl, km, 1. )
925        end do
926     #endif
927     
928     !  Ad Hoc transfer of hydrophobic to hydrophilic aerosols
929     !  Following Chin's parameterization, the rate constant is
930     !  k = 4.63e-6 s-1 (.4 day-1; e-folding time = 2.5 days)
931        if(associated(BC_toHydrophilic%data2d)) &
932          BC_toHydrophilic%data2d(i1:i2,j1:j2) = 0.0
933     
934        do k = 1, km
935         do j = j1, j2
936          do i = i1, i2
937           qUpdate = w_c%qa(nbeg)%data3d(i,j,k)*exp(-4.63e-6*cdt)
938           qUpdate = max(qUpdate,1e-32)
939           delq = max(0.,w_c%qa(nbeg)%data3d(i,j,k)-qUpdate)
940           w_c%qa(nbeg)%data3d(i,j,k) = qUpdate
941           w_c%qa(nend)%data3d(i,j,k) = w_c%qa(nend)%data3d(i,j,k)+delq
942           if(associated(BC_toHydrophilic%data2d)) &
943            BC_toHydrophilic%data2d(i,j) = BC_toHydrophilic%data2d(i,j) &
944             + delq*w_c%delp(i,j,k)/grav/cdt
945          end do
946         end do
947        end do
948     
949     !  BC Deposition
950     !  -----------
951        call Chem_Deposition( i1, i2, j1, j2, km, nbeg, nend, nbins, cdt, w_c, &
952                              BC_radius, BC_rhop, tmpu, rhoa, hsurf, hghte, oro, ustar, &
953                              u10m, v10m, fraclake, gwettop, pblh, shflux, z0h, BC_dep, rc )
954     
955     !  BC Wet Removal
956     !  -----------
957        call BC_Wet_Removal ( i1, i2, j1, j2, km, nbins, cdt, rhoa, gcBC, w_c, &
958                              precc, precl, dqcond, tmpu, BC_wet, rc )
959     
960     !  Compute the desired output diagnostics here
961     !  Ideally this will go where chemout is called in fvgcm.F since that
962     !  will reflect the distributions after transport, etc.
963     !  ------------------------------------------------------------------
964        call BC_Compute_Diags(i1, i2, j1, j2, km, nbins, gcBC, w_c, tmpu, rhoa, &
965                              BC_sfcmass, BC_colmass, BC_mass, BC_exttau, &
966                              BC_scatau, rc)
967     
968     !  Clean up
969     !  --------
970        deallocate(BC_radius, BC_rhop, stat=ios)
971     
972     #ifndef GEOS5
973     
974     !\/--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---\/ 
975     
976        deallocate(BC_emis, BC_dep, BC_wet, BC_sfcmass, BC_colmass, BC_mass, &
977                   BC_emisAN, BC_emisBB, BC_emisBF, BC_toHydrophilic, &
978                   BC_exttau, BC_scatau, stat=ios)
979     
980     !/\--- cut --- --- cut --- --- cut --- --- cut --- --- cut --- --- cut ---/\
981     
982     #endif
983     
984        return
985     
986     CONTAINS
987     
988     !-------------------------------------------------------------------------
989     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
990     !-------------------------------------------------------------------------
991     !BOP
992     !
993     ! !IROUTINE:  BC_Emission - Adds Black Carbon emission for one timestep
994     !             We have emissions from 4 sources, which are distributed
995     !             differently in the vertical
996     !             1) biomass burning - uniformly mixed in PBL
997     !             2) biofuel sources - emitted into lowest 100 m
998     !             3) anthropogenic l1 - emitted into lowest 100 m
999     !             4) anthropogenic l2 - emitted into 100 - 500 m levels
1000     !
1001     ! !INTERFACE:
1002     !
1003     
1004        subroutine BC_Emission ( i1, i2, j1, j2, km, nbins, cdt, gcBC, w_c, &
1005                                 pblh, tmpu, rhoa, BC_emis, &
1006                                 BC_emisAN, BC_emisBB, BC_emisBF, rc )
1007     
1008     ! !USES:
1009     
1010       implicit NONE
1011     
1012     ! !INPUT PARAMETERS:
1013     
1014        integer, intent(in) :: i1, i2, j1, j2, km, nbins
1015        real, intent(in)    :: cdt
1016        type(BC_GridComp), intent(in)    :: gcBC       ! BC Grid Component
1017        real, pointer, dimension(:,:)    :: pblh
1018        real, pointer, dimension(:,:,:)  :: tmpu
1019        real, pointer, dimension(:,:,:)  :: rhoa
1020     
1021     ! !OUTPUT PARAMETERS:
1022     
1023        type(Chem_Bundle), intent(inout) :: w_c         ! Chemical tracer fields
1024        type(Chem_Array), intent(inout)  :: BC_emis(nbins) ! BC emissions, kg/m2/s
1025        type(Chem_Array), intent(inout)  :: BC_emisAN      ! BC emissions, kg/m2/s
1026        type(Chem_Array), intent(inout)  :: BC_emisBB      ! BC emissions, kg/m2/s
1027        type(Chem_Array), intent(inout)  :: BC_emisBF      ! BC emissions, kg/m2/s
1028        integer, intent(out)             :: rc          ! Error return code:
1029                                                        !  0 - all is well
1030                                                        !  1 - 
1031        character(len=*), parameter :: myname = 'BC_Emission'
1032     
1033     ! !DESCRIPTION: Updates the BC concentration with emissions every timestep
1034     !
1035     ! !REVISION HISTORY:
1036     !
1037     !  06Nov2003, Colarco
1038     !  Based on Ginoux
1039     !
1040     !EOP
1041     !-------------------------------------------------------------------------
1042     
1043     ! !Local Variables
1044        integer  ::  i, j, k, m, n, ios, ijl
1045        integer  ::  n1, n2
1046                                            ! pressure at 100m, 500m, & PBLH
1047        real, dimension(i1:i2,j1:j2) :: p100, p500, pPblh  
1048        real, dimension(i1:i2,j1:j2) :: p0, z0, ps
1049        real :: p1, z1, dz, delz, delp, f100, f500, fPblh
1050        real :: qmax, qmin, eBiofuel, eBiomass
1051     
1052        real, dimension(i1:i2,j1:j2) :: factor, srcHydrophobic, srcHydrophilic
1053        real, dimension(i1:i2,j1:j2) :: srcBiofuel, srcBiomass, srcAnthro
1054        real                         :: srcAll, zpbl, maxAll
1055     
1056     !  Initialize local variables
1057     !  --------------------------
1058        n1  = w_c%reg%i_BC
1059        n2  = w_c%reg%j_BC
1060        ijl = ( i2 - i1 + 1 ) * ( j2 - j1 + 1 )
1061        eBiomass = gcBC%eBiomassBurning
1062        eBiofuel = gcBC%eBiofuel
1063     
1064     !  Zero diagnostic accumulators
1065        do n = 1, nbins
1066          if( associated(BC_emis(n)%data2d) ) BC_emis(n)%data2d = 0.0
1067        end do
1068          if(associated(BC_emisAN%data2d) )   BC_emisAN%data2d  = 0.0
1069          if(associated(BC_emisBF%data2d) )   BC_emisBF%data2d  = 0.0
1070          if(associated(BC_emisBB%data2d) )   BC_emisBB%data2d  = 0.0
1071     
1072     !  Determine surface pressure
1073     !  AMS Note: pass this in
1074     !  --------------------------
1075        ps = 0.0
1076        do k = 1, km
1077         ps(i1:i2,j1:j2) = ps(i1:i2,j1:j2) + w_c%delp(i1:i2,j1:j2,k)
1078        end do
1079     
1080     !  Find the pressure of the 100m, 500m, and PBLH altitudes
1081     !  AMS Note: this could be greatly simplified by using ze/zm and having a
1082     !      generic routine from the bottom up with an early exit condition
1083     !  -----------------------------------------------------------------------
1084        p0 = ps  
1085        z0(i1:i2,j1:j2) = 0.
1086        do k = km, 1, -1
1087         do j = j1, j2
1088          do i = i1, i2
1089           p1 = p0(i,j) - w_c%delp(i,j,k)
1090           dz = w_c%delp(i,j,k)/rhoa(i,j,k)/grav
1091           z1 = z0(i,j)+dz
1092           if(z0(i,j) .lt. 100 .and. z1 .ge. 100.) then
1093            delz = z1-100.
1094            delp = delz*rhoa(i,j,k)*grav
1095            p100(i,j) = p1+delp
1096           endif
1097           if(z0(i,j) .lt. 500 .and. z1 .ge. 500.) then
1098            delz = z1-500.
1099            delp = delz*rhoa(i,j,k)*grav
1100            p500(i,j) = p1+delp
1101           endif
1102           zpbl = max ( pblh(i,j), 100. )
1103           if(z0(i,j) .lt. zpbl .and. z1 .ge. zpbl ) then
1104            delz = z1-zpbl
1105            delp = delz*rhoa(i,j,k)*grav
1106            pPblh(i,j) = p1+delp
1107           endif
1108           p0(i,j) = p1
1109           z0(i,j) = z1
1110          end do
1111         end do
1112        end do
1113     
1114     #if 0
1115        call pmaxmin ( 'BC: p100   ', p100,  qmin, qmax, ijl, 1, 1. )
1116        call pmaxmin ( 'BC: p500   ', p500,  qmin, qmax, ijl, 1, 1. )
1117        call pmaxmin ( 'BC: pPBL   ', pPBLh, qmin, qmax, ijl, 1, 1. )
1118     #endif
1119     
1120     !  Now update the tracer mixing ratios with the aerosol sources
1121     !  ------------------------------------------------------------
1122        p0 = ps
1123     K_LOOP: do k = km, 1, -1
1124     
1125     !!!    print *, 'BC_Emissions: getting emissions for layer ', k
1126     
1127     !   First determine emissions for this layer
1128     !   ----------------------------------------
1129         do j = j1, j2
1130          do i = i1, i2
1131     
1132           p1 = p0(i,j) - w_c%delp(i,j,k)
1133     
1134     !     Pressure @ 100m
1135     !     ---------------
1136           f100 = 0.
1137           if(p1 .ge. p100(i,j)) f100 = w_c%delp(i,j,k)/(ps(i,j)-p100(i,j))
1138           if(p1 .lt. p100(i,j) .and. p0(i,j) .ge. p100(i,j)) &
1139            f100 = (p0(i,j)-p100(i,j))/(ps(i,j)-p100(i,j))
1140     
1141     !     Pressure @ 500m
1142     !     ---------------
1143           f500 = 0.
1144           if ( p0(i,j) .ge. p100(i,j) .and. p1 .lt. p100(i,j) .and. p1 .ge. p500(i,j)) &
1145            f500 = (p100(i,j)-p1)/(p100(i,j)-p500(i,j))
1146           if(p0(i,j) .lt. p100(i,j) .and. p1 .ge. p500(i,j)) &
1147            f500 = w_c%delp(i,j,k)/(p100(i,j)-p500(i,j))
1148           if(p0(i,j) .ge. p500(i,j) .and. p1 .lt. p500(i,j)) &
1149            f500 = (p0(i,j)-p500(i,j))/(p100(i,j)-p500(i,j))
1150     
1151     !     Pressure @ PBL height
1152     !     ---------------------
1153           fPblh = 0.
1154           if(p1 .ge. pPblh(i,j)) fPblh = w_c%delp(i,j,k)/(ps(i,j)-pPblh(i,j))
1155           if(p1 .lt. pPblh(i,j) .and. p0(i,j) .ge. pPblh(i,j)) &
1156            fPblh = (p0(i,j)-pPblh(i,j))/(ps(i,j)-pPblh(i,j))
1157     
1158     !     Sources by class in kg m-2 s-1
1159     !     ------------------------------
1160           srcBiofuel(i,j) = f100 *eBiofuel*gcBC%biofuel_src(i,j)
1161           srcAnthro(i,j)  = f100 *         gcBC%ebcant1_src(i,j) &
1162                           + f500 *         gcBC%ebcant2_src(i,j) &
1163     		      + f100 *         gcBC%bc_ship_src(i,j)
1164           srcBiomass(i,j) = fPblh*eBiomass*gcBC%biomass_src(i,j)
1165     
1166           srcAll = srcBiofuel(i,j) + srcAnthro(i,j) + srcBiomass(i,j)
1167           srcHydrophobic(i,j) =     gcBC%fHydrophobic  * srcAll
1168           srcHydrophilic(i,j) = (1.-gcBC%fHydrophobic) * srcAll
1169     
1170     !     Update pressure of lower level
1171     !     ------------------------------
1172           p0(i,j) = p1
1173     
1174     #ifndef GEOS5
1175           maxAll = max( maxAll, max( srcHydrophobic(i,j), srcHydrophilic(i,j)))
1176     #endif
1177     
1178          end do ! i
1179         end do  ! j
1180     
1181     #ifdef GEOS5
1182     !   Determine global max/min
1183     !   ------------------------
1184         call pmaxmin ( 'BC: Phobic ', srcHydrophobic, qmin, qmax, ijl, 1, 0. )
1185         maxAll = abs(qmax) + abs(qmin)
1186         call pmaxmin ( 'BC: Philic ', srcHydrophilic, qmin, qmax, ijl, 1, 0. )
1187         maxAll = max ( maxAll, abs(qmax) + abs(qmin) )
1188     #endif
1189     
1190     !   If emissions are zero at this level (globally), we are done
1191     !   -----------------------------------------------------------
1192         if ( maxAll .eq. 0.0 ) exit K_LOOP
1193     
1194     !   Update concentrations at this layer
1195     !   The "1" element is hydrophobic 
1196     !   The "2" element is hydrophilic
1197     !   -----------------------------------    
1198         factor = cdt * grav / w_c%delp(:,:,k)
1199     
1200         w_c%qa(n1)%data3d(:,:,k) = w_c%qa(n1)%data3d(:,:,k) & 
1201                                  + factor * srcHydrophobic 
1202     
1203         w_c%qa(n2)%data3d(:,:,k) = w_c%qa(n2)%data3d(:,:,k) & 
1204                                  + factor * srcHydrophilic
1205     
1206     !   Fill in diagnostics if requested
1207     !   --------------------------------
1208         if ( associated(BC_emis(1)%data2d)) &
1209                         BC_emis(1)%data2d = BC_emis(1)%data2d + srcHydrophobic
1210     
1211         if ( associated(BC_emis(2)%data2d)) &
1212                         BC_emis(2)%data2d = BC_emis(2)%data2d + srcHydrophilic
1213     
1214         if ( associated(BC_emisBF%data2d)) &
1215                         BC_emisBF%data2d  = BC_emisBF%data2d  + srcBiofuel
1216     
1217         if ( associated(BC_emisBB%data2d)) &
1218                         BC_emisBB%data2d  = BC_emisBB%data2d  + srcBiomass
1219     
1220         if ( associated(BC_emisAN%data2d)) &
1221                         BC_emisAN%data2d  = BC_emisAN%data2d  + srcAnthro
1222     
1223        end do K_LOOP
1224     
1225        rc = 0
1226     
1227        end subroutine BC_Emission
1228     
1229     !-------------------------------------------------------------------------
1230     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
1231     !-------------------------------------------------------------------------
1232     !BOP
1233     !
1234     ! !IROUTINE:  BC_Wet_Removal - Removal of dust by precipitation
1235     !  NOTE: For the removal term, fluxout is the sum of the in-cloud
1236     !        convective and large-scale washout and the total flux across
1237     !        the surface due to below-cloud (rainout) convective and
1238     !        large-scale precipitation reaching the surface.  The fluxout
1239     !        is initialized to zero at the beginning and then at each i, j
1240     !        grid point it is added to.
1241     !        
1242     !
1243     ! !INTERFACE:
1244     !
1245     
1246        subroutine BC_Wet_Removal ( i1, i2, j1, j2, km, nbins, cdt, rhoa,gcBC, w_c, &
1247                                    precc, precl, dqcond, tmpu, fluxout, rc )
1248     
1249     ! !USES:
1250     
1251       implicit NONE
1252     
1253     ! !INPUT PARAMETERS:
1254     
1255        integer, intent(in) :: i1, i2, j1, j2, km, nbins
1256        real, intent(in)    :: cdt
1257        type(BC_GridComp), intent(in)   :: gcBC  ! BC Grid Component
1258        real, pointer, dimension(:,:)   :: precc ! total convective precip [mm day-1]
1259        real, pointer, dimension(:,:)   :: precl ! total large-scale prec. [mm day-1]
1260        real, pointer, dimension(:,:,:) :: dqcond  ! change in q due to moist
1261                                                   ! processes [kg kg-1 s-1] 
1262        real, pointer, dimension(:,:,:) :: tmpu    ! temperature [K]
1263        real, pointer, dimension(:,:,:) :: rhoa    ! air density [kg m-3]
1264     
1265     ! !OUTPUT PARAMETERS:
1266     
1267        type(Chem_Bundle), intent(inout) :: w_c        ! Chemical tracer fields
1268        type(Chem_Array), intent(inout)  :: fluxout(nbins) ! Mass lost by wet dep
1269                                                           ! to surface, kg/m2/s
1270        integer, intent(out)             :: rc         ! Error return code:
1271                                                       !  0 - all is well
1272                                                       !  1 - 
1273        character(len=*), parameter :: myname = 'BC_Wet_Removal'
1274     
1275     ! !DESCRIPTION: Updates the dust concentration in each vertical layer
1276     !               due to wet removal
1277     !
1278     ! !REVISION HISTORY:
1279     !
1280     !  17Nov2003, Colarco
1281     !  Based on Ginoux
1282     !
1283     !EOP
1284     !-------------------------------------------------------------------------
1285     
1286     ! !Local Variables
1287        integer  ::  i, j, k, iit, n, LH, kk, ios
1288        integer  ::  nbeg, nend
1289        integer  ::  nHydrophilic
1290        real :: pdog(i1:i2,j1:j2,km)      ! air mass factor dp/g [kg m-2]
1291        real :: Td_ls, Td_cv              ! ls and cv timescales [s]
1292        real :: pls, pcv, pac             ! ls, cv, tot precip [mm day-1]
1293        real :: qls(km), qcv(km)          ! ls, cv portion dqcond [kg m-3 s-1]
1294        real :: qmx, qd, A                ! temporary variables on moisture
1295        real :: F, B, BT                  ! temporary variables on cloud, freq.
1296        real, allocatable :: fd(:,:)      ! flux across layers [kg m-2]
1297        real, allocatable :: DC(:)        ! scavenge change in mass mixing ratio
1298     
1299     !  Rain parameters (from where?)
1300        real, parameter :: B0_ls = 1.0e-4
1301        real, parameter :: F0_ls = 1.0
1302        real, parameter :: XL_ls = 5.0e-4
1303        real, parameter :: B0_cv = 1.5e-3
1304        real, parameter :: F0_cv = 0.3
1305        real, parameter :: XL_cv = 2.0e-3
1306     
1307        rc=0
1308     
1309     !  Initialize local variables
1310     !  --------------------------
1311        do n = 1, nbins
1312         if( associated(fluxout(n)%data2d) ) fluxout(n)%data2d(i1:i2,j1:j2) = 0.0
1313        end do
1314     
1315        nbeg  = w_c%reg%i_BC
1316        nend  = w_c%reg%j_BC
1317        nHydrophilic = 2
1318     
1319     !  Allocate the dynamic arrays
1320        allocate(fd(km,nbins),stat=ios)
1321        if(ios .ne. 0) stop
1322        allocate(dc(nbins),stat=ios)
1323        if(ios .ne. 0) stop
1324     
1325     !  Duration of rain: ls = model timestep, cv = 1800 s (<= cdt)
1326        Td_ls = cdt
1327        Td_cv = 1800.
1328     
1329     !  Accumulate the 3-dimensional arrays of rhoa and pdog
1330        pdog = w_c%delp/grav
1331     
1332     !  Loop over spatial indices
1333        do j = j1, j2
1334         do i = i1, i2
1335     
1336     !    Check for total precipitation amount
1337     !    Assume no precip in column if precl+precc = 0
1338          pac = precl(i,j) + precc(i,j)
1339          if(pac .le. 0.) goto 100
1340          pls = precl(i,j)
1341          pcv = precc(i,j)
1342     
1343     !    Initialize the precipitation fields
1344          qls(:)  = 0.
1345          qcv(:)  = 0.
1346          fd(:,:) = 0.
1347     
1348     !    Find the highest model layer experiencing rainout.  Assumes no
1349     !    scavenging if T < 258 K
1350          LH = 0
1351          do k = 1, km
1352           if(dqcond(i,j,k) .lt. 0. .and. tmpu(i,j,k) .gt. 258.) then
1353            LH = k
1354            goto 15
1355           endif
1356          end do
1357      15  continue
1358          if(LH .lt. 1) goto 100
1359     
1360     !    convert dqcond from kg water/kg air/s to kg water/m3/s and reverse
1361     !    sign so that dqcond < 0. (positive precip) means qls and qcv > 0.
1362          do k = LH, km
1363           qls(k) = -dqcond(i,j,k)*pls/pac*rhoa(i,j,k)
1364     !      qcv(k) = -dqcond(i,j,k)*pcv/pac*rhoa(i,j,k)
1365          end do
1366     
1367     !    Loop over vertical to do the scavenging!
1368          do k = LH, km
1369     
1370     !-----------------------------------------------------------------------------
1371     !   (1) LARGE-SCALE RAINOUT:             
1372     !       Tracer loss by rainout = TC0 * F * exp(-B*dt)
1373     !         where B = precipitation frequency,
1374     !               F = fraction of grid box covered by precipitating clouds.
1375     !       We assume that tracer scavenged by rain is falling down to the
1376     !       next level, where a fraction could be re-evaporated to gas phase
1377     !       if Qls is less then 0 in that level.
1378     !-----------------------------------------------------------------------------
1379           if (qls(k) .gt. 0.) then
1380            F  = F0_ls / (1. + F0_ls*B0_ls*XL_ls/(qls(k)*cdt/Td_ls))
1381            B  = B0_ls/F0_ls +1./(F0_ls*XL_ls/qls(k)) 
1382            BT = B * Td_ls
1383            if (BT.gt.10.) BT = 10.               !< Avoid overflow >
1384     !      Adjust du level:
1385            do n = nHydrophilic, nHydrophilic
1386             DC(n) = w_c%qa(nbeg+n-1)%data3d(i,j,k) * F * (1.-exp(-BT))
1387             if (DC(n).lt.0.) DC(n) = 0.
1388             w_c%qa(nbeg+n-1)%data3d(i,j,k) = w_c%qa(nbeg+n-1)%data3d(i,j,k)-DC(n)
1389             w_c%qa(nbeg+n-1)%data3d(i,j,k) = max(w_c%qa(nbeg+n-1)%data3d(i,j,k),1.e-32)
1390            end do
1391     !      Flux down:  unit is kg m-2
1392     !      Formulated in terms of production in the layer.  In the revaporation step
1393     !      we consider possibly adding flux from above...
1394            do n = nHydrophilic, nHydrophilic
1395             Fd(k,n) = DC(n) * pdog(i,j,k)
1396            end do
1397     
1398           end if                                    ! if Qls > 0  >>>
1399     
1400     !-----------------------------------------------------------------------------
1401     ! * (2) LARGE-SCALE WASHOUT:
1402     ! *     Occurs when rain at this level is less than above.
1403     !-----------------------------------------------------------------------------
1404           if(k .gt. LH .and. qls(k) .ge. 0.) then
1405            if(qls(k) .lt. qls(k-1)) then
1406     !       Find a maximum F overhead until the level where Qls<0.
1407             Qmx   = 0.
1408             do kk = k-1,LH,-1
1409              if (Qls(kk).gt.0.) then
1410               Qmx = max(Qmx,Qls(kk))
1411              else
1412               goto 333
1413              end if
1414             end do
1415     
1416      333    continue
1417             F = F0_ls / (1. + F0_ls*B0_ls*XL_ls/(Qmx*cdt/Td_ls))
1418             if (F.lt.0.01) F = 0.01
1419     !-----------------------------------------------------------------------------
1420     !  The following is to convert Q(k) from kgH2O/m3/sec to mm/sec in order
1421     !  to use the Harvard formula.  Convert back to mixing ratio by multiplying
1422     !  by rhoa.  Multiply by pdog gives kg/m2/s of precip.  Divide by density
1423     !  of water (=1000 kg/m3) gives m/s of precip and multiply by 1000 gives
1424     !  units of mm/s (omit the multiply and divide by 1000).
1425     !-----------------------------------------------------------------------------
1426     
1427             Qd = Qmx /rhoa(i,j,k)*pdog(i,j,k)
1428             if (Qd.ge.50.) then
1429              B = 0.
1430             else
1431              B = Qd * 0.1
1432             end if
1433             BT = B * cdt
1434             if (BT.gt.10.) BT = 10.
1435     
1436     !       Adjust du level:
1437             do n = nHydrophilic, nHydrophilic
1438              DC(n) = w_c%qa(nbeg+n-1)%data3d(i,j,k) * F * (1.-exp(-BT))
1439              if (DC(n).lt.0.) DC(n) = 0.
1440              w_c%qa(nbeg+n-1)%data3d(i,j,k) = w_c%qa(nbeg+n-1)%data3d(i,j,k)-DC(n)
1441              w_c%qa(nbeg+n-1)%data3d(i,j,k) = max(w_c%qa(nbeg+n-1)%data3d(i,j,k),1.e-32)
1442              if( associated(fluxout(n)%data2d) ) then
1443               fluxout(n)%data2d(i,j) = fluxout(n)%data2d(i,j)+DC(n)*pdog(i,j,k)/cdt
1444              endif
1445             end do
1446     
1447            end if
1448           end if                                    ! if ls washout  >>>
1449     
1450     !-----------------------------------------------------------------------------
1451     !  (3) CONVECTIVE RAINOUT:
1452     !      Tracer loss by rainout = dd0 * F * exp(-B*dt)
1453     !        where B = precipitation frequency,
1454     !              F = fraction of grid box covered by precipitating clouds.
1455     !-----------------------------------------------------------------------------
1456     
1457           if (qcv(k) .gt. 0.) then
1458            F  = F0_cv / (1. + F0_cv*B0_cv*XL_cv/(Qcv(k)*cdt/Td_cv))
1459            B  = B0_cv
1460            BT = B * Td_cv
1461            if (BT.gt.10.) BT = 10.               !< Avoid overflow >
1462     
1463     !      Adjust du level: 
1464            do n = nHydrophilic, nHydrophilic
1465             DC(n) = w_c%qa(nbeg+n-1)%data3d(i,j,k) * F * (1.-exp(-BT))
1466             if (DC(n).lt.0.) DC(n) = 0.
1467             w_c%qa(nbeg+n-1)%data3d(i,j,k) = w_c%qa(nbeg+n-1)%data3d(i,j,k)-DC(n)
1468             w_c%qa(nbeg+n-1)%data3d(i,j,k) = max(w_c%qa(nbeg+n-1)%data3d(i,j,k),1.e-32)
1469            end do
1470     
1471     !------  Flux down:  unit is kg. Including both ls and cv.
1472            do n = nHydrophilic, nHydrophilic
1473             Fd(k,n) = Fd(k,n) + DC(n)*pdog(i,j,k)
1474            end do
1475     
1476           end if                                  ! if Qcv > 0   >>>
1477     
1478     !-----------------------------------------------------------------------------
1479     !  (4) CONVECTIVE WASHOUT:
1480     !      Occurs when rain at this level is less than above.
1481     !-----------------------------------------------------------------------------
1482     
1483           if (k.gt.LH .and. Qcv(k).ge.0.) then
1484            if (Qcv(k).lt.Qcv(k-1)) then
1485     !-----  Find a maximum F overhead until the level where Qls<0.
1486             Qmx   = 0.
1487             do kk = k-1, LH, -1
1488              if (Qcv(kk).gt.0.) then
1489               Qmx = max(Qmx,Qcv(kk))
1490              else
1491               goto 444
1492              end if
1493             end do
1494     
1495      444    continue
1496             F = F0_cv / (1. + F0_cv*B0_cv*XL_cv/(Qmx*cdt/Td_cv))
1497             if (F.lt.0.01) F = 0.01
1498     !-----------------------------------------------------------------------------
1499     !  The following is to convert Q(k) from kgH2O/m3/sec to mm/sec in order
1500     !  to use the Harvard formula.  Convert back to mixing ratio by multiplying
1501     !  by rhoa.  Multiply by pdog gives kg/m2/s of precip.  Divide by density
1502     !  of water (=1000 kg/m3) gives m/s of precip and multiply by 1000 gives
1503     !  units of mm/s (omit the multiply and divide by 1000).
1504     !-----------------------------------------------------------------------------
1505     
1506             Qd = Qmx / rhoa(i,j,k)*pdog(i,j,k)
1507             if (Qd.ge.50.) then
1508              B = 0.
1509             else
1510              B = Qd * 0.1
1511             end if
1512             BT = B * cdt
1513             if (BT.gt.10.) BT = 10.
1514     
1515     !       Adjust du level:
1516             do n = nHydrophilic, nHydrophilic
1517              DC(n) = w_c%qa(nbeg+n-1)%data3d(i,j,k) * F * (1.-exp(-BT))
1518              if (DC(n).lt.0.) DC(n) = 0.
1519              w_c%qa(nbeg+n-1)%data3d(i,j,k) = w_c%qa(nbeg+n-1)%data3d(i,j,k)-DC(n)
1520              w_c%qa(nbeg+n-1)%data3d(i,j,k) = max(w_c%qa(nbeg+n-1)%data3d(i,j,k),1.e-32)
1521              if( associated(fluxout(n)%data2d) ) then
1522               fluxout(n)%data2d(i,j) = fluxout(n)%data2d(i,j)+DC(n)*pdog(i,j,k)/cdt
1523              endif
1524             end do
1525     
1526            end if
1527           end if                                    ! if cv washout  >>>
1528     
1529     !-----------------------------------------------------------------------------
1530     !  (5) RE-EVAPORATION.  Assume that SO2 is re-evaporated as SO4 since it
1531     !      has been oxidized by H2O2 at the level above. 
1532     !-----------------------------------------------------------------------------
1533     !     Add in the flux from above, which will be subtracted if reevaporation occurs
1534           if(k .gt. LH) then
1535            do n = 1, nbins
1536             Fd(k,n) = Fd(k,n) + Fd(k-1,n)
1537            end do
1538     
1539     !      Is there evaporation in the currect layer?
1540            if (-dqcond(i,j,k) .lt. 0.) then
1541     !       Fraction evaporated = H2O(k)evap / H2O(next condensation level).
1542             if (-dqcond(i,j,k-1) .gt. 0.) then
1543     
1544               A =  abs(  dqcond(i,j,k) * pdog(i,j,k)    &
1545                 /      ( dqcond(i,j,k-1) * pdog(i,j,k-1))  )
1546               if (A .gt. 1.) A = 1.
1547     
1548     !         Adjust tracer in the level
1549               do n = nHydrophilic, nHydrophilic
1550                DC(n) =  Fd(k-1,n) / pdog(i,j,k) * A
1551                w_c%qa(nbeg+n-1)%data3d(i,j,k) = w_c%qa(nbeg+n-1)%data3d(i,j,k) + DC(n)
1552                w_c%qa(nbeg+n-1)%data3d(i,j,k) = max(w_c%qa(nbeg+n-1)%data3d(i,j,k),1.e-32)
1553     !          Adjust the flux out of the bottom of the layer
1554                Fd(k,n)  = Fd(k,n) - DC(n)*pdog(i,j,k)
1555               end do
1556     
1557             endif
1558            endif                                   ! if -moistq < 0
1559           endif
1560          end do  ! k
1561     
1562          do n = nHydrophilic, nHydrophilic
1563           if( associated(fluxout(n)%data2d) ) then
1564            fluxout(n)%data2d(i,j) = fluxout(n)%data2d(i,j)+Fd(km,n)/cdt
1565           endif
1566          end do
1567     
1568      100 continue
1569         end do   ! i
1570        end do    ! j
1571     
1572        deallocate(fd,DC,stat=ios)
1573     
1574        end subroutine BC_Wet_Removal
1575     
1576     !-------------------------------------------------------------------------
1577     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
1578     !-------------------------------------------------------------------------
1579     !BOP
1580     !
1581     ! !IROUTINE:  BC_Compute_Diags - Calculate dust 2D diagnostics
1582     !
1583     ! !INTERFACE:
1584     !
1585     
1586        subroutine BC_Compute_Diags ( i1, i2, j1, j2, km, nbins, gcBC, w_c, tmpu, &
1587                                      rhoa, sfcmass, colmass, mass, exttau, scatau, &
1588                                      rc )
1589     
1590     ! !USES:
1591     
1592       implicit NONE
1593     
1594     ! !INPUT PARAMETERS:
1595        integer, intent(in) :: i1, i2, j1, j2, km, nbins
1596        type(BC_GridComp), intent(inout):: gcBC     ! BC Grid Component
1597        type(Chem_Bundle), intent(in)   :: w_c      ! Chem Bundle
1598        real, pointer, dimension(:,:,:) :: tmpu     ! temperature [K]
1599        real, pointer, dimension(:,:,:) :: rhoa     ! air density [kg m-3]
1600     
1601     ! !OUTPUT PARAMETERS:
1602        type(Chem_Array), intent(inout)  :: sfcmass ! sfc mass concentration kg/m3
1603        type(Chem_Array), intent(inout)  :: colmass ! col mass density kg/m2
1604        type(Chem_Array), intent(inout)  :: mass    ! 3d mass mixing ratio kg/kg
1605        type(Chem_Array), intent(inout)  :: exttau  ! ext. AOT at 550 nm
1606        type(Chem_Array), intent(inout)  :: scatau  ! sct. AOT at 550 nm
1607        integer, intent(out)             :: rc      ! Error return code:
1608                                                    !  0 - all is well
1609                                                    !  1 - 
1610     
1611     ! !DESCRIPTION: Calculates some simple 2d diagnostics from the BC fields
1612     !               Surface concentration (dry)
1613     !               Column mass load (dry)
1614     !               Extinction aot 550 (wet)
1615     !               Scattering aot 550 (wet)
1616     !               For the moment, this is hardwired.
1617     !
1618     ! !REVISION HISTORY:
1619     !
1620     !  16APR2004, Colarco
1621     !
1622     !EOP
1623     !-------------------------------------------------------------------------
1624     
1625     ! !Local Variables
1626        character(len=*), parameter :: myname = 'BC_Compute_Diags'
1627        integer :: i, j, k, n, nbeg, nend, ios, nch, idx
1628        real :: tau, ssa
1629        character(len=255) :: qname
1630     
1631     
1632     !  Initialize local variables
1633     !  --------------------------
1634        nbeg  = w_c%reg%i_BC
1635        nend  = w_c%reg%j_BC
1636        nch   = gcBC%mie_tables%nch
1637     
1638     !  Calculate the diagnostic variables if requested
1639     !  -----------------------------------------------
1640     
1641     !  Calculate the surface mass concentration
1642        if( associated(sfcmass%data2d) ) then
1643           sfcmass%data2d(i1:i2,j1:j2) = 0.
1644           do n = 1, nbins
1645              sfcmass%data2d(i1:i2,j1:j2) &
1646                   =   sfcmass%data2d(i1:i2,j1:j2) &
1647                   + w_c%qa(n+nbeg-1)%data3d(i1:i2,j1:j2,km)*rhoa(i1:i2,j1:j2,km)
1648           end do
1649        endif
1650     
1651     !  Calculate the dust column loading
1652        if( associated(colmass%data2d) ) then
1653           colmass%data2d(i1:i2,j1:j2) = 0.
1654           do n = 1, nbins
1655            do k = 1, km
1656             colmass%data2d(i1:i2,j1:j2) &
1657              =   colmass%data2d(i1:i2,j1:j2) &
1658                + w_c%qa(n+nbeg-1)%data3d(i1:i2,j1:j2,k)*w_c%delp(i1:i2,j1:j2,k)/grav
1659            end do
1660           end do
1661        endif
1662     
1663     !  Calculate the total mass mixing ratio
1664        if( associated(mass%data3d) ) then
1665           mass%data3d(i1:i2,j1:j2,1:km) = 0.
1666           do n = 1, nbins
1667            mass%data3d(i1:i2,j1:j2,1:km) &
1668              =   mass%data3d(i1:i2,j1:j2,1:km) &
1669                + w_c%qa(n+nbeg-1)%data3d(i1:i2,j1:j2,1:km)
1670           end do
1671        endif
1672     
1673     !  Calculate the extinction and/or scattering AOD
1674        if( associated(exttau%data2d) .or. associated(scatau%data2d) ) then
1675     
1676           if( associated(exttau%data2d) ) then
1677            exttau%data2d(i1:i2,j1:j2) = 0.
1678           endif
1679           if( associated(scatau%data2d) ) then
1680            scatau%data2d(i1:i2,j1:j2) = 0.
1681           endif
1682     
1683           do n = 1, nbins
1684     
1685     !      Select the name for species
1686            qname = trim(w_c%reg%vname(w_c%reg%i_BC+n-1))
1687            idx = Chem_MieQueryIdx(gcBC%mie_tables,qname,rc)
1688            if(rc .ne. 0) call die(myname, 'cannot find proper Mie table index')
1689     
1690     !      Recall -- at present need to divide RH by 100 to get to tables
1691            do k = 1, km
1692             do j = j1, j2
1693              do i = i1, i2
1694               call Chem_MieQuery(gcBC%mie_tables, idx, 1., &
1695                   w_c%qa(nbeg+n-1)%data3d(i,j,k)*w_c%delp(i,j,k)/grav, &
1696                   w_c%rh(i,j,k)/100., tau=tau, ssa=ssa)
1697     
1698     !         Integrate in the vertical
1699               if( associated(exttau%data2d) ) then
1700                exttau%data2d(i,j) = exttau%data2d(i,j) + tau
1701               endif
1702               if( associated(scatau%data2d) ) then
1703                scatau%data2d(i,j) = scatau%data2d(i,j) + tau*ssa
1704               endif
1705     
1706              enddo
1707             enddo
1708            enddo
1709     
1710           enddo  ! nbins
1711     
1712        endif
1713     
1714     
1715        rc = 0
1716     
1717        end subroutine BC_Compute_Diags
1718     
1719      end subroutine BC_GridCompRun
1720     
1721     !-------------------------------------------------------------------------
1722     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
1723     !-------------------------------------------------------------------------
1724     !BOP
1725     !
1726     ! !IROUTINE:  BC_GridCompFinalize --- The Chem Driver 
1727     !
1728     ! !INTERFACE:
1729     !
1730     
1731        subroutine BC_GridCompFinalize ( gcBC, w_c, impChem, expChem, &
1732                                         nymd, nhms, cdt, rc )
1733     
1734     ! !USES:
1735     
1736       implicit NONE
1737     
1738     ! !INPUT/OUTPUT PARAMETERS:
1739     
1740        type(BC_GridComp), intent(inout) :: gcBC   ! Grid Component
1741     
1742     ! !INPUT PARAMETERS:
1743     
1744        type(Chem_Bundle), intent(in)  :: w_c      ! Chemical tracer fields   
1745        integer, intent(in) :: nymd, nhms          ! time
1746        real,    intent(in) :: cdt                 ! chemical timestep (secs)
1747     
1748     
1749     ! !OUTPUT PARAMETERS:
1750     
1751        type(ESMF_State), intent(inout) :: impChem   ! Import State
1752        type(ESMF_State), intent(inout) :: expChem   ! Import State
1753        integer, intent(out) ::  rc                  ! Error return code:
1754                                                     !  0 - all is well
1755                                                     !  1 -
1756      
1757     ! !DESCRIPTION: This routine finalizes this Grid Component.
1758     !
1759     ! !REVISION HISTORY:
1760     !
1761     !  18Sep2003 da Silva  First crack.
1762     !
1763     !EOP
1764     !-------------------------------------------------------------------------
1765     
1766        character(len=*), parameter :: myname = 'BC_GridCompFinalize'
1767        rc=0
1768        return
1769     
1770      end subroutine BC_GridCompFinalize
1771     
1772      end module BC_GridCompMod
1773     
1774