!### macro's ##################################################### ! #define TRACEBACK write (gol,'("in ",a," (",a,", line",i5,")")') rname, __FILE__, __LINE__; call goErr #define IF_NOTOK_RETURN(action) if (status/=0) then; TRACEBACK; action; return; end if #define IF_ERROR_RETURN(action) if (status> 0) then; TRACEBACK; action; return; end if ! #include "tm5.inc" ! !################################################################# module photolysis use GO, only : gol, goPr, goErr ! use binas, only : pi ! use photolysis_data, only : phot_dat ! use dims implicit none ! --- in/out ------------------------------------- private !public :: get_theta, nbands, fparam, photolysis_rates, phot_dat public :: init_photolysis, end_photolysis !public :: cloud_info, ozone_info ! --- const ----------------------- character(len=*), parameter :: mname = 'photolysis' ! ! vars ! ! integer,parameter :: nbands=7 ! real ,parameter :: todu = 3.767e-17 !from #/cm2 --> du ! real ,parameter :: to_m = 3.767e-22 !from #/cm2 --> m (for o2) ! real ,parameter :: sp = 6.022e23*1.e-4*0.2095/(28.964*1e-3*9.81) ! ! sp from pa ---> #o2/cm2 ! ! real,dimension(nbands),parameter :: crs_o3 = & ! (/ 3.6448e-19, 1.7115e-18, 2.8218e-19, 1.1091e-19, & ! 2.8145e-20, 1.9000e-23, 4.5500e-21/) ! ! crs_o2: cross-section o2 205 nm ! real,parameter :: crs_o2= 7.63e-24 ! ! ! reference tables (pressure,mu0,band) ! real,dimension(0:102,0:19,nbands) :: fint ! ! extraterrestial fluxes at selected wavelengths ! real,dimension(nbands) :: fe contains subroutine init_photolysis !(iunit,o3du) use dims , only : nregions, im, jm use dims , only : newsrun use photolysis_data, only : phot_dat !#ifdef MPI ! use mpi_const, only: root,localComm, ierr, my_real, myid, lmloc !#endif ! use chem_param, only : d23_data ! --- in/out -------------------------------------- ! ! input ! integer,intent(in) :: iunit ! type(d23_data),dimension(nregions),intent(in) :: o3du ! --- local ---------------------------------------- integer :: region, imr, jmr !integer :: i,m,lmr !real :: lat ! --- begin ------------------------------------------ ! ! routine to read some reference tables for photolysis calculations ! ! to be called at the start of any photolysis calculation ! ! if(newsrun) then !#ifdef MPI ! if(myid == root) then !#endif ! open(unit=iunit,file='fintwb.data',status='old') ! ! read(iunit,*) ! read(iunit,*) ! read(iunit,*) (fe(i),i=1,7) ! read(iunit,*) fint ! ! close(iunit) !#ifdef MPI ! end if ! call mpi_bcast(fint, 103*20*nbands , my_real, root, localComm, ierr) ! call mpi_bcast(fe, 7 , my_real, root, localComm, ierr) !#endif ! end if ! ! allocate memory photolysis and initialise: do region=1,nregions imr = im(region) ; jmr = jm(region) !#ifdef MPI ! lmr = lmloc !#else ! lmr = lm(region) !#endif ! if(newsrun) allocate(phot_dat(region)%o3klim_top(jmr)) ! phot_dat(region)%o3klim_top(:) = o3du(region)%d23(:,lm(region)) if ( newsrun ) then ! allocate(phot_dat(region)%vo3_surf (imr,jmr)) ! allocate(phot_dat(region)%taus (imr,jmr)) ! allocate(phot_dat(region)%cloudfrac (imr,jmr)) allocate(phot_dat(region)%albedo (imr,jmr)) ! allocate(phot_dat(region)%vo3s_av (imr,jmr)) ! allocate(phot_dat(region)%taus_av (imr,jmr)) !JEW: allocate(phot_dat(region)%ags_av (imr,jmr)) ! allocate(phot_dat(region)%cfr_av (imr,jmr)) ! phot_dat(region)%albedo = 0.05 ! phot_dat(region)%cfr_av = 0.0 ! phot_dat(region)%taus_av = 0.0 ! phot_dat(region)%ags_av = 0.0 ! phot_dat(region)%ncloud_av = 0 ! phot_dat(region)%vo3s_av = 0.0 ! phot_dat(region)%nvo3s_av= 0 !JEW: phot_dat(region)%nalb_av= 0 ! ! allocate(phot_dat(region)%vo3 (imr,jmr,lmr)) ! allocate(phot_dat(region)%tauc(imr,jmr,lmr)) ! allocate(phot_dat(region)%jo3_av(imr,jmr,lm(region))) ! allocate(phot_dat(region)%jno2_av(imr,jmr,lm(region))) ! phot_dat(region)%jo3_av = 0.0 ! phot_dat(region)%jno2_av = 0.0 ! phot_dat(region)%nj_av = 0 end if end do end subroutine init_photolysis ! real function get_theta(daynr,lon_deg,lat_deg,time) ! implicit none ! ! ! input ! integer,intent(in) :: daynr ! real,intent(in) :: lat_deg,lon_deg ! real,intent(in) :: time ! ! ! local ! real :: lon_rad,lat_rad,obliq, deday, delta,houra ! ! lat_rad=lat_deg*pi/180. ! lon_rad=lon_deg*pi/180. ! obliq = 23.45 * pi/180. ! deday = 4.88 + 2.*pi/365. * daynr ! delta = asin(sin(obliq)*sin(deday)) ! houra = lon_rad -pi + time * (2.*pi/24.) ! ! get_theta= acos(sin(delta)*sin(lat_rad) + & ! cos(delta)*cos(lat_rad)*cos(houra))*180./pi ! ! end function get_theta ! ! ! ! real function sundis(imonth,iday) ! !-----------------------------------------------------------------------------* ! != purpose: =* ! != calculate earth-sun distance variation for a given date. based on =* ! != fourier coefficients originally from: spencer, j.w., 1971, fourier =* ! != series representation of the position of the sun, search, 2:172 =* ! !-----------------------------------------------------------------------------* ! != parameters: =* ! != idate - integer, specification of the date, from yymmdd (i)=* ! != esrm2 - real, variation of the earth-sun distance (o)=* ! != esrm2 = (average e/s dist)^2 / (e/s dist on day idate)^2 =* ! !-----------------------------------------------------------------------------* ! != edit history: =* ! != 01/95 changed computation of trig function values =* ! !-----------------------------------------------------------------------------* ! != this program is free software; you can redistribute it and/or modify =* ! != it under the terms of the gnu general public license as published by the =* ! != free software foundation; either version 2 of the license, or (at your =* ! != option) any later version. =* ! != the tuv package is distributed in the hope that it will be useful, but =* ! != without any warranty; without even the implied warranty of merchantibi- =* ! != lity or fitness for a particular purpose. see the gnu general public =* ! != license for more details. =* ! != to obtain a copy of the gnu general public license, write to: =* ! != free software foundation, inc., 675 mass ave, cambridge, ma 02139, usa. =* ! !-----------------------------------------------------------------------------* ! != to contact the authors, please mail to: =* ! != sasha madronich, ncar/acd, p.o.box 3000, boulder, co, 80307-3000, usa or =* ! != send email to: sasha@ucar.edu =* ! !-----------------------------------------------------------------------------* ! != copyright (c) 1994,95,96 university corporation for atmospheric research =* ! !-----------------------------------------------------------------------------* ! ! implicit none ! integer, intent(in) ::iday,imonth ! ! ! internal: ! integer mday, month, jday ! real dayn, thet0 ! real sinth, costh, sin2th, cos2th ! ! integer,dimension(12) ::imn=(/ 31,28,31,30,31,30,31,31,30,31,30,31 /) ! ! ! start ! ! ! parse date to find day number (julian day) ! ! mday = 0 ! do month = 1, imonth-1 ! mday = mday + imn(month) ! end do ! jday = mday + iday ! dayn = float(jday - 1) + 0.5 ! ! ! define angular day number and compute esrm2: ! ! thet0 = 2.*pi*dayn/365. ! ! ! calculate sin(2*thet0), cos(2*thet0) from ! ! addition theoremes for trig functions for better ! ! performance; the computation of sin2th, cos2th ! ! is about 5-6 times faster than the evaluation ! ! of the intrinsic functions sin and cos ! ! ! sinth = sin(thet0) ! costh = cos(thet0) ! sin2th = 2.*sinth*costh ! cos2th = costh*costh - sinth*sinth ! sundis = 1.000110 + 0.034221*costh + 0.001280*sinth + & ! 0.000719*cos2th + 0.000077*sin2th ! ! ! ! end function sundis ! ! ! ! subroutine ozone_info(region) ! ! --------------------------------------------------------------------- ! ! subroutine to calculate, from an ozone field, ! ! the overhead ozone at all grid points. ! ! upper values are given by the ozone climatology above the highest ! ! model mid-level ! ! MK 01/2003 adapted for new TM5 memory structure. ! ! --------------------------------------------------------------------- ! use binas, only: Avog ! use dims ! use photolysis_data, only: phot_dat ! use global_data, only: region_dat, mass_dat ! use chem_param, only: xmo3,io3 !#ifdef MPI ! use mpi_const, only: tracer_active, tracer_loc, my_real ! use mpi_const, only: localComm, ierr, proc_tracer ! use mpi_comm, only: scatter_after_read_k !#endif ! implicit none ! ! ! input ! integer, intent(in) :: region ! ! ! local ! real,dimension(:,:,:),allocatable :: o3_overhead ! in #cm-2 ! real,dimension(:,:),pointer :: o3_overhead_surf !in #cm-2 ! real,dimension(:),pointer :: ozone_klimtop !in #cm-2 ! real,dimension(:,:,:),pointer :: o3 !in kg/gridbox ! real,dimension(:),pointer :: area !in m2 ! integer :: i,j,l ! real,parameter :: conv = 0.1*Avog/xmo3 ! ! ! conv: from kg/m2 --> #/cm2 ! ! ! start ! ! o3_overhead_surf => phot_dat(region)%vo3_surf ! allocate(o3_overhead(im(region),jm(region),lm(region))) ! !#ifdef MPI ! if (tracer_active(io3) ) then ! o3 => mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lm(region),tracer_loc(io3)) !#else ! o3 => mass_dat(region)%rm_t(1:im(region),1:jm(region),1:lm(region),io3) !#endif ! ozone_klimtop => phot_dat(region)%o3klim_top ! area => region_dat(region)%dxyp ! ! do j=1,jm(region) ! do i=1,im(region) ! o3_overhead(i,j,lm(region)) = ozone_klimtop(j) ! end do ! end do ! do l=lm(region)-1,1,-1 ! do j=1,jm(region) ! do i=1,im(region) ! o3_overhead(i,j,l) = o3_overhead(i,j,l+1) + & ! 0.5*conv*(o3(i,j,l)+o3(i,j,l+1))/area(j) ! end do ! end do ! end do ! ! do j=1,jm(region) ! do i=1,im(region) ! o3_overhead_surf(i,j) = o3_overhead(i,j,1) + & ! 0.5*conv*o3(i,j,1)/area(j) ! end do ! end do ! !#ifdef MPI ! end if ! tracer_active(io3) !#endif ! !#ifdef MPI ! call mpi_bcast(o3_overhead_surf, im(region)*jm(region), & ! my_real, proc_tracer(io3) , localComm, ierr) ! call scatter_after_read_k(o3_overhead,& ! im(region),jm(region),lm(region),0,0,0,1,& ! phot_dat(region)%vo3,proc_tracer(io3)) !#else ! phot_dat(region)%vo3 = o3_overhead !#endif ! ! phot_dat(region)%vo3s_av = phot_dat(region)%vo3s_av + o3_overhead_surf ! phot_dat(region)%nvo3s_av = phot_dat(region)%nvo3s_av + 1 ! !#ifdef MPI ! if (tracer_active(io3) ) then !#endif ! nullify(o3) ! nullify(ozone_klimtop) ! nullify(area) !#ifdef MPI ! end if !#endif ! deallocate (o3_overhead) ! nullify(o3_overhead_surf) ! ! end subroutine ozone_info ! ! ! ! subroutine cloud_info(region) ! ! ! ! routine to convert lwc, iwc, and cf to tauc and taus, and cfr ! ! to be called each time fresh cloud info becomes available ! ! ADAPTED for the new data structure, MK, 01/2003 ! ! ! use binas, only : xmair, Avog ! use photolysis_data, only : phot_dat ! use global_data, only : cloud_dat, mass_dat !#ifdef MPI ! use mpi_const, only : lmar, myid, lmloc !#endif ! implicit none ! ! ! input ! integer, intent(in) :: region ! ! ! local ! real,dimension(:,:,:),pointer :: cloudliqw,cloudicew,cloudcov ! real,dimension(:,:,:),pointer :: pressure ! real,dimension(:,:,:),allocatable :: tau_c ! ! tau_c: optical thickness clouds above level ! real,dimension(:,:), pointer :: tau_surf ! ! tau_surf: optical thickness clouds at surface ! real,dimension(:,:), pointer :: cloud_frac !cloud fraction at the surface ! real,parameter :: reff = 8.0 !fd effective cloud radius ! real,parameter :: factor = 7.24e16*1e6*xmair*29.2605/Avog ! real,dimension(:,:,:),allocatable :: wc !ice + lwc ! real :: airp,rhodz,maxtau ! integer :: i,j,l,lsave,offsetl ! ! ! start ! ! ! calculates cloud optical thickness and cloud fraction ! ! each time a new cloud field is read in (6 hours) ! ! cloudliqw => cloud_dat(region)%lwc ! cloudicew => cloud_dat(region)%iwc ! cloudcov => cloud_dat(region)%cc ! pressure => mass_dat(region)%phlb_t ! tau_surf => phot_dat(region)%taus ! cloud_frac => phot_dat(region)%cloudfrac ! ! allocate(tau_c(im(region),jm(region),lm(region))) ! allocate(wc(im(region),jm(region),lm(region))) ! ! wc = cloudliqw + cloudicew ! no differentiation between water and ice.... ! ! do l=1,lm(region) ! do j=1,jm(region) ! do i=1,im(region) ! ! threshold to account for the cloud ! if ( cloudcov(i,j,l) > 0.02) then ! airp = (pressure(i,j,l)+pressure(i,j,l+1))/2. ! ! aird = 7.24e16*airp/t(i,j,k) !#/cm3 ! ! rho = aird*1e6*xmair/avo !g/m3 ! ! dz = 29.2605*t(i,j,k)* alog(phlb(i,j,l)/(1.0+phlb(i,j,l+1))) ! rhodz = factor*airp*alog(pressure(i,j,l)/(1.0+pressure(i,j,l+1))) ! !CMK for security ! tau_c(i,j,l) = max(0.0,wc(i,j,l)*rhodz*1.5/reff) ! else ! tau_c(i,j,l) = 0.0 ! end if ! end do ! end do ! end do ! ! ! ! use optical depth to find maximum cloudcover ! !do j=1,jm(region) ! !do i=1,im(region) ! ! cfr(i,j) = maxval(cf(i,j,:)) !maximum cloud cover in the column ! !end do ! !end do ! !mktest use cloudcover from box with maximum tau ! ! ! do j=1,jm(region) ! do i=1,im(region) ! maxtau = 0.0 ! lsave = -1 ! do l=1,lm(region) ! if(tau_c(i,j,l) > maxtau) then ! maxtau = tau_c(i,j,l) ! lsave = l ! end if ! end do ! if (lsave.ne.-1) then ! cloud_frac(i,j) = cloudcov(i,j,lsave) ! else ! cloud_frac(i,j) = 0.0 ! end if ! end do ! end do ! ! ! cloud_frac=max(1e-7,cloud_frac)! avoid division by zero ! do l=1,lm(region) ! do j=1,jm(region) ! do i=1,im(region) ! !scale to cloudy fraction (temp in wc) ! wc(i,j,l) = tau_c(i,j,l)/cloud_frac(i,j) ! end do ! end do ! end do ! ! ! for radiative transfer: calculate clouds above the level!! ! ! do j=1,jm(region) ! do i=1,im(region) ! tau_c(i,j,lm(region)) = 0.5*wc(i,j,lm(region)) ! end do ! end do ! ! !cmk check.....should air mass be used? ! do l=lm(region)-1,1,-1 ! do j=1,jm(region) ! do i=1,im(region) ! tau_c(i,j,l) = tau_c(i,j,l+1) + 0.5*(wc(i,j,l)+wc(i,j,l+1)) ! end do ! end do ! end do ! ! ! ! cloud optical thickness above surface ! ! do j=1,jm(region) ! do i=1,im(region) ! tau_surf(i,j) = tau_c(i,j,1) + 0.5*wc(i,j,1) ! end do ! end do ! ! tau_c = max(tau_c,0.0) ! tau_surf = max(tau_surf,0.0) ! ! !...copy tau_c to tau_c which is parallel over the levels (for photolysis) ! offsetl=0 !#ifdef MPI ! if(myid>0)offsetl=sum(lmar(0:myid-1) ) ! offset for global value of l !#endif !#ifdef MPI ! phot_dat(region)%tauc = tau_c(:,:,offsetl+1: offsetl+lmloc) !#else ! phot_dat(region)%tauc = tau_c(:,:,1:lm(region)) !#endif ! ! ! ! phot_dat(region)%taus_av = phot_dat(region)%taus_av + tau_surf ! phot_dat(region)%cfr_av = phot_dat(region)%cfr_av + cloud_frac ! phot_dat(region)%ncloud_av = phot_dat(region)%ncloud_av + 1 ! ! ! deallocate(wc) ! deallocate(tau_c) ! ! nullify(cloudliqw) ! nullify(cloudicew) ! nullify(cloudcov) ! nullify(pressure) ! nullify(tau_surf) ! nullify(cloud_frac) ! ! end subroutine cloud_info ! ! ! ! subroutine fparam(region,level,th,mu,fact) ! ! ! use dims ! use photolysis_data, only: phot_dat ! use global_data, only: region_dat, mass_dat !#ifdef MPI ! use advect_tools, only : calc_phlb_k !#endif ! ! implicit none ! ! ! input/output ! integer, intent(in) :: region ! integer, intent(in) :: level ! the level ! real,dimension(im(region),jm(region)) :: th,mu ! sza, cos(sza) ! real,dimension(im(region),jm(region),nbands),intent(out) :: fact ! ! ! local ! real,dimension(:,:,:),pointer :: phlb ! real,dimension(:,:),pointer :: ps ! surface pressure ! real,dimension(:,:,:),pointer :: vo3 ! overhead ozone at surface ! real,dimension(:,:),pointer :: vo3s ! overhead ozone at the surface ! real,dimension(:,:,:),pointer :: tauc ! optical depth ! real,dimension(:,:),pointer :: taus ! optical depth at surface ! real,dimension(:,:),pointer :: cfr ! cloud fraction in column ! real,dimension(:,:),pointer :: ags ! surface albedo of column ! integer,dimension(:,:),pointer :: zoomed ! ! real :: tu,tuc,tl,tlc !optical thicknesses lower/upper ! real :: etu,etuc,etl,etlc !exp(tau/mu) ! real :: eto3 !for ozon or tau ! real :: ala0,alad ! real :: ala0c,aladc,al0,ald,al0c,aldc ! real :: eto3u ! real :: au,au0,au0c,auc,ag ! real :: f0,f0c,fs,fsc ! real :: fu,fuc,fd,fdc ! real :: cref ! real,parameter :: one_g = 1.-0.87 ! real,parameter :: one_g2 = 1.-0.87*0.87 ! integer :: rd=2 ! ! integer,dimension(4),parameter :: band= (/ 3,4,5,7 /) ! real,dimension(3:7),parameter :: taur= & ! (/ 1.177e-5,1.067e-5,0.919e-5,0.498e-5,0.0782e-5 /) ! ! integer :: i,j,ilon,ilat,indi,ip,ith ,ipres,k,jb ! real :: o3c_surf,mu0,pp,dp,dth,tth,taui,vo3i,vo2,p ! ! ! start ! !#ifdef MPI ! phlb => mass_dat(region)%phlb_k !#else ! phlb => mass_dat(region)%phlb_t !#endif ! ps => mass_dat(region)%p ! vo3 => phot_dat(region)%vo3 ! vo3s => phot_dat(region)%vo3_surf ! tauc => phot_dat(region)%tauc ! taus => phot_dat(region)%taus ! cfr => phot_dat(region)%cloudfrac ! ags => phot_dat(region)%albedo ! zoomed => region_dat(region)%zoomed ! !#ifdef MPI ! call calc_phlb_k(region) !#endif ! ! jb = 1 ! first Wavelength band ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5 !pressure ! if(th(i,j).lt.90.0.and.p.lt.2e4) then ! ! vo2 = p*sp ! !ozone and oxygen ! eto3 = exp(-(crs_o2*vo2 + crs_o3(jb)*vo3(i,j,level))/mu(i,j)) ! ! pp = max(102-1e-3*p,0.0) ! ip = int(pp) ! dp = pp-ip ! tth = min(th(i,j),89.9999)*0.2 ! ith = int(tth) ! dth = tth-ith ! cref =( (1-dp)*((1-dth)*fint(ip ,ith ,jb) & ! + dth *fint(ip ,ith+1,jb)) & ! + dp *((1-dth)*fint(ip+1,ith ,jb) & ! + dth *fint(ip+1,ith+1,jb)) ) ! ! ! actinic flux only depends on ozone absorption and ozygen absorption ! ! fact(i,j,jb) = cref*eto3 ! ! else ! ! fact(i,j,jb) = 0.0 ! ! end if ! ! end do ! end do ! ! jb=2 ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! ! p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5 !pressure ! if(th(i,j).lt.90.0.and.p.lt.2e4) then ! ! ! eto3 = exp(-crs_o3(jb)*vo3(i,j,level)/mu(i,j)) ! ! pp = max(102-1e-3*p,0.0) ! ip = int(pp) ! dp = pp-ip ! tth = min(th(i,j),89.9999)*0.2 ! ith = int(tth) ! dth = tth-ith ! cref =( (1-dp)*((1-dth)*fint(ip ,ith ,jb) & ! + dth *fint(ip ,ith+1,jb)) & ! + dp *((1-dth)*fint(ip+1,ith ,jb) & ! + dth *fint(ip+1,ith+1,jb)) ) ! ! ! actinic flux only depends on ozone absorption ! ! fact(i,j,jb) = cref*eto3 ! ! else ! ! fact(i,j,jb) = 0.0 ! ! end if ! end do ! end do !_____jb = 2 ! ! do k=1,4 !over wavelength intervals ! jb = band(k) ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! ! if(th(i,j).lt.90.0) then ! ! ! taui = tauc(i,j,level) ! vo3i = vo3(i,j,level) ! mu0 = mu(i,j) ! ag = ags(i,j) ! p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5 !pressure ! tu = taur(jb)*p !optical thickness upper atmosphere clear ! tuc = tu + one_g*taui ! tl = taur(jb)*(ps(i,j) - p) !lower atmosphere ! tlc = tl + one_g*(taus(i,j)-taui) !clouds included ! eto3 = exp(-crs_o3(jb)*vo3i/mu0) ! eto3u = exp(-crs_o3(jb)*(vo3s(i,j)-vo3i)*2.5) ! etu = exp(-tu/mu0) ! etuc = etu*exp(-taui*one_g2/mu0) !includes cloud effects ! etl = exp(-tl/mu0) ! etlc = etl*exp(-(taus(i,j)-taui)*one_g2/mu0) ! ala0 = 1.-(2.+3.*mu0+(2.-3.*mu0)*etl)/(4.+3.*tl) ! alad = 1.-4./(4.+3.*tl) ! ala0c = 1.-(2.+3.*mu0+(2.-3.*mu0)*etlc)/(4.+3.*tlc) ! aladc = 1.-4./(4.+3.*tlc) ! al0 = ala0 + (1.-ala0) *ag* (1.-alad) /(1.-ag*alad) ! ald = alad + (1.-alad) *ag* (1.-alad) /(1.-ag*alad) ! al0c = ala0c + (1.-ala0c) *ag* (1.-aladc) /(1.-ag*aladc) ! aldc = aladc + (1.-aladc) *ag* (1.-aladc) /(1.-ag*aladc) ! au = 1.-4./(4.+3*tu) ! au0 = 1.-(2.+3.*mu0+(2.-3.*mu0)*etu) /(4.+3.*tu) ! au0c = 1.-(2.+3.*mu0+(2.-3.*mu0)*etuc)/(4.+3.*tuc) ! auc = 1.-4./(4.+3*tuc) ! ! ! ! direct radiation ! ! ! f0 = eto3*etu ! f0c = eto3*etuc ! ! ! ! downward diffuse, only from events above p(i) ! ! ! ! fs = rd *mu0*(1.-au0-etu) *eto3 ! fsc = rd *mu0*(1.-au0c-etuc)*eto3 ! ! ! !correct lower albedo for ozone absorbtion ! ! ! ! al0 = al0 *eto3u ! ald = ald *eto3u ! al0c = al0c *eto3u ! aldc = aldc *eto3u ! ! ! ! upward diffuse radiation ! ! ! fu = (rd *mu0*f0 *al0 + ald *fs) /(1.-ald *au) ! fuc = (rd *mu0*f0c *al0c + aldc *fsc) /(1.-aldc *auc) ! ! ! ! downward direct + diffuse radiation ! ! ! fd = f0 + fs + au *fu ! fdc = f0c + fsc + auc *fuc ! ! ! ! interpolate actinic flux table ! ! ! pp = max(102-1e-3*p,0.0) ! ip = int(pp) ! dp = pp-ip ! tth = min(th(i,j),89.9999)*0.2 ! ith = int(tth) ! dth = tth-ith ! cref =( (1-dp)*((1-dth)*fint(ip ,ith ,jb) & ! + dth *fint(ip ,ith+1,jb)) & ! + dp *((1-dth)*fint(ip+1,ith ,jb) & ! + dth *fint(ip+1,ith+1,jb)) ) ! ! ! ! calculate actinic flux, weighted with the cloud fraction ! ! ! fact(i,j,jb) = (cfr(i,j)*(fdc + fuc) + & ! (1.-cfr(i,j))*(fu+fd)) *cref ! else ! ! fact(i,j,jb) = 0.0 ! ! end if ! end do ! end do ! end do ! ! jb = 6 !no ozone absorption ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! ! p = (phlb(i,j,level+1) + phlb(i,j,level))*0.5 !pressure ! if(th(i,j).lt.90.0) then ! ! taui = tauc(i,j,level) ! mu0 = mu(i,j) ! ag = ags(i,j) ! tu = taur(jb)*p !optical thickness upper atmosphere clear ! tuc = tu + one_g*taui ! tl = taur(jb)*(ps(i,j) - p) !lower atmosphere ! tlc = tl + one_g*(taus(i,j)-taui) !clouds included ! etu = exp(-tu/mu0) ! etuc = etu*exp(-taui*one_g2/mu0) !includes cloud effects ! etl = exp(-tl/mu0) ! etlc = etl*exp(-(taus(i,j)-taui)*one_g2/mu0) ! ala0 = 1.-(2.+3.*mu0+(2.-3.*mu0)*etl)/(4.+3.*tl) ! alad = 1.-4./(4.+3.*tl) ! ala0c = 1.-(2.+3.*mu0+(2.-3.*mu0)*etlc)/(4.+3.*tlc) ! aladc = 1.-4./(4.+3.*tlc) ! al0 = ala0 + (1.-ala0) *ag* (1.-alad) /(1.-ag*alad) ! ald = alad + (1.-alad) *ag* (1.-alad) /(1.-ag*alad) ! al0c = ala0c + (1.-ala0c) *ag* (1.-aladc) /(1.-ag*aladc) ! aldc = aladc + (1.-aladc) *ag* (1.-aladc) /(1.-ag*aladc) ! au = 1.-4./(4.+3*tu) ! au0 = 1.-(2.+3.*mu0+(2.-3.*mu0)*etu) /(4.+3.*tu) ! au0c = 1.-(2.+3.*mu0+(2.-3.*mu0)*etuc)/(4.+3.*tuc) ! auc = 1.-4./(4.+3*tuc) ! ! ! ! direct radiation ! ! ! f0 = etu ! f0c = etuc ! ! ! ! downward diffuse, only from events above p(i) ! ! ! fs = rd *mu0*(1.-au0-etu) ! fsc = rd *mu0*(1.-au0c-etuc) ! ! ! ! upward diffuse radiation ! ! ! fu = (rd *mu0*f0 *al0 + ald *fs) /(1.-ald *au) ! fuc = (rd *mu0*f0c *al0c + aldc *fsc) /(1.-aldc *auc) ! ! ! ! downward direct + diffuse radiation ! ! ! fd = f0 + fs + au *fu ! fdc = f0c + fsc + auc *fuc ! ! ! ! interpolate actinic flux table ! ! ! pp = max(102-1e-3*p,0.0) ! ip = int(pp) ! dp = pp-ip ! tth = min(th(i,j),89.9999)*0.2 ! ith = int(tth) ! dth = tth-ith ! cref =( (1-dp)*((1-dth)*fint(ip ,ith ,jb) & ! + dth *fint(ip ,ith+1,jb)) & ! + dp *((1-dth)*fint(ip+1,ith ,jb) & ! + dth *fint(ip+1,ith+1,jb)) ) ! ! ! ! calculate actinic flux, weighted with the cloud fraction ! ! ! fact(i,j,jb) = (cfr(i,j)*(fdc + fuc) + (1.-cfr(i,j))*(fu+fd)) *cref ! !cmk fact(i,j) = (fu+fd) *cref ! ! else ! ! fact(i,j,jb) = 0.0 ! ! end if ! ! end do ! end do !___jb=6 ! ! ! nullify(phlb) ! nullify(ps) ! nullify(vo3) ! nullify(vo3s) ! nullify(tauc) ! nullify(taus) ! nullify(cfr) ! nullify(ags) ! nullify(zoomed) ! ! end subroutine fparam ! ! ! ! subroutine photolysis_rates(region,level,th,mu,fact,rjx) ! ! ! use dims !#ifdef MPI ! use mpi_const, only : lmar, myid, lmloc !#endif ! use global_data, only: mass_dat, region_dat, meteo_dat ! use photolysis_data, only: nj ,jo3d,jno2,jh2o2,jhno3,jhno4,jn2o5,jach2o, & ! jbch2o,jmepe,jano3,jbno3,jpan,jorgn,j45,j74,jrooh, phot_dat ! ! implicit none ! ! ! input/output ! !region & model level ! integer,intent(in) :: region ! integer,intent(in) :: level ! ! fact: actinic flux calculated by fparam... ! real,dimension(im(region),jm(region),nbands),intent(in) :: fact ! ! th: sza (deg) mu: cos(sza) ! real,dimension(im(region),jm(region)),intent(in) :: th ! real,dimension(im(region),jm(region)),intent(in) :: mu ! ! rjx: photolysis rates ! real,dimension(im(region),jm(region),nj ),intent(out) :: rjx ! ! ! local ! real,dimension(:,:,:), pointer :: phlb ! integer,dimension(:,:) , pointer :: zoomed ! ! t: temprature (K) NOTE: this if in parallel full 3D array ! real,dimension(:,:,:), pointer :: t ! ! vo3: overhead ozone (#/cm2) ! real,dimension(:,:,:), pointer :: vo3 ! ! ! ratio incoming f0 (band)/ incoming f0(li) ! ! depends on extraterrestial flux chosen ! ! ! test works with fint2.data real,dimension(1:7),parameter ::f0 = & ! ! (/ 70.61, 4.850, 19.81, 9.303, 24.86, 18.28, 48.04/) ! ! f0: bruhl wavelength works with fintwb.data ! real,dimension(nbands),parameter :: f0 = & ! (/ 73.44, 4.977, 16.13, 8.017, 27.9, 19.00, 48.12/) ! ! real, parameter :: densref = 1.3e19 ! ! real :: dens,vo2 ! real :: tau_2,tau_3,tau_4,tau_5,tau_7 !optical depth ! real :: tau_1,sig_1,fint_1 ! real :: sig_2,sig_3,sig_4,sig_5,sig_7,sig_6 !effective cross-sections ! real :: fint_2,fint_3,fint_4 !integrated actinic fluxes ! real :: fint_5,fint_6,fint_7 !integrated actinic fluxes ! real :: vo3s,v3_du2,v3_du1,v3_du,v3s1,v3s2 !ozone column variables ! real :: vo2s,v2s_m,v2s2 !oxygen column variables ! real :: tscale,tscale2,pscale !temperature and pressure scales ! real :: p2s,p3s,p4,p5,p6,p7 !tempvar for crosssection calculations ! real :: cc0,cc1,cc2 ! real :: pres,esrm2 ! integer :: i1,i2,i3,ind !indices lookup tables ! integer :: i,j !count ! integer :: ilon,ilat,ipres,offsetl ! ! include 'photolysis.i' ! ! ! start ! !#ifdef MPI ! phlb => mass_dat(region)%phlb_k !#else ! phlb => mass_dat(region)%phlb_t !#endif ! zoomed => region_dat(region)%zoomed ! ! offsetl=0 !#ifdef MPI ! if(myid>0)offsetl=sum(lmar(0:myid-1) ) ! offset for global value of l !#endif ! !#ifdef MPI ! t => meteo_dat(region)%t(1:im(region),1:jm(region),offsetl+1:offsetl+lmloc) !#else ! t => meteo_dat(region)%t(1:im(region),1:jm(region),1:lm(region)) !#endif ! !print *, 'CHECK in photolysis:', myid, lbound(t),ubound(t) ! vo3 => phot_dat(region)%vo3 ! ! rjx = 0.0 ! ! do j=jsr(region),jer(region) ! do i=isr(region),ier(region) ! if(zoomed(i,j)/=region) cycle ! ! if(th(i,j).lt.90.) then ! ! pres = 0.5*(phlb(i,j,level+1)+phlb(i,j,level)) ! vo2= pres*sp !make oxygen conc. in #/cm3 ! vo3s = vo3(i,j,level)/mu(i,j) ! ! vo2s = vo2/mu(i,j) ! v3_du = vo3s*todu ! v2s_m = vo2s*to_m ! ! ! interval: 1 ! ! if (v2s_m.ge.500.) then ! v2s_m = 500. ! v2s2 = 1.3602e+24 ! else ! v2s2 = vo2s ! end if ! ! ! interval: 1+2 ! if (v3_du.ge.300.) then ! v3_du1 = 300. ! v3s1 = 8.161e+18 ! else ! v3_du1 = v3_du ! v3s1 = vo3s ! end if ! ! ! interval: 3 - 7 ! if (v3_du.ge.3000.) then ! v3_du2 = 3000. ! v3s2 = 8.161e+19 ! else ! v3_du2 = v3_du ! v3s2 = vo3s ! end if ! ! ! get indices for lookup-tables ! ind = min(int((v3_du-0.5)/2.5) +1,115) ! i2 = ifil(ind) ! i1 = i2 ! ind = min(int((v3_du-5)/25) +1,115) ! i3 = ifil(ind) ! ! !---calculate optical depth due to ozone absorbtion ! ! tau_1 = p1(taub1(i1,1),taua1(i1,1),v3_du1) + & ! p1(taub1(i1,2),taua1(i1,2),v3_du1) * v2s_m + & ! p1(taub1(i1,3),taua1(i1,3),v3_du1) * v2s_m**2 ! tau_2 = p1(taub2(i2),taua2(i2),v3_du1) ! tau_3 = p1(taub3(i3),taua3(i3),v3_du2) ! tau_4 = p3(c4_tau(1),c4_tau(2),c4_tau(3),c4_tau(4),v3_du2) ! tau_5 = p2(c5_tau(1),c5_tau(2),c5_tau(3),v3_du2) ! tau_7 = p1(c7_tau(1),c7_tau(2),v3_du2) ! ! !---get effective actinic fluxes in units quanta/s/cm2 ! ! fint_1 = exp(-tau_1 + crs_o3(1)*v3s1 + & ! crs_o2 *v2s2)*f0(1)*fact(i,j,1)*fe(1) ! fint_2 = exp(-tau_2 + crs_o3(2)*v3s1)*f0(2)*fact(i,j,2)*fe(2) ! fint_3 = exp(-tau_3 + crs_o3(3)*v3s2)*f0(3)*fact(i,j,3)*fe(3) ! fint_4 = exp(-tau_4 + crs_o3(4)*v3s2)*f0(4)*fact(i,j,4)*fe(4) ! fint_5 = exp(-tau_5 + crs_o3(5)*v3s2)*f0(5)*fact(i,j,5)*fe(5) ! fint_6 = f0(6)*fact(i,j,6)*fe(6) ! fint_7 = exp(-tau_7 + crs_o3(7)*v3s2)*f0(7)*fact(i,j,7)*fe(7) ! ! tscale = (t(i,j,level)-250.)/250. ! dens = 7.24e16*pres/t(i,j,level) ! pscale = (dens-densref)/densref ! ! !-------------------- ! ! o3 -> o(1d) + o2 ! !-------------------- ! p2s = p1(b2_o1d(i2),a2_o1d(i2),v3_du1) ! p3s = p1(b3_o1d(i3),a3_o1d(i3),v3_du2) ! p4 = p3(c4_o1d(1),c4_o1d(2),c4_o1d(3),c4_o1d(4),v3_du2) ! p5 = p3(c5_o1d(1),c5_o1d(2),c5_o1d(3),c5_o1d(4),v3_du2) ! ! sig_2 = p3(tj_o1d(1,2),tj_o1d(2,2),tj_o1d(3,2),tj_o1d(4,2),tscale) & ! * p2s ! sig_3 = p3(tj_o1d(1,3),tj_o1d(2,3),tj_o1d(3,3),tj_o1d(4,3),tscale) & ! * p3s ! sig_4 = p3(tj_o1d(1,4),tj_o1d(2,4),tj_o1d(3,4),tj_o1d(4,4),tscale) & ! * p4 ! sig_5 = p3(tj_o1d(1,5),tj_o1d(2,5),tj_o1d(3,5),tj_o1d(4,5),tscale) & ! * p5 ! ! rjx(i,j,jo3d) = sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 ! !-------------------- ! ! no2 ! !-------------------- ! ! sig_2 =p1(tj_no2(1,2),tj_no2(2,2),tscale) * & ! p1(b2_no2(i2),a2_no2(i2),v3_du1) ! sig_3 =p1(tj_no2(1,3),tj_no2(2,3),tscale) * & ! p1(b3_no2(i3),a3_no2(i3),v3_du2) ! sig_4 =p1(tj_no2(1,4),tj_no2(2,4),tscale) * & ! p2(c4_no2(1) ,c4_no2(2),c4_no2(3),v3_du2) ! sig_5 =p1(tj_no2(1,5),tj_no2(2,5),tscale) * & ! p2(c5_no2(1) ,c5_no2(2),c5_no2(3),v3_du2) ! sig_6 =p1(tj_no2(1,6),tj_no2(2,6),tscale) * & ! p1(c6_no2(1) ,c6_no2(2),v3_du2) ! ! rjx(i,j,jno2) = sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! !-------------------- ! ! h2o2 ! !-------------------- ! ! sig_1 = p1(b1_h2o2(i1,1),a1_h2o2(i1,1),v3_du1) + & ! p1(b1_h2o2(i1,2),a1_h2o2(i1,2),v3_du1)* v2s_m ! sig_2 = p2(tj_h2o2(1,2),tj_h2o2(2,2),tj_h2o2(3,2),tscale) * & ! p1(b2_h2o2(i2),a2_h2o2(i2),v3_du1) ! sig_3 = p2(tj_h2o2(1,3),tj_h2o2(2,3),tj_h2o2(3,3),tscale) * & ! p1(b3_h2o2(i3),a3_h2o2(i3),v3_du2) ! sig_4 = p2(tj_h2o2(1,4),tj_h2o2(2,4),tj_h2o2(3,4),tscale) * & ! p2(c4_h2o2(1),c4_h2o2(2),c4_h2o2(3),v3_du2) ! sig_5 = p2(tj_h2o2(1,5),tj_h2o2(2,5),tj_h2o2(3,5),tscale) * & ! p2(c5_h2o2(1),c5_h2o2(2),c5_h2o2(3),v3_du2) ! sig_6 = p2(tj_h2o2(1,6),tj_h2o2(2,6),tj_h2o2(3,6),tscale) * & ! p1(c6_h2o2(1),c6_h2o2(2),v3_du2) ! ! rjx(i,j,jh2o2) = sig_1 * fint_1 & ! + sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! !-------------------- ! ! hno3 ! !-------------------- ! ! sig_1 =(p1(b1_hno3(i1,1),a1_hno3(i1,1), v3_du1) + & ! p1(b1_hno3(i1,2),a1_hno3(i1,2), v3_du1)*v2s_m)* & ! p3(tj_hno3(1,1),tj_hno3(2,1),tj_hno3(3,1),tj_hno3(4,1),tscale) ! sig_2 = & ! p3(tj_hno3(1,2),tj_hno3(2,2),tj_hno3(3,2),tj_hno3(4,2),tscale) & ! * p1(b2_hno3(i2),a2_hno3(i2),v3_du1) ! sig_3 = & ! p3(tj_hno3(1,3),tj_hno3(2,3),tj_hno3(3,3),tj_hno3(4,3),tscale) & ! * p1(b3_hno3(i3),a3_hno3(i3),v3_du2) ! sig_4 = & ! p3(tj_hno3(1,4),tj_hno3(2,4),tj_hno3(3,4),tj_hno3(4,4),tscale) & ! * p2(c4_hno3(1),c4_hno3(2),c4_hno3(3),v3_du2) ! sig_5 = & ! p3(tj_hno3(1,5),tj_hno3(2,5),tj_hno3(3,5),tj_hno3(4,5),tscale) & ! * p2(c5_hno3(1),c5_hno3(2),c5_hno3(3),v3_du2) ! sig_6 = & ! p3(tj_hno3(1,6),tj_hno3(2,6),tj_hno3(3,6),tj_hno3(4,6),tscale) & ! * p1(c6_hno3(1),c6_hno3(2),v3_du2) ! ! rjx(i,j,jhno3) = sig_1 * fint_1 & ! + sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! !-------------------- ! ! hno4 ! !-------------------- ! ! sig_1 = p1(b1_hno4(i1,1),a1_hno4(i1,1), v3_du1) + & ! p1(b1_hno4(i1,2),a1_hno4(i1,2), v3_du1) *v2s_m ! sig_2 = p1(b2_hno4(i2), a2_hno4(i2),v3_du1) ! sig_3 = p1(b3_hno4(i3), a3_hno4(i3),v3_du2) ! sig_4 = p2(c4_hno4(1),c4_hno4(2),c4_hno4(3),v3_du2) ! sig_5 = p2(c5_hno4(1),c5_hno4(2),c5_hno4(3),v3_du2) ! ! rjx(i,j,jhno4) = sig_1 * fint_1 & ! + sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 ! ! !-------------------- ! ! n2o5 ! !-------------------- ! ! tscale2 = (max(min(t(i,j,level),300.0),225.0)-250.0)/250.0 ! sig_1 = p1(b1_n2o5(i1,1),a1_n2o5(i1,1), v3_du1) + & ! p1(b1_n2o5(i1,2),a1_n2o5(i1,2), v3_du1) *v2s_m ! sig_2 = p3(tj_n2o5(1,2),tj_n2o5(2,2),tj_n2o5(3,2), & ! tj_n2o5(4,2),tscale2) * p1(b2_n2o5(i2),a2_n2o5(i2),v3_du1) ! sig_3 = p3(tj_n2o5(1,3),tj_n2o5(2,3),tj_n2o5(3,3), & ! tj_n2o5(4,3),tscale2) * p1(b3_n2o5(i3),a3_n2o5(i3),v3_du2) ! sig_4 = p3(tj_n2o5(1,4),tj_n2o5(2,4),tj_n2o5(3,4), & ! tj_n2o5(4,4),tscale2) * & ! p2(c4_n2o5(1),c4_n2o5(2),c4_n2o5(3),v3_du2) ! sig_5 = p3(tj_n2o5(1,5),tj_n2o5(2,5),tj_n2o5(3,5), & ! tj_n2o5(4,5),tscale2) * & ! p2(c5_n2o5(1),c5_n2o5(2),c5_n2o5(3),v3_du2) ! sig_6 = p3(tj_n2o5(1,6),tj_n2o5(2,6),tj_n2o5(3,6), & ! tj_n2o5(4,6),tscale2) * p1(c6_n2o5(1),c6_n2o5(2),v3_du2) ! ! rjx(i,j,jn2o5) = sig_1 * fint_1 & ! + sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! !-------------------- ! ! ch2o -> co + h2 ! !-------------------- ! ! sig_3 = p1(tj_coh2(1,3),tj_coh2(2,3),tscale) * & ! p1(b3_coh2(i3),a3_coh2(i3),v3_du2) ! sig_4 = p3(c4_coh2(1),c4_coh2(2),c4_coh2(3),c4_coh2(4),v3_du2) ! sig_5 = p1(tj_coh2(1,5),tj_coh2(2,5),tscale) * & ! p2(c5_coh2(1),c5_coh2(2),c5_coh2(3),v3_du2) ! cc0 = p2( 1.00997 , -0.398144 , -0.265157, tscale) ! cc1 = p2( 0.387158 , -0.434852 , -0.221940, tscale) ! cc2 = p2( -0.0388784 , 0.0571797 , 0.0876116, tscale) ! sig_6 = p1(c6_coh2(1),c6_coh2(2),v3_du2)/ p2(cc0,cc1,cc2,pscale) ! ! rjx(i,j,jach2o) = sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! !-------------------- ! ! ch2o -> cho + h ! !-------------------- ! sig_3 = p1(tj_choh(1,3),tj_choh(2,3),tscale) * & ! p1(b3_choh(i3),a3_choh(i3),v3_du2) ! sig_4 = p3(c4_choh(1),c4_choh(2),c4_choh(3), c4_choh(4),v3_du2) ! sig_5 = p1(tj_choh(1,5),tj_choh(2,5),tscale) * & ! p2(c5_choh(1),c5_choh(2),c5_choh(3),v3_du2) ! sig_6 = p1(tj_choh(1,6),tj_choh(2,6),tscale) * & ! p1(c6_choh(1),c6_choh(2),v3_du2) ! ! rjx(i,j,jbch2o) = sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! !-------------------- ! ! ch3ooh --> oh + ch3o ! !-------------------- ! sig_1 =p1(b1_ch3ooh(i1,1),a1_ch3ooh(i1,1),v3_du1) + & ! p1(b1_ch3ooh(i1,2),a1_ch3ooh(i1,2),v3_du1) * v2s_m ! sig_2 =p1(b2_ch3ooh(i2),a2_ch3ooh(i2),v3_du1) ! sig_3 =p1(b3_ch3ooh(i3),a3_ch3ooh(i3),v3_du2) ! sig_4 =p2(c4_ch3ooh(1),c4_ch3ooh(2),c4_ch3ooh(3),v3_du2) ! sig_5 =p2(c5_ch3ooh(1),c5_ch3ooh(2),c5_ch3ooh(3),v3_du2) ! sig_6 =p1(c6_ch3ooh(1),c6_ch3ooh(2),v3_du2) ! ! rjx(i,j,jmepe) = sig_1 * fint_1 & ! + sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! ! !-------------------- ! ! no3 -> no2 + o ! !-------------------- ! ! ! sig_6 = p1(tj_no2o(1,6),tj_no2o(2,6),tscale) * & ! p1(c6_no2o(1),c6_no2o(2),v3_du2) ! sig_7 = p1(tj_no2o(1,7),tj_no2o(2,7),tscale) * & ! p1(c7_no2o(1),c7_no2o(2),v3_du2) ! ! rjx(i,j,jano3) = sig_6 * fint_6 + & ! sig_7 * fint_7 ! ! !-------------------- ! ! no3 -> no + o2 ! !-------------------- ! ! sig_7 =p1(tj_noo2(1,7),tj_noo2(2,7),tscale)* & ! p1(c7_noo2(1),c7_noo2(2),v3_du2) ! ! rjx(i,j,jbno3) = sig_7 * fint_7 ! ! ! --- ! ! pan ! ! --- ! ! sig_1 = p3(tj_pan(1,1),tj_pan(2,1),tj_pan(3,1),tj_pan(4,1),tscale)* & ! (p1(b1_pan(i1,1),a1_pan(i1,1), v3_du1) + & ! p1(b1_pan(i1,2),a1_pan(i1,2), v3_du1) *v2s_m) ! sig_2 = p3(tj_pan(1,2),tj_pan(2,2),tj_pan(3,2),tj_pan(4,2),tscale)* & ! p1(b2_pan(i2),a2_pan(i2),v3_du1) ! sig_3 = p3(tj_pan(1,3),tj_pan(2,3),tj_pan(3,3),tj_pan(4,3),tscale)* & ! p1(b3_pan(i3),a3_pan(i3),v3_du2) ! sig_4 = p3(tj_pan(1,4),tj_pan(2,4),tj_pan(3,4),tj_pan(4,4),tscale)* & ! p2(c4_pan(1),c4_pan(2),c4_pan(3),v3_du2) ! sig_5 = p3(tj_pan(1,5),tj_pan(2,5),tj_pan(3,5),tj_pan(4,5),tscale)* & ! p2(c5_pan(1),c5_pan(2),c5_pan(3),v3_du2) ! sig_6 = p3(tj_pan(1,6),tj_pan(2,6),tj_pan(3,6),tj_pan(4,6),tscale)* & ! p1(c6_pan(1),c6_pan(2),v3_du2) ! ! rjx(i,j,jpan) = sig_1 * fint_1 & ! + sig_2 * fint_2 & ! + sig_3 * fint_3 & ! + sig_4 * fint_4 & ! + sig_5 * fint_5 & ! + sig_6 * fint_6 ! ! ! --- ! ! ch3ono2 ! ! --- ! ! sig_2 = p1(b2_ch3ono2(i2),a2_ch3ono2(i2),v3_du1) ! sig_3 = p1(b3_ch3ono2(i2),a3_ch3ono2(i2),v3_du2) ! sig_4 = p2(c4_ch3ono2(1), c4_ch3ono2(2), c4_ch3ono2(3), v3_du2) ! sig_5 = p2(c5_ch3ono2(1), c5_ch3ono2(2), c5_ch3ono2(3), v3_du2) ! ! rjx(i,j,jorgn) = sig_2 * fint_2 + & ! sig_3 * fint_3 + & ! sig_4 * fint_4 + & ! sig_5 * fint_5 ! ! ! --- ! ! ch3cho ! ! --- ! ! sig_2 = p1(b2_ch3cho(i2),a2_ch3cho(i2),v3_du1) & ! / p1(tp_ch3cho(1,2),tp_ch3cho(2,2),pscale) ! sig_3 = p1(b3_ch3cho(i2),a3_ch3cho(i2),v3_du2) & ! / p1(tp_ch3cho(1,3),tp_ch3cho(2,3),pscale) ! sig_4 = p2(c4_ch3cho(1), c4_ch3cho(2), c4_ch3cho(3), v3_du2) & ! / p1(tp_ch3cho(1,4),tp_ch3cho(2,4),pscale) ! sig_5 = p2(c5_ch3cho(1), c5_ch3cho(2), c5_ch3cho(3), v3_du2) & ! / p1(tp_ch3cho(1,5),tp_ch3cho(2,5),pscale) ! ! rjx(i,j,j45) = sig_2 * fint_2 + & ! sig_3 * fint_3 + & ! sig_4 * fint_4 + & ! sig_5 * fint_5 ! ! ! multiply methyl-nitrate to get organic nitrate ! rjx(i,j,jorgn) = 2.5 * rjx(i,j,jorgn) ! rjx(i,j,j74) = 50.0* rjx(i,j,jbch2o) !methylglyoxal ! rjx(i,j,jrooh) = rjx(i,j,jmepe) !other peroxides... ! end if ! ! end do ! end do ! i,j loop ! ! ! correction photolysis for dsitance sun-earth ! esrm2 = sundis(idate(2),idate(3)) ! ! multiply with seasonal correction factor ! rjx = rjx*esrm2 ! ! phot_dat(region)%jo3_av (:,:,offsetl+level) = & ! phot_dat(region)%jo3_av (:,:,offsetl+level) + rjx(:,:,jo3d) ! phot_dat(region)%jno2_av(:,:,offsetl+level) = & ! phot_dat(region)%jno2_av(:,:,offsetl+level) + rjx(:,:,jno2) !#ifdef MPI ! if ( level == lmloc ) phot_dat(region)%nj_av = & ! phot_dat(region)%nj_av + 1 !proceed only after last level.... !#else ! if ( level == lm(region) ) phot_dat(region)%nj_av = & ! phot_dat(region)%nj_av + 1 !proceed only after last level.... !#endif ! ! nullify(phlb) ! nullify(zoomed) ! nullify(t) ! nullify(vo3) ! ! contains ! ! real function p1(c0,c1,x) ! implicit none ! real, intent(in) :: c0, c1, x ! p1 = c0 + c1*x ! end function p1 ! ! real function p2(c0,c1,c2,x) ! implicit none ! real, intent(in) :: c0, c1, c2, x ! p2 = c0 + (c1 + c2*x)*x ! end function p2 ! ! real function p3(c0,c1,c2,c3,x) ! implicit none ! real, intent(in) :: c0, c1, c2, c3, x ! p3 = c0 + (c1 + (c2 + c3*x)*x)*x ! end function p3 ! ! end subroutine photolysis_rates ! ! subroutine that deallocates photolysis memory ! and writes photolysis statistics ! called from trace_end ! subroutine end_photolysis use dims , only : nregions use photolysis_data, only : phot_dat ! use global_data, only : outdir ! use io_hdf !#ifdef MPI ! use mpi_const, only : lmloc, my_real, myid, root_k, mpi_sum, com_lev, ierr !#endif ! ! ! called by trace_end ! ! output of (e.g. monthly) averaged photolysis jo3 jno2 ! ! and surface data of: albedo, cloud-cover, cloud optical thickness ! ! total ozone ! ! Deallocation of memory ! ! integer :: io,n ,sfstart,sfend,sfsattr ! integer :: k,iret,imr,jmr,lmr ! character(len=22) :: file_name ! character(len=7) :: rg integer :: region ! real,dimension(:,:,:),allocatable :: average_field ! ! file_name = trim(outdir)//'/j_statistics.hdf' !#ifdef MPI ! if (myid == root_k) io = sfstart(file_name,DFACC_CREATE) !#else ! io = sfstart(file_name,DFACC_CREATE) !#endif ! ! do region=1,nregions ! rg = '_'//trim(region_name(region)) ! imr=im(region) ! jmr=jm(region) ! lmr=lm(region) ! ! print *, 'end_photolysis: j_statistics ', & ! '(region, nj_av, nvo3s_av, ncloud_av)', region, & ! phot_dat(region)%nj_av, phot_dat(region)%nvo3s_av, & ! phot_dat(region)%ncloud_av ! if ( phot_dat(region)%nj_av > 0 ) then ! allocate(average_field(imr,jmr,lmr)) !#ifdef MPI ! if ( lmloc /= 0 ) then ! call mpi_allreduce(phot_dat(region)%jo3_av,average_field, & ! imr*jmr*lmr,my_real,mpi_sum,com_lev,ierr) ! end if !#else ! average_field = phot_dat(region)%jo3_av !#endif !#ifdef MPI ! if ( myid == root_k ) then !#endif ! average_field = average_field/phot_dat(region)%nj_av ! call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, & ! lmr,'HYBRID'//rg,average_field,'jo3_av') !#ifdef MPI ! end if !#endif !#ifdef MPI ! if ( lmloc /= 0 ) then ! call mpi_allreduce(phot_dat(region)%jno2_av,average_field, & ! imr*jmr*lmr,my_real,mpi_sum,com_lev,ierr) ! end if !#else ! average_field = phot_dat(region)%jno2_av ! !#endif !#ifdef MPI ! if ( myid == root_k ) then !#endif ! average_field = average_field/phot_dat(region)%nj_av ! call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, & ! lmr,'HYBRID'//rg,average_field,'jno2_av') !#ifdef MPI ! end if !#endif ! deallocate(average_field) ! end if !#ifdef MPI ! if(myid == root_k) then !#endif ! if (phot_dat(region)%nvo3s_av > 0) then ! phot_dat(region)%vo3s_av = & ! todu*phot_dat(region)%vo3s_av/phot_dat(region)%nvo3s_av ! call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, & ! phot_dat(region)%vo3s_av,'vo3s_av') ! end if ! if ( phot_dat(region)%nalb_av > 0 ) then ! phot_dat(region)%ags_av = & ! phot_dat(region)%ags_av/phot_dat(region)%nalb_av ! call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, & ! phot_dat(region)%ags_av,'ags_av') ! end if !#ifdef MPI ! end if !#endif ! !#ifdef MPI ! if ( myid == root_k ) then !#endif ! if ( phot_dat(region)%ncloud_av > 0 ) then ! phot_dat(region)%taus_av = & ! phot_dat(region)%taus_av/phot_dat(region)%ncloud_av ! call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, & ! phot_dat(region)%taus_av,'taus_av') ! phot_dat(region)%cfr_av = & ! phot_dat(region)%cfr_av/phot_dat(region)%ncloud_av ! call io_write(io,imr,'LON'//rg,jmr,'LAT'//rg, & ! phot_dat(region)%cfr_av,'cfr_av') ! end if !#ifdef MPI ! end if !#endif ! end do ! !#ifdef MPI ! if ( myid == root_k ) iret = sfend(io) !#else ! iret = sfend(io) !#endif do region = 1,nregions ! deallocate(phot_dat(region)%taus ) ! deallocate(phot_dat(region)%cloudfrac ) deallocate(phot_dat(region)%albedo ) ! deallocate(phot_dat(region)%vo3s_av ) ! deallocate(phot_dat(region)%taus_av ) ! deallocate(phot_dat(region)%cfr_av ) ! deallocate(phot_dat(region)%vo3 ) ! deallocate(phot_dat(region)%tauc ) ! deallocate(phot_dat(region)%jo3_av ) ! deallocate(phot_dat(region)%jno2_av ) !JEW: deallocate(phot_dat(region)%ags_av ) end do end subroutine end_photolysis end module photolysis