Context Navigation

readaerosol.F90 @ 3414

Last change on this file since 3414 was 3413, checked in by Laurent Fairhead, 6 years ago
Inclusion of Yann's latest (summer/fall 2018) modifications for convergence of DYNAMICO/LMDZ physics 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

Rev	Line
[1179]	1	! $Id: readaerosol.F90 3413 2018-11-12 12:52:29Z fairhead $
[524]	2	!
[1179]	3	MODULE readaerosol_mod
	4
	5	REAL, SAVE :: not_valid=-333.
[3413]	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)
[1179]	13
	14	CONTAINS
	15
[1492]	16	SUBROUTINE readaerosol(name_aero, type, filename, iyr_in, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1179]	17
	18	!****************************************************************************************
	19	! This routine will read the aersosol from file.
[524]	20	!
[1179]	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
[1150]	25	!
[1179]	26	! Return pointer, pt_out, to the year read or result from interpolation
	27	!****************************************************************************************
	28	USE dimphy
[2311]	29	USE print_control_mod, ONLY: lunout
[524]	30
[1179]	31	IMPLICIT NONE
[766]	32
[1179]	33	! Input arguments
	34	CHARACTER(len=7), INTENT(IN) :: name_aero
[1492]	35	CHARACTER(len=*), INTENT(IN) :: type ! actuel, annuel, scenario or preind
	36	CHARACTER(len=8), INTENT(IN) :: filename
[1179]	37	INTEGER, INTENT(IN) :: iyr_in
[1150]	38
[1179]	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
[1270]	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
[1150]	46
[1179]	47	! Local variables
[1270]	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.
[1150]	54
[1179]	55	!****************************************************************************************
	56	! Read data depending on aer_type
	57	!
	58	!****************************************************************************************
[524]	59
[1179]	60	IF (type == 'actuel') THEN
	61	! Read and return data for year 1980
	62	!****************************************************************************************
	63	cyear='1980'
[1270]	64	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
	65	! pt_out has dimensions (klon, klev_src, 12)
[2841]	66	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1179]	67
[640]	68
[1179]	69	ELSE IF (type == 'preind') THEN
	70	! Read and return data from file with suffix .nat
	71	!****************************************************************************************
	72	cyear='.nat'
[1270]	73	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
	74	! pt_out has dimensions (klon, klev_src, 12)
[2841]	75	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1179]	76
[1321]	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)
[2841]	84	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1321]	85
[1179]	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
[1270]	92	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
	93	! pt_out has dimensions (klon, klev_src, 12)
[2841]	94	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1179]	95
	96	ELSE IF (iyr_in .GE. 2100) THEN
	97	cyear='2100'
	98	WRITE(lunout,*) 'get_aero 2 iyr_in=', iyr_in,' ',cyear
[1270]	99	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
	100	! pt_out has dimensions (klon, klev_src, 12)
[2841]	101	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1179]	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
[1270]	119	! get_aero_fromfile returns pt_out allocated and initialized with data for 12 month
	120	! pt_out has dimensions (klon, klev_src, 12)
[2841]	121	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src, pt_ap, pt_b, pt_out, psurf, load)
[1179]	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)
[1270]	131	! get_aero_fromfile returns pt_2 allocated and initialized with data for 12 month
	132	! pt_2 has dimensions (klon, klev_src, 12)
[2841]	133	CALL get_aero_fromfile(name_aero, cyear, filename, klev_src2, pt_ap, pt_b, pt_2, psurf2, load2)
[1179]	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'
[2311]	137	CALL abort_physic('readaersosol','Error in number of vertical levels',1)
[1179]	138	END IF
	139
	140	! Linare interpolate to the actual year:
[1270]	141	DO it=1,12
[1179]	142	DO k=1,klev_src
	143	DO i = 1, klon
	144	pt_out(i,k,it) = &
