Context Navigation

readaerosol.F90 @ 2339

Last change on this file since 2339 was 2317, checked in by jyg, 9 years ago
Bug fixing after physics/dynamic separation enforcing in svn 2311 (corrections in readaerosol.F90 and in dyn1d/lmdz1d.F90).
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: 24.5 KB

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

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

Download in other formats:

Original Format