Context Navigation

readaerosol.F90 @ 1321

Last change on this file since 1321 was 1321, checked in by Laurent Fairhead, 15 years ago

The ability to read annual aerosol forcing files that had disappeared in
recent versions of the code is put back in

On ré-intégre la possibilité de lire des fichiers de forçage aérosols de
type annuel

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

File size: 23.6 KB

Rev	Line
[1179]	1	! $Id: readaerosol.F90 1321 2010-03-09 14:42:32Z fairhead $
[524]	2	!
[1179]	3	MODULE readaerosol_mod
	4
	5	REAL, SAVE :: not_valid=-333.
	6
	7	CONTAINS
	8
[1270]	9	SUBROUTINE readaerosol(name_aero, type, 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
	29	CHARACTER(len=*), INTENT(IN) :: type ! correspond to aer_type in clesphys.h
	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)
	60	CALL get_aero_fromfile(name_aero, cyear, 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)
	69	CALL get_aero_fromfile(name_aero, cyear, 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)
	78	CALL get_aero_fromfile(name_aero, cyear, klev_src, pt_ap, pt_b, p0, pt_out, psurf, load)
	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)
	88	CALL get_aero_fromfile(name_aero, cyear, 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)
	95	CALL get_aero_fromfile(name_aero, cyear, 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)
	115	CALL get_aero_fromfile(name_aero, cyear, 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)
	127	CALL get_aero_fromfile(name_aero, cyear, 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'
	131	CALL abort_gcm('readaersosol','Error in number of vertical levels',1)
	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) = &
	139	pt_out(i,k,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
	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) = &
	146	psurf(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
	147	(psurf(i,it) - psurf2(i,it))
[640]	148
[1270]	149	load(i,it) = &
	150	load(i,it) - FLOAT(iyr_in-iyr1)/FLOAT(iyr2-iyr1) * &
	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
	162	WRITE(lunout,*)'This option is not implemented : aer_type = ', type
	163	CALL abort_gcm('readaerosol','Error : aer_type parameter not accepted',1)
	164	END IF ! type
[640]	165
[1151]	166
[1179]	167	END SUBROUTINE readaerosol
[1143]	168
	169
[1270]	170	SUBROUTINE get_aero_fromfile(varname, cyr, 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
[1179]	193
	194	IMPLICIT NONE
[524]	195
[1179]	196	INCLUDE "dimensions.h"
	197	INCLUDE "iniprint.h"
[524]	198
[1179]	199	! Input argumets
	200	CHARACTER(len=7), INTENT(IN) :: varname
	201	CHARACTER(len=4), INTENT(IN) :: cyr
[1150]	202
[1179]	203	! Output arguments
	204	INTEGER, INTENT(OUT) :: klev_src ! Number of vertical levels in file
	205	REAL, POINTER, DIMENSION(:) :: pt_ap ! Pointer for describing the vertical levels
	206	REAL, POINTER, DIMENSION(:) :: pt_b ! Pointer for describing the vertical levels
	207	REAL :: p0 ! Reference pressure value
[1270]	208	REAL, POINTER, DIMENSION(:,:,:) :: pt_year ! Pointer-variabale from file, 12 month, grid : klon,klev_src
	209	REAL, DIMENSION(klon,12), INTENT(OUT) :: psurf_out ! Surface pression for 12 months
	210	REAL, DIMENSION(klon,12), INTENT(OUT) :: load_out ! Aerosol mass load in each column
[1319]	211	INTEGER :: nbr_tsteps ! number of month in file read
[1150]	212
[1179]	213	! Local variables
	214	CHARACTER(len=30) :: fname
[1223]	215	CHARACTER(len=8) :: filename='aerosols'
[1179]	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	!****************************************************************************************
[1223]	244	fname = filename//cyr//'.nc'
[1150]	245
[1179]	246	WRITE(lunout,*) 'reading ', TRIM(fname)
	247	CALL check_err( nf90_open(TRIM(fname), NF90_NOWRITE, ncid) )
[1150]	248
[1183]	249	! Test for equal longitudes and latitudes in file and model
	250	!****************************************************************************************
	251	! Read and test longitudes
	252	CALL check_err( nf90_inq_varid(ncid, 'lon', varid) )
	253	CALL check_err( nf90_get_var(ncid, varid, lon_src(:)) )
	254
[1202]	255	IF (maxval(ABS(lon_src - io_lon)) > 0.001) THEN
[1183]	256	WRITE(lunout,*) 'Problem in longitudes read from file : ',TRIM(fname)
	257	WRITE(lunout,*) 'longitudes in file ', TRIM(fname),' : ', lon_src
	258	WRITE(lunout,*) 'longitudes in model :', io_lon
	259
	260	CALL abort_gcm('get_aero_fromfile', 'longitudes are not the same in file and model',1)
	261	END IF
	262
	263	! Read and test latitudes
	264	CALL check_err( nf90_inq_varid(ncid, 'lat', varid) )
	265	CALL check_err( nf90_get_var(ncid, varid, lat_src(:)) )
	266
	267	! Invert source latitudes
	268	DO j = 1, jjm+1
	269	lat_src_inv(j) = lat_src(jjm+1 +1 -j)
	270	END DO
	271
[1202]	272	IF (maxval(ABS(lat_src - io_lat)) < 0.001) THEN
[1183]	273	! Latitudes are the same
	274	invert_lat=.FALSE.
[1202]	275	ELSE IF (maxval(ABS(lat_src_inv - io_lat)) < 0.001) THEN
[1183]	276	! Inverted source latitudes correspond to model latitudes
	277	WRITE(lunout,*) 'latitudes will be inverted for file : ',TRIM(fname)
	278	invert_lat=.TRUE.
	279	ELSE
	280	WRITE(lunout,*) 'Problem in latitudes read from file : ',TRIM(fname)
	281	WRITE(lunout,*) 'latitudes in file ', TRIM(fname),' : ', lat_src
	282	WRITE(lunout,*) 'latitudes in model :', io_lat
	283	CALL abort_gcm('get_aero_fromfile', 'latitudes do not correspond between file and model',1)
	284	END IF
	285
[1319]	286	! 1.5) Check number of month in file opened
	287	!
	288	!**************************************************************************************************
	289	ierr = nf90_inq_dimid(ncid, 'TIME',dimid)
	290	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = nbr_tsteps) )
	291	! IF (nbr_tsteps /= 12 .AND. nbr_tsteps /= 14) THEN
	292	IF (nbr_tsteps /= 12 ) THEN
	293	CALL abort_gcm('get_aero_fromfile', 'not the right number of months in aerosol file read (should be 12 for the moment)',1)
	294	ENDIF
	295
	296
