Ignore:
Timestamp:
Feb 9, 2013, 5:48:51 PM (12 years ago)
Author:
aslmd
Message:

LMDZ.GENERIC (and LMDZ.UNIVERSAL)

  • Optimized calculations for continuum (done for H2 and He, to be done for others)
    • new common bilinear interpolation routine (bilinearbig)
    • optimization: only one calculation is actually needed

to find indexes of wavelength for bilinear interpolation
... because this will not change with level and integration step!

  • optimization: use while loop in bilinearbig
  • completely similar results obtained (test case for a gas giant, many simulated days)

NB: those changes really improve gcm speed (factor 2.2 for whole model!)

continuum was very expensive, now very cheap
--> e.g. 1 day, 25 dyn ts, 5 phys ts
--> before: 243 seconds (including 120 seconds for continuum bilinear interpolation)
--> after: 108 seconds

  • Corrected a bug: Continuum in inifis instead of continuum

... until now, most users (unbeknownst to them) were running with the continuum by default!

  • Cosmetic changes in optcv (mostly spaces and line breaks)

... so that comparisons with optci are easy e.g. through vimdiff

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

Legend:

Unmodified
Added
Removed

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

    r728 r873  
    88     &   , iaervar,iddist,topdustref,callstats,calleofdump              &
    99     &   , enertest                                                     &
    10      &   , callgasvis,Continuum,H2Ocont_simple,graybody                 &
     10     &   , callgasvis,continuum,H2Ocont_simple,graybody                 &
    1111     &   , Nmix_co2, radfixed, dusttau                                  &
    1212     &   , meanOLR, specOLR                                             &
    … …  
    3535     &   , season,diurnal,tlocked,lwrite                                &
    3636     &   , callstats,calleofdump                                        &
    37      &   , callgasvis,Continuum,H2Ocont_simple,graybody
     37     &   , callgasvis,continuum,H2Ocont_simple,graybody
    3838
    3939      logical enertest

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

    r804 r873  
    203203         write(*,*) "call continuum opacities in radiative transfer ?",
    204204     &              "(matters only if callrad=T)"
    205          Continuum=.true. ! default value
    206          call getin("Continuum",Continuum)
    207          write(*,*) " Continuum = ",Continuum
     205         continuum=.true. ! default value
     206         call getin("continuum",continuum)
     207         write(*,*) " continuum = ",continuum
    208208
    209209         write(*,*) "use analytic function for H2O continuum ?"

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

    r716 r873  
    1      subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall)
     1     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,a)
    22
    33!==================================================================
    … …  
    1515
    1616      use datafile_mod, only: datadir
     17
    1718      implicit none
    1819
    … …  
    5253      double precision blah, Ttemp
    5354      integer nres
     55
     56      integer a
    5457
    5558      if(temp.gt.400)then
    … …  
    104107         print*,'   pressure ',pres,' Pa'
    105108
    106          call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
     109      endif
    107110
    108          print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
    109          print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
     111         call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
     112
     113         !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
     114         !print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
    110115
    111116         abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
    112117
    113          print*,'We have ',amagat,' amagats of H2'
    114          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)
    135 
    136       else
    137  
    138          call bilinearH2H2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    139          abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
     118         !print*,'We have ',amagat,' amagats of H2'
     119         !print*,'So the absorption is ',abcoef,' m^-1'
    140120
    141121         ! unlike for Rayleigh scattering, we do not currently weight by the BB function
    142122         ! however our bands are normally thin, so this is no big deal.
    143123
    144       endif
    145 
    146124      return
    147125    end subroutine interpolateH2H2
    148 
    149 
    150 !-------------------------------------------------------------------------
    151       subroutine bilinearH2H2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
    152 !     Necessary for interpolation of continuum data
    153 
    154       implicit none
    155 
    156       integer nX,nY,i,j,a,b
    157       parameter(nX=2428)
    158       parameter(nY=10)
    159 
    160       real*8 x_in,y_in,x,y,x1,x2,y1,y2
    161       real*8 f,f11,f12,f21,f22,fA,fB
    162       real*8 x_arr(nX)
    163       real*8 y_arr(nY)
    164       real*8 f2d_arr(nX,nY)
    165      
    166       integer strlen
    167       character*100 label
    168       label='subroutine bilinear'
    169 
    170       x=x_in
    171       y=y_in
    172 
    173 
    174 !     1st check we're within the wavenumber range
    175       if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
    176          f=0.0D+0
    177          return
    178       else
    179          
    180 !     in the x (wavenumber) direction 1st
    181          i=1
    182  10      if (i.lt.(nX+1)) then
    183             if (x_arr(i).gt.x) then
    184                x1=x_arr(i-1)
    185                x2=x_arr(i)
    186                a=i-1
    187                i=9999
    188             endif
    189             i=i+1
    190             goto 10
    191          endif
    192       endif
    193      
    194       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
    195          write(*,*) 'Warning from bilinearH2H2:'
    196          write(*,*) 'Outside continuum temperature range!'
    197          if(y.lt.y_arr(1))then
    198             y=y_arr(1)+0.01
    199          endif
    200          if(y.gt.y_arr(nY))then
    201             y=y_arr(nY)-0.01
    202          endif
    203       else
    204 
    205 !     in the y (temperature) direction 2nd
    206          j=1
    207  20      if (j.lt.(nY+1)) then
    208             if (y_arr(j).gt.y) then
    209                y1=y_arr(j-1)
    210                y2=y_arr(j)
    211                b=j-1
    212                j=9999
    213             endif
    214             j=j+1
    215             goto 20
    216          endif
    217       endif
    218      
    219       f11=f2d_arr(a,b)
    220       f21=f2d_arr(a+1,b)
    221       f12=f2d_arr(a,b+1)
    222       f22=f2d_arr(a+1,b+1)
    223 
    224       call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
    225 
    226       return
    227     end subroutine bilinearH2H2

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

    r716 r873  
    1      subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall)
     1     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,a)
    22
    33!==================================================================
    … …  
    1515
    1616      use datafile_mod, only: datadir
     17
    1718      implicit none
    1819
    … …  
    5556      integer nres
    5657
    57 
     58      integer a
     59 
    5860      if(temp.gt.400)then
    5961         print*,'Your temperatures are too high for this H2-He CIA dataset.'
    … …  
    108110         print*,'   and He partial pressure     ',presHe,' Pa'
    109111
    110          call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
     112      endif
    111113
    112          print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
    113          print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
     114         call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
     115
     116         !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
     117         !print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
    114118
    115119         abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1
    116120
    117          print*,'We have ',amagatH2,' amagats of H2'
    118          print*,'and     ',amagatHe,' amagats of He'
    119          print*,'So the absorption is ',abcoef,' m^-1'
    120 
    121       else
    122  
    123          call bilinearH2He(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    124          abcoef=abcoef*losch**2*100.0*amagatH2*amagatHe ! convert to m^-1
     121         !print*,'We have ',amagatH2,' amagats of H2'
     122         !print*,'and     ',amagatHe,' amagats of He'
     123         !print*,'So the absorption is ',abcoef,' m^-1'
    125124
    126125         ! unlike for Rayleigh scattering, we do not currently weight by the BB function
    127126         ! however our bands are normally thin, so this is no big deal.
    128127
    129       endif
    130 
    131128      return
    132129    end subroutine interpolateH2He
    133 
    134 
    135 !-------------------------------------------------------------------------
    136       subroutine bilinearH2He(x_arr,y_arr,f2d_arr,x_in,y_in,f)
    137 !     Necessary for interpolation of continuum data
    138 
    139       implicit none
    140 
    141       integer nX,nY,i,j,a,b
    142       parameter(nX=2428)
    143       parameter(nY=10)
    144 
    145       real*8 x_in,y_in,x,y,x1,x2,y1,y2
    146       real*8 f,f11,f12,f21,f22,fA,fB
    147       real*8 x_arr(nX)
    148       real*8 y_arr(nY)
    149       real*8 f2d_arr(nX,nY)
    150      
    151       integer strlen
    152       character*100 label
    153       label='subroutine bilinear'
    154 
    155       x=x_in
    156       y=y_in
    157 
    158 !     1st check we're within the wavenumber range
    159       if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
    160          f=0.0D+0
    161          return
    162       else
    163          
    164 !     in the x (wavenumber) direction 1st
    165          i=1
    166  10      if (i.lt.(nX+1)) then
    167             if (x_arr(i).gt.x) then
    168                x1=x_arr(i-1)
    169                x2=x_arr(i)
    170                a=i-1
    171                i=9999
    172             endif
    173             i=i+1
    174             goto 10
    175          endif
    176       endif
    177      
    178       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
    179          write(*,*) 'Warning from bilinearH2He:'
    180          write(*,*) 'Outside continuum temperature range!'
    181          if(y.lt.y_arr(1))then
    182             y=y_arr(1)+0.01
    183          endif
    184          if(y.gt.y_arr(nY))then
    185             y=y_arr(nY)-0.01
    186          endif
    187       else
    188 
    189 !     in the y (temperature) direction 2nd
    190          j=1
    191  20      if (j.lt.(nY+1)) then
    192             if (y_arr(j).gt.y) then
    193                y1=y_arr(j-1)
    194                y2=y_arr(j)
    195                b=j-1
    196                j=9999
    197             endif
    198             j=j+1
    199             goto 20
    200          endif
    201       endif
    202      
    203       f11=f2d_arr(a,b)
    204       f21=f2d_arr(a+1,b)
    205       f12=f2d_arr(a,b+1)
    206       f22=f2d_arr(a+1,b+1)
    207 
    208       call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
    209 
    210      
    211       return
    212     end subroutine bilinearH2He

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

    r716 r873  
    44
    55  use radinc_h
    6   use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep
     6  use radcommon_h, only: gasi,tlimit,wrefVAR,Cmk,tgasref,pfgasref,wnoi,scalep,indi
    77  use gases_h
    88  implicit none
    … …  
    7070
    7171  ! variables for k in units m^-1
    72   real*8 rho, dz(L_LEVELS)
     72  real*8 dz(L_LEVELS)
     73  !real*8 rho !! see test below
    7374
    7475  integer igas, jgas
     76
     77  integer interm
    7578
    7679  !--- Kasting's CIA ----------------------------------------
    … …  
    8790  !----------------------------------------------------------
    8891
     92  !! AS: to save time in computing continuum (see bilinearbig)
     93  IF (.not.ALLOCATED(indi)) THEN
     94      ALLOCATE(indi(L_NSPECTI,ngasmx))
     95      indi = -9999  ! this initial value means "to be calculated"
     96  ENDIF
     97
    8998  !=======================================================================
    9099  !     Determine the total gas opacity throughout the column, for each
    … …  
    133142
    134143  do K=2,L_LEVELS
    135      do nw=1,L_NSPECTI
     144
     145     do NW=1,L_NSPECTI
    136146
    137147        DCONT = 0.0 ! continuum absorption
    138148
    139         if(Continuum.and.(.not.graybody))then
     149        if(continuum.and.(.not.graybody))then
    140150           ! include continua if necessary
    141151           wn_cont = dble(wnoi(nw))
    … …  
    157167
    158168                 ! first do self-induced absorption
    159                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
     169                 interm = indi(nw,igas)
     170                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     171                 indi(nw,igas) = interm
    160172
    161173                 ! then cross-interactions with other gases
    … …  
    166178                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
    167179                    elseif(jgas.eq.igas_He)then
    168                        call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
     180                       interm = indi(nw,jgas)
     181                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     182                       indi(nw,jgas) = interm
    169183                    endif
    170184                    dtemp = dtemp + dtempc
    … …  
    228242
    229243              KCOEF(1) = GASI(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*     &
    230                    (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   & 
     244                   (GASI(MT(K),MP(K),NVAR(K)+1,NW,NG) -                   &
    231245                   GASI(MT(K),MP(K),NVAR(K),NW,NG))
    232246
    … …  
    242256                   (GASI(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                 &
    243257                   GASI(MT(K)+1,MP(K),NVAR(K),NW,NG))
     258
    244259           endif
    245260
    … …  
    252267
    253268           TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS + DCONT
    254            DTAUKI(K,nw,ng) = TAUGAS + DCONT ! For parameterized continuum absorption
     269           DTAUKI(K,nw,ng) = TAUGAS &
     270                             + DCONT ! For parameterized continuum absorption
    255271
    256272           do iaer=1,naerkind
    … …  
    278294  !     we need to calculate the scattering albedo and asymmetry factors
    279295
     296  do iaer=1,naerkind
    280297  DO NW=1,L_NSPECTI
    281298     DO K=2,L_LEVELS+1
    282         do iaer=1,naerkind
    283299           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER)*QSIAER(K,NW,IAER)
    284         end do
    285300     ENDDO
    286301  ENDDO
     302  end do
    287303
    288304  DO NW=1,L_NSPECTI

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

    r716 r873  
    44
    55  use radinc_h
    6   use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep
     6  use radcommon_h, only: gasv, tlimit, wrefVAR, Cmk, tgasref, pfgasref,wnov,scalep,indv
    77  use gases_h
    88
    … …  
    3535  !-------------------------------------------------------------------
    3636
     37#include "comcstfi.h"
    3738#include "callkeys.h"
    38 #include "comcstfi.h"
     39
    3940
    4041  real*8 DTAUV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    … …  
    4647  real*8 COSBV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    4748  real*8 WBARV(L_NLAYRAD,L_NSPECTV,L_NGAUSS)
    48   real*8 TAURAY(L_NSPECTV)
    4949
    5050  ! 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
     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, LK, IAER
    5959  integer MT(L_LEVELS), MP(L_LEVELS), NP(L_LEVELS)
    6060  real*8  ANS, TAUGAS
     61  real*8  TAURAY(L_NSPECTV)
    6162  real*8  TRAY(L_LEVELS,L_NSPECTV)
     63  real*8  TRAYAER
    6264  real*8  DPR(L_LEVELS), U(L_LEVELS)
    6365  real*8  LCOEF(4), LKCOEF(L_LEVELS,4)
    6466
    65   real*8 taugsurf(L_NSPECTV,L_NGAUSS-1), TRAYAER
     67  real*8 taugsurf(L_NSPECTV,L_NGAUSS-1)
     68  real*8 DCONT
     69  double precision wn_cont, p_cont, p_air, T_cont, dtemp, dtempc
     70  double precision p_cross
    6671
    6772  ! variable species mixing ratio variables
    68   real*8 QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
    69   real*8 KCOEF(4)
     73  real*8  QVAR(L_LEVELS), WRATIO(L_LEVELS), MUVAR(L_LEVELS)
     74  real*8  KCOEF(4)
    7075  integer NVAR(L_LEVELS)
    7176
    … …  
    7479
    7580  ! 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
     81  real*8 dz(L_LEVELS)
    7982
    8083  integer igas, jgas
     84
     85  integer interm
     86
     87  !! AS: to save time in computing continuum (see bilinearbig)
     88  IF (.not.ALLOCATED(indv)) THEN
     89      ALLOCATE(indv(L_NSPECTV,ngasmx))
     90      indv = -9999 ! this initial value means "to be calculated"
     91  ENDIF
    8192
    8293  !=======================================================================
    … …  
    95106     ! if we have continuum opacities, we need dz
    96107     if(kastprof)then
    97         dz(k) = dpr(k)*(8314.5/muvar(k))*TMID(K)/(g*PMID(K))
     108        dz(k) = dpr(k)*(1000.0*8.3145/muvar(k))*TMID(K)/(g*PMID(K))
    98109        U(k)  = (Cmk*mugaz/(muvar(k)))*DPR(k)
    99110     else
    … …  
    109120     end do
    110121
     122
    111123     DO NW=1,L_NSPECTV
     124        do iaer=1,naerkind
     125           TAEROS(K,NW,IAER) = TAUAERO(K,IAER) * QXVAER(K,NW,IAER)
     126        end do
    112127        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 
    118128     END DO
    119   end do
    120 
    121   !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
     129  end do                    ! levels
    122130
    123131  !     we ignore K=1...
    124132  do K=2,L_LEVELS
     133
    125134     do NW=1,L_NSPECTV
    126135
    127136        TRAYAER = TRAY(K,NW)
     137        !     TRAYAER is Tau RAYleigh scattering, plus AERosol opacity
    128138        do iaer=1,naerkind
    129139           TRAYAER = TRAYAER + TAEROS(K,NW,IAER)
    … …  
    132142        DCONT = 0.0 ! continuum absorption
    133143
    134         if(callgasvis.and.continuum.and.(.not.graybody))then
     144        if(continuum.and.(.not.graybody).and.callgasvis)then
    135145           ! include continua if necessary
    136146           wn_cont = dble(wnov(nw))
    … …  
    153163
    154164                 ! first do self-induced absorption
    155                  call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.)
     165                 interm = indv(nw,igas)
     166                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     167                 indv(nw,igas) = interm
    156168
    157169                 ! then cross-interactions with other gases
    158170                 do jgas=1,ngasmx
    159171                    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.)
     172                    dtempc  = 0.0
     173                    if(jgas.eq.igas_N2)then
     174                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
    162175                       ! should be irrelevant in the visible
    163176                    elseif(jgas.eq.igas_He)then
    164                        call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtemp,.false.)
     177                       interm = indv(nw,jgas)
     178                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     179                       indv(nw,jgas) = interm
    165180                    endif
     181                    dtemp = dtemp + dtempc
    166182                 enddo
    167183
    … …  
    181197           enddo
    182198
    183            DCONT = DCONT*dz(k)
     199           DCONT = DCONT*dz(k)
     200
    184201        endif
    185202
    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
     203        do ng=1,L_NGAUSS-1
     204
     205           ! Now compute TAUGAS
     206
     207           ! Interpolate between water mixing ratios
     208           ! WRATIO = 0.0 if the requested water amount is equal to, or outside the
     209           ! the water data range
     210
     211           if(L_REFVAR.eq.1)then ! added by RW for special no variable case
    195212              KCOEF(1) = GASV(MT(K),MP(K),1,NW,NG)
    196213              KCOEF(2) = GASV(MT(K),MP(K)+1,1,NW,NG)
    … …  
    198215              KCOEF(4) = GASV(MT(K)+1,MP(K),1,NW,NG)
    199216           else
     217
    200218              KCOEF(1) = GASV(MT(K),MP(K),NVAR(K),NW,NG) + WRATIO(K)*    &
    201219                   (GASV(MT(K),MP(K),NVAR(K)+1,NW,NG) -                  &
    … …  
    213231                   (GASV(MT(K)+1,MP(K),NVAR(K)+1,NW,NG) -                &
    214232                   GASV(MT(K)+1,MP(K),NVAR(K),NW,NG))
     233
    215234           endif
    216235
    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) +           &
     236           ! Interpolate the gaseous k-coefficients to the requested T,P values
     237
     238           ANS = LKCOEF(K,1)*KCOEF(1) + LKCOEF(K,2)*KCOEF(2) +            &
    220239                LKCOEF(K,3)*KCOEF(3) + LKCOEF(K,4)*KCOEF(4)
    221240
    222            TAUGAS = U(k)*ANS
    223 
    224            !TAUGSURF(NW,NG) = TAUGSURF(NW,NG) + TAUGAS
     241           TAUGAS  = U(k)*ANS
     242
    225243           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
     244           DTAUKV(K,nw,ng) = TAUGAS &
     245                             + TRAYAER & ! TRAYAER includes all scattering contributions
     246                             + DCONT ! For parameterized continuum aborption
    228247
    229248        end do
    230249
    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
     250        ! Now fill in the "clear" part of the spectrum (NG = L_NGAUSS),
     251        ! which holds continuum opacity only
     252
     253        NG              = L_NGAUSS
    236254        DTAUKV(K,nw,ng) = TRAY(K,NW) + DCONT ! For parameterized continuum absorption
     255
    237256        do iaer=1,naerkind
    238            DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) + TAEROS(K,NW,IAER)
    239            !     &                           + DCONT          ! For parameterized continuum absorption
     257           DTAUKV(K,nw,ng) = DTAUKV(K,nw,ng) +  TAEROS(K,NW,IAER)
    240258        end do ! a bug was here!
    241259
    … …  
    248266  !     we need to calculate the scattering albedo and asymmetry factors
    249267
     268  do iaer=1,naerkind
    250269  DO NW=1,L_NSPECTV
    251      DO K=2,L_LEVELS
    252         do iaer=1,naerkind
     270     DO K=2,L_LEVELS   ! AS: shouldn't this be L_LEVELS+1 ? (see optci)
    253271           TAUAEROLK(K,NW,IAER) = TAUAERO(K,IAER) * QSVAER(K,NW,IAER)
    254         end do
    255272     ENDDO
    256273  ENDDO
    257 
     274  end do
    258275
    259276  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))
     277    DO NG=1,L_NGAUSS
     278     DO L=1,L_NLAYRAD-1
     279
     280        K              = 2*L+1
     281        DTAUV(L,nw,ng) = DTAUKV(K,NW,NG) + DTAUKV(K+1,NW,NG)
     282
     283        atemp = 0.
     284        btemp = TRAY(K,NW) + TRAY(K+1,NW)
     285        ctemp=0.9999*(TRAY(K,NW) + TRAY(K+1,NW))
    269286           do iaer=1,naerkind
    270287              atemp = atemp +                                     &
    271288                   GVAER(K,NW,IAER)   * TAUAEROLK(K,NW,IAER) +    &
    272289                   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)
     290              btemp = btemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
     291              ctemp = ctemp + TAUAEROLK(K,NW,IAER) + TAUAEROLK(K+1,NW,IAER)
    277292           end do
    278 
     293           WBARV(L,nw,ng) = ctemp / DTAUV(L,nw,ng)
    279294           COSBV(L,NW,NG) = atemp/btemp
    280            WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
    281 
    282         END DO
     295
     296      END DO ! L vertical loop
    283297
    284298        !     No vertical averaging on bottom layer
    … …  
    299313        WBARV(L,nw,ng) = ctemp/DTAUV(L,nw,ng)
    300314
    301      END DO                 ! NG gauss point loop
     315     END DO                 ! NG Gauss loop
    302316  END DO                    ! NW spectral loop
    303317
    304 
    305 
    306318  ! Total extinction optical depths
    307319
    308   DO NW=1,L_NSPECTV        
     320  DO NW=1,L_NSPECTV       
    309321     DO NG=1,L_NGAUSS       ! full gauss loop
    310322        TAUV(1,NW,NG)=0.0D0
    … …  
    321333
    322334
    323   RETURN
    324 END SUBROUTINE OPTCV
     335  return
     336
     337
     338end subroutine optcv

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

    r726 r873  
    6666      REAL*8 blamv(L_NSPECTV+1) ! these are needed by suaer.F90
    6767
     68      !! AS: introduced to avoid doing same computations again for continuum
     69      INTEGER, DIMENSION(:,:), ALLOCATABLE :: indi
     70      INTEGER, DIMENSION(:,:), ALLOCATABLE :: indv
     71
    6872      !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 
    6973      REAL*8, DIMENSION(:,:,:,:,:), ALLOCATABLE :: gasi, gasv

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

    r869 r873  
    5555
    5656      double precision testcont ! for continuum absorption initialisation
     57
     58      integer :: dummy
    5759
    5860!=======================================================================
    … …  
    596598
    597599            ! first do self-induced absorption
    598             call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.)
     600            dummy = -9999
     601            call interpolateH2H2(500.D+0,250.D+0,17500.D+0,testcont,.true.,dummy)
    599602            ! then cross-interactions with other gases
    600603            do jgas=1,ngasmx
    … …  
    602605                  call interpolateN2H2(592.D+0,278.15D+0,200000.D+0,10000.D+0,testcont,.true.)
    603606               elseif (jgas .eq. igas_He) then
    604                   call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.)
     607                  dummy = -9999
     608                  call interpolateH2He(500.D+0,250.D+0,200000.D+0,10000.D+0,testcont,.true.,dummy)
    605609               endif
    606610            enddo
Note: See TracChangeset for help on using the changeset viewer.