Context Navigation

readaerosol_mod.F90 @ 5080

Last change on this file since 5080 was 5075, checked in by abarral, 4 months ago

[continued & end] replace netcdf by lmdz_netcdf.F90 wrapper
"use netcdf" is now only used in lmdz_netcdf.F90 (except ecrad and obsolete/)
<include "netcdf.inc"> is now likewise only used in lmdz_netcdf.F90.

systematically specify explicitely <USE lmdz_netcdf, ONLY:> (probably left some missing, to correct later on)

Further replacement of nf_put_* by nf90_put_* (same for _get_)

[minor] replace deprecated boolean operators along the way

File size: 30.6 KB

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