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)

File:
1 edited

Legend:

Unmodified
Added
Removed

  • 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
Note: See TracChangeset for help on using the changeset viewer.