Context Navigation

readaerosol.F90 @ 3435

Last change on this file since 3435 was 3435, checked in by Laurent Fairhead, 5 years ago

"Historic" :-) commit merging the physics branch used for DYNAMICO with the LMDZ trunk.
The same physics branch can now be used seamlessly with the traditional lon-lat LMDZ
dynamical core and DYNAMICO.
Testing consisted in running a lon-lat LMDZ bucket simulation with the NPv6.1 physics package
with the original trunk sources and the merged sources. Tests were succesful in the sense that
numeric continuity was preserved in the restart files from both simulation. Further tests
included running both versions of the physics codes for one year in a LMDZOR setting in which
the restart files also came out identical.

Caution:

as the physics package now manages unstructured grids, grid information needs to be transmitted

to the surface scheme ORCHIDEE. This means that the interface defined in surf_land_orchidee_mod.F90
is only compatible with ORCHIDEE version orchidee2.1 and later versions. If previous versions of
ORCHIDEE need to be used, the CPP key ORCHIDEE_NOUNSTRUCT needs to be set at compilation time.
This is done automatically if makelmdz/makelmdz_fcm are called with the veget orchidee2.0 switch

due to a limitation in XIOS, the time at which limit conditions will be read in by DYNAMICO will be

delayed by one physic timestep with respect to the time it is read in by the lon-lat model. This is caused
by the line

IF (MOD(itime-1, lmt_pas) == 0 .OR. (jour_lu /= jour .AND. grid_type /= unstructured)) THEN ! time to read

in limit_read_mod.F90

Work still needed on COSP integration and XML files for DYNAMICO

EM, 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 3435 2019-01-22 15:21:59Z fairhead $
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