Changeset 3968

Timestamp:

Nov 20, 2025, 4:23:54 PM (35 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

Location:

trunk/LMDZ.GENERIC

Files:

• changelog.txt (modified) (1 diff)
• libf/phystd/calc_rayleigh.F90 (modified) (9 diffs)

trunk/LMDZ.GENERIC/changelog.txt

@@ -2129 +2129 @@
 Minor fix in start2archive: the dimension for albedoS (albedo recasted on the
 dynamics scalar grid) shoudld be ip1jmp1.
+
+== 20/11/2025 == GM
+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).

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

@@ -102 +102 @@
         
         if ((BWNV(L_NSPECTV+1).gt.60000.).and.(strictboundrayleigh)) then
-          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."
+          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."
           call abort_physic(subname,message,1)
         elseif ((BWNV(L_NSPECTV+1).gt.60000.).and..not.(strictboundrayleigh)) then
@@ -160 +160 @@
                      T0(:) = 288.15
                      P0(:) = 1.01325e5
-                     ! ng -> valid range of the measurements : 0.1807 - 1.8172 um
-                     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
+                     if (wn .lt. 55331) then
+                       ! Sneep et al, 2005
+                       ! doi:10.1016/j.jqsrt.2004.07.025
+                       ! ng -> valid range of the measurements : 0.1807 - 1.8172 um
+                       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
+                     else
+                       ! Cuthbertson and Cuthbertson, 1920 (extrapolation)
+                       ! doi:10.1098/rspa.1920.0020
+                       ng(:) = 1. + (6914.45/(156.85 - (wn*1e-4)**2))*1e-5
+                     endif
                      Fk = 1.1364 + 25.3e-12*wn**2
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4 ! wn*100 -> cm-1 to m-1
+                     ! N=P/(kB*T) and muvar/1000 -> g/mol to kg/mol
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_N2)then
-                     ! Thalman et al, 2014
-                     ! doi:10.1016/j.jqsrt.2014.05.030
+                     ! Sneep et al, 2005
+                     ! doi:10.1016/j.jqsrt.2004.07.025
                      T0(:) = 288.15
                      P0(:) = 1.01325e5
                      if(wn.gt.21360)then !between 21360 and 39370 cm-1. We extrapolate above.
-                       ng(:) = 1. + (5677.465 + 318.81874e13/(14.4e9 - wn**2))*1e-8  !there is an error on the paper e12 -> e13
+                       ng(:) = 1. + (5677.465 + 318.81874e12/(14.4e9 - wn**2))*1e-8  !there is an error on the paper e12 -> e13
                      else !between 4860 and 21360 cm-1. We extrapolate below.
-                       ng(:) = 1. + (6498.2 + 307.4335e13/(14.4e9 - wn**2))*1e-8
-                     endif
-                     Fk = 1.034 + 3.17e-12*wn
-                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                       ng(:) = 1. + (6498.2 + 307.4335e12/(14.4e9 - wn**2))*1e-8
+                     endif
+                     Fk = 1.034 + 3.17e-12*wn**2
+                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_H2O)then
-                     ! Harvey et al, 1998
-                     ! doi:10.1063/1.556029
-                     ! doi:10.1364/AO.35.001566 for wn<4840 cm-1
-                     a0 =  0.244257733
-                     a1 = 9.74634476e-3
-                     a2 = -3.73234996e-3
-                     a3 = 2.68678472e-4
-                     a4 = 1.58920570e-3
-                     a5 = 2.45934259e-3
-                     a6 = 0.900704920
-                     a7 = -1.66626219e-2
-                     luv = 0.2292020
-                     lir = 5.432937
-                     Tr = 273.15
-                     rhor = 1000.
-                     T0(:) = tmid(:)
-                     P0(:) = pmid(:)*scalep
-                     lr = 0.589
-                     rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314463*T0(:)/(muvar(:)/1000.))
-                     rhos(:) = rho(:)/rhor
-                     ts(:) = T0(:)/Tr
-                     
-                     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(:)
-                     ! ng -> valid range of the measurements : 0.2 - 1.1 um
-                     ng(:) = sqrt(2.*b(:)+1.)/sqrt(1.-b(:))
+                     Fk = (6.+3.*3e-4)/(6.-7.*3e-4) ! delta=3e-4 Murphy 1977 doi:10.1063/1.434794
                      if(wn<4840.) then ! necessary to prevent a singularity at 3230 cm-1
+                       ! Ciddor, 1996
+                       ! doi:10.1364/AO.35.001566 for wn<4840 cm-1
+                       T0(:)=293.15
+                       P0(:)=1333.
                        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)
-                     endif
-                     Fk = (6.+3.*3e-4)/(6.-7.*3e-4) ! delta=3e-4 Murphy 1977 doi:10.1063/1.434794
-                     tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                       tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                       tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
+                     elseif(wn>50000.) then
+                       ! Barrell and Sears, 1939 (extrapolation)
+                       ! doi:10.1098/rsta.1939.0004
+                       T0(:)=273.15
+                       P0(:)=101325.
+                       ng(:) = 1. + (245.40+2.187*(1e4/wn)**(-2))*1e-6
+                       tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                       tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
+                     else
+                       ! Harvey et al, 1998
+                       ! doi:10.1063/1.556029
+                       ! ng -> valid range of the measurements : 0.2 - 1.1 um
+                       a0 =  0.244257733
+                       a1 = 9.74634476e-3
+                       a2 = -3.73234996e-3
+                       a3 = 2.68678472e-4
+                       a4 = 1.58920570e-3
+                       a5 = 2.45934259e-3
+                       a6 = 0.900704920
+                       a7 = -1.66626219e-2
+                       luv = 0.2292020
+                       lir = 5.432937
+                       Tr = 273.15
+                       rhor = 1000.
+                       T0(:) = tmid(:)
+                       P0(:) = pmid(:)*scalep
+                       lr = 0.589
+                       rho(:) = mass_frac(igas,:)*muvar(:)/massmol(igas)*P0(:)/(8.314463*T0(:)/(muvar(:)/1000.))
+                       rhos(:) = rho(:)/rhor
+                       ts(:) = T0(:)/Tr
+                       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(:)
+                       ng(:) = sqrt(2.*b(:)+1.)/sqrt(1.-b(:))
+                       tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                       tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
+                     endif
                  elseif(igas.eq.igas_H2)then
                      ! Peck and Hung, 1977
