Ignore:
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:
1 edited

Legend:

Unmodified
Added
Removed
  • LMDZ5/trunk/libf/phylmd/StratAer/miecalc_aer.F90

    r2704 r2948  
    6060  COMPLEX nn
    6161  REAL Q_ext, Q_abs, Q_sca, g, omega   !--parameters for radius r
    62   REAL x  !--size parameter
     62  REAL x, x_old  !--size parameter
    6363  REAL r, r_lower, r_upper  !--radius
    6464  REAL sigma_sca, sigma_ext, sigma_abs
     
    7070  INTEGER bin, Nbin, it
    7171  PARAMETER (Nbin=10)
    72 
     72  LOGICAL smallx
    7373
    7474!---wavelengths STREAMER
     
    347347      r_lower=exp(log(rmin)+FLOAT(bin-1)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
    348348      r_upper=exp(log(rmin)+FLOAT(bin)/FLOAT(Nbin)*(log(rmax)-log(rmin)))
     349      deltar=r_upper-r_lower
     350
    349351      r=sqrt(r_lower*r_upper)
    350352      x=2.*RPI*r/lambda_int(Nwv)
    351       deltar=r_upper-r_lower
     353
     354!we impose a minimum value for x and extrapolate quantities for small x values
     355      smallx = .FALSE.
     356      IF (x.LT.0.001) THEN
     357        smallx = .TRUE.
     358        x_old = x
     359        x = 0.001
     360      ENDIF
    352361
    353362      number=Ntot*deltar/(rmax-rmin) !dN/dr constant over tracer bin
     
    373382      ENDIF
    374383
    375       Nstart=Nmax+10
     384      Nstart=Nmax+100
    376385
    377386    !-----------loop for nu1z1, nu1z2
     
    417426      Q_sca=0.0
    418427      g=0.0
     428
    419429      DO n=Nmax-1,1,-1
    420430        nnn=FLOAT(n)
     
    426436              (2.*nnn+1.)/nnn/(nnn+1.) * REAL(a(n)*CONJG(b(n)))
    427437      ENDDO
     438
    428439      Q_ext=2./x**2 * Q_ext
    429440      Q_sca=2./x**2 * Q_sca
     441    !--extrapolation in case of small x values
     442      IF (smallx) THEN
     443        Q_ext = x_old/x * Q_ext
     444        Q_sca = x_old/x * Q_sca
     445      ENDIF
     446
    430447      Q_abs=Q_ext-Q_sca
     448
    431449      IF (AIMAG(m).EQ.0.0) Q_abs=0.0
    432450      omega=Q_sca/Q_ext
     451
     452    ! g is wrong in the smallx case (but that does not matter as long as we ignore LW scattering)
    433453      g=g*4./x**2/Q_sca
    434454
Note: See TracChangeset for help on using the changeset viewer.