[1403]	145	pt_out(i,k,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
[1179]	146	(pt_out(i,k,it) - pt_2(i,k,it))
	147	END DO
	148	END DO
[1143]	149
[1179]	150	DO i = 1, klon
[1270]	151	psurf(i,it) = &
[1403]	152	psurf(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
[1270]	153	(psurf(i,it) - psurf2(i,it))
[640]	154
[1270]	155	load(i,it) = &
[1403]	156	load(i,it) - REAL(iyr_in-iyr1)/REAL(iyr2-iyr1) * &
[1270]	157	(load(i,it) - load2(i,it))
[1179]	158	END DO
	159	END DO
[524]	160
[1179]	161	! Deallocate pt_2 no more needed
	162	DEALLOCATE(pt_2)
	163
	164	END IF ! lonlyone
	165	END IF ! iyr_in .LT. 1850
[524]	166
[1179]	167	ELSE
[1492]	168	WRITE(lunout,*)'This option is not implemented : aer_type = ', type,' name_aero=',name_aero
[2311]	169	CALL abort_physic('readaerosol','Error : aer_type parameter not accepted',1)
[1179]	170	END IF ! type
[640]	171
[1151]	172
[1179]	173	END SUBROUTINE readaerosol
[1143]	174
	175
[3413]	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
[2841]	246	SUBROUTINE get_aero_fromfile(varname, cyr, filename, klev_src, pt_ap, pt_b, pt_year, psurf_out, load_out)
[1179]	247	!****************************************************************************************
[1270]	248	! Read 12 month aerosol from file and distribute to local process on physical grid.
[1179]	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
[2346]	256	! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
[1179]	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	!****************************************************************************************
[1150]	263
[1179]	264	USE netcdf
	265	USE dimphy
[3413]	266	USE mod_grid_phy_lmdz, ONLY: nbp_lon_=>nbp_lon, nbp_lat_=>nbp_lat, klon_glo, &
	267	grid2Dto1D_glo, grid_type, unstructured
[1179]	268	USE mod_phys_lmdz_para
[1183]	269	USE iophy, ONLY : io_lon, io_lat
[2311]	270	USE print_control_mod, ONLY: lunout
[3413]	271	USE xios
[1179]	272
	273	IMPLICIT NONE
[524]	274
[1179]	275	! Input argumets
	276	CHARACTER(len=7), INTENT(IN) :: varname
	277	CHARACTER(len=4), INTENT(IN) :: cyr
[1492]	278	CHARACTER(len=8), INTENT(IN) :: filename
[1150]	279
[1179]	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
[1270]	284	REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src
[3413]	285	REAL, POINTER, DIMENSION(:,:,:) :: pt_year_mpi ! Pointer-variabale from file, 12 month, grid : klon,klev_src
[1270]	286	REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months
[3413]	287	REAL, DIMENSION(klon_mpi,12) :: psurf_out_mpi ! Surface pression for 12 months
[1270]	288	REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column
[3413]	289	REAL, DIMENSION(klon_mpi,12) :: load_out_mpi ! Aerosol mass load in each column
[1319]	290	INTEGER :: nbr_tsteps ! number of month in file read
[1150]	291
[1179]	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
[1270]	299	REAL, ALLOCATABLE, DIMENSION(:,:,:,:) :: varyear ! Global variable read from file, 12 month
	300	REAL, ALLOCATABLE, DIMENSION(:,:,:) :: varyear_glo1D !(klon_glo, klev_src, 12)
[1179]	301	REAL, ALLOCATABLE, DIMENSION(:) :: varktmp
[524]	302
[3413]	303	REAL, ALLOCATABLE :: psurf_glo2D(:,:,:) ! Surface pression for 12 months on dynamics global grid
[1270]	304	REAL, DIMENSION(klon_glo,12) :: psurf_glo1D ! -"- on physical global grid
[3413]	305	REAL, ALLOCATABLE :: load_glo2D(:,:,:) ! Load for 12 months on dynamics global grid
[1270]	306	REAL, DIMENSION(klon_glo,12) :: load_glo1D ! -"- on physical global grid
[3413]	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
[1183]	310	LOGICAL :: new_file ! true if new file format detected
	311	LOGICAL :: invert_lat ! true if the field has to be inverted for latitudes
[3413]	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
[1179]	338	! Deallocate pointers
	339	IF (ASSOCIATED(pt_ap)) DEALLOCATE(pt_ap)
	340	IF (ASSOCIATED(pt_b)) DEALLOCATE(pt_b)
[766]	341
[1249]	342	IF (is_mpi_root .AND. is_omp_root) THEN
[524]	343
[1179]	344	! 1) Open file
	345	!****************************************************************************************
[1492]	346	! Add suffix to filename
	347	fname = trim(filename)//cyr//'.nc'
[1150]	348
[1492]	349	WRITE(lunout,*) 'reading variable ',TRIM(varname),' in file ', TRIM(fname)
[2823]	350	CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid), "pb open "//trim(fname) )
[1150]	351
[3413]	352
	353	IF (grid_type/=unstructured) THEN
	354
