Ignore:
Timestamp:
Nov 30, 2018, 12:55:49 AM (6 years ago)
Author:
jvatant
Message:

Major radiative transfer contribution : Add the 'corrk_recombin' option that allows to use
correlated-k for single species instead of pre-mix and enables more flexiblity for variable species.
-> Algorithm inspired from Lacis and Oinas 1991 and Amundsen et al 2016

+ Added 'recombin_corrk.F90' - Important improvements in 'sugas_corrk.90' and 'callcork.F90'

  • Must have the desired variable species as tracers -> TBD : Enable to force composition even if no tracers
  • To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values who match the atmospheric conditions ) which is then processed as a standard pre-mix in optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.


READ CAREFULY :

  • In case of 'corrk_recombin', the variable L_NREFVAR doesn't have the same meaning as before and doesn't stand for the different mixing ratios but the different species.
  • Input corr-k should be found in corrkdir within 'corrk_gcm_IR/VI_XXX.dat' and can contain a 'fixed' specie ( compulsory if you include self-broadening ) that MUST have been created with correct mixing ratios, or a variable specie for which mixing ratio MUST have been set to 1 ( no self-broadening then, assume it's a trace sepecie ) -> You can't neither have CIA of variable species included upstream in the corr-k
Location:
trunk/LMDZ.TITAN/libf/phytitan
Files:
1 added
9 edited

Legend:

Unmodified
Added
Removed

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

    r2046 r2050  
    1515      use callkeys_mod, only: global1d, szangle
    1616      use comcstfi_mod, only: pi, mugaz, cpp
    17       use callkeys_mod, only: diurnal,tracer,seashaze,         &
     17      use callkeys_mod, only: diurnal,tracer,seashaze,corrk_recombin,   &
    1818                              strictboundcorrk,specOLR
    1919      use geometry_mod, only: latitude
    … …  
    3333!     Emmanuel 01/2001, Forget 09/2001
    3434!     Robin Wordsworth (2009)
     35!     Jan Vatant d'Ollone (2018) -> corrk recombining case
    3536!
    3637!==================================================================
    … …  
    4647      INTEGER,INTENT(IN) :: ngrid                  ! Number of atmospheric columns.
    4748      INTEGER,INTENT(IN) :: nlayer                 ! Number of atmospheric layers.
    48       REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (X/m2).
     49      REAL,INTENT(IN) :: pq(ngrid,nlayer,nq)       ! Tracers (X/kg).
    4950      INTEGER,INTENT(IN) :: nq                     ! Number of tracers.
    5051      REAL,INTENT(IN) :: qsurf(ngrid,nq)           ! Tracers on surface (kg.m-2).
    … …  
    105106      REAL*8 albv(L_NSPECTV)                         ! Spectral Visible Albedo.
    106107
    107       INTEGER ig,l,k,nw
     108      INTEGER ig,l,k,nw,iq,ip,ilay,it
     109     
     110      LOGICAL found
    108111
    109112      real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
    … …  
    120123      REAL*8 seashazefact(L_LEVELS)
    121124
     125      ! For muphys optics
     126      REAL*8 pqmo(ngrid,nlayer,nmicro)  ! Tracers for microphysics optics (X/m2).
     127      REAL*8 i2e(ngrid,nlayer)          ! int 2 ext factor ( X.kg-1 -> X.m-2 for optics )
     128     
     129      ! For corr-k recombining
     130      REAL*8 pqr(ngrid,L_PINT,L_REFVAR)  ! Tracers for corr-k recombining (mol/mol).
     131      REAL*8 fact, tmin, tmax
     132     
     133      LOGICAL usept(L_PINT,L_NTREF)  ! mask if pfref grid point will be used
     134      INTEGER inflay(L_PINT)         ! nearest inferior GCM layer for pfgasref grid points
     135     
    122136!=======================================================================
    123137!             I.  Initialization on every call   
    … …  
    130144      end do
    131145
     146      ! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2
     147      ! NOTE: it should be moved somewhere else: calmufi performs the same kind of
     148      ! computations... waste of time...
     149      i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
     150      pqmo(:,:,:) = 0.0
     151      DO iq=1,nmicro
     152        pqmo(:,:,iq) = pq(:,:,iq)*i2e(:,:)
     153      ENDDO
     154
     155      ! Default value for fixed species for whom vmr has been taken
     156      ! into account while computing high-resolution spectra
     157      if (corrk_recombin) pqr(:,:,:) = 1.0
    132158
    133159!-----------------------------------------------------------------------   
    134160      do ig=1,ngrid ! Starting Big Loop over every GCM column
    135161!-----------------------------------------------------------------------
    136 
     162         
     163         ! Recombine reference corr-k if needed
     164         if (corrk_recombin) then
     165         
     166         ! NB : To have decent CPU time recombining is not done on all gridpoints and wavelenghts but we
     167         ! calculate a gasi/v_recomb variable on the reference corrk-k T,P grid (only for T,P values
     168         ! who match the atmospheric conditions ) which is then processed as a standard pre-mix in
     169         ! optci/v routines, but updated every time tracers on the ref P grid have varied > 1%.
     170
     171            ! Extract tracers for variable radiative species
     172            ! Also find the nearest GCM layer under each ref pressure
     173            do ip=1,L_PINT
     174               
     175               ilay=0
     176               found = .false.
     177               do l=1,nlayer
     178                  if ( pplay(ig,l) .gt. 10.0**(pfgasref(ip)+2.0) ) then ! pfgasref=log(p[mbar])
     179                     found=.true.
     180                     ilay=l
     181                  endif
     182               enddo       
     183
     184               if (.not. found ) then ! set to min
     185                  do iq=1,L_REFVAR
     186                     if ( radvar_mask(iq) ) then
     187                        pqr(ig,ip,iq) = pq(ig,1,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
     188                     endif
     189                  enddo
     190               else
     191                  if (ilay==nlayer) then ! set to max
     192                     do iq=1,L_REFVAR
     193                        if ( radvar_mask(iq) ) then
     194                           pqr(ig,ip,iq) = pq(ig,nlayer,radvar_indx(iq)) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
     195                        endif
     196                     enddo
     197                  else ! standard
     198                     fact = ( 10.0**(pfgasref(ip)+2.0) - pplay(1,ilay+1) ) / ( pplay(1,ilay) - pplay(1,ilay+1) ) ! pfgasref=log(p[mbar])
     199                     do iq=1,L_REFVAR
     200                        if ( radvar_mask(iq) ) then
     201                           pqr(ig,ip,iq) = pq(ig,ilay,radvar_indx(iq))**fact * pq(ig,ilay+1,radvar_indx(iq))**(1.0-fact)
     202                           pqr(ig,ip,iq) = pqr(ig,ip,iq) / rat_mmol(radvar_indx(iq)-nmicro) ! mol/mol
     203                        endif
     204                     enddo
     205                  endif ! if ilay==nlayer
     206               endif ! if not found
     207
     208               inflay(ip) = ilay
     209
     210            enddo ! ip=1,L_PINT
     211
     212            ! NB : The following usept is a trick to call recombine only for the reference T-P
     213            ! grid points that are useful given the temperature range at this altitude
     214            ! It saves a looot of time - JVO 18
     215            usept(:,:) = .true.
     216            do ip=1,L_PINT-1
     217              if ( inflay(ip+1)==nlayer ) then
     218                usept(ip,:) = .false.
     219              endif
     220              if ( inflay(ip) == 0 ) then
     221                usept(ip+1:,:) = .false.
     222              endif
     223              if ( usept(ip,1) ) then ! if not all false
     224                tmin = minval(pt(ig,min(inflay(ip+1)+1,nlayer):max(inflay(ip),1)))
     225                tmax = maxval(pt(ig,min(inflay(ip+1)+1,nlayer):max(inflay(ip),1)))
     226                do it=1,L_NTREF-1
     227                  if ( tgasref(it+1) .lt. tmin ) then
     228                    usept(ip,it) = .false.
     229                  endif
     230                enddo
     231                do it=2,L_NTREF
     232                  if ( tgasref(it-1) .gt. tmax ) then
     233                    usept(ip,it) = .false.
     234                  endif
     235                enddo
     236                ! in case of out-of-bounds
     237                if ( tgasref(1)         .lt. tmin ) usept(ip,1) = .true.
     238                if ( tgasref(L_NTREF)   .gt. tmax ) usept(ip,L_NTREF) = .true.
     239              endif
     240            enddo ! ip=1,L_PINT-1
     241            ! deal with last bound
     242            if ( inflay(L_PINT-1).ne.0 ) usept(L_PINT,:) = usept(L_PINT-1,:)
     243
     244
     245            do ip=1,L_PINT
     246
     247               ! Recombine k at (useful only!) reference T-P values if tracers or T have enough varied
     248               do it=1,L_NTREF
     249
     250                 if ( usept(ip,it) .eqv. .false. ) cycle
     251
     252                 do l=1,L_REFVAR
     253                    if ( abs( (pqr(ig,ip,l) - pqrold(ip,l)) / max(1.0e-30,pqrold(ip,l))) .GT. 0.01  & ! +- 1%
     254                         .or. ( useptold(ip,it) .eqv. .false.  ) ) then ! in case T change but not the tracers
     255                       call recombin_corrk( pqr(ig,ip,:),ip,it )
     256                       exit ! one is enough to trigger the update
     257                    endif
     258                 enddo
     259                 
     260               enddo
     261
     262            enddo ! ip=1,L_PINT
     263
     264            useptold(:,:)=usept(:,:)
     265
     266       endif ! if corrk_recombin
    137267
    138268!=======================================================================
    … …  
    304434            endif
    305435           
    306             call optcv(pq(ig,:,1:nmicro),nlayer,plevrad,tmid,pmid,        &
     436            call optcv(pqmo(ig,:,:),nlayer,plevrad,tmid,pmid,   &
    307437                 dtauv,tauv,taucumv,wbarv,cosbv,tauray,taugsurf,seashazefact)
    308438
    … …  
    346476!-----------------------------------------------------------------------
    347477
    348          call optci(pq(ig,:,1:nmicro),nlayer,plevrad,tlevrad,tmid,pmid,    &
     478         call optci(pqmo(ig,:,:),nlayer,plevrad,tlevrad,tmid,pmid,   &
    349479              dtaui,taucumi,cosbi,wbari,taugsurfi,seashazefact)
    350480
    … …  
    459589      endif
    460590
    461       ! See physiq.F for explanations about CLFvarying. This is temporary.
    462591      if (lastcall) then
    463592        IF( ALLOCATED( gasi ) ) DEALLOCATE( gasi )
    464593        IF( ALLOCATED( gasv ) ) DEALLOCATE( gasv )
     594        IF( ALLOCATED( gasi_recomb ) ) DEALLOCATE( gasi_recomb )
     595        IF( ALLOCATED( gasv_recomb ) ) DEALLOCATE( gasv_recomb )
     596        IF( ALLOCATED( pqrold ) ) DEALLOCATE( pqrold )
     597        IF( ALLOCATED( useptold ) ) DEALLOCATE( useptold )
    465598!$OMP BARRIER
    466599!$OMP MASTER

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

    r2046 r2050  
    1414      logical,save :: strictboundcorrk
    1515!$OMP THREADPRIVATE(strictboundcorrk)
     16      logical,save :: corrk_recombin
     17!$OMP_THREADPRIVATE(corrk_recombin)
    1618      logical,save :: seashaze,uncoupl_optic_haze
    1719!$OMP THREADPRIVATE(seashaze,uncoupl_optic_haze)

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

    r2046 r2050  
    234234       call getin_p("corrkdir",corrkdir)
    235235       write(*,*) " corrkdir = ",corrkdir
     236       
     237       write(*,*) "use correlated-k recombination instead of pre-mixed values ?"
     238       corrk_recombin=.true. ! default value
     239       call getin_p("corrk_recombin",corrk_recombin)
     240       write(*,*) " corrk_recombin = ",corrk_recombin
    236241     endif
    237242     

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

    r2046 r2050  
    1 subroutine optci(PQO,NLAY,PLEV,TLEV,TMID,PMID,      &
     1subroutine optci(PQMO,NLAY,PLEV,TLEV,TMID,PMID,      &
    22     DTAUI,TAUCUMI,COSBI,WBARI,TAUGSURF,SEASHAZEFACT)
    33
    44  use radinc_h
    5   use radcommon_h, only: gasi,tlimit,Cmk,gzlat_ig,tgasref,pfgasref,wnoi,scalep,indi,gweight
     5  use radcommon_h, only: gasi,gasi_recomb,tlimit,Cmk,gzlat_ig, &
     6                         tgasref,pfgasref,wnoi,scalep,indi,gweight
    67  use gases_h
    78  use comcstfi_mod, only: r
    8   use callkeys_mod, only: continuum,graybody,callclouds,callmufi,seashaze,uncoupl_optic_haze
     9  use callkeys_mod, only: continuum,graybody,corrk_recombin,      &
     10                          callclouds,callmufi,seashaze,uncoupl_optic_haze
    911  use tracer_h, only : nmicro,nice
    1012  use MMP_OPTICS
    … …  
    3739  ! Input/Output
    3840  !==========================================================
    39   REAL*8, INTENT(IN)  :: PQO(nlay,nmicro) ! Tracers (X/m2).
    40   INTEGER, INTENT(IN) :: NLAY             ! Number of pressure layers (for pqo)
     41  REAL*8, INTENT(IN)  :: PQMO(nlay,nmicro)  ! Tracers for microphysics optics (X/m2).
     42  INTEGER, INTENT(IN) :: NLAY               ! Number of pressure layers (for pqmo)
    4143  REAL*8, INTENT(IN)  :: PLEV(L_LEVELS), TLEV(L_LEVELS)
    4244  REAL*8, INTENT(IN)  :: TMID(L_LEVELS), PMID(L_LEVELS)
    … …  
    144146        ! as long as the microphysics only isn't fully debugged -- JVO 01/18
    145147        IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    146           m0as = pqo(ilay,1)
    147           m3as = pqo(ilay,2)
    148           m0af = pqo(ilay,3)
    149           m3af = pqo(ilay,4)
     148          m0as = pqmo(ilay,1)
     149          m3as = pqmo(ilay,2)
     150          m0af = pqmo(ilay,3)
     151          m3af = pqmo(ilay,4)
    150152
    151153          IF (.NOT.mmp_sph_optics_ir(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
    … …  
    238240           ! transfer on the tested simulations !
    239241
    240            tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     242           if (corrk_recombin) then
     243             tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
     244           else
     245             tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     246           endif
    241247
    242248           KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)

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

    r2046 r2050  
    1 SUBROUTINE OPTCV(PQO,NLAY,PLEV,TMID,PMID,  &
     1SUBROUTINE OPTCV(PQMO,NLAY,PLEV,TMID,PMID,  &
    22     DTAUV,TAUV,TAUCUMV,WBARV,COSBV,TAURAY,TAUGSURF,SEASHAZEFACT)
    33
    44  use radinc_h
    5   use radcommon_h, only: gasv,tlimit,Cmk,gzlat_ig,tgasref,pfgasref,wnov,scalep,indv,gweight
     5  use radcommon_h, only: gasv,gasv_recomb,tlimit,Cmk,gzlat_ig, &
     6                         tgasref,pfgasref,wnov,scalep,indv,gweight
    67  use gases_h
    78  use comcstfi_mod, only: r
    8   use callkeys_mod, only: continuum,graybody,callgasvis,callclouds,callmufi,seashaze,uncoupl_optic_haze
     9  use callkeys_mod, only: continuum,graybody,callgasvis,corrk_recombin,  &
     10                          callclouds,callmufi,seashaze,uncoupl_optic_haze
    911  use tracer_h, only: nmicro,nice
    1012  use MMP_OPTICS
    … …  
    4345  ! Input/Output
    4446  !==========================================================
    45   REAL*8, INTENT(IN)  :: PQO(nlay,nmicro) ! Tracers (X/m2).
    46   INTEGER, INTENT(IN) :: NLAY             ! Number of pressure layers (for pqo)
     47  REAL*8, INTENT(IN)  :: PQMO(nlay,nmicro)  ! Tracers for microphysics optics (X/m2).
     48  INTEGER, INTENT(IN) :: NLAY               ! Number of pressure layers (for pqmo)
    4749  REAL*8, INTENT(IN)  :: PLEV(L_LEVELS)
    4850  REAL*8, INTENT(IN)  :: TMID(L_LEVELS), PMID(L_LEVELS)
    … …  
    167169        ! as long as the microphysics only isn't fully debugged -- JVO 01/18
    168170        IF (callmufi .AND. (.NOT. uncoupl_optic_haze)) THEN
    169           m0as = pqo(ilay,1)
    170           m3as = pqo(ilay,2)
    171           m0af = pqo(ilay,3)
    172           m3af = pqo(ilay,4)
     171          m0as = pqmo(ilay,1)
     172          m3as = pqmo(ilay,2)
     173          m0af = pqmo(ilay,3)
     174          m3af = pqmo(ilay,4)
    173175
    174176          IF (.NOT.mmp_sph_optics_vis(m0as,m3as,nw,ext_s,sca_s,ssa_s,asf_s)) &
    … …  
    267269           ! transfer on the tested simulations !
    268270
    269            tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     271           if (corrk_recombin) then
     272             tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NW,NG)
     273           else
     274             tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     275           endif
    270276             
    271277           KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
    … …  
    302308  !     Now the full treatment for the layers, where besides the opacity
    303309  !     we need to calculate the scattering albedo and asymmetry factors
     310  ! ======================================================================
    304311
    305312  ! Haze scattering

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

    r2046 r2050  
    391391!$OMP THREADPRIVATE(tankCH4)
    392392
    393 
    394       ! -----******----- FOR MUPHYS OPTICS -----******-----
    395       integer :: i,j
    396       real :: pqo(ngrid,nlayer,nq)   ! Tracers for the optics (X/m2).
    397       ! -----******----- END FOR MUPHYS OPTICS -----******-----
    398 
    399       real :: i2e(ngrid,nlayer)      ! int 2 ext factor ( X.kg-1 -> X.m-3 or X.m-2 , depending if optics or diags )
     393      real :: i2e(ngrid,nlayer)      ! int 2 ext factor ( X.kg-1 -> X.m-3 for diags )
    400394
    401395      real,save,dimension(:,:,:), allocatable :: tpq ! Tracers for decoupled microphysical tests ( temporary in 01/18 )
    … …  
    844838               call call_profilgases(nlayer)
    845839
    846                ! Convert (microphysical) tracers for optics: X.kg-1 --> X.m-2_
    847                ! NOTE: it should be moved somewhere else: calmufi performs the same kind of
    848                ! computations... waste of time...
    849                i2e(:,1:nlayer) = ( pplev(:,1:nlayer)-pplev(:,2:nlayer+1) ) / gzlat(:,1:nlayer)
    850                pqo(:,:,:) = 0.0
    851                DO j=1,nmicro-nice
    852                  pqo(:,:,j) = pq(:,:,j)*i2e(:,:)
    853                ENDDO
    854 
    855840               ! standard callcorrk
    856                call callcorrk(ngrid,nlayer,pqo,nq,qsurf,zday,                     &
     841               call callcorrk(ngrid,nlayer,pq,nq,qsurf,zday,                     &
    857842                              albedo,albedo_equivalent,emis,mu0,pplev,pplay,pt,   &
    858843                              tsurf,fract,dist_star,                              &

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

    r2026 r2050  
    5656      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv
    5757      REAL*8, DIMENSION(:), ALLOCATABLE :: PGASREF, TGASREF, PFGASREF, GWEIGHT
     58 
     59      ! For corrk recombining - JVO 18
     60      REAL*8,  DIMENSION(:,:,:,:), ALLOCATABLE :: gasi_recomb, gasv_recomb
     61      REAL*8,  DIMENSION(:,:),     ALLOCATABLE :: pqrold
     62      REAL*8,  DIMENSION(:),       ALLOCATABLE :: w_cum
     63      LOGICAL, DIMENSION(:,:),     ALLOCATABLE :: useptold
     64      INTEGER, DIMENSION(:),       ALLOCATABLE :: permut_idx
     65!$OMP THREADPRIVATE(gasi_recomb,gasv_recomb,pqrold,w_cum,useptold,permut_idx)
     66     
     67      LOGICAL, DIMENSION(:), ALLOCATABLE :: RADVAR_MASK ! for corrk recombin -> read by master in sugas_corrk
     68      INTEGER, DIMENSION(:), ALLOCATABLE :: RADVAR_INDX ! for corrk recombin -> read by master in sugas_corrk
     69     
    5870      real*8 FZEROI(L_NSPECTI)
    5971      real*8 FZEROV(L_NSPECTV)
    6072      real*8 pgasmin, pgasmax
    6173      real*8 tgasmin, tgasmax
    62 !$OMP THREADPRIVATE(gasi,gasv,&  !pgasref,tgasref,pfgasref read by master in sugas_corrk
     74!$OMP THREADPRIVATE(gasi,gasv,&  !pgasref,tgasref,pfgasref, gweight read by master in sugas_corrk
    6375        !$OMP FZEROI,FZEROV)     !pgasmin,pgasmax,tgasmin,tgasmax read by master in sugas_corrk
    6476
    … …  
    7385      real*8,parameter :: UBARI = 0.5D0
    7486
    75 !$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFvis,omegaREFir,&     ! gweight read by master in sugas_corrk
     87!$OMP THREADPRIVATE(QREFvis,QREFir,omegaREFvis,omegaREFir,&
    7688                !$OMP tstellar,planckir,PTOP)
    7789

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

    r2026 r2050  
    4040!     L_PINT    is the number of Lagrange interpolated reference
    4141!               pressures for the gas k-coefficients - now for a
    42 !               smaller p-grid than before
     42!               finer p-grid than before
    4343!     L_NTREF   is the number of reference temperatures for the
    4444!               k-coefficients
    … …  
    4646!               to this value
    4747!
    48 !     L_REFVAR  The number of different mixing ratio values for
    49 !               the k-coefficients. Variable component of the mixture
    50 !               can in princple be anything: currently it's H2O.
     48!     L_REFVAR 
     49!               EITHER (corrk_recombin = false)
     50!      _          The number of different mixing ratio values for
     51!     / \         the k-coefficients. Variable component of the mixture
     52!    /   \        can in princple be anything: currently it's H2O.
     53!   /  !  \
     54!  /_______\    OR (corrk_recombin = true)
     55!                 The TOTAL number of species to recombine when we
     56!                 recombine corr-k instead of pre-mixed values.
    5157!----------------------------------------------------------------------
    5258

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

    r2026 r2050  
    2323      use radcommon_h, only : pgasref,pfgasref,pgasmin,pgasmax
    2424      use radcommon_h, only : tgasref,tgasmin,tgasmax
    25       use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight
     25      use radcommon_h, only : gasv,gasi,FZEROI,FZEROV,gweight, w_cum
    2626      use radcommon_h, only : WNOI,WNOV
     27      use radcommon_h, only : gasi_recomb, gasv_recomb, pqrold, useptold
     28      use radcommon_h, only : radvar_mask, radvar_indx, permut_idx
    2729      use datafile_mod, only: datadir, corrkdir, banddir
    2830      use comcstfi_mod, only: mugaz
    2931      use gases_h
    3032      use ioipsl_getin_p_mod, only: getin_p
    31       use callkeys_mod, only: graybody,callgasvis, continuum
     33      use callkeys_mod, only: graybody,callgasvis, continuum, tracer, corrk_recombin
     34      use tracer_h, only: noms, nqtot_p
    3235      implicit none
    3336
    … …  
    4346      ! ALLOCATABLE ARRAYS -- AS 12/2011
    4447      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE,SAVE :: gasi8, gasv8    !read by master
     48      character*20,allocatable,DIMENSION(:),SAVE :: gastype ! for name check, read by master
    4549
    4650      real*8 x, xi(4), yi(4), ans, p
    … …  
    5458
    5559      integer :: dummy
     60     
     61      integer ngasvar
     62      logical found
    5663
    5764!=======================================================================
    … …  
    7784      read(111,*) ngas
    7885
    79       if(ngas.gt.5 .or. ngas.lt.1)then
     86      if(ngas.lt.1)then
    8087         print*,ngas,' species in database [',               &
    8188                     corrkdir(1:LEN_TRIM(corrkdir)),           &
    … …  
    8491      endif
    8592     
    86       L_REFVAR = 1 ! JVO 2017 : set to 1 to keep the code running until the new variable species treatment
     93      if ( corrk_recombin ) then
     94         L_REFVAR =  ngas
     95         
     96         ! dynamical allocations and read from Q.dat
     97         IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( L_REFVAR ) )
     98         IF ( .NOT. ALLOCATED( radvar_mask ) ) ALLOCATE( radvar_mask( L_REFVAR ) )         
     99         IF ( .NOT. ALLOCATED( radvar_indx ) ) ALLOCATE( radvar_indx( L_REFVAR ) )
     100         
     101         ngasvar = 0
     102         
     103         write(*,*) ' '
     104         do igas=1,L_REFVAR
     105           read(111,*) gastype(igas)
     106           read(111,*) radvar_mask(igas) ! boolean
     107           print *, 'Found radiatively active specie : ',trim(gastype(igas)), ' : Variable =', radvar_mask(igas)
     108           if ( radvar_mask(igas) ) ngasvar = ngasvar+1
     109         enddo
     110         
     111         print *, ' '
     112         print *, 'TOTAL number of species found for radiative transfer = ', L_REFVAR
     113         print *, '      including : ', ngasvar, 'variable species.'
     114         print *, ' '
     115       
     116         if (.not. tracer .and. ngasvar .gt. 0) then
     117            write(*,*) 'sugas_corrk:error: You try to run variable species for corrk recombining'
     118            write(*,*) 'in the radiative transfer but without tracers : I will stop here ! '
     119            call abort
     120         endif
     121         
     122         ! NB (JVO 18) We should enable here to force composition even if no associated tracers
     123     
     124         ! Check that we have matching tracers for the variable radiative active species
     125         if ( ngasvar .GT. 0 ) then
     126
     127           do igas=1,L_REFVAR
     128            if ( radvar_mask(igas) ) then
     129           
     130               found = .false.
     131               do jgas=1,nqtot_p
     132                 if ( trim(noms(jgas)) == trim(gastype(igas)) .OR. &
     133                      trim(noms(jgas)) == trim(gastype(igas))//'_') then
     134                    radvar_indx(igas) = jgas
     135                    found=.true.
     136                    exit
     137                 endif
     138               enddo
     139               
     140               if (.not. found) then
     141                 write(*,*) "sugas_corrk:error: "//trim(gastype(igas))//" is missing from tracers."
     142                 call abort
     143               endif
     144               
     145             endif
     146           enddo ! igas=1,L_REFVAR
     147           
     148         endif ! if ngasvar .gt. 0
     149     
     150      else
     151        L_REFVAR = 1 ! JVO 2017 : those will not be used, just set to 1 to keep the code running
     152      endif ! if corrk_recombin
    87153
    88154!=======================================================================
    … …  
    118184      print*,''
    119185
     186      IF (.NOT. ALLOCATED(permut_idx)) ALLOCATE(permut_idx(L_NGAUSS*L_NGAUSS)) ! for corr-k recombining
     187      permut_idx = (/(i, i=1,L_NGAUSS*L_NGAUSS)/) ! for the recombin_corrk firstcall
     188
     189      IF (.NOT. ALLOCATED(w_cum)) ALLOCATE(w_cum(L_NGAUSS)) ! for corr-k recombining     
     190      w_cum(1)= gweight(1)
     191      DO n=2,L_NGAUSS
     192        w_cum(n) = w_cum(n-1)+gweight(n)
     193      ENDDO
     194
    120195!=======================================================================
    121196!     Set the reference pressure and temperature arrays.  These are
    … …  
    213288      IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
    214289
     290      if (corrk_recombin) then
     291        IF( .NOT. ALLOCATED( gasi_recomb ) ) ALLOCATE( gasi_recomb(L_NTREF,L_PINT,L_NSPECTI,L_NGAUSS) )
     292        IF( .NOT. ALLOCATED( gasv_recomb ) ) ALLOCATE( gasv_recomb(L_NTREF,L_PINT,L_NSPECTV,L_NGAUSS) )
     293        IF( .NOT. ALLOCATED( pqrold ) ) ALLOCATE( pqrold(L_PINT,L_REFVAR) )
     294        IF( .NOT. ALLOCATED( useptold ) ) ALLOCATE( useptold(L_PINT,L_NTREF) )
     295        pqrold(:,:) = 0.0
     296        useptold(:,:) = .false.
     297      endif
     298     
    215299      ! display the values
    216300      print*,''
    … …  
    272356            print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
    273357         else
    274             file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
    275             file_path=TRIM(datadir)//TRIM(file_id)
     358             
     359            if (corrk_recombin .eqv. .false. ) then
     360              file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
     361              file_path=TRIM(datadir)//TRIM(file_id)
     362             
     363              ! check that the file exists
     364              inquire(FILE=file_path,EXIST=file_ok)
     365              if(.not.file_ok) then
     366                 write(*,*)'The file ',TRIM(file_path)
     367                 write(*,*)'was not found by sugas_corrk.F90.'
     368                 write(*,*)'Are you sure you have absorption data for these bands?'
     369                 call abort
     370              endif
     371           
     372              open(111,file=TRIM(file_path),form='formatted')
     373              read(111,*) gasv8
     374              close(111)
     375             
     376            else
     377             
     378              do igas=1,L_REFVAR
     379                ! JVO 18 : To comply with Q.dat watch out for 2-letters gas trailing underscore :
     380                ! e.g. if ever you have N2_ (symmetric so shouldn't happen) file should be ..._N2_.dat
     381                file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI_'//trim(adjustl(gastype(igas)))//'.dat'
     382                file_path=TRIM(datadir)//TRIM(file_id)
     383               
     384                ! check that the file exists
     385                inquire(FILE=file_path,EXIST=file_ok)
     386                if(.not.file_ok) then
     387                   write(*,*)'The file ',TRIM(file_path)
     388                   write(*,*)'was not found by sugas_corrk.F90.'
     389                   write(*,*)'Are you sure you have absorption data for these bands and this specie?'
     390                   call abort
     391                 endif
     392           
     393                open(111,file=TRIM(file_path),form='formatted')
     394                read(111,*) gasv8(:,:,igas,:,:)
     395                close(111)
     396              enddo             
     397             
     398            endif !  if corrk_recombin
    276399           
    277             ! check that the file exists
    278             inquire(FILE=file_path,EXIST=file_ok)
    279             if(.not.file_ok) then
    280                write(*,*)'The file ',TRIM(file_path)
    281                write(*,*)'was not found by sugas_corrk.F90.'
    282                write(*,*)'Are you sure you have absorption data for these bands?'
    283                call abort
    284             endif
    285          
    286             open(111,file=TRIM(file_path),form='formatted')
    287             read(111,*) gasv8
    288             close(111)
    289400         end if
    290401
    … …  
    316427!$OMP BARRIER
    317428      else
    318          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
    319          file_path=TRIM(datadir)//TRIM(file_id)
    320429     
    321          ! check that the file exists
    322          inquire(FILE=file_path,EXIST=file_ok)
    323          if(.not.file_ok) then
    324             write(*,*)'The file ',TRIM(file_path)
    325             write(*,*)'was not found by sugas_corrk.F90.'
    326             write(*,*)'Are you sure you have absorption data for these bands?'
    327             call abort
    328          endif
    329          
     430         if ( corrk_recombin .eqv. .false. ) then
     431           file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
     432           file_path=TRIM(datadir)//TRIM(file_id)
     433       
     434           ! check that the file exists
     435           inquire(FILE=file_path,EXIST=file_ok)
     436           if(.not.file_ok) then
     437              write(*,*)'The file ',TRIM(file_path)
     438              write(*,*)'was not found by sugas_corrk.F90.'
     439              write(*,*)'Are you sure you have absorption data for these bands?'
     440              call abort
     441           endif
     442           
    330443!$OMP MASTER                   
    331          open(111,file=TRIM(file_path),form='formatted')
    332          read(111,*) gasi8
    333          close(111)
     444           open(111,file=TRIM(file_path),form='formatted')
     445           read(111,*) gasi8
     446           close(111)
    334447!$OMP END MASTER
    335448!$OMP BARRIER
     449
     450         else
     451         
     452           do igas=1,L_REFVAR
     453             ! JVO 18 : To comply with Q.dat watch out for 2-letters gas trailing underscore :
     454             ! e.g. if ever tou have N2_ (symmetric so shouldn't happen) file should be ..._N2_.dat
     455             file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR_'//trim(adjustl(gastype(igas)))//'.dat'
     456             file_path=TRIM(datadir)//TRIM(file_id)
     457         
     458             ! check that the file exists
     459             inquire(FILE=file_path,EXIST=file_ok)
     460             if(.not.file_ok) then
     461                write(*,*)'The file ',TRIM(file_path)
     462                write(*,*)'was not found by sugas_corrk.F90.'
     463                write(*,*)'Are you sure you have absorption data for these bands and this specie?'
     464                call abort
     465             endif
     466             
     467!$OMP MASTER                   
     468             open(111,file=TRIM(file_path),form='formatted')
     469             read(111,*) gasi8(:,:,igas,:,:)
     470             close(111)
     471!$OMP END MASTER
     472!$OMP BARRIER
     473           enddo ! do igas=1,L_REFVAR
     474           
     475         endif! if corrk_recombin
     476         
    336477     
    337478         ! 'fzero' is a currently unused feature that allows optimisation
    … …  
    559700      end do
    560701
     702      gasi_recomb(:,:,:,:) = kappa_IR ! non-zero init
     703      gasv_recomb(:,:,:,:) = kappa_VI ! non-zero init
    561704
    562705!=======================================================================
    … …  
    613756      IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
    614757      IF( ALLOCATED( pgasref ) ) DEALLOCATE( pgasref )
     758      IF( ALLOCATED( gastype ) ) DEALLOCATE( gastype )
    615759!$OMP END MASTER
    616760!$OMP BARRIER
Note: See TracChangeset for help on using the changeset viewer.