subroutine sw_venus_corrk(ngrid,nlayer, & mu0,pplev,pplay,pt,tsurf,fract,dist_star, & dtsw,nfluxsurf,nfluxtop,netflux,firstcall) use mod_phys_lmdz_para, only : is_master use radinc_h, only: NBinfrared, NBvisible, L_NSPECTV, naerkind,& L_LEVELS, L_NGAUSS, L_NLEVRAD, L_NLAYRAD, L_REFVAR use radcommon_h, only: fzerov, gasv, & gweight, pfgasref, pgasmax, pgasmin, & pgasref, tgasmax, tgasmin, tgasref, scalep, & stellarf, dwnv, tauray use datafile_mod, only: datadir,banddir,corrkdir use ioipsl_getin_p_mod, only: getin_p use gases_h, only: ngasmx use optcv_mod, only: optcv use cpdet_phy_mod, only: cpdet implicit none #include "YOMCST.h" #include "clesphys.h" !================================================================== ! ! Purpose ! ------- ! Solve the radiative transfer using the correlated-k method for ! the gaseous absorption and the Toon et al. (1989) method for ! scatttering due to aerosols. ! ! Based on callcorrk (Generic GCM) ! adapted for the SW in the Venus GCM (S. Lebonnois, summer 2020) !================================================================== !----------------------------------------------------------------------- ! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid ! Layer #1 is the layer near the ground. ! Layer #nlayer is the layer at the top. !----------------------------------------------------------------------- ! INPUT INTEGER,INTENT(IN) :: ngrid ! Number of atmospheric columns. INTEGER,INTENT(IN) :: nlayer ! Number of atmospheric layers. REAL,INTENT(IN) :: mu0(ngrid) ! Cosine of sun incident angle. REAL,INTENT(IN) :: pplev(ngrid,nlayer+1) ! Inter-layer pressure (Pa). REAL,INTENT(IN) :: pplay(ngrid,nlayer) ! Mid-layer pressure (Pa). REAL,INTENT(IN) :: pt(ngrid,nlayer) ! Air temperature (K). REAL,INTENT(IN) :: tsurf(ngrid) ! Surface temperature (K). REAL,INTENT(IN) :: fract(ngrid) ! Fraction of day. REAL,INTENT(IN) :: dist_star ! Distance star-planet (AU). logical,intent(in) :: firstcall ! Signals first call to physics. ! OUTPUT REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation. REAL,INTENT(OUT) :: nfluxsurf(ngrid) ! Net SW flux absorbed at surface (W/m2). REAL,INTENT(OUT) :: nfluxtop(ngrid) ! Net incident top of atmosphere SW flux (W/m2). REAL,INTENT(OUT) :: netflux(ngrid,nlayer) ! net SW flux (W/m2). ! interesting variables ! albedo could be an input map... For future adaptation REAL :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015 ! potential outputs REAL :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2). REAL :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2). REAL :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2) REAL :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1). ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h. REAL :: QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind) REAL :: omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) REAL :: gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind) REAL :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght. REAL :: tau_col(ngrid) ! Diagnostic from aeropacity. REAL,ALLOCATABLE,SAVE :: reffrad(:,:,:) ! aerosol effective radius (m) REAL,ALLOCATABLE,SAVE :: nueffrad(:,:,:) ! aerosol effective variance !$OMP THREADPRIVATE(reffrad,nueffrad) !----------------------------------------------------------------------- ! Declaration of the variables required by correlated-k subroutines ! Numbered from top to bottom (unlike in the GCM) !----------------------------------------------------------------------- REAL*8 tmid(L_LEVELS),pmid(L_LEVELS) REAL*8 tlevrad(L_LEVELS),plevrad(L_LEVELS) ! Optical values for the optci/cv subroutines REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV) ! NB: Arrays below are "save" to avoid reallocating them at every call ! not because their content needs be reused from call to the next REAL*8,allocatable,save :: dtauv(:,:,:) REAL*8,allocatable,save :: cosbv(:,:,:) REAL*8,allocatable,save :: wbarv(:,:,:) !$OMP THREADPRIVATE(dtauv,cosbv,wbarv) REAL*8,allocatable,save :: tauv(:,:,:) REAL*8,allocatable,save :: taucumv(:,:,:) !$OMP THREADPRIVATE(tauv,taucumv) REAL*8 tauaero(L_LEVELS,naerkind) REAL*8 nfluxtopv,fluxtopvdn REAL*8 nfluxoutv_nu(L_NSPECTV) ! Outgoing band-resolved VI flux at TOA (W/m2). REAL*8 fmnetv(L_NLAYRAD) REAL*8 fluxupv(L_NLAYRAD) REAL*8 fluxdnv(L_NLAYRAD) REAL*8 acosz REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. INTEGER ig,l,k,nw,iaer real*8,allocatable,save :: taugsurf(:,:) !$OMP THREADPRIVATE(taugsurf) real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol). ! Local aerosol optical properties for each column on RADIATIVE grid. real*8,save,allocatable :: QXVAER(:,:,:) ! Extinction coeff (QVISsQREF*QREFvis) real*8,save,allocatable :: QSVAER(:,:,:) real*8,save,allocatable :: GVAER(:,:,:) !$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER) real, dimension(:,:,:), save, allocatable :: QREFvis3d !$OMP THREADPRIVATE(QREFvis3d) ! Miscellaneous : real*8 temp,temp1,temp2,pweight character(len=10) :: tmp1 character(len=10) :: tmp2 character(len=100) :: message character(len=10),parameter :: subname="sw_corrk" ! For fixed water vapour profiles. integer i_var real RH real psat,qsat logical OLRz real*8 NFLUXGNDV_nu(L_NSPECTV) ! Included by RW for runaway greenhouse 1D study. real vtmp(nlayer) REAL*8 muvarrad(L_LEVELS) !=============================================================== ! I.a Initialization on first call !=============================================================== if(firstcall) then call iniaerosol() allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind)) allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind)) allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind)) ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind)) ! Effective radius and variance of the aerosols allocate(reffrad(ngrid,nlayer,naerkind)) allocate(nueffrad(ngrid,nlayer,naerkind)) !-------------------------------------------------- ! Set up correlated k !-------------------------------------------------- print*, "callcorrk: Correlated-k data base folder:",trim(datadir) print*, "corrkdir = ",corrkdir write( tmp1, '(i3)' ) NBinfrared write( tmp2, '(i3)' ) NBvisible banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)) banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)) call su_gases ! reading of gases.def bypassed in this, hardcoded. cf su_gases.F90 call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) call setspv ! Basic visible properties. call sugas_corrk ! Set up gaseous absorption properties. call suaer_corrk ! Set up aerosol optical properties. ! now that L_NGAUSS has been initialized (by sugas_corrk) ! allocate related arrays ALLOCATE(dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)) ALLOCATE(cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)) ALLOCATE(wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)) ALLOCATE(tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)) ALLOCATE(taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)) ALLOCATE(taugsurf(L_NSPECTV,L_NGAUSS-1)) OSR_nu(:,:) = 0. ! Surface albedo in the solar range ! example of data: Cutler et al 2020 ! reflectance of basalts in UV and visible varies from a few to 10% ! it also depends on mineralogy ! in the absence of further data, I take 5% as a mean value... Sensitivity could be tested... albedo(:,:) = 0.05 end if ! of if (firstcall) !======================================================================= ! I.b Initialization on every call !======================================================================= qxvaer(:,:,:)=0.0 qsvaer(:,:,:)=0.0 gvaer(:,:,:) =0.0 ! How much light do we get ? do nw=1,L_NSPECTV stel(nw)=stellarf(nw)/(dist_star**2) end do ! Get 3D aerosol optical properties. call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & QVISsQREF3d,omegaVIS3d,gVIS3d,QREFvis3d) ! Get aerosol optical depths. call aeropacity(ngrid,nlayer,pplay,pplev,pt,aerosol, & reffrad,nueffrad,QREFvis3d,tau_col) !----------------------------------------------------------------------- !----------------------------------------------------------------------- do ig=1,ngrid ! Starting Big Loop over every GCM column !----------------------------------------------------------------------- !----------------------------------------------------------------------- !======================================================================= ! II. Transformation of the GCM variables !======================================================================= !----------------------------------------------------------------------- ! Aerosol optical properties Qext, Qscat and g. ! The transformation in the vertical is the same as for temperature. !----------------------------------------------------------------------- do iaer=1,naerkind ! Shortwave. do nw=1,L_NSPECTV do l=1,nlayer temp1=QVISsQREF3d(ig,nlayer+1-l,nw,iaer) & *QREFvis3d(ig,nlayer+1-l,iaer) temp2=QVISsQREF3d(ig,max(nlayer-l,1),nw,iaer) & *QREFvis3d(ig,max(nlayer-l,1),iaer) qxvaer(2*l,nw,iaer) = temp1 qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=temp1*omegavis3d(ig,nlayer+1-l,nw,iaer) temp2=temp2*omegavis3d(ig,max(nlayer-l,1),nw,iaer) qsvaer(2*l,nw,iaer) = temp1 qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=gvis3d(ig,nlayer+1-l,nw,iaer) temp2=gvis3d(ig,max(nlayer-l,1),nw,iaer) gvaer(2*l,nw,iaer) = temp1 gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 end do ! nlayer qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer) qxvaer(2*nlayer+1,nw,iaer)=0. qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer) qsvaer(2*nlayer+1,nw,iaer)=0. gvaer(1,nw,iaer)=gvaer(2,nw,iaer) gvaer(2*nlayer+1,nw,iaer)=0. end do ! L_NSPECTV end do ! naerkind ! Test / Correct for freaky s. s. albedo values. do iaer=1,naerkind do k=1,L_LEVELS do nw=1,L_NSPECTV if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then message='Serious problems with qsvaer values' call abort_physic(subname,message,1) endif if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer) endif end do end do ! L_LEVELS end do ! naerkind !----------------------------------------------------------------------- ! Aerosol optical depths !----------------------------------------------------------------------- do iaer=1,naerkind do k=0,nlayer-1 pweight=(pplay(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1))/ & (pplev(ig,L_NLAYRAD-k)-pplev(ig,L_NLAYRAD-k+1)) ! As 'aerosol' is at reference (visible) wavelenght we scale it as ! it will be multplied by qxi/v in optci/v temp=aerosol(ig,L_NLAYRAD-k,iaer)/QREFvis3d(ig,L_NLAYRAD-k,iaer) tauaero(2*k+2,iaer)=max(temp*pweight,0.d0) tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0) end do ! boundary conditions tauaero(1,iaer) = tauaero(2,iaer) !tauaero(1,iaer) = 0. !JL18 at time of testing, the two above conditions gave the same results bit for bit. end do ! naerkind ! Albedo DO nw=1,L_NSPECTV ! Short Wave loop. albv(nw)=albedo(ig,nw) ENDDO acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude. !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !!! Note by JL13 : In the following, some indices were changed in the interpolations, !!! so that the model results are less dependent on the number of layers ! !!! !!! --- The older versions are commented with the comment !JL13index --- !!! !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! !----------------------------------------------------------------------- ! Variable species... Not used for Venus !----------------------------------------------------------------------- do k=1,L_LEVELS qvar(k) = 1.0D-7 end do DO l=1,nlayer muvarrad(2*l) = RMD muvarrad(2*l+1) = RMD END DO muvarrad(1) = muvarrad(2) muvarrad(2*nlayer+1)=RMD ! Keep values inside limits for which we have radiative transfer coefficients !!! if(L_REFVAR.gt.1)then ! (there was a bug here) message='no variable species for Venus yet' call abort_physic(subname,message,1) endif !----------------------------------------------------------------------- ! Pressure and temperature !----------------------------------------------------------------------- DO l=1,nlayer plevrad(2*l) = pplay(ig,nlayer+1-l)/scalep plevrad(2*l+1) = pplev(ig,nlayer+1-l)/scalep tlevrad(2*l) = pt(ig,nlayer+1-l) tlevrad(2*l+1) = (pt(ig,nlayer+1-l)+pt(ig,max(nlayer-l,1)))/2 END DO plevrad(1) = 0. ! plevrad(2) = 0. !! JL18 enabling this line puts the radiative top at p=0 which was the idea before, but does not seem to perform best after all. tlevrad(1) = tlevrad(2) tlevrad(2*nlayer+1)=tsurf(ig) pmid(1) = pplay(ig,nlayer)/scalep pmid(2) = pmid(1) tmid(1) = tlevrad(2) tmid(2) = tmid(1) DO l=1,L_NLAYRAD-1 tmid(2*l+1) = tlevrad(2*l+1) tmid(2*l+2) = tlevrad(2*l+1) pmid(2*l+1) = plevrad(2*l+1) pmid(2*l+2) = plevrad(2*l+1) END DO pmid(L_LEVELS) = plevrad(L_LEVELS) tmid(L_LEVELS) = tlevrad(L_LEVELS) !!Alternative interpolation: ! pmid(3) = pmid(1) ! pmid(4) = pmid(1) ! tmid(3) = tmid(1) ! tmid(4) = tmid(1) ! DO l=2,L_NLAYRAD-1 ! tmid(2*l+1) = tlevrad(2*l) ! tmid(2*l+2) = tlevrad(2*l) ! pmid(2*l+1) = plevrad(2*l) ! pmid(2*l+2) = plevrad(2*l) ! END DO ! pmid(L_LEVELS) = plevrad(L_LEVELS-1) ! tmid(L_LEVELS) = tlevrad(L_LEVELS-1) ! Test for out-of-bounds pressure. if(plevrad(3).lt.pgasmin)then ! print*,'Minimum pressure is outside the radiative' ! print*,'transfer kmatrix bounds, exiting.' ! message="Minimum pressure outside of kmatrix bounds" ! call abort_physic(subname,message,1) elseif(plevrad(L_LEVELS).gt.pgasmax)then ! print*,'Maximum pressure is outside the radiative' ! print*,'transfer kmatrix bounds, exiting.' ! message="Maximum pressure outside of kmatrix bounds" ! call abort_physic(subname,message,1) endif ! Test for out-of-bounds temperature. do k=1,L_LEVELS if(tlevrad(k).lt.tgasmin)then print*,'Minimum temperature is outside the radiative' print*,'transfer kmatrix bounds' print*,"k=",k," tlevrad(k)=",tlevrad(k) print*,"tgasmin=",tgasmin ! if (strictboundcorrk) then ! message="Minimum temperature outside of kmatrix bounds" ! call abort_physic(subname,message,1) ! else print*,'***********************************************' print*,'we allow model to continue with tlevradtgasmax' print*,' ... we assume we know what you are doing ... ' print*,' ... but do not let this happen too often ... ' print*,'***********************************************' !tlevrad(k)=tgasmax ! Used in the source function ! ! endif endif enddo do k=1,L_NLAYRAD+1 if(tmid(k).lt.tgasmin)then print*,'Minimum temperature is outside the radiative' print*,'transfer kmatrix bounds, exiting.' print*,"k=",k," tmid(k)=",tmid(k) print*,"tgasmin=",tgasmin ! if (strictboundcorrk) then ! message="Minimum temperature outside of kmatrix bounds" ! call abort_physic(subname,message,1) ! else print*,'***********************************************' print*,'we allow model to continue but with tmid=tgasmin' print*,' ... we assume we know what you are doing ... ' print*,' ... but do not let this happen too often ... ' print*,'***********************************************' tmid(k)=tgasmin ! endif elseif(tmid(k).gt.tgasmax)then print*,'Maximum temperature is outside the radiative' print*,'transfer kmatrix bounds, exiting.' print*,"k=",k," tmid(k)=",tmid(k) print*,"tgasmax=",tgasmax ! if (strictboundcorrk) then ! message="Maximum temperature outside of kmatrix bounds" ! call abort_physic(subname,message,1) ! else print*,'***********************************************' print*,'we allow model to continue but with tmid=tgasmax' print*,' ... we assume we know what you are doing ... ' print*,' ... but do not let this happen too often ... ' print*,'***********************************************' tmid(k)=tgasmax ! endif endif enddo !======================================================================= ! III. Calling the main radiative transfer subroutines !======================================================================= !----------------------------------------------------------------------- ! Short Wave Part !----------------------------------------------------------------------- if(fract(ig) .ge. 1.0e-4) then ! Only during daylight. do nw=1,L_NSPECTV stel_fract(nw)= stel(nw) * fract(ig) end do call optcv(dtauv,tauv,taucumv,plevrad, & qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, & tmid,pmid,taugsurf,qvar,muvarrad) call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, & acosz,stel_fract, & nfluxtopv,fluxtopvdn,nfluxoutv_nu,nfluxgndv_nu, & fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) else ! During the night, fluxes = 0. nfluxtopv = 0.0d0 fluxtopvdn = 0.0d0 nfluxoutv_nu(:) = 0.0d0 nfluxgndv_nu(:) = 0.0d0 do l=1,L_NLAYRAD fmnetv(l)=0.0d0 fluxupv(l)=0.0d0 fluxdnv(l)=0.0d0 end do end if !----------------------------------------------------------------------- ! Transformation of the correlated-k code outputs ! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw) ! Net incident flux tops, positive downward nfluxtop(ig) = -1.*nfluxtopv ! Net incident flux at surface, positive downward nfluxsurf(ig) = -1.*fmnetv(L_NLAYRAD) ! Flux incident at the top of the atmosphere fluxtop_dn(ig)= fluxtopvdn fluxabs_sw(ig) = real(-nfluxtopv) fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD)) if(fluxtop_dn(ig).lt.0.0)then print*,'Achtung! fluxtop_dn has lost the plot!' print*,'fluxtop_dn=',fluxtop_dn(ig) print*,'acosz=',acosz print*,'aerosol=',aerosol(ig,:,:) print*,'temp= ',pt(ig,:) print*,'pplay= ',pplay(ig,:) message="Achtung! fluxtop_dn has lost the plot!" call abort_physic(subname,message,1) endif ! Net solar flux, positive downward do l=1,L_NLAYRAD netflux(ig,L_NLAYRAD+1-l)= -1.*fmnetv(l) enddo netflux(ig,L_NLAYRAD+1) = -1.*nfluxtopv ! Finally, the heating rates DO l=2,L_NLAYRAD dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) & *RG/(cpdet(tmid(l))*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) END DO ! These are values at top of atmosphere dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) & *RG/(cpdet(tmid(1))*scalep*(plevrad(3)-plevrad(2))) !----------------------------------------------------------------------- !----------------------------------------------------------------------- end do ! End of big loop over every GCM column. !----------------------------------------------------------------------- !----------------------------------------------------------------------- end subroutine sw_venus_corrk