Ignore:
Timestamp:
Dec 6, 2013, 1:01:19 AM (11 years ago)
Author:
jleconte
Message:

05/12/2013 == JL

  • corrected sugascorrk to work in the two band gray aproximation with -b 1x1 and NGAUSS=2
File:
1 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/libf/phystd/sugas_corrk.F90

    r961 r1122  
    5050      real*8 x, xi(4), yi(4), ans, p
    5151!     For gray case (JL12)
    52       real kappa_IR, kappa_VI, IR_VI_wnlimit, nVI_limit,nIR_limit
     52      real kappa_IR, kappa_VI, IR_VI_wnlimit
     53      integer nVI_limit,nIR_limit
    5354
    5455      integer ngas, igas, jgas
    … …  
    102103      do igas=1,ngas
    103104         read(111,*) gastype(igas)
    104          print*,'Gas ',igas,' is ',gastype(igas)
     105         print*,'Gas ',igas,' is ',gastype(igas+      real kappa_IR, kappa_VI, IR_VI_wnlimit
     106)
    105107      enddo
    106108
    … …  
    276278!=======================================================================
    277279!     Get gaseous k-coefficients and interpolate onto finer pressure grid
     280
     281
     282!        wavelength used to separate IR from VI in graybody. We will need that anyway
     283         IR_VI_wnlimit=3000.
     284         write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
     285         
     286         nVI_limit=0
     287         Do nw=1,L_NSPECTV
     288            if ((WNOV(nw).gt.IR_VI_wnlimit).and.(L_NSPECTV.gt.1)) then
     289               nVI_limit=nw-1
     290               exit
     291            endif
     292         End do
     293         nIR_limit=L_NSPECTI
     294         Do nw=1,L_NSPECTI
     295            if ((WNOI(nw).gt.IR_VI_wnlimit).and.(L_NSPECTI.gt.1)) then
     296               nIR_limit=nw-1
     297               exit
     298            endif
     299         End do
    278300
    279301      if (graybody) then
    … …  
    291313         write(*,*)" kappa_IR = ",kappa_IR       
    292314         kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule
     315
     316         write(*,*)"graybody: Visible / Infrared separation set at band: IR=",nIR_limit,", VI=",nVI_limit
    293317               
    294 !        wavelength used to separate IR from VI
    295          IR_VI_wnlimit=3000.
    296          write(*,*)"graybody: Visible / Infrared separation set at",10000./IR_VI_wnlimit,"um"
    297          
    298          nVI_limit=L_NSPECTV
    299          Do nw=1,L_NSPECTV
    300             if (WNOV(nw).gt.IR_VI_wnlimit) then
    301                nVI_limit=nw-1
    302                exit
    303             endif
    304          End do
    305          nIR_limit=L_NSPECTI
    306          Do nw=1,L_NSPECTI
    307             if (WNOI(nw).gt.IR_VI_wnlimit) then
    308                nIR_limit=nw-1
    309                exit
    310             endif
    311          End do
    312          write(*,*)"graybody: Visible / Infrared separation set at band ",nIR_limit,nVI_limit
    313              
     318      Else
     319         kappa_VI=1.e-30     
     320         kappa_IR=1.e-30       
    314321      End if
    315322
    316 
     323!      print*,corrkdir(1:4)
    317324      ! VISIBLE
    318       if (callgasvis.and.(.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
    319          file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
    320          file_path=TRIM(datadir)//TRIM(file_id)
     325      if (callgasvis) then
     326         if ((corrkdir(1:4).eq.'null'))then   !(TRIM(corrkdir).eq.'null_LowTeffStar')) then
     327            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=0.0
     328            print*,'using no corrk data'
     329            print*,'Visible corrk gaseous absorption is set to zero if graybody=F'
     330         else
     331            file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_VI.dat'
     332            file_path=TRIM(datadir)//TRIM(file_id)
     333           
     334            ! check that the file exists
     335            inquire(FILE=file_path,EXIST=file_ok)
     336            if(.not.file_ok) then
     337               write(*,*)'The file ',TRIM(file_path)
     338               write(*,*)'was not found by sugas_corrk.F90.'
     339               write(*,*)'Are you sure you have absorption data for these bands?'
     340               call abort
     341            endif
    321342         
    322          ! check that the file exists
    323          inquire(FILE=file_path,EXIST=file_ok)
    324          if(.not.file_ok) then
    325             write(*,*)'The file ',TRIM(file_path)
    326             write(*,*)'was not found by sugas_corrk.F90.'
    327             write(*,*)'Are you sure you have absorption data for these bands?'
    328             call abort
    329          endif
    330          
    331          open(111,file=TRIM(file_path),form='formatted')
    332          read(111,*) gasv8
    333          close(111)
    334 
    335       else if (callgasvis.and.graybody) then
     343            open(111,file=TRIM(file_path),form='formatted')
     344            read(111,*) gasv8
     345            close(111)
     346         end if
     347
    336348         if(nVI_limit.eq.0) then
    337             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
     349            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=   &
     350                  gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
    338351         else if (nVI_limit.eq.L_NSPECTV) then
    339             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=kappa_IR
     352            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)=   &
     353                  gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTV,1:L_NGAUSS)+kappa_IR
    340354         else
    341             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=kappa_IR
    342             gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=kappa_VI
     355            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)=   &
     356                  gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nVI_limit,1:L_NGAUSS)+kappa_IR
     357            gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)=   &
     358                  gasv8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nVI_limit+1:L_NSPECTV,1:L_NGAUSS)+kappa_VI
    343359         end if
    344360      else
    … …  
    348364
    349365      ! INFRA-RED
    350       if ((.not.TRIM(corrkdir).eq.'null').and.(.not.TRIM(corrkdir).eq.'null_LowTeffStar').and.(.not.graybody)) then
     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'
     368         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
     369      else 
    351370         file_id='/corrk_data/'//trim(adjustl(banddir))//'/corrk_gcm_IR.dat'
    352371         file_path=TRIM(datadir)//TRIM(file_id)
    … …  
    403422         end do
    404423
    405       else if (graybody) then
    406          if(nIR_limit.eq.0) then
    407             gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
    408          else if (nIR_limit.eq.L_NSPECTI) then
    409             gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=kappa_IR
    410          else
    411             gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=kappa_IR
    412             gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=kappa_VI
    413          end if
     424      endif
     425
     426      if(nIR_limit.eq.0) then
     427         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=   &
     428                  gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
     429      else if (nIR_limit.eq.L_NSPECTI) then
     430         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=   &
     431                  gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)+kappa_IR
    414432      else
    415          print*,'Infrared corrk gaseous absorption is set to zero.'
    416          gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:L_NSPECTI,1:L_NGAUSS)=0.0
    417       endif
     433         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)=   &
     434                  gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,1:nIR_limit,1:L_NGAUSS)+kappa_IR
     435         gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)=   &
     436                  gasi8(1:L_NTREF,1:L_NPREF,1:L_REFVAR,nIR_limit+1:L_NSPECTI,1:L_NGAUSS)+kappa_VI
     437      end if
    418438
    419439
Note: See TracChangeset for help on using the changeset viewer.