GCC Code Coverage Report


Directory: ./
File: rad/readaerosolstrato2_rrtm.f90
Date: 2022-01-11 19:19:34
Exec Total Coverage
Lines: 0 159 0.0%
Branches: 0 602 0.0%

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