File: C:\NOAA\NEMS_11731\src\chem\gocart\src\GMAO_Shared\Chem_Base\Chem_Aod3d.F90

1     ! Colarco, June 2006
2     ! Offline calculator of 3D aerosol optical properties
3     ! 1) Read a registry determining aerosol optical properties to compute
4     ! 2) Read a chem bundle of aerosol distribution from a file
5     ! 3) Read an output registry naming the output chem bundle
6     ! 4) Compute
7     
8       program chem_aodcalc
9     
10       use m_die, only: die
11       use Chem_MieMod
12       use Chem_RegistryMod
13       use Chem_BundleMod
14     
15       implicit none
16     
17       real, parameter :: grav = 9.80616
18       real, parameter :: rgas = 287.04
19     
20       character(len=*), parameter :: myname = 'chem_aod'
21       type(Chem_Mie)      :: mie_tables
22       type(Chem_Registry) :: regInp      ! chemistry registry
23       type(Chem_Registry) :: regOut      ! chemistry registry
24       type(Chem_Bundle)   :: w_c         ! input aerosol chemistry bundle 
25       type(Chem_Bundle)   :: w_tau       ! output total tau chemistry bundle
26       type(Chem_Bundle)   :: w_taudu     ! output du tau chemistry bundle
27       type(Chem_Bundle)   :: w_tauant    ! output anthro tau chemistry bundle
28       type(Chem_Bundle)   :: w_tauss     ! output seasalt tau chemistry bundle
29       type(Chem_Bundle)   :: w_tauoc     ! output organic carbon tau chemistry bundle
30       type(Chem_Bundle)   :: w_taubc     ! output black carbon tau chemistry bundle
31       type(Chem_Bundle)   :: w_taucc     ! output total carbon tau chemistry bundle
32       type(Chem_Bundle)   :: w_tausu     ! output sulfate tau chemistry bundle
33       integer :: i, j, k, im, jm, km, idx, n
34       integer :: i1, i2, ig, j1, j2, jg, ik, iq, iz, kk
35       integer :: nymd, nhms, timidx, freq, rc, ier, fid
36       integer :: idxTable
37       integer :: iarg, iargc, argc, lenfile
38       logical :: doing_dust, doing_anthro, doing_ss, doing_oc, doing_bc, doing_su, doing_cc
39       logical :: doing_geos4
40       logical :: doing_dry
41       logical :: new
42       logical :: found_airdensfile, found_hghtefile
43       real :: channel, tau_, ssa_, bbck_, bext_, taulev, gasym_, scaleRH, maxRH, qMass
44       integer :: itau, iext, issa, immr, ibbck, iabck0, iabck1, ietob, igasym
45       integer, parameter :: READ_ONLY = 1
46       real, pointer :: delz(:,:,:), t(:,:,:), q(:,:,:), hghte(:,:,:)
47       real, pointer :: airdens(:,:,:) => null()
48       character(len=255) :: infile, outfile, filename, airdensfile, hghtefile, rcfile, argv
49       character(len=8)   :: datestr
50     
51     ! Parse the command line (see usage() below)
52       argc = iargc()
53       if(argc .lt. 1) call usage()
54       iarg = 0
55       outfile = 'tau3d.nc4'
56       rcfile  = 'Aod_CALIPSO.rc'
57       doing_dust   = .false.
58       doing_ss     = .false.
59       doing_bc     = .false.
60       doing_oc     = .false.
61       doing_cc     = .false.
62       doing_su     = .false.
63       doing_anthro = .false.
64       doing_geos4  = .false.
65       doing_dry    = .false.
66       found_airdensfile = .false.
67       found_hghtefile = .false.
68       do i = 0, 32767
69        iarg = iarg+1
70        if(iarg .gt. argc) exit
71        call GetArg(iarg, argv)
72        select case(argv)
73         case ("-geos4")
74          doing_geos4 = .true.
75         case ("-dryaer")
76          doing_dry = .true.
77         case ("-dust")
78          doing_dust = .true.
79         case ("-ss")
80          doing_ss = .true.
81         case ("-su")
82          doing_su = .true.
83         case ("-bc")
84          doing_bc = .true.
85         case ("-oc")
86          doing_oc = .true.
87         case ("-cc")
88          doing_cc = .true.
89         case ("-anthro")
90          doing_anthro = .true.
91         case ("-airdensfile")
92          if(iarg+1 .gt. argc) call usage()
93          iarg = iarg+1
94          call GetArg(iarg, airdensfile)
95          found_airdensfile = .true.
96         case ("-hghtefile")
97          if(iarg+1 .gt. argc) call usage()
98          iarg = iarg+1
99          call GetArg(iarg, hghtefile)
100          found_hghtefile = .true.
101         case ("-o")
102          if(iarg+1 .gt. argc) call usage()
103          iarg = iarg+1
104          call GetArg(iarg, outfile)
105         case ("-t")
106          if(iarg+1 .gt. argc) call usage()
107          iarg = iarg+1
108          call GetArg(iarg, rcfile)
109         case default
110          infile = argv
111        end select
112       end do
113       rcfile = trim(rcfile)
114       infile = trim(infile)
115       outfile = trim(outfile)
116       lenfile = len(trim(outfile))
117       if(found_hghtefile) hghtefile = trim(hghtefile)
118       if(found_airdensfile) then
119        airdensfile = trim(airdensfile)
120       else
121        airdensfile = infile
122       endif
123     
124     ! Scaling of Relative Humidity
125     ! Input optics files (e.g., optics_XX.nc4) have a fractional RH
126     ! coordinate (that is, RH varies 0 - 1 in the file, as in GEOS-5).
127     ! In GEOS-4, however, RH in the chem.eta file is represented as
128     ! a percentage, so it varies 0 - 100%.  For compatibility we
129     ! introduce a flag "-geos4" on execution.  If present, the RH
130     ! from the chem.eta file is divided by 100 on input to the Mie
131     ! calculator.
132     ! If you requested to do the calculation like the aerosols were
133     ! dry, then we set the RH like it is 0%
134       scaleRH = 1.
135       if(doing_geos4) scaleRH = 1. / 100.
136       if(doing_dry)   scaleRH = 0.
137     
138     ! Hardwired: Create the output bundle
139     ! -----------------------------------
140       regOut = Chem_RegistryCreate(ier, rcfile)
141       if(ier /= 0) call die(myname, 'cannot create output registry')
142     
143     ! Check the output registry for desired output fields
144     ! ---------------------------------------------------
145       itau   = -1
146       iext   = -1
147       issa   = -1
148       immr   = -1
149       ietob  = -1
150       ibbck  = -1
151       iabck0 = -1
152       iabck1 = -1
153       igasym = -1
154       do iq = 1, regOut%nq
155        if(trim(regOut%vname(iq)) .eq. 'tau')        itau   = iq
156        if(trim(regOut%vname(iq)) .eq. 'extinction') iext   = iq
157        if(trim(regOut%vname(iq)) .eq. 'ssa')        issa   = iq
158        if(trim(regOut%vname(iq)) .eq. 'mmr')        immr   = iq
159        if(trim(regOut%vname(iq)) .eq. 'ext2back')   ietob  = iq
160        if(trim(regOut%vname(iq)) .eq. 'attback0')   iabck0 = iq
161        if(trim(regOut%vname(iq)) .eq. 'aback_sfc')  iabck0 = iq ! alias
162        if(trim(regOut%vname(iq)) .eq. 'attback1')   iabck1 = iq
163        if(trim(regOut%vname(iq)) .eq. 'aback_toa')  iabck1 = iq ! alias
164        if(trim(regOut%vname(iq)) .eq. 'backscat')   ibbck  = iq
165        if(trim(regOut%vname(iq)) .eq. 'gasym')      igasym = iq
166       enddo
167     
168     ! At this point should do some checking that certain fields are
169     ! provided.
170       if(issa .gt. 0   .and. itau .lt. 0)                      ier = 1
171       if(igasym .gt. 0 .and. (itau .lt. 0 .or. issa .lt. 0))   ier = 2
172       if(ibbck .gt. 0  .and. itau .lt. 0)                      ier = 3
173       if(ietob .gt. 0  .and. (itau .lt. 0 .or. ibbck .lt. 0))  ier = 4
174       if((iabck0 .gt. 0 .or. iabck1 .gt. 0) .and. &
175          (itau .lt. 0 .or.ibbck .lt. 0) )                      ier = 5
176     
177       if(ier /= 0) then
178          print *,'----------------------------------------------------'
179          print *,'ier = ', ier
180          print *,'issa, itau, igasym, ibbck, ietob, iabck0, iabck1 = ', &
181                 issa, itau, igasym, ibbck, ietob, iabck0, iabck1
182          call die(myname,'inconsistency in output registry')
183       end if
184     
185     ! Hardwired: Read the input chemistry registry from Chem_Registry.rc
186     ! This registry file describes the input chemistry bundle being 
187     ! operated on.
188     ! ------------------------------------------------------------------
189       regInp = Chem_RegistryCreate(ier,'Chem_MieRegistry.rc')
190       if(ier /= 0) call die(myname, 'cannot create registry')
191     
192     ! Hardwired: use the Chem_MieMod function to create the Mie tables
193     ! ------------------------------------------------------------------
194       mie_tables = Chem_MieCreate(rcfile,ier)
195     
196     ! This part needs some work: in GEOS-4 the chemistry bundle input 
197     ! files generally had 4 times per file, and in GEOS-5 they generally
198     ! have 1 time per file.  The loop over IDX below is the number of
199     ! times on the file.  In the future this should come from the
200     ! chemistry bundle and not have to be hardwired.
201     ! ------------------------------------------------------------------
202       new = .true.
203       do idx = 1, 1
204     
205     !  Read the input chemistry bundle from the infile
206     !  -------------------------------------------------
207        call Chem_BundleRead(infile, nymd, nhms, w_c, rc, freq=freq, &
208                             ChemReg=regInp, timidx=idx)
209     
210        i1 = w_c%grid%i1
211        i2 = w_c%grid%i2
212        ig = w_c%grid%ig
213        j1 = w_c%grid%j1
214        j2 = w_c%grid%j2
215        jg = w_c%grid%jg
216        im = w_c%grid%im
217        jm = w_c%grid%jm
218        km = w_c%grid%km
219     
220     !  Layer thickness in cartesian sense is needed for certain extrinsic
221     !  properties, specifically any calculations requiring backscatter
222     !  or extinction.  If these fields are requested then we need to allocate
223     !  space for the layer thickness ("delz").  The value of delz is
224     !  filled in either by (1) reading the layer edge heights from a file (the
225     !  "hghtefile") or (2) by reading the air density from a file (the
226     !  "airdensfile") and computing via hydrostatic relation 
227     !  (i.e., dz = dp/g/rhoa) or (3) by computing directly from provided 
228     !  thermodynamic variables (e.g., t, qv).
229        if(iext .gt. 0 .or. ibbck .gt. 0) then 
230     
231     !   Allocate space to contain the thickness of the layer
232         allocate(delz(i1:i2,j1:j2,1:km), stat=ier)
233         if(ier /= 0) call die(myname,'could not allocate space for delz')
234     
235     !   If "hghtefile" is provided let's get delz that way
236         if(found_hghtefile) then
237          allocate(hghte(i1:i2,j1:j2,1:km+1), stat=ier)
238          if(ier /= 0) call die(myname,'could not allocate space for HGHTE')
239          call GFIO_Open ( hghtefile, READ_ONLY, fid, ier )
240          if(ier /= 0) call die(myname,'could not open '//hghtefile//' to read HGHTE')
241          call GFIO_GetVar ( fid, 'HGHTE', nymd, nhms,     &
242                             im, jm, 1, km+1, hghte, ier )
243          if(ier /= 0) call die(myname,'could not read HGHTE from '//hghtefile)
244          call GFIO_Close ( fid, ier)
245          if(ier /= 0) call die(myname,'could not close infile for HGHTE read')
246          do ik = 1, km
247           delz(:,:,ik) = hghte(:,:,ik)-hghte(:,:,ik+1)
248          end do
249          deallocate(hghte,stat=ier)
250          if(ier /= 0) call die(myname,'could not deallocate space for HGHTE')
251     
252         else
253     !    Get the air density, either from aerosol file or externally
254          if(trim(airdensfile) .eq. trim(infile)) then
255           airdens => delz
256     !     Get the air density from the aerosol file
257           call GFIO_Open ( airdensfile, READ_ONLY, fid, ier )
258           if(ier /= 0) call die(myname,'could not open '//airdensfile//' to read airdens')
259           call GFIO_GetVar ( fid, 'AIRDENS', nymd, nhms,     &
260                              im, jm, 1, km, airdens, ier )
261           if(ier /= 0) then
262              call GFIO_GetVar ( fid, 'airdens', nymd, nhms,     &
263                                 im, jm, 1, km, airdens, ier )
264              if(ier /= 0) call die(myname,'could not read airdens from '//airdensfile)
265           end if
266           call GFIO_Close ( fid, ier)
267           if(ier /= 0) call die(myname,'could not close infile for airdens read')
268          else
269     !     Get air density from an external file
270           call GFIO_Open ( airdensfile, READ_ONLY, fid, ier )
271           if(ier /= 0) call die(myname,'could not open '//airdensfile//' to read airdens')
272           call GFIO_GetVar ( fid, 'AIRDENS', nymd, nhms,     &
273                              im, jm, 1, km, delz, ier )
274           if ( ier /= 0 ) then
275                 call GFIO_GetVar ( fid, 'airdens', nymd, nhms,     &
276                                    im, jm, 1, km, delz, ier )
277           end if
278           if(ier .eq. 0) then
279            airdens => delz
280            call GFIO_Close ( fid, ier)
281            if(ier /= 0) call die(myname,'could not close infile for airdens read')
282           else
283     !      If the air density is not on the file try to compute delz directly
284     !      from the thermodynamic variables (need T, QV)
285            allocate(t(i1:i2,j1:j2,1:km), q(i1:i2,j1:j2,1:km), stat=ier)
286            if(ier /= 0) call die(myname,'could not allocate space to compute airdens')
287            call GFIO_GetVar ( fid, 'T', nymd, nhms,     &
288                               im, jm, 1, km, t, ier )
289            if(ier /= 0) call die(myname,'could not read temperature '//airdensfile)
290            call GFIO_GetVar ( fid, 'QV', nymd, nhms,     &
291                               im, jm, 1, km, q, ier )
292            if(ier /= 0) call die(myname,'could not read specific humidity from '//airdensfile)
293            call GFIO_Close ( fid, ier)
294            if(ier /= 0) call die(myname,'could not close infile for airdens read')
295            do j = j1, j2
296             do i = i1, i2
297              call delz_
298             enddo
299            enddo
300            deallocate(t,q,stat=ier)
301            if(ier /= 0) call die(myname,'cloud not deallocate t,q')
302           endif
303          endif
304     
305     !    If you actually read air density, compute dz
306          if(associated(airdens)) then
307           delz = w_c%delp/grav/airdens
308          endif
309         endif
310     
311     !   Turn delz from m to km
312         delz = delz/1000.
313     
314        endif
315     
316        print *, 'Computing AOD for ', nymd, nhms
317     
318     !  Check that the RH seems sane
319        if(.not. doing_geos4 .and. .not. doing_dry) then
320         maxrh = maxval(w_c%rh)
321         if(maxrh .gt. 2.) then
322          print *, 'Maximum RH value = ', maxRH
323          print *, 'Should you have chosen "-geos4" as a command line option?'
324         endif
325        endif
326     
327     !  ==================================================================================
328     !  The enclosed bundle of code selects on what calculation we run and what is written
329     
330     !  Create and initialize output chemistry bundle
331        call create_species_bundle(w_tau)
332     
333     !  If doing species
334        if(doing_dust) call create_species_bundle(w_taudu)
335        if(doing_anthro) call create_species_bundle(w_tauant)
336        if(doing_ss) call create_species_bundle(w_tauss)
337        if(doing_su) call create_species_bundle(w_tausu)
338        if(doing_oc) call create_species_bundle(w_tauoc)
339        if(doing_bc) call create_species_bundle(w_taubc)
340        if(doing_cc) call create_species_bundle(w_taucc)
341     
342        do ik = 1, 1
343         channel = mie_tables%channels(ik)
344     
345         do iq = 1, mie_tables%nq
346     
347     !    Sanity check: Does the mie_table name = chem_registry name?
348          if(trim(mie_tables%vname(iq)) .ne. trim(w_c%reg%vname(iq)) ) &
349            call die(myname, 'mie_tables and chem_registry vname mismatch')
350     !    Cases
351          idxTable = Chem_MieQueryIdx(mie_tables,mie_tables%vname(iq),rc)
352     
353          if(idxTable .ne. -1) then
354           do k = 1, km
355           do j = 1, jm
356           do i = 1, im
357     
358     !     fix in case NaN on input of qa
359     #ifndef sysAIX
360           if(isnan(w_c%qa(iq)%data3d(i,j,k))) &
361             w_c%qa(iq)%data3d(i,j,k) = tiny(w_c%qa(iq)%data3d(i,j,k))
362     #endif
363           qMass = w_c%qa(iq)%data3d(i,j,k)*w_c%delp(i,j,k)/9.81
364           call Chem_MieQuery(mie_tables, idxTable, 1.*ik, &
365                              qMass, &
366                              w_c%rh(i,j,k) * scaleRH, tau=tau_, ssa=ssa_, &
367                              bbck=bbck_, bext=bext_, gasym=gasym_)
368     
369     !     Fill in the total values
370     !     Note the weighting of the ssa, backscatter, and e_to_b ratio
371           call fill(w_tau)
372     
373           if(doing_dust .and. (iq .ge. w_c%reg%i_du .and. iq .le. w_c%reg%j_du)) then
374            call fill(w_taudu)
375           endif
376     
377           if(doing_su .and. (iq .ge. w_c%reg%i_su .and. iq .le. w_c%reg%j_su)) then
378            call fill(w_tausu)
379           endif
380     
381           if(doing_bc .and. (iq .ge. w_c%reg%i_bc .and. iq .le. w_c%reg%j_bc)) then
382            call fill(w_taubc)
383           endif
384     
385           if(doing_oc .and. (iq .ge. w_c%reg%i_oc .and. iq .le. w_c%reg%j_oc)) then
386            call fill(w_tauoc)
387           endif
388     
389           if(doing_ss .and. (iq .ge. w_c%reg%i_ss .and. iq .le. w_c%reg%j_ss)) then
390            call fill(w_tauss)
391           endif
392     
393           if(doing_anthro .and. &
394              ( (iq .ge. w_c%reg%i_oc .and. iq .le. w_c%reg%j_oc) .or. &
395                (iq .ge. w_c%reg%i_bc .and. iq .le. w_c%reg%j_bc) .or. &
396                (iq .ge. w_c%reg%i_su .and. iq .le. w_c%reg%j_su) )      ) then
397            call fill(w_tauant)
398           endif
399     
400           if(doing_cc .and. &
401              ( (iq .ge. w_c%reg%i_oc .and. iq .le. w_c%reg%j_oc) .or. &
402                (iq .ge. w_c%reg%i_bc .and. iq .le. w_c%reg%j_bc) )      ) then
403            call fill(w_taucc)
404           endif
405     
406           enddo  ! i
407           enddo  ! j
408           enddo  ! k
409          endif   
410         enddo    ! iq
411     
412     call Chem_RegistryPrint(regout)
413     
414     !   Normalize the ssa calculation
415         call normal(w_tau)
416         if(doing_dust)   call normal(w_taudu)
417         if(doing_ss)     call normal(w_tauss)
418         if(doing_su)     call normal(w_tausu)
419         if(doing_bc)     call normal(w_taubc)
420         if(doing_oc)     call normal(w_tauoc)
421         if(doing_cc)     call normal(w_taucc)
422         if(doing_anthro) call normal(w_tauant)
423     
424     !   Compute the attentuated backscatter terms
425         call backscatter(w_tau)
426         if(doing_dust) call backscatter(w_taudu)
427         if(doing_ss)   call backscatter(w_tauss)
428         if(doing_su)   call backscatter(w_tausu)
429         if(doing_oc)   call backscatter(w_tauoc)
430         if(doing_bc)   call backscatter(w_taubc)
431         if(doing_cc)   call backscatter(w_taucc)
432         if(doing_anthro) call backscatter(w_tauant)
433     
434        enddo
435     
436     !  Clean up allocation of space for delz if doing extinction
437        if(iext .gt. 0 .or. ibbck .gt. 0) deallocate(delz)
438     
439     !  Write the Chem_Bundle out
440        write(datestr,'(i8.8)') nymd
441        filename = trim(outfile(1:lenfile))
442        call Chem_BundleWrite( filename, nymd, nhms, 0, w_tau, rc, &
443                               verbose=.true., new=new)
444        if(doing_dust) then 
445         filename = trim(outfile(1:lenfile)//'.dust')
446         call Chem_BundleWrite( filename, nymd, nhms, 0, w_taudu, rc, &
447                                verbose=.true., new=new)
448        endif
449        if(doing_ss) then 
450         filename = trim(outfile(1:lenfile)//'.ss')
451         call Chem_BundleWrite( filename, nymd, nhms, 0, w_tauss, rc, &
452                                verbose=.true., new=new)
453        endif
454        if(doing_su) then 
455         filename = trim(outfile(1:lenfile)//'.su')
456         call Chem_BundleWrite( filename, nymd, nhms, 0, w_tausu, rc, &
457                                verbose=.true., new=new)
458        endif
459        if(doing_oc) then 
460         filename = trim(outfile(1:lenfile)//'.oc')
461         call Chem_BundleWrite( filename, nymd, nhms, 0, w_tauoc, rc, &
462                                verbose=.true., new=new)
463        endif
464        if(doing_bc) then 
465         filename = trim(outfile(1:lenfile)//'.bc')
466         call Chem_BundleWrite( filename, nymd, nhms, 0, w_taubc, rc, &
467                                verbose=.true., new=new)
468        endif
469        if(doing_cc) then 
470         filename = trim(outfile(1:lenfile)//'.cc')
471         call Chem_BundleWrite( filename, nymd, nhms, 0, w_taucc, rc, &
472                                verbose=.true., new=new)
473        endif
474        if(doing_anthro) then 
475         filename = trim(outfile(1:lenfile)//'.anthro')
476         call Chem_BundleWrite( filename, nymd, nhms, 0, w_tauant, rc, &
477                                verbose=.true., new=new)
478        endif
479     
480     
481     !  ==================================================================================
482     
483     !  Don't overwrite the file
484     !  ------------------------
485        new = .false.
486     
487     !  Clean up pointers
488     !  -----------------
489        call Chem_BundleDestroy(w_tau,rc)
490        if(doing_anthro) call Chem_BundleDestroy(w_tauant,rc)
491        if(doing_dust)   call Chem_BundleDestroy(w_taudu,rc)
492        if(doing_ss)     call Chem_BundleDestroy(w_tauss,rc)
493        if(doing_su)     call Chem_BundleDestroy(w_tausu,rc)
494        if(doing_oc)     call Chem_BundleDestroy(w_tauoc,rc)
495        if(doing_bc)     call Chem_BundleDestroy(w_taubc,rc)
496        if(doing_cc)     call Chem_BundleDestroy(w_taucc,rc)
497     
498        call Chem_BundleDestroy(w_c,rc)
499     
500       enddo   ! idx (time increment in input file)
501     
502     ! Destroy Mie tables
503       call Chem_MieDestroy(mie_tables,ier)
504       call Chem_RegistryDestroy(regInp, rc)
505       call Chem_RegistryDestroy(regOut, rc)
506     
507     ! ----------------------------------------------------------------------------
508       contains
509     
510       subroutine usage()
511       print *
512       print *,'Usage: '
513       print *,'  Chem_Aod3d.x [-dust -anthro '
514       print *,'                -o outfile -t rcfile ] infile'
515       print *
516       print *, 'where'
517       print *
518       print *, '-geos4       to specify that the relative humidity of input file'
519       print *, '             varies 0 - 100 instead of 0 - 1 as in GEOS-5'
520       print *, '-dryaer      to specify to ignore the relative humidity in the'
521       print *, '             input file; compute all properties like RH = 0%'
522       print *, '-dust        compute for dust separately'
523       print *, '-ss          compute for seasalt separately'
524       print *, '-su          compute for sulfate separately'
525       print *, '-bc          compute for black carbon separately'
526       print *, '-oc          compute for organic carbon separately'
527       print *, '-cc          compute for total carbon separately'
528       print *, '-anthro      compute for SU, BC, and OC separately'
529       print *, '-o outfile   output file header'
530       print *, '-t rcfile    resource file specifying channels for AOD calc'
531       print *, '-airdensfile airdensfilename   filename providing variable AIRDENS'
532       print *, '                               if omitted, use infile'
533       print *, '-hghtefile   hghtefilename     filename providing variable HGHTE'
534       print *, '                               supercedes airdensfile if both provided'
535       print *, 'infile       mandatory c_rst file'
536       print *
537       call exit(1)
538       end subroutine usage
539     
540     
541       subroutine create_species_bundle(this)
542       type(Chem_Bundle) :: this
543       integer i
544       call Chem_BundleCreate(regOut, &
545                              i1, i2, ig, im, &
546                              j1, j2, jg, jm, km, &
547                              this, ier)
548       if(ier /= 0) call die(myname, 'cannot create bundle')
549       this%delp = w_c%delp
550       this%rh   = w_c%rh
551       do n = 1, regOut%nq
552        this%qa(n)%data3d = 0.0
553       end do
554       end subroutine create_species_bundle
555     
556     
557       subroutine fill(this)
558       type(Chem_Bundle) :: this
559           if(itau .gt. 0)   this%qa(itau)%data3d(i,j,k) &
560                              = this%qa(itau)%data3d(i,j,k) + tau_
561           if(iext .gt. 0)   this%qa(iext)%data3d(i,j,k) &
562                              = this%qa(iext)%data3d(i,j,k) + tau_/delz(i,j,k)
563           if(issa .gt. 0)   this%qa(issa)%data3d(i,j,k) &
564                              = this%qa(issa)%data3d(i,j,k) + ssa_*tau_
565           if(igasym .gt. 0) this%qa(igasym)%data3d(i,j,k) &
566                              = this%qa(igasym)%data3d(i,j,k) + gasym_*ssa_*tau_
567           if(immr .gt. 0)   this%qa(immr)%data3d(i,j,k) &
568                              = this%qa(immr)%data3d(i,j,k) + w_c%qa(iq)%data3d(i,j,k)
569           if(ibbck .gt. 0)  this%qa(ibbck)%data3d(i,j,k) &
570                              = this%qa(ibbck)%data3d(i,j,k) + bbck_*qMass/delz(i,j,k)
571       end subroutine fill
572     
573     
574       subroutine normal(this)
575       type(Chem_Bundle) :: this
576     !     Check floating underflow in tau
577           do k = 1, km
578            do j = 1, jm
579             do i = 1, im
580              if(itau .gt. 0) this%qa(itau)%data3d(i,j,k) = max(     this%qa(itau)%data3d(i,j,k), &
581                                                                tiny(this%qa(itau)%data3d(i,j,k)))
582              if(iext .gt. 0) this%qa(iext)%data3d(i,j,k) = max(     this%qa(iext)%data3d(i,j,k), &
583                                                                tiny(this%qa(iext)%data3d(i,j,k)))
584              if(ibbck .gt. 0) this%qa(ibbck)%data3d(i,j,k) = max(     this%qa(ibbck)%data3d(i,j,k), &
585                                                                tiny(this%qa(ibbck)%data3d(i,j,k)))
586              if(issa .gt. 0) this%qa(issa)%data3d(i,j,k) = max(     this%qa(issa)%data3d(i,j,k), &
587                                                                tiny(this%qa(issa)%data3d(i,j,k)))
588              if(issa .gt. 0)   this%qa(issa)%data3d(i,j,k) = this%qa(issa)%data3d(i,j,k)/this%qa(itau)%data3d(i,j,k)
589              if(igasym .gt. 0) this%qa(igasym)%data3d(i,j,k) = max(     this%qa(igasym)%data3d(i,j,k), &
590                                                                tiny(this%qa(igasym)%data3d(i,j,k)))
591              if(igasym .gt. 0) this%qa(igasym)%data3d(i,j,k) = this%qa(igasym)%data3d(i,j,k) &
592                                                         /(this%qa(issa)%data3d(i,j,k)*this%qa(itau)%data3d(i,j,k))
593              if(ietob .gt. 0)  this%qa(ietob)%data3d(i,j,k) = this%qa(itau)%data3d(i,j,k)/delz(i,j,k)/this%qa(ibbck)%data3d(i,j,k)
594              if(ietob .gt. 0) this%qa(ietob)%data3d(i,j,k) = max(     this%qa(ietob)%data3d(i,j,k), &
595                                                                tiny(this%qa(ietob)%data3d(i,j,k)))
596             enddo
597            enddo
598           enddo
599       end subroutine normal
600     
601     
602       subroutine backscatter(this)
603       type(Chem_Bundle) :: this
604     
605     !  Attenuated backscatter from space (I think, right, k = 1 is TOA)
606        if(iabck1 .gt. 0) then 
607         this%qa(iabck1)%data3d(:,:,1) = this%qa(ibbck)%data3d(:,:,1)*exp(-this%qa(itau)%data3d(:,:,1))
608         do k = 2, km
609          do j = 1, jm
610           do i = 1, im
611            taulev = 0.
612            do kk = 1, k-1
613             taulev = taulev + this%qa(itau)%data3d(i,j,kk)
614            enddo
615            taulev = taulev + 0.5 * this%qa(itau)%data3d(i,j,k)
616            this%qa(iabck1)%data3d(i,j,k) = this%qa(ibbck)%data3d(i,j,k)*exp(-2.*taulev)
617           enddo
618          enddo
619         enddo
620        endif
621     
622     !  Attenuated backscatter from surface (I think, right, k = 1 is TOA)
623        if(iabck0 .gt. 0) then 
624         this%qa(iabck0)%data3d(:,:,km) = this%qa(ibbck)%data3d(:,:,km)*exp(-this%qa(itau)%data3d(:,:,km))
625         do k = km-1, 1, -1
626          do j = 1, jm
627           do i = 1, im
628            taulev = 0.
629            do kk = km, k+1, -1
630             taulev = taulev + this%qa(itau)%data3d(i,j,kk)
631            enddo
632            taulev = taulev + 0.5 * this%qa(itau)%data3d(i,j,k)
633            this%qa(iabck0)%data3d(i,j,k) = this%qa(ibbck)%data3d(i,j,k)*exp(-2.*taulev)
634           enddo
635          enddo
636         enddo
637        endif
638       end subroutine backscatter
639     
640     
641       subroutine delz_
642        implicit none
643     
644        real :: Tv, mixr, pe(km+1)
645        real, parameter :: ptop = 1.
646     !  Construct edge pressures
647     !  ------------------------
648        pe(1) = ptop
649        do k = 2, km + 1
650           pe(k) = pe(k-1) + w_c%delp(i,j,k-1)
651        end do
652     
653     !  Construct mid-layer pressures and layer thickness
654     !  -------------------------------------------------
655        do k = 1, km
656           mixr = q(i,j,k) / ( 1.0 - q(i,j,k) ) ! Mixing ratio from specific humidity
657           Tv = T(i,j,k) * ( 1 + 0.61 * mixr )
658           delz(i,j,k)  = Rgas * Tv * log(pe(k+1)/pe(k)) / grav
659        end do
660     
661       end subroutine delz_
662     
663     
664     end
665