Changeset 1648 for trunk/LMDZ.TITAN/libf


Ignore:
Timestamp:
Jan 19, 2017, 2:46:53 PM (8 years ago)
Author:
jvatant
Message:

Modifications to custom radiative transfer to Titan
+ Enables an altitude dependant gfrac for CIA computations

-> many radical changes in su_gases and co ..
-> read vertical CH4 profile with call_profilgases
-> Now you need a 'profile.def' that I will add in the deftank

+ Added interpolate CIA routines for CH4
+ Added temporary mean aerosol profile opacity routine (disr_haze)

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

Legend:

Unmodified
Added
Removed

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

    r1647 r1648  
    3434      do igas=1,ngasmx
    3535
    36          if(igas.eq.vgas)then
    37             ! ignore variable gas in cpp calculation
     36         ! all values at 300 K from Engineering Toolbox
     37         if(igas.eq.igas_N2)then
     38            mugaz_c = mugaz_c + 28.01*gfrac(igas,nivref)
     39         elseif(igas.eq.igas_H2)then
     40            mugaz_c = mugaz_c + 2.01*gfrac(igas,nivref)
     41         elseif(igas.eq.igas_CH4)then
     42            mugaz_c = mugaz_c + 16.04*gfrac(igas,nivref)
    3843         else
    39             ! all values at 300 K from Engineering Toolbox
    40             if(igas.eq.igas_N2)then
    41                mugaz_c = mugaz_c + 28.01*gfrac(igas)
    42             elseif(igas.eq.igas_H2)then
    43                mugaz_c = mugaz_c + 2.01*gfrac(igas)
    44             elseif(igas.eq.igas_CH4)then
    45                mugaz_c = mugaz_c + 16.04*gfrac(igas)
    46             elseif(igas.eq.igas_C2H6)then
    47                ! C2H6 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
    48                mugaz_c = mugaz_c + 30.07*gfrac(igas)
    49             elseif(igas.eq.igas_C2H2)then
    50                ! C2H2 http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
    51                mugaz_c = mugaz_c + 26.04*gfrac(igas)
    52             else
    53                print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
    54                call abort
    55             endif
     44            print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
     45            call abort
    5646         endif
    57 
    5847      enddo
    5948
    … …  
    6150      do igas=1,ngasmx
    6251
    63          if(igas.eq.vgas)then
    64             ! ignore variable gas in cpp calculation
     52         ! all values at 300 K from Engineering Toolbox
     53         if(igas.eq.igas_N2)then
     54            cpp_c   = cpp_c   + 1.040*gfrac(igas,nivref)*28.01/mugaz_c
     55         elseif(igas.eq.igas_H2)then
     56            cpp_c   = cpp_c   + 14.31*gfrac(igas,nivref)*2.01/mugaz_c
     57         elseif(igas.eq.igas_CH4)then
     58            cpp_c   = cpp_c   + 2.226*gfrac(igas,nivref)*16.04/mugaz_c
    6559         else
    66             ! all values at 300 K from Engineering Toolbox
    67             if(igas.eq.igas_N2)then
    68                cpp_c   = cpp_c   + 1.040*gfrac(igas)*28.01/mugaz_c
    69             elseif(igas.eq.igas_H2)then
    70                cpp_c   = cpp_c   + 14.31*gfrac(igas)*2.01/mugaz_c
    71             elseif(igas.eq.igas_CH4)then
    72                cpp_c   = cpp_c   + 2.226*gfrac(igas)*16.04/mugaz_c
    73             elseif(igas.eq.igas_C2H6)then
    74                ! C2H6
    75                ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=28
    76                cpp_c   = cpp_c   + 1.763*gfrac(igas)*30.07/mugaz_c
    77             elseif(igas.eq.igas_C2H2)then
    78                ! C2H2
    79                ! http://encyclopedia.airliquide.com/Encyclopedia.asp?GasID=1
    80                cpp_c   = cpp_c   + 1.575*gfrac(igas)*26.04/mugaz_c
    81             else
    82                print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
    83                call abort
    84             endif
     60            print*,'Error in calc_cpp_mugaz: Gas species not recognised!'
     61            call abort
    8562         endif
    86 
    8763      enddo
    88 
    89 
    90 
    91 
    9264
    9365      cpp_c = 1000.0*cpp_c

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

    r1647 r1648  
    2626
    2727      use radinc_h, only: L_NSPECTV
    28       use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, taurayvar, scalep
     28      use radcommon_h, only: WAVEV, BWNV, DWNV, tstellar, tauray, scalep
    2929      use gases_h
    3030      use comcstfi_mod, only: g, mugaz
    … …  
    3838
    3939      logical typeknown
    40       real*8 tauvar,tauvarvar,tausum,tausumvar,tauwei,tauweivar,bwidth,bstart
     40      real*8 tauvar,tausum,tauwei,bwidth,bstart
    4141      double precision df
    4242
    … …  
    5252      typeknown=.false.
    5353      do igas=1,ngasmx
    54          if(igas.eq.vgas)then
    55             print*,'variable gas is ',trim(gnom(igas)),' in Rayleigh scattering '
    56          endif
    57          if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
     54         if(gfrac(igas,nivref).lt.5.e-2)then
    5855            print*,'Ignoring ',trim(gnom(igas)),' in Rayleigh scattering '// &
    5956            'as its mixing ratio is less than 0.05.'
    … …  
    7471            endif
    7572
    76             if((gfrac(igas).eq.1.0).and.(vgas.eq.0))then
     73            if(gfrac(igas,nivref).eq.1.0)then
    7774               print*,'Rayleigh scattering is for a pure ',trim(gnom(igas)),' atmosphere.'
    7875               typeknown=.true.
    … …  
    9289         tausum = 0.0
    9390         tauwei = 0.0
    94          tausumvar = 0.0
    95          tauweivar = 0.0
    9691         bstart = 10000.0/BWNV(N+1)
    9792         bwidth = (10000.0/BWNV(N)) - (10000.0/BWNV(N+1))
    … …  
    10095
    10196            tauvar=0.0
    102             tauvarvar=0.0
    10397            do igas=1,ngasmx
    104                if((igas/=vgas).and.(gfrac(igas).lt.5.e-2))then
     98               if(gfrac(igas,nivref).lt.5.e-2)then
    10599                  ! ignore variable gas in Rayleigh calculation
    106100                  ! ignore gases of mixing ratio < 0.05 in Rayleigh calculation
    … …  
    117111               endif
    118112
    119                if(igas.eq.vgas) then
    120                   tauvarvar=tauvarvar+tauconsti(igas)*tauvari(igas)
    121                   tauvar=tauvar+0.0*0.0*gfrac(igas)
    122                else
    123                   tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas)
    124                endif
     113               tauvar=tauvar+tauconsti(igas)*tauvari(igas)*gfrac(igas,nivref)
    125114
    126115            enddo
    … …  
    130119            tauwei=tauwei+df
    131120            tausum=tausum+tauvar*df
    132             tauweivar=tauweivar+df
    133             tausumvar=tausumvar+tauvarvar*df
    134121         
    135122         enddo
    136123         TAURAY(N)=tausum/tauwei
    137          TAURAYVAR(N)=tausumvar/tauweivar
    138124         ! we multiply by scalep here (100) because plev, which is used in optcv,
    139125         ! is in units of mBar, so we need to convert.

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

    r1647 r1648  
    11      subroutine callcorrk(ngrid,nlayer,pq,nq,qsurf,           &
    22          albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,    &
    3           tsurf,fract,dist_star,aerosol,muvar,                 &
     3          tsurf,fract,dist_star,aerosol,                       &
    44          dtlw,dtsw,fluxsurf_lw,                               &
    55          fluxsurf_sw,fluxsurfabs_sw,fluxtop_lw,               &
    … …  
    1717      USE tracer_h
    1818      use comcstfi_mod, only: pi, mugaz, cpp
    19       use callkeys_mod, only: diurnal,tracer,nosurf,satval,        &
     19      use callkeys_mod, only: diurnal,tracer,nosurf,        &
    2020                              strictboundcorrk,specOLR
    2121
    … …  
    5959      REAL,INTENT(IN) :: fract(ngrid)              ! Fraction of day.
    6060      REAL,INTENT(IN) :: dist_star                 ! Distance star-planet (AU).
    61       REAL,INTENT(IN) :: muvar(ngrid,nlayer+1)
    6261      logical,intent(in) :: firstcall              ! Signals first call to physics.
    6362      logical,intent(in) :: lastcall               ! Signals last call to physics.
    … …  
    133132      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
    134133      real*8 taugsurfi(L_NSPECTI,L_NGAUSS-1)
    135       real*8 qvar(L_LEVELS)                    ! Mixing ratio of variable component (mol/mol).
    136134
    137135      ! Local aerosol optical properties for each column on RADIATIVE grid.
    … …  
    153151      character(len=10) :: tmp2
    154152
    155       ! For fixed water vapour profiles.
    156       integer i_var
    157       real RH
    158       real*8 pq_temp(nlayer)
    159 ! real(KIND=r8) :: pq_temp(nlayer) ! better F90 way.. DOESNT PORT TO F77!!!
    160       real ptemp, Ttemp, qsat
    161 
    162153      logical OLRz
    163154      real*8 NFLUXGNDV_nu(L_NSPECTV)
    164 
    165       ! Included by RW for runaway greenhouse 1D study.
    166       real vtmp(nlayer)
    167       REAL*8 muvarrad(L_LEVELS)
    168      
     155   
    169156      ! Included by MT for albedo calculations.     
    170157      REAL*8 albedo_temp(L_NSPECTV) ! For equivalent albedo calculation.
    … …  
    459446      else
    460447         acosz=mu0(ig) ! Cosine of sun incident angle : 3D simulations or local 1D simulations using latitude.
    461       endif
    462 
    463 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    464 !!! Note by JL13 : In the following, some indices were changed in the interpolations,
    465 !!!                so that the model results are less dependent on the number of layers !
    466 !!!
    467 !!!           ---  The older versions are commented with the comment !JL13index  ---
    468 !!!
    469 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!
    470 
    471 
    472      do k=1,L_LEVELS
    473          qvar(k) = 1.0D-7
    474      end do
    475 
    476 
    477        ! IMPORTANT: Now convert from kg/kg to mol/mol.
    478 !       do k=1,L_LEVELS
    479 !          qvar(k) = qvar(k)/(epsi+qvar(k)*(1.-epsi))
    480 !       end do
    481 
    482        DO l=1,nlayer
    483           muvarrad(2*l)   = muvar(ig,nlayer+2-l)
    484           muvarrad(2*l+1) = (muvar(ig,nlayer+2-l)+muvar(ig,max(nlayer+1-l,1)))/2
    485        END DO
    486      
    487        muvarrad(1) = muvarrad(2)
    488        muvarrad(2*nlayer+1)=muvar(ig,1)         
    489      
    490       ! Keep values inside limits for which we have radiative transfer coefficients !!!
    491       if(L_REFVAR.gt.1)then ! (there was a bug here)
    492          do k=1,L_LEVELS
    493             if(qvar(k).lt.wrefvar(1))then
    494                qvar(k)=wrefvar(1)+1.0e-8
    495             elseif(qvar(k).gt.wrefvar(L_REFVAR))then
    496                qvar(k)=wrefvar(L_REFVAR)-1.0e-8
    497             endif
    498          end do
    499448      endif
    500449
    … …  
    574523            endif
    575524         elseif(tlevrad(k).gt.tgasmax)then
    576             print*,'Maximum temperature is outside the radiative'
    577             print*,'transfer kmatrix bounds, exiting.'
    578             print*,"k=",k," tlevrad(k)=",tlevrad(k)
    579             print*,"tgasmax=",tgasmax
     525!            print*,'Maximum temperature is outside the radiative'
     526!            print*,'transfer kmatrix bounds, exiting.'
     527!            print*,"k=",k," tlevrad(k)=",tlevrad(k)
     528!            print*,"tgasmax=",tgasmax
    580529            if (strictboundcorrk) then
    581530              call abort
    582531            else
    583               print*,'***********************************************'
    584               print*,'we allow model to continue with tlevrad=tgasmax' 
    585               print*,'  ... we assume we know what you are doing ... '
    586               print*,'  ... but do not let this happen too often ... '
    587               print*,'***********************************************'
     532!              print*,'***********************************************'
     533!              print*,'we allow model to continue with tlevrad=tgasmax' 
     534!              print*,'  ... we assume we know what you are doing ... '
     535!              print*,'  ... but do not let this happen too often ... '
     536!              print*,'***********************************************'
    588537              !tlevrad(k)=tgasmax
    589538            endif
    … …  
    607556            endif
    608557         elseif(tmid(k).gt.tgasmax)then
    609             print*,'Maximum temperature is outside the radiative'
    610             print*,'transfer kmatrix bounds, exiting.'
    611             print*,"k=",k," tmid(k)=",tmid(k)
    612             print*,"tgasmax=",tgasmax
     558!            print*,'Maximum temperature is outside the radiative'
     559!            print*,'transfer kmatrix bounds, exiting.'
     560!            print*,"k=",k," tmid(k)=",tmid(k)
     561!            print*,"tgasmax=",tgasmax
    613562            if (strictboundcorrk) then
    614563              call abort
    615564            else
    616               print*,'***********************************************'
    617               print*,'we allow model to continue with tmid=tgasmin'
    618               print*,'  ... we assume we know what you are doing ... '
    619               print*,'  ... but do not let this happen too often ... '
    620               print*,'***********************************************'
     565!              print*,'***********************************************'
     566!              print*,'we allow model to continue with tmid=tgasmin'
     567!              print*,'  ... we assume we know what you are doing ... '
     568!              print*,'  ... but do not let this happen too often ... '
     569!              print*,'***********************************************'
    621570              tmid(k)=tgasmax
    622571            endif
    … …  
    649598            call optcv(dtauv,tauv,taucumv,plevrad,                 &
    650599                 qxvaer,qsvaer,gvaer,wbarv,cosbv,tauray,tauaero,   &
    651                  tmid,pmid,taugsurf,qvar,muvarrad)
     600                 tmid,pmid,taugsurf,gweight)
    652601
    653602            call sfluxv(dtauv,tauv,taucumv,albv,dwnv,wbarv,cosbv,  &
    … …  
    692641         call optci(plevrad,tlevrad,dtaui,taucumi,                  &
    693642              qxiaer,qsiaer,giaer,cosbi,wbari,tauaero,tmid,pmid,    &
    694               taugsurfi,qvar,muvarrad)
     643              taugsurfi,gweight)
    695644
    696645         call sfluxi(plevrad,tlevrad,dtaui,taucumi,ubari,albi,      &
    … …  
    813762        IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
    814763        IF( ALLOCATED( tgasref ) ) DEALLOCATE( tgasref )
    815         IF( ALLOCATED( wrefvar ) ) DEALLOCATE( wrefvar )
    816764        IF( ALLOCATED( pfgasref ) ) DEALLOCATE( pfgasref )
    817765!$OMP END MASTER

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

    r1647 r1648  
    6060      real,save :: size_tropo
    6161      real,save :: size_strato
    62       real,save :: satval
    63 !$OMP THREADPRIVATE(size_tropo,size_strato,satval)
     62!$OMP THREADPRIVATE(size_tropo,size_strato)
    6463      real,save :: pceil
    6564!$OMP THREADPRIVATE(pceil)

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

    r1647 r1648  
    33      implicit none
    44
    5 !======================================================================C
     5!============================================================================================C
    66!
    77!     gases_h   
    88!
    9 !     A F90-allocatable version for gases.h -- AS 12/2011
     9!     * A F90-allocatable version for gases.h -- AS 12/2011
    1010!
    11 !======================================================================C
     11!     * Titan's version : J. Vatant d'Ollone (2017)
     12!                  * gfrac is now 2D (altitude-dependent for the CIA-calculations)
     13!                  * Reference level (nivref) is added for the cpp_mu and rayleigh routines
     14!===========================================================================================C
    1215
    1316      ! Set and allocated in su_gases.F90
    1417      integer :: ngasmx
    15       integer :: vgas
    16       character*20,allocatable,DIMENSION(:) :: gnom ! name of the gas, read by master
    17       real,allocatable,DIMENSION(:) :: gfrac
     18      integer :: nivref
     19      character*20,allocatable,DIMENSION(:) :: gnom ! name of the gas, set by master
     20      real,allocatable,DIMENSION(:,:) :: gfrac
    1821
    1922      ! in analogy with tracer.h ...
    20       integer :: igas_H2
     23     
    2124      integer :: igas_N2
    2225      integer :: igas_CH4
    23       integer :: igas_C2H2
    24       integer :: igas_C2H6
    25 !!$OMP THREADPRIVATE(ngasmx,vgas,gnom,gfrac,&
    26 !       !$OMP igas_H2,igas_N2,igas_CH4,igas_C2H2,igas_C2H6)
     26      integer :: igas_H2
     27     
     28!!$OMP THREADPRIVATE(ngasmx,nivref,gnom,gfrac,&
     29!       !$OMP igas_H2,igas_N2,igas_CH4)
    2730
    2831      end module gases_h

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

    r1647 r1648  
    439439
    440440!=================================
    441 
    442      write(*,*)"What is the saturation % of the variable species?"
    443      satval=0.8
    444      call getin_p("satval",satval)
    445      write(*,*)" satval = ",satval
    446441
    447442     write(*,*) "Gravitationnal sedimentation ?"
    … …  
    477472!           mugaz=8.314*1000./pr
    478473     endif ! of if (force_cpp)
    479      call su_gases
     474     
     475     
     476     call su_gases(nlayer,tracer)     
    480477     call calc_cpp_mugaz
    481 
     478     
     479     
    482480     PRINT*,'--------------------------------------------'
    483481     PRINT*
    … …  
    519517  ENDDO
    520518
    521   ! initialize variables in radinc_h
    522519  call ini_radinc_h(nlayer)
    523520 

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

    r1315 r1648  
    7878        write(*,*) 'is correct. You can change it in callphys.def with:'
    7979        write(*,*) 'datadir = /absolute/path/to/datagcm'
    80         write(*,*) 'Also check that the continuum data continuum_data/N2-N2_norm_2011.cia is there.'
     80        write(*,*) 'Also check that the continuum data continuum_data/N2-N2_2011.cia is there.'
    8181        call abort
    8282     else

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

    r1647 r1648  
    11subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
    22     QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
    3      TMID,PMID,TAUGSURF,QVAR,MUVAR)
     3     TMID,PMID,TAUGSURF,GWEIGHT)
    44
    55  use radinc_h
    6   use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
     6  use radcommon_h, only: gasi,tlimit,Cmk,tgasref,pfgasref,wnoi,scalep,indi,glat_ig
    77  use gases_h
    8   use comcstfi_mod, only: g, r, mugaz
     8  use comcstfi_mod, only: g, r
    99  use callkeys_mod, only: continuum,graybody
    1010  implicit none
    … …  
    4949  real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
    5050
     51  ! Titan customisation
     52  ! J. Vatant d'Ollone (2016)
     53  real*8 GWEIGHT(L_NGAUSS)
     54  real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
     55  real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
     56  real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
     57  real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
     58  real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
     59  real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
     60
     61  CHARACTER*2  str2
     62  ! ==========================
     63
    5164  integer L, NW, NG, K, LK, IAER
    5265  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
    … …  
    6073  double precision p_cross
    6174
    62   ! variable species mixing ratio variables
    63   real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    6475  real*8  KCOEF(4)
    65   integer NVAR(L_LEVELS)
    66 
     76 
    6777  ! temporary variables for multiple aerosol calculation
    6878  real*8 atemp
    … …  
    7383  !real*8 rho !! see test below
    7484
    75   integer igas, jgas
     85  integer igas, jgas, ilay
    7686
    7787  integer interm
    … …  
    95105
    96106     ! if we have continuum opacities, we need dz
    97       dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
    98       U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
    99             !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
    100 
    101      call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
    102           LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
     107      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))
     108      U(k)  = Cmk*DPR(k)     ! only Cmk line in optci.F
     109     
     110     call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
    103111
    104112     do LK=1,4
    … …  
    120128     do K=2,L_LEVELS
    121129
     130        ilay = k / 2 ! int. arithmetic => gives the gcm layer index
     131       
    122132! continuum absorption
    123133        DCONT = 0.0d0
    … …  
    129139           do igas=1,ngasmx
    130140
    131               if(gfrac(igas).eq.-1)then ! variable
    132                  p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
    133               else
    134                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
    135               endif
     141              p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
    136142
    137143              dtemp=0.0d0
    … …  
    151157                 ! then cross-interactions with other gases
    152158                 do jgas=1,ngasmx
    153                     p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
     159                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
    154160                    dtempc  = 0.0d0
    155161                    if(jgas.eq.igas_N2)then
    156162                       interm = indi(nw,igas,jgas)
    157163                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     164                       indi(nw,igas,jgas) = interm
     165                    endif
     166                    dtemp = dtemp + dtempc
     167                 enddo
     168
     169               elseif(igas.eq.igas_CH4)then
     170
     171                 ! first do self-induced absorption
     172                 interm = indi(nw,igas,igas)
     173                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     174                 indi(nw,igas,igas) = interm
     175
     176                 ! then cross-interactions with other gases
     177                 do jgas=1,ngasmx
     178                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
     179                    dtempc  = 0.0d0
     180                    if(jgas.eq.igas_N2)then
     181                       interm = indi(nw,igas,jgas)
     182                       call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    158183                       indi(nw,igas,jgas) = interm
    159184                    endif
    … …  
    185210        DAERO=SUM(TAEROS(K,NW,1:naerkind))
    186211
     212!================= Titan customisation ========================================
     213        call disr_haze(dz(k),plev(k),wnoi(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
     214! =============================================================================
     215
     216
    187217        do ng=1,L_NGAUSS-1
    188218
    189219           ! Now compute TAUGAS
    190220
    191            ! Interpolate between water mixing ratios
    192            ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
    193            ! the water data range
    194 
    195            if(L_REFVAR.eq.1)then ! added by RW for special no variable case
    196               KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
    197               KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
    198               KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
    199               KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
    200            else
    201 
    202               KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*     &
    203                    (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   &
    204                    GASI(MT(K),MP(K),NVAR(K),NW,NG))
    205 
    206               KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*   &
    207                    (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                 &
    208                    GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
    209 
    210               KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
    211                    (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -               &
    212                    GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
    213 
    214               KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*   &
    215                    (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
    216                    GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
    217 
    218            endif
     221           KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
     222           KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
     223           KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
     224           KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
    219225
    220226           ! Interpolate the gaseous k-coefficients to the requested T,P values
    … …  
    228234           DTAUKI(K,nw,ng) = TAUGAS    &
    229235                             + DCONT   & ! For parameterized continuum absorption
    230                              + DAERO     ! For aerosol absorption
     236                             + DAERO   & ! For aerosol absorption
     237                             + DHAZE_T(K,NW)  ! For Titan haze
    231238
    232239        end do
    … …  
    238245        DTAUKI(K,nw,ng) = 0.d0      &
    239246                          + DCONT   & ! For parameterized continuum absorption
    240                           + DAERO     ! For aerosol absorption
     247                          + DAERO   & ! For aerosol absorption
     248                          + DHAZE_T(K,NW)     ! For Titan Haze
    241249
    242250     end do
    … …  
    244252
    245253  DTAUKI(L_LEVELS+1,1:L_NSPECTI,1:L_NGAUSS)=0.d0
     254 
     255 
    246256  !=======================================================================
    247257  !     Now the full treatment for the layers, where besides the opacity
    248258  !     we need to calculate the scattering albedo and asymmetry factors
     259  ! ======================================================================
    249260
    250261  do iaer=1,naerkind
    … …  
    256267  end do
    257268 
     269  ! Haze scattering
    258270  DO NW=1,L_NSPECTI
    259      DO L=1,L_NLAYRAD
     271    DO K=2,L_LEVELS+1
     272      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
     273    ENDDO
     274  ENDDO
     275
     276  DO NW=1,L_NSPECTI
     277     DO L=1,L_NLAYRAD-1
    260278        K              = 2*L+1
    261         btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
     279        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
     280                      + DHAZES_T(K,NW) + DHAZES_T(K+1,NW)
    262281     END DO ! L vertical loop
     282     
     283     !last level
     284     L              = L_NLAYRAD
     285     K              = 2*L+1
     286     
     287     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + DHAZES_T(K,NW)
     288     
    263289  END DO                    ! NW spectral loop
    264290 
     291! ======================================================================
    265292
    266293  DO NW=1,L_NSPECTI
    267294     NG = L_NGAUSS
    268      DO L=1,L_NLAYRAD
     295     DO L=1,L_NLAYRAD-1
    269296
    270297        K              = 2*L+1
     298
    271299        DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
    272300
    … …  
    278306                   GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
    279307           end do
     308           atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
    280309           WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
    281310        else
    … …  
    291320
    292321     END DO ! L vertical loop
     322     
     323     !     No vertical averaging on bottom layer
     324
     325     L = L_NLAYRAD
     326     K = 2*L+1
     327     DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
     328     
     329        atemp = 0.
     330        if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
     331           do iaer=1,naerkind
     332              atemp = atemp +                                     &
     333                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER)
     334           end do
     335           atemp = atemp + ASF_T(K,NW)*DHAZES_T(K,NW)
     336           WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
     337        else
     338           WBARI(L,nw,ng) = 0.0D0
     339           DTAUI(L,NW,NG) = 1.0D-9
     340        endif
     341
     342        if(btemp(L,nw) .GT. 0.0d0) then
     343           cosbi(L,NW,NG) = atemp/btemp(L,nw)
     344        else
     345           cosbi(L,NW,NG) = 0.0D0
     346        end if
    293347
    294348     ! Now the other Gauss points, if needed.
    295349
    296350     DO NG=1,L_NGAUSS-1
     351     
    297352        IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
    298353
    299            DO L=1,L_NLAYRAD
     354           DO L=1,L_NLAYRAD-1
    300355              K              = 2*L+1
    301356              DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
    302357
    303358              if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
    304 
    305                  WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
    306 
     359                WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
    307360              else
    308361                 WBARI(L,nw,ng) = 0.0D0
    … …  
    312365              cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
    313366           END DO ! L vertical loop
     367           
     368            !     No vertical averaging on bottom layer
     369
     370            L = L_NLAYRAD
     371            K = 2*L+1
     372            DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)
     373            if(DTAUI(L,NW,NG) .GT. 1.0D-9) then
     374              WBARI(L,nw,ng) = btemp(L,nw)  / DTAUI(L,NW,NG)
     375            else
     376               WBARI(L,nw,ng) = 0.0D0
     377               DTAUI(L,NW,NG) = 1.0D-9
     378            endif
     379            cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
     380           
    314381        END IF
    315382
    … …  
    334401  ! ascending ray with angle theta = 0.
    335402
    336   !open(127,file='taucum.out')
    337   !do nw=1,L_NSPECTI
    338   !   write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
    339   !enddo
    340   !close(127)
    341  
    342 !  print*,'WBARI'
    343 !  print*,WBARI
    344 !  print*,'DTAUI'
    345 !  print*,DTAUI
    346 !  call abort
    347  
     403
     404!  Titan's outputs (J.V.O, 2016)===============================================
     405!      do l=1,L_NLAYRAD
     406!         do nw=1,L_NSPECTI
     407!          INT_DTAU(L,NW) = 0.0d+0
     408!            DO NG=1,L_NGAUSS
     409!               INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtaui(L,nw,ng)*gweight(NG)
     410!            enddo
     411!         enddo
     412!      enddo
     413
     414!       do nw=1,L_NSPECTI
     415!          write(str2,'(i2.2)') nw
     416!         call writediagfi(1,'kgi'//str2,'Gaz extinction coefficient IR band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))
     417!          call writediagfi(1,'khi'//str2,'Haze extinction coefficient IR band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))       
     418!       enddo 
     419 
     420! ============================================================================== 
    348421
    349422  return

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

    r1647 r1648  
    11SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
    22     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
    3      TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
     3     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,GWEIGHT)
    44
    55  use radinc_h
    6   use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
     6  use radcommon_h, only: gasv, tlimit, Cmk, tgasref, pfgasref,wnov,scalep,indv,glat_ig
    77  use gases_h
    8   use comcstfi_mod, only: g, r, mugaz
     8  use comcstfi_mod, only: g, r
    99  use callkeys_mod, only: continuum,graybody,callgasvis
    1010
    … …  
    5555  real*8  TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
    5656
     57  ! Titan customisation
     58  ! J. Vatant d'Ollone (2016)
     59  real*8 GWEIGHT(L_NGAUSS)
     60  real*8 DHAZE_T(L_LEVELS+1,L_NSPECTI)
     61  real*8 DHAZES_T(L_LEVELS+1,L_NSPECTI)
     62  real*8 SSA_T(L_LEVELS+1,L_NSPECTI)
     63  real*8 ASF_T(L_LEVELS+1,L_NSPECTI)
     64  real*8 INT_DTAU(L_NLAYRAD,L_NSPECTI)
     65  real*8 K_HAZE(L_NLAYRAD,L_NSPECTI)
     66 
     67  CHARACTER*2  str2
     68  ! ==========================
     69
    5770  integer L, NW, NG, K, LK, IAER
    5871  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
    … …  
    6982  double precision p_cross
    7083
    71   ! variable species mixing ratio variables
    72   real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    7384  real*8  KCOEF(4)
    74   integer NVAR(L_LEVELS)
    7585
    7686  ! temporary variables for multiple aerosol calculation
    … …  
    8292  real*8 dz(L_LEVELS)
    8393
    84   integer igas, jgas
     94  integer igas, jgas, ilay
    8595
    8696  integer interm
    … …  
    107117     ! if we have continuum opacities, we need dz
    108118
    109       dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))*mugaz/muvar(k)
    110       U(k)  = Cmk*DPR(k)*mugaz/muvar(k)     ! only Cmk line in optci.F 
    111           !JL13 the mugaz/muvar factor takes into account water meanmolecular weight if water is present
    112      
    113 
    114      call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
    115           LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
     119      dz(k) = dpr(k)*R*TMID(K)/(glat_ig*PMID(K))
     120      U(k)  = Cmk*DPR(k)     ! only Cmk line in optcv.F     
     121
     122     call tpindex(PMID(K),TMID(K),pfgasref,tgasref,LCOEF,MT(K),MP(K))
    116123
    117124     do LK=1,4
    … …  
    133140     end do                    ! levels
    134141  end do
    135  
     142
    136143  !     we ignore K=1...
    137144  do K=2,L_LEVELS
     145 
     146     ilay = k / 2 ! int. arithmetic => gives the gcm layer index
    138147
    139148     do NW=1,L_NSPECTV
    … …  
    153162           do igas=1,ngasmx
    154163
    155               if(gfrac(igas).eq.-1)then ! variable
    156                  p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
    157               else
    158                  p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
    159               endif
     164              p_cont  = dble(PMID(k)*scalep*gfrac(igas,ilay))
    160165
    161166              dtemp=0.0
    … …  
    176181                 ! then cross-interactions with other gases
    177182                 do jgas=1,ngasmx
    178                     p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
     183                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
    179184                    dtempc  = 0.0
    180185                    if(jgas.eq.igas_N2)then
    … …  
    187192                 enddo
    188193
     194               elseif(igas.eq.igas_CH4)then
     195
     196                 ! first do self-induced absorption
     197                 interm = indv(nw,igas,igas)
     198                 call interpolateCH4CH4(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     199                 indv(nw,igas,igas) = interm
     200
     201                 ! then cross-interactions with other gases
     202                 do jgas=1,ngasmx
     203                    p_cross = dble(PMID(k)*scalep*gfrac(jgas,ilay))
     204                    dtempc  = 0.0
     205                    if(jgas.eq.igas_N2)then
     206                       interm = indv(nw,igas,jgas)
     207                       call interpolateN2CH4(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     208                       indv(nw,igas,jgas) = interm
     209                    endif
     210                    dtemp = dtemp + dtempc
     211                 enddo
     212
    189213              endif
    190214
    … …  
    197221        endif
    198222
     223!================= Titan customisation ========================================
     224        call disr_haze(dz(k),plev(k),wnov(nw),dhaze_T(k,nw),SSA_T(k,nw),ASF_T(k,nw))
     225! =============================================================================
     226
    199227        do ng=1,L_NGAUSS-1
    200228
    201229           ! Now compute TAUGAS
    202230
    203            ! Interpolate between water mixing ratios
    204            ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
    205            ! the water data range
    206 
    207            if(L_REFVAR.eq.1)then ! added by RW for special no variable case
    208               KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
    209               KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
    210               KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
    211               KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
    212            else
    213 
    214               KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
    215                    (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
    216                    GASV(MT(K),MP(K),NVAR(K),NW,NG))
    217 
    218               KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
    219                    (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
    220                    GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
    221 
    222               KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
    223                    (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
    224                    GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
    225 
    226               KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
    227                    (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
    228                    GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
    229 
    230            endif
     231           KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
     232           KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
     233           KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
     234           KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
    231235
    232236           ! Interpolate the gaseous k-coefficients to the requested T,P values
    … …  
    240244           DTAUKV(K,nw,ng) = TAUGAS &
    241245                             + TRAYAER & ! TRAYAER includes all scattering contributions
    242                              + DCONT ! For parameterized continuum aborption
     246                             + DCONT    & ! For parameterized continuum aborption
     247                             + DHAZE_T(K,NW)  ! For Titan haze
    243248
    244249        end do
    … …  
    248253
    249254        NG              = L_NGAUSS
    250         DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
     255        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT  & ! For parameterized continuum absorption
     256                              + DHAZE_T(K,NW)  ! For Titan haze
    251257
    252258        do iaer=1,naerkind
    … …  
    261267  !     Now the full treatment for the layers, where besides the opacity
    262268  !     we need to calculate the scattering albedo and asymmetry factors
    263 
     269  ! ======================================================================
     270 
    264271  do iaer=1,naerkind
    265272    DO NW=1,L_NSPECTV
    … …  
    269276    ENDDO
    270277  end do
     278 
     279  ! Haze scattering
     280  DO NW=1,L_NSPECTV
     281    DO K=2,L_LEVELS+1
     282      DHAZES_T(K,NW) = DHAZE_T(K,NW) * SSA_T(K,NW)
     283    ENDDO
     284  ENDDO
     285
    271286
    272287  DO NW=1,L_NSPECTV
    273288     DO L=1,L_NLAYRAD-1
    274289        K              = 2*L+1
    275         atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))+SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind))
    276         btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind))
     290       
     291        atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
     292                    + SUM(GVAER(K+1,NW,1:naerkind) * TAUAEROLK(K+1,NW,1:naerkind)) &
     293                    + ASF_T(K,NW)*DHAZES_T(K,NW)
     294                   
     295        btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) + SUM(TAUAEROLK(K+1,NW,1:naerkind)) &
     296                    + DHAZES_T(K,NW)
     297       
    277298        ctemp(L,NW) = btemp(L,NW) + 0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
    278299        btemp(L,NW) = btemp(L,NW) + TRAY(K,NW) + TRAY(K+1,NW)
    … …  
    283304     L              = L_NLAYRAD
    284305     K              = 2*L+1
    285      atemp(L,NW)    = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind))
    286      btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind))
     306     
     307     atemp(L,NW) = SUM(GVAER(K,NW,1:naerkind) * TAUAEROLK(K,NW,1:naerkind)) &
     308                 + ASF_T(K,NW)*DHAZES_T(K,NW)
     309     btemp(L,NW) = SUM(TAUAEROLK(K,NW,1:naerkind)) &
     310                 + DHAZES_T(K,NW)
    287311     ctemp(L,NW) = btemp(L,NW) + 0.9999*TRAY(K,NW)
    288312     btemp(L,NW) = btemp(L,NW) + TRAY(K,NW)
    … …  
    292316  END DO                    ! NW spectral loop
    293317
     318! ===========================================================================================
     319
    294320  DO NG=1,L_NGAUSS
    295321    DO NW=1,L_NSPECTV
    … …  
    309335
    310336        WBARV(L,NW,NG) = ctemp(L,NW) / DTAUV(L,NW,NG)
     337
    311338     END DO                 ! NW spectral loop
    312339  END DO                    ! NG Gauss loop
    … …  
    329356
    330357
     358!  Titan's outputs (J.V.O, 2016)===============================================
     359!      do l=1,L_NLAYRAD
     360!         do nw=1,L_NSPECTV
     361!          INT_DTAU(L,NW) = 0.0d+0
     362!            DO NG=1,L_NGAUSS
     363!               INT_DTAU(L,NW)= INT_DTAU(L,NW) + dtauv(L,nw,ng)*gweight(NG)
     364!            enddo
     365!         enddo
     366!      enddo
     367
     368!       do nw=1,L_NSPECTV
     369!          write(str2,'(i2.2)') nw
     370!         call writediagfi(1,'kgv'//str2,'Gaz extinction coefficient VI band '//str2,'m-1',1,int_dtau(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))
     371!          call writediagfi(1,'khv'//str2,'Haze extinction coefficient VI band '//str2,'m-1',1,k_haze(L_NLAYRAD:1:-1,nw)/dz_lay(L_NLAYRAD:1:-1))       
     372!       enddo 
     373
     374! ============================================================================== 
     375
     376
    331377  return
    332378

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

    r1647 r1648  
    1515 
    1616      use radinc_h, only : L_NSPECTI,L_NSPECTV,naerkind
    17       use gases_h, only: gnom, gfrac
    1817      use radcommon_h, only: sigma, glat, grav, BWNV
    1918      use surfdat_h, only: phisfi, zmea, zstd, zsig, zgam, zthe
    … …  
    2322      use geometry_mod, only: latitude, longitude, cell_area
    2423      USE comgeomfi_h, only: totarea, totarea_planet
    25       USE tracer_h, only: noms, mmol, radius, rho_q, qext, &
     24      USE tracer_h, only: noms, mmol, radius, qext, &
    2625                          alpha_lift, alpha_devil, qextrhor, &
    2726                          igcm_dustbin
    … …  
    353352!$OMP THREADPRIVATE(qsurf_hist)
    354353
    355       real muvar(ngrid,nlayer+1) ! For Runaway Greenhouse 1D study. By RW
    356 
    357354      real,allocatable,dimension(:,:,:),save :: reffrad  ! Aerosol effective radius (m). By RW
    358355      real,allocatable,dimension(:,:,:),save :: nueffrad ! Aerosol effective radius variance. By RW
    … …  
    660657! II.a Call correlated-k radiative transfer scheme
    661658! ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    662  
    663                muvar(1:ngrid,1:nlayer+1)=mugaz
     659
     660               call call_profilgases(nlayer)
    664661
    665662               ! standard callcorrk
    666663               call callcorrk(ngrid,nlayer,pq,nq,qsurf,                           &
    667664                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
    668                               tsurf,fract,dist_star,aerosol,muvar,                &
     665                              tsurf,fract,dist_star,aerosol,                      &
    669666                              zdtlw,zdtsw,fluxsurf_lw,fluxsurf_sw,                &
    670667                              fluxsurfabs_sw,fluxtop_lw,                          &

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

    r1529 r1648  
    77!
    88!                             radcommon.h
    9 !
     9!v
    1010!----------------------------------------------------------------------C
    1111!
    … …  
    3636!     TAURAY     - Array (NSPECTV elements) of the pressure-independent
    3737!                  part of Rayleigh scattering optical depth.
    38 !     TAURAYVAR  - Array (NSPECTV elements) of the pressure-independent
    39 !                  part of Rayleigh scattering optical depth for the variable gas.
    4038!     FZEROI     - Fraction of zeros in the IR CO2 k-coefficients, for
    4139!                  each temperature, pressure, and spectral interval
    … …  
    6361      REAL*8 BWNI(L_NSPECTI+1), WNOI(L_NSPECTI), DWNI(L_NSPECTI), WAVEI(L_NSPECTI) !BWNI read by master in setspi
    6462      REAL*8 BWNV(L_NSPECTV+1), WNOV(L_NSPECTV), DWNV(L_NSPECTV), WAVEV(L_NSPECTV) !BWNV read by master in setspv
    65       REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV), TAURAYVAR(L_NSPECTV)
     63      REAL*8 STELLARF(L_NSPECTV), TAURAY(L_NSPECTV)
    6664!$OMP THREADPRIVATE(WNOI,DWNI,WAVEI,&
    6765        !$OMP WNOV,DWNV,WAVEV,&
    68         !$OMP STELLARF,TAURAY,TAURAYVAR)
     66        !$OMP STELLARF,TAURAY)
    6967
    7068      REAL*8 blami(L_NSPECTI+1)
    … …  
    7977      !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 
    8078      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
    81       REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, WREFVAR, PFGASREF
     79      REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF
    8280      real*8 FZEROI(L_NSPECTI)
    8381      real*8 FZEROV(L_NSPECTV)
    8482      real*8 pgasmin, pgasmax
    8583      real*8 tgasmin, tgasmax
    86 !$OMP THREADPRIVATE(gasi,gasv,&  !wrefvar,pgasref,tgasref,pfgasref read by master in sugas_corrk
     84!$OMP THREADPRIVATE(gasi,gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
    8785        !$OMP FZEROI,FZEROV)     !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
    8886

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

    r1647 r1648  
    1 subroutine su_gases
     1subroutine su_gases(nlayer,tracer)
    22
    33  use gases_h
    … …  
    55  implicit none
    66
    7   integer igas, ierr, count
    8 
     7  integer,intent(in) :: nlayer
     8  logical,intent(in) :: tracer
     9 
     10  integer igas, ierr
     11 
    912  !==================================================================
    1013  !     
    … …  
    1720  !     R. Wordsworth (2011)
    1821  !     Allocatable arrays by A. Spiga (2011)
    19   !     
     22  !     Titan's version : J. Vatant d'Ollone (2017)
     23  !              * force igas_X labels and def gnom for N2,CH4,H2
     24  !              * get rid of variable species ( enrichment dimension will be set in corrk routines )
     25  !
    2026  !==================================================================
    2127
    2228!$OMP MASTER
    23   ! load gas names from file 'gases.def'
     29
     30
     31  ngasmx   = 3 ! N2, CH4, H2
     32
     33
     34  ! load reference level and reference molar fractions from file 'gases.def'
    2435  open(90,file='gases.def',status='old',form='formatted',iostat=ierr)
     36 
    2537  if (ierr.eq.0) then
    2638     write(*,*) "sugases.F90: reading file gases.def"
    27      read(90,*)
    28      read(90,*,iostat=ierr) ngasmx
     39     read(90,*) ! header
     40
     41     ! We allocate gfrac and we set gas molar fractions for the reference level only.
     42     ! This will be useful for the cpp_mu and rayleigh routines
     43     ! Other gas molar fractions are now set in routine callprofilgases
     44     
     45     write(*,*) 'sugases.F90: allocating and reading gas molar fractions from reference level in gases.def...'
     46     
     47     if(.not.allocated(gfrac)) allocate(gfrac(ngasmx,nlayer))
     48     
     49     read(90,*,iostat=ierr) nivref     
    2950     if (ierr.ne.0) then
    30         write(*,*) "sugases.F90: error reading number of gases"
     51        write(*,*) "sugases.F90: error reading reference level"
    3152        write(*,*) "   (first line of gases.def) "
    3253        call abort
    3354     endif
    34 
    35      print*,ngasmx, " gases found in gases.def. Allocating names and molar fractions..."
    36 
    37      if (.not.allocated(gnom)) allocate(gnom(ngasmx))
     55     
     56     print*, "layer", nivref, "is reference level found in gases.def ..."
     57     
    3858     do igas=1,ngasmx
    39         read(90,*,iostat=ierr) gnom(igas)
     59        read(90,*,iostat=ierr) gfrac(igas,nivref)
    4060        if (ierr.ne.0) then
    41            write(*,*) 'sugases.F90: error reading gas names in gases.def...'
     61           write(*,*) 'sugases.F90: error reading reference gas molar fractions in gases.def... aborting'
    4262           call abort
    4363        endif
    4464     enddo                  !of do igas=1,ngasmx
    4565
    46      vgas=0
    47      if(.not.allocated(gfrac)) allocate(gfrac(ngasmx))
    48      do igas=1,ngasmx
    49         read(90,*,iostat=ierr) gfrac(igas)
    50         if (ierr.ne.0) then
    51            write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...'
    52            call abort
    53         endif
    5466
    55         ! find variable gas (if any)
    56         if(gfrac(igas).eq.-1.0)then
    57            if(vgas.eq.0)then
    58               vgas=igas
    59            else
    60               print*,'You seem to be choosing two variable gases'
    61               print*,'Check that gases.def is correct'
    62               call abort
    63            endif
    64         endif
    65 
    66      enddo                  !of do igas=1,ngasmx
     67     ! We force gnom = (N2, CH4, H2) and igas_X for Titan
     68     
     69     igas_N2  = 1
     70     igas_CH4 = 2
     71     igas_H2  = 3
    6772
    6873
    69      ! assign the 'igas_X' labels
    70      count=0
    71      do igas=1,ngasmx
    72         if (trim(gnom(igas)).eq."H2_" .or. trim(gnom(igas)).eq."H2") then
    73            igas_H2=igas
    74            count=count+1
    75         elseif (trim(gnom(igas)).eq."N2_" .or. trim(gnom(igas)).eq."N2") then
    76            igas_N2=igas
    77            count=count+1
    78         elseif (trim(gnom(igas)).eq."CH4") then
    79            igas_CH4=igas
    80            count=count+1
    81         elseif (trim(gnom(igas)).eq."C2H6") then
    82            igas_C2H6=igas
    83            count=count+1
    84         elseif (trim(gnom(igas)).eq."C2H2") then
    85            igas_C2H2=igas
    86            count=count+1
    87         endif
    88      enddo
    89 
    90      if(count.ne.ngasmx)then
    91         print*,'Mismatch between ngas and number of recognised gases in sugas_corrk.F90.'
    92         print*,'Either we haven`t managed to assign all the gases, or there are duplicates.'
    93         print*,'Please try again.'
    94      endif
    95 
     74     if (.not.allocated(gnom)) allocate(gnom(ngasmx))
     75     gnom(igas_N2)  = "N2_"
     76     gnom(igas_CH4) = "CH4"
     77     gnom(igas_H2)  = "H2_"
     78     
     79     
    9680  else
    9781     write(*,*) 'Cannot find required file "gases.def"'
    9882     call abort
    9983  endif
     84 
    10085  close(90)
    10186!$OMP END MASTER

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

    r1647 r1648  
    2424      use radcommon_h, only : tgasref,tgasmin,tgasmax
    2525      use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
    26       use radcommon_h, only : wrefvar,WNOI,WNOV
     26      use radcommon_h, only : WNOI,WNOV
    2727      use datafile_mod, only: datadir
    2828      use comcstfi_mod, only: mugaz
    … …  
    4444      ! ALLOCATABLE ARRAYS -- AS 12/2011
    4545      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8    !read by master
    46       character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for check with gnom, read by master
    4746
    4847      real*8 x, xi(4), yi(4), ans, p
    … …  
    7978      read(111,*) ngas
    8079
    81       if(ngas.ne.ngasmx)then
    82          print*,'Number of gases in radiative transfer data (',ngas,') does not', &
    83                 'match that in gases.def (',ngasmx,'), exiting.'
    84          call abort
    85       endif
    86 
    8780      if(ngas.gt.5 .or. ngas.lt.1)then
    8881         print*,ngas,' species in database [',               &
    … …  
    9083                '], radiative code cannot handle this.'
    9184         call abort
    92       endif
    93 
    94       ! dynamically allocate gastype and read from Q.dat
    95       IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
    96 
    97       do igas=1,ngas
    98          read(111,*) gastype(igas)
    99          print*,'Gas ',igas,' is ',gastype(igas)
    100       enddo
    101 
    102       ! get array size, load the coefficients
    103       open(111,file=TRIM(file_path),form='formatted')
    104       read(111,*) L_REFVAR
    105       IF( .NOT. ALLOCATED( wrefvar ) ) ALLOCATE( WREFVAR(L_REFVAR) )
    106       read(111,*) wrefvar
    107       close(111)
    108 
    109       ! Check that gastype and gnom match
    110       do igas=1,ngas
    111          print*,'Gas ',igas,' is ',trim(gnom(igas))
    112          if (trim(gnom(igas)).ne.trim(gastype(igas))) then
    113             print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
    114                  'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', &
    115                  'gases.def with Q.dat in your radiative transfer directory.'
    116             call abort
    117          endif
    118       enddo
    119       print*,'Confirmed gas match in radiative transfer and gases.def!'
    120 
    121       ! display the values
    122       print*,'Variable gas volume mixing ratios:'
    123       do n=1,L_REFVAR
    124          !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
    125          print*,n,'.',wrefvar(n),' mol/mol'
    126       end do
    127       print*,''
     85      endif
     86     
     87      L_REFVAR = 1 ! JVO 2017 : set to 1 to keep the code running until the new variable species treatment
    12888
    12989!=======================================================================
    … …  
    628588            enddo
    629589
    630          endif
    631 
     590         elseif (igas .eq. igas_CH4) then
     591
     592            ! first do self-induced absorption
     593            dummy = -9999
     594            call interpolateCH4CH4(200.D+0,200.D+0,7500.D+0,testcont,.true.,dummy)
     595            ! then cross-interactions with other gases
     596            do jgas=1,ngasmx
     597               if (jgas .eq. igas_N2) then
     598                  dummy = -9999
     599                  call interpolateN2CH4(200.D+0,250.0D+0,100000.D+0,5000.D+0,testcont,.true.,dummy)
     600               endif
     601            enddo
     602
     603         endif 
     604         
    632605      enddo
    633606      endif
    … …  
    645618      IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
    646619      IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
    647       IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
    648620!$OMP END MASTER
    649621!$OMP BARRIER

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

    r1616 r1648  
    1       subroutine tpindex(pw,tw,qvar,pref,tref,wrefvar,LCOEF,MT,MP,
    2      &     NVAR,wratio)
     1      subroutine tpindex(pw,tw,pref,tref,LCOEF,MT,MP)
    32
    43!==================================================================
    … …  
    65!     Purpose
    76!     -------
    8 !     Interpolate K-coefficients to the given P,T and Qvar values.
     7!     Interpolate K-coefficients to the given P,T.
    98!
    109!     Notes
    … …  
    5857      real*8 Tref(L_NTREF)
    5958      real*8 pref(L_PINT)
    60       real*8 wrefvar(L_REFVAR)
    6159
    62       integer MT, MP, N, M, NP, NVAR
    63       real*8  PW, TW, Qvar, wratio
     60      integer MT, MP, N, M, NP
     61      real*8  PW, TW
    6462      real*8  PWL, LCOEF(4), T, U
    6563
    … …  
    143141      LCOEF(4) = (1.0-T)*U
    144142
    145 !  Get the indicies for abundance of the varying species. There are 10 sets of
    146 !  k-coefficients with differing amounts of variable vs. constant gas.
    147 
    148       IF(QVAR.le.WREFVAR(1)) then
    149         NVAR   = 1
    150         WRATIO = 0.0D0      ! put all the weight on the first point
    151       ELSEIF(QVAR.ge.WREFVAR(L_REFVAR)) then
    152         NVAR   = L_REFVAR-1 ! TB16 in order to not oversize NVAr when doing
    153                                   !NVAR+1
    154         WRATIO = 1.00D0     ! put all the weight on the last point
    155       ELSE
    156         DO N=2,L_REFVAR
    157           IF(QVAR.GE.WREFVAR(N-1) .and. QVAR.lt.WREFVAR(N)) then
    158             NVAR   = N-1
    159             WRATIO = (QVAR - WREFVAR(N-1))/(WREFVAR(N) - WREFVAR(N-1))
    160             GOTO 30
    161           END IF
    162         END DO
    163       END IF
    164 
    165143   30 CONTINUE
    166144

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

    r1647 r1648  
    1919      real varian      ! Characteristic variance of log-normal distribution
    2020      real r3n_q     ! used to compute r0 from number and mass mixing ratio
    21       real rho_dust     ! Mars dust density (kg.m-3)
    2221      real ref_r0        ! for computing reff=ref_r0*r0 (in log.n. distribution)
    2322!$OMP THREADPRIVATE(noms,mmol,radius,rho_q,qext,alpha_lift,alpha_devil,qextrhor, &
    24         !$OMP varian,r3n_q,rho_dust,rho_ice,ref_r0)
     23        !$OMP varian,r3n_q,ref_r0)
    2524
    2625! tracer indexes: these are initialized in initracer and should be 0 if the
Note: See TracChangeset for help on using the changeset viewer.