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

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

Introducing stratomask diagnostic on where the stratosphere is
flag_aerosol_strat = 2 for CMIP6 strat aerosol forcing

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