Changeset 3667


Ignore:
Timestamp:
Apr 21, 2020, 9:10:37 AM (4 years ago)
Author:
oboucher
Message:

Correction of a big big bug on refractive indiex and Mie calc in the StratAer? module

File:
1 edited

Legend:

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

    r3526 r3667  
    1 !
    2 ! $Id$
    3 !
    41SUBROUTINE MIECALC_AER(tau_strat, piz_strat, cg_strat, tau_strat_wave, tau_lw_abs_rrtm, paprs, debut)
    52
     
    157154  REAL ss_g(nbands_sw_rrtm+nbands_lw_rrtm+nwave_sw+nwave_lw)
    158155
     156  !-- fichier h2so4_0.75_300.00_hummel_1988_p_q.dat
     157  ! -- wavenumber (cm-1), wavelength (um), n_r, n_i
    159158  INTEGER, PARAMETER :: nb_lambda_h2so4=62
    160159  REAL, DIMENSION (nb_lambda_h2so4,4) :: ref_ind
    161   !-- fichier h2so4_0.75_300.00_hummel_1988_p_q.dat
    162   ! -- wavenumber (cm-1), wavelength (um), n_r, n_i
    163   DATA ref_ind /                                &
    164    200.000,   50.0000,   2.01000,   6.5000E-01, &
    165    250.000,   40.0000,   1.94000,   6.3000E-01, &
    166    285.714,   35.0000,   1.72000,   5.2000E-01, &
    167    333.333,   30.0000,   1.73000,   2.9000E-01, &
    168    358.423,   27.9000,   1.78000,   2.5000E-01, &
    169    400.000,   25.0000,   1.84000,   2.4000E-01, &
    170    444.444,   22.5000,   1.82000,   2.9000E-01, &
    171    469.484,   21.3000,   1.79000,   2.5000E-01, &
    172    500.000,   20.0000,   1.81000,   2.3000E-01, &
    173    540.541,   18.5000,   1.92700,   3.0200E-01, &
    174    555.556,   18.0000,   1.95000,   4.1000E-01, &
    175    581.395,   17.2000,   1.72400,   5.9000E-01, &
    176    609.756,   16.4000,   1.52000,   4.1400E-01, &
    177    666.667,   15.0000,   1.59000,   2.1100E-01, &
    178    675.676,   14.8000,   1.61000,   2.0500E-01, &
    179    714.286,   14.0000,   1.64000,   1.9500E-01, &
    180    769.231,   13.0000,   1.69000,   1.9500E-01, &
    181    800.000,   12.5000,   1.74000,   1.9800E-01, &
    182    869.565,   11.5000,   1.89000,   3.7400E-01, &
    183    909.091,   11.0000,   1.67000,   4.8500E-01, &
    184    944.198,   10.5910,   1.72000,   3.4000E-01, &
    185   1000.000,   10.0000,   1.89000,   4.5500E-01, &
    186   1020.408,    9.8000,   1.91000,   6.8000E-01, &
    187   1052.632,    9.5000,   1.67000,   7.5000E-01, &
    188   1086.957,    9.2000,   1.60000,   5.8600E-01, &
    189   1111.111,    9.0000,   1.65000,   6.3300E-01, &
    190   1149.425,    8.7000,   1.53000,   7.7200E-01, &
    191   1176.471,    8.5000,   1.37000,   7.5500E-01, &
    192   1219.512,    8.2000,   1.20000,   6.4500E-01, &
    193   1265.823,    7.9000,   1.14000,   4.8800E-01, &
    194   1388.889,    7.2000,   1.21000,   1.7600E-01, &
    195   1538.462,    6.5000,   1.37000,   1.2800E-01, &
    196   1612.903,    6.2000,   1.42400,   1.6500E-01, &
    197   1666.667,    6.0000,   1.42500,   1.9500E-01, &
    198   1818.182,    5.5000,   1.33700,   1.8300E-01, &
    199   2000.000,    5.0000,   1.36000,   1.2100E-01, &
    200   2222.222,    4.5000,   1.38500,   1.2000E-01, &
    201   2500.000,    4.0000,   1.39800,   1.2600E-01, &
    202   2666.667,    3.7500,   1.39600,   1.3100E-01, &
    203   2857.143,    3.5000,   1.37600,   1.5800E-01, &
    204   2948.113,    3.3920,   1.35200,   1.5900E-01, &
    205   3125.000,    3.2000,   1.31100,   1.3500E-01, &
    206   3333.333,    3.0000,   1.29300,   9.5500E-02, &
    207   3703.704,    2.7000,   1.30300,   5.7000E-03, &
    208   4000.000,    2.5000,   1.34400,   3.7600E-03, &
    209   4444.444,    2.2500,   1.37000,   1.8000E-03, &
    210   5000.000,    2.0000,   1.38400,   1.2600E-03, &
    211   5555.556,    1.8000,   1.39000,   5.5000E-04, &
    212   6510.417,    1.5360,   1.40300,   1.3700E-04, &
    213   7692.308,    1.3000,   1.41000,   1.0000E-05, &
    214   9433.962,    1.0600,   1.42000,   1.5000E-06, &
    215  11627.907,    0.8600,   1.42500,   1.7900E-07, &
    216  14409.222,    0.6940,   1.42800,   1.9900E-08, &
    217  15797.788,    0.6330,   1.42900,   1.4700E-08, &
    218  18181.818,    0.5500,   1.43000,   1.0000E-08, &
    219  19417.476,    0.5150,   1.43100,   1.0000E-08, &
    220  20491.803,    0.4880,   1.43200,   1.0000E-08, &
    221  25000.000,    0.4000,   1.44000,   1.0000E-08, &
    222  29673.591,    0.3370,   1.45900,   1.0000E-08, &
    223  33333.333,    0.3000,   1.46900,   1.0000E-08, &
    224  40000.000,    0.2500,   1.48400,   1.0000E-08, &
    225  50000.000,    0.2000,   1.49800,   1.0000E-08 /
     160
    226161!---------------------------------------------------------
    227162
    228163  IF (debut) THEN   
    229164
     165     ref_ind = RESHAPE( (/ &
     166      200.000,   50.0000,   2.01000,   6.5000E-01, &
     167      250.000,   40.0000,   1.94000,   6.3000E-01, &
     168      285.714,   35.0000,   1.72000,   5.2000E-01, &
     169      333.333,   30.0000,   1.73000,   2.9000E-01, &
     170      358.423,   27.9000,   1.78000,   2.5000E-01, &
     171      400.000,   25.0000,   1.84000,   2.4000E-01, &
     172      444.444,   22.5000,   1.82000,   2.9000E-01, &
     173      469.484,   21.3000,   1.79000,   2.5000E-01, &
     174      500.000,   20.0000,   1.81000,   2.3000E-01, &
     175      540.541,   18.5000,   1.92700,   3.0200E-01, &
     176      555.556,   18.0000,   1.95000,   4.1000E-01, &
     177      581.395,   17.2000,   1.72400,   5.9000E-01, &
     178      609.756,   16.4000,   1.52000,   4.1400E-01, &
     179      666.667,   15.0000,   1.59000,   2.1100E-01, &
     180      675.676,   14.8000,   1.61000,   2.0500E-01, &
     181      714.286,   14.0000,   1.64000,   1.9500E-01, &
     182      769.231,   13.0000,   1.69000,   1.9500E-01, &
     183      800.000,   12.5000,   1.74000,   1.9800E-01, &
     184      869.565,   11.5000,   1.89000,   3.7400E-01, &
     185      909.091,   11.0000,   1.67000,   4.8500E-01, &
     186      944.198,   10.5910,   1.72000,   3.4000E-01, &
     187     1000.000,   10.0000,   1.89000,   4.5500E-01, &
     188     1020.408,    9.8000,   1.91000,   6.8000E-01, &
     189     1052.632,    9.5000,   1.67000,   7.5000E-01, &
     190     1086.957,    9.2000,   1.60000,   5.8600E-01, &
     191     1111.111,    9.0000,   1.65000,   6.3300E-01, &
     192     1149.425,    8.7000,   1.53000,   7.7200E-01, &
     193     1176.471,    8.5000,   1.37000,   7.5500E-01, &
     194     1219.512,    8.2000,   1.20000,   6.4500E-01, &
     195     1265.823,    7.9000,   1.14000,   4.8800E-01, &
     196     1388.889,    7.2000,   1.21000,   1.7600E-01, &
     197     1538.462,    6.5000,   1.37000,   1.2800E-01, &
     198     1612.903,    6.2000,   1.42400,   1.6500E-01, &
     199     1666.667,    6.0000,   1.42500,   1.9500E-01, &
     200     1818.182,    5.5000,   1.33700,   1.8300E-01, &
     201     2000.000,    5.0000,   1.36000,   1.2100E-01, &
     202     2222.222,    4.5000,   1.38500,   1.2000E-01, &
     203     2500.000,    4.0000,   1.39800,   1.2600E-01, &
     204     2666.667,    3.7500,   1.39600,   1.3100E-01, &
     205     2857.143,    3.5000,   1.37600,   1.5800E-01, &
     206     2948.113,    3.3920,   1.35200,   1.5900E-01, &
     207     3125.000,    3.2000,   1.31100,   1.3500E-01, &
     208     3333.333,    3.0000,   1.29300,   9.5500E-02, &
     209     3703.704,    2.7000,   1.30300,   5.7000E-03, &
     210     4000.000,    2.5000,   1.34400,   3.7600E-03, &
     211     4444.444,    2.2500,   1.37000,   1.8000E-03, &
     212     5000.000,    2.0000,   1.38400,   1.2600E-03, &
     213     5555.556,    1.8000,   1.39000,   5.5000E-04, &
     214     6510.417,    1.5360,   1.40300,   1.3700E-04, &
     215     7692.308,    1.3000,   1.41000,   1.0000E-05, &
     216     9433.962,    1.0600,   1.42000,   1.5000E-06, &
     217    11627.907,    0.8600,   1.42500,   1.7900E-07, &
     218    14409.222,    0.6940,   1.42800,   1.9900E-08, &
     219    15797.788,    0.6330,   1.42900,   1.4700E-08, &
     220    18181.818,    0.5500,   1.43000,   1.0000E-08, &
     221    19417.476,    0.5150,   1.43100,   1.0000E-08, &
     222    20491.803,    0.4880,   1.43200,   1.0000E-08, &
     223    25000.000,    0.4000,   1.44000,   1.0000E-08, &
     224    29673.591,    0.3370,   1.45900,   1.0000E-08, &
     225    33333.333,    0.3000,   1.46900,   1.0000E-08, &
     226    40000.000,    0.2500,   1.48400,   1.0000E-08, &
     227    50000.000,    0.2000,   1.49800,   1.0000E-08 /), (/nb_lambda_h2so4,4/), order=(/2,1/) )
     228
    230229  !--initialising dry diameters to geometrically spaced mass/volume (see Jacobson 1994)
    231       mdw(1)=mdwmin
    232       IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
    233         mdw(2)=mdw(1)*2.**(1./3.)
    234         DO it=3, nbtr_bin
    235           mdw(it)=mdw(it-1)*V_rat**(1./3.)
    236         ENDDO
    237       ELSE
    238         DO it=2, nbtr_bin
    239           mdw(it)=mdw(it-1)*V_rat**(1./3.)
    240         ENDDO
    241       ENDIF
    242       PRINT *,'init mdw=', mdw
     230    mdw(1)=mdwmin
     231    IF (V_rat.LT.1.62) THEN ! compensate for dip in second bin for lower volume ratio
     232      mdw(2)=mdw(1)*2.**(1./3.)
     233      DO it=3, nbtr_bin
     234        mdw(it)=mdw(it-1)*V_rat**(1./3.)
     235      ENDDO
     236    ELSE
     237      DO it=2, nbtr_bin
     238        mdw(it)=mdw(it-1)*V_rat**(1./3.)
     239      ENDDO
     240    ENDIF
     241    WRITE(lunout,*) 'init mdw=', mdw
    243242
    244243    !--compute particle radius for a composition of 75% H2SO4 / 25% H2O at T=293K
     
    286285
    287286      IF (refr_ind_interpol) THEN
    288 
     287 
    289288        ilambda_max=ref_ind(1,2)/1.e6 !--in m
    290289        ilambda_min=ref_ind(nb_lambda_h2so4,2)/1.e6 !--in m
     
    381380        Nmax=INT(x+4*x**(1./3.)+2.)+1
    382381      ELSE
    383         PRINT *, 'x out of bound, x=', x
     382        WRITE(lunout,*) 'x out of bound, x=', x
    384383        STOP
    385384      ENDIF
     
    461460      omegatot=omegatot+r**2*Q_ext*omega*number
    462461      gtot    =gtot+r**2*Q_sca*g*number
    463 
     462 
    464463      ENDDO   !---bin
     464
    465465    !------------------------------------------------------------------
    466466
Note: See TracChangeset for help on using the changeset viewer.