Ignore:
Timestamp:
Nov 20, 2025, 4:23:54 PM (41 hours ago)
Author:
gmilcareck
Message:

Generic PCM:
Fix bug and incorrect data and cleanup calc_rayleigh.F90.

  • [minor] Fix a very minor bug for reference pressure and temperature.

The same reference temperature and pressure was used for all gases
instead of use the one for each gas. But no impact on simulations is expected
because all gases have almost the same reference T/P.

  • [major] Data for H2O and O2 refractive index was incorrect on the reference paper. It is now fixed.
  • To remove singularities at wavenumber > 60000 cm-1 on some gases, we add extrapolated formula of refractive index at these wavenumbers. However, we keep the safety check strictboundrayleigh to remember that data is extrapolated after 60000 cm-1. So now, you can use rayleigh scattering at wavenumber > 60000 cm-1 but keep on mind that this is extrapolated data. Nobody knows what is going on at these wavenumbers (but Nobody is lost at sea with his crew).

GM

File:
1 edited

Legend:

Unmodified
Added
Removed

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

    r3863 r3968  
    102102       
    103103        if ((BWNV(L_NSPECTV+1).gt.60000.).and.(strictboundrayleigh)) then
    104           message="Rayleigh scattering is unknown for wn>60000 cm-1 - some singularities are present for higher wavenumber - if you know what you are doing, use strictboundrayleigh=.false."
     104          message="Rayleigh scattering is unknown for wn>60000 cm-1 - all data is extrapolated for higher wavenumber - if you know what you are doing, use strictboundrayleigh=.false."
    105105          call abort_physic(subname,message,1)
    106106        elseif ((BWNV(L_NSPECTV+1).gt.60000.).and..not.(strictboundrayleigh)) then
    … …  
    160160                     T0(:) = 288.15
    161161                     P0(:) = 1.01325e5
    162                      ! ng -> valid range of the measurements : 0.1807 - 1.8172 um
    163                      ng(:) = 1. + 1.1427e3*(5799.25/(128908.9**2 - wn**2) + 120.05/(89223.8**2 - wn**2) + 5.3334/(75037.5**2 - wn**2) + 4.3244/(67837.7**2 - wn**2) + 0.1218145e-4/(2418.136**2 - wn**2)) ! there is an error on the paper 1.1427e6 -> 1.1427e3
     162                     if (wn .lt. 55331) then
     163                       ! Sneep et al, 2005
     164                       ! doi:10.1016/j.jqsrt.2004.07.025
     165                       ! ng -> valid range of the measurements : 0.1807 - 1.8172 um
     166                       ng(:) = 1. + 1.1427e3*(5799.25/(128908.9**2 - wn**2) + 120.05/(89223.8**2 - wn**2) + 5.3334/(75037.5**2 - wn**2) + 4.3244/(67837.7**2 - wn**2) + 0.1218145e-4/(2418.136**2 - wn**2)) ! there is an error on the paper 1.1427e6 -> 1.1427e3
     167                     else
     168                       ! Cuthbertson and Cuthbertson, 1920 (extrapolation)
     169                       ! doi:10.1098/rspa.1920.0020
     170                       ng(:) = 1. + (6914.45/(156.85 - (wn*1e-4)**2))*1e-5
     171                     endif
    164172                     Fk = 1.1364 + 25.3e-12*wn**2
    165173                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4 ! wn*100 -> cm-1 to m-1
     174                     ! N=P/(kB*T) and muvar/1000 -> g/mol to kg/mol
     175                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    166176                 elseif(igas.eq.igas_N2)then
    167                      ! Thalman et al, 2014
    168                      ! doi:10.1016/j.jqsrt.2014.05.030
     177                     ! Sneep et al, 2005
     178                     ! doi:10.1016/j.jqsrt.2004.07.025
    169179                     T0(:) = 288.15
    170180                     P0(:) = 1.01325e5
    171181                     if(wn.gt.21360)then !between 21360 and 39370 cm-1. We extrapolate above.
    172                        ng(:) = 1. + (5677.465 + 318.81874e13/(14.4e9 - wn**2))*1e-8  !there is an error on the paper e12 -> e13
     182                       ng(:) = 1. + (5677.465 + 318.81874e12/(14.4e9 - wn**2))*1e-8  !there is an error on the paper e12 -> e13
    173183                     else !between 4860 and 21360 cm-1. We extrapolate below.
    174                        ng(:) = 1. + (6498.2 + 307.4335e13/(14.4e9 - wn**2))*1e-8
    175                      endif
    176                      Fk = 1.034 + 3.17e-12*wn
    177                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     184                       ng(:) = 1. + (6498.2 + 307.4335e12/(14.4e9 - wn**2))*1e-8
     185                     endif
     186                     Fk = 1.034 + 3.17e-12*wn**2
     187                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     188                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    178189                 elseif(igas.eq.igas_H2O)then
    179                      ! Harvey et al, 1998
    180                      ! doi:10.1063/1.556029
    181                      ! doi:10.1364/AO.35.001566 for wn<4840 cm-1
    182                      a0 =  0.244257733
    183                      a1 = 9.74634476e-3
    184                      a2 = -3.73234996e-3
    185                      a3 = 2.68678472e-4
    186                      a4 = 1.58920570e-3
    187                      a5 = 2.45934259e-3
    188                      a6 = 0.900704920
    189                      a7 = -1.66626219e-2
    190                      luv = 0.2292020
    191                      lir = 5.432937
    192                      Tr = 273.15
    193                      rhor = 1000.
    194                      T0(:) = tmid(:)
    195                      P0(:) = pmid(:)*scalep
    196                      lr = 0.589
    197                      rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314463*T0(:)/(muvar(:)/1000.))
    198                      rhos(:) = rho(:)/rhor
    199                      ts(:) = T0(:)/Tr
    200                      
    201                      b(:) = (a0 + a1*rhos(:) + a2*ts(:) + a3*ts(:)*(10000./wn/lr)**2 + a4/(10000./wn/lr)**2 + a5/((10000./wn/lr)**2 - luv**2) + a6/((10000./wn/lr)**2 - lir**2) + a7*rhos(:)**2)*rhos(:)
    202                      ! ng -> valid range of the measurements : 0.2 - 1.1 um
    203                      ng(:) = sqrt(2.*b(:)+1.)/sqrt(1.-b(:))
     190                     Fk = (6.+3.*3e-4)/(6.-7.*3e-4) ! delta=3e-4 Murphy 1977 doi:10.1063/1.434794
    204191                     if(wn<4840.) then ! necessary to prevent a singularity at 3230 cm-1
     192                       ! Ciddor, 1996
     193                       ! doi:10.1364/AO.35.001566 for wn<4840 cm-1
     194                       T0(:)=293.15
     195                       P0(:)=1333.
    205196                       ng(:) = 1. + 1.022e-8*(295.235 + 2.6422*(wn*1e-4)**2 - 0.032380*(wn*1e-4)**4 + 0.004028*(wn*1e-4)**6)
    206                      endif
    207                      Fk = (6.+3.*3e-4)/(6.-7.*3e-4) ! delta=3e-4 Murphy 1977 doi:10.1063/1.434794
    208                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     197                       tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     198                       tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
     199                     elseif(wn>50000.) then
     200                       ! Barrell and Sears, 1939 (extrapolation)
     201                       ! doi:10.1098/rsta.1939.0004
     202                       T0(:)=273.15
     203                       P0(:)=101325.
     204                       ng(:) = 1. + (245.40+2.187*(1e4/wn)**(-2))*1e-6
     205                       tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     206                       tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
     207                     else
     208                       ! Harvey et al, 1998
     209                       ! doi:10.1063/1.556029
     210                       ! ng -> valid range of the measurements : 0.2 - 1.1 um
     211                       a0 =  0.244257733
     212                       a1 = 9.74634476e-3
     213                       a2 = -3.73234996e-3
     214                       a3 = 2.68678472e-4
     215                       a4 = 1.58920570e-3
     216                       a5 = 2.45934259e-3
     217                       a6 = 0.900704920
     218                       a7 = -1.66626219e-2
     219                       luv = 0.2292020
     220                       lir = 5.432937
     221                       Tr = 273.15
     222                       rhor = 1000.
     223                       T0(:) = tmid(:)
     224                       P0(:) = pmid(:)*scalep
     225                       lr = 0.589
     226                       rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314463*T0(:)/(muvar(:)/1000.))
     227                       rhos(:) = rho(:)/rhor
     228                       ts(:) = T0(:)/Tr
     229                       b(:) = (a0 + a1*rhos(:) + a2*ts(:) + a3*ts(:)*(10000./wn/lr)**2 + a4/(10000./wn/lr)**2 + a5/((10000./wn/lr)**2 - luv**2) + a6/((10000./wn/lr)**2 - lir**2) + a7*rhos(:)**2)*rhos(:)
     230                       ng(:) = sqrt(2.*b(:)+1.)/sqrt(1.-b(:))
     231                       tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     232                       tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
     233                     endif
    209234                 elseif(igas.eq.igas_H2)then
    210235                     ! Peck and Hung, 1977
    … …  
    213238                     P0(:) = 1.01325e5
    214239                     ! ng -> valid range of the measurements : 0.1680 - 1.6945 um
    215                      ng(:) = 1. + (14895.6/(180.7 - (wn*1e-4)**2) + 4903.7/(92.-(wn*1e-4)**2))*1e-6
     240                     if(wn<59534.) then
     241                       ng(:) = 1. + (14895.6/(180.7 - (wn*1e-4)**2) + 4903.7/(92.-(wn*1e-4)**2))*1e-6
     242                     else
     243                       ng(:) = 1. + (23.79 + 12307.2/(109.832-(wn*1e-4)**2))*1e-6 ! extrapolation
     244                     endif
    216245                     Fk = (6.+3.*0.02)/(6.-7.*0.02) ! delta=0.02 Hansen 1974
    217246                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     247                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    218248                 elseif(igas.eq.igas_He)then
    219249                     ! Thalman et al, 2014
    … …  
    225255                     Fk = 1.
    226256                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     257                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    227258                 elseif(igas.eq.igas_CH4)then
    228259                     ! Sneep et al, 2005
    … …  
    234265                     Fk = 1.
    235266                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     267                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    236268                 elseif(igas.eq.igas_CO)then
    237269                     ! Sneep et al, 2005
    … …  
    240272                     P0(:) = 1.01325e5
    241273                     ! ng -> valid range of the measurements : 0.168 - 0.288 um
    242                      ng(:) = 1. + 22851e-8 + 0.456e4/(71427.**2 - wn**2)
     274                     if(wn<59809.) then
     275                       ng(:) = 1. + 22851e-8 + 0.456e4/(71427.**2 - wn**2)
     276                     else
     277                       ng(:) = 1.00028476-2.01518666e-9*wn+1.88043553e-14*wn**2 ! extrapolation from previous data
     278                     endif
    243279                     Fk = 1.016
    244280                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     281                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    245282                 elseif(igas.eq.igas_Ar)then
    246                      ! Thalman et al, 2014
    247                      ! doi:10.1016/j.jqsrt.2014.05.030
     283                     ! Sneep et al, 2005
     284                     ! doi:10.1016/j.jqsrt.2004.07.025
    248285                     T0(:) = 288.15
    249286                     P0(:) = 1.01325e5
    … …  
    252289                     Fk = 1.
    253290                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     291                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    254292                 elseif(igas.eq.igas_O2)then
    255                      ! Thalman et al, 2014
    256                      ! doi:10.1016/j.jqsrt.2014.05.030
     293                     ! Sneep et al, 2005
     294                     ! doi:10.1016/j.jqsrt.2004.07.025
    257295                     T0(:) = 273.15
    258296                     P0(:) = 1.01325e5
    259                      ! ng -> valid range of the measurements : 0.303 - 2.0 um
    260                      ng(:) = 1. + (20564.8 + 2.480899e13/(4.09e9 - wn**2))*1e-8
     297                     if (wn .lt. 18315) then
     298                       ! ng -> valid range of the measurements : > 0.546 um
     299                       ng(:) = 1. + (21351.3 + 21.85670/(4.09e9 - wn**2))*1e-8
     300                     elseif ((18315 .le. wn) .and. (wn .lt. 34722)) then
     301                       ! ng -> valid range of the measurements : 0.288 - 0.546 um
     302                       ng(:) = 1. + (20564.8 + 24.80899/(4.09e9 - wn**2))*1e-8
     303                     elseif ((34722 .le. wn) .and. (wn .lt. 45248)) then
     304                       ! ng -> valid range of the measurements : 0.288 - 0.221 um
     305                       ng(:) = 1. + (22120.4 + 20.31876/(4.09e9 - wn**2))*1e-8
     306                     else
     307                       ! ng -> valid range of the measurements : < 0.221 um
     308                       ng(:) = 1. + (23796.7 + 16.89884/(4.09e9 - wn**2))*1e-8
     309                     endif
    261310                     Fk = 1.09 + 1.385e-11*wn**2 + 1.488e-20*wn**4
    262311                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
     312                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
    263313                 else
    264314                     print*,'No rayleigh scattering for ',trim(gnom(igas)),'. No data found.'
    265315                 endif
    266                  
    267                  tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
     316               
    268317                 ! N=P/(kB*T)
    269318                 ! pmid*scalep -> mbar to Pa
    … …  
    276325            enddo !ngasmx
    277326
    278             call blackn(dble(wn*1e2),dble(tstellar),df)
     327            call blackl(dble(wn),dble(tstellar),df)
    279328            df=df*(BWNV(N+1)-BWNV(N))/Nfine
    280329            tauwei=tauwei+df
    … …  
    282331         
    283332         enddo !Nfine
     333         ! We add a scalep because pressure in radiative transfer is in mbar
    284334         TAURAY(:,N)=tausum(:)*scalep/tauwei
    285          ! scalep because tausum/tauwei is in SI units and dp on optcv is in hPa unit
    286335
    287336      end do !L_NSPECTV
Note: See TracChangeset for help on using the changeset viewer.