subroutine callcorrk(icount,ngrid,nlayer,pq,nq,qsurf, & albedo,emis,mu0,pplev,pplay,pt, & tsurf,fract,dist_star,igout,aerosol,cpp3D, & dtlw,dtsw,fluxsurf_lw, & fluxsurf_sw,fluxtop_lw,fluxtop_sw,fluxtop_dn, & reffrad,tau_col,ptime,pday,firstcall,lastcall,zzlay) use radinc_h use radcommon_h use ioipsl_getincom use radii_mod use aerosol_mod use datafile_mod, only: datadir implicit none !================================================================== ! ! 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. ! ! Authors ! ------- ! Emmanuel 01/2001, Forget 09/2001 ! Robin Wordsworth (2009) ! Modif Pluton Tanguy Bertrand 2017 !================================================================== #include "dimphys.h" #include "comcstfi.h" #include "callkeys.h" #include "tracer.h" !----------------------------------------------------------------------- ! Declaration of the arguments (INPUT - OUTPUT) on the LMD GCM grid ! Layer #1 is the layer near the ground. ! Layer #nlayermx is the layer at the top. ! INPUT INTEGER icount INTEGER ngrid,nlayer INTEGER igout REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol opacity tau REAL albedo(ngrid) ! SW albedo REAL emis(ngrid) ! LW emissivity REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer REAL zzlay(ngrid,nlayermx) ! altitude at the middle of the layers REAL pplev(ngrid,nlayermx+1) ! pres. level at GCM layer boundaries REAL pt(ngrid,nlayermx) ! air temperature (K) REAL tsurf(ngrid) ! surface temperature (K) REAL dist_star,mu0(ngrid) ! distance star-planet (AU) REAL fract(ngrid) ! fraction of day ! Globally varying aerosol optical properties on GCM grid ! Not needed everywhere so not in radcommon_h REAL :: QVISsQREF3d(ngridmx,nlayermx,L_NSPECTV,naerkind) REAL :: omegaVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind) REAL :: gVIS3d(ngridmx,nlayermx,L_NSPECTV,naerkind) REAL :: QIRsQREF3d(ngridmx,nlayermx,L_NSPECTI,naerkind) REAL :: omegaIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind) REAL :: gIR3d(ngridmx,nlayermx,L_NSPECTI,naerkind) REAL :: QREFvis3d(ngridmx,nlayermx,naerkind) REAL :: QREFir3d(ngridmx,nlayermx,naerkind) ! REAL :: omegaREFvis3d(ngridmx,nlayermx,naerkind) ! REAL :: omegaREFir3d(ngridmx,nlayermx,naerkind) ! not sure of the point of these... REAL reffrad(ngrid,nlayer,naerkind) REAL nueffrad(ngrid,nlayer,naerkind) REAL profhaze(ngrid,nlayer) ! TB17 fixed profile of haze mmr ! OUTPUT REAL dtsw(ngridmx,nlayermx) ! heating rate (K/s) due to SW REAL dtlw(ngridmx,nlayermx) ! heating rate (K/s) due to LW REAL fluxsurf_lw(ngridmx) ! incident LW flux to surf (W/m2) REAL fluxtop_lw(ngridmx) ! outgoing LW flux to space (W/m2) REAL fluxsurf_sw(ngridmx) ! incident SW flux to surf (W/m2) REAL fluxtop_sw(ngridmx) ! outgoing LW flux to space (W/m2) REAL fluxtop_dn(ngridmx) ! incident top of atmosphere SW flux (W/m2) !----------------------------------------------------------------------- ! 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) REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS) REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS) REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS) REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS) REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS) REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS) REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS) REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS) REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS) REAL*8 tauaero(L_LEVELS+1,naerkind) REAL*8 nfluxtopv,nfluxtopi,nfluxtop real*8 NFLUXTOPV_nu(L_NSPECTV) real*8 NFLUXTOPI_nu(L_NSPECTI) real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD) real*8 fmneti_nu(L_NLAYRAD,L_NSPECTI) ! real*8 fmnetv_nu(L_NLAYRAD,L_NSPECTV) ! REAL*8 fluxupv(L_NLAYRAD),fluxupi(L_NLAYRAD) REAL*8 fluxdnv(L_NLAYRAD),fluxdni(L_NLAYRAD) REAL*8 albi,albv,acosz INTEGER ig,l,k,nw,iaer,irad real fluxtoplanet real*8 taugsurf(L_NSPECTV,L_NGAUSS-1) real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1) real*8 qvar(L_LEVELS) ! mixing ratio of variable component REAL pq(ngridmx,nlayer,nq) REAL qsurf(ngridmx,nqmx) ! tracer on surface (e.g. kg.m-2) integer nq ! Local aerosol optical properties for each column on RADIATIVE grid real*8 QXVAER(L_LEVELS+1,L_NSPECTV,naerkind) real*8 QSVAER(L_LEVELS+1,L_NSPECTV,naerkind) real*8 GVAER(L_LEVELS+1,L_NSPECTV,naerkind) real*8 QXIAER(L_LEVELS+1,L_NSPECTI,naerkind) real*8 QSIAER(L_LEVELS+1,L_NSPECTI,naerkind) real*8 GIAER(L_LEVELS+1,L_NSPECTI,naerkind) save qxvaer, qsvaer, gvaer save qxiaer, qsiaer, giaer save QREFvis3d, QREFir3d REAL tau_col(ngrid) ! diagnostic from aeropacity ! Misc. logical firstcall, lastcall, nantest real*8 tempv(L_NSPECTV) real*8 tempi(L_NSPECTI) real*8 temp,temp1,temp2,pweight character(len=10) :: tmp1 character(len=10) :: tmp2 ! for fixed vapour profiles real RH real*8 pq_temp(nlayer) real ptemp, Ttemp, qsat ! for OLR spec integer OLRcount save OLRcount integer OLRcount2 save OLRcount2 character(len=2) :: tempOLR character(len=30) :: filenomOLR real ptime, pday REAL epsi_ch4 SAVE epsi_ch4 logical diagrad_OLRz,diagrad_OLR,diagrad_surf,diagrad_rates real OLR_nu(ngrid,L_NSPECTI) ! Allow variations in cp with location REAL cpp3D(ngridmx,nlayermx) ! specific heat capacity at const. pressure ! NLTE factor for CH4 real eps_nlte_sw23(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw real eps_nlte_sw33(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw real eps_nlte_lw(ngridmx,nlayermx) ! CH4 NLTE efficiency factor for zdtsw REAL dtlw_nu(nlayermx,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands REAL dtsw_nu(nlayermx,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands REAL dpp ! intermediate REAL dtlw_co(ngridmx, nlayermx) ! cooling rate (K/s) due to CO (diagnostic) REAL dtlw_hcn_c2h2(ngridmx, nlayermx) ! cooling rate (K/s) due to C2H2/HCN (diagnostic) !!read altitudes and vmrch4 integer Nfine,ifine parameter(Nfine=701) character(len=100) :: file_path real,save :: levdat(Nfine),vmrdat(Nfine) real :: vmrch4(ngridmx,nlayermx) ! vmr ch4 from vmrch4_proffix !======================================================================= ! Initialization on first call CALL zerophys((L_LEVELS+1)*L_NSPECTV*naerkind,qxvaer) CALL zerophys((L_LEVELS+1)*L_NSPECTV*naerkind,qsvaer) CALL zerophys((L_LEVELS+1)*L_NSPECTV*naerkind,gvaer) CALL zerophys((L_LEVELS+1)*L_NSPECTI*naerkind,qxiaer) CALL zerophys((L_LEVELS+1)*L_NSPECTI*naerkind,qsiaer) CALL zerophys((L_LEVELS+1)*L_NSPECTI*naerkind,giaer) if(firstcall) then print*, "callcorrk: Correlated-k data folder:",trim(datadir) call getin("corrkdir",corrkdir) print*, "corrkdir = ",corrkdir write( tmp1, '(i3)' ) L_NSPECTI write( tmp2, '(i3)' ) L_NSPECTV banddir=trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2)) banddir=trim(adjustl(corrkdir))//'/'//trim(adjustl(banddir)) print*,'starting sugas' call sugas_corrk ! set up gaseous absorption properties print*,'starting setspi' call setspi ! basic infrared properties print*,'starting setspv' call setspv ! basic visible properties ! Radiative Hazes if (aerohaze) then print*,'aerohaze: starting suaer_corrk' call suaer_corrk ! set up aerosol optical properties print*,'ending suaer_corrk' !-------------------------------------------------- ! Effective radius and variance of the aerosols !-------------------------------------------------- do iaer=1,naerkind if ((iaer.eq.iaero_haze)) then call haze_reffrad(ngrid,nlayer,reffrad(1,1,iaer), & nueffrad(1,1,iaer)) endif end do !iaer=1,naerkind. if (haze_radproffix) then print*, 'haze_radproffix=T : fixed profile for haze rad' else print*,'reffrad haze:',reffrad(1,1,iaero_haze) print*,'nueff haze',nueffrad(1,1,iaero_haze) endif endif ! radiative haze Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed epsi_ch4=mmol(igcm_ch4_gas)/mmol(igcm_n2) ! If fixed profile of CH4 gas IF (vmrch4_proffix) then file_path=trim(datadir)//'/gas_prop/vmr_ch4.txt' open(115,file=file_path,form='formatted') do ifine=1,Nfine read(115,*) levdat(ifine), vmrdat(ifine) enddo close(115) ENDIF end if ! firstcall !======================================================================= ! L_NSPECTV is the number of Visual(or Solar) spectral intervals ! how much light we get fluxtoplanet=0 DO nw=1,L_NSPECTV stel(nw)=stellarf(nw)/(dist_star**2) !flux fluxtoplanet=fluxtoplanet + stel(nw) END DO !----------------------------------------------------------------------- ! Get 3D aerosol optical properties. ! ici on selectionne les proprietes opt correspondant a reffrad if (aerohaze) then !-------------------------------------------------- ! Effective radius and variance of the aerosols if profil non ! uniform. Called only if aerohaze=true. !-------------------------------------------------- if (haze_radproffix) then call haze_reffrad_fix(ngrid,nlayer,zzlay, & reffrad,nueffrad) endif call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & QVISsQREF3d,omegaVIS3d,gVIS3d, & QIRsQREF3d,omegaIR3d,gIR3d, & QREFvis3d,QREFir3d) ! Get aerosol optical depths. call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, & reffrad,QREFvis3d,QREFir3d, & tau_col) endif !----------------------------------------------------------------------- ! Prepare CH4 mixing ratio for radiative transfer IF (methane) then vmrch4(:,:)=0. if (ch4fix) then if (vmrch4_proffix) then !! Interpolate on the model vertical grid do ig=1,ngridmx ! CALL interp_line(levdat,vmrdat,Nfine, ! & zzlay(ig,:)/1000.,vmrch4(ig,:),nlayer) enddo else vmrch4(:,:)=vmrch4fix endif else vmrch4(:,:)=pq(:,:,igcm_ch4_gas)*100.* & mmol(igcm_n2)/mmol(igcm_ch4_gas) endif ENDIF ! Prepare NON LTE correction in Pluto atmosphere IF (nlte) then CALL nlte_ch4(ngrid,nlayer,nq,pplay,pplev,pt,vmrch4, & eps_nlte_sw23,eps_nlte_sw33,eps_nlte_lw) ENDIF c Net atmospheric radiative cooling rate from C2H2 (K.s-1): c ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ dtlw_hcn_c2h2=0. if (cooling) then CALL cooling_hcn_c2h2(ngrid,nlayer,pplay, & pt,dtlw_hcn_c2h2) endif !----------------------------------------------------------------------- !======================================================================= !----------------------------------------------------------------------- ! starting big loop over every GCM column do ig=1,ngridmx !======================================================================= ! Transformation of the GCM variables !----------------------------------------------------------------------- ! Aerosol optical properties Qext, Qscat and g on each band ! The transformation in the vertical is the same as for temperature if (aerohaze) then ! shortwave do iaer=1,naerkind DO nw=1,L_NSPECTV do l=1,nlayermx temp1=QVISsQREF3d(ig,nlayermx+1-l,nw,iaer) $ *QREFvis3d(ig,nlayermx+1-l,iaer) temp2=QVISsQREF3d(ig,max(nlayermx-l,1),nw,iaer) $ *QREFvis3d(ig,max(nlayermx-l,1),iaer) qxvaer(2*l,nw,iaer) = temp1 qxvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=temp1*omegavis3d(ig,nlayermx+1-l,nw,iaer) temp2=temp2*omegavis3d(ig,max(nlayermx-l,1),nw,iaer) qsvaer(2*l,nw,iaer) = temp1 qsvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=gvis3d(ig,nlayermx+1-l,nw,iaer) temp2=gvis3d(ig,max(nlayermx-l,1),nw,iaer) gvaer(2*l,nw,iaer) = temp1 gvaer(2*l+1,nw,iaer)=(temp1+temp2)/2 end do qxvaer(1,nw,iaer)=qxvaer(2,nw,iaer) qxvaer(2*nlayermx+1,nw,iaer)=0. qsvaer(1,nw,iaer)=qsvaer(2,nw,iaer) qsvaer(2*nlayermx+1,nw,iaer)=0. gvaer(1,nw,iaer)=gvaer(2,nw,iaer) gvaer(2*nlayermx+1,nw,iaer)=0. end do ! longwave DO nw=1,L_NSPECTI do l=1,nlayermx temp1=QIRsQREF3d(ig,nlayermx+1-l,nw,iaer) $ *QREFir3d(ig,nlayermx+1-l,iaer) temp2=QIRsQREF3d(ig,max(nlayermx-l,1),nw,iaer) $ *QREFir3d(ig,max(nlayermx-l,1),iaer) qxiaer(2*l,nw,iaer) = temp1 qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=temp1*omegair3d(ig,nlayermx+1-l,nw,iaer) temp2=temp2*omegair3d(ig,max(nlayermx-l,1),nw,iaer) qsiaer(2*l,nw,iaer) = temp1 qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=gir3d(ig,nlayermx+1-l,nw,iaer) temp2=gir3d(ig,max(nlayermx-l,1),nw,iaer) giaer(2*l,nw,iaer) = temp1 giaer(2*l+1,nw,iaer)=(temp1+temp2)/2 end do qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer) qxiaer(2*nlayermx+1,nw,iaer)=0. qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer) qsiaer(2*nlayermx+1,nw,iaer)=0. giaer(1,nw,iaer)=giaer(2,nw,iaer) giaer(2*nlayermx+1,nw,iaer)=0. end do 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 print*,'Serious problems with qsvaer values' print*,'in callcorrk' call abort endif if(qsvaer(k,nw,iaer).gt.qxvaer(k,nw,iaer))then qsvaer(k,nw,iaer)=qxvaer(k,nw,iaer) endif end do do nw=1,L_NSPECTI if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then print*,'Serious problems with qsiaer values' print*,'in callcorrk' call abort endif if(qsiaer(k,nw,iaer).gt.qxiaer(k,nw,iaer))then qsiaer(k,nw,iaer)=qxiaer(k,nw,iaer) endif end do end do ! levels end do ! naerkind endif ! aerohaze !----------------------------------------------------------------------- ! Aerosol optical depths IF (aerohaze) THEN do iaer=1,naerkind ! heritage generic 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)) if (QREFvis3d(ig,L_NLAYRAD-k,iaer).ne.0) then temp=aerosol(ig,L_NLAYRAD-k,iaer)/ $ QREFvis3d(ig,L_NLAYRAD-k,iaer) else print*, 'stop corrk',k,QREFvis3d(ig,L_NLAYRAD-k,iaer) stop end if tauaero(2*k+2,iaer)=max(temp*pweight,0.d0) tauaero(2*k+3,iaer)=max(temp-tauaero(2*k+2,iaer),0.d0) ! tauaero en L_LEVELS soit deux fois plus que nlayer end do ! generic New boundary conditions tauaero(1,iaer) = tauaero(2,iaer) !tauaero(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer) !tauaero(1,iaer) = 0. !tauaero(L_LEVELS+1,iaer) = 0. end do ! naerkind ELSE tauaero(:,:)=0 ENDIF !----------------------------------------------------------------------- ! Albedo and emissivity albi=1-emis(ig) ! longwave albv=albedo(ig) ! shortwave acosz=mu0(ig) ! cosine of sun incident angle !----------------------------------------------------------------------- ! Methane vapour c qvar = mixing ratio c L_LEVELS (51) different de GCM levels (25) . L_LEVELS = 2*(llm-1)+3=2*(ngridmx-1)+3 c L_REFVAR The number of different mixing ratio values in c datagcm/composition.in for the k-coefficients. qvar(:)=0. IF (methane) then do l=1,nlayer qvar(2*l) = vmrch4(ig,nlayer+1-l)/100.* & mmol(igcm_ch4_gas)/mmol(igcm_n2) qvar(2*l+1) = ((vmrch4(ig,nlayer+1-l)+vmrch4(ig, & max(nlayer-l,1)))/2.)/100.* & mmol(igcm_ch4_gas)/mmol(igcm_n2) end do qvar(1)=qvar(2) ! Keep values inside limits for which we have radiative transfer coefficients if(L_REFVAR.gt.1)then ! there was a bug here! do k=1,L_LEVELS if(qvar(k).lt.wrefvar(1))then qvar(k)=wrefvar(1)+1.0e-8 elseif(qvar(k).gt.wrefvar(L_REFVAR))then qvar(k)=wrefvar(L_REFVAR)-1.0e-8 endif end do endif ! IMPORTANT: Now convert from kg/kg to mol/mol do k=1,L_LEVELS qvar(k) = qvar(k)/(epsi_ch4+qvar(k)*(1.-epsi_ch4)) end do ENDIF ! methane !----------------------------------------------------------------------- ! Pressure and temperature ! generic updated: 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. !! Trick to have correct calculations of fluxes ! in gflux(i/v).F, but the pmid levels are not impacted by this change. tlevrad(1) = tlevrad(2) tlevrad(2*nlayer+1)=tsurf(ig) pmid(1) = max(pgasmin,0.0001*plevrad(3)) pmid(2) = pmid(1) tmid(1) = tlevrad(2) tmid(2) = tmid(1) ! INI ! 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 !! following lines changed in 03/2015 to solve upper atmosphere bug ! plevrad(1) = 0. ! plevrad(2) = max(pgasmin,0.0001*plevrad(3)) ! ! tlevrad(1) = tlevrad(2) ! tlevrad(2*nlayermx+1)=tsurf(ig) ! ! tmid(1) = tlevrad(2) ! tmid(2) = tlevrad(2) ! ! pmid(1) = plevrad(2) ! pmid(2) = plevrad(2) 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 ! end of changes pmid(L_LEVELS) = plevrad(L_LEVELS) tmid(L_LEVELS) = tlevrad(L_LEVELS) !TB if ((PMID(2).le.1.e-5).and.(ig.eq.1)) then if ((TMID(2).le.30.).and.(ig.eq.1)) then write(*,*) 'Caution! Pres/temp of upper levels lower than' write(*,*) 'ref pressure/temp: kcoef fixed for upper levels' !! cf tpindex.F endif endif ! test for out-of-bounds pressure if(plevrad(3).lt.pgasmin)then print*,'Minimum pressure is outside the radiative' print*,'transfer kmatrix bounds, exiting.' ! call abort elseif(plevrad(L_LEVELS).gt.pgasmax)then print*,'Maximum pressure is outside the radiative' print*,'transfer kmatrix bounds, exiting.' ! call abort 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, exiting.' print*,tlevrad(k),k,tgasmin ! call abort elseif(tlevrad(k).gt.tgasmax)then print*,'Maximum temperature is outside the radiative' print*,'transfer kmatrix bounds, exiting.' ! call abort endif enddo !======================================================================= ! Calling the main radiative transfer subroutines !----------------------------------------------------------------------- ! Shortwave IF(fract(ig) .GE. 1.0e-4) THEN ! only during daylight IPM?! flux UV... fluxtoplanet=0. DO nw=1,L_NSPECTV stel_fract(nw)= stel(nw) * fract(ig) fluxtoplanet=fluxtoplanet + stel_fract(nw) END DO !print*, 'starting optcv' call optcv(dtauv,tauv,taucumv,plevrad, $ qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero, $ tmid,pmid,taugsurf,qvar) call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv, $ acosz,stel_fract,gweight,nfluxtopv,nfluxtopv_nu, $ fmnetv,fmnetv_nu,fluxupv,fluxdnv,fzerov,taugsurf) ELSE ! during the night, fluxes = 0 nfluxtopv=0.0 DO l=1,L_NLAYRAD fmnetv(l)=0.0 fluxupv(l)=0.0 fluxdnv(l)=0.0 END DO END IF !----------------------------------------------------------------------- ! Longwave call optci(plevrad,tlevrad,dtaui,taucumi, $ qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid, $ taugsurfi,qvar) call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, $ wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, $ fmneti,fmneti_nu,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi) !----------------------------------------------------------------------- ! Transformation of the correlated-k code outputs ! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw) fluxtop_lw(ig) = fluxupi(1) fluxsurf_lw(ig) = fluxdni(L_NLAYRAD) fluxtop_sw(ig) = fluxupv(1) fluxsurf_sw(ig) = fluxdnv(L_NLAYRAD) ! Flux incident at the top of the atmosphere fluxtop_dn(ig)=fluxdnv(1) ! IR spectral output from top of the atmosphere if(specOLR)then do nw=1,L_NSPECTI OLR_nu(ig,nw)=nfluxtopi_nu(nw) end do endif ! ********************************************************** ! Finally, the heating rates ! g/cp*DF/DP ! ********************************************************** DO l=2,L_NLAYRAD dpp = g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) ! DTSW : !dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1))*dpp !averaged dtlw on each wavelength do nw=1,L_NSPECTV dtsw_nu(L_NLAYRAD+1-l,nw)= & (fmnetv_nu(l,nw)-fmnetv_nu(l-1,nw))*dpp end do ! DTLW : detailed with wavelength so that we can apply NLTE !dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1))*dpp !averaged dtlw on each wavelength do nw=1,L_NSPECTI dtlw_nu(L_NLAYRAD+1-l,nw)= & (fmneti_nu(l,nw)-fmneti_nu(l-1,nw))*dpp end do END DO ! values at top of atmosphere dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1))) ! SW !dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv)*dpp do nw=1,L_NSPECTV dtsw_nu(L_NLAYRAD,nw)= & (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp end do ! LW c dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) *dpp do nw=1,L_NSPECTI dtlw_nu(L_NLAYRAD,nw)= & (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp end do ! ********************************************************** ! NON NLTE correction in Pluto atmosphere ! And conversion of LW spectral heating rates to total rates ! ********************************************************** if (.not.nlte) then eps_nlte_sw23(ig,:) =1. ! IF no NLTE eps_nlte_sw33(ig,:) =1. ! IF no NLTE eps_nlte_lw(ig,:) =1. ! IF no NLTE endif do l=1,nlayer !LW dtlw(ig,l) =0. ! dtlw_co(ig,l) =0. ! only for diagnostic do nw=1,L_NSPECTI ! wewei : wavelength in micrometer if ((wavei(nw).gt.6.).and.(wavei(nw).lt.9)) then dtlw_nu(l,nw)=dtlw_nu(l,nw)*eps_nlte_lw(ig,l) else !dtlw_nu(l,nw)=1.*dtlw_nu(l,nw) ! no CO correction (Strobbel 1996) dtlw_nu(l,nw)=0.33*dtlw_nu(l,nw) ! CO correction (Strobbel 1996) ! dtlw_co(ig,l)=dtlw_co(ig,l)+ dtlw_nu(l,nw) ! diagnostic end if dtlw(ig,l)=dtlw(ig,l)+ dtlw_nu(l,nw) !average now on each wavelength end do ! adding c2h2 if cooling active dtlw(ig,l)=dtlw(ig,l)+dtlw_hcn_c2h2(ig,l) !SW dtsw(ig,l) =0. if (strobel) then do nw=1,L_NSPECTV if ((wavev(nw).gt.2).and.(wavev(nw).lt.2.6)) then dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l) elseif ((wavev(nw).gt.3).and.(wavev(nw).lt.3.6)) then dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw33(ig,l) else dtsw_nu(l,nw)=dtsw_nu(l,nw) end if dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw) end do else ! total heating rates multiplied by nlte coef do nw=1,L_NSPECTV dtsw_nu(l,nw)=dtsw_nu(l,nw)*eps_nlte_sw23(ig,l) dtsw(ig,l)=dtsw(ig,l)+ dtsw_nu(l,nw) enddo endif end do ! ********************************************************** ! Diagnotics for last call for each grid point !if (lastcall) then !print*,'albedi vis=',albv !print*,'albedo ir=',albi !print*,'fluxup ir (:)=',fluxupi !print*,'flux ir net (:)=',fluxdni-fluxupi !print*,'cumulative flux net ir (:)=',fmneti !print*,'cumulative flux net vis (:)=',fmnetv !print*,'fluxdn vis (:)=',fluxdnv !print*,'fluxtop vis=',fluxtop_sw !print*,'fluxsurf vis=',fluxsurf_sw !print*,'fluxtop ir=',fluxtop_lw !print*,'fluxsurf ir=',fluxsurf_lw c write(*,*) 'pressure, eps_nlte_sw, eps_nlte_lw' c do l=1,nlayer c write(*,*)pplay(1,l),eps_nlte_sw(1,l),eps_nlte_lw(1,l) c end do !endif ! --------------------------------------------------------------- end do ! end of big loop over every GCM column (ig = 1:ngrid) !----------------------------------------------------------------------- ! Additional diagnostics ! IR spectral output, for exoplanet observational comparison if(specOLR)then if(ngrid.ne.1)then call writediagspec(ngrid,"OLR3D", & "OLR(lon,lat,band)","W m^-2",3,OLR_nu) endif endif if(lastcall)then ! 1D Output if(ngrid.eq.1)then ! surface diagnotics diagrad_surf=.true. if(diagrad_surf)then open(116,file='surf_vals.out') write(116,*) tsurf(1),pplev(1,1), & fluxtop_dn(1) - fluxtop_sw(1),fluxtop_lw(1) do nw=1,L_NSPECTV write(116,*) wavev(nw),fmnetv_nu(L_NLAYRAD,nw) enddo do nw=1,L_NSPECTI write(116,*) wavei(nw),fmneti_nu(L_NLAYRAD,nw) enddo close(116) endif ! OLR by band diagrad_OLR=.true. if(diagrad_OLR)then open(117,file='OLRnu.out') write(117,*) 'IR wavel - band width - OLR' do nw=1,L_NSPECTI write(117,*) wavei(nw), & abs(1.e4/bwnv(nw)-1.e4/bwnv(nw+1)),OLR_nu(1,nw) enddo close(117) endif ! OLR vs altitude: in a .txt file diagrad_OLRz=.true. if(diagrad_OLRz)then open(118,file='OLRz_plevs.out') open(119,file='OLRz.out') do l=1,L_NLAYRAD write(118,*) plevrad(2*l) do nw=1,L_NSPECTI write(119,*) fluxupi_nu(l,nw) enddo enddo close(118) close(119) endif ! Heating rates vs altitude in a .txt file diagrad_rates=.true. if(diagrad_rates)then open(120,file='heating_rates.out') write(120,*) "Pressure - Alt - HR tot - Rates (wavel SW)" do l=1,nlayermx write(120,*) pplay(1,l),zzlay(1,l),dtsw(1,l),dtsw_nu(l,:) enddo close(120) open(121,file='cooling_rates.out') write(121,*) "Pressure - Alt - CR tot - Rates (wavel LW)" do l=1,nlayermx write(121,*) pplay(1,l),zzlay(1,l),dtlw(1,l),dtlw_nu(l,:) enddo close(121) open(122,file='bands.out') write(122,*) "wavel - bands boundaries (microns)" do nw=1,L_NSPECTV write(122,*) wavev(nw),1.e4/bwnv(nw+1),1.e4/bwnv(nw) enddo do nw=1,L_NSPECTI write(122,*) wavei(nw),1.e4/bwni(nw+1),1.e4/bwni(nw) enddo close(122) open(123,file='c2h2_rates.out') write(123,*) "Pressure - Alt - CR c2h2" do l=1,nlayermx write(123,*) pplay(1,l),zzlay(1,l),dtlw_hcn_c2h2(1,l) enddo close(123) endif endif ! ngrid.eq.1 endif ! lastcall end