#ifdef AIX #define USE_ESSL #endif #define USE_BDE module mo_jlong use shr_kind_mod, only : r4 => shr_kind_r4 use shr_kind_mod, only : r8 => shr_kind_r8 use cam_logfile, only : iulog use cam_abortutils, only : endrun #ifdef SPMD use mpishorthand, only : mpicom,mpiint,mpir8, mpilog, mpir4 #endif use spmd_utils, only : masterproc implicit none interface jlong module procedure jlong_photo module procedure jlong_hrates end interface private public :: jlong_init public :: jlong_timestep_init public :: jlong public :: numj save real(r8), parameter :: hc = 6.62608e-34_r8 * 2.9979e8_r8 / 1.e-9_r8 real(r8), parameter :: wc_o2_b = 242.37_r8 ! (nm) real(r8), parameter :: wc_o3_a = 310.32_r8 ! (nm) real(r8), parameter :: wc_o3_b = 1179.87_r8 ! (nm) integer :: nw ! wavelengths >200nm integer :: nt ! number of temperatures in xsection table integer :: np_xs ! number of pressure levels in xsection table integer :: numj ! number of photorates in xsqy, rsf integer :: nump ! number of altitudes in rsf integer :: numsza ! number of zen angles in rsf integer :: numalb ! number of albedos in rsf integer :: numcolo3 ! number of o3 columns in rsf real(r4), allocatable :: xsqy(:,:,:,:) real(r8), allocatable :: wc(:) real(r8), allocatable :: we(:) real(r8), allocatable :: wlintv(:) real(r8), allocatable :: etfphot(:) real(r8), allocatable :: bde_o2_b(:) real(r8), allocatable :: bde_o3_a(:) real(r8), allocatable :: bde_o3_b(:) real(r8), allocatable :: xs_o2b(:,:,:) real(r8), allocatable :: xs_o3a(:,:,:) real(r8), allocatable :: xs_o3b(:,:,:) real(r8), allocatable :: p(:) real(r8), allocatable :: del_p(:) real(r8), allocatable :: prs(:) real(r8), allocatable :: dprs(:) real(r8), allocatable :: sza(:) real(r8), allocatable :: del_sza(:) real(r8), allocatable :: alb(:) real(r8), allocatable :: del_alb(:) real(r8), allocatable :: o3rat(:) real(r8), allocatable :: del_o3rat(:) real(r8), allocatable :: colo3(:) real(r4), allocatable :: rsf_tab(:,:,:,:,:) logical :: jlong_used = .false. contains subroutine jlong_init( xs_long_file, rsf_file, lng_indexer ) use ppgrid, only : pver use time_manager, only : is_end_curr_day use mo_util, only : rebin use solar_data, only : data_nw => nbins, data_we => we, data_etf => sol_etf implicit none !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ integer, intent(inout) :: lng_indexer(:) character(len=*), intent(in) :: xs_long_file, rsf_file !------------------------------------------------------------------------------ ! ... read Cross Section * QY NetCDF file ! find temperature index for given altitude ! derive cross*QY results, returns xsqy(nj,nz,nw) !------------------------------------------------------------------------------ call get_xsqy( xs_long_file, lng_indexer ) !------------------------------------------------------------------------------ ! ... read radiative source function NetCDF file !------------------------------------------------------------------------------ if(masterproc) write(iulog,*) 'jlong_init: before get_rsf' call get_rsf(rsf_file) if(masterproc) write(iulog,*) 'jlong_init: after get_rsf' we(:nw) = wc(:nw) - .5_r8*wlintv(:nw) we(nw+1) = wc(nw) + .5_r8*wlintv(nw) if (masterproc) then write(iulog,*) ' ' write(iulog,*) '--------------------------------------------------' endif call rebin( data_nw, nw, data_we, we, data_etf, etfphot ) if (masterproc) then write(iulog,*) 'jlong_init: etfphot after data rebin' write(iulog,'(1p,5g15.7)') etfphot(:) write(iulog,*) '--------------------------------------------------' write(iulog,*) ' ' endif jlong_used = .true. end subroutine jlong_init subroutine get_xsqy( xs_long_file, lng_indexer ) !=============================================================================! ! PURPOSE: ! ! Reads a NetCDF file that contains: ! ! cross section * QY temperature dependence, >200nm ! ! ! !=============================================================================! ! PARAMETERS: ! ! Input: ! ! filepath.... NetCDF filepath that contains the "cross sections" ! !=============================================================================! ! EDIT HISTORY: ! ! Created by Doug Kinnison, 3/14/2002 ! !=============================================================================! use ioFileMod, only : getfil use error_messages, only : alloc_err use chem_mods, only : phtcnt, pht_alias_lst, rxt_tag_lst use netcdf, only: & nf90_open, & nf90_nowrite, & nf90_inq_dimid, & nf90_inquire_dimension, & nf90_inq_varid, & nf90_noerr, & nf90_get_var, & nf90_close !------------------------------------------------------------------------------ ! ... Dummy arguments !------------------------------------------------------------------------------ integer, intent(inout) :: lng_indexer(:) character(len=*) :: xs_long_file !------------------------------------------------------------------------------ ! ... Local variables !------------------------------------------------------------------------------ integer :: varid, dimid, ndx integer :: ncid integer :: iret integer :: i, k, m, n integer :: wrk_ndx(phtcnt) character(len=256) :: locfn Masterproc_only : if( masterproc ) then !------------------------------------------------------------------------------ ! ... open NetCDF File !------------------------------------------------------------------------------ call getfil(xs_long_file, locfn, 0) iret = nf90_open(trim(locfn), NF90_NOWRITE, ncid) !------------------------------------------------------------------------------ ! ... get dimensions !------------------------------------------------------------------------------ iret = nf90_inq_dimid( ncid, 'numtemp', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= nt ) iret = nf90_inq_dimid( ncid, 'numwl', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= nw ) iret = nf90_inq_dimid( ncid, 'numprs', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= np_xs ) !------------------------------------------------------------------------------ ! ... check for cross section in dataset !------------------------------------------------------------------------------ do m = 1,phtcnt if( pht_alias_lst(m,2) == ' ' ) then iret = nf90_inq_varid( ncid, rxt_tag_lst(m), varid ) if( iret == nf90_noerr ) then lng_indexer(m) = varid end if else if( pht_alias_lst(m,2) == 'userdefined' ) then lng_indexer(m) = -1 else iret = nf90_inq_varid( ncid, pht_alias_lst(m,2), varid ) if( iret == nf90_noerr ) then lng_indexer(m) = varid else write(iulog,*) 'get_xsqy : ',rxt_tag_lst(m)(:len_trim(rxt_tag_lst(m))),' alias ', & pht_alias_lst(m,2)(:len_trim(pht_alias_lst(m,2))),' not in dataset' call endrun end if end if end do numj = 0 do m = 1,phtcnt if( lng_indexer(m) > 0 ) then if( any( lng_indexer(:m-1) == lng_indexer(m) ) ) then cycle end if numj = numj + 1 end if end do !------------------------------------------------------------------------------ ! ... allocate arrays !------------------------------------------------------------------------------ allocate( xsqy(numj,nw,nt,np_xs),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_xsqy', 'xsqy', numj*nt*np_xs*nw ) end if allocate( prs(np_xs),dprs(np_xs-1),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_xsqy', 'prs,dprs', np_xs ) end if allocate( xs_o2b(nw,nt,np_xs),xs_o3a(nw,nt,np_xs),xs_o3b(nw,nt,np_xs),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_xsqy', 'xs_o2b ... xs_o3b', np_xs ) end if !------------------------------------------------------------------------------ ! ... read cross sections !------------------------------------------------------------------------------ ndx = 0 do m = 1,phtcnt if( lng_indexer(m) > 0 ) then if( any( lng_indexer(:m-1) == lng_indexer(m) ) ) then cycle end if ndx = ndx + 1 iret = nf90_get_var( ncid, lng_indexer(m), xsqy(ndx,:,:,:) ) end if end do if( ndx /= numj ) then write(iulog,*) 'get_xsqy : ndx count /= cross section count' call endrun end if iret = nf90_inq_varid( ncid, 'jo2_b', varid ) iret = nf90_get_var( ncid, varid, xs_o2b ) iret = nf90_inq_varid( ncid, 'jo3_a', varid ) iret = nf90_get_var( ncid, varid, xs_o3a ) iret = nf90_inq_varid( ncid, 'jo3_b', varid ) iret = nf90_get_var( ncid, varid, xs_o3b ) !------------------------------------------------------------------------------ ! ... setup final lng_indexer !------------------------------------------------------------------------------ ndx = 0 wrk_ndx(:) = lng_indexer(:) do m = 1,phtcnt if( wrk_ndx(m) > 0 ) then ndx = ndx + 1 i = wrk_ndx(m) where( wrk_ndx(:) == i ) lng_indexer(:) = ndx wrk_ndx(:) = -100000 end where end if end do iret = nf90_inq_varid( ncid, 'pressure', varid ) iret = nf90_get_var( ncid, varid, prs ) iret = nf90_close( ncid ) end if Masterproc_only #ifdef SPMD ! call mpibarrier( mpicom ) call mpibcast( numj, 1, mpiint, 0, mpicom ) call mpibcast( nt, 1, mpiint, 0, mpicom ) call mpibcast( nw, 1, mpiint, 0, mpicom ) call mpibcast( np_xs, 1, mpiint, 0, mpicom ) call mpibcast( lng_indexer, phtcnt, mpiint, 0, mpicom ) #endif if( .not. masterproc ) then !------------------------------------------------------------------------------ ! ... allocate arrays !------------------------------------------------------------------------------ allocate( xsqy(numj,nw,nt,np_xs),stat=iret ) if( iret /= nf90_noerr) then write(iulog,*) 'get_xsqy : failed to allocate xsqy ; error = ',iret call endrun end if allocate( prs(np_xs),dprs(np_xs-1),stat=iret ) if( iret /= nf90_noerr) then write(iulog,*) 'get_xsqy : failed to allocate prs,dprs ; error = ',iret call endrun end if allocate( xs_o2b(nw,nt,np_xs),xs_o3a(nw,nt,np_xs),xs_o3b(nw,nt,np_xs),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_xsqy', 'xs_o2b ... xs_o3b', np_xs ) end if end if #ifdef SPMD ! call mpibarrier( mpicom ) call mpibcast( prs, np_xs, mpir8, 0, mpicom ) call mpibcast( xsqy, numj*nt*np_xs*nw, mpir4, 0, mpicom ) call mpibcast( xs_o2b, nw*nt*np_xs, mpir8, 0, mpicom ) call mpibcast( xs_o3a, nw*nt*np_xs, mpir8, 0, mpicom ) call mpibcast( xs_o3b, nw*nt*np_xs, mpir8, 0, mpicom ) #endif dprs(:np_xs-1) = 1._r8/(prs(1:np_xs-1) - prs(2:np_xs)) end subroutine get_xsqy subroutine get_rsf(rsf_file) !=============================================================================! ! PURPOSE: ! ! Reads a NetCDF file that contains: ! Radiative Souce function ! !=============================================================================! ! PARAMETERS: ! ! Input: ! ! filepath.... NetCDF file that contains the RSF ! ! ! ! Output: ! ! rsf ........ Radiative Source Function (quanta cm-2 sec-1 ! ! ! ! EDIT HISTORY: ! ! Created by Doug Kinnison, 3/14/2002 ! !=============================================================================! use ioFileMod, only : getfil use error_messages, only : alloc_err use netcdf, only: & nf90_open, & nf90_nowrite, & nf90_inq_dimid, & nf90_inquire_dimension, & nf90_inq_varid, & nf90_noerr, & nf90_get_var, & nf90_close !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ character(len=*) :: rsf_file !------------------------------------------------------------------------------ ! ... local variables !------------------------------------------------------------------------------ integer :: varid, dimid integer :: ncid integer :: i, j, k, l, w integer :: iret integer :: count(5) integer :: start(5) real(r8) :: wrk character(len=256) :: locfn Masterproc_only : if( masterproc ) then !------------------------------------------------------------------------------ ! ... open NetCDF File !------------------------------------------------------------------------------ call getfil(rsf_file, locfn, 0) iret = nf90_open(trim(locfn), NF90_NOWRITE, ncid) !------------------------------------------------------------------------------ ! ... get dimensions !------------------------------------------------------------------------------ iret = nf90_inq_dimid( ncid, 'numz', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= nump ) iret = nf90_inq_dimid( ncid, 'numsza', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= numsza ) iret = nf90_inq_dimid( ncid, 'numalb', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= numalb ) iret = nf90_inq_dimid( ncid, 'numcolo3fact', dimid ) iret = nf90_inquire_dimension( ncid, dimid,len= numcolo3 ) end if Masterproc_only #ifdef SPMD ! call mpibarrier( mpicom ) call mpibcast( nump, 1, mpiint, 0, mpicom ) call mpibcast( numsza, 1, mpiint, 0, mpicom ) call mpibcast( numalb, 1, mpiint, 0, mpicom ) call mpibcast( numcolo3, 1, mpiint, 0, mpicom ) #endif !------------------------------------------------------------------------------ ! ... allocate arrays !------------------------------------------------------------------------------ allocate( wc(nw),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'wc', nw ) end if allocate( wlintv(nw),we(nw+1),etfphot(nw),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'wlintv,etfphot', nw ) end if allocate( bde_o2_b(nw),bde_o3_a(nw),bde_o3_b(nw),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'bde', nw ) end if allocate( p(nump),del_p(nump-1),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'p,delp', nump ) end if allocate( sza(numsza),del_sza(numsza-1),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'sza,del_sza', numsza ) end if allocate( alb(numalb),del_alb(numalb-1),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'alb,del_alb', numalb ) end if allocate( o3rat(numcolo3),del_o3rat(numcolo3-1),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'o3rat,del_o3rat', numcolo3 ) end if allocate( colo3(nump),stat=iret ) if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'colo3', nump ) end if allocate( rsf_tab(nw,nump,numsza,numcolo3,numalb),stat=iret ) if( iret /= 0 ) then write(iulog,*) 'get_rsf : dimensions = ',nw,nump,numsza,numcolo3,numalb call alloc_err( iret, 'get_rsf', 'rsf_tab', numalb*numcolo3*numsza*nump ) end if Masterproc_only2 : if( masterproc ) then !------------------------------------------------------------------------------ ! ... read variables !------------------------------------------------------------------------------ iret = nf90_inq_varid( ncid, 'wc', varid ) iret = nf90_get_var( ncid, varid, wc ) iret = nf90_inq_varid( ncid, 'wlintv', varid ) iret = nf90_get_var( ncid, varid, wlintv ) iret = nf90_inq_varid( ncid, 'pm', varid ) iret = nf90_get_var( ncid, varid, p ) iret = nf90_inq_varid( ncid, 'sza', varid ) iret = nf90_get_var( ncid, varid, sza ) iret = nf90_inq_varid( ncid, 'alb', varid ) iret = nf90_get_var( ncid, varid, alb ) iret = nf90_inq_varid( ncid, 'colo3fact', varid ) iret = nf90_get_var( ncid, varid, o3rat ) iret = nf90_inq_varid( ncid, 'colo3', varid ) iret = nf90_get_var( ncid, varid, colo3 ) iret = nf90_inq_varid( ncid, 'RSF', varid ) if (masterproc) then write(iulog,*) ' ' write(iulog,*) '----------------------------------------------' write(iulog,*) 'get_rsf: numalb, numcolo3, numsza, nump = ', & numalb, numcolo3, numsza, nump write(iulog,*) 'get_rsf: size of rsf_tab = ',size(rsf_tab,dim=1),size(rsf_tab,dim=2), & size(rsf_tab,dim=3),size(rsf_tab,dim=4),size(rsf_tab,dim=5) write(iulog,*) '----------------------------------------------' write(iulog,*) ' ' endif iret = nf90_get_var( ncid, varid, rsf_tab ) iret = nf90_close( ncid ) do w = 1,nw wrk = wlintv(w) rsf_tab(w,:,:,:,:) = wrk*rsf_tab(w,:,:,:,:) enddo end if Masterproc_only2 #ifdef SPMD call mpibcast( wc, nw, mpir8, 0, mpicom ) call mpibcast( wlintv, nw, mpir8, 0, mpicom ) call mpibcast( p, nump, mpir8, 0, mpicom ) call mpibcast( sza, numsza, mpir8, 0, mpicom ) call mpibcast( alb, numalb, mpir8, 0, mpicom ) call mpibcast( o3rat, numcolo3, mpir8, 0, mpicom ) call mpibcast( colo3, nump, mpir8, 0, mpicom ) call mpibcast( rsf_tab, nw*numalb*numcolo3*numsza*nump, mpir4, 0, mpicom ) #endif #ifdef USE_BDE if (masterproc) write(iulog,*) 'Jlong using bdes' bde_o2_b(:) = max( 0._r8, hc*(wc_o2_b - wc(:))/(wc_o2_b*wc(:)) ) bde_o3_a(:) = max( 0._r8, hc*(wc_o3_a - wc(:))/(wc_o3_a*wc(:)) ) bde_o3_b(:) = max( 0._r8, hc*(wc_o3_b - wc(:))/(wc_o3_b*wc(:)) ) #else if (masterproc) write(iulog,*) 'Jlong not using bdes' bde_o2_b(:) = hc/wc(:) bde_o3_a(:) = hc/wc(:) bde_o3_b(:) = hc/wc(:) #endif del_p(:nump-1) = 1._r8/abs(p(1:nump-1)- p(2:nump)) del_sza(:numsza-1) = 1._r8/(sza(2:numsza) - sza(1:numsza-1)) del_alb(:numalb-1) = 1._r8/(alb(2:numalb) - alb(1:numalb-1)) del_o3rat(:numcolo3-1) = 1._r8/(o3rat(2:numcolo3) - o3rat(1:numcolo3-1)) end subroutine get_rsf subroutine jlong_timestep_init !--------------------------------------------------------------- ! ... set etfphot if required !--------------------------------------------------------------- use time_manager, only : is_end_curr_day use mo_util, only : rebin use solar_data, only : data_nw => nbins, data_we => we, data_etf => sol_etf implicit none if (.not.jlong_used) return call rebin( data_nw, nw, data_we, we, data_etf, etfphot ) end subroutine jlong_timestep_init subroutine jlong_hrates( nlev, sza_in, alb_in, p_in, t_in, & mw, o2_vmr, o3_vmr, colo3_in, qrl_col, & cparg, kbot ) !============================================================================== ! Purpose: ! To calculate the thermal heating rates longward of 200nm. !============================================================================== ! Approach: ! 1) Reads the Cross Section*QY NetCDF file ! 2) Given a temperature profile, derives the appropriate XS*QY ! ! 3) Reads the Radiative Source function (RSF) NetCDF file ! Units = quanta cm-2 sec-1 ! ! 4) Indices are supplied to select a RSF that is consistent with ! the reference atmosphere in TUV (for direct comparision of J's). ! This approach will be replaced in the global model. Here colo3, zenith ! angle, and altitude will be inputed and the correct entry in the table ! will be derived. !============================================================================== use physconst, only : avogad use error_messages, only : alloc_err implicit none !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ integer, intent (in) :: nlev ! number vertical levels integer, intent (in) :: kbot ! heating levels real(r8), intent(in) :: o2_vmr(nlev) ! o2 conc (mol/mol) real(r8), intent(in) :: o3_vmr(nlev) ! o3 conc (mol/mol) real(r8), intent(in) :: sza_in ! solar zenith angle (degrees) real(r8), intent(in) :: alb_in(nlev) ! albedo real(r8), intent(in) :: p_in(nlev) ! midpoint pressure (hPa) real(r8), intent(in) :: t_in(nlev) ! Temperature profile (K) real(r8), intent(in) :: colo3_in(nlev) ! o3 column density (molecules/cm^3) real(r8), intent(in) :: mw(nlev) ! atms molecular weight real(r8), intent(in) :: cparg(nlev) ! specific heat capacity real(r8), intent(inout) :: qrl_col(:,:) ! heating rates !---------------------------------------------------------------------- ! ... local variables !---------------------------------------------------------------------- integer :: astat integer :: k, km, m integer :: t_index ! Temperature index integer :: pndx real(r8) :: ptarget real(r8) :: delp real(r8) :: hfactor real(r8), allocatable :: rsf(:,:) ! Radiative source function real(r8), allocatable :: xswk(:,:) ! working xsection array real(r8), allocatable :: wrk(:) ! work vector !---------------------------------------------------------------------- ! ... allocate variables !---------------------------------------------------------------------- allocate( rsf(nw,kbot),stat=astat ) if( astat /= 0 ) then call alloc_err( astat, 'jlong_hrates', 'rsf', nw*nlev ) end if allocate( xswk(nw,3),wrk(nw),stat=astat ) if( astat /= 0 ) then call alloc_err( astat, 'jlong_hrates', 'xswk,wrk', 3*nw ) end if !---------------------------------------------------------------------- ! ... interpolate table rsf to model variables !---------------------------------------------------------------------- call interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf ) !------------------------------------------------------------------------------ ! ... calculate thermal heating rates for wavelengths >200nm !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! LLNL LUT approach to finding temperature index... ! Calculate the temperature index into the cross section ! data which lists coss sections for temperatures from ! 150 to 350 degrees K. Make sure the index is a value ! between 1 and 201. !------------------------------------------------------------------------------ qrl_col(kbot+1:nlev,:) = 0._r8 level_loop_1 : & do k = 1,kbot !---------------------------------------------------------------------- ! ... get index into xsqy !---------------------------------------------------------------------- t_index = t_in(k) - 148.5_r8 t_index = min( 201,max( t_index,1) ) !---------------------------------------------------------------------- ! ... find pressure level !---------------------------------------------------------------------- ptarget = p_in(k) if( ptarget >= prs(1) ) then xswk(:,1) = xs_o2b(:,t_index,1) xswk(:,2) = xs_o3a(:,t_index,1) xswk(:,3) = xs_o3b(:,t_index,1) else if( ptarget <= prs(np_xs) ) then xswk(:,1) = xs_o2b(:,t_index,np_xs) xswk(:,2) = xs_o3a(:,t_index,np_xs) xswk(:,3) = xs_o3b(:,t_index,np_xs) else do km = 2,np_xs if( ptarget >= prs(km) ) then pndx = km - 1 delp = (prs(pndx) - ptarget)*dprs(pndx) exit end if end do xswk(:,1) = xs_o2b(:,t_index,pndx) & + delp*(xs_o2b(:,t_index,pndx+1) - xs_o2b(:,t_index,pndx)) xswk(:,2) = xs_o3a(:,t_index,pndx) & + delp*(xs_o3a(:,t_index,pndx+1) - xs_o3a(:,t_index,pndx)) xswk(:,3) = xs_o3b(:,t_index,pndx) & + delp*(xs_o3b(:,t_index,pndx+1) - xs_o3b(:,t_index,pndx)) end if hfactor = avogad/(cparg(k)*mw(k)) wrk(:) = xswk(:,1)*bde_o2_b(:) qrl_col(k,2) = dot_product( wrk(:),rsf(:,k) ) * o2_vmr(k) * hfactor wrk(:) = xswk(:,2)*bde_o3_a(:) qrl_col(k,3) = dot_product( wrk(:),rsf(:,k) ) * o3_vmr(k) * hfactor wrk(:) = xswk(:,3)*bde_o3_b(:) qrl_col(k,4) = dot_product( wrk(:),rsf(:,k) ) * o3_vmr(k) * hfactor end do level_loop_1 deallocate( rsf, xswk, wrk ) end subroutine jlong_hrates subroutine jlong_photo( nlev, sza_in, alb_in, p_in, t_in, & colo3_in, j_long ) !============================================================================== ! Purpose: ! To calculate the total J for selective species longward of 200nm. !============================================================================== ! Approach: ! 1) Reads the Cross Section*QY NetCDF file ! 2) Given a temperature profile, derives the appropriate XS*QY ! ! 3) Reads the Radiative Source function (RSF) NetCDF file ! Units = quanta cm-2 sec-1 ! ! 4) Indices are supplied to select a RSF that is consistent with ! the reference atmosphere in TUV (for direct comparision of J's). ! This approach will be replaced in the global model. Here colo3, zenith ! angle, and altitude will be inputed and the correct entry in the table ! will be derived. !============================================================================== use spmd_utils, only : masterproc use error_messages, only : alloc_err implicit none !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ integer, intent (in) :: nlev ! number vertical levels real(r8), intent(in) :: sza_in ! solar zenith angle (degrees) real(r8), intent(in) :: alb_in(nlev) ! albedo real(r8), intent(in) :: p_in(nlev) ! midpoint pressure (hPa) real(r8), intent(in) :: t_in(nlev) ! Temperature profile (K) real(r8), intent(in) :: colo3_in(nlev) ! o3 column density (molecules/cm^3) real(r8), intent(out) :: j_long(:,:) ! photo rates (1/s) !---------------------------------------------------------------------- ! ... local variables !---------------------------------------------------------------------- integer :: astat integer :: k, km, m integer :: wn integer :: t_index ! Temperature index integer :: pndx real(r8) :: ptarget real(r8) :: delp real(r8) :: hfactor real(r8), allocatable :: rsf(:,:) ! Radiative source function real(r8), allocatable :: xswk(:,:) ! working xsection array !---------------------------------------------------------------------- ! ... allocate variables !---------------------------------------------------------------------- allocate( rsf(nw,nlev),stat=astat ) if( astat /= 0 ) then call alloc_err( astat, 'jlong_photo', 'rsf', nw*nlev ) end if allocate( xswk(numj,nw),stat=astat ) if( astat /= 0 ) then call alloc_err( astat, 'jlong_photo', 'xswk', numj*nw ) end if !---------------------------------------------------------------------- ! ... interpolate table rsf to model variables !---------------------------------------------------------------------- call interpolate_rsf( alb_in, sza_in, p_in, colo3_in, nlev, rsf ) !------------------------------------------------------------------------------ ! ... calculate total Jlong for wavelengths >200nm !------------------------------------------------------------------------------ !------------------------------------------------------------------------------ ! LLNL LUT approach to finding temperature index... ! Calculate the temperature index into the cross section ! data which lists coss sections for temperatures from ! 150 to 350 degrees K. Make sure the index is a value ! between 1 and 201. !------------------------------------------------------------------------------ level_loop_1 : & do k = 1,nlev !---------------------------------------------------------------------- ! ... get index into xsqy !---------------------------------------------------------------------- t_index = t_in(k) - 148.5_r8 t_index = min( 201,max( t_index,1) ) !---------------------------------------------------------------------- ! ... find pressure level !---------------------------------------------------------------------- ptarget = p_in(k) if( ptarget >= prs(1) ) then do wn = 1,nw xswk(:,wn) = xsqy(:,wn,t_index,1) end do else if( ptarget <= prs(np_xs) ) then do wn = 1,nw xswk(:,wn) = xsqy(:,wn,t_index,np_xs) end do else do km = 2,np_xs if( ptarget >= prs(km) ) then pndx = km - 1 delp = (prs(pndx) - ptarget)*dprs(pndx) exit end if end do do wn = 1,nw xswk(:,wn) = xsqy(:,wn,t_index,pndx) & + delp*(xsqy(:,wn,t_index,pndx+1) - xsqy(:,wn,t_index,pndx)) end do end if #ifdef USE_ESSL call dgemm( 'N', 'N', numj, 1, nw, & 1._r8, xswk, numj, rsf(1,k), nw, & 0._r8, j_long(1,k), numj ) #else j_long(:,k) = matmul( xswk(:,:),rsf(:,k) ) #endif end do level_loop_1 deallocate( rsf, xswk ) end subroutine jlong_photo !---------------------------------------------------------------------- !---------------------------------------------------------------------- ! ... interpolate table rsf to model variables !---------------------------------------------------------------------- !---------------------------------------------------------------------- subroutine interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf ) use error_messages, only : alloc_err implicit none !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ real(r8), intent(in) :: alb_in(:) ! albedo real(r8), intent(in) :: sza_in ! solar zenith angle (degrees) integer, intent(in) :: kbot ! heating levels real(r8), intent(in) :: p_in(:) ! midpoint pressure (hPa) real(r8), intent(in) :: colo3_in(:) ! o3 column density (molecules/cm^3) real(r8), intent(out) :: rsf(:,:) !---------------------------------------------------------------------- ! ... local variables !---------------------------------------------------------------------- integer :: astat integer :: is, iv, ial integer :: isp1, ivp1, ialp1 real(r8), dimension(3) :: dels real(r8), dimension(0:1,0:1,0:1) :: wghtl, wghtu real(r8) :: psum_u real(r8), allocatable :: psum_l(:) real(r8) :: v3ratl, v3ratu integer :: pind, albind real(r8) :: wrk0, wrk1, wght1 integer :: iz, k, m integer :: izl integer :: ratindl, ratindu integer :: wn !---------------------------------------------------------------------- ! ... allocate variables !---------------------------------------------------------------------- allocate( psum_l(nw),stat=astat ) if( astat /= 0 ) then call alloc_err( astat, 'jlong_hrates', 'psum_l', nw ) end if !---------------------------------------------------------------------- ! ... find the zenith angle index ( same for all levels ) !---------------------------------------------------------------------- do is = 1,numsza if( sza(is) > sza_in ) then exit end if end do is = max( min( is,numsza ) - 1,1 ) isp1 = is + 1 dels(1) = max( 0._r8,min( 1._r8,(sza_in - sza(is)) * del_sza(is) ) ) wrk0 = 1._r8 - dels(1) izl = 2 Level_loop : & do k = kbot,1,-1 !---------------------------------------------------------------------- ! ... find albedo indicies !---------------------------------------------------------------------- do ial = 1,numalb if( alb(ial) > alb_in(k) ) then exit end if end do albind = max( min( ial,numalb ) - 1,1 ) !---------------------------------------------------------------------- ! ... find pressure level indicies !---------------------------------------------------------------------- if( p_in(k) > p(1) ) then pind = 2 wght1 = 1._r8 else if( p_in(k) <= p(nump) ) then pind = nump wght1 = 0._r8 else do iz = izl,nump if( p(iz) < p_in(k) ) then izl = iz exit end if end do pind = max( min( iz,nump ),2 ) wght1 = max( 0._r8,min( 1._r8,(p_in(k) - p(pind)) * del_p(pind-1) ) ) end if !---------------------------------------------------------------------- ! ... find "o3 ratios" !---------------------------------------------------------------------- v3ratu = colo3_in(k) / colo3(pind-1) do iv = 1,numcolo3 if( o3rat(iv) > v3ratu ) then exit end if end do ratindu = max( min( iv,numcolo3 ) - 1,1 ) if( colo3(pind) /= 0._r8 ) then v3ratl = colo3_in(k) / colo3(pind) do iv = 1,numcolo3 if( o3rat(iv) > v3ratl ) then exit end if end do ratindl = max( min( iv,numcolo3 ) - 1,1 ) else ratindl = ratindu v3ratl = o3rat(ratindu) end if !---------------------------------------------------------------------- ! ... compute the weigths !---------------------------------------------------------------------- ial = albind ialp1 = ial + 1 iv = ratindl dels(2) = max( 0._r8,min( 1._r8,(v3ratl - o3rat(iv)) * del_o3rat(iv) ) ) dels(3) = max( 0._r8,min( 1._r8,(alb_in(k) - alb(ial)) * del_alb(ial) ) ) wrk1 = (1._r8 - dels(2))*(1._r8 - dels(3)) wghtl(0,0,0) = wrk0*wrk1 wghtl(1,0,0) = dels(1)*wrk1 wrk1 = (1._r8 - dels(2))*dels(3) wghtl(0,0,1) = wrk0*wrk1 wghtl(1,0,1) = dels(1)*wrk1 wrk1 = dels(2)*(1._r8 - dels(3)) wghtl(0,1,0) = wrk0*wrk1 wghtl(1,1,0) = dels(1)*wrk1 wrk1 = dels(2)*dels(3) wghtl(0,1,1) = wrk0*wrk1 wghtl(1,1,1) = dels(1)*wrk1 iv = ratindu dels(2) = max( 0._r8,min( 1._r8,(v3ratu - o3rat(iv)) * del_o3rat(iv) ) ) wrk1 = (1._r8 - dels(2))*(1._r8 - dels(3)) wghtu(0,0,0) = wrk0*wrk1 wghtu(1,0,0) = dels(1)*wrk1 wrk1 = (1._r8 - dels(2))*dels(3) wghtu(0,0,1) = wrk0*wrk1 wghtu(1,0,1) = dels(1)*wrk1 wrk1 = dels(2)*(1._r8 - dels(3)) wghtu(0,1,0) = wrk0*wrk1 wghtu(1,1,0) = dels(1)*wrk1 wrk1 = dels(2)*dels(3) wghtu(0,1,1) = wrk0*wrk1 wghtu(1,1,1) = dels(1)*wrk1 iz = pind iv = ratindl ivp1 = iv + 1 do wn = 1,nw psum_l(wn) = wghtl(0,0,0) * rsf_tab(wn,iz,is,iv,ial) & + wghtl(0,0,1) * rsf_tab(wn,iz,is,iv,ialp1) & + wghtl(0,1,0) * rsf_tab(wn,iz,is,ivp1,ial) & + wghtl(0,1,1) * rsf_tab(wn,iz,is,ivp1,ialp1) & + wghtl(1,0,0) * rsf_tab(wn,iz,isp1,iv,ial) & + wghtl(1,0,1) * rsf_tab(wn,iz,isp1,iv,ialp1) & + wghtl(1,1,0) * rsf_tab(wn,iz,isp1,ivp1,ial) & + wghtl(1,1,1) * rsf_tab(wn,iz,isp1,ivp1,ialp1) end do iz = iz - 1 iv = ratindu ivp1 = iv + 1 do wn = 1,nw psum_u = wghtu(0,0,0) * rsf_tab(wn,iz,is,iv,ial) & + wghtu(0,0,1) * rsf_tab(wn,iz,is,iv,ialp1) & + wghtu(0,1,0) * rsf_tab(wn,iz,is,ivp1,ial) & + wghtu(0,1,1) * rsf_tab(wn,iz,is,ivp1,ialp1) & + wghtu(1,0,0) * rsf_tab(wn,iz,isp1,iv,ial) & + wghtu(1,0,1) * rsf_tab(wn,iz,isp1,iv,ialp1) & + wghtu(1,1,0) * rsf_tab(wn,iz,isp1,ivp1,ial) & + wghtu(1,1,1) * rsf_tab(wn,iz,isp1,ivp1,ialp1) rsf(wn,k) = (psum_l(wn) + wght1*(psum_u - psum_l(wn))) end do !------------------------------------------------------------------------------ ! etfphot comes in as photons/cm^2/sec/nm (rsf includes the wlintv factor -- nm) ! ... --> convert to photons/cm^2/s !------------------------------------------------------------------------------ rsf(:,k) = etfphot(:) * rsf(:,k) end do Level_loop deallocate( psum_l ) end subroutine interpolate_rsf end module mo_jlong