!************************************************************************ ! This computer software was prepared by Battelle Memorial Institute, ! hereinafter the Contractor, under Contract No. DE-AC05-76RL0 1830 with ! the Department of Energy (DOE). NEITHER THE GOVERNMENT NOR THE ! CONTRACTOR MAKES ANY WARRANTY, EXPRESS OR IMPLIED, OR ASSUMES ANY ! LIABILITY FOR THE USE OF THIS SOFTWARE. ! ! Photolysis Option: Fast-J ! * Primary investigators: Elaine G. Chapman and James C. Barnard ! * Co-investigators: Jerome D. Fast, William I. Gustafson Jr. ! Last update: February 2009 ! ! Contact: ! Jerome D. Fast, PhD ! Staff Scientist ! Pacific Northwest National Laboratory ! P.O. Box 999, MSIN K9-30 ! Richland, WA, 99352 ! Phone: (509) 372-6116 ! Email: Jerome.Fast@pnl.gov ! ! The original Fast-J code was provided by Oliver Wild (Univ. of Calif. ! Irvine); however, substantial modifications were necessary to make it ! compatible with WRF-Chem and to include the effect of prognostic ! aerosols on photolysis rates. ! ! Please report any bugs or problems to Jerome Fast, the WRF-chem ! implmentation team leader for PNNL. ! ! References: ! * Wild, O., X. Zhu, M.J. Prather, (2000), Accurate simulation of in- ! and below cloud photolysis in tropospheric chemical models, J. Atmos. ! Chem., 37, 245-282. ! * Barnard, J.C., E.G. Chapman, J.D. Fast, J.R. Schmelzer, J.R. ! Schulsser, and R.E. Shetter (2004), An evaluation of the FAST-J ! photolysis model for predicting nitrogen dioxide photolysis rates ! under clear and cloudy sky conditions, Atmos. Environ., 38, ! 3393-3403. ! * Fast, J.D., W.I. Gustafson Jr., R.C. Easter, R.A. Zaveri, J.C. ! Barnard, E.G. Chapman, G.A. Grell, and S.E. Peckham (2005), Evolution ! of ozone, particulates, and aerosol direct radiative forcing in the ! vicinity of Houston using a fully-coupled meteorology-chemistry- ! aerosol model, J. Geophys. Res., 111, D21305, ! doi:10.1029/2005JD006721. ! * Barnard, J.C., J.D. Fast, G. Paredes-Miranda, W.P. Arnott, and ! A. Laskin (2010), Technical note: evaluation of the WRF-Chem ! "aerosol chemical to aerosol optical properties" module using data ! from the MILAGRO campaign, Atmos. Chem. Phys., 10, 7325-7340, ! doi:10.5194/acp-10-7325-2010. ! ! Contact Jerome Fast for updates on the status of manuscripts under ! review. ! ! Additional information: ! * www.pnl.gov/atmospheric/research/wrf-chem ! ! Support: ! Funding for adapting Fast-J was provided by the U.S. Department of ! Energy under the auspices of Atmospheric Science Program of the Office ! of Biological and Environmental Research the PNNL Laboratory Research ! and Directed Research and Development program. !************************************************************************ !WRF:MODEL_LAYER:CHEMISTRY ! module module_phot_fastj integer, parameter :: lunerr = -1 contains !*********************************************************************** subroutine fastj_driver(id,curr_secs,dtstep,config_flags, & gmt,julday,t_phy,moist,p8w,p_phy, & chem,rho_phy,dz8w,xlat,xlong,z_at_w, & ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, & ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, & ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, & ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,& ph_n2o5, & tauaer1,tauaer2,tauaer3,tauaer4, & gaer1,gaer2,gaer3,gaer4, & waer1,waer2,waer3,waer4, & bscoef1,bscoef2,bscoef3,bscoef4, & l2aer,l3aer,l4aer,l5aer,l6aer,l7aer, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte ) USE module_configure USE module_state_description USE module_data_mosaic_therm, only: nbin_a, nbin_a_maxd ! USE module_data_mosaic_asect USE module_data_mosaic_other, only: kmaxd, nsubareas USE module_fastj_mie ! USE module_mosaic_therm, only: aerosol_optical_properties ! USE module_mosaic_driver, only: mapaer_tofrom_host USE module_fastj_data, only: nb, nc implicit none !jdf ! integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer,save :: lpar !Number of levels in CTM integer,save :: jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1),save :: sizeaer,extaer,waer,gaer,tauaer,bscoef real, dimension (nspint, kmaxd+1),save :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / integer nphoto_fastj parameter (nphoto_fastj = 14) integer & lfastj_no2, lfastj_o3a, lfastj_o3b, lfastj_h2o2, & lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x, & lfastj_no3l, lfastj_hono, lfastj_n2o5, lfastj_hno3, & lfastj_hno4 parameter( lfastj_no2 = 1 ) parameter( lfastj_o3a = 2 ) parameter( lfastj_o3b = 3 ) parameter( lfastj_h2o2 = 4 ) parameter( lfastj_hchoa = 5 ) parameter( lfastj_hchob = 6 ) parameter( lfastj_ch3ooh= 7 ) parameter( lfastj_no3x = 8 ) parameter( lfastj_no3l = 9 ) parameter( lfastj_hono = 10 ) parameter( lfastj_n2o5 = 11 ) parameter( lfastj_hno3 = 12 ) parameter( lfastj_hno4 = 13 ) !jdf INTEGER, INTENT(IN ) :: id,julday, & ids,ide, jds,jde, kds,kde, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte REAL(KIND=8), INTENT(IN ) :: curr_secs REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_moist ), & INTENT(IN ) :: moist REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(INOUT ) :: & ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, & ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, & ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, & ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, & ph_n2o5 REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(IN ) :: & tauaer1,tauaer2,tauaer3,tauaer4, & gaer1,gaer2,gaer3,gaer4, & waer1,waer2,waer3,waer4, & bscoef1,bscoef2,bscoef3,bscoef4 REAL, DIMENSION( ims:ime, kms:kme, jms:jme, 1:4 ), & INTENT(IN ) :: & l2aer,l3aer,l4aer,l5aer,l6aer,l7aer REAL, DIMENSION( ims:ime, kms:kme, jms:jme, num_chem ), & INTENT(INOUT ) :: chem REAL, DIMENSION( ims:ime , kms:kme , jms:jme ) , & INTENT(IN ) :: & t_phy, & p_phy, & dz8w, & p8w, & rho_phy, & z_at_w REAL, DIMENSION( ims:ime , jms:jme ) , & INTENT(IN ) :: xlat, & xlong REAL, INTENT(IN ) :: & dtstep,gmt TYPE(grid_config_rec_type), INTENT(IN ) :: config_flags !local variables integer iclm, jclm real number_bin(nbin_a_maxd,kmaxd) !These have to be the full kmaxd to real radius_wet(nbin_a_maxd,kmaxd) !match arrays in MOSAIC routines. complex refindx(nbin_a_maxd,kmaxd) integer(kind=8) :: ixhour integer i,j, k, nsub real(kind=8) :: xtime, xhour real xmin, gmtp, tmidh real sla, slo real psfc real cos_sza real, dimension(kts:kte) :: temp, ozone, dz real, dimension(0:kte) :: pbnd real, dimension(kts:kte) :: cloudmr, airdensity, relhum real, dimension(kts:kte+1) :: zatw real valuej(kte,nphoto_fastj) logical processingAerosols ! set "pegasus" grid size variables ! itot = ite ! jtot = jte lpar = kte jpnl = kte nb = lpar + 1 !for module module_fastj_data nc = 2*nb !ditto, and don't confuse this with nc in module_fastj_mie nsubareas = 1 if ((kte > kmaxd) .or. (lpar <= 0)) then write( wrf_err_message, '(a,4i5)' ) & '*** subr fastj_driver -- ' // & 'lpar, kmaxd, kts, kte', lpar, kmaxd, kts, kte call wrf_message( trim(wrf_err_message) ) wrf_err_message = '*** subr fastj_driver -- ' // & 'kte>kmaxd OR lpar<=0' call wrf_error_fatal( wrf_err_message ) end if ! Determine if aerosol data is provided in the chem array. Currently, ! only MOSAIC will work. The Mie routine does not know how to handle ! SORGAM aerosols. select case (config_flags%chem_opt) case ( CBMZ_MOSAIC_4BIN, CBMZ_MOSAIC_8BIN, CBMZ_MOSAIC_KPP, & CBMZ_MOSAIC_4BIN_AQ, CBMZ_MOSAIC_8BIN_AQ, & CBMZ_MOSAIC_DMS_4BIN, CBMZ_MOSAIC_DMS_8BIN, & CBMZ_MOSAIC_DMS_4BIN_AQ, CBMZ_MOSAIC_DMS_8BIN_AQ, & MOZART_MOSAIC_4BIN_KPP, MOZART_MOSAIC_4BIN_AQ_KPP, & CRI_MOSAIC_8BIN_AQ_KPP, CRI_MOSAIC_4BIN_AQ_KPP, SAPRC99_MOSAIC_8BIN_VBS2_AQ_KPP, &!BSINGH(12/05/2013): Added for SAPRC 8 bin vbs SAPRC99_MOSAIC_8BIN_VBS2_KPP )!BSINGH(04/07/2014): Added for SAPRC 8 bin vbs non-aq processingAerosols = .true. case default processingAerosols = .false. end select ! Set nbin_a = ntot_amode and check nbin_a <= nbin_a_maxd. ! This duplicates the same assignment and check in module_mosaic_therm.F, ! but photolysis is called before aeorosols so this must be set too. ! ! rce 2004-dec-07 -- nbin_a is initialized elsewhere !!$ nbin_a = ntot_amode !!$ if ((nbin_a .gt. nbin_a_maxd) .or. (nbin_a .le. 0)) then !!$ write( wrf_err_message, '(a,2(1x,i4))' ) & !!$ '*** subr fastj_driver -- nbin_a, nbin_a_maxd =', & !!$ nbin_a, nbin_a_maxd !!$ call wrf_message( wrf_err_message ) !!$ call wrf_error_fatal( & !!$ '*** subr fastj_driver -- BAD VALUE for nbin_a' ) !!$ end if ! determine current time of day in Greenwich Mean Time at middle ! of current time step, tmidh. do this by computing the number of minutes ! from beginning of run to middle of current time step xtime = curr_secs/60._8 + real(dtstep,8)/120._8 ixhour = int(gmt + 0.01,8) + int(xtime/60._8,8) xhour = real(ixhour,8) !current hour xmin = 60.*gmt + real(xtime-xhour*60_8,8) gmtp = mod(xhour,24._8) tmidh = gmtp + xmin/60. ! execute for each i,j column and each nsub subarea do nsub = 1, nsubareas do jclm = jts, jte do iclm = its, ite do k = kts, lpar dz(k) = dz8w(iclm, k, jclm) ! cell depth (m) end do if( processingAerosols ) then do k = kts, lpar l2(1,k)=l2aer(iclm,k,jclm,1) l2(2,k)=l2aer(iclm,k,jclm,2) l2(3,k)=l2aer(iclm,k,jclm,3) l2(4,k)=l2aer(iclm,k,jclm,4) l3(1,k)=l3aer(iclm,k,jclm,1) l3(2,k)=l3aer(iclm,k,jclm,2) l3(3,k)=l3aer(iclm,k,jclm,3) l3(4,k)=l3aer(iclm,k,jclm,4) l4(1,k)=l4aer(iclm,k,jclm,1) l4(2,k)=l4aer(iclm,k,jclm,2) l4(3,k)=l4aer(iclm,k,jclm,3) l4(4,k)=l4aer(iclm,k,jclm,4) l5(1,k)=l5aer(iclm,k,jclm,1) l5(2,k)=l5aer(iclm,k,jclm,2) l5(3,k)=l5aer(iclm,k,jclm,3) l5(4,k)=l5aer(iclm,k,jclm,4) l6(1,k)=l6aer(iclm,k,jclm,1) l6(2,k)=l6aer(iclm,k,jclm,2) l6(3,k)=l6aer(iclm,k,jclm,3) l6(4,k)=l6aer(iclm,k,jclm,4) l7(1,k)=l7aer(iclm,k,jclm,1) l7(2,k)=l7aer(iclm,k,jclm,2) l7(3,k)=l7aer(iclm,k,jclm,3) l7(4,k)=l7aer(iclm,k,jclm,4) enddo ! take chem data and extract 1 column to create rsub(l,k,m) array ! call mapaer_tofrom_host( 0, & ! ims,ime, jms,jme, kms,kme, & ! its,ite, jts,jte, kts,kte, & ! iclm, jclm, kts,lpar, & ! num_moist, num_chem, moist, chem, & ! t_phy, p_phy, rho_phy ) ! generate aerosol optical properties for cells in this column ! subroutine is located in file module_mosaic_therm ! call aerosol_optical_properties(iclm, jclm, lpar, refindx, & ! radius_wet, number_bin) ! execute mie code , located in file module_fastj_mie ! CALL wrf_debug(250,'fastj_driver: calling mieaer') ! call mieaer( & ! id, iclm, jclm, nbin_a, & ! number_bin, radius_wet, refindx, & ! dz, curr_secs, lpar, & ! sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7,bscoef) end if ! set lat, long sla = xlat(iclm,jclm) slo = xlong(iclm,jclm) ! set column pressures, temperature, and ozone psfc = p8w(iclm,1,jclm) * 10. ! convert pascals to dynes/cm2 do k = kts, lpar pbnd(k) = p8w(iclm,k+1,jclm) *10. ! convert pascals to dynes/cm2 temp(k) = t_phy(iclm,k,jclm) ozone(k) = chem(iclm,k,jclm,p_o3) / 1.0e6 ! ppm->mol/mol air cloudmr(k) = moist(iclm,k,jclm,p_qc)/0.622 airdensity(k) = rho_phy(iclm,k,jclm)/28.966e3 relhum(k) = MIN( .95, moist(iclm,k,jclm,p_qv) / & (3.80*exp(17.27*(t_phy(iclm,k,jclm)-273.)/ & (t_phy(iclm,k,jclm)-36.))/(.01*p_phy(iclm,k,jclm)))) relhum(k) = MAX(.001,relhum(k)) zatw(k)=z_at_w(iclm,k,jclm) end do zatw(lpar+1)=z_at_w(iclm,lpar+1,jclm) ! call interface_fastj CALL wrf_debug(250,'fastj_driver: calling interface_fastj') call interface_fastj(tmidh,sla,slo,julday, & pbnd, psfc, temp, ozone, & dz, cloudmr, airdensity, relhum, zatw, & iclm, jclm, lpar, jpnl, & curr_secs, valuej, cos_sza, processingAerosols, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) ! put column photolysis rates (valuej) into wrf photolysis (i,k,j) arrays CALL wrf_debug(250,'fastj_driver: calling mapJrates_tofrom_host') call mapJrates_tofrom_host( 0, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & iclm, jclm, kts,lpar, & valuej, & ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, & ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, & ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, & ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, & ph_n2o5 ) ! put the aerosol optical properties into the wrf arrays (this is hard- ! coded to 4 spectral bins, nspint) ! do k=kts,kte ! tauaer1(iclm,k,jclm) = tauaer(1,k) ! tauaer2(iclm,k,jclm) = tauaer(2,k) ! tauaer3(iclm,k,jclm) = tauaer(3,k) ! tauaer4(iclm,k,jclm) = tauaer(4,k) ! gaer1(iclm,k,jclm) = gaer(1,k) ! gaer2(iclm,k,jclm) = gaer(2,k) ! gaer3(iclm,k,jclm) = gaer(3,k) ! gaer4(iclm,k,jclm) = gaer(4,k) ! waer1(iclm,k,jclm) = waer(1,k) ! waer2(iclm,k,jclm) = waer(2,k) ! waer3(iclm,k,jclm) = waer(3,k) ! waer4(iclm,k,jclm) = waer(4,k) ! bscoef1(iclm,k,jclm) = bscoef(1,k) ! bscoef2(iclm,k,jclm) = bscoef(2,k) ! bscoef3(iclm,k,jclm) = bscoef(3,k) ! bscoef4(iclm,k,jclm) = bscoef(4,k) ! end do end do end do end do return end subroutine fastj_driver !----------------------------------------------------------------------- subroutine mapJrates_tofrom_host( iflag, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & iclm, jclm, ktmaps,ktmape, & valuej, & ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, & ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, & ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, & ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob, & ph_n2o5 ) USE module_data_cbmz implicit none !jdf integer nphoto_fastj parameter (nphoto_fastj = 14) integer & lfastj_no2, lfastj_o3a, lfastj_o3b, lfastj_h2o2, & lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x, & lfastj_no3l, lfastj_hono, lfastj_n2o5, lfastj_hno3, & lfastj_hno4 parameter( lfastj_no2 = 1 ) parameter( lfastj_o3a = 2 ) parameter( lfastj_o3b = 3 ) parameter( lfastj_h2o2 = 4 ) parameter( lfastj_hchoa = 5 ) parameter( lfastj_hchob = 6 ) parameter( lfastj_ch3ooh= 7 ) parameter( lfastj_no3x = 8 ) parameter( lfastj_no3l = 9 ) parameter( lfastj_hono = 10 ) parameter( lfastj_n2o5 = 11 ) parameter( lfastj_hno3 = 12 ) parameter( lfastj_hno4 = 13 ) !jdf INTEGER, INTENT(IN ) :: iflag, & ims,ime, jms,jme, kms,kme, & its,ite, jts,jte, kts,kte, & iclm, jclm, ktmaps, ktmape REAL, DIMENSION( ims:ime, kms:kme, jms:jme ), & INTENT(INOUT ) :: & ph_o2,ph_o31d,ph_o33p,ph_no2,ph_no3o2,ph_no3o,ph_hno2, & ph_hno3,ph_hno4,ph_h2o2,ph_ch2or,ph_ch2om,ph_ch3cho, & ph_ch3coch3,ph_ch3coc2h5,ph_hcocho,ph_ch3cocho, & ph_hcochest,ph_ch3o2h,ph_ch3coo2h,ph_ch3ono2,ph_hcochob,& ph_n2o5 REAL, DIMENSION( kte,nphoto_fastj ), INTENT(INOUT) :: valuej ! local variables real ft integer kt ft = 60. if (iflag .gt. 0) go to 2000 ! flag is <=0, put pegasus column J rates (in 1/sec) into WRF arrays (in 1/min) do kt = ktmaps, ktmape ph_no2(iclm,kt,jclm) = valuej(kt,lfastj_no2) * ft ph_no3o(iclm,kt,jclm) = valuej(kt,lfastj_no3x) * ft ph_no3o2(iclm,kt,jclm) = valuej(kt,lfastj_no3l) * ft ph_o33p(iclm,kt,jclm) = valuej(kt,lfastj_o3a) * ft ph_o31d(iclm,kt,jclm) = valuej(kt,lfastj_o3b) * ft ph_hno2(iclm,kt,jclm) = valuej(kt,lfastj_hono) * ft ph_hno3(iclm,kt,jclm) = valuej(kt,lfastj_hno3) * ft ph_hno4(iclm,kt,jclm) = valuej(kt,lfastj_hno4) * ft ph_h2o2(iclm,kt,jclm) = valuej(kt,lfastj_h2o2) * ft ph_ch3o2h(iclm,kt,jclm) = valuej(kt,lfastj_ch3ooh) * ft ph_ch2or(iclm,kt,jclm) = valuej(kt,lfastj_hchoa) * ft ph_ch2om(iclm,kt,jclm) = valuej(kt,lfastj_hchob) * ft ph_n2o5(iclm,kt,jclm) = valuej(kt,lfastj_n2o5) * ft ph_o2(iclm,kt,jclm) = 0.0 ph_ch3cho(iclm,kt,jclm) = 0.0 ph_ch3coch3(iclm,kt,jclm) = 0.0 ph_ch3coc2h5(iclm,kt,jclm) = 0.0 ph_hcocho(iclm,kt,jclm) = 0.0 ph_ch3cocho(iclm,kt,jclm) = 0.0 ph_hcochest(iclm,kt,jclm) = 0.0 ph_ch3coo2h(iclm,kt,jclm) = 0.0 ph_ch3ono2(iclm,kt,jclm) = 0.0 ph_hcochob(iclm,kt,jclm) = 0.0 end do return !finished peg-> wrf mapping 2000 continue ! iflag > 0 ; put wrf ph_xxx Jrates (1/min) into pegasus column valuej (1/sec) do kt = ktmaps, ktmape valuej(kt,lfastj_no2) = ph_no2(iclm,kt,jclm) / ft valuej(kt,lfastj_no3x) = ph_no3o(iclm,kt,jclm) / ft valuej(kt,lfastj_no3l) = ph_no3o2(iclm,kt,jclm)/ ft valuej(kt,lfastj_o3a) = ph_o33p(iclm,kt,jclm) / ft valuej(kt,lfastj_o3b) = ph_o31d(iclm,kt,jclm) / ft valuej(kt,lfastj_hono) = ph_hno2(iclm,kt,jclm) / ft valuej(kt,lfastj_hno3) = ph_hno3(iclm,kt,jclm) / ft valuej(kt,lfastj_hno4) = ph_hno4(iclm,kt,jclm) / ft valuej(kt,lfastj_h2o2) = ph_h2o2(iclm,kt,jclm) / ft valuej(kt,lfastj_ch3ooh) = ph_ch3o2h(iclm,kt,jclm)/ft valuej(kt,lfastj_hchoa) = ph_ch2or(iclm,kt,jclm)/ ft valuej(kt,lfastj_hchob) = ph_ch2om(iclm,kt,jclm)/ ft valuej(kt,lfastj_n2o5) = ph_n2o5(iclm,kt,jclm) / ft end do return !finished wrf->peg mapping end subroutine mapJrates_tofrom_host !----------------------------------------------------------------------- !*********************************************************************** subroutine interface_fastj(tmidh,sla,slo,julian_day, & pbnd, psfc, temp, ozone, & dz, cloudmr, airdensity, relhum, zatw, & isvode, jsvode, lpar, jpnl, & curr_secs, valuej, cos_sza, processingAerosols, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !----------------------------------------------------------------- ! sets parameters for fastj call. ! inputs ! tmidh -- GMT time in decimal hours at which to calculate ! photolysis rates ! sla -- latitude, decimal degrees in real*8 ! slo -- negative of the longitude, decimal degrees in real*8 ! julian_day -- day of the year in julian days ! pbnd(0:lpar) = pressure at top boundary of cell k (dynes/cm^2). ! psfc = surface pressure (dynes/cm^2). ! temp(lpar)= mid-cell temperature values (deg K) ! ozone(lpar) = mid-cell ozone mixing ratios ! surface_albedo -- broadband albedo (dimensionless) ! curr_secs -- elapsed time from start of simulation in seconds. ! isvode,jsvode -- current column i,j. ! ! lpar -- vertical extent of column (from module_fastj_cmnh) ! ! outputs ! cos_sza -- cosine of solar zenith angle. ! valuej(lpar,nphoto_fastj-1) -- array of photolysis rates, s-1. ! ! ! local variables ! surface_pressure_mb -- surface pressure (mb). equal to col_press_mb(1). ! col_press_mb(lpar) -- for the column, grid cell boundary pressures ! (not at cell centers) up until the bottom pressure for the ! top cell (mb). ! col_temp_K(lpar+1) -- for the column, grid cell center temperature (deg K) ! col_ozone(lpar+1) -- for the column, grid cell center ozone mixing ! ratios (dimensionless) ! col_optical_depth(lpar+1) -- for the column, grid cell center cloud ! optical depths (dimensionless).SET TO ZERO IN THIS VERSION ! tauaer_550 -- aerosol optical thickness at 550 nm. ! note: photolysis rates are calculated at centers of model layers ! the pressures are given at the boundaries defining ! the top and bottom of the layers ! so the number of pressure values is equal ! to the (number of layers) + 1 ; the last pressure is set = 0 in fastj code. ! pressures from the surface up to the bottom of the top (lpar<=kmaxd) cell ! ******** pressure 2 ! layer 1 - temperature,optical depth, and O3 given here ! ******** pressure 1 ! the optical depth is appropriate for the layer depth ! conversion factor: 1 dyne/cm2 = 0.001 mb !----------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_peg_util, only: peg_message, peg_error_fatal IMPLICIT NONE !jdf integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nphoto_fastj parameter (nphoto_fastj = 14) integer & lfastj_no2, lfastj_o3a, lfastj_o3b, lfastj_h2o2, & lfastj_hchoa, lfastj_hchob, lfastj_ch3ooh, lfastj_no3x, & lfastj_no3l, lfastj_hono, lfastj_n2o5, lfastj_hno3, & lfastj_hno4 parameter( lfastj_no2 = 1 ) parameter( lfastj_o3a = 2 ) parameter( lfastj_o3b = 3 ) parameter( lfastj_h2o2 = 4 ) parameter( lfastj_hchoa = 5 ) parameter( lfastj_hchob = 6 ) parameter( lfastj_ch3ooh= 7 ) parameter( lfastj_no3x = 8 ) parameter( lfastj_no3l = 9 ) parameter( lfastj_hono = 10 ) parameter( lfastj_n2o5 = 11 ) parameter( lfastj_hno3 = 12 ) parameter( lfastj_hno4 = 13 ) integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / !jdf real pbnd(0:lpar), psfc real temp(lpar), ozone(lpar), surface_albedo real dz(lpar), cloudmr(lpar), airdensity(lpar), relhum(lpar), zatw(lpar+1) real(kind=8) :: curr_secs integer isvode, jsvode real cos_sza integer,parameter :: lunout=41 real valuej(lpar,nphoto_fastj) real hl,rhl,factor1,part1,part2,cfrac,rhfrac real emziohl(lpar+1),clwp(lpar) !ec material to check output real valuej_no3rate(lpar) real*8 lat,lon real*8 jvalue(lpar,nphoto_fastj) real sza real tau1 real tmidh, sla, slo integer julian_day,iozone1 integer,parameter :: nfastj_rxns = 14 integer k, l real surface_pressure_mb, tauaer_550, & col_press_mb,col_temp_K,col_ozone,col_optical_depth dimension col_press_mb(lpar+2),col_temp_K(lpar+1), & col_ozone(lpar+1),col_optical_depth(lpar+1) character*80 msg ! define logical processingAerosols ! if processingAerosols = true, uses values calculated in subroutine ! mieaer for variables & arrays in common block mie. ! if processingAerosols = false, sets all variables & arrays in common ! block mie to zero. (JCB-revised Fast-J requires common block mie info, ! regardless of whether aerosols are present or not. Original Wild Fast-J ! did not use common block mie info.) logical processingAerosols ! set lat and longitude as real*8 for consistency with fastj code. ! variables lat and lon previously declared as reals lat = sla lon = slo ! ! cloud optical depths currently treated by using fractional cloudiness ! based on relative humidity. cloudmr set up to use cloud liquid water ! but hooks into microphysics need to be tested - for now set cloudmr=0 ! ! parameters to calculate 'typical' liquid cloudwater path values for ! non convective clouds based on approximations in NCAR's CCM2 ! 0.18 = reference liquid water concentration (gh2o/m3) ! hl = liquid water scale height (m) ! hl=1080.+2000.0*cos(lat*0.017454329) rhl=1.0/hl do k =1, lpar+1 emziohl(k)=exp(-zatw(k)*rhl) enddo do k =1, lpar clwp(k)=0.18*hl*(emziohl(k)-emziohl(k+1)) enddo ! assume radius of cloud droplets is constant at 10 microns (0.001 cm) and ! that density of water is constant at 1 g2ho/cm3 ! factor1=3./2./0.001/1. factor1=1500.0 do k =1, lpar col_optical_depth(k) = 0.0 cfrac=0.0 cloudmr(k)=0.0 if(cloudmr(k).gt.0.0) cfrac=1.0 ! 18.0*airdensity converts mole h2o/mole air to g h2o/cm3, part1 is in g h2o/cm2 part1=cloudmr(k)*cfrac*18.0*airdensity(k)*dz(k)*100.0 if(relhum(k).lt.0.8) then rhfrac=0.0 elseif(relhum(k).le.1.0.and.relhum(k).ge.0.8) then ! rhfrac=(relhum(k)-0.8)/(1.0-0.8) rhfrac=(relhum(k)-0.8)/0.2 else rhfrac=1.0 endif if(rhfrac.ge.0.01) then ! factor 1.0e4 converts clwp of g h2o/m2 to g h2o/cm2 part2=rhfrac*clwp(k)/1.e4 else part2=0.0 endif if(cfrac.gt.0) part2=0.0 col_optical_depth(k) = factor1*(part1+part2) ! col_optical_depth(k) = 0.0 ! if(isvode.eq.33.and.jsvode.eq.34) & ! print*,'jdf opt',isvode,jsvode,k,col_optical_depth(k), & ! cfrac,rhfrac,relhum(k),clwp(k) end do col_optical_depth(lpar+1) = 0.0 if (.not.processingAerosols) then ! set common block mie variables to 0 if call set_common_mie(lpar, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) end if ! processingAerosols=false ! set pressure, temperature, ozone of each cell in the column ! set iozone1 = lpar to allow replacement of climatological ozone with model ! predicted ozone to top of chemistry column; standard fastj climatological o3 ! thereafter. surface_pressure_mb = psfc * 0.001 tau1 = tmidh col_press_mb(1) = psfc * 0.001 iozone1 = lpar do k =1, lpar col_press_mb(k+1) = pbnd(k) * 0.001 col_temp_K(k) = temp(k) col_ozone(k) = ozone(k) end do surface_albedo=0.055 ! set aerosol parameters needed by Fast-J if (processingAerosols) then tauaer_550 = 0.0 ! needed parameters already calculated by subroutine ! mieaer and passed into proper parts of fastj code ! via module_fastj_cmnmie else tauaer_550 = 0.05 ! no aerosols, assume typical constant aerosol optical thickness end if CALL wrf_debug(250,'interface_fastj: calling fastj') call fastj(isvode,jsvode,lat,lon,surface_pressure_mb,surface_albedo, & julian_day, tau1, & col_press_mb, col_temp_K, col_optical_depth, col_ozone, & iozone1,tauaer_550,jvalue,sza,lpar,jpnl, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !cos_sza = cosd(sza) cos_sza = cos(sza*3.141592653/180.) ! array jvalue (real*8) is returned from fastj. array valuej(unspecified ! real, default of real*4) is sent on to ! other chemistry subroutines do k = 1, lpar valuej(k, lfastj_no2) = jvalue(k,lfastj_no2) valuej(k, lfastj_o3a) = jvalue(k,lfastj_o3a) valuej(k, lfastj_o3b) = jvalue(k,lfastj_o3b) valuej(k, lfastj_h2o2) = jvalue(k,lfastj_h2o2) valuej(k, lfastj_hchoa) = jvalue(k,lfastj_hchoa) valuej(k, lfastj_hchob) = jvalue(k,lfastj_hchob) valuej(k, lfastj_ch3ooh) = jvalue(k,lfastj_ch3ooh) valuej(k, lfastj_no3x) = jvalue(k,lfastj_no3x) valuej(k, lfastj_no3l) = jvalue(k,lfastj_no3l) valuej(k, lfastj_hono) = jvalue(k,lfastj_hono) valuej(k, lfastj_n2o5) = jvalue(k,lfastj_n2o5) valuej(k, lfastj_hno3) = jvalue(k,lfastj_hno3) valuej(k, lfastj_hno4) = jvalue(k,lfastj_hno4) end do ! diagnostic output and zeroed value if negative photolysis rates returned do k = 1, lpar valuej(k,nphoto_fastj)=0.0 do l = 1, nphoto_fastj-1 if (valuej(k,l) .lt. 0) then write( msg, '(a,f14.2,4i4,1x,e11.4)' ) & 'FASTJ negative Jrate ' // & 'tsec i j k l J(k,l)', curr_secs,isvode,jsvode,k,l,valuej(k,l) call peg_message( lunerr, msg ) valuej(k,l) = 0.0 ! following code used if want run stopped with negative Jrate ! msg = '*** subr interface_fastj -- ' // & ! 'Negative J rate returned from Fast-J' ! call peg_error_fatal( lunerr, msg ) end if end do end do ! compute overall no3 photolysis rate ! wig: commented out since it is not used anywhere ! do k = 1, lpar ! valuej_no3rate(k) = & ! valuej(k, lfastj_no3x) + valuej(k,lfastj_no3l) ! end do end subroutine interface_fastj !*********************************************************************** subroutine set_common_mie(lpar, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) ! for use when aerosols are not included in a model run. sets variables ! in common block mie to zero, except for wavelengths. ! OUTPUT: in module_fastj_cmnmie ! wavmid ! fast-J wavelengths (cm) ! tauaer ! aerosol optical depth ! waer ! aerosol single scattering albedo ! gaer ! aerosol asymmetery factor ! extaer ! aerosol extinction ! l2,l3,l4,l5,l6,l7 ! Legendre coefficients, numbered 0,1,2,...... ! sizeaer ! average wet radius ! INPUTS ! lpar = total number of vertical layers in chemistry model. Passed ! via module_fastf_cmnh ! NB = total vertical layers + 1 considered by FastJ (lpar+1=kmaxd+1). ! passed via module_fast j_data !------------------------------------------------------------------------ USE module_data_mosaic_other, only : kmaxd IMPLICIT NONE !jdf integer lpar ! Number of levels in CTM integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / !jdf ! LOCAL VARIABLES integer klevel ! vertical level index integer ns ! spectral loop index ! aerosol optical properties: set everything = 0 when no aerosol do 1000 ns=1,nspint do 1000 klevel = 1, lpar tauaer(ns,klevel)=0. waer(ns,klevel)=0. gaer(ns,klevel)=0. sizeaer(ns,klevel)=0.0 extaer(ns,klevel)=0.0 l2(ns,klevel)=0.0 l3(ns,klevel)=0.0 l4(ns,klevel)=0.0 l5(ns,klevel)=0.0 l6(ns,klevel)=0.0 l7(ns,klevel)=0.0 1000 continue return end subroutine set_common_mie !*********************************************************************** subroutine fastj(isvode,jsvode,lat,lon,surface_pressure,surface_albedo, & julian_day,tau1, pressure, temperature, optical_depth, my_ozone1, & iozone1,tauaer_550_1,jvalue,SZA_dum,lpar,jpnl, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) ! input: ! lat = latitute; must be real*8 ! lon = longitude; must be real*8 ! surface_pressure (mb); real*4 ! surface_albedo (broadband albedo); real*4 ! julian_day; integer ! tau1 = time of calculation (GMT); real*4 ! pressure (mb) = vector of pressure values, pressure(NB); ! real*4; NB is the number of model layers; ! pressure (NB+1) is defined as 0 mb in model ! temperature (degree K)= vector of temperature values, temperature(NB); ! real*4 ! optical_depth (dimensionless) = vector of cloud optical depths, ! optical_depth(NB); real*4 ! my_ozone1 (volume mixing ratio) = ozone at layer center ! ozone(iozone); real*4; if iozone1 <= NB-1, then climatology is ! used in the upper layers ! tauaer_550; real*4 aerosol optical thickness @ 550 nm ! input note: NB is the number of model layers -- photolysis rates are calculated ! at layer centers while pressures are given at the boundaries defining ! the top and bottom of the layers. The number of pressure values = ! (number of layers) + 1 ; see below ! ******** pressure 2 ! layer 1 - optical depth, O3, and temperature given here ! ******** pressure 1 ! temperature and o3 are defined at the layer center. optical depth is ! appropriate for the layer depth. ! output: ! jvalue = photolysis rates, an array of dimension jvalue(jpnl,jppj) where ! jpnl = # of models level at which photolysis rates are calculated ! note: level 1 = first level of model (adjacent to ground) ! jppj = # of chemical species for which photolysis rates are calculated; ! this is fixed and is not easy to change on the fly ! jpnl land jppl are defined in the common block "cmn_h.f" ! SZA_dum = solar zenith angle !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 ! following specific for fastj ! jppj will be gas phase mechanism dependent integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution ! The vertical level variables are set in fastj_driver. integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry ! following should be available from other wrf modules and passed into ! photodriver real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year real(kind=double), dimension(ipar_fastj,jpar) :: P , SA real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD real(kind=double) my_ozone(kmaxd) !real*8 version of ozone mixing ratios real tauaer_550 integer iozone integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point save :: nslat, nslon integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / !jdf integer julian_day real surface_pressure,surface_albedo,pressure(lpar+2), & temperature(lpar+1) real optical_depth(lpar+1) real tau1 real*8 pi_fastj,lat,lon,timej,jvalue(lpar,jppj) integer isvode,jsvode integer iozone1,i real my_ozone1(lpar+1) real tauaer_550_1 real sza_dum integer ientryno_fastj save ientryno_fastj data ientryno_fastj / 0 / ! Just focus on one column nslat = 1 nslon = 1 pi_fastj=3.141592653589793D0 ! ! JCB - note that pj(NB+1) = p and is defined such elsewhere ! wig: v3.0.0 release has do loop from 1:lpar+1 but T is never ! initialized at lpar+1 leading to model instabilities. T and OD ! are ultimately never used at lpar+1. So, I changed this loop to ! 1:lpar and then just assign pj a value at lpar+1. do i=1,lpar pj(i)=pressure(i) T(nslon,nslat,i)=temperature(i) OD(nslon,nslat,i)=optical_depth(i) enddo pj(lpar+1) = pressure(i) ! surface albedo SA(nslon,nslat)=surface_albedo ! iozone=iozone1 do i=1,iozone1 my_ozone(i)=my_ozone1(i) enddo ! tau_fastj=tau1 ! fix time iday_fastj=julian_day ! fix optical depth for situations where aerosols not considered tauaer_550=tauaer_550_1 ! month_fastj=int(dble(iday_fastj)*12.d0/365.d0)+1 ! Approximately xgrd(nslon)=lon*pi_fastj/180.d0 ygrd(nslat)=lat*pi_fastj/180.d0 ydgrd(nslat)=lat ! Initial call to Fast-J to set things up--done only once if (ientryno_fastj .eq. 0) then call inphot2 ientryno_fastj = 1 end if ! ! Now call fastj as appropriate timej=0.0 ! manually set offset to zero (JCB: 14 November 2001) call photoj(isvode,jsvode,jvalue,timej,nslat,nslon,iozone,tauaer_550, & my_ozone,p,t,od,sa,lpar,jpnl, & xgrd,ygrd,tau_fastj,month_fastj,iday_fastj,ydgrd, & sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) sza_dum=SZA return end subroutine fastj !********************************************************************** !---(pphot.f)-------generic CTM shell from UCIrvine (p-code 4.0, 7/99) !---------PPHOT calculates photolysis rates with the Fast-J scheme !---------subroutines: inphot, photoj, Fast-J schemes... !----------------------------------------------------------------------- ! subroutine inphot2 !----------------------------------------------------------------------- ! Routine to initialise photolysis rate data, called directly from the ! cinit routine in ASAD. Currently use it to read the JPL spectral data ! and standard O3 and T profiles and to set the appropriate reaction index. !----------------------------------------------------------------------- ! ! iph Channel number for reading all data files ! rad Radius of Earth (cm) ! zzht Effective scale height above top of atmosphere (cm) ! dtaumax Maximum opt.depth above which a new level should be inserted ! dtausub No. of opt.depths at top of cloud requiring subdivision ! dsubdiv Number of additional levels to add at top of cloud ! szamax Solar zenith angle cut-off, above which to skip calculation ! !----------------------------------------------------------------------- ! USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year !jdf ! ! Set labels of photolysis rates required !ec032504 CALL RD_JS(iph,path_fastj_ratjd) ! call rd_js2 ! ! Read in JPL spectral data set !ec032504 CALL RD_TJPL(iph,path_fastj_jvspec) call rd_tjpl2 ! ! Read in T & O3 climatology !ec032504 CALL RD_PROF(iph,path_fastj_jvatms) ! call rd_prof2 ! ! Select Aerosol/Cloud types to be used ! call set_aer2 return end subroutine inphot2 !************************************************************************* subroutine photoj(isvode,jsvode,zpj,timej,nslat,nslon,iozone,tauaer_550_1, & my_ozone,p,t,od,sa,lpar,jpnl,xgrd,ygrd,tau_fastj,month_fastj,iday_fastj, & ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !----------------------------------------------------------------------- !----jv_trop.f: new FAST J-Value code, troposphere only (mjprather 6/96) !---- uses special wavelength quadrature spectral data (jv_spec.dat) !--- that includes only 289 nm - 800 nm (later a single 205 nm add-on) !--- uses special compact Mie code based on Feautrier/Auer/Prather vers. !----------------------------------------------------------------------- ! ! zpj External array providing J-values to main CTM code ! timej Offset in hours from start of timestep to time J-values ! required for - take as half timestep for mid-step Js. ! solf Solar distance factor, for scaling; normally given by: ! 1.0-(0.034*cos(real(iday_fastj-186)*2.0*pi_fastj/365.)) ! !----------basic common blocks:----------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year real(kind=double), dimension(ipar_fastj,jpar) :: P , SA real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD real(kind=double) my_ozone(kmaxd) !real*8 version of ozone mixing ratios real tauaer_550_1 integer iozone integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / !jdf real*8 zpj(lpar,jppj),timej,solf real*8 pi_fastj integer i,j integer isvode,jsvode !----------------------------------------------------------------------- ! do i=1,jpnl do j=1,jppj zj(i,j)=0.D0 zpj(i,j)=0.D0 enddo enddo ! !---Calculate new solar zenith angle CALL SOLAR2(timej,nslat,nslon, & xgrd,ygrd,tau_fastj,month_fastj,iday_fastj) if(SZA.gt.szamax) go to 10 ! !---Set up profiles on model levels CALL SET_PROF(isvode,jsvode,nslat,nslon,iozone,tauaer_550_1, & my_ozone,p,t,od,sa,lpar,jpnl,month_fastj,ydgrd) ! !---Print out atmosphere if(iprint.ne.0)CALL PRTATM(3,nslat,nslon,tau_fastj,month_fastj,ydgrd) ! code change jcb ! call prtatm(0) ! !----------------------------------------------------------------------- CALL JVALUE(isvode,jsvode,lpar,jpnl, & ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !----------------------------------------------------------------------- !---Print solar flux terms ! WRITE(6,'(A16,I5,20I9)') ' wave (beg/end)',(i,i=1,jpnl) ! DO j=NW1,NW2 ! WRITE(6,'(2F8.2,20F9.6)') WBIN(j),WBIN(j+1), ! $ (FFF(j,i)/FL(j),i=1,jpnl) ! ENDDO ! !---Include variation in distance to sun pi_fastj=3.1415926536d0 solf=1.d0-(0.034d0*cos(dble(iday_fastj-186)*2.d0 & *pi_fastj/365.d0)) if(iprint.ne.0)then ! write(6,'('' solf = '', f10.5)')solf write(*,'('' solf = '', f10.5)')solf endif ! solf=1.d0 ! code change jcb !----------------------------------------------------------------------- CALL JRATET(solf,nslat,nslon,p,t,od,sa,lpar,jpnl) !----------------------------------------------------------------------- ! ! "zj" updated in JRATET - pass this back to ASAD as "zpj" do i=1,jpnl do j=1,jppj zpj(i,j)= zj(i,j) enddo enddo ! !---Output selected values 10 if((.not.ldeg45.and.nslon.eq.37.and.nslat.eq.36).or. & (ldeg45.and.nslon.eq.19.and.nslat.eq.18)) then i=min(jppj,8) ! write(6,1000)iday_fastj,tau_fastj+timej,sza,jlabel(i),zpj(1,i) endif ! return ! 1000 format(' Photolysis on day ',i4,' at ',f4.1,' hrs: SZA = ',f7.3, & ! ' J',a7,'= ',1pE10.3) end subroutine photoj !************************************************************************* subroutine set_prof(isvode,jsvode,nslat,nslon,iozone,tauaer_550, & my_ozone,p,t,od,sa,lpar,jpnl,month_fastj,ydgrd) !----------------------------------------------------------------------- ! Routine to set up atmospheric profiles required by Fast-J using a ! doubled version of the level scheme used in the CTM. First pressure ! and z* altitude are defined, then O3 and T are taken from the supplied ! climatology and integrated to the CTM levels (may be overwritten with ! values directly from the CTM, if desired) and then black carbon and ! aerosol profiles are constructed. ! Oliver (04/07/99) !----------------------------------------------------------------------- ! ! pj Pressure at boundaries of model levels (hPa) ! z Altitude of boundaries of model levels (cm) ! odcol Optical depth at each model level ! masfac Conversion factor for pressure to column density ! ! TJ Temperature profile on model grid ! DM Air column for each model level (molecules.cm-2) ! DO3 Ozone column for each model level (molecules.cm-2) ! DBC Mass of Black Carbon at each model level (g.cm-3) ! .....! ! PSTD Approximate pressures of levels for supplied climatology ! !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year real(kind=double), dimension(ipar_fastj,jpar) :: P , SA real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD real(kind=double) my_ozone(kmaxd) !real*8 version of ozone mixing ratios real tauaer_550 integer iozone integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point !jdf real*8 pstd(52),oref2(51),tref2(51),bref2(51) real*8 odcol(lpar),dlogp,f0,t0,b0,pb,pc,xc,masfac,scaleh real vis, aerd1, aerd2 integer i, k, l, m integer isvode,jsvode pj(NB+1) = 0.d0 ! define top level ! ! Set up cloud and surface properties call CLDSRF(isvode,jsvode,odcol,nslat,nslon,p,t,od,sa,lpar,jpnl) ! Mass factor - delta-Pressure (mbars) to delta-Column (molecules.cm-2) masfac=100.d0*6.022d+23/(28.97d0*9.8d0*10.d0) ! ! Set up pressure levels for O3/T climatology - assume that value ! given for each 2 km z* level applies from 1 km below to 1 km above, ! so select pressures at these boundaries. Surface level values at ! 1000 mb are assumed to extend down to the actual P(nslon,nslat). ! pstd(1) = max(pj(1),1000.d0) pstd(2) = 1000.d0*10.d0**(-1.d0/16.d0) dlogp = 10.d0**(-2.d0/16.d0) do i=3,51 pstd(i) = pstd(i-1)*dlogp enddo pstd(52) = 0.d0 ! ! Select appropriate monthly and latitudinal profiles m = max(1,min(12,month_fastj)) l = max(1,min(18,(int(ydgrd(nslat))+99)/10)) ! ! Temporary arrays for climatology data do i=1,51 oref2(i)=oref(i,l,m) tref2(i)=tref(i,l,m) bref2(i)=bref(i) enddo ! ! Apportion O3 and T on supplied climatology z* levels onto CTM levels ! with mass (pressure) weighting, assuming constant mixing ratio and ! temperature half a layer on either side of the point supplied. ! do i = 1,NB F0 = 0.d0 T0 = 0.d0 B0 = 0.d0 do k = 1,51 PC = min(pj(i),pstd(k)) PB = max(pj(i+1),pstd(k+1)) if(PC.gt.PB) then XC = (PC-PB)/(pj(i)-pj(i+1)) F0 = F0 + oref2(k)*XC T0 = T0 + tref2(k)*XC B0 = B0 + bref2(k)*XC endif enddo TJ(i) = T0 DO3(i)= F0*1.d-6 DBC(i)= B0 enddo ! ! Insert model values here to replace or supplement climatology. ! Note that CTM temperature is always used in x-section calculations ! (see JRATET); TJ is used in actinic flux calculation only. ! do i=1,lpar ! JCB code change; just use climatlogy for upper levels if(i.le.iozone)DO3(i) = my_ozone(i) ! Volume Mixing Ratio ! TJ(i) = T(nslon,nslat,I) ! Kelvin ! JCB - overwrite climatology ! TJ(i) = (T(nslon,nslat,i)+T(nslon,nslat,i+1))/2. ! JCB - take midpoint ! code change - now take temperature as appropriate for midpoint of layer TJ(i)=T(nslon,nslat,i) enddo if(lpar+1.le.iozone)then DO3(lpar+1) = my_ozone(lpar+1) ! Above top of model (or use climatology) endif ! TJ(lpar+1) = my_temp(lpar) ! Above top of model (or use climatology) !wig 26-Aug-2000: Comment out following line so that climatology is used for ! above the model top. ! TJ(lpar+1) = T(nslon,nslat,NB) ! JCB - just use climatology or given temperature ! JCB read in O3 ! ! ! Calculate effective altitudes using scale height at each level z(1) = 0.d0 do i=1,lpar scaleh=1.3806d-19*masfac*TJ(i) z(i+1) = z(i)-(log(pj(i+1)/pj(i))*scaleh) enddo ! ! Add Aerosol Column - include aerosol types here. Currently use soot ! water and ice; assume black carbon x-section of 10 m2/g, independent ! of wavelength; assume limiting temperature for ice of -40 deg C. do i=1,lpar ! AER(1,i) = DBC(i)*10.d0*(z(i+1)-z(i)) ! DBC must be g/m^3 ! calculate AER(1,i) according to aerosol density - use trap rule vis=23.0 call aeroden(z(i)/100000.,vis,aerd1) ! convert cm to km call aeroden(z(i+1)/100000.,vis,aerd2) ! trap rule used here; convert cm to km; divide by 100000. AER(1,i)=(z(i+1)-z(i))/100000.*(aerd1+aerd2)/2./4287.55*tauaer_550 ! write(6,*)i,z(i)/100000.,aerd1,aerd2,tauaer_550,AER(1,i) ! if(T(nslon,nslat,I).gt.233.d0) then AER(2,i) = odcol(i) AER(3,i) = 0.d0 else AER(2,i) = 0.d0 AER(3,i) = odcol(i) endif enddo do k=1,MX AER(k,lpar+1) = 0.d0 ! just set equal to zero enddo AER(1,lpar+1)=2.0*AER(1,lpar) ! kludge ! ! Calculate column quantities for Fast-J do i=1,NB DM(i) = (PJ(i)-PJ(i+1))*masfac DO3(i) = DO3(i)*DM(i) enddo ! end subroutine set_prof !****************************************************************** ! SUBROUTINE CLDSRF(odcol) SUBROUTINE CLDSRF(isvode,jsvode,odcol,nslat,nslon,p,t,od,sa, & lpar,jpnl) !----------------------------------------------------------------------- ! Routine to set cloud and surface properties !----------------------------------------------------------------------- ! rflect Surface albedo (Lambertian) ! odmax Maximum allowed optical depth, above which they are scaled ! odcol Optical depth at each model level ! odsum Column optical depth ! nlbatm Level of lower photolysis boundary - usually surface !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point real(kind=double), dimension(ipar_fastj,jpar) :: P , SA real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD !jdf integer i, j, k integer isvode, jsvode real*8 odcol(lpar), odsum, odmax, odtot ! ! Default lower photolysis boundary as bottom of level 1 nlbatm = 1 ! ! Set surface albedo RFLECT = dble(SA(nslon,nslat)) RFLECT = max(0.d0,min(1.d0,RFLECT)) ! ! Zero aerosol column do k=1,MX do i=1,NB AER(k,i) = 0.d0 enddo enddo ! ! Scale optical depths as appropriate - limit column to 'odmax' odmax = 200.d0 odsum = 0.d0 do i=1,lpar odcol(i) = dble(OD(nslon,nslat,i)) odsum = odsum + odcol(i) ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf odcol',odcol(i),odsum enddo if(odsum.gt.odmax) then odsum = odmax/odsum do i=1,lpar odcol(i) = odcol(i)*odsum enddo odsum = odmax endif ! Set sub-division switch if appropriate odtot=0.d0 jadsub(nb)=0 jadsub(nb-1)=0 do i=nb-1,1,-1 k=2*i jadsub(k)=0 jadsub(k-1)=0 odtot=odtot+odcol(i) if(odcol(i).gt.0.d0.and.dtausub.gt.0.d0) then if(odtot.le.dtausub) then jadsub(k)=1 jadsub(k-1)=1 ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in cldsf1',i,k,jadsub(k) elseif(odtot.gt.dtausub) then jadsub(k)=1 jadsub(k-1)=0 do j=1,2*(i-1) jadsub(j)=0 enddo ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in cldsf2',i,k,jadsub(k) go to 20 endif endif enddo 20 continue ! return end SUBROUTINE CLDSRF !******************************************************************** subroutine solar2(timej,nslat,nslon, & xgrd,ygrd,tau_fastj,month_fastj,iday_fastj) !----------------------------------------------------------------------- ! Routine to set up SZA for given lat, lon and time !----------------------------------------------------------------------- ! timej Offset in hours from start of timestep to time J-values ! required for - take as half timestep for mid-step Js. !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point !jdf real*8 pi_fastj, pi180, loct, timej real*8 sindec, soldek, cosdec, sinlat, sollat, coslat, cosz ! pi_fastj=3.141592653589793D0 pi180=pi_fastj/180.d0 sindec=0.3978d0*sin(0.9863d0*(dble(iday_fastj)-80.d0)*pi180) soldek=asin(sindec) cosdec=cos(soldek) sinlat=sin(ygrd(nslat)) sollat=asin(sinlat) coslat=cos(sollat) ! loct = (((tau_fastj+timej)*15.d0)-180.d0)*pi180 + xgrd(nslon) cosz = cosdec*coslat*cos(loct) + sindec*sinlat sza = acos(cosz)/pi180 U0 = cos(SZA*pi180) ! return end subroutine solar2 !********************************************************************** SUBROUTINE JRATET(SOLF,nslat,nslon,p,t,od,sa,lpar,jpnl) !----------------------------------------------------------------------- ! Calculate and print J-values. Note that the loop in this routine ! only covers the jpnl levels actually needed by the CTM. !----------------------------------------------------------------------- ! ! FFF Actinic flux at each level for each wavelength bin ! QQQ Cross sections for species (read in in RD_TJPL) ! SOLF Solar distance factor, for scaling; normally given by: ! 1.0-(0.034*cos(real(iday_fastj-186)*2.0*pi_fastj/365.)) ! Assumes aphelion day 186, perihelion day 3. ! TQQ Temperatures at which QQQ cross sections supplied ! !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point real(kind=double), dimension(ipar_fastj,jpar) :: P , SA real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD !jdf integer i, j, k real*8 qo2tot, qo3tot, qo31d, qo33p, qqqt ! real*8 xseco2, xseco3, xsec1d, solf, tfact real*8 solf, tfact ! do I=1,jpnl VALJ(1) = 0.d0 VALJ(2) = 0.d0 VALJ(3) = 0.d0 do K=NW1,NW2 ! Using model 'T's here QO2TOT= xseco2(K,dble(T(nslon,nslat,I))) VALJ(1) = VALJ(1) + QO2TOT*FFF(K,I) QO3TOT= xseco3(K,dble(T(nslon,nslat,I))) QO31D = xsec1d(K,dble(T(nslon,nslat,I)))*QO3TOT QO33P = QO3TOT - QO31D VALJ(2) = VALJ(2) + QO33P*FFF(K,I) VALJ(3) = VALJ(3) + QO31D*FFF(K,I) enddo !------Calculate remaining J-values with T-dep X-sections do J=4,NJVAL VALJ(J) = 0.d0 TFACT = 0.d0 if(TQQ(2,J).gt.TQQ(1,J)) TFACT = max(0.d0,min(1.d0, & (T(nslon,nslat,I)-TQQ(1,J))/(TQQ(2,J)-TQQ(1,J)) )) do K=NW1,NW2 QQQT = QQQ(K,1,J-3) + (QQQ(K,2,J-3) - QQQ(K,1,J-3))*TFACT VALJ(J) = VALJ(J) + QQQT*FFF(K,I) !------Additional code for pressure dependencies ! if(jpdep(J).ne.0) then ! VALJ(J) = VALJ(J) + QQQT*FFF(K,I)* ! $ (zpdep(K,L)*(pj(i)+pj(i+1))*0.5d0) ! endif enddo enddo do j=1,jppj zj(i,j)=VALJ(jind(j))*jfacta(j)*SOLF enddo ! Herzberg bin do j=1,nhz zj(i,hzind(j))=hztoa(j)*fhz(i)*SOLF enddo enddo return end SUBROUTINE JRATET !********************************************************************* SUBROUTINE PRTATM(N,nslat,nslon,tau_fastj,month_fastj,ydgrd) !----------------------------------------------------------------------- ! Print out the atmosphere and calculate appropriate columns ! N=1 Print out column totals only ! N=2 Print out full columns ! N=3 Print out full columns and climatology !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nslat ! Latitude of current profile point integer nslon ! Longitude of current profile point !jdf integer n, i, k, l, m real*8 COLO3(NB),COLO2(NB),COLAX(MX,NB),ZKM,ZSTAR,PJC real*8 climat(9),masfac,dlogp if(N.eq.0) return !---Calculate columns, for diagnostic output only: COLO3(NB) = DO3(NB) COLO2(NB) = DM(NB)*0.20948d0 do K=1,MX COLAX(K,NB) = AER(K,NB) enddo do I=NB-1,1,-1 COLO3(i) = COLO3(i+1)+DO3(i) COLO2(i) = COLO2(i+1)+DM(i)*0.20948d0 do K=1,MX COLAX(k,i) = COLAX(k,i+1)+AER(k,i) enddo enddo write(*,1200) ' Tau=',tau_fastj,' SZA=',sza write(*,1200) ' O3-column(DU)=',COLO3(1)/2.687d16, & ' column aerosol @1000nm=',(COLAX(K,1),K=1,MX) !---Print out atmosphere if(N.gt.1) then write(*,1000) (' AER-X ','col-AER',k=1,mx) do I=NB,1,-1 PJC = PJ(I) ZKM =1.d-5*Z(I) ZSTAR = 16.d0*DLOG10(1013.d0/PJC) write(*,1100) I,ZKM,ZSTAR,DM(I),DO3(I),1.d6*DO3(I)/DM(I), & TJ(I),PJC,COLO3(I),COLO2(I),(AER(K,I),COLAX(K,I),K=1,MX) enddo endif ! !---Print out climatology if(N.gt.2) then do i=1,9 climat(i)=0.d0 enddo m = max(1,min(12,month_fastj)) l = max(1,min(18,(int(ydgrd(nslat))+99)/10)) masfac=100.d0*6.022d+23/(28.97d0*9.8d0*10.d0) write(*,*) 'Specified Climatology' write(*,1000) do i=51,1,-1 dlogp=10.d0**(-1.d0/16.d0) PJC = 1000.d0*dlogp**(2*i-2) climat(1) = 16.d0*DLOG10(1000.D0/PJC) climat(2) = climat(1) climat(3) = PJC*(1.d0/dlogp-dlogp)*masfac if(i.eq.1) climat(3)=PJC*(1.d0-dlogp)*masfac climat(4)=climat(3)*oref(i,l,m)*1.d-6 climat(5)=oref(i,l,m) climat(6)=tref(i,l,m) climat(7)=PJC climat(8)=climat(8)+climat(4) climat(9)=climat(9)+climat(3)*0.20948d0 write(*,1100) I,(climat(k),k=1,9) enddo write(*,1200) ' O3-column(DU)=',climat(8)/2.687d16 endif return 1000 format(5X,'Zkm',3X,'Z*',8X,'M',8X,'O3',6X,'f-O3',5X,'T',7X,'P',6x, & 'col-O3',3X,'col-O2',2X,10(a7,2x)) 1100 format(1X,I2,0P,2F6.2,1P,2E10.3,0P,F7.3,F8.2,F10.4,1P,10E9.2) 1200 format(A,F8.1,A,10(1pE10.3)) end SUBROUTINE PRTATM SUBROUTINE JVALUE(isvode,jsvode,lpar,jpnl, & ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !----------------------------------------------------------------------- ! Calculate the actinic flux at each level for the current SZA value. ! quit when SZA > 98.0 deg ==> tangent height = 63 km ! or 99. 80 km !----------------------------------------------------------------------- ! ! AVGF Attenuation of beam at each level for each wavelength ! FFF Actinic flux at each desired level ! FHZ Actinic flux in Herzberg bin ! WAVE Effective wavelength of each wavelength bin ! XQO2 Absorption cross-section of O2 ! XQO3 Absorption cross-section of O3 ! !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data USE module_peg_util, only: peg_message !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year real(kind=double), dimension(ipar_fastj,jpar) :: P , SA real(kind=double), dimension(ipar_fastj,jpar,kmaxd+1) :: T, OD integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / !jdf integer j, k ! real*8 wave, xseco3, xseco2 real*8 wave real*8 AVGF(lpar),XQO3(NB),XQO2(NB) ! diagnostics for error situations ! integer lunout ! parameter (lunout = 41) integer isvode,jsvode character*80 msg ! do J=1,jpnl do K=NW1,NW2 FFF(K,J) = 0.d0 enddo FHZ(J) = 0.d0 enddo ! !---SZA check if(SZA.gt.szamax) GOTO 99 ! !---Calculate spherical weighting functions CALL SPHERE(lpar) ! !---Loop over all wavelength bins do K=NW1,NW2 WAVE = WL(K) do J=1,NB XQO3(J) = xseco3(K,dble(TJ(J))) enddo do J=1,NB XQO2(J) = xseco2(K,dble(TJ(J))) enddo !----------------------------------------- CALL OPMIE(isvode,jsvode,K,WAVE,XQO2,XQO3,AVGF,lpar,jpnl, & ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !----------------------------------------- do J=1,jpnl FFF(K,J) = FFF(K,J) + FL(K)*AVGF(J) ! diagnostic if ( FFF(K,J) .lt. 0) then write( msg, '(a,2i4,e14.6)' ) & 'FASTJ neg actinic flux ' // & 'k j FFF(K,J) ', k, j, fff(k,j) call peg_message( lunerr, msg ) end if ! end diagnostic enddo enddo ! !---Herzberg continuum bin above 10 km, if required if(NHZ.gt.0) then K=NW2+1 WAVE = 204.d0 do J=1,NB XQO3(J) = HZO3 XQO2(J) = HZO2 enddo CALL OPMIE(isvode,jsvode,K,WAVE,XQO2,XQO3,AVGF,lpar,jpnl, & ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) do J=1,jpnl if(z(j).gt.1.d6) FHZ(J)=AVGF(J) enddo endif ! 99 continue 1000 format(' SZA=',f6.1,' Reflectvty=',f6.3,' OD=',10(1pe10.3)) return end SUBROUTINE JVALUE FUNCTION xseco3(K,TTT) !----------------------------------------------------------------------- ! Cross-sections for O3 for all processes interpolated across 3 temps !----------------------------------------------------------------------- USE module_fastj_data integer k ! real*8 ttt, flint, xseco3 real*8 ttt, xseco3 xseco3 = & flint(TTT,TQQ(1,2),TQQ(2,2),TQQ(3,2),QO3(K,1),QO3(K,2),QO3(K,3)) return end FUNCTION xseco3 FUNCTION xsec1d(K,TTT) !----------------------------------------------------------------------- ! Quantum yields for O3 --> O2 + O(1D) interpolated across 3 temps !----------------------------------------------------------------------- USE module_fastj_data integer k ! real*8 ttt, flint, xsec1d real*8 ttt, xsec1d xsec1d = & flint(TTT,TQQ(1,3),TQQ(2,3),TQQ(3,3),Q1D(K,1),Q1D(K,2),Q1D(K,3)) return end FUNCTION xsec1d FUNCTION xseco2(K,TTT) !----------------------------------------------------------------------- ! Cross-sections for O2 interpolated across 3 temps; No S_R Bands yet! !----------------------------------------------------------------------- USE module_fastj_data integer k ! real*8 ttt, flint, xseco2 real*8 ttt, xseco2 xseco2 = & flint(TTT,TQQ(1,1),TQQ(2,1),TQQ(3,1),QO2(K,1),QO2(K,2),QO2(K,3)) return end FUNCTION xseco2 REAL*8 FUNCTION flint (TINT,T1,T2,T3,F1,F2,F3) !----------------------------------------------------------------------- ! Three-point linear interpolation function !----------------------------------------------------------------------- real*8 TINT,T1,T2,T3,F1,F2,F3 IF (TINT .LE. T2) THEN IF (TINT .LE. T1) THEN flint = F1 ELSE flint = F1 + (F2 - F1)*(TINT -T1)/(T2 -T1) ENDIF ELSE IF (TINT .GE. T3) THEN flint = F3 ELSE flint = F2 + (F3 - F2)*(TINT -T2)/(T3 -T2) ENDIF ENDIF return end FUNCTION flint SUBROUTINE SPHERE(lpar) !----------------------------------------------------------------------- ! Calculation of spherical geometry; derive tangent heights, slant path ! lengths and air mass factor for each layer. Not called when ! SZA > 98 degrees. Beyond 90 degrees, include treatment of emergent ! beam (where tangent height is below altitude J-value desired at). !----------------------------------------------------------------------- ! ! GMU MU, cos(solar zenith angle) ! RZ Distance from centre of Earth to each point (cm) ! RQ Square of radius ratios ! TANHT Tangent height for the current SZA ! XL Slant path between points ! AMF Air mass factor for slab between level and level above ! !----------------------------------------------------------------------- USE module_fastj_data !jdf integer lpar !jdf integer i, j, k, ii real*8 airmas, gmu, xmu1, xmu2, xl, diff REAL*8 Ux,H,RZ(NB),RQ(NB),ZBYR ! ! Inlined air mass factor function for top of atmosphere AIRMAS(Ux,H) = (1.0d0+H)/SQRT(Ux*Ux+2.0d0*H*(1.0d0- & 0.6817d0*EXP(-57.3d0*ABS(Ux)/SQRT(1.0d0+5500.d0*H))/ & (1.0d0+0.625d0*H))) ! GMU = U0 RZ(1)=RAD+Z(1) ZBYR = ZZHT/RAD DO 2 II=2,NB RZ(II) = RAD + Z(II) RQ(II-1) = (RZ(II-1)/RZ(II))**2 2 CONTINUE IF (GMU.LT.0.0D0) THEN TANHT = RZ(nlbatm)/DSQRT(1.0D0-GMU**2) ELSE TANHT = RZ(nlbatm) ENDIF ! ! Go up from the surface calculating the slant paths between each level ! and the level above, and deriving the appropriate Air Mass Factor DO 16 J=1,NB DO K=1,NB AMF(K,J)=0.D0 ENDDO ! ! Air Mass Factors all zero if below the tangent height IF (RZ(J).LT.TANHT) GOTO 16 ! Ascend from layer J calculating AMFs XMU1=ABS(GMU) DO 12 I=J,lpar XMU2=DSQRT(1.0D0-RQ(I)*(1.0D0-XMU1**2)) XL=RZ(I+1)*XMU2-RZ(I)*XMU1 AMF(I,J)=XL/(RZ(I+1)-RZ(I)) XMU1=XMU2 12 CONTINUE ! Use function and scale height to provide AMF above top of model AMF(NB,J)=AIRMAS(XMU1,ZBYR) ! ! Twilight case - Emergent Beam IF (GMU.GE.0.0D0) GOTO 16 XMU1=ABS(GMU) ! Descend from layer J DO 14 II=J-1,1,-1 DIFF=RZ(II+1)*DSQRT(1.0D0-XMU1**2)-RZ(II) if(II.eq.1) DIFF=max(DIFF,0.d0) ! filter ! Tangent height below current level - beam passes through twice IF (DIFF.LT.0.0D0) THEN XMU2=DSQRT(1.0D0-(1.0D0-XMU1**2)/RQ(II)) XL=ABS(RZ(II+1)*XMU1-RZ(II)*XMU2) AMF(II,J)=2.d0*XL/(RZ(II+1)-RZ(II)) XMU1=XMU2 ! Lowest level intersected by emergent beam ELSE XL=RZ(II+1)*XMU1*2.0D0 ! WTING=DIFF/(RZ(II+1)-RZ(II)) ! AMF(II,J)=(1.0D0-WTING)*2.D0**XL/(RZ(II+1)-RZ(II)) AMF(II,J)=XL/(RZ(II+1)-RZ(II)) GOTO 16 ENDIF 14 CONTINUE ! 16 CONTINUE RETURN END SUBROUTINE SPHERE SUBROUTINE OPMIE(isvode,jsvode,KW,WAVEL,XQO2,XQO3,FMEAN,lpar,jpnl, & ydgrd,sizeaer,extaer,waer,gaer,tauaer,l2,l3,l4,l5,l6,l7) !----------------------------------------------------------------------- ! NEW Mie code for J's, only uses 8-term expansion, 4-Gauss pts ! Currently allow up to NP aerosol phase functions (at all altitudes) to ! be associated with optical depth AER(1:NC) = aerosol opt.depth @ 1000 nm ! ! Pick Mie-wavelength with phase function and Qext: ! ! 01 RAYLE = Rayleigh phase ! 02 ISOTR = isotropic ! 03 ABSRB = fully absorbing 'soot', wavelength indep. ! 04 S_Bkg = backgrnd stratospheric sulfate (n=1.46,log-norm:r=.09um/sigma=.6) ! 05 S_Vol = volcanic stratospheric sulfate (n=1.46,log-norm:r=.08um/sigma=.8) ! 06 W_H01 = water haze (H1/Deirm.) (n=1.335, gamma: r-mode=0.1um /alpha=2) ! 07 W_H04 = water haze (H1/Deirm.) (n=1.335, gamma: r-mode=0.4um /alpha=2) ! 08 W_C02 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=2.0um /alpha=6) ! 09 W_C04 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=4.0um /alpha=6) ! 10 W_C08 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=8.0um /alpha=6) ! 11 W_C13 = water cloud (C1/Deirm.) (n=1.335, gamma: r-mode=13.3um /alpha=6) ! 12 W_L06 = water cloud (Lacis) (n=1.335, r-mode=5.5um / alpha=11/3) ! 13 Ice-H = hexagonal ice cloud (Mishchenko) ! 14 Ice-I = irregular ice cloud (Mishchenko) ! ! Choice of aerosol index MIEDX is made in SET_AER; optical depths are ! apportioned to the AER array in SET_PROF ! !----------------------------------------------------------------------- ! FUNCTION RAYLAY(WAVE)---RAYLEIGH CROSS-SECTION for wave > 170 nm ! WSQI = 1.E6/(WAVE*WAVE) ! REFRM1 = 1.0E-6*(64.328+29498.1/(146.-WSQI)+255.4/(41.-WSQI)) ! RAYLAY = 5.40E-21*(REFRM1*WSQI)**2 !----------------------------------------------------------------------- ! ! DTAUX Local optical depth of each CTM level ! PIRAY Contribution of Rayleigh scattering to extinction ! PIAER Contribution of Aerosol scattering to extinction ! TTAU Optical depth of air vertically above each point (to top of atm) ! FTAU Attenuation of solar beam ! POMEGA Scattering phase function ! FMEAN Mean actinic flux at desired levels ! !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data USE module_peg_util, only: peg_message, peg_error_fatal !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year integer nspint ! Num of spectral intervals across solar parameter ( nspint = 4 ) ! spectrum for FAST-J real, dimension (nspint),save :: wavmid !cm real, dimension (nspint, kmaxd+1) :: sizeaer,extaer,waer,gaer,tauaer real, dimension (nspint, kmaxd+1) :: l2,l3,l4,l5,l6,l7 data wavmid & / 0.30e-4, 0.40e-4, 0.60e-4 ,0.999e-04 / INTEGER NL, N__, M__ PARAMETER (NL=500, N__=2*NL, M__=4) !wig increased nl from 350 to 500, 31-Oct-2005 REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1 REAL(kind=double), dimension(M__,2*M__) :: PM REAL(kind=double), dimension(2*M__) :: PM0 REAL(kind=double), dimension(2*M__,N__) :: POMEGA REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ REAL(kind=double), dimension(M__,M__,N__) :: DD REAL(kind=double), dimension(M__,N__) :: RR REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0 INTEGER ND,N,M,MFIT !jdf integer jndlev(lpar),jaddlv(nc),jaddto(nc+1) integer KW,km,i,j,k,l,ix,j1 integer isvode,jsvode real*8 QXMIE(MX),XLAER(MX),SSALB(MX) real*8 xlo2,xlo3,xlray,xltau,zk,taudn,tauup,zk2 real*8 WAVEL,XQO2(NB),XQO3(NB),FMEAN(lpar),POMEGAJ(2*M__,NC+1) real*8 DTAUX(NB),PIRAY(NB),PIAER(MX,NB),TTAU(NC+1),FTAU(NC+1) real*8 ftaulog,dttau,dpomega(2*M__) real*8 ftaulog2,dttau2,dpomega2(2*M__) ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ real*8 PIAER_MX1(NB) ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB character*80 msg ! !---Pick nearest Mie wavelength, no interpolation-------------- KM=1 if( WAVEL .gt. 355.d0 ) KM=2 if( WAVEL .gt. 500.d0 ) KM=3 ! if( WAVEL .gt. 800.d0 ) KM=4 !drop the 1000 nm wavelength ! !---For Mie code scale extinction at 1000 nm to wavelength WAVEL (QXMIE) ! define angstrom exponent ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ ang=log(QAA(1,MIEDX(1))/QAA(4,MIEDX(1)))/log(300./999.) do I=1,MX ! QAA is extinction efficiency QXMIE(I) = QAA(KM,MIEDX(I))/QAA(4,MIEDX(I)) ! scale to 550 nm using angstrom relationship ! note that this gives QXMIE at 550.0 nm = 1.0, aerosol optical thickness ! is defined at 550 nm ! convention -- I = 1 is aerosol, I > 1 are clouds if(I.eq.1) QXMIE(I) = (WAVEL/550.0)**ang SSALB(I) = SSA(KM,MIEDX(I)) ! single scattering albedo enddo ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB ! !---Reinitialize arrays do j=1,nc+1 ttau(j)=0.d0 ftau(j)=0.d0 enddo ! !---Set up total optical depth over each CTM level, DTAUX J1 = NLBATM do J=J1,NB XLO3=DO3(J)*XQO3(J) XLO2=DM(J)*XQO2(J)*0.20948d0 XLRAY=DM(J)*QRAYL(KW) ! Zero absorption for testing purposes ! call NOABS(XLO3,XLO2,XLRAY,AER(1,j),RFLECT) do I=1,MX ! I is aerosol type, j is level, AER(I,J)*QXMIE(I) is the layer aerosol optical thickness ! at 1000 nm , AER(I,J), times extinction efficiency for the layer (normalized to be one at 1000 nm) ! therefore xlaer(i) is the layer optical depth at the wavelength index KM XLAER(I)=AER(I,J)*QXMIE(I) enddo ! Total optical depth from all elements DTAUX(J)=XLO3+XLO2+XLRAY do I=1,MX DTAUX(J)=DTAUX(J)+XLAER(I) ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf xlaer',& ! i,j,km,dtaux(j),xlaer(i),aer(i,j),qxmie(i),xlo3,xlo2,xlray enddo ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ ! add in new aerosol information from Mie code ! layer aerosol optical thickness at wavelength index KM, layer j ! tauaer(km,j)=0.0 dtaux(j)=dtaux(j)+tauaer(km,j) ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf dtaux',& ! j,km,dtaux(j),tauaer(km,j) ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB ! Fractional extinction for Rayleigh scattering and each aerosol type PIRAY(J)=XLRAY/DTAUX(J) do I=1,MX PIAER(I,J)=SSALB(I)*XLAER(I)/DTAUX(J) enddo ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ ! note the level is now important PIAER_MX1(J)=waer(km,j)*tauaer(km,j)/DTAUX(J) ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB enddo ! of the level "j" loop ! !---Define the scattering phase fn. with mix of Rayleigh(1) & Mie(MIEDX) ! No. of quadrature pts fixed at 4 (M__), expansion of phase fn @ 8 N = M__ MFIT = 2*M__ do j=j1,NB ! jcb: layer index do i=1,MFIT pomegaj(i,j) = PIRAY(J)*PAA(I,KM,1) ! jcb: paa are the expansion coefficients of the phase function do k=1,MX ! jcb: mx is # of aerosols pomegaj(i,j) = pomegaj(i,j) + PIAER(K,J)*PAA(I,KM,MIEDX(K)) enddo enddo ! JJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJJ ! i is the # of coefficients, KM is the wavelength index, j is the level ! note the level in now relevant because we allow aerosol properties to ! vary by level pomegaj(1,j) = pomegaj(1,j) + PIAER_MX1(J)*1.0 ! 1.0 is l0 pomegaj(2,j) = pomegaj(2,j) + PIAER_MX1(J)*gaer(KM,j)*3.0 ! the three converts gear to l1 pomegaj(3,j) = pomegaj(3,j) + PIAER_MX1(J)*l2(KM,j) pomegaj(4,j) = pomegaj(4,j) + PIAER_MX1(J)*l3(KM,j) pomegaj(5,j) = pomegaj(5,j) + PIAER_MX1(J)*l4(KM,j) pomegaj(6,j) = pomegaj(6,j) + PIAER_MX1(J)*l5(KM,j) pomegaj(7,j) = pomegaj(7,j) + PIAER_MX1(J)*l6(KM,j) pomegaj(8,j) = pomegaj(7,j) + PIAER_MX1(J)*l7(KM,j) ! BBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBBB enddo ! !---Calculate attenuated incident beam EXP(-TTAU/U0) and flux on surface do J=J1,NB if(AMF(J,J).gt.0.0D0) then XLTAU=0.0D0 do I=1,NB XLTAU=XLTAU + DTAUX(I)*AMF(I,J) enddo if(XLTAU.gt.450.d0) then ! for compilers with no underflow trapping FTAU(j)=0.d0 else FTAU(J)=DEXP(-XLTAU) endif else FTAU(J)=0.0D0 endif enddo if(U0.gt.0.D0) then ZFLUX = U0*FTAU(J1)*RFLECT/(1.d0+RFLECT) else ZFLUX = 0.d0 endif ! !------------------------------------------------------------------------ ! Take optical properties on CTM layers and convert to a photolysis ! level grid corresponding to layer centres and boundaries. This is ! required so that J-values can be calculated for the centre of CTM ! layers; the index of these layers is kept in the jndlev array. !------------------------------------------------------------------------ ! ! Set lower boundary and levels to calculate J-values at J1=2*J1-1 do j=1,lpar jndlev(j)=2*j enddo ! ! Calculate column optical depths above each level, TTAU TTAU(NC+1)=0.0D0 do J=NC,J1,-1 I=(J+1)/2 TTAU(J)=TTAU(J+1) + 0.5d0*DTAUX(I) jaddlv(j)=int(0.5d0*DTAUX(I)/dtaumax) ! Subdivide cloud-top levels if required if(jadsub(j).gt.0) then jadsub(j)=min(jaddlv(j)+1,nint(dtausub))*(nint(dsubdiv)-1) jaddlv(j)=jaddlv(j)+jadsub(j) endif ! if(isvode.eq.8.and.jsvode.eq.1) print*,'jdf in opmie',& ! j,jadsub(j),jaddlv(j),ttau(j),dtaux(i) enddo ! ! Calculate attenuated beam, FTAU, level boundaries then level centres FTAU(NC+1)=1.0D0 do J=NC-1,J1,-2 I=(J+1)/2 FTAU(J)=FTAU(I) enddo do J=NC,J1,-2 FTAU(J)=sqrt(FTAU(J+1)*FTAU(J-1)) enddo ! ! Calculate scattering properties, level centres then level boundaries ! using an inverse interpolation to give correctly-weighted values do j=NC,J1,-2 do i=1,MFIT pomegaj(i,j) = pomegaj(i,j/2) enddo enddo do j=J1+2,nc,2 taudn = ttau(j-1)-ttau(j) tauup = ttau(j)-ttau(j+1) do i=1,MFIT pomegaj(i,j) = (pomegaj(i,j-1)*taudn + & pomegaj(i,j+1)*tauup) / (taudn+tauup) enddo enddo ! Define lower and upper boundaries do i=1,MFIT pomegaj(i,J1) = pomegaj(i,J1+1) pomegaj(i,nc+1) = pomegaj(i,nc) enddo ! !------------------------------------------------------------------------ ! Calculate cumulative total and define levels we want J-values at. ! Sum upwards for levels, and then downwards for Mie code readjustments. ! ! jaddlv(i) Number of new levels to add between (i) and (i+1) ! jaddto(i) Total number of new levels to add to and above level (i) ! jndlev(j) Level needed for J-value for CTM layer (j) ! !------------------------------------------------------------------------ ! ! Reinitialize level arrays do j=1,nc+1 jaddto(j)=0 enddo ! jaddto(J1)=jaddlv(J1) do j=J1+1,nc jaddto(j)=jaddto(j-1)+jaddlv(j) enddo if((jaddto(nc)+nc).gt.nl) then ! print*,'jdf mie',isvode,jsvode,jaddto(nc),nc,nl write ( msg, '(a, 2i6)' ) & 'FASTJ Max NL exceeded ' // & 'jaddto(nc)+nc NL', jaddto(nc)+nc,NL call peg_message( lunerr, msg ) msg = 'FASTJ subr OPMIE error. Max NL exceeded' call peg_error_fatal( lunerr, msg ) ! write(6,1500) jaddto(nc)+nc, 'NL',NL ! stop endif do i=1,lpar jndlev(i)=jndlev(i)+jaddto(jndlev(i)-1) enddo jaddto(nc)=jaddlv(nc) do j=nc-1,J1,-1 jaddto(j)=jaddto(j+1)+jaddlv(j) enddo ! !---------------------SET UP FOR MIE CODE------------------------------- ! ! Transpose the ascending TTAU grid to a descending ZTAU grid. ! Double the resolution - TTAU points become the odd points on the ! ZTAU grid, even points needed for asymm phase fn soln, contain 'h'. ! Odd point added at top of grid for unattenuated beam (Z='inf') ! ! Surface: TTAU(1) now use ZTAU(2*NC+1) ! Top: TTAU(NC) now use ZTAU(3) ! Infinity: now use ZTAU(1) ! ! Mie scattering code only used from surface to level NC !------------------------------------------------------------------------ ! ! Initialise all Fast-J optical property arrays do k=1,N__ do i=1,MFIT pomega(i,k) = 0.d0 enddo ztau(k) = 0.d0 fz(k) = 0.d0 enddo ! ! Ascend through atmosphere transposing grid and adding extra points do j=J1,nc+1 k = 2*(nc+1-j)+2*jaddto(j)+1 ztau(k)= ttau(j) fz(k) = ftau(j) do i=1,MFIT pomega(i,k) = pomegaj(i,j) enddo enddo ! ! Check profiles if desired ! ND = 2*(NC+jaddto(J1)-J1) + 3 ! if(kw.eq.1) call CH_PROF ! !------------------------------------------------------------------------ ! Insert new levels, working downwards from the top of the atmosphere ! to the surface (down in 'j', up in 'k'). This allows ztau and pomega ! to be incremented linearly (in a +ve sense), and the flux fz to be ! attenuated top-down (avoiding problems where lower level fluxes are ! zero). ! ! zk fractional increment in level ! dttau change in ttau per increment (linear, positive) ! dpomega change in pomega per increment (linear) ! ftaulog change in ftau per increment (exponential, normally < 1) ! !------------------------------------------------------------------------ ! do j=nc,J1,-1 zk = 0.5d0/(1.d0+dble(jaddlv(j)-jadsub(j))) dttau = (ttau(j)-ttau(j+1))*zk do i=1,MFIT dpomega(i) = (pomegaj(i,j)-pomegaj(i,j+1))*zk enddo ! Filter attenuation factor - set minimum at 1.0d-05 if(ftau(j+1).eq.0.d0) then ftaulog=0.d0 else ftaulog = ftau(j)/ftau(j+1) if(ftaulog.lt.1.d-150) then ftaulog=1.0d-05 else ftaulog=exp(log(ftaulog)*zk) endif endif k = 2*(nc-j+jaddto(j)-jaddlv(j))+1 ! k at level j+1 l = 0 ! Additional subdivision of first level if required if(jadsub(j).ne.0) then l=jadsub(j)/nint(dsubdiv-1) zk2=1.d0/dsubdiv dttau2=dttau*zk2 ftaulog2=ftaulog**zk2 do i=1,MFIT dpomega2(i)=dpomega(i)*zk2 enddo do ix=1,2*(jadsub(j)+l) ztau(k+1) = ztau(k) + dttau2 fz(k+1) = fz(k)*ftaulog2 do i=1,MFIT pomega(i,k+1) = pomega(i,k) + dpomega2(i) enddo k = k+1 enddo endif l = 2*(jaddlv(j)-jadsub(j)-l)+1 ! ! Add values at all intermediate levels do ix=1,l ztau(k+1) = ztau(k) + dttau fz(k+1) = fz(k)*ftaulog do i=1,MFIT pomega(i,k+1) = pomega(i,k) + dpomega(i) enddo k = k+1 enddo ! ! Alternate method to attenuate fluxes, fz, using 2nd-order finite ! difference scheme - just need to comment in section below ! ix = 2*(jaddlv(j)-jadsub(j))+1 ! if(l.le.0) then ! l=k-ix-1 ! else ! l=k-ix ! endif ! call efold(ftau(j+1),ftau(j),ix+1,fz(l)) ! if(jadsub(j).ne.0) then ! k = 2*(nc-j+jaddto(j)-jaddlv(j))+1 ! k at level j+1 ! ix=2*(jadsub(j)+(jadsub(j)/nint(dsubdiv-1))) ! call efold(ftau(j+1),fz(k+ix),ix,fz(k)) ! endif ! enddo ! !---Update total number of levels and check doesn't exceed N__ ND = 2*(NC+jaddto(J1)-J1) + 3 if(nd.gt.N__) then write ( msg, '(a, 2i6)' ) & 'FASTJ Max N__ exceeded ' // & 'ND N__', ND, N__ call peg_message( lunerr, msg ) msg = 'FASTJ subr OPMIE error. Max N__ exceeded' call peg_error_fatal( lunerr, msg ) ! write(6,1500) ND, 'N__',N__ ! stop endif ! !---Add boundary/ground layer to ensure no negative J's caused by !---too large a TTAU-step in the 2nd-order lower b.c. ZTAU(ND+1) = ZTAU(ND)*1.000005d0 ZTAU(ND+2) = ZTAU(ND)*1.000010d0 zk=max(abs(U0),0.01d0) zk=dexp(-ZTAU(ND)*5.d-6/zk) FZ(ND+1) = FZ(ND)*zk FZ(ND+2) = FZ(ND+1)*zk do I=1,MFIT POMEGA(I,ND+1) = POMEGA(I,ND) POMEGA(I,ND+2) = POMEGA(I,ND) enddo ND = ND+2 ! ZU0 = U0 ZREFL = RFLECT ! !----------------------------------------- CALL MIESCT(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, & ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) !----------------------------------------- ! Accumulate attenuation for selected levels l=2*(NC+jaddto(J1))+3 do j=1,lpar k=l-(2*jndlev(j)) if(k.gt.ND-2) then FMEAN(j) = 0.d0 else FMEAN(j) = FJ(k) endif enddo ! return 1000 format(1x,i3,3(2x,1pe10.4),1x,i3) 1300 format(1x,50(i3)) 1500 format(' Too many levels in photolysis code: need ',i3,' but ',a, & ' dimensioned as ',i3) END SUBROUTINE OPMIE !********************************************************************* subroutine EFOLD (F0, F1, N, F) !----------------------------------------------------------------------- !--- calculate the e-fold between two boundaries, given the value !--- at both boundaries F0(x=0) = top, F1(x=1) = bottom. !--- presume that F(x) proportional to exp[-A*x] for x=0 to x=1 !--- d2F/dx2 = A*A*F and thus expect F1 = F0 * exp[-A] !--- alternatively, could define A = ln[F0/F1] !--- let X = A*x, d2F/dX2 = F !--- assume equal spacing (not necessary, but makes this easier) !--- with N-1 intermediate points (and N layers of thickness dX = A/N) !--- !--- 2nd-order finite difference: (F(i-1) - 2F(i) + F(i+1)) / dX*dX = F(i) !--- let D = 1 / dX*dX: ! ! 1 | 1 0 0 0 0 0 | | F0 | ! | | | 0 | ! 2 | -D 2D+1 -D 0 0 0 | | 0 | ! | | | 0 | ! 3 | 0 -D 2D+1 -D 0 0 | | 0 | ! | | | 0 | ! | 0 0 -D 2D+1 -D 0 | | 0 | ! | | | 0 | ! N | 0 0 0 -D 2D+1 -D | | 0 | ! | | | 0 | ! N+1 | 0 0 0 0 0 1 | | F1 | ! !----------------------------------------------------------------------- ! Advantage of scheme over simple attenuation factor: conserves total ! number of photons - very useful when using scheme for heating rates. ! Disadvantage: although reproduces e-folds very well for small flux ! differences, starts to drift off when many orders of magnitude are ! involved. !----------------------------------------------------------------------- implicit none real*8 F0,F1,F(250) !F(N+1) integer N integer I real*8 A,DX,D,DSQ,DDP1, B(101),R(101) ! if(F0.eq.0.d0) then do I=1,N F(I)=0.d0 enddo return elseif(F1.eq.0.d0) then A = DLOG(F0/1.d-250) else A = DLOG(F0/F1) endif ! DX = float(N)/A D = DX*DX DSQ = D*D DDP1 = D+D+1.d0 ! B(2) = DDP1 R(2) = +D*F0 do I=3,N B(I) = DDP1 - DSQ/B(I-1) R(I) = +D*R(I-1)/B(I-1) enddo F(N+1) = F1 do I=N,2,-1 F(I) = (R(I) + D*F(I+1))/B(I) enddo F(1) = F0 return end subroutine EFOLD subroutine CH_PROF(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, & ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) !----------------------------------------------------------------------- ! Check profiles to be passed to MIESCT !----------------------------------------------------------------------- USE module_peg_util, only: peg_message implicit none !jdf integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 INTEGER NL, N__, M__ PARAMETER (NL=350, N__=2*NL, M__=4) REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1 REAL(kind=double), dimension(M__,2*M__) :: PM REAL(kind=double), dimension(2*M__) :: PM0 REAL(kind=double), dimension(2*M__,N__) :: POMEGA REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ REAL(kind=double), dimension(M__,M__,N__) :: DD REAL(kind=double), dimension(M__,N__) :: RR REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0 INTEGER ND,N,M,MFIT !jdf integer i,j character*80 msg ! write(6,1100) 'lev','ztau','fz ','pomega( )' do i=1,ND if(ztau(i).ne.0.d0) then write ( msg, '(a, i3, 2(1x,1pe9.3))' ) & 'FASTJ subr CH_PROF ztau ne 0. check pomega. ' // & 'k ztau(k) fz(k) ', i,ztau(i),fz(i) call peg_message( lunerr, msg ) ! write(6,1200) i,ztau(i),fz(i),(pomega(j,i),j=1,8) endif enddo return 1100 format(1x,a3,4(a9,2x)) 1200 format(1x,i3,11(1x,1pe9.3)) end subroutine CH_PROF SUBROUTINE MIESCT(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, & ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) ! SUBROUTINE MIESCT !----------------------------------------------------------------------- ! This is an adaption of the Prather radiative transfer code, (mjp, 10/95) ! Prather, 1974, Astrophys. J. 192, 787-792. ! Sol'n of inhomogeneous Rayleigh scattering atmosphere. ! (original Rayleigh w/ polarization) ! Cochran and Trafton, 1978, Ap.J., 219, 756-762. ! Raman scattering in the atmospheres of the major planets. ! (first use of anisotropic code) ! Jacob, Gottlieb and Prather, 1989, J.Geophys.Res., 94, 12975-13002. ! Chemistry of a polluted cloudy boundary layer, ! (documentation of extension to anisotropic scattering) ! ! takes atmospheric structure and source terms from std J-code ! ALSO limited to 4 Gauss points, only calculates mean field! ! ! mean rad. field ONLY (M=1) ! initialize variables FIXED/UNUSED in this special version: ! FTOP = 1.0 = astrophysical flux (unit of pi) at SZA, -ZU0, use for scaling ! FBOT = 0.0 = external isotropic flux on lower boundary ! SISOTP = 0.0 = Specific Intensity of isotropic radiation incident from top ! ! SUBROUTINES: MIESCT needs 'jv_mie.cmn' ! BLKSLV needs 'jv_mie.cmn' ! GEN (ID) needs 'jv_mie.cmn' ! LEGND0 (X,PL,N) ! MATIN4 (A) ! GAUSSP (N,XPT,XWT) !----------------------------------------------------------------------- ! INCLUDE 'jv_mie.h' implicit none !jdf integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 INTEGER NL, N__, M__ PARAMETER (NL=350, N__=2*NL, M__=4) REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1 REAL(kind=double), dimension(M__,2*M__) :: PM REAL(kind=double), dimension(2*M__) :: PM0 REAL(kind=double), dimension(2*M__,N__) :: POMEGA REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ REAL(kind=double), dimension(M__,M__,N__) :: DD REAL(kind=double), dimension(M__,N__) :: RR REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0 INTEGER ND,N,M,MFIT !jdf integer i, id, im real*8 cmeq1 !----------------------------------------------------------------------- !---fix scattering to 4 Gauss pts = 8-stream CALL GAUSSP (N,EMU,WT) !---solve eqn of R.T. only for first-order M=1 ! ZFLUX = (ZU0*FZ(ND)*ZREFL+FBOT)/(1.0d0+ZREFL) ZFLUX = (ZU0*FZ(ND)*ZREFL)/(1.0d0+ZREFL) M=1 DO I=1,N CALL LEGND0 (EMU(I),PM0,MFIT) DO IM=M,MFIT PM(I,IM) = PM0(IM) ENDDO ENDDO ! CMEQ1 = 0.25D0 CALL LEGND0 (-ZU0,PM0,MFIT) DO IM=M,MFIT PM0(IM) = CMEQ1*PM0(IM) ENDDO ! ! CALL BLKSLV CALL BLKSLV(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, & ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) ! DO ID=1,ND,2 FJ(ID) = 4.0d0*FJ(ID) + FZ(ID) ENDDO RETURN END SUBROUTINE MIESCT SUBROUTINE BLKSLV(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, & ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) !----------------------------------------------------------------------- ! Solves the block tri-diagonal system: ! A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = H(I) !----------------------------------------------------------------------- ! INCLUDE 'jv_mie.h' implicit none !jdf integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 INTEGER NL, N__, M__ PARAMETER (NL=350, N__=2*NL, M__=4) REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1 REAL(kind=double), dimension(M__,2*M__) :: PM REAL(kind=double), dimension(2*M__) :: PM0 REAL(kind=double), dimension(2*M__,N__) :: POMEGA REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ REAL(kind=double), dimension(M__,M__,N__) :: DD REAL(kind=double), dimension(M__,N__) :: RR REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0 INTEGER ND,N,M,MFIT !jdf integer i, j, k, id real*8 thesum !-----------UPPER BOUNDARY ID=1 CALL GEN(1,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, & ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) CALL MATIN4 (B) DO I=1,N RR(I,1) = 0.0d0 DO J=1,N THESUM = 0.0d0 DO K=1,N THESUM = THESUM - B(I,K)*CC(K,J) ENDDO DD(I,J,1) = THESUM RR(I,1) = RR(I,1) + B(I,J)*H(J) ENDDO ENDDO !----------CONTINUE THROUGH ALL DEPTH POINTS ID=2 TO ID=ND-1 DO ID=2,ND-1 CALL GEN(ID,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, & ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) DO I=1,N DO J=1,N B(I,J) = B(I,J) + A(I)*DD(I,J,ID-1) ENDDO H(I) = H(I) - A(I)*RR(I,ID-1) ENDDO CALL MATIN4 (B) DO I=1,N RR(I,ID) = 0.0d0 DO J=1,N RR(I,ID) = RR(I,ID) + B(I,J)*H(J) DD(I,J,ID) = - B(I,J)*C1(J) ENDDO ENDDO ENDDO !---------FINAL DEPTH POINT: ND CALL GEN(ND,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU,ZFLUX, & ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) DO I=1,N DO J=1,N THESUM = 0.0d0 DO K=1,N THESUM = THESUM + AA(I,K)*DD(K,J,ND-1) ENDDO B(I,J) = B(I,J) + THESUM H(I) = H(I) - AA(I,J)*RR(J,ND-1) ENDDO ENDDO CALL MATIN4 (B) DO I=1,N RR(I,ND) = 0.0d0 DO J=1,N RR(I,ND) = RR(I,ND) + B(I,J)*H(J) ENDDO ENDDO !-----------BACK SOLUTION DO ID=ND-1,1,-1 DO I=1,N DO J=1,N RR(I,ID) = RR(I,ID) + DD(I,J,ID)*RR(J,ID+1) ENDDO ENDDO ENDDO !----------MEAN J & H DO ID=1,ND,2 FJ(ID) = 0.0d0 DO I=1,N FJ(ID) = FJ(ID) + RR(I,ID)*WT(I) ENDDO ENDDO DO ID=2,ND,2 FJ(ID) = 0.0d0 DO I=1,N FJ(ID) = FJ(ID) + RR(I,ID)*WT(I)*EMU(I) ENDDO ENDDO ! Output fluxes for testing purposes ! CALL CH_FLUX ! CALL CH_FLUX(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, & ! ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) ! RETURN END SUBROUTINE BLKSLV SUBROUTINE CH_FLUX(ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, & ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) !----------------------------------------------------------------------- ! Diagnostic routine to check fluxes at each level - makes most sense ! when running a conservative atmosphere (zero out absorption in ! OPMIE by calling the NOABS routine below) !----------------------------------------------------------------------- ! INCLUDE 'jv_mie.h' IMPLICIT NONE !jdf integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 INTEGER NL, N__, M__ PARAMETER (NL=350, N__=2*NL, M__=4) REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1 REAL(kind=double), dimension(M__,2*M__) :: PM REAL(kind=double), dimension(2*M__) :: PM0 REAL(kind=double), dimension(2*M__,N__) :: POMEGA REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ REAL(kind=double), dimension(M__,M__,N__) :: DD REAL(kind=double), dimension(M__,N__) :: RR REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0 INTEGER ND,N,M,MFIT !jdf integer I,ID real*8 FJCHEK(N__),FZMEAN ! ! Odd (h) levels held as actinic flux, so recalculate irradiances DO ID=1,ND,2 FJCHEK(ID) = 0.0d0 DO I=1,N FJCHEK(ID) = FJCHEK(ID) + RR(I,ID)*WT(I)*EMU(i) ENDDO ENDDO ! ! Even (j) levels are already held as irradiances DO ID=2,ND,2 DO I=1,N FJCHEK(ID) = FJ(ID) ENDDO ENDDO ! ! Output Downward and Upward fluxes down through atmosphere ! WRITE(6,1200) WRITE(34,1200) DO ID=2,ND,2 FZMEAN=sqrt(FZ(ID)*FZ(ID-1)) ! WRITE(6,1000) ID, ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)), & WRITE(34,1000) ID, ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1)), & 2.0*(FJCHEK(id)+FJCHEK(id-1)), & 2.0*(FJCHEK(id)+FJCHEK(id-1))/ & (ZU0*FZMEAN-2.0*(FJCHEK(id)-FJCHEK(id-1))) ENDDO RETURN 1000 FORMAT(1x,i3,1p,2E12.4,1x,0p,f9.4) 1200 FORMAT(1x,'Lev',3x,'Downward',4x,'Upward',7x,'Ratio') END SUBROUTINE CH_FLUX SUBROUTINE NOABS(XLO3,XLO2,XLRAY,BCAER,RFLECT) !----------------------------------------------------------------------- ! Zero out absorption terms to check scattering code. Leave a little ! Rayleigh to provide a minimal optical depth, and set surface albedo ! to unity. !----------------------------------------------------------------------- IMPLICIT NONE real*8 XLO3,XLO2,XLRAY,BCAER,RFLECT XLO3=0.d0 XLO2=0.d0 XLRAY=XLRAY*1.d-10 BCAER=0.d0 RFLECT=1.d0 RETURN END SUBROUTINE NOABS SUBROUTINE GEN(ID,ND,N,M,MFIT,POMEGA,PM,PM0,FZ,WT,EMU,ZTAU, & ZFLUX,ZREFL,FJ,A,C1,H,V1,B,AA,CC,S,W,U1,DD,RR,RADIUS,ZU0) !----------------------------------------------------------------------- ! Generates coefficient matrices for the block tri-diagonal system: ! A(I)*X(I-1) + B(I)*X(I) + C(I)*X(I+1) = H(I) !----------------------------------------------------------------------- ! INCLUDE 'jv_mie.h' IMPLICIT NONE !jdf integer, parameter :: single = 4 !compiler dependent value real*4 integer, parameter :: double = 8 !compiler dependent value real*8 INTEGER NL, N__, M__ PARAMETER (NL=350, N__=2*NL, M__=4) REAL(kind=double), dimension(M__) :: A,C1,H,V1,WT,EMU REAL(kind=double), dimension(M__,M__) :: B,AA,CC,S,W,U1 REAL(kind=double), dimension(M__,2*M__) :: PM REAL(kind=double), dimension(2*M__) :: PM0 REAL(kind=double), dimension(2*M__,N__) :: POMEGA REAL(kind=double), dimension(N__) :: ZTAU, FZ, FJ REAL(kind=double), dimension(M__,M__,N__) :: DD REAL(kind=double), dimension(M__,N__) :: RR REAL(kind=double) :: ZREFL,ZFLUX,RADIUS,ZU0 INTEGER ND,N,M,MFIT !jdf integer id, id0, id1, im, i, j, k, mstart real*8 sum0, sum1, sum2, sum3 real*8 deltau, d1, d2, surfac !--------------------------------------------- IF(ID.EQ.1 .OR. ID.EQ.ND) THEN !---------calculate generic 2nd-order terms for boundaries ID0 = ID ID1 = ID+1 IF(ID.GE.ND) ID1 = ID-1 DO 10 I=1,N SUM0 = 0.0d0 SUM1 = 0.0d0 SUM2 = 0.0d0 SUM3 = 0.0d0 DO IM=M,MFIT,2 SUM0 = SUM0 + POMEGA(IM,ID0)*PM(I,IM)*PM0(IM) SUM2 = SUM2 + POMEGA(IM,ID1)*PM(I,IM)*PM0(IM) ENDDO DO IM=M+1,MFIT,2 SUM1 = SUM1 + POMEGA(IM,ID0)*PM(I,IM)*PM0(IM) SUM3 = SUM3 + POMEGA(IM,ID1)*PM(I,IM)*PM0(IM) ENDDO H(I) = 0.5d0*(SUM0*FZ(ID0) + SUM2*FZ(ID1)) A(I) = 0.5d0*(SUM1*FZ(ID0) + SUM3*FZ(ID1)) DO J=1,I SUM0 = 0.0d0 SUM1 = 0.0d0 SUM2 = 0.0d0 SUM3 = 0.0d0 DO IM=M,MFIT,2 SUM0 = SUM0 + POMEGA(IM,ID0)*PM(I,IM)*PM(J,IM) SUM2 = SUM2 + POMEGA(IM,ID1)*PM(I,IM)*PM(J,IM) ENDDO DO IM=M+1,MFIT,2 SUM1 = SUM1 + POMEGA(IM,ID0)*PM(I,IM)*PM(J,IM) SUM3 = SUM3 + POMEGA(IM,ID1)*PM(I,IM)*PM(J,IM) ENDDO S(I,J) = - SUM2*WT(J) S(J,I) = - SUM2*WT(I) W(I,J) = - SUM1*WT(J) W(J,I) = - SUM1*WT(I) U1(I,J) = - SUM3*WT(J) U1(J,I) = - SUM3*WT(I) SUM0 = 0.5d0*(SUM0 + SUM2) B(I,J) = - SUM0*WT(J) B(J,I) = - SUM0*WT(I) ENDDO S(I,I) = S(I,I) + 1.0d0 W(I,I) = W(I,I) + 1.0d0 U1(I,I) = U1(I,I) + 1.0d0 B(I,I) = B(I,I) + 1.0d0 10 CONTINUE DO I=1,N SUM0 = 0.0d0 DO J=1,N SUM0 = SUM0 + S(I,J)*A(J)/EMU(J) ENDDO C1(I) = SUM0 ENDDO DO I=1,N DO J=1,N SUM0 = 0.0d0 SUM2 = 0.0d0 DO K=1,N SUM0 = SUM0 + S(J,K)*W(K,I)/EMU(K) SUM2 = SUM2 + S(J,K)*U1(K,I)/EMU(K) ENDDO A(J) = SUM0 V1(J) = SUM2 ENDDO DO J=1,N W(J,I) = A(J) U1(J,I) = V1(J) ENDDO ENDDO IF (ID.EQ.1) THEN !-------------upper boundary, 2nd-order, C-matrix is full (CC) DELTAU = ZTAU(2) - ZTAU(1) D2 = 0.25d0*DELTAU DO I=1,N D1 = EMU(I)/DELTAU DO J=1,N B(I,J) = B(I,J) + D2*W(I,J) CC(I,J) = D2*U1(I,J) ENDDO B(I,I) = B(I,I) + D1 CC(I,I) = CC(I,I) - D1 ! H(I) = H(I) + 2.0d0*D2*C1(I) + D1*SISOTP H(I) = H(I) + 2.0d0*D2*C1(I) A(I) = 0.0d0 ENDDO ELSE !-------------lower boundary, 2nd-order, A-matrix is full (AA) DELTAU = ZTAU(ND) - ZTAU(ND-1) D2 = 0.25d0*DELTAU SURFAC = 4.0d0*ZREFL/(1.0d0 + ZREFL) DO I=1,N D1 = EMU(I)/DELTAU H(I) = H(I) - 2.0d0*D2*C1(I) SUM0 = 0.0d0 DO J=1,N SUM0 = SUM0 + W(I,J) ENDDO SUM0 = D1 + D2*SUM0 SUM1 = SURFAC*SUM0 DO J=1,N B(I,J) = B(I,J) + D2*W(I,J) - SUM1*EMU(J)*WT(J) ENDDO B(I,I) = B(I,I) + D1 H(I) = H(I) + SUM0*ZFLUX DO J=1,N AA(I,J) = - D2*U1(I,J) ENDDO AA(I,I) = AA(I,I) + D1 C1(I) = 0.0d0 ENDDO ENDIF !------------intermediate points: can be even or odd, A & C diagonal ELSE DELTAU = ZTAU(ID+1) - ZTAU(ID-1) MSTART = M + MOD(ID+1,2) DO I=1,N A(I) = EMU(I)/DELTAU C1(I) = -A(I) SUM0 = 0.0d0 DO IM=MSTART,MFIT,2 SUM0 = SUM0 + POMEGA(IM,ID)*PM(I,IM)*PM0(IM) ENDDO H(I) = SUM0*FZ(ID) DO J=1,I SUM0 = 0.0d0 DO IM=MSTART,MFIT,2 SUM0 = SUM0 + POMEGA(IM,ID)*PM(I,IM)*PM(J,IM) ENDDO B(I,J) = - SUM0*WT(J) B(J,I) = - SUM0*WT(I) ENDDO B(I,I) = B(I,I) + 1.0d0 ENDDO ENDIF RETURN END SUBROUTINE GEN SUBROUTINE LEGND0 (X,PL,N) !---Calculates ORDINARY LEGENDRE fns of X (real) !--- from P[0] = PL(1) = 1, P[1] = X, .... P[N-1] = PL(N) IMPLICIT NONE INTEGER N,I REAL*8 X,PL(N),DEN !---Always does PL(2) = P[1] PL(1) = 1.D0 PL(2) = X DO I=3,N DEN = (I-1) PL(I) = PL(I-1)*X*(2.d0-1.D0/DEN) - PL(I-2)*(1.d0-1.D0/DEN) ENDDO RETURN END SUBROUTINE LEGND0 SUBROUTINE MATIN4 (A) !----------------------------------------------------------------------- ! invert 4x4 matrix A(4,4) in place with L-U decomposition (mjp, old...) !----------------------------------------------------------------------- IMPLICIT NONE REAL*8 A(4,4) !---SETUP L AND U A(2,1) = A(2,1)/A(1,1) A(2,2) = A(2,2)-A(2,1)*A(1,2) A(2,3) = A(2,3)-A(2,1)*A(1,3) A(2,4) = A(2,4)-A(2,1)*A(1,4) A(3,1) = A(3,1)/A(1,1) A(3,2) = (A(3,2)-A(3,1)*A(1,2))/A(2,2) A(3,3) = A(3,3)-A(3,1)*A(1,3)-A(3,2)*A(2,3) A(3,4) = A(3,4)-A(3,1)*A(1,4)-A(3,2)*A(2,4) A(4,1) = A(4,1)/A(1,1) A(4,2) = (A(4,2)-A(4,1)*A(1,2))/A(2,2) A(4,3) = (A(4,3)-A(4,1)*A(1,3)-A(4,2)*A(2,3))/A(3,3) A(4,4) = A(4,4)-A(4,1)*A(1,4)-A(4,2)*A(2,4)-A(4,3)*A(3,4) !---INVERT L A(4,3) = -A(4,3) A(4,2) = -A(4,2)-A(4,3)*A(3,2) A(4,1) = -A(4,1)-A(4,2)*A(2,1)-A(4,3)*A(3,1) A(3,2) = -A(3,2) A(3,1) = -A(3,1)-A(3,2)*A(2,1) A(2,1) = -A(2,1) !---INVERT U A(4,4) = 1.D0/A(4,4) A(3,4) = -A(3,4)*A(4,4)/A(3,3) A(3,3) = 1.D0/A(3,3) A(2,4) = -(A(2,3)*A(3,4)+A(2,4)*A(4,4))/A(2,2) A(2,3) = -A(2,3)*A(3,3)/A(2,2) A(2,2) = 1.D0/A(2,2) A(1,4) = -(A(1,2)*A(2,4)+A(1,3)*A(3,4)+A(1,4)*A(4,4))/A(1,1) A(1,3) = -(A(1,2)*A(2,3)+A(1,3)*A(3,3))/A(1,1) A(1,2) = -A(1,2)*A(2,2)/A(1,1) A(1,1) = 1.D0/A(1,1) !---MULTIPLY (U-INVERSE)*(L-INVERSE) A(1,1) = A(1,1)+A(1,2)*A(2,1)+A(1,3)*A(3,1)+A(1,4)*A(4,1) A(1,2) = A(1,2)+A(1,3)*A(3,2)+A(1,4)*A(4,2) A(1,3) = A(1,3)+A(1,4)*A(4,3) A(2,1) = A(2,2)*A(2,1)+A(2,3)*A(3,1)+A(2,4)*A(4,1) A(2,2) = A(2,2)+A(2,3)*A(3,2)+A(2,4)*A(4,2) A(2,3) = A(2,3)+A(2,4)*A(4,3) A(3,1) = A(3,3)*A(3,1)+A(3,4)*A(4,1) A(3,2) = A(3,3)*A(3,2)+A(3,4)*A(4,2) A(3,3) = A(3,3)+A(3,4)*A(4,3) A(4,1) = A(4,4)*A(4,1) A(4,2) = A(4,4)*A(4,2) A(4,3) = A(4,4)*A(4,3) RETURN END SUBROUTINE MATIN4 SUBROUTINE GAUSSP (N,XPT,XWT) !----------------------------------------------------------------------- ! Loads in pre-set Gauss points for 4 angles from 0 to +1 in cos(theta)=mu !----------------------------------------------------------------------- IMPLICIT NONE INTEGER N,I REAL*8 XPT(N),XWT(N) REAL*8 GPT4(4),GWT4(4) DATA GPT4/.06943184420297D0,.33000947820757D0,.66999052179243D0, & .93056815579703D0/ DATA GWT4/.17392742256873D0,.32607257743127D0,.32607257743127D0, & .17392742256873D0/ N = 4 DO I=1,N XPT(I) = GPT4(I) XWT(I) = GWT4(I) ENDDO RETURN END SUBROUTINE GAUSSP ! subroutine aeroden(zz,v,aerd) ! purpose: find number density of boundary layer aerosols, aerd, ! at a given altitude, zz, and for a specified visibility ! input: ! zz altitude (km) ! v visibility for a horizontal surface path (km) ! output: ! aerd aerosol density at altitude z ! the vertical distribution of the boundary layer aerosol density is ! based on the 5s vertical profile models for 5 and 23 km visibility. ! above 5 km, the aden05 and aden23 models are the same ! below 5 km, the models differ as follows; ! aden05 0.99 km scale height (94% of extinction occurs below 5 km) ! aden23 1.45 km scale heigth (80% of extinction occurs below 5 km) ! implicit none integer mz, nz parameter (mz=33) real v,aerd real*8 zz ! compatability with fastj real alt, aden05, aden23, aer05,aer23 dimension alt(mz),aden05(mz),aden23(mz) !jdf dimension zbaer(*),dbaer(*) real z, f, wth integer k, kp save alt,aden05,aden23,nz data nz/mz/ data alt/ & 0.0, 1.0, 2.0, 3.0, 4.0, & 5.0, 6.0, 7.0, 8.0, 9.0, & 10.0, 11.0, 12.0, 13.0, 14.0, & 15.0, 16.0, 17.0, 18.0, 19.0, & 20.0, 21.0, 22.0, 23.0, 24.0, & 25.0, 30.0, 35.0, 40.0, 45.0, & 50.0, 70.0, 100.0/ data aden05/ & 1.378E+04, 5.030E+03, 1.844E+03, 6.731E+02, 2.453E+02, & 8.987E+01, 6.337E+01, 5.890E+01, 6.069E+01, 5.818E+01, & 5.675E+01, 5.317E+01, 5.585E+01, 5.156E+01, 5.048E+01, & 4.744E+01, 4.511E+01, 4.458E+01, 4.314E+01, 3.634E+01, & 2.667E+01, 1.933E+01, 1.455E+01, 1.113E+01, 8.826E+00, & 7.429E+00, 2.238E+00, 5.890E-01, 1.550E-01, 4.082E-02, & 1.078E-02, 5.550E-05, 1.969E-08/ data aden23/ & 2.828E+03, 1.244E+03, 5.371E+02, 2.256E+02, 1.192E+02, & 8.987E+01, 6.337E+01, 5.890E+01, 6.069E+01, 5.818E+01, & 5.675E+01, 5.317E+01, 5.585E+01, 5.156E+01, 5.048E+01, & 4.744E+01, 4.511E+01, 4.458E+01, 4.314E+01, 3.634E+01, & 2.667E+01, 1.933E+01, 1.455E+01, 1.113E+01, 8.826E+00, & 7.429E+00, 2.238E+00, 5.890E-01, 1.550E-01, 4.082E-02, & 1.078E-02, 5.550E-05, 1.969E-08/ ! z=max(0.,min(100.,real(zz))) aerd=0. if(z.gt.alt(nz)) return call locate(alt,nz,z,k) kp=k+1 f=(z-alt(k))/(alt(kp)-alt(k)) if(min(aden05(k),aden05(kp),aden23(k),aden23(kp)).le.0.) then aer05=aden05(k)*(1.-f)+aden05(kp)*f aer23=aden23(k)*(1.-f)+aden23(kp)*f else aer05=aden05(k)*(aden05(kp)/aden05(k))**f aer23=aden23(k)*(aden23(kp)/aden23(k))**f endif wth=(1./v-1/5.)/(1./23.-1./5.) wth=max(0.,min(1.,wth)) aerd=(1.-wth)*aer05+wth*aer23 ! write(*,*) 'aeroden k,kp,z,aer05(k),aer05(kp),f,aerd' ! write(*,'(2i5,1p5e11.3)') k,kp,z,aden05(k),aden05(kp),f,aerd return end subroutine aeroden !======================================================================= subroutine locate(xx,n,x,j) ! ! purpose: given an array xx of length n, and given a value X, returns ! a value J such that X is between xx(j) and xx(j+1). xx must ! be monotonic, either increasing of decreasing. this function ! returns j=1 or j=n-1 if x is out of range. !c ! input: ! xx monitonic table ! n size of xx ! x single floating point value perhaps within the range of xx ! ! output: ! function returns index value j, such that ! ! for an increasing table ! ! xx(j) .lt. x .le. xx(j+1), ! j=1 for x .lt. xx(1) ! j=n-1 for x .gt. xx(n) ! ! for a decreasing table ! xx(j) .le. x .lt. xx(j+1) ! j=n-1 for x .lt. xx(n) ! j=1 for x .gt. xx(1) ! implicit none integer j,n real x,xx(n) integer jl,jm,ju ! . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . if(x.eq.xx(1)) then j=1 return endif if(x.eq.xx(n)) then j=n-1 return endif jl=1 ju=n 10 if(ju-jl.gt.1) then jm=(ju+jl)/2 if((xx(n).gt.xx(1)).eqv.(x.gt.xx(jm)))then jl=jm else ju=jm endif goto 10 endif j=jl return end subroutine locate !************************************************************************ subroutine rd_tjpl2 !----------------------------------------------------------------------- ! set wavelength bins, solar fluxes, Rayleigh parameters, temperature- ! dependent cross sections and Rayleigh/aerosol scattering phase functions ! with temperature dependences. Current data originates from JPL 2000 !----------------------------------------------------------------------- ! ! NJVAL Number of species to calculate J-values for ! NWWW Number of wavelength bins, from NW1:NW2 ! WBIN Boundaries of wavelength bins ! WL Centres of wavelength bins - 'effective wavelength' ! FL Solar flux incident on top of atmosphere (cm-2.s-1) ! QRAYL Rayleigh parameters (effective cross-section) (cm2) ! QBC Black Carbon absorption extinct. (specific cross-sect.) (m2/g) ! QO2 O2 cross-sections ! QO3 O3 cross-sections ! Q1D O3 => O(1D) quantum yield ! TQQ Temperature for supplied cross sections ! QQQ Supplied cross sections in each wavelength bin (cm2) ! NAA Number of categories for scattering phase functions ! QAA Aerosol scattering phase functions ! NK Number of wavelengths at which functions supplied (set as 4) ! WAA Wavelengths for the NK supplied phase functions ! PAA Phase function: first 8 terms of expansion ! RAA Effective radius associated with aerosol type ! SSA Single scattering albedo ! ! npdep Number of pressure dependencies ! zpdep Pressure dependencies by wavelength bin ! jpdep Index of cross sections requiring pressure dependence ! lpdep Label for pressure dependence ! !----------------------------------------------------------------------- USE module_data_mosaic_other, only : kmaxd USE module_fastj_data USE module_peg_util, only: peg_message, peg_error_fatal IMPLICIT NONE !jdf ! Print Fast-J diagnostics if iprint /= 0 integer, parameter :: iprint = 0 integer, parameter :: single = 4 !compiler dependent value real*4 ! integer, parameter :: double = 8 !compiler dependent value real*8 integer,parameter :: ipar_fastj=1,jpar=1 ! integer,parameter :: jppj=14 !Number of photolytic reactions supplied logical,parameter :: ldeg45=.false. !Logical flag for degraded CTM resolution integer lpar !Number of levels in CTM integer jpnl !Number of levels requiring chemistry real(kind=double), dimension(ipar_fastj) :: xgrd ! Longitude (midpoint, radians) real(kind=double), dimension(jpar) :: ygrd ! Latitude (midpoint, radians) real(kind=double), dimension(jpar) :: ydgrd ! Latitude (midpoint, degrees) real(kind=double), dimension(kmaxd+1) :: etaa ! Eta(a) value for level boundaries real(kind=double), dimension(kmaxd+1) :: etab ! Eta(b) value for level boundaries real(kind=double) :: tau_fastj ! Time of Day (hours, GMT) integer month_fastj ! Number of month (1-12) integer iday_fastj ! Day of year !jdf integer i, j, k character*7 lpdep(3) character*80 msg if(NJVAL.gt.NS) then ! fastj input files are not set up for current situation write ( msg, '(a, 2i6)' ) & 'FASTJ # xsect supplied > max allowed ' // & 'NJVAL NS ', NJVAL, NS call peg_message( lunerr, msg ) msg = & 'FASTJ Setup Error: # xsect supplied > max allowed. Increase NS' call peg_error_fatal( lunerr, msg ) ! write(6,300) NJVAL,NS ! stop endif if(NAA.gt.NP) then write ( msg, '(a, 2i6)' ) & 'FASTJ # aerosol/cloud types > NP ' // & 'NAA NP ', NAA ,NP call peg_message( lunerr, msg ) msg = & 'FASTJ Setup Error: Too many phase functions supplied. Increase NP' call peg_error_fatal( lunerr, msg ) ! write(6,350) NAA ! stop endif !---Zero index arrays do j=1,jppj jind(j)=0 enddo do j=1,NJVAL jpdep(j)=0 enddo do j=1,nh hzind(j)=0 enddo ! !---Set mapping index do j=1,NJVAL do k=1,jppj if(jlabel(k).eq.titlej(1,j)) jind(k)=j ! write(6,*)k,jind(k) ! jcb ! write(6,*)jlabel(k),titlej(1,j) ! jcb enddo do k=1,npdep if(lpdep(k).eq.titlej(1,j)) jpdep(j)=k enddo enddo do k=1,jppj if(jfacta(k).eq.0.d0) then ! write(6,*) 'Not using photolysis reaction ',k write ( msg, '(a, i6)' ) & 'FASTJ Not using photolysis reaction ' , k call peg_message( lunerr, msg ) end if if(jind(k).eq.0) then if(jfacta(k).eq.0.d0) then jind(k)=1 else write ( msg, '(a, i6)' ) & 'FASTJ Which J-rate for photolysis reaction ' , k call peg_message( lunerr, msg ) ! write(6,*) 'Which J-rate for photolysis reaction ',k,' ?' ! stop msg = 'FASTJ subr rd_tjpl2 Unknown Jrate. Incorrect FASTJ setup' call peg_error_fatal( lunerr, msg ) endif endif enddo ! Herzberg index i=0 do j=1,nhz do k=1,jppj if(jlabel(k).eq.hzlab(j)) then i=i+1 hzind(i)=k hztoa(i)=hztmp(j)*jfacta(k) endif enddo enddo nhz=i if(nhz.eq.0) then if(iprint.ne.0) then write ( msg, '(a)' ) & 'FASTJ Not using Herzberg bin ' call peg_message( lunerr, msg ) ! write(6,400) end if else if(iprint.ne.0) then write ( msg, '(a)' ) & 'FASTJ Using Herzberg bin for: ' call peg_message( lunerr, msg ) write( msg, '(a,10a7)' ) & 'FASTJ ', (jlabel(hzind(i)),i=1,nhz) ! write(6,420) (jlabel(hzind(i)),i=1,nhz) end if endif ! 300 format(' Number of x-sections supplied to Fast-J: ',i3,/, & ! ' Maximum number allowed (NS) only set to: ',i3, & ! ' - increase in jv_cmn.h',/, & ! 'RESULTS WILL BE IN ERROR' ) ! 350 format(' Too many phase functions supplied; increase NP to ',i2, & ! /,'RESULTS WILL BE IN ERROR' ) ! 400 format(' Not using Herzberg bin') ! 420 format(' Using Herzberg bin for: ',10a7) return end subroutine rd_tjpl2 !******************************************************************** end module module_phot_fastj