source: LMDZ6/trunk/libf/phylmd/ecrad/lmdz/readaerosolstrato_ecrad.F90 @ 5151

Last change on this file since 5151 was 4853, checked in by idelkadi, 8 months ago

Ecrad update in LMDZ / Implementation of Ecrad double call in LMDZ

  • version 1.6.1 (commit 210d7911380f53a788c3cad73b3cf9b4e022ef87)
  • interface routines between lmdz and ecrad grouped in a new "lmdz" directory
  • double call method redesigned so as to go through the Ecrad initialization and configuration part only once for the entire simulation
  • clean-up in the read.F routine: delete unitules arguments
  • modification of default gas model in namelist (default: ECCKD)
File size: 24.6 KB
RevLine 
[4853]1!
2! $Id: readaerosolstrato_ecrad.F90 tlurton $
3!
4SUBROUTINE readaerosolstrato_ecrad(config, 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 physiq_mod, ONLY : ok_volcan
18    USE aero_mod
19    USE dimphy
20!    USE YOERAD
21    USE radiation_config, ONLY : config_type
22    USE YOMCST
23#ifdef CPP_XIOS
24    USE xios
25#endif
26
27    IMPLICIT NONE
28
29    INCLUDE "clesphys.h"
30
31    CHARACTER (len = 80) :: abort_message
32    CHARACTER (LEN=20) :: modname = 'readaerosolstrato_ecrad'
33
34! Variable input
35    LOGICAL, INTENT(IN) ::  debut
36    LOGICAL, INTENT(IN) ::  ok_volcan !activate volcanic diags
37
38! Variables locales
39    INTEGER n_lat   ! number of latitudes in the input data
40    INTEGER n_lon   ! number of longitudes
41    INTEGER n_lev   ! number of levels in the input data
42    INTEGER n_month ! number of months in the input data
43    INTEGER n_wav   ! number of wavelengths in the input data
44    REAL, ALLOCATABLE:: latitude(:)
45    REAL, ALLOCATABLE:: time(:)
46    REAL, ALLOCATABLE:: lev(:)
47    REAL, ALLOCATABLE:: wav(:)
48    INTEGER i,k,wave,band
49    INTEGER, SAVE :: mth_pre=1
50!$OMP THREADPRIVATE(mth_pre)
51
52    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: tau_aer_strat
53    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: piz_aer_strat
54    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cg_aer_strat
55    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: taulw_aer_strat
56! add ThL
57    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: pizlw_aer_strat
58    REAL, ALLOCATABLE, DIMENSION(:,:,:), SAVE :: cglw_aer_strat
59! end add ThL
60! mod ThL
61!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat,pizlw_aer_strat,cglw_aer_strat)
62! end mod ThL
63
64! Champs reconstitues
65    REAL, ALLOCATABLE:: tauaerstrat(:, :, :, :)
66    REAL, ALLOCATABLE:: pizaerstrat(:, :, :, :)
67    REAL, ALLOCATABLE:: cgaerstrat(:, :, :, :)
68    REAL, ALLOCATABLE:: taulwaerstrat(:, :, :, :)
69! add ThL
70    REAL, ALLOCATABLE:: pizlwaerstrat(:, :, :, :)
71    REAL, ALLOCATABLE:: cglwaerstrat(:, :, :, :)
72! end add ThL
73
74    REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :, :)
75    REAL, ALLOCATABLE:: pizaerstrat_mois(:, :, :, :)
76    REAL, ALLOCATABLE:: cgaerstrat_mois(:, :, :, :)
77    REAL, ALLOCATABLE:: taulwaerstrat_mois(:, :, :, :)
78! add ThL
79    REAL, ALLOCATABLE:: pizlwaerstrat_mois(:, :, :, :)
80    REAL, ALLOCATABLE:: cglwaerstrat_mois(:, :, :, :)
81! end add ThL
82
83    REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :, :)
84    REAL, ALLOCATABLE:: pizaerstrat_mois_glo(:, :, :)
85    REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
86    REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
87! add ThL
88    REAL, ALLOCATABLE:: pizlwaerstrat_mois_glo(:, :, :)
89    REAL, ALLOCATABLE:: cglwaerstrat_mois_glo(:, :, :)
90! end add ThL
91    REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
92    REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
93    REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
94    REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
95! add ThL
96    REAL, ALLOCATABLE:: pizlwaerstrat_mpi(:, :, :)
97    REAL, ALLOCATABLE:: cglwaerstrat_mpi(:, :, :)
98! end add ThL
99
100! For NetCDF:
101    INTEGER ncid_in  ! IDs for input files
102    INTEGER varid, ncerr
103
104    TYPE(config_type), INTENT(IN) :: config
105
106!--------------------------------------------------------
107
108    IF (.not.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon,klev,config%n_bands_sw))
109    IF (.not.ALLOCATED(piz_aer_strat)) ALLOCATE(piz_aer_strat(klon,klev,config%n_bands_sw))
110    IF (.not.ALLOCATED(cg_aer_strat))  ALLOCATE(cg_aer_strat(klon,klev,config%n_bands_sw))
111
112    IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,config%n_bands_lw))
113! add ThL   
114    IF (.not.ALLOCATED(pizlw_aer_strat)) ALLOCATE(pizlw_aer_strat(klon,klev,config%n_bands_lw))
115    IF (.not.ALLOCATED(cglw_aer_strat)) ALLOCATE(cglw_aer_strat(klon,klev,config%n_bands_lw))
116! end add ThL
117
118!--we only read monthly strat aerosol data
119    IF (debut.OR.mth_cur.NE.mth_pre) THEN
120
121!--only root reads the data
122      IF (is_mpi_root.AND.is_omp_root) THEN
123
124!--check mth_cur
125        IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
126          print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
127        ENDIF
128
129!--initialize n_lon as input data is 2D (lat-alt) only
130        n_lon = nbp_lon
131
132!--Starts with SW optical properties
133
134        CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
135
136        CALL nf95_inq_varid(ncid_in, "LEV", varid)
137        CALL nf95_gw_var(ncid_in, varid, lev)
138        n_lev = size(lev)
139        IF (n_lev.NE.klev) THEN
140           abort_message='Le nombre de niveaux n est pas egal a klev'
141           CALL abort_physic(modname,abort_message,1)
142        ENDIF
143
144        CALL nf95_inq_varid(ncid_in, "LAT", varid)
145        CALL nf95_gw_var(ncid_in, varid, latitude)
146        n_lat = size(latitude)
147
148        IF (grid_type/=unstructured) THEN
149           IF (n_lat.NE.nbp_lat) THEN
150             print *, 'latitude=', n_lat, nbp_lat
151             abort_message='Le nombre de lat n est pas egal a nbp_lat'
152             CALL abort_physic(modname,abort_message,1)
153           ENDIF
154        ENDIF
155
156        CALL nf95_inq_varid(ncid_in, "TIME", varid)
157        CALL nf95_gw_var(ncid_in, varid, time)
158        n_month = size(time)
159        IF (n_month.NE.12) THEN
160           abort_message='Le nombre de month n est pas egal a 12'
161           CALL abort_physic(modname,abort_message,1)
162        ENDIF
163
164        CALL nf95_inq_varid(ncid_in, "WAV", varid)
165        CALL nf95_gw_var(ncid_in, varid, wav)
166        n_wav = size(wav)
167        print *, 'WAV aerosol strato=', n_wav, wav
168!        IF (n_wav.NE.config%n_bands_sw) THEN
169!           abort_message='Le nombre de wav n est pas egal a config%n_bands_sw'
170!           CALL abort_physic(modname,abort_message,1)
171!        ENDIF
172
173        ALLOCATE(tauaerstrat(n_lat, n_lev, n_wav, n_month))
174        ALLOCATE(pizaerstrat(n_lat, n_lev, n_wav, n_month))
175        ALLOCATE(cgaerstrat(n_lat, n_lev, n_wav, n_month))
176
177!--reading stratospheric aerosol tau per layer
178        CALL nf95_inq_varid(ncid_in, "TAU_SUN", varid)
179        ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
180        print *,'code erreur readaerosolstrato=', ncerr, varid
181
182!--reading stratospheric aerosol omega per layer
183        CALL nf95_inq_varid(ncid_in, "OME_SUN", varid)
184        ncerr = nf90_get_var(ncid_in, varid, pizaerstrat)
185        print *,'code erreur readaerosolstrato=', ncerr, varid
186
187!--reading stratospheric aerosol g per layer
188        CALL nf95_inq_varid(ncid_in, "GGG_SUN", varid)
189        ncerr = nf90_get_var(ncid_in, varid, cgaerstrat)
190        print *,'code erreur readaerosolstrato sw=', ncerr, varid
191
192        CALL nf95_close(ncid_in)
193
194       
195        IF (grid_type/=unstructured) THEN
196          ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
197          ALLOCATE(pizaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
198          ALLOCATE(cgaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
199
200          ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev, n_wav))
201          ALLOCATE(pizaerstrat_mois_glo(klon_glo, n_lev, n_wav))
202          ALLOCATE(cgaerstrat_mois_glo(klon_glo, n_lev, n_wav))
203!--select the correct month
204!--and copy into 1st longitude
205          tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
206          pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
207          cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
208
209!--copy longitudes
210          DO i=2, n_lon
211           tauaerstrat_mois(i,:,:,:) = tauaerstrat_mois(1,:,:,:)
212           pizaerstrat_mois(i,:,:,:) = pizaerstrat_mois(1,:,:,:)
213           cgaerstrat_mois(i,:,:,:)  = cgaerstrat_mois(1,:,:,:)
214          ENDDO
215
216!---reduce to a klon_glo grid
217          DO band=1, config%n_bands_sw
218            CALL grid2dTo1d_glo(tauaerstrat_mois(:,:,:,band),tauaerstrat_mois_glo(:,:,band))
219            CALL grid2dTo1d_glo(pizaerstrat_mois(:,:,:,band),pizaerstrat_mois_glo(:,:,band))
220            CALL grid2dTo1d_glo(cgaerstrat_mois(:,:,:,band),cgaerstrat_mois_glo(:,:,band))
221          ENDDO
222        ENDIF
223!--Now LW optical properties
224!
225
226        CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
227
228        CALL nf95_inq_varid(ncid_in, "LEV", varid)
229        CALL nf95_gw_var(ncid_in, varid, lev)
230        n_lev = size(lev)
231        IF (n_lev.NE.klev) THEN
232           abort_message='Le nombre de niveaux n est pas egal a klev'
233           CALL abort_physic(modname,abort_message,1)
234        ENDIF
235
236        CALL nf95_inq_varid(ncid_in, "LAT", varid)
237        CALL nf95_gw_var(ncid_in, varid, latitude)
238        n_lat = size(latitude)
239
240        IF (grid_type/=unstructured) THEN
241          IF (n_lat.NE.nbp_lat) THEN
242             abort_message='Le nombre de lat n est pas egal a nbp_lat'
243             CALL abort_physic(modname,abort_message,1)
244          ENDIF
245        ENDIF
246       
247        CALL nf95_inq_varid(ncid_in, "TIME", varid)
248        CALL nf95_gw_var(ncid_in, varid, time)
249        n_month = size(time)
250        IF (n_month.NE.12) THEN
251           abort_message='Le nombre de month n est pas egal a 12'
252           CALL abort_physic(modname,abort_message,1)
253        ENDIF
254
255        CALL nf95_inq_varid(ncid_in, "WAV", varid)
256        CALL nf95_gw_var(ncid_in, varid, wav)
257        n_wav = size(wav)
258        print *, 'WAV aerosol strato=', n_wav, wav
259!        IF (n_wav.NE.config%n_bands_lw) THEN
260!           abort_message='Le nombre de wav n est pas egal a config%n_bands_lw'
261!           CALL abort_physic(modname,abort_message,1)
262!        ENDIF
263
264        ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
265! add ThL
266        ALLOCATE(pizlwaerstrat(n_lat, n_lev, n_wav, n_month))
267        ALLOCATE(cglwaerstrat(n_lat, n_lev, n_wav, n_month))
268! end add ThL
269
270!--reading stratospheric aerosol lw tau per layer
271        CALL nf95_inq_varid(ncid_in, "TAU_EAR", varid)
272        ncerr = nf90_get_var(ncid_in, varid, taulwaerstrat)
273        print *,'code erreur readaerosolstrato lw=', ncerr, varid
274
275! add/mod ThL
276!--with ECRAd, we also read omega and g:
277       !--reading stratospheric aerosol omega per layer
278        CALL nf95_inq_varid(ncid_in, "OME_EAR", varid)
279        ncerr = nf90_get_var(ncid_in, varid, pizlwaerstrat)
280        print *,'code erreur readaerosolstrato lw=', ncerr, varid
281
282!--reading stratospheric aerosol g per layer
283        CALL nf95_inq_varid(ncid_in, "GGG_EAR", varid)
284        ncerr = nf90_get_var(ncid_in, varid, cglwaerstrat)
285        print *,'code erreur readaerosolstrato lw=', ncerr, varid
286
287        CALL nf95_close(ncid_in)
288
289        IF (grid_type/=unstructured) THEN
290
291          ALLOCATE(taulwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
292          ALLOCATE(pizlwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
293          ALLOCATE(cglwaerstrat_mois(n_lon, n_lat, n_lev, n_wav))
294
295          ALLOCATE(taulwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
296          ALLOCATE(pizlwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
297          ALLOCATE(cglwaerstrat_mois_glo(klon_glo, n_lev, n_wav))
298!--select the correct month
299!--and copy into 1st longitude
300          taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
301          pizlwaerstrat_mois(1,:,:,:) = pizlwaerstrat(:,:,:,mth_cur)
302          cglwaerstrat_mois(1,:,:,:) = cglwaerstrat(:,:,:,mth_cur)
303!--copy longitudes
304          DO i=2, n_lon
305            taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
306            pizlwaerstrat_mois(i,:,:,:) = pizlwaerstrat_mois(1,:,:,:)
307            cglwaerstrat_mois(i,:,:,:) = cglwaerstrat_mois(1,:,:,:)
308          ENDDO
309
310!---reduce to a klon_glo grid
311          DO band=1, config%n_bands_lw
312            CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
313            CALL grid2dTo1d_glo(pizlwaerstrat_mois(:,:,:,band),pizlwaerstrat_mois_glo(:,:,band))
314            CALL grid2dTo1d_glo(cglwaerstrat_mois(:,:,:,band),cglwaerstrat_mois_glo(:,:,band))
315          ENDDO
316        ENDIF
317     
318      ELSE !--proc other than mpi_root and omp_root
319           !--dummy allocation needed for debug mode
320
321        ALLOCATE(tauaerstrat_mois_glo(1,1,1))
322        ALLOCATE(pizaerstrat_mois_glo(1,1,1))
323        ALLOCATE(cgaerstrat_mois_glo(1,1,1))
324        ALLOCATE(taulwaerstrat_mois_glo(1,1,1))
325        ALLOCATE(pizlwaerstrat_mois_glo(1,1,1))
326        ALLOCATE(cglwaerstrat_mois_glo(1,1,1))
327
328        ALLOCATE(tauaerstrat(0,0,0,12))
329        ALLOCATE(pizaerstrat(0,0,0,12))
330        ALLOCATE(cgaerstrat(0,0,0,12))
331        ALLOCATE(taulwaerstrat(0,0,0,12))
332        ALLOCATE(pizlwaerstrat(0,0,0,12))
333        ALLOCATE(cglwaerstrat(0,0,0,12))
334! end add/mod ThL
335
336      ENDIF !--is_mpi_root and is_omp_root
337
338!$OMP BARRIER
339
340!--keep memory of previous month
341      mth_pre=mth_cur
342
343      IF (grid_type==unstructured) THEN
344
345#ifdef CPP_XIOS
346
347        IF (is_omp_master) THEN
348          ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, config%n_bands_sw))
349          ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, config%n_bands_sw))
350          ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, config%n_bands_sw))       
351          ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, config%n_bands_lw))
352! add ThL
353          ALLOCATE(pizlwaerstrat_mpi(klon_mpi, klev, config%n_bands_lw))
354          ALLOCATE(cglwaerstrat_mpi(klon_mpi, klev, config%n_bands_lw))
355! end add ThL
356
357          CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
358          CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
359          CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
360          CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
361          CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
362          CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
363          CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
364          CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
365! add ThL
366          CALL xios_send_field("pizlwaerstrat_in",SPREAD(pizlwaerstrat(:,:,:,mth_cur),1,8))
367          CALL xios_recv_field("pizlwaerstrat_out",pizlwaerstrat_mpi)
368          CALL xios_send_field("cglwaerstrat_in",SPREAD(cglwaerstrat(:,:,:,mth_cur),1,8))
369          CALL xios_recv_field("cglwaerstrat_out",cglwaerstrat_mpi)
370! end add ThL
371        ELSE
372          ALLOCATE(tauaerstrat_mpi(0, 0, 0))
373          ALLOCATE(pizaerstrat_mpi(0, 0, 0))
374          ALLOCATE(cgaerstrat_mpi(0, 0, 0))       
375          ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
376! add ThL
377          ALLOCATE(pizlwaerstrat_mpi(0, 0, 0))
378          ALLOCATE(cglwaerstrat_mpi(0, 0, 0))
379! end add ThL
380        ENDIF 
381       
382        CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
383        CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
384        CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
385        CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
386! add ThL
387        CALL scatter_omp(pizlwaerstrat_mpi,pizlw_aer_strat)
388        CALL scatter_omp(cglwaerstrat_mpi,cglw_aer_strat)
389! end add ThL
390#endif
391      ELSE 
392       
393!--scatter on all proc
394        CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
395        CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
396        CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
397        CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
398! add ThL
399        CALL scatter(pizlwaerstrat_mois_glo,pizlw_aer_strat)
400        CALL scatter(cglwaerstrat_mois_glo,cglw_aer_strat)
401! end add ThL
402        IF (is_mpi_root.AND.is_omp_root)  DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois, taulwaerstrat_mois, pizlwaerstrat_mois, cglwaerstrat_mois) ! mod ThL
403
404      ENDIF
405
406      IF (is_mpi_root.AND.is_omp_root) THEN
407        DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat, taulwaerstrat, pizlwaerstrat, cglwaerstrat) ! mod ThL
408      ENDIF
409     
410
411!$OMP BARRIER
412
413    ENDIF !--debut ou nouveau mois
414
415!--total vertical aod at the SW wavelengths
416!--for now use ECRAD band #10 AOD into all wavelengths
417!--it is only a reasonable approximation for 550 nm
418!--WAV(10)=0.53 µm
419!--<!!!> This is considering that wavebands are in their original order,
420!--and have NOT been re-ordered monotonously, i.e. with
421!--WAV(1) = 3.46 µm
422!--WAV(2) = 2.79 µm
423!--(...)
424!--WAV(last-1) = 0.23 µm
425!--WAV(last)= 8.02 µm
426!--If bands were to be re-ordered, i.e. last band be put in 1st position,
427!--then we should consider band=11 to fetch WAV=0.53 µm.
428    band=10
429    DO i=1, klon
430    DO k=1, klev
431      IF (stratomask(i,k).GT.0.999999) THEN
432        DO wave=1, config%n_bands_sw
433          tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
434        ENDDO
435      ENDIF
436    ENDDO
437    ENDDO
438
439    IF (.NOT. ok_volcan) THEN
440!
441!--this is the default case
442!--stratospheric aerosols are added to both index 2 and 1 for double radiation calls
443!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
444    DO band=1, config%n_bands_sw
445      WHERE (stratomask.GT.0.999999)
446!--strat aerosols are added to index 2: natural and anthropogenic aerosols for bands 1 to config%n_bands_sw
447        cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
448                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
449                                    MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
450                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
451        piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
452                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
453                                    MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
454        tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
455!--strat aerosols are added to index 1: natural aerosols only for bands 1 to config%n_bands_sw
456        cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
457                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
458                                    MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
459                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
460        piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
461                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
462                                    MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
463        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
464    ENDWHERE
465    ENDDO
466!
467    ELSE
468!
469!--this is the VOLMIP case
470!--stratospheric aerosols are only added to index 2 in this case
471!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
472    DO band=1, config%n_bands_sw
473      WHERE (stratomask.GT.0.999999)
474!--strat aerosols are added to index 2 : natural and anthropogenic aerosols for bands 1 to config%n_bands_sw
475        cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
476                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
477                                    MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
478                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
479        piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
480                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
481                                    MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
482        tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
483     ENDWHERE
484  ENDDO
485  ENDIF
486
487!--total vertical aod at 10 um
488!--this is approximated from band 7 of ECRAD (9.73 µm)
489    band=7
490    DO i=1, klon
491    DO k=1, klev
492      IF (stratomask(i,k).GT.0.999999) THEN
493        DO wave=1, config%n_bands_lw
494          tausum_aero(i,config%n_bands_sw+wave,id_STRAT_phy)=tausum_aero(i,config%n_bands_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
495        ENDDO
496      ENDIF
497    ENDDO
498    ENDDO
499
500    IF (.NOT. ok_volcan) THEN
501!--this is the default case
502!--stratospheric aerosols are added to both index 2 and 1
503    DO band=1, config%n_bands_lw
504      WHERE (stratomask.GT.0.999999)
505! add ThL
506        cg_aero_lw_rrtm(:,:,2,band)  = ( cg_aero_lw_rrtm(:,:,2,band)*piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) + &
507                                         cglw_aer_strat(:,:,band)*pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) /        &
508                                    MAX( piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) +                             &
509                                         pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band), 1.e-15 )
510        piz_aero_lw_rrtm(:,:,2,band) = ( piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) +                             &
511                                         pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) /                                 &
512                                    MAX( tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band), 1.e-15 )
513        cg_aero_lw_rrtm(:,:,1,band)  = ( cg_aero_lw_rrtm(:,:,1,band)*piz_aero_lw_rrtm(:,:,1,band)*tau_aero_lw_rrtm(:,:,1,band) + &
514                                         cglw_aer_strat(:,:,band)*pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) /        &
515                                    MAX( piz_aero_lw_rrtm(:,:,1,band)*tau_aero_lw_rrtm(:,:,1,band) +                             &
516                                         pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band), 1.e-15 )
517        piz_aero_lw_rrtm(:,:,1,band) = ( piz_aero_lw_rrtm(:,:,1,band)*tau_aero_lw_rrtm(:,:,1,band) +                             &
518                                         pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) /                                 &
519                                    MAX( tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band), 1.e-15 )
520! end add ThL
521        tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
522        tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
523      ENDWHERE
524    ENDDO
525!
526    ELSE
527!
528! mod ThL
529!--this is the VOLMIP case
530    DO band=1, config%n_bands_lw
531      WHERE (stratomask.GT.0.999999)
532!--stratospheric aerosols are not added to index 1
533!--and we copy index 2 in index 1 because we want the same dust aerosol LW properties as above
534        cg_aero_lw_rrtm(:,:,1,band)  = cg_aero_lw_rrtm(:,:,2,band)
535        piz_aero_lw_rrtm(:,:,1,band) = piz_aero_lw_rrtm(:,:,2,band)
536        tau_aero_lw_rrtm(:,:,1,band) = tau_aero_lw_rrtm(:,:,2,band)
537
538!--stratospheric aerosols are only added to index 2
539        cg_aero_lw_rrtm(:,:,2,band)  = (cg_aero_lw_rrtm(:,:,2,band)*piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) + &
540                                        cglw_aer_strat(:,:,band)*pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band)) /         &
541                                        MAX( piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) +                        &
542                                        pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band), 1.e-15 )
543        piz_aero_lw_rrtm(:,:,2,band) = (piz_aero_lw_rrtm(:,:,2,band)*tau_aero_lw_rrtm(:,:,2,band) +                             &
544                                        pizlw_aer_strat(:,:,band)*taulw_aer_strat(:,:,band) ) /                                 &
545                                        MAX( tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band), 1.e-15 )
546        tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
547      ENDWHERE
548    ENDDO
549    ENDIF
550! end mod ThL
551
552!--default SSA value if there is no aerosol
553!--to avoid 0 values that seems to cause some problem to RRTM
554    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
555      piz_aero_sw_rrtm = 1.0
556    ENDWHERE
557
558!--in principle this should not be necessary
559!--as these variables have min values already but just in case
560!--put 1e-15 min value to both SW and LW AOD
561    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
562    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
563
564END SUBROUTINE readaerosolstrato_ecrad
Note: See TracBrowser for help on using the repository browser.