source: LMDZ5/trunk/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 @ 2639

Last change on this file since 2639 was 2550, checked in by oboucher, 8 years ago

Correcting a bug in the calculation of od550_strat
Computing od_10um_strat in readaerosolstrato*.F90

File size: 12.9 KB
RevLine 
[2536]1!
2! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
3!
4subroutine readaerosolstrato2_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, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo
12    USE mod_phys_lmdz_mpi_data, ONLY :  is_mpi_root
13    USE mod_phys_lmdz_omp_data, ONLY :  is_omp_root
14    USE mod_phys_lmdz_para
15    USE phys_state_var_mod
16    USE phys_local_var_mod
17    USE aero_mod
18    USE dimphy
19    USE YOERAD   , ONLY : NLW
20
21    implicit none
22
23    include "YOMCST.h"
[2539]24    INCLUDE "clesphys.h"
[2536]25
26    CHARACTER (len = 80) :: abort_message
27    CHARACTER (LEN=20) :: modname = 'readaerosolstrato2'
28
29! Variable input
30    logical, intent(in) ::  debut
31
32! Variables locales
33    integer n_lat   ! number of latitudes in the input data
34    integer n_lon   ! number of longitudes
35    integer n_lev   ! number of levels in the input data
36    integer n_month ! number of months in the input data
37    integer n_wav   ! number of wavelengths in the input data
38    real, pointer:: latitude(:)
39    real, pointer:: longitude(:)
40    real, pointer:: time(:)
41    real, pointer:: lev(:)
42    real, pointer:: wav(:)
43    integer i,k,wave,band
44    integer, save :: mth_pre
45
46    real, allocatable, dimension(:,:,:), save :: tau_aer_strat
47    real, allocatable, dimension(:,:,:), save :: piz_aer_strat
48    real, allocatable, dimension(:,:,:), save :: cg_aer_strat
49    real, allocatable, dimension(:,:,:), save :: taulw_aer_strat
50!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat)
51
52! Champs reconstitues
53    real, allocatable:: tauaerstrat(:, :, :, :)
54    real, allocatable:: pizaerstrat(:, :, :, :)
55    real, allocatable:: cgaerstrat(:, :, :, :)
56    real, allocatable:: taulwaerstrat(:, :, :, :)
57
58    real, allocatable:: tauaerstrat_mois(:, :, :, :)
59    real, allocatable:: pizaerstrat_mois(:, :, :, :)
60    real, allocatable:: cgaerstrat_mois(:, :, :, :)
61    real, allocatable:: taulwaerstrat_mois(:, :, :, :)
62
63    real, allocatable:: tauaerstrat_mois_glo(:, :, :)
64    real, allocatable:: pizaerstrat_mois_glo(:, :, :)
65    real, allocatable:: cgaerstrat_mois_glo(:, :, :)
66    real, allocatable:: taulwaerstrat_mois_glo(:, :, :)
67
68! For NetCDF:
69    integer ncid_in  ! IDs for input files
70    integer varid, ncerr
71
72!--------------------------------------------------------
73
[2539]74    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev,NSW))
75    IF (.not.ALLOCATED(piz_aer_strat)) ALLOCATE(piz_aer_strat(klon,klev,NSW))
76    IF (.not.ALLOCATED(cg_aer_strat))  ALLOCATE(cg_aer_strat(klon,klev,NSW))
[2536]77
[2539]78    IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,NLW))
[2536]79
80!--we only read monthly strat aerosol data
81    IF (debut.OR.mth_cur.NE.mth_pre) THEN
82
83!--only root reads the data
84    IF (is_mpi_root.AND.is_omp_root) THEN
85
86!--check mth_cur
87    IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
88      print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
89    ENDIF
90
91!--initialize n_lon as input data is 2D (lat-alt) only
92    n_lon = nbp_lon
93
94!--Starts with SW optical properties
95
96    CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
97
98    CALL nf95_inq_varid(ncid_in, "LEV", varid)
99    CALL nf95_gw_var(ncid_in, varid, lev)
100    n_lev = size(lev)
101    IF (n_lev.NE.klev) THEN
102       abort_message='Le nombre de niveaux n est pas egal a klev'
103       CALL abort_physic(modname,abort_message,1)
104    ENDIF
105
106    CALL nf95_inq_varid(ncid_in, "LAT", varid)
107    CALL nf95_gw_var(ncid_in, varid, latitude)
108    n_lat = size(latitude)
109    IF (n_lat.NE.nbp_lat) THEN
110       print *, 'latitude=', n_lat, nbp_lat
111       abort_message='Le nombre de lat n est pas egal a nbp_lat'
112       CALL abort_physic(modname,abort_message,1)
113    ENDIF
114
115    CALL nf95_inq_varid(ncid_in, "TIME", varid)
116    CALL nf95_gw_var(ncid_in, varid, time)
117    n_month = size(time)
118    IF (n_month.NE.12) THEN
119       abort_message='Le nombre de month n est pas egal a 12'
120       CALL abort_physic(modname,abort_message,1)
121    ENDIF
122
123    CALL nf95_inq_varid(ncid_in, "WAV", varid)
124    CALL nf95_gw_var(ncid_in, varid, wav)
125    n_wav = size(wav)
126    print *, 'WAV aerosol strato=', n_wav, wav
[2539]127    IF (n_wav.NE.NSW) THEN
[2536]128       abort_message='Le nombre de wav n est pas egal a NSW'
129       CALL abort_physic(modname,abort_message,1)
130    ENDIF
131
132    ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
133    ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
134    ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
135
136    ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
137    ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
138    ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
139
140    ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
141    ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
142    ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
143
144!--reading stratospheric aerosol tau per layer
145    CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
146    ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
147    print *,'code erreur readaerosolstrato=', ncerr, varid
148
149!--reading stratospheric aerosol omega per layer
150    CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
151    ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
152    print *,'code erreur readaerosolstrato=', ncerr, varid
153
154!--reading stratospheric aerosol g per layer
155    CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
156    ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
157    print *,'code erreur readaerosolstrato sw=', ncerr, varid
158
159    CALL nf95_close(ncid_in)
160
161!--select the correct month
162!--and copy into 1st longitude
163    tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
164    pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
165    cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
166
167!--copy longitudes
168    DO i=2, n_lon
169     tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
170     pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
171      cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
172    ENDDO
173
174!---reduce to a klon_glo grid
[2539]175    DO band=1, NSW
[2536]176      CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
177      CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
178      CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
179    ENDDO
180
181!--Now LW optical properties
182!
183    CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
184
185    CALL nf95_inq_varid(ncid_in, "LEV", varid)
186    CALL nf95_gw_var(ncid_in, varid, lev)
187    n_lev = size(lev)
188    IF (n_lev.NE.klev) THEN
189       abort_message='Le nombre de niveaux n est pas egal a klev'
190       CALL abort_physic(modname,abort_message,1)
191    ENDIF
192
193    CALL nf95_inq_varid(ncid_in, "LAT", varid)
194    CALL nf95_gw_var(ncid_in, varid, latitude)
195    n_lat = size(latitude)
196    IF (n_lat.NE.nbp_lat) THEN
197       abort_message='Le nombre de lat n est pas egal a nbp_lat'
198       CALL abort_physic(modname,abort_message,1)
199    ENDIF
200
201    CALL nf95_inq_varid(ncid_in, "TIME", varid)
202    CALL nf95_gw_var(ncid_in, varid, time)
203    n_month = size(time)
204    IF (n_month.NE.12) THEN
205       abort_message='Le nombre de month n est pas egal a 12'
206       CALL abort_physic(modname,abort_message,1)
207    ENDIF
208
209    CALL nf95_inq_varid(ncid_in, "WAV", varid)
210    CALL nf95_gw_var(ncid_in, varid, wav)
211    n_wav = size(wav)
212    print *, 'WAV aerosol strato=', n_wav, wav
[2539]213    IF (n_wav.NE.NLW) THEN
[2536]214       abort_message='Le nombre de wav n est pas egal a NLW'
215       CALL abort_physic(modname,abort_message,1)
216    ENDIF
217
218    ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
219    ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
220    ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
221
222!--reading stratospheric aerosol lw tau per layer
223    CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
224    ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
225    print *,'code erreur readaerosolstrato lw=', ncerr, varid
226
227    CALL nf95_close(ncid_in)
228
229!--select the correct month
230!--and copy into 1st longitude
231    taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
232!--copy longitudes
233    DO i=2, n_lon
234      taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
235    ENDDO
236
237!---reduce to a klon_glo grid
[2539]238    DO band=1, NLW
[2536]239      CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
240    ENDDO
241
242    ENDIF !--is_mpi_root and is_omp_root
243
244!$OMP BARRIER
245
246!--keep memory of previous month
247    mth_pre=mth_cur
248
249!--scatter on all proc
250    CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
251    CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
252    CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
253    CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
254
255    IF (is_mpi_root.AND.is_omp_root) THEN
256!
257    DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat)
258    DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
259    DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
260
261    DEALLOCATE(taulwaerstrat,taulwaerstrat_mois,taulwaerstrat_mois_glo)
262!
263    ENDIF !--is_mpi_root and is_omp_root
264
265!$OMP BARRIER
266
267    ENDIF !--debut ou nouveau mois
268
269!--total vertical aod at the 5 SW wavelengths
270!--for now use band 3 AOD into all 5 wavelengths
[2550]271!--it is only a reasonable approximation for 550 nm (wave=2)
[2536]272    band=3
273    DO i=1, klon
274    DO k=1, klev
275      IF (stratomask(i,k).GT.0.999999) THEN
[2550]276        DO wave=1, nwave_sw
277          tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
[2536]278        ENDDO
279      ENDIF
280    ENDDO
281    ENDDO
282
283!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
[2539]284    DO band=1, NSW
[2536]285    WHERE (stratomask.GT.0.999999)
[2539]286!--anthropogenic aerosols bands 1 to NSW
[2536]287    cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
288                                     cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
289                                MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
290                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
291    piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
292                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
293                                MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
294    tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
[2539]295!--natural aerosols bands 1 to NSW
[2536]296    cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
297                                     cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
298                                MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
299                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
300    piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
301                                     piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
302                                MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
303    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
304!--no stratospheric aerosol in index 1 for these tests
305!    cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
306!    piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
307!    tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
308    ENDWHERE
309    ENDDO
310
[2550]311!--total vertical aod at 10 um
312!--this is approximated from band 7 of RRTM
313    band=7
314    DO i=1, klon
315    DO k=1, klev
316      IF (stratomask(i,k).GT.0.999999) THEN
317        DO wave=1, nwave_lw
318          tausum_aero(i,nwave_sw+wave,id_STRAT_phy)=tausum_aero(i,nwave_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
319        ENDDO
320      ENDIF
321    ENDDO
322    ENDDO
323
[2539]324    DO band=1, NLW
[2536]325    WHERE (stratomask.GT.0.999999)
326    tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
327    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
328!--no stratospheric aerosols in index 1 for these tests
329!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band)
330    ENDWHERE
331    ENDDO
332
333!--default SSA value if there is no aerosol
334!--to avoid 0 values that seems to cause some problem to RRTM
335    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
336      piz_aero_sw_rrtm = 1.0
337    ENDWHERE
338
339!--in principle this should not be necessary
340!--as these variables have min values already but just in case
341!--put 1e-15 min value to both SW and LW AOD
342    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
343    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
344
345end subroutine readaerosolstrato2_rrtm
Note: See TracBrowser for help on using the repository browser.