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

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

Adding threadprivate for mth_pre variable
This should fix some reproducibility effects for runs longer than 1 months....

File size: 13.8 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    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
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:: time(:)
40    REAL, POINTER:: lev(:)
41    REAL, POINTER:: wav(:)
42    INTEGER i,k,wave,band
43    INTEGER, SAVE :: mth_pre=1
44!$OMP THREADPRIVATE(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
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))
77
78    IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,NLW))
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
127        IF (n_wav.NE.NSW) THEN
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
175        DO band=1, NSW
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
213        IF (n_wav.NE.NLW) THEN
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
238        DO band=1, NLW
239          CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
240        ENDDO
241
242      ELSE !--proc other than mpi_root and omp_root
243           !--dummy allocation needed for debug mode
244
245        ALLOCATE(tauaerstrat_mois_glo(1,1,1))
246        ALLOCATE(pizaerstrat_mois_glo(1,1,1))
247        ALLOCATE(cgaerstrat_mois_glo(1,1,1))
248        ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
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(taulwaerstrat,taulwaerstrat_mois)
268!
269      ENDIF !--is_mpi_root and is_omp_root
270
271      DEALLOCATE(tauaerstrat_mois_glo,pizaerstrat_mois_glo,cgaerstrat_mois_glo)
272      DEALLOCATE(taulwaerstrat_mois_glo)
273
274!$OMP BARRIER
275
276    ENDIF !--debut ou nouveau mois
277
278!--total vertical aod at the 5 SW wavelengths
279!--for now use band 3 AOD into all 5 wavelengths
280!--it is only a reasonable approximation for 550 nm (wave=2)
281    band=3
282    DO i=1, klon
283    DO k=1, klev
284      IF (stratomask(i,k).GT.0.999999) THEN
285        DO wave=1, nwave_sw
286          tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
287        ENDDO
288      ENDIF
289    ENDDO
290    ENDDO
291
292!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
293    DO band=1, NSW
294      WHERE (stratomask.GT.0.999999)
295!--anthropogenic aerosols bands 1 to NSW
296        cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
297                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
298                                    MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
299                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
300        piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
301                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
302                                    MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
303        tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
304!--natural aerosols bands 1 to NSW
305        cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
306                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
307                                    MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
308                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
309        piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
310                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
311                                    MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
312        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
313!--no stratospheric aerosol in index 1 for these tests
314!        cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
315!        piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
316!        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
317    ENDWHERE
318    ENDDO
319
320!--total vertical aod at 10 um
321!--this is approximated from band 7 of RRTM
322    band=7
323    DO i=1, klon
324    DO k=1, klev
325      IF (stratomask(i,k).GT.0.999999) THEN
326        DO wave=1, nwave_lw
327          tausum_aero(i,nwave_sw+wave,id_STRAT_phy)=tausum_aero(i,nwave_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
328        ENDDO
329      ENDIF
330    ENDDO
331    ENDDO
332
333    DO band=1, NLW
334      WHERE (stratomask.GT.0.999999)
335        tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
336        tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
337!--no stratospheric aerosols in index 1 for these tests
338!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band)
339      ENDWHERE
340    ENDDO
341
342!--default SSA value if there is no aerosol
343!--to avoid 0 values that seems to cause some problem to RRTM
344    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
345      piz_aero_sw_rrtm = 1.0
346    ENDWHERE
347
348!--in principle this should not be necessary
349!--as these variables have min values already but just in case
350!--put 1e-15 min value to both SW and LW AOD
351    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
352    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
353
354END SUBROUTINE readaerosolstrato2_rrtm
Note: See TracBrowser for help on using the repository browser.