File: C:\NOAA\NEMS_11731\src\chem\gocart\src\Components\GOCART_GridComp\CFC_GridComp\CFC_GridCompMod.F90

1     #include "MAPL_Generic.h"
2     !-------------------------------------------------------------------------
3     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
4     !-------------------------------------------------------------------------
5     !BOP
6     !
7     ! !MODULE:  CFC_GridCompMod --- CFC Grid Component Class
8     !
9     ! !INTERFACE:
10     !
11     
12        MODULE  CFC_GridCompMod
13     
14     ! !USES:
15     
16        USE ESMF_Mod
17        USE MAPL_Mod
18        USE Chem_Mod 	     ! Chemistry Base Class
19        USE Chem_StateMod	     ! Chemistry State
20        USE Chem_ConstMod, ONLY: grav
21        USE Chem_UtilMod	     ! I/O
22        USE m_inpak90	     ! Resource file management
23        USE m_die, ONLY: die
24     
25        IMPLICIT NONE
26     
27     ! !PUBLIC TYPES:
28     !
29        PRIVATE
30        PUBLIC  CFC_GridComp       ! The CFC object 
31     
32     !
33     ! !PUBLIC MEMBER FUNCTIONS:
34     !
35     
36        PUBLIC  CFC_GridCompInitialize
37        PUBLIC  CFC_GridCompRun
38        PUBLIC  CFC_GridCompFinalize
39     
40     !
41     ! !DESCRIPTION:
42     !
43     !  This module implements the CFC Grid Component. 
44     !
45     ! !REVISION HISTORY:
46     !
47     !  16Sep2003 da Silva  First crack.
48     !  01Aug2006 da Silva  Extensions for GEOS-5.
49     !   1Jan2008  Nielsen  CFC-12 configuration for ARCTAS.
50     !   8Feb2008  Nielsen  Standard configuration call(s) from AeroChem.
51     !
52     !EOP
53     !-------------------------------------------------------------------------
54     
55       TYPE CFC_GridComp
56     
57         CHARACTER(LEN=255) :: name
58     
59     ! For CFC-12 photolysis
60     ! ---------------------
61         INTEGER :: nlam
62         INTEGER :: nsza
63         INTEGER :: numo3
64         INTEGER :: nx
65         INTEGER :: nxdo
66         INTEGER :: nts
67     
68         REAL(KIND=4), POINTER :: sdat(:,:,:,:)
69         REAL(KIND=4), POINTER :: xtab(:,:)
70         REAL(KIND=4), POINTER :: o3_tab(:,:)
71         REAL(KIND=4), POINTER :: sza_tab(:)
72     
73         REAL, POINTER :: CFCsfcFlux(:,:)	 ! CFC-12 surface flux kg m^-2 s^-1
74         REAL, POINTER :: CFCloss(:,:,:,:)	 ! CFC loss due to photolysis m^-3 s^-1
75     
76       END TYPE CFC_GridComp
77     
78     CONTAINS
79     
80     !-------------------------------------------------------------------------
81     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
82     !-------------------------------------------------------------------------
83     !BOP
84     !
85     ! !IROUTINE:  CFC_GridCompInitialize --- Initialize CFC_GridComp
86     !
87     ! !INTERFACE:
88     !
89     
90        SUBROUTINE CFC_GridCompInitialize( gcCFC, w_c, impChem, expChem, &
91                                           nymd, nhms, cdt, rc )
92     
93     ! !USES:
94     
95       IMPLICIT NONE
96     
97     ! !INPUT PARAMETERS:
98     
99        TYPE(Chem_Bundle), intent(in) :: w_c        ! Chemical tracer fields      
100        INTEGER, INTENT(IN) :: nymd, nhms	       ! time
101        REAL,    INTENT(IN) :: cdt		       ! chemical timestep (secs)
102     
103     
104     ! !OUTPUT PARAMETERS:
105     
106        TYPE(CFC_GridComp), INTENT(INOUT) :: gcCFC   ! Grid Component
107        TYPE(ESMF_State), INTENT(INOUT)  :: impChem  ! Import State
108        TYPE(ESMF_State), INTENT(INOUT)  :: expChem  ! Export State
109        INTEGER, INTENT(OUT) ::  rc                  ! Error return code:
110                                                     !  0 - all is well
111                                                     !  1 - 
112     
113     ! !DESCRIPTION: Initializes the CFC Grid Component. It primarily sets
114     !               the import state for each active constituent package.
115     !
116     ! !REVISION HISTORY:
117     !
118     !  18Sep2003 da Silva  First crack.
119     !  31May2005  Nielsen  Mods for 7 CO bins, 5 region masks
120     !  04Nov2005     Bian  CO tagged to 4 regions 
121     !                      (global, North America, South America, and Africa)
122     !                      for CR-AVE
123     !  12Feb2005  Nielsen  8 regions for INTEX-B 2006
124     !   1Jan2008  Nielsen  CFC-12 configuration for ARCTAS
125     !
126     !EOP
127     !-------------------------------------------------------------------------
128     
129        CHARACTER(LEN=*), PARAMETER :: myname = 'CFC_GridCompInitialize'
130        
131        CHARACTER(LEN=255) :: rcfilen = 'CFC_GridComp.rc'
132     
133        CHARACTER(LEN=255) :: dir4files
134        CHARACTER(LEN=255) :: fnO2Jdat
135        CHARACTER(LEN=255) :: fnO3SZA
136        CHARACTER(LEN=255) :: fnXsectJPL
137        CHARACTER(LEN=255) :: eFileName
138     
139        INTEGER :: ier(128)
140        INTEGER :: i, i1, i2, im, j1, j2, jm, km, nbins
141     
142        gcCFC%name = 'CFC-12 Chemistry for ARCTAS'
143     
144        rc = 0
145     !  Initialize local variables
146     !  --------------------------
147        rc = 0
148        i1 = w_c%grid%i1
149        i2 = w_c%grid%i2
150        im = w_c%grid%im
151        
152        j1 = w_c%grid%j1
153        j2 = w_c%grid%j2
154        jm = w_c%grid%jm
155        
156        km = w_c%grid%km
157        
158        nbins = w_c%reg%n_CFC
159        ier(:)=0
160     
161     !  Initialize photolysis variables
162     !  -------------------------------
163        gcCFC%nlam  =  79
164        gcCFC%nsza  =  20
165        gcCFC%numo3 =  12
166        gcCFC%nx    =  35
167        gcCFC%nxdo  =  33
168        gcCFC%nts   = 200
169     
170     !  Load resource file
171     !  ------------------
172        CALL I90_loadf ( TRIM(rcfilen), ier(1) )
173        IF ( ier(1) .NE. 0 ) THEN
174         rc = 10
175         RETURN
176        END IF
177     
178        CALL I90_label ( 'directory:', ier(30) )
179        CALL I90_Gtoken( dir4files, ier(31) )
180     
181        CALL I90_label ( 'O2Jtable:', ier(32) )
182        CALL I90_Gtoken( fnO2Jdat, ier(33) )
183     
184        CALL I90_label ( 'O3&SZAtables:', ier(36) )
185        CALL I90_Gtoken( fnO3SZA, ier(37) )
186     
187        CALL I90_label ( 'JPLXsections:', ier(40) )
188        CALL I90_Gtoken( fnXsectJPL, ier(41) )
189     
190        CALL I90_label ( 'CFC_emission_filename:', ier(50) )
191        CALL I90_Gtoken( eFileName, ier(51) )
192     
193        IF( ANY( ier(1:128) /= 0 ) ) THEN
194         rc = 12
195         RETURN
196        END IF
197        ier(:)=0
198     
199     ! Allocate space for grid-size dependent arrays
200     ! ---------------------------------------------
201        ALLOCATE(gcCFC%sdat(gcCFC%nsza,gcCFC%numo3,km,gcCFC%nlam), stat=ier( 9) )
202        ALLOCATE(	            gcCFC%xtab(gcCFC%nlam,gcCFC%nts), stat=ier(10) )
203        ALLOCATE(	                gcCFC%o3_tab(gcCFC%numo3,km), stat=ier(11) )
204        ALLOCATE(	                   gcCFC%sza_tab(gcCFC%nsza), stat=ier(12) )
205        ALLOCATE(                   gcCFC%CFCsfcFlux(i1:i2,j1:j2), stat=ier(13) )
206        ALLOCATE(           gcCFC%CFCloss(i1:i2,j1:j2,1:km,nbins), stat=ier(14) )
207     
208        IF( ANY( ier(1:128) /= 0 ) ) THEN
209         rc = 14
210         RETURN
211        END IF
212        ier(:)=0
213     
214     !  Acquire the CFC-12 emissions
215     !  ----------------------------
216        CALL Chem_UtilMPread ( TRIM(eFileName), 'CFC-12_EMISSION', 20080101, &
217        			  120000, i1, i2, 0, im, j1, j2, 0, jm, 0, &
218        			  var2d=gcCFC%CFCsfcFlux, cyclic=.true., &
219        			  grid=w_c%grid_esmf )
220     
221     !  Read the tables
222     !  ---------------
223        CALL rdPhotFiles(km,dir4files,fnO2Jdat,fnO3SZA,fnXsectJPL)
224     
225        RETURN
226        CONTAINS
227        SUBROUTINE rdPhotFiles(km,dir,fnO2Jdat,fnO3SZA,fnXsectJPL)
228     !---------------------------------------------------------------------------
229     !
230     ! Read several external files for the photolysis.
231     !
232     ! Input parameters:
233     !
234     ! km    Grid dimensions
235     ! dir	Directory on which these files reside
236     ! fns   File names
237     !
238     ! Output parameters:
239     !
240     !  None.
241     !
242     ! Restrictions:
243     !
244     !  This runs on each processor
245     !
246     !  This version requires input data at jnp latitudes.
247     !
248     !-----------------------------------------------------------------------
249           IMPLICIT NONE
250     
251           INTEGER, INTENT(IN) :: km
252           
253           CHARACTER(LEN=*), INTENT(IN) :: dir
254           CHARACTER(LEN=*), INTENT(IN) :: fnO2Jdat
255           CHARACTER(LEN=*), INTENT(IN) :: fnO3SZA
256           CHARACTER(LEN=*), INTENT(IN) :: fnXsectJPL
257     
258           REAL(KIND=4), ALLOCATABLE :: dxtab(:,:,:)
259     
260           INTEGER :: i, j, k, l, ierr, iunit, iuchem, kReverse
261           INTEGER :: npr_in, nlam_in, nsza_in, no3_in
262           REAL :: deg2Rad, pi
263           REAL (KIND=4) :: pr_tab(km)
264           REAL (KIND=4) :: rlam(gcCFC%nlam)
265     
266           LOGICAL :: exists,open,found
267     
268           pi = 4.00*ATAN(1.00)
269           deg2Rad = pi/180.00
270     
271     ! Find an available logical unit 
272     ! ------------------------------
273           found=.FALSE.
274           iunit=11
275     
276           DO WHILE (.NOT. found .AND. iunit <= 99)
277            INQUIRE(UNIT=iunit,EXIST=exists,OPENED=open)
278            IF(exists .AND. .NOT. open) THEN
279             found=.TRUE.
280             iuchem=iunit
281            END IF
282            iunit=iunit+1
283           END DO
284     
285           IF(.NOT. found) THEN
286            WRITE(*,FMT="(/,'rdPhotFiles: No available logical units.')")
287            STOP
288           ELSE
289            IF(MAPL_AM_I_ROOT()) THEN
290             WRITE(*,FMT="(' ')")
291             WRITE(*,FMT="(' ','rdPhotFiles: Reading from UNIT ',I3)") iuchem
292     	END IF
293           END IF
294           
295     ! Read in sdat(nsza,numo3,levels,nlam)
296     ! ------------------------------------
297           OPEN(iuchem,FILE=TRIM(dir)//'/'//TRIM(fnO2Jdat),STATUS='old', &
298                 FORM='unformatted',ACTION='read')
299           READ(iuchem) gcCFC%sdat
300           CLOSE(iuchem)
301       
302     ! Read solar zenith angle and O3 references after checking sizes
303     ! --------------------------------------------------------------
304           OPEN(iuchem,FILE = TRIM(dir)//'/'//TRIM(fnO3SZA),STATUS='old', &
305                FORM='unformatted',ACTION='read')
306           READ(iuchem) npr_in,nlam_in,nsza_in,no3_in
307           READ(iuchem) pr_tab
308           READ(iuchem) rlam
309           READ(iuchem) gcCFC%sza_tab
310           READ(iuchem) gcCFC%o3_tab
311           CLOSE(iuchem)
312           IF(( npr_in .NE.         km) .OR. (nlam_in .NE.  gcCFC%nlam) .OR. &
313              (nsza_in .NE. gcCFC%nsza) .OR. ( no3_in .NE. gcCFC%numo3)) THEN
314              PRINT *,'rdPhotFiles: Array sizes of table do not match ', &
315                      ' those expected:',npr_in,km,nlam_in,gcCFC%nlam, &
316                      nsza_in,gcCFC%nsza,no3_in,gcCFC%numo3
317          	 STOP
318           END IF
319      
320     ! Convert sza_tab(nsza) to radians
321     ! --------------------------------
322           DO i=1,gcCFC%nsza
323            gcCFC%sza_tab(i) = gcCFC%sza_tab(i)*deg2Rad
324           END DO
325     
326     ! Reverse sdat(nsza,numo3,levels,nlam) in the vertical to accomodate GEOS-5
327     ! -------------------------------------------------------------------------
328           DO l=1,gcCFC%nlam				       
329            DO j=1,gcCFC%numo3				       
330             DO i=1,gcCFC%nsza				       
331              pr_tab(1:km) = gcCFC%sdat(i,j,1:km,l)		       
332              DO k=1,km					       
333               kReverse = km-k+1				       
334               gcCFC%sdat(i,j,k,l) = pr_tab(kReverse)
335              END DO						       
336             END DO						       
337            END DO						       
338           END DO						       
339     
340     ! Reverse o3_tab(numo3,km) in the vertical to accomodate GEOS-5
341     ! -------------------------------------------------------------
342           DO j=1,gcCFC%numo3					 
343            pr_tab(1:km) = gcCFC%o3_tab(j,1:km)			 
344            DO k=1,km						 
345             kReverse = km-k+1					 
346             gcCFC%o3_tab(j,k) = pr_tab(kReverse)
347            END DO 						 
348           END DO  						 
349     
350           ALLOCATE(dxtab(gcCFC%nlam,gcCFC%nts,gcCFC%nx),STAT=ierr)
351     
352     ! JPL cross sections
353     ! ------------------ 
354           OPEN(iuchem,FILE=TRIM(dir)//'/'//TRIM(fnXsectJPL), &
355                STATUS='old',ACTION='read',FORM='unformatted')
356           READ(iuchem) dxtab
357           CLOSE(iuchem)
358     
359     ! Need only #25
360     ! -------------
361           k=25
362           DO j=1,gcCFC%nts
363            DO i=1,gcCFC%nlam
364             gcCFC%xtab(i,j) = dxtab(i,j,k)
365            END DO
366           END DO
367     
368           DEALLOCATE(dxtab,STAT=ierr)
369     
370           IF(MAPL_AM_I_ROOT()) THEN
371            print *,'rdPhotFiles: Done'
372            print *,' '
373           END IF
374     
375           RETURN
376           END SUBROUTINE rdPhotFiles
377     
378        END SUBROUTINE CFC_GridCompInitialize
379     
380     !-------------------------------------------------------------------------
381     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
382     !-------------------------------------------------------------------------
383     !BOP
384     !
385     ! !IROUTINE:  CFC_GridCompRun --- The CFC Driver 
386     !
387     ! !INTERFACE:
388     !
389     
390        SUBROUTINE CFC_GridCompRun( gcCFC, w_c, impChem, expChem, nymd, nhms, &
391                                    cdt, rc)
392     
393     ! !USES:
394     
395       IMPLICIT NONE
396     
397     ! !INPUT/OUTPUT PARAMETERS:
398     
399        TYPE(CFC_GridComp), INTENT(INOUT) :: gcCFC   ! Grid Component
400        TYPE(Chem_Bundle), INTENT(INOUT) :: w_c	! Chemical tracer fields   
401     
402     ! !INPUT PARAMETERS:
403     
404        TYPE(ESMF_State), INTENT(INOUT) :: impChem   ! Import State
405        INTEGER, INTENT(IN) :: nymd, nhms	        ! time
406        REAL,    INTENT(IN) :: cdt		        ! chemical timestep (secs)
407     
408     ! !OUTPUT PARAMETERS:
409     
410        TYPE(ESMF_State), INTENT(INOUT) :: expChem   ! Export State
411        INTEGER, INTENT(OUT) ::  rc                  ! Error return code:
412                                                     !  0 - all is well
413                                                     !  1 -
414     
415        CHARACTER(LEN=*), PARAMETER :: myname = 'CFC_GridCompRun'
416        CHARACTER(LEN=*), PARAMETER :: Iam = myname
417     
418        INTEGER :: ier(128)
419        INTEGER ::  i1, i2, im, iXj, j1, j2, jm, km, status
420        INTEGER ::  i, indt, j, k, m, n, nbeg, nbins, nend
421        REAL :: o3c, qmin, qmax, r, rg, szan
422     
423     !  Imports
424     !  -------
425        REAL, POINTER, DIMENSION(:,:,:) ::  T  => null()
426        REAL, POINTER, DIMENSION(:,:,:) ::  O3 => null()
427        REAL, POINTER, DIMENSION(:,:)   ::  tropp => null()
428     
429     !  Local Variables
430     !  ---------------
431        REAL, PARAMETER :: badVal=2.00E+05
432        REAL, PARAMETER :: grav=9.80
433        REAL, PARAMETER :: mwtAir=28.97
434        REAL, PARAMETER :: mwtCFC12=120.917
435        REAL, PARAMETER :: Nsuba=6.022E+26
436        REAL, PARAMETER :: rstar=8.3143E+03
437        REAL, PARAMETER :: O3abv80km = 1.10E+15 !m^{-2}
438     
439        REAL, ALLOCATABLE :: emit2vmr(:,:)
440        REAL, ALLOCATABLE :: tropPa(:,:)
441        REAL, ALLOCATABLE :: pPa(:,:,:)
442        REAL, ALLOCATABLE :: nd(:,:,:)
443        REAL, ALLOCATABLE :: O3Col(:,:,:)
444        REAL, ALLOCATABLE :: photoRate(:,:,:)
445        REAL, ALLOCATABLE :: s(:,:,:,:)
446     
447     ! Disable the ACG'd CFC_GetPointer___.h for now. [Maybe fix it soon.]
448     ! -------------------------------------------------------------------
449     #define EXPORT     expChem
450     #define ptrCFCEM   CFC_emis
451     #define ptrCFCLS   CFC_loss
452     #define ptrCFCCL   CFC_column
453     
454     !JEN#include "CFC_GetPointer___.h"
455     
456     !  Bin sizes
457     !  ---------
458        integer, parameter		   :: NBIN_CFCEM = 1 ! CFC Emission
459        integer, parameter		   :: NBIN_CFCLS = 2 ! CFC Loss due to photolysis
460        integer, parameter		   :: NBIN_CFCCL = 2 ! CFC Column
461     
462     !  Bin-indexed Chem Arrays
463     !  -----------------------
464        type(Chem_Array), target	   ::	 CFCEM(NBIN_CFCEM) ! Export: CFC Surface flux
465        type(Chem_Array), pointer	   :: ptrCFCEM(:)	   ! Export: CFC Surface flux
466        type(Chem_Array), target	   ::	 CFCLS(NBIN_CFCLS) ! Export: CFC Loss due to photolysis
467        type(Chem_Array), pointer	   :: ptrCFCLS(:)	   ! Export: CFC Loss due to photolysis
468        type(Chem_Array), target	   ::	 CFCCL(NBIN_CFCCL) ! Export: CFC Column
469        type(Chem_Array), pointer	   :: ptrCFCCL(:)	   ! Export: CFC Column
470     
471     !  Local array referencing the Import/Export states
472     !  ------------------------------------------------
473        type(Chem_Array), target	   ::	 CFC12S ! Export: Stratospheric CFC-12 (CCl2F2)
474        type(Chem_Array), pointer	   :: ptrCFC12S ! Export: Stratospheric CFC-12 (CCl2F2)
475        type(Chem_Array), target	   ::	 CFC12T ! Export: Tropospheric CFC-12 (CCl2F2)
476        type(Chem_Array), pointer	   :: ptrCFC12T ! Export: Tropospheric CFC-12 (CCl2F2)
477     
478     !  Get pointers to data in state
479     !  -----------------------------
480        ptrCFC12S => CFC12S   ! Stratospheric CFC-12 (CCl2F2)
481        call MAPL_GetPointer ( EXPORT, CFC12S%data3d,  'CFC12S', RC=STATUS )
482        VERIFY_(STATUS)
483     
484        ptrCFC12T => CFC12T   ! Tropospheric CFC-12 (CCl2F2)
485        call MAPL_GetPointer ( EXPORT, CFC12T%data3d,  'CFC12T', RC=STATUS )
486        VERIFY_(STATUS)
487     
488        ptrCFCEM => CFCEM	 ! CFC-12 Surface flux
489        call MAPL_GetPointer ( EXPORT, CFCEM(1)%data2d,  'CFC12EM', RC=STATUS )
490        VERIFY_(STATUS)
491     
492        ptrCFCLS => CFCLS	 ! CFC-12 Loss due to photolysis
493        call MAPL_GetPointer ( EXPORT, CFCLS(1)%data3d,  'CFC12SLS', RC=STATUS )
494        VERIFY_(STATUS)
495        call MAPL_GetPointer ( EXPORT, CFCLS(2)%data3d,  'CFC12TLS', RC=STATUS )
496        VERIFY_(STATUS)
497     
498        ptrCFCCL => CFCCL	 ! CFC-12 Column mass density
499        call MAPL_GetPointer ( EXPORT, CFCCL(1)%data2d,  'CFC12SCL', RC=STATUS )
500        VERIFY_(STATUS)
501        call MAPL_GetPointer ( EXPORT, CFCCL(2)%data2d,  'CFC12TCL', RC=STATUS )
502        VERIFY_(STATUS)
503     
504     !  Initialize local variables
505     !  --------------------------
506        rc = 0
507        i1 = w_c%grid%i1
508        i2 = w_c%grid%i2
509        im = w_c%grid%im
510     
511        j1 = w_c%grid%j1
512        j2 = w_c%grid%j2
513        jm = w_c%grid%jm
514     
515        km = w_c%grid%km
516     
517        iXj = ( i2 - i1 + 1 ) * ( j2 - j1 + 1 )
518     
519        nbins = w_c%reg%n_CFC
520        nbeg  = w_c%reg%i_CFC
521        nend  = w_c%reg%j_CFC
522        
523     !  Imports
524     !  -------
525        call MAPL_GetPointer( impChem,     T,     'T', rc=ier(1) ) 
526        call MAPL_GetPointer( impChem,    O3,    'O3', rc=ier(2) ) 
527        call MAPL_GetPointer( impChem, tropp, 'TROPP', rc=ier(3) ) 
528     
529        IF(ANY(ier(:) /= 0 )) THEN
530         rc = 1
531         RETURN
532        END IF
533        ier(:)=0
534     
535     #ifdef DEBUG
536        CALL pmaxmin('    T',     T, qmin, qmax, iXj, km, 1. )
537        CALL pmaxmin('   O3',    O3, qmin, qmax, iXj, km, 1. )
538        CALL pmaxmin('TROPP', tropp, qmin, qmax, iXj,  1, 1. )
539     #endif
540     
541     !  Allocate temporary workspace
542     !  ----------------------------
543        ALLOCATE(    emit2vmr(i1:i2,j1:j2), STAT=ier(1))
544        ALLOCATE(      tropPa(i1:i2,j1:j2), STAT=ier(1))
545        ALLOCATE(      pPa(i1:i2,j1:j2,km), STAT=ier(2))
546        ALLOCATE(       nd(i1:i2,j1:j2,km), STAT=ier(3))
547        ALLOCATE(    O3Col(i1:i2,j1:j2,km), STAT=ier(4))
548        ALLOCATE(photoRate(i1:i2,j1:j2,km), STAT=ier(5))
549     
550        IF(ANY(ier(:) /= 0 )) THEN
551         rc = 10
552         RETURN
553        END IF
554        ier(:)=0
555     
556     !  Fix bad tropopause pressure values if they exist.
557     !  -------------------------------------------------
558        CALL Chem_UtilTroppFixer(i2, j2, tropp, VERBOSE=.TRUE., &
559                                 NEWTROPP=tropPa, RC=STATUS)
560        VERIFY_(STATUS)
561     
562     !  Find the pressure at mid-layer
563     !  ------------------------------
564        pPa(i1:i2,j1:j2,1) = w_c%grid%ptop + 0.50*w_c%delp(i1:i2,j1:j2,1)
565        DO k = 2, km
566         pPa(i1:i2,j1:j2,k) = pPa(i1:i2,j1:j2,k-1) + 0.50* &
567                              (w_c%delp(i1:i2,j1:j2,k-1)+w_c%delp(i1:i2,j1:j2,k))
568        END DO
569     
570     !  Number density
571     !  --------------
572        nd(i1:i2,j1:j2,1:km)= nsuba*pPa(i1:i2,j1:j2,1:km)/(rstar*T(i1:i2,j1:j2,1:km))
573     
574     !  Compute the overlying ozone from mole fraction.  Result: m^{-2}
575     !  ---------------------------------------------------------------
576        r = Nsuba*0.50/(mwtAir*grav)
577        O3col(i1:i2,j1:j2,1) = O3abv80km + O3(i1:i2,j1:j2,1)*w_c%delp(i1:i2,j1:j2,1)*r
578        DO k=2,km
579         O3col(i1:i2,j1:j2,k) = O3col(i1:i2,j1:j2,k-1) + &
580                                  (O3(i1:i2,j1:j2,k-1) * w_c%delp(i1:i2,j1:j2,k-1) + &
581                                   O3(i1:i2,j1:j2,  k) * w_c%delp(i1:i2,j1:j2,  k))*r
582        END DO
583     
584     !  Enable the conversion from emission [kg CFC m^{-2} s^{-1}] 
585     !  to an incremental change in the mixing ratio [s^{-1}].
586     !  ----------------------------------------------------------
587        emit2vmr(i1:i2,j1:j2) = mwtAir*grav/(mwtCFC12*w_c%delp(i1:i2,j1:j2,km))
588     
589     !  Increment mixing ratio in surface layer of tropospheric CFC-12
590     !  --------------------------------------------------------------
591        w_c%qa(nbeg+1)%data3d(i1:i2,j1:j2,km) = w_c%qa(nbeg+1)%data3d(i1:i2,j1:j2,km)+cdt* &
592                                                gcCFC%CFCsfcFlux(i1:i2,j1:j2)*emit2vmr(i1:i2,j1:j2)
593     
594     !  When tropospheric CFC-12 migrates to the stratosphere, reassign it
595     !  ------------------------------------------------------------------
596        DO k = 1, km
597         WHERE(pPa(i1:i2,j1:j2,k) < tropPa(i1:i2,j1:j2) .AND. &
598               w_c%qa(nbeg+1)%data3d(i1:i2,j1:j2,k) > 0.00 )
599          w_c%qa(nbeg)%data3d(i1:i2,j1:j2,k) = w_c%qa(nbeg)%data3d(i1:i2,j1:j2,k) + &
600                                               w_c%qa(nbeg+1)%data3d(i1:i2,j1:j2,k)
601          w_c%qa(nbeg+1)%data3d(i1:i2,j1:j2,k) = 0.00
602         END WHERE
603        END DO
604     
605     !  Convert CFC-12 to number density
606     !  --------------------------------
607        DO n=nbeg,nend
608         w_c%qa(n)%data3d(i1:i2,j1:j2,1:km) = w_c%qa(n)%data3d(i1:i2,j1:j2,1:km)* &
609                                              nd(i1:i2,j1:j2,1:km)
610        END DO
611     
612        ALLOCATE(s(gcCFC%nlam,i1:i2,j1:j2,1:km), STAT=ier(1))
613     
614     !  Photolysis:  Loop over horizontal domain
615     !  ----------------------------------------
616        DO j=j1,j2
617         DO i=i1,i2
618     
619     !  Solar zenith angle (radians).  w_c%cosz has no negative values,
620     !  which are required for correct interpolation in the S-dat tables.
621     !  -----------------------------------------------------------------
622          IF(w_c%cosz(i,j) <= 1.00E-06) THEN
623           szan = ACOS(-0.50)
624          ELSE
625           szan = ACOS(w_c%cosz(i,j))
626          END IF
627     
628          DO k=1,km
629           o3c = O3Col(i,j,k)*1.00E-04 !to cm^{-2}
630     
631     !  Interpolate radiative flux function values.
632     !  Call getS even when sun is below the horizon.
633     !  ---------------------------------------------
634           CALL getS(k,km,szan,o3c,s(:,i,j,k))
635           indt = T(i,j,k)-148.5
636           indt = MAX(1,indt)
637           indt = MIN(indt,200)
638     
639     !  Rate constant is sum over wavelengths
640     !  -------------------------------------
641           photoRate(i,j,k) = SUM(s(1:gcCFC%nlam,i,j,k)*gcCFC%xtab(1:gcCFC%nlam,indt))
642     
643          END DO ! Layer
644         END DO  ! Longitude
645        END DO   ! Latitude
646     
647        DEALLOCATE(s, STAT=ier(2))
648        m = 0
649     
650     !  Apply photolysis
651     !  ----------------
652        DO n=nbeg,nend
653         m = m+1
654         w_c%qa(n)%data3d(i1:i2,j1:j2,1:km) = w_c%qa(n)%data3d(i1:i2,j1:j2,1:km) - cdt * &
655                                              w_c%qa(n)%data3d(i1:i2,j1:j2,1:km) * &
656                                              photoRate(i1:i2,j1:j2,1:km)
657         gcCFC%CFCloss(i1:i2,j1:j2,1:km,m) = w_c%qa(n)%data3d(i1:i2,j1:j2,1:km) * &
658                                             photoRate(i1:i2,j1:j2,1:km)
659        END DO
660      
661     !  Return CFC-12 to mole fraction
662     !  ------------------------------
663        DO n=nbeg,nend
664         w_c%qa(n)%data3d(i1:i2,j1:j2,1:km) = w_c%qa(n)%data3d(i1:i2,j1:j2,1:km)/ &
665                                              nd(i1:i2,j1:j2,1:km)
666        END DO
667     
668     !  Fill the export states.
669     
670     !  CFC-12 Surface emission in kg m^{-2} s^{-1}
671     !  -------------------------------------------
672        IF(ASSOCIATED(CFC_emis(1)%data2d)) &
673         CFC_emis(1)%data2d(i1:i2,j1:j2) = gcCFC%CFCsfcFlux(i1:i2,j1:j2)
674     
675     !  Loss due to photolysis: Currently m^{-3} s^(-1), and positive for loss.
676     !  -----------------------------------------------------------------------
677        DO n = 1, nbins
678         IF(ASSOCIATED(CFC_loss(n)%data3d)) &
679          CFC_loss(n)%data3d(i1:i2,j1:j2,1:km) = gcCFC%CFCloss(i1:i2,j1:j2,1:km,n)
680        END DO
681     
682     !  Column burden in kg m(^-2)
683     !  --------------------------
684        DO n = 1, nbins
685          IF(ASSOCIATED(CFC_column(n)%data2d)) THEN
686           CFC_column(n)%data2d(i1:i2,j1:j2) = 0.
687           DO k = 1, km
688            CFC_column(n)%data2d(i1:i2,j1:j2) &
689             =   CFC_column(n)%data2d(i1:i2,j1:j2) &
690               +   w_c%qa(nbeg+n-1)%data3d(i1:i2,j1:j2,k)*mwtCFC12/mwtAir &
691                 * w_c%delp(i1:i2,j1:j2,k)/grav
692          END DO
693         END IF
694        END DO
695     
696     !  Clean up
697     !  --------
698        DEALLOCATE( emit2vmr, STAT=ier(1))
699        DEALLOCATE(      pPa, STAT=ier(2))
700        DEALLOCATE(       nd, STAT=ier(3))
701        DEALLOCATE(    O3Col, STAT=ier(4))
702        DEALLOCATE(photoRate, STAT=ier(5))
703     
704        IF(ANY(ier(:) /= 0 )) THEN
705         rc = 99
706         RETURN
707        END IF
708        ier(:)=0
709     
710        RETURN
711     
712        CONTAINS
713     
714        SUBROUTINE getS(ik,levels,sza,o3column,s)
715     ! --------------------------------------------------------------------------
716     ! NAME:
717     !   interp_s
718     ! PURPOSE:
719     !   Interpolate s values for each wavelength in table to specified O3
720     !   col and zenith angles
721     ! CATEGORY:
722     ! CALLING SEQUENCE:
723     !   Call interp_s(nlam,sza,o3column,s)
724     ! INPUTS:
725     !   nlam     --  number of wavelength intervals used
726     !   sza      --  zenith angle
727     !   o3column --  overhead o3 column value 
728     ! OPTIONAL INPUT PARAMETERS:
729     ! OUTPUTS:
730     !   s -- array of s values (nlam) for each wavelength 
731     !	 at model p-level interpolated to o3column and sza values
732     ! INTERNAL VARIABLES
733     !   sza_tab -- values of sza corresponding to sdat table values
734     !   o3_tab  -- array of overhead O3 values at each p-level (numo3s,np_ctm)
735     !		used to index sdat
736     !   sdat    -- input array of values of radiative source function 
737     !	       (nzens,numo3,np_ctm,nlam) gridded to ctm p layers
738     ! COMMON BLOCKS:
739     ! SIDE EFFECTS:
740     ! PROCEDURE:
741     !   bi-linear interpolation, for sza>94 s=0, for o3 out of range use min/max
742     ! RESTRICTIONS:
743     ! REQUIRED ROUTINES:
744     ! MODIFICATION HISTORY: 
745     !   Created 930825 - SR Kawa
746     !   Modified 960710 for 28 levels and to handle J(O2) separately
747     !   1Jan2008  Nielsen  CFC-12 configuration for ARCTAS.
748     ! --------------------------------------------------------------------------
749     
750        IMPLICIT NONE
751     
752        INTEGER, INTENT(IN) :: ik,levels
753        REAL, INTENT(IN) :: sza,o3column 
754        REAL, INTENT(OUT) :: s(gcCFC%nlam)
755     
756        INTEGER :: ijj,ikk,ikkm,il,is
757        REAL :: omt,omu,t,u
758            
759     ! For each input solar zenith angle, find the first element of
760     ! tabled sza_tab values that is greater than it.  Use this
761     ! table element and previous table element to determined
762     ! interpolated value.
763     ! ------------------------------------------------------------
764        DO is=1,gcCFC%nsza
765           ijj = is 
766           if(gcCFC%sza_tab(is) > sza) EXIT
767        END DO
768           
769     ! Location is dark, set s/jo2=0        
770     ! -----------------------------
771        IF(sza > gcCFC%sza_tab(gcCFC%nsza)) THEN
772           s(1:gcCFC%nlam) = 0.
773        ELSE
774           t = (sza-gcCFC%sza_tab(ijj-1))/(gcCFC%sza_tab(ijj)-gcCFC%sza_tab(ijj-1))
775           omt = 1.-t
776              
777     ! For each input overhead o3 column find the first element
778     ! of tabled o3_tab values that is > than it.  Use this
779     ! table element and previous table element to determine
780     ! interpolated value
781     ! --------------------------------------------------------
782           DO is=1,gcCFC%numo3
783       	 ikk = is 
784       	 IF (gcCFC%o3_tab(is,ik) > o3column) EXIT
785           END DO 
786     
787           ikkm = ikk-1
788     
789           IF(ikk > 1 .AND. o3column <= gcCFC%o3_tab(gcCFC%numo3,ik)) THEN 
790       	 u = (o3column-gcCFC%o3_tab(ikkm,ik))/ &
791       	     (gcCFC%o3_tab(ikk,ik)-gcCFC%o3_tab(ikkm,ik))
792       	 omu = 1.-u
793     
794     ! Bilinear interpolation at ik for each wavelength
795     ! ------------------------------------------------
796         	 DO il=1,gcCFC%nlam	    
797         	    s(il) = omt*omu*gcCFC%sdat(ijj-1,ikkm,ik,il) &
798         		 +t*omu*gcCFC%sdat(ijj,ikkm,ik,il) &
799         		 +t*u*gcCFC%sdat(ijj,ikk,ik,il) &
800         		 +omt*u*gcCFC%sdat(ijj-1,ikk,ik,il)
801         	 END DO
802         
803     ! Extrapolate before table
804     ! ------------------------
805           ELSE IF (ikk == 1) THEN
806        	 DO il=1,gcCFC%nlam
807        	    s(il) = omt*gcCFC%sdat(ijj-1,1,ik,il)+t*gcCFC%sdat(ijj,1,ik,il)
808        	 END DO
809     
810     ! Extrapolate past table
811     ! ----------------------
812           ELSE 
813       	 DO il=1,gcCFC%nlam
814       	    s(il) = omt*gcCFC%sdat(ijj-1,gcCFC%numo3,ik,il)+ &
815       		    t*gcCFC%sdat(ijj,gcCFC%numo3,ik,il)
816       	 END DO 
817           END IF  
818        END IF
819        
820        RETURN
821        END SUBROUTINE getS
822     
823        END SUBROUTINE CFC_GridCompRun
824     
825     !-------------------------------------------------------------------------
826     !     NASA/GSFC, Global Modeling and Assimilation Office, Code 900.3     !
827     !-------------------------------------------------------------------------
828     !BOP
829     !
830     ! !IROUTINE:  CFC_GridCompFinalize
831     !
832     ! !INTERFACE:
833     !
834     
835        SUBROUTINE CFC_GridCompFinalize( gcCFC, w_c, impChem, expChem, &
836                                         nymd, nhms, cdt, rc )
837     
838     ! !USES:
839     
840       IMPLICIT NONE
841     
842     ! !INPUT/OUTPUT PARAMETERS:
843     
844        TYPE(CFC_GridComp), INTENT(INOUT) :: gcCFC ! Grid Component
845     
846     ! !INPUT PARAMETERS:
847     
848        TYPE(Chem_Bundle), INTENT(IN)  :: w_c      ! Chemical tracer fields   
849        INTEGER, INTENT(IN) :: nymd, nhms	      ! time
850        REAL,    INTENT(IN) :: cdt		      ! chemical timestep (secs)
851     
852     
853     ! !OUTPUT PARAMETERS:
854     
855        TYPE(ESMF_State), INTENT(INOUT) :: impChem	! Import State
856        TYPE(ESMF_State), INTENT(INOUT) :: expChem	! Import State
857        INTEGER, INTENT(OUT) ::  rc                  ! Error return code:
858                                                     !  0 - all is well
859                                                     !  1 -
860      
861     ! !DESCRIPTION: This routine finalizes this Grid Component.
862     !
863     ! !REVISION HISTORY:
864     !
865     !  18Sep2003 da Silva  First crack.
866     !
867     !EOP
868     !-------------------------------------------------------------------------
869     
870        CHARACTER(LEN=*), PARAMETER :: myname = 'CFC_GridCompFinalize'
871        INTEGER :: ios
872     
873        rc = 0
874     
875        DEALLOCATE(gcCFC%sdat, gcCFC%xtab, gcCFC%o3_tab, gcCFC%sza_tab, &
876                   gcCFC%CFCloss, gcCFC%CFCsfcFlux, STAT=ios )
877     
878        IF( ios /= 0 ) THEN
879         rc = 1
880         IF(MAPL_AM_I_ROOT()) PRINT *,myname,': DEALLOCATE return code is ',ios
881        END IF
882     
883        RETURN
884        END SUBROUTINE CFC_GridCompFinalize
885     
886      END MODULE CFC_GridCompMod
887     
888