Context Navigation

readaerosol_mod.F90 @ 4894

Last change on this file since 4894 was 4627, checked in by acozic, 16 months ago
add a test on aerosol_couple before reading aerosol files in init_aero_fromfile. If we are using inca aerosols, we don't need to read these files
File size: 30.5 KB

Rev	Line
[3440]	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
[4627]	175	SUBROUTINE init_aero_fromfile(flag_aerosol, aerosol_couple)
[3440]	176	USE netcdf
	177	USE mod_phys_lmdz_para
	178	USE mod_grid_phy_lmdz, ONLY: grid_type, unstructured
[4619]	179	USE lmdz_xios
[3440]	180	IMPLICIT NONE
[4627]	181
[3440]	182	INTEGER, INTENT(IN) :: flag_aerosol
[4627]	183	LOGICAL, INTENT(IN) :: aerosol_couple
	184
[3440]	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
[4619]	193	IF (using_xios) THEN
[4627]	194	IF (flag_aerosol>0 .AND. grid_type==unstructured .AND. (.NOT. aerosol_couple) ) THEN
[3440]	195
[4619]	196	IF (is_omp_root) THEN
[3440]	197
[4619]	198	IF (is_mpi_root) THEN
[3440]	199
[4619]	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
[3440]	203
[4619]	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" )
[3440]	210
[4619]	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
[3440]	221	ENDIF
[4619]	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" )
[3440]	224	ENDIF
	225
[4619]	226	CALL bcast_mpi(nbp_lat_src)
	227	CALL bcast_mpi(nbp_lon_src)
	228	CALL bcast_mpi(klev_src)
[3440]	229
[4619]	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
[3440]	240	ENDIF
	241
[4619]	242	ENDIF
	243	ENDIF !using_xios
	244	END SUBROUTINE init_aero_fromfile
[3440]	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
[4619]	274	USE lmdz_xios
[3440]	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)
[4103]	691	IF (.not. allocated(psurf_interp)) THEN
	692	! psurf_interp is a shared array
[3440]	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