source: LMDZ6/branches/Amaury_dev/libf/phylmd/Dust/splaeropt_6bands_rrtm.f90 @ 5202

Last change on this file since 5202 was 5137, checked in by abarral, 4 months ago

Put gradsdef.h, tracstoke.h, clesphys.h into modules
File size: 11.2 KB
Line 
1! $Id: splaeropt_6bands_rrtm.F90 2644 2016-10-02 16:55:08Z oboucher $
2
3SUBROUTINE SPLAEROPT_6BANDS_RRTM(&
4        zdm, tr_seri, RHcl, &
5        tau_allaer, piz_allaer, cg_allaer)
6
7  USE dimphy
8  USE aero_mod
9  USE infotrac_phy, ONLY: nqtot, nbtr, tracers
10  USE phys_local_var_mod, ONLY: abs550aer
11  USE lmdz_abort_physic, ONLY: abort_physic
12  USE lmdz_clesphys
13
14  ! Olivier Boucher Jan 2017
15  ! based on Mie routines on ciclad CMIP6
16
17  IMPLICIT NONE
18
19  ! Input arguments:
20
21  REAL, DIMENSION(klon, klev), INTENT(IN) :: zdm        !--hauteur des couches en kg/m2
22  REAL, DIMENSION(klon, klev), INTENT(IN) :: RHcl       ! humidite relative ciel clair
23  REAL, DIMENSION(klon, klev, nbtr), INTENT(IN) :: tr_seri
24
25  ! Output arguments:
26  ! 2= total aerosols
27  ! 1= natural aerosols
28
29  REAL, DIMENSION(klon, klev, 2, nbands_sw_rrtm), INTENT(OUT) :: tau_allaer ! epaisseur optique aerosol
30  REAL, DIMENSION(klon, klev, 2, nbands_sw_rrtm), INTENT(OUT) :: piz_allaer ! single scattering albedo aerosol
31  REAL, DIMENSION(klon, klev, 2, nbands_sw_rrtm), INTENT(OUT) :: cg_allaer  ! asymmetry parameter aerosol
32
33  ! Local
34
35  LOGICAL :: soluble
36  INTEGER :: i, k, irh, iq, itr, inu
37  INTEGER :: aerindex, spsol, spinsol
38  INTEGER :: RH_num(klon, klev)
39
40  INTEGER, PARAMETER :: naero_soluble = 2    ! 1- accumulation soluble; 2- coarse soluble
41  INTEGER, PARAMETER :: naero_insoluble = 2  ! 1- coarse dust; 2- supercoarse dust
42  INTEGER, PARAMETER :: naero = naero_soluble + naero_insoluble
43
44  INTEGER, PARAMETER :: nbre_RH = 12
45  REAL, PARAMETER :: RH_tab(nbre_RH) = (/0., 10., 20., 30., 40., 50., 60., 70., 80., 85., 90., 95./)
46  REAL, PARAMETER :: RH_MAX = 95.
47  REAL :: delta(klon, klev), rh(klon, klev)
48  REAL :: tau_ae2b_int   ! Intermediate computation of epaisseur optique aerosol
49  REAL :: piz_ae2b_int   ! Intermediate computation of Single scattering albedo
50  REAL :: cg_ae2b_int    ! Intermediate computation of Assymetry parameter
51  REAL :: fact_RH(nbre_RH), tmp_var
52
53  REAL, DIMENSION(klon, klev, naero_tot, nbands_sw_rrtm) :: tau_ae
54  REAL, DIMENSION(klon, klev, naero_tot, nbands_sw_rrtm) :: piz_ae
55  REAL, DIMENSION(klon, klev, naero_tot, nbands_sw_rrtm) :: cg_ae
56
57  ! Proprietes optiques
58
59  REAL :: alpha_aers_6bands(nbre_RH, nbands_sw_rrtm, naero_soluble)   !--unit m2/g aer
60  REAL :: alpha_aeri_6bands(nbands_sw_rrtm, naero_insoluble)         !--unit m2/g aer
61  REAL :: cg_aers_6bands(nbre_RH, nbands_sw_rrtm, naero_soluble)      !--unitless
62  REAL :: cg_aeri_6bands(nbands_sw_rrtm, naero_insoluble)            !--unitless
63  REAL :: piz_aers_6bands(nbre_RH, nbands_sw_rrtm, naero_soluble)     !--unitless
64  REAL :: piz_aeri_6bands(nbands_sw_rrtm, naero_insoluble)           !--unitless
65
66  REAL, PARAMETER :: tau_min = 1.e-8
67  CHARACTER*20 modname
68
69  !***************************************************************************
70  !--the order of the soluble   species has to follow the spsol   index below
71  !--the order of the insoluble species has to follow the spinsol index below
72
73  DATA alpha_aers_6bands/  &
74          ! accumulation (sulfate+2% bc) mode soluble
75          5.212, 5.212, 5.212, 5.212, 6.973, 7.581, 8.349, 9.400, 11.078, 12.463, 14.857, 20.837, &
76          4.906, 4.906, 4.906, 4.906, 6.568, 7.195, 7.989, 9.088, 10.869, 12.354, 14.951, 21.545, &
77          3.940, 3.940, 3.940, 3.940, 5.291, 5.861, 6.591, 7.620, 9.332, 10.791, 13.410, 20.370, &
78          2.292, 2.292, 2.292, 2.292, 3.105, 3.493, 3.996, 4.724, 5.978, 7.084, 9.146, 15.067, &
79          0.762, 0.762, 0.762, 0.762, 1.050, 1.201, 1.401, 1.699, 2.232, 2.721, 3.675, 6.678, &
80          0.090, 0.090, 0.090, 0.090, 0.122, 0.141, 0.166, 0.204, 0.275, 0.344, 0.484, 0.973, &
81          ! coarse seasalt
82          0.547, 0.657, 0.705, 0.754, 0.817, 0.896, 1.008, 1.169, 1.456, 1.724, 2.199, 3.358, &
83          0.566, 0.679, 0.727, 0.776, 0.840, 0.920, 1.032, 1.196, 1.492, 1.760, 2.238, 3.416, &
84          0.596, 0.714, 0.764, 0.816, 0.882, 0.965, 1.081, 1.250, 1.552, 1.828, 2.310, 3.509, &
85          0.644, 0.771, 0.825, 0.880, 0.951, 1.040, 1.164, 1.345, 1.666, 1.957, 2.462, 3.700, &
86          0.640, 0.772, 0.829, 0.887, 0.965, 1.061, 1.198, 1.398, 1.758, 2.085, 2.658, 4.031, &
87          0.452, 0.562, 0.609, 0.659, 0.728, 0.813, 0.938, 1.125, 1.471, 1.797, 2.384, 3.855  /
88
89  DATA alpha_aeri_6bands/  &
90          ! coarse dust insoluble
91          0.594, 0.610, 0.619, 0.762, 0.791, 0.495, &
92          ! supercoarse dust insoluble
93          0.151, 0.152, 0.155, 0.159, 0.168, 0.167  /
94
95  DATA cg_aers_6bands/ &
96          ! accumulation (sulfate+2% bc) mode soluble
97          0.692, 0.692, 0.692, 0.692, 0.735, 0.739, 0.744, 0.749, 0.755, 0.759, 0.765, 0.772, &
98          0.690, 0.690, 0.690, 0.690, 0.736, 0.740, 0.746, 0.752, 0.760, 0.765, 0.771, 0.779, &
99          0.678, 0.678, 0.678, 0.678, 0.727, 0.733, 0.740, 0.748, 0.759, 0.766, 0.775, 0.787, &
100          0.641, 0.641, 0.641, 0.641, 0.692, 0.700, 0.710, 0.721, 0.736, 0.746, 0.760, 0.781, &
101          0.553, 0.553, 0.553, 0.553, 0.603, 0.615, 0.627, 0.643, 0.664, 0.678, 0.699, 0.735, &
102          0.343, 0.343, 0.343, 0.343, 0.386, 0.399, 0.414, 0.433, 0.460, 0.480, 0.510, 0.569, &
103          ! seasalt coarse Soluble
104          0.754, 0.770, 0.776, 0.781, 0.784, 0.791, 0.797, 0.805, 0.815, 0.822, 0.828, 0.840, &
105          0.736, 0.753, 0.759, 0.765, 0.771, 0.778, 0.785, 0.793, 0.804, 0.811, 0.820, 0.831, &
106          0.716, 0.735, 0.742, 0.748, 0.754, 0.762, 0.769, 0.778, 0.789, 0.796, 0.807, 0.819, &
107          0.704, 0.725, 0.733, 0.739, 0.745, 0.752, 0.759, 0.768, 0.778, 0.784, 0.792, 0.803, &
108          0.716, 0.737, 0.744, 0.751, 0.756, 0.763, 0.770, 0.777, 0.786, 0.790, 0.795, 0.800, &
109          0.688, 0.730, 0.741, 0.751, 0.761, 0.771, 0.782, 0.795, 0.810, 0.820, 0.833, 0.849  /
110
111  DATA cg_aeri_6bands/ &
112          ! coarse dust insoluble
113          0.801, 0.779, 0.709, 0.698, 0.710, 0.687, &
114          ! super coarse dust insoluble
115          0.862, 0.871, 0.852, 0.799, 0.758, 0.651  /
116
117  DATA piz_aers_6bands/&
118          ! accumulation (sulfate+2% bc) mode soluble
119          0.941, 0.941, 0.941, 0.941, 0.958, 0.961, 0.965, 0.969, 0.974, 0.977, 0.981, 0.987, &
120          0.941, 0.941, 0.941, 0.941, 0.959, 0.963, 0.967, 0.971, 0.976, 0.979, 0.983, 0.988, &
121          0.953, 0.953, 0.953, 0.953, 0.967, 0.971, 0.974, 0.978, 0.982, 0.984, 0.988, 0.992, &
122          0.955, 0.955, 0.955, 0.955, 0.969, 0.972, 0.976, 0.980, 0.984, 0.986, 0.989, 0.994, &
123          0.936, 0.936, 0.936, 0.936, 0.955, 0.961, 0.966, 0.972, 0.978, 0.982, 0.987, 0.993, &
124          0.792, 0.792, 0.792, 0.792, 0.848, 0.867, 0.887, 0.907, 0.931, 0.944, 0.960, 0.980, &
125          ! seasalt coarse soluble
126          1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.001, 1.000, &
127          1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
128          1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
129          1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, 1.000, &
130          0.994, 0.994, 0.995, 0.995, 0.995, 0.995, 0.996, 0.996, 0.996, 0.996, 0.996, 0.996, &
131          0.976, 0.867, 0.837, 0.814, 0.796, 0.774, 0.754, 0.735, 0.713, 0.702, 0.690, 0.675  /
132
133  DATA piz_aeri_6bands/ &
134          ! coarse dust insoluble
135          0.866, 0.875, 0.915, 0.977, 0.993, 0.971, &
136          ! super coarse dust insoluble
137          0.789, 0.749, 0.791, 0.918, 0.970, 0.909  /
138
139  spsol = 0
140  spinsol = 0
141
142  modname = 'splaeropt_6bands_rrt'
143
144  IF (NSW/=nbands_sw_rrtm) THEN
145    CALL abort_physic(modname, 'Erreur NSW doit etre egal a 6 pour cette routine', 1)
146  ENDIF
147
148  DO irh = 1, nbre_RH - 1
149    fact_RH(irh) = 1. / (RH_tab(irh + 1) - RH_tab(irh))
150  ENDDO
151
152  DO k = 1, klev
153    DO i = 1, klon
154      rh(i, k) = MIN(RHcl(i, k) * 100., RH_MAX)
155      RH_num(i, k) = INT(rh(i, k) / 10. + 1.)
156      IF (rh(i, k)>85.) RH_num(i, k) = 10
157      IF (rh(i, k)>90.) RH_num(i, k) = 11
158      delta(i, k) = (rh(i, k) - RH_tab(RH_num(i, k))) * fact_RH(RH_num(i, k))
159    ENDDO
160  ENDDO
161
162  tau_ae(:, :, :, :) = 0.
163  piz_ae(:, :, :, :) = 0.
164  cg_ae(:, :, :, :) = 0.
165
166  itr = 0
167  DO iq = 1, nqtot
168    IF(.NOT.tracers(iq)%isInPhysics) CYCLE
169    itr = itr + 1
170    SELECT CASE(tracers(iq)%name)
171    CASE('PREC'); CYCLE                                  !--precursor
172    CASE('FINE'); soluble = .TRUE.;  spsol = 1; aerindex = 1   !--fine mode accumulation mode
173    CASE('COSS'); soluble = .TRUE.;  spsol = 2; aerindex = 2   !--coarse mode sea salt
174    CASE('CODU'); soluble = .FALSE.; spinsol = 1; aerindex = 3   !--coarse mode dust
175    CASE('SCDU'); soluble = .FALSE.; spinsol = 2; aerindex = 4   !--super coarse mode dust
176    CASE DEFAULT; CALL abort_physic(modname, 'I cannot do aerosol optics for ' // tracers(iq)%name, 1)
177    END SELECT
178
179    IF (soluble) THEN ! For aerosol soluble components
180
181      DO k = 1, klev
182        DO i = 1, klon
183
184          tmp_var = tr_seri(i, k, itr) * zdm(i, k) !-- g/m2
185
186          DO inu = 1, NSW
187
188            tau_ae2b_int = alpha_aers_6bands(RH_num(i, k), inu, spsol) + &
189                    delta(i, k) * (alpha_aers_6bands(RH_num(i, k) + 1, inu, spsol) - &
190                            alpha_aers_6bands(RH_num(i, k), inu, spsol))
191
192            piz_ae2b_int = piz_aers_6bands(RH_num(i, k), inu, spsol) + &
193                    delta(i, k) * (piz_aers_6bands(RH_num(i, k) + 1, inu, spsol) - &
194                            piz_aers_6bands(RH_num(i, k), inu, spsol))
195
196            cg_ae2b_int = cg_aers_6bands(RH_num(i, k), inu, spsol) + &
197                    delta(i, k) * (cg_aers_6bands(RH_num(i, k) + 1, inu, spsol) - &
198                            cg_aers_6bands(RH_num(i, k), inu, spsol))
199
200            tau_ae(i, k, aerindex, inu) = tmp_var * tau_ae2b_int
201            piz_ae(i, k, aerindex, inu) = piz_ae2b_int
202            cg_ae(i, k, aerindex, inu) = cg_ae2b_int
203
204          ENDDO
205        ENDDO
206      ENDDO
207
208    ELSE    ! For all aerosol insoluble components
209
210      DO k = 1, klev
211        DO i = 1, klon
212
213          tmp_var = tr_seri(i, k, itr) * zdm(i, k) !-- g/m2
214
215          DO inu = 1, NSW
216            tau_ae2b_int = alpha_aeri_6bands(inu, spinsol)
217            piz_ae2b_int = piz_aeri_6bands(inu, spinsol)
218            cg_ae2b_int = cg_aeri_6bands(inu, spinsol)
219
220            tau_ae(i, k, aerindex, inu) = tmp_var * tau_ae2b_int
221            piz_ae(i, k, aerindex, inu) = piz_ae2b_int
222            cg_ae(i, k, aerindex, inu) = cg_ae2b_int
223          ENDDO
224        ENDDO
225      ENDDO
226
227    ENDIF ! soluble / insoluble
228
229  ENDDO  ! nbtr
230
231  !--all (natural + anthropogenic) aerosol
232  tau_allaer(:, :, 2, :) = SUM(tau_ae(:, :, 1:naero, :), dim = 3)
233  tau_allaer(:, :, 2, :) = MAX(tau_allaer(:, :, 2, :), tau_min)
234
235  piz_allaer(:, :, 2, :) = SUM(tau_ae(:, :, 1:naero, :) * piz_ae(:, :, 1:naero, :), dim = 3) / tau_allaer(:, :, 2, :)
236  piz_allaer(:, :, 2, :) = MIN(MAX(piz_allaer(:, :, 2, :), 0.01), 1.0)
237  WHERE (tau_allaer(:, :, 2, :)<=tau_min) piz_allaer(:, :, 2, :) = 1.0
238
239  cg_allaer(:, :, 2, :) = SUM(tau_ae(:, :, 1:naero, :) * piz_ae(:, :, 1:naero, :) * cg_ae(:, :, 1:naero, :), dim = 3) / &
240          (tau_allaer(:, :, 2, :) * piz_allaer(:, :, 2, :))
241  cg_allaer(:, :, 2, :) = MIN(MAX(cg_allaer(:, :, 2, :), 0.0), 1.0)
242
243  !--no aerosol
244  tau_allaer(:, :, 1, :) = tau_min
245  piz_allaer(:, :, 1, :) = 1.0
246  cg_allaer(:, :, 1, :) = 0.0
247
248  !--waveband 2 and all aerosol (third index = 2)
249  inu = 2
250  abs550aer(:) = SUM((1 - piz_allaer(:, :, 2, inu)) * tau_allaer(:, :, 2, inu), dim = 2)
251
252END SUBROUTINE SPLAEROPT_6BANDS_RRTM
Note: See TracBrowser for help on using the repository browser.