Changeset 716


Ignore:
Timestamp:
Jul 3, 2012, 8:09:54 PM (12 years ago)
Author:
rwordsworth
Message:

Mainly updates to radiative transfer and gas management scheme.
Most CIA data now read from standard HITRAN datafiles. For the H2O
continuum, two options have been added: the standard CKD continuum,
and the empirical formula in PPC (Pierrehumbert 2010). Use the toggle
'H2Ocont_simple' in callphys.def to choose.

Note to Martians: I've changed the default values of 'sedimentation' and
'co2cond' in inifis to false. Both these are defined in the standard deftank
callphys.def file, so there shouldn't be any compatibility problems.

Location:
trunk/LMDZ.GENERIC/libf/phystd
Files:
5 added
1 deleted
23 edited

Legend:

Unmodified
Added
Removed

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

    r526 r716  
    235235                          (  0.75 * QREFvis3d(ig,l,iaer) /        &
    236236                          ( rho_ice * reffrad(ig,l,iaer) )  ) *   &
    237                             pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
    238                           ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the
    239                             CLF1                                    !   clear=false/pq=0 case
     237                          !  pq(ig,l,i_h2oice) *                   & !JL I dropped the +1e-9 here to have the same
     238                          !( pplev(ig,l) - pplev(ig,l+1) ) / g /   & !   opacity in the clearsky=true and the
     239                          !  CLF1                                   !   clear=false/pq=0 case
     240                          ( pq(ig,l,i_h2oice) + 1.E-9 ) *         & ! Doing this makes the code unstable, so I have restored it (RW)
     241                          ( pplev(ig,l) - pplev(ig,l+1) ) / g /   & 
     242                            CLF1
    240243
    241244                     aerosol(ig,l,iaer) = -log(1 - CLF2 + CLF2*exp(-aerosol0))

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

    r305 r716  
    11!-------------------------------------------------------------------------
    2       subroutine bilinear(x_arr,y_arr,nX,nY,f2d_arr,x_in,y_in,f)
    3 !     Necessary for interpolation of continuum data
     2subroutine bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2) 
     3! Used for interpolation of continuum data
    44
    5       implicit none
     5  implicit none
    66
    7       integer nX,nY,i,j,a,b
     7  real*8 x,y,x1,x2,y1,y2
     8  real*8 f,f11,f12,f21,f22,fA,fB
    89
    9       real*8 x_in,y_in,x,y,x1,x2,y1,y2
    10       real*8 f,f11,f12,f21,f22,fA,fB
    11       real*8 x_arr(nX)
    12       real*8 y_arr(nY)
    13       real*8 f2d_arr(nX,nY)
    14      
    15       integer strlen
    16       character*100 label
    17       label='subroutine bilinear'
     10  ! 1st in x-direction
     11  fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
     12  fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
    1813
    19       x=x_in
    20       y=y_in
     14  ! then in y-direction
     15  f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
    2116
    22 !     1st check we're within the wavenumber range
    23       if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
    24          f=0.0D+0
    25          return
    26       else
    27          
    28 !     in the x (wavenumber) direction 1st
    29          i=1
    30  10      if (i.lt.(nX+1)) then
    31             if (x_arr(i).gt.x) then
    32                x1=x_arr(i-1)
    33                x2=x_arr(i)
    34                a=i-1
    35                i=9999
    36             endif
    37             i=i+1
    38             goto 10
    39          endif
    40       endif
    41      
    42       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
    43          write(*,*) 'Warning from bilinear.for:'
    44          write(*,*) 'Outside continuum temperature range!'
    45          if(y.lt.y_arr(1))then
    46             y=y_arr(1)+0.01
    47          endif
    48          if(y.gt.y_arr(nY))then
    49             y=y_arr(nY)-0.01
    50          endif
    51       else
    52 
    53 !     in the y (temperature) direction 2nd
    54          j=1
    55  20      if (j.lt.(nY+1)) then
    56             if (y_arr(j).gt.y) then
    57                y1=y_arr(j-1)
    58                y2=y_arr(j)
    59                b=j-1
    60                j=9999
    61             endif
    62             j=j+1
    63             goto 20
    64          endif
    65       endif
    66      
    67       f11=f2d_arr(a,b)
    68       f21=f2d_arr(a+1,b)
    69       f12=f2d_arr(a,b+1)
    70       f22=f2d_arr(a+1,b+1)
    71      
    72 !     1st in x-direction
    73       fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
    74       fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
    75      
    76 !     then in y-direction
    77       f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
    78      
    79       return
    80     end subroutine bilinear
     17  return
     18end subroutine bilinear

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

    r589 r716  
    4040            ! all values at 300 K from Engineering Toolbox
    4141            if(gnom(igas).eq.'CO2')then
     42               !cpp_c   = cpp_c   + 0.744*gfrac(igas) ! @ ~210 K (better for Mars conditions)
    4243               cpp_c   = cpp_c   + 0.846*gfrac(igas)
    4344               mugaz_c = mugaz_c + 44.01*gfrac(igas)
    … …  
    4849               cpp_c   = cpp_c   + 14.31*gfrac(igas)
    4950               mugaz_c = mugaz_c + 2.01*gfrac(igas)
     51            elseif(gnom(igas).eq.'He_')then
     52               cpp_c   = cpp_c   + 5.19*gfrac(igas)
     53               mugaz_c = mugaz_c + 4.003*gfrac(igas)
    5054            elseif(gnom(igas).eq.'H2O')then
    5155               cpp_c   = cpp_c   + 1.864*gfrac(igas)
    5256               mugaz_c = mugaz_c + 18.02*gfrac(igas)
     57            elseif(gnom(igas).eq.'SO2')then
     58               cpp_c   = cpp_c   + 0.64*gfrac(igas)
     59               mugaz_c = mugaz_c + 64.066*gfrac(igas)
     60            elseif(gnom(igas).eq.'H2S')then
     61               cpp_c   = cpp_c   + 1.003*gfrac(igas) ! from wikipedia...
     62               mugaz_c = mugaz_c + 34.08*gfrac(igas)
    5363            elseif(gnom(igas).eq.'CH4')then
    5464               cpp_c   = cpp_c   + 2.226*gfrac(igas)

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

    r471 r716  
    7474            elseif(gnom(igas).eq.'He_')then
    7575               print*,'Helium not ready yet!'
     76               tauconsti(igas) = 0.0
     77               call abort
    7678            else
    7779               print*,'Error in calc_rayleigh: Gas species not recognised!'
    … …  
    117119                  elseif(gnom(igas).eq.'H2_')then
    118120                     tauvari(igas) = 1.0/wl**4
     121                  elseif(gnom(igas).eq.'He_')then
     122                     print*,'Helium not ready yet!'
     123                     tauvari(igas) = 0.0
    119124                  else
    120125                     print*,'Error in calc_rayleigh: Gas species not recognised!'

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

    r650 r716  
    77          reffrad,tau_col,cloudfrac,totcloudfrac,              &
    88          clearsky,firstcall,lastcall)
    9 
    109
    1110      use radinc_h
    … …  
    201200      giaer(:,:,:) =0.0
    202201      radfixed=.false.
    203 
    204202
    205203      if(firstcall) then
    … …  
    275273         endif
    276274
    277          OLR_nu=0.
    278          OSR_nu=0.
     275         OLR_nu(:,:) = 0.
     276         OSR_nu(:,:) = 0.
    279277
    280278         if (ngridmx.eq.1) then
    … …  
    429427           reffrad,QREFvis3d,QREFir3d,                             &
    430428           tau_col,cloudfrac,totcloudfrac,clearsky)                ! get aerosol optical depths
    431 
    432429
    433430!-----------------------------------------------------------------------
    … …  
    933930      endif
    934931
     932
    935933    end subroutine callcorrk

  • trunk/LMDZ.GENERIC/libf/phystd/callkeys.h

    r650 r716  
    88     &   , iaervar,iddist,topdustref,callstats,calleofdump              &
    99     &   , enertest                                                     &
    10      &   , callgasvis,Continuum, graybody                               &
     10     &   , callgasvis,Continuum,H2Ocont_simple,graybody                 &
    1111     &   , Nmix_co2, Nmix_h2o                                           &
    1212     &   , dusttau                                                      &
    … …  
    5252     &   , season,diurnal,tlocked,lwrite                                &
    5353     &   , callstats,calleofdump                                        &
    54      &   , callgasvis,Continuum,graybody
     54     &   , callgasvis,Continuum,H2Ocont_simple,graybody
    5555
    5656      logical enertest

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

    r471 r716  
    88  double precision T
    99
     10
     11  ! this function has been disabled in gradients_kcm.F90 because it dont
     12  ! work if you have gaseous mixtures. need to decide whether to generalise
     13  ! it or simply remove entirely...
    1014
    1115  ! Cp_n : cf CO2 dans abe&matsui (1988)

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

    r374 r716  
    44      implicit none
    55
     6      ! Default for Berserker @ UChicago:
     7!      character(len=300) :: datadir='/home/rwordsworth/datagcm'
    68      ! Default for Gnome Idataplex:
    7       character(len=300) :: datadir='/san/home/rdword/gcm/datagcm'
     9!      character(len=300) :: datadir='/san/home/rdword/gcm/datagcm'
    810      ! Default for LMD machines:
    9 !      character(len=300) :: datadir='/u/rwlmd/datagcm'
     11      character(len=300) :: datadir='/u/rwlmd/datagcm'
    1012
    1113      end module datafile_mod

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

    r471 r716  
    55!======================================================================C
    66!
    7 !     GASES_H   
     7!     gases_h   
    88!
    9 !     THIS A F90-ALLOCATABLE VERSION FOR gases.h -- AS 12/2011
     9!     A F90-allocatable version for gases.h -- AS 12/2011
    1010!
    1111!======================================================================C
    1212
    13       !!! THOSE ARE SET AND ALLOCATED IN su_gases.F90
     13      ! Set and allocated in su_gases.F90
    1414      integer :: ngasmx
    1515      integer :: vgas
    … …  
    1717      real,allocatable,DIMENSION(:) :: gfrac
    1818
     19      ! in analogy with tracer.h ...
     20      integer :: igas_H2
     21      integer :: igas_He
     22      integer :: igas_H2O
     23      integer :: igas_CO2
     24      integer :: igas_CO
     25      integer :: igas_N2
     26      integer :: igas_O2
     27      integer :: igas_SO2
     28      integer :: igas_H2S
     29      integer :: igas_CH4
     30      integer :: igas_NH3
     31
    1932      end module gases_h

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

    r471 r716  
    1414
    1515  ! internal
    16   double precision cp_n,cp_v
     16  !double precision cp_n,cp_v
     17  double precision cp_v
    1718
    1819  double precision press, rho_plus, rho_minus, dVdT, rho_c
    … …  
    2627  double precision cp_neutral
    2728
    28   cp_n = cp_neutral(T)
    29 
     29  !cp_n = cp_neutral(T)
    3030
    3131  select case(profil_flag)

  • trunk/LMDZ.GENERIC/libf/phystd/inifis.F

    r650 r716  
    204204         write(*,*) " Continuum = ",Continuum
    205205
     206         write(*,*) "use analytic function for H2O continuum ?"
     207         H2Ocont_simple=.false. ! default value
     208         call getin("H2Ocont_simple",H2Ocont_simple)
     209         write(*,*) " H2Ocont_simple = ",H2Ocont_simple
    206210 
    207211         write(*,*) "call turbulent vertical diffusion ?"
    … …  
    221225
    222226         write(*,*) "call CO2 condensation ?"
    223          co2cond=.true. ! default value
     227         co2cond=.false. ! default value
    224228         call getin("co2cond",co2cond)
    225229         write(*,*) " co2cond = ",co2cond
    … …  
    228232            print*,'We need a CO2 ice tracer to condense CO2'
    229233            call abort
    230          endif 
    231        
     234         endif
     235
    232236         write(*,*) "CO2 supersaturation level ?"
    233237         co2supsat=1.0 ! default value
    … …  
    436440
    437441         write(*,*) "Gravitationnal sedimentation ?"
    438          sedimentation=.true. ! default value
     442         sedimentation=.false. ! default value
    439443         call getin("sedimentation",sedimentation)
    440444         write(*,*) " sedimentation = ",sedimentation

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

    r601 r716  
    66!     -------
    77!     Calculates the H2-H2 CIA opacity, using a lookup table from
    8 !     Borysow (2002)
     8!     HITRAN (2011)
    99!
    1010!     Authors
    1111!     -------
    12 !     R. Wordsworth (2009)
     12!     R. Wordsworth (2011)
    1313!     
    1414!==================================================================
    … …  
    2626
    2727      integer nS,nT
    28       parameter(nS=1649)
    29       parameter(nT=14)
     28      parameter(nS=2428)
     29      parameter(nT=10)
     30
     31      double precision, parameter :: losch = 2.6867774e19
     32      ! Loschmit's number (molecule cm^-3 at STP)
     33      ! converts cm^5 molecule^-2 --> cm^-1 amagat^-2
     34      ! see Richard et al. 2011, JQSRT for details
    3035
    3136      double precision amagat
    … …  
    3338      double precision temp_arr(nT)
    3439      double precision abs_arr(nS,nT)
    35       double precision data_tmp(nT/2+1)
    36 
    37       integer k
     40
     41      integer k,iT
    3842      logical firstcall
    3943
    … …  
    4347      integer strlen,ios
    4448
    45       amagat=(273.15/temp)*(pres/101325.0)
    46 
    47       if(firstcall)then
     49      character(LEN=*), parameter :: fmat1 = "(A20,F10.3,F10.3,I7,F7.1,E10.3,F5.3)"
     50
     51      character*20 bleh
     52      double precision blah, Ttemp
     53      integer nres
     54
     55      if(temp.gt.400)then
     56         print*,'Your temperatures are too high for this H2-H2 CIA dataset. If you '
     57         print*,'really want to run simulations with hydrogen at T > 400 K, contact'
     58         print*,'Robin Wordsworth [rwordsworth@uchicago.edu].'
     59         stop
     60      endif
     61
     62      amagat = (273.15/temp)*(pres/101325.0)
     63
     64      if(firstcall)then ! called by sugas_corrk only
     65         print*,'----------------------------------------------------'
     66         print*,'Initialising H2-H2 continuum from HITRAN database...'
    4867
    4968!     1.1 Open the ASCII files
    50 
    51          ! cold
    52          dt_file=TRIM(datadir)//'/continuum_data/H2H2_CIA_LT.dat'
     69         dt_file=TRIM(datadir)//'/continuum_data/H2-H2_norm_2011.cia'
     70
    5371         open(33,file=dt_file,form='formatted',status='old',iostat=ios)
    5472         if (ios.ne.0) then        ! file not found
    5573           write(*,*) 'Error from interpolateH2H2'
    5674           write(*,*) 'Data file ',trim(dt_file),' not found.'
    57            write(*,*)'Check that your path to datagcm:',trim(datadir)
    58            write(*,*)' is correct. You can change it in callphys.def with:'
    59            write(*,*)' datadir = /absolute/path/to/datagcm'
    60            write(*,*)' Also check that there is a continuum_data/H2H2_CIA_LT.dat there.'
     75           write(*,*) 'Check that your path to datagcm:',trim(datadir)
     76           write(*,*) 'is correct. You can change it in callphys.def with:'
     77           write(*,*) 'datadir = /absolute/path/to/datagcm'
     78           write(*,*) 'Also check that the continuum data continuum_data/H2-H2_norm_2011.cia is there.'
    6179           call abort
    6280         else
    63             do k=1,nS
    64                read(33,*) data_tmp
    65                wn_arr(k)=data_tmp(1)
    66                abs_arr(k,1:7)=data_tmp(2:8)
     81
     82            do iT=1,nT
     83
     84               read(33,fmat1) bleh,blah,blah,nres,Ttemp
     85               if(nS.ne.nres)then
     86                  print*,'Resolution given in file: ',trim(dt_file)
     87                  print*,'is ',nres,', which does not match nS.'
     88                  print*,'Please adjust nS value in interpolateH2H2.F90'
     89                  stop
     90               endif
     91               temp_arr(iT)=Ttemp
     92
     93               do k=1,nS
     94                  read(33,*) wn_arr(k),abs_arr(k,it)
     95               end do
     96
    6797            end do
     98     
    6899         endif
    69100         close(33)
    70 
    71          ! hot
    72          dt_file=TRIM(datadir)//'/continuum_data/H2H2_CIA_HT.dat'
    73          open(34,file=dt_file,form='formatted',status='old',iostat=ios)
    74          if (ios.ne.0) then        ! file not found
    75            write(*,*) 'Error from interpolateH2H2'
    76            write(*,*) 'Data file ',trim(dt_file),' not found.'
    77            write(*,*)'Check that your path to datagcm:',trim(datadir)
    78            write(*,*)' is correct. You can change it in callphys.def with:'
    79            write(*,*)' datadir = /absolute/path/to/datagcm'
    80            write(*,*)' Also check that there is a continuum_data/H2H2_CIA_HT.dat there.'
    81            call abort
    82          else
    83             do k=1,nS
    84                read(34,*) data_tmp
    85                wn_arr(k)=data_tmp(1)
    86                ! wn_arr is identical
    87                abs_arr(k,8:14)=data_tmp(2:8)
    88             end do
    89          endif
    90          close(34)
    91 
    92          temp_arr(1)  = 60.
    93          temp_arr(2)  = 100.
    94          temp_arr(3)  = 150.
    95          temp_arr(4)  = 200.
    96          temp_arr(5)  = 250.
    97          temp_arr(6)  = 300.
    98          temp_arr(7)  = 350.
    99          temp_arr(8)  = 400.
    100          temp_arr(9)  = 500.
    101          temp_arr(10) = 600.
    102          temp_arr(11) = 700.
    103          temp_arr(12) = 800.
    104          temp_arr(13) = 900.
    105          temp_arr(14) = 1000.
    106 
    107101
    108102         print*,'interpolateH2H2: At wavenumber ',wn,' cm^-1'
    … …  
    112106         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    113107
    114          print*,'the absorption is ',abcoef,' cm^-1 amg^-2'
    115 
    116          abcoef=abcoef*100.0*amagat**2 ! convert to m^-1
    117 
    118          print*,'We have ',amagat,' amagats'
     108         print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
     109         print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
     110
     111         abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
     112
     113         print*,'We have ',amagat,' amagats of H2'
    119114         print*,'So the absorption is ',abcoef,' m^-1'
     115
     116               !open(117,file='T_array.dat')
     117               !do iT=1,nT
     118               !   write(117,*), temp_arr(iT)
     119               !end do
     120               !close(117)
     121
     122               !open(115,file='nu_array.dat')
     123               !do k=1,nS
     124               !   write(115,*), wn_arr(k)
     125               !end do
     126               !close(115)
     127
     128               !open(113,file='abs_array.dat')
     129               !do iT=1,nT
     130               !   do k=1,nS
     131               !      write(113,*), abs_arr(k,iT)*losch**2
     132               !   end do
     133               !end do
     134               !close(113)
    120135
    121136      else
    122137 
    123138         call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    124          abcoef=abcoef*100.0*amagat**2 ! convert to m^-1
    125          !print*,'The absorption is ',abcoef,' m^-1'
     139         abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
    126140
    127141         ! unlike for Rayleigh scattering, we do not currently weight by the BB function
    … …  
    141155
    142156      integer nX,nY,i,j,a,b
    143       parameter(nX=1649)
    144       parameter(nY=14)
     157      parameter(nX=2428)
     158      parameter(nY=10)
    145159
    146160      real*8 x_in,y_in,x,y,x1,x2,y1,y2
    … …  
    156170      x=x_in
    157171      y=y_in
     172
    158173
    159174!     1st check we're within the wavenumber range
    … …  
    180195         write(*,*) 'Warning from bilinearH2H2:'
    181196         write(*,*) 'Outside continuum temperature range!'
    182          write(*,*) y,y_arr(1),y_arr(nY)
    183197         if(y.lt.y_arr(1))then
    184198            y=y_arr(1)+0.01
    … …  
    187201            y=y_arr(nY)-0.01
    188202         endif
    189       endif
    190       !else
     203      else
    191204
    192205!     in the y (temperature) direction 2nd
    … …  
    202215            goto 20
    203216         endif
    204       !endif
     217      endif
    205218     
    206219      f11=f2d_arr(a,b)
    … …  
    208221      f12=f2d_arr(a,b+1)
    209222      f22=f2d_arr(a+1,b+1)
    210      
    211 !     1st in x-direction
    212       fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
    213       fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
    214      
    215 !     then in y-direction
    216       f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
    217      
     223
     224      call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
     225
    218226      return
    219227    end subroutine bilinearH2H2

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

    r590 r716  
    88  implicit none
    99
    10 !==================================================================
    11 !     
    12 !     Purpose
    13 !     -------
    14 !     Run the universal model radiative transfer once in a 1D column.
    15 !     Useful for climate evolution studies etc.
    16 !     
    17 !     It can be compiled with a command like (e.g. for 25 layers):
    18 !                                  "makegcm -p std -d 25 kcm1d"
    19 !     It requires the files "callphys.def", "gases.def"
    20 !     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
    21 !
    22 !     Authors
    23 !     -------
    24 !     R. Wordsworth
    25 !
    26 !==================================================================
     10  !==================================================================
     11  !     
     12  !     Purpose
     13  !     -------
     14  !     Run the universal model radiative transfer once in a 1D column.
     15  !     Useful for climate evolution studies etc.
     16  !     
     17  !     It can be compiled with a command like (e.g. for 25 layers):
     18  !                                  "makegcm -p std -d 25 kcm1d"
     19  !     It requires the files "callphys.def", "gases.def"
     20  !     "traceur.def", and "run.def" with a line "INCLUDEDEF=callphys.def"
     21  !
     22  !     Authors
     23  !     -------
     24  !     R. Wordsworth
     25  !
     26  !==================================================================
    2727
    2828#include "dimensions.h"
    … …  
    3535#include "advtrac.h"
    3636
    37 ! --------------------------------------------------------------
    38 !  Declarations
    39 ! --------------------------------------------------------------
     37  ! --------------------------------------------------------------
     38  !  Declarations
     39  ! --------------------------------------------------------------
    4040
    4141  integer nlayer,nlevel,nq
    … …  
    6161  real fluxabs_sw    ! SW flux absorbed by planet (W/m2)
    6262  real fluxtop_dn    ! incident top of atmosphere SW flux (W/m2)
    63      
    64 ! not used
     63
     64  ! not used
    6565  real reffrad(nlayermx,naerkind)
    6666  real nueffrad(nlayermx,naerkind)
    … …  
    7171  real dTstrat
    7272  real aerosol(nlayermx,naerkind) ! aerosol tau (kg/kg)
    73   real OLR_nu(L_NSPECTI)
    74   real GSR_nu(L_NSPECTV)
     73  real OLR_nu(ngridmx,L_NSPECTI)
     74  real OSR_nu(ngridmx,L_NSPECTV)
    7575  real Eatmtot
    7676
    7777  integer ierr
    78   logical firstcall,lastcall
    79 
    80 ! --------------
    81 ! Initialisation
    82 ! --------------
     78  logical firstcall,lastcall,global1d
     79
     80  ! --------------
     81  ! Initialisation
     82  ! --------------
    8383
    8484  pi=2.E+0*asin(1.E+0)
    … …  
    8989  totcloudfrac  = 0.0
    9090
    91 !  Load parameters from "run.def"
    92 ! -------------------------------
    93 
    94 ! check if 'run1d.def' file is around (otherwise reading parameters
    95 ! from callphys.def via getin() routine won't work.)
    96   open(90,file='run.def',status='old',form='formatted',&
     91
     92  nlayer=nlayermx
     93  nlevel=nlayer+1
     94
     95
     96  !  Load parameters from "run.def"
     97  ! -------------------------------
     98
     99  ! check if 'kcm1d.def' file is around (otherwise reading parameters
     100  ! from callphys.def via getin() routine won't work.)
     101  open(90,file='kcm1d.def',status='old',form='formatted',&
    97102       iostat=ierr)
    98103  if (ierr.ne.0) then
    99      write(*,*) 'Cannot find required file "run.def"'
     104     write(*,*) 'Cannot find required file "kcm1d.def"'
    100105     write(*,*) '  (which should contain some input parameters'
    101106     write(*,*) '   along with the following line:'
    … …  
    108113  endif
    109114
    110   nlayer=nlayermx
    111   nlevel=nlayer+1
     115! now, run.def is needed anyway. so we create a dummy temporary one
     116! for ioipsl to work. if a run.def is already here, stop the
     117! program and ask the user to do a bit of cleaning
     118  open(90,file='run.def',status='old',form='formatted',&
     119       iostat=ierr)
     120  if (ierr.eq.0) then
     121     close(90)
     122     write(*,*) 'There is already a run.def file.'
     123     write(*,*) 'This is not compatible with 1D runs.'
     124     write(*,*) 'Please remove the file and restart the run.'
     125     write(*,*) 'Runtime parameters are supposed to be in kcm1d.def'
     126     stop
     127  else
     128     call system('touch run.def')
     129     call system("echo 'INCLUDEDEF=callphys.def' >> run.def")
     130     call system("echo 'INCLUDEDEF=kcm1d.def' >> run.def")
     131  endif
     132
     133  global1d = .false. ! default value
     134  call getin("global1d",global1d)
     135  if(.not.global1d)then
     136     print*,'Error, kcm1d must have global1d=.true. in kcm1d.def!'
     137     stop
     138  end if
    112139
    113140  psurf_n=100000. ! default value for psurf
    … …  
    115142  call getin("psurf",psurf_n)
    116143  write(*,*) " psurf = ",psurf_n
    117      
     144
     145! OK. now that run.def has been read once -- any variable is in memory.
     146! so we can dump the dummy run.def
     147  call system("rm -rf run.def")
     148
    118149  tsurf=300.0 ! default value for tsurf
    119150  print*,'Surface temperature (K)?'
    … …  
    125156  call getin("g",g)
    126157  write(*,*) " g = ",g
    127  
     158
    128159  periastr = 1.0
    129160  apoastr  = 1.0
    … …  
    136167  call getin("apoastr",apoastr)
    137168  write(*,*) "apoastron = ",apoastr
    138  
     169
    139170  albedo=0.2 ! default value for albedo
    140171  print*,'Albedo of bare ground?'
    … …  
    155186  cpp=0.
    156187
     188  check_cpp_match = .false.
     189  call getin("check_cpp_match",check_cpp_match)
     190  if (check_cpp_match) then
     191     print*,"In 1D modeling, check_cpp_match is supposed to be F"
     192     print*,"Please correct callphys.def"
     193     stop
     194  endif
     195
    157196  call su_gases
    158197  call calc_cpp_mugaz
    159198
    160 
    161199  call inifis(1,llm,0,86400.0,1.0,0.0,0.0,1.0,rad,g,r,cpp)
    162200
    163 ! Tracer initialisation
    164 ! ---------------------
     201  ! Tracer initialisation
     202  ! ---------------------
    165203  if (tracer) then
    166204     ! load tracer names from file 'traceur.def'
    … …  
    183221           endif
    184222        endif
    185        
     223
    186224        do iq=1,nq
    187225           ! minimal version, just read in the tracer names, 1 per line
    … …  
    193231        enddo !of do iq=1,nq
    194232     endif
    195   !endif
    196  
    197   call initracer()
     233
     234     call initracer()
    198235
    199236  endif
    … …  
    205242     enddo
    206243  enddo
    207  
     244
    208245  do iq=1,nqmx
    209246     qsurf(iq) = 0.
    … …  
    214251
    215252  iter    = 1
    216   Tstrat  = 60.0
     253  Tstrat  = 150.0
    217254  dTstrat = 1000.0
    218255
    219 ! ---------
    220 ! Run model
    221 ! ---------
    222   do
     256  ! ---------
     257  ! Run model
     258  ! ---------
     259  !do
    223260     psurf = psurf_n
    224261
    225 !    Create vertical profiles
     262     !    Create vertical profiles
    226263     call kcmprof_fn(psurf,qsurf(1),tsurf,    &
    227264          tstrat,play,plev,zlay,temp,q(:,1),muvar(1))
    228265
    229 !    Run radiative transfer
     266     !    Run radiative transfer
    230267     call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
    231268          albedo,emis,mu0,plev,play,temp,                    &
    232           tsurf,fract,dist_star,aerosol,cpp,muvar,         &
     269          tsurf,fract,dist_star,aerosol,muvar,         &
    233270          dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
    234           fluxabs_sw,fluxtop_dn,reffrad,tau_col,  &
     271          fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,reffrad,tau_col,  &
    235272          cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
    236  
    237 !    Iterate stratospheric temperature
     273
     274     !    Iterate stratospheric temperature
    238275     print*,'Tstrat = ',Tstrat
    239276     dTstrat = Tstrat
    … …  
    244281     dTstrat = dTstrat-Tstrat
    245282
    246      if(abs(dTstrat).lt.1.0)then
    247         print*,'dTstrat = ',dTstrat
    248         exit
    249      endif
    250 
    251      iter=iter+1
    252      if(iter.eq.100)then
    253         print*,'Stratosphere failed to converge after'
    254         print*,'100 iteration, aborting run.'
    255         call abort
    256      endif
    257 
    258   end do
    259 
    260 ! Calculate total atmospheric energy
     283     !if(abs(dTstrat).lt.1.0)then
     284     !   print*,'dTstrat = ',dTstrat
     285     !   exit
     286     !endif
     287
     288     !iter=iter+1
     289     !if(iter.eq.100)then
     290     !   print*,'Stratosphere failed to converge after'
     291     !   print*,'100 iteration, aborting run.'
     292     !   call abort
     293     !endif
     294
     295  !end do
     296
     297  ! Run radiative transfer one last time to get OLR,OSR
     298  firstcall=.false.
     299  lastcall=.true.
     300  call callcorrk(ngridmx,nlayer,q,nqmx,qsurf,      &
     301       albedo,emis,mu0,plev,play,temp,                    &
     302       tsurf,fract,dist_star,aerosol,muvar,         &
     303       dtlw,dtsw,fluxsurf_lw,fluxsurf_sw,fluxtop_lw,    &
     304       fluxabs_sw,fluxtop_dn,OLR_nu,OSR_nu,             &
     305       reffrad,tau_col,  &
     306       cloudfrac,totcloudfrac,.false.,firstcall,lastcall)
     307
     308
     309  ! Calculate total atmospheric energy
    261310  Eatmtot=0.0
    262 !  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
    263 !     q(:,1),muvar,Eatmtot)
    264 
    265 ! ------------------------
    266 ! Save data to ascii files
    267 ! ------------------------
     311  !  call calcenergy_kcm(tsurf,temp,play,plev,qsurf,&
     312  !     q(:,1),muvar,Eatmtot)
     313
     314  ! ------------------------
     315  ! Save data to ascii files
     316  ! ------------------------
    268317
    269318  print*,'Saving profiles...'
    … …  
    293342  close(119)
    294343
     344  print*, tsurf,psurf,fluxtop_dn,fluxabs_sw,fluxtop_lw
     345
    295346  print*,'Saving scalars...'
    296347  open(116,file='surf_vals.out')
    … …  
    305356  open(117,file='OLRnu.out')
    306357  do iw=1,L_NSPECTI
    307      write(117,*) OLR_nu(iw)
     358     write(117,*) OLR_nu(1,iw)
    308359  enddo
    309360  close(117)
    310  
    311   open(127,file='GSRnu.out')
     361
     362  open(127,file='OSRnu.out')
    312363  do iw=1,L_NSPECTV
    313      write(127,*) GSR_nu(iw)
     364     write(127,*) OSR_nu(1,iw)
    314365  enddo
    315366  close(127) 

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

    r471 r716  
    4040  double precision Dz, Dp
    4141  double precision Ptop, dlogp, Psat_max
    42   parameter (Ptop=1.0)           ! Pressure at TOA [Pa]
    43   !parameter (Psat_max=100000.0) ! Maximum vapour pressure [Pa]
    44   parameter (Psat_max=0.0) ! set to zero for dry calculations
     42  parameter (Ptop=1.0)                             ! Pressure at TOA [Pa]
    4543
    4644  double precision T(1:nlay)                       ! temperature [K]
    … …  
    6462  double precision psat_v      ! local Psat_H2O value
    6563  double precision Tcrit       ! Critical temperature [K]
    66 !  double precision Tsat        ! local Tsat_H2O value
    6764  double precision rho_vTEMP,rho_nTEMP
     65
     66  double precision TCO2cond ! for CO2 condensation quasi-hack
    6867
    6968  ! variables necessary for steam.f90
    … …  
    7473
    7574  logical verbose
    76   parameter(verbose=.false.)
     75  parameter(verbose=.true.)
     76
     77  logical add_Pvar_to_total
     78  parameter(add_Pvar_to_total=.true.)
    7779
    7880  ! initialise flags
    … …  
    8284  ! assign input variables
    8385  m_n     = dble(mugaz/1000.)
     86  cp_n    = cpp
    8487  ! modify/generalise later??
    8588
    86   if(ngasmx.gt.2)then
    87      print*,'Only two species possible at present'
    88      stop
    89   elseif(ngasmx.eq.1)then
     89  Psat_max = 1000000.0 ! maximum vapour pressure [Pa]
     90                       ! set huge until further notice
     91
     92  if(vgas.lt.1)then
    9093     if(psat_max.gt.0.0)then
    9194        print*,'Must have Psat_max=0 if no variable species'
    92         stop
     95        psat_max=0.0
     96        !stop
    9397     endif
    9498     print*, 'Assuming pure atmosphere'
    9599     m_v   = 1.0
    96100     tcrit = 1000.0
    97   elseif(gnom(ngasmx).eq.'H2O')then
     101  elseif(gnom(vgas).eq.'H2O')then
    98102     m_v   = dble(mH2O/1000.)
    99103     tcrit = 6.47d2
    100   elseif(gnom(ngasmx).eq.'NH3')then
     104  elseif(gnom(vgas).eq.'NH3')then
    101105     m_v   = 17.031/1000.
    102106     tcrit = 4.06d2
    103   elseif(gnom(ngasmx).eq.'CH4')then
     107  elseif(gnom(vgas).eq.'CH4')then
    104108     m_v   = 16.04/1000.
    105109     tcrit = 1.91d2
    … …  
    111115
    112116  rmn     = rc/m_n
    113  
    114   Psurf_n = dble(psurf_rcm)
    115117  Ttop    = dble(Tstra_rcm)
    116118  Tsurf   = dble(Tsurf_rcm)
    117119
    118 
    119   ! assume vapour saturation at surface (for now)
    120 !  if (tsurf > tcrit) then
    121 !     print*,'Above critical temperature for ',gnom(2),&
    122 !          ', at surface, assuming 10 bar pressure instead'
    123 !     Psurf_v = 10*1d5
    124 !     profil_flag(1) = 0
    125 !i!  endif
    126 
    127 
    128   psat_v=psat_max
    129   if(gnom(ngasmx).eq.'H2O')then
    130      call Psat_H2O(tsurf,psat_v)
    131   elseif(gnom(ngasmx).eq.'NH3')then
    132      call Psat_NH3(tsurf,psat_v)
     120  psat_v  = psat_max
     121  if(vgas.gt.0)then
     122     if(gnom(vgas).eq.'H2O')then
     123        call Psat_H2O(tsurf,psat_v)
     124     elseif(gnom(vgas).eq.'NH3')then
     125        call Psat_NH3(tsurf,psat_v)
     126     endif
    133127  endif
    134128
    … …  
    142136  endif
    143137
    144 
     138  if(add_Pvar_to_total)then
     139    Psurf_n = dble(psurf_rcm)
     140    psurf_rcm = real(Psurf_n+Psurf_v)
     141  else
     142    Psurf_n = dble(psurf_rcm) - Psurf_v
     143  endif
     144
     145  ! include relative humidity option
     146  !if(satval.lt.1.0)then
     147  !   Psurf_v = Psurf_v*satval
     148  !   profil_flag(1) = 0     
     149  !endif
    145150
    146151  if(verbose)then
    … …  
    155160  endif
    156161
    157 
    158 
    159162  ! define fine pressure grid
    160   psurf_rcm = real(Psurf_n+Psurf_v)
    161163  dlogp_rcm = -(log(psurf_rcm)-log(ptop))/nlayermx
    162164
    … …  
    207209  do ilay=2,nlay-1
    208210
    209 
    210 
    211 
    212211     ! calculate altitude levels (for diagnostic only)
    213212     Dz         = -Dp/(  g*(rho_n(ilay) + rho_v(ilay))  )
    … …  
    228227     ! test for moist adiabat at next level
    229228     psat_v=psat_max
    230      if(gnom(ngasmx).eq.'H2O')then
     229
     230     if(vgas.gt.0)then
     231     if(gnom(vgas).eq.'H2O')then
    231232        call Psat_H2O(T(ilay+1),psat_v)
    232      elseif(gnom(ngasmx).eq.'NH3')then
     233     elseif(gnom(vgas).eq.'NH3')then
    233234        call Psat_NH3(T(ilay+1),psat_v)
     235     endif
    234236     endif
    235237
    … …  
    253255       
    254256        psat_v=psat_max
    255         if(gnom(ngasmx).eq.'H2O')then
     257
     258        if(vgas.gt.0)then
     259        if(gnom(vgas).eq.'H2O')then
    256260           call Psat_H2O(T(ilay+1),psat_v)
    257         elseif(gnom(ngasmx).eq.'NH3')then
     261        elseif(gnom(vgas).eq.'NH3')then
    258262           call Psat_NH3(T(ilay+1),psat_v)
     263        endif
    259264        endif
    260265
    … …  
    281286     end select
    282287
    283 
    284 
    285 !         print*,'profil_flag=',profil_flag(ilay)
    286 !         print*,'T=',T(ilay)
    287 !         print*,'a=',rho_v(ilay)/rho_n(ilay)
    288 !      if(profil_flag(ilay).eq.1)then
    289 !           stop
    290 !     endif
    291 
    292 
    293 
    294 
    295 
    296288  enddo
    297289
    … …  
    378370  enddo
    379371
     372!    CO2 condensation 'haircut' of temperature profile if necessary
     373  if(co2cond)then
     374     print*,'CO2 condensation haircut - assumes CO2-dominated atmosphere!'
     375     do ilay=2,nlayermx
     376        if(P_rcm(ilay).lt.518000.)then
     377           TCO2cond = (-3167.8)/(log(.01*P_rcm(ilay))-23.23) ! Fanale's formula
     378        else
     379           TCO2cond = 684.2-92.3*log(P_rcm(ilay))+4.32*log(P_rcm(ilay))**2
     380           ! liquid-vapour transition (based on CRC handbook 2003 data)
     381        endif
     382
     383        print*,'p=',P_rcm(ilay),', T=',T_rcm(ilay),' Tcond=',TCO2cond
     384        if(T_rcm(ilay).lt.TCO2cond)then
     385           T_rcm(ilay)=TCO2cond
     386        endif
     387     enddo
     388  endif
     389
    380390  return
    381391end subroutine kcmprof_fn

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

    r596 r716  
    1       subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
    2            QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
    3            TMID,PMID,TAUGSURF,QVAR,MUVAR)
    4 
    5       use radinc_h
    6       use radcommon_h, only: gasi, tlimit, wrefVAR, Cmk,tgasref,pfgasref,wnoi,scalep
    7       use gases_h
    8       implicit none
    9 
    10 !==================================================================
    11 !     
    12 !     Purpose
    13 !     -------
    14 !     Calculates longwave optical constants at each level. For each
    15 !     layer and spectral interval in the IR it calculates WBAR, DTAU
    16 !     and COSBAR. For each level it calculates TAU.
    17 !     
    18 !     TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
    19 !     at the *bottom* of layer L), LW is the spectral wavelength interval.
    20 !     
    21 !     TLEV(L) - Temperature at the layer boundary (i.e., level)
    22 !     PLEV(L) - Pressure at the layer boundary (i.e., level)
    23 !
    24 !     Authors
    25 !     -------
    26 !     Adapted from the NASA Ames code by R. Wordsworth (2009)
    27 !     
    28 !==================================================================
    29 
     1subroutine optci(PLEV,TLEV,DTAUI,TAUCUMI,      &
     2     QXIAER,QSIAER,GIAER,COSBI,WBARI,TAUAERO,  &
     3     TMID,PMID,TAUGSURF,QVAR,MUVAR)
     4
     5  use radinc_h
     6  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep
     7  use gases_h
     8  implicit none
     9
     10  !==================================================================
     11  !     
     12  !     Purpose
     13  !     -------
     14  !     Calculates longwave optical constants at each level. For each
     15  !     layer and spectral interval in the IR it calculates WBAR, DTAU
     16  !     and COSBAR. For each level it calculates TAU.
     17  !     
     18  !     TAUI(L,LW) is the cumulative optical depth at level L (or alternatively
     19  !     at the *bottom* of layer L), LW is the spectral wavelength interval.
     20  !     
     21  !     TLEV(L) - Temperature at the layer boundary (i.e., level)
     22  !     PLEV(L) - Pressure at the layer boundary (i.e., level)
     23  !
     24  !     Authors
     25  !     -------
     26  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
     27  !     
     28  !==================================================================
    3029
    3130#include "comcstfi.h"
    … …  
    3332
    3433
    35       real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    36       real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
    37       real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
    38       real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
    39       real*8 PLEV(L_LEVELS)
    40       real*8 TLEV(L_LEVELS)
    41       real*8 TMID(L_LEVELS), PMID(L_LEVELS)
    42       real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    43       real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
    44 
    45 !     For aerosols
    46       real*8  QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
    47       real*8  QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
    48       real*8  GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
    49       real*8  TAUAERO(L_LEVELS+1,NAERKIND)
    50       real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
    51       real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
    52 
    53       integer L, NW, NG, K, LK, IAER
    54       integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
    55       real*8  ANS, TAUGAS
    56       real*8  DPR(L_LEVELS), U(L_LEVELS)
    57       real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
    58 
    59       real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
    60       real*8 DCONT
    61       double precision wn_cont, p_cont, p_air, T_cont, dtemp
    62 
    63 !     Variable species mixing ratio variables
    64       real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    65       real*8  KCOEF(4)
    66       integer NVAR(L_LEVELS)
    67 
    68 !     temporary variables for multiple aerosol calculation
    69       real*8 atemp, btemp
    70 
    71 !     variables for k in units m^-1
    72       real*8 rho, dz(L_LEVELS)
    73 
    74       integer igas
    75 
    76       !--- Kasting's CIA ----------------------------------------
    77       !real*8, parameter :: Ci(L_NSPECTI)=[                         &
    78       !     3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7,    &
    79       !     5.4E-7, 1.6E-6, 0.0,                               &
    80       !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,            &
    81       !     0.0, 4.0E-7, 4.0E-6, 1.4E-5,    &
    82       !     1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
    83       !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9,             &
    84       !     -1.7, -1.7, -1.7, -1.7, -1.7, -1.7,                &
    85       !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
    86       !     -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
    87       !----------------------------------------------------------
    88 
    89 !=======================================================================
    90 !     Determine the total gas opacity throughout the column, for each
    91 !     spectral interval, NW, and each Gauss point, NG.
    92 
    93       taugsurf(:,:) = 0.0
    94       dpr(:)        = 0.0
    95       lkcoef(:,:)   = 0.0
    96 
    97       do K=2,L_LEVELS
    98          DPR(k) = PLEV(K)-PLEV(K-1)
    99 
    100          !--- Kasting's CIA ----------------------------------------
    101          !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K))
    102          ! this is CO2 path length (in cm) as written by Francois
    103          ! delta_z = delta_p * R_specific * T / (g * P)
    104          ! But Kasting states that W is in units of _atmosphere_ cm
    105          ! So we do
    106          !dz(k)=dz(k)*(PMID(K)/1013.25)
    107          !dz(k)=dz(k)/100.0 ! in m for SI calc
    108          !----------------------------------------------------------
    109 
    110          ! if we have continuum opacities, we need dz
    111          if(kastprof)then
    112             dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
    113             U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
    114          else
    115             dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
    116             U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
    117          endif
    118 
    119          call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
    120               LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
    121          
    122          do LK=1,4
    123             LKCOEF(K,LK) = LCOEF(LK)
    124          end do
    125 
    126 
    127          DO NW=1,L_NSPECTI
    128             do iaer=1,naerkind
    129                TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
    130             end do
    131          END DO
    132       end do                    ! levels
    133 
    134       do K=2,L_LEVELS
    135          do nw=1,L_NSPECTI
    136  
    137             DCONT = 0.0 ! continuum absorption
    138 
    139             if(Continuum.and.(.not.graybody))then
    140             ! include continua if necessary
    141                wn_cont = dble(wnoi(nw))
    142                T_cont  = dble(TMID(k))
    143                do igas=1,ngasmx
    144 
    145                   if(gfrac(igas).eq.-1)then ! variable
    146                      p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
    147                   else
    148                      p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
    149                   endif
    150 
    151                   dtemp=0.0
    152                   if(gnom(igas).eq.'H2_')then
    153                      call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
    154                   elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
    155                      p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
    156                      call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
    157                   endif
    158 
    159                   DCONT = DCONT + dtemp
    160                enddo
    161 
    162                DCONT = DCONT*dz(k)
    163             endif
    164 
    165             !--- Kasting's CIA ----------------------------------------
    166             !DCO2   = dz(k)*Ci(nw)*(1.2859*PMID(k)/1000.0)*(TMID(k)/300.)**Ti(nw)
    167             !DCO2 = 130*Ci(nw)*(pmid(k)/1013.25)**2*(tmid(k)/300.)**Ti(nw) * dz(k)
    168             ! these two have been verified to give the same results
    169             !----------------------------------------------------------
    170 
    171 
    172             do ng=1,L_NGAUSS-1
    173 
    174 !     Now compute TAUGAS
    175 
    176 !     Interpolate between water mixing ratios
    177 !     WRATIO = 0.0 if the requested water amount is equal to, or outside the
    178 !     the water data range
    179 
    180                if(L_REFVAR.eq.1)then ! added by RW for special no variable case
    181                   KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
    182                   KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
    183                   KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
    184                   KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
    185                else
    186 
    187                   KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*     &
    188                        (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   &
    189                        GASI(MT(K),MP(K),NVAR(K),NW,NG))
    190 
    191                   KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*   &
    192                        (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                 &
    193                        GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
    194 
    195                   KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
    196                        (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -               &
    197                        GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
    198 
    199                   KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*   &
    200                        (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
    201                        GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
    202                endif
    203 
    204 !     Interpolate the gaseous k-coefficients to the requested T,P values
    205 
    206                ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
    207                     LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
    208                
    209                TAUGAS  = U(k)*ANS
    210 
    211                TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    212                !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
    213                DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
    214 
    215                do iaer=1,naerkind
    216                   DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
    217                end do ! a bug was here!
    218 
    219             end do
    220 
    221 !     Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
    222 !     which holds continuum opacity only
    223 
    224             NG              = L_NGAUSS
    225             DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
    226 
    227             do iaer=1,naerkind
    228                DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) +  TAEROS(K,NW,IAER)
    229             end do ! a bug was here!
    230 
    231          end do
    232       end do
    233 
    234 
    235 !=======================================================================
    236 !     Now the full treatment for the layers, where besides the opacity
    237 !     we need to calculate the scattering albedo and asymmetry factors
    238 
    239       DO NW=1,L_NSPECTI
    240          DO K=2,L_LEVELS+1
    241             do iaer=1,naerkind
    242                TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
    243             end do
    244          ENDDO
    245       ENDDO
    246 
    247       DO NW=1,L_NSPECTI
    248          NG = L_NGAUSS
    249          DO L=1,L_NLAYRAD
    250 
    251             K              = 2*L+1
    252             DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
    253 
    254             atemp = 0.
    255             btemp = 0.
    256             if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
    257                do iaer=1,naerkind
    258                   atemp = atemp +                                     &
    259                       GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
    260                       GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
    261                   btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
    262 !     *                    + 1.e-10
    263                end do
    264                WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
    265             else
    266                WBARI(L,nw,ng) = 0.0D0
    267                DTAUI(L,NW,NG) = 1.0E-9
    268             endif
    269 
    270             if(btemp .GT. 0.0) then
    271                cosbi(L,NW,NG) = atemp/btemp
    272             else
    273                cosbi(L,NW,NG) = 0.0D0
    274             end if
    275 
    276          END DO ! L vertical loop
    277 
    278 !     Now the other Gauss points, if needed.
    279 
    280          DO NG=1,L_NGAUSS-1
    281             IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
    282 
    283                DO L=1,L_NLAYRAD
    284                   K              = 2*L+1
    285                   DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
    286 
    287                   btemp = 0.
    288                   if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
    289 
    290                      do iaer=1,naerkind
    291                         btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
    292                      end do
    293                      WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
    294 
    295                   else
    296                      WBARI(L,nw,ng) = 0.0D0
    297                      DTAUI(L,NW,NG) = 1.0E-9
    298                   endif
    299 
    300                   cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
    301                END DO ! L vertical loop
    302             END IF
    303            
    304          END DO                 ! NG Gauss loop
    305       END DO                    ! NW spectral loop
    306 
    307 !     Total extinction optical depths
    308 
    309       DO NW=1,L_NSPECTI       
    310          DO NG=1,L_NGAUSS       ! full gauss loop
    311             TAUI(1,NW,NG)=0.0D0
    312             DO L=1,L_NLAYRAD
    313                TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
    314             END DO
    315 
    316             TAUCUMI(1,NW,NG)=0.0D0
    317             DO K=2,L_LEVELS
    318                TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
    319             END DO
    320          END DO                 ! end full gauss loop
    321       END DO
    322 
    323       return
    324 
    325 
    326     end subroutine optci
    327 
    328 
    329 
     34  real*8 DTAUI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     35  real*8 DTAUKI(L_LEVELS+1,L_NSPECTI,L_NGAUSS)
     36  real*8 TAUI(L_NLEVRAD,L_NSPECTI,L_NGAUSS)
     37  real*8 TAUCUMI(L_LEVELS,L_NSPECTI,L_NGAUSS)
     38  real*8 PLEV(L_LEVELS)
     39  real*8 TLEV(L_LEVELS)
     40  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
     41  real*8 COSBI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     42  real*8 WBARI(L_NLAYRAD,L_NSPECTI,L_NGAUSS)
     43
     44  ! for aerosols
     45  real*8  QXIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
     46  real*8  QSIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
     47  real*8  GIAER(L_LEVELS+1,L_NSPECTI,NAERKIND)
     48  real*8  TAUAERO(L_LEVELS+1,NAERKIND)
     49  real*8  TAUAEROLK(L_LEVELS+1,L_NSPECTI,NAERKIND)
     50  real*8  TAEROS(L_LEVELS,L_NSPECTI,NAERKIND)
     51
     52  integer L, NW, NG, K, LK, IAER
     53  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
     54  real*8  ANS, TAUGAS
     55  real*8  DPR(L_LEVELS), U(L_LEVELS)
     56  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
     57
     58  real*8 taugsurf(L_NSPECTI,L_NGAUSS-1)
     59  real*8 DCONT
     60  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
     61  double precision p_cross
     62
     63  ! variable species mixing ratio variables
     64  real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
     65  real*8  KCOEF(4)
     66  integer NVAR(L_LEVELS)
     67
     68  ! temporary variables for multiple aerosol calculation
     69  real*8 atemp, btemp
     70
     71  ! variables for k in units m^-1
     72  real*8 rho, dz(L_LEVELS)
     73
     74  integer igas, jgas
     75
     76  !--- Kasting's CIA ----------------------------------------
     77  !real*8, parameter :: Ci(L_NSPECTI)=[                         &
     78  !     3.8E-5, 1.2E-5, 2.8E-6, 7.6E-7, 4.5E-7, 2.3E-7,    &
     79  !     5.4E-7, 1.6E-6, 0.0,                               &
     80  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,            &
     81  !     0.0, 4.0E-7, 4.0E-6, 1.4E-5,    &
     82  !     1.0E-5, 1.2E-6, 2.0E-7, 5.0E-8, 3.0E-8, 0.0 ]
     83  !real*8, parameter :: Ti(L_NSPECTI)=[ -2.2, -1.9,             &
     84  !     -1.7, -1.7, -1.7, -1.7, -1.7, -1.7,                &
     85  !     0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, &
     86  !     -1.7,-1.7,-1.7,-1.7,-1.7,-1.7,-1.7, -1.7,0.0 ]
     87  !----------------------------------------------------------
     88
     89  !=======================================================================
     90  !     Determine the total gas opacity throughout the column, for each
     91  !     spectral interval, NW, and each Gauss point, NG.
     92
     93  taugsurf(:,:) = 0.0
     94  dpr(:)        = 0.0
     95  lkcoef(:,:)   = 0.0
     96
     97  do K=2,L_LEVELS
     98     DPR(k) = PLEV(K)-PLEV(K-1)
     99
     100     !--- Kasting's CIA ----------------------------------------
     101     !dz(k)=dpr(k)*189.02*TMID(K)/(0.03720*PMID(K))
     102     ! this is CO2 path length (in cm) as written by Francois
     103     ! delta_z = delta_p * R_specific * T / (g * P)
     104     ! But Kasting states that W is in units of _atmosphere_ cm
     105     ! So we do
     106     !dz(k)=dz(k)*(PMID(K)/1013.25)
     107     !dz(k)=dz(k)/100.0 ! in m for SI calc
     108     !----------------------------------------------------------
     109
     110     ! if we have continuum opacities, we need dz
     111     if(kastprof)then
     112        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
     113        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
     114     else
     115        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
     116        U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
     117     endif
     118
     119     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
     120          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
     121
     122     do LK=1,4
     123        LKCOEF(K,LK) = LCOEF(LK)
     124     end do
     125
     126
     127     DO NW=1,L_NSPECTI
     128        do iaer=1,naerkind
     129           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXIAER(K,NW,IAER)
     130        end do
     131     END DO
     132  end do                    ! levels
     133
     134  do K=2,L_LEVELS
     135     do nw=1,L_NSPECTI
     136
     137        DCONT = 0.0 ! continuum absorption
     138
     139        if(Continuum.and.(.not.graybody))then
     140           ! include continua if necessary
     141           wn_cont = dble(wnoi(nw))
     142           T_cont  = dble(TMID(k))
     143           do igas=1,ngasmx
     144
     145              if(gfrac(igas).eq.-1)then ! variable
     146                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
     147              else
     148                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
     149              endif
     150
     151              dtemp=0.0
     152              if(igas.eq.igas_N2)then
     153
     154                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
     155
     156              elseif(igas.eq.igas_H2)then
     157
     158                 ! first do self-induced absorption
     159                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
     160
     161                 ! then cross-interactions with other gases
     162                 do jgas=1,ngasmx
     163                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
     164                    dtempc  = 0.0
     165                    if(jgas.eq.igas_N2)then
     166                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
     167                    elseif(jgas.eq.igas_He)then
     168                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
     169                    endif
     170                    dtemp = dtemp + dtempc
     171                 enddo
     172
     173              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
     174
     175                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
     176                 if(H2Ocont_simple)then
     177                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     178                 else
     179                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     180                 endif
     181
     182              endif
     183
     184              DCONT = DCONT + dtemp
     185
     186           enddo
     187
     188           ! Oobleck test
     189           !rho = PMID(k)*scalep / (TMID(k)*286.99)
     190           !if(WNOI(nw).gt.300.0 .and. WNOI(nw).lt.500.0)then
     191           !   DCONT = rho * 0.125 * 4.6e-4
     192           !elseif(WNOI(nw).gt.500.0 .and. WNOI(nw).lt.700.0)then
     193           !   DCONT = 1000*dpr(k) * 1.0 * 4.6e-4 / g
     194           !   DCONT = rho * 1.0 * 4.6e-4
     195           !elseif(WNOI(nw).gt.700.0 .and. WNOI(nw).lt.900.0)then
     196           !   DCONT = rho * 0.125 * 4.6e-4
     197           !endif
     198
     199           DCONT = DCONT*dz(k)
     200
     201        endif
     202
     203        ! RW 7/3/12: already done above
     204        !if(.not.Continuum)then
     205        !   DCONT=0.0
     206        !endif
     207
     208        !--- Kasting's CIA ----------------------------------------
     209        !DCO2   = dz(k)*Ci(nw)*(1.2859*PMID(k)/1000.0)*(TMID(k)/300.)**Ti(nw)
     210        !DCO2 = 130*Ci(nw)*(pmid(k)/1013.25)**2*(tmid(k)/300.)**Ti(nw) * dz(k)
     211        ! these two have been verified to give the same results
     212        !----------------------------------------------------------
     213
     214        do ng=1,L_NGAUSS-1
     215
     216           ! Now compute TAUGAS
     217
     218           ! Interpolate between water mixing ratios
     219           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
     220           ! the water data range
     221
     222           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
     223              KCOEF(1) = GASI(MT(K),MP(K),1,NW,NG)
     224              KCOEF(2) = GASI(MT(K),MP(K)+1,1,NW,NG)
     225              KCOEF(3) = GASI(MT(K)+1,MP(K)+1,1,NW,NG)
     226              KCOEF(4) = GASI(MT(K)+1,MP(K),1,NW,NG)
     227           else
     228
     229              KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*     &
     230                   (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   &
     231                   GASI(MT(K),MP(K),NVAR(K),NW,NG))
     232
     233              KCOEF(2) = GASI(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*   &
     234                   (GASI(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                 &
     235                   GASI(MT(K),MP(K)+1,NVAR(K),NW,NG))
     236
     237              KCOEF(3) = GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)* &
     238                   (GASI(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -               &
     239                   GASI(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
     240
     241              KCOEF(4) = GASI(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*   &
     242                   (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
     243                   GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
     244           endif
     245
     246           ! Interpolate the gaseous k-coefficients to the requested T,P values
     247
     248           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
     249                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
     250
     251           TAUGAS  = U(k)*ANS
     252
     253           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
     254           DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
     255
     256           do iaer=1,naerkind
     257              DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) + TAEROS(K,NW,IAER)
     258           end do ! a bug was here!
     259
     260        end do
     261
     262        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
     263        ! which holds continuum opacity only
     264
     265        NG              = L_NGAUSS
     266        DTAUKI(K,nw,ng) = 0.0 + DCONT ! For parameterized continuum absorption
     267
     268        do iaer=1,naerkind
     269           DTAUKI(K,nw,ng) = DTAUKI(K,nw,ng) +  TAEROS(K,NW,IAER)
     270        end do ! a bug was here!
     271
     272     end do
     273  end do
     274
     275
     276  !=======================================================================
     277  !     Now the full treatment for the layers, where besides the opacity
     278  !     we need to calculate the scattering albedo and asymmetry factors
     279
     280  DO NW=1,L_NSPECTI
     281     DO K=2,L_LEVELS+1
     282        do iaer=1,naerkind
     283           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
     284        end do
     285     ENDDO
     286  ENDDO
     287
     288  DO NW=1,L_NSPECTI
     289     NG = L_NGAUSS
     290     DO L=1,L_NLAYRAD
     291
     292        K              = 2*L+1
     293        DTAUI(L,nw,ng) = DTAUKI(K,NW,NG) + DTAUKI(K+1,NW,NG)! + 1.e-50
     294
     295        atemp = 0.
     296        btemp = 0.
     297        if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
     298           do iaer=1,naerkind
     299              atemp = atemp +                                     &
     300                   GIAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
     301                   GIAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
     302              btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
     303              !     *                    + 1.e-10
     304           end do
     305           WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
     306        else
     307           WBARI(L,nw,ng) = 0.0D0
     308           DTAUI(L,NW,NG) = 1.0E-9
     309        endif
     310
     311        if(btemp .GT. 0.0) then
     312           cosbi(L,NW,NG) = atemp/btemp
     313        else
     314           cosbi(L,NW,NG) = 0.0D0
     315        end if
     316
     317     END DO ! L vertical loop
     318
     319     ! Now the other Gauss points, if needed.
     320
     321     DO NG=1,L_NGAUSS-1
     322        IF(TAUGSURF(NW,NG) .gt. TLIMIT) THEN
     323
     324           DO L=1,L_NLAYRAD
     325              K              = 2*L+1
     326              DTAUI(L,nw,ng) = DTAUKI(K,NW,NG)+DTAUKI(K+1,NW,NG)! + 1.e-50
     327
     328              btemp = 0.
     329              if(DTAUI(L,NW,NG) .GT. 1.0E-9) then
     330
     331                 do iaer=1,naerkind
     332                    btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
     333                 end do
     334                 WBARI(L,nw,ng) = btemp  / DTAUI(L,NW,NG)
     335
     336              else
     337                 WBARI(L,nw,ng) = 0.0D0
     338                 DTAUI(L,NW,NG) = 1.0E-9
     339              endif
     340
     341              cosbi(L,NW,NG) = cosbi(L,NW,L_NGAUSS)
     342           END DO ! L vertical loop
     343        END IF
     344
     345     END DO                 ! NG Gauss loop
     346  END DO                    ! NW spectral loop
     347
     348  ! Total extinction optical depths
     349
     350  DO NW=1,L_NSPECTI       
     351     DO NG=1,L_NGAUSS       ! full gauss loop
     352        TAUI(1,NW,NG)=0.0D0
     353        DO L=1,L_NLAYRAD
     354           TAUI(L+1,NW,NG)=TAUI(L,NW,NG)+DTAUI(L,NW,NG)
     355        END DO
     356
     357        TAUCUMI(1,NW,NG)=0.0D0
     358        DO K=2,L_LEVELS
     359           TAUCUMI(K,NW,NG)=TAUCUMI(K-1,NW,NG)+DTAUKI(K,NW,NG)
     360        END DO
     361     END DO                 ! end full gauss loop
     362  END DO
     363
     364  ! be aware when comparing with textbook results
     365  ! (e.g. Pierrehumbert p. 218) that
     366  ! taucumi does not take the <cos theta>=0.5 factor into
     367  ! account. It is the optical depth for a vertically
     368  ! ascending ray with angle theta = 0.
     369
     370  !open(127,file='taucum.out')
     371  !do nw=1,L_NSPECTI
     372  !   write(127,*) taucumi(L_LEVELS,nw,L_NGAUSS)
     373  !enddo
     374  !close(127)
     375
     376  return
     377
     378
     379end subroutine optci
     380
     381
     382

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

    r596 r716  
    1       SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
    2           QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
    3           TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
    4 
    5       use radinc_h
    6       use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
    7       use gases_h
    8 
    9       implicit none
    10 
    11 !==================================================================
    12 !     
    13 !     Purpose
    14 !     -------
    15 !     Calculates shortwave optical constants at each level.
    16 !     
    17 !     Authors
    18 !     -------
    19 !     Adapted from the NASA Ames code by R. Wordsworth (2009)
    20 !     
    21 !==================================================================
    22 !     
    23 !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL 
    24 !     IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
    25 !     LAYER: WBAR, DTAU, COSBAR
    26 !     LEVEL: TAU
    27 !     
    28 !     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
    29 !     layer L. NW is spectral wavelength interval, ng the Gauss point index.
    30 !     
    31 !     TLEV(L) - Temperature at the layer boundary
    32 !     PLEV(L) - Pressure at the layer boundary (i.e. level)
    33 !     GASV(NT,NPS,NW,NG) - Visible k-coefficients
    34 !     
    35 !-------------------------------------------------------------------
     1SUBROUTINE OPTCV(DTAUV,TAUV,TAUCUMV,PLEV,  &
     2     QXVAER,QSVAER,GVAER,WBARV,COSBV,       &
     3     TAURAY,TAUAERO,TMID,PMID,TAUGSURF,QVAR,MUVAR)
     4
     5  use radinc_h
     6  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
     7  use gases_h
     8
     9  implicit none
     10
     11  !==================================================================
     12  !     
     13  !     Purpose
     14  !     -------
     15  !     Calculates shortwave optical constants at each level.
     16  !     
     17  !     Authors
     18  !     -------
     19  !     Adapted from the NASA Ames code by R. Wordsworth (2009)
     20  !     
     21  !==================================================================
     22  !     
     23  !     THIS SUBROUTINE SETS THE OPTICAL CONSTANTS IN THE VISUAL 
     24  !     IT CALCUALTES FOR EACH LAYER, FOR EACH SPECRAL INTERVAL IN THE VISUAL
     25  !     LAYER: WBAR, DTAU, COSBAR
     26  !     LEVEL: TAU
     27  !     
     28  !     TAUV(L,NW,NG) is the cumulative optical depth at the top of radiation code
     29  !     layer L. NW is spectral wavelength interval, ng the Gauss point index.
     30  !     
     31  !     TLEV(L) - Temperature at the layer boundary
     32  !     PLEV(L) - Pressure at the layer boundary (i.e. level)
     33  !     GASV(NT,NPS,NW,NG) - Visible k-coefficients
     34  !     
     35  !-------------------------------------------------------------------
    3636
    3737#include "callkeys.h"
    3838#include "comcstfi.h"
    3939
    40       real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    41       real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
    42       real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
    43       real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
    44       real*8 PLEV(L_LEVELS)
    45       real*8 TMID(L_LEVELS), PMID(L_LEVELS)
    46       real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    47       real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    48       real*8 TAURAY(L_NSPECTV)
    49 
    50 !     For aerosols
    51       real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
    52       real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
    53       real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
    54       real*8 TAUAERO(L_LEVELS+1,NAERKIND)
    55       real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
    56       real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
    57 
    58       integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER
    59       integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
    60       real*8  ANS, TAUGAS
    61       real*8  TRAY(L_LEVELS,L_NSPECTV)
    62       real*8  DPR(L_LEVELS), U(L_LEVELS)
    63       real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
    64 
    65       real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
    66 
    67 !     Variable species mixing ratio variables
    68       real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    69       real*8 KCOEF(4)
    70       integer NVAR(L_LEVELS)
    71 
    72 !     temporary variables for multiple aerosol calculation
    73       real*8 atemp, btemp, ctemp
    74 
    75 !     variables for k in units m^-1
    76       double precision wn_cont, p_cont, p_air, T_cont, dtemp
    77       real*8 dz(L_LEVELS), DCONT
    78 
    79       integer igas
    80 
    81 !=======================================================================
    82 !     Determine the total gas opacity throughout the column, for each
    83 !     spectral interval, NW, and each Gauss point, NG.
    84 !     Calculate the continuum opacities, i.e., those that do not depend on
    85 !     NG, the Gauss index.
    86 
    87       taugsurf(:,:) = 0.0
    88       dpr(:)        = 0.0
    89       lkcoef(:,:)   = 0.0
    90 
    91       do K=2,L_LEVELS
    92          DPR(k) = PLEV(K)-PLEV(K-1)
    93 
    94          ! if we have continuum opacities, we need dz
    95          if(kastprof)then
    96             dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
    97             U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
    98          else
    99             dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
    100             U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
    101          endif
    102 
    103          call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
    104               LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
    105 
    106          do LK=1,4
    107             LKCOEF(K,LK) = LCOEF(LK)
    108          end do
    109 
    110          DO NW=1,L_NSPECTV
    111             TRAY(K,NW)   = TAURAY(NW) * DPR(K)
    112 
    113             do iaer=1,naerkind
    114                TAEROS(K,NW,IAER) = TAUAERO(K,IAER)  * QXVAER(K,NW,IAER)
    115             end do
    116 !     
    117 
    118          END DO
    119       end do
    120 
    121 !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
    122 
    123 !     we ignore K=1...
    124       do K=2,L_LEVELS
    125          do NW=1,L_NSPECTV
    126 
    127             TRAYAER = TRAY(K,NW)
    128             do iaer=1,naerkind
    129                TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
    130             end do
    131 
    132             DCONT = 0.0 ! continuum absorption
    133 
    134             if(callgasvis.and.Continuum.and.(.not.graybody))then
    135             ! include continua if necessary
    136                wn_cont = dble(wnov(nw))
    137                T_cont  = dble(TMID(k))
    138                do igas=1,ngasmx
    139    
    140                   if(gfrac(igas).eq.-1)then ! variable
    141                      p_cont  = dble(PMID(k)*scalep*QVAR(K)) ! qvar = mol/mol
    142                   else
    143                      p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
    144                   endif
    145 
    146                   dtemp=0.0
    147                   if(gnom(igas).eq.'H2_')then
    148                      call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
    149                   elseif(gnom(igas).eq.'H2O'.and.T_cont.gt.200.0)then
    150                      p_air = dble(PMID(k)*scalep) - p_cont ! note assumes air!!
    151                      call interpolateH2Ocont(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
    152                   endif
    153 
    154                   DCONT = DCONT + dtemp
    155 
    156                enddo
    157 
    158                DCONT = DCONT*dz(k)
    159             endif
    160 
    161             do NG=1,L_NGAUSS-1
    162 
    163 !=======================================================================
    164 !     Now compute TAUGAS
    165 !     Interpolate between water mixing ratios
    166 !     WRATIO = 0.0 if the requested water amount is equal to, or outside the
    167 !     the water data range
    168 
    169                if (L_REFVAR.eq.1)then ! added by RW for special no variable case
    170                   KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
    171                   KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
    172                   KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
    173                   KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
    174                else
    175                   KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
    176                        (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
    177                        GASV(MT(K),MP(K),NVAR(K),NW,NG))
    178 
    179                   KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
    180                        (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
    181                        GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
    182 
    183                   KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
    184                        (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
    185                        GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
    186 
    187                   KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
    188                        (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
    189                        GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
    190                endif
    191 
    192 !     Interpolate the gaseous k-coefficients to the requested T,P values
    193                
    194                ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +           &
    195                     LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
    196 
    197                TAUGAS          = U(k)*ANS
    198 
    199 
    200                !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
    201                TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    202                DTAUKV(K,nw,ng) = TAUGAS + TRAYAER  & ! TRAYAER includes all scattering contributions
    203                                  + DCONT             ! for continuum absorption
    204 
    205             end do
    206 
    207 
    208 !     Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
    209 !     which holds continuum opacity only
    210 
    211             NG = L_NGAUSS
    212             DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
    213             do iaer=1,naerkind
    214                DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
    215 !     &                           + DCONT          ! For parameterized continuum absorption
    216             end do ! a bug was here!
    217 
    218          end do
    219       end do
    220 
    221 
    222 !=======================================================================
    223 !     Now the full treatment for the layers, where besides the opacity
    224 !     we need to calculate the scattering albedo and asymmetry factors
    225 
    226       DO NW=1,L_NSPECTV
    227          DO K=2,L_LEVELS
    228             do iaer=1,naerkind
    229               TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
    230             end do
    231          ENDDO
    232       ENDDO
    233      
    234 
    235       DO NW=1,L_NSPECTV
    236          DO NG=1,L_NGAUSS
    237             DO L=1,L_NLAYRAD-1
    238                K              = 2*L+1
    239 
    240                DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
    241 
    242                atemp=0.
    243                btemp=TRAY(K,NW) + TRAY(K+1,NW)
    244                ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
    245                do iaer=1,naerkind
    246                   atemp = atemp +                                     &
    247                        GVAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
    248                        GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
    249                   btemp = btemp +                                     &
    250                        TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
    251                   ctemp = ctemp +                                     &
    252                        TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
    253                end do
    254 
    255                COSBV(L,NW,NG) = atemp/btemp
    256                WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
    257 
    258             END DO
    259 
    260 !     No vertical averaging on bottom layer
    261 
    262             L              = L_NLAYRAD
    263             K              = 2*L+1
    264             DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
    265 
    266             atemp=0.
    267             btemp=TRAY(K,NW)
    268             ctemp=0.9999*TRAY(K,NW)
    269             do iaer=1,naerkind
    270                atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
    271                btemp = btemp + TAUAEROLK(K,NW,IAER)
    272                ctemp = ctemp + TAUAEROLK(K,NW,IAER)
    273             end do
    274             COSBV(L,NW,NG) = atemp/btemp
    275             WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
    276 
    277          END DO                 ! NG gauss point loop
    278       END DO                    ! NW spectral loop
    279 
    280 
    281 
    282 !     Total extinction optical depths
    283 
    284       DO NW=1,L_NSPECTV       
    285          DO NG=1,L_NGAUSS       ! full gauss loop
    286             TAUV(1,NW,NG)=0.0D0
    287             DO L=1,L_NLAYRAD
    288                TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
    289             END DO
    290 
    291             TAUCUMV(1,NW,NG)=0.0D0
    292             DO K=2,L_LEVELS
    293                TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
    294             END DO
    295          END DO                 ! end full gauss loop
    296       END DO
    297 
    298 
    299       RETURN
    300     END SUBROUTINE OPTCV
     40  real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     41  real*8 DTAUKV(L_LEVELS+1,L_NSPECTV,L_NGAUSS)
     42  real*8 TAUV(L_NLEVRAD,L_NSPECTV,L_NGAUSS)
     43  real*8 TAUCUMV(L_LEVELS,L_NSPECTV,L_NGAUSS)
     44  real*8 PLEV(L_LEVELS)
     45  real*8 TMID(L_LEVELS), PMID(L_LEVELS)
     46  real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     47  real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
     48  real*8 TAURAY(L_NSPECTV)
     49
     50  ! for aerosols
     51  real*8 QXVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
     52  real*8 QSVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
     53  real*8 GVAER(L_LEVELS+1,L_NSPECTV,NAERKIND)
     54  real*8 TAUAERO(L_LEVELS+1,NAERKIND)
     55  real*8 TAUAEROLK(L_LEVELS+1,L_NSPECTV,NAERKIND)
     56  real*8 TAEROS(L_LEVELS,L_NSPECTV,NAERKIND)
     57
     58  integer L, NW, NG, K, NG1(L_NSPECTV), LK, IAER
     59  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
     60  real*8  ANS, TAUGAS
     61  real*8  TRAY(L_LEVELS,L_NSPECTV)
     62  real*8  DPR(L_LEVELS), U(L_LEVELS)
     63  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
     64
     65  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
     66
     67  ! variable species mixing ratio variables
     68  real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
     69  real*8 KCOEF(4)
     70  integer NVAR(L_LEVELS)
     71
     72  ! temporary variables for multiple aerosol calculation
     73  real*8 atemp, btemp, ctemp
     74
     75  ! variables for k in units m^-1
     76  double precision wn_cont, p_cont, p_air, T_cont, dtemp
     77  double precision p_cross
     78  real*8 dz(L_LEVELS), DCONT
     79
     80  integer igas, jgas
     81
     82  !=======================================================================
     83  !     Determine the total gas opacity throughout the column, for each
     84  !     spectral interval, NW, and each Gauss point, NG.
     85  !     Calculate the continuum opacities, i.e., those that do not depend on
     86  !     NG, the Gauss index.
     87
     88  taugsurf(:,:) = 0.0
     89  dpr(:)        = 0.0
     90  lkcoef(:,:)   = 0.0
     91
     92  do K=2,L_LEVELS
     93     DPR(k) = PLEV(K)-PLEV(K-1)
     94
     95     ! if we have continuum opacities, we need dz
     96     if(kastprof)then
     97        dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
     98        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
     99     else
     100        dz(k) = dpr(k)*R*TMID(K)/(g*PMID(K))
     101        U(k)  = Cmk*DPR(k)    ! only Cmk line in optci.F
     102     endif
     103
     104     call tpindex(PMID(K),TMID(K),QVAR(K),pfgasref,tgasref,WREFVAR, &
     105          LCOEF,MT(K),MP(K),NVAR(K),WRATIO(K))
     106
     107     do LK=1,4
     108        LKCOEF(K,LK) = LCOEF(LK)
     109     end do
     110
     111     DO NW=1,L_NSPECTV
     112        TRAY(K,NW)   = TAURAY(NW) * DPR(K)
     113
     114        do iaer=1,naerkind
     115           TAEROS(K,NW,IAER) = TAUAERO(K,IAER)  * QXVAER(K,NW,IAER)
     116        end do
     117
     118     END DO
     119  end do
     120
     121  !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
     122
     123  !     we ignore K=1...
     124  do K=2,L_LEVELS
     125     do NW=1,L_NSPECTV
     126
     127        TRAYAER = TRAY(K,NW)
     128        do iaer=1,naerkind
     129           TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
     130        end do
     131
     132        DCONT = 0.0 ! continuum absorption
     133
     134        if(callgasvis.and.continuum.and.(.not.graybody))then
     135           ! include continua if necessary
     136           wn_cont = dble(wnov(nw))
     137           T_cont  = dble(TMID(k))
     138           do igas=1,ngasmx
     139
     140              if(gfrac(igas).eq.-1)then ! variable
     141                 p_cont  = dble(PMID(k)*scalep*QVAR(k)) ! qvar = mol/mol
     142              else
     143                 p_cont  = dble(PMID(k)*scalep*gfrac(igas)*(1.-QVAR(k)))
     144              endif
     145
     146              dtemp=0.0
     147              if(igas.eq.igas_N2)then
     148
     149                 !call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
     150                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
     151
     152              elseif(igas.eq.igas_H2)then
     153
     154                 ! first do self-induced absorption
     155                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
     156
     157                 ! then cross-interactions with other gases
     158                 do jgas=1,ngasmx
     159                    p_cross = dble(PMID(k)*scalep*gfrac(jgas)*(1.-QVAR(k)))
     160                    if(jgas.eq.igas_N2)then
     161                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
     162                       ! should be irrelevant in the visible
     163                    elseif(jgas.eq.igas_He)then
     164                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
     165                    endif
     166                 enddo
     167
     168              elseif(igas.eq.igas_H2O.and.T_cont.gt.200.0)then
     169
     170                 p_air = dble(PMID(k)*scalep) - p_cont ! note assumes background is air!
     171                 if(H2Ocont_simple)then
     172                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     173                 else
     174                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     175                 endif
     176
     177              endif
     178
     179              DCONT = DCONT + dtemp
     180
     181           enddo
     182
     183           DCONT = DCONT*dz(k)
     184        endif
     185
     186        do NG=1,L_NGAUSS-1
     187
     188           !=======================================================================
     189           !     Now compute TAUGAS
     190           !     Interpolate between water mixing ratios
     191           !     WRATIO = 0.0 if the requested water amount is equal to, or outside the
     192           !     the water data range
     193
     194           if (L_REFVAR.eq.1)then ! added by RW for special no variable case
     195              KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
     196              KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
     197              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,1,NW,NG)
     198              KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
     199           else
     200              KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
     201                   (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
     202                   GASV(MT(K),MP(K),NVAR(K),NW,NG))
     203
     204              KCOEF(2) = GASV(MT(K),MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*  &
     205                   (GASV(MT(K),MP(K)+1,NVAR(K)+1,NW,NG) -                &
     206                   GASV(MT(K),MP(K)+1,NVAR(K),NW,NG))
     207
     208              KCOEF(3) = GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG) + WRATIO(K)*&
     209                   (GASV(MT(K)+1,MP(K)+1,NVAR(K)+1,NW,NG) -              &
     210                   GASV(MT(K)+1,MP(K)+1,NVAR(K),NW,NG))
     211
     212              KCOEF(4) = GASV(MT(K)+1,MP(K),NVAR(K),NW,NG) + WRATIO(K)*  &
     213                   (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
     214                   GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
     215           endif
     216
     217           !     Interpolate the gaseous k-coefficients to the requested T,P values
     218
     219           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +           &
     220                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
     221
     222           TAUGAS = U(k)*ANS
     223
     224           !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
     225           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
     226           DTAUKV(K,nw,ng) = TAUGAS + TRAYAER  & ! TRAYAER includes all scattering contributions
     227                + DCONT             ! for continuum absorption
     228
     229        end do
     230
     231
     232        !     Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
     233        !     which holds continuum opacity only
     234
     235        NG = L_NGAUSS
     236        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
     237        do iaer=1,naerkind
     238           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
     239           !     &                           + DCONT          ! For parameterized continuum absorption
     240        end do ! a bug was here!
     241
     242     end do
     243  end do
     244
     245
     246  !=======================================================================
     247  !     Now the full treatment for the layers, where besides the opacity
     248  !     we need to calculate the scattering albedo and asymmetry factors
     249
     250  DO NW=1,L_NSPECTV
     251     DO K=2,L_LEVELS
     252        do iaer=1,naerkind
     253           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
     254        end do
     255     ENDDO
     256  ENDDO
     257
     258
     259  DO NW=1,L_NSPECTV
     260     DO NG=1,L_NGAUSS
     261        DO L=1,L_NLAYRAD-1
     262           K              = 2*L+1
     263
     264           DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)+DTAUKV(K+1,NW,NG)
     265
     266           atemp=0.
     267           btemp=TRAY(K,NW) + TRAY(K+1,NW)
     268           ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
     269           do iaer=1,naerkind
     270              atemp = atemp +                                     &
     271                   GVAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
     272                   GVAER(K+1,NW,IAER) * TAUAEROLK(K+1,NW,IAER)
     273              btemp = btemp +                                     &
     274                   TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
     275              ctemp = ctemp +                                     &
     276                   TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
     277           end do
     278
     279           COSBV(L,NW,NG) = atemp/btemp
     280           WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
     281
     282        END DO
     283
     284        !     No vertical averaging on bottom layer
     285
     286        L              = L_NLAYRAD
     287        K              = 2*L+1
     288        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG)
     289
     290        atemp=0.
     291        btemp=TRAY(K,NW)
     292        ctemp=0.9999*TRAY(K,NW)
     293        do iaer=1,naerkind
     294           atemp = atemp + GVAER(K,NW,IAER) * TAUAEROLK(K,NW,IAER)
     295           btemp = btemp + TAUAEROLK(K,NW,IAER)
     296           ctemp = ctemp + TAUAEROLK(K,NW,IAER)
     297        end do
     298        COSBV(L,NW,NG) = atemp/btemp
     299        WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
     300
     301     END DO                 ! NG gauss point loop
     302  END DO                    ! NW spectral loop
     303
     304
     305
     306  ! Total extinction optical depths
     307
     308  DO NW=1,L_NSPECTV       
     309     DO NG=1,L_NGAUSS       ! full gauss loop
     310        TAUV(1,NW,NG)=0.0D0
     311        DO L=1,L_NLAYRAD
     312           TAUV(L+1,NW,NG)=TAUV(L,NW,NG)+DTAUV(L,NW,NG)
     313        END DO
     314
     315        TAUCUMV(1,NW,NG)=0.0D0
     316        DO K=2,L_LEVELS
     317           TAUCUMV(K,NW,NG)=TAUCUMV(K-1,NW,NG)+DTAUKV(K,NW,NG)
     318        END DO
     319     END DO                 ! end full gauss loop
     320  END DO
     321
     322
     323  RETURN
     324END SUBROUTINE OPTCV

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

    r305 r716  
    66      implicit none
    77
    8       double precision rc,m_n,m_v,rmn
    9       save m_n,m_v,rmn
     8      double precision rc,cp_n,m_n,m_v,rmn
     9      save cp_n,m_n,m_v,rmn
    1010
    1111      logical ideal_v

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

    r651 r716  
    100100!           To F90: R. Wordsworth (2010)
    101101!           New turbulent diffusion scheme: J. Leconte (2012)
     102!           Loops converted to F90 matrix format: J. Leconte (2012)
    102103!     
    103104!==================================================================
    … …  
    690691              if(CLFvarying)then
    691692
    692                  !!! ---> PROBLEMS WITH ALLOCATED ARRAYS
    693                  !!! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
     693                 ! ---> PROBLEMS WITH ALLOCATED ARRAYS
     694                 ! (temporary solution in callcorrk: do not deallocate if CLFvarying ...)
    694695                 clearsky=.true.
    695696                 call callcorrk(ngrid,nlayer,pq,nq,qsurf,                  &
    … …  
    700701                      reffrad,tau_col1,cloudfrac,totcloudfrac,             &
    701702                      clearsky,firstcall,lastcall)
    702                  clearsky = .false.  !! just in case.     
     703                 clearsky = .false.  ! just in case.     
    703704
    704705                 ! Sum the fluxes and heating rates from cloudy/clear cases
    … …  
    18661867      icount=icount+1
    18671868
    1868       !!! DEALLOCATE STUFF
     1869      ! deallocate gas variables
    18691870      if (lastcall) then
    1870         IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )    !! this was allocated in su_gases.F90
    1871         IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  !! this was allocated in su_gases.F90
     1871        IF ( ALLOCATED( gnom ) ) DEALLOCATE( gnom )
     1872        IF ( ALLOCATED( gfrac ) ) DEALLOCATE( gfrac )  ! both allocated in su_gases.F90
    18721873      endif
    18731874

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

    r671 r716  
    9090      integer, parameter :: nsizemax = 60
    9191
    92       character (len=20) :: corrkdir
     92      character (len=100) :: corrkdir
    9393      save corrkdir
    9494
    95       character (len=30) :: banddir
     95      character (len=100) :: banddir
    9696      save banddir
    9797

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

    r543 r716  
    3636
    3737      character(len=30)  :: temp1
    38       character(len=100) :: file_id
    39       character(len=100) :: file_path
     38      character(len=200) :: file_id
     39      character(len=200) :: file_path
    4040
    4141!     C1 and C2 values from Goody and Yung (2nd edition)  MKS units

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

    r484 r716  
    3838
    3939      character(len=30)  :: temp1
    40       character(len=100) :: file_id
    41       character(len=100) :: file_path
     40      character(len=200) :: file_id
     41      character(len=200) :: file_path
    4242
    4343      real*8 :: lastband(2)

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

    r471 r716  
    55  implicit none
    66
    7        integer igas, ierr
     7  integer igas, ierr, count
    88
    9 !==================================================================
    10 !     
    11 !     Purpose
    12 !     -------
    13 !     Load atmospheric composition info
    14 !     
    15 !     Authors
    16 !     -------
    17 !     R. Wordsworth (2011)
    18 !     
    19 !==================================================================
     9  !==================================================================
     10  !     
     11  !     Purpose
     12  !     -------
     13  !     Load atmospheric composition info
     14  !     
     15  !     Authors
     16  !     -------
     17  !     R. Wordsworth (2011)
     18  !     Allocatable arrays by A. Spiga (2011)
     19  !     
     20  !==================================================================
    2021
    21        ! load gas names from file 'gases.def'
    22        open(90,file='gases.def',status='old',form='formatted',iostat=ierr)
    23        if (ierr.eq.0) then
    24           write(*,*) "sugases.F90: reading file gases.def"
    25           read(90,*)
    26           read(90,*,iostat=ierr) ngasmx
    27           if (ierr.ne.0) then
    28              write(*,*) "sugases.F90: error reading number of gases"
    29              write(*,*) "   (first line of gases.def) "
    30              call abort
    31           endif
     22  ! load gas names from file 'gases.def'
     23  open(90,file='gases.def',status='old',form='formatted',iostat=ierr)
     24  if (ierr.eq.0) then
     25     write(*,*) "sugases.F90: reading file gases.def"
     26     read(90,*)
     27     read(90,*,iostat=ierr) ngasmx
     28     if (ierr.ne.0) then
     29        write(*,*) "sugases.F90: error reading number of gases"
     30        write(*,*) "   (first line of gases.def) "
     31        call abort
     32     endif
    3233
    33           PRINT *, "OK I HAVE ",ngasmx, " GASES. NOW I ALLOCATE ARRAYS AND READ NAMES AND FRAC."
     34     print*,ngasmx, " gases found in gases.def. Allocating names and molar fractions..."
    3435
    35           IF ( .NOT. ALLOCATED( gnom ) ) ALLOCATE( gnom( ngasmx ) )
    36           do igas=1,ngasmx
    37              read(90,*,iostat=ierr) gnom(igas)
    38              if (ierr.ne.0) then
    39                 write(*,*) 'sugases.F90: error reading gas names in gases.def...'
    40                 call abort
    41              endif
    42           enddo                  !of do igas=1,ngasmx
    43          
    44           vgas=0
    45           IF ( .NOT. ALLOCATED( gfrac ) ) ALLOCATE( gfrac( ngasmx ) )
    46           do igas=1,ngasmx
    47              read(90,*,iostat=ierr) gfrac(igas)
    48              if (ierr.ne.0) then
    49                 write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...'
    50                 call abort
    51              endif
     36     if (.not.allocated(gnom)) allocate(gnom(ngasmx))
     37     do igas=1,ngasmx
     38        read(90,*,iostat=ierr) gnom(igas)
     39        if (ierr.ne.0) then
     40           write(*,*) 'sugases.F90: error reading gas names in gases.def...'
     41           call abort
     42        endif
     43     enddo                  !of do igas=1,ngasmx
    5244
    53              ! find variable gas (if any)
    54              if(gfrac(igas).eq.-1.0)then
    55                 if(vgas.eq.0)then
    56                    vgas=igas
    57                 else
    58                    print*,'You seem to be choosing two variable gases'
    59                    print*,'Check that gases.def is correct'
    60                    call abort
    61                 endif
    62              endif
    63              
    64           enddo                  !of do igas=1,ngasmx
     45     vgas=0
     46     if(.not.allocated(gfrac)) allocate(gfrac(ngasmx))
     47     do igas=1,ngasmx
     48        read(90,*,iostat=ierr) gfrac(igas)
     49        if (ierr.ne.0) then
     50           write(*,*) 'sugases.F90: error reading gas molar fractions in gases.def...'
     51           call abort
     52        endif
    6553
    66        else
    67           write(*,*) 'Cannot find required file "gases.def"'
    68           call abort
    69        endif
    70        close(90)
     54        ! find variable gas (if any)
     55        if(gfrac(igas).eq.-1.0)then
     56           if(vgas.eq.0)then
     57              vgas=igas
     58           else
     59              print*,'You seem to be choosing two variable gases'
     60              print*,'Check that gases.def is correct'
     61              call abort
     62           endif
     63        endif
     64
     65     enddo                  !of do igas=1,ngasmx
     66
     67
     68     ! assign the 'igas_X' labels
     69     count=0
     70     do igas=1,ngasmx
     71        if (gnom(igas).eq."H2_") then
     72           igas_H2=igas
     73           count=count+1
     74        elseif (gnom(igas).eq."He_") then
     75           igas_He=igas
     76           count=count+1
     77        elseif (gnom(igas).eq."H2O") then
     78           igas_H2O=igas
     79           count=count+1
     80        elseif (gnom(igas).eq."CO2") then
     81           igas_CO2=igas
     82           count=count+1
     83        elseif (gnom(igas).eq."CO_") then
     84           igas_CO=igas
     85           count=count+1
     86        elseif (gnom(igas).eq."N2_") then
     87           igas_N2=igas
     88           count=count+1
     89        elseif (gnom(igas).eq."O2_") then
     90           igas_O2=igas
     91           count=count+1
     92        elseif (gnom(igas).eq."SO2") then
     93           igas_SO2=igas
     94           count=count+1
     95        elseif (gnom(igas).eq."H2S") then
     96           igas_H2S=igas
     97           count=count+1
     98        elseif (gnom(igas).eq."CH4") then
     99           igas_CH4=igas
     100           count=count+1
     101        elseif (gnom(igas).eq."NH3") then
     102           igas_NH3=igas
     103           count=count+1
     104        endif
     105     enddo
     106
     107     if(count.ne.ngasmx)then
     108        print*,'Mismatch between ngas and number of recognised gases in sugas_corrk.F90.'
     109        print*,'Either we haven`t managed to assign all the gases, or there are duplicates.'
     110        print*,'Please try again.'
     111     endif
     112
     113  else
     114     write(*,*) 'Cannot find required file "gases.def"'
     115     call abort
     116  endif
     117  close(90)
    71118
    72119end subroutine su_gases

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

    r596 r716  
    1313!     Adapted and generalised from the NASA Ames code by Robin Wordsworth (2009)
    1414!     Added double gray case by Jeremy Leconte (2012)
     15!     New HITRAN continuum data section by RW (2012)
    1516!
    1617!     Summary
    … …  
    2526      use radcommon_h, only : wrefvar
    2627      use datafile_mod, only: datadir
     28
    2729      use gases_h
    2830      use ioipsl_getincom
    … …  
    3133#include "callkeys.h"
    3234#include "comcstfi.h"
     35
    3336!==================================================================
    3437
    … …  
    3841      integer L_NGAUSScheck
    3942
    40       character(len=100) :: file_id
    41       character(len=300) :: file_path
    42 
    43       !! ALLOCATABLE ARRAYS -- AS 12/2011
     43      character(len=200) :: file_id
     44      character(len=500) :: file_path
     45
     46      ! ALLOCATABLE ARRAYS -- AS 12/2011
    4447      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi8, gasv8
    4548      character*3,allocatable,DIMENSION(:) :: gastype ! for check with gnom
    … …  
    4851      real kappa_IR, kappa_VI
    4952
    50       integer ngas, igas
     53      integer ngas, igas, jgas
    5154
    5255      double precision testcont ! for continuum absorption initialisation
    … …  
    8487                '] has no variable species, exiting.'
    8588         call abort
    86       elseif(ngas.eq.2 .and. (.not.varactive) .and. (.not.varfixed))then
    87          print*,'You have varactive and varfixed=.false. and the database [', &
    88                      corrkdir(1:LEN_TRIM(corrkdir)),           &
    89                 '] has a variable species.'
    90          call abort
    91       elseif(ngas.gt.4 .or. ngas.lt.1)then
     89      elseif(ngas.gt.5 .or. ngas.lt.1)then
    9290         print*,ngas,' species in database [',               &
    9391                     corrkdir(1:LEN_TRIM(corrkdir)),           &
    … …  
    9694      endif
    9795
    98       if(ngas.gt.3)then
    99          print*,'ngas>3, EXPERIMENTAL!'
    100       endif
    101 
     96      ! dynamically allocate gastype and read from Q.dat
    10297      IF ( .NOT. ALLOCATED( gastype ) ) ALLOCATE( gastype( ngas ) )
    103       do n=1,ngas
    104          read(111,*) gastype(n)
    105          print*,'Gas ',n,' is ',gastype(n)
     98
     99      do igas=1,ngas
     100         read(111,*) gastype(igas)
     101         print*,'Gas ',igas,' is ',gastype(igas)
    106102      enddo
    107103
    108       ! GET array size, load the coefficients
     104      ! get array size, load the coefficients
    109105      open(111,file=TRIM(file_path),form='formatted')
    110106      read(111,*) L_REFVAR
    … …  
    113109      close(111)
    114110
     111      if(L_REFVAR.gt.1 .and. (.not.varactive) .and. (.not.varfixed))then
     112         print*,'You have varactive and varfixed=.false. and the database [', &
     113                     corrkdir(1:LEN_TRIM(corrkdir)),           &
     114                '] has a variable species.'
     115         call abort
     116      endif
     117
    115118      ! Check that gastype and gnom match
    116       do n=1,ngas
    117          print*,'Gas ',n,' is ',gnom(n)
    118          if(gnom(n).ne.gastype(n))then
    119             print*,'Name of a gas in radiative transfer data (',gastype(n),') does not ', &
    120                  'match that in gases.def (',gnom(n),'), exiting. You should compare ', &
     119      do igas=1,ngas
     120         print*,'Gas ',igas,' is ',gnom(igas)
     121         if(gnom(igas).ne.gastype(igas))then
     122            print*,'Name of a gas in radiative transfer data (',gastype(igas),') does not ', &
     123                 'match that in gases.def (',gnom(igas),'), exiting. You should compare ', &
    121124                 'gases.def with Q.dat in your radiative transfer directory.'
    122125            call abort
    … …  
    126129
    127130      ! display the values
    128       print*,'Variable gas mixing ratios:'
     131      print*,'Variable gas volume mixing ratios:'
    129132      do n=1,L_REFVAR
    130133         !print*,n,'.',wrefvar(n),' kg/kg' ! pay attention!
    … …  
    191194      endif
    192195     
    193       ! GET array size, load the coefficients
     196      ! get array size, load the coefficients
    194197      open(111,file=TRIM(file_path),form='formatted')
    195198      read(111,*) L_NPREF
    … …  
    200203      IF( .NOT. ALLOCATED( pfgasref ) ) ALLOCATE( PFGASREF(L_PINT) )
    201204
    202 
    203205      ! display the values
    204206      print*,'Correlated-k pressure grid (mBar):'
    … …  
    219221      end do
    220222      pfgasref(L_PINT) = pgasref(L_NPREF)
    221       print*,'Warning, pfgasref needs generalisation to uneven grids!!'
     223      ! Warning, pfgasref needs generalisation to uneven grids!
    222224
    223225!-----------------------------------------------------------------------
    … …  
    239241      endif
    240242
    241       ! GET array size, load the coefficients
     243      ! get array size, load the coefficients
    242244      open(111,file=TRIM(file_path),form='formatted')
    243245      read(111,*) L_NTREF
    … …  
    256258      tgasmax = tgasref(L_NTREF)
    257259
    258 
    259260!-----------------------------------------------------------------------
    260 !-----------------------------------------------------------------------
    261 ! ALLOCATE THE MULTIDIM ARRAYS IN radcommon_h
    262       PRINT *, L_NTREF,L_NPREF,L_REFVAR
     261! allocate the multidimensional arrays in radcommon_h
     262
    263263      IF( .NOT. ALLOCATED( gasi8 ) ) ALLOCATE( gasi8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTI,L_NGAUSS) )
    264264      IF( .NOT. ALLOCATED( gasv8 ) ) ALLOCATE( gasv8(L_NTREF,L_NPREF,L_REFVAR,L_NSPECTV,L_NGAUSS) )
    265265      IF( .NOT. ALLOCATED( gasi ) ) ALLOCATE( gasi(L_NTREF,L_PINT,L_REFVAR,L_NSPECTI,L_NGAUSS) )
    266266      IF( .NOT. ALLOCATED( gasv ) ) ALLOCATE( gasv(L_NTREF,L_PINT,L_REFVAR,L_NSPECTV,L_NGAUSS) )
    267 !-----------------------------------------------------------------------
    268 !-----------------------------------------------------------------------
     267
     268      ! display the values
     269      print*,''
     270      print*,'Correlated-k matrix size:'
     271      print*,'[',L_NTREF,',',L_NPREF,',',L_REFVAR,',',L_NGAUSS,']'
    269272
    270273!=======================================================================
    … …  
    288291         read(111,*) gasv8
    289292         close(111)
    290          
     293
    291294      else if (callgasvis.and.graybody) then
    292 !    constant absorption coefficient in visible
     295         ! constant absorption coefficient in visible
    293296         write(*,*)"constant absorption coefficient in visible:"
    294          kappa_VI=-100000.
     297         kappa_VI = -100000.
    295298         call getin("kappa_VI",kappa_VI)
    296299         write(*,*)" kappa_VI = ",kappa_VI
    297          kappa_VI=kappa_VI*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule         
    298          gasv8(:,:,:,:,:)=kappa_VI
     300         kappa_VI = kappa_VI * 1.e4 * mugaz * 1.672621e-27       ! conversion from m^2/kg to cm^2/molecule         
     301         gasv8(:,:,:,:,:) = kappa_VI
     302
    299303      else
    300304         print*,'Visible corrk gaseous absorption is set to zero.'
    301          gasv8(:,:,:,:,:)=0.0
     305         gasv8(:,:,:,:,:) = 0.0
     306
    302307      endif
    303308
    … …  
    319324         read(111,*) gasi8
    320325         close(111)
     326     
     327         ! 'fzero' is a currently unused feature that allows optimisation
     328         ! of the radiative transfer by neglecting bands where absorption
     329         ! is close to zero. As it could be useful in the future, this
     330         ! section of the code has been kept commented and not erased.
     331         ! RW 7/3/12.
     332
     333         do nw=1,L_NSPECTI
     334            fzeroi(nw) = 0.
     335!            do nt=1,L_NTREF
     336!               do np=1,L_NPREF
     337!                  do nh=1,L_REFVAR
     338!                     do ng = 1,L_NGAUSS
     339!                        if(gasi8(nt,np,nh,nw,ng).lt.1.0e-25)then
     340!                           fzeroi(nw)=fzeroi(nw)+1.
     341!                        endif
     342!                     end do
     343!                  end do
     344!               end do
     345!            end do
     346!            fzeroi(nw)=fzeroi(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
     347         end do
     348
     349         do nw=1,L_NSPECTV
     350            fzerov(nw) = 0.
     351!            do nt=1,L_NTREF
     352!               do np=1,L_NPREF
     353!                  do nh=1,L_REFVAR
     354!                     do ng = 1,L_NGAUSS
     355!                        if(gasv8(nt,np,nh,nw,ng).lt.1.0e-25)then
     356!                           fzerov(nw)=fzerov(nw)+1.
     357!                        endif
     358!                     end do
     359!                  end do
     360!               end do
     361!            end do
     362!            fzerov(nw)=fzerov(nw)/dble(L_NTREF*L_NPREF*L_REFVAR*L_NGAUSS)
     363         end do
    321364
    322365      else if (graybody) then
    323 !    constant absorption coefficient in IR
     366         ! constant absorption coefficient in IR
    324367         write(*,*)"constant absorption coefficient in InfraRed:"
    325          kappa_IR=-100000.
     368         kappa_IR = -100000.
    326369         call getin("kappa_IR",kappa_IR)
    327          write(*,*)" kappa_IR = ",kappa_IR       
    328          kappa_IR=kappa_IR*1.e4* mugaz * 1.672621e-27    ! conversion from m^2/kg to cm^2/molecule       
    329          gasi8(:,:,:,:,:)=kappa_IR
     370         write(*,*)" kappa_IR = ",kappa_IR
     371         kappa_IR = kappa_IR * 1.e4 * mugaz * 1.672621e-27       ! conversion from m^2/kg to cm^2/molecule       
     372         gasi8(:,:,:,:,:) = kappa_IR
     373
    330374      else
    331375         print*,'Infrared corrk gaseous absorption is set to zero.'
    332          gasi8(:,:,:,:,:)=0.0
    333       endif
    334 
    335       do nw=1,L_NSPECTI
    336          fzeroi(nw) = 0.
    337       end do
    338 
    339       do nw=1,L_NSPECTV
    340          fzerov(nw) = 0.
    341       end do
     376         gasi8(:,:,:,:,:) = 0.0
     377
     378      endif
    342379
    343380
    … …  
    512549
    513550
     551!=======================================================================
     552!     Initialise the continuum absorption data
    514553
    515554      do igas=1,ngasmx
    516          if(gnom(igas).eq.'H2_')then
     555
     556         if(gnom(igas).eq.'N2_')then
     557
     558            call interpolateN2N2(100.D+0,250.D+0,17500.D+0,testcont,.true.)
     559
     560         elseif(gnom(igas).eq.'H2_')then
     561
     562            ! first do self-induced absorption
    517563            call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
    518          elseif(gnom(igas).eq.'H2O')then
    519             call interpolateH2Ocont(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
     564            ! then cross-interactions with other gases
     565            do jgas=1,ngasmx
     566               if(gnom(jgas).eq.'N2_')then
     567                  call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.)
     568               elseif(gnom(jgas).eq.'He_')then
     569                  call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.)
     570               endif
     571            enddo
     572
     573         elseif(gnom(igas).eq.'H2O')then
     574
     575            ! H2O is special
     576            if(H2Ocont_simple)then
     577               call interpolateH2Ocont_PPC(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
     578            else
     579               call interpolateH2Ocont_CKD(990.D+0,296.D+0,683.2D+0*2,0.D+0,testcont,.true.)
     580            endif
     581
    520582         endif
     583
    521584      enddo
    522585
    523       !!! DEALLOCATE LOCAL ARRAYS
     586      print*,'----------------------------------------------------'
     587      print*,'And that`s all we have. It`s possible that other'
     588      print*,'continuum absorption may be present, but if it is we'
     589      print*,'don`t yet have data for it...'
     590      print*,''
     591
     592!     Deallocate local arrays
    524593      IF( ALLOCATED( gasi8 ) ) DEALLOCATE( gasi8 )
    525594      IF( ALLOCATED( gasv8 ) ) DEALLOCATE( gasv8 )
Note: See TracChangeset for help on using the changeset viewer.