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