Context Navigation

readaerosol_mod.F90 @ 5451

Last change on this file since 5451 was 4727, checked in by idelkadi, 15 months ago
Merged trunk changes -r4488:4726 LMDZ_ECRad branch
File size: 30.5 KB

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

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

Download in other formats:

Original Format