Context Navigation

readaerosol_mod.F90 @ 5442

Last change on this file since 5442 was 5160, checked in by abarral, 5 months ago
Put .h into modules
File size: 30.3 KB

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

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

Download in other formats:

Original Format