[1183]	355	! Test for equal longitudes and latitudes in file and model
	356	!****************************************************************************************
	357	! Read and test longitudes
[3413]	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" )
[1183]	360
[3413]	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
[1183]	365
[3413]	366	CALL abort_physic('get_aero_fromfile', 'longitudes are not the same in file and model',1)
	367	END IF
[1183]	368
[3413]	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" )
[1183]	372
[3413]	373	! Invert source latitudes
	374	DO j = 1, nbp_lat
	375	lat_src_inv(j) = lat_src(nbp_lat +1 -j)
	376	END DO
[1183]	377
[3413]	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
[1183]	392
[1179]	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.
[2841]	400	ierr = nf90_inq_dimid(ncid, 'presnivs', dimid)
[1179]	401	IF (ierr /= NF90_NOERR) THEN
[2841]	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
[1179]	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	!****************************************************************************************
[1609]	424	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src),"pb inq dim for PRESNIVS or lev" )
[1179]	425
	426	! Allocate variables depending on the number of vertical levels
[2346]	427	ALLOCATE(varmth(nbp_lon,nbp_lat, klev_src), varyear(nbp_lon,nbp_lat, klev_src, 12), stat=ierr)
[2311]	428	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 1',1)
[1150]	429
[1179]	430	ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
[2311]	431	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 2',1)
[1150]	432
[1179]	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	!****************************************************************************************
[1150]	438
[1179]	439	IF (new_file) THEN
[1492]	440	! ++) Check number of month in file opened
	441	!**************************************************************************************************
	442	ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
[2841]	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" )
[1492]	447	! IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
	448	IF (nbr_tsteps /= 12 ) THEN
[2317]	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)
[1492]	451	ENDIF
[1150]	452
[1179]	453	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
	454	!****************************************************************************************
	455	! Get variable id
[2823]	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
[1179]	466
[1270]	467	! ++) Read surface pression, 12 month in one variable
[1179]	468	!****************************************************************************************
	469	! Get variable id
[1609]	470	CALL check_err( nf90_inq_varid(ncid, "ps", varid),"pb inq var ps" )
[1179]	471	! Get the variable
[1609]	472	CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D),"pb get var ps" )
[1179]	473
	474	! ++) Read mass load, 12 month in one variable
	475	!****************************************************************************************
	476	! Get variable id
[2823]	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
[1179]	485
	486	! ++) Read ap
	487	!****************************************************************************************
	488	! Get variable id
[1609]	489	CALL check_err( nf90_inq_varid(ncid, "ap", varid),"pb inq var ap" )
[1179]	490	! Get the variable
[1609]	491	CALL check_err( nf90_get_var(ncid, varid, pt_ap),"pb get var ap" )
[1150]	492
[1179]	493	! ++) Read b
	494	!****************************************************************************************
	495	! Get variable id
[1609]	496	CALL check_err( nf90_inq_varid(ncid, "b", varid),"pb inq var b" )
[1179]	497	! Get the variable
[1609]	498	CALL check_err( nf90_get_var(ncid, varid, pt_b),"pb get var b" )
[1150]	499
	500
