Ignore:
Timestamp:
Oct 12, 2023, 10:30:22 AM (15 months ago)
Author:
slebonnois
Message:

BBT : Update for the titan microphysics and cloud model

Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
4 added
16 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.TITAN/libf/phytitan/calc_rayleigh.F90

    r1648 r3083  
    1414!     Robin Wordsworth (2010)
    1515!     Jeremy Leconte (2012): Added option for variable gas. Improved water rayleigh (Bucholtz 1995).
     16!     Bruno de Batz de Trenquelleon (2023) : Added CH4 rayleigh.
    1617!
    1718!     Called by
    … …  
    2829      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep
    2930      use gases_h
    30       use comcstfi_mod, only: g, mugaz
     31      use comcstfi_mod, only: g, mugaz, pi
    3132
    3233      implicit none
    … …  
    4647      integer icantbewrong
    4748
    48       ! tau0/p0=tau/p (Hansen 1974)
    49       ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4)
    50       ! where delta=n-1 and N0 is an amagat
     49       ! we multiply by scalep here (100) because plev, which is used in optcv,
     50       ! is in units of mBar, so we need to convert.
     51       ! we calculate here TAURAY which is in m2/mBar
     52
     53       ! tau0/p0=tau/p (Hansen 1974 : equation (2.31) )
     54       ! Then calculate tau0 = (8*pi^3*p_1atm)/(3*N0^2) * 4*delta^2/(g*mugaz*lambda^4) (Hansen 1974 : equation (2.29) )
     55       ! where delta=n-1 and N0 is an amagat
     56
     57       ! In exo_k : (8*pi^3)/(3*N0^2) * 4 * delta^2 ------- is written : (24*pi^3)/(N0^2) * (4/9) * delta^2
     58
     59       ! (tauconsti * tauvari) = scalep * cross_section_molecule_exo_k /  ( gravity * molecular_mass )
     60       ! (tauconsti * tauvari) in m2/mBar because of scalep factor
     61       ! cross_section_molecule in m2/molecule
     62       ! molecular_mass is the mass of th considered molecule
    5163
    5264      typeknown=.false.
    5365      do igas=1,ngasmx
    54          if(gfrac(igas,nivref).lt.5.e-2)then
     66         if(gfrac(igas,nivref).lt.1.e-2)then
    5567            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
    56             'as its mixing ratio is less than 0.05.'
     68            'as its mixing ratio is less than 0.01.'
    5769            ! ignore variable gas in Rayleigh calculation
    58             ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
     70            ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
    5971            tauconsti(igas) = 0.0
    6072         else
    … …  
    6678               ! uses n=1.000132 from Optics, Fourth Edition.
    6779               ! was out by a factor 2!
     80            elseif(igas.eq.igas_CH4)then
     81                ! For CH4 we use the formalism of exo_k (developed by Jeremy Leconte)
     82                ! the relation between exo_k formalism and LMDZ formalism is : (tauconsti*tauvari)=scalep*sigma_exok/(gravity*molecular_mass)
     83                ! Values are taken from Caldas (2019), equation 12 & appendix D
     84                ! 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 is a constant
     85                ! 4/9=(2/3)**2 is approximately ((n+1)/(n**2+2))**2
     86                ! 1E24 = (1E6)**4 because wavelength is in micron
     87                ! 16.04*1.6726*1E-27 is CH4 molecular mass
     88                tauconsti(igas) = 24.*pi**3/(1E5/(1.380648813E-23*273.15))**2 * 4./9. * 1E24 * (1/( g *16.04*1.6726*1E-27))  * scalep             
    6889            else
    6990               print*,'Error in calc_rayleigh: Gas species not recognised!'
    … …  
    96117            tauvar=0.0
    97118            do igas=1,ngasmx
    98                if(gfrac(igas,nivref).lt.5.e-2)then
     119               if(gfrac(igas,nivref).lt.1.e-2)then
    99120                  ! ignore variable gas in Rayleigh calculation
    100                   ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
     121                  ! ignore gases of mixing ratio < 0.01 in Rayleigh calculation
    101122                  tauvari(igas)   = 0.0
    102123               else
    … …  
    105126                  elseif(igas.eq.igas_H2)then
    106127                     tauvari(igas) = 1.0/wl**4
     128                  elseif(igas.eq.igas_CH4)then
     129                      tauvari(igas) = (4.6662E-4+4.02E-6/wl**2)**2 / wl**4
    107130                  else
    108131                     print*,'Error in calc_rayleigh: Gas species not recognised!'

  • trunk/LMDZ.TITAN/libf/phytitan/calc_ysat.F90

    r2326 r3083  
    2828  INTEGER, INTENT(IN)                           :: nlon  ! # of grid points in the chunk
    2929  INTEGER, INTENT(IN)                           :: nlay  ! # of vertical layes
    30   REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: press ! Mid-layers pressure (Pa)
     30  REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: press ! Mid-layers pressure (!!! mbar !!!)
    3131  REAL, DIMENSION(nlon,nlay),      INTENT(IN)   :: temp  ! Mid-layers temperature (K)
    3232  REAL, DIMENSION(nlon,nlay,nkim), INTENT(OUT)  :: ysat  ! Saturation profiles (mol/mol)
    … …  
    4949     if(trim(cnames(ic)).eq."CH4") then
    5050
    51         where (temp(:,:).lt.90.65)                                     
    52            ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                  &
    53                 -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                    &
    54                 + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
    55         elsewhere
    56            ysat(:,:,ic) = 10.0**(3.901408e0 - ( ( 154567.02e0 / temp(:,:) - 1598.8512e0 )   &
    57                 / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
    58         endwhere
     51        !where (temp(:,:).lt.90.65)                                     
     52        !   ysat(:,:,ic) = 10.0**(4.42507e0 - ( ( ( 1165560.7e0 / temp(:,:)                  &
     53        !        -  115352.19e0 ) / temp(:,:) + 4055.6016e0 ) / temp(:,:)                    &
     54        !        + 453.92414e0 ) / temp(:,:) ) / press(:,:) * 1013.25e0
     55        !elsewhere
     56        !   ysat(:,:,ic) = 10.0**(3.901408e0 - ( ( 154567.02e0 / temp(:,:) - 1598.8512e0 )   &
     57        !        / temp(:,:) + 437.54809e0 ) / temp(:,:)) / press(:,:) * 1013.25e0
     58        !endwhere
     59
     60        ! Fray and Schmidt (2009)
     61        ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.051e1 - 1.110e3/temp(:,:)              &
     62                       - 4.341e3/temp(:,:)**2 + 1.035e5/temp(:,:)**3 - 7.910e5/temp(:,:)**4)
    5963
    6064        ! Forcing CH4 to 1.4% minimum               
    … …  
    6670!             * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
    6771
    68          ysat(:,:,ic) = 1.E5 * exp(13.4-2536./temp(:,:)) / press(:,:) ! Fray and Schmidt (2009)
     72         ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.340e1 - 2.536e3/temp(:,:)) ! Fray and Schmidt (2009)
    6973
    7074     else if(trim(cnames(ic)).eq."C2H4") then
    … …  
    8791     else if(trim(cnames(ic)).eq."C2H6") then
    8892
    89         where (temp(:,:).lt.90.)
    90 !           ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                     &
     93!        where (temp(:,:).lt.90.)
     94!           ysat(:,:,ic) = 10.0**(10.01e0-1085.0e0/(temp(:,:)-0.561e0) )                    &
    9195!                / press(:,:) * 1013.25e0 / 760.0e0
    92            ysat(:,:,ic) = 1.E5 * exp ( 15.11 - 2207./temp(:,:) - 24110./(temp(:,:)**2)       &
    93                         + 7.744E5/(temp(:,:)**3) - 1.161E7/(temp(:,:)**4)                    &
    94                         + 6.763E7/(temp(:,:)**5) ) / press(:,:) ! Fray and Schmidt (2009)
    95         elsewhere
    96            ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                 &
    97                 * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
    98         endwhere
     96!           ysat(:,:,ic) = 1.E5 * exp ( 15.11 - 2207./temp(:,:) - 24110./(temp(:,:)**2)     &
     97!                        + 7.744E5/(temp(:,:)**3) - 1.161E7/(temp(:,:)**4)                  &
     98!                        + 6.763E7/(temp(:,:)**5) ) / press(:,:) ! Fray and Schmidt (2009)
     99!        elsewhere
     100!           ysat(:,:,ic) = 10.0**(5.9366e0 - 1086.17e0/temp(:,:) + 3.83464e0                &
     101!                * alog10(1.0e3/temp(:,:)) ) / press(:,:)*1013.25e0/760.0
     102!        endwhere
     103         ysat(:,:,ic) = (1000. / press(:,:)) * exp(1.511e1 - 2.207e3/temp(:,:)            &
     104                      - 2.411e4/temp(:,:)**2 + 7.744e5/temp(:,:)**3 - 1.161e7/temp(:,:)**4  &
     105                      + 6.763e7/temp(:,:)**5)
    99106
    100107     else if((trim(cnames(ic)).eq."CH3CCH") .or. (trim(cnames(ic)).eq."CH2CCH2")) then
    … …  
    150157
    151158        !ysat(:,:,ic)= 10.0**(8.6165e0 - 1516.5e0/(temp(:,:) - 26.2e0) ) / press(:,:) * 1013.25e0 / 760.0e0
    152        
    153         ysat(:,:,ic) = 1.E5 * exp ( 13.93 - 3624./temp(:,:) - 1.325E5/(temp(:,:)**2)       &
    154                      + 6.314E6/(temp(:,:)**3) - 1.128E8/(temp(:,:)**4) ) / press(:,:) ! Fray and Schmidt (2009)
     159
     160        ysat(:,:,ic) = (1000. / press(:,:)) * exp (1.393e1 - 3.624e3/temp(:,:)             &
     161                     - 1.325e5/temp(:,:)**2 + 6.314e6/temp(:,:)**3 - 1.128e8/temp(:,:)**4) ! Fray and Schmidt (2009)
    155162
    156163     else if(trim(cnames(ic)).eq."CH3CN")  then

  • trunk/LMDZ.TITAN/libf/phytitan/callcorrk.F90

    r2408 r3083  
    11      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,      &
    2           albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
    3           tsurf,fract,dist_star,                               &
     2          albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev, &
     3          pt,tsurf,fract,dist_star,                            &
    44          dtlw,dtsw,fluxsurf_lw,                               &
    55          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
    66          fluxabs_sw,fluxtop_dn,                               &
    77          OLR_nu,OSR_nu,                                       &
    8           int_dtaui,int_dtauv,                                 &
     8          int_dtaui,int_dtauv,popthi,popthv,poptci,poptcv,     &
    99          lastcall)
    1010
    … …  
    1818      use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin,   &
    1919                              strictboundcorrk,specOLR,diagdtau,        &
    20                               tplanckmin,tplanckmax
     20                              tplanckmin,tplanckmax,callclouds,Fcloudy
    2121      use geometry_mod, only: latitude
    2222
    … …  
    3636!     Robin Wordsworth (2009)
    3737!     Jan Vatant d'Ollone (2018) -> corrk recombining case
     38!     Bruno de Batz de Trenquelléon (2023) --> Optics of Titan's haze and coulds
    3839!
    3940!==================================================================
    … …  
    5859      REAL,INTENT(IN) :: pplev(ngrid,nlayer+1)     ! Inter-layer pressure (Pa).
    5960      REAL,INTENT(IN) :: pplay(ngrid,nlayer)       ! Mid-layer pressure (Pa).
     61      REAL,INTENT(IN) :: zzlev(ngrid,nlayer+1)     ! Altitude at the atmospheric layer boundaries (ref : local surf).
    6062      REAL,INTENT(IN) :: pt(ngrid,nlayer)          ! Air temperature (K).
    6163      REAL,INTENT(IN) :: tsurf(ngrid)              ! Surface temperature (K).
    … …  
    7678      REAL,INTENT(OUT) :: OSR_nu(ngrid,L_NSPECTV)        ! Outgoing SW radition in each band (Normalized to the band width (W/m2/cm-1).
    7779      REAL,INTENT(OUT) :: albedo_equivalent(ngrid)       ! Spectrally Integrated Albedo. For Diagnostic. By MT2015
    78       REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! VI optical thickness of layers within narrowbands for diags ().
    79       REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! IR optical thickness of layers within narrowbands for diags ().
     80      REAL,INTENT(OUT) :: int_dtaui(ngrid,nlayer,L_NSPECTI) ! IR optical thickness of layers within narrowbands for diags ().
     81      REAL,INTENT(OUT) :: int_dtauv(ngrid,nlayer,L_NSPECTV) ! VI optical thickness of layers within narrowbands for diags ().
     82      REAL,INTENT(OUT) :: popthv(ngrid,nlayer,L_NSPECTV,5)  ! VI optical properties of haze within narrowbands (dtau,tau,k,w,g).
     83      REAL,INTENT(OUT) :: poptcv(ngrid,nlayer,L_NSPECTV,5)  ! VI optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
     84      REAL,INTENT(OUT) :: popthi(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze within narrowbands (dtau,tau,k,w,g).
     85      REAL,INTENT(OUT) :: poptci(ngrid,nlayer,L_NSPECTI,5)  ! IR optical properties of haze and clouds within narrowbands (dtau,tau,k,w,g).
     86     
    8087     
    8188     
    … …  
    9097      ! Optical values for the optci/cv subroutines
    9198      REAL*8 stel(L_NSPECTV),stel_fract(L_NSPECTV)
     99     
    92100      REAL*8 dtaui(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     101      REAL*8 dtaui_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     102      REAL*8 dtaui_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    93103      REAL*8 dtauv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     104      REAL*8 dtauv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     105      REAL*8 dtauv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     106     
    94107      REAL*8 cosbv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     108      REAL*8 cosbv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     109      REAL*8 cosbv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    95110      REAL*8 cosbi(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     111      REAL*8 cosbi_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     112      REAL*8 cosbi_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     113     
    96114      REAL*8 wbari(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     115      REAL*8 wbari_cc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     116      REAL*8 wbari_dc(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    97117      REAL*8 wbarv(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     118      REAL*8 wbarv_cc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     119      REAL*8 wbarv_dc(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     120     
    98121      REAL*8 tauv(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     122      REAL*8 tauv_cc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     123      REAL*8 tauv_dc(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     124     
    99125      REAL*8 taucumv(L_LEVELS,L_NSPECTV,L_NGAUSS)
     126      REAL*8 taucumv_cc(L_LEVELS,L_NSPECTV,L_NGAUSS)
     127      REAL*8 taucumv_dc(L_LEVELS,L_NSPECTV,L_NGAUSS)
    100128      REAL*8 taucumi(L_LEVELS,L_NSPECTI,L_NGAUSS)
     129      REAL*8 taucumi_cc(L_LEVELS,L_NSPECTI,L_NGAUSS)
     130      REAL*8 taucumi_dc(L_LEVELS,L_NSPECTI,L_NGAUSS)
     131     
     132      REAL*8 popt_hazev(L_LEVELS,L_NSPECTV,3)
     133      REAL*8 popt_hazev_cc(L_LEVELS,L_NSPECTV,3)
     134      REAL*8 popt_hazev_dc(L_LEVELS,L_NSPECTV,3)
     135      REAL*8 popt_hazei(L_LEVELS,L_NSPECTI,3)
     136      REAL*8 popt_hazei_cc(L_LEVELS,L_NSPECTI,3)
     137      REAL*8 popt_hazei_dc(L_LEVELS,L_NSPECTI,3)
     138
     139      REAL*8 popt_cloudsv(L_LEVELS,L_NSPECTV,3)
     140      REAL*8 popt_cloudsv_cc(L_LEVELS,L_NSPECTV,3)
     141      REAL*8 popt_cloudsv_dc(L_LEVELS,L_NSPECTV,3)
     142      REAL*8 popt_cloudsi(L_LEVELS,L_NSPECTI,3)
     143      REAL*8 popt_cloudsi_cc(L_LEVELS,L_NSPECTI,3)
     144      REAL*8 popt_cloudsi_dc(L_LEVELS,L_NSPECTI,3)
    101145
    102146      REAL*8 nfluxtopv,nfluxtopi,nfluxtop,fluxtopvdn
    … …  
    110154      REAL*8 albv(L_NSPECTV)                         ! Spectral Visible Albedo.
    111155
    112       INTEGER ig,l,k,nw,iq,ip,ilay,it
     156      INTEGER ig,l,k,nw,iq,ip,ilay,it,lev2lay,cdcolumn,ng
    113157     
    114158      LOGICAL found
    … …  
    449493!-----------------------------------------------------------------------
    450494
     495         ! Clear column :
     496         cdcolumn = 0
     497         call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
     498              dtauv_cc,tauv_cc,taucumv_cc,wbarv_cc,cosbv_cc,tauray,taugsurf,seashazefact,&
     499              popt_hazev_cc,popt_cloudsv_cc,cdcolumn)
     500         ! Dark column :
     501         cdcolumn = 1
     502         call optcv(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tmid,pmid,    &
     503              dtauv_dc,tauv_dc,taucumv_dc,wbarv_dc,cosbv_dc,tauray,taugsurf,seashazefact,&
     504              popt_hazev_dc,popt_cloudsv_dc,cdcolumn)
     505         
     506         ! Mean opacity, ssa and asf :
     507         where ((exp(-dtauv_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtauv_dc(:,:,:)).ge.1.d-40))
     508            dtauv(:,:,:) = -log((1-Fcloudy)*exp(-dtauv_cc(:,:,:)) + Fcloudy*exp(-dtauv_dc(:,:,:)))
     509         elsewhere
     510            dtauv(:,:,:) = dtauv_dc(:,:,:) ! No need to average...
     511         endwhere
     512         do ng = 1, L_NGAUSS
     513            do nw = 1, L_NSPECTV
     514               taucumv(1,nw,ng) = 0.0d0
     515               do k = 2, L_LEVELS
     516                  if ((exp(-taucumv_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumv_dc(k,nw,ng)).ge.1.d-40)) then
     517                     taucumv(k,nw,ng) = taucumv(k-1,nw,ng) - log((1-Fcloudy)*exp(-taucumv_cc(k,nw,ng)) + Fcloudy*exp(-taucumv_dc(k,nw,ng)))
     518                  else
     519                     taucumv(k,nw,ng) = taucumv(k-1,nw,ng) + taucumv_dc(k,nw,ng) ! No need to average...
     520                  end if
     521               end do
     522               do l = 1, L_NLAYRAD
     523                  tauv(l,nw,ng) = taucumv(2*l,nw,ng)
     524               end do
     525               tauv(l,nw,ng) = taucumv(2*L_NLAYRAD+1,nw,ng)
     526            end do
     527         end do
     528         wbarv = (dtauv_cc*wbarv_cc + dtauv_dc*wbarv_dc) / (dtauv_cc + dtauv_dc)
     529         cosbv = (dtauv_cc*wbarv_cc*cosbv_cc + dtauv_dc*wbarv_dc*cosbv_dc) / (dtauv_cc*wbarv_cc + dtauv_dc*wbarv_dc)
     530         
     531         ! Diagnostics for haze :
     532         where (popt_hazev_cc(:,:,1) .lt. 1.d-30)
     533            popt_hazev_cc(:,:,1) = 1.d-30
     534         endwhere
     535         where (popt_hazev_dc(:,:,1) .lt. 1.d-30)
     536            popt_hazev_dc(:,:,1) = 1.d-30
     537         endwhere
     538         popt_hazev(:,:,1) = -log( (1-Fcloudy)*exp(-popt_hazev_cc(:,:,1)) + Fcloudy*exp(-popt_hazev_dc(:,:,1)) )
     539         popt_hazev(:,:,2) = ((popt_hazev_cc(:,:,1)*popt_hazev_cc(:,:,2) + popt_hazev_dc(:,:,1)*popt_hazev_dc(:,:,2)) &
     540                             / (popt_hazev_cc(:,:,1) + popt_hazev_dc(:,:,1)))
     541         where ((popt_hazev_cc(:,:,2).gt.0.0D0) .or. (popt_hazev_dc(:,:,2).gt.0.0D0))
     542            popt_hazev(:,:,3) = (popt_hazev_cc(:,:,1)*popt_hazev_cc(:,:,2)*popt_hazev_cc(:,:,3) + popt_hazev_dc(:,:,1)*popt_hazev_dc(:,:,2)*popt_hazev_dc(:,:,3)) &
     543                             / (popt_hazev_cc(:,:,1)*popt_hazev_cc(:,:,2) + popt_hazev_dc(:,:,1)*popt_hazev_dc(:,:,2))
     544         elsewhere
     545            popt_hazev(:,:,3) = 0.0D0
     546         endwhere
     547         ! Diagnostics for clouds :
     548         where (popt_cloudsv_cc(:,:,1) .lt. 1.d-30)
     549            popt_cloudsv_cc(:,:,1) = 1.d-30
     550         endwhere
     551         where (popt_cloudsv_dc(:,:,1) .lt. 1.d-30)
     552            popt_cloudsv_dc(:,:,1) = 1.d-30
     553         endwhere
     554         popt_cloudsv(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsv_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsv_dc(:,:,1)) )
     555         popt_cloudsv(:,:,2) = ((popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)) &
     556                             / (popt_cloudsv_cc(:,:,1) + popt_cloudsv_dc(:,:,1)))
     557         where ((popt_cloudsv_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsv_dc(:,:,2).gt.0.0D0))
     558            popt_cloudsv(:,:,3) = (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2)*popt_cloudsv_cc(:,:,3) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2)*popt_cloudsv_dc(:,:,3)) &
     559                             / (popt_cloudsv_cc(:,:,1)*popt_cloudsv_cc(:,:,2) + popt_cloudsv_dc(:,:,1)*popt_cloudsv_dc(:,:,2))
     560         elsewhere
     561            popt_cloudsv(:,:,3) = 0.0D0
     562         endwhere
     563         
    451564         if(fract(ig) .ge. 1.0e-4) then ! Only during daylight.
    452565            if((ngrid.eq.1).and.(global1d))then
    … …  
    459572               end do
    460573            endif
    461            
    462             call optcv(pqmo(ig,:,:),nlayer,plevrad,tmid,pmid,   &
    463                  dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf,seashazefact)
    464 
    465574            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
    466575                 acosz,stel_fract,                                 &
    … …  
    469578
    470579         else ! During the night, fluxes = 0.
    471             nfluxtopv       = 0.0d0
    472             fluxtopvdn      = 0.0d0
    473             nfluxoutv_nu(:) = 0.0d0
    474             nfluxgndv_nu(:) = 0.0d0
     580            nfluxtopv        = 0.0d0
     581            fluxtopvdn       = 0.0d0
     582            nfluxoutv_nu(:)  = 0.0d0
     583            nfluxgndv_nu(:)  = 0.0d0
    475584            do l=1,L_NLAYRAD
    476585               fmnetv(l)=0.0d0
    … …  
    502611!-----------------------------------------------------------------------
    503612
    504          call optci(pqmo(ig,:,:),nlayer,plevrad,tlevrad,tmid,pmid,   &
    505               dtaui,taucumi,cosbi,wbari,taugsurfi,seashazefact)
     613
     614         ! Clear column :
     615         cdcolumn = 0
     616         call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
     617              dtaui_cc,taucumi_cc,cosbi_cc,wbari_cc,taugsurfi,seashazefact,   &
     618              popt_hazei_cc,popt_cloudsi_cc,cdcolumn)
     619         ! Dark column :
     620         cdcolumn = 1
     621         call optci(pqmo(ig,:,:),nlayer,zzlev(ig,:),plevrad,tlevrad,tmid,pmid,&
     622              dtaui_dc,taucumi_dc,cosbi_dc,wbari_dc,taugsurfi,seashazefact,   &
     623              popt_hazei_dc,popt_cloudsi_dc,cdcolumn)
     624
     625         ! Mean opacity, ssa and asf :
     626         where ((exp(-dtaui_cc(:,:,:)).ge.1.d-40) .and. (exp(-dtaui_dc(:,:,:)).ge.1.d-40))
     627            dtaui(:,:,:) = -log((1-Fcloudy)*exp(-dtaui_cc(:,:,:)) + Fcloudy*exp(-dtaui_dc(:,:,:)))
     628         elsewhere
     629            dtaui(:,:,:) = dtaui_dc(:,:,:) ! No need to average...
     630         endwhere
     631         do ng = 1, L_NGAUSS
     632            do nw = 1, L_NSPECTI
     633               taucumi(1,nw,ng) = 0.0d0
     634               do k = 2, L_LEVELS
     635                  if ((exp(-taucumi_cc(k,nw,ng)).ge.1.d-40) .and. (exp(-taucumi_dc(k,nw,ng)).ge.1.d-40)) then
     636                     taucumi(k,nw,ng) = taucumi(k-1,nw,ng) - log((1-Fcloudy)*exp(-taucumi_cc(k,nw,ng)) + Fcloudy*exp(-taucumi_dc(k,nw,ng)))
     637                  else
     638                     taucumi(k,nw,ng) = taucumi(k-1,nw,ng) + taucumi_dc(k,nw,ng) ! No need to average...
     639                  end if
     640               end do
     641            end do
     642         end do
     643         where ((wbari_cc.gt.0.0D0) .or. (wbari_dc.gt.0.0D0))
     644            wbari = (dtaui_cc*wbari_cc + dtaui_dc*wbari_dc) / (dtaui_cc + dtaui_dc)
     645            cosbi = (dtaui_cc*wbari_cc*cosbi_cc + dtaui_dc*wbari_dc*cosbi_dc) / (dtaui_cc*wbari_cc + dtaui_dc*wbari_dc)
     646         elsewhere
     647            wbari = 0.0D0
     648            cosbi = 0.0D0
     649         endwhere
     650
     651         ! Diagnostics for haze :
     652         where (popt_hazei_cc(:,:,1) .lt. 1.d-30)
     653            popt_hazei_cc(:,:,1) = 1.d-30
     654         endwhere
     655         where (popt_hazei_dc(:,:,1) .lt. 1.d-30)
     656            popt_hazei_dc(:,:,1) = 1.d-30
     657         endwhere
     658         popt_hazei(:,:,1) = -log( (1-Fcloudy)*exp(-popt_hazei_cc(:,:,1)) + Fcloudy*exp(-popt_hazei_dc(:,:,1)) )
     659         popt_hazei(:,:,2) = (popt_hazei_cc(:,:,1)*popt_hazei_cc(:,:,2) + popt_hazei_dc(:,:,1)*popt_hazei_dc(:,:,2)) &
     660                             / (popt_hazei_cc(:,:,1) + popt_hazei_dc(:,:,1))
     661         where ((popt_hazei_cc(:,:,2).gt.0.0D0) .or. (popt_hazei_dc(:,:,2).gt.0.0D0))
     662            popt_hazei(:,:,3) = (popt_hazei_cc(:,:,1)*popt_hazei_cc(:,:,2)*popt_hazei_cc(:,:,3) + popt_hazei_dc(:,:,1)*popt_hazei_dc(:,:,2)*popt_hazei_dc(:,:,3)) &
     663                             / (popt_hazei_cc(:,:,1)*popt_hazei_cc(:,:,2) + popt_hazei_dc(:,:,1)*popt_hazei_dc(:,:,2))
     664         elsewhere
     665            popt_hazei(:,:,3) = 0.0D0
     666         endwhere
     667         ! Diagnostics for clouds :
     668         where (popt_cloudsi_cc(:,:,1) .lt. 1.d-30)
     669            popt_cloudsi_cc(:,:,1) = 1.d-30
     670         endwhere
     671         where (popt_cloudsi_dc(:,:,1) .lt. 1.d-30)
     672            popt_cloudsi_dc(:,:,1) = 1.d-30
     673         endwhere
     674         popt_cloudsi(:,:,1) = -log( (1-Fcloudy)*exp(-popt_cloudsi_cc(:,:,1)) + Fcloudy*exp(-popt_cloudsi_dc(:,:,1)) )
     675         popt_cloudsi(:,:,2) = ((popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)) &
     676                             / (popt_cloudsi_cc(:,:,1) + popt_cloudsi_dc(:,:,1)))
     677         where ((popt_cloudsi_cc(:,:,2).gt.0.0D0) .or. (popt_cloudsi_dc(:,:,2).gt.0.0D0))
     678            popt_cloudsi(:,:,3) = (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2)*popt_cloudsi_cc(:,:,3) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2)*popt_cloudsi_dc(:,:,3)) &
     679                             / (popt_cloudsi_cc(:,:,1)*popt_cloudsi_cc(:,:,2) + popt_cloudsi_dc(:,:,1)*popt_cloudsi_dc(:,:,2))
     680         elsewhere
     681            popt_cloudsi(:,:,3) = 0.0D0
     682         endwhere
    506683
    507684         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
    … …  
    562739      !  Optical thickness diagnostics (added by JVO)
    563740      if (diagdtau) then
    564         do l=1,L_NLAYRAD
    565           do nw=1,L_NSPECTV
    566             int_dtauv(ig,l,nw) = 0.0d0
    567              DO k=1,L_NGAUSS
    568               ! Output exp(-tau) because gweight ponderates exp and not tau itself
    569               int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
    570              ENDDO
    571           enddo
    572           do nw=1,L_NSPECTI
    573            int_dtaui(ig,l,nw) = 0.0d0
    574              DO k=1,L_NGAUSS
    575               ! Output exp(-tau) because gweight ponderates exp and not tau itself
    576               int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
    577              ENDDO
    578           enddo
    579         enddo
    580       endif       
     741         do l=1,L_NLAYRAD
     742            do nw=1,L_NSPECTV
     743               int_dtauv(ig,l,nw) = 0.0d0
     744               DO k=1,L_NGAUSS
     745               ! Output exp(-tau) because gweight ponderates exp and not tau itself
     746               int_dtauv(ig,l,nw)= int_dtauv(ig,l,nw) + exp(-dtauv(l,nw,k))*gweight(k)
     747               ENDDO
     748            enddo
     749            do nw=1,L_NSPECTI
     750            int_dtaui(ig,l,nw) = 0.0d0
     751               DO k=1,L_NGAUSS
     752               ! Output exp(-tau) because gweight ponderates exp and not tau itself
     753               int_dtaui(ig,l,nw)= int_dtaui(ig,l,nw) + exp(-dtaui(l,nw,k))*gweight(k)
     754               ENDDO
     755            enddo
     756         enddo
     757      endif
     758     
     759     
     760      ! Optical properties of haze and clouds (dtau, tau, k, w, g) :
     761      do l = 2, L_LEVELS, 2
     762         lev2lay = L_NLAYRAD+1 - l/2
     763         
     764         ! Visible :
     765         do nw = 1, L_NSPECTV
     766            popthv(ig,lev2lay,nw,:) = 0.0d0
     767            poptcv(ig,lev2lay,nw,:) = 0.0d0
     768            ! Optical thickness (dtau) :
     769            popthv(ig,lev2lay,nw,1) = (popt_hazev(l,nw,1) + popt_hazev(l+1,nw,1)) / 2.0
     770            poptcv(ig,lev2lay,nw,1) = (popt_cloudsv(l,nw,1) + popt_cloudsv(l+1,nw,1)) / 2.0
     771            ! Opacity (tau) :
     772            do k = L_NLAYRAD, lev2lay, -1
     773               popthv(ig,lev2lay,nw,2) = popthv(ig,lev2lay,nw,2) + popthv(ig,k,nw,1)
     774               poptcv(ig,lev2lay,nw,2) = poptcv(ig,lev2lay,nw,2) + poptcv(ig,k,nw,1)
     775            enddo
     776            ! Extinction (k) :
     777            popthv(ig,lev2lay,nw,3) = popthv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     778            poptcv(ig,lev2lay,nw,3) = poptcv(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     779            ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
     780            popthv(ig,lev2lay,nw,4) = (popt_hazev(l,nw,2) + popt_hazev(l+1,nw,2)) / 2.0
     781            poptcv(ig,lev2lay,nw,4) = (popt_cloudsv(l,nw,2) + popt_cloudsv(l+1,nw,2)) / 2.0
     782            popthv(ig,lev2lay,nw,5) = (popt_hazev(l,nw,3) + popt_hazev(l+1,nw,3)) / 2.0
     783            poptcv(ig,lev2lay,nw,5) = (popt_cloudsv(l,nw,3) + popt_cloudsv(l+1,nw,3)) / 2.0         
     784         enddo
     785
     786         ! Infra-Red
     787         do nw=1,L_NSPECTI
     788            popthi(ig,lev2lay,nw,:) = 0.0d0
     789            poptci(ig,lev2lay,nw,:) = 0.0d0
     790            ! Optical thickness (dtau) :
     791            popthi(ig,lev2lay,nw,1) = (popt_hazei(l,nw,1) + popt_hazei(l+1,nw,1)) / 2.0
     792            poptci(ig,lev2lay,nw,1) = (popt_cloudsi(l,nw,1) + popt_cloudsi(l+1,nw,1)) / 2.0
     793            ! Opacity (tau) :
     794            do k = L_NLAYRAD, lev2lay, -1
     795               popthi(ig,lev2lay,nw,2) = popthi(ig,lev2lay,nw,2) + popthi(ig,k,nw,1)
     796               poptci(ig,lev2lay,nw,2) = poptci(ig,lev2lay,nw,2) + poptci(ig,k,nw,1)
     797            enddo
     798            ! Extinction (k) :
     799            popthi(ig,lev2lay,nw,3) = popthi(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     800            poptci(ig,lev2lay,nw,3) = poptci(ig,lev2lay,nw,1) / (zzlev(ig,lev2lay+1) - zzlev(ig,lev2lay))
     801            ! Simple Scattering Albedo (w) and Asymmetry Parameter (g) :
     802            popthi(ig,lev2lay,nw,4) = (popt_hazei(l,nw,2) + popt_hazei(l+1,nw,2)) / 2.0
     803            poptci(ig,lev2lay,nw,4) = (popt_cloudsi(l,nw,2) + popt_cloudsi(l+1,nw,2)) / 2.0
     804            popthi(ig,lev2lay,nw,5) = (popt_hazei(l,nw,3) + popt_hazei(l+1,nw,3)) / 2.0
     805            poptci(ig,lev2lay,nw,5) = (popt_cloudsi(l,nw,3) + popt_cloudsi(l+1,nw,3)) / 2.0
     806         enddo
     807      enddo
    581808
    582809

  • trunk/LMDZ.TITAN/libf/phytitan/callkeys_mod.F90

    r2793 r3083  
    4545      logical,save :: eff_gz
    4646!$OMP THREADPRIVATE(eff_gz)
    47      
     47
     48      logical,save :: nudging_u
     49!$OMP THREADPRIVATE(nudging_u)
     50      real,save    :: nudging_dt
     51!$OMP THREADPRIVATE(nudging_dt)
     52
     53      logical,save :: opt4clouds
     54!$OMP THREADPRIVATE(opt4clouds)
     55      real,save    :: Fcloudy
     56      real,save    :: Fssa
     57      real,save    :: FHVIS
     58      real,save    :: FHIR
     59!$OMP THREADPRIVATE(Fcloudy,Fssa,FHVIS,FHIR)
     60      logical,save :: resCH4_inf
     61!$OMP THREADPRIVATE(resCH4_inf)
     62       
    4863      integer,save :: ichim
    4964!$OMP THREADPRIVATE(ichim)

  • trunk/LMDZ.TITAN/libf/phytitan/calmufi.F90

    r2369 r3083  
    1818  !! done elsewhere.
    1919  !!
    20   !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017
     20  !! Authors : J.Burgalat, J.Vatant d'Ollone - 2017, B. de Batz de Trenquelléon (2023)
    2121  !!
    2222  USE MMP_GCM
    … …  
    6161
    6262  REAL(kind=8), DIMENSION(:,:), ALLOCATABLE ::  int2ext !! (\(m^{-2}\)).
     63  REAL(kind=8), DIMENSION(:), ALLOCATABLE   :: tmp
    6364  TYPE(error) :: err
    6465
    … …  
    131132      m3n(:) = zq(ilon,:,6) * int2ext(ilon,:)
    132133      do i=1,nices
    133         m3i(:,nices) = zq(ilon,:,6+i) * int2ext(ilon,:)
    134         gazs(:,i)    = zq(ilon,:,ices_indx(i)) / rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
     134        m3i(:,i)  = zq(ilon,:,ices_indx(i)) * int2ext(ilon,:)
     135        gazs(:,i) = zq(ilon,:,gazs_indx(i)) / rat_mmol(gazs_indx(i)) ! For gazs we work on the full tracer array !!
    135136        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
    136137      enddo
    137138    endif
    138139
     140    ! Hackin the pressure level
     141    tmp = plev(ilon,:)
     142    if (tmp(nlay+1) == 0.0) then
     143      tmp(nlay+1) = 2*tmp(nlay) - tmp(nlay-1)
     144    endif
     145
    139146    ! Initialize YAMMS atmospheric column
    140     err = mm_column_init(plev(ilon,:),zlev(ilon,:),play(ilon,:),zlay(ilon,:),temp(ilon,:)) ; IF (err /= 0) call abort_program(err)
     147    err = mm_column_init(tmp,zlev(ilon,:),play(ilon,:),zlay(ilon,:),temp(ilon,:)) ; IF (err /= 0) call abort_program(err)
    141148    ! Initialize YAMMS aerosols moments column
    142149    err = mm_aerosols_init(m0as,m3as,m0af,m3af) ; IF (err /= 0) call abort_program(err)
    … …  
    148155
    149156    ! initializes tendencies:
    150     !dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp
    151     !dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp
    152 
    153     dm0as(:) = 0.D0 ; dm3as(:) = 0.D0 ; dm0af(:) = 0.D0 ; dm3af(:) = 0.D0
    154     dm0n(:) = 0.D0 ; dm3n(:) = 0.D0 ; dm3i(:,:) = 0.D0 ; dgazs(:,:) = 0.D0
     157    dm0as(:) = 0._mm_wp ; dm3as(:) = 0._mm_wp ; dm0af(:) = 0._mm_wp ; dm3af(:) = 0._mm_wp
     158    dm0n(:) = 0._mm_wp ; dm3n(:) = 0._mm_wp ; dm3i(:,:) = 0._mm_wp ; dgazs(:,:) = 0._mm_wp
     159   
    155160    ! call microphysics
    156 
    157161    IF (callclouds) THEN ! call clouds
    158162      IF(.NOT.mm_muphys(dm0as,dm3as,dm0af,dm3af,dm0n,dm3n,dm3i,dgazs)) &
    … …  
    163167    ENDIF
    164168    ! save diags (if no clouds, relevant arrays will be set to 0 !)
    165     call mm_diagnostics(mmd_aer_prec(ilon),mmd_aer_s_flux(ilon,:),mmd_aer_f_flux(ilon,:),  &
    166                         mmd_ccn_prec(ilon),mmd_ccn_flux(ilon,:), mmd_ice_prec(ilon,:),   &
     169    call mm_diagnostics(mmd_aer_prec(ilon),mmd_aer_s_w(ilon,:),mmd_aer_f_w(ilon,:),mmd_aer_s_flux(ilon,:),mmd_aer_f_flux(ilon,:),  &
     170                        mmd_ccn_prec(ilon),mmd_ccn_w(ilon,:),mmd_ccn_flux(ilon,:),mmd_ice_prec(ilon,:),   &
    167171                        mmd_ice_fluxes(ilon,:,:),mmd_gazs_sat(ilon,:,:))
    168172    call mm_get_radii(mmd_rc_sph(ilon,:),mmd_rc_fra(ilon,:),mmd_rc_cld(ilon,:))
    … …  
    180184      zdq(ilon,:,6) = dm3n(:) / int2ext(ilon,:)
    181185      do i=1,nices
    182         zdq(ilon,:,6+i) = dm3i(:,nices) / int2ext(ilon,:)
    183         zdq(ilon,:,ices_indx(i)) = dgazs(:,i) * rat_mmol(ices_indx(i)) ! For gazs we work on the full tracer array !!
     186        zdq(ilon,:,ices_indx(i)) = dm3i(:,i) / int2ext(ilon,:)
     187        zdq(ilon,:,gazs_indx(i)) = dgazs(:,i) * rat_mmol(gazs_indx(i)) ! For gazs we work on the full tracer array !!
    184188        ! We use the molar mass ratio from GCM in case there is discrepancy with the mm one
    185189      enddo

  • trunk/LMDZ.TITAN/libf/phytitan/datafile_mod.F90

    r2242 r3083  
    2727      ! Default file for coupled microphysics optical properties
    2828      ! Set in physiq_mod
    29       character(LEN=100),save :: haze_opt_file ='datagcm/HAZE_OPT_23x23.DAT'
     29      character(LEN=100),save :: haze_opt_file ='datagcm/optical_tables/HAZE_OPT_23x23.DAT'
    3030!$OMP THREADPRIVATE(haze_opt_file)
     31
     32      ! Default file for nudging of zonal wind
     33      ! Set in physiq_mod
     34      character(LEN=100),save :: nudging_file ='datagcm/SuperNudging.dat'
     35!$OMP THREADPRIVATE(nudging_file)
    3136
    3237      end module datafile_mod

  • trunk/LMDZ.TITAN/libf/phytitan/disr_haze.F90

    r2040 r3083  
    4545! wl_PL in nm
    4646! press_PL in Pa
    47    open(11,file=TRIM(datadir)//'/hazetable_PL_original.dat',status="old",iostat=ierr)
     47   open(11,file=TRIM(datadir)//'/optical_tables/hazetable_PL_original.dat',status="old",iostat=ierr)
    4848   read(11,*) dummy,wl_PL
    4949   do i=1,nblev_PL

  • trunk/LMDZ.TITAN/libf/phytitan/inifis_mod.F90

    r2366 r3083  
    194194     call getin_p("eff_gz", eff_gz)
    195195     write(*,*) "eff_gz = ", eff_gz
     196
     197     write(*,*) "Nudging of zonal wind to keep the super-rotation"
     198     nudging_u = .true.
     199     call getin_p("nudging_u", nudging_u)
     200     write(*,*) "nudging_u = ", nudging_u
     201
     202     write(*,*) "Nudging relaxation time (here ~ 1 titanian year in secondes)"
     203     nudging_dt = 9.3e8
     204     call getin_p("nudging_dt", nudging_dt)
     205     write(*,*) "nudging_dt = ", nudging_dt
     206
     207     write(*,*) "New optics for clouds"
     208     opt4clouds = .true.
     209     call getin_p("opt4clouds", opt4clouds)
     210     write(*,*) "opt4clouds = ", opt4clouds
     211
     212     write(*,*) "Cloudy fraction in the Clear / Dark column method"
     213     Fcloudy = 1.0
     214     call getin_p("Fcloudy", Fcloudy)
     215     write(*,*) "Fcloudy = ", Fcloudy
     216
     217     write(*,*) "Adjustment of single scattering albedo for cloudy dropplet"
     218     Fssa = 1.0
     219     call getin_p("Fssa", Fssa)
     220     write(*,*) "Fssa = ", Fssa
     221     
     222     write(*,*) "Adjustment of aerosol properties in the VIS"
     223     FHVIS = 1.0
     224     call getin_p("FHVIS", FHVIS)
     225     write(*,*) "FHVIS = ", FHVIS
     226
     227     write(*,*) "Adjustment of aerosol properties in the IR"
     228     FHIR = 1.0
     229     call getin_p("FHIR", FHIR)
     230     write(*,*) "FHIR = ", FHIR
     231
     232     write(*,*) "Infinite tank of CH4 ?"
     233     resCH4_inf = .false.
     234     call getin_p("resCH4_inf", resCH4_inf)
     235     write(*,*) "resCH4_inf = ", resCH4_inf
    196236     
    197237     ! sanity check warning

  • trunk/LMDZ.TITAN/libf/phytitan/inimufi.F90

    r2369 r3083  
    2323  !     -------
    2424  !     J. Vatant d'Ollone, J.Burgalat, S.Lebonnois (09/2017)
     25  !     B. de Batz de Trenquelléon (02/2023)
    2526  !
    2627  !============================================================================
    … …  
    7374       ENDIF
    7475     ENDDO
     76     
     77     nice = mm_nesp
     78
    7579     ALLOCATE(ices_indx(mm_nesp))
     80     ALLOCATE(gazs_indx(mm_nesp))
    7681     ices_indx(:) = -1
     82     gazs_indx(:) = -1
    7783     DO i=1,mm_nesp
    7884       idx = indexOfTracer("mu_m3"//TRIM(mm_spcname(i)),.false.)
    … …  
    8389         ices_indx(i) = idx
    8490       ENDIF
     91
     92       idx = indexOfTracer(TRIM(mm_spcname(i)),.false.)
     93       IF (idx <= 0) THEN
     94         WRITE(*,*) "inimufi: '"//TRIM(mm_spcname(i))//"' not found in tracers table."
     95         err = .true.
     96       ELSE
     97         gazs_indx(i) = idx
     98       ENDIF
    8599     ENDDO
     100
    86101     IF (err) THEN
    87102       WRITE(*,*) "You loose.... Try again"
    88103       STOP
    89104     ENDIF
     105     
     106     write(*,*) 'DUMPTRACERS - MICRO_INDX'
     107     call dumpTracers(micro_indx)
     108     write(*,*) 'DUMPTRACERS - ICES_INDX'
     109     call dumpTracers(ices_indx)
     110     write(*,*) 'DUMPTRACERS - GAZS_INDX'
     111     call dumpTracers(gazs_indx)
    90112   ENDIF
    91 
     113   
    92114#endif
    93115

  • trunk/LMDZ.TITAN/libf/phytitan/muphy_diag.F90

    r2793 r3083  
    1414  REAL(kind=8), ALLOCATABLE, DIMENSION(:)     :: mmd_aer_prec   !! Aerosols precipitations (both modes) (m).
    1515  REAL(kind=8), ALLOCATABLE, DIMENSION(:)     :: mmd_ccn_prec   !! CCN precipitations (m).
     16  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mmd_aer_s_w    !! Spherical aerosol settling velocity (\(m.s^{-1}\)).
     17  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mmd_aer_f_w    !! Fractal aerosol settling velocity (\(m.s^{-1}\)).
     18  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mmd_ccn_w      !! CCN settling velocity (\(m.s^{-1}\)).
    1619  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mmd_aer_s_flux !! Spherical aerosol mass flux (\(kg.m^{-2}.s^{-1}\)).
    1720  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mmd_aer_f_flux !! Fractal aerosol mass flux (\(kg.m^{-2}.s^{-1}\)).
    … …  
    2427  REAL(kind=8), ALLOCATABLE, DIMENSION(:,:)   :: mmd_rc_cld     !! Cloud drop radius (m).
    2528
    26   !$OMP THREADPRIVATE(mmd_aer_prec,mmd_ccn_prec,mmd_aer_s_flux,mmd_aer_f_flux,mmd_ccn_flux,mmd_ice_fluxes)
     29  !$OMP THREADPRIVATE(mmd_aer_prec,mmd_ccn_prec,mmd_aer_s_w,mmd_aer_f_w,mmd_ccn_w,mmd_aer_s_flux,mmd_aer_f_flux,mmd_ccn_flux,mmd_ice_fluxes)
    2730  !$OMP THREADPRIVATE(mmd_gazs_sat,mmd_ice_prec,mmd_rc_sph,mmd_rc_fra,mmd_rc_cld)
    2831
    … …  
    3639    ALLOCATE(mmd_aer_prec(ngrid))
    3740    ALLOCATE(mmd_ccn_prec(ngrid))
     41    ALLOCATE(mmd_aer_s_w(ngrid,nlayer))
     42    ALLOCATE(mmd_aer_f_w(ngrid,nlayer))
     43    ALLOCATE(mmd_ccn_w(ngrid,nlayer))
    3844    ALLOCATE(mmd_aer_s_flux(ngrid,nlayer))
    3945    ALLOCATE(mmd_aer_f_flux(ngrid,nlayer))
    … …  
    4854    mmd_aer_prec(:)       = 0d0
    4955    mmd_ccn_prec(:)       = 0d0
     56    mmd_aer_s_w(:,:)      = 0d0
     57    mmd_aer_f_w(:,:)      = 0d0
     58    mmd_ccn_w(:,:)        = 0d0
    5059    mmd_aer_s_flux(:,:)   = 0d0
    5160    mmd_aer_f_flux(:,:)   = 0d0
    … …  
    6473    IF (ALLOCATED(mmd_aer_prec))   DEALLOCATE(mmd_aer_prec)
    6574    IF (ALLOCATED(mmd_ccn_prec))   DEALLOCATE(mmd_ccn_prec)
     75    IF (ALLOCATED(mmd_aer_s_w))    DEALLOCATE(mmd_aer_s_w)
     76    IF (ALLOCATED(mmd_aer_f_w))    DEALLOCATE(mmd_aer_f_w)
     77    IF (ALLOCATED(mmd_ccn_w))      DEALLOCATE(mmd_ccn_w)
    6678    IF (ALLOCATED(mmd_aer_s_flux)) DEALLOCATE(mmd_aer_s_flux)
    6779    IF (ALLOCATED(mmd_aer_f_flux)) DEALLOCATE(mmd_aer_f_flux)

  • trunk/LMDZ.TITAN/libf/phytitan/optci.F90

    r2242 r3083  
    1 subroutine optci(PQMO,NLAY,PLEV,TLEV,TMID,PMID,      &
    2      DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT)
     1subroutine optci(PQMO,NLAY,ZLEV,PLEV,TLEV,TMID,PMID, &
     2     DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT,&
     3     POPT_HAZE,POPT_CLOUDS,CDCOLUMN)
    34
    45  use radinc_h
    … …  
    78  use gases_h
    89  use datafile_mod, only: haze_opt_file
    9   use comcstfi_mod, only: r
    10   use callkeys_mod, only: continuum,graybody,corrk_recombin,               &
    11                           callclouds,callmufi,seashaze,uncoupl_optic_haze
    12   use tracer_h, only : nmicro,nice
     10  use comcstfi_mod, only: pi,r
     11  use callkeys_mod, only: continuum,graybody,corrk_recombin,              &
     12                          callclouds,callmufi,seashaze,uncoupl_optic_haze,&
     13                          opt4clouds,FHIR
     14  use tracer_h, only: nmicro,nice,ices_indx
    1315
    1416  implicit none
    … …  
    3234  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
    3335  !     Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
     36  !     Clean and correction to Titan by B. de Batz de Trenquelléon (2022)
     37  !     New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023)
     38  !     The whole routine needs to be redone...
    3439  !     
    3540  !==================================================================
    … …  
    4146  REAL*8, INTENT(IN)  :: PQMO(nlay,nmicro)  ! Tracers for microphysics optics (X/m2).
    4247  INTEGER, INTENT(IN) :: NLAY               ! Number of pressure layers (for pqmo)
     48  REAL*8, INTENT(IN)  :: ZLEV(NLAY+1)
    4349  REAL*8, INTENT(IN)  :: PLEV(L_LEVELS), TLEV(L_LEVELS)
    4450  REAL*8, INTENT(IN)  :: TMID(L_LEVELS), PMID(L_LEVELS)
    4551  REAL*8, INTENT(IN)  :: SEASHAZEFACT(L_LEVELS)
     52  INTEGER, INTENT(IN) :: CDCOLUMN
    4653 
    4754  REAL*8, INTENT(OUT) :: DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    … …  
    5057  REAL*8, INTENT(OUT) :: WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    5158  REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTI,L_NGAUSS-1)
     59  REAL*8, INTENT(OUT) :: POPT_HAZE(L_LEVELS,L_NSPECTI,3)
     60  REAL*8, INTENT(OUT) :: POPT_CLOUDS(L_LEVELS,L_NSPECTI,3)
    5261  ! ==========================================================
    5362 
    … …  
    94103  real*8  dumwvl
    95104
    96   real*8 m3as,m3af
    97   real*8 dtauaer_s,dtauaer_f
     105  ! Variables for new optics
     106  integer iq, iw, FTYPE, CTYPE
     107  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld
     108  real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
    98109  real*8,save :: rhoaer_s(L_NSPECTI),ssa_s(L_NSPECTI),asf_s(L_NSPECTI)
    99110  real*8,save :: rhoaer_f(L_NSPECTI),ssa_f(L_NSPECTI),asf_f(L_NSPECTI)
    100 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
     111  real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI)
     112  real*8,save :: ssa_cld(L_NSPECTI),asf_cld(L_NSPECTI)
     113!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld)
    101114
    102115  logical,save :: firstcall=.true.
    … …  
    127140  ENDIF
    128141
    129   !=======================================================================
    130   !     Determine the total gas opacity throughout the column, for each
    131   !     spectral interval, NW, and each Gauss point, NG.
    132 
    133   taugsurf(:,:) = 0.0
    134   dpr(:)        = 0.0
    135   lkcoef(:,:)   = 0.0
    136 
    137   do K=2,L_LEVELS
    138  
    139      ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    140      
    141      DPR(k) = PLEV(K)-PLEV(K-1)
    142 
    143      ! if we have continuum opacities, we need dz
     142   !=======================================================================
     143   !     Determine the total gas opacity throughout the column, for each
     144   !     spectral interval, NW, and each Gauss point, NG.
     145
     146   taugsurf(:,:) = 0.0
     147   dpr(:)        = 0.0
     148   lkcoef(:,:)   = 0.0
     149
     150   do K=2,L_LEVELS
     151      ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
     152      DPR(k) = PLEV(K)-PLEV(K-1)
     153
     154      ! if we have continuum opacities, we need dz
    144155      dz(k) = dpr(k)*R*TMID(K)/(gzlat_ig(ilay)*PMID(K))
    145156      U(k)  = Cmk(ilay)*DPR(k)     ! only Cmk line in optci.F
    146      
    147      call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
    148 
    149      do LK=1,4
    150         LKCOEF(K,LK) = LCOEF(LK)
    151      end do
    152   end do                    ! levels
    153 
    154   do NW=1,L_NSPECTI
    155 
    156      do K=2,L_LEVELS
    157      
    158         ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    159        
    160         IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    161           m3as = pqmo(ilay,2) / 2.0
    162           m3af = pqmo(ilay,4) / 2.0
    163          
    164           IF ( ilay .lt. 18 ) THEN
    165             m3as = pqmo(18,2) / 2.0 *dz(k)/dz(18)
    166             m3af = pqmo(18,4) / 2.0 *dz(k)/dz(18)
    167           ENDIF
    168 
    169           dtauaer_s     = m3as*rhoaer_s(nw)
    170           dtauaer_f     = m3af*rhoaer_f(nw)
    171           DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
    172          
    173           IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
    174             SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
    175             ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
    176                             / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
    177           ELSE
    178              DHAZE_T(k,nw) = 0.D0
    179              SSA_T(k,nw)   = 1.0
    180              ASF_T(k,nw)   = 1.0
    181           ENDIF
    182          
    183           IF (callclouds.and.firstcall) &
    184             WRITE(*,*) 'WARNING: In optci, optical properties &
    185                        &calculations are not implemented yet'
    186         ELSE
    187           ! Call fixed vertical haze profile of extinction - same for all columns
    188           call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
    189           if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
    190         ENDIF
    191 
    192         DCONT = 0.0d0 ! continuum absorption
    193 
    194         if(continuum.and.(.not.graybody))then
    195            ! include continua if necessary
    196            wn_cont = dble(wnoi(nw))
    197            T_cont  = dble(TMID(k))
    198            do igas=1,ngasmx
    199 
    200               p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
    201 
    202               dtemp=0.0d0
    203               if(igas.eq.igas_N2)then
    204 
    205                  interm = indi(nw,igas,igas)
    206                  call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    207                  indi(nw,igas,igas) = interm
    208 
    209               elseif(igas.eq.igas_H2)then
    210 
    211                  ! first do self-induced absorption
    212                  interm = indi(nw,igas,igas)
    213                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    214                  indi(nw,igas,igas) = interm
    215 
    216                  ! then cross-interactions with other gases
    217                  do jgas=1,ngasmx
    218                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
    219                     dtempc  = 0.0d0
    220                     if(jgas.eq.igas_N2)then
    221                        interm = indi(nw,igas,jgas)
    222                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    223                        indi(nw,igas,jgas) = interm
    224                     endif
    225                     dtemp = dtemp + dtempc
    226                  enddo
    227 
    228                elseif(igas.eq.igas_CH4)then
    229 
    230                  ! first do self-induced absorption
    231                  interm = indi(nw,igas,igas)
    232                  call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    233                  indi(nw,igas,igas) = interm
    234 
    235                  ! then cross-interactions with other gases
    236                  do jgas=1,ngasmx
    237                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
    238                     dtempc  = 0.0d0
    239                     if(jgas.eq.igas_N2)then
    240                        interm = indi(nw,igas,jgas)
    241                        call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    242                        indi(nw,igas,jgas) = interm
    243                     endif
    244                     dtemp = dtemp + dtempc
    245                  enddo
    246 
    247               endif
    248 
    249               DCONT = DCONT + dtemp
    250 
    251            enddo
    252 
    253            DCONT = DCONT*dz(k)
    254 
    255         endif
    256 
    257         do ng=1,L_NGAUSS-1
    258 
    259            ! Now compute TAUGAS
    260 
    261            ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
    262            ! the execution time of optci/v -> ~ factor 2 on the whole radiative
    263            ! transfer on the tested simulations !
    264 
    265            if (corrk_recombin) then
    266              tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
    267            else
    268              tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
    269            endif
    270 
    271            KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
    272            KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
    273            KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
    274            KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
    275 
    276 
    277            ! Interpolate the gaseous k-coefficients to the requested T,P values
    278 
    279            ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
    280                 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
    281 
    282            TAUGAS  = U(k)*ANS
    283 
    284            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    285            DTAUKI(K,nw,ng) = TAUGAS    &
    286                              + DCONT   & ! For parameterized continuum absorption
    287                              + DHAZE_T(K,NW)  ! For Titan haze
    288 
    289         end do
    290 
    291         ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
    292         ! which holds continuum opacity only
    293 
    294         NG              = L_NGAUSS
    295         DTAUKI(K,nw,ng) = 0.d0      &
    296                           + DCONT   & ! For parameterized continuum absorption
    297                           + DHAZE_T(K,NW)     ! For Titan Haze
    298 
    299      end do
    300   end do
     157         
     158      call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
     159
     160      do LK=1,4
     161         LKCOEF(K,LK) = LCOEF(LK)
     162      end do
     163   end do ! L_LEVELS
     164
     165   do NW=1,L_NSPECTI
     166      ! We ignore K=1...
     167      do K=2,L_LEVELS
     168         ! int. arithmetic => gives the gcm layer index (reversed)
     169         ilay = L_NLAYRAD+1 - k/2
     170
     171         ! Optics coupled with the microphysics :
     172         IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
     173
     174            !==========================================================================================
     175            ! Old optics (must have callclouds = .false.):
     176            !==========================================================================================
     177            IF (.NOT. opt4clouds) THEN
     178               m3as = pqmo(ilay,2) / 2.0
     179               m3af = pqmo(ilay,4) / 2.0
     180               ! Cut-off
     181               IF (ilay .lt. 19) THEN
     182                  m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     183                  m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     184               ENDIF
     185
     186               dtauaer_s = m3as*rhoaer_s(nw)
     187               dtauaer_f = m3af*rhoaer_f(nw)
     188           
     189            !==========================================================================================
     190            ! New optics :
     191            !==========================================================================================
     192            ELSE
     193               iw = (L_NSPECTI + 1) - NW + L_NSPECTV ! Visible first and return
     194               !-----------------------------
     195               ! HAZE (Spherical + Fractal) :
     196               !-----------------------------
     197               FTYPE = 1
     198
     199               ! Spherical aerosols :
     200               !---------------------
     201               CTYPE = 5
     202               m0as  = pqmo(ilay,1) / 2.0
     203               m3as  = pqmo(ilay,2) / 2.0
     204               ! If not callclouds : must have a cut-off
     205               IF (.NOT. callclouds) THEN
     206                  IF (ilay .lt. 19) THEN
     207                     m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     208                     m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     209                  ENDIF
     210               ENDIF
     211               call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
     212               
     213               ! Fractal aerosols :
     214               !-------------------
     215               CTYPE = FTYPE
     216               m0af  = pqmo(ilay,3) / 2.0
     217               m3af  = pqmo(ilay,4) / 2.0
     218               ! If not callclouds : must have a cut-off
     219               IF (.NOT. callclouds) THEN
     220                  IF (ilay .lt. 19) THEN
     221                     m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     222                     m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     223                  ENDIF
     224               ENDIF
     225               call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
     226            ENDIF
     227
     228            ! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
     229            DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
     230            IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
     231               SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
     232               ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
     233                             / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
     234            ELSE
     235               DHAZE_T(k,nw) = 0.D0
     236               SSA_T(k,nw)   = 1.0
     237               ASF_T(k,nw)   = 1.0
     238            ENDIF
     239            ! Diagnostics for the haze :
     240            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
     241            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
     242            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
     243
     244            !---------------------
     245            ! CLOUDS (Spherical) :
     246            !---------------------
     247            IF (callclouds) THEN
     248               CTYPE = 0
     249               m0ccn = pqmo(ilay,5) / 2.0
     250               m3ccn = pqmo(ilay,6) / 2.0
     251               m3ices = 0.0d0
     252               m3cld  = m3ccn
     253               
     254               ! Clear / Dark column method :
     255               !-----------------------------
     256
     257               ! CCN's SSA :
     258               call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
     259               
     260               ! Clear column (CCN, C2H2, C2H6, HCN) :
     261               IF (CDCOLUMN == 0) THEN
     262                  DO iq = 2, nice
     263                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
     264                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
     265                  ENDDO
     266                  call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
     267               
     268               ! Dark column (CCN, CH4, C2H2, C2H6, HCN) :
     269               ELSEIF (CDCOLUMN == 1) THEN
     270                  DO iq = 1, nice
     271                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
     272                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
     273                  ENDDO
     274                  call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
     275               
     276               ELSE
     277                  WRITE(*,*) 'WARNING OPTCI.F90 : CDCOLUMN MUST BE 0 OR 1'
     278                  WRITE(*,*) 'WE USE DARK COLUMN ...'
     279                  DO iq = 1, nice
     280                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
     281                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
     282                  ENDDO
     283                  call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
     284               ENDIF
     285               
     286               ! For small dropplets, opacity of nucleus dominates...
     287               !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices)
     288               ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
     289
     290               ! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
     291               DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
     292               IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
     293                  SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
     294                  ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) + dtau_cld*ssa_cld(nw)*asf_cld(nw) ) &
     295                                / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f + ssa_cld(nw)*dtau_cld )
     296               ELSE
     297                  DHAZE_T(k,nw) = 0.D0
     298                  SSA_T(k,nw)   = 1.0
     299                  ASF_T(k,nw)   = 1.0
     300               ENDIF
     301
     302               ! Tuning of optical properties (now useless...) :
     303               DHAZE_T(k,nw) = DHAZE_T(k,nw) * (FHIR * (1-SSA_T(k,nw)) + SSA_T(k,nw))
     304               SSA_T(k,nw)   = SSA_T(k,nw) / (FHIR * (1-SSA_T(k,nw)) + SSA_T(k,nw))
     305
     306               ! Diagnostics for clouds :
     307               POPT_CLOUDS(k,nw,1) = DHAZE_T(k,nw) ! dtau
     308               POPT_CLOUDS(k,nw,2) = SSA_T(k,nw)   ! wbar
     309               POPT_CLOUDS(k,nw,3) = ASF_T(k,nw)   ! gbar
     310
     311            ELSE
     312               ! Diagnostics for clouds :
     313               POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
     314               POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
     315               POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
     316            ENDIF
     317           
     318         ! Optics and microphysics no coupled :
     319         ELSE
     320            ! Call fixed vertical haze profile of extinction - same for all columns
     321            call disr_haze(dz(k),plev(k),wnoi(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
     322            if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
     323            ! Diagnostics for the haze :
     324            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
     325            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
     326            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
     327            ! Diagnostics for clouds :
     328            POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
     329            POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
     330            POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
     331         ENDIF ! ENDIF callmufi
     332         
     333         DCONT = 0.0d0 ! continuum absorption
     334
     335         if(continuum.and.(.not.graybody))then
     336            ! include continua if necessary
     337            wn_cont = dble(wnoi(nw))
     338            T_cont  = dble(TMID(k))
     339            do igas=1,ngasmx
     340
     341               p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
     342
     343               dtemp=0.0d0
     344               if(igas.eq.igas_N2)then
     345
     346                  interm = indi(nw,igas,igas)
     347                  call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     348                  indi(nw,igas,igas) = interm
     349
     350               elseif(igas.eq.igas_H2)then
     351
     352                  ! first do self-induced absorption
     353                  interm = indi(nw,igas,igas)
     354                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     355                  indi(nw,igas,igas) = interm
     356
     357                  ! then cross-interactions with other gases
     358                  do jgas=1,ngasmx
     359                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
     360                     dtempc  = 0.0d0
     361                     if(jgas.eq.igas_N2)then
     362                        interm = indi(nw,igas,jgas)
     363                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     364                        indi(nw,igas,jgas) = interm
     365                     endif
     366                     dtemp = dtemp + dtempc
     367                  enddo
     368
     369                  elseif(igas.eq.igas_CH4)then
     370
     371                  ! first do self-induced absorption
     372                  interm = indi(nw,igas,igas)
     373                  call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     374                  indi(nw,igas,igas) = interm
     375
     376                  ! then cross-interactions with other gases
     377                  do jgas=1,ngasmx
     378                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
     379                     dtempc  = 0.0d0
     380                     if(jgas.eq.igas_N2)then
     381                        interm = indi(nw,igas,jgas)
     382                        call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     383                        indi(nw,igas,jgas) = interm
     384                     endif
     385                     dtemp = dtemp + dtempc
     386                  enddo
     387
     388               endif
     389
     390               DCONT = DCONT + dtemp
     391
     392            enddo
     393
     394            DCONT = DCONT*dz(k)
     395
     396         endif
     397
     398         do ng=1,L_NGAUSS-1
     399
     400            ! Now compute TAUGAS
     401
     402            ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
     403            ! the execution time of optci/v -> ~ factor 2 on the whole radiative
     404            ! transfer on the tested simulations !
     405
     406            if (corrk_recombin) then
     407               tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
     408            else
     409               tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     410            endif
     411
     412            KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
     413            KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
     414            KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
     415            KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
     416
     417
     418            ! Interpolate the gaseous k-coefficients to the requested T,P values
     419
     420            ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
     421                  LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
     422
     423            TAUGAS  = U(k)*ANS
     424
     425            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
     426            DTAUKI(K,nw,ng) = TAUGAS    &
     427                              + DCONT   & ! For parameterized continuum absorption
     428               + DHAZE_T(K,NW)  ! For Titan haze
     429
     430         end do
     431
     432         ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
     433         ! which holds continuum opacity only
     434
     435         NG              = L_NGAUSS
     436         DTAUKI(K,nw,ng) = 0.d0      &
     437                           + DCONT   & ! For parameterized continuum absorption
     438                        + DHAZE_T(K,NW)     ! For Titan Haze
     439
     440      end do ! k = L_LEVELS
     441   end do ! nw = L_NSPECTI
    301442
    302443  !=======================================================================
    … …  
    418559
    419560  ! Total extinction optical depths
    420 
    421   DO NG=1,L_NGAUSS       ! full gauss loop
    422      DO NW=1,L_NSPECTI       
    423         TAUCUMI(1,NW,NG)=0.0D0
    424         DO K=2,L_LEVELS
    425            TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
    426         END DO
    427      END DO                 ! end full gauss loop
    428   END DO
     561  !DO NG=1,L_NGAUSS       ! full gauss loop
     562  !   DO NW=1,L_NSPECTI       
     563  !      TAUCUMI(1,NW,NG)=0.0D0
     564  !      DO K=2,L_LEVELS
     565  !         TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
     566  !      END DO
     567  !   END DO                 ! end full gauss loop
     568  !END DO
     569 
     570  TAUCUMI(:,:,:) = DTAUKI(:,:,:)
    429571
    430572  ! be aware when comparing with textbook results

  • trunk/LMDZ.TITAN/libf/phytitan/optcv.F90

    r2366 r3083  
    1 SUBROUTINE OPTCV(PQMO,NLAY,PLEV,TMID,PMID,  &
    2      DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT)
     1SUBROUTINE OPTCV(PQMO,NLAY,ZLEV,PLEV,TMID,PMID,  &
     2     DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT,&
     3     POPT_HAZE,POPT_CLOUDS,CDCOLUMN)
    34
    45  use radinc_h
    … …  
    78  use gases_h
    89  use datafile_mod, only: haze_opt_file
    9   use comcstfi_mod, only: r
    10   use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin,     &
    11                           callclouds,callmufi,seashaze,uncoupl_optic_haze
    12   use tracer_h, only: nmicro,nice
     10  use comcstfi_mod, only: pi,r
     11  use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin,   &
     12                          callclouds,callmufi,seashaze,uncoupl_optic_haze,&
     13                          opt4clouds,FHVIS, Fssa
     14  use tracer_h, only: nmicro,nice,ices_indx
    1315
    1416  implicit none
    … …  
    2426  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
    2527  !     Clean and adaptation to Titan by J. Vatant d'Ollone (2016-17)
     28  !     Clean and correction to Titan by B. de Batz de Trenquelléon (2022)
     29  !     New optics added for Titan's clouds by B. de Batz de Trenquelléon (2023)
     30  !     The whole routine needs to be redone...
    2631  !     
    2732  !==================================================================
    … …  
    4752  REAL*8, INTENT(IN)  :: PQMO(nlay,nmicro)  ! Tracers for microphysics optics (X/m2).
    4853  INTEGER, INTENT(IN) :: NLAY               ! Number of pressure layers (for pqmo)
     54  REAL*8, INTENT(IN)  :: ZLEV(NLAY+1)
    4955  REAL*8, INTENT(IN)  :: PLEV(L_LEVELS)
    5056  REAL*8, INTENT(IN)  :: TMID(L_LEVELS), PMID(L_LEVELS)
    5157  REAL*8, INTENT(IN)  :: TAURAY(L_NSPECTV)
    5258  REAL*8, INTENT(IN)  :: SEASHAZEFACT(L_LEVELS)
     59  INTEGER, INTENT(IN) :: CDCOLUMN
    5360 
    5461  REAL*8, INTENT(OUT) :: DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    … …  
    5865  REAL*8, INTENT(OUT) :: WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    5966  REAL*8, INTENT(OUT) :: TAUGSURF(L_NSPECTV,L_NGAUSS-1)
     67  REAL*8, INTENT(OUT) :: POPT_HAZE(L_LEVELS,L_NSPECTI,3)
     68  REAL*8, INTENT(OUT) :: POPT_CLOUDS(L_LEVELS,L_NSPECTI,3)
    6069  ! ==========================================================
    6170 
    … …  
    6473  ! Titan customisation
    6574  ! J. Vatant d'Ollone (2016)
    66   real*8 DHAZE_T(L_LEVELS,L_NSPECTI)
    67   real*8 DHAZES_T(L_LEVELS,L_NSPECTI)
    68   real*8 SSA_T(L_LEVELS,L_NSPECTI)
    69   real*8 ASF_T(L_LEVELS,L_NSPECTI)
     75  real*8 DHAZE_T(L_LEVELS,L_NSPECTV)
     76  real*8 DHAZES_T(L_LEVELS,L_NSPECTV)
     77  real*8 SSA_T(L_LEVELS,L_NSPECTV)
     78  real*8 ASF_T(L_LEVELS,L_NSPECTV)
    7079  ! ==========================
    7180
    … …  
    105114  real*8  dumwvl
    106115
    107   real*8 m3as,m3af
    108   real*8 dtauaer_s,dtauaer_f
     116  ! Variables for new optics
     117  integer iq, iw, FTYPE, CTYPE
     118  real*8 m0as,m0af,m0ccn,m3as,m3af,m3ccn,m3ices,m3cld
     119  real*8 dtauaer_s,dtauaer_f,dtau_ccn,dtau_cld
    109120  real*8,save :: rhoaer_s(L_NSPECTV),ssa_s(L_NSPECTV),asf_s(L_NSPECTV)
    110121  real*8,save :: rhoaer_f(L_NSPECTV),ssa_f(L_NSPECTV),asf_f(L_NSPECTV)
    111 !$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,asf_s,asf_f)
     122  real*8,save :: ssa_ccn(L_NSPECTI),asf_ccn(L_NSPECTI)
     123  real*8,save :: ssa_cld(L_NSPECTV),asf_cld(L_NSPECTV)
     124!$OMP THREADPRIVATE(rhoaer_s,rhoaer_f,ssa_s,ssa_f,ssa_cld,asf_s,asf_f,asf_cld)
    112125 
    113126  logical,save :: firstcall=.true.
    … …  
    141154  ENDIF
    142155
    143   !=======================================================================
    144   !     Determine the total gas opacity throughout the column, for each
    145   !     spectral interval, NW, and each Gauss point, NG.
    146   !     Calculate the continuum opacities, i.e., those that do not depend on
    147   !     NG, the Gauss index.
    148 
    149   taugsurf(:,:) = 0.0
    150   dpr(:)        = 0.0
    151   lkcoef(:,:)   = 0.0
    152 
    153   do K=2,L_LEVELS
    154  
    155      ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    156  
    157      DPR(k) = PLEV(K)-PLEV(K-1)
    158 
    159      ! if we have continuum opacities, we need dz
    160 
     156   !=======================================================================
     157   !     Determine the total gas opacity throughout the column, for each
     158   !     spectral interval, NW, and each Gauss point, NG.
     159   !     Calculate the continuum opacities, i.e., those that do not depend on
     160   !     NG, the Gauss index.
     161
     162   taugsurf(:,:) = 0.0
     163   dpr(:)        = 0.0
     164   lkcoef(:,:)   = 0.0
     165
     166   do K=2,L_LEVELS
     167      ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
     168      DPR(k) = PLEV(K)-PLEV(K-1)
     169
     170      ! if we have continuum opacities, we need dz
    161171      dz(k) = dpr(k)*R*TMID(K)/(gzlat_ig(ilay)*PMID(K))
    162172      U(k)  = Cmk(ilay)*DPR(k)     ! only Cmk line in optcv.F     
    163173
    164      call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
    165 
    166      do LK=1,4
    167         LKCOEF(K,LK) = LCOEF(LK)
    168      end do
    169   end do                    ! levels
    170 
    171   ! Rayleigh scattering
    172   do NW=1,L_NSPECTV
    173      TRAY(1:4,NW)   = 1.d-30
    174      do K=5,L_LEVELS
    175         TRAY(K,NW)   = TAURAY(NW) * DPR(K)
    176      end do                    ! levels
    177   end do
    178  
    179   !     we ignore K=1...
    180   do K=2,L_LEVELS
    181  
    182      ilay = L_NLAYRAD+1 - k/2 ! int. arithmetic => gives the gcm layer index (reversed)
    183 
    184      do NW=1,L_NSPECTV
    185      
    186         IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    187           m3as = pqmo(ilay,2) / 2.0
    188           m3af = pqmo(ilay,4) / 2.0
    189          
    190           IF ( ilay .lt. 18 ) THEN
    191             m3as = pqmo(18,2) / 2.0 * dz(k) / dz(18)
    192             m3af = pqmo(18,4) / 2.0 * dz(k) / dz(18)
    193           ENDIF
    194 
    195           dtauaer_s     = m3as*rhoaer_s(nw)
    196           dtauaer_f     = m3af*rhoaer_f(nw)
    197           DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
    198 
    199           IF ( dtauaer_s + dtauaer_f .GT. 1.D-30 ) THEN
    200             SSA_T(k,nw)   = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
    201             ASF_T(k,nw)   = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
    202                             / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
    203           ELSE
    204              DHAZE_T(k,nw) = 0.D0
    205              SSA_T(k,nw)   = 1.0
    206              ASF_T(k,nw)   = 1.0
    207           ENDIF
    208          
    209           IF (callclouds.and.firstcall) &
    210             WRITE(*,*) 'WARNING: In optcv, optical properties &
    211                        &calculations are not implemented yet'
    212         ELSE
    213           ! Call fixed vertical haze profile of extinction - same for all columns
    214           call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
    215           if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
    216         ENDIF
    217        
    218         !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
    219         !   but visible does not handle very well diffusion in first layer.
    220         !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
    221         !   in the 4 first semilayers in optcv, but not optci.
    222         !   This solves random variations of the sw heating at the model top.
    223         if (k<5)  DHAZE_T(K,:) = 0.0
     174      call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
     175
     176      do LK=1,4
     177         LKCOEF(K,LK) = LCOEF(LK)
     178      end do
     179   end do ! L_LEVELS
     180
     181   ! Rayleigh scattering
     182   do NW=1,L_NSPECTV
     183      TRAY(1:4,NW)   = 1.d-30
     184      do K=5,L_LEVELS
     185         TRAY(K,NW)   = TAURAY(NW) * DPR(K)
     186      end do ! L_LEVELS
     187   end do
     188 
     189   do NW=1,L_NSPECTV
     190      ! We ignore K=1...
     191      do K=2,L_LEVELS
     192         ! int. arithmetic => gives the gcm layer index (reversed)
     193         ilay = L_NLAYRAD+1 - k/2
    224194         
    225         DRAYAER = TRAY(K,NW)
    226         !     DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
    227         DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol
    228 
    229         DCONT = 0.0 ! continuum absorption
    230 
    231         if(continuum.and.(.not.graybody).and.callgasvis)then
    232            ! include continua if necessary
    233            wn_cont = dble(wnov(nw))
    234            T_cont  = dble(TMID(k))
    235            do igas=1,ngasmx
    236 
    237               p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
    238 
    239               dtemp=0.0
    240               if(igas.eq.igas_N2)then
    241 
    242                  interm = indv(nw,igas,igas)
    243 !                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    244                  indv(nw,igas,igas) = interm
    245                  ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
    246 
    247               elseif(igas.eq.igas_H2)then
    248 
    249                  ! first do self-induced absorption
    250                  interm = indv(nw,igas,igas)
    251                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    252                  indv(nw,igas,igas) = interm
    253 
    254                  ! then cross-interactions with other gases
    255                  do jgas=1,ngasmx
    256                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
    257                     dtempc  = 0.0
    258                     if(jgas.eq.igas_N2)then
    259                        interm = indv(nw,igas,jgas)
    260                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    261                        indv(nw,igas,jgas) = interm
    262                        ! should be irrelevant in the visible
    263                     endif
    264                     dtemp = dtemp + dtempc
    265                  enddo
    266 
    267                elseif(igas.eq.igas_CH4)then
    268 
    269                  ! first do self-induced absorption
    270                  interm = indv(nw,igas,igas)
    271                  call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    272                  indv(nw,igas,igas) = interm
    273 
    274                  ! then cross-interactions with other gases
    275                  do jgas=1,ngasmx
    276                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
    277                     dtempc  = 0.0
    278                     if(jgas.eq.igas_N2)then
    279                        interm = indv(nw,igas,jgas)
    280                        call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    281                        indv(nw,igas,jgas) = interm
    282                     endif
    283                     dtemp = dtemp + dtempc
    284                  enddo
    285 
    286               endif
    287 
    288               DCONT = DCONT + dtemp
    289 
    290            enddo
    291 
    292            DCONT = DCONT*dz(k)
    293 
    294         endif
    295 
    296         do ng=1,L_NGAUSS-1
    297 
    298            ! Now compute TAUGAS
    299 
    300            ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
    301            ! the execution time of optci/v -> ~ factor 2 on the whole radiative
    302            ! transfer on the tested simulations !
    303 
    304            if (corrk_recombin) then
    305              tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
    306            else
    307              tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
    308            endif
    309              
    310            KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
    311            KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
    312            KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
    313            KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
    314 
    315            ! Interpolate the gaseous k-coefficients to the requested T,P values
    316 
    317            ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
    318                 LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
    319 
    320 
    321            TAUGAS  = U(k)*ANS
    322 
    323            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    324            DTAUKV(K,nw,ng) = TAUGAS &
    325                              + DRAYAER & ! DRAYAER includes all scattering contributions
    326                              + DCONT ! For parameterized continuum aborption
    327 
    328         end do
    329 
    330         ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
    331         ! which holds continuum opacity only
    332 
    333         NG              = L_NGAUSS
    334         DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze
    335 
    336      end do
    337   end do
    338 
     195         ! Optics coupled with the microphysics :
     196         IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
     197
     198            !==========================================================================================
     199            ! Old optics (must have callclouds = .false.):
     200            !==========================================================================================
     201            IF (.NOT. opt4clouds) THEN
     202               m3as = pqmo(ilay,2) / 2.0
     203               m3af = pqmo(ilay,4) / 2.0
     204               ! Cut-off
     205               IF (ilay .lt. 19) THEN
     206                  m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     207                  m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     208               ENDIF
     209               
     210               dtauaer_s = m3as*rhoaer_s(nw)
     211               dtauaer_f = m3af*rhoaer_f(nw)
     212
     213            !==========================================================================================
     214            ! New optics :
     215            !==========================================================================================
     216            ELSE
     217               iw = (L_NSPECTV + 1) - NW ! Visible first and return
     218               !-----------------------------
     219               ! HAZE (Spherical + Fractal) :
     220               !-----------------------------
     221               FTYPE = 1
     222
     223               ! Spherical aerosols :
     224               !---------------------
     225               CTYPE = 5
     226               m0as  = pqmo(ilay,1) / 2.0
     227               m3as  = pqmo(ilay,2) / 2.0
     228               ! If not callclouds : must have a cut-off
     229               IF (.NOT. callclouds) THEN
     230                  IF (ilay .lt. 19) THEN
     231                     m0as = pqmo(19,1) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     232                     m3as = pqmo(19,2) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     233                  ENDIF
     234               ENDIF
     235               call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0as,m3as,iw,dtauaer_s,ssa_s(nw),asf_s(nw))
     236
     237               ! Fractal aerosols :
     238               !-------------------
     239               CTYPE = FTYPE
     240               m0af  = pqmo(ilay,3) / 2.0
     241               m3af  = pqmo(ilay,4) / 2.0
     242               ! If not callclouds : must have a cut-off
     243               IF (.NOT. callclouds) THEN
     244                  IF (ilay .lt. 19) THEN
     245                     m0af = pqmo(19,3) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     246                     m3af = pqmo(19,4) / 2.0 * (zlev(ilay+1)-zlev(ilay)) / (zlev(20)-zlev(19))
     247                  ENDIF
     248               ENDIF
     249               call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0af,m3af,iw,dtauaer_f,ssa_f(nw),asf_f(nw))
     250            ENDIF
     251
     252            ! Total of Haze opacity (dtau), SSA (w) and ASF (COS) :
     253            DHAZE_T(k,nw) = dtauaer_s + dtauaer_f
     254            IF (dtauaer_s + dtauaer_f .GT. 1.D-30) THEN
     255               SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) ) / ( dtauaer_s+dtauaer_f )
     256               ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) )  &
     257                             / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f )
     258            ELSE
     259               DHAZE_T(k,nw) = 0.D0
     260               SSA_T(k,nw)   = 1.0
     261               ASF_T(k,nw)   = 1.0
     262            ENDIF
     263            ! Diagnostics for the haze :
     264            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
     265            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
     266            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
     267
     268            !---------------------
     269            ! CLOUDS (Spherical) :
     270            !---------------------
     271            IF (callclouds) THEN
     272               CTYPE = 0
     273               m0ccn = pqmo(ilay,5) / 2.0
     274               m3ccn = pqmo(ilay,6) / 2.0
     275               m3ices = 0.0d0
     276               m3cld  = m3ccn
     277               
     278               ! Clear / Dark column method :
     279               !-----------------------------
     280
     281               ! CCN's SSA :
     282               call get_haze_and_cloud_opacity(FTYPE,FTYPE,m0ccn,m3ccn,iw,dtau_ccn,ssa_ccn(nw),asf_ccn(nw))
     283
     284               ! Clear column (CCN, C2H2, C2H6, HCN) :
     285               IF (CDCOLUMN == 0) THEN
     286                  DO iq = 2, nice
     287                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
     288                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
     289                  ENDDO
     290                  call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
     291               
     292               ! Dark column (CCN, CH4, C2H2, C2H6, HCN) :
     293               ELSEIF (CDCOLUMN == 1) THEN
     294                  DO iq = 1, nice
     295                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
     296                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
     297                  ENDDO
     298                  call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
     299               
     300               ELSE
     301                  WRITE(*,*) 'WARNING OPTCV.F90 : CDCOLUMN MUST BE 0 OR 1'
     302                  WRITE(*,*) 'WE USE DARK COLUMN ...'
     303                  DO iq = 1, nice
     304                     m3ices = m3ices + (pqmo(ilay,ices_indx(iq)) / 2.0)
     305                     m3cld  = m3cld + (pqmo(ilay,ices_indx(iq)) / 2.0)
     306                  ENDDO
     307                  call get_haze_and_cloud_opacity(FTYPE,CTYPE,m0ccn,m3cld,iw,dtau_cld,ssa_cld(nw),asf_cld(nw))
     308               ENDIF
     309
     310               ! For small dropplets, opacity of nucleus dominates
     311               !ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3ices) / (m3ccn + m3ices)
     312               ssa_cld(nw) = (ssa_ccn(nw)*m3ccn + ssa_cld(nw)*m3cld) / (m3ccn + m3cld)
     313               ssa_cld(nw) = ssa_cld(nw) * Fssa
     314
     315               ! Total of Haze + Clouds opacity (dtau), SSA (w) and ASF (COS) :
     316               DHAZE_T(k,nw) = dtauaer_s + dtauaer_f + dtau_cld
     317               IF (DHAZE_T(k,nw) .GT. 1.D-30) THEN
     318                  SSA_T(k,nw) = ( dtauaer_s*ssa_s(nw) + dtauaer_f*ssa_f(nw) + dtau_cld*ssa_cld(nw) ) / ( dtauaer_s+dtauaer_f+dtau_cld )
     319                  ASF_T(k,nw) = ( dtauaer_s*ssa_s(nw)*asf_s(nw) + dtauaer_f*ssa_f(nw)*asf_f(nw) + dtau_cld*ssa_cld(nw)*asf_cld(nw) ) &
     320                                / ( ssa_s(nw)*dtauaer_s + ssa_f(nw)*dtauaer_f + ssa_cld(nw)*dtau_cld )
     321               ELSE
     322                  DHAZE_T(k,nw) = 0.D0
     323                  SSA_T(k,nw)   = 1.0
     324                  ASF_T(k,nw)   = 1.0
     325               ENDIF
     326               
     327               ! Tuning of optical properties (now useless...) :
     328               DHAZE_T(k,nw) = DHAZE_T(k,nw) * (FHVIS * (1-SSA_T(k,nw)) + SSA_T(k,nw))
     329               SSA_T(k,nw)   = SSA_T(k,nw) / (FHVIS * (1-SSA_T(k,nw)) + SSA_T(k,nw))
     330
     331               ! Diagnostics for clouds :
     332               POPT_CLOUDS(k,nw,1) = DHAZE_T(k,nw) ! dtau
     333               POPT_CLOUDS(k,nw,2) = SSA_T(k,nw)   ! wbar
     334               POPT_CLOUDS(k,nw,3) = ASF_T(k,nw)   ! gbar
     335           
     336            ELSE
     337               ! Diagnostics for clouds :
     338               POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
     339               POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
     340               POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
     341            ENDIF
     342           
     343         ! Optics and microphysics no coupled :
     344         ELSE
     345            ! Call fixed vertical haze profile of extinction - same for all columns
     346            call disr_haze(dz(k),plev(k),wnov(nw),DHAZE_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
     347            if (seashaze) DHAZE_T(k,nw) = DHAZE_T(k,nw)*seashazefact(k)
     348            ! Diagnostics for the haze :
     349            POPT_HAZE(k,nw,1) = DHAZE_T(k,nw) ! dtau
     350            POPT_HAZE(k,nw,2) = SSA_T(k,nw)   ! wbar
     351            POPT_HAZE(k,nw,3) = ASF_T(k,nw)   ! gbar
     352            ! Diagnostics for clouds :
     353            POPT_CLOUDS(k,nw,1) = 0.D0 ! dtau
     354            POPT_CLOUDS(k,nw,2) = 1.0  ! wbar
     355            POPT_CLOUDS(k,nw,3) = 1.0  ! gbar
     356         ENDIF ! ENDIF callmufi
     357         
     358         !JL18 It seems to be good to have aerosols in the first "radiative layer" of the gcm in the IR
     359         !   but visible does not handle very well diffusion in first layer.
     360         !   The tauaero and tauray are thus set to 0 (a small value for rayleigh because the code crashes otherwise)
     361         !   in the 4 first semilayers in optcv, but not optci.
     362         !   This solves random variations of the sw heating at the model top.
     363         if (k<5)  DHAZE_T(K,:) = 0.0
     364         
     365         DRAYAER = TRAY(K,NW)
     366         ! DRAYAER is Tau RAYleigh scattering, plus AERosol opacity
     367         DRAYAER = DRAYAER + DHAZE_T(K,NW) ! Titan's aerosol
     368
     369         DCONT = 0.0 ! continuum absorption
     370
     371         if(continuum.and.(.not.graybody).and.callgasvis)then
     372            ! include continua if necessary
     373            wn_cont = dble(wnov(nw))
     374            T_cont  = dble(TMID(k))
     375            do igas=1,ngasmx
     376
     377               p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
     378
     379               dtemp=0.0
     380               if(igas.eq.igas_N2)then
     381
     382                  interm = indv(nw,igas,igas)
     383   !                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     384                  indv(nw,igas,igas) = interm
     385                  ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
     386
     387               elseif(igas.eq.igas_H2)then
     388
     389                  ! first do self-induced absorption
     390                  interm = indv(nw,igas,igas)
     391                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     392                  indv(nw,igas,igas) = interm
     393
     394                  ! then cross-interactions with other gases
     395                  do jgas=1,ngasmx
     396                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
     397                     dtempc  = 0.0
     398                     if(jgas.eq.igas_N2)then
     399                        interm = indv(nw,igas,jgas)
     400                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     401                        indv(nw,igas,jgas) = interm
     402                        ! should be irrelevant in the visible
     403                     endif
     404                     dtemp = dtemp + dtempc
     405                  enddo
     406
     407                  elseif(igas.eq.igas_CH4)then
     408
     409                  ! first do self-induced absorption
     410                  interm = indv(nw,igas,igas)
     411                  call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     412                  indv(nw,igas,igas) = interm
     413
     414                  ! then cross-interactions with other gases
     415                  do jgas=1,ngasmx
     416                     p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
     417                     dtempc  = 0.0
     418                     if(jgas.eq.igas_N2)then
     419                        interm = indv(nw,igas,jgas)
     420                        call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     421                        indv(nw,igas,jgas) = interm
     422                     endif
     423                     dtemp = dtemp + dtempc
     424                  enddo
     425
     426               endif
     427
     428               DCONT = DCONT + dtemp
     429
     430            enddo
     431
     432            DCONT = DCONT*dz(k)
     433
     434         endif
     435
     436         do ng=1,L_NGAUSS-1
     437
     438            ! Now compute TAUGAS
     439
     440            ! JVO 2017 : added tmpk because the repeated calls to gasi/v increased dramatically
     441            ! the execution time of optci/v -> ~ factor 2 on the whole radiative
     442            ! transfer on the tested simulations !
     443
     444            if (corrk_recombin) then
     445               tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
     446            else
     447               tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     448            endif
     449               
     450            KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
     451            KCOEF(2) = tmpk(1,2) ! KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
     452            KCOEF(3) = tmpk(2,2) ! KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
     453            KCOEF(4) = tmpk(2,1) ! KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
     454
     455            ! Interpolate the gaseous k-coefficients to the requested T,P values
     456
     457            ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
     458                  LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
     459
     460
     461            TAUGAS  = U(k)*ANS
     462
     463            TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
     464            DTAUKV(K,nw,ng) = TAUGAS &
     465                              + DRAYAER & ! DRAYAER includes all scattering contributions
     466                              + DCONT ! For parameterized continuum aborption
     467         end do
     468
     469         ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
     470         ! which holds continuum opacity only
     471
     472         NG              = L_NGAUSS
     473         DTAUKV(K,nw,ng) = DRAYAER + DCONT ! Scattering + parameterized continuum absorption, including Titan's haze
     474
     475      end do ! k = L_LEVELS
     476   end do ! nw = L_NSPECTI
    339477
    340478  !=======================================================================
    … …  
    355493  ENDDO
    356494
    357 
     495  ! NW spectral loop
    358496  DO NW=1,L_NSPECTV
    359      DO L=1,L_NLAYRAD-1
    360         K              = 2*L+1
    361         atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
    362         btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
    363         ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
    364         btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
    365         COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
    366      END DO ! L vertical loop
     497   ! L vertical loop
     498   DO L=1,L_NLAYRAD-1
     499      K = 2*L+1
     500           atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW) + ASF_T(K+1,NW)*DHAZES_T(K+1,NW)
     501      btemp(L,NW) = DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
     502           ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW)) ! JVO 2017 : does this 0.999 is really meaningful ?
     503           btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
     504      COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
     505   END DO
    367506     
    368      ! Last level
    369      L           = L_NLAYRAD
    370      K           = 2*L+1
    371      atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW)
    372      btemp(L,NW) = DHAZES_T(K,NW)
    373      ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
    374      btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
    375      COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
    376      
    377      
    378   END DO                    ! NW spectral loop
    379 
     507   ! Last level
     508   L = L_NLAYRAD
     509   K = 2*L+1
     510   atemp(L,NW) = ASF_T(K,NW)*DHAZES_T(K,NW)
     511   btemp(L,NW) = DHAZES_T(K,NW)
     512   ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW) ! JVO 2017 : does this 0.999 is really meaningful ?
     513   btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
     514   COSBV(L,NW,1:L_NGAUSS) = atemp(L,NW)/btemp(L,NW)
     515  END DO
     516
     517  ! NG Gauss loop
    380518  DO NG=1,L_NGAUSS
    381     DO NW=1,L_NSPECTV
    382      DO L=1,L_NLAYRAD-1
    383 
    384         K              = 2*L+1
    385         DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
    386         WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
    387 
    388       END DO ! L vertical loop
    389 
    390         ! Last level
    391 
    392         L              = L_NLAYRAD
    393         K              = 2*L+1
    394         DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
    395 
    396         WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
    397 
    398      END DO                 ! NW spectral loop
    399   END DO                    ! NG Gauss loop
     519   ! NW spectral loop
     520   DO NW=1,L_NSPECTV
     521      ! L vertical loop
     522      DO L=1,L_NLAYRAD-1
     523         K = 2*L+1
     524         DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
     525         WBARV(L,nw,ng) = ctemp(L,NW) / DTAUV(L,nw,ng)
     526      END DO
     527
     528      ! Last level
     529      L = L_NLAYRAD
     530      K = 2*L+1
     531           DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
     532      WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
     533   END DO
     534  END DO
    400535
    401536  ! Total extinction optical depths
    402 
    403   DO NG=1,L_NGAUSS       ! full gauss loop
    404      DO NW=1,L_NSPECTV       
    405         TAUCUMV(1,NW,NG)=0.0D0
    406         DO K=2,L_LEVELS
    407            TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
    408         END DO
    409 
    410         DO L=1,L_NLAYRAD
    411            TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
    412         END DO
    413         TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
    414      END DO           
    415   END DO                 ! end full gauss loop
     537  !DO NG=1,L_NGAUSS       ! full gauss loop
     538  !   DO NW=1,L_NSPECTV       
     539  !      TAUCUMV(1,NW,NG)=0.0D0
     540  !      DO K=2,L_LEVELS
     541  !         TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
     542  !      END DO
     543  !      DO L=1,L_NLAYRAD
     544  !         TAUV(L,NW,NG)=TAUCUMV(2*L,NW,NG)
     545  !      END DO
     546  !      TAUV(L,NW,NG)=TAUCUMV(2*L_NLAYRAD+1,NW,NG)
     547  !   END DO           
     548  !END DO                 ! end full gauss loop
     549 
     550  TAUCUMV(:,:,:) = DTAUKV(:,:,:)
     551  DO L=1,L_NLAYRAD
     552   TAUV(L,:,:)=TAUCUMV(2*L,:,:)
     553  END DO
     554  TAUV(L,:,:)=TAUCUMV(2*L_NLAYRAD+1,:,:)
    416555
    417556  if(firstcall) firstcall = .false.

  • trunk/LMDZ.TITAN/libf/phytitan/phys_state_var_mod.F90

    r2793 r3083  
    7171      real,dimension(:,:,:),allocatable,save :: int_dtaui   ! IR optical thickness of layers within narrowbands for diags ().
    7272!$OMP THREADPRIVATE(int_dtaui,int_dtauv)
     73
     74      real,dimension(:,:,:,:),allocatable,save :: zpopthi ! IR optical properties of haze within narrowbands for diags (dtau,tau,k,wbar,gbar).
     75      real,dimension(:,:,:,:),allocatable,save :: zpopthv ! VI optical properties of haze within narrowbands for diags (dtau,tau,k,wbar,gbar).
     76      real,dimension(:,:,:,:),allocatable,save :: zpoptci ! IR optical properties of clouds within narrowbands for diags (dtau,tau,k,wbar,gbar).
     77      real,dimension(:,:,:,:),allocatable,save :: zpoptcv ! VI optical properties of clouds within narrowbands for diags (dtau,tau,k,wbar,gbar).
     78!$OMP THREADPRIVATE(zpopthi,zpopthv,zpoptci,zpoptcv)
     79
     80      real,dimension(:,:),allocatable,save :: u_ref ! Reference profile for the zonal wind nudging (m.s-1).
     81!$OMP THREADPRIVATE(u_ref)
    7382
    7483      real,allocatable,dimension(:,:),save :: qsurf_hist
    … …  
    130139        ALLOCATE(int_dtaui(klon,klev,L_NSPECTI))
    131140        ALLOCATE(int_dtauv(klon,klev,L_NSPECTV))
     141        ALLOCATE(zpopthi(klon,klev,L_NSPECTI,5))
     142        ALLOCATE(zpopthv(klon,klev,L_NSPECTV,5))
     143        ALLOCATE(zpoptci(klon,klev,L_NSPECTI,5))
     144        ALLOCATE(zpoptcv(klon,klev,L_NSPECTV,5))
     145        ALLOCATE(u_ref(49,24))
    132146        allocate(dycchi(klon,klev,nkim))
    133147        ! This is defined in comsaison_h
    … …  
    191205        DEALLOCATE(int_dtaui)
    192206        DEALLOCATE(int_dtauv)
     207        DEALLOCATE(zpopthi)
     208        DEALLOCATE(zpopthv)
     209        DEALLOCATE(zpoptci)
     210        DEALLOCATE(zpoptcv)
     211        DEALLOCATE(u_ref)
    193212        DEALLOCATE(dycchi)
    194213        DEALLOCATE(mu0)

  • trunk/LMDZ.TITAN/libf/phytitan/physiq_mod.F90

    r2742 r3083  
    2020      use comsaison_h, only: mu0, fract, dist_star, declin, right_ascen
    2121      use comsoil_h, only: nsoilmx, layer, mlayer, inertiedat
    22       use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file
    23       use geometry_mod, only: latitude, longitude, cell_area
     22      use datafile_mod, only: datadir, corrkdir, banddir, haze_opt_file, nudging_file
     23      use geometry_mod, only: latitude, latitude_deg, longitude, cell_area
    2424      USE comgeomfi_h, only: totarea, totarea_planet
    2525      USE tracer_h
    … …  
    162162!                + Titan's chemistry
    163163!           Microphysical moment model - J.Burgalat / J.Vatant d'Ollone (2017-2018)
     164!           Improvement of microphysical moment model - B. de Batz de Trenquelléon (2022-2023)
     165!           Optics ofnn Titan's coulds - B. de Batz de Trenquelléon (2022-2023)
    164166!============================================================================================
    165167
    166 
    167 !    0.  Declarations :
    168 !    ------------------
     168! ---------------
     169! 0. DECLARATIONS
     170! ---------------
    169171
    170172include "netcdf.inc"
    … …  
    219221      real pw(ngrid,nlayer)               ! Vertical velocity (m/s). (NOTE : >0 WHEN DOWNWARDS !!)
    220222
    221       integer l,ig,ierr,iq,nw,isoil
    222      
     223      integer l,ig,ierr,iq,nw,isoil,ilat,lat_idx,i,j
     224
    223225      ! FOR DIAGNOSTIC :
    224226
    … …  
    236238      real zdtsurf(ngrid)     ! Cumulated tendencies.
    237239      real zdtsurfmr(ngrid)   ! Mass_redistribution routine.
     240      real zdtsurfevap(ngrid) ! Evaporation.
    238241      real zdtsdif(ngrid)     ! Turbdiff/vdifc routines.
    239242           
    240243      ! For Atmospheric Temperatures : (K/s)   
    241       real zdtdif(ngrid,nlayer)                               ! Turbdiff/vdifc routines.
    242       real zdtmr(ngrid,nlayer)                                ! Mass_redistribution routine.
    243       real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer)         ! Callcorrk routine.
     244      real zdtdif(ngrid,nlayer)                       ! Turbdiff/vdifc routines.
     245      real zdtmr(ngrid,nlayer)                        ! Mass_redistribution routine.
     246      real zdtsw1(ngrid,nlayer), zdtlw1(ngrid,nlayer) ! Callcorrk routine.
     247      real zdtlc(ngrid,nlayer)                        ! Condensation heating rate.
    244248                             
    245249      ! For Surface Tracers : (kg/m2/s)
    … …  
    267271      real zdhdif(ngrid,nlayer)                      ! Turbdiff/vdifc routines.
    268272      real zdhadj(ngrid,nlayer)                      ! Convadj routine.
     273      real zdundg(ngrid,nlayer)                      ! Nudging for zonal wind (m/s/s).
    269274     
    270275      ! For Pressure and Mass :
    … …  
    296301      real zdudyn(ngrid,nlayer)                          ! Dynamical Zonal Wind tendency (m.s-2).
    297302
    298       real zhorizwind(ngrid,nlayer) ! Horizontal Wind ( sqrt(u**+v*v))
     303      real zhorizwind(ngrid,nlayer) ! Horizontal Wind ( sqrt(u*u+v*v))
    299304
    300305      real vmr(ngrid,nlayer)                        ! volume mixing ratio
    … …  
    303308      real ISR,ASR,OLR,GND,DYN,GSR,Ts1,Ts2,Ts3,TsS ! for Diagnostic.
    304309     
    305 !     to test energy conservation (RW+JL)
     310      ! To test energy conservation (RW+JL)
    306311      real mass(ngrid,nlayer),massarea(ngrid,nlayer)
    307312      real dEtot, dEtots, AtmToSurf_TurbFlux
    … …  
    311316      real dEdiffs(ngrid),dEdiff(ngrid)
    312317     
    313 !JL12 conservation test for mean flow kinetic energy has been disabled temporarily
    314 
     318      ! JL12 : conservation test for mean flow kinetic energy has been disabled temporarily
    315319      real dItot, dItot_tmp, dVtot, dVtot_tmp
    316      
    317320      real dWtot, dWtot_tmp, dWtots, dWtots_tmp
    318      
    319321     
    320322      ! For Clear Sky Case.
    … …  
    327329      character(len=10) :: tmp1
    328330      character(len=10) :: tmp2
    329      
    330331      character*2 :: str2
    331332
    … …  
    338339!$OMP THREADPRIVATE(ctimestep)
    339340 
    340       ! Chemical tracers in molar fraction
    341       real, dimension(ngrid,nlayer,nkim)          :: ychim ! (mol/mol)
    342       real, dimension(ngrid,nlayer,nkim)          :: ychimbar ! For 2D chemistry
    343 
    344       ! Molar fraction tendencies ( chemistry, condensation and evaporation ) for tracers (mol/mol/s)
    345       real, dimension(ngrid,nlayer,nq)            :: dyccond        ! Condensation rate. NB : for all tracers, as we want to use indx on it.
    346       real, dimension(ngrid,nlayer,nq)            :: dyccondbar     ! For 2D chemistry
    347       real, dimension(ngrid)                      :: dycevapCH4     ! Surface "pseudo-evaporation" rate (forcing constant surface humidity).
    348 
    349       ! Saturation profiles
    350       real, dimension(ngrid,nlayer,nkim)          :: ysat ! (mol/mol)
    351 
    352 
    353       real :: i2e(ngrid,nlayer)      ! int 2 ext factor ( X.kg-1 -> X.m-3 for diags )
     341      ! Chemical tracers in molar fraction [mol/mol]
     342      real, dimension(ngrid,nlayer,nkim) :: ychim
     343      real, dimension(ngrid,nlayer,nkim) :: ychimbar ! For 2D chemistry
     344
     345      ! Molar fraction tendencies (chemistry, condensation and evaporation) for tracers [mol/mol/s]
     346      real, dimension(ngrid,nlayer,nq)              :: dyccond    ! Condensation rate. NB : for all tracers, as we want to use indx on it.
     347      real, dimension(ngrid,nlayer,size(ices_indx)) :: dmuficond  ! Condensation rate from microphysics [kg/kg/s].
     348      real, dimension(ngrid,nlayer,nq)              :: dyccondbar ! For 2D chemistry
     349      real, dimension(ngrid)                        :: dycevapCH4 ! Surface "pseudo-evaporation" rate for CH4.
     350
     351      ! Saturation profiles [mol/mol]
     352      real, dimension(ngrid,nlayer,nkim) :: ysat
     353     
     354      ! Fraction of CH4 [mol/mol]
     355      real, dimension(ngrid,nlayer) :: tpq_CH4
     356
     357      real :: i2e(ngrid,nlayer) ! int 2 ext factor (X.kg-1 -> X.m-3)
    354358
    355359#ifdef USE_QTEST
    … …  
    411415!        ~~~~~~~~~~~~~~~~~~
    412416         dtrad(:,:) = 0.D0
     417         zdtlc(:,:) = 0.D0
    413418         fluxrad(:) = 0.D0
    414419         zdtsw(:,:) = 0.D0
    415420         zdtlw(:,:) = 0.D0
     421         zpopthi(:,:,:,:) = 0.D0
     422         zpopthv(:,:,:,:) = 0.D0
     423         zpoptci(:,:,:,:) = 0.D0
     424         zpoptcv(:,:,:,:) = 0.D0
    416425
    417426!        Initialize setup for correlated-k radiative transfer
    … …  
    440449#ifndef MESOSCALE
    441450           IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    442              haze_opt_file=trim(datadir)//'/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
     451             haze_opt_file=trim(datadir)//'/optical_tables/HAZE_OPTIC_'//trim(adjustl(tmp1))//'x'//trim(adjustl(tmp2))//'.DAT'
    443452             inquire(file=trim(haze_opt_file),exist=file_ok)
    444453             if(.not.file_ok) then
    … …  
    579588         endif
    580589
     590         ! Read SuperNudging.dat and initialize the nudging for pu
     591         ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
     592         if (nudging_u) then
     593            nudging_file=trim(datadir)//'/SuperNudging.dat'
     594            inquire(file=trim(nudging_file),exist=file_ok)
     595            if(.not.file_ok) then
     596               write(*,*) 'ERROR : The file ',TRIM(nudging_file),' was not found by physiq_mod.F90 ! Check in ', TRIM(datadir)
     597               write(*,*) 'Meanwhile I abort ...'
     598               call abort
     599            endif
     600           
     601            open(88,file=nudging_file,form='formatted')
     602            do ilat = 1, 49
     603               read(88,*) u_ref(ilat,:)
     604            enddo
     605            close(88)
     606         endif
     607
    581608#ifndef MESOSCALE
    582609         if (ngrid.ne.1) then ! Note : no need to create a restart file in 1d.
    … …  
    609636      ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    610637
    611       pdt(:,:)    = 0.D0     
    612       zdtsurf(:)  = 0.D0
    613       pdq(:,:,:)  = 0.D0
    614       dqsurf(:,:) = 0.D0
    615       pdu(:,:)    = 0.D0
    616       pdv(:,:)    = 0.D0
    617       pdpsrf(:)   = 0.D0
    618       zflubid(:)  = 0.D0     
    619       taux(:)     = 0.D0
    620       tauy(:)     = 0.D0
     638      pdt(:,:)       = 0.D0     
     639      zdtsurf(:)     = 0.D0
     640      zdtsurfevap(:) = 0.D0
     641      pdq(:,:,:)     = 0.D0
     642      dqsurf(:,:)    = 0.D0
     643      pdu(:,:)       = 0.D0
     644      pdv(:,:)       = 0.D0
     645      pdpsrf(:)      = 0.D0
     646      zflubid(:)     = 0.D0     
     647      taux(:)        = 0.D0
     648      tauy(:)        = 0.D0
    621649
    622650      zday=pday+ptime ! Compute time, in sols (and fraction thereof).
    … …  
    788816               call stelang(ngrid,sinlon,coslon,sinlat,coslat,    &
    789817                            ztim1,ztim2,ztim3,mu0,fract, flatten)
     818           
    790819            else if(diurnal .eqv. .false.) then
    791820
    … …  
    810839
    811840               ! standard callcorrk
    812                call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,                     &
    813                               albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
    814                               tsurf,fract,dist_star,                              &
     841               call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,                      &
     842                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,zzlev,&
     843                              pt,tsurf,fract,dist_star,                           &
    815844                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
    816845                              fluxsurfabs_sw,fluxtop_lw,                          &
    817846                              fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,                &
    818                               int_dtaui,int_dtauv,lastcall)
     847                              int_dtaui,int_dtauv,zpopthi,zpopthv,zpoptci,zpoptcv,&
     848                              lastcall)
    819849
    820850               ! Radiative flux from the sky absorbed by the surface (W.m-2).
    … …  
    941971           pdt(:,:)=pdt(:,:)+zdtdif(:,:)
    942972         endif
    943            zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
     973         
     974         zdtsurf(:)=zdtsurf(:)+zdtsdif(:)
    944975
    945976         if (tracer) then
    … …  
    10661097            !        - Weird unphysical ratio of m0 and m3, ok for now, but take care of them if
    10671098            !          you want to compute optics from radii.
    1068             WHERE ( (pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0) )
    1069                     pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
    1070                     pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
     1099            WHERE ((pq(:,:,1)+pdq(:,:,1)*ptimestep < 0.D0) .OR. (pq(:,:,2)+pdq(:,:,2)*ptimestep < 0.D0))
     1100               pdq(:,:,1) = (epsilon(1.0)-1.D0)*pq(:,:,1)/ptimestep
     1101               pdq(:,:,2) = (epsilon(1.0)-1.D0)*pq(:,:,2)/ptimestep
    10711102            ENDWHERE
    1072             WHERE ( (pq(:,:,3)+pdq(:,:,3)*ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)*ptimestep < 0.D0) )
    1073                     pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
    1074                     pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
     1103            WHERE ((pq(:,:,3)+pdq(:,:,3)*ptimestep < 0.D0) .OR. (pq(:,:,4)+pdq(:,:,4)*ptimestep < 0.D0))
     1104               pdq(:,:,3) = (epsilon(1.0)-1.D0)*pq(:,:,3)/ptimestep
     1105               pdq(:,:,4) = (epsilon(1.0)-1.D0)*pq(:,:,4)/ptimestep
    10751106            ENDWHERE
     1107            IF (callclouds) THEN
     1108               WHERE ((pq(:,:,5)+pdq(:,:,5)*ptimestep < 0.D0) .OR. (pq(:,:,6)+pdq(:,:,6)*ptimestep < 0.D0))
     1109                  pdq(:,:,5) = (epsilon(1.0)-1.D0)*pq(:,:,5)/ptimestep
     1110                  pdq(:,:,6) = (epsilon(1.0)-1.D0)*pq(:,:,6)/ptimestep
     1111               ENDWHERE
     1112               DO iq = 1, size(ices_indx)
     1113                  ! For ices :
     1114                  WHERE (pq(:,:,ices_indx(iq))+pdq(:,:,ices_indx(iq))*ptimestep < 0.D0)
     1115                     pdq(:,:,ices_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,ices_indx(iq))/ptimestep
     1116                  ENDWHERE
     1117                  ! For gases :
     1118                  WHERE (pq(:,:,gazs_indx(iq))+pdq(:,:,gazs_indx(iq))*ptimestep < 0.D0)
     1119                     pdq(:,:,gazs_indx(iq)) = (epsilon(1.0)-1.D0)*pq(:,:,gazs_indx(iq))/ptimestep
     1120                  ENDWHERE
     1121               ENDDO
     1122            ENDIF
    10761123#endif
    10771124
    … …  
    10831130              ! TODO : Add a sanity check here !
    10841131            endif
    1085 
     1132         
     1133            ! >>> TEMPO : BBT
     1134            ! Default value -> no condensation [kg/kg_air/s] :
     1135            dmuficond(:,:,:) = 0.D0
     1136            do iq = 1, size(ices_indx)
     1137               dmuficond(:,:,iq) = zdqmufi(:,:,gazs_indx(iq))
     1138            enddo
     1139            call cond_muphy(ngrid,nlayer,pt,dmuficond,zdtlc)
     1140            !pdt(:,:) = pdt(:,:) + zdtlc(:,:)
     1141            ! <<< TEMPO : BBT
    10861142         endif
    10871143     
    … …  
    10991155            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    11001156            do iq = 1,nkim
    1101                ychim(:,:,iq) = ( pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro) ) / rat_mmol(iq+nmicro)
     1157               ychim(:,:,iq) = ( pq(:,:,iq+nmicro) + pdq(:,:,iq+nmicro)*ptimestep ) / rat_mmol(iq+nmicro)
    11021158            enddo
    11031159
    … …  
    11101166                ychimbar(:,:,iq) =  zqfibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro)
    11111167                  if ( callclouds ) then
    1112                     ychimbar(:,:,iq) =  ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro) / rat_mmol(iq+nmicro) )
     1168                    ychimbar(:,:,iq) =  ychimbar(:,:,iq) + ( zdqmufibar(:,:,iq+nmicro)*ptimestep / rat_mmol(iq+nmicro) )
    11131169                  endif
    11141170              enddo
    … …  
    11271183            enddo
    11281184
    1129             if ( callclouds ) dyccond(:,:,ices_indx) = 0.D0 ! Condensation have been calculated in the cloud microphysics
     1185            if (callclouds) then
     1186               do iq = 1, size(ices_indx)
     1187                  dyccond(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
     1188               enddo
     1189            endif
    11301190
    11311191            do iq=1,nkim
    1132               ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro) ! update molar ychim for following calchim
     1192              ychim(:,:,iq) = ychim(:,:,iq) + dyccond(:,:,iq+nmicro)*ptimestep ! update molar ychim for following calchim
    11331193
    11341194              pdq(:,:,iq+nmicro) = pdq(:,:,iq+nmicro) + dyccond(:,:,iq+nmicro)*rat_mmol(iq+nmicro) ! convert tendencies to mass mixing ratio
    … …  
    11491209              enddo
    11501210
    1151               if ( callclouds ) dyccondbar(:,:,ices_indx) = 0.D0 ! Condensation have been calculated in the cloud microphysics
     1211               if (callclouds) then
     1212                  do iq = 1, size(ices_indx)
     1213                     dyccondbar(:,:,gazs_indx(iq)) = 0.D0 ! Condensation have been calculated in the cloud microphysics
     1214                  enddo
     1215               endif
    11521216
    11531217              do iq=1,nkim
    1154                 ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)
     1218                ychimbar(:,:,iq) = ychimbar(:,:,iq) + dyccondbar(:,:,iq+nmicro)*ptimestep
    11551219              enddo
    11561220
    11571221              ! Pseudo-evap ( forcing constant surface humidity )
    1158               do ig=1,ngrid
    1159                  if ( ychimbar(ig,1,7) .lt. botCH4 ) ychimbar(ig,1,7) = botCH4
    1160               enddo
     1222              !do ig=1,ngrid
     1223              !   if ( ychimbar(ig,1,7) .lt. botCH4 ) ychimbar(ig,1,7) = botCH4
     1224              !enddo
    11611225
    11621226            endif ! if ( moyzon_ch .and. mod(icount-1,ichim).eq. 0 )
    … …  
    11821246            ! iv. Surface pseudo-evaporation
    11831247            ! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    1184             do ig=1,ngrid
    1185                if ( (ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4 ) then ! +dycchi because ychim not yet updated
    1186                   dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
    1187                else
    1188                   dycevapCH4(ig) = 0.D0
    1189                endif
    1190             enddo
     1248            if (resCH4_inf) then
     1249                do ig=1,ngrid
     1250                        if ( (ychim(ig,1,7)+dycchi(ig,1,7)*ptimestep) .lt. botCH4 ) then ! + dycchi because ychim not yet updated
     1251                                dycevapCH4(ig) = ( -ychim(ig,1,7)+botCH4 ) / ptimestep - dycchi(ig,1,7)
     1252                        else
     1253                                dycevapCH4(ig) = 0.D0
     1254                        endif
     1255                enddo
     1256
     1257            ! >>> New evaporation flux
     1258            else
     1259                tankCH4(:) = 2. ! Pour le moment reservoir infini...
     1260                if (moyzon_ch) then
     1261                    tpq_CH4(:,:) = ychimbar(:,:,7) + dycchi(:,:,7)*ptimestep ! + dycchi because ychim not yet updated [mol/mol]
     1262                else
     1263                    tpq_CH4(:,:) = ychim(:,:,7) + dycchi(:,:,7)*ptimestep    ! + dycchi because ychim not yet updated [mol/mol]
     1264                endif
     1265                call evapCH4(ngrid,nlayer,ptimestep,pplev,zzlay,zzlev,&
     1266                             pu,pv,tsurf,tpq_CH4,tankCH4,dycevapCH4,zdtsurfevap)
     1267                zdtsurf(:) = zdtsurf(:) + zdtsurfevap(:)
     1268            endif
     1269            ! <<< New evaporation flux
    11911270
    11921271            pdq(:,1,7+nmicro) = pdq(:,1,7+nmicro) + dycevapCH4(:)*rat_mmol(7+nmicro)
    … …  
    12041283
    12051284            pdq(:,:,:) = pdq(:,:,:) + zdqchi(:,:,:)
    1206            
     1285
    12071286         endif ! end of 'callchim'
    12081287
    … …  
    12371316            pdv(:,:)    = pdv(:,:)    + zdvmr(:,:)
    12381317            pdpsrf(:)   = pdpsrf(:)   + zdpsrfmr(:)
    1239             zdtsurf(:)  = zdtsurf(:)  + zdtsurfmr(:)
    1240            
    12411318         endif
    12421319
    … …  
    12801357!   VII. Perform diagnostics and write output files
    12811358!---------------------------------------------------
    1282 
     1359     
     1360      ! Nudging of zonal wind !
     1361      ! ~~~~~~~~~~~~~~~~~~~~~~~
     1362      if (nudging_u) then
     1363         zdundg(:,:) = 0.D0
     1364         
     1365         ! boucle sur les points de grille :
     1366         do i = 1, ngrid
     1367            ! boucle sur les latitudes :
     1368            do j = 1, 49
     1369
     1370               ! Nudging of the first 23 layers only !!
     1371               if (ABS(REAL(u_ref(j,1)) - REAL(latitude_deg(i))) .lt. 1e-2) then
     1372                  zdundg(i,1:23) = (u_ref(j,2:24) - (pu(i,1:23)+pdu(i,1:23)*ptimestep)) / nudging_dt
     1373               endif
     1374           
     1375            enddo
     1376         enddo
     1377
     1378         pdu(:,:) = pdu(:,:) + zdundg(:,:)
     1379      endif
     1380       
    12831381      ! Note : For output only: the actual model integration is performed in the dynamics.
    12841382 
    … …  
    13091407
    13101408      ! Tracers.
     1409      ! >>> TEMPO : BBT
     1410      pdq(:,1:13,gazs_indx(2)) = (2e-8 - pq(:,1:13,gazs_indx(2))) / (ptimestep * 100)  ! C2H6
     1411      pdq(:,1:12,gazs_indx(3)) = (1e-9 - pq(:,1:12,gazs_indx(3))) / (ptimestep * 100)  ! C2H2
     1412      pdq(:,1:11,gazs_indx(4)) = (1e-20 - pq(:,1:11,gazs_indx(4))) / (ptimestep * 100) ! HCN
     1413      ! <<< TEMPO : BBT
    13111414      zq(:,:,:) = pq(:,:,:) + pdq(:,:,:)*ptimestep
    13121415
    13131416      ! Surface pressure.
    13141417      ps(:) = pplev(:,1) + pdpsrf(:)*ptimestep
    1315 
    13161418
    13171419
    … …  
    14541556      call writediagfi(ngrid,"p","Pressure","Pa",3,pplay)
    14551557
    1456 !     Subsurface temperatures
    1457 !        call writediagsoil(ngrid,"tsurf","Surface temperature","K",2,tsurf)
    1458 !        call writediagsoil(ngrid,"temp","temperature","K",3,tsoil)
    1459 
    1460      
    1461 
    14621558      ! Total energy balance diagnostics
    14631559      if(callrad.and.(.not.newtonian))then
    1464      
    1465          !call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
     1560         call writediagfi(ngrid,"ALB","Surface albedo"," ",2,albedo_equivalent)
    14661561         call writediagfi(ngrid,"ISR","incoming stellar rad.","W m-2",2,fluxtop_dn)
    14671562         call writediagfi(ngrid,"ASR","absorbed stellar rad.","W m-2",2,fluxabs_sw)
    14681563         call writediagfi(ngrid,"OLR","outgoing longwave rad.","W m-2",2,fluxtop_lw)
    14691564         
    1470 !           call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
    1471 !           call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
    1472 !           call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
    1473 !           call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
    1474 !           call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
    1475 !           call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
    1476 
    1477 
    1478          call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)
    1479          
     1565         call writediagfi(ngrid,"ASRcs","absorbed stellar rad (cs).","W m-2",2,fluxabs_sw1)
     1566         call writediagfi(ngrid,"OLRcs","outgoing longwave rad (cs).","W m-2",2,fluxtop_lw1)
     1567         call writediagfi(ngrid,"fluxsurfsw","sw surface flux.","W m-2",2,fluxsurf_sw)
     1568         call writediagfi(ngrid,"fluxsurflw","lw back radiation.","W m-2",2,fluxsurf_lw)
     1569         call writediagfi(ngrid,"fluxsurfswcs","sw surface flux (cs).","W m-2",2,fluxsurf_sw1)
     1570         call writediagfi(ngrid,"fluxsurflwcs","lw back radiation (cs).","W m-2",2,fluxsurf_lw1)
     1571
     1572         call writediagfi(ngrid,"GND","heat flux from ground","W m-2",2,fluxgrd)         
    14801573         call writediagfi(ngrid,"DYN","dynamical heat input","W m-2",2,fluxdyn)
    1481          
    14821574      endif ! end of 'callrad'
    14831575       
    14841576      if(enertest) then
    1485      
    14861577         if (calldifv) then
    1487          
    14881578            call writediagfi(ngrid,"q2","turbulent kinetic energy","J.kg^-1",3,q2)
    14891579            call writediagfi(ngrid,"sensibFlux","sensible heat flux","w.m^-2",2,sensibFlux)
    1490            
    1491 !             call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
    1492 !             call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
    1493 !             call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
    1494              
     1580            call writediagfi(ngrid,"dEzdiff","turbulent diffusion heating (-sensible flux)","w.m^-2",3,dEzdiff)
     1581            call writediagfi(ngrid,"dEdiff","integrated turbulent diffusion heating (-sensible flux)","w.m^-2",2,dEdiff)
     1582            call writediagfi(ngrid,"dEdiffs","In TurbDiff (correc rad+latent heat) surf nrj change","w.m^-2",2,dEdiffs)
    14951583         endif
    14961584         
    … …  
    14991587            call writediagfi(ngrid,"dEzradlw","radiative heating","w.m^-2",3,dEzradlw)
    15001588         endif
    1501          
    15021589      endif ! end of 'enertest'
    15031590
    1504         ! Diagnostics of optical thickness
    1505         ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
    1506         if (diagdtau) then               
    1507           do nw=1,L_NSPECTV
     1591      ! Diagnostics of optical thickness
     1592      ! Warning this is exp(-tau), I let you postproc with -log to have tau itself - JVO 19
     1593      if (diagdtau) then               
     1594         do nw=1,L_NSPECTV
    15081595            write(str2,'(i2.2)') nw
    15091596            call writediagfi(ngrid,'dtauv'//str2,'Layer optical thickness attenuation in VI band '//str2,'',1,int_dtauv(:,nlayer:1:-1,nw))
    1510           enddo
    1511           do nw=1,L_NSPECTI
     1597         enddo
     1598         do nw=1,L_NSPECTI
    15121599            write(str2,'(i2.2)') nw
    15131600            call writediagfi(ngrid,'dtaui'//str2,'Layer optical thickness attenuation in IR band '//str2,'',1,int_dtaui(:,nlayer:1:-1,nw))
    1514           enddo
    1515         endif
    1516 
    1517         ! Temporary inclusions for winds diagnostics.
    1518         call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
    1519         call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
    1520 
    1521         ! Temporary inclusions for heating diagnostics.
    1522         call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
    1523         call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
    1524         call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
    1525         call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
     1601         enddo
     1602      endif
     1603
     1604      ! Temporary inclusions for winds diagnostics.
     1605      call writediagfi(ngrid,"zdudif","Turbdiff tend. zon. wind","m s-2",3,zdudif)
     1606      call writediagfi(ngrid,"zdudyn","Dyn. tend. zon. wind","m s-2",3,zdudyn)
     1607
     1608      ! Temporary inclusions for heating diagnostics.
     1609      call writediagfi(ngrid,"zdtsw","SW heating","T s-1",3,zdtsw)
     1610      call writediagfi(ngrid,"zdtlw","LW heating","T s-1",3,zdtlw)
     1611      call writediagfi(ngrid,"dtrad","radiative heating","K s-1",3,dtrad)
     1612      call writediagfi(ngrid,"zdtdyn","Dyn. heating","T s-1",3,zdtdyn)
    15261613       
    1527         ! For Debugging.
    1528         call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
    1529        
     1614      ! For Debugging.
     1615      call writediagfi(ngrid,'pphi','Geopotential',' ',3,pphi)
    15301616
    15311617      ! Output tracers.
    … …  
    15611647            call writediagfi(ngrid,"mmd_rc_sph","Spherical mode caracteristic radius",'m',3,mmd_rc_sph)
    15621648            call writediagfi(ngrid,"mmd_rc_fra","Fractal mode caracteristic radius",'m',3,mmd_rc_fra)
    1563 
    15641649         endif ! end of 'callmufi'
    15651650
    … …  
    15681653           do iq=1,nkim
    15691654             call writediagfi(ngrid,cnames(iq),cnames(iq),'mol/mol',3,zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro))
     1655             call writediagfi(ngrid,'dqcond_'//cnames(iq),'dqcond_'//cnames(iq),'mol/mol/s',3,dyccond(:,:,iq+nmicro))
    15701656           enddo
    15711657           call writediagfi(ngrid,"evapCH4","Surface CH4 pseudo-evaporation rate",'mol/mol/s',2,dycevapCH4)
    … …  
    15741660       endif ! end of 'tracer'
    15751661
     1662#ifdef CPP_XIOS
     1663!-----------------------------------------------------------------------------------------------------
    15761664! XIOS outputs
    1577 #ifdef CPP_XIOS     
     1665!-----------------------------------------------------------------------------------------------------
    15781666      ! Send fields to XIOS: (NB these fields must also be defined as
    15791667      ! <field id="..." /> in context_lmdz_physics.xml to be correctly used)
     1668     
     1669      !--------------------------------------------------------
     1670      ! General diagnostics :
     1671      !--------------------------------------------------------
    15801672      CALL send_xios_field("ls",zls*180./pi)
    15811673      CALL send_xios_field("lss",zlss*180./pi)
    … …  
    15831675      CALL send_xios_field("Declin",declin*180./pi)
    15841676     
    1585       ! Total energy balance diagnostics
    1586       if (callrad.and.(.not.newtonian)) then
    1587          CALL send_xios_field("ISR_TOA",fluxtop_dn)
    1588          CALL send_xios_field("OLR_TOA",fluxtop_lw)
    1589       endif
    1590      
    1591       CALL send_xios_field("area",cell_area)
    1592       CALL send_xios_field("pphi",pphi)
    1593       CALL send_xios_field("pphis",phisfi)
    1594      
    1595       CALL send_xios_field("ps",ps)
    1596       CALL send_xios_field("tsurf",tsurf)
    1597 
    1598       if(enertest) then
    1599          if (calldifv) then
    1600             CALL send_xios_field("sensibFlux",sensibFlux)
    1601          endif
    1602       endif
    1603 
     1677      ! Atmosphere (3D) :
    16041678      CALL send_xios_field("temp",zt)
    16051679      CALL send_xios_field("teta",zh)
     1680      CALL send_xios_field("p",pplay)
    16061681      CALL send_xios_field("u",zu)
    16071682      CALL send_xios_field("v",zv)
    16081683      CALL send_xios_field("w",pw)
    1609       CALL send_xios_field("p",pplay)
    1610      
    1611       ! Winds diagnostics.
     1684
     1685      CALL send_xios_field("area",cell_area)
     1686      CALL send_xios_field("pphi",pphi)
     1687
     1688      ! Surface (2D) :
     1689      CALL send_xios_field("ps",ps)
     1690      CALL send_xios_field("tsurf",tsurf)
     1691      CALL send_xios_field("pphis",phisfi)
     1692
     1693      ! Total energy balance diagnostics (2D) :
     1694      IF (callrad.and.(.not.newtonian)) THEN
     1695         CALL send_xios_field("ISR_TOA",fluxtop_dn)
     1696         CALL send_xios_field("OLR_TOA",fluxtop_lw)
     1697      ENDIF
     1698
     1699
     1700      !--------------------------------------------------------
     1701      ! Winds trends :
     1702      !--------------------------------------------------------
     1703      ! Atmosphere (3D) :
     1704      ! du_tot = zdudyn + pdu
     1705      CALL send_xios_field("dudyn",zdudyn)
     1706      CALL send_xios_field("pdu",pdu)
     1707     
     1708      ! pdu = zdudif + zduadj + zdundg
    16121709      CALL send_xios_field("dudif",zdudif)
    1613       CALL send_xios_field("dudyn",zdudyn)
    1614 
     1710      CALL send_xios_field("duadj",zduadj)
     1711      CALL send_xios_field("dundg",zdundg)
     1712     
     1713      ! zhorizwind = sqrt(u*u + v*v)
    16151714      CALL send_xios_field("horizwind",zhorizwind)
    16161715
    1617       ! Heating diagnostics.
     1716
     1717      !--------------------------------------------------------
     1718      ! Heating trends :
     1719      !--------------------------------------------------------
     1720      ! Atmosphere (3D) :
     1721      ! dt_tot = dtdyn + pdt
     1722      CALL send_xios_field("dtdyn",zdtdyn)
     1723      CALL send_xios_field("pdt",pdt)
     1724
     1725      ! pdt = dtrad + zdtdif + dtadj + zdtlc
     1726      CALL send_xios_field("dtrad",dtrad)
     1727      CALL send_xios_field("dtdif",zdtdif)
     1728      CALL send_xios_field("dtadj",zdtadj(:,:))
     1729      CALL send_xios_field("dtlc",zdtlc)
     1730
     1731      ! dtrad = zdtsw + zdtlw
    16181732      CALL send_xios_field("dtsw",zdtsw)
    16191733      CALL send_xios_field("dtlw",zdtlw)
    1620       CALL send_xios_field("dtrad",dtrad)
    1621       CALL send_xios_field("dtdyn",zdtdyn)
    1622       CALL send_xios_field("dtdif",zdtdif)
    1623 
    1624       ! Chemical tracers
    1625       if (callchim) then
    1626      
    1627         ! Advected fields
    1628         do iq=1,nkim
    1629           CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
    1630         enddo
    1631        
    1632         ! Upper chemistry fields
    1633         do iq=1,nkim
    1634           CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
    1635         enddo
    1636        
    1637         ! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
    1638         do iq=1,nkim
    1639          
    1640           ! GCM levels
    1641           do l=1,nlayer
    1642             ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
    1643           enddo
    1644           ! Upper levels
    1645           do l=1,nlaykim_up
    1646             ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
    1647           enddo
    1648          
    1649           CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
    1650          
    1651         enddo
    1652        
    1653         ! Condensation tendencies ( mol/mol/s )
    1654         CALL send_xios_field("dqcond_CH4",dyccond(:,:,7+nmicro))
    1655         CALL send_xios_field("dqcond_C2H2",dyccond(:,:,10+nmicro))
    1656         CALL send_xios_field("dqcond_C2H4",dyccond(:,:,12+nmicro))
    1657         CALL send_xios_field("dqcond_C2H6",dyccond(:,:,14+nmicro))
    1658         CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
    1659         CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
    1660         CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
    1661         CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
    1662         CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
    1663         CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
    1664         CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
    1665         CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
    1666         CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
    1667         CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
    1668         CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
    1669         CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
    1670         CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
    1671 
    1672         ! Pseudo-evaporation flux (mol/mol/s)
    1673         CALL send_xios_field("evapCH4",dycevapCH4(:))
    1674 
    1675       endif ! of 'if callchim'
    1676 
    1677       ! Microphysical tracers
    1678       if (callmufi) then
    1679         CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
    1680         CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
    1681         CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
    1682         CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
    1683       endif       
     1734
     1735      ! Surface (2D) :
     1736      IF(enertest) THEN
     1737         IF (calldifv) THEN
     1738            CALL send_xios_field("sensibFlux",sensibFlux)
     1739         ENDIF
     1740      ENDIF
     1741      CALL send_xios_field("fluxsurf_lw",fluxsurf_lw(:))
     1742      CALL send_xios_field("fluxsurfabs_sw",fluxsurfabs_sw(:))
     1743      CALL send_xios_field("emis",emis(:))
     1744     
     1745      ! dtsurf = dtsdif + dtsurfevap
     1746      CALL send_xios_field("dtsurf",zdtsurf(:))
     1747      CALL send_xios_field("dtsdif",zdtsdif(:))
     1748      CALL send_xios_field("dtsurfevap",zdtsurfevap(:))
     1749     
     1750
     1751      !--------------------------------------------------------
     1752      ! Optical diagnostics :
     1753      !--------------------------------------------------------
     1754      ! Optics diagnostics for haze (Channel VI04 : 2.80um).
     1755      CALL send_xios_field('tauhv_04',zpopthv(:,:,4,2))
     1756      CALL send_xios_field('khv_04',zpopthv(:,:,4,3))
     1757      CALL send_xios_field('wbarhv_04',zpopthv(:,:,4,4))
     1758      CALL send_xios_field('gbarhv_04',zpopthv(:,:,4,5))
     1759      ! Optics diagnostics for haze (Channel VI18 : 0.85um).
     1760      CALL send_xios_field('tauhv_18',zpopthv(:,:,18,2))
     1761      CALL send_xios_field('khv_18',zpopthv(:,:,18,3))
     1762      CALL send_xios_field('wbarhv_18',zpopthv(:,:,18,4))
     1763      CALL send_xios_field('gbarhv_18',zpopthv(:,:,18,5))
     1764      ! Optics diagnostics for haze (Channel IR07 : 29.5um).
     1765      CALL send_xios_field('tauhi_07',zpopthi(:,:,7,2))
     1766      CALL send_xios_field('khi_07',zpopthi(:,:,7,3))
     1767      CALL send_xios_field('wbarhi_07',zpopthi(:,:,7,4))
     1768      CALL send_xios_field('gbarhi_07',zpopthi(:,:,7,5))
     1769      ! Optics diagnostics for haze (Channel IR22 : 5.66um).
     1770      CALL send_xios_field('tauhi_22',zpopthi(:,:,22,2))
     1771      CALL send_xios_field('khi_22',zpopthi(:,:,22,3))
     1772      CALL send_xios_field('wbarhi_22',zpopthi(:,:,22,4))
     1773      CALL send_xios_field('gbarhi_22',zpopthi(:,:,22,5))
     1774
     1775      ! Optics diagnostics for haze + clouds (Channel VI04 : 2.80um).
     1776      CALL send_xios_field('taucv_04',zpoptcv(:,:,4,2))
     1777      CALL send_xios_field('kcv_04',zpoptcv(:,:,4,3))
     1778      CALL send_xios_field('wbarcv_04',zpoptcv(:,:,4,4))
     1779      CALL send_xios_field('gbarcv_04',zpoptcv(:,:,4,5))
     1780      ! Optics diagnostics for haze + clouds (Channel VI18 : 0.85um).
     1781      CALL send_xios_field('taucv_18',zpoptcv(:,:,18,2))
     1782      CALL send_xios_field('kcv_18',zpoptcv(:,:,18,3))
     1783      CALL send_xios_field('wbarcv_18',zpoptcv(:,:,18,4))
     1784      CALL send_xios_field('gbarcv_18',zpoptcv(:,:,18,5))
     1785      ! Optics diagnostics for haze + clouds (Channel IR07 : 29.5um).
     1786      CALL send_xios_field('tauci_07',zpoptci(:,:,7,2))
     1787      CALL send_xios_field('kci_07',zpoptci(:,:,7,3))
     1788      CALL send_xios_field('wbarci_07',zpoptci(:,:,7,4))
     1789      CALL send_xios_field('gbarci_07',zpoptci(:,:,7,5))
     1790      ! Optics diagnostics for haze + clouds (Channel IR22 : 5.66um).
     1791      CALL send_xios_field('tauci_22',zpoptci(:,:,22,2))
     1792      CALL send_xios_field('kci_22',zpoptci(:,:,22,3))
     1793      CALL send_xios_field('wbarci_22',zpoptci(:,:,22,4))
     1794      CALL send_xios_field('gbarci_22',zpoptci(:,:,22,5))
     1795
     1796      ! >>> [TEMPO : BBT]
     1797      !DO nw = 1, L_NSPECTV
     1798      !   WRITE(str2,'(i2.2)') nw
     1799      !   CALL send_xios_field('khv'//TRIM(nw),zpopthv(:,:,nw,3))
     1800      !ENDDO
     1801      ! <<< [TEMPO : BBT]
     1802     
     1803
     1804      !--------------------------------------------------------
     1805      ! Microphysical tracers :
     1806      !--------------------------------------------------------
     1807      IF (callmufi) THEN
     1808         ! Atmosphere (3D) :
     1809         ! Moments M0 and M3 :
     1810         CALL send_xios_field("mu_m0as",zq(:,:,micro_indx(1))*i2e)
     1811         CALL send_xios_field("mu_m3as",zq(:,:,micro_indx(2))*i2e)
     1812         CALL send_xios_field("mu_m0af",zq(:,:,micro_indx(3))*i2e)
     1813         CALL send_xios_field("mu_m3af",zq(:,:,micro_indx(4))*i2e)
     1814         IF (callclouds) THEN
     1815            CALL send_xios_field("mu_m0n",zq(:,:,micro_indx(5))*i2e)
     1816            CALL send_xios_field("mu_m3n",zq(:,:,micro_indx(6))*i2e)
     1817            DO iq = 1, size(ices_indx)
     1818               CALL send_xios_field(TRIM(nameOfTracer(ices_indx(iq))),zq(:,:,ices_indx(iq))*i2e)
     1819            ENDDO
     1820         ENDIF
     1821
     1822         ! Microphysical diagnostics :
     1823         CALL send_xios_field("rc_sph",mmd_rc_sph(:,:))
     1824         CALL send_xios_field("rc_fra",mmd_rc_fra(:,:))
     1825         CALL send_xios_field("vsed_aers",mmd_aer_s_w(:,:))
     1826         CALL send_xios_field("vsed_aerf",mmd_aer_f_w(:,:))
     1827         CALL send_xios_field("flux_aers",mmd_aer_s_flux(:,:))
     1828         CALL send_xios_field("flux_aerf",mmd_aer_f_flux(:,:))
     1829         IF (callclouds) THEN
     1830            CALL send_xios_field("rc_cld",mmd_rc_cld(:,:))
     1831            CALL send_xios_field("vsed_ccn",mmd_ccn_w(:,:))
     1832            CALL send_xios_field("flux_ccn",mmd_ccn_flux(:,:))
     1833            CALL send_xios_field("flux_iCH4",mmd_ice_fluxes(:,:,1))
     1834            CALL send_xios_field("flux_iC2H2",mmd_ice_fluxes(:,:,2))
     1835            CALL send_xios_field("flux_iC2H6",mmd_ice_fluxes(:,:,3))
     1836            CALL send_xios_field("flux_iHCN",mmd_ice_fluxes(:,:,3))
     1837            DO iq = 1, size(ices_indx)
     1838               CALL send_xios_field(TRIM(nameOfTracer(gazs_indx(iq)))//'_sat',mmd_gazs_sat(:,:,iq))
     1839            ENDDO
     1840         ENDIF
     1841         
     1842         ! Surface (2D) :
     1843         CALL send_xios_field("aer_prec",mmd_aer_prec(:))
     1844         CALL send_xios_field("ccn_prec",mmd_ccn_prec(:))
     1845         CALL send_xios_field("iCH4_prec",mmd_ice_prec(:,1))
     1846         CALL send_xios_field("iC2H2_prec",mmd_ice_prec(:,2))
     1847         CALL send_xios_field("iC2H6_prec",mmd_ice_prec(:,3))
     1848         CALL send_xios_field("iHCN_prec",mmd_ice_prec(:,3))
     1849      ENDIF ! of 'if callmufi'
     1850
     1851      !--------------------------------------------------------
     1852      ! Chemical tracers :
     1853      !--------------------------------------------------------
     1854      IF (callchim) THEN
     1855         ! Atmosphere (3D) :
     1856         ! Chemical species :
     1857         DO iq = 1, nkim
     1858            CALL send_xios_field(trim(cnames(iq)),zq(:,:,iq+nmicro)/rat_mmol(iq+nmicro)) ! kg/kg -> mol/mol
     1859         ENDDO
     1860
     1861         ! Condensation tendencies from microphysics (mol/mol/s) :
     1862         DO iq = 1, size(ices_indx)
     1863            CALL send_xios_field('dmuficond_'//trim(nameOfTracer(gazs_indx(iq))),dmuficond(:,:,iq)/rat_mmol(gazs_indx(iq))) ! kg/kg/s -> mol/mol/s
     1864         ENDDO
     1865         
     1866         ! Condensation tendencies (mol/mol/s) :
     1867         CALL send_xios_field("dqcond_CH4",dyccond(:,:,7+nmicro))
     1868         CALL send_xios_field("dqcond_C2H2",dyccond(:,:,10+nmicro))
     1869         CALL send_xios_field("dqcond_C2H4",dyccond(:,:,12+nmicro))
     1870         CALL send_xios_field("dqcond_C2H6",dyccond(:,:,14+nmicro))
     1871         CALL send_xios_field("dqcond_C3H6",dyccond(:,:,17+nmicro))
     1872         CALL send_xios_field("dqcond_C4H4",dyccond(:,:,21+nmicro))
     1873         CALL send_xios_field("dqcond_CH3CCH",dyccond(:,:,23+nmicro))
     1874         CALL send_xios_field("dqcond_C3H8",dyccond(:,:,24+nmicro))
     1875         CALL send_xios_field("dqcond_C4H2",dyccond(:,:,25+nmicro))
     1876         CALL send_xios_field("dqcond_C4H6",dyccond(:,:,26+nmicro))
     1877         CALL send_xios_field("dqcond_C4H10",dyccond(:,:,27+nmicro))
     1878         CALL send_xios_field("dqcond_AC6H6",dyccond(:,:,28+nmicro))
     1879         CALL send_xios_field("dqcond_HCN",dyccond(:,:,35+nmicro))
     1880         CALL send_xios_field("dqcond_CH3CN",dyccond(:,:,39+nmicro))
     1881         CALL send_xios_field("dqcond_HC3N",dyccond(:,:,41+nmicro))
     1882         CALL send_xios_field("dqcond_NCCN",dyccond(:,:,42+nmicro))
     1883         CALL send_xios_field("dqcond_C4N2",dyccond(:,:,43+nmicro))
     1884     
     1885         ! Upper atmosphere chemistry variables (3D) :
     1886         DO iq = 1, nkim
     1887            CALL send_xios_field(trim(cnames(iq))//"_up",ykim_up(iq,:,:)) ! mol/mol
     1888         ENDDO
     1889         
     1890         ! Total atmosphere chemistry variables (3D) :
     1891         ! Append fields in ykim_tot for output on the total vertical grid (0->1300km)
     1892         DO iq = 1, nkim   
     1893            ! GCM levels
     1894            DO l = 1, nlayer
     1895               ykim_tot(iq,:,l) = zq(:,l,iq+nmicro)/rat_mmol(iq+nmicro)
     1896            ENDDO
     1897            ! Upper levels
     1898            DO l = 1, nlaykim_up
     1899               ykim_tot(iq,:,nlayer+l) = ykim_up(iq,:,l)
     1900            ENDDO
     1901            CALL send_xios_field(trim(cnames(iq))//"_tot",ykim_tot(iq,:,:)) ! mol/mol
     1902         ENDDO
     1903
     1904         ! Surface (2D) :
     1905         CALL send_xios_field("evapCH4",dycevapCH4(:)) ! Pseudo-evaporation flux (mol/mol/s)
     1906         CALL send_xios_field("tankCH4",tankCH4(:))    ! CH4 tank at the surface (m)
     1907      ENDIF ! of 'if callchim'
     1908
    16841909
    16851910      if (lastcall.and.is_omp_master) then

  • trunk/LMDZ.TITAN/libf/phytitan/tpindex.F

    r1648 r3083  
    4848!     -------
    4949!     Adapted from the NASA Ames code by R. Wordsworth (2009)
     50!     Last update : B. de Batz de Trenquelleon (2023)
    5051!     
    5152!==================================================================
    … …  
    5556      implicit none
    5657
    57       real*8 Tref(L_NTREF)
    58       real*8 pref(L_PINT)
     58      ! INPUTS :
     59      !=========
     60      REAL,INTENT(IN) :: Tref(L_NTREF)
     61      REAL,INTENT(IN) :: pref(L_PINT)
     62      REAL,INTENT(IN) :: PW, TW
    5963
    60       integer MT, MP, N, M, NP
    61       real*8  PW, TW
    62       real*8  PWL, LCOEF(4), T, U
     64      ! OUTPUTS :
     65      !==========
     66      INTEGER,INTENT(OUT) :: MT, MP
     67      REAL,INTENT(OUT) :: LCOEF(4)
    6368
    64 C======================================================================C
    65  
    66 !     Get the upper and lower temperature grid indicies that bound the
    67 !     requested temperature. If the requested temperature is outside
    68 !     the T-grid, set up to extrapolate from the appropriate end.
    69 !     TW : temperature to be interpolated
    70 !     TREF : grid array
    71 !     MT : index of TREF for bounding new temperature
    72 !     U : new index (real) for temperature interpolated
     69      ! LOCAL :
     70      !========
     71      INTEGER :: N
     72      REAL :: PWL ! local value for log10(PW)
     73      REAL :: TWL ! local value for TW
     74      REAL :: T ! linear interpolation weight along log(pressure)
     75      REAL :: U ! linear interpolation weight along temperature
    7376
    74       IF(TW.LE.TREF(1)) THEN
     77
     78!     1. For the interpolation along temperature
     79!        Get the upper and lower temperature grid indicies that bound the
     80!        requested temperature. If the sought temperature is below
     81!        the T-grid minimum, assume that lower value. If the sought temperature
     82!        is above the T-grid maximum, assume that maximum value.
     83!        TW : temperature to interpolate to
     84!        TREF : reference temperature array
     85!        MT : index of TREF for (lower) bounding temperature
     86!        U : weight for interpolation in temperature
     87
     88      TWL = TW
     89      IF(TWL.LT.TREF(1)) THEN
     90
     91        ! Below lowest tabulated temperature
    7592        MT = 1
    76         IF (TW.LT.TREF(1)) THEN
    77          write(*,*) 'tpindex: Caution! Temperature of upper levels lower
    78      $ than ref temperature for k-coef: k-coeff fixed for upper levels'
    79          write(*,*) "         TW=",TW
    80          write(*,*) "         TREF(1)=",TREF(1)
    81         ENDIF
    82       ELSE
     93        TWL=TREF(1)*1.00
     94      ELSEIF(TWL.GE.TREF(L_NTREF)) THEN
     95       ! Above highest tabulated temperature
     96        MT=L_NTREF-1
     97        TWL = TREF(L_NTREF)*1.00
     98      ELSE
     99        ! Regular case, find encompassing tabulated temperatures
    83100        do n=1,L_NTREF-1
    84           if(tw.gt.Tref(n) .and. TW.LE.TREF(N+1)) then
     101          IF (TWL.GE.TREF(N) .and. TWL.LT.TREF(N+1)) THEN
    85102            MT = n
    86             goto 10
    87           end if
    88         end do
     103            EXIT
     104          END IF
     105        END DO
     106     
     107      END IF
     108!     weight for the temperature interpolation
     109      U = (TWL-TREF(MT))/(TREF(MT+1)-TREF(MT))
    89110
    90         MT = L_NTREF-1
    91      
    92    10   continue
    93       END IF
    94111
    95       !TB15 : case low temp : MT=1: fixed TW right above tref(1)
    96       IF (MT.eq.1) THEN
    97          TW=tref(1)*1.00
    98 !         write(*,*) 'tpindex: Caution! Temperature of upper levels lower
    99 !     $than ref temperature for k-coef: k-coeff fixed for upper levels'
    100 !         write(*,*) "         TW=",TW
    101 !         write(*,*) "         TREF(1)=",TREF(1)
     112!     2. For interpolation along pressure
     113!        Get the upper and lower pressure grid indicies that bound the
     114!        requested pressure. If the sought pressure is below
     115!        the P-grid minimum, assume that lower value. If the sought pressure
     116!        is above the P-grid maximum, assume that maximum value.
     117
     118      PWL = log10(PW)
     119      IF(PWL.LT.PREF(1)) THEN
     120        ! Below lowest tabulated pressure
     121        MP = 1
     122        PWL=Pref(1)*1.00
     123      ELSEIF(PWL.GE.PREF(L_PINT)) THEN
     124        ! Above highest tabulated pressure
     125        MP = L_PINT-1
     126        PWL = PREF(L_PINT)*1.00
     127      ELSE
     128        ! Regular case, find encompassing tabulated pressures
     129        do n=1,L_PINT-1
     130          IF (PWL.GE.PREF(N) .and. PWL.LT.PREF(N+1)) THEN
     131            MP = n
     132            EXIT
     133          END IF
     134        END DO
    102135      ENDIF
    103136
    104       U = (TW-TREF(MT))/(TREF(MT+1)-TREF(MT))
    105 
    106 !     Get the upper and lower pressure grid indicies that bound the
    107 !     requested pressure. If the requested pressure is outside
    108 !     the P-grid, set up to extrapolate from the appropriate end.
    109 
    110       pwl = log10(pw)
    111 
    112       do n=2,L_PINT-1
    113         if(pwl.le.Pref(n)) then
    114           MP = n-1
    115           goto 20
    116         end if
    117       end do
    118 
    119       MP = L_PINT-1
    120 
    121    20 continue
    122      
    123       !TB15 : case low pressure : n=2 : fixed pwl, right above pref(1)
    124       IF (MP.eq.1) THEN
    125         IF (PWL.LT.PREF(1)) THEN
    126          write(*,*) 'tpindex: Caution! Pressure of upper levels lower
    127      $than ref pressure for k-coef: k-coeff fixed for upper levels'
    128          write(*,*) "         PWL=",PWL
    129          write(*,*) "         PREF(1)=",PREF(1)
    130         ENDIF
    131         PWL=Pref(1)*1.00
    132       ENDIF
    133 
    134 !     interpolated pressure
     137!     weight for the pressure interpolation
    135138      T = (PWL-PREF(MP))/(PREF(MP+1)-PREF(MP))
    136139
    137 !  Fill in the interpolation coefficients
     140!  Build the interpolation coefficients (bilinear in temperature-log(pressure))
    138141      LCOEF(1) = (1.0-T)*(1.0-U)
    139142      LCOEF(2) = T*(1.0-U)
    … …  
    141144      LCOEF(4) = (1.0-T)*U
    142145
    143    30 CONTINUE
    144 
    145       return
    146       end
     146      end subroutine tpindex

  • trunk/LMDZ.TITAN/libf/phytitan/tracer_h.F90

    r2373 r3083  
    1 ! WARNING (JB 27/03/2018)
     1! AUTHOR : J. Burgalat (27/03/2018)
     2! NOTES : Add gazs_indx by B. de Batz de Trenquelléon (31/05/2022)
    23!
    34! OpenMP directives in this module are inconsistent:
    … …  
    3839  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: micro_indx  !! Indexes of all microphysical tracers
    3940  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: ices_indx   !! Indexes of all ice microphysical tracers
     41  INTEGER, DIMENSION(:), ALLOCATABLE, SAVE :: gazs_indx   !! Indexes of all gazs microphysical tracers
    4042
    4143  CHARACTER(len=30), DIMENSION(:), ALLOCATABLE, SAVE :: noms      !! name of the tracer
    … …  
    294296          IF (ANY(ices_indx == idx)) suffix=suffix//", ice"
    295297        ENDIF
     298        IF (ALLOCATED(gazs_indx)) THEN
     299          IF (ANY(gazs_indx == idx)) suffix=suffix//", gaz"
     300        ENDIF
    296301        suffix=suffix//")"
    297302      ELSE
    298303        suffix=" ()"
    299304      ENDIF
    300       WRITE(*,'(I5,(a),I6,(a))') i," -> ",idx ," : "//TRIM(noms(i))//suffix
     305      WRITE(*,'(I5,(a),I6,(a))') i," -> ",idx ," : "//TRIM(noms(idx))//suffix
    301306    ENDDO
    302307  END SUBROUTINE dumpTracers
Note: See TracChangeset for help on using the changeset viewer.