module MO_FASTJX_MOD !---------------------------------------------------------------------- ! FAST-JX Photolysis rate ! Revised by Junfeng Liu ! ! Feb 19, 2010 ! Junfeng.Liu !---------------------------------------------------------------------- ! Below is the original information from FASTJX code ! ! >>>>>>>>>>>>>>>>current code revised to JX ver 6.4 (8/08)<<<<<<<<<<<< ! ! 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 ! !<<<<<<<<<<<<<<<<<<<<<< O(1D) quantum yield character*7 :: TITLEJ(X_) !Species name in FJX_spec.dat character*7 :: TITLEJ2 !Temp var character*7 :: TITLEJ3 !Temp var real*8 :: TQQ(3,X_) !Temperature for supplied cross sections(K) real*8 :: QQQ(WX_,2,X_) !Supplied cross sections in each wavelength bin (cm2) !-----RD_MIE <- FJX_scat.dat------------------------------------------- integer :: NAA !Number of categories for scattering phase functions, is 30 integer :: NAA_AM3 !Number of categories for scattering phase functions for am3, is 880 real*8 :: MEE_AM3(5,A_AM3) !Aerosol mass extinction efficiency, MEE*colume mass =od real*8 :: SAA_AM3(5,A_AM3) !Single scattering albedo real*8 :: PAA_AM3(8,5,A_AM3) !Phase function: first 8 terms of expansion real*8 :: WAA_AM3(5,A_AM3) !Wavelengths for the NK supplied phase functions(nm) character*78 :: TITLE0 !Title of the first line in FJX_spec.dat character*78 :: TITLE0_AM3 !Title of the first line in am3_scat.dat integer :: JTAUMX !Max # is 250 real*8 :: ATAU ! is 1.120 real*8 :: ATAU0 ! is 0.010 character*20 :: TITLAA(A_) ! A_ =40 real*8 :: RAA(A_) !Effective radius associated with aerosol type real*8 :: DAA(A_) ! rho real*8 :: WAA(5,A_) !Wavelengths for the NK supplied phase functions(nm) real*8 :: QAA(5,A_) !Aerosol scattering phase functions real*8 :: SAA(5,A_) !Single scattering albedo real*8 :: PAA(8,5,A_) !Phase function: first 8 terms of expansion !-----RD_UM <- FJX_UMaer.dat------------------------------------------ character*20 :: TITLUM(33) !UM aerosol type names real*8 :: UMAER(3,5,21,33) !UM aerosol data !-----RD_JS <- chem_Js.dat-------------------------------------------- real*8 :: JFACTA(JVN_) !Quantum yield (or multiplication factor) for photolysis character*7 :: JLABEL(JVN_) !Reference label identifying appropriate J-value to use integer :: JIND(JVN_) !Set index arrays that map Jvalue(j) onto rates integer :: NRATJ !number of J species logical :: module_is_initialized = .false. real*8, parameter :: MASFAC = 100.d0*6.022d+23/(28.97d0*9.8d0*10.d0) !used to calc ZH real*8, parameter :: pfactor1 = 1.e3/9.8/28.9644*6.02e23/1.e4 !pa => molec/cm2 real*8, dimension(L1_) :: ZH ! CONTAINS ! ! ! Initialize data for fastjx calculation ! ! ! This subroutine initializes the calculation of photolysis rates ! for FAST-JX calculation ! ! !----------------------------------------------------------------------- subroutine fastjx_init !----------------------------------------------------------------------- ! 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 in !----------------------------------------------------------------------- ! ! IPH Channel number for reading all data files ! RAD Radius of Earth (cm) ! ZZHT Effective scale height above top of atmosphere (cm) ! DATUMX Maximum opt.depth above which sub layers should be inser ! SZAMAX Solar zenith angle cut-off, above which to skip calculat ! !----------------------------------------------------------------------- ! implicit none ! include 'parm_CTM.f' ! include 'parm_MIE.f' ! include 'cmn_metdat.f' ! include 'cmn_JVdat.f' integer I, J, K ! real*8, parameter:: PI = 3.141592653589793d0 ! real*8, parameter:: PI180=3.141592653589793d0/180.d0 if (module_is_initialized) return !---------------------------------------------------------------------- ! output version number and tagname to logfile. !---------------------------------------------------------------------- call write_version_number(version, tagname) ! ! Defaults & constants RAD = 6375.d5 !cm earth radii ZZHT = 5.d5 ! STT(:,:,:,:) = 1.d6 ! TCNAME(:) = 'CO2' ! ! Read in annual mean pressure field !RSH if this code is ever reactivated, change to use mpp_open ! open (IPH,file='ECT42_grid.dat',form='formatted',status='old') ! read (IPH,*) ! read (IPH,'(16F8.2)') (XDGRD(I),I=1,I_) ! read (IPH,*) ! read (IPH,'(16F8.2)') (YDGRD(J),J=1,J_) ! read (IPH,*) ! read (IPH,'(16F8.2)') (AREAXY(1,J),J=1,J_) ! read (IPH,*) ! read (IPH,'(16F8.1)') ((PMEAN(I,J),I=1,I_),J=1,J_) ! close(IPH) !RSH call mpp_close (IPH) ! do J=1,J_ ! AREAXY(1,J) = AREAXY(1,J)*1.d9 ! do I=2,I_ ! AREAXY(I,J) = AREAXY(1,J) ! enddo ! enddo ! do I = 1,I_ ! XGRD(I) = XDGRD(I)*PI180 ! enddo ! do J = 1,J_ ! YGRD(J) = YDGRD(J)*PI180 ! enddo ! ! Read in fast-J X-sections (spectral data) <<<<<<<<<<<<<< new fast-JX call RD_XXX('INPUT/FJX_spec.dat') ! Read in aerosol/cloud scattering data <<<<<<<<<<<<<<<<<< new fast-JX call RD_MIE('INPUT/FJX_scat.dat') ! Read in aerosol/cloud scattering data <<<<<<<<<<<<<<<<<< new fast-JX call RD_MIE_AM3('INPUT/am3_scat.dat') ! Read in UMich aerosol optical data <<<<<<<<<<<<<<<<<< new fast-JX 6 ! call RD_UM (IPH,'INPUT/FJX_UMaer.dat') ! Read in labels of photolysis rates required >>>>> 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('INPUT/chem_Js.dat') ! ! Read in T & O3 climatology >>>> general backup clim ! call RD_PROF(IPH,'atmos_std.dat') !----------------------------------------------------------------------- ! mark module as initialized. !----------------------------------------------------------------------- module_is_initialized = .true. end subroutine fastjx_init !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! include 'cmn_metdat.f' for fast-JX code v5.3+ (prather 6/05) ! ! needs 'parm_ctm.f' for dimensions ! delivers p, T, Surf Albedo, and Optical Depth from CTM to fastJX ! >>>>this is for standalone fast-JX ver 5.3 (6.05) !----------------------------------------------------------------------- ! ! real*8 P(I_,J_) ! Surface pressure ! real*8 T(I_,J_,L_) ! Temperature profile ! ! real*8 XGRD(I_) ! Longitude (midpoint, radians) ! real*8 XDGRD(I_) ! real*8 YGRD(J_) ! Latitude (midpoint, radians) !! real*8 YDGRD(J_) ! real*8 ETAA(L1_) ! Eta(a) value for level boundaries ! real*8 ETAB(L1_) ! Eta(b) value for level boundaries ! real*8 AREAXY(I_,J_) ! area (m^2) ! integer MONTH ! integer NSLAT ! Latitude(J) index of current column ! integer NSLON ! Longitude(I) index of current column ! Note that the climatology for O3, T, clouds, etc include a layer above ! the CTM (L1_=L_+1) for calculating the J's ! real*8, dimension(L1_) :: RELH ! Rel Hum (0.00 to 1.00) ! real*8, dimension(I_,J_,L1_) :: TJ, DM, DO3, ZH ! real*8, dimension(I_,J_,L1_) :: CLDLWP,AER1P, AER2P ! integer,dimension(I_,J_,L1_) :: CLDNDX,AER1N, AER2N ! real*8, dimension(I_,J_) :: PMEAN, SA ! 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 ! ! real*8 STT(I_,J_,L_,NTR_) ! real*8 TREF(51,18,12),OREF(51,18,12) ! character*10 TCNAME(NTR_) ! ! common/metdat/P,T, STT, & ! & XGRD,YGRD,XDGRD,YDGRD,ETAA,ETAB,AREAXY, & ! & TREF,OREF,PMEAN,SA, MONTH,NSLAT,NSLON, TCNAME ! ! common /jvdatIJ/TJ,DM,DO3,ZH,RELH, CLDLWP,AER1P,AER2P, & ! & CLDNDX,AER1N,AER2N !----------------------------------------------------------------------- !-------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 !----------------------------------------------------------------------- ! 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(18,3) O2 cross-sections ! QO3(18,3) O3 cross-sections ! Q1D(18,3) O3 => O(1D) quantum yield ! TQQ(3,64) Temperature for supplied cross sections ! QQQ Supplied cross sections in each wavelength bin (cm2) ! !----------------------------------------------------------------------- ! common /jvchem/JFACTA,JIND,NRATJ,JLABEL ! ! common /jvdat/WBIN,WL,FL,QO2,QO3,Q1D,QQQ,QRAYL,TQQ, & ! & ZZHT,ATAU,ATAU0, WAA,QAA,PAA,SAA,RAA,DAA,RAD, UMAER, & ! & JTAUMX, NJVAL,NW1,NW2,NAA ,TITLE0,TITLEJ,TITLAA,TITLUM !----------------------------------------------------------------------- ! fastjx_photo should be called by fphoto, its function is to provide ! the calculated J-values to fphoto, based on clould and ! aerosol information ! subroutine fastjx_photo( U0, SOLF, phalf1, zhalf1, pfull1, tfull, & XO3, pwt, lwc, cloudsf, clouds_ndx ,qfld, & srf_alb, aerop, aeron, o3_column_top, ZPJ) !----------------------------------------------------------------------- ! ! PHOTOJ is the gateway to fast-JX calculations: ! only access to CTM 3-D GLOBAL arrays ! sets up the 1-D column arrays for calculating J's ! v6.1 calculates the optical properties for each wavelength ! and passes single column, single wavelength data to OPMIE ! OPMIE no longer needs aerosol data, just the overall properties ! DTAUX() = opt depth of each layer, ! POMEGAX(8,) = scat phase fn, includes s-s albedo factor. !----------------------------------------------------------------------- implicit none ! include 'parm_CTM.f' ! include 'parm_MIE.f' ! include 'cmn_metdat.f' ! include 'cmn_JVdat.f' real*8, intent(in) :: U0, & !cos(SZA) SOLF !solar-earth distance factor !In AM3, the first layer is top of ATM, !but in FASTJ, the first layer is surface real*8, intent(in) :: phalf1(L1_),& !Pressure at boundaries (Pa) 49 layer zhalf1(L1_),& !height at layer boundaries pfull1(L_),& !Pressure at mid-layer (Pa) 48 tfull(L_) !Temperature at mid-Layer (K) real*8, intent(in) :: XO3(L_) !O3 mixing ratio (VMR) real*8, intent(in) :: pwt(L_) !Column air density (Kg/m2) real*8, intent(in) :: lwc(:,:) !(L_,ncloud=2)Liquid water content (Kg/Kg) real*8, intent(in) :: cloudsf(:,:) !Clouds fraction integer, intent(in) :: clouds_ndx(:,:) ! clouds type index =2: 1 water, 2 ice real*8, intent(in) :: qfld(L_) !Specific humidity (kg/kg) real*8, intent(in) :: srf_alb !Surface albedo real*8, intent(in) :: aerop(:,:) !Aerosol mass path (g/m2) (L_,AERONUM) integer, intent(in) :: aeron(:,:) !(AERONUM) Aerosol type index, see FJX_scat.dat !for aerosol type info real*8, intent(in) :: o3_column_top !O3 column density, a small number real*8, intent(out) :: ZPJ(JVL_,JVN_) !ZPJ(JVL_,JVN_) 2-D array of J-values !----------------------------------------------------------------------- real*8 :: SZA !solar zenith angle real*8 :: FREFL !fraction of energy refle real*8 :: ZH(L1_) !height of each layer (cm) real*8 :: SCALEH !scale height (cm) real*8 :: pz, qs, wrk, es !parameters to calc RH real*8 :: phalf(L1_), pfull(L_) !pressures at boundaries and mid-layer (hPa) real*8, dimension(L1_+1) :: TTJ,& !Temperature used by FAST-J(K) DDJ,& !Air column density (molecs/cm2) for each layer ZZJ,& !O3 column density (molecs/cm2) ZHL !Distance of each layer to surface (cm), same to ZH, !but have another top layer, used to calculate cloud density real*8, dimension(L1_+1) :: PPJ !pressure at boundaries (hPa) !--------key amtospheric data needed to solve plane-parallel J--------- integer,dimension(L2_+1) :: JXTRA real*8, dimension(W_) :: FJTOP,& FJBOT,& FSBOT,& FLXD0,& RFL real*8, dimension(L_, W_) :: AVGF,& FJFLX real*8, dimension(L1_,W_) :: DTAUX,& FLXD real*8, dimension(8,L1_,W_) :: POMEGAX real*8, dimension(W_,L1_) :: FFX real*8, dimension(W_,8) :: FFXNET !---flux/heating arrays (along with FJFLX,FLXD,FLXD0) real*8 :: FLXJ(L1_),& FFX0,& FXBOT,& FABOT real*8 :: ODABS,& ODRAY real*8 :: RFLECT,& FREFS,& FREFI !SOLF real*8 :: AMF2(2*L1_+1,2*L1_+1) !------------key SCATTERING arrays for clouds+aerosols------------------ real*8 :: OPTX(5),& SSAX(5),& SLEGX(8,5) real*8 :: OD(5,L1_),& SSA(5,L1_),& SLEG(8,5,L1_) real*8 :: OD600(L1_) real*8 :: PATH,PATH1,PATH2,& DENSWP,& DENS,& RH,& RHO,& REFF,& XTINCT,& ODCLD,ODCLD1,ODCLD2 !optical depth of cloud !------------key arrays AFTER solving for J's--------------------------- real*8 :: FFF(W_,JVL_),& VALJ(X_) real*8 :: FLXUP(W_),& FLXDN(W_),& DIRUP(W_),& DIRDN(W_) !------------2-D array of J_s returned by JRATET------------------------ real*8 :: VALJL(JVL_,NJVAL) integer :: I,J,K,& L,M,& KMIE,KW,& NAER,NDCLD,NDCLD1,NDCLD2, & RATIO(W_) real*8 :: XQO3,XQO2,& DTAUC,& WAVE, & ! FLINT, & TTT logical :: test_J = .false. !-------------check model initialization ------------------------------- if (.not. module_is_initialized) & call error_mesg ('ATMOS: fastjx_init:','fastjx_init must be called first.', FATAL) !----------------------------------------------------------------------- ZPJ(:,:) = 0.d0 FFF(:,:) = 0.d0 FREFI = 0.d0 FREFL = 0.d0 FREFS = 0.d0 !----------------------------------------------------------------------- ! call SOLARZ(UTIME,IDAY,YGRD,XGRD1, SZA,U0,SOLF) !----------------------------------------------------------------------- !-----SOLF = 1.d0 ! this needs to be dropped to include 6.7% annual cy SZA = acos(U0)*180./ 3.141592653589793d0 !-----check for dark conditions SZA > 98.0 deg => tan ht = 63 km ! or 99. 80 km if (SZA .gt. SZAMAX) goto 99 !-----load the amtospheric column data phalf=phalf1/100. pfull=pfull1/100. !-----calculate effective altitude of each CTM level edge (cm) ! ZH(1) = 16d5*log10(1000.d0/phalf(L1_)) ! do L = 1,L_ ! SCALEH = 1.3806d-19*MASFAC*tfull(L1_-L) ! ZH(L+1) = ZH(L) -(log(phalf(L1_-L)/phalf(L1_-L+1))*SCALEH) ! enddo do L=1,L1_ ZH(L)= zhalf1(L1_+1-L)*100. ! height at boundaries (cm) end do ! if (mpp_pe() == mpp_root_pe() ) then ! write(*,*) 'fastjx: zh(cm)=', zh ! end if !---calculate spherical weighting functions (AMF: Air Mass Factor) do L = 1,L1_ ZHL(L) = ZH(L) enddo ZHL(L1_+1) = ZHL(L1_) + ZZHT do L = 1,L1_ !L=1 surface AM3 L=1 TOA ! PPJ(L) = ETAA(L) + ETAB(L)*P(ILNG,JLAT) ! TTJ(L) = TJ(ILNG,JLAT,L) ! DDJ(L) = DM(ILNG,JLAT,L) ! ZZJ(L) = DO3(ILNG,JLAT,L) PPJ(L) = phalf(L1_+1-L) ! pressure at boundaries (hPa) if(L molec/cm2 ZZJ(L) = o3_column_top !molec/cm2 end if enddo PPJ(L1_+1) = 0.d0 !----------------------------------------------------------------------- call SPHERE2 (U0,RAD,ZHL,ZZHT,AMF2,L1_) !calculate how much sunlight is blocked by earth ! !at a few lowest layers when zenith angle is high !----------------------------------------------------------------------- !---calculate the optical properties (opt-depth, single-scat-alb, phase-fn(1:8) !--- 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 OD(K,L) = 0.d0 SSA(K,L) = 0.d0 do I=1,8 SLEG(I,K,L) = 0.d0 enddo enddo enddo !-----sst and phase-fn are weighted by optical depth do L = 1,L1_ ! do M = 1,size(lwc,2) !-----cloud in layer: if valid cloud index (now 4:13) if(L90, 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 fastjx_photo !<<<<<<<<<<<<<<<<<<<<<<>>>>>>>>>>>>>>>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 do J = 1,NJVAL T_JX = TITLEJ(J) do K = 1,NRATJ T_CHEM = JLABEL(K) if (T_CHEM(1:6) .eq. T_JX(1:6)) JIND(K)=J enddo enddo ! 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 return END subroutine RD_JS !----------------------------------------------------------------------- ! subroutine RD_PROF(NJ2,NAMFIL) !----------------------------------------------------------------------- ! Routine to input T and O3 reference profiles !----------------------------------------------------------------------- ! implicit none ! include 'parm_CTM.f' ! include 'cmn_metdat.f' ! integer, intent(in) :: NJ2 ! character(*), intent(in) :: NAMFIL ! ! integer IA, I, M, L, LAT, MON, NTLATS, NTMONS, N216 ! real*8 OFAC, OFAK ! character*72 TITLE0 ! !RSH if this code is ever reactivated, change to use mpp_open: ! open (NJ2,file=NAMFIL,status='old',form='formatted') !! read (NJ2,'(A)') TITLE0 ! read (NJ2,'(2I5)') NTLATS,NTMONS ! write(6,'(1X,A)') TITLE0 ! write(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) !RSH call mpp_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 ! return ! 1000 format(1x,'std atmos profiles: ',i3,' lat x ',i2,' mon') ! END !----------------------------------------------------------------------- ! subroutine SET_CLD0(TINIT,CLDP,NCLD,ALBEDO) !----------------------------------------------------------------------- ! Routine to set cloud and surface properties: now loads the input the ! input profiles (FJX_test.dat) INCLUDING Temperatures and Surf Albe ! ! >>>>>> this subroutine will need to be customized ! >>>>>> it is separate from aerosols since it comes from met fields !----------------------------------------------------------------------- ! implicit none !! include 'parm_CTM.f' ! include 'parm_MIE.f' ! include 'cmn_metdat.f' ! include 'cmn_JVdat.f' ! integer, intent(in) :: NCLD(L_) ! real*8, intent(in) :: TINIT(L_),CLDP(L_),ALBEDO ! integer I, J, L ! ! do J = 1,J_ ! do I = 1,I_ ! do L = 1,L_ ! CLDLWP(I,J,L) = CLDP(L) ! CLDNDX(I,J,L) = NCLD(L) ! T(I,J,L) = TINIT(L) ! enddo ! SA(I,J) = ALBEDO ! enddo ! enddo ! ! return ! END !----------------------------------------------------------------------- ! subroutine SET_AER0 (AER1,AER2, NAA1,NAA2) !----------------------------------------------------------------------- ! Set up aerosols >>>>customize for climatology or CTM as source. ! this needs to load in the aerosol mass/column and type in each lay ! In a CTM, this should be done as a climatology or in real time using ! computed STT(kg) of the aerosol and calculating PATH=1000.*STT/ARE !----------------------------------------------------------------------- ! implicit none ! include 'parm_CTM.f' ! include 'parm_MIE.f' ! include 'cmn_metdat.f' ! include 'cmn_JVdat.f' ! integer, intent(in) :: NAA1(L_),NAA2(L_) ! real*8, intent(in) :: AER1(L_),AER2(L_) ! integer I, J, L, K ! AER1P(:,:,:) = 0.d0 !! AER2P(:,:,:) = 0.d0 ! AER1N(:,:,:) = 0 ! AER2N(:,:,:) = 0 ! do J = 1,J_ ! do I = 1,I_ ! do L = 1,L_ ! AER1P(I,J,L) = AER1(L) ! AER2P(I,J,L) = AER2(L) ! AER1N(I,J,L) = NAA1(L) ! AER2N(I,J,L) = NAA2(L) ! enddo ! enddo ! enddo ! END !----------------------------------------------------------------------- ! 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.f' ! include 'cmn_metdat.f' ! integer I, J, K, L, M, N ! real*8 PSTD(52),OREF2(51),TREF2(51),PJ(L1_+1) ! 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 ! PJ(L1_+1) = 0.d0 ! do K = 1,L1_ ! PJ(K) = ETAA(K) + ETAB(K)*PMEAN(I,J) ! enddo ! 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 !<<<<<<<<<<<<<<<<<<<<<<<1.0 ! real*8, intent(in):: RELH ! index of cloud/aerosols ! integer,intent(inout):: L ! integer I,J ! real*8 XTINCT, REFF,RHO ! if (L .gt. NAA .or. L .lt. 4) then ! write(6,*) 'ATMOS:fastjx_photo: aerosol index out-of-range: L/NAA',L,NAA ! L = 18 ! endif ! if (L .lt. 14) then ! write(6,*) 'ATMOS:fastjx_photo: aerosol as cloud: L/NAA',L,NAA ! endif ! REFF = RAA(L) ! RHO = DAA(L) ! do J=1,5 !---extinction K(m2/g) = Q(wvl) / [4/3 * Reff(micron) * aerosol-density( ! XTINCT = 0.75d0*QAA(J,L)/(REFF*RHO) ! OPTD(J) = PATH*XTINCT ! 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.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' ! 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,*) 'ATMOS:fastjx_init: 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 (hPa) ! 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.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' ! real*8, intent(in) :: PPJ(L1_+1),TTJ(L1_) real*8, intent(in) :: PPJ(L_),TTJ(L_) !set to pfull (hPa) and tfull from AM3 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(L1_-L) ! if (L .eq. 1) then PP = PPJ(L1_-L) !PPJ(1) ! else ! PP = (PPJ(L)+PPJ(L+1))*0.5d0 ! endif DD = 7.24e18*PP/TT !molec/cm3 do J = 1,NJVAL VALJ(J) = 0.d0 enddo !j 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) -QO31D) *FFF(K,L) !changed by jul VALJ(3) = VALJ(3) + QO31D*FFF(K,L) enddo do J = 4,NJVAL if (TQQ(2,J) .gt. TQQ(1,J)) then TFACT = max(0.0,min(1.0,(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.0, min(1.0, 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.0, min(1.0, TFAC0)) !---IV=NJVAL+1 = Q1A for Acetone => (1), allow full range of T = 200K-30 IV = NJVAL+1 TT200 = min(300.0, max(200.0, 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 do J=1,NJVAL VALJL(L,J) = VALJ(J) ! if(VALJ(J) .lt. 0.D0)then ! write(*,*)'fastjx: !!!! negative J: J=',J, ' TITLE=',TITLEJ(J),' VALJ(J)=',VALJ(J),' L=',L ! endif enddo enddo ! master loop over L=1,JVL_ return END subroutine JRATET !----------------------------------------------------------------------- ! subroutine JP_ATM(PPJ,TTJ,DDJ,ZZJ,ZHL,ZZHT,DTAUX,POMEGAX,JXTRA) !----------------------------------------------------------------------- ! implicit none ! include 'parm_CTM.f' !----------------------------------------------------------------------- !--------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_XXX(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(18,3) O2 cross-sections ! QO3(18,3) O3 cross-sections ! Q1D(18,3) O3 => O(1D) quantum yield ! TQQ(3,64) Temperature for supplied cross sections ! QQQ Supplied cross sections in each wavelength bin (cm2) ! !----------------------------------------------------------------------- !jul++ ! RD_XXX ! will access WL(18), FL(18), QRAYL(18) ! QO2(18,3) O2 cross-sections ! QO3(18,3) O3 cross-sections ! Q1D(18,3) O3 => O(1D) quantum yield ! TITLEJ(64), SPECIES NAME O2 O3 ! TQQ(3, 64), TEMPERATURE 180k, 260k, 300K ! QQQ(18,2,64) !jul-- implicit none ! include 'parm_CTM.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' character(*), intent(in) :: NAMFIL integer NJ1 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 call mpp_open (NJ1, trim(NAMFIL), MPP_RDONLY, MPP_ASCII, & MPP_SEQUENTIAL, MPP_MULTI, MPP_SINGLE) 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) 'ATMOS:fastjx_init: RD_XXX stop',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 if (mpp_pe() == mpp_root_pe()) then write(*,*)'ATMOS:fastjx_init: RD_XXX:' write(*,200) J,TITLEJ(J),(TQQ(I,J),I=1,3) end if 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!) if (mpp_pe() == mpp_root_pe()) then write(*,200) JJ,TITLEJ(JJ),(TQQ(I,JJ),I=1,2) end if 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 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 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 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 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,*) 'ATMOS:fastjx_init: 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' call mpp_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(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.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' character(*), intent(in) :: NAMFIL integer I, J, K, NJ1 !jul++ ! NAA,TITLE0 ! JTAUMX,ATAU,ATAU0 ! TITLAA(J),RAA(J),DAA(J) ! WAA(K,J),QAA(K,J),SAA(K,J),PAA(I,K,J) ! !jul-- call mpp_open (NJ1, trim(NAMFIL), MPP_RDONLY, MPP_ASCII, & MPP_SEQUENTIAL, MPP_MULTI, MPP_SINGLE) read (NJ1,'(i2,a78)') NAA,TITLE0 if (NAA .gt. A_) then write(*,*) 'ATMOS:fastjx_init: too many scat-data sets:', NAA, A_ stop endif read (NJ1,'(5x,i5,2f10.5)') JTAUMX,ATAU,ATAU0 if (mpp_pe() == mpp_root_pe()) then write(*,*)'ATMOS:fastjx_init: RD_MIE:' write(*,'(a,2f9.5,i5)') ' ATAU/ATAU0/JMX',ATAU,ATAU0,JTAUMX end if 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 call mpp_close (NJ1) if (mpp_pe() == mpp_root_pe()) then write(*,'(a,9f8.1)') ' ATMOS:fastjx_init: RD_MIE: Aerosol optical: r-eff/rho/Q(@wavel):' & ,(WAA(K,1),K=1,5) write(*,*) TITLE0 end if do J=1,NAA if (mpp_pe() == mpp_root_pe()) then write(*,*) ' ATMOS:fastjx_init: RD_MIE:' write(*,'(i3,1x,a8,7f8.3)') & J,TITLAA(J),RAA(J),DAA(J),(QAA(K,J),K=1,5) end if enddo return END subroutine RD_MIE !----------------------------------------------------------------------- subroutine RD_MIE_AM3(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.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' character(*), intent(in) :: NAMFIL integer I, J, K, NJ1 character*90 :: str1 !jul++ ! NAA,TITLE0 ! JTAUMX,ATAU,ATAU0 ! TITLAA(J),RAA(J),DAA(J) ! WAA(K,J),QAA(K,J),SAA(K,J),PAA(I,K,J) ! !jul-- call mpp_open (NJ1, trim(NAMFIL), MPP_RDONLY, MPP_ASCII, & MPP_SEQUENTIAL, MPP_MULTI, MPP_SINGLE) read (NJ1,'(i4,a78)') NAA_AM3, TITLE0_AM3 if (NAA_AM3 .gt. A_AM3) then write(*,*) 'ATMOS:fastjx_init: too many scat-data sets for AM3:', NAA_AM3, A_AM3 stop endif ! read (NJ1,'(5x,i5,2f10.5)') JTAUMX,ATAU,ATAU0 ! if (mpp_pe() == mpp_root_pe()) then ! write(*,*)'ATMOS:fastjx_init: RD_MIE_AM3:' ! write(*,'(a,2f9.5,i5)') ' ATAU/ATAU0/JMX',ATAU,ATAU0,JTAUMX ! end if ! read (NJ1,*) do J = 1,NAA_AM3 read (NJ1,'(a90)') str1 if (mpp_pe() == mpp_root_pe()) then write (*,'(a12,a90)') & 'fastjx_init:',str1 endif ! 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,2e14.6,7f6.3,1x,f7.3,f8.4)') & WAA_AM3(K,J),MEE_AM3(K,J),SAA_AM3(K,J),(PAA_AM3(I,K,J),I=2,8) PAA_AM3(1,K,J) = 1.d0 if (mpp_pe() == mpp_root_pe()) then !!! write (*,'(a12,i4,1x,f4.0,2e14.6,8f6.3,1x,f7.3,f8.4)') & !!! 'fastjx_init:',J,WAA_AM3(K,J),MEE_AM3(K,J),SAA_AM3(K,J),(PAA_AM3(I,K,J),I=1,8) endif enddo enddo call mpp_close (NJ1) ! if (mpp_pe() == mpp_root_pe()) then ! write(*,'(a,9f8.1)') ' ATMOS:fastjx_init: RD_MIE: Aerosol optical: r-eff/rho/Q(@wavel):' & ! ,(WAA_AM3(K,1),K=1,5) !! write(*,*) TITLE0_AM3 ! end if ! if (mpp_pe() == mpp_root_pe()) then ! do J=1,NAA_AM3 ! write(*,*) ' ATMOS:fastjx_init: RD_MIE_AM3:' ! write(*,'(i3,1x,a8,7f8.3)') & ! J,TITLAA_AM3(J),(QAA(K,J),K=1,5) ! end if ! enddo return END subroutine RD_MIE_AM3 !----------------------------------------------------------------------- subroutine RD_UM(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.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' !jul++ ! TITLUM(L) ! UMAER !jul-- character(*), intent(in) :: NAMFIL integer I, J, K, L, NJ1 call mpp_open (NJ1, trim(NAMFIL), MPP_RDONLY, MPP_ASCII, & MPP_SEQUENTIAL, MPP_MULTI, MPP_SINGLE) 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 call mpp_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(out) :: SZA,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 ! if (mpp_pe() == mpp_root_pe()) then ! write(*,*) 'FASTJ:SOLARZ: XGRDI,YGRDJ',XGRDI,YGRDJ ! write(*,*) 'FASTJ:SOLARZ: LOCT (rad)',LOCT ! write(*,*) 'FASTJ:SOLARZ: SINDEC,COSDEC', SINDEC,COSDEC ! write(*,*) 'FASTJ:SOLARZ: SINLAT,COSLAT', SINLAT,COSLAT ! write(*,*) 'FASTJ:SOLARZ: COS, SZA',COSSZA,SZA ! end if ! 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.0) ! 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 continue 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 ATAU ! ! 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 chosen !--- 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=100 !--- ATAU = 1.12 now recommended for more -accurate heating rates (not J's) !----------------------------------------------------------------------- ! implicit none integer, intent(in) :: JTAUMX,L1X,L2X !index of cloud/aerosol integer, intent(in) :: NX !Mie scattering array size real*8, intent(in) :: DTAUX(L1X) !cloud+3aerosol OD in each layer real*8, intent(in) :: ATAU,ATAU0 integer, intent(out):: JXTRA(L2X+1) !number of sub-layers to be added ! 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 code !---In the original fast-J, equal sub-layers were chosen, this is wasteful !---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 into !--- 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 ATAUMN !---These are hardwired below, can be changed, but have been tested/optimized 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 ATAU>1 !---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 levels JTOTL = L2X + 2 do L2 = L2X,1,-1 JTOTL = JTOTL + JXTRA(L2) if (JTOTL .gt. NX/2) then write(6,'(A,2I5,F9.2)') 'ATMOS:fastjx_photo: N_/L2_/L2-cutoff JXTRA:',NX,L2X,L2 do L = L2,1,-1 JXTRA(L) = 0 enddo go to 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.f' ! include 'parm_MIE.f' ! include 'cmn_JVdat.f' real*8, intent(in) :: DTAUX(L1_,W_),& ! optical depth POMEGAX(8,L1_,W_) ! scattering phase fn 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 abledo) ! 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 layer-L ! AMF2 now does both edges and middle of CTM layers ! JXTRA(1:L1_) = number 0:J = no. of additional levels to be inserted ! out: ! FJACT(1:L_,1:W_) = mean actinic flux(diff+direct) at std CTM levels(mid-lyr) ! (new ver 5.7 diagnostics for fluxes, deposition) fluxes 'down' are <0 ! FJTOP(1:W_) = diffuse flux out top-of-atmosphere (TAU=0 above top model lyr) ! 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 dim!!) ! FLXD(1:L_+1,1:W_) = solar flux deposited in layer L (includes lyr above CTM) ! 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 of atm) ! 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 pomega ! 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 code. ! there are twice the number of layers in the LZ arrays (2*L2_ + 2*JADDTO + 1) ! 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 OD_s ! 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 (L2) ! JNDLEV(L) = L2 index that maps on CTM mid-layer L ! !---JADDLV(L2=1:L2_) = number of levels to add between TTAU2(L2) and TTAU(L2+1) !--- JADDLV is taken from JXTRA, which is based on visible OD. !--- JADDTO(L2=1:L2_+1) is the cumulative number of levels to be added !---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)') 'ATMOS:fastjx_photo: 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 * AirMassFactor) !--- 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 above 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) !!!!!!!!good place to 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 layer !--- note FLXD0 .ne. (1.d0 - FTAU(L_+1))/AMF(L_+1,L_+1) with spherical atmos 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 ! Ftau2 indicate how much of light is still alive at that layer do L2 = 1,L2_+1 ! L2 = index of CTM edge- and mid-layers L2L = L2LEV(L2) ! L2L = index for L2 in expanded scale(JADD) LZ = ND + 2 - 2*L2L ! LZ = index for L2 in scatt arrays 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 !Loop L endif enddo ! Loop L2 ! Now fill in the even points with simple interpolation in scatter arrays 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 enddo ! wavelength loop! !----------------------------------------------------------------------- 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: 4* + solar at mid-layer 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 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 !---NB if one needs the mean intensity throughout layer L (instead of mi !--- then average (tau-weighted) the odd-points from: NELEV(L-1) to NE !---NB This is NOT now used. ! real*8 FJACT2(L_,W_) !-----LZ are indices to the 1:ND array LZ1 > LZMID > LZ0 ! LZ1 = ND ! do L = 1,L_ ! LZMID = ND+2-2*JNDLEV(L) ! LZ0 = ND+2-2*JNELEV(L) ! SUMT = 0.d0 ! SUMJ = 0.d0 ! do L2 = LZ0,LZ1-2,2 ! SUMT = SUMT + ZTAU(L2+2,K)-ZTAU(L2,K) ! SUMJ = SUMJ + (ZTAU(L2+2,K)-ZTAU(L2,K))*(FJ(L2,K)+FJ(L2+2,K)) ! enddo ! FJACT2(L,K) = 2.d0*SUMJ/SUMT + FZ(LZMID,K) ! LZ1 = LZ0 ! enddo !---diffuse fluxes reflected at top, incident at bottom FJTOP(K) = FJT(K) FJBOT(K) = FJB(K) ! wavelength loop! enddo !if( any(FJACT(:,:) < 0.0D0 )) then ! write(*,*) 'Mie: find negative FJACT' ! write(*,*) 'Mie: U0=', U0, ' ND=',ND ! if (mpp_pe() == mpp_root_pe() ) then ! write(*,*) 'Mie: FJ diffusive', FJ ! write(*,*) 'Mie: FZ solar at mid-layer', FZ ! write(*,*) 'Mie: ZTAU', ZTAU ! write(*,*) 'Mie: RFL', RFL ! write(*,*) 'Mie: FJT', FJT ! write(*,*) 'Mie: FJB', FJB ! write(*,*) 'Mie: POMEGA', POMEGA ! write(*,*) 'Mie: ZFLUX', ZFLUX ! write(*,*) 'Mie: AMF2', AMF2 ! endif !end if return END subroutine OPMIE !<<<<<<<<<<<<<<<<<<<<<<