source: LMDZ6/branches/IPSLCM6.0.15/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 @ 3607

Last change on this file since 3607 was 3427, checked in by oboucher, 6 years ago

Correcting a bug for strat aerosol forcing diagnostics

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