Changeset 2948 for LMDZ5/trunk/libf/phylmd/StratAer

Timestamp:

Jul 13, 2017, 9:15:00 PM (7 years ago)

Author:

oboucher

Message:

Nstart is increased for the downward recurrence for more accurate results.
The calculations for the very small x (small particles and large wavelengths) has accuracy problem

File:

• LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90 (modified) (6 diffs)

LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90

@@ -60 +60 @@
   COMPLEX nn
   REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
-  REAL x  !--size parameter
+  REAL x, x_old  !--size parameter
   REAL r, r_lower, r_upper  !--radius
   REAL sigma_sca, sigma_ext, sigma_abs
@@ -70 +70 @@
   INTEGER bin, Nbin, it
   PARAMETER (Nbin=10)
-
+  LOGICAL smallx
 
 !---wavelengths STREAMER
@@ -347 +347 @@
       r_lower=exp(log(rmin)+FLOAT(bin-1)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
       r_upper=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
+      deltar=r_upper-r_lower
+
       r=sqrt(r_lower*r_upper)
       x=2.*RPI*r/lambda_int(Nwv)
-      deltar=r_upper-r_lower
+
+!we impose a minimum value for x and extrapolate quantities for small x values
+      smallx = .FALSE.
+      IF (x.LT.0.001) THEN
+        smallx = .TRUE.
+        x_old = x
+        x = 0.001
+      ENDIF
 
       number=Ntot*deltar/(rmax-rmin) !dN/dr constant over tracer bin
@@ -373 +382 @@
       ENDIF
 
-      Nstart=Nmax+10
+      Nstart=Nmax+100
 
     !-----------loop for nu1z1, nu1z2
@@ -417 +426 @@
       Q_sca=0.0
       g=0.0
+
       DO n=Nmax-1,1,-1
         nnn=FLOAT(n)
@@ -426 +436 @@
               (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
       ENDDO
+
       Q_ext=2./x**2 * Q_ext
       Q_sca=2./x**2 * Q_sca
+    !--extrapolation in case of small x values
+      IF (smallx) THEN
+        Q_ext = x_old/x * Q_ext
+        Q_sca = x_old/x * Q_sca
+      ENDIF
+
       Q_abs=Q_ext-Q_sca
+
       IF (AIMAG(m).EQ.0.0) Q_abs=0.0
       omega=Q_sca/Q_ext
+
+    ! g is wrong in the smallx case (but that does not matter as long as we ignore LW scattering)
       g=g*4./x**2/Q_sca