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