source: LMDZ6/branches/Amaury_dev/libf/phylmd/readaerosolstrato.F90 @ 5151

Last change on this file since 5151 was 5144, checked in by abarral, 3 months ago

Put YOMCST.h into modules
  • Property copyright set to
    Name of program: LMDZ
    Creation date: 1984
    Version: LMDZ5
    License: CeCILL version 2
    Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539
    See the license file in the root directory
File size: 7.7 KB
Line 
1SUBROUTINE readaerosolstrato(debut)
2
3  USE netcdf95, ONLY: nf95_close, nf95_gw_var, nf95_inq_dimid, &
4          nf95_inq_varid, nf95_open
5  USE netcdf, ONLY: nf90_get_var, nf90_noerr, nf90_nowrite
6
7  USE phys_cal_mod, ONLY: mth_cur
8  USE lmdz_grid_phy, ONLY: nbp_lon, nbp_lat, klon_glo, &
9          grid2dto1d_glo, grid_type, unstructured
10  USE lmdz_phys_mpi_data, ONLY: is_mpi_root
11  USE lmdz_phys_omp_data, ONLY: is_omp_root
12  USE lmdz_phys_para
13  USE phys_state_var_mod
14  USE phys_local_var_mod
15  USE aero_mod
16  USE dimphy
17  USE lmdz_print_control, ONLY: prt_level, lunout
18  USE lmdz_xios
19  USE lmdz_abort_physic, ONLY: abort_physic
20  USE lmdz_yomcst
21
22  IMPLICIT NONE
23
24  ! Variable input
25  LOGICAL debut
26
27  ! Variables locales
28  INTEGER n_lat   ! number of latitudes in the input data
29  INTEGER n_lon   ! number of longitudes in the input data
30  INTEGER n_lev   ! number of levels in the input data
31  INTEGER n_month ! number of months in the input data
32  REAL, ALLOCATABLE :: latitude(:)
33  REAL, ALLOCATABLE :: longitude(:)
34  REAL, ALLOCATABLE :: time(:)
35  REAL, ALLOCATABLE :: lev(:)
36  INTEGER i, k, band, wave
37  INTEGER, save :: mth_pre = 1
38  !$OMP THREADPRIVATE(mth_pre)
39
40  REAL, ALLOCATABLE, DIMENSION(:, :), save :: tau_aer_strat
41  !$OMP THREADPRIVATE(tau_aer_strat)
42
43  ! Champs reconstitues
44  REAL, ALLOCATABLE :: tauaerstrat(:, :, :, :)
45  REAL, ALLOCATABLE :: tauaerstrat_mois(:, :, :)
46  REAL, ALLOCATABLE :: tauaerstrat_mois_glo(:, :)
47  REAL, ALLOCATABLE :: tau_aer_strat_mpi(:, :)
48
49  ! For NetCDF:
50  INTEGER ncid_in  ! IDs for input files
51  INTEGER varid, ncerr
52
53  ! Stratospheric aerosols optical properties
54  ! alpha_strat over the 2 bands is normalised by the 550 nm extinction coefficient
55  ! alpha_strat_wave is *not* normalised by the 550 nm extinction coefficient
56  REAL, DIMENSION(nbands) :: alpha_strat, piz_strat, cg_strat
57  data alpha_strat/0.9922547, 0.7114912 /
58  data piz_strat  /0.9999998, 0.99762493/
59  data cg_strat   /0.73107845, 0.73229635/
60  REAL, DIMENSION(nwave_sw) :: alpha_strat_wave
61  data alpha_strat_wave/3.36780953, 3.34667683, 3.20444202, 3.0293026, 2.82108808/
62
63  CHARACTER (len = 20) :: modname = 'readaerosolstrato'
64  CHARACTER (len = 80) :: abort_message
65
66  !--------------------------------------------------------
67
68  IF (.NOT.ALLOCATED(tau_aer_strat)) ALLOCATE(tau_aer_strat(klon, klev))
69
70  !--only read file if beginning of run or start of new month
71  IF (debut.OR.mth_cur/=mth_pre) THEN
72
73    !--only root reads
74    IF (is_mpi_root.AND.is_omp_root) THEN
75
76      IF (nbands/=2) THEN
77        abort_message = 'nbands doit etre egal a 2 dans readaerosolstrat'
78        CALL abort_physic(modname, abort_message, 1)
79      ENDIF
80
81      CALL nf95_open("taustrat.nc", nf90_nowrite, ncid_in)
82
83      CALL nf95_inq_varid(ncid_in, "LEV", varid)
84      CALL nf95_gw_var(ncid_in, varid, lev)
85      n_lev = size(lev)
86      IF (n_lev/=klev) THEN
87        abort_message = 'Le nombre de niveaux n est pas egal a klev'
88        CALL abort_physic(modname, abort_message, 1)
89      ENDIF
90
91      CALL nf95_inq_varid(ncid_in, "LAT", varid)
92      CALL nf95_gw_var(ncid_in, varid, latitude)
93      n_lat = size(latitude)
94      WRITE(lunout, *) 'LAT aerosol strato=', n_lat, latitude
95      IF (grid_type/=unstructured) THEN
96        IF (n_lat/=nbp_lat) THEN
97          abort_message = 'Le nombre de lat n est pas egal a nbp_lat'
98          CALL abort_physic(modname, abort_message, 1)
99        ENDIF
100      ENDIF
101
102      CALL nf95_inq_varid(ncid_in, "LON", varid)
103      CALL nf95_gw_var(ncid_in, varid, longitude)
104      n_lon = size(longitude)
105      IF (grid_type/=unstructured) THEN
106        WRITE(lunout, *) 'LON aerosol strato=', n_lon, longitude
107        IF (n_lon/=nbp_lon) THEN
108          abort_message = 'Le nombre de lon n est pas egal a nbp_lon'
109          CALL abort_physic(modname, abort_message, 1)
110        ENDIF
111      ENDIF
112
113      CALL nf95_inq_varid(ncid_in, "TIME", varid)
114      CALL nf95_gw_var(ncid_in, varid, time)
115      n_month = size(time)
116      WRITE(lunout, *) 'TIME aerosol strato=', n_month, time
117      IF (n_month/=12) THEN
118        abort_message = 'Le nombre de month n est pas egal a 12'
119        CALL abort_physic(modname, abort_message, 1)
120      ENDIF
121
122      IF (.NOT.ALLOCATED(tauaerstrat))          ALLOCATE(tauaerstrat(n_lon, n_lat, n_lev, n_month))
123      IF (.NOT.ALLOCATED(tauaerstrat_mois))     ALLOCATE(tauaerstrat_mois(n_lon, n_lat, n_lev))
124      IF (.NOT.ALLOCATED(tauaerstrat_mois_glo)) ALLOCATE(tauaerstrat_mois_glo(klon_glo, n_lev))
125
126      !--reading stratospheric AOD at 550 nm
127      CALL nf95_inq_varid(ncid_in, "TAUSTRAT", varid)
128      ncerr = nf90_get_var(ncid_in, varid, tauaerstrat)
129      WRITE(lunout, *) 'code erreur readaerosolstrato=', ncerr, varid
130
131      CALL nf95_close(ncid_in)
132
133      !---select the correct month
134      IF (mth_cur<1.OR.mth_cur>12) THEN
135        WRITE(lunout, *) 'probleme avec le mois dans readaerosolstrat =', mth_cur
136      ENDIF
137      tauaerstrat_mois(:, :, :) = tauaerstrat(:, :, :, mth_cur)
138
139      !---reduce to a klon_glo grid
140      CALL grid2dTo1d_glo(tauaerstrat_mois, tauaerstrat_mois_glo)
141
142    ELSE
143      ALLOCATE(tauaerstrat_mois(0, 0, 0))
144    ENDIF !--is_mpi_root and is_omp_root
145
146    !$OMP BARRIER
147
148    IF (grid_type==unstructured) THEN
149      IF (is_omp_master) THEN
150        CALL xios_send_field("taustrat_in", tauaerstrat_mois)
151        ALLOCATE(tau_aer_strat_mpi(klon_mpi, klev))
152        CALL xios_recv_field("taustrat_out", tau_aer_strat_mpi)
153      ELSE
154        ALLOCATE(tau_aer_strat_mpi(0, 0))
155      ENDIF
156      CALL scatter_omp(tau_aer_strat_mpi, tau_aer_strat)
157    ELSE
158      !--scatter on all proc
159      CALL scatter(tauaerstrat_mois_glo, tau_aer_strat)
160    ENDIF
161
162    !--keep memory of previous month
163    mth_pre = mth_cur
164
165    IF (is_mpi_root.AND.is_omp_root) THEN
166
167      DEALLOCATE(tauaerstrat)
168      DEALLOCATE(tauaerstrat_mois)
169      DEALLOCATE(tauaerstrat_mois_glo)
170
171    ENDIF !-is_mpi_root and is_omp_root
172
173    !$OMP BARRIER
174
175  ENDIF !--debut ou nouveau mois
176
177  !--total vertical aod at the 6 wavelengths
178  DO wave = 1, nwave_sw
179    DO k = 1, klev
180      tausum_aero(:, wave, id_STRAT_phy) = tausum_aero(:, wave, id_STRAT_phy) + tau_aer_strat(:, k) * alpha_strat_wave(wave) / alpha_strat_wave(2)
181    ENDDO
182  ENDDO
183
184  !--weighted average for cg, piz and tau, adding strat aerosols on top of tropospheric ones
185  DO band = 1, nbands
186    !--anthropogenic aerosols bands 1 and 2
187    cg_aero(:, :, 3, band) = (cg_aero(:, :, 3, band) * piz_aero(:, :, 3, band) * tau_aero(:, :, 3, band) + &
188            cg_strat(band) * piz_strat(band) * alpha_strat(band) * tau_aer_strat(:, :)) / &
189            MAX(piz_aero(:, :, 3, band) * tau_aero(:, :, 3, band) + &
190                    piz_strat(band) * alpha_strat(band) * tau_aer_strat(:, :), 1.e-15)
191    piz_aero(:, :, 3, band) = (piz_aero(:, :, 3, band) * tau_aero(:, :, 3, band) + &
192            piz_strat(band) * alpha_strat(band) * tau_aer_strat(:, :)) / &
193            MAX(tau_aero(:, :, 3, band) + alpha_strat(band) * tau_aer_strat(:, :), 1.e-15)
194    tau_aero(:, :, 3, band) = tau_aero(:, :, 3, band) + alpha_strat(band) * tau_aer_strat(:, :)
195    !--natural aerosols bands 1 and 2
196    cg_aero(:, :, 2, band) = (cg_aero(:, :, 2, band) * piz_aero(:, :, 2, band) * tau_aero(:, :, 2, band) + &
197            cg_strat(band) * piz_strat(band) * alpha_strat(band) * tau_aer_strat(:, :)) / &
198            MAX(piz_aero(:, :, 2, band) * tau_aero(:, :, 2, band) + &
199                    piz_strat(band) * alpha_strat(band) * tau_aer_strat(:, :), 1.e-15)
200    piz_aero(:, :, 2, band) = (piz_aero(:, :, 2, band) * tau_aero(:, :, 2, band) + &
201            piz_strat(band) * alpha_strat(band) * tau_aer_strat(:, :)) / &
202            MAX(tau_aero(:, :, 2, band) + alpha_strat(band) * tau_aer_strat(:, :), 1.e-15)
203    tau_aero(:, :, 2, band) = tau_aero(:, :, 2, band) + alpha_strat(band) * tau_aer_strat(:, :)
204  ENDDO
205
206END SUBROUTINE  readaerosolstrato
Note: See TracBrowser for help on using the repository browser.