Changeset 878


Ignore:
Timestamp:
Feb 12, 2013, 12:26:32 PM (12 years ago)
Author:
jleconte
Message:

12/02/2013 == JL

  • Follows previous commit by Aymeric about bilinear interpolations:
    • Extended to all existing continua
    • generalized bilinearbig to work for various size inputs
  • because N2 and H2O continua databases are smaller, improvement around 15% for

an earth case.

Location:
trunk/LMDZ.GENERIC
Files:
9 edited

Legend:

Unmodified
Added
Removed

  • trunk/LMDZ.GENERIC/README

    r875 r878  
    901901- Added a keyword to enable ocean runoff in callphys.def (activerunoff)
    902902
     903== 12/02/2013 == JL
     904- Follows previous commit by Aymeric about bilinear interpolations:
     905   - Extended to all existing continua
     906   - generalized bilinearbig to work for various size inputs
     907- because N2 and H2O continua databases are smaller, improvement around 15% for
     908an earth case.
    903909 
    904910
     911

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

    r873 r878  
    1       subroutine bilinearbig(x_arr,y_arr,f2d_arr,x_in,y_in,f,a)
     1      subroutine bilinearbig(nX,nY,x_arr,y_arr,f2d_arr,x_in,y_in,f,ind)
    22
    33!     Necessary for interpolation of continuum data
    … …  
    66      implicit none
    77
    8       integer nX,nY,i,j,a,b
    9       parameter(nX=2428)
    10       parameter(nY=10)
     8      integer nX,nY,i,j,ind,b
    119
    1210      real*8 x_in,y_in,x1,x2,y1,y2
    … …  
    2826   !! ... and actually calculations only need to be done once
    2927   !! --> Case 1 : we have not calculated yet
    30    if ( a == -9999) then
     28   if ( ind == -9999) then
    3129      !1st check we're within the wavenumber range
    3230      if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
    3331         f=0.0D+0
    34          a=-1
     32         ind=-1
    3533      else
    3634        i=1
    … …  
    4038          i=i+1
    4139          x2=x_arr(i)
    42           a=i-1
     40          ind=i-1
    4341        end do
    4442      endif
    4543   !! --> case 2 : we already saw we are out of wavenumber range
    46    else if ( a == -1) then
     44   else if ( ind == -1) then
    4745      f=0.0D+0
    4846      return
    4947   !! --> case 3 : we already determined a -- so we just proceed
    5048   else
    51       x1=x_arr(a)
    52       x2=x_arr(a+1)
     49      x1=x_arr(ind)
     50      x2=x_arr(ind+1)
    5351   endif
    5452
    … …  
    7472      endif
    7573
    76       f11=f2d_arr(a,b)
    77       f21=f2d_arr(a+1,b)
    78       f12=f2d_arr(a,b+1)
    79       f22=f2d_arr(a+1,b+1)
     74      f11=f2d_arr(ind,b)
     75      f21=f2d_arr(ind+1,b)
     76      f12=f2d_arr(ind,b+1)
     77      f22=f2d_arr(ind+1,b+1)
    8078
    8179      call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)

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

    r873 r878  
    1      subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,a)
     1     subroutine interpolateH2H2(wn,temp,pres,abcoef,firstcall,ind)
    22
    33!==================================================================
    … …  
    5454      integer nres
    5555
    56       integer a
     56      integer ind
    5757
    5858      if(temp.gt.400)then
    … …  
    109109      endif
    110110
    111          call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
     111         call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
    112112
    113113         !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'

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

    r873 r878  
    1      subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,a)
     1     subroutine interpolateH2He(wn,temp,presH2,presHe,abcoef,firstcall,ind)
    22
    33!==================================================================
    … …  
    5656      integer nres
    5757
    58       integer a
     58      integer ind
    5959 
    6060      if(temp.gt.400)then
    … …  
    112112      endif
    113113
    114          call bilinearbig(wn_arr,temp_arr,abs_arr,wn,temp,abcoef,a)
     114         call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
    115115
    116116         !print*,'the absorption is ',abcoef,' cm^5 molecule^-2'

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

    r716 r878  
    1      subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall)
     1     subroutine interpolateH2Ocont_CKD(wn,temp,presS,presF,abcoef,firstcall,ind)
    22
    33!==================================================================
    … …  
    4040      double precision data_tmp(nT)
    4141
    42       integer k
     42      integer k,ind
    4343      logical firstcall
    4444
    … …  
    129129         print*,'   H2O pressure ',presS,' Pa'
    130130         print*,'   air pressure ',presF,' Pa'
     131         
     132      endif
    131133
    132          call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
    133          print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
     134      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS,ind)
     135!      print*,'the self absorption is ',abcoefS,' cm^2 molecule^-1'
    134136
    135          call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
    136          print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
     137      call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF,ind)
     138!      print*,'the foreign absorption is ',abcoefF,' cm^2 molecule^-1'
    137139
    138          print*,'We have ',amagatS,' amagats of H2O vapour'
    139          print*,'and ',amagatF,' amagats of air'
     140!      print*,'We have ',amagatS,' amagats of H2O vapour'
     141!      print*,'and ',amagatF,' amagats of air'
    140142
    141          abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
    142          abcoef = abcoef*(presS/(presF+presS))      ! take H2O mixing ratio into account
     143      abcoef = abcoefS*amagatS + abcoefF*amagatF ! Eq. (15) in Clough (1989)
     144      abcoef = abcoef*(presS/(presF+presS))      ! take H2O mixing ratio into account
    143145                                                    ! abs coeffs are given per molecule of H2O
    144146
    145          Nmolec = (presS+presF)/(kB*temp)           ! assume ideal gas
    146          print*,'Total number of molecules per m^3 is',Nmolec
     147      Nmolec = (presS+presF)/(kB*temp)           ! assume ideal gas
     148!      print*,'Total number of molecules per m^3 is',Nmolec
    147149
    148          abcoef = abcoef*Nmolec/(100.0**2)          ! convert to m^-1
    149          print*,'So the total absorption is ',abcoef,' m^-1'
    150          print*,'And optical depth / km : ',1000.0*abcoef
     150      abcoef = abcoef*Nmolec/(100.0**2)          ! convert to m^-1
     151!      print*,'So the total absorption is ',abcoef,' m^-1'
     152!      print*,'And optical depth / km : ',1000.0*abcoef
    151153
     154
     155      if(wn.gt.500 .and. wn.lt.1400)then
     156      elseif(wn.gt.2100 .and. wn.lt.3000)then
    152157      else
     158         abcoef = 0.0
     159      endif
    153160
    154          call bilinearH2Ocont(wn_arr,temp_arr,abs_arrS,wn,temp,abcoefS)
    155          call bilinearH2Ocont(wn_arr,temp_arr,abs_arrF,wn,temp,abcoefF)
     161      ! unlike for Rayleigh scattering, we do not currently weight by the BB function
     162      ! however our bands are normally thin, so this is no big deal.
    156163
    157          abcoef = abcoefS*amagatS + abcoefF*amagatF
    158          abcoef = abcoef*(presS/(presF+presS))
    159 
    160          Nmolec = (presS+presF)/(kB*temp)
    161          abcoef = abcoef*Nmolec/(100.0**2)
    162 
    163          if(wn.gt.500 .and. wn.lt.1400)then
    164          elseif(wn.gt.2100 .and. wn.lt.3000)then
    165          else
    166             abcoef = 0.0
    167          endif
    168 
    169          ! unlike for Rayleigh scattering, we do not currently weight by the BB function
    170          ! however our bands are normally thin, so this is no big deal.
    171 
    172       endif
    173164
    174165      return
    175166    end subroutine interpolateH2Ocont_CKD
    176167
    177 
    178 !-------------------------------------------------------------------------
    179       subroutine bilinearH2Ocont(x_arr,y_arr,f2d_arr,x_in,y_in,f)
    180 !     Necessary for interpolation of continuum data
    181 
    182       implicit none
    183 
    184       integer nX,nY,i,j,a,b
    185       parameter(nX=1001)
    186       parameter(nY=11)
    187 
    188       real*8 x_in,y_in,x,y,x1,x2,y1,y2
    189       real*8 f,f11,f12,f21,f22,fA,fB
    190       real*8 x_arr(nX)
    191       real*8 y_arr(nY)
    192       real*8 f2d_arr(nX,nY)
    193      
    194       integer strlen
    195       character*100 label
    196       label='subroutine bilinear'
    197 
    198       x=x_in
    199       y=y_in
    200 
    201 !     1st check we're within the wavenumber range
    202       if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
    203          f=0.0D+0
    204          return
    205       else
    206          
    207 !     in the x (wavenumber) direction 1st
    208          i=1
    209  10      if (i.lt.(nX+1)) then
    210             if (x_arr(i).gt.x) then
    211                x1=x_arr(i-1)
    212                x2=x_arr(i)
    213                a=i-1
    214                i=9999
    215             endif
    216             i=i+1
    217             goto 10
    218          endif
    219       endif
    220      
    221       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
    222          write(*,*) 'Warning from bilinearH2Ocont:'
    223          write(*,*) 'Outside continuum temperature range!'
    224          if(y.lt.y_arr(1))then
    225             y=y_arr(1)+0.01
    226          endif
    227          if(y.gt.y_arr(nY))then
    228             y=y_arr(nY)-0.01
    229          endif
    230       else
    231 
    232 !     in the y (temperature) direction 2nd
    233          j=1
    234  20      if (j.lt.(nY+1)) then
    235             if (y_arr(j).gt.y) then
    236                y1=y_arr(j-1)
    237                y2=y_arr(j)
    238                b=j-1
    239                j=9999
    240             endif
    241             j=j+1
    242             goto 20
    243          endif
    244       endif
    245      
    246       f11=f2d_arr(a,b)
    247       f21=f2d_arr(a+1,b)
    248       f12=f2d_arr(a,b+1)
    249       f22=f2d_arr(a+1,b+1)
    250      
    251 !     1st in x-direction
    252       fA=f11*(x2-x)/(x2-x1)+f21*(x-x1)/(x2-x1)
    253       fB=f12*(x2-x)/(x2-x1)+f22*(x-x1)/(x2-x1)
    254      
    255 !     then in y-direction
    256       f=fA*(y2-y)/(y2-y1)+fB*(y-y1)/(y2-y1)
    257 
    258       return
    259     end subroutine bilinearH2Ocont

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

    r716 r878  
    1 subroutine interpolateN2N2(wn,temp,pres,abcoef,firstcall)
     1subroutine interpolateN2N2(wn,temp,pres,abcoef,firstcall,ind)
    22
    33  !==================================================================
    … …  
    2121  double precision temp               ! temperature            (Kelvin)
    2222  double precision pres               ! pressure               (Pascals)
     23  integer :: ind
    2324
    2425  ! output
    … …  
    104105     print*,'   pressure ',pres,' Pa'
    105106
    106      call bilinearN2N2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
     107  endif
     108     call bilinearbig(nS,nT,wn_arr,temp_arr,abs_arr,wn,temp,abcoef,ind)
    107109
    108      print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
    109      print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
     110!     print*,'the absorption is ',abcoef,' cm^5 molecule^-2'
     111!     print*,'or ',abcoef*losch**2,' cm^-1 amagat^-2'
    110112
    111113     abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
     114!     abcoef=0.
    112115
    113      print*,'We have ',amagat,' amagats of N2'
    114      print*,'So the absorption is ',abcoef,' m^-1'
     116!     print*,'We have ',amagat,' amagats of N2'
     117!     print*,'So the absorption is ',abcoef,' m^-1'
    115118
    116   else
     119!     unlike for Rayleigh scattering, we do not currently weight by the BB function
     120!     however our bands are normally thin, so this is no big deal.
    117121
    118      call bilinearN2N2(wn_arr,temp_arr,abs_arr,wn,temp,abcoef)
    119      abcoef=abcoef*losch**2*100.0*amagat**2 ! convert to m^-1
    120 
    121      ! unlike for Rayleigh scattering, we do not currently weight by the BB function
    122      ! however our bands are normally thin, so this is no big deal.
    123 
    124   endif
    125122
    126123  return
    127124end subroutine interpolateN2N2
    128125
    129 
    130 !-------------------------------------------------------------------------
    131       subroutine bilinearN2N2(x_arr,y_arr,f2d_arr,x_in,y_in,f)
    132 
    133       implicit none
    134 
    135       integer nX,nY,i,j,a,b
    136       parameter(nX=582)
    137       parameter(nY=10)
    138 
    139       real*8 x_in,y_in,x,y,x1,x2,y1,y2
    140       real*8 f,f11,f12,f21,f22,fA,fB
    141       real*8 x_arr(nX)
    142       real*8 y_arr(nY)
    143       real*8 f2d_arr(nX,nY)
    144      
    145       integer strlen
    146       character*100 label
    147       label='subroutine bilinear'
    148 
    149       x=x_in
    150       y=y_in
    151 
    152 !     1st check we're within the wavenumber range
    153       if ((x.lt.x_arr(2)).or.(x.gt.x_arr(nX-2))) then
    154          f=0.0D+0
    155          return
    156       else
    157          
    158 !     in the x (wavenumber) direction 1st
    159          i=1
    160  10      if (i.lt.(nX+1)) then
    161             if (x_arr(i).gt.x) then
    162                x1=x_arr(i-1)
    163                x2=x_arr(i)
    164                a=i-1
    165                i=9999
    166             endif
    167             i=i+1
    168             goto 10
    169          endif
    170       endif
    171      
    172       if ((y.lt.y_arr(1)).or.(y.gt.y_arr(nY))) then
    173          write(*,*) 'Warning from bilinearN2N2:'
    174          write(*,*) 'Outside continuum temperature range!'
    175          if(y.lt.y_arr(1))then
    176             y=y_arr(1)+0.01
    177          endif
    178          if(y.gt.y_arr(nY))then
    179             y=y_arr(nY)-0.01
    180          endif
    181       else
    182 
    183 !     in the y (temperature) direction 2nd
    184          j=1
    185  20      if (j.lt.(nY+1)) then
    186             if (y_arr(j).gt.y) then
    187                y1=y_arr(j-1)
    188                y2=y_arr(j)
    189                b=j-1
    190                j=9999
    191             endif
    192             j=j+1
    193             goto 20
    194          endif
    195       endif
    196      
    197       f11=f2d_arr(a,b)
    198       f21=f2d_arr(a+1,b)
    199       f12=f2d_arr(a,b+1)
    200       f22=f2d_arr(a+1,b+1)
    201 
    202       call bilinear(f,f11,f21,f12,f22,x,x1,x2,y,y1,y2)
    203 
    204       return
    205     end subroutine bilinearN2N2

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

    r873 r878  
    9292  !! AS: to save time in computing continuum (see bilinearbig)
    9393  IF (.not.ALLOCATED(indi)) THEN
    94       ALLOCATE(indi(L_NSPECTI,ngasmx))
     94      ALLOCATE(indi(L_NSPECTI,ngasmx,ngasmx))
    9595      indi = -9999  ! this initial value means "to be calculated"
    9696  ENDIF
    … …  
    162162              if(igas.eq.igas_N2)then
    163163
    164                  call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
     164                 interm = indi(nw,igas,igas)
     165                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     166                 indi(nw,igas,igas) = interm
    165167
    166168              elseif(igas.eq.igas_H2)then
    167169
    168170                 ! first do self-induced absorption
    169                  interm = indi(nw,igas)
     171                 interm = indi(nw,igas,igas)
    170172                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    171                  indi(nw,igas) = interm
     173                 indi(nw,igas,igas) = interm
    172174
    173175                 ! then cross-interactions with other gases
    … …  
    176178                    dtempc  = 0.0
    177179                    if(jgas.eq.igas_N2)then
    178                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
     180                       interm = indi(nw,igas,jgas)
     181                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     182                       indi(nw,igas,jgas) = interm
    179183                    elseif(jgas.eq.igas_He)then
    180                        interm = indi(nw,jgas)
     184                       interm = indi(nw,igas,jgas)
    181185                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    182                        indi(nw,jgas) = interm
     186                       indi(nw,igas,jgas) = interm
    183187                    endif
    184188                    dtemp = dtemp + dtempc
    … …  
    191195                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
    192196                 else
    193                     call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     197                    interm = indi(nw,igas,igas)
     198                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
     199                    indi(nw,igas,igas) = interm
    194200                 endif
    195201

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

    r873 r878  
    8787  !! AS: to save time in computing continuum (see bilinearbig)
    8888  IF (.not.ALLOCATED(indv)) THEN
    89       ALLOCATE(indv(L_NSPECTV,ngasmx))
     89      ALLOCATE(indv(L_NSPECTV,ngasmx,ngasmx))
    9090      indv = -9999 ! this initial value means "to be calculated"
    9191  ENDIF
    … …  
    157157              if(igas.eq.igas_N2)then
    158158
    159                  !call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.)
     159                 interm = indv(nw,igas,igas)
     160!                 call interpolateN2N2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
     161                 indv(nw,igas,igas) = interm
    160162                 ! only goes to 500 cm^-1, so unless we're around a cold brown dwarf, this is irrelevant in the visible
    161163
    … …  
    163165
    164166                 ! first do self-induced absorption
    165                  interm = indv(nw,igas)
     167                 interm = indv(nw,igas,igas)
    166168                 call interpolateH2H2(wn_cont,T_cont,p_cont,dtemp,.false.,interm)
    167                  indv(nw,igas) = interm
     169                 indv(nw,igas,igas) = interm
    168170
    169171                 ! then cross-interactions with other gases
    … …  
    172174                    dtempc  = 0.0
    173175                    if(jgas.eq.igas_N2)then
    174                        call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.)
     176                       interm = indv(nw,igas,jgas)
     177                       call interpolateN2H2(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
     178                       indv(nw,igas,jgas) = interm
    175179                       ! should be irrelevant in the visible
    176180                    elseif(jgas.eq.igas_He)then
    177                        interm = indv(nw,jgas)
     181                       interm = indv(nw,igas,jgas)
    178182                       call interpolateH2He(wn_cont,T_cont,p_cross,p_cont,dtempc,.false.,interm)
    179                        indv(nw,jgas) = interm
     183                       indv(nw,igas,jgas) = interm
    180184                    endif
    181185                    dtemp = dtemp + dtempc
    … …  
    188192                    call interpolateH2Ocont_PPC(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
    189193                 else
    190                     call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.)
     194                    interm = indv(nw,igas,igas)
     195                    call interpolateH2Ocont_CKD(wn_cont,T_cont,p_cont,p_air,dtemp,.false.,interm)
     196                    indv(nw,igas,igas) = interm
    191197                 endif
    192198

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

    r873 r878  
    6767
    6868      !! AS: introduced to avoid doing same computations again for continuum
    69       INTEGER, DIMENSION(:,:), ALLOCATABLE :: indi
    70       INTEGER, DIMENSION(:,:), ALLOCATABLE :: indv
     69      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indi
     70      INTEGER, DIMENSION(:,:,:), ALLOCATABLE :: indv
    7171
    7272      !!! ALLOCATABLE STUFF SO THAT DIMENSIONS ARE READ in *.dat FILES -- AS 12/2011 
Note: See TracChangeset for help on using the changeset viewer.