Context Navigation

readaerosol.F90 @ 4799

Last change on this file since 4799 was 2841, checked in by acozic, 8 years ago
remove the reading of p0 in aerosol files because it's not use and add the possibility to have new axes in these files (presnivs and time_counter)
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: 25.1 KB

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

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

Download in other formats:

Original Format