[1179]	297	! 2) Check if old or new file is avalabale.
	298	! New type of file should contain the dimension 'lev'
	299	! Old type of file should contain the dimension 'PRESNIVS'
	300	!****************************************************************************************
	301	ierr = nf90_inq_dimid(ncid, 'lev', dimid)
	302	IF (ierr /= NF90_NOERR) THEN
	303	! Coordinate axe lev not found. Check for presnivs.
	304	ierr = nf90_inq_dimid(ncid, 'PRESNIVS', dimid)
	305	IF (ierr /= NF90_NOERR) THEN
	306	! Dimension PRESNIVS not found either
	307	CALL abort_gcm('get_aero_fromfile', 'dimension lev or presnivs not in file',1)
	308	ELSE
	309	! Old file found
	310	new_file=.FALSE.
	311	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will not be done'
	312	END IF
	313	ELSE
	314	! New file found
	315	new_file=.TRUE.
	316	WRITE(lunout,*) 'Vertical interpolation for ',TRIM(varname),' will be done'
	317	END IF
	318
	319	! 2) Find vertical dimension klev_src
	320	!****************************************************************************************
	321	CALL check_err( nf90_inquire_dimension(ncid, dimid, len = klev_src) )
	322
	323	! Allocate variables depending on the number of vertical levels
