Context Navigation

readaerosol.F90 @ 5428

Last change on this file since 5428 was 2841, checked in by acozic, 8 years ago
remove the reading of p0 in aerosol files because it's not use and add the possibility to have new axes in these files (presnivs and time_counter)
Property copyright set to Name of program: LMDZ Creation date: 1984 Version: LMDZ5 License: CeCILL version 2 Holder: Laboratoire de m\'et\'eorologie dynamique, CNRS, UMR 8539 See the license file in the root directory Property svn:eol-style set to `native` Property svn:keywords set to `Author Date Id Revision`
File size: 25.1 KB

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

Note: See TracBrowser for help on using the repository browser.

Download in other formats:

Original Format