Ignore:
Timestamp:
Jul 5, 2021, 4:04:28 PM (3 years ago)
Author:
aslmd
Message:

Generic GCM:

Adding k-coefficients mixing on the fly
Working with MordernTrac?

JVO + YJ

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
1 added
6 edited

Legend:

Unmodified
Added
Removed
  • trunk/LMDZ.GENERIC/libf/phystd/callcorrk.F90

    r2537 r2543  
    3636      use optcv_mod, only: optcv
    3737      use optci_mod, only: optci
     38      use recombin_corrk_mod, only: corrk_recombin, call_recombin
    3839      implicit none
    3940
     
    811812!=======================================================================
    812813
     814! ----------------------------------------------------------------
     815! Recombine reference corrk tables if needed - Added by JVO, 2020.
     816         if (corrk_recombin) then
     817           call call_recombin(ig,nlayer,pq(ig,:,:),pplay(ig,:),pt(ig,:),qvar(:),tmid(:),pmid(:))
     818         endif
     819! ----------------------------------------------------------------
    813820
    814821         Cmk= 0.01 * 1.0 / (glat(ig) * mugaz * 1.672621e-27) ! q_main=1.0 assumed.
  • trunk/LMDZ.GENERIC/libf/phystd/initracer.F

    r2542 r2543  
    44      USE tracer_h
    55      USE callkeys_mod, only: water
     6      USE recombin_corrk_mod, ONLY: ini_recombin
    67      IMPLICIT NONE
    78c=======================================================================
     
    8182       !! we allocate once for all arrays in common in tracer_h.F90
    8283       !! (supposedly those are not used before call to initracer)
    83        IF (.NOT.ALLOCATED(noms))         ALLOCATE(noms(nq))
    84        IF (.NOT.ALLOCATED(mmol))         ALLOCATE(mmol(nq))
    85        IF (.NOT.ALLOCATED(aki))          ALLOCATE(aki(nqtot))
    86        IF (.NOT.ALLOCATED(cpi))          ALLOCATE(cpi(nqtot))
    87        IF (.NOT.ALLOCATED(radius))       ALLOCATE(radius(nq))
    88        IF (.NOT.ALLOCATED(rho_q))        ALLOCATE(rho_q(nq))
    89        IF (.NOT.ALLOCATED(qext))         ALLOCATE(qext(nq))
    90        IF (.NOT.ALLOCATED(alpha_lift))   ALLOCATE(alpha_lift(nq))
    91        IF (.NOT.ALLOCATED(alpha_devil))  ALLOCATE(alpha_devil(nq))
    92        IF (.NOT.ALLOCATED(qextrhor))     ALLOCATE(qextrhor(nq))
    93        IF (.NOT.ALLOCATED(igcm_dustbin)) ALLOCATE(igcm_dustbin(nq))
    94        IF (.NOT.ALLOCATED(is_chim))      ALLOCATE(is_chim(nqtot))
     84       IF (.NOT.ALLOCATED(noms))           ALLOCATE(noms(nq))
     85       IF (.NOT.ALLOCATED(mmol))           ALLOCATE(mmol(nq))
     86       IF (.NOT.ALLOCATED(aki))            ALLOCATE(aki(nqtot))
     87       IF (.NOT.ALLOCATED(cpi))            ALLOCATE(cpi(nqtot))
     88       IF (.NOT.ALLOCATED(radius))         ALLOCATE(radius(nq))
     89       IF (.NOT.ALLOCATED(rho_q))          ALLOCATE(rho_q(nq))
     90       IF (.NOT.ALLOCATED(qext))           ALLOCATE(qext(nq))
     91       IF (.NOT.ALLOCATED(alpha_lift))     ALLOCATE(alpha_lift(nq))
     92       IF (.NOT.ALLOCATED(alpha_devil))    ALLOCATE(alpha_devil(nq))
     93       IF (.NOT.ALLOCATED(qextrhor))       ALLOCATE(qextrhor(nq))
     94       IF (.NOT.ALLOCATED(igcm_dustbin))   ALLOCATE(igcm_dustbin(nq))
     95       IF (.NOT.ALLOCATED(is_chim))        ALLOCATE(is_chim(nqtot))
     96       IF (.NOT.ALLOCATED(is_rad))         ALLOCATE(is_rad(nqtot))
     97       IF (.NOT.ALLOCATED(is_recomb))      ALLOCATE(is_recomb(nqtot))
     98       IF (.NOT.ALLOCATED(is_recomb_qset)) THEN
     99         ALLOCATE(is_recomb_qset(nqtot))
     100       ENDIF
     101       IF (.NOT.ALLOCATED(is_recomb_qotf)) THEN
     102         ALLOCATE(is_recomb_qotf(nqtot))
     103       ENDIF
    95104       !! initialization
    96        alpha_lift(:)  = 0.
    97        alpha_devil(:) = 0.
    98        mmol(:)        = 0.
    99        aki(:)         = 0.
    100        cpi(:)         = 0.
    101        is_chim(:)     = 0
     105       alpha_lift(:)     = 0.
     106       alpha_devil(:)    = 0.
     107       mmol(:)           = 0.
     108       aki(:)            = 0.
     109       cpi(:)            = 0.
     110       is_chim(:)        = 0
     111       is_rad(:)         = 0
     112       is_recomb(:)      = 0
     113       is_recomb_qset(:) = 0
     114       is_recomb_qotf(:) = 0
    102115       
    103116       ! Added by JVO 2017 : these arrays are handled later
     
    393406      write(*,*) 'Number of species in the chemistry nesp = ',nesp
    394407
     408!     Processing modern traceur options
     409      if(moderntracdef) then
     410        call ini_recombin
     411      endif
     412
    395413c------------------------------------------------------------
    396414c     Initialisation tracers ....
     
    520538     $             is_chim(iq)
    521539          end if
     540          ! option is_rad
     541          if (index(tracline,'is_rad=') /= 0) then
     542              read(tracline(index(tracline,'is_rad=')
     543     $         +len('is_rad='):),*) is_rad(iq)
     544              write(*,*) ' Parameter value (traceur.def) : is_rad=',
     545     $         is_rad(iq)
     546          else
     547              write(*,*) ' Parameter value (default)     : is_rad=',
     548     $         is_rad(iq)
     549          end if
     550          ! option is_recomb
     551          if (index(tracline,'is_recomb=') /= 0) then
     552              read(tracline(index(tracline,'is_recomb=')
     553     $         +len('is_recomb='):),*) is_recomb(iq)
     554              write(*,*) ' Parameter value (traceur.def) : is_recomb=',
     555     $         is_recomb(iq)
     556          else
     557              write(*,*) ' Parameter value (default)     : is_recomb=',
     558     $         is_recomb(iq)
     559          end if
     560          ! option is_recomb_qset
     561          if (index(tracline,'is_recomb_qset=') /= 0) then
     562              read(tracline(index(tracline,'is_recomb_qset=')
     563     $         +len('is_recomb_qset='):),*) is_recomb_qset(iq)
     564              write(*,*) ' Parameter value (traceur.def) :'//
     565     $         ' is_recomb_qset=',
     566     $         is_recomb_qset(iq)
     567          else
     568              write(*,*) ' Parameter value (default)     :'//
     569     $         ' is_recomb_qset=',
     570     $         is_recomb_qset(iq)
     571          endif
     572          ! option is_recomb_qotf
     573          if (index(tracline,'is_recomb_qotf=') /= 0) then
     574              read(tracline(index(tracline,'is_recomb_qotf=')
     575     $         +len('is_recomb_qotf='):),*) is_recomb_qotf(iq)
     576              write(*,*) ' Parameter value (traceur.def) :'//
     577     $         ' is_recomb_qotf=',
     578     $         is_recomb_qotf(iq)
     579          else
     580              write(*,*) ' Parameter value (default)     :'//
     581     $         ' is_recomb_qotf=',
     582     $         is_recomb_qotf(iq)
     583          end if
    522584      end subroutine get_tracdat
    523585
  • trunk/LMDZ.GENERIC/libf/phystd/optci.F90

    r2520 r2543  
    1515  use comcstfi_mod, only: g, r, mugaz
    1616  use callkeys_mod, only: kastprof,continuum,graybody
     17  use recombin_corrk_mod, only: corrk_recombin, gasi_recomb
    1718  implicit none
    1819
     
    247248              ! transfer on the tested simulations !
    248249
    249               tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     250              IF (corrk_recombin) THEN ! added by JVO
     251                tmpk = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
     252              ELSE
     253                tmpk = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     254              ENDIF
    250255
    251256              KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
     
    256261           else
    257262
    258               tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
     263              IF (corrk_recombin) THEN ! added by JVO
     264                tmpkvar = GASI_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
     265              ELSE
     266                tmpkvar = GASI(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
     267              ENDIF
    259268
    260269              KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) *  &
  • trunk/LMDZ.GENERIC/libf/phystd/optcv.F90

    r2520 r2543  
    1414  use comcstfi_mod, only: g, r, mugaz
    1515  use callkeys_mod, only: kastprof,continuum,graybody,callgasvis
     16  use recombin_corrk_mod, only: corrk_recombin, gasv_recomb
    1617
    1718  implicit none
     
    245246              ! transfer on the tested simulations !
    246247
    247               tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     248              IF (corrk_recombin) THEN ! Added by JVO
     249                tmpk = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG) ! contains the mix of recombined species
     250              ELSE
     251                tmpk = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,1,NW,NG)
     252              ENDIF
    248253             
    249254              KCOEF(1) = tmpk(1,1) ! KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
     
    254259           else
    255260
    256               tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
     261              IF (corrk_recombin) THEN
     262                tmpkvar = GASV_RECOMB(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
     263              ELSE
     264                tmpkvar = GASV(MT(K):MT(K)+1,MP(K):MP(K)+1,NVAR(K):NVAR(K)+1,NW,NG)
     265              ENDIF
    257266
    258267              KCOEF(1) = tmpkvar(1,1,1) + WRATIO(K) *  &
  • trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90

    r2520 r2543  
    1414!     Added double gray case by Jeremy Leconte (2012)
    1515!     New HITRAN continuum data section by RW (2012)
     16!     Modern traceur.def & corrk recombing scheme by J.Vatant d'Ollone (2020)
    1617!
    1718!     Summary
     
    3334      use callkeys_mod, only: varactive,varfixed,graybody,callgasvis,&
    3435                continuum
     36      use tracer_h, only : nqtot, moderntracdef, is_recomb, noms
     37      use recombin_corrk_mod, only: su_recombin,        &
     38                corrk_recombin, use_premix, nrecomb_tot, rcb2tot_idx
    3539      implicit none
    3640
     
    5862
    5963      integer :: dummy
     64
     65      if (.not. moderntracdef) use_premix=.true. ! Added by JVO for compatibility with 'old' traceur.def
     66     
     67      if (use_premix) then ! use_premix flag added by JVO, thus if pure recombining then premix is skipped
    6068
    6169!=======================================================================
     
    8189      read(111,*) ngas
    8290
    83       if(ngas.ne.ngasmx)then
    84          print*,'Number of gases in radiative transfer data (',ngas,') does not', &
    85                 'match that in gases.def (',ngasmx,'), exiting.'
    86          call abort
    87       endif
     91      if(moderntracdef) then
     92           ! JVO 20 - TODO : Sanity check with nspcrad !
     93           print*, 'Warning : Sanity check on # of gases still not implemented !!'
     94      else
     95        if(ngas.ne.ngasmx)then
     96           print*,'Number of gases in radiative transfer data (',ngas,') does not', &
     97                  'match that in gases.def (',ngasmx,'), exiting.'
     98           call abort
     99        endif
     100      endif
    88101
    89102      if(ngas.eq.1 .and. (varactive.or.varfixed))then
     
    104117      do igas=1,ngas
    105118         read(111,*) gastype(igas)
    106          print*,'Gas ',igas,' is ',gastype(igas)
     119         if(corrk_recombin) then
     120            print*,'Premix gas ',igas,' is ',gastype(igas)
     121         else
     122           print*,'Gas ',igas,' is ',gastype(igas)
     123         endif
    107124      enddo
    108125
     
    121138      endif
    122139
    123       ! Check that gastype and gnom match
    124       do igas=1,ngas
    125          print*,'Gas ',igas,' is ',trim(gnom(igas))
    126          if (trim(gnom(igas)).ne.trim(gastype(igas))) then
    127             print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
    128                  'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', &
    129                  'gases.def with Q.dat in your radiative transfer directory.'
    130             call abort
    131          endif
    132       enddo
    133       print*,'Confirmed gas match in radiative transfer and gases.def!'
     140      if (moderntracdef) then
     141          ! JVO 20 - TODO : Sanity check with nspcrad !
     142          print*, 'Warning : Sanity check on name of gases still not implemented !!'
     143      else
     144        ! Check that gastype and gnom match
     145        do igas=1,ngas
     146           print*,'Gas ',igas,' is ',trim(gnom(igas))
     147           if (trim(gnom(igas)).ne.trim(gastype(igas))) then
     148              print*,'Name of a gas in radiative transfer data (',trim(gastype(igas)),') does not ', &
     149                   'match that in gases.def (',trim(gnom(igas)),'), exiting. You should compare ', &
     150                   'gases.def with Q.dat in your radiative transfer directory.'
     151              call abort
     152           endif
     153        enddo
     154        print*,'Confirmed gas match in radiative transfer and gases.def!'
     155      endif
    134156
    135157      ! display the values
     
    140162      end do
    141163      print*,''
     164     
     165      else
     166        L_REFVAR = 1
     167      endif ! use_premix
    142168
    143169!=======================================================================
     
    258284      tgasmax = tgasref(L_NTREF)
    259285
    260       IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
    261       IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
     286      if(corrk_recombin) then ! even if no premix we keep a L_REFVAR=1 to store precombined at firstcall (see
     287         IF(.NOT. ALLOCATED(gasi8)) ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
     288         IF(.NOT. ALLOCATED(gasv8)) ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
     289      else
     290         IF(.NOT. ALLOCATED(gasi8)) ALLOCATE(gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS))
     291         IF(.NOT. ALLOCATED(gasv8)) ALLOCATE(gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS))
     292      endif
    262293!$OMP END MASTER
    263294!$OMP BARRIER
    264295
     296! JVO 20 : In these gasi/gasi8 matrix  we stack in same dim. variable specie and species to recombine (to keep code small)
     297
    265298!-----------------------------------------------------------------------
    266299! allocate the multidimensional arrays in radcommon_h
    267       IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
    268       IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
    269 
    270       ! display the values
    271       print*,''
    272       print*,'Correlated-k matrix size:'
    273       print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
     300     if(corrk_recombin) then
     301        IF(.NOT. ALLOCATED(gasi)) ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTI,L_NGAUSS))
     302        IF(.NOT. ALLOCATED(gasv)) ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR+nrecomb_tot,L_NSPECTV,L_NGAUSS))
     303        ! display the values
     304        print*,''
     305        print*,'Correlated-k matrix size:'
     306        print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR+nrecomb_tot,' (',L_REFVAR,'+',nrecomb_tot,'),',L_NGAUSS,']'
     307      else
     308        IF(.NOT. ALLOCATED(gasi)) ALLOCATE(gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS))
     309        IF(.NOT. ALLOCATED(gasv)) ALLOCATE(gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS))
     310        ! display the values
     311        print*,''
     312        print*,'Correlated-k matrix size:'
     313        print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
     314      endif
    274315
    275316!=======================================================================
     
    319360
    320361!$OMP MASTER         
    321 !      print*,corrkdir(1:4)
     362
    322363      ! VISIBLE
    323364      if (callgasvis) then
    324          if ((corrkdir(1:4).eq.'null'))then   !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
    325             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
    326             print*,'using no corrk data'
    327             print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
    328          else
    329             file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
     365           
     366        ! Looking for premixed corrk files corrk_gcm_VI.dat if needed
     367        if (use_premix) then
     368           ! print*,corrkdir(1:4)
     369           if ((corrkdir(1:4).eq.'null'))then   !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
     370              gasv8(:,:,:,:,:)=0.0
     371              print*,'using no corrk data'
     372              print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
     373           else
     374              file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
     375              file_path=TRIM(datadir)//TRIM(file_id)
     376             
     377              ! check that the file exists
     378              inquire(FILE=file_path,EXIST=file_ok)
     379              if(.not.file_ok) then
     380                 write(*,*)'The file ',TRIM(file_path)
     381                 write(*,*)'was not found by sugas_corrk.F90.'
     382                 write(*,*)'Are you sure you have absorption data for these bands?'
     383                 call abort
     384              endif
     385           
     386              open(111,file=TRIM(file_path),form='formatted')
     387              read(111,*) gasv8(:,:,:L_REFVAR,:,:)
     388              close(111)
     389           end if
     390        else
     391           gasv8(:,:,1,:,:)=0.0 ! dummy but needs to be initialized
     392        endif ! use_premix
     393       
     394        ! Looking for others files corrk_gcm_VI_XXX.dat if needed
     395        if (corrk_recombin) then
     396          do igas=1,nrecomb_tot
     397
     398            file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI_'//trim(adjustl(noms(rcb2tot_idx(igas))))//'.dat'
    330399            file_path=TRIM(datadir)//TRIM(file_id)
    331            
     400             
    332401            ! check that the file exists
    333402            inquire(FILE=file_path,EXIST=file_ok)
     
    335404               write(*,*)'The file ',TRIM(file_path)
    336405               write(*,*)'was not found by sugas_corrk.F90.'
    337                write(*,*)'Are you sure you have absorption data for these bands?'
     406               write(*,*)'Are you sure you have absorption data for this specie at these bands?'
    338407               call abort
    339408            endif
    340409         
    341410            open(111,file=TRIM(file_path),form='formatted')
    342             read(111,*) gasv8
     411            read(111,*) gasv8(:,:,L_REFVAR+igas,:,:)
    343412            close(111)
    344          end if
    345 
    346          if(nVI_limit.eq.0) then
    347             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=   &
    348                   gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
    349          else if (nVI_limit.eq.L_NSPECTV) then
    350             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=   &
    351                   gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_IR
    352          else
    353             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=   &
    354                   gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)+kappa_IR
    355             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=   &
    356                   gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
    357          end if
     413          enddo                           
     414        endif ! corrk_recombin
     415
     416        if(nVI_limit.eq.0) then
     417         gasv8(:,:,:,:,:)= gasv8(:,:,:,:,:)+kappa_VI
     418           else if (nVI_limit.eq.L_NSPECTV) then
     419         gasv8(:,:,:,:,:)= gasv8(:,:,:,:,:)+kappa_IR
    358420      else
    359          print*,'Visible corrk gaseous absorption is set to zero.'
    360          gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
    361       endif
     421         gasv8(:,:,:,1:nVI_limit,:)= gasv8(:,:,:,1:nVI_limit,:)+kappa_IR
     422         gasv8(:,:,:,nVI_limit+1:L_NSPECTV,:)= gasv8(:,:,:,nVI_limit+1:L_NSPECTV,:)+kappa_VI
     423      end if
     424           
     425         else
     426            print*,'Visible corrk gaseous absorption is set to zero.'
     427            gasv8(:,:,:,:,:)=0.0
     428         endif ! callgasvis
     429         
    362430!$OMP END MASTER
    363431!$OMP BARRIER
    364432
    365433      ! INFRA-RED
    366       if ((corrkdir(1:4).eq.'null'))then       !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
    367          print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
     434     
     435      ! Looking for premixed corrk files corrk_gcm_IR.dat if needed
     436      if (use_premix) then
     437        if ((corrkdir(1:4).eq.'null'))then       !.or.(TRIM(corrkdir).eq.'null_LowTeffStar')) then
     438           print*,'Infrared corrk gaseous absorption is set to zero if graybody=F'
    368439!$OMP MASTER         
    369          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
     440           gasi8(:,:,:,:,:)=0.0
    370441!$OMP END MASTER
    371442!$OMP BARRIER
    372       else
    373          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
    374          file_path=TRIM(datadir)//TRIM(file_id)
    375      
    376          ! check that the file exists
    377          inquire(FILE=file_path,EXIST=file_ok)
    378          if(.not.file_ok) then
    379             write(*,*)'The file ',TRIM(file_path)
    380             write(*,*)'was not found by sugas_corrk.F90.'
    381             write(*,*)'Are you sure you have absorption data for these bands?'
    382             call abort
    383          endif
     443        else
     444           file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
     445           file_path=TRIM(datadir)//TRIM(file_id)
     446       
     447           ! check that the file exists
     448           inquire(FILE=file_path,EXIST=file_ok)
     449           if(.not.file_ok) then
     450              write(*,*)'The file ',TRIM(file_path)
     451              write(*,*)'was not found by sugas_corrk.F90.'
     452              write(*,*)'Are you sure you have absorption data for these bands?'
     453              call abort
     454           endif
    384455         
    385456!$OMP MASTER                   
    386          open(111,file=TRIM(file_path),form='formatted')
    387          read(111,*) gasi8
    388          close(111)
     457           open(111,file=TRIM(file_path),form='formatted')
     458           read(111,*) gasi8(:,:,:L_REFVAR,:,:)
     459           close(111)
    389460!$OMP END MASTER
    390461!$OMP BARRIER
    391462     
    392          ! 'fzero' is a currently unused feature that allows optimisation
    393          ! of the radiative transfer by neglecting bands where absorption
    394          ! is close to zero. As it could be useful in the future, this
    395          ! section of the code has been kept commented and not erased.
    396          ! RW 7/3/12.
    397 
    398          do nw=1,L_NSPECTI
    399             fzeroi(nw) = 0.d0
    400 !            do nt=1,L_NTREF
    401 !               do np=1,L_NPREF
    402 !                  do nh=1,L_REFVAR
    403 !                     do ng = 1,L_NGAUSS
    404 !                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
    405 !                           fzeroi(nw)=fzeroi(nw)+1.d0
    406 !                        endif
    407 !                     end do
    408 !                  end do
    409 !               end do
    410 !            end do
    411 !            fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
    412          end do
    413 
    414          do nw=1,L_NSPECTV
    415             fzerov(nw) = 0.d0
    416 !            do nt=1,L_NTREF
    417 !               do np=1,L_NPREF
    418 !                  do nh=1,L_REFVAR
    419 !                     do ng = 1,L_NGAUSS
    420 !                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
    421 !                           fzerov(nw)=fzerov(nw)+1.d0
    422 !                        endif
    423 !                     end do
    424 !                  end do
    425 !               end do
    426 !            end do
    427 !            fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
    428          end do
    429 
    430       endif
     463           ! 'fzero' is a currently unused feature that allows optimisation
     464           ! of the radiative transfer by neglecting bands where absorption
     465           ! is close to zero. As it could be useful in the future, this
     466           ! section of the code has been kept commented and not erased.
     467           ! RW 7/3/12.
     468
     469           do nw=1,L_NSPECTI
     470              fzeroi(nw) = 0.d0
     471!              do nt=1,L_NTREF
     472!                 do np=1,L_NPREF
     473!                    do nh=1,L_REFVAR
     474!                       do ng = 1,L_NGAUSS
     475!                          if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
     476!                             fzeroi(nw)=fzeroi(nw)+1.d0
     477!                          endif
     478!                       end do
     479!                    end do
     480!                 end do
     481!              end do
     482!              fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
     483           end do
     484
     485           do nw=1,L_NSPECTV
     486              fzerov(nw) = 0.d0
     487!              do nt=1,L_NTREF
     488!                 do np=1,L_NPREF
     489!                    do nh=1,L_REFVAR
     490!                       do ng = 1,L_NGAUSS
     491!                          if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
     492!                             fzerov(nw)=fzerov(nw)+1.d0
     493!                          endif
     494!                       end do
     495!                    end do
     496!                 end do
     497!              end do
     498!              fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
     499           end do
     500
     501        endif ! if corrkdir=null
     502
     503      else
     504         gasi8(:,:,1,:,:)=0.0 ! dummy but needs to be initialized
     505      endif ! use_premix
     506   
     507      ! Looking for others files corrk_gcm_IR_XXX.dat if needed
     508      if (corrk_recombin) then
     509        do igas=1,nrecomb_tot
     510
     511          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR_'//trim(adjustl(noms(rcb2tot_idx(igas))))//'.dat'
     512          file_path=TRIM(datadir)//TRIM(file_id)
     513           
     514          ! check that the file exists
     515          inquire(FILE=file_path,EXIST=file_ok)
     516          if(.not.file_ok) then
     517             write(*,*)'The file ',TRIM(file_path)
     518             write(*,*)'was not found by sugas_corrk.F90.'
     519             write(*,*)'Are you sure you have absorption data for this specie at these bands?'
     520             call abort
     521          endif
     522!$OMP MASTER           
     523          open(111,file=TRIM(file_path),form='formatted')
     524          read(111,*) gasi8(:,:,L_REFVAR+igas,:,:)
     525          close(111)
     526!$OMP END MASTER
     527!$OMP BARRIER
     528        enddo                           
     529      endif ! corrk_recombin
    431530
    432531!$OMP MASTER                 
    433532      if(nIR_limit.eq.0) then
    434          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=   &
    435                   gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
     533         gasi8(:,:,:,:,:)= gasi8(:,:,:,:,:)+kappa_VI
    436534      else if (nIR_limit.eq.L_NSPECTI) then
    437          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=   &
    438                   gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_IR
     535         gasi8(:,:,:,:,:)= gasi8(:,:,:,:,:)+kappa_IR
    439536      else
    440          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=   &
    441                   gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)+kappa_IR
    442          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=   &
    443                   gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
     537         gasi8(:,:,:,1:nIR_limit,:)= gasi8(:,:,:,1:nIR_limit,:)+kappa_IR
     538         gasi8(:,:,:,nIR_limit+1:L_NSPECTI,:)= gasi8(:,:,:,nIR_limit+1:L_NSPECTI,:)+kappa_VI
    444539      end if
    445540
     
    614709      end do
    615710
     711! Allocate and initialise arrays for corrk_recombin
     712      if (corrk_recombin) call su_recombin
    616713
    617714!=======================================================================
  • trunk/LMDZ.GENERIC/libf/phystd/tracer_h.F90

    r2542 r2543  
    3232
    3333       integer, save, allocatable :: is_chim(:) ! 1 if tracer used in chemistry, else 0
    34 !$OMP THREADPRIVATE(is_chim)
     34       integer, save, allocatable :: is_rad(:)  ! 1 if   ""    ""  in radiative transfer, else 0
     35!$OMP THREADPRIVATE(is_chim,is_rad)
     36
     37       integer, save, allocatable :: is_recomb(:)      ! 1 if tracer used in recombining scheme, else 0 (if 1, must have is_rad=1)
     38       integer, save, allocatable :: is_recomb_qset(:) ! 1 if tracer k-corr assume predefined vmr, else 0 (if 1, must have is_recomb=1)
     39       integer, save, allocatable :: is_recomb_qotf(:) ! 1 if tracer recombination is done on-the-fly, else 0 (if 1, must have is_recomb_qset=0)
     40!$OMP THREADPRIVATE(is_recomb,is_recomb_qset,is_recomb_qotf)
    3541
    3642! tracer indexes: these are initialized in initracer and should be 0 if the
Note: See TracChangeset for help on using the changeset viewer.