Context Navigation

readaerosol.F90 @ 1483

Last change on this file since 1483 was 1483, checked in by jghattas, 13 years ago

Added option aer_type=mix1 and mix2

mix1: as aer_type=scenario for SO4 and aer_type=annuel for the other aerosols.
mix2: as aer_type=scenario for SO4 and aer_type=preind for the other aerosols.
For these 2 options 2 files are read : so4.runxxxx.nc for SO4 and aerosolsxxxx.nc the other aerosols.

Bug fix for using old sulfate data(files AR4 from O Boucher) : moved test for number of time steps read from file but not available in old file format.

Property svn:eol-style set to native
Property svn:keywords set to Author Date Id Revision

File size: 23.9 KB

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

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

Download in other formats:

Original Format