!--------------------------------------------------------------------------- !BOP module fastJX65_mod ! ! !USES: !use GmiASCIIoperations_mod, only : AsciiOpenRead implicit none private ! ! !PUBLIC MEMBER FUNCTIONS: public :: controlFastJX65 public :: initializeFastJX65 public :: getQAA_RAAinFastJX65 ! INTEGER, PARAMETER :: numScat = 5 real*8, parameter :: PI=3.141592653589793d0 real*8, parameter :: Deg2Rad = PI / 180.0d0 real*8, parameter :: SecPerHour = 3600.0d0 !!! testing only real*8, allocatable, public :: ODdiag(:,:,:,:) real*8, allocatable, public :: OOJdiag(:,:,:), DDJdiag(:,:,:) ! # include "gmi_AerDust_const.h" # include "setkin_par.h" ! ! !DESCRIPTION: ! ! !HISTORY ! 2016.05.26 Luke Oman & Huisheng Bian ! Finished the FastJX 6.5 implementation ! !EOP !--------------------------------------------------------------------------- contains !--------------------------------------------------------------------------- !BOP ! ! !IROUTINE: controlFastJX65 ! ! !INTERFACE: ! subroutine controlFastJX65 (k1, k2, chem_mask_khi, num_qjs, & & month_gmi, jday, time_sec, SZA_ij, & & tau_clw_ij, tau_cli_ij, & & press3e_ij, pctm_ij, kel_ij, optdepth_ij, & & surf_alb_ij, qjgmi_ij, relHumidity_ij, & & overheadO3col_ij, ODAER_ij, ODMDUST_ij, ozone_ij) ! # include "parm_CTM_fastJX65.h" # include "cmn_JVdat_fastJX65.h" # include "cmn_metdat_fastJX65.h" # include "parm_MIE_fastJX65.h" ! ! !INPUT PARAMETERS: integer, intent(in) :: k1 ! first interface index integer, intent(in) :: k2 ! last interface index integer, intent(in) :: chem_mask_khi ! number of chemistry levels [jpnl] integer, intent(in) :: num_qjs ! number of photolysis reactions [jppj] integer, intent(in) :: month_gmi ! number of month (1- 12) [month] integer, intent(in) :: jday real*8, intent(in) :: kel_ij(k1:k2) ! temperature at box centers (degK) [t] real*8, intent(in) :: press3e_ij(k1:k2) real*8, intent(in) :: time_sec ! time of day in model (s, GMT) [gmtau (hrs)] real*8, intent(in) :: SZA_ij ! solar zenith angle real*8 , intent(in) :: tau_clw_ij(k1:k2) ! optical thickness for liquid cloud real*8 , intent(in) :: tau_cli_ij(k1:k2) ! optical thickness for ice cloud real*8, intent(in) :: pctm_ij ! surface pressure (mb) [p] real*8, intent(in) :: surf_alb_ij ! surface albedo (fraction 0-1) [sa] real*8, intent(in) :: relHumidity_ij(k1:k2) ! relative humidity (fraction 0-1) [relh] ! real*8, intent(in) :: optdepth_ij(k1:k2) ! optical depth in box (unitless) [od] real*8, intent(in) :: ODAER_ij(k1:k2,NSADaer*nrh_b) ! optical depth for aerosol real*8, intent(in) :: ODMDUST_ij(k1:k2,NSADdust) ! optical depth for mineral dust ! ! real*8, optional, intent(in) :: ozone_ij(k1:k2) ! mixing ratio for ozone real*8, intent(in), optional :: ozone_ij(k1:k2) ! mixing ratio for ozone ! ! !OUTPUT PARAMETERS: real*8, intent(out) :: overheadO3col_ij(k1:k2) real*8, intent(out) :: qjgmi_ij(k1:chem_mask_khi, num_qjs) ! ! !LOCAL VARIABLES: real*8 ZPJ(chem_mask_khi-k1+1, num_qjs) !2-D array of J's indexed to CTM chemistry! real*8 GMTAU,DELTAU, ALBEDO,SZA real*8 SOLF, U0, FREFL integer I,J,K,L,ILNG,JLAT,ODINDX, IDAY, il, ik, ic logical, save :: first = .true. logical :: rootProc !EOP !--------------------------------------------------------------------------- !BOC ILNG = I_ JLAT = J_ L_ = k2 - k1 + 1 JVL_ = chem_mask_khi JVN_ = num_qjs L1_ = L_ + 1 L2_ = 2 * L1_ MONTH = month_gmi IDAY = jday GMTAU = time_sec / SecPerHour !XDGRD(I_) = londeg_i !YDGRD(J_) = latdeg_j ! Set surface pressure ! PMEAN(I_, J_) = pctm_ij ! Set column temperature T(I_, J_, 1:L_) = kel_ij(k1:k2) ! Set column relative humidity RELH(1:L_) = relHumidity_ij(k1:k2) ! set cloud optical depth ! begBian, 11/21/2011 ! CLDLWP(I_, J_, 1:L_) = optdepth_ij(k1:k2) ! CLDLWP(I_, J_, L1_) = 0.0 ! ! set cloud index !! start-by-hbian ! do ik = k1, k2 ! if( tau_clw_ij(ik) .ge. tau_cli_ij(ik)) then ! CLDNDX(I_, J_, ik) = 9 ! else ! CLDNDX(I_, J_, ik) = 13 ! endif ! enddo ! please define these variables and input their values !CLDW(I_, J_, 1:L_) = 0.0d0 ! tau_clw_ij(k1:k2) !CLDW(I_, J_, L1_) = 0.d0 !CLDI(I_, J_, 1:L_) = 0.0d0 ! tau_cli_ij(k1:k2) !CLDI(I_, J_, L1_) = 0.d0 CLDW(I_, J_, 1:L_) = tau_clw_ij(k1:k2) CLDW(I_, J_, L1_) = 0.d0 CLDI(I_, J_, 1:L_) = tau_cli_ij(k1:k2) CLDI(I_, J_, L1_) = 0.d0 !if (londeg_i .ge. 70.d0 .and. londeg_i .lt. 72.d0 .and. & ! latdeg_j .ge. 44.0d0 .and. latdeg_j .lt. 46.d0) then !CLDWNDX(I_, J_, k1:k2) = 9 !CLDINDX(I_, J_, k1:k2) = 13 CLDWNDX(I_, J_, 1:L1_) = 9 CLDINDX(I_, J_, 1:L1_) = 13 ! fixend-by-hbian ! endBian, 11/21/2011 ! Set surface albedo SA(I_, J_) = surf_alb_ij ! Set Aerosol properties OPTDUST(1:L_,:) = ODMDUST_ij(k1:k2,:) OPTAER (1:L_,:) = ODAER_ij (k1:k2,:) OPTDUST(L1_,:) = 0.0d0 ! ODMDUST_ij(k1:k2,:) OPTAER (L1_,:) = 0.0d0 ! ODAER_ij (k1:k2,:) !!!!! testing only OPTDUST(:,:) = 0.0d0 OPTAER (:,:) = 0.0d0 ! Set pressure at the center of grid box PJ(1) = pctm_ij PJ(2:L1_) = press3e_ij(k1:k2) PJ(L1_+1) = 0.0d0 !XGRD(I_) = XDGRD(I_) *PI/180.d0 !YGRD(J_) = YDGRD(J_) *PI/180.d0 !--------------- ! Sets up clouds !--------------- !call SET_CLD (GMTAU) <------- This is empty at this time if (.not. present(ozone_ij)) then call SET_ATM0 endif !---reset the atmosphere, aerosols, clouds for the time step (now = dummy) !----------------------------------------------------------------------- if (present(ozone_ij)) then DO3(I_,J_,1:L_) = ozone_ij(k1:k2) DO3(I_,J_,L1_ ) = ozone_ij(k2) call SET_ATM(GMTAU) endif ! call SET_PROF <--- does not exist. Do we have to write it? !-------------------------------------- ! Compute the photolysis rate constants !-------------------------------------- ! GMTAU: Number days of simulation (in) ! IDAY: Current day (in) ! ILNG: (in) ! JLAT: (in) ! SOLF: Solar flux (out) ! SZA_ij: Solar zenith angle in degree (in) ! U0: cosine of the solar zenith angle or cos(SZA) (out) ! FREFL: fraction of flux (energy-wtd) reflected (out) ! ZPJ: 2D array of J's indexed to CTM chemistry (out) call PHOTOJ(GMTAU, IDAY, ILNG, JLAT, SOLF, SZA_ij, U0, FREFL, ZPJ) ! Store the photolysis rate constants qjgmi_ij(:,:) = ZPJ(:,:) !---------------------------------------------- ! Compute the Overhead ozone column diagnostics !---------------------------------------------- overheadO3col_ij = 0.0d0 overheadO3col_ij(k2) = DO3(I_,J_,k2) do il = k1, k2-1 do ik = il+1, k2 overheadO3col_ij(il) = overheadO3col_ij(il) + DO3(I_,J_,ik) end do end do RETURN end subroutine controlFastJX65 !EOC !--------------------------------------------------------------------------- !BOP ! ! !IROUTINE: initializeFastJX65 ! ! !INTERFACE: ! subroutine initializeFastJX65 (k1, k2, chem_mask_khi, num_qjs, & cross_section_file, rate_file, T_O3_climatology_file, rootProc) ! aerosolOpticalData_file, & ! cross_section_file , scattering_data_file, & ! rate_file, T_O3_climatology_file) ! # include "gmi_AerDust_const.h" # include "setkin_par.h" # include "parm_CTM_fastJX65.h" # include "cmn_JVdat_fastJX65.h" # include "cmn_metdat_fastJX65.h" ! ! !INPUT PARAMETERS: logical , intent(in) :: rootProc integer , intent(in) :: k1, k2 integer , intent(in) :: chem_mask_khi ! number of chemistry levels [JVL_] integer , intent(in) :: num_qjs ! number of photolysis reactions [JVN_] ! UMichigan aerosol optical data file ! character (len=128), intent(in) :: aerosolOpticalData_file ! fast-J X-sections (spectral data) input file name character (len=128), intent(in) :: cross_section_file ! Aerosol/cloud scattering data input file name ! character (len=128), intent(in) :: scattering_data_file ! Labels of photolysis rates required input file name ! keyed to chem code character (len=128), intent(in) :: rate_file ! Read in T & O3 climatology input file name ! general backup clim. character (len=128), intent(in) :: T_O3_climatology_file ! !EOP !---------------------------------------------------------------------- !BOC ! total number of cloud + aerosols + Rayleigh NCTA = 2 + NSADdust + NSADaer L_ = k2 - k1 + 1 L1_ = L_ + 1 L2_ = 2 * L_ + 2 JVL_ = chem_mask_khi JVN_ = num_qjs ! ! allocate(OPTDUSt(1:L_,NSADdust)) allocate(optDust(1:L1_,NSADdust)) optDust = 0.0d0 ! allocate(OPTAER (1:L_,NSADaer*nrh_b)) allocate(optAer (1:L1_,NSADaer*nrh_b)) optAer = 0.0d0 allocate (CLDW(I_,J_,L1_)) allocate (CLDI(I_,J_,L1_)) CLDW = 0.0d0 CLDI = 0.0d0 allocate (CLDWNDX(I_,J_,L1_)) allocate (CLDINDX(I_,J_,L1_)) !! testing only allocate (OOJdiag(I_,J_,L1_)) allocate (DDJdiag(I_,J_,L1_)) allocate (ODdiag(I_,J_,L1_,5)) OOJdiag = 0.0d0 DDJdiag = 0.0d0 ODdiag = 0.0d0 ! begBian, 11/21/2011 ! allocate(AER (NCTA,L1_)) ! AER = 0.0d0 ! allocate(ADX (NCTA,L1_)) ! allocate(AER1P_aero(1:L_,1:NSADaer)) ! AER1P_aero = 0.0d0 ! allocate(AER1P_dust(1:L_,1:NSADdust)) ! AER1P_dust = 0.0d0 ! Read in the input files ! call INPHOT (cross_section_file, & ! aerosolOpticalData_file, & ! scattering_data_file, & ! rate_file, & ! T_O3_climatology_file) call INPHOT (cross_section_file, & rate_file, T_O3_climatology_file, rootProc) ! bian (11/16/06) start ! ! Set up aerosol ! call SET_AER0 ! bian end return end subroutine initializeFastJX65 !EOC !--------------------------------------------------------------------------- !BOP ! ! !IROUTINE: getQAA_RAAinFastJX65 ! ! !INTERFACE: ! subroutine GetQAA_RAAinFastJX65 (bRAA, bQAA) ! # include "parm_CTM_fastJX65.h" # include "cmn_JVdat_fastJX65.h" ! ! !OUTPUT PARAMETERS: real*8 , intent(out) :: bRAA(:,:) real*8 , intent(out) :: bQAA(:,:) !EOP !--------------------------------------------------------------------------- !BOC ! Get effective radius associated with aerosol type bRAA(:,:) = RAA(:,:) !<----- Check the dimensions. used to be 2D now 1D ! Get aerosol scattering phase functions bQAA(:,:) = QAA(:,:) RETURN end subroutine getQAA_RAAinFastJX65 !EOC !--------------------------------------------------------------------------- ! >>>>>>>>>>>>>>>>current code revised to JX ver 6.5 (8/09)<<<<<<<<<<<< ! ! version 6.5 ! >>>> redo J-values to be average over CTM layer (not mid-layer) ! ! version 6.4 ! >>>> allows for shortened, speeded up troposphere versions<<<<<< ! STD: W_=18 ! identical results to v-6.2 if cloud OD is consistent ! TROP-ONLY: W_=12 ! collapses the wavelength bins from 18 to 12 (5-6-7-8 & 11-18) ! drops many 'stratospheric' cross-sections (denoted by 'x' in 2n ! allows use of single standard spectral data set: FJX_spec.dat ! results close to W_=18, largest difference is J-O2 (<1% in 13-1 ! This is recommended as accurate for troposphere only calculatio ! TROP-QUICK: W_=8 ! reverts to original fast-J 7-bins (12-18) plus 1 scaled UV (5) ! errors in 12-18 km range for J-O2, high sun are 10%, worse abov ! ***Photolysis of O2 in the upper tropical troposphere is an import ! source of O3. It needs to be included in tropospheric runs. ! TROP-ONLY is recommended, W_=8 is a quick fix if speed essentia ! ! Major rewrite of code to minimize calls and allow better vector-ty ! loop over wavelengths internal to Mie soln. ! Driven by profiling of CTM code, may still need optimization. ! Wavelengths can be contracted to W_=12 (trop only) and strat-only ! X-sections are dropped. With parm W_=18, the std fast-JX is re ! Many call eliminated and summations in BLKSLV and GEN_ID are expli ! GEN_ID replaces GEN and calculates all matrix coeff's (1:L_) at on ! RD_XXX changed to collapse wavelengths & x-sections to Trop-only: ! WX_ = 18 (parm_CTM.f) should match the JX_spec.dat wavelengt ! W_ = 12 (Trop-only) or 18 (std) is set in (parm_MIE.f). ! if W_=12 then drop strat wavels, and drop x-sects (e.g. N2O, ... ! ! version 6.3 ! revise cloud/aerosol OD & wavelength properties for CTM link: ! OPTICL is new sub for cloud optical properties, but it ! now starts with cloud OD @ 600 nm and cloud NDX ! fast-JX now uses cloud NDX to scale OD to other wavelengt ! OPTICA & OPTICM are new subs to convert aerosol path (g/m2) to ! A is std UCI scat data ! M is U Michigan data tables for aerosols, includes Rel Hu ! drop sub GAUSSP and put into Parameter statement (parm_MIE.f) ! ! version 6.2 ! corrects a long-standing problem at SZA > 89 degrees. ! In prior versions the ray-tracing of the path (and air-mass functi ! back to the sun was done at the edges of the CTM layers (it was de ! for the grid-point J-value code at Harvard/GISS/UCI). This left t ! interpolation to the mid-layer (needed for J's) open. The prior m ! gave irregular fluctuations in the direct solar beam at mid-layer ! large SZA > 88. This is now corrected with exact ray-tracing from ! the mid-pt of each CTM layer. For small SZA, there is no effectiv ! difference, for large SZA, results could be erratic. ! ! v-6.2 fix should be easy if you have migrated to v6.1, else some min ! caution may be needed: ! replace sub SPHERE with SPHERE2, AMF2 report factors for mid and ! replace sub OPMIE with new OPMIE, this uses the new AMF2 correctl ! replace sub PHOTOJ with new PHOTOJ, this just hands off AMF2 from ! SPHERE2 to OPMIE. ! ! version 6.1 adds ! 6.1b simplifies calling sequences feeds solar factor, albedo, to ! and read LAT, LNG directly. No substantive changes. ! new read-in of scat data for clouds/aerosols to allow for UMich d ! This has required substantial rewrite of some of the core subrout ! OPMIE is now called for each wavelength and without aersol/clo ! all subs below OPMIE are unchanged ! OPTICD & OPTICM are new subs to convert path (g/m2) to OD and ! D is std UCI scat data (re-ordered for clouds 1st) ! M is U Michigan data tables for aerosols, includes Rel Hu ! PHOTOJ now assembles the aerosol data (better for CTM implemen ! This version can reproduce earlier versions exactly, but the test ! is changed from OD and NDX to PATH (g/m2) and NDX. ! ! version 6.0 adds ! new 200-nm scattering data so that stratospheric aerosols can be ! version 5.7 ! adds the new flux diagnostics (including heating rates) ! accurate fluxes for spherical atmos and SZA > 90 ! ! recommend geometric delta-tau factor from 1.18 to 1.12 for more ac ! heating rates (but more layers!) ! tuned and corrected to be almost flux conserving (1.e-5), except ! deep clouds, where diffusive flux is created (1.e-4) ! still needs to return to the original 1970-code for the block-tri ! after extensive profiling with F95 and 'modern' versions ! it was found that they are much more expensive!!! ! corrects typo in JAC(2000) fast-J paper on I+ (reflected from l.b. ! I+(lb) = refl/(1+refl) * (4*Integ[j(lb)*mu*dmu] + mu0*Fdirect(l ! version 5.6 adds ! clean up problems with thick clouds does correct solar attenuatio ! into cloud sub-layers and into the mid-point of the CTM level ! New calculated upward and downward FLUXES at each wavelength at T ! Correct deposition of solar flux in each CTM layer (spherical) ! awaits new diagnostics of the h's for heating rates. ! back to old matrix solver (UCI blocksolver and matinv-4) ! version 5.5 adds ! new code for generating and solving the block tri-diagonal scatte ! problem. Uses single call to GEM and general 4x4 block-tri ! version 5.3c adds ! calculates reflected UV-vis solar energy (relative to 1.0) ! new solar spectrum (J-O2 increases in strat by 10%, J-NO by 15+%) ! ! version 5.3b changes include: ! new data files for specral Xsection and mie-scattering. ! add sub-layers (JXTRA) to thick cloud/aerosol layers, ! sets up log-spaced sub-layers of increasing thickness ATAU ! correction 'b' does massive clean up of the linking code, ! now the only subroutine that has access to CTM arrays is PHO ! Also, the access to the cmn_JVdat.f is 'read-only' after ini ! This should enable safe openMP/MPI coding. ! ! common files and what they mean: ! parm_CTM.f dimensions & params for code (CTM and fast-JX) ! parm_MIE.f dimensions for mie code variables. ! cmn_metdat.f CTM 3-D arrays, time of day, grid, etc. ! cmn_JVdat.f Xsects, Mie, etc., (initialized and then read-only) ! !<<<<<<<<<<<<<<<<<<<<>>will be superseded by CTM routines ! COMMON BLOCKS: cmn_JVdat.f ! ! RD_XXX(NJ1,NAMFIL) ! Read wavelength bins, solar fluxes, Rayleigh, Temp-dep X-sec ! >>>v-6.4 now collapse to Trop-only Wl's and Xs's if W_=12 ! called once by INPHOT. ! COMMON BLOCKS: cmn_JVdat.f ! Input files: FJX_spec.dat ! ! RD_MIE(NJ1,NAMFIL) ! Set aerosols/cloud scattering, called once by INPHOT ! COMMON BLOCKS: cmn_JVdat.f ! Input files: FJX_scat.dat ! ! RD_UM(NJ1,NAMFIL) ! UMich aerosol optical data, called once by INPHOT ! COMMON BLOCKS: cmn_JVdat.f ! Input files: FJX_UMaer.dat ! ! SOLARZ(GMTIME,NDAY,YGRDJ,XGRDI, SZA,COSSZA,SOLFX) ! SOLARZ(GMTIME,NDAY,YGRDJ,XGRDI, SZA,COSSZA,SOLFX) ! calc SZA and Solar Flux factor for given lat/lon/UT ! ! SPHERE2(GMU,RAD,ZHL,ZZHT,AMF2,L1_) ! calculate spherical geometry, air-mass factors (v 6.2) ! ! EXTRAL(DTAUX,L1X,L2X,NX,JTAUMX,ATAU,ATAU0, JXTRA) ! add sub-layers (JXTRA) to thick cloud/aerosol layers ! ! OPMIE (DTAUX,POMEGAX,U0,RFL,AMF2,JXTRA,FJACT,FJTOP,FJBOT,FSBOT,FJFLX ! calculate mean intensity (actinic) at each CTM levels ! calculate fluxes and deposition (heating rates) ! COMMON BLOCKS: cmn_JVdat.f ! !<<<<<<<<<<<<<<<<<<<<<<>>>> keyed to users che ! this is a tranfer map from the J's automatically calculated in fast- ! onto the names and order in the users chemistry code ! call RD_JS (IPH, 'chem_Js.dat') call RD_JS(IPH, rate_file, rootProc) ! ! Read in T & O3 climatology >>>> general backup clim ! call RD_PROF (IPH, 'atmos_std.dat') call RD_PROF(IPH, T_O3_climatology_file) ! return end subroutine INPHOT !###################################################################### subroutine RD_JS (NJ1, NAMFIL, rootProc) !----------------------------------------------------------------------- ! Reread the chem_Js.dat file to map photolysis rate to reaction ! Read in quantum yield 'jfacta' and fastj2 label 'jlabel' !----------------------------------------------------------------------- ! ! jfacta Quantum yield (or multiplication factor) for photolysis ! jlabel Reference label identifying appropriate J-value to use ! ipr Photolysis reaction counter - should total 'JVN_' ! !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" # include "parm_MIE_fastJX65.h" # include "cmn_metdat_fastJX65.h" # include "cmn_JVdat_fastJX65.h" ! integer , intent (in) ::NJ1 logical , intent (in) :: rootProc character ( * ), intent (in) ::NAMFIL integer :: IPR, I, J, K character (len=120) :: CLINE character (len=7) :: T_JX, T_CHEM ! Reread the chem_Js.dat file to map photolysis rate to reaction ! Read in quantum yield jfacta and fastj2 label jlab IF (rootProc) PRINT*," Reading the file: ", TRIM(NAMFIL) IPR = 0 open (NJ1, file = NAMFIL, status = 'old', form = 'formatted') 10 read (NJ1, '(A)', err = 20) CLINE if (IPR.eq.JVN_) goto 20 if (CLINE (2:5) .eq.'9999') then goto 20 elseif (CLINE (1:1) .eq.'#') then goto 10 elseif (CLINE (5:5) .eq.'$') then goto 10 else IPR = IPR + 1 read (CLINE (79:83) , '(F5.1)') JFACTA (IPR) read (CLINE (86:92) , '(A7)') JLABEL (IPR) JFACTA (IPR) = JFACTA (IPR) / 100.d0 goto 10 endif 20 close (NJ1) NRATJ = IPR !----------------------------------------------------------------------- ! compare Xsections titles with J-values listed in chem code (jratd.dat ! map J-values needed for chemistry (chem_Js.dat) onto the fast-JX rate ! >>>>>>>>>>>>>>>>current code revised to JPL-02 ver 8.5 (5/05)<<<<<<<< ! >>>this must now follow the read in of Xsects, etc<<< !----------------------------------------------------------------------- !---Zero / Set index arrays that map Jvalue(j) onto rates do J = 1, JVN_ JIND (J) = 0 enddo if(rootProc) then ! print *,'RD_JS JVN_ =',JVN_ ! print *,'RD_JS NJVAL =',NJVAL ! print *,'RD_JS TITLEJ =',TITLEJ(1:NJVAL) endif do J = 1, NJVAL T_JX = TITLEJ (J) ! Note: C3H6Ot renamed as Acet-a already !if (TITLEJ(J) .eq. 'Acet-a' .or. TITLEJ(J) .eq. 'Acet-b') T_JX = 'C3H6O' if (TITLEJ(J) .eq. 'Acet-a') T_JX = 'C3H6O' !!! don't delete comment for following line !if (TITLEJ(J) .eq. 'Acet-b') T_JX = 'C3H6O' do K = 1, NRATJ T_CHEM = JLABEL (K) if (T_CHEM (1:6) .eq.T_JX (1:6) ) JIND (K) = J enddo enddo IF (rootProc) THEN write (6, '(a,i4,a)') ' Photochemistry Scheme with ', IPR, & ' J-values' do K = 1, NRATJ J = JIND (K) if (J.eq.0) then write(6,'(i5,a9,f6.2,a,i4,a9)') K, JLABEL(K), JFACTA(K) & , ' has no mapping onto onto fast-JX' else write(6,'(i5,a9,f6.2,a,i4,a9)') K, JLABEL(K), JFACTA(K) & , ' mapped onto fast-JX:', J, TITLEJ (J) endif enddo END IF return end subroutine RD_JS !###################################################################### !----------------------------------------------------------------------- subroutine RD_PROF (NJ2, NAMFIL) !----------------------------------------------------------------------- ! Routine to input T and O3 reference profiles !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" # include "cmn_metdat_fastJX65.h" integer , intent (in) ::NJ2 character ( * ), intent (in) ::NAMFIL ! integer :: IA, I, M, L, LAT, MON, NTLATS, NTMONS, N216 real (8) :: OFAC, OFAK character (len=72) :: TITLE0 ! open (NJ2, file = NAMFIL, status = 'old', form = 'formatted') read (NJ2, '(A)') TITLE0 read (NJ2, '(2I5)') NTLATS, NTMONS ! write(6,'(1X,A)') TITLE0 !JENwrite (6, 1000) NTLATS, NTMONS N216 = min (216, NTLATS * NTMONS) do IA = 1, N216 read (NJ2, '(1X,I3,3X,I2)') LAT, MON M = min (12, max (1, MON) ) L = min (18, max (1, (LAT + 95) / 10) ) read (NJ2, '(3X,11F7.1)') (TREF (I, L, M) , I = 1, 41) read (NJ2, '(3X,11F7.4)') (OREF (I, L, M) , I = 1, 31) enddo close (NJ2) ! Extend climatology to 100 km OFAC = exp ( - 2.d5 / 5.d5) do I = 32, 51 OFAK = OFAC** (I - 31) do M = 1, NTMONS do L = 1, NTLATS OREF (I, L, M) = OREF (31, L, M) * OFAK enddo enddo enddo do L = 1, NTLATS do M = 1, NTMONS do I = 42, 51 TREF (I, L, M) = TREF (41, L, M) enddo enddo enddo 1000 format(1x,'std atmos profiles: ',i3,' lat x ',i2,' mon') return end subroutine RD_PROF !######################################################################## !----------------------------------------------------------------------- subroutine SET_ATM0 !----------------------------------------------------------------------- ! Routine to set up atmospheric profiles required by Fast-J2 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 supplie ! climatology and integrated to the CTM levels (may be overwritten with ! values directly from the CTM, if desired). ! Oliver (04/07/99) & MJP (7/05) !----------------------------------------------------------------------- ! ! ! PJ Pressure at boundaries of model levels (hPa) ! 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) ! ZH Altitude of boundaries of model levels (cm) ! PSTD Approximate pressures of levels for supplied climatology !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" # include "cmn_metdat_fastJX65.h" integer :: I, J, K, L, M, N real (8) :: PSTD (52), OREF2 (51), TREF2 (51) real (8) :: DLOGP, F0, T0, PB, PC, XC, MASFAC, SCALEH ! ! Select appropriate month M = max (1, min (12, MONTH) ) ! Mass factor - delta-Pressure (mbars) to delta-Column (molecules.cm-2) MASFAC = 100.d0 * 6.022d+23 / (28.97d0 * 9.8d0 * 10.d0) !----------------------------------------------------------------------- do J = 1, J_ ! Select appropriate latitudinal profiles N = max (1, min (18, (int (YDGRD (J) ) + 99) / 10) ) ! Temporary zonal arrays for climatology data do K = 1, 51 OREF2 (K) = OREF (K, N, M) TREF2 (K) = TREF (K, N, M) enddo do I = 1, I_ ! 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. ! L1_ = L_+1: ! PJ(L=1:L1_) = pressure at CTM layer edge, PJ(L1_+1)=0 (top-of-at !begJules ! PJ (L1_ + 1) = 0.d0 ! do K = 1, L1_ ! PJ (K) = ETAA (K) + ETAB (K) * PMEAN (I, J) ! enddo !endJules ! 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(ILNG,JLAT). ! 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 K = 3, 51 PSTD (K) = PSTD (K - 1) * DLOGP enddo PSTD (52) = 0.d0 do L = 1, L1_ F0 = 0.d0 T0 = 0.d0 do K = 1, 51 PC = min (PJ (L), PSTD (K) ) PB = max (PJ (L + 1), PSTD (K + 1) ) if (PC.gt.PB) then XC = (PC - PB) / (PJ (L) - PJ (L + 1) ) F0 = F0 + OREF2 (K) * XC T0 = T0 + TREF2 (K) * XC endif enddo TJ (I, J, L) = T0 DM (I, J, L) = (PJ (L) - PJ (L + 1) ) * MASFAC DO3 (I, J, L) = F0 * 1.d-6 * DM (I, J, L) enddo ! Calculate effective altitudes using scale height at each level ZH (I, J, 1) = 0.d0 do L = 1, L_ SCALEH = 1.3806d-19 * MASFAC * TJ (I, J, L) ZH (I, J, L + 1) = ZH (I, J, L) - (LOG (PJ (L + 1) / PJ (L) ) & * SCALEH) enddo enddo enddo return end subroutine SET_ATM0 !###################################################################### !<<<<<<<<<<<<<<<<<<<<<<< 98.0 deg => tan ht = 63 km ! or 99. 80 km if (SZA.gt.SZAMAX) goto 99 !---load the amtospheric column data ! begBian, 11/21/2011 do L = 1, L1_ ! PPJ (L) = ETAA (L) + ETAB (L) * P (ILNG, JLAT) PPJ (L) = PJ (L) TTJ (L) = TJ (ILNG, JLAT, L) DDJ (L) = DM (ILNG, JLAT, L) ZZJ (L) = DO3 (ILNG, JLAT, L) enddo ! endBian, 11/21/2011 PPJ (L1_ + 1) = 0.d0 !---calculate spherical weighting functions (AMF: Air Mass Factor) do L = 1, L1_ ZHL (L) = ZH (ILNG, JLAT, L) enddo ZHL (L1_ + 1) = ZHL (L1_) + ZZHT !----------------------------------------------------------------------- call SPHERE2 (U0, RAD, ZHL, ZZHT, AMF2, L1_) !----------------------------------------------------------------------- !---calculate the optical properties (opt-depth, single-scat-alb, phase- !--- at the 5 std wavelengths 200-300-400-600-999 nm for cloud+aerosols do L = 1, L1_ OD600 (L) = 0.D0 !do K = 1, 5 do K = 1, numScat-1 OD (K, L) = 0.d0 SSA (K, L) = 0.d0 do I = 1, 8 SLEG (I, K, L) = 0.d0 enddo enddo enddo do L = 1, L1_ !---cloud in layer: if valid cloud index (now 4:13) ! begBian, 11/22/2011 do M = 1, 2 ! NDCLD = CLDNDX (ILNG, JLAT, L) ! PATH = CLDLWP (ILNG, JLAT, L) ! DENS = PATH / ZH (ILNG, JLAT, L) if (M .eq. 1) then NDCLD = CLDWNDX (ILNG, JLAT, L) ODCLD = CLDW (ILNG, JLAT, L) else NDCLD = CLDINDX (ILNG, JLAT, L) ODCLD = CLDI (ILNG, JLAT, L) endif RH = RELH (L) ! endBian, 11/22/2011 if (ODCLD.gt.0.d0) then if (NDCLD.lt.4.or.NDCLD.gt.13) then NDCLD = 9 endif !---now calculate cloud OD at 600 nm from mixed formulae: !REFF = RAA (NDCLD) !RHO = DAA (NDCLD) ! begBian, 11/22/2011 ! if (NDCLD.eq.12.or.NDCLD.eq.13) then !!---Ice Water Content (density at g/m3) varies from 0.0001 to 0.1 !!---correct R-eff for ice clouds, increase Reff (microns) with density ! REFF = 50.d0 * (1.d0 + 8.333d0 * DENS) ! endif !!---extinction K(m2/g) = Q(wvl) / [4/3 * Reff(micron) * aerosol-density( ! ODCLD = PATH * 0.75d0 * QAA (4, NDCLD) / (REFF * RHO) ! endBian, 11/22/2011 call OPTICL (OPTX, SSAX, SLEGX, ODCLD, NDCLD) ! begBian, 11/22/2011 ! do K = 1, 5 do K = 1, numScat-1 ! endBian, 11/22/2011 OD (K, L) = OD (K, L) + OPTX (K) SSA (K, L) = SSA (K, L) + SSAX (K) * OPTX (K) do I = 1, 8 SLEG(I,K,L) = SLEG(I,K,L) + SLEGX(I,K) * SSAX(K) * OPTX(K) enddo enddo endif !---use OD of clouds (not aerosols) at 600 nm to determine added layers ! begBian, 11/22/2011 enddo ! m ! OD600 (L) = OD (4, L) OD600 (L) = OD (3, L) ! endBian, 11/22/2011 !---aerosols in layer: check aerosol index !---this uses data from climatology OR from current CTM (STT of aerosols !---FIND useful way to sum over different aerosol types! ! !============================================================== ! Loop over the no. of aerosol/cloud types supplied from CTM ! begBian, 11/22/2011 ! do M = 1, MX do M = 1, NSADdust+NSADaer ! NAER = AERindex (M) ! it = 3 ! if (M .ge. it) then ! PATH = 0.0d0 ! end if ! ! For mineral dust ! if ((M .gt. it) .and. (M .ge. it + NSADdust)) then if (M .le. NSADdust) then ODaer = OPTdust (L, M) NAER = 14 + M !NAER = M call OPTICA (OPTX, SSAX, SLEGX, ODaer, RH, NAER) ! do K = 1, 5 do K = 1, numScat-1 OD (K, L) = OD (K, L) + OPTX (K) SSA(K,L) = SSA(K,L) + SSAX(K) * OPTX(K) do I = 1, 8 SLEG(I,K,L) = SLEG(I,K,L) + SLEGX(I,K) * SSAX(K) * OPTX (K) enddo enddo else do MS = 1, NRH_b MSOD = (M-NSADdust-1)*5+MS MSN = 14+NSADdust+(M-NSADdust-1)*7+MS !MSOD = (M-NSADdust-1)*NRH_b + MS !MSN = NSADdust+(M-NSADdust-1)*NRH_b+MS ODaer = OPTaer (L, MSOD) NAER = MSN call OPTICA (OPTX, SSAX, SLEGX, ODaer, RH, NAER) ! do K = 1, 5 do K = 1, numScat-1 OD (K, L) = OD (K, L) + OPTX (K) SSA (K, L) = SSA (K, L) + SSAX (K) * OPTX (K) do I = 1, 8 SLEG(I,K,L) = SLEG(I,K,L) + SLEGX(I,K) * SSAX(K) * OPTX (K) enddo enddo enddo endif ! PATH = AER1P_dust (L, M-it) ! call OPTICA (OPTX, SSAX, SLEGX, PATH, RH, NAER) ! end if ! ! For other aerosols (sulfate, sea salt, etc.) ! if ((M .gt. it + NSADdust) .and. (M .ge. it + NSADdust+NSADaer)) then ! PATH = AER1P_aero (L, M-(it+NSADdust)) ! call OPTICM (OPTX, SSAX, SLEGX, PATH, RH, NAERnew) ! end if !!end JulesKouatchou !!============================================================== ! !!---subroutines OPTICA & OPTICM return the same information: !!--- optical depth (OPTX), single-scat albedo (SSAX) and phase fn (SLEG !!---subs have slightly different inputs: !!--- PATH is the g/m2 in the layer, NAER in the cloud/aerosol index !!--- UMich aerosols use relative humidity (RH) ! if (PATH.gt.0.d0) then !!KNJR 30Sept2011 if (NAER.gt.0) then !!KNJR 30Sept2011 call OPTICA (OPTX, SSAX, SLEGX, PATH, RH, NAER) !!KNJR 30Sept2011 else !!KNJR 30Sept2011 NAERnew = -NAER !!KNJR 30Sept2011 call OPTICM (OPTX, SSAX, SLEGX, PATH, RH, NAERnew) !!KNJR 30Sept2011 endif ! do K = 1, 5 ! OD (K, L) = OD (K, L) + OPTX (K) ! SSA (K, L) = SSA (K, L) + SSAX (K) * OPTX (K) ! do I = 1, 8 ! SLEG (I, K, L) = SLEG (I, K, L) + SLEGX (I, K) * SSAX (K) * OPTX (K) ! enddo ! enddo ! endif enddo ! M ! endBian, 11/22/2011 ! begBian, 11/22/2011 ! do K = 1, 5 do K = 1, numScat-1 ! endBian, 11/22/2011 if (OD (K, L) .gt.0.d0) then SSA (K, L) = SSA (K, L) / OD (K, L) do I = 1, 8 SLEG (I, K, L) = SLEG (I, K, L) / OD (K, L) enddo endif enddo enddo ! L !---can add aerosol OD at 600 nm to determine added layers, but not done !--- OD600(L) = OD(4,L) + aerosol OD's !---when combining with Rayleigh and O2-O3 abs, remember the SSA and !--- phase fn SLEG are weighted by OD and OD*SSA, respectively. !---Given the aerosol+cloud OD/layer in visible (600 nm) calculate how t ! additonal levels at top of clouds (now uses log spacing) !----------------------------------------------------------------------- call EXTRAL (OD600, L1_, L2_, N_, JTAUMX, ATAU, ATAU0, JXTRA) !----------------------------------------------------------------------- !---set surface reflectance RFLECT = SA (ILNG, JLAT) RFL (:) = max (0.d0, min (1.d0, RFLECT) ) !----------------------------------------------------------------------- !---Loop over all wavelength bins to calc mean actinic flux AVGF(L) !----------------------------------------------------------------------- do K = 1, W_ WAVE = WL (K) !---Pick nearest Mie wavelength to get scattering properites------------ ! use 200 nm prop for <255 nm KMIE = 1 ! use 300 nm prop for 255-355 nm ! if (WAVE.gt.255.d0) KMIE = 2 ! use 400 nm prop for 355-500 nm ! if (WAVE.gt.355.d0) KMIE = 3 ! if (WAVE.gt.500.d0) KMIE = 4 ! if (WAVE.gt.800.d0) KMIE = 5 if (WAVE.gt.355.d0) KMIE = 2 if (WAVE.gt.500.d0) KMIE = 3 if (WAVE.gt.800.d0) KMIE = 4 !---Combine: Rayleigh scatters & O2 & O3 absorbers to get optical proper !---values at L1_=L_+1 are a pseudo/climatol layer above the top CTM lay do L = 1, L1_ TTT = TTJ (L) XQO3 = FLINT (TTT, TQQ (1, 2), TQQ (2, 2), TQQ (3, 2), QO3 (K, 1), & QO3 (K, 2), QO3 (K, 3) ) XQO2 = FLINT (TTT, TQQ (1, 1), TQQ (2, 1), TQQ (3, 1), QO2 (K, 1), & QO2 (K, 2), QO2 (K, 3) ) ! beg biandix, Dec19, 2011 ! this line was missed in previous implementation ODABS = XQO3*ZZJ(L) + XQO2*DDJ(L)*0.20948d0 ! end biandix, Dec19, 2011 ODRAY = DDJ (L) * QRAYL (K) DTAUX (L, K) = OD (KMIE, L) + ODABS + ODRAY do I = 1, 8 POMEGAX (I, L, K) = SLEG (I, KMIE, L) * OD (KMIE, L) enddo POMEGAX (1, L, K) = POMEGAX (1, L, K) + 1.0d0 * ODRAY POMEGAX (3, L, K) = POMEGAX (3, L, K) + 0.5d0 * ODRAY do I = 1, 8 POMEGAX (I, L, K) = POMEGAX (I, L, K) / DTAUX (L, K) enddo enddo enddo !!! testing only: diagnostics !write(*,*) 'now filling FAST-JX diagnostics A: ', sum(od),sum(zzj),sum(ddj) !write(*,*) 'now filling FAST-JX diagnostics B: ', size(od), size(oddiag) do L=1,L1_ do k=1,4 if(allocated(ODdiag)) ODdiag(ilng,jlat,l,k)=od(k,l)/1.0d18 enddo if(allocated(OOJdiag)) OOJdiag(ilng,jlat,l)=zzj(l)/1.0d18 if(allocated(DDJdiag)) DDJdiag(ilng,jlat,l)=ddj(l)/1.0d18 enddo !write(*,*) 'FAST-JX diagnostics filled!' !----------------------------------------------------------------------- call OPMIE (DTAUX, POMEGAX, U0, RFL, AMF2, JXTRA, AVGF, FJTOP, & FJBOT, FSBOT, FJFLX, FLXD, FLXD0) !----------------------------------------------------------------------- do K = 1, W_ !----direct(DIR) and diffuse(FLX) fluxes at top(UP) (solar = negative by !---- also at bottom (DN), does not include diffuse reflected flux. FLXUP (K) = FJTOP (K) DIRUP (K) = - FLXD0 (K) FLXDN (K) = - FJBOT (K) DIRDN (K) = - FSBOT (K) do L = 1, JVL_ FFF (K, L) = FFF (K, L) + SOLF * FL (K) * AVGF (L, K) enddo FREFI = FREFI + SOLF * FL (K) * FLXD0 (K) / WL (K) FREFL = FREFL + SOLF * FL (K) * FJTOP (K) / WL (K) FREFS = FREFS + SOLF * FL (K) / WL (K) !---for each wavelength calculate the flux budget/heating rates: ! FLXD(L) = direct flux deposited in layer L [approx = MU0*(F(L+1) -F( ! but for spherical atmosphere! ! FJFLX(L) = diffuse flux across top of layer L !---calculate divergence of diffuse flux in each CTM layer (& t-o-a) !--- need special fix at top and bottom: !---FABOT = total abs at L.B. & FXBOT = net diffusive flux at L.B. FABOT = (1.d0 - RFL (K) ) * (FJBOT (K) + FSBOT (K) ) FXBOT = - FJBOT (K) + RFL (K) * (FJBOT (K) + FSBOT (K) ) FLXJ (1) = FJFLX (1, K) - FXBOT do L = 2, L_ FLXJ (L) = FJFLX (L, K) - FJFLX (L - 1, K) enddo FLXJ (L_ + 1) = FJTOP (K) - FJFLX (L_, K) !---calculate net flux deposited in each CTM layer (direct & diffuse): FFX0 = 0.d0 do L = 1, L1_ FFX (K, L) = FLXD (L, K) - FLXJ (L) FFX0 = FFX0 + FFX (K, L) enddo ! NB: the radiation level ABOVE the top CTM level is included in these ! these are the flux budget/heating terms for the column: ! FFXNET(K,1) = FLXD0 direct(solar) flux dep into atmos (spheric ! FFXNET(K,2) = FSBOT direct(solar) flux dep onto LB (surface) ! FFXNET(K,3) = FLXD0+FSBOT TOTAL solar into atmopshere+surface ! FFXNET(K,4) = FJTOP diffuse flux leaving top-of-atmos ! FFXNET(K,5) = FFX0 diffuse flux absorbed in atmos ! FFXNET(K,6) = FABOT total (dir+dif) absorbed at LB (surface) ! these are surface fluxes to compare direct vs. diffuse: ! FFXNET(K,7) = FSBOT direct flux dep onto LB (surface) - for sr ! FFXNET(K,8) = FJBOT diffuse flux dep onto LB (surface) FFXNET (K, 1) = FLXD0 (K) FFXNET (K, 2) = FSBOT (K) FFXNET (K, 3) = FLXD0 (K) + FSBOT (K) FFXNET (K, 4) = FJTOP (K) FFXNET (K, 5) = FFX0 FFXNET (K, 6) = FABOT FFXNET (K, 7) = FSBOT (K) FFXNET (K, 8) = FJBOT (K) !----------------------------------------------------------------------- ! end loop over wavelength K enddo !----------------------------------------------------------------------- !calculate reflected flux (energy wei FREFL = FREFL / FREFS FREFI = FREFI / FREFS !---NB UVB = 280-320 = bins 12:15, UVA = 320-400 = bins 16:17, VIS = bin !----------------------------------------------------------------------- call JRATET (PPJ, TTJ, FFF, VALJL) !----------------------------------------------------------------------- !---map the J-values from fast-JX onto ASAD ones (use JIND & JFACTA) do L = 1, JVL_ do J = 1, NRATJ if (JIND (J) .gt.0) then ZPJ (L, J) = VALJL (L, JIND (J) ) * JFACTA (J) else ZPJ (L, J) = 0.d0 endif enddo enddo !---diagnostics that are NOT returned to the CTM code !----------------------------------------------------------------------- ! write (6, * ) 'fast-JX-(6.5)----PHOTOJ internal print: Atmosphere & !&----' !!---used last called values of DTAUX and POMEGAX, should be 600 nm ! call JP_ATM (PPJ, TTJ, DDJ, ZZJ, ZHL, ZZHT, DTAUX (1, W_), & POMEGAX (1, 1, W_), JXTRA) !!---PRINT SUMMARY of mean intensity, flux, heating rates: ! write (6, * ) ! write (6, * ) 'fast-JX(6.5)----PHOTOJ internal print: Mean Intens & !&----' ! ! write (6, '(a,5f10.4)') ' SUMMARY fast-JX: albedo/SZA/u0/F-incd/F- & ! &refl/', RFLECT, SZA, U0, FREFI, FREFL ! write (6, '(a5,18i8)') ' bin:', (K, K = NW2, NW1, - 1) ! write (6, '(a5,18f8.1)') ' wvl:', (WL (K) , K = NW2, NW1, - 1) ! write (6, '(a)') ' ---- 100000=Fsolar MEAN INTENSITY per wvl bin' do L = JVL_, 1, - 1 do K = NW1, NW2 RATIO (K) = (1.d5 * FFF (K, L) / FL (K) ) enddo ! write (6, '(i3,2x,18i8)') L, (RATIO (K) , K = NW2, NW1, - 1) enddo ! write (6, * ) ! write (6, * ) 'fast-JX(6.5)----PHOTOJ internal print: Net Fluxes----' ! write (6, '(a11,18i8)') ' bin:', (K, K = NW2, NW1, - 1) ! write (6, '(a11,18f8.1)') ' wvl:', (WL (K) , K = NW2, NW1, - 1) !! write(6,'(a11,18f8.4)') ' sol in atm',(FFXNET(K,1), K=NW2,NW1,-1) !! write(6,'(a11,18f8.4)') ' sol at srf',(FFXNET(K,2), K=NW2,NW1,-1) ! write (6, * ) ' ---NET FLUXES--- ' ! write (6, '(a11,18f8.4)') ' sol TOTAL ', (FFXNET (K, 3) , K = NW2, NW1, - 1) ! write (6, '(a11,18f8.4)') ' dif outtop', (FFXNET (K, 4) , K = NW2, NW1, - 1) ! write (6, '(a11,18f8.4)') ' abs in atm', (FFXNET (K, 5) , K = NW2, NW1, - 1) ! write (6, '(a11,18f8.4)') ' abs at srf', (FFXNET (K, 6) , K = NW2, NW1, - 1) ! write (6, * ) ' ---SRF FLUXES--- ' ! write (6, '(a11,18f8.4)') ' srf direct', (FFXNET (K, 7) , K = NW2, NW1, - 1) ! write (6, '(a11,18f8.4)') ' srf diffus', (FFXNET (K, 8) , K = NW2, NW1, - 1) ! write (6, '(2a)') ' ---NET ABS per layer: 10000=Fsolar', ' & ! & [NB: values <0 = numerical error w/clouds or SZA>90, colm OK]' do L = JVL_, 1, - 1 do K = NW1, NW2 RATIO (K) = 1.d5 * FFX (K, L) enddo ! write (6, '(i9,2x,18i8)') L, (RATIO (K) , K = NW2, NW1, - 1) enddo !----------------------------------------------------------------------- 99 continue return end subroutine PHOTOJ !###################################################################### !<<<<<<<<<<<<<<<<<<<<<<1.0 real (8) , intent (in) ::RELH ! index of cloud/aerosols integer , intent (inout) :: L integer :: I, J, LL real (8) :: XTINCT, REFF, RHO !if (L.gt.NAA.or.L.lt.4) then ! !write (6, * ) ' aerosol index out-of-range: L/NAA', L, NAA ! L = 18 !endif if (L.gt.NAA) then !write (6, * ) ' aerosol index out-of-range: L/NAA', L, NAA LL = L-(L-50)/7*7 L = LL endif ! if (L.lt.14) then ! write (6, * ) ' aerosol as cloud: L/NAA', L, NAA ! endif !REFF = RAA (L) !RHO = DAA (L) ! begBian, 11/22/2011 ! do J = 1, 5 do J = 1, numScat-1 !---extinction K(m2/g) = Q(wvl) / [4/3 * Reff(micron) * aerosol-density( ! XTINCT = 0.75d0 * QAA (J, L) / (REFF * RHO) ! OPTD (J) = PATH * XTINCT ! note: ODaer was calculated at k = 4 OPTD(J) = ODaer * QAA(J,L) / QAA(4, L) ! endBian, 11/22/2011 SSALB (J) = SAA (J, L) do I = 1, 8 SLEG (I, J) = PAA (I, J, L) enddo enddo return end subroutine OPTICA !###################################################################### !----------------------------------------------------------------------- subroutine OPTICM (OPTD, SSALB, SLEG, PATH, RELH, L) !----------------------------------------------------------------------- !---for the U Michigan aerosol data sets, this generate fast-JX data for !---NB Approximates the Legendre expansion(L) of the scattering phase fn !--- as (2*L+1)*g**L !---UMAER(I,J,K,L): ! I=1:3 = [SSAbldeo, g, k-ext(m2/g)] ! J=1:6 = [200, 300, 400, 550, 600 , 1000 nm] ! K=1:21= [0, 5, 10, 15, ..., 90, 95, 99 %RelHum] ! L=1:33= UM aerosol types [SULF, SS-1,-2,-3,-4, DD-1,-2,-3,-4, FF00(0 ! FF02, ...FF14(14%BC), BB00, BB02, ...BB30(30%BC)] implicit none # include "parm_CTM_fastJX65.h" # include "parm_MIE_fastJX65.h" # include "cmn_JVdat_fastJX65.h" ! optical depth of layer real (8) , intent (out) ::OPTD (5) ! single-scattering albedo real (8) , intent (out) ::SSALB (5) ! scatt phase fn (Leg coeffs) real (8) , intent (out) ::SLEG (8, 5) ! path (g/m2) of aerosol/cloud real (8) , intent (in) ::PATH ! relative humidity (0.00->1.0 real (8) , intent (in) ::RELH ! index of cloud/aerosols integer , intent (inout) ::L integer :: KR, J real (8) :: R, FRH, GCOS, XTINCT !---calculate fast-JX properties at the std 5 wavelengths:200-300-400-60 !---interpolate in Relative Humidity !---extrapolate pahse fn from first term (g) if (L.gt.33) then write (6, * ) ' UM aer index too large: L', L L = 1 endif R = 100.d0 * min (1.d0, max (.01d0, RELH) ) KR = (R / 5.d0) KR = max (0, min (19, KR) ) + 1 if (KR.lt.20) then FRH = 0.20d0 * (R - 5.d0 * float (KR - 1) ) else FRH = 0.25d0 * (R - 5.d0 * float (KR - 1) ) endif do J = 1, 5 SSALB (J) = UMAER (1, J, KR, L) * (1.d0 - FRH) + UMAER (1, J, KR + & 1, L) * FRH XTINCT = UMAER (3, J, KR, L) * (1.d0 - FRH) + UMAER (3, J, KR + 1, & L) * FRH OPTD (J) = PATH * XTINCT GCOS = UMAER (2, J, KR, L) * (1.d0 - FRH) + UMAER (2, J, KR + 1, & L) * FRH SLEG (1, J) = 1.d0 SLEG (2, J) = 3.d0 * GCOS SLEG (3, J) = 5.d0 * GCOS**2 SLEG (4, J) = 7.d0 * GCOS**3 SLEG (5, J) = 9.d0 * GCOS**4 SLEG (6, J) = 11.d0 * GCOS**5 SLEG (7, J) = 13.d0 * GCOS**6 SLEG (8, J) = 15.d0 * GCOS**7 enddo return end subroutine OPTICM !###################################################################### !----------------------------------------------------------------------- subroutine JRATET (PPJ, TTJ, FFF, VALJL) !----------------------------------------------------------------------- ! in: ! PPJ(L1_+1) = pressure profile at edges ! TTJ(L1_) = = temperatures at mid-level ! FFF(K=1:NW, L=1:JVL_) = mean actinic flux ! out: ! VALJ(JVL_,NJVAL) JVL_ = no of levels !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" # include "parm_MIE_fastJX65.h" # include "cmn_JVdat_fastJX65.h" real (8) , intent (in) ::PPJ (L1_ + 1), TTJ (L1_) real (8) , intent (in) ::FFF (W_, JVL_) real (8) , intent (out) ::VALJL (JVL_, NJVAL) ! temp for call J's at one L real (8) :: VALJ (X_) real (8) :: QO2TOT (W_), QO3TOT (W_), QO31DY (W_), QO31D, QQQT, & TFACT real (8) :: TT, PP, DD, TT200, TFACA, TFAC0, TFAC1, TFAC2, QQQA, & QQ2, QQ1A, QQ1B integer :: J, K, L, IV ! master loop over layer = L do L = 1, JVL_ !---need temperature, pressure, and density at mid-layer (for some quant TT = TTJ (L) if (L.eq.1) then PP = PPJ (1) else PP = (PPJ (L) + PPJ (L + 1) ) * 0.5d0 endif DD = 7.24e18 * PP / TT do J = 1, NJVAL VALJ (J) = 0.d0 enddo if (TT.le.TQQ (2, 1) ) then if (TT.le.TQQ (1, 1) ) then TFACT = 0.d0 else TFACT = (TT - TQQ (1, 1) ) / (TQQ (2, 1) - TQQ (1, 1) ) endif do K = 1, W_ QO2TOT (K) = QO2 (K, 1) + (QO2 (K, 2) - QO2 (K, 1) ) * TFACT enddo else if (TT.ge.TQQ (3, 1) ) then TFACT = 1.d0 else TFACT = (TT - TQQ (2, 1) ) / (TQQ (3, 1) - TQQ (2, 1) ) endif do K = 1, W_ QO2TOT (K) = QO2 (K, 2) + (QO2 (K, 3) - QO2 (K, 2) ) * TFACT enddo endif if (TT.le.TQQ (2, 2) ) then if (TT.le.TQQ (1, 2) ) then TFACT = 0.d0 else TFACT = (TT - TQQ (1, 2) ) / (TQQ (2, 2) - TQQ (1, 2) ) endif do K = 1, W_ QO3TOT (K) = QO3 (K, 1) + (QO3 (K, 2) - QO3 (K, 1) ) * TFACT enddo else if (TT.ge.TQQ (3, 2) ) then TFACT = 1.d0 else TFACT = (TT - TQQ (2, 2) ) / (TQQ (3, 2) - TQQ (2, 2) ) endif do K = 1, W_ QO3TOT (K) = QO3 (K, 2) + (QO3 (K, 3) - QO3 (K, 2) ) * TFACT enddo endif if (TT.le.TQQ (2, 3) ) then if (TT.le.TQQ (1, 3) ) then TFACT = 0.d0 else TFACT = (TT - TQQ (1, 3) ) / (TQQ (2, 3) - TQQ (1, 3) ) endif do K = 1, W_ QO31DY (K) = Q1D (K, 1) + (Q1D (K, 2) - Q1D (K, 1) ) * TFACT enddo else if (TT.ge.TQQ (3, 3) ) then TFACT = 1.d0 else TFACT = (TT - TQQ (2, 3) ) / (TQQ (3, 3) - TQQ (2, 3) ) endif do K = 1, W_ QO31DY (K) = Q1D (K, 2) + (Q1D (K, 3) - Q1D (K, 2) ) * TFACT enddo endif do K = 1, W_ QO31D = QO31DY (K) * QO3TOT (K) VALJ (1) = VALJ (1) + QO2TOT (K) * FFF (K, L) VALJ (2) = VALJ (2) + QO3TOT (K) * FFF (K, L) VALJ (3) = VALJ (3) + QO31D * FFF (K, L) enddo ! begBian, 1/07/2012 !cc ***NB*** original fast-J/2 code had VALJ(2) as J[O3 -> O(3P)+O2] !cc which is different from std pratmo code - fixed 7/04 MP !cc QO33P = QO3TOT - QO31D VALJ(2) = VALJ(2) - VALJ(3) !AOO for GMI per the above comment ! endBian, 1/07/2012 do J = 4, NJVAL if (TQQ (2, J) .gt.TQQ (1, J) ) then TFACT = max (0.d0, min (1.d0, (TT - TQQ (1, J) ) / (TQQ (2, & J) - TQQ (1, J) ) ) ) else TFACT = 0.d0 endif do K = 1, W_ QQQT = QQQ (K, 1, J) + (QQQ (K, 2, J) - QQQ (K, 1, J) ) * TFACT VALJ (J) = VALJ (J) + QQQT * FFF (K, L) enddo !>>>>Methylvinyl ketone 'MeVK ' q(M) = 1/(1 + 1.67e-19*[M]) if (TITLEJ (J) .eq.'MeVK ') then VALJ (J) = VALJ (J) / (1.0 + 1.67e-19 * DD) endif !>>>>Methylethyl ketone MEKeto q(M) = 1/(1 + 2.0*[M/2.5e19]) if (TITLEJ (J) .eq.'MEKeto') then VALJ (J) = VALJ (J) / (1.0 + 0.80E-19 * DD) endif !>>>>Methyl glyoxal MGlyxl q(M) = 1/(1 + 4.15*[M/2.5E19]) if (TITLEJ (J) .eq.'MGlyxl') then VALJ (J) = VALJ (J) / (1.0 + 1.66e-19 * DD) endif enddo if (TITLEJ (NJVAL - 1) .eq.'Acet-a') then !--------------J-ref v8.3 includes Blitz ACETONE q-yields-------------- !---Acetone is a special case: (as per Blitz et al GRL, 2004) !--- 61 = NJVAL-1 = J1(acetone-a) ==> CH3CO + CH3 !--- 62 = NJVAL = J2(acetone-b) ==> CH3 + CO + CH3 VALJ (NJVAL - 1) = 0.d0 VALJ (NJVAL) = 0.d0 !---IV=NJVAL-1 = Xsect (total abs) for Acetone - pre-calc Temp interp fa IV = NJVAL - 1 TFACA = (TT - TQQ (1, IV) ) / (TQQ (2, IV) - TQQ (1, IV) ) TFACA = max (0.d0, min (1.d0, TFACA) ) !---IV=NJVAL = Q2 for Acetone=>(2), specifically designed for quadratic !--- but force to Q2=0 by 210K IV = NJVAL TFAC0 = ( (TT - TQQ (1, IV) ) / (TQQ (2, IV) - TQQ (1, IV) ) & ) **2 if (TT.lt.TQQ (1, IV) ) then TFAC0 = (TT - 210.d0) / (TQQ (1, IV) - 210.d0) endif TFAC0 = max (0.d0, min (1.d0, TFAC0) ) !---IV=NJVAL+1 = Q1A for Acetone => (1), allow full range of T = 200K-30 IV = NJVAL + 1 TT200 = min (300.d0, max (200.d0, TT) ) TFAC1 = (TT200 - TQQ (1, IV) ) / (TQQ (2, IV) - TQQ (1, IV) ) !---IV=NJVAL+2 = Q1B for Acetone => (1) IV = NJVAL + 2 TFAC2 = (TT200 - TQQ (1, IV) ) / (TQQ (2, IV) - TQQ (1, IV) ) !---now integrate over wavelengths do K = 1, W_ !---NJVAL-1 = Xsect (total abs) for Acetone IV = NJVAL - 1 QQQA = QQQ (K, 1, IV) + (QQQ (K, 2, IV) - QQQ (K, 1, IV) ) & * TFACA !---NJVAL = Q2 for Acetone=>(2), specifically designed for quadratic i IV = NJVAL QQ2 = QQQ (K, 1, IV) + (QQQ (K, 2, IV) - QQQ (K, 1, IV) ) & * TFAC0 if (TT.lt.TQQ (1, IV) ) then QQ2 = QQQ (K, 1, IV) * TFAC0 endif !---NJVAL+1 = Q1A for Acetone => (1), allow full range of T = 200K-300K IV = NJVAL + 1 QQ1A = QQQ (K, 1, IV) + (QQQ (K, 2, IV) - QQQ (K, 1, IV) ) & * TFAC1 !---NJVAL+2 = Q1B for Acetone => (1) ! scaled to [M]=2.5e19 IV = NJVAL + 2 QQ1B = QQQ (K, 1, IV) + (QQQ (K, 2, IV) - QQQ (K, 1, IV) ) & * TFAC2 QQ1B = QQ1B * 4.d-20 !---J(61) VALJ (NJVAL - 1) = VALJ (NJVAL - 1) + FFF (K, L) * QQQA * & (1.d0 - QQ2) / (QQ1A + QQ1B * DD) !---J(62) VALJ (NJVAL) = VALJ (NJVAL) + FFF (K, L) * QQQA * QQ2 !K enddo !-----------end v-8.3 includes Blitz ACETONE q-yields-------------- endif !----Load array of J-values in native order, need to be indexed/scaled ! by ASAD-related code later: ZPJ(L,JJ) = VALJL(L,JIND(JJ))*JFACTA(JJ ! begBian, 2/08/2012 ! add two acetone reaction branches together VALJ(NJVAL-1) = VALJ(NJVAL-1)+VALJ(NJVAL) ! endBian, 2/08/2012 do J = 1, NJVAL VALJL (L, J) = VALJ (J) enddo ! master loop over L=1,JVL_ enddo return end subroutine JRATET !###################################################################### !----------------------------------------------------------------------- subroutine JP_ATM (PPJ, TTJ, DDJ, ZZJ, ZHL, ZZHT, DTAUX, POMEGAX, & JXTRA) !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" !----------------------------------------------------------------------- !--------key amtospheric data needed to solve plane-parallel J---------- real (8) , dimension (L1_ + 1) ::TTJ, DDJ, ZZJ, ZHL real (8) , dimension (L1_ + 1) ::PPJ integer , dimension (L2_ + 1) ::JXTRA real (8) ::ZZHT real (8) ::DTAUX (L1_), POMEGAX (8, L1_) !----------------------------------------------------------------------- integer :: I, J, K, L real (8) :: COLO2, COLO3, ZKM, DELZ, ZTOP ! write(6,'(4a)') ' L z(km) p T ', & ! & ' d(air) d(O3)',' col(O2) col(O3) d-TAU SS-alb', & ! & ' g(cos) CTM lyr=>' COLO2 = 0.d0 COLO3 = 0.d0 ZTOP = ZHL (L1_) + ZZHT do L = L1_, 1, - 1 COLO2 = COLO2 + DDJ (L) * 0.20948d0 COLO3 = COLO3 + ZZJ (L) DELZ = ZTOP - ZHL (L) ZTOP = ZHL (L) ZKM = ZHL (L) * 1.d-5 ! write(6,'(1x,i3,0p,f6.2,f10.3,f7.2,1p,4e9.2,0p,f10.4,2f8.5,2i3)') & ! & L,ZKM,PPJ(L),TTJ(L),DDJ(L)/DELZ,ZZJ(L)/DELZ, & ! & COLO2,COLO3,DTAUX(L),POMEGAX(1,L),POMEGAX(2,L)/3.d0, & ! & JXTRA(L+L),JXTRA(L+L-1) enddo return end subroutine JP_ATM !###################################################################### subroutine RD_TJPL(rootProc,NJ1,NAMFIL) !----------------------------------------------------------------------- ! Read in wavelength bins, solar fluxes, Rayleigh parameters, temperature- ! dependent cross sections and Rayleigh/aerosol scattering phase functions ! with temperature dependences. Current data originates from JPL'97 !----------------------------------------------------------------------- ! ! NAMFIL Name of spectral data file (j2_spec.dat) >> j2 for fast-J2 ! NJ1 Channel number for reading data file ! 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 ! !----------------------------------------------------------------------- # include "parm_CTM_fastJX65.h" # include "parm_MIE_fastJX65.h" # include "cmn_JVdat_fastJX65.h" character (len=128), intent(in) :: NAMFIL logical, intent(in) :: rootProc INTEGER, intent(in) :: NJ1 integer :: I, J, K, IW, NK, NQQQ, NWWW IF (rootProc) PRINT*," Reading the file: ", TRIM(NAMFIL) do J=1,X_ do K=1,3 TQQ(K,J) = 0.d0 enddo enddo !-------------spectral data--------------------------------------------- !c open(NJ1, FILE=NAMFIL) !call AsciiOpenRead (NJ1, NAMFIL) open (NJ1, FILE = NAMFIL, status = 'old', form = 'formatted') read(NJ1,'(A)') TITLE0 read(NJ1,'(10X,14I5)') NJVAL,NWWW,NW1,NW2 if (NJVAL.gt.X_) then write(6,300) NJVAL,X_ STOP "Fast_JX 6.5" endif !------------NQQQ = no. additional J-values from X-sects (O2,O3P,O3D+NQQQ) ! begBian, 1/07/2012 !NQQQ = NJVAL-3 NQQQ = NJVAL-1 ! endBian, 1/07/2012 if (rootProc) then print *,'RD_TJPL NJVAL NQQQ =',NJVAL,NQQQ endif read(NJ1,102) (WL(IW),IW=1,NWWW) read(NJ1,102) (FL(IW),IW=1,NWWW) read(NJ1,102) (QRAYL(IW),IW=1,NWWW) read(NJ1,102) (QBC(IW),IW=1,NWWW) ! From Loiusse et al. [JGR, 1996] ! !---Read O2 X-sects, O3 X-sects, O3=>O(1D) quant yields (each at 3 temps) read (NJ1, 103) TITLEJ(1), TQQ (1, 1), (QO2 (IW, 1), IW = 1, NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, 1), (QO2 (IW, 2), IW = 1, NWWW) read (NJ1, 103) TITLEJ3, TQQ (3, 1), (QO2 (IW, 3), IW = 1, NWWW) read (NJ1, 103) TITLEJ(2), TQQ (1, 2), (QO3 (IW, 1), IW = 1, NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, 2), (QO3 (IW, 2), IW = 1, NWWW) read (NJ1, 103) TITLEJ3, TQQ (3, 2), (QO3 (IW, 3), IW = 1, NWWW) read (NJ1, 103) TITLEJ (3), TQQ (1, 3), (Q1D (IW, 1), IW = 1, NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, 3), (Q1D (IW, 2), IW = 1, NWWW) read (NJ1, 103) TITLEJ3, TQQ (3, 3), (Q1D (IW, 3), IW = 1, NWWW) ! do K=1,3 ! read(NJ1,103) TITLEJ(K,1),TQQ(K,1), (QO2(IW,K),IW=1,NWWW) ! enddo ! do K=1,3 ! read(NJ1,103) TITLEJ(K,2),TQQ(K,2), (QO3(IW,K),IW=1,NWWW) ! enddo ! do K=1,3 ! read(NJ1,103) TITLEJ(K,3),TQQ(K,3), (Q1D(IW,K),IW=1,NWWW) ! enddo ! !---Read remaining species: X-sections at 2 T's do J=1,NQQQ ! begBian, bugfix, 11/30/2011 ! read(NJ1,103) TITLEJ(J+3),TQQ(1,J+3),(QQQ(IW,1,J),IW=1,NWWW) ! read(NJ1,103) TITLEJ2,TQQ(2,J+3),(QQQ(IW,2,J),IW=1,NWWW) read(NJ1,103) TITLEJ(J+3),TQQ(1,J+3),(QQQ(IW,1,J+3),IW=1,NWWW) read(NJ1,103) TITLEJ2, TQQ(2,J+3),(QQQ(IW,2,J+3),IW=1,NWWW) ! endBian !if(TRIM(titlej(j+3)).eq."C3H6Ot") ncat_acetone_trop = J+3 ! Added for GMI to catch trop acetone {AOO, 8/04} !if(TRIM(titlej(j+3)).eq."C3H6Os") ncat_acetone_stra = J+3 ! Added for GMI to catch strat acetone {AOO, 8/04} if(TRIM(titlej(j+3)).eq."MeVK") ncat_met_vinyl_ketone = J+3 ! Added for GMI to catch methylvinylketone {AOO, 8/04} if(TRIM(titlej(j+3)).eq."MEKeto") ncat_met_ethyl_ketone = J+3 ! Added for GMI to catch methylethylketone {AOO, 8/04} if(TRIM(titlej(j+3)).eq."MGlyxl") ncat_methyl_glyoxal = J+3 ! Added for GMI to catch methy glyoxal {AOO, 8/04} enddo read(NJ1,'(A)') TITLE0 ! begBian, 11/30/2011 read (NJ1, '(5x,i5,2f10.5)') JTAUMX, ATAU, ATAU0 !write (6, '(a,2f9.5,i5)') ' ATAU/ATAU0/JMX', ATAU, ATAU0, JTAUMX ! endBian, 11/30/2011 ! ! From Fast JX !---Read aerosol phase functions: read(NJ1,'(A10,I5,/)') TITLE0,NAA if(NAA.gt.A_) then ! write(6,350) NAA STOP " too many scat-data for Fast_JX" endif ! begBian, 11/30/2011 ! note: change NK when using new table NK=4 ! Fix number of wavelengths at 4 ! endBian do j=1,NAA read(NJ1,110) TITLAA(j) do k=1,NK read(NJ1,*) WAA(k,j),QAA(k,j),RAA(k,j),SAA(k,j), & & (PAA(i,k,j),i=1,8) enddo enddo ! ! !! ----------- !! FROM RD_MIE !! ----------- ! read (NJ1, '(i2,a78)') NAA, TITLE0 ! if (NAA.gt.A_) then ! write (6, * ) ' too many scat-data sets:', NAA, A_ ! stop ! endif ! read (NJ1, '(5x,i5,2f10.5)') JTAUMX, ATAU, ATAU0 ! write (6, '(a,2f9.5,i5)') ' ATAU/ATAU0/JMX', ATAU, ATAU0, JTAUMX ! read (NJ1, * ) ! do J = 1, NAA ! read (NJ1, '(3x,a20,32x,f5.3,15x,f5.3)') & ! TITLAA (J) , RAA (J) , DAA (J) ! ! ver 6.0 extend to 5 ref wavelengths for mie-sca ! do K = 1, 5 ! read (NJ1, '(f4.0,f7.4,f7.4,7f6.3,1x,f7.3,f8.4)') & ! WAA (K, J) , QAA (K, J), SAA (K, J), & ! (PAA (I, K, J) , I = 2, 8) ! PAA (1, K, J) = 1.d0 ! enddo ! enddo !----------------------------- ! Begin Tropospheric Treatment !----------------------------- !---truncate number of wavelengths to do troposphere-only if (W_.ne.WX_) then !---TROP-ONLY if (W_.eq.12) then IF (rootProc) write (6,'(a)') & '>TROP-ONLY reduce wavelengths to 12, drop strat X-sects' NW2 = 12 do IW = 1, 4 WL (IW) = WL (IW + 4) FL (IW) = FL (IW + 4) QRAYL (IW) = QRAYL (IW + 4) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 4, K) QO3 (IW, K) = QO3 (IW + 4, K) Q1D (IW, K) = Q1D (IW + 4, K) enddo ! begBian, bugfix, 11/30/2011 ! do J = 4, NQQQ do J = 4, NJVAL ! endBian QQQ (IW, 1, J) = QQQ (IW + 4, 1, J) QQQ (IW, 2, J) = QQQ (IW + 4, 2, J) enddo enddo do IW = 5, 12 WL (IW) = WL (IW + 6) FL (IW) = FL (IW + 6) QRAYL (IW) = QRAYL (IW + 6) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 6, K) QO3 (IW, K) = QO3 (IW + 6, K) Q1D (IW, K) = Q1D (IW + 6, K) enddo ! begBian, bugfix, 11/30/2011 ! do J = 4, NQQQ do J = 4, NJVAL ! endBian QQQ (IW, 1, J) = QQQ (IW + 6, 1, J) QQQ (IW, 2, J) = QQQ (IW + 6, 2, J) enddo enddo !---TROP-QUICK (must scale solar flux for W=5) elseif (W_.eq.8) then IF (rootProc) write (6, '(a)') & '>TROP-QUICK reduce wavelengths to 8, drop strat X-sects' NW2 = 8 do IW = 1, 1 WL (IW) = WL (IW + 4) FL (IW) = FL (IW + 4) * 2.d0 QRAYL (IW) = QRAYL (IW + 4) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 4, K) QO3 (IW, K) = QO3 (IW + 4, K) Q1D (IW, K) = Q1D (IW + 4, K) enddo ! begBian, bugfix, 11/30/2011 ! do J = 4, NQQQ do J = 4, NJVAL ! endBian QQQ (IW, 1, J) = QQQ (IW + 4, 1, J) QQQ (IW, 2, J) = QQQ (IW + 4, 2, J) enddo enddo do IW = 2, 8 WL (IW) = WL (IW + 10) FL (IW) = FL (IW + 10) QRAYL (IW) = QRAYL (IW + 10) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 10, K) QO3 (IW, K) = QO3 (IW + 10, K) Q1D (IW, K) = Q1D (IW + 10, K) enddo ! begBian, bugfix, 11/30/2011 ! do J = 4, NQQQ do J = 4, NJVAL ! endBian QQQ (IW, 1, J) = QQQ (IW + 10, 1, J) QQQ (IW, 2, J) = QQQ (IW + 10, 2, J) enddo enddo else IF (rootProc) write(6,*)'Number of used wavelengths wrong:', W_ stop endif endif !----------------------------- ! End Tropospheric Treatment !----------------------------- ! !!---Zero index arrays ! do j=1,JVN_ ! jind(j)=0 ! enddo ! !! !!---Set mapping index ! do j=1,NJVAL ! do k=1,JVN_ ! if(TRIM(JLABEL(k)) .eq. TRIM(TITLEJ(j))) jind(k)=j ! enddo ! enddo ! ! do k=1,JVN_ ! if(JFACTA(k).eq.0.d0 .AND. rootProc) & ! & write(6,*) 'RD_TJPL: Not using photolysis reaction ',k ! if(jind(k).eq.0) then ! if(JFACTA(k).eq.0.d0) then ! jind(k)=1 ! else ! write(6,*) 'RD_TJPL: Which J-rate for photolysis reaction ',k,' ?' ! STOP "Fast_JX" ! endif ! endif ! enddo ! !-------- ! Reset the titles for NJVAL-1 & NJVAL to be the two acetone J_s ! 61: C3H6O = Acet-a (CH3CO + CH3) ! 62: Q2-Ac = Acet-b (CH3 + CO + CH3) TITLEJ (NJVAL - 1) = 'Acet-a' TITLEJ (NJVAL) = 'Acet-b' 101 FORMAT(8E10.3) 102 FORMAT((10X,6E10.3)/(10X,6E10.3)/(10X,6E10.3)) 103 FORMAT(A7,F3.0,6E10.3/(10X,6E10.3)/(10X,6E10.3)) ! 103 FORMAT(A7,F3.0,7E10.3/(10X,7E10.3)) 104 FORMAT(13x,i2) 105 FORMAT(A7,3x,7E10.3) 110 format(3x,a20) 200 format(1x,' x-sect:',a10,3(3x,f6.2)) 201 format(1x,' pr.dep:',a10,7(1pE10.3)) 300 format(' Number of x-sections supplied to Fast-J2: ',i3,/, & & ' Maximum number allowed (NS) only set to: ',i3, & & ' - increase in fast_JX_jv_cmn.h') 350 format(' Too many phase functions supplied; increase NP to ',i2) 400 format(' Not using Herzberg bin') close(NJ1) return end subroutine RD_TJPL !###################################################################### !----------------------------------------------------------------------- subroutine RD_XXX (NJ1, NAMFIL) !----------------------------------------------------------------------- ! Read in wavelength bins, solar fluxes, Rayleigh, T-dep X-sections. ! !>>>>NEW v-6.4 changed to collapse wavelengths & x-sections to Trop-onl ! WX_ = 18 (parm_CTM.f) should match the JX_spec.dat wavelengt ! W_ = 12 (Trop-only) or 18 (std) is set in (parm_MIE.f). ! if W_=12 then drop strat wavels, and drop x-sects (e.g. N2O, ... ! W_ = 8, reverts to quick fix: fast-J (12-18) plus bin (5) s ! !----------------------------------------------------------------------- ! NAMFIL Name of spectral data file (JX_spec.dat) >> j2 for fast-J ! NJ1 Channel number for reading data file ! ! NJVAL Number of species to calculate J-values for ! NWWW Number of wavelength bins, from 1:NWWW ! 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) ! 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) !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" # include "parm_MIE_fastJX65.h" # include "cmn_JVdat_fastJX65.h" integer , intent (in) ::NJ1 character ( * ), intent (in) ::NAMFIL integer :: I, J, JJ, K, IW, NQQQ, NWWW, NQRD TQQ (:, :) = 0.d0 !----------spectral data----set for new format data J-ver8.3------------ ! note that NJVAL = # J-values, but NQQQ (>NJVAL) = # Xsects rea ! for 2005a data, NJVAL = 62 (including a spare XXXX) and ! NQQQ = 64 so that 4 wavelength datasets read in for aceto ! note NQQQ is not used outside this subroutine! ! >>>> W_ = 12 <<<< means trop-only, discard WL #1-4 and #9-10, some X-s open (NJ1, FILE = NAMFIL, status = 'old', form = 'formatted') read (NJ1, 100) TITLE0 read (NJ1, 101) NJVAL, NQRD, NWWW NW1 = 1 NW2 = NWWW if (NJVAL.gt.X_.or.NQRD.gt.X_) then write (6, 201) NJVAL, X_ stop endif write (6, '(1X,A)') TITLE0 !----J-values: 1=O2, 2=O3P,3=O3D 4=readin Xsects read (NJ1, 102) (WL (IW), IW = 1, NWWW) read (NJ1, 102) (FL (IW), IW = 1, NWWW) read (NJ1, 102) (QRAYL (IW), IW = 1, NWWW) !---Read O2 X-sects, O3 X-sects, O3=>O(1D) quant yields (each at 3 temps read (NJ1, 103) TITLEJ (1), TQQ (1, 1), (QO2 (IW, 1), IW = 1, & NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, 1), (QO2 (IW, 2), IW = 1, NWWW) read (NJ1, 103) TITLEJ3, TQQ (3, 1), (QO2 (IW, 3), IW = 1, NWWW) read (NJ1, 103) TITLEJ (2), TQQ (1, 2), (QO3 (IW, 1), IW = 1, & NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, 2), (QO3 (IW, 2), IW = 1, NWWW) read (NJ1, 103) TITLEJ3, TQQ (3, 2), (QO3 (IW, 3), IW = 1, NWWW) read (NJ1, 103) TITLEJ (3), TQQ (1, 3), (Q1D (IW, 1), IW = 1, & NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, 3), (Q1D (IW, 2), IW = 1, NWWW) read (NJ1, 103) TITLEJ3, TQQ (3, 3), (Q1D (IW, 3), IW = 1, NWWW) do J = 1, 3 write (6, 200) J, TITLEJ (J), (TQQ (I, J), I = 1, 3) enddo !---Read remaining species: X-sections at 2 T_s JJ = 4 do J = 4, NQRD read (NJ1, 103) TITLEJ (JJ), TQQ (1, JJ), (QQQ (IW, 1, JJ), & IW = 1, NWWW) read (NJ1, 103) TITLEJ2, TQQ (2, JJ), (QQQ (IW, 2, JJ), IW = 1, & NWWW) if (W_.eq.18.or.TITLEJ2 (7:7) .ne.'x') then !---include stratospheric J's (this also includes Cl and Br compounds!) write (6, 200) JJ, TITLEJ (JJ), (TQQ (I, JJ), I = 1, 2) JJ = JJ + 1 endif enddo NQQQ = JJ - 1 NJVAL = NJVAL + (NQQQ - NQRD) !---truncate number of wavelengths to do troposphere-only if (W_.ne.WX_) then !---TROP-ONLY if (W_.eq.12) then write (6, '(a)') ' >>>TROP-ONLY reduce wavelengths to 12, drop str & &at X-sects' NW2 = 12 do IW = 1, 4 WL (IW) = WL (IW + 4) FL (IW) = FL (IW + 4) QRAYL (IW) = QRAYL (IW + 4) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 4, K) QO3 (IW, K) = QO3 (IW + 4, K) Q1D (IW, K) = Q1D (IW + 4, K) enddo do J = 4, NQQQ QQQ (IW, 1, J) = QQQ (IW + 4, 1, J) QQQ (IW, 2, J) = QQQ (IW + 4, 2, J) enddo enddo do IW = 5, 12 WL (IW) = WL (IW + 6) FL (IW) = FL (IW + 6) QRAYL (IW) = QRAYL (IW + 6) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 6, K) QO3 (IW, K) = QO3 (IW + 6, K) Q1D (IW, K) = Q1D (IW + 6, K) enddo do J = 4, NQQQ QQQ (IW, 1, J) = QQQ (IW + 6, 1, J) QQQ (IW, 2, J) = QQQ (IW + 6, 2, J) enddo enddo !---TROP-QUICK (must scale solar flux for W=5) elseif (W_.eq.8) then write (6, '(a)') ' >>>TROP-QUICK reduce wavelengths to 8, drop str & &at X-sects' NW2 = 8 do IW = 1, 1 WL (IW) = WL (IW + 4) FL (IW) = FL (IW + 4) * 2.d0 QRAYL (IW) = QRAYL (IW + 4) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 4, K) QO3 (IW, K) = QO3 (IW + 4, K) Q1D (IW, K) = Q1D (IW + 4, K) enddo do J = 4, NQQQ QQQ (IW, 1, J) = QQQ (IW + 4, 1, J) QQQ (IW, 2, J) = QQQ (IW + 4, 2, J) enddo enddo do IW = 2, 8 WL (IW) = WL (IW + 10) FL (IW) = FL (IW + 10) QRAYL (IW) = QRAYL (IW + 10) do K = 1, 3 QO2 (IW, K) = QO2 (IW + 10, K) QO3 (IW, K) = QO3 (IW + 10, K) Q1D (IW, K) = Q1D (IW + 10, K) enddo do J = 4, NQQQ QQQ (IW, 1, J) = QQQ (IW + 10, 1, J) QQQ (IW, 2, J) = QQQ (IW + 10, 2, J) enddo enddo else write (6, * ) ' number of used wavelengths wrong:', W_ stop endif endif ! Reset the titles for NJVAL-1 & NJVAL to be the two acetone J_s ! 61: C3H6O = Acet-a (CH3CO + CH3) ! 62: Q2-Ac = Acet-b (CH3 + CO + CH3) TITLEJ (NJVAL - 1) = 'Acet-a' TITLEJ (NJVAL) = 'Acet-b' close (NJ1) 100 format(a) 101 format(10x,5i5) 102 format(10x, 6e10.3/(10x,6e10.3)/(10x,6e10.3)) 103 format(a7,f3.0,6e10.3/(10x,6e10.3)/(10x,6e10.3)) 200 format(1x,' x-sect:',i3,a10,3(3x,f6.2)) 201 format(' Number of x-sections supplied to Fast-J2: ',i3,/, & & ' Maximum number allowed (X_) only set to: ',i3, & & ' - increase in cmn_jv.f') return end subroutine RD_XXX !###################################################################### !----------------------------------------------------------------------- !subroutine RD_MIE (NJ1, NAMFIL) !----------------------------------------------------------------------- !-------aerosols/cloud scattering data set for fast-JX (ver 5.3+) ! >>>>>>>>>>>>>>>>spectral data rev to J-ref ver8.5 (5/05)<<<<<<<<<<<< !----------------------------------------------------------------------- ! NAMFIL Name of scattering data file (e.g., FJX_scat.dat) ! NJ1 Channel number for reading data file ! NAA Number of categories for scattering phase functions ! QAA Aerosol scattering phase functions ! NK Number of wavelengths at which functions supplied (set as ! WAA Wavelengths for the NK supplied phase functions ! PAA Phase function: first 8 terms of expansion ! RAA Effective radius associated with aerosol type ! SAA Single scattering albedo !----------------------------------------------------------------------- !implicit none !# include "parm_CTM_fastJX65.h" !# include "parm_MIE_fastJX65.h" ! !# include "cmn_JVdat_fastJX65.h" !integer , intent (in) ::NJ1 !character ( * ), intent (in) ::NAMFIL ! !integer :: I, J, K ! !open (NJ1, FILE = NAMFIL, status = 'old', form = 'formatted') !read (NJ1, '(i2,a78)') NAA, TITLE0 !if (NAA.gt.A_) then ! write (6, * ) ' too many scat-data sets:', NAA, A_ ! stop !endif !read (NJ1, '(5x,i5,2f10.5)') JTAUMX, ATAU, ATAU0 !write (6, '(a,2f9.5,i5)') ' ATAU/ATAU0/JMX', ATAU, ATAU0, JTAUMX !read (NJ1, * ) !do J = 1, NAA !read (NJ1, '(3x,a20,32x,f5.3,15x,f5.3)') TITLAA (J) , RAA (J) , & ! DAA (J) ! ! ver 6.0 extend to 5 ref wavelengths for mie-sca !do K = 1, 5 !read (NJ1, '(f4.0,f7.4,f7.4,7f6.3,1x,f7.3,f8.4)') WAA (K, J) , & ! QAA (K, J) , SAA (K, J) , (PAA (I, K, J) , I = 2, 8) !PAA (1, K, J) = 1.d0 !enddo ! !enddo ! !close (NJ1) !write (6, '(a,9f8.1)') ' Aerosol optical: r-eff/rho/Q(@wavel):', & ! (WAA (K, 1) , K = 1, 5) !write (6, * ) TITLE0 !do J = 1, NAA !write (6, '(i3,1x,a8,7f8.3)') J, TITLAA (J) , RAA (J) , DAA (J) , & ! (QAA (K, J) , K = 1, 5) ! !enddo !return ! ! !end subroutine RD_MIE !###################################################################### !----------------------------------------------------------------------- !subroutine RD_UM (NJ1, NAMFIL) !----------------------------------------------------------------------- !-------UMich aerosol optical data for fast-JX (ver 6.1+) !----------------------------------------------------------------------- ! NAMFIL Name of scattering data file (e.g., FJX_scat.dat) ! NJ1 Channel number for reading data file !----------------------------------------------------------------------- !implicit none !# include "parm_CTM_fastJX65.h" !# include "parm_MIE_fastJX65.h" ! !# include "cmn_JVdat_fastJX65.h" !integer , intent (in) ::NJ1 ! !character ( * ), intent (in) ::NAMFIL ! !integer :: I, J, K, L ! !open (NJ1, FILE = NAMFIL, status = 'old', form = 'formatted') !read (NJ1, '(a78)') TITLE0 !write (6, * ) 'UMichigan Aerosol optical data' ! !write (6, * ) TITLE0 !---33 Different UM Aerosol Types: SULF, SS-1,-2,-3,-4, DD-1,-2,-3,-4, !--- FF00(0%BC), FF02, ...FF14(14%BC), BB00, BB02, ...BB30(30%BC) !do L = 1, 33 !read (NJ1, '(a4)') TITLUM (L) !---21 Rel Hum: K=1=0%, =2=5%, ... =20=95%, =21=99% !do K = 1, 21 !---6 wavelengths: J=1=200nm, 2=300nm, 3=400nm, (4'=550nm) 4=600nm, 5=10 !---3 optic vars: I=1=SSAlbedo, =2=g, =3=k-ext !read (NJ1, '(9f9.5,27x,6f9.5)') ( (UMAER (I, J, K, L) , I = 1, 3) & ! , J = 1, 5) !enddo ! !enddo ! !close (NJ1) ! !write (6, '(7(i5,1x,a4))') (L, TITLUM (L) , L = 1, 33) !return ! ! !end subroutine RD_UM !###################################################################### !----------------------------------------------------------------------- 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 SOLARZ (GMTIME, NDAY, YGRDJ, XGRDI, SZA, COSSZA, SOLFX) !----------------------------------------------------------------------- ! GMTIME = UT for when J-values are wanted ! (for implicit solver this is at the end of the time step) ! NDAY = integer day of the year (used for solar lat and declin) ! YGRDJ = laitude (radians) for grid (I,J) ! XGDRI = longitude (radians) for grid (I,J) ! ! SZA = solar zenith angle in degrees ! COSSZA = U0 = cos(SZA) !----------------------------------------------------------------------- implicit none real (8) , intent (in) ::GMTIME, YGRDJ, XGRDI integer , intent (in) ::NDAY real (8) , intent (in) :: SZA real (8) , intent (out) :: COSSZA, SOLFX real (8) :: PI, PI180, LOCT real (8) :: SINDEC, SOLDEK, COSDEC, SINLAT, SOLLAT, COSLAT, COSZ ! PI = 3.141592653589793d0 PI180 = PI / 180.d0 !SINDEC = 0.3978d0 * sin (0.9863d0 * (dble (NDAY) - 80.d0) * PI180) !SOLDEK = asin (SINDEC) !COSDEC = cos (SOLDEK) !SINLAT = sin (YGRDJ) !SOLLAT = asin (SINLAT) !COSLAT = cos (SOLLAT) ! !LOCT = ( ( (GMTIME) * 15.d0) - 180.d0) * PI180 + XGRDI !COSSZA = COSDEC * COSLAT * cos (LOCT) + SINDEC * SINLAT !SZA = acos (COSSZA) / PI180 COSSZA = COS(SZA*PI180) !write (6, * ) ' XGRDI,YGRDJ', XGRDI, YGRDJ !write (6, * ) ' LOCT (rad)', LOCT !write (6, * ) ' SINDEC,COSDEC', SINDEC, COSDEC !write (6, * ) ' SINLAT,COSLAT', SINLAT, COSLAT !write (6, * ) ' COS, SZA', COSSZA, SZA SOLFX = 1.d0 - (0.034d0 * cos (dble (NDAY - 186) * 2.d0 * PI / 365.d0) ) return end subroutine SOLARZ !###################################################################### !----------------------------------------------------------------------- subroutine SPHERE2 (GMU, RAD, ZHL, ZZHT, AMF2, L1_) !----------------------------------------------------------------------- !----new v6.2: does AirMassFactors for mid-layer, needed for SZA ~ 90 ! This new AMF2 does each of the half-layers of the CTM separately, ! whereas the original, based on the pratmo code did the whole layer ! and thus calculated the ray-path to the CTM layre edges, NOT the m ! Since fast-JX is meant to calculate the intensity at the mid-layer, t ! solar beam at low sun (interpolated between layer edges) was incor ! This new model does make some approximations of the geometry of the l ! the CTM layer is split evenly in mass (good) and in height (approx ! ! 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). !----------------------------------------------------------------------- ! in: ! GMU = MU0 = cos(solar zenith angle) ! RAD radius of Earth mean sea level (cm) ! ZHL(L) height (cm) of the bottome edge of CTM level L ! ZZHT scale height (cm) used above top of CTM (ZHL(L_+1) ! L1_ dimension of CTM = levels +1 (L+1 = above-CTM level) ! out: ! AMF2(I,J) = air mass factor for CTM level I for sunlight reaching !----------------------------------------------------------------------- implicit none integer , intent (in) ::L1_ real (8) , intent (in) ::GMU, RAD, ZHL (L1_ + 1), ZZHT real (8) , intent (out) ::AMF2 (2 * L1_ + 1, 2 * L1_ + 1) ! RZ Distance from centre of Earth to each point (cm) ! RQ Square of radius ratios ! SHADHT Shadow height for the current SZA ! XL Slant path between points integer :: I, J, K, II, L2 real (8) :: XMU1, XMU2, XL, DIFF, SHADHT, RZ (L1_ + 1) real (8) :: RZ2 (2 * L1_ + 1), RQ2 (2 * L1_ + 1) ! !--- must have top-of-atmos (NOT top-of-CTM) defined ! ZHL(L1_+1) = ZHL(L1_) + ZZHT RZ (1) = RAD+ZHL (1) do II = 2, L1_ + 1 RZ (II) = RAD+ZHL (II) enddo !---calculate heights for edges of split CTM-layers L2 = 2 * L1_ do II = 2, L2, 2 I = II / 2 RZ2 (II - 1) = RZ (I) RZ2 (II) = 0.5d0 * (RZ (I) + RZ (I + 1) ) enddo RZ2 (L2 + 1) = RZ (L1_ + 1) do II = 1, L2 RQ2 (II) = (RZ2 (II) / RZ2 (II + 1) ) **2 enddo !---shadow height for SZA > 90 if (GMU.lt.0.0d0) then SHADHT = RZ2 (1) / dsqrt (1.0d0 - GMU**2) else SHADHT = 0.d0 endif !---up from the surface calculating the slant paths between each level !--- and the level above, and deriving the appropriate Air Mass Factor AMF2 (:, :) = 0.d0 do 16 J = 1, 2 * L1_ + 1 ! Air Mass Factors all zero if below the tangent height if (RZ2 (J) .lt.SHADHT) goto 16 ! Ascend from layer J calculating AMF2s XMU1 = abs (GMU) do I = J, 2 * L1_ XMU2 = dsqrt (1.0d0 - RQ2 (I) * (1.0d0 - XMU1**2) ) XL = RZ2 (I + 1) * XMU2 - RZ2 (I) * XMU1 AMF2 (I, J) = XL / (RZ2 (I + 1) - RZ2 (I) ) XMU1 = XMU2 enddo !--fix above top-of-atmos (L=L1_+1), must set DTAU(L1_+1)=0 AMF2 (2 * L1_ + 1, J) = 1.d0 ! ! Twilight case - Emergent Beam, calc air mass factors below layer if (GMU.ge.0.0d0) goto 16 ! Descend from layer J XMU1 = abs (GMU) do II = J - 1, 1, - 1 DIFF = RZ2 (II + 1) * sqrt (1.0d0 - XMU1**2) - RZ2 (II) ! filter if (II.eq.1) DIFF = max (DIFF, 0.d0) ! Tangent height below current level - beam passes through twice if (DIFF.lt.0.0d0) then XMU2 = sqrt (1.0d0 - (1.0d0 - XMU1**2) / RQ2 (II) ) XL = abs (RZ2 (II + 1) * XMU1 - RZ2 (II) * XMU2) AMF2 (II, J) = 2.d0 * XL / (RZ2 (II + 1) - RZ2 (II) ) XMU1 = XMU2 ! Lowest level intersected by emergent beam else XL = RZ2 (II + 1) * XMU1 * 2.0d0 AMF2 (II, J) = XL / (RZ2 (II + 1) - RZ2 (II) ) goto 16 endif enddo 16 end do return end subroutine SPHERE2 !###################################################################### !----------------------------------------------------------------------- subroutine EXTRAL (DTAUX, L1X, L2X, NX, JTAUMX, ATAU, ATAU0, & JXTRA) !----------------------------------------------------------------------- ! ! new version 6.1, add sub-layers (JXTRA) to thick cloud/aerosol laye ! this version sets up log-spaced sub-layers of increasing thickness ! ! DTAUX(L=1:L1X) = Optical Depth in layer L (generally 600 nm OD) ! This can be just cloud or cloud+aerosol, it is used only to set ! the number in levels to insert in each layer L ! Set for log-spacing of tau levels, increasing top-down. ! ! N.B. the TTAU, etc calculated here are NOT used elsewhere !---The log-spacing parameters have been tested for convergence and chos !--- to be within 0.5% for ranges OD=1-500, rflect=0-100%, mu0=0.1-1.0 !--- use of ATAU = 1.18 and min = 0.01, gives at most +135 pts for OD=1 !--- ATAU = 1.12 now recommended for more -accurate heating rates (not !----------------------------------------------------------------------- ! implicit none !index of cloud/aerosol integer , intent (in) ::JTAUMX, L1X, L2X !Mie scattering array size integer , intent (in) ::NX !cloud+3aerosol OD in each real (8) , intent (in) ::DTAUX (L1X) real (8) , intent (in) ::ATAU, ATAU0 !number of sub-layers to b integer , intent (out) ::JXTRA (L2X + 1) ! integer :: JTOTL, I, L, L2 real (8) :: TTAU (L2X + 1), DTAUJ, ATAU1, ATAULN, ATAUM, ATAUN1 ! !---Reinitialize arrays TTAU (:) = 0.d0 JXTRA (:) = 0 ! !---combine these edge- and mid-layer points into grid of size: !--- L2X+1 = 2*L1X+1 = 2*L_+3 !---calculate column optical depths above each level, TTAU(1:L2X+1) !--- note that TTAU(L2X+1)=0 and TTAU(1)=total OD ! !---Divide thick layers to achieve better accuracy in the scattering cod !---In the original fast-J, equal sub-layers were chosen, this is wastef !---and this new code (ver 5.3) uses log-scale: !--- Each succesive layer (down) increase thickness by ATAU > 1 !--- e.g., if ATAU = 2, a layer with OD = 15 could be divided int !--- 4 sub-layers with ODs = 1 - 2 - 4 - 8 !---The key parameters are: !--- ATAU = factor increase from one layer to the next !--- ATAUMN = the smallest OD layer desired !--- JTAUMX = maximum number of divisions (i.e., may not get to A !---These are hardwired below, can be changed, but have been tested/opti ATAU1 = ATAU - 1.d0 ATAULN = log (ATAU) TTAU (L2X + 1) = 0.0d0 do L2 = L2X, 1, - 1 L = (L2 + 1) / 2 DTAUJ = 0.5d0 * DTAUX (L) TTAU (L2) = TTAU (L2 + 1) + DTAUJ !---Now compute the number of log-spaced sub-layers to be added in !--- the interval TTAU(L2) > TTAU(L2+1) !---The objective is to have successive TAU-layers increasing by factor !---the number of sub-layers + 1 if (TTAU (L2) .lt.ATAU0) then JXTRA (L2) = 0 else ATAUM = max (ATAU0, TTAU (L2 + 1) ) ATAUN1 = log (TTAU (L2) / ATAUM) / ATAULN JXTRA (L2) = min (JTAUMX, max (0, int (ATAUN1 - 0.5d0) ) ) endif enddo !---check on overflow of arrays, cut off JXTRA at lower L if too many le JTOTL = L2X + 2 do L2 = L2X, 1, - 1 JTOTL = JTOTL + JXTRA (L2) if (JTOTL.gt.NX / 2) then ! write (6, '(A,2I5,F9.2)') 'N_/L2_/L2-cutoff JXTRA:', NX, L2X, L2 do L = L2, 1, - 1 JXTRA (L) = 0 enddo goto 10 endif enddo 10 continue return end subroutine EXTRAL !###################################################################### !----------------------------------------------------------------------- subroutine OPMIE (DTAUX, POMEGAX, U0, RFL, AMF2, JXTRA, FJACT, & FJTOP, FJBOT, FSBOT, FJFLX, FLXD, FLXD0) !----------------------------------------------------------------------- implicit none # include "parm_CTM_fastJX65.h" # include "parm_MIE_fastJX65.h" # include "cmn_JVdat_fastJX65.h" real (8) , intent (in) ::DTAUX (L1_, W_), POMEGAX (8, L1_, W_) real (8) , intent (in) ::AMF2 (2 * L1_ + 1, 2 * L1_ + 1) real (8) , intent (in) ::U0, RFL (W_) integer , intent (in) ::JXTRA (L2_ + 1) real (8) , intent (out) ::FJACT (L_, W_), FJTOP (W_), FJBOT (W_), FSBOT (W_) real (8) , intent (out) ::FJFLX (L_, W_), FLXD (L1_, W_), FLXD0 (W_) integer :: JNDLEV (L_), JNELEV (L1_) integer :: JADDLV (L2_ + 1), JADDTO (L2_ + 1), L2LEV (L2_ + 1) integer :: JTOTL, I, II, J, K, L, LL, IX, JK, L2, L2L, L22, LZ, LZZ, ND integer :: LZ0, LZ1, LZMID real (8) :: SUMT, SUMJ real (8) :: DTAU (L1_ + 1, W_), POMEGAJ (M2_, L2_ + 1, W_), TTAU (L2_ + 1, W_) real (8) :: FTAU2 (L2_ + 1, W_), POMEGAB (M2_, W_) real (8) :: ATAUA, ATAUZ, XLTAU, TAUDN, TAUUP, DTAUJ, FJFLX0 real (8) , dimension (W_) ::TAUBTM, TAUTOP, FBTM, FTOP, ZFLUX !--- variables used in mie code----------------------------------------- real (8) , dimension (W_) ::FJT, FJB real (8) , dimension (N_, W_) ::FJ, FZ, ZTAU real (8) , dimension (M2_, N_, W_) ::POMEGA real (8) , dimension (2 * L1_, W_) ::FLXD2 ! fast-J Mie code for J_s, only uses 8-term expansion, 4-Gauss pts ! ! in: ! DTAUX(1:L1_,1:W_) = optical depth of each layer ! POMEGAX(1:8,1:L1_,1:W_) = scattering phase fn (multiplied by s-s a ! U0 = cos (SZA) ! RFL(1:W_) = Lambertian albedo of surface ! AMF2(1:2*L1_+1,1:2*L1_+1) = air mass factor (I,L)=wt of layer-I to ! AMF2 now does both edges and middle of CTM layers ! JXTRA(1:L1_) = number 0:J = no. of additional levels to be inserte ! out: ! FJACT(1:L_,1:W_) = mean actinic flux(diff+direct) at std CTM level ! (new ver 5.7 diagnostics for fluxes, deposition) fluxes 'down' are < ! FJTOP(1:W_) = diffuse flux out top-of-atmosphere (TAU=0 above top ! FJBOT(1:W_) = diffuse flux onto surface (<0 by definition) ! FSBOT(1:W_) = direct/solar flux onto surface (<0 by definition) ! FJFLX(1:L_,1:W_) = diffuse flux across top of model layer L ! this connects with FJBOT = FJFLX(0) & FJTOP = FJFLX(L_+1) (not ! FLXD(1:L_+1,1:W_) = solar flux deposited in layer L (includes lyr ! this should take into account sphericity, and is not just = mu0 ! FLXD0(1:W_) = sum of solar flux deposited in atmos ! does NOT include flux on lower surface, does NOT mean absorbed! !----------------------------------------------------------------------- ! ! DTAU Local optical depth of each CTM level ! TTAU Optical depth of air vertically above each point (to top ! FTAU2 Attenuation of solar beam ! POMEGAJ Scattering phase function ! !---new ver 5.3 code adds sub-layers (# = JXTRA(L2)) using ATAU as the ! factor increase from sub-layer to sub-layer ! !---------------------SET UP FOR MIE CODE------------------------------- !-----------------wavelength independent-------------------------------- ! ! 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') ! ! The following mapping holds for JADDLV=0 ! Surface: TTAU(1) ==> ZTAU(2*L2_+1) ! Top: TTAU(L2_) ==> ZTAU(3) ! Infinity: 0.0 ==> ZTAU(1) ! index: 2*(L2_+1-L2)+1 ==> LZ ! ! Mie scattering code only used from surface to level L2_ !----------------------------------------------------------------------- ! !----------------------------------------------------------------------- ! Insert new levels, working downwards from the top of the atmosphere ! to the surface (down in 'LZ', up in 'L2'). This allows ztau and pomeg ! to be incremented linearly, and the flux fz to be attenuated top-down ! (avoiding problems where lower level fluxes are zero). !----------------------------------------------------------------------- ! ! Ascend through atmosphere transposing grid and adding extra points ! remember L2=1 is surface of CTM, but last layer (LZ) in scattering co ! there are twice the number of layers in the LZ arrays (2*L2_ + 2*JADD ! because we need to insert the intermediate layers (even LZ) for the ! asymmetric scattering code. ! Transfer the L2=1:L2_+1 values (TTAU,FTAU2,POMEGAJ) onto the reverse ! order, expanded, doubled-level scatter grid. ! Note that we need to deal with the expansion by JADD levels (L2L). ! These JADDLV levels are skipped and need to be interpolated later ! Note that only odd LZ levels are filled, !----------------------re-grid data------------------------------------- ! Calculate cumulative total and define levels we want J-values at. ! Sum upwards for levels, and then downwards for Mie code readjustments ! ! JXTRA(L2) Number of new levels to add between (L2) and (L2+1) ! ***JXTRA(1:L2_+1) is calculated based on the aerosol+cloud O ! JADDLV(L2) Number of new levels actually added at each wavelength ! where JADDLV = 0 when there is effectively no FTAU2 ! JADDTO(L2) Total number of new levels to add to and above level ! JNDLEV(L) = L2 index that maps on CTM mid-layer L ! !---JADDLV(L2=1:L2_) = number of levels to add between TTAU2(L2) and TTA !--- JADDLV is taken from JXTRA, which is based on visible OD. !--- JADDTO(L2=1:L2_+1) is the cumulative number of levels to be adde !---these should be fixed for all wavelengths to lock-in the array sizes do L2 = 1, L2_, 1 JADDLV (L2) = JXTRA (L2) enddo JADDTO (L2_ + 1) = 0 do L2 = L2_, 1, - 1 JADDTO (L2) = JADDTO (L2 + 1) + JADDLV (L2) enddo !---expanded grid now included CTM edge and mid layers plus expanded !--- grid to allow for finer delta-tau at tops of clouds. !--- DIM of new grid = L2_ + JADDTO(1) + 1 !---L2LEV(L2) = L2-index for old level L2 in expanded J-grid (w/JADDLV) ! in absence of JADDLV, L2LEV(L2) = L2 L2LEV (1) = 1 do L2 = 2, L2_ + 1 L2LEV (L2) = L2LEV (L2 - 1) + 1 + JADDLV (L2 - 1) enddo !---JNDLEV(L=1:L_) = L2-index in expanded grid for CTM mid-layer L !---JNELEV(L=1:L_) = L2-index for top of layer L do L = 1, L_ JNDLEV (L) = L2LEV (2 * L) JNELEV (L) = L2LEV (2 * L + 1) enddo !need to set this to top-of-atmosphere JNELEV (L_ + 1) = 0 ND = 2 * L2_ + 2 * JADDTO (1) + 1 if (ND.gt.N_) then write (6, '(a,2i9)') ' overflow of scatter arrays:', ND, N_ stop endif !----------------begin wavelength dependent set up---------------------- !---Reinitialize arrays ZTAU (:, :) = 0.d0 FZ (:, :) = 0.d0 POMEGA (:, :, :) = 0.d0 do K = 1, W_ !---Set up optical depth DTAU(L) do L = 1, L1_ DTAU (L, K) = DTAUX (L, K) enddo DTAU (L1_ + 1, K) = 0.d0 !---Define the total scattering phase fn for each CTM layer L=1:L_+1 !--- from a DTAU-wt_d mix of aerosols, cloud & Rayleigh !---No. of quadrature pts fixed at 4(M_), expansion of phase fn @ 8 do L = 1, L1_ do I = 1, M2_ POMEGAJ (I, L, K) = POMEGAX (I, L, K) enddo enddo !---Calculate attenuated incident beam exp(-TTAU/U0 = DTAU * AirMassFact !--- at the middle & edges of the CTM layers L=1:2*L1_+1 !--- L1_ is top-edge of CTM (ie, L=38 = 2 hPa) which has TAU > 0 !--- note that DTAU(L1_) is optical depth in the FULL CTM layer just ab FTAU2 (:, :) = 0.d0 FTAU2 (L2_ + 1, :) = 1.0d0 do LL = 1, 2 * L1_ + 1 L = (LL + 1) / 2 if (AMF2 (LL, LL) .gt.0.0d0) then XLTAU = 0.0d0 do II = 1, 2 * L1_ + 1 I = (II + 1) / 2 XLTAU = XLTAU + 0.5d0 * DTAU (I, K) * AMF2 (II, LL) enddo ! zero out flux at 1e-33 if (XLTAU.lt.76.d0) then FTAU2 (LL, K) = exp ( - XLTAU) endif endif enddo !---calculate direct solar flux deposited in each CTM half-layer: L=1:L2 !--- use FSBOT for surface flux, cannot do layer above CTM (L_+1) FLXD2 (:, :) = 0.d0 do LL = 1, 2 * L1_ if (AMF2 (LL, LL) .gt.0.d0) then FLXD2 (LL, K) = (FTAU2 (LL + 1, K) - FTAU2 (LL, K) ) / AMF2 ( & LL, LL) endif enddo if (AMF2 (1, 1) .gt.0.d0) then FSBOT (K) = FTAU2 (1, K) / AMF2 (1, 1) else FSBOT (K) = 0.d0 endif do LL = 2, 2 * L1_, 2 L = LL / 2 FLXD (L, K) = FLXD2 (LL, K) + FLXD2 (LL - 1, K) enddo !---integrate solar flux depositied in CTM layers L=1:L_, cannot do top !--- note FLXD0 .ne. (1.d0 - FTAU(L_+1))/AMF(L_+1,L_+1) with spherical FLXD0 (K) = 0.d0 if (AMF2 (2 * L1_, 2 * L1_) .gt.0.d0) then do L = 1, L1_ FLXD0 (K) = FLXD0 (K) + FLXD (L, K) enddo 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. !----------------------------------------------------------------------- !---Now combine the CTM layer edges (1:L_+2) with the CTM mid-layer !--- points (1:L_) plus 1 for the mid point of added top layer. !---combine these edge- and mid-layer points into grid of size: !--- L2_+1 = 2*L1_+1 = 2*L_+3 !---calculate column optical depths above each level, TTAU(1:L2_+1) !--- note that TTAU(L2_+1)=0 and TTAU(1)=total OD TTAU (L2_ + 1, K) = 0.0d0 do L2 = L2_, 1, - 1 L = (L2 + 1) / 2 DTAUJ = 0.5d0 * DTAU (L, K) TTAU (L2, K) = TTAU (L2 + 1, K) + DTAUJ enddo !----solar flux incident on lower boundary & Lambertian reflect factor: if (FSBOT (K) .gt.0.d0) then ZFLUX (K) = FSBOT (K) * RFL (K) / (1.d0 + RFL (K) ) else ZFLUX (K) = 0.d0 endif ! Calculate scattering properties, level centres then level boundaries !>>>>>be careful of order, we are overwriting/shifting the 'POMEGAJ' upw do L2 = L2_, 2, - 2 L = L2 / 2 do I = 1, M2_ POMEGAJ (I, L2, K) = POMEGAJ (I, L, K) enddo enddo !---lower boundary value is set (POMEGAJ(I,1), but set upper: do I = 1, M2_ POMEGAJ (I, L2_ + 1, K) = POMEGAJ (I, L2_, K) enddo !---now have POMEGAJ filled at even points from L2=3:L2_-1 !---use inverse interpolation for correct tau-weighted values at edges do L2 = 3, L2_ - 1, 2 TAUDN = TTAU (L2 - 1, K) - TTAU (L2, K) TAUUP = TTAU (L2, K) - TTAU (L2 + 1, K) do I = 1, M2_ POMEGAJ (I, L2, K) = (POMEGAJ (I, L2 - 1, K) * TAUDN + POMEGAJ (I, & L2 + 1, K) * TAUUP) / (TAUDN + TAUUP) enddo enddo !---at this point FTAU2(1:L2_+1) and POMEAGJ(1:8, 1:L2_+1) !--- where FTAU2(L2_+1) = 1.0 = top-of-atmos, FTAU2(1) = surface ! L2 = index of CTM edge- and mid-layers do L2 = 1, L2_ + 1 ! L2L = index for L2 in expanded scale(JA L2L = L2LEV (L2) ! LZ = index for L2 in scatt arrays LZ = ND+2 - 2 * L2L ZTAU (LZ, K) = TTAU (L2, K) FZ (LZ, K) = FTAU2 (L2, K) do I = 1, M2_ POMEGA (I, LZ, K) = POMEGAJ (I, L2, K) enddo enddo ! Now go thru the pairs of L2 levels to see if we need JADD levels ! L2 = index of CTM edge- and mid-layer do L2 = 1, L2_ ! L2L = index for L2 in expanded scale( L2L = L2LEV (L2) ! LZ = index for L2 in scatt arrays LZ = ND+2 - 2 * L2L ! L22 = 0 if no added level L22 = L2LEV (L2 + 1) - L2LEV (L2) - 1 if (L22.gt.0) then TAUBTM (K) = TTAU (L2, K) TAUTOP (K) = TTAU (L2 + 1, K) FBTM (K) = FTAU2 (L2, K) FTOP (K) = FTAU2 (L2 + 1, K) do I = 1, M2_ POMEGAB (I, K) = POMEGAJ (I, L2, K) enddo !---to fit L22 new layers between TAUBOT > TAUTOP, calculate new 1/ATAU !--- such that TAU(just above TAU-btm) = ATUAZ * TAUBTM < TAUBTM ATAUZ = exp ( - log (TAUBTM (K) / max (TAUTOP (K), ATAU0) ) & / float (L22 + 1) ) ! add odd levels between L2LEV(L2) & L2L do L = 1, L22 ! LZZ = index(odd) of added level in scat LZZ = LZ - 2 * L ZTAU (LZZ, K) = TAUBTM (K) * ATAUZ !---fraction from TAUBTM=>TAUTOP ATAUA = (TAUBTM (K) - ZTAU (LZZ, K) ) / (TAUBTM (K) - TAUTOP & (K) ) !---solar flux at interp-levels: use exp(TAU/U0) if U0>0.02 (89 deg), !---else scale by TAU if (U0.gt.0.02d0) then FZ (LZZ, K) = FTOP (K) * exp ( (TAUTOP (K) - ZTAU (LZZ, & K) ) / U0) else if (FBTM (K) .lt.1.d-32) then FZ (LZZ, K) = 0.d0 else FZ (LZZ, K) = FBTM (K) * (FTOP (K) / FBTM (K) ) ** & ATAUA endif endif do I = 1, M2_ POMEGA (I, LZZ, K) = POMEGAB (I, K) + ATAUA * (POMEGAJ (I, & L2 + 1, K) - POMEGAB (I, K) ) enddo TAUBTM (K) = ZTAU (LZZ, K) FBTM (K) = FZ (LZZ, K) do I = 1, M2_ POMEGAB (I, K) = POMEGA (I, LZZ, K) enddo enddo endif enddo ! Now fill in the even points with simple interpolation in scatter arr do LZ = 2, ND-1, 2 ZTAU (LZ, K) = 0.5d0 * (ZTAU (LZ - 1, K) + ZTAU (LZ + 1, K) ) FZ (LZ, K) = sqrt (FZ (LZ - 1, K) * FZ (LZ + 1, K) ) do I = 1, M2_ POMEGA (I, LZ, K) = 0.5d0 * (POMEGA (I, LZ - 1, K) + POMEGA (I, & LZ + 1, K) ) enddo enddo ! wavelength loop! enddo !----------------------------------------------------------------------- call MIESCT (FJ, FJT, FJB, POMEGA, FZ, ZTAU, ZFLUX, RFL, U0, ND) !----------------------------------------------------------------------- !---Move mean intensity from scatter array FJ(LZ=1:ND) !--- to CTM mid-level array FJACT(L=1:L_) do K = 1, W_ !---mean intensity at mid-layer: 4* + solar ! do L = 1,L_ ! L2L = JNDLEV(L) ! LZ = ND+2 - 2*L2L ! FJACT(L,K) = 4.d0*FJ(LZ,K) + FZ(LZ,K) ! enddo !---mean intensity averaged throughout layer: do L = 1, L_ LZ0 = ND+2 - 2 * JNELEV (L) if (L.gt.1) then LZ1 = ND+2 - 2 * JNELEV (L - 1) else LZ1 = ND endif SUMJ = (4.d0 * FJ (LZ0, K) + FZ (LZ0, K) ) * (ZTAU (LZ0 + 2, K) & - ZTAU (LZ0, K) ) + (4.d0 * FJ (LZ1, K) + FZ (LZ1, K) ) * (ZTAU ( & LZ1, K) - ZTAU (LZ1 - 2, K) ) SUMT = ZTAU (LZ0 + 2, K) - ZTAU (LZ0, K) + ZTAU (LZ1, K) - ZTAU ( & LZ1 - 2, K) do LZ = LZ0 + 2, LZ1 - 2, 2 SUMJ = SUMJ + (4.d0 * FJ (LZ, K) + FZ (LZ, K) ) * (ZTAU (LZ + 2, & K) - ZTAU (LZ - 2, K) ) SUMT = SUMT + ZTAU (LZ + 2, K) - ZTAU (LZ - 2, K) enddo FJACT (L, K) = SUMJ / SUMT enddo !---mean diffuse flux: 4 (not solar) at top of layer L !--- average (tau-wtd) the h's just above and below the L-edge do L = 1, L_ L2L = JNELEV (L) LZ = ND+2 - 2 * L2L FJFLX0 = (ZTAU (LZ + 1, K) - ZTAU (LZ, K) ) / (ZTAU (LZ + 1, K) & - ZTAU (LZ - 1, K) ) FJFLX (L, K) = 4.d0 * (FJ (LZ - 1, K) * FJFLX0 + FJ (LZ + 1, K) & * (1.d0 - FJFLX0) ) enddo !---diffuse fluxes reflected at top, incident at bottom FJTOP (K) = FJT (K) FJBOT (K) = FJB (K) ! wavelength loop! enddo return end subroutine OPMIE !###################################################################### !<<<<<<<<<<<<<<<<<<<<<<