Context Navigation

readaerosol.F90 @ 2408

Last change on this file since 2408 was 2408, checked in by Laurent Fairhead, 9 years ago
Merged trunk changes r2298:2396 into testing branch
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: 24.7 KB

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

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

Download in other formats:

Original Format