From 67c2e26ed271aa2bdbe32548f443f6a916464f95 Mon Sep 17 00:00:00 2001 From: Dustin Swales Date: Wed, 24 Apr 2019 10:53:16 -0600 Subject: [PATCH] Working./gmtb_scm twpice_control_RRTMGP_cloud --- physics/rrtmgp_lw_cloud_optics.F90 | 4 +- physics/rrtmgp_lw_main.F90 | 3 +- physics/rrtmgp_sw_cloud_optics.F90 | 191 ++++++++++++++--------------- physics/rrtmgp_sw_main.F90 | 73 ++++++----- 4 files changed, 141 insertions(+), 130 deletions(-) diff --git a/physics/rrtmgp_lw_cloud_optics.F90 b/physics/rrtmgp_lw_cloud_optics.F90 index d21a9221d..7b8c5f6f5 100644 --- a/physics/rrtmgp_lw_cloud_optics.F90 +++ b/physics/rrtmgp_lw_cloud_optics.F90 @@ -267,7 +267,7 @@ module mo_rrtmgp_lw_cloud_optics 1.202445e-02, 1.142838e-02, 1.086257e-02, 1.032445e-02, 9.811791e-03, & !3 9.322587e-03, 8.855053e-03, 8.407591e-03, 7.978763e-03, 7.567273e-03, & !3 7.171949e-03, 6.791728e-03, 6.425642e-03, 6.072809e-03, 5.732424e-03, & !3 - 5.403748e-03, 5.086103e-03, 4.778865e-03, & !3 + 5.403748e-03, 5.086103e-03, 4.778865e-03, & !3 1.804566e-01, 1.168987e-01, 8.680442e-02, 6.910060e-02, 5.738174e-02, & !4 4.902332e-02, 4.274585e-02, 3.784923e-02, 3.391734e-02, 3.068690e-02, & !4 2.798301e-02, 2.568480e-02, 2.370600e-02, 2.198337e-02, 2.046940e-02, & !4 @@ -321,7 +321,7 @@ module mo_rrtmgp_lw_cloud_optics 1.077117e-02, 1.035369e-02, 9.962643e-03, 9.595509e-03, 9.250088e-03, & !9 8.924447e-03, 8.616876e-03, 8.325862e-03, 8.050057e-03, 7.788258e-03, & !9 7.539388e-03, 7.302478e-03, 7.076656e-03, 6.861134e-03, 6.655197e-03, & !9 - 6.458197e-03, 6.269543e-03, 6.088697e-03, & !9 + 6.458197e-03, 6.269543e-03, 6.088697e-03, & !9 1.593628e-01, 1.014552e-01, 7.458955e-02, 5.903571e-02, 4.887582e-02, & !10 4.171159e-02, 3.638480e-02, 3.226692e-02, 2.898717e-02, 2.631256e-02, & !10 2.408925e-02, 2.221156e-02, 2.060448e-02, 1.921325e-02, 1.799699e-02, & !10 diff --git a/physics/rrtmgp_lw_main.F90 b/physics/rrtmgp_lw_main.F90 index c2837d294..09a092510 100644 --- a/physics/rrtmgp_lw_main.F90 +++ b/physics/rrtmgp_lw_main.F90 @@ -1058,6 +1058,7 @@ subroutine rrtmgp_lw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr do iCol=1,nCol do iGpt=1,nGptsLW iBand = kdist_lw_clr%convert_gpt2band(iGpt) + optical_props_clr%tau(iCol,1:nlay,iGpt) = optical_props_clr%tau(iCol,1:nlay,iGpt) * secdiff(iBand,iCol) optical_props_aer%tau(iCol,1:nlay,iGpt) = tau_aer(iCol,1:nlay,iBand) * (1. - ssa_aer(iCol,1:nlay,iBand)) * secdiff(iBand,iCol) enddo enddo @@ -1190,7 +1191,7 @@ subroutine rrtmgp_lw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr sum(fluxClrSky%flux_dn(1,iLay:iLay+1))/2.,sum(fluxAllSky%flux_up(1,iLay:iLay+1))/2.,& sum(fluxAllSky%flux_dn(1,iLay:iLay+1))/2. write(60,*) optical_props_clr%tau(1,iLay,:) - write(61,*) tau_cld(:,1,iLay)*secdiff(:,1) + write(61,'(16f12.3)') tau_cld(:,1,iLay)!*secdiff(:,1) enddo ! ####################################################################################### ! Copy fluxes from RRTGMP types into model radiation types. diff --git a/physics/rrtmgp_sw_cloud_optics.F90 b/physics/rrtmgp_sw_cloud_optics.F90 index cc353af66..10d4e8e4b 100644 --- a/physics/rrtmgp_sw_cloud_optics.F90 +++ b/physics/rrtmgp_sw_cloud_optics.F90 @@ -2070,25 +2070,35 @@ subroutine rrtmgp_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cl ! Local variables integer :: iCol, iLay, iBand, index, ia, istr - real(kind_phys) :: tau_rain, tau_snow, factor, fint, cld_ref_iceTemp + real(kind_phys) :: tau_rain, tau_snow, factor, fint, cld_ref_iceTemp,asyw,ssaw,za1,za2 + real(kind_phys), dimension(nBandsSW) :: ssa_rain, ssa_snow, asy_rain, asy_snow, & tau_liq, ssa_liq, asy_liq, tau_ice, ssa_ice, asy_ice, forwliq, asycoliq, & - forwice, extcoice, asycoice, ssacoice, fdelta, extcoliq, ssacoliq, & - tau_iceORG, tau_liqORG + forwice, extcoice, asycoice, ssacoice, fdelta, extcoliq, ssacoliq ! Initialize tau_cld(:,:,:) = 0._kind_phys - ssa_cld(:,:,:) = 0._kind_phys + ssa_cld(:,:,:) = 1._kind_phys asy_cld(:,:,:) = 0._kind_phys ! Compute cloud radiative properties for cloud. - print*,'RRTMGP_CLOUD_OPTICS_CONFIG: ' - print*,'iswcliq: ',iswcliq - print*,'iswcice: ',iswcice if (iswcliq > 0) then do iCol=1,ncol do iLay=1,nlay - if (cld_frac(iCol,iLay) .gt. 0.) then + ! Initialize + tau_liq(:) = 0._kind_phys + tau_ice(:) = 0._kind_phys + tau_rain = 0._kind_phys + tau_snow = 0._kind_phys + ssa_liq(:) = 0._kind_phys + ssa_ice(:) = 0._kind_phys + ssa_rain(:) = 0._kind_phys + ssa_snow(:) = 0._kind_phys + asy_liq(:) = 0._kind_phys + asy_ice(:) = 0._kind_phys + asy_rain(:) = 0._kind_phys + asy_snow(:) = 0._kind_phys + if (cld_frac(iCol,iLay) .gt. 0._kind_phys) then ! ########################################################################### ! Rain clouds ! ########################################################################### @@ -2121,44 +2131,36 @@ subroutine rrtmgp_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cl ! Liquid clouds ! ########################################################################### if (cld_lwp(iCol,iLay) .gt. 0) then + ! Find index in coefficient LUT for corresponding partice size. factor = cld_ref_liq(iCol,iLay) - 1.5 index = max( 1, min( 57, int( factor ) )) fint = factor - float(index) - if ( iswcliq == 1 ) then - do iBand=1,nBandsSW - extcoliq(iBand) = extliq1(index,iBand) + fint*(extliq1(index+1,iBand)-extliq1(index,iBand)) - ssacoliq(iBand) = ssaliq1(index,iBand) + fint*(ssaliq1(index+1,iBand)-ssaliq1(index,iBand)) - asycoliq(iBand) = asyliq1(index,iBand) + fint*(asyliq1(index+1,iBand)-asyliq1(index,iBand)) - if (fint .lt. 0._kind_phys .and. ssacoliq(iBand) .gt. 1._kind_phys) then - ssacoliq(iBand) = ssaliq1(index,iBand) - endif - forwliq(iBand) = asycoliq(iBand)*asycoliq(iBand) - tau_liqORG(iBand) = cld_lwp(iCol,iLay) * extcoliq(iBand) - tau_liq(iBand) = (1._kind_phys - forwliq(iBand) * ssacoliq(iBand)) * tau_liqORG(iBand) - ssa_liq(iBand) = ssacoliq(iBand) * (1._kind_phys - forwliq(iBand) ) / & - (1._kind_phys - forwliq(iBand)*ssacoliq(iBand)) - asy_liq(iBand) = (asycoliq(iBand) - forwliq(iBand)) / (1._kind_phys - forwliq(iBand)) - enddo - elseif ( iswcliq == 2 ) then ! use updated coeffs - do iBand=1,nBandsSW - extcoliq(iBand) = extliq2(index,iBand) + fint*(extliq2(index+1,iBand)-extliq2(index,iBand)) - ssacoliq(iBand) = ssaliq2(index,iBand) + fint*(ssaliq2(index+1,iBand)-ssaliq2(index,iBand)) - asycoliq(iBand) = asyliq2(index,iBand) + fint*(asyliq2(index+1,iBand)-asyliq2(index,iBand)) - if (fint .lt. 0._kind_phys .and. ssacoliq(iBand) .gt. 1._kind_phys) then - ssacoliq(iBand) = ssaliq2(index,iBand) - endif - forwliq(iBand) = asycoliq(iBand)*asycoliq(iBand) - tau_liqORG(iBand) = cld_lwp(iCol,iLay) * extcoliq(iBand) - tau_liq(iBand) = (1._kind_phys - forwliq(iBand) * ssacoliq(iBand)) * tau_liqORG(iBand) - ssa_liq(iBand) = ssacoliq(iBand) * (1._kind_phys - forwliq(iBand) ) / & - (1._kind_phys - forwliq(iBand)*ssacoliq(iBand)) - asy_liq(iBand) = (asycoliq(iBand) - forwliq(iBand)) / (1._kind_phys - forwliq(iBand)) - enddo - endif - else ! Cloudy, but no liquid condensate present - tau_liq(iBand) = 0._kind_phys - ssa_liq(iBand) = 0._kind_phys - asy_liq(iBand) = 0._kind_phys + + ! Extract coefficents for all bands and compute radiative properties + do iBand=1,nBandsSW + ! Interpolate coefficients + if ( iswcliq == 1 ) then + extcoliq(iBand) = max(0._kind_phys, extliq1(index,iBand) + & + fint*(extliq1(index+1,iBand)-extliq1(index,iBand))) + ssacoliq(iBand) = max(0._kind_phys, min(1._kind_phys, ssaliq1(index,iBand) + & + fint*(ssaliq1(index+1,iBand)-ssaliq1(index,iBand)))) + asycoliq(iBand) = max(0._kind_phys, min(1._kind_phys, asyliq1(index,iBand) + & + fint*(asyliq1(index+1,iBand)-asyliq1(index,iBand)))) + elseif ( iswcliq == 2 ) then ! use updated coeffs + extcoliq(iBand) = max(0._kind_phys, extliq2(index,iBand) + & + fint*(extliq2(index+1,iBand)-extliq2(index,iBand))) + ssacoliq(iBand) = max(0._kind_phys, min(1._kind_phys, ssaliq2(index,iBand) + & + fint*(ssaliq2(index+1,iBand)-ssaliq2(index,iBand)))) + asycoliq(iBand) = max(0._kind_phys, min(1._kind_phys, asyliq2(index,iBand) + & + fint*(asyliq2(index+1,iBand)-asyliq2(index,iBand)))) + endif + if (fint .lt. 0._kind_phys .and. ssacoliq(iBand) .gt. 1._kind_phys) then + ssacoliq(iBand) = ssaliq1(index,iBand) + endif + tau_liq(iBand) = cld_lwp(iCol,iLay) * extcoliq(iBand) + ssa_liq(iBand) = tau_liq(iBand) * ssacoliq(iBand) + asy_liq(iBand) = ssa_liq(iBand) * asycoliq(iBand) + enddo endif ! IF cloudy with liquid condensate ! ########################################################################### @@ -2186,9 +2188,12 @@ subroutine rrtmgp_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cl index = max( 1, min( 42, int( factor ) )) fint = factor - float(index) do iBand = 1,nBandsSW - extcoice(iBand) = extice2(index,iBand) + fint*(extice2(index+1,iBand)-extice2(index,iBand)) - ssacoice(iBand) = ssaice2(index,iBand) + fint*(ssaice2(index+1,iBand)-ssaice2(index,iBand)) - asycoice(iBand) = asyice2(index,iBand) + fint*(asyice2(index+1,iBand)-asyice2(index,iBand)) + extcoice(iBand) = extice2(index,iBand) + & + fint*(extice2(index+1,iBand)-extice2(index,iBand)) + ssacoice(iBand) = ssaice2(index,iBand) + & + fint*(ssaice2(index+1,iBand)-ssaice2(index,iBand)) + asycoice(iBand) = asyice2(index,iBand) + & + fint*(asyice2(index+1,iBand)-asyice2(index,iBand)) tau_ice(iBand) = cld_iwp(iCol,iLay) * extcoice(iBand) ssa_ice(iBand) = tau_ice(iBand) * ssacoice(iBand) asy_ice(iBand) = ssa_ice(iBand) * asycoice(iBand) @@ -2198,67 +2203,56 @@ subroutine rrtmgp_sw_cloud_optics(ncol, nlay, nBandsSW, cld_lwp, cld_ref_liq, cl ! (generalized effective size from 5 to 140 microns). ! https://doi.org/10.1175/1520-0442(1996)009<2058:AAPOTS>2.0.CO;2 elseif ( iswcice == 3 ) then + cld_ref_iceTemp = max( 5.0, min( 140.0, 1.0315*cld_ref_ice(iCol,iLay) )) ! Determine indices for table interpolation. - factor = (cld_ref_ice(iCol,iLay) - 2._kind_phys) / 3._kind_phys - index = int( factor ) + factor = (cld_ref_iceTemp - 2._kind_phys) / 3._kind_phys + index = max( 1, min( 45, int( factor ) )) fint = factor - float(index) do iBand = 1,nBandsSW - ! Interpolate coefficient tables to appropriate ice-particle size and band. - extcoice(iBand) = extice3(index,iBand) + fint*(extice3(index+1,iBand)-extice3(index,iBand)) ! eq (3.9a) - ssacoice(iBand) = ssaice3(index,iBand) + fint*(ssaice3(index+1,iBand)-ssaice3(index,iBand)) ! eq (3.9b) - asycoice(iBand) = asyice3(index,iBand) + fint*(asyice3(index+1,iBand)-asyice3(index,iBand)) ! eq (3.9c) - fdelta(iBand) = fdlice3(index,iBand) + fint*(fdlice3(index+1,iBand)-fdlice3(index,iBand)) ! eq (3.9d) - forwice(iBand) = fdelta(iBand) + 0.5_kind_phys/ssacoice(iBand) ! eq (3.8) + ! Interpolate coefficient tables to appropriate ice-particle size. + extcoice(iBand) = max(0._kind_phys, extice3(index,iBand) + & + fint*(extice3(index+1,iBand)-extice3(index,iBand))) ! eq (3.9a) + ssacoice(iBand) = max(0._kind_phys, min(1._kind_phys, ssaice3(index,iBand) + & + fint*(ssaice3(index+1,iBand)-ssaice3(index,iBand)))) ! eq (3.9b) + asycoice(iBand) = max(0._kind_phys, min(1._kind_phys, asyice3(index,iBand) + & + fint*(asyice3(index+1,iBand)-asyice3(index,iBand)))) ! eq (3.9c) + fdelta(iBand) = fdlice3(index,iBand) + & + fint*(fdlice3(index+1,iBand)-fdlice3(index,iBand)) ! eq (3.9d) + forwice(iBand) = fdelta(iBand) + 0.5_kind_phys / ssacoice(iBand) if (forwice(iBand) .gt. asycoice(iBand)) forwice(iBand) = asycoice(iBand) - - ! Calculatice radiative properties - tau_iceORG(iBand) = cld_iwp(iCol,iLay) * extcoice(iBand) - tau_ice(iBand) = (1._kind_phys - forwice(iBand)*ssacoice(iBand)) * tau_iceORG(iBand) ! eq (A.2a) - ssa_ice(iBand) = ssacoice(iBand) * (1._kind_phys - forwice(iBand) ) / & ! eq (A.2b) - (1._kind_phys - forwice(iBand)*ssacoice(iBand)) - asy_ice(iBand) = (asycoice(iBand) - forwice(iBand)) / & ! eq (A.2c) - (1._kind_phys - forwice(iBand)) + tau_ice(iBand) = cld_iwp(iCol,iLay) * extcoice(iBand) + ssa_ice(iBand) = tau_ice(iBand) * ssacoice(iBand) + asy_ice(iBand) = ssa_ice(iBand) * asycoice(iBand) enddo endif - else ! Cloudy, but no ice condensate present - tau_ice(iBand) = 0._kind_phys - ssa_ice(iBand) = 0._kind_phys - asy_ice(iBand) = 0._kind_phys endif ! IF cloudy column with ice condensate - else ! No clouds - tau_liq(:) = 0._kind_phys - tau_ice(:) = 0._kind_phys - tau_rain = 0._kind_phys - tau_snow = 0._kind_phys - ssa_liq(:) = 0._kind_phys - ssa_ice(:) = 0._kind_phys - ssa_rain(:) = 0._kind_phys - ssa_snow(:) = 0._kind_phys - asy_liq(:) = 0._kind_phys - asy_ice(:) = 0._kind_phys - asy_rain(:) = 0._kind_phys - asy_snow(:) = 0._kind_phys endif ! IF cloudy column ! ########################################################################### - ! Compute total cloud optical depth + ! Compute total cloud radiative properties (tau, omega, and g) ! ########################################################################### - do iBand = 1,nBandsSW - tau_cld(iBand,iCol,iLay) = tau_liq(iBand) + tau_ice(iBand) + tau_rain + tau_snow - ssa_cld(iBand,iCol,iLay) = (tau_liq(iBand) * ssa_liq(iBand) + & - tau_ice(iBand) * ssa_ice(iBand) + & - tau_rain * ssa_rain(iBand) + & - tau_snow * ssa_snow(iBand)) / & - (tau_liq(iBand) + tau_ice(iBand) + tau_rain + tau_snow) - asy_cld(iBand,iCol,iLay) = (tau_liq(iBand) * ssa_liq(iBand) * asy_liq(iBand) + & - tau_ice(iBand) * ssa_ice(iBand) * asy_ice(iBand) + & - tau_rain * ssa_rain(iBand) * asy_rain(iBand) + & - tau_snow * ssa_snow(iBand) * asy_snow(iBand)) / & - (tau_liq(iBand) * ssa_liq(iBand) + & - tau_ice(iBand) * ssa_ice(iBand) + & - tau_rain * ssa_rain(iBand) + & - tau_snow * ssa_snow(iBand)) - enddo ! Loop over G-points + if (cld_frac(iCol,iLay) .gt. 0._kind_phys) then + do iBand = 1,nBandsSW + ! Sum up radiative properties by type. + tau_cld(iBand,iCol,iLay) = tau_liq(iBand) + tau_ice(iBand) + tau_rain + tau_snow + ssa_cld(iBand,iCol,iLay) = ssa_liq(iBand) + ssa_ice(iBand) + ssa_rain(iBand) + ssa_snow(iBand) + asy_cld(iBand,iCol,iLay) = asy_liq(iBand) + asy_ice(iBand) + asy_rain(iBand) + asy_snow(iBand) + ! Delta-scale + asyw = asy_cld(iband,iCol,iLay)/max(0._kind_phys, ssa_cld(iBand,iCol,iLay)) + ssaw = min(1._kind_phys-0.000001, ssa_cld(iBand,iCol,iLay)/tau_cld(iBand,iCol,iLay)) + za1 = asyw * asyw + za2 = ssaw * za1 + tau_cld(iBand,iCol,iLay) = (1._kind_phys - za2) * tau_cld(iband,iCol,iLay) + ssa_cld(iBand,iCol,iLay) = (ssaw - za2) / (1._kind_phys - za2) + asy_cld(iBand,iCol,iLay) = (asyw - za2/ssaw)/(1-za2/ssaw) + !asy_cld(iBand,iCol,iLay) = (tau_liq(iBand)*ssa_liq(iBand)*(asycoliq(iBand)-asycoliq(iBand)**2)/& + ! (1 - asycoliq(iBand)**2) + & + ! tau_ice(iBand)*ssa_ice(iBand)*(asycoice(iBand)-forwice(iBand))/& + ! (1 - forwice(iBand)))/& + ! (tau_liq(iBand)*ssa_liq(iBand) + tau_ice(iBand)*ssa_ice(iBand)) + enddo ! Loop over SW bands + endif ! END sum cloudy properties + ! enddo ! Loop over layers enddo ! Loop over columns endif @@ -2266,8 +2260,9 @@ end subroutine rrtmgp_sw_cloud_optics ! ####################################################################################### ! SUBROUTINE mcica_subcol_sw - ! ####################################################################################### - subroutine mcica_subcol_sw(ncol, nlay, ngpts, cld_frac, icseed, dzlyr, de_lgth, cld_frac_mcica) + ! ###################################################################################### + subroutine mcica_subcol_sw(ncol, nlay, ngpts, cld_frac, icseed, dzlyr, de_lgth, & + cld_frac_mcica) ! Inputs integer,intent(in) :: & ncol, & ! Number of horizontal gridpoints diff --git a/physics/rrtmgp_sw_main.F90 b/physics/rrtmgp_sw_main.F90 index 1554d5699..53f1ae716 100644 --- a/physics/rrtmgp_sw_main.F90 +++ b/physics/rrtmgp_sw_main.F90 @@ -154,6 +154,9 @@ subroutine rrtmgp_sw_init(Model,mpicomm, mpirank, mpiroot, errmsg, errflg) integer,dimension(:),allocatable :: temp1,temp2,temp3,temp4,temp_log_array1,& temp_log_array2, temp_log_array3, temp_log_array4 + open(69,file='rrtmgp_sw_aux_dump.txt',status='unknown') + open(70,file='rrtmgp_sw_aux_taucld.txt',status='unknown') + ! Initialize errmsg = '' errflg = 0 @@ -644,9 +647,10 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr ! RTE+RRTMGP classes type(ty_optical_props_2str) :: & - optical_props_clr, & ! Optical properties for gaseous atmosphere - optical_props_aer, & ! Optical properties for aerosols - optical_props_cldy ! Optical properties for clouds + optical_props_clr, & ! Optical properties for gaseous atmosphere + optical_props_aer, & ! Optical properties for aerosols + optical_props_cldy ! Optical properties for clouds + type(ty_fluxes_broadband) :: & fluxAllSky, & ! All-sky flux (W/m2) fluxClrSky ! Clear-sky flux (W/m2) @@ -654,12 +658,12 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr ! Local variables integer :: iCol, iBand, iGpt, iDay, iLay integer,dimension(ncol) :: ipseed - real(kind_phys) :: solAdjFac,cfrac + real(kind_phys) :: solAdjFac, cfrac, asyw, ssaw, za1, za2 logical :: top_at_1=.false. real(kind_phys), dimension(ncol) :: clrfracSFC, cldfracSFC real(kind_phys), dimension(ncol,nlay) :: vmr_o3, vmr_h2o, thetaTendClrSky, & thetaTendAllSky,coldry,tem0 - real(kind_phys), dimension(nday,nlay) :: cldfrac2 + real(kind_phys), dimension(nday,nlay) :: cldfrac2, cld_lwp2 real(kind_phys), dimension(nday,nlay+1),target :: & flux_up_allSky, flux_up_clrSky, flux_dn_allSky, flux_dn_clrSky, p_lev2 real(kind_phys), dimension(nday,nGptsSW) :: toa_flux @@ -821,7 +825,7 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr optical_props_clr, & toa_flux)) - ! 1c) Add contribution from aerosols. Also, apply diffusivity angle + ! 1c) Add contribution from aerosols. print*,'Clear-Sky(SW): Increment Aerosol' do iDay=1,nDay do iGpt=1,nGptsSW @@ -831,7 +835,7 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr optical_props_aer%g(iDay,1:nlay,iGpt) = asy_aer(idxday(iDay),1:nlay,iBand) enddo enddo - call check_error_msg(optical_props_aer%increment(optical_props_clr)) + !call check_error_msg(optical_props_aer%increment(optical_props_clr)) ! 1d) Compute the clear-sky broadband fluxes print*,'Clear-Sky(SW): Fluxes' @@ -853,9 +857,11 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr ! #################################################################################### ! 2a) Compute in-cloud optics print*,'All-Sky(SW): Optics ' - tau_cld(1:nBandsSW,1:nDay,1:nLay) = 0._kind_phys - asy_cld(1:nBandsSW,1:nDay,1:nLay) = 0._kind_phys - ssa_cld(1:nBandsSW,1:nDay,1:nLay) = 0._kind_phys + tau_cld(:,:,:) = 0._kind_phys + asy_cld(:,:,:) = 0._kind_phys + ssa_cld(:,:,:) = 0._kind_phys + + cldfrac2 = merge(cldfrac(idxday,:), 0., cldfrac(idxday,:) .gt. 0._kind_phys) if (any(cldfrac(idxday,:) .gt. 0)) then ! Cloud-optical properties by type provided. Compute optical-depth, single- ! scattering albedo, and asymmetry parameter @@ -866,7 +872,7 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr cld_rwp(idxday,1:nLay), cld_ref_rain(idxday,1:nLay), & cld_swp(idxday,1:nLay), cld_ref_snow(idxday,1:nLay), & cldfrac(idxday,1:nLay), & - tau_cld(:,:,:), ssa_cld(:,:,:), asy_cld(:,:,:)) + tau_cld, ssa_cld, asy_cld) else ! Cloud-optical depth, single scattering albedo, and asymmetry parameter provided. do iDay=1,nDay @@ -877,7 +883,7 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr asy_cld(:,iDay,iLay) = cld_asy(idxday(iDay),iLay) else tau_cld(:,iDay,iLay) = 0. - ssa_cld(:,iDay,iLay) = 0. + ssa_cld(:,iDay,iLay) = 1. asy_cld(:,iDay,iLay) = 0. endif end do @@ -885,11 +891,21 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr endif endif + write(69,'(a20)') "RRTMGP TAUs" + write(70,*) "#" + do iDay=1,nDay + do iLay=1,nlay + write(69,'(a5,i2,4f12.3)') '',iLay,p_lay(idxday(iDay),iLay),sum(optical_props_clr%tau(iDay,iLay,:)) + write(70,'(16f12.3)') tau_cld(:,1,iLay) + write(70,'(16f12.3)') ssa_cld(:,1,iLay) + write(70,'(16f12.3)') asy_cld(:,1,iLay) + enddo + enddo + ! 2b) Call McICA to sample clouds. if (isubcsw .gt. 0) then print*,'All-Sky(SW): McICA' allocate(tau(nDay,nLay,nGptsSW),asy(nDay,nLay,nGptsSW),ssa(nDay,nLay,nGptsSW)) - cldfrac2 = merge(cldfrac(idxday,:), 0., cldfrac(idxday,:) .gt. 0._kind_phys) call mcica_subcol_sw(nday, nlay, nGptsSW, cldfrac2, ipseed(idxday), dzlyr(idxday,:), & de_lgth(idxday), cldfracMCICA) @@ -899,13 +915,13 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr do iGpt=1,nGptsSW iBand = kdist_sw_clr%convert_gpt2band(iGpt) if (cldfracMCICA(iGpt,iDay,iLay) .gt. 0.) then - tau(iDay,iLay,iGpt) = tau_cld(iband,iDay,iLay) - asy(iDay,iLay,iGpt) = asy_cld(iband,iDay,iLay) - ssa(iDay,iLay,iGpt) = ssa_cld(iband,iDay,iLay) + optical_props_cldy%tau(iDay,iLay,iGpt) = tau_cld(iband,iDay,iLay) + optical_props_cldy%ssa(iDay,iLay,iGpt) = ssa_cld(iBand,iDay,iLay) + optical_props_cldy%g(iDay,iLay,iGpt) = asy_cld(iBand,iDay,iLay) else - tau(iDay,iLay,iGpt) = 0. - ssa(iDay,iLay,iGpt) = 0. - asy(iDay,iLay,iGpt) = 0. + optical_props_cldy%tau(iDay,iLay,iGpt) = 0._kind_phys + optical_props_cldy%ssa(iDay,iLay,iGpt) = 1._kind_phys + optical_props_cldy%g(iDay,iLay,iGpt) = 0._kind_phys endif enddo enddo @@ -919,21 +935,18 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr cfrac = cldfrac(idxday(iDay),iLay)/cldfracSFC(idxday(iDay)) if (cfrac .gt. 0.) then do iBand=1,nBandsSW - tau(iDay,iLay,iBand) = tau_cld(iBand,iDay,iLay) - asy(iDay,iLay,iBand) = asy_cld(iBand,iDay,iLay) - ssa(iDay,iLay,iBand) = ssa_cld(iBand,iDay,iLay) + optical_props_cldy%tau(iDay,iLay,iBand) = tau_cld(iBand,iDay,iLay) + optical_props_cldy%g(iDay,iLay,iBand) = asy_cld(iBand,iDay,iLay) + optical_props_cldy%ssa(iDay,iLay,iBand) = ssa_cld(iBand,iDay,iLay) enddo else - tau(iDay,iLay,:) = 0. - asy(iDay,iLay,:) = 0. - ssa(iDay,iLay,:) = 0. + optical_props_cldy%tau(iDay,iLay,:) = 0._kind_phys + optical_props_cldy%g(iDay,iLay,:) = 1._kind_phys + optical_props_cldy%ssa(iDay,iLay,:) = 0._kind_phys endif enddo enddo endif - optical_props_cldy%tau = tau - optical_props_cldy%ssa = ssa - optical_props_cldy%g = asy ! 2c) Add cloud contribution from the gaseous (clear-sky) atmosphere. print*,'All-Sky(SW): Increment' @@ -966,7 +979,7 @@ subroutine rrtmgp_sw_run(p_lay, p_lev, t_lay, t_lev, q_lay, o3_lay, vmr_co2, vmr sfcflx(idxday)%upfx0 = fluxClrSky%flux_up(:,1) sfcflx(idxday)%dnfxc = fluxAllSky%flux_dn(:,1) sfcflx(idxday)%dnfx0 = fluxClrSky%flux_dn(:,1) - cldtau(idxday,:) = 0. + cldtau(idxday,:) = tau(:,:,10) hswc(idxday,:) = thetaTendAllSky ! Optional output @@ -989,6 +1002,8 @@ end subroutine rrtmgp_sw_run ! ######################################################################################### ! ######################################################################################### subroutine rrtmgp_sw_finalize() + close(69) + close(70) end subroutine rrtmgp_sw_finalize ! #########################################################################################