1 |
|
|
! $Id: readaerosol.F90 3436 2019-01-22 16:26:21Z emillour $ |
2 |
|
|
! |
3 |
|
|
MODULE readaerosol_mod |
4 |
|
|
|
5 |
|
|
REAL, SAVE :: not_valid=-333. |
6 |
|
|
|
7 |
|
|
INTEGER, SAVE :: nbp_lon_src |
8 |
|
|
!$OMP THREADPRIVATE(nbp_lon_src) |
9 |
|
|
INTEGER, SAVE :: nbp_lat_src |
10 |
|
|
!$OMP THREADPRIVATE(nbp_lat_src) |
11 |
|
|
REAL, ALLOCATABLE, SAVE :: psurf_interp(:,:) |
12 |
|
|
|
13 |
|
|
CONTAINS |
14 |
|
|
|
15 |
|
|
SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
16 |
|
|
|
17 |
|
|
!**************************************************************************************** |
18 |
|
|
! This routine will read the aersosol from file. |
19 |
|
|
! |
20 |
|
|
! Read a year data with get_aero_fromfile depending on aer_type : |
21 |
|
|
! - actuel : read year 1980 |
22 |
|
|
! - preind : read natural data |
23 |
|
|
! - scenario : read one or two years and do eventually linare time interpolation |
24 |
|
|
! |
25 |
|
|
! Return pointer, pt_out, to the year read or result from interpolation |
26 |
|
|
!**************************************************************************************** |
27 |
|
|
USE dimphy |
28 |
|
|
USE print_control_mod, ONLY: lunout |
29 |
|
|
|
30 |
|
|
IMPLICIT NONE |
31 |
|
|
|
32 |
|
|
! Input arguments |
33 |
|
|
CHARACTER(len=7), INTENT(IN) :: name_aero |
34 |
|
|
CHARACTER(len=*), INTENT(IN) :: type ! actuel, annuel, scenario or preind |
35 |
|
|
CHARACTER(len=8), INTENT(IN) :: filename |
36 |
|
|
INTEGER, INTENT(IN) :: iyr_in |
37 |
|
|
|
38 |
|
|
! Output |
39 |
|
|
INTEGER, INTENT(OUT) :: klev_src |
40 |
|
|
REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels |
41 |
|
|
REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels |
42 |
|
|
REAL, POINTER, DIMENSION(:,:,:) :: pt_out ! The massvar distributions, DIMENSION(klon, klev_src, 12) |
43 |
|
|
REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf ! Surface pression for 12 months |
44 |
|
|
REAL, DIMENSION(klon,12), INTENT(OUT) :: load ! Aerosol mass load in each column for 12 months |
45 |
|
|
|
46 |
|
|
! Local variables |
47 |
|
|
CHARACTER(len=4) :: cyear |
48 |
|
|
REAL, POINTER, DIMENSION(:,:,:) :: pt_2 |
49 |
|
|
REAL, DIMENSION(klon,12) :: psurf2, load2 |
50 |
|
|
INTEGER :: iyr1, iyr2, klev_src2 |
51 |
|
|
INTEGER :: it, k, i |
52 |
|
|
LOGICAL, PARAMETER :: lonlyone=.FALSE. |
53 |
|
|
|
54 |
|
|
!**************************************************************************************** |
55 |
|
|
! Read data depending on aer_type |
56 |
|
|
! |
57 |
|
|
!**************************************************************************************** |
58 |
|
|
|
59 |
|
|
IF (type == 'actuel') THEN |
60 |
|
|
! Read and return data for year 1980 |
61 |
|
|
!**************************************************************************************** |
62 |
|
|
cyear='1980' |
63 |
|
|
! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
64 |
|
|
! pt_out has dimensions (klon, klev_src, 12) |
65 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
66 |
|
|
|
67 |
|
|
|
68 |
|
|
ELSE IF (type == 'preind') THEN |
69 |
|
|
! Read and return data from file with suffix .nat |
70 |
|
|
!**************************************************************************************** |
71 |
|
|
cyear='.nat' |
72 |
|
|
! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
73 |
|
|
! pt_out has dimensions (klon, klev_src, 12) |
74 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
75 |
|
|
|
76 |
|
|
ELSE IF (type == 'annuel') THEN |
77 |
|
|
! Read and return data from scenario annual files |
78 |
|
|
!**************************************************************************************** |
79 |
|
|
WRITE(cyear,'(I4)') iyr_in |
80 |
|
|
WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear |
81 |
|
|
! get_aero_fromfile returns pt_out allocated and initialized with data for nbr_tsteps month |
82 |
|
|
! pt_out has dimensions (klon, klev_src, 12) |
83 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
84 |
|
|
|
85 |
|
|
ELSE IF (type == 'scenario') THEN |
86 |
|
|
! Read data depending on actual year and interpolate if necessary |
87 |
|
|
!**************************************************************************************** |
88 |
|
|
IF (iyr_in .LT. 1850) THEN |
89 |
|
|
cyear='.nat' |
90 |
|
|
WRITE(lunout,*) 'get_aero 1 iyr_in=', iyr_in,' ',cyear |
91 |
|
|
! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
92 |
|
|
! pt_out has dimensions (klon, klev_src, 12) |
93 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
94 |
|
|
|
95 |
|
|
ELSE IF (iyr_in .GE. 2100) THEN |
96 |
|
|
cyear='2100' |
97 |
|
|
WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,' ',cyear |
98 |
|
|
! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
99 |
|
|
! pt_out has dimensions (klon, klev_src, 12) |
100 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
101 |
|
|
|
102 |
|
|
ELSE |
103 |
|
|
! Read data from 2 decades and interpolate to actual year |
104 |
|
|
! a) from actual 10-yr-period |
105 |
|
|
IF (iyr_in.LT.1900) THEN |
106 |
|
|
iyr1 = 1850 |
107 |
|
|
iyr2 = 1900 |
108 |
|
|
ELSE IF (iyr_in.GE.1900.AND.iyr_in.LT.1920) THEN |
109 |
|
|
iyr1 = 1900 |
110 |
|
|
iyr2 = 1920 |
111 |
|
|
ELSE |
112 |
|
|
iyr1 = INT(iyr_in/10)*10 |
113 |
|
|
iyr2 = INT(1+iyr_in/10)*10 |
114 |
|
|
ENDIF |
115 |
|
|
|
116 |
|
|
WRITE(cyear,'(I4)') iyr1 |
117 |
|
|
WRITE(lunout,*) 'get_aero 3 iyr_in=', iyr_in,' ',cyear |
118 |
|
|
! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month |
119 |
|
|
! pt_out has dimensions (klon, klev_src, 12) |
120 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load) |
121 |
|
|
|
122 |
|
|
! If to read two decades: |
123 |
|
|
IF (.NOT.lonlyone) THEN |
124 |
|
|
|
125 |
|
|
! b) from the next following one |
126 |
|
|
WRITE(cyear,'(I4)') iyr2 |
127 |
|
|
WRITE(lunout,*) 'get_aero 4 iyr_in=', iyr_in,' ',cyear |
128 |
|
|
|
129 |
|
|
NULLIFY(pt_2) |
130 |
|
|
! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month |
131 |
|
|
! pt_2 has dimensions (klon, klev_src, 12) |
132 |
|
|
CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, pt_2, psurf2, load2) |
133 |
|
|
! Test for same number of vertical levels |
134 |
|
|
IF (klev_src /= klev_src2) THEN |
135 |
|
|
WRITE(lunout,*) 'Two aerosols files with different number of vertical levels is not allowded' |
136 |
|
|
CALL abort_physic('readaersosol','Error in number of vertical levels',1) |
137 |
|
|
END IF |
138 |
|
|
|
139 |
|
|
! Linare interpolate to the actual year: |
140 |
|
|
DO it=1,12 |
141 |
|
|
DO k=1,klev_src |
142 |
|
|
DO i = 1, klon |
143 |
|
|
pt_out(i,k,it) = & |
144 |
|
|
pt_out(i,k,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * & |
145 |
|
|
(pt_out(i,k,it) - pt_2(i,k,it)) |
146 |
|
|
END DO |
147 |
|
|
END DO |
148 |
|
|
|
149 |
|
|
DO i = 1, klon |
150 |
|
|
psurf(i,it) = & |
151 |
|
|
psurf(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * & |
152 |
|
|
(psurf(i,it) - psurf2(i,it)) |
153 |
|
|
|
154 |
|
|
load(i,it) = & |
155 |
|
|
load(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * & |
156 |
|
|
(load(i,it) - load2(i,it)) |
157 |
|
|
END DO |
158 |
|
|
END DO |
159 |
|
|
|
160 |
|
|
! Deallocate pt_2 no more needed |
161 |
|
|
DEALLOCATE(pt_2) |
162 |
|
|
|
163 |
|
|
END IF ! lonlyone |
164 |
|
|
END IF ! iyr_in .LT. 1850 |
165 |
|
|
|
166 |
|
|
ELSE |
167 |
|
|
WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero |
168 |
|
|
CALL abort_physic('readaerosol','Error : aer_type parameter not accepted',1) |
169 |
|
|
END IF ! type |
170 |
|
|
|
171 |
|
|
|
172 |
|
|
END SUBROUTINE readaerosol |
173 |
|
|
|
174 |
|
|
|
175 |
|
1 |
SUBROUTINE init_aero_fromfile(flag_aerosol) |
176 |
|
|
USE netcdf |
177 |
|
|
USE mod_phys_lmdz_para |
178 |
|
|
USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured |
179 |
|
|
#ifdef CPP_XIOS |
180 |
|
|
USE xios |
181 |
|
|
#endif |
182 |
|
|
IMPLICIT NONE |
183 |
|
|
INTEGER, INTENT(IN) :: flag_aerosol |
184 |
|
|
#ifdef CPP_XIOS |
185 |
|
|
REAL,ALLOCATABLE :: lat_src(:) |
186 |
|
|
REAL,ALLOCATABLE :: lon_src(:) |
187 |
|
|
CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc' |
188 |
|
|
CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc' |
189 |
|
|
INTEGER :: klev_src |
190 |
|
|
INTEGER :: ierr ,ncid, dimID, varid |
191 |
|
|
REAL :: null_array(0) |
192 |
|
|
|
193 |
|
|
IF (flag_aerosol>0 .AND. grid_type==unstructured) THEN |
194 |
|
|
|
195 |
|
|
IF (is_omp_root) THEN |
196 |
|
|
|
197 |
|
|
IF (is_mpi_root) THEN |
198 |
|
|
|
199 |
|
|
IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN |
200 |
|
|
CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) ) |
201 |
|
|
ENDIF |
202 |
|
|
|
203 |
|
|
! Read and test longitudes |
204 |
|
|
CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon") |
205 |
|
|
CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon") |
206 |
|
|
CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" ) |
207 |
|
|
ALLOCATE(lon_src(nbp_lon_src)) |
208 |
|
|
CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" ) |
209 |
|
|
|
210 |
|
|
! Read and test latitudes |
211 |
|
|
CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat") |
212 |
|
|
CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat") |
213 |
|
|
CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" ) |
214 |
|
|
ALLOCATE(lat_src(nbp_lat_src)) |
215 |
|
|
CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" ) |
216 |
|
|
IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN |
217 |
|
|
IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN |
218 |
|
|
CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file') |
219 |
|
|
ENDIF |
220 |
|
|
ENDIF |
221 |
|
|
CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" ) |
222 |
|
|
CALL check_err( nf90_close(ncid),"pb in close" ) |
223 |
|
|
ENDIF |
224 |
|
|
|
225 |
|
|
CALL bcast_mpi(nbp_lat_src) |
226 |
|
|
CALL bcast_mpi(nbp_lon_src) |
227 |
|
|
CALL bcast_mpi(klev_src) |
228 |
|
|
|
229 |
|
|
IF (is_mpi_root ) THEN |
230 |
|
|
CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src) |
231 |
|
|
CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src) |
232 |
|
|
ELSE |
233 |
|
|
CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array ) |
234 |
|
|
CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array) |
235 |
|
|
ENDIF |
236 |
|
|
CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src) |
237 |
|
|
CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.) |
238 |
|
|
|
239 |
|
|
ENDIF |
240 |
|
|
|
241 |
|
|
ENDIF |
242 |
|
|
#endif |
243 |
|
1 |
END SUBROUTINE init_aero_fromfile |
244 |
|
|
|
245 |
|
|
|
246 |
|
|
|
247 |
|
|
|
248 |
|
|
SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out) |
249 |
|
|
!**************************************************************************************** |
250 |
|
|
! Read 12 month aerosol from file and distribute to local process on physical grid. |
251 |
|
|
! Vertical levels, klev_src, may differ from model levels if new file format. |
252 |
|
|
! |
253 |
|
|
! For mpi_root and master thread : |
254 |
|
|
! 1) Open file |
255 |
|
|
! 2) Find vertical dimension klev_src |
256 |
|
|
! 3) Read field month by month |
257 |
|
|
! 4) Close file |
258 |
|
|
! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo) |
259 |
|
|
! - Also the levels and the latitudes have to be inversed |
260 |
|
|
! |
261 |
|
|
! For all processes and threads : |
262 |
|
|
! 6) Scatter global field(klon_glo) to local process domain(klon) |
263 |
|
|
! 7) Test for negative values |
264 |
|
|
!**************************************************************************************** |
265 |
|
|
|
266 |
|
|
USE netcdf |
267 |
|
|
USE dimphy |
268 |
|
|
USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, & |
269 |
|
|
grid2Dto1D_glo, grid_type, unstructured |
270 |
|
|
USE mod_phys_lmdz_para |
271 |
|
|
USE iophy, ONLY : io_lon, io_lat |
272 |
|
|
USE print_control_mod, ONLY: lunout |
273 |
|
|
#ifdef CPP_XIOS |
274 |
|
|
USE xios |
275 |
|
|
#endif |
276 |
|
|
IMPLICIT NONE |
277 |
|
|
|
278 |
|
|
! Input argumets |
279 |
|
|
CHARACTER(len=7), INTENT(IN) :: varname |
280 |
|
|
CHARACTER(len=4), INTENT(IN) :: cyr |
281 |
|
|
CHARACTER(len=8), INTENT(IN) :: filename |
282 |
|
|
|
283 |
|
|
! Output arguments |
284 |
|
|
INTEGER, INTENT(OUT) :: klev_src ! Number of vertical levels in file |
285 |
|
|
REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels |
286 |
|
|
REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels |
287 |
|
|
REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src |
288 |
|
|
REAL, POINTER, DIMENSION(:,:,:) :: pt_year_mpi ! Pointer-variabale from file, 12 month, grid : klon,klev_src |
289 |
|
|
REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months |
290 |
|
|
REAL, DIMENSION(klon_mpi,12) :: psurf_out_mpi ! Surface pression for 12 months |
291 |
|
|
REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column |
292 |
|
|
REAL, DIMENSION(klon_mpi,12) :: load_out_mpi ! Aerosol mass load in each column |
293 |
|
|
INTEGER :: nbr_tsteps ! number of month in file read |
294 |
|
|
|
295 |
|
|
! Local variables |
296 |
|
|
CHARACTER(len=30) :: fname |
297 |
|
|
CHARACTER(len=30) :: cvar |
298 |
|
|
INTEGER :: ncid, dimid, varid |
299 |
|
|
INTEGER :: imth, i, j, k, ierr |
300 |
|
|
REAL :: npole, spole |
301 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varmth |
302 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12 month |
303 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12) |
304 |
|
|
REAL, ALLOCATABLE, DIMENSION(:) :: varktmp |
305 |
|
|
|
306 |
|
|
REAL, ALLOCATABLE :: psurf_glo2D(:,:,:) ! Surface pression for 12 months on dynamics global grid |
307 |
|
|
REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid |
308 |
|
|
REAL, ALLOCATABLE :: load_glo2D(:,:,:) ! Load for 12 months on dynamics global grid |
309 |
|
|
REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid |
310 |
|
|
REAL, ALLOCATABLE, DIMENSION(:,:) :: vartmp |
311 |
|
|
REAL, ALLOCATABLE,DIMENSION(:) :: lon_src ! longitudes in file |
312 |
|
|
REAL, ALLOCATABLE,DIMENSION(:) :: lat_src, lat_src_inv ! latitudes in file |
313 |
|
|
LOGICAL :: new_file ! true if new file format detected |
314 |
|
|
LOGICAL :: invert_lat ! true if the field has to be inverted for latitudes |
315 |
|
|
INTEGER :: nbp_lon, nbp_lat |
316 |
|
|
LOGICAL,SAVE :: first=.TRUE. |
317 |
|
|
!$OMP THREADPRIVATE(first) |
318 |
|
|
|
319 |
|
|
IF (grid_type==unstructured) THEN |
320 |
|
|
nbp_lon=nbp_lon_src |
321 |
|
|
nbp_lat=nbp_lat_src |
322 |
|
|
ELSE |
323 |
|
|
nbp_lon=nbp_lon_ |
324 |
|
|
nbp_lat=nbp_lat_ |
325 |
|
|
ENDIF |
326 |
|
|
|
327 |
|
|
IF (is_mpi_root) THEN |
328 |
|
|
|
329 |
|
|
ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12)) |
330 |
|
|
ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12)) |
331 |
|
|
ALLOCATE(vartmp(nbp_lon,nbp_lat)) |
332 |
|
|
ALLOCATE(lon_src(nbp_lon)) |
333 |
|
|
ALLOCATE(lat_src(nbp_lat)) |
334 |
|
|
ALLOCATE(lat_src_inv(nbp_lat)) |
335 |
|
|
ELSE |
336 |
|
|
ALLOCATE(varyear(0,0,0,0)) |
337 |
|
|
ALLOCATE(psurf_glo2D(0,0,0)) |
338 |
|
|
ALLOCATE(load_glo2D(0,0,0)) |
339 |
|
|
ENDIF |
340 |
|
|
|
341 |
|
|
! Deallocate pointers |
342 |
|
|
IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap) |
343 |
|
|
IF (ASSOCIATED(pt_b)) DEALLOCATE(pt_b) |
344 |
|
|
|
345 |
|
|
IF (is_mpi_root .AND. is_omp_root) THEN |
346 |
|
|
|
347 |
|
|
! 1) Open file |
348 |
|
|
!**************************************************************************************** |
349 |
|
|
! Add suffix to filename |
350 |
|
|
fname = trim(filename)//cyr//'.nc' |
351 |
|
|
|
352 |
|
|
WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname) |
353 |
|
|
CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) ) |
354 |
|
|
|
355 |
|
|
|
356 |
|
|
IF (grid_type/=unstructured) THEN |
357 |
|
|
|
358 |
|
|
! Test for equal longitudes and latitudes in file and model |
359 |
|
|
!**************************************************************************************** |
360 |
|
|
! Read and test longitudes |
361 |
|
|
CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" ) |
362 |
|
|
CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" ) |
363 |
|
|
|
364 |
|
|
IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN |
365 |
|
|
WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname) |
366 |
|
|
WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src |
367 |
|
|
WRITE(lunout,*) 'longitudes in model :', io_lon |
368 |
|
|
|
369 |
|
|
CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1) |
370 |
|
|
END IF |
371 |
|
|
|
372 |
|
|
! Read and test latitudes |
373 |
|
|
CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" ) |
374 |
|
|
CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" ) |
375 |
|
|
|
376 |
|
|
! Invert source latitudes |
377 |
|
|
DO j = 1, nbp_lat |
378 |
|
|
lat_src_inv(j) = lat_src(nbp_lat +1 -j) |
379 |
|
|
END DO |
380 |
|
|
|
381 |
|
|
IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN |
382 |
|
|
! Latitudes are the same |
383 |
|
|
invert_lat=.FALSE. |
384 |
|
|
ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN |
385 |
|
|
! Inverted source latitudes correspond to model latitudes |
386 |
|
|
WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname) |
387 |
|
|
invert_lat=.TRUE. |
388 |
|
|
ELSE |
389 |
|
|
WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname) |
390 |
|
|
WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src |
391 |
|
|
WRITE(lunout,*) 'latitudes in model :', io_lat |
392 |
|
|
CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1) |
393 |
|
|
END IF |
394 |
|
|
ENDIF |
395 |
|
|
|
396 |
|
|
! 2) Check if old or new file is avalabale. |
397 |
|
|
! New type of file should contain the dimension 'lev' |
398 |
|
|
! Old type of file should contain the dimension 'PRESNIVS' |
399 |
|
|
!**************************************************************************************** |
400 |
|
|
ierr = nf90_inq_dimid(ncid, 'lev', dimid) |
401 |
|
|
IF (ierr /= NF90_NOERR) THEN |
402 |
|
|
! Coordinate axe lev not found. Check for presnivs. |
403 |
|
|
ierr = nf90_inq_dimid(ncid, 'presnivs', dimid) |
404 |
|
|
IF (ierr /= NF90_NOERR) THEN |
405 |
|
|
ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid) |
406 |
|
|
IF (ierr /= NF90_NOERR) THEN |
407 |
|
|
! Dimension PRESNIVS not found either |
408 |
|
|
CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file',1) |
409 |
|
|
ELSE |
410 |
|
|
! Old file found |
411 |
|
|
new_file=.FALSE. |
412 |
|
|
WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done' |
413 |
|
|
END IF |
414 |
|
|
ELSE |
415 |
|
|
! New file found |
416 |
|
|
new_file=.TRUE. |
417 |
|
|
WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done' |
418 |
|
|
ENDIF |
419 |
|
|
ELSE |
420 |
|
|
! New file found |
421 |
|
|
new_file=.TRUE. |
422 |
|
|
WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done' |
423 |
|
|
END IF |
424 |
|
|
|
425 |
|
|
! 2) Find vertical dimension klev_src |
426 |
|
|
!**************************************************************************************** |
427 |
|
|
CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" ) |
428 |
|
|
|
429 |
|
|
! Allocate variables depending on the number of vertical levels |
430 |
|
|
ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr) |
431 |
|
|
IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1) |
432 |
|
|
|
433 |
|
|
ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr) |
434 |
|
|
IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1) |
435 |
|
|
|
436 |
|
|
! 3) Read all variables from file |
437 |
|
|
! There is 2 options for the file structure : |
438 |
|
|
! new_file=TRUE : read varyear, ps, pt_ap and pt_b |
439 |
|
|
! new_file=FALSE : read varyear month by month |
440 |
|
|
!**************************************************************************************** |
441 |
|
|
|
442 |
|
|
IF (new_file) THEN |
443 |
|
|
! ++) Check number of month in file opened |
444 |
|
|
!************************************************************************************************** |
445 |
|
|
ierr = nf90_inq_dimid(ncid, 'TIME',dimid) |
446 |
|
|
if (ierr /= NF90_NOERR) THEN |
447 |
|
|
ierr = nf90_inq_dimid(ncid, 'time_counter', dimid) |
448 |
|
|
ENDIF |
449 |
|
|
CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME or time_counter" ) |
450 |
|
|
! IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN |
451 |
|
|
IF (nbr_tsteps /= 12 ) THEN |
452 |
|
|
CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' & |
453 |
|
|
,1) |
454 |
|
|
ENDIF |
455 |
|
|
|
456 |
|
|
! ++) Read the aerosol concentration month by month and concatenate to total variable varyear |
457 |
|
|
!**************************************************************************************** |
458 |
|
|
! Get variable id |
459 |
|
|
!CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) ) |
460 |
|
|
print *,'readaerosol ', TRIM(varname) |
461 |
|
|
IF ( nf90_inq_varid(ncid, TRIM(varname), varid) /= NF90_NOERR ) THEN |
462 |
|
|
! Variable is not there |
463 |
|
|
WRITE(lunout,*) 'Attention '//TRIM(varname)//' is not in aerosol input file' |
464 |
|
|
varyear(:,:,:,:)=0.0 |
465 |
|
|
ELSE |
466 |
|
|
! Get the variable |
467 |
|
|
CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) ) |
468 |
|
|
ENDIF |
469 |
|
|
|
470 |
|
|
! ++) Read surface pression, 12 month in one variable |
471 |
|
|
!**************************************************************************************** |
472 |
|
|
! Get variable id |
473 |
|
|
CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" ) |
474 |
|
|
! Get the variable |
475 |
|
|
CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" ) |
476 |
|
|
|
477 |
|
|
! ++) Read mass load, 12 month in one variable |
478 |
|
|
!**************************************************************************************** |
479 |
|
|
! Get variable id |
480 |
|
|
!CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname)) |
481 |
|
|
IF ( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) /= NF90_NOERR ) THEN |
482 |
|
|
WRITE(lunout,*) 'Attention load_'//TRIM(varname)//' is not in aerosol input file' |
483 |
|
|
load_glo2D(:,:,:)=0.0 |
484 |
|
|
ELSE |
485 |
|
|
! Get the variable |
486 |
|
|
CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) ) |
487 |
|
|
ENDIF |
488 |
|
|
|
489 |
|
|
! ++) Read ap |
490 |
|
|
!**************************************************************************************** |
491 |
|
|
! Get variable id |
492 |
|
|
CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" ) |
493 |
|
|
! Get the variable |
494 |
|
|
CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" ) |
495 |
|
|
|
496 |
|
|
! ++) Read b |
497 |
|
|
!**************************************************************************************** |
498 |
|
|
! Get variable id |
499 |
|
|
CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" ) |
500 |
|
|
! Get the variable |
501 |
|
|
CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" ) |
502 |
|
|
|
503 |
|
|
|
504 |
|
|
|
505 |
|
|
ELSE ! old file |
506 |
|
|
|
507 |
|
|
! ++) Read the aerosol concentration month by month and concatenate to total variable varyear |
508 |
|
|
!**************************************************************************************** |
509 |
|
|
DO imth=1, 12 |
510 |
|
|
IF (imth.EQ.1) THEN |
511 |
|
|
cvar=TRIM(varname)//'JAN' |
512 |
|
|
ELSE IF (imth.EQ.2) THEN |
513 |
|
|
cvar=TRIM(varname)//'FEB' |
514 |
|
|
ELSE IF (imth.EQ.3) THEN |
515 |
|
|
cvar=TRIM(varname)//'MAR' |
516 |
|
|
ELSE IF (imth.EQ.4) THEN |
517 |
|
|
cvar=TRIM(varname)//'APR' |
518 |
|
|
ELSE IF (imth.EQ.5) THEN |
519 |
|
|
cvar=TRIM(varname)//'MAY' |
520 |
|
|
ELSE IF (imth.EQ.6) THEN |
521 |
|
|
cvar=TRIM(varname)//'JUN' |
522 |
|
|
ELSE IF (imth.EQ.7) THEN |
523 |
|
|
cvar=TRIM(varname)//'JUL' |
524 |
|
|
ELSE IF (imth.EQ.8) THEN |
525 |
|
|
cvar=TRIM(varname)//'AUG' |
526 |
|
|
ELSE IF (imth.EQ.9) THEN |
527 |
|
|
cvar=TRIM(varname)//'SEP' |
528 |
|
|
ELSE IF (imth.EQ.10) THEN |
529 |
|
|
cvar=TRIM(varname)//'OCT' |
530 |
|
|
ELSE IF (imth.EQ.11) THEN |
531 |
|
|
cvar=TRIM(varname)//'NOV' |
532 |
|
|
ELSE IF (imth.EQ.12) THEN |
533 |
|
|
cvar=TRIM(varname)//'DEC' |
534 |
|
|
END IF |
535 |
|
|
|
536 |
|
|
! Get variable id |
537 |
|
|
CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) ) |
538 |
|
|
|
539 |
|
|
! Get the variable |
540 |
|
|
CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) ) |
541 |
|
|
|
542 |
|
|
! Store in variable for the whole year |
543 |
|
|
varyear(:,:,:,imth)=varmth(:,:,:) |
544 |
|
|
|
545 |
|
|
END DO |
546 |
|
|
|
547 |
|
|
! Putting dummy |
548 |
|
|
psurf_glo2D(:,:,:) = not_valid |
549 |
|
|
load_glo2D(:,:,:) = not_valid |
550 |
|
|
pt_ap(:) = not_valid |
551 |
|
|
pt_b(:) = not_valid |
552 |
|
|
|
553 |
|
|
END IF |
554 |
|
|
|
555 |
|
|
! 4) Close file |
556 |
|
|
!**************************************************************************************** |
557 |
|
|
CALL check_err( nf90_close(ncid),"pb in close" ) |
558 |
|
|
|
559 |
|
|
|
560 |
|
|
! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo) |
561 |
|
|
!**************************************************************************************** |
562 |
|
|
! Test if vertical levels have to be inversed |
563 |
|
|
|
564 |
|
|
IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN |
565 |
|
|
! WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted' |
566 |
|
|
! WRITE(lunout,*) 'before pt_ap = ', pt_ap |
567 |
|
|
! WRITE(lunout,*) 'before pt_b = ', pt_b |
568 |
|
|
|
569 |
|
|
! Inverse vertical levels for varyear |
570 |
|
|
DO imth=1, 12 |
571 |
|
|
varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly |
572 |
|
|
DO k=1, klev_src |
573 |
|
|
DO j=1, nbp_lat |
574 |
|
|
DO i=1, nbp_lon |
575 |
|
|
varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k) |
576 |
|
|
END DO |
577 |
|
|
END DO |
578 |
|
|
END DO |
579 |
|
|
END DO |
580 |
|
|
|
581 |
|
|
! Inverte vertical axes for pt_ap and pt_b |
582 |
|
|
varktmp(:) = pt_ap(:) |
583 |
|
|
DO k=1, klev_src |
584 |
|
|
pt_ap(k) = varktmp(klev_src+1-k) |
585 |
|
|
END DO |
586 |
|
|
|
587 |
|
|
varktmp(:) = pt_b(:) |
588 |
|
|
DO k=1, klev_src |
589 |
|
|
pt_b(k) = varktmp(klev_src+1-k) |
590 |
|
|
END DO |
591 |
|
|
WRITE(lunout,*) 'after pt_ap = ', pt_ap |
592 |
|
|
WRITE(lunout,*) 'after pt_b = ', pt_b |
593 |
|
|
|
594 |
|
|
ELSE |
595 |
|
|
WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done' |
596 |
|
|
WRITE(lunout,*) 'pt_ap = ', pt_ap |
597 |
|
|
WRITE(lunout,*) 'pt_b = ', pt_b |
598 |
|
|
END IF |
599 |
|
|
|
600 |
|
|
|
601 |
|
|
IF (grid_type/=unstructured) THEN |
602 |
|
|
! - Invert latitudes if necessary |
603 |
|
|
DO imth=1, 12 |
604 |
|
|
IF (invert_lat) THEN |
605 |
|
|
|
606 |
|
|
! Invert latitudes for the variable |
607 |
|
|
varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly |
608 |
|
|
DO k=1,klev_src |
609 |
|
|
DO j=1,nbp_lat |
610 |
|
|
DO i=1,nbp_lon |
611 |
|
|
varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k) |
612 |
|
|
END DO |
613 |
|
|
END DO |
614 |
|
|
END DO |
615 |
|
|
|
616 |
|
|
! Invert latitudes for surface pressure |
617 |
|
|
vartmp(:,:) = psurf_glo2D(:,:,imth) |
618 |
|
|
DO j=1,nbp_lat |
619 |
|
|
DO i=1,nbp_lon |
620 |
|
|
psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j) |
621 |
|
|
END DO |
622 |
|
|
END DO |
623 |
|
|
|
624 |
|
|
! Invert latitudes for the load |
625 |
|
|
vartmp(:,:) = load_glo2D(:,:,imth) |
626 |
|
|
DO j=1,nbp_lat |
627 |
|
|
DO i=1,nbp_lon |
628 |
|
|
load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j) |
629 |
|
|
END DO |
630 |
|
|
END DO |
631 |
|
|
END IF ! invert_lat |
632 |
|
|
|
633 |
|
|
! Do zonal mead at poles and distribut at whole first and last latitude |
634 |
|
|
DO k=1, klev_src |
635 |
|
|
npole=0. ! North pole, j=1 |
636 |
|
|
spole=0. ! South pole, j=nbp_lat |
637 |
|
|
DO i=1,nbp_lon |
638 |
|
|
npole = npole + varyear(i,1,k,imth) |
639 |
|
|
spole = spole + varyear(i,nbp_lat,k,imth) |
640 |
|
|
END DO |
641 |
|
|
npole = npole/REAL(nbp_lon) |
642 |
|
|
spole = spole/REAL(nbp_lon) |
643 |
|
|
varyear(:,1, k,imth) = npole |
644 |
|
|
varyear(:,nbp_lat,k,imth) = spole |
645 |
|
|
END DO |
646 |
|
|
END DO ! imth |
647 |
|
|
|
648 |
|
|
ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr) |
649 |
|
|
IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1) |
650 |
|
|
|
651 |
|
|
! Transform from 2D to 1D field |
652 |
|
|
CALL grid2Dto1D_glo(varyear,varyear_glo1D) |
653 |
|
|
CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D) |
654 |
|
|
CALL grid2Dto1D_glo(load_glo2D,load_glo1D) |
655 |
|
|
|
656 |
|
|
ENDIF |
657 |
|
|
|
658 |
|
|
ELSE |
659 |
|
|
ALLOCATE(varyear_glo1D(0,0,0)) |
660 |
|
|
END IF ! is_mpi_root .AND. is_omp_root |
661 |
|
|
|
662 |
|
|
!$OMP BARRIER |
663 |
|
|
|
664 |
|
|
! 6) Distribute to all processes |
665 |
|
|
! Scatter global field(klon_glo) to local process domain(klon) |
666 |
|
|
! and distribute klev_src to all processes |
667 |
|
|
!**************************************************************************************** |
668 |
|
|
|
669 |
|
|
! Distribute klev_src |
670 |
|
|
CALL bcast(klev_src) |
671 |
|
|
|
672 |
|
|
! Allocate and distribute pt_ap and pt_b |
673 |
|
|
IF (.NOT. ASSOCIATED(pt_ap)) THEN ! if pt_ap is allocated also pt_b is allocated |
674 |
|
|
ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr) |
675 |
|
|
IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1) |
676 |
|
|
END IF |
677 |
|
|
CALL bcast(pt_ap) |
678 |
|
|
CALL bcast(pt_b) |
679 |
|
|
|
680 |
|
|
! Allocate space for output pointer variable at local process |
681 |
|
|
IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year) |
682 |
|
|
ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr) |
683 |
|
|
ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr) |
684 |
|
|
IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1) |
685 |
|
|
|
686 |
|
|
IF (grid_type==unstructured) THEN |
687 |
|
|
#ifdef CPP_XIOS |
688 |
|
|
IF (is_omp_master) THEN |
689 |
|
|
CALL xios_send_field(TRIM(varname)//"_in",varyear) |
690 |
|
|
CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi) |
691 |
|
|
CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D) |
692 |
|
|
CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi) |
693 |
|
|
IF (.not. allocated(psurf_interp)) THEN |
694 |
|
|
! psurf_interp is a shared array |
695 |
|
|
ALLOCATE(psurf_interp(klon_mpi,12)) |
696 |
|
|
CALL xios_send_field("psurf_aerosol_in",psurf_glo2D) |
697 |
|
|
CALL xios_recv_field("psurf_aerosol_out",psurf_interp) |
698 |
|
|
ENDIF |
699 |
|
|
ENDIF |
700 |
|
|
CALL scatter_omp(pt_year_mpi,pt_year) |
701 |
|
|
CALL scatter_omp(load_out_mpi,load_out) |
702 |
|
|
CALL scatter_omp(psurf_interp,psurf_out) |
703 |
|
|
first=.FALSE. |
704 |
|
|
#endif |
705 |
|
|
ELSE |
706 |
|
|
! Scatter global field to local domain at local process |
707 |
|
|
CALL scatter(varyear_glo1D, pt_year) |
708 |
|
|
CALL scatter(psurf_glo1D, psurf_out) |
709 |
|
|
CALL scatter(load_glo1D, load_out) |
710 |
|
|
ENDIF |
711 |
|
|
! 7) Test for negative values |
712 |
|
|
!**************************************************************************************** |
713 |
|
|
IF (MINVAL(pt_year) < 0.) THEN |
714 |
|
|
WRITE(lunout,*) 'Warning! Negative values read from file :', fname |
715 |
|
|
END IF |
716 |
|
|
|
717 |
|
|
END SUBROUTINE get_aero_fromfile |
718 |
|
|
|
719 |
|
|
|
720 |
|
|
SUBROUTINE check_err(status,text) |
721 |
|
|
USE netcdf |
722 |
|
|
USE print_control_mod, ONLY: lunout |
723 |
|
|
IMPLICIT NONE |
724 |
|
|
|
725 |
|
|
INTEGER, INTENT (IN) :: status |
726 |
|
|
CHARACTER(len=*), INTENT (IN), OPTIONAL :: text |
727 |
|
|
|
728 |
|
|
IF (status /= NF90_NOERR) THEN |
729 |
|
|
WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status |
730 |
|
|
IF (PRESENT(text)) THEN |
731 |
|
|
WRITE(lunout,*) 'Error in get_aero_fromfile : ',text |
732 |
|
|
END IF |
733 |
|
|
CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1) |
734 |
|
|
END IF |
735 |
|
|
|
736 |
|
|
END SUBROUTINE check_err |
737 |
|
|
|
738 |
|
|
|
739 |
|
|
END MODULE readaerosol_mod |