[1270]	324	ALLOCATE(varmth(iim, jjm+1, klev_src), varyear(iim, jjm+1, klev_src, 12), stat=ierr)
[1179]	325	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 1',1)
[1150]	326
[1179]	327	ALLOCATE(pt_ap(klev_src), pt_b(klev_src), varktmp(klev_src), stat=ierr)
	328	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 2',1)
[1150]	329
[1179]	330	! 3) Read all variables from file
	331	! There is 2 options for the file structure :
	332	! new_file=TRUE : read varyear, ps, pt_ap and pt_b
	333	! new_file=FALSE : read varyear month by month
	334	!****************************************************************************************
[1150]	335
[1179]	336	IF (new_file) THEN
[1150]	337
[1179]	338	! ++) Read the aerosol concentration month by month and concatenate to total variable varyear
	339	!****************************************************************************************
	340	! Get variable id
	341	CALL check_err( nf90_inq_varid(ncid, TRIM(varname), varid) )
	342
	343	! Get the variable
	344	CALL check_err( nf90_get_var(ncid, varid, varyear(:,:,:,:)) )
	345
[1270]	346	! ++) Read surface pression, 12 month in one variable
[1179]	347	!****************************************************************************************
	348	! Get variable id
	349	CALL check_err( nf90_inq_varid(ncid, "ps", varid) )
	350	! Get the variable
	351	CALL check_err( nf90_get_var(ncid, varid, psurf_glo2D) )
	352
	353	! ++) Read mass load, 12 month in one variable
	354	!****************************************************************************************
	355	! Get variable id
	356	CALL check_err( nf90_inq_varid(ncid, "load_"//TRIM(varname), varid) )
	357	! Get the variable
	358	CALL check_err( nf90_get_var(ncid, varid, load_glo2D) )
	359
	360	! ++) Read ap
	361	!****************************************************************************************
	362	! Get variable id
	363	CALL check_err( nf90_inq_varid(ncid, "ap", varid) )
	364	! Get the variable
	365	CALL check_err( nf90_get_var(ncid, varid, pt_ap) )
[1150]	366
[1179]	367	! ++) Read b
	368	!****************************************************************************************
	369	! Get variable id
	370	CALL check_err( nf90_inq_varid(ncid, "b", varid) )
	371	! Get the variable
	372	CALL check_err( nf90_get_var(ncid, varid, pt_b) )
[1150]	373
[1179]	374	! ++) Read p0 : reference pressure
	375	!****************************************************************************************
	376	! Get variable id
	377	CALL check_err( nf90_inq_varid(ncid, "p0", varid) )
	378	! Get the variable
	379	CALL check_err( nf90_get_var(ncid, varid, 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
	414	CALL check_err( nf90_inq_varid(ncid, TRIM(cvar), varid) )
	415
	416	! Get the variable
	417	CALL check_err( nf90_get_var(ncid, varid, varmth) )
	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	!****************************************************************************************
	434	CALL check_err( nf90_close(ncid) )
	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
	516	npole = npole/FLOAT(iim)
	517	spole = spole/FLOAT(iim)
	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)
[1179]	524	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 3',1)
	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)
	548	IF (ierr /= 0) CALL abort_gcm('get_aero_fromfile', 'pb in allocation 4',1)
	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)
	556	IF (ierr /= 0) CALL abort_gcm('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
[1179]	572	SUBROUTINE check_err(status)
	573	USE netcdf
	574	IMPLICIT NONE
[524]	575
[1179]	576	INCLUDE "iniprint.h"
	577	INTEGER, INTENT (IN) :: status
[524]	578
[1179]	579	IF (status /= NF90_NOERR) THEN
	580	WRITE(lunout,*) 'Error in get_aero_fromfile ',status
	581	CALL abort_gcm('get_aero_fromfile',trim(nf90_strerror(status)),1)
	582	END IF
[524]	583
[1179]	584	END SUBROUTINE check_err
	585
	586
	587	END MODULE readaerosol_mod

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

Download in other formats:

Original Format