source: LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato_rrtm.F90 @ 2231

Last change on this file since 2231 was 2231, checked in by oboucher, 9 years ago

Putting minimum values for SW and LW aerosol optical depth values
Putting default value for aerosol single scattering albedo if no aerosol
Doing this consistently in all cases covered by flag_aerosol and flag_aerosol_strat
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
  • Property svn:keywords set to Author Date Id Revi
File size: 9.1 KB
Line 
1!
2! $Id: readaerosolstrato_rrtm.F90 2231 2015-03-12 16:46:56Z oboucher $
3!
4subroutine readaerosolstrato_rrtm(debut)
5
6    use netcdf95, only: nf95_close, nf95_gw_var, nf95_inq_dimid, &
7                        nf95_inq_varid, nf95_open
8    use netcdf, only: nf90_get_var, nf90_noerr, nf90_nowrite
9
10    USE phys_cal_mod, ONLY : mth_cur
11    USE mod_grid_phy_lmdz
12    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
13    USE mod_phys_lmdz_para
14    USE phys_state_var_mod
15    USE phys_local_var_mod
16    USE aero_mod
17    USE dimphy
18    USE YOERAD   , ONLY : NLW
19
20    implicit none
21
22    include "YOMCST.h"
23    include "dimensions.h"
24
25! Variable input
26    logical debut
27
28! Variables locales
29    integer n_lat   ! number of latitudes in the input data
30    integer n_lon   ! number of longitudes in the input data
31    integer n_lev   ! number of levels in the input data
32    integer n_month ! number of months in the input data
33    real, pointer:: latitude(:)
34    real, pointer:: longitude(:)
35    real, pointer:: time(:)
36    real, pointer:: lev(:)
37    integer k, band, wave, i
38    integer, save :: mth_pre
39
40    real, allocatable, dimension(:,:), save :: tau_aer_strat
41!$OMP THREADPRIVATE(tau_aer_strat)
42
43! Champs reconstitues
44    real, allocatable:: tauaerstrat(:, :, :, :)
45    real, allocatable:: tauaerstrat_mois(:, :, :)
46    real, allocatable:: tauaerstrat_mois_glo(:, :)
47
48    real, allocatable:: sum_tau_aer_strat(:)
49
50! For NetCDF:
51    integer ncid_in  ! IDs for input files
52    integer varid, ncerr
53
54! Stratospheric aerosols optical properties
55! alpha_sw_strat over the 6 bands is normalised by the 550 nm extinction coefficient
56    real, dimension(nbands_sw_rrtm) :: alpha_sw_strat, piz_sw_strat, cg_sw_strat
57    data alpha_sw_strat/0.8545564, 0.8451642, 0.9821724, 0.8145110, 0.3073565, 7.7966176E-02/
58    data cg_sw_strat   /0.6997170, 0.6810035, 0.7403592, 0.7562674, 0.6676504, 0.3478689/
59    data piz_sw_strat  /0.9999998, 0.9999998, 1.000000000, 0.9999958, 0.9977155, 0.4510679/
60!
61!--diagnostics AOD in the SW
62! alpha_sw_strat_wave is *not* normalised by the 550 nm extinction coefficient
63    real, dimension(nwave) :: alpha_sw_strat_wave
64    data alpha_sw_strat_wave/3.708007,4.125824,4.136584,3.887478,3.507738/
65!
66!--diagnostics AOD in the LW at 10 um
67    real :: alpha_lw_strat_wave
68    data alpha_lw_strat_wave/0.2746812/
69!
70    real, dimension(nbands_lw_rrtm) :: alpha_lw_abs_rrtm
71    data alpha_lw_abs_rrtm/   8.8340312E-02, 6.9856711E-02, 6.2652975E-02, 5.7188231E-02, &
72                              6.3157059E-02, 5.5072524E-02, 5.0571125E-02, 0.1349073, &   
73                              0.1381676, 9.6506312E-02, 5.1312990E-02, 2.4256418E-02, &
74                              2.7191756E-02, 3.3862915E-02, 1.6132960E-02, 1.4275438E-02/ ! calculated with Mie_SW_LW_RRTM_V2.4 (bimodal, corrected)
75                                                                                          ! for r_0=/0.13E-6, 0.41E-6/ m, sigma_g=/1.26, 1.30/
76                                                                                          ! order: increasing wavelength!
77!--------------------------------------------------------
78
79    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev))
80    IF (.not.ALLOCATED(sum_tau_aer_strat)) ALLOCATE(sum_tau_aer_strat(klon))
81
82    IF (debut.OR.mth_cur.NE.mth_pre) THEN
83
84    IF (is_mpi_root) THEN
85
86    IF (nbands_sw_rrtm.NE.6) THEN
87        print *,'nbands_sw_rrtm doit etre egal a 6 dans readaerosolstrat_rrtm'
88        STOP
89    ENDIF
90
91    CALL nf95_open("taustrat.nc", nf90_nowrite, ncid_in)
92
93    CALL nf95_inq_varid(ncid_in, "LEV", varid)
94    CALL nf95_gw_var(ncid_in, varid, lev)
95    n_lev = size(lev)
96    IF (n_lev.NE.klev) THEN
97       print *,'Le nombre de niveaux n est pas egal a klev'
98       STOP
99    ENDIF
100
101    CALL nf95_inq_varid(ncid_in, "LAT", varid)
102    CALL nf95_gw_var(ncid_in, varid, latitude)
103    n_lat = size(latitude)
104    print *, 'LAT aerosol strato=', n_lat, latitude
105    IF (n_lat.NE.jjm+1) THEN
106       print *,'Le nombre de lat n est pas egal a jjm+1'
107       STOP
108    ENDIF
109
110    CALL nf95_inq_varid(ncid_in, "LON", varid)
111    CALL nf95_gw_var(ncid_in, varid, longitude)
112    n_lon = size(longitude)
113    print *, 'LON aerosol strato=', n_lon, longitude
114    IF (n_lon.NE.iim) THEN
115       print *,'Le nombre de lon n est pas egal a iim'
116       STOP
117    ENDIF
118
119    CALL nf95_inq_varid(ncid_in, "TIME", varid)
120    CALL nf95_gw_var(ncid_in, varid, time)
121    n_month = size(time)
122    print *, 'TIME aerosol strato=', n_month, time
123    IF (n_month.NE.12) THEN
124       print *,'Le nombre de month n est pas egal a 12'
125       STOP
126    ENDIF
127
128    ALLOCATE(tauaerstrat(n_lon, n_lat, n_lev, n_month))
129    ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev))
130    ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev))
131
132!--reading stratospheric AOD at 550 nm
133    CALL nf95_inq_varid(ncid_in, "TAUSTRAT", varid)
134    ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
135    print *,'code erreur readaerosolstrato=', ncerr, varid
136
137    CALL nf95_close(ncid_in)
138
139!---select the correct month
140    IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
141      print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
142    ENDIF
143    tauaerstrat_mois(:,:,:) = tauaerstrat(:,:,:,mth_cur)
144
145!---reduce to a klon_glo grid
146    CALL grid2dTo1d_glo(tauaerstrat_mois,tauaerstrat_mois_glo)
147
148    ENDIF !--is_mpi_root
149
150!--keep memory of previous month
151    mth_pre=mth_cur
152
153!--scatter on all proc
154    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
155
156    IF (is_mpi_root) THEN
157!
158    DEALLOCATE(tauaerstrat)
159    DEALLOCATE(tauaerstrat_mois)
160    DEALLOCATE(tauaerstrat_mois_glo)
161!
162    ENDIF !--is_mpi_root
163
164    ENDIF !--debut ou nouveau mois
165
166!--total vertical aod at the 5 SW wavelengths
167    DO wave=1, nwave
168    DO k=1, klev
169    tausum_aero(:,wave,id_STRAT_phy)=tausum_aero(:,wave,id_STRAT_phy)+ &
170       tau_aer_strat(:,k)*alpha_sw_strat_wave(wave)/alpha_sw_strat_wave(2)
171    ENDDO
172    ENDDO
173
174!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
175    DO band=1, nbands_sw_rrtm
176!--anthropogenic aerosols bands 1 to nbands_sw_rrtm
177    cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
178                                  cg_sw_strat(band)*piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /           &
179                             MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                                &
180                                  piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
181    piz_aero_sw_rrtm(:,:,2,band)  = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                            &
182                              piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /                                 &
183                              MAX( tau_aero_sw_rrtm(:,:,2,band) + alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
184    tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + alpha_sw_strat(band)*tau_aer_strat(:,:)
185!--natural aerosols bands 1 to nbands_sw_rrtm
186    cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
187                             cg_sw_strat(band)*piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /                &
188                             MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                                &
189                                  piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:), 1.e-15 )
190    piz_aero_sw_rrtm(:,:,1,band)  = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                            &
191                              piz_sw_strat(band)*alpha_sw_strat(band)*tau_aer_strat(:,:) ) /                                 &
192                              MAX( tau_aero_sw_rrtm(:,:,1,band) + alpha_sw_strat(band)*tau_aer_strat(:,:),1.e-15 )
193    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + alpha_sw_strat(band)*tau_aer_strat(:,:)
194!--no stratospheric aerosol in index 1 for these tests
195!    cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
196!    piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
197!    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
198    ENDDO
199
200!--stratospheric AOD in LW
201    IF (nbands_lw_rrtm .NE. NLW) then
202      print*, 'different values for NLW (=',NLW,') and nbands_lw_rrtm (=', nbands_lw_rrtm, ')'
203      STOP
204    ENDIF
205
206    DO band=1, nbands_lw_rrtm
207    tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:)
208    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + alpha_lw_abs_rrtm(band)*tau_aer_strat(:,:)
209!--no stratospheric aerosols in index 1 for these tests
210!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band)
211    ENDDO
212
213!--default SSA value if there is no aerosol
214!--to avoid 0 values that seems to cause some problem to RRTM
215    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
216      piz_aero_sw_rrtm = 1.0
217    ENDWHERE
218
219!--in principle this should not be necessary
220!--as these variables have min values already but just in case
221!--put 1e-15 min value to both SW and LW AOD
222    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
223    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
224
225end subroutine readaerosolstrato_rrtm
Note: See TracBrowser for help on using the repository browser.