@@ -213 +238 @@
                      P0(:) = 1.01325e5
                      ! ng -> valid range of the measurements : 0.1680 - 1.6945 um
-                     ng(:) = 1. + (14895.6/(180.7 - (wn*1e-4)**2) + 4903.7/(92.-(wn*1e-4)**2))*1e-6
+                     if(wn<59534.) then
+                       ng(:) = 1. + (14895.6/(180.7 - (wn*1e-4)**2) + 4903.7/(92.-(wn*1e-4)**2))*1e-6
+                     else
+                       ng(:) = 1. + (23.79 + 12307.2/(109.832-(wn*1e-4)**2))*1e-6 ! extrapolation
+                     endif
                      Fk = (6.+3.*0.02)/(6.-7.*0.02) ! delta=0.02 Hansen 1974
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_He)then
                      ! Thalman et al, 2014
@@ -225 +255 @@
                      Fk = 1.
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_CH4)then
                      ! Sneep et al, 2005
@@ -234 +265 @@
                      Fk = 1. 
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_CO)then
                      ! Sneep et al, 2005
@@ -240 +272 @@
                      P0(:) = 1.01325e5
                      ! ng -> valid range of the measurements : 0.168 - 0.288 um
-                     ng(:) = 1. + 22851e-8 + 0.456e4/(71427.**2 - wn**2)
+                     if(wn<59809.) then
+                       ng(:) = 1. + 22851e-8 + 0.456e4/(71427.**2 - wn**2)
+                     else
+                       ng(:) = 1.00028476-2.01518666e-9*wn+1.88043553e-14*wn**2 ! extrapolation from previous data
+                     endif
                      Fk = 1.016
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_Ar)then
-                     ! Thalman et al, 2014
-                     ! doi:10.1016/j.jqsrt.2014.05.030
+                     ! Sneep et al, 2005
+                     ! doi:10.1016/j.jqsrt.2004.07.025
                      T0(:) = 288.15
                      P0(:) = 1.01325e5
@@ -252 +289 @@
                      Fk = 1.
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  elseif(igas.eq.igas_O2)then
-                     ! Thalman et al, 2014
-                     ! doi:10.1016/j.jqsrt.2014.05.030
+                     ! Sneep et al, 2005
+                     ! doi:10.1016/j.jqsrt.2004.07.025
                      T0(:) = 273.15
                      P0(:) = 1.01325e5
-                     ! ng -> valid range of the measurements : 0.303 - 2.0 um
-                     ng(:) = 1. + (20564.8 + 2.480899e13/(4.09e9 - wn**2))*1e-8
+                     if (wn .lt. 18315) then
+                       ! ng -> valid range of the measurements : > 0.546 um
+                       ng(:) = 1. + (21351.3 + 21.85670/(4.09e9 - wn**2))*1e-8
+                     elseif ((18315 .le. wn) .and. (wn .lt. 34722)) then
+                       ! ng -> valid range of the measurements : 0.288 - 0.546 um
+                       ng(:) = 1. + (20564.8 + 24.80899/(4.09e9 - wn**2))*1e-8
+                     elseif ((34722 .le. wn) .and. (wn .lt. 45248)) then
+                       ! ng -> valid range of the measurements : 0.288 - 0.221 um
+                       ng(:) = 1. + (22120.4 + 20.31876/(4.09e9 - wn**2))*1e-8
+                     else
+                       ! ng -> valid range of the measurements : < 0.221 um
+                       ng(:) = 1. + (23796.7 + 16.89884/(4.09e9 - wn**2))*1e-8
+                     endif
                      Fk = 1.09 + 1.385e-11*wn**2 + 1.488e-20*wn**4
                      tauvari(igas,:) = mass_frac(igas,:)*((ng(:)**2-1.)/(ng(:)**2+2.))**2 * Fk * (wn*100.)**4
+                     tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
                  else
                      print*,'No rayleigh scattering for ',trim(gnom(igas)),'. No data found.'
                  endif
-                 
-                 tauconsti(igas,:) = 24.*pi**3 *6.022141E+023 / (g*(muvar(:)/1000.)*(P0(:)/(1.380649E-23*T0(:)))**2)
+                
                  ! N=P/(kB*T)
                  ! pmid*scalep -> mbar to Pa
@@ -276 +325 @@
             enddo !ngasmx
 
-            call blackn(dble(wn*1e2),dble(tstellar),df)
+            call blackl(dble(wn),dble(tstellar),df)
             df=df*(BWNV(N+1)-BWNV(N))/Nfine
             tauwei=tauwei+df
@@ -282 +331 @@
          
          enddo !Nfine
+         ! We add a scalep because pressure in radiative transfer is in mbar 
          TAURAY(:,N)=tausum(:)*scalep/tauwei
-         ! scalep because tausum/tauwei is in SI units and dp on optcv is in hPa unit
 
       end do !L_NSPECTV