MODULE callcorrk_mod IMPLICIT NONE CONTAINS subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf, & albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt, & zzlay,tsurf,fract,dist_star,aerosol,muvar, & dtlw,dtsw,fluxsurf_lw, & fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw, & fluxabs_sw,fluxtop_dn, & OLR_nu,OSR_nu,GSR_nu, & int_dtaui,int_dtauv, & tau_col,cloudfrac,totcloudfrac, & clearsky,firstcall,lastcall) use mod_phys_lmdz_para, only : is_master use radinc_h, only: L_NSPECTV, L_NSPECTI, naerkind, banddir, corrkdir,& L_LEVELS, L_NGAUSS, L_NLEVRAD, L_NLAYRAD, L_REFVAR use radcommon_h, only: wrefvar, Cmk, fzeroi, fzerov, gasi, gasv, & glat_ig, gweight, pfgasref, pgasmax, pgasmin, & pgasref, tgasmax, tgasmin, tgasref, scalep, & ubari, wnoi, stellarf, glat, dwnv, dwni, tauray use datafile_mod, only: datadir use ioipsl_getin_p_mod, only: getin_p use gases_h, only: ngasmx use radii_mod, only : su_aer_radii, haze_reffrad_fix use aerosol_mod, only : iaero_haze use aeropacity_mod, only: aeropacity use aeroptproperties_mod, only: aeroptproperties use tracer_h, only: constants_epsi_generic,igcm_ch4_gas,igcm_n2,mmol use comcstfi_mod, only: pi, mugaz, cpp, r, g use callkeys_mod, only: varactive,diurnal,tracer,varfixed,satval, & diagdtau,kastprof,strictboundcorrk,specOLR, & tplanckmin,tplanckmax,global1d, & generic_condensation,aerohaze,haze_radproffix,& methane,carbox,cooling,nlte,strobel,& ch4fix,vmrch4_proffix,vmrch4fix use optcv_mod, only: optcv use optci_mod, only: optci use sfluxi_mod, only: sfluxi use sfluxv_mod, only: sfluxv use recombin_corrk_mod, only: corrk_recombin, call_recombin use generic_cloud_common_h, only: Psat_generic, epsi_generic use generic_tracer_index_mod, only: generic_tracer_index use planetwide_mod, only: planetwide_maxval, planetwide_minval use radcommon_h, only: wavev,wavei 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) ! !================================================================== !----------------------------------------------------------------------- ! 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) :: pq(ngrid,nlayer,nq) ! Tracers (kg/kg_of_air). INTEGER,INTENT(IN) :: nq ! Number of tracers. REAL,INTENT(IN) :: qsurf(ngrid,nq) ! Tracers on surface (kg.m-2). REAL,INTENT(IN) :: albedo(ngrid,L_NSPECTV) ! Spectral Short Wavelengths Albedo. By MT2015 REAL,INTENT(IN) :: emis(ngrid) ! Long Wave emissivity. 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) :: zzlay(ngrid,nlayer) ! Mid-layer altitude 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). REAL,INTENT(IN) :: muvar(ngrid,nlayer+1) REAL,INTENT(IN) :: cloudfrac(ngrid,nlayer) ! Fraction of clouds (%). logical,intent(in) :: clearsky logical,intent(in) :: firstcall ! Signals first call to physics. logical,intent(in) :: lastcall ! Signals last call to physics. ! OUTPUT REAL,INTENT(OUT) :: aerosol(ngrid,nlayer,naerkind) ! Aerosol tau at reference wavelenght. REAL,INTENT(OUT) :: dtlw(ngrid,nlayer) ! Heating rate (K/s) due to LW radiation. REAL,INTENT(OUT) :: dtsw(ngrid,nlayer) ! Heating rate (K/s) due to SW radiation. REAL,INTENT(OUT) :: fluxsurf_lw(ngrid) ! Incident LW flux to surf (W/m2). REAL,INTENT(OUT) :: fluxsurf_sw(ngrid) ! Incident SW flux to surf (W/m2) REAL,INTENT(OUT) :: fluxsurfabs_sw(ngrid) ! Absorbed SW flux by the surface (W/m2). By MT2015. REAL,INTENT(OUT) :: fluxtop_lw(ngrid) ! Outgoing LW flux to space (W/m2). REAL,INTENT(OUT) :: fluxabs_sw(ngrid) ! SW flux absorbed by the planet (W/m2). REAL,INTENT(OUT) :: fluxtop_dn(ngrid) ! Incident top of atmosphere SW flux (W/m2). REAL,INTENT(OUT) :: OLR_nu(ngrid,L_NSPECTI) ! Outgoing LW radiation in each band (Normalized to the band width (W/m2/cm-1). REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV) ! Outgoing SW radiation in each band (Normalized to the band width (W/m2/cm-1). REAL,INTENT(OUT) :: GSR_nu(ngrid,L_NSPECTV) ! Surface SW radiation in each band (Normalized to the band width (W/m2/cm-1). REAL,INTENT(OUT) :: tau_col(ngrid) ! Diagnostic from aeropacity. REAL,INTENT(OUT) :: albedo_equivalent(ngrid) ! Spectrally Integrated Albedo. For Diagnostic. By MT2015 REAL,INTENT(OUT) :: totcloudfrac(ngrid) ! Column Fraction of clouds (%). REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags (). REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags (). ! Globally varying aerosol optical properties on GCM grid ; not needed everywhere so not in radcommon_h. ! made "save" variables so they are allocated once in for all, not because ! the values need be saved from a time step to the next REAL,SAVE,ALLOCATABLE :: QVISsQREF3d(:,:,:,:) REAL,SAVE,ALLOCATABLE :: omegaVIS3d(:,:,:,:) REAL,SAVE,ALLOCATABLE :: gVIS3d(:,:,:,:) !$OMP THREADPRIVATE(QVISsQREF3d,omegaVIS3d,gVIS3d) REAL,SAVE,ALLOCATABLE :: QIRsQREF3d(:,:,:,:) REAL,SAVE,ALLOCATABLE :: omegaIR3d(:,:,:,:) REAL,SAVE,ALLOCATABLE :: gIR3d(:,:,:,:) !$OMP THREADPRIVATE(QIRsQREF3d,omegaIR3d,gIR3d) ! REAL :: omegaREFvis3d(ngrid,nlayer,naerkind) ! REAL :: omegaREFir3d(ngrid,nlayer,naerkind) ! not sure of the point of these... 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 :: dtaui(:,:,:) REAL*8,allocatable,save :: dtauv(:,:,:) REAL*8,allocatable,save :: cosbv(:,:,:) REAL*8,allocatable,save :: cosbi(:,:,:) REAL*8,allocatable,save :: wbari(:,:,:) REAL*8,allocatable,save :: wbarv(:,:,:) !$OMP THREADPRIVATE(dtaui,dtauv,cosbv,cosbi,wbari,wbarv) REAL*8,allocatable,save :: tauv(:,:,:) REAL*8,allocatable,save :: taucumv(:,:,:) REAL*8,allocatable,save :: taucumi(:,:,:) !$OMP THREADPRIVATE(tauv,taucumv,taucumi) REAL*8,allocatable,save :: tauaero(:,:) !$OMP THREADPRIVATE(tauaero) REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn real*8 nfluxtopv_nu(L_NSPECTV) 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 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,acosz REAL*8 albv(L_NSPECTV) ! Spectral Visible Albedo. INTEGER ig,l,k,nw,iaer,iq real*8,allocatable,save :: taugsurf(:,:) real*8,allocatable,save :: taugsurfi(:,:) !$OMP THREADPRIVATE(taugsurf,taugsurfi) real*8 qvar(L_LEVELS) ! Mixing ratio of variable component (mol/mol). index 1 is the top of the atmosphere, index L_LEVELS is the bottom ! 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(:,:,:) real*8,save,allocatable :: QXIAER(:,:,:) ! Extinction coeff (QIRsQREF*QREFir) real*8,save,allocatable :: QSIAER(:,:,:) real*8,save,allocatable :: GIAER(:,:,:) !$OMP THREADPRIVATE(QXVAER,QSVAER,GVAER,QXIAER,QSIAER,GIAER) real, dimension(:,:,:), save, allocatable :: QREFvis3d real, dimension(:,:,:), save, allocatable :: QREFir3d !$OMP THREADPRIVATE(QREFvis3d,QREFir3d) ! Miscellaneous : real*8 temp,temp1,temp2,pweight character(len=10) :: tmp1 character(len=10) :: tmp2 character(len=100) :: message character(len=10),parameter :: subname="callcorrk" ! For fixed water vapour profiles. integer i_var real RH real*8 pq_temp(nlayer) ! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!! 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) ! Included by MT for albedo calculations. REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation. REAL*8 surface_stellar_flux ! Stellar flux reaching the surface. Useful for equivalent albedo calculation. ! NLTE factor for CH4 real eps_nlte_sw23(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw real eps_nlte_sw33(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw real eps_nlte_lw(ngrid,nlayer) ! CH4 NLTE efficiency factor for zdtsw integer Nfine,ifine parameter(Nfine=701) real,save :: levdat(Nfine),vmrdat(Nfine) REAL dtlw_hcn_c2h2(ngrid, nlayer) ! cooling rate (K/s) due to C2H2/HCN (diagnostic) real :: vmrch4(ngrid,nlayer) ! vmr ch4 from vmrch4_proffix REAL dtlw_nu(nlayer,L_NSPECTI) ! heating rate (K/s) due to LW in spectral bands REAL dtsw_nu(nlayer,L_NSPECTV) ! heating rate (K/s) due to SW in spectral bands ! local variable REAL dpp ! intermediate integer ok ! status (returned by NetCDF functions) integer igcm_generic_gas, igcm_generic_ice! index of the vap and ice of generic_tracer logical call_ice_gas_generic ! to call only one time the ice/vap pair of a tracer real, save :: metallicity ! metallicity of planet --- is not used here, but necessary to call function Psat_generic !$OMP THREADPRIVATE(metallicity) REAL, SAVE :: qvap_deep ! deep mixing ratio of water vapor when simulating bottom less planets !$OMP THREADPRIVATE(qvap_deep) REAL :: maxvalue,minvalue !=============================================================== ! I.a Initialization on first call !=============================================================== if(firstcall) then ! test on allocated necessary because of CLFvarying (two calls to callcorrk in physiq) if(.not.allocated(QVISsQREF3d)) then allocate(QVISsQREF3d(ngrid,nlayer,L_NSPECTV,naerkind)) endif if(.not.allocated(omegaVIS3d)) then allocate(omegaVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)) endif if(.not.allocated(gVIS3d)) then allocate(gVIS3d(ngrid,nlayer,L_NSPECTV,naerkind)) endif if (.not.allocated(QIRsQREF3d)) then allocate(QIRsQREF3d(ngrid,nlayer,L_NSPECTI,naerkind)) endif if (.not.allocated(omegaIR3d)) then allocate(omegaIR3d(ngrid,nlayer,L_NSPECTI,naerkind)) endif if (.not.allocated(gIR3d)) then allocate(gIR3d(ngrid,nlayer,L_NSPECTI,naerkind)) endif if (.not.allocated(tauaero)) then allocate(tauaero(L_LEVELS,naerkind)) endif if(.not.allocated(QXVAER)) then allocate(QXVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok) if (ok /= 0) then write(*,*) "memory allocation failed for QXVAER!" call abort_physic(subname,'allocation failure for QXVAER',1) endif endif if(.not.allocated(QSVAER)) then allocate(QSVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok) if (ok /= 0) then write(*,*) "memory allocation failed for QSVAER!" call abort_physic(subname,'allocation failure for QSVAER',1) endif endif if(.not.allocated(GVAER)) then allocate(GVAER(L_LEVELS,L_NSPECTV,naerkind), stat=ok) if (ok /= 0) then write(*,*) "memory allocation failed for GVAER!" call abort_physic(subname,'allocation failure for GVAER',1) endif endif if(.not.allocated(QXIAER)) then allocate(QXIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok) if (ok /= 0) then write(*,*) "memory allocation failed for QXIAER!" call abort_physic(subname,'allocation failure for QXIAER',1) endif endif if(.not.allocated(QSIAER)) then allocate(QSIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok) if (ok /= 0) then write(*,*) "memory allocation failed for QSIAER!" call abort_physic(subname,'allocation failure for QSIAER',1) endif endif if(.not.allocated(GIAER)) then allocate(GIAER(L_LEVELS,L_NSPECTI,naerkind), stat=ok) if (ok /= 0) then write(*,*) "memory allocation failed for GIAER!" call abort_physic(subname,'allocation failure for GIAER',1) endif endif !!! ALLOCATED instances are necessary because of CLFvarying (strategy to call callcorrk twice in physiq...) IF(.not.ALLOCATED(QREFvis3d))THEN ALLOCATE(QREFvis3d(ngrid,nlayer,naerkind), stat=ok) IF (ok/=0) THEN write(*,*) "memory allocation failed for QREFvis3d!" call abort_physic(subname,'allocation failure for QREFvis3d',1) ENDIF ENDIF IF(.not.ALLOCATED(QREFir3d)) THEN ALLOCATE(QREFir3d(ngrid,nlayer,naerkind), stat=ok) IF (ok/=0) THEN write(*,*) "memory allocation failed for QREFir3d!" call abort_physic(subname,'allocation failure for QREFir3d',1) ENDIF ENDIF ! Effective radius and variance of the aerosols IF(.not.ALLOCATED(reffrad)) THEN allocate(reffrad(ngrid,nlayer,naerkind), stat=ok) IF (ok/=0) THEN write(*,*) "memory allocation failed for reffrad!" call abort_physic(subname,'allocation failure for reffrad',1) ENDIF ENDIF IF(.not.ALLOCATED(nueffrad)) THEN allocate(nueffrad(ngrid,nlayer,naerkind), stat=ok) IF (ok/=0) THEN write(*,*) "memory allocation failed for nueffrad!" call abort_physic(subname,'allocation failure for nueffrad',1) ENDIF ENDIF if (is_master) call system('rm -f surf_vals_long.out') call su_aer_radii(ngrid,nlayer,reffrad,nueffrad) !-------------------------------------------------- ! Set up correlated k !-------------------------------------------------- !this block is now done at firstcall of physiq_mod ! print*, "callcorrk: Correlated-k data base folder:",trim(datadir) ! call getin_p("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 setspi ! Basic infrared properties. ! 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 if(.not.allocated(dtaui)) then ALLOCATE(dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for dtaui!" call abort_physic(subname,'allocation failure for dtaui',1) endif endif if(.not.allocated(dtauv)) then ALLOCATE(dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for dtauv!" call abort_physic(subname,'allocation failure for dtauv',1) endif endif if(.not.allocated(cosbv)) then ALLOCATE(cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for cosbv!" call abort_physic(subname,'allocation failure for cobsv',1) endif endif if(.not.allocated(cosbi)) then ALLOCATE(cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for cosbi!" call abort_physic(subname,'allocation failure for cobsi',1) endif endif if(.not.allocated(wbari)) then ALLOCATE(wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for wbari!" call abort_physic(subname,'allocation failure for wbari',1) endif endif if(.not.allocated(wbarv)) then ALLOCATE(wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for wbarv!" call abort_physic(subname,'allocation failure for wbarv',1) endif endif if(.not.allocated(tauv)) then ALLOCATE(tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for tauv!" call abort_physic(subname,'allocation failure for tauv',1) endif endif if(.not.allocated(taucumv)) then ALLOCATE(taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for taucumv!" call abort_physic(subname,'allocation failure for taucumv',1) endif endif if(.not.allocated(taucumi)) then ALLOCATE(taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for taucumi!" call abort_physic(subname,'allocation failure for taucumi',1) endif endif if(.not.allocated(taugsurf)) then ALLOCATE(taugsurf(L_NSPECTV,L_NGAUSS-1), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for taugsurf!" call abort_physic(subname,'allocation failure for taugsurf',1) endif endif if(.not.allocated(taugsurfi)) then ALLOCATE(taugsurfi(L_NSPECTI,L_NGAUSS-1), stat=ok) if (ok/=0) then write(*,*) "memory allocation failed for taugsurfi!" call abort_physic(subname,'allocation failure for taugsurfi',1) endif endif if(varfixed .and. (generic_condensation .or. methane .or. carbox))then write(*,*) "Deep generic tracer vapor mixing ratio ? (no effect if negative) " qvap_deep=-1. ! default value call getin_p("qvap_deep",qvap_deep) write(*,*) " qvap_deep = ",qvap_deep metallicity=0.0 ! default value --- is not used here but necessary to call function Psat_generic call getin_p("metallicity",metallicity) ! --- is not used here but necessary to call function Psat_generic endif end if ! of if (firstcall) !======================================================================= ! I.b Initialization on every call !======================================================================= qxvaer(:,:,:)=0.0 qsvaer(:,:,:)=0.0 gvaer(:,:,:) =0.0 qxiaer(:,:,:)=0.0 qsiaer(:,:,:)=0.0 giaer(:,:,:) =0.0 OLR_nu(:,:) = 0. OSR_nu(:,:) = 0. GSR_nu(:,:) = 0. !-------------------------------------------------- ! Effective radius and variance of the aerosols !-------------------------------------------------- if (aerohaze) then do iaer=1,naerkind if ((iaer.eq.iaero_haze)) then call su_aer_radii(ngrid,nlayer,reffrad(1,1,iaer), & nueffrad(1,1,iaer)) endif end do !iaer=1,naerkind. if (haze_radproffix) then if (haze_radproffix) then call haze_reffrad_fix(ngrid,nlayer,zzlay,& reffrad,nueffrad) endif 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 ! How much light do we get ? do nw=1,L_NSPECTV stel(nw)=stellarf(nw)/(dist_star**2) end do if (aerohaze) then ! Get 3D aerosol optical properties. 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, pt,pq,aerosol, & reffrad,nueffrad,QREFvis3d,QREFir3d, & tau_col,cloudfrac,totcloudfrac,clearsky) 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,ngrid 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 ! Net atmospheric radiative cooling rate from C2H2 (K.s-1): ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~ ! dtlw_hcn_c2h2=0. if (cooling) then CALL cooling_hcn_c2h2(ngrid,nlayer,pplay,& pt,dtlw_hcn_c2h2) endif !----------------------------------------------------------------------- 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. !----------------------------------------------------------------------- ! AF24: for now only consider one aerosol (=haze) if (aerohaze) then 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 do nw=1,L_NSPECTI ! Longwave do l=1,nlayer temp1=QIRsQREF3d(ig,nlayer+1-l,nw,iaer) & *QREFir3d(ig,nlayer+1-l,iaer) temp2=QIRsQREF3d(ig,max(nlayer-l,1),nw,iaer) & *QREFir3d(ig,max(nlayer-l,1),iaer) qxiaer(2*l,nw,iaer) = temp1 qxiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=temp1*omegair3d(ig,nlayer+1-l,nw,iaer) temp2=temp2*omegair3d(ig,max(nlayer-l,1),nw,iaer) qsiaer(2*l,nw,iaer) = temp1 qsiaer(2*l+1,nw,iaer)=(temp1+temp2)/2 temp1=gir3d(ig,nlayer+1-l,nw,iaer) temp2=gir3d(ig,max(nlayer-l,1),nw,iaer) giaer(2*l,nw,iaer) = temp1 giaer(2*l+1,nw,iaer)=(temp1+temp2)/2 end do ! nlayer qxiaer(1,nw,iaer)=qxiaer(2,nw,iaer) qxiaer(2*nlayer+1,nw,iaer)=0. qsiaer(1,nw,iaer)=qsiaer(2,nw,iaer) qsiaer(2*nlayer+1,nw,iaer)=0. giaer(1,nw,iaer)=giaer(2,nw,iaer) giaer(2*nlayer+1,nw,iaer)=0. end do ! L_NSPECTI 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 do nw=1,L_NSPECTI if(qsiaer(k,nw,iaer).gt.1.05*qxiaer(k,nw,iaer))then message='Serious problems with qsvaer values' call abort_physic(subname,message,1) 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 ! L_LEVELS end do ! naerkind end if ! aerohaze !----------------------------------------------------------------------- ! Aerosol optical depths !----------------------------------------------------------------------- if (aerohaze) then 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)) ! 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 else tauaero(:,:)=0 end if ! aerohaze ! Albedo and Emissivity. albi=1-emis(ig) ! Long Wave. 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. !----------------------------------------------------------------------- ! GCS (Generic Condensable Specie) Vapor ! If you have GCS tracers and they are : variable & radiatively active ! ! NC22 !----------------------------------------------------------------------- if (generic_condensation .or. methane .or. carbox) then ! 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) ! ELSE ! For now, only one GCS tracer can be both variable and radiatively active ! If you set two GCS tracers, that are variable and radiatively active, ! the last one in tracer.def will be chosen as the one that will be vadiatively active do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer if(varactive)then i_var=igcm_generic_gas 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) !JL13index qvar(2*l+1) = (pq(ig,nlayer+1-l,i_var)+pq(ig,max(nlayer-l,1),i_var))/2 !JL13index ! Average approximation as for temperature... end do qvar(1)=qvar(2) elseif(varfixed .and. (qvap_deep .ge. 0))then do l=1,nlayer ! Here we will assign fixed water vapour profiles globally. call Psat_generic(pt(ig,l),pplay(ig,l),metallicity,psat,qsat) if (qsat .lt. qvap_deep) then pq_temp(l) = qsat ! fully saturated everywhere else pq_temp(l) = qvap_deep end if 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) else do k=1,L_LEVELS qvar(k) = 1.0D-7 end do end if ! varactive/varfixed endif end do ! do iq=1,nq loop on tracers end if ! if (generic_condensation) !----------------------------------------------------------------------- ! No Water vapor and No GCS (Generic Condensable Specie) vapor !----------------------------------------------------------------------- if (.not. (generic_condensation .or. methane .or. carbox)) then do k=1,L_LEVELS qvar(k) = 1.0D-7 end do end if ! if (.not. generic_condensation) if(.not.kastprof)then ! IMPORTANT: Now convert from kg/kg to mol/mol. do k=1,L_LEVELS if (generic_condensation .or. methane .or. carbox) then do iq=1,nq call generic_tracer_index(nq,iq,igcm_generic_gas,igcm_generic_ice,call_ice_gas_generic) if (call_ice_gas_generic) then ! to call only one time the ice/vap pair of a tracer if(.not. varactive .or. i_var.eq.iq)then epsi_generic=constants_epsi_generic(iq) qvar(k) = qvar(k)/(epsi_generic+qvar(k)*(1.-epsi_generic)) endif endif end do ! do iq=1,nq loop on tracers endif end do end if !----------------------------------------------------------------------- ! kcm mode only ! !----------------------------------------------------------------------- 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*nlayer+1)=muvar(ig,1) ! 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) = 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="Minimum pressure outside of kmatrix bounds" call abort_physic(subname,message,1) endif ! Test for out-of-bounds temperature. ! -- JVO 20 : Also add a sanity test checking that tlevrad is ! within Planck function temperature boundaries, ! which would cause gfluxi/sfluxi to crash. 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 if (tlevrad(k).lt.tplanckmin) then print*,'Minimum temperature is outside the boundaries for' print*,'Planck function integration set in callphys.def, aborting.' print*,"k=",k," tlevrad(k)=",tlevrad(k) print*,"tplanckmin=",tplanckmin message="Minimum temperature outside Planck function bounds - Change tplanckmin in callphys.def" call abort_physic(subname,message,1) else if (tlevrad(k).gt.tplanckmax) then print*,'Maximum temperature is outside the boundaries for' print*,'Planck function integration set in callphys.def, aborting.' print*,"k=",k," tlevrad(k)=",tlevrad(k) print*,"tplanckmax=",tplanckmax message="Maximum temperature outside Planck function bounds - Change tplanckmax in callphys.def" call abort_physic(subname,message,1) 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 !======================================================================= ! ---------------------------------------------------------------- ! Recombine reference corrk tables if needed - Added by JVO, 2020. if (corrk_recombin) then call call_recombin(ig,nlayer,pq(ig,:,:),pplay(ig,:),pt(ig,:),qvar(:),tmid(:),pmid(:)) endif ! ---------------------------------------------------------------- Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed. glat_ig=glat(ig) !----------------------------------------------------------------------- ! Short Wave Part !----------------------------------------------------------------------- if(fract(ig) .ge. 1.0e-4) then ! Only during daylight. if((ngrid.eq.1).and.(global1d))then do nw=1,L_NSPECTV stel_fract(nw)= stel(nw)* 0.25 / acosz ! globally averaged = divide by 4, and we correct for solar zenith angle end do else do nw=1,L_NSPECTV stel_fract(nw)= stel(nw) * fract(ig) 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,nfluxtopv,fluxtopvdn,nfluxoutv_nu,& nfluxgndv_nu,nfluxtopv_nu, & fmnetv,fmnetv_nu,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 fmnetv_nu(l,:)=0.0d0 fluxupv(l)=0.0d0 fluxdnv(l)=0.0d0 end do end if ! Equivalent Albedo Calculation (for OUTPUT). MT2015 if(fract(ig) .ge. 1.0e-4) then ! equivalent albedo makes sense only during daylight. surface_stellar_flux=sum(nfluxgndv_nu(1:L_NSPECTV)) if(surface_stellar_flux .gt. 1.0e-3) then ! equivalent albedo makes sense only if the stellar flux received by the surface is positive. DO nw=1,L_NSPECTV albedo_temp(nw)=albedo(ig,nw)*nfluxgndv_nu(nw) ENDDO albedo_temp(1:L_NSPECTV)=albedo_temp(1:L_NSPECTV)/surface_stellar_flux albedo_equivalent(ig)=sum(albedo_temp(1:L_NSPECTV)) else albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0. endif else albedo_equivalent(ig)=0.0 ! Spectrally Integrated Albedo not defined for non-irradiated grid points. So we arbitrary set the equivalent albedo to 0. endif !----------------------------------------------------------------------- ! Long Wave Part !----------------------------------------------------------------------- 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,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) ! Flux incident at the top of the atmosphere fluxtop_dn(ig)=fluxtopvdn fluxtop_lw(ig) = real(nfluxtopi) fluxabs_sw(ig) = real(-nfluxtopv) fluxsurf_lw(ig) = real(fluxdni(L_NLAYRAD)) fluxsurf_sw(ig) = real(fluxdnv(L_NLAYRAD)) ! Flux absorbed by the surface. By MT2015. fluxsurfabs_sw(ig) = fluxsurf_sw(ig)*(1.-albedo_equivalent(ig)) 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 ! 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)/DWNV(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)) & ! *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) dpp = glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 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(ig,L_NLAYRAD+1-l)=(fmneti(l)-fmneti(l-1)) & ! *glat(ig)/(cpp*scalep*(plevrad(2*l+1)-plevrad(2*l-1))) 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 ! These are values at top of atmosphere ! dtsw(ig,L_NLAYRAD)=(fmnetv(1)-nfluxtopv) & ! *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2))) ! dtlw(ig,L_NLAYRAD)=(fmneti(1)-nfluxtopi) & ! *glat(ig)/(cpp*scalep*(plevrad(3)-plevrad(2))) dpp = g/(cpp*scalep*(plevrad(3)-plevrad(1))) do nw=1,L_NSPECTV dtsw_nu(L_NLAYRAD,nw)= & (fmnetv_nu(1,nw)-nfluxtopv_nu(nw))*dpp end do do nw=1,L_NSPECTI dtlw_nu(L_NLAYRAD,nw)= & (fmneti_nu(1,nw)-nfluxtopi_nu(nw))*dpp end do ! Optical thickness diagnostics (added by JVO) if (diagdtau) then do l=1,L_NLAYRAD do nw=1,L_NSPECTV int_dtauv(ig,l,nw) = 0.0d0 DO k=1,L_NGAUSS ! Output exp(-tau) because gweight ponderates exp and not tau itself int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k) ENDDO enddo do nw=1,L_NSPECTI int_dtaui(ig,l,nw) = 0.0d0 DO k=1,L_NGAUSS ! Output exp(-tau) because gweight ponderates exp and not tau itself int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k) ENDDO enddo enddo endif ! ********************************************************** ! 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 ! ********************************************************** !----------------------------------------------------------------------- end do ! End of big loop over every GCM column. !----------------------------------------------------------------------- !----------------------------------------------------------------------- ! 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) ! USEFUL COMMENT - Do Not Remove. ! ! 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) then IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi ) IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv ) !$OMP BARRIER !$OMP MASTER IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref ) IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref ) IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar ) IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref ) IF( ALLOCATED( gweight ) ) DEALLOCATE( gweight ) !$OMP END MASTER !$OMP BARRIER IF ( ALLOCATED(reffrad)) DEALLOCATE(reffrad) IF ( ALLOCATED(nueffrad)) DEALLOCATE(nueffrad) endif end subroutine callcorrk END MODULE callcorrk_mod