source: LMDZ6/branches/DYNAMICO-conv/libf/phylmd/rrtm/readaerosolstrato2_rrtm.F90 @ 3863

Last change on this file since 3863 was 3413, checked in by Laurent Fairhead, 6 years ago

Inclusion of Yann's latest (summer/fall 2018) modifications for
convergence of DYNAMICO/LMDZ physics
YM/LF

File size: 15.7 KB
RevLine 
[2536]1!
2! $Id: readaerosolstrato2_rrtm.F90 2526 2016-05-26 22:13:40Z oboucher $
3!
[3413]4
[2694]5SUBROUTINE readaerosolstrato2_rrtm(debut)
[2536]6
[2694]7    USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
[2536]8                        nf95_inq_varid, nf95_open
[2694]9    USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
[2536]10
11    USE phys_cal_mod, ONLY : mth_cur
[3413]12    USE mod_grid_phy_lmdz, ONLY: nbp_lon, nbp_lat, klon_glo, grid2dTo1d_glo, grid_type, unstructured
13    USE mod_phys_lmdz_mpi_data
14    USE mod_phys_lmdz_omp_data
[2536]15    USE mod_phys_lmdz_para
16    USE phys_state_var_mod
17    USE phys_local_var_mod
18    USE aero_mod
19    USE dimphy
[2694]20    USE YOERAD, ONLY : NLW
21    USE YOMCST
[3413]22    USE xios
[2694]23    IMPLICIT NONE
[2536]24
[2539]25    INCLUDE "clesphys.h"
[2536]26
27    CHARACTER (len = 80) :: abort_message
28    CHARACTER (LEN=20) :: modname = 'readaerosolstrato2'
29
30! Variable input
[2694]31    LOGICAL, INTENT(IN) ::  debut
[2536]32
33! Variables locales
[2694]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
[2744]44    INTEGER, SAVE :: mth_pre=1
45!$OMP THREADPRIVATE(mth_pre)
[2536]46
[2694]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
[2536]51!$OMP THREADPRIVATE(tau_aer_strat,piz_aer_strat,cg_aer_strat,taulw_aer_strat)
52
53! Champs reconstitues
[2694]54    REAL, ALLOCATABLE:: tauaerstrat(:, :, :, :)
55    REAL, ALLOCATABLE:: pizaerstrat(:, :, :, :)
56    REAL, ALLOCATABLE:: cgaerstrat(:, :, :, :)
57    REAL, ALLOCATABLE:: taulwaerstrat(:, :, :, :)
[2536]58
[2694]59    REAL, ALLOCATABLE:: tauaerstrat_mois(:, :, :, :)
60    REAL, ALLOCATABLE:: pizaerstrat_mois(:, :, :, :)
61    REAL, ALLOCATABLE:: cgaerstrat_mois(:, :, :, :)
62    REAL, ALLOCATABLE:: taulwaerstrat_mois(:, :, :, :)
[2536]63
[2694]64    REAL, ALLOCATABLE:: tauaerstrat_mois_glo(:, :, :)
65    REAL, ALLOCATABLE:: pizaerstrat_mois_glo(:, :, :)
66    REAL, ALLOCATABLE:: cgaerstrat_mois_glo(:, :, :)
67    REAL, ALLOCATABLE:: taulwaerstrat_mois_glo(:, :, :)
[3413]68    REAL, ALLOCATABLE:: tauaerstrat_mpi(:, :, :)
69    REAL, ALLOCATABLE:: pizaerstrat_mpi(:, :, :)
70    REAL, ALLOCATABLE:: cgaerstrat_mpi(:, :, :)
71    REAL, ALLOCATABLE:: taulwaerstrat_mpi(:, :, :)
[2536]72
73! For NetCDF:
[2694]74    INTEGER ncid_in  ! IDs for input files
75    INTEGER varid, ncerr
[2536]76
77!--------------------------------------------------------
78
[2539]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))
[2536]82
[2539]83    IF (.not.ALLOCATED(taulw_aer_strat)) ALLOCATE(taulw_aer_strat(klon,klev,NLW))
[2536]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
[2725]89      IF (is_mpi_root.AND.is_omp_root) THEN
[2536]90
91!--check mth_cur
[2725]92        IF (mth_cur.LT.1.OR.mth_cur.GT.12) THEN
93          print *,'probleme avec le mois dans readaerosolstrat =', mth_cur
94        ENDIF
[2536]95
96!--initialize n_lon as input data is 2D (lat-alt) only
[2725]97        n_lon = nbp_lon
[2536]98
99!--Starts with SW optical properties
100
[2725]101        CALL nf95_open("tauswstrat.2D.nc", nf90_nowrite, ncid_in)
[2536]102
[2725]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
[2536]110
[2725]111        CALL nf95_inq_varid(ncid_in, "LAT", varid)
112        CALL nf95_gw_var(ncid_in, varid, latitude)
113        n_lat = size(latitude)
[3413]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
[2725]121        ENDIF
[2536]122
[2725]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
[2536]130
[2725]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
[2536]139
[2725]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))
[2536]143
144!--reading stratospheric aerosol tau per layer
[2725]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
[2536]148
149!--reading stratospheric aerosol omega per layer
[2725]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
[2536]153
154!--reading stratospheric aerosol g per layer
[2725]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
[2536]158
[2725]159        CALL nf95_close(ncid_in)
[2536]160
[3413]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))
[2536]170!--select the correct month
171!--and copy into 1st longitude
[3413]172          tauaerstrat_mois(1,:,:,:) = tauaerstrat(:,:,:,mth_cur)
173          pizaerstrat_mois(1,:,:,:) = pizaerstrat(:,:,:,mth_cur)
174          cgaerstrat_mois(1,:,:,:)  = cgaerstrat(:,:,:,mth_cur)
[2536]175
176!--copy longitudes
[3413]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
[2536]182
183!---reduce to a klon_glo grid
[3413]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
[2536]190!--Now LW optical properties
191!
[3413]192
[2725]193        CALL nf95_open("taulwstrat.2D.nc", nf90_nowrite, ncid_in)
[2536]194
[2725]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
[2536]202
[2725]203        CALL nf95_inq_varid(ncid_in, "LAT", varid)
204        CALL nf95_gw_var(ncid_in, varid, latitude)
205        n_lat = size(latitude)
[3413]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
[2725]212        ENDIF
[3413]213       
[2725]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
[2536]221
[2725]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
[2536]230
[2725]231        ALLOCATE(taulwaerstrat(n_lat, n_lev, n_wav, n_month))
[2536]232
233!--reading stratospheric aerosol lw tau per layer
[2725]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
[2536]237
[2725]238        CALL nf95_close(ncid_in)
[2536]239
[3413]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
[2536]245!--select the correct month
246!--and copy into 1st longitude
[3413]247          taulwaerstrat_mois(1,:,:,:) = taulwaerstrat(:,:,:,mth_cur)
[2536]248!--copy longitudes
[3413]249          DO i=2, n_lon
250            taulwaerstrat_mois(i,:,:,:) = taulwaerstrat_mois(1,:,:,:)
251          ENDDO
[2536]252
253!---reduce to a klon_glo grid
[3413]254          DO band=1, NLW
255            CALL grid2dTo1d_glo(taulwaerstrat_mois(:,:,:,band),taulwaerstrat_mois_glo(:,:,band))
256          ENDDO
257        ENDIF
258     
[2725]259      ELSE !--proc other than mpi_root and omp_root
260           !--dummy allocation needed for debug mode
[2536]261
[2725]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
[3413]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
[2725]273      ENDIF !--is_mpi_root and is_omp_root
274
[2536]275!$OMP BARRIER
276
277!--keep memory of previous month
[2725]278      mth_pre=mth_cur
[2536]279
[3413]280      IF (grid_type==unstructured) THEN
281
282        IF (is_omp_master) THEN
283          ALLOCATE(tauaerstrat_mpi(klon_mpi, klev, NSW))
284          ALLOCATE(pizaerstrat_mpi(klon_mpi, klev, NSW))
285          ALLOCATE(cgaerstrat_mpi(klon_mpi, klev, NSW))       
286          ALLOCATE(taulwaerstrat_mpi(klon_mpi, klev, NLW))
287         
288          CALL xios_send_field("tauaerstrat_in",SPREAD(tauaerstrat(:,:,:,mth_cur),1,8))
289          CALL xios_recv_field("tauaerstrat_out",tauaerstrat_mpi)
290          CALL xios_send_field("pizaerstrat_in",SPREAD(pizaerstrat(:,:,:,mth_cur),1,8))
291          CALL xios_recv_field("pizaerstrat_out",pizaerstrat_mpi)
292          CALL xios_send_field("cgaerstrat_in",SPREAD(cgaerstrat(:,:,:,mth_cur),1,8))
293          CALL xios_recv_field("cgaerstrat_out",cgaerstrat_mpi)
294          CALL xios_send_field("taulwaerstrat_in",SPREAD(taulwaerstrat(:,:,:,mth_cur),1,8))
295          CALL xios_recv_field("taulwaerstrat_out",taulwaerstrat_mpi)
296        ELSE
297          ALLOCATE(tauaerstrat_mpi(0, 0, 0))
298          ALLOCATE(pizaerstrat_mpi(0, 0, 0))
299          ALLOCATE(cgaerstrat_mpi(0, 0, 0))       
300          ALLOCATE(taulwaerstrat_mpi(0, 0, 0))
301        ENDIF 
302       
303        CALL scatter_omp(tauaerstrat_mpi,tau_aer_strat)
304        CALL scatter_omp(pizaerstrat_mpi,piz_aer_strat)
305        CALL scatter_omp(cgaerstrat_mpi,cg_aer_strat)
306        CALL scatter_omp(taulwaerstrat_mpi,taulw_aer_strat)
307      ELSE 
308       
[2536]309!--scatter on all proc
[3413]310        CALL scatter(tauaerstrat_mois_glo,tau_aer_strat)
311        CALL scatter(pizaerstrat_mois_glo,piz_aer_strat)
312        CALL scatter(cgaerstrat_mois_glo,cg_aer_strat)
313        CALL scatter(taulwaerstrat_mois_glo,taulw_aer_strat)
[2536]314
[2725]315        DEALLOCATE(tauaerstrat_mois, pizaerstrat_mois, cgaerstrat_mois)
[3413]316        DEALLOCATE(taulwaerstrat_mois)
317      ENDIF
[2536]318
[3413]319      IF (is_mpi_root.AND.is_omp_root) THEN
320        DEALLOCATE(tauaerstrat, pizaerstrat, cgaerstrat,taulwaerstrat)
321      ENDIF
322     
[2725]323
[2536]324!$OMP BARRIER
325
326    ENDIF !--debut ou nouveau mois
327
328!--total vertical aod at the 5 SW wavelengths
329!--for now use band 3 AOD into all 5 wavelengths
[2550]330!--it is only a reasonable approximation for 550 nm (wave=2)
[2536]331    band=3
332    DO i=1, klon
333    DO k=1, klev
334      IF (stratomask(i,k).GT.0.999999) THEN
[2550]335        DO wave=1, nwave_sw
336          tausum_aero(i,wave,id_STRAT_phy)=tausum_aero(i,wave,id_STRAT_phy)+tau_aer_strat(i,k,band)
[2536]337        ENDDO
338      ENDIF
339    ENDDO
340    ENDDO
341
342!--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
[2539]343    DO band=1, NSW
[2725]344      WHERE (stratomask.GT.0.999999)
[2539]345!--anthropogenic aerosols bands 1 to NSW
[2725]346        cg_aero_sw_rrtm(:,:,2,band)  = ( cg_aero_sw_rrtm(:,:,2,band)*piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) + &
347                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
348                                    MAX( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
349                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
350        piz_aero_sw_rrtm(:,:,2,band) = ( piz_aero_sw_rrtm(:,:,2,band)*tau_aero_sw_rrtm(:,:,2,band) +                             &
351                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
352                                    MAX( tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band), 1.e-15 )
353        tau_aero_sw_rrtm(:,:,2,band)  = tau_aero_sw_rrtm(:,:,2,band) + tau_aer_strat(:,:,band)
[2539]354!--natural aerosols bands 1 to NSW
[2725]355        cg_aero_sw_rrtm(:,:,1,band)  = ( cg_aero_sw_rrtm(:,:,1,band)*piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) + &
356                                         cg_aer_strat(:,:,band)*piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /              &
357                                    MAX( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
358                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band), 1.e-15 )
359        piz_aero_sw_rrtm(:,:,1,band) = ( piz_aero_sw_rrtm(:,:,1,band)*tau_aero_sw_rrtm(:,:,1,band) +                             &
360                                         piz_aer_strat(:,:,band)*tau_aer_strat(:,:,band) ) /                                     &
361                                    MAX( tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band), 1.e-15 )
362        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band) + tau_aer_strat(:,:,band)
[2536]363!--no stratospheric aerosol in index 1 for these tests
[2725]364!        cg_aero_sw_rrtm(:,:,1,band)  =  cg_aero_sw_rrtm(:,:,1,band)
365!        piz_aero_sw_rrtm(:,:,1,band)  = piz_aero_sw_rrtm(:,:,1,band)
366!        tau_aero_sw_rrtm(:,:,1,band)  = tau_aero_sw_rrtm(:,:,1,band)
[2536]367    ENDWHERE
368    ENDDO
369
[2550]370!--total vertical aod at 10 um
371!--this is approximated from band 7 of RRTM
372    band=7
373    DO i=1, klon
374    DO k=1, klev
375      IF (stratomask(i,k).GT.0.999999) THEN
376        DO wave=1, nwave_lw
377          tausum_aero(i,nwave_sw+wave,id_STRAT_phy)=tausum_aero(i,nwave_sw+wave,id_STRAT_phy)+taulw_aer_strat(i,k,band)
378        ENDDO
379      ENDIF
380    ENDDO
381    ENDDO
382
[2539]383    DO band=1, NLW
[2725]384      WHERE (stratomask.GT.0.999999)
385        tau_aero_lw_rrtm(:,:,2,band)  = tau_aero_lw_rrtm(:,:,2,band) + taulw_aer_strat(:,:,band)
386        tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band) + taulw_aer_strat(:,:,band)
[2536]387!--no stratospheric aerosols in index 1 for these tests
388!    tau_aero_lw_rrtm(:,:,1,band)  = tau_aero_lw_rrtm(:,:,1,band)
[2725]389      ENDWHERE
[2536]390    ENDDO
391
392!--default SSA value if there is no aerosol
393!--to avoid 0 values that seems to cause some problem to RRTM
394    WHERE (tau_aero_sw_rrtm.LT.1.e-14)
395      piz_aero_sw_rrtm = 1.0
396    ENDWHERE
397
398!--in principle this should not be necessary
399!--as these variables have min values already but just in case
400!--put 1e-15 min value to both SW and LW AOD
401    tau_aero_sw_rrtm = MAX(tau_aero_sw_rrtm,1.e-15)
402    tau_aero_lw_rrtm = MAX(tau_aero_lw_rrtm,1.e-15)
403
[2694]404END SUBROUTINE readaerosolstrato2_rrtm
Note: See TracBrowser for help on using the repository browser.