Context Navigation

readaerosol.F90 @ 5177

Last change on this file since 5177 was 3413, checked in by Laurent Fairhead, 6 years ago
Inclusion of Yann's latest (summer/fall 2018) modifications for convergence of DYNAMICO/LMDZ physics YM/LF
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: 30.3 KB

Line
1	! $Id: readaerosol.F90 3413 2018-11-12 12:52:29Z abarral $
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	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	USE xios
181	IMPLICIT NONE
182	INTEGER, INTENT(IN) :: flag_aerosol
183	REAL,ALLOCATABLE :: lat_src(:)
184	REAL,ALLOCATABLE :: lon_src(:)
185	CHARACTER(LEN=*),PARAMETER :: file_aerosol='aerosols.nat.nc'
186	CHARACTER(LEN=*),PARAMETER :: file_so4='so4.nat.nc'
187	INTEGER :: klev_src
188	INTEGER :: ierr ,ncid, dimID, varid
189	REAL :: null_array(0)
190
191	IF (flag_aerosol>0 .AND. grid_type==unstructured) THEN
192
193	IF (is_omp_root) THEN
194
195	IF (is_mpi_root) THEN
196
197	IF (nf90_open(TRIM(file_aerosol), NF90_NOWRITE, ncid) /= NF90_NOERR) THEN
198	CALL check_err( nf90_open(TRIM(file_so4), NF90_NOWRITE, ncid), "pb open "//trim(file_so4) )
199	ENDIF
200
201	! Read and test longitudes
202	CALL check_err( nf90_inq_dimid(ncid, "lon", dimID),"pb inq dim lon")
203	CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lon_src),"pb inq dim lon")
204	CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
205	ALLOCATE(lon_src(nbp_lon_src))
206	CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
207
208	! Read and test latitudes
209	CALL check_err( nf90_inq_dimid(ncid, "lat", dimID),"pb inq dim lat")
210	CALL check_err( nf90_inquire_dimension(ncid, dimID, len = nbp_lat_src),"pb inq dim lat")
211	CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
212	ALLOCATE(lat_src(nbp_lat_src))
213	CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
214	IF (nf90_inq_dimid(ncid, 'lev', dimid) /= NF90_NOERR) THEN
215	IF (nf90_inq_dimid(ncid, 'presnivs', dimid)/= NF90_NOERR) THEN
216	CALL check_err(nf90_inq_dimid(ncid, 'PRESNIVS', dimid),'dimension lev,PRESNIVS or presnivs not in file')
217	ENDIF
218	ENDIF
219	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
220	CALL check_err( nf90_close(ncid),"pb in close" )
221	ENDIF
222
223	CALL bcast_mpi(nbp_lat_src)
224	CALL bcast_mpi(nbp_lon_src)
225	CALL bcast_mpi(klev_src)
226
227	IF (is_mpi_root ) THEN
228	CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=nbp_lat_src, jbegin=0, latvalue_1d=lat_src)
229	CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=nbp_lon_src, ibegin=0, lonvalue_1d=lon_src)
230	ELSE
231	CALL xios_set_domain_attr("domain_aerosol",nj_glo=nbp_lat_src, nj=0, jbegin=0, latvalue_1d=null_array )
232	CALL xios_set_domain_attr("domain_aerosol",ni_glo=nbp_lon_src, ni=0, ibegin=0, lonvalue_1d=null_array)
233	ENDIF
234	CALL xios_set_axis_attr("axis_aerosol",n_glo=klev_src)
235	CALL xios_set_fieldgroup_attr("aerosols", enabled=.TRUE.)
236
237	ENDIF
238
239	ENDIF
240
241	END SUBROUTINE init_aero_fromfile
242
243
244
245
246	SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
247	!****************************************************************************************
248	! Read 12 month aerosol from file and distribute to local process on physical grid.
249	! Vertical levels, klev_src, may differ from model levels if new file format.
250	!
251	! For mpi_root and master thread :
252	! 1) Open file
253	! 2) Find vertical dimension klev_src
254	! 3) Read field month by month
255	! 4) Close file
256	! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
257	! - Also the levels and the latitudes have to be inversed
258	!
259	! For all processes and threads :
260	! 6) Scatter global field(klon_glo) to local process domain(klon)
261	! 7) Test for negative values
262	!****************************************************************************************
263
264	USE netcdf
265	USE dimphy
266	USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
267	grid2Dto1D_glo, grid_type, unstructured
268	USE mod_phys_lmdz_para
269	USE iophy, ONLY : io_lon, io_lat
270	USE print_control_mod, ONLY: lunout
271	USE xios
272
273	IMPLICIT NONE
274
275	! Input argumets
276	CHARACTER(len=7), INTENT(IN) :: varname
277	CHARACTER(len=4), INTENT(IN) :: cyr
278	CHARACTER(len=8), INTENT(IN) :: filename
279
280	! Output arguments
281	INTEGER, INTENT(OUT) :: klev_src ! Number of vertical levels in file
282	REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels
283	REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels
284	REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src
285	REAL, POINTER, DIMENSION(:,:,:) :: pt_year_mpi ! Pointer-variabale from file, 12 month, grid : klon,klev_src
286	REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months
287	REAL, DIMENSION(klon_mpi,12) :: psurf_out_mpi ! Surface pression for 12 months
288	REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column
289	REAL, DIMENSION(klon_mpi,12) :: load_out_mpi ! Aerosol mass load in each column
290	INTEGER :: nbr_tsteps ! number of month in file read
291
292	! Local variables
293	CHARACTER(len=30) :: fname
294	CHARACTER(len=30) :: cvar
295	INTEGER :: ncid, dimid, varid
296	INTEGER :: imth, i, j, k, ierr
297	REAL :: npole, spole
298	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varmth
299	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12 month
300	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12)
301	REAL, ALLOCATABLE, DIMENSION(:) :: varktmp
302
303	REAL, ALLOCATABLE :: psurf_glo2D(:,:,:) ! Surface pression for 12 months on dynamics global grid
304	REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid
305	REAL, ALLOCATABLE :: load_glo2D(:,:,:) ! Load for 12 months on dynamics global grid
306	REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid
307	REAL, ALLOCATABLE, DIMENSION(:,:) :: vartmp
308	REAL, ALLOCATABLE,DIMENSION(:) :: lon_src ! longitudes in file
309	REAL, ALLOCATABLE,DIMENSION(:) :: lat_src, lat_src_inv ! latitudes in file
310	LOGICAL :: new_file ! true if new file format detected
311	LOGICAL :: invert_lat ! true if the field has to be inverted for latitudes
312	INTEGER :: nbp_lon, nbp_lat
313	LOGICAL,SAVE :: first=.TRUE.
314	!$OMP THREADPRIVATE(first)
315
316	IF (grid_type==unstructured) THEN
317	nbp_lon=nbp_lon_src
318	nbp_lat=nbp_lat_src
319	ELSE
320	nbp_lon=nbp_lon_
321	nbp_lat=nbp_lat_
322	ENDIF
323
324	IF (is_mpi_root) THEN
325
326	ALLOCATE(psurf_glo2D(nbp_lon,nbp_lat,12))
327	ALLOCATE(load_glo2D(nbp_lon,nbp_lat,12))
328	ALLOCATE(vartmp(nbp_lon,nbp_lat))
329	ALLOCATE(lon_src(nbp_lon))
330	ALLOCATE(lat_src(nbp_lat))
331	ALLOCATE(lat_src_inv(nbp_lat))
332	ELSE
333	ALLOCATE(varyear(0,0,0,0))
334	ALLOCATE(psurf_glo2D(0,0,0))
335	ALLOCATE(load_glo2D(0,0,0))
336	ENDIF
337
338	! Deallocate pointers
339	IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
340	IF (ASSOCIATED(pt_b)) DEALLOCATE(pt_b)
341
342	IF (is_mpi_root .AND. is_omp_root) THEN
343
344	! 1) Open file
345	!****************************************************************************************
346	! Add suffix to filename
347	fname = trim(filename)//cyr//'.nc'
348
349	WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
350	CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
351
352
353	IF (grid_type/=unstructured) THEN
354
355	! Test for equal longitudes and latitudes in file and model
356	!****************************************************************************************
357	! Read and test longitudes
358	CALL check_err( nf90_inq_varid(ncid, 'lon', varid),"pb inq lon" )
359	CALL check_err( nf90_get_var(ncid, varid, lon_src(:)),"pb get lon" )
360
361	IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
362	WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
363	WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
364	WRITE(lunout,*) 'longitudes in model :', io_lon
365
366	CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
367	END IF
368
369	! Read and test latitudes
370	CALL check_err( nf90_inq_varid(ncid, 'lat', varid),"pb inq lat" )
371	CALL check_err( nf90_get_var(ncid, varid, lat_src(:)),"pb get lat" )
372
373	! Invert source latitudes
374	DO j = 1, nbp_lat
375	lat_src_inv(j) = lat_src(nbp_lat +1 -j)
376	END DO
377
378	IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
379	! Latitudes are the same
380	invert_lat=.FALSE.
381	ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
382	! Inverted source latitudes correspond to model latitudes
383	WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
384	invert_lat=.TRUE.
385	ELSE
386	WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
387	WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src
388	WRITE(lunout,*) 'latitudes in model :', io_lat
389	CALL abort_physic('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
390	END IF
391	ENDIF
392
393	! 2) Check if old or new file is avalabale.
394	! New type of file should contain the dimension 'lev'
395	! Old type of file should contain the dimension 'PRESNIVS'
396	!****************************************************************************************
397	ierr = nf90_inq_dimid(ncid, 'lev', dimid)
398	IF (ierr /= NF90_NOERR) THEN
399	! Coordinate axe lev not found. Check for presnivs.
400	ierr = nf90_inq_dimid(ncid, 'presnivs', dimid)
401	IF (ierr /= NF90_NOERR) THEN
402	ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
403	IF (ierr /= NF90_NOERR) THEN
404	! Dimension PRESNIVS not found either
405	CALL abort_physic('get_aero_fromfile', 'dimension lev,PRESNIVS or presnivs not in file',1)
406	ELSE
407	! Old file found
408	new_file=.FALSE.
409	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
410	END IF
411	ELSE
412	! New file found
413	new_file=.TRUE.
414	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
415	ENDIF
416	ELSE
417	! New file found
418	new_file=.TRUE.
419	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
420	END IF
421
422	! 2) Find vertical dimension klev_src
423	!****************************************************************************************
424	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
425
426	! Allocate variables depending on the number of vertical levels
427	ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr)
428	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1)
429
430	ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
431	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1)
432
433	! 3) Read all variables from file
434	! There is 2 options for the file structure :
435	! new_file=TRUE : read varyear, ps, pt_ap and pt_b
436	! new_file=FALSE : read varyear month by month
437	!****************************************************************************************
438
439	IF (new_file) THEN
440	! ++) Check number of month in file opened
441	!**************************************************************************************************
442	ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
443	if (ierr /= NF90_NOERR) THEN
444	ierr = nf90_inq_dimid(ncid, 'time_counter', dimid)
445	ENDIF
446	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps),"pb inq dim TIME or time_counter" )
447	! IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
448	IF (nbr_tsteps /= 12 ) THEN
449	CALL abort_physic('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)' &
450	,1)
451	ENDIF
452
453	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
454	!****************************************************************************************
455	! Get variable id
456	!CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid),"pb inq var "//TRIM(varname) )
457	print *,'readaerosol ', TRIM(varname)
458	IF ( nf90_inq_varid(ncid, TRIM(varname), varid) /= NF90_NOERR ) THEN
459	! Variable is not there
460	WRITE(lunout,*) 'Attention '//TRIM(varname)//' is not in aerosol input file'
461	varyear(:,:,:,:)=0.0
462	ELSE
463	! Get the variable
464	CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)),"pb get var "//TRIM(varname) )
465	ENDIF
466
467	! ++) Read surface pression, 12 month in one variable
468	!****************************************************************************************
469	! Get variable id
470	CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
471	! Get the variable
472	CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
473
474	! ++) Read mass load, 12 month in one variable
475	!****************************************************************************************
476	! Get variable id
477	!CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) ,"pb inq var load_"//TRIM(varname))
478	IF ( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) /= NF90_NOERR ) THEN
479	WRITE(lunout,*) 'Attention load_'//TRIM(varname)//' is not in aerosol input file'
480	load_glo2D(:,:,:)=0.0
481	ELSE
482	! Get the variable
483	CALL check_err( nf90_get_var(ncid, varid, load_glo2D),"pb get var load_"//TRIM(varname) )
484	ENDIF
485
486	! ++) Read ap
487	!****************************************************************************************
488	! Get variable id
489	CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
490	! Get the variable
491	CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
492
493	! ++) Read b
494	!****************************************************************************************
495	! Get variable id
496	CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
497	! Get the variable
498	CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
499
500
501
502	ELSE ! old file
503
504	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
505	!****************************************************************************************
506	DO imth=1, 12
507	IF (imth.EQ.1) THEN
508	cvar=TRIM(varname)//'JAN'
509	ELSE IF (imth.EQ.2) THEN
510	cvar=TRIM(varname)//'FEB'
511	ELSE IF (imth.EQ.3) THEN
512	cvar=TRIM(varname)//'MAR'
513	ELSE IF (imth.EQ.4) THEN
514	cvar=TRIM(varname)//'APR'
515	ELSE IF (imth.EQ.5) THEN
516	cvar=TRIM(varname)//'MAY'
517	ELSE IF (imth.EQ.6) THEN
518	cvar=TRIM(varname)//'JUN'
519	ELSE IF (imth.EQ.7) THEN
520	cvar=TRIM(varname)//'JUL'
521	ELSE IF (imth.EQ.8) THEN
522	cvar=TRIM(varname)//'AUG'
523	ELSE IF (imth.EQ.9) THEN
524	cvar=TRIM(varname)//'SEP'
525	ELSE IF (imth.EQ.10) THEN
526	cvar=TRIM(varname)//'OCT'
527	ELSE IF (imth.EQ.11) THEN
528	cvar=TRIM(varname)//'NOV'
529	ELSE IF (imth.EQ.12) THEN
530	cvar=TRIM(varname)//'DEC'
531	END IF
532
533	! Get variable id
534	CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
535
536	! Get the variable
537	CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
538
539	! Store in variable for the whole year
540	varyear(:,:,:,imth)=varmth(:,:,:)
541
542	END DO
543
544	! Putting dummy
545	psurf_glo2D(:,:,:) = not_valid
546	load_glo2D(:,:,:) = not_valid
547	pt_ap(:) = not_valid
548	pt_b(:) = not_valid
549
550	END IF
551
552	! 4) Close file
553	!****************************************************************************************
554	CALL check_err( nf90_close(ncid),"pb in close" )
555
556
557	! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
558	!****************************************************************************************
559	! Test if vertical levels have to be inversed
560
561	IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
562	! WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' needs to be inverted'
563	! WRITE(lunout,*) 'before pt_ap = ', pt_ap
564	! WRITE(lunout,*) 'before pt_b = ', pt_b
565
566	! Inverse vertical levels for varyear
567	DO imth=1, 12
568	varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
569	DO k=1, klev_src
570	DO j=1, nbp_lat
571	DO i=1, nbp_lon
572	varyear(i,j,k,imth) = varmth(i,j,klev_src+1-k)
573	END DO
574	END DO
575	END DO
576	END DO
577
578	! Inverte vertical axes for pt_ap and pt_b
579	varktmp(:) = pt_ap(:)
580	DO k=1, klev_src
581	pt_ap(k) = varktmp(klev_src+1-k)
582	END DO
583
584	varktmp(:) = pt_b(:)
585	DO k=1, klev_src
586	pt_b(k) = varktmp(klev_src+1-k)
587	END DO
588	WRITE(lunout,*) 'after pt_ap = ', pt_ap
589	WRITE(lunout,*) 'after pt_b = ', pt_b
590
591	ELSE
592	WRITE(lunout,*) 'Vertical axis in file ',TRIM(fname), ' is ok, no vertical inversion is done'
593	WRITE(lunout,*) 'pt_ap = ', pt_ap
594	WRITE(lunout,*) 'pt_b = ', pt_b
595	END IF
596
597
598	IF (grid_type/=unstructured) THEN
599	! - Invert latitudes if necessary
600	DO imth=1, 12
601	IF (invert_lat) THEN
602
603	! Invert latitudes for the variable
604	varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
605	DO k=1,klev_src
606	DO j=1,nbp_lat
607	DO i=1,nbp_lon
608	varyear(i,j,k,imth) = varmth(i,nbp_lat+1-j,k)
609	END DO
610	END DO
611	END DO
612
613	! Invert latitudes for surface pressure
614	vartmp(:,:) = psurf_glo2D(:,:,imth)
615	DO j=1,nbp_lat
616	DO i=1,nbp_lon
617	psurf_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
618	END DO
619	END DO
620
621	! Invert latitudes for the load
622	vartmp(:,:) = load_glo2D(:,:,imth)
623	DO j=1,nbp_lat
624	DO i=1,nbp_lon
625	load_glo2D(i,j,imth)= vartmp(i,nbp_lat+1-j)
626	END DO
627	END DO
628	END IF ! invert_lat
629
630	! Do zonal mead at poles and distribut at whole first and last latitude
631	DO k=1, klev_src
632	npole=0. ! North pole, j=1
633	spole=0. ! South pole, j=nbp_lat
634	DO i=1,nbp_lon
635	npole = npole + varyear(i,1,k,imth)
636	spole = spole + varyear(i,nbp_lat,k,imth)
637	END DO
638	npole = npole/REAL(nbp_lon)
639	spole = spole/REAL(nbp_lon)
640	varyear(:,1, k,imth) = npole
641	varyear(:,nbp_lat,k,imth) = spole
642	END DO
643	END DO ! imth
644
645	ALLOCATE(varyear_glo1D(klon_glo, klev_src, 12), stat=ierr)
646	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 3',1)
647
648	! Transform from 2D to 1D field
649	CALL grid2Dto1D_glo(varyear,varyear_glo1D)
650	CALL grid2Dto1D_glo(psurf_glo2D,psurf_glo1D)
651	CALL grid2Dto1D_glo(load_glo2D,load_glo1D)
652
653	ENDIF
654
655	ELSE
656	ALLOCATE(varyear_glo1D(0,0,0))
657	END IF ! is_mpi_root .AND. is_omp_root
658
659	!$OMP BARRIER
660
661	! 6) Distribute to all processes
662	! Scatter global field(klon_glo) to local process domain(klon)
663	! and distribute klev_src to all processes
664	!****************************************************************************************
665
666	! Distribute klev_src
667	CALL bcast(klev_src)
668
669	! Allocate and distribute pt_ap and pt_b
670	IF (.NOT. ASSOCIATED(pt_ap)) THEN ! if pt_ap is allocated also pt_b is allocated
671	ALLOCATE(pt_ap(klev_src), pt_b(klev_src), stat=ierr)
672	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1)
673	END IF
674	CALL bcast(pt_ap)
675	CALL bcast(pt_b)
676
677	! Allocate space for output pointer variable at local process
678	IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
679	ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
680	ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
681	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
682
683	IF (grid_type==unstructured) THEN
684	IF (is_omp_master) THEN
685	CALL xios_send_field(TRIM(varname)//"_in",varyear)
686	CALL xios_recv_field(TRIM(varname)//"_out",pt_year_mpi)
687	CALL xios_send_field("load_"//TRIM(varname)//"_in",load_glo2D)
688	CALL xios_recv_field("load_"//TRIM(varname)//"_out",load_out_mpi)
689	IF (first) THEN
690	ALLOCATE(psurf_interp(klon_mpi,12))
691	CALL xios_send_field("psurf_aerosol_in",psurf_glo2D)
692	CALL xios_recv_field("psurf_aerosol_out",psurf_interp)
693	ENDIF
694	ENDIF
695	CALL scatter_omp(pt_year_mpi,pt_year)
696	CALL scatter_omp(load_out_mpi,load_out)
697	CALL scatter_omp(psurf_interp,psurf_out)
698	first=.FALSE.
699	ELSE
700	! Scatter global field to local domain at local process
701	CALL scatter(varyear_glo1D, pt_year)
702	CALL scatter(psurf_glo1D, psurf_out)
703	CALL scatter(load_glo1D, load_out)
704	ENDIF
705	! 7) Test for negative values
706	!****************************************************************************************
707	IF (MINVAL(pt_year) < 0.) THEN
708	WRITE(lunout,*) 'Warning! Negative values read from file :', fname
709	END IF
710
711	END SUBROUTINE get_aero_fromfile
712
713
714	SUBROUTINE check_err(status,text)
715	USE netcdf
716	USE print_control_mod, ONLY: lunout
717	IMPLICIT NONE
718
719	INTEGER, INTENT (IN) :: status
720	CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
721
722	IF (status /= NF90_NOERR) THEN
723	WRITE(lunout,*) 'Error in get_aero_fromfile, netcdf error code = ',status
724	IF (PRESENT(text)) THEN
725	WRITE(lunout,*) 'Error in get_aero_fromfile : ',text
726	END IF
727	CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1)
728	END IF
729
730	END SUBROUTINE check_err
731
732
733	END MODULE readaerosol_mod

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

Download in other formats:

Original Format