diff --git a/bld/namelist_files/namelist_definition.xml b/bld/namelist_files/namelist_definition.xml index 37a569022a..a80f3560b4 100644 --- a/bld/namelist_files/namelist_definition.xml +++ b/bld/namelist_files/namelist_definition.xml @@ -4609,6 +4609,21 @@ Use spectral scaling in the radiation heating Default: set by build-namelist + +Solar shade insolation factor + + + +Solar shade insolation factor + + + +Solar shade insolation factor + + nbins, data_we => we, data_etf => sol_etf implicit none @@ -100,20 +101,8 @@ subroutine jlong_init( xs_long_file, rsf_file, lng_indexer ) 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 ) @@ -182,18 +171,18 @@ subroutine get_xsqy( xs_long_file, lng_indexer ) 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 + 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 + 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' + pht_alias_lst(m,2)(:len_trim(pht_alias_lst(m,2))),' not in dataset' call endrun end if end if @@ -213,15 +202,15 @@ subroutine get_xsqy( xs_long_file, lng_indexer ) !------------------------------------------------------------------------------ allocate( xsqy(numj,nw,nt,np_xs),stat=iret ) - if( iret /= 0 ) then + 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 + 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 + if( iret /= 0 ) then call alloc_err( iret, 'get_xsqy', 'xs_o2b ... xs_o3b', np_xs ) end if !------------------------------------------------------------------------------ @@ -281,17 +270,17 @@ subroutine get_xsqy( xs_long_file, lng_indexer ) ! ... allocate arrays !------------------------------------------------------------------------------ allocate( xsqy(numj,nw,nt,np_xs),stat=iret ) - if( iret /= nf90_noerr) then + 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 + 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 + if( iret /= 0 ) then call alloc_err( iret, 'get_xsqy', 'xs_o2b ... xs_o3b', np_xs ) end if end if @@ -384,39 +373,41 @@ subroutine get_rsf(rsf_file) ! ... allocate arrays !------------------------------------------------------------------------------ allocate( wc(nw),stat=iret ) - if( iret /= 0 ) then + 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 + allocate( wlintv(nw),we(nw+1),etfphot(nw,pcols,begchunk:endchunk),stat=iret ) + if( iret /= 0 ) then call alloc_err( iret, 'get_rsf', 'wlintv,etfphot', nw ) end if + etfphot = nan + allocate( bde_o2_b(nw),bde_o3_a(nw),bde_o3_b(nw),stat=iret ) - if( iret /= 0 ) then + 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 + 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 + 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 + 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 + 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 + 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 + 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 @@ -441,7 +432,7 @@ subroutine get_rsf(rsf_file) iret = nf90_get_var( ncid, varid, colo3 ) iret = nf90_inq_varid( ncid, 'RSF', varid ) - + if (masterproc) then write(iulog,*) ' ' write(iulog,*) '----------------------------------------------' @@ -501,21 +492,32 @@ subroutine jlong_timestep_init use mo_util, only : rebin use solar_irrad_data,only : data_nw => nbins, data_we => we, data_etf => sol_etf + use solar_shade, only: sun_shade + use phys_grid, only : get_ncols_p implicit none + integer :: i, c, ncols + if (.not.jlong_used) return - call rebin( data_nw, nw, data_we, we, data_etf, etfphot ) + do c = begchunk,endchunk + ncols = get_ncols_p(c) + do i = 1,ncols + ! apply sun shade factor to input ETF + call rebin( data_nw, nw, data_we, we, data_etf(:)*sun_shade(:,i,c), etfphot(:,i,c) ) + end do + end do + end subroutine jlong_timestep_init - subroutine jlong_hrates( nlev, sza_in, alb_in, p_in, t_in, & + subroutine jlong_hrates( icol, lchnk, 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. +! Purpose: +! To calculate the thermal heating rates longward of 200nm. !============================================================================== ! Approach: ! 1) Reads the Cross Section*QY NetCDF file @@ -539,6 +541,7 @@ subroutine jlong_hrates( nlev, sza_in, alb_in, p_in, t_in, & !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ + integer, intent (in) :: icol, lchnk integer, intent (in) :: nlev ! number vertical levels integer, intent (in) :: kbot ! heating levels real(r8), intent(in) :: o2_vmr(nlev) ! o2 conc (mol/mol) @@ -581,7 +584,7 @@ subroutine jlong_hrates( nlev, sza_in, alb_in, p_in, t_in, & !---------------------------------------------------------------------- ! ... interpolate table rsf to model variables !---------------------------------------------------------------------- - call interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf ) + call interpolate_rsf( icol, lchnk, alb_in, sza_in, p_in, colo3_in, kbot, rsf ) !------------------------------------------------------------------------------ ! ... calculate thermal heating rates for wavelengths >200nm @@ -641,11 +644,11 @@ subroutine jlong_hrates( nlev, sza_in, alb_in, p_in, t_in, & end subroutine jlong_hrates - subroutine jlong_photo( nlev, sza_in, alb_in, p_in, t_in, & + subroutine jlong_photo( icol, lchnk, 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. +! Purpose: +! To calculate the total J for selective species longward of 200nm. !============================================================================== ! Approach: ! 1) Reads the Cross Section*QY NetCDF file @@ -669,6 +672,7 @@ subroutine jlong_photo( nlev, sza_in, alb_in, p_in, t_in, & !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ + integer, intent (in) :: icol, lchnk 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 @@ -706,7 +710,7 @@ subroutine jlong_photo( nlev, sza_in, alb_in, p_in, t_in, & !---------------------------------------------------------------------- ! ... interpolate table rsf to model variables !---------------------------------------------------------------------- - call interpolate_rsf( alb_in, sza_in, p_in, colo3_in, nlev, rsf ) + call interpolate_rsf( icol, lchnk, alb_in, sza_in, p_in, colo3_in, nlev, rsf ) !------------------------------------------------------------------------------ ! ... calculate total Jlong for wavelengths >200nm @@ -768,7 +772,7 @@ end subroutine jlong_photo ! ... interpolate table rsf to model variables !---------------------------------------------------------------------- !---------------------------------------------------------------------- - subroutine interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf ) + subroutine interpolate_rsf( icol, lchnk, alb_in, sza_in, p_in, colo3_in, kbot, rsf ) use error_messages, only : alloc_err @@ -777,6 +781,7 @@ subroutine interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf ) !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ + integer , intent(in) :: icol, lchnk real(r8), intent(in) :: alb_in(:) ! albedo real(r8), intent(in) :: sza_in ! solar zenith angle (degrees) integer, intent(in) :: kbot ! heating levels @@ -947,9 +952,9 @@ subroutine interpolate_rsf( alb_in, sza_in, p_in, colo3_in, kbot, rsf ) end do !------------------------------------------------------------------------------ ! etfphot comes in as photons/cm^2/sec/nm (rsf includes the wlintv factor -- nm) -! ... --> convert to photons/cm^2/s +! ... --> convert to photons/cm^2/s !------------------------------------------------------------------------------ - rsf(:,k) = etfphot(:) * rsf(:,k) + rsf(:,k) = etfphot(:,icol,lchnk) * rsf(:,k) end do Level_loop diff --git a/src/chemistry/mozart/mo_jshort.F90 b/src/chemistry/mozart/mo_jshort.F90 index e6b635c7be..6611d46b25 100644 --- a/src/chemistry/mozart/mo_jshort.F90 +++ b/src/chemistry/mozart/mo_jshort.F90 @@ -13,6 +13,8 @@ module mo_jshort use spmd_utils, only : masterproc use ppgrid, only : pver use phys_control, only : waccmx_is + use infnan, only : nan, assignment(=) + use ppgrid, only : pcols, begchunk, endchunk implicit none @@ -64,8 +66,8 @@ module mo_jshort real(r8), allocatable :: bde_o2_b(:) real(r8), allocatable :: bde_o3_a(:) real(r8), allocatable :: bde_o3_b(:) - real(r8), allocatable :: etfphot(:) - real(r8), allocatable :: etfphot_ms93(:) + real(r8), allocatable :: etfphot(:,:,:) + real(r8), allocatable :: etfphot_ms93(:,:,:) real(r8), allocatable :: xs_o2src(:) real(r8), allocatable :: xs_o3a(:) real(r8), allocatable :: xs_o3b(:) @@ -76,7 +78,7 @@ module mo_jshort subroutine jshort_init( xs_coef_file, xs_short_file, sht_indexer ) !------------------------------------------------------------------------------ ! ... initialize photorates for 120nm <= lambda <= 200nm -!------------------------------------------------------------------------------ +!------------------------------------------------------------------------------ use mo_util, only : rebin use solar_irrad_data, only : data_nbins=>nbins, data_we => we, data_etf => sol_etf @@ -85,7 +87,7 @@ subroutine jshort_init( xs_coef_file, xs_short_file, sht_indexer ) !------------------------------------------------------------------------------ ! ... dummy arguments -!------------------------------------------------------------------------------ +!------------------------------------------------------------------------------ character(len=*), intent(in) :: xs_coef_file, xs_short_file integer, intent(inout) :: sht_indexer(:) @@ -117,18 +119,6 @@ subroutine jshort_init( xs_coef_file, xs_short_file, sht_indexer ) write(iulog,*) 'jshort_init: we range' write(iulog,'(1p,5g15.7)') minval(we(:)),maxval(we(:)) end if - call rebin( data_nbins, nw, data_we, we, data_etf, etfphot ) - if(masterproc) then - write(iulog,*) 'jshort_init: etfphot' - write(iulog,'(1p,5g15.7)') etfphot(:) - write(iulog,*) '-------------------------------------------' - write(iulog,*) ' ' - write(iulog,*) 'jshort_init: diagnostics for ms93' - call rebin( data_nbins, nw_ms93, data_we, we_ms, data_etf, etfphot_ms93 ) - write(iulog,'(1p,5g15.7)') etfphot_ms93(:) - write(iulog,*) '-------------------------------------------' - write(iulog,*) ' ' - end if !------------------------------------------------------------------------------ ! ... loads Chebyshev polynomial Coeff !------------------------------------------------------------------------------ @@ -211,18 +201,18 @@ subroutine get_crs( xs_short_file, sht_indexer ) do m = 1,phtcnt if( pht_alias_lst(m,1) == ' ' ) then ierr = pio_inq_varid( ncid, rxt_tag_lst(m), varid ) - if( ierr == PIO_noerr ) then + if( ierr == PIO_noerr ) then sht_indexer(m) = varid end if else if( pht_alias_lst(m,1) == 'userdefined' ) then sht_indexer(m) = -1 else ierr = pio_inq_varid( ncid, pht_alias_lst(m,1), varid ) - if( ierr == PIO_noerr ) then + if( ierr == PIO_noerr ) then sht_indexer(m) = varid else write(iulog,*) 'get_crs : ',rxt_tag_lst(m)(:len_trim(rxt_tag_lst(m))),' alias ', & - pht_alias_lst(m,1)(:len_trim(pht_alias_lst(m,1))),' not in dataset' + pht_alias_lst(m,1)(:len_trim(pht_alias_lst(m,1))),' not in dataset' call endrun end if end if @@ -252,39 +242,43 @@ subroutine get_crs( xs_short_file, sht_indexer ) ! ... allocate arrays !------------------------------------------------------------------------------ allocate( wc(nw),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'wc', nw ) end if allocate( we(nw+1),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'we', nw+1 ) end if allocate( wlintv(nw),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'wlintv', nw ) end if - allocate( etfphot(nw),stat=ierr ) - if( ierr /= 0 ) then + allocate( etfphot(nw,pcols,begchunk:endchunk),stat=ierr ) + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'etfphot', nw ) end if + etfphot = nan + allocate( bde_o2_a(nw),bde_o2_b(nw),bde_o3_a(nw),bde_o3_b(nw),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'bde_o2_a ... bde_o3_b', nw ) end if - allocate( etfphot_ms93(nw_ms93),stat=ierr ) - if( ierr /= 0 ) then + allocate( etfphot_ms93(nw_ms93,pcols,begchunk:endchunk),stat=ierr ) + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'etfphot_ms93', nw_ms93 ) end if + etfphot_ms93 = nan + allocate( xs_o2src(nw),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'xs_o2src', nw ) end if allocate( xs_o3a(nw),xs_o3b(nw),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'xs_o3a,xs_o3b', nw ) end if allocate( xs_species(nw),xs_wl(nw,nj),stat=ierr ) - if( ierr /= 0 ) then + if( ierr /= 0 ) then call alloc_err( ierr, 'get_crs', 'xs_species,xs_wl', nw*nj ) end if @@ -371,7 +365,7 @@ subroutine xs_init(xs_coef_file) !------------------------------------------------------------------------------ ! ... Dummy arguments -!------------------------------------------------------------------------------ +!------------------------------------------------------------------------------ character(len=*), intent(in) :: xs_coef_file !------------------------------------------------------------- @@ -499,17 +493,27 @@ subroutine jshort_timestep_init use mo_util, only : rebin use solar_irrad_data, only : data_nbins=>nbins, data_we => we, data_etf => sol_etf + use solar_shade, only: sun_shade + use phys_grid, only : get_ncols_p implicit none - call rebin( data_nbins, nw, data_we, we, data_etf, etfphot ) - call rebin( data_nbins, nw_ms93, data_we, we_ms, data_etf, etfphot_ms93 ) + integer :: i, c, ncols + + do c = begchunk,endchunk + ncols = get_ncols_p(c) + do i = 1,ncols + ! apply sun shade factor to input ETF + call rebin( data_nbins, nw, data_we, we, data_etf(:)*sun_shade(:,i,c), etfphot(:,i,c) ) + call rebin( data_nbins, nw_ms93, data_we, we_ms, data_etf(:)*sun_shade(:,i,c), etfphot_ms93(:,i,c) ) + end do + end do end subroutine jshort_timestep_init subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & o3cc, tlev, zkm, mw, qrs, cparg, & - lchnk, long, co2cc, scco2, do_diag ) + lchnk, icol, co2cc, scco2, do_diag ) !==============================================================================! ! Subroutine Jshort ! !==============================================================================! @@ -569,7 +573,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & !------------------------------------------------------------------------------ integer, intent(in) :: nlev ! model vertical levels integer, intent(in) :: lchnk ! chunk index - integer, intent(in) :: long ! chunk index + integer, intent(in) :: icol ! column index real(r8), intent(in) :: zen ! Zenith angle (degrees) real(r8), intent(in) :: o2_vmr(nlev) ! o2 conc (mol/mol) real(r8), intent(in) :: o3_vmr(nlev) ! o3 conc (mol/mol) @@ -719,7 +723,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & ! corrected for O2 and O3 absorption !------------------------------------------------------------------------------ do wn = 1,nw ! nw = 33 (nsrb_tot+nsrc_tot) - fnorm(:,wn) = etfphot(wn)*trans_o2(:,wn)*trans_o3(:,wn) + fnorm(:,wn) = etfphot(wn,icol,lchnk)*trans_o2(:,wn)*trans_o3(:,wn) end do !------------------------------------------------------------------------------ @@ -737,7 +741,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & !------------------------------------------------------------------------------ ! ... lyman alpha !------------------------------------------------------------------------------ - jo2_lya(:) = etfphot(1)*ro2la(:)*wlintv(1) + jo2_lya(:) = etfphot(1,icol,lchnk)*ro2la(:)*wlintv(1) wrk(1:nsrc_tot) = xs_o2src(1:nsrc_tot)*wlintv(1:nsrc_tot) & *bde_o2_a(1:nsrc_tot) @@ -747,7 +751,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & !------------------------------------------------------------------------------ if( do_diag ) then write(iulog,*) '-------------------------------------------------' - write(iulog,*) 'jshort_hrates: fnorm,wrk at long,lchnk = ',long,lchnk + write(iulog,*) 'jshort_hrates: fnorm,wrk at icol,lchnk = ',icol,lchnk write(iulog,'(1p,5g12.5)') fnorm(nlev,1:nsrc_tot) write(iulog,*) ' ' write(iulog,'(1p,5g12.5)') wrk(1:nsrc_tot) @@ -771,7 +775,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & if( do_diag ) then write(iulog,*) '-------------------------------------------------' - write(iulog,*) 'jshort_hrates: lya,bde_o2_a,qrs(nlev) at long,lchnk = ',long,lchnk + write(iulog,*) 'jshort_hrates: lya,bde_o2_a,qrs(nlev) at icol,lchnk = ',icol,lchnk write(iulog,'(1p,5g12.5)') jo2_lya(nlev),bde_o2_a(2),qrs(nlev,1) write(iulog,*) '-------------------------------------------------' end if @@ -788,7 +792,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & qrs(:,1) = qrs(:,1) + jo2_lya(:)*.53_r8*bde_o2_a(2) if( do_diag ) then write(iulog,*) '-------------------------------------------------' - write(iulog,*) 'jshort_hrates: o2(1),qrs(nlev) at long,lchnk = ',long,lchnk + write(iulog,*) 'jshort_hrates: o2(1),qrs(nlev) at icol,lchnk = ',icol,lchnk write(iulog,'(1p,5g12.5)') o2_vmr(1),qrs(nlev,1) write(iulog,*) '-------------------------------------------------' end if @@ -829,7 +833,7 @@ subroutine jshort_hrates( nlev, zen, o2_vmr, o3_vmr, o2cc, & end subroutine jshort_hrates - subroutine jshort_photo( nlev, zen, n2cc, o2cc, o3cc, & + subroutine jshort_photo( icol, lchnk, nlev, zen, n2cc, o2cc, o3cc, & nocc, tlev, zkm, jo2_sht, jno_sht, jsht ) !==============================================================================! ! Subroutine Jshort ! @@ -896,6 +900,7 @@ subroutine jshort_photo( nlev, zen, n2cc, o2cc, o3cc, & !------------------------------------------------------------------------------ ! ... dummy arguments !------------------------------------------------------------------------------ + integer, intent(in) :: icol, lchnk ! column, chunk index integer, intent(in) :: nlev ! model vertical levels real(r8), intent(in) :: zen ! Zenith angle (degrees) real(r8), intent(in) :: n2cc(nlev) ! Molecular Nitrogen conc (mol/cm^3) @@ -1045,7 +1050,7 @@ subroutine jshort_photo( nlev, zen, n2cc, o2cc, o3cc, & ! corrected for O2 and O3 absorption !------------------------------------------------------------------------------ do wn = 1,nw ! nw = 33 (nsrb_tot+nsrc_tot) - fnorm(:,wn) = etfphot(wn)*trans_o2(:,wn)*trans_o3(:,wn) + fnorm(:,wn) = etfphot(wn,icol,lchnk)*trans_o2(:,wn)*trans_o3(:,wn) end do !------------------------------------------------------------------------------ @@ -1063,7 +1068,7 @@ subroutine jshort_photo( nlev, zen, n2cc, o2cc, o3cc, & !------------------------------------------------------------------------------ ! ... Lyman Alpha !------------------------------------------------------------------------------ - jo2_lya(:) = etfphot(1)*ro2la(:)*wlintv(1) + jo2_lya(:) = etfphot(1,icol,lchnk)*ro2la(:)*wlintv(1) wrk(1:nsrc_tot) = xs_o2src(1:nsrc_tot)*wlintv(1:nsrc_tot) wrk(1) = 0._r8 @@ -1100,7 +1105,7 @@ subroutine jshort_photo( nlev, zen, n2cc, o2cc, o3cc, & !------------------------------------------------------------------------------ ! ... Derive the NO rate constant Minsch. and Siskind, JGR, 98, 20401, 1993 !------------------------------------------------------------------------------ - call calc_jno( nlev, etfphot_ms93, n2cc, o2scol, o3scol, & + call calc_jno( nlev, etfphot_ms93(:,icol,lchnk), n2cc, o2scol, o3scol, & noscol, jno_sht ) !------------------------------------------------------------------------------ @@ -1227,7 +1232,7 @@ subroutine sphers( nlev, z, zenith_angle, dsdh, nid ) !------------------------------------------------------------------------------ ! Find index of layer in which the screening height lies !------------------------------------------------------------------------------ - id = i + id = i if( zenith_angle > 90._r8 ) then do j = 1,nlayer if( rpsinz < (zd(j-1) + re) .and. rpsinz >= (zd(j) + re) ) then @@ -1236,7 +1241,7 @@ subroutine sphers( nlev, z, zenith_angle, dsdh, nid ) end if end do end if - + do j = 1,id sm = 1._r8 if( j == id .and. id == i .and. zenith_angle > 90._r8 ) then @@ -1329,7 +1334,7 @@ subroutine slant_col( nlev, delz, dsdh, nid, absden, scol ) ! height needs to be increased for higher model top !------------------------------------------------------------------------------ if (nlev==pver) then - if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then hscale = 20.e5_r8 else hscale = 10.e5_r8 @@ -1384,7 +1389,7 @@ subroutine lymana( nlev, o2scol, rm, ro2 ) ! xso2la - REAL, molecular absorption cross section in LA bands (O) ! !-----------------------------------------------------------------------------! ! EDIT HISTORY: ! -! 01/15/2002 Taken from Sasha Madronich's TUV Version 4.1a, Doug Kinnison ! ! +! 01/15/2002 Taken from Sasha Madronich's TUV Version 4.1a, Doug Kinnison ! ! 01/15/2002 Upgraded to F90, DK ! !-----------------------------------------------------------------------------! @@ -1501,7 +1506,7 @@ subroutine calc_o2srb( nlev, nid, o2col, tlev, tsrb, xscho2 ) kbot = k end if end do - + if ( kbot == nlev ) then tsrb(:,:) = 0._r8 xscho2(:,:) = 0._r8 @@ -1511,9 +1516,9 @@ subroutine calc_o2srb( nlev, nid, o2col, tlev, tsrb, xscho2 ) ! ... Fill in cross section where X is out of range ! by repeating edge table values !------------------------------------------------------- - if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then - ! Need to be careful with nlev values for kbot and ktop. + ! Need to be careful with nlev values for kbot and ktop. ! This was handled by Hanli Liu fix. if ( kbot < nlev ) then do k = 1,kbot @@ -1544,7 +1549,7 @@ subroutine calc_o2srb( nlev, nid, o2col, tlev, tsrb, xscho2 ) !------------------------------------------------------- ! ... Calculate incremental optical depths !------------------------------------------------------- - do i = 1,nsrbtuv + do i = 1,nsrbtuv do k = 1,nlev-1 if( nid(nlev-k) /= -1 ) then !------------------------------------------------------- diff --git a/src/chemistry/mozart/mo_photo.F90 b/src/chemistry/mozart/mo_photo.F90 index a953c2db74..636e1ac5ab 100644 --- a/src/chemistry/mozart/mo_photo.F90 +++ b/src/chemistry/mozart/mo_photo.F90 @@ -7,7 +7,7 @@ module mo_photo use ppgrid, only : pcols, pver, pverp, begchunk, endchunk use cam_abortutils, only : endrun use mo_constants, only : pi,r2d,boltz,d2r - use ref_pres, only : num_pr_lev, ptop_ref + use ref_pres, only : num_pr_lev, ptop_ref use pio use cam_pio_utils, only : cam_pio_openfile use spmd_utils, only : masterproc @@ -21,7 +21,7 @@ module mo_photo public :: photo_inti, table_photo, xactive_photo public :: set_ub_col - public :: setcol + public :: setcol public :: photo_timestep_init public :: photo_register @@ -71,7 +71,7 @@ module mo_photo integer :: jhno3_ndx, jno3_ndx, jpan_ndx, jmpan_ndx integer :: jo1da_ndx, jo3pa_ndx, jno2a_ndx, jn2o5a_ndx, jn2o5b_ndx - integer :: jhno3a_ndx, jno3a_ndx, jpana_ndx, jmpana_ndx, jho2no2a_ndx + integer :: jhno3a_ndx, jno3a_ndx, jpana_ndx, jmpana_ndx, jho2no2a_ndx integer :: jonitra_ndx integer :: jppi_ndx, jepn1_ndx, jepn2_ndx, jepn3_ndx, jepn4_ndx, jepn6_ndx @@ -85,7 +85,7 @@ module mo_photo contains - + !---------------------------------------------------------------------- !---------------------------------------------------------------------- subroutine photo_register @@ -94,7 +94,7 @@ subroutine photo_register ! add photo-ionization rates to phys buffer for waccmx ionosphere module - call pbuf_add_field('IonRates' , 'physpkg', dtype_r8, (/pcols,pver,nIonRates/), ion_rates_idx) ! Ionization rates for O+,O2+,N+,N2+,NO+ + call pbuf_add_field('IonRates' , 'physpkg', dtype_r8, (/pcols,pver,nIonRates/), ion_rates_idx) ! Ionization rates for O+,O2+,N+,N2+,NO+ endsubroutine photo_register @@ -110,7 +110,7 @@ subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & use mo_photoin, only : photoin_inti use mo_tuv_inti, only : tuv_inti use mo_tuv_inti, only : nlng - use mo_seto2, only : o2_xsect_inti + use mo_seto2, only : o2_xsect_inti use interpolate_data, only: lininterp_init, lininterp, lininterp_finish, interp_type use chem_mods, only : phtcnt use chem_mods, only : ncol_abs => nabscol @@ -122,7 +122,7 @@ subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & use seasalt_model, only : sslt_names=>seasalt_names, sslt_ncnst=>seasalt_nbin use mo_jshort, only : jshort_init use mo_jeuv, only : jeuv_init, neuv - use phys_grid, only : get_ncols_p, get_rlat_all_p + use phys_grid, only : get_ncols_p, get_rlat_all_p use solar_irrad_data,only : has_spectrum use photo_bkgrnd, only : photo_bkgrnd_init use cam_history, only : addfld @@ -137,7 +137,7 @@ subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & character(len=*), intent(in) :: tuv_xsect_file character(len=*), intent(in) :: o2_xsect_file logical, intent(in) :: xactive_prates - ! waccm + ! waccm character(len=*), intent(in) :: xs_coef_file character(len=*), intent(in) :: xs_short_file character(len=*), intent(in) :: photon_file @@ -176,7 +176,7 @@ subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & !---------------------------------------------------------------------------- ! Need a larger maximum zenith angle for WACCM-X extended to high altitudes !---------------------------------------------------------------------------- - if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then max_zen_angle = 116._r8 else if ( ptop_ref < 10._r8 ) then max_zen_angle = 97.01_r8 ! degrees @@ -184,14 +184,14 @@ subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & max_zen_angle = 88.85_r8 ! degrees endif - ! jeuv_1,,, jeuv_25 --> need euv calculations --> must be waccm + ! jeuv_1,,, jeuv_25 --> need euv calculations --> must be waccm ! how to determine if shrt calc is needed ?? -- use top level pressure => waccm = true ? false if ( .not. has_spectrum ) then write(iulog,*) 'photo_inti: solar_irrad_data file needs to contain irradiance spectrum' call endrun('photo_inti: ERROR -- solar irradiance spectrum is missing') endif - + !---------------------------------------------------------------------- ! ... allocate indexers !---------------------------------------------------------------------- @@ -287,7 +287,7 @@ subroutine photo_inti( xs_coef_file, xs_short_file, xs_long_file, rsf_file, & end if do_jshort = o_ndx>0 .and. o2_ndx>0 .and. (o3_ndx>0.or.o3_inv_ndx>0) .and. n2_ndx>0 .and. no_ndx>0 - + call jeuv_init( photon_file, electron_file, euv_indexer ) do_jeuv = any(euv_indexer(:)>0) @@ -639,7 +639,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & real(r8), allocatable :: lng_prates(:,:) ! photorates matrix (1/s) real(r8), allocatable :: sht_prates(:,:) ! photorates matrix (1/s) real(r8), allocatable :: euv_prates(:,:) ! photorates matrix (1/s) - + real(r8), allocatable :: zarg(:) real(r8), allocatable :: tline(:) ! vertical temperature array @@ -676,7 +676,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & if ((.not.do_jshort) .or. (ptop_ref < 10._r8)) then n_jshrt_levs = pver - p1 = 1 + p1 = 1 p2 = pver else n_jshrt_levs = pver+1 @@ -708,7 +708,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & end if endif endif - + if (nsht>0) then allocate( sht_prates(n_jshrt_levs,nsht),stat=astat ) if( astat /= 0 ) then @@ -781,7 +781,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & lwc_line(:) = lwc(i,:) cld_line(:) = clouds(i,:) - + tline(p1:p2) = temper(i,:pver) zarg(p1:p2) = zmid(i,:pver) @@ -799,7 +799,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & if (jpni3_ndx > 0 ) photos(i,:,jpni3_ndx) = photos(i,:,jpni3_ndx) + esfact * 0.15_r8 if (jpni4_ndx > 0 ) photos(i,:,jpni4_ndx) = photos(i,:,jpni4_ndx) + esfact * 6.2e-3_r8 ! added to v02 - if (jpni5_ndx > 0 ) photos(i,:,jpni5_ndx) = photos(i,:,jpni5_ndx) + esfact * 1.0_r8 + if (jpni5_ndx > 0 ) photos(i,:,jpni5_ndx) = photos(i,:,jpni5_ndx) + esfact * 1.0_r8 endif if (do_jshort) then if ( ptop_ref > 10._r8 ) then @@ -833,8 +833,8 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & !----------------------------------------------------------------- ! ... short wave length component !----------------------------------------------------------------- - call jshort( n_jshrt_levs, sza, n2_den, o2_den, o3_den, & - no_den, tline, zarg, jo2_sht, jno_sht, sht_prates ) + call jshort( i, lchnk, n_jshrt_levs, sza, n2_den, o2_den, o3_den, & + no_den, tline, zarg, jo2_sht, jno_sht, sht_prates ) do m = 1,phtcnt if( sht_indexer(m) > 0 ) then @@ -874,7 +874,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & !----------------------------------------------------------------- ! ... long wave length component !----------------------------------------------------------------- - call jlong( pver, sza, eff_alb, parg, tline, colo3, lng_prates ) + call jlong( i, lchnk, pver, sza, eff_alb, parg, tline, colo3, lng_prates ) do m = 1,phtcnt if( lng_indexer(m) > 0 ) then alias_factor = pht_alias_mult(m,2) @@ -906,9 +906,9 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & ! Save photo-ionization rates to physics buffer accessed in ionosphere module for WACCMX if (ion_rates_idx>0) then - + ionRates(i,1:pver,1:nIonRates) = esfact * euv_prates(1:pver,1:nIonRates) - + endif end if daylight @@ -924,7 +924,7 @@ subroutine table_photo( photos, pmid, pdel, temper, zmid, zint, & endif end do col_loop - + if ( do_jeuv ) then qbktot(:ncol,:) = qbko1(:ncol,:)+qbko2(:ncol,:)+qbkn2(:ncol,:)+qbkn1(:ncol,:)+qbkno(:ncol,:) call outfld( 'Qbkgndtot', qbktot(:ncol,:),ncol, lchnk ) @@ -1026,15 +1026,15 @@ subroutine xactive_photo( photos, vmr, temper, cwat, cldfr, & real(r8) :: xfrc(pverp) ! cloud fraction xuexi real(r8) :: airdens(pverp) ! atmospheric density real(r8) :: o3line(pverp) ! vertical o3 vmr - real(r8) :: aerocs1(pverp) - real(r8) :: aerocs2(pverp) - real(r8) :: aercbs1(pverp) - real(r8) :: aercbs2(pverp) - real(r8) :: aersoa(pverp) - real(r8) :: aerant(pverp) - real(r8) :: aerso4(pverp) + real(r8) :: aerocs1(pverp) + real(r8) :: aerocs2(pverp) + real(r8) :: aercbs1(pverp) + real(r8) :: aercbs2(pverp) + real(r8) :: aersoa(pverp) + real(r8) :: aerant(pverp) + real(r8) :: aerso4(pverp) real(r8) :: aerds(4,pverp) - real(r8) :: rh(pverp) + real(r8) :: rh(pverp) real(r8) :: zarg(pverp) ! vertical height array real(r8) :: aersal(pverp,4) real(r8) :: albedo(kw) ! wavelenght dependent albedo @@ -1169,7 +1169,7 @@ subroutine xactive_photo( photos, vmr, temper, cwat, cldfr, & do ndx = 1,4 aerds(ndx,pverp:2:-1) = dust_vmr(i,:,ndx) end do - else + else do ndx = 1,4 aerds(ndx,pverp:2:-1) = 0._r8 end do @@ -1205,8 +1205,8 @@ subroutine xactive_photo( photos, vmr, temper, cwat, cldfr, & airdens, aerocs1, aerocs2, aercbs1, aercbs2, & aersoa, aerant, aerso4, aersal, aerds, o3line, rh, & prates, sza, nw, dt_xdiag ) - dt_diag(i,:) = dt_xdiag(:) - + dt_diag(i,:) = dt_xdiag(:) + do m = 1,phtcnt if( lng_indexer(m) > 0 ) then alias_factor = pht_alias_mult(m,2) @@ -1436,7 +1436,7 @@ subroutine set_ub_col( col_delta, vmr, invariants, ptop, pdel, ncol, lchnk ) real(r8) :: o2_exo_col(ncol) real(r8) :: o3_exo_col(ncol) integer :: i - + !--------------------------------------------------------------- ! ... assign column density at the upper boundary ! the first column is o3 and the second is o2. @@ -1458,8 +1458,8 @@ subroutine set_ub_col( col_delta, vmr, invariants, ptop, pdel, ncol, lchnk ) + delp * (o2_exo_coldens(ki,i,lchnk,next) & - o2_exo_coldens(kl,i,lchnk,next)) else - tint_vals(1) = o2_exo_coldens( 1,i,lchnk,last) - tint_vals(2) = o2_exo_coldens( 1,i,lchnk,next) + tint_vals(1) = o2_exo_coldens( 1,i,lchnk,last) + tint_vals(2) = o2_exo_coldens( 1,i,lchnk,next) endif o2_exo_col(i) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1)) end do @@ -1476,8 +1476,8 @@ subroutine set_ub_col( col_delta, vmr, invariants, ptop, pdel, ncol, lchnk ) + delp * (o3_exo_coldens(ki,i,lchnk,next) & - o3_exo_coldens(kl,i,lchnk,next)) else - tint_vals(1) = o3_exo_coldens( 1,i,lchnk,last) - tint_vals(2) = o3_exo_coldens( 1,i,lchnk,next) + tint_vals(1) = o3_exo_coldens( 1,i,lchnk,last) + tint_vals(2) = o3_exo_coldens( 1,i,lchnk,next) endif o3_exo_col(i) = tint_vals(1) + dels * (tint_vals(2) - tint_vals(1)) end do @@ -1711,7 +1711,7 @@ subroutine photo_timestep_init( calday ) ! Set jlong etf !----------------------------------------------------------------------- call jlong_timestep_init - + if ( do_jshort ) then !----------------------------------------------------------------------- ! Set jshort etf diff --git a/src/chemistry/mozart/mo_waccm_hrates.F90 b/src/chemistry/mozart/mo_waccm_hrates.F90 index 5e7b8a8405..b0cf007162 100644 --- a/src/chemistry/mozart/mo_waccm_hrates.F90 +++ b/src/chemistry/mozart/mo_waccm_hrates.F90 @@ -26,7 +26,7 @@ module mo_waccm_hrates integer :: ele_temp_ndx, ion_temp_ndx contains - + subroutine init_hrates( ) use mo_chem_utls, only : get_spc_ndx use cam_history, only : addfld @@ -81,7 +81,7 @@ subroutine init_hrates( ) attr = 'total jo2 euv photolysis rate' call addfld( 'JO2_EUV', (/ 'lev' /), 'I', '/s', trim(attr) ) - ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index + ele_temp_ndx = pbuf_get_index('TElec',errcode=err)! electron temperature index ion_temp_ndx = pbuf_get_index('TIon',errcode=err) ! ion temperature index end subroutine init_hrates @@ -208,19 +208,19 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) qrs_tot(:ncol,:) = 0._r8 if (.not. has_hrates) return - -!------------------------------------------------------------------------- + +!------------------------------------------------------------------------- ! ... set maximum zenith angle - higher value for higher top model -!------------------------------------------------------------------------- - if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then +!------------------------------------------------------------------------- + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then max_zen_angle = 116._r8 else max_zen_angle = 97.01_r8 ! degrees endif -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... get chunk latitudes and longitudes -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- lchnk = state%lchnk call get_lat_all_p( lchnk, ncol, latndx ) @@ -228,23 +228,23 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) call get_rlat_all_p( lchnk, ncol, rlats ) call get_rlon_all_p( lchnk, ncol, rlons ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set lower limit for heating rates which is now dictated by radheat module -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- kbot_hrates = bot_mlt_lev kbot_hrates = min( kbot_hrates,pver ) ! write(iulog,*) 'hrates: kbot_hrates = ',kbot_hrates -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... calculate cosine of zenith angle then cast back to angle -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- calday = get_curr_calday() call zenith( calday, rlats, rlons, zen_angle, ncol ) zen_angle(:) = acos( zen_angle(:) ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... map incoming concentrations to working array -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do m = 1,pcnst n = map2chm(m) if( n > 0 ) then @@ -255,49 +255,49 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) end do call get_short_lived_species( mmr, lchnk, ncol, pbuf ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set atmosphere mean mass -!----------------------------------------------------------------------- - if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then +!----------------------------------------------------------------------- + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then do k = 1,pver mbar(:ncol,k) = mbarv(:ncol,k,lchnk) enddo - else + else call set_mean_mass( ncol, mmr, mbar ) endif ! -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... xform from mmr to vmr -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call mmr2vmr( mmr(:ncol,:,:), vmr(:ncol,:,:), mbar(:ncol,:), ncol ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... xform water vapor from mmr to vmr -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do k = 1,pver h2ovmr(:ncol,k) = vmr(:ncol,k,id_h2o) end do -!----------------------------------------------------------------------- -! ... xform geopotential height from m to km +!----------------------------------------------------------------------- +! ... xform geopotential height from m to km ! and pressure from Pa to mb -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- zsurf(:ncol) = rga * state%phis(:ncol) do k = 1,pver zmid(:ncol,k) = m2km * (state%zm(:ncol,k) + zsurf(:ncol)) end do -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set the "invariants" -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call setinv( invariants, state%t, h2ovmr, vmr, state%pmid, ncol, lchnk, pbuf ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set the column densities at the upper boundary -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call set_ub_col( col_delta, vmr, invariants, state%pint(:,1), state%pdel, ncol, lchnk ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set rates for "tabular" and user specified reactions -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do m = 1,rxntot do k = 1,pver reaction_rates(:,k,m) = 0._r8 @@ -307,15 +307,15 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) call usrrxt_hrates( reaction_rates, state%t, ele_temp_fld, ion_temp_fld, & h2ovmr, invariants(:,:,indexm), ncol, kbot_hrates ) call adjrxt( reaction_rates, invariants, invariants(1,1,indexm), ncol,pver ) - -!----------------------------------------------------------------------- + +!----------------------------------------------------------------------- ! ... set cp array -!----------------------------------------------------------------------- - if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then +!----------------------------------------------------------------------- + if ( waccmx_is('ionosphere') .or. waccmx_is('neutral') ) then do k = 1, pver cpair(:ncol,k) = cpairv(:ncol,k,lchnk) enddo - else + else call calc_cp( ncol, vmr, cpair ) endif @@ -329,18 +329,18 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) write(iulog,*) ' ' #endif -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set the earth-sun distance factor -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call shr_orb_decl( calday, eccen, mvelpp, lambm0, obliqr , & delta, esfact ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... set the column densities -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call setcol( col_delta, col_dens, vmr, state%pdel, ncol ) !----------------------------------------------------------------------- ! ... compute the thermal heating rates -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do m = 1,4 do k = 1,pver qrs(:,k,m) = 0._r8 @@ -376,7 +376,7 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) call jshort( pver, sza, o2_line, o3_line, o2cc, & o3cc, tline, zarg, mw, qrs_col, & cparg, lchnk, i, co2cc, scco2, do_diag ) - call jlong( pver, sza, eff_alb, parg, tline, & + call jlong( i, lchnk, pver, sza, eff_alb, parg, tline, & mw, o2_line, o3_line, colo3, qrl_col, & cparg, kbot_hrates ) do m = 1,4 @@ -419,15 +419,15 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) call outfld( 'QRS_EUV', euv_hrate(:,:), ncol, lchnk ) call outfld( 'QRS_CO2NIR', co2_hrate(:,:), ncol, lchnk ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... chemical pot heating rate -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call cph( cphrate, vmr, reaction_rates, cpair, mbar, & kbot_hrates, ncol, lchnk ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... auroral ion production -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call aurora( state%t, mbar, rlats, & aur_hrate, cpair, state%pmid, lchnk, calday, & ncol, rlons, pbuf ) @@ -436,15 +436,15 @@ subroutine waccm_hrates(ncol, state, asdir, bot_mlt_lev, qrs_tot, pbuf ) end do call outfld( 'QRS_AUR', aur_hrate(:,:), ncol, lchnk ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... airglow heating rate -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- call airglow( aghrate, vmr(1,1,id_o2_1s), vmr(1,1,id_o2_1d), vmr(1,1,id_o1d), reaction_rates, cpair, & ncol, lchnk ) -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- ! ... form total heating rate -!----------------------------------------------------------------------- +!----------------------------------------------------------------------- do k = 1,kbot_hrates qrs_tot(:ncol,k) = qrs(:,k,1) + qrs(:,k,2) + qrs(:,k,3) + qrs(:,k,4) & + qrl(:,k,1) + qrl(:,k,2) + qrl(:,k,3) + qrl(:,k,4) diff --git a/src/chemistry/utils/solar_shade.F90 b/src/chemistry/utils/solar_shade.F90 new file mode 100644 index 0000000000..ad416fb077 --- /dev/null +++ b/src/chemistry/utils/solar_shade.F90 @@ -0,0 +1,99 @@ +module solar_shade + use shr_kind_mod, only: r8 => shr_kind_r8 + use solar_irrad_data, only: nbins, we + use physics_types,only : physics_state + use ppgrid, only : pcols, begchunk, endchunk + use cam_logfile, only: iulog + use spmd_utils, only: masterproc + use cam_abortutils, only: endrun + + implicit none + private + public :: solar_shade_readnl + public :: solar_shade_init + + real(r8), public, protected, allocatable :: sun_shade(:,:,:) ! chunk, col, wavelen dependent + + real(r8) :: soll0 = 1._r8 ! namelist variable + real(r8) :: soll1 = 1._r8 ! namelist variable + real(r8) :: soll2 = 1._r8 ! namelist variable + +contains + + !============================================================================ + subroutine solar_shade_readnl(nlfile) + use namelist_utils, only: find_group_name + use units, only: getunit, freeunit + use spmd_utils, only: mpicom, masterprocid, mpi_real8, mpi_success + + ! arguments + character(len=*), intent(in) :: nlfile ! filepath for file containing namelist input + + ! local vars + integer :: unitn, ierr + + character(len=*), parameter :: prefix = 'solar_shade_readnl: ' + + namelist /solar_shade_opts/ soll0, soll1, soll2 + + if (masterproc) then + unitn = getunit() + open( unitn, file=trim(nlfile), status='old' ) + call find_group_name(unitn, 'solar_shade_opts', status=ierr) + if (ierr == 0) then + read(unitn, solar_shade_opts, iostat=ierr) + if (ierr /= 0) then + call endrun(prefix//'ERROR reading namelist') + end if + end if + close(unitn) + call freeunit(unitn) + end if + + call mpi_bcast(soll0, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr/=mpi_success) then + call endrun(prefix//'ERROR bcast soll0') + end if + call mpi_bcast(soll1, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr/=mpi_success) then + call endrun(prefix//'ERROR bcast soll1') + end if + call mpi_bcast(soll2, 1, mpi_real8, masterprocid, mpicom, ierr) + if (ierr/=mpi_success) then + call endrun(prefix//'ERROR bcast soll2') + end if + + if (masterproc) then + write(iulog,*) prefix,'factor soll0 = ', soll0 + write(iulog,*) prefix,'factor soll1 = ', soll1 + write(iulog,*) prefix,'factor soll2 = ', soll2 + end if + + end subroutine solar_shade_readnl + + !============================================================================ + subroutine solar_shade_init() + use phys_grid, only : get_ncols_p, get_rlat_all_p, get_rlon_all_p + + integer :: lchnk, ncol, icol + real(r8) :: clat(pcols) ! current latitudes(radians) + real(r8) :: modval + + allocate(sun_shade(nbins, pcols, begchunk:endchunk)) + sun_shade = 1._r8 + + ! copy of Ben's insolation function + do lchnk=begchunk,endchunk + ncol = get_ncols_p(lchnk) + call get_rlat_all_p(lchnk, ncol, clat) + do icol = 1,ncol + modval = (1._r8-soll0) + modval = modval + (1._r8-soll1)*sin(clat(icol)) + modval = modval + (1._r8-soll2)*0.5_r8*(3._r8*(sin(clat(icol)))**2-1._r8) + sun_shade(:nbins, icol, lchnk) = 1._r8-modval + end do + end do + + end subroutine solar_shade_init + +end module solar_shade diff --git a/src/control/runtime_opts.F90 b/src/control/runtime_opts.F90 index 8018c3a1a3..b8b73e7753 100644 --- a/src/control/runtime_opts.F90 +++ b/src/control/runtime_opts.F90 @@ -70,6 +70,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) use prescribed_strataero,only: prescribed_strataero_readnl use aerodep_flx, only: aerodep_flx_readnl use solar_data, only: solar_data_readnl + use solar_shade, only: solar_shade_readnl use tropopause, only: tropopause_readnl use aoa_tracers, only: aoa_tracers_readnl use prescribed_ozone, only: prescribed_ozone_readnl @@ -162,6 +163,7 @@ subroutine read_namelist(nlfilename, single_column, scmlat, scmlon) call prescribed_volcaero_readnl(nlfilename) call prescribed_strataero_readnl(nlfilename) call solar_data_readnl(nlfilename) + call solar_shade_readnl(nlfilename) call carma_readnl(nlfilename) call tropopause_readnl(nlfilename) call aoa_tracers_readnl(nlfilename) diff --git a/src/physics/cam/physpkg.F90 b/src/physics/cam/physpkg.F90 index 218935f6aa..588aea4ac0 100644 --- a/src/physics/cam/physpkg.F90 +++ b/src/physics/cam/physpkg.F90 @@ -742,6 +742,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out ) use sslt_rebin, only: sslt_rebin_init use tropopause, only: tropopause_init use solar_data, only: solar_data_init + use solar_shade, only: solar_shade_init use dadadj_cam, only: dadadj_init use cam_abortutils, only: endrun use nudging, only: Nudge_Model, nudging_init @@ -820,6 +821,7 @@ subroutine phys_init( phys_state, phys_tend, pbuf2d, cam_out ) ! solar irradiance data modules call solar_data_init() + call solar_shade_init() ! Prognostic chemistry. call chem_init(phys_state,pbuf2d) diff --git a/src/physics/rrtmg/aer_src/rrtmg_sw_rad.f90 b/src/physics/rrtmg/aer_src/rrtmg_sw_rad.f90 index 450c64fac6..b2e31a160e 100644 --- a/src/physics/rrtmg/aer_src/rrtmg_sw_rad.f90 +++ b/src/physics/rrtmg/aer_src/rrtmg_sw_rad.f90 @@ -72,7 +72,7 @@ module rrtmg_sw_rad ! Set iaer to select aerosol option ! iaer = 0, no aerosols -! iaer = 10, input total aerosol optical depth, single scattering albedo +! iaer = 10, input total aerosol optical depth, single scattering albedo ! and asymmetry parameter (tauaer, ssaaer, asmaer) directly integer, parameter :: iaer = 10 @@ -104,28 +104,28 @@ subroutine rrtmg_sw & ! ------- Description ------- -! This program is the driver for RRTMG_SW, the AER SW radiation model for +! This program is the driver for RRTMG_SW, the AER SW radiation model for ! application to GCMs, that has been adapted from RRTM_SW for improved ! efficiency and to provide fractional cloudiness and cloud overlap ! capability using McICA. ! ! This routine ! b) calls INATM_SW to read in the atmospheric profile; -! all layering in RRTMG is ordered from surface to toa. +! all layering in RRTMG is ordered from surface to toa. ! c) calls CLDPRMC_SW to set cloud optical depth for McICA based ! on input cloud properties -! d) calls SETCOEF_SW to calculate various quantities needed for +! d) calls SETCOEF_SW to calculate various quantities needed for ! the radiative transfer algorithm -! e) calls SPCVMC to call the two-stream model that in turn -! calls TAUMOL to calculate gaseous optical depths for each +! e) calls SPCVMC to call the two-stream model that in turn +! calls TAUMOL to calculate gaseous optical depths for each ! of the 16 spectral bands and to perform the radiative transfer ! using McICA, the Monte-Carlo Independent Column Approximation, ! to represent sub-grid scale cloud variability ! f) passes the calculated fluxes and cooling rates back to GCM ! ! *** This version uses McICA *** -! Monte Carlo Independent Column Approximation (McICA, Pincus et al., -! JC, 2003) method is applied to the forward model calculation +! Monte Carlo Independent Column Approximation (McICA, Pincus et al., +! JC, 2003) method is applied to the forward model calculation ! This method is valid for clear sky or partial cloud conditions. ! ! This call to RRTMG_SW must be preceeded by a call to the module @@ -134,7 +134,7 @@ subroutine rrtmg_sw & ! on the RRTMG quadrature point (ngptsw) dimension. ! ! *** This version only allows input of cloud optical properties *** -! Input cloud fraction, cloud optical depth, single scattering albedo +! Input cloud fraction, cloud optical depth, single scattering albedo ! and asymmetry parameter directly (inflg = 0) ! ! *** This version only allows input of aerosol optical properties *** @@ -145,7 +145,7 @@ subroutine rrtmg_sw & ! ------- Modifications ------- ! ! This version of RRTMG_SW has been modified from RRTM_SW to use a reduced -! set of g-point intervals and a two-stream model for application to GCMs. +! set of g-point intervals and a two-stream model for application to GCMs. ! !-- Original version (derived from RRTM_SW) ! 2002: AER. Inc. @@ -154,10 +154,10 @@ subroutine rrtmg_sw & !-- Additional modifications for GCM application ! Aug 2003: M. J. Iacono, AER Inc. !-- Total number of g-points reduced from 224 to 112. Original -! set of 224 can be restored by exchanging code in module parrrsw.f90 +! set of 224 can be restored by exchanging code in module parrrsw.f90 ! and in file rrtmg_sw_init.f90. ! Apr 2004: M. J. Iacono, AER, Inc. -!-- Modifications to include output for direct and diffuse +!-- Modifications to include output for direct and diffuse ! downward fluxes. There are output as "true" fluxes without ! any delta scaling applied. Code can be commented to exclude ! this calculation in source file rrtmg_sw_spcvrt.f90. @@ -166,7 +166,7 @@ subroutine rrtmg_sw & ! Nov 2005: M. J. Iacono, AER, Inc. !-- Reformatted for consistency with rrtmg_lw. ! Feb 2007: M. J. Iacono, AER, Inc. -!-- Modifications to formatting to use assumed-shape arrays. +!-- Modifications to formatting to use assumed-shape arrays. ! Aug 2007: M. J. Iacono, AER, Inc. !-- Modified to output direct and diffuse fluxes either with or without ! delta scaling based on setting of idelm flag @@ -178,7 +178,7 @@ subroutine rrtmg_sw & ! ----- Input ----- integer, intent(in) :: lchnk ! chunk identifier - integer, intent(in) :: ncol ! Number of horizontal columns + integer, intent(in) :: ncol ! Number of horizontal columns integer, intent(in) :: nlay ! Number of model layers integer, intent(in) :: icld ! Cloud overlap method ! 0: Clear only @@ -221,7 +221,7 @@ subroutine rrtmg_sw & real(kind=r8), intent(in) :: adjes ! Flux adjustment for Earth/Sun distance real(kind=r8), intent(in) :: coszen(:) ! Cosine of solar zenith angle ! Dimensions: (ncol) - real(kind=r8), intent(in) :: solvar(1:nbndsw) ! Solar constant (Wm-2) scaling per band + real(kind=r8), intent(in) :: solvar(1:nbndsw,1:ncol) ! Solar constant (Wm-2) scaling per band real(kind=r8), intent(in) :: cldfmcl(:,:,:) ! Cloud fraction ! Dimensions: (ngptsw,ncol,nlay) @@ -243,13 +243,13 @@ subroutine rrtmg_sw & ! Dimensions: (ncol,nlay) real(kind=r8), intent(in) :: tauaer(:,:,:) ! Aerosol optical depth (iaer=10 only) ! Dimensions: (ncol,nlay,nbndsw) - ! (non-delta scaled) + ! (non-delta scaled) real(kind=r8), intent(in) :: ssaaer(:,:,:) ! Aerosol single scattering albedo (iaer=10 only) ! Dimensions: (ncol,nlay,nbndsw) - ! (non-delta scaled) + ! (non-delta scaled) real(kind=r8), intent(in) :: asmaer(:,:,:) ! Aerosol asymmetry parameter (iaer=10 only) ! Dimensions: (ncol,nlay,nbndsw) - ! (non-delta scaled) + ! (non-delta scaled) ! ----- Output ----- @@ -371,13 +371,13 @@ subroutine rrtmg_sw & real(kind=r8) :: znifu(ncol,nlay+2) ! temporary near-IR downward shortwave flux (w/m2) real(kind=r8) :: znicu(ncol,nlay+2) ! temporary clear sky near-IR downward shortwave flux (w/m2) - ! Optional output fields + ! Optional output fields real(kind=r8) :: swnflx(nlay+2) ! Total sky shortwave net flux (W/m2) real(kind=r8) :: swnflxc(nlay+2) ! Clear sky shortwave net flux (W/m2) real(kind=r8) :: dirdflux(nlay+2) ! Direct downward shortwave surface flux real(kind=r8) :: difdflux(nlay+2) ! Diffuse downward shortwave surface flux - real(kind=r8) :: uvdflx(nlay+2) ! Total sky downward shortwave flux, UV/vis - real(kind=r8) :: nidflx(nlay+2) ! Total sky downward shortwave flux, near-IR + real(kind=r8) :: uvdflx(nlay+2) ! Total sky downward shortwave flux, UV/vis + real(kind=r8) :: nidflx(nlay+2) ! Total sky downward shortwave flux, near-IR ! Initializations @@ -588,7 +588,7 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & taua, ssaa, asma) ! Input atmospheric profile from GCM, and prepare it for use in RRTMG_SW. - ! Set other RRTMG_SW input parameters. + ! Set other RRTMG_SW input parameters. use parrrsw, only: nbndsw, ngptsw, nmol, mxmol, & jpband, jpb1, jpb2 @@ -626,7 +626,7 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & integer, intent(in) :: dyofyr ! Day of the year (used to get Earth/Sun ! distance if adjflx not provided) real(kind=r8), intent(in) :: adjes ! Flux adjustment for Earth/Sun distance - real(kind=r8), intent(in) :: solvar(jpb1:jpb2) ! Solar constant (Wm-2) scaling per band + real(kind=r8), intent(in) :: solvar(jpb1:jpb2,1:ncol) ! Solar constant (Wm-2) scaling per band real(kind=r8), intent(in) :: cldfmcl(:,:,:) ! Cloud fraction ! Dimensions: (ngptsw,ncol,nlay) @@ -719,7 +719,7 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & real(kind=r8), parameter :: sbc = 5.67e-08_r8 ! Stefan-Boltzmann constant (W/m2K4) integer :: isp, l, ix, n, imol, ib, ig, iplon ! Loop indices - real(kind=r8) :: amm, summol ! + real(kind=r8) :: amm, summol ! real(kind=r8) :: adjflx ! flux adjustment for Earth/Sun distance !----------------------------------------------------------------------------------------- @@ -728,16 +728,16 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & adjflx = adjes ! 2) Calculate Earth/Sun distance from DYOFYR, the cumulative day of the year. - ! (Set adjflx to 1. to use constant Earth/Sun distance of 1 AU). + ! (Set adjflx to 1. to use constant Earth/Sun distance of 1 AU). if (dyofyr .gt. 0) then adjflx = earth_sun(dyofyr) endif ! Set incoming solar flux adjustment to include adjustment for ! current Earth/Sun distance (ADJFLX) and scaling of default internal - ! solar constant (rrsw_scon = 1368.22 Wm-2) by band (SOLVAR). SOLVAR can be set - ! to a single scaling factor as needed, or to a different value in each - ! band, which may be necessary for paleoclimate simulations. + ! solar constant (rrsw_scon = 1368.22 Wm-2) by band (SOLVAR). SOLVAR can be set + ! to a single scaling factor as needed, or to a different value in each + ! band, which may be necessary for paleoclimate simulations. do iplon = 1 ,ncol adjflux(iplon,:) = 0._r8 @@ -745,7 +745,7 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & do ib = jpb1,jpb2 do iplon = 1, ncol - adjflux(iplon,ib) = adjflx * solvar(ib) + adjflux(iplon,ib) = adjflx * solvar(ib,iplon) end do end do @@ -754,15 +754,15 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & tbound(iplon) = tsfc(iplon) ! Install input GCM arrays into RRTMG_SW arrays for pressure, temperature, - ! and molecular amounts. + ! and molecular amounts. ! Pressures are input in mb, or are converted to mb here. - ! Molecular amounts are input in volume mixing ratio, or are converted from + ! Molecular amounts are input in volume mixing ratio, or are converted from ! mass mixing ratio (or specific humidity for h2o) to volume mixing ratio - ! here. These are then converted to molecular amount (molec/cm2) below. - ! The dry air column COLDRY (in molec/cm2) is calculated from the level - ! pressures, pz (in mb), based on the hydrostatic equation and includes a - ! correction to account for h2o in the layer. The molecular weight of moist - ! air (amm) is calculated for each layer. + ! here. These are then converted to molecular amount (molec/cm2) below. + ! The dry air column COLDRY (in molec/cm2) is calculated from the level + ! pressures, pz (in mb), based on the hydrostatic equation and includes a + ! correction to account for h2o in the layer. The molecular weight of moist + ! air (amm) is calculated for each layer. ! Note: In RRTMG, layer indexing goes from bottom to top, and coding below ! assumes GCM input fields are also bottom to top. Input layer indexing ! from GCM fields should be reversed here if necessary. @@ -790,8 +790,8 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & wkl(iplon,4,l) = n2ovmr(iplon,nlay-l+1) wkl(iplon,5,l) = 0._r8 wkl(iplon,6,l) = ch4vmr(iplon,nlay-l+1) - wkl(iplon,7,l) = o2vmr(iplon,nlay-l+1) - amm = (1._r8 - wkl(iplon,1,l)) * amd + wkl(iplon,1,l) * amw + wkl(iplon,7,l) = o2vmr(iplon,nlay-l+1) + amm = (1._r8 - wkl(iplon,1,l)) * amd + wkl(iplon,1,l) * amw coldry(iplon,l) = (pz(iplon,l-1)-pz(iplon,l)) * 1.e3_r8 * avogad / & (1.e2_r8 * grav * amm * (1._r8 + wkl(iplon,1,l))) end do @@ -812,7 +812,7 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & ! Transfer aerosol optical properties to RRTM variables; ! modify to reverse layer indexing here if necessary. - if (iaer .ge. 1) then + if (iaer .ge. 1) then do l = 1, nlay-1 do ib = 1, nbndsw do iplon = 1, ncol @@ -827,7 +827,7 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & ! Transfer cloud fraction and cloud optical properties to RRTM variables; ! modify to reverse layer indexing here if necessary. - if (icld .ge. 1) then + if (icld .ge. 1) then ! Move incoming GCM cloud arrays to RRTMG cloud arrays. ! For GCM input, incoming reice is in effective radius; for Fu parameterization (iceflag = 3) ! convert effective radius to generalized effective size using method of Mitchell, JAS, 2002: @@ -878,5 +878,3 @@ subroutine inatm_sw (ncol, nlay, icld, iaer, & end subroutine inatm_sw end module rrtmg_sw_rad - - diff --git a/src/physics/rrtmg/rad_solar_var.F90 b/src/physics/rrtmg/rad_solar_var.F90 index 10b98f10eb..81635f3373 100644 --- a/src/physics/rrtmg/rad_solar_var.F90 +++ b/src/physics/rrtmg/rad_solar_var.F90 @@ -1,12 +1,13 @@ !------------------------------------------------------------------------------- ! This module uses the Lean solar irradiance data to provide a solar cycle -! scaling factor used in heating rate calculations +! scaling factor used in heating rate calculations !------------------------------------------------------------------------------- module rad_solar_var use shr_kind_mod , only : r8 => shr_kind_r8 use solar_irrad_data, only : sol_irrad, we, nbins, has_spectrum, sol_tsi use solar_irrad_data, only : do_spctrl_scaling + use solar_shade, only : sun_shade use cam_abortutils, only : endrun implicit none @@ -18,7 +19,7 @@ module rad_solar_var real(r8), allocatable :: ref_band_irrad(:) ! scaling will be relative to ref_band_irrad in each band real(r8), allocatable :: irrad(:) ! solar irradiance at model timestep in each band - real(r8) :: tsi_ref ! total solar irradiance assumed by rrtmg + real(r8) :: tsi_ref ! total solar irradiance assumed by rrtmg real(r8), allocatable :: radbinmax(:) real(r8), allocatable :: radbinmin(:) @@ -89,21 +90,31 @@ end subroutine rad_solar_var_init !------------------------------------------------------------------------------- !------------------------------------------------------------------------------- - subroutine get_variability( sfac ) + subroutine get_variability( lchnk, sfac ) + use ppgrid, only : pcols + use phys_grid, only: get_ncols_p - real(r8), intent(out) :: sfac(nradbins) ! scaling factors for CAM heating + integer, intent(in) :: lchnk + real(r8), intent(out) :: sfac(nradbins,pcols) ! scaling factors for CAM heating integer :: yr, mon, day, tod + integer :: icol, ncols + + ncols = get_ncols_p(lchnk) if ( do_spctrl_scaling ) then - call integrate_spectrum( nbins, nradbins, we, radbinmin, radbinmax, sol_irrad, irrad) + sfac = 0._r8 - sfac(:nradbins) = irrad(:nradbins)/ref_band_irrad(:nradbins) + ! apply the column-dependent sun shade factor before integrating to RRTMG bands + do icol = 1,ncols + call integrate_spectrum( nbins, nradbins, we, radbinmin, radbinmax, sol_irrad(:nbins)*sun_shade(:nbins,icol,lchnk), irrad) + sfac(:nradbins,icol) = irrad(:nradbins)/ref_band_irrad(:nradbins) + end do else - sfac(:nradbins) = sol_tsi/tsi_ref + sfac(:nradbins,:ncols) = sol_tsi/tsi_ref endif @@ -129,7 +140,7 @@ subroutine integrate_spectrum( nsrc, ntrg, src_x, min_trg, max_trg, src, trg ) real(r8), intent(in) :: min_trg(ntrg) ! target coordinates real(r8), intent(in) :: src(nsrc) ! source array real(r8), intent(out) :: trg(ntrg) ! target array - + !--------------------------------------------------------------- ! ... local variables !--------------------------------------------------------------- diff --git a/src/physics/rrtmg/radiation.F90 b/src/physics/rrtmg/radiation.F90 index 9f9cfa40d1..5ef2c860b1 100644 --- a/src/physics/rrtmg/radiation.F90 +++ b/src/physics/rrtmg/radiation.F90 @@ -107,7 +107,7 @@ module radiation real(r8) :: aer_tau400(pcols,0:pver) real(r8) :: aer_tau550(pcols,0:pver) real(r8) :: aer_tau700(pcols,0:pver) - + end type rad_out_t ! Namelist variables @@ -126,19 +126,19 @@ module radiation logical :: use_rad_uniform_angle = .false. ! if true, use the namelist rad_uniform_angle for the coszrs calculation ! Physics buffer indices -integer :: qrs_idx = 0 -integer :: qrl_idx = 0 -integer :: su_idx = 0 -integer :: sd_idx = 0 -integer :: lu_idx = 0 -integer :: ld_idx = 0 +integer :: qrs_idx = 0 +integer :: qrl_idx = 0 +integer :: su_idx = 0 +integer :: sd_idx = 0 +integer :: lu_idx = 0 +integer :: ld_idx = 0 integer :: fsds_idx = 0 integer :: fsns_idx = 0 integer :: fsnt_idx = 0 integer :: flns_idx = 0 integer :: flnt_idx = 0 -integer :: cldfsnow_idx = 0 -integer :: cld_idx = 0 +integer :: cldfsnow_idx = 0 +integer :: cld_idx = 0 character(len=4) :: diag(0:N_DIAG) =(/' ','_d1 ','_d2 ','_d3 ','_d4 ','_d5 ','_d6 ','_d7 ','_d8 ','_d9 ','_d10'/) @@ -213,7 +213,7 @@ subroutine radiation_readnl(nlfile) if (iradlw < 0) iradlw = nint((-iradlw *3600._r8)/dtime) if (irad_always < 0) irad_always = nint((-irad_always*3600._r8)/dtime) - !----------------------------------------------------------------------- + !----------------------------------------------------------------------- ! Print runtime options to log. !----------------------------------------------------------------------- @@ -239,8 +239,8 @@ subroutine radiation_register use physics_buffer, only: pbuf_add_field, dtype_r8 use radiation_data, only: rad_data_register - call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate - call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate + call pbuf_add_field('QRS' , 'global',dtype_r8,(/pcols,pver/), qrs_idx) ! shortwave radiative heating rate + call pbuf_add_field('QRL' , 'global',dtype_r8,(/pcols,pver/), qrl_idx) ! longwave radiative heating rate call pbuf_add_field('FSDS' , 'global',dtype_r8,(/pcols/), fsds_idx) ! Surface solar downward flux call pbuf_add_field('FSNS' , 'global',dtype_r8,(/pcols/), fsns_idx) ! Surface net shortwave flux @@ -303,15 +303,15 @@ end function radiation_do !================================================================================================ real(r8) function radiation_nextsw_cday() - + ! Return calendar day of next sw radiation calculation ! Local variables integer :: nstep ! timestep counter - logical :: dosw ! true => do shosrtwave calc + logical :: dosw ! true => do shosrtwave calc integer :: offset ! offset for calendar day calculation integer :: dTime ! integer timestep size - real(r8):: calday ! calendar day of + real(r8):: calday ! calendar day of !----------------------------------------------------------------------- radiation_nextsw_cday = -1._r8 @@ -323,14 +323,14 @@ real(r8) function radiation_nextsw_cday() nstep = nstep + 1 offset = offset + dtime if (radiation_do('sw', nstep)) then - radiation_nextsw_cday = get_curr_calday(offset=offset) + radiation_nextsw_cday = get_curr_calday(offset=offset) dosw = .true. end if end do if(radiation_nextsw_cday == -1._r8) then call endrun('error in radiation_nextsw_cday') end if - + end function radiation_nextsw_cday !================================================================================================ @@ -368,7 +368,7 @@ subroutine radiation_init(pbuf2d) integer :: dtime !----------------------------------------------------------------------- - + call rad_solar_var_init() call rrtmg_state_init() call rad_data_init(pbuf2d) ! initialize output fields for offline driver @@ -408,12 +408,12 @@ subroutine radiation_init(pbuf2d) end if if (docosp) call cospsimulator_intr_init - + allocate(cosp_cnt(begchunk:endchunk)) if (is_first_restart_step()) then cosp_cnt(begchunk:endchunk) = cosp_cnt_init else - cosp_cnt(begchunk:endchunk) = 0 + cosp_cnt(begchunk:endchunk) = 0 end if call addfld('O3colAbove', horiz_only, 'A', 'DU', 'Column O3 above model top', sampling_seq='rad_lwsw') @@ -625,7 +625,7 @@ subroutine radiation_define_restart(file) end if end subroutine radiation_define_restart - + !=============================================================================== subroutine radiation_write_restart(file) @@ -644,7 +644,7 @@ subroutine radiation_write_restart(file) end if end subroutine radiation_write_restart - + !=============================================================================== subroutine radiation_read_restart(file) @@ -674,21 +674,21 @@ subroutine radiation_read_restart(file) end if end subroutine radiation_read_restart - + !=============================================================================== subroutine radiation_tend( & state, ptend, pbuf, cam_out, cam_in, net_flx, rd_out) - !----------------------------------------------------------------------- - ! + !----------------------------------------------------------------------- + ! ! Driver for radiation computation. - ! + ! ! Revision history: ! 2007-11-05 M. Iacono Install rrtmg_lw and sw as radiation model. ! 2007-12-27 M. Iacono Modify to use CAM cloud optical properties with rrtmg. !----------------------------------------------------------------------- - + use phys_grid, only: get_rlat_all_p, get_rlon_all_p use cam_control_mod, only: eccen, mvelpp, lambm0, obliqr use shr_orb_mod, only: shr_orb_decl, shr_orb_cosz @@ -718,7 +718,7 @@ subroutine radiation_tend( & ! Arguments type(physics_state), intent(in), target :: state type(physics_ptend), intent(out) :: ptend - + type(physics_buffer_desc), pointer :: pbuf(:) type(cam_out_t), intent(inout) :: cam_out type(cam_in_t), intent(in) :: cam_in @@ -731,7 +731,7 @@ subroutine radiation_tend( & type(rad_out_t), pointer :: rd ! allow rd_out to be optional by allocating a local object ! if the argument is not present logical :: write_output - + integer :: i, k integer :: lchnk, ncol logical :: dosw, dolw @@ -743,7 +743,7 @@ subroutine radiation_tend( & real(r8) :: clon(pcols) ! current longitudes(radians) real(r8) :: coszrs(pcols) ! Cosine solar zenith angle - ! Gathered indices of day and night columns + ! Gathered indices of day and night columns ! chunk_column_index = IdxDay(daylight_column_index) integer :: Nday ! Number of daylight columns integer :: Nnite ! Number of night columns @@ -754,8 +754,8 @@ subroutine radiation_tend( & real(r8), pointer :: cld(:,:) ! cloud fraction real(r8), pointer :: cldfsnow(:,:) ! cloud fraction of just "snow clouds- whatever they are" - real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate - real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate + real(r8), pointer :: qrs(:,:) ! shortwave radiative heating rate + real(r8), pointer :: qrl(:,:) ! longwave radiative heating rate real(r8), pointer :: fsds(:) ! Surface solar down flux real(r8), pointer :: fsns(:) ! Surface solar absorbed flux real(r8), pointer :: fsnt(:) ! Net column abs solar flux at model top @@ -809,7 +809,7 @@ subroutine radiation_tend( & real(r8) :: c_cld_tau_w_f(nswbands,pcols,pver) ! combined cloud forward scattered fraction * w * tau real(r8) :: c_cld_lw_abs (nlwbands,pcols,pver) ! combined cloud absorption optics depth (LW) - real(r8) :: sfac(1:nswbands) ! time varying scaling factors due to Solar Spectral Irrad at 1 A.U. per band + real(r8) :: sfac(1:nswbands,pcols) ! time varying scaling factors due to Solar Spectral Irrad at 1 A.U. per band integer :: icall ! index through climate/diagnostic radiation calls logical :: active_calls(0:N_DIAG) @@ -825,7 +825,7 @@ subroutine radiation_tend( & real(r8) :: fcns(pcols,pverp) ! net clear-sky shortwave flux real(r8) :: fnl(pcols,pverp) ! net longwave flux real(r8) :: fcnl(pcols,pverp) ! net clear-sky longwave flux - + ! for COSP real(r8) :: emis(pcols,pver) ! Cloud longwave emissivity real(r8) :: gb_snow_tau(pcols,pver) ! grid-box mean snow_tau @@ -934,8 +934,8 @@ subroutine radiation_tend( & else cldfprime(:ncol,:) = cld(:ncol,:) end if - - + + if (dosw) then if (oldcldoptics) then @@ -965,7 +965,7 @@ subroutine radiation_tend( & cld_tau_w(:,:ncol,:) = liq_tau_w(:,:ncol,:) + ice_tau_w(:,:ncol,:) cld_tau_w_g(:,:ncol,:) = liq_tau_w_g(:,:ncol,:) + ice_tau_w_g(:,:ncol,:) cld_tau_w_f(:,:ncol,:) = liq_tau_w_f(:,:ncol,:) + ice_tau_w_f(:,:ncol,:) - + if (cldfsnow_idx > 0) then ! add in snow call get_snow_optics_sw(state, pbuf, snow_tau, snow_tau_w, snow_tau_w_g, snow_tau_w_f) @@ -1081,7 +1081,7 @@ subroutine radiation_tend( & if (dosw) then - call get_variability(sfac) + call get_variability(lchnk,sfac) ! Get the active climate/diagnostic shortwave calculations call rad_cnst_get_call_list(active_calls) @@ -1096,7 +1096,7 @@ subroutine radiation_tend( & call aer_rad_props_sw(icall, state, pbuf, nnite, idxnite, & aer_tau, aer_tau_w, aer_tau_w_g, aer_tau_w_f) - + rd%cld_tau_cloudsim(:ncol,:) = cld_tau(rrtmg_sw_cloudsim_band,:ncol,:) rd%aer_tau550(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag) rd%aer_tau400(:ncol,:) = aer_tau(:ncol,:,idx_sw_diag+1) @@ -1133,7 +1133,7 @@ subroutine radiation_tend( & ! Output aerosol mmr call rad_cnst_out(0, state, pbuf) - + ! Longwave radiation computation if (dolw) then @@ -1149,7 +1149,7 @@ subroutine radiation_tend( & call rrtmg_state_update( state, pbuf, icall, r_state) call aer_rad_props_lw(icall, state, pbuf, aer_lw_abs) - + call rad_rrtmg_lw( & lchnk, ncol, num_rrtmg_levs, r_state, state%pmid, & aer_lw_abs, cldfprime, c_cld_lw_abs, qrl, rd%qrlc, & @@ -1390,7 +1390,7 @@ subroutine radiation_output_lw(lchnk, ncol, icall, rd, pbuf, cam_out, freqclr, f call outfld('FLUT'//diag(icall), rd%flut, pcols, lchnk) call outfld('FLUTC'//diag(icall), rd%flutc, pcols, lchnk) - + ftem(:ncol) = rd%flutc(:ncol) - rd%flut(:ncol) call outfld('LWCF'//diag(icall), ftem, pcols, lchnk) @@ -1411,7 +1411,7 @@ end subroutine radiation_output_lw subroutine calc_col_mean(state, mmr_pointer, mean_value) - ! Compute the column mean mass mixing ratio. + ! Compute the column mean mass mixing ratio. type(physics_state), intent(in) :: state real(r8), dimension(:,:), pointer :: mmr_pointer ! mass mixing ratio (lev) @@ -1440,4 +1440,3 @@ end subroutine calc_col_mean !=============================================================================== end module radiation - diff --git a/src/physics/rrtmg/radsw.F90 b/src/physics/rrtmg/radsw.F90 index df222557dd..5e0edd6443 100644 --- a/src/physics/rrtmg/radsw.F90 +++ b/src/physics/rrtmg/radsw.F90 @@ -1,7 +1,7 @@ module radsw -!----------------------------------------------------------------------- -! +!----------------------------------------------------------------------- +! ! Purpose: Solar radiation calculations. ! !----------------------------------------------------------------------- @@ -44,7 +44,7 @@ module radsw subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & E_pmid ,E_cld , & E_aer_tau,E_aer_tau_w,E_aer_tau_w_g,E_aer_tau_w_f, & - eccf ,E_coszrs ,solin ,sfac , & + eccf ,E_coszrs ,solin ,E_sfac , & E_asdir ,E_asdif ,E_aldir ,E_aldif , & qrs ,qrsc ,fsnt ,fsntc ,fsntoa,fsutoa, & fsntoac ,fsnirtoa ,fsnrtoac ,fsnrtoaq ,fsns , & @@ -57,33 +57,33 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & !----------------------------------------------------------------------- -! -! Purpose: +! +! Purpose: ! Solar radiation code -! -! Method: +! +! Method: ! mji/rrtmg ! RRTMG, two-stream, with McICA -! +! ! Divides solar spectrum into 14 intervals from 0.2-12.2 micro-meters. ! solar flux fractions specified for each interval. allows for ! seasonally and diurnally varying solar input. Includes molecular, -! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud, +! cloud, aerosol, and surface scattering, along with h2o,o3,co2,o2,cloud, ! and surface absorption. Computes delta-eddington reflections and -! transmissions assuming homogeneously mixed layers. Adds the layers -! assuming scattering between layers to be isotropic, and distinguishes +! transmissions assuming homogeneously mixed layers. Adds the layers +! assuming scattering between layers to be isotropic, and distinguishes ! direct solar beam from scattered radiation. -! +! ! Longitude loops are broken into 1 or 2 sections, so that only daylight ! (i.e. coszrs > 0) computations are done. -! +! ! Note that an extra layer above the model top layer is added. -! +! ! mks units are used. -! +! ! Special diagnostic calculation of the clear sky surface and total column ! absorbed flux is also done for cloud forcing diagnostics. -! +! !----------------------------------------------------------------------- use cmparray_mod, only: CmpDayNite, ExpDayNite @@ -91,8 +91,8 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & use mcica_subcol_gen_sw, only: mcica_subcol_sw use physconst, only: cpair use rrtmg_state, only: rrtmg_state_t - - ! Minimum cloud amount (as a fraction of the grid-box area) to + + ! Minimum cloud amount (as a fraction of the grid-box area) to ! distinguish from clear sky real(r8), parameter :: cldmin = 1.0e-80_r8 @@ -126,12 +126,12 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8), intent(in) :: E_aldir(pcols) ! 0.7-5.0 micro-meter srfc alb: direct rad real(r8), intent(in) :: E_asdif(pcols) ! 0.2-0.7 micro-meter srfc alb: diffuse rad real(r8), intent(in) :: E_aldif(pcols) ! 0.7-5.0 micro-meter srfc alb: diffuse rad - real(r8), intent(in) :: sfac(nbndsw) ! factor to account for solar variability in each band + real(r8), intent(in) :: E_sfac(nbndsw,pcols) ! factor to account for solar variability in each band real(r8), optional, intent(in) :: E_cld_tau (nbndsw, pcols, pver) ! cloud optical depth - real(r8), optional, intent(in) :: E_cld_tau_w (nbndsw, pcols, pver) ! cloud optical - real(r8), optional, intent(in) :: E_cld_tau_w_g(nbndsw, pcols, pver) ! cloud optical - real(r8), optional, intent(in) :: E_cld_tau_w_f(nbndsw, pcols, pver) ! cloud optical + real(r8), optional, intent(in) :: E_cld_tau_w (nbndsw, pcols, pver) ! cloud optical + real(r8), optional, intent(in) :: E_cld_tau_w_g(nbndsw, pcols, pver) ! cloud optical + real(r8), optional, intent(in) :: E_cld_tau_w_f(nbndsw, pcols, pver) ! cloud optical logical, optional, intent(in) :: old_convert ! Output arguments @@ -175,6 +175,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: rel(pcols,rrtmg_levs-1) ! Liquid effective drop size (microns) real(r8) :: rei(pcols,rrtmg_levs-1) ! Ice effective drop size (microns) + real(r8) :: sfac(nbndsw,pcols) ! factor to account for solar variability in each band real(r8) :: coszrs(pcols) ! Cosine solar zenith angle real(r8) :: asdir(pcols) ! 0.2-0.7 micro-meter srfc alb: direct rad real(r8) :: aldir(pcols) ! 0.7-5.0 micro-meter srfc alb: direct rad @@ -183,16 +184,16 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: h2ovmr(pcols,rrtmg_levs) ! h2o volume mixing ratio real(r8) :: o3vmr(pcols,rrtmg_levs) ! o3 volume mixing ratio - real(r8) :: co2vmr(pcols,rrtmg_levs) ! co2 volume mixing ratio - real(r8) :: ch4vmr(pcols,rrtmg_levs) ! ch4 volume mixing ratio - real(r8) :: o2vmr(pcols,rrtmg_levs) ! o2 volume mixing ratio - real(r8) :: n2ovmr(pcols,rrtmg_levs) ! n2o volume mixing ratio + real(r8) :: co2vmr(pcols,rrtmg_levs) ! co2 volume mixing ratio + real(r8) :: ch4vmr(pcols,rrtmg_levs) ! ch4 volume mixing ratio + real(r8) :: o2vmr(pcols,rrtmg_levs) ! o2 volume mixing ratio + real(r8) :: n2ovmr(pcols,rrtmg_levs) ! n2o volume mixing ratio real(r8) :: tsfc(pcols) ! surface temperature integer :: dyofyr ! Set to day of year for Earth/Sun distance calculation in ! rrtmg_sw, or pass in adjustment directly into adjes - real(r8) :: solvar(nbndsw) ! solar irradiance variability in each band + real(r8) :: solvar(nbndsw,pcols) ! solar irradiance variability in each band integer, parameter :: nsubcsw = ngptsw ! rrtmg_sw g-point (quadrature point) dimension integer :: permuteseed ! permute seed for sub-column generator @@ -209,7 +210,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: asm_aer_sw(pcols, rrtmg_levs-1, nbndsw) ! aer asymmetry parameter real(r8) :: cld_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud fraction - real(r8) :: rei_stosw(pcols, rrtmg_levs-1) ! stochastic ice particle size + real(r8) :: rei_stosw(pcols, rrtmg_levs-1) ! stochastic ice particle size real(r8) :: rel_stosw(pcols, rrtmg_levs-1) ! stochastic liquid particle size real(r8) :: cicewp_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud ice water path real(r8) :: cliqwp_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud liquid wter path @@ -219,7 +220,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & real(r8) :: fsfc_stosw(nsubcsw, pcols, rrtmg_levs-1) ! stochastic cloud forward scattering fraction (optional) real(r8), parameter :: dps = 1._r8/86400._r8 ! Inverse of seconds per day - + real(r8) :: swuflx(pcols,rrtmg_levs+1) ! Total sky shortwave upward flux (W/m2) real(r8) :: swdflx(pcols,rrtmg_levs+1) ! Total sky shortwave downward flux (W/m2) real(r8) :: swhr(pcols,rrtmg_levs) ! Total sky shortwave radiative heating rate (K/d) @@ -255,7 +256,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Aerosol radiative property arrays real(r8) :: tauxar(pcols,0:pver) ! aerosol extinction optical depth real(r8) :: wa(pcols,0:pver) ! aerosol single scattering albedo - real(r8) :: ga(pcols,0:pver) ! aerosol assymetry parameter + real(r8) :: ga(pcols,0:pver) ! aerosol asymmetry parameter real(r8) :: fa(pcols,0:pver) ! aerosol forward scattered fraction ! CRM @@ -304,7 +305,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & qrsc(1:ncol,1:pver) = 0.0_r8 fns(1:ncol,1:pverp) = 0.0_r8 fcns(1:ncol,1:pverp) = 0.0_r8 - if (single_column.and.scm_crm_mode) then + if (single_column.and.scm_crm_mode) then fus(1:ncol,1:pverp) = 0.0_r8 fds(1:ncol,1:pverp) = 0.0_r8 fusc(:ncol,:pverp) = 0.0_r8 @@ -331,6 +332,10 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & call CmpDayNite(r_state%o3vmr, o3vmr, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, rrtmg_levs) call CmpDayNite(r_state%co2vmr, co2vmr, Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, rrtmg_levs) + do i = 1,nbndsw + call CmpDayNite(E_sfac(i,:), sfac(i,:), Nday, IdxDay, Nnite, IdxNite, 1, pcols) + end do + call CmpDayNite(E_coszrs, coszrs, Nday, IdxDay, Nnite, IdxNite, 1, pcols) call CmpDayNite(E_asdir, asdir, Nday, IdxDay, Nnite, IdxNite, 1, pcols) call CmpDayNite(E_aldir, aldir, Nday, IdxDay, Nnite, IdxNite, 1, pcols) @@ -385,14 +390,14 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Define solar incident radiation do i = 1, Nday - solin(i) = sum(sfac(:)*solar_band_irrad(:)) * eccf * coszrs(i) + solin(i) = sum(sfac(:,i)*solar_band_irrad(:)) * eccf * coszrs(i) end do ! Calculate cloud optical properties here if using CAM method, or if using one of the - ! methods in RRTMG_SW, then pass in cloud physical properties and zero out cloud optical + ! methods in RRTMG_SW, then pass in cloud physical properties and zero out cloud optical ! properties here - ! Zero optional cloud optical property input arrays tauc_sw, ssac_sw, asmc_sw, + ! Zero optional cloud optical property input arrays tauc_sw, ssac_sw, asmc_sw, ! if inputting cloud physical properties to RRTMG_SW !tauc_sw(:,:,:) = 0.0_r8 !ssac_sw(:,:,:) = 1.0_r8 @@ -415,7 +420,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & fsfc_sw(ns,i,k) = 0._r8 asmc_sw(ns,i,k) = 0._r8 endif - + tauc_sw(ns,i,k)=E_cld_tau(ns,IdxDay(i),kk) if (tauc_sw(ns,i,k) > 0._r8) then ssac_sw(ns,i,k)=E_cld_tau_w(ns,IdxDay(i),kk)/tauc_sw(ns,i,k) @@ -441,7 +446,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & fsfc_sw(ns,i,k) = 0._r8 asmc_sw(ns,i,k) = 0._r8 endif - + tauc_sw(ns,i,k)=E_cld_tau(ns,IdxDay(i),kk) if (tauc_sw(ns,i,k) > 0._r8) then ssac_sw(ns,i,k)=max(E_cld_tau_w(ns,IdxDay(i),kk),1.e-80_r8)/max(tauc_sw(ns,i,k),1.e-80_r8) @@ -487,7 +492,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Call sub-column generator for McICA in radiation call t_startf('mcica_subcol_sw') - ! Set permute seed (must be offset between LW and SW by at least 140 to insure + ! Set permute seed (must be offset between LW and SW by at least 140 to insure ! effective randomization) permuteseed = 1 @@ -510,13 +515,13 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & tsfc(:ncol) = tlev(:ncol,rrtmg_levs+1) - solvar(1:nbndsw) = sfac(1:nbndsw) + solvar(1:nbndsw,1:ncol) = sfac(1:nbndsw,1:ncol) call rrtmg_sw(lchnk, Nday, rrtmg_levs, icld, & pmidmb, pintmb, tlay, tlev, tsfc, & h2ovmr, o3vmr, co2vmr, ch4vmr, o2vmr, n2ovmr, & asdir, asdif, aldir, aldif, & - coszrs, eccf, dyofyr, solvar, & + coszrs, eccf, dyofyr, solvar(1:nbndsw,1:Nday), & cld_stosw, tauc_stosw, ssac_stosw, asmc_stosw, fsfc_stosw, & cicewp_stosw, cliqwp_stosw, rei, rel, & tau_aer_sw, ssa_aer_sw, asm_aer_sw, & @@ -526,8 +531,8 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & ! Flux units are in W/m2 on output from rrtmg_sw and contain output for ! extra layer above model top with vertical indexing from bottom to top. ! - ! Heating units are in J/kg/s on output from rrtmg_sw and contain output - ! for extra layer above model top with vertical indexing from bottom to top. + ! Heating units are in J/kg/s on output from rrtmg_sw and contain output + ! for extra layer above model top with vertical indexing from bottom to top. ! ! Reverse vertical indexing to go from top to bottom for CAM output. @@ -545,7 +550,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & fsnt(1:Nday) = swdflx(1:Nday,rrtmg_levs) - swuflx(1:Nday,rrtmg_levs) fsntc(1:Nday) = swdflxc(1:Nday,rrtmg_levs) - swuflxc(1:Nday,rrtmg_levs) - ! Set the downwelling flux at the surface + ! Set the downwelling flux at the surface fsds(1:Nday) = swdflx(1:Nday,1) fsdsc(1:Nday) = swdflxc(1:Nday,1) @@ -622,7 +627,7 @@ subroutine rad_rrtmg_sw(lchnk,ncol ,rrtmg_levs ,r_state , & end if ! these outfld calls don't work for spmd only outfield in scm mode (nonspmd) - if (single_column .and. scm_crm_mode) then + if (single_column .and. scm_crm_mode) then ! Following outputs added for CRM call ExpDayNite(fus,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) call ExpDayNite(fds,Nday, IdxDay, Nnite, IdxNite, 1, pcols, 1, pverp) @@ -639,9 +644,9 @@ end subroutine rad_rrtmg_sw !------------------------------------------------------------------------------- subroutine radsw_init() -!----------------------------------------------------------------------- -! -! Purpose: +!----------------------------------------------------------------------- +! +! Purpose: ! Initialize various constants for radiation scheme. ! !----------------------------------------------------------------------- @@ -654,7 +659,7 @@ subroutine radsw_init() ! Initialize rrtmg_sw call rrtmg_sw_ini - + end subroutine radsw_init