subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,emis,mu0,pplev,pplay,pt, & tsurf,fract,dist_star,aerosol,muvar, & dtlw,dtsw,fluxsurf_lw, & fluxsurf_sw,fluxtop_lw,fluxabs_sw,fluxtop_dn, & OLR_nu,OSR_nu, & reffrad,tau_col,cloudfrac,totcloudfrac, & clearsky,firstcall,lastcall) use radinc_h use radcommon_h use watercommon_h use datafile_mod, only: datadir use ioipsl_getincom use gases_h 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) ! !================================================================== #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 REAL aerosol(ngrid,nlayermx,naerkind) ! aerosol tau (kg/kg) REAL albedo(ngrid) ! SW albedo REAL emis(ngrid) ! LW emissivity REAL pplay(ngrid,nlayermx) ! pres. level in GCM mid of layer 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) ! 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 fluxabs_sw(ngridmx) ! SW flux absorbed by planet (W/m2) REAL fluxtop_dn(ngridmx) ! incident top of atmosphere SW flux (W/m2) REAL OLR_nu(ngrid,L_NSPECTI)! Outgoing LW radition in each band (Normalized to the band width (W/m2/cm-1) REAL OSR_nu(ngrid,L_NSPECTV)! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1) !----------------------------------------------------------------------- ! 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 nfluxoutv_nu(L_NSPECTV) ! outgoing band-resolved VI flux at TOA (W/m2) real*8 nfluxtopi_nu(L_NSPECTI) ! net band-resolved IR flux at TOA (W/m2) real*8 fluxupi_nu(L_NLAYRAD,L_NSPECTI) ! for 1D diagnostic REAL*8 fmneti(L_NLAYRAD),fmnetv(L_NLAYRAD) 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 szangle logical global1d save szangle,global1d 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 (mol/mol) 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 water vapour profiles integer i_var real RH real*8 pq_temp(nlayer) real ptemp, Ttemp, qsat ! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!! !real ptime, pday logical OLRz real*8 NFLUXGNDV_nu(L_NSPECTV) ! for H2O cloud fraction in aeropacity real cloudfrac(ngridmx,nlayermx) real totcloudfrac(ngridmx) logical clearsky ! for weird cloud test real pqtest(ngridmx,nlayer,nq) ! are particle radii fixed? logical radfixed real maxrad, minrad real CBRT external CBRT ! included by RW for runaway greenhouse 1D study real muvar(ngridmx,nlayermx+1) real vtmp(nlayermx) REAL*8 muvarrad(L_LEVELS) radfixed=.false. !======================================================================= ! Initialization on first call qxvaer(:,:,:)=0.0 qsvaer(:,:,:)=0.0 gvaer(:,:,:) =0.0 qxiaer(:,:,:)=0.0 qsiaer(:,:,:)=0.0 giaer(:,:,:) =0.0 if(firstcall) then call system('rm -f surf_vals_long.out') !-------------------------------------------------- ! Effective radius and variance of the aerosols do iaer=1,naerkind ! these values will change once the microphysics gets to work ! UNLESS tracer=.false., in which case we should be working with ! a fixed aerosol layer, and be able to define reffrad in a ! .def file. To be improved! if(iaer.eq.1)then ! CO2 ice do l=1,nlayer do ig=1,ngrid reffrad(ig,l,iaer) = 1.e-4 nueffrad(ig,l,iaer) = 0.1 enddo enddo endif if(iaer.eq.2)then ! H2O ice do l=1,nlayer do ig=1,ngrid reffrad(ig,l,iaer) = 1.e-5 nueffrad(ig,l,iaer) = 0.1 enddo enddo endif if(iaer.eq.3)then ! dust do l=1,nlayer do ig=1,ngrid reffrad(ig,l,iaer) = 1.e-5 nueffrad(ig,l,iaer) = 0.1 enddo enddo endif if(iaer.gt.3)then print*,'Error in callcorrk, naerkind is too high.' print*,'The code still needs generalisation to arbitrary' print*,'aerosol kinds and number.' call abort endif enddo print*, "callcorrk: Correlated-k data base 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)) call sugas_corrk ! set up gaseous absorption properties call setspi ! basic infrared properties call setspv ! basic visible properties call suaer_corrk ! set up aerosol optical properties Cmk= 0.01 * 1.0 / (g * mugaz * 1.672621e-27) ! q_main=1.0 assumed if((igcm_h2o_vap.eq.0) .and. varactive)then print*,'varactive in callcorrk but no h2o_vap tracer.' stop endif OLR_nu=0. OSR_nu=0. if (ngridmx.eq.1) then PRINT*, 'Simulate global averaged conditions ?' global1d = .false. ! default value call getin("global1d",global1d) write(*,*) "global1d = ",global1d if (global1d) then PRINT *,'Solar Zenith angle (deg.) ?' PRINT *,'(assumed for averaged solar flux S/4)' szangle=60.0 ! default value call getin("szangle",szangle) write(*,*) "szangle = ",szangle endif endif firstcall=.false. end if !======================================================================= ! Initialization on every call do l=1,nlayer do ig=1,ngrid do iaer=1,naerkind nueffrad(ig,l,iaer) = 0.1 ! stays at 0.1 enddo enddo enddo if(kastprof)then radfixed=.true. endif if(radfixed)then do l=1,nlayer do ig=1,ngrid reffrad(ig,l,1) = 5.e-5 ! CO2 ice enddo enddo print*,'CO2 ice particle size = ',reffrad(1,1,1)/1.e-6,' um' if(naerkind.ge.2)then do l=1,nlayer do ig=1,ngrid !reffrad(ig,l,2) = 2.e-5 ! H2O ice reffrad(ig,l,2) = 5.e-6 ! H2O ice enddo enddo print*,'H2O ice particle size = ',reffrad(1,1,2)/1.e-6,' um' endif if(naerkind.eq.3)then do l=1,nlayer do ig=1,ngrid reffrad(ig,l,naerkind) = 2.e-6 ! dust enddo enddo print*,'Dust particle size = ',reffrad(1,1,naerkind)/1.e-6,' um' endif if(naerkind.gt.3)then print*,'Code not general enough to deal with naerkind > 3 yet.' call abort endif else maxrad=0.0 !averad=0.0 minrad=1.0 do l=1,nlayer !masse = (pplev(ig,l) - pplev(ig,l+1))/g do ig=1,ngrid !if(tracer)then if(tracer.and.igcm_co2_ice.gt.0)then if(igcm_co2_ice.lt.1)then print*,'Tracers but no CO2 ice still seems to be a problem...' print*,'Aborting in callcorrk.' stop endif reffrad(ig,l,1) = CBRT( 3*pq(ig,l,igcm_co2_ice)/ & (4*Nmix_co2*pi*rho_co2) ) endif reffrad(ig,l,1) = max(reffrad(ig,l,1),1.e-6) reffrad(ig,l,1) = min(reffrad(ig,l,1),500.e-6) !averad = averad + reffrad(ig,l,1)*area(ig) maxrad = max(reffrad(ig,l,1),maxrad) minrad = min(reffrad(ig,l,1),minrad) enddo enddo if(igcm_co2_ice.gt.0)then print*,'Max. CO2 ice particle size = ',maxrad/1.e-6,' um' print*,'Min. CO2 ice particle size = ',minrad/1.e-6,' um' endif if((naerkind.ge.2).and.water)then maxrad=0.0 minrad=1.0 do l=1,nlayer do ig=1,ngrid reffrad(ig,l,2) = CBRT( 3*pq(ig,l,igcm_h2o_ice)/ & (4*Nmix_h2o*pi*rho_ice) ) reffrad(ig,l,2) = max(reffrad(ig,l,2),1.e-6) reffrad(ig,l,2) = min(reffrad(ig,l,2),100.e-6) maxrad = max(reffrad(ig,l,2),maxrad) minrad = min(reffrad(ig,l,2),minrad) enddo enddo print*,'Max. H2O ice particle size = ',maxrad/1.e-6,' um' print*,'Min. H2O ice particle size = ',minrad/1.e-6,' um' endif if(naerkind.eq.3)then do l=1,nlayer do ig=1,ngrid reffrad(ig,l,naerkind) = 2.e-6 ! dust enddo enddo endif endif ! how much light we get fluxtoplanet=0 do nw=1,L_NSPECTV stel(nw)=stellarf(nw)/(dist_star**2) fluxtoplanet=fluxtoplanet + stel(nw) end do call aeroptproperties(ngrid,nlayer,reffrad,nueffrad, & QVISsQREF3d,omegaVIS3d,gVIS3d, & QIRsQREF3d,omegaIR3d,gIR3d, & QREFvis3d,QREFir3d) ! get 3D aerosol optical properties call aeropacity(ngrid,nlayer,nq,pplay,pplev,pq,aerosol, & reffrad,QREFvis3d,QREFir3d, & tau_col,cloudfrac,totcloudfrac,clearsky) ! get aerosol optical depths !----------------------------------------------------------------------- ! Starting Big Loop over every GCM column do ig=1,ngridmx !======================================================================= ! Transformation of the GCM variables !----------------------------------------------------------------------- ! Aerosol optical properties Qext, Qscat and g ! The transformation in the vertical is the same as for temperature ! 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 ! test / correct for freaky s. s. albedo values do iaer=1,naerkind do k=1,L_LEVELS+1 do nw=1,L_NSPECTV if(qsvaer(k,nw,iaer).gt.1.05*qxvaer(k,nw,iaer))then print*,'Serious problems with qsvaer values 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 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 end do !----------------------------------------------------------------------- ! Aerosol optical depths do iaer=1,naerkind ! a bug was here 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)) 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(L_LEVELS+1,iaer) = tauaero(L_LEVELS,iaer) !tauaero(1,iaer) = 0. !tauaero(L_LEVELS+1,iaer) = 0. end do ! Albedo and emissivity albi=1-emis(ig) ! longwave albv=albedo(ig) ! shortwave if(noradsurf.and.(albv.gt.0.0))then print*,'For open lower boundary in callcorrk must' print*,'have surface albedo set to zero!' call abort endif if ((ngridmx.eq.1).and.(global1d)) then ! fixed zenith angle 'szangle' in 1D simulations w/ globally-averaged sunlight acosz = cos(pi*szangle/180.0) print*,'acosz=',acosz,', szangle=',szangle else acosz=mu0(ig) ! cosine of sun incident angle : 3D simulations or local 1D simulations using latitude endif !----------------------------------------------------------------------- ! Water vapour (to be generalised for other gases eventually) if(varactive)then i_var=igcm_h2o_vap do l=1,nlayer qvar(2*l) = pq(ig,nlayer+1-l,i_var) qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2 ! Average approximation as for temperature... end do qvar(1)=qvar(2) elseif(varfixed)then do l=1,nlayermx ! here we will assign fixed water vapour profiles globally RH = satval * ((pplay(ig,l)/pplev(ig,1) - 0.02) / 0.98) if(RH.lt.0.0) RH=0.0 ptemp=pplay(ig,l) Ttemp=pt(ig,l) call watersat(Ttemp,ptemp,qsat) !pq_temp(l) = qsat ! fully saturated everywhere pq_temp(l) = RH * qsat ! ~realistic profile (e.g. 80% saturation at ground) end do do l=1,nlayer qvar(2*l) = pq_temp(nlayer+1-l) qvar(2*l+1) = (pq_temp(nlayer+1-l)+pq_temp(max(nlayer-l,1)))/2 end do qvar(1)=qvar(2) ! Lowest layer of atmosphere RH = satval * (1 - 0.02) / 0.98 if(RH.lt.0.0) RH=0.0 ptemp = pplev(ig,1) Ttemp = tsurf(ig) call watersat(Ttemp,ptemp,qsat) !qvar(2*nlayermx+1)=qsat ! fully saturated everywhere qvar(2*nlayermx+1)= RH * qsat ! ~realistic profile (e.g. 80% saturation at ground) !qvar=0.005 ! completely constant profile (JL) else do k=1,L_LEVELS qvar(k) = 1.0D-7 end do end if if(.not.kastprof)then ! IMPORTANT: Now convert from kg/kg to mol/mol do k=1,L_LEVELS qvar(k) = qvar(k)/epsi end do end if !----------------------------------------------------------------------- ! kcm mode only if(kastprof)then ! initial values equivalent to mugaz DO l=1,nlayer muvarrad(2*l) = mugaz muvarrad(2*l+1) = mugaz END DO !do k=1,L_LEVELS ! qvar(k) = 0.0 !end do !print*,'ASSUMING qH2O=0 EVERYWHERE IN CALLCORRK!' endif if(kastprof.and.(ngasmx.gt.1))then DO l=1,nlayer muvarrad(2*l) = muvar(ig,nlayer+2-l) muvarrad(2*l+1) = (muvar(ig,nlayer+2-l) + & muvar(ig,max(nlayer+1-l,1)))/2 END DO muvarrad(1) = muvarrad(2) muvarrad(2*nlayermx+1)=muvar(ig,1) print*,'Recalculating qvar with VARIABLE epsi for kastprof' print*,'Assumes that the variable gas is H2O!!!' print*,'Assumes that there is only one tracer' !i_var=igcm_h2o_vap i_var=1 if(nqmx.gt.1)then print*,'Need 1 tracer only to run kcm1d.e' stop endif do l=1,nlayer vtmp(l)=pq(ig,l,i_var)*muvar(ig,l+1)/mH2O end do do l=1,nlayer qvar(2*l) = vtmp(nlayer+1-l) qvar(2*l+1) = ( vtmp(nlayer+1-l) + vtmp(max(nlayer-l,1)) )/2 end do qvar(1)=qvar(2) print*,'Warning: reducing qvar in callcorrk.' print*,'Temperature profile no longer consistent ', & 'with saturated H2O.' do k=1,L_LEVELS qvar(k) = qvar(k)*satval end do endif ! 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 !----------------------------------------------------------------------- ! 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) = 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 pmid(L_LEVELS) = plevrad(L_LEVELS) tmid(L_LEVELS) = tlevrad(L_LEVELS) ! 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*,'WARNING, OVERRIDING FOR TEST' call abort elseif(tlevrad(k).gt.tgasmax)then print*,'Maximum temperature is outside the radiative' print*,'transfer kmatrix bounds, exiting.' !print*,'WARNING, OVERRIDING FOR TEST' call abort endif enddo !======================================================================= ! Calling the main radiative transfer subroutines !----------------------------------------------------------------------- ! Shortwave if(fract(ig) .ge. 1.0e-4) then ! only during daylight! fluxtoplanet=0. if((ngridmx.eq.1).and.(global1d))then do nw=1,L_NSPECTV stel_fract(nw)= stel(nw) * 0.25 / acosz fluxtoplanet=fluxtoplanet + stel_fract(nw) ! globally averaged = divide by 4 ! but we correct for solar zenith angle end do else do nw=1,L_NSPECTV stel_fract(nw)= stel(nw) * fract(ig) fluxtoplanet=fluxtoplanet + stel_fract(nw) end do endif 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,gweight, & nfluxtopv,nfluxoutv_nu,nfluxgndv_nu, & fmnetv,fluxupv,fluxdnv,fzerov,taugsurf) else ! during the night, fluxes = 0 nfluxtopv = 0.0 nfluxoutv_nu(:) = 0.0 nfluxgndv_nu(:) = 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,muvarrad) call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi, & wnoi,dwni,cosbi,wbari,gweight,nfluxtopi,nfluxtopi_nu, & fmneti,fluxupi,fluxdni,fluxupi_nu,fzeroi,taugsurfi) !----------------------------------------------------------------------- ! Transformation of the correlated-k code outputs ! (into dtlw, dtsw, fluxsurf_lw, fluxsurf_sw, fluxtop_lw, fluxtop_sw) ! Flux incident at the top of the atmosphere fluxtop_dn(ig)=fluxdnv(1) fluxtop_lw(ig) = real(nfluxtopi) fluxabs_sw(ig) = real(-nfluxtopv) fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD)) 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,:) call abort endif ! Spectral output, for exoplanet observational comparison if(specOLR)then do nw=1,L_NSPECTI OLR_nu(ig,nw)=nfluxtopi_nu(nw)/DWNI(nw) !JL Normalize to the bandwidth end do do nw=1,L_NSPECTV !GSR_nu(ig,nw)=nfluxgndv_nu(nw) OSR_nu(ig,nw)=nfluxoutv_nu(nw)/DWNV(nw) !JL Normalize to the bandwidth end do endif ! Finally, the heating rates DO l=2,L_NLAYRAD dtsw(ig,L_NLAYRAD+1-l)=(fmnetv(l)-fmnetv(l-1)) & *g/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) dtlw(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) & *g/(cpp*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) & *g/(cpp*scalep*(plevrad(3)-plevrad(1))) dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) & *g/(cpp*scalep*(plevrad(3)-plevrad(1))) ! --------------------------------------------------------------- end do ! end of big loop over every GCM column (ig = 1:ngrid) !----------------------------------------------------------------------- ! Additional diagnostics ! IR spectral output, for exoplanet observational comparison if(lastcall.and.(ngrid.eq.1))then ! could disable the 1D output, they are in the diagfi and diagspec... JL12 print*,'Saving scalar quantities in surf_vals.out...' print*,'psurf = ', pplev(1,1),' Pa' open(116,file='surf_vals.out') write(116,*) tsurf(1),pplev(1,1),fluxtop_dn(1), & real(-nfluxtopv),real(nfluxtopi) close(116) ! I am useful, please don`t remove me! ! if(specOLR)then ! open(117,file='OLRnu.out') ! do nw=1,L_NSPECTI ! write(117,*) OLR_nu(1,nw) ! enddo ! close(117) ! ! open(127,file='OSRnu.out') ! do nw=1,L_NSPECTV ! write(127,*) OSR_nu(1,nw) ! enddo ! close(127) ! endif ! OLR vs altitude: do it as a .txt file OLRz=.false. if(OLRz)then print*,'saving IR vertical flux for OLRz...' 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 endif ! see physiq.F for explanations about CLFvarying. This is temporary. if (lastcall .and. .not.CLFvarying) then IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi ) IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv ) IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref ) IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref ) IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar ) IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref ) endif end subroutine callcorrk