[524]	501
[1179]	502	ELSE ! old file
[524]	503
[1179]	504	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
	505	!****************************************************************************************
[1270]	506	DO imth=1, 12
[1179]	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
[1609]	534	CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid),"pb inq var "//TRIM(cvar) )
[1179]	535
	536	! Get the variable
[1609]	537	CALL check_err( nf90_get_var(ncid, varid, varmth),"pb get var "//TRIM(cvar) )
[1179]	538
	539	! Store in variable for the whole year
	540	varyear(:,:,:,imth)=varmth(:,:,:)
	541
	542	END DO
[1150]	543
[1179]	544	! Putting dummy
	545	psurf_glo2D(:,:,:) = not_valid
	546	load_glo2D(:,:,:) = not_valid
	547	pt_ap(:) = not_valid
	548	pt_b(:) = not_valid
[524]	549
[1179]	550	END IF
[524]	551
[1179]	552	! 4) Close file
	553	!****************************************************************************************
[1609]	554	CALL check_err( nf90_close(ncid),"pb in close" )
[1179]	555
[524]	556
[2346]	557	! 5) Transform the global field from 2D(nbp_lon,nbp_lat) to 1D(klon_glo)
[1179]	558	!****************************************************************************************
	559	! Test if vertical levels have to be inversed
[782]	560
[1179]	561	IF ((pt_b(1) < pt_b(klev_src)) .OR. .NOT. new_file) THEN
[1321]	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
[1179]	565
	566	! Inverse vertical levels for varyear
[1270]	567	DO imth=1, 12
[1179]	568	varmth(:,:,:) = varyear(:,:,:,imth) ! use varmth temporarly
	569	DO k=1, klev_src
[2346]	570	DO j=1, nbp_lat
	571	DO i=1, nbp_lon
[1179]	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
[782]	583
[1179]	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
[524]	590
[1179]	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
[524]	596
[1183]	597
[3413]	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
[1183]	612
[3413]	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
[1183]	620
[3413]	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
[1183]	629
[3413]	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
[1183]	644
[3413]	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)
[1179]	647
[3413]	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
[1249]	655	ELSE
[3413]	656	ALLOCATE(varyear_glo1D(0,0,0))
[1249]	657	END IF ! is_mpi_root .AND. is_omp_root
	658
[1189]	659	!$OMP BARRIER
[1179]	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	!****************************************************************************************
[524]	665
[1179]	666	! Distribute klev_src
	667	CALL bcast(klev_src)
[524]	668
[1179]	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)
[2311]	672	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 4',1)
[1179]	673	END IF
	674	CALL bcast(pt_ap)
	675	CALL bcast(pt_b)
[524]	676
[1179]	677	! Allocate space for output pointer variable at local process
	678	IF (ASSOCIATED(pt_year)) DEALLOCATE(pt_year)
[1270]	679	ALLOCATE(pt_year(klon, klev_src, 12), stat=ierr)
[3413]	680	ALLOCATE(pt_year_mpi(klon_mpi, klev_src, 12), stat=ierr)
[2311]	681	IF (ierr /= 0) CALL abort_physic('get_aero_fromfile', 'pb in allocation 5',1)
[524]	682
[3413]	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
[1179]	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
[524]	710
[1179]	711	END SUBROUTINE get_aero_fromfile
[524]	712
	713
[1609]	714	SUBROUTINE check_err(status,text)
[1179]	715	USE netcdf
[2311]	716	USE print_control_mod, ONLY: lunout
[1179]	717	IMPLICIT NONE
[524]	718
[1179]	719	INTEGER, INTENT (IN) :: status
[1609]	720	CHARACTER(len=*), INTENT (IN), OPTIONAL :: text
[524]	721
[1179]	722	IF (status /= NF90_NOERR) THEN
[1609]	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
[2311]	727	CALL abort_physic('get_aero_fromfile',trim(nf90_strerror(status)),1)
[1179]	728	END IF
[524]	729
[1179]	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